How to add a new PDE type to Inciter

This page describes how to add a new PDE type to Inciter.

Inciter supports multiple partial differential equation (PDE) types. At the time this was written, there were two PDE types: Transport and CompFlow. The former is specialized to solving transport equations of single or multiple independent scalars, governed by a set of independent advection-diffusion equations, while the latter is specialized to solving the Euler and Navier-Stokes equations for single-material compressible ideal gases. The simplest way to learn about the currently implemented PDE types and how their source is organized into directories and files is to do

cd <quinoa>/src/PDE
tree

which at this time gives

.
├── CGPDE.h
├── CMakeLists.txt
├── CompFlow
│   ├── CGCompFlow.h
│   ├── DGCompFlow.h
│   ├── Physics
│   │   ├── CGEuler.h
│   │   ├── CG.h
│   │   ├── CGNavierStokes.h
│   │   ├── DGEuler.h
│   │   ├── DG.h
│   │   └── DGNavierStokes.h
│   ├── Problem
│   │   ├── NLEnergyGrowth.h
│   │   ├── RayleighTaylor.h
│   │   ├── SedovBlastwave.h
│   │   ├── SodShocktube.h
│   │   ├── TaylorGreen.h
│   │   ├── UserDefined.h
│   │   └── VorticalFlow.h
│   └── Problem.h
├── ConfigureCompFlow.C
├── ConfigureCompFlow.h
├── ConfigureTransport.C
├── ConfigureTransport.h
├── DGPDE.h
├── FunctionPrototypes.h
├── Integrate
│   ├── Boundary.C
│   ├── Boundary.h
│   ├── Initialize.C
│   ├── Initialize.h
│   ├── Mass.C
│   ├── Mass.h
│   ├── Quadrature.C
│   ├── Quadrature.h
│   ├── Riemann
│   │   ├── HLLC.h
│   │   ├── LaxFriedrichs.h
│   │   ├── RiemannFactory.C
│   │   ├── RiemannFactory.h
│   │   ├── RiemannSolver.h
│   │   └── Upwind.h
│   ├── Source.C
│   ├── Source.h
│   ├── Surface.C
│   ├── Surface.h
│   ├── Volume.C
│   └── Volume.h
├── Limiter.C
├── Limiter.h
├── PDEFactory.h
├── PDEStack.C
├── PDEStack.h
└── Transport
    ├── CGTransport.h
    ├── DGTransport.h
    ├── Physics
    │   ├── CGAdvDiff.h
    │   ├── CGAdvection.h
    │   ├── CG.h
    │   ├── DGAdvection.h
    │   └── DG.h
    ├── Problem
    │   ├── CylAdvect.h
    │   ├── GaussHump.h
    │   ├── ShearDiff.h
    │   └── SlotCyl.h
    └── Problem.h

8 directories, 61 files

Once can see separate directories for CompFlow and Transport, each with subdirectories Physics and Problem. For example, two PDE classes are implemented under CompFlow: CGCompFlow and DGCompFlow, using continuous and discoontinuous Galerkin finite element (CG, DG) discretizations, respectively,

This page describes how to add a new PDE, such as Transport or CompFlow.

Rationale and plan

This page will discuss how to add a new PDE type by adding one for multi-material compressible flow using a DG finite element discretization. We will explain how and where to create the new files for the PDE class, its Physics and Problem policy classes, and how to define these so that the existing discontinuous Galerkin integration driver class, DG, can drive it , potentially coupled to other existing and future DG systems. The new PDE type will be instantiated from a factory based on user input, similar to the existing PDE types. What will not be described are the details of the implementation of the discretization of the newly added equation system and its algoritm. We describe how to create a skeleton PDE type in the existing framework and how to build on the existing software infrastructure that supports such a new algorithm, e.g., user input, mesh I/O, partitioning using over-decomposition, automatic load balancing via migration, etc.

1. Specification of a PDE in the input file

A PDE type, or more precisely, a PDE system is selected by the user via specifying a PDE-block in the Inciter input file. An example, configuring the existing CompFlow PDE is given by a compflow...end block,

inciter
  ...
  scheme dg
  ...
  compflow
    depvar u    # assign dependent variable 'u' to this PDE
    physics euler       # select 'euler' as compressible flow physics policy
    problem vortical_flow       # select 'vortical flow` as problem policy
    ...
  end
  ...
end

This is a snippet, configuring inciter to compute a single PDE system governing compressible flow. Every PDE type is configured via a physics and a problem policy. Here the physics is euler, configuring the solution of the Euler equations for inviscid flow, and the problem is vortical_flow, which is a test case derived using the method of manufactured solutions, specifying a prescribed velocity field with vortices. Since the analytical solution is known for this specific case, it can be used to verify the correctness of the implementation of the algorithm and the order of accuracy of the numerical method.

