navier stokes test.

dune/fem-dg/test/navierstokes/Makefile.am

+# install these headers 
+nseqdir=$(includedir)/src/problems/nseq
+nseq_HEADERS= ns_model.hh problemcreator.hh nswaves.hh \
+ns_model_spec.hh constant.hh problemtype.hh nssmooth.hh thermodynamics.hh passtraits.hh
+
+#USE_OMP=-fopenmp
+#USE_OMP=-DUSE_PTHREADS=1
+
+LDADD=$(ALL_PKG_LDFLAGS) $(ALL_PKG_LIBS) $(LOCAL_LIBS) $(DUNEMPILDFLAGS) $(DUNEMPILIBS)
+
+BASEDIR=$(DUNE_FEM_DG_ROOT)/dune/fem-dg/main
+
+SOURCES = $(BASEDIR)/main.cc $(BASEDIR)/main_0.cc $(BASEDIR)/main_1.cc \
+						$(BASEDIR)/main_2.cc $(BASEDIR)/main_3.cc $(BASEDIR)/main_4.cc
+SOURCES12 = $(BASEDIR)/main.cc $(BASEDIR)/main_1.cc $(BASEDIR)/main_2.cc
+SOURCESONEP = $(BASEDIR)/main.cc $(BASEDIR)/main_pol.cc 
+
+#GRIDTYPE = ALUGRID_SIMPLEX
+#GRIDTYPE=PARALLELGRID_ALUGRID_CUBE
+#GRIDTYPE=PARALLELGRID_ALUGRID_SIMPLEX
+GRIDTYPE=SPGRID
+#GRIDTYPE=CARTESIANGRID_ALUGRID_CUBE
+#GRIDTYPE = ALUGRID_CUBE
+#GRIDTYPE=YASPGRID
+#GRIDTYPE = YASPGRID
+GRIDDIM=2
+POLORDER=2
+
+PROBLEM = 2 # check ./problemtype.hh
+FLUX = 1 # check ./problemcreator.hh
+
+DIFFFLUX=PRIMALDG
+#DIFFFLUX=FLUXDG
+
+# INFO SPACE OPERATOR:
+# 		1. define PRIMALDG for various space operators in primal formulation
+# 	  2. define DUALDG for the LDG space operator in dual formulation
+# INFO TRACK LIFTING:
+# 		1. define LOCALDEBUG to calculate \sum_e\int_Omega(r_e*l_e) and 
+# 			 \sum_e\int_Omega(r_e*l_e). They will be output to std::cout from the Stepper
+# INFO LIMITER
+# 		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) -D$(DIFFFLUX) -DLAGRANGEBASIS # -DLEGENDREBASIS # -DWBPROBLEM -DWELLBALANCE
+AM_LDFLAGS = $(ALL_PKG_LDFLAGS) $(LAPACK_LDFLAGS) $(USE_OMP) 
+
+
+check_PROGRAMS = constant constantonep constant12 \
+								 nseq nseqonep nseq12
+nodist_constant_SOURCES = $(SOURCES)
+nodist_constantonep_SOURCES = $(SOURCESONEP)
+nodist_constant12_SOURCES = $(SOURCES12) 
+nodist_nseq_SOURCES = $(SOURCES) 
+nodist_nseqonep_SOURCES = $(SOURCESONEP) 
+nodist_nseq12_SOURCES = $(SOURCES12) 
+
+constant_CPPFLAGS = $(AM_CPPFLAGS)
+constantonep_CPPFLAGS = $(constant_CPPFLAGS) -DONLY_ONE_P -DPOLORDER=$(POLORDER)
+constant12_CPPFLAGS = $(constant_CPPFLAGS) -DONLY_P1_P2
+nseq_CPPFLAGS = $(AM_CPPFLAGS)
+nseqonep_CPPFLAGS = $(nseq_CPPFLAGS) -DONLY_ONE_P -DPOLORDER=$(POLORDER)
+nseq12_CPPFLAGS = $(nseq_CPPFLAGS) -DONLY_P1_P2
+
+DISTCHECK_CONFIGURE_FLAGS = DOXYGEN="true"
+
+EXTRA_DIST = paramFiles/param*
+
+CLEANFILES = manager.*.log 
+
+PROG=nseqonep
+# codegeneration 
+generatecode:
+	$(MAKE) -i clean
+	$(MAKE) CODEGEN_OBJFILE=  CXXFLAGS="-g -Wall -Wfatal-errors -DBASEFUNCTIONSET_CODEGEN_GENERATE" $(PROG)
+	$(MAKE) generate
+
+generate:
+	./$(PROG) femhowto.eocSteps:1 femhowto.maximaltimesteps:1 fem.io.outputformat:none fem.ode.order:1 paramFiles/paramNSWaves 
+
+compilecode:
+	$(MAKE) clean 
+	$(MAKE) CXXFLAGS="$(CXXFLAGS) -DUSE_BASEFUNCTIONSET_CODEGEN" $(PROG)
+
+codegen:
+	$(MAKE) generatecode
+	$(MAKE) compilecode
+
+include $(top_srcdir)/am/global-rules

dune/fem-dg/test/navierstokes/autogeneratedcode.hh

+#if CODEDIM == 1
+#include "./autogeneratedcode/autogeneratedcode_1d.hh"
+#endif
+#if CODEDIM == 2
+#include "./autogeneratedcode/autogeneratedcode_2d.hh"
+#endif
+#if CODEDIM == 3
+#include "./autogeneratedcode/autogeneratedcode_3d.hh"
+#endif

dune/fem-dg/test/navierstokes/constant.hh

