porepy.numerics.fv.fvutils module

Various FV specific utility functions.

Implementation note: This could perhaps have been implemneted as a superclass for MPxA discertizations, however, due to the somewhat intricate inheritance relation between these methods, the current structure with multiple auxiliary methods emerged.

class ExcludeBoundaries(subcell_topology, bound, nd)

Bases: object

Wrapper class to store mappings needed in the finite volume discretizations. The original use for this class was for exclusion of equations that are redundant due to the presence of boundary conditions, hence the name. The use of this class has increased to also include linear transformation that can be applied to the subfaces before the exclusion operator. This will typically be a rotation matrix, so that the boundary conditions can be specified in an arbitary coordinate system.

The systems being set up in mpfa (and mpsa) describe continuity of flux and potential (respectively stress and displacement) on all sub-faces. For the boundary faces, unless a Robin condition is specified, only one of the two should be included (e.g. for a Dirichlet boundary condition, there is no concept of continuity of flux/stresses). This class contains mappings to eliminate the necessary fields in order set the correct boundary conditions.

exclude_boundary(other, transform=False)

Mapping to exclude faces/component with any boundary condition from local systems.

Parameters

other (scipy.sparse matrix) – Matrix of local equations for continuity of flux and pressure.

Returns

scipy.sparse matrix, with rows corresponding to faces/components with

Neumann conditions eliminated.

exclude_dirichlet(other, transform=True)

Mapping to exclude faces/components with Dirichlet boundary conditions from local systems.

Parameters

other (scipy.sparse matrix) – Matrix of local equations for continuity of flux and pressure.

Returns

scipy.sparse matrix, with rows corresponding to faces/components with

Dirichlet conditions eliminated.

exclude_neumann(other, transform=True)

Mapping to exclude faces/components with Neumann boundary conditions from local systems.

Parameters

other (scipy.sparse matrix) – Matrix of local equations for continuity of flux and pressure.

Returns

scipy.sparse matrix, with rows corresponding to faces/components with

Neumann conditions eliminated.

exclude_neumann_dirichlet(other, transform=True)

Mapping to exclude faces/components with Neumann and Dirichlet boundary conditions from local systems.

Parameters

other (scipy.sparse matrix) – Matrix of local equations for continuity of flux and pressure.

Returns

scipy.sparse matrix, with rows corresponding to faces/components with

Neumann conditions eliminated.

exclude_neumann_robin(other, transform=True)

Mapping to exclude faces/components with Neumann and Robin boundary conditions from local systems.

Parameters

other (scipy.sparse matrix) – Matrix of local equations for continuity of flux and pressure.

Returns

scipy.sparse matrix, with rows corresponding to faces/components with

Neumann conditions eliminated.

exclude_robin_dirichlet(other, transform=True)

Mapping to exclude faces/components with Robin and Dirichlet boundary conditions from local systems.

Parameters

other (scipy.sparse matrix) – Matrix of local equations for continuity of flux and pressure.

Returns

scipy.sparse matrix, with rows corresponding to faces/components with

Neumann conditions eliminated.

keep_neumann(other, transform=True)

Mapping to exclude faces/components that is not on the Neumann boundary conditions from local systems.

Parameters

other (scipy.sparse matrix) – Matrix of local equations for continuity of flux and pressure.

Returns

scipy.sparse matrix, with rows corresponding to faces/components with

all but Neumann conditions eliminated.

keep_robin(other, transform=True)

Mapping to exclude faces/components that is not on the Robin boundary conditions from local systems.

Parameters

other (scipy.sparse matrix) – Matrix of local equations for continuity of flux and pressure.

Returns

scipy.sparse matrix, with rows corresponding to faces/components with

all but Robin conditions eliminated.

class SubcellTopology(sd)

Bases: object

Class to represent data of subcell topology (interaction regions) for mpsa/mpfa.

g - the grid

Subcell topology, all cells seen apart

nno - node numbers

fno - face numbers

cno - cell numbers

subfno - subface numbers. Has exactly two duplicates for internal faces

subhfno - subface numbers. No duplicates

num_cno - cno.max() + 1

num_subfno - subfno.max() + 1

Subcell topology, after cells sharing a face have been joined. Use

terminology _unique, since subfno is unique after

nno_unique - node numbers after pairing sub-faces

fno_unique - face numbers after pairing sub-faces

cno_unique - cell numbers after pairing sub-faces

subfno_unique - subface numbers  after pairing sub-faces

unique_subfno - index of those elements in subfno that are included

in subfno_unique, that is subfno_unique = subfno[unique_subfno], and similarly cno_unique = cno[subfno_unique] etc.

num_subfno_unique = subfno_unique.max() + 1

pair_over_subfaces(other)

Transfer quantities from a cell-face base (cells sharing a face have their own expressions) to a face-based. The quantities should live on sub-faces (not faces)

