From f4f0e599231ca03d8d5f1b286c5404d80e382745 Mon Sep 17 00:00:00 2001
From: Robert K <robertk@posteo.org>
Date: Mon, 5 Sep 2022 15:30:43 +0200
Subject: [PATCH] [feature][fluxes] Added HLL for potential temperature
formulation.
---
dune/fem-dg/models/defaultmodel.hh | 10 ++
dune/fem-dg/models/modelwrapper.hh | 28 ++++-
dune/fem-dg/operator/dg/dgpyoperator.hh | 1 +
.../operator/fluxes/advection/parameters.hh | 6 +-
.../operator/fluxes/euler/eulerflux_impl.hh | 110 ++++++++++--------
dune/fem-dg/operator/fluxes/euler/fluxes.hh | 19 +++
python/dune/femdg/patch.py | 10 ++
7 files changed, 130 insertions(+), 54 deletions(-)
diff --git a/dune/fem-dg/models/defaultmodel.hh b/dune/fem-dg/models/defaultmodel.hh
index ec7de0f6..eb7c38bf 100644
--- a/dune/fem-dg/models/defaultmodel.hh
+++ b/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$
*
diff --git a/dune/fem-dg/models/modelwrapper.hh b/dune/fem-dg/models/modelwrapper.hh
index 81293776..8a6cc25c 100644
--- a/dune/fem-dg/models/modelwrapper.hh
+++ b/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 );
}
diff --git a/dune/fem-dg/operator/dg/dgpyoperator.hh b/dune/fem-dg/operator/dg/dgpyoperator.hh
index 6b3aaf9e..5b465ab4 100644
--- a/dune/fem-dg/operator/dg/dgpyoperator.hh
+++ b/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>
diff --git a/dune/fem-dg/operator/fluxes/advection/parameters.hh b/dune/fem-dg/operator/fluxes/advection/parameters.hh
index c839bf0c..de4313d5 100644
--- a/dune/fem-dg/operator/fluxes/advection/parameters.hh
+++ b/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.
diff --git a/dune/fem-dg/operator/fluxes/euler/eulerflux_impl.hh b/dune/fem-dg/operator/fluxes/euler/eulerflux_impl.hh
index eea74349..4524100a 100644
--- a/dune/fem-dg/operator/fluxes/euler/eulerflux_impl.hh
+++ b/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 )
{}
diff --git a/dune/fem-dg/operator/fluxes/euler/fluxes.hh b/dune/fem-dg/operator/fluxes/euler/fluxes.hh
index 55ea1fbb..4561e1cc 100644
--- a/dune/fem-dg/operator/fluxes/euler/fluxes.hh
+++ b/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.
*
diff --git a/python/dune/femdg/patch.py b/python/dune/femdg/patch.py
index 04336da7..5a590e24 100644
--- a/python/dune/femdg/patch.py
+++ b/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:
--
GitLab