diff --git a/src/porepy/numerics/ad/functions.py b/src/porepy/numerics/ad/functions.py index 45366cb474..f8c149906e 100644 --- a/src/porepy/numerics/ad/functions.py +++ b/src/porepy/numerics/ad/functions.py @@ -1,22 +1,26 @@ -""" -This module contains functions to be wrapped in a pp.ad.Function and used as part -of compound pp.ad.Operators, i.e. as (terms of) equations. +"""This module contains functions to be wrapped in a +:class:`~porepy.numerics.ad.operator_functions.Function` and used as part +of compound :class:`~porepy.numerics.ad.operators.Operator`, i.e. as (terms of) equations. + +Some functions depend on non-ad objects. This requires that the function ``f`` be wrapped +in an ``ad.Function`` using partial evaluation: -Some functions depend on non-ad objects. This requires that the function (f) be wrapped -in an ad Function using partial evaluation: +Examples: + >>> from functools import partial + >>> AdFunction = pp.ad.Function(partial(f, other_parameter), "name") + >>> equation: pp.ad.Operator = AdFunction(var) - 2 * var - from functools import partial - AdFunction = pp.ad.Function(partial(f, other_parameter), "name") - equation: pp.ad.Operator = AdFunction(var) - 2 * var + with ``var`` being some AD variable. -with var being some ad variable. + Note that while the argument to ``AdFunction`` is a + :class:`~porepy.numerics.ad.operators.Operator, the wrapping in + ``pp.ad.Function`` implies that upon parsing, + the argument passed to ``f`` will be an Ad_array. -Note that while the argument to AdFunction is a pp.ad.Operator, the wrapping in -pp.ad.Function implies that upon parsing, the argument passed to f will be an Ad_array. """ from __future__ import annotations -from typing import Callable, Union +from typing import Callable import numpy as np import scipy.sparse as sps @@ -90,25 +94,18 @@ def l2_norm(dim: int, var: pp.ad.Ad_array) -> pp.ad.Ad_array: """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] + 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 - Usage note: + 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 : pp.ad.Ad_array - Ad operator (variable or expression) which is argument of the norm - function. - - Returns - ------- - pp.ad.Ad_array + Parameters: + dim: Dimension, i.e. number of vector components. + var: Ad operator (variable or expression) which is argument of the norm function. + + Returns: The norm of var with appropriate val and jac attributes. """ @@ -265,30 +262,24 @@ def heaviside(var, zerovalue: float = 0.5): def heaviside_smooth(var, eps: float = 1e-3): - """ - Smooth (regularized) version of the Heaviside function. - - Parameters - ---------- - var : Ad_array or ndarray - Input array. - eps : float, optional - Regularization parameter. The default is 1E-3. The function will - convergence to the Heaviside function in the limit when eps --> 0 - - Returns - ------- - Ad_array or ndarray (depending on the input) - Regularized heaviside function (and its Jacobian if applicable). - - 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 + """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). """ if isinstance(var, Ad_array): @@ -315,26 +306,19 @@ def __call__(self, var, zerovalue: float = 0.5): return np.heaviside(var) # type: ignore -def maximum( - var0: pp.ad.Ad_array, var1: Union[pp.ad.Ad_array, np.ndarray] -) -> pp.ad.Ad_array: +def maximum(var0: pp.ad.Ad_array, var1: pp.ad.Ad_array | np.ndarray) -> pp.ad.Ad_array: """Ad maximum function represented as an Ad_array. The second argument is allowed to be constant, with a numpy array originally wrapped in a pp.ad.Array, whereas the first argument is expected to be an Ad_array originating from a pp.ad.Operator. - Parameters - ---------- - var0 : pp.ad.Ad_array - Ad operator (variable or expression). - var1 : Union[pp.ad.Ad_array, pp.ad.Array] - Ad operator (variable or expression) OR ad Array. + Parameters: + var0: First argument to the maximum function. + var1: Second argument. - Returns - ------- - pp.ad.Ad_array - The maximum of var0 and var1 with appropriate val and jac attributes. + Returns: + The maximum of ``var0`` and ``var1`` with appropriate val and jac attributes. """ vals = [var0.val.copy()] @@ -357,22 +341,17 @@ def maximum( def characteristic_function(tol: float, var: pp.ad.Ad_array): """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. + Returns 1 if ``var.val`` is within absolute tolerance = ``tol`` of zero. + The derivative is set to zero independent of ``var.val``. - Usage note: - See module level documentation on how to wrap functions like this in ad.Function. + 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 : pp.ad.Ad_array - Ad operator (variable or expression). + Parameters: + tol: Absolute tolerance for comparison with 0 using np.isclose. + var: Ad operator (variable or expression). - Returns - ------- - pp.ad.Ad_array + Returns: The characteristic function of var with appropriate val and jac attributes. """ diff --git a/src/porepy/numerics/ad/operator_functions.py b/src/porepy/numerics/ad/operator_functions.py index 4442e5981b..3c38ca5b2d 100644 --- a/src/porepy/numerics/ad/operator_functions.py +++ b/src/porepy/numerics/ad/operator_functions.py @@ -1,14 +1,14 @@ -""" -Contains callable operators representing functions to be called with other +"""This module contains callable operators representing functions to be called with other operators as input arguments. Contains also a decorator class for callables, which transforms them automatically in the specified operator function type. + """ from __future__ import annotations import abc from functools import partial -from typing import Callable, List, Optional, Type, Union +from typing import Callable, Optional, Type, Union import numpy as np import scipy.sparse as sps @@ -38,7 +38,7 @@ class AbstractFunction(Operator): The abstraction intends to provide means for approximating operators, where values are e.g. interpolated and the Jacobian is approximated using FD. - IMPLEMENTATION NOTE: + Note: One can flag the operator as ``ad_compatible``. If flagged, the AD framework passes AD arrays directly to the callable ``func`` and will **not** call the abstract methods for values and the Jacobian during operator parsing. @@ -47,6 +47,25 @@ class AbstractFunction(Operator): For now only one child class, porepy.ad.Function, flags itself always as AD compatible. + Parameters: + func: callable Python object representing a (numeric) function. + Expected to take numerical information in some form and return numerical + information in the same form. + name: name of this instance as an AD operator + array_compatible (optional): If true, the callable ``func`` will be called + using arrays (numpy.typing.ArrayLike). Flagging this true, the user ensures + that the callable can work with arrays and return respectively + formatted output. If false, the function will be evaluated element-wise + (scalar input). Defaults to False. + ad_compatible (Optional): If true, the callable ``func`` will be called using + the porepy.ad.Ad_array. + + Note that as of now, this will effectively bypass the abstract methods + for generating values and the Jacobian, assuming both will be provided + correctly by the return value of ``func``. + + Defaults to False. + """ def __init__( @@ -56,26 +75,6 @@ def __init__( array_compatible: bool = False, ad_compatible: bool = False, ): - """Configures this instance as a valid AD operator. - The passed callable is expected to take numerical information in some form and - return numerical information in the same form form. - - Args: - func (Callable): callable Python object representing a (numeric) function - name (str): name of this instance as an AD operator - array_compatible (default=False): If true, the callable ``func`` will be called - using arrays (numpy.typing.ArrayLike). Flagging this true, the user ensures - that the callable can work with arrays and return respectively - formatted output. If false, the function will be evaluated element-wise - (scalar input) - ad_compatible (default=False): If true, the callable ``func`` will be called using - the porepy.ad.Ad_array. - - Important NOTE: As of now, this will effectively bypass the abstract methods - for generating values and the Jacobian, assuming both will be provided - correctly by the return value of ``func`` - - """ ### PUBLIC # Reference to callable passed at instantiation self.func: Callable = func @@ -89,16 +88,16 @@ def __init__( self._operation: Operation = Operation.evaluate self._set_tree() - def __call__(self, *args: Union[pp.ad.Operator, Ad_array]): + def __call__(self, *args: pp.ad.Operator | Ad_array) -> pp.ad.Operator: """Renders this function operator callable, fulfilling its notion as 'function'. - Args: + Parameters: *args: AD operators passed as symbolic arguments for the callable passed at instantiation. Returns: - porepy.ad.Operator: Operator with call-arguments as children in the operator tree. - The assigned operation is ``evaluate``. + Operator with call-arguments as children in the operator tree. + The assigned operation is ``evaluate``. """ children = [self, *args] @@ -142,18 +141,18 @@ def parse(self, md: pp.MixedDimensionalGrid): The real work will be done by combining the function with arguments, during parsing of an operator tree. - Args: + Parameters: md: Mixed-dimensional grid. Returns: - The instance itself + The instance itself. + """ return self @abc.abstractmethod def get_values(self, *args: Ad_array) -> np.ndarray: - """ - Abstract method for evaluating the callable passed at instantiation. + """Abstract method for evaluating the callable passed at instantiation. This method will be called during the operator parsing. The AD arrays passed as arguments will be in the same order as the operators passed to @@ -162,12 +161,13 @@ def get_values(self, *args: Ad_array) -> np.ndarray: The returned numpy array will be set as 'val' argument for the AD array representing this instance. - Args: + Parameters: *args: Ad_array representation of the operators passed during the call to this instance Returns: - numpy.ndarray: numeric function values + Function values in numerical format. + """ pass @@ -183,28 +183,35 @@ def get_jacobian(self, *args: Ad_array) -> sps.spmatrix: The returned numpy array will be be set as 'jac' argument for the AD array representing this instance. - NOTE: the necessary dimensions for the jacobian can be extracted from the dimensions - of the Jacobians of passed Ad_array instances. + Note: + The necessary dimensions for the jacobian can be extracted from the dimensions + of the Jacobians of passed Ad_array instances. - Args: + Parameters: *args: Ad_array representation of the operators passed during the call to this instance Returns: - numpy.ndarray: numeric representation of the Jacobian of this function + Numeric representation of the Jacobian of this function. + """ pass class AbstractJacobianFunction(AbstractFunction): - """'Half'-abstract base class, providing a call to the callable ``func`` in order to obtain - numeric function values. + """Partially abstract base class, providing a call to the callable ``func`` in order to + obtain numeric function values. What remains abstract is the Jacobian. """ def get_values(self, *args: Ad_array) -> np.ndarray: + """ + Returns: + The direct evaluation of the callable using ``val`` of passed AD arrays. + + """ # get values of argument Ad_arrays. vals = (arg.val for arg in args) @@ -231,14 +238,15 @@ class Function(AbstractFunction): 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. + 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 (see ``ad_compatible``) + 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 (see ``ad_compatible``) """ @@ -259,14 +267,13 @@ class ConstantFunction(AbstractFunction): zero Jacobian. It still has to be called though since it fulfills the notion of a 'function'. + + Parameters: + values: constant values per cell. + """ def __init__(self, name: str, values: np.ndarray): - """ - Args: - values: constant values per cell - - """ # dummy function, takes whatever and returns only the pre-set values def func(*args): return values @@ -284,13 +291,13 @@ def get_values(self, *args: Ad_array) -> np.ndarray: def get_jacobian(self, *args: Ad_array) -> sps.spmatrix: """ - IMPLEMENTATION 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. + 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. @@ -300,25 +307,23 @@ def get_jacobian(self, *args: Ad_array) -> sps.spmatrix: class DiagonalJacobianFunction(AbstractJacobianFunction): - """ - Approximates the Jacobian of the function using identities and scalar - multipliers per dependency. + """Approximates the Jacobian of the function using identities and scalar multipliers + per dependency. + + Parameters: + multipliers: 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. + """ def __init__( self, func: Callable, name: str, - multipliers: Union[List[float], float], + multipliers: float | list[float], array_compatible: bool = False, ): - """ - Args: - multipliers: 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. - - """ super().__init__(func, name, array_compatible) # check and format input for further use if isinstance(multipliers, list): @@ -351,6 +356,22 @@ class InterpolatedFunction(AbstractFunction): 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: lower bounds for the domain of ``func``. + max_val: upper bound for the domain. + npt: number of interpolation points per dimension of the domain. + order: 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. + """ def __init__( @@ -364,20 +385,6 @@ def __init__( preval: bool = False, array_compatible: bool = False, ): - """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. - - Args: - min_val: lower bounds for the domain of ``func``. - max_val: upper bound for the domain. - npt: number of interpolation points per dimension of the domain. - order: Order of interpolation. Supports currently only linear order. - preval (default=False): 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. - """ super().__init__(func, name, array_compatible) ### PUBLIC @@ -432,59 +439,53 @@ def get_jacobian(self, *args: Ad_array) -> sps.spmatrix: class ADmethod: - """ - Automatic-Differentiation method/function. - - (Decorator) Class for methods representing e.g., physical properties. + """(Decorator) Class for methods representing e.g., physical properties. The decorated function is expected to take scalars/vectors and return a scalar/vector. - See example usage below. The return value will be an AD operator of a type passed to the decorator. - EXAMPLE USAGE: - .. code-block:: python - import porepy as pp + Examples: + .. code:: python + + import porepy as pp + + # decorating class methods + class IdealGas: - # 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 - @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 - # 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 + :class:`~porepy.numerics.ad.operators.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: decorated function object + ad_function_type: type reference to an AD operator function to be instantiated. + When instantiated, that instance will effectively replace ``func``. + operator_kwargs: keyword arguments to be passed when instantiating an operator + of type ``ad_function_type``. - With above code, the density of an instance of 'IdealGas' can be called using - :class:`~porepy.numerics.ad.operators.MergedVariable` representing - pressure and temperature. - Analogously, `dummy_rel_perm` can be called with one representing the saturation. """ def __init__( self, func: Optional[Callable] = None, - ad_function_type: Type["AbstractFunction"] = Function, + ad_function_type: Type[AbstractFunction] = Function, operator_kwargs: dict = {}, ) -> None: - """ - Decorator class constructor. - Saves information about the requested AD Function type and keyword arguments necessary - for its instantiation. - - NOTE: If used as decorator WITHOUT explicit instantiation, the instantiation will be - done implicitly with above default arguments (that's how Python decorators work). - - Args: - func: decorated function object - ad_function_type: type reference to an AD operator function to be instantiated. - When instantiated, that instance will effectively replace ``func``. - operator_kwargs: keyword arguments to be passed when instantiating an operator - of type ``ad_function_type``. - """ # reference to decorated function object self._func = func # mark if decoration without explicit call to constructor @@ -498,8 +499,7 @@ def __init__( self._op_kwargs = operator_kwargs def __call__(self, *args, **kwargs) -> Union[ADmethod, pp.ad.Operator]: - """ - Wrapper factory. + """Wrapper factory. The decorated object is wrapped and/or evaluated here. Dependent on whether the decorated function is a method belonging to a class, @@ -508,7 +508,7 @@ def __call__(self, *args, **kwargs) -> Union[ADmethod, pp.ad.Operator]: If bound to a class instance, the wrapper will include a partial function, where the instance of the class was already passed beforehand. - IMPLEMENTATION NOTE: + Note: If the decorator was explicitly instantiated during decoration, that instance will effectively be replaced by another decorator instance created here in the call. @@ -516,6 +516,7 @@ def __call__(self, *args, **kwargs) -> Union[ADmethod, pp.ad.Operator]: used as a decorator, hence properly dereferencing the original instance. If used differently or if another reference is saved between explicit instantiation and call, this is a potential memory leak. + """ # If decorated without explicit init, the function is passed during a call to # the decorator as first argument. @@ -572,11 +573,10 @@ def __get__(self, binding_instance: object, binding_type: type) -> Callable: are always passed in unbound form to the decorator when the code is evaluated (i.e. they don't have a reference to the `self` argument, contrary to bound - methods) + methods) - Args: - binding_instance: instance, whose method/attribute has been decorated - by this class. + Parameters: + binding_instance: instance, whose method has been decorated by this class. binding_type: type instance of the decorated method's class/owner. """ @@ -586,14 +586,15 @@ def __get__(self, binding_instance: object, binding_type: type) -> Callable: # This will trigger the function evaluation. return partial(self.__call__, binding_instance) - def ad_wrapper(self, *args, **kwargs) -> AbstractFunction: + def ad_wrapper(self, *args, **kwargs) -> Operator: """Actual wrapper function. Constructs the necessary AD-Operator class wrapping the decorated callable and performs the evaluation/call. - Args: + Parameters: *args: arguments for the call to the wrapping AD operator function **kwargs: keyword argument for the call to the wrapping Ad operator function + """ # Make sure proper assignment of callable was made assert self._func is not None diff --git a/src/porepy/numerics/ad/operators.py b/src/porepy/numerics/ad/operators.py index 8b3b044322..3f96b5fbd2 100644 --- a/src/porepy/numerics/ad/operators.py +++ b/src/porepy/numerics/ad/operators.py @@ -43,18 +43,25 @@ def _get_shape(mat): class Operator: - """Superclass for all Ad operators. + """Parent class for all AD operators. - Objects of this class is not meant to be initiated directly; rather the various + This class is 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 operations. - Attributes: - subdomains: List of subdomains (subdomains) on which the operator is defined. Will be - empty for operators not associated with specific subdomains. - interfaces: List of edges (tuple of subdomains) in the mixed-dimensional grid on which - the operator is defined. Will be empty for operators not associated with - specific subdomains. + Contains a tree structure of child operators for the recursive forward evaluation. + + Provides overload functions for basic arithmetic operations. + + Parameters: + name: 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). """ @@ -68,9 +75,23 @@ def __init__( if name is None: name = "" self._name = name + + self._set_tree(tree) + self.subdomains: list[pp.Grid] = [] if subdomains is None else subdomains + """List of subdomains on which the operator is defined, passed at instantiation. + + Will be empty for operators not associated with specific subdomains. + + """ + self.interfaces: list[pp.MortarGrid] = [] if interfaces is None else interfaces - self._set_tree(tree) + """List of mortar grids in the mixed-dimensional grid on which the operator is defined, + passed at instantiation. + + Will be empty for operators not associated with specific subdomains. + + """ def _set_tree(self, tree=None): if tree is None: @@ -230,11 +251,8 @@ def _identify_discretizations( return _ad_utils.uniquify_discretization_list(all_discr) def discretize(self, mdg: pp.MixedDimensionalGrid) -> None: - """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 + """Perform the ``discretize`` function of all child operators which are discretizations + using data from mdg. """ unique_discretizations: dict[ @@ -243,15 +261,23 @@ def discretize(self, mdg: pp.MixedDimensionalGrid) -> None: _ad_utils.discretize_from_list(unique_discretizations, mdg) def is_leaf(self) -> bool: - """Check if this operator is a leaf in the tree-representation of an object. + """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: - bool: True if the operator has no children. Note that this implies that the - method parse() is expected to be implemented. + True if the operator has no children. + """ return len(self.tree.children) == 0 def set_name(self, name: str) -> None: + """Reset this object's name originally passed at instantiation. + + Parameters: + name: the new name to be assigned. + + """ self._name = name def previous_timestep(self) -> pp.ad.Operator: @@ -264,7 +290,7 @@ def previous_timestep(self) -> pp.ad.Operator: Returns: A copy of self, with all time dependent operators evaluated at the previous - time step. + time step. """ # Create a copy of the operator tree evaluated at a previous time step. This is done @@ -335,11 +361,19 @@ def _traverse_tree(op: Operator) -> Operator: def parse(self, mdg: pp.MixedDimensionalGrid) -> Any: """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(). + of atomic operators; such operators should be evaluated by the method :meth:`evaluate`. + + Parameters: + mdg: Mixed-dimensional grid on which this operator is to be parsed. + + Returns: + A numerical format representing this operator;s values on given domain. + """ raise NotImplementedError("This type of operator cannot be parsed right away") @@ -650,7 +684,7 @@ def __str__(self) -> str: return self._name if self._name is not None else "" def viz(self): - """Give a visualization of the operator tree that has this operator at the top.""" + """Draws a visualization of the operator tree that has this operator as its root.""" G = nx.Graph() def parse_subgraph(node): @@ -709,16 +743,17 @@ def evaluate( """Evaluate the residual and Jacobian matrix for a given state. Parameters: - dof_manager (pp.DofManager): used to represent the problem. Will be used + dof_manager: Used to represent the problem. Will be used to parse the sub-operators that combine to form this operator. - state (np.ndarray, optional): State vector for which the residual and its + 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: - An Ad-array representation of the residual and Jacobian. Note that the Jacobian - matrix need not be invertible, or ever square; this depends on the operator. + 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. """ # Get the mixed-dimensional grid used for the dof-manager. @@ -870,27 +905,24 @@ def _parse_other(self, other): class Matrix(Operator): """Ad representation of a sparse matrix. - For dense matrices, use an Array instead. + For dense matrices, use :class:`Array` instead. This is a shallow wrapper around the real matrix; it is needed to combine the matrix with other types of Ad objects. - Attributes: - shape (Tuple of ints): Shape of the wrapped matrix. + Parameters: + mat: Sparse matrix to be wrapped as an AD operator. + name: Name of this operator """ def __init__(self, mat: sps.spmatrix, name: Optional[str] = None) -> None: - """Construct an Ad representation of a matrix. - - Parameters: - mat (sps.spmatrix): Sparse matrix to be represented. - - """ super().__init__(name=name) self._mat = mat self._set_tree() + self.shape = mat.shape + """Shape of the wrapped matrix.""" def __repr__(self) -> str: return f"Matrix with shape {self._mat.shape} and {self._mat.data.size} elements" @@ -901,32 +933,31 @@ def __str__(self) -> str: s += self._name return s - def parse(self, mdg) -> sps.spmatrix: - """Convert the Ad matrix into an actual matrix. - - Parameters: - mdg (pp.MixedDimensionalGrid): Mixed-dimensional grid. Not used, but it is - needed as input to be compatible with parse methods for other operators. + def parse(self, mdg: pp.MixedDimensionalGrid) -> sps.spmatrix: + """See :meth:`Operator.parse`. Returns: - sps.spmatrix: The wrapped matrix. + The wrapped matrix. """ return self._mat def transpose(self) -> "Matrix": + """Returns an AD operator representing the transposed matrix.""" return Matrix(self._mat.transpose()) class Array(Operator): - """Ad representation of a numpy array. + """AD representation of a constant numpy array. - For sparse matrices, use a Matrix instead. + For sparse matrices, use :class:`Matrix` instead. + For time-dependent arrays see :class:`TimeDependentArray`. This is a shallow wrapper around the real array; it is needed to combine the array - with other types of Ad objects. + with other types of AD operators. - See also TimeDependentArray. + Parameters: + values: Numpy array to be represented. """ @@ -934,7 +965,7 @@ def __init__(self, values: np.ndarray, name: Optional[str] = None) -> None: """Construct an Ad representation of a numpy array. Parameters: - values (np.ndarray): Numpy array to be represented. + values: Numpy array to be represented. """ super().__init__(name=name) @@ -951,14 +982,10 @@ def __str__(self) -> str: return s def parse(self, mdg: pp.MixedDimensionalGrid) -> np.ndarray: - """Convert the Ad Array into an actual array. - - Parameters: - mdg (pp.MixedDimensionalGrid): Mixed-dimensional grid. Not used, but it is - needed as input to be compatible with parse methods for other operators. + """See :meth:`Operator.parse`. Returns: - np.ndarray: The wrapped array. + The wrapped array. """ return self._values @@ -968,20 +995,27 @@ class TimeDependentArray(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. + 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 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. - Attributes: - prev_time (boolean): 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. + Parameters: + name: Name of the variable. Should correspond to items in data[pp.STATE]. + subdomains: Subdomains on which the array is defined. Defaults to None. + interfaces: Interfaces on which the array is defined. Defaults to None. + Exactly one of subdomains and interfaces must be non-empty. + previous_timestep: Flag indicating if the array should be evaluated at the + previous time step. + + Raises: + ValueError: If either none of, or both of, subdomains and interfaces are empty. """ @@ -992,23 +1026,6 @@ def __init__( interfaces: Optional[list[pp.MortarGrid]] = None, previous_timestep: bool = False, ): - """Initialize a TimeDependentArray. - - Args: - name: Name of the variable. Should correspond to items in data[pp.STATE]. - subdomains: Subdomains on which the array is defined. Defaults to None. - interfaces: Interfaces on which the array is defined. Defaults to None. - - Exactly one of subdomains and interfaces must be non-empty. - - previous_timestep: Flag indicating if the array should be evaluated at the - previous time step. - - Raises: - ValueError: If either none of, or both of, subdomains and interfaces are - empty. - - """ self._name: str = name @@ -1047,12 +1064,17 @@ def __init__( self._is_interface_array = True self.prev_time: bool = previous_timestep + """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. + + """ self._set_tree() def previous_timestep(self) -> TimeDependentArray: - """Return a representation of this variable on the previous time step. - + """ Returns: This array represented at the previous time step. @@ -1070,10 +1092,11 @@ def parse(self, mdg: pp.MixedDimensionalGrid) -> np.ndarray: """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]. + ``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]``. - Args: + Parameters: mdg: Mixed-dimensional grid. Returns: @@ -1113,7 +1136,7 @@ def __repr__(self) -> str: class Scalar(Operator): """Ad representation of a scalar. - This is a shallow wrapper around the real scalar; it may be useful to combine + 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 @@ -1124,12 +1147,6 @@ class Scalar(Operator): """ def __init__(self, value: float, name: Optional[str] = None) -> None: - """Construct an Ad representation of a float. - - Parameters: - value (float): Number to be represented - - """ super().__init__(name=name) self._value = value self._set_tree() @@ -1144,40 +1161,35 @@ def __str__(self) -> str: return s def parse(self, mdg: pp.MixedDimensionalGrid) -> float: - """Convert the Ad Scalar into an actual number. - - Parameters: - mdg (pp.MixedDimensionalGrid): Mixed-dimensional grid. Not used, but it is - needed as input to be compatible with parse methods for other operators. + """See :meth:`Operator.parse`. Returns: - float: The wrapped number. + The wrapped number. """ return self._value class Variable(Operator): - """Ad representation of a variable defined on a single Grid or MortarGrid. + """AD operator representing a variable defined on a single grid or mortar grid. - For combinations of variables on different subdomains, see MergedVariable. + For combinations of variables on different subdomains, see :class:`MergedVariable`. - Conversion of the variable into numerical value should be done with respect to the - state of an array; see the method evaluate(). Therefore, the variable does not - implement a parse() method. + A conversion of the variable into numerical value should be done with respect to the + state of an array; see :meth:`Operator.evaluate`. Therefore, the variable does not + implement the method :meth:`Operator.parse`. - Attributes: - id (int): Unique identifier of this variable. - prev_iter (boolean): Whether the variable represents the state at the - previous iteration. - prev_time (boolean): Whether the variable represents the state at the - previous time step. - subdomains: List with one item, giving the single grid on which the operator is - defined. - interfaces: List with one item, giving the single edge (tuple of subdomains) on - which the operator is defined. + 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. - It is assumed that exactly one of subdomains and interfaces is defined. + Parameters: + name: Variable name. + ndof: 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. """ @@ -1196,20 +1208,6 @@ def __init__( previous_timestep: bool = False, previous_iteration: bool = False, ): - """Initiate an Ad representation of a variable associated with a grid or edge. - - It is assumed that exactly one of subdomains and interfaces is defined. - Parameters: - name (str): Variable name. - ndof (dict): Number of dofs per grid element. - subdomains (optional list of pp.Grid ): List with length one containing a grid. - interfaces (optional list of pp.MortarGrid): List with length one containing - an interface. - num_cells (int): Number of cells in the grid. Only sued if the variable - is on an interface. - - """ - self._name: str = name self._cells: int = ndof.get("cells", 0) self._faces: int = ndof.get("faces", 0) @@ -1230,21 +1228,25 @@ def __init__( self._g = self.interfaces[0] self._is_edge_var = True - self.prev_time: bool = previous_timestep - self.prev_iter: bool = previous_iteration - # The number of cells in the grid. Will only be used if grid_like is a tuple # that is, if this is a mortar variable self._num_cells = num_cells - self.id = next(self._ids) self._set_tree() - def size(self) -> int: - """Get the number of dofs for this grid. + self.prev_time: bool = previous_timestep + """Indicates whether the variable represents the state at the previous time step.""" + + self.prev_iter: bool = previous_iteration + """Indicates whether the variable represents the state at the previous iteration.""" + + self.id = next(self._ids) + """Unique identifier of this variable.""" + def size(self) -> int: + """ Returns: - int: Number of dofs. + The number of dofs of this variable on its grid. """ if self._is_edge_var: @@ -1260,10 +1262,10 @@ def size(self) -> int: ) def previous_timestep(self) -> Variable: - """Return a representation of this variable on the previous time step. - + """ Returns: - Variable: A representation of this variable, with self.prev_time=True. + A representation of this variable at the previous time step, + with its ``prev_time`` attribute set to ``True``. """ ndof = {"cells": self._cells, "faces": self._faces, "nodes": self._nodes} @@ -1277,10 +1279,10 @@ def previous_timestep(self) -> Variable: ) def previous_iteration(self) -> Variable: - """Return a representation of this variable on the previous time iteration. - + """ Returns: - Variable: A representation of this variable, with self.prev_iter=True. + A representation of this variable on the previous time iteration, + with its ``prev_iter`` attribute set to ``True``. """ ndof = {"cells": self._cells, "faces": self._faces, "nodes": self._nodes} @@ -1308,35 +1310,23 @@ def __repr__(self) -> str: class MergedVariable(Variable): - """Ad representation of a collection of variables that individually live on separate - subdomains of interfaces, but which it is useful to treat jointly. + """AD operator representing 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 the method evaluate(). Therefore, the class does not implement - a parse() method. - - Attributes: - sub_vars (List of Variable): List of variable on different subdomains or interfaces. - id (int): Counter of all variables. Used to identify variables. Usage of this - term is not clear, it may change. - prev_iter (boolean): Whether the variable represents the state at the - previous iteration. - prev_time (boolean): Whether the variable represents the state at the - previous time step. + state of an array; see :meth:`Operator.evaluate`. Therefore, the MergedVariable does not + implement the method :meth:`Operator.parse`. - It is assumed that exactly one of subdomains and interfaces is defined. + Parameters: + variables: List of variables to be merged. Should all have the same name. """ def __init__(self, variables: list[Variable]) -> None: - """Create a merged representation of variables. - - Parameters: - variables (list of Variable): Variables to be merged. Should all have the - same name. - """ self.sub_vars = variables + """List of single-grid variables which are merged into this merged variable, + passed at instantiation.""" # Use counter from superclass to ensure unique Variable ids self.id = next(Variable._ids) @@ -1367,19 +1357,18 @@ def __init__(self, variables: list[Variable]) -> None: self.prev_iter: bool = False def size(self) -> int: - """Get total size of the merged variable. - + """ Returns: - int: Total size of this merged variable. + Total size of this merged variable (sum of sizes of respective sub-variables). """ return sum([v.size() for v in self.sub_vars]) def previous_timestep(self) -> MergedVariable: - """Return a representation of this merged variable on the previous time step. - + """ Returns: - Variable: A representation of this variable, with self.prev_time=True. + A representation of this merged variable on the previous time + iteration, with its ``prev_iter`` attribute set to ``True``. """ new_subs = [var.previous_timestep() for var in self.sub_vars] @@ -1388,10 +1377,10 @@ def previous_timestep(self) -> MergedVariable: return new_var def previous_iteration(self) -> MergedVariable: - """Return a representation of this merged variable on the previous iteration. - + """ Returns: - Variable: A representation of this variable, with self.prev_iter=True. + A representation of this merged variable on the previous + iteration, with its ``prev_iter`` attribute set to ``True`` """ new_subs = [var.previous_iteration() for var in self.sub_vars] @@ -1400,8 +1389,12 @@ def previous_iteration(self) -> MergedVariable: return new_var def copy(self) -> "MergedVariable": - # A shallow copy should be sufficient here; the attributes are not expected to - # change. + """ + Returns: + A shallow copy should be sufficient here; the attributes are not expected to + change. + + """ return copy.deepcopy(self) def __repr__(self) -> str: @@ -1435,9 +1428,17 @@ def __repr__(self) -> str: class Tree: """Simple implementation of a Tree class. Used to represent combinations of Ad operators. + + References: + https://stackoverflow.com/questions/2358045/how-can-i-implement-a-tree-in-python + + + Parameters: + operation: See :data:`Operation` + children: List of children, either as Ad arrays or other :class:`Operator`. + """ - # https://stackoverflow.com/questions/2358045/how-can-i-implement-a-tree-in-python def __init__( self, operation: Operation, @@ -1452,5 +1453,5 @@ def __init__( self.add_child(child) def add_child(self, node: Union[Operator, Ad_array]) -> None: - # assert isinstance(node, (Operator, "pp.ad.Operator")) + """Adds a child to this instance.""" self.children.append(node)