The normal vector is honored, so that the here and there side get different signs when paired up.

The method is intended used for combining forces, fluxes, displacements and pressures, as used in MPSA / MPFA.

Parameters

other – sps.matrix, size (self.subhfno.size x something)

Returns

sps.matrix, size (self.subfno_unique.size x something)

pair_over_subfaces_nd(other)

nd-version of pair_over_subfaces, see above.

append_dofs_of_discretization(sd, d, kw1, kw2, k_dof)

Appends rows to existing discretizations stored as ‘stress’ and ‘bound_stress’ in the data dictionary on the nodes. Only applies to the highest dimension (for now, at least). The newly added faces are found from ‘new_faces’ in the data dictionary. Assumes all new rows and columns should be appended, not inserted to the “interior” of the discretization matrices. sd - grid object d - corresponding data dictionary kw1 - keyword for the stored discretization in the data dictionary,

e.g. ‘flux’

kw2 - keyword for the stored boundary discretization in the data

dictionary, e.g. ‘bound_flux’

boundary_to_sub_boundary(bound, subcell_topology)

Convert a boundary condition defined for faces to a boundary condition defined by subfaces.

Args: bound (pp.BoundaryCondition/pp.BoundarConditionVectorial):

Boundary condition given for faces.

subcell_topology (pp.fvutils.SubcellTopology):

The subcell topology defining the finite volume subgrid.

Returns

An instance of the

BoundaryCondition/BoundaryConditionVectorial class, where all face values of bound has been copied to the subfaces as defined by subcell_topology.

Return type

pp.BoundarCondition/pp.BoundarConditionVectorial

cell_ind_for_partial_update(sd, cells=None, faces=None, nodes=None)

Obtain indices of cells and faces needed for a partial update of the discretization stencil.

Implementation note: This function should really be split into three parts, one for each of the modes (cell, face, node).

The subgrid can be specified in terms of cells, faces and nodes to be updated. The method will then define a sufficiently large subgrid to account for changes in the flux discretization. The idea is that cells are used to account for updates in material parameters (or geometry), faces when faces are split (fracture growth), while the parameter nodes is mainly aimed at a gradual build of the discretization of the entire grid (for memory conservation, see comments in mpfa.mpfa()). For more details, see the implementations and comments below.

Cautionary note: The option to combine cells, faces and nodes in one go has not been tested. Problems may arise for grid configurations where separate regions are close to touching. This is however speculation at the time of writing.

Parameters

• g (core.grids.grid) – grid to be discretized

• cells (np.array, int, optional) – Index of cells on which to base the subgrid computation. Defaults to None.

• faces (np.array, int, optional) – Index of faces on which to base the subgrid computation. Defaults to None.

• nodes (np.array, int, optional) – Index of faces on which to base the subgrid computation. Defaults to None.

• sd (Grid) –

Returns

Cell indexes of the subgrid. No guarantee that they form

a connected grid.

np.array, int: Indexes of faces to have their discretization updated.

Return type

np.array, int

cell_scalar_to_subcell_vector(nd, sub_cell_index, cell_index)

Create mapping from sub-cells to cells for vector problems. For example, discretization of grad_p-term in Biot, where p is a cell-center scalar

Parameters

• nd – dimension

• sub_cell_index – sub-cell indices

• cell_index – cell indices

Return type

scipy.sparse.csr_matrix (shape num_subcells * nd, num_cells)

cell_vector_to_subcell(nd, sub_cell_index, cell_index)

Create mapping from sub-cells to cells for scalar problems. For example, discretization of div_g-term in mpfa with gravity, where g is a cell-center vector (dim nd)

Parameters

• nd – dimension

• sub_cell_index – sub-cell indices

• cell_index – cell indices

Return type

scipy.sparse.csr_matrix (shape num_subcells * nd, num_cells * nd)

compute_darcy_flux(mdg, keyword='flow', keyword_store=None, d_name='darcy_flux', p_name='pressure', lam_name='mortar_solution', well_name='well_flux', data=None, from_iterate=False)

Computes darcy_flux over all faces in the entire mixed-dimensional grid given pressures for all nodes, provided as node properties.

Parameter: mgd: mixed-dimensional grid with the following data fields for all nodes/grids:

‘flux’: Internal discretization of fluxes. ‘bound_flux’: Discretization of boundary fluxes. ‘pressure’: Pressure values for each cell of the grid (overwritten by p_name). ‘bc_val’: Boundary condition values.

and the following edge property field for all connected grids:

‘coupling_flux’: Discretization of the coupling fluxes.

keyword (str): defaults to ‘flow’. The parameter keyword used to obtain the

data necessary to compute the fluxes.

keyword_store (str): defaults to keyword. The parameter keyword determining

where the data will be stored.