+#ifndef DUNE_CONSTANT_HH
+#define DUNE_CONSTANT_HH
+
+// dune-fem include
+#include <dune/fem/io/parameter.hh>
+#include <dune/fem/space/common/functionspace.hh>
+
+// local includes
+#include <dune/fem/probleminterfaces.hh>
+#include "thermodynamics.hh"
+
+
+/*****************************************************
+ *                                                   *
+ * 2d problem                                        *
+ * NSConstant solution problem                       *
+ *   with exact analytical solution for Euler and    *
+ *   NS eqns with source term equal 0                *
+ *                                                   *
+ ****************************************************/
+
+namespace Dune {
+
+template <class GridType>
+class NSConstant : public EvolutionProblemInterface<
+                  Dune::FunctionSpace< double, double, GridType::dimension, GridType::dimension + 2 >,
+                  true >,
+                public Thermodynamics< GridType::dimensionworld >
+{
+  NSConstant( const NSConstant& );
+ public:
+  typedef FunctionSpace<typename GridType::ctype, 
+                        double, GridType::dimensionworld,
+                        GridType::dimensionworld + 2 > FunctionSpaceType ;
+
+  enum{ dimension = GridType::dimensionworld };
+  typedef typename FunctionSpaceType :: DomainFieldType DomainFieldType;
+  typedef typename FunctionSpaceType :: DomainType DomainType;
+  typedef typename FunctionSpaceType :: RangeFieldType  RangeFieldType;
+  typedef typename FunctionSpaceType :: RangeType  RangeType;
+  typedef Thermodynamics< dimension >                     ThermodynamicsType;
+
+
+  NSConstant() 
+    : myName_( "NSConstant" )
+
+    , Re_( Parameter::getValue< double >( "reynold" ) )
+    , Re_inv_( 1. / Re_ )
+    , Pr_( Parameter::getValue< double >( "prandtl" ) )
+    , Pr_inv_( 1. / Pr_ )
+  {
+    init();
+  }
+
+
+  void init();
+
+
+  void printInitInfo() const;
+
+  // source implementations 
+  inline bool hasStiffSource() { return false; }
+  inline bool hasNonStiffSource() { return false; }
+  inline double stiffSource( const double t, const DomainType& x, RangeType& res ) const;
+  inline double nonStiffSource( const double t, const DomainType& x, RangeType& res ) const;
+
+  // this is U0
+  inline void evaluate( const DomainType& arg, RangeType& res ) const
+  {
+    evaluate( 0, arg, res );
+  }
+
+
+  inline void evaluate( const DomainType& arg, const double t, RangeType& res ) const
+  {
+    evaluate( t, arg, res );
+  }
+
+
+  inline void evaluate( const double t, const DomainType& arg, RangeType& res ) const
+  {
+    res[0] = 1.;
+    res[1] = 0.;
+    res[2] = 0.;
+#if GRIDDIM==3    
+    res[3] = 0.;
+    res[4] = 3.;
+#elif GRIDDIM==2
+    res[3] = 3.;
+#endif
+  }
+
+
+  template <class RangeImp>
+  inline void pressAndTemp( const RangeImp& u, double& p, double& T ) const;
+
+
+  const ThermodynamicsType& thermodynamics() const { return *this; }
+  inline double Re() const { return Re_; }
+  inline double Re_inv() const { return Re_inv_; }
+  inline double Pr() const { return Pr_; }
+  inline double Pr_inv() const { return Pr_inv_; }
+  inline double c_p() const { return c_p_; }
+  inline double c_v_inv() const { return c_v_inv_; }
+  inline double R_d() const { return R_d_; }
+  inline double R_v() const { return R_v_; }
+  inline double endTime() const { return endTime_; }
+  inline double gamma() const { return gamma_; }
+  inline double g() const { return g_; }
+  inline std::string myName() const { return myName_; }
+  inline double lambda() const { return -2./3.*mu_; }
+  inline double lambda( const double T ) const { return -2./3.*mu_; }
+  inline double mu() const { return mu_; }
+  inline double mu( const double T ) const { return mu_; }
+  inline double k() const { return c_p_*mu()*Pr_inv_; }
+  inline double k( const double T ) const { return c_p_*mu(T)*Pr_inv_; }
+  void printmyInfo( const std::string& filename )const {}
+  std::string description() const { return ""; }
+  void paraview_conv2prim() const {}
+
+private:
+  const ThermodynamicsType& thermodynamics_;
+  const std::string myName_;
+  const double Re_;
+  const double Re_inv_;
+  const double Pr_;
+  const double Pr_inv_;
+  double gamma_;
+  double g_;
+  double endTime_;
+  double c_v_;
+  double c_v_inv_;
+  double c_p_;
+  double c_p_inv_;
+  double R_d_;
+  double R_d_inv_;
+  double R_v_;
+  double R_v_inv_;
+  double lambda_;
+  double mu_;
+};
+
+
+
+
+template <class GridType> 
+inline void NSConstant<GridType>
+:: init()
+{
+  g_ = Dune::Parameter::getValue< double >( "g" );
+  endTime_ = Dune::Parameter::getValue< double >( "femhowto.endTime" );
+  mu_ = Dune::Parameter :: getValue< double >( "mu" ); 
+  c_v_ = Dune::Parameter :: getValue< double >( "c_v" );
+  c_p_ = Dune::Parameter :: getValue< double >( "c_p" );
+
+  // additional coefficients
+  assert( g_ == 0. );
+  c_v_inv_ = 1. / c_v_;
+  c_p_inv_ = 1. / c_p_;
+  R_d_ = c_p_ - c_v_;
+  R_d_inv_ = 1. / R_d_;
+  gamma_ = c_p_ * c_v_inv_;
+}
+
+
+
+template <class GridType> 
+inline void NSConstant<GridType>
+:: printInitInfo() const
+{
+  std::cout <<R_d_ <<" = R_d.\n";
+  std::cout <<R_v_ <<" = R_v.\n";
+  std::cout <<c_v_inv_ <<" = c_v_inv.\n";
+  std::cout <<mu_ <<" = mu\n";
+  std::cout <<c_p_ <<" = c_p_.\n";
+  std::cout <<gamma_ <<" = gamma.\n";
+  std::cout <<"DONE.\n";
+}
+
+
+template <class GridType>
+inline double NSConstant<GridType>
+:: stiffSource( const double t, const DomainType& x, RangeType& res ) const
+{
+  res = 0.;
+  return 0.;
+}
+
+
+template <class GridType>
+inline double NSConstant<GridType>
+:: nonStiffSource( const double t, const DomainType& x, RangeType& res ) const
+{
+  res = 0.;
+  return 0.;
+}
+
+
+template <class GridType> 
+template <class RangeImp>
+inline void NSConstant<GridType>
+:: pressAndTemp( const RangeImp& u, double& p, double& T ) const
+{
+  const double rho_inv = 1./u[0];
+#if GRIDDIM==2
+  const double rhoeps = u[3] - 0.5*rho_inv*(u[1]*u[1] + u[2]*u[2]);
+#elif GRIDDIM==3
+  const double rhoeps = u[4] - 0.5*rho_inv*(u[1]*u[1] + u[2]*u[2]+u[3]*u[3]);
+#else
+#error "Grid dimension must be 2 or 3"
+#endif
+  assert( rhoeps > 1e-15 );
+  p = (gamma_ - 1.)*rhoeps;
+  T = p * rho_inv * R_d_inv_;
+}
+
+}
+#endif

dune/fem-dg/test/navierstokes/ns_model.hh

