new Setup class base for better encapsulation

src/dataBlock/coarsen.cpp

@@ -41,15 +41,15 @@ void DataBlock::EnrollGridCoarseningLevels(GridCoarseningFunc func) {
 void DataBlock::ComputeGridCoarseningLevels() {
   idfx::pushRegion("DataBlock::ComputeGridCoarseningLevels");
   static bool levelsHaveBeenComputedOnce = false;
-  if((gridCoarseningFunc == NULL) && (haveGridCoarsening == GridCoarsening::dynamic)) {
+  if(!gridCoarseningFunc && (haveGridCoarsening == GridCoarsening::dynamic)) {
     IDEFIX_ERROR("Dynamic grid Coarsening is enabled, "
                  "but no function has been enrolled to compute coarsening levels");
   }
   // if grid coarsening is enabled(=static), we compute the levels once
   // levels can be either initialised with the initial conditions, or with a dedicated
   // Coarsening function (if Enrollment has been called)
   if((haveGridCoarsening == GridCoarsening::enabled) && (!levelsHaveBeenComputedOnce)) {
-    if(gridCoarseningFunc != NULL) {
+    if(gridCoarseningFunc) {
       idfx::pushRegion("User-defined Coarsening function");
         gridCoarseningFunc(*this);
       idfx::popRegion();

src/dataBlock/dataBlock.cpp

@@ -386,31 +386,15 @@ void DataBlock::DeriveVectorPotential() {
 void DataBlock::LaunchUserStepLast() {
   if(haveUserStepLast) {
     idfx::pushRegion("User::UserStepLast");
-    if(userStepLast != nullptr)
-      userStepLast(*this, this->t, this->dt);
-    else
-      IDEFIX_ERROR("UserStepLast not properly initialized");
+    userStepLast(*this, this->t, this->dt);
     idfx::popRegion();
   }
 }
 
 void DataBlock::LaunchUserStepFirst() {
   if(haveUserStepFirst) {
     idfx::pushRegion("User::UserStepFirst");
-    if(userStepFirst != nullptr)
-      userStepFirst(*this, this->t, this->dt);
-    else
-      IDEFIX_ERROR("UserStepLast not properly initialized");
+    userStepFirst(*this, this->t, this->dt);
     idfx::popRegion();
   }
 }
-
-void DataBlock::EnrollUserStepLast(StepFunc func) {
-  haveUserStepLast = true;
-  userStepLast = func;
-}
-
-void DataBlock::EnrollUserStepFirst(StepFunc func) {
-  haveUserStepFirst = true;
-  userStepFirst = func;
-}

src/dataBlock/dataBlock.hpp

@@ -17,11 +17,12 @@
 #ifndef DATABLOCK_DATABLOCK_HPP_
 #define DATABLOCK_DATABLOCK_HPP_
 
-#include <vector>
-#include <string>
+#include <cstdio>
+#include <functional>
 #include <map>
 #include <memory>
-#include <cstdio>
+#include <string>
+#include <vector>
 
 #include "idefix.hpp"
 #include "grid.hpp"
@@ -50,15 +51,13 @@ class PlanetarySystem;
 template<typename Phys>
 class Fluid;
 class SubGrid;
-class Vtk;
-class Dump;
 
 // Forward class hydro declaration
 #include "physics.hpp"
 
-using GridCoarseningFunc = void(*) (DataBlock &);
+using GridCoarseningFunc = std::function<void(DataBlock &)>;
 
-using StepFunc = void (*) (DataBlock &, const real t, const real dt);
+using StepFunc = std::function<void(DataBlock &, const real, const real)>;
 
 class DataBlock {
  public:
@@ -101,7 +100,7 @@ class DataBlock {
 
   GridCoarsening haveGridCoarsening{GridCoarsening::disabled};
                                 ///< Is grid coarsening enabled?
-  GridCoarseningFunc gridCoarseningFunc{NULL};
+  GridCoarseningFunc gridCoarseningFunc;
                                ///< The user-defined grid coarsening level computation function
 
 
@@ -176,8 +175,16 @@ class DataBlock {
   void LaunchUserStepFirst();     ///< perform user-defined step before main integration step
   void LaunchUserStepLast();      ///< Perform user-defined step after main integration step
 
-  void EnrollUserStepFirst(StepFunc);
-  void EnrollUserStepLast(StepFunc);
+  template<typename Func>
+  void EnrollUserStepFirst(Func func) {
+    haveUserStepFirst = true;
+    userStepFirst = func;
+  }
+  template<typename Func>
+  void EnrollUserStepLast(Func func) {
+    haveUserStepLast = true;
+    userStepLast = func;
+  }
 
  private:
   void WriteVariable(FILE* , int , int *, char *, void*);
@@ -187,8 +194,8 @@ class DataBlock {
   bool haveUserStepFirst{false};
   bool haveUserStepLast{false};
 
-  StepFunc userStepFirst{nullptr};
-  StepFunc userStepLast{nullptr};
+  StepFunc userStepFirst;
+  StepFunc userStepLast;
 };
 
 #endif // DATABLOCK_DATABLOCK_HPP_

src/dataBlock/fargo.cpp

@@ -201,7 +201,7 @@ void Fargo::GetFargoVelocity(real t) {
   if(type != userdef)
     IDEFIX_ERROR("Fargo::GetFargoVelocity should not be called when fargo != userdef");
   if(velocityHasBeenComputed == false) {
-    if(fargoVelocityFunc== NULL) {
+    if(!fargoVelocityFunc) {
       IDEFIX_ERROR("No Fargo velocity function has been defined");
     }
     fargoVelocityFunc(*data, meanVelocity);

src/dataBlock/fargo.hpp

@@ -83,7 +83,7 @@ class Fargo {
   bool velocityHasBeenComputed{false};
   bool haveDomainDecomposition{false};
 
-  FargoVelocityFunc fargoVelocityFunc{NULL};  // The user-defined fargo velocity function
+  FargoVelocityFunc fargoVelocityFunc;  // The user-defined fargo velocity function
 };
 
 

src/fluid/RiemannSolver/CMakeLists.txt

@@ -22,6 +22,7 @@ target_sources(idefix
   PUBLIC ${CMAKE_CURRENT_LIST_DIR}/extrapolateToFaces.hpp
   PUBLIC ${CMAKE_CURRENT_LIST_DIR}/flux.hpp
   PUBLIC ${CMAKE_CURRENT_LIST_DIR}/riemannSolver.hpp
+  PUBLIC ${CMAKE_CURRENT_LIST_DIR}/shockFlag.hpp
   PUBLIC ${CMAKE_CURRENT_LIST_DIR}/shockFlattening.hpp
   PUBLIC ${CMAKE_CURRENT_LIST_DIR}/slopeLimiter.hpp
   )

src/fluid/RiemannSolver/HDsolvers/hllHD.hpp

@@ -36,7 +36,7 @@ void RiemannSolver<Phys>::HllHD(IdefixArray4D<real> &Flux) {
   IdefixArray4D<real> Vc = this->Vc;
   IdefixArray4D<real> Vs = this->Vs;
   IdefixArray3D<real> cMax = this->cMax;
-  EquationOfState eos = *(hydro->eos.get());
+  EquationOfStateDevice eos = *(hydro->eos.get());
 
   // Required for high order interpolations
   IdefixArray1D<real> dx = this->data->dx[DIR];

src/fluid/RiemannSolver/HDsolvers/hllcHD.hpp

@@ -37,7 +37,7 @@ void RiemannSolver<Phys>::HllcHD(IdefixArray4D<real> &Flux) {
   IdefixArray4D<real> Vs = this->Vs;
   IdefixArray3D<real> cMax = this->cMax;
 
-  EquationOfState eos = *(hydro->eos.get());
+  EquationOfStateDevice eos = *(hydro->eos.get());
 
   ExtrapolateToFaces<Phys,DIR> extrapol = *this->GetExtrapolator<DIR>();
 

src/fluid/RiemannSolver/HDsolvers/roeHD.hpp

@@ -58,7 +58,7 @@ void RiemannSolver<Phys>::RoeHD(IdefixArray4D<real> &Flux) {
   IdefixArray4D<real> Vs = this->Vs;
   IdefixArray3D<real> cMax = this->cMax;
 
-  EquationOfState eos = *(hydro->eos.get());
+  EquationOfStateDevice eos = *(hydro->eos.get());
 
   ExtrapolateToFaces<Phys,DIR> extrapol = *this->GetExtrapolator<DIR>();
 

src/fluid/RiemannSolver/HDsolvers/tvdlfHD.hpp

@@ -36,7 +36,7 @@ void RiemannSolver<Phys>::TvdlfHD(IdefixArray4D<real> &Flux) {
   IdefixArray4D<real> Vc = this->Vc;
   IdefixArray4D<real> Vs = this->Vs;
   IdefixArray3D<real> cMax = this->cMax;
-  EquationOfState eos = *(hydro->eos.get());
+  EquationOfStateDevice eos = *(hydro->eos.get());
 
   // Required for high order interpolations
   IdefixArray1D<real> dx = this->data->dx[DIR];

src/fluid/RiemannSolver/MHDsolvers/hllMHD.hpp

@@ -83,7 +83,7 @@ void RiemannSolver<Phys>::HllMHD(IdefixArray4D<real> &Flux) {
   IdefixArray3D<real> dL;
   IdefixArray3D<real> dR;
 
-  EquationOfState eos = *(hydro->eos.get());
+  EquationOfStateDevice eos = *(hydro->eos.get());
 
   [[maybe_unused]] real xHConstant = hydro->xH;
 

src/fluid/RiemannSolver/MHDsolvers/hlldMHD.hpp

@@ -80,7 +80,7 @@ void RiemannSolver<Phys>::HlldMHD(IdefixArray4D<real> &Flux) {
   IdefixArray3D<real> dL;
   IdefixArray3D<real> dR;
 
-  EquationOfState eos = *(hydro->eos.get());
+  EquationOfStateDevice eos = *(hydro->eos.get());
 
   ExtrapolateToFaces<Phys,DIR> extrapol = *this->GetExtrapolator<DIR>();
 

src/fluid/RiemannSolver/MHDsolvers/roeMHD.hpp

@@ -104,7 +104,7 @@ void RiemannSolver<Phys>::RoeMHD(IdefixArray4D<real> &Flux) {
   IdefixArray3D<real> dL;
   IdefixArray3D<real> dR;
 
-  EquationOfState eos = *(hydro->eos.get());
+  EquationOfStateDevice eos = *(hydro->eos.get());
 
   // TODO(baghdads) what is this delta?
   real delta    = 1.e-6;

src/fluid/RiemannSolver/MHDsolvers/tvdlfMHD.hpp

@@ -71,7 +71,7 @@ void RiemannSolver<Phys>::TvdlfMHD(IdefixArray4D<real> &Flux) {
   IdefixArray3D<real> dL;
   IdefixArray3D<real> dR;
 
-  EquationOfState eos = *(hydro->eos.get());
+  EquationOfStateDevice eos = *(hydro->eos.get());
 
   ExtrapolateToFaces<Phys,DIR> extrapol = *this->GetExtrapolator<DIR>();
   // Define normal, tangent and bi-tanget indices