ns works with flop count

dune/fem-dg/main/main.cc

@@ -8,6 +8,7 @@
 #include <config.h>
 
 #if defined USE_BASEFUNCTIONSET_CODEGEN || defined BASEFUNCTIONSET_CODEGEN_GENERATE
+#warning "Codegen enabled!"
 #define USE_FEMDG_BASISFUNCTIONSET
 #endif
 

dune/fem-dg/main/main.hh

 #ifndef DUNE_FEM_DG_MAINHEADER_HH
 #define DUNE_FEM_DG_MAINHEADER_HH
 
+//#ifdef COUNT_FLOPS
+#include <dune/fem/misc/double.hh>
+//#endif
+
 // streams for backup
 #include <dune/fem-dg/misc/streams.hh>
 

dune/fem-dg/misc/cons2prim.hh

@@ -27,6 +27,7 @@ void setupAdditionalVariables( const TimeProvider& tp,
                                const ModelType& model,
                                PrimDiscreteFunctionType& primDF )
 {
+#if 0
   typedef typename ConsDiscreteFunctionType::Traits::DiscreteFunctionSpaceType
     ConsDiscreteFunctionSpaceType;
   typedef typename PrimDiscreteFunctionType::Traits::DiscreteFunctionSpaceType
@@ -78,6 +79,7 @@ void setupAdditionalVariables( const TimeProvider& tp,
       primLF.axpy( quad[qP] , prim );
     }
   }
+#endif
 }
 
 #endif

dune/fem-dg/misc/crs.hh

