DG solver extension

josie/frac_mom/fields.py

@@ -0,0 +1,55 @@
+# josiepy
+# Copyright © 2021 Ruben Di Battista
+# All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+# 1. Redistributions of source code must retain the above copyright
+# notice, this list of conditions and the following disclaimer.
+# 2. Redistributions in binary form must reproduce the above copyright
+# notice, this list of conditions and the following disclaimer in the
+# documentation and/or other materials provided with the distribution.
+#
+# THIS SOFTWARE IS PROVIDED BY Ruben Di Battista ''AS IS'' AND ANY
+# EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+# DISCLAIMED. IN NO EVENT SHALL Ruben Di Battista BE LIABLE FOR ANY
+# DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+#
+# The views and conclusions contained in the software and documentation
+# are those of the authors and should not be interpreted as representing
+# official policies, either expressed or implied, of Ruben Di Battista.
+
+from __future__ import annotations
+
+from josie.fluid.fields import FluidFields
+from josie.state import Fields
+
+
+class ConsFields(Fields):
+    """Indexing enum for the conservative state variables of the problem"""
+
+    m0 = 0
+    m12 = 1
+    m1 = 2
+    m32 = 3
+    m1U = 4
+    m1V = 5
+
+
+class FracMomFields(FluidFields):
+    """Indexing enum for the state variables of the problem"""
+
+    m0 = 0
+    m12 = 1
+    m1 = 2
+    m32 = 3
+    m1U = 4
+    m1V = 5
+    U = 6
+    V = 7

josie/frac_mom/problem.py

@@ -0,0 +1,130 @@
+# josiepy
+# Copyright © 2019 Ruben Di Battista
+# All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+# 1. Redistributions of source code must retain the above copyright
+# notice, this list of conditions and the following disclaimer.
+# 2. Redistributions in binary form must reproduce the above copyright
+# notice, this list of conditions and the following disclaimer in the
+# documentation and/or other materials provided with the distribution.
+#
+# THIS SOFTWARE IS PROVIDED BY Ruben Di Battista ''AS IS'' AND ANY
+# EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+# DISCLAIMED. IN NO EVENT SHALL Ruben Di Battista BE LIABLE FOR ANY
+# DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+#
+# The views and conclusions contained in the software and documentation
+# are those of the authors and should not be interpreted as representing
+# official policies, either expressed or implied, of Ruben Di Battista.
+from __future__ import annotations
+
+import numpy as np
+
+
+from josie.dimension import MAX_DIMENSIONALITY
+from josie.math import Direction
+from josie.mesh.cellset import CellSet
+from josie.problem import ConvectiveProblem
+
+from .state import Q
+from .fields import ConsFields
+
+
+class FracMomProblem(ConvectiveProblem):
+    """A class representing an PGD system problem
+
+    Attributes
+    ---------
+    eos
+        An instance of :class:`~.EOS`, an equation of state for the fluid
+    """
+
+    def F(self, cells: CellSet) -> np.ndarray:
+        r""" This returns the tensor representing the flux for an PGD model
+
+        A general problem can be written in a compact way:
+
+        .. math::
+
+            \pdv{\vb{q}}{t} + \div{\vb{F\qty(\vb{q})}} + \vb{B}\qty(\vb{q})
+            \cdot \gradient{\vb{q}} = \vb{s\qty(\vb{q})}
+
+        This function needs to return :math:`\vb{F}\qty(\vb{q})`
+
+        Parameters
+        ----------
+        cells
+            A :class:`MeshCellSet` that contains the cell data
+
+        Returns
+        ---------
+        F
+            An array of dimension :math:`Nx \times Ny \times 4 \times 2`, i.e.
+            an array that of each :math:`x` cell and :math:`y` cell stores the
+            :math:`4 \times 2` flux tensor
+
+            The flux tensor is:
+
+            .. math::
+
+                \pdeConvective =
+                \begin{bmatrix}
+                    \rho u & \rho v \\
+                    \rho u^2 + p & \rho uv \\
+                    \rho vu & \rho v^ 2 + p \\
+                    (\rho E + p)U & (\rho E + p)V
+                \end{bmatrix}
+        """
+        values: Q = cells.values.view(Q)
+        fields = values.fields
+
+        num_cells_x, num_cells_y, num_dofs, _ = values.shape
+
+        # Flux tensor
+        F = np.empty(
+            (
+                num_cells_x,
+                num_cells_y,
+                num_dofs,
+                len(ConsFields),
+                MAX_DIMENSIONALITY,
+            )
+        )
+        U = values[..., fields.U]
+        V = values[..., fields.V]
+        m0U = np.multiply(values[..., fields.m0], U)
+        m0V = np.multiply(values[..., fields.m0], V)
+        m12U = np.multiply(values[..., fields.m12], U)
+        m12V = np.multiply(values[..., fields.m12], V)
+        m1U = values[..., fields.m1U]
+        m1V = values[..., fields.m1V]
+        m32U = np.multiply(values[..., fields.m32], U)
+        m32V = np.multiply(values[..., fields.m32], V)
+
+        m1UU = np.multiply(m1U, U)
+        m1UV = np.multiply(m1U, V)
+        m1VV = np.multiply(m1V, V)
+        m1VU = m1UV
+
+        F[..., fields.m0, Direction.X] = m0U
+        F[..., fields.m0, Direction.Y] = m0V
+        F[..., fields.m12, Direction.X] = m12U
+        F[..., fields.m12, Direction.Y] = m12V
+        F[..., fields.m1, Direction.X] = m1U
+        F[..., fields.m1, Direction.Y] = m1V
+        F[..., fields.m32, Direction.X] = m32U
+        F[..., fields.m32, Direction.Y] = m32V
+        F[..., fields.m1U, Direction.X] = m1UU
+        F[..., fields.m1U, Direction.Y] = m1UV
+        F[..., fields.m1V, Direction.X] = m1VU
+        F[..., fields.m1V, Direction.Y] = m1VV
+
+        return F

