porepy.numerics.ad package

Init file for all AD functionality.

They should all be accessible through a calling
>>> import porepy as pp
>>> pp.ad.Matrix???

etc.

class ADmethod(func=None, ad_function_type=Function, operator_kwargs={})[source]

Bases: object

(Decorator) Class for methods representing e.g., physical properties. The decorated function is expected to take scalars/vectors and return a scalar/vector.

The return value will be an AD operator of a type passed to the decorator.

Examples

import porepy as pp

# decorating class methods
class IdealGas:

    @ADmethod(ad_operator=pp.ad.DiagonalJacobianFunction,
            operators_args={"multipliers"=[1,1]})
    def density(self, p: float, T: float) -> float:
        return p/T

# decorating function
@ADmethod(ad_operator=pp.ad.Function)
def dummy_rel_perm(s):
    return s**2

With above code, the density of an instance of IdealGas can be called using MergedVariable representing pressure and temperature. Analogously, dummy_rel_perm can be called with one representing the saturation.

Note

If used as decorator WITHOUT explicit instantiation, the instantiation will be done implicitly with above default arguments (that’s how Python decorators work).

Parameters

  • func (Optional[Callable]) – decorated function object

  • ad_function_type (Type[AbstractFunction]) – type reference to an AD operator function to be instantiated. When instantiated, that instance will effectively replace func.

  • operator_kwargs (dict) – keyword arguments to be passed when instantiating an operator of type ad_function_type.

ad_wrapper(*args, **kwargs)[source]

Actual wrapper function. Constructs the necessary AD-Operator class wrapping the decorated callable and performs the evaluation/call.

Parameters

  • *args – arguments for the call to the wrapping AD operator function

  • **kwargs – keyword argument for the call to the wrapping Ad operator function

Return type

Operator

class Ad_array(val=1.0, jac=0.0)[source]

Bases: object

copy()[source]
diagvec_mul_jac(a)[source]
full_jac()[source]
jac_mul_diagvec(a)[source]
class Array(values, name=None)[source]

Bases: Operator

AD representation of a constant numpy array.

For sparse matrices, use Matrix instead. For time-dependent arrays see TimeDependentArray.

This is a shallow wrapper around the real array; it is needed to combine the array with other types of AD operators.

Parameters

  • values (np.ndarray) – Numpy array to be represented.

  • name (Optional[str]) –

parse(mdg)[source]

See Operator.parse().

Returns

The wrapped array.

Parameters

mdg (MixedDimensionalGrid) –

Return type

ndarray

class BiotAd(keyword, subdomains, flow_keyword='flow')[source]

Bases: Discretization

Ad wrapper around the Biot discretization class.

For description of the method, we refer to the standard Biot class.

Parameters

  • keyword (str) –

  • subdomains (List[pp.Grid]) –

  • flow_keyword (str) –

class BiotStabilizationAd(keyword, subdomains)[source]

Bases: Discretization

Parameters

  • keyword (str) –

  • subdomains (List[pp.Grid]) –

class BoundaryCondition(keyword, subdomains, name=None)[source]

Bases: Operator

Wrapper class for Ad representations of boundary conditions for a given keyword.

Parameters

  • keyword (str) –

  • subdomains (list[pp.Grid]) –

  • name (Optional[str]) –

parse(mdg)[source]

Convert the Ad expression into numerical values for the boundary conditions, in the form of an np.ndarray concatenated for all subdomains.

Parameters

mdg (pp.MixedDimensionalGrid) – Mixed-dimensional grid. The boundary condition will be taken from the data dictionaries with the relevant keyword.

Returns

Value of boundary conditions.

Return type

np.ndarray

class ColoumbContactAd(keyword, interfaces)[source]

Bases: Discretization

Parameters

  • keyword (str) –

  • interfaces (List[pp.MortarGrid]) –

class ConstantFunction(name, values)[source]

Bases: AbstractFunction

Function representing constant, scalar values with no dependencies and ergo a zero Jacobian.

It still has to be called though since it fulfills the notion of a ‘function’.

Parameters

  • values (np.ndarray) – constant values per cell.

  • name (str) –

get_jacobian(*args)[source]

Note

The return value is not a sparse matrix as imposed by the parent method signature, but a zero. Numerical operations with a zero always works with any numeric formats in numpy, scipy and PorePy’s AD framework. Since the constant function (most likely) gets no arguments passed, we have no way of knowing the necessary shape for a zero matrix. Hence scalar.

Returns: the trivial derivative of a constant.

Parameters

args (Ad_array) –

Return type

spmatrix

get_values(*args)[source]
Returns

The values passed at instantiation.

Parameters

args (Ad_array) –

Return type

ndarray

class ContactTractionAd(keyword, interfaces, low_dim_subdomains)[source]

Bases: Discretization

Parameters

  • keyword (str) –

  • interfaces (List[pp.MortarGrid]) –

  • low_dim_subdomains (List[pp.Grid]) –

class DiagonalJacobianFunction(func, name, multipliers, array_compatible=False)[source]

Bases: AbstractJacobianFunction

Approximates the Jacobian of the function using identities and scalar multipliers per dependency.

Parameters

  • multipliers (float | list[float]) – scalar multipliers for the identity blocks in the Jacobian, per dependency of func. The order in multipliers is expected to match the order of AD operators passed to the call of this function.

  • func (Callable) –

  • name (str) –

  • array_compatible (bool) –

get_jacobian(*args)[source]

