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)