Inaccurate solve on unstructured mesh

[Andrew McRae]

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()

[Mikael Mortensen (mikael-mortensen)]

Hi Andrew,

You need to provide the mesh for us to test it. You could probably just copy-paste it into the top of you code or something since I don't think you can attach files to bugs?

Mikael

Den Nov 23, 2012 kl. 2:28 PM skrev Andrew McRae:

> Public bug reported:
>
> 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()
>
> ** Affects: dolfin
> Importance: Undecided
> Status: New
>
> --
> 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...

[Andrew McRae (andymc)]

Sorry! I forgot to attach it when relisting the bug. Attached to this comment.

[Mikael Mortensen (mikael-mortensen)]

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

[Andrew McRae (andymc)]

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)
...

[Mikael Mortensen (mikael-mortensen)]

Try

mesh = refine(mesh)

This will make the result exact. Perhaps there is something wrong with the mesh you have that refine fixes. I tried mesh.order(), but that changes nothing.

Mikael

Den 23. nov. 2012 kl. 18:24 skrev Andrew McRae <email address hidden>:

> 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)]
>>>
...

[Andrew McRae (andymc)]

I briefly tried it and it indeed seemed to work. Thanks! May I ask what
"refine" actually does? Naively, I'd expect it to require both a mesh and
a function for the mesh to be refined wrt., but here it's being called on
the mesh only.

Andrew

On 23 November 2012 17:59, Mikael Mortensen <email address hidden>wrote:

> Try
>
> mesh = refine(mesh)
>
> This will make the result exact. Perhaps there is something wrong with
> the mesh you have that refine fixes. I tried mesh.order(), but that
> changes nothing.
>
> Mikael
>
> Den 23. nov. 2012 kl. 18:24 skrev Andrew McRae
> <email address hidden>:
>
> > 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], ...