Unequal results for linear and nonlinear solver with dolfinx_mpc

dolfinx
1

Hi Community,

I have a problem solving nonlinear problems with periodic boundary conditions using dolfinx_mpc. I solved a 2D linear variational problem with a linear and a nonlinear solver and the results are not equal. How is that possible?

I am solving the following Poisson equation: \kappa \nabla^2 (\overline{\nabla T} \cdot \vec{x} + \tilde{T} ) = 0
The weak formulation is
\kappa \nabla (\overline{\nabla T} \cdot \vec{x} + \tilde{T} ) \cdot \nabla \delta v = 0 with zero Neumann boundary conditions.

As you can see in the following picture, the results between linear and nonlinear solver are only qualitatively the same, but not quantitatively.

Comparison
Comparison4849×2538 732 KB

Do you have any advice what is wrong?
I am not sure, but since I am using Docker and version 0.5.0 of dolfinx and dolfinx_mpc, I can only assume that something is wrong in this version. I appreciate any help!
Here is my MWE:

# Load required libraries
#########################
import gmsh
import dolfinx
from dolfinx import fem, io
from dolfinx.io.gmshio import model_to_mesh
import numpy as np
from mpi4py import MPI
import ufl
from petsc4py import PETSc
from typing import List
import dolfinx_mpc

from Nonlinear_assembly import NonlinearMPCProblem, NewtonSolverMPC
# from: https://github.com/jorgensd/dolfinx_mpc/blob/main/python/tests/test_nonlinear_assembly.py


##############################################################################
#                                   Generate mesh                            #
##############################################################################
gmsh.initialize()
comm = MPI.COMM_WORLD
rank = MPI.COMM_WORLD.rank

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                    # dimension

if rank == 0:
    inclusion = gmsh.model.occ.addDisk(L / 2, L / 2, 0, r, r)
    gmsh.model.occ.synchronize()
    matrix = gmsh.model.occ.addRectangle(0,0,0, L, L)
    gmsh.model.occ.synchronize()
    whole_domain = gmsh.model.occ.fragment([(2, matrix)], [(2, inclusion)])
    gmsh.model.occ.synchronize()

    gmsh.model.addPhysicalGroup(gdim, [matrix], tag = 1)
    gmsh.model.addPhysicalGroup(gdim, [inclusion], tag = 2)

    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.05)
    gmsh.model.mesh.field.setNumber(2, "DistMin", 0.025)
    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, _ = model_to_mesh(gmsh.model, comm, rank, gdim=gdim)
gmsh.finalize()

with io.XDMFFile(domain.comm, "MeshWithGmsh_v2.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_meshtags(ct)

dx = ufl.Measure('dx', domain=domain, subdomain_data=ct)
x_coor = ufl.SpatialCoordinate(domain) # coordinates
num_cells = domain.topology.index_map(domain.topology.dim).size_local
##############################################################################
#                          Function space and elements                       #
##############################################################################
P1 = ufl.FiniteElement("CG", domain.ufl_cell(), 1)  
V = fem.FunctionSpace(domain, P1)                                              

### for linear problem:
T_lin  = ufl.TrialFunction(V)                                                       # fluctuation: primary variable linear problem

### for nonlinear problem:
T_nonl  = fem.Function(V)                                                           # fluctuation: primary variable nonlinear problem
dT_nonl  = ufl.TrialFunction(V)  

# testfunction      
v = ufl.TestFunction(V)
##############################################################################
#                         periodic boundary conditions                       #
##############################################################################
boundary_condition: List[str] = ["periodic", "periodic"]
fdim = domain.topology.dim - 1

bcs = []
pbc_directions = []
pbc_slave_tags = []
pbc_is_slave = []
pbc_is_master = []
pbc_meshtags = []
pbc_slave_to_master_maps = []

# Create MultiPointConstraint object
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] - domain.geometry.x[:, i].max()
        return out_x
    return pbc_slave_to_master_map

def generate_pbc_is_slave(i):
    return lambda x: np.isclose(x[i], domain.geometry.x[:, i].max())

