Comment 3 for bug 1082370

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