Dear all,
I am wondering whether mpc.backsubstitution also works in parallel. I solve the Stokes equations with periodic boundary conditions in 2D and the solution with a single rank is true. However, errors (only) occur at the slave nodes in parallel. These are the lines after solving the problems:
ksp.solve(b, func_sol.vector)
func_sol.x.scatter_forward()
mpc.backsubstitution(func_sol.vector)
ve, pr = func_sol.sub(0).collapse(), func_sol.sub(1).collapse()
ve.x.scatter_forward()
pr.x.scatter_forward()
Or am I missing something else after the solution of the problem? I appreciate any help! Thank you.
Here is the MWE:
import gmsh
import dolfinx
import dolfinx_mpc
from dolfinx.fem import FunctionSpace, VectorFunctionSpace, TensorFunctionSpace, Function, Constant, Expression, assemble_scalar, form
from dolfinx.io import XDMFFile
from dolfinx.mesh import locate_entities_boundary, meshtags
from dolfinx.io.gmshio import model_to_mesh
from ufl import grad, div, inner, Measure, TrialFunctions, TestFunctions, FiniteElement, VectorElement
import ufl
from petsc4py import PETSc
import numpy as np
from mpi4py import MPI
# MPI initialization
comm = MPI.COMM_WORLD
rank = MPI.COMM_WORLD.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, 0, gdim=gdim)
gmsh.finalize()
#############################################
# Mesh generation finished #
#############################################
dx = Measure('dx', domain=domain, subdomain_data=ct)
num_cells = domain.topology.index_map(domain.topology.dim).size_global
# elements and funcionspace
###########################
# Taylor-Hood elements and Mixed Functionspace
P1 = FiniteElement("CG", domain.ufl_cell(), 1)
P2 = VectorElement("CG", domain.ufl_cell(), 2)
TH = ufl.MixedElement([P2, P1])
V = FunctionSpace(domain, TH)
# define function, trial function and test function
(ve, pr) = TrialFunctions(V)
(dve, dpr) = TestFunctions(V)
# define material parameter
S = FunctionSpace(domain, ("DG", 0))
vis = Function(S)
vis = Function(S)
fluid_visc = ct.find(1)
vis.x.array[fluid_visc] = np.full_like(fluid_visc, 1, dtype=PETSc.ScalarType)
inclusion_visc = ct.find(2)
vis.x.array[inclusion_visc] = np.full_like(inclusion_visc, 1e+5, dtype=PETSc.ScalarType)
vis.x.scatter_forward()
# boundary conditions
#####################
## set all fluid velocities at the interface to zero (no slip boundary condition)
# create collapsed functionspace
V_coll0, _ = V.sub(0).collapse()
v_bc = Function(V_coll0)
v_bc.x.set(0)
v_bc.x.scatter_forward()
# dofs of interface and dirichlet bcs
dofs_v_interf = dolfinx.fem.locate_dofs_topological((V.sub(0), V_coll0) , ft.dim, ft.find(9)) # interface tag: 9
bcs_v_interf = dolfinx.fem.dirichletbc(v_bc, dofs_v_interf, V.sub(0))
## set pressures at the corners to zero
def boundary_node(x):
return np.logical_or(
np.logical_and(
np.isclose(x[0], 0.),
np.logical_or(
np.isclose(x[1], .0),
np.isclose(x[1], L),
)
),
np.logical_and(
np.isclose(x[0], L),
np.logical_or(
np.isclose(x[1], 0.),
np.isclose(x[1], L),
)
),
)
boundary_facets = locate_entities_boundary(domain, 0, boundary_node) # it is zero dimension for a point
dofs_corners = dolfinx.fem.locate_dofs_topological(V.sub(1), 0, boundary_facets) # it is zero dimension for a point
bc_corners = dolfinx.fem.dirichletbc(PETSc.ScalarType(0.0), dofs_corners, V.sub(1))
bcs = [bcs_v_interf, bc_corners]
## periodic boundary conditions
pbc_directions = []
pbc_slave_tags = []
pbc_is_slave = []
pbc_is_master = []
pbc_meshtags = []
pbc_slave_to_master_maps = []
mpc = 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)))
N_pbc = len(pbc_directions)
for i in range(N_pbc):
if N_pbc > 1:
def pbc_slave_to_master_map(x):
out_x = pbc_slave_to_master_maps[i](x)
idx = pbc_is_slave[(i + 1) % N_pbc](x)
out_x[pbc_directions[i]][idx] = np.nan
return out_x
else:
pbc_slave_to_master_map = pbc_slave_to_master_maps[i]
for ij in range(gdim):
mpc.create_periodic_constraint_topological(V.sub(0).sub(ij), pbc_meshtags[i],
pbc_slave_tags[i],
pbc_slave_to_master_map,
bcs)
mpc.create_periodic_constraint_topological(V.sub(1), pbc_meshtags[i],
pbc_slave_tags[i],
pbc_slave_to_master_map,
bcs)
if len(pbc_directions) > 1:
# 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(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.create_periodic_constraint_topological(V.sub(0).sub(ij), pbc_meshtags[1],
pbc_slave_tags[1],
pbc_slave_to_master_map,
bcs)
mpc.create_periodic_constraint_topological(V.sub(1), pbc_meshtags[1],
pbc_slave_tags[1],
pbc_slave_to_master_map,
bcs)
mpc.finalize()
# Variational problem
#####################
f_constant = Constant(domain, PETSc.ScalarType((1,0)))
## weak formulation
a = vis * inner(grad(ve), grad(dve)) * dx + pr * div(dve) * dx + div(ve) * dpr * dx
L = inner(f_constant, dve) * dx
a_f = form(a)
L_f = form(L)
# solution function
func_sol = Function(V)
# assemble matrices
A = dolfinx_mpc.assemble_matrix(a_f, mpc, bcs=bcs)
A.assemble()
b = dolfinx_mpc.assemble_vector(L_f, mpc)
dolfinx_mpc.apply_lifting(b, [a_f], [bcs], mpc)
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
dolfinx.fem.petsc.set_bc(b, bcs)
b.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
# Preconditioner matrix form
Pf = dolfinx.fem.form(ufl.inner(ufl.grad(ve), ufl.grad(dve)) * dx + pr * dpr * dx)
P = dolfinx_mpc.assemble_matrix(Pf, mpc, bcs=bcs)
P.assemble()
# solver settings
ksp = PETSc.KSP().create(domain.comm)
ksp.setOperators(A, P)
ksp.setType(PETSc.KSP.Type.MINRES)
pc = ksp.getPC()
pc.setType("hypre")
pc.setHYPREType("boomeramg")
pc.setUp()
# solving and backsubstitution to obtain the values at the slave degrees of freedom
ksp.solve(b, func_sol.vector)
func_sol.x.scatter_forward()
mpc.backsubstitution(func_sol.vector)
# collapse mixed functionspace
ve, pr = func_sol.sub(0).collapse(), func_sol.sub(1).collapse()
ve.x.scatter_forward()
pr.x.scatter_forward()
# average velocity to see error immediately
v_avg = domain.comm.allreduce(assemble_scalar(form(ve[0] * dx)), op=MPI.SUM)
if rank == 0:
print("Average velocity in x direction:", v_avg)
file_ve = XDMFFile(domain.comm, "results_2D_iterative/velocity.xdmf", "w")
file_ve.write_mesh(domain)
file_ve.write_function(ve)
file_ve.close()