josie/frac_mom/schemes/LF.py

@@ -0,0 +1,131 @@
+# josiepy
+# Copyright © 2020 Ruben Di Battista
+# All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+# 1. Redistributions of source code must retain the above copyright
+# notice, this list of conditions and the following disclaimer.
+# 2. Redistributions in binary form must reproduce the above copyright
+# notice, this list of conditions and the following disclaimer in the
+# documentation and/or other materials provided with the distribution.
+#
+# THIS SOFTWARE IS PROVIDED BY Ruben Di Battista ''AS IS'' AND ANY
+# EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+# DISCLAIMED. IN NO EVENT SHALL Ruben Di Battista BE LIABLE FOR ANY
+# DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+#
+# The views and conclusions contained in the software and documentation
+# are those of the authors and should not be interpreted as representing
+# official policies, either expressed or implied, of Ruben Di Battista.
+from __future__ import annotations
+
+import numpy as np
+
+from typing import TYPE_CHECKING
+from josie.frac_mom.state import Q
+
+from .scheme import FracMomScheme
+
+if TYPE_CHECKING:
+    from josie.mesh.cellset import NeighboursCellSet, MeshCellSet
+
+
+class LF(FracMomScheme):
+    r"""This class implements the Rusanov scheme. See
+    :cite:`toro_riemann_2009` for a detailed view on compressible schemes.
+    The Rusanov scheme is discretized by:
+
+    .. math::
+
+        \numConvective  =
+            \frac{1}{2} \qty[%
+            \qty|\pdeConvective|_{i+1} + \qty|\pdeConvective|_{i}
+            - \sigma \qty(\pdeState_{i+1} - \pdeState_{i})
+            ] S_f
+    """
+
+    def F(self, cells: MeshCellSet, neighs: NeighboursCellSet):
+        Q_L: Q = cells.values.view(Q)
+        Q_R: Q = neighs.values.view(Q)
+
+        fields = Q.fields
+
+        FS = np.zeros_like(Q_L).view(Q)
+        DeltaF = np.zeros_like(self.problem.F(cells)).view(Q)
+
+        u_L = Q_L[..., fields.U]
+        u_R = Q_R[..., fields.U]
+        v_L = Q_L[..., fields.V]
+        v_R = Q_R[..., fields.V]
+
+        sigma = max(
+            np.amax(np.abs(u_L)),
+            np.amax(np.abs(u_R)),
+            np.amax(np.abs(v_R)),
+            np.amax(np.abs(v_L)),
+        )
+
+        # First four variables of the total state are the conservative
+        # variables (rho, rhoU, rhoV, rhoE)
+        Q_L_cons = Q_L.get_conservative()
+        Q_R_cons = Q_R.get_conservative()
+        DeltaQ = np.zeros_like(Q_L_cons).view(Q)
+
+        if neighs.direction == 0:
+            DeltaF[..., 0:2, :, :] = 0.5 * (
+                self.problem.F(cells)[..., 0:2, :, :]
+                + self.problem.F(neighs)[..., 2:4, :, :]
+            )
+            DeltaQ[..., 0:2, :] = (
+                0.5 * sigma * (Q_R_cons[..., 2:4, :] - Q_L_cons[..., 0:2, :])
+            )
+        if neighs.direction == 1:
+            DeltaF[..., 0, :, :] = 0.5 * (
+                self.problem.F(cells)[..., 0, :, :]
+                + self.problem.F(neighs)[..., 1, :, :]
+            )
+            DeltaQ[..., 0, :] = (
+                0.5 * sigma * (Q_R_cons[..., 1, :] - Q_L_cons[..., 0, :])
+            )
+            DeltaF[..., 2, :, :] = 0.5 * (
+                self.problem.F(cells)[..., 2, :, :]
+                + self.problem.F(neighs)[..., 3, :, :]
+            )
+            DeltaQ[..., 2, :] = (
+                0.5 * sigma * (Q_R_cons[..., 3, :] - Q_L_cons[..., 2, :])
+            )
+        if neighs.direction == 2:
+            DeltaF[..., 2:4, :, :] = 0.5 * (
+                self.problem.F(cells)[..., 2:4, :, :]
+                + self.problem.F(neighs)[..., 0:2, :, :]
+            )
+
+            DeltaQ[..., 2:4, :] = (
+                0.5 * sigma * (Q_R_cons[..., 0:2, :] - Q_L_cons[..., 2:4, :])
+            )
+        if neighs.direction == 3:
+            DeltaF[..., 1, :, :] = 0.5 * (
+                self.problem.F(cells)[..., 1, :, :]
+                + self.problem.F(neighs)[..., 0, :, :]
+            )
+            DeltaQ[..., 1, :] = (
+                0.5 * sigma * (Q_R_cons[..., 0, :] - Q_L_cons[..., 1, :])
+            )
+            DeltaF[..., 3, :, :] = 0.5 * (
+                self.problem.F(cells)[..., 3, :, :]
+                + self.problem.F(neighs)[..., 2, :, :]
+            )
+            DeltaQ[..., 3, :] = (
+                0.5 * sigma * (Q_R_cons[..., 2, :] - Q_L_cons[..., 3, :])
+            )
+
+        Delta = np.einsum("...mkl,...l->...mk", DeltaF, neighs.normals)
+        FS.set_conservative(Delta - DeltaQ)
+        return FS

josie/frac_mom/schemes/__init__.py

@@ -0,0 +1,2 @@
+from .LF import LF
+from .scheme import FracMomScheme