porepy.numerics.ad.equation_system module

Contains the EquationSystem, managing variables and equations for a system modelled using the AD framework.

class EquationSystem(mdg)[source]

Bases: object

Represents an equation system, modelled by AD variables and equations in AD form.

This class provides functionalities to create and manage variables, as well as managing equations on the form of Operator. It further provides functions to assemble subsystems and using subsets of equations and variables.

Note

As of now, the system matrix (Jacobian) is assembled with respect to ALL variables and then the columns belonging to the requested subset of variables and grids are sliced out and returned. This will be optimized with minor changes to the AD operator class and its recursive forward AD mode in the future.

Currently, this class optimizes the block structure of the Jacobian only regarding the subdomains and interfaces. A more localized optimization (e.g. cell-wise for equations without spatial differential operators) is not performed.

Parameters

mdg (MixedDimensionalGrid) –

SubSystem(equation_names=None, variable_names=None)[source]

Creates an EquationSystem for a given subset of equations and variables.

Currently only subsystems containing whole equations and variables in the mixed-dimensional sense can be created. Restrictions of equations to subdomains is not supported.

Parameters
Returns

A new instance of EquationSystem. The subsystem equations and variables are ordered as imposed by this systems’s order.

Raises

ValueError – if passed names are not among created variables and set equations.

Return type

EquationSystem

assemble(state=None)[source]

Assemble Jacobian matrix and residual vector of the whole system.

This is a shallow wrapper of assemble_subsystem(). Here, the subsystem is the complete set of equations, variables and grids.

Parameters

state (optional) – see assemble_subsystem(). Defaults to None.

Returns

Tuple containing

sps.spmatrix: Jacobian matrix corresponding to the targeted state. The ordering of the equations (rows) is determined by the order the equations were added. The DOFs (columns) are ordered according to the global order.

np.ndarray: Residual vector corresponding to the targeted state, scaled by -1 (moved to rhs).

Return type

tuple[scipy.sparse._base.spmatrix, numpy.ndarray]

assemble_schur_complement_system(primary_equations, primary_variables, inverter=None, state=None)[source]

Assemble Jacobian matrix and residual vector using a Schur complement elimination of the variables and equations not to be included.

The specified equations and variables will define blocks of the linearized system as

\[\begin{split}\left [ \begin{matrix} A_{pp} & A_{ps} \\ A_{sp} & A_{ss} \end{matrix} \right] \left [ \begin{matrix} x_p \\ x_s \end{matrix}\right] = \left [ \begin{matrix} b_p \\ b_s \end{matrix}\right]\end{split}\]

where subscripts p and s define primary and secondary blocks. The Schur complement system is then given by

\[\left( A_{pp} - A_{ps} * A_{ss}^{-1} * A_{sp}\right) * x_p = b_p - A_{ps} * A_{ss} * b_s\]

The Schur complement is well-defined only if the inverse of \(A_{ss}\) exists, and the efficiency of the approach assumes that an efficient inverter for \(A_{ss}\) can be found. The user must ensure both requirements are fulfilled.

Note

The optional arguments defining the secondary block, and the flag excl_loc_prim_to_sec are meant for nested Schur-complements and splitting solvers. This is an advanced usage and requires the user to be careful, since the resulting blocks \(A_{pp}\) and \(A_{ss}\) might end up to be not square. This will result in errors.

Examples

The default inverter can be defined by

import scipy.sparse as sps
inverter = lambda A: sps.csr_matrix(sps.linalg.inv(A.A))

It is costly in terms of computational time and memory, though.

TODO: We should rather use the block inverter in pp.matrix_operations. This will require some work on ensuring the system is block-diagonal.

Parameters

  • primary_equations (EquationList | EquationRestriction) – a subset of equations specifying the primary subspace in row-sense.

  • primary_variables (VariableList) – VariableType input specifying the primary subspace in column-sense.

  • inverter (optional) – callable object to compute the inverse of the matrix \(A_{ss}\). By default, the scipy direct sparse inverter is used.

  • state (optional) – see assemble_subsystem(). Defaults to None.

Returns

Tuple containing

sps.spmatrix: Jacobian matrix representing the Schur complement with respect to the targeted state. np.ndarray: Residual vector for the Schur complement with respect to the targeted state. Scaled with -1 (moved to rhs).

Raises

  • AssertionError –

    • If the primary block would have 0 rows or columns. - If the secondary block would have 0 rows or columns. - If the secondary block is not square.

  • ValueError – If primary and secondary columns overlap.

Return type

tuple[sps.spmatrix, np.ndarray]

assemble_subsystem(equations=None, variables=None, state=None)[source]

Assemble Jacobian matrix and residual vector using a specified subset of equations, variables and grids.

The method is intended for use in splitting algorithms. Matrix blocks not included will simply be sliced out.

Note

The ordering of columns in the returned system are defined by the global DOF order. The row blocks are in the same order as equations were added to this system. If an equation is defined on multiple grids, the respective row-block is internally ordered as given by the mixed-dimensional grid (for sd in subdomains, for intf in interfaces).

