discretize.TreeMesh.cell_gradient_weak_form_robin#

TreeMesh.cell_gradient_weak_form_robin(alpha=0.0, beta=1.0, gamma=0.0)[source]#

Create Robin conditions pieces for weak form of the cell gradient operator (cell centers to faces).

This method returns the pieces required to impose Robin boundary conditions for the discrete weak form gradient operator that maps from cell centers to faces. These pieces are needed when constructing the discrete representation of the inner product \(\langle \vec{u} , \nabla \phi \rangle\) according to the finite volume method.

Where general boundary conditions are defined on \(\phi\) according to the Robin condition:

\[\alpha \phi + \beta \frac{\partial \phi}{\partial n} = \gamma\]

the user supplies values for \(\alpha\), \(\beta\) and \(\gamma\) for all boundary faces. cell_gradient_weak_form_robin returns the matrix \(\mathbf{B}\) and vector \(\mathbf{b}\) required for the discrete representation of \(\langle \vec{u} , \nabla \phi \rangle\). See the notes section for a comprehensive description.

Parameters:

alphascalar or (n_boundary_faces) array_like: Defines \(\alpha\) for Robin boundary condition. Can be defined as a scalar or array_like. If array_like, the length of the array must be equal to the number of boundary faces.
betascalar or (n_boundary_faces) array_like: Defines \(\beta\) for Robin boundary condition. Can be defined as a scalar or array_like. If array_like, must have the same length as alpha .
gamma: scalar or (n_boundary_faces) array_like or (n_boundary_faces, n_rhs) array_like: Defines \(\gamma\) for Robin boundary condition. If array like, gamma can have shape (n_boundary_face,). Can also have shape (n_boundary_faces, n_rhs) if multiple systems have the same alpha and beta parameters.

Returns:

B(n_faces, n_cells) scipy.sparse.csr_matrix: A sparse matrix dependent on the values of alpha, beta and gamma supplied
b(n_faces) numpy.ndarray or (n_faces, n_rhs) numpy.ndarray: A vector dependent on the values of alpha, beta and gamma supplied

Notes

For the gradient of a scalar \(\phi\), the weak form is implemented by taking the inner product with a piecewise-constant test function \(\vec{u}\) and integrating over the domain:

\[\langle \vec{u} , \nabla \phi \rangle \; = \int_\Omega \vec{u} \cdot (\nabla \phi) \, dv\]

For a discrete representation of \(\phi\) at cell centers, the gradient operator maps from cell centers to faces. To implement the boundary conditions in this case, we must use integration by parts and re-express the inner product as:

\[\langle \vec{u} , \nabla \phi \rangle \; = - \int_V \phi \, (\nabla \cdot \vec{u} ) \, dV + \oint_{\partial \Omega} \phi \hat{n} \cdot \vec{u} \, da\]

where the Robin condition is applied to \(\phi\) on the boundary. The discrete approximation to the above expression is given by:

\[\langle \vec{u} , \nabla \phi \rangle \; \approx - \boldsymbol{u^T \big ( D^T M_c - B \big ) \phi + u^T b}\]

where \(\mathbf{D}\) is the face_divergence and \(\mathbf{M_c}\) is the cell center inner product matrix (just a diagonal matrix comprised of the cell volumes). cell_gradient_weak_form_robin returns the matrix \(\mathbf{B}\) and vector \(\mathbf{b}\) based on the parameters alpha , beta and gamma provided.

Examples

Here we form all of pieces required to construct the discrete representation of the inner product between \(\mathbf{u}\) for specified Robin boundary conditions. We define \(\boldsymbol{\phi}\) at cell centers and \(\mathbf{u}\) on the faces. We begin by creating a small 2D tensor mesh:

>>> from discretize import TensorMesh
>>> import numpy as np
>>> import scipy.sparse as sp

>>> h = np.ones(32)
>>> mesh = TensorMesh([h, h])

We then define alpha, beta, and gamma parameters for a zero Dirichlet condition on the boundary faces. This corresponds to setting:

>>> alpha = 1.0
>>> beta = 0.0
>>> gamma = 0.0

Next, we construct all of the necessary pieces required to take the discrete inner product:

>>> B, b = mesh.cell_gradient_weak_form_robin(alpha, beta, gamma)
>>> Mc = sp.diags(mesh.cell_volumes)
>>> Df = mesh.face_divergence

In practice, these pieces are usually re-arranged when used to solve PDEs with the finite volume method. However, if you wanted to create a function which computes the discrete inner product for any \(\mathbf{u}\) and \(\boldsymbol{\phi}\):

>>> def inner_product(u, phi):
...     return u @ (-Df.T @ Mc + B) @ phi + u @ b