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
>