The approximate Jacobian consists of identity blocks times scalar multiplier per every function dependency.

Parameters

args (Ad_array) –

Return type

spmatrix

class DifferentiableFVAd(subdomains, mdg, base_discr, system_manager, permeability_function, permeability_argument, potential, keyword)[source]

Bases: object

This class represents the application of the product and chain rule of the flux expression

q = T(k(u)) * p

Where the transmissibility matrix T is a function of the cell permeability k, which again is a function of a primary variable u, while p is the potential (pressure). The chain rule applied to this expression reads

dq = p * dT/dk * dk/du * du + T * dp

The transmissibility matrix can be computed from a Tpfa or Mpfa discretization, or (in principle) any other finite volume method. The derivative of the transmissibilities, dT/dk, is approximated with a two-point flux stencil.

If vector sources are included, p should be replaced by (p - dist * vector_source), with dist the distance roughly corresponding to the inverse of the gradient.

Parameters

  • subdomains (List[pp.Grid]) –

  • mdg (pp.MixedDimensionalGrid) –

  • base_discr (Union[pp.ad.MpfaAd, pp.ad.TpfaAd]) –

  • system_manager (pp.ad.EquationSystem) –

  • permeability_function (Callable[[pp.ad.Variable], pp.ad.Ad_array]) –

  • permeability_argument (pp.ad.Variable) –

  • potential (pp.ad.Variable) –

  • keyword (str) –

bound_pressure_face()[source]

Boundary face contribution to pressure reconstruction.

Returns

pp.ad.Operator representing the value and jacobian of bound_pressure_face,

which is used to reconstruct the pressure trace on the boundary (see FV discretizations).

Return type

Operator

flux()[source]

Flux from potential, BCs and vector sources.

Returns

pp.ad.Operator representing the flux.

Return type

Operator

class Discretization[source]

Bases: ABC

General/utility methods for AD discretization classes.

The init of the children classes below typically calls wrap_discretization and has arguments including subdomains or interfaces and keywords for parameter and possibly matrix storage.

class DivUAd(keyword, subdomains, mat_dict_keyword)[source]

Bases: Discretization

Parameters

  • keyword (str) –

  • subdomains (List[pp.Grid]) –

  • mat_dict_keyword (str) –

class Divergence(subdomains, dim=1, name=None)[source]

Bases: Operator

Wrapper class for Ad representations of divergence operators.

Parameters

  • subdomains (list[pp.Grid]) –

  • dim (int) –

  • name (Optional[str]) –

parse(mdg)[source]

Convert the Ad expression into a divergence operators on all relevant subdomains, represented as a sparse block matrix.

Parameters

mdg (pp.MixedDimensionalGrid) – Not used, but needed for compatibility with the general parsing method for Operators.

Returns

Block matrix representation of a divergence operator on

multiple subdomains.

Return type

sps.spmatrix

class EquationManager(mdg, dof_manager, equations=None, secondary_variables=None)[source]

Bases: object

Representation of a set of equations specified on Ad form.

The equations are tied to a specific MixedDimensionalGrid, with variables fixed in a corresponding DofManager.

Parameters

  • mdg (pp.MixedDimensionalGrid) –

  • dof_manager (pp.DofManager) –

  • equations (Optional[Dict[str, 'pp.ad.Operator']]) –

  • secondary_variables (Optional[Sequence['pp.ad.Variable']]) –

mdg

Mixed-dimensional grid on which this EquationManager operates.

Type

pp.MixedDimensionalGrid

dof_manager

Degree of freedom manager used for this EquationManager.

Type

pp.DofManager

equations

Equations assigned to this EquationManager. can be expanded by direct addition to the list.

Type

List of Ad Operators

variables

Mapping from subdomains or grid tuples (interfaces) to Ad variables. These are set at initialization from the MixedDimensionalGrid, and should not be changed later.

Type

Dict

secondary_variables

List of variables that are secondary, that is, their derivatives will not be included in the Jacobian matrix. Variables will be represented on atomic form, that is, mixed-dimensional variables are unravelled. Secondary variables act as a filter during assembly, that is, they do not impact the ordering or treatment of variables.

Type

List of Ad Variables

row_block_indices_last_assembled

Row indices for the start of blocks corresponding to different equations in the last assembled system. The last item in the array is the total number of rows, so that row indices for block i can be recovered by np.arange(row_bl..[i], row_bl..[i+1]). The user must relate the indices to equations (either in self.equations or the equation list given to the relevant assembly method). This information is intended for diagnostic usage.

Type

np.ndarray

assemble(state=None)[source]

Assemble Jacobian matrix and residual vector with respect to the current state represented in self.mdg.

Derivatives for secondary variables are not included in the Jacobian matrix.

Parameters

state (np.ndarray, optional) – State vector to assemble from. If not provided, the default behavior of pp.ad.Expression.to_ad() will be followed.

Returns

Jacobian matrix corresponding to the current variable state,

as found in self.mdg. The ordering of the equations is determined by the ordering in self.equations (for rows) and self.dof_manager (for columns).

np.ndarray: Residual vector corresponding to the current variable state,

as found in self.mdg. Scaled with -1 (moved to rhs).

Return type

sps.spmatrix

