4.2. Diffusion#

mirgecom.diffusion computes the diffusion operator.

Diffusion equation:

\[\frac{\partial u}{\partial t} = \nabla \cdot (\boldsymbol{\kappa} \nabla u)\]

where:

  • conserved variable \(u\)

  • scalar (isotropic) or diagonal-tensor (orthotropic) \(\boldsymbol{\kappa}\)

In the orthotropic case, the product \(\boldsymbol{\kappa} \nabla u\) is treated as a Hadamard product of arrays rather than a matrix-array multiplication. Fully anisotropic materials are not currently handled.

Due to the possibility of big differences in the magnitude of diffusivity coefficients, such as thermal conductivity in air and a solid, an harmonic average can be used to increase robustness and avoid numerical instabilities. In this case, to ensure flux continuity,

\[F = \kappa^- \nabla T^- = \kappa^+ \nabla T^+ = \bar{\kappa} \left( \frac{\nabla T^- + \nabla T^+}{2} \right)\]

where

\[\bar{\kappa} = \frac{2 \kappa^-_{ii} \kappa^+_{ii}} {\kappa^-_{ii} + \kappa^+_{ii}}\]

with \(\kappa_{ii}\) being either the individual components of the diffusivity array (orthotropic material) or a single scalar (isotropic).

4.2.1. Flux functions#

mirgecom.diffusion.diffusion_flux(kappa, grad_u)[source]#

Compute the diffusive flux \(-\kappa \nabla u\).

Parameters:
Returns:

The diffusive flux.

Return type:

meshmode.dof_array.DOFArray

mirgecom.diffusion.grad_facial_flux_central(kappa_tpair, u_tpair, normal)[source]#

Compute the numerical flux for \(\nabla u\).

Uses a simple average of the two sides’ values:

\[F = -\frac{u^- + u^+}{2} \hat{n}.\]
mirgecom.diffusion.grad_facial_flux_weighted(kappa_tpair, u_tpair, normal)[source]#

Compute the numerical flux for \(\nabla u\) using a weighted average.

Weights each side’s value by the corresponding thermal conductivity \(\kappa\):

\[F = -\frac{\kappa^- u^- + \kappa^+ u^+}{\kappa^- + \kappa^+} \hat{n}\]

For an orthotropic material, uses an averaging wrt to the normal component according to

\[\kappa = n^T \cdot \kappa \cdot n\]
mirgecom.diffusion.diffusion_facial_flux_central(kappa_tpair, u_tpair, grad_u_tpair, lengthscales_tpair, normal, *, penalty_amount=None)[source]#

Compute the numerical flux for \(\nabla \cdot (\kappa \nabla u)\).

Uses a simple average of the two sides’ values:

\[F = -\frac{\kappa^- \nabla u^- + \kappa^+ \nabla u^+}{2} \cdot \hat{n} - \tau (u^+ - u^-).\]

The amount of penalization \(\tau\) is given by

\[\tau = \frac{\alpha \bar{\kappa}_{avg}}{l},\]

where \(\alpha\) is a user-definied value, \(l\) is the element characteristic lengthscale and \(\bar{\kappa}_{avg}\) is the averaged value, considering both isotropic (scalar) or orthotropic (array) cases. In the latter, the normal value is used for the penalization (see [Ern_2008]).

mirgecom.diffusion.diffusion_facial_flux_harmonic(kappa_tpair, u_tpair, grad_u_tpair, lengthscales_tpair, normal, *, penalty_amount=None)[source]#

Compute the numerical flux for \(\nabla \cdot (\kappa \nabla u)\).

Uses a modified average of the two sides’ values that replaces \(\kappa^-\) and \(\kappa^+\) with their harmonic mean, plus a penalization term

\[F = -\frac{2 \kappa_{ii}^- \kappa_{ii}^+}{\kappa_{ii}^- + \kappa_{ii}^+} \frac{\nabla u^- + \nabla u^+}{2} \cdot \hat{n} - \tau (u^+ - u^-).\]

The amount of penalization \(\tau\) is given by

\[\tau = \frac{\alpha \bar{\kappa}_{harm}}{l},\]

where \(\alpha\) is a user-defined value, \(l\) is the element characteristic lengthscale and \(\bar{\kappa}_{harm}\) is the harmonic mean, considering both isotropic (scalar) or orthotropic (array) cases. In the latter, the normal value is used for the penalization (see [Ern_2008]).

4.2.2. RHS Evaluation#

