RusanovWithExplicitEuler.py

Go to the documentation of this file.

# This file is part of the ExaHyPE2 project. For conditions of distribution and
# use, please see the copyright notice at www.peano-framework.org
import os
import re
 
import peano4.api
import peano4.api.datamodel
import peano4.api.dastgen2
import toolbox.blockstructured.api
import toolbox.loadbalancing.api.actionsets
import exahype2.api
from exahype2.api import Target
import exahype2.api.dslhype as dslhype
 
from ..Solver import Solver
from ..blockstructured.BoundaryConditions import get_access_to_persistent_face_data
from ..blockstructured.create_dynamic_mesh_refinement_action_set import (
    create_dynamic_mesh_refinement_action_set,
)
from .actionsets import ApplyDirichletBoundaryConditions
from .actionsets import ProjectPatchOntoFaces
from .actionsets import RollOverUpdatedFaceWithSommerfeldBoundaryConditions
from .actionsets import UpdateCell
from .actionsets import ComputeGlobalCFLCondition
 
 
 
 
 
 
def timeStepWithRusanov2D(patchData: "::exahype2::CellData<double, double>&") -> "void":
    N = Integer("patchData.numberOfCells")
    t = DataBlock([[0, N]], "patchData.t", "double")
    dt = DataBlock([[0, N]], "patchData.dt", "double")
    h = DataBlock([[0, N]], "patchData.cellSize", "tarch::la::Vector<2, double>")
    x = DataBlock([[0, N]], "patchData.cellCentre", "tarch::la::Vector<2, double>")
    maxEigenvalues = DataBlock([[0, N]], "patchData.maxEigenvalue", "double")
 
    input = DataBlock([[0, NumberOfUnknowns + NumberOfAuxiliaryVariables], [-H, NV + H], [-H, NV + H], [0, N]], "patchData.QIn")
    output = DataBlock([[0, NumberOfUnknowns + NumberOfAuxiliaryVariables], [0, NV], [0, NV], [0, N]], "patchData.QOut")
    output = input[:, 0:NV, 0:NV]
 
    if volumetricSource:
       source = DataBlock([[0, NumberOfUnknowns], [0, NV], [0, NV], [0, N]], eval("volumetricSource", input[:, 0:NV, 0:NV], x, t))
       output[0:NumberOfUnknowns, :, :] = output[0:NumberOfUnknowns, :, :] + dt * source
 
    # The Rusanov flux at a face needs the eigenvalues of both adjacent
    # volumes, i.e. the patch volumes plus the one halo layer across each
    # face along the normal. We never evaluate the remaining halo (edges,
    # corners), which holds no valid data.
    maxEigenvaluesX = DataBlock([[0, 1], [-1, NV + 1], [0, NV], [0, N]], eval("maxEigenvalue", input[:, -1:NV+1, 0:NV], x, 0))
    maxEigenvaluesY = DataBlock([[0, 1], [0, NV], [-1, NV + 1], [0, N]], eval("maxEigenvalue", input[:, 0:NV, -1:NV+1], x, 1))
 
    lambdaLeft = max(maxEigenvaluesX[0:1, -1:NV-1, 0:NV], maxEigenvaluesX[0:1, 0:NV, 0:NV])
    lambdaRight = max(maxEigenvaluesX[0:1, 0:NV, 0:NV], maxEigenvaluesX[0:1, 1:NV + 1, 0:NV])
    lambdaBottom = max(maxEigenvaluesY[0:1, 0:NV, -1:NV-1], maxEigenvaluesY[0:1, 0:NV, 0:NV])
    lambdaTop = max(maxEigenvaluesY[0:1, 0:NV, 0:NV], maxEigenvaluesY[0:1, 0:NV, 1:NV + 1])
 
    fluxLeft = -1.0 * lambdaLeft * (input[0:NumberOfUnknowns, 0:NV, 0:NV] - input[0:NumberOfUnknowns, -1:NV-1, 0:NV])
    fluxRight = -1.0 * lambdaRight * (input[0:NumberOfUnknowns, 1:NV+1, 0:NV] - input[0:NumberOfUnknowns, 0:NV, 0:NV])
    output[0:NumberOfUnknowns, :, :] = output[0:NumberOfUnknowns, :, :] + dt / h(0) * (fluxLeft - fluxRight) * NV * 0.5
 
    fluxBottom = -1.0 * lambdaBottom * (input[0:NumberOfUnknowns, 0:NV, 0:NV] - input[0:NumberOfUnknowns, 0:NV, -1:NV-1])
    fluxTop = -1.0 * lambdaTop * (input[0:NumberOfUnknowns, 0:NV, 1:NV+1] - input[0:NumberOfUnknowns, 0:NV, 0:NV])
    output[0:NumberOfUnknowns, :, :] = output[0:NumberOfUnknowns, :, :] + dt / h(1) * (fluxBottom - fluxTop) * NV * 0.5
 
    if flux:
      flux_x = DataBlock([[0, NumberOfUnknowns], [-H, NV + H], [0, NV], [0, N]], eval("flux", input[:, -1:NV+1, 0:NV], x, 0))
      flux_y = DataBlock([[0, NumberOfUnknowns], [0, NV], [-H, NV + H], [0, N]], eval("flux", input[:, 0:NV, -1:NV+1], x, 1))
      
      output[0:NumberOfUnknowns, :, :] = output[0:NumberOfUnknowns, :, :] + 0.5 * NV * dt / h(0) *  (flux_x[:, -1:NV-1, 0:NV] - flux_x[:, 1:NV+1, 0:NV])
      output[0:NumberOfUnknowns, :, :] = output[0:NumberOfUnknowns, :, :] + 0.5 * NV * dt / h(1) * (flux_y[:, 0:NV, -1:NV-1] - flux_y[:, 0:NV, 1:NV+1])
 
    if nonconservativeProduct:
       average_x = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, -1:NV, 0:NV] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV + 1, 0:NV])
       delta_x = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV + H, 0:NV] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, -H:NV, 0:NV]
       nonconservative_x = DataBlock([[0, NumberOfUnknowns], [0, NV + H], [0, NV], [0, N]], eval("nonconservativeProduct", average_x, delta_x, x, h, t, dt, 0))
       output[0:NumberOfUnknowns, :, :] = output[0:NumberOfUnknowns, :, :] + 0.5 * NV * dt / h(0) * (-nonconservative_x[:, 0:NV, 0:NV] - nonconservative_x[:, 1:NV+1, 0:NV])
       
       average_y = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, -1:NV] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV + 1])
       delta_y = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV + H] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, -H:NV]
       nonconservative_y = DataBlock([[0, NumberOfUnknowns], [0, NV], [0, NV + H], [0, N]], eval("nonconservativeProduct", average_y, delta_y, x, h, t, dt, 1))
       output[0:NumberOfUnknowns, :, :] = output[0:NumberOfUnknowns, :, :] + 0.5 * NV * dt / h(1) * (-nonconservative_y[:, 0:NV, 0:NV] - nonconservative_y[:, 0:NV, 1:NV+1])
 
    if maxEigenvalue:
      maxEigenvaluesX2 = DataBlock([[0, 1], [0, NV], [0, NV], [0, N]], eval("maxEigenvalue", output, x, 0))
      maxEigenvaluesY2 = DataBlock([[0, 1], [0, NV], [0, NV], [0, N]], eval("maxEigenvalue", output, x, 1))
      maxEigenvaluesPerVolume = max(maxEigenvaluesX2, maxEigenvaluesY2)
      maxEigenvalues = max(maxEigenvaluesPerVolume)
 
 
 