assemble_schur_complement_system(primary_equations, primary_variables, inverter)[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 a reordering of the linearized system into

[A_pp, A_ps [x_p = [b_p

A_sp, A_ss] x_s] b_s]

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

(A_pp - A_ps * inv(A_ss) * A_sp) * x_p = b_p - A_ps * inv(A_pp) * 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. The simplest option is a lambda function on the form:

inverter = lambda A: sps.csr_matrix(np.linalg.inv(A.A))

but depending on A (size and sparsity pattern), this can be costly in terms of computational time and memory.

The method can be used e.g. for splitting between primary and secondary variables, where the latter can be efficiently eliminated (for instance, they contain no spatial derivatives).

Parameters

  • primary_equations (Sequence of str) – Equations to be assembled, specified as keys in self.equations. Should have length > 0.

  • primary_variables (Sequence of Variables) – Variables to be assembled. Should have length > 0.

  • inverter (Callable) – Method to compute the inverse of the matrix A_ss.

Returns

Jacobian matrix corresponding to the current variable state,

as found in self.mdg, for the specified equations and variables.

np.ndarray: Residual vector corresponding to the current variable state,

as found in self.mdg, for the specified equations and variables. Scaled with -1 (moved to rhs).

Return type

sps.spmatrix

assemble_subsystem(eq_names=None, variables=None)[source]

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

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

Parameters

  • eq_names (Sequence of str, optional) – Equations to be assembled, specified as keys in self.equations. If not provided (None), all equations known to this EquationManager will be included.

  • variables (Sequence of Variables, optional) – Variables to be assembled. If not provided (None), all variables known to this EquationManager will be included. If a secondary variable is specified, this will be included in the returned system.

Returns

Jacobian matrix corresponding to the current variable state,

as found in self.mdg, for the specified equations and variables.

np.ndarray: Residual vector corresponding to the current variable state,

as found in self.mdg, for the specified equations and variables. Scaled with -1 (moved to rhs).

NOTE: The ordering of columns in the system are defined by the order of the

variables specified in DofManager. For the rows, no corresponding global ordering of equations exists, and the rows will therefore be organized by the ordering in the parameter eq_names.

Return type

sps.spmatrix

discretize(mdg)[source]

Loop over all discretizations in self.equations, find all unique discretizations 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

mdg (pp.MixedDimensionalGrid) – Mixed-dimensional grid from which parameters etc. will be taken and where discretization matrices will be stored.

Return type

None

merge_variables(grid_var)[source]

Concatenate a variable defined over several subdomains or interfaces.

The mixed-dimensional variable can be used to define mathematical operations on multiple subdomains simultaneously (provided it is combined with other operators defined on the same subdomains).

NOTE: mixed-dimensional variables are assigned unique ids (see documentation of Variable and MixedDimensionalVariable), thus two MixedDimensionalVariables will have different ids even if they represent the same combination of subdomains and variables. This does not impact the parsing of the variables into numerical values.

Parameters

grid_var (Sequence[Tuple[Union[Grid, MortarGrid], str]]) – Tuple containing first a (Mortar)grid representing the subdomain or interface and second the name of the variable.

Returns

Joint representation of the variable on the

specified subdomains.

Return type

pp.ad.MixedDimensionalVariable

name_and_assign_equations(equation_dictionary)[source]

Convenience method for assigning and naming equations.

Parameters

equation_dictionary (Dict) – Dictionary containing name: equation pairs.

Return type

None

subsystem_equation_manager(eq_names, variables)[source]

Extract an EquationManager for a subset of variables and equations. In effect, this produce a nonlinear subsystem.

Parameters

  • eq_names (Sequence of str) – Equations assigned to the new EquationManager, specified as keys in self.equations.

  • variables (Sequence of Variables) – Variables for which the new EquationManager is defined.

Returns

System of nonlinear equations. The ordering of the

equations in the subsystem will be the same as in the original EquationManager (i.e. self), disregarding equations not included in the subset. Variables that were excluded are added to the set of secondary_variables in the new EquationManager.

Return type

EquationManager

variable(grid_like, variable)[source]

Get a variable for a specified grid or interface between grids, that is a mortar grid.

Subsequent calls of this method with the same grid and variable will return references to the same variable.

Returns

Ad representation of a variable.

Return type

pp.ad.Variable

Parameters
variable_state(grid_var, state)[source]
Parameters
Return type

List[ndarray]

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 (pp.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.

class Function(func, name, array_compatible=True)[source]

Bases: AbstractFunction

Ad representation of an analytically given function, where it is expected that passing Ad_arrays directly to func will return the proper result.

Here the values and the Jacobian are obtained exactly by the AD framework.

The intended use is as a wrapper for operations on pp.ad.Ad_array objects, in forms which are not directly or easily expressed by the rest of the Ad framework.

Note

This is a special case where the abstract methods for getting values and the Jacobian are formally implemented but never used by the AD framework. A separate operation called evaluate is implemented instead, which simply feeds the AD arrays to func.

Parameters

  • func (Callable) –

  • name (str) –

  • array_compatible (bool) –

get_jacobian(*args)[source]
Parameters

args (Ad_array) –

Return type

ndarray

get_values(*args)[source]
Parameters

args (Ad_array) –

Return type

ndarray

class Geometry(subdomains, nd, name=None, matrix_names=None)[source]

Bases: Operator

Wrapper class for Ad representations of grids.

Parameters

  • subdomains (list[pp.Grid]) –

  • nd (int) –

  • name (Optional[str]) –

  • matrix_names (Optional[list[str]]) –

cell_volumes

Diagonal ad matrix of cell volumes.

Type

pp.ad.Matrix

face_areas

Diagonal ad matrix of face areas.

Type

pp.ad.Matrix

nd

Ambient/highest dimension of the mixed-dimensional grid.

Type

int

FIXME: Implement parse??

basis(dim=None)[source]

Return a cell-wise basis for all subdomains.

Parameters

dim (int, optional) – Dimension of the base. Defaults to the dimension of the Geometry.

Returns

Array of dim pp.ad.Matrix, each of which is represents a basis function.

Return type

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

e_i(i, dim=None)[source]

Return a cell-wise basis function for all subdomains.

Parameters

  • dim (int) – Dimension of the functions.

  • i (int) – Index of the basis function. Note: Counts from 0.

Returns

Ad representation of a matrix with the basis functions as

columns.

Return type

pp.ad.Matrix

class GradPAd(keyword, subdomains)[source]

Bases: Discretization

Parameters

  • keyword (str) –

  • subdomains (List[pp.Grid]) –

class InterpolatedFunction(func, name, min_val, max_val, npt, order=1, preval=False, array_compatible=False)[source]

Bases: AbstractFunction

Represents the passed function as an interpolation of chosen order on a cartesian, uniform grid.

The image of the function is expected to be of dimension 1, while the domain can be multidimensional.

Note

All vector-valued ndarray arguments are assumed to be column vectors. Each row-entry represents a value for an argument of func in respective order.

Parameters

  • min_val (np.ndarray) – lower bounds for the domain of func.

  • max_val (np.ndarray) – upper bound for the domain.

  • npt (np.ndarray) – number of interpolation points per dimension of the domain.

  • order (int) – Order of interpolation. Supports currently only linear order.

  • preval (optional) – If True, pre-evaluates the values of the function at the points of interpolation and stores them. If False, evaluates them if necessary. Influences the runtime. Defaults to False.

  • func (Callable) –

  • name (str) –

  • array_compatible (bool) –

get_jacobian(*args)[source]
Parameters

args (Ad_array) –

Return type

spmatrix

get_values(*args)[source]
Parameters

args (Ad_array) –

Return type

ndarray

class MassMatrixAd(keyword, subdomains)[source]

Bases: Discretization

Parameters

  • keyword (str) –

  • subdomains (List[pp.Grid]) –

class Matrix(mat, name=None)[source]

Bases: Operator

Ad representation of a sparse matrix.

For dense matrices, use Array instead.

This is a shallow wrapper around the real matrix; it is needed to combine the matrix with other types of Ad objects.

Parameters

  • mat (sps.spmatrix) – Sparse matrix to be wrapped as an AD operator.

  • name (Optional[str]) – Name of this operator

parse(mdg)[source]

See Operator.parse().

Returns

The wrapped matrix.

Parameters

mdg (MixedDimensionalGrid) –

Return type

spmatrix

transpose()[source]

Returns an AD operator representing the transposed matrix.

Return type

Matrix

property T: Matrix

Shorthand for transpose.

shape

Shape of the wrapped matrix.

class MixedDimensionalVariable(variables)[source]

Bases: Variable

Ad representation of a collection of variables that individually live on separate subdomains or interfaces, but treated jointly in the mixed-dimensional sense.

Conversion of the variables into numerical value should be done with respect to the state of an array; see Operator.evaluate(). Therefore, the MergedVariable does not implement the method Operator.parse().

Parameters

variables (list[Variable]) – List of variables to be merged. Should all have the same name.

copy()[source]

Copy the mixed-dimensional variable.

Returns

A shallow copy should be sufficient here; the attributes are not expected to change.

Return type

MixedDimensionalVariable

copy_common_sub_tags()[source]

Copy any shared tags from the sub variables to this variable.

Only tags with identical values are copied. Thus, the md variable can “trust” that its tags are consistent with all sub variables.

Return type

None

previous_iteration()[source]

Return a representation of this mixed-dimensional variable on the previous iteration.

Returns

A representation of this merged variable on the previous iteration, with its prev_iter attribute set to True

Return type

MixedDimensionalVariable

previous_timestep()[source]

Return a representation of this mixed-dimensional variable on the previous time step.

If the md-variable is already defined on the previous time step, return itself.

Returns

A representation of this merged variable on the previous time iteration, with its prev_iter attribute set to True.

Return type

MixedDimensionalVariable

size()[source]

Returns the total size of the mixed-dimensional variable by summing the sizes of sub-variables.

Return type

int

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

A tuple of all domains on which the atomic sub-variables are defined.

id: int

ID counter. Used to identify variables during operator parsing.

original_variable: MixedDimensionalVariable

The original variable, if this variable is a copy of another variable.

This attribute is used by the methods Variable.previous_timestep() and Variable.previous_iteration() to keep a link to the original variable.

prev_iter: bool

Flag indicating if the variable represents the state at the previous iteration.

prev_time: bool

Flag indicating if the variable represents the state at the previous time step.

sub_vars

List of sub-variables passed at instantiation, each defined on a separate domain.

class MortarProjections(mdg, subdomains, interfaces, dim=1)[source]

Bases: Operator

Wrapper class to generate projections to and from MortarGrids.

Parameters

  • mdg (pp.MixedDimensionalGrid) –

  • subdomains (list[pp.Grid]) –

  • interfaces (list[pp.MortarGrid]) –

  • dim (int) –

mortar_to_primary_int

Matrix of projections from the mortar grid to the primary grid. Intended for extensive quantities (so fluxes). Represented as an Ad Matrix operator.

Type

pp.ad.Matrix

mortar_to_primary_avg

Matrix of projections from the mortar grid to the primary grid. Intended for intensive quantities (so pressures). Represented as an Ad Matrix operator.

Type

pp.ad.Matrix

primary_to_mortar_int

Matrix of projections from the primary grid to the mortar grid. Intended for extensive quantities (so fluxes). Represented as an Ad Matrix operator.

Type

pp.ad.Matrix

primary_to_mortar_avg

Matrix of projections from the primary grid to the mortar grid. Intended for intensive quantities (so pressures). Represented as an Ad Matrix operator.

Type

pp.ad.Matrix

mortar_to_secondary_int

Matrix of projections from the mortar grid to the secondary grid. Intended for extensive quantities (so fluxes). Represented as an Ad Matrix operator.

Type

pp.ad.Matrix

mortar_to_secondary_avg

Matrix of projections from the mortar grid to the secondary grid. Intended for intensive quantities (so pressures). Represented as an Ad Matrix operator.

Type

pp.ad.Matrix

secondary_to_mortar_int

Matrix of projections from the secondary grid to the mortar grid. Intended for extensive quantities (so fluxes). Represented as an Ad Matrix operator.

Type

pp.ad.Matrix

secondary_to_mortar_avg

Matrix of projections from the secondary grid to the mortar grid. Intended for intensive quantities (so pressures). Represented as an Ad Matrix operator.

Type

pp.ad.Matrix

sign_of_mortar_sides

Matrix representation that assigns signs to two mortar sides. Needed to implement a jump operator in contact mechanics.

Type

pp.Ad.Matrix

class MpfaAd(keyword, subdomains)[source]

Bases: Discretization

Parameters

  • keyword (str) –

  • subdomains (List[pp.Grid]) –

class MpsaAd(keyword, subdomains)[source]

Bases: Discretization

Parameters

  • keyword (str) –

  • subdomains (List[pp.Grid]) –

class Operator(name=None, subdomains=None, interfaces=None, tree=None)[source]

Bases: object

Parent class for all AD operators.

Objects of this class are not meant to be initiated directly, rather the various subclasses should be used. Instances of this class will still be created when subclasses are combined by Operator.Operations.

Contains a tree structure of child operators for the recursive forward evaluation.

Provides overload functions for basic arithmetic operations.

Parameters

  • name (Optional[str]) – Name of this operator. Used for string representations

  • subdomains (optional) – List of subdomains on which the operator is defined. Will be empty for operators not associated with any subdomains. Defaults to None (converted to empty list).

  • interfaces (optional) – List of interfaces in the mixed-dimensional grid on which the operator is defined. Will be empty for operators not associated with any interface. Defaults to None (converted to empty list).

  • tree (Optional[Tree]) –

class Operations(value)

Bases: Enum

Object representing all supported operations by the operator class.

Used to construct the operator tree and identify Operator.Operations.

add
approximate
div
evaluate
mul
pow
sub
void
discretize(mdg)[source]

Perform discretization operation on all discretizations identified in the tree of this operator, using data from mdg.

IMPLEMENTATION NOTE: The discretizations was identified at initialization of Expression - it is now done here to accommodate updates (?) and

Parameters

mdg (MixedDimensionalGrid) –

Return type

None

evaluate(system_manager, state=None)[source]

Evaluate the residual and Jacobian matrix for a given state.

Parameters

  • system_manager (pp.ad.EquationSystem | pp.DofManager) – Used to represent the problem. Will be used to parse the sub-operators that combine to form this operator.

  • state (optional) – State vector for which the residual and its derivatives should be formed. If not provided, the state will be pulled from the previous iterate (if this exists), or alternatively from the state at the previous time step.

Returns

A representation of the residual and Jacobian in form of an AD Array. Note that the Jacobian matrix need not be invertible, or even square; this depends on the operator.

is_leaf()[source]

Check if this operator is a leaf in the tree-representation of an expression.

Note that this implies that the method parse() is expected to be implemented.

Returns

True if the operator has no children.

Return type

bool

parse(mdg)[source]

Translate the operator into a numerical expression.

Subclasses that represent atomic operators (leaves in a tree-representation of an operator) should override this method to return e.g. a number, an array or a matrix. This method should not be called on operators that are formed as combinations of atomic operators; such operators should be evaluated by the method evaluate().

Parameters

mdg (MixedDimensionalGrid) – Mixed-dimensional grid on which this operator is to be parsed.

Returns

A numerical format representing this operator;s values on given domain.

Return type

Any

previous_timestep()[source]

Return an operator that represents the value of this operator at the previous timestep.

The operator tree at the previous time step is created as a shallow copy, and will thus be identical to the original operator, except that all time dependent operators are evaluated at the previous time step.

Returns

A copy of self, with all time dependent operators evaluated at the previous time step.

Return type

Operator

set_name(name)[source]

Reset this object’s name originally passed at instantiation.

Parameters

name (str) – the new name to be assigned.

Return type

None

viz()[source]

Draws a visualization of the operator tree that has this operator as its root.

interfaces: list[porepy.grids.mortar_grid.MortarGrid]

List of interfaces on which the operator is defined, passed at instantiation.

Will be empty for operators not associated with specific interfaces.

property name: str

The name given to this variable.

subdomains: list[porepy.grids.grid.Grid]

List of subdomains on which the operator is defined, passed at instantiation.

Will be empty for operators not associated with specific subdomains.

class ParameterArray(param_keyword, array_keyword, subdomains=None, interfaces=None, name=None)[source]

Bases: Operator

Extract an array from the parameter dictionaries for a given set of subdomains.

Can be used to implement sources, and general arrays to be picked from the parameter array (and thereby could be changed during the simulation, without having to redefine the abstract Ad representation of the equations).

Parameters
parse(mdg)[source]

Convert the Ad expression into numerical values for the scalar sources, in the form of an np.ndarray concatenated for all subdomains.

Parameters

mdg (pp.MixedDimensionalGrid) – Mixed-dimensional grid. The boundary condition will be taken from the data dictionaries with the relevant keyword.

Returns

Value of boundary conditions.

Return type

np.ndarray

interfaces: list[pp.MortarGrid]

List of interfaces on which the operator is defined, passed at instantiation.

Will be empty for operators not associated with specific interfaces.

subdomains: list[pp.Grid]

List of subdomains on which the operator is defined, passed at instantiation.

Will be empty for operators not associated with specific subdomains.

class ParameterMatrix(param_keyword, array_keyword, subdomains=None, interfaces=None, name=None)[source]

Bases: ParameterArray

Extract a matrix from the parameter dictionaries for a given set of subdomains.

Typical use: Parameters which are left multiplied with an ad expression. Note that

array parameters are represented by one diagonal matrix for each grid.

Parameters

  • param_keyword (str) –

  • array_keyword (str) –

  • subdomains (Optional[list[pp.Grid]]) –

  • interfaces (Optional[list[pp.MortarGrid]]) –

  • name (Optional[str]) –

parse(mdg)[source]

Convert the Ad expression into numerical values for the scalar sources, in the form of an np.ndarray concatenated for all subdomains.

Parameters

mdg (pp.MixedDimensionalGrid) – Mixed-dimensional grid. The boundary condition will be taken from the data dictionaries with the relevant keyword.

Returns

Value of boundary conditions.

Return type

sps.spmatrix

interfaces: list[pp.MortarGrid]

List of interfaces on which the operator is defined, passed at instantiation.

Will be empty for operators not associated with specific interfaces.

subdomains: list[pp.Grid]

List of subdomains on which the operator is defined, passed at instantiation.

Will be empty for operators not associated with specific subdomains.

class RegularizedHeaviside(regularization)[source]

Bases: object

Parameters

regularization (Callable) –

class RobinCouplingAd(keyword, interfaces)[source]

Bases: Discretization

Parameters

  • keyword (str) –

  • interfaces (List[pp.MortarGrid]) –

class Scalar(value, name=None)[source]

Bases: Operator

Ad representation of a scalar.

This is a shallow wrapper around a real scalar. It may be useful to combine the scalar with other types of Ad objects.

NOTE: Since this is a wrapper around a Python immutable, copying a Scalar will effectively create a deep copy, i.e., changes in the value of one Scalar will not be reflected in the other. This is in contrast to the behavior of the other Ad objects.

Parameters

  • value (float) –

  • name (Optional[str]) –

parse(mdg)[source]

See Operator.parse().

Returns

The wrapped number.

Parameters

mdg (MixedDimensionalGrid) –

Return type

float

interfaces: list[porepy.grids.mortar_grid.MortarGrid]

List of interfaces on which the operator is defined, passed at instantiation.

Will be empty for operators not associated with specific interfaces.

subdomains: list[porepy.grids.grid.Grid]

List of subdomains on which the operator is defined, passed at instantiation.

Will be empty for operators not associated with specific subdomains.

class SubdomainProjections(subdomains, dim=1)[source]

Bases: Operator

Wrapper class for generating projection to and from subdomains.

One use case in when variables are defined on only some subdomains.

The class should be used through the methods {cell, face}_{projection, restriction}.

See also MortarProjections for projections to and from mortar subdomains.

Parameters
cell_prolongation(subdomains)[source]

Construct prolongation from subdomain to global cell quantities.

Parameters

subdomains (List of pp.Grid) – One or several subdomains to which the prolongation should apply.

Returns

Matrix operator (in the Ad sense) that represent the

prolongation.

Return type

pp.ad.Matrix

cell_restriction(subdomains)[source]

Construct restrictions from global to subdomain cell quantities.

Parameters

subdomains (List of pp.Grid) – One or several subdomains to which the projection should apply.

Returns

Matrix operator (in the Ad sense) that represents the

projection.

Return type

pp.ad.Matrix

face_prolongation(subdomains)[source]

Construct prolongation from subdomain to global face quantities.

Parameters

subdomains (List of pp.Grid) – One or several subdomains to which the prolongation should apply.

Returns

Matrix operator (in the Ad sense) that represent the

prolongation.

Return type

pp.ad.Matrix

face_restriction(subdomains)[source]

Construct restrictions from global to subdomain face quantities.

Parameters

subdomains (List of pp.Grid) – One or several subdomains to which the projection should apply.

Returns

Matrix operator (in the Ad sense) that represent the

projection.

Return type

pp.ad.Matrix

interfaces: list[pp.MortarGrid]

List of interfaces on which the operator is defined, passed at instantiation.

Will be empty for operators not associated with specific interfaces.

subdomains: list[pp.Grid]

List of subdomains on which the operator is defined, passed at instantiation.

Will be empty for operators not associated with specific subdomains.

class TimeDependentArray(name, subdomains=None, interfaces=None, previous_timestep=False)[source]

Bases: Array

An Ad-wrapper around a time-dependent numpy array.

The array is tied to a MixedDimensionalGrid, and is distributed among the data dictionaries associated with subdomains and interfaces. The array values are stored in data[pp.STATE][pp.ITERATE][self._name] for the current time and data[pp.STATE][self._name] for the previous time.

The array can be differentiated in time using pp.ad.dt().

The intended use is to represent time-varying quantities in equations, e.g., source terms. Future use will also include numerical values of boundary conditions, however, this is pending an update to the model classes.

Parameters

  • name (str) – Name of the variable. Should correspond to items in data[pp.STATE].

  • subdomains (list[porepy.grids.grid.Grid]) – Subdomains on which the array is defined. Defaults to None.

  • interfaces (list[porepy.grids.mortar_grid.MortarGrid]) – Interfaces on which the array is defined. Defaults to None. Exactly one of subdomains and interfaces must be non-empty.

  • previous_timestep (bool) – Flag indicating if the array should be evaluated at the previous time step.

previous_timestep[source]

If True, the array will be evaluated using data[pp.STATE] (data being the data dictionaries for subdomains and interfaces), if False, data[pp.STATE][pp.ITERATE] is used.

Return type

TimeDependentArray

Raises

ValueError – If either none of, or both of, subdomains and interfaces are empty.

Parameters
parse(mdg)[source]

Convert this array into numerical values.

The numerical values will be picked from the representation of the array in data[pp.STATE][pp.ITERATE] (where data is the data dictionary of the subdomains or interfaces of this Array), or, if self.prev_time = True, from data[pp.STATE].

Parameters

mdg (MixedDimensionalGrid) – Mixed-dimensional grid.

Returns

A numpy ndarray containing the numerical values of this array.

Return type

ndarray

previous_timestep()[source]
Returns

This array represented at the previous time step.

Return type

TimeDependentArray

interfaces: list[porepy.grids.mortar_grid.MortarGrid]

List of interfaces on which the operator is defined, passed at instantiation.

Will be empty for operators not associated with specific interfaces.

prev_time: bool

If True, the array will be evaluated using data[pp.STATE] (data being the data dictionaries for subdomains and interfaces).

If False, data[pp.STATE][pp.ITERATE] is used.

subdomains: list[porepy.grids.grid.Grid]

List of subdomains on which the operator is defined, passed at instantiation.

Will be empty for operators not associated with specific subdomains.

class TpfaAd(keyword, subdomains)[source]

Bases: Discretization

Parameters

  • keyword (str) –

  • subdomains (List[pp.Grid]) –

class Trace(subdomains, dim=1, name=None)[source]

Bases: Operator

Wrapper class for Ad representations of trace operators and their inverse, that is, mappings between grid cells and faces.

NOTE: The mapping will hit both boundary and interior faces, so the values to be mapped should be carefully filtered (e.g. by combining it with a mortar mapping).

The mapping does not alter signs of variables, that is, the direction of face normal vectors is not accounted for.

Parameters
trace

Matrix of trace projections from cells to faces.

Type

pp.ad.Matrix

inv_trace

Matrix of trace projections from faces to cells.

Type

pp.ad.Matrix

interfaces: list[pp.MortarGrid]

List of interfaces on which the operator is defined, passed at instantiation.

Will be empty for operators not associated with specific interfaces.

subdomains: list[pp.Grid]

List of subdomains on which the operator is defined, passed at instantiation.

Will be empty for operators not associated with specific subdomains.

class UpwindAd(keyword, subdomains)[source]

Bases: Discretization

Parameters

  • keyword (str) –

  • subdomains (List[pp.Grid]) –

class UpwindCouplingAd(keyword, interfaces)[source]

Bases: Discretization

Parameters

  • keyword (str) –

  • interfaces (List[pp.MortarGrid]) –

class Variable(name, ndof, domain, tags=None, previous_timestep=False, previous_iteration=False)[source]

Bases: Operator

AD operator representing a variable defined on a single grid or mortar grid.

For combinations of variables on different subdomains, see MergedVariable.

Conversion of the variable into numerical value should be done with respect to the state of an array; see Operator.evaluate(). Therefore, the variable does not implement the method Operator.parse().

A variable is associated with either a grid or an interface. Therefore it is assumed that either subdomains or interfaces is passed as an argument.

Parameters

  • name (str) – Variable name.

  • ndof (dict[Literal['cells', 'faces', 'nodes'], int]) – Number of dofs per grid element. Valid keys are cells, faces and nodes.

  • subdomains (length=1) – List containing a single grid.

  • interfaces (length=1) – List containing a single mortar grid.

  • num_cells – Number of cells in the grid. Only relevant if this is an interface variable.

  • domain (GridLike) –

  • tags (Optional[dict[str, Any]]) –

  • previous_timestep (bool) –

  • previous_iteration (bool) –

previous_iteration()[source]
Returns

A representation of this variable on the previous time iteration, with its prev_iter attribute set to True.

Return type

Variable

previous_timestep()[source]

Return a representation of this variable on the previous time step.

If the variable is already a representation of the previous time step, the method returns itself.

Returns

A representation of this variable at the previous time step, with its prev_time attribute set to True.

Return type

Variable

set_name(name)[source]
Raises

RuntimeError – Variables must not be re-named once defined, since the name is used as an identifier.

Parameters

name (str) –

Return type

None

size()[source]

Returns the total number of dofs this variable has.

Return type

int

property domain: Union[Grid, MortarGrid]

The grid or mortar grid on which this variable is defined.

id: int

ID counter. Used to identify variables during operator parsing.

interfaces: list[porepy.grids.mortar_grid.MortarGrid]

List of interfaces on which the operator is defined, passed at instantiation.

Will be empty for operators not associated with specific interfaces.

original_variable: Variable

The original variable, if this variable is a copy of another variable.

This attribute is used by the methods Variable.previous_timestep() and Variable.previous_iteration() to keep a link to the original variable.

prev_iter: bool

Flag indicating if the variable represents the state at the previous iteration.

prev_time: bool

Flag indicating if the variable represents the state at the previous time step.

subdomains: list[porepy.grids.grid.Grid]

List of subdomains on which the operator is defined, passed at instantiation.

Will be empty for operators not associated with specific subdomains.

property tags: dict[str, Any]

A dictionary of tags associated with this variable.

class WellCouplingAd(keyword, interfaces)[source]

Bases: Discretization

Parameters

  • keyword (str) –

  • interfaces (List[pp.MortarGrid]) –

abs(var)[source]
arccos(var)[source]
arccosh(var)[source]
arcsin(var)[source]
arcsinh(var)[source]
arctan(var)[source]
arctanh(var)[source]
characteristic_function(tol, var)[source]

Characteristic function of an ad variable.

Returns 1 if var.val is within absolute tolerance = tol of zero. The derivative is set to zero independent of var.val.

Note

See module level documentation on how to wrap functions like this in ad.Function.

Parameters

  • tol (float) – Absolute tolerance for comparison with 0 using np.isclose.

  • var (Ad_array) – Ad operator (variable or expression).

Returns

The characteristic function of var with appropriate val and jac attributes.

cos(var)[source]
cosh(var)[source]
dt(op, time_step)[source]

Approximate the time-derivative of an operator tree.

This time derivative is obtained by taking a first-order finite difference.

The operator tree at the previous time step is created as a shallow copy, and will thus be identical to the original operator, except that all time dependent operators are evaluated at the previous time step.

If the time-dependent quantity q is already evaluated at the previous time step, its derivative will be defined as (q(time=n-1) - q(time=n-1)) / dt = 0.

Parameters

  • op (Operator) – Operator tree to be differentiated.

  • dt – Size of time step.

  • time_step (Scalar) –

Returns

A first-order approximation of the time derivative of op.

Return type

Operator

exp(var)[source]
heaviside(var, zerovalue=0.5)[source]
Parameters

zerovalue (float) –

heaviside_smooth(var, eps=1e-3)[source]

Smooth (regularized) version of the Heaviside function.

Note

The analytical expression for the smooth version Heaviside function reads:

H_eps(x) = (1/2) * (1 + (2/pi) * arctan(x/eps)),

with its derivative smoothly approximating the Dirac delta function:

d(H(x))/dx = delta_eps = (1/pi) * (eps / (eps^2 + x^2)).

Reference: https://ieeexplore.ieee.org/document/902291

Parameters

  • var – Input array.

  • eps (optional) – Regularization parameter. The function will converge to the Heaviside function in the limit when eps --> 0. The default is 1e-3.

Returns

Regularized heaviside function (and its Jacobian if applicable) in form of a Ad_array or ndarray (depending on the input).

initAdArrays(variables)[source]
l2_norm(dim, var)[source]

L2 norm of a vector variable.

For the example of dim=3 components and n vectors, the ordering is assumed to be [u0, v0, w0, u1, v1, w1, ..., un, vn, wn]

Vectors satisfying ui=vi=wi=0 are assigned zero entries in the jacobi matrix

Note

See module level documentation on how to wrap functions like this in ad.Function.

Parameters

  • dim (int) – Dimension, i.e. number of vector components.

  • var (Ad_array) – Ad operator (variable or expression) which is argument of the norm function.

Returns

The norm of var with appropriate val and jac attributes.

Return type

Ad_array

log(var)[source]
maximum(var_0, var_1)[source]

Ad maximum function represented as an Ad_array.

The arguments can be either Ad_arrays or ndarrays, this duality is needed to allow for parsing of operators that can be taken at the current iteration (in which case it will parse as an Ad_array) or at the previous iteration or time step (in which case it will parse as a numpy array).

Parameters

  • var_0 (pp.ad.Ad_array) – First argument to the maximum function.

  • var_1 (pp.ad.Ad_array | np.ndarray) – Second argument.

  • scalar (If one of the input arguments is) –

  • used. (broadcasting will be) –

Returns

The maximum of the two arguments, taken element-wise in the arrays. The return type is Ad_array if at least one of the arguments is an Ad_array, otherwise it is an ndarray. If an Ad_array is returned, the Jacobian is computed according to the maximum values of the Ad_arrays (so if element i of the maximum is picked from var_0, row i of the Jacobian is also picked from the Jacobian of var_0). If var_0 is a ndarray, its Jacobian is set to zero.

Return type

pp.ad.Ad_array

sign(var)[source]
sin(var)[source]
sinh(var)[source]
tan(var)[source]
tanh(var)[source]
time_increment(op)[source]

Find the time increment of an operator tree.

Parameters

op (Operator) – Operator tree for which we need the time increment.

Returns

The difference of the operator tree now and at the previous time step.

Return type

Operator

Submodules