Solve equation in only one phase of a two-phase problem

General
1

I have a two-phase system for FSI simulation. With everything built on a single big mesh, I need to compute variables on the solid phase only as shown in the attached figure.

two_phase_FSI
two_phase_FSI1981×511 433 KB

(a) mask and velocity vector arrows, (b) velocity magnitude. (c) result of \omega=0.5\nabla \times u.

The test in the figure shows that when the velocity field u is given, the rotation velocity curlU of the solid is computed via projection of \omega=0.5\nabla \times u. Interpolation cannot do this as the velocity is continuous but not differentiable on the interface.

To do so, I have 2 scenarios:

  1. Do projection on the big mesh with the help of a function mask for solid and apply Dirichlet BC for fluid-only DOFs as
a  = fem.form(ufl.inner(Pv, w) * mask * ufl.dx)
L  = fem.form(ufl.inner( v, w) * mask * ufl.dx)

problem = LinearProblem(a, L, bcs, u=target_func,
        petsc_options={"ksp_type": "preonly", 
                       "pc_type": "lu", 
                       "pc_factor_mat_solver_type": "mumps"})
problem.solve()

The full MWE is attached in the end.

  1. Create a submesh for solid with submesh = mesh.create_submesh(big_mesh, tdim, solid_cells)[0]. First, I interpolate u from the big mesh to the submesh. Then, solve the projection on the submesh to get curlU. Finally, interpolate curlU to the function on the big mesh. The whole process involves 2 non-matching meshes interpolation.

In my application, I do two-phase FSI simulations in 3D with 2nd order iso-parametric meshes. Scenario 1 gets really slow as it involves too many unnecessary DOFs. Scenario 2 is fast but it often crashes with interpolation related GJK or Newton method failed to converge runtime errors. I sense the scenario 1 could be tuned for this problem, but still I am looking for a hand on this. Thanks.

MWE with DOLFINx-v0.8.0:

from dolfinx import mesh, io, fem
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import numpy as np
import ufl

comm = MPI.COMM_WORLD
eps = 1e-8

def project(v, target_func, mask, bcs=[]) -> None:

    V = target_func.function_space
    w  = ufl.TestFunction(V)
    Pv = ufl.TrialFunction(V)
    a  = fem.form(ufl.inner(Pv, w) * mask * ufl.dx)
    L  = fem.form(ufl.inner( v, w) * mask * ufl.dx)

    problem = LinearProblem(a, L, bcs, u=target_func,
            petsc_options={"ksp_type": "preonly", 
                           "pc_type": "lu", 
                           "pc_factor_mat_solver_type": "mumps"})
    problem.solve()
    target_func.x.scatter_forward()
    
def solid_domain(x):
    return np.logical_and(
        np.logical_and(x[0] > 0.25 - eps, x[0] < 0.75 + eps),
        np.logical_and(x[1] > 0.25 - eps, x[1] < 0.75 + eps))
    
_mesh = mesh.create_unit_square(comm, 32, 32)
tdim = _mesh.topology.dim

# Mask on a DG-0 function space
DG0 = fem.functionspace(_mesh, ("DG", 0))
mask = fem.Function(DG0)
_mesh.topology.create_connectivity(tdim, tdim)

solid_cells = mesh.locate_entities(_mesh, tdim, solid_domain)
solid_dofs = fem.locate_dofs_topological(DG0, tdim, solid_cells)
mask.x.array[solid_dofs] = 1.0

with io.XDMFFile(comm, "output/mask.xdmf", "w") as file:
    file.write_mesh(_mesh)
    file.write_function(mask)

# Build velocity field for 2 phases
omega = 1.0      # rotation velocity
center = ufl.as_vector([0.5, 0.5])
x   = ufl.SpatialCoordinate(_mesh) - center
r   = ufl.sqrt((x[0] - 0.5)**2 + (x[1] - 0.5)**2)
tri = ufl.as_vector([- x[1] / (r + eps), x[0] / (r + eps)])
V = fem.functionspace(_mesh, ("Lagrange", 1, (tdim,)))
Q = fem.functionspace(_mesh, ("Lagrange", 1))
u_fluid = fem.Function(V)
u_solid = fem.Function(V)
u       = fem.Function(V)
mask_V  = fem.Function(V)

solid_dofs_V = fem.locate_dofs_topological(V, tdim, solid_cells)
for d in range(tdim):
    mask_V.x.array[solid_dofs_V * tdim + d] = 1.0