def generate_pbc_is_master(i):
    return lambda x: np.isclose(x[i], domain.geometry.x[:, i].min())

for i, bc_type in enumerate(boundary_condition):

    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 = dolfinx.mesh.locate_entities_boundary(domain, fdim, pbc_is_slave[-1])
    arg_sort = np.argsort(facets)
    pbc_meshtags.append(dolfinx.mesh.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]

    mpc.create_periodic_constraint_topological(V, 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) to the master dof at (0, 0)
if len(pbc_directions) > 1:
    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

    mpc.create_periodic_constraint_topological(V, pbc_meshtags[1],
                                               pbc_slave_tags[1],
                                               pbc_slave_to_master_map,
                                               bcs)
mpc.finalize()
##############################################################################
#                           Definition of parameters                         #
##############################################################################
### material parameters
S = fem.FunctionSpace(domain, ("DG", 0))
kappa = fem.Function(S)

for i in range(num_cells):
    if ct.values[i] == 1: # matrix
        kappa.x.array[i] = 1
    elif ct.values[i] == 2: # inclusions
        kappa.x.array[i] = 5

### prescribed Gradient \bar{\nabla T}
T_gradient = fem.Constant(domain, PETSc.ScalarType((10,0)))
##############################################################################
#                                variational problem                         #
##############################################################################
# linear form
eqn_lin  = kappa * ufl.dot(ufl.grad(ufl.dot(T_gradient , x_coor) + T_lin), ufl.grad(v)) * dx

# nonlinear form
eqn_nonl = kappa * ufl.dot(ufl.grad(ufl.dot(T_gradient , x_coor) + T_nonl), ufl.grad(v)) * dx
Jac = ufl.derivative(eqn_nonl, T_nonl, dT_nonl)
##############################################################################
#                                  Solve the problem                         #
##############################################################################
### linear solver
a_form = ufl.lhs(eqn_lin)
L_form = ufl.rhs(eqn_lin)
problem_lin = dolfinx_mpc.LinearProblem(a_form, L_form, mpc, bcs=bcs)
T_lin = problem_lin.solve()

### nonlinear solver
problem_nonl = NonlinearMPCProblem(eqn_nonl, T_nonl, mpc, bcs=bcs, J = Jac)
solver_nonl  = NewtonSolverMPC(comm, problem_nonl, mpc)
solver_nonl.solve(T_nonl)

