diff --git a/src/dataBlock/coarsen.cpp b/src/dataBlock/coarsen.cpp index f4b75b9d..be0ccddb 100644 --- a/src/dataBlock/coarsen.cpp +++ b/src/dataBlock/coarsen.cpp @@ -41,7 +41,7 @@ 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"); } @@ -49,7 +49,7 @@ void DataBlock::ComputeGridCoarseningLevels() { // 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(); diff --git a/src/dataBlock/dataBlock.cpp b/src/dataBlock/dataBlock.cpp index 3440faca..463f148c 100644 --- a/src/dataBlock/dataBlock.cpp +++ b/src/dataBlock/dataBlock.cpp @@ -386,10 +386,7 @@ 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(); } } @@ -397,20 +394,7 @@ void DataBlock::LaunchUserStepLast() { 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; -} diff --git a/src/dataBlock/dataBlock.hpp b/src/dataBlock/dataBlock.hpp index 6d322e7e..3acc4432 100644 --- a/src/dataBlock/dataBlock.hpp +++ b/src/dataBlock/dataBlock.hpp @@ -17,11 +17,12 @@ #ifndef DATABLOCK_DATABLOCK_HPP_ #define DATABLOCK_DATABLOCK_HPP_ -#include -#include +#include +#include #include #include -#include +#include +#include #include "idefix.hpp" #include "grid.hpp" @@ -50,15 +51,13 @@ class PlanetarySystem; template class Fluid; class SubGrid; -class Vtk; -class Dump; // Forward class hydro declaration #include "physics.hpp" -using GridCoarseningFunc = void(*) (DataBlock &); +using GridCoarseningFunc = std::function; -using StepFunc = void (*) (DataBlock &, const real t, const real dt); +using StepFunc = std::function; 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 + void EnrollUserStepFirst(Func func) { + haveUserStepFirst = true; + userStepFirst = func; + } + template + 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_ diff --git a/src/dataBlock/fargo.cpp b/src/dataBlock/fargo.cpp index 046a11d5..682b5340 100644 --- a/src/dataBlock/fargo.cpp +++ b/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); diff --git a/src/dataBlock/fargo.hpp b/src/dataBlock/fargo.hpp index 9347f716..fb62a220 100644 --- a/src/dataBlock/fargo.hpp +++ b/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 }; diff --git a/src/fluid/RiemannSolver/CMakeLists.txt b/src/fluid/RiemannSolver/CMakeLists.txt index 9099fea2..646c4bdd 100644 --- a/src/fluid/RiemannSolver/CMakeLists.txt +++ b/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 ) diff --git a/src/fluid/RiemannSolver/HDsolvers/hllHD.hpp b/src/fluid/RiemannSolver/HDsolvers/hllHD.hpp index 276f4e4c..cbd2322f 100644 --- a/src/fluid/RiemannSolver/HDsolvers/hllHD.hpp +++ b/src/fluid/RiemannSolver/HDsolvers/hllHD.hpp @@ -36,7 +36,7 @@ void RiemannSolver::HllHD(IdefixArray4D &Flux) { IdefixArray4D Vc = this->Vc; IdefixArray4D Vs = this->Vs; IdefixArray3D cMax = this->cMax; - EquationOfState eos = *(hydro->eos.get()); + EquationOfStateDevice eos = *(hydro->eos.get()); // Required for high order interpolations IdefixArray1D dx = this->data->dx[DIR]; diff --git a/src/fluid/RiemannSolver/HDsolvers/hllcHD.hpp b/src/fluid/RiemannSolver/HDsolvers/hllcHD.hpp index 2585fd7a..a25b3832 100644 --- a/src/fluid/RiemannSolver/HDsolvers/hllcHD.hpp +++ b/src/fluid/RiemannSolver/HDsolvers/hllcHD.hpp @@ -37,7 +37,7 @@ void RiemannSolver::HllcHD(IdefixArray4D &Flux) { IdefixArray4D Vs = this->Vs; IdefixArray3D cMax = this->cMax; - EquationOfState eos = *(hydro->eos.get()); + EquationOfStateDevice eos = *(hydro->eos.get()); ExtrapolateToFaces extrapol = *this->GetExtrapolator(); diff --git a/src/fluid/RiemannSolver/HDsolvers/roeHD.hpp b/src/fluid/RiemannSolver/HDsolvers/roeHD.hpp index b985b40a..3b770f13 100644 --- a/src/fluid/RiemannSolver/HDsolvers/roeHD.hpp +++ b/src/fluid/RiemannSolver/HDsolvers/roeHD.hpp @@ -58,7 +58,7 @@ void RiemannSolver::RoeHD(IdefixArray4D &Flux) { IdefixArray4D Vs = this->Vs; IdefixArray3D cMax = this->cMax; - EquationOfState eos = *(hydro->eos.get()); + EquationOfStateDevice eos = *(hydro->eos.get()); ExtrapolateToFaces extrapol = *this->GetExtrapolator(); diff --git a/src/fluid/RiemannSolver/HDsolvers/tvdlfHD.hpp b/src/fluid/RiemannSolver/HDsolvers/tvdlfHD.hpp index 7a389dd1..4fd030d8 100644 --- a/src/fluid/RiemannSolver/HDsolvers/tvdlfHD.hpp +++ b/src/fluid/RiemannSolver/HDsolvers/tvdlfHD.hpp @@ -36,7 +36,7 @@ void RiemannSolver::TvdlfHD(IdefixArray4D &Flux) { IdefixArray4D Vc = this->Vc; IdefixArray4D Vs = this->Vs; IdefixArray3D cMax = this->cMax; - EquationOfState eos = *(hydro->eos.get()); + EquationOfStateDevice eos = *(hydro->eos.get()); // Required for high order interpolations IdefixArray1D dx = this->data->dx[DIR]; diff --git a/src/fluid/RiemannSolver/MHDsolvers/hllMHD.hpp b/src/fluid/RiemannSolver/MHDsolvers/hllMHD.hpp index ae4c4b6c..1a4f4fc4 100644 --- a/src/fluid/RiemannSolver/MHDsolvers/hllMHD.hpp +++ b/src/fluid/RiemannSolver/MHDsolvers/hllMHD.hpp @@ -83,7 +83,7 @@ void RiemannSolver::HllMHD(IdefixArray4D &Flux) { IdefixArray3D dL; IdefixArray3D dR; - EquationOfState eos = *(hydro->eos.get()); + EquationOfStateDevice eos = *(hydro->eos.get()); [[maybe_unused]] real xHConstant = hydro->xH; diff --git a/src/fluid/RiemannSolver/MHDsolvers/hlldMHD.hpp b/src/fluid/RiemannSolver/MHDsolvers/hlldMHD.hpp index 7d0237a6..c95f6d5c 100644 --- a/src/fluid/RiemannSolver/MHDsolvers/hlldMHD.hpp +++ b/src/fluid/RiemannSolver/MHDsolvers/hlldMHD.hpp @@ -80,7 +80,7 @@ void RiemannSolver::HlldMHD(IdefixArray4D &Flux) { IdefixArray3D dL; IdefixArray3D dR; - EquationOfState eos = *(hydro->eos.get()); + EquationOfStateDevice eos = *(hydro->eos.get()); ExtrapolateToFaces extrapol = *this->GetExtrapolator(); diff --git a/src/fluid/RiemannSolver/MHDsolvers/roeMHD.hpp b/src/fluid/RiemannSolver/MHDsolvers/roeMHD.hpp index 4395c202..607cd9c3 100644 --- a/src/fluid/RiemannSolver/MHDsolvers/roeMHD.hpp +++ b/src/fluid/RiemannSolver/MHDsolvers/roeMHD.hpp @@ -104,7 +104,7 @@ void RiemannSolver::RoeMHD(IdefixArray4D &Flux) { IdefixArray3D dL; IdefixArray3D dR; - EquationOfState eos = *(hydro->eos.get()); + EquationOfStateDevice eos = *(hydro->eos.get()); // TODO(baghdads) what is this delta? real delta = 1.e-6; diff --git a/src/fluid/RiemannSolver/MHDsolvers/tvdlfMHD.hpp b/src/fluid/RiemannSolver/MHDsolvers/tvdlfMHD.hpp index 743c6ada..3e4fc849 100644 --- a/src/fluid/RiemannSolver/MHDsolvers/tvdlfMHD.hpp +++ b/src/fluid/RiemannSolver/MHDsolvers/tvdlfMHD.hpp @@ -71,7 +71,7 @@ void RiemannSolver::TvdlfMHD(IdefixArray4D &Flux) { IdefixArray3D dL; IdefixArray3D dR; - EquationOfState eos = *(hydro->eos.get()); + EquationOfStateDevice eos = *(hydro->eos.get()); ExtrapolateToFaces extrapol = *this->GetExtrapolator(); // Define normal, tangent and bi-tanget indices diff --git a/src/fluid/RiemannSolver/extrapolateToFaces.hpp b/src/fluid/RiemannSolver/extrapolateToFaces.hpp index c832bad9..60911454 100644 --- a/src/fluid/RiemannSolver/extrapolateToFaces.hpp +++ b/src/fluid/RiemannSolver/extrapolateToFaces.hpp @@ -3,7 +3,7 @@ // // Source file src/fluid/RiemannSolver/extrapolateToFaces.hpp // -// Last modified : 10/2023 +// Last modified : 10/2026 // // Copyright(C) by : // - Geoffroy Lesur (IPAG/UGA/CNRS - 2023) @@ -14,9 +14,8 @@ #ifndef FLUID_RIEMANNSOLVER_EXTRAPOLATETOFACES_HPP_ #define FLUID_RIEMANNSOLVER_EXTRAPOLATETOFACES_HPP_ -#include "fluid.hpp" #include "dataBlock.hpp" -#include "shockFlattening.hpp" +#include "shockFlag.hpp" #include "slopeLimiter.hpp" // Build a left and right extrapolation of the primitive variables along direction dir @@ -33,302 +32,259 @@ template; + // PLM weights for irregular grids (trivially copyable, hence device-safe) + struct GridWeights { + IdefixArray1D cp; ///< (xgc(i+1)-xgc(i)) / (xr(i)-xgc(i)) + IdefixArray1D cm; ///< (xgc(i)-xgc(i-1)) / (xgc(i)-xr(i-1)) + IdefixArray1D dp; ///< (xr(i)-xgc(i)) / dx(i) + IdefixArray1D dm; ///< (xgc(i)-xr(i-1)) / dx(i) + IdefixArray1D wp; ///< dx(i) / (xgc(i+1)-xgc(i)) + IdefixArray1D wm; ///< dx(i) / (xgc(i)-xgc(i-1)) + }; + + // Index offsets along the reconstruction direction + static constexpr int ioffset = (dir == IDIR ? 1 : 0); + static constexpr int joffset = (dir == JDIR ? 1 : 0); + static constexpr int koffset = (dir == KDIR ? 1 : 0); + public: - explicit ExtrapolateToFaces(RiemannSolver *rSolver): - Vc(rSolver->hydro->Vc), - dx(rSolver->hydro->data->dx[dir]), - shockFlattening(rSolver->haveShockFlattening), - isRegularGrid(rSolver->hydro->data->mygrid->isRegularCartesian) { - if(shockFlattening) { - flags = rSolver->shockFlattening->flagArray; - } - if(!isRegularGrid) { - ComputePLMweights(rSolver->hydro->data); - } + // data: datablock holding the grid information + // VcArray: primitive variables in cell centers + // shockFlags: shock flag array (only used when hasShockFlattening is true) + explicit ExtrapolateToFaces(DataBlock *data, + IdefixArray4D VcArray, + IdefixArray3D shockFlags = {}, + bool hasShockFlattening = false): + Vc(VcArray), + dx(data->dx[dir]), + flags(shockFlags), + isRegularGrid(data->mygrid->isRegularCartesian), + shockFlattening(hasShockFlattening) { + if(!isRegularGrid) { + ComputeGridWeights(data); + } + } + + // Extrapolate the primitive variables to the left (vL) and right (vR) states + // of the interface located between cell i-1 and cell i + KOKKOS_FORCEINLINE_FUNCTION void ExtrapolatePrimVar(const int i, + const int j, + const int k, + real vL[], real vR[]) const { + if constexpr(order == 1) { + ExtrapolateOrder1(i,j,k,vL,vR); + } else if constexpr(order == 2) { + if(isRegularGrid) { + ExtrapolateOrder2Regular(i,j,k,vL,vR); + } else { + ExtrapolateOrder2Irregular(i,j,k,vL,vR); + } + } else if constexpr(order == 3) { + ExtrapolateOrder3(i,j,k,vL,vR); + } else if constexpr(order == 4) { + ExtrapolateOrder4(i,j,k,vL,vR); + } } - void ComputePLMweights(DataBlock *data) { - // Allocate weight arrays - cpArray = IdefixArray1D("ExtrapolateToFaces_cp",data->np_tot[dir]); - cmArray = IdefixArray1D("ExtrapolateToFaces_cm",data->np_tot[dir]); - dpArray = IdefixArray1D("ExtrapolateToFaces_dp",data->np_tot[dir]); - dmArray = IdefixArray1D("ExtrapolateToFaces_dm",data->np_tot[dir]); - wpArray = IdefixArray1D("ExtrapolateToFaces_wp",data->np_tot[dir]); - wmArray = IdefixArray1D("ExtrapolateToFaces_wm",data->np_tot[dir]); + // Allocate and compute the PLM weights along direction dir (irregular grids only) + void ComputeGridWeights(DataBlock *data) { + gridWeights.cp = IdefixArray1D("ExtrapolateToFaces_cp",data->np_tot[dir]); + gridWeights.cm = IdefixArray1D("ExtrapolateToFaces_cm",data->np_tot[dir]); + gridWeights.dp = IdefixArray1D("ExtrapolateToFaces_dp",data->np_tot[dir]); + gridWeights.dm = IdefixArray1D("ExtrapolateToFaces_dm",data->np_tot[dir]); + gridWeights.wp = IdefixArray1D("ExtrapolateToFaces_wp",data->np_tot[dir]); + gridWeights.wm = IdefixArray1D("ExtrapolateToFaces_wm",data->np_tot[dir]); - auto dx = data->dx[dir]; + auto dxGrid = data->dx[dir]; auto xgc = data->xgc[dir]; auto xr = data->xr[dir]; - auto wp = wpArray; - auto wm = wmArray; - auto cp = cpArray; - auto cm = cmArray; - auto dp = dpArray; - auto dm = dmArray; - - idefix_for("ComputePLMweights",1,data->np_tot[dir]-1, - KOKKOS_LAMBDA(const int i) { - wp(i) = dx(i) / (xgc(i+1) - xgc(i)); - wm(i) = dx(i) / (xgc(i) - xgc(i-1)); - cp(i) = (xgc(i+1) - xgc(i)) / (xr(i) - xgc(i)); - cm(i) = (xgc(i) - xgc(i-1)) / (xgc(i) - xr(i-1)); - dp(i) = (xr(i) - xgc(i)) / dx(i); - dm(i) = (xgc(i) - xr(i-1)) / dx(i); - }); + auto wp = gridWeights.wp; + auto wm = gridWeights.wm; + auto cp = gridWeights.cp; + auto cm = gridWeights.cm; + auto dp = gridWeights.dp; + auto dm = gridWeights.dm; + + idefix_for("ComputeGridWeights",1,data->np_tot[dir]-1, + KOKKOS_LAMBDA(const int i) { + wp(i) = dxGrid(i) / (xgc(i+1) - xgc(i)); + wm(i) = dxGrid(i) / (xgc(i) - xgc(i-1)); + cp(i) = (xgc(i+1) - xgc(i)) / (xr(i) - xgc(i)); + cm(i) = (xgc(i) - xgc(i-1)) / (xgc(i) - xr(i-1)); + dp(i) = (xr(i) - xgc(i)) / dxGrid(i); + dm(i) = (xgc(i) - xr(i-1)) / dxGrid(i); + }); } + // Shock flag at cell (k,j,i), or FlagShock::None when shock flattening is disabled + KOKKOS_FORCEINLINE_FUNCTION FlagShock GetFlag(const int k, const int j, const int i) const { + return shockFlattening ? flags(k,j,i) : FlagShock::None; + } + // PLM slope with a minmod fallback in shocked cells (regular grid) + KOKKOS_FORCEINLINE_FUNCTION real LimitSlope(const real dvp, const real dvm, + const FlagShock flag) const { + if(shockFlattening && flag == FlagShock::Shock) { + return SL::MinModLim(dvp,dvm); + } + return SL::PLMLim(dvp,dvm); + } - KOKKOS_FORCEINLINE_FUNCTION void ExtrapolatePrimVar(const int i, - const int j, - const int k, - real vL[], real vR[]) const { - // 1-- Store the primitive variables on the left, right, and averaged states - constexpr int ioffset = (dir==IDIR ? 1 : 0); - constexpr int joffset = (dir==JDIR ? 1 : 0); - constexpr int koffset = (dir==KDIR ? 1 : 0); - - for(int nv = 0 ; nv < Phys::nvar ; nv++) { - if constexpr(order == 1) { - vL[nv] = Vc(nv,k-koffset,j-joffset,i-ioffset); - vR[nv] = Vc(nv,k,j,i); - } else if constexpr(order == 2) { - if(isRegularGrid) { - ///////////////////////////////////// - // Regular Grid, PLM reconstruction - ///////////////////////////////////// - real dvm = Vc(nv,k-koffset,j-joffset,i-ioffset) - -Vc(nv,k-2*koffset,j-2*joffset,i-2*ioffset); - real dvp = Vc(nv,k,j,i)-Vc(nv,k-koffset,j-joffset,i-ioffset); - - real dv; - if(shockFlattening) { - if(flags(k-koffset,j-joffset,i-ioffset) == FlagShock::Shock) { - // Force slope limiter to minmod - dv = SL::MinModLim(dvp,dvm); - } else { - dv = SL::PLMLim(dvp,dvm); - } - } else { // No shock flattening - dv = SL::PLMLim(dvp,dvm); - } - - vL[nv] = Vc(nv,k-koffset,j-joffset,i-ioffset) + HALF_F*dv; - - dvm = dvp; - dvp = Vc(nv,k+koffset,j+joffset,i+ioffset) - Vc(nv,k,j,i); - - if(shockFlattening) { - if(flags(k,j,i) == FlagShock::Shock) { - dv = SL::MinModLim(dvp,dvm); - } else { - dv = SL::PLMLim(dvp,dvm); - } - } else { // No shock flattening - dv = SL::PLMLim(dvp,dvm); - } - - vR[nv] = Vc(nv,k,j,i) - HALF_F*dv; - } else { - ///////////////////////////////////// - // Irregular Grid, PLM reconstruction - ///////////////////////////////////// - const int index = ioffset*i + joffset*j + koffset*k; - - real dvm = Vc(nv,k-koffset,j-joffset,i-ioffset) - -Vc(nv,k-2*koffset,j-2*joffset,i-2*ioffset); - real dvp = Vc(nv,k,j,i)-Vc(nv,k-koffset,j-joffset,i-ioffset); - - dvm *= wmArray(index-1); - dvp *= wpArray(index-1); - real cp = cpArray(index-1); - real cm = cmArray(index-1); - - real dv; - if(shockFlattening) { - if(flags(k-koffset,j-joffset,i-ioffset) == FlagShock::Shock) { - // Force slope limiter to minmod - dv = SL::MinModLim(dvp,dvm); - } else { - dv = SL::PLMLim(dvp,dvm,cp,cm); - } - } else { // No shock flattening - dv = SL::PLMLim(dvp,dvm,cp,cm); - } - - vL[nv] = Vc(nv,k-koffset,j-joffset,i-ioffset) + dpArray(index-1)*dv; - - dvm = Vc(nv,k,j,i)-Vc(nv,k-koffset,j-joffset,i-ioffset); - dvp = Vc(nv,k+koffset,j+joffset,i+ioffset) - Vc(nv,k,j,i); - dvm *= wmArray(index); - dvp *= wpArray(index); - cp = cpArray(index); - cm = cmArray(index); - - if(shockFlattening) { - if(flags(k,j,i) == FlagShock::Shock) { - dv = SL::MinModLim(dvp,dvm); - } else { - dv = SL::PLMLim(dvp,dvm,cp,cm); - } - } else { // No shock flattening - dv = SL::PLMLim(dvp,dvm,cp,cm); - } - vR[nv] = Vc(nv,k,j,i) - dmArray(index)*dv; - } // Regular grid - - } else if constexpr(order == 3) { - // 1D index along the chosen direction - const int index = ioffset*i + joffset*j + koffset*k; - real dvm = Vc(nv,k-koffset,j-joffset,i-ioffset) - -Vc(nv,k-2*koffset,j-2*joffset,i-2*ioffset); - real dvp = Vc(nv,k,j,i)-Vc(nv,k-koffset,j-joffset,i-ioffset); - - // Limo3 limiter - real dv; - if(shockFlattening) { - if(flags(k-koffset,j-joffset,i-ioffset) == FlagShock::Shock) { - // Force slope limiter to minmod - dv = SL::MinModLim(dvp,dvm); - } else { - dv = dvp * SL::LimO3Lim(dvp, dvm, dx(index-1)); - } - } else { // No shock flattening - dv = dvp * SL::LimO3Lim(dvp, dvm, dx(index-1)); - } - - vL[nv] = Vc(nv,k-koffset,j-joffset,i-ioffset) + HALF_F*dv; - - // Check positivity - if(nv==RHO) { - // If face element is negative, revert to minmod - if(vL[nv] <= 0.0) { - dv = SL::MinModLim(dvp,dvm); - vL[nv] = Vc(nv,k-koffset,j-joffset,i-ioffset) + HALF_F*dv; - } - } - if constexpr(Phys::pressure) { - if(nv==PRS) { - // If face element is negative, revert to minmod - if(vL[nv] <= 0.0) { - dv = SL::MinModLim(dvp,dvm); - vL[nv] = Vc(nv,k-koffset,j-joffset,i-ioffset) + HALF_F*dv; - } - } - } - - dvm = dvp; - dvp = Vc(nv,k+koffset,j+joffset,i+ioffset) - Vc(nv,k,j,i); - - // Limo3 limiter - if(shockFlattening) { - if(flags(k,j,i) == FlagShock::Shock) { - // Force slope limiter to minmod - dv = SL::MinModLim(dvp,dvm); - } else { - dv = dvm * SL::LimO3Lim(dvm, dvp, dx(index)); - } - } else { // No shock flattening - dv = dvm * SL::LimO3Lim(dvm, dvp, dx(index)); - } - - vR[nv] = Vc(nv,k,j,i) - HALF_F*dv; + // PLM slope with a minmod fallback in shocked cells (irregular grid) + KOKKOS_FORCEINLINE_FUNCTION real LimitSlope(const real dvp, const real dvm, + const FlagShock flag, + const real cp, const real cm) const { + if(shockFlattening && flag == FlagShock::Shock) { + return SL::MinModLim(dvp,dvm); + } + return SL::PLMLim(dvp,dvm,cp,cm); + } - // Check positivity - if(nv==RHO) { - // If face element is negative, revert to vanleer - if(vR[nv] <= 0.0) { - dv = SL::MinModLim(dvp,dvm); - vR[nv] = Vc(nv,k,j,i) - HALF_F*dv; - } - } - if constexpr(Phys::pressure) { - if(nv==PRS) { - // If face element is negative, revert to vanleer - if(vR[nv] <= 0.0) { - dv = SL::MinModLim(dvp,dvm); - vR[nv] = Vc(nv,k,j,i) - HALF_F*dv; - } - } - } - } else if constexpr(order == 4) { - // Reconstruction in cell i-1 - real vm2 = Vc(nv,k-3*koffset,j-3*joffset,i-3*ioffset);; - real vm1 = Vc(nv,k-2*koffset,j-2*joffset,i-2*ioffset); - real v0 = Vc(nv,k-koffset,j-joffset,i-ioffset); - real vp1 = Vc(nv,k,j,i); - real vp2 = Vc(nv,k+koffset,j+joffset,i+ioffset); + // LimO3 slope with a minmod fallback in shocked cells. + // dvU multiplies the limiter, (dvA,dvB) are the limiter arguments. + KOKKOS_FORCEINLINE_FUNCTION real LimitSlopeO3(const real dvU, const real dvA, const real dvB, + const FlagShock flag, const real dxm) const { + if(shockFlattening && flag == FlagShock::Shock) { + return SL::MinModLim(dvA,dvB); + } + return dvU * SL::LimO3Lim(dvA,dvB,dxm); + } - real vr,vl; - SL::getPPMStates(vm2, vm1, v0, vp1, vp2, vl, vr); - // vL= left side of current interface (i-1/2)= right side of cell i-1 + // Enforce positivity of the density and pressure face values by reverting to the + // fallback slope dvFallback (sign = +1 for a left state, -1 for a right state) + KOKKOS_FORCEINLINE_FUNCTION void EnforcePositivity(const int nv, real &vFace, + const real v0, const real dvFallback, + const real sign) const { + bool enforce = (nv == RHO); + if constexpr(Phys::pressure) { + enforce = enforce || (nv == PRS); + } + if(enforce && vFace <= 0.0) { + vFace = v0 + sign*HALF_F*dvFallback; + } + } - // Check positivity - if(nv==RHO) { - // If face element is negative, revert to vanleer - if(vr <= 0.0) { - real dv = SL::PLMLim(vp1-v0,v0-vm1); - vr = v0+HALF_F*dv; - } - } - if constexpr(Phys::pressure) { - if(nv==PRS) { - // If face element is negative, revert to vanleer - if(vr <= 0.0) { - real dv = SL::PLMLim(vp1-v0,v0-vm1); - vr = v0+HALF_F*dv; - } - } - } + // First order: piecewise constant reconstruction + KOKKOS_FORCEINLINE_FUNCTION void ExtrapolateOrder1(const int i, const int j, const int k, + real vL[], real vR[]) const { + for(int nv = 0 ; nv < Phys::nvar ; nv++) { + vL[nv] = Vc(nv,k-koffset,j-joffset,i-ioffset); + vR[nv] = Vc(nv,k,j,i); + } + } - vL[nv] = vr; - // Reconstruction in cell i + // Second order (PLM) reconstruction on a regular grid + KOKKOS_FORCEINLINE_FUNCTION void ExtrapolateOrder2Regular(const int i, const int j, const int k, + real vL[], real vR[]) const { + for(int nv = 0 ; nv < Phys::nvar ; nv++) { + // Left state: reconstruction in cell i-1 + real dvm = Vc(nv,k-koffset,j-joffset,i-ioffset) + -Vc(nv,k-2*koffset,j-2*joffset,i-2*ioffset); + real dvp = Vc(nv,k,j,i)-Vc(nv,k-koffset,j-joffset,i-ioffset); + real dv = LimitSlope(dvp,dvm,GetFlag(k-koffset,j-joffset,i-ioffset)); + vL[nv] = Vc(nv,k-koffset,j-joffset,i-ioffset) + HALF_F*dv; + + // Right state: reconstruction in cell i + dvm = dvp; + dvp = Vc(nv,k+koffset,j+joffset,i+ioffset) - Vc(nv,k,j,i); + dv = LimitSlope(dvp,dvm,GetFlag(k,j,i)); + vR[nv] = Vc(nv,k,j,i) - HALF_F*dv; + } + } - vm2 = vm1; - vm1 = v0; - v0 = vp1; - vp1 = vp2; - vp2 = Vc(nv,k+2*koffset,j+2*joffset,i+2*ioffset); + // Second order (PLM) reconstruction on an irregular grid + KOKKOS_FORCEINLINE_FUNCTION void ExtrapolateOrder2Irregular(const int i, const int j, const int k, + real vL[], real vR[]) const { + // 1D index along the chosen direction + const int index = ioffset*i + joffset*j + koffset*k; - SL::getPPMStates(vm2, vm1, v0, vp1, vp2, vl, vr); + for(int nv = 0 ; nv < Phys::nvar ; nv++) { + // Left state: reconstruction in cell i-1 + real dvm = Vc(nv,k-koffset,j-joffset,i-ioffset) + -Vc(nv,k-2*koffset,j-2*joffset,i-2*ioffset); + real dvp = Vc(nv,k,j,i)-Vc(nv,k-koffset,j-joffset,i-ioffset); + dvm *= gridWeights.wm(index-1); + dvp *= gridWeights.wp(index-1); + real dv = LimitSlope(dvp,dvm,GetFlag(k-koffset,j-joffset,i-ioffset), + gridWeights.cp(index-1),gridWeights.cm(index-1)); + vL[nv] = Vc(nv,k-koffset,j-joffset,i-ioffset) + gridWeights.dp(index-1)*dv; + + // Right state: reconstruction in cell i + dvm = Vc(nv,k,j,i)-Vc(nv,k-koffset,j-joffset,i-ioffset); + dvp = Vc(nv,k+koffset,j+joffset,i+ioffset) - Vc(nv,k,j,i); + dvm *= gridWeights.wm(index); + dvp *= gridWeights.wp(index); + dv = LimitSlope(dvp,dvm,GetFlag(k,j,i), + gridWeights.cp(index),gridWeights.cm(index)); + vR[nv] = Vc(nv,k,j,i) - gridWeights.dm(index)*dv; + } + } - // Check positivity - if(nv==RHO) { - // If face element is negative, revert to vanleer - if(vl <= 0.0) { - real dv = SL::PLMLim(vp1-v0,v0-vm1); - vl = v0-HALF_F*dv; - } - } - if constexpr(Phys::pressure) { - if(nv==PRS) { - // If face element is negative, revert to vanleer - if(vl <= 0.0) { - real dv = SL::PLMLim(vp1-v0,v0-vm1); - vl = v0-HALF_F*dv; - } - } - } + // Third order (LimO3) reconstruction + KOKKOS_FORCEINLINE_FUNCTION void ExtrapolateOrder3(const int i, const int j, const int k, + real vL[], real vR[]) const { + // 1D index along the chosen direction + const int index = ioffset*i + joffset*j + koffset*k; - vR[nv] = vl; - } + for(int nv = 0 ; nv < Phys::nvar ; nv++) { + // Left state: reconstruction in cell i-1 + real dvm = Vc(nv,k-koffset,j-joffset,i-ioffset) + -Vc(nv,k-2*koffset,j-2*joffset,i-2*ioffset); + real dvp = Vc(nv,k,j,i)-Vc(nv,k-koffset,j-joffset,i-ioffset); + real dv = LimitSlopeO3(dvp,dvp,dvm,GetFlag(k-koffset,j-joffset,i-ioffset),dx(index-1)); + vL[nv] = Vc(nv,k-koffset,j-joffset,i-ioffset) + HALF_F*dv; + // Check positivity of density and pressure + EnforcePositivity(nv,vL[nv],Vc(nv,k-koffset,j-joffset,i-ioffset), + SL::MinModLim(dvp,dvm),ONE_F); + + // Right state: reconstruction in cell i + dvm = dvp; + dvp = Vc(nv,k+koffset,j+joffset,i+ioffset) - Vc(nv,k,j,i); + dv = LimitSlopeO3(dvm,dvm,dvp,GetFlag(k,j,i),dx(index)); + vR[nv] = Vc(nv,k,j,i) - HALF_F*dv; + // Check positivity of density and pressure + EnforcePositivity(nv,vR[nv],Vc(nv,k,j,i),SL::MinModLim(dvp,dvm),-ONE_F); } } - IdefixArray4D Vc; - IdefixArray1D dx; - IdefixArray3D flags; + // Fourth order (PPM) reconstruction + KOKKOS_FORCEINLINE_FUNCTION void ExtrapolateOrder4(const int i, const int j, const int k, + real vL[], real vR[]) const { + for(int nv = 0 ; nv < Phys::nvar ; nv++) { + // Reconstruction in cell i-1 + real vm2 = Vc(nv,k-3*koffset,j-3*joffset,i-3*ioffset); + real vm1 = Vc(nv,k-2*koffset,j-2*joffset,i-2*ioffset); + real v0 = Vc(nv,k-koffset,j-joffset,i-ioffset); + real vp1 = Vc(nv,k,j,i); + real vp2 = Vc(nv,k+koffset,j+joffset,i+ioffset); + + real vr, vl; + SL::getPPMStates(vm2,vm1,v0,vp1,vp2,vl,vr); + // vL is the left side of the current interface, i.e. the right side of cell i-1 + EnforcePositivity(nv,vr,v0,SL::PLMLim(vp1-v0,v0-vm1),ONE_F); + vL[nv] = vr; + + // Reconstruction in cell i + vm2 = vm1; + vm1 = v0; + v0 = vp1; + vp1 = vp2; + vp2 = Vc(nv,k+2*koffset,j+2*joffset,i+2*ioffset); + + SL::getPPMStates(vm2,vm1,v0,vp1,vp2,vl,vr); + // vR is the right side of the current interface, i.e. the left side of cell i + EnforcePositivity(nv,vl,v0,SL::PLMLim(vp1-v0,v0-vm1),-ONE_F); + vR[nv] = vl; + } + } - IdefixArray1D cpArray; - IdefixArray1D cmArray; - IdefixArray1D dpArray; - IdefixArray1D dmArray; - IdefixArray1D wpArray; - IdefixArray1D wmArray; + IdefixArray4D Vc; ///< primitive variables in cell centers + IdefixArray1D dx; ///< cell width along direction dir + IdefixArray3D flags; ///< shock flags (empty when shock flattening is disabled) + GridWeights gridWeights; ///< PLM weights (irregular grids only) bool isRegularGrid{true}; bool shockFlattening{false}; }; - #endif // FLUID_RIEMANNSOLVER_EXTRAPOLATETOFACES_HPP_ diff --git a/src/fluid/RiemannSolver/riemannSolver.hpp b/src/fluid/RiemannSolver/riemannSolver.hpp index bffaaa43..f8fffc14 100644 --- a/src/fluid/RiemannSolver/riemannSolver.hpp +++ b/src/fluid/RiemannSolver/riemannSolver.hpp @@ -73,9 +73,6 @@ class RiemannSolver { std::unique_ptr> shockFlattening; private: - template - friend class ExtrapolateToFaces; - IdefixArray4D Vc; IdefixArray4D Vs; IdefixArray3D cMax; @@ -163,12 +160,19 @@ RiemannSolver::RiemannSolver(Input &input, Fluid* hydro) : Vc{hydro- } // init slope limiters - slopeLimIDIR = std::make_unique>(this); + IdefixArray3D shockFlags{}; + if(haveShockFlattening) { + shockFlags = shockFlattening->flagArray; + } + slopeLimIDIR = std::make_unique>(data, Vc, shockFlags, + haveShockFlattening); #if DIMENSIONS >= 2 - slopeLimJDIR = std::make_unique>(this); + slopeLimJDIR = std::make_unique>(data, Vc, shockFlags, + haveShockFlattening); #endif #if DIMENSIONS == 3 - slopeLimKDIR = std::make_unique>(this); + slopeLimKDIR = std::make_unique>(data, Vc, shockFlags, + haveShockFlattening); #endif } diff --git a/src/fluid/RiemannSolver/shockFlag.hpp b/src/fluid/RiemannSolver/shockFlag.hpp new file mode 100644 index 00000000..33e35c59 --- /dev/null +++ b/src/fluid/RiemannSolver/shockFlag.hpp @@ -0,0 +1,33 @@ +// *********************************************************************************** +// Idefix MHD astrophysical code +// +// Source file src/fluid/RiemannSolver/shockFlag.hpp +// +// Last modified : 10/2026 +// +// Copyright(C) by : +// - Geoffroy Lesur (IPAG/UGA/CNRS - 2023) +// and other code contributors +// +// Licensed under CeCILL 2.1 License, see COPYING for more information +// *********************************************************************************** + +#ifndef FLUID_RIEMANNSOLVER_SHOCKFLAG_HPP_ +#define FLUID_RIEMANNSOLVER_SHOCKFLAG_HPP_ + +#include + +#include "arrays.hpp" + +// Forward declaration (used in the enrollment function signature) +class DataBlock; + +// Flag attached to each cell: whether shock flattening should be applied there +enum class FlagShock{None, Shock}; + +// User-defined function used to flag additional cells for shock flattening. +// It is called on the host and receives the datablock, the current time and +// the (writable) shock flag array. +using UserShockFunc = std::function&)>; + +#endif // FLUID_RIEMANNSOLVER_SHOCKFLAG_HPP_ diff --git a/src/fluid/RiemannSolver/shockFlattening.hpp b/src/fluid/RiemannSolver/shockFlattening.hpp index 454b7db0..b3fc0c9d 100644 --- a/src/fluid/RiemannSolver/shockFlattening.hpp +++ b/src/fluid/RiemannSolver/shockFlattening.hpp @@ -3,7 +3,7 @@ // // Source file src/fluid/RiemannSolver/shockFlattening.hpp // -// Last modified : 10/2023 +// Last modified : 10/2026 // // Copyright(C) by : // - Geoffroy Lesur (IPAG/UGA/CNRS - 2023) @@ -15,31 +15,50 @@ #ifndef FLUID_RIEMANNSOLVER_SHOCKFLATTENING_HPP_ #define FLUID_RIEMANNSOLVER_SHOCKFLATTENING_HPP_ +#include + #include "fluid.hpp" #include "../physics.hpp" #include "eos.hpp" +#include "shockFlag.hpp" -enum class FlagShock{None, Shock}; - -using UserShockFunc = void (*) (DataBlock &, const real t, IdefixArray3D&); +// Forward declaration of the functor performing the actual shock detection +template +struct ShockFlattening_FindShockFunctor; template class ShockFlattening { + // The find-shock functor reads the private members of this class + template + friend struct ShockFlattening_FindShockFunctor; + public: ShockFlattening(Fluid*, real); ShockFlattening() {} + // Compute the shock flag array from the current primitive state void FindShock(); - Fluid *hydro; - EquationOfState *eos; + // Enroll a user-defined function flagging additional cells for shock flattening. + // The function is called on the host at the end of FindShock() and receives the + // datablock, the current time and the flag array (see UserShockFunc). + template + void EnrollUserShockFlag(Func myFunc); + + bool IsActive() const {return isActive;} + real Smoothing() const {return smoothing;} + + // Shock flag array. Kept public so that it can be passed to (and modified by) + // user-enrolled shock flag functions. IdefixArray3D flagArray; + + private: + Fluid *hydro{nullptr}; bool isActive{false}; real smoothing{0}; - void EnrollUserShockFlag(UserShockFunc); bool haveUserShockFlag{false}; - UserShockFunc userShockFunc{NULL}; + UserShockFunc userShockFunc; }; template @@ -47,7 +66,7 @@ struct ShockFlattening_FindShockFunctor { //***************************************************************** // Functor constructor //***************************************************************** - explicit ShockFlattening_FindShockFunctor(ShockFlattening *sf) { + explicit ShockFlattening_FindShockFunctor(ShockFlattening *sf) { flags = sf->flagArray; Vc = sf->hydro->Vc; smoothing = sf->smoothing; @@ -66,8 +85,9 @@ struct ShockFlattening_FindShockFunctor { eos = *(sf->hydro->eos.get()); } } + //***************************************************************** - // Functor Variables + // Functor variables //***************************************************************** real smoothing; IdefixArray3D flags; @@ -75,20 +95,19 @@ struct ShockFlattening_FindShockFunctor { #if GEOMETRY == CARTESIAN IdefixArray1D dx1, dx2, dx3; #else - IdefixArray3D Ax1,Ax2,Ax3, dV; + IdefixArray3D Ax1, Ax2, Ax3, dV; #endif - EquationOfState eos; + EquationOfStateDevice eos; //***************************************************************** - // Functor Operator + // Divergence of the velocity field at cell (k,j,i) //***************************************************************** - KOKKOS_INLINE_FUNCTION void operator() (const int k, const int j, const int i) const { - flags(k,j,i) = FlagShock::None; + KOKKOS_INLINE_FUNCTION real ComputeDivV(const int k, const int j, const int i) const { #if GEOMETRY == CARTESIAN - real divV = D_EXPAND( (Vc(VX1,k,j,i+1) - Vc(VX1,k,j,i-1))/(2*dx1(i)) , - + (Vc(VX2,k,j+1,i) - Vc(VX2,k,j-1,i))/(2*dx2(j)) , - + (Vc(VX3,k+1,j,i) - Vc(VX3,k-1,j,i))/(2*dx3(k)) ); + return D_EXPAND( (Vc(VX1,k,j,i+1) - Vc(VX1,k,j,i-1))/(2*dx1(i)) , + + (Vc(VX2,k,j+1,i) - Vc(VX2,k,j-1,i))/(2*dx2(j)) , + + (Vc(VX3,k+1,j,i) - Vc(VX3,k-1,j,i))/(2*dx3(k)) ); #else real divV = D_EXPAND( (Vc(VX1,k,j,i+1) + Vc(VX1,k,j,i))*Ax1(k,j,i+1) - (Vc(VX1,k,j,i-1) + Vc(VX1,k,j,i))*Ax1(k,j,i) , @@ -96,62 +115,80 @@ struct ShockFlattening_FindShockFunctor { (Vc(VX2,k,j-1,i) + Vc(VX2,k,j,i))*Ax2(k,j,i) , + (Vc(VX3,k+1,j,i) + Vc(VX3,k,j,i))*Ax3(k+1,j,i) - (Vc(VX3,k-1,j,i) + Vc(VX3,k,j,i))*Ax3(k,j,i) ); - divV = 0.5*divV/dV(k,j,i); + return 0.5*divV/dV(k,j,i); #endif + } - if(divV= 2 + cs = eos.GetWaveSpeed(k,j+1,i); + pR = Vc(RHO,k,j+1,i)*cs*cs; + cs = eos.GetWaveSpeed(k,j-1,i); + pL = Vc(RHO,k,j-1,i)*cs*cs; pmin = FMIN(pmin,pL); pmin = FMIN(pmin,pR); - gradP = FABS(pR - pL); - #if DIMENSIONS >= 2 - cs = eos.GetWaveSpeed(k,j+1,i); - pR = Vc(RHO,k,j+1,i)*cs*cs; - cs = eos.GetWaveSpeed(k,j-1,i); - pL = Vc(RHO,k,j-1,i)*cs*cs; - - pmin = FMIN(pmin,pL); - pmin = FMIN(pmin,pR); - - gradP += FABS(pR - pL); - #endif - #if DIMENSIONS == 3 - cs = eos.GetWaveSpeed(k+1,j,i); - pR = Vc(RHO,k+1,j,i)*cs*cs; - cs = eos.GetWaveSpeed(k-1,j,i); - pL = Vc(RHO,k-1,j,i)*cs*cs; - - pmin = FMIN(pmin,pL); - pmin = FMIN(pmin,pR); - - gradP += FABS(pR - pL); - #endif - } else if constexpr(Phys::pressure) { - pmin = Vc(PRS,k,j,i); - pmin = FMIN(pmin,Vc(PRS,k,j,i+1)); - pmin = FMIN(pmin,Vc(PRS,k,j,i-1)); - gradP = FABS(Vc(PRS,k,j,i+1) - Vc(PRS,k,j,i-1)); - #if DIMENSIONS >= 2 - pmin = FMIN(pmin,Vc(PRS,k,j+1,i)); - pmin = FMIN(pmin,Vc(PRS,k,j-1,i)); - gradP += FABS(Vc(PRS,k,j+1,i) - Vc(PRS,k,j-1,i)); - #endif - #if DIMENSIONS == 3 - pmin = FMIN(pmin,Vc(PRS,k+1,j,i)); - pmin = FMIN(pmin,Vc(PRS,k-1,j,i)); - gradP += FABS(Vc(PRS,k+1,j,i) - Vc(PRS,k-1,j,i)); - #endif - } + gradP += FABS(pR - pL); + #endif + #if DIMENSIONS == 3 + cs = eos.GetWaveSpeed(k+1,j,i); + pR = Vc(RHO,k+1,j,i)*cs*cs; + cs = eos.GetWaveSpeed(k-1,j,i); + pL = Vc(RHO,k-1,j,i)*cs*cs; + + pmin = FMIN(pmin,pL); + pmin = FMIN(pmin,pR); + + gradP += FABS(pR - pL); + #endif + } else if constexpr(Phys::pressure) { + pmin = Vc(PRS,k,j,i); + pmin = FMIN(pmin,Vc(PRS,k,j,i+1)); + pmin = FMIN(pmin,Vc(PRS,k,j,i-1)); + gradP = FABS(Vc(PRS,k,j,i+1) - Vc(PRS,k,j,i-1)); + #if DIMENSIONS >= 2 + pmin = FMIN(pmin,Vc(PRS,k,j+1,i)); + pmin = FMIN(pmin,Vc(PRS,k,j-1,i)); + gradP += FABS(Vc(PRS,k,j+1,i) - Vc(PRS,k,j-1,i)); + #endif + #if DIMENSIONS == 3 + pmin = FMIN(pmin,Vc(PRS,k+1,j,i)); + pmin = FMIN(pmin,Vc(PRS,k-1,j,i)); + gradP += FABS(Vc(PRS,k+1,j,i) - Vc(PRS,k-1,j,i)); + #endif + } + } + + //***************************************************************** + // Functor operator + //***************************************************************** + KOKKOS_INLINE_FUNCTION void operator() (const int k, const int j, const int i) const { + flags(k,j,i) = FlagShock::None; + + const real divV = ComputeDivV(k,j,i); + if(divV < ZERO_F) { + // Compressing flow: look for a strong pressure jump if constexpr(Phys::pressure || Phys::isothermal) { + real pmin, gradP; + ComputePressureCriterion(k,j,i,pmin,gradP); if(gradP > smoothing*pmin) { flags(k,j,i) = FlagShock::Shock; } @@ -192,7 +229,8 @@ void ShockFlattening::FindShock() { } template -void ShockFlattening::EnrollUserShockFlag(UserShockFunc myFunc) { +template +void ShockFlattening::EnrollUserShockFlag(Func myFunc) { this->haveUserShockFlag = true; this->userShockFunc = myFunc; } diff --git a/src/fluid/addSourceTerms.hpp b/src/fluid/addSourceTerms.hpp index 0350fab0..d5b7a8a1 100644 --- a/src/fluid/addSourceTerms.hpp +++ b/src/fluid/addSourceTerms.hpp @@ -67,7 +67,7 @@ struct Fluid_AddSourceTermsFunctor { IdefixArray1D rt; #endif - EquationOfState eos; + EquationOfStateDevice eos; bool haveRotation; real OmegaZ; @@ -208,12 +208,7 @@ void Fluid::AddSourceTerms(real t, real dt) { idfx::pushRegion("Fluid::AddSourceTerms"); if(haveUserSourceTerm) { - if(userSourceTerm != NULL) { - userSourceTerm(this, t, dt); - } else { - // Deprecated version - userSourceTermOld(*data, t, dt); - } + userSourceTerm(this, t, dt); } auto func = Fluid_AddSourceTermsFunctor(this,dt); diff --git a/src/fluid/boundary/boundary.hpp b/src/fluid/boundary/boundary.hpp index 56c7cac2..236e1b6c 100644 --- a/src/fluid/boundary/boundary.hpp +++ b/src/fluid/boundary/boundary.hpp @@ -18,6 +18,7 @@ #include #include #include +#include #include "idefix.hpp" #include "fluid_defs.hpp" #include "grid.hpp" @@ -34,13 +35,13 @@ class Axis; // User-defined boundary conditions signature template -using UserDefBoundaryFunc = void (*) (Fluid *, int dir, BoundarySide side, const real t); -using UserDefBoundaryFuncOld = void (*) (DataBlock &, int dir, BoundarySide side, const real t); +using UserDefBoundaryFunc = std::function *, int, BoundarySide, const real)>; +using UserDefBoundaryFuncOld = std::function; // This last one is deprecated template -using InternalBoundaryFunc = void (*) (Fluid *, const real t); -using InternalBoundaryFuncOld = void (*) (DataBlock &, const real t); // DEPRECATED +using InternalBoundaryFunc = std::function *, const real)>; +using InternalBoundaryFuncOld = std::function; // DEPRECATED using BoundingBox = std::array,3>; @@ -55,12 +56,12 @@ class Boundary { void EnforceFluxBoundaries(int,real); ///< Apply boundary condition conditions to the fluxes - void EnrollUserDefBoundary(UserDefBoundaryFuncOld); ///< Deprecated - void EnrollUserDefBoundary(UserDefBoundaryFunc); ///< User-defined boundary condition - void EnrollInternalBoundary(InternalBoundaryFuncOld); ///< Deprecated - void EnrollInternalBoundary(InternalBoundaryFunc); ///< User-defined internal boundary - void EnrollFluxBoundary(UserDefBoundaryFuncOld); ///< Deprecated - void EnrollFluxBoundary(UserDefBoundaryFunc); ///< Flux boundary condition + template + void EnrollUserDefBoundary(Func myFunc); ///< User-defined boundary condition + template + void EnrollInternalBoundary(Func myFunc); ///< User-defined internal boundary + template + void EnrollFluxBoundary(Func myFunc); ///< Flux boundary condition void EnforcePeriodic(int, BoundarySide ); ///< Enforce periodic BC in direction and side void EnforceReflective(int, BoundarySide ); ///< Enforce reflective BC in direction and side @@ -72,19 +73,16 @@ class Boundary { #endif // User defined Boundary conditions - UserDefBoundaryFuncOld userDefBoundaryFuncOld{NULL}; - UserDefBoundaryFunc userDefBoundaryFunc{NULL}; + UserDefBoundaryFunc userDefBoundaryFunc; bool haveUserDefBoundary{false}; // Internal boundary function bool haveInternalBoundary{false}; - InternalBoundaryFuncOld internalBoundaryFuncOld{NULL}; - InternalBoundaryFunc internalBoundaryFunc{NULL}; + InternalBoundaryFunc internalBoundaryFunc; // Flux boundary function bool haveFluxBoundary{false}; - UserDefBoundaryFuncOld fluxBoundaryFuncOld{NULL}; - UserDefBoundaryFunc fluxBoundaryFunc{NULL}; + UserDefBoundaryFunc fluxBoundaryFunc; // specific for loops on ghost cells template @@ -261,21 +259,8 @@ Boundary::Boundary(Fluid* fluid) { template -void Boundary::EnrollFluxBoundary(UserDefBoundaryFuncOld myFunc) { - std::stringstream msg; - msg << "The old signature for flux boundary condition " << std::endl - << "(DataBlock &, int dir, BoundarySide side, const real t)" << std::endl - << "is deprecated. You should now use "<< std::endl - << "(Fluid *, int dir, BoundarySide side, const real t)" << std::endl - << "With the Phys of your choice (DefaultPhysics, DustPhysics...)" << std::endl; - - IDEFIX_WARNING(msg); - this->fluxBoundaryFuncOld = myFunc; - this->haveFluxBoundary = true; -} - -template -void Boundary::EnrollFluxBoundary(UserDefBoundaryFunc myFunc) { +template +void Boundary::EnrollFluxBoundary(Func myFunc) { this->haveFluxBoundary = true; this->fluxBoundaryFunc = myFunc; } @@ -285,18 +270,10 @@ void Boundary::EnforceFluxBoundaries(int dir,const real t) { idfx::pushRegion("Boundary::EnforceFluxBoundaries"); if(haveFluxBoundary) { if(data->lbound[dir] != internal) { - if(fluxBoundaryFunc != NULL) { - this->fluxBoundaryFunc(fluid, dir, left, t); - } else { - this->fluxBoundaryFuncOld(*data, dir, left, t); - } + this->fluxBoundaryFunc(fluid, dir, left, t); } if(data->rbound[dir] != internal) { - if(fluxBoundaryFunc != NULL) { - this->fluxBoundaryFunc(fluid, dir, right, t); - } else { - this->fluxBoundaryFuncOld(*data, dir, right, t); - } + this->fluxBoundaryFunc(fluid, dir, right, t); } } else { IDEFIX_ERROR("Cannot enforce flux boundary conditions without enrolling a specific function"); @@ -310,11 +287,7 @@ void Boundary::SetBoundaries(real t) { // set internal boundary conditions if(haveInternalBoundary) { idfx::pushRegion("Boundary::UserDefInternalBoundary"); - if(internalBoundaryFunc != NULL) { - internalBoundaryFunc(fluid, t); - } else { - internalBoundaryFuncOld(*data, t); - } + internalBoundaryFunc(fluid, t); idfx::popRegion(); } for(int dir=0 ; dir < DIMENSIONS ; dir++ ) { @@ -384,12 +357,7 @@ void Boundary::EnforceBoundaryDir(real t, int dir) { case BoundaryType::userdef: if(this->haveUserDefBoundary) { idfx::pushRegion("Boundary::UserDefBoundary"); - if(this->userDefBoundaryFunc != NULL) { - this->userDefBoundaryFunc(fluid, dir, left, t); - } else { - // Revert to deprecated Boundary condition - this->userDefBoundaryFuncOld(*data, dir, left, t); - } + this->userDefBoundaryFunc(fluid, dir, left, t); idfx::popRegion(); } else { std::stringstream msg; @@ -430,12 +398,7 @@ void Boundary::EnforceBoundaryDir(real t, int dir) { case BoundaryType::userdef: if(this->haveUserDefBoundary) { idfx::pushRegion("Boundary::UserDefBoundary"); - if(this->userDefBoundaryFunc != NULL) { - this->userDefBoundaryFunc(fluid, dir, right, t); - } else { - // Revert to deprecated Boundary condition - this->userDefBoundaryFuncOld(*data, dir, right, t); - } + this->userDefBoundaryFunc(fluid, dir, right, t); idfx::popRegion(); } else { std::stringstream msg; @@ -599,44 +562,16 @@ void Boundary::ReconstructNormalField(int dir) { template -void Boundary::EnrollUserDefBoundary(UserDefBoundaryFuncOld myFunc) { - std::stringstream msg; - msg << "The old signature for user-defined boundary condition " << std::endl - << "(DataBlock &, int dir, BoundarySide side, const real t)" << std::endl - << "is deprecated. You should now use "<< std::endl - << "(Fluid *, int dir, BoundarySide side, const real t)" << std::endl - << "With the Phys of your choice (DefaultPhysics, DustPhysics...)" << std::endl; - - IDEFIX_WARNING(msg); - this->userDefBoundaryFuncOld = myFunc; - this->haveUserDefBoundary = true; -} - -template -void Boundary::EnrollUserDefBoundary(UserDefBoundaryFunc myFunc) { - std::stringstream msg; - +template +void Boundary::EnrollUserDefBoundary(Func myFunc) { this->userDefBoundaryFunc = myFunc; this->haveUserDefBoundary = true; } template -void Boundary::EnrollInternalBoundary(InternalBoundaryFuncOld myFunc) { - std::stringstream msg; - msg << "The old signature for internal boundary condition " << std::endl - << "(DataBlock &, const real t)" << std::endl - << "is deprecated. You should now use "<< std::endl - << "(Fluid *, const real t)" << std::endl - << "With the Phys of your choice (DefaultPhysics, DustPhysics...)" << std::endl; - - IDEFIX_WARNING(msg); - this->internalBoundaryFuncOld = myFunc; - this->haveInternalBoundary = true; -} - -template -void Boundary::EnrollInternalBoundary(InternalBoundaryFunc myFunc) { +template +void Boundary::EnrollInternalBoundary(Func myFunc) { this->internalBoundaryFunc = myFunc; this->haveInternalBoundary = true; } diff --git a/src/fluid/braginskii/bragThermalDiffusion.hpp b/src/fluid/braginskii/bragThermalDiffusion.hpp index 06597661..b6dae331 100644 --- a/src/fluid/braginskii/bragThermalDiffusion.hpp +++ b/src/fluid/braginskii/bragThermalDiffusion.hpp @@ -81,7 +81,7 @@ class BragThermalDiffusion { real clessAlpha, clessBeta; // equation of state (required to get the heat capacity) - EquationOfState *eos; + EquationOfStateDevice *eos; PLMLimiter limiter{PLMLimiter::VanLeer}; }; @@ -260,7 +260,7 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, IdefixArray4D Vc = this->Vc; IdefixArray4D Vs = this->Vs; IdefixArray3D dMax = this->dMax; - EquationOfState eos = *(this->eos); + EquationOfStateDevice eos = *(this->eos); HydroModuleStatus haveThermalDiffusion = this->status.status; bool haveSlopeLimiter = this->haveSlopeLimiter; diff --git a/src/fluid/convertConsToPrim.hpp b/src/fluid/convertConsToPrim.hpp index 3a201b61..46af1acf 100644 --- a/src/fluid/convertConsToPrim.hpp +++ b/src/fluid/convertConsToPrim.hpp @@ -20,7 +20,7 @@ #include "tracer.hpp" template -KOKKOS_INLINE_FUNCTION void K_ConsToPrim(real Vc[], real Uc[], const EquationOfState *eos) { +KOKKOS_INLINE_FUNCTION void K_ConsToPrim(real Vc[], real Uc[], const EquationOfStateDevice *eos) { Vc[RHO] = Uc[RHO]; EXPAND( Vc[VX1] = Uc[MX1]/Uc[RHO]; , @@ -74,7 +74,7 @@ KOKKOS_INLINE_FUNCTION void K_ConsToPrim(real Vc[], real Uc[], const EquationOfS } template -KOKKOS_INLINE_FUNCTION void K_PrimToCons(real Uc[], real Vc[], const EquationOfState *eos) { +KOKKOS_INLINE_FUNCTION void K_PrimToCons(real Uc[], real Vc[], const EquationOfStateDevice *eos) { Uc[RHO] = Vc[RHO]; EXPAND( Uc[MX1] = Vc[VX1]*Vc[RHO]; , @@ -115,7 +115,7 @@ void Fluid::ConvertConsToPrim() { IdefixArray4D Vc = this->Vc; IdefixArray4D Uc = this->Uc; - EquationOfState eos; + EquationOfStateDevice eos; if constexpr(Phys::eos) { eos = *(this->eos.get()); } @@ -162,7 +162,7 @@ void Fluid::ConvertPrimToCons() { IdefixArray4D Vc = this->Vc; IdefixArray4D Uc = this->Uc; - EquationOfState eos; + EquationOfStateDevice eos; if constexpr(Phys::eos) { eos = *(this->eos.get()); } diff --git a/src/fluid/drag.cpp b/src/fluid/drag.cpp index 79464cdf..db96af8a 100644 --- a/src/fluid/drag.cpp +++ b/src/fluid/drag.cpp @@ -31,7 +31,7 @@ void Drag::AddDragForce(const real dt) { } auto gammaDrag = this->gammaDrag; - gammaDrag.RefreshUserDrag(data); + RefreshUserDrag(gammaDrag); // Compute a drag force fd = - gamma*rhod*rhog*(vd-vg) // Where gamma is computed according to the choice of drag type @@ -88,7 +88,7 @@ void Drag::AddImplicitBackReaction(const real dt, IdefixArray3D preFactor) } auto gammaDrag = this->gammaDrag; - gammaDrag.RefreshUserDrag(data); + RefreshUserDrag(gammaDrag); // Compute a drag force fd = - gamma*rhod*rhog*(vd-vg) // Where gamma is computed according to the choice of drag type @@ -154,7 +154,7 @@ void Drag::AddImplicitFluidMomentum(const real dt) { } auto gammaDrag = this->gammaDrag; - if(!feedback) gammaDrag.RefreshUserDrag(data); + if(!feedback) RefreshUserDrag(gammaDrag); // Compute a drag force fd = - gamma*rhod*rhog*(vd-vg) // Where gamma is computed according to the choice of drag type @@ -182,6 +182,18 @@ void Drag::AddImplicitFluidMomentum(const real dt) { idfx::popRegion(); } +void Drag::RefreshUserDrag(GammaDrag &gammaDrag) { + if(dragType == GammaDrag::Type::Userdef) { + if(userDrag) { + idfx::pushRegion("Drag::UserDrag"); + userDrag(data, gammaDrag.instanceNumber, gammaDrag.dragCoeff, gammaDrag.gammai); + idfx::popRegion(); + } else { + IDEFIX_ERROR("No User-defined drag function has been enrolled"); + } + } +} + void Drag::ShowConfig() { idfx::cout << "Drag: Using "; if(implicit) { @@ -212,32 +224,14 @@ void Drag::ShowConfig() { } } -void Drag::EnrollUserDrag(UserDefDragFunc func) { - gammaDrag.EnrollUserDrag(func); -} + //////////////////////////////////////////// // GammaDrag function definitions //////////////////////////////////////////// -void GammaDrag::RefreshUserDrag(DataBlock *data) { - if(type == Type::Userdef) { - if(userDrag != NULL) { - idfx::pushRegion("GammaDrag::UserDrag"); - userDrag(data, instanceNumber, dragCoeff, gammai); - idfx::popRegion(); - } else { - IDEFIX_ERROR("No User-defined drag function has been enrolled"); - } - } -} -void GammaDrag::EnrollUserDrag(UserDefDragFunc func) { - if(type != Type::Userdef) { - IDEFIX_ERROR("User-defined drag function requires drag entry to be set to \"userdef\""); - } - this->userDrag = func; -} + GammaDrag::GammaDrag(Input &input, std::string BlockName, int instanceNumber, DataBlock *data) { if(input.CheckEntry(BlockName,"drag")>=0) { diff --git a/src/fluid/drag.hpp b/src/fluid/drag.hpp index 786d48d6..bcdae56d 100644 --- a/src/fluid/drag.hpp +++ b/src/fluid/drag.hpp @@ -15,12 +15,14 @@ #define FLUID_DRAG_HPP_ #include +#include + #include "idefix.hpp" #include "input.hpp" #include "fluid_defs.hpp" #include "eos.hpp" -using UserDefDragFunc = void (*) (DataBlock *, int n, real beta, IdefixArray3D &gammai); +using UserDefDragFunc = std::function &)>; /* A class that holds the kind of drag law we intend to use */ class GammaDrag { @@ -28,8 +30,6 @@ class GammaDrag { enum class Type{Gamma, Tau, Size, Userdef}; GammaDrag() = default; GammaDrag(Input &, std::string BlockName, int instanceNumber, DataBlock *data); - void EnrollUserDrag(UserDefDragFunc); - void RefreshUserDrag(DataBlock *); KOKKOS_INLINE_FUNCTION real GetGamma(const int k, const int j, const int i) const { real gamma; // The drag coefficient @@ -55,15 +55,13 @@ class GammaDrag { } return gamma; } - Type type; real dragCoeff; - EquationOfState eos; + EquationOfStateDevice eos; IdefixArray3D gammai; IdefixArray4D VcGas; int instanceNumber; - UserDefDragFunc userDrag{NULL}; }; @@ -75,6 +73,7 @@ class Drag { void ShowConfig(); // print configuration void AddDragForce(const real); + void RefreshUserDrag(GammaDrag &); ////////////////////////// // Implicit functions void AddImplicitBackReaction(const real, IdefixArray3D); // Add the back reaction @@ -82,7 +81,13 @@ class Drag { void AddImplicitFluidMomentum(const real); // Add the implicit drag force on dust grains ///////////////////////// - void EnrollUserDrag(UserDefDragFunc); // User defined drag function enrollment + template + void EnrollUserDrag(Func func) { + if(dragType != GammaDrag::Type::Userdef) { + IDEFIX_ERROR("User-defined drag function requires drag entry to be set to \"userdef\""); + } + this->userDrag = func; + } bool IsImplicit() const { return implicit; } // Check if the drag is implicit IdefixArray4D UcDust; // Dust conservative quantities @@ -99,6 +104,8 @@ class Drag { bool feedback{false}; bool implicit{false}; int instanceNumber; + UserDefDragFunc userDrag; + GammaDrag::Type dragType; }; @@ -127,7 +134,8 @@ Drag::Drag(Input &input, Fluid *hydroin): if(input.CheckEntry(blockName,"drag")>=0) { this->instanceNumber = hydroin->instanceNumber; this->gammaDrag = GammaDrag(input,blockName,instanceNumber,data); - + this->gammaDrag.eos = *(data->hydro->eos.get()); + this->dragType = gammaDrag.type; // Feedback is true by default, but can be switched off. this->feedback = input.GetOrSet(blockName,"drag_feedback",0,true); this->implicit = input.GetOrSet(blockName,"drag_implicit",0,false); diff --git a/src/fluid/enroll.hpp b/src/fluid/enroll.hpp index 68406d56..d550e1ae 100644 --- a/src/fluid/enroll.hpp +++ b/src/fluid/enroll.hpp @@ -50,29 +50,13 @@ void Fluid::EnrollFluxBoundary(T myFunc) { } template -void Fluid::EnrollUserSourceTerm(SrcTermFunc myFunc) { +template +void Fluid::EnrollUserSourceTerm(Func myFunc) { this->userSourceTerm = myFunc; this->haveUserSourceTerm = true; this->haveSourceTerms = true; } -// Deprecated enrollment function -template -void Fluid::EnrollUserSourceTerm(SrcTermFuncOld myFunc) { - std::stringstream msg; - msg << "The old signature for user-defined source terms " << std::endl - << "(DataBlock &, const real t, const real dt)" << std::endl - << "is deprecated. You should now use "<< std::endl - << "(Fluid *, const real t, const real dt)" << std::endl - << "With the Phys of your choice (DefaultPhysics, DustPhysics...)" << std::endl; - - IDEFIX_WARNING(msg); - - this->userSourceTermOld = myFunc; - this->haveUserSourceTerm = true; - this->haveSourceTerms = true; -} - template void Fluid::EnrollEmfBoundary(EmfBoundaryFunc myFunc) { diff --git a/src/fluid/eos/eos_adiabatic.hpp b/src/fluid/eos/eos_adiabatic.hpp index 49b35a4e..21b53483 100644 --- a/src/fluid/eos/eos_adiabatic.hpp +++ b/src/fluid/eos/eos_adiabatic.hpp @@ -19,23 +19,19 @@ #include "idefix.hpp" #include "input.hpp" -// This is the adiabatic (constant gamma) implementation of the equation of state -class EquationOfState { +// This is the adiabatic (constant gamma) implementation +// of the equation of state (device-compatible subset) +class EquationOfStateDevice { public: - EquationOfState() = default; + EquationOfStateDevice() = default; - EquationOfState(Input & input, DataBlock *, std::string prefix) { + EquationOfStateDevice(Input & input, DataBlock *, std::string prefix) { this->gamma = input.GetOrSet(prefix,"gamma",0, 5.0/3.0); } - void ShowConfig() { - idfx::cout << "EquationOfState: ideal with gamma=" << this->gamma << std::endl; - } - // First adiabatic exponent. In the ideal EOS, gamma does not depend on the gas state, // So we add default values to 0 here so that GetGamma can be called without any argument KOKKOS_INLINE_FUNCTION real GetGamma(real P = 0.0, real rho = 0.0) const {return gamma;} - void Refresh(DataBlock &, real) {} // Refresh the eos (recompute coefficients and tables) KOKKOS_INLINE_FUNCTION real GetWaveSpeed(int k, int j, int i) const { @@ -51,7 +47,22 @@ class EquationOfState { return Eint * (gamma-1.0); } - private: + protected: real gamma; }; + +// This is the adiabatic (constant gamma) implementation of the equation of state (full class) +class EquationOfState : public EquationOfStateDevice { + public: + EquationOfState() = default; + + EquationOfState(Input & input, DataBlock *data, std::string prefix) + : EquationOfStateDevice(input, data, prefix) {} + + void ShowConfig() { + idfx::cout << "EquationOfState: ideal with gamma=" << this->gamma << std::endl; + } + + void Refresh(DataBlock &, real) {} // Refresh the eos (recompute coefficients and tables) +}; #endif // FLUID_EOS_EOS_ADIABATIC_HPP_ diff --git a/src/fluid/eos/eos_isothermal.hpp b/src/fluid/eos/eos_isothermal.hpp index 72ca0604..6667f82d 100644 --- a/src/fluid/eos/eos_isothermal.hpp +++ b/src/fluid/eos/eos_isothermal.hpp @@ -22,24 +22,55 @@ #include "dataBlock.hpp" -// This is the isothermal implementation of the equation of state -class EquationOfState { +// This is the isothermal implementation of the equation of state (device-compatible subset) +class EquationOfStateDevice { public: - EquationOfState() = default; + EquationOfStateDevice() = default; + + // In the isothermal EOS, gamma does not depend on the gas state, + // So we add default values to 0 here so that GetGamma can be called without any argument + KOKKOS_INLINE_FUNCTION real GetGamma(real P = 0.0, real rho = 0.0) const { + return 1.0; + } - EquationOfState(Input & input, DataBlock *data, std::string prefix) { - std::string isoString = input.Get(prefix,"csiso",0); - if(isoString.compare("constant") == 0) { - this->haveIsoSoundSpeed = Constant; - this->isoSoundSpeed = input.Get(prefix,"csiso",1); - } else if(isoString.compare("userdef") == 0) { - this->haveIsoSoundSpeed = UserDefFunction; - this->isoSoundSpeedArray = IdefixArray3D(prefix+"_csIso", - data->np_tot[KDIR], data->np_tot[JDIR], data->np_tot[IDIR]); + KOKKOS_INLINE_FUNCTION real GetWaveSpeed(const int k, const int j, const int i) const { + if(haveIsoSoundSpeed == UserDefFunction) { + return isoSoundSpeedArray(k,j,i); } else { - IDEFIX_ERROR("csiso admits only constant or userdef entries"); + return isoSoundSpeed; } } + KOKKOS_INLINE_FUNCTION real GetInternalEnergy(real P, real rho) const { + Kokkos::abort("Internal energy is not defined in the isothermal EOS"); + } + KOKKOS_INLINE_FUNCTION real GetPressure(real Eint, real rho) const { + Kokkos::abort("Pressure is not defined in the isothermal EOS"); + } + + protected: + real isoSoundSpeed; + HydroModuleStatus haveIsoSoundSpeed{Disabled}; + IdefixArray3D isoSoundSpeedArray; +}; + +// This is the isothermal implementation of the equation of state (full class) +class EquationOfState : public EquationOfStateDevice { + public: + EquationOfState() = default; + EquationOfState(Input & input, DataBlock *data, std::string prefix) + : EquationOfStateDevice() { + std::string isoString = input.Get(prefix,"csiso",0); + if(isoString.compare("constant") == 0) { + this->haveIsoSoundSpeed = Constant; + this->isoSoundSpeed = input.Get(prefix,"csiso",1); + } else if(isoString.compare("userdef") == 0) { + this->haveIsoSoundSpeed = UserDefFunction; + this->isoSoundSpeedArray = IdefixArray3D(prefix+"_csIso", + data->np_tot[KDIR], data->np_tot[JDIR], data->np_tot[IDIR]); + } else { + IDEFIX_ERROR("csiso admits only constant or userdef entries"); + } + } void ShowConfig() { if(haveIsoSoundSpeed == Constant) { @@ -54,12 +85,6 @@ class EquationOfState { } } - // In the isothermal EOS, gamma does not depend on the gas state, - // So we add default values to 0 here so that GetGamma can be called without any argument - KOKKOS_INLINE_FUNCTION real GetGamma(real P = 0.0, real rho = 0.0) const { - return 1.0; - } - void Refresh(DataBlock &data, real t) { // Refresh the coefficients (and tables) idfx::pushRegion("EquationOfState::Refresh"); if(haveIsoSoundSpeed == UserDefFunction) { @@ -74,22 +99,9 @@ class EquationOfState { idfx::popRegion(); } - KOKKOS_INLINE_FUNCTION real GetWaveSpeed(const int k, const int j, const int i) const { - if(haveIsoSoundSpeed == UserDefFunction) { - return isoSoundSpeedArray(k,j,i); - } else { - return isoSoundSpeed; - } - } - KOKKOS_INLINE_FUNCTION real GetInternalEnergy(real P, real rho) const { - Kokkos::abort("Internal energy is not defined in the isothermal EOS"); - } - KOKKOS_INLINE_FUNCTION real GetPressure(real Eint, real rho) const { - Kokkos::abort("Pressure is not defined in the isothermal EOS"); - } - // Enroll user-defined isothermal sound speed - void EnrollIsoSoundSpeed(IsoSoundSpeedFunc func) { + template + void EnrollIsoSoundSpeed(Func func) { if(this->haveIsoSoundSpeed != UserDefFunction) { IDEFIX_WARNING("Isothermal sound speed enrollment requires Hydro/csiso " " to be set to userdef in .ini file"); @@ -98,10 +110,7 @@ class EquationOfState { } private: - real isoSoundSpeed; - HydroModuleStatus haveIsoSoundSpeed{Disabled}; - IdefixArray3D isoSoundSpeedArray; - IsoSoundSpeedFunc isoSoundSpeedFunc{NULL}; + IsoSoundSpeedFunc isoSoundSpeedFunc; }; #endif // FLUID_EOS_EOS_ISOTHERMAL_HPP_ diff --git a/src/fluid/eos/eos_template.hpp b/src/fluid/eos/eos_template.hpp index df5e1c06..84eb7099 100644 --- a/src/fluid/eos/eos_template.hpp +++ b/src/fluid/eos/eos_template.hpp @@ -19,21 +19,16 @@ #include "idefix.hpp" #include "input.hpp" -// This is a template for a generalized custom equation of state -class EquationOfState { +// This is a template for a generalized custom equation of state (device-compatible subset) +class EquationOfStateDevice { public: - EquationOfState() = default; + EquationOfStateDevice() = default; // Add here the required steps to initialize your equation of state - EquationOfState(Input & input, DataBlock *, std::string prefix) { + EquationOfStateDevice(Input & input, DataBlock *, std::string prefix) { // .... } - // Information message displayed before entering the main loop - void ShowConfig() { - idfx::cout << "EquationOfState: Custom." << std::endl; - } - // First adiabatic exponent. KOKKOS_INLINE_FUNCTION real GetGamma(real P , real rho ) const { real gamma; @@ -41,11 +36,6 @@ class EquationOfState { return gamma; } - // Refresh the eos (recompute coefficients and tables) - void Refresh(DataBlock &data, real t) { - // .... - } - // This function is used only when the isothermal approximation is enabled. Not needed here KOKKOS_INLINE_FUNCTION real GetWaveSpeed(int k, int j, int i) const { @@ -67,10 +57,30 @@ class EquationOfState { return P; } - private: + protected: real gamma; // Add here the internal variables required by your equation of state. }; +// This is a template for a generalized custom equation of state (full class) +class EquationOfState : public EquationOfStateDevice { + public: + EquationOfState() = default; + + // Add here the required steps to initialize your equation of state + EquationOfState(Input & input, DataBlock *data, std::string prefix) + : EquationOfStateDevice(input, data, prefix) {} + + // Information message displayed before entering the main loop + void ShowConfig() { + idfx::cout << "EquationOfState: Custom." << std::endl; + } + + // Refresh the eos (recompute coefficients and tables) + void Refresh(DataBlock &data, real t) { + // .... + } +}; + #endif // FLUID_EOS_EOS_TEMPLATE_HPP_ diff --git a/src/fluid/fluid.hpp b/src/fluid/fluid.hpp index d32b074f..1313039f 100644 --- a/src/fluid/fluid.hpp +++ b/src/fluid/fluid.hpp @@ -158,8 +158,8 @@ class Fluid { void EnrollEmfBoundary(EmfBoundaryFunc); // Add some user source terms - void EnrollUserSourceTerm(SrcTermFunc); - void EnrollUserSourceTerm(SrcTermFuncOld); // Deprecated + template + void EnrollUserSourceTerm(Func myFunc); // Enroll user-defined ohmic, ambipolar and Hall diffusivities void EnrollOhmicDiffusivity(DiffusivityFunc); @@ -227,19 +227,18 @@ class Fluid { // Emf boundary conditions bool haveEmfBoundary{false}; - EmfBoundaryFunc emfBoundaryFunc{NULL}; + EmfBoundaryFunc emfBoundaryFunc; // User defined source term - SrcTermFunc userSourceTerm{NULL}; - SrcTermFuncOld userSourceTermOld{NULL}; + SrcTermFunc userSourceTerm; bool haveUserSourceTerm{false}; real etaO, xH, xA; // Ohmic resistivity, Hall, ambipolar (when constant) // Ohmic, Hall and ambipolar diffusivity (when function-defined) - DiffusivityFunc ohmicDiffusivityFunc{NULL}; - DiffusivityFunc ambipolarDiffusivityFunc{NULL}; - DiffusivityFunc hallDiffusivityFunc{NULL}; + DiffusivityFunc ohmicDiffusivityFunc; + DiffusivityFunc ambipolarDiffusivityFunc; + DiffusivityFunc hallDiffusivityFunc; IdefixArray3D cMax; // Maximum propagation speed diff --git a/src/fluid/fluid_defs.hpp b/src/fluid/fluid_defs.hpp index 9acf98fb..e22f4e8d 100644 --- a/src/fluid/fluid_defs.hpp +++ b/src/fluid/fluid_defs.hpp @@ -18,6 +18,7 @@ #define FLUID_FLUID_DEFS_HPP_ #include +#include #include "../idefix.hpp" @@ -33,7 +34,6 @@ class Fluid; #define eps_UCT_CONTACT (1.0e-6) - // Parabolic terms can have different status enum HydroModuleStatus {Disabled, Constant, UserDefFunction}; @@ -45,26 +45,26 @@ struct ParabolicModuleStatus { }; -using GravPotentialFunc = void (*) (DataBlock &, const real t, IdefixArray1D&, - IdefixArray1D&, IdefixArray1D&, - IdefixArray3D &); +using GravPotentialFunc = std::function&, + IdefixArray1D&, IdefixArray1D&, + IdefixArray3D &)>; -using BodyForceFunc = void (*) (DataBlock &, const real t, IdefixArray4D&); +using BodyForceFunc = std::function&)>; -using IsoSoundSpeedFunc = void (*) (DataBlock &, const real t, IdefixArray3D &); +using IsoSoundSpeedFunc = std::function &)>; template -using SrcTermFunc = void (*) (Fluid *, const real t, const real dt); +using SrcTermFunc = std::function *, const real, const real)>; -using EmfBoundaryFunc = void (*) (DataBlock &, const real t); -using DiffusivityFunc = void (*) (DataBlock &, const real t, IdefixArray3D &); +using EmfBoundaryFunc = std::function; +using DiffusivityFunc = std::function &)>; -using BragDiffusivityFunc = void (*) (DataBlock &, const real, - std::vector> &); +using BragDiffusivityFunc = std::function> &)>; -// Deprecated signatures -using SrcTermFuncOld = void (*) (DataBlock &, const real t, const real dt); +// Deprecated signatures (kept as aliases for backward compatibility) +using SrcTermFuncOld = std::function; #endif //FLUID_FLUID_DEFS_HPP_ diff --git a/src/fluid/thermalDiffusion.cpp b/src/fluid/thermalDiffusion.cpp index b7abecb7..4ad08174 100644 --- a/src/fluid/thermalDiffusion.cpp +++ b/src/fluid/thermalDiffusion.cpp @@ -65,7 +65,7 @@ void ThermalDiffusion::AddDiffusiveFlux(int dir, const real t, const IdefixArray IdefixArray3D kappaArr = this->kappaArr; IdefixArray1D dx = this->data->dx[dir]; - EquationOfState eos = *(this->eos); + EquationOfStateDevice eos = *(this->eos); #if GEOMETRY == POLAR IdefixArray1D x1 = this->data->x[IDIR]; #endif diff --git a/src/fluid/thermalDiffusion.hpp b/src/fluid/thermalDiffusion.hpp index 9f45221e..7f1cfbc1 100644 --- a/src/fluid/thermalDiffusion.hpp +++ b/src/fluid/thermalDiffusion.hpp @@ -63,7 +63,7 @@ class ThermalDiffusion { real kappa; // equation of state (required to get the heat capacity) - EquationOfState *eos; + EquationOfStateDevice *eos; }; #include "fluid.hpp" diff --git a/src/gravity/gravity.hpp b/src/gravity/gravity.hpp index 40c5732d..c747c30b 100644 --- a/src/gravity/gravity.hpp +++ b/src/gravity/gravity.hpp @@ -15,6 +15,8 @@ #ifndef GRAVITY_GRAVITY_HPP_ #define GRAVITY_GRAVITY_HPP_ +#include + #include "idefix.hpp" #include "input.hpp" #include "selfGravity.hpp" @@ -23,11 +25,11 @@ class DataBlock; // Enrolled functions signature -using GravPotentialFunc = void (*) (DataBlock &, const real t, IdefixArray1D&, +using GravPotentialFunc = std::function&, IdefixArray1D&, IdefixArray1D&, - IdefixArray3D &); + IdefixArray3D &)>; -using BodyForceFunc = void (*) (DataBlock &, const real t, IdefixArray4D&); +using BodyForceFunc = std::function&)>; class Gravity { @@ -81,10 +83,10 @@ class Gravity { DataBlock *data; // User defined gravitational potential - GravPotentialFunc gravPotentialFunc{NULL}; + GravPotentialFunc gravPotentialFunc; // Body force - BodyForceFunc bodyForceFunc{NULL}; + BodyForceFunc bodyForceFunc; #ifdef DEBUG_GRAVITY // Used to get fields usefull for debugging diff --git a/src/gravity/laplacian.cpp b/src/gravity/laplacian.cpp index 2c7dafb7..ccc24098 100644 --- a/src/gravity/laplacian.cpp +++ b/src/gravity/laplacian.cpp @@ -632,10 +632,6 @@ void Laplacian::EnforceBoundary(int dir, BoundarySide side, LaplacianBoundaryTyp } -void Laplacian::EnrollUserDefBoundary(UserDefBoundaryFunc myFunc) { - this->userDefBoundaryFunc = myFunc; - this->haveUserDefBoundary = true; -} void Laplacian::SetBoundaries(IdefixArray3D &arr) { idfx::pushRegion("Laplacian::SetBoundaries"); diff --git a/src/gravity/laplacian.hpp b/src/gravity/laplacian.hpp index f3cd6a82..bd6908cd 100644 --- a/src/gravity/laplacian.hpp +++ b/src/gravity/laplacian.hpp @@ -16,6 +16,7 @@ #define GRAVITY_LAPLACIAN_HPP_ #include +#include #include "idefix.hpp" #ifdef WITH_MPI #include "mpi.hpp" @@ -52,16 +53,20 @@ class Laplacian { void operator() (IdefixArray3D in, IdefixArray3D laplacian); // Handling userdef boundary. - using UserDefBoundaryFunc = void (*) (DataBlock &, int dir, BoundarySide side, - const real t, IdefixArray3D &arr); + using UserDefBoundaryFunc = std::function &)>; void EnforceBoundary(int dir, BoundarySide side, LaplacianBoundaryType type, IdefixArray3D &); bool haveUserDefBoundary{false}; - void EnrollUserDefBoundary(UserDefBoundaryFunc); // Enroll user-defined boundary conditions + template + void EnrollUserDefBoundary(Func myFunc) { + this->userDefBoundaryFunc = myFunc; + this->haveUserDefBoundary = true; + } // User defined Boundary conditions - UserDefBoundaryFunc userDefBoundaryFunc{NULL}; + UserDefBoundaryFunc userDefBoundaryFunc; // Local potential array size std::array np_tot; diff --git a/src/main.cpp b/src/main.cpp index b46c8a43..38c3c0ac 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -34,6 +34,7 @@ #include #include #include +#include #include @@ -47,10 +48,14 @@ #include "dataBlock.hpp" #include "timeIntegrator.hpp" #include "setup.hpp" +#include "setup_registry.hpp" #include "output.hpp" #ifdef WITH_MPI #include "mpi.hpp" #endif +#ifdef WITH_PYTHON + #include "pydefix.hpp" +#endif @@ -109,7 +114,7 @@ int main( int argc, char* argv[] ) { Pydefix pydefix(input); #endif Output output(input, data); - Setup mysetup(input, grid, data, output); + auto mysetup = SetupRegistry::Create(input, grid, data, output); idfx::cout << "Main: initialisation finished." << std::endl; @@ -146,13 +151,13 @@ int main( int argc, char* argv[] ) { pydefix.InitFlow(data); } else { idfx::pushRegion("Setup::Initflow"); - mysetup.InitFlow(data); + mysetup->InitFlow(data); } data.DeriveVectorPotential(); idfx::popRegion(); #else idfx::pushRegion("Setup::Initflow"); - mysetup.InitFlow(data); + mysetup->InitFlow(data); data.DeriveVectorPotential(); idfx::popRegion(); #endif @@ -174,11 +179,11 @@ int main( int argc, char* argv[] ) { pydefix.InitFlow(data); } else { idfx::pushRegion("Setup::Initflow"); - mysetup.InitFlow(data); + mysetup->InitFlow(data); } #else idfx::pushRegion("Setup::Initflow"); - mysetup.InitFlow(data); + mysetup->InitFlow(data); #endif idfx::popRegion(); data.DeriveVectorPotential(); // This does something only when evolveVectorPotential is on diff --git a/src/output/output.cpp b/src/output/output.cpp index 30e84b46..11ce490c 100644 --- a/src/output/output.cpp +++ b/src/output/output.cpp @@ -399,34 +399,7 @@ void Output::ForceWriteXdmf(DataBlock &data) { } #endif -void Output::EnrollAnalysis(AnalysisFunc myFunc) { - idfx::pushRegion("Output::EnrollAnalysis"); - if(!analysisEnabled) { - IDEFIX_WARNING("You are enrolling an analysis function " - "but analysis are not enabled in the input file," - "so no analysis will be performed"); - } - analysisFunc = myFunc; - haveAnalysisFunc = true; - - idfx::popRegion(); -} -void Output::EnrollUserDefVariables(UserDefVariablesFunc myFunc) { - idfx::pushRegion("Output::EnrollUserDefVariable"); - if(!userDefVariablesEnabled) { - IDEFIX_ERROR("You are enrolling a user-defined variables function " - "but the userdef variables are not set in the input file"); - } - userDefVariablesFunc = myFunc; - if(haveSlices) { - for(int i = 0 ; i < slices.size() ; i++) { - slices[i]->EnrollUserDefFunc(myFunc); - } - } - haveUserDefVariablesFunc = true; - idfx::popRegion(); -} void Output::ResetTimer() { elapsedTime = 0.0; diff --git a/src/output/output.hpp b/src/output/output.hpp index 1f76aa94..2741c72b 100644 --- a/src/output/output.hpp +++ b/src/output/output.hpp @@ -6,8 +6,7 @@ // Last modified : 10/2024 // // Copyright(C) by : -// - Geoffroy Lesur (IPAG/UGA/CNRS - 2021 - 2024) -// - Alankar Dutta (MPA Garching - 2023) +// - Geoffroy Lesur (IPAG/UGA/CNRS - 2023 - 2024) // - Gaylor Wafflard (IPAG/UGA/CNRS - 2024) // and other code contributors // @@ -16,27 +15,29 @@ #ifndef OUTPUT_OUTPUT_HPP_ #define OUTPUT_OUTPUT_HPP_ -#include + #include #include +#include #include +#include #include "idefix.hpp" -#include "input.hpp" #include "dataBlock.hpp" +#include "input.hpp" #include "vtk.hpp" +#include "dump.hpp" #ifdef WITH_HDF5 #include "xdmf.hpp" #endif +#include "slice.hpp" #ifdef WITH_PYTHON -#include "pydefix.hpp" + #include "pydefix.hpp" #endif -#include "dump.hpp" -#include "slice.hpp" -using AnalysisFunc = void (*) (DataBlock &); +using AnalysisFunc = std::function; using UserDefVariablesContainer = std::map>; -using UserDefVariablesFunc = void (*) (DataBlock &, UserDefVariablesContainer &); +using UserDefVariablesFunc = std::function; class Output { @@ -57,8 +58,37 @@ class Output { #endif void ResetTimer(); // Reset internal timer double GetTimer(); - void EnrollAnalysis(AnalysisFunc); - void EnrollUserDefVariables(UserDefVariablesFunc); + + // Template enrollment methods (definitions inline) + template + void EnrollAnalysis(Func myFunc) { + idfx::pushRegion("Output::EnrollAnalysis"); + if(!analysisEnabled) { + IDEFIX_WARNING("You are enrolling an analysis function " + "but analysis are not enabled in the input file," + "so no analysis will be performed"); + } + analysisFunc = myFunc; + haveAnalysisFunc = true; + idfx::popRegion(); + } + + template + void EnrollUserDefVariables(Func myFunc) { + idfx::pushRegion("Output::EnrollUserDefVariable"); + if(!userDefVariablesEnabled) { + IDEFIX_ERROR("You are enrolling a user-defined variables function " + "but the userdef variables are not set in the input file"); + } + userDefVariablesFunc = myFunc; + if(haveSlices) { + for(int i = 0 ; i < slices.size() ; i++) { + slices[i]->EnrollUserDefFunc(myFunc); + } + } + haveUserDefVariablesFunc = true; + idfx::popRegion(); + } private: bool forceNoWrite = false; //< explicitely disable all writes @@ -102,4 +132,5 @@ class Output { Kokkos::Timer timer; double elapsedTime{0.0}; }; + #endif // OUTPUT_OUTPUT_HPP_ diff --git a/src/output/slice.cpp b/src/output/slice.cpp index eb83b650..2fe27e60 100644 --- a/src/output/slice.cpp +++ b/src/output/slice.cpp @@ -82,9 +82,6 @@ void Slice::EnrollUserDefVariables(std::map> this->haveUserDefinedVariables = true; } -void Slice::EnrollUserDefFunc(UserDefVariablesFunc myFunc) { - userDefVariablesFunc = myFunc; -} void Slice::CheckForWrite(DataBlock &data, bool force) { idfx::pushRegion("Slice:CheckForWrite"); diff --git a/src/output/slice.hpp b/src/output/slice.hpp index 2f72c7be..bb2f999a 100644 --- a/src/output/slice.hpp +++ b/src/output/slice.hpp @@ -27,17 +27,21 @@ class Grid; class SubGrid; class Vtk; -using UserDefVariablesFunc = void (*) (DataBlock &, - std::map> &); +using UserDefVariablesFunc = std::function> &)>; class Slice { public: Slice(Input &, DataBlock &, int, SliceType, int, real, real); void CheckForWrite(DataBlock &, bool = false); void EnrollUserDefVariables(std::map>); - void EnrollUserDefFunc(UserDefVariablesFunc); + template + void EnrollUserDefFunc(Func myFunc) { + userDefVariablesFunc = myFunc; + } real slicePeriod = 0.0; real sliceLast = 0.0; + private: bool containsX0; SliceType type; @@ -48,7 +52,7 @@ class Slice { bool haveUserDefinedVariables{false}; std::map> variableMap; std::map> userDefVariableMap; - UserDefVariablesFunc userDefVariablesFunc{NULL}; + UserDefVariablesFunc userDefVariablesFunc; #ifdef WITH_MPI MPI_Comm avgComm; // Communicator for averages #endif diff --git a/src/setup.cpp b/src/setup.cpp index 6713749b..ea848f3e 100644 --- a/src/setup.cpp +++ b/src/setup.cpp @@ -3,7 +3,7 @@ // // Source file src/setup.cpp // -// Last modified : 03/2026 +// Last modified : 10/2026 // // Copyright(C) by : // - Geoffroy Lesur (IPAG/UGA/CNRS - 2021 - 2026) @@ -12,6 +12,7 @@ // Licensed under CeCILL 2.1 License, see COPYING for more information // *********************************************************************************** +#include #include "setup.hpp" #include "idefix.hpp" @@ -28,5 +29,4 @@ __attribute__((weak)) void Setup::InitFlow(DataBlock &data) { IDEFIX_ERROR("Please create your own setup.cpp following the documentation"); } -__attribute__((weak)) Setup::~Setup() { -} +__attribute__((weak)) Setup::~Setup() {} diff --git a/src/setup.hpp b/src/setup.hpp index c2532260..46d6a80b 100644 --- a/src/setup.hpp +++ b/src/setup.hpp @@ -17,6 +17,7 @@ #ifndef SETUP_HPP_ #define SETUP_HPP_ +#include #include "idefix.hpp" #include "input.hpp" #include "grid.hpp" @@ -29,11 +30,46 @@ #include "dataBlockHost.hpp" #include "boundary.hpp" +// Forward declarations +class Setup; + +// Base class for user-defined setup with OOP style +class SetupBase { + public: + SetupBase(Input &input, Grid &grid, DataBlock &data, Output &output) + : myInput(input), myGrid(grid), myData(data), myOutput(output) {} + + virtual ~SetupBase() = default; + virtual void InitFlow(DataBlock &data) = 0; + + protected: + Input &myInput; + Grid &myGrid; + DataBlock &myData; + Output &myOutput; +}; + +// Default Setup class for backward compatibility (unchanged) class Setup { public: - Setup(Input &, Grid &, DataBlock &, Output&); - ~Setup(); - void InitFlow(DataBlock &); + Setup(Input &, Grid &, DataBlock &, Output&); + ~Setup(); + void InitFlow(DataBlock &); +}; + +// Default Setup wrapper that inherits from SetupBase +// This allows the factory to return SetupBase* for the default case +class DefaultSetup : public SetupBase { + public: + DefaultSetup(Input &input, Grid &grid, DataBlock &data, Output &output) + : SetupBase(input, grid, data, output), mySetup(input, grid, data, output) {} + + void InitFlow(DataBlock &data) override { + mySetup.InitFlow(data); + } + + private: + Setup mySetup; }; #endif // SETUP_HPP_ diff --git a/src/setup_registry.hpp b/src/setup_registry.hpp new file mode 100644 index 00000000..5b753f23 --- /dev/null +++ b/src/setup_registry.hpp @@ -0,0 +1,77 @@ +// *********************************************************************************** +// Idefix MHD astrophysical code +// +// Source file src/setup_registry.hpp +// +// Last modified : 10/2026 +// +// Copyright(C) by : +// - Geoffroy Lesur (IPAG/UGA/CNRS - 2021 - 2026) +// and other code contributors +// +// Licensed under CeCILL 2.1 License, see COPYING for more information +// *********************************************************************************** + +#ifndef SETUP_REGISTRY_HPP_ +#define SETUP_REGISTRY_HPP_ + +#include +#include +#include // for std::move +#include "setup.hpp" + +class SetupRegistry { + public: + using SetupFactory = std::function( + Input&, Grid&, DataBlock&, Output&)>; + + // Register the user's setup factory (only one allowed per binary) + static void Register(SetupFactory factory) { + auto& instance = GetInstance(); + if (instance.factory) { + IDEFIX_ERROR("Multiple Setup registrations detected. " + "Only one REGISTER_SETUP() call is allowed per binary."); + } + instance.factory = std::move(factory); + } + + // Create the setup: uses registered factory or falls back to DefaultSetup + static std::unique_ptr Create( + Input& input, Grid& grid, DataBlock& data, Output& output) { + auto& instance = GetInstance(); + if (instance.factory) { + return instance.factory(input, grid, data, output); + } + // Fallback to default (backward compatible with old Setup class) + return std::make_unique(input, grid, data, output); + } + + private: + SetupFactory factory; + + static SetupRegistry& GetInstance() { + static SetupRegistry instance; // Meyer's singleton + return instance; + } + + SetupRegistry() = default; + ~SetupRegistry() = default; + SetupRegistry(const SetupRegistry&) = delete; + SetupRegistry& operator=(const SetupRegistry&) = delete; +}; + +// Macro for one-line registration in user's setup.cpp +#define REGISTER_SETUP(SetupClass) \ + namespace { \ + struct SetupRegistrar { \ + SetupRegistrar() { \ + SetupRegistry::Register( \ + [](Input& input, Grid& grid, DataBlock& data, Output& output) { \ + return std::make_unique( \ + input, grid, data, output); \ + }); \ + } \ + } setup_registrar; \ + } // namespace + +#endif // SETUP_REGISTRY_HPP_ diff --git a/test/HD/MachReflection_newsetup/definitions.hpp b/test/HD/MachReflection_newsetup/definitions.hpp new file mode 100644 index 00000000..bfe30684 --- /dev/null +++ b/test/HD/MachReflection_newsetup/definitions.hpp @@ -0,0 +1,18 @@ +// *********************************************************************************** +// Idefix MHD astrophysical code +// +// Source file test/HD/MachReflection/definitions.hpp +// +// Last modified : 11/2020 +// +// Copyright(C) by : +// - Geoffroy Lesur (IPAG/UGA/CNRS - 2020) +// and other code contributors +// +// Licensed under CeCILL 2.1 License, see COPYING for more information +// *********************************************************************************** + +#define COMPONENTS 2 +#define DIMENSIONS 2 + +#define GEOMETRY CARTESIAN diff --git a/test/HD/MachReflection_newsetup/idefix-hll.ini b/test/HD/MachReflection_newsetup/idefix-hll.ini new file mode 100644 index 00000000..8572dded --- /dev/null +++ b/test/HD/MachReflection_newsetup/idefix-hll.ini @@ -0,0 +1,26 @@ +[Grid] +X1-grid 1 0.0 480 u 4.0 +X2-grid 1 0.0 120 u 1.0 +X3-grid 1 0.0 1 u 1.0 + +[TimeIntegrator] +CFL 0.8 +tstop 0.2 +first_dt 1.e-5 +nstages 2 + +[Hydro] +solver hll +gamma 1.4 + +[Boundary] +X1-beg userdef +X1-end outflow +X2-beg userdef +X2-end userdef +X3-beg outflow +X3-end outflow + +[Output] +vtk 0.2 +dmp 0.2 diff --git a/test/HD/MachReflection_newsetup/idefix-hllc.ini b/test/HD/MachReflection_newsetup/idefix-hllc.ini new file mode 100644 index 00000000..9dee90ee --- /dev/null +++ b/test/HD/MachReflection_newsetup/idefix-hllc.ini @@ -0,0 +1,26 @@ +[Grid] +X1-grid 1 0.0 480 u 4.0 +X2-grid 1 0.0 120 u 1.0 +X3-grid 1 0.0 1 u 1.0 + +[TimeIntegrator] +CFL 0.8 +tstop 0.2 +first_dt 1.e-5 +nstages 2 + +[Hydro] +solver hllc +gamma 1.4 + +[Boundary] +X1-beg userdef +X1-end outflow +X2-beg userdef +X2-end userdef +X3-beg outflow +X3-end outflow + +[Output] +vtk 0.2 +dmp 0.2 diff --git a/test/HD/MachReflection_newsetup/idefix-tvdlf.ini b/test/HD/MachReflection_newsetup/idefix-tvdlf.ini new file mode 100644 index 00000000..dfee317a --- /dev/null +++ b/test/HD/MachReflection_newsetup/idefix-tvdlf.ini @@ -0,0 +1,26 @@ +[Grid] +X1-grid 1 0.0 480 u 4.0 +X2-grid 1 0.0 120 u 1.0 +X3-grid 1 0.0 1 u 1.0 + +[TimeIntegrator] +CFL 0.8 +tstop 0.2 +first_dt 1.e-5 +nstages 2 + +[Hydro] +solver tvdlf +gamma 1.4 + +[Boundary] +X1-beg userdef +X1-end outflow +X2-beg userdef +X2-end userdef +X3-beg outflow +X3-end outflow + +[Output] +vtk 0.2 +dmp 0.2 diff --git a/test/HD/MachReflection_newsetup/idefix.ini b/test/HD/MachReflection_newsetup/idefix.ini new file mode 100644 index 00000000..caccba34 --- /dev/null +++ b/test/HD/MachReflection_newsetup/idefix.ini @@ -0,0 +1,26 @@ +[Grid] +X1-grid 1 0.0 480 u 4.0 +X2-grid 1 0.0 120 u 1.0 +X3-grid 1 0.0 1 u 1.0 + +[TimeIntegrator] +CFL 0.8 +tstop 0.2 +first_dt 1.e-5 +nstages 2 + +[Hydro] +solver roe +gamma 1.4 + +[Boundary] +X1-beg userdef +X1-end outflow +X2-beg userdef +X2-end userdef +X3-beg outflow +X3-end outflow + +[Output] +vtk 0.2 +dmp 0.2 diff --git a/test/HD/MachReflection_newsetup/setup.cpp b/test/HD/MachReflection_newsetup/setup.cpp new file mode 100644 index 00000000..a92564e6 --- /dev/null +++ b/test/HD/MachReflection_newsetup/setup.cpp @@ -0,0 +1,141 @@ +// *********************************************************************************** +// Idefix MHD astrophysical code +// +// Source file test/HD/MachReflection_newsetup/setup.cpp +// +// Last modified : 10/2026 +// +// Copyright(C) by : +// - Geoffroy Lesur (IPAG/UGA/CNRS - 2026) +// and other code contributors +// +// Licensed under CeCILL 2.1 License, see COPYING for more information +// *********************************************************************************** + +#include +#include "idefix.hpp" +#include "setup.hpp" +#include "setup_registry.hpp" + +// User's Setup class derived from SetupBase +class MachReflectionSetup : public SetupBase { +public: + using SetupBase::SetupBase; // Inherit constructor + + // ===== Member functions (previously global) ===== + + // User-defined boundary - now a member function! + void MyUserdefBoundary(Fluid* hydro, int dir, BoundarySide side, + const real t) { + const real alpha = 60./180.*M_PI; + + if ((dir==IDIR) && (side == left)) { + IdefixArray4D Vc = hydro->Vc; + int ighost = hydro->data->nghost[IDIR]; + hydro->boundary->BoundaryFor("UserDefBoundaryX1Beg", dir, side, + KOKKOS_LAMBDA (int k, int j, int i) { + + Vc(RHO,k,j,i) = 8.0; + Vc(VX1,k,j,i) = 8.25*sin(alpha); + Vc(VX2,k,j,i) = - 8.25*cos(alpha); + Vc(PRS,k,j,i) = 116.5; + }); + } + + if(dir==JDIR) { + IdefixArray4D Vc = hydro->Vc; + IdefixArray1D x1 = hydro->data->x[IDIR]; + IdefixArray1D x2 = hydro->data->x[JDIR]; + if(side == left) { + const int jend = hydro->data->beg[JDIR]; + hydro->boundary->BoundaryFor("UserDefBoundaryX2Beg", dir, side, + KOKKOS_LAMBDA (int k, int j, int i) { + if(x1(i) < 1.0/6.0) { + Vc(RHO,k,j,i) = 8.0; + Vc(VX1,k,j,i) = 8.25*sin(alpha); + Vc(VX2,k,j,i) = - 8.25*cos(alpha); + Vc(PRS,k,j,i) = 116.5; + } else { + Vc(RHO,k,j,i) = Vc(RHO,k,2*jend-j-1,i); + Vc(VX1,k,j,i) = Vc(VX1,k,2*jend-j-1,i); + Vc(VX2,k,j,i) = -Vc(VX2,k,2*jend-j-1,i); + Vc(PRS,k,j,i) = Vc(PRS,k,2*jend-j-1,i); + } + + }); + //return; + } else if(side==right) { + real xs = 10.0*t/sin(alpha) + 1.0/6.0 + 1.0/tan(alpha); + + hydro->boundary->BoundaryFor("UserDefBoundaryX2End", dir, side, + KOKKOS_LAMBDA (int k, int j, int i) { + if(x1(i) < xs) { + Vc(RHO,k,j,i) = 8.0; + Vc(VX1,k,j,i) = 8.25*sin(alpha); + Vc(VX2,k,j,i) = - 8.25*cos(alpha); + Vc(PRS,k,j,i) = 116.5; + } else { + Vc(RHO,k,j,i) = 1.4; + Vc(VX1,k,j,i) = 0.0; + Vc(VX2,k,j,i) = 0.0; + Vc(PRS,k,j,i) = 1.0; + } + + }); + } + } + + } + + // ===== Constructor - enroll member functions ===== + MachReflectionSetup(Input &input, Grid &grid, DataBlock &data, Output &output) + : SetupBase(input, grid, data, output) { + + // Enroll boundary function using lambda that captures 'this' + myData.hydro->EnrollUserDefBoundary( + [this](Fluid* hydro, int dir, BoundarySide side, + const real t) { + this->MyUserdefBoundary(hydro, dir, side, t); + } + ); + } + + // ===== InitFlow ===== + void InitFlow(DataBlock &data) override { + // Create a host copy + DataBlockHost d(data); + real alpha, xs, x1, x2; + + alpha = 1.0/3.0*M_PI; + + for(int k = 0; k < d.np_tot[KDIR] ; k++) { + for(int j = 0; j < d.np_tot[JDIR] ; j++) { + for(int i = 0; i < d.np_tot[IDIR] ; i++) { + x1 = d.x[IDIR](i); + x2 = d.x[JDIR](j); + xs = 1.0/6.0 + x2/tan(alpha); + if(x1>xs) { + d.Vc(RHO,k,j,i) = 1.4; + d.Vc(VX1,k,j,i) = 0.0; + d.Vc(VX2,k,j,i) = 0.0; + d.Vc(PRS,k,j,i) = 1.0; + } else { + d.Vc(RHO,k,j,i) = 8.0; + d.Vc(VX1,k,j,i) = 8.25*sin(alpha); + d.Vc(VX2,k,j,i) = -8.25*cos(alpha); + d.Vc(PRS,k,j,i) = 116.5; + } + } + } + } + + // Send it all, if needed + d.SyncToDevice(); + } + +private: + // Can add member variables here if needed +}; + +// ===== Registration ===== +REGISTER_SETUP(MachReflectionSetup) diff --git a/test/HD/MachReflection_newsetup/testme.json b/test/HD/MachReflection_newsetup/testme.json new file mode 100644 index 00000000..c701471f --- /dev/null +++ b/test/HD/MachReflection_newsetup/testme.json @@ -0,0 +1,17 @@ +{ + "namings": "ini,mpi", + "default": { + "dumpname": "dump.0001.dmp", + "noplot": true, + "reconstruction": 2, + "single": false, + "dec": [2, 2], + "tolerance": 0 + }, + "variants": [ + { + "ini": ["idefix.ini","idefix-hll.ini","idefix-hllc.ini","idefix-tvdlf.ini"], + "mpi": [false, true] + } + ] +}