[feature][fluxes] Added HLL for potential temperature formulation.

dune/fem-dg/models/defaultmodel.hh

@@ -218,6 +218,16 @@ namespace Fem
       globalMax = -1e308;
     }
 
+    /**
+     *  \brief compute pressure and temperature for compressible flow
+     */
+    template <class State>
+    void pressureTemperatur( const State& u, double& p, double& T ) const
+    {
+      p = 0.0;
+      T = 0.0;
+    }
+
     /**
      * \brief returns the mass factor term \f$ R \f$
      *

dune/fem-dg/models/modelwrapper.hh

@@ -226,9 +226,9 @@ namespace Fem
     template <class Entity>
     void setEntity( const Entity& entity ) const
     {
-      if( hasAdvection )
+      if constexpr ( hasAdvection )
         advection().init( entity );
-      if( hasDiffusion )
+      if constexpr ( hasDiffusion )
         diffusion().init( entity );
     }
 
@@ -249,6 +249,28 @@ namespace Fem
       return (min - max).infinity_norm() < 1e-10;
     }
 
+    template <class State> // double[dim+2]
+    void pressureTemperature( const State& u, double& p, double& T ) const
+    {
+      if constexpr ( hasAdvection )
+      {
+        const auto pT = advection().pressureTemperature( u );
+        p = pT[0];
+        T = pT[1];
+        return ;
+      }
+
+      if constexpr ( hasDiffusion )
+      {
+        const auto pT = diffusion().pressureTemperature( u );
+        p = pT[0];
+        T = pT[1];
+        return ;
+      }
+
+      p = 0.0;
+      T = 0.0;
+    }
 
     template <class LocalEvaluation>
     inline double stiffSource( const LocalEvaluation& local,
@@ -288,7 +310,7 @@ namespace Fem
                             const RangeType& u,
                             RangeType& maxValue) const
     {
-      if( hasDiffusion )
+      if constexpr ( hasDiffusion )
       {
         maxValue = diffusion().diffusionTimeStep( local.entity(), local.quadraturePoint(), u );
       }

dune/fem-dg/operator/dg/dgpyoperator.hh

@@ -18,6 +18,7 @@
 // dune-fem-dg includes
 #include <dune/fem-dg/algorithm/evolution.hh>
 #include <dune/fem-dg/operator/fluxes/advection/fluxes.hh>
+#include <dune/fem-dg/operator/fluxes/euler/fluxes.hh>
 #include <dune/fem-dg/operator/dg/operatortraits.hh>
 #include <dune/fem-dg/operator/dg/primaloperator.hh>
 #include <dune/fem-dg/models/modelwrapper.hh>

dune/fem-dg/operator/fluxes/advection/parameters.hh

@@ -38,6 +38,8 @@ namespace Fem
       euler_llf,
       //! the Harten, Lax and van Leer (HLL) flux
       euler_hll,
+      //! the Harten, Lax and van Leer (HLL) flux with potential temperature formulation
+      euler_hll_pt,
       //! the HLLC flux
       euler_hllc,
       //! general flux: Parameter selection is done via parameter file!
@@ -55,12 +57,12 @@ namespace Fem
 
     //! Contains all known enums for advection fluxes which can be chosen via parameter file.
     const Enum        _enums[] = { Enum::none, Enum::upwind, Enum::llf,
-                                   Enum::euler_llf, Enum::euler_hll, Enum::euler_hllc,
+                                   Enum::euler_llf, Enum::euler_hll, Enum::euler_hll_pt, Enum::euler_hllc,
                                    Enum::mhd_dw, Enum::mhd_hllem
                                    };
     //! Contains all known names of advection fluxes which can be chosen via parameter file.
     const std::string _strings[] = { "NONE", "UPWIND" , "LLF",
-                                     "EULER-LLF", "EULER-HLL" , "EULER-HLLC",
+                                     "EULER-LLF", "EULER-HLL" , "EULER-HLL-PT", "EULER-HLLC",
                                      "MHD-DW", "MHD-HLLEM"
                                      };
     //! Number of known advection fluxes which can be chosen via parameter file.

dune/fem-dg/operator/fluxes/euler/eulerflux_impl.hh

@@ -24,7 +24,7 @@ namespace EulerNumFlux
   //
   ////////////////////////////////////////////////////////
 
-  typedef enum {LLF, HLL, HLLC} EulerFluxType;
+  typedef enum {LLF, HLL, HLL_PT, HLLC} EulerFluxType;
   const EulerFluxType std_flux_type = HLL;
 
   /**
@@ -74,11 +74,14 @@ namespace EulerNumFlux
     FieldType num_flux(const FieldType Uj[dim+2], const FieldType Un[dim+2],
                        const FieldType normal[dim], FieldType gj[dim+2]) const
     {
-      if (flux_type == LLF) return num_flux_LLF(Uj, Un, normal, gj);
+      if constexpr (flux_type == LLF) return num_flux_LLF(Uj, Un, normal, gj);
 
-      if (flux_type == HLL) return num_flux_HLL(Uj, Un, normal, gj);
+      if constexpr (flux_type == HLL) return num_flux_HLL(Uj, Un, normal, gj, std::false_type());
 
-      if (flux_type == HLLC) return num_flux_HLLC(Uj, Un, normal, gj);
+      // potential temperature
+      if constexpr (flux_type == HLL_PT) return num_flux_HLL(Uj, Un, normal, gj, std::true_type());
+
+      if constexpr (flux_type == HLLC) return num_flux_HLLC(Uj, Un, normal, gj);
 
       DUNE_THROW( Dune::NotImplemented, "Numerical flux not implemented" );
     }
@@ -142,28 +145,34 @@ namespace EulerNumFlux
     }
 
 
+    template <bool pottemp>
     FieldType num_flux_HLL(const FieldType Uj[dim+2], const FieldType Un[dim+2],
-                           const FieldType normal[dim], FieldType gj[dim+2]) const
+                           const FieldType normal[dim], FieldType gj[dim+2],
+                           const std::integral_constant<bool, pottemp> ) const
     {
+      static constexpr bool potentialTemperature = pottemp;
+
       const FieldType rhoj = Uj[0];
       FieldType Ej = Uj[dim+1];
       const FieldType rhon = Un[0];
       FieldType En = Un[dim+1];
-#ifdef POTTEMP
-      FieldType pressj, tempj;
-      FieldType pressn, tempn;
-      model_.pressAndTemp( Uj, pressj, tempj );
-      model_.pressAndTemp( Un, pressn, tempn );
-
-      FieldType Ekinj = 0;
-      FieldType Ekinn = 0;
-      for(int i=1; i<dim+1; i++){
-        Ekinj += (0.5/rhoj) * Uj[i] * Uj[i];
-        Ekinn += (0.5/rhon) * Un[i] * Un[i];
+
+      if constexpr (potentialTemperature)
+      {
+        FieldType pressj, tempj;
+        FieldType pressn, tempn;
+        model_.pressureTemperature( Uj, pressj, tempj );
+        model_.pressureTemperature( Un, pressn, tempn );
+
+        FieldType Ekinj = 0;
+        FieldType Ekinn = 0;
+        for(int i=1; i<dim+1; i++){
+          Ekinj += (0.5/rhoj) * Uj[i] * Uj[i];
+          Ekinn += (0.5/rhon) * Un[i] * Un[i];
+        }
+        Ej = pressj/(_gamma-1.) + Ekinj;
+        En = pressn/(_gamma-1.) + Ekinn;
       }
-      Ej = pressj/(_gamma-1.) + Ekinj;
-      En = pressn/(_gamma-1.) + Ekinn;
-#endif
 
       FieldType rho_uj[dim], rho_un[dim], uj[dim], un[dim];
       FieldType Ekin2j=0.0, Ekin2n=0.0;
@@ -204,11 +213,10 @@ namespace EulerNumFlux
           for(int i=0; i<dim; i++) guj[i] = rho_uj[i]*uj[0];
           guj[0] += pj;
 
-#ifdef POTTEMP
-          gj[dim+1] = Uj[dim+1]*uj[0];
-#else
-          gj[dim+1] = (Ej+pj)*uj[0];
-#endif
+          if constexpr (potentialTemperature)
+            gj[dim+1] = Uj[dim+1]*uj[0];
+          else
+            gj[dim+1] = (Ej+pj)*uj[0];
         }
         else{
           const FieldType tmp1 = sj * sn;
@@ -220,15 +228,18 @@ namespace EulerNumFlux
           }
           guj[0] += tmp2 * (sn*pj - sj*pn);
 
-#ifdef POTTEMP
-          const FieldType Etmpj = Uj[dim+1]*uj[0];
-          const FieldType Etmpn = Un[dim+1]*un[0];
-          gj[dim+1] = tmp2 * (sn*Etmpj-sj*Etmpn + tmp1*(Un[dim+1] - Uj[dim+1]));
-#else
-          const FieldType Etmpj = (Ej+pj)*uj[0];
-          const FieldType Etmpn = (En+pn)*un[0];
-          gj[dim+1] = tmp2 * (sn*Etmpj-sj*Etmpn + tmp1*(En - Ej));
-#endif
+          if constexpr(potentialTemperature)
+          {
+            const FieldType Etmpj = Uj[dim+1]*uj[0];
+            const FieldType Etmpn = Un[dim+1]*un[0];
+            gj[dim+1] = tmp2 * (sn*Etmpj-sj*Etmpn + tmp1*(Un[dim+1] - Uj[dim+1]));
+          }
+          else
+          {
+            const FieldType Etmpj = (Ej+pj)*uj[0];
+            const FieldType Etmpn = (En+pn)*un[0];
+            gj[dim+1] = tmp2 * (sn*Etmpj-sj*Etmpn + tmp1*(En - Ej));
+          }
         }
       }
       else{
@@ -238,11 +249,10 @@ namespace EulerNumFlux
           for(int i=0; i<dim; i++) guj[i] = rho_un[i]*un[0];
           guj[0] += pn;
 
-#ifdef POTTEMP
-          gj[dim+1] = Uj[dim+1]*un[0];
-#else
-          gj[dim+1] = (En+pn)*un[0];
-#endif
+          if constexpr(potentialTemperature)
+            gj[dim+1] = Uj[dim+1]*un[0];
+          else
+            gj[dim+1] = (En+pn)*un[0];
         }
         else{
           const FieldType tmp1 = sj * sn;
@@ -254,16 +264,18 @@ namespace EulerNumFlux
           }
           guj[0] += tmp2 * (sn*pj - sj*pn);
 
-#ifdef POTTEMP
-          const FieldType Etmpj = Uj[dim+1]*uj[0];
-          const FieldType Etmpn = Un[dim+1]*un[0];
-          gj[dim+1] = tmp2 * (sn*Etmpj-sj*Etmpn + tmp1*(Un[dim+1] - Uj[dim+1]));
-#else
-          const FieldType Etmpj = (Ej+pj)*uj[0];
-          const FieldType Etmpn = (En+pn)*un[0];
-          gj[dim+1] = tmp2 * (sn*Etmpj-sj*Etmpn + tmp1*(En - Ej));
-#endif
-
+          if constexpr(potentialTemperature)
+          {
+            const FieldType Etmpj = Uj[dim+1]*uj[0];
+            const FieldType Etmpn = Un[dim+1]*un[0];
+            gj[dim+1] = tmp2 * (sn*Etmpj-sj*Etmpn + tmp1*(Un[dim+1] - Uj[dim+1]));
+          }
+          else
+          {
+            const FieldType Etmpj = (Ej+pj)*uj[0];
+            const FieldType Etmpn = (En+pn)*un[0];
+            gj[dim+1] = tmp2 * (sn*Etmpj-sj*Etmpn + tmp1*(En - Ej));
+          }
         }
       }
 
@@ -482,7 +494,7 @@ namespace Fem
     /**
      * \copydoc DGAdvectionFluxBase::DGAdvectionFluxBase()
      */
