solve_system#
This submodule contains all the basic building blocks which are needed to
create the solver in mfv2d.solve_system_2d. Some of these should
be replaced or removed, as they are outdated, but it works for now.
Computing the Explicit Terms#
To compute mfv2d.kform.KExplicit terms of the
mfv2d.kform.KFormSystem the function compute_element_rhs()
is used, which is to be passed to mfv2d.element.call_per_leaf_flex()
to compute forcing terms. Under the hood it just calls _extract_rhs_2d()
for each element and then joins them together. This function in turn calls
rhs_2d_element_projection(), since boundary projections are evaluated
separately.
- mfv2d.solve_system.compute_element_rhs(ie: int, system: KFormSystem, element_caches: ObjectElementArray[ElementMassMatrixCache]) ndarray[tuple[Any, ...], dtype[float64]][source]#
Compute rhs for an element.
This basically means just concatenating the projections of the functions on the element for each of the equations in the system.
- Parameters:
ie (int) – Index of the element.
system (KFormSystem) – System for which to compute the rhs.
element_caches (FemCache) – Cache from which to get the basis from.
- Returns:
Array with the resulting rhs.
- Return type:
array
- mfv2d.solve_system._extract_rhs_2d(proj: Sequence[tuple[float, KExplicit]], weight: KWeight, element_cache: ElementMassMatrixCache) ndarray[tuple[Any, ...], dtype[float64]][source]#
Extract the rhs resulting from element projections.
Combines the sequence of
KExplicitterms together.- Parameters:
- Returns:
Array of the resulting projection degrees of freedom.
- Return type:
array
- mfv2d.solve_system.rhs_2d_element_projection(right: KElementProjection, element_cache: ElementMassMatrixCache) ndarray[tuple[Any, ...], dtype[float64]][source]#
Evaluate the differential form projections on the 1D element.
- Parameters:
right (KElementProjection) – The projection of a function on the element.
corners ((4, 2) array) – Array with corners of the element.
basis (Basis2D) – Basis to use for computing the projection.
- Returns:
The resulting projection vector.
- Return type:
array of
numpy.float64
Matrix Assembly#
Global system matrix is created by assemble_matrix(), which creates
individual merges leaf element matrices. This used to be more involved,
but now we are in better times.
- mfv2d.solve_system.assemble_matrix(leaf_matrices: FlexibleElementArray[float64, uint32]) csr_array[source]#
Assemble global matrix.
Vector Assembly#
Global system vector is similarly created by assemble_vector(). It used
to be slow, complicated, and painful, but it’s all ogre now.
Forcing Assembly#
Last of the assembly routines is the assemble_forcing() function. This
assembles together the leaf element forcing (value of the given expression
given solution).
Global Reconstruction#
Based on mfv2d.element.reconstruct() a global reconstruction
of the solution for all leaf elements is implemented with the
reconstruct_mesh_from_solution() function. It returns a
pyvista.UnstructuredGrid, which can be added to the list
of outputs.
- mfv2d.solve_system.reconstruct_mesh_from_solution(system: KFormSystem, recon_order: int | None, element_collection: ElementCollection, caches: FemCache, dof_offsets: FixedElementArray[uint32], solution: FlexibleElementArray[float64, uint32]) UnstructuredGrid[source]#
Reconstruct the unknown differential forms.
Time March Support#
Since for time marching certain quantities must be extracted from
the non-constraint equations, supporting functions are provided here. Namely,
_extract_time_carry() is called on each element, which
extract_carry() wraps. To determine what are the indices of degrees
of freedom to carry find_time_carry_indices() is used.
- mfv2d.solve_system.extract_carry(element_collection: ElementCollection, time_carry_index_array: FlexibleElementArray[uint32, uint32], initial_solution: FlexibleElementArray[float64, uint32])[source]#
Extract carry terms from solution.
- mfv2d.solve_system._extract_time_carry(ie: int, time_carry_index_array: FlexibleElementArray[uint32, uint32], solution: FlexibleElementArray[float64, uint32]) ndarray[tuple[Any, ...], dtype[float64]][source]#
Element function for extraction.
Leaf Calculations#
For computing leaf element matrices and forcing, compute_leaf_matrix()
and compute_leaf_vector() are available. These pretty much just
coordinate calls to mfv2d._mfv2d.compute_element_matrix() and
mfv2d._mfv2d.compute_element_vector().
- mfv2d.solve_system.compute_leaf_matrix(ie: int, expressions: Sequence[Sequence[Sequence[MatOpCode | int | float] | None]], unknowns: UnknownOrderings, element_caches: ObjectElementArray[ElementMassMatrixCache], fields: FixedElementArray[object_]) ndarray[tuple[Any, ...], dtype[float64]][source]#
Compute the element matrix.
- mfv2d.solve_system.compute_leaf_vector(ie: int, expressions: Sequence[Sequence[Sequence[MatOpCode | int | float] | None]], unknowns: UnknownOrderings, element_caches: ObjectElementArray[ElementMassMatrixCache], fields: FixedElementArray[object_], solution: FlexibleElementArray[float64, uint32]) ndarray[tuple[Any, ...], dtype[float64]][source]#
Compute the element vector.
Since these also need vector field information for any interior products,
functions to compute are also given:
compute_element_vector_fields_nonlin() and its wrapper
compute_element_vector_fields().
- mfv2d.solve_system.compute_element_vector_fields_nonlin(unknown_forms: Sequence[KFormUnknown], element_basis: Basis2D, output_basis: Basis2D, vector_fields: Sequence[Function2D | KFormUnknown], element_corners: ndarray[tuple[Any, ...], dtype[float64]], unknown_offsets: ndarray[tuple[Any, ...], dtype[uint32]], solution: ndarray[tuple[Any, ...], dtype[float64]] | None) tuple[ndarray[tuple[Any, ...], dtype[float64]], ...][source]#
Evaluate vector fields which may be non-linear.
- Parameters:
unknown_forms (Sequence of KFormUnknown) – Unknown forms in the order they appear in the system. This is used to determine what degrees of freedom are needed for the vector fields based on differential forms.
element_basis (Basis2D) – Basis functions that the element uses. This needs to match with the number of degrees of freedom of the element.
output_basis (Basis2D) – Basis onto which the result is to be computed.
vector_fields (Sequence of Function2D or KFormUnknown) – Description of the vector fields. Can be a callable which gives its value at a point, or instead it can be an unknown in the system.
element_corners ((4, 2) array) – Array of the element corner points.
unknown_offsets (array) – Array with offsets of the degrees of freedom within the element. This is used to pick correct degrees of freedom from the element DoFs vector.
solution (array, optional) – Array of the element degrees of freedom. If not provided, all are assumed to be zero instead.
- Returns:
Tuple with arrays with values of the vector field at each point of the 2D basis integration rules.
- Return type:
tuple of array
- mfv2d.solve_system.compute_element_vector_fields(ie: int, system: KFormSystem, child_count_array: FixedElementArray[uint32], orders_in: FixedElementArray[uint32], orders_out: FixedElementArray[uint32], basis_cache: FemCache, vector_fields: Sequence[Function2D | KFormUnknown], corners: FixedElementArray[float64], dof_offsets: FixedElementArray[uint32], solution: FlexibleElementArray[float64, uint32]) tuple[source]#
Wrap the vector fields func.
The Actual Solver#
The actual solver that runs to solve the (potentially non-linear) system
is implemented in the non_linear_solve_run() function.
- mfv2d.solve_system.non_linear_solve_run(system: KFormSystem, max_iterations: int, relax: float, atol: float, rtol: float, print_residual: bool, unknown_ordering: UnknownOrderings, element_collection: ElementCollection, leaf_elements: ndarray[tuple[Any, ...], dtype[integer]], cache_2d: FemCache, element_caches: ObjectElementArray[ElementMassMatrixCache], compiled_system: CompiledSystem, explicit_vec: ndarray[tuple[Any, ...], dtype[float64]], dof_offsets: FixedElementArray[uint32], element_offsets: ndarray[tuple[Any, ...], dtype[uint32]], linear_element_matrices: FlexibleElementArray[float64, uint32], time_carry_index_array: FlexibleElementArray[uint32, uint32] | None, time_carry_term: FlexibleElementArray[float64, uint32] | None, solution: FlexibleElementArray[float64, uint32], global_lagrange: ndarray[tuple[Any, ...], dtype[float64]], max_mag: float, vector_fields: Sequence[Function2D | KFormUnknown], system_decomp: SuperLU, lagrange_mat: csr_array | None, return_all_residuals: bool = False) tuple[FlexibleElementArray[float64, uint32], ndarray[tuple[Any, ...], dtype[float64]], int, ndarray[tuple[Any, ...], dtype[float64]]][source]#
Run the iterative non-linear solver.
Based on how the compiled system looks, this may only take a single iteration, otherwise, it may run for as long as it needs to converge.
Solver Settings#
To configure the solver, several settings dataclasses are provided. These each handle a different aspect of the solver and most have some default values.
- class mfv2d.solve_system.SystemSettings(system: ~mfv2d.kform.KFormSystem, boundary_conditions: ~collections.abc.Sequence[~mfv2d.boundary.BoundaryCondition2DSteady] = <factory>, constrained_forms: ~collections.abc.Sequence[tuple[float, ~mfv2d.kform.KFormUnknown]] = <factory>, initial_conditions: ~collections.abc.Mapping[~mfv2d.kform.KFormUnknown, ~collections.abc.Callable[[~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.float64]], ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.float64]]], ~collections.abc.Buffer | ~numpy._typing._array_like._SupportsArray[~numpy.dtype[~typing.Any]] | ~numpy._typing._nested_sequence._NestedSequence[~numpy._typing._array_like._SupportsArray[~numpy.dtype[~typing.Any]]] | complex | bytes | str | ~numpy._typing._nested_sequence._NestedSequence[complex | bytes | str]]] = <factory>)[source]#
Type used to hold system information for solving.
- Parameters:
system (KFormSystem) – System of equations to solve.
boundaray_conditions (Sequence of BoundaryCondition2DSteady, optional) – Sequence of boundary conditions to be applied to the system.
constrained_forms (Sequence of (float, KFormUnknown), optional) –
Sequence of 2-form unknowns which must be constrained. These can arrise form cases where a continuous variable acts as a Lagrange multiplier on the continuous level and only appears in the PDE as a gradient. In that case it will result in a singular system if not constrained manually to a fixed value.
An example of such a case is pressure in Stokes flow or incompressible Navier-Stokes equations.
intial_conditions (Mapping of (KFormUnknown, Callable), optional) – Functions which give initial conditions for different forms.
- boundary_conditions: Sequence[BoundaryCondition2DSteady]#
- constrained_forms: Sequence[tuple[float, KFormUnknown]]#
- initial_conditions: Mapping[KFormUnknown, Callable[[ndarray[tuple[Any, ...], dtype[float64]], ndarray[tuple[Any, ...], dtype[float64]]], Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str]]]#
- system: KFormSystem#
- class mfv2d.solve_system.SolverSettings(maximum_iterations: int = 100, relaxation: float = 1.0, absolute_tolerance: float = 1e-06, relative_tolerance: float = 1e-05)[source]#
Settings used by the non-linear solver.
Solver finds a solution to the system (1), where \(\mathbb{I}\) is the implicit part of the system, \(\vec{E}\) is the explicit part of the system, and \(\vec{F}\) is the constant part of the system.
This is done by computing updates to the state \(\Delta u\) as per (2). The iterations stop once the residual \(\vec{E}\left({\vec{u}}^i\right) + \vec{F} + \vec{I}\left({\vec{u}}^i\right)\) falls under the specified criterion.
The stopping criterion has two tolerances - absolute and relative - which both need to be met in order for the solver to stop. The absolute criterion checks if the highest absolute value of the residual elements is bellow the value specified. The relative first scales it by the maximum absolute value of the constant forcing. Depending on the system and the solution, these may need to be adjusted. If the system is linear, meaning that \(\vec{E} = 0\) and \(\mathbb{I}\) is not dependant on math:vec{u}, then the solver will terminate in a single iteration.
Last thing to consider is the relaxation. Sometimes, the system is very stiff and does not converge nicely. This can be the case for steady-state calculations of non-linear systems with bad initial guesses. In such cases, the correction can be too large and overshoot the solution. In such cases, convergence may still be achieved if the update is scaled by a “relaxation factor”. Conversely, convergence may be slightly sped up for some very stable problems if the increment is amplified, meaning the relaxation factor is greater than 1.
(1)#\[\mathbb{I}\left({\vec{u}}\right) {\vec{u}} = \vec{E}\left({\vec{u}}^i \right) + \vec{F}\](2)#\[\Delta {\vec{u}}^i = \left(\mathbb{I}\left({\vec{u}}^i\right)\right)^{-1} \left( \vec{E}\left({\vec{u}}^i\right) + \vec{F} + \vec{I}\left({\vec{u}}^i\right)\right)\]- Parameters:
maximum_iterations (int, default: 100) – Maximum number of iterations the solver performs.
relaxation (float, default: 1.0) – Fraction of solution increment to use.
absolute_tolerance (float, default: 1e-6) – Maximum value of the residual must meet in order for the solution to be considered converged.
relative_tolerance (float, default: 1e-5) – Maximum fraction of the maximum of the right side of the equation the residual must meet in order for the solution to be considered converged.
- class mfv2d.solve_system.RefinementSettings(refinement_levels: int, division_predicate: ~collections.abc.Callable[[~mfv2d.mimetic2d.ElementLeaf2D, int], bool] | None = None, division_function: ~collections.abc.Callable[[int, int, int], tuple[int, int, int, int]] = <function divide_old>)[source]#
Type used to hold settings related to refinement information.
- Parameters:
refinement_levels (int) – Number of mesh refinement levels which can be done. When zero, no refinement is done.
division_predicate (Callable (Element2D, int) -> bool, optional) – Callable used to determine if an element should be divided further. If not specified, no refinement is done.
division_function (OrderDivisionFunction, optional) – Function which determines order of the parent and child elements resulting from the division of the element. When not specified, the “old” method is used.
- class mfv2d.solve_system.TimeSettings(dt: float, nt: int, time_march_relations: Mapping[KWeight, KFormUnknown], sample_rate: int = 1)[source]#
Type for defining time settings of the solver.
- Parameters:
dt (float) – Time step to take.
nt (int) – Number of time steps to simulate.
time_march_relations (dict of (KWeight, KFormUnknown)) – Pairs of weights and unknowns, which determine what equations are treated as time marching equations for which unknowns. At least one should be present.
sample_rate (int, optional) – How often the output is saved. If not specified, every time step is saved. First and last steps are always saved.
- time_march_relations: Mapping[KWeight, KFormUnknown]#
Solver Statistics#
To return feedback on how the solver ran and more information about the system
it just solved, SolutionStatistics type is provided. As such, an
object of this type is returned when the solver runs.
- class mfv2d.solve_system.SolutionStatistics(element_orders: dict[int, int], n_total_dofs: int, n_leaf_dofs: int, n_lagrange: int, n_elems: int, n_leaves: int, iter_history: ndarray[tuple[Any, ...], dtype[uint32]], residual_history: ndarray[tuple[Any, ...], dtype[float64]])[source]#
Information about the solution.