diff --git a/doc/wiki/Runtime-Parameters-related/Runtime-Parameters:-All.md b/doc/wiki/Runtime-Parameters-related/Runtime-Parameters:-All.md index 9129bf3665..833a1ff590 100644 --- a/doc/wiki/Runtime-Parameters-related/Runtime-Parameters:-All.md +++ b/doc/wiki/Runtime-Parameters-related/Runtime-Parameters:-All.md @@ -54,6 +54,7 @@ For variables with `Default/Min/Max` labeled as `Depend`, click the parameter na | [[ DT__CR_DIFFUSION \| Runtime-Parameters:-Timestep#DT__CR_DIFFUSION ]] | 3.0e-1 | 0.0 | None | dt criterion: cosmic-ray diffusion CFL factor [0.3] ##CR_DIFFUSION only## | | [[ DT__FLUID \| Runtime-Parameters:-Timestep#DT__FLUID ]] | -1.0 | None | None | dt criterion: fluid solver CFL factor (<0=auto) [-1.0] | | [[ DT__FLUID_INIT \| Runtime-Parameters:-Timestep#DT__FLUID_INIT ]] | -1.0 | None | None | dt criterion: DT__FLUID at the first step (<0=auto) [-1.0] | +| [[ DT__GRACKLE_COOLING \| Runtime-Parameters:-Timestep#DT__GRACKLE_COOLING ]] | -1.0 | None | None | dt criterion: factor for Grackle cooling time (<0=off) [-1.0] ##SUPPORT_GRACKLE only## | | [[ DT__GRAVITY \| Runtime-Parameters:-Timestep#DT__GRAVITY ]] | -1.0 | None | None | dt criterion: gravity solver safety factor (<0=auto) [-1.0] | | DT__HYBRID_CFL | -1.0 | None | None | dt criterion: hybrid solver CFL factor (<0=auto) (diffusion) [-1.0] ## ELBDM_HYBRID ONLY## | | DT__HYBRID_CFL_INIT | -1.0 | None | None | dt criterion: DT__HYBRID_CFL in the first step (<0=auto) [-1.0] ## ELBDM_HYBRID ONLY## | @@ -122,16 +123,19 @@ For variables with `Default/Min/Max` labeled as `Depend`, click the parameter na | [[ GFUNC_COEFF0 \| Runtime-Parameters:-Gravity#GFUNC_COEFF0 ]] | -1.0 | None | None | Green's function coefficient at the origin for the isolated BC (<0=auto) [-1.0] | | [[ GPU_NSTREAM \| Runtime-Parameters:-GPU#GPU_NSTREAM ]] | -1 | None | None | number of CUDA streams for the asynchronous memory copy in GPU (<=0=auto) [-1] | | [[ GRACKLE_ACTIVATE \| Runtime-Parameters:-Chemistry-and-Radiation#GRACKLE_ACTIVATE ]] | 1 | None | None | enable Grackle [1] | -| GRACKLE_CIE_COOLING | 0 | None | None | 0: off; 1:on | +| [[ GRACKLE_CIE_COOLING \| Runtime-Parameters:-Chemistry-and-Radiation#GRACKLE_CIE_COOLING ]] | 0 | None | None | map to "cie_cooling" [0] | | [[ GRACKLE_CLOUDY_TABLE \| Runtime-Parameters:-Chemistry-and-Radiation#GRACKLE_CLOUDY_TABLE ]] | None | None | None | "grackle_data_file" | | [[ GRACKLE_CMB_FLOOR \| Runtime-Parameters:-Chemistry-and-Radiation#GRACKLE_CMB_FLOOR ]] | 1 | None | None | map to "cmb_temperature_floor" [1] | | [[ GRACKLE_COOLING \| Runtime-Parameters:-Chemistry-and-Radiation#GRACKLE_COOLING ]] | 1 | None | None | map to "with_radiative_cooling" [1] | -| GRACKLE_H2_OPA_APPROX | 0 | 0 | 1 | H2 opacity from Ripamonti+04; 0:off, 1:Ripomonti+04 | +| [[ GRACKLE_H2_OPA_APPROX \| Runtime-Parameters:-Chemistry-and-Radiation#GRACKLE_H2_OPA_APPROX ]] | 0 | 0 | 1 | map to "h2_optical_depth_approximation" (0=off, 1=Ripomonti+04) [0] | +| [[ GRACKLE_HYDROGEN_MFRAC \| Runtime-Parameters:-Chemistry-and-Radiation#GRACKLE_HYDROGEN_MFRAC ]] | 0.76 | 0.0 | 1.0 | map to "HydrogenFractionByMass" (NOT set to Grackle when GRACKLE_PRIMORDIAL==0) [0.76] | | [[ GRACKLE_METAL \| Runtime-Parameters:-Chemistry-and-Radiation#GRACKLE_METAL ]] | 0 | None | None | map to "metal_cooling" (must increase NCOMP_PASSIVE_USER by 1) [0] | | [[ GRACKLE_PE_HEATING \| Runtime-Parameters:-Chemistry-and-Radiation#GRACKLE_PE_HEATING ]] | 0 | None | None | map to "photoelectric_heating" [0] | -| [[ GRACKLE_PE_HEATING_RATE \| Runtime-Parameters:-Chemistry-and-Radiation#GRACKLE_PE_HEATING_RATE ]] | 8.5e-26 | 0.0 | None | map to "photoelectric_heating_rate (in erg/cm^3/s)" [8.5e-26] | -| [[ GRACKLE_PRIMORDIAL \| Runtime-Parameters:-Chemistry-and-Radiation#GRACKLE_PRIMORDIAL ]] | 0 | 0 | 3 | map to "primordial_chemistry" (0=Cloudy, 1/2/3=6-/9-/12-species) [0] (must increase NCOMP_PASSIVE_USER by 6/9/12, respectively) | -| GRACKLE_THREE_BODY_RATE | 0 | 0 | 5 | used Glover+08 rate | +| [[ GRACKLE_PE_HEATING_RATE \| Runtime-Parameters:-Chemistry-and-Radiation#GRACKLE_PE_HEATING_RATE ]] | 8.5e-26 | 0.0 | None | map to "photoelectric_heating_rate" (in erg/cm^3/s/n_H) [8.5e-26] | +| [[ GRACKLE_PRIMORDIAL \| Runtime-Parameters:-Chemistry-and-Radiation#GRACKLE_PRIMORDIAL ]] | 0 | 0 | 3 | map to "primordial_chemistry" (0=Cloudy, 1/2/3=6-/9-/12-species) (must increase NCOMP_PASSIVE_USER by 6/9/12, respectively) [0] | +| [[ GRACKLE_THREE_BODY_RATE \| Runtime-Parameters:-Chemistry-and-Radiation#GRACKLE_THREE_BODY_RATE ]] | 0 | 0 | 5 | map to "three_body_rate" (0=Abel+02, 1=Palla+83, 2=Cohen+83, 3=Flower+07, 4=Glover+08, 5=Forrey+13) [0] | +| [[ GRACKLE_USE_S_HEATING_RATE \| Runtime-Parameters:-Chemistry-and-Radiation#GRACKLE_USE_S_HEATING_RATE ]] | 0 | None | None | map to "use_specific_heating_rate" [0] | +| [[ GRACKLE_USE_V_HEATING_RATE \| Runtime-Parameters:-Chemistry-and-Radiation#GRACKLE_USE_V_HEATING_RATE ]] | 0 | None | None | map to "use_volumetric_heating_rate" [0] | | [[ GRACKLE_UV \| Runtime-Parameters:-Chemistry-and-Radiation#GRACKLE_UV ]] | 0 | None | None | map to "UVbackground" [0] | | [[ GRACKLE_VERBOSE \| Runtime-Parameters:-Chemistry-and-Radiation#GRACKLE_VERBOSE ]] | 1 | None | None | map to "grackle_verbose" [1] | @@ -229,6 +233,7 @@ For variables with `Default/Min/Max` labeled as `Depend`, click the parameter na | [[ OPT__FIXUP_FLUX \| Runtime-Parameters:-Hydro#OPT__FIXUP_FLUX ]] | Depend | Depend | Depend | correct coarse grids by the fine-grid boundary fluxes [1] ##HYDRO and ELBDM ONLY## | | [[ OPT__FIXUP_RESTRICT \| Runtime-Parameters:-Hydro#OPT__FIXUP_RESTRICT ]] | 1 | None | None | correct coarse grids by averaging the fine-grid data [1] | | [[ OPT__FLAG_ANGULAR \| Runtime-Parameters:-Refinement#OPT__FLAG_ANGULAR ]] | 0 | None | None | flag: angular resolution (Input__Flag_AngularResolution) [0] | +| [[ OPT__FLAG_COOLING_LEN \| Runtime-Parameters:-Refinement#OPT__FLAG_COOLING_LEN ]] | 0 | None | None | flag: cooling length (Input__Flag_CoolingLen) [0] ##HYDRO+SUPPORT_GRACKLE ONLY## | | [[ OPT__FLAG_CRAY \| Runtime-Parameters:-Refinement#OPT__FLAG_CRAY ]] | 0 | None | None | flag: cosmic-ray energy (Input__Flag_CRay) [0] ##COSMIC_RAY ONLY## | | [[ OPT__FLAG_CURRENT \| Runtime-Parameters:-Refinement#OPT__FLAG_CURRENT ]] | 0 | None | None | flag: current density in MHD (Input__Flag_Current) [0] ##MHD ONLY## | | OPT__FLAG_ENGY_DENSITY | 0 | None | None | flag: energy density (Input_Flag_EngyDensity) [0] ##ELBDM ONLY## | @@ -289,6 +294,9 @@ For variables with `Default/Min/Max` labeled as `Depend`, click the parameter na | [[ OPT__OUTPUT_DIVVEL \| Runtime-Parameters:-Outputs#OPT__OUTPUT_DIVVEL ]] | 0 | None | None | output divergence(velocity) [0] ##HYDRO ONLY## | | [[ OPT__OUTPUT_ENTHALPY \| Runtime-Parameters:-Outputs#OPT__OUTPUT_ENTHALPY ]] | 1 | None | None | output reduced enthalpy [1] ##SRHD ONLY## | | [[ OPT__OUTPUT_ENTR \| Runtime-Parameters:-Outputs#OPT__OUTPUT_ENTR ]] | 0 | None | None | output gas entropy [0] ##HYDRO ONLY## | +| [[ OPT__OUTPUT_GRACKLE_MU \| Runtime-Parameters:-Outputs#OPT__OUTPUT_GRACKLE_MU ]] | 0 | None | None | output the mean molecular weight calculated by Grackle [0] ##SUPPORT_GRACKLE only## | +| [[ OPT__OUTPUT_GRACKLE_TCOOL \| Runtime-Parameters:-Outputs#OPT__OUTPUT_GRACKLE_TCOOL ]] | 0 | None | None | output the cooling time calculated by Grackle [0] ##SUPPORT_GRACKLE only## | +| [[ OPT__OUTPUT_GRACKLE_TEMP \| Runtime-Parameters:-Outputs#OPT__OUTPUT_GRACKLE_TEMP ]] | 0 | None | None | output the temperature calculated by Grackle [0] ##SUPPORT_GRACKLE only## | | [[ OPT__OUTPUT_LORENTZ \| Runtime-Parameters:-Outputs#OPT__OUTPUT_LORENTZ ]] | 0 | None | None | output Lorentz factor [0] ##SRHD ONLY## | | [[ OPT__OUTPUT_MACH \| Runtime-Parameters:-Outputs#OPT__OUTPUT_MACH ]] | 0 | None | None | output mach number [0] ##HYDRO ONLY## | | [[ OPT__OUTPUT_MODE \| Runtime-Parameters:-Outputs#OPT__OUTPUT_MODE ]] | -1 | 1 | 3 | (1=const step, 2=const dt, 3=dump table) -> edit "Input__DumpTable" for 3 | @@ -341,6 +349,7 @@ For variables with `Default/Min/Max` labeled as `Depend`, click the parameter na | [[ OPT__UM_IC_NLEVEL \| Runtime-Parameters:-Initial-Conditions#OPT__UM_IC_NLEVEL ]] | 1 | 1 | None | number of AMR levels UM_IC [1] --> edit "Input__UM_IC_RefineRegion" if >1 | | [[ OPT__UM_IC_NVAR \| Runtime-Parameters:-Initial-Conditions#OPT__UM_IC_NVAR ]] | Depend | Depend | Depend | number of variables in UM_IC: (1~NCOMP_TOTAL; <=0=auto) [HYDRO=5+passive/ELBDM=2] | | [[ OPT__UM_IC_REFINE \| Runtime-Parameters:-Initial-Conditions#OPT__UM_IC_REFINE ]] | 1 | None | None | refine UM_IC from level OPT__UM_IC_LEVEL to MAX_LEVEL [1] | +| [[ OPT__UNFREEZE_GRACKLE \| Runtime-Parameters:-Chemistry-and-Radiation#OPT__UNFREEZE_GRACKLE ]] | 0 | None | None | allow the evolution by Grackle solver when the fluid is frozen (for OPT__FREEZE_FLUID==1 only) [0] | | [[ OPT__UNIT \| Runtime-Parameters:-Units#OPT__UNIT ]] | 0 | None | None | specify code units -> must set exactly 3 basic units below [0] ##USELESS FOR COMOVING## | | [[ OPT__VERBOSE \| Runtime-Parameters:-Miscellaneous#OPT__VERBOSE ]] | 0 | None | None | output the simulation progress in detail [0] | | [[ OUTPUT_DIR \| Runtime-Parameters:-Outputs#OUTPUT_DIR ]] | "." | None | None | set the output directory [.] | @@ -404,7 +413,7 @@ For variables with `Default/Min/Max` labeled as `Depend`, click the parameter na # T | Name | Default | Min | Max | Short description | | :--- | :--- | :--- | :--- | :--- | -| [[ TESTPROB_ID \| Runtime-Parameters:-General#TESTPROB_ID ]] | 0 | 0 | None | 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
20: HYDRO MHD Cosmic Ray Soundwave
21: HYDRO MHD Cosmic Ray Shocktube
23: HYDRO MHD Cosmic Ray Diffusion
100: HYDRO CDM cosmological simulation (+GRAVITY & COMOVING & PARTICLE)
101: HYDRO Zeldovich pancake collapse (+GRAVITY & COMOVING & PARTICLE)
1000: ELBDM external potential (+GRAVITY)
1001: ELBDM Jeans instability in the comoving frame (+GRAVITY, +COMOVING)
1002: ELBDM Jeans instability in the physical frame (+GRAVITY)
1003: ELBDM soliton merger (+GRAVITY)
1004: ELBDM self-similar halo (+GRAVITY, +COMOVING)
1005: ELBDM rotating vortex pair
1006: ELBDM vortex pair in linear motion
1007: ELBDM halo extracted from a large-scale structure simulation (+GRAVITY)
1008: ELBDM 1D Gaussian wave packet
1009: ELBDM large-scale structure simulation (+GRAVITY, +COMOVING)
1010: ELBDM plane wave
1011: ELBDM small wave perturbations on homogeneous background | +| [[ TESTPROB_ID \| Runtime-Parameters:-General#TESTPROB_ID ]] | 0 | 0 | None | 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
20: HYDRO MHD Cosmic Ray Soundwave
21: HYDRO MHD Cosmic Ray Shocktube
23: HYDRO MHD Cosmic Ray Diffusion
24: HYDRO Grackle Test (+SUPPORT_GRACKLE)
100: HYDRO CDM cosmological simulation (+GRAVITY & COMOVING & PARTICLE)
101: HYDRO Zeldovich pancake collapse (+GRAVITY & COMOVING & PARTICLE)
1000: ELBDM external potential (+GRAVITY)
1001: ELBDM Jeans instability in the comoving frame (+GRAVITY, +COMOVING)
1002: ELBDM Jeans instability in the physical frame (+GRAVITY)
1003: ELBDM soliton merger (+GRAVITY)
1004: ELBDM self-similar halo (+GRAVITY, +COMOVING)
1005: ELBDM rotating vortex pair
1006: ELBDM vortex pair in linear motion
1007: ELBDM halo extracted from a large-scale structure simulation (+GRAVITY)
1008: ELBDM 1D Gaussian wave packet
1009: ELBDM large-scale structure simulation (+GRAVITY, +COMOVING)
1010: ELBDM plane wave
1011: ELBDM small wave perturbations on homogeneous background | # U | Name | Default | Min | Max | Short description | diff --git a/doc/wiki/Runtime-Parameters-related/Runtime-Parameters:-Chemistry-and-Radiation.md b/doc/wiki/Runtime-Parameters-related/Runtime-Parameters:-Chemistry-and-Radiation.md index e097fd3a77..8588f2e2b1 100644 --- a/doc/wiki/Runtime-Parameters-related/Runtime-Parameters:-Chemistry-and-Radiation.md +++ b/doc/wiki/Runtime-Parameters-related/Runtime-Parameters:-Chemistry-and-Radiation.md @@ -9,6 +9,13 @@ Parameters described on this page: [GRACKLE_PE_HEATING](#GRACKLE_PE_HEATING),   [GRACKLE_PE_HEATING_RATE](#GRACKLE_PE_HEATING_RATE),   [GRACKLE_CLOUDY_TABLE](#GRACKLE_CLOUDY_TABLE)   +[GRACKLE_THREE_BODY_RATE](#GRACKLE_THREE_BODY_RATE)   +[GRACKLE_CIE_COOLING](#GRACKLE_CIE_COOLING)   +[GRACKLE_H2_OPA_APPROX](#GRACKLE_H2_OPA_APPROX)   +[GRACKLE_USE_V_HEATING_RATE](#GRACKLE_USE_V_HEATING_RATE)   +[GRACKLE_USE_S_HEATING_RATE](#GRACKLE_USE_S_HEATING_RATE)   +[GRACKLE_HYDROGEN_MFRAC](#GRACKLE_HYDROGEN_MFRAC)   +[OPT__UNFREEZE_GRACKLE](#OPT__UNFREEZE_GRACKLE)   Parameters below are shown in the format:   **`Name`   (Valid Values)   [Default Value]** @@ -24,28 +31,28 @@ Only applicable when enabling the compilation option * #### `GRACKLE_VERBOSE`   (0=off, 1=on)   [1] * **Description:** -Map to the "grackle_verbose" option in GRACKLE. +Map to the ["grackle_verbose" option in GRACKLE](https://grackle.readthedocs.io/en/latest/Interaction.html#enabling-output). * **Restriction:** * #### `GRACKLE_COOLING`   (0=off, 1=on)   [1] * **Description:** -Map to the "with_radiative_cooling" runtime parameter in GRACKLE. +Map to the ["with_radiative_cooling" runtime parameter in GRACKLE](https://grackle.readthedocs.io/en/latest/Parameters.html#c.with_radiative_cooling). * **Restriction:** * #### `GRACKLE_PRIMORDIAL`   (0=Cloudy, 1=6-species, 2=9-species, 3=12-species)   [0] * **Description:** -Map to the "primordial_chemistry" runtime parameter in GRACKLE. +Map to the ["primordial_chemistry" runtime parameter in GRACKLE](https://grackle.readthedocs.io/en/latest/Parameters.html#c.primordial_chemistry). One must increase [[--passive | Installation:-Option-List#--passive]] -by 3, 6, or 9 for GRACKLE_PRIMORDIAL=1, 2, or 3, respectively. +by 6, 9, or 12 for GRACKLE_PRIMORDIAL=1, 2, or 3, respectively. * **Restriction:** * #### `GRACKLE_METAL`   (0=off, 1=on)   [0] * **Description:** -Map to the "metal_cooling" runtime parameter in GRACKLE. One must increase +Map to the ["metal_cooling" runtime parameter in GRACKLE](https://grackle.readthedocs.io/en/latest/Parameters.html#c.metal_cooling). One must increase [[--passive | Installation:-Option-List#--passive]] by 1 and initialize the field `Metal` using the field index `Idx_Metal` properly. * **Restriction:** @@ -53,33 +60,81 @@ by 1 and initialize the field `Metal` using the field index `Idx_Metal` properly * #### `GRACKLE_UV`   (0=off, 1=on)   [0] * **Description:** -Map to the "UVbackground" runtime parameter in GRACKLE. +Map to the ["UVbackground" runtime parameter in GRACKLE](https://grackle.readthedocs.io/en/latest/Parameters.html#c.UVbackground). * **Restriction:** * #### `GRACKLE_CMB_FLOOR`   (0=off, 1=on)   [1] * **Description:** -Map to the "cmb_temperature_floor" runtime parameter in GRACKLE. +Map to the ["cmb_temperature_floor" runtime parameter in GRACKLE](https://grackle.readthedocs.io/en/latest/Parameters.html#c.cmb_temperature_floor). * **Restriction:** * #### `GRACKLE_PE_HEATING`   (0=off, 1=on)   [0] * **Description:** -Map to the "photoelectric_heating" runtime parameter in GRACKLE. +Map to the ["photoelectric_heating" runtime parameter in GRACKLE](https://grackle.readthedocs.io/en/latest/Parameters.html#c.photoelectric_heating). * **Restriction:** * #### `GRACKLE_PE_HEATING_RATE`   (≥0.0)   [8.5e-26] * **Description:** -Map to the "photoelectric_heating_rate" runtime parameter in GRACKLE. +Map to the ["photoelectric_heating_rate" runtime parameter in GRACKLE](https://grackle.readthedocs.io/en/latest/Parameters.html#c.photoelectric_heating_rate). Note that the input value should always be in units of -erg cm-3 s-1. +erg cm-3 s-1 nH-1. * **Restriction:** * #### `GRACKLE_CLOUDY_TABLE`   (string)   [none] * **Description:** -Map to the "grackle_data_file" runtime parameter in GRACKLE. +Map to the ["grackle_data_file" runtime parameter in GRACKLE](https://grackle.readthedocs.io/en/latest/Parameters.html#c.grackle_data_file). + * **Restriction:** + + +* #### `GRACKLE_THREE_BODY_RATE`   (0=Abel+02, 1=Palla+83, 2=Cohen+83, 3=Flower+07, 4=Glover+08, 5=Forrey+13)   [0] + * **Description:** +Map to the ["three_body_rate" runtime parameter in GRACKLE](https://grackle.readthedocs.io/en/latest/Parameters.html#c.three_body_rate). + * **Restriction:** + + +* #### `GRACKLE_CIE_COOLING`   (0=off, 1=on)   [0] + * **Description:** +Map to the ["cie_cooling" runtime parameter in GRACKLE](https://grackle.readthedocs.io/en/latest/Parameters.html#c.cie_cooling). + * **Restriction:** + + +* #### `GRACKLE_H2_OPA_APPROX`   (0=off, 1=Ripomonti+04)   [0] + * **Description:** +Map to the ["h2_optical_depth_approximation" runtime parameter in GRACKLE](https://grackle.readthedocs.io/en/latest/Parameters.html#c.h2_optical_depth_approximation). + * **Restriction:** + + +* #### `GRACKLE_USE_V_HEATING_RATE`   (0=off, 1=on)   [0] + * **Description:** +Map to the ["use_volumetric_heating_rate" runtime parameter in GRACKLE](https://grackle.readthedocs.io/en/latest/Parameters.html#c.use_volumetric_heating_rate). + * **Restriction:** + + +* #### `GRACKLE_USE_S_HEATING_RATE`   (0=off, 1=on)   [0] + * **Description:** +Map to the ["use_specific_heating_rate" runtime parameter in GRACKLE](https://grackle.readthedocs.io/en/latest/Parameters.html#c.use_specific_heating_rate). + * **Restriction:** + + +* #### `GRACKLE_PE_HEATING_RATE`   (≥0.0)   [0.76] + * **Description:** +Map to the ["HydrogenFractionByMass" runtime parameter in GRACKLE](https://grackle.readthedocs.io/en/latest/Parameters.html#c.HydrogenFractionByMass). +Note that the input value will only set to Grackle when it is in non-equilibrium mode [GRACKLE_PRIMORDIAL](#GRACKLE_PRIMORDIAL)>0, +because the tables for tabulated mode were created assuming hydrogen mass fraction of about 0.716. +See the Grackle document for the details. +When in tabulated mode ([GRACKLE_PRIMORDIAL](#GRACKLE_PRIMORDIAL)==0), this value may still be used elsewhere in GAMER and +it may be different from HydrogenFractionByMass inside Grackle. + * **Restriction:** + + +* #### `OPT__UNFREEZE_GRACKLE`   (0=off, 1=on)   [0] + * **Description:** +Allow the evolution by Grackle solver when the fluid is frozen +(for [[OPT__FREEZE_FLUID | Runtime-Parameters:-Hydro#OPT__FREEZE_FLUID]]==1 only). * **Restriction:** diff --git a/doc/wiki/Runtime-Parameters-related/Runtime-Parameters:-Outputs.md b/doc/wiki/Runtime-Parameters-related/Runtime-Parameters:-Outputs.md index 71dbedb1c9..7eb26c61b2 100644 --- a/doc/wiki/Runtime-Parameters-related/Runtime-Parameters:-Outputs.md +++ b/doc/wiki/Runtime-Parameters-related/Runtime-Parameters:-Outputs.md @@ -21,6 +21,9 @@ Parameters described on this page: [OPT__OUTPUT_LORENTZ](#OPT__OUTPUT_LORENTZ),   [OPT__OUTPUT_3VELOCITY](#OPT__OUTPUT_3VELOCITY),   [OPT__OUTPUT_ENTHALPY](#OPT__OUTPUT_ENTHALPY),   +[OPT__OUTPUT_GRACKLE_TEMP](#OPT__OUTPUT_GRACKLE_TEMP),   +[OPT__OUTPUT_GRACKLE_MU](#OPT__OUTPUT_GRACKLE_MU),   +[OPT__OUTPUT_GRACKLE_TCOOL](#OPT__OUTPUT_GRACKLE_TCOOL),   [OPT__OUTPUT_USER_FIELD](#OPT__OUTPUT_USER_FIELD),   [OPT__OUTPUT_MODE](#OPT__OUTPUT_MODE),   [OPT__OUTPUT_RESTART](#OPT__OUTPUT_RESTART),   @@ -211,6 +214,27 @@ Output the SRHD reduced enthalpy on grids. * **Restriction:** For [[--srhd | Installation:-Option-List#--srhd]] only. + +* #### `OPT__OUTPUT_GRACKLE_TEMP`   (0=off, 1=on)   [0] + * **Description:** +Output the temperature calculated by Grackle. + * **Restriction:** +For [[--grackle | Installation:-Option-List#--grackle]] only. + + +* #### `OPT__OUTPUT_GRACKLE_MU`   (0=off, 1=on)   [0] + * **Description:** +Output the mean molecular weight calculated by Grackle. + * **Restriction:** +For [[--grackle | Installation:-Option-List#--grackle]] only. + + +* #### `OPT__OUTPUT_GRACKLE_TCOOL`   (0=off, 1=on)   [0] + * **Description:** +Output the cooling time calculated by Grackle. + * **Restriction:** +For [[--grackle | Installation:-Option-List#--grackle]] only. + * #### `OPT__OUTPUT_MODE`   (1=const step, 2=const dt, 3=dump table)   [none] * **Description:** diff --git a/doc/wiki/Runtime-Parameters-related/Runtime-Parameters:-Refinement.md b/doc/wiki/Runtime-Parameters-related/Runtime-Parameters:-Refinement.md index 9c8ff3f3b8..c9161b4dfe 100644 --- a/doc/wiki/Runtime-Parameters-related/Runtime-Parameters:-Refinement.md +++ b/doc/wiki/Runtime-Parameters-related/Runtime-Parameters:-Refinement.md @@ -34,6 +34,7 @@ Parameters described on this page: [OPT__FLAG_LRTZ_GRADIENT](#OPT__FLAG_LRTZ_GRADIENT),   [OPT__FLAG_VORTICITY](#OPT__FLAG_VORTICITY),   [OPT__FLAG_JEANS](#OPT__FLAG_JEANS),   +[OPT__FLAG_COOLING_LEN](#OPT__FLAG_COOLING_LEN),   [OPT__FLAG_CURRENT](#OPT__FLAG_CURRENT),   [OPT__FLAG_CRAY](#OPT__FLAG_CRAY),   [OPT__FLAG_LOHNER_DENS](#OPT__FLAG_LOHNER_DENS),   @@ -225,6 +226,23 @@ An example file can be found at `example/input/Input__Flag_Jeans`. Recommended values: ≥4. * **Restriction:** + +* #### `OPT__FLAG_COOLING_LEN`   (0=off, 1=on)   [0] + * **Description:** +Refinement criterion: gas cooling length. It ensures that the cooling length +is resolved by at least N cells. Specifically, a cell +on level l will be flagged for refinement if its estimated +cooling length Lcool satisfies +Lcool ≡tcoolcs  < NlΔξl, +where tcool is the cooling time (currently calculated by Grackle), cs is sound speed of gas, Δξl is the cell width along ξ on level l, and Nl is the refinement threshold on level l. +Specify the refinement +thresholds on different levels in the input file `Input__Flag_CoolingLen` +with the [[specific format | Runtime-Parameters:-Input__Flag_*]]. +An example file can be found at `example/input/Input__Flag_CoolingLen`. +Recommended values: ≥1. + * **Restriction:** +Must compile with [[--grackle | Installation:-Option-List#--grackle]]. + * #### `OPT__FLAG_CURRENT`   (0=off, 1=on)   [0] * **Description:** diff --git a/doc/wiki/Runtime-Parameters-related/Runtime-Parameters:-Timestep.md b/doc/wiki/Runtime-Parameters-related/Runtime-Parameters:-Timestep.md index ba8000e787..094a5cba4e 100644 --- a/doc/wiki/Runtime-Parameters-related/Runtime-Parameters:-Timestep.md +++ b/doc/wiki/Runtime-Parameters-related/Runtime-Parameters:-Timestep.md @@ -18,6 +18,7 @@ Parameters described on this page: [DT__PARVEL](#DT__PARVEL),   [DT__PARVEL_MAX](#DT__PARVEL_MAX),   [DT__PARACC](#DT__PARACC),   +[DT__GRACKLE_COOLING](#DT__GRACKLE_COOLING),   [DT__CR_DIFFUSION](#DT__CR_DIFFUSION),   [DT__MAX_DELTA_A](#DT__MAX_DELTA_A),   [DT__SYNC_PARENT_LV](#DT__SYNC_PARENT_LV),   @@ -108,6 +109,14 @@ for the exact formula. Only applicable when adopting the compilation option [[--store_par_acc | Installation:-Option-List#--store_par_acc]]. + +* #### `DT__GRACKLE_COOLING`   (≥0.0; <0.0 → off)   [-1.0] + * **Description:** +Factor for Grackle cooling time. + * **Restriction:** +Only applicable when adopting the compilation option +[[--grackle | Installation:-Option-List#--grackle]]. + * #### `DT__CR_DIFFUSION`   (≥0.0)   [0.3] * **Description:** diff --git a/example/input/Input__Flag_CoolingLen b/example/input/Input__Flag_CoolingLen new file mode 100644 index 0000000000..9e1eeeca26 --- /dev/null +++ b/example/input/Input__Flag_CoolingLen @@ -0,0 +1,13 @@ +# Level cooling_length/cell_size + 0 1.0 + 1 1.0 + 2 1.0 + 3 1.0 + 4 1.0 + 5 1.0 + 6 1.0 + 7 1.0 + 8 1.0 + 9 1.0 + 10 1.0 + 11 1.0 diff --git a/example/test_problem/Hydro/AGORA_IsolatedGalaxy/Input_Options/Input__Parameter.high-res b/example/test_problem/Hydro/AGORA_IsolatedGalaxy/Input_Options/Input__Parameter.high-res index 40ea1be1a9..deabb8c93a 100644 --- a/example/test_problem/Hydro/AGORA_IsolatedGalaxy/Input_Options/Input__Parameter.high-res +++ b/example/test_problem/Hydro/AGORA_IsolatedGalaxy/Input_Options/Input__Parameter.high-res @@ -131,8 +131,15 @@ GRACKLE_METAL 1 # ... "metal_cooling" (must increas GRACKLE_UV 1 # ... "UVbackground" [0] GRACKLE_CMB_FLOOR 1 # ... "cmb_temperature_floor" [1] GRACKLE_PE_HEATING 1 # ... "photoelectric_heating" [0] -GRACKLE_PE_HEATING_RATE 8.5e-26 # ... "photoelectric_heating_rate (in erg/cm^3/s)" [8.5e-26] +GRACKLE_PE_HEATING_RATE 8.5e-26 # ... "photoelectric_heating_rate" (in erg/cm^3/s/n_H) [8.5e-26] GRACKLE_CLOUDY_TABLE CloudyData_UVB=HM2012.h5 # "grackle_data_file" +GRACKLE_THREE_BODY_RATE 0 # map to "three_body_rate" (0=Abel+02, 1=Palla+83, 2=Cohen+83, 3=Flower+07, 4=Glover+08, 5=Forrey+13) [0] +GRACKLE_CIE_COOLING 0 # map to "cie_cooling" [0] +GRACKLE_H2_OPA_APPROX 0 # map to "h2_optical_depth_approximation" (0=off, 1=Ripomonti+04) [0] +GRACKLE_USE_V_HEATING_RATE 0 # map to "use_volumetric_heating_rate" [0] +GRACKLE_USE_S_HEATING_RATE 0 # map to "use_specific_heating_rate" [0] +GRACKLE_HYDROGEN_MFRAC 0.76 # map to "HydrogenFractionByMass" (NOT set to Grackle when GRACKLE_PRIMORDIAL==0) [0.76] +OPT__UNFREEZE_GRACKLE 0 # allow the evolution by Grackle solver when the fluid is frozen (for OPT__FREEZE_FLUID==1 only) [0] CHE_GPU_NPGROUP -1 # number of patch groups sent into the CPU/GPU Grackle solver (<=0=auto) [-1] diff --git a/example/test_problem/Hydro/AGORA_IsolatedGalaxy/Input_Options/Input__Parameter.low-res b/example/test_problem/Hydro/AGORA_IsolatedGalaxy/Input_Options/Input__Parameter.low-res index 41980091c2..2e256ec991 100644 --- a/example/test_problem/Hydro/AGORA_IsolatedGalaxy/Input_Options/Input__Parameter.low-res +++ b/example/test_problem/Hydro/AGORA_IsolatedGalaxy/Input_Options/Input__Parameter.low-res @@ -131,8 +131,15 @@ GRACKLE_METAL 1 # ... "metal_cooling" (must increas GRACKLE_UV 1 # ... "UVbackground" [0] GRACKLE_CMB_FLOOR 1 # ... "cmb_temperature_floor" [1] GRACKLE_PE_HEATING 1 # ... "photoelectric_heating" [0] -GRACKLE_PE_HEATING_RATE 8.5e-26 # ... "photoelectric_heating_rate (in erg/cm^3/s)" [8.5e-26] +GRACKLE_PE_HEATING_RATE 8.5e-26 # ... "photoelectric_heating_rate" (in erg/cm^3/s/n_H) [8.5e-26] GRACKLE_CLOUDY_TABLE CloudyData_UVB=HM2012.h5 # "grackle_data_file" +GRACKLE_THREE_BODY_RATE 0 # map to "three_body_rate" (0=Abel+02, 1=Palla+83, 2=Cohen+83, 3=Flower+07, 4=Glover+08, 5=Forrey+13) [0] +GRACKLE_CIE_COOLING 0 # map to "cie_cooling" [0] +GRACKLE_H2_OPA_APPROX 0 # map to "h2_optical_depth_approximation" (0=off, 1=Ripomonti+04) [0] +GRACKLE_USE_V_HEATING_RATE 0 # map to "use_volumetric_heating_rate" [0] +GRACKLE_USE_S_HEATING_RATE 0 # map to "use_specific_heating_rate" [0] +GRACKLE_HYDROGEN_MFRAC 0.76 # map to "HydrogenFractionByMass" (NOT set to Grackle when GRACKLE_PRIMORDIAL==0) [0.76] +OPT__UNFREEZE_GRACKLE 0 # allow the evolution by Grackle solver when the fluid is frozen (for OPT__FREEZE_FLUID==1 only) [0] CHE_GPU_NPGROUP -1 # number of patch groups sent into the CPU/GPU Grackle solver (<=0=auto) [-1] diff --git a/example/test_problem/Hydro/CMZ/Input__Parameter b/example/test_problem/Hydro/CMZ/Input__Parameter index d6ff47aa8e..11d7c374d9 100644 --- a/example/test_problem/Hydro/CMZ/Input__Parameter +++ b/example/test_problem/Hydro/CMZ/Input__Parameter @@ -137,11 +137,15 @@ GRACKLE_METAL 1 # ... "metal_cooling" (must increas GRACKLE_UV 1 # ... "UVbackground" [0] GRACKLE_CMB_FLOOR 1 # ... "cmb_temperature_floor" [1] GRACKLE_PE_HEATING 1 # ... "photoelectric_heating" [0] -GRACKLE_PE_HEATING_RATE 8.5e-26 # ... "photoelectric_heating_rate (in erg/cm^3/s)" [8.5e-26] +GRACKLE_PE_HEATING_RATE 8.5e-26 # ... "photoelectric_heating_rate" (in erg/cm^3/s/n_H) [8.5e-26] GRACKLE_CLOUDY_TABLE CloudyData_UVB=HM2012.h5 # "grackle_data_file" -GRACKLE_THREE_BODY_RATE 0 # used Glover+08 rate -GRACKLE_CIE_COOLING 0 # 0: off; 1:on -GRACKLE_H2_OPA_APPROX 0 # H2 opacity from Ripamonti+04; 0:off, 1:Ripomonti+04 +GRACKLE_THREE_BODY_RATE 0 # map to "three_body_rate" (0=Abel+02, 1=Palla+83, 2=Cohen+83, 3=Flower+07, 4=Glover+08, 5=Forrey+13) [0] +GRACKLE_CIE_COOLING 0 # map to "cie_cooling" [0] +GRACKLE_H2_OPA_APPROX 0 # map to "h2_optical_depth_approximation" (0=off, 1=Ripomonti+04) [0] +GRACKLE_USE_V_HEATING_RATE 0 # map to "use_volumetric_heating_rate" [0] +GRACKLE_USE_S_HEATING_RATE 0 # map to "use_specific_heating_rate" [0] +GRACKLE_HYDROGEN_MFRAC 0.76 # map to "HydrogenFractionByMass" (NOT set to Grackle when GRACKLE_PRIMORDIAL==0) [0.76] +OPT__UNFREEZE_GRACKLE 0 # allow the evolution by Grackle solver when the fluid is frozen (for OPT__FREEZE_FLUID==1 only) [0] CHE_GPU_NPGROUP -1 # number of patch groups sent into the CPU/GPU Grackle solver (<=0=auto) [-1] diff --git a/example/test_problem/Hydro/CR_Diffusion/Input__Parameter b/example/test_problem/Hydro/CR_Diffusion/Input__Parameter index 1d4e672346..a5f483f1fb 100644 --- a/example/test_problem/Hydro/CR_Diffusion/Input__Parameter +++ b/example/test_problem/Hydro/CR_Diffusion/Input__Parameter @@ -160,11 +160,15 @@ GRACKLE_METAL 0 # ... "metal_cooling" (must increas 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_PE_HEATING_RATE 8.5e-26 # ... "photoelectric_heating_rate" (in erg/cm^3/s/n_H) [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 +GRACKLE_THREE_BODY_RATE 0 # map to "three_body_rate" (0=Abel+02, 1=Palla+83, 2=Cohen+83, 3=Flower+07, 4=Glover+08, 5=Forrey+13) [0] +GRACKLE_CIE_COOLING 0 # map to "cie_cooling" [0] +GRACKLE_H2_OPA_APPROX 0 # map to "h2_optical_depth_approximation" (0=off, 1=Ripomonti+04) [0] +GRACKLE_USE_V_HEATING_RATE 0 # map to "use_volumetric_heating_rate" [0] +GRACKLE_USE_S_HEATING_RATE 0 # map to "use_specific_heating_rate" [0] +GRACKLE_HYDROGEN_MFRAC 0.76 # map to "HydrogenFractionByMass" (NOT set to Grackle when GRACKLE_PRIMORDIAL==0) [0.76] +OPT__UNFREEZE_GRACKLE 0 # allow the evolution by Grackle solver when the fluid is frozen (for OPT__FREEZE_FLUID==1 only) [0] CHE_GPU_NPGROUP -1 # number of patch groups sent into the CPU/GPU Grackle solver (<=0=auto) [-1] diff --git a/example/test_problem/Hydro/CR_ShockTube/Input__Parameter b/example/test_problem/Hydro/CR_ShockTube/Input__Parameter index 0360331616..3d6daf305d 100644 --- a/example/test_problem/Hydro/CR_ShockTube/Input__Parameter +++ b/example/test_problem/Hydro/CR_ShockTube/Input__Parameter @@ -159,11 +159,15 @@ GRACKLE_METAL 0 # ... "metal_cooling" (must increas 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_PE_HEATING_RATE 8.5e-26 # ... "photoelectric_heating_rate" (in erg/cm^3/s/n_H) [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 +GRACKLE_THREE_BODY_RATE 0 # map to "three_body_rate" (0=Abel+02, 1=Palla+83, 2=Cohen+83, 3=Flower+07, 4=Glover+08, 5=Forrey+13) [0] +GRACKLE_CIE_COOLING 0 # map to "cie_cooling" [0] +GRACKLE_H2_OPA_APPROX 0 # map to "h2_optical_depth_approximation" (0=off, 1=Ripomonti+04) [0] +GRACKLE_USE_V_HEATING_RATE 0 # map to "use_volumetric_heating_rate" [0] +GRACKLE_USE_S_HEATING_RATE 0 # map to "use_specific_heating_rate" [0] +GRACKLE_HYDROGEN_MFRAC 0.76 # map to "HydrogenFractionByMass" (NOT set to Grackle when GRACKLE_PRIMORDIAL==0) [0.76] +OPT__UNFREEZE_GRACKLE 0 # allow the evolution by Grackle solver when the fluid is frozen (for OPT__FREEZE_FLUID==1 only) [0] CHE_GPU_NPGROUP -1 # number of patch groups sent into the CPU/GPU Grackle solver (<=0=auto) [-1] diff --git a/example/test_problem/Hydro/CR_SoundWave/Input__Parameter b/example/test_problem/Hydro/CR_SoundWave/Input__Parameter index acc884ac52..6bd277bb8c 100644 --- a/example/test_problem/Hydro/CR_SoundWave/Input__Parameter +++ b/example/test_problem/Hydro/CR_SoundWave/Input__Parameter @@ -158,11 +158,15 @@ GRACKLE_METAL 0 # ... "metal_cooling" (must increas 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_PE_HEATING_RATE 8.5e-26 # ... "photoelectric_heating_rate" (in erg/cm^3/s/n_H) [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 +GRACKLE_THREE_BODY_RATE 0 # map to "three_body_rate" (0=Abel+02, 1=Palla+83, 2=Cohen+83, 3=Flower+07, 4=Glover+08, 5=Forrey+13) [0] +GRACKLE_CIE_COOLING 0 # map to "cie_cooling" [0] +GRACKLE_H2_OPA_APPROX 0 # map to "h2_optical_depth_approximation" (0=off, 1=Ripomonti+04) [0] +GRACKLE_USE_V_HEATING_RATE 0 # map to "use_volumetric_heating_rate" [0] +GRACKLE_USE_S_HEATING_RATE 0 # map to "use_specific_heating_rate" [0] +GRACKLE_HYDROGEN_MFRAC 0.76 # map to "HydrogenFractionByMass" (NOT set to Grackle when GRACKLE_PRIMORDIAL==0) [0.76] +OPT__UNFREEZE_GRACKLE 0 # allow the evolution by Grackle solver when the fluid is frozen (for OPT__FREEZE_FLUID==1 only) [0] CHE_GPU_NPGROUP -1 # number of patch groups sent into the CPU/GPU Grackle solver (<=0=auto) [-1] diff --git a/example/test_problem/Hydro/GrackleTest/Input__Flag_CoolingLen b/example/test_problem/Hydro/GrackleTest/Input__Flag_CoolingLen new file mode 100644 index 0000000000..3836537ce2 --- /dev/null +++ b/example/test_problem/Hydro/GrackleTest/Input__Flag_CoolingLen @@ -0,0 +1,13 @@ +# Level cooling_length/cell_size + 0 8.0 + 1 8.0 + 2 8.0 + 3 8.0 + 4 8.0 + 5 8.0 + 6 8.0 + 7 8.0 + 8 8.0 + 9 8.0 + 10 8.0 + 11 8.0 diff --git a/example/test_problem/Hydro/GrackleTest/Input__Parameter b/example/test_problem/Hydro/GrackleTest/Input__Parameter new file mode 100644 index 0000000000..3429a6051d --- /dev/null +++ b/example/test_problem/Hydro/GrackleTest/Input__Parameter @@ -0,0 +1,239 @@ + + +# ================================================================================================================= +# 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) +# 6. For boolean options: 0/1 -> off/on +# ================================================================================================================= + + +# simulation scale +BOX_SIZE 8.0 # box size along the longest side (in Mpc/h if COMOVING is adopted) +NX0_TOT_X 128 # number of base-level cells along x +NX0_TOT_Y 128 # number of base-level cells along y +NX0_TOT_Z 32 # number of base-level cells along z +OMP_NTHREAD -1 # number of OpenMP threads (<=0=auto) [-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 24 # test problem ID [0] + # 24: HYDRO Grackle Test (+SUPPORT_GRACKLE) + + +# code units (in cgs) +OPT__UNIT 1 # specify code units -> must set exactly 3 basic units below [0] ##USELESS FOR COMOVING## +UNIT_L 3.08567758149e21 # length unit (<=0 -> set to UNIT_V*UNIT_T or (UNIT_M/UNIT_D)^(1/3)) [-1.0] +UNIT_M 1.9885e42 # mass unit (<=0 -> set to UNIT_D*UNIT_L^3) [-1.0] +UNIT_T 3.15569252e13 # 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) +OPT__BC_FLU_XP 1 # fluid boundary condition at the +x face: (1=periodic, 2=outflow, 3=reflecting, 4=user, 5=diode) +OPT__BC_FLU_YM 1 # fluid boundary condition at the -y face: (1=periodic, 2=outflow, 3=reflecting, 4=user, 5=diode) +OPT__BC_FLU_YP 1 # fluid boundary condition at the +y face: (1=periodic, 2=outflow, 3=reflecting, 4=user, 5=diode) +OPT__BC_FLU_ZM 1 # fluid boundary condition at the -z face: (1=periodic, 2=outflow, 3=reflecting, 4=user, 5=diode) +OPT__BC_FLU_ZP 1 # fluid boundary condition at the +z face: (1=periodic, 2=outflow, 3=reflecting, 4=user, 5=diode) +OPT__BC_POT 1 # gravity boundary condition: (1=periodic, 2=isolated) + + +# time-step +DT__MAX 1.0 # dt criterion: maximum allowed dt (<0=off) [-1.0] +DT__FLUID 0.5 # 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 0.35355 # dt criterion: gravity solver safety factor (<0=auto) [-1.0] +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.0 # dt criterion: particle acceleration safety factor (0=off) [0.5] ##STORE_PAR_ACC ONLY## +DT__GRACKLE_COOLING 0.25 # dt criterion: factor for Grackle cooling time (<0=off) [-1.0] ##SUPPORT_GRACKLE only## +DT__SYNC_PARENT_LV 0.05 # 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 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## +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 1 # refine every REGRID_COUNT sub-step [4] +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 2 # 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_COOLING_LEN 1 # flag: cooling length (Input__Flag_CoolingLen) [0] ##HYDRO+SUPPORT_GRACKLE 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_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 0 # minimize MPI barriers to improve load balance, especially with particles [0] + # (STORE_POT_GHOST, PAR_IMPROVE_ACC=1, OPT__TIMING_BARRIER=0 only; recommend AUTO_REDUCE_DT=0) + + +# 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 # map to "with_radiative_cooling" [1] +GRACKLE_PRIMORDIAL 3 # map to "primordial_chemistry" (0=Cloudy, 1/2/3=6-/9-/12-species) (must increase NCOMP_PASSIVE_USER by 6/9/12, respectively) [0] +GRACKLE_METAL 1 # map to "metal_cooling" (must increase NCOMP_PASSIVE_USER by 1) [0] +GRACKLE_UV 1 # map to "UVbackground" [0] +GRACKLE_CMB_FLOOR 1 # map to "cmb_temperature_floor" [1] +GRACKLE_PE_HEATING 1 # map to "photoelectric_heating" [0] +GRACKLE_PE_HEATING_RATE 8.5e-26 # map to "photoelectric_heating_rate" (in erg/cm^3/s/n_H) [8.5e-26] +GRACKLE_CLOUDY_TABLE CloudyData_UVB=HM2012.h5 # "grackle_data_file" +GRACKLE_THREE_BODY_RATE 0 # map to "three_body_rate" (0=Abel+02, 1=Palla+83, 2=Cohen+83, 3=Flower+07, 4=Glover+08, 5=Forrey+13) [0] +GRACKLE_CIE_COOLING 0 # map to "cie_cooling" [0] +GRACKLE_H2_OPA_APPROX 0 # map to "h2_optical_depth_approximation" (0=off, 1=Ripomonti+04) [0] +GRACKLE_USE_V_HEATING_RATE 0 # map to "use_volumetric_heating_rate" [0] +GRACKLE_USE_S_HEATING_RATE 0 # map to "use_specific_heating_rate" [0] +GRACKLE_HYDROGEN_MFRAC 0.76 # map to "HydrogenFractionByMass" (NOT set to Grackle when GRACKLE_PRIMORDIAL==0) [0.76] +OPT__UNFREEZE_GRACKLE 1 # allow the evolution by Grackle solver when the fluid is frozen (for OPT__FREEZE_FLUID==1 only) [0] +CHE_GPU_NPGROUP -1 # number of patch groups sent into the CPU/GPU Grackle solver (<=0=auto) [-1] + + +# fluid solver in HYDRO (MODEL==HYDRO only) +GAMMA 1.666667 # ratio of specific heats (i.e., adiabatic index) [5.0/3.0] +MOLECULAR_WEIGHT 0.6 # mean molecular weight -> currently only for post-processing [0.6] +MINMOD_COEFF 1.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: + # (-1=auto, 0=none, 1=vanLeer, 2=generalized MinMod, 3=vanAlbada, 4=vanLeer+generalized MinMod, 6=central, 7=Athena) [-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 (<0=auto, 0=none, 1=Roe, 2=HLLC, 3=HLLE, 4=HLLD) [-1] +DUAL_ENERGY_SWITCH 2.0e-1 # apply dual-energy if E_int/E_kin < DUAL_ENERGY_SWITCH [2.0e-2] ##DUAL_ENERGY ONLY## + + +# fluid solvers in all models +FLU_GPU_NPGROUP -1 # number of patch groups sent into the 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_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 0 # ensure "sum(passive_scalar_density) == gas_density" [1] +OPT__RESET_FLUID 0 # reset fluid variables after each update -> edit "Flu_ResetByUser.cpp" [0] +OPT__FREEZE_FLUID 1 # do not evolve fluid at all [0] +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_ENTR 0.0 # minimum entropy (must >= 0.0) [0.0] ##HYDRO and MHD ONLY## +JEANS_MIN_PRES 1 # 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] + + +# initialization +OPT__INIT 1 # initialization option: (1=FUNCTION, 2=RESTART, 3=FILE->"UM_IC") +RESTART_LOAD_NRANK 256 # number of parallel I/O (i.e., number of MPI ranks) for restart [1] +OPT__INIT_RESTRICT 1 # restrict all data during the initialization [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__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__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] +MONO_MAX_ITER 10 # maximum number of iterations to reduce INT_MONO_COEFF when interpolation fails (0=off) [10] + + +# data dump +OPT__OUTPUT_TOTAL 1 # output the simulation snapshot: (0=off, 1=HDF5, 2=C-binary) [1] +OPT__OUTPUT_PART 1 # output a single line or slice: (0=off, 1=xy, 2=yz, 3=xz, 4=x, 5=y, 6=z, 7=diag, 8=entire box) [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_POT 1 # output gravitational potential [1] ##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_PRES 1 # output gas pressure [0] ##HYDRO ONLY## +OPT__OUTPUT_TEMP 1 # output gas temperature [0] ##HYDRO ONLY## +OPT__OUTPUT_ENTR 0 # output gas entropy [0] ##HYDRO ONLY## +OPT__OUTPUT_CS 1 # 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_LORENTZ 0 # output Lorentz factor [0] ##SRHD ONLY## +OPT__OUTPUT_3VELOCITY 0 # output 3-velocities [0] ##SRHD ONLY## +OPT__OUTPUT_GRACKLE_TEMP 1 # output the temperature calculated by Grackle [0] ##SUPPORT_GRACKLE only## +OPT__OUTPUT_GRACKLE_MU 1 # output the mean molecular weight calculated by Grackle [0] ##SUPPORT_GRACKLE only## +OPT__OUTPUT_GRACKLE_TCOOL 1 # output the cooling time calculated by Grackle [0] ##SUPPORT_GRACKLE only## +OPT__OUTPUT_USER_FIELD 0 # output user-defined derived fields [0] -> edit "Flu_DerivedField_User.cpp" +OPT__OUTPUT_MODE 1 # (1=const step, 2=const dt, 3=dump table) -> edit "Input__DumpTable" for 3 +OUTPUT_STEP 1 # output data every OUTPUT_STEP step ##OPT__OUTPUT_MODE==1 ONLY## +OUTPUT_DT 1.0 # output data every OUTPUT_DT time interval ##OPT__OUTPUT_MODE==2 ONLY## +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] + + +# 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_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 1 # 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 1 # 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 ##HYDRO and ELBDM ONLY## [0] +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] diff --git a/example/test_problem/Hydro/GrackleTest/Input__TestProb b/example/test_problem/Hydro/GrackleTest/Input__TestProb new file mode 100644 index 0000000000..a132b595a8 --- /dev/null +++ b/example/test_problem/Hydro/GrackleTest/Input__TestProb @@ -0,0 +1,21 @@ +# problem-specific runtime parameters +GrackleTest_DefaultTestMode 0 # Default test mode (0=none, >=1: see README) [0] +GrackleTest_MassDensity_Min 1.0e-29 # Minimum total mass density in the box (in g cm^-3) [1.0e-29] +GrackleTest_MassDensity_Max 1.0e-21 # Maximum total mass density in the box (in g cm^-3) [1.0e-21] +GrackleTest_TempOverMMW_Min 1.0e+00 # Minimum temperature over mean molecular weight (T/mu) in the box (in K) [1.0e+00] +GrackleTest_TempOverMMW_Max 1.0e+08 # Maximum temperature over mean molecular weight (T/mu) in the box (in K) [1.0e+08] +GrackleTest_MFrac_Metal 0.01295 # Metal mass fraction (GRACKLE_METAL only) [0.01295] +GrackleTest_MFrac_e 0.0 # Electron mass fraction (GRACKLE_PRIMORDIAL >= 1 only) [0.0] +GrackleTest_MFrac_HI 0.750158 # HI mass fraction (GRACKLE_PRIMORDIAL >= 1 only) [0.750158] +GrackleTest_MFrac_HII 0.0 # HII mass fraction (GRACKLE_PRIMORDIAL >= 1 only) [0.0] +GrackleTest_MFrac_HeI 0.236892 # HeI mass fraction (GRACKLE_PRIMORDIAL >= 1 only) [0.236892] +GrackleTest_MFrac_HeII 0.0 # HeII mass fraction (GRACKLE_PRIMORDIAL >= 1 only) [0.0] +GrackleTest_MFrac_HeIII 0.0 # HeIII mass fraction (GRACKLE_PRIMORDIAL >= 1 only) [0.0] +GrackleTest_MFrac_HM 0.0 # HM mass fraction (GRACKLE_PRIMORDIAL >= 2 only) [0.0] +GrackleTest_MFrac_H2I 0.0 # H2I mass fraction (GRACKLE_PRIMORDIAL >= 2 only) [0.0] +GrackleTest_MFrac_H2II 0.0 # H2II mass fraction (GRACKLE_PRIMORDIAL >= 2 only) [0.0] +GrackleTest_MFrac_DI 0.0 # DI mass fraction (GRACKLE_PRIMORDIAL >= 3 only) [0.0] +GrackleTest_MFrac_DII 0.0 # DII mass fraction (GRACKLE_PRIMORDIAL >= 3 only) [0.0] +GrackleTest_MFrac_HDI 0.0 # HDI mass fraction (GRACKLE_PRIMORDIAL >= 3 only) [0.0] +GrackleTest_HeatingRate 0.0 # User-provided heating rate (in erg cm^-3 s^-1 n_H^-1) [0.0] +GrackleTest_CoolingRate 0.0 # User-provided cooling rate (in erg cm^-3 s^-1 n_H^-2) [0.0] diff --git a/example/test_problem/Hydro/GrackleTest/README b/example/test_problem/Hydro/GrackleTest/README new file mode 100644 index 0000000000..94e11b9f34 --- /dev/null +++ b/example/test_problem/Hydro/GrackleTest/README @@ -0,0 +1,119 @@ +Compilation flags: +======================================== +Enable : MODEL=HYDRO, SUPPORT_GRACKLE +Disable: COMOVING, PARTICLE, GRAVITY + + +Default setup: +======================================== +1. Units: + (1) External (for Input__TestProb only): + See the comment of each runtime variable + + (2) Internal (for all other input files and internal usage): + [L] = kpc + [M] = 1.0e9 Msun + [T] = Myr + -->[V] = [L]/[T] ~ 9.8e2 km/s + -->[D] = [M]/[L]^3 ~ 6.8e-23 g/cm^3 + +2. BOX_SIZE 8.0 kpc + NX0_TOT_X 128 + NX0_TOT_Y 128 + NX0_TOT_Z 32 + +3. MAX_LEVEL 0 + OPT__FLAG_COOLING_LEN 1 + REGRID_COUNT 1 + +4. DT__MAX 1.00 Myr + DT__GRACKLE_COOLING 0.25 + +5. GRACKLE_ACTIVATE 1 + OPT__FREEZE_FLUID 1 + OPT__UNFREEZE_GRACKLE 1 + +6. OPT__OUTPUT_TEMP 1 + OPT__OUTPUT_GRACKLE_TEMP 1 + OPT__OUTPUT_GRACKLE_MU 1 + OPT__OUTPUT_GRACKLE_TCOOL 1 + OUTPUT_STEP 1 + +7. END_STEP 10 Steps + END_T 10.0 Myr or 10 * DT__GRACKLE_COOLING * cooling time + + +Note: +======================================== +1. Test the functionality and the interface to the Grackle library + +2. The initial conditions of varying phases + - In the x direction, linearly varying log( mass density ) + from `GrackleTest_MassDensity_Min` to `GrackleTest_MassDensity_Max` + - In the y direction, linearly varying log( temperature / mean molecular weight ) + from `GrackleTest_TempOverMMW_Min to `GrackleTest_TempOverMMW_Max` + - No variation in the z direction + +3. The chemical composition + - During the compilation, one must set `NCOMP_PASSIVE_USER` (`--passive` in `generate_make.sh`) + to the total number of chemical species (e.g., 13 for metal + 12 primordial species) + - One must turn on the corresponding `GRACKLE_METAL` and `GRACKLE_PRIMORDIAL` in `Input__Parameter` + - The mass fraction of each species can be set by `GrackleTest_MFrac_*` in `Input__TestProb` + - Currently, the mass fractions are assumed to be spatially uniform initially + --> Modify SetGridIC() if they are not + +4. The user-provided heating and cooling + - Turn on `GRACKLE_USE_V_HEATING_RATE` in `Input__Parameter` + - Set `GrackleTest_HeatingRate` and `GrackleTest_CoolingRate` in `Input__TestProb` + - Edit `Grackle_vHeatingRate_GrackleTest()` to set the distribution + +5. The default test mode can be set by `GrackleTest_DefaultTestMode` in `Input__TestProb` + - 0: + --> Non-default, one can choose the parameters freely + - 1: + - The density range is from 1.0e-29 g cm^-3 to 1.0e-21 g cm^-3 + - The T/mu range is from 1 K to 1e8 K + --> To see the cooling rates across the phase diagram + - 2: + - The density is uniformly 1.0e-24 g cm^-3 + - The T/mu is uniformly 1e4 K + --> To see the time evolution of temperature and composition for a given initial phase + - 3: + - The density is uniformly 1.0e-28 g cm^-3 + - The T/mu is uniformly 1e6 K + - The metal mass fraction is 0 + - `GRACKLE_PRIMORDIAL` should be 0 + - `GRACKLE_UV`, `GRACKLE_CMB_FLOOR`, `GRACKLE_PE_HEATING` should all be 0 + - `GRACKLE_USE_V_HEATING_RATE` should be 1 + - `GrackleTest_HeatingRate` = 1.0e-24 erg cm^-3 s^-1 n_H^-1 + - `GrackleTest_CoolingRate` = 1.6e-20 erg cm^-3 s^-1 n_H^-2 + --> To see the time evolution and spatial distribution of temperature for the given + user-provided one heating and one cooling 2D Gaussian distributions (as a demo) + +6. The evolution of fluid + - The fluid is not evolved by hydrodynamics (`OPT__FREEZE_FLUID`) + - But it is evolved by Grackle solver (`OPT__UNFREEZE_GRACKLE`) + +7. The default time step is set to min( 0.25x the minimum cooling time, 0.1 Myr ) + - Set `DT__GRACKLE_COOLING` and `DT__MAX` to change it + +8. Output the Grackle-calculated fields in the snapshots for analysis + - Turn on `OPT__OUTPUT_GRACKLE_TEMP`, `OPT__OUTPUT_GRACKLE_MU`, and `OPT__OUTPUT_GRACKLE_TCOOL` + +9. Refine the grids by the cooling length (= cooling time x sound speed) + - Turn on `OPT__FLAG_COOLING_LEN` + - Set the critera in `Input__Flag_CoolingLen` + - Set the target `MAX_LEVEL` + +10. Some handy yt analysis scripts are put at "yt_script" + - `plot_projection.py`: to plot the spatial distribution of gas quantities + - `plot_time_evolution.py`: to plot the time evolution of gas quantities + +11. CloudyData_UVB table + - Source: https://github.com/grackle-project/grackle_data_files/tree/main/input + - The default table `CloudyData_UVB=HM2012.h5` can be downloaded by `sh download_CloudyData_UVB.sh` + + +First-time GRACKLE installation guide (on NTU clusters as an example): +======================================== +Please refer to the README file at "./example/grackle/" diff --git a/example/test_problem/Hydro/GrackleTest/clean.sh b/example/test_problem/Hydro/GrackleTest/clean.sh new file mode 100644 index 0000000000..6a330bb25b --- /dev/null +++ b/example/test_problem/Hydro/GrackleTest/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 diff --git a/example/test_problem/Hydro/GrackleTest/download_CloudyData_UVB.sh b/example/test_problem/Hydro/GrackleTest/download_CloudyData_UVB.sh new file mode 100644 index 0000000000..4edd38919e --- /dev/null +++ b/example/test_problem/Hydro/GrackleTest/download_CloudyData_UVB.sh @@ -0,0 +1,10 @@ +#!/bin/bash +LOCAL_FILENAME="CloudyData_UVB=HM2012.h5" +FILE_ID="677ca211999605c485c8de6c" +FILE_SHA256="8715f1b39e90a7296ec2adcd442fa13a3d45d2ad021c6fa2fae9e4ab7a4700b2" + +# file download +curl -L https://hub.yt/api/v1/item/${FILE_ID}/download -o "${LOCAL_FILENAME}" + +# compare sha256sum +! [ `sha256sum ${LOCAL_FILENAME} | awk '{print $1}'` = "${FILE_SHA256}" ] && echo "File broken: ${LOCAL_FILENAME}" diff --git a/example/test_problem/Hydro/GrackleTest/generate_make.sh b/example/test_problem/Hydro/GrackleTest/generate_make.sh new file mode 100644 index 0000000000..bc933ccf2d --- /dev/null +++ b/example/test_problem/Hydro/GrackleTest/generate_make.sh @@ -0,0 +1,5 @@ +# This script should run in the same directory as configure.py + +PYTHON=python3 + +${PYTHON} configure.py --mpi=true --hdf5=true --gpu=true --model=HYDRO --grackle=true --passive=13 "$@" diff --git a/example/test_problem/Hydro/GrackleTest/yt_script/plot_projection.py b/example/test_problem/Hydro/GrackleTest/yt_script/plot_projection.py new file mode 100644 index 0000000000..b4eda9629e --- /dev/null +++ b/example/test_problem/Hydro/GrackleTest/yt_script/plot_projection.py @@ -0,0 +1,124 @@ +import argparse +import sys +import yt +import numpy as np + +# load the command-line parameters +parser = argparse.ArgumentParser( description='Plot the gas projections' ) + +parser.add_argument( '-p', action='store', required=False, type=str, dest='prefix', + help='path prefix [%(default)s]', default='../' ) +parser.add_argument( '-s', action='store', required=True, type=int, dest='idx_start', + help='first data index' ) +parser.add_argument( '-e', action='store', required=True, type=int, dest='idx_end', + help='last data index' ) +parser.add_argument( '-d', action='store', required=False, type=int, dest='didx', + help='delta data index [%(default)d]', default=1 ) + +args=parser.parse_args() + +# take note +print( '\nCommand-line arguments:' ) +print( '-------------------------------------------------------------------' ) +print( ' '.join(map(str, sys.argv)) ) +print( '-------------------------------------------------------------------\n' ) + + +idx_start = args.idx_start +idx_end = args.idx_end +didx = args.didx +prefix = args.prefix + + +yt.enable_parallelism() + +ts = yt.DatasetSeries( [ prefix+'/Data_%06d'%idx for idx in range(idx_start, idx_end+1, didx) ] ) + + +# target fields and their colorbar limits +zlim = { + 'density' : ( 1.0e-29, 1.0e-21 ), + 'temperature' : ( 1.0e+00, 1.0e+08 ), + 'thermal_energy_density' : ( 1.0e-21, 1.0e-05 ), + 'sound_speed' : ( 1.0e+04, 1.0e+08 ), + 'electron_density' : ( 1.0e-29, 1.0e-21 ), + 'DI_density' : ( 1.0e-29, 1.0e-21 ), + 'DII_density' : ( 1.0e-29, 1.0e-21 ), + 'H2I_density' : ( 1.0e-29, 1.0e-21 ), + 'H2II_density' : ( 1.0e-29, 1.0e-21 ), + 'HDI_density' : ( 1.0e-29, 1.0e-21 ), + 'HI_density' : ( 1.0e-29, 1.0e-21 ), + 'HII_density' : ( 1.0e-29, 1.0e-21 ), + 'HM_density' : ( 1.0e-29, 1.0e-21 ), + 'HeI_density' : ( 1.0e-29, 1.0e-21 ), + 'HeII_density' : ( 1.0e-29, 1.0e-21 ), + 'HeIII_density' : ( 1.0e-29, 1.0e-21 ), + 'metal_density' : ( 1.0e-29, 1.0e-21 ), + 'grackle_T_over_mu' : ( 1.0e+00, 1.0e+08 ), + 'grackle_temperature' : ( 1.0e+00, 1.0e+08 ), + 'grackle_mu' : ( 5.0e-01, 1.5e+00 ), + 'grackle_cooling_rate' : ( -1.0e-24, 1.0e-24 ), + 'grackle_cooling_strength' : ( -1.0e+03, 1.0e+03 ), + 'grackle_cooling_length' : ( 1.0e-04, 1.0e+00 ), + 'grackle_cooling_length_over_dh' : ( 5.0e-02, 2.0e+01 ), + } + + +def set_derived_fields( ds ): + + def _grackle_T_over_mu( field, data ): + return data['grackle_temperature'] / data['grackle_mu'] + ds.add_field( ('gas', 'grackle_T_over_mu'), function=_grackle_T_over_mu, sampling_type='cell', units='K', display_name='Grackle $T/\mu$' ) + + def _grackle_cooling_rate( field, data ): + return data['thermal_energy_density'] / data['grackle_cooling_time'] + ds.add_field( ('gas', 'grackle_cooling_rate'), function=_grackle_cooling_rate, sampling_type='cell', units='erg/cm**3/s' ) + + def _grackle_cooling_strength( field, data ): + return 1.0 / data['grackle_cooling_time'] + ds.add_field( ('gas', 'grackle_cooling_strength'), function=_grackle_cooling_strength, sampling_type='cell', units='Myr**-1' ) + + def _grackle_cooling_length( field, data ): + factor = np.where( data['grackle_cooling_time'] < 0.0, 1.0, np.inf ) + return factor * np.abs( data['grackle_cooling_time'] ) * data['sound_speed'] + ds.add_field( ('gas', 'grackle_cooling_length'), function=_grackle_cooling_length, sampling_type='cell', units='kpc' ) + + def _grackle_cooling_length_over_dh( field, data ): + return data['grackle_cooling_length'] / data['dx'] + ds.add_field( ('gas', 'grackle_cooling_length_over_dh'), function=_grackle_cooling_length_over_dh, sampling_type='cell', units='dimensionless', display_name='Grackle Cooling Length / dh' ) + + +for ds in ts.piter(): + + set_derived_fields( ds ) + + ad = ds.all_data() + + # plot the projections + for field in zlim.keys(): + + if ds.parameters["Grackle_Primordial"] < 3 and (field == 'DI_density' or field == 'DII_density' or field == 'HDI_density'): + continue + if ds.parameters["Grackle_Primordial"] < 2 and (field == 'H2I_density' or field == 'H2II_density' or field == 'HM_density'): + continue + if ds.parameters["Grackle_Primordial"] < 1 and (field == 'HI_density' or field == 'HII_density' or field == 'HeI_density' or field == 'HeII_density' or field == 'HeIII_density' or field == 'electron_density'): + continue + if ds.parameters["Grackle_Metal"] == 0 and field == 'metal_density': + continue + + # Projection weighted by 1 -> get the average + pz = yt.ProjectionPlot( ds, 'z', field, weight_field=('index','ones'), center='c', data_source=ad ) + pz.set_axes_unit( 'kpc' ) + pz.set_zlim( field, zlim[field][0], zlim[field][1] ) + + pz.annotate_text( (0.1,0.9), + "min = {: >11.3e}\nmax = {: >11.3e}".format(ad[field].min(),ad[field].max()), + coord_system='axis', text_args={"color": "grey"} ) + + pz.annotate_timestamp( time_unit='Myr', corner='upper_right' ) + + pz.save( 'fig_%s_Projection_z_%s.png'%(ds, field), mpl_kwargs={'dpi':150} ) + + if field == 'grackle_cooling_length' or field == 'grackle_cooling_length_over_dh': + pz.annotate_grids() + pz.save( 'fig_%s_Projection_z_%s_withgrids.png'%(ds, field), mpl_kwargs={'dpi':150} ) diff --git a/example/test_problem/Hydro/GrackleTest/yt_script/plot_time_evolution.py b/example/test_problem/Hydro/GrackleTest/yt_script/plot_time_evolution.py new file mode 100644 index 0000000000..afac47176c --- /dev/null +++ b/example/test_problem/Hydro/GrackleTest/yt_script/plot_time_evolution.py @@ -0,0 +1,161 @@ +import argparse +import sys +import yt +import matplotlib.pyplot as plt +import numpy as np + +# load the command-line parameters +parser = argparse.ArgumentParser( description='Plot the gas time evolution' ) + +parser.add_argument( '-p', action='store', required=False, type=str, dest='prefix', + help='path prefix [%(default)s]', default='../' ) +parser.add_argument( '-s', action='store', required=True, type=int, dest='idx_start', + help='first data index' ) +parser.add_argument( '-e', action='store', required=True, type=int, dest='idx_end', + help='last data index' ) +parser.add_argument( '-d', action='store', required=False, type=int, dest='didx', + help='delta data index [%(default)d]', default=1 ) + +args=parser.parse_args() + +# take note +print( '\nCommand-line arguments:' ) +print( '-------------------------------------------------------------------' ) +print( ' '.join(map(str, sys.argv)) ) +print( '-------------------------------------------------------------------\n' ) + + +idx_start = args.idx_start +idx_end = args.idx_end +didx = args.didx +prefix = args.prefix + + +yt.enable_parallelism() + +ts = yt.DatasetSeries( [ prefix+'/Data_%06d'%idx for idx in range(idx_start, idx_end+1, didx) ] ) + + +time = np.array([]) +density = np.array([]) +metal_density = np.array([]) +electron_density = np.array([]) +HI_density = np.array([]) +HII_density = np.array([]) +HeI_density = np.array([]) +HeII_density = np.array([]) +HeIII_density = np.array([]) +HM_density = np.array([]) +H2I_density = np.array([]) +H2II_density = np.array([]) +DI_density = np.array([]) +DII_density = np.array([]) +HDI_density = np.array([]) +temperature = np.array([]) +gr_temperature = np.array([]) +gr_mu = np.array([]) +temperature_max = np.array([]) +temperature_min = np.array([]) +gr_temperature_max = np.array([]) +gr_temperature_min = np.array([]) + + +# collect the results +for ds in ts.piter(): + + time = np.append( time , ds.current_time.in_units('Myr').d ) + density = np.append( density , np.average(ds.all_data()['density' ].in_units('g/cm**3').d) ) + + if ds.parameters['Grackle_Metal'] == 1: + metal_density = np.append( metal_density , np.average(ds.all_data()['metal_density' ].in_units('g/cm**3').d) ) + + if ds.parameters['Grackle_Primordial'] >= 1: + electron_density = np.append( electron_density , np.average(ds.all_data()['electron_density' ].in_units('g/cm**3').d) ) + HI_density = np.append( HI_density , np.average(ds.all_data()['HI_density' ].in_units('g/cm**3').d) ) + HII_density = np.append( HII_density , np.average(ds.all_data()['HII_density' ].in_units('g/cm**3').d) ) + HeI_density = np.append( HeI_density , np.average(ds.all_data()['HeI_density' ].in_units('g/cm**3').d) ) + HeII_density = np.append( HeII_density , np.average(ds.all_data()['HeII_density' ].in_units('g/cm**3').d) ) + HeIII_density = np.append( HeIII_density , np.average(ds.all_data()['HeIII_density' ].in_units('g/cm**3').d) ) + + if ds.parameters['Grackle_Primordial'] >= 2: + HM_density = np.append( HM_density , np.average(ds.all_data()['HM_density' ].in_units('g/cm**3').d) ) + H2I_density = np.append( H2I_density , np.average(ds.all_data()['H2I_density' ].in_units('g/cm**3').d) ) + H2II_density = np.append( H2II_density , np.average(ds.all_data()['H2II_density' ].in_units('g/cm**3').d) ) + + if ds.parameters['Grackle_Primordial'] >= 3: + DI_density = np.append( DI_density , np.average(ds.all_data()['DI_density' ].in_units('g/cm**3').d) ) + DII_density = np.append( DII_density , np.average(ds.all_data()['DII_density' ].in_units('g/cm**3').d) ) + HDI_density = np.append( HDI_density , np.average(ds.all_data()['HDI_density' ].in_units('g/cm**3').d) ) + + temperature = np.append( temperature , np.average(ds.all_data()['temperature' ].in_units('K').d) ) + gr_temperature = np.append( gr_temperature , np.average(ds.all_data()['grackle_temperature'].in_units('K').d) ) + gr_mu = np.append( gr_mu , np.average(ds.all_data()['grackle_mu' ].d) ) + + temperature_max = np.append( temperature_max , np.max( ds.all_data()['temperature' ].in_units('K').d) ) + temperature_min = np.append( temperature_min , np.min( ds.all_data()['temperature' ].in_units('K').d) ) + gr_temperature_max = np.append( gr_temperature_max, np.max( ds.all_data()['grackle_temperature'].in_units('K').d) ) + gr_temperature_min = np.append( gr_temperature_min, np.min( ds.all_data()['grackle_temperature'].in_units('K').d) ) + + +# plot the results +f, ax = plt.subplots( 1, 2, figsize=(12.8, 4.8) ) + + +# plot the temperature evolution +ax[0].plot( time, gr_temperature, 'k-s', label='Grackle T' ) +ax[0].plot( time, temperature, 'k--', label='GAMER T' ) +ax[0].plot( time, gr_temperature/gr_mu, 'k:', label='Grackle T/mu' ) + +if ds.parameters['GrackleTest_DefaultTestMode'] == 3: + # plot the expected temperature evolution for the given heating and cooling rates + r_heat = ds.parameters['GrackleTest_HeatingRate'] * ds.parameters['Grackle_HydrogenMFrac'] * (ds.parameters['Gamma']-1) * ds.parameters['MolecularWeight'] / ds.units.boltzmann_constant.in_units('erg/K').v *(1.*ds.units.Myr).in_units('s').v + r_cool = (density[0]*ds.parameters['Grackle_HydrogenMFrac']/ds.units.hydrogen_mass.in_units('g').v*ds.parameters['GrackleTest_CoolingRate']) * ds.parameters['Grackle_HydrogenMFrac'] * (ds.parameters['Gamma']-1) * ds.parameters['MolecularWeight'] / ds.units.boltzmann_constant.in_units('erg/K').v *(1.*ds.units.Myr).in_units('s').v + ax[0].plot( time, temperature_max[0]+time*r_heat, 'g-.', label='Analytical max(T)' ) + ax[0].plot( time, temperature_min[0]-time*r_cool, 'y-.', label='Analytical min(T)' ) + + # plot the maximum and minimum temperature evolution as the results of the given heating and cooling rates + ax[0].plot( time, temperature_max, 'r--', label='max(T)' ) + ax[0].plot( time, temperature_min, 'b--', label='min(T)' ) + +ax[0].set_yscale( 'log' ) +ax[0].set_xlim( 0.0, 1.33*np.max(time) ) +ax[0].set_ylim( 0.1*np.min(gr_temperature_min), 10.0*np.max(gr_temperature_max) ) +ax[0].set_xlabel( 'Time (Myr)', fontsize='large' ) +ax[0].set_ylabel( 'Temperature (K)', fontsize='large' ) +ax[0].legend() + + +# plot the density evolution +ax[1].plot( time, density, 'k-s', label='' ) + +if ds.parameters['Grackle_Metal'] == 1: + ax[1].plot( time, metal_density, 'g-', label='Metal' ) + +if ds.parameters['Grackle_Primordial'] >= 1: + ax[1].plot( time, electron_density, 'g--', label='e' ) + ax[1].plot( time, HI_density, 'r-', label='HI' ) + ax[1].plot( time, HII_density, 'r--', label='HII' ) + ax[1].plot( time, HeI_density, 'b-', label='HeI' ) + ax[1].plot( time, HeII_density, 'b--', label='HeII' ) + ax[1].plot( time, HeIII_density, 'b-.', label='HEIII' ) + +if ds.parameters['Grackle_Primordial'] >= 2: + ax[1].plot( time, HM_density, 'r:', label='HM' ) + ax[1].plot( time, H2I_density, 'm--', label='H2I' ) + ax[1].plot( time, H2II_density, 'm-', label='H2II' ) + +if ds.parameters['Grackle_Primordial'] >= 3: + ax[1].plot( time, DI_density, 'c-', label='DI' ) + ax[1].plot( time, DII_density, 'c--', label='DII' ) + ax[1].plot( time, HDI_density, 'y-', label='HDI' ) + +ax[1].set_yscale( 'log' ) +ax[1].set_xlim( 0.0, 1.33*np.max(time) ) +ax[1].set_ylim( 3.0e-6*np.min(density), 3.0*np.max(density) ) +ax[1].set_xlabel( 'Time (Myr)', fontsize='large' ) +ax[1].set_ylabel( 'Density (g/cm$^3$)', fontsize='large' ) +ax[1].legend() + + +f.subplots_adjust( wspace=0.2 ) +f.savefig( 'fig_time_evolution.png', bbox_inches='tight', pad_inches=0.05, dpi=150 ) diff --git a/example/test_problem/Hydro/GrackleTest/yt_script/yt.toml b/example/test_problem/Hydro/GrackleTest/yt_script/yt.toml new file mode 100644 index 0000000000..db8ca3d31d --- /dev/null +++ b/example/test_problem/Hydro/GrackleTest/yt_script/yt.toml @@ -0,0 +1,201 @@ +[fields.gamer.GrackleTemp] +aliases = ["grackle_temperature"] +units = "code_temperature" +display_name = "Grackle Temperature" + +[fields.gamer.GrackleMu] +aliases = ["grackle_mu"] +units = "dimensionless" +display_name = "Grackle Mean Molecular Weight" + +[fields.gamer.GrackleTCool] +aliases = ["grackle_cooling_time"] +units = "code_time" +display_name = "Grackle Cooling Time" + +[fields.gamer.Metal] +aliases = ["metal_density"] +units = "code_density" +display_name = "Metal Mass Density" + +[fields.gamer.Electron] +aliases = ["electron_density"] +units = "code_density" +display_name = "Electron Mass Density" + +[fields.gamer.HI] +aliases = ["HI_density"] +units = "code_density" +display_name = "HI Mass Density" + +[fields.gamer.HII] +aliases = ["HII_density"] +units = "code_density" +display_name = "HII Mass Density" + +[fields.gamer.HeI] +aliases = ["HeI_density"] +units = "code_density" +display_name = "HeI Mass Density" + +[fields.gamer.HeII] +aliases = ["HeII_density"] +units = "code_density" +display_name = "HeII Mass Density" + +[fields.gamer.HeIII] +aliases = ["HeIII_density"] +units = "code_density" +display_name = "HeIII Mass Density" + +[fields.gamer.HM] +aliases = ["HM_density"] +units = "code_density" +display_name = "HM Mass Density" + +[fields.gamer.H2I] +aliases = ["H2I_density"] +units = "code_density" +display_name = "H2I Mass Density" + +[fields.gamer.H2II] +aliases = ["H2II_density"] +units = "code_density" +display_name = "H2II Mass Density" + +[fields.gamer.DI] +aliases = ["DI_density"] +units = "code_density" +display_name = "DI Mass Density" + +[fields.gamer.DII] +aliases = ["DII_density"] +units = "code_density" +display_name = "DII Mass Density" + +[fields.gamer.HDI] +aliases = ["HDI_density"] +units = "code_density" +display_name = "HDI Mass Density" + +[plot.gas.density] +cmap = "viridis" +log = true +units = "g/cm**3" + +[plot.gas.metal_density] +cmap = "viridis" +log = true +units = "g/cm**3" + +[plot.gas.HI_density] +cmap = "viridis" +log = true +units = "g/cm**3" + +[plot.gas.HII_density] +cmap = "viridis" +log = true +units = "g/cm**3" + +[plot.gas.HeI_density] +cmap = "viridis" +log = true +units = "g/cm**3" + +[plot.gas.HeII_density] +cmap = "viridis" +log = true +units = "g/cm**3" + +[plot.gas.HeIII_density] +cmap = "viridis" +log = true +units = "g/cm**3" + +[plot.gas.electron_density] +cmap = "viridis" +log = true +units = "g/cm**3" + +[plot.gas.H2I_density] +cmap = "viridis" +log = true +units = "g/cm**3" + +[plot.gas.H2II_density] +cmap = "viridis" +log = true +units = "g/cm**3" + +[plot.gas.HM_density] +cmap = "viridis" +log = true +units = "g/cm**3" + +[plot.gas.DI_density] +cmap = "viridis" +log = true +units = "g/cm**3" + +[plot.gas.DII_density] +cmap = "viridis" +log = true +units = "g/cm**3" + +[plot.gas.HDI_density] +cmap = "viridis" +log = true +units = "g/cm**3" + +[plot.gas.temperature] +cmap = "magma" +log = true +units = "K" + +[plot.gas.grackle_temperature] +cmap = "magma" +log = true +units = "K" + +[plot.gas.grackle_T_over_mu] +cmap = "magma" +log = true +units = "K" + +[plot.gas.thermal_energy_density] +cmap = "magma" +log = true +units = "erg/cm**3" + +[plot.gas.sound_speed] +cmap = "magma" +log = true +units = "cm/s" + +[plot.gas.grackle_mu] +cmap = "cividis" +log = false +units = "dimensionless" + +[plot.gas.grackle_cooling_rate] +cmap = "bwr" +log = true +linthresh=1.0e-29 +units = "erg/cm**3/s" + +[plot.gas.grackle_cooling_strength] +cmap = "bwr" +log = true +linthresh=1.0e-3 +units = "Myr**-1" + +[plot.gas.grackle_cooling_length] +cmap = "Blues_r" +log = true +units = "kpc" + +[plot.gas.grackle_cooling_length_over_dh] +cmap = "PiYG" +log = true +units = "dimensionless" diff --git a/example/test_problem/Template/Input__Parameter b/example/test_problem/Template/Input__Parameter index d7742ab30f..eb3c429009 100644 --- a/example/test_problem/Template/Input__Parameter +++ b/example/test_problem/Template/Input__Parameter @@ -47,6 +47,7 @@ TESTPROB_ID 0 # test problem ID [0] # 20: HYDRO MHD Cosmic Ray Soundwave # 21: HYDRO MHD Cosmic Ray Shocktube # 23: HYDRO MHD Cosmic Ray Diffusion + # 24: HYDRO Grackle Test (+SUPPORT_GRACKLE) # 100: HYDRO CDM cosmological simulation (+GRAVITY & COMOVING & PARTICLE) # 101: HYDRO Zeldovich pancake collapse (+GRAVITY & COMOVING & PARTICLE) # 1000: ELBDM external potential (+GRAVITY) @@ -121,6 +122,7 @@ DT__HYBRID_VELOCITY_INIT -1.0 # dt criterion: DT__HYBRID_VELOCITY in 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__GRACKLE_COOLING -1.0 # dt criterion: factor for Grackle cooling time (<0=off) [-1.0] ##SUPPORT_GRACKLE only## DT__CR_DIFFUSION 0.3 # dt criterion: cosmic-ray diffusion CFL factor [0.3] ##CR_DIFFUSION 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 @@ -152,6 +154,7 @@ OPT__FLAG_PRES_GRADIENT 0 # flag: pressure gradient (Input__Flag OPT__FLAG_LRTZ_GRADIENT 0 # flag: Lorentz factor gradient (Input__Flag_LrtzGradient) [0] ##SRHD 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_COOLING_LEN 0 # flag: cooling length (Input__Flag_CoolingLen) [0] ##HYDRO+SUPPORT_GRACKLE ONLY## OPT__FLAG_CURRENT 0 # flag: current density in MHD (Input__Flag_Current) [0] ##MHD ONLY## OPT__FLAG_CRAY 0 # flag: cosmic-ray energy (Input__Flag_CRay) [0] ##COSMIC_RAY ONLY## OPT__FLAG_ENGY_DENSITY 0 # flag: energy density (Input_Flag_EngyDensity) [0] ##ELBDM ONLY## @@ -205,17 +208,20 @@ SRC_GPU_NPGROUP -1 # number of patch groups sent into the GRACKLE_ACTIVATE 1 # enable Grackle [1] GRACKLE_VERBOSE 1 # map to "grackle_verbose" [1] GRACKLE_COOLING 1 # map to "with_radiative_cooling" [1] -GRACKLE_PRIMORDIAL 0 # map to "primordial_chemistry" (0=Cloudy, 1/2/3=6-/9-/12-species) [0] - # (must increase NCOMP_PASSIVE_USER by 6/9/12, respectively) +GRACKLE_PRIMORDIAL 0 # map to "primordial_chemistry" (0=Cloudy, 1/2/3=6-/9-/12-species) (must increase NCOMP_PASSIVE_USER by 6/9/12, respectively) [0] GRACKLE_METAL 0 # map to "metal_cooling" (must increase NCOMP_PASSIVE_USER by 1) [0] GRACKLE_UV 0 # map to "UVbackground" [0] GRACKLE_CMB_FLOOR 1 # map to "cmb_temperature_floor" [1] GRACKLE_PE_HEATING 0 # map to "photoelectric_heating" [0] -GRACKLE_PE_HEATING_RATE 8.5e-26 # map to "photoelectric_heating_rate (in erg/cm^3/s)" [8.5e-26] +GRACKLE_PE_HEATING_RATE 8.5e-26 # map to "photoelectric_heating_rate" (in erg/cm^3/s/n_H) [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 +GRACKLE_THREE_BODY_RATE 0 # map to "three_body_rate" (0=Abel+02, 1=Palla+83, 2=Cohen+83, 3=Flower+07, 4=Glover+08, 5=Forrey+13) [0] +GRACKLE_CIE_COOLING 0 # map to "cie_cooling" [0] +GRACKLE_H2_OPA_APPROX 0 # map to "h2_optical_depth_approximation" (0=off, 1=Ripomonti+04) [0] +GRACKLE_USE_V_HEATING_RATE 0 # map to "use_volumetric_heating_rate" [0] +GRACKLE_USE_S_HEATING_RATE 0 # map to "use_specific_heating_rate" [0] +GRACKLE_HYDROGEN_MFRAC 0.76 # map to "HydrogenFractionByMass" (NOT set to Grackle when GRACKLE_PRIMORDIAL==0) [0.76] +OPT__UNFREEZE_GRACKLE 0 # allow the evolution by Grackle solver when the fluid is frozen (for OPT__FREEZE_FLUID==1 only) [0] CHE_GPU_NPGROUP -1 # number of patch groups sent into the CPU/GPU Grackle solver (<=0=auto) [-1] @@ -401,6 +407,9 @@ OPT__OUTPUT_DIVMAG 0 # output |divergence(B)*dh/|B|| [0] ## OPT__OUTPUT_LORENTZ 0 # output Lorentz factor [0] ##SRHD ONLY## OPT__OUTPUT_3VELOCITY 0 # output 3-velocities [0] ##SRHD ONLY## OPT__OUTPUT_ENTHALPY 1 # output reduced enthalpy [1] ##SRHD ONLY## +OPT__OUTPUT_GRACKLE_TEMP 0 # output the temperature calculated by Grackle [0] ##SUPPORT_GRACKLE only## +OPT__OUTPUT_GRACKLE_MU 0 # output the mean molecular weight calculated by Grackle [0] ##SUPPORT_GRACKLE only## +OPT__OUTPUT_GRACKLE_TCOOL 0 # output the cooling time calculated by Grackle [0] ##SUPPORT_GRACKLE only## OPT__OUTPUT_USER_FIELD 0 # output user-defined derived fields [0] -> edit "Flu_DerivedField_User.cpp" OPT__OUTPUT_MODE 1 # (1=const step, 2=const dt, 3=dump table) -> edit "Input__DumpTable" for 3 OPT__OUTPUT_RESTART 0 # output data immediately after restart [0] diff --git a/include/Global.h b/include/Global.h index aabea2f5b6..c8e8834ce5 100644 --- a/include/Global.h +++ b/include/Global.h @@ -278,6 +278,16 @@ extern char GRACKLE_CLOUDY_TABLE[MAX_STRING]; extern int GRACKLE_THREE_BODY_RATE; extern bool GRACKLE_CIE_COOLING; extern int GRACKLE_H2_OPA_APPROX; +extern bool GRACKLE_USE_V_HEATING_RATE; +extern bool GRACKLE_USE_S_HEATING_RATE; +extern double GRACKLE_HYDROGEN_MFRAC; +extern bool OPT__UNFREEZE_GRACKLE; +extern bool OPT__OUTPUT_GRACKLE_TEMP; +extern bool OPT__OUTPUT_GRACKLE_MU; +extern bool OPT__OUTPUT_GRACKLE_TCOOL; +extern bool OPT__FLAG_COOLING_LEN; +extern double FlagTable_CoolingLen[NLEVEL-1]; +extern double DT__GRACKLE_COOLING; extern int CHE_GPU_NPGROUP; #endif diff --git a/include/HDF5_Typedef.h b/include/HDF5_Typedef.h index edde083e3b..aeb60bfd99 100644 --- a/include/HDF5_Typedef.h +++ b/include/HDF5_Typedef.h @@ -479,6 +479,9 @@ struct InputPara_t # ifdef CR_DIFFUSION double Dt__CR_Diffusion; # endif +# ifdef SUPPORT_GRACKLE + double Dt__GrackleCooling; +# endif # ifdef COMOVING double Dt__MaxDeltaA; # endif @@ -527,6 +530,9 @@ struct InputPara_t # ifdef SRHD int Opt__Flag_LrtzGradient; # endif +# ifdef SUPPORT_GRACKLE + int Opt__Flag_CoolingLen; +# endif # ifdef COSMIC_RAY int Opt__Flag_CRay; # endif @@ -699,6 +705,10 @@ struct InputPara_t int Grackle_ThreeBodyRate; int Grackle_CIE_Cooling; int Grackle_H2_OpaApprox; + int Grackle_UseVHeatingRate; + int Grackle_UseSHeatingRate; + double Grackle_HydrogenMFrac; + int Opt__UnfreezeGrackle; int Che_GPU_NPGroup; # endif @@ -831,6 +841,11 @@ struct InputPara_t int Opt__Output_3Velocity; int Opt__Output_Enthalpy; # endif +# ifdef SUPPORT_GRACKLE + int Opt__Output_GrackleTemp; + int Opt__Output_GrackleMu; + int Opt__Output_GrackleTCool; +# endif # endif // #if ( MODEL == HYDRO ) int Opt__Output_UserField; int Opt__Output_Mode; @@ -913,6 +928,9 @@ struct InputPara_t # ifdef SRHD double FlagTable_LrtzGradient[NLEVEL-1]; # endif +# ifdef SUPPORT_GRACKLE + double FlagTable_CoolingLen [NLEVEL-1]; +# endif # ifdef COSMIC_RAY double FlagTable_CRay [NLEVEL-1]; # endif diff --git a/include/Prototype.h b/include/Prototype.h index 767fd31d06..d7ca6d52cb 100644 --- a/include/Prototype.h +++ b/include/Prototype.h @@ -368,7 +368,8 @@ void FindFather( const int lv, const int Mode ); void Flag_Real( const int lv, const UseLBFunc_t UseLBFunc ); bool Flag_Check( const int lv, const int PID, const int i, const int j, const int k, const real dv, const real Fluid[][PS1][PS1][PS1], const real Pot[][PS1][PS1], const real MagCC[][PS1][PS1][PS1], - const real Vel[][PS1][PS1][PS1], const real Pres[][PS1][PS1], const real Lrtz[][PS1][PS1], + const real Vel[][PS1][PS1][PS1], const real Pres[][PS1][PS1], const real Lrtz[][PS1][PS1], + const real Lcool[][PS1][PS1], const real *Lohner_Var, const real *Lohner_Ave, const real *Lohner_Slope, const int Lohner_NVar, const real ParCount[][PS1][PS1], const real ParDens[][PS1][PS1], const real JeansCoeff, const real *Interf_Var, const real Spectral_Cond ); @@ -814,16 +815,18 @@ void CPU_SrcSolver( const real h_Flu_Array_In [][FLU_NIN_S ][ CUBE(SRC_NXT) // Grackle #ifdef SUPPORT_GRACKLE -void Grackle_Init(); -void Grackle_Init_FieldData(); -void Grackle_End(); -void Init_MemAllocate_Grackle( const int Che_NPG ); -void End_MemFree_Grackle(); -void Grackle_Prepare( const int lv, real_che h_Che_Array[], const int NPG, const int *PID0_List ); -void Grackle_Close( const int lv, const int SaveSg, const real_che h_Che_Array[], const int NPG, const int *PID0_List ); -void Grackle_AdvanceDt( const int lv, const double TimeNew, const double TimeOld, const double dt, const int SaveSg, - const bool OverlapMPI, const bool Overlap_Sync ); -void CPU_GrackleSolver( grackle_field_data *Che_FieldData, code_units Che_Units, const int NPatchGroup, const real dt ); +void Grackle_Init(); +void Grackle_Init_FieldData(); +void Grackle_End(); +void Init_MemAllocate_Grackle( const int Che_NPG ); +void End_MemFree_Grackle(); +void Grackle_Prepare( const int lv, real_che h_Che_Array[], const int NPG, const int *PID0_List ); +void Grackle_Close( const int lv, const int SaveSg, const real_che h_Che_Array[], const int NPG, const int *PID0_List ); +void Grackle_AdvanceDt( const int lv, const double TimeNew, const double TimeOld, const double dt, const int SaveSg, + const bool OverlapMPI, const bool Overlap_Sync ); +void CPU_GrackleSolver( grackle_field_data *Che_FieldData, code_units Che_Units, const int NPatchGroup, const real dt ); +void Grackle_Calculate( real Out[], const GrackleFieldBIdx_t TFields, const int lv, const int NPG, const int *PID0_List ); +double Grackle_GetTimeStep_CoolingTime( const int lv ); #endif // #ifdef SUPPORT_GRACKLE diff --git a/include/TestProb.h b/include/TestProb.h index e9d4a43afe..17b384c1ff 100644 --- a/include/TestProb.h +++ b/include/TestProb.h @@ -120,6 +120,10 @@ extern int (*FB_User_Ptr)( const int lv, const double TimeNew, const double Tim real (*Fluid)[FB_NXT][FB_NXT][FB_NXT], const double EdgeL[], const double dh, bool CoarseFine[], const int TID, RandomNumber_t *RNG ); #endif +#ifdef SUPPORT_GRACKLE +extern real_che (*Grackle_vHeatingRate_User_Ptr)( const double x, const double y, const double z, const double Time, const double n_H ); +extern real_che (*Grackle_sHeatingRate_User_Ptr)( const double x, const double y, const double z, const double Time ); +#endif diff --git a/include/Typedef.h b/include/Typedef.h index c2d767a47e..f6053aa92c 100644 --- a/include/Typedef.h +++ b/include/Typedef.h @@ -93,6 +93,7 @@ const TestProbID_t TESTPROB_HYDRO_CR_SOUNDWAVE = 20, TESTPROB_HYDRO_CR_SHOCKTUBE = 21, TESTPROB_HYDRO_CR_DIFFUSION = 23, + TESTPROB_HYDRO_GRACKLE_TEST = 24, TESTPROB_HYDRO_BARRED_POT = 51, TESTPROB_HYDRO_JET_ICM_WALL = 52, TESTPROB_HYDRO_CDM_LSS = 100, @@ -521,6 +522,13 @@ const GracklePriChe_t GRACKLE_PRI_CHE_NSPE6 = 1, GRACKLE_PRI_CHE_NSPE9 = 2, GRACKLE_PRI_CHE_NSPE12 = 3; + +// bitwise field indices used by Grackle +typedef long GrackleFieldBIdx_t; +const GrackleFieldBIdx_t + _GRACKLE_TEMP = ( 1L << 0 ), + _GRACKLE_MU = ( 1L << 1 ), + _GRACKLE_TCOOL = ( 1L << 2 ); #endif diff --git a/src/Auxiliary/Aux_Check_Parameter.cpp b/src/Auxiliary/Aux_Check_Parameter.cpp index a018831b4b..ed2bc78cb5 100644 --- a/src/Auxiliary/Aux_Check_Parameter.cpp +++ b/src/Auxiliary/Aux_Check_Parameter.cpp @@ -310,8 +310,13 @@ void Aux_Check_Parameter() { int NDerField = UserDerField_Num; # if ( MODEL == HYDRO ) - if ( OPT__OUTPUT_DIVVEL ) NDerField ++; - if ( OPT__OUTPUT_MACH ) NDerField ++; + if ( OPT__OUTPUT_DIVVEL ) NDerField ++; + if ( OPT__OUTPUT_MACH ) NDerField ++; +# ifdef SUPPORT_GRACKLE + if ( OPT__OUTPUT_GRACKLE_TEMP ) NDerField ++; + if ( OPT__OUTPUT_GRACKLE_MU ) NDerField ++; + if ( OPT__OUTPUT_GRACKLE_TCOOL ) NDerField ++; +# endif # endif if ( NDerField > DER_NOUT_MAX ) @@ -428,6 +433,9 @@ void Aux_Check_Parameter() # ifdef SRHD Flag |= OPT__FLAG_LRTZ_GRADIENT; # endif +# ifdef SUPPORT_GRACKLE + Flag |= OPT__FLAG_COOLING_LEN; +# endif # ifdef COSMIC_RAY Flag |= OPT__FLAG_CRAY; Flag |= OPT__FLAG_LOHNER_CRAY; @@ -1793,6 +1801,21 @@ void Aux_Check_Parameter() # error : ERROR : SUPPORT_GRACKLE must work with EOS_GAMMA/EOS_COSMIC_RAY !! # endif + if ( OPT__OUTPUT_GRACKLE_TEMP && ! GRACKLE_ACTIVATE ) + Aux_Error( ERROR_INFO, "OPT__OUTPUT_GRACKLE_TEMP requires the option GRACKLE_ACTIVATE !!\n"); + + if ( OPT__OUTPUT_GRACKLE_MU && ! GRACKLE_ACTIVATE ) + Aux_Error( ERROR_INFO, "OPT__OUTPUT_GRACKLE_MU requires the option GRACKLE_ACTIVATE !!\n"); + + if ( OPT__OUTPUT_GRACKLE_TCOOL && ! GRACKLE_ACTIVATE ) + Aux_Error( ERROR_INFO, "OPT__OUTPUT_GRACKLE_TCOOL requires the option GRACKLE_ACTIVATE !!\n"); + + if ( OPT__FLAG_COOLING_LEN && ! GRACKLE_ACTIVATE ) + Aux_Error( ERROR_INFO, "OPT__FLAG_COOLING_LEN requires the option GRACKLE_ACTIVATE !!\n"); + + if ( DT__GRACKLE_COOLING >= 0.0 && ! GRACKLE_ACTIVATE ) + Aux_Error( ERROR_INFO, "DT__GRACKLE_COOLING requires the option GRACKLE_ACTIVATE !!\n"); + // warning // ------------------------------ if ( MPI_Rank == 0 ) { @@ -1803,6 +1826,25 @@ void Aux_Check_Parameter() if ( GRACKLE_PRIMORDIAL > 0 ) Aux_Message( stderr, "WARNING : adiabatic index gamma is currently fixed to %13.7e for Grackle !!\n", GAMMA ); + if ( OPT__UNFREEZE_GRACKLE && ! OPT__FREEZE_FLUID ) + Aux_Message( stderr, "WARNING : OPT__UNFREEZE_GRACKLE is useless when OPT__FREEZE_FLUID is off !!\n" ); + + if ( OPT__UNFREEZE_GRACKLE && OPT__FREEZE_FLUID ) + Aux_Message( stderr, "REMINDER : OPT__UNFREEZE_GRACKLE will allow fluid variables to be updated by Grackle solver even with OPT__FREEZE_FLUID enabled !!\n" ); + + if ( OPT__OUTPUT_GRACKLE_TEMP && OPT__OUTPUT_TEMP ) + { + Aux_Message( stderr, "WARNING : temperature field calculated by Grackle is named \"GrackleTemp\"\n" ); + Aux_Message( stderr, " and it is different from the default temperature field \"Temp\" !!\n" ); + } + + if ( GRACKLE_PRIMORDIAL == GRACKLE_PRI_CHE_CLOUDY ) + { + Aux_Message( stderr, "WARNING : For GRACKLE_PRIMORDIAL == %d, GRACKLE_HYDROGEN_MFRAC = %13.7e is only used outside Grackle\n", + GRACKLE_PRI_CHE_CLOUDY, GRACKLE_HYDROGEN_MFRAC ); + Aux_Message( stderr, " and it can be different from HydrogenFractionByMass determined by and used inside Grackle !!\n" ); + } + } // if ( MPI_Rank == 0 ) #endif // #ifdef SUPPORT_GRACKLE diff --git a/src/Auxiliary/Aux_TakeNote.cpp b/src/Auxiliary/Aux_TakeNote.cpp index ef833f01fb..7e800e2966 100644 --- a/src/Auxiliary/Aux_TakeNote.cpp +++ b/src/Auxiliary/Aux_TakeNote.cpp @@ -980,6 +980,9 @@ void Aux_TakeNote() # ifdef CR_DIFFUSION fprintf( Note, "DT__CR_DIFFUSION % 14.7e\n", DT__CR_DIFFUSION ); # endif +# ifdef SUPPORT_GRACKLE + fprintf( Note, "DT__GRACKLE_COOLING % 14.7e\n", DT__GRACKLE_COOLING ); +# endif # ifdef COMOVING fprintf( Note, "DT__MAX_DELTA_A % 14.7e\n", DT__MAX_DELTA_A ); # endif @@ -1022,6 +1025,9 @@ void Aux_TakeNote() # ifdef SRHD fprintf( Note, "OPT__FLAG_LRTZ_GRADIENT % d\n", OPT__FLAG_LRTZ_GRADIENT ); # endif +# ifdef SUPPORT_GRACKLE + fprintf( Note, "OPT__FLAG_COOLING_LEN % d\n", OPT__FLAG_COOLING_LEN ); +# endif # endif # if ( MODEL == ELBDM ) fprintf( Note, "OPT__FLAG_ENGY_DENSITY % d\n", OPT__FLAG_ENGY_DENSITY ); @@ -1123,21 +1129,25 @@ void Aux_TakeNote() # ifdef SUPPORT_GRACKLE fprintf( Note, "Parameters of Grackle\n" ); fprintf( Note, "***********************************************************************************\n" ); - fprintf( Note, "GRACKLE_ACTIVATE % d\n", GRACKLE_ACTIVATE ); + fprintf( Note, "GRACKLE_ACTIVATE % d\n", GRACKLE_ACTIVATE ); if ( GRACKLE_ACTIVATE ) { - fprintf( Note, "GRACKLE_VERBOSE % d\n", GRACKLE_VERBOSE ); - fprintf( Note, "GRACKLE_COOLING % d\n", GRACKLE_COOLING ); - fprintf( Note, "GRACKLE_PRIMORDIAL % d\n", GRACKLE_PRIMORDIAL ); - fprintf( Note, "GRACKLE_METAL % d\n", GRACKLE_METAL ); - fprintf( Note, "GRACKLE_UV % d\n", GRACKLE_UV ); - fprintf( Note, "GRACKLE_CMB_FLOOR % d\n", GRACKLE_CMB_FLOOR ); - fprintf( Note, "GRACKLE_PE_HEATING % d\n", GRACKLE_PE_HEATING ); - fprintf( Note, "GRACKLE_PE_HEATING_RATE % 14.7e\n", GRACKLE_PE_HEATING_RATE ); - fprintf( Note, "GRACKLE_CLOUDY_TABLE %s\n", GRACKLE_CLOUDY_TABLE ); - fprintf( Note, "GRACKLE_THREE_BODY_RATE % d\n", GRACKLE_THREE_BODY_RATE ); - fprintf( Note, "GRACKLE_CIE_COOLING % d\n", GRACKLE_CIE_COOLING ); - fprintf( Note, "GRACKLE_H2_OPA_APPROX % d\n", GRACKLE_H2_OPA_APPROX ); - fprintf( Note, "CHE_GPU_NPGROUP % d\n", CHE_GPU_NPGROUP ); } + fprintf( Note, "GRACKLE_VERBOSE % d\n", GRACKLE_VERBOSE ); + fprintf( Note, "GRACKLE_COOLING % d\n", GRACKLE_COOLING ); + fprintf( Note, "GRACKLE_PRIMORDIAL % d\n", GRACKLE_PRIMORDIAL ); + fprintf( Note, "GRACKLE_METAL % d\n", GRACKLE_METAL ); + fprintf( Note, "GRACKLE_UV % d\n", GRACKLE_UV ); + fprintf( Note, "GRACKLE_CMB_FLOOR % d\n", GRACKLE_CMB_FLOOR ); + fprintf( Note, "GRACKLE_PE_HEATING % d\n", GRACKLE_PE_HEATING ); + fprintf( Note, "GRACKLE_PE_HEATING_RATE % 14.7e\n", GRACKLE_PE_HEATING_RATE ); + fprintf( Note, "GRACKLE_CLOUDY_TABLE %s\n", GRACKLE_CLOUDY_TABLE ); + fprintf( Note, "GRACKLE_THREE_BODY_RATE % d\n", GRACKLE_THREE_BODY_RATE ); + fprintf( Note, "GRACKLE_CIE_COOLING % d\n", GRACKLE_CIE_COOLING ); + fprintf( Note, "GRACKLE_H2_OPA_APPROX % d\n", GRACKLE_H2_OPA_APPROX ); + fprintf( Note, "GRACKLE_USE_V_HEATING_RATE % d\n", GRACKLE_USE_V_HEATING_RATE ); + fprintf( Note, "GRACKLE_USE_S_HEATING_RATE % d\n", GRACKLE_USE_S_HEATING_RATE ); + fprintf( Note, "GRACKLE_HYDROGEN_MFRAC % 14.7e\n", GRACKLE_HYDROGEN_MFRAC ); + fprintf( Note, "OPT__UNFREEZE_GRACKLE % d\n", OPT__UNFREEZE_GRACKLE ); + fprintf( Note, "CHE_GPU_NPGROUP % d\n", CHE_GPU_NPGROUP ); } fprintf( Note, "***********************************************************************************\n" ); fprintf( Note, "\n\n" ); # endif // #ifdef SUPPORT_GRACKLE @@ -1619,7 +1629,6 @@ void Aux_TakeNote() fprintf( Note, "OPT__OUTPUT_CS % d\n", OPT__OUTPUT_CS ); fprintf( Note, "OPT__OUTPUT_DIVVEL % d\n", OPT__OUTPUT_DIVVEL ); fprintf( Note, "OPT__OUTPUT_MACH % d\n", OPT__OUTPUT_MACH ); -# endif # ifdef MHD fprintf( Note, "OPT__OUTPUT_DIVMAG % d\n", OPT__OUTPUT_DIVMAG ); # endif @@ -1629,6 +1638,12 @@ void Aux_TakeNote() fprintf( Note, "OPT__OUTPUT_LORENTZ % d\n", OPT__OUTPUT_LORENTZ ); fprintf( Note, "OPT__OUTPUT_ENTHALPY % d\n", OPT__OUTPUT_ENTHALPY ); # endif +# ifdef SUPPORT_GRACKLE + fprintf( Note, "OPT__OUTPUT_GRACKLE_TEMP % d\n", OPT__OUTPUT_GRACKLE_TEMP ); + fprintf( Note, "OPT__OUTPUT_GRACKLE_MU % d\n", OPT__OUTPUT_GRACKLE_MU ); + fprintf( Note, "OPT__OUTPUT_GRACKLE_TCOOL % d\n", OPT__OUTPUT_GRACKLE_TCOOL ); +# endif +# endif // user-defined derived fields if ( OPT__OUTPUT_USER_FIELD ) { @@ -1832,6 +1847,18 @@ void Aux_TakeNote() } # endif +# ifdef SUPPORT_GRACKLE + if ( OPT__FLAG_COOLING_LEN ) + { + fprintf( Note, "Flag Criterion (Cooling Length over Cell Size in HYDRO+GRACKLE)\n" ); + fprintf( Note, "***********************************************************************************\n" ); + fprintf( Note, " Level l_cool / dh\n" ); + for (int lv=0; lv declared in Init_MemAllocate_Grackle.cpp +extern int Che_NField; +extern int CheIdx_sEint; + + + + +//------------------------------------------------------------------------------------------------------- +// Function : Grackle_Calculate +// Description : Calculate the fields by calling Grackle's API functions +// +// Note : 1. Ref: https://grackle.readthedocs.io/en/latest/Interaction.html#calling-the-available-functions +// 2. Currently, the fields have to be calculated in units of patch group +// 3. The output field data type is "real" instead of "real_che" +// --> to be consistent with the usage of the output +// 4. This function can be invoked by +// e.g., Output_DumpData_*(), Flag_Real() and Grackle_GetTimeStep_CoolingTime() +// +// Parameter : Out : Output array, array size = NFieldOut*(NPG*PS2*PS2*PS2) +// TFields : Target fields to be calculated: +// _GRACKLE_TEMP, _GRACKLE_MU, _GRACKLE_TCOOL +// --> defined in include/Typedef.h +// lv : Target refinement level +// NPG : Number of patch groups calculated at a time +// PID0_List : List recording the patch indices with LocalID==0 to be calculated +// +// Return : Out[] +//------------------------------------------------------------------------------------------------------- +void Grackle_Calculate( real Out[], const GrackleFieldBIdx_t TFields, + const int lv, const int NPG, const int *PID0_List ) +{ + +// 0. nothing to do if there is no target patch group + if ( NPG == 0 ) return; + + const int Size1v = NPG * CUBE(PS2); // size of one field + + +// 1. check the number of output fields + int NFieldOut = 0; + const int IdxOut_Undefined = -1; + const int IdxOut_temp = ( TFields & _GRACKLE_TEMP ) ? NFieldOut++ : IdxOut_Undefined; + const int IdxOut_mu = ( TFields & _GRACKLE_MU ) ? NFieldOut++ : IdxOut_Undefined; + const int IdxOut_tcool = ( TFields & _GRACKLE_TCOOL ) ? NFieldOut++ : IdxOut_Undefined; + + if ( NFieldOut == 0 ) return; + + +// 2. allocate and prepare the input array for the Grackle solver + real_che *gr_fields_input = new real_che[ (long)Che_NField*(long)Size1v ]; + + Grackle_Prepare( lv, gr_fields_input, NPG, PID0_List ); + + +// 3. allocate the output array for the Grackle solver + real_che *gr_fields_temperature = NULL; + real_che *gr_fields_gamma = NULL; + real_che *gr_fields_cooling_time = NULL; + + if ( IdxOut_temp != IdxOut_Undefined || IdxOut_mu != IdxOut_Undefined ) + gr_fields_temperature = new real_che[Size1v]; + + if ( IdxOut_mu != IdxOut_Undefined ) + gr_fields_gamma = new real_che[Size1v]; + + if ( IdxOut_tcool != IdxOut_Undefined ) + gr_fields_cooling_time = new real_che[Size1v]; + + typedef real (*vla_out)[Size1v]; + vla_out Out1D = ( vla_out )Out; + + +// 4. set grid_dimension, grid_start, and grid_end + const int OptFac = 16; // optimization factor + if ( SQR(PS2)%OptFac != 0 ) Aux_Error( ERROR_INFO, "SQR(PS2) %% OptFac != 0 !!\n" ); + + Che_FieldData->grid_dimension[0] = PS2*OptFac; + Che_FieldData->grid_dimension[1] = 1; + Che_FieldData->grid_dimension[2] = SQR(PS2)*NPG/OptFac; + + for (int d=0; d<3; d++) + { + Che_FieldData->grid_start[d] = 0; + Che_FieldData->grid_end [d] = Che_FieldData->grid_dimension[d] - 1; + } + + +// 5. invoke Grackle +// --> note that we use the OpenMP implementation in Grackle directly, which applies the parallelization to the first two +// dimensiones of the input grid +// --> this approach is found to be much more efficient than parallelizing different patches or patch groups here + if ( IdxOut_temp != IdxOut_Undefined || IdxOut_mu != IdxOut_Undefined ) + { + if ( calculate_temperature( &Che_Units, Che_FieldData, gr_fields_temperature ) == 0 ) + Aux_Error( ERROR_INFO, "Grackle calculate_temperature() failed !!\n" ); + } + + if ( IdxOut_mu != IdxOut_Undefined ) + { + if ( calculate_gamma( &Che_Units, Che_FieldData, gr_fields_gamma ) == 0 ) + Aux_Error( ERROR_INFO, "Grackle calculate_gamma() failed !!\n" ); + } + + if ( IdxOut_tcool != IdxOut_Undefined ) + { + if ( calculate_cooling_time( &Che_Units, Che_FieldData, gr_fields_cooling_time ) == 0 ) + Aux_Error( ERROR_INFO, "Grackle calculate_cooling_time() failed !!\n" ); + } + + +// 6. fill in the output array + const real_che *gr_fields_sEint = gr_fields_input + CheIdx_sEint*Size1v; + const real_che temperature_units = (real_che) get_temperature_units( &Che_Units ); + + for (int i=0; i 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. Time-step is set to restrict the cooling process +// --> dt = DT__GRACKLE_COOLING * Min(CoolingTime) +// 3. This is usually for the purpose of outputting data for analysis; +// Inside the Grackle solver, subcycling time-steps will be adopted for the correctness of evolution +// +// Parameter : lv : Target refinement level +// +// Return : dt +//------------------------------------------------------------------------------------------------------- +double Grackle_GetTimeStep_CoolingTime( const int lv ) +{ + +// get the minimum cooling time + const real MinCoolingTime = GetMinCoolingTime( lv ); + +// get the time-step + const double dt = DT__GRACKLE_COOLING * MinCoolingTime; + + return dt; + +} // FUNCTION : Grackle_GetTimeStep_CoolingTime + + + +//------------------------------------------------------------------------------------------------------- +// Function : GetMinCoolingTime +// Description : Get the minimum Grackle cooling time at the target level among all MPI ranks +// +// Note : 1. Invoked by Grackle_GetTimeStep_CoolingTime() +// +// Parameter : lv : Target refinement level +// +// Return : MinCoolingTime_AllRank +//------------------------------------------------------------------------------------------------------- +real GetMinCoolingTime( const int lv ) +{ + + real MinCoolingTime = INFINITY; + bool AnyCell = false; + +# pragma omp parallel + { + + real *Grackle_TCool = new real [ CUBE(PS2) ]; // array storing ONE patch group of grackle cooling time + +// get the minimum Grckle cooling time in this rank +# pragma omp for reduction( min:MinCoolingTime ) reduction( ||:AnyCell ) schedule( runtime ) + for (int PID0=0; PID0NPatchComma[lv][1]; PID0+=8) + { + +// calculate the Grackle cooling time + Grackle_Calculate( Grackle_TCool, _GRACKLE_TCOOL, lv, 1, &PID0 ); + + for (int LocalID=0; LocalID<8; LocalID++) + { + const int PID = PID0 + LocalID; + +// if OPT__FIXUP_RESTRICT is enabled, skip all non-leaf patches +// becausethey are later overwritten by the refined patches +// note that this leads to the timestep being "inf" when a level is completely refined + if ( OPT__FIXUP_RESTRICT && amr->patch[0][lv][PID]->son != -1 ) continue; + + AnyCell = true; + +// calculate the minimum + for (int k=0; kNPatchComma[lv][1]; PID0+=8) + + delete [] Grackle_TCool; + + } // OpenMP parallel region + + +// get the minimum Grackle cooling time in all ranks + real MinCoolingTime_AllRank; + bool AnyCell_AllRank; + + MPI_Allreduce( &MinCoolingTime, &MinCoolingTime_AllRank, 1, MPI_GAMER_REAL, MPI_MIN, MPI_COMM_WORLD ); + MPI_Reduce( &AnyCell, &AnyCell_AllRank, 1, MPI_C_BOOL, MPI_LOR, 0, MPI_COMM_WORLD ); + + +// check + if ( MinCoolingTime_AllRank == INFINITY && AnyCell_AllRank && MPI_Rank == 0 ) + Aux_Error( ERROR_INFO, "MinCoolingTime == INFINITY at lv %d !!\n", lv ); + + + return MinCoolingTime_AllRank; + +} // FUNCTION : GetMinCoolingTime + + + +#endif // #ifdef SUPPORT_GRACKLE diff --git a/src/Grackle/Grackle_Init.cpp b/src/Grackle/Grackle_Init.cpp index a6843c9fc6..4a4e8b8ec3 100644 --- a/src/Grackle/Grackle_Init.cpp +++ b/src/Grackle/Grackle_Init.cpp @@ -99,6 +99,13 @@ void Grackle_Init() grackle_data->three_body_rate = GRACKLE_THREE_BODY_RATE; grackle_data->cie_cooling = GRACKLE_CIE_COOLING; grackle_data->h2_optical_depth_approximation = GRACKLE_H2_OPA_APPROX; + grackle_data->use_volumetric_heating_rate = GRACKLE_USE_V_HEATING_RATE; + grackle_data->use_specific_heating_rate = GRACKLE_USE_S_HEATING_RATE; +// hydrogen mass fraction is only set when it is in non-equilibrium mode +// because the tables for tabulated mode were created assuming hydrogen mass fraction of about 0.716 +// see https://grackle.readthedocs.io/en/latest/Parameters.html#c.HydrogenFractionByMass + if ( GRACKLE_PRIMORDIAL != GRACKLE_PRI_CHE_CLOUDY ) + grackle_data->HydrogenFractionByMass = GRACKLE_HYDROGEN_MFRAC; # ifdef OPENMP // currently we adopt the OpenMP implementation in Grackle directly, which applies the parallelization to diff --git a/src/Grackle/Grackle_Init_FieldData.cpp b/src/Grackle/Grackle_Init_FieldData.cpp index d8b0d8eafc..7c11268b11 100644 --- a/src/Grackle/Grackle_Init_FieldData.cpp +++ b/src/Grackle/Grackle_Init_FieldData.cpp @@ -64,13 +64,13 @@ void Grackle_Init_FieldData() Che_FieldData->DII_density = NULL; Che_FieldData->HDI_density = NULL; Che_FieldData->metal_density = NULL; + Che_FieldData->volumetric_heating_rate = NULL; + Che_FieldData->specific_heating_rate = NULL; // fields not supported yet Che_FieldData->x_velocity = NULL; Che_FieldData->y_velocity = NULL; Che_FieldData->z_velocity = NULL; - Che_FieldData->volumetric_heating_rate = NULL; - Che_FieldData->specific_heating_rate = NULL; Che_FieldData->RT_HI_ionization_rate = NULL; Che_FieldData->RT_HeI_ionization_rate = NULL; Che_FieldData->RT_HeII_ionization_rate = NULL; diff --git a/src/Grackle/Grackle_Prepare.cpp b/src/Grackle/Grackle_Prepare.cpp index dd71c474ff..9f7a7a241f 100644 --- a/src/Grackle/Grackle_Prepare.cpp +++ b/src/Grackle/Grackle_Prepare.cpp @@ -22,8 +22,19 @@ extern int CheIdx_DI; extern int CheIdx_DII; extern int CheIdx_HDI; extern int CheIdx_Metal; +extern int CheIdx_vHeatingRate; +extern int CheIdx_sHeatingRate; +// declare as static so that other functions cannot invoke it directly and must use the function pointer +static real_che Grackle_vHeatingRate_User_Template( const double x, const double y, const double z, const double Time, const double n_H ); +static real_che Grackle_sHeatingRate_User_Template( const double x, const double y, const double z, const double Time ); + + +// these function pointers must be set by a test problem initializer +real_che (*Grackle_vHeatingRate_User_Ptr)( const double x, const double y, const double z, const double Time, const double n_H ) = NULL; +real_che (*Grackle_sHeatingRate_User_Ptr)( const double x, const double y, const double z, const double Time ) = NULL; + //------------------------------------------------------------------------------------------------------- @@ -89,6 +100,20 @@ void Grackle_Prepare( const int lv, real_che h_Che_Array[], const int NPG, const if ( Idx_Metal == Idx_Undefined || CheIdx_Metal == Idx_Undefined ) Aux_Error( ERROR_INFO, "[Che]Idx_Metal is undefined for \"GRACKLE_METAL\" !!\n" ); } + + if ( GRACKLE_USE_V_HEATING_RATE ) { + if ( Grackle_vHeatingRate_User_Ptr == NULL ) + Aux_Error( ERROR_INFO, "Grackle_vHeatingRate_User_Ptr == NULL !!\n" ); + if ( CheIdx_vHeatingRate == Idx_Undefined ) + Aux_Error( ERROR_INFO, "CheIdx_vHeatingRate is undefined for \"GRACKLE_USE_V_HEATING_RATE\" !!\n" ); + } + + if ( GRACKLE_USE_S_HEATING_RATE ) { + if ( Grackle_sHeatingRate_User_Ptr == NULL ) + Aux_Error( ERROR_INFO, "Grackle_sHeatingRate_User_Ptr == NULL !!\n" ); + if ( CheIdx_sHeatingRate == Idx_Undefined ) + Aux_Error( ERROR_INFO, "CheIdx_sHeatingRate is undefined for \"GRACKLE_USE_S_HEATING_RATE\" !!\n" ); + } # endif // #ifdef GAMER_DEBUG @@ -100,23 +125,26 @@ void Grackle_Prepare( const int lv, real_che h_Che_Array[], const int NPG, const # else const bool CheckMinEint_Yes = true; # endif - - real_che *Ptr_Dens0 = h_Che_Array + CheIdx_Dens *Size1v; - real_che *Ptr_sEint0 = h_Che_Array + CheIdx_sEint*Size1v; - real_che *Ptr_Ent0 = h_Che_Array + CheIdx_Ent *Size1v; - real_che *Ptr_e0 = h_Che_Array + CheIdx_e *Size1v; - real_che *Ptr_HI0 = h_Che_Array + CheIdx_HI *Size1v; - real_che *Ptr_HII0 = h_Che_Array + CheIdx_HII *Size1v; - real_che *Ptr_HeI0 = h_Che_Array + CheIdx_HeI *Size1v; - real_che *Ptr_HeII0 = h_Che_Array + CheIdx_HeII *Size1v; - real_che *Ptr_HeIII0 = h_Che_Array + CheIdx_HeIII*Size1v; - real_che *Ptr_HM0 = h_Che_Array + CheIdx_HM *Size1v; - real_che *Ptr_H2I0 = h_Che_Array + CheIdx_H2I *Size1v; - real_che *Ptr_H2II0 = h_Che_Array + CheIdx_H2II *Size1v; - real_che *Ptr_DI0 = h_Che_Array + CheIdx_DI *Size1v; - real_che *Ptr_DII0 = h_Che_Array + CheIdx_DII *Size1v; - real_che *Ptr_HDI0 = h_Che_Array + CheIdx_HDI *Size1v; - real_che *Ptr_Metal0 = h_Che_Array + CheIdx_Metal*Size1v; + const double dh = amr->dh[lv]; + + real_che *Ptr_Dens0 = h_Che_Array + CheIdx_Dens *Size1v; + real_che *Ptr_sEint0 = h_Che_Array + CheIdx_sEint *Size1v; + real_che *Ptr_Ent0 = h_Che_Array + CheIdx_Ent *Size1v; + real_che *Ptr_e0 = h_Che_Array + CheIdx_e *Size1v; + real_che *Ptr_HI0 = h_Che_Array + CheIdx_HI *Size1v; + real_che *Ptr_HII0 = h_Che_Array + CheIdx_HII *Size1v; + real_che *Ptr_HeI0 = h_Che_Array + CheIdx_HeI *Size1v; + real_che *Ptr_HeII0 = h_Che_Array + CheIdx_HeII *Size1v; + real_che *Ptr_HeIII0 = h_Che_Array + CheIdx_HeIII *Size1v; + real_che *Ptr_HM0 = h_Che_Array + CheIdx_HM *Size1v; + real_che *Ptr_H2I0 = h_Che_Array + CheIdx_H2I *Size1v; + real_che *Ptr_H2II0 = h_Che_Array + CheIdx_H2II *Size1v; + real_che *Ptr_DI0 = h_Che_Array + CheIdx_DI *Size1v; + real_che *Ptr_DII0 = h_Che_Array + CheIdx_DII *Size1v; + real_che *Ptr_HDI0 = h_Che_Array + CheIdx_HDI *Size1v; + real_che *Ptr_Metal0 = h_Che_Array + CheIdx_Metal *Size1v; + real_che *Ptr_vHeatingRate0 = h_Che_Array + CheIdx_vHeatingRate*Size1v; + real_che *Ptr_sHeatingRate0 = h_Che_Array + CheIdx_sHeatingRate*Size1v; # pragma omp parallel @@ -135,30 +163,33 @@ void Grackle_Prepare( const int lv, real_che h_Che_Array[], const int NPG, const real_che *Ptr_Dens=NULL, *Ptr_sEint=NULL, *Ptr_Ent=NULL, *Ptr_e=NULL, *Ptr_HI=NULL, *Ptr_HII=NULL; real_che *Ptr_HeI=NULL, *Ptr_HeII=NULL, *Ptr_HeIII=NULL, *Ptr_HM=NULL, *Ptr_H2I=NULL, *Ptr_H2II=NULL; real_che *Ptr_DI=NULL, *Ptr_DII=NULL, *Ptr_HDI=NULL, *Ptr_Metal=NULL; + real_che *Ptr_vHeatingRate=NULL, *Ptr_sHeatingRate=NULL; # pragma omp for schedule( static ) for (int TID=0; TIDpatch[ amr->FluSg[lv] ][lv][PID]->fluid; + const double x0 = amr->patch[0][lv][PID]->EdgeL[0] + 0.5*dh; + const double y0 = amr->patch[0][lv][PID]->EdgeL[1] + 0.5*dh; + const double z0 = amr->patch[0][lv][PID]->EdgeL[2] + 0.5*dh; + for (int k=0; kmetal_density = Ptr_Metal0; + if ( GRACKLE_USE_V_HEATING_RATE ) + Che_FieldData->volumetric_heating_rate = Ptr_vHeatingRate0; + + if ( GRACKLE_USE_S_HEATING_RATE ) + Che_FieldData->specific_heating_rate = Ptr_sHeatingRate0; + } // FUNCTION : Grackle_Prepare +//------------------------------------------------------------------------------------------------------- +// Function : Grackle_vHeatingRate_User_Template +// Description : Function template to set Grackle's volumetric heating rate +// +// Note : 1. Invoked by Grackle_Prepare() using the function pointer +// "Grackle_vHeatingRate_User_Ptr", which must be set by a test problem initializer +// 2. This function will be invoked by multiple OpenMP threads when OPENMP is enabled +// --> Please ensure that everything here is thread-safe +// 3. Returned rate should be in the unit of erg s^-1 cm^-3 +// +// Parameter : x/y/z : Target physical coordinates +// Time : Target physical time +// n_H : Hydrogen number density, should be in the unit of cm^-3 +// +// Return : volumetric_heating_rate +//------------------------------------------------------------------------------------------------------- +real_che Grackle_vHeatingRate_User_Template( const double x, const double y, const double z, const double Time, const double n_H ) +{ + + const double Center[3] = { amr->BoxCenter[0], amr->BoxCenter[1], amr->BoxCenter[2] }; + const double radius_to_center = sqrt( SQR(x-Center[0]) + SQR(y-Center[1]) ); + const double ScaleLength = 0.125*amr->BoxSize[0]; + const double R_core = 0.125*amr->BoxSize[0]; + const double volumetric_heating_rate_0 = n_H * GRACKLE_PE_HEATING_RATE; // GRACKLE_PE_HEATING_RATE is in units of erg cm^-3 s^-1 n_H^-1 + const real_che volumetric_heating_rate = volumetric_heating_rate_0 * exp( (R_core - MAX(radius_to_center, R_core)) / ScaleLength ); + + return volumetric_heating_rate; + +} // FUNCTION : Grackle_vHeatingRate_User_Template + + + +//------------------------------------------------------------------------------------------------------- +// Function : Grackle_sHeatingRate_User_Template +// Description : Function template to set Grackle's specific heating rate +// +// Note : 1. Invoked by Grackle_Prepare() using the function pointer +// "Grackle_sHeatingRate_User_Ptr", which must be set by a test problem initializer +// 2. This function will be invoked by multiple OpenMP threads when OPENMP is enabled +// --> Please ensure that everything here is thread-safe +// 3. Returned rate should be in the units of erg s^-1 g^-1 +// +// Parameter : x/y/z : Target physical coordinates +// Time : Target physical time +// +// Return : specific_heating_rate +//------------------------------------------------------------------------------------------------------- +real_che Grackle_sHeatingRate_User_Template( const double x, const double y, const double z, const double Time ) +{ + + const double Center[3] = { amr->BoxCenter[0], amr->BoxCenter[1], amr->BoxCenter[2] }; + const double radius_to_center = sqrt( SQR(x-Center[0]) + SQR(y-Center[1]) ); + const double ScaleLength = 0.125*amr->BoxSize[0]; + const double R_core = 0.125*amr->BoxSize[0]; + const double specific_heating_rate_0 = GRACKLE_PE_HEATING_RATE * GRACKLE_HYDROGEN_MFRAC / Const_mH; // GRACKLE_PE_HEATING_RATE is in units of erg cm^-3 s^-1 n_H^-1 + const real_che specific_heating_rate = specific_heating_rate_0 * exp( (R_core - MAX(radius_to_center, R_core)) / ScaleLength ); + + return specific_heating_rate; + +} // FUNCTION : Grackle_sHeatingRate_User_Template + + + #endif // #ifdef SUPPORT_GRACKLE diff --git a/src/Grackle/Init_MemAllocate_Grackle.cpp b/src/Grackle/Init_MemAllocate_Grackle.cpp index 36a64b2164..1045898631 100644 --- a/src/Grackle/Init_MemAllocate_Grackle.cpp +++ b/src/Grackle/Init_MemAllocate_Grackle.cpp @@ -6,23 +6,25 @@ // global variables for accessing h_Che_Array[] // --> also used by Grackle_Prepare.cpp and Grackle_Close.cpp // --> they are not declared in "Global.h" simply because they are only used by a few Grackle routines -int Che_NField = NULL_INT; -int CheIdx_Dens = Idx_Undefined; -int CheIdx_sEint = Idx_Undefined; -int CheIdx_Ent = Idx_Undefined; -int CheIdx_e = Idx_Undefined; -int CheIdx_HI = Idx_Undefined; -int CheIdx_HII = Idx_Undefined; -int CheIdx_HeI = Idx_Undefined; -int CheIdx_HeII = Idx_Undefined; -int CheIdx_HeIII = Idx_Undefined; -int CheIdx_HM = Idx_Undefined; -int CheIdx_H2I = Idx_Undefined; -int CheIdx_H2II = Idx_Undefined; -int CheIdx_DI = Idx_Undefined; -int CheIdx_DII = Idx_Undefined; -int CheIdx_HDI = Idx_Undefined; -int CheIdx_Metal = Idx_Undefined; +int Che_NField = NULL_INT; +int CheIdx_Dens = Idx_Undefined; +int CheIdx_sEint = Idx_Undefined; +int CheIdx_Ent = Idx_Undefined; +int CheIdx_e = Idx_Undefined; +int CheIdx_HI = Idx_Undefined; +int CheIdx_HII = Idx_Undefined; +int CheIdx_HeI = Idx_Undefined; +int CheIdx_HeII = Idx_Undefined; +int CheIdx_HeIII = Idx_Undefined; +int CheIdx_HM = Idx_Undefined; +int CheIdx_H2I = Idx_Undefined; +int CheIdx_H2II = Idx_Undefined; +int CheIdx_DI = Idx_Undefined; +int CheIdx_DII = Idx_Undefined; +int CheIdx_HDI = Idx_Undefined; +int CheIdx_Metal = Idx_Undefined; +int CheIdx_vHeatingRate = Idx_Undefined; +int CheIdx_sHeatingRate = Idx_Undefined; @@ -76,6 +78,12 @@ void Init_MemAllocate_Grackle( const int Che_NPG ) if ( GRACKLE_METAL ) CheIdx_Metal = Che_NField ++; + if ( GRACKLE_USE_V_HEATING_RATE ) + CheIdx_vHeatingRate = Che_NField ++; + + if ( GRACKLE_USE_S_HEATING_RATE ) + CheIdx_sHeatingRate = Che_NField ++; + // allocate the input/output array for the Grackle solver for (int t=0; t<2; t++) diff --git a/src/Init/Init_ByRestart_HDF5.cpp b/src/Init/Init_ByRestart_HDF5.cpp index e08c16a93f..dfbdb4f58d 100644 --- a/src/Init/Init_ByRestart_HDF5.cpp +++ b/src/Init/Init_ByRestart_HDF5.cpp @@ -2015,6 +2015,9 @@ void Check_InputPara( const char *FileName, const int FormatVersion ) # ifdef CR_DIFFUSION LoadField( "Dt__CR_Diffusion", &RS.Dt__CR_Diffusion, SID, TID, NonFatal, &RT.Dt__CR_Diffusion, 1, NonFatal ); # endif +# ifdef SUPPORT_GRACKLE + LoadField( "Dt__GrackleCooling", &RS.Dt__GrackleCooling, SID, TID, NonFatal, &RT.Dt__GrackleCooling, 1, NonFatal ); +# endif # ifdef COMOVING LoadField( "Dt__MaxDeltaA", &RS.Dt__MaxDeltaA, SID, TID, NonFatal, &RT.Dt__MaxDeltaA, 1, NonFatal ); # endif @@ -2055,6 +2058,9 @@ void Check_InputPara( const char *FileName, const int FormatVersion ) # ifdef SRHD LoadField( "Opt__Flag_LrtzGradient", &RS.Opt__Flag_LrtzGradient, SID, TID, NonFatal, &RT.Opt__Flag_LrtzGradient, 1, NonFatal ); # endif +# ifdef SUPPORT_GRACKLE + LoadField( "Opt__Flag_CoolingLen", &RS.Opt__Flag_CoolingLen, SID, TID, NonFatal, &RT.Opt__Flag_CoolingLen, 1, NonFatal ); +# endif # ifdef COSMIC_RAY LoadField( "Opt__Flag_CRay", &RS.Opt__Flag_CRay, SID, TID, NonFatal, &RT.Opt__Flag_CRay, 1, NonFatal ); # endif @@ -2225,6 +2231,10 @@ void Check_InputPara( const char *FileName, const int FormatVersion ) LoadField( "Grackle_ThreeBodyRate", &RS.Grackle_ThreeBodyRate, SID, TID, NonFatal, &RT.Grackle_ThreeBodyRate, 1, NonFatal ); LoadField( "Grackle_CIE_Cooling", &RS.Grackle_CIE_Cooling, SID, TID, NonFatal, &RT.Grackle_CIE_Cooling, 1, NonFatal ); LoadField( "Grackle_H2_OpaApprox", &RS.Grackle_H2_OpaApprox, SID, TID, NonFatal, &RT.Grackle_H2_OpaApprox, 1, NonFatal ); + LoadField( "Grackle_UseVHeatingRate", &RS.Grackle_UseVHeatingRate, SID, TID, NonFatal, &RT.Grackle_UseVHeatingRate, 1, NonFatal ); + LoadField( "Grackle_UseSHeatingRate", &RS.Grackle_UseSHeatingRate, SID, TID, NonFatal, &RT.Grackle_UseSHeatingRate, 1, NonFatal ); + LoadField( "Grackle_HydrogenMFrac", &RS.Grackle_HydrogenMFrac, SID, TID, NonFatal, &RT.Grackle_HydrogenMFrac, 1, NonFatal ); + LoadField( "Opt__UnfreezeGrackle", &RS.Opt__UnfreezeGrackle, SID, TID, NonFatal, &RT.Opt__UnfreezeGrackle, 1, NonFatal ); LoadField( "Che_GPU_NPGroup", &RS.Che_GPU_NPGroup, SID, TID, NonFatal, &RT.Che_GPU_NPGroup, 1, NonFatal ); # endif @@ -2355,6 +2365,11 @@ void Check_InputPara( const char *FileName, const int FormatVersion ) LoadField( "Opt__Output_3Velocity", &RS.Opt__Output_3Velocity, SID, TID, NonFatal, &RT.Opt__Output_3Velocity, 1, NonFatal ); LoadField( "Opt__Output_Enthalpy", &RS.Opt__Output_Enthalpy, SID, TID, NonFatal, &RT.Opt__Output_Enthalpy, 1, NonFatal ); # endif +# ifdef SUPPORT_GRACKLE + LoadField( "Opt__Output_GrackleTemp", &RS.Opt__Output_GrackleTemp, SID, TID, NonFatal, &RT.Opt__Output_GrackleTemp, 1, NonFatal ); + LoadField( "Opt__Output_GrackleMu", &RS.Opt__Output_GrackleMu, SID, TID, NonFatal, &RT.Opt__Output_GrackleMu, 1, NonFatal ); + LoadField( "Opt__Output_GrackleTCool", &RS.Opt__Output_GrackleTCool, SID, TID, NonFatal, &RT.Opt__Output_GrackleTCool, 1, NonFatal ); +# endif # endif // #if ( MODEL == HYDRO ) LoadField( "Opt__Output_UserField", &RS.Opt__Output_UserField, SID, TID, NonFatal, &RT.Opt__Output_UserField, 1, NonFatal ); # ifdef PARTICLE @@ -2474,6 +2489,9 @@ void Check_InputPara( const char *FileName, const int FormatVersion ) # ifdef SRHD RS.FlagTable_LrtzGradient[lv] = -1.0; # endif +# ifdef SRHD + RS.FlagTable_CoolingLen [lv] = -1.0; +# endif # elif ( MODEL == ELBDM ) for (int t=0; t<2; t++) { @@ -2573,6 +2591,11 @@ void Check_InputPara( const char *FileName, const int FormatVersion ) LoadField( "FlagTable_LrtzGradient", RS.FlagTable_LrtzGradient, SID, TID, NonFatal, RT.FlagTable_LrtzGradient, N1, NonFatal ); # endif +# ifdef SUPPORT_GRACKLE + if ( OPT__FLAG_COOLING_LEN ) + LoadField( "FlagTable_CoolingLen", RS.FlagTable_CoolingLen, SID, TID, NonFatal, RT.FlagTable_CoolingLen, N1, NonFatal ); +# endif + # elif ( MODEL == ELBDM ) if ( OPT__FLAG_ENGY_DENSITY ) { LoadField( "FlagTable_EngyDensity", RS.FlagTable_EngyDensity, SID, TID, NonFatal, NullPtr, -1, NonFatal ); diff --git a/src/Init/Init_Load_FlagCriteria.cpp b/src/Init/Init_Load_FlagCriteria.cpp index a856bf55bc..826ca709ef 100644 --- a/src/Init/Init_Load_FlagCriteria.cpp +++ b/src/Init/Init_Load_FlagCriteria.cpp @@ -35,6 +35,12 @@ void Init_Load_FlagCriteria() const bool OPT__FLAG_LRTZ_GRADIENT = false; double *FlagTable_LrtzGradient = NULL; # endif + +# ifndef SUPPORT_GRACKLE + const bool OPT__FLAG_COOLING_LEN = false; + double *FlagTable_CoolingLen = NULL; +# endif + # ifndef COSMIC_RAY const bool OPT__FLAG_LOHNER_CRAY = false; const bool OPT__FLAG_CRAY = false; @@ -74,32 +80,35 @@ void Init_Load_FlagCriteria() # error : unsupported MODEL !! # endif - const int NFlagMode = 18; + const int NFlagMode = 19; const bool Flag[NFlagMode] = { OPT__FLAG_RHO, OPT__FLAG_RHO_GRADIENT, OPT__FLAG_PRES_GRADIENT, OPT__FLAG_ENGY_DENSITY, OPT__FLAG_LOHNER, OPT__FLAG_USER, (bool)OPT__FLAG_NPAR_PATCH, OPT__FLAG_NPAR_CELL, OPT__FLAG_PAR_MASS_CELL, OPT__FLAG_VORTICITY, OPT__FLAG_JEANS, OPT__FLAG_CURRENT, OPT__FLAG_CRAY, OPT__FLAG_LRTZ_GRADIENT, OPT__FLAG_INTERFERENCE, - OPT__FLAG_SPECTRAL, OPT__FLAG_ANGULAR, OPT__FLAG_RADIAL }; + OPT__FLAG_SPECTRAL, OPT__FLAG_ANGULAR, OPT__FLAG_RADIAL, + OPT__FLAG_COOLING_LEN }; const char ModeName[][100] = { "OPT__FLAG_RHO", "OPT__FLAG_RHO_GRADIENT", "OPT__FLAG_PRES_GRADIENT", "OPT__FLAG_ENGY_DENSITY", "OPT__FLAG_LOHNER", "OPT__FLAG_USER", "OPT__FLAG_NPAR_PATCH", "OPT__FLAG_NPAR_CELL", "OPT__FLAG_PAR_MASS_CELL", "OPT__FLAG_VORTICITY", "OPT__FLAG_JEANS", "OPT__FLAG_CURRENT", "OPT__FLAG_CRAY", "OPT__FLAG_LRTZ_GRADIENT", "OPT__FLAG_INTERFERENCE", - "OPT__FLAG_SPECTRAL", "OPT__FLAG_ANGULAR", "OPT__FLAG_RADIAL" }; + "OPT__FLAG_SPECTRAL", "OPT__FLAG_ANGULAR", "OPT__FLAG_RADIAL", + "OPT__FLAG_COOLING_LEN" }; const char FileName[][100] = { "Input__Flag_Rho", "Input__Flag_RhoGradient", "Input__Flag_PresGradient", "Input__Flag_EngyDensity", "Input__Flag_Lohner", "Input__Flag_User", "Input__Flag_NParPatch", "Input__Flag_NParCell", "Input__Flag_ParMassCell", "Input__Flag_Vorticity", "Input__Flag_Jeans", "Input__Flag_Current", "Input__Flag_CRay", "Input__Flag_LrtzGradient", "Input__Flag_Interference", "Input__Flag_Spectral", "Input__Flag_AngularResolution", - "Input__Flag_RadialResolution"}; + "Input__Flag_RadialResolution", "Input__Flag_CoolingLen" }; + double *FlagTable[NFlagMode] = { FlagTable_Rho, FlagTable_RhoGradient, FlagTable_PresGradient, NULL, NULL, NULL, NULL, NULL, FlagTable_ParMassCell, FlagTable_Vorticity, FlagTable_Jeans, FlagTable_Current, FlagTable_CRay, FlagTable_LrtzGradient, NULL, NULL, NULL, - FlagTable_Radial }; + FlagTable_Radial, FlagTable_CoolingLen }; FILE *File; char *input_line = NULL, TargetName[100]; @@ -138,6 +147,9 @@ void Init_Load_FlagCriteria() # ifdef SRHD FlagTable_LrtzGradient[lv] = -1.0; # endif +# ifdef SUPPORT_GRACKLE + FlagTable_CoolingLen [lv] = -1.0; +# endif # elif ( MODEL == ELBDM ) for (int t=0; t<2; t++) { diff --git a/src/Init/Init_Load_Parameter.cpp b/src/Init/Init_Load_Parameter.cpp index 7149d02b3a..9f31d7c863 100644 --- a/src/Init/Init_Load_Parameter.cpp +++ b/src/Init/Init_Load_Parameter.cpp @@ -130,6 +130,9 @@ void Init_Load_Parameter() # ifdef CR_DIFFUSION ReadPara->Add( "DT__CR_DIFFUSION", &DT__CR_DIFFUSION, 3.0e-1, 0.0, NoMax_double ); # endif +# ifdef SUPPORT_GRACKLE + ReadPara->Add( "DT__GRACKLE_COOLING", &DT__GRACKLE_COOLING, -1.0, NoMin_double, NoMax_double ); +# endif # ifdef COMOVING ReadPara->Add( "DT__MAX_DELTA_A", &DT__MAX_DELTA_A, 0.01, 0.0, NoMax_double ); # endif @@ -166,6 +169,9 @@ void Init_Load_Parameter() # ifdef SRHD ReadPara->Add( "OPT__FLAG_LRTZ_GRADIENT", &OPT__FLAG_LRTZ_GRADIENT, false, Useless_bool, Useless_bool ); # endif +# ifdef SUPPORT_GRACKLE + ReadPara->Add( "OPT__FLAG_COOLING_LEN", &OPT__FLAG_COOLING_LEN, false, Useless_bool, Useless_bool ); +# endif # ifdef MHD ReadPara->Add( "OPT__FLAG_CURRENT", &OPT__FLAG_CURRENT, false, Useless_bool, Useless_bool ); # endif @@ -253,6 +259,10 @@ void Init_Load_Parameter() ReadPara->Add( "GRACKLE_THREE_BODY_RATE", &GRACKLE_THREE_BODY_RATE, 0, 0, 5 ); ReadPara->Add( "GRACKLE_CIE_COOLING", &GRACKLE_CIE_COOLING, false, Useless_bool, Useless_bool ); ReadPara->Add( "GRACKLE_H2_OPA_APPROX", &GRACKLE_H2_OPA_APPROX, 0, 0, 1 ); + ReadPara->Add( "GRACKLE_USE_V_HEATING_RATE", &GRACKLE_USE_V_HEATING_RATE, false, Useless_bool, Useless_bool ); + ReadPara->Add( "GRACKLE_USE_S_HEATING_RATE", &GRACKLE_USE_S_HEATING_RATE, false, Useless_bool, Useless_bool ); + ReadPara->Add( "GRACKLE_HYDROGEN_MFRAC", &GRACKLE_HYDROGEN_MFRAC, 0.76, 0.0, 1.0 ); + ReadPara->Add( "OPT__UNFREEZE_GRACKLE", &OPT__UNFREEZE_GRACKLE, false, Useless_bool, Useless_bool ); // do not check CHE_GPU_NPGROUP since it may be reset by either Init_ResetParameter() or CUAPI_SetMemSize() ReadPara->Add( "CHE_GPU_NPGROUP", &CHE_GPU_NPGROUP, -1, NoMin_int, NoMax_int ); # endif @@ -532,7 +542,6 @@ void Init_Load_Parameter() ReadPara->Add( "OPT__OUTPUT_CS", &OPT__OUTPUT_CS, false, Useless_bool, Useless_bool ); ReadPara->Add( "OPT__OUTPUT_DIVVEL", &OPT__OUTPUT_DIVVEL, false, Useless_bool, Useless_bool ); ReadPara->Add( "OPT__OUTPUT_MACH", &OPT__OUTPUT_MACH, false, Useless_bool, Useless_bool ); -# endif # ifdef MHD ReadPara->Add( "OPT__OUTPUT_DIVMAG", &OPT__OUTPUT_DIVMAG, false, Useless_bool, Useless_bool ); # endif @@ -541,6 +550,12 @@ void Init_Load_Parameter() ReadPara->Add( "OPT__OUTPUT_3VELOCITY", &OPT__OUTPUT_3VELOCITY, false, Useless_bool, Useless_bool ); ReadPara->Add( "OPT__OUTPUT_ENTHALPY", &OPT__OUTPUT_ENTHALPY, true, Useless_bool, Useless_bool ); # endif +# ifdef SUPPORT_GRACKLE + ReadPara->Add( "OPT__OUTPUT_GRACKLE_TEMP", &OPT__OUTPUT_GRACKLE_TEMP, false, Useless_bool, Useless_bool ); + ReadPara->Add( "OPT__OUTPUT_GRACKLE_MU", &OPT__OUTPUT_GRACKLE_MU, false, Useless_bool, Useless_bool ); + ReadPara->Add( "OPT__OUTPUT_GRACKLE_TCOOL", &OPT__OUTPUT_GRACKLE_TCOOL, false, Useless_bool, Useless_bool ); +# endif +# endif // #if ( MODEL == HYDRO ) ReadPara->Add( "OPT__OUTPUT_USER_FIELD", &OPT__OUTPUT_USER_FIELD, false, Useless_bool, Useless_bool ); ReadPara->Add( "OPT__OUTPUT_MODE", &OPT__OUTPUT_MODE, -1, 1, 3 ); ReadPara->Add( "OPT__OUTPUT_RESTART", &OPT__OUTPUT_RESTART, false, Useless_bool, Useless_bool ); diff --git a/src/Init/Init_TestProb.cpp b/src/Init/Init_TestProb.cpp index 819bc5cd4a..b108f187ca 100644 --- a/src/Init/Init_TestProb.cpp +++ b/src/Init/Init_TestProb.cpp @@ -30,6 +30,7 @@ void Init_TestProb_Hydro_EnergyPowerSpectrum(); void Init_TestProb_Hydro_CR_SoundWave(); void Init_TestProb_Hydro_CR_ShockTube(); void Init_TestProb_Hydro_CR_Diffusion(); +void Init_TestProb_Hydro_GrackleTest(); void Init_TestProb_ELBDM_ExtPot(); void Init_TestProb_ELBDM_JeansInstabilityComoving(); @@ -96,6 +97,7 @@ void Init_TestProb() case TESTPROB_HYDRO_CR_SOUNDWAVE : Init_TestProb_Hydro_CR_SoundWave(); break; case TESTPROB_HYDRO_CR_SHOCKTUBE : Init_TestProb_Hydro_CR_ShockTube(); break; case TESTPROB_HYDRO_CR_DIFFUSION : Init_TestProb_Hydro_CR_Diffusion(); break; + case TESTPROB_HYDRO_GRACKLE_TEST : Init_TestProb_Hydro_GrackleTest(); break; case TESTPROB_ELBDM_EXTPOT : Init_TestProb_ELBDM_ExtPot(); break; case TESTPROB_ELBDM_JEANS_INSTABILITY_COMOVING : Init_TestProb_ELBDM_JeansInstabilityComoving(); break; diff --git a/src/Main/InvokeSolver.cpp b/src/Main/InvokeSolver.cpp index e256197cc3..3398f8ac0a 100644 --- a/src/Main/InvokeSolver.cpp +++ b/src/Main/InvokeSolver.cpp @@ -102,7 +102,12 @@ void InvokeSolver( const Solver_t TSolver, const int lv, const double TimeNew, c // reset the time-step actually adopted to zero for OPT__FREEZE_FLUID - const double dt = ( OPT__FREEZE_FLUID ) ? 0.0 : dt_in; +# ifdef SUPPORT_GRACKLE + const bool UnFreeze = ( TSolver == GRACKLE_SOLVER && OPT__UNFREEZE_GRACKLE ); +# else + const bool UnFreeze = false; +# endif + const double dt = ( OPT__FREEZE_FLUID && ! UnFreeze ) ? 0.0 : dt_in; // set the maximum number of patch groups to be updated at a time diff --git a/src/Main/Main.cpp b/src/Main/Main.cpp index 70358820de..690ab0b315 100644 --- a/src/Main/Main.cpp +++ b/src/Main/Main.cpp @@ -260,6 +260,16 @@ char GRACKLE_CLOUDY_TABLE[MAX_STRING]; int GRACKLE_THREE_BODY_RATE; bool GRACKLE_CIE_COOLING; int GRACKLE_H2_OPA_APPROX; +bool GRACKLE_USE_V_HEATING_RATE; +bool GRACKLE_USE_S_HEATING_RATE; +double GRACKLE_HYDROGEN_MFRAC; +bool OPT__UNFREEZE_GRACKLE; +bool OPT__OUTPUT_GRACKLE_TEMP; +bool OPT__OUTPUT_GRACKLE_MU; +bool OPT__OUTPUT_GRACKLE_TCOOL; +bool OPT__FLAG_COOLING_LEN; +double FlagTable_CoolingLen[NLEVEL-1]; +double DT__GRACKLE_COOLING; int CHE_GPU_NPGROUP; #endif diff --git a/src/Makefile_base b/src/Makefile_base index 05a5ff17b0..994b2cf366 100644 --- a/src/Makefile_base +++ b/src/Makefile_base @@ -324,7 +324,8 @@ ifeq "$(filter -DSUPPORT_GRACKLE, $(SIMU_OPTION))" "-DSUPPORT_GRACKLE" CPU_FILE += CPU_GrackleSolver.cpp CPU_FILE += Grackle_Init.cpp Grackle_End.cpp Init_MemAllocate_Grackle.cpp End_MemFree_Grackle.cpp \ - Grackle_Prepare.cpp Grackle_Close.cpp Grackle_Init_FieldData.cpp Grackle_AdvanceDt.cpp + Grackle_Prepare.cpp Grackle_Close.cpp Grackle_Init_FieldData.cpp Grackle_AdvanceDt.cpp \ + Grackle_Calculate.cpp Grackle_GetTimeStep_CoolingTime.cpp vpath %.cpp Grackle Grackle/CPU_Grackle endif # SUPPORT_GRACKLE diff --git a/src/Miscellaneous/Mis_GetTimeStep.cpp b/src/Miscellaneous/Mis_GetTimeStep.cpp index 322d9e545e..7e031d54b8 100644 --- a/src/Miscellaneous/Mis_GetTimeStep.cpp +++ b/src/Miscellaneous/Mis_GetTimeStep.cpp @@ -236,6 +236,20 @@ double Mis_GetTimeStep( const int lv, const double dTime_SyncFaLv, const double # endif +// 1.10 CRITERION TEN : Grackle cooling time +// ============================================================================================================= +# ifdef SUPPORT_GRACKLE + if ( DT__GRACKLE_COOLING >= 0.0 ) + { + dTime[NdTime] = dTime_dt * Grackle_GetTimeStep_CoolingTime( lv ); + sprintf( dTime_Name[NdTime++], "%s", "Grackle_TCool" ); + +// when fluid is freezed, disable this criterion by resetting it to a huge value + if ( OPT__FREEZE_FLUID && ! OPT__UNFREEZE_GRACKLE ) dTime[NdTime-1] = HUGE_NUMBER; + } +# endif + + // 2. get the minimum time-step from all criteria // ============================================================================================================= diff --git a/src/Output/Output_DumpData_Part.cpp b/src/Output/Output_DumpData_Part.cpp index 7bed3b7093..2118db4a28 100644 --- a/src/Output/Output_DumpData_Part.cpp +++ b/src/Output/Output_DumpData_Part.cpp @@ -4,9 +4,10 @@ static void WriteFile( FILE *File, const int lv, const int PID, const int i, con const int ii, const int jj, const int kk, const real (*DerField)[ CUBE(PS1) ] ); static void GetDerivedField( real (*Der_FluIn)[NCOMP_TOTAL][ CUBE(DER_NXT) ], real (*Der_Out ) [ CUBE(PS1) ], + real (*Der_Out_1PG) [ CUBE(PS2) ], real (*Der_MagFC)[NCOMP_MAG ][ (DER_NXT+1)*SQR(DER_NXT) ], real (*Der_MagCC) [ CUBE(DER_NXT) ], - const int lv, const int PID, const bool PrepFluIn ); + const int lv, const int PID, const bool PrepFluIn, const bool PrepOutPG ); @@ -111,9 +112,11 @@ void Output_DumpData_Part( const OptOutputPart_t Part, const bool BaseOnly, cons // for the derived fields const int Der_NP = 8; bool Der_PrepFluIn; + bool Der_PrepOutPG; real (*Der_FluIn)[NCOMP_TOTAL][ CUBE(DER_NXT) ] = new real [Der_NP][NCOMP_TOTAL ][ CUBE(DER_NXT) ]; real (*Der_Out ) [ CUBE(PS1) ] = new real [DER_NOUT_MAX][ CUBE(PS1) ]; + real (*Der_Out_1PG ) [ CUBE(PS2) ] = new real [DER_NOUT_MAX][ CUBE(PS2) ]; # ifdef MHD real (*Der_MagFC)[NCOMP_MAG ][ (DER_NXT+1)*SQR(DER_NXT) ] = new real [Der_NP][NCOMP_MAG ][ (DER_NXT+1)*SQR(DER_NXT) ]; real (*Der_MagCC) [ CUBE(DER_NXT) ] = new real [NCOMP_MAG ][ CUBE(DER_NXT) ]; @@ -170,6 +173,14 @@ void Output_DumpData_Part( const OptOutputPart_t Part, const bool BaseOnly, cons } if ( OPT__OUTPUT_ENTHALPY ) fprintf( File, " %*s", StrLen_Flt, "Reduced enthalpy" ); +# endif +# ifdef SUPPORT_GRACKLE + if ( OPT__OUTPUT_GRACKLE_TEMP ) + fprintf( File, " %*s", StrLen_Flt, "Grackle temperature" ); + if ( OPT__OUTPUT_GRACKLE_MU ) + fprintf( File, " %*s", StrLen_Flt, "Grackle mu" ); + if ( OPT__OUTPUT_GRACKLE_TCOOL ) + fprintf( File, " %*s", StrLen_Flt, "Grackle cooling time" ); # endif if ( OPT__OUTPUT_USER_FIELD ) { for (int v=0; vNPatchComma[lv][1]; PID++) { if ( PID % 8 == 0 ) Der_PrepFluIn = true; + if ( PID % 8 == 0 ) Der_PrepOutPG = true; // output the patch data only if it has no son (if the option "BaseOnly" is turned off) if ( amr->patch[0][lv][PID]->son == -1 || BaseOnly ) @@ -203,8 +215,9 @@ void Output_DumpData_Part( const OptOutputPart_t Part, const bool BaseOnly, cons if ( Corner[0] == Corner[1] && Corner[0] == Corner[2] ) { // compute the derived fields - GetDerivedField( Der_FluIn, Der_Out, Der_MagFC, Der_MagCC, lv, PID, Der_PrepFluIn ); + GetDerivedField( Der_FluIn, Der_Out, Der_Out_1PG, Der_MagFC, Der_MagCC, lv, PID, Der_PrepFluIn, Der_PrepOutPG ); Der_PrepFluIn = false; + Der_PrepOutPG = false; // write data for (int k=0; kz ) ) { // compute the derived fields - GetDerivedField( Der_FluIn, Der_Out, Der_MagFC, Der_MagCC, lv, PID, Der_PrepFluIn ); + GetDerivedField( Der_FluIn, Der_Out, Der_Out_1PG, Der_MagFC, Der_MagCC, lv, PID, Der_PrepFluIn, Der_PrepOutPG ); Der_PrepFluIn = false; + Der_PrepOutPG = false; // write data // --> check whether the cell is within the target range @@ -260,6 +274,7 @@ void Output_DumpData_Part( const OptOutputPart_t Part, const bool BaseOnly, cons delete [] Der_FluIn; delete [] Der_Out; + delete [] Der_Out_1PG; # ifdef MHD delete [] Der_MagFC; delete [] Der_MagCC; @@ -414,6 +429,18 @@ void WriteFile( FILE *File, const int lv, const int PID, const int i, const int } # endif +# ifdef SUPPORT_GRACKLE + if ( OPT__OUTPUT_GRACKLE_TEMP ) + fprintf( File, BlankPlusFormat_Flt, DerField[ Der_FieldIdx ++ ][Der_CellIdx] ); + + if ( OPT__OUTPUT_GRACKLE_MU ) + fprintf( File, BlankPlusFormat_Flt, DerField[ Der_FieldIdx ++ ][Der_CellIdx] ); + + if ( OPT__OUTPUT_GRACKLE_TCOOL ) + fprintf( File, BlankPlusFormat_Flt, DerField[ Der_FieldIdx ++ ][Der_CellIdx] ); + +# endif + if ( OPT__OUTPUT_USER_FIELD ) { for (int v=0; v To prepare patches within the same patch group just once +// PrepOutPG : Whether to fill in Out_1PG[] +// --> To prepare patches within the same patch group just once // -// Return : FluIn[], Out[], MagFC[] (MagCC[] is not useful outside this function) +// Return : FluIn[], Out[], Out_1PG[], MagFC[] (MagCC[] is not useful outside this function) //------------------------------------------------------------------------------------------------------- void GetDerivedField( real (*FluIn)[NCOMP_TOTAL][ CUBE(DER_NXT) ], real (*Out ) [ CUBE(PS1) ], + real (*Out_1PG) [ CUBE(PS2) ], real (*MagFC)[NCOMP_MAG ][ (DER_NXT+1)*SQR(DER_NXT) ], real (*MagCC) [ CUBE(DER_NXT) ], - const int lv, const int PID, const bool PrepFluIn ) + const int lv, const int PID, const bool PrepFluIn, const bool PrepOutPG ) { const double dh = amr->dh[lv]; @@ -526,6 +558,65 @@ void GetDerivedField( real (*FluIn)[NCOMP_TOTAL][ CUBE(DER_NXT) ], OutFieldIdx += NFieldOut; } + +# ifdef SUPPORT_GRACKLE + if ( OPT__OUTPUT_GRACKLE_TEMP ) + { + const int NFieldOut = 1; + + if ( OutFieldIdx + NFieldOut > DER_NOUT_MAX ) + Aux_Error( ERROR_INFO, "OutFieldIdx (%d) + NFieldOut (%d) > DER_NOUT_MAX (%d) !!\n", + OutFieldIdx, NFieldOut, DER_NOUT_MAX ); + +// Grackle_Calculate has to prepare the field with one patch group each time + const int PID0 = PID - LocalID; + if ( PrepOutPG ) Grackle_Calculate( Out_1PG[OutFieldIdx], _GRACKLE_TEMP, lv, 1, &PID0 ); + +// copy the corresponding patch from the already calculated patch group + const real *Out_1P = Out_1PG[OutFieldIdx] + LocalID*CUBE(PS1); + memcpy( Out[OutFieldIdx], Out_1P, sizeof(real)*CUBE(PS1) ); + + OutFieldIdx += NFieldOut; + } + + if ( OPT__OUTPUT_GRACKLE_MU ) + { + const int NFieldOut = 1; + + if ( OutFieldIdx + NFieldOut > DER_NOUT_MAX ) + Aux_Error( ERROR_INFO, "OutFieldIdx (%d) + NFieldOut (%d) > DER_NOUT_MAX (%d) !!\n", + OutFieldIdx, NFieldOut, DER_NOUT_MAX ); + +// Grackle_Calculate has to prepare the field with one patch group each time + const int PID0 = PID - LocalID; + if ( PrepOutPG ) Grackle_Calculate( Out_1PG[OutFieldIdx], _GRACKLE_MU, lv, 1, &PID0 ); + +// copy the corresponding one patch from the already calculated patch group + const real *Out_1P = Out_1PG[OutFieldIdx] + LocalID*CUBE(PS1); + memcpy( Out[OutFieldIdx], Out_1P, sizeof(real)*CUBE(PS1) ); + + OutFieldIdx += NFieldOut; + } + + if ( OPT__OUTPUT_GRACKLE_TCOOL ) + { + const int NFieldOut = 1; + + if ( OutFieldIdx + NFieldOut > DER_NOUT_MAX ) + Aux_Error( ERROR_INFO, "OutFieldIdx (%d) + NFieldOut (%d) > DER_NOUT_MAX (%d) !!\n", + OutFieldIdx, NFieldOut, DER_NOUT_MAX ); + +// Grackle_Calculate has to prepare the field with one patch group each time + const int PID0 = PID - LocalID; + if ( PrepOutPG ) Grackle_Calculate( Out_1PG[OutFieldIdx], _GRACKLE_TCOOL, lv, 1, &PID0 ); + +// copy the corresponding one patch from the already calculated patch group + const real *Out_1P = Out_1PG[OutFieldIdx] + LocalID*CUBE(PS1); + memcpy( Out[OutFieldIdx], Out_1P, sizeof(real)*CUBE(PS1) ); + + OutFieldIdx += NFieldOut; + } +# endif // #ifdef SUPPORT_GRACKLE # endif // #if ( MODEL == HYDRO ) if ( OPT__OUTPUT_USER_FIELD ) diff --git a/src/Output/Output_DumpData_Total_HDF5.cpp b/src/Output/Output_DumpData_Total_HDF5.cpp index 75462ea3e7..dc09db964e 100644 --- a/src/Output/Output_DumpData_Total_HDF5.cpp +++ b/src/Output/Output_DumpData_Total_HDF5.cpp @@ -79,7 +79,7 @@ Procedure for outputting new variables: //------------------------------------------------------------------------------------------------------- -// Function : Output_DumpData_Total_HDF5 (FormatVersion = 2505) +// Function : Output_DumpData_Total_HDF5 (FormatVersion = 2506) // Description : Output all simulation data in the HDF5 format, which can be used as a restart file // or loaded by YT // @@ -279,6 +279,10 @@ Procedure for outputting new variables: // Input__TestProb parameters in "Info/InputTest" // 2504 : 2025/04/29 --> output OPT__PAR_INIT_CHECK // 2505 : 2025/05/07 --> output PassiveFloor_Var +// 2506 : 2025/12/07 --> output GRACKLE_USE_V_HEATING_RATE, GRACKLE_USE_S_HEATING_RATE, +// GRACKLE_HYDROGEN_MFRAC, OPT__UNFREEZE_GRACKLE, +// OPT__OUTPUT_GRACKLE_TEMP, OPT__OUTPUT_GRACKLE_MU, OPT__OUTPUT_GRACKLE_TCOOL, +// DT__GRACKLE_COOLING, OPT__FLAG_COOLING_LEN, FlagTable_CoolingLen //------------------------------------------------------------------------------------------------------- void Output_DumpData_Total_HDF5( const char *FileName ) { @@ -381,7 +385,6 @@ void Output_DumpData_Total_HDF5( const char *FileName ) if ( MachDumpIdx >= NFIELD_STORED_MAX ) Aux_Error( ERROR_INFO, "exceed NFIELD_STORED_MAX (%d) !!\n", NFIELD_STORED_MAX ); if ( OPT__OUTPUT_MACH ) sprintf( FieldLabelOut[MachDumpIdx ], "%s", "Mach" ); -# endif # ifdef MHD const int DivMagDumpIdx = ( OPT__OUTPUT_DIVMAG ) ? NFieldStored++ : NoDump; @@ -413,6 +416,24 @@ void Output_DumpData_Total_HDF5( const char *FileName ) if ( OPT__OUTPUT_ENTHALPY ) sprintf( FieldLabelOut[EnthalpyDumpIdx], "%s", "Enth" ); # endif // #ifdef SRHD +# ifdef SUPPORT_GRACKLE + const int GrackleTempDumpIdx = ( OPT__OUTPUT_GRACKLE_TEMP ) ? NFieldStored++ : NoDump; + if ( GrackleTempDumpIdx >= NFIELD_STORED_MAX ) + Aux_Error( ERROR_INFO, "exceed NFIELD_STORED_MAX (%d) !!\n", NFIELD_STORED_MAX ); + if ( OPT__OUTPUT_GRACKLE_TEMP ) sprintf( FieldLabelOut[GrackleTempDumpIdx], "%s", "GrackleTemp" ); + + const int GrackleMuDumpIdx = ( OPT__OUTPUT_GRACKLE_MU ) ? NFieldStored++ : NoDump; + if ( GrackleMuDumpIdx >= NFIELD_STORED_MAX ) + Aux_Error( ERROR_INFO, "exceed NFIELD_STORED_MAX (%d) !!\n", NFIELD_STORED_MAX ); + if ( OPT__OUTPUT_GRACKLE_MU ) sprintf( FieldLabelOut[GrackleMuDumpIdx], "%s", "GrackleMu" ); + + const int GrackleTCoolDumpIdx = ( OPT__OUTPUT_GRACKLE_TCOOL ) ? NFieldStored++ : NoDump; + if ( GrackleTCoolDumpIdx >= NFIELD_STORED_MAX ) + Aux_Error( ERROR_INFO, "exceed NFIELD_STORED_MAX (%d) !!\n", NFIELD_STORED_MAX ); + if ( OPT__OUTPUT_GRACKLE_TCOOL ) sprintf( FieldLabelOut[GrackleTCoolDumpIdx], "%s", "GrackleTCool" ); +# endif // ifdef SUPPORT_GRACKLE +# endif // if ( MODEL == HYDRO ) + const int UserDumpIdx0 = ( OPT__OUTPUT_USER_FIELD ) ? NFieldStored : NoDump; if ( UserDumpIdx0+UserDerField_Num-1 >= NFIELD_STORED_MAX ) Aux_Error( ERROR_INFO, "exceed NFIELD_STORED_MAX (%d) !!\n", NFIELD_STORED_MAX ); @@ -777,6 +798,7 @@ void Output_DumpData_Total_HDF5( const char *FileName ) real (*Der_FluIn)[NCOMP_TOTAL][ CUBE(DER_NXT) ] = new real [Der_NP][NCOMP_TOTAL ][ CUBE(DER_NXT) ]; real (*Der_Out ) [ CUBE(PS1) ] = new real [DER_NOUT_MAX][ CUBE(PS1) ]; + real (*Der_Out_1PG) [ CUBE(PS2) ] = new real [DER_NOUT_MAX][ CUBE(PS2) ]; # ifdef MHD real (*Der_MagFC)[NCOMP_MAG ][ (DER_NXT+1)*SQR(DER_NXT) ] = new real [Der_NP][NCOMP_MAG ][ (DER_NXT+1)*SQR(DER_NXT) ]; real (*Der_MagCC) [ CUBE(DER_NXT) ] = new real [NCOMP_MAG ][ CUBE(DER_NXT) ]; @@ -1133,9 +1155,68 @@ void Output_DumpData_Total_HDF5( const char *FileName ) } } # endif // #ifdef SRHD + +# ifdef SUPPORT_GRACKLE +// d-11. Grackle's temperature + else if ( v == GrackleTempDumpIdx ) + { + for (int PID0=0; PID0NPatchComma[lv][1]; PID0+=8) + { +// compute the target derived field with one patch group each time + Grackle_Calculate( Der_Out_1PG[0], _GRACKLE_TEMP, lv, 1, &PID0 ); + +// store the target derived field patch by patch in a patch group + for (int LocalID=0; LocalID<8; LocalID++) + { + const int PID = PID0 + LocalID; + const real *Der_Out_1P = Der_Out_1PG[0] + LocalID*CUBE(PS1); +// copy data of one patch from Der_Out_1P[] to FieldData[] + memcpy( FieldData[PID], Der_Out_1P, FieldSizeOnePatch ); + } // for (int LocalID=0; LocalID<8; LocalID++) + } // for (int PID0=0; PID0NPatchComma[lv][1]; PID0+=8) + } // if ( v == GrackleTempDumpIdx ) + +// d-12. Grackle's mu (mean molecular weight) + else if ( v == GrackleMuDumpIdx ) + { + for (int PID0=0; PID0NPatchComma[lv][1]; PID0+=8) + { +// compute the target derived field with one patch group each time + Grackle_Calculate( Der_Out_1PG[0], _GRACKLE_MU, lv, 1, &PID0 ); + +// store the target derived field patch by patch in a patch group + for (int LocalID=0; LocalID<8; LocalID++) + { + const int PID = PID0 + LocalID; + const real *Der_Out_1P = Der_Out_1PG[0] + LocalID*CUBE(PS1); +// copy data of one patch from Der_Out_1P[] to FieldData[] + memcpy( FieldData[PID], Der_Out_1P, FieldSizeOnePatch ); + } // for (int LocalID=0; LocalID<8; LocalID++) + } // for (int PID0=0; PID0NPatchComma[lv][1]; PID0+=8) + } // if ( v == GrackleMuDumpIdx ) + +// d-13. Grackle's cooling time + else if ( v == GrackleTCoolDumpIdx ) + { + for (int PID0=0; PID0NPatchComma[lv][1]; PID0+=8) + { +// compute the target derived field with one patch group each time + Grackle_Calculate( Der_Out_1PG[0], _GRACKLE_TCOOL, lv, 1, &PID0 ); + +// store the target derived field patch by patch in a patch group + for (int LocalID=0; LocalID<8; LocalID++) + { + const int PID = PID0 + LocalID; + const real *Der_Out_1P = Der_Out_1PG[0] + LocalID*CUBE(PS1); +// copy data of one patch from Der_Out_1P[] to FieldData[] + memcpy( FieldData[PID], Der_Out_1P, FieldSizeOnePatch ); + } // for (int LocalID=0; LocalID<8; LocalID++) + } // for (int PID0=0; PID0NPatchComma[lv][1]; PID0+=8) + } // if ( v == GrackleTCoolDumpIdx ) +# endif // #ifdef SUPPORT_GRACKLE # endif // #if ( MODEL == HYDRO ) -// d-11. user-defined derived fields +// d-14. user-defined derived fields // the following check also works for OPT__OUTPUT_USER_FIELD==false since UserDerField_Num is initialized as 0 else if ( v >= UserDumpIdx0 && v < UserDumpIdx0 + UserDerField_Num ) { @@ -1324,6 +1405,7 @@ void Output_DumpData_Total_HDF5( const char *FileName ) delete [] Der_FluIn; delete [] Der_Out; + delete [] Der_Out_1PG; # ifdef MHD delete [] Der_MagFC; delete [] Der_MagCC; @@ -1667,7 +1749,7 @@ void FillIn_KeyInfo( KeyInfo_t &KeyInfo, const int NFieldStored ) const time_t CalTime = time( NULL ); // calendar time - KeyInfo.FormatVersion = 2505; + KeyInfo.FormatVersion = 2506; KeyInfo.Model = MODEL; KeyInfo.NLevel = NLEVEL; KeyInfo.NCompFluid = NCOMP_FLUID; @@ -2455,6 +2537,9 @@ void FillIn_InputPara( InputPara_t &InputPara, const int NFieldStored, char Fiel # ifdef CR_DIFFUSION InputPara.Dt__CR_Diffusion = DT__CR_DIFFUSION; # endif +# ifdef SUPPORT_GRACKLE + InputPara.Dt__GrackleCooling = DT__GRACKLE_COOLING; +# endif # ifdef COMOVING InputPara.Dt__MaxDeltaA = DT__MAX_DELTA_A; # endif @@ -2501,6 +2586,9 @@ void FillIn_InputPara( InputPara_t &InputPara, const int NFieldStored, char Fiel InputPara.Dt__SpeedOfLight = DT__SPEED_OF_LIGHT; InputPara.Opt__Flag_LrtzGradient = OPT__FLAG_LRTZ_GRADIENT; # endif +# ifdef SUPPORT_GRACKLE + InputPara.Opt__Flag_CoolingLen = OPT__FLAG_COOLING_LEN; +# endif # ifdef COSMIC_RAY InputPara.Opt__Flag_CRay = OPT__FLAG_CRAY; # endif @@ -2688,6 +2776,10 @@ void FillIn_InputPara( InputPara_t &InputPara, const int NFieldStored, char Fiel InputPara.Grackle_ThreeBodyRate = GRACKLE_THREE_BODY_RATE; InputPara.Grackle_CIE_Cooling = GRACKLE_CIE_COOLING; InputPara.Grackle_H2_OpaApprox = GRACKLE_H2_OPA_APPROX; + InputPara.Grackle_UseVHeatingRate = GRACKLE_USE_V_HEATING_RATE; + InputPara.Grackle_UseSHeatingRate = GRACKLE_USE_S_HEATING_RATE; + InputPara.Grackle_HydrogenMFrac = GRACKLE_HYDROGEN_MFRAC; + InputPara.Opt__UnfreezeGrackle = OPT__UNFREEZE_GRACKLE; InputPara.Che_GPU_NPGroup = CHE_GPU_NPGROUP; # endif @@ -2837,6 +2929,11 @@ void FillIn_InputPara( InputPara_t &InputPara, const int NFieldStored, char Fiel InputPara.Opt__Output_Lorentz = OPT__OUTPUT_LORENTZ; InputPara.Opt__Output_Enthalpy = OPT__OUTPUT_ENTHALPY; # endif +# ifdef SUPPORT_GRACKLE + InputPara.Opt__Output_GrackleTemp = OPT__OUTPUT_GRACKLE_TEMP; + InputPara.Opt__Output_GrackleMu = OPT__OUTPUT_GRACKLE_MU; + InputPara.Opt__Output_GrackleTCool = OPT__OUTPUT_GRACKLE_TCOOL; +# endif # endif // #if ( MODEL == HYDRO ) InputPara.Opt__Output_UserField = OPT__OUTPUT_USER_FIELD; InputPara.Opt__Output_Mode = OPT__OUTPUT_MODE; @@ -2931,6 +3028,9 @@ void FillIn_InputPara( InputPara_t &InputPara, const int NFieldStored, char Fiel # ifdef SRHD InputPara.FlagTable_LrtzGradient[lv] = FlagTable_LrtzGradient[lv]; # endif +# ifdef SUPPORT_GRACKLE + InputPara.FlagTable_CoolingLen [lv] = FlagTable_CoolingLen [lv]; +# endif # ifdef COSMIC_RAY InputPara.FlagTable_CRay [lv] = FlagTable_CRay [lv]; # endif @@ -3514,6 +3614,9 @@ void GetCompound_InputPara( hid_t &H5_TypeID, const int NFieldStored ) # ifdef CR_DIFFUSION H5Tinsert( H5_TypeID, "Dt__CR_Diffusion", HOFFSET(InputPara_t,Dt__CR_Diffusion ), H5T_NATIVE_DOUBLE ); # endif +# ifdef SUPPORT_GRACKLE + H5Tinsert( H5_TypeID, "Dt__GrackleCooling", HOFFSET(InputPara_t,Dt__GrackleCooling ), H5T_NATIVE_DOUBLE ); +# endif # ifdef COMOVING H5Tinsert( H5_TypeID, "Dt__MaxDeltaA", HOFFSET(InputPara_t,Dt__MaxDeltaA ), H5T_NATIVE_DOUBLE ); # endif @@ -3552,6 +3655,9 @@ void GetCompound_InputPara( hid_t &H5_TypeID, const int NFieldStored ) # ifdef SRHD H5Tinsert( H5_TypeID, "Opt__Flag_LrtzGradient", HOFFSET(InputPara_t,Opt__Flag_LrtzGradient ), H5T_NATIVE_INT ); # endif +# ifdef SUPPORT_GRACKLE + H5Tinsert( H5_TypeID, "Opt__Flag_CoolingLen", HOFFSET(InputPara_t,Opt__Flag_CoolingLen ), H5T_NATIVE_INT ); +# endif # ifdef COSMIC_RAY H5Tinsert( H5_TypeID, "Opt__Flag_CRay", HOFFSET(InputPara_t,Opt__Flag_CRay ), H5T_NATIVE_INT ); # endif @@ -3755,6 +3861,10 @@ void GetCompound_InputPara( hid_t &H5_TypeID, const int NFieldStored ) H5Tinsert( H5_TypeID, "Grackle_ThreeBodyRate", HOFFSET(InputPara_t,Grackle_ThreeBodyRate ), H5T_NATIVE_INT ); H5Tinsert( H5_TypeID, "Grackle_CIE_Cooling", HOFFSET(InputPara_t,Grackle_CIE_Cooling ), H5T_NATIVE_INT ); H5Tinsert( H5_TypeID, "Grackle_H2_OpaApprox", HOFFSET(InputPara_t,Grackle_H2_OpaApprox ), H5T_NATIVE_INT ); + H5Tinsert( H5_TypeID, "Grackle_UseVHeatingRate", HOFFSET(InputPara_t,Grackle_UseVHeatingRate), H5T_NATIVE_INT ); + H5Tinsert( H5_TypeID, "Grackle_UseSHeatingRate", HOFFSET(InputPara_t,Grackle_UseSHeatingRate), H5T_NATIVE_INT ); + H5Tinsert( H5_TypeID, "Grackle_HydrogenMFrac", HOFFSET(InputPara_t,Grackle_HydrogenMFrac ), H5T_NATIVE_DOUBLE ); + H5Tinsert( H5_TypeID, "Opt__UnfreezeGrackle", HOFFSET(InputPara_t,Opt__UnfreezeGrackle ), H5T_NATIVE_INT ); H5Tinsert( H5_TypeID, "Che_GPU_NPGroup", HOFFSET(InputPara_t,Che_GPU_NPGroup ), H5T_NATIVE_INT ); # endif @@ -3887,6 +3997,11 @@ void GetCompound_InputPara( hid_t &H5_TypeID, const int NFieldStored ) H5Tinsert( H5_TypeID, "Opt__Output_Lorentz", HOFFSET(InputPara_t,Opt__Output_Lorentz ), H5T_NATIVE_INT ); H5Tinsert( H5_TypeID, "Opt__Output_Enthalpy", HOFFSET(InputPara_t,Opt__Output_Enthalpy ), H5T_NATIVE_INT ); # endif +# ifdef SUPPORT_GRACKLE + H5Tinsert( H5_TypeID, "Opt__Output_GrackleTemp", HOFFSET(InputPara_t,Opt__Output_GrackleTemp ), H5T_NATIVE_INT ); + H5Tinsert( H5_TypeID, "Opt__Output_GrackleMu", HOFFSET(InputPara_t,Opt__Output_GrackleMu ), H5T_NATIVE_INT ); + H5Tinsert( H5_TypeID, "Opt__Output_GrackleTCool", HOFFSET(InputPara_t,Opt__Output_GrackleTCool ), H5T_NATIVE_INT ); +# endif # endif // #if ( MODEL == HYDRO ) H5Tinsert( H5_TypeID, "Opt__Output_UserField", HOFFSET(InputPara_t,Opt__Output_UserField ), H5T_NATIVE_INT ); H5Tinsert( H5_TypeID, "Opt__Output_Mode", HOFFSET(InputPara_t,Opt__Output_Mode ), H5T_NATIVE_INT ); @@ -3980,6 +4095,9 @@ void GetCompound_InputPara( hid_t &H5_TypeID, const int NFieldStored ) # ifdef SRHD H5Tinsert( H5_TypeID, "FlagTable_LrtzGradient", HOFFSET(InputPara_t,FlagTable_LrtzGradient ), H5_TypeID_Arr_NLvM1Double ); # endif +# ifdef SUPPORT_GRACKLE + H5Tinsert( H5_TypeID, "FlagTable_CoolingLen", HOFFSET(InputPara_t,FlagTable_CoolingLen ), H5_TypeID_Arr_NLvM1Double ); +# endif # ifdef COSMIC_RAY H5Tinsert( H5_TypeID, "FlagTable_CRay", HOFFSET(InputPara_t,FlagTable_CRay ), H5_TypeID_Arr_NLvM1Double ); # endif diff --git a/src/Refine/Flag_Check.cpp b/src/Refine/Flag_Check.cpp index ebd6801851..2ab947ae8b 100644 --- a/src/Refine/Flag_Check.cpp +++ b/src/Refine/Flag_Check.cpp @@ -34,6 +34,7 @@ static bool Check_Radial( const int i, const int j, const int k, const int lv, c // Vel : Input velocity array // Pres : Input pressure array // Lrtz : Input Lorentz factor array +// LCool : Input cooling length array // Lohner_Ave : Input array storing the averages for the Lohner error estimator // Lohner_Slope : Input array storing the slopes for the Lohner error estimator // Lohner_NVar : Number of variables stored in Lohner_Ave and Lohner_Slope @@ -52,6 +53,7 @@ static bool Check_Radial( const int i, const int j, const int k, const int lv, c bool Flag_Check( const int lv, const int PID, const int i, const int j, const int k, const real dv, const real Fluid[][PS1][PS1][PS1], const real Pot[][PS1][PS1], const real MagCC[][PS1][PS1][PS1], const real Vel[][PS1][PS1][PS1], const real Pres[][PS1][PS1], const real Lrtz[][PS1][PS1], + const real LCool[][PS1][PS1], const real *Lohner_Var, const real *Lohner_Ave, const real *Lohner_Slope, const int Lohner_NVar, const real ParCount[][PS1][PS1], const real ParDens[][PS1][PS1], const real JeansCoeff, const real *Interf_Var, const real Spectral_Cond ) @@ -181,6 +183,17 @@ bool Flag_Check( const int lv, const int PID, const int i, const int j, const in # endif +// check cooling length +// =========================================================================================== +# if ( MODEL == HYDRO && defined SUPPORT_GRACKLE ) + if ( OPT__FLAG_COOLING_LEN ) + { + Flag |= ( amr->dh[lv]*FlagTable_CoolingLen[lv] > LCool[k][j][i] ); + if ( Flag ) return Flag; + } +# endif + + // check vorticity // =========================================================================================== # if ( MODEL == HYDRO ) diff --git a/src/Refine/Flag_Real.cpp b/src/Refine/Flag_Real.cpp index 2531fdfa0c..2348be9783 100644 --- a/src/Refine/Flag_Real.cpp +++ b/src/Refine/Flag_Real.cpp @@ -192,6 +192,7 @@ void Flag_Real( const int lv, const UseLBFunc_t UseLBFunc ) real (*Pres)[PS1][PS1] = NULL; real (*Cs2)[PS1][PS1] = NULL; real (*Lrtz)[PS1][PS1] = NULL; + real (*LCool)[PS1][PS1] = NULL; real (*ParCount)[PS1][PS1] = NULL; // declare as **real** to be consistent with Par_MassAssignment() real (*ParDens )[PS1][PS1] = NULL; real *Lohner_Var = NULL; // array storing the variables for Lohner @@ -200,6 +201,7 @@ void Flag_Real( const int lv, const UseLBFunc_t UseLBFunc ) real *Interf_Var = NULL; // array storing the density and phase for the interference criterion real *Spectral_Var = NULL; // array storing a patch group of real and imaginary parts for the spectral criterion real Spectral_Cond = 0.0; // variable storing the magnitude of the largest coefficient for the spectral criterion + real *Grackle_TCool = NULL; // array storing a patch group of grackle cooling time int i_start, i_end, j_start, j_end, k_start, k_end, SibID, SibPID, PID; bool ProperNesting, NextPatch; @@ -212,12 +214,19 @@ void Flag_Real( const int lv, const UseLBFunc_t UseLBFunc ) if ( OPT__FLAG_JEANS ) NeedPres = true; if ( OPT__FLAG_JEANS ) NeedCs2 = true; # endif +# ifdef SUPPORT_GRACKLE + if ( OPT__FLAG_COOLING_LEN ) NeedCs2 = true; +# endif + if ( NeedCs2 ) NeedPres = true; # ifdef MHD if ( OPT__FLAG_CURRENT || NeedPres ) MagCC = new real [3][PS1][PS1][PS1]; # endif # ifdef SRHD if ( OPT__FLAG_LRTZ_GRADIENT ) Lrtz = new real [PS1][PS1][PS1]; +# endif +# ifdef SUPPORT_GRACKLE + if ( OPT__FLAG_COOLING_LEN ) LCool = new real [PS1][PS1][PS1]; # endif if ( OPT__FLAG_VORTICITY ) Vel = new real [3][PS1][PS1][PS1]; if ( NeedPres ) Pres = new real [PS1][PS1][PS1]; @@ -246,6 +255,11 @@ void Flag_Real( const int lv, const UseLBFunc_t UseLBFunc ) Lohner_Slope = new real [ 3*Lohner_NVar*CUBE(Lohner_NSlope) ]; // 3: X/Y/Z of 1 patch } +# ifdef SUPPORT_GRACKLE + if ( OPT__FLAG_COOLING_LEN ) + Grackle_TCool = new real [ CUBE(PS2) ]; // prepare one patch group +# endif + // loop over all REAL patches (the buffer patches will be flagged only due to the FlagBuf // extension or the grandson check) @@ -278,6 +292,11 @@ void Flag_Real( const int lv, const UseLBFunc_t UseLBFunc ) MinDens, MinPres, MinTemp, MinEntr, DE_Consistency_No ); # endif +# ifdef SUPPORT_GRACKLE + if ( OPT__FLAG_COOLING_LEN ) + Grackle_Calculate( Grackle_TCool, _GRACKLE_TCOOL, lv, NPG, &PID0 ); +# endif + // loop over all local patches within the same patch group for (int LocalID=0; LocalID<8; LocalID++) { @@ -454,6 +473,22 @@ void Flag_Real( const int lv, const UseLBFunc_t UseLBFunc ) } // i,j,k } // if ( OPT__FLAG_LRTZ_GRADIENT ) # endif // #ifdef SRHD + +# ifdef SUPPORT_GRACKLE +// evaluate cooling length + if ( OPT__FLAG_COOLING_LEN ) + { + for (int k=0; k= 1 only) [0.0] +static double GrackleTest_MFrac_HI; // HI mass fraction (GRACKLE_PRIMORDIAL >= 1 only) [0.750158] +static double GrackleTest_MFrac_HII; // HII mass fraction (GRACKLE_PRIMORDIAL >= 1 only) [0.0] +static double GrackleTest_MFrac_HeI; // HeI mass fraction (GRACKLE_PRIMORDIAL >= 1 only) [0.236892] +static double GrackleTest_MFrac_HeII; // HeII mass fraction (GRACKLE_PRIMORDIAL >= 1 only) [0.0] +static double GrackleTest_MFrac_HeIII; // HeIII mass fraction (GRACKLE_PRIMORDIAL >= 1 only) [0.0] +static double GrackleTest_MFrac_HM; // HM mass fraction (GRACKLE_PRIMORDIAL >= 2 only) [0.0] +static double GrackleTest_MFrac_H2I; // H2I mass fraction (GRACKLE_PRIMORDIAL >= 2 only) [0.0] +static double GrackleTest_MFrac_H2II; // H2II mass fraction (GRACKLE_PRIMORDIAL >= 2 only) [0.0] +static double GrackleTest_MFrac_DI; // DI mass fraction (GRACKLE_PRIMORDIAL >= 3 only) [0.0] +static double GrackleTest_MFrac_DII; // DII mass fraction (GRACKLE_PRIMORDIAL >= 3 only) [0.0] +static double GrackleTest_MFrac_HDI; // HDI mass fraction (GRACKLE_PRIMORDIAL >= 3 only) [0.0] +static double GrackleTest_HeatingRate; // User-provided heating rate (in erg cm^-3 s^-1 n_H^-1) [0.0] +static double GrackleTest_CoolingRate; // User-provided cooling rate (in erg cm^-3 s^-1 n_H^-2) [0.0] + +static double GrackleTest_logDens_Min; // Minimum log( mass density ) in the box +static double GrackleTest_logDens_Max; // Maximum log( mass density ) in the box +static double GrackleTest_logDens_Range; // Range of log ( mass density ) +static double GrackleTest_logTemp_Min; // Minimum log( temperature ) in the box +static double GrackleTest_logTemp_Max; // Maximum log( temperature ) in the box +static double GrackleTest_logTemp_Range; // Range of log ( temperature ) +// ======================================================================================= + + + + +//------------------------------------------------------------------------------------------------------- +// 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 ); + + +// errors +# if ( MODEL != HYDRO ) + Aux_Error( ERROR_INFO, "MODEL != HYDRO !!\n" ); +# endif + +# ifndef SUPPORT_GRACKLE + Aux_Error( ERROR_INFO, "SUPPORT_GRACKLE must be enabled !!\n" ); +# endif + +# ifdef PARTICLE + Aux_Error( ERROR_INFO, "PARTICLE must be disabled !!\n" ); +# endif + +# ifdef GRAVITY + Aux_Error( ERROR_INFO, "GRAVITY must be disabled !!\n" ); +# endif + +# ifdef COMOVING + Aux_Error( ERROR_INFO, "COMOVING must be disabled !!\n" ); +# endif + + if ( !OPT__UNIT ) + Aux_Error( ERROR_INFO, "OPT__UNIT must be enabled for this test !!\n" ); + + if ( !OPT__FREEZE_FLUID ) + Aux_Error( ERROR_INFO, "OPT__FREEZE_FLUID must be enabled for this test !!\n" ); + +# ifdef SUPPORT_GRACKLE + if ( !OPT__UNFREEZE_GRACKLE ) + Aux_Error( ERROR_INFO, "OPT__UNFREEZE_GRACKLE must be enabled for this test !!\n" ); + + if ( !GRACKLE_ACTIVATE ) + Aux_Error( ERROR_INFO, "GRACKLE_ACTIVATE must be enabled for this test !!\n" ); + + if ( !OPT__OUTPUT_GRACKLE_TEMP ) + Aux_Error( ERROR_INFO, "OPT__OUTPUT_GRACKLE_TEMP must be enabled for this test !!\n" ); + + if ( !OPT__OUTPUT_GRACKLE_MU ) + Aux_Error( ERROR_INFO, "OPT__OUTPUT_GRACKLE_MU must be enabled for this test !!\n" ); + + if ( !OPT__OUTPUT_GRACKLE_TCOOL ) + Aux_Error( ERROR_INFO, "OPT__OUTPUT_GRACKLE_TCOOL must be enabled for this test !!\n" ); + + if ( END_STEP != 0 && END_T != 0.0 && !GRACKLE_COOLING ) + Aux_Error( ERROR_INFO, "GRACKLE_COOLING must be enabled for time evolution in this test !!\n" ); + + + if ( MPI_Rank == 0 ) + { + if ( DT__GRACKLE_COOLING < 0.0 ) + Aux_Message( stderr, "WARNING : it's recommended to set DT__GRACKLE_COOLING for this test !!\n" ); + + if ( !OPT__FLAG_COOLING_LEN ) + Aux_Message( stderr, "WARNING : it's recommended to set OPT__FLAG_COOLING_LEN for this test !!\n" ); + } // if ( MPI_Rank == 0 ) +# endif // #ifdef SUPPORT_GRACKLE + + + if ( MPI_Rank == 0 ) Aux_Message( stdout, " Validating test problem %d ... done\n", TESTPROB_ID ); + +} // FUNCTION : Validate + + + +#if ( MODEL == HYDRO && defined SUPPORT_GRACKLE ) +//------------------------------------------------------------------------------------------------------- +// Function : LoadInputTestProb +// Description : Read problem-specific runtime parameters from Input__TestProb and store them in HDF5 snapshots (Data_*) +// +// Note : 1. Invoked by SetParameter() to read parameters +// 2. Invoked by Output_DumpData_Total_HDF5() using the function pointer Output_HDF5_InputTest_Ptr to store parameters +// 3. If there is no problem-specific runtime parameter to load, add at least one parameter +// to prevent an empty structure in HDF5_Output_t +// --> Example: +// LOAD_PARA( load_mode, "TestProb_ID", &TESTPROB_ID, TESTPROB_ID, TESTPROB_ID, TESTPROB_ID ); +// +// Parameter : load_mode : Mode for loading parameters +// --> LOAD_READPARA : Read parameters from Input__TestProb +// LOAD_HDF5_OUTPUT : Store parameters in HDF5 snapshots +// ReadPara : Data structure for reading parameters (used with LOAD_READPARA) +// HDF5_InputTest : Data structure for storing parameters in HDF5 snapshots (used with LOAD_HDF5_OUTPUT) +// +// Return : None +//------------------------------------------------------------------------------------------------------- +void LoadInputTestProb( const LoadParaMode_t load_mode, ReadPara_t *ReadPara, HDF5_Output_t *HDF5_InputTest ) +{ + +# ifndef SUPPORT_HDF5 + if ( load_mode == LOAD_HDF5_OUTPUT ) Aux_Error( ERROR_INFO, "please turn on SUPPORT_HDF5 in the Makefile for load_mode == LOAD_HDF5_OUTPUT !!\n" ); +# endif + + if ( load_mode == LOAD_READPARA && ReadPara == NULL ) Aux_Error( ERROR_INFO, "load_mode == LOAD_READPARA and ReadPara == NULL !!\n" ); + if ( load_mode == LOAD_HDF5_OUTPUT && HDF5_InputTest == NULL ) Aux_Error( ERROR_INFO, "load_mode == LOAD_HDF5_OUTPUT and HDF5_InputTest == NULL !!\n" ); + +// add parameters in the following format: +// --> note that VARIABLE, DEFAULT, MIN, and MAX must have the same data type +// --> some handy constants (e.g., NoMin_int, Eps_float, ...) are defined in "include/ReadPara.h" +// --> LOAD_PARA() is defined in "include/TestProb.h" +// ******************************************************************************************************************************** +// LOAD_PARA( load_mode, "KEY_IN_THE_FILE", &VARIABLE, DEFAULT, MIN, MAX ); +// ******************************************************************************************************************************** + LOAD_PARA( load_mode, "GrackleTest_DefaultTestMode", &GrackleTest_DefaultTestMode, 0, 0, 3 ); + LOAD_PARA( load_mode, "GrackleTest_MassDensity_Min", &GrackleTest_MassDensity_Min, 1.0e-29, Eps_double, NoMax_double ); + LOAD_PARA( load_mode, "GrackleTest_MassDensity_Max", &GrackleTest_MassDensity_Max, 1.0e-21, Eps_double, NoMax_double ); + LOAD_PARA( load_mode, "GrackleTest_TempOverMMW_Min", &GrackleTest_TempOverMMW_Min, 1.0e+00, Eps_double, NoMax_double ); + LOAD_PARA( load_mode, "GrackleTest_TempOverMMW_Max", &GrackleTest_TempOverMMW_Max, 1.0e+08, Eps_double, NoMax_double ); + LOAD_PARA( load_mode, "GrackleTest_MFrac_Metal", &GrackleTest_MFrac_Metal, 0.01295, 0.0, 1.0 ); + LOAD_PARA( load_mode, "GrackleTest_MFrac_e", &GrackleTest_MFrac_e, 0.0, 0.0, 1.0 ); + LOAD_PARA( load_mode, "GrackleTest_MFrac_HI", &GrackleTest_MFrac_HI, 0.750158, 0.0, 1.0 ); + LOAD_PARA( load_mode, "GrackleTest_MFrac_HII", &GrackleTest_MFrac_HII, 0.0, 0.0, 1.0 ); + LOAD_PARA( load_mode, "GrackleTest_MFrac_HeI", &GrackleTest_MFrac_HeI, 0.236892, 0.0, 1.0 ); + LOAD_PARA( load_mode, "GrackleTest_MFrac_HeII", &GrackleTest_MFrac_HeII, 0.0, 0.0, 1.0 ); + LOAD_PARA( load_mode, "GrackleTest_MFrac_HeIII", &GrackleTest_MFrac_HeIII, 0.0, 0.0, 1.0 ); + LOAD_PARA( load_mode, "GrackleTest_MFrac_HM", &GrackleTest_MFrac_HM, 0.0, 0.0, 1.0 ); + LOAD_PARA( load_mode, "GrackleTest_MFrac_H2I", &GrackleTest_MFrac_H2I, 0.0, 0.0, 1.0 ); + LOAD_PARA( load_mode, "GrackleTest_MFrac_H2II", &GrackleTest_MFrac_H2II, 0.0, 0.0, 1.0 ); + LOAD_PARA( load_mode, "GrackleTest_MFrac_DI", &GrackleTest_MFrac_DI, 0.0, 0.0, 1.0 ); + LOAD_PARA( load_mode, "GrackleTest_MFrac_DII", &GrackleTest_MFrac_DII, 0.0, 0.0, 1.0 ); + LOAD_PARA( load_mode, "GrackleTest_MFrac_HDI", &GrackleTest_MFrac_HDI, 0.0, 0.0, 1.0 ); + LOAD_PARA( load_mode, "GrackleTest_HeatingRate", &GrackleTest_HeatingRate, 0.0, 0.0, NoMax_double ); + LOAD_PARA( load_mode, "GrackleTest_CoolingRate", &GrackleTest_CoolingRate, 0.0, 0.0, NoMax_double ); + +} // FUNCITON : LoadInputTestProb + + + +//------------------------------------------------------------------------------------------------------- +// 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 +// (1-1) read parameters from Input__TestProb + const char FileName[] = "Input__TestProb"; + ReadPara_t *ReadPara = new ReadPara_t; + + LoadInputTestProb( LOAD_READPARA, ReadPara, NULL ); + + ReadPara->Read( FileName ); + + delete ReadPara; + +// (1-2) set the default values + if ( ! GRACKLE_METAL ) + { + GrackleTest_MFrac_Metal = 0.0; PRINT_RESET_PARA( GrackleTest_MFrac_Metal, FORMAT_REAL, "for GRACKLE_METAL disabled" ); + } + + if ( GRACKLE_PRIMORDIAL < GRACKLE_PRI_CHE_NSPE12 ) + { + GrackleTest_MFrac_DI = 0.0; PRINT_RESET_PARA( GrackleTest_MFrac_DI, FORMAT_REAL, "for GRACKLE_PRIMORDIAL < GRACKLE_PRI_CHE_NSPE12" ); + GrackleTest_MFrac_DII = 0.0; PRINT_RESET_PARA( GrackleTest_MFrac_DII, FORMAT_REAL, "for GRACKLE_PRIMORDIAL < GRACKLE_PRI_CHE_NSPE12" ); + GrackleTest_MFrac_HDI = 0.0; PRINT_RESET_PARA( GrackleTest_MFrac_HDI, FORMAT_REAL, "for GRACKLE_PRIMORDIAL < GRACKLE_PRI_CHE_NSPE12" ); + } + + if ( GRACKLE_PRIMORDIAL < GRACKLE_PRI_CHE_NSPE9 ) + { + GrackleTest_MFrac_HM = 0.0; PRINT_RESET_PARA( GrackleTest_MFrac_HM, FORMAT_REAL, "for GRACKLE_PRIMORDIAL < GRACKLE_PRI_CHE_NSPE9" ); + GrackleTest_MFrac_H2I = 0.0; PRINT_RESET_PARA( GrackleTest_MFrac_H2I, FORMAT_REAL, "for GRACKLE_PRIMORDIAL < GRACKLE_PRI_CHE_NSPE9" ); + GrackleTest_MFrac_H2II = 0.0; PRINT_RESET_PARA( GrackleTest_MFrac_H2II, FORMAT_REAL, "for GRACKLE_PRIMORDIAL < GRACKLE_PRI_CHE_NSPE9" ); + } + + if ( GRACKLE_PRIMORDIAL < GRACKLE_PRI_CHE_NSPE6 ) + { + GrackleTest_MFrac_e = 0.0; PRINT_RESET_PARA( GrackleTest_MFrac_e, FORMAT_REAL, "for GRACKLE_PRIMORDIAL < GRACKLE_PRI_CHE_NSPE6" ); + GrackleTest_MFrac_HI = 0.0; PRINT_RESET_PARA( GrackleTest_MFrac_HI, FORMAT_REAL, "for GRACKLE_PRIMORDIAL < GRACKLE_PRI_CHE_NSPE6" ); + GrackleTest_MFrac_HII = 0.0; PRINT_RESET_PARA( GrackleTest_MFrac_HII, FORMAT_REAL, "for GRACKLE_PRIMORDIAL < GRACKLE_PRI_CHE_NSPE6" ); + GrackleTest_MFrac_HeI = 0.0; PRINT_RESET_PARA( GrackleTest_MFrac_HeI, FORMAT_REAL, "for GRACKLE_PRIMORDIAL < GRACKLE_PRI_CHE_NSPE6" ); + GrackleTest_MFrac_HeII = 0.0; PRINT_RESET_PARA( GrackleTest_MFrac_HeII, FORMAT_REAL, "for GRACKLE_PRIMORDIAL < GRACKLE_PRI_CHE_NSPE6" ); + GrackleTest_MFrac_HeIII = 0.0; PRINT_RESET_PARA( GrackleTest_MFrac_HeIII, FORMAT_REAL, "for GRACKLE_PRIMORDIAL < GRACKLE_PRI_CHE_NSPE6" ); + } + + if ( !GRACKLE_USE_V_HEATING_RATE ) + { + GrackleTest_HeatingRate = 0.0; PRINT_RESET_PARA( GrackleTest_HeatingRate, FORMAT_REAL, "for GRACKLE_USE_V_HEATING_RATE disabled" ); + GrackleTest_CoolingRate = 0.0; PRINT_RESET_PARA( GrackleTest_CoolingRate, FORMAT_REAL, "for GRACKLE_USE_V_HEATING_RATE disabled" ); + } + + if ( GrackleTest_DefaultTestMode == 1 ) + { + GrackleTest_MassDensity_Min = 1.0e-29; PRINT_RESET_PARA( GrackleTest_MassDensity_Min, FORMAT_REAL, "for GrackleTest_DefaultTestMode == 1" ); + GrackleTest_MassDensity_Max = 1.0e-21; PRINT_RESET_PARA( GrackleTest_MassDensity_Max, FORMAT_REAL, "for GrackleTest_DefaultTestMode == 1" ); + GrackleTest_TempOverMMW_Min = 1.0e+00; PRINT_RESET_PARA( GrackleTest_TempOverMMW_Min, FORMAT_REAL, "for GrackleTest_DefaultTestMode == 1" ); + GrackleTest_TempOverMMW_Max = 1.0e+08; PRINT_RESET_PARA( GrackleTest_TempOverMMW_Max, FORMAT_REAL, "for GrackleTest_DefaultTestMode == 1" ); + GrackleTest_HeatingRate = 0.0; PRINT_RESET_PARA( GrackleTest_HeatingRate, FORMAT_REAL, "for GrackleTest_DefaultTestMode == 1" ); + GrackleTest_CoolingRate = 0.0; PRINT_RESET_PARA( GrackleTest_CoolingRate, FORMAT_REAL, "for GrackleTest_DefaultTestMode == 1" ); + } + else if ( GrackleTest_DefaultTestMode == 2 ) + { + GrackleTest_MassDensity_Min = 1.0e-24; PRINT_RESET_PARA( GrackleTest_MassDensity_Min, FORMAT_REAL, "for GrackleTest_DefaultTestMode == 2" ); + GrackleTest_MassDensity_Max = 1.0e-24; PRINT_RESET_PARA( GrackleTest_MassDensity_Max, FORMAT_REAL, "for GrackleTest_DefaultTestMode == 2" ); + GrackleTest_TempOverMMW_Min = 1.0e+04; PRINT_RESET_PARA( GrackleTest_TempOverMMW_Min, FORMAT_REAL, "for GrackleTest_DefaultTestMode == 2" ); + GrackleTest_TempOverMMW_Max = 1.0e+04; PRINT_RESET_PARA( GrackleTest_TempOverMMW_Max, FORMAT_REAL, "for GrackleTest_DefaultTestMode == 2" ); + GrackleTest_HeatingRate = 0.0; PRINT_RESET_PARA( GrackleTest_HeatingRate, FORMAT_REAL, "for GrackleTest_DefaultTestMode == 2" ); + GrackleTest_CoolingRate = 0.0; PRINT_RESET_PARA( GrackleTest_CoolingRate, FORMAT_REAL, "for GrackleTest_DefaultTestMode == 2" ); + } + else if ( GrackleTest_DefaultTestMode == 3 ) + { + if ( GRACKLE_PRIMORDIAL != GRACKLE_PRI_CHE_CLOUDY ) + Aux_Error( ERROR_INFO, "GRACKLE_PRIMORDIAL must be 0 for GrackleTest_DefaultTestMode = %d !!\n", + GrackleTest_DefaultTestMode ); + + if ( GRACKLE_COOLING != 1 ) + Aux_Error( ERROR_INFO, "GRACKLE_COOLING must be 1 for GrackleTest_DefaultTestMode = %d !!\n", + GrackleTest_DefaultTestMode ); + + if ( GRACKLE_USE_V_HEATING_RATE != 1 ) + Aux_Error( ERROR_INFO, "GRACKLE_USE_V_HEATING_RATE must be 1 for GrackleTest_DefaultTestMode = %d !!\n", + GrackleTest_DefaultTestMode ); + + if ( GRACKLE_UV != 0 ) + Aux_Error( ERROR_INFO, "GRACKLE_UV must be 0 for GrackleTest_DefaultTestMode = %d !!\n", + GrackleTest_DefaultTestMode ); + + if ( GRACKLE_CMB_FLOOR != 0 ) + Aux_Error( ERROR_INFO, "GRACKLE_CMB_FLOOR must be 0 for GrackleTest_DefaultTestMode = %d !!\n", + GrackleTest_DefaultTestMode ); + + if ( GRACKLE_PE_HEATING != 0 ) + Aux_Error( ERROR_INFO, "GRACKLE_PE_HEATING must be 0 for GrackleTest_DefaultTestMode = %d !!\n", + GrackleTest_DefaultTestMode ); + + GrackleTest_MassDensity_Min = 1.0e-28; PRINT_RESET_PARA( GrackleTest_MassDensity_Min, FORMAT_REAL, "for GrackleTest_DefaultTestMode == 3" ); + GrackleTest_MassDensity_Max = 1.0e-28; PRINT_RESET_PARA( GrackleTest_MassDensity_Max, FORMAT_REAL, "for GrackleTest_DefaultTestMode == 3" ); + GrackleTest_TempOverMMW_Min = 1.0e+06; PRINT_RESET_PARA( GrackleTest_TempOverMMW_Min, FORMAT_REAL, "for GrackleTest_DefaultTestMode == 3" ); + GrackleTest_TempOverMMW_Max = 1.0e+06; PRINT_RESET_PARA( GrackleTest_TempOverMMW_Max, FORMAT_REAL, "for GrackleTest_DefaultTestMode == 3" ); + GrackleTest_MFrac_Metal = 0.0; PRINT_RESET_PARA( GrackleTest_MFrac_Metal, FORMAT_REAL, "for GrackleTest_DefaultTestMode == 3" ); + GrackleTest_HeatingRate = 1.0e-24; PRINT_RESET_PARA( GrackleTest_HeatingRate, FORMAT_REAL, "for GrackleTest_DefaultTestMode == 3" ); + GrackleTest_CoolingRate = 1.6e-20; PRINT_RESET_PARA( GrackleTest_CoolingRate, FORMAT_REAL, "for GrackleTest_DefaultTestMode == 3" ); + } + else + { + Aux_Error( ERROR_INFO, "Unknown GrackleTest_DefaultTestMode = %d !!\n", GrackleTest_DefaultTestMode ); + } + +// (1-3) check the runtime parameters + if ( GrackleTest_MassDensity_Min > GrackleTest_MassDensity_Max ) + Aux_Error( ERROR_INFO, "MassDensity_Min = %14.7e > MassDensity_Max = %14.7e !!\n", + GrackleTest_MassDensity_Min, GrackleTest_MassDensity_Max ); + + if ( GrackleTest_TempOverMMW_Min > GrackleTest_TempOverMMW_Max ) + Aux_Error( ERROR_INFO, "TempOverMMW_Min = %14.7e > TempOverMMW_Max = %14.7e !!\n", + GrackleTest_TempOverMMW_Min, GrackleTest_TempOverMMW_Max ); + + if ( GRACKLE_PRIMORDIAL >= GRACKLE_PRI_CHE_NSPE6 ) + { + const double MFrac_Total = GrackleTest_MFrac_e + GrackleTest_MFrac_HI + GrackleTest_MFrac_HII + + GrackleTest_MFrac_HeI + GrackleTest_MFrac_HeII + GrackleTest_MFrac_HeIII + + GrackleTest_MFrac_HM + GrackleTest_MFrac_H2I + GrackleTest_MFrac_H2II + + GrackleTest_MFrac_DI + GrackleTest_MFrac_DII + GrackleTest_MFrac_HDI + + GrackleTest_MFrac_Metal; + + if ( ! Mis_CompareRealValue( MFrac_Total, 1.0, NULL, false ) ) + Aux_Error( ERROR_INFO, "Sum of mass fraction = %14.7e != 1.0 !!\n", MFrac_Total ); + } + + + +// (2) set the problem-specific derived parameters +// convert to code units + GrackleTest_MassDensity_Min /= UNIT_D; + GrackleTest_MassDensity_Max /= UNIT_D; + + GrackleTest_logDens_Min = log10( GrackleTest_MassDensity_Min ); + GrackleTest_logDens_Max = log10( GrackleTest_MassDensity_Max ); + GrackleTest_logDens_Range = GrackleTest_logDens_Max - GrackleTest_logDens_Min; + GrackleTest_logTemp_Min = log10( GrackleTest_TempOverMMW_Min * MOLECULAR_WEIGHT ); // from Temp/mu to built-in Temp in GAMER + GrackleTest_logTemp_Max = log10( GrackleTest_TempOverMMW_Max * MOLECULAR_WEIGHT ); + GrackleTest_logTemp_Range = GrackleTest_logTemp_Max - GrackleTest_logTemp_Min; + + +// (3) reset other general-purpose parameters +// --> a helper macro PRINT_RESET_PARA is defined in Macro.h + const long End_Step_Default = 10; // 10 * DT__GRACKLE_COOLING * cooling time + const double End_T_Default = 10.0*Const_Myr/UNIT_T; // 10 Myr + + if ( END_STEP < 0 ) { + END_STEP = End_Step_Default; + PRINT_RESET_PARA( END_STEP, FORMAT_LONG, "" ); + } + + if ( END_T < 0.0 ) { + END_T = End_T_Default; + PRINT_RESET_PARA( 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, " GrackleTest_DefaultTestMode = %d\n", GrackleTest_DefaultTestMode ); + Aux_Message( stdout, " GrackleTest_MassDensity_Min = %13.7e UNIT_D\n", GrackleTest_MassDensity_Min ); + Aux_Message( stdout, " = %13.7e g/cm^3\n", GrackleTest_MassDensity_Min*UNIT_D ); + Aux_Message( stdout, " = %13.7e mH/cm^3\n", GrackleTest_MassDensity_Min*UNIT_D/Const_mH ); + Aux_Message( stdout, " GrackleTest_MassDensity_Max = %13.7e UNIT_D\n", GrackleTest_MassDensity_Max ); + Aux_Message( stdout, " = %13.7e g/cm^3\n", GrackleTest_MassDensity_Max*UNIT_D ); + Aux_Message( stdout, " = %13.7e mH/cm^3\n", GrackleTest_MassDensity_Max*UNIT_D/Const_mH ); + Aux_Message( stdout, " GrackleTest_TempOverMMW_Min = %13.7e K\n", GrackleTest_TempOverMMW_Min ); + Aux_Message( stdout, " GrackleTest_TempOverMMW_Max = %13.7e K\n", GrackleTest_TempOverMMW_Max ); + Aux_Message( stdout, " GrackleTest_MFrac_Metal = %13.7e\n", GrackleTest_MFrac_Metal ); + Aux_Message( stdout, " GrackleTest_MFrac_e = %13.7e\n", GrackleTest_MFrac_e ); + Aux_Message( stdout, " GrackleTest_MFrac_HI = %13.7e\n", GrackleTest_MFrac_HI ); + Aux_Message( stdout, " GrackleTest_MFrac_HII = %13.7e\n", GrackleTest_MFrac_HII ); + Aux_Message( stdout, " GrackleTest_MFrac_HeI = %13.7e\n", GrackleTest_MFrac_HeI ); + Aux_Message( stdout, " GrackleTest_MFrac_HeII = %13.7e\n", GrackleTest_MFrac_HeII ); + Aux_Message( stdout, " GrackleTest_MFrac_HeIII = %13.7e\n", GrackleTest_MFrac_HeIII ); + Aux_Message( stdout, " GrackleTest_MFrac_HM = %13.7e\n", GrackleTest_MFrac_HM ); + Aux_Message( stdout, " GrackleTest_MFrac_H2I = %13.7e\n", GrackleTest_MFrac_H2I ); + Aux_Message( stdout, " GrackleTest_MFrac_H2II = %13.7e\n", GrackleTest_MFrac_H2II ); + Aux_Message( stdout, " GrackleTest_MFrac_DI = %13.7e\n", GrackleTest_MFrac_DI ); + Aux_Message( stdout, " GrackleTest_MFrac_DII = %13.7e\n", GrackleTest_MFrac_DII ); + Aux_Message( stdout, " GrackleTest_MFrac_HDI = %13.7e\n", GrackleTest_MFrac_HDI ); + Aux_Message( stdout, " GrackleTest_HeatingRate = %13.7e erg cm^-3 s^-1 n_H^-1\n", GrackleTest_HeatingRate ); + Aux_Message( stdout, " GrackleTest_CoolingRate = %13.7e erg cm^-3 s^-1 n_H^-2\n", GrackleTest_CoolingRate ); + 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[] ) +{ + + +// check +# ifdef GAMER_DEBUG + if ( GRACKLE_PRIMORDIAL >= GRACKLE_PRI_CHE_NSPE6 && Idx_e == Idx_Undefined ) Aux_Error( ERROR_INFO, "Idx_e is undefined !!\n" ); + if ( GRACKLE_PRIMORDIAL >= GRACKLE_PRI_CHE_NSPE6 && Idx_HI == Idx_Undefined ) Aux_Error( ERROR_INFO, "Idx_HI is undefined !!\n" ); + if ( GRACKLE_PRIMORDIAL >= GRACKLE_PRI_CHE_NSPE6 && Idx_HII == Idx_Undefined ) Aux_Error( ERROR_INFO, "Idx_HII is undefined !!\n" ); + if ( GRACKLE_PRIMORDIAL >= GRACKLE_PRI_CHE_NSPE6 && Idx_HeI == Idx_Undefined ) Aux_Error( ERROR_INFO, "Idx_HeI is undefined !!\n" ); + if ( GRACKLE_PRIMORDIAL >= GRACKLE_PRI_CHE_NSPE6 && Idx_HeII == Idx_Undefined ) Aux_Error( ERROR_INFO, "Idx_HeII is undefined !!\n" ); + if ( GRACKLE_PRIMORDIAL >= GRACKLE_PRI_CHE_NSPE6 && Idx_HeIII == Idx_Undefined ) Aux_Error( ERROR_INFO, "Idx_HeIII is undefined !!\n" ); + if ( GRACKLE_PRIMORDIAL >= GRACKLE_PRI_CHE_NSPE9 && Idx_HM == Idx_Undefined ) Aux_Error( ERROR_INFO, "Idx_HM is undefined !!\n" ); + if ( GRACKLE_PRIMORDIAL >= GRACKLE_PRI_CHE_NSPE9 && Idx_H2I == Idx_Undefined ) Aux_Error( ERROR_INFO, "Idx_H2I is undefined !!\n" ); + if ( GRACKLE_PRIMORDIAL >= GRACKLE_PRI_CHE_NSPE9 && Idx_H2II == Idx_Undefined ) Aux_Error( ERROR_INFO, "Idx_H2II is undefined !!\n" ); + if ( GRACKLE_PRIMORDIAL >= GRACKLE_PRI_CHE_NSPE12 && Idx_DI == Idx_Undefined ) Aux_Error( ERROR_INFO, "Idx_DI is undefined !!\n" ); + if ( GRACKLE_PRIMORDIAL >= GRACKLE_PRI_CHE_NSPE12 && Idx_DII == Idx_Undefined ) Aux_Error( ERROR_INFO, "Idx_DII is undefined !!\n" ); + if ( GRACKLE_PRIMORDIAL >= GRACKLE_PRI_CHE_NSPE12 && Idx_HDI == Idx_Undefined ) Aux_Error( ERROR_INFO, "Idx_HDI is undefined !!\n" ); + if ( GRACKLE_METAL && Idx_Metal == Idx_Undefined ) Aux_Error( ERROR_INFO, "Idx_Metal is undefined !!\n" ); + if ( EoS_DensTemp2Pres_CPUPtr == NULL ) Aux_Error( ERROR_INFO, "EoS_DensTemp2Pres_CPUPtr == NULL !!\n" ); +# endif + +// compute the gas log density by linear interpolation in x-direction + const double logDens = GrackleTest_logDens_Min + + GrackleTest_logDens_Range*( x / amr->BoxSize[0] ); + const double Dens = pow( 10, logDens ); + +// compute the gas log temperature by linear interpolation in y-direction + const double logTemp = GrackleTest_logTemp_Min + + GrackleTest_logTemp_Range*( y / amr->BoxSize[1] ); + const double Temp = pow( 10, logTemp ); + +// compute the gas pressure + const double Pres = EoS_DensTemp2Pres_CPUPtr( Dens, Temp, NULL, EoS_AuxArray_Flt, EoS_AuxArray_Int, + h_EoS_Table ); // assuming EoS requires no passive scalars + +// assume no momentum + const double MomX = 0.0; + const double MomY = 0.0; + const double MomZ = 0.0; + +// compute the total gas energy + const double Eint = EoS_DensPres2Eint_CPUPtr( Dens, Pres, NULL, EoS_AuxArray_Flt, + EoS_AuxArray_Int, h_EoS_Table ); // assuming EoS requires no passive scalars + const double 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; + + +// set passive scalars +// 6-species network + if ( GRACKLE_PRIMORDIAL >= GRACKLE_PRI_CHE_NSPE6 ) + { + fluid[Idx_e ] = Dens * GrackleTest_MFrac_e; + fluid[Idx_HI ] = Dens * GrackleTest_MFrac_HI; + fluid[Idx_HII ] = Dens * GrackleTest_MFrac_HII; + fluid[Idx_HeI ] = Dens * GrackleTest_MFrac_HeI; + fluid[Idx_HeII ] = Dens * GrackleTest_MFrac_HeII; + fluid[Idx_HeIII] = Dens * GrackleTest_MFrac_HeIII; + } + +// 9-species network + if ( GRACKLE_PRIMORDIAL >= GRACKLE_PRI_CHE_NSPE9 ) + { + fluid[Idx_HM ] = Dens * GrackleTest_MFrac_HM; + fluid[Idx_H2I ] = Dens * GrackleTest_MFrac_H2I; + fluid[Idx_H2II ] = Dens * GrackleTest_MFrac_H2II; + } + +// 12-species network + if ( GRACKLE_PRIMORDIAL >= GRACKLE_PRI_CHE_NSPE12 ) + { + fluid[Idx_DI ] = Dens * GrackleTest_MFrac_DI; + fluid[Idx_DII ] = Dens * GrackleTest_MFrac_DII; + fluid[Idx_HDI ] = Dens * GrackleTest_MFrac_HDI; + } + +// metallicity for metal cooling + if ( GRACKLE_METAL ) + fluid[Idx_Metal] = Dens * GrackleTest_MFrac_Metal; + +} // FUNCTION : SetGridIC + + + +//------------------------------------------------------------------------------------------------------- +// Function : Grackle_vHeatingRate_GrackleTest +// Description : Function to set Grackle's volumetric heating rate for GrackleTest +// +// Note : 1. Invoked by Grackle_Prepare() using the function pointer +// "Grackle_vHeatingRate_User_Ptr", which must be set by a test problem initializer +// 2. This function will be invoked by multiple OpenMP threads when OPENMP is enabled +// --> Please ensure that everything here is thread-safe +// 3. Returned rate should be in the unit of erg s^-1 cm^-3 +// +// Parameter : x/y/z : Target physical coordinates +// Time : Target physical time +// n_H : Hydrogen number density, should be in the unit of cm^-3 +// +// Return : volumetric_heating_rate +//------------------------------------------------------------------------------------------------------- +real_che Grackle_vHeatingRate_GrackleTest( const double x, const double y, const double z, const double Time, const double n_H ) +{ + const double volumetric_heating_rate_0 = n_H * GrackleTest_HeatingRate; // GrackleTest_HeatingRate is in the unit of erg cm^-3 s^-1 n_H^-1 + const double volumetric_cooling_rate_0 = n_H * n_H * GrackleTest_CoolingRate; // GrackleTest_CoolingRate is in the unit of erg cm^-3 s^-1 n_H^-2 + +// assume uniform distribution if it is not specified + if ( GrackleTest_DefaultTestMode != 3 ) + return volumetric_heating_rate_0 - volumetric_cooling_rate_0; + +// an example of spatial distribution only for GrackleTest_DefaultTestMode == 3 +// two 2D Gaussian distributions: one heating and one cooling + const double Center_H[3] = { amr->BoxCenter[0]+0.25*amr->BoxSize[0], amr->BoxCenter[1], amr->BoxCenter[2] }; // center of heating + const double Center_C[3] = { amr->BoxCenter[0]-0.25*amr->BoxSize[0], amr->BoxCenter[1], amr->BoxCenter[2] }; // center of cooling + const double r_to_center_h = sqrt( SQR(x-Center_H[0]) + SQR(y-Center_H[1]) ); // radius to heating + const double r_to_center_c = sqrt( SQR(x-Center_C[0]) + SQR(y-Center_C[1]) ); // radius to cooling + const double distr_width = 0.0625*amr->BoxSize[0]; // width of Gaussian + const double distr_r_max = 0.1250*amr->BoxSize[0]; // cutoff radius + const real_che volumetric_heating_rate_r = ( r_to_center_h <= distr_r_max ) ? volumetric_heating_rate_0 * exp( -0.5*SQR(r_to_center_h/distr_width) ) + : (real_che)0.0; + const real_che volumetric_cooling_rate_r = ( r_to_center_c <= distr_r_max ) ? volumetric_cooling_rate_0 * exp( -0.5*SQR(r_to_center_c/distr_width) ) + : (real_che)0.0; + + return volumetric_heating_rate_r - volumetric_cooling_rate_r; + +} // FUNCTION : Grackle_vHeatingRate_GrackleTest +#endif // #if ( MODEL == HYDRO && defined SUPPORT_GRACKLE ) + + + +//------------------------------------------------------------------------------------------------------- +// Function : Init_TestProb_Hydro_GrackleTest +// Description : Test problem initializer +// +// Note : None +// +// Parameter : None +// +// Return : None +//------------------------------------------------------------------------------------------------------- +void Init_TestProb_Hydro_GrackleTest() +{ + + if ( MPI_Rank == 0 ) Aux_Message( stdout, "%s ...\n", __FUNCTION__ ); + + +// validate the compilation flags and runtime parameters + Validate(); + + +# if ( MODEL == HYDRO && defined SUPPORT_GRACKLE ) +// set the problem-specific runtime parameters + SetParameter(); + + +// set the function pointers of various problem-specific routines + Init_Function_User_Ptr = SetGridIC; + Grackle_vHeatingRate_User_Ptr = Grackle_vHeatingRate_GrackleTest; +# ifdef SUPPORT_HDF5 + Output_HDF5_InputTest_Ptr = LoadInputTestProb; +# endif +# endif // #if ( MODEL == HYDRO && defined SUPPORT_GRACKLE ) + + + if ( MPI_Rank == 0 ) Aux_Message( stdout, "%s ... done\n", __FUNCTION__ ); + +} // FUNCTION : Init_TestProb_Hydro_GrackleTest