-    EulerFluxImpl (const ModelImp& mod )
+    EulerFluxImpl (const ModelImp& mod, const ParameterType& parameter = ParameterType() )
       : BaseType( mod ),
         numFlux_( mod )
     {}

dune/fem-dg/operator/fluxes/euler/fluxes.hh

@@ -53,6 +53,25 @@ namespace Fem
     static std::string name () { return "HLL (Dennis)"; }
   };
 
+  /**
+   * \brief class specialization for the HLLC flux.
+   *
+   * The purpose of this class is to allow the selection of an Euler flux
+   * via an enum given in AdvectionFlux::Enum.
+   */
+  template< class Model >
+  class DGAdvectionFlux< Model, AdvectionFlux::Enum::euler_hll_pt >
+    : public EulerFluxImpl< Model, EulerNumFlux::EulerFlux<Model,EulerNumFlux::EulerFluxType::HLL_PT > >
+  {
+    typedef EulerFluxImpl< Model, EulerNumFlux::EulerFlux<Model,EulerNumFlux::EulerFluxType::HLL_PT > > BaseType ;
+  public:
+    template< class ... Args>
+    DGAdvectionFlux(  Args&&... args )
+      : BaseType( std::forward<Args>(args)... )
+    {}
+    static std::string name () { return "HLL_PT (Dennis)"; }
+  };
+
   /**
    * \brief class specialization for the HLLC flux.
    *

python/dune/femdg/patch.py

@@ -188,6 +188,16 @@ def codeFemDg(self, *args, **kwargs):
     predefined.update( {t: arg_t} )
     self.predefineCoefficients(predefined,'x')
 
+    # pressure and temperature for potential temperature formulation
+    pressureTemperature = getattr(Model,"pressureTemperature",None)
+    if pressureTemperature is not None:
+        pressureTemperature = pressureTemperature(u)
+    self.generateMethod(code, pressureTemperature,
+          'Dune::FieldVector<double, 2>', 'pressureTemperature',
+          args=['const State &u'],
+          targs=['class State'], const=True, inline=True,
+          predefined=predefined)
+
     maxWaveSpeed = getattr(Model,"maxWaveSpeed",None)
     # check for deprecated maxLambda
     if maxWaveSpeed is None: