Hi, happy new year!
I have a problem about the scifem of BlockedNewtionSolver() in dolfinx v0.9. Link: BlockedNewtonSolver() and test_blocked_newton_solver. Code:
petsc_options = {
"pc_factor_mat_solver_type": "mumps",
}
def NewtonSolver_Option(E_u, variables, bcs, entity_maps):
F = list(ufl.extract_blocks(E_u))
blocked_solver = BlockedNewtonSolver(F, variables, bcs=bcs, petsc_options=petsc_options, entity_maps=entity_maps)
blocked_solver.convergence_criterion = "residual"
blocked_solver.rtol = 1e-6
blocked_solver.atol = 1e-6
blocked_solver.max_it = 500
blocked_solver.relaxation_parameter = 1
return blocked_solver
# 引用自scifem库
# https://github.com/scientificcomputing/scifem/blob/main/src/scifem/solvers.py
from typing import Callable
from packaging.version import parse as _v
_alpha_kw: str = "alpha"
class BlockedNewtonSolver(df.cpp.nls.petsc.NewtonSolver):
def __init__(
self,
F: list[ufl.form.Form],
u: list[df.fem.Function],
bcs: list[df.fem.DirichletBC] = [],
J: list[list[ufl.form.Form]] | None = None,
form_compiler_options: dict | None = None,
jit_options: dict | None = None,
petsc_options: dict | None = None,
entity_maps: dict | None = None,
):
"""Initialize solver for solving a non-linear problem using Newton's method.
Args:
F: List of PDE residuals [F_0(u, v_0), F_1(u, v_1), ...]
u: List of unknown functions u=[u_0, u_1, ...]
bcs: List of Dirichlet boundary conditions
J: UFL representation of the Jacobian (Optional)
Note:
If not provided, the Jacobian will be computed using the
assumption that the test functions come from a ``ufl.MixedFunctionSpace``
form_compiler_options: Options used in FFCx
compilation of this form. Run ``ffcx --help`` at the
command line to see all available options.
jit_options: Options used in CFFI JIT compilation of C
code generated by FFCx. See ``python/dolfinx/jit.py``
for all available options. Takes priority over all
other option values.
petsc_options:
Options passed to the PETSc Krylov solver.
entity_maps: Maps used to map entities between different meshes.
Only needed if the forms have not been compiled a priori,
and has coefficients, test, or trial functions that are defined on different meshes.
"""
# Initialize base class
super().__init__(u[0].function_space.mesh.comm)
# Set PETSc options for Krylov solver
prefix = self.krylov_solver.getOptionsPrefix()
if prefix is None:
prefix = ""
if petsc_options is not None:
# Set PETSc options
opts = PETSc.Options()
opts.prefixPush(prefix)
for k, v in petsc_options.items():
opts[k] = v
opts.prefixPop()
self.krylov_solver.setFromOptions()
self._F = df.fem.form(
F,
form_compiler_options=form_compiler_options,
jit_options=jit_options,
entity_maps=entity_maps,
)
# Create the Jacobian matrix, dF/du
if J is None:
if _v(df.__version__) < _v("0.9"):
raise RuntimeError(
"Automatic computation of Jacobian for blocked problem is only"
+ "supported in DOLFINx 0.9 and later"
)
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._a = df.fem.form(
J,
form_compiler_options=form_compiler_options,
jit_options=jit_options,
entity_maps=entity_maps,
)
self._bcs = bcs
self._u = u
self._pre_solve_callback: Callable[["BlockedNewtonSolver"], None] | None = None
self._post_solve_callback: Callable[["BlockedNewtonSolver"], None] | None = None
# Create structures for holding arrays and matrix
self._b = df.fem.petsc.create_vector_block(self._F)
self._J = df.fem.petsc.create_matrix_block(self._a)
self._dx = df.fem.petsc.create_vector_block(self._F)
self._x = df.fem.petsc.create_vector_block(self._F)
self._J.setOptionsPrefix(prefix)
self._J.setFromOptions()
self.setJ(self._assemble_jacobian, self._J)
self.setF(self._assemble_residual, self._b)
self.set_form(self._pre_newton_iteration)
self.set_update(self._update_function)
def set_pre_solve_callback(self, callback: Callable[["BlockedNewtonSolver"], None]):
"""Set a callback function that is called before each Newton iteration."""
self._pre_solve_callback = callback
def set_post_solve_callback(self, callback: Callable[["BlockedNewtonSolver"], None]):
"""Set a callback function that is called after each Newton iteration."""
self._post_solve_callback = callback
@property
def L(self) -> list[df.fem.Form]:
"""Compiled linear form (the residual form)"""
return self._F
@property
def a(self) -> list[list[df.fem.Form]]:
"""Compiled bilinear form (the Jacobian form)"""
return self._a
@property
def u(self):
return self._u
def __del__(self):
self._J.destroy()
self._b.destroy()
self._dx.destroy()
self._x.destroy()
def _pre_newton_iteration(self, x: PETSc.Vec) -> None:
"""Function called before the residual or Jacobian is
computed.
Args:
x: The vector containing the latest solution
"""
if self._pre_solve_callback is not None:
self._pre_solve_callback(self)
# Scatter previous solution `u=[u_0, ..., u_N]` to `x`; the
# blocked version used for lifting
df.cpp.la.petsc.scatter_local_vectors(
x,
[ui.x.petsc_vec.array_r for ui in self._u],
[
(
ui.function_space.dofmap.index_map,
ui.function_space.dofmap.index_map_bs,
)
for ui in self._u
],
)
x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
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
"""
# Assemble F(u_{i-1}) - J(u_D - u_{i-1}) and set du|_bc= u_D - u_{i-1}
with b.localForm() as b_local:
b_local.set(0.0)
df.fem.petsc.assemble_vector_block(
b,
self._F,
self._a,
bcs=self._bcs,
x0=x,
# dolfinx 0.8 compatibility
# this is called 'scale' in 0.8, 'alpha' in 0.9
**{_alpha_kw: -1.0},
)
b.ghostUpdate(PETSc.InsertMode.INSERT_VALUES, PETSc.ScatterMode.FORWARD)
def _assemble_jacobian(self, x: PETSc.Vec, A: PETSc.Mat) -> None:
"""Assemble the Jacobian matrix.
Args:
x: The vector containing the latest solution
"""
# Assemble Jacobian
A.zeroEntries()
df.fem.petsc.assemble_matrix_block(A, self._a, bcs=self._bcs)
A.assemble()
def _update_function(self, solver, dx: PETSc.Vec, x: PETSc.Vec):
if self._post_solve_callback is not None:
self._post_solve_callback(self)
# Update solution
offset_start = 0
for ui in self._u:
Vi = ui.function_space
num_sub_dofs = Vi.dofmap.index_map.size_local * Vi.dofmap.index_map_bs
ui.x.petsc_vec.array_w[:num_sub_dofs] -= (
self.relaxation_parameter * dx.array_r[offset_start : offset_start + num_sub_dofs]
)
ui.x.petsc_vec.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
offset_start += num_sub_dofs
def solve(self):
"""Solve non-linear problem into function. Returns the number
of iterations and if the solver converged."""
df.log.set_log_level(df.log.LogLevel.INFO)
n, converged = super().solve(self._x)
self._x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
return n, converged
In my problem, I need to use the ufl.MixedFunctionSpace() to solve the problem. I completed the program using BlockedNewtionSolver() and tested it with essentially no errors. However, in my problem, the loads are loaded progressively with each step of my computation and I need to update the equations and boundary conditions at each computation step. The tests have found that there is an unacceptable huge memory consumption, and I found that it was caused by creating the function at each step (I needed to update the F and bcs parameters in BlockedNewtionSolver() at each step) but the resource could not be released at the end. I try to use BlockedNewtionSolver().__del__() to release but there’s no noticeable effect.
A simple loop structure is as follows:
def Displacement_update(step):
bcs = BC_update(step)
E_u = DisplacementEquation()
return NewtonSolver_Option(E_u, [u_sub1, u_sub2], bcs, entity_maps)
while n <= Load_num:
print('+-' * 20, " Begin of Step {} ".format(n), '-+' * 20, '')
Pre_Process(n)
Iter_num = 1
solver_u = Displacement_update(n)
while Iter_num <= Iter_max:
print("> Solver: Step {} Iteration {}".format(n, Iter_num))
print("> Solver: PETSc Krylov Solver & PETSc SNES Solver")
solver_u.solve() # displacement
# Here is the SNES solver for d_sub1 and d_sub2
# Here is the computation of d_czm
# Here are the residuals for d_sub1, d_sub2 and d_czm
if (d_czm_res <= Iter_tol) and (d_sub1_error <= Iter_tol) and (d_sub2_error <= Iter_tol):
break
Iter_num += 1
if Iter_num > Iter_max:
NonConvergence_list.append([n, d_sub1_error, d_sub2_error, d_czm_res])
print('>' * 20, ' Not Converge ', '<' * 20)
Post_Process(n, u_sub1, u_sub2, d_sub1, d_sub2)
print('+-' * 20, " End of Step {} ".format(n), '-+' * 20, '\n\n')
solver_u.__del__()
n += 1