The important parts of this specification are the following.

  1. A new block configures a particular PDE system, such as compflow...end. There can be multiple such PDE-blocks specified, which then will be integrated together in an unsplit strongly-coupled fashion, e.g, compflow + transport with, say, 12 scalar components, which then together will comprise a system of 17 coupled scalar equations.
  2. The PDE must be specialized with a selected physics type. In the example above this is euler, but this is only the simplest choice. Every PDE has a default physics which is selected if no physics is specified.
  3. The PDE must be specialized with a selected problem type. In the example above this is vortical_flow, but this is only one choice. Every PDE has a default problem which is selected if no problem is specified.
  4. The depvar keyword specifies the dependent variable by a single character that is used within the input file to refer to a particular PDE system. Note that there can be multiple PDE blocks specified with the same or different types, but their depvar must be unique. This depvar then may be used downstream of its specification to configure a scalar component of a given PDE to be used as a refinement variable as error indicator for adaptive refinement. Another way depvars are used is to specify coupling among PDEs. An example in Walker, which uses the same coupling mechanism in its input files and where systems of stochastic differential equations (SDE) are specified analogously as PDEs in Inciter, is given in glm_homogeneous_shear.q, which couples three systems of SDEs to integrate a joint Fokker-Planck equation of particle positions, velocities, and turbulent kinetic energy dissipation rate in a homogeneous shear-driven turbulent flow.

2. Add new keywords

To add a new PDE type, we start by adding a new keyword that will be used to designate it in the input file. Since the new PDE is for methods specifialized for multi-material compressible flow, we add a new keyword multimat. We also add a new keyword for a new physics policy, veleq, to designate an important feature of the new method to be implemented, velocity equilibrium, assuming a single velocity field for all materials, while solving for separate material fields for mass and internal energy. Another new keyword we add is nmat which will be used to specify the number of materials for a multi-material flow. New keyword definitions with their documentation go to src/Control/Keywords.h.

Control/Keywords.h

git diff src/Control/Keywords.h
diff --git a/src/Control/Keywords.h b/src/Control/Keywords.h
index b3fbcebd..c57c65ce 100644
--- a/src/Control/Keywords.h
+++ b/src/Control/Keywords.h
@@ -1727,6 +1727,23 @@ struct ncomp_info {
 };
 using ncomp = keyword< ncomp_info,  TAOCPP_PEGTL_STRING("ncomp") >;
 
