"""Module for base finite volume discretization classes."""
import abc
import warnings
from typing import Type
import numpy as np
import porepy as pp
import scipy.sparse as sps
import pygeon as pg
[docs]
class FiniteVolumeDiscretization(abc.ABC):
"""
Abstract class for PyGeoN finite volume discretization methods.
"""
[docs]
def __init__(self, keyword=pg.UNITARY_DATA) -> None:
"""
Initialize the FiniteVolumeDiscretization object.
Args:
keyword (str): The keyword used to identify the discretization method.
Default is pg.UNITARY_DATA.
Returns:
None
"""
self.keyword = keyword
self.bc_type: Type[pg.FiniteVolumeBC]
[docs]
def ndof(self, sd) -> int:
"""
Returns the number of degrees of freedom on a given grid.
Args:
sd (pg.Grid): The grid object.
Returns:
int: The number of degrees of freedom.
"""
return self.ndof_per_cell(sd) * sd.num_cells
[docs]
def assemble_system_matrix(
self, sd: pg.Grid, data: dict | None = None
) -> sps.csc_array:
"""
Assemble the system matrix, using the material parameters in the data dictionary
Args:
sd (pg.Grid): Grid, or a subclass.
data (dict): The data dictionary.
Returns:
sps.csc_array: The discretization matrix.
"""
# Assemble the zero'th order accumulation terms
M = self.assemble_accumulation_terms(sd, data)
# Assemble the second order terms, given by div(dual variables)
A = self.div(sd) @ self.assemble_dual_var_map(sd, data)
# Assemble the matrix
return (M + A).tocsc()
[docs]
def div(self, sd) -> sps.csc_array:
"""
Assembles the block-diagonal divergence operator acting on all the face-based
dual variables.
Args:
sd (pg.Grid): The grid object.
Returns:
sps.csc_array: The divergence operator
"""
return sps.kron(sps.eye_array(self.ndof_per_cell(sd)), pg.div(sd), format="csc")
[docs]
def face_area_scaling(self, sd) -> np.ndarray:
"""
Assembles the scaling vector with the face areas, for the face-based dual
variables.
Args:
sd (pg.Grid): The grid object.
Returns:
np.ndarray: The scaling vector
"""
return np.tile(sd.face_areas, self.ndof_per_cell(sd))
[docs]
def check_nonnegative_weights(self, weight: np.ndarray) -> None:
"""
Check whether all weighted distances are nonnegative.
Args:
weight (np.ndarray): The physical parameter weights.
"""
if np.any(weight < 0):
warnings.warn(
f"There are {(weight < 0).sum()} extra-cellular \
cell centers."
)
[docs]
def compute_harmonic_avg(self, faces: np.ndarray, dists: np.ndarray) -> np.ndarray:
"""
Compute $(1 / \delta_i + 1 / \delta_j)^{-1}$ at each face, between cells i and
j. This is used to compute the effective permeability at the face.
Args:
faces (np.ndarray): The extended array of faces.
dists (np.ndarray): The extended array of weighted distances.
Returns:
np.ndarray: The face-wise harmonic averages.
"""
return np.array(1 / np.bincount(faces, weights=dists))
[docs]
def get_bcs_from_data(self, sd: pg.Grid, data: dict | None) -> pg.FiniteVolumeBC:
"""
Extracts the FiniteVolumeBC object from the data dictionary, if it exists.
Else, it creates a new one and inserts it in data.
Args:
sd (pg.Grid): The grid object.
data (dict): The data dictionary
Returns:
pg.FiniteVolumeBC: The boundary condition object
"""
if not data:
data = pp.initialize_data({}, self.keyword)
if "bc" in data[pp.PARAMETERS][self.keyword]:
bcs = data[pp.PARAMETERS][self.keyword]["bc"]
else:
bcs = self.bc_type(sd, data, self.keyword)
return bcs
[docs]
def assemble_rhs_boundary_vector(
self, sd: pg.Grid, data: dict | None = None
) -> np.ndarray:
"""
Assembles the right-hand side vector related to the boundary conditions.
Args:
sd (pg.Grid): The grid object.
data (dict): The data dictionary
Returns:
np.ndarray: The right-hand side vector
"""
A_rhs = self.assemble_bdry_dual_var_map(sd, data)
bcs = self.get_bcs_from_data(sd, data)
dual = bcs.dual_var / sd.face_areas
prim = bcs.primary_var
g = np.hstack((dual.ravel(), prim.ravel()))
return -self.div(sd) @ A_rhs @ g
[docs]
@abc.abstractmethod
def ndof_per_cell(self, sd: pg.Grid) -> int:
"""
Returns the number of degrees of freedom per cell.
Args:
sd (pg.Grid): The grid object.
Returns:
int: The number of degrees of freedom per cell.
"""
[docs]
@abc.abstractmethod
def assemble_dual_var_map(self, sd: pg.Grid, data: dict | None) -> sps.csc_array:
"""
Assemble the mapping from cell-based primary variables to face-based dual
variables.
Args:
sd (pg.Grid): Grid, or a subclass.
data (dict): The data dictionary
Returns:
sps.csc_array: The matrix mapping primary to dual variables
"""
[docs]
@abc.abstractmethod
def assemble_accumulation_terms(
self, sd: pg.Grid, data: dict | None
) -> sps.csc_array:
"""
Assemble the zeroth-order terms for the primary variables.
Args:
sd (pg.Grid): Grid, or a subclass.
data (dict): The data dictionary
Returns:
sps.csc_array: The diagonal mass matrix
"""
[docs]
@abc.abstractmethod
def assemble_bdry_dual_var_map(
self, sd: pg.Grid, data: dict | None
) -> sps.csc_array:
"""
Assembles the matrix that maps from the boundary condition values to the dual
variables on the boundary faces.
Args:
sd (pg.Grid): Grid, or a subclass.
data (dict): The data dictionary
Returns:
sps.csc_array: The matrix to be multiplied with the boundary data g.
"""