Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Update the CDM halo construction in HaloMerger test problem #134

Open
wants to merge 4 commits into
base: psidm
Choose a base branch
from
Open

Update the CDM halo construction in HaloMerger test problem #134

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
f56cc3c
Use the FDM density profile to construct CDM halo
hsinhaoHHuang Sep 28, 2024
c14daab
Update README
hsinhaoHHuang Sep 28, 2024
c5ee852
Merge branch 'psidm' into psidm-halo-merger_CDM
hsinhaoHHuang Oct 30, 2024
da0b8aa
Support constructing CDM halos with external potential in HaloMerger
hsinhaoHHuang Oct 31, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
82 changes: 82 additions & 0 deletions example/test_problem/ELBDM/HaloMerger/FDMHaloDensityProfile
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,82 @@
# r (code_length) density (code_density)
5.90000000e-04 1.13660000e+05
6.01141644e-04 1.13656279e+05
6.32624146e-04 1.06501823e+05
6.65755423e-04 9.82635653e+04
7.00621825e-04 9.05168097e+04
7.37314222e-04 8.34638603e+04
7.75928243e-04 7.64868480e+04
8.16564526e-04 6.96644797e+04
8.59328980e-04 6.35464195e+04
9.04333059e-04 5.82486348e+04
9.51694055e-04 5.33252240e+04
1.00153540e-03 4.93928376e+04
1.05398700e-03 4.60739413e+04
1.10918555e-03 4.35503984e+04
1.16727492e-03 4.15915695e+04
1.22840649e-03 4.01334779e+04
1.29273959e-03 3.94419839e+04
1.36044190e-03 3.87196510e+04
1.43168986e-03 3.79971660e+04
1.50666916e-03 3.76434488e+04
1.58557521e-03 3.67997220e+04
1.66861367e-03 3.60365889e+04
1.75600094e-03 3.51860751e+04
1.84796480e-03 3.41940573e+04
1.94474491e-03 3.30525659e+04
2.04659351e-03 3.19651276e+04
2.15377604e-03 3.07071955e+04
2.26657185e-03 2.93033596e+04
2.38527491e-03 2.78406016e+04
2.51019459e-03 2.62849468e+04
2.64165646e-03 2.50520799e+04
2.78000314e-03 2.43616745e+04
2.92559521e-03 2.43881516e+04
3.07881210e-03 2.47977486e+04
3.24005314e-03 2.49677406e+04
3.40973857e-03 2.40270939e+04
3.58831063e-03 2.17110865e+04
3.77623472e-03 1.84180138e+04
3.97400061e-03 1.52040351e+04
4.18212374e-03 1.28395146e+04
4.40114653e-03 1.14651875e+04
4.63163980e-03 1.06861057e+04
4.87420427e-03 1.01002103e+04
5.12947214e-03 9.50788131e+03
5.39810868e-03 8.73582821e+03
5.68081404e-03 7.75677212e+03
5.97832502e-03 6.78824550e+03
6.29141700e-03 6.03614732e+03
6.62090598e-03 5.41074813e+03
6.96765068e-03 4.76378247e+03
7.33255482e-03 4.08847618e+03
7.71656942e-03 3.46692604e+03
8.12069532e-03 2.98810224e+03
8.54598578e-03 2.62300084e+03
8.99354920e-03 2.27291464e+03
9.46455205e-03 1.93016437e+03
9.96022187e-03 1.63331875e+03
1.04818505e-02 1.39506325e+03
1.10307975e-02 1.17154780e+03
1.16084934e-02 9.62725766e+02
1.22164440e-02 7.93051154e+02
1.28562337e-02 6.62872300e+02
1.35295299e-02 5.72070436e+02
1.42380875e-02 4.79246216e+02
1.49837530e-02 3.87451856e+02
1.57684699e-02 3.08119262e+02
1.65942834e-02 2.39956006e+02
1.74633458e-02 1.92108585e+02
1.83779219e-02 1.54373954e+02
1.93403955e-02 1.15859677e+02
2.03532750e-02 8.50437637e+01
2.14192001e-02 6.41426070e+01
2.25409491e-02 4.62022768e+01
2.37214453e-02 3.16797249e+01
2.49637656e-02 1.97108242e+01
2.62711476e-02 1.02333643e+01
2.76469988e-02 4.05611318e+00
2.90949050e-02 1.10750601e+00
3.06186397e-02 1.94123040e-01
3.22221742e-02 2.11390042e-02
3.39096877e-02 1.38335859e-03
0 ...oblem/ELBDM/HaloMerger/HaloDensityProfile → ...em/ELBDM/HaloMerger/FitHaloDensityProfile
View file
Open in desktop
File renamed without changes.
8 changes: 6 additions & 2 deletions example/test_problem/ELBDM/HaloMerger/Input__TestProb_ParCloud
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,18 +6,22 @@ HaloMerger_ParCloud_1_CenCoordZ 0.12938422
HaloMerger_ParCloud_1_VelocityX 0.5 # x/y/z-component of the bulk velocity of the 1st particle cloud [0.0]
HaloMerger_ParCloud_1_VelocityY 0.0
HaloMerger_ParCloud_1_VelocityZ 0.0
HaloMerger_ParCloud_1_DensProf_Filename HaloDensityProfile # filename of the density profile table for the 1st particle cloud (HaloMerger_ParCloud_InitMode == 1 only)
HaloMerger_ParCloud_1_DensProf_Filename FDMHaloDensityProfile # filename of the density profile table for the 1st particle cloud (HaloMerger_ParCloud_InitMode == 1 only)
HaloMerger_ParCloud_1_DensProf_MaxR 0.03234605475 # maximum radius for particles for the 1st particle cloud (must > 0.0) (HaloMerger_ParCloud_InitMode == 1 only) [0.5*amr->BoxSize[0]]
HaloMerger_ParCloud_1_RSeed 123 # random seed for setting particle position and velocity for the 1st particle cloud (must >= 0) (HaloMerger_ParCloud_InitMode == 1 only) [123]
HaloMerger_ParCloud_1_NPar 5000000 # number of particles for the 1st particle cloud (must >= 0) (HaloMerger_ParCloud_InitMode == 1 only) [0]
HaloMerger_ParCloud_1_BuiltWithExtPot 0 # whether to consider an external potential for the velocity dispersion for the 1st particle cloud (HaloMerger_ParCloud_InitMode == 1 only) [0]
HaloMerger_ParCloud_1_ExtPot_Filename ParCloud_ExtPotTable # filename of the external potential table for the 1st particle cloud (HaloMerger_ParCloud_InitMode == 1 only) (HaloMerger_ParCloud_1_BuiltWithExtPot == 1 only)
# 2nd particle cloud
HaloMerger_ParCloud_2_CenCoordX 0.19407633
HaloMerger_ParCloud_2_CenCoordY 0.12938422
HaloMerger_ParCloud_2_CenCoordZ 0.12938422
HaloMerger_ParCloud_2_VelocityX -0.5
HaloMerger_ParCloud_2_VelocityY 0.0
HaloMerger_ParCloud_2_VelocityZ 0.0
HaloMerger_ParCloud_2_DensProf_Filename HaloDensityProfile
HaloMerger_ParCloud_2_DensProf_Filename FDMHaloDensityProfile
HaloMerger_ParCloud_2_DensProf_MaxR 0.03234605475
HaloMerger_ParCloud_2_RSeed 123
HaloMerger_ParCloud_2_NPar 5000000
HaloMerger_ParCloud_2_BuiltWithExtPot 0
HaloMerger_ParCloud_2_ExtPot_Filename ParCloud_ExtPotTable
4 changes: 2 additions & 2 deletions example/test_problem/ELBDM/HaloMerger/Make_DensityProfile.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -99,7 +99,7 @@ def Soliton_fitting_analytical_dens(r, m22, r_c):

