Stokes equations with periodic boundary conditions

General
1

Hi Community,

I am solving the 2D linear Stokes-equations with periodic boundary conditions (from here) for the pressure and all velocity components in x- and y-direction using dolfinx_mpc (version 0.6.1).

Problem description: The domain is a unit square with an inclusion in the center that has a high dynamic viscosity (rigid particle), so the fluid has to flow around this particle. Therefore, the velocities at the interface between the matrix and the particle are set to zero (non-slip boundary conditions). In order to make the solution unique, the pressures at the corners are set to zero (because of the periodic boundary conditions).

Direct Solver: The solution is verified for computing the linear system of equations with a direct solver. The implementation is correct.

Iterative Solver: I would like to use an iterative solver for upcoming 3D simulations (according to this forum post), but the solution is not periodic anymore, because the solution at two edges are strange (see the attached figure). Since the solution for the direct solver is correct, I assume that there is something wrong in the implementation for the iterative solver. Do you have an idea what the problem may be?

ComparisonDirect_IterativeSolver
ComparisonDirect_IterativeSolver2710×1543 225 KB

Thank you very much!

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, rank, 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)
for i in range(num_cells):
    if ct.values[i] == 1: # matrix
        vis.x.array[i] = 1
    elif ct.values[i] == 2: # inclusions
        vis.x.array[i] = 1e+5
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
            print(out_x[pbc_directions[i]][idx])
            print(idx)
            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

## solver - can be easily switched
solver = ["direct", "iterative"][1]
func_sol = Function(V)