+#ifndef NS_MODEL_HH
+#define NS_MODEL_HH
+
+// DUNE includes
+#include <dune/common/version.hh>
+#include <dune/fem/misc/fmatrixconverter.hh>
+
+// local includes
+#include <dune/fem-dg/models/defaultmodel.hh>
+#include <dune/fem-dg/operator/limiter/limitpass.hh>
+#include "ns_model_spec.hh"
+
+
+namespace Dune {
+
+//////////////////////////////////////////////////////
+//
+//                 NAVIER-STOKES EQNS 
+//
+//////////////////////////////////////////////////////
+template< class GridPart > 
+class NSModelTraits 
+{
+  public:
+  typedef GridPart GridPartType;
+  typedef typename GridPart :: GridType                     GridType;
+  enum { dimDomain = GridType :: dimensionworld };
+  enum { dimRange = dimDomain + 2 };
+  enum { dimGradRange = dimRange * dimDomain };
+  enum { dimGradient = dimDomain + 1 };
+
+  typedef Fem :: FunctionSpace< double, double, dimDomain, dimRange > FunctionSpaceType ;
+  
+  typedef typename FunctionSpaceType :: DomainType         DomainType ;
+  typedef typename FunctionSpaceType :: DomainFieldType    DomainFieldType ;
+  typedef FieldVector< DomainFieldType, dimDomain - 1 >    FaceDomainType;
+  typedef typename FunctionSpaceType :: RangeType          RangeType ;
+
+  typedef FieldVector< double, dimGradRange >               GradientType;
+  typedef typename FunctionSpaceType :: JacobianRangeType  JacobianRangeType;
+  typedef JacobianRangeType                                FluxRangeType;
+
+  typedef FieldVector< double, dimGradRange >               GradientRangeType;
+  typedef FieldMatrix< double, dimGradRange, dimDomain >    JacobianFluxRangeType;
+
+  typedef typename GridPart :: IntersectionIteratorType     IntersectionIterator;
+  typedef typename IntersectionIterator::Intersection       Intersection;
+  typedef Intersection       IntersectionType;
+  typedef typename GridPartType :: template Codim<0> :: EntityType  EntityType;
+  typedef typename EntityType :: EntityPointer              EntityPointerType;
+
+  typedef Thermodynamics< dimDomain >                       ThermodynamicsType;
+
+  typedef MinModLimiter< DomainFieldType > LimiterFunctionType;
+};
+
+
+
+template< class GridPartType , class ProblemImp >
+class NSModel : public DefaultModel < NSModelTraits< GridPartType > >
+{
+  public:
+  typedef ProblemImp                                        ProblemType;
+  typedef typename GridPartType::GridType                   GridType;
+  typedef NSModelTraits< GridPartType >                     Traits;
+
+  typedef NSFlux< Traits >  FluxType ;
+
+  enum { dimDomain = Traits :: dimDomain };
+  enum { dimRange  = Traits :: dimRange  };
+  enum { dimGradRange = Traits::dimGradRange } ;
+
+  typedef typename Traits :: EntityType                     EntityType;
+  typedef typename Traits :: EntityPointerType              EntityPointerType;
+  typedef typename Traits :: IntersectionIterator           IntersectionIterator;
+  typedef typename Traits :: Intersection                   IntersectionType;
+  typedef typename Traits :: FaceDomainType                 FaceDomainType;
+
+  typedef typename Traits :: DomainType                     DomainType;
+  typedef typename Traits :: RangeType                      RangeType;
+  typedef typename Traits :: GradientRangeType              GradientRangeType;
+  typedef typename Traits :: JacobianRangeType              JacobianRangeType;
+  typedef typename Traits :: JacobianFluxRangeType          JacobianFluxRangeType;
+  typedef typename Traits :: ThermodynamicsType             ThermodynamicsType;
+
+ public:
+  NSModel( const ProblemType& problem ) 
+    : thermodynamics_( problem.thermodynamics() )
+    , problem_( problem )
+    , nsFlux_( problem )
+    , alpha_( std::pow( problem.gamma(), 1.5 ) * (problem.Re_inv() * problem.Pr_inv()) )
+  {
+  }
+
+  double gamma() const { return problem_.gamma() ; }
+
+  inline bool hasStiffSource() const { return problem_.hasStiffSource(); }
+  inline bool hasNonStiffSource() const { return problem_.hasNonStiffSource(); }
+  inline bool hasFlux() const { return true ; }
+
+  inline double stiffSource( const EntityType& en,
+                        const double time,
+                        const DomainType& x,
+                        const RangeType& u,
+                        const GradientRangeType& du,
+                        RangeType & s) const
+  {
+    return stiffSource( en, time, x, u, s );
+  }
+
+
+  inline double stiffSource( const EntityType& en,
+                        const double time,
+                        const DomainType& x,
+                        const RangeType& u,
+                        const JacobianRangeType& jac,
+                        RangeType & s) const
+  {
+    return stiffSource( en, time, x, u, s ); 
+  }
+
+
+  inline double stiffSource( const EntityType& en
+                      , const double time
+                      , const DomainType& x
+                      , const RangeType& u
+                      , RangeType& s ) const 
+  {
+    // some special RHS for testcases/NSWaves
+    const DomainType& xgl = en.geometry().global(x);
+    return problem_.stiffSource( time, xgl, u, s );
+  }
+
+
+  inline double nonStiffSource( const EntityType& en,
+                                const double time,
+                                const DomainType& x,
+                                const RangeType& u,
+                                const GradientRangeType& du,
+                                RangeType & s) const
+  {
+    Fem::FieldMatrixConverter< GradientRangeType, JacobianRangeType > jac( du );
+    return nonStiffSource( en, time, x, u, jac, s );
+  }
+
+
+
+  template< class JacobianRangeTypeImp >
+  inline double nonStiffSource( const EntityType& en,
+                                const double time,
+                                const DomainType& x,
+                                const RangeType& u,
+                                const JacobianRangeTypeImp& jac,
+                                RangeType& s) const
+  {
+    const DomainType& xgl = en.geometry().global(x);
+    return problem_.nonStiffSource( time, xgl, u, s );
+  }
+
+  inline void advection( const EntityType& en,
+                         const double time,
+                         const DomainType& x,
+                         const RangeType& u,
+                         JacobianRangeType& f ) const 
+  {
+    nsFlux_.analyticalFlux( u, f );
+  }
+
+  /////////////////////////////////////////////////////////////////
+  // Limiter section 
+  ////////////////////////////////////////////////////////////////
+  inline void velocity( const EntityType& en,
+                        const double time,
+                        const DomainType& x,
+                        const RangeType& u,
+                        DomainType& velocity) const
+  {
+    for(int i=0; i<dimDomain; ++i)
+    {
+      // we store \rho u but do not need to divide by \rho here since only
+      // sign is needed.
+      velocity[i] = u[i+1];
+    }
+  }
+
+  // adjust average value if necessary 
+  // (e.g. transform from conservative to primitive variables )
+  void adjustAverageValue( const EntityType& entity,
+                           const DomainType& xLocal,
+                           RangeType& u ) const
+  {
+    // nothing to be done here for this test case 
+  }
+
+  // calculate jump between left and right value 
+  inline void jump(const IntersectionType& it,
+                   const double time,
+                   const FaceDomainType& x,
+                   const RangeType& uLeft,
+                   const RangeType& uRight,
+                   RangeType& jump) const
+  {
+    const double pl = pressure( uLeft );
+    const double pr = pressure( uRight );
+
+    // take pressure as shock detection values 
+    jump  = (pl-pr)/(0.5*(pl+pr));
+  }
+
+  // calculate jump between left and right value 
+  inline void adaptationIndicator(
+                   const IntersectionType& it,
+                   const double time,
+                   const FaceDomainType& x,
+                   const RangeType& uLeft,
+                   const RangeType& uRight,
+                   RangeType& indicator) const
+  {
+    jump( it, time, x, uLeft, uRight, indicator );
+  }
+
+
+  //! \brief returns true if physical check does something useful */
+  bool hasPhysical () const { return true; }
+
+  template <class LocalDomainType>
+  bool physical( const EntityType& en,
+                 const LocalDomainType& x,
+                 const RangeType& u) const
+  {
+    if (u[0]<1e-8)
+      return false;
+    else
+    {
+      const double p = pressure( u );
+      return p > 1e-8;
+    }
+
+
+    return u[ 0 ] > 1e-8 ;
+  }
+
+
+  inline double diffusionTimeStep( const IntersectionType& it,
+                                   const double enVolume,
+                                   const double circumEstimate,
+                                   const double time,
+                                   const FaceDomainType& x,
+                                   const RangeType& u ) const
+  {
+    // look at Ch. Merkle Diplom thesis, pg. 38
+    // get value of mu at temperature T
+
+    // get mu 
+    const double mu = problem_.mu( u );
+
+    // ksi = 0.25 
+    return mu * circumEstimate * alpha_ / (0.25 * u[0] * enVolume);
+  }
+
+
+  //! return analyticalFlux for 1st pass 
+  inline void jacobian( const EntityType& en,
+                        const double time,
+                        const DomainType& x,
+                        const RangeType& u, 
+                        JacobianFluxRangeType& a ) const 
+  {
+    nsFlux_.jacobian( u, a );
+  }
+  
+
+  inline bool hasBoundaryValue( const IntersectionType& it,
+                                const double time,
+                                const FaceDomainType& x ) const 
+  { 
+    return true;
+  }
+ 
+
+  inline double boundaryFlux( const IntersectionType& it,
+                              const double time,
+                              const FaceDomainType& x,
+                              const RangeType& uLeft,
+                              RangeType& gLeft ) const
+  {
+    abort();
+    gLeft = 0;
+    return 0.0;
+  }
+
+  inline double boundaryFlux( const IntersectionType& it,
+                              const double time,
+                              const FaceDomainType& x,
+                              const RangeType& uLeft,
+                              const GradientRangeType& duLeft,
+                              RangeType& gLeft ) const
+  {
+    abort();
+    DomainType xgl=it.intersectionGlobal().global(x);
+    const typename Traits :: DomainType normal = it.integrationOuterNormal(x); 
+    double p;
+    double T;
+    pressAndTemp( uLeft, p, T );
+    gLeft = 0;
+
+    // bnd. cond. from euler part
+    for (int i=0;i<dimDomain; ++i)
+      gLeft[i+1] = normal[i]*p;
+
+    return 0.0;
+  } 
+
+  inline double diffusionBoundaryFlux( const IntersectionType& it,
+                                       const double time,
+                                       const FaceDomainType& x,
+                                       const RangeType& uLeft,
+                                       const GradientRangeType& gradLeft,
+                                       RangeType& gLeft ) const  
+  {
+    Fem::FieldMatrixConverter< GradientRangeType, JacobianRangeType> jacLeft( gradLeft );
+    return diffusionBoundaryFlux( it, time, x, uLeft, jacLeft, gLeft );
+  }
+
+  /** \brief boundary flux for the diffusion part
+   */
+  template <class JacobianRangeImp>
+  inline double diffusionBoundaryFlux( const IntersectionType& it,
+                                       const double time,
+                                       const FaceDomainType& x,
+                                       const RangeType& uLeft,
+                                       const JacobianRangeImp& jacLeft,
+                                       RangeType& gLeft ) const  
+  {
+    std::cerr <<"diffusionBoundaryFlux shouldn't be used for this testcase" <<std::endl;
+    abort();
+  }
+
+
+  inline void boundaryValue( const IntersectionType& it,
+                             const double time,
+                             const FaceDomainType& x,
+                             const RangeType& uLeft,
+                             RangeType& uRight ) const
+  {
+    const DomainType xgl = it.geometry().global( x );
+    problem_.evaluate( time , xgl , uRight );
+  }
+  
+  // here x is in global coordinates
+  inline void maxSpeed( const EntityType& entity,
+                        const double time,
+                        const DomainType& x,
+                        const DomainType& normal,
+                        const RangeType& u,
+                        double& advspeed,
+                        double& totalspeed ) const 
+  {
+    advspeed = nsFlux_.maxSpeed( normal , u );
+    totalspeed=advspeed;
+  }
+
+
+  void diffusion( const EntityType& en,
+                  const double time,
+                  const DomainType& x,
+                  const RangeType& u,
+                  const GradientRangeType& v,
+                  JacobianRangeType& diff ) const
+  {
+    Fem::FieldMatrixConverter< GradientRangeType, JacobianRangeType> jac( v );
+    diffusion( en, time, x, u, jac, diff );
+  }
+
+
+  template <class JacobianRangeImp>
+  void diffusion( const EntityType& en,
+                  const double time,
+                  const DomainType& x,
+                  const RangeType& u,
+                  const JacobianRangeImp& jac,
+                  JacobianRangeType& diff ) const
+  {
+    nsFlux_.diffusion( u, jac, diff );
+  }
+
+  inline void pressAndTemp( const RangeType& u, double& p, double& T ) const
+  {
+    thermodynamics_.pressAndTempEnergyForm( u, p, T );
+  }
+
+  inline void conservativeToPrimitive( const double time,
+                                       const DomainType& xgl,
+                                       const RangeType& cons, 
+                                       RangeType& result,
+                                       bool ) const
+  {
+    problem_.evaluate( time, xgl, result );
+    //thermodynamics_.conservativeToPrimitiveEnergyForm( cons, result );
+  }
+
+
+  inline const ProblemType& problem() const { return problem_; }
+
+ protected:
+  const ThermodynamicsType& thermodynamics_;
+  const ProblemType& problem_;
+  const FluxType  nsFlux_;
+  const double alpha_;
+};
+
+} // end namespace Dune
+
+#endif // file definition

dune/fem-dg/test/navierstokes/ns_model_spec.hh