The columns of the subsystem are assumed to be properly defined by variables, otherwise a matrix of shape (M,) is returned. This happens if grid variables are passed which are unknown to this EquationSystem.

Parameters

  • equations (optional) –

    a subset of equations to which the subsystem should be restricted. If not provided (None), all equations known to this EquationSystem will be included.

    The user can specify grids per equation (name) to which the subsystem should be restricted in the row-sense. Grids not belonging to the domain of an equation will raise an error.

  • variables (optional) – VariableType input specifying the subspace in column-sense. If not provided (None), all variables will be included.

  • state (optional) – State vector to assemble from. By default, the stored pp.ITERATE or pp.STATE are used, in that order.

Returns

Tuple with two elements

spmatrix: (Part of the) Jacobian matrix corresponding to the targeted variable state, for the specified equations and variables. ndarray: Residual vector corresponding to the targeted variable state, for the specified equations. Scaled with -1 (moved to rhs).

Return type

tuple[sps.spmatrix, np.ndarray]

create_variables(name, dof_info=None, subdomains=None, interfaces=None, tags=None)[source]

Creates new variables according to specifications.

This method does not assign any values to the variable. This has to be done in a subsequent step (e.g. using set_var_values()).

Examples

An example on how to define a pressure variable with cell-wise one DOF (default) on all subdomains and no interfaces would be

p = ad_system.create_variables('pressure', subdomains=mdg.subdomains())
Parameters

  • name (str) – Name of the variable.

  • dof_info (Optional[dict[Literal['cells', 'faces', 'nodes'], int]]) – Dictionary containing information about number of DOFs per admissible type. Defaults to {'cells':1}.

  • subdomains (optional) – List of subdomains on which the variable is defined. If None, then it will not be defined on any subdomain.

  • interfaces (optional) – list of interfaces on which the variable is defined. If None, then it will not be defined on any interface.

  • tags (optional) – dictionary containing tags for the variables. The tags are assigned to all variables created by this method and can be updated using update_variable_tags().

Returns

A mixed-dimensional variable with above specifications.

Raises

  • ValueError – If non-admissible DOF types are used as local DOFs.

  • ValueError – If one attempts to create a variable not defined on any grid.

  • KeyError – If a variable with given name is already defined.

Return type

MixedDimensionalVariable

discretize(equations=None)[source]

Find and loop over all discretizations in the equation operators, extract unique references and discretize.

This is more efficient than discretizing on the Operator level, since discretizations which occur more than once in a set of equations will be identified and only discretized once.

Parameters

equations (optional) – A subset of equations. If not provided (None), all known equations will be discretized.

Return type

None

dofs_of(variables)[source]

Get the indices in the global vector of unknowns belonging to the variable(s).

Parameters

variables (Union[list[str], list[porepy.numerics.ad.operators.Variable], list[porepy.numerics.ad.operators.MixedDimensionalVariable]]) – VariableType input for which the indices are requested.

Returns

an order-preserving array of indices of DOFs belonging to the VariableType input.

Raises

ValueError – if unknown VariableType arguments are passed.

Return type

ndarray

expand_schur_complement_solution(reduced_solution)[source]

Expands the solution of the last assembled Schur complement system to the whole solution.

With reduced_solution as \(x_p\) from

\[\begin{split}\left [ \begin{matrix} A_{pp} & A_{ps} \\ A_{sp} & A_{ss} \end{matrix} \right] \left [ \begin{matrix} x_p \\ x_s \end{matrix}\right] = \left [ \begin{matrix} b_p \\ b_s \end{matrix}\right],\end{split}\]

the method returns the whole vector \([x_p, x_s]\), where

\[x_s = A_{ss}^{-1} * (b_s - A_{sp} * x_p).\]

Note

Independent of how the primary and secondary blocks were chosen, this method always returns a vector of size num_dofs. Especially when the primary and secondary variables did not constitute the whole vector of unknowns, the result is still of size num_dofs. The entries corresponding to the excluded grid variables are zero.

Parameters

reduced_solution (ndarray) – Solution to the linear system returned by assemble_schur_complement_system().

Returns

The expanded Schur solution in global size.

Raises

ValueError – If the Schur complement system was not assembled before.

Return type

ndarray

get_variable_values(variables=None, from_iterate=False)[source]

Assembles an array containing values for the passed variable-like argument.

The global order is preserved and independent of the order of the argument.

Parameters

  • variables (optional) – VariableType input for which the values are requested. If None (default), the global vector of unknowns is returned.

  • from_iterate (optional) – flag to return values stored as ITERATE, instead of STATE (default).

Returns

The respective (sub) vector in numerical format, size anywhere between 0 and

num_dofs().

Raises

  • KeyError – If no values are stored for the VariableType input.

  • ValueError – If unknown VariableType arguments are passed.

Return type

ndarray

get_variables(variables=None, grids=None, tag_name=None, tag_value=None)[source]

Filter variables based on grid, tag name and tag value.

Particular usage: calling without arguments will return all variables in the system.

Parameters
Returns

List of filtered variables.

Return type

list[porepy.numerics.ad.operators.Variable]

