Thanks to both of you for giving me the advice to create a simple example first. I try to provide a minimal working example:
In this example, we have a quadratic domain with one inclusion in the center. These two subdomains have different material parameters.
I want to solve the steady state heat equation (which is actually in this case the poisson equation). However, in the context of computational homogenization, we split the primary variable (the temperature) into a constant part (a prescribed temperature gradient * x) and a fluctuation part, which is the new primary variable. We can write:
\kappa \nabla^2(\overline{\nabla T} \cdot x + \widetilde{T}) = 0
I want to use periodic boundary conditions (function modified from here). Therefore, I adapted this function. When I solve this linear problem with multipoint constraints, the result seems to be right. However, I tried to solve the same problem with the nonlinear MPC classes. I do not understand two things:
- When I try to use a boundary condition (in addition to periodic boundary conditions) in order to fix one degree of freedom (as you prevent “rigid body movements” in the mechanical case), the solver diverges. This behavior cannot be seen for the linear solver.
- When I do not use an additional boundary condition for one DOF, the nonlinear solver converges with one iteration. This should be the case for linear problems. However, the results between the nonlinear and linear solver are totally different.
The code is provided below. The mesh files can be found here.
# Load required libraries
#########################
import dolfinx
from dolfinx import fem, mesh, io, log
import dolfinx_mpc
import numpy as np
from mpi4py import MPI
import ufl
from petsc4py import PETSc
from typing import List
##############################################################################
# Define mesh #
##############################################################################
comm = MPI.COMM_WORLD
rank = comm.rank
meshFile = "2D_OneInclusion_v4" # mesh file name
with io.XDMFFile(comm, "mesh/"+meshFile+"_mesh.xdmf", "r") as xdmf:
domain = xdmf.read_mesh(name="Grid") # reads the mesh
ct = xdmf.read_meshtags(domain, name="Grid") # cell tags
domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim-1)
with io.XDMFFile(comm, "mesh/"+meshFile+"_facets.xdmf", "r") as xdmf:
ft = xdmf.read_meshtags(domain, name="Grid") # facet tags
num_cells = domain.topology.index_map(domain.topology.dim).size_local # number of cells
dx = ufl.Measure('dx', domain=domain, subdomain_data=ct)
x_coor = ufl.SpatialCoordinate(domain)
##############################################################################
# Function space #
##############################################################################
# Elementformulierung
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)
##############################################################################
# Randbedinungen #
##############################################################################
### periodic boundary conditions
if(1):
mpc, bcs = dirichlet_and_periodic_bcs(domain, V, ["periodic", "periodic"])
### Gradient
T_gradient = fem.Constant(domain, PETSc.ScalarType((1,0)))
### fix fluctuation on one point
def boundary_node(x):
return np.logical_and( np.isclose(x[0], domain.geometry.x.min()), np.isclose(x[1], domain.geometry.x.min()))
dofs_node = fem.locate_dofs_geometrical(V, boundary_node)
bc_node = fem.dirichletbc(PETSc.ScalarType(0.0), dofs_node, V)
# this cause divergence for nonlinear solving ... Why?
if(1):
bcs.append(bc_node)
#############################################################################
# Definition von Parametern #
##############################################################################
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
##############################################################################
# Variationsproblem #
##############################################################################
# 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)
##############################################################################
# Lösungsschleife #
##############################################################################
### 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(domain.comm, problem_nonl, mpc)
log.set_log_level(log.LogLevel.INFO)
number, converged = solver_nonl.solve(T_nonl)
assert(converged)
log.set_log_level(log.LogLevel.OFF)
if(0):
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)
This function is used to create periodic boundary conditions. It should work since reasonable results appear for linear solving.
# adapted from https://github.com/jorgensd/dolfinx_mpc/blob/main/python/demos/demo_periodic_gep.py
def dirichlet_and_periodic_bcs(domain, funcspace, boundary_condition: List[str] = ["dirichlet", "periodic"], dbc_value = 0):
"""
- Function to set either dirichlet or periodic boundary conditions.
- This function is was taken from https://github.com/jorgensd/dolfinx_mpc/blob/main/python/demos/demo_periodic_gep.py and adapted.
----------
boundary_condition
First item describes b.c. on {x=0} and {x=1}
Second item describes b.c. on {y=0} and {y=1}
"""
fdim = domain.topology.dim - 1
bcs = []
pbc_directions = []
pbc_slave_tags = []
pbc_is_slave = []
pbc_is_master = []
pbc_meshtags = []
pbc_slave_to_master_maps = []
# number of subspaces in functionspace
N_subspaces = funcspace.num_sub_spaces
# Create MultiPointConstraint object
mpc = dolfinx_mpc.MultiPointConstraint(funcspace)
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())
# Parse boundary conditions
def parse_bcs():
for i, bc_type in enumerate(boundary_condition):
if bc_type == "dirichlet":
def dirichletboundary(x):
return np.logical_or(np.isclose(x[i], domain.geometry.x[:,i].min()), np.isclose(x[i], domain.geometry.x[:,i].max()))
if N_subspaces == 0:
u_bc = fem.Function(functionspace)
u_bc.x.array[:] = dbc_value # value of dirichlet bc needs to be passed into this function!
facets = mesh.locate_entities_boundary(domain, fdim, dirichletboundary)
topological_dofs = fem.locate_dofs_topological(functionspace, fdim, facets)
bcs.append(fem.dirichletbc(u_bc, topological_dofs))
else:
geometrical_dofs = fem.locate_dofs_geometrical(functionspace.collapse()[0], dirichletboundary)
bcs.append(fem.dirichletbc(PETSc.ScalarType(dbc_value), geometrical_dofs, functionspace))
elif bc_type == "periodic":
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 = mesh.locate_entities_boundary(domain, fdim, pbc_is_slave[-1])
arg_sort = np.argsort(facets)
pbc_meshtags.append(mesh.meshtags(domain,
fdim,
facets[arg_sort],
np.full(len(facets), pbc_slave_tags[-1], dtype=np.int32)))
elif bc_type == "neumann":
pass
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(functionspace, 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) to the master dof at (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
mpc.create_periodic_constraint_topological(functionspace, pbc_meshtags[1],
pbc_slave_tags[1],
pbc_slave_to_master_map,
bcs)
if N_subspaces == 0:
functionspace = funcspace
parse_bcs()
else:
for j in range(N_subspaces):
functionspace = funcspace.sub(j)
parse_bcs()
mpc.finalize()
return mpc, bcs
These classes are used for NonlinearMPC problems.
# copied from https://github.com/jorgensd/dolfinx_mpc/blob/main/python/tests/test_nonlinear_assembly.py
class NonlinearMPCProblem(dolfinx.fem.petsc.NonlinearProblem):
def __init__(self, F, u, mpc, bcs=[], J=None, form_compiler_params={},
jit_params={}):
self.mpc = mpc
super().__init__(F, u, bcs=bcs, J=J,
form_compiler_params=form_compiler_params,
jit_params=jit_params)
def F(self, x: PETSc.Vec, F: PETSc.Vec):
with F.localForm() as F_local:
F_local.set(0.0)
dolfinx_mpc.assemble_vector(self._L, self.mpc, b=F)
# Apply boundary condition
dolfinx_mpc.apply_lifting(F, [self._a], bcs=[self.bcs],
constraint=self.mpc, x0=[x], scale=-1.0)
F.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
dolfinx.fem.petsc.set_bc(F, self.bcs, x, -1.0)
def J(self, x: PETSc.Vec, A: PETSc.Mat):
A.zeroEntries()
dolfinx_mpc.assemble_matrix(self._a, self.mpc, bcs=self.bcs, A=A)
A.assemble()
class NewtonSolverMPC(dolfinx.cpp.nls.petsc.NewtonSolver):
def __init__(self, comm: MPI.Intracomm, problem: NonlinearMPCProblem,
mpc: dolfinx_mpc.MultiPointConstraint):
"""A Newton solver for non-linear MPC problems."""
super().__init__(comm)
self.mpc = mpc
self.u_mpc = dolfinx.fem.Function(mpc.function_space)
# Create matrix and vector to be used for assembly of the non-linear
# MPC problem
self._A = dolfinx_mpc.cpp.mpc.create_matrix(
problem.a, mpc._cpp_object)
self._b = dolfinx.cpp.la.petsc.create_vector(
mpc.function_space.dofmap.index_map,
mpc.function_space.dofmap.index_map_bs)
self.setF(problem.F, self._b)
self.setJ(problem.J, self._A)
self.set_form(problem.form)
self.set_update(self.update)
def update(self, solver: dolfinx.nls.petsc.NewtonSolver,
dx: PETSc.Vec, x: PETSc.Vec):
# We need to use a vector created on the MPC's space to update ghosts
self.u_mpc.vector.array = x.array_r
self.u_mpc.vector.axpy(-1.0, dx)
self.u_mpc.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT,
mode=PETSc.ScatterMode.FORWARD)
self.mpc.homogenize(self.u_mpc.vector)
self.mpc.backsubstitution(self.u_mpc.vector)
x.array = self.u_mpc.vector.array_r
x.ghostUpdate(addv=PETSc.InsertMode.INSERT,
mode=PETSc.ScatterMode.FORWARD)
def solve(self, u: dolfinx.fem.Function):
"""Solve non-linear problem into function u. Returns the number
of iterations and if the solver converged."""
n, converged = super().solve(u.vector)
u.x.scatter_forward()
return n, converged
@property
def A(self) -> PETSc.Mat:
"""Jacobian matrix"""
return self._A
@property
def b(self) -> PETSc.Vec:
"""Residual vector"""
return self._b