Dear all,
I am trying to solve a non linear problem with periodic bcs and a Langrange multiplier using ufl.MixedFunctionSpace.
For that, I am trying to adapt the BlockedNewtonSolver from scifem to include a list of mpc from dolfinx_mpc using dolfinx_mpc.assemble_matrix_nest and dolfinx_mpc.create_matrix_nest.
I am still learning and debugging my code, but i run into an error before the very first iteration, when assemble the first residual:
[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[0]PETSC ERROR: or see https://petsc.org/release/faq/#valgrind and https://petsc.org/release/faq/
[0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run
[0]PETSC ERROR: to get more information on the crash.
[0]PETSC ERROR: Run with -malloc_debug to check if memory corruption is causing the crash.
Abort(59) on node 0 (rank 0 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 59) - process 0
Here is my draft code
from __future__ import annotations
import logging
from packaging.version import parse as _v
from petsc4py import PETSc
import ufl
import dolfinx
import dolfinx_mpc
logger = logging.getLogger(__name__)
import numpy as np
#Original code avaiable at https://github.com/scientificcomputing/scifem/blob/main/src/scifem/solvers.py
def cast_seq_to_nest(dx: PETSc.Vec, x: PETSc.Vec) -> PETSc.Vec:
"""
Cast a sequential Vec (dx) into a VecNest using the structure of x.
Parameters:
dx (PETSc.Vec): Sequential vector to be converted.
x (PETSc.Vec): Nest vector providing the structure.
Returns:
PETSc.Vec: A new VecNest constructed from dx.
"""
if x.getType() != 'nest':
raise ValueError("x must be a VecNest type.")
# Extract sub-vector sizes from x
sub_vecs_x = x.getNestSubVecs()
sizes = [sub_vec.getSize() for sub_vec in sub_vecs_x]
# Validate total size of dx
if dx.getSize() != sum(sizes):
raise ValueError("The size of dx does not match the total size of x's sub-vectors.")
# Split dx into sub-vectors based on sizes
offsets = np.concatenate(([0], np.cumsum(sizes)[:-1])) # Calculate the offsets
sub_vecs_dx = []
for offset, size in zip(offsets, sizes):
sub_vec = PETSc.Vec().createSeq(size)
sub_vec.array = dx.array[offset : offset + size]
sub_vecs_dx.append(sub_vec)
# Create a new VecNest
dx_nest = PETSc.Vec().createNest(sub_vecs_dx)
dx_nest.assemble()
return dx_nest
class BlockedNewtonSolverMPC(dolfinx.cpp.nls.petsc.NewtonSolver):
def __init__(
self,
F: list[ufl.form.Form],
u: list[dolfinx.fem.Function],
mpc: list[dolfinx_mpc.MultiPointConstraint],
bcs: list[dolfinx.fem.DirichletBC] = [],
J: list[list[ufl.form.Form]] | None = None,
):
# Initialize base class
super().__init__(u[0].function_space.mesh.comm)
self._F = dolfinx.fem.form(F)
# Create the Jacobian matrix, dF/du
if J is None:
du = ufl.TrialFunctions(ufl.MixedFunctionSpace(*[ui.function_space for ui in u]))
J = ufl.extract_blocks(sum(ufl.derivative(sum(F), u[i], du[i]) for i in range(len(u))))
self.jacobian = J
self._a = dolfinx.fem.form(J)
self.mpc = mpc
self.u_mpc = [dolfinx.fem.Function(m.function_space) for m in self.mpc]
self._bcs = bcs
self._u = u
self._b = dolfinx_mpc.create_vector_nest(self._F, mpc)
self._J = dolfinx_mpc.create_matrix_nest(self._a, mpc)
self._dx = dolfinx.fem.petsc.create_vector_nest(self._F)
self._x = dolfinx.fem.petsc.create_vector_nest(self._F)
self.setF(self._assemble_residual, self._b)
self.setJ(self._assemble_jacobian, self._J)
self.set_update(self._update_function)
print(type(self._a))
print("Everything set in solver ")
@property
def L(self) -> list[dolfinx.fem.Form]:
"""Compiled linear form (the residual form)"""
return self._F
@property
def a(self) -> list[list[dolfinx.fem.Form]]:
"""Compiled bilinear form (the Jacobian form)"""
return self._a
@property
def u(self):
return self._u
def __del__(self):
self._A.destroy()
self._b.destroy()
self._dx.destroy()
self._x.destroy()
def _assemble_residual(self, x: PETSc.Vec, b: PETSc.Vec) -> None:
"""Assemble the residual F into the vector b.
Args:
x: The vector containing the latest solution
b: Vector to assemble the residual into
"""
print("Assembling residual", b.getArray(),x.getArray())
dolfinx_mpc.assemble_vector_nest(b, self._F, self.mpc)
b_local.array[:] = 0.0
print("Setting done ")
with b.localForm() as b_local:
b_local.set(0.0)
# x0 = [] if x is None else x.getNestSubVecs()
# bcs1 = dolfinx.fem.bcs_by_block(dolfinx.fem.extract_function_spaces(self._F), self._bcs)
for i,(b_sub,x_sub) in enumerate(zip(b.getNestSubVecs(),x.getNestSubVecs())):
suba=self._J
print(x_sub.getSize(),self._J.getSize())
print("self.jacobian[j][i] ")
dolfinx_mpc.apply_lifting(
b_sub,
# self.jacobian,
[self.jacobian[j][i] for j in range(2)],
bcs=[],
constraint=self.mpc[i],
x0=[x_sub],
scale=dolfinx.default_scalar_type(-1.0),
)
b_sub.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
bcs0 = dolfinx.fem.bcs_by_block(dolfinx.fem.extract_function_spaces(self._F), self._bcs)
dolfinx.fem.petsc.set_bc_nest(b, bcs0)
def _assemble_jacobian(self, x: PETSc.Vec, A: PETSc.Mat) -> None:
print("ASSEMBLING Jacobian attempt")
"""Assemble the Jacobian matrix.
Args:
x: The vector containing the latest solution
"""
A.zeroEntries()
dolfinx_mpc.assemble_matrix_nest(A, self._a, self.mpc, self._bcs)
A.assemble()
# print(A.getArray())
print("ASSEMBLING Jacobian done")
def _update_function(self, solver, dx: PETSc.Vec, x: PETSc.Vec):
"""
Update and add MPC TODO
"""
print("Updating...")
newdx = cast_seq_to_nest(dx,x)
for i,(subdx,subx) in enumerate(zip(newdx.getNestSubVecs(),x.getNestSubVecs())):
self.u_mpc[i].x.petsc_vec.array = subx.array_r
self.u_mpc[i].x.petsc_vec.axpy(-1.0, subdx)
self.u_mpc[i].x.petsc_vec.ghostUpdate(
addv=PETSc.InsertMode.INSERT, # type: ignore
mode=PETSc.ScatterMode.FORWARD, # type: ignore
) # type: ignore
self.mpc[i].homogenize(self.u_mpc[i])
self.mpc[i].backsubstitution(self.u_mpc[i])
subx.array = self.u_mpc[i].x.petsc_vec.array_r
subx.ghostUpdate(
addv=PETSc.InsertMode.INSERT, # type: ignore
mode=PETSc.ScatterMode.FORWARD, # type: ignore
)
print("Updating Done")
def solve(self):
print("Starting solver")
"""Solve non-linear problem into function. Returns the number
of iterations and if the solver converged."""
# dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO)
print("SOlving with ", self._x,self._x.getArray())
print("Starting solver")
print("\n")
n, converged = super().solve(self._x)
self._x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
return n, converged
Thank you very much for your help