Unexpected results for periodic constraints

[tuderic]

Dear all,

I am working with periodic boundary conditions and this is a follow-up topic of this post. The aim is to set periodic BCs for large meshes, but the solution of my test cases are not periodic, unfortunately.

The following MWE consists of 4 test cases with simple linear Lagrange elements, and the thermal homogenization method is applied.
The result: Applying periodic BCs to a unit cube works (see the Figure;I also tested it with e.g. 200**3 elements on the HPC) and it also works for a box with the length of 20 in each direction, but not for 40 (see the Figure; the element size is 1 for each box mesh):

import dolfinx
import dolfinx_mpc
from dolfinx.fem import FunctionSpace, Function, Constant, assemble_scalar, form
from dolfinx.io import XDMFFile
from dolfinx.mesh import locate_entities_boundary, meshtags, locate_entities


from ufl import grad, dot, lhs, rhs, Measure, TrialFunction, TestFunction, FiniteElement
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
##################################################
n_size_1 = 20
n_size_2 = 40
unit_cube20 = dolfinx.mesh.create_unit_cube(
        MPI.COMM_WORLD, n_size_1, n_size_1, n_size_1, cell_type=dolfinx.mesh.CellType.hexahedron)
unit_cube40 = dolfinx.mesh.create_unit_cube(
        MPI.COMM_WORLD, n_size_2, n_size_2, n_size_2, cell_type=dolfinx.mesh.CellType.hexahedron)
box20 = dolfinx.mesh.create_box(MPI.COMM_WORLD,
                                     [np.array([0, 0, 0]), np.array([n_size_1, n_size_1, n_size_1])],
                                     [n_size_1, n_size_1, n_size_1], cell_type=dolfinx.mesh.CellType.hexahedron)
box40 = dolfinx.mesh.create_box(MPI.COMM_WORLD,
                                     [np.array([0, 0, 0]), np.array([n_size_2, n_size_2, n_size_2])],
                                     [n_size_2, n_size_2, n_size_2], cell_type=dolfinx.mesh.CellType.hexahedron)