mirgecom.diffusion.grad_operator(dcoll, kappa, boundaries, u, *, quadrature_tag=<class 'grudge.dof_desc.DISCR_TAG_BASE'>, dd=DOFDesc(domain_tag=VolumeDomainTag(tag=<class 'grudge.dof_desc.VTAG_ALL'>), discretization_tag=<class 'grudge.dof_desc.DISCR_TAG_BASE'>), comm_tag=None, numerical_flux_func=<function grad_facial_flux_weighted>, kappa_tpairs=None, u_tpairs=None)[source]#

Compute the gradient of u.

Parameters:
Returns:

grad_u – the gradient of u

Return type:

numpy.ndarray

mirgecom.diffusion.diffusion_operator(dcoll, kappa, boundaries, u, *, return_grad_u=False, penalty_amount=None, gradient_numerical_flux_func=<function grad_facial_flux_weighted>, diffusion_numerical_flux_func=<function diffusion_facial_flux_harmonic>, quadrature_tag=<class 'grudge.dof_desc.DISCR_TAG_BASE'>, dd=DOFDesc(domain_tag=VolumeDomainTag(tag=<class 'grudge.dof_desc.VTAG_ALL'>), discretization_tag=<class 'grudge.dof_desc.DISCR_TAG_BASE'>), comm_tag=None, grad_u=None)[source]#

Compute the diffusion operator.

The diffusion operator is defined as \(\nabla\cdot(\kappa\nabla u)\), where \(\kappa\) is the conductivity and \(u\) is a scalar field.

Parameters:
  • dcoll (grudge.discretization.DiscretizationCollection) – the discretization collection to use

  • kappa (numbers.Number or meshmode.dof_array.DOFArray or numpy.ndarray) – the conductivity value(s)

  • boundaries – dictionary (or list of dictionaries) mapping boundary domain tags to DiffusionBoundary instances

  • u (meshmode.dof_array.DOFArray or numpy.ndarray) – the DOF array (or object array of DOF arrays) to which the operator should be applied

  • return_grad_u (bool) – an optional flag indicating whether \(\nabla u\) should also be returned

  • penalty_amount (float) – strength parameter for the diffusion flux interior penalty (temporary?); the default value is 0.05

  • gradient_numerical_flux_func – function that computes the numerical gradient flux

  • diffusion_numerical_flux_func – function that computes the numerical diffusive flux

  • quadrature_tag – quadrature tag indicating which discretization in dcoll to use for overintegration

  • dd (grudge.dof_desc.DOFDesc) – the DOF descriptor of the discretization on which u lives. Must be a volume on the base discretization.

  • comm_tag (Hashable) – Tag for distributed communication

Returns:

  • diff_u (meshmode.dof_array.DOFArray or numpy.ndarray) – the diffusion operator applied to u

  • grad_u (numpy.ndarray) – the gradient of u; only returned if return_grad_u is True

4.2.3. Boundary conditions#

class mirgecom.diffusion.DiffusionBoundary[source]#

Diffusion boundary base class.

abstract get_grad_flux(dcoll, dd_bdry, kappa_minus, u_minus, *, numerical_flux_func=<function grad_facial_flux_weighted>)[source]#

Compute the flux for grad(u) on the boundary dd_bdry.

abstract get_diffusion_flux(dcoll, dd_bdry, kappa_minus, u_minus, grad_u_minus, lengthscales_minus, *, penalty_amount=None, numerical_flux_func=<function diffusion_facial_flux_harmonic>)[source]#

Compute the flux for diff(u) on the boundary dd_bdry.

class mirgecom.diffusion.DirichletDiffusionBoundary(value)[source]#

Dirichlet boundary condition for the diffusion operator.

For the boundary condition \(u|_\Gamma = f\), uses external data

\[ \begin{align}\begin{aligned} u^+ &= 2 f - u^-\\(\nabla u)^+ &= (\nabla u)^-\end{aligned}\end{align} \]

to compute boundary fluxes as shown in [Hesthaven_2008], Section 7.1.

__init__(value)[source]#

Initialize the boundary condition.

Parameters:

value (float or meshmode.dof_array.DOFArray) – the value(s) of \(f\) along the boundary

get_grad_flux(dcoll, dd_bdry, kappa_minus, u_minus, *, numerical_flux_func=<function grad_facial_flux_weighted>)[source]#

Compute the flux for grad(u) on the boundary dd_bdry.

get_diffusion_flux(dcoll, dd_bdry, kappa_minus, u_minus, grad_u_minus, lengthscales_minus, *, penalty_amount=None, numerical_flux_func=<function diffusion_facial_flux_harmonic>)[source]#

