Inaccurate solve on unstructured mesh
Affects | Status | Importance | Assigned to | Milestone | |
---|---|---|---|---|---|
DOLFIN |
Invalid
|
Undecided
|
Unassigned |
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_
# 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 PeriodicBoundar
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 PeriodicBoundar
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(
D = Function(
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_
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()
Changed in dolfin: | |
status: | New → Invalid |
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: level(ERROR) yX(SubDomain) : yY(SubDomain) : "1+10*pow( x[0],2) *pow(x[ 1],3)*( 1-x[0]) *(1-x[1] )") interpolate( Dexpr,V) ) vector( [gamma. dx(1),- gamma.dx( 0)]),du) *dx /bugs.launchpad .net/bugs/ 1082370 level(ERROR)
>
> 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_
> # 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 PeriodicBoundar
> 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 PeriodicBoundar
> 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(
> D = Function(
>
> 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_
> 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:/
>
> 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_
> # 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...