def timeStepWithRusanov3D(patchData: "::exahype2::CellData<double, double>&") -> "void":
    N = Integer("patchData.numberOfCells")
    t = DataBlock([[0, N]], "patchData.t", "double")
    dt = DataBlock([[0, N]], "patchData.dt", "double")
    h = DataBlock([[0, N]], "patchData.cellSize", "tarch::la::Vector<3, double>")
    x = DataBlock([[0, N]], "patchData.cellCentre", "tarch::la::Vector<3, double>")
    maxEigenvalues = DataBlock([[0, N]], "patchData.maxEigenvalue", "double")
 
    input = DataBlock([[0, NumberOfUnknowns + NumberOfAuxiliaryVariables], [-H, NV + H], [-H, NV + H], [-H, NV + H], [0, N]], "patchData.QIn")
    output = DataBlock([[0, NumberOfUnknowns + NumberOfAuxiliaryVariables], [0, NV], [0, NV], [0, NV], [0, N]], "patchData.QOut")
    output = input[:, 0:NV, 0:NV, 0:NV]
 
    if volumetricSource:
       source = DataBlock([[0, NumberOfUnknowns], [0, NV], [0, NV], [0, NV], [0, N]], eval("volumetricSource", input[:, 0:NV, 0:NV, 0:NV], x, t))
       output[0:NumberOfUnknowns, :, :, :] = output[0:NumberOfUnknowns, :, :, :] + dt * source
 
    # As in 2d: patch volumes plus the one halo layer across each face along
    # the normal only, never the edges and corners of the halo
    maxEigenvaluesX = DataBlock([[0, 1], [-1, NV + 1], [0, NV], [0, NV], [0, N]], eval("maxEigenvalue", input[:, -1:NV+1, 0:NV, 0:NV], x, 0))
    maxEigenvaluesY = DataBlock([[0, 1], [0, NV], [-1, NV + 1], [0, NV], [0, N]], eval("maxEigenvalue", input[:, 0:NV, -1:NV+1, 0:NV], x, 1))
    maxEigenvaluesZ = DataBlock([[0, 1], [0, NV], [0, NV], [-1, NV + 1], [0, N]], eval("maxEigenvalue", input[:, 0:NV, 0:NV, -1:NV+1], x, 2))
    
    lambdaLeft = max(maxEigenvaluesX[0:1, -1:NV-1, 0:NV, 0:NV], maxEigenvaluesX[0:1, 0:NV, 0:NV, 0:NV])
    lambdaRight = max(maxEigenvaluesX[0:1, 0:NV, 0:NV, 0:NV], maxEigenvaluesX[0:1, 1:NV + 1, 0:NV, 0:NV])
    lambdaBottom = max(maxEigenvaluesY[0:1, 0:NV, -1:NV-1, 0:NV], maxEigenvaluesY[0:1, 0:NV, 0:NV, 0:NV])
    lambdaTop = max(maxEigenvaluesY[0:1, 0:NV, 0:NV, 0:NV], maxEigenvaluesY[0:1, 0:NV, 1:NV + 1, 0:NV])
    lambdaNear = max(maxEigenvaluesZ[0:1, 0:NV, 0:NV, -1:NV-1], maxEigenvaluesZ[0:1, 0:NV, 0:NV, 0:NV])
    lambdaFar = max(maxEigenvaluesZ[0:1, 0:NV, 0:NV, 0:NV], maxEigenvaluesZ[0:1, 0:NV, 0:NV, 1:NV + 1])
 
    fluxLeft = -1.0 * lambdaLeft * (input[0:NumberOfUnknowns, 0:NV, 0:NV, 0:NV] - input[0:NumberOfUnknowns, -1:NV-1, 0:NV, 0:NV])
    fluxRight = -1.0 * lambdaRight * (input[0:NumberOfUnknowns, 1:NV+1, 0:NV, 0:NV] - input[0:NumberOfUnknowns, 0:NV, 0:NV, 0:NV])
    output[0:NumberOfUnknowns, :, :, :] = output[0:NumberOfUnknowns, :, :, :] + dt / h(0) * (fluxLeft - fluxRight) * NV * 0.5
    fluxBottom = -1.0 * lambdaBottom * (input[0:NumberOfUnknowns, 0:NV, 0:NV, 0:NV] - input[0:NumberOfUnknowns, 0:NV, -1:NV-1, 0:NV])
    fluxTop = -1.0 * lambdaTop * (input[0:NumberOfUnknowns, 0:NV, 1:NV+1, 0:NV] - input[0:NumberOfUnknowns, 0:NV, 0:NV, 0:NV])
    output[0:NumberOfUnknowns, :, :, :] = output[0:NumberOfUnknowns, :, :, :] + dt / h(1) * (fluxBottom - fluxTop) * NV * 0.5
 
    fluxNear = -1.0 * lambdaNear * (input[0:NumberOfUnknowns, 0:NV, 0:NV, 0:NV] - input[0:NumberOfUnknowns, 0:NV, 0:NV, -1:NV-1])
    fluxFar = -1.0 * lambdaFar * (input[0:NumberOfUnknowns, 0:NV, 0:NV, 1:NV+1] - input[0:NumberOfUnknowns, 0:NV, 0:NV, 0:NV])
    output[0:NumberOfUnknowns, :, :, :] = output[0:NumberOfUnknowns, :, :, :] + dt / h(2) * (fluxNear - fluxFar) * NV * 0.5
 
    if flux:
      flux_x = DataBlock([[0, NumberOfUnknowns], [-H, NV + H], [0, NV], [0, NV], [0, N]], eval("flux", input[:, -1:NV+1, 0:NV, 0:NV], x, 0))
      flux_y = DataBlock([[0, NumberOfUnknowns], [0, NV], [-H, NV + H], [0, NV], [0, N]], eval("flux", input[:, 0:NV, -1:NV+1, 0:NV], x, 1))
      flux_z = DataBlock([[0, NumberOfUnknowns], [0, NV], [0, NV], [-H, NV + H], [0, N]], eval("flux", input[:, 0:NV, 0:NV, -1:NV+1], x, 2))
      output[0:NumberOfUnknowns, :, :, :] = output[0:NumberOfUnknowns, :, :, :] + 0.5 * NV * dt / h(0) *  (flux_x[:, -1:NV-1, 0:NV, 0:NV] - flux_x[:, 1:NV+1, 0:NV, 0:NV])
      output[0:NumberOfUnknowns, :, :, :] = output[0:NumberOfUnknowns, :, :, :] + 0.5 * NV * dt / h(1) * (flux_y[:, 0:NV, -1:NV-1, 0:NV] - flux_y[:, 0:NV, 1:NV+1, 0:NV])
      output[0:NumberOfUnknowns, :, :, :] = output[0:NumberOfUnknowns, :, :, :] + 0.5 * NV * dt / h(2) * (flux_z[:, 0:NV, 0:NV, -1:NV-1] - flux_z[:, 0:NV, 0:NV, 1:NV+1])
 
    if nonconservativeProduct:
       average_x = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, -1:NV, 0:NV, 0:NV] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV + 1, 0:NV, 0:NV])
       delta_x = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV + H, 0:NV, 0:NV] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, -H:NV, 0:NV, 0:NV]
       nonconservative_x = DataBlock([[0, NumberOfUnknowns], [0, NV + H], [0, NV], [0, NV], [0, N]], eval("nonconservativeProduct", average_x, delta_x, x, h, t, dt, 0))
       output[0:NumberOfUnknowns, :, :, :] = output[0:NumberOfUnknowns, :, :, :] + 0.5 * NV * dt / h(0) * (-nonconservative_x[:, 0:NV, 0:NV, 0:NV] - nonconservative_x[:, 1:NV+1, 0:NV, 0:NV])
       
       average_y = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, -1:NV, 0:NV] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV + 1, 0:NV])
       delta_y = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV + H, 0:NV] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, -H:NV, 0:NV]
       nonconservative_y = DataBlock([[0, NumberOfUnknowns], [0, NV], [0, NV + H], [0, NV], [0, N]], eval("nonconservativeProduct", average_y, delta_y, x, h, t, dt, 1))
       output[0:NumberOfUnknowns, :, :, :] = output[0:NumberOfUnknowns, :, :, :] + 0.5 * NV * dt / h(1) * (-nonconservative_y[:, 0:NV, 0:NV, 0:NV] - nonconservative_y[:, 0:NV, 1:NV+1, 0:NV])
 
       average_z = 0.5 * (input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV, -1:NV] + input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV, 0:NV + 1])
       delta_z = input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV, 0:NV + H] - input[0:NumberOfUnknowns + NumberOfAuxiliaryVariables, 0:NV, 0:NV, -H:NV]
       nonconservative_z = DataBlock([[0, NumberOfUnknowns], [0, NV], [0, NV], [0, NV + H], [0, N]], eval("nonconservativeProduct", average_z, delta_z, x, h, t, dt, 2))
       output[0:NumberOfUnknowns, :, :, :] = output[0:NumberOfUnknowns, :, :, :] + 0.5 * NV * dt / h(2) * (-nonconservative_z[:, 0:NV, 0:NV, 0:NV] - nonconservative_z[:, 0:NV, 0:NV, 1:NV+1])
 
    if maxEigenvalue:
      maxEigenvaluesX2 = DataBlock([[0, 1], [0, NV], [0, NV], [0, NV], [0, N]], eval("maxEigenvalue", output, x, 0))
      maxEigenvaluesY2 = DataBlock([[0, 1], [0, NV], [0, NV], [0, NV], [0, N]], eval("maxEigenvalue", output, x, 1))
      maxEigenvaluesZ2 = DataBlock([[0, 1], [0, NV], [0, NV], [0, NV], [0, N]], eval("maxEigenvalue", output, x, 2))
      maxEigenvaluesXY = max(maxEigenvaluesX2, maxEigenvaluesY2)
      maxEigenvaluesPerVolume = max(maxEigenvaluesXY, maxEigenvaluesZ2)
      maxEigenvalues = max(maxEigenvaluesPerVolume)
 
 
