diff --git a/doc/wiki/Runtime-Parameters-related/Runtime-Parameters:-All.md b/doc/wiki/Runtime-Parameters-related/Runtime-Parameters:-All.md index 9129bf366..c5811b6d3 100644 --- a/doc/wiki/Runtime-Parameters-related/Runtime-Parameters:-All.md +++ b/doc/wiki/Runtime-Parameters-related/Runtime-Parameters:-All.md @@ -79,7 +79,9 @@ For variables with `Default/Min/Max` labeled as `Depend`, click the parameter na | ELBDM_MASS | -1.0 | 2.22507386e-308 | None | particle mass in ev/c^2 (input unit is fixed even when OPT__UNIT or COMOVING is on) | | ELBDM_MATCH_PHASE | 1 | None | None | match child phases with father phases during data restriction [1] ##ELBDM_HYBRID ONLY## | | ELBDM_PLANCK_CONST | -1.0 | None | None | reduced Planck constant (will be overwritten if OPT__UNIT or COMOVING is on) | -| ELBDM_REMOVE_MOTION_CM | ELBDM_REMOVE_MOTION_CM_NONE | 0 | 2 | remove the motion of center-of-mass (must enable OPT__CK_CONSERVATION): (0=off, 1=init, 2=every step) [0] | +| ELBDM_REMOVE_MOTION_CM | ELBDM_REMOVE_MOTION_CM_NONE | 0 | 2 | remove the motion of center-of-mass (must enable OPT__CK_CONSERVATION): (0=off, 1=init, 2=every step) [0] | +| ELBDM_RESCALE_MASS_ERROR | 0 | None | None | Rescale the total ELBDM mass to initial total ELBDM mass [0] | +| ELBDM_RESCALE_MASS_STEPS | 100 | 1 | None | Per number of steps to rescale the total ELBDM mass [100] | | ELBDM_TAYLOR3_AUTO | 0 | None | None | Optimize ELBDM_TAYLOR3_COEFF automatically to minimize the damping at kmax [0] | | ELBDM_TAYLOR3_COEFF | 1.0/6.0 | None | None | 3rd Taylor expansion coefficient [1.0/6.0] ##USELESS if ELBDM_TAYLOR3_AUTO is on## | | [[ END_STEP \| Runtime-Parameters:-General#END_STEP ]] | -1L | None | None | end step (<0=auto -> must be set by test problems or restart) [-1] | diff --git a/example/test_problem/Template/Input__Parameter b/example/test_problem/Template/Input__Parameter index d7742ab30..0ad70ab1e 100644 --- a/example/test_problem/Template/Input__Parameter +++ b/example/test_problem/Template/Input__Parameter @@ -268,6 +268,8 @@ ELBDM_TAYLOR3_COEFF 0.166666667 # 3rd Taylor expansion coefficient [1. ELBDM_TAYLOR3_AUTO 0 # Optimize ELBDM_TAYLOR3_COEFF automatically to minimize the damping at kmax [0] ELBDM_REMOVE_MOTION_CM 0 # remove the motion of center-of-mass (must enable OPT__CK_CONSERVATION): # (0=off, 1=init, 2=every step) [0] +ELBDM_RESCALE_MASS_ERROR 0 # rescale the total ELBDM mass to initial total ELBDM mass (must enable OPT__CK_CONSERVATION) [0] +ELBDM_RESCALE_MASS_STEPS 100 # Per number of steps to rescale the total ELBDM mass (must enable OPT__CK_CONSERVATION) [100] ELBDM_BASE_SPECTRAL 0 # adopt the spectral method to evolve base-level wave function (must enable SUPPORT_FFTW) [0] ELBDM_MATCH_PHASE 1 # match child phases with father phases during data restriction [1] ##ELBDM_HYBRID ONLY## ELBDM_FIRST_WAVE_LEVEL -1 # level at which to switch to the wave solver (must >=1) [-1] ##ELBDM_HYBRID ONLY## diff --git a/include/Global.h b/include/Global.h index aabea2f5b..846713754 100644 --- a/include/Global.h +++ b/include/Global.h @@ -167,6 +167,8 @@ extern int OPT__FLAG_SPECTRAL_N; extern double FlagTable_Spectral[NLEVEL-1][2]; extern ELBDMRemoveMotionCM_t ELBDM_REMOVE_MOTION_CM; +extern bool ELBDM_RESCALE_MASS_ERROR; +extern int ELBDM_RESCALE_MASS_STEPS; extern bool ELBDM_BASE_SPECTRAL; #else diff --git a/include/HDF5_Typedef.h b/include/HDF5_Typedef.h index edde083e3..3de936ea7 100644 --- a/include/HDF5_Typedef.h +++ b/include/HDF5_Typedef.h @@ -606,6 +606,8 @@ struct InputPara_t double ELBDM_Taylor3_Coeff; int ELBDM_Taylor3_Auto; int ELBDM_RemoveMotionCM; + int ELBDM_RescaleMassError; + int ELBDM_RescaleMassSteps; int ELBDM_BaseSpectral; # if ( ELBDM_SCHEME == ELBDM_HYBRID ) int ELBDM_FirstWaveLevel; diff --git a/include/Prototype.h b/include/Prototype.h index 767fd31d0..e034adebf 100644 --- a/include/Prototype.h +++ b/include/Prototype.h @@ -608,6 +608,7 @@ bool ELBDM_Flag_Interference( const int i, const int j, const int k, const rea real ELBDM_UnwrapPhase( const real Phase_Ref, const real Phase_Wrapped ); real ELBDM_SetTaylor3Coeff( const real dt, const real dh, const real Eta ); void ELBDM_RemoveMotionCM(); +void ELBDM_RescaleMassError(); #ifdef SUPPORT_FFTW void CPU_ELBDMSolver_FFT( const real dt, const double PrepTime, const int SaveSg ); #endif diff --git a/src/Auxiliary/Aux_Check_Conservation.cpp b/src/Auxiliary/Aux_Check_Conservation.cpp index 850090e7b..3979ba492 100644 --- a/src/Auxiliary/Aux_Check_Conservation.cpp +++ b/src/Auxiliary/Aux_Check_Conservation.cpp @@ -5,6 +5,7 @@ // --> declared in "Model_ELBDM/ELBDM_RemoveMotionCM.cpp" #if ( MODEL == ELBDM ) extern double ELBDM_Vcm[3]; +extern double ELBDM_MassPsi; #endif @@ -784,4 +785,17 @@ void Aux_Check_Conservation( const char *comment ) } # endif +// calculate the ELBDM total mass for ELBDM_RescaleMassError.cpp +# if ( MODEL == ELBDM ) + if ( ELBDM_RESCALE_MASS_ERROR == true ) + { + if ( MPI_Rank == 0 ) + { + ELBDM_MassPsi = Fluid_AllRank[0]; + } +// broadcast + MPI_Bcast( &ELBDM_MassPsi, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD ); + } +# endif + } // FUNCTION : Aux_Check_Conservation diff --git a/src/Auxiliary/Aux_Check_Parameter.cpp b/src/Auxiliary/Aux_Check_Parameter.cpp index a018831b4..10b209cc4 100644 --- a/src/Auxiliary/Aux_Check_Parameter.cpp +++ b/src/Auxiliary/Aux_Check_Parameter.cpp @@ -1272,6 +1272,9 @@ void Aux_Check_Parameter() Aux_Error( ERROR_INFO, "\"%s\" does NOT support \"%s\" !!\n", "ELBDM_REMOVE_MOTION_CM", "BITWISE_REPRODUCIBILITY" ); # endif + if ( ELBDM_RESCALE_MASS_ERROR == true && !OPT__CK_CONSERVATION ) + Aux_Error( ERROR_INFO, "\"%s\" must work with \"%s\" !!\n", "ELBDM_RESCALE_MASS_ERROR", "OPT__CK_CONSERVATION" ); + for (int f=0; f<6; f++) if ( ELBDM_BASE_SPECTRAL && OPT__BC_FLU[f] != BC_FLU_PERIODIC ) Aux_Error( ERROR_INFO, "ELBDM_BASE_SPECTRAL only works with periodic boundary condition (OPT__BC_FLU=1) !!\n" ); diff --git a/src/Auxiliary/Aux_TakeNote.cpp b/src/Auxiliary/Aux_TakeNote.cpp index ef833f01f..eb18cbb3b 100644 --- a/src/Auxiliary/Aux_TakeNote.cpp +++ b/src/Auxiliary/Aux_TakeNote.cpp @@ -1247,20 +1247,22 @@ void Aux_TakeNote() # endif } else { - fprintf( Note, "ELBDM_MASS % 14.7e\n", ELBDM_MASS ); + fprintf( Note, "ELBDM_MASS % 14.7e\n", ELBDM_MASS ); } - fprintf( Note, "ELBDM_PLANCK_CONST % 14.7e\n", ELBDM_PLANCK_CONST ); - fprintf( Note, "ELBDM_ETA % 14.7e\n", ELBDM_ETA ); + fprintf( Note, "ELBDM_PLANCK_CONST % 14.7e\n", ELBDM_PLANCK_CONST ); + fprintf( Note, "ELBDM_ETA % 14.7e\n", ELBDM_ETA ); # ifdef QUARTIC_SELF_INTERACTION - fprintf( Note, "ELBDM_LAMBDA % 14.7e\n", ELBDM_LAMBDA ); -# endif - fprintf( Note, "ELBDM_TAYLOR3_COEFF % 14.7e\n", ELBDM_TAYLOR3_COEFF ); - fprintf( Note, "ELBDM_TAYLOR3_AUTO % d\n", ELBDM_TAYLOR3_AUTO ); - fprintf( Note, "ELBDM_REMOVE_MOTION_CM % d\n", ELBDM_REMOVE_MOTION_CM ); - fprintf( Note, "ELBDM_BASE_SPECTRAL % d\n", ELBDM_BASE_SPECTRAL ); + fprintf( Note, "ELBDM_LAMBDA % 14.7e\n", ELBDM_LAMBDA ); +# endif + fprintf( Note, "ELBDM_TAYLOR3_COEFF % 14.7e\n", ELBDM_TAYLOR3_COEFF ); + fprintf( Note, "ELBDM_TAYLOR3_AUTO % d\n", ELBDM_TAYLOR3_AUTO ); + fprintf( Note, "ELBDM_REMOVE_MOTION_CM % d\n", ELBDM_REMOVE_MOTION_CM ); + fprintf( Note, "ELBDM_RESCALE_MASS_ERROR % d\n", ELBDM_RESCALE_MASS_ERROR); + fprintf( Note, "ELBDM_RESCALE_MASS_STEPS % d\n", ELBDM_RESCALE_MASS_STEPS); + fprintf( Note, "ELBDM_BASE_SPECTRAL % d\n", ELBDM_BASE_SPECTRAL ); # if ( ELBDM_SCHEME == ELBDM_HYBRID ) - fprintf( Note, "ELBDM_MATCH_PHASE % d\n", ELBDM_MATCH_PHASE ); - fprintf( Note, "ELBDM_FIRST_WAVE_LEVEL % d\n", ELBDM_FIRST_WAVE_LEVEL ); + fprintf( Note, "ELBDM_MATCH_PHASE % d\n", ELBDM_MATCH_PHASE ); + fprintf( Note, "ELBDM_FIRST_WAVE_LEVEL % d\n", ELBDM_FIRST_WAVE_LEVEL ); # endif # else diff --git a/src/Init/Init_ByRestart_HDF5.cpp b/src/Init/Init_ByRestart_HDF5.cpp index e08c16a93..2507df867 100644 --- a/src/Init/Init_ByRestart_HDF5.cpp +++ b/src/Init/Init_ByRestart_HDF5.cpp @@ -2134,6 +2134,8 @@ void Check_InputPara( const char *FileName, const int FormatVersion ) LoadField( "ELBDM_Taylor3_Coeff", &RS.ELBDM_Taylor3_Coeff, SID, TID, NonFatal, &RT.ELBDM_Taylor3_Coeff, 1, NonFatal ); LoadField( "ELBDM_Taylor3_Auto", &RS.ELBDM_Taylor3_Auto, SID, TID, NonFatal, &RT.ELBDM_Taylor3_Auto, 1, NonFatal ); LoadField( "ELBDM_RemoveMotionCM", &RS.ELBDM_RemoveMotionCM, SID, TID, NonFatal, &RT.ELBDM_RemoveMotionCM, 1, NonFatal ); + LoadField( "ELBDM_RescaleMassError", &RS.ELBDM_RescaleMassError, SID, TID, NonFatal, &RT.ELBDM_RescaleMassError, 1, NonFatal ); + LoadField( "ELBDM_RescaleMassSteps", &RS.ELBDM_RescaleMassSteps, SID, TID, NonFatal, &RT.ELBDM_RescaleMassSteps, 1, NonFatal ); LoadField( "ELBDM_BaseSpectral", &RS.ELBDM_BaseSpectral, SID, TID, NonFatal, &RT.ELBDM_BaseSpectral, 1, NonFatal ); # if ( ELBDM_SCHEME == ELBDM_HYBRID ) // ELBDM_FIRST_WAVE_LEVEL currently cannot be changed upon restart because the code cannot robustly handle the conversion diff --git a/src/Init/Init_Load_Parameter.cpp b/src/Init/Init_Load_Parameter.cpp index 7149d02b3..7a565a4b0 100644 --- a/src/Init/Init_Load_Parameter.cpp +++ b/src/Init/Init_Load_Parameter.cpp @@ -336,6 +336,8 @@ void Init_Load_Parameter() ReadPara->Add( "ELBDM_TAYLOR3_COEFF", &ELBDM_TAYLOR3_COEFF, 1.0/6.0, NoMin_double, NoMax_double ); ReadPara->Add( "ELBDM_TAYLOR3_AUTO", &ELBDM_TAYLOR3_AUTO, false, Useless_bool, Useless_bool ); ReadPara->Add( "ELBDM_REMOVE_MOTION_CM", &ELBDM_REMOVE_MOTION_CM, ELBDM_REMOVE_MOTION_CM_NONE, 0, 2 ); + ReadPara->Add( "ELBDM_RESCALE_MASS_ERROR", &ELBDM_RESCALE_MASS_ERROR, false, Useless_bool, Useless_bool ); + ReadPara->Add( "ELBDM_RESCALE_MASS_STEPS", &ELBDM_RESCALE_MASS_STEPS, 100, 1, NoMax_int ); ReadPara->Add( "ELBDM_BASE_SPECTRAL", &ELBDM_BASE_SPECTRAL, false, Useless_bool, Useless_bool ); # if ( ELBDM_SCHEME == ELBDM_HYBRID ) ReadPara->Add( "ELBDM_MATCH_PHASE", &ELBDM_MATCH_PHASE, true, Useless_bool, Useless_bool ); diff --git a/src/Main/Main.cpp b/src/Main/Main.cpp index 70358820d..620b05c70 100644 --- a/src/Main/Main.cpp +++ b/src/Main/Main.cpp @@ -159,6 +159,8 @@ double DT__HYBRID_CFL, DT__HYBRID_CFL_INIT, DT__HYBRID_VELOCITY, D double ELBDM_LAMBDA; #endif ELBDMRemoveMotionCM_t ELBDM_REMOVE_MOTION_CM; +bool ELBDM_RESCALE_MASS_ERROR; +int ELBDM_RESCALE_MASS_STEPS; bool ELBDM_BASE_SPECTRAL; #else @@ -746,6 +748,9 @@ int main( int argc, char *argv[] ) if ( ELBDM_REMOVE_MOTION_CM == ELBDM_REMOVE_MOTION_CM_EVERY_STEP ) TIMING_FUNC( ELBDM_RemoveMotionCM(), Timer_Main[4], TIMER_ON ); + + if ( ELBDM_RESCALE_MASS_ERROR && Step % ELBDM_RESCALE_MASS_STEPS == 0 ) + TIMING_FUNC( ELBDM_RescaleMassError(), Timer_Main[4], TIMER_ON ); # endif // #if ( MODEL == ELBDM ) // --------------------------------------------------------------------------------------------------- diff --git a/src/Makefile_base b/src/Makefile_base index 05a5ff17b..25d2f356e 100644 --- a/src/Makefile_base +++ b/src/Makefile_base @@ -174,7 +174,7 @@ CPU_FILE += CPU_ELBDMSolver_FD.cpp CPU_ELBDMSolver_FFT.cpp CPU_ELBDMSolver_ CPU_ELBDMSolver_HJ.cpp ELBDM_Init_ByFunction_AssignData.cpp ELBDM_GetTimeStep_Fluid.cpp ELBDM_GetTimeStep_Hybrid_CFL.cpp \ ELBDM_Flag_EngyDensity.cpp ELBDM_Flag_Interference.cpp ELBDM_Flag_Spectral.cpp ELBDM_UnwrapPhase.cpp \ ELBDM_GetTimeStep_Phase.cpp ELBDM_GetTimeStep_Hybrid_Velocity.cpp ELBDM_HasWaveCounterpart.cpp ELBDM_SetTaylor3Coeff.cpp \ - ELBDM_GramFE_EvolutionMatrix.cpp ELBDM_RemoveMotionCM.cpp ELBDM_Aux_Record_Hybrid.cpp + ELBDM_GramFE_EvolutionMatrix.cpp ELBDM_RemoveMotionCM.cpp ELBDM_RescaleMassError.cpp ELBDM_Aux_Record_Hybrid.cpp vpath %.cu Model_ELBDM/GPU_ELBDM vpath %.cpp Model_ELBDM/CPU_ELBDM Model_ELBDM diff --git a/src/Model_ELBDM/ELBDM_RescaleMassError.cpp b/src/Model_ELBDM/ELBDM_RescaleMassError.cpp new file mode 100644 index 000000000..92634228f --- /dev/null +++ b/src/Model_ELBDM/ELBDM_RescaleMassError.cpp @@ -0,0 +1,90 @@ +#include "GAMER.h" + +#if ( MODEL == ELBDM ) + +// global variable to store the ELBDM total mass + double ELBDM_MassPsi = NULL_REAL; +static double ELBDM_InitMassPsi = NULL_REAL; + +//------------------------------------------------------------------------------------------------------- +// Function : ELBDM_RescaleMassError +// Description : Remove the mass error created bo floating numerical error, only for base level. +// +// Note : 1. Work with the option ELBDM_RESCALE_MASS_ERROR +// --> Must also enable OPT__CK_CONSERVATION since it relies on Aux_Check_Conservation() +// to calculate the total ELBDM mass (ELBDM_MassPsi) +// 2. Invoked by Main() +// +// Parameter : None +// +// Return : amr->fluid[DENS/REAL/IMAG] +//------------------------------------------------------------------------------------------------------- +void ELBDM_RescaleMassError() +{ + if ( ELBDM_InitMassPsi == NULL_REAL ) + { + if ( MPI_Rank == 0 ) + { + ELBDM_InitMassPsi = ConRef[1]; + } + MPI_Bcast( &ELBDM_InitMassPsi, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD ); + } +// check + + if ( ELBDM_MassPsi == NULL_REAL ) + Aux_Error( ERROR_INFO, "ELBDM_MassPsi == NULL_REAL !!\n"); + + if ( ! Aux_IsFinite(ELBDM_MassPsi) ) + Aux_Error( ERROR_INFO, "ELBDM_MassPsi = %14.7e !!\n", ELBDM_MassPsi ); + +// Rescale the total ELBDM mass + for (int lv=0; lvNPatchComma[lv][1]; PID++) + { + real (*const fluid)[PS1][PS1][PS1] = amr->patch[ amr->FluSg[lv] ][lv][PID]->fluid; + + for (int k=0; kuse_wave_flag[lv] ) { +# endif + + fluid[REAL][k][j][i] *= SQRT(ELBDM_InitMassPsi/ELBDM_MassPsi); + fluid[IMAG][k][j][i] *= SQRT(ELBDM_InitMassPsi/ELBDM_MassPsi); + fluid[DENS][k][j][i] = SQR(fluid[REAL][k][j][i]) + SQR(fluid[IMAG][k][j][i]); + +# if ( ELBDM_SCHEME == ELBDM_HYBRID ) + } else { + fluid[DENS][k][j][i] *= (ELBDM_InitMassPsi/ELBDM_MassPsi); + } +# endif + + }}} + } // for (int PID=0; PIDNPatchComma[lv][1]; PID++) + +// update the data on MPI buffer patches +# if ( ELBDM_SCHEME == ELBDM_HYBRID ) + if ( amr->use_wave_flag[lv] ) { +# endif + Buf_GetBufferData( lv, amr->FluSg[lv], NULL_INT, NULL_INT, DATA_GENERAL, _REAL|_IMAG|_DENS, _NONE, Flu_ParaBuf, USELB_YES ); +# if ( ELBDM_SCHEME == ELBDM_HYBRID ) + } else { + Buf_GetBufferData( lv, amr->FluSg[lv], NULL_INT, NULL_INT, DATA_GENERAL, _DENS, _NONE, Flu_ParaBuf, USELB_YES ); + } +# endif + + } // for (int lv=0; lv output OPT__PAR_INIT_CHECK // 2505 : 2025/05/07 --> output PassiveFloor_Var +// 2506 : 2025/12/01 --> output ELBDM_RESCALE_MASS_ERROR, ELBDM_RESCALE_MASS_STEPS //------------------------------------------------------------------------------------------------------- void Output_DumpData_Total_HDF5( const char *FileName ) { @@ -1667,7 +1668,7 @@ void FillIn_KeyInfo( KeyInfo_t &KeyInfo, const int NFieldStored ) const time_t CalTime = time( NULL ); // calendar time - KeyInfo.FormatVersion = 2505; + KeyInfo.FormatVersion = 2506; KeyInfo.Model = MODEL; KeyInfo.NLevel = NLEVEL; KeyInfo.NCompFluid = NCOMP_FLUID; @@ -2580,6 +2581,8 @@ void FillIn_InputPara( InputPara_t &InputPara, const int NFieldStored, char Fiel InputPara.ELBDM_Taylor3_Coeff = ELBDM_TAYLOR3_COEFF; InputPara.ELBDM_Taylor3_Auto = ELBDM_TAYLOR3_AUTO; InputPara.ELBDM_RemoveMotionCM = ELBDM_REMOVE_MOTION_CM; + InputPara.ELBDM_RescaleMassError = ELBDM_RESCALE_MASS_ERROR; + InputPara.ELBDM_RescaleMassSteps = ELBDM_RESCALE_MASS_STEPS; InputPara.ELBDM_BaseSpectral = ELBDM_BASE_SPECTRAL; # if ( ELBDM_SCHEME == ELBDM_HYBRID ) InputPara.ELBDM_FirstWaveLevel = ELBDM_FIRST_WAVE_LEVEL; @@ -3639,6 +3642,8 @@ void GetCompound_InputPara( hid_t &H5_TypeID, const int NFieldStored ) H5Tinsert( H5_TypeID, "ELBDM_Taylor3_Coeff", HOFFSET(InputPara_t,ELBDM_Taylor3_Coeff ), H5T_NATIVE_DOUBLE ); H5Tinsert( H5_TypeID, "ELBDM_Taylor3_Auto", HOFFSET(InputPara_t,ELBDM_Taylor3_Auto ), H5T_NATIVE_INT ); H5Tinsert( H5_TypeID, "ELBDM_RemoveMotionCM", HOFFSET(InputPara_t,ELBDM_RemoveMotionCM ), H5T_NATIVE_INT ); + H5Tinsert( H5_TypeID, "ELBDM_RescaleMassError", HOFFSET(InputPara_t,ELBDM_RescaleMassError ), H5T_NATIVE_INT ); + H5Tinsert( H5_TypeID, "ELBDM_RescaleMassSteps", HOFFSET(InputPara_t,ELBDM_RescaleMassSteps ), H5T_NATIVE_INT ); H5Tinsert( H5_TypeID, "ELBDM_BaseSpectral", HOFFSET(InputPara_t,ELBDM_BaseSpectral ), H5T_NATIVE_INT ); # if ( ELBDM_SCHEME == ELBDM_HYBRID )