Domain Interface for Stokes

General
,
1

Hi, I am new in Fenics. I want to implement the Stokes problem in stenosis, which has two domains. I want to compute the solution in whole domain then along with this not able to find the solution in the separate domain. Moreover, to compute the error with the whole domain solution for u in domain 1 and domain 2.
Code is here, which works properly but does not satisfy the solution because the parabolic inlet but not sure about the solution. I was also not able to find the solution for the separate domains.

FB removed code

2

What do you mean by this?

Secondly please attach plots of current flow field to highlight what is wrong.

Finally, please format your code with 3x` encapsulation, ie

```python
from mpi4py import MPI
#add code here
```

and please list what version of FEniCSx and multiphenicsx you are using, along with how you have installed the packages.

3

Thank you for your response. I am using google colab for the simulation. Multiphenicsx library direct linked to Colab.

steno
steno1280×568 108 KB

This is the geometry with two domains.
Firstly interested to find the solution in the whole domain.
At the next step I already define the mesh for the both domains, find the solution in the domain 1 separately and domain 2 as well.
pressure
pressure1599×945 230 KB

vel1
vel11600×801 172 KB

4 
from mpi4py import MPI
import gmsh
import numpy as np
from dolfinx import fem, plot
import ufl
import dolfinx.io
from dolfinx.io import gmshio
from dolfinx.io import XDMFFile
from mpi4py import MPI
from ufl import div, dx, grad, inner, SpatialCoordinate, FacetNormal
from dolfinx.fem import (Constant, Function, FunctionSpace, locate_dofs_geometrical)
import mpi4py.MPI
import petsc4py.PETSc
import multiphenicsx.fem
from petsc4py import PETSc
import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, fem, la
from dolfinx.fem import (
Constant,
Function,
dirichletbc,
extract_function_spaces,
form,
functionspace,
locate_dofs_topological,
)
from dolfinx.fem.petsc import assemble_matrix_block, assemble_vector_block
from dolfinx.io import XDMFFile
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from ufl import div, dx, grad, inner
import viskex

Initialize Gmsh
gmsh.initialize()
gmsh.model.add(“mortar”)

Define parameters
ls = 5
Xi = 100 # um
Xo = 100 # um
L = 100.0 # um
x0 = Xi + L / 2.0
R = 50.0 # um
f0 = 0.5 # 0–1

Add points
p1 = gmsh.model.geo.addPoint(0, 0, 0, ls)
p2 = gmsh.model.geo.addPoint(Xi, 0, 0, ls)
p3 = gmsh.model.geo.addPoint(Xi, R, 0, ls)
p4 = gmsh.model.geo.addPoint(0, R, 0, ls)
p5 = gmsh.model.geo.addPoint(Xi + L, 0, 0, ls)
p6 = gmsh.model.geo.addPoint(Xi + L + Xo, 0, 0, ls)
p7 = gmsh.model.geo.addPoint(Xi + L + Xo, R, 0, ls)
p8 = gmsh.model.geo.addPoint(Xi + L, R, 0, ls)

Define lines
l1 = gmsh.model.geo.addLine(p1, p2)
l2 = gmsh.model.geo.addLine(p2, p3)
l3 = gmsh.model.geo.addLine(p3, p4)
l4 = gmsh.model.geo.addLine(p4, p1)
l5 = gmsh.model.geo.addLine(p6, p5)
l6 = gmsh.model.geo.addLine(p7, p6)
l7 = gmsh.model.geo.addLine(p8, p7)
#l8 = gmsh.model.geo.addLine(p8, p5)

Define the lower stenosis profile as a B-spline
nPoints = 21
pListLower = [p3]
for i in range(1, nPoints + 1):
x = Xi + L * i / (nPoints + 1)
y = R * (1 - f0 / 2.5 * (1 + np.cos(2.0 * np.pi * (x - x0) / L)))
pListLower.append(gmsh.model.geo.addPoint(x, y, 0, ls))
pListLower.append(p8)
l8 = gmsh.model.geo.addBSpline(pListLower)

Define the upper stenosis profile as a B-spline
pListUpper = [p2]
for i in range(1, nPoints + 1):
x = Xi + L * i / (nPoints + 1)
y = R - (R * (1 - f0 / 2.5 * (1 + np.cos(2.0 * np.pi * (x - x0) / L))))
pListUpper.append(gmsh.model.geo.addPoint(x, y, 0, ls))
pListUpper.append(p5)
l9 = gmsh.model.geo.addBSpline(pListUpper)

loop_rect = gmsh.model.geo.addCurveLoop([l4, l1, l2, l3])
rectangle = gmsh.model.geo.addPlaneSurface([loop_rect])

loop_stenosis = gmsh.model.geo.addCurveLoop([l2, l8, l7, l6, l5, -l9])
stenosis = gmsh.model.geo.addPlaneSurface([loop_stenosis])

gmsh.model.geo.synchronize()

gmsh.model.addPhysicalGroup(1, [l4], 1, “Inlet”)
gmsh.model.addPhysicalGroup(1, [l2], 2, “Interface”)
gmsh.model.addPhysicalGroup(1, [l6], 3, “Outlet”)
gmsh.model.addPhysicalGroup(1, [l3, l8, l7, l5, l9, l1], 4, “Wall”)
gmsh.model.addPhysicalGroup(2, [rectangle], name = “left_rectangle”)
gmsh.model.addPhysicalGroup(2, [stenosis], name = “right_stenosis”)

gmsh.model.mesh.generate(2)

Save mesh to file
gmsh.write(“new_geometry.msh”)

partitioner = dolfinx.mesh.create_cell_partitioner(dolfinx.mesh.GhostMode.shared_facet)
mesh, subdomains, boundaries = dolfinx.io.gmshio.model_to_mesh(
gmsh.model, comm=mpi4py.MPI.COMM_WORLD, rank=0, gdim=2, partitioner=partitioner)
gmsh.finalize()

mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim)

cells_Omega1 = subdomains.indices[subdomains.values == 1]
cells_Omega2 = subdomains.indices[subdomains.values == 2]
facets_partial_Omega = boundaries.indices[boundaries.values == 1]
facets_Gamma = boundaries.indices[boundaries.values == 2]

Define associated measures
dx = ufl.Measure(“dx”)(subdomain_data=subdomains)
dS = ufl.Measure(“dS”)(subdomain_data=boundaries)
dS = dS(2) # restrict to the interface, which has facet ID equal to 2

viskex.dolfinx.plot_mesh(mesh)
viskex.dolfinx.plot_mesh_tags(mesh, subdomains, “subdomains”)
viskex.dolfinx.plot_mesh_tags(mesh, boundaries, “boundaries and interfaces”)

Define associated measures
#dx = ufl.Measure(“dx”, subdomain_data=subdomains)
#dS = ufl.Measure(“dS”, subdomain_data=boundaries)
#ds = ufl.Measure(“ds”, subdomain_data=boundaries)
#n = FacetNormal(mesh)
#locate cells in the two subdomains

cells_Omega1 = subdomains.indices[subdomains.values == 1]

cells_Omega2 = subdomains.indices[subdomains.values == 2]

#locate boundary facets

facets_inlet_boundary = boundaries.indices[boundaries.values == 1]

facets_interface = boundaries.indices[boundaries.values == 2]

facets_outlet_boundary = boundaries.indices[boundaries.values == 3]

facets_wall_boundary = boundaries.indices[boundaries.values == 4]

multiphenicsx plot
#DEFINITION OF THE MEASURES

Define associated measures
dx = ufl.Measure(“dx”)(subdomain_data=subdomains)

dS = ufl.Measure(“dS”)(subdomain_data=boundaries)

dS = dS(2) # restrict to the interface, which has facet ID equal to 2

n= FacetNormal(mesh)

#DEFINITION OF THE SPACES

Define function spaces
#We firstly define vector and finite element for upper and lower half of the domain

V1_element = basix.ufl.element(“Lagrange”, mesh.basix_cell(), 2, shape=(mesh.geometry.dim, ))
Q1_element = basix.ufl.element(“Lagrange”, mesh.basix_cell(), 1)
W1_element = basix.ufl.mixed_element([V1_element, Q1_element])
W1 = dolfinx.fem.functionspace(mesh, W1_element)
V1, _ = W1.sub(0).collapse()
Q1, _ = W1.sub(1).collapse()

V2_element = basix.ufl.element(“Lagrange”, mesh.basix_cell(), 2, shape=(mesh.geometry.dim, ))
Q2_element = basix.ufl.element(“Lagrange”, mesh.basix_cell(), 1)
W2_element = basix.ufl.mixed_element([V2_element, Q2_element])
W2 = dolfinx.fem.functionspace(mesh, W1_element)
V2, _ = W2.sub(0).collapse()
Q2, _ = W2.sub(1).collapse()

#space of the lagrange multiplie

M_element = basix.ufl.element(“Lagrange”, mesh.basix_cell(), 2, shape=(mesh.geometry.dim, ))

M = dolfinx.fem.functionspace(mesh, M_element)

mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)

mesh.topology.create_connectivity(subdomains.dim, subdomains.dim)

#define the restrictions

dofs_V1_Omega1 = dolfinx.fem.locate_dofs_topological(V1, subdomains.dim, cells_Omega1) #locate dofs of the V1space

dofs_V2_Omega2 = dolfinx.fem.locate_dofs_topological(V2, subdomains.dim, cells_Omega2) #locate dofs of the V2space

dofs_M_Gamma = dolfinx.fem.locate_dofs_topological(M, boundaries.dim, facets_interface) #locate dofs of the LM space

dofs_Q1_Omega1 = dolfinx.fem.locate_dofs_topological(Q1, subdomains.dim, cells_Omega1) #locate dofs of the Q1space

dofs_Q2_Omega2 = dolfinx.fem.locate_dofs_topological(Q2, subdomains.dim, cells_Omega2) #locate dofs of the Q2space

restriction_V1_Omega1 = multiphenicsx.fem.DofMapRestriction(V1.dofmap, dofs_V1_Omega1)

restriction_V2_Omega2 = multiphenicsx.fem.DofMapRestriction(V2.dofmap, dofs_V2_Omega2)

restriction_Q1_Omega1 = multiphenicsx.fem.DofMapRestriction(Q1.dofmap, dofs_Q1_Omega1)

restriction_Q2_Omega2 = multiphenicsx.fem.DofMapRestriction(Q2.dofmap, dofs_Q2_Omega2)

restriction_M_Gamma = multiphenicsx.fem.DofMapRestriction(M.dofmap, dofs_M_Gamma)

restriction = [restriction_V1_Omega1, restriction_V2_Omega2, restriction_Q1_Omega1, restriction_Q2_Omega2, restriction_M_Gamma]

Define function spaces
V_el = basix.ufl.element(“Lagrange”, mesh.basix_cell(), 2, shape=(mesh.geometry.dim,))

Q_el = basix.ufl.element(“Lagrange”, mesh.basix_cell(), 1)

W_el = basix.ufl.mixed_element([V_el, Q_el])

W = dolfinx.fem.functionspace(mesh, W_el)

V, _ = W.sub(0).collapse()

Q, _ = W.sub(1).collapse()

Define boundary conditions
Define parabolic velocity profile for the inlet
def parabolic_inlet(x):

Assume the inlet height (vertical distance) is 50 units
R = 50.0

y = x[1] # y-coordinate (vertical direction)

The parabolic profile is zero at y = 0 and y = R, and maximum at y = R/2
return np.stack((0.5 * y * (R - y) / R**2, np.zeros_like(y)))

u_inlet = Function(V)

u_inlet.interpolate(parabolic_inlet)

Apply the inlet boundary condition
inlet_dofs = locate_dofs_topological((W.sub(0), V), mesh.topology.dim - 1, facets_inlet_boundary)

bc_inlet = dirichletbc(u_inlet, inlet_dofs, W.sub(0))

Pressure boundary condition at the outlet (zero pressure)
p_outlet = dolfinx.fem.Function(Q) # Create a function in the pressure space Q

p_outlet.interpolate(lambda x: np.zeros_like(x[0])) # Set pressure to zero

dofs_outlet = dolfinx.fem.locate_dofs_topological(Q, mesh.topology.dim - 1, facets_outlet_boundary)

bc_outlet = dolfinx.fem.dirichletbc(p_outlet, dofs_outlet)

Zero normal velocity for the outlet
#u_outlet = Function(V)

#u_outlet.interpolate(lambda x: np.zeros((mesh.geometry.dim, x.shape[1])))

#outlet_dofs = locate_dofs_topological((W.sub(0), V), mesh.topology.dim - 1, facets_outlet_boundary)

#bc_outlet = dirichletbc(u_outlet, outlet_dofs, W.sub(0))

No-slip condition on the walls
u_walls = Function(V)

u_walls.interpolate(lambda x: np.zeros((mesh.geometry.dim, x.shape[1])))

wall_dofs = locate_dofs_topological((W.sub(0), V), mesh.topology.dim - 1, facets_wall_boundary)

bc_walls = dirichletbc(u_walls, wall_dofs, W.sub(0))

If there’s an interface boundary that requires specific handling, define it here:
Example:
interface_dofs = locate_dofs_topological((W.sub(0), V), mesh.topology.dim - 1, facets_interface)
bc_interface = dirichletbc(, interface_dofs, W.sub(0))
Combine boundary conditions
bcs = [bc_inlet, bc_walls, bc_outlet]

Define variational problem
u, p = ufl.TrialFunctions(W)

v, q = ufl.TestFunctions(W)

mu = fem.Constant(mesh, PETSc.ScalarType(1.0))

f = fem.Constant(mesh, PETSc.ScalarType((0.0, 0.0)))

a = mu * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx - p * ufl.div(v) * ufl.dx + q * ufl.div(u) * ufl.dx

L = ufl.inner(v, f) * ufl.dx

Solve the linear problem
uh = fem.Function(W)

problem = fem.petsc.LinearProblem(a, L, bcs=bcs, u=uh, petsc_options={

“ksp_type”: “preonly”,

“pc_type”: “lu”,

“pc_factor_mat_solver_type”: “mumps”,

})

problem.solve()

uh = problem.solve()

us = uh.split()[0].collapse()

ps = uh.split()[1].collapse()

us.name = “vel”

ps.name = “pres”

Plotting using viskex.dolfinx
viskex.dolfinx.plot_vector_field(us)

viskex.dolfinx.plot_scalar_field(ps)

u_ex1_norm = np.sqrt(mesh.comm.allreduce(

dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(us, us) * dx(1))), op=mpi4py.MPI.SUM))

u_ex2_norm = np.sqrt(mesh.comm.allreduce(

dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(us, us) * dx(2))), op=mpi4py.MPI.SUM))

p_ex1_norm = np.sqrt(mesh.comm.allreduce(

dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(ps, ps) * dx(1))), op=mpi4py.MPI.SUM))

p_ex2_norm = np.sqrt(mesh.comm.allreduce(

dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(ps, ps) * dx(2))), op=mpi4py.MPI.SUM))
5

This is probably a question for @francesco-ballarin as I am not very familiar with the multiphenics API.

Secondly, I would advice you to try to minimize the code.
Currently you use a mesh built by GMSH, but you could illustrate the problem on a unit square (which you want to divide in two pieces after solving for the Stokes flow on the whole domain). This would massively shorten your code.

It seems like you have two different issues:

  1. The standard solution of Stokes problem on the whole domain doesn’t work.

Here I would as I already mention start with a simpler domain to ensure that you set boundary conditions correctly. I also guess that you do not need multiphenicsx for this first step.

  1. How to use multiphenicsx to compute the error in two subdomains separately.

For this I do not understand why you need to use multiphenicsx. You could simply restrict your volume integrals during error computations by creating a restricted integration measure.

dx = ufl.Measure("dx", domain=mesh)
error_1 = dolfinx.fem.form(ufl.inner(u-u_ex,u-u_ex)*dx(1))
error_2 = dolfinx.fem.form(ufl.inner(u-u_ex,u-u_ex)*dx(2))
print(mesh.comm.allreduce(dolfinx.fem.assemble_scalar(error_1), op=MPI.SUM))
print(mesh.comm.allreduce(dolfinx.fem.assemble_scalar(error_2), op=MPI.SUM))
6

I can confirm that you are never using the defined restriction, so it seems to me that you are never really using multiphenicsx. Apart from that, I couldn’t really understand the description of the problem.

7

without multiphenicsx how we can do the domain decompositon and then find the solution in the seperate domain. thank you

8

Again, you seem to be switching back and forth between notions.

In your previous post, you stated that you wanted to compute the error in a subdomain, you did not specify that you want to solve an equation in each of the subdomains.

However, you can do that with either multiphenicsx, or the dolfinx.mesh.create_submesh function.

However, you would have to specify how your domain decomposition method is supposed to work. What are the boundary conditions on the interface between the two sub-domains that you want to solve on individually.

And as has already been stated, you haven’t even tried to use the multiphenics restriction to solve a problem on a restricted domain in the code you have presented to us.