Merge branch 'walker' into develop

src/Control/Keywords.hpp

@@ -3132,7 +3132,7 @@ struct fluctuation_info {
     "Select fluctuation (as the dependent variable) to solve for"; }
   static std::string longDescription() { return
     R"(This keyword is used to select the fluctuation of a random variable as
-    what quantity to solve for, i.e., use asthe dependent variable, e.g., in a
+    what quantity to solve for, i.e., use as the dependent variable, e.g., in a
     position or velocity model for a stochastic particle. This configures how
     statistics must be interpreted.)"; }
   struct expect {
@@ -3142,6 +3142,23 @@ struct fluctuation_info {
 using fluctuation =
   keyword< fluctuation_info, TAOCPP_PEGTL_STRING("fluctuation") >;
 
+struct product_info {
+  static std::string name() { return "product"; }
+  static std::string shortDescription() { return
+    "Select product (as the dependent variable) to solve for"; }
+  static std::string longDescription() { return
+    R"(This keyword is used to select the product of multiple random variables
+    as what quantity to solve for, i.e., use as the dependent variable, e.g., in
+    a velocity model, solve for the product of the density and velocity, i.e.,
+    the momentum, for a stochastic particle. This configures how
+    statistics must be interpreted.)"; }
+  struct expect {
+    static std::string description() { return "string"; }
+  };
+};
+using product =
+  keyword< product_info, TAOCPP_PEGTL_STRING("product") >;
+
 struct solve_info {
   static std::string name() { return "solve for"; }
   static std::string shortDescription() { return
@@ -3154,7 +3171,8 @@ struct solve_info {
     static std::string description() { return "string"; }
     static std::string choices() {
       return '\'' + fullvar::string() + "\' | \'"
-                  + fluctuation::string() + '\'';
+                  + fluctuation::string() + "\' | \'"
+                  + product::string() + '\'';
     }
   };
 };

src/Control/Walker/InputDeck/Grammar.hpp

@@ -277,7 +277,7 @@ namespace grm {
   //!   coupled to eq among other DiffEqs of type coupledeq
   //! \tparam depvar_msg Error message key to use on missing coupled depvar
   //! \param[in] in Parser input
-  //! \param[in,out] stack Grammar stack to wrok with
+  //! \param[in,out] stack Grammar stack to work with
   //! \param[in] missing Error message key to use on missing coupled equation
   //!   if the coupling is required. Pass MsgKey::OPTIONAL as missing if the
   //!   coupling is optional.
@@ -359,6 +359,34 @@ namespace grm {
     }
   }
 
+  //! Query if equation 'eq' has been coupled to equation 'coupledeq'
+  //! \tparam eq Tag of the equation to query
+  //! \tparam coupledeq Tag of the equation that is potentially coupled to
+  //!   equation 'eq'
+  //! \param[in] stack Grammar stack to work with
+  //! \return True if equation 'eq' is coupled to equation 'coupledeq'
+  //! \note Always the eq system that is parsed last is interrogated.
+  template< typename eq, typename coupledeq, typename Stack >
+  static bool coupled( const Stack& stack ) {
+    return stack.template get< tag::param, eq, coupledeq >().back() != '-';
+  }
+
+  //! Query number of components of coupled equation
+  //! \tparam eq Tag of the equation that is coupled
+  //! \tparam coupledeq Tag of the equation that is coupled to equation 'eq'
+  //! \tparam id Tag to access the coupled equation 'eq' (relative) ids, see
+  //!   tk::grm::couple.
+  //! \param[in] stack Grammar stack to work with
+  //! \return Number of scalar components of coupled equation
+  //! \note Always the eq system that is parsed last is interrogated.
+  template< typename eq, typename coupledeq, typename id, typename Stack >
+  static std::size_t ncomp_coupled( const Stack& stack ) {
+    Assert( (coupled< eq, coupledeq >( stack )), "Eq must be coupled" );
+    // Query relative id of coupled eq
+    auto cid = stack.template get< tag::param, eq, id >().back();
+    return stack.template get< tag::component, coupledeq >().at( cid );
+  }
+
   //! Rule used to trigger action
   struct check_velocity : pegtl::success {};
   //! \brief Do error checking on the velocity eq block
@@ -386,6 +414,31 @@ namespace grm {
           stack.template get< tag::param, tag::velocity, tag::solve >();
         if (solve.size() != neq.get< tag::velocity >())
           Message< Stack, ERROR, MsgKey::NOSOLVE >( stack, in );
+
+        // Increase number of components by the number of particle densities if
+        // we solve for particle momentum, coupled to mass fractions. This is
+        // to allocate storage for particle velocity as variables derived from
+        // momentum.
+        if (!solve.empty() && solve.back() == walker::ctr::DepvarType::PRODUCT)
+        {
+          //! Error out if not coupled to mixmassfracbeta
+          if (!coupled< tag::velocity, tag::mixmassfracbeta >( stack )) {
+            Message< Stack, ERROR, MsgKey::MIXMASSFRACBETA_DEPVAR >(stack,in);
+          } else {
+            // access number of components of velocity eq just parsed
+            auto& ncomp = stack.template get< tag::component, tag::velocity >();
+            // query number of components of coupled mixmassfracbeta model
+            auto nc = ncomp_coupled< tag::velocity, tag::mixmassfracbeta,
+                                     tag::mixmassfracbeta_id >( stack );
+            // Augment storage of velocity equation, solving for momentum by
+            // the number of scalar components the coupled mixmassfracbeta mix
+            // model solves for. The magic number, 4, below is
+            // MixMassFractionBeta::NUMDERIVED + 1, and the 3 is the number of
+            // velocity components derived from momentum.
+            ncomp.back() += (nc/4)*3;
+          }
+        }
+
         // Set C0 = 2.1 if not specified
         auto& C0 = stack.template get< tag::param, tag::velocity, tag::c0 >();
         if (C0.size() != neq.get< tag::velocity >()) C0.push_back( 2.1 );

src/Control/Walker/InputDeck/InputDeck.hpp

@@ -218,6 +218,7 @@ class InputDeck : public tk::TaggedTuple< InputDeckMembers > {
                                  , kw::sde_com2
                                  , kw::fullvar
                                  , kw::fluctuation
+                                 , kw::product
                                  , kw::solve
                                  , kw::variant
                                  , kw::slm

src/Control/Walker/Options/Depvar.hpp

@@ -23,7 +23,8 @@ namespace ctr {
 
 //! Dependent variable options types
 enum class DepvarType : uint8_t { FULLVAR=0
-                                , FLUCTUATION };
+                                , FLUCTUATION
+                                , PRODUCT };
 
 //! Pack/Unpack DepvarType: forward overload to generic enum class packer
 inline void operator|( PUP::er& p, DepvarType& e ) { PUP::pup( p, e ); }
@@ -35,6 +36,7 @@ class Depvar : public tk::Toggle< DepvarType > {
     //! Valid expected choices to make them also available at compile-time
     using keywords = brigand::list< kw::fullvar
                                   , kw::fluctuation
+                                  , kw::product
                                   >;
 
     //! \brief Options constructor
@@ -43,13 +45,15 @@ class Depvar : public tk::Toggle< DepvarType > {
     explicit Depvar() :
       tk::Toggle< DepvarType >(
         //! Group, i.e., options, name
-        "Solve for",
+        "solve for",
         //! Enums -> names
         { { DepvarType::FULLVAR, kw::fullvar::name() },
-          { DepvarType::FLUCTUATION, kw::fluctuation::name() } },
+          { DepvarType::FLUCTUATION, kw::fluctuation::name() },
+          { DepvarType::PRODUCT, kw::product::name() } },
         //! keywords -> Enums
         { { kw::fullvar::string(), DepvarType::FULLVAR },
-          { kw::fluctuation::string(), DepvarType::FLUCTUATION } } )
+          { kw::fluctuation::string(), DepvarType::FLUCTUATION },
+          { kw::product::string(), DepvarType::PRODUCT } } )
     {}
 };
 

src/Control/Walker/Options/VelocityVariant.hpp

@@ -45,7 +45,7 @@ class VelocityVariant : public tk::Toggle< VelocityVariantType > {
     explicit VelocityVariant() :
       tk::Toggle< VelocityVariantType >(
         //! Group, i.e., options, name
-        "Model variant",
+        "model variant",
         //! Enums -> names
         { { VelocityVariantType::SLM, kw::slm::name() },
           { VelocityVariantType::GLM, kw::glm::name() } },

src/DiffEq/Beta/ConfigureMixMassFractionBeta.cpp

@@ -75,8 +75,12 @@ infoMixMassFractionBeta( std::map< ctr::DiffEqType, tk::ctr::ncomp_t >& cnt )
   nfo.emplace_back( "start offset in particle array", std::to_string(
     g_inputdeck.get< tag::component >().offset< eq >(c) ) );
   auto ncomp = g_inputdeck.get< tag::component, eq >()[c];
-  nfo.emplace_back( "number of components",
-    std::to_string( ncomp ) + " (" + std::to_string(ncomp/4) + "*4) " );
+
+  auto numderived =
+    MixMassFractionBeta<InitZero,MixMassFracBetaCoeffInstVel>::NUMDERIVED;
+  nfo.emplace_back( "number of components", std::to_string(ncomp) + " (=" +
+                    std::to_string(ncomp/(numderived+1)) + '*' +
+                    std::to_string(numderived+1) + ") " );
 
   coupledInfo< eq, tag::velocity, tag::velocity_id >
              ( c, "velocity", nfo );

src/DiffEq/Beta/MixMassFractionBeta.hpp

@@ -111,6 +111,13 @@ class MixMassFractionBeta {
     using eq = tag::mixmassfracbeta;
 
   public:
+    //! Number of derived variables computed by the SDE
+    //! \details Derived variables: density, specific volume, 1 - mass fraction
+    //! \warning If you change this, you must also change the constant in
+    //!   Velocity::m_mixmassfracbeta_ncomp in DiffEq/Velocity/Velocity.hpp.
+    //!   To see where, search for NUMDERIVED there.
+    static const std::size_t NUMDERIVED = 3;
+
     //! \brief Constructor
     //! \param[in] c Index specifying which system of mix mass-fraction beta
     //!   SDEs to construct. There can be multiple mixmassfracbeta ... end blocks
@@ -121,7 +128,8 @@ class MixMassFractionBeta {
       m_c( c ),
       m_depvar( g_inputdeck.get< tag::param, eq, tag::depvar >().at(c) ),
       // divide by the number of derived variables computed, see derived()
-      m_ncomp( g_inputdeck.get< tag::component >().get< eq >().at(c) / 4 ),
+      m_ncomp( g_inputdeck.get< tag::component >().get< eq >().at(c) /
+               (NUMDERIVED + 1) ),
       m_offset( g_inputdeck.get< tag::component >().offset< eq >(c) ),
       m_rng( g_rng.at( tk::ctr::raw(
         g_inputdeck.get< tag::param, eq, tag::rng >().at(c) ) ) ),
@@ -213,6 +221,8 @@ class MixMassFractionBeta {
                     m_velocity_solve, m_solve, m_ncomp, moments, m_bprime,
                     m_kprime, m_rho2, m_r, m_hts, m_hp, m_b, m_k, m_S, t );
 
+      const auto eps = std::numeric_limits< tk::real >::epsilon();
+
       // Advance particles
       const auto npar = particles.nunk();
       for (auto p=decltype(npar){0}; p<npar; ++p) {
@@ -222,7 +232,8 @@ class MixMassFractionBeta {
 
         // Access coupled particle velocity
         tk::real u = 0.0, v = 0.0, w = 0.0;
-        if (m_velocity_coupled) {
+        using std::abs;
+        if (abs(m_dY[0]) > eps || abs(m_dY[1]) > eps || abs(m_dY[2]) > eps) {
           u = particles( p, 0, m_velocity_offset );
           v = particles( p, 1, m_velocity_offset );
           w = particles( p, 2, m_velocity_offset );
@@ -259,7 +270,7 @@ class MixMassFractionBeta {
     const char m_dissipation_depvar;    //!< Depvar of coupled dissipation eq
     const ncomp_t m_dissipation_offset; //!< Offset of coupled dissipation eq
 
-    std::array<tk::real, 3> m_dY;       //! Prescribed mean scalar gradient
+    std::array< tk::real, 3 > m_dY;     //! Prescribed mean scalar gradient
 
     //! Coefficients
     std::vector< kw::sde_bprime::info::expect::type > m_bprime;

src/DiffEq/CoupledEq.hpp

@@ -131,13 +131,9 @@ std::size_t offset( std::size_t system )
 //!   tk::grm::couple in Control/Walker/InputDeck/Grammar.h
 //! \param[in] system Relative equation system id of equation 'eq'
 //! \return System offset of coupled equation in tk::Data array of all systems
-//! \note If equation 'eq' is not coupled to equation 'coupledeq', we return a
-//!   large number which will hopefully trigger some problems in client code if
-//!   used. This is used to indicate "no coupling" so that client code can still
-//!   call this function even for an equation that is not coupled, without
-//!   chaning client code, compared to equations that are coupled. In other
-//!   words, calling this function on a coupledeq equation that is not coupled
-//!   is not an error.
+//! \note If equation 'eq' is not coupled to equation 'coupledeq', we return the
+//!   id of the coupled equation. In other words, calling this function on a
+//!   coupledeq equation that is not coupled is not an error.
 // *****************************************************************************
 {
   static_assert( !std::is_same< eq, coupledeq >::value,
@@ -153,6 +149,34 @@ std::size_t offset( std::size_t system )
     return cid;
 }
 
+template< typename eq, typename coupledeq, typename id >
+std::size_t ncomp( std::size_t system )
+// *****************************************************************************
+//  Query number of components of coupled equation
+//! \tparam eq Tag of the equation that is coupled
+//! \tparam coupledeq Tag of the equation that is coupled to equation 'eq'
+//! \tparam id Tag to access the coupled equation 'eq' (relative) ids, see
+//!   tk::grm::couple in Control/Walker/InputDeck/Grammar.h
+//! \param[in] system Relative equation system id of equation 'eq'
+//! \return Number of scalar components of coupled equation
+//! \note If equation 'eq' is not coupled to equation 'coupledeq', we return the
+//!   id of the coupled equation. In other words, calling this function on a
+//!   coupledeq equation that is not coupled is not an error.
+// *****************************************************************************
+{
+  static_assert( !std::is_same< eq, coupledeq >::value,
+                 "Eq and coupled eq must differ" );
+
+  // Query relative id of coupled eq
+  auto cid = system_id< eq, coupledeq, id >( system );
+
+  // Return number of scalar components of coupled eq
+  if (coupled< eq, coupledeq >( system ))
+    return g_inputdeck.get< tag::component >().get< coupledeq >().at( cid );
+  else
+    return cid;
+}
+
 } // walker::
 
 #endif // CoupledEq_h

src/DiffEq/Velocity/ConfigureVelocity.cpp

@@ -65,10 +65,23 @@ infoVelocity( std::map< ctr::DiffEqType, tk::ctr::ncomp_t >& cnt )
 
   nfo.emplace_back( ctr::DiffEq().name( ctr::DiffEqType::VELOCITY ), "" );
 
+  auto solve = g_inputdeck.get< tag::param, eq, tag::solve >()[c];
+
   nfo.emplace_back( "start offset in particle array", std::to_string(
     g_inputdeck.get< tag::component >().offset< eq >(c) ) );
   auto ncomp = g_inputdeck.get< tag::component, eq >()[c];
-  nfo.emplace_back( "number of components", std::to_string( ncomp ) );
+
+  if (solve == ctr::DepvarType::PRODUCT &&
+      coupled< eq, tag::mixmassfracbeta >(c))
+  {
+    auto numderived =
+      Velocity<InitZero,VelocityCoeffStationary>(c).numderived();
+    nfo.emplace_back( "number of components", std::to_string(ncomp) + " (=" +
+                      std::to_string(ncomp - numderived) + '+' +
+                      std::to_string(numderived) + ") " );
+  } else {
+    nfo.emplace_back( "number of components", std::to_string(ncomp) );
+  }
 
   coupledInfo< eq, tag::position, tag::position_id >
              ( c, "position", nfo );
@@ -85,7 +98,6 @@ infoVelocity( std::map< ctr::DiffEqType, tk::ctr::ncomp_t >& cnt )
   auto cp = g_inputdeck.get< tag::param, eq, tag::coeffpolicy >()[c];
   nfo.emplace_back( "coefficients policy", ctr::CoeffPolicy().name( cp ) );
 
-  auto solve = g_inputdeck.get< tag::param, eq, tag::solve >()[c];
   auto dv = ctr::Depvar();
   nfo.emplace_back( dv.group(), dv.name(solve) );