if solver == "direct":
    problem = dolfinx_mpc.LinearProblem(a, L, mpc, bcs=bcs, petsc_options={
        "ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})  # define linear problem
    func_sol = problem.solve()

elif solver == "iterative":

    a_f = form(a)
    L_f = form(L)

    ## iterative solver
    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()

    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()
    ksp.solve(b, func_sol.vector)

ve, pr = func_sol.sub(0).collapse(), func_sol.sub(1).collapse()

# intialize files for the results
file_results1 = XDMFFile(domain.comm, "results_2D_"+solver+"/pressure.xdmf", "w")
file_results1.write_mesh(domain)
file_results4 = XDMFFile(domain.comm, "results_2D_"+solver+"/velocity.xdmf", "w")
file_results4.write_mesh(domain)
file_results1.write_function(pr)
file_results4.write_function(ve)
file_results1.close()
file_results4.close()
2

You are missing the backsubstitution
https://jsdokken.com/dolfinx_mpc/api.html?highlight=backsubstitute#dolfinx_mpc.MultiPointConstraint.backsubstitution
after solving the problem.

Note that you can use LinearProblem with iterative solvers, see for instance: https://github.com/jorgensd/dolfinx_mpc/blob/main/python/demos/demo_contact_2D.py

You can also access the petsc ksp solver if that is more what you’d like to do.

2 Likes
3

Thank you! This works perfectly in 2D.

However, I have some issues in 3D. The problem description is the same, only the circle and the rectangle are replaced by a sphere in the center of a box.

Executing the file with the previous defined solver settings does not work, as

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

only gives infinity values and the iterative solver (used in 2D) does not converge.

I tried another iterative solver (see MWE or the link above) which gives reasonable results at first glance. However, it needs a very long time for one CPU and one corner and one edge at the cube are not periodic. Therefore, I am wondering, if I did any mistake within the definition of 3D periodic bcs (PBCs) or if I need to use different solver settings/preconditioning steps.

For the PBCs: At first, I am mapping the master and slave surfaces of opposite sides without the slave-to-slave intersections. After this, I am mapping the corner of the slave_x-slave_y-slave_z intersection to the master_x-master_y-master_z intersection. And finally, I map all the intersections of two slaves to the intersections of the opposite master edges. Therefore, I think that there is neither an over-constraint nor one constrained missing.

I appreciate any help!

Here is the actual MWE in 3D:

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, MixedElement
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 Start
##################################################
gmsh.initialize()
r = 0.1     # radius of sphere
L = 1.0     # length of box

gdim = 3         # Geometric dimension of the mesh
fdim = gdim - 1  # facet dimension

if comm.rank == 0:
    # Define geometry for RVE
    sphere = gmsh.model.occ.add_sphere(L/2, L/2, L/2, r)
    box = gmsh.model.occ.add_box(0.0, 0.0, 0.0, L, L, L)
    whole_domain = gmsh.model.occ.fragment([(3, box)], [(3, sphere)])
    gmsh.model.occ.synchronize()

    matrix_physical = gmsh.model.addPhysicalGroup(gdim, [box], tag=1)         # tag for matrix
    inclusion_physical = gmsh.model.addPhysicalGroup(gdim, [sphere], tag=2)      # tag for inclusion

    # Get the interface elements between rectangle and circle and tag
    interface_elements = gmsh.model.getBoundary([(gdim, inclusion_physical)])
    interface_physical = gmsh.model.addPhysicalGroup(fdim, (1, 1), 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, 4/3 * np.pi * (r ** 3)):  # 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, "CurvesList", [e[0] for e in edges])
    gmsh.model.mesh.field.setNumber(1, "Sampling", 300)

    gmsh.model.mesh.field.add("Threshold", 2)
    gmsh.model.mesh.field.setNumber(2, "IField", 1)
    gmsh.model.mesh.field.setNumber(2, "LcMin", 0.010)
    gmsh.model.mesh.field.setNumber(2, "LcMax", 0.030)
    gmsh.model.mesh.field.setNumber(2, "DistMin", 0.00)
    gmsh.model.mesh.field.setNumber(2, "DistMax", 0.075)
    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, 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
x = ufl.SpatialCoordinate(domain)

# elements and funcionspace
###########################

# Taylor-Hood elements and Mixed Functionspace
P1 = FiniteElement("CG", domain.ufl_cell(), 1)
P2 = VectorElement("CG", domain.ufl_cell(), 2)
TH = 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)
for i in range(num_cells):
    if ct.values[i] == 1: # matrix
        vis.x.array[i] = 1
    elif ct.values[i] == 2: # inclusions
        vis.x.array[i] = 1e+5
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_v_interf = dolfinx.fem.locate_dofs_topological((V.sub(0), V_coll0)  , ft.dim, ft.find(9))  #dofs of interface tag: 9
bcs_v_interf = dolfinx.fem.dirichletbc(v_bc, dofs_v_interf, V.sub(0)) # dirichlet bc for nonslip bc

## set pressures at the corners to zero
def boundary_node(x):
    return np.logical_or(
        np.logical_and(
            np.isclose(x[0], 0.0),
            np.logical_or(
                np.logical_and(
                    np.isclose(x[1], 0.0),
                    np.logical_or(
                        np.isclose(x[2], 0.0),
                        np.isclose(x[2], L),
                    )
                ),
                np.logical_and(
                    np.isclose(x[1], L),
                    np.logical_or(
                        np.isclose(x[2], 0.0),
                        np.isclose(x[2], L),
                    )
                )
            )
        ),
        np.logical_and(
            np.isclose(x[0], L),
            np.logical_or(
                np.logical_and(
                    np.isclose(x[1], 0.0),
                    np.logical_or(
                        np.isclose(x[2], 0.0),
                        np.isclose(x[2], L),
                    )
                ),
                np.logical_and(
                    np.isclose(x[1], L),
                    np.logical_or(
                        np.isclose(x[2], 0.0),
                        np.isclose(x[2], 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] # combine dirichlet boundary conditions

## 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 + 3) 
    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):
    # slave/master mapping of opposite surfaces (without slave-slave intersections)
 
    def pbc_slave_to_master_map(x):
        out_x = pbc_slave_to_master_maps[i](x)
        idx = np.logical_or(pbc_is_slave[(i + 1) % N_pbc](x),pbc_is_slave[(i + 2) % N_pbc](x))
        out_x[pbc_directions[i]][idx] = np.nan
        return out_x

    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:
    def pbc_slave_to_master_map_corner(x):
        '''
        Maps the slave corner dof [intersection of slave_x, slave_y, slave_z] (1,1,1) to
        master corner dof [master_x, master_y, master_z] (0,0,0)
        '''
        out_x = x.copy()
        out_x[0] = x[0] - L
        out_x[1] = x[1] - L
        out_x[2] = x[2] - L
        idx_corner = np.logical_and(pbc_is_slave[0](x), np.logical_and(pbc_is_slave[1](x), pbc_is_slave[2](x)))
        out_x[0][~idx_corner] = np.nan
        out_x[1][~idx_corner] = np.nan
        out_x[2][~idx_corner] = np.nan
        print(out_x[0][idx_corner])
        return out_x

    def generate_pbc_slave_to_master_map_edges(dir_i, dir_j):
        def pbc_slave_to_master_map_edges(x):
            '''
            Maps the slave edge dofs [intersection of slave_i, slave_j] (i=1,j=1) to
            master corner dof [master_i, master_j] (i=0,j=0)
            (1) map slave_x, slave_y to master_x, master_y: i=0, j=1
            (2) map slave_x, slave_z to master_x, master_z: i=0, j=2
            (3) map slave_y, slave_z to master_y, master_z: i=1, j=2
            '''
            out_x = x.copy()
            out_x[dir_i] = x[dir_i] - L
            out_x[dir_j] = x[dir_j] - L
            idx = np.logical_and(pbc_is_slave[dir_i](x), pbc_is_slave[j](x))
            out_x[dir_i][~idx] = np.nan
            out_x[dir_j][~idx] = np.nan
            idx_corner = np.logical_and(pbc_is_slave[0](x), np.logical_and(pbc_is_slave[1](x), pbc_is_slave[2](x)))
            out_x[dir_i][idx_corner] = np.nan
            out_x[dir_j][idx_corner] = np.nan
            print(len(out_x[dir_i][idx]))
            return out_x
        return pbc_slave_to_master_map_edges

    mapping_slave_intersections = [(0,1),(0,2),(1,2)] # pairs of slave intersections

    for ij in range(gdim):
        # mapping slave intersection node (corner) to master intersection node (opposite corner)
        mpc.create_periodic_constraint_topological(V.sub(0).sub(ij), pbc_meshtags[0],
                                                   pbc_slave_tags[0],
                                                   pbc_slave_to_master_map_corner,
                                                   bcs)

        for inters in mapping_slave_intersections:
            # mapping slave intersections to opposite master intersections
            mpc.create_periodic_constraint_topological(V.sub(0).sub(ij), pbc_meshtags[inters[0]],
                                                       pbc_slave_tags[inters[0]],
                                                       generate_pbc_slave_to_master_map_edges(inters[0], inters[1]),
                                                       bcs)

            if ij == 0: # only one additional constraint
                # mapping slave intersection node (corner) to master intersection node (opposite corner) for PRESSURE
                # not necessary, as the pressure at the corners is set to 0
                # mapping slave intersections to opposite master intersections
                mpc.create_periodic_constraint_topological(V.sub(1), pbc_meshtags[inters[0]],
                                                           pbc_slave_tags[inters[0]],
                                                           generate_pbc_slave_to_master_map_edges(inters[0], inters[1]),
                                                           bcs)

mpc.finalize()

# Variational problem
#####################
f_constant = Constant(domain, PETSc.ScalarType((1,0,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

## solving
func_sol = Function(V)
petsc_options = {"ksp_rtol": 1e-9, "ksp_atol": 1e-9, "pc_type": "gamg",
                 "pc_gamg_type": "agg", "pc_gamg_square_graph": 2,
                 "pc_gamg_threshold": 0.02, "pc_gamg_coarse_eq_limit": 1000, "pc_gamg_sym_graph": True,
                 "mg_levels_ksp_type": "chebyshev", "mg_levels_pc_type": "jacobi",
                 "mg_levels_esteig_ksp_type": "cg"}
problem = dolfinx_mpc.LinearProblem(a, L, mpc, bcs=bcs, petsc_options=petsc_options)
func_sol = problem.solve()

ve = func_sol.sub(0).collapse()
pr = func_sol.sub(1).collapse()

# write results to file
file_results1 = XDMFFile(domain.comm, "results_3D_iterative/pressure.xdmf", "w")
file_results1.write_mesh(domain)
file_results4 = XDMFFile(domain.comm, "results_3D_iterative/velocity.xdmf", "w")
file_results4.write_mesh(domain)
file_results4.write_function(ve)
file_results4.close()
file_results1.write_function(pr)
file_results1.close()
4

I would like to come back to the topic again as I do not get the correct results for parallel computing.

After using the backsubstitution after solving the problem, I get correct results for the computation on one processor. However, using multiple processors, the results are not equal.

For 1 processor:
Average velocity in x direction: 0.009726822408585722
For 2 processors:
Average velocity in x direction: 0.009613288690081394

I think the problem can be related to the backsubstitution as the results for 2 processors looks like this (only the DOFs of one slave surface are obtained from the backsubstitution):

Velocity_field
Velocity_field2956×1873 189 KB

Am I missing something?

Here is the slightly modified 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)

## iterative solver
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)
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
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)

# intialize files for the results
file_results1 = XDMFFile(domain.comm, "results_2D_iterative/pressure.xdmf", "w")
file_results1.write_mesh(domain)
file_results4 = XDMFFile(domain.comm, "results_2D_iterative/velocity.xdmf", "w")
file_results4.write_mesh(domain)
file_results1.write_function(pr)
file_results4.write_function(ve)
file_results1.close()
file_results4.close()