+#ifndef NS_MODEL_SPEC_HH
+#define NS_MODEL_SPEC_HH
+
+#include <dune/common/fvector.hh>
+
+#include <dune/fem-dg/operator/fluxes/eulerfluxes.hh>
+#include <dune/fem/misc/fmatrixconverter.hh>
+
+#include "problemtype.hh"
+
+namespace Dune {
+
+template< class Traits >
+class NSFlux
+{
+  enum { dimDomain = Traits :: dimDomain }; 
+
+  enum { e = dimDomain + 1 };
+  enum { e_ = e };
+  enum { dimRange = Traits :: dimRange };
+  enum { dimGradRange = dimRange * dimDomain };
+  
+public:
+  typedef typename Traits::DomainType  DomainType;
+  typedef typename Traits::RangeType   RangeType;
+  typedef Dune::FieldVector< double, dimGradRange > GradientRangeType;
+  typedef typename Traits::JacobianRangeType  JacobianRangeType;
+
+  typedef Dune::FieldMatrix< double, dimGradRange, dimDomain > JacobianFluxRangeType;
+  typedef Dune::Fem::FieldMatrixConverter< GradientRangeType, JacobianRangeType > 
+    ConvertedJacobianRangeType;
+  typedef Dune::FieldMatrix< double, dimDomain, dimDomain > VelocityGradientType;
+
+  NSFlux( const NSProblemType& problem )
+    : eulerFlux_() 
+    , problem_( problem )
+    , gamma_( problem.gamma() )
+    , R_d_inv_( problem.thermodynamics().R_d_inv() )
+    , Re_inv_( problem.Re_inv() )
+    , c_v_inv_( problem.thermodynamics().c_vd_inv() )
+    , c_v_( problem.thermodynamics().c_vd() )
+  {}
+
+  inline void analyticalFlux( const RangeType& u, JacobianRangeType& f ) const
+  {
+    eulerFlux_.analyticalFlux( gamma_ , u , f );
+  }
+
+  inline void jacobian( const RangeType& u, JacobianFluxRangeType& du ) const
+  {
+    du = 0;
+
+    assert( int(du.rows) == int(dimRange * dimDomain) );
+    assert( int(du.cols) == int(dimDomain) );
+
+    for (int r=0;r<dimRange;r++)
+      for (int d=0;d<dimDomain;d++)
+        du[dimDomain*r+d][d] = u[r];
+  }
+
+  inline double maxSpeed( const DomainType& n, const RangeType& u ) const
+  {
+    return eulerFlux_.maxSpeed( gamma_, n, u );
+  }
+
+  inline void diffusion( const RangeType& u,
+                         const GradientRangeType& du,
+                         JacobianRangeType& f ) const;
+
+  template <class JacobianRangeImp>
+  inline void diffusion( const RangeType& u,
+                         const JacobianRangeImp& du,
+                         JacobianRangeType& f ) const;
+
+  inline double mu( const double T ) const { return problem_.mu(T); }
+  inline double lambda( const double T ) const { return problem_.lambda(T); }
+  inline double k( const double T ) const { return problem_.k(T); }
+
+protected:
+  const EulerAnalyticalFlux<dimDomain> eulerFlux_;
+  const NSProblemType& problem_;
+  const double gamma_;
+  const double R_d_inv_;
+  const double Re_inv_;
+  const double c_v_inv_;
+  const double c_v_;
+};
+
+
+/////////////////////////////////////////////
+// Implementation of flux functions 
+/////////////////////////////////////////////
+template< class Traits >
+template< class JacobianRangeImp >
+void NSFlux< Traits > 
+:: diffusion( const RangeType& u,
+              const JacobianRangeImp& du,
+              JacobianRangeType& diff ) const 
+{
+  //std::cout << du << " du " << std::endl;
+  assert( u[0] > 1e-10 );
+  const double rho_inv = 1. / u[0];
+
+  // get velocity field
+  DomainType v;
+  for( int i = 0; i < dimDomain; ++i )
+    v[i] = u[ i+1 ] * rho_inv;
+
+  // | v |^2 
+  const double vTwoNorm = v.two_norm2(); 
+
+  // get all partial derivatives of all velocities
+  VelocityGradientType dVel;
+  for( int j = 0; j < dimDomain; ++j ) // v components 
+  {
+    for( int i = 0; i < dimDomain; ++i ) // space derivatives 
+    {
+      // substract d_x rho * v from the derivative of the conservative variables 
+      dVel[ j ][ i ] = rho_inv * (du[ j+1 ][ i ] - v[ j ]*du[ 0 ][ i ]);
+    }
+  } 
+
+  DomainType vGradVel;
+  for( int i = 0; i < dimDomain; ++i )
+  {
+    vGradVel[ i ] = 0.;
+    for( int j = 0; j < dimDomain; ++j )
+      vGradVel[ i ] += v[ j ] * dVel[ j ][ i ];
+  }
+
+  // get the absolute temperature
+  const double T = problem_.temperature( u );
+
+  // get all derivatives of the absolute temperature
+  DomainType dTemp;
+  for( int i = 0; i < dimDomain; ++i )
+  {
+    dTemp[ i ]  = du[ e_ ][ i ] - T * c_v_ * du[ 0 ][ i ] - 0.5*du[ 0 ][ i ]*vTwoNorm - u[ 0 ] * vGradVel[ i ];
+    dTemp[ i ] *= rho_inv * c_v_inv_;
+  }
+
+  // get temperature dependend viscosity coefficients
+  const double muLoc = mu( T );
+  const double kLoc = k( T );
+
+  // get divergence of the velocity field
+  double divVel = 0;
+  for( int i = 0; i < dimDomain; ++i )
+    divVel += dVel[ i ][ i ];
+
+  // apply lambda to divergence 
+  divVel *= lambda( T );
+
+  for( int i = 0; i < dimDomain; ++i )
+  {
+    // assemble the diffusion part for the density equation
+    // this equation is purely hyperbolic, so no diffusion 
+    diff[ 0 ][ i ] = 0;
+
+    // assemble the diffusion part for the momentum equation
+    // the stesstensor tau 
+    for( int j = 0; j < dimDomain; ++j )
+    {
+      diff[ j+1 ][ i ] = muLoc*(dVel[ j ][ i ] + dVel[ i ][ j ]);
+    }
+
+    // add   lambda * div v * 1I 
+    diff[ i+1 ][ i ] += divVel;
+  }
+
+  for( int i = 0; i < dimDomain; ++i )
+  {
+    // assemble the diffusion part for the energy (internal+kinetic) equation
+    for( int j = 0; j < dimDomain; ++j )
+      diff[e_][ i ] = v[ j ] * diff[ i+1 ][ j ];
+
+    //diff[e_][i] += problem_.Pr_inv() * kLoc * dTemp[i];
+    diff[e_][i] += kLoc * dTemp[ i ];
+  }
+
+  // scale diffusion part with Reynold's number
+  diff *= Re_inv_;
+}
+
+} // end namespace Dune
+#endif

dune/fem-dg/test/navierstokes/nssmooth.hh

