solving

solving#

This sub-module contains functions and types that are used to implement different solvers for the system of equations. They rely on the types implemented in C to speed the most time consuming parts of the calculations.

The key idea behind all these is to write it in a way that would make it straight forward to port to a distributed memory architecture. As such, data and operations are made in a way which requires minimal communication between the different element data.

Support Types#

A wrapper around mfv2d._mfv2d.LinearSystem is provided by LinearSystem. This provides an easier to deal with constructor, as well as a few convenience methods.

class mfv2d.solving.LinearSystem(n_elem: int, form_spec: ElementFormSpecification, orders: npt.NDArray[np.uint32], element_matrices: Sequence[npt.NDArray[np.float64]], constraints: Sequence[Constraint])[source]#

Class used to represent a linear system with element equations and constraints.

Parameters:

  • n_elem (int) – Number of elements in the system.

  • form_spec (ElementFormSpecification) – Specification of forms in the element system.

  • orders (array) – Array of element orders.

  • element_matrices (Sequence of arrays) – Sequence of element matrices.

  • constraints (Sequence of Constraint) – Sequences of constraints the element systems are subject to.

static __new__(cls, n_elem: int, form_spec: ElementFormSpecification, orders: npt.NDArray[np.uint32], element_matrices: Sequence[npt.NDArray[np.float64]], constraints: Sequence[Constraint]) → Self[source]#

Create a new LinearSystem.

apply_diagonal(x: DenseVector, out: DenseVector, /) → None#

Apply multiplication by the diagonal part of the system.

Parameters:

  • DenseVector – Dense vector to which this is applied.

  • DenseVector – Dense vector to which the output is returned.

apply_diagonal_inverse(x: DenseVector, out: DenseVector, /) → None#

Apply inverse of the diagonal part of the system.

Parameters:

  • DenseVector – Dense vector to which this is applied.

  • DenseVector – Dense vector to which the output is returned.

apply_full_trace_system(x: TraceVector, /, out: TraceVector, tmp1: DenseVector, tmp2: DenseVector) → None[source]#

Apply the Schur’s trace system.

Parameters:

  • TraceVector – Trace vector to which this is applied to.

  • out (TraceVector) – Trace vector, which receives output. Can be the same as input.

  • tmp1 (DenseVector) – Dense vector used to store intermediate result. Must be unique.

  • tmp2 (DenseVector) – Dense vector used to store intermediate result. Must be unique.

apply_trace(x: DenseVector, out: TraceVector, /) → None#

Apply the trace constraints to the dense vector.

Parameters:

  • DenseVector – Dense vector to which this is applied.

  • TraceVector – Trace vector to which the output is returned.

apply_trace_transpose(x: TraceVector, out: DenseVector, /) → None#

Apply the transpose of the constraints to the trace vector.

Parameters:

  • TraceVector – Trace vector to which this is applied.

  • DenseVector – Dense vector to which the output is returned.

combined_system_matrix() → csr_array[source]#

Combine the system matrix into a scipy CSR array.

get_constraint_blocks()#

get_system_constraint_blocks() -> tuple[MatrixCRS, …] Get the constraint blocks of the system.

get_dense_blocks()#

get_system_dense_blocks() -> tuple[numpy.typing.NDArray[numpy.float64], …] Get the dense blocks of the system.

To make writing some solvers easier FullVector is provided. This type just packs the dense parts represented by mfv2d._mfv2d.DenseVector and the trace parts mfv2d._mfv2d.TraceVector together and provides methods these two provide.

class mfv2d.solving.FullVector(dense: DenseVector, trace: TraceVector)[source]#

Class contaning a dense and trace vector.

__post_init__() → None[source]#

Check both are from the same parent.

static add(v1: FullVector, v2: FullVector, v_out: FullVector, k: float, /) → None[source]#

Add two vectors.

copy() → FullVector[source]#

Create a copy of itself.

dense: DenseVector#
static dot(v1: FullVector, v2: FullVector) → float[source]#

Compute dot product of two vectors.

classmethod make_empty(system: LinearSystem) → Self[source]#

Create an empty vector.

static scale(v: FullVector, k: float, v_out: FullVector, /) → None[source]#

Scale the vector by a constant.

set_from(other: FullVector) → None[source]#

Set the vector from another.

static subtract(v1: FullVector, v2: FullVector, v_out: FullVector, k: float, /) → None[source]#

Subtracts two vectors.

trace: TraceVector#

Core Solvers#

The solvers implemented in this sub-module all follow the phylosophy that they should be type-agnostic. As such, they all take the system along with all functions needed to operate on it. This allows the same solver to be used for different vector and system types by just writing a few short wrapper functions.

mfv2d.solving.gmres_general(mat: _Mat, rhs: _Vec, initial_guess: _Vec, convergence: ConvergenceSettings, system_application_function: Callable[[_Mat, _Vec, _Vec], None], vec_dot_function: Callable[[_Vec, _Vec], float], vec_add_to_function: Callable[[_Vec, _Vec, _Vec, float], None], vec_sub_from_scaled_function: Callable[[_Vec, _Vec, _Vec, float], None], vec_scale_by_function: Callable[[_Vec, float, _Vec], None], vec_copy_function: Callable[[_Vec], _Vec]) → tuple[_Vec, float, int][source]#

General implementation of GMRES to use for any data type with operators.

Returns:

  • _Vec – Computed solution.

  • float – Estimated residual.

  • int – Iterations done.

mfv2d.solving.cg_general(mat: _Mat, rhs: _Vec, initial_guess: _Vec, convergence: ConvergenceSettings, system_application_function: Callable[[_Mat, _Vec], None], vec_dot_function: Callable[[_Vec, _Vec], float], vec_add_to_scaled_function: Callable[[_Vec, _Vec, float], None], vec_sub_from_scaled_function: Callable[[_Vec, _Vec, float], None], vec_copy_function: Callable[[_Vec], _Vec], vec_set_function: Callable[[_Vec, _Vec], None]) → tuple[_Vec, float, int][source]#

General implementation of GMRES to use for any data type with operators.

Parameters:

  • m (int) – Number of basis to use for the solution.

  • mat (_Mat) – System matrix state.

  • rhs (_Vec) – Right side of the system.

  • convergence (ConvergenceSettings) – Settings to use for convergence.

  • system_application_function ((_Mat, _Vec) -> _Vec) – Function that applies the system to the input vector.

  • vec_dot_function ((_Vec, _Vec) -> float) – Function that computes the dot product of system vectors.

  • vec_add_function ((_Vec, _Vec) -> _Vec) – Function that adds two vectors together.

  • vec_sub_function ((_Vec, _Vec) -> _Vec) – Function which subtracts two vectors.

  • vec_scale_function ((float, _Vec) -> _Vec) – Function which scales a vector by a constant

Returns:

  • _Vec – Computed solution.

  • float – Estimated residual.

  • int – Iterations done.

mfv2d.solving.pcg_general(mat: _Mat, rhs: _Vec, initial_guess: _Vec, convergence: ConvergenceSettings, system_application_function: Callable[[_Mat, _Vec, _Vec], None], precondition_function: Callable[[_Mat, _Vec, _Vec], None], vec_dot_function: Callable[[_Vec, _Vec], float], vec_add_to_scaled_function: Callable[[_Vec, _Vec, float, _Vec], None], vec_sub_from_scaled_function: Callable[[_Vec, _Vec, float, _Vec], None], vec_copy_function: Callable[[_Vec], _Vec], degen_limit: float = 1e-12) → tuple[_Vec, float, int][source]#

General implementation of GMRES to use for any data type with operators.

Returns:

  • _Vec – Computed solution.

  • float – Estimated residual.

  • int – Iterations done.

System Solvers#

Solvers which actually solve the system are based on the core solvers and use them at some point, though the specific way the do it can differ.

GMRES#

GMRES may be slow and memory hungry, but it is guaranteed to give you the best solution possible from a Krylov method solver. It solves the system by using the FullVector and solve it all at once.

mfv2d.solving.solve_gmres_iterative(system: LinearSystem, rhs: DenseVector, constraints: TraceVector, convergence: ConvergenceSettings) → tuple[DenseVector, TraceVector, float, int][source]#

Solve the system using GMRES.

CG#

Conjugate gradient is the optimal method for symmetric positive definite (SPD) matrices. Unfortunately, mixed problems (mixed Poisson, Navier-Stokes) give rise to indefinete problems, which means that CG will likely terminate due to the degeneration of the Krylov subspace.

mfv2d.solving.solve_cg_iterative(system: LinearSystem, rhs: DenseVector, constraints: TraceVector, convergence: ConvergenceSettings) → tuple[DenseVector, TraceVector, float, int][source]#

Solve the system using CG.

PCG#

Preconditioned conjugate gradient can deal with the problems of CG solver by applying a preconditioner to the system. This allows for both faster convergence and nicer behavior with indefinete systems. It uses block-Jacobi preconditioner.

mfv2d.solving.solve_pcg_iterative(system: LinearSystem, rhs: DenseVector, constraints: TraceVector, convergence: ConvergenceSettings) → tuple[DenseVector, TraceVector, float, int][source]#

Solve the system using preconditioned CG.

Schur’s Complement#

Schur’s complement is the best performing solver I have written thus far. The idea is to use Gaussian eleimination to obtain the system for Lagrange multipliers using the equation (1), where \(\vec{\lambda}\) are the Lagrange multipliers, \(\mathbf{N}\) the matrix that enforces the constraints, \(\mathbf{A}\) the block diagonal part, and the \(\vec{y}\) and \(\vec{\phi}\) being the dense and trace forcing respectively.

While the matrix \(\mathbf{N} \mathbf{A}^{-1} \mathbf{N}^T\) is dense, the is solved iteratively using PCG in a matrix-free manner. This is because the dense solve would scale with the cube of the degrees of freedom.

(1)#\[\mathbf{N} \mathbf{A}^{-1} \mathbf{N}^T \vec{\lambda} = \mathbf{N} \mathbf{A}^{-1} \vec{y} - \phi\]
Contents