class RusanovWithExplicitEuler(Solver.NumericalScheme):
    """!
 
    Finite Volumes with a Rusanov Riemann solver and explicit Euler time
    stepping, realised over a block-structured mesh.
 
    """
 
    def __init__(
        self,
        solver: Solver,
        patch_size: int,
        overlap: int = 1,
        cfl_constant: float = 0.9,
        interpolation: "toolbox.blockstructured.api.InterpolateHangingFaces.Scheme" = toolbox.blockstructured.api.InterpolateHangingFaces.Scheme.Linear,
        restriction: "toolbox.blockstructured.api.RestrictHangingFaces.Scheme" = toolbox.blockstructured.api.RestrictHangingFaces.Scheme.Averaging,
        min_volume_size: float = None,
        max_volume_size: float = None,
    ):
        """!
 
        The constructor memorises the attributes. However, it does not do this
        within this class, but it adds the important properties straight to the
        solver, so anybody else working against this solver knows about the
        properties introduced by the numerical scheme.
 
        @see @ref __create_data_structures()
 
        patch_size: int
          Number of Finite Volumes per axis within one patch (i.e. within
          one mesh cell).
 
        overlap: int
          Number of Finite Volumes the halo layer reconstructed around a
          patch - from the neighbouring patches' faces - is wide.
 
        cfl_constant: float
          Safety factor the admissible time step size (@see
          exahype2.api.fv.actionsets.ComputeGlobalCFLCondition,
          which this is passed straight through to) is scaled down by,
          relative to the theoretical CFL=1 stability limit
          h / (Dimensions * lambda_max). 0.9 is a
          conservative default for an explicit-Euler Finite Volume scheme;
          this codebase has no prior precedent value to match (grep for
          "cfl" outside this class finds only exaparty's unrelated SPH
          solvers, which default to 0.1).
 
        interpolation: toolbox.blockstructured.api.InterpolateHangingFaces.Scheme
          How the scheme initialises the halo of a fine patch next to a
          coarser patch on an adaptive mesh. Whenever the time step
          reconstructs the halo of a fine cell, it interpolates the data of
          the adjacent coarse face onto the fine cell's hanging face
          (toolbox.blockstructured.api.InterpolateHangingFaces). The
          default is linear interpolation.
 
        restriction: toolbox.blockstructured.api.RestrictHangingFaces.Scheme
          How the scheme hands the data of fine patches over to an adjacent
          coarser patch. Whenever the fine cells project their data onto
          their faces, the scheme restricts the data of the hanging faces
          onto the coarse face (toolbox.blockstructured.api.RestrictHangingFaces).
          The default is averaging.
 
        min_volume_size: float
          The time step asks for a refinement of a cell only if its volumes
          are at least of this size (@ref exahype2.api.blockstructured.DynamicMeshRefinement).
          None (the default) does not bound the refinement.
 
        max_volume_size: float
          The time step asks for an erase, i.e. coarsening, of a cell only
          if the volumes of the coarser cell would be strictly smaller than
          this size. None (the default) does not bound the coarsening.
 
        On a regular mesh, there are no hanging faces, and the two schemes
        are never used. The fluxes along a face between a coarse and a fine
        patch are not corrected, i.e. the scheme is not exactly conservative
        on adaptive meshes.
 
        """
        super().__init__(solver)
        assert (
            min_volume_size is None
            or max_volume_size is None
            or min_volume_size <= max_volume_size
        ), f"min_volume_size={min_volume_size} has to be smaller than or equal to max_volume_size={max_volume_size}"
        self.min_volume_size = min_volume_size
        self.max_volume_size = max_volume_size
        self.interpolation = interpolation
        self.restriction = restriction
        self.patch_size = patch_size
        self.overlap = overlap
        # Set in construct_output(): the generated kernels live in the project
        # namespace, as their names are solver-specific
        self.kernel_namespace = None
        self.cfl_constant = cfl_constant
        self.__create_data_structures()
 
    @property
    def kernel_function_name(self) -> str:
        """!
        Name of the generated time step kernel. It carries the solver name,
        as do the kernels/<SolverName>_rusanov*.h files holding it, so two
        solvers' kernels do not clash.
        """
        return self._solver_solver.name + "_timeStepWithRusanov"
 
    def __create_data_structures(self):
        """!
 
        Create the mesh data structures (a Finite Volume patch per cell,
        plus the face-based overlap patches used to reconstruct its halo)
        this numerical scheme requires.
 
        Ported from the legacy exahype2.solvers.fv.FV.create_data_structures():
        one @ref peano4.api.datamodel.Patch of patch_size^d Finite Volumes
        per mesh cell (_patch), plus a halo-augmented copy of it
        (_patch_with_halo, (patch_size+2*overlap)^d, @see
        toolbox.blockstructured.api.ReconstructPatchWithHalo) and two overlap
        patches of width 2*overlap attached to the mesh faces:
 
        - _patch_overlap_current_time_step: the persistent halo data
          toolbox.blockstructured.api.ReconstructPatchWithHalo reads from to
          rebuild _patch_with_halo (@see construct_output(), which hands it
          _patch/_patch_overlap_current_time_step/_patch_with_halo directly
          rather than the whole solver).
        - _patch_overlap_update: the projection of the just-updated _patch
          accumulated by the current sweep (@see
          exahype2.api.fv.actionsets.ProjectPatchOntoFaces), rolled over
          into _patch_overlap_current_time_step by
          exahype2.api.fv.actionsets.RollOverUpdatedFace once the update is
          complete.
        - _max_eigenvalue: one double per cell, named after the solver
          followed by MaxEigenvalue (e.g. EulerFVSmoothMaxEigenvalue),
          which holds the maximum eigenvalue of
          the patch after its last update
          (@see exahype2.api.fv.actionsets.UpdateCell).
          exahype2.api.fv.actionsets.ComputeGlobalCFLCondition reduces these
          values into the next time step size. It is modelled via DaStGen2
          as an array of cardinality one. The array is a
          peano4.api.dastgen2.Peano4SmartPointerDoubleArray, i.e. the value
          resides on the heap.
 
        All five are assigned onto self._solver (the shared
        exahype2.api.Solver instance both this numerical scheme and its
        surrounding solver's Solver.GridConstruction/Solver.Plotter
        counterparts hold a reference to as self._solver), not onto self,
        since the remaining FV action sets that need them - @see
        exahype2.api.fv.actionsets.ApplyDirichletBoundaryConditions,
        exahype2.api.fv.actionsets.UpdateCell,
        exahype2.api.fv.actionsets.ComputeGlobalCFLCondition - are handed that
        shared solver and look the patches up as solver._patch /
        solver._patch_overlap_*.
 
        Storage is hardwired to smart pointers (peano4.api.datamodel's
        PatchToDoubleArrayWithSmartPointer), matching the legacy default
        (Storage.SmartPointers for both cell and face data) - the legacy
        Storage enum's other two backends (CallStack, Heap) are not
        exposed here.
 
        Two pieces of the legacy method are deliberately not ported yet:
        - The cell/face label creation
          (exahype2.grid.create_cell_label()/create_face_label()) - that
          module has not been ported to the new API.
        - The load/store/compute guards
          (FV._provide_cell_data_to_compute_kernels_default_guard() and
          friends) - they are keyed off the legacy SolverState state
          machine, which the task-graph architecture has no equivalent of:
          a task only ever runs within its own phase's task graph in the
          first place, so "are we still constructing the grid" does not
          need to be re-checked at runtime the same way.
 
        """
        total_unknowns = self._solver_solver.unknowns + self._solver_solver.auxiliary_variables
        identifier = self._solver_solver._unknown_identifier()
 
        self._solver_solver._patch = peano4.api.datamodel.Patch(
            (self.patch_size, self.patch_size, self.patch_size),
            total_unknowns,
            identifier,
        )
        self._solver_solver._patch_with_halo = peano4.api.datamodel.Patch(
            (self.patch_size+2*self.overlap, self.patch_size+2*self.overlap, self.patch_size+2*self.overlap),
            total_unknowns,
            identifier + "WithHalo",
        )
        self._solver_solver._patch_overlap_current_time_step = peano4.api.datamodel.Patch(
            (2 * self.overlap, self.patch_size, self.patch_size),
            total_unknowns,
            identifier + "CurrentTimeStep",
        )
        self._solver_solver._patch_overlap_update = peano4.api.datamodel.Patch(
            (2 * self.overlap, self.patch_size, self.patch_size),
            total_unknowns,
            identifier + "Update",
        )
        self._solver_solver._max_eigenvalue = peano4.api.datamodel.DaStGen2(self._solver_solver.name + "MaxEigenvalue")
        self._solver_solver._max_eigenvalue.data.add_attribute(
            peano4.api.dastgen2.Peano4SmartPointerDoubleArray("maxEigenvalue", 1)
        )
 
 
    def __add_data_to_repository(self, peano_project):
        """!
        
        Consult @ref page_exahype_programmers_guide_new_solver for some information
        how this routine fits to the typical workflow when we introduce a new
        solver.
        
        """
        peano4.api.Logger().log("""
Patch data
----------
{}
Patch overlap data
------------------
{}
{}""".format( str(self._solver_solver._patch), str(self._solver_solver._patch_overlap_current_time_step), str(self._solver_solver._patch_overlap_update) ) )
 
        peano_project.datamodel.add_cell(self._solver_solver._patch)
        peano_project.datamodel.add_cell(self._solver_solver._patch_with_halo)
        peano_project.datamodel.add_face(self._solver_solver._patch_overlap_current_time_step)
        peano_project.datamodel.add_face(self._solver_solver._patch_overlap_update)
        peano_project.datamodel.add_cell(self._solver_solver._max_eigenvalue)
 
        # The task graph compiler only uses, loads and stores data where the
        # solver's grid selector holds, too (@see peano4.api.datamodel.DoF).
        for patch in (
            self._solver_solver._patch,
            self._solver_solver._patch_with_halo,
            self._solver_solver._patch_overlap_current_time_step,
            self._solver_solver._patch_overlap_update,
            self._solver_solver._max_eigenvalue,
        ):
            patch.grid_selector = self._solver_solver.grid_selector
 
        self._solver_solver._patch.generator =peano4.api.datamodel.PatchToDoubleArrayWithSmartPointer(
            self._solver_solver._patch, "double"
        )
        self._solver_solver._patch_with_halo.generator = peano4.api.datamodel.PatchToDoubleArrayWithSmartPointer(
            self._solver_solver._patch_with_halo, "double"
        )
        self._solver_solver._patch_overlap_current_time_step.generator = peano4.api.datamodel.PatchToDoubleArrayWithSmartPointer(
            self._solver_solver._patch_overlap_current_time_step, "double"
        )
        self._solver_solver._patch_overlap_update.generator = peano4.api.datamodel.PatchToDoubleArrayWithSmartPointer(
            self._solver_solver._patch_overlap_update, "double"
        )
        
        self._solver_solver._patch_overlap_current_time_step.generator.merge_method_definition = (
            toolbox.blockstructured.api.get_face_merge_implementation(
                self._solver_solver._patch_overlap_current_time_step
            )
        )
        self._solver_solver._patch_overlap_current_time_step.generator.includes += """
#include "peano4/utils/Loop.h"
"""
 
        # Only faces selected by the solver's grid selector hold valid data,
        # so only those are exchanged.
        self._solver_solver._patch_overlap_current_time_step.generator.send_condition = self._solver_solver.grid_selector.createFaceGuard()
        self._solver_solver._patch_overlap_current_time_step.generator.receive_and_merge_condition = self._solver_solver.grid_selector.createFaceGuard()
 
        self._solver_solver._patch_overlap_update.generator.send_condition = "false"
        self._solver_solver._patch_overlap_update.generator.receive_and_merge_condition = "false"
 
    def create_default_mesh_adaptation(self) -> Solver.MeshAdaptation:
        """!
 
        The Finite Volume schemes adapt the mesh and balance the load
        through exahype2.api.fv.AdaptMeshAndBalanceLoad by default
 
        """
        from .AdaptMeshAndBalanceLoad import AdaptMeshAndBalanceLoad
 
        return AdaptMeshAndBalanceLoad(self._solver_solver)
 
    def construct_output(self,
                         exahype2_project,
                         peano_project,
                         time_step):
        """!
 
        """
        peano4.api.Logger().log("Create time stepping task graph for Finite Volume scheme over block-structured meshes with Rusanov flux")
 
        self.__add_data_to_repository(peano_project)
 
        self.kernel_namespace = peano_project.namespace
        # Construct the call first, as it validates the PDE terms' variants
        self._solver_solver.compute_kernel_call = self.__construct_kernel_call(Target.Serial, False)
        self.__create_time_step_functions(peano_project.directory)
        
        reconstruction_task         = peano4.api.taskgraph.Task(name=f"{self._solver.name}_Reconstruct", entity=peano4.api.taskgraph.Task.GridEntity.CELL)
        update_cell_task            = peano4.api.taskgraph.Task(name=f"{self._solver.name}_Update cell", entity=peano4.api.taskgraph.Task.GridEntity.CELL)
        project_on_faces_task       = peano4.api.taskgraph.Task(name=f"{self._solver.name}_Project on faces", entity=peano4.api.taskgraph.Task.GridEntity.CELL)
        compute_cfl_condition_task  = peano4.api.taskgraph.Task(name=f"{self._solver.name}_Compute global CFL condition", entity=peano4.api.taskgraph.Task.GridEntity.CELL, may_spawn_background_tasks=False)
        
        time_step.add_task(reconstruction_task)
        time_step.add_task(update_cell_task)
        time_step.add_task(project_on_faces_task)
        time_step.add_task(compute_cfl_condition_task)
 
        # The load balancing is switched off throughout the time step, and
        # switched on again only at its very end, so it may then decide on
        # splits. Peano carries these splits out in the subsequent mesh sweep.
        reconstruction_task.add_action_set(
            toolbox.loadbalancing.api.actionsets.ToggleLoadBalancing(
                enable_load_balancing=False, plug_into_prepare_traversal=True
            )
        )
 
        # The very last task evaluates the adaptivity criterion if the
        # solver has one (create_dynamic_mesh_refinement_action_set()), i.e. it
        # collects refine and erase requests for the mesh adaptation phase,
        # and it switches the load balancing on again and informs it that
        # the time step has finished.
        dynamic_mesh_refinement_task = peano4.api.taskgraph.Task(
            name=f"{self._solver.name}_Dynamic mesh refinement",
            entity=peano4.api.taskgraph.Task.GridEntity.CELL,
            may_spawn_background_tasks=False,
        )
        dynamic_mesh_refinement_task.use_cell(self._solver_solver._patch, peano4.api.taskgraph.Task.DataAccess.READ_AND_WRITE)
        dynamic_mesh_refinement = create_dynamic_mesh_refinement_action_set(
            self._solver_solver,
            self.patch_size,
            guard=self._solver_solver.grid_selector.createCellGuard(),
            min_volume_size=self.min_volume_size,
            max_volume_size=self.max_volume_size,
        )
        if dynamic_mesh_refinement is not None:
            dynamic_mesh_refinement_task.add_action_set(dynamic_mesh_refinement)
        dynamic_mesh_refinement_task.add_action_set(
            toolbox.loadbalancing.api.actionsets.LoadBalancing(enable_load_balancing=True)
        )
        dynamic_mesh_refinement_task.add_dependency(compute_cfl_condition_task)
        dynamic_mesh_refinement_task.add_dependency(project_on_faces_task)
        time_step.add_task(dynamic_mesh_refinement_task)
 
        
        reconstruction_task.use_cell(self._solver_solver._patch, peano4.api.taskgraph.Task.DataAccess.READ_AND_WRITE)
        reconstruction_task.use_cell(self._solver_solver._patch_with_halo, peano4.api.taskgraph.Task.DataAccess.WRITE)
        
        reconstruction_task.use_face(
            self._solver_solver._patch_overlap_current_time_step,
            get_access_to_persistent_face_data(self._solver_solver, exahype2_project, peano4.api.taskgraph.Task.DataAccess.READ_AND_DISCARD),
        )
        reconstruction_task.add_action_set(
            toolbox.blockstructured.api.ReconstructPatchWithHalo(
                patch=self._solver_solver._patch,
                patch_overlap=self._solver_solver._patch_overlap_current_time_step,
                reconstructed_patch=self._solver_solver._patch_with_halo,
                guard=self._solver_solver.grid_selector.createCellGuard(),
                add_assertions_to_halo_exchange=False,
                additional_includes=self._solver_solver._get_default_includes() + self._solver_solver.user_solver_includes,
            )
        )
        
        reconstruction_task.add_action_set(
            toolbox.blockstructured.api.InterpolateHangingFaces(
                self._solver_solver._patch_overlap_current_time_step,
                scheme=self.interpolation,
                additional_includes=self._solver_solver._get_default_includes() + self._solver_solver.user_solver_includes,
            )
        )
        
        ApplyDirichletBoundaryConditions.add_task_to_time_step(self._solver_solver, exahype2_project, time_step, reconstruction_task)
 
        
        update_cell_task.add_dependency(reconstruction_task)
        update_cell_task.use_cell(self._solver_solver._patch_with_halo, peano4.api.taskgraph.Task.DataAccess.READ_AND_DISCARD)
        update_cell_task.use_cell(self._solver_solver._patch,           peano4.api.taskgraph.Task.DataAccess.WRITE)
        update_cell_task.use_cell(self._solver_solver._max_eigenvalue,  peano4.api.taskgraph.Task.DataAccess.WRITE)
        update_cell_task.add_action_set(
            UpdateCell(
                self._solver_solver,
                guard=self._solver_solver.grid_selector.createCellGuard(),
            )
        )
        update_cell_task.tags.add(Solver.SolutionUpdateTaskTag)
 
        
        compute_cfl_condition_task.add_dependency(update_cell_task)
        compute_cfl_condition_task.use_cell(self._solver_solver._max_eigenvalue, peano4.api.taskgraph.Task.DataAccess.READ_AND_DISCARD)
        compute_cfl_condition_task.add_action_set(
            ComputeGlobalCFLCondition(
                self._solver_solver,
                self.cfl_constant,
                self.patch_size,
                guard=self._solver_solver.grid_selector.createCellGuard(),
            )
        )
        
        
        project_on_faces_task.add_dependency(update_cell_task)
        project_on_faces_task.use_cell(self._solver_solver._patch, peano4.api.taskgraph.Task.DataAccess.READ_AND_WRITE)
        project_on_faces_task.use_face(self._solver_solver._patch_overlap_update, peano4.api.taskgraph.Task.DataAccess.TEMPORARY)
        project_on_faces_task.use_face(
            self._solver_solver._patch_overlap_current_time_step,
            get_access_to_persistent_face_data(self._solver_solver, exahype2_project, peano4.api.taskgraph.Task.DataAccess.WRITE),
        )
        project_on_faces_task.add_action_set(
            ProjectPatchOntoFaces(
                self._solver_solver._patch,
                self._solver_solver._patch_overlap_update,
                guard=self._solver_solver.grid_selector.createCellGuard(),
                includes=self._solver_solver._get_default_includes() + self._solver_solver.user_solver_includes,
            )
        )
        
        project_on_faces_task.add_action_set(self.create_restriction_of_hanging_faces())
        
        project_on_faces_task.add_action_set(
            RollOverUpdatedFaceWithSommerfeldBoundaryConditions.create_roll_over(
                self._solver_solver,
                exahype2_project,
                guard=self._solver_solver.grid_selector.createFaceGuard(),
                includes=self._solver_solver._get_default_includes() + self._solver_solver.user_solver_includes,
                time_step_size=f"repositories::{self._solver.get_name_of_global_instance()}.getMinTimeStepSize()",
            )
        )
        
        if peano4.api.Logger().is_visualisation:
            time_step.visualise(title="Time stepping task graph for {}".format(self._solver_solver.name))
        if peano4.api.Logger().is_text:
            time_step.log_as_text(title="Time stepping task graph for {}".format(self._solver_solver.name))
 
    def create_restriction_of_hanging_faces(self) -> "toolbox.blockstructured.api.RestrictHangingFaces":
        """!
 
        Action set which restricts the projected fine data from the hanging
        faces onto the coarse faces
 
        Every task which projects the solution onto the faces has to host
        this action set, so the faces between coarse and fine patches hold
        consistent data afterwards. Besides the time step, this is the
        projection at the end of the setup
        (@ref exahype2.api.fv.SetupSimulationOverRegularGrid). The task has
        to write all of _patch_overlap_update, as averaging clears it first.
 
        """
        return toolbox.blockstructured.api.RestrictHangingFaces(
            self._solver_solver._patch_overlap_update,
            scheme=self.restriction,
            clear_guard=self._solver_solver.grid_selector.createFaceGuard(),
            additional_includes=self._solver_solver._get_default_includes() + self._solver_solver.user_solver_includes,
        )
 
    def __create_time_step_functions(self, output_directory):
        """!
        
        Calling this function is delayed until all PDE term information is provided.
        In this way we can avoid generating code that will not be used. The following
        information is required to translate a Python function using DSLHyPE:
        - template_parameters: the variables that are used within the Python function
          to represent external information
        - functor_parameters: the PDE terms that can be used by the solver
        - function: the Python function itself
 
        TODO: Currently the MLIR types are passed manually here, this should ideally
        be inferred
 
        """
 
        template_parameters = [
            dslhype.SyntaxTree.Argument("NV", dslhype.SyntaxTree.TInteger()),
            dslhype.SyntaxTree.Argument("H", dslhype.SyntaxTree.TInteger()),
            dslhype.SyntaxTree.Argument("NumberOfUnknowns", dslhype.SyntaxTree.TInteger()),
            dslhype.SyntaxTree.Argument("NumberOfAuxiliaryVariables", dslhype.SyntaxTree.TInteger())
        ]
    
        functor_parameters = [
            dslhype.SyntaxTree.Argument("flux", dslhype.SyntaxTree.TCustom("const ::exahype2::fv::rusanov::FluxFunctor&"), "(!llvm.ptr, !llvm.ptr, !llvm.ptr, f64, f64, i32, !llvm.ptr) -> ()"),
            dslhype.SyntaxTree.Argument("nonconservativeProduct", dslhype.SyntaxTree.TCustom("const ::exahype2::fv::rusanov::NonconservativeProductFunctor&"), "(!llvm.ptr, !llvm.ptr, !llvm.ptr, !llvm.ptr, f64, f64, i32, !llvm.ptr) -> ()"),
            dslhype.SyntaxTree.Argument("volumetricSource", dslhype.SyntaxTree.TCustom("const ::exahype2::fv::rusanov::VolumetricSourceFunctor&"), "(!llvm.ptr, !llvm.ptr, f64, !llvm.ptr) -> ()"),
            dslhype.SyntaxTree.Argument("maxEigenvalue", dslhype.SyntaxTree.TCustom("const ::exahype2::fv::rusanov::MaxEigenvalueFunctor&"), "(!llvm.ptr, !llvm.ptr, !llvm.ptr, f64, f64, i32, !llvm.ptr) -> ()")
        ]
 
        time_step_dsl_function_2d = timeStepWithRusanov2D
        time_step_dsl_function_3d = timeStepWithRusanov3D
        dslhype.translate_dslhype_function(
            function = time_step_dsl_function_2d,
            function_namespace = self.kernel_namespace,
            template_parameters=template_parameters,
            functor_parameters=functor_parameters,
            pde_terms=self._solver_solver.pde_terms,
            output_file_name=self._solver_solver.name + "_rusanov2d",
            target=dslhype.Target.Serial,
            stateless=False,
            function_name=self.kernel_function_name,
            output_directory=output_directory,
        )
        dslhype.translate_dslhype_function(
            function = time_step_dsl_function_3d,
            function_namespace = self.kernel_namespace,
            template_parameters=template_parameters,
            functor_parameters=functor_parameters,
            pde_terms=self._solver_solver.pde_terms,
            output_file_name=self._solver_solver.name + "_rusanov3d",
            target=dslhype.Target.Serial,
            stateless=False,
            function_name=self.kernel_function_name,
            output_directory=output_directory,
        )
        file = open(os.path.join(output_directory, "kernels", f"{self._solver.name}_rusanov.h"), "w")
        file.write(f"""#if Dimensions == 2
    #include "{self._solver.name}_rusanov2d.h"
#else
    #include "{self._solver.name}_rusanov3d.h"
#endif""")
        file.close()
 
    def __construct_kernel_call(self, target: Target, stateless: bool):
        """!
 
        This function constructs the calls to the compute kernels based on the
        application PDE terms and the version chosen.
 
        """
        available_terms = ["flux", "nonconservativeProduct", "volumetricSource", "maxEigenvalue"]
 
        if stateless:
            return f"""::{"::".join(self.kernel_namespace)}::{self.kernel_function_name}<
        {self.patch_size},
        {self.overlap},
        {self._solver.unknowns},
        {self._solver.auxiliary_variables}
    >(
        patchData,
        {self._solver.name}
    );"""
        else:
            # Variants per term the kernel takes: {name: {stateless: PDETerm}}
            variants = {}
            for application_term in self._solver_solver.pde_terms:
                if application_term.NAME in available_terms:
                    variants.setdefault(application_term.NAME, {})[application_term.stateless] = application_term
 
            terms_with_both_variants = [name for name in variants if len(variants[name]) == 2]
            terms_without_stateless_variant = [name for name in variants if True not in variants[name]]
            if terms_with_both_variants and terms_without_stateless_variant:
                raise Exception(
                    f"Solver {self._solver.name}: the PDE terms {terms_with_both_variants} are "
                    f"registered both with and without state. In this case, the generated code "
                    f"asks cellCanUseStatelessPDETerms() per cell before it invokes the kernel: "
                    f"if it returns true, the kernel is called with the stateless variants of all "
                    f"PDE terms, otherwise with the variants with state. However, the PDE terms "
                    f"{terms_without_stateless_variant} have no stateless variant, so a kernel "
                    f"call would mix stateless and stateful routines. Either add stateless "
                    f"variants of {terms_without_stateless_variant}, or register "
                    f"{terms_with_both_variants} only once."
                )
 
            def kernel_call(stateless: bool) -> str:
                pde_term_calls = []
                for term in available_terms:
                    if term not in variants:
                        continue
                    # A term registered only once is used in either case
                    application_term = variants[term].get(stateless, next(iter(variants[term].values())))
                    # Commas inside template brackets, e.g. Vector<Dimensions, double>,
                    # do not separate parameters
                    arguments = [
                        parameter.split()[-1]
                        for parameter in re.split(r",(?![^<]*>)", application_term.SIGNATURE)
                    ]
                    # Stateless terms are static and take a dummy tag (@see PDETerm)
                    if application_term.stateless:
                        arguments.append("::exahype2::Solver::Offloadable::Yes")
                    pde_term_calls.append(f"""
    [&](
        {application_term.SIGNATURE}
    )->void {{
        repositories::{self._solver.get_name_of_global_instance()}.{application_term.name}({", ".join(arguments)});
    }}""")
                return f"""::{"::".join(self.kernel_namespace)}::{self.kernel_function_name}<
    {self.patch_size},
    {self.overlap},
    {self._solver.unknowns},
    {self._solver.auxiliary_variables}
>(
    patchData,
{",\n".join(pde_term_calls)}
);"""
 
            if not terms_with_both_variants:
                return kernel_call(stateless=False)
            return f"""if (repositories::{self._solver.get_name_of_global_instance()}.cellCanUseStatelessPDETerms(
    patchData.cellCentre[0],
    patchData.cellSize[0],
    patchData.t[0],
    patchData.dt[0]
)) {{
{kernel_call(stateless=True)}
}} else {{
{kernel_call(stateless=False)}
}}"""