From ab340ed54c34892826716dfa12d57958f76fb102 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Mon, 10 Feb 2025 17:55:33 +0100 Subject: [PATCH 01/17] implements a 1st order IMEX1 solver for polydisperse dust species --- src/dataBlock/evolveStage.cpp | 10 ++ src/fluid/drag.cpp | 194 +++++++++++++++++++++++ src/fluid/drag.hpp | 28 +++- src/fluid/evolveStage.hpp | 6 +- test/Dust/DustEnergy/idefix-implicit.ini | 36 +++++ test/Dust/DustEnergy/testme.py | 2 +- test/Dust/DustyWave/idefix-implicit.ini | 36 +++++ test/Dust/DustyWave/testme.py | 2 +- 8 files changed, 309 insertions(+), 5 deletions(-) create mode 100644 test/Dust/DustEnergy/idefix-implicit.ini create mode 100644 test/Dust/DustyWave/idefix-implicit.ini diff --git a/src/dataBlock/evolveStage.cpp b/src/dataBlock/evolveStage.cpp index 981c060a..c97c9f71 100644 --- a/src/dataBlock/evolveStage.cpp +++ b/src/dataBlock/evolveStage.cpp @@ -19,6 +19,16 @@ void DataBlock::EvolveStage() { for(int i = 0 ; i < dust.size() ; i++) { dust[i]->EvolveStage(this->t,this->dt); } + // Add implicit term for dust drag + if(dust[0]->drag->IsImplicit()) { + for(int i = 0 ; i < dust.size() ; i++) { + dust[i]->drag->AddImplicitBackReaction(this->dt,dust[0]->drag->implicitFactor); + } + dust[0]->drag->NormalizeImplicitBackReaction(this->dt); + for(int i = 0 ; i < dust.size() ; i++) { + dust[i]->drag->AddImplicitFluidMomentum(this->dt); + } + } } idfx::popRegion(); diff --git a/src/fluid/drag.cpp b/src/fluid/drag.cpp index d24e0bea..46542147 100644 --- a/src/fluid/drag.cpp +++ b/src/fluid/drag.cpp @@ -21,6 +21,10 @@ void Drag::AddDragForce(const real dt) { EquationOfState eos = *(this->eos); + if(implicit) { + IDEFIX_ERROR("Add DragForce should not be called when drag is implicit"); + } + auto userGammai = this->gammai; if(type == Type::Userdef) { if(userDrag != NULL) { @@ -78,9 +82,199 @@ void Drag::AddDragForce(const real dt) { }); idfx::popRegion(); } +// +// $ p_g^{(n+1)}=p_g^{(n)}+\sum_i\frac{\rho_g\gamma_i dt}{1+\rho_g \gamma_i dt}p_i^{(n)} $ +// We accumulate in the array "prefactor" +// $ \sum_i\frac{\rho_i\gamma_i dt}{1+\rho_g\gamma_i dt} +// + +void Drag::AddImplicitBackReaction(const real dt, IdefixArray3D preFactor) { + if(!feedback) { + // no feedback, no need for anything + return; + } + idfx::pushRegion("AddImplicitFluidMomentum"); + + auto UcGas = this->UcGas; + auto VcGas = this->VcGas; + auto UcDust = this->UcDust; + auto VcDust = this->VcDust; + + const Type type = this->type; + real dragCoeff = this->dragCoeff; + bool feedback = this->feedback; + + bool isFirst = this->instanceNumber == 0; + + EquationOfState eos = *(this->eos); + + if(!implicit) { + IDEFIX_ERROR("AddImplicitGasMomentum should not be called when drag is explicit"); + } + + auto userGammai = this->gammai; + if(type == Type::Userdef) { // If feedback enabled, + // the coefficients were already computed in the gas + // momentum pass + if(userDrag != NULL) { + idfx::pushRegion("Drag::UserDrag"); + userDrag(data, dragCoeff, userGammai); + idfx::popRegion(); + } else { + IDEFIX_ERROR("No User-defined drag function has been enrolled"); + } + } + // Compute a drag force fd = - gamma*rhod*rhog*(vd-vg) + // Where gamma is computed according to the choice of drag type + idefix_for("DragForce",0,data->np_tot[KDIR],0,data->np_tot[JDIR],0,data->np_tot[IDIR], + KOKKOS_LAMBDA (int k, int j, int i) { + real gamma; // The drag coefficient + if(type == Type::Gamma) { + gamma = dragCoeff; + + } else if(type == Type::Tau) { + // In this case, the coefficient is the stopping time (assumed constant) + gamma = 1/(dragCoeff*VcGas(RHO,k,j,i)); + } else if(type == Type::Size) { + real cs; + // Assume a fixed size, hence for both Epstein or Stokes, gamma~1/rho_g/cs + // Get the sound speed + #if HAVE_ENERGY == 1 + cs = std::sqrt(eos.GetGamma(VcGas(PRS,k,j,i),VcGas(RHO,k,j,i) + *VcGas(PRS,k,j,i)/VcGas(RHO,k,j,i))); + #else + cs = eos.GetWaveSpeed(k,j,i); + #endif + gamma = cs/dragCoeff; + } else if(type == Type::Userdef) { + gamma = userGammai(k,j,i); + } + + const real factor = VcDust(RHO,k,j,i)*gamma*dt/(1+VcGas(RHO,k,j,i)*gamma*dt); + if(isFirst) { + preFactor(k,j,i) = factor; + } else { + preFactor(k,j,i) += factor; + } + + for(int n = MX1 ; n < MX1+COMPONENTS ; n++) { + UcGas(n,k,j,i) += dt * gamma * VcGas(RHO,k,j,i) * UcDust(n,k,j,i) / + (1 + VcGas(RHO,k,j,i)*dt*gamma); + } + }); + idfx::popRegion(); +} + +void Drag::NormalizeImplicitBackReaction(const real dt) { + if(!feedback) { + // no feedback, no need for anything + return; + } + idfx::pushRegion("AddImplicitFluidMomentum"); + + auto UcGas = this->UcGas; + auto preFactor = this->implicitFactor; + + if(!implicit) { + IDEFIX_ERROR("AddImplicitGasMomentum should not be called when drag is explicit"); + } + + // Compute a drag force fd = - gamma*rhod*rhog*(vd-vg) + // Where gamma is computed according to the choice of drag type + idefix_for("DragForce",0,data->np_tot[KDIR],0,data->np_tot[JDIR],0,data->np_tot[IDIR], + KOKKOS_LAMBDA (int k, int j, int i) { + const real factor = 1+preFactor(k,j,i); + for(int n = MX1 ; n < MX1+COMPONENTS ; n++) { + UcGas(n,k,j,i) /= factor; + } + }); + idfx::popRegion(); +} + +void Drag::AddImplicitFluidMomentum(const real dt) { + idfx::pushRegion("AddImplicitFluidMomentum"); + + auto UcGas = this->UcGas; + auto VcGas = this->VcGas; + auto UcDust = this->UcDust; + auto VcDust = this->VcDust; + auto InvDt = this->InvDt; + + const Type type = this->type; + real dragCoeff = this->dragCoeff; + bool feedback = this->feedback; + + EquationOfState eos = *(this->eos); + + if(!implicit) { + IDEFIX_ERROR("AddImplicitGasMomentum should not be called when drag is explicit"); + } + + auto userGammai = this->gammai; + if(type == Type::Userdef && !feedback) { // If feedback enabled, + // the coefficients were already computed in the gas + // momentum pass + if(userDrag != NULL) { + idfx::pushRegion("Drag::UserDrag"); + userDrag(data, dragCoeff, userGammai); + idfx::popRegion(); + } else { + IDEFIX_ERROR("No User-defined drag function has been enrolled"); + } + } + // Compute a drag force fd = - gamma*rhod*rhog*(vd-vg) + // Where gamma is computed according to the choice of drag type + idefix_for("DragForce",0,data->np_tot[KDIR],0,data->np_tot[JDIR],0,data->np_tot[IDIR], + KOKKOS_LAMBDA (int k, int j, int i) { + real gamma; // The drag coefficient + if(type == Type::Gamma) { + gamma = dragCoeff; + + } else if(type == Type::Tau) { + // In this case, the coefficient is the stopping time (assumed constant) + gamma = 1/(dragCoeff*VcGas(RHO,k,j,i)); + } else if(type == Type::Size) { + real cs; + // Assume a fixed size, hence for both Epstein or Stokes, gamma~1/rho_g/cs + // Get the sound speed + #if HAVE_ENERGY == 1 + cs = std::sqrt(eos.GetGamma(VcGas(PRS,k,j,i),VcGas(RHO,k,j,i) + *VcGas(PRS,k,j,i)/VcGas(RHO,k,j,i))); + #else + cs = eos.GetWaveSpeed(k,j,i); + #endif + gamma = cs/dragCoeff; + } else if(type == Type::Userdef) { + gamma = userGammai(k,j,i); + } + + for(int n = MX1 ; n < MX1+COMPONENTS ; n++) { + real oldUc = UcDust(n,k,j,i); + UcDust(n,k,j,i) = (oldUc + dt * gamma * VcDust(RHO,k,j,i) * UcGas(n,k,j,i)) / + (1 + VcGas(RHO,k,j,i)*dt*gamma); + + #if HAVE_ENERGY == 1 + real dp = UcDust(n,k,j,i) - oldUc; + + // We add back the energy dissipated for the dust which is not accounted for + // (since there is no energy equation for dust grains) + + // TODO(GL): this should be disabled in the case of a true multifluid system where + // both fluids have a proper energy equation + UcGas(ENG,k,j,i) -= dp*VcDust(n,k,j,i); + #endif + } + }); + idfx::popRegion(); +} void Drag::ShowConfig() { idfx::cout << "Drag: Using "; + if(implicit) { + idfx::cout << " IMPLICIT "; + } else { + idfx::cout << " EXPLICIT "; + } switch(type) { case Type::Gamma: idfx::cout << "constant gamma"; diff --git a/src/fluid/drag.hpp b/src/fluid/drag.hpp index d5f7a5a7..c676a88e 100644 --- a/src/fluid/drag.hpp +++ b/src/fluid/drag.hpp @@ -24,7 +24,16 @@ class Drag { Drag(Input &, Fluid *); void ShowConfig(); // print configuration void AddDragForce(const real); + + ////////////////////////// + // Implicit functions + void AddImplicitBackReaction(const real, IdefixArray3D); // Add the back reaction + void NormalizeImplicitBackReaction(const real); // Normalize the implicit back reaction + void AddImplicitFluidMomentum(const real); // Add the implicit drag force on dust grains + ///////////////////////// + void EnrollUserDrag(UserDefDragFunc); // User defined drag function enrollment + bool IsImplicit() const { return implicit; } // Check if the drag is implicit IdefixArray4D UcDust; // Dust conservative quantities IdefixArray4D UcGas; // Gas conservative quantities @@ -32,12 +41,17 @@ class Drag { IdefixArray4D VcGas; // Gas primitive quantities IdefixArray3D InvDt; // The InvDt of current dust specie IdefixArray3D gammai; // the drag coefficient (only used for user-defined dust grains) + IdefixArray3D implicitFactor; // The prefactor used by the implicit timestepping + Type type; private: DataBlock* data; real dragCoeff; bool feedback{false}; + bool implicit{false}; + int instanceNumber; + UserDefDragFunc userDrag{NULL}; @@ -92,11 +106,21 @@ Drag::Drag(Input &input, Fluid *hydroin): IDEFIX_ERROR(msg); } // Fetch the drag coefficient for the current specie. - const int n = hydroin->instanceNumber; - this->dragCoeff = input.Get(BlockName,"drag",n+1); + this->instanceNumber = hydroin->instanceNumber; + this->dragCoeff = input.Get(BlockName,"drag",instanceNumber+1); // 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); + + if(implicit && instanceNumber == 0) { + this->implicitFactor = IdefixArray3D("ImplicitFactor", + data->np_tot[KDIR], + data->np_tot[JDIR], + data->np_tot[IDIR]); + } + } else { IDEFIX_ERROR("A [Drag] block is required in your input file to define the drag force."); } diff --git a/src/fluid/evolveStage.hpp b/src/fluid/evolveStage.hpp index 06754d0b..31dae5b7 100644 --- a/src/fluid/evolveStage.hpp +++ b/src/fluid/evolveStage.hpp @@ -61,7 +61,11 @@ void Fluid::EvolveStage(const real t, const real dt) { if(haveSourceTerms) AddSourceTerms(t, dt); // Step 5: add drag when needed - if(haveDrag) drag->AddDragForce(dt); + if(haveDrag) { + if(!drag->IsImplicit()) { + drag->AddDragForce(dt); + } + } if constexpr(Phys::mhd) { #if DIMENSIONS >= 2 diff --git a/test/Dust/DustEnergy/idefix-implicit.ini b/test/Dust/DustEnergy/idefix-implicit.ini new file mode 100644 index 00000000..9ba3c781 --- /dev/null +++ b/test/Dust/DustEnergy/idefix-implicit.ini @@ -0,0 +1,36 @@ +# This test checks that the total energy (thermal+dust kinetic+gas kinetic) +# is effectively conserved when drag is present + +[Grid] +X1-grid 1 0.0 500 u 1.0 +X2-grid 1 0.0 1 u 1.0 +X3-grid 1 0.0 1 u 1.0 + +[TimeIntegrator] +CFL 0.8 +tstop 1.0 +first_dt 1.e-4 +nstages 2 + +[Hydro] +solver hllc +gamma 1.4 + +[Dust] +nSpecies 1 +drag tau 1.0 +drag_feedback yes +drag_implicit yes + +[Boundary] +X1-beg periodic +X1-end periodic +X2-beg outflow +X2-end outflow +X3-beg outflow +X3-end outflow + +[Output] +dmp 1.0 +analysis 0.01 +log 1000 diff --git a/test/Dust/DustEnergy/testme.py b/test/Dust/DustEnergy/testme.py index f3b413e8..7187d299 100755 --- a/test/Dust/DustEnergy/testme.py +++ b/test/Dust/DustEnergy/testme.py @@ -15,7 +15,7 @@ def testMe(test): test.configure() test.compile() - inifiles=["idefix.ini"] + inifiles=["idefix.ini","idefix-implicit.ini"] # loop on all the ini files for this test for ini in inifiles: diff --git a/test/Dust/DustyWave/idefix-implicit.ini b/test/Dust/DustyWave/idefix-implicit.ini new file mode 100644 index 00000000..05da6ccc --- /dev/null +++ b/test/Dust/DustyWave/idefix-implicit.ini @@ -0,0 +1,36 @@ +# This test checks the dissipation of a sound wave by a dust grains +# partially coupled to the gas (Riols & Lesur 2018, appendix A) + +[Grid] +X1-grid 1 0.0 500 u 1.0 +X2-grid 1 0.0 1 u 1.0 +X3-grid 1 0.0 1 u 1.0 + +[TimeIntegrator] +CFL 0.8 +tstop 10.0 +first_dt 1.e-4 +nstages 2 + +[Hydro] +solver hllc +csiso constant 1.0 + +[Dust] +nSpecies 1 +drag tau 1.0e-3 +drag_feedback yes +drag_implicit yes + +[Boundary] +X1-beg periodic +X1-end periodic +X2-beg outflow +X2-end outflow +X3-beg outflow +X3-end outflow + +[Output] +dmp 10.0 +analysis 0.01 +log 1000 diff --git a/test/Dust/DustyWave/testme.py b/test/Dust/DustyWave/testme.py index b6a75ce3..2b03d163 100755 --- a/test/Dust/DustyWave/testme.py +++ b/test/Dust/DustyWave/testme.py @@ -15,7 +15,7 @@ def testMe(test): test.configure() test.compile() - inifiles=["idefix.ini"] + inifiles=["idefix.ini","idefix-implicit.ini"] # loop on all the ini files for this test for ini in inifiles: From 22eaf5b872f26e8ec2f4a329497e2ef4805ed315 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Tue, 11 Feb 2025 10:33:14 +0100 Subject: [PATCH 02/17] refactor dust drag law --- src/fluid/drag.cpp | 188 ++++++++++-------------- src/fluid/drag.hpp | 91 +++++++----- test/Dust/DustyWave/idefix-implicit.ini | 2 +- 3 files changed, 127 insertions(+), 154 deletions(-) diff --git a/src/fluid/drag.cpp b/src/fluid/drag.cpp index 46542147..7739d84d 100644 --- a/src/fluid/drag.cpp +++ b/src/fluid/drag.cpp @@ -5,7 +5,12 @@ // Licensed under CeCILL 2.1 License, see COPYING for more information // *********************************************************************************** #include "drag.hpp" +#include + #include "physics.hpp" + + + void Drag::AddDragForce(const real dt) { idfx::pushRegion("Drag::AddDragForce"); @@ -15,51 +20,20 @@ void Drag::AddDragForce(const real dt) { auto VcDust = this->VcDust; auto InvDt = this->InvDt; - const Type type = this->type; - real dragCoeff = this->dragCoeff; bool feedback = this->feedback; - EquationOfState eos = *(this->eos); - if(implicit) { IDEFIX_ERROR("Add DragForce should not be called when drag is implicit"); } - auto userGammai = this->gammai; - if(type == Type::Userdef) { - if(userDrag != NULL) { - idfx::pushRegion("Drag::UserDrag"); - userDrag(data, dragCoeff, userGammai); - idfx::popRegion(); - } else { - IDEFIX_ERROR("No User-defined drag function has been enrolled"); - } - } + auto gammaDrag = this->gammaDrag; + gammaDrag.RefreshUserDrag(data); + // Compute a drag force fd = - gamma*rhod*rhog*(vd-vg) // Where gamma is computed according to the choice of drag type idefix_for("DragForce",0,data->np_tot[KDIR],0,data->np_tot[JDIR],0,data->np_tot[IDIR], KOKKOS_LAMBDA (int k, int j, int i) { - real gamma; // The drag coefficient - if(type == Type::Gamma) { - gamma = dragCoeff; - - } else if(type == Type::Tau) { - // In this case, the coefficient is the stopping time (assumed constant) - gamma = 1/(dragCoeff*VcGas(RHO,k,j,i)); - } else if(type == Type::Size) { - real cs; - // Assume a fixed size, hence for both Epstein or Stokes, gamma~1/rho_g/cs - // Get the sound speed - #if HAVE_ENERGY == 1 - cs = std::sqrt(eos.GetGamma(VcGas(PRS,k,j,i),VcGas(RHO,k,j,i) - *VcGas(PRS,k,j,i)/VcGas(RHO,k,j,i))); - #else - cs = eos.GetWaveSpeed(k,j,i); - #endif - gamma = cs/dragCoeff; - } else if(type == Type::Userdef) { - gamma = userGammai(k,j,i); - } + real gamma = gammaDrag.GetGamma(k,j,i); // The drag coefficient real dp = dt * gamma * VcDust(RHO,k,j,i) * VcGas(RHO,k,j,i); for(int n = MX1 ; n < MX1+COMPONENTS ; n++) { @@ -97,58 +71,26 @@ void Drag::AddImplicitBackReaction(const real dt, IdefixArray3D preFactor) auto UcGas = this->UcGas; auto VcGas = this->VcGas; - auto UcDust = this->UcDust; auto VcDust = this->VcDust; - const Type type = this->type; - real dragCoeff = this->dragCoeff; bool feedback = this->feedback; bool isFirst = this->instanceNumber == 0; - EquationOfState eos = *(this->eos); if(!implicit) { IDEFIX_ERROR("AddImplicitGasMomentum should not be called when drag is explicit"); } - auto userGammai = this->gammai; - if(type == Type::Userdef) { // If feedback enabled, - // the coefficients were already computed in the gas - // momentum pass - if(userDrag != NULL) { - idfx::pushRegion("Drag::UserDrag"); - userDrag(data, dragCoeff, userGammai); - idfx::popRegion(); - } else { - IDEFIX_ERROR("No User-defined drag function has been enrolled"); - } - } + auto gammaDrag = this->gammaDrag; + gammaDrag.RefreshUserDrag(data); + // Compute a drag force fd = - gamma*rhod*rhog*(vd-vg) // Where gamma is computed according to the choice of drag type idefix_for("DragForce",0,data->np_tot[KDIR],0,data->np_tot[JDIR],0,data->np_tot[IDIR], KOKKOS_LAMBDA (int k, int j, int i) { - real gamma; // The drag coefficient - if(type == Type::Gamma) { - gamma = dragCoeff; - - } else if(type == Type::Tau) { - // In this case, the coefficient is the stopping time (assumed constant) - gamma = 1/(dragCoeff*VcGas(RHO,k,j,i)); - } else if(type == Type::Size) { - real cs; - // Assume a fixed size, hence for both Epstein or Stokes, gamma~1/rho_g/cs - // Get the sound speed - #if HAVE_ENERGY == 1 - cs = std::sqrt(eos.GetGamma(VcGas(PRS,k,j,i),VcGas(RHO,k,j,i) - *VcGas(PRS,k,j,i)/VcGas(RHO,k,j,i))); - #else - cs = eos.GetWaveSpeed(k,j,i); - #endif - gamma = cs/dragCoeff; - } else if(type == Type::Userdef) { - gamma = userGammai(k,j,i); - } + real gamma = gammaDrag.GetGamma(k,j,i); // The drag coefficient + const real factor = VcDust(RHO,k,j,i)*gamma*dt/(1+VcGas(RHO,k,j,i)*gamma*dt); if(isFirst) { @@ -200,53 +142,20 @@ void Drag::AddImplicitFluidMomentum(const real dt) { auto VcDust = this->VcDust; auto InvDt = this->InvDt; - const Type type = this->type; - real dragCoeff = this->dragCoeff; bool feedback = this->feedback; - EquationOfState eos = *(this->eos); - if(!implicit) { IDEFIX_ERROR("AddImplicitGasMomentum should not be called when drag is explicit"); } - auto userGammai = this->gammai; - if(type == Type::Userdef && !feedback) { // If feedback enabled, - // the coefficients were already computed in the gas - // momentum pass - if(userDrag != NULL) { - idfx::pushRegion("Drag::UserDrag"); - userDrag(data, dragCoeff, userGammai); - idfx::popRegion(); - } else { - IDEFIX_ERROR("No User-defined drag function has been enrolled"); - } - } + auto gammaDrag = this->gammaDrag; + if(!feedback) gammaDrag.RefreshUserDrag(data); + // Compute a drag force fd = - gamma*rhod*rhog*(vd-vg) // Where gamma is computed according to the choice of drag type idefix_for("DragForce",0,data->np_tot[KDIR],0,data->np_tot[JDIR],0,data->np_tot[IDIR], KOKKOS_LAMBDA (int k, int j, int i) { - real gamma; // The drag coefficient - if(type == Type::Gamma) { - gamma = dragCoeff; - - } else if(type == Type::Tau) { - // In this case, the coefficient is the stopping time (assumed constant) - gamma = 1/(dragCoeff*VcGas(RHO,k,j,i)); - } else if(type == Type::Size) { - real cs; - // Assume a fixed size, hence for both Epstein or Stokes, gamma~1/rho_g/cs - // Get the sound speed - #if HAVE_ENERGY == 1 - cs = std::sqrt(eos.GetGamma(VcGas(PRS,k,j,i),VcGas(RHO,k,j,i) - *VcGas(PRS,k,j,i)/VcGas(RHO,k,j,i))); - #else - cs = eos.GetWaveSpeed(k,j,i); - #endif - gamma = cs/dragCoeff; - } else if(type == Type::Userdef) { - gamma = userGammai(k,j,i); - } + real gamma = gammaDrag.GetGamma(k,j,i); // The drag coefficient for(int n = MX1 ; n < MX1+COMPONENTS ; n++) { real oldUc = UcDust(n,k,j,i); @@ -275,17 +184,17 @@ void Drag::ShowConfig() { } else { idfx::cout << " EXPLICIT "; } - switch(type) { - case Type::Gamma: + switch(gammaDrag.type) { + case GammaDrag::Type::Gamma: idfx::cout << "constant gamma"; break; - case Type::Tau: + case GammaDrag::Type::Tau: idfx::cout << "constant stopping time"; break; - case Type::Size: + case GammaDrag::Type::Size: idfx::cout << "constant dust size"; break; - case Type::Userdef: + case GammaDrag::Type::Userdef: idfx::cout << "user-defined"; break; } @@ -299,8 +208,59 @@ 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, 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) { + std::string dragType = input.Get(BlockName,"drag",0); + if(dragType.compare("gamma") == 0) { + this->type = Type::Gamma; + } else if(dragType.compare("tau") == 0) { + this->type = Type::Tau; + this->VcGas = data->hydro->Vc; + } else if(dragType.compare("size") == 0) { + this->type = Type::Size; + this->eos = *(data->hydro->eos.get()); + this->VcGas = data->hydro->Vc; + } else if(dragType.compare("userdef") == 0) { + this->type = Type::Userdef; + this->gammai = IdefixArray3D("UserDrag", + data->np_tot[KDIR], + data->np_tot[JDIR], + data->np_tot[IDIR]); + } else { + std::stringstream msg; + msg << "Unknown drag type \"" << dragType + << "\" in your input file." << std::endl + << "Allowed values are: gamma, tau, epstein, stokes, userdef." << std::endl; + + IDEFIX_ERROR(msg); + } + } + // Fetch the drag coefficient for the current specie. + this->dragCoeff = input.Get(BlockName,"drag",instanceNumber+1); +} diff --git a/src/fluid/drag.hpp b/src/fluid/drag.hpp index c676a88e..26332eff 100644 --- a/src/fluid/drag.hpp +++ b/src/fluid/drag.hpp @@ -15,10 +15,52 @@ using UserDefDragFunc = void (*) (DataBlock *, real beta, IdefixArray3D &gammai); +/* A class that holds the kind of drag law we intend to use */ +class GammaDrag { + public: + 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 + if(type == Type::Gamma) { + gamma = dragCoeff; + + } else if(type == Type::Tau) { + // In this case, the coefficient is the stopping time (assumed constant) + gamma = 1/(dragCoeff*VcGas(RHO,k,j,i)); + } else if(type == Type::Size) { + real cs; + // Assume a fixed size, hence for both Epstein or Stokes, gamma~1/rho_g/cs + // Get the sound speed + #if HAVE_ENERGY == 1 + cs = std::sqrt(eos.GetGamma(VcGas(PRS,k,j,i),VcGas(RHO,k,j,i) + *VcGas(PRS,k,j,i)/VcGas(RHO,k,j,i))); + #else + cs = eos.GetWaveSpeed(k,j,i); + #endif + gamma = cs/dragCoeff; + } else if(type == Type::Userdef) { + gamma = gammai(k,j,i); + } + return gamma; + } + + Type type; + real dragCoeff; + EquationOfState eos; + IdefixArray3D gammai; + IdefixArray4D VcGas; + + UserDefDragFunc userDrag{NULL}; +}; + class Drag { public: - enum class Type{Gamma, Tau, Size, Userdef}; // Different types of implementation for the drag force. template Drag(Input &, Fluid *); @@ -40,24 +82,20 @@ class Drag { IdefixArray4D VcDust; // Gas primitive quantities IdefixArray4D VcGas; // Gas primitive quantities IdefixArray3D InvDt; // The InvDt of current dust specie - IdefixArray3D gammai; // the drag coefficient (only used for user-defined dust grains) IdefixArray3D implicitFactor; // The prefactor used by the implicit timestepping - Type type; + GammaDrag gammaDrag; // The drag law private: DataBlock* data; - real dragCoeff; bool feedback{false}; bool implicit{false}; int instanceNumber; +}; + - UserDefDragFunc userDrag{NULL}; - // Sound speed computation - EquationOfState *eos; -}; #include "fluid.hpp" @@ -72,47 +110,22 @@ Drag::Drag(Input &input, Fluid *hydroin): // Save the parent hydro object this->data = hydroin->data; - // We use the EOS of the basic fluid instance, as there is no EOS for dust! - this->eos = hydroin->data->hydro->eos.get(); // Check in which block we should fetch our information - std::string BlockName; + std::string blockName; if(Phys::dust) { - BlockName = "Dust"; + blockName = "Dust"; } else { IDEFIX_ERROR("Drag is currently implemented only for dusty fuids"); } - if(input.CheckEntry(BlockName,"drag")>=0) { - std::string dragType = input.Get(BlockName,"drag",0); - if(dragType.compare("gamma") == 0) { - this->type = Type::Gamma; - } else if(dragType.compare("tau") == 0) { - this->type = Type::Tau; - } else if(dragType.compare("size") == 0) { - this->type = Type::Size; - } else if(dragType.compare("userdef") == 0) { - this->type = Type::Userdef; - this->gammai = IdefixArray3D("UserDrag", - data->np_tot[KDIR], - data->np_tot[JDIR], - data->np_tot[IDIR]); - } else { - std::stringstream msg; - msg << "Unknown drag type \"" << dragType - << "\" in your input file." << std::endl - << "Allowed values are: gamma, tau, epstein, stokes, userdef." << std::endl; - - IDEFIX_ERROR(msg); - } - // Fetch the drag coefficient for the current specie. + if(input.CheckEntry(blockName,"drag")>=0) { this->instanceNumber = hydroin->instanceNumber; - this->dragCoeff = input.Get(BlockName,"drag",instanceNumber+1); + this->gammaDrag = GammaDrag(input,blockName,instanceNumber,data); // 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); + this->feedback = input.GetOrSet(blockName,"drag_feedback",0,true); + this->implicit = input.GetOrSet(blockName,"drag_implicit",0,false); if(implicit && instanceNumber == 0) { this->implicitFactor = IdefixArray3D("ImplicitFactor", diff --git a/test/Dust/DustyWave/idefix-implicit.ini b/test/Dust/DustyWave/idefix-implicit.ini index 05da6ccc..726585e0 100644 --- a/test/Dust/DustyWave/idefix-implicit.ini +++ b/test/Dust/DustyWave/idefix-implicit.ini @@ -18,7 +18,7 @@ csiso constant 1.0 [Dust] nSpecies 1 -drag tau 1.0e-3 +drag tau 1.0 drag_feedback yes drag_implicit yes From cbcc63fffd13e8f3ec7cbc58bcd88a2ec0112e31 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Tue, 11 Feb 2025 10:35:35 +0100 Subject: [PATCH 03/17] clean up --- src/fluid/drag.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/fluid/drag.cpp b/src/fluid/drag.cpp index 7739d84d..791cdf88 100644 --- a/src/fluid/drag.cpp +++ b/src/fluid/drag.cpp @@ -6,11 +6,8 @@ // *********************************************************************************** #include "drag.hpp" #include - #include "physics.hpp" - - void Drag::AddDragForce(const real dt) { idfx::pushRegion("Drag::AddDragForce"); From 05ca488b9b069397532f5a9f3970b5a45c0d29d6 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Tue, 11 Feb 2025 14:11:04 +0100 Subject: [PATCH 04/17] change signature of userdef drag function --- src/fluid/drag.cpp | 3 ++- src/fluid/drag.hpp | 6 ++---- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/src/fluid/drag.cpp b/src/fluid/drag.cpp index 791cdf88..bbbf11b2 100644 --- a/src/fluid/drag.cpp +++ b/src/fluid/drag.cpp @@ -216,7 +216,7 @@ void GammaDrag::RefreshUserDrag(DataBlock *data) { if(type == Type::Userdef) { if(userDrag != NULL) { idfx::pushRegion("GammaDrag::UserDrag"); - userDrag(data, dragCoeff, gammai); + userDrag(data, instanceNumber, dragCoeff, gammai); idfx::popRegion(); } else { IDEFIX_ERROR("No User-defined drag function has been enrolled"); @@ -260,4 +260,5 @@ GammaDrag::GammaDrag(Input &input, std::string BlockName, int instanceNumber, Da } // Fetch the drag coefficient for the current specie. this->dragCoeff = input.Get(BlockName,"drag",instanceNumber+1); + this->instanceNumber = instanceNumber; } diff --git a/src/fluid/drag.hpp b/src/fluid/drag.hpp index 26332eff..29eca309 100644 --- a/src/fluid/drag.hpp +++ b/src/fluid/drag.hpp @@ -13,7 +13,7 @@ #include "fluid_defs.hpp" #include "eos.hpp" -using UserDefDragFunc = void (*) (DataBlock *, real beta, IdefixArray3D &gammai); +using UserDefDragFunc = void (*) (DataBlock *, int n, real beta, IdefixArray3D &gammai); /* A class that holds the kind of drag law we intend to use */ class GammaDrag { @@ -55,6 +55,7 @@ class GammaDrag { IdefixArray3D gammai; IdefixArray4D VcGas; + int instanceNumber; UserDefDragFunc userDrag{NULL}; }; @@ -94,9 +95,6 @@ class Drag { }; - - - #include "fluid.hpp" template From 6fbd722439179c67bf22518098d1b210ebede845 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Tue, 11 Feb 2025 14:11:57 +0100 Subject: [PATCH 05/17] add the dustyShock test from Benitez-Llambay 2019 --- test/Dust/DustyShock/definitions.hpp | 5 ++ test/Dust/DustyShock/idefix-implicit.ini | 36 +++++++++ test/Dust/DustyShock/idefix.ini | 35 ++++++++ test/Dust/DustyShock/python/testidefix.py | 97 ++++++++++++++++++++++ test/Dust/DustyShock/setup.cpp | 99 +++++++++++++++++++++++ test/Dust/DustyShock/testme.py | 39 +++++++++ 6 files changed, 311 insertions(+) create mode 100644 test/Dust/DustyShock/definitions.hpp create mode 100644 test/Dust/DustyShock/idefix-implicit.ini create mode 100644 test/Dust/DustyShock/idefix.ini create mode 100755 test/Dust/DustyShock/python/testidefix.py create mode 100644 test/Dust/DustyShock/setup.cpp create mode 100755 test/Dust/DustyShock/testme.py diff --git a/test/Dust/DustyShock/definitions.hpp b/test/Dust/DustyShock/definitions.hpp new file mode 100644 index 00000000..cdbe23d6 --- /dev/null +++ b/test/Dust/DustyShock/definitions.hpp @@ -0,0 +1,5 @@ +#define COMPONENTS 1 +#define DIMENSIONS 1 +#define ISOTHERMAL + +#define GEOMETRY CARTESIAN diff --git a/test/Dust/DustyShock/idefix-implicit.ini b/test/Dust/DustyShock/idefix-implicit.ini new file mode 100644 index 00000000..e83d5bf4 --- /dev/null +++ b/test/Dust/DustyShock/idefix-implicit.ini @@ -0,0 +1,36 @@ +# This test checks the behaviour of a dust sound shock +# following the 4 fluids test of Benitez-Llambay+ 2019 + +[Grid] +X1-grid 1 0.0 400 u 40.0 +X2-grid 1 0.0 1 u 1.0 +X3-grid 1 0.0 1 u 1.0 + +[TimeIntegrator] +CFL 0.8 +tstop 500.0 +first_dt 1.e-4 +nstages 2 + +[Hydro] +solver hllc +csiso constant 1.0 + +[Dust] +nSpecies 3 +drag userdef 1.0 3.0 5.0 +drag_feedback yes +drag_implicit yes + +[Boundary] +X1-beg userdef +X1-end userdef +X2-beg outflow +X2-end outflow +X3-beg outflow +X3-end outflow + +[Output] +# dmp 10.0 +vtk 500.0 +log 1000 diff --git a/test/Dust/DustyShock/idefix.ini b/test/Dust/DustyShock/idefix.ini new file mode 100644 index 00000000..1b9c9ea2 --- /dev/null +++ b/test/Dust/DustyShock/idefix.ini @@ -0,0 +1,35 @@ +# This test checks the behaviour of a dust sound shock +# following the 4 fluids test of Benitez-Llambay+ 2019 + +[Grid] +X1-grid 1 0.0 400 u 40.0 +X2-grid 1 0.0 1 u 1.0 +X3-grid 1 0.0 1 u 1.0 + +[TimeIntegrator] +CFL 0.8 +tstop 500.0 +first_dt 1.e-4 +nstages 2 + +[Hydro] +solver hllc +csiso constant 1.0 + +[Dust] +nSpecies 3 +drag userdef 1.0 3.0 5.0 +drag_feedback yes + +[Boundary] +X1-beg userdef +X1-end userdef +X2-beg outflow +X2-end outflow +X3-beg outflow +X3-end outflow + +[Output] +# dmp 10.0 +vtk 500.0 +log 1000 diff --git a/test/Dust/DustyShock/python/testidefix.py b/test/Dust/DustyShock/python/testidefix.py new file mode 100755 index 00000000..8fca15a3 --- /dev/null +++ b/test/Dust/DustyShock/python/testidefix.py @@ -0,0 +1,97 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Thu Mar 5 11:29:41 2020 + +@author: glesur +""" + +import os +import sys +sys.path.append(os.getenv("IDEFIX_DIR")) +from pytools.vtk_io import readVTK +import argparse +import numpy as np +import matplotlib.pyplot as plt +from scipy.integrate import solve_ivp + +parser = argparse.ArgumentParser() +parser.add_argument("-noplot", + default=False, + help="disable plotting", + action="store_true") + + +args, unknown=parser.parse_known_args() + +V=readVTK('../data.0001.vtk', geometry='cartesian') + +# Find where the shock starts +istart = np.argwhere(V.data['RHO'][:,0,0]>2.0)[0][0]-1 + +# Compute analytical solution +M=2 # Mach number +K=np.asarray([1,3,5])/M # Drag coefficients + +# (56) from Benitez-Llambay 2018 +def get_wg(Mach, wi): + coeff=[1,np.sum(wi-1)-1/Mach**2-1,1/Mach**2] + solution=np.roots(coeff) + return(np.min(solution)) + +# (55) from Benitez-Llambay 2018 +def rhs(t,wi,Mach,K): + wg = get_wg(Mach,wi) + return K*(wg-wi) + +wi0=np.asarray([1.0,1.0,1.0]) +x=V.x[istart:] +arguments=M,K + +sol = solve_ivp(rhs, t_span=[x[0],x[-1]], y0=wi0, t_eval=x,args=arguments) + +# compute gas velocity +wg=np.zeros(len(x)) +for i in range(0,len(wg)): + wg[i]=get_wg(M,sol.y[:,i]) + + +if(not args.noplot): + plt.figure(1) + plt.plot(V.x,V.data['RHO'][:,0,0],'o',markersize=4,color='tab:purple') + plt.plot(x,1/wg,color='tab:purple') + plt.plot(V.x,V.data['Dust0_RHO'][:,0,0],'o',markersize=4,color='tab:orange') + plt.plot(x,1/sol.y[0,:],color='tab:orange') + plt.plot(V.x,V.data['Dust1_RHO'][:,0,0],'o',markersize=4,color='tab:blue') + plt.plot(x,1/sol.y[1,:],color='tab:blue') + plt.plot(V.x,V.data['Dust2_RHO'][:,0,0],'o',markersize=4,color='tab:green') + plt.plot(x,1/sol.y[2,:],color='tab:green') + plt.xlim([0,10]) + plt.title('Density') + + plt.figure(2) + plt.plot(V.x,V.data['VX1'][:,0,0],'o',markersize=4,color='tab:purple') + plt.plot(x,M*wg,color='tab:purple') + plt.plot(V.x[::4],V.data['Dust0_VX1'][::4,0,0],'o',markersize=4,color='tab:orange') + plt.plot(x,M*sol.y[0,:],color='tab:orange') + plt.plot(V.x[::4],V.data['Dust1_VX1'][::4,0,0],'o',markersize=4,color='tab:blue') + plt.plot(x,M*sol.y[1,:],color='tab:blue') + plt.plot(V.x[::4],V.data['Dust2_VX1'][::4,0,0],'o',markersize=4,color='tab:green') + plt.plot(x,M*sol.y[2,:],color='tab:green') + + plt.xlim([0,10]) + plt.title('Velocity') + + plt.ioff() + plt.show() + +# Compute the error on the dust densities + +error=np.mean(np.fabs(V.data['Dust0_RHO'][istart:,0,0]-1/sol.y[0,:]) + np.fabs(V.data['Dust1_RHO'][istart:,0,0]-1/sol.y[1,:]) + np.fabs(V.data['Dust2_RHO'][istart:,0,0]-1/sol.y[2,:])) / 3 +print("Error=%e"%error) +if error<3.6e-2: + print("SUCCESS!") + sys.exit(0) +else: + print("FAILURE!") + sys.exit(1) diff --git a/test/Dust/DustyShock/setup.cpp b/test/Dust/DustyShock/setup.cpp new file mode 100644 index 00000000..f8632248 --- /dev/null +++ b/test/Dust/DustyShock/setup.cpp @@ -0,0 +1,99 @@ +#include "idefix.hpp" +#include "setup.hpp" + +#define FILENAME "timevol.dat" + +void MyDrag(DataBlock *data, int n, real beta, IdefixArray3D &gamma) { + // Compute the drag coefficient gamma from the input beta + auto VcGas = data->hydro->Vc; + auto VcDust = data->dust[n]->Vc; + + idefix_for("MyDrag",0,data->np_tot[KDIR],0,data->np_tot[JDIR],0,data->np_tot[IDIR], + KOKKOS_LAMBDA (int k, int j, int i) { + gamma(k,j,i) = beta/VcGas(RHO,k,j,i)/VcDust(RHO,k,j,i); + }); +} + + +void ApplyBoundary(DataBlock *data, IdefixArray4D Vc, int dir, BoundarySide side) { + if(dir==IDIR) { + int iref,ibeg,iend; + if(side == left) { + iref = data->beg[IDIR]; + ibeg = 0; + iend = data->beg[IDIR]; + } else { + iref =data->end[IDIR] - 1; + ibeg=data->end[IDIR]; + iend=data->np_tot[IDIR]; + } + int nvar = Vc.extent(0); + idefix_for("UserDefBoundary",0,data->np_tot[KDIR],0,data->np_tot[JDIR],ibeg,iend, + KOKKOS_LAMBDA (int k, int j, int i) { + for(int n = 0 ; n < nvar ; n++) { + Vc(n,k,j,i) = Vc(n,k,j,iref ); + } + }); + } +} + +void UserdefBoundary(Hydro *hydro, int dir, BoundarySide side, real t) { + ApplyBoundary(hydro->data, hydro->Vc, dir, side); +} + +void UserdefBoundaryDust(Fluid *dust, int dir, BoundarySide side, real t) { + ApplyBoundary(dust->data, dust->Vc, dir, side); +} + + +// Initialisation routine. Can be used to allocate +// Arrays or variables which are used later on +Setup::Setup(Input &input, Grid &grid, DataBlock &data, Output &output) { + + data.hydro->EnrollUserDefBoundary(&UserdefBoundary); + if(data.haveDust) { + int nSpecies = data.dust.size(); + for(int n = 0 ; n < nSpecies ; n++) { + data.dust[n]->EnrollUserDefBoundary(&UserdefBoundaryDust); + data.dust[n]->drag->EnrollUserDrag(&MyDrag); + } + } + +} + +// This routine initialize the flow +// Note that data is on the device. +// One can therefore define locally +// a datahost and sync it, if needed +void Setup::InitFlow(DataBlock &data) { + // Create a host copy + DataBlockHost d(data); + + const real x0 = 4.0; + + 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++) { + real x = d.x[IDIR](i); + + d.Vc(RHO,k,j,i) = (x < x0) ? 1.0 : 16.0; + d.Vc(VX1,k,j,i) = (x < x0) ? 2.0 : 0.125; + + for(int n = 0 ; n < data.dust.size(); n++) { + d.dustVc[n](RHO,k,j,i) = (x < x0) ? 1.0 : 16.0; + d.dustVc[n](VX1,k,j,i) = (x < x0) ? 2.0 : 0.125; + } + + + } + } + } + + // Send it all, if needed + d.SyncToDevice(); +} + +// Analyse data to produce an output +void MakeAnalysis(DataBlock & data) { + +} diff --git a/test/Dust/DustyShock/testme.py b/test/Dust/DustyShock/testme.py new file mode 100755 index 00000000..2b03d163 --- /dev/null +++ b/test/Dust/DustyShock/testme.py @@ -0,0 +1,39 @@ +#!/usr/bin/env python3 + +""" + +@author: glesur +""" +import os +import sys +sys.path.append(os.getenv("IDEFIX_DIR")) + +import pytools.idfx_test as tst + +name="dump.0001.dmp" + +def testMe(test): + test.configure() + test.compile() + inifiles=["idefix.ini","idefix-implicit.ini"] + + # loop on all the ini files for this test + for ini in inifiles: + test.run(inputFile=ini) + if test.init: + test.makeReference(filename=name) + test.standardTest() + test.nonRegressionTest(filename=name,tolerance=1e-14) + + +test=tst.idfxTest() + +if not test.all: + if(test.check): + test.checkOnly(filename=name) + else: + testMe(test) +else: + test.noplot = True + test.reconstruction=2 + testMe(test) From 60740e78eecd2e70ea5bd0768ec100f017a8356d Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Tue, 11 Feb 2025 14:20:18 +0100 Subject: [PATCH 06/17] add back dump files for reference --- test/Dust/DustyShock/idefix-implicit.ini | 2 +- test/Dust/DustyShock/idefix.ini | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/test/Dust/DustyShock/idefix-implicit.ini b/test/Dust/DustyShock/idefix-implicit.ini index e83d5bf4..2d92a0a5 100644 --- a/test/Dust/DustyShock/idefix-implicit.ini +++ b/test/Dust/DustyShock/idefix-implicit.ini @@ -31,6 +31,6 @@ X3-beg outflow X3-end outflow [Output] -# dmp 10.0 +dmp 500.0 vtk 500.0 log 1000 diff --git a/test/Dust/DustyShock/idefix.ini b/test/Dust/DustyShock/idefix.ini index 1b9c9ea2..a98ebd09 100644 --- a/test/Dust/DustyShock/idefix.ini +++ b/test/Dust/DustyShock/idefix.ini @@ -30,6 +30,6 @@ X3-beg outflow X3-end outflow [Output] -# dmp 10.0 +dmp 500.0 vtk 500.0 log 1000 From 4b8e277e51e78709642cb8cdd3a2f2d510afd240 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Tue, 11 Feb 2025 14:21:02 +0100 Subject: [PATCH 07/17] update reference files for new dusyShock --- reference | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/reference b/reference index b675bcea..3e155509 160000 --- a/reference +++ b/reference @@ -1 +1 @@ -Subproject commit b675bceaa6aabc01dded346e2d631857f698dc76 +Subproject commit 3e1555090230ceb21ea6428d82c032734fef7488 From e844bc19c14b497b4e590b4a32feeff38d605308 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Tue, 11 Feb 2025 17:03:09 +0100 Subject: [PATCH 08/17] add documentation for implicit dust module --- doc/source/modules/dust.rst | 41 +++++++++++++++++++++++++++++++++---- 1 file changed, 37 insertions(+), 4 deletions(-) diff --git a/doc/source/modules/dust.rst b/doc/source/modules/dust.rst index ca371dad..0cdac49b 100644 --- a/doc/source/modules/dust.rst +++ b/doc/source/modules/dust.rst @@ -39,15 +39,45 @@ which guarantees total momentum conservation. -Drag CFL condition -------------------- -*Idefix* computes the drag terms with a time-explicit scheme. Hence, an addition CFL constraint arrises because of the drag: +Time Integration +---------------- + +*Idefix* can compute the drag with a 2nd order time-explicit (default) or 1st order time-implicit scheme. Each scheme has its pros and cons. The choice +is made by the user in the input file. + +Explicit scheme: drag CFL condition ++++++++++++++++++++++++++++++++++++ + +When using the time explicit scheme, an addition CFL constraint arrises because of the drag: .. math:: dt < \min(\frac{1}{\sum_i\gamma_i(\rho_i+\rho)}) -*Idefix* automatically adjusts the CFL to satisfy this inequality, in addition to the usual CFL condition. +*Idefix* automatically adjusts the CFL to satisfy this inequality, in addition to the usual CFL condition. This can severly limit the time step, +especially when the gas density is high. + +Implicit scheme ++++++++++++++++ + +When the drag is applied with a 1st order implicit scheme, the drag force is applied at the end of each step. In order to avoid +the complete inversion of the system, we follow a simplified inversion procedure, where we first update the gas momentum as: + +.. math:: + + v_g^{(n+1)}=\left(v_g^{(n)}+\sum_i\frac{\rho_i\gamma_i dt}{1+\rho \gamma_i dt}v_i^{(n)}\right)/\left(1+\sum_i\frac{\rho_i\gamma_i dt}{1+\rho\gamma_idt}\right) + +And then update the dust momentum as: + +.. math:: + + v_i^{(n+1)}=\left(v_i^{(n)}+\rho\gamma_i v_g^{(n+1)}dt\right)/(1+\rho\gamma_i dt) + +Note that the latter equation relies on the *updated* gas velocity. + +.. warning:: + While the implicit scheme is more stable than the explicit one, and it does not require any additional CFL condition, it is less accurate and + possibly lead to inacurrate dust velocities when :math:`dt\gg (\gamma_i\rho)^{-1}`. Use it at your own risk. Dust parameters --------------- @@ -66,6 +96,9 @@ The dust module can be enabled adding a block `[Dust]` in your input .ini file. +----------------+-------------------------+---------------------------------------------------------------------------------------------+ | drag_feedback | bool | | (optionnal) whether the gas feedback is enabled (default true). | +----------------+-------------------------+---------------------------------------------------------------------------------------------+ +| drag_implicit | bool | | (optionnal) whether the drag uses a 1st order implicit method. Otherwise use the | +| | | | 2nd order time-explicit scheme (default is false=time explicit) | ++----------------+-------------------------+---------------------------------------------------------------------------------------------+ The drag parameter :math:`\beta_i` above sets the functional form of :math:`\gamma_i(\rho, \rho_i, c_s)` depending on the drag type: From 78fae024a76ee95b94940ffacb8ad36ec0fd2a90 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Thu, 13 Feb 2025 13:57:46 +0100 Subject: [PATCH 09/17] use proper update variables istead of pre-stage ones --- src/fluid/drag.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/fluid/drag.cpp b/src/fluid/drag.cpp index bbbf11b2..b6eb1b3d 100644 --- a/src/fluid/drag.cpp +++ b/src/fluid/drag.cpp @@ -89,7 +89,7 @@ void Drag::AddImplicitBackReaction(const real dt, IdefixArray3D preFactor) real gamma = gammaDrag.GetGamma(k,j,i); // The drag coefficient - const real factor = VcDust(RHO,k,j,i)*gamma*dt/(1+VcGas(RHO,k,j,i)*gamma*dt); + const real factor = UcDust(RHO,k,j,i)*gamma*dt/(1+UcGas(RHO,k,j,i)*gamma*dt); if(isFirst) { preFactor(k,j,i) = factor; } else { @@ -97,8 +97,8 @@ void Drag::AddImplicitBackReaction(const real dt, IdefixArray3D preFactor) } for(int n = MX1 ; n < MX1+COMPONENTS ; n++) { - UcGas(n,k,j,i) += dt * gamma * VcGas(RHO,k,j,i) * UcDust(n,k,j,i) / - (1 + VcGas(RHO,k,j,i)*dt*gamma); + UcGas(n,k,j,i) += dt * gamma * UcGas(RHO,k,j,i) * UcDust(n,k,j,i) / + (1 + UcGas(RHO,k,j,i)*dt*gamma); } }); idfx::popRegion(); @@ -156,8 +156,8 @@ void Drag::AddImplicitFluidMomentum(const real dt) { for(int n = MX1 ; n < MX1+COMPONENTS ; n++) { real oldUc = UcDust(n,k,j,i); - UcDust(n,k,j,i) = (oldUc + dt * gamma * VcDust(RHO,k,j,i) * UcGas(n,k,j,i)) / - (1 + VcGas(RHO,k,j,i)*dt*gamma); + UcDust(n,k,j,i) = (oldUc + dt * gamma * UcDust(RHO,k,j,i) * UcGas(n,k,j,i)) / + (1 + UcGas(RHO,k,j,i)*dt*gamma); #if HAVE_ENERGY == 1 real dp = UcDust(n,k,j,i) - oldUc; @@ -167,7 +167,7 @@ void Drag::AddImplicitFluidMomentum(const real dt) { // TODO(GL): this should be disabled in the case of a true multifluid system where // both fluids have a proper energy equation - UcGas(ENG,k,j,i) -= dp*VcDust(n,k,j,i); + UcGas(ENG,k,j,i) -= dp*UcDust(n,k,j,i)/UcDust(RHO,k,j,i); #endif } }); From 59a91d3cea976f8de290ee461a609cc865ee0fa2 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Thu, 13 Feb 2025 14:07:17 +0100 Subject: [PATCH 10/17] switch off heating terms when feedback is disabled (feedback assumes rho_d is negligible, so this makes sense anyway) --- src/fluid/drag.cpp | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/src/fluid/drag.cpp b/src/fluid/drag.cpp index b6eb1b3d..a1aae441 100644 --- a/src/fluid/drag.cpp +++ b/src/fluid/drag.cpp @@ -36,16 +36,17 @@ void Drag::AddDragForce(const real dt) { for(int n = MX1 ; n < MX1+COMPONENTS ; n++) { real dv = VcDust(n,k,j,i) - VcGas(n,k,j,i); UcDust(n,k,j,i) -= dp*dv; - if(feedback) UcGas(n,k,j,i) += dp*dv; - #if HAVE_ENERGY == 1 - // We add back the energy dissipated for the dust which is not accounted for - // (since there is no energy equation for dust grains) - - // TODO(GL): this should be disabled in the case of a true multifluid system where - // both fluids have a proper energy equation - UcGas(ENG,k,j,i) += dp*dv*VcDust(n,k,j,i); - #endif - } + if(feedback) { + UcGas(n,k,j,i) += dp*dv; + #if HAVE_ENERGY == 1 + // We add back the energy dissipated for the dust which is not accounted for + // (since there is no energy equation for dust grains) + + // TODO(GL): this should be disabled in the case of a true multifluid system where + // both fluids have a proper energy equation + UcGas(ENG,k,j,i) += dp*dv*VcDust(n,k,j,i); + #endif + } // feedback // Cfl constraint real idt = gamma*VcGas(RHO,k,j,i); if(feedback) idt += gamma*VcDust(RHO,k,j,i); @@ -167,7 +168,7 @@ void Drag::AddImplicitFluidMomentum(const real dt) { // TODO(GL): this should be disabled in the case of a true multifluid system where // both fluids have a proper energy equation - UcGas(ENG,k,j,i) -= dp*UcDust(n,k,j,i)/UcDust(RHO,k,j,i); + if(feedback) UcGas(ENG,k,j,i) -= dp*UcDust(n,k,j,i)/UcDust(RHO,k,j,i); #endif } }); From f1141eef04e8955cfa9a0d336b389ab87f9d2790 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Thu, 13 Feb 2025 14:09:21 +0100 Subject: [PATCH 11/17] typo --- src/fluid/drag.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/fluid/drag.cpp b/src/fluid/drag.cpp index a1aae441..8bc76516 100644 --- a/src/fluid/drag.cpp +++ b/src/fluid/drag.cpp @@ -47,6 +47,7 @@ void Drag::AddDragForce(const real dt) { UcGas(ENG,k,j,i) += dp*dv*VcDust(n,k,j,i); #endif } // feedback + } // Cfl constraint real idt = gamma*VcGas(RHO,k,j,i); if(feedback) idt += gamma*VcDust(RHO,k,j,i); From 9b2afc546fae579c34be795b4e26ab4fbb6d2c1d Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Thu, 13 Feb 2025 14:13:00 +0100 Subject: [PATCH 12/17] restore energy conservation in implicit scheme with feedback --- src/fluid/drag.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fluid/drag.cpp b/src/fluid/drag.cpp index 8bc76516..d1b5de95 100644 --- a/src/fluid/drag.cpp +++ b/src/fluid/drag.cpp @@ -169,7 +169,7 @@ void Drag::AddImplicitFluidMomentum(const real dt) { // TODO(GL): this should be disabled in the case of a true multifluid system where // both fluids have a proper energy equation - if(feedback) UcGas(ENG,k,j,i) -= dp*UcDust(n,k,j,i)/UcDust(RHO,k,j,i); + if(feedback) UcGas(ENG,k,j,i) -= dp*VcDust(n,k,j,i); #endif } }); From e8d34bb2b16671598f3fb20914368a3f53a1baea Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Thu, 13 Feb 2025 21:02:41 +0100 Subject: [PATCH 13/17] fix label of dust tracers with non-isothermal equation of state fix #320 --- src/fluid/fluid.hpp | 91 +++++++++++++++++++++++---------------------- 1 file changed, 47 insertions(+), 44 deletions(-) diff --git a/src/fluid/fluid.hpp b/src/fluid/fluid.hpp index ae809b48..44cace6b 100644 --- a/src/fluid/fluid.hpp +++ b/src/fluid/fluid.hpp @@ -585,52 +585,55 @@ Fluid::Fluid(Grid &grid, Input &input, DataBlock *datain, int n) { } for(int i = 0 ; i < Phys::nvar+nTracer ; i++) { - switch(i) { - case RHO: - VcName.push_back("RHO"); - data->dump->RegisterVariable(Vc, outputPrefix+"Vc-RHO", RHO); - break; - case VX1: - VcName.push_back("VX1"); - data->dump->RegisterVariable(Vc, outputPrefix+"Vc-VX1", VX1, IDIR); - break; - case VX2: - VcName.push_back("VX2"); - data->dump->RegisterVariable(Vc, outputPrefix+"Vc-VX2", VX2, JDIR); - break; - case VX3: - VcName.push_back("VX3"); - data->dump->RegisterVariable(Vc, outputPrefix+"Vc-VX3", VX3, KDIR); - break; - case BX1: - VcName.push_back("BX1"); - // never save cell-centered BX1 in dumps - break; - case BX2: - VcName.push_back("BX2"); - #if DIMENSIONS < 2 - data->dump->RegisterVariable(Vc, outputPrefix+"Vc-BX2", BX2, JDIR); - #endif - break; - case BX3: - VcName.push_back("BX3"); - #if DIMENSIONS < 3 - data->dump->RegisterVariable(Vc, outputPrefix+"Vc-BX3", BX3, KDIR); - #endif - break; - case PRS: - VcName.push_back("PRS"); - data->dump->RegisterVariable(Vc, outputPrefix+"Vc-PRS", PRS); - break; - default: - if(i>=Phys::nvar) { - std::string tracerLabel = std::string("TR")+std::to_string(i-Phys::nvar); // ="TRn" - VcName.push_back(tracerLabel); - data->dump->RegisterVariable(Vc, outputPrefix+"Vc-"+tracerLabel, i); - } else { + if(i < Phys::nvar) { + // These are standard variable names + switch(i) { + case RHO: + VcName.push_back("RHO"); + data->dump->RegisterVariable(Vc, outputPrefix+"Vc-RHO", RHO); + break; + case VX1: + VcName.push_back("VX1"); + data->dump->RegisterVariable(Vc, outputPrefix+"Vc-VX1", VX1, IDIR); + break; + case VX2: + VcName.push_back("VX2"); + data->dump->RegisterVariable(Vc, outputPrefix+"Vc-VX2", VX2, JDIR); + break; + case VX3: + VcName.push_back("VX3"); + data->dump->RegisterVariable(Vc, outputPrefix+"Vc-VX3", VX3, KDIR); + break; + case BX1: + VcName.push_back("BX1"); + // never save cell-centered BX1 in dumps + break; + case BX2: + VcName.push_back("BX2"); + #if DIMENSIONS < 2 + data->dump->RegisterVariable(Vc, outputPrefix+"Vc-BX2", BX2, JDIR); + #endif + break; + case BX3: + VcName.push_back("BX3"); + #if DIMENSIONS < 3 + data->dump->RegisterVariable(Vc, outputPrefix+"Vc-BX3", BX3, KDIR); + #endif + break; + case PRS: + VcName.push_back("PRS"); + data->dump->RegisterVariable(Vc, outputPrefix+"Vc-PRS", PRS); + break; + default: + // unknown variable name (this should never happen, but who knows...) VcName.push_back("Vc-"+std::to_string(i)); data->vtk->RegisterVariable(Vc, outputPrefix+"Vc-"+std::to_string(i), i); - } + break; + } + } else { //if(i>=Phys::nvar) + std::string tracerLabel = std::string("TR")+std::to_string(i-Phys::nvar); // ="TRn" + VcName.push_back(tracerLabel); + data->dump->RegisterVariable(Vc, outputPrefix+"Vc-"+tracerLabel, i); } data->vtk->RegisterVariable(Vc, outputPrefix+VcName[i], i); #ifdef WITH_HDF5 From af6b3b7064000de48d69c5034635c520186e4670 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Thu, 20 Feb 2025 21:08:11 +0100 Subject: [PATCH 14/17] update reference test change tolerance for dustywave to cope with the implicit drag --- reference | 2 +- test/Dust/DustyWave/python/testidefix.py | 7 ++----- 2 files changed, 3 insertions(+), 6 deletions(-) diff --git a/reference b/reference index 3e155509..5ecb9588 160000 --- a/reference +++ b/reference @@ -1 +1 @@ -Subproject commit 3e1555090230ceb21ea6428d82c032734fef7488 +Subproject commit 5ecb95884f845a23fa0baf45492770c30e22b4fc diff --git a/test/Dust/DustyWave/python/testidefix.py b/test/Dust/DustyWave/python/testidefix.py index c98bc821..ad52d330 100755 --- a/test/Dust/DustyWave/python/testidefix.py +++ b/test/Dust/DustyWave/python/testidefix.py @@ -35,9 +35,6 @@ # Get the minimal decay rate (this should be the one that pops up) tau=np.amax(np.imag(sol)) - - - # load the dat file produced by the setup raw=np.loadtxt('../timevol.dat',skiprows=1) t=raw[:,0] @@ -58,10 +55,10 @@ # Compute decay rate: tau_measured=t[-1]/(2*np.log(etot[-1]/etot[0])) # error on the decay rate: -error=(tau_measured-tau)/tau +error=np.abs((tau_measured-tau)/tau) print("error=%f"%error) -if(error<0.065): +if(error<0.07): print("Success!") else: print("Failure!") From 9ea85b0edcd1d4d7dd83ce9194dd79ba96f72d45 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Thu, 20 Feb 2025 21:15:35 +0100 Subject: [PATCH 15/17] do not perform implicit capture --- src/fluid/drag.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/fluid/drag.cpp b/src/fluid/drag.cpp index d1b5de95..06a88f02 100644 --- a/src/fluid/drag.cpp +++ b/src/fluid/drag.cpp @@ -71,6 +71,7 @@ void Drag::AddImplicitBackReaction(const real dt, IdefixArray3D preFactor) auto UcGas = this->UcGas; auto VcGas = this->VcGas; auto VcDust = this->VcDust; + auto UcDust = this->UcDust; bool feedback = this->feedback; From ff459c0d8691caf55d222779e5b8face02a5f243 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Thu, 20 Feb 2025 21:21:37 +0100 Subject: [PATCH 16/17] remove reference to useless feedback --- src/fluid/drag.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/fluid/drag.cpp b/src/fluid/drag.cpp index 06a88f02..4cf879cd 100644 --- a/src/fluid/drag.cpp +++ b/src/fluid/drag.cpp @@ -73,8 +73,6 @@ void Drag::AddImplicitBackReaction(const real dt, IdefixArray3D preFactor) auto VcDust = this->VcDust; auto UcDust = this->UcDust; - bool feedback = this->feedback; - bool isFirst = this->instanceNumber == 0; From 3c9f26df2dddbe6d5c7f92b2a388d4b595e3950e Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Fri, 21 Feb 2025 10:42:14 +0100 Subject: [PATCH 17/17] update signature of fargo+planet+dust drag function --- test/Dust/FargoPlanet/setup.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/Dust/FargoPlanet/setup.cpp b/test/Dust/FargoPlanet/setup.cpp index af611911..8447efe0 100644 --- a/test/Dust/FargoPlanet/setup.cpp +++ b/test/Dust/FargoPlanet/setup.cpp @@ -42,7 +42,7 @@ void MySoundSpeed(DataBlock &data, const real t, IdefixArray3D &cs) { // a = 20 cm * beta * (Sigma_phys/100 g.cm^-2) / (rho_s / 2 g.cm^-3) // NB: checked against A. Johansen+ (2007) -void MyDrag(DataBlock *data, real beta, IdefixArray3D &gamma) { +void MyDrag(DataBlock *data, int nSpecie, real beta, IdefixArray3D &gamma) { // Compute the drag coefficient gamma from the input beta auto x1 = data->x[IDIR]; real sigma0 = sigma0Glob;