ADERDG.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
from enum import IntEnum
from abc import abstractmethod
 
import jinja2
 
import peano4.output
import peano4.datamodel
import peano4.toolbox.blockstructured
 
import exahype2.grid
import exahype2.solvers.aderdg.kernels
import exahype2.solvers.aderdg.actionsets
 
from exahype2.solvers.Storage import Storage
from exahype2.solvers.PDETerms import PDETerms
from exahype2.solvers.aderdg.Architectures import Architectures
from exahype2.solvers.aderdg.LagrangeBasis import GaussLegendreBasis, GaussLobattoBasis
 
 
"""
Used to have a central way for all of this class' descendants to differentiate between Gauss-Legendre
and Gauss-Lobatto nodes, helps to avoid having to constantly pass a string and ensure that this string is
the same everywhere.
"""
class Polynomials(IntEnum):
    Gauss_Legendre = 0,
    Gauss_Lobatto = 1
 
"""
Denotes valid precisions for kernels and storage but also is supposed to allow a little more flexibility
in the naming. If a user wants to name it float because they're more used to C-based languages that's
fine but if they'd rather use single we allow it as well.
"""
PrecisionType = {
    "SP":     "float",
    "float":  "float",
    "single": "float",
    "fp32":   "float",
    "DP":     "double",
    "double": "double",
    "fp64":   "double",
    "quadruple":   "long double",
    "long double": "long double",
    "fp128":       "long double",
    "float16_t":        "std::float16_t",
    "std::float16_t":   "std::float16_t",
    "_Float16":         "_Float16",
    "fp16":             "_Float16",
    "bfloat16_t":       "std::bfloat16_t",
    "std::bfloat16_t":  "std::bfloat16_t",
    "__bf16":           "__bf16",
    "bfloat":           "__bf16",
    "bf16":             "__bf16"
}
 
