diff --git a/doc/wiki/Adding-New-Simulations.md b/doc/wiki/Adding-New-Simulations.md index be3294b524..b8d6feee15 100644 --- a/doc/wiki/Adding-New-Simulations.md +++ b/doc/wiki/Adding-New-Simulations.md @@ -244,7 +244,7 @@ and get the field index. For example, void AddNewField_NewProblem() { if ( NewFieldIdx == Idx_Undefined ) - NewFieldIdx = AddField( "NewFieldLabel", NORMALIZE_YES, INTERP_FRAC_YES ); + NewFieldIdx = AddField( "NewFieldLabel", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, NORMALIZE_YES, INTERP_FRAC_YES ); } ``` @@ -253,12 +253,23 @@ in the output files, and the field index `NewFieldIdx` can be used to access the field data (see the next step). The check `if ( NewFieldIdx == Idx_Undefined )` is just to avoid redundant assignments to the same field index variable. - The second parameter should be set to either `NORMALIZE_YES` or `NORMALIZE_NO`. + The second parameter should be set to either `FIXUP_FLUX_YES` or `FIXUP_FLUX_NO`. +It controls whether the new field will be corrected by the fluxes across the coarse-fine boundaries +when enabling [[ OPT__FIXUP_FLUX | Runtime-Parameters:-Hydro#OPT__FIXUP_FLUX ]]. + + The third parameter should be set to either `FIXUP_REST_YES` or `FIXUP_REST_NO`. +It controls whether the new field will be corrected by the volume-weighted average of the fine-grid data +when enabling [[ OPT__FIXUP_RESTRICT | Runtime-Parameters:-Hydro#OPT__FIXUP_RESTRICT ]]. + + The fourth parameter should be set to either `FLOOR_YES` or `FLOOR_NO`. +It controls whether the new field will be floored to `TINY_NUMBER` to ensure positivity. + + The fifth parameter should be set to either `NORMALIZE_YES` or `NORMALIZE_NO`. It controls whether the new field will be renormalized by the total gas density after every update when enabling [[ OPT__NORMALIZE_PASSIVE | Runtime-Parameters:-Hydro#OPT__NORMALIZE_PASSIVE ]]. - The third parameter should be set to either `INTERP_FRAC_YES` or `INTERP_FRAC_NO`. + The sixth parameter should be set to either `INTERP_FRAC_YES` or `INTERP_FRAC_NO`. It controls whether the new field will be converted to mass fraction during interpolation when enabling [[ OPT__INT_FRAC_PASSIVE_LR | Runtime-Parameters:-Hydro#OPT__INT_FRAC_PASSIVE_LR ]]. diff --git a/include/Global.h b/include/Global.h index ac1d5eb66e..aabea2f5b6 100644 --- a/include/Global.h +++ b/include/Global.h @@ -50,7 +50,7 @@ extern int *BaseP; // table recording the IDs extern int Flu_ParaBuf; // number of parallel buffers to exchange all fluid // variables for the fluid solver and fluid refinement -extern long FixUpVar_Flux, FixUpVar_Restrict; +extern long FixUpVar_Flux, FixUpVar_Restrict, PassiveFloorMask; extern int PassiveNorm_NVar, PassiveNorm_VarIdx[NCOMP_PASSIVE]; extern int PassiveIntFrac_NVar, PassiveIntFrac_VarIdx[NCOMP_PASSIVE]; diff --git a/include/HDF5_Typedef.h b/include/HDF5_Typedef.h index 8b60e10636..edde083e3b 100644 --- a/include/HDF5_Typedef.h +++ b/include/HDF5_Typedef.h @@ -387,7 +387,7 @@ struct SymConst_t // Structure : InputPara_t // Description : Data structure for outputting the run-time parameters // -// Note : 1. All run-time parameters are loaded from the files "Input__XXX" +// Note : 1. Most of the run-time parameters are loaded from the files "Input__XXX" //------------------------------------------------------------------------------------------------------- struct InputPara_t { @@ -623,6 +623,7 @@ struct InputPara_t int Opt__FixUp_Restrict; long FixUpRestrict_Var; int Opt__CorrAfterAllSync; + long PassiveFloor_Var; int Opt__NormalizePassive; int NormalizePassive_NVar; int NormalizePassive_VarIdx[NCOMP_PASSIVE]; diff --git a/include/Prototype.h b/include/Prototype.h index 2bf9474a87..767fd31d06 100644 --- a/include/Prototype.h +++ b/include/Prototype.h @@ -103,7 +103,7 @@ void CPU_FluidSolver( real h_Flu_Array_In[][FLU_NIN][ CUBE(FLU_NXT) ], const real ELBDM_Eta, real ELBDM_Taylor3_Coeff, const bool ELBDM_Taylor3_Auto, const double Time, const bool UsePot, const OptExtAcc_t ExtAcc, const MicroPhy_t MicroPhy, const real MinDens, const real MinPres, const real MinEint, - const real DualEnergySwitch, + const real DualEnergySwitch, const long PassiveFloor, const bool NormPassive, const int NNorm, const int NormIdx[], const bool FracPassive, const int NFrac, const int FracIdx[], const bool JeansMinPres, const real JeansMinPres_Coeff, @@ -111,23 +111,23 @@ void CPU_FluidSolver( real h_Flu_Array_In[][FLU_NIN][ CUBE(FLU_NXT) ], void Hydro_NormalizePassive( const real GasDens, real Passive[], const int NNorm, const int NormIdx[] ); #if ( MODEL == HYDRO ) real Hydro_Con2Pres( const real Dens, const real MomX, const real MomY, const real MomZ, const real Engy, - const real Passive[], const bool CheckMinPres, const real MinPres, const real Emag, + const real Passive[], const bool CheckMinPres, const real MinPres, const long PassiveFloor, const real Emag, const EoS_DE2P_t EoS_DensEint2Pres, const EoS_GUESS_t EoS_GuessHTilde, const EoS_H2TEM_t EoS_HTilde2Temp, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], const real *const EoS_Table[EOS_NTABLE_MAX], real *EintOut ); real Hydro_Con2Eint( const real Dens, const real MomX, const real MomY, const real MomZ, const real Engy, - const bool CheckMinEint, const real MinEint, const real Emag, + const bool CheckMinEint, const real MinEint, const long PassiveFloor, const real Emag, const EoS_GUESS_t EoS_GuessHTilde, const EoS_H2TEM_t EoS_HTilde2Temp, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], const real *const EoS_Table[EOS_NTABLE_MAX] ); real Hydro_ConEint2Etot( const real Dens, const real MomX, const real MomY, const real MomZ, const real Eint, const real Emag ); real Hydro_Con2Temp( const real Dens, const real MomX, const real MomY, const real MomZ, const real Engy, - const real Passive[], const bool CheckMinTemp, const real MinTemp, const real Emag, + const real Passive[], const bool CheckMinTemp, const real MinTemp, const long PassiveFloor, const real Emag, const EoS_DE2T_t EoS_DensEint2Temp, const EoS_GUESS_t EoS_GuessHTilde, const EoS_H2TEM_t EoS_HTilde2Temp, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], const real *const EoS_Table[EOS_NTABLE_MAX] ); real Hydro_Con2Entr( const real Dens, const real MomX, const real MomY, const real MomZ, const real Engy, - const real Passive[], const bool CheckMinEntr, const real MinEntr, const real Emag, + const real Passive[], const bool CheckMinEntr, const real MinEntr, const long PassiveFloor, const real Emag, const EoS_DE2S_t EoS_DensEint2Entr, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], const real *const EoS_Table[EOS_NTABLE_MAX] ); real Hydro_CheckMinPres( const real InPres, const real MinPres ); @@ -135,23 +135,23 @@ real Hydro_CheckMinEint( const real InEint, const real MinEint ); real Hydro_CheckMinTemp( const real InTemp, const real MinTemp ); real Hydro_CheckMinEntr( const real InEntr, const real MinEntr ); real Hydro_CheckMinEintInEngy( const real Dens, const real MomX, const real MomY, const real MomZ, const real InEngy, - const real MinEint, const real Emag ); + const real MinEint, const long PassiveFloor, const real Emag ); bool Hydro_IsUnphysical( const IsUnphyMode_t Mode, const real Fields[], const real Emag, const EoS_DE2P_t EoS_DensEint2Pres, const EoS_GUESS_t EoS_GuessHTilde, const EoS_H2TEM_t EoS_HTilde2Temp, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], - const real *const EoS_Table[EOS_NTABLE_MAX], + const real *const EoS_Table[EOS_NTABLE_MAX], const long PassiveFloor, const char File[], const int Line, const char Function[], const IsUnphVerb_t Verbose ); bool Hydro_IsUnphysical_Single( const real Field, const char SingleFieldName[], const real Min, const real Max, const char File[], const int Line, const char Function[], const IsUnphVerb_t Verbose ); #ifdef DUAL_ENERGY void Hydro_DualEnergyFix( const real Dens, const real MomX, const real MomY, const real MomZ, real &Etot, real &Dual, char &DE_Status, const real Gamma_m1, const real _Gamma_m1, - const bool CheckMinPres, const real MinPres, const real DualEnergySwitch, + const bool CheckMinPres, const real MinPres, const long PassiveFloor, const real DualEnergySwitch, const real Emag ); real Hydro_Con2Dual( const real Dens, const real MomX, const real MomY, const real MomZ, const real Engy, const real Emag, const EoS_DE2P_t EoS_DensEint2Pres, const double EoS_AuxArray_Flt[], - const int EoS_AuxArray_Int[], const real *const EoS_Table[EOS_NTABLE_MAX] ); + const int EoS_AuxArray_Int[], const real *const EoS_Table[EOS_NTABLE_MAX], const long PassiveFloor ); real Hydro_DensPres2Dual( const real Dens, const real Pres, const real Gamma_m1 ); real Hydro_DensDual2Pres( const real Dens, const real Dual, const real Gamma_m1, const bool CheckMinPres, const real MinPres ); @@ -250,7 +250,7 @@ void Init_ByFile(); void Init_UniformGrid( const int lv, const bool FindHomePatchForPar ); void Init_Field(); FieldIdx_t AddField( const char *InputLabel, const FixUpFlux_t FixUp_Flux, const FixUpRestrict_t FixUp_Restrict, - const NormPassive_t Norm, const IntFracPassive_t IntFrac ); + const FloorPassive_t Floor, const NormPassive_t Norm, const IntFracPassive_t IntFrac ); FieldIdx_t GetFieldIndex( const char *InputLabel, const Check_t Check ); #ifdef OPENMP void Init_OpenMP(); @@ -313,7 +313,7 @@ void dt_Close( const real h_dt_Array_T[], const int NPG ); void CPU_dtSolver( const Solver_t TSolver, real dt_Array[], const real Flu_Array[][FLU_NIN_T][ CUBE(PS1) ], const real Mag_Array[][NCOMP_MAG][ PS1P1*SQR(PS1) ], const real Pot_Array[][ CUBE(GRA_NXT) ], const double Corner_Array[][3], const int NPatchGroup, const real dh, const real Safety, - const MicroPhy_t MicroPhy, const real MinPres, const bool P5_Gradient, + const MicroPhy_t MicroPhy, const real MinPres, const long PassiveFloor, const bool P5_Gradient, const bool UsePot, const OptExtAcc_t ExtAcc, const double TargetTime ); @@ -528,7 +528,7 @@ void Hydro_BoundaryCondition_Diode( real *Array, const int BC_Face, const int NV const int ArraySizeX, const int ArraySizeY, const int ArraySizeZ, const int Idx_Start[], const int Idx_End[], const int TFluVarIdxList[], const int NVar_Der, const long TDerVarList[] ); -void Hydro_Con2Pri( const real In[], real Out[], const real MinPres, +void Hydro_Con2Pri( const real In[], real Out[], const real MinPres, const long PassiveFloor, const bool FracPassive, const int NFrac, const int FracIdx[], const bool JeansMinPres, const real JeansMinPres_Coeff, const EoS_DE2P_t EoS_DensEint2Pres, const EoS_DP2E_t EoS_DensPres2Eint, @@ -641,6 +641,7 @@ void CUAPI_Asyn_FluidSolver( real h_Flu_Array_In[][FLU_NIN ][ CUBE(FLU_NXT) ], const double Time, const bool UsePot, const OptExtAcc_t ExtAcc, const MicroPhy_t MicroPhy, const real MinDens, const real MinPres, const real MinEint, const real DualEnergySwitch, + const long PassiveFloor, const bool NormPassive, const int NNorm, const bool FracPassive, const int NFrac, const bool JeansMinPres, const real JeansMinPres_Coeff, @@ -648,7 +649,8 @@ void CUAPI_Asyn_FluidSolver( real h_Flu_Array_In[][FLU_NIN ][ CUBE(FLU_NXT) ], void CUAPI_Asyn_dtSolver( const Solver_t TSolver, real h_dt_Array[], const real h_Flu_Array[][FLU_NIN_T][ CUBE(PS1) ], const real h_Mag_Array[][NCOMP_MAG][ PS1P1*SQR(PS1) ], const real h_Pot_Array[][ CUBE(GRA_NXT) ], const double h_Corner_Array[][3], const int NPatchGroup, const real dh, const real Safety, - const MicroPhy_t MicroPhy, const real MinPres, const bool P5_Gradient, const bool UsePot, + const MicroPhy_t MicroPhy, const real MinPres, const long PassiveFloor, + const bool P5_Gradient, const bool UsePot, const OptExtAcc_t ExtAcc, const double TargetTime, const int GPU_NStream ); void CUAPI_Asyn_SrcSolver( const real h_Flu_Array_In [][FLU_NIN_S ][ CUBE(SRC_NXT) ], real h_Flu_Array_Out[][FLU_NOUT_S][ CUBE(PS1) ], @@ -656,7 +658,7 @@ void CUAPI_Asyn_SrcSolver( const real h_Flu_Array_In [][FLU_NIN_S ][ CUBE(SRC_NX const double h_Corner_Array[][3], const SrcTerms_t SrcTerms, const int NPatchGroup, const real dt, const real dh, const double TimeNew, const double TimeOld, - const real MinDens, const real MinPres, const real MinEint, + const real MinDens, const real MinPres, const real MinEint, const long PassiveFloor, const int GPU_NStream ); void CUAPI_DiagnoseDevice(); void CUAPI_MemAllocate(); @@ -807,7 +809,7 @@ void CPU_SrcSolver( const real h_Flu_Array_In [][FLU_NIN_S ][ CUBE(SRC_NXT) const double h_Corner_Array[][3], const SrcTerms_t SrcTerms, const int NPatchGroup, const real dt, const real dh, const double TimeNew, const double TimeOld, - const real MinDens, const real MinPres, const real MinEint ); + const real MinDens, const real MinPres, const real MinEint, const long PassiveFloor ); // Grackle diff --git a/include/SrcTerms.h b/include/SrcTerms.h index 2fa74df839..173e90c514 100644 --- a/include/SrcTerms.h +++ b/include/SrcTerms.h @@ -14,7 +14,7 @@ typedef void (*SrcFunc_t)( real fluid[], const real B[], const SrcTerms_t *SrcTerms, const real dt, const real dh, const double x, const double y, const double z, const double TimeNew, const double TimeOld, - const real MinDens, const real MinPres, const real MinEint, + const real MinDens, const real MinPres, const real MinEint, const long PassiveFloor, const EoS_t *EoS, const double AuxArray_Flt[], const int AuxArray_Int[] ); diff --git a/include/Typedef.h b/include/Typedef.h index fe1670b400..c2d767a47e 100644 --- a/include/Typedef.h +++ b/include/Typedef.h @@ -491,6 +491,12 @@ const FixUpRestrict_t FIXUP_REST_NO = 0, FIXUP_REST_YES = 1; +typedef int FloorPassive_t; +const FloorPassive_t + FLOOR_NULL = -1, + FLOOR_NO = 0, + FLOOR_YES = 1; + typedef int NormPassive_t; const NormPassive_t NORMALIZE_NO = 0, diff --git a/src/Auxiliary/Aux_Check_Conservation.cpp b/src/Auxiliary/Aux_Check_Conservation.cpp index 293eb1176f..850090e7b8 100644 --- a/src/Auxiliary/Aux_Check_Conservation.cpp +++ b/src/Auxiliary/Aux_Check_Conservation.cpp @@ -234,7 +234,7 @@ void Aux_Check_Conservation( const char *comment ) # endif # ifndef SRHD // Hydro_Con2Eint() calculates Eint for both HD and SRHD but we disable SRHD for now - Eint = Hydro_Con2Eint( Dens, MomX, MomY, MomZ, Etot, CheckMinEint_No, NULL_REAL, Emag, + Eint = Hydro_Con2Eint( Dens, MomX, MomY, MomZ, Etot, CheckMinEint_No, NULL_REAL, PassiveFloorMask, Emag, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); # else @@ -253,7 +253,7 @@ void Aux_Check_Conservation( const char *comment ) Cons[3] = MomZ; Cons[4] = Etot; for ( int v = NCOMP_FLUID; v < NCOMP_TOTAL; v++ ) Cons[v] = 0.0; - Hydro_Con2Pri( Cons, Prim, (real)-HUGE_NUMBER, NULL_BOOL, NULL_INT, NULL, + Hydro_Con2Pri( Cons, Prim, (real)-HUGE_NUMBER, PassiveFloorMask, NULL_BOOL, NULL_INT, NULL, NULL_BOOL, NULL_REAL, EoS_DensEint2Pres_CPUPtr, EoS_DensPres2Eint_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, NULL, &Lrtz ); HTilde = Hydro_Con2HTilde( Cons, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, diff --git a/src/Auxiliary/Aux_Check_Finite.cpp b/src/Auxiliary/Aux_Check_Finite.cpp index 70b88ad104..1b4240e625 100644 --- a/src/Auxiliary/Aux_Check_Finite.cpp +++ b/src/Auxiliary/Aux_Check_Finite.cpp @@ -63,7 +63,7 @@ void Aux_Check_Finite( const int lv, const char *comment ) const real Emag = NULL_REAL; # endif const real Pres = Hydro_Con2Pres( Data[DENS], Data[MOMX], Data[MOMY], Data[MOMZ], Data[ENGY], Data+NCOMP_FLUID, - false, NULL_REAL, Emag, + false, NULL_REAL, PassiveFloorMask, Emag, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, NULL ); diff --git a/src/Auxiliary/Aux_TakeNote.cpp b/src/Auxiliary/Aux_TakeNote.cpp index 97feb0a839..ef833f01fb 100644 --- a/src/Auxiliary/Aux_TakeNote.cpp +++ b/src/Auxiliary/Aux_TakeNote.cpp @@ -1299,6 +1299,15 @@ void Aux_TakeNote() fprintf( Note, "\n" ); } fprintf( Note, "OPT__CORR_AFTER_ALL_SYNC % d\n", OPT__CORR_AFTER_ALL_SYNC ); + fprintf( Note, "Passive_Floor_Off % d\n", -1 ); + +// target passive scalars to NOT be applied floor operations + fprintf( Note, " Target fields " ); + for (int v=0; v Should be set to the global variable "PassiveFloorMask" // NormPassive : true --> normalize passive scalars so that the sum of their mass density // is equal to the gas mass density // NNorm : Number of passive scalars to be normalized @@ -221,7 +225,7 @@ void CPU_FluidSolver( real h_Flu_Array_In[][FLU_NIN][ CUBE(FLU_NXT) ], const real ELBDM_Eta, real ELBDM_Taylor3_Coeff, const bool ELBDM_Taylor3_Auto, const double Time, const bool UsePot, const OptExtAcc_t ExtAcc, const MicroPhy_t MicroPhy, const real MinDens, const real MinPres, const real MinEint, - const real DualEnergySwitch, + const real DualEnergySwitch, const long PassiveFloor, const bool NormPassive, const int NNorm, const int NormIdx[], const bool FracPassive, const int NFrac, const int FracIdx[], const bool JeansMinPres, const real JeansMinPres_Coeff, @@ -249,7 +253,7 @@ void CPU_FluidSolver( real h_Flu_Array_In[][FLU_NIN][ CUBE(FLU_NXT) ], # if ( FLU_SCHEME == RTVD ) CPU_FluidSolver_RTVD( h_Flu_Array_In, h_Flu_Array_Out, h_Flux_Array, h_Corner_Array, h_Pot_Array_USG, - NPatchGroup, dt, dh, StoreFlux, XYZ, MinDens, MinPres, MinEint, EoS ); + NPatchGroup, dt, dh, StoreFlux, XYZ, MinDens, MinPres, MinEint, PassiveFloor, EoS ); # elif ( FLU_SCHEME == MHM || FLU_SCHEME == MHM_RP ) @@ -258,7 +262,7 @@ void CPU_FluidSolver( real h_Flu_Array_In[][FLU_NIN][ CUBE(FLU_NXT) ], h_PriVar, h_Slope_PPM, h_FC_Var, h_FC_Flux, h_FC_Mag_Half, h_EC_Ele, NPatchGroup, dt, dh, StoreFlux, StoreElectric, LR_Limiter, MinMod_Coeff, MinMod_MaxIter, Time, UsePot, ExtAcc, CPUExtAcc_Ptr, ExtAcc_AuxArray, MinDens, MinPres, MinEint, - DualEnergySwitch, NormPassive, NNorm, NormIdx, FracPassive, NFrac, FracIdx, + DualEnergySwitch, PassiveFloor, NormPassive, NNorm, NormIdx, FracPassive, NFrac, FracIdx, JeansMinPres, JeansMinPres_Coeff, EoS, MicroPhy ); # elif ( FLU_SCHEME == CTU ) @@ -268,7 +272,7 @@ void CPU_FluidSolver( real h_Flu_Array_In[][FLU_NIN][ CUBE(FLU_NXT) ], h_PriVar, h_Slope_PPM, h_FC_Var, h_FC_Flux, h_FC_Mag_Half, h_EC_Ele, NPatchGroup, dt, dh, StoreFlux, StoreElectric, LR_Limiter, MinMod_Coeff, Time, UsePot, ExtAcc, CPUExtAcc_Ptr, ExtAcc_AuxArray, MinDens, MinPres, MinEint, - DualEnergySwitch, NormPassive, NNorm, NormIdx, FracPassive, NFrac, FracIdx, + DualEnergySwitch, PassiveFloor, NormPassive, NNorm, NormIdx, FracPassive, NFrac, FracIdx, JeansMinPres, JeansMinPres_Coeff, EoS ); # else diff --git a/src/Fluid/Flu_BoundaryCondition_User.cpp b/src/Fluid/Flu_BoundaryCondition_User.cpp index fe93649785..ab6c1f4b00 100644 --- a/src/Fluid/Flu_BoundaryCondition_User.cpp +++ b/src/Fluid/Flu_BoundaryCondition_User.cpp @@ -226,14 +226,14 @@ void Flu_BoundaryCondition_User( real *Array, const int NVar_Flu, const int Ghos if ( PrepVz ) Array3D[ v2 ++ ][k][j][i] = BVal[MOMZ] / BVal[DENS]; if ( PrepPres ) Array3D[ v2 ++ ][k][j][i] = Hydro_Con2Pres( BVal[DENS], BVal[MOMX], BVal[MOMY], BVal[MOMZ], BVal[ENGY], BVal+NCOMP_FLUID, - CheckMinPres_Yes, MIN_PRES, Emag, + CheckMinPres_Yes, MIN_PRES, PassiveFloorMask, Emag, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, NULL ); if ( PrepTemp ) Array3D[ v2 ++ ][k][j][i] = Hydro_Con2Temp( BVal[DENS], BVal[MOMX], BVal[MOMY], BVal[MOMZ], BVal[ENGY], BVal+NCOMP_FLUID, - CheckMinTemp_Yes, MIN_TEMP, Emag, + CheckMinTemp_Yes, MIN_TEMP, PassiveFloorMask, Emag, EoS_DensEint2Temp_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, diff --git a/src/Fluid/Flu_Close.cpp b/src/Fluid/Flu_Close.cpp index 7cc0cd81bd..87540b8357 100644 --- a/src/Fluid/Flu_Close.cpp +++ b/src/Fluid/Flu_Close.cpp @@ -33,24 +33,24 @@ void CorrectElectric( const int SonLv, const real h_Ele_Array[][9][NCOMP_ELE][ P void ResetLongB( real L[], real R[], const real FC_B, const int d ); #endif extern void Hydro_RiemannSolver_Roe ( const int XYZ, real Flux_Out[], const real L_In[], const real R_In[], - const real MinDens, const real MinPres, const EoS_DE2P_t EoS_DensEint2Pres, + const real MinDens, const real MinPres, const long PassiveFloor, const EoS_DE2P_t EoS_DensEint2Pres, const EoS_DP2C_t EoS_DensPres2CSqr, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], const real* const EoS_Table[EOS_NTABLE_MAX] ); extern void Hydro_RiemannSolver_HLLC( const int XYZ, real Flux_Out[], const real L_In[], const real R_In[], - const real MinDens, const real MinPres, const EoS_DE2P_t EoS_DensEint2Pres, + const real MinDens, const real MinPres, const long PassiveFloor, const EoS_DE2P_t EoS_DensEint2Pres, const EoS_DP2C_t EoS_DensPres2CSqr, const EoS_GUESS_t EoS_GuessHTilde, const EoS_H2TEM_t EoS_HTilde2Temp, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], const real* const EoS_Table[EOS_NTABLE_MAX] ); extern void Hydro_RiemannSolver_HLLE( const int XYZ, real Flux_Out[], const real L_In[], const real R_In[], - const real MinDens, const real MinPres, const EoS_DE2P_t EoS_DensEint2Pres, + const real MinDens, const real MinPres, const long PassiveFloor, const EoS_DE2P_t EoS_DensEint2Pres, const EoS_DP2C_t EoS_DensPres2CSqr, const EoS_GUESS_t EoS_GuessHTilde, const EoS_H2TEM_t EoS_HTilde2Temp, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], const real* const EoS_Table[EOS_NTABLE_MAX] ); #ifdef MHD extern void Hydro_RiemannSolver_HLLD( const int XYZ, real Flux_Out[], const real L_In[], const real R_In[], - const real MinDens, const real MinPres, const EoS_DE2P_t EoS_DensEint2Pres, + const real MinDens, const real MinPres, const long PassiveFloor, const EoS_DE2P_t EoS_DensEint2Pres, const EoS_DP2C_t EoS_DensPres2CSqr, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], const real* const EoS_Table[EOS_NTABLE_MAX] ); #endif @@ -453,7 +453,7 @@ bool Unphysical( const real Fluid[], const int CheckMode, const real Emag ) # else // without DUAL_ENERGY const real Eint = Hydro_Con2Eint( Fluid[DENS], Fluid[MOMX], Fluid[MOMY], Fluid[MOMZ], Fluid[ENGY], - NoFloor, NULL_REAL, Emag, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, + NoFloor, NULL_REAL, PassiveFloorMask, Emag, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); if ( Eint < (real)MIN_EINT || !Aux_IsFinite(Eint) ) return true; @@ -464,7 +464,7 @@ bool Unphysical( const real Fluid[], const int CheckMode, const real Emag ) { const real Pres = Hydro_Con2Pres( Fluid[DENS], Fluid[MOMX], Fluid[MOMY], Fluid[MOMZ], Fluid[ENGY], Fluid+NCOMP_FLUID, - CheckMinPres_No, NULL_REAL, Emag, + CheckMinPres_No, NULL_REAL, PassiveFloorMask, Emag, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, @@ -678,6 +678,7 @@ void CorrectUnphysical( const int lv, const int NPG, const int *PID0_List, ResetLongB( VarL[d], VarC, FC_B[0], d ); // reset the longitudinal B field # endif Hydro_RiemannSolver_Roe ( d, FluxL[d], VarL[d], VarC, MIN_DENS, MIN_PRES, + PassiveFloorMask, EoS_DensEint2Pres_CPUPtr, EoS_DensPres2CSqr_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); @@ -685,6 +686,7 @@ void CorrectUnphysical( const int lv, const int NPG, const int *PID0_List, ResetLongB( VarC, VarR[d], FC_B[1], d ); // reset the longitudinal B field # endif Hydro_RiemannSolver_Roe ( d, FluxR[d], VarC, VarR[d], MIN_DENS, MIN_PRES, + PassiveFloorMask, EoS_DensEint2Pres_CPUPtr, EoS_DensPres2CSqr_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); @@ -698,10 +700,12 @@ void CorrectUnphysical( const int lv, const int NPG, const int *PID0_List, # ifndef MHD case RSOLVER_1ST_HLLC: Hydro_RiemannSolver_HLLC( d, FluxL[d], VarL[d], VarC, MIN_DENS, MIN_PRES, + PassiveFloorMask, EoS_DensEint2Pres_CPUPtr, EoS_DensPres2CSqr_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); Hydro_RiemannSolver_HLLC( d, FluxR[d], VarC, VarR[d], MIN_DENS, MIN_PRES, + PassiveFloorMask, EoS_DensEint2Pres_CPUPtr, EoS_DensPres2CSqr_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); @@ -713,6 +717,7 @@ void CorrectUnphysical( const int lv, const int NPG, const int *PID0_List, ResetLongB( VarL[d], VarC, FC_B[0], d ); // reset the longitudinal B field # endif Hydro_RiemannSolver_HLLE( d, FluxL[d], VarL[d], VarC, MIN_DENS, MIN_PRES, + PassiveFloorMask, EoS_DensEint2Pres_CPUPtr, EoS_DensPres2CSqr_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); @@ -721,6 +726,7 @@ void CorrectUnphysical( const int lv, const int NPG, const int *PID0_List, ResetLongB( VarC, VarR[d], FC_B[1], d ); // reset the longitudinal B field # endif Hydro_RiemannSolver_HLLE( d, FluxR[d], VarC, VarR[d], MIN_DENS, MIN_PRES, + PassiveFloorMask, EoS_DensEint2Pres_CPUPtr, EoS_DensPres2CSqr_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); @@ -738,6 +744,7 @@ void CorrectUnphysical( const int lv, const int NPG, const int *PID0_List, ResetLongB( VarL[d], VarC, FC_B[0], d ); // reset the longitudinal B field # endif Hydro_RiemannSolver_HLLD( d, FluxL[d], VarL[d], VarC, MIN_DENS, MIN_PRES, + PassiveFloorMask, EoS_DensEint2Pres_CPUPtr, EoS_DensPres2CSqr_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); @@ -745,6 +752,7 @@ void CorrectUnphysical( const int lv, const int NPG, const int *PID0_List, ResetLongB( VarC, VarR[d], FC_B[1], d ); // reset the longitudinal B field # endif Hydro_RiemannSolver_HLLD( d, FluxR[d], VarC, VarR[d], MIN_DENS, MIN_PRES, + PassiveFloorMask, EoS_DensEint2Pres_CPUPtr, EoS_DensPres2CSqr_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); @@ -790,7 +798,7 @@ void CorrectUnphysical( const int lv, const int NPG, const int *PID0_List, # ifdef DUAL_ENERGY Hydro_DualEnergyFix( Update[DENS], Update[MOMX], Update[MOMY], Update[MOMZ], Update[ENGY], Update[DUAL], h_DE_Array_F_Out[TID][idx_out], EoS_AuxArray_Flt[1], EoS_AuxArray_Flt[2], - CorrPres_No, NULL_REAL, DUAL_ENERGY_SWITCH, Emag_Out ); + CorrPres_No, NULL_REAL, PassiveFloorMask, DUAL_ENERGY_SWITCH, Emag_Out ); # endif if ( Unphysical(Update, CheckMinEint, Emag_Out) ) @@ -824,9 +832,11 @@ void CorrectUnphysical( const int lv, const int NPG, const int *PID0_List, { case RSOLVER_1ST_ROE: Hydro_RiemannSolver_Roe ( d, FluxL_1D, Corr1D_InOut_PtrL, Corr1D_InOut_PtrC, MIN_DENS, MIN_PRES, + PassiveFloorMask, EoS_DensEint2Pres_CPUPtr, EoS_DensPres2CSqr_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); Hydro_RiemannSolver_Roe ( d, FluxR_1D, Corr1D_InOut_PtrC, Corr1D_InOut_PtrR, MIN_DENS, MIN_PRES, + PassiveFloorMask, EoS_DensEint2Pres_CPUPtr, EoS_DensPres2CSqr_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); break; @@ -834,10 +844,12 @@ void CorrectUnphysical( const int lv, const int NPG, const int *PID0_List, # ifndef MHD case RSOLVER_1ST_HLLC: Hydro_RiemannSolver_HLLC( d, FluxL_1D, Corr1D_InOut_PtrL, Corr1D_InOut_PtrC, MIN_DENS, MIN_PRES, + PassiveFloorMask, EoS_DensEint2Pres_CPUPtr, EoS_DensPres2CSqr_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); Hydro_RiemannSolver_HLLC( d, FluxR_1D, Corr1D_InOut_PtrC, Corr1D_InOut_PtrR, MIN_DENS, MIN_PRES, + PassiveFloorMask, EoS_DensEint2Pres_CPUPtr, EoS_DensPres2CSqr_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); @@ -846,10 +858,12 @@ void CorrectUnphysical( const int lv, const int NPG, const int *PID0_List, case RSOLVER_1ST_HLLE: Hydro_RiemannSolver_HLLE( d, FluxL_1D, Corr1D_InOut_PtrL, Corr1D_InOut_PtrC, MIN_DENS, MIN_PRES, + PassiveFloorMask, EoS_DensEint2Pres_CPUPtr, EoS_DensPres2CSqr_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); Hydro_RiemannSolver_HLLE( d, FluxR_1D, Corr1D_InOut_PtrC, Corr1D_InOut_PtrR, MIN_DENS, MIN_PRES, + PassiveFloorMask, EoS_DensEint2Pres_CPUPtr, EoS_DensPres2CSqr_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); @@ -860,9 +874,11 @@ void CorrectUnphysical( const int lv, const int NPG, const int *PID0_List, Aux_Error( ERROR_INFO, "RSOLVER_1ST_HLLD in MHD is NOT supported yet !!\n" ); /* Hydro_RiemannSolver_HLLD( d, FluxL_1D, Corr1D_InOut_PtrL, Corr1D_InOut_PtrC, MIN_DENS, MIN_PRES, + PassiveFloorMask, EoS_DensEint2Pres_CPUPtr, EoS_DensPres2CSqr_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); Hydro_RiemannSolver_HLLD( d, FluxR_1D, Corr1D_InOut_PtrC, Corr1D_InOut_PtrR, MIN_DENS, MIN_PRES, + PassiveFloorMask, EoS_DensEint2Pres_CPUPtr, EoS_DensPres2CSqr_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); */ @@ -919,7 +935,8 @@ void CorrectUnphysical( const int lv, const int NPG, const int *PID0_List, // floor and normalize passive scalars # if ( NCOMP_PASSIVE > 0 ) - for (int v=NCOMP_FLUID; v apply it only when AutoReduceDt_Continue is false @@ -946,7 +963,7 @@ void CorrectUnphysical( const int lv, const int NPG, const int *PID0_List, # else if ( ! AutoReduceDt_Continue && OPT__LAST_RESORT_FLOOR ) Update[ENGY] = Hydro_CheckMinEintInEngy( Update[DENS], Update[MOMX], Update[MOMY], Update[MOMZ], Update[ENGY], - MIN_EINT, Emag_Out ); + MIN_EINT, PassiveFloorMask, Emag_Out ); # endif @@ -1025,11 +1042,11 @@ void CorrectUnphysical( const int lv, const int NPG, const int *PID0_List, fprintf( File, "%14.7e, %14.7e", Hydro_Con2Eint( In[DENS], In[MOMX], In[MOMY], In[MOMZ], In[ENGY], - CheckMinEint_No, NULL_REAL, Emag_In, + CheckMinEint_No, NULL_REAL, PassiveFloorMask, Emag_In, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ), Hydro_Con2Pres( In[DENS], In[MOMX], In[MOMY], In[MOMZ], In[ENGY], In+NCOMP_FLUID, - CheckMinPres_No, NULL_REAL, Emag_In, + CheckMinPres_No, NULL_REAL, PassiveFloorMask, Emag_In, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, NULL ) ); # if ( DUAL_ENERGY == DE_ENPY ) @@ -1052,11 +1069,11 @@ void CorrectUnphysical( const int lv, const int NPG, const int *PID0_List, fprintf( File, "%14.7e, %14.7e", Hydro_Con2Eint( Out[DENS], Out[MOMX], Out[MOMY], Out[MOMZ], Out[ENGY], - CheckMinEint_No, NULL_REAL, Emag_Out, + CheckMinEint_No, NULL_REAL, PassiveFloorMask, Emag_Out, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ), Hydro_Con2Pres( Out[DENS], Out[MOMX], Out[MOMY], Out[MOMZ], Out[ENGY], Out+NCOMP_FLUID, - CheckMinPres_No, NULL_REAL, Emag_Out, + CheckMinPres_No, NULL_REAL, PassiveFloorMask, Emag_Out, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, NULL ) ); # if ( DUAL_ENERGY == DE_ENPY ) @@ -1079,11 +1096,11 @@ void CorrectUnphysical( const int lv, const int NPG, const int *PID0_List, fprintf( File, "%14.7e, %14.7e", Hydro_Con2Eint( Update[DENS], Update[MOMX], Update[MOMY], Update[MOMZ], Update[ENGY], - CheckMinEint_No, NULL_REAL, Emag_Update, + CheckMinEint_No, NULL_REAL, PassiveFloorMask, Emag_Update, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ), Hydro_Con2Pres( Update[DENS], Update[MOMX], Update[MOMY], Update[MOMZ], Update[ENGY], Update+NCOMP_FLUID, - CheckMinPres_No, NULL_REAL, Emag_Update, + CheckMinPres_No, NULL_REAL, PassiveFloorMask, Emag_Update, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, NULL ) ); # if ( DUAL_ENERGY == DE_ENPY ) @@ -1133,7 +1150,7 @@ void CorrectUnphysical( const int lv, const int NPG, const int *PID0_List, # endif fprintf( File, " %14.7e", Hydro_Con2Eint(tmp[0], tmp[1], tmp[2], tmp[3], tmp[4], - CheckMinEint_No, NULL_REAL, Emag_tmp, + CheckMinEint_No, NULL_REAL, PassiveFloorMask, Emag_tmp, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table) ); # ifdef MHD diff --git a/src/Fluid/Flu_DerivedField_BuiltIn.cpp b/src/Fluid/Flu_DerivedField_BuiltIn.cpp index 52204124bb..8d57641d49 100644 --- a/src/Fluid/Flu_DerivedField_BuiltIn.cpp +++ b/src/Fluid/Flu_DerivedField_BuiltIn.cpp @@ -175,7 +175,7 @@ void Flu_DerivedField_Mach( real Out[], const real FluIn[], const real MagIn[], # ifdef SRHD real Prim[NCOMP_TOTAL], FourCs2, U2; - Hydro_Con2Pri( fluid, Prim, MIN_PRES, NULL_BOOL, NULL_INT, NULL, + Hydro_Con2Pri( fluid, Prim, MIN_PRES, PassiveFloorMask, NULL_BOOL, NULL_INT, NULL, NULL_BOOL, NULL_REAL, EoS_DensEint2Pres_CPUPtr, EoS_DensPres2Eint_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, NULL, NULL ); @@ -194,7 +194,7 @@ void Flu_DerivedField_Mach( real Out[], const real FluIn[], const real MagIn[], V2 = SQR( Vx ) + SQR( Vy ) + SQR( Vz ); Pres = Hydro_Con2Pres( fluid[DENS], fluid[MOMX], fluid[MOMY], fluid[MOMZ], fluid[ENGY], fluid+NCOMP_FLUID, - CheckMinPres_Yes, MIN_PRES, Emag, + CheckMinPres_Yes, MIN_PRES, PassiveFloorMask, Emag, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, NULL ); Cs2 = EoS_DensPres2CSqr_CPUPtr( fluid[DENS], Pres, fluid+NCOMP_FLUID, EoS_AuxArray_Flt, EoS_AuxArray_Int, diff --git a/src/Fluid/Flu_FixUp_Flux.cpp b/src/Fluid/Flu_FixUp_Flux.cpp index 702e5de768..4e985d48dd 100644 --- a/src/Fluid/Flu_FixUp_Flux.cpp +++ b/src/Fluid/Flu_FixUp_Flux.cpp @@ -208,7 +208,7 @@ void Flu_FixUp_Flux( const int lv, const long TVar ) # endif { Pres = Hydro_Con2Pres( ForEint[DENS], ForEint[MOMX], ForEint[MOMY], ForEint[MOMZ], ForEint[ENGY], - ForEint+NCOMP_FLUID, CheckMinPres_No, NULL_REAL, Emag, + ForEint+NCOMP_FLUID, CheckMinPres_No, NULL_REAL, PassiveFloorMask, Emag, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, &Eint ); @@ -245,7 +245,7 @@ void Flu_FixUp_Flux( const int lv, const long TVar ) if ( Hydro_IsUnphysical( UNPHY_MODE_CONS, CorrVal, NULL_REAL, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, - ERROR_INFO, UNPHY_VERBOSE ) ) + PassiveFloorMask, ERROR_INFO, UNPHY_VERBOSE ) ) # else if ( CorrVal[DENS] <= MIN_DENS # ifndef BAROTROPIC_EOS @@ -289,7 +289,8 @@ void Flu_FixUp_Flux( const int lv, const long TVar ) // floor and normalize the passive scalars # if ( NCOMP_PASSIVE > 0 && MODEL == HYDRO ) for (int v=NCOMP_FLUID; vpatch[FaFluSg][FaLv][FaPID]->fluid[ENGY][k][j][i], amr->patch[FaFluSg][FaLv][FaPID]->fluid[DUAL][k][j][i], dummy, EoS_AuxArray_Flt[1], EoS_AuxArray_Flt[2], CheckMinPres_Yes, MIN_PRES, - UseDual2FixEngy, Emag ); + PassiveFloorMask, UseDual2FixEngy, Emag ); # else // #ifdef DUAL_ENERGY @@ -438,7 +438,7 @@ void Flu_FixUp_Restrict( const int FaLv, const int SonFluSg, const int FaFluSg, amr->patch[FaFluSg][FaLv][FaPID]->fluid[MOMY][k][j][i], amr->patch[FaFluSg][FaLv][FaPID]->fluid[MOMZ][k][j][i], amr->patch[FaFluSg][FaLv][FaPID]->fluid[ENGY][k][j][i], - MIN_EINT, Emag ); + MIN_EINT, PassiveFloorMask, Emag ); # endif // #ifdef DUAL_ENERGY ... else ... } // i,j,k # endif // #if ( MODEL == HYDRO ) diff --git a/src/Fluid/Flu_Prepare.cpp b/src/Fluid/Flu_Prepare.cpp index 70e605f29c..85aa100c61 100644 --- a/src/Fluid/Flu_Prepare.cpp +++ b/src/Fluid/Flu_Prepare.cpp @@ -211,7 +211,7 @@ void Flu_Prepare( const int lv, const double PrepTime, if ( Hydro_IsUnphysical( UNPHY_MODE_CONS, fluid, Emag, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, - ERROR_INFO, UNPHY_VERBOSE ) ) + PassiveFloorMask, ERROR_INFO, UNPHY_VERBOSE ) ) CheckFailed = true; // generic diff --git a/src/Fluid/Flu_ResetByUser.cpp b/src/Fluid/Flu_ResetByUser.cpp index a0b447ee39..c83beabb11 100644 --- a/src/Fluid/Flu_ResetByUser.cpp +++ b/src/Fluid/Flu_ResetByUser.cpp @@ -287,18 +287,20 @@ void Flu_ResetByUser_API_Default( const int lv, const int FluSg, const int MagSg fluid[DENS] = FMAX( fluid[DENS], (real)MIN_DENS ); # ifndef SRHD fluid[ENGY] = Hydro_CheckMinEintInEngy( fluid[DENS], fluid[MOMX], fluid[MOMY], fluid[MOMZ], fluid[ENGY], - MIN_EINT, Emag ); + MIN_EINT, PassiveFloorMask, Emag ); # endif // calculate the dual-energy variable (entropy or internal energy) # ifdef DUAL_ENERGY fluid[DUAL] = Hydro_Con2Dual( fluid[DENS], fluid[MOMX], fluid[MOMY], fluid[MOMZ], fluid[ENGY], Emag, - EoS_DensEint2Pres_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); + EoS_DensEint2Pres_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, + PassiveFloorMask ); # endif // floor and normalize passive scalars # if ( NCOMP_PASSIVE > 0 ) - for (int v=NCOMP_FLUID; v Should be set to the global variable "PassiveFloorMask" // NormPassive : true --> normalize passive scalars so that the sum of their mass density // is equal to the gas mass density // NNorm : Number of passive scalars to be normalized @@ -277,6 +281,7 @@ void CUAPI_Asyn_FluidSolver( real h_Flu_Array_In[][FLU_NIN ][ CUBE(FLU_NXT) ], const double Time, const bool UsePot, const OptExtAcc_t ExtAcc, const MicroPhy_t MicroPhy, const real MinDens, const real MinPres, const real MinEint, const real DualEnergySwitch, + const long PassiveFloor, const bool NormPassive, const int NNorm, const bool FracPassive, const int NFrac, const bool JeansMinPres, const real JeansMinPres_Coeff, @@ -545,7 +550,7 @@ void CUAPI_Asyn_FluidSolver( real h_Flu_Array_In[][FLU_NIN ][ CUBE(FLU_NXT) ], d_Flux_Array + UsedPatch[s], d_Corner_Array_F + UsedPatch[s], d_Pot_Array_USG_F + UsedPatch[s], - dt, 1.0/dh, StoreFlux, XYZ, MinDens, MinPres, MinEint, EoS ); + dt, 1.0/dh, StoreFlux, XYZ, MinDens, MinPres, MinEint, PassiveFloor, EoS ); # elif ( FLU_SCHEME == MHM || FLU_SCHEME == MHM_RP ) @@ -567,7 +572,7 @@ void CUAPI_Asyn_FluidSolver( real h_Flu_Array_In[][FLU_NIN ][ CUBE(FLU_NXT) ], d_EC_Ele + UsedPatch[s], dt, dh, StoreFlux, StoreElectric, LR_Limiter, MinMod_Coeff, MinMod_MaxIter, Time, UsePot, ExtAcc, GPUExtAcc_Ptr, MinDens, MinPres, MinEint, - DualEnergySwitch, NormPassive, NNorm, FracPassive, NFrac, + DualEnergySwitch, PassiveFloor, NormPassive, NNorm, FracPassive, NFrac, JeansMinPres, JeansMinPres_Coeff, EoS, MicroPhy ); # elif ( FLU_SCHEME == CTU ) @@ -590,7 +595,7 @@ void CUAPI_Asyn_FluidSolver( real h_Flu_Array_In[][FLU_NIN ][ CUBE(FLU_NXT) ], d_EC_Ele + UsedPatch[s], dt, dh, StoreFlux, StoreElectric, LR_Limiter, MinMod_Coeff, Time, UsePot, ExtAcc, GPUExtAcc_Ptr, MinDens, MinPres, MinEint, - DualEnergySwitch, NormPassive, NNorm, FracPassive, NFrac, + DualEnergySwitch, PassiveFloor, NormPassive, NNorm, FracPassive, NFrac, JeansMinPres, JeansMinPres_Coeff, EoS ); # else diff --git a/src/GPU_API/CUAPI_Asyn_SrcSolver.cu b/src/GPU_API/CUAPI_Asyn_SrcSolver.cu index a5a293d9e6..4130e66946 100644 --- a/src/GPU_API/CUAPI_Asyn_SrcSolver.cu +++ b/src/GPU_API/CUAPI_Asyn_SrcSolver.cu @@ -13,7 +13,7 @@ void CUSRC_SrcSolver_IterateAllCells( const double g_Corner_Array[][3], const SrcTerms_t SrcTerms, const int NPatchGroup, const real dt, const real dh, const double TimeNew, const double TimeOld, - const real MinDens, const real MinPres, const real MinEint, const EoS_t EoS ); + const real MinDens, const real MinPres, const real MinEint, const long PassiveFloor, const EoS_t EoS ); // device pointers extern real (*d_Flu_Array_S_In )[FLU_NIN_S ][ CUBE(SRC_NXT) ]; @@ -56,6 +56,7 @@ extern cudaStream_t *Stream; // TimeOld : Physical time before update // --> This function updates physical time from TimeOld to TimeNew // MinDens/Pres/Eint : Density, pressure, and internal energy floors +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored // GPU_NStream : Number of CUDA streams for the asynchronous memory copy // // Return : h_Flu_Array_Out[] @@ -66,7 +67,7 @@ void CUAPI_Asyn_SrcSolver( const real h_Flu_Array_In [][FLU_NIN_S ][ CUBE(SRC_NX const double h_Corner_Array[][3], const SrcTerms_t SrcTerms, const int NPatchGroup, const real dt, const real dh, const double TimeNew, const double TimeOld, - const real MinDens, const real MinPres, const real MinEint, + const real MinDens, const real MinPres, const real MinEint, const long PassiveFloor, const int GPU_NStream ) { @@ -161,7 +162,7 @@ void CUAPI_Asyn_SrcSolver( const real h_Flu_Array_In [][FLU_NIN_S ][ CUBE(SRC_NX d_Mag_Array_S_In + UsedPatch[s], d_Corner_Array_S + UsedPatch[s], SrcTerms, NPatchGroup, dt, dh, TimeNew, TimeOld, - MinDens, MinPres, MinEint, EoS ); + MinDens, MinPres, MinEint, PassiveFloor, EoS ); CUDA_CHECK_ERROR( cudaGetLastError() ); } // for (int s=0; sMagSg[lv] ); # endif - Eint = Hydro_Con2Eint( Dens, Px, Py, Pz, Etot, CheckMinEint_Yes, MIN_EINT, Emag, + Eint = Hydro_Con2Eint( Dens, Px, Py, Pz, Etot, CheckMinEint_Yes, MIN_EINT, PassiveFloorMask, Emag, NULL, NULL, NULL, NULL, NULL ); # endif // #ifdef DUAL_ENERGY ... else diff --git a/src/Init/Init_ByFile.cpp b/src/Init/Init_ByFile.cpp index 858c469689..caf67174fd 100644 --- a/src/Init/Init_ByFile.cpp +++ b/src/Init/Init_ByFile.cpp @@ -704,7 +704,8 @@ void Init_ByFile_Default( real fluid_out[], const real fluid_in[], const int nva # ifdef DUAL_ENERGY fluid_out[DUAL] = Hydro_Con2Dual( fluid_in[DENS], fluid_in[MOMX], fluid_in[MOMY], fluid_in[MOMZ], fluid_in[ENGY], Emag, - EoS_DensEint2Pres_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); + EoS_DensEint2Pres_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, + PassiveFloorMask ); # endif # elif ( MODEL == ELBDM ) diff --git a/src/Init/Init_ByRestart_HDF5.cpp b/src/Init/Init_ByRestart_HDF5.cpp index f3e5efa4e5..e08c16a93f 100644 --- a/src/Init/Init_ByRestart_HDF5.cpp +++ b/src/Init/Init_ByRestart_HDF5.cpp @@ -2153,6 +2153,7 @@ void Check_InputPara( const char *FileName, const int FormatVersion ) LoadField( "Opt__FixUp_Restrict", &RS.Opt__FixUp_Restrict, SID, TID, NonFatal, &RT.Opt__FixUp_Restrict, 1, NonFatal ); LoadField( "FixUpRestrict_Var", &RS.FixUpRestrict_Var, SID, TID, NonFatal, &RT.FixUpRestrict_Var, 1, NonFatal ); LoadField( "Opt__CorrAfterAllSync", &RS.Opt__CorrAfterAllSync, SID, TID, NonFatal, &RT.Opt__CorrAfterAllSync, 1, NonFatal ); + LoadField( "PassiveFloor_Var", &RS.PassiveFloor_Var, SID, TID, NonFatal, &RT.PassiveFloor_Var, 1, NonFatal ); LoadField( "Opt__NormalizePassive", &RS.Opt__NormalizePassive, SID, TID, NonFatal, &RT.Opt__NormalizePassive, 1, NonFatal ); LoadField( "NormalizePassive_NVar", &RS.NormalizePassive_NVar, SID, TID, NonFatal, &RT.NormalizePassive_NVar, 1, NonFatal ); LoadField( "NormalizePassive_VarIdx", RS.NormalizePassive_VarIdx, SID, TID, NonFatal, RT.NormalizePassive_VarIdx, NP, NonFatal ); diff --git a/src/Init/Init_Field.cpp b/src/Init/Init_Field.cpp index bba374d938..f48d1aa263 100644 --- a/src/Init/Init_Field.cpp +++ b/src/Init/Init_Field.cpp @@ -35,6 +35,7 @@ void Init_Field() NDefinedField = 0; FixUpVar_Flux = 0; FixUpVar_Restrict = 0; + PassiveFloorMask = 0; PassiveNorm_NVar = 0; PassiveIntFrac_NVar = 0; @@ -49,11 +50,11 @@ void Init_Field() // --> must not change the following order of declaration since they must be consistent // with the symbolic constants defined in Macro.h (e.g., DENS) # if ( MODEL == HYDRO ) - Idx_Dens = AddField( "Dens", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_NO, INTERP_FRAC_NO ); - Idx_MomX = AddField( "MomX", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_NO, INTERP_FRAC_NO ); - Idx_MomY = AddField( "MomY", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_NO, INTERP_FRAC_NO ); - Idx_MomZ = AddField( "MomZ", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_NO, INTERP_FRAC_NO ); - Idx_Engy = AddField( "Engy", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_NO, INTERP_FRAC_NO ); + Idx_Dens = AddField( "Dens", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_NULL, NORMALIZE_NO, INTERP_FRAC_NO ); + Idx_MomX = AddField( "MomX", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_NULL, NORMALIZE_NO, INTERP_FRAC_NO ); + Idx_MomY = AddField( "MomY", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_NULL, NORMALIZE_NO, INTERP_FRAC_NO ); + Idx_MomZ = AddField( "MomZ", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_NULL, NORMALIZE_NO, INTERP_FRAC_NO ); + Idx_Engy = AddField( "Engy", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_NULL, NORMALIZE_NO, INTERP_FRAC_NO ); if ( Idx_Dens != DENS ) Aux_Error( ERROR_INFO, "inconsistent Idx_Dens (%d != %d) !!\n", Idx_Dens, DENS ); if ( Idx_MomX != MOMX ) Aux_Error( ERROR_INFO, "inconsistent Idx_MomX (%d != %d) !!\n", Idx_MomX, MOMX ); @@ -69,18 +70,18 @@ void Init_Field() # elif ( MODEL == ELBDM ) # if ( ELBDM_SCHEME == ELBDM_WAVE ) - Idx_Dens = AddField( "Dens", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_NO, INTERP_FRAC_NO ); - Idx_Real = AddField( "Real", FIXUP_FLUX_NO, FIXUP_REST_YES, NORMALIZE_NO, INTERP_FRAC_NO ); - Idx_Imag = AddField( "Imag", FIXUP_FLUX_NO, FIXUP_REST_YES, NORMALIZE_NO, INTERP_FRAC_NO ); + Idx_Dens = AddField( "Dens", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_NULL, NORMALIZE_NO, INTERP_FRAC_NO ); + Idx_Real = AddField( "Real", FIXUP_FLUX_NO, FIXUP_REST_YES, FLOOR_NULL, NORMALIZE_NO, INTERP_FRAC_NO ); + Idx_Imag = AddField( "Imag", FIXUP_FLUX_NO, FIXUP_REST_YES, FLOOR_NULL, NORMALIZE_NO, INTERP_FRAC_NO ); if ( Idx_Dens != DENS ) Aux_Error( ERROR_INFO, "inconsistent Idx_Dens (%d != %d) !!\n", Idx_Dens, DENS ); if ( Idx_Real != REAL ) Aux_Error( ERROR_INFO, "inconsistent Idx_Real (%d != %d) !!\n", Idx_Real, REAL ); if ( Idx_Imag != IMAG ) Aux_Error( ERROR_INFO, "inconsistent Idx_Imag (%d != %d) !!\n", Idx_Imag, IMAG ); # elif ( ELBDM_SCHEME == ELBDM_HYBRID ) - Idx_Dens = AddField( "Dens", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_NO, INTERP_FRAC_NO ); - Idx_Phas = AddField( "Phase", FIXUP_FLUX_NO, FIXUP_REST_YES, NORMALIZE_NO, INTERP_FRAC_NO ); - Idx_Stub = AddField( "Stub", FIXUP_FLUX_NO, FIXUP_REST_YES, NORMALIZE_NO, INTERP_FRAC_NO ); + Idx_Dens = AddField( "Dens", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_NULL, NORMALIZE_NO, INTERP_FRAC_NO ); + Idx_Phas = AddField( "Phase", FIXUP_FLUX_NO, FIXUP_REST_YES, FLOOR_NULL, NORMALIZE_NO, INTERP_FRAC_NO ); + Idx_Stub = AddField( "Stub", FIXUP_FLUX_NO, FIXUP_REST_YES, FLOOR_NULL, NORMALIZE_NO, INTERP_FRAC_NO ); Idx_Real = Idx_Phas; Idx_Imag = Idx_Stub; @@ -101,31 +102,31 @@ void Init_Field() // 2. add other predefined fields # ifdef SUPPORT_GRACKLE if ( GRACKLE_PRIMORDIAL >= GRACKLE_PRI_CHE_NSPE6 ) { - Idx_e = AddField( "Electron", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_YES, INTERP_FRAC_YES ); - Idx_HI = AddField( "HI", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_YES, INTERP_FRAC_YES ); - Idx_HII = AddField( "HII", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_YES, INTERP_FRAC_YES ); - Idx_HeI = AddField( "HeI", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_YES, INTERP_FRAC_YES ); - Idx_HeII = AddField( "HeII", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_YES, INTERP_FRAC_YES ); - Idx_HeIII = AddField( "HeIII", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_YES, INTERP_FRAC_YES ); + Idx_e = AddField( "Electron", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, NORMALIZE_YES, INTERP_FRAC_YES ); + Idx_HI = AddField( "HI", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, NORMALIZE_YES, INTERP_FRAC_YES ); + Idx_HII = AddField( "HII", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, NORMALIZE_YES, INTERP_FRAC_YES ); + Idx_HeI = AddField( "HeI", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, NORMALIZE_YES, INTERP_FRAC_YES ); + Idx_HeII = AddField( "HeII", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, NORMALIZE_YES, INTERP_FRAC_YES ); + Idx_HeIII = AddField( "HeIII", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, NORMALIZE_YES, INTERP_FRAC_YES ); } if ( GRACKLE_PRIMORDIAL >= GRACKLE_PRI_CHE_NSPE9 ) { - Idx_HM = AddField( "HM", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_YES, INTERP_FRAC_YES ); - Idx_H2I = AddField( "H2I", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_YES, INTERP_FRAC_YES ); - Idx_H2II = AddField( "H2II", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_YES, INTERP_FRAC_YES ); + Idx_HM = AddField( "HM", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, NORMALIZE_YES, INTERP_FRAC_YES ); + Idx_H2I = AddField( "H2I", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, NORMALIZE_YES, INTERP_FRAC_YES ); + Idx_H2II = AddField( "H2II", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, NORMALIZE_YES, INTERP_FRAC_YES ); } if ( GRACKLE_PRIMORDIAL >= GRACKLE_PRI_CHE_NSPE12 ) { - Idx_DI = AddField( "DI", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_YES, INTERP_FRAC_YES ); - Idx_DII = AddField( "DII", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_YES, INTERP_FRAC_YES ); - Idx_HDI = AddField( "HDI", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_YES, INTERP_FRAC_YES ); + Idx_DI = AddField( "DI", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, NORMALIZE_YES, INTERP_FRAC_YES ); + Idx_DII = AddField( "DII", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, NORMALIZE_YES, INTERP_FRAC_YES ); + Idx_HDI = AddField( "HDI", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, NORMALIZE_YES, INTERP_FRAC_YES ); } // normalize the metallicity field only when adopting the non-equilibrium chemistry // --> may need a machanism to allow users to overwrite this default setup if ( GRACKLE_METAL ) - Idx_Metal = AddField( "Metal", FIXUP_FLUX_YES, FIXUP_REST_YES, (GRACKLE_PRIMORDIAL==GRACKLE_PRI_CHE_CLOUDY)?NORMALIZE_NO:NORMALIZE_YES, - INTERP_FRAC_YES ); + Idx_Metal = AddField( "Metal", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, + (GRACKLE_PRIMORDIAL==GRACKLE_PRI_CHE_CLOUDY)?NORMALIZE_NO:NORMALIZE_YES, INTERP_FRAC_YES ); # endif // #ifdef SUPPORT_GRACKLE @@ -137,12 +138,12 @@ void Init_Field() // corresponding symbolic constants (e.g., DUAL/CRAY) defined in Macro.h // --> as we still rely on these constants (e.g., DENS, DUAL) in the fluid solvers # ifdef COSMIC_RAY - Idx_CRay = AddField( "CRay", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_NO, INTERP_FRAC_NO ); + Idx_CRay = AddField( "CRay", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, NORMALIZE_NO, INTERP_FRAC_NO ); if ( Idx_CRay != CRAY ) Aux_Error( ERROR_INFO, "inconsistent Idx_CRay (%d != %d) !!\n", Idx_CRay, CRAY ); # endif # ifdef DUAL_ENERGY - Idx_Dual = AddField( "Dual", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_NO, INTERP_FRAC_NO ); + Idx_Dual = AddField( "Dual", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, NORMALIZE_NO, INTERP_FRAC_NO ); if ( Idx_Dual != DUAL ) Aux_Error( ERROR_INFO, "inconsistent Idx_Dual (%d != %d) !!\n", Idx_Dual, DUAL ); # endif @@ -195,28 +196,31 @@ void Init_Field() // --> The fix-up flux operations will be applied to these scalars // (4) add the new field to FixUpVar_Restrict if FixUp_Restrict==true // --> The fix-up restrict operations will be applied to these scalars -// (5) add the new field to PassiveNorm_VarIdx[] if Norm==true +// (5) add the new field to PassiveFloorMask if Floor==true +// --> The floor operations will be applied to these passive scalars +// (6) add the new field to PassiveNorm_VarIdx[] if Norm==true // --> The sum of passive scalars in this list will be normalized to the gas density // --> sum(passive_scalar_density) == gas_density -// (6) add the new field to PassiveIntFrac_VarIdx[] if IntFrac==true -// --> These passive scalars will be converted to fracion form during interpolation +// (7) add the new field to PassiveIntFrac_VarIdx[] if IntFrac==true +// --> These passive scalars will be converted to fraction form during interpolation // 2. One must invoke AddField() exactly NCOMP_TOTAL times to set the labels of all fields // 3. Invoked by Init_Field() and various test problem initializers // // Parameter : InputLabel : Label (i.e., name) of the new field // FixUp_Flux : whether or not to apply fix-up flux operations // FixUp_Restrict : whether or not to apply fix-up restrict operations +// Floor : whether or not to floor the field to TINY_NUMBER (passive scalars only) // Norm : whether or not to normalize the new field // IntFrac : whether or not to convert the new field to fraction form during interpolation // // Return : (1) FieldLabel[] // (2) Index of the newly added field -// (3) FixUpVar_Flux and FixUpVar_Restrict +// (3) FixUpVar_Flux, FixUpVar_Restrict, and PassiveFloorMask // (4) PassiveNorm_NVar & PassiveNorm_VarIdx[] // (5) PassiveIntFrac_NVar & PassiveIntFrac_VarIdx[] //------------------------------------------------------------------------------------------------------- FieldIdx_t AddField( const char *InputLabel, const FixUpFlux_t FixUp_Flux, const FixUpRestrict_t FixUp_Restrict, - const NormPassive_t Norm, const IntFracPassive_t IntFrac ) + const FloorPassive_t Floor, const NormPassive_t Norm, const IntFracPassive_t IntFrac ) { const FieldIdx_t FieldIdx = NDefinedField ++; @@ -235,6 +239,13 @@ FieldIdx_t AddField( const char *InputLabel, const FixUpFlux_t FixUp_Flux, const if ( strcmp( FieldLabel[v], InputLabel ) == 0 ) Aux_Error( ERROR_INFO, "duplicate field label \"%s\" !!\n", InputLabel ); + if ( FieldIdx < NCOMP_FLUID && Floor != FLOOR_NULL ) + Aux_Error( ERROR_INFO, "Field index (%d) with passive floor option is not a passive scalar !!\n" + " --> Set FLOOR_NULL when adding a built-in fluid field by AddField()\n", FieldIdx ); + if ( FieldIdx >= NCOMP_FLUID && Floor == FLOOR_NULL ) + Aux_Error( ERROR_INFO, "Field index (%d) with Floor == FLOOR_NULL is a passive scalar !!\n" + " --> Set FLOOR_NO or FLOOR_YES when adding a passive scalar by AddField()\n", FieldIdx ); + // set field label strcpy( FieldLabel[FieldIdx], InputLabel ); @@ -245,6 +256,12 @@ FieldIdx_t AddField( const char *InputLabel, const FixUpFlux_t FixUp_Flux, const if ( FixUp_Restrict ) FixUpVar_Restrict |= (1L< note that PassiveFloorMask is written 1 for fields with FLOOR_NULL (not passive scalars) +// but this value should never be used as non-passive scalars should not be affected by it + if ( Floor ) PassiveFloorMask |= (1L< note that PassiveNorm_VarIdx[] starts from 0 instead of NCOMP_FLUID // --> currently we set PassiveNorm_VarIdx[] no mater OPT__NORMALIZE_PASSIVE is on or off @@ -263,6 +280,11 @@ FieldIdx_t AddField( const char *InputLabel, const FixUpFlux_t FixUp_Flux, const " --> Set NORMALIZE_NO when adding a built-in fluid field by AddField()\n", FieldIdx, NCOMP_FLUID ); + if ( Floor != FLOOR_YES ) + Aux_Error( ERROR_INFO, "Field index (%d) with normalization option is not floored !!\n" + " --> Set FLOOR_YES when adding a normalized passive scalar by AddField()\n", + FieldIdx ); + PassiveNorm_VarIdx[NormIdx] = FieldIdx - NCOMP_FLUID; } # endif @@ -348,6 +370,6 @@ void Init_Field_User_Template() { // example -// Idx_NewField = AddField( "NewFieldLabel", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_YES, INTERP_FRAC_YES ); +// Idx_NewField = AddField( "NewFieldLabel", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, NORMALIZE_YES, INTERP_FRAC_YES ); } // FUNCTION : Init_Field_User_Template diff --git a/src/Interpolation/Interpolate.cpp b/src/Interpolation/Interpolate.cpp index aa988f35bd..adaed12b54 100644 --- a/src/Interpolation/Interpolate.cpp +++ b/src/Interpolation/Interpolate.cpp @@ -277,7 +277,7 @@ void Interpolate_Iterate( real CData[], const int CSize[3], const int CStart[3], } # endif - Hydro_Con2Pri( Cons, Temp, MIN_PRES, + Hydro_Con2Pri( Cons, Temp, MIN_PRES, PassiveFloorMask, OPT__INT_FRAC_PASSIVE_LR, PassiveIntFrac_NVar, PassiveIntFrac_VarIdx, JeansMinPres_No, NULL_REAL, EoS_DensEint2Pres_CPUPtr, EoS_DensPres2Eint_CPUPtr, @@ -348,7 +348,7 @@ void Interpolate_Iterate( real CData[], const int CSize[3], const int CStart[3], if ( !FData_is_Prim ) Hydro_DualEnergyFix( Temp[DENS], Temp[MOMX], Temp[MOMY], Temp[MOMZ], Temp[ENGY], Temp[DUAL], dummy, EoS_AuxArray_Flt[1], EoS_AuxArray_Flt[2], - CheckMinPres_No, NULL_REAL, UseDual2FixEngy, Emag ); + CheckMinPres_No, NULL_REAL, PassiveFloorMask, UseDual2FixEngy, Emag ); # endif @@ -357,7 +357,7 @@ void Interpolate_Iterate( real CData[], const int CSize[3], const int CStart[3], = Hydro_IsUnphysical( (FData_is_Prim)?UNPHY_MODE_PRIM:UNPHY_MODE_CONS, Temp, Emag, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, - ERROR_INFO, UNPHY_SILENCE ); + PassiveFloorMask, ERROR_INFO, UNPHY_SILENCE ); // 5-3. additional check @@ -431,7 +431,7 @@ void Interpolate_Iterate( real CData[], const int CSize[3], const int CStart[3], else { const real CheckMinPres_No = false; const real Pres = Hydro_Con2Pres( Temp[DENS], Temp[MOMX], Temp[MOMY], Temp[MOMZ], Temp[ENGY], Temp+NCOMP_FLUID, - CheckMinPres_No, NULL_REAL, Emag, + CheckMinPres_No, NULL_REAL, PassiveFloorMask, Emag, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, &Eint ); Aux_Message( stderr, "Eint=%14.7e, Pres=%14.7e\n", Eint, Pres ); @@ -460,7 +460,7 @@ void Interpolate_Iterate( real CData[], const int CSize[3], const int CStart[3], if ( Hydro_IsUnphysical( UNPHY_MODE_CONS, Cons, Emag, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, - ERROR_INFO, UNPHY_VERBOSE ) ) + PassiveFloorMask, ERROR_INFO, UNPHY_VERBOSE ) ) Aux_Error( ERROR_INFO, "unphysical interpolated energy in %s() !!\n", __FUNCTION__ ); # endif } diff --git a/src/LoadBalance/LB_Refine_AllocateNewPatch.cpp b/src/LoadBalance/LB_Refine_AllocateNewPatch.cpp index 76a808b6d7..a5ecfcb3c3 100644 --- a/src/LoadBalance/LB_Refine_AllocateNewPatch.cpp +++ b/src/LoadBalance/LB_Refine_AllocateNewPatch.cpp @@ -1109,14 +1109,14 @@ int AllocateSonPatch( const int FaLv, const int *Cr, const int PScale, const int Hydro_DualEnergyFix( FData_Flu[DENS][k][j][i], FData_Flu[MOMX][k][j][i], FData_Flu[MOMY][k][j][i], FData_Flu[MOMZ][k][j][i], FData_Flu[ENGY][k][j][i], FData_Flu[DUAL][k][j][i], dummy, EoS_AuxArray_Flt[1], EoS_AuxArray_Flt[2], CheckMinPres_Yes, MIN_PRES, - UseDual2FixEngy, Emag ); + PassiveFloorMask, UseDual2FixEngy, Emag ); # else // #ifdef DUAL_ENERGY // apply internal energy floor FData_Flu[ENGY][k][j][i] = Hydro_CheckMinEintInEngy( FData_Flu[DENS][k][j][i], FData_Flu[MOMX][k][j][i], FData_Flu[MOMY][k][j][i], - FData_Flu[MOMZ][k][j][i], FData_Flu[ENGY][k][j][i], MIN_EINT, Emag ); + FData_Flu[MOMZ][k][j][i], FData_Flu[ENGY][k][j][i], MIN_EINT, PassiveFloorMask, Emag ); # endif // #ifdef DUAL_ENERGY ... else ... # endif // #if ( MODEL == HYDRO ) diff --git a/src/Main/InterpolateGhostZone.cpp b/src/Main/InterpolateGhostZone.cpp index 7506672dac..5e3f122419 100644 --- a/src/Main/InterpolateGhostZone.cpp +++ b/src/Main/InterpolateGhostZone.cpp @@ -438,7 +438,7 @@ void InterpolateGhostZone( const int lv, const int PID, real IntData_CC[], real # endif CData_CC_Ptr[Idx] = Hydro_Con2Pres( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], FluidForEoS+NCOMP_FLUID, - (MinPres>=(real)0.0), MinPres, Emag, + (MinPres>=(real)0.0), MinPres, PassiveFloorMask, Emag, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, NULL ); @@ -456,7 +456,7 @@ void InterpolateGhostZone( const int lv, const int PID, real IntData_CC[], real FluWeighting *CData_CC_Ptr[Idx] + FluWeighting_IntT*Hydro_Con2Pres( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], FluidForEoS+NCOMP_FLUID, - (MinPres>=(real)0.0), MinPres, Emag, + (MinPres>=(real)0.0), MinPres, PassiveFloorMask, Emag, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, NULL ); @@ -485,7 +485,7 @@ void InterpolateGhostZone( const int lv, const int PID, real IntData_CC[], real # endif CData_CC_Ptr[Idx] = Hydro_Con2Temp( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], FluidForEoS+NCOMP_FLUID, - (MinTemp>=(real)0.0), MinTemp, Emag, + (MinTemp>=(real)0.0), MinTemp, PassiveFloorMask, Emag, EoS_DensEint2Temp_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); @@ -502,7 +502,7 @@ void InterpolateGhostZone( const int lv, const int PID, real IntData_CC[], real FluWeighting *CData_CC_Ptr[Idx] + FluWeighting_IntT*Hydro_Con2Temp( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], FluidForEoS+NCOMP_FLUID, - (MinTemp>=(real)0.0), MinTemp, Emag, + (MinTemp>=(real)0.0), MinTemp, PassiveFloorMask, Emag, EoS_DensEint2Temp_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); } @@ -530,7 +530,7 @@ void InterpolateGhostZone( const int lv, const int PID, real IntData_CC[], real # endif CData_CC_Ptr[Idx] = Hydro_Con2Entr( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], FluidForEoS+NCOMP_FLUID, - (MinEntr>=(real)0.0), MinEntr, Emag, + (MinEntr>=(real)0.0), MinEntr, PassiveFloorMask, Emag, EoS_DensEint2Entr_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); @@ -547,7 +547,7 @@ void InterpolateGhostZone( const int lv, const int PID, real IntData_CC[], real FluWeighting *CData_CC_Ptr[Idx] + FluWeighting_IntT*Hydro_Con2Entr( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], FluidForEoS+NCOMP_FLUID, - (MinEntr>=(real)0.0), MinEntr, Emag, + (MinEntr>=(real)0.0), MinEntr, PassiveFloorMask, Emag, EoS_DensEint2Entr_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); } @@ -577,7 +577,7 @@ void InterpolateGhostZone( const int lv, const int PID, real IntData_CC[], real const bool CheckMinEint_No = false; // floor value is not supported for now CData_CC_Ptr[Idx] = Hydro_Con2Eint( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], - CheckMinEint_No, NULL_REAL, Emag, + CheckMinEint_No, NULL_REAL, PassiveFloorMask, Emag, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); @@ -594,7 +594,7 @@ void InterpolateGhostZone( const int lv, const int PID, real IntData_CC[], real FluWeighting *CData_CC_Ptr[Idx] + FluWeighting_IntT*Hydro_Con2Eint( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], - CheckMinEint_No, NULL_REAL, Emag, + CheckMinEint_No, NULL_REAL, PassiveFloorMask, Emag, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); } @@ -891,7 +891,7 @@ void InterpolateGhostZone( const int lv, const int PID, real IntData_CC[], real # endif CData_CC_Ptr[Idx] = Hydro_Con2Pres( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], FluidForEoS+NCOMP_FLUID, - (MinPres>=(real)0.0), MinPres, Emag, + (MinPres>=(real)0.0), MinPres, PassiveFloorMask, Emag, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, NULL ); @@ -909,7 +909,7 @@ void InterpolateGhostZone( const int lv, const int PID, real IntData_CC[], real FluWeighting *CData_CC_Ptr[Idx] + FluWeighting_IntT*Hydro_Con2Pres( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], FluidForEoS+NCOMP_FLUID, - (MinPres>=(real)0.0), MinPres, Emag, + (MinPres>=(real)0.0), MinPres, PassiveFloorMask, Emag, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, NULL ); @@ -938,7 +938,7 @@ void InterpolateGhostZone( const int lv, const int PID, real IntData_CC[], real # endif CData_CC_Ptr[Idx] = Hydro_Con2Temp( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], FluidForEoS+NCOMP_FLUID, - (MinTemp>=(real)0.0), MinTemp, Emag, + (MinTemp>=(real)0.0), MinTemp, PassiveFloorMask, Emag, EoS_DensEint2Temp_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); @@ -955,7 +955,7 @@ void InterpolateGhostZone( const int lv, const int PID, real IntData_CC[], real FluWeighting *CData_CC_Ptr[Idx] + FluWeighting_IntT*Hydro_Con2Temp( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], FluidForEoS+NCOMP_FLUID, - (MinTemp>=(real)0.0), MinTemp, Emag, + (MinTemp>=(real)0.0), MinTemp, PassiveFloorMask, Emag, EoS_DensEint2Temp_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); } @@ -983,7 +983,7 @@ void InterpolateGhostZone( const int lv, const int PID, real IntData_CC[], real # endif CData_CC_Ptr[Idx] = Hydro_Con2Entr( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], FluidForEoS+NCOMP_FLUID, - (MinEntr>=(real)0.0), MinEntr, Emag, + (MinEntr>=(real)0.0), MinEntr, PassiveFloorMask, Emag, EoS_DensEint2Entr_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); @@ -1000,7 +1000,7 @@ void InterpolateGhostZone( const int lv, const int PID, real IntData_CC[], real FluWeighting *CData_CC_Ptr[Idx] + FluWeighting_IntT*Hydro_Con2Entr( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], FluidForEoS+NCOMP_FLUID, - (MinEntr>=(real)0.0), MinEntr, Emag, + (MinEntr>=(real)0.0), MinEntr, PassiveFloorMask, Emag, EoS_DensEint2Entr_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); } @@ -1030,7 +1030,7 @@ void InterpolateGhostZone( const int lv, const int PID, real IntData_CC[], real const bool CheckMinEint_No = false; // floor value is not supported for now CData_CC_Ptr[Idx] = Hydro_Con2Eint( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], - CheckMinEint_No, NULL_REAL, Emag, + CheckMinEint_No, NULL_REAL, PassiveFloorMask, Emag, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); @@ -1047,7 +1047,7 @@ void InterpolateGhostZone( const int lv, const int PID, real IntData_CC[], real FluWeighting *CData_CC_Ptr[Idx] + FluWeighting_IntT*Hydro_Con2Eint( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], - CheckMinEint_No, NULL_REAL, Emag, + CheckMinEint_No, NULL_REAL, PassiveFloorMask, Emag, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); } @@ -1978,7 +1978,7 @@ void InterpolateGhostZone( const int lv, const int PID, real IntData_CC[], real // the runtime parameter DUAL_ENERGY_SWITCH here Hydro_DualEnergyFix( FData_Dens[t], FData_MomX[t], FData_MomY[t], FData_MomZ[t], FData_Engy[t], FData_Dual[t], dummy, EoS_AuxArray_Flt[1], EoS_AuxArray_Flt[2], (MinPres>=(real)0.0), MinPres, - UseDual2FixEngy, Emag ); + PassiveFloorMask, UseDual2FixEngy, Emag ); } } // if ( DE_Consistency && ( TVarCC & _TOTAL ) == _TOTAL && TVarFC == _MAG ) # endif // if ( MODEL == HYDRO && defined DUAL_ENERGY ) diff --git a/src/Main/InvokeSolver.cpp b/src/Main/InvokeSolver.cpp index 32a9d7b168..fc29ddd131 100644 --- a/src/Main/InvokeSolver.cpp +++ b/src/Main/InvokeSolver.cpp @@ -579,7 +579,7 @@ void Solver( const Solver_t TSolver, const int lv, const double TimeNew, const d OPT__LR_LIMITER, MINMOD_COEFF, MINMOD_MAX_ITER, ELBDM_ETA, ELBDM_TAYLOR3_COEFF, ELBDM_TAYLOR3_AUTO, TimeOld, (OPT__SELF_GRAVITY || OPT__EXT_POT), OPT__EXT_ACC, MicroPhy, - MIN_DENS, MIN_PRES, MIN_EINT, DUAL_ENERGY_SWITCH, + MIN_DENS, MIN_PRES, MIN_EINT, DUAL_ENERGY_SWITCH, PassiveFloorMask, OPT__NORMALIZE_PASSIVE, PassiveNorm_NVar, OPT__INT_FRAC_PASSIVE_LR, PassiveIntFrac_NVar, JEANS_MIN_PRES, JeansMinPres_Coeff, @@ -596,7 +596,7 @@ void Solver( const Solver_t TSolver, const int lv, const double TimeNew, const d OPT__LR_LIMITER, MINMOD_COEFF, MINMOD_MAX_ITER, ELBDM_ETA, ELBDM_TAYLOR3_COEFF, ELBDM_TAYLOR3_AUTO, TimeOld, (OPT__SELF_GRAVITY || OPT__EXT_POT), OPT__EXT_ACC, MicroPhy, - MIN_DENS, MIN_PRES, MIN_EINT, DUAL_ENERGY_SWITCH, + MIN_DENS, MIN_PRES, MIN_EINT, DUAL_ENERGY_SWITCH, PassiveFloorMask, OPT__NORMALIZE_PASSIVE, PassiveNorm_NVar, PassiveNorm_VarIdx, OPT__INT_FRAC_PASSIVE_LR, PassiveIntFrac_NVar, PassiveIntFrac_VarIdx, JEANS_MIN_PRES, JeansMinPres_Coeff, UseWaveFlag ); @@ -705,14 +705,14 @@ void Solver( const Solver_t TSolver, const int lv, const double TimeNew, const d CUAPI_Asyn_dtSolver( TSolver, h_dt_Array_T[ArrayID], h_Flu_Array_T[ArrayID], h_Mag_Array_T[ArrayID], NULL, NULL, NPG, dh, (Step==0)?DT__FLUID_INIT:DT__FLUID, - MicroPhy, MIN_PRES, NULL_BOOL, + MicroPhy, MIN_PRES, PassiveFloorMask, NULL_BOOL, EXT_POT_NONE, EXT_ACC_NONE, NULL_REAL, GPU_NSTREAM ); # else CPU_dtSolver ( TSolver, h_dt_Array_T[ArrayID], h_Flu_Array_T[ArrayID], h_Mag_Array_T[ArrayID], NULL, NULL, NPG, dh, (Step==0)?DT__FLUID_INIT:DT__FLUID, - MicroPhy, MIN_PRES, NULL_BOOL, + MicroPhy, MIN_PRES, PassiveFloorMask, NULL_BOOL, EXT_POT_NONE, EXT_ACC_NONE, NULL_REAL ); # endif break; @@ -723,14 +723,14 @@ void Solver( const Solver_t TSolver, const int lv, const double TimeNew, const d CUAPI_Asyn_dtSolver( TSolver, h_dt_Array_T[ArrayID], NULL, NULL, h_Pot_Array_T[ArrayID], h_Corner_Array_PGT[ArrayID], NPG, dh, DT__GRAVITY, - MicroPhy, NULL_REAL, OPT__GRA_P5_GRADIENT, + MicroPhy, NULL_REAL, PassiveFloorMask, OPT__GRA_P5_GRADIENT, (OPT__SELF_GRAVITY || OPT__EXT_POT), OPT__EXT_ACC, TimeNew, GPU_NSTREAM ); # else CPU_dtSolver ( TSolver, h_dt_Array_T[ArrayID], NULL, NULL, h_Pot_Array_T[ArrayID], h_Corner_Array_PGT[ArrayID], NPG, dh, DT__GRAVITY, - MicroPhy, NULL_REAL, OPT__GRA_P5_GRADIENT, + MicroPhy, NULL_REAL, PassiveFloorMask, OPT__GRA_P5_GRADIENT, (OPT__SELF_GRAVITY || OPT__EXT_POT), OPT__EXT_ACC, TimeNew ); # endif break; @@ -753,14 +753,14 @@ void Solver( const Solver_t TSolver, const int lv, const double TimeNew, const d h_Flu_Array_S_Out[ArrayID], h_Mag_Array_S_In [ArrayID], h_Corner_Array_S [ArrayID], - SrcTerms, NPG, dt, dh, TimeNew, TimeOld, MIN_DENS, MIN_PRES, MIN_EINT, + SrcTerms, NPG, dt, dh, TimeNew, TimeOld, MIN_DENS, MIN_PRES, MIN_EINT, PassiveFloorMask, GPU_NSTREAM ); # else CPU_SrcSolver ( h_Flu_Array_S_In [ArrayID], h_Flu_Array_S_Out[ArrayID], h_Mag_Array_S_In [ArrayID], h_Corner_Array_S [ArrayID], - SrcTerms, NPG, dt, dh, TimeNew, TimeOld, MIN_DENS, MIN_PRES, MIN_EINT ); + SrcTerms, NPG, dt, dh, TimeNew, TimeOld, MIN_DENS, MIN_PRES, MIN_EINT, PassiveFloorMask ); # endif break; diff --git a/src/Main/Main.cpp b/src/Main/Main.cpp index a495c688a5..70358820de 100644 --- a/src/Main/Main.cpp +++ b/src/Main/Main.cpp @@ -38,7 +38,7 @@ double *FlagTable_User [NLEVEL-1]; double *DumpTable = NULL; int DumpTable_NDump; int *UM_IC_RefineRegion = NULL; -long FixUpVar_Flux, FixUpVar_Restrict; +long FixUpVar_Flux, FixUpVar_Restrict, PassiveFloorMask; int PassiveNorm_NVar, PassiveNorm_VarIdx[NCOMP_PASSIVE]; int PassiveIntFrac_NVar, PassiveIntFrac_VarIdx[NCOMP_PASSIVE]; int StrLen_Flt; diff --git a/src/Main/Prepare_PatchData.cpp b/src/Main/Prepare_PatchData.cpp index 174c1dffbf..1f6fe8d121 100644 --- a/src/Main/Prepare_PatchData.cpp +++ b/src/Main/Prepare_PatchData.cpp @@ -761,7 +761,7 @@ void Prepare_PatchData( const int lv, const double PrepTime, real *OutputCC, rea # endif Data1PG_CC_Ptr[Idx1] = Hydro_Con2Pres( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], FluidForEoS+NCOMP_FLUID, - (MinPres>=(real)0.0), MinPres, Emag, + (MinPres>=(real)0.0), MinPres, PassiveFloorMask, Emag, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, NULL ); @@ -779,7 +779,7 @@ void Prepare_PatchData( const int lv, const double PrepTime, real *OutputCC, rea FluWeighting *Data1PG_CC_Ptr[Idx1] + FluWeighting_IntT*Hydro_Con2Pres( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], FluidForEoS+NCOMP_FLUID, - (MinPres>=(real)0.0), MinPres, Emag, + (MinPres>=(real)0.0), MinPres, PassiveFloorMask, Emag, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, NULL ); @@ -808,7 +808,7 @@ void Prepare_PatchData( const int lv, const double PrepTime, real *OutputCC, rea # endif Data1PG_CC_Ptr[Idx1] = Hydro_Con2Temp( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], FluidForEoS+NCOMP_FLUID, - (MinTemp>=(real)0.0), MinTemp, Emag, + (MinTemp>=(real)0.0), MinTemp, PassiveFloorMask, Emag, EoS_DensEint2Temp_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); @@ -825,7 +825,7 @@ void Prepare_PatchData( const int lv, const double PrepTime, real *OutputCC, rea FluWeighting *Data1PG_CC_Ptr[Idx1] + FluWeighting_IntT*Hydro_Con2Temp( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], FluidForEoS+NCOMP_FLUID, - (MinTemp>=(real)0.0), MinTemp, Emag, + (MinTemp>=(real)0.0), MinTemp, PassiveFloorMask, Emag, EoS_DensEint2Temp_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); } @@ -853,7 +853,7 @@ void Prepare_PatchData( const int lv, const double PrepTime, real *OutputCC, rea # endif Data1PG_CC_Ptr[Idx1] = Hydro_Con2Entr( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], FluidForEoS+NCOMP_FLUID, - (MinEntr>=(real)0.0), MinEntr, Emag, + (MinEntr>=(real)0.0), MinEntr, PassiveFloorMask, Emag, EoS_DensEint2Entr_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); @@ -870,7 +870,7 @@ void Prepare_PatchData( const int lv, const double PrepTime, real *OutputCC, rea FluWeighting *Data1PG_CC_Ptr[Idx1] + FluWeighting_IntT*Hydro_Con2Entr( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], FluidForEoS+NCOMP_FLUID, - (MinEntr>=(real)0.0), MinEntr, Emag, + (MinEntr>=(real)0.0), MinEntr, PassiveFloorMask, Emag, EoS_DensEint2Entr_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); } @@ -900,7 +900,7 @@ void Prepare_PatchData( const int lv, const double PrepTime, real *OutputCC, rea const bool CheckMinEint_No = false; // floor value is not supported for now Data1PG_CC_Ptr[Idx1] = Hydro_Con2Eint( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], - CheckMinEint_No, NULL_REAL, Emag, + CheckMinEint_No, NULL_REAL, PassiveFloorMask, Emag, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); @@ -917,7 +917,7 @@ void Prepare_PatchData( const int lv, const double PrepTime, real *OutputCC, rea FluWeighting *Data1PG_CC_Ptr[Idx1] + FluWeighting_IntT*Hydro_Con2Eint( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], - CheckMinEint_No, NULL_REAL, Emag, + CheckMinEint_No, NULL_REAL, PassiveFloorMask, Emag, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); } @@ -1213,7 +1213,7 @@ void Prepare_PatchData( const int lv, const double PrepTime, real *OutputCC, rea # endif Data1PG_CC_Ptr[Idx1] = Hydro_Con2Pres( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], FluidForEoS+NCOMP_FLUID, - (MinPres>=(real)0.0), MinPres, Emag, + (MinPres>=(real)0.0), MinPres, PassiveFloorMask, Emag, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, NULL ); @@ -1231,7 +1231,7 @@ void Prepare_PatchData( const int lv, const double PrepTime, real *OutputCC, rea FluWeighting *Data1PG_CC_Ptr[Idx1] + FluWeighting_IntT*Hydro_Con2Pres( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], FluidForEoS+NCOMP_FLUID, - (MinPres>=(real)0.0), MinPres, Emag, + (MinPres>=(real)0.0), MinPres, PassiveFloorMask, Emag, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, NULL ); @@ -1260,7 +1260,7 @@ void Prepare_PatchData( const int lv, const double PrepTime, real *OutputCC, rea # endif Data1PG_CC_Ptr[Idx1] = Hydro_Con2Temp( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], FluidForEoS+NCOMP_FLUID, - (MinTemp>=(real)0.0), MinTemp, Emag, + (MinTemp>=(real)0.0), MinTemp, PassiveFloorMask, Emag, EoS_DensEint2Temp_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); @@ -1277,7 +1277,7 @@ void Prepare_PatchData( const int lv, const double PrepTime, real *OutputCC, rea FluWeighting *Data1PG_CC_Ptr[Idx1] + FluWeighting_IntT*Hydro_Con2Temp( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], FluidForEoS+NCOMP_FLUID, - (MinTemp>=(real)0.0), MinTemp, Emag, + (MinTemp>=(real)0.0), MinTemp, PassiveFloorMask, Emag, EoS_DensEint2Temp_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); } @@ -1305,7 +1305,7 @@ void Prepare_PatchData( const int lv, const double PrepTime, real *OutputCC, rea # endif Data1PG_CC_Ptr[Idx1] = Hydro_Con2Entr( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], FluidForEoS+NCOMP_FLUID, - (MinEntr>=(real)0.0), MinEntr, Emag, + (MinEntr>=(real)0.0), MinEntr, PassiveFloorMask, Emag, EoS_DensEint2Entr_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); @@ -1322,7 +1322,7 @@ void Prepare_PatchData( const int lv, const double PrepTime, real *OutputCC, rea FluWeighting *Data1PG_CC_Ptr[Idx1] + FluWeighting_IntT*Hydro_Con2Entr( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], FluidForEoS+NCOMP_FLUID, - (MinEntr>=(real)0.0), MinEntr, Emag, + (MinEntr>=(real)0.0), MinEntr, PassiveFloorMask, Emag, EoS_DensEint2Entr_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); } @@ -1352,7 +1352,7 @@ void Prepare_PatchData( const int lv, const double PrepTime, real *OutputCC, rea const bool CheckMinEint_No = false; // floor value is not supported for now Data1PG_CC_Ptr[Idx1] = Hydro_Con2Eint( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], - CheckMinEint_No, NULL_REAL, Emag, + CheckMinEint_No, NULL_REAL, PassiveFloorMask, Emag, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); @@ -1369,7 +1369,7 @@ void Prepare_PatchData( const int lv, const double PrepTime, real *OutputCC, rea FluWeighting *Data1PG_CC_Ptr[Idx1] + FluWeighting_IntT*Hydro_Con2Eint( FluidForEoS[DENS], FluidForEoS[MOMX], FluidForEoS[MOMY], FluidForEoS[MOMZ], FluidForEoS[ENGY], - CheckMinEint_No, NULL_REAL, Emag, + CheckMinEint_No, NULL_REAL, PassiveFloorMask, Emag, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); } diff --git a/src/Miscellaneous/CPU_dtSolver.cpp b/src/Miscellaneous/CPU_dtSolver.cpp index 1f3ecf9d0e..a0a8fc95c1 100644 --- a/src/Miscellaneous/CPU_dtSolver.cpp +++ b/src/Miscellaneous/CPU_dtSolver.cpp @@ -7,7 +7,7 @@ void CPU_dtSolver_HydroCFL( real g_dt_Array[], const real g_Flu_Array[][FLU_NIN_T][ CUBE(PS1) ], const real g_Mag_Array[][NCOMP_MAG][ PS1P1*SQR(PS1) ], const int NPG, const real dh, const real Safety, const real MinPres, - const EoS_t EoS, const MicroPhy_t MicroPhy ); + const long PassiveFloor, const EoS_t EoS, const MicroPhy_t MicroPhy ); #ifdef GRAVITY void CPU_dtSolver_HydroGravity( real g_dt_Array[], const real g_Pot_Array[][ CUBE(GRA_NXT) ], @@ -47,6 +47,7 @@ void CPU_dtSolver_HydroGravity( real g_dt_Array[], // Safety : dt safety factor // MicroPhy : Microphysics object // MinPres : Minimum allowed pressure +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored // P5_Gradient : Use 5-points stencil to evaluate the potential gradient // UsePot : Add self-gravity and/or external potential // ExtAcc : Add external acceleration @@ -57,7 +58,7 @@ void CPU_dtSolver_HydroGravity( real g_dt_Array[], void CPU_dtSolver( const Solver_t TSolver, real dt_Array[], const real Flu_Array[][FLU_NIN_T][ CUBE(PS1) ], const real Mag_Array[][NCOMP_MAG][ PS1P1*SQR(PS1) ], const real Pot_Array[][ CUBE(GRA_NXT) ], const double Corner_Array[][3], const int NPatchGroup, const real dh, const real Safety, - const MicroPhy_t MicroPhy, const real MinPres, const bool P5_Gradient, + const MicroPhy_t MicroPhy, const real MinPres, const long PassiveFloor, const bool P5_Gradient, const bool UsePot, const OptExtAcc_t ExtAcc, const double TargetTime ) { @@ -65,7 +66,7 @@ void CPU_dtSolver( const Solver_t TSolver, real dt_Array[], const real Flu_Array { # if ( MODEL == HYDRO ) case DT_FLU_SOLVER: - CPU_dtSolver_HydroCFL( dt_Array, Flu_Array, Mag_Array, NPatchGroup, dh, Safety, MinPres, EoS, MicroPhy ); + CPU_dtSolver_HydroCFL( dt_Array, Flu_Array, Mag_Array, NPatchGroup, dh, Safety, MinPres, PassiveFloor, EoS, MicroPhy ); break; # ifdef GRAVITY diff --git a/src/Model_ELBDM/ELBDM_Init_ByFunction_AssignData.cpp b/src/Model_ELBDM/ELBDM_Init_ByFunction_AssignData.cpp index 0f4c2e50ea..ed0c473d2e 100644 --- a/src/Model_ELBDM/ELBDM_Init_ByFunction_AssignData.cpp +++ b/src/Model_ELBDM/ELBDM_Init_ByFunction_AssignData.cpp @@ -172,7 +172,8 @@ void ELBDM_Init_ByFunction_AssignData( const int lv ) // floor and normalize passive scalars (actually passive scalars are NOT supported by ELBDM yet) /* # if ( NCOMP_PASSIVE > 0 ) - for (int v=NCOMP_FLUID; v normalize passive scalars so that the sum of their mass density // is equal to the gas mass density // NNorm : Number of passive scalars to be normalized @@ -174,6 +176,7 @@ void CUFLU_FluidSolver_CTU( const bool UsePot, const OptExtAcc_t ExtAcc, const ExtAcc_t ExtAcc_Func, const real MinDens, const real MinPres, const real MinEint, const real DualEnergySwitch, + const long PassiveFloor, const bool NormPassive, const int NNorm, const bool FracPassive, const int NFrac, const bool JeansMinPres, const real JeansMinPres_Coeff, @@ -203,6 +206,7 @@ void CPU_FluidSolver_CTU( const double c_ExtAcc_AuxArray[], const real MinDens, const real MinPres, const real MinEint, const real DualEnergySwitch, + const long PassiveFloor, const bool NormPassive, const int NNorm, const int c_NormIdx[], const bool FracPassive, const int NFrac, const int c_FracIdx[], const bool JeansMinPres, const real JeansMinPres_Coeff, @@ -260,7 +264,7 @@ void CPU_FluidSolver_CTU( // 1. evaluate the face-centered values at the half time-step Hydro_DataReconstruction( g_Flu_Array_In[P], g_Mag_Array_In[P], g_PriVar_1PG, g_FC_Var_1PG, g_Slope_PPM_1PG, NULL, Con2Pri_Yes, LR_Limiter, MinMod_Coeff, dt, dh, - MinDens, MinPres, MinEint, FracPassive, NFrac, c_FracIdx, + MinDens, MinPres, MinEint, PassiveFloor, FracPassive, NFrac, c_FracIdx, JeansMinPres, JeansMinPres_Coeff, &EoS ); @@ -268,7 +272,7 @@ void CPU_FluidSolver_CTU( Hydro_ComputeFlux( g_FC_Var_1PG, g_FC_Flux_1PG, N_HF_FLUX, 0, 0, CorrHalfVel_No, NULL, NULL, NULL_REAL, NULL_REAL, NULL_REAL, EXT_POT_NONE, EXT_ACC_NONE, NULL, NULL, - MinDens, MinPres, &EoS ); + MinDens, MinPres, PassiveFloor, &EoS ); // 3. evaluate electric field and update B field at the half time-step @@ -285,7 +289,7 @@ void CPU_FluidSolver_CTU( // 4. correct the face-centered variables by the transverse flux gradients Hydro_TGradientCorrection( g_FC_Var_1PG, g_FC_Flux_1PG, g_Mag_Array_In[P], g_FC_Mag_Half_1PG, g_EC_Ele_1PG, g_PriVar_1PG, - dt, dh, MinDens, MinEint ); + dt, dh, MinDens, MinEint, PassiveFloor ); // 5. evaluate the cell-centered primitive variables at the half time-step @@ -306,7 +310,7 @@ void CPU_FluidSolver_CTU( Hydro_ComputeFlux( g_FC_Var_1PG, g_FC_Flux_1PG, N_FL_FLUX, NSkip_N, NSkip_T, CorrHalfVel, g_Pot_Array_USG[P], g_Corner_Array[P], dt, dh, Time, UsePot, ExtAcc, ExtAcc_Func, c_ExtAcc_AuxArray, - MinDens, MinPres, &EoS ); + MinDens, MinPres, PassiveFloor, &EoS ); if ( StoreFlux ) Hydro_StoreIntFlux( g_FC_Flux_1PG, g_Flux_Array[P], N_FL_FLUX ); @@ -330,7 +334,7 @@ void CPU_FluidSolver_CTU( // --> CTU does not support reducing the min-mod coefficient Hydro_FullStepUpdate( g_Flu_Array_In[P], g_Flu_Array_Out[P], g_DE_Array_Out[P], g_Mag_Array_Out[P], g_FC_Flux_1PG, dt, dh, MinDens, MinEint, DualEnergySwitch, - NormPassive, NNorm, c_NormIdx, &EoS, NULL, NULL_INT, NULL_INT ); + PassiveFloor, NormPassive, NNorm, c_NormIdx, &EoS, NULL, NULL_INT, NULL_INT ); } // loop over all patch groups } // OpenMP parallel region @@ -358,6 +362,7 @@ void CPU_FluidSolver_CTU( // dt : Time interval to advance solution // dh : Cell size // MinDens/Eint : Density and internal energy floors +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored //------------------------------------------------------------------------------------------------------- GPU_DEVICE void Hydro_TGradientCorrection( real g_FC_Var [][NCOMP_TOTAL_PLUS_MAG][ CUBE(N_FC_VAR) ], @@ -366,7 +371,8 @@ void Hydro_TGradientCorrection( real g_FC_Var [][NCOMP_TOTAL_PLUS_MAG][ const real g_FC_B_Half[][ FLU_NXT_P1*SQR(FLU_NXT) ], const real g_EC_Ele [][ CUBE(N_EC_ELE) ], const real g_PriVar [][ CUBE(FLU_NXT) ], - const real dt, const real dh, const real MinDens, const real MinEint ) + const real dt, const real dh, const real MinDens, const real MinEint, + const long PassiveFloor ) { const int didx_flux[3] = { 1, N_HF_FLUX, SQR(N_HF_FLUX) }; @@ -567,9 +573,10 @@ void Hydro_TGradientCorrection( real g_FC_Var [][NCOMP_TOTAL_PLUS_MAG][ # endif fc_var[f][0] = FMAX( fc_var[f][0], MinDens ); fc_var[f][4] = Hydro_CheckMinEintInEngy( fc_var[f][0], fc_var[f][1], fc_var[f][2], fc_var[f][3], fc_var[f][4], - MinEint, Emag ); + MinEint, PassiveFloor, Emag ); # if ( NCOMP_PASSIVE > 0 ) for (int v=NCOMP_FLUID; v normalize passive scalars so that the sum of their mass density // is equal to the gas mass density // NNorm : Number of passive scalars to be normalized @@ -283,6 +284,7 @@ void CUFLU_FluidSolver_MHM( const bool UsePot, const OptExtAcc_t ExtAcc, const ExtAcc_t ExtAcc_Func, const real MinDens, const real MinPres, const real MinEint, const real DualEnergySwitch, + const long PassiveFloor, const bool NormPassive, const int NNorm, const bool FracPassive, const int NFrac, const bool JeansMinPres, const real JeansMinPres_Coeff, @@ -312,6 +314,7 @@ void CPU_FluidSolver_MHM( const double c_ExtAcc_AuxArray[], const real MinDens, const real MinPres, const real MinEint, const real DualEnergySwitch, + const long PassiveFloor, const bool NormPassive, const int NNorm, const int c_NormIdx[], const bool FracPassive, const int NFrac, const int c_FracIdx[], const bool JeansMinPres, const real JeansMinPres_Coeff, @@ -435,7 +438,7 @@ void CPU_FluidSolver_MHM( // 1-a-2. evaluate the half-step first-order fluxes by Riemann solver // hydrodynamic fluxes Hydro_RiemannPredict_Flux( g_Flu_Array_In[P], g_Flux_Half_1PG, g_Mag_Array_In[P], g_PriVar_1PG+MAG_OFFSET, - MinDens, MinPres, &EoS ); + MinDens, MinPres, PassiveFloor, &EoS ); // add cosmic-ray fluxes # ifdef CR_DIFFUSION @@ -457,7 +460,7 @@ void CPU_FluidSolver_MHM( // 1-a-4. evaluate the half-step solutions Hydro_RiemannPredict( g_Flu_Array_In[P], g_FC_Mag_Half_1PG, g_Flux_Half_1PG, g_PriVar_Half_1PG, - dt, dh, MinDens, MinPres, MinEint, FracPassive, NFrac, c_FracIdx, + dt, dh, MinDens, MinPres, MinEint, PassiveFloor, FracPassive, NFrac, c_FracIdx, JeansMinPres, JeansMinPres_Coeff, &EoS ); @@ -483,7 +486,7 @@ void CPU_FluidSolver_MHM( // --> note that g_PriVar_Half_1PG[] returned by Hydro_RiemannPredict() stores the primitive variables Hydro_DataReconstruction( NULL, g_FC_Mag_Half_1PG, g_PriVar_Half_1PG, g_FC_Var_1PG, g_Slope_PPM_1PG, NULL, Con2Pri_No, LR_Limiter, AdaptiveMinModCoeff, dt, dh, - MinDens, MinPres, MinEint, FracPassive, NFrac, c_FracIdx, + MinDens, MinPres, MinEint, PassiveFloor, FracPassive, NFrac, c_FracIdx, JeansMinPres, JeansMinPres_Coeff, &EoS ); @@ -511,7 +514,7 @@ void CPU_FluidSolver_MHM( // evaluate the face-centered values by data reconstruction Hydro_DataReconstruction( g_Flu_Array_In[P], g_Mag_Array_In[P], g_PriVar_Half_1PG, g_FC_Var_1PG, g_Slope_PPM_1PG, g_EC_Ele_1PG, Con2Pri_Yes, LR_Limiter, AdaptiveMinModCoeff, dt, dh, - MinDens, MinPres, MinEint, FracPassive, NFrac, c_FracIdx, + MinDens, MinPres, MinEint, PassiveFloor, FracPassive, NFrac, c_FracIdx, JeansMinPres, JeansMinPres_Coeff, &EoS ); # endif // #if ( FLU_SCHEME == MHM_RP ) ... else ... @@ -530,7 +533,7 @@ void CPU_FluidSolver_MHM( Hydro_ComputeFlux( g_FC_Var_1PG, g_FC_Flux_1PG, N_FL_FLUX, NSkip_N, NSkip_T, CorrHalfVel, g_Pot_Array_USG[P], g_Corner_Array[P], dt, dh, Time, UsePot, ExtAcc, ExtAcc_Func, c_ExtAcc_AuxArray, - MinDens, MinPres, &EoS ); + MinDens, MinPres, PassiveFloor, &EoS ); // add cosmic-ray fluxes # ifdef CR_DIFFUSION @@ -565,7 +568,8 @@ void CPU_FluidSolver_MHM( // 4. full-step evolution Hydro_FullStepUpdate( g_Flu_Array_In[P], g_Flu_Array_Out[P], g_DE_Array_Out[P], g_Mag_Array_Out[P], g_FC_Flux_1PG, dt, dh, MinDens, MinEint, DualEnergySwitch, - NormPassive, NNorm, c_NormIdx, &EoS, &s_FullStepFailure, Iteration, MinMod_MaxIter ); + PassiveFloor, NormPassive, NNorm, c_NormIdx, &EoS, &s_FullStepFailure, + Iteration, MinMod_MaxIter ); // add the cosmic-ray source term of adiabatic work # ifdef COSMIC_RAY @@ -605,6 +609,7 @@ void CPU_FluidSolver_MHM( // g_CC_B : Array storing the input cell-centered magnetic field (for MHD only) // --> Accessed with a stride FLU_NXT // MinDens/Pres : Density and pressure floors +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored // EoS : EoS object //------------------------------------------------------------------------------------------------------- GPU_DEVICE @@ -612,7 +617,7 @@ void Hydro_RiemannPredict_Flux( const real g_ConVar[][ CUBE(FLU_NXT) ], real g_Flux_Half[][NCOMP_TOTAL_PLUS_MAG][ CUBE(N_FC_FLUX) ], const real g_FC_B[][ SQR(FLU_NXT)*FLU_NXT_P1 ], const real g_CC_B[][ CUBE(FLU_NXT) ], - const real MinDens, const real MinPres, + const real MinDens, const real MinPres, const long PassiveFloor, const EoS_t *EoS ) { @@ -703,25 +708,25 @@ void Hydro_RiemannPredict_Flux( const real g_ConVar[][ CUBE(FLU_NXT) ], // invoke the Riemann solver # if ( RSOLVER == EXACT && !defined MHD ) - Hydro_RiemannSolver_Exact( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, + Hydro_RiemannSolver_Exact( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, PassiveFloor, EoS->DensEint2Pres_FuncPtr, EoS->DensPres2CSqr_FuncPtr, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table ); # elif ( RSOLVER == ROE ) - Hydro_RiemannSolver_Roe ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, + Hydro_RiemannSolver_Roe ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, PassiveFloor, EoS->DensEint2Pres_FuncPtr, EoS->DensPres2CSqr_FuncPtr, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table ); # elif ( RSOLVER == HLLE ) - Hydro_RiemannSolver_HLLE ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, + Hydro_RiemannSolver_HLLE ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, PassiveFloor, EoS->DensEint2Pres_FuncPtr, EoS->DensPres2CSqr_FuncPtr, EoS->GuessHTilde_FuncPtr, EoS->HTilde2Temp_FuncPtr, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table ); # elif ( RSOLVER == HLLC && !defined MHD ) - Hydro_RiemannSolver_HLLC ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, + Hydro_RiemannSolver_HLLC ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, PassiveFloor, EoS->DensEint2Pres_FuncPtr, EoS->DensPres2CSqr_FuncPtr, EoS->GuessHTilde_FuncPtr, EoS->HTilde2Temp_FuncPtr, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table ); # elif ( RSOLVER == HLLD && defined MHD ) - Hydro_RiemannSolver_HLLD ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, + Hydro_RiemannSolver_HLLD ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, PassiveFloor, EoS->DensEint2Pres_FuncPtr, EoS->DensPres2CSqr_FuncPtr, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table ); # else @@ -740,25 +745,25 @@ void Hydro_RiemannPredict_Flux( const real g_ConVar[][ CUBE(FLU_NXT) ], # endif # if ( RSOLVER_RESCUE == EXACT && !defined MHD ) - Hydro_RiemannSolver_Exact( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, + Hydro_RiemannSolver_Exact( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, PassiveFloor, EoS->DensEint2Pres_FuncPtr, EoS->DensPres2CSqr_FuncPtr, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table ); # elif ( RSOLVER_RESCUE == ROE ) - Hydro_RiemannSolver_Roe ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, + Hydro_RiemannSolver_Roe ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, PassiveFloor, EoS->DensEint2Pres_FuncPtr, EoS->DensPres2CSqr_FuncPtr, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table ); # elif ( RSOLVER_RESCUE == HLLE ) - Hydro_RiemannSolver_HLLE ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, + Hydro_RiemannSolver_HLLE ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, PassiveFloor, EoS->DensEint2Pres_FuncPtr, EoS->DensPres2CSqr_FuncPtr, EoS->GuessHTilde_FuncPtr, EoS->HTilde2Temp_FuncPtr, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table ); # elif ( RSOLVER_RESCUE == HLLC && !defined MHD ) - Hydro_RiemannSolver_HLLC ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, + Hydro_RiemannSolver_HLLC ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, PassiveFloor, EoS->DensEint2Pres_FuncPtr, EoS->DensPres2CSqr_FuncPtr, EoS->GuessHTilde_FuncPtr, EoS->HTilde2Temp_FuncPtr, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table ); # elif ( RSOLVER_RESCUE == HLLD && defined MHD ) - Hydro_RiemannSolver_HLLD ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, + Hydro_RiemannSolver_HLLD ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, PassiveFloor, EoS->DensEint2Pres_FuncPtr, EoS->DensPres2CSqr_FuncPtr, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table ); # else @@ -814,6 +819,7 @@ void Hydro_RiemannPredict_Flux( const real g_ConVar[][ CUBE(FLU_NXT) ], // dt : Time interval to advance solution // dh : Cell size // MinDens/Pres/Eint : Density, pressure, and internal energy floors +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored // FracPassive : true --> convert passive scalars to mass fraction during data reconstruction // NFrac : Number of passive scalars for the option "FracPassive" // FracIdx : Target variable indices for the option "FracPassive" @@ -828,7 +834,7 @@ void Hydro_RiemannPredict( const real g_ConVar_In[][ CUBE(FLU_NXT) ], real g_PriVar_Half[][ CUBE(FLU_NXT) ], const real dt, const real dh, const real MinDens, const real MinPres, const real MinEint, - const bool FracPassive, const int NFrac, const int FracIdx[], + const long PassiveFloor, const bool FracPassive, const int NFrac, const int FracIdx[], const bool JeansMinPres, const real JeansMinPres_Coeff, const EoS_t *EoS ) { @@ -909,16 +915,18 @@ void Hydro_RiemannPredict( const real g_ConVar_In[][ CUBE(FLU_NXT) ], # else const real Emag = NULL_REAL; # endif - out_con[4] = Hydro_CheckMinEintInEngy( out_con[0], out_con[1], out_con[2], out_con[3], out_con[4], MinEint, Emag ); + out_con[4] = Hydro_CheckMinEintInEngy( out_con[0], out_con[1], out_con[2], out_con[3], out_con[4], MinEint, PassiveFloor, Emag ); # endif // #ifndef BAROTROPIC_EOS # endif // #ifndef SRHD # if ( NCOMP_PASSIVE > 0 ) for (int v=NCOMP_FLUID; v primitive variables - Hydro_Con2Pri( out_con, out_pri, MinPres, FracPassive, NFrac, FracIdx, JeansMinPres, JeansMinPres_Coeff, + Hydro_Con2Pri( out_con, out_pri, MinPres, PassiveFloor, + FracPassive, NFrac, FracIdx, JeansMinPres, JeansMinPres_Coeff, EoS->DensEint2Pres_FuncPtr, EoS->DensPres2Eint_FuncPtr, EoS->GuessHTilde_FuncPtr, EoS->HTilde2Temp_FuncPtr, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, diff --git a/src/Model_Hydro/CPU_Hydro/CPU_FluidSolver_RTVD.cpp b/src/Model_Hydro/CPU_Hydro/CPU_FluidSolver_RTVD.cpp index 114c90fa84..8abdfa36d6 100644 --- a/src/Model_Hydro/CPU_Hydro/CPU_FluidSolver_RTVD.cpp +++ b/src/Model_Hydro/CPU_Hydro/CPU_FluidSolver_RTVD.cpp @@ -15,7 +15,7 @@ static void CPU_AdvanceX( real u[][ CUBE(FLU_NXT) ], const real dt, const real dx, const bool StoreFlux, const int j_skip, const int k_skip, const real MinDens, const real MinPres, const real MinEint, - const EoS_t *EoS ); + const long PassiveFloor, const EoS_t *EoS ); static void TransposeXY( real u[][ FLU_NXT*FLU_NXT*FLU_NXT ] ); static void TransposeXZ( real u[][ FLU_NXT*FLU_NXT*FLU_NXT ] ); @@ -43,6 +43,7 @@ static void TransposeXZ( real u[][ FLU_NXT*FLU_NXT*FLU_NXT ] ); // MinDens : Density floor // MinPres : Pressure floor // MinEint : Internal energy floor +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored // EoS : EoS object //------------------------------------------------------------------------------------------------------- void CPU_FluidSolver_RTVD( @@ -54,7 +55,7 @@ void CPU_FluidSolver_RTVD( const int NPatchGroup, const real dt, const real dh, const bool StoreFlux, const bool XYZ, const real MinDens, const real MinPres, const real MinEint, - const EoS_t EoS ) + const long PassiveFloor, const EoS_t EoS ) { if ( XYZ ) @@ -62,15 +63,15 @@ void CPU_FluidSolver_RTVD( # pragma omp parallel for schedule( runtime ) for (int P=0; P store the coarse-fine fluxes -// j_skip : Number of cells that can be skipped on each side in the y direction -// k_skip : Number of cells that can be skipped on each side in the z direction -// MinDens : Density floor -// MinPres : Pressure floor -// MinEint : Internal energy floor -// EoS : EoS object +// Parameter : u : Input fluid array +// dt : Time interval to advance solution +// dx : Grid size +// StoreFlux : true --> store the coarse-fine fluxes +// j_skip : Number of cells that can be skipped on each side in the y direction +// k_skip : Number of cells that can be skipped on each side in the z direction +// MinDens : Density floor +// MinPres : Pressure floor +// MinEint : Internal energy floor +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored +// EoS : EoS object //------------------------------------------------------------------------------------------------------- void CPU_AdvanceX( real u[][ CUBE(FLU_NXT) ], const real dt, const real dx, const bool StoreFlux, const int j_skip, const int k_skip, const real MinDens, const real MinPres, const real MinEint, - const EoS_t *EoS ) + const long PassiveFloor, const EoS_t *EoS ) { const bool CheckMinPres_Yes = true; @@ -227,7 +229,7 @@ void CPU_AdvanceX( real u[][ CUBE(FLU_NXT) ], const real dt, const real dx, _rho = (real)1.0 / ux[0][i]; vx = _rho * ux[1][i]; p = Hydro_Con2Pres( ux[0][i], ux[1][i], ux[2][i], ux[3][i], ux[4][i], Passive, - CheckMinPres_Yes, MinPres, NULL_REAL, EoS->DensEint2Pres_FuncPtr, NULL, NULL, + CheckMinPres_Yes, MinPres, PassiveFloor, NULL_REAL, EoS->DensEint2Pres_FuncPtr, NULL, NULL, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table, NULL ); # ifdef CHECK_UNPHYSICAL_IN_FLUID @@ -275,7 +277,7 @@ void CPU_AdvanceX( real u[][ CUBE(FLU_NXT) ], const real dt, const real dx, { u_half[0][i] = FMAX( u_half[0][i], MinDens ); u_half[4][i] = Hydro_CheckMinEintInEngy( u_half[0][i], u_half[1][i], u_half[2][i], u_half[3][i], u_half[4][i], - MinEint, NULL_REAL ); + MinEint, PassiveFloor, NULL_REAL ); } @@ -289,7 +291,7 @@ void CPU_AdvanceX( real u[][ CUBE(FLU_NXT) ], const real dt, const real dx, _rho = (real)1.0 / u_half[0][i]; vx = _rho * u_half[1][i]; p = Hydro_Con2Pres( u_half[0][i], u_half[1][i], u_half[2][i], u_half[3][i], u_half[4][i], Passive, - CheckMinPres_Yes, MinPres, NULL_REAL, EoS->DensEint2Pres_FuncPtr, NULL, NULL, + CheckMinPres_Yes, MinPres, PassiveFloor, NULL_REAL, EoS->DensEint2Pres_FuncPtr, NULL, NULL, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table, NULL ); # ifdef CHECK_UNPHYSICAL_IN_FLUID @@ -369,7 +371,7 @@ void CPU_AdvanceX( real u[][ CUBE(FLU_NXT) ], const real dt, const real dx, { ux[0][i] = FMAX( ux[0][i], MinDens ); ux[4][i] = Hydro_CheckMinEintInEngy( ux[0][i], ux[1][i], ux[2][i], ux[3][i], ux[4][i], - MinEint, NULL_REAL ); + MinEint, PassiveFloor, NULL_REAL ); } diff --git a/src/Model_Hydro/CPU_Hydro/CPU_Shared_ComputeFlux.cpp b/src/Model_Hydro/CPU_Hydro/CPU_Shared_ComputeFlux.cpp index 3a01ed6834..e153de8a05 100644 --- a/src/Model_Hydro/CPU_Hydro/CPU_Shared_ComputeFlux.cpp +++ b/src/Model_Hydro/CPU_Hydro/CPU_Shared_ComputeFlux.cpp @@ -32,19 +32,19 @@ #if ( RSOLVER == EXACT || RSOLVER_RESCUE == EXACT ) void Hydro_RiemannSolver_Exact( const int XYZ, real Flux_Out[], const real L_In[], const real R_In[], - const real MinDens, const real MinPres, const EoS_DE2P_t EoS_DensEint2Pres, + const real MinDens, const real MinPres, const long PassiveFloor, const EoS_DE2P_t EoS_DensEint2Pres, const EoS_DP2C_t EoS_DensPres2CSqr, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], const real* const EoS_Table[EOS_NTABLE_MAX] ); #endif #if ( RSOLVER == ROE || RSOLVER_RESCUE == ROE ) void Hydro_RiemannSolver_Roe( const int XYZ, real Flux_Out[], const real L_In[], const real R_In[], - const real MinDens, const real MinPres, const EoS_DE2P_t EoS_DensEint2Pres, + const real MinDens, const real MinPres, const long PassiveFloor, const EoS_DE2P_t EoS_DensEint2Pres, const EoS_DP2C_t EoS_DensPres2CSqr, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], const real* const EoS_Table[EOS_NTABLE_MAX] ); #endif #if ( RSOLVER == HLLE || RSOLVER_RESCUE == HLLE ) void Hydro_RiemannSolver_HLLE( const int XYZ, real Flux_Out[], const real L_In[], const real R_In[], - const real MinDens, const real MinPres, const EoS_DE2P_t EoS_DensEint2Pres, + const real MinDens, const real MinPres, const long PassiveFloor, const EoS_DE2P_t EoS_DensEint2Pres, const EoS_DP2C_t EoS_DensPres2CSqr, const EoS_GUESS_t EoS_GuessHTilde, const EoS_H2TEM_t EoS_HTilde2Temp, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], @@ -52,7 +52,7 @@ void Hydro_RiemannSolver_HLLE( const int XYZ, real Flux_Out[], const real L_In[] #endif #if ( RSOLVER == HLLC || RSOLVER_RESCUE == HLLC ) void Hydro_RiemannSolver_HLLC( const int XYZ, real Flux_Out[], const real L_In[], const real R_In[], - const real MinDens, const real MinPres, const EoS_DE2P_t EoS_DensEint2Pres, + const real MinDens, const real MinPres, const long PassiveFloor, const EoS_DE2P_t EoS_DensEint2Pres, const EoS_DP2C_t EoS_DensPres2CSqr, const EoS_GUESS_t EoS_GuessHTilde, const EoS_H2TEM_t EoS_HTilde2Temp, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], @@ -60,7 +60,7 @@ void Hydro_RiemannSolver_HLLC( const int XYZ, real Flux_Out[], const real L_In[] #endif #if ( RSOLVER == HLLD || RSOLVER_RESCUE == HLLD ) void Hydro_RiemannSolver_HLLD( const int XYZ, real Flux_Out[], const real L_In[], const real R_In[], - const real MinDens, const real MinPres, const EoS_DE2P_t EoS_DensEint2Pres, + const real MinDens, const real MinPres, const long PassiveFloor, const EoS_DE2P_t EoS_DensEint2Pres, const EoS_DP2C_t EoS_DensPres2CSqr, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], const real* const EoS_Table[EOS_NTABLE_MAX] ); #endif @@ -106,6 +106,7 @@ void Hydro_RiemannSolver_HLLD( const int XYZ, real Flux_Out[], const real L_In[] // ExtAcc_Func : Function pointer to the external acceleration routine (for UNSPLIT_GRAVITY only) // ExtAcc_AuxArray : Auxiliary array for external acceleration (for UNSPLIT_GRAVITY only) // MinDens/Pres : Density and pressure floors +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored // EoS : EoS object //------------------------------------------------------------------------------------------------------- GPU_DEVICE @@ -115,7 +116,7 @@ void Hydro_ComputeFlux( const real g_FC_Var [][NCOMP_TOTAL_PLUS_MAG][ CUBE(N_FC_ const bool CorrHalfVel, const real g_Pot_USG[], const double g_Corner[], const real dt, const real dh, const double Time, const bool UsePot, const OptExtAcc_t ExtAcc, const ExtAcc_t ExtAcc_Func, const double ExtAcc_AuxArray[], - const real MinDens, const real MinPres, const EoS_t *EoS ) + const real MinDens, const real MinPres, const long PassiveFloor, const EoS_t *EoS ) { // check @@ -282,25 +283,25 @@ void Hydro_ComputeFlux( const real g_FC_Var [][NCOMP_TOTAL_PLUS_MAG][ CUBE(N_FC_ // 2. invoke Riemann solver # if ( RSOLVER == EXACT && !defined MHD ) - Hydro_RiemannSolver_Exact( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, + Hydro_RiemannSolver_Exact( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, PassiveFloor, EoS->DensEint2Pres_FuncPtr, EoS->DensPres2CSqr_FuncPtr, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table ); # elif ( RSOLVER == ROE ) - Hydro_RiemannSolver_Roe ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, + Hydro_RiemannSolver_Roe ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, PassiveFloor, EoS->DensEint2Pres_FuncPtr, EoS->DensPres2CSqr_FuncPtr, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table ); # elif ( RSOLVER == HLLE ) - Hydro_RiemannSolver_HLLE ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, + Hydro_RiemannSolver_HLLE ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, PassiveFloor, EoS->DensEint2Pres_FuncPtr, EoS->DensPres2CSqr_FuncPtr, EoS->GuessHTilde_FuncPtr, EoS->HTilde2Temp_FuncPtr, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table ); # elif ( RSOLVER == HLLC && !defined MHD ) - Hydro_RiemannSolver_HLLC ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, + Hydro_RiemannSolver_HLLC ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, PassiveFloor, EoS->DensEint2Pres_FuncPtr, EoS->DensPres2CSqr_FuncPtr, EoS->GuessHTilde_FuncPtr, EoS->HTilde2Temp_FuncPtr, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table ); # elif ( RSOLVER == HLLD && defined MHD ) - Hydro_RiemannSolver_HLLD ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, + Hydro_RiemannSolver_HLLD ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, PassiveFloor, EoS->DensEint2Pres_FuncPtr, EoS->DensPres2CSqr_FuncPtr, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table ); # else @@ -320,25 +321,25 @@ void Hydro_ComputeFlux( const real g_FC_Var [][NCOMP_TOTAL_PLUS_MAG][ CUBE(N_FC_ # endif # if ( RSOLVER_RESCUE == EXACT && !defined MHD ) - Hydro_RiemannSolver_Exact( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, + Hydro_RiemannSolver_Exact( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, PassiveFloor, EoS->DensEint2Pres_FuncPtr, EoS->DensPres2CSqr_FuncPtr, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table ); # elif ( RSOLVER_RESCUE == ROE ) - Hydro_RiemannSolver_Roe ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, + Hydro_RiemannSolver_Roe ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, PassiveFloor, EoS->DensEint2Pres_FuncPtr, EoS->DensPres2CSqr_FuncPtr, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table ); # elif ( RSOLVER_RESCUE == HLLE ) - Hydro_RiemannSolver_HLLE ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, + Hydro_RiemannSolver_HLLE ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, PassiveFloor, EoS->DensEint2Pres_FuncPtr, EoS->DensPres2CSqr_FuncPtr, EoS->GuessHTilde_FuncPtr, EoS->HTilde2Temp_FuncPtr, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table ); # elif ( RSOLVER_RESCUE == HLLC && !defined MHD ) - Hydro_RiemannSolver_HLLC ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, + Hydro_RiemannSolver_HLLC ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, PassiveFloor, EoS->DensEint2Pres_FuncPtr, EoS->DensPres2CSqr_FuncPtr, EoS->GuessHTilde_FuncPtr, EoS->HTilde2Temp_FuncPtr, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table ); # elif ( RSOLVER_RESCUE == HLLD && defined MHD ) - Hydro_RiemannSolver_HLLD ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, + Hydro_RiemannSolver_HLLD ( d, Flux_1Face, ConVar_L, ConVar_R, MinDens, MinPres, PassiveFloor, EoS->DensEint2Pres_FuncPtr, EoS->DensPres2CSqr_FuncPtr, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table ); # else diff --git a/src/Model_Hydro/CPU_Hydro/CPU_Shared_DataReconstruction.cpp b/src/Model_Hydro/CPU_Hydro/CPU_Shared_DataReconstruction.cpp index ce95f70fc1..a93fb329a7 100644 --- a/src/Model_Hydro/CPU_Hydro/CPU_Shared_DataReconstruction.cpp +++ b/src/Model_Hydro/CPU_Hydro/CPU_Shared_DataReconstruction.cpp @@ -21,7 +21,7 @@ #else void Hydro_Rotate3D( real InOut[], const int XYZ, const bool Forward, const int Mag_Offset ); -void Hydro_Con2Pri( const real In[], real Out[], const real MinPres, +void Hydro_Con2Pri( const real In[], real Out[], const real MinPres, const long PassiveFloor, const bool FracPassive, const int NFrac, const int FracIdx[], const bool JeansMinPres, const real JeansMinPres_Coeff, const EoS_DE2P_t EoS_DensEint2Pres, const EoS_DP2E_t EoS_DensPres2Eint, @@ -33,7 +33,7 @@ void Hydro_Pri2Con( const real In[], real Out[], const bool FracPassive, const i const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], const real *const EoS_Table[EOS_NTABLE_MAX], const real* const EintIn ); #if ( FLU_SCHEME == MHM ) -void Hydro_Con2Flux( const int XYZ, real Flux[], const real In[], const real MinPres, +void Hydro_Con2Flux( const int XYZ, real Flux[], const real In[], const real MinPres, const long PassiveFloor, const EoS_DE2P_t EoS_DensEint2Pres, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], const real *const EoS_Table[EOS_NTABLE_MAX], const real* const PresIn ); #ifdef MHD @@ -80,12 +80,12 @@ static void Hydro_HancockPredict( real fcCon[][NCOMP_LR], const real fcPri[][NCO const real g_EC_Ele[][ CUBE(N_EC_ELE) ], const int NGhost, const int NEle, const real MinDens, const real MinPres, const real MinEint, - const EoS_t *EoS ); + const EoS_t *EoS, const long PassiveFloor ); #ifdef MHD GPU_DEVICE void Hydro_ConFC2PriCC_MHM( real g_PriVar[][ CUBE(FLU_NXT) ], const real g_FC_Var [][NCOMP_TOTAL_PLUS_MAG][ CUBE(N_FC_VAR) ], - const real MinDens, const real MinPres, const real MinEint, + const real MinDens, const real MinPres, const real MinEint, const long PassiveFloor, const bool FracPassive, const int NFrac, const int FracIdx[], const bool JeansMinPres, const real JeansMinPres_Coeff, const EoS_t *EoS ); @@ -177,6 +177,7 @@ static void Hydro_Char2Pri( real InOut[], const real Dens, const real Pres, cons // dt : Time interval to advance solution (for the CTU scheme) // dh : Cell size // MinDens/Pres/Eint : Density, pressure, and internal energy floors +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored // FracPassive : true --> convert passive scalars to mass fraction during data reconstruction // NFrac : Number of passive scalars for the option "FracPassive" // FracIdx : Target variable indices for the option "FracPassive" @@ -194,7 +195,7 @@ void Hydro_DataReconstruction( const real g_ConVar [][ CUBE(FLU_NXT) ], const bool Con2Pri, const LR_Limiter_t LR_Limiter, const real MinMod_Coeff, const real dt, const real dh, const real MinDens, const real MinPres, const real MinEint, - const bool FracPassive, const int NFrac, const int FracIdx[], + const long PassiveFloor, const bool FracPassive, const int NFrac, const int FracIdx[], const bool JeansMinPres, const real JeansMinPres_Coeff, const EoS_t *EoS ) { @@ -319,7 +320,7 @@ void Hydro_DataReconstruction( const real g_ConVar [][ CUBE(FLU_NXT) ], MHD_GetCellCenteredBField( ConVar_1Cell+NCOMP_TOTAL, g_FC_B[0], g_FC_B[1], g_FC_B[2], NIn, NIn, NIn, i, j, k ); # endif - Hydro_Con2Pri( ConVar_1Cell, PriVar_1Cell, MinPres, FracPassive, NFrac, FracIdx, + Hydro_Con2Pri( ConVar_1Cell, PriVar_1Cell, MinPres, PassiveFloor, FracPassive, NFrac, FracIdx, JeansMinPres, JeansMinPres_Coeff, EoS->DensEint2Pres_FuncPtr, EoS->DensPres2Eint_FuncPtr, EoS->GuessHTilde_FuncPtr, EoS->HTilde2Temp_FuncPtr, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table, EintPtr, NULL ); @@ -612,9 +613,10 @@ void Hydro_DataReconstruction( const real g_ConVar [][ CUBE(FLU_NXT) ], fcPri[faceR][4] = Hydro_CheckMinPres( fcPri[faceR][4], MinPres ); # if ( NCOMP_PASSIVE > 0 ) - for (int v=NCOMP_FLUID; v must be done after the CGPU_LOOP( idx_fc, CUBE(N_FC_VAR) ) loop since it will update g_PriVar[] - Hydro_ConFC2PriCC_MHM( g_PriVar, g_FC_Var, MinDens, MinPres, MinEint, FracPassive, NFrac, FracIdx, + Hydro_ConFC2PriCC_MHM( g_PriVar, g_FC_Var, MinDens, MinPres, MinEint, PassiveFloor, FracPassive, NFrac, FracIdx, JeansMinPres, JeansMinPres_Coeff, EoS ); # endif // #if ( FLU_SCHEME == MHM && defined MHD ) @@ -708,7 +710,7 @@ void Hydro_DataReconstruction( const real g_ConVar [][ CUBE(FLU_NXT) ], const bool Con2Pri, const LR_Limiter_t LR_Limiter, const real MinMod_Coeff, const real dt, const real dh, const real MinDens, const real MinPres, const real MinEint, - const bool FracPassive, const int NFrac, const int FracIdx[], + const long PassiveFloor, const bool FracPassive, const int NFrac, const int FracIdx[], const bool JeansMinPres, const real JeansMinPres_Coeff, const EoS_t *EoS ) { @@ -847,7 +849,7 @@ void Hydro_DataReconstruction( const real g_ConVar [][ CUBE(FLU_NXT) ], MHD_GetCellCenteredBField( ConVar_1Cell+NCOMP_TOTAL, g_FC_B[0], g_FC_B[1], g_FC_B[2], NIn, NIn, NIn, i, j, k ); # endif - Hydro_Con2Pri( ConVar_1Cell, PriVar_1Cell, MinPres, FracPassive, NFrac, FracIdx, + Hydro_Con2Pri( ConVar_1Cell, PriVar_1Cell, MinPres, PassiveFloor, FracPassive, NFrac, FracIdx, JeansMinPres, JeansMinPres_Coeff, EoS->DensEint2Pres_FuncPtr, EoS->DensPres2Eint_FuncPtr, EoS->GuessHTilde_FuncPtr, EoS->HTilde2Temp_FuncPtr, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table, EintPtr, NULL ); @@ -1319,9 +1321,10 @@ void Hydro_DataReconstruction( const real g_ConVar [][ CUBE(FLU_NXT) ], fcPri[faceR][4] = Hydro_CheckMinPres( fcPri[faceR][4], MinPres ); # if ( NCOMP_PASSIVE > 0 ) - for (int v=NCOMP_FLUID; v must be done after the CGPU_LOOP( idx_fc, CUBE(N_FC_VAR) ) loop since it will update g_PriVar[] - Hydro_ConFC2PriCC_MHM( g_PriVar, g_FC_Var, MinDens, MinPres, MinEint, FracPassive, NFrac, FracIdx, + Hydro_ConFC2PriCC_MHM( g_PriVar, g_FC_Var, MinDens, MinPres, MinEint, PassiveFloor, FracPassive, NFrac, FracIdx, JeansMinPres, JeansMinPres_Coeff, EoS ); # endif // #if ( FLU_SCHEME == MHM && defined MHD ) @@ -2055,6 +2058,7 @@ void Hydro_LimitSlope( const real L[], const real C[], const real R[], const LR_ // NEle : Stride for accessing g_EC_Ele[] // MinDens/Pres/Eint : Density, pressure, and internal energy floors // EoS : EoS object +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored //------------------------------------------------------------------------------------------------------- GPU_DEVICE void Hydro_HancockPredict( real fcCon[][NCOMP_LR], const real fcPri[][NCOMP_LR], const real dt, @@ -2064,7 +2068,7 @@ void Hydro_HancockPredict( real fcCon[][NCOMP_LR], const real fcPri[][NCOMP_LR], const real g_EC_Ele[][ CUBE(N_EC_ELE) ], const int NGhost, const int NEle, const real MinDens, const real MinPres, const real MinEint, - const EoS_t *EoS ) + const EoS_t *EoS, const long PassiveFloor ) { const real dt_dh2 = (real)0.5*dt/dh; @@ -2074,10 +2078,10 @@ void Hydro_HancockPredict( real fcCon[][NCOMP_LR], const real fcPri[][NCOMP_LR], // calculate flux for (int f=0; f<6; f++) # ifdef SRHD - Hydro_Con2Flux( f/2, Flux[f], fcCon[f], MinPres, EoS->DensEint2Pres_FuncPtr, + Hydro_Con2Flux( f/2, Flux[f], fcCon[f], MinPres, PassiveFloor, EoS->DensEint2Pres_FuncPtr, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table, fcPri[f] ); # else - Hydro_Con2Flux( f/2, Flux[f], fcCon[f], MinPres, EoS->DensEint2Pres_FuncPtr, + Hydro_Con2Flux( f/2, Flux[f], fcCon[f], MinPres, PassiveFloor, EoS->DensEint2Pres_FuncPtr, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table, NULL ); # endif @@ -2106,7 +2110,7 @@ void Hydro_HancockPredict( real fcCon[][NCOMP_LR], const real fcPri[][NCOMP_LR], EoS->DensEint2Pres_FuncPtr, EoS->GuessHTilde_FuncPtr, EoS->HTilde2Temp_FuncPtr, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table, - ERROR_INFO, UNPHY_SILENCE ) ) + PassiveFloor, ERROR_INFO, UNPHY_SILENCE ) ) reset_cell = true; # else @@ -2143,11 +2147,12 @@ void Hydro_HancockPredict( real fcCon[][NCOMP_LR], const real fcPri[][NCOMP_LR], const real Emag = NULL_REAL; # endif // MHD fcCon[f][4] = Hydro_CheckMinEintInEngy( fcCon[f][0], fcCon[f][1], fcCon[f][2], fcCon[f][3], fcCon[f][4], - MinEint, Emag ); + MinEint, PassiveFloor, Emag ); # endif // #ifndef BAROTROPIC_EOS # endif // #ifndef SRHD # if ( NCOMP_PASSIVE > 0 ) - for (int v=NCOMP_FLUID; v Should contain NCOMP_TOTAL_PLUS_MAG variables // MinDens/Pres/Eint : Density, pressure, and internal energy floors +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored // FracPassive : true --> convert passive scalars to mass fraction during data reconstruction // NFrac : Number of passive scalars for the option "FracPassive" // FracIdx : Target variable indices for the option "FracPassive" @@ -2182,7 +2188,7 @@ void Hydro_HancockPredict( real fcCon[][NCOMP_LR], const real fcPri[][NCOMP_LR], GPU_DEVICE void Hydro_ConFC2PriCC_MHM( real g_PriVar[][ CUBE(FLU_NXT) ], const real g_FC_Var [][NCOMP_TOTAL_PLUS_MAG][ CUBE(N_FC_VAR) ], - const real MinDens, const real MinPres, const real MinEint, + const real MinDens, const real MinPres, const real MinEint, const long PassiveFloor, const bool FracPassive, const int NFrac, const int FracIdx[], const bool JeansMinPres, const real JeansMinPres_Coeff, const EoS_t *EoS ) @@ -2207,7 +2213,7 @@ void Hydro_ConFC2PriCC_MHM( real g_PriVar[][ CUBE(FLU_NXT) ], g_FC_Var[faceR][MAG_OFFSET+d][idx_fc] ); } - Hydro_Con2Pri( ConCC, PriCC, MinPres, FracPassive, NFrac, FracIdx, + Hydro_Con2Pri( ConCC, PriCC, MinPres, PassiveFloor, FracPassive, NFrac, FracIdx, JeansMinPres, JeansMinPres_Coeff, EoS->DensEint2Pres_FuncPtr, EoS->DensPres2Eint_FuncPtr, EoS->GuessHTilde_FuncPtr, EoS->HTilde2Temp_FuncPtr, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table, NULL, NULL ); diff --git a/src/Model_Hydro/CPU_Hydro/CPU_Shared_DualEnergy.cpp b/src/Model_Hydro/CPU_Hydro/CPU_Shared_DualEnergy.cpp index 9cda80175d..6e187ebc9d 100644 --- a/src/Model_Hydro/CPU_Hydro/CPU_Shared_DualEnergy.cpp +++ b/src/Model_Hydro/CPU_Hydro/CPU_Shared_DualEnergy.cpp @@ -53,6 +53,7 @@ static real Hydro_DensDual2Pres( const real Dens, const real Dual, const real Ga // --> In some cases we actually want to check if pressure becomes unphysical, // for which we don't want to enable this option // MinPres : Minimum allowed pressure +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored // DualEnergySwitch : if ( Eint/(Ekin+Emag) < DualEnergySwitch ) ==> correct Eint and Etot // else ==> correct Dual // Emag : Magnetic energy density (0.5*B^2) --> for MHD only @@ -62,7 +63,7 @@ static real Hydro_DensDual2Pres( const real Dens, const real Dual, const real Ga GPU_DEVICE void Hydro_DualEnergyFix( const real Dens, const real MomX, const real MomY, const real MomZ, real &Etot, real &Dual, char &DE_Status, const real Gamma_m1, const real _Gamma_m1, - const bool CheckMinPres, const real MinPres, const real DualEnergySwitch, + const bool CheckMinPres, const real MinPres, const long PassiveFloor, const real DualEnergySwitch, const real Emag ) { @@ -78,7 +79,7 @@ void Hydro_DualEnergyFix( const real Dens, const real MomX, const real MomY, con // --> Enth (i.e., non-thermal energy) includes both kinetic and magnetic energies real Enth, Eint, Pres; - Eint = Hydro_Con2Eint( Dens, MomX, MomY, MomZ, Etot, CheckMinEint_No, NULL_REAL, Emag, + Eint = Hydro_Con2Eint( Dens, MomX, MomY, MomZ, Etot, CheckMinEint_No, NULL_REAL, PassiveFloor, Emag, NULL, NULL, NULL, NULL, NULL ); Enth = Etot - Eint; @@ -138,12 +139,14 @@ void Hydro_DualEnergyFix( const real Dens, const real MomX, const real MomY, con // EoS_DensEint2Pres : EoS routine to compute the gas pressure // EoS_AuxArray_* : Auxiliary arrays for EoS_DensEint2Pres() // EoS_Table : EoS tables +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored // // Return : Dual //------------------------------------------------------------------------------------------------------- real Hydro_Con2Dual( const real Dens, const real MomX, const real MomY, const real MomZ, const real Engy, const real Emag, const EoS_DE2P_t EoS_DensEint2Pres, const double EoS_AuxArray_Flt[], - const int EoS_AuxArray_Int[], const real *const EoS_Table[EOS_NTABLE_MAX] ) + const int EoS_AuxArray_Int[], const real *const EoS_Table[EOS_NTABLE_MAX], + const long PassiveFloor ) { // currently this function does NOT apply pressure floor when calling Hydro_Con2Pres() @@ -153,7 +156,7 @@ real Hydro_Con2Dual( const real Dens, const real MomX, const real MomY, const re // calculate pressure and convert it to the dual-energy variable // --> note that DE_ENPY only works with EOS_GAMMA, which does not involve passive scalars - Pres = Hydro_Con2Pres( Dens, MomX, MomY, MomZ, Engy, NULL, CheckMinPres_No, NULL_REAL, Emag, + Pres = Hydro_Con2Pres( Dens, MomX, MomY, MomZ, Engy, NULL, CheckMinPres_No, NULL_REAL, PassiveFloor, Emag, EoS_DensEint2Pres, NULL, NULL, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table, NULL ); Dual = Hydro_DensPres2Dual( Dens, Pres, EoS_AuxArray_Flt[1] ); diff --git a/src/Model_Hydro/CPU_Hydro/CPU_Shared_FluUtility.cpp b/src/Model_Hydro/CPU_Hydro/CPU_Shared_FluUtility.cpp index 51cbfe378f..d2e7907be7 100644 --- a/src/Model_Hydro/CPU_Hydro/CPU_Shared_FluUtility.cpp +++ b/src/Model_Hydro/CPU_Hydro/CPU_Shared_FluUtility.cpp @@ -14,13 +14,13 @@ #ifdef __CUDACC__ GPU_DEVICE static real Hydro_Con2Pres( const real Dens, const real MomX, const real MomY, const real MomZ, const real Engy, - const real Passive[], const bool CheckMinPres, const real MinPres, const real Emag, + const real Passive[], const bool CheckMinPres, const real MinPres, const long PassiveFloor, const real Emag, const EoS_DE2P_t EoS_DensEint2Pres, const EoS_GUESS_t EoS_GuessHTilde, const EoS_H2TEM_t EoS_HTilde2Temp, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], const real *const EoS_Table[EOS_NTABLE_MAX], real *EintOut ); GPU_DEVICE static real Hydro_Con2Eint( const real Dens, const real MomX, const real MomY, const real MomZ, const real Engy, - const bool CheckMinEint, const real MinEint, const real Emag, + const bool CheckMinEint, const real MinEint, const long PassiveFloor, const real Emag, const EoS_GUESS_t EoS_GuessHTilde, const EoS_H2TEM_t EoS_HTilde2Temp, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], const real *const EoS_Table[EOS_NTABLE_MAX] ); @@ -40,7 +40,7 @@ static bool Hydro_IsUnphysical( const IsUnphyMode_t Mode, const real Fields[], const real Emag, const EoS_DE2P_t EoS_DensEint2Pres, const EoS_GUESS_t EoS_GuessHTilde, const EoS_H2TEM_t EoS_HTilde2Temp, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], - const real *const EoS_Table[EOS_NTABLE_MAX], + const real *const EoS_Table[EOS_NTABLE_MAX], const long PassiveFloor, const char File[], const int Line, const char Function[], const IsUnphVerb_t Verbose ); GPU_DEVICE static bool Hydro_IsUnphysical_Single( const real Field, const char SingleFieldName[], const real Min, const real Max, @@ -184,6 +184,7 @@ void Hydro_Rotate3D( real InOut[], const int XYZ, const bool Forward, const int // Parameter : In : Input conserved variables // Out : Output primitive variables // MinPres : Minimum allowed pressure +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored // FracPassive : true --> convert passive scalars to mass fraction // NFrac : Number of passive scalars for the option "FracPassive" // FracIdx : Target variable indices for the option "FracPassive" @@ -205,7 +206,7 @@ void Hydro_Rotate3D( real InOut[], const int XYZ, const bool Forward, const int // Return : Out[], EintOut (optional), LorentzFactorPtr (optional) //------------------------------------------------------------------------------------------------------- GPU_DEVICE -void Hydro_Con2Pri( const real In[], real Out[], const real MinPres, +void Hydro_Con2Pri( const real In[], real Out[], const real MinPres, const long PassiveFloor, const bool FracPassive, const int NFrac, const int FracIdx[], const bool JeansMinPres, const real JeansMinPres_Coeff, const EoS_DE2P_t EoS_DensEint2Pres, const EoS_DP2E_t EoS_DensPres2Eint, @@ -222,7 +223,7 @@ void Hydro_Con2Pri( const real In[], real Out[], const real MinPres, Hydro_IsUnphysical( UNPHY_MODE_CONS, In, NULL_REAL, EoS_DensEint2Pres, EoS_GuessHTilde, EoS_HTilde2Temp, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table, - ERROR_INFO, UNPHY_VERBOSE ); + PassiveFloor, ERROR_INFO, UNPHY_VERBOSE ); # endif HTilde = Hydro_Con2HTilde( In, EoS_GuessHTilde, EoS_HTilde2Temp, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table ); @@ -264,7 +265,7 @@ void Hydro_Con2Pri( const real In[], real Out[], const real MinPres, Out[1] = In[1]*_Rho; Out[2] = In[2]*_Rho; Out[3] = In[3]*_Rho; - Out[4] = Hydro_Con2Pres( In[0], In[1], In[2], In[3], In[4], In+NCOMP_FLUID, CheckMinPres_Yes, MinPres, Emag, + Out[4] = Hydro_Con2Pres( In[0], In[1], In[2], In[3], In[4], In+NCOMP_FLUID, CheckMinPres_Yes, MinPres, PassiveFloor, Emag, EoS_DensEint2Pres, EoS_GuessHTilde, EoS_HTilde2Temp, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table, EintOut ); @@ -442,6 +443,7 @@ void Hydro_Pri2Con( const real In[], real Out[], const bool FracPassive, // Flux : Array to store the output fluxes // In : Array storing the input conserved variables // MinPres : Minimum allowed pressure +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored // EoS_DensEint2Pres : EoS routine to compute the gas pressure // EoS_AuxArray_* : Auxiliary arrays for EoS_DensEint2Pres() // EoS_Table : EoS tables for EoS_DensEint2Pres() @@ -453,7 +455,7 @@ void Hydro_Pri2Con( const real In[], real Out[], const bool FracPassive, // Return : Flux[] //------------------------------------------------------------------------------------------------------- GPU_DEVICE -void Hydro_Con2Flux( const int XYZ, real Flux[], const real In[], const real MinPres, +void Hydro_Con2Flux( const int XYZ, real Flux[], const real In[], const real MinPres, const long PassiveFloor, const EoS_DE2P_t EoS_DensEint2Pres, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], const real *const EoS_Table[EOS_NTABLE_MAX], const real* const AuxArray ) { @@ -507,7 +509,7 @@ void Hydro_Con2Flux( const int XYZ, real Flux[], const real In[], const real Min # else // #ifdef SRHD const real Pres = ( AuxArray == NULL ) ? Hydro_Con2Pres( InRot[0], InRot[1], InRot[2], InRot[3], InRot[4], In+NCOMP_FLUID, - CheckMinPres_Yes, MinPres, Emag, EoS_DensEint2Pres, + CheckMinPres_Yes, MinPres, PassiveFloor, Emag, EoS_DensEint2Pres, NULL, NULL, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table, NULL ) : AuxArray[0]; @@ -772,23 +774,24 @@ real Hydro_CheckMinEntr( const real InEntr, const real MinEntr ) // 2. Input conserved instead of primitive variables // 3. For MHD, one must provide the magnetic energy density Emag (i.e., 0.5*B^2) // -// Parameter : Dens : Mass density -// MomX/Y/Z : Momentum density -// InEngy : Energy density -// MinEint : Internal energy density floor -// Emag : Magnetic energy density (0.5*B^2) --> For MHD only +// Parameter : Dens : Mass density +// MomX/Y/Z : Momentum density +// InEngy : Energy density +// MinEint : Internal energy density floor +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored +// Emag : Magnetic energy density (0.5*B^2) --> For MHD only // // Return : Total energy density with internal energy density greater than a given threshold //------------------------------------------------------------------------------------------------------- GPU_DEVICE real Hydro_CheckMinEintInEngy( const real Dens, const real MomX, const real MomY, const real MomZ, const real InEngy, - const real MinEint, const real Emag ) + const real MinEint, const long PassiveFloor, const real Emag ) { const bool CheckMinEint_No = false; real InEint, OutEint, OutEngy; - InEint = Hydro_Con2Eint( Dens, MomX, MomY, MomZ, InEngy, CheckMinEint_No, NULL_REAL, Emag, + InEint = Hydro_Con2Eint( Dens, MomX, MomY, MomZ, InEngy, CheckMinEint_No, NULL_REAL, PassiveFloor, Emag, NULL, NULL, NULL, NULL, NULL ); OutEint = Hydro_CheckMinEint( InEint, MinEint ); @@ -826,6 +829,7 @@ real Hydro_CheckMinEintInEngy( const real Dens, const real MomX, const real MomY // Fields : Field data to be checked // Emag : Magnetic energy density (0.5*B^2) --> For MHD only // EoS_* : EoS parameters +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored // File : __FILE__ // Line : __LINE__ // Function : __FUNCTION__ @@ -840,7 +844,7 @@ bool Hydro_IsUnphysical( const IsUnphyMode_t Mode, const real Fields[], const real Emag, const EoS_DE2P_t EoS_DensEint2Pres, const EoS_GUESS_t EoS_GuessHTilde, const EoS_H2TEM_t EoS_HTilde2Temp, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], - const real *const EoS_Table[EOS_NTABLE_MAX], + const real *const EoS_Table[EOS_NTABLE_MAX], const long PassiveFloor, const char File[], const int Line, const char Function[], const IsUnphVerb_t Verbose ) { @@ -881,7 +885,9 @@ bool Hydro_IsUnphysical( const IsUnphyMode_t Mode, const real Fields[], // check passive scalars (which can be zero) else { - if ( Fields[v] < (real)0.0 || Fields[v] > HUGE_NUMBER ) + if ( Fields[v] < (real)0.0 && PassiveFloor & BIDX(v) ) + UnphyCell = true; + if ( Fields[v] < -HUGE_NUMBER || Fields[v] > HUGE_NUMBER ) UnphyCell = true; } } // for (int v=0; v HUGE_NUMBER || Eint != Eint ) @@ -996,7 +1002,9 @@ bool Hydro_IsUnphysical( const IsUnphyMode_t Mode, const real Fields[], // check passive scalars (which can be zero) else { - if ( Fields[v] < (real)0.0 || Fields[v] > HUGE_NUMBER ) + if ( Fields[v] < (real)0.0 && PassiveFloor & BIDX(v) ) + UnphyCell = true; + if ( Fields[v] < -HUGE_NUMBER || Fields[v] > HUGE_NUMBER ) UnphyCell = true; } } // for (int v=0; v HUGE_NUMBER ) +// check negative (passive scalars can be zero) + if ( Fields[v] < (real)0.0 && PassiveFloor & BIDX(v+NCOMP_FLUID) ) + UnphyCell = true; + +// check infinity + if ( Fields[v] < -HUGE_NUMBER || Fields[v] > HUGE_NUMBER ) UnphyCell = true; } @@ -1136,6 +1148,7 @@ bool Hydro_IsUnphysical_Single( const real Field, const char SingleFieldName[], // for which this option should be disabled // --> For example: Flu_FixUp(), Flu_Close(), Hydro_Aux_Check_Negative() // MinPres : Pressure floor +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored // Emag : Magnetic energy density (0.5*B^2) --> For MHD only // EoS_DensEint2Pres : EoS routine to compute the gas pressure // EoS_GuessHTilde : EoS routine to compute guessed reduced enthalpy @@ -1150,7 +1163,7 @@ bool Hydro_IsUnphysical_Single( const real Field, const char SingleFieldName[], //------------------------------------------------------------------------------------------------------- GPU_DEVICE real Hydro_Con2Pres( const real Dens, const real MomX, const real MomY, const real MomZ, const real Engy, - const real Passive[], const bool CheckMinPres, const real MinPres, const real Emag, + const real Passive[], const bool CheckMinPres, const real MinPres, const long PassiveFloor, const real Emag, const EoS_DE2P_t EoS_DensEint2Pres, const EoS_GUESS_t EoS_GuessHTilde, const EoS_H2TEM_t EoS_HTilde2Temp, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], const real *const EoS_Table[EOS_NTABLE_MAX], real *EintOut ) @@ -1167,7 +1180,7 @@ real Hydro_Con2Pres( const real Dens, const real MomX, const real MomY, const re Cons[3] = MomZ; Cons[4] = Engy; - Hydro_Con2Pri( Cons, Prim, (CheckMinPres)?MinPres:-HUGE_NUMBER, false, NULL_INT, NULL, + Hydro_Con2Pri( Cons, Prim, (CheckMinPres)?MinPres:-HUGE_NUMBER, PassiveFloor, false, NULL_INT, NULL, NULL_BOOL, NULL_REAL, NULL, NULL, EoS_GuessHTilde, EoS_HTilde2Temp, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table, NULL, NULL ); Pres = Prim[4]; @@ -1177,7 +1190,7 @@ real Hydro_Con2Pres( const real Dens, const real MomX, const real MomY, const re const bool CheckMinEint_No = false; real Eint; - Eint = Hydro_Con2Eint( Dens, MomX, MomY, MomZ, Engy, CheckMinEint_No, NULL_REAL, Emag, + Eint = Hydro_Con2Eint( Dens, MomX, MomY, MomZ, Engy, CheckMinEint_No, NULL_REAL, PassiveFloor, Emag, NULL, NULL, NULL, NULL, NULL ); Pres = EoS_DensEint2Pres( Dens, Eint, Passive, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table ); @@ -1233,6 +1246,7 @@ real Hydro_Con2Pres( const real Dens, const real MomX, const real MomY, const re // --> In some cases we actually want to check if internal energy // becomes unphysical, for which this option should be disabled // MinEint : Internal energy floor +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored // Emag : Magnetic energy density (0.5*B^2) --> For MHD only // EoS_GuessHTilde : EoS routine to compute guessed reduced enthalpy // EoS_HTilde2Temp : EoS routine to compute temperature @@ -1243,7 +1257,7 @@ real Hydro_Con2Pres( const real Dens, const real MomX, const real MomY, const re //------------------------------------------------------------------------------------------------------- GPU_DEVICE real Hydro_Con2Eint( const real Dens, const real MomX, const real MomY, const real MomZ, const real Engy, - const bool CheckMinEint, const real MinEint, const real Emag, + const bool CheckMinEint, const real MinEint, const long PassiveFloor, const real Emag, const EoS_GUESS_t EoS_GuessHTilde, const EoS_H2TEM_t EoS_HTilde2Temp, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], const real *const EoS_Table[EOS_NTABLE_MAX] ) @@ -1261,7 +1275,7 @@ real Hydro_Con2Eint( const real Dens, const real MomX, const real MomY, const re Cons[3] = MomZ; Cons[4] = Engy; - Hydro_Con2Pri( Cons, Prim, -HUGE_NUMBER, false, NULL_INT, NULL, + Hydro_Con2Pri( Cons, Prim, -HUGE_NUMBER, PassiveFloor, false, NULL_INT, NULL, NULL_BOOL, NULL_REAL, NULL, NULL, EoS_GuessHTilde, EoS_HTilde2Temp, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table, NULL, NULL ); @@ -1345,6 +1359,7 @@ real Hydro_ConEint2Etot( const real Dens, const real MomX, const real MomY, cons // --> In some cases we actually want to check if temperature becomes unphysical, // for which we don't want to enable this option // MinTemp : Temperature floor +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored // Emag : Magnetic energy density (0.5*B^2) --> For MHD only // EoS_DensEint2Temp : EoS routine to compute the gas temperature // EoS_GuessHTilde : EoS routine to compute guessed reduced enthalpy @@ -1356,7 +1371,7 @@ real Hydro_ConEint2Etot( const real Dens, const real MomX, const real MomY, cons //------------------------------------------------------------------------------------------------------- GPU_DEVICE real Hydro_Con2Temp( const real Dens, const real MomX, const real MomY, const real MomZ, const real Engy, - const real Passive[], const bool CheckMinTemp, const real MinTemp, const real Emag, + const real Passive[], const bool CheckMinTemp, const real MinTemp, const long PassiveFloor, const real Emag, const EoS_DE2T_t EoS_DensEint2Temp, const EoS_GUESS_t EoS_GuessHTilde, const EoS_H2TEM_t EoS_HTilde2Temp, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], const real *const EoS_Table[EOS_NTABLE_MAX] ) @@ -1401,7 +1416,7 @@ real Hydro_Con2Temp( const real Dens, const real MomX, const real MomY, const re Cons[3] = MomZ; Cons[4] = Engy; - Hydro_Con2Pri( Cons, Prim, -HUGE_NUMBER, false, NULL_INT, NULL, + Hydro_Con2Pri( Cons, Prim, -HUGE_NUMBER, PassiveFloor, false, NULL_INT, NULL, NULL_BOOL, NULL_REAL, NULL, NULL, EoS_GuessHTilde, EoS_HTilde2Temp, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table, NULL, NULL ); Temp = Prim[4]/Prim[0]; @@ -1412,7 +1427,7 @@ real Hydro_Con2Temp( const real Dens, const real MomX, const real MomY, const re const bool CheckMinEint_No = false; real Eint; - Eint = Hydro_Con2Eint( Dens, MomX, MomY, MomZ, Engy, CheckMinEint_No, NULL_REAL, Emag, + Eint = Hydro_Con2Eint( Dens, MomX, MomY, MomZ, Engy, CheckMinEint_No, NULL_REAL, PassiveFloor, Emag, NULL, NULL, NULL, NULL, NULL ); Temp = EoS_DensEint2Temp( Dens, Eint, Passive, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table ); # endif // #ifdef SRHD ... else ... @@ -1445,6 +1460,7 @@ real Hydro_Con2Temp( const real Dens, const real MomX, const real MomY, const re // --> In some cases we actually want to check if entropy becomes unphysical, // for which we don't want to enable this option // MinEntr : Entropy floor +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored // Emag : Magnetic energy density (0.5*B^2) --> For MHD only // EoS_DensEint2Entr : EoS routine to compute the gas entropy // EoS_AuxArray_* : Auxiliary arrays for EoS_DensEint2Entr() @@ -1454,7 +1470,7 @@ real Hydro_Con2Temp( const real Dens, const real MomX, const real MomY, const re //------------------------------------------------------------------------------------------------------- GPU_DEVICE real Hydro_Con2Entr( const real Dens, const real MomX, const real MomY, const real MomZ, const real Engy, - const real Passive[], const bool CheckMinEntr, const real MinEntr, const real Emag, + const real Passive[], const bool CheckMinEntr, const real MinEntr, const long PassiveFloor, const real Emag, const EoS_DE2S_t EoS_DensEint2Entr, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], const real *const EoS_Table[EOS_NTABLE_MAX] ) { @@ -1485,7 +1501,7 @@ real Hydro_Con2Entr( const real Dens, const real MomX, const real MomY, const re const bool CheckMinEint_No = false; real Eint, Entr; - Eint = Hydro_Con2Eint( Dens, MomX, MomY, MomZ, Engy, CheckMinEint_No, NULL_REAL, Emag, + Eint = Hydro_Con2Eint( Dens, MomX, MomY, MomZ, Engy, CheckMinEint_No, NULL_REAL, PassiveFloor, Emag, NULL, NULL, NULL, NULL, NULL ); Entr = EoS_DensEint2Entr( Dens, Eint, Passive, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table ); diff --git a/src/Model_Hydro/CPU_Hydro/CPU_Shared_FullStepUpdate.cpp b/src/Model_Hydro/CPU_Hydro/CPU_Shared_FullStepUpdate.cpp index d48a10faee..83dce5da58 100644 --- a/src/Model_Hydro/CPU_Hydro/CPU_Shared_FullStepUpdate.cpp +++ b/src/Model_Hydro/CPU_Hydro/CPU_Shared_FullStepUpdate.cpp @@ -46,6 +46,7 @@ // dh : Cell size // MinDens/Eint : Density and internal energy floors // DualEnergySwitch : Use the dual-energy formalism if E_int/E_kin < DualEnergySwitch +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored // NormPassive : true --> normalize passive scalars so that the sum of their mass density // is equal to the gas mass density // NNorm : Number of passive scalars to be normalized @@ -62,8 +63,8 @@ GPU_DEVICE void Hydro_FullStepUpdate( const real g_Input[][ CUBE(FLU_NXT) ], real g_Output[][ CUBE(PS2) ], char g_DE_Status[], const real g_FC_B[][ PS2P1*SQR(PS2) ], const real g_Flux[][NCOMP_TOTAL_PLUS_MAG][ CUBE(N_FC_FLUX) ], - const real dt, const real dh, const real MinDens, const real MinEint, - const real DualEnergySwitch, const bool NormPassive, const int NNorm, const int NormIdx[], + const real dt, const real dh, const real MinDens, const real MinEint, const real DualEnergySwitch, + const long PassiveFloor, const bool NormPassive, const int NNorm, const int NormIdx[], const EoS_t *EoS, int *s_FullStepFailure, const int Iteration, const int MinMod_MaxIter ) { @@ -136,13 +137,14 @@ void Hydro_FullStepUpdate( const real g_Input[][ CUBE(FLU_NXT) ], real g_Output[ // Output_1Cell[DENS] = FMAX( Output_1Cell[DENS], MinDens ); Output_1Cell[ENGY] = Hydro_CheckMinEintInEngy( Output_1Cell[DENS], Output_1Cell[MOMX], Output_1Cell[MOMY], Output_1Cell[MOMZ], - Output_1Cell[ENGY], MinEint, Emag ); + Output_1Cell[ENGY], MinEint, PassiveFloor, Emag ); # endif // #ifdef BAROTROPIC_EOS // 2. floor and normalize passive scalars # if ( NCOMP_PASSIVE > 0 ) - for (int v=NCOMP_FLUID; vAuxArrayDevPtr_Flt[1], EoS->AuxArrayDevPtr_Flt[2], CheckMinPres_No, NULL_REAL, - DualEnergySwitch, Emag ); + PassiveFloor, DualEnergySwitch, Emag ); # endif // #ifdef DUAL_ENERGY @@ -182,7 +184,7 @@ void Hydro_FullStepUpdate( const real g_Input[][ CUBE(FLU_NXT) ], real g_Output[ EoS->DensEint2Pres_FuncPtr, EoS->GuessHTilde_FuncPtr, EoS->HTilde2Temp_FuncPtr, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table, - ERROR_INFO, UNPHY_SILENCE ) ) + PassiveFloor, ERROR_INFO, UNPHY_SILENCE ) ) { # ifdef __CUDACC__ // GPU // use atomicExch_block() on Pascal (or later) GPUs to avoid inter-block synchronization for better performance @@ -207,7 +209,7 @@ void Hydro_FullStepUpdate( const real g_Input[][ CUBE(FLU_NXT) ], real g_Output[ { // get pressure const real Pres = Hydro_Con2Pres( Output_1Cell[DENS], Output_1Cell[MOMX], Output_1Cell[MOMY], Output_1Cell[MOMZ], - Output_1Cell[ENGY], Output_1Cell+NCOMP_FLUID, CheckMinPres_No, NULL_REAL, Emag, + Output_1Cell[ENGY], Output_1Cell+NCOMP_FLUID, CheckMinPres_No, NULL_REAL, PassiveFloor, Emag, EoS->DensEint2Pres_FuncPtr, EoS->GuessHTilde_FuncPtr, EoS->HTilde2Temp_FuncPtr, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, diff --git a/src/Model_Hydro/CPU_Hydro/CPU_Shared_RiemannSolver_Exact.cpp b/src/Model_Hydro/CPU_Hydro/CPU_Shared_RiemannSolver_Exact.cpp index 88345b70e4..430d7637f7 100644 --- a/src/Model_Hydro/CPU_Hydro/CPU_Shared_RiemannSolver_Exact.cpp +++ b/src/Model_Hydro/CPU_Hydro/CPU_Shared_RiemannSolver_Exact.cpp @@ -18,7 +18,7 @@ #else // #ifdef __CUDACC__ -void Hydro_Con2Pri( const real In[], real Out[], const real MinPres, +void Hydro_Con2Pri( const real In[], real Out[], const real MinPres, const long PassiveFloor, const bool FracPassive, const int NFrac, const int FracIdx[], const bool JeansMinPres, const real JeansMinPres_Coeff, const EoS_DE2P_t EoS_DensEint2Pres, const EoS_DP2E_t EoS_DensPres2Eint, @@ -51,6 +51,7 @@ GPU_DEVICE static void Set_Flux( real flux[], const real val[], const real Gamma // Flux_Out : Output array to store the average flux along t axis // L/R_In : Input left/right states (conserved variables) // MinDens/Pres : Density and pressure floors +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored // EoS_DensEint2Pres : EoS routine to compute the gas pressure // EoS_DensPres2CSqr : EoS routine to compute the sound speed squared // EoS_AuxArray_* : Auxiliary arrays for the EoS routines @@ -60,7 +61,7 @@ GPU_DEVICE static void Set_Flux( real flux[], const real val[], const real Gamma //------------------------------------------------------------------------------------------------------ GPU_DEVICE void Hydro_RiemannSolver_Exact( const int XYZ, real Flux_Out[], const real L_In[], const real R_In[], - const real MinDens, const real MinPres, const EoS_DE2P_t EoS_DensEint2Pres, + const real MinDens, const real MinPres, const long PassiveFloor, const EoS_DE2P_t EoS_DensEint2Pres, const EoS_DP2C_t EoS_DensPres2CSqr, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], const real* const EoS_Table[EOS_NTABLE_MAX] ) { @@ -90,9 +91,9 @@ void Hydro_RiemannSolver_Exact( const int XYZ, real Flux_Out[], const real L_In[ real L[NCOMP_TOTAL], R[NCOMP_TOTAL], Temp; // convert conserved variables to primitive variables - Hydro_Con2Pri( L_In, L, MinPres, FracPassive_No, NULL_INT, NULL, JeansMinPres_No, NULL_REAL, + Hydro_Con2Pri( L_In, L, MinPres, PassiveFloor, FracPassive_No, NULL_INT, NULL, JeansMinPres_No, NULL_REAL, EoS_DensEint2Pres, NULL, NULL, NULL, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table, NULL, NULL ); - Hydro_Con2Pri( R_In, R, MinPres, FracPassive_No, NULL_INT, NULL, JeansMinPres_No, NULL_REAL, + Hydro_Con2Pri( R_In, R, MinPres, PassiveFloor, FracPassive_No, NULL_INT, NULL, JeansMinPres_No, NULL_REAL, EoS_DensEint2Pres, NULL, NULL, NULL, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table, NULL, NULL ); diff --git a/src/Model_Hydro/CPU_Hydro/CPU_Shared_RiemannSolver_HLLC.cpp b/src/Model_Hydro/CPU_Hydro/CPU_Shared_RiemannSolver_HLLC.cpp index e6c1c1c34b..010b59d872 100644 --- a/src/Model_Hydro/CPU_Hydro/CPU_Shared_RiemannSolver_HLLC.cpp +++ b/src/Model_Hydro/CPU_Hydro/CPU_Shared_RiemannSolver_HLLC.cpp @@ -17,7 +17,7 @@ #else // #ifdef __CUDACC__ void Hydro_Rotate3D( real InOut[], const int XYZ, const bool Forward, const int Mag_Offset ); -void Hydro_Con2Flux( const int XYZ, real Flux[], const real In[], const real MinPres, +void Hydro_Con2Flux( const int XYZ, real Flux[], const real In[], const real MinPres, const long PassiveFloor, const EoS_DE2P_t EoS_DensEint2Pres, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], const real *const EoS_Table[EOS_NTABLE_MAX], const real* const PresIn ); @@ -48,6 +48,7 @@ void Hydro_Con2Flux( const int XYZ, real Flux[], const real In[], const real Min // Flux_Out : Array to store the output flux // L/R_In : Input left/right states (conserved variables) // MinDens/Pres : Density and pressure floors +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored // EoS_DensEint2Pres : EoS routine to compute the gas pressure // EoS_DensPres2CSqr : EoS routine to compute the sound speed squared // EoS_AuxArray_* : Auxiliary arrays for the EoS routines @@ -57,7 +58,7 @@ void Hydro_Con2Flux( const int XYZ, real Flux[], const real In[], const real Min //------------------------------------------------------------------------------------------------------- GPU_DEVICE void Hydro_RiemannSolver_HLLC( const int XYZ, real Flux_Out[], const real L_In[], const real R_In[], - const real MinDens, const real MinPres, const EoS_DE2P_t EoS_DensEint2Pres, + const real MinDens, const real MinPres, const long PassiveFloor, const EoS_DE2P_t EoS_DensEint2Pres, const EoS_DP2C_t EoS_DensPres2CSqr, const EoS_GUESS_t EoS_GuessHTilde, const EoS_H2TEM_t EoS_HTilde2Temp, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], @@ -80,11 +81,11 @@ void Hydro_RiemannSolver_HLLC( const int XYZ, real Flux_Out[], const real L_In[] Hydro_IsUnphysical( UNPHY_MODE_CONS, L, NULL_REAL, EoS_DensEint2Pres, EoS_GuessHTilde, EoS_HTilde2Temp, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table, - ERROR_INFO, UNPHY_VERBOSE ); + PassiveFloor, ERROR_INFO, UNPHY_VERBOSE ); Hydro_IsUnphysical( UNPHY_MODE_CONS, R, NULL_REAL, EoS_DensEint2Pres, EoS_GuessHTilde, EoS_HTilde2Temp, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table, - ERROR_INFO, UNPHY_VERBOSE ); + PassiveFloor, ERROR_INFO, UNPHY_VERBOSE ); # endif @@ -101,11 +102,11 @@ void Hydro_RiemannSolver_HLLC( const int XYZ, real Flux_Out[], const real L_In[] real lV1, rV1, lV2, rV2, lV3, rV3; real lFactor,rFactor; - Hydro_Con2Pri( L, PL, MinPres, NULL_BOOL, NULL_INT, NULL, NULL_BOOL, + Hydro_Con2Pri( L, PL, MinPres, PassiveFloor, NULL_BOOL, NULL_INT, NULL, NULL_BOOL, (real)NULL_REAL, NULL, NULL, EoS_GuessHTilde, EoS_HTilde2Temp, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table, NULL, &lFactor ); - Hydro_Con2Pri( R, PR, MinPres, NULL_BOOL, NULL_INT, NULL, NULL_BOOL, + Hydro_Con2Pri( R, PR, MinPres, PassiveFloor, NULL_BOOL, NULL_INT, NULL, NULL_BOOL, (real)NULL_REAL, NULL, NULL, EoS_GuessHTilde, EoS_HTilde2Temp, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table, NULL, &rFactor ); @@ -113,11 +114,11 @@ void Hydro_RiemannSolver_HLLC( const int XYZ, real Flux_Out[], const real L_In[] Hydro_IsUnphysical( UNPHY_MODE_PRIM, PL, NULL_REAL, EoS_DensEint2Pres, EoS_GuessHTilde, EoS_HTilde2Temp, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table, - ERROR_INFO, UNPHY_VERBOSE ); + PassiveFloor, ERROR_INFO, UNPHY_VERBOSE ); Hydro_IsUnphysical( UNPHY_MODE_PRIM, PR, NULL_REAL, EoS_DensEint2Pres, EoS_GuessHTilde, EoS_HTilde2Temp, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table, - ERROR_INFO, UNPHY_VERBOSE ); + PassiveFloor, ERROR_INFO, UNPHY_VERBOSE ); # endif @@ -389,7 +390,7 @@ void Hydro_RiemannSolver_HLLC( const int XYZ, real Flux_Out[], const real L_In[] //------------------------------------------------------------------------------------------------------- GPU_DEVICE void Hydro_RiemannSolver_HLLC( const int XYZ, real Flux_Out[], const real L_In[], const real R_In[], - const real MinDens, const real MinPres, const EoS_DE2P_t EoS_DensEint2Pres, + const real MinDens, const real MinPres, const long PassiveFloor, const EoS_DE2P_t EoS_DensEint2Pres, const EoS_DP2C_t EoS_DensPres2CSqr, const EoS_GUESS_t EoS_GuessHTilde, const EoS_H2TEM_t EoS_HTilde2Temp, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], @@ -430,10 +431,10 @@ void Hydro_RiemannSolver_HLLC( const int XYZ, real Flux_Out[], const real L_In[] _RhoR = ONE / R[0]; u_L = _RhoL*L[1]; u_R = _RhoR*R[1]; - P_L = Hydro_Con2Pres( L[0], L[1], L[2], L[3], L[4], L+NCOMP_FLUID, CheckMinPres_Yes, MinPres, Emag, + P_L = Hydro_Con2Pres( L[0], L[1], L[2], L[3], L[4], L+NCOMP_FLUID, CheckMinPres_Yes, MinPres, PassiveFloor, Emag, EoS_DensEint2Pres, EoS_GuessHTilde, EoS_HTilde2Temp, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table, NULL ); - P_R = Hydro_Con2Pres( R[0], R[1], R[2], R[3], R[4], R+NCOMP_FLUID, CheckMinPres_Yes, MinPres, Emag, + P_R = Hydro_Con2Pres( R[0], R[1], R[2], R[3], R[4], R+NCOMP_FLUID, CheckMinPres_Yes, MinPres, PassiveFloor, Emag, EoS_DensEint2Pres, EoS_GuessHTilde, EoS_HTilde2Temp, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table, NULL ); Cs_L = SQRT( EoS_DensPres2CSqr( L[0], P_L, L+NCOMP_FLUID, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table ) ); @@ -589,7 +590,7 @@ void Hydro_RiemannSolver_HLLC( const int XYZ, real Flux_Out[], const real L_In[] { const real MaxV_L = FMIN( W_L, ZERO ); - Hydro_Con2Flux( 0, Flux_LR, L, MinPres, NULL, NULL, NULL, NULL, &P_L ); + Hydro_Con2Flux( 0, Flux_LR, L, MinPres, PassiveFloor, NULL, NULL, NULL, NULL, &P_L ); for (int v=0; v= ZERO ) diff --git a/src/Model_Hydro/CPU_Hydro/CPU_Shared_RiemannSolver_HLLE.cpp b/src/Model_Hydro/CPU_Hydro/CPU_Shared_RiemannSolver_HLLE.cpp index f5f77d1405..78deae0dc6 100644 --- a/src/Model_Hydro/CPU_Hydro/CPU_Shared_RiemannSolver_HLLE.cpp +++ b/src/Model_Hydro/CPU_Hydro/CPU_Shared_RiemannSolver_HLLE.cpp @@ -17,7 +17,7 @@ #else // #ifdef __CUDACC__ void Hydro_Rotate3D( real InOut[], const int XYZ, const bool Forward, const int Mag_Offset ); -void Hydro_Con2Flux( const int XYZ, real Flux[], const real In[], const real MinPres, +void Hydro_Con2Flux( const int XYZ, real Flux[], const real In[], const real MinPres, const long PassiveFloor, const EoS_DE2P_t EoS_DensEint2Pres, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], const real *const EoS_Table[EOS_NTABLE_MAX], const real* const PresIn ); @@ -48,6 +48,7 @@ void Hydro_Con2Flux( const int XYZ, real Flux[], const real In[], const real Min // Flux_Out : Array to store the output flux // L/R_In : Input left/right states (conserved variables) // MinDens/Pres : Density and pressure floors +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored // EoS_DensEint2Pres : EoS routine to compute the gas pressure // EoS_DensPres2CSqr : EoS routine to compute the sound speed squared // EoS_AuxArray_* : Auxiliary arrays for the EoS routines @@ -57,7 +58,7 @@ void Hydro_Con2Flux( const int XYZ, real Flux[], const real In[], const real Min //------------------------------------------------------------------------------------------------------- GPU_DEVICE void Hydro_RiemannSolver_HLLE( const int XYZ, real Flux_Out[], const real L_In[], const real R_In[], - const real MinDens, const real MinPres, const EoS_DE2P_t EoS_DensEint2Pres, + const real MinDens, const real MinPres, const long PassiveFloor, const EoS_DE2P_t EoS_DensEint2Pres, const EoS_DP2C_t EoS_DensPres2CSqr, const EoS_GUESS_t EoS_GuessHTilde, const EoS_H2TEM_t EoS_HTilde2Temp, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], @@ -80,11 +81,11 @@ void Hydro_RiemannSolver_HLLE( const int XYZ, real Flux_Out[], const real L_In[] Hydro_IsUnphysical( UNPHY_MODE_CONS, L, NULL_REAL, EoS_DensEint2Pres, EoS_GuessHTilde, EoS_HTilde2Temp, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table, - ERROR_INFO, UNPHY_VERBOSE ); + PassiveFloor, ERROR_INFO, UNPHY_VERBOSE ); Hydro_IsUnphysical( UNPHY_MODE_CONS, R, NULL_REAL, EoS_DensEint2Pres, EoS_GuessHTilde, EoS_HTilde2Temp, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table, - ERROR_INFO, UNPHY_VERBOSE ); + PassiveFloor, ERROR_INFO, UNPHY_VERBOSE ); # endif @@ -100,11 +101,11 @@ void Hydro_RiemannSolver_HLLE( const int XYZ, real Flux_Out[], const real L_In[] real lFactor,rFactor; // Lorentz factor - Hydro_Con2Pri( L, PL, MinPres, NULL_BOOL, NULL_INT, NULL, NULL_BOOL, + Hydro_Con2Pri( L, PL, MinPres, PassiveFloor, NULL_BOOL, NULL_INT, NULL, NULL_BOOL, (real)NULL_REAL, NULL, NULL, EoS_GuessHTilde, EoS_HTilde2Temp, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table, NULL, &lFactor ); - Hydro_Con2Pri( R, PR, MinPres, NULL_BOOL, NULL_INT, NULL, NULL_BOOL, + Hydro_Con2Pri( R, PR, MinPres, PassiveFloor, NULL_BOOL, NULL_INT, NULL, NULL_BOOL, (real)NULL_REAL, NULL, NULL, EoS_GuessHTilde, EoS_HTilde2Temp, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table, NULL, &rFactor ); @@ -299,7 +300,7 @@ void Hydro_RiemannSolver_HLLE( const int XYZ, real Flux_Out[], const real L_In[] //------------------------------------------------------------------------------------------------------- GPU_DEVICE void Hydro_RiemannSolver_HLLE( const int XYZ, real Flux_Out[], const real L_In[], const real R_In[], - const real MinDens, const real MinPres, const EoS_DE2P_t EoS_DensEint2Pres, + const real MinDens, const real MinPres, const long PassiveFloor, const EoS_DE2P_t EoS_DensEint2Pres, const EoS_DP2C_t EoS_DensPres2CSqr, const EoS_GUESS_t EoS_GuessHTilde, const EoS_H2TEM_t EoS_HTilde2Temp, const double EoS_AuxArray_Flt[], const int EoS_AuxArray_Int[], @@ -365,10 +366,10 @@ void Hydro_RiemannSolver_HLLE( const int XYZ, real Flux_Out[], const real L_In[] Emag_R = NULL_REAL; # endif - P_L = Hydro_Con2Pres( L[0], L[1], L[2], L[3], L[4], L+NCOMP_FLUID, CheckMinPres_Yes, MinPres, Emag_L, + P_L = Hydro_Con2Pres( L[0], L[1], L[2], L[3], L[4], L+NCOMP_FLUID, CheckMinPres_Yes, MinPres, PassiveFloor, Emag_L, EoS_DensEint2Pres, EoS_GuessHTilde, EoS_HTilde2Temp, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table, NULL ); - P_R = Hydro_Con2Pres( R[0], R[1], R[2], R[3], R[4], R+NCOMP_FLUID, CheckMinPres_Yes, MinPres, Emag_R, + P_R = Hydro_Con2Pres( R[0], R[1], R[2], R[3], R[4], R+NCOMP_FLUID, CheckMinPres_Yes, MinPres, PassiveFloor, Emag_R, EoS_DensEint2Pres, EoS_GuessHTilde, EoS_HTilde2Temp, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table, NULL ); a2_L = EoS_DensPres2CSqr( L[0], P_L, L+NCOMP_FLUID, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table ); @@ -640,8 +641,8 @@ void Hydro_RiemannSolver_HLLE( const int XYZ, real Flux_Out[], const real L_In[] # endif real Flux_L[NCOMP_TOTAL_PLUS_MAG], Flux_R[NCOMP_TOTAL_PLUS_MAG]; // use NCOMP_TOTAL_PLUS_MAG for Hydro_Con2Flux() - Hydro_Con2Flux( 0, Flux_L, L, MinPres, NULL, NULL, NULL, NULL, &P_L ); - Hydro_Con2Flux( 0, Flux_R, R, MinPres, NULL, NULL, NULL, NULL, &P_R ); + Hydro_Con2Flux( 0, Flux_L, L, MinPres, PassiveFloor, NULL, NULL, NULL, NULL, &P_L ); + Hydro_Con2Flux( 0, Flux_R, R, MinPres, PassiveFloor, NULL, NULL, NULL, NULL, &P_R ); for (int v=0; v= ZERO ) @@ -629,7 +630,7 @@ void Hydro_RiemannSolver_Roe( const int XYZ, real Flux_Out[], const real L_In[], const real Emag = NULL_REAL; # endif I_Pres = Hydro_Con2Pres( I_States[0], I_States[1], I_States[2], I_States[3], I_States[4], Passive, - CheckMinPres_No, NULL_REAL, Emag, EoS_DensEint2Pres, NULL, NULL, + CheckMinPres_No, NULL_REAL, PassiveFloor, Emag, EoS_DensEint2Pres, NULL, NULL, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table, NULL ); // if unphysical results occur, recalculate fluxes by a substitute Riemann solver @@ -641,16 +642,16 @@ void Hydro_RiemannSolver_Roe( const int XYZ, real Flux_Out[], const real L_In[], # endif # if ( CHECK_INTERMEDIATE == EXACT && !defined MHD ) - Hydro_RiemannSolver_Exact( 0, Flux_Out, L, R, MinDens, MinPres, EoS_DensEint2Pres, EoS_DensPres2CSqr, + Hydro_RiemannSolver_Exact( 0, Flux_Out, L, R, MinDens, MinPres, PassiveFloor, EoS_DensEint2Pres, EoS_DensPres2CSqr, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table ); # elif ( CHECK_INTERMEDIATE == HLLE ) - Hydro_RiemannSolver_HLLE ( 0, Flux_Out, L, R, MinDens, MinPres, EoS_DensEint2Pres, EoS_DensPres2CSqr, + Hydro_RiemannSolver_HLLE ( 0, Flux_Out, L, R, MinDens, MinPres, PassiveFloor, EoS_DensEint2Pres, EoS_DensPres2CSqr, NULL, NULL, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table ); # elif ( CHECK_INTERMEDIATE == HLLC && !defined MHD ) - Hydro_RiemannSolver_HLLC ( 0, Flux_Out, L, R, MinDens, MinPres, EoS_DensEint2Pres, EoS_DensPres2CSqr, + Hydro_RiemannSolver_HLLC ( 0, Flux_Out, L, R, MinDens, MinPres, PassiveFloor, EoS_DensEint2Pres, EoS_DensPres2CSqr, NULL, NULL, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table ); # elif ( CHECK_INTERMEDIATE == HLLD && defined MHD ) - Hydro_RiemannSolver_HLLD ( 0, Flux_Out, L, R, MinDens, MinPres, EoS_DensEint2Pres, EoS_DensPres2CSqr, + Hydro_RiemannSolver_HLLD ( 0, Flux_Out, L, R, MinDens, MinPres, PassiveFloor, EoS_DensEint2Pres, EoS_DensPres2CSqr, EoS_AuxArray_Flt, EoS_AuxArray_Int, EoS_Table ); # else diff --git a/src/Model_Hydro/CPU_Hydro/CPU_dtSolver_HydroCFL.cpp b/src/Model_Hydro/CPU_Hydro/CPU_dtSolver_HydroCFL.cpp index 76cfe0e75f..f4a9e87746 100644 --- a/src/Model_Hydro/CPU_Hydro/CPU_dtSolver_HydroCFL.cpp +++ b/src/Model_Hydro/CPU_Hydro/CPU_dtSolver_HydroCFL.cpp @@ -39,15 +39,16 @@ // --> CFL condition // 3. Arrays with a prefix "g_" are stored in the global memory of GPU // -// Parameter : g_dt_Array : Array to store the minimum dt in each target patch -// g_Flu_Array : Array storing the prepared fluid data of each target patch -// g_Mag_Array : Array storing the prepared B field data of each target patch -// NPG : Number of target patch groups (for CPU only) -// dh : Cell size -// Safety : dt safety factor -// MinPres : Minimum allowed pressure -// EoS : EoS object -// MicroPhy : Microphysics object +// Parameter : g_dt_Array : Array to store the minimum dt in each target patch +// g_Flu_Array : Array storing the prepared fluid data of each target patch +// g_Mag_Array : Array storing the prepared B field data of each target patch +// NPG : Number of target patch groups (for CPU only) +// dh : Cell size +// Safety : dt safety factor +// MinPres : Minimum allowed pressure +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored +// EoS : EoS object +// MicroPhy : Microphysics object // // Return : g_dt_Array //----------------------------------------------------------------------------------------- @@ -55,13 +56,13 @@ __global__ void CUFLU_dtSolver_HydroCFL( real g_dt_Array[], const real g_Flu_Array[][FLU_NIN_T][ CUBE(PS1) ], const real g_Mag_Array[][NCOMP_MAG][ PS1P1*SQR(PS1) ], - const real dh, const real Safety, const real MinPres, const EoS_t EoS, - const MicroPhy_t MicroPhy ) + const real dh, const real Safety, const real MinPres, + const long PassiveFloor, const EoS_t EoS, const MicroPhy_t MicroPhy ) #else void CPU_dtSolver_HydroCFL ( real g_dt_Array[], const real g_Flu_Array[][FLU_NIN_T][ CUBE(PS1) ], const real g_Mag_Array[][NCOMP_MAG][ PS1P1*SQR(PS1) ], const int NPG, - const real dh, const real Safety, const real MinPres, const EoS_t EoS, - const MicroPhy_t MicroPhy ) + const real dh, const real Safety, const real MinPres, + const long PassiveFloor, const EoS_t EoS, const MicroPhy_t MicroPhy ) #endif { @@ -94,7 +95,7 @@ void CPU_dtSolver_HydroCFL ( real g_dt_Array[], const real g_Flu_Array[][FLU_NI # ifdef SRHD real Pri[FLU_NIN_T], LorentzFactor, U_Max, Us_Max, LorentzFactor_Max, LorentzFactor_s_Max, Us, Rho; - Hydro_Con2Pri( fluid, Pri, MinPres, NULL_BOOL, NULL_INT, NULL, NULL_BOOL, + Hydro_Con2Pri( fluid, Pri, MinPres, PassiveFloor, NULL_BOOL, NULL_INT, NULL, NULL_BOOL, (real)NULL_REAL, EoS.DensEint2Pres_FuncPtr, EoS.DensPres2Eint_FuncPtr, EoS.GuessHTilde_FuncPtr, EoS.HTilde2Temp_FuncPtr, EoS.AuxArrayDevPtr_Flt, EoS.AuxArrayDevPtr_Int, EoS.Table, NULL, &LorentzFactor ); @@ -129,7 +130,7 @@ void CPU_dtSolver_HydroCFL ( real g_dt_Array[], const real g_Flu_Array[][FLU_NI Vy = FABS( fluid[MOMY] )*_Rho; Vz = FABS( fluid[MOMZ] )*_Rho; Pres = Hydro_Con2Pres( fluid[DENS], fluid[MOMX], fluid[MOMY], fluid[MOMZ], fluid[ENGY], fluid+NCOMP_FLUID, - CheckMinPres_Yes, MinPres, Emag, + CheckMinPres_Yes, MinPres, PassiveFloor, Emag, EoS.DensEint2Pres_FuncPtr, EoS.GuessHTilde_FuncPtr, EoS.HTilde2Temp_FuncPtr, EoS.AuxArrayDevPtr_Flt, EoS.AuxArrayDevPtr_Int, EoS.Table, NULL ); a2 = EoS.DensPres2CSqr_FuncPtr( fluid[DENS], Pres, fluid+NCOMP_FLUID, EoS.AuxArrayDevPtr_Flt, EoS.AuxArrayDevPtr_Int, diff --git a/src/Model_Hydro/GPU_Hydro/CUFLU_FluidSolver_RTVD.cu b/src/Model_Hydro/GPU_Hydro/CUFLU_FluidSolver_RTVD.cu index db325d9595..773eff4f4d 100644 --- a/src/Model_Hydro/GPU_Hydro/CUFLU_FluidSolver_RTVD.cu +++ b/src/Model_Hydro/GPU_Hydro/CUFLU_FluidSolver_RTVD.cu @@ -22,7 +22,7 @@ static __device__ void CUFLU_Advance( real g_Fluid_In [][5][ CUBE(FLU_NXT) ], const int j_gap, const int k_gap, real s_cu[][5][FLU_NXT], real s_cw[][5][FLU_NXT], real s_flux[][5][FLU_NXT], real s_RLflux[][5][FLU_NXT], const bool FinalOut, const int XYZ, - const real MinDens, const real MinPres, const real MinEint, + const real MinDens, const real MinPres, const real MinEint, const long PassiveFloor, const EoS_t *EoS ); @@ -36,20 +36,21 @@ static __device__ void CUFLU_Advance( real g_Fluid_In [][5][ CUBE(FLU_NXT) ], // Prefix "s" for pointers pointing to the "Shared" memory space // b. The three-dimensional evolution is achieved by using the dimensional-split method // -// Parameter : g_Fluid_In : Global memory array to store the input fluid variables -// g_Fluid_Out : Global memory array to store the output fluid variables -// g_Flux : Global memory array to store the output fluxes -// g_Corner : Global memory array storing the physical corner coordinates of each patch group (USELESS CURRENTLY) -// g_Pot_USG : Global memory array storing the input potential for UNSPLIT_GRAVITY (NOT SUPPORTED in RTVD) -// dt : Time interval to advance solution -// _dh : 1 / grid size -// StoreFlux : true --> store the coarse-fine fluxes -// XYZ : true : x->y->z ( forward sweep) -// false : z->y->x (backward sweep) -// MinDens : Density floor -// MinPres : Pressure floor -// MinEint : Internal energy floor -// EoS : EoS object +// Parameter : g_Fluid_In : Global memory array to store the input fluid variables +// g_Fluid_Out : Global memory array to store the output fluid variables +// g_Flux : Global memory array to store the output fluxes +// g_Corner : Global memory array storing the physical corner coordinates of each patch group (USELESS CURRENTLY) +// g_Pot_USG : Global memory array storing the input potential for UNSPLIT_GRAVITY (NOT SUPPORTED in RTVD) +// dt : Time interval to advance solution +// _dh : 1 / grid size +// StoreFlux : true --> store the coarse-fine fluxes +// XYZ : true : x->y->z ( forward sweep) +// false : z->y->x (backward sweep) +// MinDens : Density floor +// MinPres : Pressure floor +// MinEint : Internal energy floor +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored +// EoS : EoS object //------------------------------------------------------------------------------------------------------- __global__ void CUFLU_FluidSolver_RTVD( real g_Fluid_In [][NCOMP_TOTAL][ CUBE(FLU_NXT) ], @@ -58,7 +59,7 @@ __global__ void CUFLU_FluidSolver_RTVD( const double g_Corner[][3], const real g_Pot_USG[][ CUBE(USG_NXT_F) ], const real dt, const real _dh, const bool StoreFlux, - const bool XYZ, const real MinDens, const real MinPres, const real MinEint, + const bool XYZ, const real MinDens, const real MinPres, const real MinEint, const long PassiveFloor, const EoS_t EoS ) { @@ -70,25 +71,25 @@ __global__ void CUFLU_FluidSolver_RTVD( if ( XYZ ) { CUFLU_Advance( g_Fluid_In, g_Fluid_Out, g_Flux, dt, _dh, StoreFlux, 0, 0, - s_cu, s_cw, s_flux, s_RLflux, false, 0, MinDens, MinPres, MinEint, &EoS ); + s_cu, s_cw, s_flux, s_RLflux, false, 0, MinDens, MinPres, MinEint, PassiveFloor, &EoS ); CUFLU_Advance( g_Fluid_In, g_Fluid_Out, g_Flux, dt, _dh, StoreFlux, FLU_GHOST_SIZE, 0, - s_cu, s_cw, s_flux, s_RLflux, false, 3, MinDens, MinPres, MinEint, &EoS ); + s_cu, s_cw, s_flux, s_RLflux, false, 3, MinDens, MinPres, MinEint, PassiveFloor, &EoS ); CUFLU_Advance( g_Fluid_In, g_Fluid_Out, g_Flux, dt, _dh, StoreFlux, FLU_GHOST_SIZE, FLU_GHOST_SIZE, - s_cu, s_cw, s_flux, s_RLflux, true, 6, MinDens, MinPres, MinEint, &EoS ); + s_cu, s_cw, s_flux, s_RLflux, true, 6, MinDens, MinPres, MinEint, PassiveFloor, &EoS ); } else { CUFLU_Advance( g_Fluid_In, g_Fluid_Out, g_Flux, dt, _dh, StoreFlux, 0, 0, - s_cu, s_cw, s_flux, s_RLflux, false, 6, MinDens, MinPres, MinEint, &EoS ); + s_cu, s_cw, s_flux, s_RLflux, false, 6, MinDens, MinPres, MinEint, PassiveFloor, &EoS ); CUFLU_Advance( g_Fluid_In, g_Fluid_Out, g_Flux, dt, _dh, StoreFlux, 0, FLU_GHOST_SIZE, - s_cu, s_cw, s_flux, s_RLflux, false, 3, MinDens, MinPres, MinEint, &EoS ); + s_cu, s_cw, s_flux, s_RLflux, false, 3, MinDens, MinPres, MinEint, PassiveFloor, &EoS ); CUFLU_Advance( g_Fluid_In, g_Fluid_Out, g_Flux, dt, _dh, StoreFlux, FLU_GHOST_SIZE, FLU_GHOST_SIZE, - s_cu, s_cw, s_flux, s_RLflux, true, 0, MinDens, MinPres, MinEint, &EoS ); + s_cu, s_cw, s_flux, s_RLflux, true, 0, MinDens, MinPres, MinEint, PassiveFloor, &EoS ); } } // FUNCTION : CUFLU_FluidSolver_RTVD @@ -103,26 +104,27 @@ __global__ void CUFLU_FluidSolver_RTVD( // Prefix "s" for pointers pointing to the "Shared" memory space // b. The direction of the one dimensional sweep is determined by the input parameter "XYZ" // -// Parameter : g_Fluid_In : Global memory array to store the input fluid variables -// g_Fluid_Out : Global memory array to store the output fluid variables -// g_Flux : Global memory array to store the output fluxes -// dt : Time interval to advance solution -// _dh : 1 / grid size -// StoreFlux : true --> store the coarse-fine fluxes -// j_gap : Number of useless grids in each side in the j direction (j may not be equal to y) -// k_gap : Number of useless grids in each side in the k direction (k mya not be equal to z) -// s_cu : Shared memory array storing the normal flux -// s_cw : Shared memory array storing the auxiliary flux -// s_flux : Shared memory array storing the final flux used to update the fluid variables -// s_RLflux : Shared memory array storing the left/right-moving flux -// XYZ : 0 : Update the solution in the x direction -// 3 : Update the solution in the y direction -// 6 : Update the solution in the z direction -// --> This parameter is also used to determine the place to store the output fluxes -// MinDens : Density floor -// MinPres : Pressure floor -// MinEint : Internal energy floor -// EoS : EoS object +// Parameter : g_Fluid_In : Global memory array to store the input fluid variables +// g_Fluid_Out : Global memory array to store the output fluid variables +// g_Flux : Global memory array to store the output fluxes +// dt : Time interval to advance solution +// _dh : 1 / grid size +// StoreFlux : true --> store the coarse-fine fluxes +// j_gap : Number of useless grids in each side in the j direction (j may not be equal to y) +// k_gap : Number of useless grids in each side in the k direction (k mya not be equal to z) +// s_cu : Shared memory array storing the normal flux +// s_cw : Shared memory array storing the auxiliary flux +// s_flux : Shared memory array storing the final flux used to update the fluid variables +// s_RLflux : Shared memory array storing the left/right-moving flux +// XYZ : 0 : Update the solution in the x direction +// 3 : Update the solution in the y direction +// 6 : Update the solution in the z direction +// --> This parameter is also used to determine the place to store the output fluxes +// MinDens : Density floor +// MinPres : Pressure floor +// MinEint : Internal energy floor +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored +// EoS : EoS object //------------------------------------------------------------------------------------------------------- __device__ void CUFLU_Advance( real g_Fluid_In [][5][ CUBE(FLU_NXT) ], real g_Fluid_Out[][5][ CUBE(PS2) ], @@ -131,7 +133,7 @@ __device__ void CUFLU_Advance( real g_Fluid_In [][5][ CUBE(FLU_NXT) ], const int j_gap, const int k_gap, real s_cu[][5][FLU_NXT], real s_cw[][5][FLU_NXT], real s_flux[][5][FLU_NXT], real s_RLflux[][5][FLU_NXT], const bool FinalOut, const int XYZ, - const real MinDens, const real MinPres, const real MinEint, + const real MinDens, const real MinPres, const real MinEint, const long PassiveFloor, const EoS_t *EoS ) { @@ -192,7 +194,7 @@ __device__ void CUFLU_Advance( real g_Fluid_In [][5][ CUBE(FLU_NXT) ], _rho = (real)1.0 / Fluid[0]; vx = _rho * Fluid[1]; p = Hydro_Con2Pres( Fluid[0], Fluid[1], Fluid[2], Fluid[3], Fluid[4], Passive, - CheckMinPres_Yes, MinPres, NULL_REAL, EoS->DensEint2Pres_FuncPtr, NULL, NULL, + CheckMinPres_Yes, MinPres, PassiveFloor, NULL_REAL, EoS->DensEint2Pres_FuncPtr, NULL, NULL, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table, NULL ); # ifdef CHECK_UNPHYSICAL_IN_FLUID @@ -237,7 +239,7 @@ __device__ void CUFLU_Advance( real g_Fluid_In [][5][ CUBE(FLU_NXT) ], // apply density and internal energy floors Fluid_half[0] = FMAX( Fluid_half[0], MinDens ); Fluid_half[4] = Hydro_CheckMinEintInEngy( Fluid_half[0], Fluid_half[1], Fluid_half[2], Fluid_half[3], Fluid_half[4], - MinEint, NULL_REAL ); + MinEint, PassiveFloor, NULL_REAL ); } @@ -251,7 +253,7 @@ __device__ void CUFLU_Advance( real g_Fluid_In [][5][ CUBE(FLU_NXT) ], _rho = (real)1.0 / Fluid_half[0]; vx = _rho * Fluid_half[1]; p = Hydro_Con2Pres( Fluid_half[0], Fluid_half[1], Fluid_half[2], Fluid_half[3], Fluid_half[4], Passive, - CheckMinPres_Yes, MinPres, NULL_REAL, EoS->DensEint2Pres_FuncPtr, NULL, NULL, + CheckMinPres_Yes, MinPres, PassiveFloor, NULL_REAL, EoS->DensEint2Pres_FuncPtr, NULL, NULL, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table, NULL ); # ifdef CHECK_UNPHYSICAL_IN_FLUID @@ -339,7 +341,7 @@ __device__ void CUFLU_Advance( real g_Fluid_In [][5][ CUBE(FLU_NXT) ], // apply density and internal energy floors Fluid[0] = FMAX( Fluid[0], MinDens ); Fluid[4] = Hydro_CheckMinEintInEngy( Fluid[0], Fluid[1], Fluid[2], Fluid[3], Fluid[4], - MinEint, NULL_REAL ); + MinEint, PassiveFloor, NULL_REAL ); // check negative density and energy diff --git a/src/Model_Hydro/Hydro_Aux_Check_Negative.cpp b/src/Model_Hydro/Hydro_Aux_Check_Negative.cpp index 811e9f18d2..3e62ec33ea 100644 --- a/src/Model_Hydro/Hydro_Aux_Check_Negative.cpp +++ b/src/Model_Hydro/Hydro_Aux_Check_Negative.cpp @@ -74,7 +74,7 @@ void Hydro_Aux_Check_Negative( const int lv, const int Mode, const char *comment const real Emag = NULL_REAL; # endif // MHD Pres = Hydro_Con2Pres( Fluid[DENS], Fluid[MOMX], Fluid[MOMY], Fluid[MOMZ], Fluid[ENGY], Fluid+NCOMP_FLUID, - CheckMinPres_No, NULL_REAL, Emag, + CheckMinPres_No, NULL_REAL, PassiveFloorMask, Emag, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, NULL ); # endif // DUAL_ENERGY diff --git a/src/Model_Hydro/Hydro_Init_ByFunction_AssignData.cpp b/src/Model_Hydro/Hydro_Init_ByFunction_AssignData.cpp index b871f36b4a..56f3d92fa1 100644 --- a/src/Model_Hydro/Hydro_Init_ByFunction_AssignData.cpp +++ b/src/Model_Hydro/Hydro_Init_ByFunction_AssignData.cpp @@ -365,18 +365,20 @@ void Hydro_Init_ByFunction_AssignData( const int lv ) fluid[DENS] = FMAX( fluid[DENS], (real)MIN_DENS ); # ifndef SRHD fluid[ENGY] = Hydro_CheckMinEintInEngy( fluid[DENS], fluid[MOMX], fluid[MOMY], fluid[MOMZ], fluid[ENGY], - MIN_EINT, Emag ); + MIN_EINT, PassiveFloorMask, Emag ); # endif // calculate the dual-energy variable (entropy or internal energy) # ifdef DUAL_ENERGY fluid[DUAL] = Hydro_Con2Dual( fluid[DENS], fluid[MOMX], fluid[MOMY], fluid[MOMZ], fluid[ENGY], Emag, - EoS_DensEint2Pres_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); + EoS_DensEint2Pres_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, + PassiveFloorMask ); # endif // floor and normalize passive scalars # if ( NCOMP_PASSIVE > 0 ) - for (int v=NCOMP_FLUID; v output user-defined parameters in "User/UserPara" and // Input__TestProb parameters in "Info/InputTest" // 2504 : 2025/04/29 --> output OPT__PAR_INIT_CHECK +// 2505 : 2025/05/07 --> output PassiveFloor_Var //------------------------------------------------------------------------------------------------------- void Output_DumpData_Total_HDF5( const char *FileName ) { @@ -925,7 +926,7 @@ void Output_DumpData_Total_HDF5( const char *FileName ) Emag = MHD_GetCellCenteredBEnergyInPatch( lv, PID, i, j, k, amr->MagSg[lv] ); # endif Temp = Hydro_Con2Temp( u[DENS], u[MOMX], u[MOMY], u[MOMZ], u[ENGY], u+NCOMP_FLUID, - CheckMinTemp_No, NULL_REAL, Emag, EoS_DensEint2Temp_CPUPtr, + CheckMinTemp_No, NULL_REAL, PassiveFloorMask, Emag, EoS_DensEint2Temp_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); FieldData[PID][k][j][i] = Temp; @@ -951,7 +952,7 @@ void Output_DumpData_Total_HDF5( const char *FileName ) Emag = MHD_GetCellCenteredBEnergyInPatch( lv, PID, i, j, k, amr->MagSg[lv] ); # endif Entr = Hydro_Con2Entr( u[DENS], u[MOMX], u[MOMY], u[MOMZ], u[ENGY], u+NCOMP_FLUID, - CheckMinEntr_No, NULL_REAL, Emag, EoS_DensEint2Entr_CPUPtr, + CheckMinEntr_No, NULL_REAL, PassiveFloorMask, Emag, EoS_DensEint2Entr_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); FieldData[PID][k][j][i] = Entr; } @@ -978,7 +979,7 @@ void Output_DumpData_Total_HDF5( const char *FileName ) # ifdef SRHD real Prim[NCOMP_TOTAL]; - Hydro_Con2Pri( u, Prim, (real)-HUGE_NUMBER, NULL_BOOL, NULL_INT, NULL, + Hydro_Con2Pri( u, Prim, (real)-HUGE_NUMBER, PassiveFloorMask, NULL_BOOL, NULL_INT, NULL, NULL_BOOL, NULL_REAL, EoS_DensEint2Pres_CPUPtr, EoS_DensPres2Eint_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, NULL, NULL ); @@ -986,7 +987,7 @@ void Output_DumpData_Total_HDF5( const char *FileName ) Cs2 = EoS_DensPres2CSqr_CPUPtr( Prim[0], Prim[4], NULL, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); # else Pres = Hydro_Con2Pres( u[DENS], u[MOMX], u[MOMY], u[MOMZ], u[ENGY], u+NCOMP_FLUID, - CheckMinPres_No, NULL_REAL, Emag, EoS_DensEint2Pres_CPUPtr, + CheckMinPres_No, NULL_REAL, PassiveFloorMask, Emag, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, NULL ); Cs2 = EoS_DensPres2CSqr_CPUPtr( u[DENS], Pres, u+NCOMP_FLUID, EoS_AuxArray_Flt, EoS_AuxArray_Int, @@ -1100,7 +1101,7 @@ void Output_DumpData_Total_HDF5( const char *FileName ) { for (int fv=0; fvpatch[ amr->FluSg[lv] ][lv][PID]->fluid[fv][k][j][i]; - Hydro_Con2Pri( Cons, Prim, (real)-HUGE_NUMBER, false, NULL_INT, NULL, + Hydro_Con2Pri( Cons, Prim, (real)-HUGE_NUMBER, PassiveFloorMask, false, NULL_INT, NULL, NULL_BOOL, (real)NULL_REAL, NULL, NULL, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, NULL, &LorentzFactor ); @@ -1666,7 +1667,7 @@ void FillIn_KeyInfo( KeyInfo_t &KeyInfo, const int NFieldStored ) const time_t CalTime = time( NULL ); // calendar time - KeyInfo.FormatVersion = 2504; + KeyInfo.FormatVersion = 2505; KeyInfo.Model = MODEL; KeyInfo.NLevel = NLEVEL; KeyInfo.NCompFluid = NCOMP_FLUID; @@ -2596,6 +2597,7 @@ void FillIn_InputPara( InputPara_t &InputPara, const int NFieldStored, char Fiel InputPara.Opt__FixUp_Restrict = OPT__FIXUP_RESTRICT; InputPara.FixUpRestrict_Var = FixUpVar_Restrict; InputPara.Opt__CorrAfterAllSync = OPT__CORR_AFTER_ALL_SYNC; + InputPara.PassiveFloor_Var = PassiveFloorMask; InputPara.Opt__NormalizePassive = OPT__NORMALIZE_PASSIVE; InputPara.NormalizePassive_NVar = PassiveNorm_NVar; @@ -3655,6 +3657,7 @@ void GetCompound_InputPara( hid_t &H5_TypeID, const int NFieldStored ) H5Tinsert( H5_TypeID, "Opt__FixUp_Restrict", HOFFSET(InputPara_t,Opt__FixUp_Restrict ), H5T_NATIVE_INT ); H5Tinsert( H5_TypeID, "FixUpRestrict_Var", HOFFSET(InputPara_t,FixUpRestrict_Var ), H5T_NATIVE_LONG ); H5Tinsert( H5_TypeID, "Opt__CorrAfterAllSync", HOFFSET(InputPara_t,Opt__CorrAfterAllSync ), H5T_NATIVE_INT ); + H5Tinsert( H5_TypeID, "PassiveFloor_Var", HOFFSET(InputPara_t,PassiveFloor_Var ), H5T_NATIVE_LONG ); H5Tinsert( H5_TypeID, "Opt__NormalizePassive", HOFFSET(InputPara_t,Opt__NormalizePassive ), H5T_NATIVE_INT ); H5Tinsert( H5_TypeID, "NormalizePassive_NVar", HOFFSET(InputPara_t,NormalizePassive_NVar ), H5T_NATIVE_INT ); # if ( NCOMP_PASSIVE > 0 ) diff --git a/src/Output/Output_L1Error.cpp b/src/Output/Output_L1Error.cpp index c284fea1a0..e336f2565c 100644 --- a/src/Output/Output_L1Error.cpp +++ b/src/Output/Output_L1Error.cpp @@ -404,13 +404,13 @@ void WriteFile( void (*AnalFunc_Flu)( real fluid[], const double x, const double const real Pres_Nume = Hydro_Con2Pres( Nume[DENS], Nume[MOMX], Nume[MOMY], Nume[MOMZ], Nume[ENGY], Nume+NCOMP_FLUID, - CheckMinPres_No, NULL_REAL, Emag_Nume, + CheckMinPres_No, NULL_REAL, PassiveFloorMask, Emag_Nume, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, NULL ); const real Temp_Nume = Hydro_Con2Temp( Nume[DENS], Nume[MOMX], Nume[MOMY], Nume[MOMZ], Nume[ENGY], Nume+NCOMP_FLUID, - CheckMinTemp_No, NULL_REAL, Emag_Nume, + CheckMinTemp_No, NULL_REAL, PassiveFloorMask, Emag_Nume, EoS_DensEint2Temp_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); @@ -434,12 +434,12 @@ void WriteFile( void (*AnalFunc_Flu)( real fluid[], const double x, const double # if ( MODEL == HYDRO ) const real Emag_Zero = 0.0; // Anal[ENGY] set by AnalFunc_Flu() does NOT include magnetic energy const real Pres_Anal = Hydro_Con2Pres( Anal[DENS], Anal[MOMX], Anal[MOMY], Anal[MOMZ], Anal[ENGY], - Anal+NCOMP_FLUID, CheckMinPres_No, NULL_REAL, Emag_Zero, + Anal+NCOMP_FLUID, CheckMinPres_No, NULL_REAL, PassiveFloorMask, Emag_Zero, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, NULL ); const real Temp_Anal = Hydro_Con2Temp( Anal[DENS], Anal[MOMX], Anal[MOMY], Anal[MOMZ], Anal[ENGY], - Anal+NCOMP_FLUID, CheckMinTemp_No, NULL_REAL, Emag_Zero, + Anal+NCOMP_FLUID, CheckMinTemp_No, NULL_REAL, PassiveFloorMask, Emag_Zero, EoS_DensEint2Temp_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); diff --git a/src/Output/Output_Patch.cpp b/src/Output/Output_Patch.cpp index a46dae5c65..c315d57b45 100644 --- a/src/Output/Output_Patch.cpp +++ b/src/Output/Output_Patch.cpp @@ -192,12 +192,12 @@ void Output_Patch( const int lv, const int PID, const int FluSg, const int MagSg const real Pres = ( magnetic == NULL ) ? NULL_REAL : Hydro_Con2Pres( u[DENS], u[MOMX], u[MOMY], u[MOMZ], u[ENGY], u+NCOMP_FLUID, - CheckMinPres_No, NULL_REAL, Emag, + CheckMinPres_No, NULL_REAL, PassiveFloorMask, Emag, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, NULL ); # else const real Pres = Hydro_Con2Pres( u[DENS], u[MOMX], u[MOMY], u[MOMZ], u[ENGY], u+NCOMP_FLUID, - CheckMinPres_No, NULL_REAL, NULL_REAL, + CheckMinPres_No, NULL_REAL, PassiveFloorMask, NULL_REAL, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, NULL ); # endif diff --git a/src/Output/Output_PreparedPatch_Fluid.cpp b/src/Output/Output_PreparedPatch_Fluid.cpp index 3b92b7ccba..a4c436bb05 100644 --- a/src/Output/Output_PreparedPatch_Fluid.cpp +++ b/src/Output/Output_PreparedPatch_Fluid.cpp @@ -166,7 +166,7 @@ void Output_PreparedPatch_Fluid( const int TLv, const int TPID, const real Emag = NULL_REAL; # endif fprintf( File, BlankPlusFormat_Flt, Hydro_Con2Pres(u[DENS],u[MOMX],u[MOMY],u[MOMZ],u[ENGY],u+NCOMP_FLUID, - CheckMinPres_No,NULL_REAL,Emag, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, + CheckMinPres_No,NULL_REAL,PassiveFloorMask,Emag, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, NULL) ); # endif // #if ( MODEL == HYDRO ) diff --git a/src/Refine/Flag_Real.cpp b/src/Refine/Flag_Real.cpp index 7ccd838aa0..2531fdfa0c 100644 --- a/src/Refine/Flag_Real.cpp +++ b/src/Refine/Flag_Real.cpp @@ -397,7 +397,7 @@ void Flag_Real( const int lv, const UseLBFunc_t UseLBFunc ) Pres[k][j][i] = Hydro_Con2Pres( Fluid[DENS][k][j][i], Fluid[MOMX][k][j][i], Fluid[MOMY][k][j][i], Fluid[MOMZ][k][j][i], Fluid[ENGY][k][j][i], Passive, - CheckMinPres_Yes, MIN_PRES, Emag, + CheckMinPres_Yes, MIN_PRES, PassiveFloorMask, Emag, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, NULL ); @@ -440,7 +440,7 @@ void Flag_Real( const int lv, const UseLBFunc_t UseLBFunc ) Hydro_IsUnphysical( UNPHY_MODE_CONS, Cons, NULL_REAL, EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, - ERROR_INFO, UNPHY_VERBOSE ); + PassiveFloorMask, ERROR_INFO, UNPHY_VERBOSE ); # endif HTilde = Hydro_Con2HTilde( Cons, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, diff --git a/src/Refine/Refine.cpp b/src/Refine/Refine.cpp index b34c89c0a2..a6f6ff643b 100644 --- a/src/Refine/Refine.cpp +++ b/src/Refine/Refine.cpp @@ -939,14 +939,14 @@ void Refine( const int lv, const UseLBFunc_t UseLBFunc ) Hydro_DualEnergyFix( Flu_FData[DENS][k][j][i], Flu_FData[MOMX][k][j][i], Flu_FData[MOMY][k][j][i], Flu_FData[MOMZ][k][j][i], Flu_FData[ENGY][k][j][i], Flu_FData[DUAL][k][j][i], dummy, EoS_AuxArray_Flt[1], EoS_AuxArray_Flt[2], CheckMinPres_Yes, MIN_PRES, - UseDual2FixEngy, Emag ); + PassiveFloorMask, UseDual2FixEngy, Emag ); # else // #ifdef DUAL_ENERGY // apply internal energy floor Flu_FData[ENGY][k][j][i] = Hydro_CheckMinEintInEngy( Flu_FData[DENS][k][j][i], Flu_FData[MOMX][k][j][i], Flu_FData[MOMY][k][j][i], - Flu_FData[MOMZ][k][j][i], Flu_FData[ENGY][k][j][i], MIN_EINT, Emag ); + Flu_FData[MOMZ][k][j][i], Flu_FData[ENGY][k][j][i], MIN_EINT, PassiveFloorMask, Emag ); # endif // #ifdef DUAL_ENERGY ... else ... # endif // #if ( MODEL == HYDRO ) diff --git a/src/SelfGravity/Gra_Close.cpp b/src/SelfGravity/Gra_Close.cpp index d48043f60c..a997a002a4 100644 --- a/src/SelfGravity/Gra_Close.cpp +++ b/src/SelfGravity/Gra_Close.cpp @@ -76,7 +76,8 @@ void Gra_Close( const int lv, const int SaveSg, const real h_Flu_Array_G[][GRA_N amr->patch[SaveSg][lv][PID]->fluid[MOMY][k][j][i], amr->patch[SaveSg][lv][PID]->fluid[MOMZ][k][j][i], amr->patch[SaveSg][lv][PID]->fluid[ENGY][k][j][i], - Emag, EoS_DensEint2Pres_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); + Emag, EoS_DensEint2Pres_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, + PassiveFloorMask ); # endif } # endif // #ifdef UNSPLIT_GRAVITY diff --git a/src/SourceTerms/CPU_SrcSolver.cpp b/src/SourceTerms/CPU_SrcSolver.cpp index c213c2854f..2ac2c382af 100644 --- a/src/SourceTerms/CPU_SrcSolver.cpp +++ b/src/SourceTerms/CPU_SrcSolver.cpp @@ -10,7 +10,7 @@ void CPU_SrcSolver_IterateAllCells( const double g_Corner_Array[][3], const SrcTerms_t SrcTerms, const int NPatchGroup, const real dt, const real dh, const double TimeNew, const double TimeOld, - const real MinDens, const real MinPres, const real MinEint, const EoS_t EoS ); + const real MinDens, const real MinPres, const real MinEint, const long PassiveFloor, const EoS_t EoS ); @@ -38,6 +38,7 @@ void CPU_SrcSolver_IterateAllCells( // TimeOld : Physical time before update // --> This function updates physical time from TimeOld to TimeNew // MinDens/Pres/Eint : Density, pressure, and internal energy floors +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored //------------------------------------------------------------------------------------------------------- void CPU_SrcSolver( const real h_Flu_Array_In [][FLU_NIN_S ][ CUBE(SRC_NXT) ], real h_Flu_Array_Out[][FLU_NOUT_S][ CUBE(PS1) ], @@ -45,7 +46,7 @@ void CPU_SrcSolver( const real h_Flu_Array_In [][FLU_NIN_S ][ CUBE(SRC_NXT) const double h_Corner_Array[][3], const SrcTerms_t SrcTerms, const int NPatchGroup, const real dt, const real dh, const double TimeNew, const double TimeOld, - const real MinDens, const real MinPres, const real MinEint ) + const real MinDens, const real MinPres, const real MinEint, const long PassiveFloor ) { // check @@ -67,7 +68,7 @@ void CPU_SrcSolver( const real h_Flu_Array_In [][FLU_NIN_S ][ CUBE(SRC_NXT) CPU_SrcSolver_IterateAllCells( h_Flu_Array_In, h_Flu_Array_Out, h_Mag_Array_In, h_Corner_Array, SrcTerms, NPatchGroup, dt, dh, TimeNew, TimeOld, - MinDens, MinPres, MinEint, EoS ); + MinDens, MinPres, MinEint, PassiveFloor, EoS ); } // FUNCTION : CPU_SrcSolver diff --git a/src/SourceTerms/CPU_SrcSolver_IterateAllCells.cpp b/src/SourceTerms/CPU_SrcSolver_IterateAllCells.cpp index 789d57a38e..728e1d33c7 100644 --- a/src/SourceTerms/CPU_SrcSolver_IterateAllCells.cpp +++ b/src/SourceTerms/CPU_SrcSolver_IterateAllCells.cpp @@ -48,7 +48,7 @@ void CUSRC_SrcSolver_IterateAllCells( const double g_Corner_Array[][3], const SrcTerms_t SrcTerms, const int NPatchGroup, const real dt, const real dh, const double TimeNew, const double TimeOld, - const real MinDens, const real MinPres, const real MinEint, const EoS_t EoS ) + const real MinDens, const real MinPres, const real MinEint, const long PassiveFloor, const EoS_t EoS ) #else void CPU_SrcSolver_IterateAllCells( const real g_Flu_Array_In [][FLU_NIN_S ][ CUBE(SRC_NXT) ], @@ -57,7 +57,7 @@ void CPU_SrcSolver_IterateAllCells( const double g_Corner_Array[][3], const SrcTerms_t SrcTerms, const int NPatchGroup, const real dt, const real dh, const double TimeNew, const double TimeOld, - const real MinDens, const real MinPres, const real MinEint, const EoS_t EoS ) + const real MinDens, const real MinPres, const real MinEint, const long PassiveFloor, const EoS_t EoS ) #endif { @@ -110,13 +110,13 @@ void CPU_SrcSolver_IterateAllCells( // (1) deleptonization # if ( MODEL == HYDRO ) if ( SrcTerms.Deleptonization ) - SrcTerms.Dlep_FuncPtr( fluid, B, &SrcTerms, dt, dh, x, y, z, TimeNew, TimeOld, MinDens, MinPres, MinEint, &EoS, + SrcTerms.Dlep_FuncPtr( fluid, B, &SrcTerms, dt, dh, x, y, z, TimeNew, TimeOld, MinDens, MinPres, MinEint, PassiveFloor, &EoS, SrcTerms.Dlep_AuxArrayDevPtr_Flt, SrcTerms.Dlep_AuxArrayDevPtr_Int ); # endif // (2) user-defined if ( SrcTerms.User ) - SrcTerms.User_FuncPtr( fluid, B, &SrcTerms, dt, dh, x, y, z, TimeNew, TimeOld, MinDens, MinPres, MinEint, &EoS, + SrcTerms.User_FuncPtr( fluid, B, &SrcTerms, dt, dh, x, y, z, TimeNew, TimeOld, MinDens, MinPres, MinEint, PassiveFloor, &EoS, SrcTerms.User_AuxArrayDevPtr_Flt, SrcTerms.User_AuxArrayDevPtr_Int ); // store the updated results diff --git a/src/SourceTerms/Deleptonization/CPU_Src_Deleptonization.cpp b/src/SourceTerms/Deleptonization/CPU_Src_Deleptonization.cpp index 0f6289667d..4b9637b244 100644 --- a/src/SourceTerms/Deleptonization/CPU_Src_Deleptonization.cpp +++ b/src/SourceTerms/Deleptonization/CPU_Src_Deleptonization.cpp @@ -106,6 +106,7 @@ void Src_SetAuxArray_Deleptonization( double AuxArray_Flt[], int AuxArray_Int[] // TimeOld : Physical time before update // --> This function updates physical time from TimeOld to TimeNew // MinDens/Pres/Eint : Density, pressure, and internal energy floors +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored // EoS : EoS object // AuxArray_* : Auxiliary arrays (see the Note above) // @@ -116,7 +117,7 @@ static void Src_Deleptonization( real fluid[], const real B[], const SrcTerms_t *SrcTerms, const real dt, const real dh, const double x, const double y, const double z, const double TimeNew, const double TimeOld, - const real MinDens, const real MinPres, const real MinEint, + const real MinDens, const real MinPres, const real MinEint, const long PassiveFloor, const EoS_t *EoS, const double AuxArray_Flt[], const int AuxArray_Int[] ) { diff --git a/src/SourceTerms/User_Template/CPU_Src_User_Template.cpp b/src/SourceTerms/User_Template/CPU_Src_User_Template.cpp index 561d883afb..25d579c846 100644 --- a/src/SourceTerms/User_Template/CPU_Src_User_Template.cpp +++ b/src/SourceTerms/User_Template/CPU_Src_User_Template.cpp @@ -110,6 +110,7 @@ void Src_SetAuxArray_User_Template( double AuxArray_Flt[], int AuxArray_Int[] ) // TimeOld : Physical time before update // --> This function updates physical time from TimeOld to TimeNew // MinDens/Pres/Eint : Density, pressure, and internal energy floors +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored // EoS : EoS object // AuxArray_* : Auxiliary arrays (see the Note above) // @@ -120,7 +121,7 @@ static void Src_User_Template( real fluid[], const real B[], const SrcTerms_t *SrcTerms, const real dt, const real dh, const double x, const double y, const double z, const double TimeNew, const double TimeOld, - const real MinDens, const real MinPres, const real MinEint, + const real MinDens, const real MinPres, const real MinEint, const long PassiveFloor, const EoS_t *EoS, const double AuxArray_Flt[], const int AuxArray_Int[] ) { @@ -143,7 +144,7 @@ static void Src_User_Template( real fluid[], const real B[], Emag = (real)0.0; # endif Eint = Hydro_Con2Eint( fluid[DENS], fluid[MOMX], fluid[MOMY], fluid[MOMZ], fluid[ENGY], - CheckMinEint_Yes, MinEint, Emag, EoS->GuessHTilde_CPUPtr, + CheckMinEint_Yes, MinEint, PassiveFloor, Emag, EoS->GuessHTilde_CPUPtr, EoS->HTilde2Temp_CPUPtr, EoS->AuxArray_Flt, EoS->AuxArray_Int, EoS->Table ); Enth = fluid[ENGY] - Eint; diff --git a/src/TestProblem/ELBDM/PlaneWave/Init_TestProb_ELBDM_PlaneWave.cpp b/src/TestProblem/ELBDM/PlaneWave/Init_TestProb_ELBDM_PlaneWave.cpp index 5a6adae4c8..b482cd38e4 100644 --- a/src/TestProblem/ELBDM/PlaneWave/Init_TestProb_ELBDM_PlaneWave.cpp +++ b/src/TestProblem/ELBDM/PlaneWave/Init_TestProb_ELBDM_PlaneWave.cpp @@ -357,7 +357,7 @@ void AddNewField_PlaneWave() { if ( PWave_Idx_WrappedPhase == Idx_Undefined ) - PWave_Idx_WrappedPhase = AddField( "WrappedPhase", FIXUP_FLUX_NO, FIXUP_REST_NO, NORMALIZE_NO, INTERP_FRAC_NO ); + PWave_Idx_WrappedPhase = AddField( "WrappedPhase", FIXUP_FLUX_NO, FIXUP_REST_NO, FLOOR_NO, NORMALIZE_NO, INTERP_FRAC_NO ); } // FUNCTION : AddNewField_PlaneWave diff --git a/src/TestProblem/Hydro/AGORA_IsolatedGalaxy/Init_TestProb_Hydro_AGORA_IsolatedGalaxy.cpp b/src/TestProblem/Hydro/AGORA_IsolatedGalaxy/Init_TestProb_Hydro_AGORA_IsolatedGalaxy.cpp index 068de8dc29..94b338b49e 100644 --- a/src/TestProblem/Hydro/AGORA_IsolatedGalaxy/Init_TestProb_Hydro_AGORA_IsolatedGalaxy.cpp +++ b/src/TestProblem/Hydro/AGORA_IsolatedGalaxy/Init_TestProb_Hydro_AGORA_IsolatedGalaxy.cpp @@ -505,7 +505,7 @@ void AddNewField_AGORA() // --> since Grackle may already add this field automatically when GRACKLE_METAL is enabled // --> also note that "Idx_Metal" has been predefined in Field.h if ( AGORA_UseMetal && Idx_Metal == Idx_Undefined ) - Idx_Metal = AddField( "Metal", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_NO, INTERP_FRAC_YES ); + Idx_Metal = AddField( "Metal", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, NORMALIZE_NO, INTERP_FRAC_YES ); } // FUNCTION : AddNewField_AGORA diff --git a/src/TestProblem/Hydro/Bondi/Flu_ResetByUser_Bondi.cpp b/src/TestProblem/Hydro/Bondi/Flu_ResetByUser_Bondi.cpp index c8a6a2bb32..1af81bdbb1 100644 --- a/src/TestProblem/Hydro/Bondi/Flu_ResetByUser_Bondi.cpp +++ b/src/TestProblem/Hydro/Bondi/Flu_ResetByUser_Bondi.cpp @@ -163,17 +163,19 @@ void Flu_ResetByUser_API_Bondi( const int lv, const int FluSg, const int MagSg, fluid[DENS] = FMAX( fluid[DENS], (real)MIN_DENS ); fluid[ENGY] = Hydro_CheckMinEintInEngy( fluid[DENS], fluid[MOMX], fluid[MOMY], fluid[MOMZ], fluid[ENGY], - (real)MIN_EINT, Emag ); + (real)MIN_EINT, PassiveFloorMask, Emag ); // calculate the dual-energy variable # ifdef DUAL_ENERGY fluid[DUAL] = Hydro_Con2Dual( fluid[DENS], fluid[MOMX], fluid[MOMY], fluid[MOMZ], fluid[ENGY], Emag, - EoS_DensEint2Pres_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); + EoS_DensEint2Pres_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, + PassiveFloorMask ); # endif // floor and normalize passive scalars # if ( NCOMP_PASSIVE > 0 ) - for (int v=NCOMP_FLUID; v since Grackle may already add this field automatically when GRACKLE_METAL is enabled // --> also note that "Idx_Metal" has been predefined in Field.h if ( Idx_Metal == Idx_Undefined ) - Idx_Metal = AddField( "Metal", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_NO, INTERP_FRAC_YES ); + Idx_Metal = AddField( "Metal", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, NORMALIZE_NO, INTERP_FRAC_YES ); } // FUNCTION : AddNewField_BarredPot diff --git a/src/TestProblem/Hydro/ClusterMerger/Init_TestProb_ClusterMerger.cpp b/src/TestProblem/Hydro/ClusterMerger/Init_TestProb_ClusterMerger.cpp index eff80b2fac..c4a8b6ed55 100644 --- a/src/TestProblem/Hydro/ClusterMerger/Init_TestProb_ClusterMerger.cpp +++ b/src/TestProblem/Hydro/ClusterMerger/Init_TestProb_ClusterMerger.cpp @@ -829,13 +829,13 @@ void AddNewField_ClusterMerger() { if ( Merger_Coll_UseMetals ) - Idx_Metal = AddField( "Metal", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_NO, INTERP_FRAC_NO ); + Idx_Metal = AddField( "Metal", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, NORMALIZE_NO, INTERP_FRAC_NO ); if ( ColorField1Idx == Idx_Undefined ) - ColorField1Idx = AddField( "ColorField1", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_NO, INTERP_FRAC_NO ); + ColorField1Idx = AddField( "ColorField1", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, NORMALIZE_NO, INTERP_FRAC_NO ); if ( Merger_Coll_NumHalos > 1 && ColorField2Idx == Idx_Undefined ) - ColorField2Idx = AddField( "ColorField2", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_NO, INTERP_FRAC_NO ); + ColorField2Idx = AddField( "ColorField2", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, NORMALIZE_NO, INTERP_FRAC_NO ); if ( Merger_Coll_NumHalos > 2 && ColorField3Idx == Idx_Undefined ) - ColorField3Idx = AddField( "ColorField3", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_NO, INTERP_FRAC_NO ); + ColorField3Idx = AddField( "ColorField3", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, NORMALIZE_NO, INTERP_FRAC_NO ); } // FUNCTION : AddNewField_ClusterMerger #endif diff --git a/src/TestProblem/Hydro/JetICMWall/Init_TestProb_Hydro_JetICMWall.cpp b/src/TestProblem/Hydro/JetICMWall/Init_TestProb_Hydro_JetICMWall.cpp index 7eb60d93dd..2b0d87516b 100644 --- a/src/TestProblem/Hydro/JetICMWall/Init_TestProb_Hydro_JetICMWall.cpp +++ b/src/TestProblem/Hydro/JetICMWall/Init_TestProb_Hydro_JetICMWall.cpp @@ -437,16 +437,16 @@ void AddNewField_JetICMWall() if ( JetFieldIdx == Idx_Undefined ) JetFieldIdx = AddField( "JetField", FIXUP_FLUX_YES, FIXUP_REST_YES, - NORMALIZE_YES, INTERP_FRAC_NO ); + FLOOR_YES, NORMALIZE_YES, INTERP_FRAC_NO ); if ( ICMFieldIdx == Idx_Undefined ) ICMFieldIdx = AddField( "ICMField", FIXUP_FLUX_YES, FIXUP_REST_YES, - NORMALIZE_YES, INTERP_FRAC_NO ); + FLOOR_YES, NORMALIZE_YES, INTERP_FRAC_NO ); if ( LobeFieldIdx == Idx_Undefined ) LobeFieldIdx = AddField( "LobeField", FIXUP_FLUX_YES, FIXUP_REST_YES, - NORMALIZE_YES, INTERP_FRAC_NO ); + FLOOR_YES, NORMALIZE_YES, INTERP_FRAC_NO ); if ( IntFieldIdx == Idx_Undefined ) IntFieldIdx = AddField( "IntField", FIXUP_FLUX_YES, FIXUP_REST_YES, - NORMALIZE_YES, INTERP_FRAC_NO ); + FLOOR_YES, NORMALIZE_YES, INTERP_FRAC_NO ); } // FUNCTION : AddNewField_JetICMWall #endif // #if ( MODEL == HYDRO ) diff --git a/src/TestProblem/Hydro/Plummer/Init_TestProb_Hydro_Plummer.cpp b/src/TestProblem/Hydro/Plummer/Init_TestProb_Hydro_Plummer.cpp index a7d9a49917..400eb703ef 100644 --- a/src/TestProblem/Hydro/Plummer/Init_TestProb_Hydro_Plummer.cpp +++ b/src/TestProblem/Hydro/Plummer/Init_TestProb_Hydro_Plummer.cpp @@ -529,8 +529,8 @@ void AddNewField_Plummer() if ( Plummer_AddColor ) { - Plummer_Idx_Cloud0 = AddField( "Cloud0", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_YES, INTERP_FRAC_YES ); - Plummer_Idx_Cloud1 = AddField( "Cloud1", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_YES, INTERP_FRAC_YES ); + Plummer_Idx_Cloud0 = AddField( "Cloud0", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, NORMALIZE_YES, INTERP_FRAC_YES ); + Plummer_Idx_Cloud1 = AddField( "Cloud1", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, NORMALIZE_YES, INTERP_FRAC_YES ); } } // FUNCTION : AddNewField_Plummer diff --git a/src/YT/YT_DerivedFunction.cpp b/src/YT/YT_DerivedFunction.cpp index 108eb9b276..97917ec64e 100644 --- a/src/YT/YT_DerivedFunction.cpp +++ b/src/YT/YT_DerivedFunction.cpp @@ -259,7 +259,7 @@ void Temperature_DerivedFunc( const int list_len, const long *list_gid, const ch // get temperature ((real *) data_array[lid].data_ptr)[idx] = Hydro_Con2Temp( Data[DENS][idx], Data[MOMX][idx], Data[MOMY][idx], Data[MOMZ][idx], Data[ENGY][idx], - Passive, false, NULL_REAL, Emag, + Passive, false, NULL_REAL, PassiveFloorMask, Emag, EoS_DensEint2Temp_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); }}} // i, j, k diff --git a/test_problem_deprecated/Model_ParticleOnly/TwoParOrbit/Flu_ResetByUser.cpp b/test_problem_deprecated/Model_ParticleOnly/TwoParOrbit/Flu_ResetByUser.cpp index 7689deaab2..f842d8de87 100644 --- a/test_problem_deprecated/Model_ParticleOnly/TwoParOrbit/Flu_ResetByUser.cpp +++ b/test_problem_deprecated/Model_ParticleOnly/TwoParOrbit/Flu_ResetByUser.cpp @@ -100,7 +100,8 @@ void Flu_ResetByUser( const int lv, const int FluSg, const double TTime ) // floor and normalize passive scalars # if ( NCOMP_PASSIVE > 0 ) - for (int v=NCOMP_FLUID; v