@@ -20,7 +20,7 @@ namespace Dune
   template<class K, int ROWS, int COLS>
   class BlockCRSMatrix
   {
-    std::vector< double > data_;
+    std::vector< K > data_;
     std::vector< int > rows_;
     std::vector< int > cols_;
   public:

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

@@ -95,7 +95,7 @@ inline
 EulerAnalyticalFlux<2>::FieldType EulerAnalyticalFlux<2>::rhoeps(const RangeType& u) const
 {
   assert( u[0] >= 1e-10 );
-  return u[e]-0.5*(u[1]*u[1]+u[2]*u[2])/u[0];
+  return u[e] -(0.5*((u[1]*u[1])+(u[2]*u[2]))/u[0]);
 }
 
 template <>
@@ -105,7 +105,7 @@ EulerAnalyticalFlux<2>::FieldType EulerAnalyticalFlux<2>::pressure(const FieldTy
 {
   const FieldType re = rhoeps(u);
   assert(re>=1e-10);
-  return (gamma-1)*re;
+  return (gamma-1.0)*re;
 }
 
 template <>
@@ -222,12 +222,14 @@ EulerAnalyticalFlux<2>::maxSpeed(const FieldType gamma,
                                  const DomainType& n,
                                  const RangeType& u) const
 {
-  assert( u[0] > 1e-10 );
-  FieldType u_normal = (u[1]*n[0]+u[2]*n[1]) / u[0];
+  assert( u[0] > FieldType( 1e-10 ) );
+  FieldType u_normal = FieldType(u[1]*n[0]+u[2]*n[1]) / u[0];
   FieldType p = pressure(gamma,u);
   FieldType c2 = gamma * p/ u[0] * n.two_norm2();
-  assert( c2 > 1e-10 );
-  return std::abs(u_normal) + std::sqrt(c2);
+  assert( c2 > FieldType( 1e-10 ) );
+  FieldType maxspd = std::abs(double(u_normal)) + std::sqrt(double(c2));
+  //std::cout << maxspd << " maxspd" << std::endl;
+  return maxspd;
 }
 
 template <>

dune/fem-dg/solver/smartodesolver.hh

@@ -12,7 +12,7 @@ namespace Dune {
 struct SmartOdeSolverParameters : public DuneODE :: ODEParameters
 {
   SmartOdeSolverParameters( const std::string keyPrefix = "fem.ode." )
-    : DuneODE :: ODEParameters( keyPrefix )
+    : DuneODE :: ODEParameters()
   {}
 
   using DuneODE :: ODEParameters :: keyPrefix_;

dune/fem-dg/stepper/advectiondiffusionstepper.hh

@@ -71,16 +71,18 @@ struct AdvectionDiffusionStepper
   using BaseType :: space;
   using BaseType :: problem;
   using BaseType :: adaptationHandler_ ;
+  using BaseType :: adaptParam_;
   using BaseType :: adaptive ;
   using BaseType :: doEstimateMarkAdapt ;
+  using BaseType :: name ;
 
   AdvectionDiffusionStepper( GridType& grid, const std::string name = "" ) :
     BaseType( grid, name ),
     dgOperator_( gridPart_, problem(), name  ),
     dgAdvectionOperator_( gridPart_, problem(), name ),
     dgDiffusionOperator_( gridPart_, problem(), name ),
-    dgIndicator_( gridPart_, problem(), name )
-    // gradientIndicator_( space(), problem(), adaptParam_ )
+    dgIndicator_( gridPart_, problem(), name ),
+    gradientIndicator_( space(), problem(), adaptParam_ )
   {
   }
 
@@ -102,8 +104,11 @@ struct AdvectionDiffusionStepper
     }
 
     // one of them is not zero,
+    size_t advSize     = dgAdvectionOperator_.numberOfElements();
+    size_t diffSize    = dgDiffusionOperator_.numberOfElements();
+    size_t dgIndSize   = gradientIndicator_.numberOfElements();
     size_t dgSize      = dgOperator_.numberOfElements();
-    UInt64Type grSize  = dgSize; // std::max( std::max(advSize, dgSize ), std::max( diffSize, dgIndSize ) );
+    UInt64Type grSize  = std::max( std::max(advSize, dgSize ), std::max( diffSize, dgIndSize ) );
     double minMax[ 2 ] = { double(grSize), 1.0/double(grSize) } ;
     grid_.comm().max( &minMax[ 0 ], 2 );
     if( Dune::Fem::Parameter :: verbose () )
@@ -118,7 +123,7 @@ struct AdvectionDiffusionStepper
     // create adaptation handler in case of apost indicator
     if( adaptive() )
     {
-      if( ! adaptationHandler_ ) // && adaptParam_.aposterioriIndicator() )
+      if( ! adaptationHandler_ && adaptParam_.aposterioriIndicator() )
       {
         adaptationHandler_ = new AdaptationHandlerType( grid_, tp );
         dgIndicator_.setAdaptation( *adaptationHandler_ );
@@ -130,19 +135,20 @@ struct AdvectionDiffusionStepper
                               LinearInverseOperatorType > OdeSolverImpl;
     return new OdeSolverImpl( tp, dgOperator_,
                               dgAdvectionOperator_,
-                              dgDiffusionOperator_ );
+                              dgDiffusionOperator_,
+                              name() );
   }
 
   //! estimate and mark solution
   virtual void initialEstimateMarkAdapt( )
   {
-    // doEstimateMarkAdapt( dgIndicator_, gradientIndicator_, true );
+    doEstimateMarkAdapt( dgIndicator_, gradientIndicator_, true );
   }
 
   //! estimate and mark solution
   virtual void estimateMarkAdapt( )
   {
-    // doEstimateMarkAdapt( dgIndicator_, gradientIndicator_, false );
+    doEstimateMarkAdapt( dgIndicator_, gradientIndicator_, false );
   }
 
   const ModelType& model() const { return dgOperator_.model(); }
@@ -152,6 +158,6 @@ protected:
   ExplicitOperatorType    dgAdvectionOperator_;
   ImplicitOperatorType    dgDiffusionOperator_;
   DGIndicatorType         dgIndicator_;
-  //GradientIndicatorType   gradientIndicator_;
+  GradientIndicatorType   gradientIndicator_;
 };
 #endif // FEMHOWTO_STEPPER_HH

dune/fem-dg/stepper/stepperbase.hh

@@ -122,6 +122,9 @@ struct StepperBase
 
   using BaseType :: grid_ ;
   using BaseType :: limitSolution ;
+  using BaseType :: adaptParam_ ;
+  using BaseType :: eocParam_;
+  using BaseType :: param_;
 
   // constructor taking grid
   StepperBase(GridType& grid, const std::string name = "" ) :
@@ -183,12 +186,10 @@ struct StepperBase
     return IOTupleType( &solution_, additionalVariables_, indicator() );
   }
 
-  /*
   bool checkDofsValid ( TimeProviderType& tp, const int loop ) const
   {
     return solution_.dofsValid();
   }
-  */
 
   void writeData( TimeProviderType& tp, const bool writeAnyway = false )
   {
@@ -207,7 +208,7 @@ struct StepperBase
     problem().postProcessTimeStep( grid_, solution(), tp.time() );
   }
 
-  bool adaptive () const { return false; }
+  bool adaptive () const { return adaptParam_.adaptive(); }
 
   // function creating the ode solvers
   virtual OdeSolverType* createOdeSolver( TimeProviderType& ) = 0;
@@ -290,10 +291,10 @@ struct StepperBase
     {
       // copy data tuple
       dataTuple_ = dataTuple () ;
-      //dataWriter_ = new DataWriterType( grid_, dataTuple_, tp ); //,
-      //  eocParam_.dataOutputParameters( loop, problem_->dataPrefix() ) );
-      //assert( dataWriter_ );
+      dataWriter_ = new DataWriterType( grid_, dataTuple_, tp,
+        eocParam_.dataOutputParameters( loop, problem_->dataPrefix() ) );
     }
+    assert( dataWriter_ );
 
     typedef typename InitialDataType :: TimeDependentFunctionType
       TimeDependentFunctionType;
@@ -357,8 +358,7 @@ struct StepperBase
 
   inline double error(TimeProviderType& tp, DiscreteFunctionType& u)
   {
-    // const int order = eocParam_.quadOrder() < 0 ? 2*u.space().order()+4 : eocParam_.quadOrder();
-    const int order = 2*u.space().order()+4 ;
+    const int order = eocParam_.quadOrder() < 0 ? 2*u.space().order()+4 : eocParam_.quadOrder();
     Fem :: L2Norm< GridPartType > l2norm( u.space().gridPart(), order );
     return l2norm.distance( problem().fixedTimeFunction( tp.time() ), u );
   }
@@ -421,7 +421,7 @@ protected:
                             GradientIndicator& gradientIndicator,
                             const bool initialAdaptation = false )
   {
-    if( adaptive() )
+    if( adaptParam_.adaptive() )
     {
       // get grid sequence before adaptation
       const int sequence = space().sequence();
@@ -434,7 +434,6 @@ protected:
         // do marking and adaptation
         adaptationHandler_->adapt( adaptationManager(), initialAdaptation );
       }
-/*
       else if( adaptParam_.gradientBasedIndicator() )
       {
         gradientIndicator.estimateAndMark( solution_ );
@@ -445,7 +444,7 @@ protected:
         // marking has been done by limiter
         adaptationManager().adapt();
       }
-*/
+
       // if grid has changed then limit solution again
       if( sequence != space().sequence() )
       {

dune/fem-dg/test/grids/nswaves_2d.dgf

@@ -3,7 +3,7 @@ DGF
 Interval
 0 0
 2 2
-16 16 
+32 32 
 #
 
 GridParameter

dune/fem-dg/test/navierstokes/parameter

@@ -11,7 +11,7 @@
 
 fem.prefix: data
 fem.io.datafileprefix: NS
-fem.io.outputformat: vtk-cell
+fem.io.outputformat: none
 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

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

@@ -143,7 +143,7 @@ class Thermodynamics
     assert( p > 1e-12 );
     assert( theta > 1e-12 );
 
-    const Field rho = std::pow( p/p0_ , gamma_inv_ ) * p0_ * R_d_inv_ / theta ;
+    const Field rho = std::pow( double(p/p0_) , double(gamma_inv_) ) * p0_ * R_d_inv_ / theta ;
 
     assert( rho > 0.0 );
     return rho;