Solving a system of coupled PDEs with different dimensions

dolfinx
, ,
1

Hello,

I have a question regarding this tutorial:

Coupling PDEs of multiple dimensions

I have used the same approach outlined in this tutorial to address my problem of solving partial differential equations (PDEs) in both 2D and 1D domains and coupling them together. Specifically, I need to solve an equation that closely resembles the Poisson equation on a 1D domain, which serves as one of the boundaries of the 2D domain. Here is what I have:

∇2ψ=f(x)

Then, in variational form, I have:

ufl.inner(ufl.grad(psi), ufl.grad(w)) * ds

I attempted to follow the tutorial mentioned earlier by adding the above terms into the variational form. However, I am encountering an error when I try to form and extract the residuals.

Code:

F = alpha * ufl.inner(sigma(u, mu, lmbda), epsilon(v)) * dx
F -= alpha * ufl.inner(f, v) * dx
F += -ufl.inner(psi - psi_k, ufl.dot(v, n)) * ds
F += ufl.inner(ufl.dot(u, n_g), w) * ds
F += ufl.inner(ufl.exp(psi), w) * ds - ufl.inner(g, w) * ds
F += ufl.inner(ufl.grad(psi), ufl.grad(w)) * ds # Addiotinal term for Poisson equation
residual = dolfinx.fem.form(ufl.extract_blocks(F), entity_maps=entity_maps)

Error when running in debug mode on the last line of the previous code

Exception has occurred: VerificationError
CompileError: command '/usr/bin/gcc' failed with exit code 1
subprocess.CalledProcessError: Command '['/usr/bin/gcc', '-pthread', '-B', '/home/phavaej/anaconda3/envs/fenicsx-env/compiler_compat', '-fno-strict-overflow', '-Wsign-compare', '-DNDEBUG', '-O2', '-Wall', '-fPIC', '-O2', '-isystem', '/home/phavaej/anaconda3/envs/fenicsx-env/include', '-fPIC', '-O2', '-isystem', '/home/phavaej/anaconda3/envs/fenicsx-env/include', '-fPIC', '-I/home/phavaej/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ffcx/codegeneration', '-I/home/phavaej/anaconda3/envs/fenicsx-env/include/python3.12', '-c', 'libffcx_forms_bd87454e19278e453f4db79408a381bd7f579eaf.c', '-o', './libffcx_forms_bd87454e19278e453f4db79408a381bd7f579eaf.o', '-std=c17', '-O2', '-g0']' returned non-zero exit status 1.

The above exception was the direct cause of the following exception:

distutils.errors.DistutilsExecError: command '/usr/bin/gcc' failed with exit code 1

During handling of the above exception, another exception occurred:

distutils.errors.CompileError: command '/usr/bin/gcc' failed with exit code 1

During handling of the above exception, another exception occurred:

  File "/home/phavaej/Fenics/test_interpolation_block_coupling/coupling_example_org.py", line 132, in <module>
    residual = dolfinx.fem.form(ufl.extract_blocks(F), entity_maps=entity_maps)
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
cffi.VerificationError: CompileError: command '/usr/bin/gcc' failed with exit code 1

I think this issue may be related to the dimensions or accessing the correct subspace of the mixed domain, but I struggled to resolve it. I am using dolphinx v 0.9.

I would greatly appreciate any help and guidance on this.

Kind regards,
Peyman

2

Please provide a minimal reproducible example, as it is not possible to reproduce your error with only this snippet.

3

Here is the MWE:

from mpi4py import MPI
import dolfinx
import gmsh
import numpy as np
import ufl

gmsh_model_rank = 0
mesh_comm = MPI.COMM_WORLD
c_y = 1
R = 0.5
potential_contact_marker = 2
displacement_marker = 1
res = 0.2
order = 2
refinement_level = 2

