- changed status to invalid
Error in mixed hybrid Poisson equation with discontinuous Lagrange trace
Hi.
I tried to implement hybridized mixed method for the Poisson equation using discontinuous Lagrange trace space. However, in contrast to standard mixed method, the method does not give a correct answer. I used Dolfin 1.6.0. Please take a look.
Added: Standard mixed method and the hybridized one in this code have to give same answer for sigma and u but errors of the numerical solutions with this code are "nan".
Regards, Jeonghun J. Lee
# hybridized mixed method for Poisson equation
from dolfin import *
import numpy
#exact solution
u_ex = Expression('x[0]*(1-x[0])*sin(pi*x[1])', degree = 3)
sigma_ex = Expression(('(1 - 2*x[0])*sin(pi*x[1])', 'pi*x[0]*(1-x[0])*cos(pi*x[1])'), degree = 3)
f = Expression('-sin(pi*x[1])*(2+pi*pi*x[0]*(1-x[0]))', degree = 3)
def compute(nx):
mesh = UnitSquareMesh(nx,nx)
# define function spaces
Mh = VectorElement("DG", mesh.ufl_cell(), 1)
Vh = FiniteElement("DG", mesh.ufl_cell(), 0)
Qh = FiniteElement("Discontinuous Lagrange Trace", mesh.ufl_cell(), 1)
element = MixedElement([Mh,Vh,Qh])
V = FunctionSpace(mesh, element)
#def uex_bdy(x, on_boundary):
# return on_boundary
bc = DirichletBC(V.sub(2), u_ex, "on_boundary")
# test and trial functions
sigma, u, p = TrialFunctions(V)
tau, v, q = TestFunctions(V)
# bilinear form
n = FacetNormal(mesh)
a0 = inner(sigma, tau)*dx + u*div(tau)*dx + div(sigma)*v*dx
a1 = avg(p)*jump(tau,n)*dS + jump(sigma,n)*avg(q)*dS + p*dot(tau, n)*ds + dot(sigma,n)*q*ds
### To verify convergence of standard mixed method
# Mh = FiniteElement('BDM', mesh.ufl_cell(), 1)
# Vh = FiniteElement('DG', mesh.ufl_cell(), 0)
# element = MixedElement([Mh, Vh])
# V = FunctionSpace(mesh, element)
#
# (tau, v) = TestFunctions(V)
# (sigma, u) = TrialFunctions(V)
#
# a = inner(sigma, tau)*dx + inner(u, div(tau))*dx + inner(div(sigma), v)*dx
###
L = f*v*dx
#A = assemble(a)
A0 = assemble(a0)
A1 = assemble(a1)
A = A0 + A1
b = assemble(L)
bc.apply(A, b)
U = Function(V)
solve(A, U.vector(), b)
(sigma, u, p) = U.split(deepcopy = True)
# Compute errors
V1 = FunctionSpace(mesh, 'BDM', 2)
sigma_i = interpolate(sigma_ex, V1)
V2 = FunctionSpace(mesh, 'DG', 1)
u_i = interpolate(u_ex, V2)
sigma_err = sqrt(abs(assemble(inner(sigma_i - sigma, sigma_i - sigma)*dx)))
u_err = sqrt(assemble((u_i - u)*(u_i - u)*dx))
errors = {'u_L2_err': u_err, 'sigma_L2_err': sigma_err}
print 'u_err', u_err
print 'sigma_err', sigma_err
return errors
h = [] # element sizes
E = [] # errors
for nx in [4, 8, 16, 32, 64]:
h.append(1.0/nx)
E.append(compute(nx)) # list of dicts
# Compute convergence rates
from math import log as ln # log is a dolfin name too
error_types = E[0].keys()
for error_type in sorted(error_types):
l = len(E)
print '\nError norm based on', error_type
for i in range(0, l):
if i==0:
print 'Mesh size=%.10f Error=%.10f Convergence rate= -- ' % (h[i], E[i][error_type])
else:
r = ln(E[i][error_type]/E[i-1][error_type])/ln(0.5)
print 'Mesh size=%.10f Error=%.10f Convergence rate=%.2f' % (h[i], E[i][error_type], r)
Comments (15)
-
-
-
assigned issue to
John, could you please add code formatting to your example.
Garth, I am unsure if this is a bug and asked John to report it here to keep track of it.
-
assigned issue to
-
- edited description
-
The report doesn't say what the issue is, just 'does not give a correct answer'. It needs say what the expected answer is. Also, the code should be paired back.
-
reporter @Johannes : I think you changed the code format. Thank you for doing that.
-
reporter - edited description
-
@johnlee04 - Yes, I added the missing ``` at the top. You are welcome.
-
-
- changed status to open
-
reporter - edited description
-
reporter Updated the code for 2017.1.0
-
- marked as blocker
-
@johnlee04 FIAT with the fix to this bug https://bitbucket.org/fenics-project/dolfin/issues/871/interpolation-onto-discontinuous-lagrange, in e.g. this branch: https://bitbucket.org/fenics-project/fiat/branch/meg/fix-dolfin-issue-871, yields convergence rates of 2 for sigma in L2 and 1 for u in L2 for your test case above. I will push the issue 871 fix to master asap.
-
Testing these issues more closely, this seems to run as expected with current master. Interesting... Marking as resolved and removing blocker.
-
- changed status to resolved
Example converges as expected with current master, resolving.
- Log in to comment
Please use fenics-support@googlegroups.com