Comment 4 for bug 1082370

Thanks for the tip regarding curl.

I believe the choice of function spaces used means the solve should be
exact. I.e. the numerical error should be at the level of machine
precision, as with the UnitSquare case.

Here are three more meshes, with char. mesh length set to 1/8, 1/16, 1/32
in the obvious way. I suggest changing the last few lines to

#print max(q_.vector())
print assemble(q_*q_*dx)

This gives ~3.5 for test4.xml, ~2.5 for test8.xml, but then 600 and 1800
for test16.xml and test32.xml. So refinement isn't helping!

I plan to try some other, more structured, meshes this weekend.

On 23 November 2012 16:14, Mikael Mortensen <email address hidden>wrote:

> Andrew,
>
> I don't think there's anything wrong here, probably just numerical error
> due to poor mesh. Refine your mesh once and see what happens.
>
> Also, you can use
>
> L = inner(curl(gamma),du)*dx
>
> instead of what you have.
>
> Best regards
>
> Mikael
>
>
> Den Nov 23, 2012 kl. 3:43 PM skrev Andrew McRae:
>
> > Sorry! I forgot to attach it when relisting the bug. Attached to this
> > comment.
> >
> > ** Attachment added: "Mesh"
> >
> https://bugs.launchpad.net/dolfin/+bug/1082370/+attachment/3442192/+files/test4.xml
> >
> > --
> > You received this bug notification because you are a member of DOLFIN
> > Team, which is subscribed to DOLFIN.
> > https://bugs.launchpad.net/bugs/1082370
> >
> > Title:
> > Inaccurate solve on unstructured mesh
> >
> > Status in DOLFIN:
> > New
> >
> > Bug description:
> > Mesh was generated in gmsh (simple unit square with characteristic
> > mesh length 0.25 everywhere), converted to xml using dolfin-convert.
> >
> > Based on the linear SW equations. The code below is equivalent to
> > starting with a given height field, and a zero velocity field, then
> > updating the velocity field by a single Euler timestep. The vorticity
> > is calculated. This should be numerically 0 (change the mesh to
> > UnitSquare to see this).
> >
> > from dolfin import *
> > set_log_level(ERROR)
> > # mesh and function spaces
> > mesh = Mesh("test4.xml")
> > #mesh = UnitSquare(4,4)
> > plot(mesh)
> > E = FunctionSpace(mesh, 'CG', 2)
> > S = FunctionSpace(mesh, 'BDM', 1)
> > V = FunctionSpace(mesh, 'DG', 0)
> >
> > class PeriodicBoundaryX(SubDomain):
> > def inside(self, x, on_boundary):
> > return on_boundary and near(x[0], 0)
> > def map(self, x, y):
> > y[0] = x[0] - 1
> > y[1] = x[1]
> >
> > class PeriodicBoundaryY(SubDomain):
> > def inside(self, x, on_boundary):
> > return on_boundary and near(x[1], 0)
> > def map(self, x, y):
> > y[0] = x[0]
> > y[1] = x[1] - 1
> >
> > left = PeriodicBoundaryX()
> > bottom = PeriodicBoundaryY()
> > bcE = [PeriodicBC(E, left), PeriodicBC(E, bottom)]
> > bcS = [PeriodicBC(S, left), PeriodicBC(S, bottom)]
> >
> > Dexpr = Expression("1+10*pow(x[0],2)*pow(x[1],3)*(1-x[0])*(1-x[1])")
> > D = Function(interpolate(Dexpr,V))
> >
> > dui = TrialFunction(S)
> > w = TestFunction(S)
> > a = inner(w,dui)*dx
> > L = div(w)*D*dx
> > A = assemble(a)
> > b = assemble(L)
> > [bc.apply(A, b) for bc in bcS]
> > du = Function(S)
> > solve(A, du.vector(), b)
> >
> > q = TrialFunction(E)
> > gamma = TestFunction(E)
> > a = gamma*q*dx
> > L = inner(as_vector([gamma.dx(1),-gamma.dx(0)]),du)*dx
> > A = assemble(a)
> > b = assemble(L)
> > [bc.apply(A, b) for bc in bcE]
> > q_ = Function(E)
> > solve(A, q_.vector(), b)
> >
> > print q_.vector().array() # should be 0
> >
> > interactive()
> >
> > To manage notifications about this bug go to:
> > https://bugs.launchpad.net/dolfin/+bug/1082370/+subscriptions
>
> --
> You received this bug notification because you are subscribed to the bug
> report.
> https://bugs.launchpad.net/bugs/1082370
>
> Title:
> Inaccurate solve on unstructured mesh
>
> Status in DOLFIN:
> New
>
> Bug description:
> Mesh was generated in gmsh (simple unit square with characteristic
> mesh length 0.25 everywhere), converted to xml using dolfin-convert.
>
> Based on the linear SW equations. The code below is equivalent to
> starting with a given height field, and a zero velocity field, then
> updating the velocity field by a single Euler timestep. The vorticity
> is calculated. This should be numerically 0 (change the mesh to
> UnitSquare to see this).
>
> from dolfin import *
> set_log_level(ERROR)
> # mesh and function spaces
> mesh = Mesh("test4.xml")
> #mesh = UnitSquare(4,4)
> plot(mesh)
> E = FunctionSpace(mesh, 'CG', 2)
> S = FunctionSpace(mesh, 'BDM', 1)
> V = FunctionSpace(mesh, 'DG', 0)
>
> class PeriodicBoundaryX(SubDomain):
> def inside(self, x, on_boundary):
> return on_boundary and near(x[0], 0)
> def map(self, x, y):
> y[0] = x[0] - 1
> y[1] = x[1]
>
> class PeriodicBoundaryY(SubDomain):
> def inside(self, x, on_boundary):
> return on_boundary and near(x[1], 0)
> def map(self, x, y):
> y[0] = x[0]
> y[1] = x[1] - 1
>
> left = PeriodicBoundaryX()
> bottom = PeriodicBoundaryY()
> bcE = [PeriodicBC(E, left), PeriodicBC(E, bottom)]
> bcS = [PeriodicBC(S, left), PeriodicBC(S, bottom)]
>
> Dexpr = Expression("1+10*pow(x[0],2)*pow(x[1],3)*(1-x[0])*(1-x[1])")
> D = Function(interpolate(Dexpr,V))
>
> dui = TrialFunction(S)
> w = TestFunction(S)
> a = inner(w,dui)*dx
> L = div(w)*D*dx
> A = assemble(a)
> b = assemble(L)
> [bc.apply(A, b) for bc in bcS]
> du = Function(S)
> solve(A, du.vector(), b)
>
> q = TrialFunction(E)
> gamma = TestFunction(E)
> a = gamma*q*dx
> L = inner(as_vector([gamma.dx(1),-gamma.dx(0)]),du)*dx
> A = assemble(a)
> b = assemble(L)
> [bc.apply(A, b) for bc in bcE]
> q_ = Function(E)
> solve(A, q_.vector(), b)
>
> print q_.vector().array() # should be 0
>
> interactive()
>
> To manage notifications about this bug go to:
> https://bugs.launchpad.net/dolfin/+bug/1082370/+subscriptions
>