Restricting equation terms to a subdomain

General
1

After some searching, I have not found a way to do the following in dolfinx (I believe it can be done in FeniCS with Expression): I am considering reaction-diffusion in 3D where the reaction term is only applied on a simple internal subdomain of a larger (e.g., cubic) domain – for instance, a small sphere, with center (d, 0, 0). How do I restrict the reaction term in this way? Here is a MWE based on some posts I have seen. It returns the error TypeError: unsupported operand type(s) for +: ‘int’ and ‘NoneType’ on the last line. The weak form is the same as for the heat equation demo, just changing f. Sorry if the issue is very obvious!

 # Parameters for reaction-diffusion
D = 800 # µm^2/s
kp = 57
Km = 0.001
d = 1

def initial_condition(x):
    return 1000.0 * (x[0]**2 + x[1]**2 + x[2]**2 <= 0.250)

nx, ny, nz = 50, 50, 50
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([-10,-10,-10]), np.array([10, 10, 10])], 
                               [nx, ny, nz], mesh.CellType.tetrahedron)

V = fem.FunctionSpace(domain, ("CG", 2))

u_n = fem.Function(V)
u_n.name = "u_n"
u_n.interpolate(initial_condition)

fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(
    domain, fdim, lambda x: np.full(x.shape[1], True, dtype=bool))
bc = fem.dirichletbc(PETSc.ScalarType(0), fem.locate_dofs_topological(V, fdim, boundary_facets), V)

uh = fem.Function(V)
uh.name = "uh"
uh.interpolate(initial_condition)

def E(x):
    return 1.0 * ((x[0]-d)**2 + x[1]**2 + x[2]**2 <= 0.250)

def E_inv(x):
    return 1.0 * ((x[0]-d)**2 + x[1]**2 + x[2]**2 > 0.250)

x = ufl.SpatialCoordinate(domain)
v = ufl.TestFunction(V)
f = - E(x) * kp * uh / Km / (1 + uh / Km) - E_inv(x) * kp * rho * uh / Km / (1 + uh / Km)
F = uh * v * ufl.dx + dt*ufl.dot(ufl.grad(uh), ufl.grad(v)) * ufl.dx - (u_n + dt * f) * v * ufl.dx

problem = fem.petsc.NonlinearProblem(F, uh, bcs=[bc])
2

You can either:

  1. Use ufl.conditional combined with x = ufl.SpatialCoordinate(domain) to generate your conditional expressions in proper ufl syntax
  2. Interpolate the expressions E and E_inv into a suitable function space.

For further questions, please make sure that all imports and parameters are defined, such that one can re-run the code easily.

3

Okay, got it – thanks!