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)