Convergence difficulties with 3D periodic BCs

General
, , ,
1

Hello everyone,

Following my attempts of mixing dolfinx, dolfinx_mpc and scifem packages all together, I came up with some issues regarding the convergence of the problem. To this end, I created a benchmark with a MWE, i.e., a 3D periodic poisson problem, to mimic the convergence issue that I encountered. In this benchmark, I compare the imposition of 3D periodic BCs and dirichlet BCs on master surfaces to the imposition of dirichlet BCs on all faces (w.o. periodic BCs). This is the script:

from __future__ import annotations
from typing import Union
from mpi4py import MPI
import dolfinx.fem as fem
import numpy as np
from dolfinx import default_scalar_type, plot
from dolfinx.common import Timer
from dolfinx.fem import Function, assemble_scalar, form
from dolfinx.mesh import CellType, create_unit_cube, locate_entities_boundary, meshtags, exterior_facet_indices
from numpy.typing import NDArray
from ufl import (
    SpatialCoordinate,
    TestFunction,
    TrialFunction,
    derivative,
    dx,
    grad,
    inner
)
from dolfinx_mpc import NonlinearProblem, LinearProblem, MultiPointConstraint

def benchmark_periodic3D(celltype: CellType, N: int, periodic: bool):

    # Create mesh and finite element
    mesh = create_unit_cube(MPI.COMM_WORLD, N, N, N, celltype)
    tdim = mesh.topology.dim
    V = fem.functionspace(mesh, ("Lagrange", 2))
    tol = float(5e2 * np.finfo(default_scalar_type).resolution)
    print(f'\nTolerance: {tol}, Celltype: {celltype}, Nº cells: {mesh.topology.index_map(mesh.topology.dim).size_local}')
    
    def dirichlet_boundary_total(x: NDArray[Union[np.float32, np.float64]]) -> NDArray[np.bool_]:
        return np.logical_or(
            np.logical_or(
                np.logical_or(np.isclose(x[0], 0.0, atol=tol), np.isclose(x[0], 1.0, atol=tol)),
                np.logical_or(np.isclose(x[1], 0.0, atol=tol), np.isclose(x[1], 1.0, atol=tol))
            ),
            np.logical_or(np.isclose(x[2], 0.0, atol=tol), np.isclose(x[2], 1.0, atol=tol))
        )
    
    def dirichlet_boundary_periodic(x: NDArray[Union[np.float32, np.float64]]) -> NDArray[np.bool_]:
        return np.logical_or(
            np.logical_or(np.isclose(x[0], 0.0, atol=tol), np.isclose(x[1], 0.0, atol=tol)),
            np.isclose(x[2], 0.0, atol=tol)
        )

    # Create Dirichlet boundary condition
    u_L = Function(V)
    u_L.interpolate(lambda x: -1.0/6.0 * ((x[0]-0.5)**2 + (x[1]-0.5)**2 + (x[2]-0.5)**2))
    dirichletboundary = dirichlet_boundary_periodic if periodic else dirichlet_boundary_total
    geometrical_dofs = fem.locate_dofs_geometrical(V, dirichletboundary)
    bc = fem.dirichletbc(u_L, geometrical_dofs)#, V)
    bcs = [bc]

    def PeriodicBoundary(axis):
        """
        Full surface minus dofs constrained by BCs
        """
        return lambda x: np.isclose(x[axis], 1.0, atol=tol)

    def create_periodic_map(axis):
        def periodic_relation(x):
            out_x = np.zeros(x.shape)
            for i in range(3):
                out_x[i] = x[i] - 1.0 if i == axis else x[i]
            return out_x
        return periodic_relation

    with Timer("~~Periodic: Compute mpc condition"):
        mpc = MultiPointConstraint(V)
        for axis in range(3):
            facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, PeriodicBoundary(axis))
            arg_sort = np.argsort(facets)
            mt = meshtags(mesh, mesh.topology.dim - 1, facets[arg_sort], np.full(len(facets), axis+1, dtype=np.int32))
            mpc.create_periodic_constraint_topological(V, mt, axis+1, create_periodic_map(axis), bcs)
        mpc.finalize()

    # Define variational problem
    u = Function(mpc.function_space, name="Solution_u")
    v = TestFunction(V)

    x = SpatialCoordinate(mesh)
    f = 1.0
    F = inner(grad(u), grad(v)) * dx - inner(f, v) * dx
    J = derivative(F, u, TrialFunction(V))

    u_ = TrialFunction(V)
    a = inner(grad(u_), grad(v)) * dx 
    L = inner(f, v) * dx

    u_soln = -1.0/6.0 * ((x[0]-0.5)**2 + (x[1]-0.5)**2 + (x[2]-0.5)**2)

    # Sanity check that the MPC class has some constraints to impose
    num_slaves_global = mesh.comm.allreduce(len(mpc.slaves), op=MPI.SUM)
    num_masters_global = mesh.comm.allreduce(len(mpc.masters.array), op=MPI.SUM)
    assert (not periodic) | ((num_slaves_global > 0) & (num_masters_global > 0))
    print(f'Periodicity: {periodic}, Nº slaves: {num_slaves_global}, Nº masters: {num_masters_global}')

    tol = np.finfo(u.x.array.dtype).resolution
    petsc_options_nonlinear = {
        "snes_type": "ksponly",
        "ksp_type": "preonly",
        "pc_type": "lu",
        "pc_factor_mat_solver_type": "mumps",
        "snes_atol": tol,
        "snes_rtol": tol,
        "ksp_error_if_not_converged": True,
        "snes_error_if_not_converged": True,
        "snes_monitor": None,
    }

    petsc_options_linear = {"ksp_type": "preonly", "pc_type": "lu"}

    print('LINEAR:')
    problem = LinearProblem(a, L, mpc, bcs=bcs, petsc_options=petsc_options_linear)
    u_h = problem.solve()
    l2_error_local = assemble_scalar(form((u_h - u_soln) ** 2 * dx))
    l2_error_global = mesh.comm.allreduce(l2_error_local, op=MPI.SUM)
    print(f' - Error L2 local: {l2_error_local}')
    print(f' - Error L2 global: {l2_error_global}')

    print('NONLINEAR:')
    problem = NonlinearProblem(F, u, mpc=mpc, bcs=bcs, J=J, petsc_options=petsc_options_nonlinear)
    _, converged, num_iterations = problem.solve()
    print(f" - Converged: {converged}, iterations: {num_iterations}")
    l2_error_local = assemble_scalar(form((u - u_soln) ** 2 * dx))
    l2_error_global = mesh.comm.allreduce(l2_error_local, op=MPI.SUM)
    print(f' - Error L2 local: {l2_error_local}')
    print(f' - Error L2 global: {l2_error_global}')

