Dear all,
I am trying to solve the Stokes equations with periodic boundary conditions similar to this post, but with a nested, iterative solver. Therefore, I tried to adapt this tutorial example.
However, the results are qualitatively not as expected. The field distribution is strange at some DOFs, and the velocity components are somehow not exactly periodic. Is it not possible to give two multipoint constraints simultaneously?
In the following, you can find the implemented code. I am sorry that the code is not minimal, but I am not sure how to present the problem differently. I am confident that everything is implemented correctly except for the nested solver, as it converges to the correct solution with other solvers.
What am I missing? Has anyone an idea?
import gmsh
import dolfinx
import dolfinx_mpc
from dolfinx.fem import functionspace, Function, Constant, form
from dolfinx.io.gmshio import model_to_mesh
from dolfinx.mesh import locate_entities_boundary, meshtags
import ufl
import numpy as np
from mpi4py import MPI
from ufl import grad, div, inner, Measure
import basix.ufl
from petsc4py import PETSc
comm = MPI.COMM_WORLD
rank = comm.rank
#############################################
# Mesh generation #
#############################################
gmsh.initialize()
L = 1 # length of rectangle
f = 0.3 # volume fraction of inclusion
r = np.sqrt(f*L**2/np.pi) # radius of inclusion
gdim = 2 # Geometric dimension of the mesh
fdim = gdim - 1 # facet dimension
if rank == 0:
inclusion = gmsh.model.occ.addDisk(L / 2, L / 2, 0, r, r)
matrix = gmsh.model.occ.addRectangle(0,0,0, L, L)
whole_domain = gmsh.model.occ.fragment([(2, matrix)], [(2, inclusion)])
gmsh.model.occ.synchronize()
matrix_physical = gmsh.model.addPhysicalGroup(gdim, [matrix], tag = 1)
inclusion_physical = gmsh.model.addPhysicalGroup(gdim, [inclusion], tag = 2)
# Get the interface elements between rectangle and circle and tag
interface_elements = gmsh.model.getBoundary([(gdim, inclusion_physical)])
interface_physical = gmsh.model.addPhysicalGroup(fdim, (5, 5), tag = 9)
background_surfaces = []
inclusion_surfaces = []
for domain in whole_domain[0]:
com = gmsh.model.occ.getCenterOfMass(domain[0], domain[1])
mass = gmsh.model.occ.getMass(domain[0], domain[1])
if np.isclose(mass, np.pi * (r ** 2)): # identify inclusion by its mass
inclusion_surfaces.append(domain)
else:
background_surfaces.append(domain)
gmsh.model.mesh.field.add("Distance", 1)
edges = gmsh.model.getBoundary(inclusion_surfaces, oriented=False)
gmsh.model.mesh.field.setNumbers(1, "EdgesList", [e[1] for e in edges])
gmsh.model.mesh.field.add("Threshold", 2)
gmsh.model.mesh.field.setNumber(2, "IField", 1)
gmsh.model.mesh.field.setNumber(2, "LcMin", 0.01)
gmsh.model.mesh.field.setNumber(2, "LcMax", 0.03)
gmsh.model.mesh.field.setNumber(2, "DistMin", 0.01)
gmsh.model.mesh.field.setNumber(2, "DistMax", 0.1)
gmsh.model.mesh.field.setAsBackgroundMesh(2)
gmsh.option.setNumber("Mesh.Algorithm", 6)
gmsh.model.mesh.generate(gdim)
domain, ct, ft = model_to_mesh(gmsh.model, comm, rank=0, gdim=gdim)
gmsh.finalize()
dx = Measure('dx', domain=domain, subdomain_data=ct)
domain.topology.create_connectivity(gdim, gdim)
#############################################
# Function spaces #
#############################################
# Element formulation
P1 = basix.ufl.element("Lagrange", domain.basix_cell(), 1) # pressure
P2 = basix.ufl.element("Lagrange", domain.basix_cell(), 2, shape=(domain.geometry.dim,)) # velocity
V, Q = functionspace(domain, P2), functionspace(domain, P1)
# define function, trial function and test function
(ve, pr) = ufl.TrialFunction(V), ufl.TrialFunction(Q)
(dve, dpr) = ufl.TestFunction(V), ufl.TestFunction(Q)
#############################################
# Boundary conditions #
#############################################
p_gradient = Constant(domain, PETSc.ScalarType((1,0)))
# set all fluid velocities within the inclusion and at the interface to zero
dofs_v_ct2 = dolfinx.fem.locate_dofs_topological(V , ct.dim, ct.find(2))
bcs_v_ct2 = dolfinx.fem.dirichletbc(np.zeros(domain.geometry.dim, dtype=PETSc.ScalarType), dofs_v_ct2, V)
# set all pressures within inclusion to zero, but not at the interface
dofs_p_ct2 = dolfinx.fem.locate_dofs_topological(Q , ct.dim, ct.find(2)) # find pressure dofs in cell tag 2
dofs_p_ft9 = dolfinx.fem.locate_dofs_topological(Q , ft.dim, ft.find(9)) # find pressure dofs for facet tag 9
dofs_p_fin = np.array(list(set(dofs_p_ct2).difference(dofs_p_ft9))) # delete pressure dofs of ft9 from ct2 for the first array # final dofs of the pressure in cell tag2 (without facet tag)
bcs_p_ct2 = dolfinx.fem.dirichletbc(PETSc.ScalarType(0.0), dofs_p_fin, Q) # set dirichlet bc
bcs = [bcs_v_ct2, bcs_p_ct2]
## periodic boundary conditions
###############################
pbc_directions = []
pbc_slave_tags = []
pbc_is_slave = []
pbc_is_master = []
pbc_meshtags = []
pbc_slave_to_master_maps = []
mpc_p = dolfinx_mpc.MultiPointConstraint(Q)
mpc_u = dolfinx_mpc.MultiPointConstraint(V)
def generate_pbc_slave_to_master_map(i):
def pbc_slave_to_master_map(x):
out_x = x.copy()
out_x[i] = x[i] - L
return out_x
return pbc_slave_to_master_map
def generate_pbc_is_slave(i):
return lambda x: np.isclose(x[i], L)
def generate_pbc_is_master(i):
return lambda x: np.isclose(x[i], 0.0)
for i in range(gdim):
pbc_directions.append(i)
pbc_slave_tags.append(i + 2)
pbc_is_slave.append(generate_pbc_is_slave(i))
pbc_is_master.append(generate_pbc_is_master(i))
pbc_slave_to_master_maps.append(generate_pbc_slave_to_master_map(i))
facets = locate_entities_boundary(domain, fdim, pbc_is_slave[-1])
arg_sort = np.argsort(facets)
pbc_meshtags.append(meshtags(domain,
fdim,
facets[arg_sort],
np.full(len(facets), pbc_slave_tags[-1], dtype=np.int32)))
for i in range(len(pbc_directions)):
def pbc_slave_to_master_map(x):
out_x = pbc_slave_to_master_maps[i](x)
idx = pbc_is_slave[(i + 1) % len(pbc_directions)](x)
out_x[pbc_directions[i]][idx] = np.nan
return out_x
for ij in range(gdim):
mpc_u.create_periodic_constraint_topological(V.sub(ij), pbc_meshtags[i],
pbc_slave_tags[i],
pbc_slave_to_master_map,
bcs)
mpc_p.create_periodic_constraint_topological(Q, pbc_meshtags[i],
pbc_slave_tags[i],
pbc_slave_to_master_map,
bcs)
# Map intersection(slaves_x, slaves_y) to intersection(masters_x, masters_y),
# i.e. map the slave dof at (1, 1)/(1,1,1) to the master dof at (0, 0)/(0,0,0)
def pbc_slave_to_master_map_corner(x):
out_x = x.copy()
out_x[0] = x[0] - domain.geometry.x[:, 0].max()
out_x[1] = x[1] - domain.geometry.x[:, 1].max()
idx = np.logical_and(pbc_is_slave[0](x), pbc_is_slave[1](x))
out_x[0][~idx] = np.nan
out_x[1][~idx] = np.nan
return out_x
for ij in range(gdim):
mpc_u.create_periodic_constraint_topological(V.sub(ij), pbc_meshtags[1],
pbc_slave_tags[1],
pbc_slave_to_master_map_corner,
bcs)
mpc_p.create_periodic_constraint_topological(Q, pbc_meshtags[1],
pbc_slave_tags[1],
pbc_slave_to_master_map_corner,
bcs)
mpc_p.finalize()
mpc_u.finalize()
#############################################
# Variational formulation #
#############################################
## weak formulation
a00 = inner(grad(ve), grad(dve)) * dx
a01 = -inner(pr, div(dve)) * dx
a10 = -inner(div(ve), dpr) * dx
a11 = None
L0 = -inner(p_gradient, dve) * dx
L1 = inner(Constant(domain, PETSc.ScalarType(0)), dpr) * dx
a = [[form(a00), form(a01)], [form(a10), form(a11)]]
L = [form(L0), form(L1)]
#############################################
# Nested iterative matrix solver #
#############################################
A = dolfinx_mpc.create_matrix_nest(a, [mpc_u, mpc_p])
dolfinx_mpc.assemble_matrix_nest(A, a, [mpc_u, mpc_p], bcs)
A.assemble()
b = dolfinx_mpc.create_vector_nest(L, [mpc_u, mpc_p])
dolfinx_mpc.assemble_vector_nest(b, L, [mpc_u, mpc_p])
# Set Dirichlet boundary condition values in the RHS
dolfinx.fem.petsc.apply_lifting_nest(b, a, bcs)
for b_sub in b.getNestSubVecs():
b_sub.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) # type: ignore
# bcs0 = dolfinx.cpp.fem.bcs_rows(dolfinx.fem.assemble._create_cpp_form(L), bcs)
bcs0 = dolfinx.fem.bcs_by_block(dolfinx.fem.extract_function_spaces(L), bcs)
dolfinx.fem.petsc.set_bc_nest(b, bcs0)
P11 = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(pr * dpr * ufl.dx))
P = PETSc.Mat().createNest([[A.getNestSubMatrix(0, 0), None], [None, P11]]) # type: ignore
P.assemble()
# Create a MINRES Krylov solver and a block-diagonal preconditioner
# using PETSc's additive fieldsplit preconditioner
ksp = PETSc.KSP().create(domain.comm)
ksp.setOperators(A, P)
ksp.setTolerances(rtol=1e-5, atol=1e-7, max_it=500)
ksp.setType("minres")
ksp.getPC().setType("fieldsplit")
ksp.getPC().setFieldSplitType(PETSc.PC.CompositeType.ADDITIVE)
nested_IS = P.getNestISs()
ksp.getPC().setFieldSplitIS(("ve", nested_IS[0][0]),
("pr", nested_IS[0][1]))
# # Configure velocity and pressure sub-solvers
ksp_u, ksp_p = ksp.getPC().getFieldSplitSubKSP()
ksp_u.setType("preonly")
ksp_u.getPC().setType("gamg")
ksp_p.setType("preonly")
ksp_p.getPC().setType("jacobi")
ksp.setFromOptions()
## solve system
Uh = b.copy()
ksp.solve(b, Uh)
print(f"Converged reason: {ksp.getConvergedReason()}")
for Uh_sub in Uh.getNestSubVecs():
Uh_sub.ghostUpdate(
addv=PETSc.InsertMode.INSERT, # type: ignore
mode=PETSc.ScatterMode.FORWARD, # type: ignore
) # type: ignore
# intialize files for the results
vel, p = Function(mpc_u.function_space), Function(mpc_p.function_space)
vtx_pr = dolfinx.io.VTXWriter(domain.comm, "pressure.bp", p, engine="BP4")
vtx_ve = dolfinx.io.VTXWriter(domain.comm, "velocity.bp", vel, engine="BP4")
vel.x.petsc_vec.setArray(Uh.getNestSubVecs()[0].array)
p.x.petsc_vec.setArray(Uh.getNestSubVecs()[1].array)
vel.x.scatter_forward()
p.x.scatter_forward()
# Backsubstitute to update slave dofs in solution vector
mpc_u.backsubstitution(vel)
mpc_p.backsubstitution(p)
# save results
vtx_pr.write(0)
vtx_ve.write(0)
# destroy
A.destroy()
b.destroy()
for Uh_sub in Uh.getNestSubVecs():
Uh_sub.destroy()
Uh.destroy()
ksp.destroy()