# Initialize gmsh
gmsh.initialize()
if MPI.COMM_WORLD.rank == 0:
    # Create disk and subtract top part
    membrane = gmsh.model.occ.addDisk(0, c_y, 0, R, R)
    square = gmsh.model.occ.addRectangle(-R, c_y, 0, 2 * R, 1.1 * R)
    gmsh.model.occ.synchronize()
    new_tags, _ = gmsh.model.occ.cut([(2, membrane)], [(2, square)])
    gmsh.model.occ.synchronize()

    # Split boundary into two components
    boundary = gmsh.model.getBoundary(new_tags, oriented=False)
    contact_boundary = []
    dirichlet_boundary = []
    for bnd in boundary:
        mass = gmsh.model.occ.getMass(bnd[0], bnd[1])
        if np.isclose(mass, np.pi * R):
            contact_boundary.append(bnd[1])
        elif np.isclose(mass, 2 * R):
            dirichlet_boundary.append(bnd[1])
        else:
            raise RuntimeError("Unknown boundary")

    # Tag physical groups for the surface
    for i, tags in enumerate(new_tags):
        gmsh.model.addPhysicalGroup(tags[0], [tags[1]], i+1)

    # Tag physical groups for the boundary
    gmsh.model.add_physical_group(1, contact_boundary, potential_contact_marker)
    gmsh.model.add_physical_group(1, dirichlet_boundary, displacement_marker)

    # Create higher resolution mesh at the contact boundary
    distance_field = gmsh.model.mesh.field.add("Distance")
    gmsh.model.mesh.field.setNumbers(distance_field, "EdgesList", contact_boundary)
    threshold = gmsh.model.mesh.field.add("Threshold")
    gmsh.model.mesh.field.setNumber(threshold, "IField", distance_field)
    gmsh.model.mesh.field.setNumber(threshold, "LcMin", res)
    gmsh.model.mesh.field.setNumber(threshold, "LcMax", 20 * res)
    gmsh.model.mesh.field.setNumber(threshold, "DistMin", 0.075 * R)
    gmsh.model.mesh.field.setNumber(threshold, "DistMax", 0.5 * R)
    gmsh.model.mesh.field.setAsBackgroundMesh(threshold)

    # Generate mesh, make second order and refine
    gmsh.model.mesh.generate(2)
    gmsh.model.mesh.setOrder(order)
    for _ in range(refinement_level):
        gmsh.model.mesh.refine()
        gmsh.model.mesh.setOrder(order)
        
    # We inspect the generated mesh and markers
omega, ct, ft = dolfinx.io.gmshio.model_to_mesh(gmsh.model, mesh_comm, gmsh_model_rank, gdim=2)
gmsh.finalize()

tdim = omega.topology.dim
fdim = tdim - 1
gdim = omega.geometry.dim
gamma, gamma_to_omega = dolfinx.mesh.create_submesh(omega, fdim, ft.find(potential_contact_marker))[
        0:2]


V = dolfinx.fem.functionspace(omega, ("Lagrange", 1, (omega.geometry.dim, )))
Q = dolfinx.fem.functionspace(gamma, ("Lagrange", 1))
W = ufl.MixedFunctionSpace(V, Q)

facet_imap = omega.topology.index_map(ft.dim)
num_facets = facet_imap.size_local + facet_imap.num_ghosts
omega_to_gamma = np.full(num_facets, -1)
omega_to_gamma[gamma_to_omega] = np.arange(len(gamma_to_omega))
entity_maps = {gamma: omega_to_gamma}
# -

dx = ufl.Measure("dx", domain=omega)
ds = ufl.Measure("ds", domain=omega, subdomain_data=ft,
                 subdomain_id=potential_contact_marker)

E = dolfinx.fem.Constant(omega, dolfinx.default_scalar_type(2e4))
nu = dolfinx.fem.Constant(omega, 0.4)
alpha = dolfinx.fem.Constant(omega, dolfinx.default_scalar_type(0.1))


h = 0.13
x, y = ufl.SpatialCoordinate(omega)
g = y + dolfinx.fem.Constant(omega, dolfinx.default_scalar_type(h))
uD = 0.72

n = ufl.FacetNormal(omega)
n_g = dolfinx.fem.Constant(omega, np.zeros(gdim, dtype=dolfinx.default_scalar_type))
n_g.value[-1] = -1
f = dolfinx.fem.Constant(omega, np.zeros(gdim, dtype=dolfinx.default_scalar_type))


v, w = ufl.TestFunctions(W)
u = dolfinx.fem.Function(V, name="Displacement")
psi = dolfinx.fem.Function(Q, name="Latent_variable")
psi_k = dolfinx.fem.Function(Q, name="Previous_Latent_variable")

mu = E / (2.0 * (1.0 + nu))
lmbda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))

def epsilon(w):
    return ufl.sym(ufl.grad(w))

def sigma(w, mu, lmbda):
    ew = epsilon(w)
    gdim = ew.ufl_shape[0]
    return 2.0 * mu * epsilon(w) + lmbda * ufl.tr(ufl.grad(w)) * ufl.Identity(gdim)


F = alpha * ufl.inner(sigma(u, mu, lmbda), epsilon(v)) * dx
F -= alpha * ufl.inner(f, v) * dx
F += -ufl.inner(psi - psi_k, ufl.dot(v, n)) * ds
F += ufl.inner(ufl.dot(u, n_g), w) * ds
F += ufl.inner(ufl.exp(psi), w) * ds - ufl.inner(g, w) * ds
F += ufl.inner(ufl.grad(psi), ufl.grad(w)) * ds # Addiotinal term for Poisson equation
residual = dolfinx.fem.form(ufl.extract_blocks(F), entity_maps=entity_maps)