Integration measures for coupled bulk-surface equations

Hi all!

I am trying to implement a solver for a system of coupled bulk-surface equations with FEniCSx version 0.8.0. The variational form of the system is given by

\begin{cases} \displaystyle \delta_\Omega \int_\Omega u_t \varphi \, \mathrm{d}V + \int_\Omega \vec{\nabla} u \cdot \vec{\nabla} \varphi \, \mathrm{d}V = -\frac{1}{\delta_k} \int_{\Gamma} u w \varphi \, \mathrm{d}S\text{,}\\ \displaystyle \int_\Gamma w_t \psi \, \mathrm{d}S + \delta_\Gamma \int_\Gamma w_x \psi' \, \mathrm{d}S = - \frac{1}{\delta_k} \int_\Gamma u w \psi \, \mathrm{d}S\text{,}\\ \end{cases}

where the unknown function u is defined on \Omega \in \mathbb{R}^n, and the unknown function w on \Gamma \subset \partial \Omega. The constants \delta_\Omega, \delta_\Gamma, \delta_k are given. I used the backward Euler approximation for the time derivatives and employed a two-step evolution to arrive at a fully discretised scheme, but I am unsure how to define the correct integration measures for each variational form. Here is a MWE (or I guess a minimal non-working example) of what I initially tried:

from mpi4py import MPI
from dolfinx.mesh import create_unit_square, create_submesh, locate_entities
from dolfinx.fem import functionspace, Function
from ufl import TestFunction, TrialFunction, dx, ds, dot, nabla_grad
import numpy as np
from pathlib import Path
from dolfinx.io import VTXWriter

delta_omega = 1
delta_gamma = 1
delta_k = 1

t = 0
T = 10
num_steps = 500
dt = T / num_steps


def on_gamma(x):
    return np.isclose(x[1], 0)


bulk_mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
surface_entities = locate_entities(bulk_mesh, bulk_mesh.topology.dim - 1, on_gamma)
surface_mesh, entity_map, vertex_map, geom_map = create_submesh(bulk_mesh, bulk_mesh.topology.dim - 1, surface_entities)

Vb = functionspace(bulk_mesh, ("Lagrange", 1))
Vs = functionspace(surface_mesh, ("Lagrange", 1))


# Create initial conditions
def u0(x):
    return x[1]


def w0(x):
    return np.square(np.sin(np.pi * x[0]))


u_n = Function(Vb)
u_n.name = "u_n"
u_n.interpolate(u0)
w_n = Function(Vs)
w_n.name = "w_n"
w_n.interpolate(w0)

uh = Function(Vb)
uh.name = "uh"
uh.interpolate(u0)
wh = Function(Vs)
wh.name = "uh"
wh.interpolate(w0)

results_folder = Path("results")
results_folder.mkdir(exist_ok=True, parents=True)
filename_u = results_folder / "bulk-surface_u"
vtx_writer_u = VTXWriter(bulk_mesh.comm, filename_u.with_suffix(".bp"), [uh], engine="BP4")
vtx_writer_u.write(0.0)
filename_w = results_folder / "bulk-surface_w"
vtx_writer_w = VTXWriter(surface_mesh.comm, filename_w.with_suffix(".bp"), [wh], engine="BP4")
vtx_writer_w.write(0.0)

u = TrialFunction(Vb)
phi = TestFunction(Vb)
w = TrialFunction(Vs)
psi = TestFunction(Vs)

a_bulk = delta_omega * u * phi * dx + dt * dot(nabla_grad(u), nabla_grad(phi)) * dx + dt / delta_k * u * w_n * phi * ds
L_bulk = delta_omega * u_n * phi * dx
a_surface = w * psi * dx + delta_gamma * dt * dot(nabla_grad(w), nabla_grad(psi)) * dx + dt / delta_k * u_n * w * phi * dx
L_surface = w_n * psi * dx

This raises the error:

Traceback (most recent call last):
File "/home/yowie/anaconda3/envs/fenicsx-env/Test.py", line 68, in <module>
    a_bulk = delta_omega * u * phi * dx + dt * dot(nabla_grad(u), nabla_grad(phi)) * dx + dt / delta_k * u * w_n * phi * ds
                                                                                        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~
File "/home/yowie/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ufl/measure.py", line 424, in __rmul__
    raise ValueError(
ValueError: Multiple domains found, making the choice of integration domain ambiguous.

I’m guessing that the problem lies in the fact that u and w_n (in the bulk bilinear form) and u_n and w (in the surface bilinear form) are from different function spaces. My best guess is that I would have to interpolate them somehow. But how can I interpolate w_n into Vb? Maybe I’m not even on the right track with these guesses though, I’m quite new to FEniCS. Any help would be greatly appreciated!

See for instance dolfinx/python/demo/demo_stokes.py at v0.8.0 · FEniCS/dolfinx · GitHub for two different options. What you are trying to do (sum equations with different test functions) isn’t allowed if you use two independent function spaces.

In future posts it would be appreciated if you could provide a clearer description of what was wrong with the code. For instance, I guess you obtained an error at some point (maybe when defining a_bulk already, or, if not, surely when trying to assemble the matrix or solve the system): in that case, post the error, since it usually provides helpful information for people who are willing to answer your question.

Thank for the quick reply! I’ve been staring at the lid-driven cavity flow demo for a while now, but I can’t pinpoint the “two different options” you mention. Could you point me in the right direction? I also looked at the “Combining Dirichlet and Neumann conditions” tutorial, but it doesn’t quite tell me what I need to know.

For future reference, I’ve added the error message to the original post. I’ll make sure to include it in the future, thanks for pointing it out!

Sorry, by two different options I meant nest vs block.

I think I’m starting to understand what’s going on in the demo you linked. But I’m not sure how to deal with the integrals of u w \varphi and u w \psi. The demo only deals with sums of different trial functions, not products. What would be the right way of dealing with this integral in FEniCSx?

it doesn’t matter much. For instance, have a look at tutorial_navier_stokes, function run_block, definition of F.

I can’t make much sense of this demo without a strong or weak formulation of the problem it’s solving. It doesn’t seem to document it anywhere. Do you happen to know the specifics or where to find them?

It’s a Navier-Stokes problem with a fluid that is moving from the left to the right. Feel free to search yourself for the strong form of the equation, and derive the corresponding weak form.

Seems like it’s the same one used in this tutorial. I just can’t see how any of this helps me. Could you give me some direction on how to split the variational form?

Indeed. In

    F = [(nu * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx + ufl.inner(ufl.grad(u) * u, v) * ufl.dx
          - ufl.inner(p, ufl.div(v)) * ufl.dx),
         ufl.inner(ufl.div(u), q) * ufl.dx]

the first entry F[0] corresponds to the weak form of the first equation in the system (42) at the page you linked, while F[1] is the second equation. Note that there is a quadratic term grad(u) * u, which plays the same role as your quadratic terms u w (it doesn’t matter that much that one quadratic term has both unknowns, and the other only one of the two)…

I’ve installed multifenicsx and am taking some time to familiarise myself with the Navier-Stokes demo. I also noticed the demo on subdomain restrictions. I’ll get back as soon as I have a working example (in which case I’ll post the solution) or run into another problem related to the variational form. Thanks for all the help so far!

If we can use both methods—nested matrices with an iterative solver and a block iterative solver — is there a way to determine which method is best for our problem? What criteria should we consider?

Thanks

I’ve spent a considerable amount of time looking into multiphenicsx and the various tutorials provided by its creator. I think I have a better understanding of what needs to be done now, but I’m still running into some trouble. Following the multiphenicsx Navier/Stokes tutorial, I whipped up this code:

from mpi4py import MPI
from dolfinx.mesh import create_unit_square, create_submesh, locate_entities, meshtags
from dolfinx.fem import functionspace, Function, form
from ufl import TestFunction, dx, dot, nabla_grad, Measure
import numpy as np


delta_omega = 10
delta_gamma = 10
delta_k = 0.05

T = 10
num_steps = 500
dt = T / num_steps


def on_gamma(x):
    return np.isclose(x[1], 0)


bulk_mesh = create_unit_square(MPI.COMM_WORLD, 100, 100)
surface_facets = locate_entities(bulk_mesh, bulk_mesh.topology.dim - 1, on_gamma)
facet_tag = meshtags(bulk_mesh, bulk_mesh.topology.dim - 1, surface_facets, 1)

surface_mesh, entity_map, vertex_map, geom_map = create_submesh(bulk_mesh, bulk_mesh.topology.dim - 1, surface_facets)

Vb = functionspace(bulk_mesh, ("Lagrange", 1))
Vs = functionspace(surface_mesh, ("Lagrange", 1))


# Create initial conditions
def u0(x):
    return 2 * x[1]


def w0(x):
    return np.square(np.sin(np.pi * x[0])) + 1


u_n = Function(Vb)
u_n.name = "u_n"
u_n.interpolate(u0)

w_n = Function(Vs)
w_n.name = "w_n"
w_n.interpolate(w0)

u, w = Function(Vb), Function(Vs)
phi, psi = TestFunction(Vb), TestFunction(Vs)

dS = Measure("dS", domain=bulk_mesh, subdomain_data=facet_tag)
ds = Measure("ds", domain=surface_mesh)

F = [delta_omega * u * phi * dx + dt * dot(nabla_grad(u), nabla_grad(phi)) * dx + dt / delta_k * u * w * phi("+") * dS(1) - delta_omega * u_n * phi * dx,
        w * psi * dx + delta_gamma * dt * dot(nabla_grad(w), nabla_grad(psi)) * dx + dt / delta_k * u * w * psi * ds - w_n * psi * dx]

index_map = bulk_mesh.topology.index_map(bulk_mesh.topology.dim)
inv_entity_map = np.full(index_map.size_local + index_map.num_ghosts, -1)
inv_entity_map[entity_map] = np.arange(len(entity_map))
entity_maps = {surface_mesh: inv_entity_map}

form(F, entity_maps=entity_maps)

However, I get the error

Traceback (most recent call last):
  File "/home/yowie/anaconda3/envs/fenicsx-env/Test.py", line 62, in <module>
	form(F, entity_maps=entity_maps)
  File "/home/yowie/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 249, in form
	return _create_form(form)
		   ^^^^^^^^^^^^^^^^^^
  File "/home/yowie/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 246, in _create_form
	return list(map(lambda sub_form: _create_form(sub_form), form))
		   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/yowie/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 246, in <lambda>
	return list(map(lambda sub_form: _create_form(sub_form), form))
									 ^^^^^^^^^^^^^^^^^^^^^^
  File "/home/yowie/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 244, in _create_form
	return _form(form)
		   ^^^^^^^^^^^
  File "/home/yowie/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 186, in _form
	ufcx_form, module, code = jit.ffcx_jit(
							  ^^^^^^^^^^^^^
  File "/home/yowie/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/jit.py", line 51, in mpi_jit
	return local_jit(*args, **kwargs)
		   ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/yowie/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/jit.py", line 201, in ffcx_jit
	r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
		^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/yowie/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ffcx/codegeneration/jit.py", line 276, in compile_forms
	raise e
  File "/home/yowie/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ffcx/codegeneration/jit.py", line 256, in compile_forms
	impl = _compile_objects(
		   ^^^^^^^^^^^^^^^^^
  File "/home/yowie/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ffcx/codegeneration/jit.py", line 383, in _compile_objects
	_, code_body = ffcx.compiler.compile_ufl_objects(
				   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/yowie/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ffcx/compiler.py", line 113, in compile_ufl_objects
	ir = compute_ir(analysis, _object_names, _prefix, options, visualise)
		 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/yowie/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ffcx/ir/representation.py", line 225, in compute_ir
	_compute_integral_ir(
  File "/home/yowie/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ffcx/ir/representation.py", line 588, in _compute_integral_ir
	integral_ir = compute_integral_ir(
				  ^^^^^^^^^^^^^^^^^^^^
  File "/home/yowie/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ffcx/ir/integral.py", line 85, in compute_integral_ir
	mt_table_reference = build_optimized_tables(
						 ^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/yowie/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ffcx/ir/elementtables.py", line 360, in build_optimized_tables
	get_ffcx_table_values(
  File "/home/yowie/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ffcx/ir/elementtables.py", line 139, in get_ffcx_table_values
	tbl = component_element.tabulate(deriv_order, entity_points)
		  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/yowie/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/basix/ufl.py", line 660, in tabulate
	tables = self._element.tabulate(nderivs, points)
			 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/yowie/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/basix/ufl.py", line 415, in tabulate
	tab = self._element.tabulate(nderivs, points)
		  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/yowie/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/basix/finite_element.py", line 137, in tabulate
	return self._e.tabulate(n, x)
		   ^^^^^^^^^^^^^^^^^^^^^^
RuntimeError: Point dim (2) does not match element dim (1).

I got this error even before I included the entity_maps argument in form, and a quick search on how to solve it told me to add the entity_maps argument. However, did this not fix the error. What am I doing wrong here?

Here you are trying to do a mixed variational form with a codim-1 mesh, which is currently under development in the following branches:

With these two PRs you should be able to create a mixed variational form with variables living in cells, and variables that only live in a function space on a subset of facets.

Thanks Jørgen, this seems to be exactly what I’m looking for! If any further questions pop up while I dive into the mixed-dimensional PRs, I’ll open a new thread.