Hi @dokken, hope you’re doing well. I tried changing all 3 lines of the formulation in @RemDelaporteMathurin’s MWE that rely on u_sub/v_sub to use ufl.dx(domain=submesh) instead of ds(1), as I believe you were indicating in your diagnosis.
Full Modified MWE
from mpi4py import MPI
from petsc4py import PETSc
import dolfinx
import dolfinx.fem.petsc
import matplotlib.pyplot as plt
import numpy as np
import pyvista
import ufl
from dolfinx import plot
dx = 1 / 5
L = 100
nx = int(L / dx)
mesh = dolfinx.mesh.create_rectangle(
MPI.COMM_WORLD,
[np.array([0, 0]), np.array([L, 1])],
[nx, 10],
cell_type=dolfinx.mesh.CellType.quadrilateral,
)
vdim = mesh.topology.dim
fdim = mesh.topology.dim - 1
mesh.topology.create_connectivity(fdim, vdim)
# facet meshtags top and bottom
tag_to_marker = {
1: lambda x: np.isclose(x[1], 0), # bottom
2: lambda x: np.isclose(x[1], 1), # top
3: lambda x: np.isclose(x[0], 0), # left
4: lambda x: np.isclose(x[0], L), # right
}
facets = np.array([], dtype=np.int64)
tags = np.array([], dtype=np.int32)
for tag, marker in tag_to_marker.items():
facet_indices = dolfinx.mesh.locate_entities(mesh, fdim, marker)
facets = np.concatenate((facets, facet_indices))
tags = np.concatenate((tags, np.full_like(facet_indices, tag, dtype=np.int32)))
facet_tags = dolfinx.mesh.meshtags(mesh, fdim, facets, tags)
cell_tags = dolfinx.mesh.meshtags(
mesh,
vdim,
np.arange(mesh.topology.index_map(vdim).size_local),
np.ones(mesh.topology.index_map(vdim).size_local, dtype=np.int32),
)
with dolfinx.io.XDMFFile(mesh.comm, "results/facet_tags.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_meshtags(facet_tags, x=mesh.geometry)
with dolfinx.io.XDMFFile(mesh.comm, "results/cell_tags.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_meshtags(cell_tags, x=mesh.geometry)
# make submesh of the bottom boundary
submesh, cmap, vmap, nmap = dolfinx.mesh.create_submesh(
mesh, dim=fdim, entities=facet_tags.find(1)
)
submesh.topology.create_connectivity(0, 1)
# Function spaces and functions
V_bulk = dolfinx.fem.functionspace(mesh, ("CG", 1))
V_sub = dolfinx.fem.functionspace(submesh, ("CG", 1))
W = ufl.MixedFunctionSpace(V_bulk, V_sub)
u = dolfinx.fem.Function(V_bulk)
u.name = "u"
u_sub = dolfinx.fem.Function(V_sub)
u_sub.name = "u_sub"
v, v_sub = ufl.TestFunctions(W)
# Formulation
dx = ufl.dx(domain=mesh, subdomain_data=cell_tags)
ds = ufl.ds(domain=mesh, subdomain_data=facet_tags)
dx_sub = ufl.dx(domain=submesh)
F = ufl.inner(ufl.grad(u), ufl.grad(v)) * dx
F += ufl.inner(ufl.grad(u_sub), ufl.grad(v_sub)) * dx_sub # <-- CHANGED
vel_x = 10
# Option 1: Full grad with 2D vector. Works but odd that we need a 2D velocity
vel = dolfinx.fem.Constant(submesh, PETSc.ScalarType([vel_x, vel_x]))
F += ufl.inner(ufl.dot(ufl.grad(u_sub), vel), v_sub) * dx_sub # <-- CHANGED
# coupling term
h_l = dolfinx.fem.Constant(mesh, 0.4)
flux = h_l * (u - u_sub)
F += flux * v * ds(1)
F += -flux * v_sub * dx_sub # <-- CHANGED
forms = ufl.extract_blocks(F)
# Dirichlet BC left
bc_top_dofs = dolfinx.fem.locate_dofs_topological(
V_bulk,
mesh.topology.dim - 1,
dolfinx.mesh.locate_entities(
mesh, mesh.topology.dim - 1, lambda x: np.isclose(x[1], 1)
),
)
bc_top = dolfinx.fem.dirichletbc(
dolfinx.default_scalar_type(0.0),
bc_top_dofs,
V_bulk,
)
bc_left_dofs = dolfinx.fem.locate_dofs_topological(
V_sub, 0, dolfinx.mesh.locate_entities(submesh, 0, lambda x: np.isclose(x[0], 0))
)
bc_left = dolfinx.fem.dirichletbc(
dolfinx.default_scalar_type(1.0),
bc_left_dofs,
V_sub,
)
# Nonlinear problem
problem = dolfinx.fem.petsc.NonlinearProblem(
forms,
[u, u_sub],
bcs=[
bc_top,
bc_left,
],
petsc_options_prefix="codim1_prob",
entity_maps=[cmap],
)
problem.solve()
# Post processing
with dolfinx.io.VTXWriter(mesh.comm, "results/u.bp", [u]) as writer:
writer.write(0.0)
with dolfinx.io.VTXWriter(submesh.comm, "results/u_sub.bp", [u_sub]) as writer:
writer.write(0.0)
topology, cell_types, geometry = plot.vtk_mesh(u.function_space)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
grid.point_data["c"] = u.x.array
grid.set_active_scalars("c")
plotter = pyvista.Plotter()
plotter.add_mesh(grid)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
plotter.show()
else:
figure = plotter.screenshot("u.png")
topology, cell_types, geometry = plot.vtk_mesh(u_sub.function_space)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
grid.point_data["c"] = u_sub.x.array
grid.set_active_scalars("c")
# Make two points to construct the line between
a = [0, 0, 0]
b = [L, 0, 0]
sample = grid.sample_over_line(a, b, resolution=100)
plt.plot(sample["Distance"], sample["c"])
plt.ylim(0, 1)
plt.xlabel("x")
plt.ylabel("u_sub")
plt.show()
File ".../lib/python3.11/site-packages/ffcx/ir/elementtables.py", line 527, in build_optimized_tables
tbl = clamp_table_small_numbers(t["array"], rtol=rtol, atol=atol)
^
UnboundLocalError: cannot access local variable 't' where it is not associated with a value
Full traceback
Traceback (most recent call last):
File ".../MWE_revised.py", line 132, in <module>
problem = dolfinx.fem.petsc.NonlinearProblem(
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File ".../lib/python3.11/site-packages/dolfinx/fem/petsc.py", line 1239, in __init__
self._F = _create_form(
^^^^^^^^^^^^^
File ".../lib/python3.11/site-packages/dolfinx/fem/forms.py", line 449, in form
return _create_form(form)
^^^^^^^^^^^^^^^^^^
File ".../lib/python3.11/site-packages/dolfinx/fem/forms.py", line 445, in _create_form
return list(map(lambda sub_form: _create_form(sub_form), form))
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File ".../lib/python3.11/site-packages/dolfinx/fem/forms.py", line 445, in <lambda>
return list(map(lambda sub_form: _create_form(sub_form), form))
^^^^^^^^^^^^^^^^^^^^^^
File ".../lib/python3.11/site-packages/dolfinx/fem/forms.py", line 441, in _create_form
return _form(form)
^^^^^^^^^^^
File ".../lib/python3.11/site-packages/dolfinx/fem/forms.py", line 361, in _form
ufcx_form, module, code = jit.ffcx_jit(
^^^^^^^^^^^^^
File ".../lib/python3.11/site-packages/dolfinx/jit.py", line 60, in mpi_jit
return local_jit(*args, **kwargs)
^^^^^^^^^^^^^^^^^^^^^^^^^^
File "..../lib/python3.11/site-packages/dolfinx/jit.py", line 215, in ffcx_jit
r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File ".../lib/python3.11/site-packages/ffcx/codegeneration/jit.py", line 244, in compile_forms
raise e
File ".../lib/python3.11/site-packages/ffcx/codegeneration/jit.py", line 224, in compile_forms
impl = _compile_objects(
^^^^^^^^^^^^^^^^^
File ".../lib/python3.11/site-packages/ffcx/codegeneration/jit.py", line 349, in _compile_objects
_, code_body = ffcx.compiler.compile_ufl_objects(
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File ".../lib/python3.11/site-packages/ffcx/compiler.py", line 113, in compile_ufl_objects
ir = compute_ir(analysis, _object_names, _prefix, options, visualise)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File ".../lib/python3.11/site-packages/ffcx/ir/representation.py", line 150, in compute_ir
irs = [
^
File ".../python3.11/site-packages/ffcx/ir/representation.py", line 151, in <listcomp>
_compute_integral_ir(
File ".../lib/python3.11/site-packages/ffcx/ir/representation.py", line 402, in _compute_integral_ir
integral_ir = compute_integral_ir(
^^^^^^^^^^^^^^^^^^^^
File ".../lib/python3.11/site-packages/ffcx/ir/integral.py", line 170, in compute_integral_ir
mt_table_reference = build_optimized_tables(
^^^^^^^^^^^^^^^^^^^^^^^
File ".../lib/python3.11/site-packages/ffcx/ir/elementtables.py", line 527, in build_optimized_tables
tbl = clamp_table_small_numbers(t["array"], rtol=rtol, atol=atol)
^
UnboundLocalError: cannot access local variable 't' where it is not associated with a value
Exception ignored in: <function NonlinearProblem.__del__ at 0x77a1def66520>
Traceback (most recent call last):
File ".../lib/python3.11/site-packages/dolfinx/fem/petsc.py", line 1361, in __del__
lambda obj: obj is not None, (self._snes, self._A, self._b, self._x, self._P_mat)
^^^^^^^^^^
AttributeError: 'NonlinearProblem' object has no attribute '_snes'
Looking at the full traceback appears to indicate some sort of form compilation error with UFL. Are my changes not the intended fix?
Thank you for any clarification!