Source code for dolfinx.fem.linearproblem
# -*- coding: utf-8 -*-
# Copyright (C) 2020 Jørgen S. Dokken
#
# This file is part of DOLFINX (https://www.fenicsproject.org)
#
# SPDX-License-Identifier: LGPL-3.0-or-later
import typing
import ufl
from dolfinx import fem
from petsc4py import PETSc
[docs]class LinearProblem():
"""Class for solving a linear variational problem of the form
a(u, v) = L(v) for all v using PETSc as a linear algebra backend.
"""
def __init__(self, a: ufl.Form, L: ufl.Form, bcs: typing.List[fem.DirichletBC] = [],
u: fem.Function = None, petsc_options={}, form_compiler_parameters={}, jit_parameters={}):
"""Initialize solver for a linear variational problem.
Parameters
----------
a
A bilinear UFL form, the left hand side of the variational problem.
L
A linear UFL form, the right hand side of the variational problem.
bcs
A list of Dirichlet boundary conditions.
u
The solution function. It will be created if not provided.
petsc_options
Parameters that is passed to the linear algebra backend PETSc.
For available choices for the 'petsc_options' kwarg, see the
`PETSc-documentation <https://www.mcs.anl.gov/petsc/documentation/index.html>`.
form_compiler_parameters
Parameters used in FFCX compilation of this form. Run `ffcx --help` at
the commandline to see all available options. Takes priority over all
other parameter values, except for `scalar_type` which is determined by
DOLFINX.
jit_parameters
Parameters used in CFFI JIT compilation of C code generated by FFCX.
See `python/dolfinx/jit.py` for all available parameters.
Takes priority over all other parameter values.
.. code-block:: python
problem = LinearProblem(a, L, [bc0, bc1], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
"""
self._a = fem.Form(a, form_compiler_parameters=form_compiler_parameters, jit_parameters=jit_parameters)
self._A = fem.create_matrix(self._a)
self._L = fem.Form(L, form_compiler_parameters=form_compiler_parameters, jit_parameters=jit_parameters)
self._b = fem.create_vector(self._L)
if u is None:
# Extract function space from TrialFunction (which is at the
# end of the argument list as it is numbered as 1, while the
# Test function is numbered as 0)
self.u = fem.Function(a.arguments()[-1].ufl_function_space())
else:
self.u = u
self.bcs = bcs
self._solver = PETSc.KSP().create(self.u.function_space.mesh.mpi_comm())
self._solver.setOperators(self._A)
# Give PETSc solver options a unique prefix
solver_prefix = "dolfinx_solve_{}".format(id(self))
self._solver.setOptionsPrefix(solver_prefix)
# Set PETSc options
opts = PETSc.Options()
opts.prefixPush(solver_prefix)
for k, v in petsc_options.items():
opts[k] = v
opts.prefixPop()
self._solver.setFromOptions()
[docs] def solve(self) -> fem.Function:
"""Solve the problem."""
# Assemble lhs
self._A.zeroEntries()
fem.assemble_matrix(self._A, self._a, bcs=self.bcs)
self._A.assemble()
# Assemble rhs
with self._b.localForm() as b_loc:
b_loc.set(0)
fem.assemble_vector(self._b, self._L)
# Apply boundary conditions to the rhs
fem.apply_lifting(self._b, [self._a], [self.bcs])
self._b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
fem.set_bc(self._b, self.bcs)
# Solve linear system and update ghost values in the solution
self._solver.solve(self._b, self.u.vector)
self.u.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
return self.u
@property
def L(self) -> fem.Form:
"""Get the compiled linear form"""
return self._L
@property
def a(self) -> fem.Form:
"""Get the compiled bilinear form"""
return self._a
@property
def A(self) -> PETSc.Mat:
"""Get the matrix operator"""
return self._A
@property
def b(self) -> PETSc.Vec:
"""Get the RHS vector"""
return self._b
@property
def solver(self) -> PETSc.KSP:
"""Get the linear solver"""
return self._solver