Coverage for python/src/dolfinx_mpc/problem.py: 78%
263 statements
« prev ^ index » next coverage.py v7.9.2, created at 2025-07-11 21:02 +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 _ghost_update, _zero_vector, 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 _ghost_update(x, PETSc.InsertMode.INSERT, PETSc.ScatterMode.FORWARD) # type: ignore
63 _fem.petsc.assign(x, u)
65 # Assemble Jacobian
66 J.zeroEntries()
67 if J.getType() == "nest":
68 assemble_matrix_nest(J, jacobian, mpc, bcs, diagval=1.0) # type: ignore
69 else:
70 assemble_matrix(jacobian, mpc, bcs, diagval=1.0, A=J) # type: ignore
71 J.assemble()
72 if preconditioner is not None:
73 P.zeroEntries()
74 if P.getType() == "nest":
75 assemble_matrix_nest(P, preconditioner, mpc, bcs, diagval=1.0) # type: ignore
76 else:
77 assemble_matrix(mpc, preconditioner, bcs, diagval=1.0, A=P) # type: ignore
79 P.assemble()
82def assemble_residual_mpc(
83 u: typing.Union[_fem.Function, Sequence[_fem.Function]],
84 residual: typing.Union[_fem.Form, typing.Iterable[_fem.Form]],
85 jacobian: typing.Union[_fem.Form, typing.Iterable[typing.Iterable[_fem.Form]]],
86 bcs: typing.Iterable[_fem.DirichletBC],
87 mpc: typing.Union[MultiPointConstraint, Sequence[MultiPointConstraint]],
88 _snes: PETSc.SNES, # type: ignore
89 x: PETSc.Vec, # type: ignore
90 F: PETSc.Vec, # type: ignore
91):
92 """Assemble the residual into the vector `F`.
94 A function conforming to the interface expected by SNES.setResidual can
95 be created by fixing the first four arguments:
97 functools.partial(assemble_residual, u, jacobian, preconditioner,
98 bcs)
100 Args:
101 u: Function(s) tied to the solution vector within the residual and
102 Jacobian.
103 residual: Form of the residual. It can be a sequence of forms.
104 jacobian: Form of the Jacobian. It can be a nested sequence of
105 forms.
106 bcs: List of Dirichlet boundary conditions.
107 _snes: The solver instance.
108 x: The vector containing the point to evaluate the residual at.
109 F: Vector to assemble the residual into.
110 """
111 # Update input vector before assigning
112 _ghost_update(x, PETSc.InsertMode.INSERT, PETSc.ScatterMode.FORWARD) # type: ignore
113 # Assign the input vector to the unknowns
114 _fem.petsc.assign(x, u)
115 if isinstance(u, Sequence):
116 assert isinstance(mpc, Sequence)
117 for i in range(len(u)):
118 mpc[i].homogenize(u[i])
119 mpc[i].backsubstitution(u[i])
120 else:
121 assert isinstance(u, _fem.Function)
122 assert isinstance(mpc, MultiPointConstraint)
123 mpc.homogenize(u)
124 mpc.backsubstitution(u)
125 # Assemble the residual
126 _zero_vector(F)
127 if x.getType() == "nest":
128 assemble_vector_nest(F, residual, mpc) # type: ignore
129 else:
130 assert isinstance(residual, _fem.Form)
131 assert isinstance(mpc, MultiPointConstraint)
132 assemble_vector(residual, mpc, F)
134 # Lift vector
135 try:
136 # Nest and blocked lifting
137 bcs1 = _fem.bcs.bcs_by_block(_fem.forms.extract_function_spaces(jacobian, 1), bcs) # type: ignore
138 _fem.petsc._assign_block_data(residual, x) # type: ignore
139 apply_lifting(F, jacobian, bcs=bcs1, constraint=mpc, x0=x, scale=-1.0) # type: ignore
140 _ghost_update(F, PETSc.InsertMode.ADD, PETSc.ScatterMode.REVERSE) # type: ignore
141 bcs0 = _fem.bcs.bcs_by_block(_fem.forms.extract_function_spaces(residual), bcs) # type: ignore
142 _fem.petsc.set_bc(F, bcs0, x0=x, alpha=-1.0)
143 except RuntimeError:
144 # Single form lifting
145 apply_lifting(F, [jacobian], bcs=[bcs], constraint=mpc, x0=[x], scale=-1.0) # type: ignore
146 _ghost_update(F, PETSc.InsertMode.ADD, PETSc.ScatterMode.REVERSE) # type: ignore
147 _fem.petsc.set_bc(F, bcs, x0=x, alpha=-1.0)
148 _ghost_update(F, PETSc.InsertMode.INSERT, PETSc.ScatterMode.FORWARD) # type: ignore
151class NonlinearProblem(dolfinx.fem.petsc.NonlinearProblem):
152 def __init__(
153 self,
154 F: typing.Union[ufl.form.Form, Sequence[ufl.form.Form]],
155 u: typing.Union[_fem.Function, Sequence[_fem.Function]],
156 mpc: typing.Union[MultiPointConstraint, Sequence[MultiPointConstraint]],
157 bcs: typing.Optional[Sequence[_fem.DirichletBC]] = None,
158 J: typing.Optional[typing.Union[ufl.form.Form, Iterable[Iterable[ufl.form.Form]]]] = None,
159 P: typing.Optional[typing.Union[ufl.form.Form, Iterable[Iterable[ufl.form.Form]]]] = None,
160 kind: typing.Optional[typing.Union[str, typing.Iterable[typing.Iterable[str]]]] = None,
161 form_compiler_options: typing.Optional[dict] = None,
162 jit_options: typing.Optional[dict] = None,
163 petsc_options: typing.Optional[dict] = None,
164 entity_maps: typing.Optional[dict[dolfinx.mesh.Mesh, npt.NDArray[np.int32]]] = None,
165 ):
166 """Class for solving nonlinear problems with SNES.
168 Solves problems of the form
169 :math:`F_i(u, v) = 0, i=0,...N\\ \\forall v \\in V` where
170 :math:`u=(u_0,...,u_N), v=(v_0,...,v_N)` using PETSc SNES as the
171 non-linear solver.
173 Note: The deprecated version of this class for use with
174 NewtonSolver has been renamed NewtonSolverNonlinearProblem.
176 Args:
177 F: UFL form(s) of residual :math:`F_i`.
178 u: Function used to define the residual and Jacobian.
179 bcs: Dirichlet boundary conditions.
180 J: UFL form(s) representing the Jacobian
181 :math:`J_ij = dF_i/du_j`.
182 P: UFL form(s) representing the preconditioner.
183 kind: The PETSc matrix type(s) for the Jacobian and
184 preconditioner (``MatType``).
185 See :func:`dolfinx.fem.petsc.create_matrix` for more
186 information.
187 form_compiler_options: Options used in FFCx compilation of all
188 forms. Run ``ffcx --help`` at the command line to see all
189 available options.
190 jit_options: Options used in CFFI JIT compilation of C code
191 generated by FFCx. See ``python/dolfinx/jit.py`` for all
192 available options. Takes priority over all other option
193 values.
194 petsc_options: Options to pass to the PETSc SNES object.
195 entity_maps: If any trial functions, test functions, or
196 coefficients in the form are not defined over the same mesh
197 as the integration domain, ``entity_maps`` must be
198 supplied. For each key (a mesh, different to the
199 integration domain mesh) a map should be provided relating
200 the entities in the integration domain mesh to the entities
201 in the key mesh e.g. for a key-value pair ``(msh, emap)``
202 in ``entity_maps``, ``emap[i]`` is the entity in ``msh``
203 corresponding to entity ``i`` in the integration domain
204 mesh.
205 """
206 # Compile residual and Jacobian forms
207 self._F = _fem.form(
208 F,
209 form_compiler_options=form_compiler_options,
210 jit_options=jit_options,
211 entity_maps=entity_maps,
212 )
214 if J is None:
215 dus = [ufl.TrialFunction(Fi.arguments()[0].ufl_function_space()) for Fi in F]
216 J = _fem.forms.derivative_block(F, u, dus)
218 self._J = _fem.form(
219 J,
220 form_compiler_options=form_compiler_options,
221 jit_options=jit_options,
222 entity_maps=entity_maps,
223 )
225 if P is not None:
226 self._preconditioner = _fem.form(
227 P,
228 form_compiler_options=form_compiler_options,
229 jit_options=jit_options,
230 entity_maps=entity_maps,
231 )
232 else:
233 self._preconditioner = None
235 self._u = u
236 # Set default values if not supplied
237 bcs = [] if bcs is None else bcs
238 self.mpc = mpc
239 # Create PETSc structures for the residual, Jacobian and solution
240 # vector
241 if kind == "nest" or isinstance(kind, Sequence):
242 assert isinstance(mpc, Sequence)
243 assert isinstance(self._J, Sequence)
244 self._A = create_matrix_nest(self._J, mpc)
245 elif kind is None:
246 assert isinstance(mpc, MultiPointConstraint)
247 self._A = _cpp_mpc.create_matrix(self._J._cpp_object, mpc._cpp_object)
248 else:
249 raise ValueError("Unsupported kind for matrix: {}".format(kind))
251 kind = "nest" if self._A.getType() == "nest" else kind
252 if kind == "nest":
253 assert isinstance(mpc, Sequence)
254 assert isinstance(self._F, Sequence)
255 self._b = create_vector_nest(self._F, mpc)
256 self._x = create_vector_nest(self._F, mpc)
257 else:
258 assert isinstance(mpc, MultiPointConstraint)
259 self._b = create_vector(mpc.function_space.dofmap.index_map, mpc.function_space.dofmap.index_map_bs)
260 self._x = create_vector(mpc.function_space.dofmap.index_map, mpc.function_space.dofmap.index_map_bs)
262 # Create PETSc structure for preconditioner if provided
263 if self.preconditioner is not None:
264 if kind == "nest":
265 assert isinstance(self.preconditioner, Sequence)
266 assert isinstance(self.mpc, Sequence)
267 self._P_mat = create_matrix_nest(self.preconditioner, self.mpc)
268 else:
269 assert isinstance(self._P, _fem.Form)
270 self._P_mat = _cpp_mpc.create_matrix(self.preconditioner, kind=kind)
271 else:
272 self._P_mat = None
274 # Create the SNES solver and attach the corresponding Jacobian and
275 # residual computation functions
276 self._snes = PETSc.SNES().create(comm=self.A.comm) # type: ignore
277 self.solver.setJacobian(
278 partial(assemble_jacobian_mpc, u, self.J, self.preconditioner, bcs, mpc), self._A, self.P_mat
279 )
280 self.solver.setFunction(partial(assemble_residual_mpc, u, self.F, self.J, bcs, mpc), self.b)
282 # Set PETSc options
283 problem_prefix = f"dolfinx_nonlinearproblem_{id(self)}"
284 self.solver.setOptionsPrefix(problem_prefix)
285 opts = PETSc.Options() # type: ignore
286 opts.prefixPush(problem_prefix)
288 if petsc_options is not None:
289 for k, v in petsc_options.items():
290 opts[k] = v
291 self.solver.setFromOptions()
293 for k in petsc_options.keys():
294 del opts[k]
296 opts.prefixPop()
298 def solve(self) -> tuple[PETSc.Vec, int, int]: # type: ignore
299 """Solve the problem and update the solution in the problem
300 instance.
302 Returns:
303 The solution, convergence reason and number of iterations.
304 """
306 # Move current iterate into the work array.
307 _fem.petsc.assign(self._u, self.x)
309 # Solve problem
310 self.solver.solve(None, self.x)
312 # Move solution back to function
313 dolfinx.fem.petsc.assign(self.x, self._u)
314 if isinstance(self.mpc, Sequence):
315 for i in range(len(self._u)):
316 self.mpc[i].backsubstitution(self._u[i])
317 self.mpc[i].backsubstitution(self._u[i])
318 else:
319 assert isinstance(self._u, _fem.Function)
320 self.mpc.homogenize(self._u)
321 self.mpc.backsubstitution(self._u)
323 converged_reason = self.solver.getConvergedReason()
324 return self._u, converged_reason, self.solver.getIterationNumber()
327class LinearProblem(dolfinx.fem.petsc.LinearProblem):
328 """
329 Class for solving a linear variational problem with multi point constraints of the form
330 a(u, v) = L(v) for all v using PETSc as a linear algebra backend.
332 Args:
333 a: A bilinear UFL form, the left hand side of the variational problem.
334 L: A linear UFL form, the right hand side of the variational problem.
335 mpc: The multi point constraint.
336 bcs: A list of Dirichlet boundary conditions.
337 u: The solution function. It will be created if not provided. The function has
338 to be based on the functionspace in the mpc, i.e.
340 .. highlight:: python
341 .. code-block:: python
343 u = dolfinx.fem.Function(mpc.function_space)
344 petsc_options: Parameters that is passed to the linear algebra backend PETSc. #type: ignore
345 For available choices for the 'petsc_options' kwarg, see the PETSc-documentation
346 https://www.mcs.anl.gov/petsc/documentation/index.html.
347 form_compiler_options: Parameters used in FFCx compilation of this form. Run `ffcx --help` at
348 the commandline to see all available options. Takes priority over all
349 other parameter values, except for `scalar_type` which is determined by DOLFINx.
350 jit_options: Parameters used in CFFI JIT compilation of C code generated by FFCx.
351 See https://github.com/FEniCS/dolfinx/blob/main/python/dolfinx/jit.py#L22-L37
352 for all available parameters. Takes priority over all other parameter values.
353 P: A preconditioner UFL form.
354 Examples:
355 Example usage:
357 .. highlight:: python
358 .. code-block:: python
360 problem = LinearProblem(a, L, mpc, [bc0, bc1],
361 petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
363 """
365 u: typing.Union[_fem.Function, list[_fem.Function]]
366 _a: typing.Union[_fem.Form, typing.Sequence[typing.Sequence[_fem.Form]]]
367 _L: typing.Union[_fem.Form, typing.Sequence[_fem.Form]]
368 _preconditioner: typing.Optional[typing.Union[_fem.Form, typing.Sequence[typing.Sequence[_fem.Form]]]]
369 _mpc: typing.Union[MultiPointConstraint, typing.Sequence[MultiPointConstraint]]
370 _A: PETSc.Mat # type: ignore
371 _P: typing.Optional[PETSc.Mat] # type: ignore
372 _b: PETSc.Vec # type: ignore
373 _solver: PETSc.KSP # type: ignore
374 _x: PETSc.Vec # type: ignore
375 bcs: typing.List[_fem.DirichletBC]
376 __slots__ = tuple(__annotations__)
378 def __init__(
379 self,
380 a: typing.Union[ufl.Form, typing.Iterable[typing.Iterable[ufl.Form]]],
381 L: typing.Union[ufl.Form, typing.Iterable[ufl.Form]],
382 mpc: typing.Union[MultiPointConstraint, typing.Sequence[MultiPointConstraint]],
383 bcs: typing.Optional[typing.List[_fem.DirichletBC]] = None,
384 u: typing.Optional[typing.Union[_fem.Function, typing.Sequence[_fem.Function]]] = None,
385 petsc_options_prefix: str = "dolfinx_mpc_linear_problem_",
386 petsc_options: typing.Optional[dict] = None,
387 form_compiler_options: typing.Optional[dict] = None,
388 jit_options: typing.Optional[dict] = None,
389 P: typing.Optional[typing.Union[ufl.Form, typing.Iterable[ufl.Form]]] = None,
390 ):
391 # Compile forms
392 form_compiler_options = {} if form_compiler_options is None else form_compiler_options
393 jit_options = {} if jit_options is None else jit_options
394 self._a = _fem.form(a, jit_options=jit_options, form_compiler_options=form_compiler_options)
395 self._L = _fem.form(L, jit_options=jit_options, form_compiler_options=form_compiler_options)
397 self._mpc = mpc
398 # Nest assembly
399 if isinstance(mpc, Sequence):
400 is_nest = True
401 # Sanity check
402 for mpc_i in mpc:
403 if not mpc_i.finalized:
404 raise RuntimeError("The multi point constraint has to be finalized before calling initializer")
405 # Create function containing solution vector
406 else:
407 is_nest = False
408 if not mpc.finalized:
409 raise RuntimeError("The multi point constraint has to be finalized before calling initializer")
411 # Create function(s) containing solution vector(s)
412 if is_nest:
413 if u is None:
414 assert isinstance(self._mpc, Sequence)
415 self.u = [_fem.Function(self._mpc[i].function_space) for i in range(len(self._mpc))]
416 else:
417 assert isinstance(self._mpc, Sequence)
418 assert isinstance(u, Sequence)
419 for i, (mpc_i, u_i) in enumerate(zip(self._mpc, u)):
420 assert isinstance(u_i, _fem.Function)
421 assert isinstance(mpc_i, MultiPointConstraint)
422 if u_i.function_space is not mpc_i.function_space:
423 raise ValueError(
424 "The input function has to be in the function space in the multi-point constraint",
425 "i.e. u = dolfinx.fem.Function(mpc.function_space)",
426 )
427 self.u = list(u)
428 else:
429 if u is None:
430 assert isinstance(self._mpc, MultiPointConstraint)
431 self.u = _fem.Function(self._mpc.function_space)
432 else:
433 assert isinstance(u, _fem.Function)
434 assert isinstance(self._mpc, MultiPointConstraint)
435 if u.function_space is self._mpc.function_space:
436 self.u = u
437 else:
438 raise ValueError(
439 "The input function has to be in the function space in the multi-point constraint",
440 "i.e. u = dolfinx.fem.Function(mpc.function_space)",
441 )
443 # Create MPC matrix and vector
444 self._preconditioner = _fem.form(P, jit_options=jit_options, form_compiler_options=form_compiler_options)
446 if is_nest:
447 assert isinstance(mpc, Sequence)
448 assert isinstance(self._L, Sequence)
449 assert isinstance(self._a, Sequence)
450 self._A = create_matrix_nest(self._a, mpc)
451 self._b = create_vector_nest(self._L, mpc)
452 self._x = create_vector_nest(self._L, mpc)
453 if self._preconditioner is None:
454 self._P_mat = None
455 else:
456 assert isinstance(self._preconditioner, Sequence)
457 self._P_mat = create_matrix_nest(self._preconditioner, mpc)
458 else:
459 assert isinstance(mpc, MultiPointConstraint)
460 assert isinstance(self._L, _fem.Form)
461 assert isinstance(self._a, _fem.Form)
462 self._A = _cpp_mpc.create_matrix(self._a._cpp_object, mpc._cpp_object)
463 self._b = create_vector(mpc.function_space.dofmap.index_map, mpc.function_space.dofmap.index_map_bs)
464 self._x = create_vector(mpc.function_space.dofmap.index_map, mpc.function_space.dofmap.index_map_bs)
465 if self._preconditioner is None:
466 self._P_mat = None
467 else:
468 assert isinstance(self._preconditioner, _fem.Form)
469 self._P_mat = _cpp_mpc.create_matrix(self._preconditioner._cpp_object, mpc._cpp_object)
471 self.bcs = [] if bcs is None else bcs
473 if is_nest:
474 assert isinstance(self.u, Sequence)
475 comm = self.u[0].function_space.mesh.comm
476 else:
477 assert isinstance(self.u, _fem.Function)
478 comm = self.u.function_space.mesh.comm
480 self._solver = PETSc.KSP().create(comm) # type: ignore
481 self._solver.setOperators(self._A, self._P_mat)
483 self._petsc_options_prefix = petsc_options_prefix
484 self.solver.setOptionsPrefix(petsc_options_prefix)
485 self.A.setOptionsPrefix(f"{petsc_options_prefix}A_")
486 self.b.setOptionsPrefix(f"{petsc_options_prefix}b_")
487 self.x.setOptionsPrefix(f"{petsc_options_prefix}x_")
488 if self.P_mat is not None:
489 self.P_mat.setOptionsPrefix(f"{petsc_options_prefix}P_mat_")
491 # Set options on KSP only
492 if petsc_options is not None:
493 opts = PETSc.Options() # type: ignore
494 opts.prefixPush(self.solver.getOptionsPrefix())
496 for k, v in petsc_options.items():
497 opts[k] = v
499 self.solver.setFromOptions()
501 # Tidy up global options
502 for k in petsc_options.keys():
503 del opts[k]
504 opts.prefixPop()
506 def solve(self) -> typing.Union[list[_fem.Function], _fem.Function]:
507 """Solve the problem.
509 Returns:
510 Function containing the solution"""
512 # Assemble lhs
513 self._A.zeroEntries()
514 if self._A.getType() == "nest":
515 assemble_matrix_nest(self._A, self._a, self._mpc, self.bcs, diagval=1.0) # type: ignore
516 else:
517 assert isinstance(self._a, _fem.Form)
518 assemble_matrix(self._a, self._mpc, bcs=self.bcs, A=self._A)
520 self._A.assemble()
521 assert self._A.assembled
523 # Assemble the preconditioner if provided
524 if self._P_mat is not None:
525 self._P_mat.zeroEntries()
526 if self._P_mat.getType() == "nest":
527 assert isinstance(self._preconditioner, typing.Sequence)
528 assemble_matrix_nest(self._P_mat, self._preconditioner, self._mpc, self.bcs) # type: ignore
529 else:
530 assert isinstance(self._preconditioner, _fem.Form)
531 assemble_matrix(self._preconditioner, self._mpc, bcs=self.bcs, A=self._P_mat)
532 self._P_mat.assemble()
534 # Assemble the residual
535 _zero_vector(self._b)
536 if self._x.getType() == "nest":
537 assemble_vector_nest(self._b, self._L, self._mpc) # type: ignore
538 else:
539 assert isinstance(self._L, _fem.Form)
540 assert isinstance(self._mpc, MultiPointConstraint)
541 assemble_vector(self._L, self._mpc, self._b)
543 # Lift vector
544 try:
545 # Nest and blocked lifting
546 bcs1 = _fem.bcs.bcs_by_block(_fem.forms.extract_function_spaces(self._a, 1), self.bcs) # type: ignore
547 apply_lifting(self._b, self._a, bcs=bcs1, constraint=self._mpc) # type: ignore
548 _ghost_update(self._b, PETSc.InsertMode.ADD, PETSc.ScatterMode.REVERSE) # type: ignore
549 bcs0 = _fem.bcs.bcs_by_block(_fem.forms.extract_function_spaces(self._L), self.bcs) # type: ignore
550 _fem.petsc.set_bc(self._b, bcs0)
551 except RuntimeError:
552 # Single form lifting
553 apply_lifting(self._b, [self._a], bcs=[self.bcs], constraint=self._mpc) # type: ignore
554 _ghost_update(self._b, PETSc.InsertMode.ADD, PETSc.ScatterMode.REVERSE) # type: ignore
555 _fem.petsc.set_bc(self._b, self.bcs)
556 _ghost_update(self._b, PETSc.InsertMode.INSERT, PETSc.ScatterMode.FORWARD) # type: ignore
558 # Solve linear system and update ghost values in the solution
559 self._solver.solve(self._b, self._x)
560 _ghost_update(self._x, PETSc.InsertMode.INSERT, PETSc.ScatterMode.FORWARD) # type: ignore
561 _fem.petsc.assign(self._x, self.u)
563 if isinstance(self.u, Sequence):
564 assert isinstance(self._mpc, Sequence)
565 for i in range(len(self.u)):
566 self._mpc[i].homogenize(self.u[i])
567 self._mpc[i].backsubstitution(self.u[i])
568 else:
569 assert isinstance(self.u, _fem.Function)
570 assert isinstance(self._mpc, MultiPointConstraint)
571 self._mpc.homogenize(self.u)
572 self._mpc.backsubstitution(self.u)
574 return self.u