class ADERDG(object):
    def __init__(
        self,
        name,
        order,
        unknowns,
        auxiliary_variables,
        min_cell_h,
        max_cell_h,
        plot_grid_properties=False,
    ):
        self._name = name
 
        """
        Polynomial order of the representation of the underlying function in each spatial dimension.
        The number of support points used to represent this polynomial is always the order + 1
        """
        self._order = order
 
        """
        self._unknowns and self._auxiliary_variables respectively hold the number of unknowns and
        auxiliary variables in the equation to be computed. Unknowns are variables that change over
        time whereas auxiliary variables can be space-dependent but don't vary over time.
        These can be specified either as simple ints or by a dictionary
        (e.g.) unknowns = {'a': 1, 'b': 1, 'c': 3}
        in which the user specifies the multiplicity of the variable (the velocity has one component
        per dimension for example.)
        If they are specified by a dictionary then the code generates a "VariableShortcuts" file which
        allows the user to specify a variable by name and automatically maps this to the right position
        in an array for better legibility. Otherwise they must manually remember the position of each
        variable manually.
 
        use_var_shortcut is used to know whether or not the user passed their variables via a dict
        variable_names and variable_pos are used internally to remember the names and respective
        positions of variables if set by a dictionary.
        """
        self._use_var_shortcut = False
        self._variable_names = []
        self._variable_pos = [0]
 
        if type(unknowns) is dict:
            self._unknowns = sum(unknowns.values())
            self._variable_names += list(unknowns.keys())
            for var in list(unknowns.values()):
                self._variable_pos.append(var + self._variable_pos[-1])
            self._use_var_shortcut = True
        elif type(unknowns) is int:
            self._unknowns = unknowns
        else:
            raise Exception(
                "not a valid type for parameter unknowns, needs to be int or dictionary."
            )
 
        if type(auxiliary_variables) is dict:
            self._auxiliary_variables = sum(auxiliary_variables.values())
            self._variable_names += list(auxiliary_variables.keys())
            for var in list(auxiliary_variables.values()):
                self._variable_pos.append(var + self._variable_pos[-1])
            self._use_var_shortcut = True
        elif type(auxiliary_variables) is int:
            self._auxiliary_variables = auxiliary_variables
        else:
            raise Exception(
                "not a valid type for parameter auxiliary_variables, needs to be int or dictionary."
            )
 
        """
        Minimal and maximal accepted sizes for the cell, domain will be subdivided until the cell size is smaller
        than the minimal size but will never be subdivided to be smaller than the minimal size.
        """
        self._min_cell_h = min_cell_h
        self._max_cell_h = max_cell_h
        if min_cell_h > max_cell_h:
            raise Exception(
                "min_cell_h ("
                + str(min_cell_h)
                + ") is bigger than max_cell_h ("
                + str(max_cell_h)
                + ")"
            )
 
        self._is_linear     = False
        self._polynomials   = Polynomials.Gauss_Legendre
        self._basis         = GaussLegendreBasis(order + 1)
        self._architecture  = "noarch"
        self._use_libxsmm   = False
 
        self._initialise_patches = False
 
        self._predictor_computation_precisions      = [PrecisionType["double"]]
        self._corrector_computation_precision       = PrecisionType["double"]
        self._solution_persistent_storage_precision = PrecisionType["double"]
        self._precompute_picard_precision           = False
        self._use_half_precision                    = False
 
        # Whether or not to plot information about the grid as well as information about the content of the patches
        self._plot_grid_properties = plot_grid_properties
        self.plot_description = ""
 
        self.plot_metadata = ""
        self.select_dofs_to_print = None
 
        # Used by children classes or setters in order to add constants or includes to the corresponding part
        self._solver_constants = ""
        self.user_action_set_includes = ""
        self.user_solver_includes = ""
 
        """
        These specify the time step, they must be set by child classes in order to advance the timestepping.
        These are the main difference between child classes as they all use the same underlying routines.
        """
        self._compute_time_step_size = "#error Not yet defined"
        self._compute_new_time_step_size = "#error Not yet defined"
 
        """
        These are all the optional and mandatory plugins that can be specified by the user.
        eigenvalues, boundary conditions, refinement criteria and initial conditions must
        all be specified, all others are optional but can be added via the set_implementation() function.
        There are four choices for these, three of which are specified in PDETerms:
        User_Defined_Implementation creates empty stubs that the user must then manually fill.
        None_implementation leaves out this optional stub and ensures the kernels don't use this function.
        Empty_implementation does... the same thing as user_defined I think?
        The final option is to fill these with a string, in which case this string is directly taken as
        C++ code and used as the implementation of this function. The user is trusted to provide code that
        compiles and runs properly, it is never checked.
        """
        self._flux_implementation           = PDETerms.None_Implementation
        self._ncp_implementation            = PDETerms.None_Implementation
        self._source_term_implementation    = PDETerms.None_Implementation
        self._point_sources_implementation  = PDETerms.None_Implementation
        self._eigenvalues_implementation    = PDETerms.None_Implementation
        self._material_param_implementation = PDETerms.None_Implementation
 
        self._refinement_criterion_implementation = PDETerms.None_Implementation
        self._initial_conditions_implementation   = PDETerms.None_Implementation
        self._boundary_conditions_implementation  = PDETerms.None_Implementation
 
        self._precision_criterion_implementation  = PDETerms.None_Implementation
        self._riemann_solver_implementation       = PDETerms.None_Implementation
        self._number_of_point_sources = 0
 
        """
        These will be used by children of this class to fill them with the declaration and definition of the
        above mandatory and optional functions.
        """
        self._abstract_solver_user_declarations = ""
        self._abstract_solver_user_definitions = ""
        self._solver_user_declarations = ""
        self._solver_user_definitions = ""
 
        """
        Used by child classes; will contain code to be executed at the start and end respectively of the
        time step. Typically specifies some timestep variable and advances the overall timestamp of the
        entire solver.
        """
        self._start_time_step_implementation  = ""
        self._finish_time_step_implementation = ""
 
        self._constructor_implementation = ""
        self._destructor_implementation = ""
 
        """
        These are ADER-DG CFL factors, they help determine the maximal stable timestep size.
        Up to order 4, computed via von Neumann analysis [1].
        [1] M. Dumbser, D. S. Balsara, E. F. Toro, and C.-D. Munz,
        ‘A unified framework for the construction of one-step finite volume and discontinuous Galerkin schemes on unstructured meshes’, Journal of Computational Physics, vol. 227, no. 18, pp. 8209–8253, Sep. 2008.
        They are necessary to compute stable timestep sizes as dt <= c*dx / (lambda_max*(2N+1)
        with c the CFL for given order, lambda_max the maximal eigenvalue, N the order and dx the minimal mesh distance.
        """
        self._cflAder = [
            1.0,
            0.33,
            0.17,
            0.1,
            0.069,
            0.045,
            0.038,
            0.03,
            0.02,
            0.015,
            0.008821428571428572,
            0.005235326086956522,
            0.0031307249999999996,
            0.0018842326388888888,
            0.001140285614224138,
            0.0006933672202620968,
        ]
 
        # Safety factor for CFL (Courant-Friedrichs-Lewy) condition
        self._cflSafetyFactor = 0.0
 
        """
        Used if the user wants to add some function to postprocessing, meaning
        it should be added after every other computation in ADER-DG is done.
        """
        self._action_set_postprocess_solution = None
 
        self._preprocess_reconstructed_patch = ""
        self._postprocess_updated_patch = ""
 
        self.switch_storage_scheme(Storage.SmartPointers, Storage.SmartPointers)
 
        # Create initial structures and action_sets, these are re-generated essentially anytime something about the solver is modified.
        self.create_data_structures()
        self.create_action_sets()
 
 
    def __str__(self):
        result = (
            """
Name:                   """
            + self._name
            + """
Type:                   """
            + self.__class__.__module__
            + """
Storage schemes:        """
            + str(self._cell_data_storage) + ", " + str(self._face_data_storage)
            + """
DG order:               """
            + str(self._order)
            + """
Polynomials:            """
            + str(self._polynomials.name)
            + """
Basis:                  """
            + str(self._basis)
            + """
Basis size:             """
            + str(self._order + 1)
            + """
Linear:                 """
            + str(self._is_linear)
            + """
Unknowns:               """
            + str(self._unknowns)
            + """
Auxiliary variables:    """
            + str(self._auxiliary_variables)
            + """
h_cell_min:             """
            + str(self._min_cell_h)
            + """
h_cell_max:             """
            + str(self._max_cell_h)
            + """
Initial conditions:     """
            + str(self._initial_conditions_implementation)
            + """
Boundary conditions:    """
            + str(self._boundary_conditions_implementation)
            + """
Refinement criterion:   """
            + str(self._refinement_criterion_implementation)
            + """
Eigenvalues:            """
            + str(self._eigenvalues_implementation)
            + """
Flux:                   """
            + str(self._flux_implementation)
            + """
NCP:                    """
            + str(self._ncp_implementation)
            + """
Source term:            """
            + str(self._source_term_implementation)
            + """
Point sources:          """
            + str(self._point_sources_implementation)
            + """
Number of point sources:"""
            + str(self._number_of_point_sources)
            + """
Riemann solver:         """
            + str(self._riemann_solver_implementation)
            + """
Material parameters:    """
            + str(self._material_param_implementation)
            + """
"""
        )
        return result
 
    __repr__ = __str__
 
    """
    Creates all the structures that are attached on the Peano grid either in a cell or on the faces.
    These are defined as Peano patches which have given sizes to contain persistent data as well as
    conditions for when they should be merged, sent, received, stored, loaded, etc.
    """
    @abstractmethod
    def create_data_structures(self):
        self._number_of_face_projections = 1
 
        # Current time step contains the solution, projections are respectively the projection
        # of the solution and of the fluxes on the faces as predicted by the space-time predictor.
        self._current_time_step = peano4.datamodel.Patch(
            (self._order + 1, self._order + 1, self._order + 1),
            self._unknowns + self._auxiliary_variables,
            self._unknown_identifier(),
            self._solution_persistent_storage_precision,
        )
        self._previous_time_step = peano4.datamodel.Patch(
            (self._order + 1, self._order + 1, self._order + 1),
            self._unknowns + self._auxiliary_variables,
            self._unknown_identifier() + "_old",
            self._solution_persistent_storage_precision,
        )
        self._rhs_estimates_projection = peano4.datamodel.Patch(
            (2, self._order+1, self._order+1),
            self._unknowns + self._auxiliary_variables,
            self._unknown_identifier() + "Estimates",
            self._corrector_computation_precision,
        )
        self._flux_estimates_projection = peano4.datamodel.Patch(
            (2, self._order+1, self._order+1),
            self._unknowns,
            self._unknown_identifier() + "FluxEstimates",
            self._corrector_computation_precision,
        )
 
        if self._cell_data_storage == Storage.CallStack:
          self._current_time_step.generator = peano4.datamodel.PatchToDoubleArray(
              self._current_time_step,
              self._solution_persistent_storage_precision
          )
          self._previous_time_step.generator = peano4.datamodel.PatchToDoubleArray(
              self._previous_time_step,
              self._solution_persistent_storage_precision
          )
        elif self._cell_data_storage == Storage.Heap:
          self._current_time_step.generator = peano4.datamodel.PatchToDoubleArrayOnHeap(
              self._current_time_step,
              self._solution_persistent_storage_precision
          )
          self._previous_time_step.generator = peano4.datamodel.PatchToDoubleArrayOnHeap(
              self._previous_time_step,
              self._solution_persistent_storage_precision
          )
        elif self._cell_data_storage == Storage.SmartPointers:
          self._current_time_step.generator = peano4.datamodel.PatchToDoubleArrayWithSmartPointer(
              self._current_time_step,
              self._solution_persistent_storage_precision
          )
          self._previous_time_step.generator = peano4.datamodel.PatchToDoubleArrayWithSmartPointer(
              self._previous_time_step,
              self._solution_persistent_storage_precision
          )
        else:
            assert False, "Storage variant {} not supported".format(
                self._cell_data_storage
            )
 
        if self._face_data_storage == Storage.CallStack:
          self._rhs_estimates_projection.generator = peano4.datamodel.PatchToDoubleArray(
              self._rhs_estimates_projection,
              self._corrector_computation_precision
          )
          self._flux_estimates_projection.generator = peano4.datamodel.PatchToDoubleArray(
              self._flux_estimates_projection,
              self._corrector_computation_precision
          )
        elif self._face_data_storage == Storage.Heap:
          self._rhs_estimates_projection.generator = peano4.datamodel.PatchToDoubleArrayOnHeap(
              self._rhs_estimates_projection,
              self._corrector_computation_precision
          )
          self._flux_estimates_projection.generator = peano4.datamodel.PatchToDoubleArrayOnHeap(
              self._flux_estimates_projection,
              self._corrector_computation_precision
          )
        elif self._face_data_storage == Storage.SmartPointers:
          self._rhs_estimates_projection.generator = peano4.datamodel.PatchToDoubleArrayWithSmartPointer(
              self._rhs_estimates_projection,
              self._corrector_computation_precision
          )
          self._flux_estimates_projection.generator = peano4.datamodel.PatchToDoubleArrayWithSmartPointer(
              self._flux_estimates_projection,
              self._corrector_computation_precision
          )
        else:
            assert False, "Storage variant {} not supported".format(
                self._face_data_storage
            )
 
        # Custom merge implementation takes into account that the faces project such that
        # the left cell projects to the left half of the face and the right cell projects to the right half.
        self._rhs_estimates_projection.generator.merge_method_definition = (
            get_face_overlap_merge_implementation(self._rhs_estimates_projection)
        )
        self._rhs_estimates_projection.generator.includes += """
#include "peano4/utils/Loop.h"
#include "repositories/SolverRepository.h"
"""
 
        self._flux_estimates_projection.generator.merge_method_definition = (
            get_face_overlap_merge_implementation(self._flux_estimates_projection)
        )
        self._flux_estimates_projection.generator.includes += """
#include "peano4/utils/Loop.h"
#include "repositories/SolverRepository.h"
"""
 
        self._current_time_step.generator.load_store_compute_flag = "::peano4::grid::constructLoadStoreComputeFlag({},{},{})".format(
            self._provide_cell_data_to_compute_kernels_default_guard(),
            self._load_cell_data_default_guard(),
            self._store_cell_data_default_guard(),
        )
 
        self._previous_time_step.generator.load_store_compute_flag = "::peano4::grid::constructLoadStoreComputeFlag({},{},{})".format(
            self._provide_cell_data_to_compute_kernels_default_guard(),
            self._load_cell_data_default_guard(),
            self._store_cell_data_default_guard(),
        )
 
        self._rhs_estimates_projection.generator.load_store_compute_flag = "::peano4::grid::constructLoadStoreComputeFlag({},{},{})".format(
            self._provide_cell_data_to_compute_kernels_default_guard(),
            self._load_face_data_default_guard(),
            self._store_face_data_default_guard()
        )
 
        # Gets overwritten in subclasses to take into account when to send or receive faces
        self._rhs_estimates_projection.generator.send_condition = "false"
        self._rhs_estimates_projection.generator.receive_and_merge_condition = "false"
 
        self._flux_estimates_projection.generator.load_store_compute_flag = "::peano4::grid::constructLoadStoreComputeFlag({},{},{})".format(
            self._provide_face_data_to_compute_kernels_default_guard(),
            self._load_face_data_default_guard(),
            self._store_face_data_default_guard()
        )
 
        self._flux_estimates_projection.generator.send_condition = "false"
        self._flux_estimates_projection.generator.receive_and_merge_condition = "false"
 
        self._current_time_step.generator.includes += """
#include "../repositories/SolverRepository.h"
"""
 
        self._previous_time_step.generator.includes += """
#include "../repositories/SolverRepository.h"
"""
 
        self._rhs_estimates_projection.generator.includes += """
#include "../repositories/SolverRepository.h"
"""
 
        self._flux_estimates_projection.generator.includes += """
#include "../repositories/SolverRepository.h"
"""
 
        self._cell_label = exahype2.grid.create_cell_label(self._name)
        self._face_label = exahype2.grid.create_face_label(self._name)
 
    """
    Creates the action sets, these are actions tied to the grid that are executed at given points of
    grid generation, construction or steps of the ADER-DG solver. Peano handles grid structures that
    contain data as well as operations that are executed on that data. These are the latter.
    The guard here are the conditions for the execution of the operation, as in they are only executed
    if the guard evaluates to true.
    """
    @abstractmethod
    def create_action_sets(self):
        # Updates faces and cells (whether at boundary, refined etc.)
        self._action_set_update_face_label = exahype2.grid.UpdateFaceLabel(self._name)
        self._action_set_update_cell_label = exahype2.grid.UpdateCellLabel(self._name)
        self._action_set_update_face_label.descend_invocation_order = 0
        self._action_set_update_cell_label.descend_invocation_order = 0
 
        # Sets initial conditions for the solution
        self._action_set_initial_conditions = exahype2.solvers.aderdg.actionsets.InitialCondition(
            self, self._store_cell_data_default_guard(), "true"
        )
        self._action_set_initial_conditions.descend_invocation_order = 1
 
        self._action_set_initial_conditions_for_grid_construction = exahype2.solvers.aderdg.actionsets.InitialCondition(
            self, self._store_cell_data_default_guard(), "false"
        )
        self._action_set_initial_conditions_for_grid_construction.descend_invocation_order = 1
 
        self._action_set_handle_boundary = exahype2.solvers.aderdg.actionsets.HandleBoundary(
            self, self._store_boundary_data_default_guard()
        )
        self._action_set_handle_boundary.descend_invocation_order = 1
 
        self._action_set_couple_resolution_transitions_and_handle_dynamic_mesh_refinement = exahype2.solvers.aderdg.actionsets.DynamicAMR(
            solver=self,
            delete_guard=self._delete_face_data_default_guard(),
            interpolate_guard=self._interpolate_face_data_default_guard(),
            restrict_guard=self._restrict_face_data_default_guard(),
        )
        self._action_set_couple_resolution_transitions_and_handle_dynamic_mesh_refinement.descend_invocation_order = 1
 
        # All the actual computation of the solution happens in the prediction and correction action sets.
        # Predictor extrapolates the solution from time t onto time-nodes, then projects this
        # extrapolated solution to the faces.
        # Correction uses the extrapolated face solutions from the predictor and the neighboring
        # cells to compute a solution at the next timestep.
        self._action_set_correction                   = exahype2.solvers.aderdg.actionsets.Correction(
            self, self._store_cell_data_default_guard())
        self._action_set_correction.descend_invocation_order = 2  # Needs to happen after handle_boundary
 
        self._action_set_prediction                   = exahype2.solvers.aderdg.actionsets.Prediction(
            self, self._store_cell_data_default_guard(), on_hanging_cells=False)
        self._action_set_prediction.descend_invocation_order = 3
 
        self._action_set_prediction_on_hanging_cells  = exahype2.solvers.aderdg.actionsets.Prediction(
            self, self._store_cell_data_default_guard(), on_hanging_cells=True)
        self._action_set_prediction_on_hanging_cells.descend_invocation_order = 3
 
        # Handles adaptive mesh refinement
        self._action_set_AMR = exahype2.solvers.aderdg.actionsets.AdaptivityCriterion(
            solver=self,
            guard=self._store_cell_data_default_guard()
            + " and (repositories::"
            + self.get_name_of_global_instance()
            + ".isFirstGridSweepOfTimeStep() or repositories::"
            + self.get_name_of_global_instance()
            + ".getSolverState()=="
            + self._name
            + "::SolverState::GridInitialisation)",
            build_up_new_refinement_instructions=True,
            called_by_grid_construction=False,
            implement_previous_refinement_instructions=True
        )
        self._action_set_AMR.descend_invocation_order = 4
 
        self._action_set_AMR_commit_without_further_analysis = exahype2.solvers.aderdg.actionsets.AdaptivityCriterion(
            solver=self,
            guard=self._store_cell_data_default_guard(),
            build_up_new_refinement_instructions=False,
            implement_previous_refinement_instructions=True,
            called_by_grid_construction=False,
        )
        self._action_set_AMR_commit_without_further_analysis.descend_invocation_order = 4
 
        self._action_set_AMR_throughout_grid_construction = exahype2.solvers.aderdg.actionsets.AdaptivityCriterion(
            solver=self,
            guard=self._store_cell_data_default_guard(),
            build_up_new_refinement_instructions=True,
            implement_previous_refinement_instructions=True,
            called_by_grid_construction=True,
        )
        self._action_set_AMR_throughout_grid_construction.descend_invocation_order = 4
 
        if self._action_set_postprocess_solution == None:
            self._action_set_postprocess_solution = exahype2.solvers.aderdg.actionsets.EmptyPostprocessSolution(solver=self)
        self._action_set_postprocess_solution.descend_invocation_order = 4
 
    def create_readme_descriptor(self, domain_offset, domain_size):
        return (
            """ ExaHyPE 2 ADER-DG solver of order: """
            + str(self._order)
            + """ with a domain size of: """
            + str(domain_size)
            + """ and a domain offset of: """
            + str(domain_offset)
        )
 
    def get_user_action_set_includes(self):
        return self.user_action_set_includes
 
    def get_user_solver_includes(self):
        return self.user_solver_includes
 
    def add_user_action_set_includes(self, value):
        """
        Add further includes to this property, if your action sets require some additional
        routines from other header files.
        """
        self.user_action_set_includes += value
 
    def add_user_solver_includes(self, value):
        """
        Add further includes to this property, if your solver requires some additional
        routines from other header files.
        """
        self.user_solver_includes += value
 
    def _store_cell_data_default_guard(self):
        return (
            "not marker.willBeRefined() "
            + "and repositories::"
            + self.get_name_of_global_instance()
            + ".getSolverState()!="
            + self._name
            + "::SolverState::GridConstruction"
        )
 
    def _provide_cell_data_to_compute_kernels_default_guard(self):
        return (
                    "not marker.willBeRefined() "
                    + "and repositories::"
                    + self.get_name_of_global_instance()
                    + ".getSolverState()!="
                    + self._name
                    + "::SolverState::GridConstruction"
                )
 
    def _provide_face_data_to_compute_kernels_default_guard(self):
        return (
                    "not marker.willBeRefined() "
                    + "and repositories::"
                    + self.get_name_of_global_instance()
                    + ".getSolverState()!="
                    + self._name
                    + "::SolverState::GridConstruction"
                )
 
    def _load_cell_data_default_guard(self):
        return (
            "not marker.hasBeenRefined() "
            + "and repositories::"
            + self.get_name_of_global_instance()
            + ".getSolverState()!="
            + self._name
            + "::SolverState::GridConstruction "
            + "and repositories::"
            + self.get_name_of_global_instance()
            + ".getSolverState()!="
            + self._name
            + "::SolverState::GridInitialisation"
        )
 
    def _store_face_data_default_guard(self):
        return (
            "not marker.willBeRefined() "
            + "and repositories::"
            + self.get_name_of_global_instance()
            + ".getSolverState()!="
            + self._name
            + "::SolverState::GridConstruction"
        )
 
    def _load_face_data_default_guard(self):
        return (
            "not marker.hasBeenRefined() "
            + "and repositories::"
            + self.get_name_of_global_instance()
            + ".getSolverState()!="
            + self._name
            + "::SolverState::GridConstruction "
            + "and repositories::"
            + self.get_name_of_global_instance()
            + ".getSolverState()!="
            + self._name
            + "::SolverState::GridInitialisation"
        )
 
    def _store_boundary_data_default_guard(self):
        return (
            self._store_face_data_default_guard()
            + " and not repositories::"
            + self.get_name_of_global_instance()
            + ".PeriodicBC[marker.getSelectedFaceNumber()%Dimensions]"
            + " and not marker.hasBeenRefined() and fineGridFace"
            + self._name
            + "FaceLabel.getBoundary()"
        )
 
    def _delete_face_data_default_guard(self):
        return "true"
 
    def _interpolate_face_data_default_guard(self):
        return (
            "repositories::"
            + self.get_name_of_global_instance()
            + ".getSolverState()!="
            + self._name
            + "::SolverState::GridInitialisation"
        )
 
    def _restrict_face_data_default_guard(self):
        return "true"
 
    def _unknown_identifier(self):
        return self._name + "Q"
 
    def get_name_of_global_instance(self):
        return "instanceOf" + self._name
 
    def _get_default_includes(self):
        return """
#include "tarch/la/Vector.h"
 
#include "peano4/utils/Globals.h"
#include "peano4/utils/Loop.h"
 
#include "repositories/SolverRepository.h"
""" + "#include <stdfloat>" if self._use_half_precision else ""
 
    """
    The following functions will be called by the Peano 4 project,
    they tell Peano which of our generated data it should attach
    to each grid element and use.
    """
    def add_to_Peano4_datamodel(self, datamodel, verbose):
        if verbose:
          print("Patch data")
          print("----------")
          print(str(self._current_time_step))
          print("Patch overlap data")
          print("----------")
          print(str(self._rhs_estimates_projection))
          print(str(self._flux_estimates_projection))
        datamodel.add_cell(self._cell_label)
        datamodel.add_cell(self._current_time_step)
        datamodel.add_cell(self._previous_time_step)
        datamodel.add_face(self._face_label)
        datamodel.add_face(self._rhs_estimates_projection)
        datamodel.add_face(self._flux_estimates_projection)
 
    def add_use_data_statements_to_Peano4_solver_step(self, step):
        step.use_cell(self._cell_label)
        step.use_cell(self._current_time_step)
        step.use_cell(self._previous_time_step)
        step.use_face(self._face_label)
        step.use_face(self._rhs_estimates_projection)
        step.use_face(self._flux_estimates_projection)
 
    def add_actions_to_init_grid(self, step):
        step.add_action_set(self._action_set_initial_conditions)
        step.add_action_set(self._action_set_update_face_label)
        step.add_action_set(self._action_set_update_cell_label)
        step.add_action_set(self._action_set_AMR)
        step.add_action_set(self._action_set_postprocess_solution)
 
    def add_actions_to_create_grid(self, step, evaluate_refinement_criterion):
        step.add_action_set(self._action_set_initial_conditions_for_grid_construction)
        step.add_action_set(self._action_set_update_face_label)
        step.add_action_set(self._action_set_update_cell_label)
        if evaluate_refinement_criterion:
            step.add_action_set(self._action_set_AMR_throughout_grid_construction)
        else:
            step.add_action_set(self._action_set_AMR_commit_without_further_analysis)
 
    def add_actions_to_plot_solution(self, step, output_path):
        d = {}
        self._init_dictionary_with_default_parameters(d)
        self.add_entries_to_text_replacement_dictionary(d)
 
        step.add_action_set(self._action_set_update_face_label)
        step.add_action_set(self._action_set_update_cell_label)
 
        mapping = []
        for z in self._basis.quadrature_points():
            for y in self._basis.quadrature_points():
                for x in self._basis.quadrature_points():
                    mapping.append((x, y, z))
 
        # Calls a generator which generates a standard plotter that uses the tarch library.
        step.add_action_set(
            peano4.toolbox.blockstructured.PlotPatchesInPeanoBlockFormat(
                filename=output_path + "solution-" + self._name,
                patch=self._previous_time_step,
                dataset_name=self._unknown_identifier()+"_old",
                description=self.plot_description,
                mapping=mapping,
                select_dofs=self.select_dofs_to_print,
                guard="repositories::plotFilter.plotPatch(marker) and "
                + self._load_cell_data_default_guard(),
                additional_includes="""
#include "exahype2/PlotFilter.h"
#include "../repositories/SolverRepository.h"
""",
                precision="PlotterPrecision",
                time_stamp_evaluation="0.5*(repositories::getMinTimeStamp()+repositories::getMaxTimeStamp())",
                plot_cell_data=False,
                dataType=self._solution_persistent_storage_precision,
            )
        )
 
        step.add_action_set(self._action_set_couple_resolution_transitions_and_handle_dynamic_mesh_refinement)
 
        if self._plot_grid_properties:
            step.add_action_set(
                peano4.toolbox.PlotGridInPeanoBlockFormat(
                    filename=output_path + "grid-" + self._name,
                    cell_unknown=None,
                    guard_guard="repositories::plotFilter.plotPatch(marker) and "
                    + self._load_cell_data_default_guard(),
                    additional_includes="""
#include "exahype2/PlotFilter.h"
#include "../repositories/SolverRepository.h"
""",
                    plot_cell_data=False,
                )
            )
        pass
 
    # These actions are executed during each time step
    def add_actions_to_perform_time_step(self, step):
        step.add_action_set(self._action_set_update_face_label)
        step.add_action_set(self._action_set_update_cell_label)
        step.add_action_set(self._action_set_couple_resolution_transitions_and_handle_dynamic_mesh_refinement)
 
        step.add_action_set(self._action_set_handle_boundary)
        step.add_action_set(self._action_set_correction)
        step.add_action_set(self._action_set_prediction)
        step.add_action_set(self._action_set_prediction_on_hanging_cells)
 
        step.add_action_set(self._action_set_AMR)
        step.add_action_set(self._action_set_postprocess_solution)
 
    def set_plot_description(self, description):
        self.plot_description = description
 
    def _generate_kernels(self, namespace, output, subdirectory, dimensions):
        if (not self._is_linear) and (
            self._material_param_implementation != PDETerms.None_Implementation
            or self._point_sources_implementation != PDETerms.None_Implementation
        ):
            raise Exception(
                "Material parameters and point sources are currently only available for linear generated kernels."
            )
 
        d = {}
        self._init_dictionary_with_default_parameters(d)
        self.add_entries_to_text_replacement_dictionary(d)
        d["DIMENSIONS"] = dimensions
        d["NUMBER_OF_DOF_3D"] = 1 if dimensions == 2 else d["NUMBER_OF_DOF"]
        d["FULL_QUALIFIED_SOLVER_NAME"] = "::".join(namespace) + "::" + self._name
 
        kernels_namespace = namespace + ["kernels", self._name]
        kernels_output_path = subdirectory + "kernels/" + self._name
        kernels_templates_prefix = os.path.dirname(os.path.realpath(__file__)) + "/kernels/"
 
        if not os.path.exists(kernels_output_path):
           os.makedirs(kernels_output_path)
 
        generated_kernels = (
            peano4.output.Jinja2TemplatedHeaderImplementationFilePair(
                headfile_template=kernels_templates_prefix + "CellErrorIntegral.template.h",
                cppfile_template=kernels_templates_prefix + "CellErrorIntegral.template.cpp",
                classname=kernels_output_path + "/CellErrorIntegral",
                namespace=kernels_namespace,
                subdirectory=subdirectory + ".",
                dictionary=d,
                default_overwrite=True,
                apply_iteratively=True,
            )
        )
        output.add(generated_kernels)
        output.makefile.add_h_file(kernels_output_path + "/CellErrorIntegral.h", generated=True)
        output.makefile.add_cpp_file(kernels_output_path + "/CellErrorIntegral.cpp", generated=True)
 
        generated_kernels = (
            peano4.output.Jinja2TemplatedHeaderImplementationFilePair(
                headfile_template=kernels_templates_prefix + "MaxScaledEigenvalue.template.h",
                cppfile_template=kernels_templates_prefix + "MaxScaledEigenvalue.template.cpp",
                classname=kernels_output_path + "/MaxScaledEigenvalue",
                namespace=kernels_namespace,
                subdirectory=subdirectory + ".",
                dictionary=d,
                default_overwrite=True,
                apply_iteratively=True,
            )
        )
        output.add(generated_kernels)
        output.makefile.add_h_file(kernels_output_path + "/MaxScaledEigenvalue.h", generated=True)
        output.makefile.add_cpp_file(kernels_output_path + "/MaxScaledEigenvalue.cpp", generated=True)
 
        generated_kernels = (
            peano4.output.Jinja2TemplatedHeaderImplementationFilePair(
                headfile_template=kernels_templates_prefix + "FaceIntegral.template.h",
                cppfile_template=kernels_templates_prefix + "FaceIntegral.template.cpp",
                classname=kernels_output_path + "/FaceIntegral",
                namespace=kernels_namespace,
                subdirectory=subdirectory + ".",
                dictionary=d,
                default_overwrite=True,
                apply_iteratively=True,
            )
        )
        output.add(generated_kernels)
        output.makefile.add_h_file(kernels_output_path + "/FaceIntegral.h", generated=True)
        output.makefile.add_cpp_file(kernels_output_path + "/FaceIntegral.cpp", generated=True)
 
        generated_kernels = (
            peano4.output.Jinja2TemplatedHeaderImplementationFilePair(
                headfile_template=kernels_templates_prefix + "RiemannSolver.template.h",
                cppfile_template=kernels_templates_prefix + "RiemannSolverLinear.template.cpp" if self._is_linear
                else kernels_templates_prefix + "RiemannSolverNonlinear.template.cpp",
                classname=kernels_output_path + "/RiemannSolver",
                namespace=kernels_namespace,
                subdirectory=subdirectory + ".",
                dictionary=d,
                default_overwrite=True,
                apply_iteratively=True,
            )
        )
        output.add(generated_kernels)
        output.makefile.add_h_file(kernels_output_path + "/RiemannSolver.h", generated=True)
        output.makefile.add_cpp_file(kernels_output_path + "/RiemannSolver.cpp", generated=True)
 
        exahype2.solvers.aderdg.kernels.AMRRoutines(d).generate_kernels(
            kernels_namespace, output, subdirectory, kernels_templates_prefix, kernels_output_path
        )
        exahype2.solvers.aderdg.kernels.CellData(d).generate_kernels(
            kernels_namespace, output, subdirectory, kernels_templates_prefix, kernels_output_path
        )
        exahype2.solvers.aderdg.kernels.DGMatrices(d).generate_kernels(
            kernels_namespace, output, subdirectory, kernels_templates_prefix, kernels_output_path
        )
        exahype2.solvers.aderdg.kernels.Quadrature(d).generate_kernels(
            kernels_namespace, output, subdirectory, kernels_templates_prefix, kernels_output_path
        )
        exahype2.solvers.aderdg.kernels.FusedSpaceTimePredictorVolumeIntegral(d).generate_kernels(
            kernels_namespace, output, subdirectory, kernels_templates_prefix, kernels_output_path
        )
 
        if self._point_sources_implementation != PDETerms.None_Implementation:
          generated_kernels = (
            peano4.output.Jinja2TemplatedHeaderImplementationFilePair(
                headfile_template=kernels_templates_prefix + "PointSources.template.h",
                cppfile_template=kernels_templates_prefix + "PointSources.template.cpp",
                classname=kernels_output_path + "/PointSources",
                namespace=kernels_namespace,
                subdirectory=subdirectory + ".",
                dictionary=d,
                default_overwrite=True,
                apply_iteratively=True,
            )
          )
          output.add(generated_kernels)
          output.makefile.add_h_file(kernels_output_path + "/PointSources.h", generated=True)
          output.makefile.add_cpp_file(kernels_output_path + "/PointSources.cpp", generated=True)
 
        if self._use_libxsmm:
          generated_kernels = (
            peano4.output.Jinja2TemplatedHeaderImplementationFilePair(
                headfile_template=kernels_templates_prefix + "Gemms.template.h",
                cppfile_template=kernels_templates_prefix + "Gemms.template.cpp",
                classname=kernels_output_path + "/Gemms",
                namespace=kernels_namespace,
                subdirectory=subdirectory + ".",
                dictionary=d,
                default_overwrite=True,
                apply_iteratively=True,
            )
          )
          output.add(generated_kernels)
          output.makefile.add_h_file(kernels_output_path + "/Gemms.h", generated=True)
          output.makefile.add_cpp_file(kernels_output_path + "/Gemms.cpp", generated=True)
 
        """"
        This is a very ugly fix that should be replaced asap.
        Essentially, the libxsmm kernel generator generates faulty code if
        'noarch' is specified. This code works fine except that it contains a
        #ifdef which is missing its corresponding #endif.
        This code just inserts that missing #endif.
        """
        if self._use_libxsmm and self._architecture == "noarch":
            file = open("kernels/" + self._name + "/FusedSpaceTimePredictorVolumeIntegral_libxsmm.c", "r")
            content = file.read()
            file.close()
            file = open("kernels/" + self._name + "/FusedSpaceTimePredictorVolumeIntegral_libxsmm.c", "w")
            content = content.replace("__FILE__)\n", "__FILE__)\n#endif\n")
            file.write(content)
            file.close()
 
            file = open("kernels/" + self._name + "/AMRRoutines_libxsmm.c", "r")
            content = file.read()
            file.close()
            file = open("kernels/" + self._name + "/AMRRoutines_libxsmm.c", "w")
            content = content.replace("__FILE__)\n", "__FILE__)\n#endif\n")
            file.write(content)
            file.close()
 
    """
    This tells Peano to add the solver files to the project. It generates any files that are not
    specifically cell or face data or grid actions such as the actionsets.
    Currently it generates solver implementation and header files as well as abstract solver
    implementation and header files.
    In addition to this if the user has specified that they would like to use shortcuts to refer
    to their variable (e.g., specified their unknowns or aux. variables through a dict) it generates
    a VariableShortcuts file which is then connected to the solver so that users can specify their
    variables through name.
    """
    def add_implementation_files_to_project(self, namespace, output, dimensions, subdirectory=""):
        """
        The ExaHyPE2 project will call this operation when it sets
        up the overall environment.
 
        This routine is typically not invoked by a user.
 
        output: peano4.output.Output
        """
        templatefile_prefix = os.path.dirname(os.path.realpath(__file__)) + "/"
 
        if self._solver_template_file_class_name is None:
            templatefile_prefix += self.__class__.__name__
        else:
            templatefile_prefix += self._solver_template_file_class_name
 
        if subdirectory:
            subdirectory += "/"
 
        abstractHeaderDictionary = {}
        implementationDictionary = {}
        self._init_dictionary_with_default_parameters(abstractHeaderDictionary)
        self._init_dictionary_with_default_parameters(implementationDictionary)
        self.add_entries_to_text_replacement_dictionary(abstractHeaderDictionary)
        self.add_entries_to_text_replacement_dictionary(implementationDictionary)
 
        generated_abstract_header_file = (
            peano4.output.Jinja2TemplatedHeaderImplementationFilePair(
                templatefile_prefix + "Abstract.template.h",
                templatefile_prefix + "Abstract.template.cpp",
                "Abstract" + self._name,
                namespace,
                subdirectory + ".",
                abstractHeaderDictionary,
                True,
                True,
            )
        )
        generated_solver_files = (
            peano4.output.Jinja2TemplatedHeaderImplementationFilePair(
                templatefile_prefix + ".template.h",
                templatefile_prefix + ".template.cpp",
                self._name,
                namespace,
                subdirectory + ".",
                implementationDictionary,
                False,
                True,
            )
        )
 
        output.add(generated_abstract_header_file)
        output.add(generated_solver_files)
        output.makefile.add_h_file(subdirectory + "Abstract" + self._name + ".h", generated=True)
        output.makefile.add_h_file(subdirectory + self._name + ".h", generated=True)
        output.makefile.add_cpp_file(subdirectory + "Abstract" + self._name + ".cpp", generated=True)
        output.makefile.add_cpp_file(subdirectory + self._name + ".cpp", generated=True)
 
        if self._use_var_shortcut:
            generated_shortcut_file = peano4.output.Jinja2TemplatedHeaderFile(
                os.path.dirname(os.path.realpath(__file__))
                + "/"
                + "../VariableShortcuts.template.h",
                "VariableShortcuts",
                namespace,
                subdirectory + ".",
                implementationDictionary,
                True,
                True,
            )
            output.add(generated_shortcut_file)
            output.makefile.add_h_file(subdirectory + "VariableShortcuts.h", generated=True)
 
        self._generate_kernels(namespace, output, subdirectory, dimensions)
 
    def set_solver_constants(self, datastring):
        self._solver_constants = datastring
 
    def add_solver_constants(self, datastring):
        self._solver_constants += datastring
 
    @property
    def order(self):
        return self._order
 
    @property
    def unknowns(self):
        return self._unknowns
 
    @property
    def auxiliary_variables(self):
        return self._auxiliary_variables
 
    @property
    def preprocess_reconstructed_patch(self):
        return self._preprocess_reconstructed_patch
 
    @preprocess_reconstructed_patch.setter
    def preprocess_reconstructed_patch(self, kernel):
        self._preprocess_reconstructed_patch = kernel
        self.create_data_structures()
        self.create_action_sets()
 
    @property
    def postprocess_updated_patch(self):
        return self._postprocess_updated_patch
 
    @postprocess_updated_patch.setter
    def postprocess_updated_patch(self, kernel):
        self._postprocess_updated_patch = kernel
        self.create_data_structures()
        self.create_action_sets()
 
    @abstractmethod
    def add_entries_to_text_replacement_dictionary(self, d):
        pass
 
    """
    Generates a dictionary of "words" that will later be used in various templates to fill these out by adding
    information from the solver or which was specified by the user.
    """
    def _init_dictionary_with_default_parameters(self, d):
        d["SOLVER_INSTANCE"]    = self.get_name_of_global_instance()
        d["SOLVER_NAME"]        = self._name
        d["UNKNOWN_IDENTIFIER"] = self._unknown_identifier()
 
        d["ALIGNMENT"]          = Architectures.get_alignment(self._architecture)
        d["SIMD_WIDTH"]         = Architectures.get_simd_width(self._architecture)
        d["ORDER"]              = self._order
        d["NUMBER_OF_DOF"]      = self._order + 1
        d["NUMBER_OF_DOF_PADDED"] = Architectures.get_size_with_padding(self._architecture, self._order + 1)
        d["NUMBER_OF_UNKNOWNS"]             = self._unknowns
        d["NUMBER_OF_UNKNOWNS_PADDED"]      = Architectures.get_size_with_padding(self._architecture, self._unknowns)
        d["NUMBER_OF_AUXILIARY_VARIABLES"]  = self._auxiliary_variables
        d["NUMBER_OF_AUXILIARY_VARIABLES_PADDED"] = Architectures.get_size_with_padding(self._architecture, self._auxiliary_variables)
        d["NUMBER_OF_DATA"]                     = self._unknowns + self._auxiliary_variables
        d["NUMBER_OF_DATA_PADDED"]              = Architectures.get_size_with_padding(self._architecture, self._unknowns + self._auxiliary_variables)
 
        d["CFL_SAFETY_FACTOR"]  = self._cflSafetyFactor
        d["CFL_ADER"]           = self._cflAder[self._order]  # /(2*(self._order+1)+1)
 
        d["USE_HALF_PRECISION"] = self._use_half_precision
        d["IS_LINEAR"] = self._is_linear
        d["LINEARITY"] = "linear" if self._is_linear else "nonlinear"
 
        d["USE_GAUSS_LOBATTO"] = ( "true" if self._polynomials is Polynomials.Gauss_Lobatto else "false" )
        d["POLYNOMIAL_TYPE"]   = ( "lobatto" if self._polynomials is Polynomials.Gauss_Lobatto else "legendre" )
 
        d["SOLUTION_STORAGE_PRECISION"] = self._solution_persistent_storage_precision
 
        d["PREDICTOR_COMPUTATION_PRECISIONS"] = self._predictor_computation_precisions
        d["CORRECTOR_COMPUTATION_PRECISION"]  = self._corrector_computation_precision
 
        d["COMPUTATION_PRECISIONS"] = self._predictor_computation_precisions[:]
        if self._corrector_computation_precision not in d["COMPUTATION_PRECISIONS"]:
            d["COMPUTATION_PRECISIONS"].append(self._corrector_computation_precision)
        if self._precompute_picard_precision!=False and self._precompute_picard_precision not in d["COMPUTATION_PRECISIONS"]:
            d["COMPUTATION_PRECISIONS"].append(self._precompute_picard_precision)
        d["NUMBER_OF_PREDICTOR_PRECISIONS"] = len(d["PREDICTOR_COMPUTATION_PRECISIONS"])
        d["NUMBER_OF_PRECISIONS"] = len(d["COMPUTATION_PRECISIONS"])
 
        if self._min_cell_h > self._max_cell_h:
            raise Exception("min/max h are inconsistent")
        d["MAX_CELL_H"] = self._max_cell_h
        d["MIN_CELL_H"] = self._min_cell_h
 
        d["SOLVER_CONSTANTS"] = self._solver_constants
        d["SOLVER_INCLUDES"]  = self.get_user_solver_includes()
 
        d["BOUNDARY_CONDITIONS_IMPLEMENTATION"] = self._boundary_conditions_implementation
        d["REFINEMENT_CRITERION_IMPLEMENTATION"] = self._refinement_criterion_implementation
        d["INITIAL_CONDITIONS_IMPLEMENTATION"] = self._initial_conditions_implementation
        d["ADAPTIVE_PRECISION_IMPLEMENTATION"] = self._precision_criterion_implementation
        d["RIEMANN_SOLVER_IMPLEMENTATION"] = self._riemann_solver_implementation
 
        d["ABSTRACT_SOLVER_USER_DECLARATIONS"] = jinja2.Template(
            self._abstract_solver_user_declarations, undefined=jinja2.DebugUndefined
        ).render(**d)
        d["ABSTRACT_SOLVER_USER_DEFINITIONS"] = jinja2.Template(
            self._abstract_solver_user_definitions, undefined=jinja2.DebugUndefined
        ).render(**d)
        d["SOLVER_USER_DECLARATIONS"] = jinja2.Template(
            self._solver_user_declarations, undefined=jinja2.DebugUndefined
        ).render(**d)
        d["SOLVER_USER_DEFINITIONS"] = jinja2.Template(
            self._solver_user_definitions, undefined=jinja2.DebugUndefined
        ).render(**d)
 
        d["START_TIME_STEP_IMPLEMENTATION"] = jinja2.Template(
            self._start_time_step_implementation, undefined=jinja2.DebugUndefined
        ).render(**d)
        d["FINISH_TIME_STEP_IMPLEMENTATION"] = jinja2.Template(
            self._finish_time_step_implementation, undefined=jinja2.DebugUndefined
        ).render(**d)
 
        d["CONSTRUCTOR_IMPLEMENTATION"] = jinja2.Template(
            self._constructor_implementation, undefined=jinja2.DebugUndefined
        ).render(**d)
        d["DESTRUCTOR_IMPLEMENTATION"] = jinja2.Template(
            self._destructor_implementation, undefined=jinja2.DebugUndefined
        ).render(**d)
        d["COMPUTE_TIME_STEP_SIZE"] = jinja2.Template(
            self._compute_time_step_size, undefined=jinja2.DebugUndefined
        ).render(**d)
        d["COMPUTE_NEW_TIME_STEP_SIZE"] = jinja2.Template(
            self._compute_new_time_step_size, undefined=jinja2.DebugUndefined
        ).render(**d)
 
        d["PREPROCESS_RECONSTRUCTED_PATCH"] = jinja2.Template(
            self._preprocess_reconstructed_patch, undefined=jinja2.DebugUndefined
        ).render(**d)
        d["POSTPROCESS_UPDATED_PATCH"] = jinja2.Template(
            self._postprocess_updated_patch, undefined=jinja2.DebugUndefined
        ).render(**d)
 
        d["USE_FLUX"]                 = self._flux_implementation != PDETerms.None_Implementation
        d["USE_NCP"]                  = self._ncp_implementation != PDETerms.None_Implementation
        d["USE_SOURCE"]               = self._source_term_implementation != PDETerms.None_Implementation
        d["USE_MATERIAL_PARAMETERS"]  = self._material_param_implementation != PDETerms.None_Implementation
        d["USE_POINT_SOURCE"]         = self._point_sources_implementation != PDETerms.None_Implementation
        d["NUMBER_OF_POINT_SOURCES"]  = self._number_of_point_sources
 
        d["USE_VARIABLE_SHORTCUT"]  = self._use_var_shortcut
        d["VARIABLE_NAMES"]         = self._variable_names
        d["VARIABLE_POSITIONS"]     = self._variable_pos
 
        d["QUADRATURE_POINTS"]      = ", ".join(self._basis.quadrature_points())
 
        d["USE_LIBXSMM"] = self._use_libxsmm
        d["ARCHITECTURE"] = self._architecture
 
        # self._basis.init_dictionary_with_default_parameters(d)
 
    def set_implementation(
        self,
        flux,
        ncp,
        eigenvalues,
        boundary_conditions,
        refinement_criterion,
        initial_conditions,
        source_term,
        material_parameters,
        point_source,
        additional_action_set_includes,
        additional_user_includes,
    ):
        """
        If you pass in User_Defined, then the generator will create C++ stubs
        that you have to befill manually. If you pass in None_Implementation, it
        will create nop, i.e., no implementation or defaults. Any other string
        is copied 1:1 into the implementation. If you pass in None, then the
        set value so far won't be overwritten.
        """
        if boundary_conditions is not None:
            self._boundary_conditions_implementation = boundary_conditions
        if refinement_criterion is not None:
            self._refinement_criterion_implementation = refinement_criterion
        if initial_conditions is not None:
            self._initial_conditions_implementation = initial_conditions
 
        if flux is not None:
            self._flux_implementation = flux
        if ncp is not None:
            self._ncp_implementation = ncp
        if eigenvalues is not None:
            self._eigenvalues_implementation = eigenvalues
        if source_term is not None:
            self._source_term_implementation = source_term
        if material_parameters is not None:
            self._material_param_implementation = material_parameters
 
        if type(point_source) != int:
            raise Exception(
                "point_source needs to be an integer, this determines the number of point sources that will be used."
            )
 
        if point_source > 0:
            self._point_sources_implementation = PDETerms.User_Defined_Implementation
            self._number_of_point_sources = point_source
 
        self.user_action_set_includes += additional_action_set_includes
        self.user_solver_includes += additional_user_includes
 
        self.create_data_structures()
        self.create_action_sets()
 
    def add_kernel_optimisations(
        self,
        is_linear                             = False,
        polynomials                           = Polynomials.Gauss_Legendre,
        use_libxsmm                           = False,
        precision                             = False,
        predictor_computation_precisions      = False,
        corrector_computation_precision       = False,
        solution_persistent_storage_precision = False,
        precompute_picard_precision           = False,
        precision_criterion_implementation    = PDETerms.User_Defined_Implementation,
        riemann_solver_implementation         = PDETerms.None_Implementation,
        architecture                          = "noarch",
        initialise_patches                    = False
    ):
        if polynomials is Polynomials.Gauss_Legendre:
            self._basis = GaussLegendreBasis(self._order + 1)
        elif polynomials is Polynomials.Gauss_Lobatto:
            self._basis = GaussLobattoBasis(self._order + 1)
        else:
            raise Exception(
                "No proper basis chosen: {}, valid options are Gauss_Legendre and Gauss_Lobatto nodes".format(
                    polynomials
                )
            )
 
        self._is_linear = is_linear
        self._polynomials = polynomials
 
        self._use_libxsmm = use_libxsmm
 
        if precision != False:
            if precision not in PrecisionType.keys():
                raise AssertionError(
                    "The chosen precision was not a valid choice. valid choices are: "
                    + ", ".join(PrecisionType.keys())
                )
            else:
                self._predictor_computation_precisions      = [PrecisionType[precision]]
                self._precompute_picard_precision           = PrecisionType[precision]
                self._solution_persistent_storage_precision = PrecisionType[precision]
                self._corrector_computation_precision       = PrecisionType[precision]
 
        if solution_persistent_storage_precision != False:
            if solution_persistent_storage_precision not in PrecisionType.keys():
                raise AssertionError(
                    "The chosen precision was not a valid choice. valid choices are: "
                    + ", ".join(PrecisionType.keys())
                )
            self._solution_persistent_storage_precision = PrecisionType[solution_persistent_storage_precision]
 
        if corrector_computation_precision != False:
            if corrector_computation_precision not in PrecisionType.keys():
                raise AssertionError(
                    "The chosen precision was not a valid choice. valid choices are: "
                    + ", ".join(PrecisionType.keys())
                )
            self._corrector_computation_precision = PrecisionType[corrector_computation_precision]
 
        if precompute_picard_precision != False:
            if precompute_picard_precision not in PrecisionType.keys():
                raise AssertionError(
                    "The chosen precision was not a valid choice. valid choices are: "
                    + ", ".join(PrecisionType.keys())
                )
            self._precompute_picard_precision = PrecisionType[precompute_picard_precision]
 
        if predictor_computation_precisions != False:
            self._predictor_computation_precisions = []
            for prec in predictor_computation_precisions:
                if prec not in PrecisionType.keys():
                    raise AssertionError(
                        "The chosen precision was not a valid choice. valid choices are: "
                        + ", ".join(PrecisionType.keys())
                    )
                self._predictor_computation_precisions.append(PrecisionType[prec])
 
        if len(self._predictor_computation_precisions) > 1:
            self._precision_criterion_implementation = precision_criterion_implementation
 
        if riemann_solver_implementation is not PDETerms.None_Implementation:
            self._riemann_solver_implementation = riemann_solver_implementation
 
        self._architecture = architecture
 
        self._initialise_patches = initialise_patches
 
        self.create_data_structures()
        self.create_action_sets()
 
    def switch_storage_scheme(
        self,
        cell_data_storage: Storage,
        face_data_storage: Storage,
    ):
        """
        By default, we hold all data on the heap using smart pointers. You can explicitly switch
        to storage on the call stack or heap using raw pointers.
 
        @see create_data_structures()
        """
        assert isinstance(cell_data_storage, Storage)
        assert isinstance(face_data_storage, Storage)
 
        self._cell_data_storage = cell_data_storage
        self._face_data_storage = face_data_storage
 
        self.create_data_structures()
        self.create_action_sets()
 
def get_face_overlap_merge_implementation(patch_overlap):
    d = {
        "OVERLAP": patch_overlap.dim[0]
        / 2,  # Assumed to always be two, but present here anyway.
        "DOFS_PER_AXIS": patch_overlap.dim[1],
        "UNKNOWNS": patch_overlap.no_of_unknowns,
    }
    Template = """
  #if Dimensions == 2
  constexpr int SpaceFaceSize = {DOFS_PER_AXIS} * {UNKNOWNS}; // Order + 1
  #else
  constexpr int SpaceFaceSize = {DOFS_PER_AXIS} * {DOFS_PER_AXIS} * {UNKNOWNS}; // (Order + 1)^2
  #endif
 
  const int faceDimension = marker.getSelectedFaceNumber() % Dimensions;
  // If the outer normal is positive, the normal points to the right so the face is on the right.
  const int faceLR = (marker.outerNormal()[faceDimension] > 0.0 ? 0 : 1);
 
  // std::copy_n(from, how_many, to);
  std::copy_n(
    &neighbour.value[(1 - faceLR) * SpaceFaceSize],
    SpaceFaceSize,
    &value[(1 - faceLR) * SpaceFaceSize]
  );
"""
    return Template.format(**d)