if __name__ == "__main__":
    celltypes = [CellType.hexahedron, CellType.tetrahedron]
    Ne = [8, 12]
    for celltype, N in zip(celltypes, Ne):
        for periodic in [False, True]:
            benchmark_periodic3D(celltype, N, periodic)
    print('\n')

And this is the output obtained:

Tolerance: 5e-13, Celltype: CellType.hexahedron, Nº cells: 512
Periodicity: False, Nº slaves: 0, Nº masters: 0
LINEAR:
 - Error L2 local: 7.650468866454569e-32
 - Error L2 global: 7.650468866454569e-32
NONLINEAR:
0 SNES Function norm 2.837973414409e+00
1 SNES Function norm 5.700934760166e-16
 - Converged: 5, iterations: 1
 - Error L2 local: 9.398975198597027e-32
 - Error L2 global: 9.398975198597027e-32

Tolerance: 5e-13, Celltype: CellType.hexahedron, Nº cells: 512
Periodicity: True, Nº slaves: 768, Nº masters: 768
LINEAR:
 - Error L2 local: 3.342374098344129e-05
 - Error L2 global: 3.342374098344129e-05
NONLINEAR:
0 SNES Function norm 2.117728887619e+00
1 SNES Function norm 2.236806980910e-01
 - Converged: 5, iterations: 1
 - Error L2 local: 3.34237409834424e-05
 - Error L2 global: 3.34237409834424e-05

Tolerance: 5e-13, Celltype: CellType.tetrahedron, Nº cells: 10368
Periodicity: False, Nº slaves: 0, Nº masters: 0
LINEAR:
 - Error L2 local: 1.1490127105419507e-32
 - Error L2 global: 1.1490127105419507e-32
NONLINEAR:
0 SNES Function norm 4.230216839918e+00
1 SNES Function norm 3.326986410762e-16
 - Converged: 5, iterations: 1
 - Error L2 local: 1.3933664178418828e-32
 - Error L2 global: 1.3933664178418828e-32

Tolerance: 5e-13, Celltype: CellType.tetrahedron, Nº cells: 10368
Periodicity: True, Nº slaves: 1728, Nº masters: 1728
LINEAR:
 - Error L2 local: 4.449450735450615e-05
 - Error L2 global: 4.449450735450615e-05
NONLINEAR:
0 SNES Function norm 3.104222328289e+00
1 SNES Function norm 4.605991868480e-01
 - Converged: 5, iterations: 1
 - Error L2 local: 4.4494507354511916e-05
 - Error L2 global: 4.4494507354511916e-05

As far as I am concerned, the multipointconstraint does not create canonical equations for master-slave relations, but, instead, simplifies the assembled matrix by substituting the slave dofs relations by the master ones and eliminating the slave dofs during computation. Shouldn’t both types of problems (with and without periodicity) be theoretically analogous? And if so, shouldn’t they reach the same level of numerical precision? Am I missing something?

Thank you very much in advance.

EDIT: I was mistaken, I did not take into account the compatibility condition of the source term, i.e., its volume integral should be 0. I was so focused on having a poisson problem with an analytical solution to compare the periodicity performance with, that I totally slipped. Sorry for the silly post :sad_but_relieved_face:

2

There are several issues in this code.

  1. Trying to use an MPC to map to a DirichletBC:

If the solution is periodic, there should be no DirichletBC, i.e.

    def dirichlet_boundary_periodic(
        x: NDArray[Union[np.float32, np.float64]],
    ) -> NDArray[np.bool_]:
        return np.full_like(x[0], False, dtype=bool)

You need to be careful with the edge mappings. I attach a code that produces a sensible (and consistent) formulation:

from __future__ import annotations
from typing import Union
from mpi4py import MPI
import dolfinx.fem as fem
import numpy as np
from dolfinx import default_scalar_type, plot
from dolfinx.common import Timer
from dolfinx.fem import Function, assemble_scalar, form
from dolfinx.mesh import (
    CellType,
    create_unit_cube,
    locate_entities_boundary,
    meshtags,
    exterior_facet_indices,
)
from numpy.typing import NDArray
from ufl import (
    SpatialCoordinate,
    TestFunction,
    TrialFunction,
    derivative,
    dx,
    grad,
    inner,
    cos,
    sin,
    pi,
    div,
)
from dolfinx_mpc import NonlinearProblem, LinearProblem, MultiPointConstraint


def benchmark_periodic3D(celltype: CellType, N: int, periodic: bool):
    # Create mesh and finite element
    mesh = create_unit_cube(MPI.COMM_WORLD, N, N, N, celltype)
    tdim = mesh.topology.dim
    V = fem.functionspace(mesh, ("Lagrange", 2))
    tol = float(5e1 * np.finfo(default_scalar_type).resolution)
    print(
        f"\nTolerance: {tol}, Celltype: {celltype}, Nº cells: {mesh.topology.index_map(mesh.topology.dim).size_local}"
    )

    def dirichlet_boundary_total(
        x: NDArray[Union[np.float32, np.float64]],
    ) -> NDArray[np.bool_]:
        return np.logical_or(
            np.logical_or(
                np.logical_or(
                    np.isclose(x[0], 0.0, atol=tol), np.isclose(x[0], 1.0, atol=tol)
                ),
                np.logical_or(
                    np.isclose(x[1], 0.0, atol=tol), np.isclose(x[1], 1.0, atol=tol)
                ),
            ),
            np.logical_or(
                np.isclose(x[2], 0.0, atol=tol), np.isclose(x[2], 1.0, atol=tol)
            ),
        )

    def dirichlet_boundary_periodic(
        x: NDArray[Union[np.float32, np.float64]],
    ) -> NDArray[np.bool_]:
        return np.full_like(x[0], False, dtype=bool)

    # Create Dirichlet boundary condition
    def u_ex(x):
        return (
            np.cos(2 * np.pi * x[0])
            * np.cos(4 * np.pi * x[1])
            * np.cos(2 * np.pi * x[2])
        )

    u_L = Function(V)
    u_L.interpolate(u_ex)
    dirichletboundary = (
        dirichlet_boundary_periodic if periodic else dirichlet_boundary_total
    )
    geometrical_dofs = fem.locate_dofs_geometrical(V, dirichletboundary)
    bc = fem.dirichletbc(u_L, geometrical_dofs)  # , V)
    bcs = [bc]
    print(f"Nº Dirichlet dofs: {len(geometrical_dofs)}")

    def PeriodicBoundary(axis):
        """
        Full surface minus dofs constrained by BCs
        """

        def complicated_mpc(x):
            condition = np.isclose(x[axis], 1.0, atol=tol)
            if axis == 0:
                return condition
            elif axis == 1:
                return np.logical_and(
                    condition, np.logical_and(x[0] > tol, x[0] < 1.0 - tol)
                )
            elif axis == 2:
                return np.logical_and(
                    condition,
                    np.logical_and(
                        np.logical_and(x[0] > tol, x[0] < 1.0 - tol),
                        np.logical_and(x[1] > tol, x[1] < 1.0 - tol),
                    ),
                )
            else:
                raise RuntimeError("Axis should be 0, 1, or 2")

        return complicated_mpc

    def create_periodic_map(axis):
        def periodic_relation(x):
            out_x = x.copy()
            out_x[axis] = 0
            return out_x

        return periodic_relation

    with Timer("~~Periodic: Compute mpc condition"):
        mpc = MultiPointConstraint(V)
        for axis in range(3):
            mpc.create_periodic_constraint_geometrical(
                V, PeriodicBoundary(axis), create_periodic_map(axis), bcs, tol=tol
            )
        mpc.finalize()

    # Define variational problem
    u = Function(mpc.function_space, name="Solution_u")
    v = TestFunction(V)

    x = SpatialCoordinate(mesh)
    u_soln = cos(2 * pi * x[0]) * cos(4 * pi * x[1]) * cos(2 * pi * x[2])
    f = -div(grad(u_soln))
    F = inner(grad(u), grad(v)) * dx - inner(f, v) * dx

    u_ = TrialFunction(V)
    a = inner(grad(u_), grad(v)) * dx
    L = inner(f, v) * dx

    # Sanity check that the MPC class has some constraints to impose
    num_slaves_global = mesh.comm.allreduce(len(mpc.slaves), op=MPI.SUM)
    num_masters_global = mesh.comm.allreduce(len(mpc.masters.array), op=MPI.SUM)
    assert (not periodic) | ((num_slaves_global > 0) & (num_masters_global > 0))
    print(
        f"Periodicity: {periodic}, Nº slaves: {num_slaves_global}, Nº masters: {num_masters_global}"
    )

    tol = np.finfo(u.x.array.dtype).resolution
    petsc_options_nonlinear = {
        "snes_type": "ksponly",
        "ksp_type": "preonly",
        "pc_type": "lu",
        "pc_factor_mat_solver_type": "mumps",
        "snes_atol": tol,
        "snes_rtol": tol,
        "ksp_error_if_not_converged": True,
        "snes_error_if_not_converged": True,
        "snes_monitor": None,
    }

    petsc_options_linear = {
        "ksp_type": "preonly",
        "mat_mumps_icntl_24": 1,
        "mat_mumps_icntl_25": 0,
        "pc_type": "lu",
        "pc_factor_mat_solver_type": "mumps",
        "ksp_error_if_not_converged": True,
    }

    problem = LinearProblem(a, L, mpc, bcs=bcs, petsc_options=petsc_options_linear)
    u_h = problem.solve()
    u_avg_analytical = mesh.comm.allreduce(
        assemble_scalar(form(u_soln * dx)), op=MPI.SUM
    )
    u_avg = mesh.comm.allreduce(assemble_scalar(form(u_h * dx)), op=MPI.SUM)
    u_h.x.array[:] -= u_avg - u_avg_analytical
    import dolfinx

    with dolfinx.io.VTXWriter(
        mesh.comm, f"linear_mpc_{periodic=}_{N}_3D.bp", [u_h]
    ) as writer:
        writer.write(0.0)
    l2_error_local = assemble_scalar(form((u_h - u_soln) ** 2 * dx))
    l2_error_global = mesh.comm.allreduce(l2_error_local, op=MPI.SUM)
    print(f" - Error L2 local: {l2_error_local}")
    print(f" - Error L2 global: {l2_error_global}")

    print("NONLINEAR:")
    J = derivative(F, u, TrialFunction(V))
    problem = NonlinearProblem(
        F, u, mpc=mpc, bcs=bcs, J=J, petsc_options=petsc_options_nonlinear
    )
    _, converged, num_iterations = problem.solve()
    print(f" - Converged: {converged}, iterations: {num_iterations}")
    l2_error_local = assemble_scalar(form((u - u_soln) ** 2 * dx))
    l2_error_global = mesh.comm.allreduce(l2_error_local, op=MPI.SUM)
    print(f" - Error L2 local: {l2_error_local}")
    print(f" - Error L2 global: {l2_error_global}")
    with dolfinx.io.VTXWriter(
        mesh.comm, f"nonlinear_mpc_{periodic=}_{N}_3D.bp", [u]
    ) as writer:
        writer.write(0.0)


if __name__ == "__main__":
    celltypes = [CellType.tetrahedron, CellType.hexahedron]
    Ne = [5, 10, 20]
    for celltype, N in zip(celltypes, Ne):
        for periodic in [True, False]:
            benchmark_periodic3D(celltype, N, periodic)
    print("\n")
NonlinearMPCProblem: error with periodic structure and nonlinear material
3

Thank you for you answer @dokken . The code I presented is a “quick and dirty” script, that relies on behind the scene checks.

With respect to the MPC mappings, I checked that the dolfinx_mpc cpp code is quite robust and it does not apply changes when imposing dirichlet BCs to slave dofs (it does seem to yield an empty MPC object, not affecting the solution). I didn’t post the more complex solution for brevity (sorry, my bad).

I appreciate this comment:

If the solution is periodic, there should be no DirichletBC

I thought fixing a point with a specific value in space was numerically equivalent in dolfinx to setting a Lagrange multiplier for the global condition \int_{\Omega}u \textrm{d}\Omega=0. I was mistaken.
Checking your code, I saw that you substract the integral afterwards for fixing the remaing constant of the field (C:=0). Now I see that your approach is the one equivalent to the lagrange multiplier method. For the sake of curiosity, have you experienced any advantages (performance, etc.) on choosing the substracting approach over the lagrange mul. method?