Skip to content

Commit

Permalink
Stab wildly at extricating HLL out of NS.
Browse files Browse the repository at this point in the history
  • Loading branch information
@MTCam
MTCam committed Apr 27, 2022
1 parent 3e1fb92 commit 50e1b6d
Show file tree
Hide file tree
Showing 8 changed files with 906 additions and 176 deletions.
3 changes: 3 additions & 0 deletions doc/misc.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,3 +71,6 @@ References
`(DOI) <https://doi.org/10.1007/978-3-642-59721-3_14>`__
.. [Ihme_2014] Yu Lv and Matthias Ihme (2014) Journal of Computationsl Physics 270 105 \
`(DOI) <http://dx.doi.org/10.1016/j.jcp.2014.03.029>`__
.. [Toro_2009] Eleuterio F. Toro (2009), Riemann Solvers and Numerical Methods for Fluid Dynamics, Springer \
`(DOI) <http://doi.org/10.1007/978-3-540-49834-6>`__
3 changes: 2 additions & 1 deletion mirgecom/boundary.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,7 +122,8 @@ def inviscid_divergence_flux(self, discr, btag, gas_model, state_minus,
boundary_state_pair = TracePair(btag, interior=state_minus,
exterior=state_plus)

return self._inviscid_div_flux_func(discr, state_tpair=boundary_state_pair)
return self._inviscid_div_flux_func(discr, state_pair=boundary_state_pair,
gas_model=gas_model)


class DummyBoundary(PrescribedFluidBoundary):
Expand Down
24 changes: 10 additions & 14 deletions mirgecom/euler.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,8 +63,10 @@

from mirgecom.inviscid import (
inviscid_flux,
inviscid_facial_flux
inviscid_flux_rusanov,
inviscid_boundary_flux_for_divergence_operator,
)

from mirgecom.operators import div_operator
from mirgecom.gas_model import (
project_fluid_state,
Expand All @@ -81,6 +83,7 @@ class _EulerTseedTag:


def euler_operator(discr, state, gas_model, boundaries, time=0.0,
inviscid_numerical_flux_func=inviscid_flux_rusanov,
quadrature_tag=None):
r"""Compute RHS of the Euler flow equations.
Expand Down Expand Up @@ -176,20 +179,13 @@ def _interp_to_surf_quad(utpair):

# Compute volume contributions
inviscid_flux_vol = inviscid_flux(vol_state_quad)

# Compute interface contributions
inviscid_flux_bnd = (

# Interior faces
sum(inviscid_facial_flux(discr, state_pair)
for state_pair in interior_states_quad)

# Domain boundary faces
+ sum(
boundaries[btag].inviscid_divergence_flux(
discr, as_dofdesc(btag).with_discr_tag(quadrature_tag), gas_model,
state_minus=boundary_states_quad[btag], time=time)
for btag in boundaries)
)
inviscid_flux_bnd = \
inviscid_boundary_flux_for_divergence_operator(
discr, gas_model, boundaries, interior_states_quad,
boundary_states_quad, quadrature_tag=quadrature_tag,
numerical_flux_func=inviscid_numerical_flux_func, time=time)

return -div_operator(discr, dd_quad_vol, dd_quad_faces,
inviscid_flux_vol, inviscid_flux_bnd)
Expand Down
243 changes: 137 additions & 106 deletions mirgecom/flux.py
View file
Original file line number Diff line number Diff line change
@@ -1,12 +1,17 @@
""":mod:`mirgecom.flux` provides inter-facial flux routines.
""":mod:`mirgecom.flux` provides generic inter-elemental flux routines.
Numerical Flux Routines
^^^^^^^^^^^^^^^^^^^^^^^
Low-level interfaces
^^^^^^^^^^^^^^^^^^^^
.. autofunction:: num_flux_lfr
.. autofunction:: num_flux_central
.. autofunction:: num_flux_hll
Flux pair interfaces for operators
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. autofunction:: gradient_flux_central
.. autofunction:: divergence_flux_central
.. autofunction:: flux_lfr
.. autofunction:: divergence_flux_lfr
"""

__copyright__ = """
Expand Down Expand Up @@ -35,179 +40,205 @@
import numpy as np # noqa


def gradient_flux_central(u_tpair, normal):
r"""Compute a central flux for the gradient operator.
def num_flux_lfr(f_minus_normal, f_plus_normal, q_minus, q_plus, lam, **kwargs):
r"""Compute Lax-Friedrichs/Rusanov flux after [Hesthaven_2008]_, Section 6.6.
The central gradient flux, $\mathbf{h}$, of a scalar quantity $u$ is calculated
as:
The Lax-Friedrichs/Rusanov flux is calculated as:
.. math::
f_{\text{LFR}}=\frac{1}{2}\left(f^-+f^+-\lambda\left(q^--q^+\right)\right)
\mathbf{h}({u}^-, {u}^+; \mathbf{n}) = \frac{1}{2}
\left({u}^{+}+{u}^{-}\right)\mathbf{\hat{n}}
where ${u}^-, {u}^+$, are the scalar function values on the interior
and exterior of the face on which the central flux is to be calculated, and
$\mathbf{\hat{n}}$ is the *normal* vector.
where $f^{\mp}$ and $q^{\mp}$ are the normal flux components and state components
on the interior and the exterior of the face over which the LFR flux is to be
calculated. $\lambda$ is the user-supplied dissipation/penalty coefficient.
*u_tpair* is the :class:`~grudge.trace_pair.TracePair` representing the scalar
quantities ${u}^-, {u}^+$. *u_tpair* may also represent a vector-quantity
:class:`~grudge.trace_pair.TracePair`, and in this case the central scalar flux
is computed on each component of the vector quantity as an independent scalar.
The $\lambda$ parameter is system-specific. Specifically, for the Rusanov flux
it is the max eigenvalue of the flux jacobian.
Parameters
----------
u_tpair: :class:`~grudge.trace_pair.TracePair`
Trace pair for the face upon which flux calculation is to be performed
normal: numpy.ndarray
object array of :class:`~meshmode.dof_array.DOFArray` with outward-pointing
normals
f_minus_normal
Normal component of physical flux interior to (left of) interface
f_plus_normal
Normal component of physical flux exterior to (right of) interface
q_minus
Physical state on interior of interface
q_plus
Physical state on exterior of interface
lam: :class:`~meshmode.dof_array.DOFArray`
State jump penalty parameter for dissipation term
Returns
-------
numpy.ndarray
object array of :class:`~meshmode.dof_array.DOFArray` with the flux for each
scalar component.
object array of :class:`~meshmode.dof_array.DOFArray` with the
Lax-Friedrichs/Rusanov numerical flux.
"""
from arraycontext import outer
return outer(u_tpair.avg, normal)
return (f_minus_normal + f_plus_normal - lam*(q_plus - q_minus))/2


def divergence_flux_central(trace_pair, normal):
r"""Compute a central flux for the divergence operator.
def num_flux_central(f_minus_normal, f_plus_normal, **kwargs):
r"""Central low-level numerical flux.
The central divergence flux, $h$, is calculated as:
The central flux is calculated as:
.. math::
h(\mathbf{v}^-, \mathbf{v}^+; \mathbf{n}) = \frac{1}{2}
\left(\mathbf{v}^{+}+\mathbf{v}^{-}\right) \cdot \hat{n}
where $\mathbf{v}^-, \mathbf{v}^+$, are the vectors on the interior and exterior
of the face across which the central flux is to be calculated, and $\hat{n}$ is
the unit normal to the face.
f_{\text{central}} = \frac{\left(f^++f^-\right)}{2}
Parameters
----------
trace_pair: :class:`~grudge.trace_pair.TracePair`
Trace pair for the face upon which flux calculation is to be performed
normal: numpy.ndarray
object array of :class:`~meshmode.dof_array.DOFArray` with outward-pointing
normals
f_minus_normal
Normal component of physical flux interior to (left of) interface
f_plus_normal
Normal component of physical flux exterior to (right of) interface
Returns
-------
numpy.ndarray
object array of :class:`~meshmode.dof_array.DOFArray` with the flux for each
scalar component.
object array of :class:`~meshmode.dof_array.DOFArray` with the
central numerical flux.
"""
return trace_pair.avg@normal
return (f_plus_normal + f_minus_normal)/2


def divergence_flux_lfr(cv_tpair, f_tpair, normal, lam):
r"""Compute Lax-Friedrichs/Rusanov flux after [Hesthaven_2008]_, Section 6.6.
def num_flux_hll(f_minus_normal, f_plus_normal, q_minus, q_plus, s_minus, s_plus):
r"""HLL low-level numerical flux.
The Lax-Friedrichs/Rusanov flux is calculated as:
The Harten, Lax, van Leer approximate Riemann numerical flux is calculated as:
.. math::
f_{\mathtt{LFR}} = \frac{1}{2}(\mathbf{F}(q^-) + \mathbf{F}(q^+)) \cdot
\hat{n} + \frac{\lambda}{2}(q^{-} - q^{+}),
f^{*}_{\text{HLL}} = \frac{\left(s^+f^--s^-f^++s^+s^-\left(q^+-q^-\right)
\right)}{\left(s^+ - s^-\right)}
where $q^-, q^+$ are the scalar solution components on the interior and the
exterior of the face on which the LFR flux is to be calculated, $\mathbf{F}$ is
the vector flux function, $\hat{n}$ is the face normal, and $\lambda$ is the
user-supplied jump term coefficient.
where $f^{\mp}$, $q^{\mp}$, and $s^{\mp}$ are the interface-normal fluxes, the
states, and the wavespeeds for the interior (-) and exterior (+) of the
interface, respectively.
The $\lambda$ parameter is system-specific. Specifically, for the Rusanov flux
it is the max eigenvalue of the flux Jacobian:
.. math::
\lambda = \text{max}\left(|\mathbb{J}_{F}(q^-)|,|\mathbb{J}_{F}(q^+)|\right)
Here, $\lambda$ is a function parameter, leaving the responsibility for the
calculation of the eigenvalues of the system-dependent flux Jacobian to the
caller.
Details about this approximate Riemann solver can be found in Section 10.3 of
[Toro_2009]_.
Parameters
----------
cv_tpair: :class:`~grudge.trace_pair.TracePair`
f_minus_normal
Normal component of physical flux interior to (left of) interface
Solution trace pair for faces for which numerical flux is to be calculated
f_plus_normal
Normal component of physical flux exterior to (right of) interface
f_tpair: :class:`~grudge.trace_pair.TracePair`
q_minus
Physical state on interior of interface
Physical flux trace pair on faces on which numerical flux is to be calculated
q_plus
Physical state on exterior of interface
normal: numpy.ndarray
q_minus
Physical state on interior of interface
object array of :class:`meshmode.dof_array.DOFArray` with outward-pointing
normals
q_plus
Physical state on exterior of interface
lam: :class:`~meshmode.dof_array.DOFArray`
s_minus: :class:`~meshmode.dof_array.DOFArray`
Interface wave-speed parameter for interior of interface
lambda parameter for Lax-Friedrichs/Rusanov flux
s_plus: :class:`~meshmode.dof_array.DOFArray`
Interface wave-speed parameter for exterior of interface
Returns
-------
numpy.ndarray
object array of :class:`~meshmode.dof_array.DOFArray` with the
Lax-Friedrichs/Rusanov flux.
HLL numerical flux.
"""
return flux_lfr(cv_tpair, f_tpair, normal, lam)@normal
actx = q_minus.array_context
f_star = (s_plus*f_minus_normal - s_minus*f_plus_normal
+ s_plus*s_minus*(q_plus - q_minus))/(s_plus - s_minus)

# choose the correct f contribution based on the wave speeds
f_check_minus = \
actx.np.greater_equal(s_minus, 0*s_minus)*(0*f_minus_normal + 1.0)
f_check_plus = \
actx.np.less_equal(s_plus, 0*s_plus)*(0*f_minus_normal + 1.0)

def flux_lfr(cv_tpair, f_tpair, normal, lam):
r"""Compute Lax-Friedrichs/Rusanov flux after [Hesthaven_2008]_, Section 6.6.
The Lax-Friedrichs/Rusanov flux is calculated as:
f = f_star
f = actx.np.where(f_check_minus, f_minus_normal, f)
f = actx.np.where(f_check_plus, f_plus_normal, f)

.. math::
return f

f_{\mathtt{LFR}} = \frac{1}{2}(\mathbf{F}(q^-) + \mathbf{F}(q^+))
+ \frac{\lambda}{2}(q^{-} - q^{+})\hat{\mathbf{n}},

where $q^-, q^+$ are the scalar solution components on the interior and the
exterior of the face on which the LFR flux is to be calculated, $\mathbf{F}$ is
the vector flux function, $\hat{\mathbf{n}}$ is the face normal, and $\lambda$
is the user-supplied jump term coefficient.
def gradient_flux_central(u_tpair, normal):
r"""Compute a central flux for the gradient operator.
The $\lambda$ parameter is system-specific. Specifically, for the Rusanov flux
it is the max eigenvalue of the flux jacobian:
The central gradient flux, $\mathbf{h}$, of a scalar quantity $u$ is calculated
as:
.. math::
\lambda = \text{max}\left(|\mathbb{J}_{F}(q^-)|,|\mathbb{J}_{F}(q^+)|\right)
Here, $\lambda$ is a function parameter, leaving the responsibility for the
calculation of the eigenvalues of the system-dependent flux Jacobian to the
caller.
\mathbf{h}({u}^-, {u}^+; \mathbf{n}) = \frac{1}{2}
\left({u}^{+}+{u}^{-}\right)\mathbf{\hat{n}}
where ${u}^-, {u}^+$, are the scalar function values on the interior
and exterior of the face on which the central flux is to be calculated, and
$\mathbf{\hat{n}}$ is the *normal* vector.
*u_tpair* is the :class:`~grudge.trace_pair.TracePair` representing the scalar
quantities ${u}^-, {u}^+$. *u_tpair* may also represent a vector-quantity
:class:`~grudge.trace_pair.TracePair`, and in this case the central scalar flux
is computed on each component of the vector quantity as an independent scalar.
Parameters
----------
cv_tpair: :class:`~grudge.trace_pair.TracePair`
u_tpair: :class:`~grudge.trace_pair.TracePair`
Trace pair for the face upon which flux calculation is to be performed
normal: numpy.ndarray
object array of :class:`~meshmode.dof_array.DOFArray` with outward-pointing
normals
Returns
-------
numpy.ndarray
object array of :class:`~meshmode.dof_array.DOFArray` with the flux for each
scalar component.
"""
from arraycontext import outer
return outer(u_tpair.avg, normal)

Solution trace pair for faces for which numerical flux is to be calculated

f_tpair: :class:`~grudge.trace_pair.TracePair`
def divergence_flux_central(trace_pair, normal):
r"""Compute a central flux for the divergence operator.
Physical flux trace pair on faces on which numerical flux is to be calculated
The central divergence flux, $h$, is calculated as:
normal: numpy.ndarray
.. math::
object array of :class:`~meshmode.dof_array.DOFArray` with outward-pointing
normals
h(\mathbf{v}^-, \mathbf{v}^+; \mathbf{n}) = \frac{1}{2}
\left(\mathbf{v}^{+}+\mathbf{v}^{-}\right) \cdot \hat{n}
lam: :class:`~meshmode.dof_array.DOFArray`
where $\mathbf{v}^-, \mathbf{v}^+$, are the vectors on the interior and exterior
of the face across which the central flux is to be calculated, and $\hat{n}$ is
the unit normal to the face.
lambda parameter for Lax-Friedrichs/Rusanov flux
Parameters
----------
trace_pair: :class:`~grudge.trace_pair.TracePair`
Trace pair for the face upon which flux calculation is to be performed
normal: numpy.ndarray
object array of :class:`~meshmode.dof_array.DOFArray` with outward-pointing
normals
Returns
-------
numpy.ndarray
object array of :class:`~meshmode.dof_array.DOFArray` with the
Lax-Friedrichs/Rusanov flux.
object array of :class:`~meshmode.dof_array.DOFArray` with the flux for each
scalar component.
"""
from arraycontext import outer
return f_tpair.avg - lam*outer(cv_tpair.diff, normal)/2
return trace_pair.avg@normal
Loading

0 comments on commit 50e1b6d

Please sign in to comment.