gdim = 3         # Geometric dimension of the mesh
fdim = gdim - 1  # facet dimension
def thermal_hom(domain, ij):

    L_max = domain.geometry.x[:, 0].max()
    L_min = domain.geometry.x[:, 0].min()
    radius = L_max/10
    midpoint = L_max/2
    def cells_sphere(x):
        return ((x[0]-midpoint)**2 + (x[1]-midpoint)**2 + (x[2]-midpoint)**2) <= radius**2
    def cells_cube(x):
        return ((x[0]-midpoint)**2 + (x[1]-midpoint)**2 + (x[2]-midpoint)**2) >= radius**2

    cells_0 = locate_entities(domain, domain.topology.dim, cells_sphere)
    cells_1 = locate_entities(domain, domain.topology.dim, cells_cube)

    ##################################################
    ## Mesh generation finished
    ##################################################

    dx = Measure('dx', domain=domain)
    x = ufl.SpatialCoordinate(domain)

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

    P1 = FiniteElement("CG", domain.ufl_cell(), 1)
    V = FunctionSpace(domain, P1)

    # define function, trial function and test function
    T_fluc = TrialFunction(V)
    dT_fluc = TestFunction(V)

    # define material parameter
    S = FunctionSpace(domain, ("DG", 0))
    kappa = Function(S)
    kappa.x.array[:] = 2.5
    kappa.x.array[cells_0] = np.full_like(cells_0, 5.0, dtype=PETSc.ScalarType)
    kappa.x.array[cells_1] = np.full_like(cells_1, 1.0, dtype=PETSc.ScalarType)
    kappa.x.scatter_forward()

    # boundary conditions
    #####################
    bcs = []
    ## 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_max
            return out_x

        return pbc_slave_to_master_map

    def generate_pbc_is_slave(i):
        return lambda x: np.isclose(x[i], L_max)

    def generate_pbc_is_master(i):
        return lambda x: np.isclose(x[i], L_min)

    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[-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[i][idx] = np.nan
            return out_x

        mpc.create_periodic_constraint_topological(V, pbc_meshtags[i],
                                                       pbc_slave_tags[i],
                                                       pbc_slave_to_master_map,
                                                       bcs)

    if len(pbc_directions) > 1:
        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_max
                out_x[dir_j] = x[dir_j] - L_max
                idx = np.logical_and(pbc_is_slave[dir_i](x), pbc_is_slave[dir_j](x))
                out_x[dir_i][~idx] = np.nan
                out_x[dir_j][~idx] = np.nan
                # remove corner DOFs at point (1,1,1)
                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
                return out_x
            return pbc_slave_to_master_map_edges

        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_max
            out_x[1] = x[1] - L_max
            out_x[2] = x[2] - L_max
            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
            return out_x

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

        # mapping slave intersection node (corner) to master intersection node (opposite corner)
        mpc.create_periodic_constraint_topological(V, 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, 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
    #####################
    T_gradient = Constant(domain, PETSc.ScalarType((1,0,0)))

    ## weak formulation
    eqn = - kappa * dot(grad(dot(T_gradient , x) + T_fluc), grad(dT_fluc)) * dx
    a_form = lhs(eqn)
    L_form = rhs(eqn)

    ## solving
    func_sol = Function(V)
    problem = dolfinx_mpc.LinearProblem(a_form, L_form, mpc, bcs=bcs)
    func_sol = problem.solve()

    # write results to file
    file_results1 = XDMFFile(domain.comm, f"results_3D_MWE_v7/T_fluc{ij}.xdmf", "w")
    file_results1.write_mesh(domain)
    file_results1.write_function(func_sol)
    file_results1.close()

if __name__ == "__main__":
    thermal_hom(unit_cube20, 0)
    thermal_hom(unit_cube40, 1)
    thermal_hom(box20, 2)
    thermal_hom(box40, 3)

Does anyone have an idea why the mapping in test case 4 is wrong?
I am using version 0.7.2 and the calculations are done in serial.

[dokken]

It would help if you could reduce the example into something simpler.

For instance, could you start with reducing the code to a periodic condition on only one face of the cube?

Secondly, you should probably add tolerances to np.isclose in both these

@tuderic

pbc_slave_to_master_map

    def generate_pbc_is_slave(i):
        return lambda x: np.isclose(x[i], L_max)

    def generate_pbc_is_master(i):
        return lambda x: np.isclose(x[i], L_min)

Could you start by check that these return similar number of indices for the two cases (unit cube and box of given size)?

[tuderic]

I hope that the following example is simple enough to discover the error.

I think the error occurs in the function

mpc.create_periodic_constraint_topological(V, pbc_meshtags[i],
                                           pbc_slave_tags[i],
                                           pbc_slave_to_master_map,
                                           bcs)

I checked the inputs of these functions: the visualization of the meshtags in Paraview was correct; the pbc_slave_tags give the same number for both meshes (unit cube and box); and the pbc_slave_to_master_map function give the correct number of DOFs that should be mapped.

So I guess the error occurs due to the implementation in C++ (probably again a rounding error?). Do you have an idea?

import dolfinx
import dolfinx_mpc
from dolfinx.fem import FunctionSpace
from dolfinx.mesh import locate_entities_boundary, meshtags

from ufl import FiniteElement
import numpy as np
from mpi4py import MPI

comm = MPI.COMM_WORLD
rank = MPI.COMM_WORLD.rank

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

n_size = 40     # number of elements in each direction
mesh1 = dolfinx.mesh.create_unit_cube(
        MPI.COMM_WORLD, n_size, n_size, n_size, cell_type=dolfinx.mesh.CellType.hexahedron)
mesh2 = dolfinx.mesh.create_box(MPI.COMM_WORLD,
                                     [np.array([0, 0, 0]), np.array([n_size, n_size, n_size])],
                                     [n_size, n_size, n_size], cell_type=dolfinx.mesh.CellType.hexahedron)

def mapping(domain):

    L_max = domain.geometry.x[:, 0].max()
    L_min = domain.geometry.x[:, 0].min()

    P1 = FiniteElement("CG", domain.ufl_cell(), 1)
    V = FunctionSpace(domain, P1)

    # start for periodic boundary conditions
    bcs = []
    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_max
            return out_x

        return pbc_slave_to_master_map

    def generate_pbc_is_slave(i):
        return lambda x: np.isclose(x[i], L_max, rtol=1e-6)

    def generate_pbc_is_master(i):
        return lambda x: np.isclose(x[i], L_min, rtol=1e-6)

    for i in range(gdim):
        def pbc_is_master_test(x):
            return np.isclose(x[i], L_min, rtol=1e-6)
        
        def pbc_is_slave_test(x):
            return np.isclose(x[i], L_max, rtol=1e-6)

        print(f"Number of indices slaves: {len(locate_entities_boundary(domain, fdim, pbc_is_slave_test))}")
        print(f"Number of indices masters: {len(locate_entities_boundary(domain, fdim, pbc_is_master_test))}")

        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)
    num_mapping_dofs = 0

    for i in range(1):

        # slave/master mapping of opposite surfaces (without slave-slave[-slave] intersections)
        def pbc_slave_to_master_map(x):
            nonlocal num_mapping_dofs
            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[i][idx] = np.nan
            print(f"direction: {i}, Number of DOFs for Surface mapping: , {len(out_x[i][~idx])}, Number of Nodes for Surface mapping that are excluded: {len(out_x[i][idx])}, Number of DOFs at one surface total mapping of: {len(out_x[i][~idx]) + len(out_x[i][idx])} =? {(n_size+1)**2}")
            num_mapping_dofs += len(out_x[i][~idx])
            return out_x

        mpc.create_periodic_constraint_topological(V, pbc_meshtags[i],
                                                       pbc_slave_tags[i],
                                                       pbc_slave_to_master_map,
                                                       bcs)

    mpc.finalize()

    print(f"Number of slaves: {len(mpc.slaves)}, total number of mapping dofs through functions: {num_mapping_dofs}")

if __name__ == "__main__":
    mapping(mesh1)
    mapping(mesh2)

[tuderic]

Update:

I tested a more simplified example with mpc.create_periodic_constraint_geometrical:

import dolfinx
import dolfinx_mpc
from mpi4py import MPI
import basix.ufl
import numpy as np

comm = MPI.COMM_WORLD

n_size = 100     # number of elements in each direction
mesh1 = dolfinx.mesh.create_unit_cube(
        MPI.COMM_WORLD, n_size, n_size, n_size, cell_type=dolfinx.mesh.CellType.hexahedron)
mesh2 = dolfinx.mesh.create_box(MPI.COMM_WORLD,
                                     [np.array([0, 0, 0]), np.array([n_size, n_size, n_size])],
                                     [n_size, n_size, n_size], cell_type=dolfinx.mesh.CellType.hexahedron)
# mesh
def mapping(domain):
    L_max = domain.geometry.x[:, 0].max()

    # functionspace
    P1 = basix.ufl.element("Lagrange", domain.topology.cell_name(), 1, gdim=domain.geometry.dim)
    V  = dolfinx.fem.functionspace(domain, P1)

    ## periodic boundary conditions
    num_mapping_dofs = 0
    mpc = dolfinx_mpc.MultiPointConstraint(V)

    def periodic_boundary_y(x):
        return np.isclose(x[1], L_max, atol=1e-6)

    def periodic_relation_y(x):
        nonlocal num_mapping_dofs
        out_x = np.zeros_like(x)
        out_x[0] = x[0]
        out_x[1] = L_max - x[1]
        out_x[2] = x[2]
        num_mapping_dofs += len(out_x[0])
        return out_x

    bcs = []
    mpc.create_periodic_constraint_geometrical(V, periodic_boundary_y, periodic_relation_y, bcs)

    mpc.finalize()

    print(f"Number of slaves: {len(mpc.slaves)}, total number of mapping DOFs through functions: {num_mapping_dofs}")

if __name__ == "__main__":
    mapping(mesh1)
    mapping(mesh2)

Locally, I get the correct result (dolfinx_mpc is used via the docker image):

Number of slaves: 10201, total number of mapping DOFs through functions: 10201
Number of slaves: 10201, total number of mapping DOFs through functions: 10201

However, at the HPC, I get this result:

Number of slaves: 10201, total number of mapping DOFs through functions: 10201
Number of slaves: 201, total number of mapping DOFs through functions: 10201

The installation at the HPC was done via spack and dolfinx_mpc was build on top from source. I changed one parameter in utils.h as discussed here. For both, I am using version 0.7.2.

My questions are:

1. Do you get the same results? I am a little bit confused, as I get different results for the execution on HPC and via docker image.
2. Why does the example here work in contrast to the example above (at least via docker image)? Is there a difference in mpc.create_periodic_constraint_geometrical and mpc.create_periodic_constraint_topological?

I am very grateful for your help. Thanks.

[tuderic]

Update: The unexpected results no longer occur locally anymore in version 0.8.0.

[dokken]

Happy to hear that the issues are resolved!