diff --git a/include/CUFLU.h b/include/CUFLU.h index 501af9b1f..2f63842c7 100644 --- a/include/CUFLU.h +++ b/include/CUFLU.h @@ -261,6 +261,14 @@ // check unphysical results in the MHM half-step prediction #if ( FLU_SCHEME == MHM ) # define MHM_CHECK_PREDICT +#ifdef MHM_CHECK_PREDICT +#ifndef SRHD +# define MHM_REPREDICT_ITER_NUM 2 +# define MHM_REPREDICT_STEPS_SAFE_FAC (real)0.4 +# define MHM_REPREDICT_SLOPE_SAFE_FAC (real)0.9 +# define MHM_REPREDICT_SUBSTEPS_MAX 4 +#endif +#endif #endif diff --git a/src/Auxiliary/Aux_Check_Parameter.cpp b/src/Auxiliary/Aux_Check_Parameter.cpp index a018831b4..5d09507d6 100644 --- a/src/Auxiliary/Aux_Check_Parameter.cpp +++ b/src/Auxiliary/Aux_Check_Parameter.cpp @@ -1107,6 +1107,31 @@ void Aux_Check_Parameter() # endif // #if ( FLU_SCHEME == MHM || FLU_SCHEME == CTU ) +// check for MHM +// ------------------------------ +# if ( FLU_SCHEME == MHM ) + +# ifdef MHM_CHECK_PREDICT +# ifndef SRHD + + if ( MHM_REPREDICT_ITER_NUM < 1 ) + Aux_Error( ERROR_INFO, "MHM_REPREDICT_ITER_NUM should be >= 1 for MHM_CHECK_PREDICT !!\n" ); + + if ( MHM_REPREDICT_STEPS_SAFE_FAC <= 0 ) + Aux_Error( ERROR_INFO, "MHM_REPREDICT_STEPS_SAFE_FAC should be > 0 for MHM_CHECK_PREDICT !!\n" ); + + if ( MHM_REPREDICT_SUBSTEPS_MAX < 1 ) + Aux_Error( ERROR_INFO, "MHM_REPREDICT_SUBSTEPS_MAX should be >= 1 for MHM_CHECK_PREDICT !!\n" ); + + if ( MHM_REPREDICT_SLOPE_SAFE_FAC <= 0 || MHM_REPREDICT_SLOPE_SAFE_FAC >= 1 ) + Aux_Error( ERROR_INFO, "MHM_REPREDICT_SLOPE_SAFE_FAC should be > 0 and < 1 for MHM_CHECK_PREDICT !!\n" ); + +# endif // #ifndef SRHD +# endif // #ifdef MHM_CHECK_PREDICT + +# endif // #if ( FLU_SCHEME == MHM ) + + // check for MHM_RP // ------------------------------ # if ( FLU_SCHEME == MHM_RP ) diff --git a/src/Auxiliary/Aux_TakeNote.cpp b/src/Auxiliary/Aux_TakeNote.cpp index ef833f01f..ae839c767 100644 --- a/src/Auxiliary/Aux_TakeNote.cpp +++ b/src/Auxiliary/Aux_TakeNote.cpp @@ -603,6 +603,12 @@ void Aux_TakeNote() # ifdef MHM_CHECK_PREDICT fprintf( Note, "MHM_CHECK_PREDICT ON\n" ); +# ifndef SRHD + fprintf( Note, "MHM_REPREDICT_ITER_NUM % d\n", MHM_REPREDICT_ITER_NUM ); + fprintf( Note, "MHM_REPREDICT_STEPS_SAFE_FAC % 14.7e\n", MHM_REPREDICT_STEPS_SAFE_FAC ); + fprintf( Note, "MHM_REPREDICT_SLOPE_SAFE_FAC % 14.7e\n", MHM_REPREDICT_SLOPE_SAFE_FAC ); + fprintf( Note, "MHM_REPREDICT_SUBSTEPS_MAX % d\n", MHM_REPREDICT_SUBSTEPS_MAX ); +# endif // #ifndef SRHD # else fprintf( Note, "MHM_CHECK_PREDICT OFF\n" ); # endif diff --git a/src/Model_Hydro/CPU_Hydro/CPU_Shared_DataReconstruction.cpp b/src/Model_Hydro/CPU_Hydro/CPU_Shared_DataReconstruction.cpp index a93fb329a..b6d7e5092 100644 --- a/src/Model_Hydro/CPU_Hydro/CPU_Shared_DataReconstruction.cpp +++ b/src/Model_Hydro/CPU_Hydro/CPU_Shared_DataReconstruction.cpp @@ -81,6 +81,34 @@ static void Hydro_HancockPredict( real fcCon[][NCOMP_LR], const real fcPri[][NCO const int NGhost, const int NEle, const real MinDens, const real MinPres, const real MinEint, const EoS_t *EoS, const long PassiveFloor ); +#ifdef MHM_CHECK_PREDICT +GPU_DEVICE +static bool Hydro_HancockPredict_CheckUnphysical( const real fcCon[][NCOMP_LR], const real dt_dh, + const EoS_t *EoS, const long PassiveFloor ); +# ifndef SRHD +GPU_DEVICE +static real Hydro_HancockPredict_GetDeplFracMax( const real fcCon[][NCOMP_LR], + const real fcCon_init[][NCOMP_LR], const real fcPri_init[][NCOMP_LR], + const real MinEint, const long PassiveFloor ); +GPU_DEVICE +static void Hydro_HancockPredict_IterateReprediction( real fcCon[][NCOMP_LR], + const real fcCon_init[][NCOMP_LR], const real dt_dh2, + const real g_cc_array[][ CUBE(FLU_NXT) ], const int cc_idx, + const real MinPres, const EoS_t *EoS, const long PassiveFloor, + const int Iter, const int IterNum, const real DeplFracMax ); +GPU_DEVICE +static void Hydro_HancockPredict_RescueByHigherSteps( real fcCon[][NCOMP_LR], + const real fcCon_init[][NCOMP_LR], const real dt_dh2, + const real MinPres, const EoS_t *EoS, const long PassiveFloor, + const int NumSubSteps ); +GPU_DEVICE +static void Hydro_HancockPredict_RescueByLowerSlopes( real fcCon[][NCOMP_LR], + const real fcCon_init[][NCOMP_LR], const real dt_dh2, + const real g_cc_array[][ CUBE(FLU_NXT) ], const int cc_idx, + const real MinPres, const EoS_t *EoS, const long PassiveFloor, + const real ReducedSlopeRatio ); +#endif // #ifndef SRHD +#endif // #ifdef MHM_CHECK_PREDICT #ifdef MHD GPU_DEVICE void Hydro_ConFC2PriCC_MHM( real g_PriVar[][ CUBE(FLU_NXT) ], @@ -2071,6 +2099,14 @@ void Hydro_HancockPredict( real fcCon[][NCOMP_LR], const real fcPri[][NCOMP_LR], const EoS_t *EoS, const long PassiveFloor ) { +# ifdef MHM_CHECK_PREDICT +// back up the input values in case they need to be recalculated + real fcCon_init[6][NCOMP_LR]; + for (int f=0; f<6; f++) + for (int v=0; v zero slope +// repredict_iter_num = 3: more steps -> lower slope -> zero slope +// repredict_iter_num = 4: more steps -> lower slope -> even lower slope -> zero slope +// repredict_iter_num = 5: more steps -> even more steps -> lower slope -> even lower slope -> zero slope + const int repredict_iter_num = MHM_REPREDICT_ITER_NUM; // maximum number of iterations of reprediction +# endif // #ifdef SRHD ... else ... + for (int repredict_iter=0; repredict_iterDensEint2Pres_FuncPtr, - EoS->GuessHTilde_FuncPtr, EoS->HTilde2Temp_FuncPtr, - EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table, - PassiveFloor, ERROR_INFO, UNPHY_SILENCE ) ) - reset_cell = true; +// check whether there are unphysical results after the update + repredict = Hydro_HancockPredict_CheckUnphysical( fcCon, dt_dh2*2, EoS, PassiveFloor ); -# else - - if ( Hydro_IsUnphysical_Single( fcCon[f][DENS], "density", TINY_NUMBER, HUGE_NUMBER, ERROR_INFO, UNPHY_SILENCE ) ) - reset_cell = true; - -# ifndef BAROTROPIC_EOS - if ( Hydro_IsUnphysical_Single( fcCon[f][4], "energy", TINY_NUMBER, HUGE_NUMBER, ERROR_INFO, UNPHY_SILENCE ) ) - reset_cell = true; - - if ( Hydro_IsUnphysical_Single( fcPri[f][4], "pressure", (real)0.0, HUGE_NUMBER, ERROR_INFO, UNPHY_SILENCE ) ) - reset_cell = true; -# endif // #ifndef BAROTROPIC_EOS -# endif // #ifdef SRHD ... else ... +// break immediately when there is no need to repredict + if ( !repredict ) break; -// set to the cell-centered values before update - if ( reset_cell ) +# ifndef SRHD +// repredict with more accurate or safer methods to avoid overshoot + if ( repredict_iter != repredict_iter_num-1 ) + { +// decide the maxium depletion fraction of density, energy, and internal energy +// from the initial states and the original updates + if ( repredict_iter == 0 ) + depl_frac_max = Hydro_HancockPredict_GetDeplFracMax( fcCon, fcCon_init, fcPri, MinEint, PassiveFloor ); + +// invoke the rescue methods to repredict + Hydro_HancockPredict_IterateReprediction( fcCon, fcCon_init, dt_dh2, + g_cc_array, cc_idx, MinPres, EoS, PassiveFloor, + repredict_iter, repredict_iter_num, depl_frac_max ); + +// reset the flag for the next check + repredict = false; + } +// repredict for the last time + else // if ( repredict_iter != repredict_iter_num-1 ) +# endif // #ifndef SRHD { - for (int face=0; face<6; face++) +// fall back to the zero-slope cell-centered values as the last resort + for (int f=0; f<6; f++) for (int v=0; vDensEint2Pres_FuncPtr, + EoS->GuessHTilde_FuncPtr, EoS->HTilde2Temp_FuncPtr, + EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table, + PassiveFloor, ERROR_INFO, UNPHY_SILENCE ) ) + { + isUnphysical = true; + break; + } + +# else +// 1. check the value of density + if ( Hydro_IsUnphysical_Single( fcCon[f][DENS], "density", TINY_NUMBER, HUGE_NUMBER, ERROR_INFO, UNPHY_SILENCE ) ) + { + isUnphysical = true; + break; + } + +// 2. check the value of velocity + if ( MAX( FABS(fcCon[f][MOMX]), MAX( FABS(fcCon[f][MOMY]), FABS(fcCon[f][MOMZ]) ) )*dt_dh > FABS(fcCon[f][DENS]) ) + { + isUnphysical = true; + break; + } + +# ifndef BAROTROPIC_EOS +// 3. check the value of total energy + if ( Hydro_IsUnphysical_Single( fcCon[f][4], "energy", TINY_NUMBER, HUGE_NUMBER, ERROR_INFO, UNPHY_SILENCE ) ) + { + isUnphysical = true; + break; + } + +// 4. check the value of pressure +# ifdef MHD + const real Emag = (real)0.5*( SQR(fcCon[f][MAG_OFFSET+0]) + SQR(fcCon[f][MAG_OFFSET+1]) + SQR(fcCon[f][MAG_OFFSET+2]) ); +# else + const real Emag = NULL_REAL; +# endif // MHD + const bool CheckMinPres_No = false; + const real Pres = Hydro_Con2Pres( fcCon[f][DENS], fcCon[f][MOMX], fcCon[f][MOMY], fcCon[f][MOMZ], fcCon[f][ENGY], + fcCon[f]+NCOMP_FLUID, CheckMinPres_No, NULL_REAL, PassiveFloor, Emag, + EoS->DensEint2Pres_FuncPtr, NULL, NULL, + EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table, NULL ); + + if ( Hydro_IsUnphysical_Single( Pres, "pressure", (real)0.0, HUGE_NUMBER, ERROR_INFO, UNPHY_SILENCE ) ) + { + isUnphysical = true; + break; + } +# endif // #ifndef BAROTROPIC_EOS +# endif // #ifdef SRHD ... else ... + } // for (int f=0; f<6; f++) + + return isUnphysical; + +} // FUNCTION : Hydro_HancockPredict_CheckUnphysical + + + +# ifndef SRHD +//------------------------------------------------------------------------------------------------------- +// Function : Hydro_HancockPredict_GetDeplFracMax +// Description : Get the maximum of depletion fraction of the face-centered variables after update among all faces +// +// Note : 1. Work for the MHM scheme +// 2. Invoked by Hydro_HancockPredict() only when unphysical result is found +// 3. The value will be used to estimate the required number of steps and the reduced slope +// to rescue the unphysical results +// +// Parameter : fcCon : Updated face-centered conserved variables +// fcCon_init : Input face-centered conserved variables before update +// fcPri_init : Input face-centered primitive variables before update +// MinEint : Internal energy floors +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored +// +// Return : depl_frac_max : Maximum depletion fraction +//------------------------------------------------------------------------------------------------------- +GPU_DEVICE +real Hydro_HancockPredict_GetDeplFracMax( const real fcCon[][NCOMP_LR], + const real fcCon_init[][NCOMP_LR], const real fcPri_init[][NCOMP_LR], + const real MinEint, const long PassiveFloor ) +{ + +// initialize it with a small but non-zero value + real depl_frac_max = MAX_ERROR; + +// depletion fraction is defined as (var_initial - var_updated)/var_initial +// positive (>0) when the variable decreased after update +// =1 leads to vacuum, >1 is unphysical, and the larger the dangerous +// negative (<0) when the variable increased after update, and it is safe + +// loop through all the faces to find the maximum + for (int f=0; f<6; f++) + { +// density depletion fraction + depl_frac_max = MAX( (real)1.0-fcCon[f][DENS]/fcCon_init[f][DENS], depl_frac_max ); + +// energy depletion fraction + depl_frac_max = MAX( (real)1.0-fcCon[f][ENGY]/fcCon_init[f][ENGY], depl_frac_max ); + +// internal energy depletion fraction +# ifndef BAROTROPIC_EOS +// compute the magnetic energy +# ifdef MHD + const real Emag = (real)0.5*( SQR( fcCon[f][MAG_OFFSET+0]) + SQR( fcCon[f][MAG_OFFSET+1]) + SQR( fcCon[f][MAG_OFFSET+2]) ); + const real initlEmag = (real)0.5*( SQR(fcCon_init[f][MAG_OFFSET+0]) + SQR(fcCon_init[f][MAG_OFFSET+1]) + SQR(fcCon_init[f][MAG_OFFSET+2]) ); +# else + const real initlEmag = NULL_REAL; +# endif // MHD + +// compute the initial internal energy; check minimum in case it is negative in the input + const bool CheckMinEint_Yes = true; + const real initlEint = Hydro_Con2Eint( fcCon_init[f][DENS], fcCon_init[f][MOMX], + fcCon_init[f][MOMY], fcCon_init[f][MOMZ], fcCon_init[f][ENGY], + CheckMinEint_Yes, MinEint, PassiveFloor, initlEmag, + NULL, NULL, NULL, NULL, NULL ); + +// estimate the decreased amount of internal energy from the initial change rate with linearization +// --> Eint = Engy - 0.5*(Mom_i**2)/Dens +// --> dEint = dEngy - Vel_i*dMom_i + 0.5*(Vel_i**2)*dDens + const real depldEint = ( (fcCon_init[f][ENGY] - fcCon[f][ENGY]) +# ifdef MHD + - ( initlEmag - Emag ) +# endif // MHD + - (fcCon_init[f][MOMX] - fcCon[f][MOMX])*fcPri_init[f][1] + - (fcCon_init[f][MOMY] - fcCon[f][MOMY])*fcPri_init[f][2] + - (fcCon_init[f][MOMZ] - fcCon[f][MOMZ])*fcPri_init[f][3] + + (fcCon_init[f][DENS] - fcCon[f][DENS])*(real)0.5*( SQR(fcPri_init[f][1])+SQR(fcPri_init[f][2])+SQR(fcPri_init[f][3]) ) ); + + depl_frac_max = MAX( depldEint/initlEint, depl_frac_max ); +# endif // #ifndef BAROTROPIC_EOS + } // for (int f=0; f<6; f++) + + return depl_frac_max; + +} // FUNCTION : Hydro_HancockPredict_GetDeplFracMax + + + +//------------------------------------------------------------------------------------------------------- +// Function : Hydro_HancockPredict_IterateReprediction +// Description : Redo the HancockPredict evolution by iterations of different trials +// +// Note : 1. Work for the MHM scheme +// 2. Invoked by Hydro_HancockPredict() only when unphysical result is found +// 3. Invoke different rescue methods accordingly to the iteration times +// 4. Input variables must be conserved variables +// +// Parameter : fcCon : Face-centered conserved variables to be updated +// fcCon_init : Input face-centered conserved variables before update +// dt_dh2 : 0.5 * Time interval to advance solution / Cell size +// g_cc_array : Array storing the cell-centered conserved variables for checking +// negative density and pressure +// --> It is just the input array Flu_Array_In[] +// cc_idx : Index for accessing g_cc_array[] +// MinPres : Pressure floor +// EoS : EoS object +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored +// Iter : Current iteration +// IterNum : Total number of iterations +// DeplFracMax : Maximum depletion fraction based on the initial prediction +//------------------------------------------------------------------------------------------------------- +GPU_DEVICE +void Hydro_HancockPredict_IterateReprediction( real fcCon[][NCOMP_LR], + const real fcCon_init[][NCOMP_LR], const real dt_dh2, + const real g_cc_array[][ CUBE(FLU_NXT) ], const int cc_idx, + const real MinPres, const EoS_t *EoS, const long PassiveFloor, + const int Iter, const int IterNum, const real DeplFracMax ) +{ + +// iteration division of rescuing methods +// Iter = [ 0, 1, ..., IterHalf-1] -> more accurate method +// Iter = [IterHalf, IterHalf+1, ..., IterNum -2] -> more diffusive method + const int IterHalf = (IterNum-1)/2; + +// for the first half, try using higher-order integration method and smaller time-step + if ( Iter < IterHalf ) + { +// estimate the required numer of sub-steps according to the maximum depletion fraction + const real safety_factor = MHM_REPREDICT_STEPS_SAFE_FAC/(1< 0.2 -> 0.1 + const int num_substeps_est = (int)ceilf( DeplFracMax/safety_factor ); +// or at least increase by iterations + const int num_substeps_min = 1< 2 -> 4 + const int num_substeps_max = MHM_REPREDICT_SUBSTEPS_MAX; // to avoid too many sub-steps, for the performance consideration + const int num_substeps = MIN( MAX( num_substeps_est, num_substeps_min ), num_substeps_max ); + + Hydro_HancockPredict_RescueByHigherSteps( fcCon, fcCon_init, dt_dh2, + MinPres, EoS, PassiveFloor, num_substeps ); + } +// for the second half, try reducing the slope for data reconstruction + else // if ( Iter < IterHalf ) + { +// reduce the slope according to the maximum depletion fraction + const real safety_factor = MHM_REPREDICT_SLOPE_SAFE_FAC/(1<<(Iter-IterHalf)); // e.g., 0.9 -> 0.45 -> 0.225 + const real reduced_slope_ratio = MIN( safety_factor/DeplFracMax, safety_factor ); + + Hydro_HancockPredict_RescueByLowerSlopes( fcCon, fcCon_init, dt_dh2, g_cc_array, cc_idx, + MinPres, EoS, PassiveFloor, reduced_slope_ratio ); + } // if ( Iter < IterHalf ) ... else ... + +} + + + +//------------------------------------------------------------------------------------------------------- +// Function : Hydro_HancockPredict_RescueByHigherSteps +// Description : Evolve the face-centered variables by half time-step by subcyling steps +// +// Note : 1. Work for the MHM scheme +// 2. Invoked by Hydro_HancockPredict_IterateReprediction() when unphysical result is found in Hydro_HancockPredict() +// 3. Input variables must be conserved variables +// 4. Warning: This is not helpful for rescuing most of the time and can waste time +// +// Parameter : fcCon : Face-centered conserved variables to be updated +// fcCon_init : Input face-centered conserved variables before update +// dt_dh2 : 0.5 * Time interval to advance solution / Cell size +// MinPres : Pressure floors +// EoS : EoS object +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored +// NumSubSteps : Number of sub-steps to evolve +//------------------------------------------------------------------------------------------------------- +GPU_DEVICE +void Hydro_HancockPredict_RescueByHigherSteps( real fcCon[][NCOMP_LR], + const real fcCon_init[][NCOMP_LR], const real dt_dh2, + const real MinPres, const EoS_t *EoS, const long PassiveFloor, + const int NumSubSteps ) +{ + +// restore the initial face-centered variables + for (int f=0; f<6; f++) + for (int v=0; vDensEint2Pres_FuncPtr, + EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table, NULL ); + } // for (int f=0; f<6; f++) + +// sub-1-2. update the face-centered variables for half sub-step by the initial flux + real fcCon_sub_half[6][NCOMP_LR]; + for (int f=0; f<6; f++) + for (int v=0; vDensEint2Pres_FuncPtr, + EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table, NULL ); + } // for (int f=0; f<6; f++) + +// sub-2-2. update the face-centered variables for full sub-step by the half sub-step flux + for (int v=0; v It is just the input array Flu_Array_In[] +// cc_idx : Index for accessing g_cc_array[] +// MinPres : Density, pressure, and internal energy floors +// EoS : EoS object +// PassiveFloor : Bitwise flag to specify the passive scalars to be floored +// ReducedSlopeRatio : Reduced slope / original slope +//------------------------------------------------------------------------------------------------------- +GPU_DEVICE +void Hydro_HancockPredict_RescueByLowerSlopes( real fcCon[][NCOMP_LR], + const real fcCon_init[][NCOMP_LR], const real dt_dh2, + const real g_cc_array[][ CUBE(FLU_NXT) ], const int cc_idx, + const real MinPres, const EoS_t *EoS, const long PassiveFloor, + const real ReducedSlopeRatio ) +{ + +// reconstruct the face-centered variables with the reduced slope + for (int f=0; f<6; f++) + for (int v=0; vDensEint2Pres_FuncPtr, + EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table, NULL ); + } // for (int f=0; f<6; f++) + +// update the face-centered variables + for (int v=0; v