with io.XDMFFile(domain.comm, "results.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_function(T_lin)

with io.XDMFFile(domain.comm, "results_nonl.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_function(T_nonl)
2

Hello Tuderic,

I am also working with periodic boundary conditions.
I was wondering, whether you were able to solve this problem?

I looked into the latest version of dolfinx and it seemed that the nonlinear class and newton solver did not change in comparison to the version 0.5.0. Maybe it is a bug?

Best regards
Martin

3

Nothing appears wrong to me. The Poisson problem with pure periodic boundary conditions does not have a unique solution, i.e. if T^h satisfies the PDE and the periodic conditions, then it is trivial to show that T^h + C (with C a constant) is also a solution. You would need to add an additional condition to obtain a unique solution, e.g. prescribing the volume average of \tilde{T} or prescribing the value of \tilde{T} at a single point.

With small modifications to your original code, you can see that changing the initial guess causes the nonlinear solver to obtain solutions with different values of the constant C (I have modified only the solution portion of your code). Note that you should use Xdmf3ReaderT to read the file results_nonl.xdmf into ParaView.

##############################################################################
#                                  Solve the problem                         #
##############################################################################
### linear solver
a_form = ufl.lhs(eqn_lin)
L_form = ufl.rhs(eqn_lin)
problem_lin = dolfinx_mpc.LinearProblem(a_form, L_form, mpc, bcs=bcs)
T_lin = problem_lin.solve()

with io.XDMFFile(domain.comm, "results.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_function(T_lin)

### nonlinear solver
problem_nonl = NonlinearMPCProblem(eqn_nonl, T_nonl, mpc, bcs=bcs, J = Jac)
solver_nonl  = NewtonSolverMPC(comm, problem_nonl, mpc)

with io.XDMFFile(domain.comm, "results_nonl.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)

    for i in range(0, 5):
        T_nonl.x.array[:] = float(i)
        solver_nonl.solve(T_nonl)
        xdmf.write_function(T_nonl, i)
3 Likes
4

Thank you very much for your answer. I really appreciate it and it helps me a lot.

However, there is still a problem since I intially wanted to prescribe the value of \tilde{T} at a single point. When I add the lines

### fix \tilde{T} on one point
def boundary_node(x):
    return np.logical_and( np.isclose(x[0], 0.0), np.isclose(x[1], 0.0))
dofs_node = fem.locate_dofs_geometrical(V, boundary_node)
bc_node = fem.dirichletbc(PETSc.ScalarType(0.0), dofs_node, V)

bcs.append(bc_node)

after the mpc is finalized, solving the linear problem still works, but the nonlinear solver no longer converges. Can you tell me if there is a reason for this as well?

5

You would need to add bc_node to bcs before finalizing mpc, as the MultiPointConstraint removes all Dirichlet dofs from the slave-to-master mapping. This is the purpose of providing bcs as the last argument to mpc.create_periodic_constraint_topological. For this problem, because of the periodicity, you would want to apply the Dirichlet condition to all four dofs on the corners of the domain. See the following, which appears to give me the desired solution for the linear and nonlinear solver:

# Load required libraries
#########################
import gmsh
import dolfinx
from dolfinx import fem, io
from dolfinx.io.gmshio import model_to_mesh
import numpy as np
from mpi4py import MPI
import ufl
from petsc4py import PETSc
from typing import List
import dolfinx_mpc

from Nonlinear_assembly import NonlinearMPCProblem, NewtonSolverMPC
# from: https://github.com/jorgensd/dolfinx_mpc/blob/main/python/tests/test_nonlinear_assembly.py


##############################################################################
#                                   Generate mesh                            #
##############################################################################
gmsh.initialize()
comm = MPI.COMM_WORLD
rank = MPI.COMM_WORLD.rank

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                    # dimension

if rank == 0:
    inclusion = gmsh.model.occ.addDisk(L / 2, L / 2, 0, r, r)
    gmsh.model.occ.synchronize()
    matrix = gmsh.model.occ.addRectangle(0,0,0, L, L)
    gmsh.model.occ.synchronize()
    whole_domain = gmsh.model.occ.fragment([(2, matrix)], [(2, inclusion)])
    gmsh.model.occ.synchronize()

    gmsh.model.addPhysicalGroup(gdim, [matrix], tag = 1)
    gmsh.model.addPhysicalGroup(gdim, [inclusion], tag = 2)

    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.05)
    gmsh.model.mesh.field.setNumber(2, "DistMin", 0.025)
    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, _ = model_to_mesh(gmsh.model, comm, rank, gdim=gdim)
gmsh.finalize()

with io.XDMFFile(domain.comm, "MeshWithGmsh_v2.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_meshtags(ct)

dx = ufl.Measure('dx', domain=domain, subdomain_data=ct)
x_coor = ufl.SpatialCoordinate(domain) # coordinates
num_cells = domain.topology.index_map(domain.topology.dim).size_local
##############################################################################
#                          Function space and elements                       #
##############################################################################
P1 = ufl.FiniteElement("CG", domain.ufl_cell(), 1)  
V = fem.FunctionSpace(domain, P1)                                              

### for linear problem:
T_lin  = ufl.TrialFunction(V)                                                       # fluctuation: primary variable linear problem

### for nonlinear problem:
T_nonl  = fem.Function(V)                                                           # fluctuation: primary variable nonlinear problem
dT_nonl  = ufl.TrialFunction(V)  

# testfunction      
v = ufl.TestFunction(V)
##############################################################################
#                         periodic boundary conditions                       #
##############################################################################
boundary_condition: List[str] = ["periodic", "periodic"]
fdim = domain.topology.dim - 1

bcs = []
### fix \tilde{T} on one point
geom_max = [domain.geometry.x[:, i].max() for i in range(domain.geometry.dim)]
geom_min = [domain.geometry.x[:, i].min() for i in range(domain.geometry.dim)]
def boundary_node(x):
    return np.logical_or(
        np.logical_and(
            np.isclose(x[0], geom_min[0]),
            np.logical_or(
                np.isclose(x[1], geom_min[1]),
                np.isclose(x[1], geom_max[1]),
            )
        ),
        np.logical_and(
            np.isclose(x[0], geom_max[0]),
            np.logical_or(
                np.isclose(x[1], geom_min[1]),
                np.isclose(x[1], geom_max[1]),
            )
        ),
    )
dofs_node = fem.locate_dofs_geometrical(V, boundary_node)
print(dofs_node)
bc_node = fem.dirichletbc(PETSc.ScalarType(0.0), dofs_node, V)
bcs = [bc_node]

pbc_directions = []
pbc_slave_tags = []
pbc_is_slave = []
pbc_is_master = []
pbc_meshtags = []
pbc_slave_to_master_maps = []

# Create MultiPointConstraint object
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] - domain.geometry.x[:, i].max()
        return out_x
    return pbc_slave_to_master_map

def generate_pbc_is_slave(i):
    return lambda x: np.isclose(x[i], domain.geometry.x[:, i].max())

def generate_pbc_is_master(i):
    return lambda x: np.isclose(x[i], domain.geometry.x[:, i].min())

for i, bc_type in enumerate(boundary_condition):

    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 = dolfinx.mesh.locate_entities_boundary(domain, fdim, pbc_is_slave[-1])
    arg_sort = np.argsort(facets)
    pbc_meshtags.append(dolfinx.mesh.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]

    mpc.create_periodic_constraint_topological(V, 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) to the master dof at (0, 0)
# if len(pbc_directions) > 1:
    # 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

    # mpc.create_periodic_constraint_topological(V, pbc_meshtags[1],
                                               # pbc_slave_tags[1],
                                               # pbc_slave_to_master_map,
                                               # bcs)
mpc.finalize()
##############################################################################
#                           Definition of parameters                         #
##############################################################################
### material parameters
S = fem.FunctionSpace(domain, ("DG", 0))
kappa = fem.Function(S)

for i in range(num_cells):
    if ct.values[i] == 1: # matrix
        kappa.x.array[i] = 1
    elif ct.values[i] == 2: # inclusions
        kappa.x.array[i] = 5

### prescribed Gradient \bar{\nabla T}
T_gradient = fem.Constant(domain, PETSc.ScalarType((10,0)))
##############################################################################
#                                variational problem                         #
##############################################################################
# linear form
eqn_lin  = kappa * ufl.dot(ufl.grad(ufl.dot(T_gradient , x_coor) + T_lin), ufl.grad(v)) * dx

# nonlinear form
eqn_nonl = kappa * ufl.dot(ufl.grad(ufl.dot(T_gradient , x_coor) + T_nonl), ufl.grad(v)) * dx
Jac = ufl.derivative(eqn_nonl, T_nonl, dT_nonl)
##############################################################################
#                                  Solve the problem                         #
##############################################################################
### linear solver
a_form = ufl.lhs(eqn_lin)
L_form = ufl.rhs(eqn_lin)
problem_lin = dolfinx_mpc.LinearProblem(a_form, L_form, mpc, bcs=bcs)
T_lin = problem_lin.solve()

with io.XDMFFile(domain.comm, "results.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_function(T_lin)

### nonlinear solver
problem_nonl = NonlinearMPCProblem(eqn_nonl, T_nonl, mpc, bcs=bcs, J = Jac)
solver_nonl  = NewtonSolverMPC(comm, problem_nonl, mpc)
solver_nonl.solve(T_nonl)

with io.XDMFFile(domain.comm, "results_nonl.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_function(T_nonl)
3 Likes
6

This ist really great. Thank you for the explanation!

1 Like