+struct nmat_info {
+  static std::string name() { return "nmat"; }
+  static std::string shortDescription() { return
+    "Set number of materials for a system of differential equations"; }
+  static std::string longDescription() { return
+    R"(This keyword is used to specify the number of materials, e.g., for
+    multi-material flow, see also the keyword 'multimat' and 'veleq'.)";
+  }
+  struct expect {
+    using type = std::size_t;
+    static constexpr type lower = 1;
+    static std::string description() { return "uint"; }
+  };
+};
+using nmat = keyword< nmat_info,  TAOCPP_PEGTL_STRING("nmat") >;
+
 struct ttyi_info {
   static std::string name() { return "ttyi"; }
   static std::string shortDescription() { return
@@ -3769,6 +3786,23 @@ struct euler_info {
 };
 using euler = keyword< euler_info, TAOCPP_PEGTL_STRING("euler") >;
 
+struct veleq_info {
+  using code = Code< V >;
+  static std::string name() { return "Velocity equilibrium"; }
+  static std::string shortDescription() { return "Specify the multi-material "
+    " compressible flow with velocity equilibrium as physics configuration"; }
+  static std::string longDescription() { return
+    R"(This keyword is used to select a compressible flow algorithm as physics
+    configuration designed for multiple materials assuming velocity equailibrium
+    (single velocity). Example: "multimat physics veleq end")";
+    }
+  struct expect {
+    static std::string description() { return "string"; }
+  };
+};
+using veleq = keyword< veleq_info, TAOCPP_PEGTL_STRING("veleq") >;
+
 struct advection_info {
   using code = Code< A >;
   static std::string name() { return "Advection"; }
@@ -4368,6 +4402,41 @@ struct compflow_info {
 };
 using compflow = keyword< compflow_info, TAOCPP_PEGTL_STRING("compflow") >;
 
+struct multimat_info {
+  static std::string name() { return "Compressible multi-material flow"; }
+  static std::string shortDescription() { return "Start configuration block "
+    "for the multi-material compressible flow equations"; }
+  static std::string longDescription() { return
+    R"(This keyword is used to introduce the multimat ... end block,
+    used to specify the configuration for a system of partial differential
+    equations, governing multi-material compressible fluid flow. Keywords
+    allowed in a multimat ... end block: )" + std::string("\'")
+    + depvar::string()+ "\', \'"
+    + physics::string() + "\', \'"
+    + problem::string() + "\', \'"
+    + material::string() + "\', \'"
+    + nmat::string() + "\', \'"
+    + pde_alpha::string() + "\', \'"
+    + pde_p0::string() + "\', \'"
+    + pde_betax::string() + "\', \'"
+    + pde_betay::string() + "\', \'"
+    + pde_betaz::string() + "\', \'"
+    + pde_beta::string() + "\', \'"
+    + pde_r0::string() + "\', \'"
+    + pde_ce::string() + "\', \'"
+    + pde_kappa::string() + "\', \'"
+    + bc_dirichlet::string() + "\', \'"
+    + bc_sym::string() + "\', \'"
+    + bc_inlet::string() + "\', \'"
+    + bc_outlet::string() + "\', \'"
+    + bc_extrapolate::string() + "\'. "
+    + R"(For an example multimat ... end block, see
+      doc/html/inicter_example_multimat.html.)";
+  }
+};
+using multimat = keyword< multimat_info, TAOCPP_PEGTL_STRING("multimat") >;
+
 struct rcb_info {
   static std::string name() { return "recursive coordinate bisection"; }
   static std::string shortDescription() { return

For more details on the various (mandatory and optional) fields of a keyword specification struct see the documentation in src/Control/Keywords.h.

We also add the new keywords to inciter's grammar's keywords pool:

Control/Inciter/InputDeck/InputDeck.h

git diff src/Control/Inciter/InputDeck/InputDeck.h 
diff --git a/src/Control/Inciter/InputDeck/InputDeck.h b/src/Control/Inciter/InputDeck/InputDeck.h
index 39eff786..155c47dc 100644
--- a/src/Control/Inciter/InputDeck/InputDeck.h
+++ b/src/Control/Inciter/InputDeck/InputDeck.h
@@ -70,12 +70,14 @@ class InputDeck :
                                    kw::phg,
                                    kw::inciter,
                                    kw::ncomp,
+                                   kw::nmat,
                                    kw::pde_diffusivity,
                                    kw::pde_lambda,
                                    kw::pde_u0,
                                    kw::bc_dirichlet,
                                    kw::sideset,
                                    kw::compflow,
+                                   kw::multimat,
                                    kw::ic,
                                    kw::txt_float_format,
                                    kw::txt_float_default,
@@ -95,6 +97,7 @@ class InputDeck :
                                    kw::advdiff,
                                    kw::navierstokes,
                                    kw::euler,
+                                   kw::veleq,
                                    kw::user_defined,
                                    kw::vortical_flow,
                                    kw::pde_alpha,

This is required so that the compiler can generate a database containing the help for all the keywords in the grammar understood by inciter's control file parser. The above changes not only add the keyword but also some documentation that gets displayed when passing the -C or -H command line arguments to the inciter executable, so quick help is available at the user's fingertips:

$ inciter -C
inciter Control File Keywords:
...
   multimat              Start configuration block for the multi-material compressible flow equations
        navierstokes     string Specify the Navier-Stokes (viscous) compressible flow physics configuration
               ncomp       uint Set number of scalar components for a system of differential equations
    nl_energy_growth     string Select the nonlinear energy growth test problem 
                nmat       uint Set number of materials for a system of differential equations
...
$ inciter -H multimat
inciter control file keyword 'multimat'

   Start configuration block for the multi-material compressible flow equations

   This keyword is used to introduce the multimat ... end block, used to specify
   the configuration for a system of partial differential equations, governing
   multi-material compressible fluid flow. Keywords allowed in a multimat ... end
   block: 'depvar', 'physics', 'problem', 'material', 'nmat', 'alpha', 'p0',
   'betax', 'betay', 'betaz', 'beta', 'r0', 'ce', 'kappa', 'bc_dirichlet',
   'bc_sym', 'bc_inlet', 'bc_outlet', 'bc_extrapolate'.For an example multimat ...
   end block, see doc/html/inicter_example_multimat.html.

3. Add new option switches

Next is to add a new state to the existing PDE option switch, ctr::PDE. This "option switch" is really only a fancy enum, used to store the user's choice of the PDE type in a type-safe manner, after parsing a PDE-block, e.g., compflow ... end, in the control file. This fancy enum is an option switch because it inherits from tk::Toggle, defined in Control/Toggle.h, which is a generic switch (or option), that helps associating enum values to keywords and querying one based on the other. Extending the existing PDE option switch is done by extending the list of PDE types in src/Control/Inciter/Options/PDE.h:

Control/Inciter/Options/PDE.h

git diff src/Control/Inciter/Options/PDE.h
diff --git a/src/Control/Inciter/Options/PDE.h b/src/Control/Inciter/Options/PDE.h
index 2b4af067..84b6d3ff 100644
--- a/src/Control/Inciter/Options/PDE.h
+++ b/src/Control/Inciter/Options/PDE.h
@@ -22,7 +22,8 @@ namespace ctr {
 
 //! Differential equation types
 enum class PDEType : uint8_t { TRANSPORT=0,
-                               COMPFLOW };
+                               COMPFLOW,
+                               MULTIMAT };
 
 //! Pack/Unpack: forward overload to generic enum class packer
 inline void operator|( PUP::er& p, PDEType& e ) { PUP::pup( p, e ); }
@@ -40,6 +41,7 @@ class PDE : public tk::Toggle< PDEType > {
     // List valid expected choices to make them also available at compile-time
     using keywords = brigand::list< kw::transport
                                   , kw::compflow
+                                  , kw::multimat
                                   >;
 
     //! Constructor: pass associations references to base, which will handle
@@ -48,10 +50,12 @@ class PDE : public tk::Toggle< PDEType > {
       tk::Toggle< PDEType >( "Partial differential equation",
         //! Enums -> names
         { { PDEType::TRANSPORT, kw::transport::name() },
-          { PDEType::COMPFLOW, kw::compflow::name() } },
+          { PDEType::COMPFLOW, kw::compflow::name() },
+          { PDEType::MULTIMAT, kw::multimat::name() } },
         //! keywords -> Enums
         { { kw::transport::string(), PDEType::TRANSPORT },
-          { kw::compflow::string(), PDEType::COMPFLOW } } ) {}
+          { kw::compflow::string(), PDEType::COMPFLOW },
+          { kw::multimat::string(), PDEType::MULTIMAT } } ) {}
 };
 
 } // ctr::