###################################################################################################################
# save to file
np.savetxt( 'HaloDensityProfile',
np.savetxt( 'FitHaloDensityProfile',
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How about printing a warning message if the file already exists?

np.column_stack( (halo_densprof_radius, halo_densprof_density) ),
fmt=' %9.8e',
header=' r density' )
Expand Down Expand Up @@ -148,7 +148,7 @@ def Soliton_fitting_analytical_dens(r, m22, r_c):

# save the figure
fig.subplots_adjust( top=0.93, bottom=0.1, left=0.1, right=0.97 )
fig.savefig( 'fig_HaloDensityProfile.png' )
fig.savefig( 'fig_FitHaloDensityProfile.png' )
plt.close()

fig = plt.figure()
Expand Down
94 changes: 94 additions & 0 deletions example/test_problem/ELBDM/HaloMerger/Make_ParCloud_ExtPotTable.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,94 @@
import numpy as np
import matplotlib.pyplot as plt


###################################################################################################################
# Note that the external potential table has to have the same radial bins as the density profile of ParCloud
ParCloud_DensProf_Filename = 'FDMHaloDensityProfile'
###################################################################################################################


###################################################################################################################
# code units (in cgs)
UNIT_L = 4.4366320345075490e+24
UNIT_D = 2.5758579724476994e-30
UNIT_T = 4.4366320345075490e+17

NEWTON_G = 6.6738e-8 / ( 1.0/UNIT_D/(UNIT_T**2) )
###################################################################################################################


###################################################################################################################
def Potential_UniDenSph(r, M, R):
return -NEWTON_G*M/r if r >= R else NEWTON_G*M/(2.0*R*R*R)*(r*r-3.0*R*R)
###################################################################################################################


###################################################################################################################
# parameters for the external potential
# the default is the potential of an uniform density sphere with 3x mass and 0.3x radius of the default CDM halo
ParCloud_ExtPot_UniDenSph_M = 3.0*3.618847000e-02 # in code_mass
ParCloud_ExtPot_UniDenSph_R = 0.3*3.234605475e-02 # in code_length
###################################################################################################################


###################################################################################################################
# output the information
print( 'Information' )
print( 'ParCloud_DensProf_Filename = '+ParCloud_DensProf_Filename )
print( 'UNIT_L = {: >16.8e} cm'.format( UNIT_L ) )
print( 'UNIT_D = {: >16.8e} g/cm**3'.format( UNIT_D ) )
print( 'UNIT_T = {: >16.8e} s'.format( UNIT_T ) )
print( 'ParCloud_ExtPot_UniDenSph_M = {: >16.8e} UNIT_M'.format( ParCloud_ExtPot_UniDenSph_M ) )
print( 'ParCloud_ExtPot_UniDenSph_R = {: >16.8e} UNIT_L'.format( ParCloud_ExtPot_UniDenSph_R ) )
###################################################################################################################


###################################################################################################################
# create the external potential table
ParCloud_DensProf_radius, _ = np.loadtxt( ParCloud_DensProf_Filename, skiprows=1, unpack=True )
ParCloud_ExtPot_table_radius = np.append( ParCloud_DensProf_radius, 2.0*ParCloud_DensProf_radius[-1] ) # probably a bug, we need to add one extra row, which will not be used
ParCloud_ExtPot_table_potential = np.array([Potential_UniDenSph( radius, ParCloud_ExtPot_UniDenSph_M, ParCloud_ExtPot_UniDenSph_R ) for radius in ParCloud_ExtPot_table_radius ])
###################################################################################################################


###################################################################################################################
# save to file
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How about printing a warning message if the file already exists?

np.savetxt( 'ParCloud_ExtPotTable',
np.column_stack( (ParCloud_ExtPot_table_radius, ParCloud_ExtPot_table_potential) ),
fmt=' %9.8e',
header=' r potential' )
Comment on lines +57 to +60
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How about using this format?

Suggested change
np.savetxt( 'ParCloud_ExtPotTable',
np.column_stack( (ParCloud_ExtPot_table_radius, ParCloud_ExtPot_table_potential) ),
fmt=' %9.8e',
header=' r potential' )
np.savetxt( 'ParCloud_ExtPotTable',
np.column_stack( (ParCloud_ExtPot_table_radius, ParCloud_ExtPot_table_potential) ),
fmt='%23.8e',
header='%21s %23s'%('r', 'potential') )

This can make it easier to maintain the text file.

###################################################################################################################


###################################################################################################################
# plot to images
fig = plt.figure()
ax = fig.add_subplot(111)

# plot some important values for reference
ax.axvline( ParCloud_ExtPot_UniDenSph_R, linestyle='--', color='grey', label=r'$R$' )

# plot the external potential
ax.plot( ParCloud_ExtPot_table_radius, ParCloud_ExtPot_table_potential, '-', color='r', label=r'$\Phi(r)$' )

# annotate the information
ax.annotate( r'$UniDenSph_M$ = %.8e'%ParCloud_ExtPot_UniDenSph_M+'\n'+
r'$UniDenSph_R$ = %.8e'%ParCloud_ExtPot_UniDenSph_R,
xy=(0.02,0.80), xycoords='axes fraction' )

# setting for the figure
ax.set_xscale('log')
ax.set_xlim( 0.5*np.min(ParCloud_ExtPot_table_radius), 2.0*np.max(ParCloud_ExtPot_table_radius) )

# set the labels
ax.set_xlabel( r'$r$'+' (code_length)' )
ax.set_ylabel( r'$\Phi$'+' (code_length$^2$/code_time$^2$)' )
fig.suptitle( 'External Potential for ParCloud' )
ax.legend( loc='lower right' )

# save the figure
fig.subplots_adjust( top=0.93, bottom=0.1, left=0.15, right=0.97 )
fig.savefig( 'fig_ParCloud_ExtPotTable.png' )
plt.close()
###################################################################################################################
25 changes: 18 additions & 7 deletions example/test_problem/ELBDM/HaloMerger/README
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,9 +59,12 @@ Note:
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
- The default particle clouds use the FDMHaloDensityProfile (see Note 7. below) to represent the CDM halos
- The initial condition can be constructed as an equilibrium with external potential by reading the provided external potential table
- The external potential table "ParCloud_ExtPotTable" can be generated by running the Python script "python3 Make_ParCloud_ExtPotTable.py",
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If the script is for Python3 only, how about adding a version check in the script?

where the default is the potential of a uniform density sphere with 3x mass and 0.3x radius of the default CDM halo
c. Enable compilation options: PARTICLE, SUPPORT_GSL
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since we do not edit Makefile directly anymore, how about changing to add options in generage_make.sh?

d. 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
Expand All @@ -71,14 +74,22 @@ Note:
b. Without soliton
c. N = 640, L = 0.0646921095 Mpc/h, single-precision

7. Generate the default HaloDensityProfile and SolitonDensityProfile: python Make_DensityProfile.py
7. The default FDMHaloDensityProfile is extracted from the default HALO_IC
a. Run the default ELBDM case to get Data_000000
b. Adjust the plot_script/plot__density_profile.py: set `r_sphere = (50.0, 'kpc')` and `nbin = 128`
c. Run `python3 plot__density_profile.py -s 0 -e 0`
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If the script is for Python3 only, how about adding a version check in the script?

b. Remove the inner region (r < 6e-4 code_length) of the output Data_000000_density_profile to achieve a more consistent total enclosed mass
-> In this way, the CDM halo would look more similar and be easier to compare to the ELBDM case

8. Generate the default FitHaloDensityProfile and SolitonDensityProfile: python Make_DensityProfile.py
a. The parameters for the halo density profile are used to fit the ELBDM halo
-> the total mass of the CDM halo would be similar to the ELBDM case, but the distribution is more extended at a large radius

8. Some examples of yt visualization scripts are put in "plot_script"
9. 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
10. 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
11. 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
2 changes: 1 addition & 1 deletion example/test_problem/ELBDM/HaloMerger/plot_script/plot__density_profile.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@
dpi = 150
nbin = 32

Ref_DensProf_filename = '../HaloDensityProfile'
Ref_DensProf_filename = '../FDMHaloDensityProfile'

yt.enable_parallelism()
ts = yt.DatasetSeries([ prefix+'Data_%06d'%idx for idx in range(idx_start, idx_end+1, didx) ])
Expand Down
Loading