Support new rescale mass runtime parameter option

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] |

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##

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

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;

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

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

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" );

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

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

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   );

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 )
 //    ---------------------------------------------------------------------------------------------------
 

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

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; lv<NLEVEL; lv++)
+   {
+#     pragma omp parallel for schedule( runtime )
+      for (int PID=0; PID<amr->NPatchComma[lv][1]; PID++)
+      {
+         real (*const fluid)[PS1][PS1][PS1] = amr->patch[ amr->FluSg[lv] ][lv][PID]->fluid;
+
+         for (int k=0; k<PS1; k++)  {
+         for (int j=0; j<PS1; j++)  {
+         for (int i=0; i<PS1; i++)  {
+
+#           if ( ELBDM_SCHEME == ELBDM_HYBRID )
+            if ( amr->use_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; PID<amr->NPatchComma[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<NLEVEL; lv++)
+
+
+// reset ELBDM_MassPsi[] to check whether it is properly recalculated by Aux_Check_Conservation()
+ELBDM_MassPsi = NULL_REAL;
+
+} // FUNCTION : ELBDM_RescaleMassError
+
+
+
+#endif // #if ( MODEL == ELBDM )

src/Output/Output_DumpData_Total_HDF5.cpp

@@ -79,7 +79,7 @@ Procedure for outputting new variables:
 
 
 //-------------------------------------------------------------------------------------------------------
-// Function    :  Output_DumpData_Total_HDF5 (FormatVersion = 2505)
+// Function    :  Output_DumpData_Total_HDF5 (FormatVersion = 2506)
 // Description :  Output all simulation data in the HDF5 format, which can be used as a restart file
 //                or loaded by YT
 //
@@ -279,6 +279,7 @@ Procedure for outputting new variables:
 //                                             Input__TestProb parameters in "Info/InputTest"
 //                2504 : 2025/04/29 --> 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 )