We also augment the switch used for selecting the physics policy:

Control/Inciter/Options/Physics.h

git diff src/Control/Inciter/Options/Physics.h
diff --git a/src/Control/Inciter/Options/Physics.h b/src/Control/Inciter/Options/Physics.h
index d24ccf1d..f59e85c1 100644
--- a/src/Control/Inciter/Options/Physics.h
+++ b/src/Control/Inciter/Options/Physics.h
@@ -23,7 +23,8 @@ namespace ctr {
 enum class PhysicsType : uint8_t { ADVECTION=0,
                                    ADVDIFF,
                                    EULER,
-                                   NAVIERSTOKES };
+                                   NAVIERSTOKES,
+                                   VELEQ };
 
 //! Pack/Unpack PhysicsType: forward overload to generic enum class packer
 inline void operator|( PUP::er& p, PhysicsType& e ) { PUP::pup( p, e ); }
@@ -35,8 +36,9 @@ class Physics : public tk::Toggle< PhysicsType > {
     //! Valid expected choices to make them also available at compile-time
     using keywords = brigand::list< kw::advection
                                   , kw::advdiff
-                                  , kw::navierstokes
                                   , kw::euler
+                                  , kw::navierstokes
+                                  , kw::veleq
                                   >;
 
     //! \brief Options constructor
@@ -50,12 +52,16 @@ class Physics : public tk::Toggle< PhysicsType > {
         { { PhysicsType::ADVECTION, kw::advection::name() },
           { PhysicsType::ADVDIFF, kw::advdiff::name() },
           { PhysicsType::EULER, kw::euler::name() },
-          { PhysicsType::NAVIERSTOKES, kw::navierstokes::name() } },
+          { PhysicsType::NAVIERSTOKES, kw::navierstokes::name() },
+          { PhysicsType::VELEQ, kw::veleq::name() }
+        },
         //! keywords -> Enums
         { { kw::advection::string(), PhysicsType::ADVECTION },
           { kw::advdiff::string(), PhysicsType::ADVDIFF },
           { kw::euler::string(), PhysicsType::EULER },
-          { kw::navierstokes::string(), PhysicsType::NAVIERSTOKES } } )
+          { kw::navierstokes::string(), PhysicsType::NAVIERSTOKES },
+          { kw::veleq::string(), PhysicsType::VELEQ } } )
     {
        brigand::for_each< keywords >( assertPolicyCodes() );
     }
@@ -89,6 +95,7 @@ class Physics : public tk::Toggle< PhysicsType > {
       , { PhysicsType::ADVDIFF, *kw::advdiff::code() }
       , { PhysicsType::EULER, *kw::euler::code() }
       , { PhysicsType::NAVIERSTOKES, *kw::navierstokes::code() }
+      , { PhysicsType::VELEQ, *kw::veleq::code() }
     };
 };

4. Add parsing/grammar for the new keywords

After adding the new keywords multimat, nmat, and veleq, we have to teach the input file parser to recognize them.

First we augment the data structures that store data parsed associated with the new keywords. First is the number of scalar components for each PDE type:

Control/Inciter/Components.h

git diff src/Control/Inciter/Components.h
diff --git a/src/Control/Inciter/Components.h b/src/Control/Inciter/Components.h
index 64f0696f..ef51a9e9 100644
--- a/src/Control/Inciter/Components.h
+++ b/src/Control/Inciter/Components.h
@@ -17,8 +17,9 @@ namespace ctr {
 
 //! Number of components of partial differential equations
 using ncomps = tk::ctr::ncomponents<
-  tag::compflow,     std::vector< tk::ctr::ncomp_type >
+  tag::compflow,     std::vector< tk::ctr::ncomp_type >,
+  tag::multimat,     std::vector< tk::ctr::ncomp_type >
 >;
 
 } // ctr::

