From 72244f8b857cc19cad88ed03fa83b6d355f615af Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Victor=20R=C3=A9ville?= Date: Tue, 3 Dec 2024 10:57:15 +0100 Subject: [PATCH 01/15] Added new branch with collisionless TC from develop --- src/fluid/braginskii/bragThermalDiffusion.cpp | 2 +- src/fluid/braginskii/bragThermalDiffusion.hpp | 80 ++++++++++++++++--- src/fluid/fluid_defs.hpp | 4 + 3 files changed, 76 insertions(+), 10 deletions(-) diff --git a/src/fluid/braginskii/bragThermalDiffusion.cpp b/src/fluid/braginskii/bragThermalDiffusion.cpp index 82d5dac9..2d340adb 100644 --- a/src/fluid/braginskii/bragThermalDiffusion.cpp +++ b/src/fluid/braginskii/bragThermalDiffusion.cpp @@ -46,7 +46,7 @@ void BragThermalDiffusion::ShowConfig() { } } -void BragThermalDiffusion::EnrollBragThermalDiffusivity(BragDiffusivityFunc myFunc) { +void BragThermalDiffusion::EnrollBragThermalDiffusivity(FourArrayDiffusivityFunc myFunc) { if(this->status.status != UserDefFunction) { IDEFIX_WARNING("Braginskii thermal diffusivity enrollment requires Hydro/BragThermalDiffusion " "to be set to userdef in .ini file"); diff --git a/src/fluid/braginskii/bragThermalDiffusion.hpp b/src/fluid/braginskii/bragThermalDiffusion.hpp index 1fe774a4..21703999 100644 --- a/src/fluid/braginskii/bragThermalDiffusion.hpp +++ b/src/fluid/braginskii/bragThermalDiffusion.hpp @@ -35,11 +35,13 @@ class BragThermalDiffusion { void AddBragDiffusiveFluxLim(int, const real, const IdefixArray4D &); // Enroll user-defined thermal conductivity - void EnrollBragThermalDiffusivity(BragDiffusivityFunc); + void EnrollBragThermalDiffusivity(FourArrayDiffusivityFunc); IdefixArray3D heatSrc; // Source terms of the thermal operator IdefixArray3D knorArr; IdefixArray3D kparArr; + IdefixArray3D alpha; // Transition from Brag to collisionless tc + IdefixArray3D clessQ; // Collisionless tc flux coefficient // pre-computed geometrical factors in non-cartesian geometry IdefixArray1D one_dmu; @@ -50,7 +52,7 @@ class BragThermalDiffusion { // status of the module ParabolicModuleStatus &status; - BragDiffusivityFunc diffusivityFunc; + FourArrayDiffusivityFunc diffusivityFunc; bool haveSlopeLimiter{false}; @@ -110,6 +112,12 @@ BragThermalDiffusion::BragThermalDiffusion(Input &input, Grid &grid, Fluid this->knorArr = IdefixArray3D("BragThermalDiffusionKnorArray",data->np_tot[KDIR], data->np_tot[JDIR], data->np_tot[IDIR]); + this->alpha = IdefixArray3D("ClessThermalDiffusionAlphaArray",data->np_tot[KDIR], + data->np_tot[JDIR], + data->np_tot[IDIR]); + this->clessQ = IdefixArray3D("ClessThermalDiffusionClessQ",data->np_tot[KDIR], + data->np_tot[JDIR], + data->np_tot[IDIR]); } else { IDEFIX_ERROR("Unknown braginskii thermal diffusion definition in idefix.ini. " "Can only be constant or userdef."); @@ -256,7 +264,7 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, if(haveThermalDiffusion == UserDefFunction && dir == IDIR) { if(diffusivityFunc) { idfx::pushRegion("UserDef::BragThermalDiffusivityFunction"); - diffusivityFunc(*this->data, t, kparArr, knorArr); + diffusivityFunc(*this->data, t, kparArr, knorArr, alpha, clessQ); idfx::popRegion(); } else { IDEFIX_ERROR("No user-defined thermal diffusion function has been enrolled"); @@ -274,6 +282,8 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, [[maybe_unused]] real dTi, dTj, dTk, dTn; dTi = dTj = dTk = dTn = ZERO_F; + [[maybe_unused]] real dvp, dvm, dpp, dpm, Pn, Vn; /* For the collisionless / saturated flux */ + real locdmax = 0; /////////////////////////////////////////// // IDIR sweep // @@ -292,9 +302,9 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, kpar = kparConstant; } - EXPAND( Bi = BX_I; , - Bj = BY_I; , - Bk = BZ_I; ) + D_EXPAND( Bi = BX_I; , + Bj = BY_I; , + Bk = BZ_I; ) Bn = BX_I; #if GEOMETRY == CARTESIAN @@ -370,7 +380,9 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, } dTi = D_DX_I_T(Vc)/dx1(i); + #if DIMENSIONS >= 2 + if (haveSlopeLimiter) { dTj = 1./x1(i-1)*SL::PLMLim(SL_DY_T(Vc,k,j,i-1)/dx2(j), SL_DY_T(Vc,k,j+1,i-1)/dx2(j+1)); @@ -399,6 +411,22 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, locdmax = FMAX(kpar,knor)/(0.5*(Vc(RHO,k,j,i) + Vc(RHO,k,j,i - 1)))*gamma_m1; dTn = dTi; + + /* For collisionless / saturated heat flux */ + + dvp = Vc(VX1,k,j,i+1) - Vc(VX1,k,j,i); + dvm = Vc(VX1,k,j,i) - Vc(VX1,k,j,i-1); + + dpp = Vc(PRS,k,j,i+1) - Vc(PRS,k,j,i); + dpm = Vc(PRS,k,j,i) - Vc(PRS,k,j,i-1); + + if(haveSlopeLimiter) { + Vn = Vc(VX1,k,j,i)+0.5*SL::PLMLim(dvp, dvm); + Pn = Vc(PRS,k,j,i)+0.5*SL::PLMLim(dpp, dpm); + } else { + Vn = 0.5*(Vc(VX1,k, j, i) + Vc(VX1,k,j,i+1)); + Pn = 0.5*(Vc(PRS,k, j, i) + Vc(PRS,k,j,i+1)); + } } else if(dir == JDIR) { ////////////// // JDIR @@ -533,6 +561,22 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, locdmax = FMAX(kpar,knor)/(0.5*(Vc(RHO,k,j,i) + Vc(RHO,k,j - 1,i)))*gamma_m1; dTn = dTj; + + /* For the collisionless / saturated flux */ + + dvp = Vc(VX2,k,j+1,i) - Vc(VX2,k,j,i); + dvm = Vc(VX2,k,j,i) - Vc(VX2,k,j-1,i); + + dpp = Vc(PRS,k,j+1,i) - Vc(PRS,k,j,i); + dpm = Vc(PRS,k,j,i) - Vc(PRS,k,j-1,i); + + if(haveSlopeLimiter) { + Vn = Vc(VX2,k,j,i)+0.5*SL::PLMLim(dvp, dvm); + Pn = Vc(PRS,k,j,i)+0.5*SL::PLMLim(dpp, dpm); + } else { + Vn = 0.5*(Vc(VX2,k,j, i) + Vc(VX2,k,j+1,i)); + Pn = 0.5*(Vc(PRS,k,j, i) + Vc(PRS,k,j+1,i)); + } } else if(dir == KDIR) { ////////////// // KDIR @@ -625,11 +669,28 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, locdmax = FMAX(kpar,knor)/(0.5*(Vc(RHO,k,j,i) + Vc(RHO,k-1,j,i)))*gamma_m1; dTn = dTk; + + /* For the collisionless / saturated flux */ + dvp = Vc(VX3,k+1,j,i) - Vc(VX3,k,j,i); + dvm = Vc(VX3,k,j,i) - Vc(VX3,k-1,j,i); + + dpp = Vc(PRS,k+1,j,i) - Vc(PRS,k,j,i); + dpm = Vc(PRS,k,j,i) - Vc(PRS,k-1,j,i); + + if(haveSlopeLimiter) { + Vn = Vc(VX3,k,j,i)+0.5*SL::PLMLim(dvp, dvm); + Pn = Vc(PRS,k,j,i)+0.5*SL::PLMLim(dpp, dpm); + } else { + Vn = 0.5*(Vc(VX3,k,j,i) + Vc(VX3,k+1,j,i)); + Pn = 0.5*(Vc(PRS,k,j,i) + Vc(PRS,k+1,j,i)); + } } // From here, gradients and normal have been computed, so we just need to get the fluxes - bgradT = EXPAND( Bi*dTi , + Bj*dTj, +Bk*dTk); - Bmag = EXPAND( Bi*Bi , + Bj*Bj, + Bk*Bk); + bgradT = D_EXPAND( Bi*dTi , + Bj*dTj, +Bk*dTk); + Bmag = D_EXPAND( Bi*Bi , + Bj*Bj, + Bk*Bk); + //printf("%f , %f\n", Bmag, EXPAND(Bi*Bi, + Bj*Bj, + Bk*Bk)); + //printf("%f , %f \n", 0.5*(Vc(VX1,k,j,i) + Vc(VX1,k,j,i+1)), Vn); Bmag = sqrt(Bmag); Bmag = FMAX(1e-6*SMALL_NUMBER,Bmag); @@ -637,10 +698,11 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, bn = Bn/Bmag; /* -- unit vector component -- */ q = kpar*bgradT*bn + knor*(dTn - bn*bgradT); + q = alpha(k,j,i)*q + (1-alpha(k,j,i))*clessQ(k,j,i)*Pn*Vn; Flux(ENG, k, j, i) -= q; - dMax(k,j,i) = FMAX(dMax(k,j,i),locdmax); + dMax(k,j,i) = FMAX(dMax(k,j,i),locdmax)*alpha(k,j,i); }); idfx::popRegion(); } diff --git a/src/fluid/fluid_defs.hpp b/src/fluid/fluid_defs.hpp index 2269fc80..90e18b87 100644 --- a/src/fluid/fluid_defs.hpp +++ b/src/fluid/fluid_defs.hpp @@ -52,6 +52,10 @@ using DiffusivityFunc = void (*) (DataBlock &, const real t, IdefixArray3D using BragDiffusivityFunc = void (*) (DataBlock &, const real t, IdefixArray3D &, IdefixArray3D &); +using FourArrayDiffusivityFunc = void (*) (DataBlock &, const real t, + IdefixArray3D &, IdefixArray3D &, + IdefixArray3D &, IdefixArray3D &); + // Deprecated signatures using SrcTermFuncOld = void (*) (DataBlock &, const real t, const real dt); From cb3bea853eff493db3d5b03f415c0aa98048c90b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Victor=20R=C3=A9ville?= Date: Thu, 12 Dec 2024 16:29:24 +0100 Subject: [PATCH 02/15] Fixed undeclared CL arrays (GPU compiling) --- reference | 2 +- src/fluid/braginskii/bragThermalDiffusion.hpp | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/reference b/reference index b675bcea..609a6016 160000 --- a/reference +++ b/reference @@ -1 +1 @@ -Subproject commit b675bceaa6aabc01dded346e2d631857f698dc76 +Subproject commit 609a60164461fb2d4fbeffa89ce148acf1191525 diff --git a/src/fluid/braginskii/bragThermalDiffusion.hpp b/src/fluid/braginskii/bragThermalDiffusion.hpp index 21703999..187209ee 100644 --- a/src/fluid/braginskii/bragThermalDiffusion.hpp +++ b/src/fluid/braginskii/bragThermalDiffusion.hpp @@ -260,6 +260,8 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, real kparConstant = this->kpar; IdefixArray3D knorArr = this->knorArr; IdefixArray3D kparArr = this->kparArr; + IdefixArray3D alpha = this->alpha; + IdefixArray3D clessQ = this->clessQ; if(haveThermalDiffusion == UserDefFunction && dir == IDIR) { if(diffusivityFunc) { From 6ec6d0db33c0f1c6ac2637fa400c832ec8a83640 Mon Sep 17 00:00:00 2001 From: Hal Bal Date: Sat, 14 Dec 2024 10:18:43 +0100 Subject: [PATCH 03/15] collisionless tc as an option of bragTDiffusion --- src/fluid/braginskii/bragThermalDiffusion.cpp | 22 ++- src/fluid/braginskii/bragThermalDiffusion.hpp | 155 +++++++++++------- src/fluid/fluid_defs.hpp | 10 +- 3 files changed, 116 insertions(+), 71 deletions(-) diff --git a/src/fluid/braginskii/bragThermalDiffusion.cpp b/src/fluid/braginskii/bragThermalDiffusion.cpp index 2d340adb..43477c9d 100644 --- a/src/fluid/braginskii/bragThermalDiffusion.cpp +++ b/src/fluid/braginskii/bragThermalDiffusion.cpp @@ -27,11 +27,17 @@ void BragThermalDiffusion::ShowConfig() { } else if (status.status==UserDefFunction) { idfx::cout << "Braginskii Thermal Diffusion: ENABLED with user-defined diffusivity function." << std::endl; - if(!diffusivityFunc) { + if(!bragDiffusivityFunc) { IDEFIX_ERROR("No braginskii thermal diffusion function has been enrolled"); } + } else if (status.status==CollisionLess) { + idfx::cout << "CollisionLess / Braginskii Thermal Diffusion: ENABLED with user-defined " + "diffusivity function."<< std::endl; + if(!clessDiffusivityFunc) { + IDEFIX_ERROR("No collisionless / Braginskii thermal diffusion function has been enrolled"); + } } else { - IDEFIX_ERROR("Unknown braginskii thermal diffusion mode"); + IDEFIX_ERROR("Unknown Braginskii thermal diffusion mode"); } if(status.isExplicit) { idfx::cout << "Braginskii Thermal Diffusion: uses an explicit time integration." << std::endl; @@ -46,12 +52,20 @@ void BragThermalDiffusion::ShowConfig() { } } -void BragThermalDiffusion::EnrollBragThermalDiffusivity(FourArrayDiffusivityFunc myFunc) { +void BragThermalDiffusion::EnrollBragThermalDiffusivity(TwoArrayDiffusivityFunc myFunc) { if(this->status.status != UserDefFunction) { IDEFIX_WARNING("Braginskii thermal diffusivity enrollment requires Hydro/BragThermalDiffusion " "to be set to userdef in .ini file"); } - this->diffusivityFunc = myFunc; + this->bragDiffusivityFunc = myFunc; +} + +void BragThermalDiffusion::EnrollClessThermalDiffusivity(FourArrayDiffusivityFunc myFunc) { + if(this->status.status != UserDefFunction) { + IDEFIX_WARNING("Collisionless/Braginskii thermal diffusivity enrollment requires " + "Hydro/BragThermalDiffusion to be set to collisionless in .ini file"); + } + this->clessDiffusivityFunc = myFunc; } void BragThermalDiffusion::AddBragDiffusiveFlux(int dir, const real t, diff --git a/src/fluid/braginskii/bragThermalDiffusion.hpp b/src/fluid/braginskii/bragThermalDiffusion.hpp index 187209ee..8cf892e5 100644 --- a/src/fluid/braginskii/bragThermalDiffusion.hpp +++ b/src/fluid/braginskii/bragThermalDiffusion.hpp @@ -35,13 +35,14 @@ class BragThermalDiffusion { void AddBragDiffusiveFluxLim(int, const real, const IdefixArray4D &); // Enroll user-defined thermal conductivity - void EnrollBragThermalDiffusivity(FourArrayDiffusivityFunc); + void EnrollBragThermalDiffusivity(TwoArrayDiffusivityFunc); + void EnrollClessThermalDiffusivity(FourArrayDiffusivityFunc); IdefixArray3D heatSrc; // Source terms of the thermal operator IdefixArray3D knorArr; IdefixArray3D kparArr; - IdefixArray3D alpha; // Transition from Brag to collisionless tc - IdefixArray3D clessQ; // Collisionless tc flux coefficient + IdefixArray3D clessAlpha; // Transition from Brag to collisionless tc + IdefixArray3D clessBeta; // Collisionless tc flux coefficient // pre-computed geometrical factors in non-cartesian geometry IdefixArray1D one_dmu; @@ -52,7 +53,10 @@ class BragThermalDiffusion { // status of the module ParabolicModuleStatus &status; - FourArrayDiffusivityFunc diffusivityFunc; + TwoArrayDiffusivityFunc bragDiffusivityFunc; + FourArrayDiffusivityFunc clessDiffusivityFunc; + + bool includeCollisionlessTD{false}; bool haveSlopeLimiter{false}; @@ -107,17 +111,26 @@ BragThermalDiffusion::BragThermalDiffusion(Input &input, Grid &grid, Fluid } else if(input.Get("Hydro","bragTDiffusion",2).compare("userdef") == 0) { this->status.status = UserDefFunction; this->kparArr = IdefixArray3D("BragThermalDiffusionKparArray",data->np_tot[KDIR], - data->np_tot[JDIR], - data->np_tot[IDIR]); + data->np_tot[JDIR], + data->np_tot[IDIR]); + this->knorArr = IdefixArray3D("BragThermalDiffusionKnorArray",data->np_tot[KDIR], + data->np_tot[JDIR], + data->np_tot[IDIR]); + } else if(input.Get("Hydro","bragTDiffusion",2).compare("collisionless") == 0) { + this->includeCollisionlessTD = true; + this->status.status = CollisionLess; + this->kparArr = IdefixArray3D("BragThermalDiffusionKparArray",data->np_tot[KDIR], + data->np_tot[JDIR], + data->np_tot[IDIR]); this->knorArr = IdefixArray3D("BragThermalDiffusionKnorArray",data->np_tot[KDIR], - data->np_tot[JDIR], - data->np_tot[IDIR]); - this->alpha = IdefixArray3D("ClessThermalDiffusionAlphaArray",data->np_tot[KDIR], - data->np_tot[JDIR], - data->np_tot[IDIR]); - this->clessQ = IdefixArray3D("ClessThermalDiffusionClessQ",data->np_tot[KDIR], - data->np_tot[JDIR], - data->np_tot[IDIR]); + data->np_tot[JDIR], + data->np_tot[IDIR]); + this->clessAlpha = IdefixArray3D("ClessThermalDiffusionAlphaArray",data->np_tot[KDIR], + data->np_tot[JDIR], + data->np_tot[IDIR]); + this->clessBeta = IdefixArray3D("ClessThermalDiffusionBetaArray",data->np_tot[KDIR], + data->np_tot[JDIR], + data->np_tot[IDIR]); } else { IDEFIX_ERROR("Unknown braginskii thermal diffusion definition in idefix.ini. " "Can only be constant or userdef."); @@ -218,6 +231,7 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, HydroModuleStatus haveThermalDiffusion = this->status.status; bool haveSlopeLimiter = this->haveSlopeLimiter; + bool includeCollisionlessTD = this->includeCollisionlessTD; using SL = SlopeLimiter; int ibeg, iend, jbeg, jend, kbeg, kend; @@ -260,17 +274,26 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, real kparConstant = this->kpar; IdefixArray3D knorArr = this->knorArr; IdefixArray3D kparArr = this->kparArr; - IdefixArray3D alpha = this->alpha; - IdefixArray3D clessQ = this->clessQ; + IdefixArray3D clessAlpha = this->clessAlpha; + IdefixArray3D clessBeta = this->clessBeta; if(haveThermalDiffusion == UserDefFunction && dir == IDIR) { - if(diffusivityFunc) { + if(bragDiffusivityFunc) { idfx::pushRegion("UserDef::BragThermalDiffusivityFunction"); - diffusivityFunc(*this->data, t, kparArr, knorArr, alpha, clessQ); + bragDiffusivityFunc(*this->data, t, kparArr, knorArr); idfx::popRegion(); - } else { + } + else { IDEFIX_ERROR("No user-defined thermal diffusion function has been enrolled"); } + } else if(haveThermalDiffusion == CollisionLess && dir == IDIR) { + if (clessDiffusivityFunc) { + idfx::pushRegion("UserDef::ClessThermalDiffusivityFunction"); + clessDiffusivityFunc(*this->data, t, kparArr, knorArr, clessAlpha, clessBeta); + idfx::popRegion(); + } else { + IDEFIX_ERROR("No user-defined collisionless thermal diffusion function has been enrolled"); + } } idefix_for("BragDiffusiveFlux",kbeg, kend, jbeg, jend, ibeg, iend, @@ -292,7 +315,7 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, /////////////////////////////////////////// if(dir == IDIR) { - if(haveThermalDiffusion == UserDefFunction) { + if(haveThermalDiffusion == UserDefFunction || haveThermalDiffusion == CollisionLess) { knor = HALF_F*(knorArr(k,j,i-1)+knorArr(k,j,i)); if(haveSlopeLimiter) { kpar = 2.*(kparArr(k,j,i-1)*kparArr(k,j,i))/(kparArr(k,j,i-1)+kparArr(k,j,i)); @@ -415,26 +438,27 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, dTn = dTi; /* For collisionless / saturated heat flux */ - - dvp = Vc(VX1,k,j,i+1) - Vc(VX1,k,j,i); - dvm = Vc(VX1,k,j,i) - Vc(VX1,k,j,i-1); - - dpp = Vc(PRS,k,j,i+1) - Vc(PRS,k,j,i); - dpm = Vc(PRS,k,j,i) - Vc(PRS,k,j,i-1); - - if(haveSlopeLimiter) { - Vn = Vc(VX1,k,j,i)+0.5*SL::PLMLim(dvp, dvm); - Pn = Vc(PRS,k,j,i)+0.5*SL::PLMLim(dpp, dpm); - } else { - Vn = 0.5*(Vc(VX1,k, j, i) + Vc(VX1,k,j,i+1)); - Pn = 0.5*(Vc(PRS,k, j, i) + Vc(PRS,k,j,i+1)); + if (includeCollisionlessTD) { + dvp = Vc(VX1,k,j,i+1) - Vc(VX1,k,j,i); + dvm = Vc(VX1,k,j,i) - Vc(VX1,k,j,i-1); + + dpp = Vc(PRS,k,j,i+1) - Vc(PRS,k,j,i); + dpm = Vc(PRS,k,j,i) - Vc(PRS,k,j,i-1); + + if(haveSlopeLimiter) { + Vn = Vc(VX1,k,j,i)+0.5*SL::PLMLim(dvp, dvm); + Pn = Vc(PRS,k,j,i)+0.5*SL::PLMLim(dpp, dpm); + } else { + Vn = 0.5*(Vc(VX1,k, j, i) + Vc(VX1,k,j,i+1)); + Pn = 0.5*(Vc(PRS,k, j, i) + Vc(PRS,k,j,i+1)); + } } } else if(dir == JDIR) { ////////////// // JDIR ///////////// - if(haveThermalDiffusion == UserDefFunction) { + if(haveThermalDiffusion == UserDefFunction || haveThermalDiffusion == CollisionLess) { knor = HALF_F*(knorArr(k,j-1,i)+knorArr(k,j,i)); if(haveSlopeLimiter) { kpar = 2.*(kparArr(k,j-1,i)*kparArr(k,j,i))/(kparArr(k,j-1,i)+kparArr(k,j,i)); @@ -565,25 +589,26 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, dTn = dTj; /* For the collisionless / saturated flux */ - - dvp = Vc(VX2,k,j+1,i) - Vc(VX2,k,j,i); - dvm = Vc(VX2,k,j,i) - Vc(VX2,k,j-1,i); - - dpp = Vc(PRS,k,j+1,i) - Vc(PRS,k,j,i); - dpm = Vc(PRS,k,j,i) - Vc(PRS,k,j-1,i); - - if(haveSlopeLimiter) { - Vn = Vc(VX2,k,j,i)+0.5*SL::PLMLim(dvp, dvm); - Pn = Vc(PRS,k,j,i)+0.5*SL::PLMLim(dpp, dpm); - } else { - Vn = 0.5*(Vc(VX2,k,j, i) + Vc(VX2,k,j+1,i)); - Pn = 0.5*(Vc(PRS,k,j, i) + Vc(PRS,k,j+1,i)); - } + if(includeCollisionlessTD) { + dvp = Vc(VX2,k,j+1,i) - Vc(VX2,k,j,i); + dvm = Vc(VX2,k,j,i) - Vc(VX2,k,j-1,i); + + dpp = Vc(PRS,k,j+1,i) - Vc(PRS,k,j,i); + dpm = Vc(PRS,k,j,i) - Vc(PRS,k,j-1,i); + + if(haveSlopeLimiter) { + Vn = Vc(VX2,k,j,i)+0.5*SL::PLMLim(dvp, dvm); + Pn = Vc(PRS,k,j,i)+0.5*SL::PLMLim(dpp, dpm); + } else { + Vn = 0.5*(Vc(VX2,k,j, i) + Vc(VX2,k,j+1,i)); + Pn = 0.5*(Vc(PRS,k,j, i) + Vc(PRS,k,j+1,i)); + } + } } else if(dir == KDIR) { ////////////// // KDIR ///////////// - if(haveThermalDiffusion == UserDefFunction) { + if(haveThermalDiffusion == UserDefFunction || haveThermalDiffusion == CollisionLess) { knor = HALF_F*(knorArr(k-1,j,i)+knorArr(k,j,i)); if(haveSlopeLimiter) { kpar = 2.*(kparArr(k-1,j,i)*kparArr(k,j,i))/(kparArr(k-1,j,i)+kparArr(k,j,i)); @@ -673,18 +698,20 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, dTn = dTk; /* For the collisionless / saturated flux */ - dvp = Vc(VX3,k+1,j,i) - Vc(VX3,k,j,i); - dvm = Vc(VX3,k,j,i) - Vc(VX3,k-1,j,i); + if(includeCollisionlessTD) { + dvp = Vc(VX3,k+1,j,i) - Vc(VX3,k,j,i); + dvm = Vc(VX3,k,j,i) - Vc(VX3,k-1,j,i); - dpp = Vc(PRS,k+1,j,i) - Vc(PRS,k,j,i); - dpm = Vc(PRS,k,j,i) - Vc(PRS,k-1,j,i); + dpp = Vc(PRS,k+1,j,i) - Vc(PRS,k,j,i); + dpm = Vc(PRS,k,j,i) - Vc(PRS,k-1,j,i); - if(haveSlopeLimiter) { - Vn = Vc(VX3,k,j,i)+0.5*SL::PLMLim(dvp, dvm); - Pn = Vc(PRS,k,j,i)+0.5*SL::PLMLim(dpp, dpm); - } else { - Vn = 0.5*(Vc(VX3,k,j,i) + Vc(VX3,k+1,j,i)); - Pn = 0.5*(Vc(PRS,k,j,i) + Vc(PRS,k+1,j,i)); + if(haveSlopeLimiter) { + Vn = Vc(VX3,k,j,i)+0.5*SL::PLMLim(dvp, dvm); + Pn = Vc(PRS,k,j,i)+0.5*SL::PLMLim(dpp, dpm); + } else { + Vn = 0.5*(Vc(VX3,k,j,i) + Vc(VX3,k+1,j,i)); + Pn = 0.5*(Vc(PRS,k,j,i) + Vc(PRS,k+1,j,i)); + } } } // From here, gradients and normal have been computed, so we just need to get the fluxes @@ -700,11 +727,15 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, bn = Bn/Bmag; /* -- unit vector component -- */ q = kpar*bgradT*bn + knor*(dTn - bn*bgradT); - q = alpha(k,j,i)*q + (1-alpha(k,j,i))*clessQ(k,j,i)*Pn*Vn; - + if(includeCollisionlessTD) { + q = clessAlpha(k,j,i)*q + (1-clessAlpha(k,j,i))*clessBeta(k,j,i)*Pn*Vn; + } Flux(ENG, k, j, i) -= q; + dMax(k,j,i) = FMAX(dMax(k,j,i),locdmax); - dMax(k,j,i) = FMAX(dMax(k,j,i),locdmax)*alpha(k,j,i); + if(includeCollisionlessTD) { + dMax(k,j,i) *= clessAlpha(k,j,i); + } }); idfx::popRegion(); } diff --git a/src/fluid/fluid_defs.hpp b/src/fluid/fluid_defs.hpp index 90e18b87..62a2099c 100644 --- a/src/fluid/fluid_defs.hpp +++ b/src/fluid/fluid_defs.hpp @@ -25,7 +25,7 @@ class Fluid; // Parabolic terms can have different status -enum HydroModuleStatus {Disabled, Constant, UserDefFunction }; +enum HydroModuleStatus {Disabled, Constant, UserDefFunction, CollisionLess}; // Structure to describe the status of parabolic modules struct ParabolicModuleStatus { @@ -49,12 +49,12 @@ using SrcTermFunc = void (*) (Fluid *, const real t, const real dt); using EmfBoundaryFunc = void (*) (DataBlock &, const real t); using DiffusivityFunc = void (*) (DataBlock &, const real t, IdefixArray3D &); -using BragDiffusivityFunc = void (*) (DataBlock &, const real t, - IdefixArray3D &, IdefixArray3D &); +using TwoArrayDiffusivityFunc = void (*) (DataBlock &, const real t, + IdefixArray3D &, IdefixArray3D &); using FourArrayDiffusivityFunc = void (*) (DataBlock &, const real t, - IdefixArray3D &, IdefixArray3D &, - IdefixArray3D &, IdefixArray3D &); + IdefixArray3D &, IdefixArray3D &, + IdefixArray3D &, IdefixArray3D &); // Deprecated signatures using SrcTermFuncOld = void (*) (DataBlock &, const real t, const real dt); From 9ac669865703c68324b6ef35999cc786deaaf8c7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Victor=20R=C3=A9ville?= Date: Mon, 16 Dec 2024 10:32:16 +0100 Subject: [PATCH 04/15] idefix error includes collisionless status --- src/fluid/braginskii/bragThermalDiffusion.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/fluid/braginskii/bragThermalDiffusion.hpp b/src/fluid/braginskii/bragThermalDiffusion.hpp index 8cf892e5..ca1cf115 100644 --- a/src/fluid/braginskii/bragThermalDiffusion.hpp +++ b/src/fluid/braginskii/bragThermalDiffusion.hpp @@ -133,7 +133,7 @@ BragThermalDiffusion::BragThermalDiffusion(Input &input, Grid &grid, Fluid data->np_tot[IDIR]); } else { IDEFIX_ERROR("Unknown braginskii thermal diffusion definition in idefix.ini. " - "Can only be constant or userdef."); + "Can only be constant, userdef or collisionless."); } } else { IDEFIX_ERROR("I cannot create a BragThermalDiffusion object without bragTDiffusion defined" @@ -728,13 +728,13 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, bn = Bn/Bmag; /* -- unit vector component -- */ q = kpar*bgradT*bn + knor*(dTn - bn*bgradT); if(includeCollisionlessTD) { - q = clessAlpha(k,j,i)*q + (1-clessAlpha(k,j,i))*clessBeta(k,j,i)*Pn*Vn; + q = clessAlpha(k,j,i)*q + (1-clessAlpha(k,j,i))*clessBeta(k,j,i)*Pn*Vn; } Flux(ENG, k, j, i) -= q; dMax(k,j,i) = FMAX(dMax(k,j,i),locdmax); if(includeCollisionlessTD) { - dMax(k,j,i) *= clessAlpha(k,j,i); + dMax(k,j,i) *= clessAlpha(k,j,i); } }); idfx::popRegion(); From 94a71d229c565fe21be263c4f06d65b4555477ed Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Victor=20R=C3=A9ville?= Date: Fri, 20 Dec 2024 15:58:55 +0100 Subject: [PATCH 05/15] Added a saturation option for collisionless formulation --- src/fluid/braginskii/bragThermalDiffusion.cpp | 15 +- src/fluid/braginskii/bragThermalDiffusion.hpp | 164 +++++++++++------- src/fluid/fluid_defs.hpp | 2 +- 3 files changed, 113 insertions(+), 68 deletions(-) diff --git a/src/fluid/braginskii/bragThermalDiffusion.cpp b/src/fluid/braginskii/bragThermalDiffusion.cpp index 43477c9d..ee1d5b99 100644 --- a/src/fluid/braginskii/bragThermalDiffusion.cpp +++ b/src/fluid/braginskii/bragThermalDiffusion.cpp @@ -27,15 +27,9 @@ void BragThermalDiffusion::ShowConfig() { } else if (status.status==UserDefFunction) { idfx::cout << "Braginskii Thermal Diffusion: ENABLED with user-defined diffusivity function." << std::endl; - if(!bragDiffusivityFunc) { + if(!bragDiffusivityFunc && !clessDiffusivityFunc) { IDEFIX_ERROR("No braginskii thermal diffusion function has been enrolled"); } - } else if (status.status==CollisionLess) { - idfx::cout << "CollisionLess / Braginskii Thermal Diffusion: ENABLED with user-defined " - "diffusivity function."<< std::endl; - if(!clessDiffusivityFunc) { - IDEFIX_ERROR("No collisionless / Braginskii thermal diffusion function has been enrolled"); - } } else { IDEFIX_ERROR("Unknown Braginskii thermal diffusion mode"); } @@ -50,12 +44,15 @@ void BragThermalDiffusion::ShowConfig() { if(haveSlopeLimiter) { idfx::cout << "Braginskii Thermal Diffusion: uses a slope limiter." << std::endl; } + if(includeCollisionlessTD) { + idfx::cout << "Saturation with collisionless flux is enabled." << std::endl; + } } void BragThermalDiffusion::EnrollBragThermalDiffusivity(TwoArrayDiffusivityFunc myFunc) { if(this->status.status != UserDefFunction) { IDEFIX_WARNING("Braginskii thermal diffusivity enrollment requires Hydro/BragThermalDiffusion " - "to be set to userdef in .ini file"); + "to be set to userdef in .ini file"); } this->bragDiffusivityFunc = myFunc; } @@ -63,7 +60,7 @@ void BragThermalDiffusion::EnrollBragThermalDiffusivity(TwoArrayDiffusivityFunc void BragThermalDiffusion::EnrollClessThermalDiffusivity(FourArrayDiffusivityFunc myFunc) { if(this->status.status != UserDefFunction) { IDEFIX_WARNING("Collisionless/Braginskii thermal diffusivity enrollment requires " - "Hydro/BragThermalDiffusion to be set to collisionless in .ini file"); + "Hydro/BragThermalDiffusion to be set to userdef in .ini file"); } this->clessDiffusivityFunc = myFunc; } diff --git a/src/fluid/braginskii/bragThermalDiffusion.hpp b/src/fluid/braginskii/bragThermalDiffusion.hpp index ca1cf115..db01497b 100644 --- a/src/fluid/braginskii/bragThermalDiffusion.hpp +++ b/src/fluid/braginskii/bragThermalDiffusion.hpp @@ -41,8 +41,8 @@ class BragThermalDiffusion { IdefixArray3D heatSrc; // Source terms of the thermal operator IdefixArray3D knorArr; IdefixArray3D kparArr; - IdefixArray3D clessAlpha; // Transition from Brag to collisionless tc - IdefixArray3D clessBeta; // Collisionless tc flux coefficient + IdefixArray3D clessAlphaArr; // Transition from Brag to collisionless tc + IdefixArray3D clessBetaArr; // Collisionless tc flux coefficient (before pv) // pre-computed geometrical factors in non-cartesian geometry IdefixArray1D one_dmu; @@ -67,6 +67,7 @@ class BragThermalDiffusion { // constant diffusion coefficient (when needed) real knor, kpar; + real clessAlpha, clessBeta; // equation of state (required to get the heat capacity) EquationOfState *eos; @@ -104,36 +105,55 @@ BragThermalDiffusion::BragThermalDiffusion(Input &input, Grid &grid, Fluid IDEFIX_ERROR("Unknown braginskii thermal diffusion limiter in idefix.ini. " "Can only be vanleer, mc or nolimiter."); } - if(input.Get("Hydro","bragTDiffusion",2).compare("constant") == 0) { - this->kpar = input.Get("Hydro","bragTDiffusion",3); - this->knor = input.GetOrSet("Hydro","bragTDiffusion",4,0.); - this->status.status = Constant; - } else if(input.Get("Hydro","bragTDiffusion",2).compare("userdef") == 0) { - this->status.status = UserDefFunction; - this->kparArr = IdefixArray3D("BragThermalDiffusionKparArray",data->np_tot[KDIR], - data->np_tot[JDIR], - data->np_tot[IDIR]); - this->knorArr = IdefixArray3D("BragThermalDiffusionKnorArray",data->np_tot[KDIR], - data->np_tot[JDIR], - data->np_tot[IDIR]); - } else if(input.Get("Hydro","bragTDiffusion",2).compare("collisionless") == 0) { - this->includeCollisionlessTD = true; - this->status.status = CollisionLess; - this->kparArr = IdefixArray3D("BragThermalDiffusionKparArray",data->np_tot[KDIR], - data->np_tot[JDIR], - data->np_tot[IDIR]); - this->knorArr = IdefixArray3D("BragThermalDiffusionKnorArray",data->np_tot[KDIR], - data->np_tot[JDIR], - data->np_tot[IDIR]); - this->clessAlpha = IdefixArray3D("ClessThermalDiffusionAlphaArray",data->np_tot[KDIR], - data->np_tot[JDIR], - data->np_tot[IDIR]); - this->clessBeta = IdefixArray3D("ClessThermalDiffusionBetaArray",data->np_tot[KDIR], - data->np_tot[JDIR], - data->np_tot[IDIR]); + if(input.Get("Hydro","bragTDiffusion",2).compare("nosat") == 0) { + if(input.Get("Hydro","bragTDiffusion",3).compare("constant") == 0) { + this->kpar = input.Get("Hydro","bragTDiffusion",4); + this->knor = input.GetOrSet("Hydro","bragTDiffusion",5,0.); + this->status.status = Constant; + } else if(input.Get("Hydro","bragTDiffusion",3).compare("userdef") == 0) { + this->status.status = UserDefFunction; + this->kparArr = IdefixArray3D("BragThermalDiffusionKparArray",data->np_tot[KDIR], + data->np_tot[JDIR], + data->np_tot[IDIR]); + this->knorArr = IdefixArray3D("BragThermalDiffusionKnorArray",data->np_tot[KDIR], + data->np_tot[JDIR], + data->np_tot[IDIR]); + } else { + IDEFIX_ERROR("Unknown braginskii thermal diffusion definition in idefix.ini. " + "Can only be constant or userdef."); + } + } else if(input.Get("Hydro","bragTDiffusion",2).compare("wcollisionless") == 0.0) { + if(input.Get("Hydro","bragTDiffusion",3).compare("constant") == 0) { + this->includeCollisionlessTD = true; + this->kpar = input.Get("Hydro","bragTDiffusion",4); + this->knor = input.GetOrSet("Hydro","bragTDiffusion",5,0.); + this->clessAlpha = input.GetOrSet("Hydro","bragTDiffusion",6,1.); + this->clessBeta = input.GetOrSet("Hydro","bragTDiffusion",7,0.); + this->status.status = Constant; + } else if(input.Get("Hydro","bragTDiffusion",3).compare("userdef") == 0) { + this->includeCollisionlessTD = true; + this->status.status = UserDefFunction; + this->kparArr = IdefixArray3D("BragThermalDiffusionKparArray",data->np_tot[KDIR], + data->np_tot[JDIR], + data->np_tot[IDIR]); + this->knorArr = IdefixArray3D("BragThermalDiffusionKnorArray",data->np_tot[KDIR], + data->np_tot[JDIR], + data->np_tot[IDIR]); + this->clessAlphaArr = IdefixArray3D("ClessThermalDiffusionAlphaArray", + data->np_tot[KDIR], + data->np_tot[JDIR], + data->np_tot[IDIR]); + this->clessBetaArr = IdefixArray3D("ClessThermalDiffusionBetaArray", + data->np_tot[KDIR], + data->np_tot[JDIR], + data->np_tot[IDIR]); + } else { + IDEFIX_ERROR("Unknown braginskii/collsiionless thermal diffusion definition in idefix.ini. " + "Can only be constant or userdef."); + } } else { - IDEFIX_ERROR("Unknown braginskii thermal diffusion definition in idefix.ini. " - "Can only be constant, userdef or collisionless."); + IDEFIX_ERROR("Unknown braginskii thermal diffusion saturation in idefix.ini. " + "Can only be nosat or wcollisionless."); } } else { IDEFIX_ERROR("I cannot create a BragThermalDiffusion object without bragTDiffusion defined" @@ -156,7 +176,7 @@ BragThermalDiffusion::BragThermalDiffusion(Input &input, Grid &grid, Fluid // and therefore not available as such in the DataBlock. // It is rather defined as PRS/RHO. // Special spatial derivative macros are therefore needed and defined here -// directly at the right cell interface according to the direciton of the flux. +// directly at the right cell interface according to the direction of the flux. #define D_DX_I_T(q) (q(PRS,k,j,i)/q(RHO,k,j,i) - q(PRS,k,j,i - 1)/q(RHO,k,j,i - 1)) #define D_DY_J_T(q) (q(PRS,k,j,i)/q(RHO,k,j,i) - q(PRS,k,j - 1,i)/q(RHO,k,j - 1,i)) #define D_DZ_K_T(q) (q(PRS,k,j,i)/q(RHO,k,j,i) - q(PRS,k - 1,j,i)/q(RHO,k - 1,j,i)) @@ -272,27 +292,31 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, #endif real knorConstant = this->knor; real kparConstant = this->kpar; + real clessAlphaConst = this->clessAlpha; + real clessBetaConst = this->clessBeta; IdefixArray3D knorArr = this->knorArr; IdefixArray3D kparArr = this->kparArr; - IdefixArray3D clessAlpha = this->clessAlpha; - IdefixArray3D clessBeta = this->clessBeta; + IdefixArray3D clessAlphaArr = this->clessAlphaArr; + IdefixArray3D clessBetaArr = this->clessBetaArr; - if(haveThermalDiffusion == UserDefFunction && dir == IDIR) { + if(includeCollisionlessTD == false && haveThermalDiffusion == UserDefFunction && dir == IDIR) { if(bragDiffusivityFunc) { idfx::pushRegion("UserDef::BragThermalDiffusivityFunction"); bragDiffusivityFunc(*this->data, t, kparArr, knorArr); idfx::popRegion(); } else { - IDEFIX_ERROR("No user-defined thermal diffusion function has been enrolled"); + IDEFIX_ERROR("No user-defined Braginskii thermal diffusion function has been enrolled"); } - } else if(haveThermalDiffusion == CollisionLess && dir == IDIR) { + } else if(includeCollisionlessTD == true && haveThermalDiffusion == UserDefFunction + && dir == IDIR) { if (clessDiffusivityFunc) { idfx::pushRegion("UserDef::ClessThermalDiffusivityFunction"); - clessDiffusivityFunc(*this->data, t, kparArr, knorArr, clessAlpha, clessBeta); + clessDiffusivityFunc(*this->data, t, kparArr, knorArr, clessAlphaArr, clessBetaArr); idfx::popRegion(); } else { - IDEFIX_ERROR("No user-defined collisionless thermal diffusion function has been enrolled"); + IDEFIX_ERROR("No user-defined Braginskii/collisionless " + "thermal diffusion function has been enrolled"); } } @@ -306,8 +330,8 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, [[maybe_unused]] real dTi, dTj, dTk, dTn; dTi = dTj = dTk = dTn = ZERO_F; - - [[maybe_unused]] real dvp, dvm, dpp, dpm, Pn, Vn; /* For the collisionless / saturated flux */ + /* For the collisionless / saturated flux */ + [[maybe_unused]] real clessAlpha, clessBeta, dvp, dvm, dpp, dpm, Pn, Vn; real locdmax = 0; /////////////////////////////////////////// @@ -315,16 +339,24 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, /////////////////////////////////////////// if(dir == IDIR) { - if(haveThermalDiffusion == UserDefFunction || haveThermalDiffusion == CollisionLess) { + if(haveThermalDiffusion == UserDefFunction) { knor = HALF_F*(knorArr(k,j,i-1)+knorArr(k,j,i)); if(haveSlopeLimiter) { kpar = 2.*(kparArr(k,j,i-1)*kparArr(k,j,i))/(kparArr(k,j,i-1)+kparArr(k,j,i)); } else { kpar = HALF_F*(kparArr(k,j,i-1)+kparArr(k,j,i)); } + if(includeCollisionlessTD) { + clessAlpha=HALF_F*(clessAlphaArr(k,j,i-1)+clessAlphaArr(k,j,i)); + clessBeta=HALF_F*(clessBetaArr(k,j,i-1)+clessAlphaArr(k,j,i)); + } } else { knor = knorConstant; kpar = kparConstant; + if(includeCollisionlessTD) { + clessAlpha=clessAlphaConst; + clessBeta=clessBetaConst; + } } D_EXPAND( Bi = BX_I; , @@ -446,11 +478,11 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, dpm = Vc(PRS,k,j,i) - Vc(PRS,k,j,i-1); if(haveSlopeLimiter) { - Vn = Vc(VX1,k,j,i)+0.5*SL::PLMLim(dvp, dvm); - Pn = Vc(PRS,k,j,i)+0.5*SL::PLMLim(dpp, dpm); + Vn = Vc(VX1,k,j,i)-HALF_F*SL::PLMLim(dvp, dvm); + Pn = Vc(PRS,k,j,i)-HALF_F*SL::PLMLim(dpp, dpm); } else { - Vn = 0.5*(Vc(VX1,k, j, i) + Vc(VX1,k,j,i+1)); - Pn = 0.5*(Vc(PRS,k, j, i) + Vc(PRS,k,j,i+1)); + Vn = HALF_F*(Vc(VX1,k,j,i) + Vc(VX1,k,j,i-1)); + Pn = HALF_F*(Vc(PRS,k,j,i) + Vc(PRS,k,j,i-1)); } } } else if(dir == JDIR) { @@ -458,16 +490,24 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, // JDIR ///////////// - if(haveThermalDiffusion == UserDefFunction || haveThermalDiffusion == CollisionLess) { + if(haveThermalDiffusion == UserDefFunction) { knor = HALF_F*(knorArr(k,j-1,i)+knorArr(k,j,i)); if(haveSlopeLimiter) { kpar = 2.*(kparArr(k,j-1,i)*kparArr(k,j,i))/(kparArr(k,j-1,i)+kparArr(k,j,i)); } else { kpar = HALF_F*(kparArr(k,j-1,i)+kparArr(k,j,i)); } + if(includeCollisionlessTD) { + clessAlpha=HALF_F*(clessAlphaArr(k,j-1,i)+clessAlphaArr(k,j,i)); + clessBeta=HALF_F*(clessBetaArr(k,j-1,i)+clessAlphaArr(k,j,i)); + } } else { knor = knorConstant; kpar = kparConstant; + if(includeCollisionlessTD) { + clessAlpha=clessAlphaConst; + clessBeta=clessBetaConst; + } } EXPAND( Bi = BX_J; , @@ -597,27 +637,35 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, dpm = Vc(PRS,k,j,i) - Vc(PRS,k,j-1,i); if(haveSlopeLimiter) { - Vn = Vc(VX2,k,j,i)+0.5*SL::PLMLim(dvp, dvm); - Pn = Vc(PRS,k,j,i)+0.5*SL::PLMLim(dpp, dpm); + Vn = Vc(VX2,k,j,i)-HALF_F*SL::PLMLim(dvp, dvm); + Pn = Vc(PRS,k,j,i)-HALF_F*SL::PLMLim(dpp, dpm); } else { - Vn = 0.5*(Vc(VX2,k,j, i) + Vc(VX2,k,j+1,i)); - Pn = 0.5*(Vc(PRS,k,j, i) + Vc(PRS,k,j+1,i)); + Vn = 0.5*(Vc(VX2,k,j-1,i) + Vc(VX2,k,j,i)); + Pn = 0.5*(Vc(PRS,k,j-1,i) + Vc(PRS,k,j,i)); } } } else if(dir == KDIR) { ////////////// // KDIR ///////////// - if(haveThermalDiffusion == UserDefFunction || haveThermalDiffusion == CollisionLess) { + if(haveThermalDiffusion == UserDefFunction) { knor = HALF_F*(knorArr(k-1,j,i)+knorArr(k,j,i)); if(haveSlopeLimiter) { kpar = 2.*(kparArr(k-1,j,i)*kparArr(k,j,i))/(kparArr(k-1,j,i)+kparArr(k,j,i)); } else { kpar = HALF_F*(kparArr(k-1,j,i)+kparArr(k,j,i)); } + if(includeCollisionlessTD) { + clessAlpha=HALF_F*(clessAlphaArr(k-1,j,i)+clessAlphaArr(k,j,i)); + clessBeta=HALF_F*(clessBetaArr(k-1,j,i)+clessBetaArr(k,j,i)); + } } else { knor = knorConstant; kpar = kparConstant; + if(includeCollisionlessTD) { + clessAlpha=clessAlphaConst; + clessBeta=clessBetaConst; + } } Bi = BX_K; @@ -706,11 +754,11 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, dpm = Vc(PRS,k,j,i) - Vc(PRS,k-1,j,i); if(haveSlopeLimiter) { - Vn = Vc(VX3,k,j,i)+0.5*SL::PLMLim(dvp, dvm); - Pn = Vc(PRS,k,j,i)+0.5*SL::PLMLim(dpp, dpm); + Vn = Vc(VX3,k,j,i)-HALF_F*SL::PLMLim(dvp, dvm); + Pn = Vc(PRS,k,j,i)-HALF_F*SL::PLMLim(dpp, dpm); } else { - Vn = 0.5*(Vc(VX3,k,j,i) + Vc(VX3,k+1,j,i)); - Pn = 0.5*(Vc(PRS,k,j,i) + Vc(PRS,k+1,j,i)); + Vn = HALF_F*(Vc(VX3,k-1,j,i) + Vc(VX3,k,j,i)); + Pn = HALF_F*(Vc(PRS,k-1,j,i) + Vc(PRS,k,j,i)); } } } @@ -728,13 +776,13 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, bn = Bn/Bmag; /* -- unit vector component -- */ q = kpar*bgradT*bn + knor*(dTn - bn*bgradT); if(includeCollisionlessTD) { - q = clessAlpha(k,j,i)*q + (1-clessAlpha(k,j,i))*clessBeta(k,j,i)*Pn*Vn; + q = clessAlpha*q + (1-clessAlpha)*clessBeta*Pn*Vn; } Flux(ENG, k, j, i) -= q; dMax(k,j,i) = FMAX(dMax(k,j,i),locdmax); if(includeCollisionlessTD) { - dMax(k,j,i) *= clessAlpha(k,j,i); + dMax(k,j,i) *= clessAlpha; } }); idfx::popRegion(); diff --git a/src/fluid/fluid_defs.hpp b/src/fluid/fluid_defs.hpp index 62a2099c..1d040071 100644 --- a/src/fluid/fluid_defs.hpp +++ b/src/fluid/fluid_defs.hpp @@ -25,7 +25,7 @@ class Fluid; // Parabolic terms can have different status -enum HydroModuleStatus {Disabled, Constant, UserDefFunction, CollisionLess}; +enum HydroModuleStatus {Disabled, Constant, UserDefFunction}; // Structure to describe the status of parabolic modules struct ParabolicModuleStatus { From d6560366a1c3a3d6c8c2dabb6c61426b56aaad71 Mon Sep 17 00:00:00 2001 From: Hal Bal Date: Wed, 15 Jan 2025 10:05:44 +0100 Subject: [PATCH 06/15] Upwind scheme for cless TC with limiters + tests --- reference | 2 +- src/fluid/braginskii/bragThermalDiffusion.hpp | 97 +++++++----- test/MHD/MTI/idefix-rkl.ini | 2 +- test/MHD/MTI/idefix-sl.ini | 2 +- test/MHD/MTI/idefix.ini | 2 +- test/MHD/clessTDiffusion/CMakeLists.txt | 2 + test/MHD/clessTDiffusion/definitions.hpp | 5 + test/MHD/clessTDiffusion/idefix.ini | 41 +++++ test/MHD/clessTDiffusion/python/testidefix.py | 40 +++++ test/MHD/clessTDiffusion/setup.cpp | 149 ++++++++++++++++++ test/MHD/clessTDiffusion/testme.py | 42 +++++ test/MHD/sphBragTDiffusion/idefix.ini | 2 +- 12 files changed, 342 insertions(+), 44 deletions(-) create mode 100644 test/MHD/clessTDiffusion/CMakeLists.txt create mode 100644 test/MHD/clessTDiffusion/definitions.hpp create mode 100644 test/MHD/clessTDiffusion/idefix.ini create mode 100644 test/MHD/clessTDiffusion/python/testidefix.py create mode 100644 test/MHD/clessTDiffusion/setup.cpp create mode 100755 test/MHD/clessTDiffusion/testme.py diff --git a/reference b/reference index 609a6016..99f339c4 160000 --- a/reference +++ b/reference @@ -1 +1 @@ -Subproject commit 609a60164461fb2d4fbeffa89ce148acf1191525 +Subproject commit 99f339c43d98274a925282cc96252c5a0fa9374f diff --git a/src/fluid/braginskii/bragThermalDiffusion.hpp b/src/fluid/braginskii/bragThermalDiffusion.hpp index db01497b..64037420 100644 --- a/src/fluid/braginskii/bragThermalDiffusion.hpp +++ b/src/fluid/braginskii/bragThermalDiffusion.hpp @@ -122,7 +122,7 @@ BragThermalDiffusion::BragThermalDiffusion(Input &input, Grid &grid, Fluid IDEFIX_ERROR("Unknown braginskii thermal diffusion definition in idefix.ini. " "Can only be constant or userdef."); } - } else if(input.Get("Hydro","bragTDiffusion",2).compare("wcollisionless") == 0.0) { + } else if(input.Get("Hydro","bragTDiffusion",2).compare("wcless") == 0.0) { if(input.Get("Hydro","bragTDiffusion",3).compare("constant") == 0) { this->includeCollisionlessTD = true; this->kpar = input.Get("Hydro","bragTDiffusion",4); @@ -148,12 +148,12 @@ BragThermalDiffusion::BragThermalDiffusion(Input &input, Grid &grid, Fluid data->np_tot[JDIR], data->np_tot[IDIR]); } else { - IDEFIX_ERROR("Unknown braginskii/collsiionless thermal diffusion definition in idefix.ini. " + IDEFIX_ERROR("Unknown braginskii/collisionless thermal diffusion definition in idefix.ini. " "Can only be constant or userdef."); } } else { IDEFIX_ERROR("Unknown braginskii thermal diffusion saturation in idefix.ini. " - "Can only be nosat or wcollisionless."); + "Can only be nosat or wcless."); } } else { IDEFIX_ERROR("I cannot create a BragThermalDiffusion object without bragTDiffusion defined" @@ -220,7 +220,7 @@ BragThermalDiffusion::BragThermalDiffusion(Input &input, Grid &grid, Fluid //We now define spatial average macros for the magnetic field. // The magnetic field appears in the expression of the Braginskii heat flux. -// It is therefore needed at the right cell interface according to the direction of the flux. +// It is therefore needed at the left cell interface according to the direction of the flux. #define BX_I Vs(BX1s,k,j,i) #define BY_J Vs(BX2s,k,j,i) #define BZ_K Vs(BX3s,k,j,i) @@ -348,7 +348,7 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, } if(includeCollisionlessTD) { clessAlpha=HALF_F*(clessAlphaArr(k,j,i-1)+clessAlphaArr(k,j,i)); - clessBeta=HALF_F*(clessBetaArr(k,j,i-1)+clessAlphaArr(k,j,i)); + clessBeta=HALF_F*(clessBetaArr(k,j,i-1)+clessBetaArr(k,j,i)); } } else { knor = knorConstant; @@ -470,21 +470,28 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, dTn = dTi; /* For collisionless / saturated heat flux */ - if (includeCollisionlessTD) { - dvp = Vc(VX1,k,j,i+1) - Vc(VX1,k,j,i); - dvm = Vc(VX1,k,j,i) - Vc(VX1,k,j,i-1); + if (includeCollisionlessTD) { + if(haveSlopeLimiter) { + dvp = Vc(VX1,k,j,i+1) - Vc(VX1,k,j,i); + dvm = Vc(VX1,k,j,i) - Vc(VX1,k,j,i-1); - dpp = Vc(PRS,k,j,i+1) - Vc(PRS,k,j,i); - dpm = Vc(PRS,k,j,i) - Vc(PRS,k,j,i-1); + dpp = Vc(PRS,k,j,i+1) - Vc(PRS,k,j,i); + dpm = Vc(PRS,k,j,i) - Vc(PRS,k,j,i-1); - if(haveSlopeLimiter) { - Vn = Vc(VX1,k,j,i)-HALF_F*SL::PLMLim(dvp, dvm); - Pn = Vc(PRS,k,j,i)-HALF_F*SL::PLMLim(dpp, dpm); - } else { - Vn = HALF_F*(Vc(VX1,k,j,i) + Vc(VX1,k,j,i-1)); - Pn = HALF_F*(Vc(PRS,k,j,i) + Vc(PRS,k,j,i-1)); - } - } + /* Upwind scheme */ + if (Vc(VX1,k,j,i) > 0.0) { + Vn = Vc(VX1,k,j,i-1)+HALF_F*SL::PLMLim(dvm, dvp); + Pn = Vc(PRS,k,j,i-1)+HALF_F*SL::PLMLim(dpm, dpp); + } + else { + Vn = Vc(VX1,k,j,i)-HALF_F*SL::PLMLim(dvm, dvp); + Pn = Vc(PRS,k,j,i)-HALF_F*SL::PLMLim(dpm, dpp); + } + } else { + Vn = HALF_F*(Vc(VX1,k,j,i) + Vc(VX1,k,j,i-1)); + Pn = HALF_F*(Vc(PRS,k,j,i) + Vc(PRS,k,j,i-1)); + } + } } else if(dir == JDIR) { ////////////// // JDIR @@ -499,7 +506,7 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, } if(includeCollisionlessTD) { clessAlpha=HALF_F*(clessAlphaArr(k,j-1,i)+clessAlphaArr(k,j,i)); - clessBeta=HALF_F*(clessBetaArr(k,j-1,i)+clessAlphaArr(k,j,i)); + clessBeta=HALF_F*(clessBetaArr(k,j-1,i)+clessBetaArr(k,j,i)); } } else { knor = knorConstant; @@ -629,21 +636,27 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, dTn = dTj; /* For the collisionless / saturated flux */ - if(includeCollisionlessTD) { - dvp = Vc(VX2,k,j+1,i) - Vc(VX2,k,j,i); - dvm = Vc(VX2,k,j,i) - Vc(VX2,k,j-1,i); + if(includeCollisionlessTD) { + if(haveSlopeLimiter) { + dvp = Vc(VX2,k,j+1,i) - Vc(VX2,k,j,i); + dvm = Vc(VX2,k,j,i) - Vc(VX2,k,j-1,i); - dpp = Vc(PRS,k,j+1,i) - Vc(PRS,k,j,i); - dpm = Vc(PRS,k,j,i) - Vc(PRS,k,j-1,i); + dpp = Vc(PRS,k,j+1,i) - Vc(PRS,k,j,i); + dpm = Vc(PRS,k,j,i) - Vc(PRS,k,j-1,i); - if(haveSlopeLimiter) { - Vn = Vc(VX2,k,j,i)-HALF_F*SL::PLMLim(dvp, dvm); - Pn = Vc(PRS,k,j,i)-HALF_F*SL::PLMLim(dpp, dpm); - } else { - Vn = 0.5*(Vc(VX2,k,j-1,i) + Vc(VX2,k,j,i)); - Pn = 0.5*(Vc(PRS,k,j-1,i) + Vc(PRS,k,j,i)); - } - } + /* Upwind scheme */ + if (Vc(VX2,k,j,i) > 0.0) { + Vn = Vc(VX2,k,j-1,i)+HALF_F*SL::PLMLim(dvm, dvp); + Pn = Vc(PRS,k,j-1,i)+HALF_F*SL::PLMLim(dpm, dpp); + } else { + Vn = Vc(VX2,k,j,i)-HALF_F*SL::PLMLim(dvm, dvp); + Pn = Vc(PRS,k,j,i)-HALF_F*SL::PLMLim(dpm, dpp); + } + } else { + Vn = HALF_F*(Vc(VX2,k,j-1,i) + Vc(VX2,k,j,i)); + Pn = HALF_F*(Vc(PRS,k,j-1,i) + Vc(PRS,k,j,i)); + } + } } else if(dir == KDIR) { ////////////// // KDIR @@ -747,15 +760,21 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, /* For the collisionless / saturated flux */ if(includeCollisionlessTD) { - dvp = Vc(VX3,k+1,j,i) - Vc(VX3,k,j,i); - dvm = Vc(VX3,k,j,i) - Vc(VX3,k-1,j,i); + if(haveSlopeLimiter) { + dvp = Vc(VX3,k+1,j,i) - Vc(VX3,k,j,i); + dvm = Vc(VX3,k,j,i) - Vc(VX3,k-1,j,i); - dpp = Vc(PRS,k+1,j,i) - Vc(PRS,k,j,i); - dpm = Vc(PRS,k,j,i) - Vc(PRS,k-1,j,i); + dpp = Vc(PRS,k+1,j,i) - Vc(PRS,k,j,i); + dpm = Vc(PRS,k,j,i) - Vc(PRS,k-1,j,i); - if(haveSlopeLimiter) { - Vn = Vc(VX3,k,j,i)-HALF_F*SL::PLMLim(dvp, dvm); - Pn = Vc(PRS,k,j,i)-HALF_F*SL::PLMLim(dpp, dpm); + /* Upwind scheme */ + if (Vc(VX3,k,j,i) > 0.0) { + Vn = Vc(VX3,k-1,j,i)+HALF_F*SL::PLMLim(dvp, dvm); + Pn = Vc(PRS,k-1,j,i)+HALF_F*SL::PLMLim(dpp, dpm); + } else { + Vn = Vc(VX3,k,j,i)-HALF_F*SL::PLMLim(dvp, dvm); + Pn = Vc(PRS,k,j,i)-HALF_F*SL::PLMLim(dpp, dpm); + } } else { Vn = HALF_F*(Vc(VX3,k-1,j,i) + Vc(VX3,k,j,i)); Pn = HALF_F*(Vc(PRS,k-1,j,i) + Vc(PRS,k,j,i)); diff --git a/test/MHD/MTI/idefix-rkl.ini b/test/MHD/MTI/idefix-rkl.ini index 1b9b11eb..09e9eb72 100644 --- a/test/MHD/MTI/idefix-rkl.ini +++ b/test/MHD/MTI/idefix-rkl.ini @@ -12,7 +12,7 @@ nstages 3 [Hydro] solver hlld gamma 1.6666666666666666 -bragTDiffusion rkl nolimiter userdef +bragTDiffusion rkl nolimiter nosat userdef bragViscosity rkl nolimiter userdef [Gravity] diff --git a/test/MHD/MTI/idefix-sl.ini b/test/MHD/MTI/idefix-sl.ini index f2c9deb2..aa707914 100644 --- a/test/MHD/MTI/idefix-sl.ini +++ b/test/MHD/MTI/idefix-sl.ini @@ -12,7 +12,7 @@ nstages 3 [Hydro] solver hlld gamma 1.6666666666666666 -bragTDiffusion rkl mc userdef +bragTDiffusion rkl mc nosat userdef bragViscosity rkl vanleer userdef [Gravity] diff --git a/test/MHD/MTI/idefix.ini b/test/MHD/MTI/idefix.ini index 7a09c406..a159c4f7 100644 --- a/test/MHD/MTI/idefix.ini +++ b/test/MHD/MTI/idefix.ini @@ -12,7 +12,7 @@ nstages 3 [Hydro] solver hlld gamma 1.6666666666666666 -bragTDiffusion explicit nolimiter userdef +bragTDiffusion explicit nolimiter nosat userdef bragViscosity explicit nolimiter userdef [Gravity] diff --git a/test/MHD/clessTDiffusion/CMakeLists.txt b/test/MHD/clessTDiffusion/CMakeLists.txt new file mode 100644 index 00000000..aa9abdd7 --- /dev/null +++ b/test/MHD/clessTDiffusion/CMakeLists.txt @@ -0,0 +1,2 @@ +enable_idefix_property(Idefix_MHD) +set_idefix_property(Idefix_RECONSTRUCTION Parabolic) diff --git a/test/MHD/clessTDiffusion/definitions.hpp b/test/MHD/clessTDiffusion/definitions.hpp new file mode 100644 index 00000000..21b1f873 --- /dev/null +++ b/test/MHD/clessTDiffusion/definitions.hpp @@ -0,0 +1,5 @@ +#define COMPONENTS 3 +#define DIMENSIONS 1 +#define GEOMETRY SPHERICAL + +#define SMALL_PRESSURE_TEMPERATURE (0.1) diff --git a/test/MHD/clessTDiffusion/idefix.ini b/test/MHD/clessTDiffusion/idefix.ini new file mode 100644 index 00000000..9c444b92 --- /dev/null +++ b/test/MHD/clessTDiffusion/idefix.ini @@ -0,0 +1,41 @@ +[Grid] +X1-grid 1 1.0 256 l 40.0 +X2-grid 1 1.4708 1 u 1.6708 +X3-grid 1 0.9892 1 u 1.0108 + +[TimeIntegrator] +CFL 0.4 +CFL_max_var 1.1 # not used +tstop 300.0 +first_dt 1.e-6 +nstages 2 + +[Hydro] +solver hlld +gamma 1.6666666666666 +bragTDiffusion rkl mc wcless userdef + +[Gravity] +potential central +Mcentral 1.0 + +[Boundary] +X1-beg userdef +X1-end outflow +X2-beg periodic +X2-end periodic +X3-beg periodic +X3-end periodic + +[Setup] +UNIT_DENSITY 1.6726e-16 +UNIT_LENGTH 6.9570e10 +UNIT_VELOCITY 4.3670e7 +cs_vesc 0.26 +va_vesc 0.3 +k0 9e-7 + +[Output] +vtk 100.0 +dmp 300.0 +log 10 diff --git a/test/MHD/clessTDiffusion/python/testidefix.py b/test/MHD/clessTDiffusion/python/testidefix.py new file mode 100644 index 00000000..444816d1 --- /dev/null +++ b/test/MHD/clessTDiffusion/python/testidefix.py @@ -0,0 +1,40 @@ +import os +import sys +sys.path.append(os.getenv("IDEFIX_DIR")) +from pytools.vtk_io import readVTK +import numpy as np +import inifix + +conf = inifix.load("../idefix.ini") +gamma = conf["Hydro"]["gamma"] +kappa = conf["Hydro"]["bragTDiffusion"][-1] + +current_VTK = readVTK('../data.0003.vtk', geometry='spherical') + +rho = current_VTK.data['RHO'].squeeze() +prs = current_VTK.data['PRS'].squeeze() +r = current_VTK.r + +idx1=np.where(r > 20)[0][0] +idx2=np.where(r > 30)[0][0] + +def gamma_prime(gamma, beta): + return (gamma+beta*(gamma-1))/(1+beta*(gamma-1)) + +success = True +eps = 1e-2 + +gamma_eff = np.gradient(prs, r)/np.gradient(rho, r)*rho/prs + +Error=np.abs(gamma_eff[idx1:idx2]-gamma_prime(gamma, 1.5)).mean() +if Error > eps: + success=False + +if success: + print("SUCCESS") + print("Error: {0}".format(Error)) + sys.exit(0) +else: + print("Failed") + print("Error: {0}".format(Error)) + sys.exit(1) diff --git a/test/MHD/clessTDiffusion/setup.cpp b/test/MHD/clessTDiffusion/setup.cpp new file mode 100644 index 00000000..624e1d72 --- /dev/null +++ b/test/MHD/clessTDiffusion/setup.cpp @@ -0,0 +1,149 @@ +#include "idefix.hpp" +#include "setup.hpp" + +real udenGlob; +real ulenGlob; +real uvelGlob; +real gammaGlob; +real cs_vescGlob; +real va_vescGlob; +real k0_Glob; +real ParkerWind(real); + +void MyClessThermalConductivity(DataBlock &data, const real t, IdefixArray3D &kparArr, IdefixArray3D &knorArr, IdefixArray3D &alpha, IdefixArray3D &clessQ) { + IdefixArray4D Vc = data.hydro->Vc; + IdefixArray1D x1 = data.x[IDIR]; + IdefixArray1D x2 = data.x[JDIR]; + real norm = 1.6726e-24*0.5/(udenGlob*uvelGlob*ulenGlob*1.3807e-16); + real uTemp=0.5*uvelGlob*uvelGlob*1.6726e-24/1.3807e-16; + real k0 = k0_Glob*norm; + idefix_for("MyThConductivity",0,data.np_tot[KDIR],0,data.np_tot[JDIR],0,data.np_tot[IDIR], + KOKKOS_LAMBDA (int k, int j, int i) { + kparArr(k,j,i) = k0*pow(Vc(PRS,k,j,i)/Vc(RHO,k,j,i)*uTemp,2.5); + knorArr(k,j,i) = 0.; + alpha(k,j,i) = (1.0-tanh(x1(i)-10))/2; + clessQ(k,j,i) = -1.5; + }); +} + +// User-defined boundaries +void UserDefBoundary(Hydro *hydro, int dir, BoundarySide side, real t) { + + DataBlock *data = hydro->data; + + if( (dir==IDIR) && (side == left)) { + IdefixArray4D Vc = hydro->Vc; + IdefixArray4D Vs = hydro->Vs; + IdefixArray1D x1 = data->x[IDIR]; + + real rc,vc,vwind0; + real cs=cs_vescGlob*sqrt(2.); + real va_vesc = va_vescGlob; + real PonRho; + + PonRho = cs*cs; + rc = 0.25 / (cs_vescGlob*cs_vescGlob); + vc = cs; + vwind0 = ParkerWind(1./rc) * cs; + + hydro->boundary->BoundaryFor("UserDefBoundary", dir, side, + KOKKOS_LAMBDA (int k, int j, int i) { + real r = x1(i); + real mu = va_vesc * sqrt(2.); + Vc(RHO,k,j,i) = vwind0/(vwind0 * r * r); + Vc(PRS,k,j,i) = PonRho * Vc(RHO, k, j, i); + Vc(VX1,k,j,i) = vwind0; + Vc(VX2,k,j,i) = 0.0; + Vc(VX3,k,j,i) = 0.0; + Vc(BX1,k,j,i) = mu / (r*r); + Vc(BX2,k,j,i) = 0.0; + Vc(BX3,k,j,i) = 0.0; + + }); + } +} + +// 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) { + // Set the function for userdefboundary + data.hydro->EnrollUserDefBoundary(&UserDefBoundary); + data.hydro->bragThermalDiffusion->EnrollClessThermalDiffusivity(&MyClessThermalConductivity); + gammaGlob=input.Get("Hydro", "gamma", 0); + udenGlob=input.Get("Setup", "UNIT_DENSITY",0); + ulenGlob=input.Get("Setup", "UNIT_LENGTH",0); + uvelGlob=input.Get("Setup", "UNIT_VELOCITY",0); + cs_vescGlob=input.Get("Setup", "cs_vesc", 0); + va_vescGlob=input.Get("Setup", "va_vesc", 0); + k0_Glob = input.Get("Setup", "k0", 0); +} + +// 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); + + real r,th,rl; + real PonRho, vwind0, rc, vc; + real cs=cs_vescGlob*sqrt(2.); + + + rc = 0.25 / (cs_vescGlob*cs_vescGlob); + vwind0 = ParkerWind(1./rc) * cs; + PonRho = cs*cs; + real mu = va_vescGlob * sqrt(2.); + + 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++) { + r=d.x[IDIR](i); + th=d.x[JDIR](j); + real vwind; + + vwind = ParkerWind(r/rc) * cs; + + d.Vc(RHO,k,j,i) = 1.0*vwind0/(vwind * r * r); + d.Vc(PRS,k,j,i) = PonRho * d.Vc(RHO, k, j, i); + d.Vc(VX1,k,j,i) = vwind; + d.Vc(VX2,k,j,i) = 0.0; + d.Vc(VX3,k,j,i) = 0.0; + + rl=d.xl[IDIR](i); // Radius on the left side of the cell + d.Vs(BX1s, k, j, i) = mu / (rl*rl); + + } + } + } + // Send it all, if needed + d.SyncToDevice(); +} + +// Analyse data to produce an output +void MakeAnalysis(DataBlock & data) { +} + +/**************************************************/ +real ParkerWind(real x) +/* Parker wind velocity in unit of iso sound speed + x = radius / critical radius. +**************************************************/ +{ + real v, f; + real vref; + + v = 1e-7; + f = v*v-2*log(v)-4/x-4*log(x)+3; + if (x>1) {v=10;} + while (fabs(f) > 1e-10){ + vref = v; + v = v - 0.5*f/(v-1/v); + while (v < 0){ + v = (v + vref)/2; + } + f = v*v-2*log(v)-4/x-4*log(x)+3; + } + return v; +} diff --git a/test/MHD/clessTDiffusion/testme.py b/test/MHD/clessTDiffusion/testme.py new file mode 100755 index 00000000..6ed5b164 --- /dev/null +++ b/test/MHD/clessTDiffusion/testme.py @@ -0,0 +1,42 @@ +#!/usr/bin/env python3 + +""" + +@author: glesur +""" +import os +import sys +sys.path.append(os.getenv("IDEFIX_DIR")) + +import pytools.idfx_test as tst + + +def testMe(test): + test.configure() + test.compile() + inifiles=["idefix.ini"] + + # loop on all the ini files for this test + name = "dump.0001.dmp" + for ini in inifiles: + test.run(inputFile=ini) + test.standardTest() + if test.init: + test.makeReference(filename=name) + test.nonRegressionTest(filename=name, tolerance=2e-15) + + +test=tst.idfxTest() + +if not test.all: + if(test.check): + test.checkOnly(filename="dump.0001.dmp") + else: + testMe(test) +else: + test.noplot = True + test.vectPot=False + test.single=False + test.reconstruction=2 + test.mpi=False + testMe(test) diff --git a/test/MHD/sphBragTDiffusion/idefix.ini b/test/MHD/sphBragTDiffusion/idefix.ini index 35ac54e8..795ae7bb 100644 --- a/test/MHD/sphBragTDiffusion/idefix.ini +++ b/test/MHD/sphBragTDiffusion/idefix.ini @@ -13,7 +13,7 @@ nstages 3 [Hydro] solver hlld gamma 1.4 -bragTDiffusion rkl nolimiter constant 50.0 +bragTDiffusion rkl nolimiter nosat constant 50.0 [Setup] amplitude 1e-4 From 79e70cf233f0ef0dec55e916b8c5ad723ba123ac Mon Sep 17 00:00:00 2001 From: Jean Kempf Date: Thu, 16 Jan 2025 14:12:10 +0100 Subject: [PATCH 07/15] updated ShowConfig functions from braginskii tc and viscosity to display which slope limiter is used, fixed typo in the braginskii doc --- doc/source/modules/braginskii.rst | 2 +- src/fluid/braginskii/bragThermalDiffusion.cpp | 12 ++++++++++-- src/fluid/braginskii/bragThermalDiffusion.hpp | 5 +++++ src/fluid/braginskii/bragViscosity.cpp | 10 +++++++++- src/fluid/braginskii/bragViscosity.hpp | 4 ++++ 5 files changed, 29 insertions(+), 4 deletions(-) diff --git a/doc/source/modules/braginskii.rst b/doc/source/modules/braginskii.rst index e657cedd..ccc841f4 100644 --- a/doc/source/modules/braginskii.rst +++ b/doc/source/modules/braginskii.rst @@ -37,7 +37,7 @@ though adapted to vector quantities. cell interface by a simple arithmetic average (Eq. (6)-(7) from Sharma & Hammett 2007). However in the same paper, the authors showed that this implementation can lead to - unphysical heat flux from high to low temperature regions. + unphysical heat flux from low to high temperature regions. So we implemented slope limiters for the computation of these transverse heat fluxes, as described in Eq. (17) from Sharma & Hammett (2007). Only the van Leer and the Monotonized Central (MC) limiters are available diff --git a/src/fluid/braginskii/bragThermalDiffusion.cpp b/src/fluid/braginskii/bragThermalDiffusion.cpp index ee1d5b99..a898ab8b 100644 --- a/src/fluid/braginskii/bragThermalDiffusion.cpp +++ b/src/fluid/braginskii/bragThermalDiffusion.cpp @@ -42,10 +42,18 @@ void BragThermalDiffusion::ShowConfig() { IDEFIX_ERROR("Unknown time integrator for braginskii thermal diffusion."); } if(haveSlopeLimiter) { - idfx::cout << "Braginskii Thermal Diffusion: uses a slope limiter." << std::endl; + if(haveMonotizedCentral) { + idfx::cout << "Braginskii Thermal Diffusion: uses the monotonized central slope limiter." << std::endl; + } + else if(haveVanLeer) { + idfx::cout << "Braginskii Thermal Diffusion: uses the van Leer slope limiter." << std::endl; + } + else { + IDEFIX_ERROR("Unknown slope limiter for braginskii thermal diffusion."); + } } if(includeCollisionlessTD) { - idfx::cout << "Saturation with collisionless flux is enabled." << std::endl; + idfx::cout << "Braginskii Thermal Diffusion: saturation with collisionless flux is enabled." << std::endl; } } diff --git a/src/fluid/braginskii/bragThermalDiffusion.hpp b/src/fluid/braginskii/bragThermalDiffusion.hpp index 64037420..2a103a49 100644 --- a/src/fluid/braginskii/bragThermalDiffusion.hpp +++ b/src/fluid/braginskii/bragThermalDiffusion.hpp @@ -59,6 +59,9 @@ class BragThermalDiffusion { bool includeCollisionlessTD{false}; bool haveSlopeLimiter{false}; + bool haveMonotizedCentral{false}; + bool haveVanLeer{false}; + bool haveMinmod{false}; // helper array IdefixArray4D &Vc; @@ -91,9 +94,11 @@ BragThermalDiffusion::BragThermalDiffusion(Input &input, Grid &grid, Fluid if(input.CheckEntry("Hydro","bragTDiffusion")>=0) { if(input.Get("Hydro","bragTDiffusion",1).compare("mc") == 0) { this->haveSlopeLimiter = true; + this->haveMonotizedCentral = true; limiter = PLMLimiter::McLim; } else if(input.Get("Hydro","bragTDiffusion",1).compare("vanleer") == 0) { this->haveSlopeLimiter = true; + this->haveVanLeer = true; limiter = PLMLimiter::VanLeer; } else if(input.Get("Hydro","bragTDiffusion",1).compare("minmod") == 0) { IDEFIX_ERROR("The minmod slope limiter is not available because it has been " diff --git a/src/fluid/braginskii/bragViscosity.cpp b/src/fluid/braginskii/bragViscosity.cpp index 6d5007b5..e05019a4 100644 --- a/src/fluid/braginskii/bragViscosity.cpp +++ b/src/fluid/braginskii/bragViscosity.cpp @@ -59,7 +59,15 @@ void BragViscosity::ShowConfig() { IDEFIX_ERROR("Unknown time integrator for braginskii viscosity."); } if(haveSlopeLimiter) { - idfx::cout << "Braginskii Viscosity: uses a slope limiter." << std::endl; + if(haveMonotizedCentral) { + idfx::cout << "Braginskii Viscosity: uses the monotonized central slope limiter." << std::endl; + } + else if(haveVanLeer) { + idfx::cout << "Braginskii Viscosity: uses the van Leer slope limiter." << std::endl; + } + else { + IDEFIX_ERROR("Unknown slope limiter for braginskii viscosity."); + } } #if GEOMETRY == CYLINDRICAL || GEOMETRY == POLAR diff --git a/src/fluid/braginskii/bragViscosity.hpp b/src/fluid/braginskii/bragViscosity.hpp index 7a8d2b75..5340110a 100644 --- a/src/fluid/braginskii/bragViscosity.hpp +++ b/src/fluid/braginskii/bragViscosity.hpp @@ -52,6 +52,8 @@ class BragViscosity { DiffusivityFunc bragViscousDiffusivityFunc; bool haveSlopeLimiter{false}; + bool haveMonotizedCentral{false}; + bool haveVanLeer{false}; IdefixArray4D &Vc; IdefixArray4D &Vs; @@ -78,9 +80,11 @@ BragViscosity::BragViscosity(Input &input, Grid &grid, Fluid *hydroin): if(input.CheckEntry("Hydro","bragViscosity")>=0) { if(input.Get("Hydro","bragViscosity",1).compare("vanleer") == 0) { this->haveSlopeLimiter = true; + this->haveVanLeer = true; limiter = PLMLimiter::VanLeer; } else if(input.Get("Hydro","bragViscosity",1).compare("mc") == 0) { this->haveSlopeLimiter = true; + this->haveMonotizedCentral = true; limiter = PLMLimiter::McLim; } else if (input.Get("Hydro","bragViscosity",1).compare("nolimiter") == 0) { this->haveSlopeLimiter = false; From a61eb61af47c10e7caf0810a0560457af1ea809d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Victor=20R=C3=A9ville?= Date: Mon, 20 Jan 2025 12:33:43 +0100 Subject: [PATCH 08/15] Updated doc for Collisionless TC --- doc/source/modules/braginskii.rst | 63 ++++++++++++++++++++++++++++--- reference | 2 +- 2 files changed, 59 insertions(+), 6 deletions(-) diff --git a/doc/source/modules/braginskii.rst b/doc/source/modules/braginskii.rst index ccc841f4..4c1c4dd6 100644 --- a/doc/source/modules/braginskii.rst +++ b/doc/source/modules/braginskii.rst @@ -4,7 +4,7 @@ Braginskii module =================== Equations solved and methods ---------------------------- +---------------------------- The ``Braginskii`` module implements the anisotropic heat and momentum fluxes specific to weakly collisional, magnetised plasma like the intracluster medium @@ -72,18 +72,69 @@ of the Braginskii heat flux and viscosity. .. _braginskiiParameterSection: + +Saturation with collisionless heat flux +--------------------------------------- + +The ``Braginskii`` module can include a collisionless saturation of the Braginskii heat flux. +The heat flux is then computed as follows: + +:math:`q = \alpha (q_B + q_\perp) + (1-\alpha)\beta*p*v`, + +where :math:`\alpha \in [0,1]` controls the transition between the Braginskii heat flux and the collisionless heat flux +and :math:`\beta` controls the amplitude of the collisionless heat flux (typically :math:`\beta \in [1,4]`, see Hollweg 1976). + +.. note:: + As a result, even with :math:`\kappa_\perp = 0`, the heat flux is no longer necessarilly strictly aligned with the magnetic field. +.. note:: + The collisionless heat flux is a hyperbolic term and the diffusion coefficient is set proportional to :math:`\alpha`. +.. note:: + If selected, slope limiters are also used in the collisionless flux, where an upwind scheme has been implemented for stability. +.. note:: + This saturation has been thought to be used mostly using the uerdef function that takes four arrays as input and is enrolled through + ``data.hydro->bragThermalDiffusion->EnrollClessThermalDiffusivity()`` + Main parameters of the module ----------------------------- The ``Braginskii`` module can be enabled adding one or two lines in the ``[Hydro]`` section -starting with the keyword -`bragTDiffusion` or/and *bragViscosity*. The following table summarises the different options +starting with the keyword `bragTDiffusion` or/and *bragViscosity*. The following tables summarise the different options associated to the activation of the Braginskii module: +--------+-----------------------+-------------------------+---------------------------------------------------------------------------------------+ | Column | Entry name | Parameter type | Comment | +========+=======================+=========================+=======================================================================================+ -| 0 | bragModule | string | | Activates Braginskii diffusion. Can be ``bragTDiffusion`` or ``bragViscosity``. | +| 0 | bragTDiffusion | string | | Activates Braginskii thermal diffusion. | ++--------+-----------------------+-------------------------+---------------------------------------------------------------------------------------+ +| 1 | integration | string | | Specifies the type of scheme to be used to integrate the parabolic term. | +| | | | | Can be ``rkl`` or ``explicit``. | ++--------+-----------------------+-------------------------+---------------------------------------------------------------------------------------+ +| 2 | slope limiter | string | | Choose the type of limiter to be used to compute anisotropic transverse flux terms. | +| | | | | Can be ``mc``, ``vanleer`` or ``nolimiter``. | ++--------+-----------------------+-------------------------+---------------------------------------------------------------------------------------+ +| 3 | saturation mode | string | | Include or not collisionless saturation. Can be ``nosat`` or ``wcless``. | ++--------+-----------------------+-------------------------+---------------------------------------------------------------------------------------+ +| 4 | diffusivity type | string | | Specifies the type of diffusivity wanted. Can be ``constant`` or ``userdef``. | ++--------+-----------------------+-------------------------+---------------------------------------------------------------------------------------+ +| 5 | parallel diffusivity | real | | Mandatory if the diffusivity type is ``constant``. Not needed otherwise. | +| | | | | Value of the parallel diffusivity. Should be a real number. | ++--------+-----------------------+-------------------------+---------------------------------------------------------------------------------------+ +| 6 | normal diffusivity | real | | When bragModule ``bragTDiffusion`` and diffusivity type ``constant``, | +| | | | | value of the normal diffusivity. Should be a real number. | ++--------+-----------------------+-------------------------+---------------------------------------------------------------------------------------+ +| 7 | alpha collisionless | real | | If the diffusivity type is ``constant`` and saturation is ``wcless``. | +| | | | | Set to 1 if not provided. | ++--------+-----------------------+-------------------------+---------------------------------------------------------------------------------------+ +| 8 | beta collisionless | real | | If the diffusivity type is ``constant`` and saturation is ``wcless``. | +| | | | | Set to 0 if not provided. | ++--------+-----------------------+-------------------------+---------------------------------------------------------------------------------------+ + +for the *bragViscosity*: + ++--------+-----------------------+-------------------------+---------------------------------------------------------------------------------------+ +| Column | Entry name | Parameter type | Comment | ++========+=======================+=========================+=======================================================================================+ +| 0 | bragViscosity | string | | Activates Braginskii viscosity. | +--------+-----------------------+-------------------------+---------------------------------------------------------------------------------------+ | 1 | integration | string | | Specifies the type of scheme to be used to integrate the parabolic term. | | | | | | Can be ``rkl`` or ``explicit``. | @@ -101,7 +152,7 @@ associated to the activation of the Braginskii module: +--------+-----------------------+-------------------------+---------------------------------------------------------------------------------------+ Numerical checks ---------------- +---------------- In Cartesian geometry, the ``Braginskii`` module has been tested with many setups and in all configurations of magnetic polarisation: @@ -119,3 +170,5 @@ The same goes for the anisotropic heat flux in Cylindrical/Polar geometry while the anisotropic viscosity has *never* been tested in this geometry. In spherical geometry, both ``Braginskii`` operators have been partially validated (diffusion along the polar axis has not been directly tested). + +The collisionless saturation has been tested in 1D and 2D spherical geometry. diff --git a/reference b/reference index 99f339c4..609a6016 160000 --- a/reference +++ b/reference @@ -1 +1 @@ -Subproject commit 99f339c43d98274a925282cc96252c5a0fa9374f +Subproject commit 609a60164461fb2d4fbeffa89ce148acf1191525 From 7bef7a974c409fea952dc75ea5f81e51f657e783 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Victor=20R=C3=A9ville?= Date: Mon, 20 Jan 2025 16:46:27 +0100 Subject: [PATCH 09/15] Corrected doc typo --- doc/source/modules/braginskii.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/source/modules/braginskii.rst b/doc/source/modules/braginskii.rst index 4c1c4dd6..d9cb4f2d 100644 --- a/doc/source/modules/braginskii.rst +++ b/doc/source/modules/braginskii.rst @@ -91,7 +91,7 @@ and :math:`\beta` controls the amplitude of the collisionless heat flux (typical .. note:: If selected, slope limiters are also used in the collisionless flux, where an upwind scheme has been implemented for stability. .. note:: - This saturation has been thought to be used mostly using the uerdef function that takes four arrays as input and is enrolled through + This saturation has been thought to be used mostly using the userdef function that takes four arrays as input and is enrolled through ``data.hydro->bragThermalDiffusion->EnrollClessThermalDiffusivity()`` Main parameters of the module From 8f48bc2b34f5223be8de1faa61d2c5bd05cf95b8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Victor=20R=C3=A9ville?= Date: Tue, 21 Jan 2025 13:22:36 +0100 Subject: [PATCH 10/15] Corrected typo --- src/fluid/braginskii/bragThermalDiffusion.cpp | 14 +++++++------- src/fluid/braginskii/bragThermalDiffusion.hpp | 4 ++-- src/fluid/braginskii/bragViscosity.cpp | 9 ++++----- 3 files changed, 13 insertions(+), 14 deletions(-) diff --git a/src/fluid/braginskii/bragThermalDiffusion.cpp b/src/fluid/braginskii/bragThermalDiffusion.cpp index a898ab8b..0dd0b033 100644 --- a/src/fluid/braginskii/bragThermalDiffusion.cpp +++ b/src/fluid/braginskii/bragThermalDiffusion.cpp @@ -42,18 +42,18 @@ void BragThermalDiffusion::ShowConfig() { IDEFIX_ERROR("Unknown time integrator for braginskii thermal diffusion."); } if(haveSlopeLimiter) { - if(haveMonotizedCentral) { - idfx::cout << "Braginskii Thermal Diffusion: uses the monotonized central slope limiter." << std::endl; - } - else if(haveVanLeer) { + if(haveMonotonizedCentral) { + idfx::cout << "Braginskii Thermal Diffusion: " + "uses the monotonized central slope limiter." << std::endl; + } else if(haveVanLeer) { idfx::cout << "Braginskii Thermal Diffusion: uses the van Leer slope limiter." << std::endl; - } - else { + } else { IDEFIX_ERROR("Unknown slope limiter for braginskii thermal diffusion."); } } if(includeCollisionlessTD) { - idfx::cout << "Braginskii Thermal Diffusion: saturation with collisionless flux is enabled." << std::endl; + idfx::cout << "Braginskii Thermal Diffusion: saturation" + " with collisionless flux is enabled." << std::endl; } } diff --git a/src/fluid/braginskii/bragThermalDiffusion.hpp b/src/fluid/braginskii/bragThermalDiffusion.hpp index 2a103a49..d05c20b4 100644 --- a/src/fluid/braginskii/bragThermalDiffusion.hpp +++ b/src/fluid/braginskii/bragThermalDiffusion.hpp @@ -59,7 +59,7 @@ class BragThermalDiffusion { bool includeCollisionlessTD{false}; bool haveSlopeLimiter{false}; - bool haveMonotizedCentral{false}; + bool haveMonotonizedCentral{false}; bool haveVanLeer{false}; bool haveMinmod{false}; @@ -94,7 +94,7 @@ BragThermalDiffusion::BragThermalDiffusion(Input &input, Grid &grid, Fluid if(input.CheckEntry("Hydro","bragTDiffusion")>=0) { if(input.Get("Hydro","bragTDiffusion",1).compare("mc") == 0) { this->haveSlopeLimiter = true; - this->haveMonotizedCentral = true; + this->haveMonotonizedCentral = true; limiter = PLMLimiter::McLim; } else if(input.Get("Hydro","bragTDiffusion",1).compare("vanleer") == 0) { this->haveSlopeLimiter = true; diff --git a/src/fluid/braginskii/bragViscosity.cpp b/src/fluid/braginskii/bragViscosity.cpp index e05019a4..e53ab2c6 100644 --- a/src/fluid/braginskii/bragViscosity.cpp +++ b/src/fluid/braginskii/bragViscosity.cpp @@ -60,12 +60,11 @@ void BragViscosity::ShowConfig() { } if(haveSlopeLimiter) { if(haveMonotizedCentral) { - idfx::cout << "Braginskii Viscosity: uses the monotonized central slope limiter." << std::endl; - } - else if(haveVanLeer) { + idfx::cout << "Braginskii Viscosity: uses the " + "monotonized central slope limiter." << std::endl; + } else if(haveVanLeer) { idfx::cout << "Braginskii Viscosity: uses the van Leer slope limiter." << std::endl; - } - else { + } else { IDEFIX_ERROR("Unknown slope limiter for braginskii viscosity."); } } From e840e2f938f05f64ad153c2fd54db457844e68b8 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Tue, 21 Jan 2025 21:23:47 +0100 Subject: [PATCH 11/15] revert back to old reference commit sha to avoid merge conflicts --- reference | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/reference b/reference index 609a6016..99f339c4 160000 --- a/reference +++ b/reference @@ -1 +1 @@ -Subproject commit 609a60164461fb2d4fbeffa89ce148acf1191525 +Subproject commit 99f339c43d98274a925282cc96252c5a0fa9374f From ac104701e4900c45888eda1b162483b7e240f5f4 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Wed, 22 Jan 2025 09:36:35 +0100 Subject: [PATCH 12/15] go back to last reference before fork --- reference | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/reference b/reference index 99f339c4..b675bcea 160000 --- a/reference +++ b/reference @@ -1 +1 @@ -Subproject commit 99f339c43d98274a925282cc96252c5a0fa9374f +Subproject commit b675bceaa6aabc01dded346e2d631857f698dc76 From 0c07a8e2f309489ccc9c7e836f16271ef1783105 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Victor=20R=C3=A9ville?= Date: Wed, 29 Jan 2025 08:35:10 +0100 Subject: [PATCH 13/15] Use a std::vector of IdefixArray as input for the userdef braginskii diffusion --- doc/source/modules/braginskii.rst | 5 +- src/fluid/braginskii/bragThermalDiffusion.cpp | 13 +--- src/fluid/braginskii/bragThermalDiffusion.hpp | 65 ++++++++++--------- src/fluid/fluid_defs.hpp | 8 +-- test/MHD/MTI/setup.cpp | 6 +- test/MHD/clessTDiffusion/CMakeLists.txt | 1 - test/MHD/clessTDiffusion/python/testidefix.py | 2 +- test/MHD/clessTDiffusion/setup.cpp | 26 +++++--- 8 files changed, 64 insertions(+), 62 deletions(-) diff --git a/doc/source/modules/braginskii.rst b/doc/source/modules/braginskii.rst index d9cb4f2d..fd5f2b8f 100644 --- a/doc/source/modules/braginskii.rst +++ b/doc/source/modules/braginskii.rst @@ -76,7 +76,7 @@ of the Braginskii heat flux and viscosity. Saturation with collisionless heat flux --------------------------------------- -The ``Braginskii`` module can include a collisionless saturation of the Braginskii heat flux. +The ``Braginskii`` module can include a collisionless saturation of the Braginskii heat flux, typically due to supra-thermal electrons. The heat flux is then computed as follows: :math:`q = \alpha (q_B + q_\perp) + (1-\alpha)\beta*p*v`, @@ -91,8 +91,7 @@ and :math:`\beta` controls the amplitude of the collisionless heat flux (typical .. note:: If selected, slope limiters are also used in the collisionless flux, where an upwind scheme has been implemented for stability. .. note:: - This saturation has been thought to be used mostly using the userdef function that takes four arrays as input and is enrolled through - ``data.hydro->bragThermalDiffusion->EnrollClessThermalDiffusivity()`` + This saturation has been thought to be used mostly using the userdef function that takes four userdef arrays as input. Main parameters of the module ----------------------------- diff --git a/src/fluid/braginskii/bragThermalDiffusion.cpp b/src/fluid/braginskii/bragThermalDiffusion.cpp index 0dd0b033..f251d0b0 100644 --- a/src/fluid/braginskii/bragThermalDiffusion.cpp +++ b/src/fluid/braginskii/bragThermalDiffusion.cpp @@ -19,7 +19,6 @@ #include "eos.hpp" - void BragThermalDiffusion::ShowConfig() { if(status.status==Constant) { idfx::cout << "Braginskii Thermal Diffusion: ENABLED with constant diffusivity kpar=" @@ -27,7 +26,7 @@ void BragThermalDiffusion::ShowConfig() { } else if (status.status==UserDefFunction) { idfx::cout << "Braginskii Thermal Diffusion: ENABLED with user-defined diffusivity function." << std::endl; - if(!bragDiffusivityFunc && !clessDiffusivityFunc) { + if(!bragDiffusivityFunc) { IDEFIX_ERROR("No braginskii thermal diffusion function has been enrolled"); } } else { @@ -57,7 +56,7 @@ void BragThermalDiffusion::ShowConfig() { } } -void BragThermalDiffusion::EnrollBragThermalDiffusivity(TwoArrayDiffusivityFunc myFunc) { +void BragThermalDiffusion::EnrollBragThermalDiffusivity(BragDiffusivityFunc myFunc) { if(this->status.status != UserDefFunction) { IDEFIX_WARNING("Braginskii thermal diffusivity enrollment requires Hydro/BragThermalDiffusion " "to be set to userdef in .ini file"); @@ -65,14 +64,6 @@ void BragThermalDiffusion::EnrollBragThermalDiffusivity(TwoArrayDiffusivityFunc this->bragDiffusivityFunc = myFunc; } -void BragThermalDiffusion::EnrollClessThermalDiffusivity(FourArrayDiffusivityFunc myFunc) { - if(this->status.status != UserDefFunction) { - IDEFIX_WARNING("Collisionless/Braginskii thermal diffusivity enrollment requires " - "Hydro/BragThermalDiffusion to be set to userdef in .ini file"); - } - this->clessDiffusivityFunc = myFunc; -} - void BragThermalDiffusion::AddBragDiffusiveFlux(int dir, const real t, const IdefixArray4D &Flux) { idfx::pushRegion("BragThermalDiffusion::AddBragDiffusiveFlux"); diff --git a/src/fluid/braginskii/bragThermalDiffusion.hpp b/src/fluid/braginskii/bragThermalDiffusion.hpp index d05c20b4..7b1875ff 100644 --- a/src/fluid/braginskii/bragThermalDiffusion.hpp +++ b/src/fluid/braginskii/bragThermalDiffusion.hpp @@ -9,6 +9,7 @@ #define FLUID_BRAGINSKII_BRAGTHERMALDIFFUSION_HPP_ #include +#include #include "idefix.hpp" #include "input.hpp" @@ -35,8 +36,7 @@ class BragThermalDiffusion { void AddBragDiffusiveFluxLim(int, const real, const IdefixArray4D &); // Enroll user-defined thermal conductivity - void EnrollBragThermalDiffusivity(TwoArrayDiffusivityFunc); - void EnrollClessThermalDiffusivity(FourArrayDiffusivityFunc); + void EnrollBragThermalDiffusivity(BragDiffusivityFunc); IdefixArray3D heatSrc; // Source terms of the thermal operator IdefixArray3D knorArr; @@ -53,8 +53,7 @@ class BragThermalDiffusion { // status of the module ParabolicModuleStatus &status; - TwoArrayDiffusivityFunc bragDiffusivityFunc; - FourArrayDiffusivityFunc clessDiffusivityFunc; + BragDiffusivityFunc bragDiffusivityFunc; bool includeCollisionlessTD{false}; @@ -307,7 +306,8 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, if(includeCollisionlessTD == false && haveThermalDiffusion == UserDefFunction && dir == IDIR) { if(bragDiffusivityFunc) { idfx::pushRegion("UserDef::BragThermalDiffusivityFunction"); - bragDiffusivityFunc(*this->data, t, kparArr, knorArr); + std::vector> userdefArr = {kparArr, knorArr}; + bragDiffusivityFunc(*this->data, t, userdefArr); idfx::popRegion(); } else { @@ -315,9 +315,10 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, } } else if(includeCollisionlessTD == true && haveThermalDiffusion == UserDefFunction && dir == IDIR) { - if (clessDiffusivityFunc) { + if (bragDiffusivityFunc) { idfx::pushRegion("UserDef::ClessThermalDiffusivityFunction"); - clessDiffusivityFunc(*this->data, t, kparArr, knorArr, clessAlphaArr, clessBetaArr); + std::vector> userdefArr = {kparArr, knorArr, clessAlphaArr, clessBetaArr}; + bragDiffusivityFunc(*this->data, t, userdefArr); idfx::popRegion(); } else { IDEFIX_ERROR("No user-defined Braginskii/collisionless " @@ -484,14 +485,14 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, dpm = Vc(PRS,k,j,i) - Vc(PRS,k,j,i-1); /* Upwind scheme */ - if (Vc(VX1,k,j,i) > 0.0) { - Vn = Vc(VX1,k,j,i-1)+HALF_F*SL::PLMLim(dvm, dvp); - Pn = Vc(PRS,k,j,i-1)+HALF_F*SL::PLMLim(dpm, dpp); - } - else { - Vn = Vc(VX1,k,j,i)-HALF_F*SL::PLMLim(dvm, dvp); - Pn = Vc(PRS,k,j,i)-HALF_F*SL::PLMLim(dpm, dpp); - } + if (Vc(VX1,k,j,i) > 0.0) { + Vn = Vc(VX1,k,j,i-1)+HALF_F*SL::PLMLim(dvm, dvp); + Pn = Vc(PRS,k,j,i-1)+HALF_F*SL::PLMLim(dpm, dpp); + } + else { + Vn = Vc(VX1,k,j,i)-HALF_F*SL::PLMLim(dvm, dvp); + Pn = Vc(PRS,k,j,i)-HALF_F*SL::PLMLim(dpm, dpp); + } } else { Vn = HALF_F*(Vc(VX1,k,j,i) + Vc(VX1,k,j,i-1)); Pn = HALF_F*(Vc(PRS,k,j,i) + Vc(PRS,k,j,i-1)); @@ -649,14 +650,14 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, dpp = Vc(PRS,k,j+1,i) - Vc(PRS,k,j,i); dpm = Vc(PRS,k,j,i) - Vc(PRS,k,j-1,i); - /* Upwind scheme */ - if (Vc(VX2,k,j,i) > 0.0) { - Vn = Vc(VX2,k,j-1,i)+HALF_F*SL::PLMLim(dvm, dvp); - Pn = Vc(PRS,k,j-1,i)+HALF_F*SL::PLMLim(dpm, dpp); - } else { - Vn = Vc(VX2,k,j,i)-HALF_F*SL::PLMLim(dvm, dvp); - Pn = Vc(PRS,k,j,i)-HALF_F*SL::PLMLim(dpm, dpp); - } + /* Upwind scheme */ + if (Vc(VX2,k,j,i) > 0.0) { + Vn = Vc(VX2,k,j-1,i)+HALF_F*SL::PLMLim(dvm, dvp); + Pn = Vc(PRS,k,j-1,i)+HALF_F*SL::PLMLim(dpm, dpp); + } else { + Vn = Vc(VX2,k,j,i)-HALF_F*SL::PLMLim(dvm, dvp); + Pn = Vc(PRS,k,j,i)-HALF_F*SL::PLMLim(dpm, dpp); + } } else { Vn = HALF_F*(Vc(VX2,k,j-1,i) + Vc(VX2,k,j,i)); Pn = HALF_F*(Vc(PRS,k,j-1,i) + Vc(PRS,k,j,i)); @@ -772,14 +773,14 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, dpp = Vc(PRS,k+1,j,i) - Vc(PRS,k,j,i); dpm = Vc(PRS,k,j,i) - Vc(PRS,k-1,j,i); - /* Upwind scheme */ - if (Vc(VX3,k,j,i) > 0.0) { - Vn = Vc(VX3,k-1,j,i)+HALF_F*SL::PLMLim(dvp, dvm); - Pn = Vc(PRS,k-1,j,i)+HALF_F*SL::PLMLim(dpp, dpm); - } else { - Vn = Vc(VX3,k,j,i)-HALF_F*SL::PLMLim(dvp, dvm); - Pn = Vc(PRS,k,j,i)-HALF_F*SL::PLMLim(dpp, dpm); - } + /* Upwind scheme */ + if (Vc(VX3,k,j,i) > 0.0) { + Vn = Vc(VX3,k-1,j,i)+HALF_F*SL::PLMLim(dvp, dvm); + Pn = Vc(PRS,k-1,j,i)+HALF_F*SL::PLMLim(dpp, dpm); + } else { + Vn = Vc(VX3,k,j,i)-HALF_F*SL::PLMLim(dvp, dvm); + Pn = Vc(PRS,k,j,i)-HALF_F*SL::PLMLim(dpp, dpm); + } } else { Vn = HALF_F*(Vc(VX3,k-1,j,i) + Vc(VX3,k,j,i)); Pn = HALF_F*(Vc(PRS,k-1,j,i) + Vc(PRS,k,j,i)); @@ -790,8 +791,8 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, bgradT = D_EXPAND( Bi*dTi , + Bj*dTj, +Bk*dTk); Bmag = D_EXPAND( Bi*Bi , + Bj*Bj, + Bk*Bk); + // EXPAND can yield uexpected behaviour when DIMENSIONS < COMPONENTS //printf("%f , %f\n", Bmag, EXPAND(Bi*Bi, + Bj*Bj, + Bk*Bk)); - //printf("%f , %f \n", 0.5*(Vc(VX1,k,j,i) + Vc(VX1,k,j,i+1)), Vn); Bmag = sqrt(Bmag); Bmag = FMAX(1e-6*SMALL_NUMBER,Bmag); diff --git a/src/fluid/fluid_defs.hpp b/src/fluid/fluid_defs.hpp index 1d040071..3f68c9b4 100644 --- a/src/fluid/fluid_defs.hpp +++ b/src/fluid/fluid_defs.hpp @@ -8,6 +8,7 @@ #ifndef FLUID_FLUID_DEFS_HPP_ #define FLUID_FLUID_DEFS_HPP_ +#include #include "../idefix.hpp" @@ -49,12 +50,9 @@ using SrcTermFunc = void (*) (Fluid *, const real t, const real dt); using EmfBoundaryFunc = void (*) (DataBlock &, const real t); using DiffusivityFunc = void (*) (DataBlock &, const real t, IdefixArray3D &); -using TwoArrayDiffusivityFunc = void (*) (DataBlock &, const real t, - IdefixArray3D &, IdefixArray3D &); -using FourArrayDiffusivityFunc = void (*) (DataBlock &, const real t, - IdefixArray3D &, IdefixArray3D &, - IdefixArray3D &, IdefixArray3D &); +using BragDiffusivityFunc = void (*) (DataBlock &, const real, + std::vector> &); // Deprecated signatures using SrcTermFuncOld = void (*) (DataBlock &, const real t, const real dt); diff --git a/test/MHD/MTI/setup.cpp b/test/MHD/MTI/setup.cpp index e631f7b2..1608b5ec 100644 --- a/test/MHD/MTI/setup.cpp +++ b/test/MHD/MTI/setup.cpp @@ -39,9 +39,13 @@ void Potential(DataBlock& data, const real t, IdefixArray1D& x1, IdefixArr } -void MyBragThermalConductivity(DataBlock &data, const real t, IdefixArray3D &kparArr, IdefixArray3D &knorArr) { +void MyBragThermalConductivity(DataBlock &data, const real t, std::vector> &userdefArr) { IdefixArray4D Vc = data.hydro->Vc; IdefixArray1D x2 = data.x[JDIR]; + + IdefixArray3D kparArr = userdefArr[0]; + IdefixArray3D knorArr = userdefArr[1]; + real ksi = ksiGlob; idefix_for("MyThConductivity",0,data.np_tot[KDIR],0,data.np_tot[JDIR],0,data.np_tot[IDIR], KOKKOS_LAMBDA (int k, int j, int i) { diff --git a/test/MHD/clessTDiffusion/CMakeLists.txt b/test/MHD/clessTDiffusion/CMakeLists.txt index aa9abdd7..629aed2d 100644 --- a/test/MHD/clessTDiffusion/CMakeLists.txt +++ b/test/MHD/clessTDiffusion/CMakeLists.txt @@ -1,2 +1 @@ enable_idefix_property(Idefix_MHD) -set_idefix_property(Idefix_RECONSTRUCTION Parabolic) diff --git a/test/MHD/clessTDiffusion/python/testidefix.py b/test/MHD/clessTDiffusion/python/testidefix.py index 444816d1..9427aca8 100644 --- a/test/MHD/clessTDiffusion/python/testidefix.py +++ b/test/MHD/clessTDiffusion/python/testidefix.py @@ -22,7 +22,7 @@ def gamma_prime(gamma, beta): return (gamma+beta*(gamma-1))/(1+beta*(gamma-1)) success = True -eps = 1e-2 +eps = 3e-4 gamma_eff = np.gradient(prs, r)/np.gradient(rho, r)*rho/prs diff --git a/test/MHD/clessTDiffusion/setup.cpp b/test/MHD/clessTDiffusion/setup.cpp index 624e1d72..7c6e7eb2 100644 --- a/test/MHD/clessTDiffusion/setup.cpp +++ b/test/MHD/clessTDiffusion/setup.cpp @@ -10,19 +10,28 @@ real va_vescGlob; real k0_Glob; real ParkerWind(real); -void MyClessThermalConductivity(DataBlock &data, const real t, IdefixArray3D &kparArr, IdefixArray3D &knorArr, IdefixArray3D &alpha, IdefixArray3D &clessQ) { +const real kB{1.3807e-16}; +const real mp{1.6726e-24}; + +void MyClessThermalConductivity(DataBlock &data, const real t, std::vector> &userdefArr) { IdefixArray4D Vc = data.hydro->Vc; IdefixArray1D x1 = data.x[IDIR]; IdefixArray1D x2 = data.x[JDIR]; - real norm = 1.6726e-24*0.5/(udenGlob*uvelGlob*ulenGlob*1.3807e-16); - real uTemp=0.5*uvelGlob*uvelGlob*1.6726e-24/1.3807e-16; + + IdefixArray3D kparArr = userdefArr[0]; + IdefixArray3D knorArr = userdefArr[1]; + IdefixArray3D clessAlpha = userdefArr[2]; + IdefixArray3D clessBeta = userdefArr[3]; + + real norm = mp*0.5/(udenGlob*uvelGlob*ulenGlob*kB); + real uTemp=0.5*uvelGlob*uvelGlob*mp/kB; real k0 = k0_Glob*norm; - idefix_for("MyThConductivity",0,data.np_tot[KDIR],0,data.np_tot[JDIR],0,data.np_tot[IDIR], + idefix_for("MyClessThConductivity",0,data.np_tot[KDIR],0,data.np_tot[JDIR],0,data.np_tot[IDIR], KOKKOS_LAMBDA (int k, int j, int i) { kparArr(k,j,i) = k0*pow(Vc(PRS,k,j,i)/Vc(RHO,k,j,i)*uTemp,2.5); knorArr(k,j,i) = 0.; - alpha(k,j,i) = (1.0-tanh(x1(i)-10))/2; - clessQ(k,j,i) = -1.5; + clessAlpha(k,j,i) = (1.0-tanh(x1(i)-10))/2; + clessBeta(k,j,i) = -1.5; }); } @@ -39,6 +48,7 @@ void UserDefBoundary(Hydro *hydro, int dir, BoundarySide side, real t) { real rc,vc,vwind0; real cs=cs_vescGlob*sqrt(2.); real va_vesc = va_vescGlob; + real mu = va_vesc * sqrt(2.); real PonRho; PonRho = cs*cs; @@ -49,7 +59,7 @@ void UserDefBoundary(Hydro *hydro, int dir, BoundarySide side, real t) { hydro->boundary->BoundaryFor("UserDefBoundary", dir, side, KOKKOS_LAMBDA (int k, int j, int i) { real r = x1(i); - real mu = va_vesc * sqrt(2.); + Vc(RHO,k,j,i) = vwind0/(vwind0 * r * r); Vc(PRS,k,j,i) = PonRho * Vc(RHO, k, j, i); Vc(VX1,k,j,i) = vwind0; @@ -68,7 +78,7 @@ void UserDefBoundary(Hydro *hydro, int dir, BoundarySide side, real t) { Setup::Setup(Input &input, Grid &grid, DataBlock &data, Output &output) { // Set the function for userdefboundary data.hydro->EnrollUserDefBoundary(&UserDefBoundary); - data.hydro->bragThermalDiffusion->EnrollClessThermalDiffusivity(&MyClessThermalConductivity); + data.hydro->bragThermalDiffusion->EnrollBragThermalDiffusivity(&MyClessThermalConductivity); gammaGlob=input.Get("Hydro", "gamma", 0); udenGlob=input.Get("Setup", "UNIT_DENSITY",0); ulenGlob=input.Get("Setup", "UNIT_LENGTH",0); From 250610a19b14c537515ef54f09a780c39bf180f7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Victor=20R=C3=A9ville?= Date: Wed, 29 Jan 2025 08:43:52 +0100 Subject: [PATCH 14/15] corrected typo --- src/fluid/braginskii/bragThermalDiffusion.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fluid/braginskii/bragThermalDiffusion.hpp b/src/fluid/braginskii/bragThermalDiffusion.hpp index 7b1875ff..116ff9a0 100644 --- a/src/fluid/braginskii/bragThermalDiffusion.hpp +++ b/src/fluid/braginskii/bragThermalDiffusion.hpp @@ -791,7 +791,7 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, bgradT = D_EXPAND( Bi*dTi , + Bj*dTj, +Bk*dTk); Bmag = D_EXPAND( Bi*Bi , + Bj*Bj, + Bk*Bk); - // EXPAND can yield uexpected behaviour when DIMENSIONS < COMPONENTS + // EXPAND can yield unexpected behaviour when DIMENSIONS < COMPONENTS //printf("%f , %f\n", Bmag, EXPAND(Bi*Bi, + Bj*Bj, + Bk*Bk)); Bmag = sqrt(Bmag); Bmag = FMAX(1e-6*SMALL_NUMBER,Bmag); From c0fbc84fe593b1ca599fa1c89a960d2f0c81b877 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Victor=20R=C3=A9ville?= Date: Tue, 11 Feb 2025 07:07:33 +0100 Subject: [PATCH 15/15] switch for std.vector.at() in userdef TD tests --- test/MHD/MTI/setup.cpp | 4 ++-- test/MHD/clessTDiffusion/setup.cpp | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/test/MHD/MTI/setup.cpp b/test/MHD/MTI/setup.cpp index 1608b5ec..552f7646 100644 --- a/test/MHD/MTI/setup.cpp +++ b/test/MHD/MTI/setup.cpp @@ -43,8 +43,8 @@ void MyBragThermalConductivity(DataBlock &data, const real t, std::vector Vc = data.hydro->Vc; IdefixArray1D x2 = data.x[JDIR]; - IdefixArray3D kparArr = userdefArr[0]; - IdefixArray3D knorArr = userdefArr[1]; + IdefixArray3D kparArr = userdefArr.at(0); + IdefixArray3D knorArr = userdefArr.at(1); real ksi = ksiGlob; idefix_for("MyThConductivity",0,data.np_tot[KDIR],0,data.np_tot[JDIR],0,data.np_tot[IDIR], diff --git a/test/MHD/clessTDiffusion/setup.cpp b/test/MHD/clessTDiffusion/setup.cpp index 7c6e7eb2..651a81f5 100644 --- a/test/MHD/clessTDiffusion/setup.cpp +++ b/test/MHD/clessTDiffusion/setup.cpp @@ -18,10 +18,10 @@ void MyClessThermalConductivity(DataBlock &data, const real t, std::vector x1 = data.x[IDIR]; IdefixArray1D x2 = data.x[JDIR]; - IdefixArray3D kparArr = userdefArr[0]; - IdefixArray3D knorArr = userdefArr[1]; - IdefixArray3D clessAlpha = userdefArr[2]; - IdefixArray3D clessBeta = userdefArr[3]; + IdefixArray3D kparArr = userdefArr.at(0); + IdefixArray3D knorArr = userdefArr.at(1); + IdefixArray3D clessAlpha = userdefArr.at(2); + IdefixArray3D clessBeta = userdefArr.at(3); real norm = mp*0.5/(udenGlob*uvelGlob*ulenGlob*kB); real uTemp=0.5*uvelGlob*uvelGlob*mp/kB;