d_name (str): defaults to ‘darcy_flux’. The parameter name which the computed

darcy_flux will be stored by in the dictionary.

p_name (str): defaults to ‘pressure’. The keyword that the pressure

field is stored by in the dictionary.

lam_name (str): defaults to ‘mortar_solution’. The keyword that the mortar flux

field is stored by in the dictionary.

data (dictionary): defaults to None. If mdg is mono-dimensional grid the data

dictionary must be given. If mdg is a multi-dimensional grid, this variable has no effect.

Returns

the same mixed-dimensional grid with the added field ‘darcy_flux’ added to all subdomain data fields. Note that the fluxes between grids will be added only at the mdg edge, not at the node fields. The signs of the darcy_flux correspond to the directions of the normals, in the edge/coupling case those of the higher grid. For edges between grids of equal dimension, there is an implicit assumption that all normals point from the second to the first of the sorted grids (mdg.sorted_nodes_of_edge(e)).

Return type

pp.MixedDimensionalGrid

compute_dist_face_cell(sd, subcell_topology, eta, return_paired=True)

Compute vectors from cell centers continuity points on each sub-face.

The location of the continuity point is given by

x_cp = (1-eta) * x_facecenter + eta * x_vertex

On the boundary, eta is set to zero, thus the continuity point is at the face center

Parameters

• sd – Grid

• subcell_topology – Of class subcell topology in this module

• eta – [0,1), eta = 0 gives cont. pt. at face midpoint, eta = 1 means at the vertex. If eta is given as a scalar this value will be applied to all subfaces except the boundary faces, where eta=0 will be enforced. If the length of eta equals the number of subfaces, eta[i] will be used in the computation of the continuity point of the subface s_t.subfno_unique[i]. Note that eta=0 at the boundary is ONLY enforced for scalar eta.

Returns

sps.csr() matrix representation of vectors. Size sd.nf x (sd.nc * sd.nd)

Raises

ValueError if the size of eta is not 1 or subcell_topology.num_subfno_unique. –

determine_eta(sd)

Set default value for the location of continuity point eta in MPFA and MSPA.

The function is intended to give a best estimate of eta in cases where the user has not specified a value.

Parameters

sd (Grid) – Grid for discretization

Returns

double. 1/3 if the grid in known to consist of simplicies (it is one of

TriangleGrid, TetrahedralGrid, or their structured versions). 0 if not.

Return type

float

expand_indices_incr(ind, dim, increment)

expand_indices_nd(ind, nd, direction='F')

Expand indices from scalar to vector form.

Examples: >>> i = np.array([0, 1, 3]) >>> expand_indices_nd(i, 2) (array([0, 1, 2, 3, 6, 7]))

>>> expand_indices_nd(i, 3, "C")
(array([0, 3, 9, 1, 4, 10, 2, 5, 11])

Parameters

• ind (ndarray) –

• nd (int) –

• direction –

Return type

ndarray

Returns

find_active_indices(parameter_dictionary, sd)

Process information in parameter dictionary on whether the discretization should consider a subgrid. Look for fields in the parameter dictionary called specified_cells, specified_faces or specified_nodes. These are then processed by the function pp.fvutils.cell-ind_for_partial_update.

If no relevant information is found, the active indices are all cells and faces in the grid.

Parameters

• parameter_dictionary (dict) – Parameters, potentially containing fields “specified_cells”, “specified_faces”, “specified_nodes”.

• sd (pp.Grid) – Grid to be discretized.

Returns

Cells to be included in the active grid. np.ndarary: Faces to have their discretization updated. NOTE: This may not

be all faces in the grid.

Return type

np.ndarray

map_hf_2_f(fno=None, subfno=None, nd=None, sd=None)

Create mapping from half-faces to faces for vector problems. Either fno, subfno and nd should be given or g (and optinally nd) should be given.

Args: EITHER:

fno (np.ndarray): face numbering in sub-cell topology based on unique subfno subfno (np.ndarrary): sub-face numbering nd (int): dimension

OR:

g (pp.Grid): If a grid is supplied the function will set:

fno = pp.fvutils.SubcellTopology(g).fno_unique subfno = pp.fvutils.SubcellTopology(g).subfno_unique

nd (int): Optinal, defaults to sd.dim. Defines the dimension of the vector.

Returns

map_subgrid_to_grid(sd, loc_faces, loc_cells, is_vector, nd=None)

Obtain mappings from the cells and faces of a subgrid back to a larger grid.

Parameters

• g (pp.Grid) – The larger grid.

• loc_faces (np.ndarray) – For each face in the subgrid, the index of the corresponding face in the larger grid.

• loc_cells (np.ndarray) – For each cell in the subgrid, the index of the corresponding cell in the larger grid.

• is_vector (bool) – If True, the returned mappings are sized to fit with vector variables, with nd elements per cell and face.

• nd (int, optional) – Dimension. Defaults to sd.dim.

• sd (Grid) –

Return type

tuple[numpy.ndarray, numpy.ndarray]

Retuns:

sps.csr_matrix, size (sd.num_faces, loc_faces.size): Mapping from local to

global faces. If is_vector is True, the size is multiplied with sd.dim.

sps.csr_matrix, size (loc_cells.size, sd.num_cells): Mapping from global to

local cells. If is_vector is True, the size is multiplied with sd.dim.

partial_discretization(sd, data, tensor, bnd, apertures, partial_discr, physics='flow')

Perform a partial (local) multi-point discretization on a grid with provided data, tensor and boundary conditions by

1) Appending the existing discretization matrices to the right size according to the added cells and faces. 2) Discretizing on the relevant subgrids by calls to the provided partial discretization function (mpfa_partial or mpsa_partial). 3) Zeroing out the rows corresponding to the updated faces. 4) Inserting the newly computed values to the just deleted rows.

partial_update_discretization(sd, data, keyword, discretize, dim=None, scalar_cell_right=None, vector_cell_right=None, scalar_face_right=None, vector_face_right=None, scalar_cell_left=None, vector_cell_left=None, scalar_face_left=None, vector_face_left=None, second_keyword=None)

Do partial update of discretization scheme.

This is intended as a helper function for the update_discretization methods of fv schemes. In particular, this method allows for a unified implementation of update in mpfa, mpsa and biot.

The implementation is somewhat heavy to cover both mpfa, mpsa and biot.

Parameters scalar_cell_right, vector_face_left etc. are lists of keys in the discretization matrix dictionary. They are used to tell whehter the matrix should be considered a cell or face quantity, and scalar of vector. Left and right are used to map rows and columns, respectively. Together these fields allows for mapping a discretization between grids.

If a term misses a right or left mapping, it will be ignored.

Parameters

• sd (Grid) –

• data (dict) –

• keyword (str) –

• discretize (Callable) –

• dim (Optional[int]) –

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

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

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

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

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

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

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

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

• second_keyword (Optional[str]) –

Return type

None

remove_nonlocal_contribution(raw_ind, nd, *args)

For a set of matrices, zero out rows associated with given faces, adjusting for the matrices being related to vector quantities if necessary.

Example: If raw_ind = np.array([2]), and nd = 2, rows 4 and 5 will be eliminated

(row 0 and 1 will in this case be associated with face 0, row 2 and 3 with face

1).

Parameters

• raw_ind (np.ndarray) – Face indices to have their rows eliminated.

• nd (int) – Spatial dimension. Needed to map face indices to rows.

• *args (sps.spmatrix) – Set of matrices. Will be eliminated in place.

Returns

DESCRIPTION.

Return type

None

scalar_divergence(sd)

Get divergence operator for a grid.

The operator is easily accessible from the grid itself, so we keep it here for completeness.

See also vector_divergence(g)

Parameters

sd (pp.Grid) – grid

Return type

csr_matrix

Returns

divergence operator

scalar_tensor_vector_prod(sd, k, subcell_topology)

Compute product of normal vectors and tensors on a sub-cell level. This is essentially defining Darcy’s law for each sub-face in terms of sub-cell gradients. Thus, we also implicitly define the global ordering of sub-cell gradient variables (via the interpretation of the columns in nk). NOTE: In the local numbering below, in particular in the variables i and j, it is tacitly assumed that sd.dim == sd.nodes.shape[0] == sd.face_normals.shape[0] etc. See implementation note in main method. :param sd: Discretization grid :type sd: pp.Grid :param k: The permeability tensor :type k: pp.Second_order_tensor :param subcell_topology: Wrapper class containing

subcell numbering.

Returns

sub-face wise product of normal vector and permeability tensor. cell_node_blocks pairings of node and cell indices, which together

define a sub-cell.

sub_cell_ind: index of all subcells

Return type

nk

Parameters

• sd (Grid) –

• k (SecondOrderTensor) –

• subcell_topology (fvutils.SubcellTopology) –

subproblems(sd, max_memory, peak_memory_estimate)

Parameters

• sd (Grid) –

• max_memory (int) –

• peak_memory_estimate (int) –

Return type

Generator[tuple[porepy.grids.grid.Grid, numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray], None, None]

vector_divergence(sd)

Get vector divergence operator for a grid g

It is assumed that the first column corresponds to the x-equation of face 0, second column is y-equation etc. (and so on in nd>2). The next column is then the x-equation for face 1. Correspondingly, the first row represents x-component in first cell etc.

Parameters

sd (pp.Grid) – grid

Return type

csr_matrix

Returns

vector_div (sparse csr matrix), dimensions: nd * (num_cells, num_faces)