tk::ctr::ncomponents is a tuple of vectors that stores the number of scalar components for each PDE type parsed. It is a vector for each PDE type because there can be multiple PDE systems of each type.

We also create a new tuple to store data during parsing the multimat...end block:

Control/Inciter/Types.h

git diff src/Control/Inciter/Types.h
index 4bc3bf64..65a3235f 100644
--- a/src/Control/Inciter/Types.h
+++ b/src/Control/Inciter/Types.h
@@ -184,10 +184,46 @@ using CompFlowPDEParameters = tk::tuple::tagged_tuple<
   tag::npar,          std::vector< kw::npar::info::expect::type >
 >;
 
+//! Compressible flow equation parameters storage
+using MultiMatPDEParameters = tk::tuple::tagged_tuple<
+  tag::depvar,        std::vector< char >,
+  tag::physics,       std::vector< PhysicsType >,
+  tag::problem,       std::vector< ProblemType >,
+  tag::bcdir,         std::vector< std::vector<
+                       kw::sideset::info::expect::type > >,
+  tag::bcsym,         std::vector< std::vector<
+                       kw::sideset::info::expect::type > >,
+  tag::bcinlet,       std::vector< std::vector<
+                        kw::sideset::info::expect::type > >,
+  tag::bcoutlet,      std::vector< std::vector<
+                        kw::sideset::info::expect::type > >,
+  tag::bcextrapolate, std::vector< std::vector<
+                         kw::sideset::info::expect::type > >,
+  //! Parameter vector (for specific, e.g., verification problems)
+  tag::alpha,         std::vector< kw::pde_alpha::info::expect::type >,
+  //! Parameter vector (for specific, e.g., verification problems)
+  tag::beta,          std::vector< kw::pde_beta::info::expect::type >,
+  //! Parameter vector (for specific, e.g., verification problems)
+  tag::p0,            std::vector< kw::pde_p0::info::expect::type >,
+  //! Material ID
+  tag::id,            std::vector< kw::id::info::expect::type >,
+  //! Ratio of spec heats
+  tag::gamma,         std::vector< kw::mat_gamma::info::expect::type >,
+  //! Dynamic viscosity
+  tag::mu,            std::vector< kw::mat_mu::info::expect::type >,
+  //! Spec. heat at const vol.
+  tag::cv,            std::vector< kw::mat_cv::info::expect::type >,
+  //! Heat conductivity
+  tag::k,             std::vector< kw::mat_k::info::expect::type >,
+  //! number of materials
+  tag::nmat,          std::vector< kw::nmat::info::expect::type >
+>;
+
 //! Parameters storage
 using parameters = tk::tuple::tagged_tuple<
-  tag::compflow,    CompFlowPDEParameters
+  tag::compflow,    CompFlowPDEParameters,
+  tag::multimat,    MultiMatPDEParameters
 >;

Note that in ctr::MultiMatPDEParameters every field is a vector, since multiple systems of this PDE type can be configured. Also note that this is where the dependent variable, physics policy (ctr::PhysicsType enum), the problem policy (ctr::ProblemType enum), boundary conditions, problem-configuration-specific parameters, and the number of materials are all stored. This is data that is either parsed from the user or can be computed immediately after and based on user input.

Next we add a couple of new tags (empty structs used as compile-time labels), tag::nmat and tag::multimat, behind which the parser stores the tuple ctr::MultiMatPDEParameters, and the number of materials, respectively,

Control/Tags.h

git diff src/Control/Tags.h
diff --git a/src/Control/Tags.h b/src/Control/Tags.h
index 188a86d0..9c0b7ec6 100644
--- a/src/Control/Tags.h
+++ b/src/Control/Tags.h
@@ -102,6 +102,7 @@ struct matched {};
 struct part {};
 struct centroid {};
 struct ncomp {};
+struct nmat {};
 struct tty {};
 struct dump {};
 struct plot {};
@@ -176,6 +177,7 @@ struct gid {};
 struct transport {};
 struct advection {};
 struct compflow {};
+struct multimat {};
 struct problem {};
 struct physics {};
 struct diffusivity {};

We are now ready to augment the grammar itself:

Control/Inciter/InputDeck/Grammar.h

git diff src/Control/Inciter/InputDeck/Grammar.h
diff --git a/src/Control/Inciter/InputDeck/Grammar.h b/src/Control/Inciter/InputDeck/Grammar.h
index 6b16b623..43a9bba7 100644
--- a/src/Control/Inciter/InputDeck/Grammar.h
+++ b/src/Control/Inciter/InputDeck/Grammar.h
@@ -35,8 +35,9 @@ namespace deck {
 
   //! \brief Number of registered equations
   //! \details Counts the number of parsed equation blocks during parsing.
-  static tk::tuple::tagged_tuple< tag::transport, std::size_t,
-                                  tag::compflow,  std::size_t > neq;
+  static tk::tuple::tagged_tuple< tag::transport, std::size_t,
+                                  tag::compflow,  std::size_t,
+                                  tag::multimat,  std::size_t > neq;
 
 } // ::deck
 } // ::inciter
@@ -201,6 +202,70 @@ namespace grm {
     }
   };

The above tagged tuple, deck::neq, stores the number of equations during parsing the inciter's input file. Since this is part of the state of inciter's parser, it is static and thus only visible locally to inciter's grammar and not the rest of the code. deck::neq is used to count the number of PDE systems parsed.

The following code is run immediately after the multimat...end block finished parsing, thus it is designed to do error checking and set/verify defaults for this particular PDE system.

+  //! Rule used to trigger action
+  template< class eq > struct check_multimat : pegtl::success {};
+  //! \brief Set defaults and do error checking on the multimaterial
+  //!    compressible flow equation block
+  //! \details This is error checking that only the multimaterial compressible
+  //!   flow equation block must satisfy. Besides error checking we also set
+  //!   defaults here as this block is called when parsing of a
+  //!   multimat...end block has just finished.
+  template< class eq >
+  struct action< check_multimat< eq > > {
+    template< typename Input, typename Stack >
+    static void apply( const Input& in, Stack& stack ) {
+      using inciter::deck::neq;
+
+      // Error out if no dependent variable has been selected
+      auto& depvar = stack.template get< tag::param, eq, tag::depvar >();
+      if (depvar.empty() || depvar.size() != neq.get< eq >())
+        Message< Stack, ERROR, MsgKey::NODEPVAR >( stack, in );
+

First we ensure that a dependent variable has been associated. This is mandatory, as this identifies the PDE system during parsing. If depvar is not specified in this equation block or the vector that stores the dependent variables for this equation block does not have the same number of entries as many PDEs have been configured so far, that means the dependent variable has not been specified, so the parser will issue an error when all parsing has been finished.

+      // If physics type is not given, default to 'veleq'
+      auto& physics = stack.template get< tag::param, eq, tag::physics >();
+      if (physics.empty() || physics.size() != neq.get< eq >())
+        physics.push_back( inciter::ctr::PhysicsType::VELEQ );
+

We also ensure that a mandatory physics policy has been selected. If the user did not select one, the default is the new veleq we are adding.

+      // Set number of scalar components based on number of materials
+      auto& nmat = stack.template get< tag::param, eq, tag::nmat >();
+      auto& ncomp = stack.template get< tag::component, eq >();
+      if (physics.back() == inciter::ctr::PhysicsType::VELEQ) {
+        // physics = veleq: m-material compressible flow
+        // scalar components: volfrac:m-1 + mass:m + momentum:3 + energy:m
+        // if nmat is unspecified, configure it be 2
+        if (nmat.empty() || nmat.size() != neq.get< eq >()) {
+          Message< Stack, WARNING, MsgKey::NONMAT >( stack, in );
+          nmat.push_back( 2 );
+        }
+        // set ncomp based on nmat
+        auto m = nmat.back();
+        ncomp.push_back( m-1 + m + 3 + m );
+      }
+

The number of materials is processed above. If the physics selected is the newly added velocity equilibrium and the user did not specify the number of materials, first we issue a warning, then we set the number of materials equal to 2. If the user has configured the number materials, we use it, and compute the number of scalar components, ncomp, for this PDE system, as appropriate for the algorithm.

+      // If problem type is not given, default to 'user_defined'
+      auto& problem = stack.template get< tag::param, eq, tag::problem >();
+      if (problem.empty() || problem.size() != neq.get< eq >())
+        problem.push_back( inciter::ctr::ProblemType::USER_DEFINED );
+      else if (problem.back() == inciter::ctr::ProblemType::VORTICAL_FLOW) {
+        const auto& alpha = stack.template get< tag::param, eq, tag::alpha >();
+        const auto& beta = stack.template get< tag::param, eq, tag::beta >();
+        const auto& p0 = stack.template get< tag::param, eq, tag::p0 >();
+        if ( alpha.size() != problem.size() ||
+             beta.size() != problem.size() ||
+             p0.size() != problem.size() )
+          Message< Stack, ERROR, MsgKey::VORTICAL_UNFINISHED >( stack, in );
+      }
+

The above is some basic sanity checking on the problem policy selected.

+      // Error check Dirichlet boundary condition block for all multimat
+      // configurations
+      for (const auto& s : stack.template get< tag::param, eq, tag::bcdir >())
+        if (s.empty())
+          Message< Stack, ERROR, MsgKey::BC_EMPTY >( stack, in );

Followed by some error checking for Dirichlet boundary conditions. Obviously, this section can and will grow in the future as more defaults and error checking will be added.

+    }
+  };
+
   //! Rule used to trigger action
   template< class Option, typename...tags >
   struct store_inciter_option : pegtl::success {};