+#ifndef NS_SMOOTH_SOLUTION_HH
+#define NS_SMOOTH_SOLUTION_HH
+
+#include <dune/common/version.hh>
+
+// dune-fem includes
+#include <dune/fem/misc/linesegmentsampler.hh>
+#include <dune/fem/io/parameter.hh>
+#include <dune/fem/space/common/functionspace.hh>
+
+// local includes
+#include "thermodynamics.hh"
+#include <dune/fem/probleminterfaces.hh>
+
+
+namespace Dune {
+
+template <class GridType>
+class NSSmoothSolution : public EvolutionProblemInterface<
+                            Dune::FunctionSpace< double, double, GridType::dimension, GridType::dimension + 2 >,
+                            true >,
+                         public Thermodynamics< GridType::dimensionworld >
+{
+  NSSmoothSolution( const NSSmoothSolution& );
+  public:
+  typedef FunctionSpace<typename GridType::ctype,
+                        double, GridType::dimensionworld,
+                        GridType::dimensionworld + 2 > FunctionSpaceType ;
+
+  enum{ dimension = GridType::dimensionworld };
+  enum { energyId = dimension + 1 };
+  typedef typename FunctionSpaceType :: DomainFieldType   DomainFieldType;
+  typedef typename FunctionSpaceType :: DomainType        DomainType;
+  typedef typename FunctionSpaceType :: RangeFieldType    RangeFieldType;
+  typedef typename FunctionSpaceType :: RangeType         RangeType;
+  typedef Thermodynamics< dimension >                     ThermodynamicsType;
+
+  NSSmoothSolution() : ThermodynamicsType(),
+    myName_( "NSSmoothSolution" ),
+    omegaGNS_( Parameter::getValue< double >( "omegaGNS" ) ),
+    kGNS_( Parameter::getValue< double >( "kGNS" ) ),
+    gammaGNS_( Parameter::getValue< double >( "gammaGNS" ) ),
+    endTime_ ( Parameter::getValue<double>( "femhowto.endTime" )), 
+    mu_( Parameter :: getValue< double >( "mu" )), 
+    //k_ ( c_pd() * Pr_inv() * mu_),
+    k_ ( 2.0 ),
+    alpha_( 1.0 ),
+    omega_( dimension ),
+    A_( init( true ) ),
+    B_( init( false ) ),
+    e0_( dimension ) // e0 > d/2
+  {
+  }
+
+
+  // initialize A and B 
+  double init(const bool returnA ) const ;
+
+  // print info 
+  void printInitInfo() const {}
+
+  // source implementations 
+  inline bool hasStiffSource() { return false; }
+  inline bool hasNonStiffSource() { return true; }
+  inline double stiffSource( const double t, const DomainType& x, const RangeType& u, RangeType& res ) const
+  {
+    res = 0;
+    return 0;
+  }
+
+  //! beta = k * sum_i x_i  - omega * t 
+  const double beta(const double t, const DomainType& x ) const
+  {
+    double sumX = 0;
+    for( int i=0; i< dimension; ++i ) sumX += x[i];
+    return k_ * sumX - omega_ * t ;
+  }
+
+  inline double nonStiffSource( const double t, const DomainType& x, const RangeType& u, RangeType& res ) const
+  {
+    const double betaa = beta( t, x );
+    const double cosBeta = std::cos( betaa );
+    res = cosBeta * alpha_ * ( dimension * k_ - omega_);
+
+    // velocity components 
+    const double add = (gamma() - 1.0) * ( e0_ - double(dimension) * 0.5 )
+                       * ( cosBeta * alpha_ * k_ );
+
+    res[ energyId ] *= e0_ ;
+
+    for( int i=1; i< RangeType :: dimension ; ++i ) res[ i ] += add ;
+    
+    // TODO: calculate time step restriction due to source term
+    return 0; 
+  }
+
+  // this is the initial data
+  inline void evaluate( const DomainType& arg , RangeType& res ) const 
+  {
+    evaluate( 0., arg, res );
+  }
+
+  // evaluate function 
+  inline void evaluate( const double t, const DomainType& x, RangeType& res ) const
+  {
+    // sinus waves 
+    res = std::sin( beta( t, x ) ) * alpha_ + 2.0 ;
+
+    // add factor for energy component 
+    res[ energyId ] *= e0_ ;
+
+    return ;
+  }
+
+  // cloned method 
+  inline void evaluate( const DomainType& x, const double t, RangeType& res ) const
+  {
+    evaluate( t, x, res );
+  }
+
+  //template< class RangeImp >
+  double pressure( const RangeType& u ) const 
+  {
+    return thermodynamics().pressureEnergyForm( u );
+  }
+
+  double 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
+  {
+    thermodynamics().pressAndTempEnergyForm( u, p, T );
+  }
+
+  /*  \brief finalize the simulation using the calculated numerical
+   *  solution u for this problem
+   *  
+   *  \param[in] variablesToOutput Numerical solution in the suitably chosen variables
+   *  \param[in] eocloop Specific EOC loop
+   */
+  template< class DiscreteFunctionType >
+  void finalizeSimulation( DiscreteFunctionType& variablesToOutput,
+                           const int eocloop) const
+  {}
+
+  const ThermodynamicsType& thermodynamics() const { return *this; }
+  using ThermodynamicsType :: Re ;
+  using ThermodynamicsType :: Re_inv;
+  using ThermodynamicsType :: Pr;
+  using ThermodynamicsType :: Pr_inv;
+  using ThermodynamicsType :: c_pd; 
+  using ThermodynamicsType :: c_pd_inv;
+  using ThermodynamicsType :: c_vd; 
+  using ThermodynamicsType :: c_vd_inv;
+  using ThermodynamicsType :: gamma;
+  using ThermodynamicsType :: g;
+  using ThermodynamicsType :: R_d_inv;
+  void printmyInfo( std::string filename ) const {}
+  inline double endtime() const { return endTime_; }
+  inline std::string myName() const { return myName_; }
+  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(); } 
+  
+protected:
+  const std::string myName_;
+  const double omegaGNS_;
+  const double kGNS_;
+  const double gammaGNS_;
+  const double endTime_;
+  const double mu_;
+  const double k_;
+  const double alpha_;
+  const double omega_ ;
+  const double A_;
+  const double B_;
+  const double e0_;
+};
+
+
+template <class GridType>
+inline double NSSmoothSolution<GridType>
+:: init(const bool returnA ) const 
+{
+  if( dimension == 1 )
+  {
+    if( returnA ) // A 
+      return (-omegaGNS_ + kGNS_*(3.5*gamma()-2.5));
+    else // B 
+      return ( kGNS_*(0.5+3.5*gamma()) - 4.*omegaGNS_ );
+  }
+  if( dimension == 2 ) 
+  {
+    if( returnA ) // A 
+      return (-omegaGNS_ + kGNS_*(3.*gamma() - 1.));
+    else // B 
+      return ( kGNS_*(2.+6.*gamma()) - 4.*omegaGNS_ );
+  }
+  else if( dimension == 3 ) 
+  {
+    if( returnA ) // A 
+      return (-omegaGNS_ + 0.5*kGNS_*(5.*gamma() + 1.));
+    else // B 
+     return (kGNS_*(4.5+7.5*gamma()) - 4.*omegaGNS_);
+  }
+
+  abort();
+  return 0;
+}
+
+
+template <class GridType>
+inline std::string NSSmoothSolution<GridType>
+:: description() const
+{
+  std::ostringstream stream;
+
+  stream <<"{\\bf Problem:}" <<myName_
+         <<", {\\bf $\\mu$:} " <<mu_
+         <<", {\\bf End time:} " <<endTime_
+         <<", {$\\gamma_{GNS}$:} " <<gammaGNS_ 
+         <<"\n"
+         <<", {$\\omega_{GNS}$:} " <<omegaGNS_
+         <<", {$k_{GNS}$:} " <<kGNS_;
+
+  std::string returnString = stream.str();
+
+  return returnString;
+}
+
+} // end namespace Dune
+#endif

dune/fem-dg/test/navierstokes/nswaves.hh