identify_dof(dof)[source]

Identifies the variable to which a specific DOF index belongs.

The intended use is to help identify entries in the global vector or the column of the Jacobian.

Parameters

dof (int) – a single index in the global vector.

Return type

Variable

Returns: the identified Variable object.

Raises

KeyError – if the dof is out of range (larger than num_dofs or smaller than 0).

Parameters

dof (int) –

Return type

Variable

md_variable(name, grids=None)[source]

Create a mixed-dimensional variable for a given name-domain list combination.

Parameters

  • name (str) – Name of the mixed-dimensional variable.

  • grids (optional) – List of grids where the variable is defined. If None (default), all grids where the variable is defined are used.

Returns

A mixed-dimensional variable.

Raises

ValueError – If variables name exist on both grids and interfaces and domain type is not specified (grids is None).

Return type

MixedDimensionalVariable

num_dofs()[source]

Returns the total number of dofs managed by this system.

Return type

int

projection_to(variables=None)[source]

Create a projection matrix from the global vector of unknowns to a specified subspace.

The transpose of the returned matrix can be used to slice respective columns out of the global Jacobian.

The projection preserves the global order defined by the system, i.e. it includes no permutation.

Parameters

variables (optional) – VariableType input for which the subspace is requested. If no subspace is specified using variables, a null-space projection is returned.

Returns

a sparse projection matrix of shape (M, num_dofs), where 0 <= M <= num_dofs.

Return type

csr_matrix

remove_equation(name)[source]

Removes a previously set equation and all related information.

Returns

A reference to the equation in operator form or None, if the equation is unknown.

Raises

ValueError – If an unknown equation is attempted removed.

Parameters

name (str) –

Return type

Operator | None

set_equation(equation, grids, equations_per_grid_entity)[source]

Sets an equation using the passed operator and uses its name as an identifier.

If an equation already exists under that name, it is overwritten.

Information about the image space must be provided for now, such that grid-wise row slicing is possible. This will hopefully be provided automatically in the future.

Note

Regarding the number of equations, this method assumes that the AD framework assembles row blocks per grid in subdomains, then per grid in interfaces, for each operator representing an equation. This is assumed to be the way PorePy AD works.

Parameters

  • equation (Operator) – An equation in AD operator form, assuming the right-hand side is zero and this instance represents the left-hand side.

  • grids (Union[list[porepy.grids.grid.Grid], list[porepy.grids.mortar_grid.MortarGrid]]) – A list of subdomain or interface grids on which the equation is defined.

  • equations_per_grid_entity (dict[Literal['cells', 'faces', 'nodes'], int]) – a dictionary describing how many equations equation_operator provides. This is a temporary work-around until operators are able to provide information on their image space. The dictionary must contain the number of equations per grid entity (cells, faces, nodes) for the operator.

Raises

  • ValueError – If the equation operator has a name already assigned to a previously set equation.

  • ValueError – If the equation is defined on both subdomains and interfaces.

  • AssertionError – If the equation is defined on an unknown grid.

  • ValueError – If indicated number of equations does not match the actual number as per evaluation of operator.

Return type

None

set_variable_values(values, variables=None, to_state=False, to_iterate=False, additive=False)[source]

Sets values for a (sub) vector of the global vector of unknowns.

The order of values is assumed to fit the global order.

Note

The vector is assumed to be of proper size and will be dissected according to the global order, starting with the index 0. Mismatches of is-size and should-be-size according to the subspace specified by variables will raise respective errors by numpy.

Parameters

  • values (ndarray) – Vector of size corresponding to number of DOFs of the specified variables.

  • variables (optional) – VariableType input for which the values are requested. If None (default), the global vector of unknowns will be set.

  • to_state (optional) – Flag to write values to pp.STATE.

  • to_iterate (optional) – Flag to write values to pp.ITERATE.

  • additive (optional) – Flag to write values additively. To be used in iterative procedures.

Raises

ValueError – If unknown VariableType arguments are passed.

Return type

None

update_variable_tags(tags, variables=None)[source]

Assigns tags to variables.

Parameters
Return type

None

admissible_dof_types: tuple[Literal['cells'], Literal['faces'], Literal['nodes']]

A set denoting admissible types of local DOFs for variables.

  • nodes: DOFs per grid node.

  • cells: DOFs per grid cell.

  • faces: DOFS per grid face.

assembled_equation_indices: dict[str, numpy.ndarray]

Contains the row indices in the last assembled (sub-) system for a given equation name (key).

This dictionary changes with every call to any assemble-method.

property equations: dict[str, porepy.numerics.ad.operators.Operator]

Dictionary containing names of operators (keys) and operators (values), which have been set as equations in this system.

mdg: MixedDimensionalGrid

Mixed-dimensional domain passed at instantiation.

property variable_domains: list[Union[porepy.grids.grid.Grid, porepy.grids.mortar_grid.MortarGrid]]

List containing all domains where at least one variable is defined.

property variables: list[porepy.numerics.ad.operators.Variable]

List containing all :class:`~porepy.numerics.ad.Variable`s known to this system.