@@ -565,6 +630,55 @@ namespace deck {
                                tag::bcextrapolate > >,
            check_errors< tag::compflow, tk::grm::check_compflow > > {};
 
+  //! compressible multi-material flow
+  struct multimat :
+         pegtl::if_must<
+           scan_eq< use< kw::multimat >, tag::multimat >,
+           tk::grm::block< use< kw::end >,
+                           tk::grm::policy< use,
+                                            use< kw::physics >,
+                                            ctr::Physics,
+                                            tag::multimat,
+                                            tag::physics >,
+                           tk::grm::policy< use,
+                                            use< kw::problem >,
+                                            ctr::Problem,
+                                            tag::multimat,
+                                            tag::problem >,
+                           tk::grm::depvar< use,
+                                            tag::multimat,
+                                            tag::depvar >,
+                           parameter< tag::multimat,
+                                      kw::nmat,
+                                      tag::nmat >,
+                           material_properties< tag::multimat >,
+                           parameter< tag::multimat,
+                                      kw::pde_alpha,
+                                      tag::alpha >,
+                           parameter< tag::multimat,
+                                      kw::pde_p0,
+                                      tag::p0 >,
+                           parameter< tag::multimat,
+                                      kw::pde_beta,
+                                      tag::beta >,
+                           bc< kw::bc_dirichlet,
+                               tag::multimat,
+                               tag::bcdir >,
+                           bc< kw::bc_sym,
+                               tag::multimat,
+                               tag::bcsym >,
+                           bc< kw::bc_inlet,
+                               tag::multimat,
+                               tag::bcinlet >,
+                           bc< kw::bc_outlet,
+                               tag::multimat,
+                               tag::bcoutlet >,
+                           bc< kw::bc_extrapolate,
+                               tag::multimat,
+                               tag::bcextrapolate > >,
+           check_errors< tag::multimat,
+                         tk::grm::check_multimat > > {};
+

The above is the code block that augments inciter's input file grammar with the multimat block. The rules of PEGTL (our parser library) are out of scope here, but even without knowing the details of how the parsing is done, one can surely infer what is (and what is not) parsed in this block. Examples of what are parsed: physics and problem policies, depvar, etc., and what is NOT parsed: number of (scalar) components, c.f., the definition of the transport equation grammar block, deck::transport, not listed, in Control/Inciter/InputDeck/Grammar.h.

   //! partitioning ... end block
   struct partitioning :
          pegtl::if_must<
@@ -580,7 +694,7 @@ namespace deck {
 
   //! equation types
   struct equations :
-         pegtl::sor< transport, compflow > {};
+         pegtl::sor< transport, compflow, multimat > {};
 
   //! refinement variable(s) (refvar) ... end block
   struct refvars :

We finally include (above) the new multimat grammar block among the equations recognized by the parser.

One last item is to augment the grammar messages data store to include our new warning on unspecified number of materials:

Control/CommonGrammar.h

git diff src/Control/CommonGrammar.h
diff --git a/src/Control/CommonGrammar.h b/src/Control/CommonGrammar.h
index 40c0b71d..c2de296e 100644
--- a/src/Control/CommonGrammar.h
+++ b/src/Control/CommonGrammar.h
@@ -98,6 +98,7 @@ namespace grm {
     WRONGGAUSSIAN,      //!< Wrong number of parameters configuring a PDF
     NEGATIVEPARAM,      //!< Negative variance given configuring a Gaussian
     NONCOMP,            //!< No number of components selected
+    NONMAT,             //!< No number of materials selected
     NORNG,              //!< No RNG selected
     NODT,               //!< No time-step-size policy selected
     MULDT,              //!< Multiple time-step-size policies selected
@@ -187,7 +188,10 @@ namespace grm {
       "variable to solve for." },
     { MsgKey::NONCOMP, "The number of components has not been specified in the "
       "block preceding this position. This is mandatory for the preceding "
       "block. Use the keyword 'ncomp' to specify the number of components." },
+    { MsgKey::NONMAT, "The number of materials has not been specified in the "
+      "block preceding this position. This is mandatory for the preceding "
+      "block. Use the keyword 'nmat' to specify the number of materials." },
     { MsgKey::NORNG, "The random number generator has not been specified in "
       "the block preceding this position. This is mandatory for the preceding "
       "block. Use the keyword 'rng' to specify the random number generator." },

With the above changes, inciter's parser can now recognize the following:

inciter
  ...
  scheme dg
  ...
  multimat    # configure our new PDE block for compressible multi-material flow
    depvar u    # assign dependent variable 'u' to this PDE
    physics veleq       # select 'euler' as compressible flow physics policy
    problem vortical_flow       # select 'vortical flow` as problem policy
    nmat 3      # configure 3 materials
    ...
  end
  ...
end

5. Augment the PDE stack

Next is to augment inciter's PDEStack, a factory from which the PDEs are instantiated according to user configuration.

Class PDEStack helps with configuring and instantiating PDE class objects from factories. There are separate factories for PDEs with continuous Galerkin (CG) and discontinuous Galerkin (DG) finite element discretizations, declared in PDE/PDEFactory.h: CGFactory and DGFactory. These factories are associative containers that assign PDE class constructors to a key that uniquely identifies a given PDE system configured to a given physics and problem policy.

First the constructors of all PDEs with all possible combinations of their valid policies are registered into their factory. This is done by PDEStack's constructor. Then when user input is known at runtime, i.e., we know what and how many PDE systems, with what discretizations, and with what policies the user wants to use, we instantiate the selected ones in PDEStack::selectedCG() and PDEStack::selectedDG().

The following diff augments PDEStack with our newly added PDE type.

PDE/PDEStack.C

git diff src/PDE/PDEStack.C
diff --git a/src/PDE/PDEStack.C b/src/PDE/PDEStack.C
index 70036aad..af1a2de7 100644
--- a/src/PDE/PDEStack.C
+++ b/src/PDE/PDEStack.C
@@ -17,6 +17,7 @@
 
 #include "ConfigureTransport.h"
 #include "ConfigureCompFlow.h"
+#include "ConfigureMultiMat.h"
 
 using inciter::PDEStack;
 
@@ -96,6 +97,7 @@ PDEStack::PDEStack() : m_cgfactory(), m_dgfactory(),
 {
   registerTransport( m_cgfactory, m_dgfactory, m_cgEqTypes, m_dgEqTypes );
   registerCompFlow( m_cgfactory, m_dgfactory, m_cgEqTypes, m_dgEqTypes );
+  registerMultiMat( m_dgfactory, m_dgEqTypes );
 }
 
 std::vector< inciter::CGPDE >
@@ -144,6 +146,8 @@ PDEStack::selectedDG() const
         pdes.push_back( createDG< tag::transport >( d, cnt ) );
       else if (d == ctr::PDEType::COMPFLOW)
         pdes.push_back( createDG< tag::compflow >( d, cnt ) );
+      else if (d == ctr::PDEType::MULTIMAT)
+        pdes.push_back( createDG< tag::multimat >( d, cnt ) );
       else Throw( "Can't find selected DGPDE" );
     }
 
@@ -169,6 +173,8 @@ PDEStack::info() const
       nfo.emplace_back( infoTransport( cnt ) );
     else if (d == ctr::PDEType::COMPFLOW)
       nfo.emplace_back( infoCompFlow( cnt ) );
+    else if (d == ctr::PDEType::MULTIMAT)
+      nfo.emplace_back( infoMultiMat( cnt ) );
     else Throw( "Can't find selected PDE" );
   }

registerMultiMat() will be defined in ConfigureMultiMat.[Ch].

6. Augment the PDE build object

Before adding the new (skeleton) files for the new PDE type (which will contain their actual implementation of the spatial discretization), we add a new object file to the PDE build object so that registration of the new PDE into PDEStack is complete:

PDE/CMakeLists.txt

git diff src/PDE/CMakeListst.txt
diff --git a/src/PDE/CMakeLists.txt b/src/PDE/CMakeLists.txt
index c2ae8d9b..76562c52 100644
--- a/src/PDE/CMakeLists.txt
+++ b/src/PDE/CMakeLists.txt
@@ -14,7 +14,8 @@ add_library(PDE
             Integrate/Riemann/RiemannFactory.C
             Limiter.C
             ConfigureTransport.C
-            ConfigureCompFlow.C)
+            ConfigureCompFlow.C
+            ConfigureMultiMat.C)
 
 target_include_directories(PDE PUBLIC
                            ${QUINOA_SOURCE_DIR}

7. New skeleton PDE classes

Finally, we add a set of new files that will implement the new PDE using DG. As an example, we also add a new physics and a couple of new problem policy classes: one for a user-defined problem and another one for a verification problem using a vortical flow setup.

Here are the new files and subdirectories under src/PDE/:

├── ConfigureMultiMat.C
├── ConfigureMultiMat.h
├── MultiMat
│   ├── DGMultiMat.h
│   ├── Physics
│   │   ├── DG.h
│   │   └── DGVelEq.h
│   ├── Problem
│   │   ├── UserDefined.h
│   │   └── VorticalFlow.h
│   └── Problem.h

Click on the filenames for the sources:

For registration into PDEStack and DGFactory:

Implementation and calls to the discretization of all DG operators for multi-material compressible flow physics:

Physics configuration:

Problem configuration: