diff --git a/doc/source/modules/braginskii.rst b/doc/source/modules/braginskii.rst index e657cedd..fd5f2b8f 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 @@ -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 @@ -72,18 +72,68 @@ 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, 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`, + +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 userdef function that takes four userdef arrays as input. + 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 +151,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 +169,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/src/fluid/braginskii/bragThermalDiffusion.cpp b/src/fluid/braginskii/bragThermalDiffusion.cpp index 82d5dac9..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,11 +26,11 @@ 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 { - 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; @@ -42,16 +41,27 @@ 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(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 { + IDEFIX_ERROR("Unknown slope limiter for braginskii thermal diffusion."); + } + } + if(includeCollisionlessTD) { + idfx::cout << "Braginskii Thermal Diffusion: saturation" + " with collisionless flux is enabled." << std::endl; } } 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"); + "to be set to userdef in .ini file"); } - this->diffusivityFunc = myFunc; + this->bragDiffusivityFunc = myFunc; } void BragThermalDiffusion::AddBragDiffusiveFlux(int dir, const real t, diff --git a/src/fluid/braginskii/bragThermalDiffusion.hpp b/src/fluid/braginskii/bragThermalDiffusion.hpp index 1fe774a4..116ff9a0 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" @@ -40,6 +41,8 @@ class BragThermalDiffusion { IdefixArray3D heatSrc; // Source terms of the thermal operator IdefixArray3D knorArr; IdefixArray3D kparArr; + 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; @@ -50,9 +53,14 @@ class BragThermalDiffusion { // status of the module ParabolicModuleStatus &status; - BragDiffusivityFunc diffusivityFunc; + BragDiffusivityFunc bragDiffusivityFunc; + + bool includeCollisionlessTD{false}; bool haveSlopeLimiter{false}; + bool haveMonotonizedCentral{false}; + bool haveVanLeer{false}; + bool haveMinmod{false}; // helper array IdefixArray4D &Vc; @@ -61,6 +69,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; @@ -84,9 +93,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->haveMonotonizedCentral = 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 " @@ -98,21 +109,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]); + 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("wcless") == 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/collisionless 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 or userdef."); + IDEFIX_ERROR("Unknown braginskii thermal diffusion saturation in idefix.ini. " + "Can only be nosat or wcless."); } } else { IDEFIX_ERROR("I cannot create a BragThermalDiffusion object without bragTDiffusion defined" @@ -135,7 +180,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)) @@ -179,7 +224,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) @@ -210,6 +255,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; @@ -250,16 +296,33 @@ 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 clessAlphaArr = this->clessAlphaArr; + IdefixArray3D clessBetaArr = this->clessBetaArr; - if(haveThermalDiffusion == UserDefFunction && dir == IDIR) { - if(diffusivityFunc) { + if(includeCollisionlessTD == false && haveThermalDiffusion == UserDefFunction && dir == IDIR) { + if(bragDiffusivityFunc) { idfx::pushRegion("UserDef::BragThermalDiffusivityFunction"); - diffusivityFunc(*this->data, t, kparArr, knorArr); + std::vector> userdefArr = {kparArr, knorArr}; + bragDiffusivityFunc(*this->data, t, userdefArr); + idfx::popRegion(); + } + else { + IDEFIX_ERROR("No user-defined Braginskii thermal diffusion function has been enrolled"); + } + } else if(includeCollisionlessTD == true && haveThermalDiffusion == UserDefFunction + && dir == IDIR) { + if (bragDiffusivityFunc) { + idfx::pushRegion("UserDef::ClessThermalDiffusivityFunction"); + std::vector> userdefArr = {kparArr, knorArr, clessAlphaArr, clessBetaArr}; + bragDiffusivityFunc(*this->data, t, userdefArr); idfx::popRegion(); } else { - IDEFIX_ERROR("No user-defined thermal diffusion function has been enrolled"); + IDEFIX_ERROR("No user-defined Braginskii/collisionless " + "thermal diffusion function has been enrolled"); } } @@ -273,6 +336,8 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, [[maybe_unused]] real dTi, dTj, dTk, dTn; dTi = dTj = dTk = dTn = ZERO_F; + /* For the collisionless / saturated flux */ + [[maybe_unused]] real clessAlpha, clessBeta, dvp, dvm, dpp, dpm, Pn, Vn; real locdmax = 0; /////////////////////////////////////////// @@ -287,14 +352,22 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, } 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)+clessBetaArr(k,j,i)); + } } else { knor = knorConstant; kpar = kparConstant; + if(includeCollisionlessTD) { + clessAlpha=clessAlphaConst; + clessBeta=clessBetaConst; + } } - 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 +443,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 +474,30 @@ 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 */ + 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); + + /* 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 @@ -411,9 +510,17 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, } 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)+clessBetaArr(k,j,i)); + } } else { knor = knorConstant; kpar = kparConstant; + if(includeCollisionlessTD) { + clessAlpha=clessAlphaConst; + clessBeta=clessBetaConst; + } } EXPAND( Bi = BX_J; , @@ -533,6 +640,29 @@ 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 */ + 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); + + /* 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 @@ -544,9 +674,17 @@ void BragThermalDiffusion::AddBragDiffusiveFluxLim(int dir, const real t, } 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; @@ -625,11 +763,36 @@ 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 */ + if(includeCollisionlessTD) { + 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); + + /* 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)); + } + } } // 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); + // 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); @@ -637,10 +800,15 @@ 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*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; + } }); idfx::popRegion(); } diff --git a/src/fluid/braginskii/bragViscosity.cpp b/src/fluid/braginskii/bragViscosity.cpp index 6d5007b5..e53ab2c6 100644 --- a/src/fluid/braginskii/bragViscosity.cpp +++ b/src/fluid/braginskii/bragViscosity.cpp @@ -59,7 +59,14 @@ 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; diff --git a/src/fluid/fluid_defs.hpp b/src/fluid/fluid_defs.hpp index 2269fc80..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" @@ -25,7 +26,7 @@ class Fluid; // Parabolic terms can have different status -enum HydroModuleStatus {Disabled, Constant, UserDefFunction }; +enum HydroModuleStatus {Disabled, Constant, UserDefFunction}; // Structure to describe the status of parabolic modules struct ParabolicModuleStatus { @@ -49,8 +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 BragDiffusivityFunc = void (*) (DataBlock &, const real t, - 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/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/MTI/setup.cpp b/test/MHD/MTI/setup.cpp index e631f7b2..552f7646 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.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], KOKKOS_LAMBDA (int k, int j, int i) { diff --git a/test/MHD/clessTDiffusion/CMakeLists.txt b/test/MHD/clessTDiffusion/CMakeLists.txt new file mode 100644 index 00000000..629aed2d --- /dev/null +++ b/test/MHD/clessTDiffusion/CMakeLists.txt @@ -0,0 +1 @@ +enable_idefix_property(Idefix_MHD) 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..9427aca8 --- /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 = 3e-4 + +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..651a81f5 --- /dev/null +++ b/test/MHD/clessTDiffusion/setup.cpp @@ -0,0 +1,159 @@ +#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); + +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]; + + 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; + real k0 = k0_Glob*norm; + 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.; + clessAlpha(k,j,i) = (1.0-tanh(x1(i)-10))/2; + clessBeta(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 mu = va_vesc * sqrt(2.); + 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); + + 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->EnrollBragThermalDiffusivity(&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