Compute the flux for diff(u) on the boundary dd_bdry.

class mirgecom.diffusion.NeumannDiffusionBoundary(value)[source]#

Neumann boundary condition for the diffusion operator.

For the boundary condition \((\nabla u \cdot \mathbf{\hat{n}})|_\Gamma = g\), uses external data

\[u^+ = u^-\]

when computing the boundary fluxes for \(\nabla u\), and uses

\[ \begin{align}\begin{aligned}(-\kappa \nabla u\cdot\mathbf{\hat{n}})|_\Gamma &= -\kappa^- (\nabla u\cdot\mathbf{\hat{n}})|_\Gamma\\ &= -\kappa^- g\end{aligned}\end{align} \]

when computing the boundary fluxes for \(\nabla \cdot (\kappa \nabla u)\).

__init__(value)[source]#

Initialize the boundary condition.

Parameters:

value (float or meshmode.dof_array.DOFArray) – the value(s) of \(g\) along the boundary

get_grad_flux(dcoll, dd_bdry, kappa_minus, u_minus, *, numerical_flux_func=<function grad_facial_flux_weighted>)[source]#

Compute the flux for grad(u) on the boundary dd_bdry.

get_diffusion_flux(dcoll, dd_bdry, kappa_minus, u_minus, grad_u_minus, lengthscales_minus, *, penalty_amount=None, numerical_flux_func=<function diffusion_facial_flux_harmonic>)[source]#

Compute the flux for diff(u) on the boundary dd_bdry.

class mirgecom.diffusion.RobinDiffusionBoundary(u_ref, alpha)[source]#

Robin boundary condition for the diffusion operator.

The non-homogeneous Robin boundary condition is a linear combination of \(u\) and its gradient \(\nabla u\), given by

\[(\alpha u - \kappa \nabla u \cdot \mathbf{\hat{n}})|_\Gamma = \alpha u_{ref}.\]

where \(u_{ref}\) is the reference value of \(u\) for \(x \to \infty\) and \(\alpha\) is the weight for \(u\). The gradient weight \(\kappa\) is the conductivity (thermal) or the diffusivity (species).

The current implementation uses external data

\[u^+ = u^-\]

when computing the boundary fluxes for \(\nabla u\), and

\[\nabla u\cdot\mathbf{\hat{n}} |_\Gamma = \frac{\alpha}{\kappa}(u_{ref} - u^-)\]

when computing the boundary fluxes for \(\nabla \cdot (\kappa \nabla u)\).

__init__(u_ref, alpha)[source]#

Initialize the boundary condition.

Parameters:
get_grad_flux(dcoll, dd_bdry, kappa_minus, u_minus, *, numerical_flux_func=<function grad_facial_flux_weighted>)[source]#

Compute the flux for grad(u) on the boundary dd_bdry.

get_diffusion_flux(dcoll, dd_bdry, kappa_minus, u_minus, grad_u_minus, lengthscales_minus, *, penalty_amount=None, numerical_flux_func=<function diffusion_facial_flux_harmonic>)[source]#

Compute the flux for diff(u) on the boundary dd_bdry.

class mirgecom.diffusion.PrescribedFluxDiffusionBoundary(value)[source]#

Prescribed flux boundary condition for the diffusion operator.

For the boundary condition \((\nabla u \cdot \mathbf{\hat{n}})|_\Gamma\), uses external data

\[u^+ = u^-\]

when computing the boundary fluxes for \(\nabla u\), and applies directly the prescribed flux \(g\) when computing \(\nabla \cdot (\kappa \nabla u)\):

\[f_{presc} \cdot \hat{n}\]
__init__(value)[source]#

Initialize the boundary condition.

Parameters:

value (float or meshmode.dof_array.DOFArray) – the value(s) of \(g\) along the boundary

get_grad_flux(dcoll, dd_bdry, kappa_minus, u_minus, *, numerical_flux_func=<function grad_facial_flux_weighted>)[source]#

Compute the flux for grad(u) on the boundary dd_bdry.

get_diffusion_flux(dcoll, dd_bdry, kappa_minus, u_minus, grad_u_minus, lengthscales_minus, *, penalty_amount=None, numerical_flux_func=<function diffusion_facial_flux_harmonic>)[source]#

Compute the flux for diff(u) on the boundary dd_bdry.

class mirgecom.diffusion.DummyDiffusionBoundary[source]#

Dummy boundary condition that duplicates the internal values.