Interpolate Constant as Initial Condition in mixed element space

Errors
1

I’m trying to set all nodes to the same initial condition, which is a constant.
However, I cannot figure out why I keep getting TypeError: interpolate(): incompatible function arguments. How exactly should I interpolate a constant to all nodes?

The same error persists even without two function spaces (which I created as only 2 of 5 variables are time-dependent. Would appreciate any advice on whether this is a wise option or whether I should just put dummy variables in place too!).

Dolfinx installed via conda on Mac M1 Monterey.

Thank you!

Below is the code:

import ufl
from dolfinx import mesh, fem, io, nls, log
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np

xmin = -5
xmax = 5
ymax = 5
domain = mesh.create_rectangle(MPI.COMM_WORLD, [(xmin,0),(xmax,ymax)], n=[200, 50])

x = ufl.SpatialCoordinate(domain)
degree = 3
P1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), degree) 
element = ufl.MixedElement([P1, P1, P1, P1, P1]) #create a mixed function space 
V = fem.FunctionSpace(domain, element)  

#we create another function space for the t dependent h1 and h2 only
element2 = ufl.MixedElement([P1, P1])
W = fem.FunctionSpace(domain, element2) 

v1, v2, v3, v4, v5 = ufl.TestFunctions(V)

u = fem.Function(V) 
u_n = fem.Function(W) #update previous solution 
h1, h2, p1, p2, ul = ufl.split(u) 
h1_n, h2_n = ufl.split(u_n) 

#Define constants
Hbath = 2.0

#ERROR
u.sub(0).interpolate(fem.Constant(domain, PETSc.ScalarType(Hbath)), V)
u.x.scatter_forward()
2 
TypeError: interpolate(): incompatible function arguments. The following argument types are supported:
    1. (self: dolfinx.cpp.fem.Function_float64, f: Callable[[numpy.ndarray[numpy.float64]], numpy.ndarray[numpy.float64]], cells: numpy.ndarray[numpy.int32]) -> None
    2. (self: dolfinx.cpp.fem.Function_float64, u: dolfinx.cpp.fem.Function_float64, cells: numpy.ndarray[numpy.int32]) -> None
    3. (self: dolfinx.cpp.fem.Function_float64, arg0: dolfinx::fem::Expression<double>, arg1: numpy.ndarray[numpy.int32]) -> None

As the error message shows, interpolate expects:

  1. a positional argument that is either a callable, Expression, or Function
  2. an optional keyword argument that specifies a list of cells over which to interpolate

You can also see the expected call signature in the source code:

To interpolate a constant, you can define a simple lambda function that returns a constant value, e.g.:

import ufl
from dolfinx import mesh, fem, io, nls, log
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np

xmin = -5
xmax = 5
ymax = 5
domain = mesh.create_rectangle(MPI.COMM_WORLD, [(xmin,0),(xmax,ymax)], n=[200, 50])

x = ufl.SpatialCoordinate(domain)
degree = 3
P1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), degree) 
element = ufl.MixedElement([P1, P1, P1, P1, P1]) #create a mixed function space 
V = fem.FunctionSpace(domain, element)  

#we create another function space for the t dependent h1 and h2 only
element2 = ufl.MixedElement([P1, P1])
W = fem.FunctionSpace(domain, element2) 

v1, v2, v3, v4, v5 = ufl.TestFunctions(V)

u = fem.Function(V) 
u_n = fem.Function(W) #update previous solution 
h1, h2, p1, p2, ul = ufl.split(u) 
h1_n, h2_n = ufl.split(u_n) 

#Define constants
Hbath = 2.0

#ERROR
u.sub(0).interpolate(lambda x: np.full((x.shape[1],), Hbath))
u.x.scatter_forward()
1 Like
3

Thanks for the quick response!
Just wanted to point out to anyone reading this in the future that it was the initial definition of Hbath to a fem.Constant that would throw a later error "Couldn’t map c_2 to a float, returning ufl object without evaluation’.

4

a bit late, but is this really needed?

5

If you want to push forward ghost values, yes.

6

If that is the case, shouldn’t the arrays be different before and after scatter forward, e.g. in the following example? If not, what would be an example to show the effect? I’m trying to understand which operations require forward/reverse scattering of the ghost values (if there is some document explaining this, please feel free to share). Maybe this would be better placed in a separate topic.

import dolfinx
from mpi4py import MPI

N = 4
degree = 1
comm = MPI.COMM_WORLD
rank = comm.rank
mesh = dolfinx.mesh.create_unit_square(comm, N, N)
V = dolfinx.fem.FunctionSpace(mesh, ("CG", degree))

u = dolfinx.fem.Function(V)
u.interpolate(lambda x: x[0] * x[1])

print(f"{rank = }: before:\n\t {u.x.array = }")
u.x.scatter_forward()
# u.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
print("SCATTER FORWARD")
print(f"{rank = }: after\n\t {u.x.array = }")

Also while I have your attention: is there any difference between u.x.scatter_forward() and u.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD) (apart from the latter being implemented in PETSc)?

7

Interpolation (when no list of cell indices is supplied) works on all cells local to the process, and will therefore update all dofs local to the process (included ghosted dofs).

If you supply a subset of cell indices (local index) you would need to use scatter forward to ensure that all data has been updated consistently.

No difference except x is the DOLFINx implementation (using MPI neighbourhood communicators), while the petsc one uses their own mpi communication system.

2 Likes