Interpolate into function, and fail in PETSc vector norm (dolfinx)

Linear Algebra
1

Hello,

I’m trying to do a fairly simple thing: Interpolate an expression into a function named “test”, and then take the vector norm via PETSc (in dolfinx from today). The underlying vector has the right entries, however the norm routine always returns the initial value after test.vector.set(0.0), while when taking the norm “manually”, the right result is obtained. Here a MWE:

#!/usr/bin/env python3
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc

import dolfinx
from dolfinx import Function, FunctionSpace, BoxMesh, cpp, solve
from dolfinx.cpp.mesh import CellType
from ufl import TrialFunction, TestFunction, FiniteElement, inner, dx, sqrt

mesh = BoxMesh(
    MPI.COMM_WORLD, [np.array([0.0, 0.0, 0.0]),
                     np.array([2.0, 1.0, 1.0])], [2, 2, 2],
    CellType.tetrahedron, dolfinx.cpp.mesh.GhostMode.none)

W0e = FiniteElement("Quadrature", mesh.ufl_cell(), degree=1, quad_scheme='default')
W0 = FunctionSpace(mesh, W0e)

dx = dx(metadata={'quadrature_degree': 1})

def project(v, V):

    w = TestFunction(V)
    Pv = TrialFunction(V)
    a = inner(w, Pv) * dx
    L = inner(w, v) * dx
    # solve linear system for projection
    function = Function(V)
    
    solve(a==L, function)
    
    return function


# a function
test = Function(W0, name="test")
test.vector.set(0.0)

expr = 1.0

test_proj = project(expr, W0)
test.interpolate(test_proj)

# with this line (that does nothing), both norms are equal - without, the PETSc norm is zero!
#test.vector[0] += 0

norm_test2 = 0
for i in range(len(test.vector[:])):
    norm_test2 += test.vector[i]**2.
    
norm_test = sqrt(norm_test2)

print("Norm           : %f" % norm_test)
print("Norm from PETSc: %f" % test.vector.norm())

The two print routines should yield the same (6.928203). However, if I insert the commented line test.vector[0] += 0 (which adds nothing on the first entry of the vector), then the norm function returns the correct value.

This seems very weird to me. Does anyone have an idea what happens? A bug in interpolate, or rather in PETSc?

Bw,
Marc

2

Add the line:

test.vector.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)

or

test.vector.assemble()

after interpolation.

3

Oh, thanks! That worked! I also tried the line ghostUpdate with ScatterMode Forward, but that didn’t work. The Reverse option does the job!