+#ifndef NSWAVES_HH
+#define NSWAVES_HH
+
+#include <dune/common/version.hh>
+
+// dune-fem includes
+#include <dune/fem/misc/linesegmentsampler.hh>
+#include <dune/fem/io/parameter.hh>
+#include <dune/fem/space/common/functionspace.hh>
+
+// local includes
+#include "thermodynamics.hh"
+#include <dune/fem/probleminterfaces.hh>
+
+
+/***********************************************************
+ *                                                   
+ * 2d problem                                        
+ * PhD thesis Gregor Gassner, pg. 99
+ *   Analytical solution for the Navier-Stokes equations
+ *                                                  
+ **********************************************************/
+
+
+namespace Dune {
+
+template <class GridType>
+class NSWaves : public EvolutionProblemInterface<
+                  Dune::Fem::FunctionSpace< double, double, GridType::dimension, GridType::dimension + 2 >,
+                  true >,
+                public Thermodynamics< GridType::dimensionworld >
+{
+  NSWaves( const NSWaves& );
+  public:
+  typedef Fem::FunctionSpace<typename GridType::ctype,
+                             double, GridType::dimensionworld,
+                             GridType::dimensionworld + 2 > FunctionSpaceType ;
+
+  typedef Fem::Parameter  ParameterType ;
+
+  enum{ dimension = GridType::dimensionworld };
+  enum { energyId = dimension + 1 };
+  typedef typename FunctionSpaceType :: DomainFieldType   DomainFieldType;
+  typedef typename FunctionSpaceType :: DomainType        DomainType;
+  typedef typename FunctionSpaceType :: RangeFieldType    RangeFieldType;
+  typedef typename FunctionSpaceType :: RangeType         RangeType;
+  typedef Thermodynamics< dimension >                     ThermodynamicsType;
+
+  NSWaves() : ThermodynamicsType(),
+    myName_( "NSWaves" ),
+    omegaGNS_( ParameterType::getValue< double >( "omegaGNS" ) ),
+    kGNS_( ParameterType::getValue< double >( "kGNS" ) ),
+    gammaGNS_( ParameterType::getValue< double >( "gammaGNS" ) ),
+    endTime_ ( ParameterType::getValue<double>( "femhowto.endTime" )), 
+    mu_( ParameterType :: getValue< double >( "mu" )), 
+    k_ ( c_pd() * Pr_inv() * mu_),
+    A_( init( true ) ),
+    B_( init( false ) )
+  {
+  }
+
+
+  // initialize A and B 
+  double init(const bool returnA ) const ;
+
+  // print info 
+  void printInitInfo() const;
+
+  // source implementations 
+  inline bool hasStiffSource() const { return false; }
+  inline bool hasNonStiffSource() const  { return true; }
+  inline double stiffSource( const double t, const DomainType& x, const RangeType& u, RangeType& res ) const;
+  inline double nonStiffSource( const double t, const DomainType& x, const RangeType& u, RangeType& res ) const;
+
+  // this is the initial data
+  inline void evaluate( const DomainType& arg , RangeType& res ) const 
+  {
+    evaluate( 0., arg, res );
+  }
+
+  // evaluate function 
+  inline void evaluate( const double t, const DomainType& x, RangeType& res ) const;
+
+  // cloned method 
+  inline void evaluate( const DomainType& x, const double t, RangeType& res ) const
+  {
+    evaluate( t, x, res );
+  }
+
+  //template< class RangeImp >
+  double pressure( const RangeType& u ) const 
+  {
+    return thermodynamics().pressureEnergyForm( u );
+  }
+
+  double 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;
+
+
+  /*  \brief finalize the simulation using the calculated numerical
+   *  solution u for this problem
+   *  
+   *  \param[in] variablesToOutput Numerical solution in the suitably chosen variables
+   *  \param[in] eocloop Specific EOC loop
+   */
+  template< class DiscreteFunctionType >
+  void finalizeSimulation( DiscreteFunctionType& variablesToOutput,
+                           const int eocloop) const
+  {}
+
+
+  const ThermodynamicsType& thermodynamics() const { return *this; }
+  using ThermodynamicsType :: Re ;
+  using ThermodynamicsType :: Re_inv;
+  using ThermodynamicsType :: Pr;
+  using ThermodynamicsType :: Pr_inv;
+  using ThermodynamicsType :: c_pd; 
+  using ThermodynamicsType :: c_pd_inv;
+  using ThermodynamicsType :: c_vd; 
+  using ThermodynamicsType :: c_vd_inv;
+  using ThermodynamicsType :: gamma;
+  using ThermodynamicsType :: g;
+  using ThermodynamicsType :: R_d_inv;
+  void printmyInfo( std::string filename ) const {}
+  inline double endtime() const { return endTime_; }
+  inline std::string myName() const { return myName_; }
+  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(); } 
+  
+protected:
+  const std::string myName_;
+  const double omegaGNS_;
+  const double kGNS_;
+  const double gammaGNS_;
+  const double endTime_;
+  const double mu_;
+  const double k_;
+  const double A_;
+  const double B_;
+};
+
+
+template <class GridType>
+inline double NSWaves<GridType>
+:: init(const bool returnA ) const 
+{
+  if( dimension == 1 )
+  {
+    if( returnA ) // A 
+      return (-omegaGNS_ + kGNS_*(3.5*gamma()-2.5));
+    else // B 
+      return ( kGNS_*(0.5+3.5*gamma()) - 4.*omegaGNS_ );
+  }
+  if( dimension == 2 ) 
+  {
+    if( returnA ) // A 
+      return (-omegaGNS_ + kGNS_*(3.*gamma() - 1.));
+    else // B 
+      return ( kGNS_*(2.+6.*gamma()) - 4.*omegaGNS_ );
+  }
+  else if( dimension == 3 ) 
+  {
+    if( returnA ) // A 
+      return (-omegaGNS_ + 0.5*kGNS_*(5.*gamma() + 1.));
+    else // B 
+     return (kGNS_*(4.5+7.5*gamma()) - 4.*omegaGNS_);
+  }
+
+  abort();
+  return 0;
+}
+
+
+
+template <class GridType>
+inline void NSWaves<GridType>
+:: printInitInfo() const
+{}
+
+
+template <class GridType>
+inline double NSWaves<GridType>
+:: stiffSource( const double t, const DomainType& x, const RangeType& u, RangeType& res ) const
+{
+  res = 0.;
+  return 0.;
+}
+
+
+template <class GridType>
+inline double NSWaves<GridType>
+:: nonStiffSource( const double t, const DomainType& x, const RangeType& u, RangeType& res ) const
+{
+  /*
+  res = 0;
+  double sumX = 0;
+  for( int i=0; i< dimension; ++i ) sumX += x[i];
+
+  double Frequency=1. ;
+  double Amplitude=0.1 ;
+  double Pi = M_PI;
+  double Omega=Pi*Frequency ;
+  double a=1.0*2.*Pi;
+  double Kappa = gamma();
+  double KappaM1 = gamma() - 1.0;
+  double mu0 = mu_ ;
+
+  res[ 0 ] = (-a+3*Omega)*std::cos(Omega*sumX-a*t);
+  res[ 1 ] = (-a+Omega/(2.)*(1.+Kappa*5.))*std::cos(Omega*sumX-a*t)
+              + Amplitude*Omega*KappaM1*std::sin(2.*(Omega*sumX-a*t));
+  for(int i=2; i<energyId; ++i )
+    res[ i ] = res[ 1 ];
+
+  res[ energyId ]  = 0.5*((9.+Kappa*15.)*Omega-8.*a)*std::cos(Omega*sumX-a*t) 
+              +Amplitude*(3.*Omega*Kappa-a)*std::sin(2.*(Omega*sumX-a*t))
+              +3.*mu0*Kappa*(Omega*Omega)*Pr_inv()*std::sin(Omega*sumX-a*t);
+
+  res *= Amplitude;
+  */
+
+  double sumX = 0;
+  for( int i=0; i< dimension; ++i ) sumX += x[i];
+  const double beta = kGNS_ * sumX - omegaGNS_*t;
+  const double cosBeta = std::cos( beta );
+  const double sinBeta = std::sin( beta );
+  const double sin2BetaGamma = std::sin( 2.*beta ) * gammaGNS_;
+  const double sinGammaKappa = sin2BetaGamma * kGNS_ * (gamma() - 1.);
+
+  res[0] = gammaGNS_*( cosBeta * ( dimension * kGNS_ - omegaGNS_) );
+  res[1] = gammaGNS_*( cosBeta * A_ + sinGammaKappa );
+
+  // set other velocity components
+  res[2] = res[ 1 ];
+  res[ energyId - 1 ] = res[ 1 ];
+
+  res[ energyId ] = gammaGNS_*( sin2BetaGamma*( dimension * gamma()*kGNS_ - omegaGNS_ ) );
+  res[ energyId ] += gammaGNS_*( cosBeta *B_ );
+  res[ energyId ] += dimension * gammaGNS_ * kGNS_ * kGNS_* k_ * c_vd_inv() 
+                     * Re_inv() * sinBeta;
+
+  // time step restriction
+  return 0.0;
+}
+
+
+
+template <class GridType>
+inline void NSWaves<GridType>
+:: evaluate( const double t, const DomainType& x, RangeType& res ) const 
+{
+#ifdef WBPROBLEM
+  // 10-0.7*(10.*x-x*x*2.) , 10.-x*4.
+  res = 0.;
+  double z = x[dimension-1];
+  res[0] = 10.-4.*z;
+  double p = thermodynamics_.p0() - g_*(10.*z - 2.*z*z);
+  res[energyId] = p/(gamma()-1.);
+  for (int i=0;i<dimension;++i)
+    res[energyId] += 0.5*res[i+1]*res[i+1]/res[0];
+#else
+  /*
+  res = 0;
+
+  double sumX = 0;
+  for( int i=0; i< dimension; ++i ) sumX += x[i];
+
+  double Frequency=1. ;
+  double Amplitude=0.1 ;
+  double Pi = M_PI;
+  double Omega=Pi*Frequency ;
+  double a=1.0*2.*Pi;
+
+  res[ 0 ] = 2.+ Amplitude*std::sin(Omega*sumX - a*t);
+  for( int i=1; i<=energyId; ++i )
+    res[ i ] = res[ 0 ];
+  res[ energyId ] *= res[ 0 ];
+  */
+
+  double sumX = 0;
+  for( int i=0; i< dimension; ++i ) sumX += x[i];
+  const double beta = kGNS_ * sumX - omegaGNS_*t;
+  const double sinBeta = std::sin( beta );
+  const double sinGamma2 = sinBeta * gammaGNS_ + 2.;
+
+  for(int i=0; i<energyId; ++i) 
+  {
+    res[ i ] = sinGamma2; // constant velocity field
+  }
+  res[ energyId ] = sinGamma2 * sinGamma2;
+#endif
+}
+
+
+
+template <class GridType>
+template< class RangeImp >
+inline void NSWaves<GridType>
+:: pressAndTemp( const RangeImp& u, double& p, double& T ) const
+{
+  thermodynamics().pressAndTempEnergyForm( u, p, T );
+}
+
+
+template <class GridType>
+inline std::string NSWaves<GridType>
+:: description() const
+{
+  std::ostringstream stream;
+
+  stream <<"{\\bf Problem:}" <<myName_
+         <<", {\\bf $\\mu$:} " <<mu_
+         <<", {\\bf End time:} " <<endTime_
+         <<", {$\\gamma_{GNS}$:} " <<gammaGNS_ 
+         <<"\n"
+         <<", {$\\omega_{GNS}$:} " <<omegaGNS_
+         <<", {$k_{GNS}$:} " <<kGNS_;
+
+  std::string returnString = stream.str();
+
+  return returnString;
+}
+
+} // end namespace Dune
+#endif

dune/fem-dg/test/navierstokes/passtraits.hh

+#ifndef DUNE_FEM_DG_PASSTRAITS_HH
+#define DUNE_FEM_DG_PASSTRAITS_HH
+
+// Dune includes
+#include <dune/fem/gridpart/common/gridpart.hh>
+#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
+
+// Dune-Fem includes
+#include <dune/fem/space/finitevolume.hh>
+#include <dune/fem/space/discontinuousgalerkin.hh>
+#include <dune/fem/function/adaptivefunction.hh>
+#include <dune/fem/quadrature/cachingquadrature.hh>
+#include <dune/fem-dg/operator/adaptation/adaptation.hh>
+#include <dune/fem/function/petscdiscretefunction.hh>
+
+
+namespace Dune {
+
+  //PassTraits
+  //----------
+
+  template <class Model,int dimRange,int polOrd>
+  class PassTraits
+  {
+  public:
+    typedef typename Model :: Traits                                 ModelTraits;
+    typedef typename ModelTraits :: GridPartType                     GridPartType;
+    typedef typename GridPartType :: GridType                        GridType;
+    typedef typename GridType :: ctype                               ctype;
+    static const int dimDomain = Model :: Traits :: dimDomain;
+
+    typedef Fem::CachingQuadrature< GridPartType, 0 >                     VolumeQuadratureType;
+    //typedef Fem::LagrangePointQuadrature< GridPartType, 2*polOrd+4 >          VolumeQuadratureType;
+    typedef Fem::CachingQuadrature< GridPartType, 1 >                     FaceQuadratureType;
+
+    //typedef Fem::ElementQuadrature< GridPartType, 0 >                     VolumeQuadratureType;
+    //typedef Fem::ElementQuadrature< GridPartType, 1 >                     FaceQuadratureType;
+
+    // Allow generalization to systems
+    typedef Fem::FunctionSpace< ctype, double, dimDomain, dimRange >      FunctionSpaceType;
+#ifdef LEGENDREBASIS
+#warning "Using Legendre basis functions"
+    typedef Fem::LegendreDiscontinuousGalerkinSpace< FunctionSpaceType,
+                                        GridPartType, polOrd,
+                                        Fem::CachingStorage >             DiscreteFunctionSpaceType;
+#elif defined LAGRANGEBASIS
+#warning "Using Lagrange basis functions"
+    typedef Fem::LagrangeDiscontinuousGalerkinSpace< FunctionSpaceType,
+                                        GridPartType, polOrd,
+                                        Fem::CachingStorage >             DiscreteFunctionSpaceType;
+#else
+#warning "Using standard modal basis functions"
+    typedef Fem::DiscontinuousGalerkinSpace< FunctionSpaceType,
+                                        GridPartType, polOrd,
+                                        Fem::CachingStorage >             DiscreteFunctionSpaceType;
+#endif
+    typedef Fem::AdaptiveDiscreteFunction< DiscreteFunctionSpaceType >    DestinationType;
+    //typedef Fem::PetscDiscreteFunction< DiscreteFunctionSpaceType >    DestinationType;
+
+    // Indicator for Limiter
+    typedef Fem::FunctionSpace< ctype, double, dimDomain, 3> FVFunctionSpaceType;
+    typedef Fem::FiniteVolumeSpace<FVFunctionSpaceType,GridPartType, 0, Fem::SimpleStorage> IndicatorSpaceType;
+    typedef Fem::AdaptiveDiscreteFunction<IndicatorSpaceType> IndicatorType;
+
+    typedef AdaptationHandler< GridType, FunctionSpaceType >  AdaptationHandlerType ;
+  };
+
+} // end namespace Dune 
+#endif

dune/fem-dg/test/navierstokes/problemcreator.hh

+#ifndef FEMHOWTO_NSEQ_RPOBLEMCREATOR_HH
+#define FEMHOWTO_NSEQ_RPOBLEMCREATOR_HH
+#include <config.h>
+
+// system includes
+#include <string>
+
+#include "passtraits.hh"
+
+// dune includes
+#include <dune/common/version.hh>
+#include <dune/fem/io/parameter.hh>
+
+// local includes
+#include <dune/fem-dg/operator/fluxes/eulerfluxes.hh>
+#include <dune/fem-dg/operator/fluxes/diffusionflux.hh>
+#include "problemtype.hh"
+#include "ns_model.hh"
+
+
+template< class GridType > 
+struct ProblemGenerator 
+{
+  typedef NSProblemType ProblemType;
+
+  template< class GridPart >
+  struct Traits
+  {
+    typedef ProblemType InitialDataType;
+    typedef Dune::NSModel< GridPart, InitialDataType > ModelType;
+    // choice of diffusion flux (see diffusionflux.hh for methods)
+    static const Dune :: DGDiffusionFluxIdentifier PrimalDiffusionFluxId 
+           = Dune :: method_general ;
+
+// ******************************** NUMERICAL FLUX *****************************
+#if (FLUX==1)
+#warning "FLUX: LLF"
+    typedef LLFFlux< ModelType > FluxType;
+#elif (FLUX==2)
+#warning "FLUX: HLL (Dennis)"
+    typedef HLLNumFlux< ModelType > FluxType;
+#elif (FLUX==3)
+#warning "FLUX: HLLC (Dennis)"
+    typedef HLLCNumFlux< ModelType > FluxType;
+#elif (FLUX==4)
+#warning "FLUX: HLL2C"
+    typedef HLL2CFlux< ModelType > FluxType;
+#elif (FLUX==5)
+#warning "FLUX: HLL2"
+    typedef HLL2Flux< ModelType > FluxType;
+#elif (FLUX==6)
+#warning "FLUX: HLLEM (Mhd)"
+    typedef HLLEMNumFlux< ModelType > FluxType;
+#else
+#error "Set the flag FLUX! See Makefile.am for details!"
+#endif
+    
+// ****************************** END NUMERICAL FLUX ***************************
+  };
+
+  static inline std::string advectionFluxName()
+  {
+#if (FLUX==1)
+    return "LLF";
+#elif (FLUX==2)
+    return "HLL(Dennis)";
+#elif (FLUX==3)
+    return "HLLC(Dennis)";
+#elif (FLUX==4)
+    return "HLL2C";
+#elif (FLUX==5)
+    return "HLL2";
+#elif (FLUX==6)
+    return "HLLEM(Mhd)";
+#endif
+  }
+
+  static inline std::string diffusionFluxName()
+  {
+#ifdef EULER
+    return "";
+#elif (defined PRIMALDG)
+    return Dune::Fem::Parameter::getValue< std::string >("dgdiffusionflux.method");
+#else
+    return "LDG";
+#endif
+  }
+
+  static inline Dune::GridPtr<GridType>               
+  initializeGrid( const std::string description )
+  {
+    // use default implementation 
+    return initialize< GridType > ( description );
+  }
+
+  static ProblemType* problem()
+  {
+    // choice of explicit or implicit ode solver
+    return new ProblemType ();
+  }
+};
+
+#endif // FEMHOWTO_NSEQ_RPOBLEMCREATOR_HH

dune/fem-dg/test/navierstokes/problemtype.hh

+#ifndef FEMHOWTO_NSEQ_PROBLEMTYPE_HH
+#define FEMHOWTO_NSEQ_PROBLEMTYPE_HH
+
+#include <dune/common/version.hh>
+#include <dune/fem/io/parameter.hh>
+
+///////////////////////////////////////
+// AVAILABLE PROBLEMS
+///////////////////////////////////////
+#if PROBLEM==1
+#error
+   // a trivial exact solution
+   // for Navier-Stokes equations
+   #include "constant.hh"
+   #warning "Problem: 2d constant Navier-Stokes solution"
+   typedef NSConstant< GridSelector :: GridType >  NSProblemType;
+   #define PROBLEM_HAS_SOLUTION
+
+#elif PROBLEM==2
+   // a nonstationary exact solution 
+   // for Navier-Stokes equations
+   #include "nswaves.hh"
+   typedef NSWaves< GridSelector :: GridType >  NSProblemType;
+   #define PROBLEM_HAS_SOLUTION
+
+#elif PROBLEM==3
+   // a nonstationary exact solution 
+   // for Navier-Stokes equations
+   #include "nssmooth.hh"
+   typedef NSSmoothSolution< GridSelector :: GridType >  NSProblemType;
+   #define PROBLEM_HAS_SOLUTION
+
+#else
+   #error "No valid problem number specified"
+#endif
+
+#endif

dune/fem-dg/test/navierstokes/thermodynamics.hh

+#ifndef DUNE_THERMODYNAMICS_HH
+#define DUNE_THERMODYNAMICS_HH
+
+// system include
+#include <cmath>
+#include <iomanip>
+#include <iostream>
+#include <fstream>
+
+// dune-fem include
+#include <dune/fem/io/parameter.hh>
+
+using namespace Dune;
+
+// Thermodynamics
+// --------------
+
+/** \class Thermodynamics
+ *  \brief deals with physics used for the atmosphere, in particular
+ *         physical constants, equation of state etc.
+ *
+ *  \tparam dimDomain dimension of the domain
+ */
+template< int dimDomain >
+class Thermodynamics
+{
+  enum{ energyId = dimDomain+1 };
+
+  public:
+  typedef Fem::Parameter  ParameterType ;
+
+  Thermodynamics() :
+    Re_( ParameterType::getValue< double >( "Re" ) ),
+    Pr_( ParameterType::getValue< double >( "Pr" ) ),
+    g_( ParameterType::getValue< double >( "g", 9.81 ) ),
+    p0_(ParameterType::getValue< double >( "p0", 100000. )),
+    p0Star_( ParameterType::getValue< double >( "p0Star", 610.7 )),
+    T0_( ParameterType::getValue< double >( "T0", 273.15 )),
+    c_pd_     ( ParameterType::getValue< double >( "c_pd", 1004. )),
+    c_vd_     ( ParameterType::getValue< double >( "c_vd", 717. )),
+    c_pl_     ( ParameterType::getValue< double >( "c_pl", 4186. )),
+    L0_       ( ParameterType::getValue< double >( "L0", 2500000. )),
+    relCloud_ ( ParameterType::getValue< double >( "relCloud", 1. )),
+    Re_inv_     ( 1. / Re_ ),
+    Pr_inv_     ( 1. / Pr_ ),
+    c_pd_inv_   ( 1. / c_pd_ ),
+    c_vd_inv_   ( 1. / c_vd_ ),
+    c_pl_inv_   ( 1. / c_pl_ ),
+    R_d_        ( c_pd_ - c_vd_ ),
+    R_d_inv_    ( 1. / R_d_ ),
+    p0_inv_     ( 1. / p0_ ),
+    T0_inv_     ( 1. / T0_ ),
+    kappa_      ( R_d_ * c_pd_inv_ ),
+    kappa_inv_  ( 1. / kappa_ ),
+    gamma_      ( c_pd_ * c_vd_inv_ ),
+    gammaM1_    ( gamma_ - 1.0 ),
+    gamma_inv_  ( 1. / gamma_ )
+  {
+    assert( gamma_ > 1. );
+    assert( R_d_ > 200. );
+  }
+
+  /** \brief calculate the pressure and the temperature assuming the energy form
+   *         for conservative variables: \f$[\rho,\rho\boldsymbol{v},\rho e]\f$
+   *
+   *  Calculate the pressure and temperature from the conservative variables
+   *  \f$[\rho,\rho\boldsymbol{v},\rho e]\f$, where \f$ e \f$ is the sum of 
+   *  the internal and the kinetic energy.
+   *
+   *  \param[in] cons Conervative variables
+   *
+   *  \return pressure in energy form 
+   */
+  template< class RangeType >
+  inline double 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.;
+    for( int i=1; i<=dimDomain; ++i )
+      kin += cons[ i ] * cons[ i ];
+    kin *= 0.5 / cons[ 0 ];
+
+    return gammaM1_ * ( cons[ energyId ] - kin );
+  }
+
+  /** \brief calculate the pressure and the temperature assuming the energy form
+   *         for conservative variables: \f$[\rho,\rho\boldsymbol{v},\rho e]\f$
+   *
+   *  Calculate the pressure and temperature from the conservative variables
+   *  \f$[\rho,\rho\boldsymbol{v},\rho e]\f$, where \f$ e \f$ is the sum of 
+   *  the internal and the kinetic energy.
+   *
+   *  \param[in] cons Conervative variables
+   *  \param[out] p Pressure
+   *  \param[out] T temperature
+   *
+   *  \tparam RangeType Type of the range value
+   */
+  template< class RangeType >
+  inline double temperatureEnergyForm( const RangeType& cons, 
+                                       const double p ) const 
+  {
+    assert( cons[0] > 1e-20 );
+    return R_d_inv_ * p / cons[ 0 ];
+  }
+
+  /** \brief calculate the pressure and the temperature assuming the energy form
+   *         for conservative variables: \f$[\rho,\rho\boldsymbol{v},\rho e]\f$
+   *
+   *  Calculate the pressure and temperature from the conservative variables
+   *  \f$[\rho,\rho\boldsymbol{v},\rho e]\f$, where \f$ e \f$ is the sum of 
+   *  the internal and the kinetic energy.
+   *
+   *  \param[in] cons Conervative variables
+   *
+   *  \return pressure in energy form 
+   */
+  template< class RangeType >
+  inline double temperatureEnergyForm( const RangeType& cons ) const
+  {
+    return temperatureEnergyForm( cons, pressureEnergyForm( cons ) );
+  }
+
+  /** \brief calculate the pressure and the temperature assuming the energy form
+   *         for conservative variables: \f$[\rho,\rho\boldsymbol{v},\rho e]\f$
+   *
+   *  \param[in] cons Conervative variables
+   *  \param[out] p Pressure
+   *  \param[out] T temperature
+   *
+   *  \tparam RangeType Type of the range value
+   */
+  template< class RangeType >
+  inline double densityThetaForm( const RangeType& prim ) const 
+  {
+    const double p     = prim[ energyId - 1 ];
+    const double 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 ;
+
+    assert( rho > 0.0 );  
+    return rho;
+  }
+
+
+  template< class RangeType >
+  void conservativeToPrimitiveThetaForm( const RangeType& cons, RangeType& prim ) const;
+
+  template< class RangeType >
+  void primitiveToConservativeThetaForm( const RangeType& prim, RangeType& cons ) const;
+
+  template< class RangeType >
+  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_; }
+
+  //! the Reynolds number Re
+  inline double Re() const { return Re_; }
+  //! the inverse Reynolds number (1/Re)
+  inline double Re_inv() const { return Re_inv_; }
+  //! the Prandtl number Pr 
+  inline double Pr() const { return Pr_; }
+  //! the inverse Prandtl number (1/Pr) 
+  inline double 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_;
+};
+
+
+/** \brief converts conservative variables in the energy form to primitive ones
+ *
+ *  Converts conservative variables \f$[\rho\boldsymbol{v},p,\rho e\f$
+ *  to primitive ones \f$[\boldsymbol{v},p,\theta]\f$, where \f$e\f$ is the sum
+ *  of internal and kinetic energy, and \f$\theta\f$ potential temperature.
+ *
+ *  \param[in] cons Conservative variables
+ *  \param[out] prim Primitive variables
+ *
+ *  \tparam dimDomain dimension of the domain
+ *  \tparam RangeType type of the range value
+ */
+template< int dimDomain >
+template< class RangeType >
+void Thermodynamics< dimDomain >
+:: conservativeToPrimitiveEnergyForm( const RangeType& cons, RangeType& prim ) const
+{
+  std::cerr <<"conservativeToPrimitiveEnergyForm not implemented" <<std::endl;
+  abort();
+
+  //const double rho_inv = 1./ cons[0];
+
+  //double p, T;
+  //pressAndTempEnergyForm( cons, p, T );
+
+  //prim[energyId-1] = p;
+  // this is not pot. temp !!!!!!!
+  //prim[energyId] = cons[energyId]/cons[0];
+}
+
+
+/** \brief converts conservative variables in the theta form to primitive ones
+ *
+ *  Converts conservative variables \f$[\rho\boldsymbol{v},p,\rho\theta]\f$
+ *  to primitive ones \f$[\boldsymbol{v},p,\theta]\f$, where \f$\theta\f$ is
+ *  potential temperature
+ *
+ *  \param[in] cons Conservative variables
+ *  \param[out] prim Primitive variables
+ *
+ *  \tparam dimDomain dimension of the domain
+ *  \tparam RangeType type of the range value
+ */
+template< int dimDomain >
+template< class RangeType >
+void Thermodynamics< dimDomain >
+:: conservativeToPrimitiveThetaForm( const RangeType& cons, RangeType& prim ) const
+{
+  assert( cons[0] > 0. );
+  assert( cons[energyId] > 0. );
+
+  double p, T;
+  pressAndTempThetaForm( cons, p, T );
+
+  for( int i = 0; i < dimDomain; ++i )
+    prim[i] = cons[i+1]/cons[0];
+
+  prim[energyId-1] = p;
+  prim[energyId] = cons[energyId] / cons[0];
+}
+
+template< int dimDomain >
+template< class RangeType >
+void Thermodynamics< dimDomain >
+:: primitiveToConservativeThetaForm( const RangeType& prim, RangeType& cons ) const
+{
+  // p,theta  --> rho 
+  cons[ 0 ] = densityThetaForm( prim );
+
+  // v_i  --> v_i+1 * rho 
+  for( int i = 0; i < dimDomain; ++i )
+    cons[ i+1 ] = prim[ i ] * cons[0];
+
+  // theta --> theta * rho 
+  cons[ energyId ] = prim[ energyId ] * cons[ 0 ];
+}
+
+#endif // file define