Skip to content

Commit

Permalink
Add ELBDM test problems
Browse files Browse the repository at this point in the history
  • Loading branch information
ChunYen-Chen committed Dec 12, 2024
1 parent 4e74710 commit 12e84db
Show file tree
Hide file tree
Showing 47 changed files with 1,521 additions and 526 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@
[Ji-hoon Kim, et al., 2016, ApJ, 833, 202](https://dx.doi.org/10.3847/1538-4357/833/2/202) [(arXiv: 1610.03066)](https://dx.doi.org/10.3847/1538-4357/833/2/202)
2. Other references
- [AGORA website](https://sites.google.com/site/santacruzcomparisonproject/)
- [AGORA initial condition](http://goo.gl/8JzbIJ)
- [AGORA initial condition](https://goo.gl/8JzbIJ)
- [Enzo setup](https://bitbucket.org/enzo/enzo-dev/src/19f4a44e06f1c386573dc77b3608ba66b64d93bc/run/Hydro/Hydro-3D/AgoraGalaxy/?at=week-of-code)
- [Goldbaum et al. 2016](https://arxiv.org/abs/1605.00646)
- [yt hub](https://girder.hub.yt/#collection/5736481ddd9119000164acf1)
Expand Down
143 changes: 143 additions & 0 deletions doc/wiki/Test-Problem-related/Test-Problems:-DiskHeating.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,143 @@
> [!CAUTION]
> Please do not edit this file(page) manually since the workflow will overwrite your changes.
>
> This file(page) is automatically generated by the workflow `Update test problem wiki page` using the script `tool/wiki/sync_test_problem_pages.py`.
>
> The workflow is triggered by push changes to any of `example/test_problem/*/*/README.md` and `tool/wiki/sync_test_problem_pages.py`.

# Compilation flags
- Must enable
- [[MODEL=ELBDM | Installation: Simulation-Options#MODEL]]
- [[GRAVITY | Installation: Simulation-Options#GRAVITY]]
- [[PARTICLE | Installation: Simulation-Options#PARTICLE]]
- [[STORE_PAR_ACC | Installation: Simulation-Options#STORE_PAR_ACC]]
- [[SUPPORT_HDF5 | Installation: Simulation-Options#SUPPORT_HDF5]]
- [[SUPPORT_GSL | Installation: Simulation-Options#SUPPORT_GSL]] (optional, only useful for thin disk)
- Must disable
- [[COMOVING | Installation: Simulation-Options#COMOVING]]
- Available options
- [[Miscellaneous Options | Installation: Simulation-Options#miscellaneous-options]]


# Quick start
0. Generate `gamer`
1. Copy `generate_make.sh` to the directory `src`
2. Generate `Makefile`
```bash
sh generate_make.sh
```
3. Compile `gamer`
```bash
make clean
make -j 4
```

1. Download initial conditions of `m_22=0.4`, `M_h=7e10` Msun halo and stellar disk with
```bash
sh download_ic.sh
```

2. Ensure [[PAR_INIT | #PAR_INIT]] = 1 and [[OPT__INIT]] = 3 in `Input__Parameter`

3. Default [[END_T | ]] is 2.5e-1 (about 3.5 Gyr) as in [Yang et al. 2023](https://doi.org/10.1093/mnras/stae793) and [[OUTPUT_DT | ]] is 1.0e-2 (about 0.14 Gyr)

4. To switch to a high-resolution run, command `ln -sf ic_files/PAR_IC_0.4_M7 PAR_IC`
Set [[PAR_NPAR | ]]=80000000, [[MAX_LEVEL | Runtime-Parameters:-Refinement#MAX_LEVEL]]=3, and change all values in `Input__Flag_NParPatch` to `800`


# General initial condition setup
1. Disk

a. Generate the disk via modified [GALIC](https://github.com/HsunYeong/GALIC.git)
The snapshots have the filenames `snap_XXX.hdf5`

b. Set the filename, units in `get_par_ic.py` to match the GALIC set-up

c. Set center to be location of the soliton in `get_par_ic.py`

d. Execute `get_par_ic.py`, it will generate `PAR_IC`

2. Halo

a. If the data is binary file `UM_IC`

* Set [[OPT__INIT | ]] = 3 and [[PAR_INIT | ]] = 1
* `Input__UM_IC_RefineRegion` is required

b. If the data is GAMER snapshot

* Command `ln -s Data_XXXXXX RESTART` to create a soft link
* Set [[OPT__INIT | ]] = 2 and [[PAR_INIT | ]] = 2 in `Input__Parameter`
* Turn on `AddParWhenRestart` and `AddParWhenRestartByFile` in `Input__TestProb`
* Set `AddParWhenRestartNPar` in `Input__TestProb`
* Turn on [[OPT__RESTART_RESET | ]] in `Input__Parameter`
# Recommand to turn off `AddParWhenRestart`, `AddParWhenRestartByFile`, [[OPT__RESTART_RESET | ]] right after the simulation starts

c. The code for [FDM halo reconstruction](https://github.com/calab-ntu/psidm-halo-reconstruction)

3. Thin disk (optional)

a. Create a soft link for restart
```bash
ln -s Data_XXXXXX RESTART
```

b. Set [[OPT__INIT | ]] = 2 and [[PAR_INIT | ]] = 2 and turn on [[OPT__RESTART_RESET | ]] in Input__Parameter

c. Turn on `AddParWhenRestart` in `Input__TestProb`

d. Set `AddParWhenRestartNPar`, `Disk_Mass`, `Disk_R`, and `DispTableFile` in `Input__TestProb`


# Analysis scripts
1. `particle_proj.py`, `plot_halo_slice.py`

* Plot the projection of disk particles or halo density slice
* Output files: `particle_proj_*.png`, `Data_*_Slice_x_Dens.png`

2. `plot_halo_density.py`, `plot_halo_potential.py`

* Compute and plot shell-averaged halo density or gravitational potential profiles
* Output files: `Data_*_1d-Profile_radius_Dens.png`, `Halo_Dens_Data_*.npy`,
`Data_*_1d-Profile_radius_Pote.png`, `Halo_Pote_Data_*.npy`

3. `data_disk.py`

* Compute the disk information (rotation speed, velocity dispersion, surface density, scale height, etc.)
* Output files: `Data_Disk_*.npy`

4. `data_halo.py`

* Compute the halo information (enclosed mass, velocity dispersion, etc.)
* Required files: `Halo_Dens_*.npy` (generated by `plot_halo_density.py`) and `Halo_Pote_*.npy` (by `plot_halo_potential.py`)
* Output files: `Data_Halo_*.npy`

5. `vel_distribution.py`

* Compute the velocity distribution in 2-kpc-wide radial bins centered on R = 4, 6, 8, 10 kpc
* Output files: `Vel_data_*.npz`, `vel_*.png`

6. `get_heating.py`

* Compute the theoretical heating rate of the stellar disk as a function of radius
* Required files: `Data_Disk_*.npy` (generated by `data_disk.py`), `Data_Halo_*.npy` (by `data_halo.py`)
* Output files: `Heating_*.npz`

7. `get_heating_rate_theory.py`

* Get the time-averaged theoretical heating rate
* Required files: `Heating_*.npz` (generated by `get_heating.py`)
* Output files: printed plain text

8. `get_heating_rate_simulation.py`

* Compute ensemble- and time-averaged disk heating rates measured from simulation data
* Required files: `Vel_data_*.npz` (generated by `vel_distribution.py`)
* Output files: `sigma_z_sqr.png` and printed plain text

9. `plot_data_example.py`

* An example script to plot the angle-averaged disk rotation curve and shell-averaged density profile
* Required files: `Data_Disk_*.npy` (generated by `data_disk.py`) and/or `Data_Halo_*.npy` (by `data_halo.py`)
* Output files: `Rotation_Curve.png`, `Halo_Density_Profile.png`
26 changes: 26 additions & 0 deletions doc/wiki/Test-Problem-related/Test-Problems:-GaussianWavePacket.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,26 @@
> [!CAUTION]
> Please do not edit this file(page) manually since the workflow will overwrite your changes.
>
> This file(page) is automatically generated by the workflow `Update test problem wiki page` using the script `tool/wiki/sync_test_problem_pages.py`.
>
> The workflow is triggered by push changes to any of `example/test_problem/*/*/README.md` and `tool/wiki/sync_test_problem_pages.py`.

# Compilation flags
- Must enable
- [[MODEL=ELBDM | Installation: Simulation-Options#MODEL]]
- Must disable
- [[GRAVITY | Installation: Simulation-Options#GRAVITY]]
- [[PARTICLE | Installation: Simulation-Options#PARTICLE]]
- Available options
- [[Miscellaneous Options | Installation: Simulation-Options#miscellaneous-options]]


# Default setup
1. Evolve the Gaussian wave packet for half of box
2. Apply the analytical solution as user-defined BC
--> Set [[OPT__BC_FLU_* | ]] = 4


# Note
1. Only support 1D --> Use `Gau_XYZ` to control the propagation direction
90 changes: 90 additions & 0 deletions doc/wiki/Test-Problem-related/Test-Problems:-HaloMerger.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,90 @@
> [!CAUTION]
> Please do not edit this file(page) manually since the workflow will overwrite your changes.
>
> This file(page) is automatically generated by the workflow `Update test problem wiki page` using the script `tool/wiki/sync_test_problem_pages.py`.
>
> The workflow is triggered by push changes to any of `example/test_problem/*/*/README.md` and `tool/wiki/sync_test_problem_pages.py`.

Compilation flags:
========================================
Enable : MODEL=ELBDM, GRAVITY (, PARTICLE, SUPPORT_GSL)
Disable: COMOVING


# Default setup
1. Code units
(1) UNIT_L = Mpc/h, where h=0.6955 is the present dimensionless Hubble parameter
= 1437.814521231748 kpc
(2) UNIT_V = 100 km/s
(3) UNIT_D = rho_bg (background matter density at z=0)
= 38.06 Msun/kpc^3
--> Mass density and wavefunction are normalized to rho_bg
(4) UNIT_T = 14068.4678922741 Myr
(5) UNIT_M = 1.131e+11 Msun

2. ELBDM_MASS 1.0e-22
ELBDM_REMOVE_MOTION_CM 0
ELBDM_TAYLOR3_AUTO 0

3. OPT__BC_FLU_* 1 (periodic)
OPT__BC_POT 2 (isolated)

4. MAX_LEVEL 3
OPT__FLAG_RHO 1

5. END_T 0.25

# Note
1. Simulate the merger of halos and solitons in an external potential

2. Edit `Input_TestProb_Halo` to specify the parameters for the halos
a. Add new parameters with increasing indexes if the number of halos is more than 2
b. For `HaloMerger_Halo_InitMode` == 1
- The initial condition of halos is constructed by reading the `HALO_IC` file as a table and performing linear interpolation
- Note that this `HALO_IC` has nothing to do with the built-in option [[OPT__INIT | ]] = 3
- The `HALO_IC` must be single AMR level
- If there is a halo `UM_IC` (for [[OPT__INIT | ]] = 3) that has multiple AMR levels (and there is `Input__UM_IC_RefineRegion`),
the Python script `Make_UM_IC_uniform.py` is provided to convert it to
a single-level `UM_IC` with the specified level, which can be used as the `HALO_IC`
(However, those levels higher than Target_lv+log_2( 2*PatchSize ) cannot be handled and will be ignored)
- If converting the halo UM_IC to be single-level is not feasible,
switching to the initialization option [[OPT__INIT | ]] = 3 to load a multi-level `UM_IC` is also an alternative
(However, only one `UM_IC` can be loaded and all the parameters in `Input_TestProb_Halo` will not be used)
- Note that `HaloMerger_Halo_*_CenCoord*` in `Input__TestProb_Halo` is the `HALO_IC` box center rather than the exact halo center

3. Edit `Input_TestProb_Soliton` to specify the parameters for the solitons
a. Add new parameters with increasing indexes if the number of solitons is more than 2
b. For `HaloMerger_Soliton_InitMode` == 1
- The initial condition of solitons is constructed by reading the table of soliton density profile
(the density profile will be rescaled to the given core radius or core density if `HaloMerger_Soliton_*_DensProf_Rescale` == 1)
c. For `HaloMerger_Soliton_InitMode` == 2
- The initial condition of solitons is constructed by using the analytical formula of soliton density profile

4. Edit `Input_TestProb_ParCloud` to specify the parameters for the particle clouds
a. Add new parameters with increasing indexes if the number of particle clouds is more than 2
b. For `HaloMerger_ParCloud_InitMode` == 1
- The initial condition of particle clouds is constructed by reading the table of density profile and using `Par_EquilibriumIC()`
- The default particle clouds use the HaloDensityProfile (see Note 7. below) to represent the CDM halos
- Enable compilation options: [[PARTICLE | ]], [[SUPPORT_GSL | ]]
- Set the parameter [[OPT__FREEZE_FLUID | ]] to 1 in `Input__Parameter` for particle-only simulations

5. Turn on [[OPT__EXT_POT | ]] == 1 to add external potential
a. The external potential of a uniform-density sphere, which is proportional to r^2 inside and proportional to 1/r outside

6. Download the default `HALO_IC` file: `sh download_ic.sh`
a. Halo mass = 4.0960e+09 Msun
b. Without soliton
c. N = 640, L = 0.0646921095 Mpc/h, single-precision

7. Generate the default HaloDensityProfile and SolitonDensityProfile: `python Make_DensityProfile.py`
a. The parameters for the halo density profile are used to fit the ELBDM halo

8. Some examples of yt visualization scripts are put in `plot_script`

9. The corresponding wavelength is 0.00083778702 Mpc/h when ELBDM_MASS = 1.0e-22 and velocity = 1.0*100 km/s
-> Make sure the resolution is high enough to resolve the wavelength

10. The input halos and solitons may overlap with each other initially
-> Their wavefunctions, instead of densities, are added directly
-> Note that the interference will cause the density distribution to be different from the sum of individual density
105 changes: 105 additions & 0 deletions doc/wiki/Test-Problem-related/Test-Problems:-IsolatedHalo.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,105 @@
> [!CAUTION]
> Please do not edit this file(page) manually since the workflow will overwrite your changes.
>
> This file(page) is automatically generated by the workflow `Update test problem wiki page` using the script `tool/wiki/sync_test_problem_pages.py`.
>
> The workflow is triggered by push changes to any of `example/test_problem/*/*/README.md` and `tool/wiki/sync_test_problem_pages.py`.

# Compilation flags
- Must enable
- [[MODEL=ELBDM | Installation: Simulation-Options#MODEL]]
- [[GRAVITY | Installation: Simulation-Options#GRAVITY]]
- Must disable
- [[COMOVING | Installation: Simulation-Options#COMOVING]]
- Available options
- [[Miscellaneous Options | Installation: Simulation-Options#miscellaneous-options]]


# Default setup
1. Code units
(1) UNIT_L = Mpc/h, where h=0.6955 is the present dimensionless Hubble parameter
(2) UNIT_V = 100 km/s
(3) UNIT_D = rho_bg (background matter density at z=0)
--> Mass density and wavefunction are normalized to rho_bg

2. ELBDM
| Parameter name | value |
|--- |--- |
| ELBDM_MASS | 8.0e-23 |
| ELBDM_REMOVE_MOTION_CM | 1 |
| ELBDM_TAYLOR3_AUTO | 0 |

3. Boundary conditions
| Parameter name | value |
|--- |--- |
| OPT__BC_FLU_* | 1 |
| OPT__BC_POT | 2 |

4. AMR
| Parameter name | value |
|--- |--- |
| MAX_LEVEL | 0 |


# Libyt covering_grid setup
1. Code units
(1) UNIT_L = Mpc/h, where h=0.6955 is the present dimensionless Hubble parameter
(2) UNIT_V = 100 km/s
(3) UNIT_D = rho_bg (background matter density at z=0)
--> Mass density and wavefunction are normalized to rho_bg

2. ELBDM
| Parameter name | value |
|--- |--- |
| ELBDM_MASS | 8.0e-23 |
| ELBDM_REMOVE_MOTION_CM | 1 |
| ELBDM_TAYLOR3_AUTO | 0 |

3. Boundary conditions
| Parameter name | value |
|--- |--- |
| OPT__BC_FLU_* | 1 |
| OPT__BC_POT | 2 |

4. AMR
| Parameter name | value |
|--- |--- |
| MAX_LEVEL | 0 |
| OPT__FLAG_RHO | 1 |

5. End simulation condition
| Parameter name | value |
|--- |--- |
| END_T | 5.7116620e-03 |
| END_STEP | 32 |

6. libyt
| Parameter name | value |
|--- |--- |
| YT_SCRIPT | libyt_script/inline_script_covering_grid |
| YT_VERBOSE | 1 |


# Note
1. Download the IC file
```bash
sh download_ic.sh
```

2. Some examples of yt visualization scripts are put in `yt_script`

3. Simulate a single isolated halo extracted from a cosmological simulation

4. About libyt covering_grid test
a. Use submit script `./libyt_script/submit_gamer.job` for job submission

b. Put `./libyt_script/inline_script_covering_grid.py` under the same folder as gamer

c. For determining the `left_edge` and `dims` in function `ds.covering_grid`:
`left_edge`: LV1 resolution is 0.175/512/2 ; region covered by LV1 box (by `ds.covering_grid`) is 0.175/512/2*512; 0.04375 = (0.175 - 0.175/512/2*512)/2
`dims`: Plan to cover region with half of the simulation box length, i.e. will have 256X256X256 base level cells -> refine to MAX_LEVEL=1 -> LV1 cells is 512X512X512

d. Use `plot_slice-dens_covering_grid.py` to extract density slices from `.npz` files

e. Use `make_movie.sh` to convert `.png` pictures into `.mp4`
Loading

0 comments on commit 12e84db

Please sign in to comment.