Skip to content
Snippets Groups Projects
Commit 62064bd5 authored by Robert K's avatar Robert K
Browse files

make field type exchangeable.

parent 84230eb3
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/fem-dg/test/navierstokes/Makefile.am
+ 9
− 7
View file @ 62064bd5
......@@ -18,11 +18,11 @@ SOURCES12 = $(BASEDIR)/main.cc $(BASEDIR)/main_1.cc $(BASEDIR)/main_2.cc
#GRIDTYPE = ALUGRID_SIMPLEX
#GRIDTYPE=PARALLELGRID_ALUGRID_CUBE
#GRIDTYPE=PARALLELGRID_ALUGRID_SIMPLEX
#GRIDTYPE=SPGRID_COUNT_FLOPS
GRIDTYPE=SPGRID_COUNT_FLOPS
#GRIDTYPE=SPGRID
#GRIDTYPE=CARTESIANGRID_ALUGRID_CUBE
# GRIDTYPE = ALUGRID_CUBE
GRIDTYPE=YASPGRID
#GRIDTYPE = ALUGRID_CUBE
#GRIDTYPE=YASPGRID
#GRIDTYPE = YASPGRID
GRIDDIM=2
POLORDER=2
......@@ -31,7 +31,6 @@ PROBLEM = 2 # check ./problemtype.hh
FLUX = 1 # check ./problemcreator.hh
DIFFFLUX=PRIMALDG
VTK=1
#DIFFFLUX=FLUXDG
# INFO SPACE OPERATOR:
......@@ -44,7 +43,7 @@ VTK=1
# 1. define LIMITER to limit the advection part of the solution (not checked)
AM_CPPFLAGS = $(USE_OMP) -I../../problems/nseq $(ALL_PKG_CPPFLAGS) -DGRIDDIM=$(GRIDDIM) \
-DPROBLEM=$(PROBLEM) -D$(GRIDTYPE) $(DUNEMPICPPFLAGS) \
-DFLUX=$(FLUX) -DUSE_VTKWRITER=$(VTK) -D$(DIFFFLUX) -DCOUNT_FLOPS # -DLEGENDREBASIS # -DWBPROBLEM -DWELLBALANCE
-DFLUX=$(FLUX) -D$(DIFFFLUX) -DLEGENDREBASIS -DCOUNT_FLOPS # -DLEGENDREBASIS # -DWBPROBLEM -DWELLBALANCE
AM_LDFLAGS = $(ALL_PKG_LDFLAGS) $(LAPACK_LDFLAGS) $(USE_OMP)
EXTRA_PROGRAMS = nseq nseq12
......@@ -65,14 +64,17 @@ EXTRA_DIST = parameter
CLEANFILES = manager.*.log
PROG=nseqonep
# width of SIMD unit (usually 2 or 4)
SIMDWIDTH=4
# codegeneration
generatecode:
$(MAKE) -i clean
$(MAKE) CODEGEN_OBJFILE= CXXFLAGS="-g -Wall -Wfatal-errors -DBASEFUNCTIONSET_CODEGEN_GENERATE" $(PROG)
$(MAKE) CXXFLAGS="-g -Wall -Wfatal-errors -DBASEFUNCTIONSET_CODEGEN_GENERATE -DCODEGEN_SIMD_WIDTH=$(SIMDWIDTH)" $(PROG)
$(MAKE) generate
generate:
./$(PROG) femhowto.eocSteps:1 femhowto.maximaltimesteps:1 fem.io.outputformat:none fem.ode.order:1 parameter
./$(PROG) fem.eoc.steps:1 fem.io.outputformat:none fem.ode.order:1 parameter
compilecode:
$(MAKE) clean
......
dune/fem-dg/test/navierstokes/ns_model.hh
+ 5
− 4
View file @ 62064bd5
......@@ -51,7 +51,7 @@ class NSModelTraits
typedef Intersection IntersectionType;
typedef typename GridPartType :: template Codim<0> :: EntityType EntityType;
typedef Thermodynamics< dimDomain > ThermodynamicsType;
typedef typename ProblemType :: ThermodynamicsType ThermodynamicsType;
typedef MinModLimiter< DomainFieldType > LimiterFunctionType;
};
......@@ -79,6 +79,7 @@ class NSModel : public DefaultModel < NSModelTraits< GridPartType, ProblemImp >
typedef typename Traits :: DomainType DomainType;
typedef typename Traits :: RangeType RangeType;
typedef typename Traits :: RangeFieldType RangeFieldType;
typedef typename Traits :: GradientRangeType GradientRangeType;
typedef typename Traits :: JacobianRangeType JacobianRangeType;
typedef typename Traits :: JacobianFluxRangeType JacobianFluxRangeType;
......@@ -300,8 +301,8 @@ class NSModel : public DefaultModel < NSModelTraits< GridPartType, ProblemImp >
abort();
DomainType xgl=it.intersectionGlobal().global(x);
const typename Traits :: DomainType normal = it.integrationOuterNormal(x);
double p;
double T;
RangeFieldType p;
RangeFieldType T;
pressAndTemp( uLeft, p, T );
gLeft = 0;
......@@ -385,7 +386,7 @@ class NSModel : public DefaultModel < NSModelTraits< GridPartType, ProblemImp >
nsFlux_.diffusion( u, jac, diff );
}
inline void pressAndTemp( const RangeType& u, double& p, double& T ) const
inline void pressAndTemp( const RangeType& u, RangeFieldType& p, RangeFieldType& T ) const
{
thermodynamics_.pressAndTempEnergyForm( u, p, T );
}
......
dune/fem-dg/test/navierstokes/nswaves.hh
+ 14
− 14
View file @ 62064bd5
......@@ -31,7 +31,7 @@ class NSWaves : public EvolutionProblemInterface<
//double,
GridType::dimension, GridType::dimension + 2 >,
true >,
public Thermodynamics< GridType::dimensionworld >
public Thermodynamics< GridType::dimensionworld, typename GridType::ctype >
{
NSWaves( const NSWaves& );
public:
......@@ -49,7 +49,7 @@ class NSWaves : public EvolutionProblemInterface<
typedef typename FunctionSpaceType :: DomainType DomainType;
typedef typename FunctionSpaceType :: RangeFieldType RangeFieldType;
typedef typename FunctionSpaceType :: RangeType RangeType;
typedef Thermodynamics< dimension > ThermodynamicsType;
typedef Thermodynamics< dimension, RangeFieldType > ThermodynamicsType;
NSWaves() : ThermodynamicsType(),
myName_( "NSWaves" ),
......@@ -66,7 +66,7 @@ class NSWaves : public EvolutionProblemInterface<
// initialize A and B
double init(const bool returnA ) const ;
RangeFieldType init(const bool returnA ) const ;
// print info
void printInitInfo() const;
......@@ -93,19 +93,19 @@ class NSWaves : public EvolutionProblemInterface<
}
//template< class RangeImp >
double pressure( const RangeType& u ) const
RangeFieldType pressure( const RangeType& u ) const
{
return thermodynamics().pressureEnergyForm( u );
}
double temperature( const RangeType& u ) const
RangeFieldType temperature( const RangeType& u ) const
{
return thermodynamics().temperatureEnergyForm( u, pressure( u ) );
}
// pressure and temperature
template< class RangeImp >
inline void pressAndTemp( const RangeImp& u, double& p, double& T ) const;
inline void pressAndTemp( const RangeImp& u, RangeFieldType& p, RangeFieldType& T ) const;
/* \brief finalize the simulation using the calculated numerical
......@@ -138,10 +138,10 @@ class NSWaves : public EvolutionProblemInterface<
void paraview_conv2prim() const {}
std::string description() const;
inline double mu( const RangeType& ) const { return mu_; }
inline double mu( const double T ) const { return mu_; }
inline double lambda( const double T ) const { return -2./3.*mu(T); }
inline double k( const double T ) const { return c_pd() *mu(T) * Pr_inv(); }
inline RangeFieldType mu( const RangeType& ) const { return mu_; }
inline RangeFieldType mu( const RangeFieldType T ) const { return mu_; }
inline RangeFieldType lambda( const RangeFieldType T ) const { return -2./3.*mu(T); }
inline RangeFieldType k( const RangeFieldType T ) const { return c_pd() *mu(T) * Pr_inv(); }
protected:
const std::string myName_;
......@@ -157,8 +157,8 @@ protected:
template <class GridType>
inline double NSWaves<GridType>
:: init(const bool returnA ) const
inline typename NSWaves<GridType> :: RangeFieldType
NSWaves<GridType>:: init(const bool returnA ) const
{
if( dimension == 1 )
{
......@@ -183,7 +183,7 @@ inline double NSWaves<GridType>
}
abort();
return 0;
return RangeFieldType(0);
}
......@@ -311,7 +311,7 @@ inline void NSWaves<GridType>
template <class GridType>
template< class RangeImp >
inline void NSWaves<GridType>
:: pressAndTemp( const RangeImp& u, double& p, double& T ) const
:: pressAndTemp( const RangeImp& u, RangeFieldType& p, RangeFieldType& T ) const
{
thermodynamics().pressAndTempEnergyForm( u, p, T );
}
......
dune/fem-dg/test/navierstokes/parameter
+ 3
− 2
View file @ 62064bd5
......@@ -11,10 +11,11 @@
fem.prefix: data
fem.io.datafileprefix: NS
fem.io.outputformat: vtk-cell
fem.io.savestep: 0.01 # SaveStep (write data every `saveStep' time period, <=0 deactivates)
fem.io.savecount: -1 # SaveCount (write data every saveCount time steps, <=0 deactivates)
femdg.stepper.printcount: 1000
femhowto.printCount: 1
femhowto.additionalvariables: false
......@@ -49,7 +50,7 @@ Pr: 0.72 # 0.72 #(Prandtl's number)
femhowto.startTime: 0. # s
femhowto.endTime: 0.1
femdg.stepper.endtime: 0.1
fem.eoc.steps: 2
fem.eoc.steps: 1
gammaGNS: 0.5
omegaGNS: 1.
......
dune/fem-dg/test/navierstokes/thermodynamics.hh
+ 68
− 68
View file @ 62064bd5
......@@ -21,7 +21,7 @@ using namespace Dune;
*
* \tparam dimDomain dimension of the domain
*/
template< int dimDomain >
template< int dimDomain, class Field = double >
class Thermodynamics
{
enum{ energyId = dimDomain+1 };
......@@ -72,14 +72,14 @@ class Thermodynamics
* \return pressure in energy form
*/
template< class RangeType >
inline double pressureEnergyForm( const RangeType& cons ) const
inline Field pressureEnergyForm( const RangeType& cons ) const
{
// cons = [rho, rho*v, rho*e]
assert( cons[0] > 1e-20 );
assert( cons[energyId] > 1e-20 );
// kinetic energy
double kin = 0.;
Field kin = 0.;
for( int i=1; i<=dimDomain; ++i )
kin += cons[ i ] * cons[ i ];
kin *= 0.5 / cons[ 0 ];
......@@ -101,8 +101,8 @@ class Thermodynamics
* \tparam RangeType Type of the range value
*/
template< class RangeType >
inline double temperatureEnergyForm( const RangeType& cons,
const double p ) const
inline Field temperatureEnergyForm( const RangeType& cons,
const Field p ) const
{
assert( cons[0] > 1e-20 );
return R_d_inv_ * p / cons[ 0 ];
......@@ -120,7 +120,7 @@ class Thermodynamics
* \return pressure in energy form
*/
template< class RangeType >
inline double temperatureEnergyForm( const RangeType& cons ) const
inline Field temperatureEnergyForm( const RangeType& cons ) const
{
return temperatureEnergyForm( cons, pressureEnergyForm( cons ) );
}
......@@ -135,15 +135,15 @@ class Thermodynamics
* \tparam RangeType Type of the range value
*/
template< class RangeType >
inline double densityThetaForm( const RangeType& prim ) const
inline Field densityThetaForm( const RangeType& prim ) const
{
const double p = prim[ energyId - 1 ];
const double theta = prim[ energyId ];
const Field p = prim[ energyId - 1 ];
const Field theta = prim[ energyId ];
assert( p > 1e-12 );
assert( theta > 1e-12 );
const double rho = std::pow( p/p0_ , gamma_inv_ ) * p0_ * R_d_inv_ / theta ;
const Field rho = std::pow( p/p0_ , gamma_inv_ ) * p0_ * R_d_inv_ / theta ;
assert( rho > 0.0 );
return rho;
......@@ -160,62 +160,62 @@ class Thermodynamics
void conservativeToPrimitiveEnergyForm( const RangeType& cons, RangeType& prim ) const;
public:
inline double g() const { return g_; }
inline double c_pd() const { return c_pd_; }
inline double c_pd_inv() const { return c_pd_inv_; }
inline double c_vd() const { return c_vd_; }
inline double c_vd_inv() const { return c_vd_inv_; }
inline double c_pl() const { return c_pl_; }
inline double c_pl_inv() const { return c_pl_inv_; }
inline double R_d() const { return R_d_; }
inline double R_d_inv() const { return R_d_inv_; }
inline double T0() const { return T0_; }
inline double T0_inv() const { return T0_inv_; }
inline double L0() const { return L0_; }
inline double p0Star() const { return p0Star_; }
inline double p0() const { return p0_; }
inline double p0_inv() const { return p0_inv_; }
inline double gamma() const { return gamma_; }
inline double gamma_inv() const { return gamma_inv_; }
inline double kappa() const { return kappa_; }
inline double kappa_inv() const { return kappa_inv_; }
inline Field g() const { return g_; }
inline Field c_pd() const { return c_pd_; }
inline Field c_pd_inv() const { return c_pd_inv_; }
inline Field c_vd() const { return c_vd_; }
inline Field c_vd_inv() const { return c_vd_inv_; }
inline Field c_pl() const { return c_pl_; }
inline Field c_pl_inv() const { return c_pl_inv_; }
inline Field R_d() const { return R_d_; }
inline Field R_d_inv() const { return R_d_inv_; }
inline Field T0() const { return T0_; }
inline Field T0_inv() const { return T0_inv_; }
inline Field L0() const { return L0_; }
inline Field p0Star() const { return p0Star_; }
inline Field p0() const { return p0_; }
inline Field p0_inv() const { return p0_inv_; }
inline Field gamma() const { return gamma_; }
inline Field gamma_inv() const { return gamma_inv_; }
inline Field kappa() const { return kappa_; }
inline Field kappa_inv() const { return kappa_inv_; }
//! the Reynolds number Re
inline double Re() const { return Re_; }
inline Field Re() const { return Re_; }
//! the inverse Reynolds number (1/Re)
inline double Re_inv() const { return Re_inv_; }
inline Field Re_inv() const { return Re_inv_; }
//! the Prandtl number Pr
inline double Pr() const { return Pr_; }
inline Field Pr() const { return Pr_; }
//! the inverse Prandtl number (1/Pr)
inline double Pr_inv() const { return Pr_inv_; }
inline Field Pr_inv() const { return Pr_inv_; }
private:
const double Re_;
const double Pr_;
const double g_;
const double p0_; // surface pressure
const double p0Star_;
const double T0_; // freezing temperature
const double c_pd_; // specific heat capacity of dry air w.r.t. pressure
const double c_vd_; // specific heat capacity of water vapour w.r.t. pressure
const double c_pl_; // specific heat capacity of liquid water w.r.t. pressure
const double L0_; // latent heat of evaporasion at 0 Celsius [J/kg]
const double relCloud_;
const double Re_inv_;
const double Pr_inv_;
const double c_pd_inv_;
const double c_vd_inv_;
const double c_pl_inv_;
const double R_d_; // gas constant for dry air
const double R_d_inv_;
const double p0_inv_;
const double T0_inv_;
const double kappa_;
const double kappa_inv_;
const double gamma_;
const double gammaM1_;
const double gamma_inv_;
const Field Re_;
const Field Pr_;
const Field g_;
const Field p0_; // surface pressure
const Field p0Star_;
const Field T0_; // freezing temperature
const Field c_pd_; // specific heat capacity of dry air w.r.t. pressure
const Field c_vd_; // specific heat capacity of water vapour w.r.t. pressure
const Field c_pl_; // specific heat capacity of liquid water w.r.t. pressure
const Field L0_; // latent heat of evaporasion at 0 Celsius [J/kg]
const Field relCloud_;
const Field Re_inv_;
const Field Pr_inv_;
const Field c_pd_inv_;
const Field c_vd_inv_;
const Field c_pl_inv_;
const Field R_d_; // gas constant for dry air
const Field R_d_inv_;
const Field p0_inv_;
const Field T0_inv_;
const Field kappa_;
const Field kappa_inv_;
const Field gamma_;
const Field gammaM1_;
const Field gamma_inv_;
};
......@@ -231,17 +231,17 @@ private:
* \tparam dimDomain dimension of the domain
* \tparam RangeType type of the range value
*/
template< int dimDomain >
template< int dimDomain, class Field >
template< class RangeType >
void Thermodynamics< dimDomain >
void Thermodynamics< dimDomain, Field >
:: conservativeToPrimitiveEnergyForm( const RangeType& cons, RangeType& prim ) const
{
std::cerr <<"conservativeToPrimitiveEnergyForm not implemented" <<std::endl;
abort();
//const double rho_inv = 1./ cons[0];
//const Field rho_inv = 1./ cons[0];
//double p, T;
//Field p, T;
//pressAndTempEnergyForm( cons, p, T );
//prim[energyId-1] = p;
......@@ -262,15 +262,15 @@ void Thermodynamics< dimDomain >
* \tparam dimDomain dimension of the domain
* \tparam RangeType type of the range value
*/
template< int dimDomain >
template< int dimDomain, class Field >
template< class RangeType >
void Thermodynamics< dimDomain >
void Thermodynamics< dimDomain, Field >
:: conservativeToPrimitiveThetaForm( const RangeType& cons, RangeType& prim ) const
{
assert( cons[0] > 0. );
assert( cons[energyId] > 0. );
double p, T;
Field p, T;
pressAndTempThetaForm( cons, p, T );
for( int i = 0; i < dimDomain; ++i )
......@@ -280,9 +280,9 @@ void Thermodynamics< dimDomain >
prim[energyId] = cons[energyId] / cons[0];
}
template< int dimDomain >
template< int dimDomain, class Field >
template< class RangeType >
void Thermodynamics< dimDomain >
void Thermodynamics< dimDomain, Field >
:: primitiveToConservativeThetaForm( const RangeType& prim, RangeType& cons ) const
{
// p,theta --> rho
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment