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!