mesh hierarchy inject function issue - possible bug
Dear team, In the following code, I call inject twice, each time on separate sets of variables {lag_f, lag_c} and {psi0, psi0_c}, where lag_f and psi0 are defined on the fine mesh, and lag_c and psi0_c are defined on the coarse mesh. It seems that the first inject call on {psi0, psi0_c} changes the result of the second inject call on {lag_f, lag_c} which is not what I would expect as psi0 and lag are unrelated. I have highlighted the relevant definitions in blue, and problem bits in red. I should add that the problem goes away if in the VectorFunctionSpace definitions of vfs_f and vfs_c I change "CG" to "DG". But I'm worried about the result of the CG inject. Many thanks in advance. Kind regards Wei from firedrake import * import numpy as np L = 1.0 # Zonal length n0 = 50 # Spatial resolution mesh_base = SquareMesh(n0, n0, L) mesh_hierarchy = MeshHierarchy(mesh_base, 1) mesh = mesh_hierarchy[1] # fine mesh mesh_c = mesh_hierarchy[0] # coarse mesh Vdg = FunctionSpace(mesh,"DG",1) # DG elements for Potential Vorticity (PV) Vcg = FunctionSpace(mesh,"CG",1) # CG elements for Streamfunction vfs_f = VectorFunctionSpace(mesh,"CG",1) vfs_c = VectorFunctionSpace(mesh_c,"CG",1) vxf, vyf = SpatialCoordinate(vfs_f) x = SpatialCoordinate(mesh) q0 = Function(Vdg).interpolate(0.1*sin(x[0])*sin(x[1])) q1 = Function(Vdg) psi0 = Function(Vcg) # Streamfunctions for different time steps psi = TrialFunction(Vcg) # Test function phi = TestFunction(Vcg) # Trial function Apsi = (inner(grad(psi),grad(phi)))*dx Lpsi = -q1*phi*dx bc1 = DirichletBC(Vcg, 0., 'on_boundary') psi_problem = LinearVariationalProblem(Apsi,Lpsi,psi0,bcs=bc1) psi_solver = LinearVariationalSolver(psi_problem, solver_parameters={ 'ksp_type':'cg', 'pc_type':'sor' }) gradperp = lambda u: as_vector((-u.dx(1), u.dx(0))) lag_f = Function(vfs_f).interpolate(as_vector((vxf, vyf))) if 1: print "-------------------------------------------------------------" q1.assign(q0) psi_solver.solve() psi0_coarse = Function(FunctionSpace(mesh_c, "CG", 1)) inject(psi0, psi0_coarse) print "psi0 norm ", norm(psi0) print "psi0_coarse norm ", norm(psi0_coarse) print "psi0_coarse relative error ", np.absolute(norm(psi0_coarse) - norm(psi0)) / norm(psi0) # relative error 0.000583086120876 print "-------------------------------------------------------------" lag_c = Function(vfs_c) inject(lag_f, lag_c) print "analytical value ", np.sqrt(2./3.) # 0.816496580928 print "lag0_f norm ", norm(lag_f) print "lag0_c norm ", norm(lag_c) # this is the correct value of 0.816496580928 if inject(psi0, psi0_coarse) is not called # but if inject(psi0, psi0_coarse) is called then # this gives 0.825399087311 which is a relative error of 0.019032990351 print "lag relative error ", np.absolute(norm(lag_c) - norm(lag_f)) / norm(lag_f)
Hmm, that does sound like a bug. I will get a chance to look at this on Tuesday. Thanks! Lawrence
On 19 May 2017, at 19:02, Pan, Wei <wei.pan@imperial.ac.uk> wrote:
Dear team,
In the following code, I call inject twice, each time on separate sets of variables {lag_f, lag_c} and {psi0, psi0_c}, where lag_f and psi0 are defined on the fine mesh, and lag_c and psi0_c are defined on the coarse mesh. It seems that the first inject call on {psi0, psi0_c} changes the result of the second inject call on {lag_f, lag_c} which is not what I would expect as psi0 and lag are unrelated. I have highlighted the relevant definitions in blue, and problem bits in red.
I should add that the problem goes away if in the VectorFunctionSpace definitions of vfs_f and vfs_c I change "CG" to "DG". But I'm worried about the result of the CG inject. Many thanks in advance.
Kind regards
Wei
from firedrake import * import numpy as np
L = 1.0 # Zonal length n0 = 50 # Spatial resolution
mesh_base = SquareMesh(n0, n0, L) mesh_hierarchy = MeshHierarchy(mesh_base, 1) mesh = mesh_hierarchy[1] # fine mesh mesh_c = mesh_hierarchy[0] # coarse mesh
Vdg = FunctionSpace(mesh,"DG",1) # DG elements for Potential Vorticity (PV) Vcg = FunctionSpace(mesh,"CG",1) # CG elements for Streamfunction
vfs_f = VectorFunctionSpace(mesh,"CG",1) vfs_c = VectorFunctionSpace(mesh_c,"CG",1) vxf, vyf = SpatialCoordinate(vfs_f)
x = SpatialCoordinate(mesh) q0 = Function(Vdg).interpolate(0.1*sin(x[0])*sin(x[1])) q1 = Function(Vdg) psi0 = Function(Vcg) # Streamfunctions for different time steps
psi = TrialFunction(Vcg) # Test function phi = TestFunction(Vcg) # Trial function
Apsi = (inner(grad(psi),grad(phi)))*dx Lpsi = -q1*phi*dx bc1 = DirichletBC(Vcg, 0., 'on_boundary') psi_problem = LinearVariationalProblem(Apsi,Lpsi,psi0,bcs=bc1) psi_solver = LinearVariationalSolver(psi_problem, solver_parameters={ 'ksp_type':'cg', 'pc_type':'sor' })
gradperp = lambda u: as_vector((-u.dx(1), u.dx(0)))
lag_f = Function(vfs_f).interpolate(as_vector((vxf, vyf)))
if 1: print "-------------------------------------------------------------"
q1.assign(q0) psi_solver.solve() psi0_coarse = Function(FunctionSpace(mesh_c, "CG", 1)) inject(psi0, psi0_coarse)
print "psi0 norm ", norm(psi0) print "psi0_coarse norm ", norm(psi0_coarse) print "psi0_coarse relative error ", np.absolute(norm(psi0_coarse) - norm(psi0)) / norm(psi0) # relative error 0.000583086120876
print "-------------------------------------------------------------"
lag_c = Function(vfs_c) inject(lag_f, lag_c) print "analytical value ", np.sqrt(2./3.) # 0.816496580928 print "lag0_f norm ", norm(lag_f) print "lag0_c norm ", norm(lag_c) # this is the correct value of 0.816496580928 if inject(psi0, psi0_coarse) is not called # but if inject(psi0, psi0_coarse) is called then # this gives 0.825399087311 which is a relative error of 0.019032990351 print "lag relative error ", np.absolute(norm(lag_c) - norm(lag_f)) / norm(lag_f)
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Dear Wei, thanks, this was a stupid bug I introduced a while ago. I fix it and add a test in this branch: https://github.com/firedrakeproject/firedrake/pull/1079 This will be merged soon, and if you then run firedrake-update things should work. Cheers, Lawrence On 19/05/17 19:02, Pan, Wei wrote:
Dear team,
In the following code, I call /inject/ twice, each time on separate sets of variables {lag_f, lag_c} and {psi0, psi0_c}, where lag_f and psi0 are defined on the fine mesh, and lag_c and psi0_c are defined on the coarse mesh. It seems that the first inject call on {psi0, psi0_c} changes the result of the second inject call on{lag_f, lag_c} which is not what I would expect as psi0 and lag are unrelated. I have highlighted the relevant definitions in blue, and problem bits in red.
I should add that the problem goes away if in the VectorFunctionSpace definitions of vfs_f and vfs_c I change "CG" to "DG". But I'm worried about the result of the CG inject. Many thanks in advance.
Kind regards
Wei
from firedrake import * import numpy as np
L = 1.0 # Zonal length n0 = 50 # Spatial resolution
mesh_base = SquareMesh(n0, n0, L) mesh_hierarchy = MeshHierarchy(mesh_base, 1) mesh = mesh_hierarchy[1] # fine mesh mesh_c = mesh_hierarchy[0] # coarse mesh
Vdg = FunctionSpace(mesh,"DG",1) # DG elements for Potential Vorticity (PV) Vcg = FunctionSpace(mesh,"CG",1) # CG elements for Streamfunction
vfs_f = VectorFunctionSpace(mesh,"CG",1) vfs_c = VectorFunctionSpace(mesh_c,"CG",1) vxf, vyf = SpatialCoordinate(vfs_f)
x = SpatialCoordinate(mesh) q0 = Function(Vdg).interpolate(0.1*sin(x[0])*sin(x[1])) q1 = Function(Vdg) psi0 = Function(Vcg) # Streamfunctions for different time steps
psi = TrialFunction(Vcg) # Test function phi = TestFunction(Vcg) # Trial function
Apsi = (inner(grad(psi),grad(phi)))*dx Lpsi = -q1*phi*dx bc1 = DirichletBC(Vcg, 0., 'on_boundary') psi_problem = LinearVariationalProblem(Apsi,Lpsi,psi0,bcs=bc1) psi_solver = LinearVariationalSolver(psi_problem, solver_parameters={ 'ksp_type':'cg', 'pc_type':'sor' })
gradperp = lambda u: as_vector((-u.dx(1), u.dx(0)))
lag_f = Function(vfs_f).interpolate(as_vector((vxf, vyf)))
if 1: print "-------------------------------------------------------------"
q1.assign(q0) psi_solver.solve() psi0_coarse = Function(FunctionSpace(mesh_c, "CG", 1)) inject(psi0, psi0_coarse)
print "psi0 norm ", norm(psi0) print "psi0_coarse norm ", norm(psi0_coarse) print "psi0_coarse relative error ", np.absolute(norm(psi0_coarse) - norm(psi0)) / norm(psi0) # relative error 0.000583086120876
print "-------------------------------------------------------------"
lag_c = Function(vfs_c) inject(lag_f, lag_c)
print "analytical value ", np.sqrt(2./3.) # 0.816496580928
print "lag0_f norm ", norm(lag_f) print "lag0_c norm ", norm(lag_c) # this is the correct value of 0.816496580928 if inject(psi0, psi0_coarse) is not called # but if inject(psi0, psi0_coarse) is called then # this gives 0.825399087311 which is a relative error of 0.019032990351 print "lag relative error ", np.absolute(norm(lag_c) - norm(lag_f)) / norm(lag_f)
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Hi Lawrence, Awesome! May I ask when will the merge happen? Thanks. Best Wei ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk> Sent: 23 May 2017 16:50:52 To: firedrake Subject: Re: [firedrake] mesh hierarchy inject function issue - possible bug Dear Wei, thanks, this was a stupid bug I introduced a while ago. I fix it and add a test in this branch: https://github.com/firedrakeproject/firedrake/pull/1079 This will be merged soon, and if you then run firedrake-update things should work. Cheers, Lawrence On 19/05/17 19:02, Pan, Wei wrote:
Dear team,
In the following code, I call /inject/ twice, each time on separate sets of variables {lag_f, lag_c} and {psi0, psi0_c}, where lag_f and psi0 are defined on the fine mesh, and lag_c and psi0_c are defined on the coarse mesh. It seems that the first inject call on {psi0, psi0_c} changes the result of the second inject call on{lag_f, lag_c} which is not what I would expect as psi0 and lag are unrelated. I have highlighted the relevant definitions in blue, and problem bits in red.
I should add that the problem goes away if in the VectorFunctionSpace definitions of vfs_f and vfs_c I change "CG" to "DG". But I'm worried about the result of the CG inject. Many thanks in advance.
Kind regards
Wei
from firedrake import * import numpy as np
L = 1.0 # Zonal length n0 = 50 # Spatial resolution
mesh_base = SquareMesh(n0, n0, L) mesh_hierarchy = MeshHierarchy(mesh_base, 1) mesh = mesh_hierarchy[1] # fine mesh mesh_c = mesh_hierarchy[0] # coarse mesh
Vdg = FunctionSpace(mesh,"DG",1) # DG elements for Potential Vorticity (PV) Vcg = FunctionSpace(mesh,"CG",1) # CG elements for Streamfunction
vfs_f = VectorFunctionSpace(mesh,"CG",1) vfs_c = VectorFunctionSpace(mesh_c,"CG",1) vxf, vyf = SpatialCoordinate(vfs_f)
x = SpatialCoordinate(mesh) q0 = Function(Vdg).interpolate(0.1*sin(x[0])*sin(x[1])) q1 = Function(Vdg) psi0 = Function(Vcg) # Streamfunctions for different time steps
psi = TrialFunction(Vcg) # Test function phi = TestFunction(Vcg) # Trial function
Apsi = (inner(grad(psi),grad(phi)))*dx Lpsi = -q1*phi*dx bc1 = DirichletBC(Vcg, 0., 'on_boundary') psi_problem = LinearVariationalProblem(Apsi,Lpsi,psi0,bcs=bc1) psi_solver = LinearVariationalSolver(psi_problem, solver_parameters={ 'ksp_type':'cg', 'pc_type':'sor' })
gradperp = lambda u: as_vector((-u.dx(1), u.dx(0)))
lag_f = Function(vfs_f).interpolate(as_vector((vxf, vyf)))
if 1: print "-------------------------------------------------------------"
q1.assign(q0) psi_solver.solve() psi0_coarse = Function(FunctionSpace(mesh_c, "CG", 1)) inject(psi0, psi0_coarse)
print "psi0 norm ", norm(psi0) print "psi0_coarse norm ", norm(psi0_coarse) print "psi0_coarse relative error ", np.absolute(norm(psi0_coarse) - norm(psi0)) / norm(psi0) # relative error 0.000583086120876
print "-------------------------------------------------------------"
lag_c = Function(vfs_c) inject(lag_f, lag_c)
print "analytical value ", np.sqrt(2./3.) # 0.816496580928
print "lag0_f norm ", norm(lag_f) print "lag0_c norm ", norm(lag_c) # this is the correct value of 0.816496580928 if inject(psi0, psi0_coarse) is not called # but if inject(psi0, psi0_coarse) is called then # this gives 0.825399087311 which is a relative error of 0.019032990351 print "lag relative error ", np.absolute(norm(lag_c) - norm(lag_f)) / norm(lag_f)
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
participants (2)
-
Lawrence Mitchell
-
Pan, Wei