[feature][BGFixFlux] Implemented HLL for potential temperature and BGFix

model.

dune/fem-dg/models/defaultmodel.hh

@@ -222,12 +222,21 @@ namespace Fem
      *  \brief compute pressure and temperature for compressible flow
      */
     template <class State>
-    void pressureTemperatur( const State& u, double& p, double& T ) const
+    void pressureTemperature( const State& u, double& p, double& T ) const
     {
       p = 0.0;
       T = 0.0;
     }
 
+    /**
+     *  \brief compute background atmosphere
+     */
+    template <class Entity, class Point>
+    RangeType background( const Entity& entity, const Point& x ) const
+    {
+      return RangeType( 0 );
+    }
+
     /**
      * \brief returns the mass factor term \f$ R \f$
      *

dune/fem-dg/models/modelwrapper.hh

@@ -272,6 +272,22 @@ namespace Fem
       T = 0.0;
     }
 
+    template <class Entity, class Point>
+    RangeType background( const Entity& entity, const Point& x ) const
+    {
+      if constexpr ( hasAdvection )
+      {
+        return advection().background( time(), entity, x );
+      }
+
+      if constexpr ( hasDiffusion )
+      {
+        return diffusion().background( time(), entity, x );
+      }
+
+      return RangeType( 0 );
+    }
+
     template <class LocalEvaluation>
     inline double stiffSource( const LocalEvaluation& local,
                                const RangeType& u,

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

@@ -42,6 +42,9 @@ namespace Fem
       euler_hll_pt,
       //! the HLLC flux
       euler_hllc,
+      //! the Harten, Lax and van Leer (HLL) flux with potential temperature formulation
+      //! and background fix
+      euler_hll_bgfix,
       //! general flux: Parameter selection is done via parameter file!
       euler_general,
 

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

@@ -175,6 +175,90 @@ namespace Fem
 
   };
 
+
+  /**
+   * \brief class specialization for a general flux chosen by a parameter file.
+   *
+   * The purpose of this class is to allow the selection of an Euler flux
+   * via an enum given in AdvectionFlux::Enum.
+   */
+  template< class ModelImp >
+  class DGAdvectionFlux< ModelImp, AdvectionFlux::Enum::euler_hll_bgfix >
+    : public DGAdvectionFluxBase< ModelImp, AdvectionFluxParameters >
+  {
+    typedef DGAdvectionFluxBase< ModelImp, AdvectionFluxParameters >        BaseType;
+    typedef DGAdvectionFlux< ModelImp, AdvectionFlux::Enum::euler_hll_pt >  NumFluxImplType;
+
+    static const int dimRange = ModelImp::dimRange;
+    typedef typename ModelImp::DomainType         DomainType;
+    typedef typename ModelImp::RangeType          RangeType;
+    typedef typename ModelImp::JacobianRangeType  JacobianRangeType;
+    typedef typename ModelImp::FluxRangeType      FluxRangeType;
+    typedef typename ModelImp::FaceDomainType     FaceDomainType;
+
+  public:
+    typedef AdvectionFlux::Enum                   IdEnum;
+    typedef typename BaseType::ModelType          ModelType;
+    typedef typename BaseType::ParameterType      ParameterType;
+
+    using BaseType::model;
+
+    /**
+     * \copydoc DGAdvectionFluxBase::DGAdvectionFluxBase()
+     */
+    template< class ... Args>
+    DGAdvectionFlux(  Args&&... args )
+      : BaseType( std::forward<Args>(args)... ),
+        numFlux_( std::forward<Args>(args)... )
+    {}
+
+    /**
+     * \copydoc DGAdvectionFluxBase::name()
+     */
+    static std::string name () { return "AdvectionFlux - Euler (via parameter file)"; }
+
+    /**
+     * \copydoc DGAdvectionFluxBase::numericalFlux()
+     */
+    template< class LocalEvaluation >
+    inline double
+    numericalFlux( const LocalEvaluation& left,
+                   const LocalEvaluation& right,
+                   const RangeType& uLeft,
+                   const RangeType& uRight,
+                   const JacobianRangeType& jacLeft,
+                   const JacobianRangeType& jacRight,
+                   RangeType& gLeft,
+                   RangeType& gRight) const
+    {
+      // Assumption: background atmosphere is continuous
+      auto bgLeft  = model().background( left.entity(), left.quadraturePoint() );
+      auto bgRight = model().background( right.entity(), right.quadraturePoint() );
+
+      RangeType uLeftTot( uLeft );
+      RangeType uRightTot( uRight );
+
+      uLeftTot  += bgLeft;
+      uRightTot += bgRight;
+
+      RangeType gBgLeft, gBgRight;
+
+      // numerical flux
+      const double ws = numFlux_.numericalFlux(left, right, uLeftTot, uRightTot,
+          jacLeft, jacRight, gLeft, gRight );
+
+      // background flux
+      numFlux_.numericalFlux(left, right, bgLeft, bgRight, jacLeft, jacRight, gBgLeft, gBgRight );
+
+      gLeft  -= gBgLeft;
+      gRight -= gBgRight;
+
+      return ws;
+    }
+
+  protected:
+    NumFluxImplType numFlux_;
+  };
 }
 }
 

python/dune/femdg/patch.py

@@ -198,6 +198,17 @@ def codeFemDg(self, *args, **kwargs):
           targs=['class State'], const=True, inline=True,
           predefined=predefined)
 
+    # background atmosphere
+    background = getattr(Model,"background",None)
+    if background is not None:
+        background = background(t,x)
+    self.generateMethod(code, background,
+          'RRangeType', 'background',
+          args=['const double &t',
+                'const Entity &entity', 'const Point &x'],
+          targs=['class Entity, class Point'], const=True, inline=True,
+          predefined=None)
+
     maxWaveSpeed = getattr(Model,"maxWaveSpeed",None)
     # check for deprecated maxLambda
     if maxWaveSpeed is None: