Function from array for rhs gives wrong result

[luisfpereira]

We (@JosquinHa and I) are trying to solve the following variational problem:

a = inner(v, u)*dx - ((
    -inner(Dx(Dx(Dx(v, 0), 0), 0), Dx(Dx(Dx(u, 0), 0), 0))
    - 3 * inner(Dx(Dx(Dx(v, 1), 0), 0), Dx(Dx(Dx(u, 1), 0), 0))
    - 3 * inner(Dx(Dx(Dx(v, 0), 1), 1), Dx(Dx(Dx(u, 0), 1), 1))
    - inner(Dx(Dx(Dx(v, 1), 1), 1), Dx(Dx(Dx(u, 1), 1), 1))
)/10)*dx

L = inner(chi, u)*dx

where chi is given by an array (representing a shape), and u belongs to "WXRobH3NC" function space of degree 7 (following some last year discussions with @andreasbock). In particular, we do:

mesh = PeriodicUnitSquareMesh(nx=n, ny=n)
WX = FunctionSpace(mesh, "WXRobH3NC", degree=7)

V = FunctionSpace(mesh, "CG", 1)

m = V.ufl_domain()
W = VectorFunctionSpace(m, V.ufl_element())
v, u = TrialFunction(WX), TestFunction(WX)

(Notice this element is not available in master yet, as far as I'm aware. Our docker setup can be found here).

If we treat chi as an interpolated function given by an expression, i.e.

x, y = SpatialCoordinate(mesh)
expr = cos(2*pi*x) + sin(2*pi*y)
chi = interpolate(expr, V)

we obtain the expected answer sol = -(cos(2*pi*x) + sin(2*pi*y)) / ((2*pi)**6).

On the other hand, if we consider chi to come from an array, i.e.

chi = Function(V)
chi.dat.data[:] = data  # where data is a 1d vector with properly ordered values

we get a wrong solution (i.e. not close at all, both quantitatively and qualitatively, to the expected result).

We solve using defaults:

res = Function(WX, name="field")
solve(a == L, res)

Does anyone have an idea on what we're doing wrong? I believe the problem may be easy to identify, if not, I can quickly provide more information on how to get data from the expression to try to emulate the known case.

[ksagiyam]

Could you print(chi.dat.data_ro) for both cases and check they are indeed identical?

[connorjward]

Assuming that this is being run in parallel what happens if you copy over data_ro_with_halos to data_with_halos?

Another thing to look at might be to check whether assemble(L) yields the same thing in both cases.

[djoharr]

Passing data_ro_with_halos does not change anything. Also, data, data_with_halos, data_ro and data_ro_with_halos are all exactly the same (for both chi).
When comparing assemble(L_1).dat.data and assemble(L_2).dat.data, I get exactly the same result. However, when printing assemble(L) in both case I don't get quite the same values:

• Case with expression: Coefficient(WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x7f30207f2da0>, FiniteElement('WuXu Robust H3 Nonconforming', triangle, 7), name=None), Mesh(VectorElement(FiniteElement('Discontinuous Lagrange', triangle, 1, variant='equispaced'), dim=2, variant='equispaced'), 73)), 179)
• Case with array: Coefficient(WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x7f3020649600>, FiniteElement('WuXu Robust H3 Nonconforming', triangle, 7), name=None), Mesh(VectorElement(FiniteElement('Discontinuous Lagrange', triangle, 1, variant='equispaced'), dim=2, variant='equispaced'), 60)), 181)
  I hope any of this is relevant.

[connorjward]

The differences you see there are trivial. The Coefficients simply have different count attributes (kind of like an ID). It just means that they are distinct objects.
Would you be able to provide a more minimal example that we could run to investigate the behaviour? For instance do you see the same behaviour with a simpler element, or with a simpler form?

[djoharr]

Hey,
I apologize for the wrong answers I gave before, it seems all I had to do was switching from python notebook to python script. When passing the data from the expression as array things are now working. It seems the notebook I was working on did not update well...
I have corrected the way of ordering the data array and everything works as intended.
Many thanks for all the instructive replies.

[connorjward]

Glad to hear it's sorted.