mask_V.x.scatter_forward()

expr = fem.Expression(omega * r * tri, V.element.interpolation_points())
u_solid.interpolate(expr)
u_solid.x.scatter_forward()

u_fluid.interpolate(lambda x: ((x[0] - 0.5), (x[1] - 0.5)))
u_fluid.x.scatter_forward()

u.x.array[:] = (u_solid.x.array[:] * mask_V.x.array[:] +
                u_fluid.x.array[:] * (1 - mask_V.x.array[:]))

with io.XDMFFile(comm, "output/u.xdmf", "w") as file:
    file.write_mesh(_mesh)
    file.write_function(u)
    
curlU = fem.Function(Q)
total_dofs_Q = (Q.dofmap.index_map.size_local + 
                Q.dofmap.index_map.num_ghosts)
solid_dofs_Q = fem.locate_dofs_topological(Q, tdim, solid_cells)
fluid_only_dofs_Q = np.setdiff1d(np.arange(total_dofs_Q), solid_dofs_Q)

expr = 1. / 2. * ufl.curl(u)
f = fem.Function(Q)
bc_f = fem.dirichletbc(f, fluid_only_dofs_Q)
project(expr, curlU, mask, bcs=[bc_f])

with io.XDMFFile(comm, "output/omega.xdmf", "w") as file:
    file.write_mesh(_mesh)
    file.write_function(curlU)
2

You don’t need to interpolate u from the parent to the sub mesh. At least if you upgrade to 0.9, you can have u in a variational form on the submesh.

See for instance: Multiphysics: Solving PDEs on subdomains — FEniCS Workshop

3

I had a go with the matched section Integration with function from parent and sub-mesh in the tutorial and the result looks correct with my case. Please correct me if I am wrong with the following idea:

u is defined on the parent mesh, while omega function, trail function and test function are defined on the submesh. It needs extra entity_maps info to explicitly formulate and assemble to A and b. Then solve projection on the submesh to get omega.

Several follow-up questions:
(1) I noticed entity_map argument of fem.form also exists in v0.8.0 with the same description, but I failed to call a_solid, L_solid = dolfinx.fem.form(ufl.system(F), entity_maps=entity_maps) as that in v0.9.0. Is it a v0.9.0-only related feature?
(2) With the use of entity_maps of fem.form, the Ax=b problem can only be solved explicitly rather than to use LinearProblem wrapper?
(3) My practical application is FSI + ALE, so the mesh moves in every step. If with a submesh, any idea to associate the movement of a submesh to the parent mesh?

4

You don’t need to project anything. You can leave u in the variational form, and it will behave as any function defined in a standard variational form.

Various fixes related to submeshes, as well as a slight api change happened from 0.8 to 0.9

You can pass in compiled forms to LinearProblem.

Then you should interpolate any displacement in the parent mesh onto the submesh and move the mesh afterwards

5

Assume I have two functions f and fs defined on parent mesh and submesh respectively. Is it able to find a DOF map from submesh to parent mesh? Say, the dof_sub_to_parent_map array in the following code:

tdim = cell_tags.dim
submesh, cell_map, vertex_map, _ = dolfinx.mesh.create_submesh(mesh, tdim, cell_tags.find(marker))

# Function space and functions
V  = dolfinx.fem.functionspace(mesh, ("Lagrange", 1, (tdim,)))
Vs = dolfinx.fem.functionspace(submesh, ("Lagrange", 1, (tdim,)))
f  = dolfinx.fem.Function(V)
fs = dolfinx.fem.Function(Vs)

# Wanted: dof_sub_to_parent_map
dof_sub_to_parent_map = ...
f.x.array[dof_sub_to_parent_map] = fs.x.array[:]
6

You can make this function if you want to.
However, it is much quicker and easier call

submesh, parent_cells, _, _ = create_submesh(mesh, tdim, cells)
f.interpolate(fs, cells0=np.arange(len(parent_cells)), cells1=parent_cells)

as done in:

7

It works with interpolation. My understanding is the interpolation between submesh and parent mesh is matching meshes interpolation. An unmatured question: for the quoted interpolation, does DOLFINx use 1-to-1 assigning on specified cells or interpolate through compute?

8

Yes it is matching interpolation.

As the spaces for both parent and submesh are Lagrange 1, it becomes a task of assigning degrees of freedom (not tabulate/compute basis functions on the input space)