How to add a new PDE type to Inciter
Contents
This page describes how to add a new PDE type to Inciter (Compressible flow solver).
Inciter (Compressible flow solver) 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.hpp
├── CMakeLists.txt
├── CompFlow
│ ├── CGCompFlow.hpp
│ ├── DGCompFlow.hpp
│ ├── Physics
│ │ ├── CGEuler.hpp
│ │ ├── CG.hpp
│ │ ├── CGNavierStokes.hpp
│ │ ├── DGEuler.hpp
│ │ ├── DG.hpp
│ │ └── DGNavierStokes.hpp
│ ├── Problem
│ │ ├── NLEnergyGrowth.hpp
│ │ ├── RayleighTaylor.hpp
│ │ ├── SedovBlastwave.hpp
│ │ ├── SodShocktube.hpp
│ │ ├── TaylorGreen.hpp
│ │ ├── UserDefined.hpp
│ │ └── VorticalFlow.hpp
│ └── Problem.hpp
├── ConfigureCompFlow.cpp
├── ConfigureCompFlow.hpp
├── ConfigureTransport.cpp
├── ConfigureTransport.hpp
├── DGPDE.hpp
├── FunctionPrototypes.hpp
├── Integrate
│ ├── Boundary.cpp
│ ├── Boundary.hpp
│ ├── Initialize.cpp
│ ├── Initialize.hpp
│ ├── Mass.cpp
│ ├── Mass.hpp
│ ├── Quadrature.cpp
│ ├── Quadrature.hpp
│ ├── Riemann
│ │ ├── HLLC.hpp
│ │ ├── LaxFriedrichs.hpp
│ │ ├── RiemannFactory.cpp
│ │ ├── RiemannFactory.hpp
│ │ ├── RiemannSolver.hpp
│ │ └── Upwind.hpp
│ ├── Source.cpp
│ ├── Source.hpp
│ ├── Surface.cpp
│ ├── Surface.hpp
│ ├── Volume.cpp
│ └── Volume.hpp
├── Limiter.cpp
├── Limiter.hpp
├── PDEFactory.hpp
├── PDEStack.cpp
├── PDEStack.hpp
└── Transport
├── CGTransport.hpp
├── DGTransport.hpp
├── Physics
│ ├── CGAdvDiff.hpp
│ ├── CGAdvection.hpp
│ ├── CG.hpp
│ ├── DGAdvection.hpp
│ └── DG.hpp
├── Problem
│ ├── CylAdvect.hpp
│ ├── GaussHump.hpp
│ ├── ShearDiff.hpp
│ └── SlotCyl.hpp
└── Problem.hpp
8 directories, 61 filesOnce 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 (Compressible flow solver) 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
...
endThis 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. 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.
- 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+transportwith, say, 12 scalar components, which then together will comprise a system of 17 coupled scalar equations. - The PDE must be specialized with a selected
physicstype. In the example above this iseuler, but this is only the simplest choice. Every PDE has a default physics which is selected if no physics is specified. - The PDE must be specialized with a selected
problemtype. In the example above this isvortical_flow, but this is only one choice. Every PDE has a default problem which is selected if no problem is specified. - The
depvarkeyword 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 theirdepvarmust be unique. Thisdepvarthen 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.
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. We also a new keyword for a problem type, interface_advection. New keyword definitions with their documentation go to src/Control/Keywords.h.pp
Control/Keywords.hpp
git diff src/Control/Keywords.hpp diff --git a/src/Control/Keywords.hpp b/src/Control/Keywords.hpp index b3fbcebd..c57c65ce 100644 --- a/src/Control/Keywords.hpp +++ b/src/Control/Keywords.hpp @@ -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 interface_advection_info { + using code = Code< I >; + static std::string name() { return "Interface advection"; } + static std::string shortDescription() { return + "Select the interface advection test problem "; } + static std::string longDescription() { return + R"(This keyword is used to select the interface advection test problem. The + purpose of this test problem is to test the well-balancedness of the + multi-material discretization and its interface capturing + capabilities. Example: "problem interface_advection".)"; } + struct expect { + static std::string description() { return "string"; } + }; +}; +using interface_advection = + keyword< interface_advection_info, + TAOCPP_PEGTL_STRING("interface_advection") >; + 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/
We also add the new keywords to inciter's grammar's keywords pool:
Control/Inciter/InputDeck/InputDeck.hpp
git diff src/Control/Inciter/InputDeck/InputDeck.hpp diff --git a/src/Control/Inciter/InputDeck/InputDeck.hpp b/src/Control/Inciter/InputDeck/InputDeck.hpp index 39eff786..155c47dc 100644 --- a/src/Control/Inciter/InputDeck/InputDeck.hpp +++ b/src/Control/Inciter/InputDeck/InputDeck.hpp @@ -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::interface_advection, 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::compflow ... end, in the control file. This fancy enum is an option switch because it inherits from tk::src/:
Control/Inciter/Options/PDE.hpp
git diff src/Control/Inciter/Options/PDE.hpp diff --git a/src/Control/Inciter/Options/PDE.hpp b/src/Control/Inciter/Options/PDE.hpp index 2b4af067..84b6d3ff 100644 --- a/src/Control/Inciter/Options/PDE.hpp +++ b/src/Control/Inciter/Options/PDE.hpp @@ -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 and problem list:
Control/Inciter/Options/Physics.hpp
git diff src/Control/Inciter/Options/Physics.hpp diff --git a/src/Control/Inciter/Options/Physics.hpp b/src/Control/Inciter/Options/Physics.hpp index d24ccf1d..f59e85c1 100644 --- a/src/Control/Inciter/Options/Physics.hpp +++ b/src/Control/Inciter/Options/Physics.hpp @@ -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() } }; };
Control/Inciter/Options/Problem.hpp
diff --git a/src/Control/Inciter/Options/Problem.hpp b/src/Control/Inciter/Options/Problem.hpp index c953f28f0..b21e2c891 100644 --- a/src/Control/Inciter/Options/Problem.hpp +++ b/src/Control/Inciter/Options/Problem.hpp @@ -34,7 +34,8 @@ enum class ProblemType : uint8_t { USER_DEFINED, CYL_ADVECT, SOD_SHOCKTUBE, ROTATED_SOD_SHOCKTUBE, - SEDOV_BLASTWAVE }; + SEDOV_BLASTWAVE, + INTERFACE_ADVECTION }; //! Pack/Unpack ProblemType: forward overload to generic enum class packer inline void operator|( PUP::er& p, ProblemType& e ) { PUP::pup( p, e ); } @@ -56,6 +57,7 @@ class Problem : public tk::Toggle< ProblemType > { , kw::sod_shocktube , kw::rotated_sod_shocktube , kw::sedov_blastwave + , kw::interface_advection >; //! \brief Options constructor @@ -78,7 +80,9 @@ class Problem : public tk::Toggle< ProblemType > { { ProblemType::SOD_SHOCKTUBE, kw::sod_shocktube::name() }, { ProblemType::ROTATED_SOD_SHOCKTUBE, kw::rotated_sod_shocktube::name() }, - { ProblemType::SEDOV_BLASTWAVE, kw::sedov_blastwave::name() } }, + { ProblemType::SEDOV_BLASTWAVE, kw::sedov_blastwave::name() }, + { ProblemType::INTERFACE_ADVECTION, + kw::interface_advection::name() } }, //! keywords -> Enums { { kw::user_defined::string(), ProblemType::USER_DEFINED }, { kw::shear_diff::string(), ProblemType::SHEAR_DIFF }, @@ -93,7 +97,9 @@ class Problem : public tk::Toggle< ProblemType > { { kw::rotated_sod_shocktube::string(), ProblemType::ROTATED_SOD_SHOCKTUBE }, { kw::sod_shocktube::string(), ProblemType::SOD_SHOCKTUBE }, - { kw::sedov_blastwave::string(), ProblemType::SEDOV_BLASTWAVE } } ) + { kw::sedov_blastwave::string(), ProblemType::SEDOV_BLASTWAVE }, + { kw::interface_advection::string(), + ProblemType::INTERFACE_ADVECTION } } ) { brigand::for_each< keywords >( assertPolicyCodes() ); } @@ -136,6 +142,7 @@ class Problem : public tk::Toggle< ProblemType > { , { ProblemType::ROTATED_SOD_SHOCKTUBE, *kw::rotated_sod_shocktube::code() } , { ProblemType::SEDOV_BLASTWAVE, *kw::sedov_blastwave::code() } + , { ProblemType::INTERFACE_ADVECTION, *kw::interface_advection::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.hpp
git diff src/Control/Inciter/Components.hpp diff --git a/src/Control/Inciter/Components.hpp b/src/Control/Inciter/Components.hpp index 64f0696f..ef51a9e9 100644 --- a/src/Control/Inciter/Components.hpp +++ b/src/Control/Inciter/Components.hpp @@ -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.hpp
git diff src/Control/Inciter/Types.hpp index 4bc3bf64..65a3235f 100644 --- a/src/Control/Inciter/Types.hpp +++ b/src/Control/Inciter/Types.hpp @@ -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::
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/Inciter/InputDeck/Grammar.hpp
git diff src/Control/Inciter/InputDeck/Grammar.hpp diff --git a/src/Control/Inciter/InputDeck/Grammar.hpp b/src/Control/Inciter/InputDeck/Grammar.hpp index 6b16b623..43a9bba7 100644 --- a/src/Control/Inciter/InputDeck/Grammar.hpp +++ b/src/Control/Inciter/InputDeck/Grammar.hpp @@ -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.hpp.
//! 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.hpp
git diff src/Control/CommonGrammar.hpp diff --git a/src/Control/CommonGrammar.hpp b/src/Control/CommonGrammar.hpp index 40c0b71d..c2de296e 100644 --- a/src/Control/CommonGrammar.hpp +++ b/src/Control/CommonGrammar.hpp @@ -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 interface_advection # select 'interface_advection` 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/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::
The following diff augments PDEStack with our newly added PDE type.
PDE/PDEStack.cpp
git diff src/PDE/PDEStack.cpp diff --git a/src/PDE/PDEStack.cpp b/src/PDE/PDEStack.cpp index 70036aad..af1a2de7 100644 --- a/src/PDE/PDEStack.cpp +++ b/src/PDE/PDEStack.cpp @@ -17,6 +17,7 @@ #include "ConfigureTransport.hpp" #include "ConfigureCompFlow.hpp" +#include "ConfigureMultiMat.hpp" 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.cpp Limiter.cpp ConfigureTransport.cpp - ConfigureCompFlow.cpp) + ConfigureCompFlow.cpp + ConfigureMultiMat.cpp) 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 interface advection.
Here are the new files and subdirectories under src/PDE/:
├── ConfigureMultiMat.cpp ├── ConfigureMultiMat.hpp ├── MultiMat │ ├── DGMultiMat.hpp │ ├── Physics │ │ ├── DG.hpp │ │ └── DGVelEq.hpp │ ├── Problem │ │ ├── UserDefined.hpp │ │ └── InterfaceAdvection.hpp │ └── 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:
- PDE/
MultiMat/ Physics/ DG.hpp - list of all MultiMat DG physics policies - PDE/
MultiMat/ Physics/ DGEuler.hpp - a new physics policy for inviscid multi-material policy
Problem configuration:
- PDE/
MultiMat/ Problem/ UserDefined.hpp - a new problem policy for multi-material compflow - PDE/
MultiMat/ Problem/ InterfaceAdvection.hpp - a new problem policy specialized for the interface advection verification problem - PDE/
MultiMat/ Problem.hpp - list of all MultiMat problem policies.