diff --git a/dune/fem-dg/test/navierstokes/Makefile.am b/dune/fem-dg/test/navierstokes/Makefile.am
index 767e992911387846a0e0c86fd48cd1d86a1ef402..cc0f3e0fdb5c486ed7c51b09c1c7a1caab345f8f 100644
--- a/dune/fem-dg/test/navierstokes/Makefile.am
+++ b/dune/fem-dg/test/navierstokes/Makefile.am
@@ -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 
diff --git a/dune/fem-dg/test/navierstokes/ns_model.hh b/dune/fem-dg/test/navierstokes/ns_model.hh
index 68fc9cd0c30d866fdba4200fd70bb9ec86cab600..c0c79e50cbe176bcdcf6a7fa31d5c6b92d7af87b 100644
--- a/dune/fem-dg/test/navierstokes/ns_model.hh
+++ b/dune/fem-dg/test/navierstokes/ns_model.hh
@@ -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 );
   }
diff --git a/dune/fem-dg/test/navierstokes/nswaves.hh b/dune/fem-dg/test/navierstokes/nswaves.hh
index c495722d4721e38eac14fe66a623d8d69cabe13e..3cfec440a8d545c8768b5ef86967e562dcd418bb 100644
--- a/dune/fem-dg/test/navierstokes/nswaves.hh
+++ b/dune/fem-dg/test/navierstokes/nswaves.hh
@@ -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 );
 }
diff --git a/dune/fem-dg/test/navierstokes/parameter b/dune/fem-dg/test/navierstokes/parameter
index 708c9fbdfa605ec8ef51cb70e6aba9866bdebc6c..cfdc9376d485dc667624c69e9c48bb619f2e8510 100644
--- a/dune/fem-dg/test/navierstokes/parameter
+++ b/dune/fem-dg/test/navierstokes/parameter
@@ -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.
diff --git a/dune/fem-dg/test/navierstokes/thermodynamics.hh b/dune/fem-dg/test/navierstokes/thermodynamics.hh
index 7edb6d54bebabfe37cfb020933343cd4119d5830..c9b4d7c3315c9caa8fc134c8106482cdda5451f5 100644
--- a/dune/fem-dg/test/navierstokes/thermodynamics.hh
+++ b/dune/fem-dg/test/navierstokes/thermodynamics.hh
@@ -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