discretize.TreeMesh.edge_divergence_weak_form_robin#

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

Create Robin conditions pieces for weak form of the edge divergence operator (edges to nodes).

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

To implement the boundary conditions, we assume

\[\vec{u} = \nabla \phi\]

for some scalar function \(\phi\). Boundary conditions are imposed on the scalar function 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 nodes or faces. For the values supplied, edge_divergence_weak_form_robin returns the matrix \(\mathbf{B}\) and vector \(\mathbf{b}\) required for the discrete representation of \(\langle \psi , \nabla \cdot \vec{u} \rangle\). See the notes section for a comprehensive description.

Parameters:

alphascalar or 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 or boundary nodes.
betascalar or 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. Cannot be zero.
gammascalar or array_like: Defines \(\gamma\) for Robin boundary condition. If array like, gamma can have shape (n_boundary_xxx,). Can also have shape (n_boundary_xxx, n_rhs) if multiple systems have the same alpha and beta parameters.

Returns:

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

Notes

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

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

For a discrete representation of the vector \(\vec{u}\) that lives on mesh edges, the divergence operator must map from edges to nodes. To implement boundary conditions in this case, we must use the integration by parts to re-express the inner product as:

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

Assuming \(\vec{u} = \nabla \phi\), the above equation becomes:

\[\langle \psi , \nabla \cdot \vec{u} \rangle \, = - \int_V \nabla \phi \cdot \nabla \psi \, dV + \oint_{\partial \Omega} \psi \, \frac{\partial \phi}{\partial n} \, da\]

Substituting in the Robin conditions:

\[\langle \psi , \nabla \cdot \vec{u} \rangle \, = - \int_V \nabla \phi \cdot \nabla \psi \, dV + \oint_{\partial \Omega} \psi \, \frac{(\gamma - \alpha \Phi)}{\beta} \, da\]

therefore, beta cannot be zero.

The discrete approximation to the above expression is given by:

\[\langle \psi , \nabla \cdot \vec{u} \rangle \, \approx - \boldsymbol{\psi^T \big ( G_n^T M_e G_n - B \big ) \phi + \psi^T b}\]

where

\[\boldsymbol{u} = \boldsymbol{G_n \, \phi}\]

\(\mathbf{G_n}\) is the nodal_gradient and \(\mathbf{M_e}\) is the edge inner product matrix (see get_edge_inner_product). edge_divergence_weak_form_robin returns the matrix \(\mathbf{B}\) and vector \(\mathbf{b}\) based on the parameters alpha , beta and gamma provided.

Examples

Here we construct all of the pieces required for the discrete representation of \(\langle \psi , \nabla \cdot \vec{u} \rangle\) for specified Robin boundary conditions. We define \(\mathbf{u}\) on the edges, and \(\boldsymbol{\psi}\) and \(\boldsymbol{\psi}\) on the nodes. 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 Neumann condition on the boundary faces. This corresponds to setting:

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

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

>>> B, b = mesh.edge_divergence_weak_form_robin(alpha, beta, gamma)
>>> Me = mesh.get_edge_inner_product()
>>> Gn = mesh.nodal_gradient

In practice, these pieces are usually re-arranged when used to solve PDEs with the finite volume method. Because the boundary conditions are applied to the scalar potential \(\phi\), we create a function which computes the discrete inner product for any \(\boldsymbol{\psi}\) and \(\boldsymbol{\phi}\) where \(\mathbf{u} = \boldsymbol{G \, \phi}\):

>>> def inner_product(psi, phi):
...     return psi @ (-Gn.T @ Me @ Gn + B) @ phi + psi @ b