diff --git a/example/test_problem/ELBDM/ExtPot/Input__Parameter b/example/test_problem/ELBDM/ExtPot/Input__Parameter index e2965934f5..a384134fb9 100644 --- a/example/test_problem/ELBDM/ExtPot/Input__Parameter +++ b/example/test_problem/ELBDM/ExtPot/Input__Parameter @@ -54,7 +54,7 @@ OPT__DT_USER 0 # dt criterion: user-defined -> edit " OPT__DT_LEVEL 3 # dt at different AMR levels (1=shared, 2=differ by two, 3=flexible) [3] OPT__RECORD_DT 1 # record info of the dt determination [1] AUTO_REDUCE_DT 1 # reduce dt automatically when the program fails (for OPT__DT_LEVEL==3 only) [1] -AUTO_REDUCE_DT_FACTOR 0.8 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [0.8] +AUTO_REDUCE_DT_FACTOR 1.0 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [1.0] AUTO_REDUCE_DT_FACTOR_MIN 0.1 # minimum allowed AUTO_REDUCE_DT_FACTOR after consecutive failures [0.1] diff --git a/example/test_problem/Hydro/AGORA_IsolatedGalaxy/Input__Parameter b/example/test_problem/Hydro/AGORA_IsolatedGalaxy/Input__Parameter index 1d202bb30b..e00e4b14d8 100644 --- a/example/test_problem/Hydro/AGORA_IsolatedGalaxy/Input__Parameter +++ b/example/test_problem/Hydro/AGORA_IsolatedGalaxy/Input__Parameter @@ -77,7 +77,7 @@ OPT__DT_USER 0 # dt criterion: user-defined -> edit " OPT__DT_LEVEL 3 # dt at different AMR levels (1=shared, 2=differ by two, 3=flexible) [3] OPT__RECORD_DT 1 # record info of the dt determination [1] AUTO_REDUCE_DT 1 # reduce dt automatically when the program fails (for OPT__DT_LEVEL==3 only) [1] -AUTO_REDUCE_DT_FACTOR 0.8 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [0.8] +AUTO_REDUCE_DT_FACTOR 1.0 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [1.0] AUTO_REDUCE_DT_FACTOR_MIN 0.1 # minimum allowed AUTO_REDUCE_DT_FACTOR after consecutive failures [0.1] AUTO_REDUCE_MINMOD_FACTOR 0.8 # reduce MINMOD_COEFF by this factor together with AUTO_REDUCE_DT (1.0=off) [0.8] ##HYDRO ONLY## AUTO_REDUCE_MINMOD_MIN 1.0e-2 # minimum allowed MINMOD_COEFF after consecutive failures [1.0e-2] ##HYDRO ONLY## diff --git a/example/test_problem/Hydro/AcousticWave/Input__Parameter b/example/test_problem/Hydro/AcousticWave/Input__Parameter index 61b6252f46..bb27ed4c1e 100644 --- a/example/test_problem/Hydro/AcousticWave/Input__Parameter +++ b/example/test_problem/Hydro/AcousticWave/Input__Parameter @@ -44,7 +44,7 @@ DT__FLUID_INIT -1.0 # dt criterion: DT__FLUID at the first OPT__DT_USER 0 # dt criterion: user-defined -> edit "Mis_GetTimeStep_UserCriteria.cpp" [0] OPT__RECORD_DT 1 # record info of the dt determination [1] AUTO_REDUCE_DT 1 # reduce dt automatically when the program fails (for OPT__DT_LEVEL==3 only) [1] -AUTO_REDUCE_DT_FACTOR 0.8 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [0.8] +AUTO_REDUCE_DT_FACTOR 1.0 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [1.0] AUTO_REDUCE_DT_FACTOR_MIN 0.1 # minimum allowed AUTO_REDUCE_DT_FACTOR after consecutive failures [0.1] AUTO_REDUCE_MINMOD_FACTOR 0.8 # reduce MINMOD_COEFF by this factor together with AUTO_REDUCE_DT (1.0=off) [0.8] ##HYDRO ONLY## AUTO_REDUCE_MINMOD_MIN 1.0e-2 # minimum allowed MINMOD_COEFF after consecutive failures [1.0e-2] ##HYDRO ONLY## diff --git a/example/test_problem/Hydro/BlastWave/Input__Parameter b/example/test_problem/Hydro/BlastWave/Input__Parameter index c2581acbf8..353f9bdbc5 100644 --- a/example/test_problem/Hydro/BlastWave/Input__Parameter +++ b/example/test_problem/Hydro/BlastWave/Input__Parameter @@ -56,7 +56,7 @@ OPT__DT_USER 0 # dt criterion: user-defined -> edit " OPT__DT_LEVEL 3 # dt at different AMR levels (1=shared, 2=differ by two, 3=flexible) [3] OPT__RECORD_DT 1 # record info of the dt determination [1] AUTO_REDUCE_DT 1 # reduce dt automatically when the program fails (for OPT__DT_LEVEL==3 only) [1] -AUTO_REDUCE_DT_FACTOR 0.8 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [0.8] +AUTO_REDUCE_DT_FACTOR 1.0 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [1.0] AUTO_REDUCE_DT_FACTOR_MIN 0.1 # minimum allowed AUTO_REDUCE_DT_FACTOR after consecutive failures [0.1] AUTO_REDUCE_MINMOD_FACTOR 0.8 # reduce MINMOD_COEFF by this factor together with AUTO_REDUCE_DT (1.0=off) [0.8] ##HYDRO ONLY## AUTO_REDUCE_MINMOD_MIN 1.0e-2 # minimum allowed MINMOD_COEFF after consecutive failures [1.0e-2] ##HYDRO ONLY## diff --git a/example/test_problem/Hydro/Bondi/Input__Parameter b/example/test_problem/Hydro/Bondi/Input__Parameter index 1fdd25e513..8ab3fde7be 100644 --- a/example/test_problem/Hydro/Bondi/Input__Parameter +++ b/example/test_problem/Hydro/Bondi/Input__Parameter @@ -63,7 +63,7 @@ OPT__DT_USER 0 # dt criterion: user-defined -> edit " OPT__DT_LEVEL 2 # dt at different AMR levels (1=shared, 2=differ by two, 3=flexible) [3] OPT__RECORD_DT 1 # record info of the dt determination [1] AUTO_REDUCE_DT 1 # reduce dt automatically when the program fails (for OPT__DT_LEVEL==3 only) [1] -AUTO_REDUCE_DT_FACTOR 0.8 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [0.8] +AUTO_REDUCE_DT_FACTOR 1.0 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [1.0] AUTO_REDUCE_DT_FACTOR_MIN 0.1 # minimum allowed AUTO_REDUCE_DT_FACTOR after consecutive failures [0.1] AUTO_REDUCE_MINMOD_FACTOR 0.8 # reduce MINMOD_COEFF by this factor together with AUTO_REDUCE_DT (1.0=off) [0.8] ##HYDRO ONLY## AUTO_REDUCE_MINMOD_MIN 1.0e-2 # minimum allowed MINMOD_COEFF after consecutive failures [1.0e-2] ##HYDRO ONLY## diff --git a/example/test_problem/Hydro/CDM_LSS/Input__Parameter b/example/test_problem/Hydro/CDM_LSS/Input__Parameter index 2a5de8db04..9adf3eed39 100644 --- a/example/test_problem/Hydro/CDM_LSS/Input__Parameter +++ b/example/test_problem/Hydro/CDM_LSS/Input__Parameter @@ -83,7 +83,7 @@ OPT__DT_USER 0 # dt criterion: user-defined -> edit " OPT__DT_LEVEL 3 # dt at different AMR levels (1=shared, 2=differ by two, 3=flexible) [3] OPT__RECORD_DT 1 # record info of the dt determination [1] AUTO_REDUCE_DT 1 # reduce dt automatically when the program fails (for OPT__DT_LEVEL==3 only) [1] -AUTO_REDUCE_DT_FACTOR 0.8 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [0.8] +AUTO_REDUCE_DT_FACTOR 1.0 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [1.0] AUTO_REDUCE_DT_FACTOR_MIN 0.1 # minimum allowed AUTO_REDUCE_DT_FACTOR after consecutive failures [0.1] AUTO_REDUCE_MINMOD_FACTOR 0.8 # reduce MINMOD_COEFF by this factor together with AUTO_REDUCE_DT (1.0=off) [0.8] ##HYDRO ONLY## AUTO_REDUCE_MINMOD_MIN 1.0e-2 # minimum allowed MINMOD_COEFF after consecutive failures [1.0e-2] ##HYDRO ONLY## diff --git a/example/test_problem/Hydro/CMZ/Input__Parameter b/example/test_problem/Hydro/CMZ/Input__Parameter index 148d9a3f52..b774ada279 100644 --- a/example/test_problem/Hydro/CMZ/Input__Parameter +++ b/example/test_problem/Hydro/CMZ/Input__Parameter @@ -81,7 +81,7 @@ OPT__DT_USER 0 # dt criterion: user-defined -> edit " OPT__DT_LEVEL 3 # dt at different AMR levels (1=shared, 2=differ by two, 3=flexible) [3] OPT__RECORD_DT 1 # record info of the dt determination [1] AUTO_REDUCE_DT 0 # reduce dt automatically when the program fails (for OPT__DT_LEVEL==3 only) [1] -AUTO_REDUCE_DT_FACTOR 0.8 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [0.8] +AUTO_REDUCE_DT_FACTOR 1.0 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [1.0] AUTO_REDUCE_DT_FACTOR_MIN 0.5 # minimum allowed AUTO_REDUCE_DT_FACTOR after consecutive failures [0.1] AUTO_REDUCE_MINMOD_FACTOR 0.8 # reduce MINMOD_COEFF by this factor together with AUTO_REDUCE_DT (1.0=off) [0.8] ##HYDRO ONLY## AUTO_REDUCE_MINMOD_MIN 1.0e-2 # minimum allowed MINMOD_COEFF after consecutive failures [1.0e-2] ##HYDRO ONLY## diff --git a/example/test_problem/Hydro/Caustic/Input__Parameter b/example/test_problem/Hydro/Caustic/Input__Parameter index 672a56c5c4..5d938c1009 100644 --- a/example/test_problem/Hydro/Caustic/Input__Parameter +++ b/example/test_problem/Hydro/Caustic/Input__Parameter @@ -63,7 +63,7 @@ OPT__DT_USER 0 # dt criterion: user-defined -> edit " OPT__DT_LEVEL 3 # dt at different AMR levels (1=shared, 2=differ by two, 3=flexible) [3] OPT__RECORD_DT 1 # record info of the dt determination [1] AUTO_REDUCE_DT 1 # reduce dt automatically when the program fails (for OPT__DT_LEVEL==3 only) [1] -AUTO_REDUCE_DT_FACTOR 0.8 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [0.8] +AUTO_REDUCE_DT_FACTOR 1.0 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [1.0] AUTO_REDUCE_DT_FACTOR_MIN 0.1 # minimum allowed AUTO_REDUCE_DT_FACTOR after consecutive failures [0.1] AUTO_REDUCE_MINMOD_FACTOR 0.8 # reduce MINMOD_COEFF by this factor together with AUTO_REDUCE_DT (1.0=off) [0.8] ##HYDRO ONLY## AUTO_REDUCE_MINMOD_MIN 1.0e-2 # minimum allowed MINMOD_COEFF after consecutive failures [1.0e-2] ##HYDRO ONLY## diff --git a/example/test_problem/Hydro/ClusterMerger/Input__Flag_User b/example/test_problem/Hydro/ClusterMerger/Input__Flag_User new file mode 100644 index 0000000000..3c83f3cc4d --- /dev/null +++ b/example/test_problem/Hydro/ClusterMerger/Input__Flag_User @@ -0,0 +1,22 @@ +# Level User-defined Criteria + 0 6.0 + 1 6.0 + 2 6.0 + 3 6.0 + 4 6.0 + 5 6.0 + 6 6.0 + 7 6.0 + 8 6.0 + 9 6.0 + 10 6.0 + 11 6.0 + 12 6.0 + 13 6.0 + 14 6.0 + 15 6.0 + 16 6.0 + 17 6.0 + 18 6.0 + 19 6.0 + 20 6.0 diff --git a/example/test_problem/Hydro/ClusterMerger/Input__Parameter b/example/test_problem/Hydro/ClusterMerger/Input__Parameter index 91f9f7d5b3..a735e03474 100644 --- a/example/test_problem/Hydro/ClusterMerger/Input__Parameter +++ b/example/test_problem/Hydro/ClusterMerger/Input__Parameter @@ -79,7 +79,7 @@ OPT__DT_USER 0 # dt criterion: user-defined -> edit " OPT__DT_LEVEL 3 # dt at different AMR levels (1=shared, 2=differ by two, 3=flexible) [3] OPT__RECORD_DT 1 # record info of the dt determination [1] AUTO_REDUCE_DT 1 # reduce dt automatically when the program fails (for OPT__DT_LEVEL==3 only) [1] -AUTO_REDUCE_DT_FACTOR 0.8 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [0.8] +AUTO_REDUCE_DT_FACTOR 1.0 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [1.0] AUTO_REDUCE_DT_FACTOR_MIN 0.1 # minimum allowed AUTO_REDUCE_DT_FACTOR after consecutive failures [0.1] AUTO_REDUCE_MINMOD_FACTOR 0.8 # reduce MINMOD_COEFF by this factor together with AUTO_REDUCE_DT (1.0=off) [0.8] ##HYDRO ONLY## AUTO_REDUCE_MINMOD_MIN 1.0e-2 # minimum allowed MINMOD_COEFF after consecutive failures [1.0e-2] ##HYDRO ONLY## @@ -124,9 +124,20 @@ OPT__MINIMIZE_MPI_BARRIER 0 # minimize MPI barriers to improve loa # (STORE_POT_GHOST, PAR_IMPROVE_ACC=1, OPT__TIMING_BARRIER=0 only; recommend AUTO_REDUCE_DT=0) +# source terms +SRC_DELEPTONIZATION 0 # deleptonization (for simulations of stellar core collapse) [0] ##HYDRO ONLY## +SRC_EXACTCOOLING 0 # exact cooling scheme from Gaspari (2009) [0] ##HYDRO ONLY## +SRC_EC_TEF_N 6001 # number of points for lambda(T) sampling in LOG +SRC_EC_SUBCYCLING 1 # perform subcycling when the cooling time step is small [0=off, 1=on] +SRC_EC_DTCOEF 0.5 # the coefficient of the cooling time step +SRC_USER 0 # user-defined source terms -> edit "Src_User.cpp" [0] +SRC_GPU_NPGROUP -1 # number of patch groups sent into the CPU/GPU source-term solver (<=0=auto) [-1] + + # fluid solver in HYDRO (MODEL==HYDRO only) GAMMA 1.666666666667 # ratio of specific heats (i.e., adiabatic index) [5.0/3.0] MOLECULAR_WEIGHT 0.59242761692650336 # mean molecular weight -> currently only for post-processing [0.6] +MU_NORM 1.6737352238051868e-24 # m_H value, fully consistent with yt MINMOD_COEFF 2.0 # coefficient of the generalized MinMod limiter (1.0~2.0) [1.5] MINMOD_MAX_ITER 0 # maximum number of iterations to reduce MINMOD_COEFF when data reconstruction fails (0=off) [0] OPT__LR_LIMITER 1 # slope limiter of data reconstruction in the MHM/MHM_RP/CTU schemes: @@ -146,9 +157,12 @@ OPT__FIXUP_RESTRICT 1 # correct coarse grids by averaging th OPT__CORR_AFTER_ALL_SYNC -1 # apply various corrections after all levels are synchronized (see "Flu_CorrAfterAllSync"): # (-1=auto, 0=off, 1=every step, 2=before dump) [-1] OPT__NORMALIZE_PASSIVE 1 # ensure "sum(passive_scalar_density) == gas_density" [1] -MIN_DENS 1.0e-07 # minimum mass density (must >= 0.0) [0.0] ##HYDRO, MHD, and ELBDM ONLY## -MIN_PRES 1.0e-09 # minimum pressure (must >= 0.0) [0.0] ##HYDRO and MHD ONLY## -MIN_EINT 1.0e-09 # minimum internal energy (must >= 0.0) [0.0] ##HYDRO and MHD ONLY## +OPT__RESET_FLUID 0 # reset fluid variables after each update -> edit "Flu_ResetByUser.cpp" [0] +OPT__RESET_FLUID_INIT 0 # reset fluid variables during initialization (<0=auto -> OPT__RESET_FLUID, 0=off, 1=on) [-1] +MIN_DENS 0.0 # minimum mass density (must >= 0.0) [0.0] ##HYDRO, MHD, and ELBDM ONLY## +MIN_PRES 0.0 # minimum pressure (must >= 0.0) [0.0] ##HYDRO and MHD ONLY## +MIN_EINT 0.0 # minimum internal energy (must >= 0.0) [0.0] ##HYDRO and MHD ONLY## +MIN_TEMP 1.0e4 # minimum temperature in K (must >= 0.0) [0.0] ##HYDRO and MHD ONLY## MIN_ENTR 0.0 # minimum entropy (must >= 0.0) [0.0] ##HYDRO and MHD ONLY## @@ -228,6 +242,7 @@ OPT__MANUAL_CONTROL 1 # support manually dump data or stop r # (by generating the file DUMP_GAMER_DUMP or STOP_GAMER_STOP) [1] OPT__RECORD_USER 0 # record the user-specified info -> edit "Aux_RecordUser.cpp" [0] OPT__OPTIMIZE_AGGRESSIVE 0 # apply aggressive optimizations (experimental) [0] +OPT__SORT_PATCH_BY_LBIDX 1 # sort patches to improve bitwise reproducibility [SERIAL:0, LOAD_BALACNE:1] # checks diff --git a/example/test_problem/Hydro/ClusterMerger/Input__TestProb b/example/test_problem/Hydro/ClusterMerger/Input__TestProb index a19023d7a8..2f56ca04c5 100644 --- a/example/test_problem/Hydro/ClusterMerger/Input__TestProb +++ b/example/test_problem/Hydro/ClusterMerger/Input__TestProb @@ -1,33 +1,68 @@ # problem-specific runtime parameters Merger_Coll_NumHalos 2 # number of halos +AGN_feedback 0 # turn on/off (1/0) AGN feedback # parameters for cluster 1 Merger_File_Prof1 profile1_gamer.h5 # profile table of cluster 1 Merger_File_Par1 1to3_b0.0_gamerp_1.h5 # particle file of cluster 1 +Unit_R1 1.0 # the unit of length in cluster 1's IC file (in cgs) +Unit_D1 1.0 # the unit of density in cluster 1's IC file (in cgs) +Unit_P1 1.0 # the unit of pressure in cluster 1's IC file (in cgs) Merger_Coll_PosX1 6000.0 # X-center of cluster 1 in kpc Merger_Coll_PosY1 7500.0 # Y-center of cluster 1 in kpc Merger_Coll_VelX1 375.0 # X-velocity of cluster 1 in km/s Merger_Coll_VelY1 0.0 # Y-velocity of cluster 1 in km/s Merger_Coll_IsGas1 1 # If cluster 1 should have gas and not have DM only +Bondi_MassBH1 3.4e8 # black hole mass (in Msun) +Mdot_BH1 0.0 # accretion rate of cluster 1 in Msun/yr +Jet_HalfHeight1 0.70 # jet 1: half height of the cylinder-shape jet source (in kpc) +Jet_Radius1 0.35 # jet 1: radius of the cylinder-shape jet source (in kpc) # parameters for cluster 2 Merger_File_Prof2 profile2_gamer.h5 # profile table of cluster 2 Merger_File_Par2 1to3_b0.0_gamerp_2.h5 # particle file of cluster 2 +Unit_R2 1.0 # the unit of length in cluster 2's IC file (in cgs) +Unit_D2 1.0 # the unit of density in cluster 2's IC file (in cgs) +Unit_P2 1.0 # the unit of pressure in cluster 2's IC file (in cgs) Merger_Coll_PosX2 9000.0 # X-center of cluster 2 in kpc Merger_Coll_PosY2 7500.0 # Y-center of cluster 2 in kpc Merger_Coll_VelX2 -1125.0000000000002 # X-velocity of cluster 2 in km/s Merger_Coll_VelY2 0.0 # Y-velocity of cluster 2 in km/s Merger_Coll_IsGas2 1 # If cluster 2 should have gas and not have DM only +Bondi_MassBH2 3.4e8 # black hole mass (in Msun) +Mdot_BH2 0.0 # accretion rate of cluster 2 in Msun/yr +Jet_HalfHeight2 0.70 # jet 2: half height of the cylinder-shape jet source (in kpc) +Jet_Radius2 0.35 # jet 2: radius of the cylinder-shape jet source (in kpc) # parameters for cluster 3 (disabled for this example) Merger_File_Prof3 none # profile table of cluster 3 Merger_File_Par3 none # particle file of cluster 3 +Unit_R3 1.0 # the unit of length in cluster 3's IC file (in cgs) +Unit_D3 1.0 # the unit of density in cluster 3's IC file (in cgs) +Unit_P3 1.0 # the unit of pressure in cluster 3's IC file (in cgs) Merger_Coll_PosX3 0 # X-center of cluster 3 in kpc Merger_Coll_PosY3 0.0 # Y-center of cluster 3 in kpc Merger_Coll_VelX3 0.0 # X-velocity of cluster 3 in km/s Merger_Coll_VelY3 0.0 # Y-velocity of cluster 3 in km/s Merger_Coll_IsGas3 0 # If cluster 3 should have gas and not have DM only +Bondi_MassBH3 3.4e8 # black hole mass (in Msun) +Mdot_BH3 0.0 # accretion rate of cluster 3 in Msun/yr +Jet_HalfHeight3 0.70 # jet 3: half height of the cylinder-shape jet source (in kpc) +Jet_Radius3 0.35 # jet 3: radius of the cylinder-shape jet source (in kpc) + +# parameters of AGN jet feedback +Accretion_Mode 1 # 1: hot mode; 2: code mode; 3: combine (hot + cold) +eta 1.0 # mass loading factor in jet feedback +eps_f 0.001 # the radiative efficiency in jet feedback +eps_m 0.0 # the fraction of total energy that goes into the thermal energy in jet feedback +R_acc 4.0 # accretion radius: compute the accretion rate (in kpc) +R_dep 4.0 # radius to deplete the accreted gas (in kpc) +JetDirection_case 1 # Methods for choosing the jet direction: 1. Fixed at x-axis; 2. Import from table (generate JetDirection.txt); 3. Align with angular momentum # other parameters Merger_Coll_UseMetals 0 # look for a metal field in the profile files [1] Merger_Coll_LabelCenter 1 # label the particle closest to the center of each cluster [1] +AdjustBHPos 1 # (true/false) --> Adjust the BH position +AdjustBHVel 0 # (true/false) --> Adjust the BH velocity +AdjustPeriod 10 # the time interval of adjustment (in Myr) +fixBH 0 # fix the BH at the simulation box center and set its velocity to be zero (1 cluster only) diff --git a/example/test_problem/Hydro/ClusterMerger/README b/example/test_problem/Hydro/ClusterMerger/README index 200f32199d..7951a5c904 100644 --- a/example/test_problem/Hydro/ClusterMerger/README +++ b/example/test_problem/Hydro/ClusterMerger/README @@ -24,6 +24,12 @@ Default setup: 2. Default resolution ~ 15 kpc (MAX_LEVEL = 3) +AGN feedback: +======================================== +1. Enable OPT__RESET_FLUID in Input__Parameter and "AGN_feedback" in Input__TestProb. +2. To refine regions around the SMBHs, enable OPT__FLAG_USER and set MAX_LEVEL accordingly. + + Note: ======================================== 1. Mimic the merging cluster simulation setup of Flash provided by John ZuHone diff --git a/example/test_problem/Hydro/ClusterMerger/clean.sh b/example/test_problem/Hydro/ClusterMerger/clean.sh index 2df8a60911..fb0539574a 100644 --- a/example/test_problem/Hydro/ClusterMerger/clean.sh +++ b/example/test_problem/Hydro/ClusterMerger/clean.sh @@ -4,4 +4,4 @@ rm -f Record__Note Record__Timing Record__TimeStep Record__PatchCount Record__Du PowerSpec_* Particle_* nohup.out Record__Performance Record__TimingMPI_* \ Record__ParticleCount Record__User Patch_* Record__NCorrUnphy FailedPatchGroup* *.pyc Record__LoadBalance -rm -f Record__Center +rm -f Record__Center BH_variable.bin diff --git a/example/test_problem/Hydro/ClusterMerger/yt_script/plot_profile.py b/example/test_problem/Hydro/ClusterMerger/yt_script/plot_profile.py index 87ed5de56b..2490c15165 100644 --- a/example/test_problem/Hydro/ClusterMerger/yt_script/plot_profile.py +++ b/example/test_problem/Hydro/ClusterMerger/yt_script/plot_profile.py @@ -40,13 +40,13 @@ # load data -ds_gamer = yt.load( "../Data_000008" ) +ds_gamer = yt.load( "../Data_000000" ) #ds_gamer = yt.load( "gamer/Data_000100" ) #ds_flash = yt.load( "flash/fiducial_1to1_b0_hdf5_plt_cnt_0100" ) # create a sphere on the max density location -sp_gamer = ds_gamer.sphere( "max", (2.5e3, "kpc") ) +sp_gamer = ds_gamer.sphere( "max", (1.0e4, "kpc") ) #sp_flash = ds_flash.sphere( "max", (2.5e3, "kpc") ) @@ -144,4 +144,4 @@ # show/save figure plt.savefig( FileOut, bbox_inches='tight', pad_inches=0.05 ) -plt.show() +#plt.show() diff --git a/example/test_problem/Hydro/EnergyPowerSpectrum/Input__Parameter b/example/test_problem/Hydro/EnergyPowerSpectrum/Input__Parameter index 0145347d2d..3570632336 100644 --- a/example/test_problem/Hydro/EnergyPowerSpectrum/Input__Parameter +++ b/example/test_problem/Hydro/EnergyPowerSpectrum/Input__Parameter @@ -42,7 +42,7 @@ DT__FLUID_INIT -1.0 # dt criterion: DT__FLUID at the first OPT__DT_USER 0 # dt criterion: user-defined -> edit "Mis_GetTimeStep_UserCriteria.cpp" [0] OPT__RECORD_DT 1 # record info of the dt determination [1] AUTO_REDUCE_DT 1 # reduce dt automatically when the program fails (for OPT__DT_LEVEL==3 only) [1] -AUTO_REDUCE_DT_FACTOR 0.8 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [0.8] +AUTO_REDUCE_DT_FACTOR 1.0 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [1.0] AUTO_REDUCE_DT_FACTOR_MIN 0.1 # minimum allowed AUTO_REDUCE_DT_FACTOR after consecutive failures [0.1] AUTO_REDUCE_MINMOD_FACTOR 0.8 # reduce MINMOD_COEFF by this factor together with AUTO_REDUCE_DT (1.0=off) [0.8] ##HYDRO ONLY## AUTO_REDUCE_MINMOD_MIN 1.0e-2 # minimum allowed MINMOD_COEFF after consecutive failures [1.0e-2] ##HYDRO ONLY## diff --git a/example/test_problem/Hydro/ExactCooling/Input__Parameter b/example/test_problem/Hydro/ExactCooling/Input__Parameter new file mode 100644 index 0000000000..34fa27ba82 --- /dev/null +++ b/example/test_problem/Hydro/ExactCooling/Input__Parameter @@ -0,0 +1,387 @@ + + +# ================================================================================================================= +# NOTE: +# 1. Comment symbol: # +# 2. [*]: defaults +# 3. Parameters set to "auto" (usually by setting to a negative value) do not have deterministic default values +# and will be set according to the adopted compilation options and/or other runtime parameters +# 4. To add new parameters, please edit "Init/Init_Load_Parameter.cpp" +# 5. All dimensional variables should be set consistently with the code units (set by UNIT_L/M/T/V/D) unless +# otherwise specified (e.g., SF_CREATE_STAR_MIN_GAS_DENS & SF_CREATE_STAR_MIN_STAR_MASS) +# 6. For boolean options: 0/1 -> off/on +# ================================================================================================================= + + +# simulation scale +BOX_SIZE 3e21 # box size along the longest side (in Mpc/h if COMOVING is adopted) +NX0_TOT_X 64 # number of base-level cells along x +NX0_TOT_Y 64 # number of base-level cells along y +NX0_TOT_Z 64 # number of base-level cells along z +OMP_NTHREAD -1 # number of OpenMP threads (<=0=auto -> omp_get_max_threads) [-1] ##OPENMP ONLY## +END_T -1.0 # end physical time (<0=auto -> must be set by test problems or restart) [-1.0] +END_STEP -1 # end step (<0=auto -> must be set by test problems or restart) [-1] + + +# test problems +TESTPROB_ID 123 # test problem ID [0] + # 0: none + # 1: HYDRO blast wave [+MHD] + # 2: HYDRO acoustic wave + # 3: HYDRO Bondi accretion (+GRAVITY) + # 4: HYDRO cluster merger vs. Flash (+GRAVITY & PARTICLE) + # 5: HYDRO AGORA isolated galaxy (+GRAVITY & PARTICLE & STAR_FORMATION & GRACKLE) + # 6: HYDRO caustic wave + # 7: HYDRO spherical collapse (+GRAVITY & COMOVING) + # 8: HYDRO Kelvin Helmholtz instability + # 9: HYDRO Riemann problems [+MHD] + # 10: HYDRO jet(s) + # 11: HYDRO Plummer cloud(s) (+GRAVITY & PARTICLE) + # 12: HYDRO gravity (+GRAVITY) + # 13: HYDRO MHD Arnold-Beltrami-Childress (ABC) flow (+MHD) + # 14: HYDRO MHD Orszag-Tang vortex (+MHD) + # 15: HYDRO MHD linear wave (+MHD) + # 16: HYDRO Jeans instability (+GRAVITY) [+MHD] + # 17: HYDRO particle in equilibrium (+GRAVITY & PARTICLE) + # 19: HYDRO energy power spectrum + # 100: HYDRO CDM cosmological simulation (+GRAVITY & COMOVING & PARTICLE) + # 101: HYDRO Zeldovich pancake collapse (+GRAVITY & COMOVING & PARTICLE) + # 1000: ELBDM external potential (+GRAVITY) + + +# code units (in cgs) +OPT__UNIT 1 # specify code units -> must set exactly 3 basic units below [0] ##USELESS FOR COMOVING## +UNIT_L 1.0 # length unit (<=0 -> set to UNIT_V*UNIT_T or (UNIT_M/UNIT_D)^(1/3)) [-1.0] +UNIT_M 1.0 # mass unit (<=0 -> set to UNIT_D*UNIT_L^3) [-1.0] +UNIT_T 1.0 # time unit (<=0 -> set to UNIT_L/UNIT_V) [-1.0] +UNIT_V -1.0 # velocity unit (<=0 -> set to UNIT_L/UNIT_T) [-1.0] +UNIT_D -1.0 # mass density unit (<=0 -> set to UNIT_M/UNIT_L^3) [-1.0] + + +# boundary conditions +OPT__BC_FLU_XM 1 # fluid boundary condition at the -x face: (1=periodic, 2=outflow, 3=reflecting, 4=user, 5=diode) ##2/3/5 for HYDRO ONLY## +OPT__BC_FLU_XP 1 # fluid boundary condition at the +x face: (1=periodic, 2=outflow, 3=reflecting, 4=user, 5=diode) ##2/3/5 for HYDRO ONLY## +OPT__BC_FLU_YM 1 # fluid boundary condition at the -y face: (1=periodic, 2=outflow, 3=reflecting, 4=user, 5=diode) ##2/3/5 for HYDRO ONLY## +OPT__BC_FLU_YP 1 # fluid boundary condition at the +y face: (1=periodic, 2=outflow, 3=reflecting, 4=user, 5=diode) ##2/3/5 for HYDRO ONLY## +OPT__BC_FLU_ZM 1 # fluid boundary condition at the -z face: (1=periodic, 2=outflow, 3=reflecting, 4=user, 5=diode) ##2/3/5 for HYDRO ONLY## +OPT__BC_FLU_ZP 1 # fluid boundary condition at the +z face: (1=periodic, 2=outflow, 3=reflecting, 4=user, 5=diode) ##2/3/5 for HYDRO ONLY## +OPT__BC_POT 1 # gravity boundary condition: (1=periodic, 2=isolated) +GFUNC_COEFF0 -1.0 # Green's function coefficient at the origin for the isolated BC (<0=auto) [-1.0] + + +# particle (PARTICLE only) +PAR_NPAR 0 # total number of particles (must be set for PAR_INIT==1/3; must be an integer) +PAR_INIT 1 # initialization option for particles: (1=FUNCTION, 2=RESTART, 3=FILE->"PAR_IC") +PAR_IC_FORMAT 1 # data format of PAR_IC: (1=[attribute][id], 2=[id][attribute]; row-major) [1] +PAR_IC_MASS -1.0 # mass of all particles for PAR_INIT==3 (<0=off) [-1.0] +PAR_IC_TYPE -1 # type of all particles for PAR_INIT==3 (<0=off) [-1] +PAR_INTERP 3 # particle interpolation scheme: (1=NGP, 2=CIC, 3=TSC) [2] +PAR_INTEG 2 # particle integration scheme: (1=Euler, 2=KDK) [2] +PAR_TR_INTERP 3 # tracer particle interpolation scheme: (1=NGP, 2=CIC, 3=TSC) [3] +PAR_TR_INTEG 2 # tracer particle integration scheme: (1=Euler, 2=RK2) [2] +PAR_IMPROVE_ACC 1 # improve force accuracy at patch boundaries [1] ##STORE_POT_GHOST and PAR_INTERP=2/3 ONLY## +PAR_PREDICT_POS 1 # predict particle position during mass assignment [1] +PAR_REMOVE_CELL -1.0 # remove particles X-root-cells from the boundaries (non-periodic BC only; <0=auto) [-1.0] +OPT__FREEZE_PAR 0 # do not update particles (except for tracers) [0] +PAR_TR_VEL_CORR 0 # correct tracer particle velocities in regions of discontinuous flow [0] + +# cosmology (COMOVING only) +A_INIT 0.01 # initial scale factor +OMEGA_M0 0.3 # omega matter at the present time +HUBBLE0 0.70 # dimensionless Hubble parameter (currently only for converting ELBDM_MASS to code units) + + +# time-step +DT__MAX -1.0 # dt criterion: maximum allowed dt (<0=off) [-1.0] +DT__FLUID -1.0 # dt criterion: fluid solver CFL factor (<0=auto) [-1.0] +DT__FLUID_INIT -1.0 # dt criterion: DT__FLUID at the first step (<0=auto) [-1.0] +DT__GRAVITY -1.0 # dt criterion: gravity solver safety factor (<0=auto) [-1.0] +DT__PHASE 0.0 # dt criterion: phase rotation safety factor (0=off) [0.0] ##ELBDM ONLY## +DT__PARVEL 0.5 # dt criterion: particle velocity safety factor [0.5] +DT__PARVEL_MAX -1.0 # dt criterion: maximum allowed dt from particle velocity (<0=off) [-1.0] +DT__PARACC 0.5 # dt criterion: particle acceleration safety factor (0=off) [0.5] ##STORE_PAR_ACC ONLY## +DT__MAX_DELTA_A 0.01 # dt criterion: maximum variation of the cosmic scale factor [0.01] +DT__SYNC_PARENT_LV 0.1 # dt criterion: allow dt to adjust by (1.0+DT__SYNC_PARENT) in order to synchronize + # with the parent level (for OPT__DT_LEVEL==3 only) [0.1] +DT__SYNC_CHILDREN_LV 0.1 # dt criterion: allow dt to adjust by (1.0-DT__SYNC_CHILDREN) in order to synchronize + # with the children level (for OPT__DT_LEVEL==3 only; 0=off) [0.1] +OPT__DT_USER 0 # dt criterion: user-defined -> edit "Mis_GetTimeStep_UserCriteria.cpp" [0] +OPT__DT_LEVEL 3 # dt at different AMR levels (1=shared, 2=differ by two, 3=flexible) [3] +OPT__RECORD_DT 1 # record info of the dt determination [1] +AUTO_REDUCE_DT 1 # reduce dt automatically when the program fails (for OPT__DT_LEVEL==3 only) [1] +AUTO_REDUCE_DT_FACTOR 0.8 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [0.8] +AUTO_REDUCE_DT_FACTOR_MIN 0.1 # minimum allowed AUTO_REDUCE_DT_FACTOR after consecutive failures [0.1] +AUTO_REDUCE_MINMOD_FACTOR 0.8 # reduce MINMOD_COEFF by this factor together with AUTO_REDUCE_DT (1.0=off) [0.8] ##HYDRO ONLY## +AUTO_REDUCE_MINMOD_MIN 1.0e-2 # minimum allowed MINMOD_COEFF after consecutive failures [1.0e-2] ##HYDRO ONLY## +AUTO_REDUCE_INT_MONO_FACTOR 0.8 # reduce INT_MONO_COEFF(_B) by this factor together with AUTO_REDUCE_DT (1.0=off) [0.8] +AUTO_REDUCE_INT_MONO_MIN 1.0e-2 # minimum allowed INT_MONO_COEFF(_B) after consecutive failures [1.0e-2] + + +# grid refinement (examples of Input__Flag_XXX tables are put at "example/input/") +REGRID_COUNT 4 # refine every REGRID_COUNT sub-step [4] +REFINE_NLEVEL 1 # number of new AMR levels to be created at once during refinement [1] +FLAG_BUFFER_SIZE -1 # number of buffer cells for the flag operation (0~PATCH_SIZE; <0=auto -> PATCH_SIZE) [-1] +FLAG_BUFFER_SIZE_MAXM1_LV -1 # FLAG_BUFFER_SIZE at the level MAX_LEVEL-1 (<0=auto -> REGRID_COUNT) [-1] +FLAG_BUFFER_SIZE_MAXM2_LV -1 # FLAG_BUFFER_SIZE at the level MAX_LEVEL-2 (<0=auto) [-1] +MAX_LEVEL 0 # maximum refinement level (0~NLEVEL-1) [NLEVEL-1] +OPT__FLAG_RHO 0 # flag: density (Input__Flag_Rho) [0] +OPT__FLAG_RHO_GRADIENT 0 # flag: density gradient (Input__Flag_RhoGradient) [0] +OPT__FLAG_PRES_GRADIENT 0 # flag: pressure gradient (Input__Flag_PresGradient) [0] ##HYDRO ONLY## +OPT__FLAG_VORTICITY 0 # flag: vorticity (Input__Flag_Vorticity) [0] ##HYDRO ONLY## +OPT__FLAG_JEANS 0 # flag: Jeans length (Input__Flag_Jeans) [0] ##HYDRO ONLY## +OPT__FLAG_CURRENT 0 # flag: current density in MHD (Input__Flag_Current) [0] ##MHD ONLY## +OPT__FLAG_ENGY_DENSITY 0 # flag: energy density (Input_Flag_EngyDensity) [0] ##ELBDM ONLY## +OPT__FLAG_LOHNER_DENS 0 # flag: Lohner for mass density (Input__Flag_Lohner) [0] ##BOTH HYDRO AND ELBDM## +OPT__FLAG_LOHNER_ENGY 0 # flag: Lohner for energy density (Input__Flag_Lohner) [0] ##HYDRO ONLY## +OPT__FLAG_LOHNER_PRES 0 # flag: Lohner for pressure (Input__Flag_Lohner) [0] ##HYDRO ONLY## +OPT__FLAG_LOHNER_TEMP 0 # flag: Lohner for temperature (Input__Flag_Lohner) [0] ##HYDRO ONLY## +OPT__FLAG_LOHNER_ENTR 0 # flag: Lohner for entropy (Input__Flag_Lohner) [0] ##HYDRO ONLY## +OPT__FLAG_LOHNER_FORM 2 # form of Lohner: (1=FLASH-1, 2=FLASH-2, 3=form-invariant-1, 4=form-invariant-2) [2] +OPT__FLAG_USER 0 # flag: user-defined (Input__Flag_User) -> edit "Flag_User.cpp" [0] +OPT__FLAG_USER_NUM 1 # number of threshold values in user-defined table (Input__Flag_User) [1] +OPT__FLAG_REGION 0 # flag: specify the regions **allowed** to be refined -> edit "Flag_Region.cpp" [0] +OPT__FLAG_NPAR_PATCH 0 # flag: # of particles per patch (Input__Flag_NParPatch): (0=off, 1=itself, 2=itself+siblings) [0] +OPT__FLAG_NPAR_CELL 0 # flag: # of particles per cell (Input__Flag_NParCell) [0] +OPT__FLAG_PAR_MASS_CELL 0 # flag: total particle mass per cell (Input__Flag_ParMassCell) [0] +OPT__NO_FLAG_NEAR_BOUNDARY 0 # flag: disallow refinement near the boundaries [0] +OPT__PATCH_COUNT 1 # record the # of patches at each level: (0=off, 1=every step, 2=every sub-step) [1] +OPT__PARTICLE_COUNT 1 # record the # of particles at each level: (0=off, 1=every step, 2=every sub-step) [1] +OPT__REUSE_MEMORY 2 # reuse patch memory to reduce memory fragmentation: (0=off, 1=on, 2=aggressive) [2] +OPT__MEMORY_POOL 0 # preallocate patches for OPT__REUSE_MEMORY=1/2 (Input__MemoryPool) [0] + + +# load balance (LOAD_BALANCE only) +LB_INPUT__WLI_MAX 0.1 # weighted-load-imbalance (WLI) threshold for redistributing all patches [0.1] +LB_INPUT__PAR_WEIGHT 0.0 # load-balance weighting of one particle over one cell [0.0] +OPT__RECORD_LOAD_BALANCE 1 # record the load-balance info [1] +OPT__MINIMIZE_MPI_BARRIER 1 # minimize MPI barriers to improve load balance, especially with particles [1] + # (STORE_POT_GHOST, PAR_IMPROVE_ACC=1, OPT__TIMING_BARRIER=0 only; recommend AUTO_REDUCE_DT=0) + + +# source terms +SRC_DELEPTONIZATION 0 # deleptonization (for simulations of stellar core collapse) [0] ##HYDRO ONLY## +SRC_EXACTCOOLING 1 # exact cooling scheme from Gaspari (2009) [0] ##HYDRO ONLY## +SRC_EC_TEF_N 1501 # number of points for lambda(T) sampling in LOG +SRC_EC_SUBCYCLING 1 # perform subcycling when the cooling time step is small [0=off, 1=on] +SRC_EC_DTCOEF 0.5 # the coefficient of the cooling time step +SRC_USER 0 # user-defined source terms -> edit "Src_User.cpp" [0] +SRC_GPU_NPGROUP -1 # number of patch groups sent into the CPU/GPU source-term solver (<=0=auto) [-1] + + +# Grackle library for chemistry and radiative cooling (SUPPORT_GRACKLE only) +GRACKLE_ACTIVATE 1 # enable Grackle [1] +GRACKLE_VERBOSE 1 # map to "grackle_verbose" [1] +GRACKLE_COOLING 1 # ... "with_radiative_cooling" [1] +GRACKLE_PRIMORDIAL 0 # ... "primordial_chemistry" (0=Cloudy, 1/2/3=6-/9-/12-species) [0] + # (must increase NCOMP_PASSIVE_USER by 6/9/12, respectively) +GRACKLE_METAL 0 # ... "metal_cooling" (must increase NCOMP_PASSIVE_USER by 1) [0] +GRACKLE_UV 0 # ... "UVbackground" [0] +GRACKLE_CMB_FLOOR 1 # ... "cmb_temperature_floor" [1] +GRACKLE_PE_HEATING 0 # ... "photoelectric_heating" [0] +GRACKLE_PE_HEATING_RATE 8.5e-26 # ... "photoelectric_heating_rate (in erg/cm^3/s)" [8.5e-26] +GRACKLE_CLOUDY_TABLE CloudyData_noUVB.h5 # "grackle_data_file" +GRACKLE_THREE_BODY_RATE 4 # used Glover+08 rate +GRACKLE_CIE_COOLING 1 # 0: off; 1:on +GRACKLE_H2_OPA_APPROX 1 # H2 opacity from Ripamonti+04; 0:off, 1:Ripomonti+04 +CHE_GPU_NPGROUP -1 # number of patch groups sent into the CPU/GPU Grackle solver (<=0=auto) [-1] + + +# star formation (STAR_FORMATION only) +SF_CREATE_STAR_SCHEME 0 # star formation schemes (0=off, 1=AGORA) [0] +SF_CREATE_STAR_RSEED 123 # random seed [123] +SF_CREATE_STAR_DET_RANDOM -1 # make random numbers deterministic (i.e., independent of OpenMP and MPI, <0=auto) [-1] +SF_CREATE_STAR_MIN_LEVEL 0 # minimum AMR level allowed to form stars (<0=auto -> MAX_LEVEL) [0] +SF_CREATE_STAR_MIN_GAS_DENS 1.0e1 # minimum gas density allowed to form stars (in HI count/cm^3) [1.0e1] +SF_CREATE_STAR_MASS_EFF 1.0e-2 # Gas-to-star mass conversion efficiency [1.0e-2] +SF_CREATE_STAR_MIN_STAR_MASS 0.0 # minimum star particle mass for the stochastical star formation (in Msun) [0.0] +SF_CREATE_STAR_MAX_STAR_MFRAC 0.5 # maximum gas mass fraction allowed to convert to stars per substep [0.5] + + +# fluid solver in HYDRO (MODEL==HYDRO only) +GAMMA 1.666666667 # ratio of specific heats (i.e., adiabatic index) [5.0/3.0] ##EOS_GAMMA ONLY## +MOLECULAR_WEIGHT 0.6 # mean molecular weight [0.6] +MU_NORM 1.6737352238051868e-24 # m_H value, fully consistent with yt +ISO_TEMP 1.0e4 # isothermal temperature in kelvin ##EOS_ISOTHERMAL ONLY## +MINMOD_COEFF 1.5 # coefficient of the generalized MinMod limiter (1.0~2.0) [1.5] +MINMOD_MAX_ITER 0 # maximum number of iterations to reduce MINMOD_COEFF when data reconstruction fails (0=off) [0] +OPT__LR_LIMITER -1 # slope limiter of data reconstruction in the MHM/MHM_RP/CTU schemes: + # (-1=auto, 0=none, 1=vanLeer, 2=generalized MinMod, 3=vanAlbada, 4=vanLeer+generalized MinMod, 6=central) [-1] +OPT__1ST_FLUX_CORR -1 # correct unphysical results (defined by MIN_DENS/PRES) by the 1st-order fluxes: + # (<0=auto, 0=off, 1=3D, 2=3D+1D) [-1] ##MHM/MHM_RP/CTU ONLY## +OPT__1ST_FLUX_CORR_SCHEME -1 # Riemann solver for OPT__1ST_FLUX_CORR (-1=auto, 0=none, 1=Roe, 2=HLLC, 3=HLLE, 4=HLLD) [-1] +DUAL_ENERGY_SWITCH 2.0e-2 # apply dual-energy if E_int/E_kin < DUAL_ENERGY_SWITCH [2.0e-2] ##DUAL_ENERGY ONLY## +OPT__SAME_INTERFACE_B 0 # ensure B field consistency on the shared interfaces between sibling patches (for debugging) [0] ##MHD ONLY## + + +# fluid solver in ELBDM (MODEL==ELBDM only) +ELBDM_MASS 1.0 # particle mass in ev/c^2 (input unit is fixed even when OPT__UNIT or COMOVING is on) +ELBDM_PLANCK_CONST 1.0 # reduced Planck constant (will be overwritten if OPT__UNIT or COMOVING is on) +ELBDM_LAMBDA 1.0 # quartic self-interaction coefficient [1.0] ##QUARTIC_SELF_INTERACTION ONLY## +ELBDM_TAYLOR3_COEFF 0.166666667 # 3rd Taylor expansion coefficient [1.0/6.0] ##USELESS if ELBDM_TAYLOR3_AUTO is on## +ELBDM_TAYLOR3_AUTO 1 # Optimize ELBDM_TAYLOR3_COEFF automatically to minimize the damping at kmax [1] + + +# fluid solvers in all models +FLU_GPU_NPGROUP -1 # number of patch groups sent into the CPU/GPU fluid solver (<=0=auto) [-1] +GPU_NSTREAM -1 # number of CUDA streams for the asynchronous memory copy in GPU (<=0=auto) [-1] +OPT__FIXUP_FLUX 1 # correct coarse grids by the fine-grid boundary fluxes [1] ##HYDRO and ELBDM ONLY## +OPT__FIXUP_ELECTRIC 1 # correct coarse grids by the fine-grid boundary electric field [1] ##MHD ONLY## +OPT__FIXUP_RESTRICT 1 # correct coarse grids by averaging the fine-grid data [1] +OPT__CORR_AFTER_ALL_SYNC -1 # apply various corrections after all levels are synchronized (see "Flu_CorrAfterAllSync"): + # (-1=auto, 0=off, 1=every step, 2=before dump) [-1] +OPT__NORMALIZE_PASSIVE 1 # ensure "sum(passive_scalar_density) == gas_density" [1] +OPT__INT_FRAC_PASSIVE_LR 1 # convert specified passive scalars to mass fraction during data reconstruction [1] +OPT__OVERLAP_MPI 0 # overlap MPI communication with CPU/GPU computations [0] ##NOT SUPPORTED YET## +OPT__RESET_FLUID 0 # reset fluid variables after each update -> edit "Flu_ResetByUser.cpp" [0] +OPT__RESET_FLUID_INIT -1 # reset fluid variables during initialization (<0=auto -> OPT__RESET_FLUID, 0=off, 1=on) [-1] +OPT__FREEZE_FLUID 0 # do not evolve fluid at all [0] +OPT__CHECK_PRES_AFTER_FLU -1 # check unphysical pressure at the end of the fluid solver (<0=auto) [-1] +OPT__LAST_RESORT_FLOOR 1 # apply floor values as the last resort when the fluid solver fails [1] ##HYDRO and MHD ONLY## +MIN_DENS 1.0e-30 # minimum mass density (must >= 0.0) [0.0] ##HYDRO, MHD, and ELBDM ONLY## +MIN_PRES 1.0e-30 # minimum pressure (must >= 0.0) [0.0] ##HYDRO and MHD ONLY## +MIN_EINT 1.0e-30 # minimum internal energy (must >= 0.0) [0.0] ##HYDRO and MHD ONLY## +MIN_TEMP 1.0e2 # minimum temperature in K (must >= 0.0) [0.0] ##HYDRO and MHD ONLY## +MIN_ENTR 0.0 # minimum entropy (must >= 0.0) [0.0] ##HYDRO and MHD ONLY## +JEANS_MIN_PRES 0 # minimum pressure estimated from the Jeans length [0] ##HYDRO/MHD and GRAVITY ONLY## +JEANS_MIN_PRES_LEVEL -1 # for JEANS_MIN_PRES; ensure Jeans length is resolved by JEANS_MIN_PRES_NCELL*dh[JEANS_MIN_PRES_LEVEL] + # (<0=auto -> MAX_LEVEL) [-1] +JEANS_MIN_PRES_NCELL 4 # for JEANS_MIN_PRES; see JEANS_MIN_PRES_LEVEL [4] + + +# gravity solvers in all models +NEWTON_G 1.0 # gravitational constant (will be overwritten if OPT__UNIT or COMOVING is on) +SOR_OMEGA -1.0 # over-relaxation parameter in SOR: (<0=auto) [-1.0] +SOR_MAX_ITER -1 # maximum number of iterations in SOR: (<0=auto) [-1] +SOR_MIN_ITER -1 # minimum number of iterations in SOR: (<0=auto) [-1] +MG_MAX_ITER -1 # maximum number of iterations in multigrid: (<0=auto) [-1] +MG_NPRE_SMOOTH -1 # number of pre-smoothing steps in multigrid: (<0=auto) [-1] +MG_NPOST_SMOOTH -1 # number of post-smoothing steps in multigrid: (<0=auto) [-1] +MG_TOLERATED_ERROR -1.0 # maximum tolerated error in multigrid (<0=auto) [-1.0] +POT_GPU_NPGROUP -1 # number of patch groups sent into the CPU/GPU Poisson solver (<=0=auto) [-1] +OPT__GRA_P5_GRADIENT 0 # 5-points gradient in the Gravity solver (must have GRA/USG_GHOST_SIZE_G>=2) [0] +OPT__SELF_GRAVITY 1 # add self-gravity [1] +OPT__EXT_ACC 0 # add external acceleration (0=off, 1=function, 2=table) [0] ##HYDRO ONLY## + # --> 2 (table) is not supported yet +OPT__EXT_POT 0 # add external potential (0=off, 1=function, 2=table) [0] + # --> for 2 (table), edit the corresponding parameters below too + +# tabular external potential +EXT_POT_TABLE_NAME ExtPotTable # external potential table: filename +EXT_POT_TABLE_NPOINT_X 129 # ... : table size (i.e., number of data points) along x/y/z +EXT_POT_TABLE_NPOINT_Y 129 # +EXT_POT_TABLE_NPOINT_Z 129 # +EXT_POT_TABLE_DH_X 0.0078125 # ... : spatial interval between adjacent data points +EXT_POT_TABLE_DH_Y 0.0078125 # +EXT_POT_TABLE_DH_Z 0.0078125 # +EXT_POT_TABLE_EDGEL_X 0.0 # ... : starting x/y/z coordinates +EXT_POT_TABLE_EDGEL_Y 0.0 # +EXT_POT_TABLE_EDGEL_Z 0.0 # +EXT_POT_TABLE_FLOAT8 -1 # ... : double precision (<0=auto -> FLOAT8, 0=off, 1=on) [-1] + # --> not supported yet; use -1 for now +OPT__GRAVITY_EXTRA_MASS 0 # add extra mass source when computing gravity [0] + + +# initialization +OPT__INIT 1 # initialization option: (1=FUNCTION, 2=RESTART, 3=FILE->"UM_IC") +OPT__INIT_BFIELD_BYFILE 0 # initialize the magnetic field from a vector potential disk file named "B_IC" + # (example python script: tool/inits/gen_vec_pot.py) [0] ##MHD ONLY## +RESTART_LOAD_NRANK 1 # number of parallel I/O (i.e., number of MPI ranks) for restart [1] +OPT__RESTART_RESET 0 # reset some simulation status parameters (e.g., current step and time) during restart [0] +OPT__UM_IC_LEVEL 0 # starting AMR level in UM_IC [0] +OPT__UM_IC_NLEVEL 1 # number of AMR levels UM_IC [1] --> edit "Input__UM_IC_RefineRegion" if >1 +OPT__UM_IC_NVAR -1 # number of variables in UM_IC: (1~NCOMP_TOTAL; <=0=auto) [HYDRO=5+passive/ELBDM=2] +OPT__UM_IC_FORMAT 1 # data format of UM_IC: (1=vzyx, 2=zyxv; row-major and v=field) [1] +OPT__UM_IC_DOWNGRADE 1 # downgrade UM_IC from level OPT__UM_IC_LEVEL to 0 [1] +OPT__UM_IC_REFINE 1 # refine UM_IC from level OPT__UM_IC_LEVEL to MAX_LEVEL [1] +OPT__UM_IC_LOAD_NRANK 1 # number of parallel I/O (i.e., number of MPI ranks) for loading UM_IC [1] +OPT__INIT_RESTRICT 1 # restrict all data during the initialization [1] +OPT__INIT_GRID_WITH_OMP 1 # enable OpenMP when assigning the initial condition of each grid patch [1] +OPT__GPUID_SELECT -1 # GPU ID selection mode: (-3=Laohu, -2=CUDA, -1=MPI rank, >=0=input) [-1] +INIT_SUBSAMPLING_NCELL 0 # perform sub-sampling during initialization: (0=off, >0=# of sub-sampling cells) [0] + + +# interpolation schemes: (-1=auto, 1=MinMod-3D, 2=MinMod-1D, 3=vanLeer, 4=CQuad, 5=Quad, 6=CQuar, 7=Quar) +OPT__INT_TIME 1 # perform "temporal" interpolation for OPT__DT_LEVEL == 2/3 [1] +OPT__INT_PRIM 1 # switch to primitive variables when the interpolation on conserved variables fails [1] ##HYDRO ONLY## +OPT__INT_PHASE 1 # interpolation on phase (does not support MinMod-1D) [1] ##ELBDM ONLY## +OPT__FLU_INT_SCHEME -1 # ghost-zone fluid variables for the fluid solver [-1] +OPT__REF_FLU_INT_SCHEME -1 # newly allocated fluid variables during grid refinement [-1] +OPT__MAG_INT_SCHEME 4 # ghost-zone magnetic field for the MHD solver (2,3,4,6 only) [4] +OPT__REF_MAG_INT_SCHEME 4 # newly allocated magnetic field during grid refinement (2,3,4,6 only) [4] +OPT__POT_INT_SCHEME 4 # ghost-zone potential for the Poisson solver (only supports 4 & 5) [4] +OPT__RHO_INT_SCHEME 4 # ghost-zone mass density for the Poisson solver [4] +OPT__GRA_INT_SCHEME 4 # ghost-zone potential for the gravity solver (for UNSPLIT_GRAVITY as well) [4] +OPT__REF_POT_INT_SCHEME 4 # newly allocated potential during grid refinement [4] +INT_MONO_COEFF 2.0 # coefficient for ensuring the interpolation monotonicity (1.0~4.0) [2.0] +INT_MONO_COEFF_B 2.0 # coefficient for ensuring the interpolation monotonicity of B field (1.0~4.0) [2.0] ##MHD ONLY## +MONO_MAX_ITER 10 # maximum number of iterations to reduce INT_MONO_COEFF when interpolation fails (0=off) [10] +INT_OPP_SIGN_0TH_ORDER 1 # switch to 0th-order interpolation if adjacent values change signs [HYDRO:1; ELBDM:0] + + +# data dump +OPT__OUTPUT_TOTAL 1 # output the simulation snapshot: (0=off, 1=HDF5, 2=C-binary) [1] +OPT__OUTPUT_PART 0 # output a single line or slice: (0=off, 1=xy, 2=yz, 3=xz, 4=x, 5=y, 6=z, 7=diag) [0] +OPT__OUTPUT_USER 0 # output the user-specified data -> edit "Output_User.cpp" [0] +OPT__OUTPUT_PAR_MODE 0 # output the particle data: (0=off, 1=text-file, 2=C-binary) [0] ##PARTICLE ONLY## +OPT__OUTPUT_BASEPS 0 # output the base-level power spectrum [0] +OPT__OUTPUT_BASE 0 # only output the base-level data [0] ##OPT__OUTPUT_PART ONLY## +OPT__OUTPUT_POT 0 # output gravitational potential [0] ##OPT__OUTPUT_TOTAL ONLY## +OPT__OUTPUT_PAR_DENS 1 # output the particle or total mass density on grids: + # (0=off, 1=particle mass density, 2=total mass density) [1] ##OPT__OUTPUT_TOTAL ONLY## +OPT__OUTPUT_CC_MAG 1 # output **cell-centered** magnetic field (necessary for yt analysis) [1] ##MHD ONLY## +OPT__OUTPUT_PRES 0 # output gas pressure [0] ##HYDRO ONLY## +OPT__OUTPUT_TEMP 0 # output gas temperature [0] ##HYDRO ONLY## +OPT__OUTPUT_ENTR 0 # output gas entropy [0] ##HYDRO ONLY## +OPT__OUTPUT_CS 0 # output sound speed [0] ##HYDRO ONLY## +OPT__OUTPUT_DIVVEL 0 # output divergence(velocity) [0] ##HYDRO ONLY## +OPT__OUTPUT_MACH 0 # output mach number [0] ##HYDRO ONLY## +OPT__OUTPUT_DIVMAG 0 # output |divergence(B)*dh/|B|| [0] ##MHD ONLY## +OPT__OUTPUT_USER_FIELD 0 # output user-defined derived fields [0] -> edit "Flu_DerivedField_User.cpp" +OPT__OUTPUT_MODE 2 # (1=const step, 2=const dt, 3=dump table) -> edit "Input__DumpTable" for 3 +OPT__OUTPUT_RESTART 0 # output data immediately after restart [0] +OUTPUT_STEP 5 # output data every OUTPUT_STEP step ##OPT__OUTPUT_MODE==1 ONLY## +OUTPUT_DT 3.15576e13 # output data every OUTPUT_DT time interval ##OPT__OUTPUT_MODE==2 ONLY## +OUTPUT_WALLTIME -1.0 # output data every OUTPUT_WALLTIME walltime (<=0.0=off) [-1.0] +OUTPUT_WALLTIME_UNIT 0 # unit of OUTPUT_WALLTIME (0=second, 1=minute, 2=hour, 3=day) [0] +OUTPUT_PART_X -1.0 # x coordinate for OPT__OUTPUT_PART [-1.0] +OUTPUT_PART_Y -1.0 # y coordinate for OPT__OUTPUT_PART [-1.0] +OUTPUT_PART_Z -1.0 # z coordinate for OPT__OUTPUT_PART [-1.0] +INIT_DUMPID -1 # set the first dump ID (<0=auto) [-1] + + +# yt inline analysis (SUPPORT_LIBYT only) +YT_SCRIPT yt_inline # yt inline analysis script (do not include the ".py" file extension) +YT_VERBOSE 1 # verbose level of yt: (0=off, 1=info, 2=warning, 3=debug) [1] +YT_FIG_BASENAME Fig # figure basename, use libyt default if not set: (default=Fig%9d) + + +# miscellaneous +OPT__VERBOSE 0 # output the simulation progress in detail [0] +OPT__TIMING_BARRIER -1 # synchronize before timing -> more accurate, but may slow down the run (<0=auto) [-1] +OPT__TIMING_BALANCE 0 # record the max/min elapsed time in various code sections for checking load balance [0] +OPT__TIMING_MPI 0 # record the MPI bandwidth achieved in various code sections [0] ##LOAD_BALANCE ONLY## +OPT__RECORD_NOTE 1 # take notes for the general simulation info [1] +OPT__RECORD_UNPHY 1 # record the number of cells with unphysical results being corrected [1] +OPT__RECORD_MEMORY 1 # record the memory consumption [1] +OPT__RECORD_PERFORMANCE 1 # record the code performance [1] +OPT__MANUAL_CONTROL 1 # support manually dump data or stop run during the runtime + # (by generating the file DUMP_GAMER_DUMP or STOP_GAMER_STOP) [1] +OPT__RECORD_USER 0 # record the user-specified info -> edit "Aux_Record_User.cpp" [0] +OPT__OPTIMIZE_AGGRESSIVE 0 # apply aggressive optimizations (experimental) [0] + + +# checks +OPT__CK_REFINE 0 # check the grid refinement [0] +OPT__CK_PROPER_NESTING 0 # check the proper-nesting condition [0] +OPT__CK_CONSERVATION 0 # check the conservation law [0] +OPT__CK_NORMALIZE_PASSIVE 0 # check the normalization of passive scalars [0] ##OPT__NORMALIZE_PASSIVE ONLY## +OPT__CK_RESTRICT 0 # check the data restriction [0] +OPT__CK_FINITE 0 # check if all variables are finite [0] +OPT__CK_PATCH_ALLOCATE 0 # check if all patches are properly allocated [0] +OPT__CK_FLUX_ALLOCATE 0 # check if all flux arrays are properly allocated [0] ##HYDRO and ELBDM ONLY## +OPT__CK_NEGATIVE 0 # check the negative values: (0=off, 1=density, 2=pressure and entropy, 3=both) [0] ##HYDRO ONLY## +OPT__CK_MEMFREE 1.0 # check the free memory in GB (0=off, >0=threshold) [1.0] +OPT__CK_PARTICLE 0 # check the particle allocation [0] +OPT__CK_INTERFACE_B 0 # check the consistency of patch interface B field [0] ##MHD ONLY## +OPT__CK_DIVERGENCE_B 0 # check the divergence-free constraint on B field (0=off, 1=on, 2=on+verbose) [0] ##MHD ONLY## +OPT__CK_INPUT_FLUID 0 # check the input data of the fluid solver [0] diff --git a/example/test_problem/Hydro/ExactCooling/Input__TestProb b/example/test_problem/Hydro/ExactCooling/Input__TestProb new file mode 100644 index 0000000000..33ff1c4b38 --- /dev/null +++ b/example/test_problem/Hydro/ExactCooling/Input__TestProb @@ -0,0 +1,3 @@ +# problem-specific runtime parameters +EC_Temp 10000000.0 # the constant temperature (in K) +EC_Dens 1.0 # the constant number density (in cm^-3) diff --git a/example/test_problem/Hydro/ExactCooling/README b/example/test_problem/Hydro/ExactCooling/README new file mode 100644 index 0000000000..515db42949 --- /dev/null +++ b/example/test_problem/Hydro/ExactCooling/README @@ -0,0 +1,24 @@ +Compilation flags: +======================================== +Enable : HYDRO +Disable: PARTICLE, GRAVITY + + +Default setup: +======================================== +1. BOX_SIZE = 1kpc +2. simulation time = 100 Myr +3. Units: in cgs +4. Enable SRC_EC + + +Note: +======================================== +1. This test problem is for testing the exact cooling scheme (reference: Farber et al. 2018) +2. Modify the source term ExactCooling to have only one branch of the cooling function +3. The analytical solution and simulation error will be recorded in Output__Error file when OPT__OUTPUT_USER is enabled. + + + + + diff --git a/example/test_problem/Hydro/ExactCooling/clean.sh b/example/test_problem/Hydro/ExactCooling/clean.sh new file mode 100644 index 0000000000..dcb84a1d4c --- /dev/null +++ b/example/test_problem/Hydro/ExactCooling/clean.sh @@ -0,0 +1,6 @@ +rm -f Record__Note Record__Timing Record__TimeStep Record__PatchCount Record__Dump Record__MemInfo Record__L1Err \ + Record__Conservation Data* stderr stdout log XYslice* YZslice* XZslice* Xline* Yline* Zline* \ + Diag* BaseXYslice* BaseYZslice* BaseXZslice* BaseXline* BaseYline* BaseZline* BaseDiag* \ + PowerSpec_* Particle_* nohup.out Record__Performance Record__TimingMPI_* \ + Record__ParticleCount Record__User Patch_* Record__NCorrUnphy FailedPatchGroup* *.pyc Record__LoadBalance \ + GRACKLE_INFO Record__DivB diff --git a/example/test_problem/Hydro/Gravity/Input__Parameter b/example/test_problem/Hydro/Gravity/Input__Parameter index bd39aeeeca..5b36d5e09e 100644 --- a/example/test_problem/Hydro/Gravity/Input__Parameter +++ b/example/test_problem/Hydro/Gravity/Input__Parameter @@ -60,7 +60,7 @@ OPT__DT_USER 0 # dt criterion: user-defined -> edit " OPT__DT_LEVEL 3 # dt at different AMR levels (1=shared, 2=differ by two, 3=flexible) [3] OPT__RECORD_DT 1 # record info of the dt determination [1] AUTO_REDUCE_DT 1 # reduce dt automatically when the program fails (for OPT__DT_LEVEL==3 only) [1] -AUTO_REDUCE_DT_FACTOR 0.8 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [0.8] +AUTO_REDUCE_DT_FACTOR 1.0 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [1.0] AUTO_REDUCE_DT_FACTOR_MIN 0.1 # minimum allowed AUTO_REDUCE_DT_FACTOR after consecutive failures [0.1] AUTO_REDUCE_MINMOD_FACTOR 0.8 # reduce MINMOD_COEFF by this factor together with AUTO_REDUCE_DT (1.0=off) [0.8] ##HYDRO ONLY## AUTO_REDUCE_MINMOD_MIN 1.0e-2 # minimum allowed MINMOD_COEFF after consecutive failures [1.0e-2] ##HYDRO ONLY## diff --git a/example/test_problem/Hydro/JeansInstability/Input__Parameter b/example/test_problem/Hydro/JeansInstability/Input__Parameter index c2eb781ab6..3be68b3153 100644 --- a/example/test_problem/Hydro/JeansInstability/Input__Parameter +++ b/example/test_problem/Hydro/JeansInstability/Input__Parameter @@ -55,7 +55,7 @@ OPT__DT_USER 0 # dt criterion: user-defined -> edit " OPT__DT_LEVEL 3 # dt at different AMR levels (1=shared, 2=differ by two, 3=flexible) [3] OPT__RECORD_DT 1 # record info of the dt determination [1] AUTO_REDUCE_DT 1 # reduce dt automatically when the program fails (for OPT__DT_LEVEL==3 only) [1] -AUTO_REDUCE_DT_FACTOR 0.8 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [0.8] +AUTO_REDUCE_DT_FACTOR 1.0 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [1.0] AUTO_REDUCE_DT_FACTOR_MIN 0.1 # minimum allowed AUTO_REDUCE_DT_FACTOR after consecutive failures [0.1] AUTO_REDUCE_MINMOD_FACTOR 0.8 # reduce MINMOD_COEFF by this factor together with AUTO_REDUCE_DT (1.0=off) [0.8] ##HYDRO ONLY## AUTO_REDUCE_MINMOD_MIN 1.0e-2 # minimum allowed MINMOD_COEFF after consecutive failures [1.0e-2] ##HYDRO ONLY## diff --git a/example/test_problem/Hydro/Jet/Input__Parameter b/example/test_problem/Hydro/Jet/Input__Parameter index 6809667f11..68082b2817 100644 --- a/example/test_problem/Hydro/Jet/Input__Parameter +++ b/example/test_problem/Hydro/Jet/Input__Parameter @@ -60,7 +60,7 @@ OPT__DT_USER 0 # dt criterion: user-defined -> edit " OPT__DT_LEVEL 3 # dt at different AMR levels (1=shared, 2=differ by two, 3=flexible) [3] OPT__RECORD_DT 1 # record info of the dt determination [1] AUTO_REDUCE_DT 1 # reduce dt automatically when the program fails (for OPT__DT_LEVEL==3 only) [1] -AUTO_REDUCE_DT_FACTOR 0.8 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [0.8] +AUTO_REDUCE_DT_FACTOR 1.0 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [1.0] AUTO_REDUCE_DT_FACTOR_MIN 0.1 # minimum allowed AUTO_REDUCE_DT_FACTOR after consecutive failures [0.1] AUTO_REDUCE_MINMOD_FACTOR 0.8 # reduce MINMOD_COEFF by this factor together with AUTO_REDUCE_DT (1.0=off) [0.8] ##HYDRO ONLY## AUTO_REDUCE_MINMOD_MIN 1.0e-2 # minimum allowed MINMOD_COEFF after consecutive failures [1.0e-2] ##HYDRO ONLY## diff --git a/example/test_problem/Hydro/KelvinHelmholtzInstability/Input__Parameter b/example/test_problem/Hydro/KelvinHelmholtzInstability/Input__Parameter index d633b1806f..3f756a0e37 100644 --- a/example/test_problem/Hydro/KelvinHelmholtzInstability/Input__Parameter +++ b/example/test_problem/Hydro/KelvinHelmholtzInstability/Input__Parameter @@ -52,7 +52,7 @@ OPT__DT_USER 0 # dt criterion: user-defined -> edit " OPT__DT_LEVEL 3 # dt at different AMR levels (1=shared, 2=differ by two, 3=flexible) [3] OPT__RECORD_DT 1 # record info of the dt determination [1] AUTO_REDUCE_DT 1 # reduce dt automatically when the program fails (for OPT__DT_LEVEL==3 only) [1] -AUTO_REDUCE_DT_FACTOR 0.8 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [0.8] +AUTO_REDUCE_DT_FACTOR 1.0 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [1.0] AUTO_REDUCE_DT_FACTOR_MIN 0.1 # minimum allowed AUTO_REDUCE_DT_FACTOR after consecutive failures [0.1] AUTO_REDUCE_MINMOD_FACTOR 0.8 # reduce MINMOD_COEFF by this factor together with AUTO_REDUCE_DT (1.0=off) [0.8] ##HYDRO ONLY## AUTO_REDUCE_MINMOD_MIN 1.0e-2 # minimum allowed MINMOD_COEFF after consecutive failures [1.0e-2] ##HYDRO ONLY## diff --git a/example/test_problem/Hydro/MHD_ABC/Input__Parameter b/example/test_problem/Hydro/MHD_ABC/Input__Parameter index 27a1882ff3..625fa27959 100644 --- a/example/test_problem/Hydro/MHD_ABC/Input__Parameter +++ b/example/test_problem/Hydro/MHD_ABC/Input__Parameter @@ -57,7 +57,7 @@ OPT__DT_USER 0 # dt criterion: user-defined -> edit " OPT__DT_LEVEL 3 # dt at different AMR levels (1=shared, 2=differ by two, 3=flexible) [3] OPT__RECORD_DT 1 # record info of the dt determination [1] AUTO_REDUCE_DT 1 # reduce dt automatically when the program fails (for OPT__DT_LEVEL==3 only) [1] -AUTO_REDUCE_DT_FACTOR 0.8 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [0.8] +AUTO_REDUCE_DT_FACTOR 1.0 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [1.0] AUTO_REDUCE_DT_FACTOR_MIN 0.1 # minimum allowed AUTO_REDUCE_DT_FACTOR after consecutive failures [0.1] AUTO_REDUCE_MINMOD_FACTOR 0.8 # reduce MINMOD_COEFF by this factor together with AUTO_REDUCE_DT (1.0=off) [0.8] ##HYDRO ONLY## AUTO_REDUCE_MINMOD_MIN 1.0e-2 # minimum allowed MINMOD_COEFF after consecutive failures [1.0e-2] ##HYDRO ONLY## diff --git a/example/test_problem/Hydro/MHD_LinearWave/Input__Parameter b/example/test_problem/Hydro/MHD_LinearWave/Input__Parameter index d05b6d3614..7b2671e186 100644 --- a/example/test_problem/Hydro/MHD_LinearWave/Input__Parameter +++ b/example/test_problem/Hydro/MHD_LinearWave/Input__Parameter @@ -52,7 +52,7 @@ OPT__DT_USER 0 # dt criterion: user-defined -> edit " OPT__DT_LEVEL 3 # dt at different AMR levels (1=shared, 2=differ by two, 3=flexible) [3] OPT__RECORD_DT 1 # record info of the dt determination [1] AUTO_REDUCE_DT 1 # reduce dt automatically when the program fails (for OPT__DT_LEVEL==3 only) [1] -AUTO_REDUCE_DT_FACTOR 0.8 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [0.8] +AUTO_REDUCE_DT_FACTOR 1.0 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [1.0] AUTO_REDUCE_DT_FACTOR_MIN 0.1 # minimum allowed AUTO_REDUCE_DT_FACTOR after consecutive failures [0.1] AUTO_REDUCE_MINMOD_FACTOR 0.8 # reduce MINMOD_COEFF by this factor together with AUTO_REDUCE_DT (1.0=off) [0.8] ##HYDRO ONLY## AUTO_REDUCE_MINMOD_MIN 1.0e-2 # minimum allowed MINMOD_COEFF after consecutive failures [1.0e-2] ##HYDRO ONLY## diff --git a/example/test_problem/Hydro/MHD_OrszagTangVortex/Input__Parameter b/example/test_problem/Hydro/MHD_OrszagTangVortex/Input__Parameter index 0d514794fd..cf3be37905 100644 --- a/example/test_problem/Hydro/MHD_OrszagTangVortex/Input__Parameter +++ b/example/test_problem/Hydro/MHD_OrszagTangVortex/Input__Parameter @@ -57,7 +57,7 @@ OPT__DT_USER 0 # dt criterion: user-defined -> edit " OPT__DT_LEVEL 3 # dt at different AMR levels (1=shared, 2=differ by two, 3=flexible) [3] OPT__RECORD_DT 1 # record info of the dt determination [1] AUTO_REDUCE_DT 1 # reduce dt automatically when the program fails (for OPT__DT_LEVEL==3 only) [1] -AUTO_REDUCE_DT_FACTOR 0.8 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [0.8] +AUTO_REDUCE_DT_FACTOR 1.0 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [1.0] AUTO_REDUCE_DT_FACTOR_MIN 0.1 # minimum allowed AUTO_REDUCE_DT_FACTOR after consecutive failures [0.1] AUTO_REDUCE_MINMOD_FACTOR 0.8 # reduce MINMOD_COEFF by this factor together with AUTO_REDUCE_DT (1.0=off) [0.8] ##HYDRO ONLY## AUTO_REDUCE_MINMOD_MIN 1.0e-2 # minimum allowed MINMOD_COEFF after consecutive failures [1.0e-2] ##HYDRO ONLY## diff --git a/example/test_problem/Hydro/ParticleEquilibriumIC/Input__Parameter b/example/test_problem/Hydro/ParticleEquilibriumIC/Input__Parameter index ab3e6e7943..f05c2a7e9b 100644 --- a/example/test_problem/Hydro/ParticleEquilibriumIC/Input__Parameter +++ b/example/test_problem/Hydro/ParticleEquilibriumIC/Input__Parameter @@ -74,7 +74,7 @@ OPT__DT_USER 0 # dt criterion: user-defined -> edit " OPT__DT_LEVEL 3 # dt at different AMR levels (1=shared, 2=differ by two, 3=flexible) [3] OPT__RECORD_DT 1 # record info of the dt determination [1] AUTO_REDUCE_DT 1 # reduce dt automatically when the program fails (for OPT__DT_LEVEL==3 only) [1] -AUTO_REDUCE_DT_FACTOR 0.8 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [0.8] +AUTO_REDUCE_DT_FACTOR 1.0 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [1.0] AUTO_REDUCE_DT_FACTOR_MIN 0.1 # minimum allowed AUTO_REDUCE_DT_FACTOR after consecutive failures [0.1] AUTO_REDUCE_MINMOD_FACTOR 0.8 # reduce MINMOD_COEFF by this factor together with AUTO_REDUCE_DT (1.0=off) [0.8] ##HYDRO ONLY## AUTO_REDUCE_MINMOD_MIN 1.0e-2 # minimum allowed MINMOD_COEFF after consecutive failures [1.0e-2] ##HYDRO ONLY## diff --git a/example/test_problem/Hydro/ParticleTest/Input__Parameter b/example/test_problem/Hydro/ParticleTest/Input__Parameter index 9648bf5dda..31822b4136 100644 --- a/example/test_problem/Hydro/ParticleTest/Input__Parameter +++ b/example/test_problem/Hydro/ParticleTest/Input__Parameter @@ -70,7 +70,7 @@ OPT__DT_USER 0 # dt criterion: user-defined -> edit " OPT__DT_LEVEL 3 # dt at different AMR levels (1=shared, 2=differ by two, 3=flexible) [3] OPT__RECORD_DT 1 # record info of the dt determination [1] AUTO_REDUCE_DT 1 # reduce dt automatically when the program fails (for OPT__DT_LEVEL==3 only) [1] -AUTO_REDUCE_DT_FACTOR 0.8 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [0.8] +AUTO_REDUCE_DT_FACTOR 1.0 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [1.0] AUTO_REDUCE_DT_FACTOR_MIN 0.1 # minimum allowed AUTO_REDUCE_DT_FACTOR after consecutive failures [0.1] AUTO_REDUCE_MINMOD_FACTOR 0.8 # reduce MINMOD_COEFF by this factor together with AUTO_REDUCE_DT (1.0=off) [0.8] ##HYDRO ONLY## AUTO_REDUCE_MINMOD_MIN 1.0e-2 # minimum allowed MINMOD_COEFF after consecutive failures [1.0e-2] ##HYDRO ONLY## diff --git a/example/test_problem/Hydro/Plummer/Input__Parameter b/example/test_problem/Hydro/Plummer/Input__Parameter index 03b87fd580..c87ce6b2d8 100644 --- a/example/test_problem/Hydro/Plummer/Input__Parameter +++ b/example/test_problem/Hydro/Plummer/Input__Parameter @@ -68,7 +68,7 @@ OPT__DT_USER 0 # dt criterion: user-defined -> edit " OPT__DT_LEVEL 3 # dt at different AMR levels (1=shared, 2=differ by two, 3=flexible) [3] OPT__RECORD_DT 1 # record info of the dt determination [1] AUTO_REDUCE_DT 1 # reduce dt automatically when the program fails (for OPT__DT_LEVEL==3 only) [1] -AUTO_REDUCE_DT_FACTOR 0.8 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [0.8] +AUTO_REDUCE_DT_FACTOR 1.0 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [1.0] AUTO_REDUCE_DT_FACTOR_MIN 0.1 # minimum allowed AUTO_REDUCE_DT_FACTOR after consecutive failures [0.1] AUTO_REDUCE_MINMOD_FACTOR 0.8 # reduce MINMOD_COEFF by this factor together with AUTO_REDUCE_DT (1.0=off) [0.8] ##HYDRO ONLY## AUTO_REDUCE_MINMOD_MIN 1.0e-2 # minimum allowed MINMOD_COEFF after consecutive failures [1.0e-2] ##HYDRO ONLY## diff --git a/example/test_problem/Hydro/Riemann/Input__Parameter b/example/test_problem/Hydro/Riemann/Input__Parameter index bb13f5ffe7..6b5076229a 100644 --- a/example/test_problem/Hydro/Riemann/Input__Parameter +++ b/example/test_problem/Hydro/Riemann/Input__Parameter @@ -52,7 +52,7 @@ OPT__DT_USER 0 # dt criterion: user-defined -> edit " OPT__DT_LEVEL 3 # dt at different AMR levels (1=shared, 2=differ by two, 3=flexible) [3] OPT__RECORD_DT 1 # record info of the dt determination [1] AUTO_REDUCE_DT 1 # reduce dt automatically when the program fails (for OPT__DT_LEVEL==3 only) [1] -AUTO_REDUCE_DT_FACTOR 0.8 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [0.8] +AUTO_REDUCE_DT_FACTOR 1.0 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [1.0] AUTO_REDUCE_DT_FACTOR_MIN 0.1 # minimum allowed AUTO_REDUCE_DT_FACTOR after consecutive failures [0.1] AUTO_REDUCE_MINMOD_FACTOR 0.8 # reduce MINMOD_COEFF by this factor together with AUTO_REDUCE_DT (1.0=off) [0.8] ##HYDRO ONLY## AUTO_REDUCE_MINMOD_MIN 1.0e-2 # minimum allowed MINMOD_COEFF after consecutive failures [1.0e-2] ##HYDRO ONLY## diff --git a/example/test_problem/Hydro/SphericalCollapse/Input__Parameter b/example/test_problem/Hydro/SphericalCollapse/Input__Parameter index c6e9c623b8..4feb821f2f 100644 --- a/example/test_problem/Hydro/SphericalCollapse/Input__Parameter +++ b/example/test_problem/Hydro/SphericalCollapse/Input__Parameter @@ -74,7 +74,7 @@ OPT__DT_USER 0 # dt criterion: user-defined -> edit " OPT__DT_LEVEL 3 # dt at different AMR levels (1=shared, 2=differ by two, 3=flexible) [3] OPT__RECORD_DT 1 # record info of the dt determination [1] AUTO_REDUCE_DT 1 # reduce dt automatically when the program fails (for OPT__DT_LEVEL==3 only) [1] -AUTO_REDUCE_DT_FACTOR 0.8 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [0.8] +AUTO_REDUCE_DT_FACTOR 1.0 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [1.0] AUTO_REDUCE_DT_FACTOR_MIN 0.1 # minimum allowed AUTO_REDUCE_DT_FACTOR after consecutive failures [0.1] AUTO_REDUCE_MINMOD_FACTOR 0.8 # reduce MINMOD_COEFF by this factor together with AUTO_REDUCE_DT (1.0=off) [0.8] ##HYDRO ONLY## AUTO_REDUCE_MINMOD_MIN 1.0e-2 # minimum allowed MINMOD_COEFF after consecutive failures [1.0e-2] ##HYDRO ONLY## diff --git a/example/test_problem/Hydro/Zeldovich/Input__Parameter b/example/test_problem/Hydro/Zeldovich/Input__Parameter index 091c2e22ab..2cc21505a3 100644 --- a/example/test_problem/Hydro/Zeldovich/Input__Parameter +++ b/example/test_problem/Hydro/Zeldovich/Input__Parameter @@ -83,7 +83,7 @@ OPT__DT_USER 0 # dt criterion: user-defined -> edit " OPT__DT_LEVEL 3 # dt at different AMR levels (1=shared, 2=differ by two, 3=flexible) [3] OPT__RECORD_DT 1 # record info of the dt determination [1] AUTO_REDUCE_DT 1 # reduce dt automatically when the program fails (for OPT__DT_LEVEL==3 only) [1] -AUTO_REDUCE_DT_FACTOR 0.8 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [0.8] +AUTO_REDUCE_DT_FACTOR 1.0 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [1.0] AUTO_REDUCE_DT_FACTOR_MIN 0.1 # minimum allowed AUTO_REDUCE_DT_FACTOR after consecutive failures [0.1] diff --git a/example/test_problem/Template/Input__Parameter b/example/test_problem/Template/Input__Parameter index a30031e5ae..bc51af8cfc 100644 --- a/example/test_problem/Template/Input__Parameter +++ b/example/test_problem/Template/Input__Parameter @@ -109,7 +109,7 @@ OPT__DT_USER 0 # dt criterion: user-defined -> edit " OPT__DT_LEVEL 3 # dt at different AMR levels (1=shared, 2=differ by two, 3=flexible) [3] OPT__RECORD_DT 1 # record info of the dt determination [1] AUTO_REDUCE_DT 1 # reduce dt automatically when the program fails (for OPT__DT_LEVEL==3 only) [1] -AUTO_REDUCE_DT_FACTOR 0.8 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [0.8] +AUTO_REDUCE_DT_FACTOR 1.0 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [1.0] AUTO_REDUCE_DT_FACTOR_MIN 0.1 # minimum allowed AUTO_REDUCE_DT_FACTOR after consecutive failures [0.1] AUTO_REDUCE_MINMOD_FACTOR 0.8 # reduce MINMOD_COEFF by this factor together with AUTO_REDUCE_DT (1.0=off) [0.8] ##HYDRO ONLY## AUTO_REDUCE_MINMOD_MIN 1.0e-2 # minimum allowed MINMOD_COEFF after consecutive failures [1.0e-2] ##HYDRO ONLY## @@ -160,6 +160,7 @@ OPT__MINIMIZE_MPI_BARRIER 0 # minimize MPI barriers to improve loa # source terms SRC_DELEPTONIZATION 0 # deleptonization (for simulations of stellar core collapse) [0] ##HYDRO ONLY## +SRC_EXACTCOOLING 0 # exact cooling scheme from Gaspari (2009) [0] ##HYDRO ONLY## SRC_USER 0 # user-defined source terms -> edit "Src_User.cpp" [0] SRC_GPU_NPGROUP -1 # number of patch groups sent into the CPU/GPU source-term solver (<=0=auto) [-1] diff --git a/include/CUDA_ConstMemory.h b/include/CUDA_ConstMemory.h index dae82dfdd2..97c446d205 100644 --- a/include/CUDA_ConstMemory.h +++ b/include/CUDA_ConstMemory.h @@ -34,6 +34,8 @@ SET_GLOBAL( __constant__ real c_Mm[3] ); #if ( MODEL == HYDRO ) SET_GLOBAL( __constant__ double c_Src_Dlep_AuxArray_Flt[SRC_NAUX_DLEP] ); SET_GLOBAL( __constant__ int c_Src_Dlep_AuxArray_Int[SRC_NAUX_DLEP] ); +SET_GLOBAL( __constant__ double c_Src_EC_AuxArray_Flt[SRC_NAUX_EC] ); +SET_GLOBAL( __constant__ int c_Src_EC_AuxArray_Int[SRC_NAUX_EC] ); #endif SET_GLOBAL( __constant__ double c_Src_User_AuxArray_Flt[SRC_NAUX_USER] ); SET_GLOBAL( __constant__ int c_Src_User_AuxArray_Int[SRC_NAUX_USER] ); diff --git a/include/CUFLU.h b/include/CUFLU.h index e5e057dc5c..307ddb0265 100644 --- a/include/CUFLU.h +++ b/include/CUFLU.h @@ -253,6 +253,13 @@ # define HLLD_WAVESPEED HLL_WAVESPEED_DAVIS +// check unphysical results in the MHM half-step prediction +#if ( FLU_SCHEME == MHM ) +# define MHM_CHECK_PREDICT +#endif + + + // 2. ELBDM macro //========================================================================================= #elif ( MODEL == ELBDM ) diff --git a/include/Field.h b/include/Field.h index bbad5d8a84..999c62b1e9 100644 --- a/include/Field.h +++ b/include/Field.h @@ -29,6 +29,7 @@ SET_GLOBAL( FieldIdx_t Idx_CRay, Idx_Undefined ); #ifdef DUAL_ENERGY SET_GLOBAL( FieldIdx_t Idx_Dual, Idx_Undefined ); #endif +SET_GLOBAL( FieldIdx_t Idx_tcool, Idx_Undefined ); // Grackle fields SET_GLOBAL( FieldIdx_t Idx_e, Idx_Undefined ); diff --git a/include/Global.h b/include/Global.h index d3fa25fa5c..8e974fa115 100644 --- a/include/Global.h +++ b/include/Global.h @@ -278,6 +278,8 @@ extern SrcTerms_t SrcTerms; #if ( MODEL == HYDRO ) extern double Src_Dlep_AuxArray_Flt[SRC_NAUX_DLEP]; extern int Src_Dlep_AuxArray_Int[SRC_NAUX_DLEP]; +extern double Src_EC_AuxArray_Flt[SRC_NAUX_EC]; +extern int Src_EC_AuxArray_Int[SRC_NAUX_EC]; #endif extern double Src_User_AuxArray_Flt[SRC_NAUX_USER]; extern int Src_User_AuxArray_Int[SRC_NAUX_USER]; @@ -376,6 +378,12 @@ extern real (*h_SrcDlepProf_Data)[SRC_DLEP_PROF_NBINMAX]; extern real *h_SrcDlepProf_Radius; #endif +#if ( MODEL == HYDRO ) +extern double *h_SrcEC_TEF_lambda; +extern double *h_SrcEC_TEF_alpha; +extern double *h_SrcEC_TEFc; +#endif + // 4/5. GPU (device) global memory arrays and timers diff --git a/include/HDF5_Typedef.h b/include/HDF5_Typedef.h index fbab5b0733..68f4ba43e9 100644 --- a/include/HDF5_Typedef.h +++ b/include/HDF5_Typedef.h @@ -282,6 +282,7 @@ struct SymConst_t # ifdef MHD int EulerY; # endif + int MHM_CheckPredict; int EoSNAuxMax; int EoSNTableMax; diff --git a/include/Macro.h b/include/Macro.h index 441a14d4fc..0fa6c32d24 100644 --- a/include/Macro.h +++ b/include/Macro.h @@ -147,8 +147,11 @@ # define NCOMP_PASSIVE_BUILTIN1 0 # endif +// ExactCooling source term +# define NCOMP_PASSIVE_BUILTIN2 1 + // total number of built-in scalars -# define NCOMP_PASSIVE_BUILTIN ( NCOMP_PASSIVE_BUILTIN0 + NCOMP_PASSIVE_BUILTIN1 ) +# define NCOMP_PASSIVE_BUILTIN ( NCOMP_PASSIVE_BUILTIN0 + NCOMP_PASSIVE_BUILTIN1 + NCOMP_PASSIVE_BUILTIN2 ) #endif // #if ( MODEL == HYDRO ) @@ -244,6 +247,8 @@ # define PASSIVE_NEXT_IDX2 ( PASSIVE_NEXT_IDX1 ) # endif +# define TCOOL ( PASSIVE_NEXT_IDX2 ) + #endif // #if ( NCOMP_PASSIVE > 0 ) // field indices of magnetic --> element of [0 ... NCOMP_MAG-1] @@ -280,6 +285,8 @@ # define FLUX_NEXT_IDX2 ( FLUX_NEXT_IDX1 ) # endif +# define FLUX_TCOOL ( FLUX_NEXT_IDX2 ) + #endif // #if ( NCOMP_PASSIVE > 0 ) // bitwise field indices @@ -302,6 +309,8 @@ # define _CRAY ( 1L << CRAY ) # endif +# define _TCOOL ( 1L << TCOOL ) + #endif // #if ( NCOMP_PASSIVE > 0 ) // magnetic field @@ -331,6 +340,8 @@ # define _FLUX_CRAY ( 1L << FLUX_CRAY ) # endif +# define _FLUX_TCOOL ( 1L << FLUX_TCOOL ) + #endif // #if ( NFLUX_PASSIVE > 0 ) // bitwise indices of derived fields @@ -506,7 +517,7 @@ // particle type macros // number of particle types (default: 4) -# define PAR_NTYPE 6 +# define PAR_NTYPE 8 // particle type indices (must be in the range 0<=index assuming the variable "Eint" is correct diff --git a/src/Fluid/Flu_FixUp_Restrict.cpp b/src/Fluid/Flu_FixUp_Restrict.cpp index 9bf3c1d7a7..73e10be156 100644 --- a/src/Fluid/Flu_FixUp_Restrict.cpp +++ b/src/Fluid/Flu_FixUp_Restrict.cpp @@ -91,6 +91,7 @@ void Flu_FixUp_Restrict( const int FaLv, const int SonFluSg, const int FaFluSg, # else const bool ResMag = false; # endif + const long EoSVar = _TOTAL - _TCOOL; const int PS1_half = PS1 / 2; @@ -300,9 +301,9 @@ void Flu_FixUp_Restrict( const int FaLv, const int SonFluSg, const int FaFluSg, # if ( MODEL == HYDRO ) // apply this correction only when preparing all fluid variables or magnetic field # ifdef MHD - if ( ( TVarCC & _TOTAL ) == _TOTAL || ResMag ) + if ( ( TVarCC & EoSVar ) == EoSVar || ResMag ) # else - if ( ( TVarCC & _TOTAL ) == _TOTAL ) + if ( ( TVarCC & EoSVar ) == EoSVar ) # endif for (int k=0; k as we still rely on these constants (e.g., DENS, DUAL) in the fluid solvers + Idx_tcool = AddField( "tcool", NORMALIZE_NO, INTERP_FRAC_NO ); + if ( Idx_tcool != TCOOL ) Aux_Error( ERROR_INFO, "inconsistent Idx_tcool (%d != %d) !!\n", Idx_tcool, TCOOL ); + # ifdef COSMIC_RAY Idx_CRay = AddField( "CRay", NORMALIZE_NO, INTERP_FRAC_NO ); if ( Idx_CRay != CRAY ) Aux_Error( ERROR_INFO, "inconsistent Idx_CRay (%d != %d) !!\n", Idx_CRay, CRAY ); @@ -121,7 +124,6 @@ void Init_Field() # endif - // 5. validate if all fields have been set properly if ( NDefinedField != NCOMP_TOTAL ) Aux_Error( ERROR_INFO, "total number of defined fields (%d) != expectation (%d) !!\n" diff --git a/src/Init/Init_GAMER.cpp b/src/Init/Init_GAMER.cpp index 273bea8d77..00d054691a 100644 --- a/src/Init/Init_GAMER.cpp +++ b/src/Init/Init_GAMER.cpp @@ -9,7 +9,7 @@ extern void (*Par_Init_ByFunction_Ptr)( const long NPar_ThisRank, const long NPa real *ParType, real *AllAttribute[PAR_NATT_TOTAL] ); #endif - +extern bool IsInit_tcool[NLEVEL]; //------------------------------------------------------------------------------------------------------- diff --git a/src/Init/Init_Load_Parameter.cpp b/src/Init/Init_Load_Parameter.cpp index 0e3cca99f2..ba5fd43556 100644 --- a/src/Init/Init_Load_Parameter.cpp +++ b/src/Init/Init_Load_Parameter.cpp @@ -124,7 +124,7 @@ void Init_Load_Parameter() ReadPara->Add( "OPT__DT_LEVEL", &OPT__DT_LEVEL, 3, 1, 3 ); ReadPara->Add( "OPT__RECORD_DT", &OPT__RECORD_DT, true, Useless_bool, Useless_bool ); ReadPara->Add( "AUTO_REDUCE_DT", &AUTO_REDUCE_DT, true, Useless_bool, Useless_bool ); - ReadPara->Add( "AUTO_REDUCE_DT_FACTOR", &AUTO_REDUCE_DT_FACTOR, 0.8, Eps_double, 1.0 ); + ReadPara->Add( "AUTO_REDUCE_DT_FACTOR", &AUTO_REDUCE_DT_FACTOR, 1.0, Eps_double, 1.0 ); ReadPara->Add( "AUTO_REDUCE_DT_FACTOR_MIN", &AUTO_REDUCE_DT_FACTOR_MIN, 0.1, 0.0, 1.0 ); # if ( MODEL == HYDRO ) ReadPara->Add( "AUTO_REDUCE_MINMOD_FACTOR", &AUTO_REDUCE_MINMOD_FACTOR, 0.8, Eps_double, 1.0 ); @@ -190,6 +190,10 @@ void Init_Load_Parameter() // source terms ReadPara->Add( "SRC_DELEPTONIZATION", &SrcTerms.Deleptonization, false, Useless_bool, Useless_bool ); + ReadPara->Add( "SRC_EXACTCOOLING", &SrcTerms.ExactCooling, false, Useless_bool, Useless_bool ); + ReadPara->Add( "SRC_EC_TEF_N", &SrcTerms.EC_TEF_N, 1501, 1, NoMax_int ); + ReadPara->Add( "SRC_EC_SUBCYCLING", &SrcTerms.EC_subcycling, false, Useless_bool, Useless_bool ); + ReadPara->Add( "SRC_EC_DTCOEF", &SrcTerms.EC_dtCoef, 0.5, Eps_double, NoMax_double ); ReadPara->Add( "SRC_USER", &SrcTerms.User, false, Useless_bool, Useless_bool ); // do not check SRC_GPU_NPGROUP since it may be reset by either Init_ResetDefaultParameter() or CUAPI_SetMemSize() ReadPara->Add( "SRC_GPU_NPGROUP", &SRC_GPU_NPGROUP, -1, NoMin_int, NoMax_int ); diff --git a/src/Init/Init_MemAllocate_Fluid.cpp b/src/Init/Init_MemAllocate_Fluid.cpp index 536dc6ec5a..f2e62473e7 100644 --- a/src/Init/Init_MemAllocate_Fluid.cpp +++ b/src/Init/Init_MemAllocate_Fluid.cpp @@ -92,6 +92,7 @@ void Init_MemAllocate_Fluid( const int Flu_NPatchGroup, const int Pot_NPatchGrou # endif # endif // FLU_SCHEME + } // FUNCTION : Init_MemAllocate_Fluid diff --git a/src/Init/Init_TestProb.cpp b/src/Init/Init_TestProb.cpp index 34818e4260..623cec2308 100644 --- a/src/Init/Init_TestProb.cpp +++ b/src/Init/Init_TestProb.cpp @@ -25,6 +25,7 @@ void Init_TestProb_Hydro_BarredPot(); void Init_TestProb_Hydro_ParticleTest(); void Init_TestProb_Hydro_CDM_LSS(); void Init_TestProb_Hydro_Zeldovich(); +void Init_TestProb_Hydro_ExactCooling(); void Init_TestProb_Hydro_EnergyPowerSpectrum(); void Init_TestProb_ELBDM_ExtPot(); @@ -74,6 +75,7 @@ void Init_TestProb() case TESTPROB_HYDRO_PARTICLE_TEST : Init_TestProb_Hydro_ParticleTest(); break; case TESTPROB_HYDRO_CDM_LSS : Init_TestProb_Hydro_CDM_LSS(); break; case TESTPROB_HYDRO_ZELDOVICH : Init_TestProb_Hydro_Zeldovich(); break; + case TESTPROB_HYDRO_EXACTCOOLING : Init_TestProb_Hydro_ExactCooling(); break; case TESTPROB_HYDRO_ENERGY_POWER_SPECTRUM : Init_TestProb_Hydro_EnergyPowerSpectrum(); break; case TESTPROB_ELBDM_EXTPOT : Init_TestProb_ELBDM_ExtPot(); break; diff --git a/src/Main/EvolveLevel.cpp b/src/Main/EvolveLevel.cpp index f19e41b3f2..17248a86d2 100644 --- a/src/Main/EvolveLevel.cpp +++ b/src/Main/EvolveLevel.cpp @@ -295,6 +295,7 @@ void EvolveLevel( const int lv, const double dTime_FaLv ) # endif // HYDRO Aux_Message( stderr, ", INT_MONO_COEFF %13.7e (min %13.7e)\n", IntMonoCoeff_Failed, AUTO_REDUCE_INT_MONO_MIN ); Aux_Message( stderr, " --> Apply floor values with the original dt and interpolation coefficients as the last resort ...\n" ); + Aux_Message( stderr, " --> Consider setting AUTO_REDUCE_DT_FACTOR < 1.0 in Input__Parameter if not done yet\n" ); } } // if ( AutoReduceDtCoeff >= AUTO_REDUCE_DT_FACTOR_MIN && ... ) ... else ... @@ -719,13 +720,15 @@ void EvolveLevel( const int lv, const double dTime_FaLv ) // 12-1. use the average data on fine grids to correct the coarse-grid data if ( OPT__FIXUP_RESTRICT ) { + const long ResVar = _TOTAL - _TCOOL; + TIMING_FUNC( Flu_FixUp_Restrict( lv, amr->FluSg[lv+1], amr->FluSg[lv], amr->MagSg[lv+1], amr->MagSg[lv], - NULL_INT, NULL_INT, _TOTAL, _MAG ), + NULL_INT, NULL_INT, ResVar, _MAG ), Timer_FixUp[lv], TIMER_ON ); # ifdef LOAD_BALANCE TIMING_FUNC( LB_GetBufferData( lv, amr->FluSg[lv], amr->MagSg[lv], NULL_INT, DATA_RESTRICT, - _TOTAL, _MAG, NULL_INT ), + ResVar, _MAG, NULL_INT ), Timer_GetBuf[lv][7], TIMER_ON ); # endif } diff --git a/src/Main/Main.cpp b/src/Main/Main.cpp index ca688119a8..54c2857c98 100644 --- a/src/Main/Main.cpp +++ b/src/Main/Main.cpp @@ -263,6 +263,8 @@ SrcTerms_t SrcTerms; #if ( MODEL == HYDRO ) double Src_Dlep_AuxArray_Flt[SRC_NAUX_DLEP]; int Src_Dlep_AuxArray_Int[SRC_NAUX_DLEP]; +double Src_EC_AuxArray_Flt[SRC_NAUX_EC]; +int Src_EC_AuxArray_Int[SRC_NAUX_EC]; #endif double Src_User_AuxArray_Flt[SRC_NAUX_USER]; int Src_User_AuxArray_Int[SRC_NAUX_USER]; @@ -364,6 +366,9 @@ double (*h_Corner_Array_S[2])[3] = { NULL, NUL #if ( MODEL == HYDRO ) real (*h_SrcDlepProf_Data)[SRC_DLEP_PROF_NBINMAX] = NULL; real *h_SrcDlepProf_Radius = NULL; +double *h_SrcEC_TEF_lambda = NULL; +double *h_SrcEC_TEF_alpha = NULL; +double *h_SrcEC_TEFc = NULL; #endif @@ -445,6 +450,9 @@ double (*d_Corner_Array_S)[3] = NULL; #if ( MODEL == HYDRO ) real (*d_SrcDlepProf_Data)[SRC_DLEP_PROF_NBINMAX] = NULL; real *d_SrcDlepProf_Radius = NULL; +double *d_SrcEC_TEF_lambda = NULL; +double *d_SrcEC_TEF_alpha = NULL; +double *d_SrcEC_TEFc = NULL; #endif #endif // #ifdef GPU diff --git a/src/Makefile b/src/Makefile index d8391f5e1e..8904684069 100644 --- a/src/Makefile +++ b/src/Makefile @@ -37,10 +37,11 @@ SIMU_OPTION += -DMODEL=HYDRO # enable gravity # --> must define SUPPORT_FFTW=FFTW2 or FFTW3 -#SIMU_OPTION += -DGRAVITY +SIMU_OPTION += -DGRAVITY # enable particles -#SIMU_OPTION += -DPARTICLE +# --> must enable GRAVITY +SIMU_OPTION += -DPARTICLE # support Grackle, a chemistry and radiative cooling library # --> must set NCOMP_PASSIVE_USER according to the primordial chemistry network set by GRACKLE_PRIMORDIAL @@ -55,7 +56,7 @@ ifeq "$(filter -DMODEL=HYDRO, $(SIMU_OPTION))" "-DMODEL=HYDRO" # hydrodynamic scheme: RTVD/MHM/MHM_RP/CTU # --> must be set when MODEL=HYDRO # --> MHD only supports MHM, MHM_RP and CTU -SIMU_OPTION += -DFLU_SCHEME=CTU +SIMU_OPTION += -DFLU_SCHEME=MHM # scheme of spatial data reconstruction: PLM/PPM (piecewise-linear/piecewise-parabolic) # --> useless for RTVD @@ -65,7 +66,7 @@ SIMU_OPTION += -DLR_SCHEME=PPM # --> pure hydro: EXACT/ROE/HLLE/HLLC(*); MHD: ROE/HLLE/HLLD(*); # (recommended solvers are highlighted by *) # --> useless for RTVD -SIMU_OPTION += -DRSOLVER=ROE +SIMU_OPTION += -DRSOLVER=HLLC # dual energy formalism: DE_ENPY/DE_EINT (evolve entropy or internal energy) # --> DE_EINT is not supported yet; useless for RTVD @@ -74,7 +75,7 @@ SIMU_OPTION += -DRSOLVER=ROE # number of user-defined passively advected scalars # --> set it to 0 or comment it out if none is required # --> useless for RTVD -SIMU_OPTION += -DNCOMP_PASSIVE_USER=0 +SIMU_OPTION += -DNCOMP_PASSIVE_USER=1 # magnetohydrodynamics #SIMU_OPTION += -DMHD @@ -170,7 +171,7 @@ SIMU_OPTION += -DPATCH_SIZE=8 # GPU acceleration # --> must set GPU_ARCH as well -#SIMU_OPTION += -DGPU +SIMU_OPTION += -DGPU # GPU architecture: FERMI/KEPLER/MAXWELL/PASCAL/VOLTA/TURING/AMPERE SIMU_OPTION += -DGPU_ARCH=TURING @@ -193,11 +194,11 @@ SIMU_OPTION += -DTIMING # serial mode (in which no MPI libraries are required) # --> must disable LOAD_BALANCE -SIMU_OPTION += -DSERIAL +#SIMU_OPTION += -DSERIAL # load-balance parallelization: HILBERT # --> must disable SERIAL -#SIMU_OPTION += -DLOAD_BALANCE=HILBERT +SIMU_OPTION += -DLOAD_BALANCE=HILBERT # overlap MPI communication with computation # --> NOT supported yet; must enable LOAD_BALANCE @@ -210,16 +211,16 @@ SIMU_OPTION += -DOPENMP #SIMU_OPTION += -DLAOHU # support HDF5 format -#SIMU_OPTION += -DSUPPORT_HDF5 +SIMU_OPTION += -DSUPPORT_HDF5 # support GNU scientific library -#SIMU_OPTION += -DSUPPORT_GSL +SIMU_OPTION += -DSUPPORT_GSL # support FFTW library # --> SUPPORT_FFTW must be defined when GRAVITY is enabled # --> use FFTW3 for fftw3 support # use FFTW2 for fftw2 support -#SIMU_OPTION += -DSUPPORT_FFTW=FFTW3 +SIMU_OPTION += -DSUPPORT_FFTW=FFTW3 # support yt inline analysis @@ -250,58 +251,57 @@ SIMU_OPTION += -DRANDOM_NUMBER=RNG_GNU_EXT # intel # ------------------------------------------------------------------------------- -# CXX = icpc # serial compiler +#CXX = icpc # serial compiler + CXX = $(MPI_PATH)/bin/mpicxx # MPI compiler +#CXX = CC # cray wrapper script + CXXFLAG = -g -O3 # general flags +#CXXFLAG = -g -O3 -std=c++11 #-gxx-name=YOUR_G++ +#CXXFLAG = -g -fast + CXXFLAG += -w1 # warning flags + OPENMPFLAG = -fopenmp # openmp flag + LIB = -limf # libraries and linker flags + +# for debug only +ifeq "$(filter -DGAMER_DEBUG, $(SIMU_OPTION))" "-DGAMER_DEBUG" +#CXXFLAG += -fstack-protector-all +#LIB += -lssp +endif + +# suppress warning when OpenMP is disabled +ifeq "$(filter -DOPENMP, $(SIMU_OPTION))" "" + CXXFLAG += -Wno-unknown-pragmas -diag-disable 3180 +endif + + +# gnu +# ------------------------------------------------------------------------------- +# CXX = g++ # serial compiler ##CXX = $(MPI_PATH)/bin/mpicxx # MPI compiler ##CXX = CC # cray wrapper script # CXXFLAG = -g -O3 # general flags -##CXXFLAG = -g -O3 -std=c++11 #-gxx-name=YOUR_G++ -##CXXFLAG = -g -fast -# CXXFLAG += -w1 # warning flags -# CXXFLAG += -fp-model precise # disable optimizations that are not value-safe +##CXXFLAG = -g -O3 -std=c++11 +##CXXFLAG = -g -Ofast +# CXXFLAG += -Wall -Wextra # warning flags +# CXXFLAG += -Wno-unused-variable -Wno-unused-parameter \ +# -Wno-maybe-uninitialized -Wno-unused-but-set-variable \ +# -Wno-unused-result -Wno-unused-function \ +# -Wno-implicit-fallthrough -Wno-parentheses # OPENMPFLAG = -fopenmp # openmp flag -# LIB = -limf # libraries and linker flags +# LIB = # libraries and linker flags # ## for debug only #ifeq "$(filter -DGAMER_DEBUG, $(SIMU_OPTION))" "-DGAMER_DEBUG" ##CXXFLAG += -fstack-protector-all -##LIB += -lssp +##CXXFLAG += -fsanitize=undefined -fsanitize=address +##LIB += -fsanitize=undefined -fsanitize=address #endif # ## suppress warning when OpenMP is disabled #ifeq "$(filter -DOPENMP, $(SIMU_OPTION))" "" -# CXXFLAG += -Wno-unknown-pragmas -diag-disable 3180 +# CXXFLAG += -Wno-unknown-pragmas #endif -# gnu -# ------------------------------------------------------------------------------- - CXX = g++ # serial compiler -#CXX = $(MPI_PATH)/bin/mpicxx # MPI compiler -#CXX = CC # cray wrapper script - CXXFLAG = -g -O3 # general flags -#CXXFLAG = -g -O3 -std=c++11 -#CXXFLAG = -g -Ofast - CXXFLAG += -Wall -Wextra # warning flags - CXXFLAG += -Wno-unused-variable -Wno-unused-parameter \ - -Wno-maybe-uninitialized -Wno-unused-but-set-variable \ - -Wno-unused-result -Wno-unused-function \ - -Wno-implicit-fallthrough -Wno-parentheses - OPENMPFLAG = -fopenmp # openmp flag - LIB = # libraries and linker flags - -# for debug only -ifeq "$(filter -DGAMER_DEBUG, $(SIMU_OPTION))" "-DGAMER_DEBUG" -#CXXFLAG += -fstack-protector-all -#CXXFLAG += -fsanitize=undefined -fsanitize=address -#LIB += -fsanitize=undefined -fsanitize=address -endif - -# suppress warning when OpenMP is disabled -ifeq "$(filter -DOPENMP, $(SIMU_OPTION))" "" - CXXFLAG += -Wno-unknown-pragmas -endif - - # CUDA # ------------------------------------------------------------------------------- NVCC = $(CUDA_PATH)/bin/nvcc # CUDA compiler @@ -316,24 +316,24 @@ endif ####################################################################################################### # template -CUDA_PATH := -FFTW2_PATH := -FFTW3_PATH := -MPI_PATH := -HDF5_PATH := -GRACKLE_PATH := -GSL_PATH := -LIBYT_PATH := +#CUDA_PATH := +#FFTW_PATH := +#MPI_PATH := +#HDF5_PATH := +#GRACKLE_PATH := +#GSL_PATH := +#LIBYT_PATH := # NTU-eureka (default: openmpi-intel) -#CUDA_PATH := /software/cuda/default +CUDA_PATH := /software/cuda/default +#FFTW_PATH := /software/fftw/default #FFTW2_PATH := /software/fftw/default -#FFTW3_PATH := /software/fftw/default3 -#MPI_PATH := /software/openmpi/default -#HDF5_PATH := /software/hdf5/default -#GRACKLE_PATH := -#GSL_PATH := /software/gsl/default -#LIBYT_PATH := +FFTW3_PATH := /software/fftw/3.3.10-intel-2022.0.1-openmpi-4.1.1-ucx_mt +MPI_PATH := /software/openmpi/default +HDF5_PATH := /software/hdf5/default +GRACKLE_PATH := +GSL_PATH := /software/gsl/default +LIBYT_PATH := # NTU-hulk (openmpi-intel-qlc) #CUDA_PATH := /opt/gpu/cuda/4.2 @@ -645,16 +645,16 @@ endif # SUPPORT_LIBYT # local source terms source files # ------------------------------------------------------------------------------------ GPU_FILE += CUAPI_Asyn_SrcSolver.cu CUSRC_SrcSolver_IterateAllCells.cu CUSRC_Src_Deleptonization.cu \ - CUSRC_Src_User_Template.cu + CUSRC_Src_User_Template.cu CUSRC_Src_ExactCooling.cu CPU_FILE += CPU_SrcSolver.cpp CPU_SrcSolver_IterateAllCells.cpp CPU_Src_Deleptonization.cpp \ - CPU_Src_User_Template.cpp + CPU_Src_User_Template.cpp CPU_Src_ExactCooling.cpp dtSolver_ExactCooling.cpp CPU_FILE += Src_AdvanceDt.cpp Src_Prepare.cpp Src_Close.cpp Src_Init.cpp Src_End.cpp \ Src_WorkBeforeMajorFunc.cpp -vpath %.cu SourceTerms SourceTerms/User_Template SourceTerms/Deleptonization -vpath %.cpp SourceTerms SourceTerms/User_Template SourceTerms/Deleptonization +vpath %.cu SourceTerms SourceTerms/User_Template SourceTerms/Deleptonization SourceTerms/ExactCooling +vpath %.cpp SourceTerms SourceTerms/User_Template SourceTerms/Deleptonization SourceTerms/ExactCooling # Grackle source files (included only if "SUPPORT_GRACKLE" is turned on) diff --git a/src/Miscellaneous/Mis_GetTimeStep.cpp b/src/Miscellaneous/Mis_GetTimeStep.cpp index 56b79213c2..d5d6c22fe5 100644 --- a/src/Miscellaneous/Mis_GetTimeStep.cpp +++ b/src/Miscellaneous/Mis_GetTimeStep.cpp @@ -1,7 +1,7 @@ #include "GAMER.h" extern double (*Mis_GetTimeStep_User_Ptr)( const int lv, const double dTime_dt ); - +extern double Mis_GetTimeStep_ExactCooling( const int lv, const double dTime_dt ); @@ -201,6 +201,21 @@ double Mis_GetTimeStep( const int lv, const double dTime_SyncFaLv, const double # endif +// 1.9 CRITERION NINE : ExactCooling source term ##HYDRO ONLY## +// ============================================================================================================= +# if ( MODEL == HYDRO ) + double EC_dtCoef = SrcTerms.EC_dtCoef; + if ( SrcTerms.EC_subcycling ) { + EC_dtCoef = HUGE_NUMBER; +# ifdef GAMER_DEBUG + Aux_Message( stderr, "WARNING : Resetting EC_dtCoef to be HUGE_NUMBER when subcycling is enabled.\n" ); +# endif + } + dTime[NdTime] = EC_dtCoef * dTime_dt * Mis_GetTimeStep_ExactCooling( lv, dTime_dt ); + sprintf( dTime_Name[NdTime++], "%s", "ExactCooling" ); +# endif + + // 2. get the minimum time-step from all criteria // ============================================================================================================= diff --git a/src/Model_Hydro/CPU_Hydro/CPU_Shared_DataReconstruction.cpp b/src/Model_Hydro/CPU_Hydro/CPU_Shared_DataReconstruction.cpp index 32f08957b2..95ab637168 100644 --- a/src/Model_Hydro/CPU_Hydro/CPU_Shared_DataReconstruction.cpp +++ b/src/Model_Hydro/CPU_Hydro/CPU_Shared_DataReconstruction.cpp @@ -1974,42 +1974,46 @@ void Hydro_HancockPredict( real fc[][NCOMP_LR], const real dt, const real dh, MHD_UpdateMagnetic_Half( fc, g_EC_Ele, dt, dh, cc_i-NGhost, cc_j-NGhost, cc_k-NGhost, NEle ); # endif -// check the negative density and energy + +// check negative, inf, and nan in density, energy, and pressure +# ifdef MHM_CHECK_PREDICT + bool reset_cell = false; for (int f=0; f<6; f++) { -# ifdef BAROTROPIC_EOS - if ( fc[f][0] <= (real)0.0 ) + if ( fc[f][0] <= (real)0.0 || fc[f][0] >= HUGE_NUMBER || fc[f][0] != fc[f][0] ) reset_cell = true; +# ifndef BAROTROPIC_EOS +# ifdef MHD + const real Emag = (real)0.5*( SQR(fc[f][MAG_OFFSET+0]) + SQR(fc[f][MAG_OFFSET+1]) + SQR(fc[f][MAG_OFFSET+2]) ); # else - if ( fc[f][0] <= (real)0.0 || fc[f][4] <= (real)0.0 ) + const real Emag = NULL_REAL; # endif + const real Pres = Hydro_Con2Pres( fc[f][DENS], fc[f][MOMX], fc[f][MOMY], fc[f][MOMZ], fc[f][ENGY], fc[f]+NCOMP_FLUID, + true, MinPres, Emag, EoS->DensEint2Pres_FuncPtr, EoS->AuxArrayDevPtr_Flt, + EoS->AuxArrayDevPtr_Int, EoS->Table, NULL ); + if ( fc[f][4] <= (real)0.0 || fc[f][4] >= HUGE_NUMBER || fc[f][4] != fc[f][4] ) reset_cell = true; + if ( Pres <= (real)0.0 || Pres >= HUGE_NUMBER || Pres != Pres ) reset_cell = true; +# endif // #ifndef BAROTROPIC_EOS + +// set to the cell-centered values before update + if ( reset_cell ) { -// set to the cell-centered values before update - for (int f=0; f<6; f++) + for (int face=0; face<6; face++) for (int v=0; v 0 ) - for (int v=NCOMP_FLUID; v output OPT__FFTW_STARTUP // 2467 : 2023/05/18 --> replace OPT__INIT_BFIELD_BYFILE by OPT__INIT_BFIELD_BYVECPOT // 2468 : 2023/06/24 --> output OPT__SORT_PATCH_BY_LBIDX +// 2469 : 2023/09/09 --> output MHM_CHECK_PREDICT //------------------------------------------------------------------------------------------------------- void Output_DumpData_Total_HDF5( const char *FileName ) { @@ -1414,7 +1415,7 @@ void FillIn_KeyInfo( KeyInfo_t &KeyInfo, const int NFieldStored ) const time_t CalTime = time( NULL ); // calendar time - KeyInfo.FormatVersion = 2468; + KeyInfo.FormatVersion = 2469; KeyInfo.Model = MODEL; KeyInfo.NLevel = NLEVEL; KeyInfo.NCompFluid = NCOMP_FLUID; @@ -1952,6 +1953,11 @@ void FillIn_SymConst( SymConst_t &SymConst ) SymConst.EulerY = 0; # endif # endif // MHD +# ifdef MHM_CHECK_PREDICT + SymConst.MHM_CheckPredict = 1; +# else + SymConst.MHM_CheckPredict = 0; +# endif SymConst.EoSNAuxMax = EOS_NAUX_MAX; SymConst.EoSNTableMax = EOS_NTABLE_MAX; @@ -2799,6 +2805,7 @@ void GetCompound_SymConst( hid_t &H5_TypeID ) # ifdef MHD H5Tinsert( H5_TypeID, "EulerY", HOFFSET(SymConst_t,EulerY ), H5T_NATIVE_INT ); # endif + H5Tinsert( H5_TypeID, "MHM_CheckPredict", HOFFSET(SymConst_t,MHM_CheckPredict ), H5T_NATIVE_INT ); H5Tinsert( H5_TypeID, "EoSNAuxMax", HOFFSET(SymConst_t,EoSNAuxMax ), H5T_NATIVE_INT ); H5Tinsert( H5_TypeID, "EoSNTableMax", HOFFSET(SymConst_t,EoSNTableMax ), H5T_NATIVE_INT ); diff --git a/src/SourceTerms/CPU_SrcSolver_IterateAllCells.cpp b/src/SourceTerms/CPU_SrcSolver_IterateAllCells.cpp index 789d57a38e..54aaaea6a8 100644 --- a/src/SourceTerms/CPU_SrcSolver_IterateAllCells.cpp +++ b/src/SourceTerms/CPU_SrcSolver_IterateAllCells.cpp @@ -112,9 +112,13 @@ void CPU_SrcSolver_IterateAllCells( if ( SrcTerms.Deleptonization ) SrcTerms.Dlep_FuncPtr( fluid, B, &SrcTerms, dt, dh, x, y, z, TimeNew, TimeOld, MinDens, MinPres, MinEint, &EoS, SrcTerms.Dlep_AuxArrayDevPtr_Flt, SrcTerms.Dlep_AuxArrayDevPtr_Int ); +// (2) exact cooling + if ( SrcTerms.ExactCooling ) + SrcTerms.EC_FuncPtr( fluid, B, &SrcTerms, dt, dh, x, y, z, TimeNew, TimeOld, MinDens, MinPres, MinEint, &EoS, + SrcTerms.EC_AuxArrayDevPtr_Flt, SrcTerms.EC_AuxArrayDevPtr_Int ); # endif -// (2) user-defined +// (3) user-defined if ( SrcTerms.User ) SrcTerms.User_FuncPtr( fluid, B, &SrcTerms, dt, dh, x, y, z, TimeNew, TimeOld, MinDens, MinPres, MinEint, &EoS, SrcTerms.User_AuxArrayDevPtr_Flt, SrcTerms.User_AuxArrayDevPtr_Int ); diff --git a/src/SourceTerms/ExactCooling/CPU_Src_ExactCooling.cpp b/src/SourceTerms/ExactCooling/CPU_Src_ExactCooling.cpp new file mode 100644 index 0000000000..7fb47d13da --- /dev/null +++ b/src/SourceTerms/ExactCooling/CPU_Src_ExactCooling.cpp @@ -0,0 +1,725 @@ +#include "CUFLU.h" +#include "Global.h" + +#if ( MODEL == HYDRO ) + + +// external functions and GPU-related set-up +#ifdef __CUDACC__ + +#include "CUDA_CheckError.h" +#include "CUFLU_Shared_FluUtility.cu" +#include "CUDA_ConstMemory.h" +#ifdef DUAL_ENERGY +#include "CUFLU_Shared_DualEnergy.cu" +#endif + +extern double *d_SrcEC_TEF_lambda; +extern double *d_SrcEC_TEF_alpha; +extern double *d_SrcEC_TEFc; + +#endif // #ifdef __CUDACC__ + + +// local function prototypes +#ifndef __CUDACC__ + +extern bool IsInit_tcool[NLEVEL]; +void Src_SetAuxArray_ExactCooling( double [], int [] ); +void Src_SetConstMemory_ExactCooling( const double AuxArray_Flt[], const int AuxArray_Int[], + double *&DevPtr_Flt, int *&DevPtr_Int ); +void Src_PassData2GPU_ExactCooling(); +void Src_SetCPUFunc_ExactCooling( SrcFunc_t & ); +#ifdef GPU +void Src_SetGPUFunc_ExactCooling( SrcFunc_t & ); +#endif +#ifdef GPU +void CUAPI_MemFree_ExactCooling(); +#endif +static bool IsInit = false; +void Src_WorkBeforeMajorFunc_ExactCooling( const int lv, const double TimeNew, const double TimeOld, const double dt, + double AuxArray_Flt[], int AuxArray_Int[] ); +void Src_End_ExactCooling(); + +void Cool_fct( double Dens, double Temp, double* Emis, double* Lambdat, double Z, double cl_moli_mole, double mp ); +#endif +GPU_DEVICE static +double TEF( double TEMP, int k, const double TEF_lambda[], const double TEF_alpha[], const double TEFc[], + const double AuxArray_Flt[], const int AuxArray_Int[] ); +GPU_DEVICE static +double TEFinv( double Y, int k, const double TEF_lambda[], const double TEF_alpha[], const double TEFc[], + const double AuxArray_Flt[], const int AuxArray_Int[] ); +real Hydro_CheckMinTemp( const real InTemp, const real MinTemp ); + + +/******************************************************** +1. Template of a user-defined source term + --> Enabled by the runtime option "SRC_USER" + +2. This file is shared by both CPU and GPU + + CUSRC_Src_ExactCooling.cu -> CPU_Src_ExactCooling.cpp + +3. Four steps are required to implement a source term + + I. Set auxiliary arrays + II. Implement the source-term function + III. [Optional] Add the work to be done every time + before calling the major source-term function + IV. Set initialization functions + +4. The source-term function must be thread-safe and + not use any global variable +********************************************************/ + + + +// ======================= +// I. Set auxiliary arrays +// ======================= + +//------------------------------------------------------------------------------------------------------- +// Function : Src_SetAuxArray_ExactCooling +// Description : Set the auxiliary arrays AuxArray_Flt/Int[] +// +// Note : 1. Invoked by Src_Init_ExactCooling() +// 2. AuxArray_Flt/Int[] have the size of SRC_NAUX_USER defined in Macro.h (default = 10) +// 3. Add "#ifndef __CUDACC__" since this routine is only useful on CPU +// +// Parameter : AuxArray_Flt/Int : Floating-point/Integer arrays to be filled up +// +// Return : AuxArray_Flt/Int[] +//------------------------------------------------------------------------------------------------------- +#ifndef __CUDACC__ +void Src_SetAuxArray_ExactCooling( double AuxArray_Flt[], int AuxArray_Int[] ) +{ + + const int TEF_N = SrcTerms.EC_TEF_N; // number of points for lambda(T) sampling in LOG + const bool subcycling = SrcTerms.EC_subcycling; // whether to use subcycling + const int TEF_int = TEF_N-1; // number of intervals + const double TEF_TN = 1.e14; // == Tref, must be high enough, but affects sampling resolution (Kelvin) + const double TEF_Tmin = MIN_TEMP; // MIN temperature +# ifdef GAMER_DEBUG + if ( TEF_Tmin <= 0.0 ){ + Aux_Error( ERROR_INFO, "TEF_Tmin invalid (can not be smaller or equal to zero)!!\n" ); + } +# endif + const double TEF_dltemp = (log10(TEF_TN) - log10(TEF_Tmin))/TEF_int; // sampling resolution (Kelvin), LOG! + + const double cl_X = 0.7; // mass-fraction of hydrogen + const double cl_Z = 0.018; // metallicity (in Zsun) + const double cl_mol = 1.0/(2*cl_X+0.75*(1-cl_X-cl_Z)+cl_Z*0.5); // mean (total) molecular weights + const double cl_mole = 2.0/(1+cl_X); // mean electron molecular weights + const double cl_moli = 1.0/cl_X; // mean proton molecular weights + const double cl_moli_mole = cl_moli*cl_mole; + +// Store them in the aux array + AuxArray_Flt[0] = 1.0/(GAMMA-1.0); + AuxArray_Flt[1] = TEF_TN; + AuxArray_Flt[2] = TEF_Tmin; + AuxArray_Flt[3] = TEF_dltemp; + AuxArray_Flt[4] = cl_Z; + AuxArray_Flt[5] = cl_moli_mole; + AuxArray_Flt[6] = cl_mol; + AuxArray_Flt[7] = MU_NORM/UNIT_M; //mp: proton mass + AuxArray_Flt[8] = (Const_kB/UNIT_E) * (MU_NORM/UNIT_M); //kB*mp + + AuxArray_Int[0] = TEF_N; + AuxArray_Int[1] = subcycling; + + +} // FUNCTION : Src_SetAuxArray_ExactCooling +#endif // #ifndef __CUDACC__ + + +// ====================================== +// II. Implement the source-term function +// ====================================== + +//------------------------------------------------------------------------------------------------------- +// Function : Src_ExactCooling +// Description : Major source-term function +// +// Note : 1. Invoked by CPU/GPU_SrcSolver_IterateAllCells() +// 2. See Src_SetAuxArray_ExactCooling() for the values stored in AuxArray_Flt/Int[] +// 3. Shared by both CPU and GPU +// +// Parameter : fluid : Fluid array storing both the input and updated values +// --> Including both active and passive variables +// B : Cell-centered magnetic field +// SrcTerms : Structure storing all source-term variables +// dt : Time interval to advance solution +// dh : Grid size +// x/y/z : Target physical coordinates +// TimeNew : Target physical time to reach +// TimeOld : Physical time before update +// --> This function updates physical time from TimeOld to TimeNew +// MinDens/Pres/Eint : Density, pressure, and internal energy floors +// EoS : EoS object +// AuxArray_* : Auxiliary arrays (see the Note above) +// +// Return : fluid[] +//----------------------------------------------------------------------------------------- +GPU_DEVICE_NOINLINE +static void Src_ExactCooling( 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 EoS_t *EoS, const double AuxArray_Flt[], const int AuxArray_Int[] ) +{ + +// check +# ifdef GAMER_DEBUG + if ( AuxArray_Flt == NULL ) printf( "ERROR : AuxArray_Flt == NULL in %s !!\n", __FUNCTION__ ); + if ( AuxArray_Int == NULL ) printf( "ERROR : AuxArray_Int == NULL in %s !!\n", __FUNCTION__ ); +# endif + + + const int TEF_N = AuxArray_Int[0]; // number of points for lambda(T) sampling in LOG + const bool subcycling = AuxArray_Int[1]; // whether to use subcycling + const double cl_CV = AuxArray_Flt[0]; // 1.0/(GAMMA-1.0) + const double TEF_TN = AuxArray_Flt[1]; // == Tref, must be high enough, but affects sampling resolution + const double TEF_Tmin = AuxArray_Flt[2]; // MIN temperature + const double TEF_dltemp = AuxArray_Flt[3]; // sampling resolution (Kelvin), LOG! + +# ifdef __CUDACC__ + const double *TEF_lambda = SrcTerms->EC_TEF_lambda_DevPtr; + const double *TEF_alpha = SrcTerms->EC_TEF_alpha_DevPtr; + const double *TEFc = SrcTerms->EC_TEFc_DevPtr; +# else + const double *TEF_lambda = h_SrcEC_TEF_lambda; + const double *TEF_alpha = h_SrcEC_TEF_alpha; + const double *TEFc = h_SrcEC_TEFc; +# endif + + + double Temp, Eint, Enth, Emag, Pres, Tini, Eintf, Tk, lambdaTini, tcool, Ynew; + int k, knew; + +// (1) Compute the internal energy and temperature + const bool CheckMinEint_No = false; + const bool CheckMinPres_No = false; +# ifdef DUAL_ENERGY +# ifdef __CUDACC__ + Pres = Hydro_DensDual2Pres( fluid[DENS], fluid[DUAL], EoS->AuxArrayDevPtr_Flt[1], CheckMinPres_No, NULL_REAL ); + Eint = EoS->DensPres2Eint_FuncPtr( fluid[DENS], Pres, NULL, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table ); +# else + Pres = Hydro_DensDual2Pres( fluid[DENS], fluid[DUAL], EoS_AuxArray_Flt[1], CheckMinPres_No, NULL_REAL ); + Eint = EoS_DensPres2Eint_CPUPtr( fluid[DENS], Pres, NULL, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); +# endif +# else +# ifdef MHD + Emag = (real)0.5*( SQR(B[MAGX]) + SQR(B[MAGY]) + SQR(B[MAGZ]) ); +# else + Emag = (real)0.0; +# endif + Eint = Hydro_Con2Eint( fluid[DENS], fluid[MOMX], fluid[MOMY], fluid[MOMZ], fluid[ENGY], + CheckMinEint_No, NULL_REAL, Emag ); +# endif + Enth = fluid[ENGY] - Eint; +# ifdef __CUDACC__ + Temp = EoS->DensEint2Temp_FuncPtr( fluid[DENS], Eint, NULL, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table ); +# else + Temp = EoS_DensEint2Temp_CPUPtr( fluid[DENS], Eint, NULL, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); +# endif + Temp = Hydro_CheckMinTemp( Temp, TEF_Tmin ); + Tini = Temp; + +// (2) Decide the index k (an interval) where Tini falls into + k = int((log10(Tini)-log10(TEF_Tmin))/TEF_dltemp); + if ( k < 0 || k > TEF_N-1 || Tini != Tini ){ +# ifdef GAMER_DEBUG + printf( "Error! Temp = %13.7e is out of range (min: %13.7e, max: %13.7e) at TimeNew = %13.7e, so the array index is invalid.\n", Tini, TEF_Tmin, TEF_TN, TimeNew ); +# endif + fluid[TCOOL] = NAN; + fluid[ENGY] = NAN; + return; + } + Tk = POW(10.0, (log10(TEF_Tmin)+k*TEF_dltemp)); + lambdaTini = TEF_lambda[k] * POW((Tini/Tk), TEF_alpha[k]); + +// Compute the cooling time and store it to the field TCOOL + tcool = cl_CV*Tini/(fluid[DENS]*lambdaTini); + fluid[TCOOL] = tcool; + +// Do NOT update the internal energy if dt = 0 + if ( dt == 0.0 ) return; + +// (3) Calculate Ynew + Ynew = TEF( Tini, k, TEF_lambda, TEF_alpha, TEFc, AuxArray_Flt, AuxArray_Int ) + (Tini/TEF_TN)*(TEF_lambda[TEF_N-1]/lambdaTini)*(dt/tcool); + +// (4) Find the new power law interval where Ynew resides + for (int i=k; i>=0; i--){ + if( Ynew < TEFc[i] ){ + knew = i; + Temp = TEFinv( Ynew, knew, TEF_lambda, TEF_alpha, TEFc, AuxArray_Flt, AuxArray_Int ); + if (Temp <= TEF_Tmin){ +# ifdef GAMER_DEBUG + printf( "Error! Temp = %13.7e has reached the floor at TimeNew = %13.7e.\n", Temp, TimeNew ); +# endif + Temp = TEF_Tmin; + } + goto label; + } + } + Temp = TEF_Tmin; // reached the floor: Tn+1 < Tfloor + knew = 0; +# ifdef GAMER_DEBUG + printf( "Error! Temp = %13.7e is reaching the floor at TimeNew = %13.7e.\n", Temp, TimeNew ); +# endif + label: // label for goto statement + +// (5) Calculate the new internal energy and update fluid[ENGY] +# ifdef __CUDACC__ + Pres = EoS->DensTemp2Pres_FuncPtr( fluid[DENS], Temp, NULL, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table ); + Eintf = EoS->DensPres2Eint_FuncPtr( fluid[DENS], Pres, NULL, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table ); +# else + Pres = EoS_DensTemp2Pres_CPUPtr( fluid[DENS], Temp, NULL, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); + Eintf = EoS_DensPres2Eint_CPUPtr( fluid[DENS], Pres, NULL, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); +# endif + + fluid[ENGY] = Enth + Eintf; + +// (6) Update fluid[DUAL] +# ifdef DUAL_ENERGY +# ifdef __CUDACC__ + fluid[DUAL] = Hydro_DensPres2Dual( fluid[DENS], Pres, EoS->AuxArrayDevPtr_Flt[1] ); +# else + fluid[DUAL] = Hydro_DensPres2Dual( fluid[DENS], Pres, EoS_AuxArray_Flt[1] ); +# endif +# endif + +# ifdef GAMER_DEBUG + const real Eintff = real(Eintf); + if ( Hydro_CheckUnphysical( UNPHY_MODE_SING, &Eintff, "output internal energy density", ERROR_INFO, UNPHY_VERBOSE ) ) + { + printf( "Dens = %13.7e, Eint = %13.7e, Eintf = %13.7e, Engy = %13.7e\n", fluid[DENS], Eint, Eintf, fluid[ENGY] ); + printf( "dt = %13.7e, Ynew = %13.7e, tcool = %13.7e, Temp = %13.7e\n", dt, Ynew, tcool, Temp ); + } +# endif // GAMER_DEBUG + +} // FUNCTION : Src_ExactCooling + + +GPU_DEVICE static +// the temporal evolution function (TEF) +double TEF( double TEMP, int k, const double TEF_lambda[], const double TEF_alpha[], const double TEFc[], + const double AuxArray_Flt[], const int AuxArray_Int[] ){ + + const int TEF_N = AuxArray_Int[0]; // number of points for lambda(T) sampling in LOG + const double TEF_TN = AuxArray_Flt[1]; // == Tref, must be high enough, but affects sampling resolution + const double TEF_Tmin = AuxArray_Flt[2]; // MIN temperature + const double TEF_dltemp = AuxArray_Flt[3]; // sampling resolution (Kelvin), LOG! + + double TEF, Tk; + Tk = POW(10.0, log10(TEF_Tmin) + k*TEF_dltemp); +// Do the integration in Gapari (2009) Eq. (24) + if ( TEF_alpha[k] != 1.0 ){ + TEF = TEFc[k] + ((1.0/(1.0-TEF_alpha[k])) * (TEF_lambda[TEF_N-1] / TEF_lambda[k]) * (Tk/TEF_TN) * (1.0 - POW((Tk/TEMP), (TEF_alpha[k]-1.0)))); + } + else TEF = TEFc[k] + ((TEF_lambda[TEF_N-1]/TEF_lambda[k]) * (Tk/TEF_TN) * log(Tk/TEMP)); + + return TEF; +} + + +GPU_DEVICE static +// the INVERSE temporal evolution function (TEF^-1) +double TEFinv( double Y, int k, const double TEF_lambda[], const double TEF_alpha[], const double TEFc[], + const double AuxArray_Flt[], const int AuxArray_Int[] ){ + + const int TEF_N = AuxArray_Int[0]; // number of points for lambda(T) sampling in LOG + const double TEF_TN = AuxArray_Flt[1]; // == Tref, must be high enough, but affects sampling resolution + const double TEF_Tmin = AuxArray_Flt[2]; // MIN temperature + const double TEF_dltemp = AuxArray_Flt[3]; // sampling resolution (Kelvin), LOG! + + double TEFinv, Tk, Yk; + Tk = POW(10.0, log10(TEF_Tmin) + k*TEF_dltemp); + Yk = TEFc[k]; + + if ( TEF_alpha[k] != 1.0 ){ + TEFinv = Tk*POW(1.0-(1.0-TEF_alpha[k])*(TEF_lambda[k]/TEF_lambda[TEF_N-1])*(TEF_TN/Tk)*(double(Y)-Yk), 1.0/(1.0-TEF_alpha[k])); + } + else TEFinv = Tk * exp(-((TEF_lambda[k]/TEF_lambda[TEF_N-1]) * (TEF_TN/Tk) * (double(Y)-Yk))); + + return TEFinv; +} + + + +// ================================================== +// III. [Optional] Add the work to be done every time +// before calling the major source-term function +// ================================================== + +//------------------------------------------------------------------------------------------------------- +// Function : Src_WorkBeforeMajorFunc_ExactCooling +// Description : Specify work to be done every time before calling the major source-term function +// +// Note : 1. Invoked by Src_WorkBeforeMajorFunc() +// --> By linking to "Src_WorkBeforeMajorFunc_EC_Ptr" in Src_Init_ExactCooling() +// 2. Add "#ifndef __CUDACC__" since this routine is only useful on CPU +// +// Parameter : lv : Target refinement level +// TimeNew : Target physical time to reach +// TimeOld : Physical time before update +// --> The major source-term function will update the system from TimeOld to TimeNew +// dt : Time interval to advance solution +// --> Physical coordinates : TimeNew - TimeOld == dt +// Comoving coordinates : TimeNew - TimeOld == delta(scale factor) != dt +// AuxArray_Flt/Int : Auxiliary arrays +// --> Can be used and/or modified here +// --> Must call Src_SetConstMemory_ExactCooling() after modification +// +// Return : AuxArray_Flt/Int[] +//------------------------------------------------------------------------------------------------------- +#ifndef __CUDACC__ +void Src_WorkBeforeMajorFunc_ExactCooling( const int lv, const double TimeNew, const double TimeOld, const double dt, + double AuxArray_Flt[], int AuxArray_Int[] ) +{ + +} // FUNCTION : Src_WorkBeforeMajorFunc_ExactCooling +#endif + + + +#ifdef __CUDACC__ +//------------------------------------------------------------------------------------------------------- +// Function : Src_PassData2GPU_ExactCooling +// Description : Transfer data to GPU +// +// Note : 1. Invoked by Src_WorkBeforeMajorFunc_ExactCooling() +// 2. Use synchronous transfer +// +// Parameter : None +// Return : None +// ------------------------------------------------------------------------------------------------------- +void Src_PassData2GPU_ExactCooling() +{ + const long EC_TEF_MemSize = sizeof(double)*SrcTerms.EC_TEF_N; + + CUDA_CHECK_ERROR( cudaMalloc( (void**) &d_SrcEC_TEF_lambda, EC_TEF_MemSize ) ); + CUDA_CHECK_ERROR( cudaMalloc( (void**) &d_SrcEC_TEF_alpha, EC_TEF_MemSize ) ); + CUDA_CHECK_ERROR( cudaMalloc( (void**) &d_SrcEC_TEFc, EC_TEF_MemSize ) ); + +// store the device pointers in SrcTerms when using GPU + SrcTerms.EC_TEF_lambda_DevPtr = d_SrcEC_TEF_lambda; + SrcTerms.EC_TEF_alpha_DevPtr = d_SrcEC_TEF_alpha; + SrcTerms.EC_TEFc_DevPtr = d_SrcEC_TEFc; + +// use synchronous transfer + CUDA_CHECK_ERROR( cudaMemcpy( d_SrcEC_TEF_lambda, h_SrcEC_TEF_lambda, EC_TEF_MemSize, cudaMemcpyHostToDevice ) ); + CUDA_CHECK_ERROR( cudaMemcpy( d_SrcEC_TEF_alpha, h_SrcEC_TEF_alpha, EC_TEF_MemSize, cudaMemcpyHostToDevice ) ); + CUDA_CHECK_ERROR( cudaMemcpy( d_SrcEC_TEFc, h_SrcEC_TEFc, EC_TEF_MemSize, cudaMemcpyHostToDevice ) ); + +} // FUNCTION : Src_PassData2GPU_ExactCooling +#endif // #ifdef __CUDACC__ + + +// ================================ +// IV. Set initialization functions +// ================================ + +#ifdef __CUDACC__ +# define FUNC_SPACE __device__ static +#else +# define FUNC_SPACE static +#endif + +FUNC_SPACE SrcFunc_t SrcFunc_Ptr = Src_ExactCooling; + +//----------------------------------------------------------------------------------------- +// Function : Src_SetCPU/GPUFunc_ExactCooling +// Description : Return the function pointer of the CPU/GPU source-term function +// +// Note : 1. Invoked by Src_Init_ExactCooling() +// 2. Call-by-reference +// +// Parameter : SrcFunc_CPU/GPUPtr : CPU/GPU function pointer to be set +// +// Return : SrcFunc_CPU/GPUPtr +//----------------------------------------------------------------------------------------- +#ifdef __CUDACC__ +__host__ +void Src_SetGPUFunc_ExactCooling( SrcFunc_t &SrcFunc_GPUPtr ) +{ + CUDA_CHECK_ERROR( cudaMemcpyFromSymbol( &SrcFunc_GPUPtr, SrcFunc_Ptr, sizeof(SrcFunc_t) ) ); +} + +#else + +void Src_SetCPUFunc_ExactCooling( SrcFunc_t &SrcFunc_CPUPtr ) +{ + SrcFunc_CPUPtr = SrcFunc_Ptr; +} + +#endif // #ifdef __CUDACC__ ... else ... + + + +#ifdef __CUDACC__ +//------------------------------------------------------------------------------------------------------- +// Function : Src_SetConstMemory_ExactCooling +// Description : Set the constant memory variables on GPU +// +// Note : 1. Adopt the suggested approach for CUDA version >= 5.0 +// 2. Invoked by Src_Init_ExactCooling() and, if necessary, Src_WorkBeforeMajorFunc_ExactCooling() +// 3. SRC_NAUX_USER is defined in Macro.h +// +// Parameter : AuxArray_Flt/Int : Auxiliary arrays to be copied to the constant memory +// DevPtr_Flt/Int : Pointers to store the addresses of constant memory arrays +// +// Return : c_Src_EC_AuxArray_Flt[], c_Src_EC_AuxArray_Int[], DevPtr_Flt, DevPtr_Int +//--------------------------------------------------------------------------------------------------- +void Src_SetConstMemory_ExactCooling( const double AuxArray_Flt[], const int AuxArray_Int[], + double *&DevPtr_Flt, int *&DevPtr_Int ) +{ + +// copy data to constant memory + CUDA_CHECK_ERROR( cudaMemcpyToSymbol( c_Src_EC_AuxArray_Flt, AuxArray_Flt, SRC_NAUX_USER*sizeof(double) ) ); + CUDA_CHECK_ERROR( cudaMemcpyToSymbol( c_Src_EC_AuxArray_Int, AuxArray_Int, SRC_NAUX_USER*sizeof(int ) ) ); + +// obtain the constant-memory pointers + CUDA_CHECK_ERROR( cudaGetSymbolAddress( (void **)&DevPtr_Flt, c_Src_EC_AuxArray_Flt) ); + CUDA_CHECK_ERROR( cudaGetSymbolAddress( (void **)&DevPtr_Int, c_Src_EC_AuxArray_Int) ); + +} // FUNCTION : Src_SetConstMemory_ExactCooling +#endif // #ifdef __CUDACC__ + + + +#ifndef __CUDACC__ + +// function pointer +//extern void (*Src_WorkBeforeMajorFunc_EC_Ptr)( const int lv, const double TimeNew, const double TimeOld, const double dt, +// double AuxArray_Flt[], int AuxArray_Int[] ); +//extern void (*Src_End_EC_Ptr)(); + +//----------------------------------------------------------------------------------------- +// Function : Src_Init_ExactCooling +// Description : Initialize a user-specified source term +// +// Note : 1. Set auxiliary arrays by invoking Src_SetAuxArray_*() +// --> Copy to the GPU constant memory and store the associated addresses +// 2. Set the source-term function by invoking Src_SetCPU/GPUFunc_*() +// 3. Set the function pointers "Src_WorkBeforeMajorFunc_EC_Ptr" and "Src_End_EC_Ptr" +// 4. Invoked by Src_Init() +// --> Enable it by linking to the function pointer "Src_Init_EC_Ptr" +// 5. Add "#ifndef __CUDACC__" since this routine is only useful on CPU +// +// Parameter : None +// +// Return : None +//----------------------------------------------------------------------------------------- +void Src_Init_ExactCooling() +{ + +// set the auxiliary arrays + Src_SetAuxArray_ExactCooling( Src_EC_AuxArray_Flt, Src_EC_AuxArray_Int ); + +// copy the auxiliary arrays to the GPU constant memory and store the associated addresses +# ifdef GPU + Src_SetConstMemory_ExactCooling( Src_EC_AuxArray_Flt, Src_EC_AuxArray_Int, + SrcTerms.EC_AuxArrayDevPtr_Flt, SrcTerms.EC_AuxArrayDevPtr_Int ); +# else + SrcTerms.EC_AuxArrayDevPtr_Flt = Src_EC_AuxArray_Flt; + SrcTerms.EC_AuxArrayDevPtr_Int = Src_EC_AuxArray_Int; +# endif + +// set the major source-term function + Src_SetCPUFunc_ExactCooling( SrcTerms.EC_CPUPtr ); + +# ifdef GPU + Src_SetGPUFunc_ExactCooling( SrcTerms.EC_GPUPtr ); + SrcTerms.EC_FuncPtr = SrcTerms.EC_GPUPtr; +# else + SrcTerms.EC_FuncPtr = SrcTerms.EC_CPUPtr; +# endif + + if ( OPT__INIT == INIT_BY_RESTART ){ + for (int i=0; i=0; i--){ + Ti = POW(10, log10(TEF_Tmin) + i*TEF_dltemp); + Tip1 = POW(10, log10(TEF_Tmin) + (i+1)*TEF_dltemp); + Cool_fct(1.0, Ti, &emis, &LAMBDAT, cl_Z, cl_moli_mole, cl_mp); + h_SrcEC_TEF_lambda[i] = LAMBDAT*cl_mol/cl_moli_mole/cl_kB_mp; +# ifdef GAMER_DEBUG + if ( h_SrcEC_TEF_lambda[i] <= 0.0 ){ + Aux_Error( ERROR_INFO, "h_SrcEC_TEF_lambda[i] invalid (can not be smaller or equal to zero)!!\n" ); + } +# endif + h_SrcEC_TEF_alpha[i] = (log10(h_SrcEC_TEF_lambda[i+1]) - log10(h_SrcEC_TEF_lambda[i])) / (log10(Tip1) - log10(Ti)); + } + +// Initialize the constant of integration + double Ti_2, Tip1_2; + h_SrcEC_TEFc[TEF_N-1] = 0.0; // TEF(Tref) + for (int i=TEF_N-2; i>=0; i--){ + Ti_2 = POW(10.0, log10(TEF_Tmin) + i*TEF_dltemp); + Tip1_2 = POW(10.0, log10(TEF_Tmin) + (i+1)*TEF_dltemp); + if (h_SrcEC_TEF_alpha[i] != 1.0){ + h_SrcEC_TEFc[i] = h_SrcEC_TEFc[i+1] - (1.0/(1.0-h_SrcEC_TEF_alpha[i]))*(h_SrcEC_TEF_lambda[TEF_N-1]/h_SrcEC_TEF_lambda[i])*(Ti_2/TEF_TN)*(1.0-POW(Ti_2/Tip1_2, h_SrcEC_TEF_alpha[i]-1.0)); + } + else h_SrcEC_TEFc[i] = h_SrcEC_TEFc[i+1] - (h_SrcEC_TEF_lambda[TEF_N-1]/h_SrcEC_TEF_lambda[i])*(Ti_2/TEF_TN)*log(Ti_2/Tip1_2); + } + +# ifdef GPU + Src_PassData2GPU_ExactCooling(); +# endif + + +// set the auxiliary functions +// Src_WorkBeforeMajorFunc_EC_Ptr = Src_WorkBeforeMajorFunc_ExactCooling; +// Src_End_EC_Ptr = Src_End_ExactCooling; + +} // FUNCTION : Src_Init_ExactCooling + + + +// Sutherland-Dopita cooling function, with optimal parmetrization over a wide range of T and Z +void Cool_fct( double Dens, double Temp, double* Emis, double* Lambdat, double Z, double cl_moli_mole, double mp ){ + + double TLOGC = 5.65; + double QLOGC = -21.566; + double QLOGINFTY = -23.1; + double PPP = 0.8; + double TLOGM = 5.1; + double QLOGM = -20.85; + double SIG = 0.65; + double Zm = Z; + double TLOG = log10(Temp); + + *Lambdat = 0.; + if (Zm < 0) Zm = 0; + double QLOG0, ARG, BUMP1RHS, BUMP2LHS, QLAMBDA0, QLOG1, QLAMBDA1, ne_ni; + + if (TLOG >= 6.1) QLOG0 = -26.39 + 0.471*log10(Temp + 3.1623e6); + else if (TLOG >= 4.9){ + ARG = pow(10.0, (-(TLOG-4.9)/0.5)) + 0.077302; + QLOG0 = -22.16 + log10(ARG); + } + else if (TLOG >= 4.25){ + BUMP1RHS = -21.98 - ((TLOG-4.25)/0.55); + BUMP2LHS = -22.16 - pow((TLOG-4.9)/0.284, 2); + QLOG0 = fmax(BUMP1RHS, BUMP2LHS); + } + else QLOG0 = -21.98 - pow((TLOG-4.25)/0.2, 2); + + if (QLOG0 < -30.0) QLOG0 = -30.0; + QLAMBDA0 = pow(10.0, QLOG0); + + if (TLOG >= 5.65){ + QLOG1 = QLOGC - PPP*(TLOG-TLOGC); + QLOG1 = fmax(QLOG1, QLOGINFTY); + } + else QLOG1 = QLOGM - pow((TLOG-TLOGM)/SIG, 2); + + if (QLOG1 < -30.0) QLOG1 = -30.0; + QLAMBDA1 = pow(10.0, QLOG1); + + *Lambdat = (QLAMBDA0 + Zm*QLAMBDA1) / (UNIT_E*pow(UNIT_L, 3)/UNIT_T); +// For testing purpose (1) +// *Lambdat = 3.2217e-27 * sqrt(Temp) / (UNIT_E*pow(UNIT_L, 3)/UNIT_T); + +// For testing purpose (2) +// if (TLOG >= 5.0) *Lambdat = 3.2217e-27 * sqrt(Temp) / (UNIT_E*pow(UNIT_L, 3)/UNIT_T); +// else *Lambdat = 3.2217e-27 * pow(Temp, 0.4) / (UNIT_E*pow(UNIT_L, 3)/UNIT_T); + + ne_ni = (Dens*Dens) / (cl_moli_mole*mp*mp); + *Emis = ne_ni * (*Lambdat); // emissivity: lum/vol + +} + + + +//----------------------------------------------------------------------------------------- +// Function : Src_End_ExactCooling +// Description : Free the resources used by a user-specified source term +// +// Note : 1. Invoked by Src_End() +// --> Enable it by linking to the function pointer "Src_End_EC_Ptr" +// 2. Add "#ifndef __CUDACC__" since this routine is only useful on CPU +// +// Parameter : None +// +// Return : None +//----------------------------------------------------------------------------------------- +void Src_End_ExactCooling() +{ + +// free the memory of the h_SrcEC_* arrays + delete [] h_SrcEC_TEF_lambda; h_SrcEC_TEF_lambda = NULL; + delete [] h_SrcEC_TEF_alpha; h_SrcEC_TEF_alpha = NULL; + delete [] h_SrcEC_TEFc; h_SrcEC_TEFc = NULL; + + SrcTerms.EC_TEF_lambda_DevPtr = NULL; + SrcTerms.EC_TEF_alpha_DevPtr = NULL; + SrcTerms.EC_TEFc_DevPtr = NULL; + +# ifdef GPU + CUAPI_MemFree_ExactCooling(); +# endif + +} // FUNCTION : Src_End_ExactCooling + +#endif // #ifndef __CUDACC__ + + +#ifdef __CUDACC__ +//------------------------------------------------------------------------------------------------------- +// Function : CUAPI_MemFree_ExactCooling +// Description : Free the GPU memory of the ExactCooling arrays +// +// Note : 1. Invoked by Src_End_ExactCooling() +// +// Parameter : None +// +// Return : None +//------------------------------------------------------------------------------------------------------- +void CUAPI_MemFree_ExactCooling() +{ + + if ( d_SrcEC_TEF_lambda != NULL ) { CUDA_CHECK_ERROR( cudaFree( d_SrcEC_TEF_lambda ) ); d_SrcEC_TEF_lambda = NULL; } + if ( d_SrcEC_TEF_alpha != NULL ) { CUDA_CHECK_ERROR( cudaFree( d_SrcEC_TEF_alpha ) ); d_SrcEC_TEF_alpha = NULL; } + if ( d_SrcEC_TEFc != NULL ) { CUDA_CHECK_ERROR( cudaFree( d_SrcEC_TEFc ) ); d_SrcEC_TEFc = NULL; } + + SrcTerms.EC_TEF_lambda_DevPtr = NULL; + SrcTerms.EC_TEF_alpha_DevPtr = NULL; + SrcTerms.EC_TEFc_DevPtr = NULL; + +} // FUNCTION : CUAPI_MemFree_ExactCooling +#endif // #ifdef __CUDACC__ + + +#endif // #if ( MODEL == HYDRO ) diff --git a/src/SourceTerms/ExactCooling/CUSRC_Src_ExactCooling.cu b/src/SourceTerms/ExactCooling/CUSRC_Src_ExactCooling.cu new file mode 120000 index 0000000000..9e4142318b --- /dev/null +++ b/src/SourceTerms/ExactCooling/CUSRC_Src_ExactCooling.cu @@ -0,0 +1 @@ +CPU_Src_ExactCooling.cpp \ No newline at end of file diff --git a/src/SourceTerms/ExactCooling/dtSolver_ExactCooling.cpp b/src/SourceTerms/ExactCooling/dtSolver_ExactCooling.cpp new file mode 100644 index 0000000000..7c319b5927 --- /dev/null +++ b/src/SourceTerms/ExactCooling/dtSolver_ExactCooling.cpp @@ -0,0 +1,117 @@ +#include "GAMER.h" + + +extern bool IsInit_tcool[NLEVEL]; + + +//------------------------------------------------------------------------------------------------------- +// Function : Mis_GetTimeStep_ExactCooling +// Description : Estimate the evolution time-step constrained by the ExactCooling source term +// +// Note : 1. This function should be applied to both physical and comoving coordinates and always +// return the evolution time-step (dt) actually used in various solvers +// --> Physical coordinates : dt = physical time interval +// Comoving coordinates : dt = delta(scale_factor) / ( Hubble_parameter*scale_factor^3 ) +// --> We convert dt back to the physical time interval, which equals "delta(scale_factor)" +// in the comoving coordinates, in Mis_GetTimeStep() +// 2. Invoked by Mis_GetTimeStep() using the function pointer "Mis_GetTimeStep_User_Ptr", +// which must be set by a test problem initializer +// 3. Enabled by the runtime option "OPT__DT_USER" +// +// Parameter : lv : Target refinement level +// dTime_dt : dTime/dt (== 1.0 if COMOVING is off) +// +// Return : dt +//------------------------------------------------------------------------------------------------------- +double Mis_GetTimeStep_ExactCooling( const int lv, const double dTime_dt ) +{ + + if ( !SrcTerms.ExactCooling ) return HUGE_NUMBER; + + +// allocate memory for per-thread arrays +# ifdef OPENMP + const int NT = OMP_NTHREAD; // number of OpenMP threads +# else + const int NT = 1; +# endif + + double dt_EC = HUGE_NUMBER; + double *OMP_dt_EC = new double [NT]; + + +# pragma omp parallel + { +# ifdef OPENMP + const int TID = omp_get_thread_num(); +# else + const int TID = 0; +# endif + +// initialize the array + OMP_dt_EC[TID] = __DBL_MAX__; + + const double dh = amr->dh[lv]; + +# pragma omp for schedule( runtime ) + for (int PID=0; PIDNPatchComma[lv][1]; PID++) + { + for (int k=0; kMagSg[lv] ); +# else + real *B = NULL; +# endif // ifdef MHD ... else ... + + real tcool_Code = NULL_REAL; + + if ( IsInit_tcool[lv] ) { +// use the stored cooling time +# ifdef TCOOL + tcool_Code = amr->patch[ amr->FluSg[lv] ][lv][PID]->fluid[TCOOL][k][j][i]; +# endif + } + else { +// call Src_ExactCooling() to compute the cooling time if not initialized yet + const double z = amr->patch[0][lv][PID]->EdgeL[2] + (k+0.5)*dh; + const double y = amr->patch[0][lv][PID]->EdgeL[1] + (j+0.5)*dh; + const double x = amr->patch[0][lv][PID]->EdgeL[0] + (i+0.5)*dh; + +// get the input arrays + real fluid[FLU_NIN_S]; + for (int v=0; vpatch[ amr->FluSg[lv] ][lv][PID]->fluid[v][k][j][i]; + SrcTerms.EC_CPUPtr( fluid, B, &SrcTerms, 0.0, NULL_REAL, x, y, z, NULL_REAL, NULL_REAL, + MIN_DENS, MIN_PRES, MIN_EINT, NULL, + Src_EC_AuxArray_Flt, Src_EC_AuxArray_Int ); +# ifdef TCOOL + tcool_Code = fluid[TCOOL]; +# endif + } // if ( IsInit_tcool[lv] ) ... else ... + +// compare the cooling time and store the minimum value + OMP_dt_EC[TID] = FMIN( OMP_dt_EC[TID], tcool_Code ); + + }}} // i,j,k + } // for (int PID=0; PIDNPatchComma[lv][1]; PID++) + } // OpenMP parallel region + +// find the minimum over all OpenMP threads + for (int TID=0; TIDpatch[SaveSg_Flu][lv][PID]->fluid[0][0][0], h_Flu_Array_S_Out[N][0], FLU_NOUT_S*CUBE(PS1)*sizeof(real) ); + +# if ( defined GAMER_DEBUG && defined TCOOL ) + for (int k=0; kpatch[SaveSg_Flu][lv][PID]->fluid[TCOOL][k][j][i]; + + if ( !Aux_IsFinite(TCool) || TCool <= (real)0.0 ) + Aux_Error( ERROR_INFO, "Unphysical cooling time %14.7e (lv=%d, PID=%d, (i,j,k)=(%d,%d,%d)) !!\n", TCool, lv, PID, i, j, k ); + } +# endif } } // for (int TID=0; TID users may not define Src_WorkBeforeMajorFunc_User_Ptr if ( SrcTerms.User && Src_WorkBeforeMajorFunc_User_Ptr != NULL ) Src_WorkBeforeMajorFunc_User_Ptr ( lv, TimeNew, TimeOld, dt, diff --git a/src/TestProblem/Hydro/ClusterMerger/Flag_ClusterMerger.cpp b/src/TestProblem/Hydro/ClusterMerger/Flag_ClusterMerger.cpp new file mode 100644 index 0000000000..f5516dcfa5 --- /dev/null +++ b/src/TestProblem/Hydro/ClusterMerger/Flag_ClusterMerger.cpp @@ -0,0 +1,64 @@ +#include "GAMER.h" + +#if ( MODEL == HYDRO && defined GRAVITY ) + + +extern int Merger_Coll_NumHalos; +extern double R_acc; // the radius to compute the accretoin rate +extern double Jet_Vec[3][3]; // jet direction +extern double ClusterCen[3][3]; + +static bool FirstTime = true; + + +//------------------------------------------------------------------------------------------------------- +// Function : Flag_ClusterMerger +// Description : Flag cells for refinement for the cluster merger test problem +// +// Note : 1. Linked to the function pointer "Flag_User_Ptr" by Init_TestProb_Hydro_ClusterMerger() +// 2. Please turn on the runtime option "OPT__FLAG_USER" +// +// Parameter : i,j,k : Indices of the targeted element in the patch ptr[ amr->FluSg[lv] ][lv][PID] +// lv : Refinement level of the targeted patch +// PID : ID of the targeted patch +// Threshold : User-provided threshold for the flag operation, which is loaded from the +// file "Input__Flag_User" +// +// Return : "true" if the flag criteria are satisfied +// "false" if the flag criteria are not satisfied +//------------------------------------------------------------------------------------------------------- +bool Flag_ClusterMerger( const int i, const int j, const int k, const int lv, const int PID, const double *Threshold ) +{ + + const double dh = amr->dh[lv]; + const double Pos[3] = { amr->patch[0][lv][PID]->EdgeL[0] + (i+0.5)*dh, + amr->patch[0][lv][PID]->EdgeL[1] + (j+0.5)*dh, + amr->patch[0][lv][PID]->EdgeL[2] + (k+0.5)*dh }; + + bool Flag = false; + +// Flag cells within the target radius, and if the radius is not resolved with a specific number (Threshold[0]) of cells + for (int c=0; cdh[MAX_LEVEL]; + if ( R_acc/dh_max <= Threshold[0] ){ + Aux_Message( stderr, "WARNING : MAX_LEVEL is less than the desired refinement level set in Input__Flag_User!! dh_max = %13.7e\n", dh_max); + } + FirstTime = false; + } + + return Flag; + +} // FUNCTION : Flag_ClusterMerger + + + +#endif // #if ( MODEL == HYDRO && defined GRAVITY ) diff --git a/src/TestProblem/Hydro/ClusterMerger/Flu_ResetByUser_ClusterMerger.cpp b/src/TestProblem/Hydro/ClusterMerger/Flu_ResetByUser_ClusterMerger.cpp new file mode 100644 index 0000000000..b18bde4115 --- /dev/null +++ b/src/TestProblem/Hydro/ClusterMerger/Flu_ResetByUser_ClusterMerger.cpp @@ -0,0 +1,666 @@ +#include "GAMER.h" + +#if ( MODEL == HYDRO && defined GRAVITY ) + + +extern int Merger_Coll_NumHalos, Accretion_Mode; +extern double eta, eps_f, eps_m, R_acc, R_dep; // parameters of jet feedback +extern bool AGN_feedback; + +extern double CM_Bondi_SinkMass[3]; +extern double CM_Bondi_SinkMomX[3]; +extern double CM_Bondi_SinkMomY[3]; +extern double CM_Bondi_SinkMomZ[3]; +extern double CM_Bondi_SinkMomXAbs[3]; +extern double CM_Bondi_SinkMomYAbs[3]; +extern double CM_Bondi_SinkMomZAbs[3]; +extern double CM_Bondi_SinkE[3]; +extern double CM_Bondi_SinkEk[3]; +extern double CM_Bondi_SinkEt[3]; +extern int CM_Bondi_SinkNCell[3]; + +extern double Bondi_MassBH1; +extern double Bondi_MassBH2; +extern double Bondi_MassBH3; +extern double Mdot_BH1; // the accretion rate +extern double Mdot_BH2; +extern double Mdot_BH3; +extern double Jet_HalfHeight1; +extern double Jet_HalfHeight2; +extern double Jet_HalfHeight3; +extern double Jet_Radius1; +extern double Jet_Radius2; +extern double Jet_Radius3; +extern double Jet_Vec[3][3]; // jet direction +extern double Mdot[3]; // the feedback injection rate +extern double Pdot[3]; +extern double Edot[3]; +extern double GasVel[3][3]; // gas velocity +extern double SoundSpeed[3]; +extern double GasDens[3]; +extern double RelativeVel[3]; // the relative velocity between BH and gas +extern double ColdGasMass[3]; +extern double GasMass[3]; +extern double ParMass[3]; +extern double ClusterCen[3][3]; +extern double BH_Pos[3][3]; // BH position (for updating ClusterCen) +extern double BH_Vel[3][3]; // BH velocity + +double Jet_WaveK[3]; // jet wavenumber used in the sin() function to have smooth bidirectional jets +double Jet_HalfHeight[3]; +double Jet_Radius[3]; +double V_cyl[3]; // the volume of jet source +double M_inj[3], P_inj[3], E_inj[3]; // the injected density +double normalize_const[3]; // The exact normalization constant + +// the variables that need to be recorded +double E_inj_exp[3] = { 0.0, 0.0, 0.0 }; // the expected amount of injected energy +double M_inj_exp[3] = { 0.0, 0.0, 0.0 }; // the expected amount of injected gas mass +double dt_base; + +extern void GetClusterCenter( int lv, bool AdjustPos, bool AdjustVel, double Cen_old[][3], double Cen_new[][3], double Cen_Vel[][3] ); + +static bool FirstTime = true; +extern int JetDirection_NBin; // number of bins of the jet direction table +extern double *Time_table; // the time table of jet direction +extern double *Theta_table[3]; // the theta table of jet direction for 3 clusters +extern double *Phi_table[3]; // the phi table of jet direction for 3 clusters + +extern bool AdjustBHPos; +extern bool AdjustBHVel; +extern double AdjustPeriod; +extern int AdjustCount; // count the number of adjustments +extern int JetDirection_case; // Methods for choosing the jet direction +int merge_index = 0; // record BH 1 merge BH 2 / BH 2 merge BH 1 + +double ang_mom_sum[3][3] = { { 1.0, 0.0, 0.0 }, + { 1.0, 0.0, 0.0 }, + { 1.0, 0.0, 0.0 } }; +bool if_overlap = false; + +#ifdef MHD +extern double (*MHD_ResetByUser_VecPot_Ptr)( const double x, const double y, const double z, const double Time, + const double dt, const int lv, const char Component, double AuxArray[] ); +extern double (*MHD_ResetByUser_BField_Ptr)( const double x, const double y, const double z, const double Time, + const double dt, const int lv, const char Component, double AuxArray[], const double B_in, + const bool UseVecPot, const real *Ax, const real *Ay, const real *Az, + const int i, const int j, const int k ); +#endif + +//------------------------------------------------------------------------------------------------------- +// Function : Flu_ResetByUser_Func_ClusterMerger +// Description : Function to reset the fluid field in the Bondi accretion problem +// +// Note : 1. Invoked by Flu_ResetByUser_API_ClusterMerger() and Hydro_Init_ByFunction_AssignData() using the +// function pointer "Flu_ResetByUser_Func_Ptr" +// --> This function pointer is reset by Init_TestProb_Hydro_ClusterMerger() +// --> Hydro_Init_ByFunction_AssignData(): constructing initial condition +// Flu_ResetByUser_API_ClusterMerger() : after each update +// 2. Input fluid[] stores the original values +// 3. Even when DUAL_ENERGY is adopted, one does NOT need to set the dual-energy variable here +// --> It will be set automatically +// 4. Enabled by the runtime option "OPT__RESET_FLUID" +// +// Parameter : fluid : Fluid array storing both the input (origial) and reset values +// --> Including both active and passive variables +// x/y/z : Target physical coordinates +// Time : Target physical time +// dt : Time interval to advance solution +// lv : Target refinement level +// AuxArray : Auxiliary array +// +// Return : true : This cell has been reset +// false : This cell has not been reset +//------------------------------------------------------------------------------------------------------- +int Flu_ResetByUser_Func_ClusterMerger( real fluid[], const double Emag, const double x, const double y, const double z, + const double Time, const double dt, const int lv, double AuxArray[] ) +{ + const double Pos[3] = { x, y, z }; + +// (1) SMBH Accretion (note: gas depletion is currently disabled) + + double dr2[3][3], r2[3]; + const double V_dep = 4.0/3.0*M_PI*pow(R_dep,3.0); // the region to remove gas +// the density need to be removed + double D_dep[3] = { Mdot_BH1*dt/V_dep, Mdot_BH2*dt/V_dep, Mdot_BH3*dt/V_dep }; + int iftrue = 0; // mark whether this cell is reset or not [0/1] + + for (int c=0; c= Edot[1] ) status = 0; // only inject cluster 1 + else status = 1; // only inject cluster 2 + n_jet = Merger_Coll_NumHalos-1; + } + else { + status = 0; + n_jet = Merger_Coll_NumHalos; + } + + for (int c=status; c<(n_jet+status); c++) + { +// distance: jet center to mesh + for (int d=0; d<3; d++) Vec_c2m[c][d] = Pos[d] - ClusterCen[c][d]; + Dis_c2m = sqrt( SQR(Vec_c2m[c][0]) + SQR(Vec_c2m[c][1]) + SQR(Vec_c2m[c][2]) ); + +// vectors for calculating the distance between cells and the jet sources + for (int d=0; d<3; d++) TempVec[d] = ClusterCen[c][d] + Jet_Vec[c][d]; + +// distance: temporary vector to mesh + for (int d=0; d<3; d++) Vec_v2m[d] = Pos[d] - TempVec[d]; + Dis_v2m = sqrt( SQR(Vec_v2m[0]) + SQR(Vec_v2m[1]) + SQR(Vec_v2m[2]) ); + +// distance: jet center to temporary vector + Dis_c2v = sqrt( SQR(Jet_Vec[c][0]) + SQR(Jet_Vec[c][1]) + SQR(Jet_Vec[c][2]) ); + +// check whether or not the target cell is within the jet source + S = 0.5*( Dis_c2m + Dis_v2m + Dis_c2v ); + Area = sqrt( S*(S-Dis_c2m)*(S-Dis_v2m)*(S-Dis_c2v) ); + Jet_dr[c] = 2.0*Area/Dis_c2v; + Jet_dh[c] = sqrt( Dis_c2m*Dis_c2m - Jet_dr[c]*Jet_dr[c] ); + + if ( Jet_dh[c] <= Jet_HalfHeight[c] && Jet_dr[c] <= Jet_Radius[c] ) + { + which_cluster += c+1; + +// Record the old momentum + double MOMX_old = fluid[MOMX]; + double MOMY_old = fluid[MOMY]; + double MOMZ_old = fluid[MOMZ]; + + fluid[DENS] += M_inj[c]; + +// Transfer into BH frame + fluid[MOMX] -= BH_Vel[c][0]*fluid[DENS]; + fluid[MOMY] -= BH_Vel[c][1]*fluid[DENS]; + fluid[MOMZ] -= BH_Vel[c][2]*fluid[DENS]; + +// use a sine function to make the velocity smooth within the jet from +Jet_Vec to -Jet_Vec + EngySin = E_inj[c]*normalize_const[c]*sin( Jet_WaveK[c]*Jet_dh[c] ); + double P_SQR = SQR(fluid[MOMX])+SQR(fluid[MOMY])+SQR(fluid[MOMZ]); + double tmp_dens = fluid[DENS]-M_inj[c]; +// the new momentum is calculated from the old density, new density, old momentum and injected energy + double P_new = sqrt(2*fluid[DENS]*(EngySin+0.5*P_SQR/tmp_dens)); + P_new *= SIGN( Vec_c2m[c][0]*Jet_Vec[c][0] + Vec_c2m[c][1]*Jet_Vec[c][1] + Vec_c2m[c][2]*Jet_Vec[c][2] ); + fluid[MOMX] = P_new*Jet_Vec[c][0]; + fluid[MOMY] = P_new*Jet_Vec[c][1]; + fluid[MOMZ] = P_new*Jet_Vec[c][2]; + +// Transfer back into the rest frame + fluid[MOMX] += BH_Vel[c][0]*fluid[DENS]; + fluid[MOMY] += BH_Vel[c][1]*fluid[DENS]; + fluid[MOMZ] += BH_Vel[c][2]*fluid[DENS]; + + fluid[ENGY] += 0.5*((SQR(fluid[MOMX])+SQR(fluid[MOMY])+SQR(fluid[MOMZ]))/fluid[DENS]-(SQR(MOMX_old)+SQR(MOMY_old)+SQR(MOMZ_old))/tmp_dens); + } // if ( Jet_dh[c] <= Jet_HalfHeight[c] && Jet_dr[c] <= Jet_Radius[c] ) + } // for (int c=status; c<(n_jet+status); c++) + + if ( which_cluster >= 3 ) Aux_Error( ERROR_INFO, "Error: which_cluster >= 3!\n" ); + return which_cluster; + +} // FUNCTION : Flu_ResetByUser_Func_ClusterMerger + + + +//------------------------------------------------------------------------------------------------------- +// Function : Flu_ResetByUser_API_ClusterMerger +// Description : API for resetting the fluid array in the Bondi accretion problem +// +// Note : 1. Enabled by the runtime option "OPT__RESET_FLUID" +// 2. Invoked using the function pointer "Flu_ResetByUser_API_Ptr" +// --> This function pointer is reset by Init_TestProb_Hydro_ClusterMerger() +// 3. Currently does not work with "OPT__OVERLAP_MPI" +// 4. Invoke Flu_ResetByUser_Func_ClusterMerger() directly +// +// Parameter : lv : Target refinement level +// FluSg : Target fluid sandglass +// MagSg : Target B field sandglass +// TimeNew : Current physical time (system has been updated from TimeOld to TimeNew in EvolveLevel()) +// dt : Time interval to advance solution (can be different from TimeNew-TimeOld in COMOVING) +//------------------------------------------------------------------------------------------------------- +void Flu_ResetByUser_API_ClusterMerger( const int lv, const int FluSg, const int MagSg, const double TimeNew, const double dt ) +{ + + const bool CurrentMaxLv = ( NPatchTotal[lv] > 0 && ( lv == MAX_LEVEL || NPatchTotal[lv+1] == 0 ) ); + const double dh = amr->dh[lv]; + const real dv = CUBE(dh); +# if ( MODEL == HYDRO || MODEL == MHD ) + const real Gamma_m1 = GAMMA - (real)1.0; + const real _Gamma_m1 = (real)1.0 / Gamma_m1; +# endif + +// (1) Get the BH position and velocity and adjust them if needed + bool AdjustPosNow = false; + bool AdjustVelNow = false; + if ( CurrentMaxLv && AdjustCount < int(TimeNew/AdjustPeriod)){ // only adjust the BHs on the current maximum level + if ( AdjustBHPos == true ) AdjustPosNow = true; + if ( AdjustBHVel == true ) AdjustVelNow = true; + AdjustCount += 1; + } + GetClusterCenter( lv, AdjustPosNow, AdjustVelNow, BH_Pos, ClusterCen, BH_Vel ); + + +// (2) Decide whether to merge BHs + double RelativeBHPos[3] = { BH_Pos[0][0]-BH_Pos[1][0], BH_Pos[0][1]-BH_Pos[1][1], BH_Pos[0][2]-BH_Pos[1][2] }; + double RelativeBHVel[3] = { BH_Vel[0][0]-BH_Vel[1][0], BH_Vel[0][1]-BH_Vel[1][1], BH_Vel[0][2]-BH_Vel[1][2] }; + double AbsRelPos = sqrt( SQR(RelativeBHPos[0])+SQR(RelativeBHPos[1])+SQR(RelativeBHPos[2]) ); + double AbsRelVel = sqrt( SQR(RelativeBHVel[0])+SQR(RelativeBHVel[1])+SQR(RelativeBHVel[2]) ); + double escape_vel[2]; + double soften = amr->dh[MAX_LEVEL]; + if ( AbsRelPos > soften ){ + escape_vel[0] = sqrt(2*NEWTON_G*Bondi_MassBH2/AbsRelPos); + escape_vel[1] = sqrt(2*NEWTON_G*Bondi_MassBH1/AbsRelPos); + } + else{ + escape_vel[0] = sqrt(2*NEWTON_G*Bondi_MassBH2/soften); + escape_vel[1] = sqrt(2*NEWTON_G*Bondi_MassBH1/soften); + } + +// Merge the two BHs if they are located within R_acc, and the relative velocity is small enough + if ( Merger_Coll_NumHalos == 2 ){ + if ( AbsRelPos < R_acc && ( AbsRelVel < 3*escape_vel[1] || AbsRelVel < 3*escape_vel[0] ) ){ + Merger_Coll_NumHalos -= 1; + if ( Bondi_MassBH1 >= Bondi_MassBH2 ) merge_index = 1; // record BH 1 merge BH 2 / BH 2 merge BH 1 + else merge_index = 2; + Bondi_MassBH1 += Bondi_MassBH2; + Bondi_MassBH2 = 0.0; +// Relabel the BH and DM particles being merged + for (long p=0; pPar->NPar_AcPlusInac; p++) { + bool if_cluster = false; + if ( amr->Par->Type[p] == real(PTYPE_CLUSTER+1) || amr->Par->Type[p] == real(PTYPE_CEN+1) ){ + if_cluster = true; + } + if ( amr->Par->Mass[p] >= (real)0.0 && if_cluster ){ + amr->Par->Type[p] = PTYPE_CLUSTER; + } + } + Aux_Message( stdout, "BHs Merge! In rank %d, TimeNew = %14.8e; merge_index = %d, BHPos1 = %14.8e, %14.8e, %14.8e; BHPos2 = %14.8e, %14.8e, %14.8e; BHVel1 = %14.8e, %14.8e, %14.8e; BHVel2 = %14.8e, %14.8e, %14.8e; AbsRelPos = %14.8e, AbsRelVel = %14.8e, escape_vel[0] = %14.8e, escape_vel[1] = %14.8e.\n", MPI_Rank, TimeNew, merge_index, BH_Pos[0][0], BH_Pos[0][1], BH_Pos[0][2], BH_Pos[1][0], BH_Pos[1][1], BH_Pos[1][2], BH_Vel[0][0], BH_Vel[0][1], BH_Vel[0][2], BH_Vel[1][0], BH_Vel[1][1], BH_Vel[1][2], AbsRelPos, AbsRelVel, escape_vel[0], escape_vel[1]); + } + } + + + if ( AGN_feedback ) { +// (3) Set the injection parameters + double Mdot_BH[3] = { Mdot_BH1, Mdot_BH2, Mdot_BH3 }; + double Bondi_MassBH[3] = { Bondi_MassBH1, Bondi_MassBH2, Bondi_MassBH3 }; + + Jet_HalfHeight[0] = Jet_HalfHeight1; + Jet_HalfHeight[1] = Jet_HalfHeight2; + Jet_HalfHeight[2] = Jet_HalfHeight3; + Jet_Radius[0] = Jet_Radius1; + Jet_Radius[1] = Jet_Radius2; + Jet_Radius[2] = Jet_Radius3; + +// Set the jet direction vector + if ( JetDirection_case == 1 ) { // Fixed at x-axis + for (int c=0; cNPatchComma[lv][1]; PID++) + { + x02 = amr->patch[0][lv][PID]->EdgeL[0] + 0.5*dh; + y02 = amr->patch[0][lv][PID]->EdgeL[1] + 0.5*dh; + z02 = amr->patch[0][lv][PID]->EdgeL[2] + 0.5*dh; + + for (int k=0; kpatch[FluSg][lv][PID]->fluid[v][k][j][i]; + } + +# ifdef MHD + const real Emag = MHD_GetCellCenteredBEnergyInPatch( lv, PID, i, j, k, MagSg ); +# else + const real Emag = (real)0.0; +# endif + + int status_overlap = 0; + + for (int c=0; c= 3 ) if_overlap_each_rank = true; + + }}} + } // for (int PID=0; PIDNPatchComma[lv][1]; PID++) + + for (int c=0; cNPatchComma[lv][1]; PID++) { + const double *EdgeL = amr->patch[0][lv][PID]->EdgeL; + const double *EdgeR = amr->patch[0][lv][PID]->EdgeR; + const double patch_pos[3] = { (EdgeL[0]+EdgeR[0])*0.5, (EdgeL[1]+EdgeR[1])*0.5, (EdgeL[2]+EdgeR[2])*0.5 }; + const double patch_d = sqrt(SQR(EdgeL[0]-EdgeR[0])+SQR(EdgeL[1]-EdgeR[1])+SQR(EdgeL[2]-EdgeR[2]))*0.5; + if (SQR(patch_pos[0]-ClusterCen[c][0])+SQR(patch_pos[1]-ClusterCen[c][1])+SQR(patch_pos[2]-ClusterCen[c][2]) <= SQR(2*R_acc+patch_d)){ + for (int p=0; ppatch[0][lv][PID]->NPar; p++) { + const long ParID = amr->patch[0][lv][PID]->ParList[p]; + const real ParX = amr->Par->PosX[ParID]; + const real ParY = amr->Par->PosY[ParID]; + const real ParZ = amr->Par->PosZ[ParID]; + const real ParM = amr->Par->Mass[ParID]; + if ( SQR(ParX-ClusterCen[c][0])+SQR(ParY-ClusterCen[c][1])+SQR(ParZ-ClusterCen[c][2]) <= SQR(R_acc) ){ + par_mass += ParM; + } } } } + MPI_Allreduce( &par_mass, &ParMass[c], 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD ); + } + } // for (int c=0; cNPatchComma[lv][1]; PID++) + { + x0 = amr->patch[0][lv][PID]->EdgeL[0] + 0.5*dh; + y0 = amr->patch[0][lv][PID]->EdgeL[1] + 0.5*dh; + z0 = amr->patch[0][lv][PID]->EdgeL[2] + 0.5*dh; + + for (int k=0; kpatch[FluSg][lv][PID]->fluid[v][k][j][i]; + +// backup the unmodified values since we want to record the amount of sunk variables removed at the maximum level + fluid_bk[v] = fluid[v]; + } + +# ifdef MHD + const real Emag = MHD_GetCellCenteredBEnergyInPatch( lv, PID, i, j, k, MagSg ); +# else + const real Emag = (real)0.0; +# endif + +// reset this cell + if (CurrentMaxLv) Reset = Flu_ResetByUser_Func_ClusterMerger( fluid, Emag, x, y, z, TimeNew, dt, lv, NULL ); + +// operations necessary only when this cell has been reset + if ( Reset != 0 ) + { +# if ( MODEL == HYDRO || MODEL == MHD ) +// abort the program instead of applying a density floor + if ( fluid[DENS] < MIN_DENS ) { + Aux_Error( ERROR_INFO, "Error: Fluid density has reached the floor!\n" ); + } +// apply an energy floor + fluid[ENGY] = Hydro_CheckMinEintInEngy( fluid[DENS], fluid[MOMX], fluid[MOMY], fluid[MOMZ], fluid[ENGY], + (real)MIN_EINT, 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 ); +# endif + +// floor and normalize passive scalars +# if ( NCOMP_PASSIVE > 0 ) + for (int v=NCOMP_FLUID; vpatch[FluSg][lv][PID]->fluid[v][k][j][i] = fluid[v]; + + }}} // i,j,k + } // for (int PID=0; PIDNPatchComma[lv][1]; PID++) + } // if ( AGN_feedback ) + +} // FUNCTION : Flu_ResetByUser_API_ClusterMerger + + +#endif // #if ( MODEL == HYDRO && defined GRAVITY ) diff --git a/src/TestProblem/Hydro/ClusterMerger/Init_TestProb_ClusterMerger.cpp b/src/TestProblem/Hydro/ClusterMerger/Init_TestProb_ClusterMerger.cpp index ab90f7f51b..b7a127c673 100644 --- a/src/TestProblem/Hydro/ClusterMerger/Init_TestProb_ClusterMerger.cpp +++ b/src/TestProblem/Hydro/ClusterMerger/Init_TestProb_ClusterMerger.cpp @@ -10,12 +10,22 @@ // problem-specific global variables // ======================================================================================= int Merger_Coll_NumHalos; // number of clusters + bool AGN_feedback; // turn on/off (1/0) AGN feedback static char Merger_File_Prof1[1000]; // profile table of cluster 1 static char Merger_File_Prof2[1000]; // profile table of cluster 2 static char Merger_File_Prof3[1000]; // profile table of cluster 3 char Merger_File_Par1 [1000]; // particle file of cluster 1 char Merger_File_Par2 [1000]; // particle file of cluster 2 char Merger_File_Par3 [1000]; // particle file of cluster 3 + double Unit_R1; // unit of length in cluster 1's IC file + double Unit_D1; // unit of density in cluster 1's IC file + double Unit_P1; // unit of pressure in cluster 1's IC file + double Unit_R2; // unit of length in cluster 2's IC file + double Unit_D2; // unit of density in cluster 2's IC file + double Unit_P2; // unit of pressure in cluster 2's IC file + double Unit_R3; // unit of length in cluster 3's IC file + double Unit_D3; // unit of density in cluster 3's IC file + double Unit_P3; // unit of pressure in cluster 3's IC file bool Merger_Coll_IsGas1; // (true/false) --> does cluster 1 have gas bool Merger_Coll_IsGas2; // (true/false) --> does cluster 2 have gas bool Merger_Coll_IsGas3; // (true/false) --> does cluster 3 have gas @@ -59,6 +69,64 @@ static FieldIdx_t ColorField1Idx = Idx_Undefined; static FieldIdx_t ColorField2Idx = Idx_Undefined; static FieldIdx_t ColorField3Idx = Idx_Undefined; + int Accretion_Mode; // 1: hot mode; 2: code mode; 3: combine (hot + cold) + double eta; // mass loading factor in jet feedback + double eps_f; // the radiative efficiency in jet feedback + double eps_m; // the fraction of total energy that goes into the thermal energy in jet feedback + double R_acc; // accretion radius: compute the accretion rate + double R_dep; // radius to deplete the accreted gas + + double CM_Bondi_SinkMass[3]; // total mass change in the feedback region in one global time-step + double CM_Bondi_SinkMomX[3]; // total x-momentum change ... + double CM_Bondi_SinkMomY[3]; // total y-momentum change ... + double CM_Bondi_SinkMomZ[3]; // total z-momentum change ... + double CM_Bondi_SinkMomXAbs[3]; // total |x-momentum| change ... + double CM_Bondi_SinkMomYAbs[3]; // total |y-momentum| change ... + double CM_Bondi_SinkMomZAbs[3]; // total |z-momentum| change ... + double CM_Bondi_SinkE[3]; // total injected energy ... + double CM_Bondi_SinkEk[3]; // total injected kinetic energy ... + double CM_Bondi_SinkEt[3]; // total injected thermal energy ... + int CM_Bondi_SinkNCell[3]; // total number of finest cells within the feedback region + + double Bondi_MassBH1; // black hole mass of cluster 1 + double Bondi_MassBH2; // black hole mass of cluster 2 + double Bondi_MassBH3; // black hole mass of cluster 3 + double Mdot_BH1; // the accretion rate of BH 1 + double Mdot_BH2; // the accretion rate of BH 2 + double Mdot_BH3; // the accretion rate of BH 3 + double Jet_HalfHeight1; // half height of the cylinder-shape jet source of cluster 1 + double Jet_HalfHeight2; // half height of the cylinder-shape jet source of cluster 2 + double Jet_HalfHeight3; // half height of the cylinder-shape jet source of cluster 3 + double Jet_Radius1; // radius of the cylinder-shape jet source of cluster 1 + double Jet_Radius2; // radius of the cylinder-shape jet source of cluster 2 + double Jet_Radius3; // radius of the cylinder-shape jet source of cluster 3 + double Mdot[3]; // the feedback injeciton rate of mass + double Pdot[3]; // the feedback injeciton rate of momentum + double Edot[3]; // the feedback injeciton rate of total energy + double Jet_Vec[3][3]; // jet direction + double GasVel[3][3]; // average gas velocity inside the accretion radius + double SoundSpeed[3]; // average sound speed inside the accreiton radius + double GasDens[3]; // average gas density inside the accreiton radius + double RelativeVel[3]; // relative velocity between BH and gas for each cluster + double ColdGasMass[3]; // cold gas mass inside the accretion radius + double GasMass[3]; // total gas mass inside the accretion radius + double ParMass[3]; // total DM mass inside the accretion radius + double ClusterCen[3][3]; // the center of each cluster + double BH_Pos[3][3]; // BH position of each cluster + double BH_Vel[3][3]; // BH velocity of each cluster + + int JetDirection_NBin; // number of bins of the jet direction table +static double *JetDirection = NULL; // jet direction[time/theta_1/phi_1/theta_2/phi_2/theta_3/phi_3] + double *Time_table; // the time table of jet direction + double *Theta_table[3]; // the theta table of jet direction for 3 clusters + double *Phi_table[3]; // the phi table of jet direction for 3 clusters + + bool AdjustBHPos; // (true/false) --> Adjust the BH position + bool AdjustBHVel; // (true/false) --> Adjust the BH velocity + double AdjustPeriod; // the time interval of adjustment + int AdjustCount = 0; // count the number of adjustments + int JetDirection_case; // Methods for choosing the jet direction: 1. Fixed at x-axis; 2. Import from table (generate JetDirection.txt); 3. Align with angular momentum + bool fixBH; // fix the BH at the simulation box center and set its velocity to be zero (1 cluster only) // ======================================================================================= // problem-specific function prototypes @@ -72,11 +140,20 @@ void Par_Init_ByFunction_ClusterMerger(const long NPar_ThisRank, void Aux_Record_ClusterMerger(); #endif +bool Flag_ClusterMerger( const int i, const int j, const int k, const int lv, const int PID, const double *Threshold ); int Read_Num_Points_ClusterMerger(std::string filename); void Read_Profile_ClusterMerger(std::string filename, std::string fieldname, double field[]); void AddNewField_ClusterMerger(); +int Flu_ResetByUser_Func_ClusterMerger( real fluid[], const double Emag, const double x, const double y, const double z, + const double Time, const double dt, const int lv, double AuxArray[] ); +void Flu_ResetByUser_API_ClusterMerger( const int lv, const int FluSg, const int MagSg, const double TimeNew, const double dt ); + +extern void (*Flu_ResetByUser_API_Ptr)( const int lv, const int FluSg, const int MagSg, const double TimeNew, const double dt ); +void Output_ClusterMerger(); +void Init_User_ClusterMerger(); + //------------------------------------------------------------------------------------------------------- // Function : Validate // Description : Validate the compilation flags and runtime parameters for this test problem @@ -169,6 +246,7 @@ void SetParameter() // ReadPara->Add( "KEY_IN_THE_FILE", &VARIABLE, DEFAULT, MIN, MAX ); // ******************************************************************************************************************************** ReadPara->Add( "Merger_Coll_NumHalos", &Merger_Coll_NumHalos, 2, 1, 3 ); + ReadPara->Add( "AGN_feedback", &AGN_feedback, false, Useless_bool, Useless_bool ); ReadPara->Add( "Merger_Coll_IsGas1", &Merger_Coll_IsGas1, true, Useless_bool, Useless_bool ); ReadPara->Add( "Merger_Coll_IsGas2", &Merger_Coll_IsGas2, true, Useless_bool, Useless_bool ); ReadPara->Add( "Merger_Coll_IsGas3", &Merger_Coll_IsGas3, true, Useless_bool, Useless_bool ); @@ -178,6 +256,15 @@ void SetParameter() ReadPara->Add( "Merger_File_Par2", Merger_File_Par2, NoDef_str, Useless_str, Useless_str ); ReadPara->Add( "Merger_File_Prof3", Merger_File_Prof3, NoDef_str, Useless_str, Useless_str ); ReadPara->Add( "Merger_File_Par3", Merger_File_Par3, NoDef_str, Useless_str, Useless_str ); + ReadPara->Add( "Unit_R1", &Unit_R1, 1.0, NoMin_double, NoMax_double ); + ReadPara->Add( "Unit_D1", &Unit_D1, 1.0, NoMin_double, NoMax_double ); + ReadPara->Add( "Unit_P1", &Unit_P1, 1.0, NoMin_double, NoMax_double ); + ReadPara->Add( "Unit_R2", &Unit_R2, 1.0, NoMin_double, NoMax_double ); + ReadPara->Add( "Unit_D2", &Unit_D2, 1.0, NoMin_double, NoMax_double ); + ReadPara->Add( "Unit_P2", &Unit_P2, 1.0, NoMin_double, NoMax_double ); + ReadPara->Add( "Unit_R3", &Unit_R3, 1.0, NoMin_double, NoMax_double ); + ReadPara->Add( "Unit_D3", &Unit_D3, 1.0, NoMin_double, NoMax_double ); + ReadPara->Add( "Unit_P3", &Unit_P3, 1.0, NoMin_double, NoMax_double ); ReadPara->Add( "Merger_Coll_PosX1", &Merger_Coll_PosX1, -1.0, NoMin_double, NoMax_double ); ReadPara->Add( "Merger_Coll_PosY1", &Merger_Coll_PosY1, -1.0, NoMin_double, NoMax_double ); ReadPara->Add( "Merger_Coll_PosX2", &Merger_Coll_PosX2, -1.0, NoMin_double, NoMax_double ); @@ -192,6 +279,29 @@ void SetParameter() ReadPara->Add( "Merger_Coll_VelY3", &Merger_Coll_VelY3, -1.0, NoMin_double, NoMax_double ); ReadPara->Add( "Merger_Coll_UseMetals", &Merger_Coll_UseMetals, true, Useless_bool, Useless_bool ); ReadPara->Add( "Merger_Coll_LabelCenter",&Merger_Coll_LabelCenter,true, Useless_bool, Useless_bool ); + ReadPara->Add( "Bondi_MassBH1", &Bondi_MassBH1, -1.0, Eps_double, NoMax_double ); + ReadPara->Add( "Bondi_MassBH2", &Bondi_MassBH2, -1.0, Eps_double, NoMax_double ); + ReadPara->Add( "Bondi_MassBH3", &Bondi_MassBH3, -1.0, Eps_double, NoMax_double ); + ReadPara->Add( "R_acc", &R_acc, -1.0, NoMin_double, NoMax_double ); + ReadPara->Add( "R_dep", &R_dep, -1.0, NoMin_double, NoMax_double ); + ReadPara->Add( "Mdot_BH1", &Mdot_BH1, -1.0, NoMin_double, NoMax_double ); + ReadPara->Add( "Mdot_BH2", &Mdot_BH2, -1.0, NoMin_double, NoMax_double ); + ReadPara->Add( "Mdot_BH3", &Mdot_BH3, -1.0, NoMin_double, NoMax_double ); + ReadPara->Add( "Jet_HalfHeight1", &Jet_HalfHeight1, -1.0, Eps_double, NoMax_double ); + ReadPara->Add( "Jet_HalfHeight2", &Jet_HalfHeight2, -1.0, Eps_double, NoMax_double ); + ReadPara->Add( "Jet_HalfHeight3", &Jet_HalfHeight3, -1.0, Eps_double, NoMax_double ); + ReadPara->Add( "Jet_Radius1", &Jet_Radius1, -1.0, Eps_double, NoMax_double ); + ReadPara->Add( "Jet_Radius2", &Jet_Radius2, -1.0, Eps_double, NoMax_double ); + ReadPara->Add( "Jet_Radius3", &Jet_Radius3, -1.0, Eps_double, NoMax_double ); + ReadPara->Add( "Accretion_Mode", &Accretion_Mode, 1, 1, 3 ); + ReadPara->Add( "eta", &eta, -1.0, NoMin_double, NoMax_double ); + ReadPara->Add( "eps_f", &eps_f, -1.0, NoMin_double, NoMax_double ); + ReadPara->Add( "eps_m", &eps_m, -1.0, NoMin_double, NoMax_double ); + ReadPara->Add( "AdjustBHPos", &AdjustBHPos, false, Useless_bool, Useless_bool ); + ReadPara->Add( "AdjustBHVel", &AdjustBHVel, false, Useless_bool, Useless_bool ); + ReadPara->Add( "AdjustPeriod", &AdjustPeriod, -1.0, NoMin_double, NoMax_double ); + ReadPara->Add( "JetDirection_case", &JetDirection_case, 1, 1, 3 ); + ReadPara->Add( "fixBH", &fixBH, false, Useless_bool, Useless_bool ); ReadPara->Read( FileName ); @@ -202,6 +312,16 @@ void SetParameter() Aux_Error( ERROR_INFO, "please set NCOMP_PASSIVE_USER (currently %d) == Merger_Coll_NumHalos + Merger_Coll_UseMetals (currently %d) in the Makefile !!\n", NCOMP_PASSIVE_USER, Merger_Coll_NumHalos + (int)Merger_Coll_UseMetals ); +// Set the correct parameters when fixing the BH + if ( Merger_Coll_NumHalos != 1 && fixBH == true ){ + fixBH = false; + Aux_Message( stdout, "WARNING! Reset fixBH to be false for multiple clusters!\n" ); + } + if ( fixBH == true ){ + AdjustBHPos = false; + AdjustBHVel = false; + } + // convert to code units Merger_Coll_PosX1 *= Const_kpc / UNIT_L; Merger_Coll_PosY1 *= Const_kpc / UNIT_L; @@ -215,6 +335,21 @@ void SetParameter() Merger_Coll_VelY2 *= (Const_km/Const_s) / UNIT_V; Merger_Coll_VelX3 *= (Const_km/Const_s) / UNIT_V; Merger_Coll_VelY3 *= (Const_km/Const_s) / UNIT_V; + Bondi_MassBH1 *= Const_Msun/UNIT_M; + Bondi_MassBH2 *= Const_Msun/UNIT_M; + Bondi_MassBH3 *= Const_Msun/UNIT_M; + R_acc *= Const_kpc / UNIT_L; + R_dep *= Const_kpc / UNIT_L; + Mdot_BH1 *= (Const_Msun/Const_yr) / (UNIT_M / UNIT_T); + Mdot_BH2 *= (Const_Msun/Const_yr) / (UNIT_M / UNIT_T); + Mdot_BH3 *= (Const_Msun/Const_yr) / (UNIT_M / UNIT_T); + Jet_HalfHeight1 *= Const_kpc / UNIT_L; + Jet_HalfHeight2 *= Const_kpc / UNIT_L; + Jet_HalfHeight3 *= Const_kpc / UNIT_L; + Jet_Radius1 *= Const_kpc / UNIT_L; + Jet_Radius2 *= Const_kpc / UNIT_L; + Jet_Radius3 *= Const_kpc / UNIT_L; + AdjustPeriod *= Const_Myr / UNIT_T; if ( OPT__INIT != INIT_BY_RESTART ) { @@ -246,12 +381,11 @@ void SetParameter() Read_Profile_ClusterMerger(filename1, "/fields/metallicity", Table_M1); else for ( int i; i < Merger_NBin1; i++ ) Table_M1[i] = 0.0; - // convert to code units (assuming the input units are cgs) for ( int b=0; b 2 && Merger_Coll_IsGas3 ) +// (2-2) Initialize the BH position and velocity + double ClusterCenter[3][3] = {{ Merger_Coll_PosX1, Merger_Coll_PosY1, amr->BoxCenter[2] }, + { Merger_Coll_PosX2, Merger_Coll_PosY2, amr->BoxCenter[2] }, + { Merger_Coll_PosX3, Merger_Coll_PosY3, amr->BoxCenter[2] }}; + double CenterVel[3][3] = {{ Merger_Coll_VelX1, Merger_Coll_VelY1, 0.0 }, + { Merger_Coll_VelX2, Merger_Coll_VelY2, 0.0 }, + { Merger_Coll_VelX3, Merger_Coll_VelY3, 0.0 }}; + + for (int c=0; c a helper macro PRINT_WARNING is defined in TestProb.h const long End_Step_Default = __INT_MAX__; const double End_T_Default = 10.0*Const_Gyr/UNIT_T; @@ -412,12 +577,13 @@ void SetParameter() } -// (4) make a note +// (6) make a note if ( MPI_Rank == 0 ) { Aux_Message( stdout, "=============================================================================\n" ); Aux_Message( stdout, " test problem ID = %d\n", TESTPROB_ID ); Aux_Message( stdout, " number of clusters = %d\n", Merger_Coll_NumHalos ); + Aux_Message( stdout, " turn on AGN feedback = %s\n", (AGN_feedback)? "yes":"no" ); if ( Merger_Coll_IsGas1 ) Aux_Message( stdout, " profile file 1 = %s\n", Merger_File_Prof1 ); Aux_Message( stdout, " particle file 1 = %s\n", Merger_File_Par1 ); @@ -452,6 +618,18 @@ void SetParameter() Aux_Message( stdout, " use metals = %s\n", (Merger_Coll_UseMetals)? "yes":"no" ); Aux_Message( stdout, " label cluster centers = %s\n", (Merger_Coll_LabelCenter)? "yes":"no" ); Aux_Message( stdout, "=============================================================================\n" ); + +// Check if the accretion region is larger than the jet cylinder + if ( R_acc < Jet_HalfHeight1 ) Aux_Message( stderr, "WARNING : R_acc is less than Jet_HalfHeight1!!\n"); + if ( R_acc < Jet_Radius1 ) Aux_Message( stderr, "WARNING : R_acc is less than Jet_Radius1!!\n"); + if ( Merger_Coll_NumHalos > 1 ) { + if ( R_acc < Jet_HalfHeight2 ) Aux_Message( stderr, "WARNING : R_acc is less than Jet_HalfHeight2!!\n"); + if ( R_acc < Jet_Radius2 ) Aux_Message( stderr, "WARNING : R_acc is less than Jet_Radius2!!\n"); + } + if ( Merger_Coll_NumHalos > 2 ) { + if ( R_acc < Jet_HalfHeight3 ) Aux_Message( stderr, "WARNING : R_acc is less than Jet_HalfHeight3!!\n"); + if ( R_acc < Jet_Radius3 ) Aux_Message( stderr, "WARNING : R_acc is less than Jet_Radius3!!\n"); + } } if ( MPI_Rank == 0 ) Aux_Message( stdout, " Setting runtime parameters ... done\n" ); @@ -694,11 +872,17 @@ void Init_TestProb_Hydro_ClusterMerger() // set the function pointers of various problem-specific routines Init_Function_User_Ptr = SetGridIC; + Flag_User_Ptr = Flag_ClusterMerger; End_User_Ptr = End_ClusterMerger; Aux_Record_User_Ptr = Aux_Record_ClusterMerger; Par_Init_ByFunction_Ptr = Par_Init_ByFunction_ClusterMerger; Init_Field_User_Ptr = AddNewField_ClusterMerger; - //Par_Init_Attribute_User_Ptr = AddNewParticleAttribute_ClusterMerger; + + Flu_ResetByUser_Func_Ptr = Flu_ResetByUser_Func_ClusterMerger; + Flu_ResetByUser_API_Ptr = Flu_ResetByUser_API_ClusterMerger; + Output_User_Ptr = Output_ClusterMerger; + Init_User_Ptr = Init_User_ClusterMerger; + # ifdef MHD Init_Function_BField_User_Ptr = SetBFieldIC; # endif @@ -792,3 +976,129 @@ void AddNewField_ClusterMerger() } #endif + + +// Restart only: initialize the BH variables (position, velocity, mass) from particles labelled with PTYPE_CEN +void Init_User_ClusterMerger() +{ + if ( OPT__INIT == INIT_BY_RESTART ){ + + if ( OPT__RESTART_RESET ) { + printf("Error! OPT__RESTART_RESET should be disabled.\n"); + return; + } + + const char FileName[] = "BH_variable.bin"; + int TargetDumpID = DumpID-1; //INIT_DUMPID; + + FILE* File_User = fopen(FileName, "rb"); + if ( File_User == NULL ) { + Aux_Error( ERROR_INFO, "Error opening the file \"%s\"\n", FileName ); + return; + } + + if ( MPI_Rank == 0 ) { + + double BH_Mass[3] = {0.0, 0.0, 0.0}; + + double Time; + int dumpID; + while (fread(&Time, sizeof(double), 1, File_User) == 1) { + fread(&dumpID, sizeof(int), 1, File_User); + printf("dumpID = %d, TargetDumpID = %d\n", dumpID, TargetDumpID); + if (dumpID == TargetDumpID) { + fread(&Merger_Coll_NumHalos, sizeof(int), 1, File_User); + for (int c=0; c 0 && ( lv == MAX_LEVEL || NPatchTotal[lv+1] == 0 ) ); + +// Initialize min_pos to be the old center + for (int c=0; cPar->PosX, amr->Par->PosY, amr->Par->PosZ }; + for (int c=0; cNPatchComma[lv][1]; PID++) { + const double *EdgeL = amr->patch[0][lv][PID]->EdgeL; + const double *EdgeR = amr->patch[0][lv][PID]->EdgeR; + const double patch_pos[3] = { (EdgeL[0]+EdgeR[0])*0.5, (EdgeL[1]+EdgeR[1])*0.5, (EdgeL[2]+EdgeR[2])*0.5 }; + const double patch_d = sqrt(SQR(EdgeL[0]-EdgeR[0])+SQR(EdgeL[1]-EdgeR[1])+SQR(EdgeL[2]-EdgeR[2]))*0.5; + if (SQR(patch_pos[0]-Cen_new_pre[c][0])+SQR(patch_pos[1]-Cen_new_pre[c][1])+SQR(patch_pos[2]-Cen_new_pre[c][2]) <= SQR(20*R_acc+patch_d)){ + for (int p=0; ppatch[0][lv][PID]->NPar; p++) { + const long ParID = amr->patch[0][lv][PID]->ParList[p]; + const real ParX_tmp = amr->Par->PosX[ParID]; + const real ParY_tmp = amr->Par->PosY[ParID]; + const real ParZ_tmp = amr->Par->PosZ[ParID]; + const real ParM_tmp = amr->Par->Mass[ParID]; + const real VelX_tmp = amr->Par->VelX[ParID]; + const real VelY_tmp = amr->Par->VelY[ParID]; + const real VelZ_tmp = amr->Par->VelZ[ParID]; + bool if_cluster = false; + if ( amr->Par->Type[ParID] == real(PTYPE_CLUSTER+c) || amr->Par->Type[ParID] == real(PTYPE_CEN+c) ){ + if_cluster = true; + } + if ( if_cluster && SQR(ParX_tmp-Cen_new_pre[c][0])+SQR(ParY_tmp-Cen_new_pre[c][1])+SQR(ParZ_tmp-Cen_new_pre[c][2]) <= SQR(10*R_acc) ){ +// Record the mass, position and velocity of this particle + ParX[c][num_par[c]] = ParX_tmp; + ParY[c][num_par[c]] = ParY_tmp; + ParZ[c][num_par[c]] = ParZ_tmp; + ParM[c][num_par[c]] = ParM_tmp; + VelX[c][num_par[c]] = VelX_tmp; + VelY[c][num_par[c]] = VelY_tmp; + VelZ[c][num_par[c]] = VelZ_tmp; + num_par[c] += 1; + } + if (num_par[c] >= N) { + N = num_par[c] + 1; // Increase the new maximum size if needed + ParX[c] = (double*)realloc(ParX[c], N * sizeof(double)); + ParY[c] = (double*)realloc(ParY[c], N * sizeof(double)); + ParZ[c] = (double*)realloc(ParZ[c], N * sizeof(double)); + ParM[c] = (double*)realloc(ParM[c], N * sizeof(double)); + VelX[c] = (double*)realloc(VelX[c], N * sizeof(double)); + VelY[c] = (double*)realloc(VelY[c], N * sizeof(double)); + VelZ[c] = (double*)realloc(VelZ[c], N * sizeof(double)); + } } } } } + +// Collect the number of target particles from each rank + MPI_Allreduce( num_par, num_par_sum, 3, MPI_INT, MPI_SUM, MPI_COMM_WORLD ); + + int num_par_eachRank[3][MPI_NRank]; + int displs[3][MPI_NRank]; + for (int c=0; cdh[MAX_LEVEL]; + for (int c=0; c soften ) pote_local[i-start] += ParM_sum[c][j]/rel_pos; + else if ( rel_pos <= soften && i != j ) pote_local[i-start] += ParM_sum[c][j]/soften; + } + pote_local[i-start] *= -NEWTON_G; + } + int recvcounts[MPI_NRank], displs[MPI_NRank]; + for (int i=0; i 1 && sqrt(dis[0]) < dis_exp && sqrt(dis[1]) < dis_exp ) IfConverge = true; + + for (int c=0; cPar->NPar_AcPlusInac; p++) { - if ( amr->Par->Type[p] == real(PTYPE_CEN+c) ) { - for (int d=0; d<3; d++) Cen_Tmp[d] = ParPos[d][p]; + if ( amr->Par->Mass[p] >= (real)0.0 && amr->Par->Type[p] == real(PTYPE_CEN+c) ){ + if ( CurrentMaxLv && AdjustPos == true ){ + amr->Par->PosX[p] = min_pos[c][0]; + amr->Par->PosY[p] = min_pos[c][1]; + amr->Par->PosZ[p] = min_pos[c][2]; + } + if ( CurrentMaxLv && AdjustVel == true ){ + amr->Par->VelX[p] = DM_Vel[c][0]; + amr->Par->VelY[p] = DM_Vel[c][1]; + amr->Par->VelZ[p] = DM_Vel[c][2]; + } + Cen_Tmp[0] = amr->Par->PosX[p]; + Cen_Tmp[1] = amr->Par->PosY[p]; + Cen_Tmp[2] = amr->Par->PosZ[p]; + Vel_Tmp[0] = amr->Par->VelX[p]; + Vel_Tmp[1] = amr->Par->VelY[p]; + Vel_Tmp[2] = amr->Par->VelZ[p]; break; } } - // use MPI_MAX since Cen_Tmp[] is initialized as -inf - MPI_Reduce( Cen_Tmp, Cen[c], 3, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD ); + +// use MPI_MAX since Cen_Tmp[] is initialized as -inf + MPI_Allreduce( Cen_Tmp, Cen_new[c], 3, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD ); + MPI_Allreduce( Vel_Tmp, Cen_Vel[c], 3, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD ); + } // for (int c=0; c + + +static void Output_ExactCooling(); +double Lambda(double Temp, double ZIRON); +double integrand(double Temp, void *params); + +// problem-specific global variables +// ======================================================================================= +static double EC_Temp; +static double EC_Dens; + int count = 0; +// ======================================================================================= + + + + +//------------------------------------------------------------------------------------------------------- +// Function : Validate +// Description : Validate the compilation flags and runtime parameters for this test problem +// +// Note : None +// +// Parameter : None +// +// Return : None +//------------------------------------------------------------------------------------------------------- +void Validate() +{ + + if ( MPI_Rank == 0 ) Aux_Message( stdout, " Validating test problem %d ...\n", TESTPROB_ID ); + +# if ( MODEL != HYDRO ) + Aux_Error( ERROR_INFO, "MODEL != HYDRO !!\n" ); +# endif + + if ( MPI_Rank == 0 ) Aux_Message( stdout, " Validating test problem %d ... done\n", TESTPROB_ID ); + +} // FUNCTION : Validate + + + +#if ( MODEL == HYDRO ) +//------------------------------------------------------------------------------------------------------- +// Function : SetParameter +// Description : Load and set the problem-specific runtime parameters +// +// Note : 1. Filename is set to "Input__TestProb" by default +// 2. Major tasks in this function: +// (1) load the problem-specific runtime parameters +// (2) set the problem-specific derived parameters +// (3) reset other general-purpose parameters if necessary +// (4) make a note of the problem-specific parameters +// 3. Must call EoS_Init() before calling any other EoS routine +// +// Parameter : None +// +// Return : None +//------------------------------------------------------------------------------------------------------- +void SetParameter() +{ + + if ( MPI_Rank == 0 ) Aux_Message( stdout, " Setting runtime parameters ...\n" ); + + +// (1) load the problem-specific runtime parameters + const char FileName[] = "Input__TestProb"; + ReadPara_t *ReadPara = new ReadPara_t; + +// (1-1) add parameters in the following format: +// --> note that VARIABLE, DEFAULT, MIN, and MAX must have the same data type +// --> some handy constants (e.g., Useless_bool, Eps_double, NoMin_int, ...) are defined in "include/ReadPara.h" +// ******************************************************************************************************************************** +// ReadPara->Add( "KEY_IN_THE_FILE", &VARIABLE, DEFAULT, MIN, MAX ); +// ******************************************************************************************************************************** + ReadPara->Add( "EC_Temp", &EC_Temp, 1000000.0, Eps_double, NoMax_double ); + ReadPara->Add( "EC_Dens", &EC_Dens, 1.0, Eps_double, NoMax_double ); + + ReadPara->Read( FileName ); + + delete ReadPara; + +// (1-2) set the default values + +// (1-3) check the runtime parameters + + +// (2) set the problem-specific derived parameters + + +// (3) reset other general-purpose parameters +// --> a helper macro PRINT_WARNING is defined in TestProb.h + const long End_Step_Default = __INT_MAX__; + const double End_T_Default = 100.0*Const_Myr/UNIT_T; + + if ( END_STEP < 0 ) { + END_STEP = End_Step_Default; + PRINT_WARNING( "END_STEP", END_STEP, FORMAT_LONG ); + } + + if ( END_T < 0.0 ) { + END_T = End_T_Default; + PRINT_WARNING( "END_T", END_T, FORMAT_REAL ); + } + + +// (4) make a note + if ( MPI_Rank == 0 ) + { + Aux_Message( stdout, "=============================================================================\n" ); + Aux_Message( stdout, " test problem ID = %d\n", TESTPROB_ID ); + Aux_Message( stdout, " EC_Temp = %13.7e\n", EC_Temp ); + Aux_Message( stdout, " EC_Dens = %13.7e\n", EC_Dens ); + Aux_Message( stdout, "=============================================================================\n" ); + } + + + if ( MPI_Rank == 0 ) Aux_Message( stdout, " Setting runtime parameters ... done\n" ); + +} // FUNCTION : SetParameter + + + +//------------------------------------------------------------------------------------------------------- +// Function : SetGridIC +// Description : Set the problem-specific initial condition on grids +// +// Note : 1. This function may also be used to estimate the numerical errors when OPT__OUTPUT_USER is enabled +// --> In this case, it should provide the analytical solution at the given "Time" +// 2. This function will be invoked by multiple OpenMP threads when OPENMP is enabled +// (unless OPT__INIT_GRID_WITH_OMP is disabled) +// --> Please ensure that everything here is thread-safe +// 3. Even when DUAL_ENERGY is adopted for HYDRO, one does NOT need to set the dual-energy variable here +// --> It will be calculated automatically +// 4. For MHD, do NOT add magnetic energy (i.e., 0.5*B^2) to fluid[ENGY] here +// --> It will be added automatically later +// +// Parameter : fluid : Fluid field to be initialized +// x/y/z : Physical coordinates +// Time : Physical time +// lv : Target refinement level +// AuxArray : Auxiliary array +// +// Return : fluid +//------------------------------------------------------------------------------------------------------- +void SetGridIC( real fluid[], const double x, const double y, const double z, const double Time, + const int lv, double AuxArray[] ) +{ + const double cl_X = 0.7; // mass-fraction of hydrogen + const double cl_Z = 0.018; // metallicity (in Zsun) + const double cl_mol = 1.0/(2*cl_X+0.75*(1-cl_X-cl_Z)+cl_Z*0.5); // mean (total) molecular weights + + double Dens, MomX, MomY, MomZ, Pres, Eint, Etot; +// Convert the input number density into mass density rho + double cl_dens = (EC_Dens*MU_NORM*cl_mol) / UNIT_D; +// double cl_dens = (EC_Dens*MU_NORM*0.61709348966) / UNIT_D; + double cl_pres = EoS_DensTemp2Pres_CPUPtr( cl_dens, EC_Temp, NULL, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); + + Dens = cl_dens; + MomX = 0.0; + MomY = 0.0; + MomZ = 0.0; + Pres = cl_pres; + Eint = EoS_DensPres2Eint_CPUPtr( Dens, Pres, NULL, EoS_AuxArray_Flt, + EoS_AuxArray_Int, h_EoS_Table ); // assuming EoS requires no passive scalars + Etot = Hydro_ConEint2Etot( Dens, MomX, MomY, MomZ, Eint, 0.0 ); // do NOT include magnetic energy here + +// set the output array + fluid[DENS] = Dens; + fluid[MOMX] = MomX; + fluid[MOMY] = MomY; + fluid[MOMZ] = MomZ; + fluid[ENGY] = Etot; + +} // FUNCTION : SetGridIC + + + +//------------------------------------------------------------------------------------------------------- +// Function : OutputExactCooling +// Description : Output the temperature relative error in the exact cooling problem +// +// Note : 1. Enabled by the runtime option "OPT__OUTPUT_USER" +// 2. Construct the analytical solution corresponding to the cooling function +// +// Parameter : None +// +// Return : None +//------------------------------------------------------------------------------------------------------- +void Output_ExactCooling() +{ + const char FileName[] = "Output__Error"; + static bool FirstTime = true; + +// header + if ( FirstTime ) { + if ( MPI_Rank == 0 ) { + if ( Aux_CheckFileExist(FileName) ) + Aux_Message( stderr, "WARNING : file \"%s\" already exists !!\n", FileName ); + + FILE *File_User = fopen( FileName, "a" ); + fprintf( File_User, "#%13s%10s ", "Time", "DumpID" ); + fprintf( File_User, "%14s %14s %14s %14s %14s", "Temp_nume", "Temp_anal", "Err", "Tcool_nume", "Tcool_anal"); + fprintf( File_User, "\n" ); + fclose( File_User ); + } + FirstTime = false; + } // if ( FirstTime ) + + + const double cl_X = 0.7; // mass-fraction of hydrogen + const double cl_Z = 0.018; // metallicity (in Zsun) + const double cl_mol = 1.0/(2*cl_X+0.75*(1-cl_X-cl_Z)+cl_Z*0.5); // mean (total) molecular weights + const double cl_mole = 2.0/(1+cl_X); // mean electron molecular weights + const double cl_moli = 1.0/cl_X; // mean proton molecular weights + const double cl_moli_mole = cl_moli*cl_mole; // Assume the molecular weights are constant, mu_e*mu_i = 1.464 + const double Temp_cut = 1e5; + +// Get the numerical result + real fluid[NCOMP_TOTAL]; + double Temp_nume = 0.0; + double Temp_nume_tmp = 0.0; + double Tcool_nume = 0.0; + int count = 0; + const int lv = 0; + + for (int k=1; kpatch[ amr->FluSg[lv] ][lv][0]->fluid[v][k][j][i]; + Temp_nume_tmp = (real) Hydro_Con2Temp( fluid[0], fluid[1], fluid[2], fluid[3], fluid[4], fluid+NCOMP_FLUID, + true, MIN_TEMP, 0.0, EoS_DensEint2Temp_CPUPtr, + EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); + Tcool_nume += 1.0/(GAMMA-1)*(Const_kB*cl_moli_mole*Temp_nume_tmp)/(fluid[0]*UNIT_D/MU_NORM*cl_mol*3.2217e-27*sqrt(Temp_nume_tmp))/Const_Myr; + Temp_nume += Temp_nume_tmp; + count += 1; + }}} + Temp_nume /= count; + Tcool_nume /= count; + +// (1) Compute the analytical solution for Sutherland-Dopita cooling function + double gsl_result, gsl_error; + double K = -(GAMMA-1)*EC_Dens*cl_mol*cl_mol/cl_moli_mole/Const_kB; + gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000); + gsl_function F; + F.function = &integrand; + gsl_integration_qags(&F, EC_Temp, Temp_nume, 0, 1e-10, 1000, w, &gsl_result, &gsl_error); + double Time_gsl = gsl_result/K; + +/* +// (2) Compute the analytical solution for single branch cooling function + double Temp_anal, Tcool_anal; + if (sqrt(EC_Temp) >= 3.2217e-27/2.0*(GAMMA-1)*EC_Dens*cl_mol*cl_mol/cl_mole/cl_moli/Const_kB*Time[0]*UNIT_T){ + Temp_anal = pow(sqrt(EC_Temp) - 3.2217e-27/2.0*(GAMMA-1)*EC_Dens*cl_mol*cl_mol/cl_mole/cl_moli/Const_kB*Time[0]*UNIT_T, 2.0); + if (Temp_anal < MIN_TEMP) Temp_anal = MIN_TEMP; + } + else Temp_anal = MIN_TEMP; + + Tcool_anal = 1.0/(GAMMA-1)*EC_Dens*Const_kB*Temp_anal/((EC_Dens*cl_mol/cl_mole)*(EC_Dens*cl_mol/cl_moli)*3.2217e-27*sqrt(Temp_anal))/Const_Myr; +*/ +/* +// (3) Compute the analytical solution for 2 branch cooling function + const int size_anal = 1001; + double time_anal[size_anal]; + double Temp_anal_arr[size_anal] = {0.0}; + double Tcool_anal_arr[size_anal] = {0.0}; + double time_cut, Temp_anal, Tcool_anal; + + for (int i=0; i= 3.2217e-27/2.0*(GAMMA-1)*EC_Dens*cl_mol*cl_mol/cl_mole/cl_moli/Const_kB*time_anal[i]*Const_Myr){ + Temp_anal_arr[i] = pow(sqrt(EC_Temp) - 3.2217e-27/2.0*(GAMMA-1)*EC_Dens*cl_mol*cl_mol/cl_mole/cl_moli/Const_kB*time_anal[i]*Const_Myr, 2.0); + if (Temp_anal_arr[i] < MIN_TEMP) Temp_anal_arr[i] = MIN_TEMP; + } + else Temp_anal_arr[i] = MIN_TEMP; + } + + for (int i=0; i= Temp_cut && Temp_anal_arr[i+1] <= Temp_cut ){ + time_cut = time_anal[i] + (time_anal[i+1]-time_anal[i])*(Temp_cut-Temp_anal_arr[i])/(Temp_anal_arr[i+1]-Temp_anal_arr[i]); + } + } + +// for (int i=0; i= 3.2217e-27/2.0*(GAMMA-1)*EC_Dens*cl_mol*cl_mol/cl_mole/cl_moli/Const_kB*Time[0]*UNIT_T){ + Temp_anal = pow(sqrt(EC_Temp) - 3.2217e-27/2.0*(GAMMA-1)*EC_Dens*cl_mol*cl_mol/cl_mole/cl_moli/Const_kB*Time[0]*UNIT_T, 2.0); + if (Temp_anal < MIN_TEMP) Temp_anal = MIN_TEMP; + } + else Temp_anal = MIN_TEMP; + + Tcool_anal = 1.0/(GAMMA-1)*EC_Dens*Const_kB*Temp_anal/((EC_Dens*cl_mol/cl_mole)*(EC_Dens*cl_mol/cl_moli)*3.2217e-27*sqrt(Temp_anal))/Const_Myr; + } + else { + if (pow(Temp_cut, 0.6) >= 3.2217e-27*0.6*(GAMMA-1)*EC_Dens*cl_mol*cl_mol/cl_mole/cl_moli/Const_kB*(Time[0]*UNIT_T-time_cut*Const_Myr)){ + Temp_anal = pow(pow(Temp_cut, 0.6) - 3.2217e-27*0.6*(GAMMA-1)*EC_Dens*cl_mol*cl_mol/cl_mole/cl_moli/Const_kB*(Time[0]*UNIT_T-time_cut*Const_Myr), 1.0/0.6); + if (Temp_anal < MIN_TEMP) Temp_anal = MIN_TEMP; + } + else Temp_anal = MIN_TEMP; + + Tcool_anal = 1.0/(GAMMA-1)*EC_Dens*Const_kB*Temp_anal/((EC_Dens*cl_mol/cl_mole)*(EC_Dens*cl_mol/cl_moli)*3.2217e-27*pow(Temp_anal, 0.4))/Const_Myr; + } +*/ + +// Record + if ( MPI_Rank == 0 ) { + FILE *File_User = fopen( FileName, "a" ); + fprintf( File_User, "%14.7e%10d ", Time[0]*UNIT_T/Const_Myr, DumpID ); +// fprintf( File_User, "%14.7e %14.7e %14.7e %14.7e %14.7e", Temp_nume, Temp_anal, (Temp_nume-Temp_anal)/Temp_anal, Tcool_nume, Tcool_anal ); + fprintf( File_User, "%14.7e %14.7e %14.7e %14.7e %14.7e", Temp_nume, Time_gsl/Const_Myr, (Time[0]*UNIT_T-Time_gsl)/Time_gsl, Tcool_nume, 0.0 ); + fprintf( File_User, "\n" ); + fclose( File_User ); + } + + +} // FUNCTION : Output_ExactCooling +#endif + + +double Lambda(double TEMP, double ZIRON){ + double TLOGC = 5.65; + double QLOGC = -21.566; + double QLOGINFTY = -23.1; + double PPP = 0.8; + double TLOGM = 5.1; + double QLOGM = -20.85; + double SIG = 0.65; + double TLOG = log10(TEMP); + + double QLOG0, QLOG1, QLAMBDA0, QLAMBDA1, ARG, BUMP1RHS, BUMP2LHS, Lambdat; + if (TLOG >= 6.1) QLOG0 = -26.39 + 0.471*log10(TEMP + 3.1623e6); + else if (TLOG >= 4.9){ + ARG = pow(10.0, (-(TLOG-4.9)/0.5)) + 0.077302; + QLOG0 = -22.16 + log10(ARG); + } + else if (TLOG >= 4.25){ + BUMP1RHS = -21.98 - ((TLOG-4.25)/0.55); + BUMP2LHS = -22.16 - pow((TLOG-4.9)/0.284, 2); + QLOG0 = fmax(BUMP1RHS, BUMP2LHS); + } + else QLOG0 = -21.98 - pow((TLOG-4.25)/0.2, 2); + + if (QLOG0 < -30.0) QLOG0 = -30.0; + QLAMBDA0 = pow(10.0, QLOG0); + + if (TLOG >= 5.65) { + QLOG1 = QLOGC - PPP * (TLOG - TLOGC); + QLOG1 = fmax(QLOG1, QLOGINFTY); + } else { + QLOG1 = QLOGM - pow((TLOG - TLOGM) / SIG, 2.0); + } + + if (QLOG1 < -30.0) QLOG1 = -30.0; + QLAMBDA1 = pow(10.0, QLOG1); + + Lambdat = QLAMBDA0 + ZIRON * QLAMBDA1; +// Lambdat = 3.2217e-27 * sqrt(TEMP); + + return Lambdat; +} + +double integrand(double Temp, void *params){ + double Lambda_T = Lambda(Temp, 0.018); + return 1.0/Lambda_T; +} + + + +//------------------------------------------------------------------------------------------------------- +// Function : Init_TestProb_Hydro_ExactCooling +// Description : Test problem initializer +// +// Note : None +// +// Parameter : None +// +// Return : None +//------------------------------------------------------------------------------------------------------- +void Init_TestProb_Hydro_ExactCooling() +{ + + if ( MPI_Rank == 0 ) Aux_Message( stdout, "%s ...\n", __FUNCTION__ ); + + +// validate the compilation flags and runtime parameters + Validate(); + + +// replace HYDRO by the target model (e.g., MHD/ELBDM) and also check other compilation flags if necessary (e.g., GRAVITY/PARTICLE) +# if ( MODEL == HYDRO ) +// set the problem-specific runtime parameters + SetParameter(); + + +// procedure to enable a problem-specific function: +// 1. define a user-specified function (example functions are given below) +// 2. declare its function prototype on the top of this file +// 3. set the corresponding function pointer below to the new problem-specific function +// 4. enable the corresponding runtime option in "Input__Parameter" +// --> for instance, enable OPT__OUTPUT_USER for Output_User_Ptr + + Init_Function_User_Ptr = SetGridIC; + Output_User_Ptr = Output_ExactCooling; +// End_User_Ptr = End_ClusterMerger; +// Aux_Record_User_Ptr = Aux_Record_ClusterMerger; + +# endif +/* + Init_Function_User_Ptr = SetGridIC; +# ifdef MHD + Init_Function_BField_User_Ptr = SetBFieldIC; +# endif +// comment out Init_ByFile_User_Ptr to use the default +// Init_ByFile_User_Ptr = NULL; // option: OPT__INIT=3; example: Init/Init_ByFile.cpp -> Init_ByFile_Default() + Init_Field_User_Ptr = NULL; // set NCOMP_PASSIVE_USER; example: TestProblem/Hydro/Plummer/Init_TestProb_Hydro_Plummer.cpp --> AddNewField() + Flag_Region_Ptr = NULL; // option: OPT__FLAG_REGION; example: Refing/Flag_Region.cpp + Flag_User_Ptr = NULL; // option: OPT__FLAG_USER; example: Refine/Flag_User.cpp + Mis_GetTimeStep_User_Ptr = NULL; // option: OPT__DT_USER; example: Miscellaneous/Mis_GetTimeStep_User.cpp + Mis_UserWorkBeforeNextLevel_Ptr = NULL; // example: Miscellaneous/Mis_UserWorkBeforeNextLevel.cpp + Mis_UserWorkBeforeNextSubstep_Ptr = NULL; // example: Miscellaneous/Mis_UserWorkBeforeNextSubstep.cpp + BC_User_Ptr = NULL; // option: OPT__BC_FLU_*=4; example: TestProblem/ELBDM/ExtPot/Init_TestProb_ELBDM_ExtPot.cpp --> BC() +# ifdef MHD + BC_BField_User_Ptr = NULL; // option: OPT__BC_FLU_*=4; +# endif + Flu_ResetByUser_Func_Ptr = NULL; // option: OPT__RESET_FLUID; example: Fluid/Flu_ResetByUser.cpp + Init_DerivedField_User_Ptr = NULL; // option: OPT__OUTPUT_USER_FIELD; example: Fluid/Flu_DerivedField_User.cpp + Output_User_Ptr = NULL; // option: OPT__OUTPUT_USER; example: TestProblem/Hydro/AcousticWave/Init_TestProb_Hydro_AcousticWave.cpp --> OutputError() + Output_UserWorkBeforeOutput_Ptr = NULL; // option: none; example: Output/Output_UserWorkBeforeOutput.cpp + Aux_Record_User_Ptr = NULL; // option: OPT__RECORD_USER; example: Auxiliary/Aux_Record_User.cpp + Init_User_Ptr = NULL; // option: none; example: none + End_User_Ptr = NULL; // option: none; example: TestProblem/Hydro/ClusterMerger_vs_Flash/Init_TestProb_ClusterMerger_vs_Flash.cpp --> End_ClusterMerger() +# ifdef GRAVITY + Init_ExtAcc_Ptr = NULL; // option: OPT__EXT_ACC; example: SelfGravity/CPU_Gravity/CPU_ExtAcc_PointMass.cpp + End_ExtAcc_Ptr = NULL; + Init_ExtPot_Ptr = NULL; // option: OPT__EXT_POT; example: SelfGravity/CPU_Poisson/CPU_ExtPot_PointMass.cpp + End_ExtPot_Ptr = NULL; + Poi_AddExtraMassForGravity_Ptr = NULL; // option: OPT__GRAVITY_EXTRA_MASS; example: none + Poi_UserWorkBeforePoisson_Ptr = NULL; // option: none; example: SelfGravity/Poi_UserWorkBeforePoisson.cpp +# endif +# ifdef PARTICLE + Par_Init_ByFunction_Ptr = NULL; // option: PAR_INIT=1; example: Particle/Par_Init_ByFunction.cpp + Par_Init_Attribute_User_Ptr = NULL; // set PAR_NATT_USER; example: TestProblem/Hydro/AGORA_IsolatedGalaxy/Init_TestProb_Hydro_AGORA_IsolatedGalaxy.cpp --> AddNewParticleAttribute() +# endif +# if ( EOS == EOS_USER ) + EoS_Init_Ptr = NULL; // option: EOS in the Makefile; example: EoS/User_Template/CPU_EoS_User_Template.cpp + EoS_End_Ptr = NULL; +# endif +# endif // #if ( MODEL == HYDRO ) + Src_Init_User_Ptr = NULL; // option: SRC_USER; example: SourceTerms/User_Template/CPU_Src_User_Template.cpp +*/ + + if ( MPI_Rank == 0 ) Aux_Message( stdout, "%s ... done\n", __FUNCTION__ ); + +} // FUNCTION : Init_TestProb_Hydro_ExactCooling diff --git a/tool/analysis/gamer_compare_data/Makefile b/tool/analysis/gamer_compare_data/Makefile index 7c8f4b37c0..c3569379bc 100644 --- a/tool/analysis/gamer_compare_data/Makefile +++ b/tool/analysis/gamer_compare_data/Makefile @@ -26,7 +26,7 @@ SIMU_OPTION += -DSUPPORT_HDF5 # siimulation parameters ####################################################################################################### NLEVEL := 10 # level : 0 ~ NLEVEL-1 -MAX_PATCH := 2000000 # maximum number of patches in each level +MAX_PATCH := 1000000 # maximum number of patches in each level NLEVEL := $(strip $(NLEVEL)) MAX_PATCH := $(strip $(MAX_PATCH))