Coverage for python/src/dolfinx_mpc/problem.py: 77%
260 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-07-02 10:36 +0000
1# -*- coding: utf-8 -*-
2# Copyright (C) 2021-2025 Jørgen S. Dokken
3#
4# This file is part of DOLFINx MPC
5#
6# SPDX-License-Identifier: MIT
7from __future__ import annotations
9import typing
10from functools import partial
11from typing import Iterable, Sequence
13from petsc4py import PETSc
15import dolfinx.fem.petsc
16import numpy as np
17import numpy.typing as npt
18import ufl
19from dolfinx import fem as _fem
20from dolfinx.la.petsc import create_vector
22from dolfinx_mpc.cpp import mpc as _cpp_mpc
24from .assemble_matrix import assemble_matrix, assemble_matrix_nest, create_matrix_nest
25from .assemble_vector import apply_lifting, assemble_vector, assemble_vector_nest, create_vector_nest
26from .multipointconstraint import MultiPointConstraint
29def assemble_jacobian_mpc(
30 u: typing.Union[Sequence[_fem.Function], _fem.Function],
31 jacobian: typing.Union[_fem.Form, typing.Iterable[typing.Iterable[_fem.Form]]],
32 preconditioner: typing.Optional[typing.Union[_fem.Form, typing.Iterable[typing.Iterable[_fem.Form]]]],
33 bcs: typing.Iterable[_fem.DirichletBC],
34 mpc: typing.Union[MultiPointConstraint, Sequence[MultiPointConstraint]],
35 _snes: PETSc.SNES, # type: ignore
36 x: PETSc.Vec, # type: ignore
37 J: PETSc.Mat, # type: ignore
38 P: PETSc.Mat, # type: ignore
39):
40 """Assemble the Jacobian matrix and preconditioner.
42 A function conforming to the interface expected by SNES.setJacobian can
43 be created by fixing the first four arguments:
45 functools.partial(assemble_jacobian, u, jacobian, preconditioner,
46 bcs)
48 Args:
49 u: Function tied to the solution vector within the residual and
50 jacobian
51 jacobian: Form of the Jacobian
52 preconditioner: Form of the preconditioner
53 bcs: List of Dirichlet boundary conditions
54 mpc: The multi point constraint or a sequence of multi point
55 _snes: The solver instance
56 x: The vector containing the point to evaluate at
57 J: Matrix to assemble the Jacobian into
58 P: Matrix to assemble the preconditioner into
59 """
60 # Copy existing soultion into the function used in the residual and
61 # Jacobian
62 try:
63 x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD) # type: ignore
64 except PETSc.Error: # type: ignore
65 for x_sub in x.getNestSubVecs():
66 x_sub.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD) # type: ignore
67 _fem.petsc.assign(x, u)
69 # Assemble Jacobian
70 J.zeroEntries()
71 if J.getType() == "nest":
72 assemble_matrix_nest(J, jacobian, mpc, bcs, diagval=1.0) # type: ignore
73 else:
74 assemble_matrix(jacobian, mpc, bcs, diagval=1.0, A=J) # type: ignore
75 J.assemble()
76 if preconditioner is not None:
77 P.zeroEntries()
78 if P.getType() == "nest":
79 assemble_matrix_nest(P, preconditioner, mpc, bcs, diagval=1.0) # type: ignore
80 else:
81 assemble_matrix(mpc, preconditioner, bcs, diagval=1.0, A=P) # type: ignore
83 P.assemble()
86def assemble_residual_mpc(
87 u: typing.Union[_fem.Function, Sequence[_fem.Function]],
88 residual: typing.Union[_fem.Form, typing.Iterable[_fem.Form]],
89 jacobian: typing.Union[_fem.Form, typing.Iterable[typing.Iterable[_fem.Form]]],
90 bcs: typing.Iterable[_fem.DirichletBC],
91 mpc: typing.Union[MultiPointConstraint, Sequence[MultiPointConstraint]],
92 _snes: PETSc.SNES, # type: ignore
93 x: PETSc.Vec, # type: ignore
94 F: PETSc.Vec, # type: ignore
95):
96 """Assemble the residual into the vector `F`.
98 A function conforming to the interface expected by SNES.setResidual can
99 be created by fixing the first four arguments:
101 functools.partial(assemble_residual, u, jacobian, preconditioner,
102 bcs)
104 Args:
105 u: Function(s) tied to the solution vector within the residual and
106 Jacobian.
107 residual: Form of the residual. It can be a sequence of forms.
108 jacobian: Form of the Jacobian. It can be a nested sequence of
109 forms.
110 bcs: List of Dirichlet boundary conditions.
111 _snes: The solver instance.
112 x: The vector containing the point to evaluate the residual at.
113 F: Vector to assemble the residual into.
114 """
115 # Update input vector before assigning
116 _fem.petsc._ghostUpdate(x, PETSc.InsertMode.INSERT, PETSc.ScatterMode.FORWARD) # type: ignore
118 # Assign the input vector to the unknowns
119 _fem.petsc.assign(x, u)
120 if isinstance(u, Sequence):
121 assert isinstance(mpc, Sequence)
122 for i in range(len(u)):
123 mpc[i].homogenize(u[i])
124 mpc[i].backsubstitution(u[i])
125 else:
126 assert isinstance(u, _fem.Function)
127 assert isinstance(mpc, MultiPointConstraint)
128 mpc.homogenize(u)
129 mpc.backsubstitution(u)
130 # Assemble the residual
131 _fem.petsc._zero_vector(F)
132 if x.getType() == "nest":
133 assemble_vector_nest(F, residual, mpc) # type: ignore
134 else:
135 assert isinstance(residual, _fem.Form)
136 assert isinstance(mpc, MultiPointConstraint)
137 assemble_vector(residual, mpc, F)
139 # Lift vector
140 try:
141 # Nest and blocked lifting
142 bcs1 = _fem.bcs.bcs_by_block(_fem.forms.extract_function_spaces(jacobian, 1), bcs) # type: ignore
143 _fem.petsc._assign_block_data(residual, x) # type: ignore
144 apply_lifting(F, jacobian, bcs=bcs1, constraint=mpc, x0=x, scale=-1.0) # type: ignore
145 _fem.petsc._ghostUpdate(F, PETSc.InsertMode.ADD, PETSc.ScatterMode.REVERSE) # type: ignore
146 bcs0 = _fem.bcs.bcs_by_block(_fem.forms.extract_function_spaces(residual), bcs) # type: ignore
147 _fem.petsc.set_bc(F, bcs0, x0=x, alpha=-1.0)
148 except RuntimeError:
149 # Single form lifting
150 apply_lifting(F, [jacobian], bcs=[bcs], constraint=mpc, x0=[x], scale=-1.0) # type: ignore
151 _fem.petsc._ghostUpdate(F, PETSc.InsertMode.ADD, PETSc.ScatterMode.REVERSE) # type: ignore
152 _fem.petsc.set_bc(F, bcs, x0=x, alpha=-1.0)
153 _fem.petsc._ghostUpdate(F, PETSc.InsertMode.INSERT, PETSc.ScatterMode.FORWARD) # type: ignore
156class NonlinearProblem(dolfinx.fem.petsc.NonlinearProblem):
157 def __init__(
158 self,
159 F: typing.Union[ufl.form.Form, Sequence[ufl.form.Form]],
160 u: typing.Union[_fem.Function, Sequence[_fem.Function]],
161 mpc: typing.Union[MultiPointConstraint, Sequence[MultiPointConstraint]],
162 bcs: typing.Optional[Sequence[_fem.DirichletBC]] = None,
163 J: typing.Optional[typing.Union[ufl.form.Form, Iterable[Iterable[ufl.form.Form]]]] = None,
164 P: typing.Optional[typing.Union[ufl.form.Form, Iterable[Iterable[ufl.form.Form]]]] = None,
165 kind: typing.Optional[typing.Union[str, typing.Iterable[typing.Iterable[str]]]] = None,
166 form_compiler_options: typing.Optional[dict] = None,
167 jit_options: typing.Optional[dict] = None,
168 petsc_options: typing.Optional[dict] = None,
169 entity_maps: typing.Optional[dict[dolfinx.mesh.Mesh, npt.NDArray[np.int32]]] = None,
170 ):
171 """Class for solving nonlinear problems with SNES.
173 Solves problems of the form
174 :math:`F_i(u, v) = 0, i=0,...N\\ \\forall v \\in V` where
175 :math:`u=(u_0,...,u_N), v=(v_0,...,v_N)` using PETSc SNES as the
176 non-linear solver.
178 Note: The deprecated version of this class for use with
179 NewtonSolver has been renamed NewtonSolverNonlinearProblem.
181 Args:
182 F: UFL form(s) of residual :math:`F_i`.
183 u: Function used to define the residual and Jacobian.
184 bcs: Dirichlet boundary conditions.
185 J: UFL form(s) representing the Jacobian
186 :math:`J_ij = dF_i/du_j`.
187 P: UFL form(s) representing the preconditioner.
188 kind: The PETSc matrix type(s) for the Jacobian and
189 preconditioner (``MatType``).
190 See :func:`dolfinx.fem.petsc.create_matrix` for more
191 information.
192 form_compiler_options: Options used in FFCx compilation of all
193 forms. Run ``ffcx --help`` at the command line to see all
194 available options.
195 jit_options: Options used in CFFI JIT compilation of C code
196 generated by FFCx. See ``python/dolfinx/jit.py`` for all
197 available options. Takes priority over all other option
198 values.
199 petsc_options: Options to pass to the PETSc SNES object.
200 entity_maps: If any trial functions, test functions, or
201 coefficients in the form are not defined over the same mesh
202 as the integration domain, ``entity_maps`` must be
203 supplied. For each key (a mesh, different to the
204 integration domain mesh) a map should be provided relating
205 the entities in the integration domain mesh to the entities
206 in the key mesh e.g. for a key-value pair ``(msh, emap)``
207 in ``entity_maps``, ``emap[i]`` is the entity in ``msh``
208 corresponding to entity ``i`` in the integration domain
209 mesh.
210 """
211 # Compile residual and Jacobian forms
212 self._F = _fem.form(
213 F,
214 form_compiler_options=form_compiler_options,
215 jit_options=jit_options,
216 entity_maps=entity_maps,
217 )
219 if J is None:
220 dus = [ufl.TrialFunction(Fi.arguments()[0].ufl_function_space()) for Fi in F]
221 J = _fem.forms.derivative_block(F, u, dus)
223 self._J = _fem.form(
224 J,
225 form_compiler_options=form_compiler_options,
226 jit_options=jit_options,
227 entity_maps=entity_maps,
228 )
230 if P is not None:
231 self._preconditioner = _fem.form(
232 P,
233 form_compiler_options=form_compiler_options,
234 jit_options=jit_options,
235 entity_maps=entity_maps,
236 )
237 else:
238 self._preconditioner = None
240 self._u = u
241 # Set default values if not supplied
242 bcs = [] if bcs is None else bcs
243 self.mpc = mpc
244 # Create PETSc structures for the residual, Jacobian and solution
245 # vector
246 if kind == "nest" or isinstance(kind, Sequence):
247 assert isinstance(mpc, Sequence)
248 assert isinstance(self._J, Sequence)
249 self._A = create_matrix_nest(self._J, mpc)
250 elif kind is None:
251 assert isinstance(mpc, MultiPointConstraint)
252 self._A = _cpp_mpc.create_matrix(self._J._cpp_object, mpc._cpp_object)
253 else:
254 raise ValueError("Unsupported kind for matrix: {}".format(kind))
256 kind = "nest" if self._A.getType() == "nest" else kind
257 if kind == "nest":
258 assert isinstance(mpc, Sequence)
259 assert isinstance(self._F, Sequence)
260 self._b = create_vector_nest(self._F, mpc)
261 self._x = create_vector_nest(self._F, mpc)
262 else:
263 assert isinstance(mpc, MultiPointConstraint)
264 self._b = create_vector(mpc.function_space.dofmap.index_map, mpc.function_space.dofmap.index_map_bs)
265 self._x = create_vector(mpc.function_space.dofmap.index_map, mpc.function_space.dofmap.index_map_bs)
267 # Create PETSc structure for preconditioner if provided
268 if self.preconditioner is not None:
269 if kind == "nest":
270 assert isinstance(self.preconditioner, Sequence)
271 assert isinstance(self.mpc, Sequence)
272 self._P_mat = create_matrix_nest(self.preconditioner, self.mpc)
273 else:
274 assert isinstance(self._P, _fem.Form)
275 self._P_mat = _cpp_mpc.create_matrix(self.preconditioner, kind=kind)
276 else:
277 self._P_mat = None
279 # Create the SNES solver and attach the corresponding Jacobian and
280 # residual computation functions
281 self._snes = PETSc.SNES().create(comm=self.A.comm) # type: ignore
282 self.solver.setJacobian(
283 partial(assemble_jacobian_mpc, u, self.J, self.preconditioner, bcs, mpc), self._A, self.P_mat
284 )
285 self.solver.setFunction(partial(assemble_residual_mpc, u, self.F, self.J, bcs, mpc), self.b)
287 # Set PETSc options
288 problem_prefix = f"dolfinx_nonlinearproblem_{id(self)}"
289 self.solver.setOptionsPrefix(problem_prefix)
290 opts = PETSc.Options() # type: ignore
291 opts.prefixPush(problem_prefix)
293 if petsc_options is not None:
294 for k, v in petsc_options.items():
295 opts[k] = v
296 self.solver.setFromOptions()
298 for k in petsc_options.keys():
299 del opts[k]
301 opts.prefixPop()
303 def solve(self) -> tuple[PETSc.Vec, int, int]: # type: ignore
304 """Solve the problem and update the solution in the problem
305 instance.
307 Returns:
308 The solution, convergence reason and number of iterations.
309 """
311 # Move current iterate into the work array.
312 _fem.petsc.assign(self._u, self.x)
314 # Solve problem
315 self.solver.solve(None, self.x)
317 # Move solution back to function
318 dolfinx.fem.petsc.assign(self.x, self._u)
319 if isinstance(self.mpc, Sequence):
320 for i in range(len(self._u)):
321 self.mpc[i].backsubstitution(self._u[i])
322 self.mpc[i].backsubstitution(self._u[i])
323 else:
324 assert isinstance(self._u, _fem.Function)
325 self.mpc.homogenize(self._u)
326 self.mpc.backsubstitution(self._u)
328 converged_reason = self.solver.getConvergedReason()
329 return self._u, converged_reason, self.solver.getIterationNumber()
332class LinearProblem(dolfinx.fem.petsc.LinearProblem):
333 """
334 Class for solving a linear variational problem with multi point constraints of the form
335 a(u, v) = L(v) for all v using PETSc as a linear algebra backend.
337 Args:
338 a: A bilinear UFL form, the left hand side of the variational problem.
339 L: A linear UFL form, the right hand side of the variational problem.
340 mpc: The multi point constraint.
341 bcs: A list of Dirichlet boundary conditions.
342 u: The solution function. It will be created if not provided. The function has
343 to be based on the functionspace in the mpc, i.e.
345 .. highlight:: python
346 .. code-block:: python
348 u = dolfinx.fem.Function(mpc.function_space)
349 petsc_options: Parameters that is passed to the linear algebra backend PETSc. #type: ignore
350 For available choices for the 'petsc_options' kwarg, see the PETSc-documentation
351 https://www.mcs.anl.gov/petsc/documentation/index.html.
352 form_compiler_options: Parameters used in FFCx compilation of this form. Run `ffcx --help` at
353 the commandline to see all available options. Takes priority over all
354 other parameter values, except for `scalar_type` which is determined by DOLFINx.
355 jit_options: Parameters used in CFFI JIT compilation of C code generated by FFCx.
356 See https://github.com/FEniCS/dolfinx/blob/main/python/dolfinx/jit.py#L22-L37
357 for all available parameters. Takes priority over all other parameter values.
358 P: A preconditioner UFL form.
359 Examples:
360 Example usage:
362 .. highlight:: python
363 .. code-block:: python
365 problem = LinearProblem(a, L, mpc, [bc0, bc1],
366 petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
368 """
370 u: typing.Union[_fem.Function, list[_fem.Function]]
371 _a: typing.Union[_fem.Form, typing.Sequence[typing.Sequence[_fem.Form]]]
372 _L: typing.Union[_fem.Form, typing.Sequence[_fem.Form]]
373 _preconditioner: typing.Optional[typing.Union[_fem.Form, typing.Sequence[typing.Sequence[_fem.Form]]]]
374 _mpc: typing.Union[MultiPointConstraint, typing.Sequence[MultiPointConstraint]]
375 _A: PETSc.Mat # type: ignore
376 _P: typing.Optional[PETSc.Mat] # type: ignore
377 _b: PETSc.Vec # type: ignore
378 _solver: PETSc.KSP # type: ignore
379 _x: PETSc.Vec # type: ignore
380 bcs: typing.List[_fem.DirichletBC]
381 __slots__ = tuple(__annotations__)
383 def __init__(
384 self,
385 a: typing.Union[ufl.Form, typing.Iterable[typing.Iterable[ufl.Form]]],
386 L: typing.Union[ufl.Form, typing.Iterable[ufl.Form]],
387 mpc: typing.Union[MultiPointConstraint, typing.Sequence[MultiPointConstraint]],
388 bcs: typing.Optional[typing.List[_fem.DirichletBC]] = None,
389 u: typing.Optional[typing.Union[_fem.Function, typing.Sequence[_fem.Function]]] = None,
390 petsc_options: typing.Optional[dict] = None,
391 form_compiler_options: typing.Optional[dict] = None,
392 jit_options: typing.Optional[dict] = None,
393 P: typing.Optional[typing.Union[ufl.Form, typing.Iterable[ufl.Form]]] = None,
394 ):
395 # Compile forms
396 form_compiler_options = {} if form_compiler_options is None else form_compiler_options
397 jit_options = {} if jit_options is None else jit_options
398 self._a = _fem.form(a, jit_options=jit_options, form_compiler_options=form_compiler_options)
399 self._L = _fem.form(L, jit_options=jit_options, form_compiler_options=form_compiler_options)
401 self._mpc = mpc
402 # Nest assembly
403 if isinstance(mpc, Sequence):
404 is_nest = True
405 # Sanity check
406 for mpc_i in mpc:
407 if not mpc_i.finalized:
408 raise RuntimeError("The multi point constraint has to be finalized before calling initializer")
409 # Create function containing solution vector
410 else:
411 is_nest = False
412 if not mpc.finalized:
413 raise RuntimeError("The multi point constraint has to be finalized before calling initializer")
415 # Create function(s) containing solution vector(s)
416 if is_nest:
417 if u is None:
418 assert isinstance(self._mpc, Sequence)
419 self.u = [_fem.Function(self._mpc[i].function_space) for i in range(len(self._mpc))]
420 else:
421 assert isinstance(self._mpc, Sequence)
422 assert isinstance(u, Sequence)
423 for i, (mpc_i, u_i) in enumerate(zip(self._mpc, u)):
424 assert isinstance(u_i, _fem.Function)
425 assert isinstance(mpc_i, MultiPointConstraint)
426 if u_i.function_space is not mpc_i.function_space:
427 raise ValueError(
428 "The input function has to be in the function space in the multi-point constraint",
429 "i.e. u = dolfinx.fem.Function(mpc.function_space)",
430 )
431 self.u = list(u)
432 else:
433 if u is None:
434 assert isinstance(self._mpc, MultiPointConstraint)
435 self.u = _fem.Function(self._mpc.function_space)
436 else:
437 assert isinstance(u, _fem.Function)
438 assert isinstance(self._mpc, MultiPointConstraint)
439 if u.function_space is self._mpc.function_space:
440 self.u = u
441 else:
442 raise ValueError(
443 "The input function has to be in the function space in the multi-point constraint",
444 "i.e. u = dolfinx.fem.Function(mpc.function_space)",
445 )
447 # Create MPC matrix and vector
448 self._preconditioner = _fem.form(P, jit_options=jit_options, form_compiler_options=form_compiler_options)
450 if is_nest:
451 assert isinstance(mpc, Sequence)
452 assert isinstance(self._L, Sequence)
453 assert isinstance(self._a, Sequence)
454 self._A = create_matrix_nest(self._a, mpc)
455 self._b = create_vector_nest(self._L, mpc)
456 self._x = create_vector_nest(self._L, mpc)
457 if self._preconditioner is None:
458 self._P_mat = None
459 else:
460 assert isinstance(self._preconditioner, Sequence)
461 self._P_mat = create_matrix_nest(self._preconditioner, mpc)
462 else:
463 assert isinstance(mpc, MultiPointConstraint)
464 assert isinstance(self._L, _fem.Form)
465 assert isinstance(self._a, _fem.Form)
466 self._A = _cpp_mpc.create_matrix(self._a._cpp_object, mpc._cpp_object)
467 self._b = create_vector(mpc.function_space.dofmap.index_map, mpc.function_space.dofmap.index_map_bs)
468 self._x = create_vector(mpc.function_space.dofmap.index_map, mpc.function_space.dofmap.index_map_bs)
469 if self._preconditioner is None:
470 self._P_mat = None
471 else:
472 assert isinstance(self._preconditioner, _fem.Form)
473 self._P_mat = _cpp_mpc.create_matrix(self._preconditioner._cpp_object, mpc._cpp_object)
475 self.bcs = [] if bcs is None else bcs
477 if is_nest:
478 assert isinstance(self.u, Sequence)
479 comm = self.u[0].function_space.mesh.comm
480 else:
481 assert isinstance(self.u, _fem.Function)
482 comm = self.u.function_space.mesh.comm
484 self._solver = PETSc.KSP().create(comm) # type: ignore
485 self._solver.setOperators(self._A, self._P_mat)
487 # Give PETSc solver options a unique prefix
488 solver_prefix = "dolfinx_mpc_solve_{}".format(id(self))
489 self._solver.setOptionsPrefix(solver_prefix)
491 # Set PETSc options
492 opts = PETSc.Options() # type: ignore
493 opts.prefixPush(solver_prefix)
494 if petsc_options is not None:
495 for k, v in petsc_options.items():
496 opts[k] = v
497 opts.prefixPop()
498 self._solver.setFromOptions()
500 def solve(self) -> typing.Union[list[_fem.Function], _fem.Function]:
501 """Solve the problem.
503 Returns:
504 Function containing the solution"""
506 # Assemble lhs
507 self._A.zeroEntries()
508 if self._A.getType() == "nest":
509 assemble_matrix_nest(self._A, self._a, self._mpc, self.bcs, diagval=1.0) # type: ignore
510 else:
511 assert isinstance(self._a, _fem.Form)
512 assemble_matrix(self._a, self._mpc, bcs=self.bcs, A=self._A)
514 self._A.assemble()
515 assert self._A.assembled
517 # Assemble the preconditioner if provided
518 if self._P_mat is not None:
519 self._P_mat.zeroEntries()
520 if self._P_mat.getType() == "nest":
521 assert isinstance(self._preconditioner, typing.Sequence)
522 assemble_matrix_nest(self._P_mat, self._preconditioner, self._mpc, self.bcs) # type: ignore
523 else:
524 assert isinstance(self._preconditioner, _fem.Form)
525 assemble_matrix(self._preconditioner, self._mpc, bcs=self.bcs, A=self._P_mat)
526 self._P_mat.assemble()
528 # Assemble the residual
529 _fem.petsc._zero_vector(self._b)
530 if self._x.getType() == "nest":
531 assemble_vector_nest(self._b, self._L, self._mpc) # type: ignore
532 else:
533 assert isinstance(self._L, _fem.Form)
534 assert isinstance(self._mpc, MultiPointConstraint)
535 assemble_vector(self._L, self._mpc, self._b)
537 # Lift vector
538 try:
539 # Nest and blocked lifting
540 bcs1 = _fem.bcs.bcs_by_block(_fem.forms.extract_function_spaces(self._a, 1), self.bcs) # type: ignore
541 apply_lifting(self._b, self._a, bcs=bcs1, constraint=self._mpc) # type: ignore
542 _fem.petsc._ghostUpdate(self._b, PETSc.InsertMode.ADD, PETSc.ScatterMode.REVERSE) # type: ignore
543 bcs0 = _fem.bcs.bcs_by_block(_fem.forms.extract_function_spaces(self._L), self.bcs) # type: ignore
544 _fem.petsc.set_bc(self._b, bcs0)
545 except RuntimeError:
546 # Single form lifting
547 apply_lifting(self._b, [self._a], bcs=[self.bcs], constraint=self._mpc) # type: ignore
548 _fem.petsc._ghostUpdate(self._b, PETSc.InsertMode.ADD, PETSc.ScatterMode.REVERSE) # type: ignore
549 _fem.petsc.set_bc(self._b, self.bcs)
550 _fem.petsc._ghostUpdate(self._b, PETSc.InsertMode.INSERT, PETSc.ScatterMode.FORWARD) # type: ignore
552 # Solve linear system and update ghost values in the solution
553 self._solver.solve(self._b, self._x)
554 _fem.petsc._ghostUpdate(self._x, PETSc.InsertMode.INSERT, PETSc.ScatterMode.FORWARD) # type: ignore
555 _fem.petsc.assign(self._x, self.u)
557 if isinstance(self.u, Sequence):
558 assert isinstance(self._mpc, Sequence)
559 for i in range(len(self.u)):
560 self._mpc[i].homogenize(self.u[i])
561 self._mpc[i].backsubstitution(self.u[i])
562 else:
563 assert isinstance(self.u, _fem.Function)
564 assert isinstance(self._mpc, MultiPointConstraint)
565 self._mpc.homogenize(self.u)
566 self._mpc.backsubstitution(self.u)
568 return self.u