diff --git a/example/test_problem/ELBDM/HaloMerger/FDMHaloDensityProfile b/example/test_problem/ELBDM/HaloMerger/FDMHaloDensityProfile new file mode 100644 index 0000000000..283c2b04ca --- /dev/null +++ b/example/test_problem/ELBDM/HaloMerger/FDMHaloDensityProfile @@ -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 diff --git a/example/test_problem/ELBDM/HaloMerger/HaloDensityProfile b/example/test_problem/ELBDM/HaloMerger/FitHaloDensityProfile similarity index 100% rename from example/test_problem/ELBDM/HaloMerger/HaloDensityProfile rename to example/test_problem/ELBDM/HaloMerger/FitHaloDensityProfile diff --git a/example/test_problem/ELBDM/HaloMerger/Input__TestProb_ParCloud b/example/test_problem/ELBDM/HaloMerger/Input__TestProb_ParCloud index 346635ff28..d65add99db 100644 --- a/example/test_problem/ELBDM/HaloMerger/Input__TestProb_ParCloud +++ b/example/test_problem/ELBDM/HaloMerger/Input__TestProb_ParCloud @@ -6,10 +6,12 @@ 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 @@ -17,7 +19,9 @@ 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 diff --git a/example/test_problem/ELBDM/HaloMerger/Make_DensityProfile.py b/example/test_problem/ELBDM/HaloMerger/Make_DensityProfile.py index ea6bea71d3..b7676624fe 100644 --- a/example/test_problem/ELBDM/HaloMerger/Make_DensityProfile.py +++ b/example/test_problem/ELBDM/HaloMerger/Make_DensityProfile.py @@ -1,5 +1,6 @@ import numpy as np import matplotlib.pyplot as plt +import os.path ################################################################################################################### @@ -99,22 +100,32 @@ def Soliton_fitting_analytical_dens(r, m22, r_c): ################################################################################################################### # save to file -np.savetxt( 'HaloDensityProfile', +filename = 'FitHaloDensityProfile' + +if os.path.exists( filename ): + print( '\nWARNING: file "%s" already exists and will be overwritten !!'%filename ) + +np.savetxt( filename, np.column_stack( (halo_densprof_radius, halo_densprof_density) ), - fmt=' %9.8e', - header=' r density' ) + fmt='%24.8e', + header='%22s %24s'%( 'r', 'density' ) ) + +filename = 'SolitonDensityProfile' -np.savetxt( 'SolitonDensityProfile', +if os.path.exists( filename ): + print( '\nWARNING: file "%s" already exists and will be overwritten !!'%filename ) + +np.savetxt( filename, np.column_stack( (soliton_densprof_radius, soliton_densprof_density) ), - fmt=' %9.8e', - header=' r density' ) + fmt='%24.8e', + header='%22s %24s'%( 'r', 'density' ) ) ################################################################################################################### ################################################################################################################### # plot to images fig = plt.figure() -ax = fig.add_subplot(111) +ax = fig.add_subplot( 111 ) # plot some important values for reference ax.plot( [halo_fitting_r_0, halo_fitting_r_0 ], [0.3*np.min(halo_densprof_density), 3.0*np.max(halo_densprof_density)], '--', color='grey', label=r'$r_0$' ) @@ -135,8 +146,8 @@ def Soliton_fitting_analytical_dens(r, m22, r_c): xy=(0.02,0.02), xycoords='axes fraction' ) # setting for the figure -ax.set_xscale('log') -ax.set_yscale('log') +ax.set_xscale( 'log' ) +ax.set_yscale( 'log' ) ax.set_xlim( 0.5*np.min(halo_densprof_radius), 2.0*np.max(halo_densprof_radius) ) ax.set_ylim( 0.5*np.min(halo_densprof_density), 2.0*np.max(halo_densprof_density) ) @@ -148,11 +159,18 @@ 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' ) + +filename_fig = 'fig_FitHaloDensityProfile.png' + +if os.path.exists( filename_fig ): + print( '\nWARNING: file "%s" already exists and will be overwritten !!'%filename_fig ) + +fig.savefig( filename_fig ) plt.close() + fig = plt.figure() -ax = fig.add_subplot(111) +ax = fig.add_subplot( 111 ) # plot some important values for reference ax.plot( [soliton_fitting_r_c, soliton_fitting_r_c ], [0.3*np.min(soliton_densprof_density), 3.0*np.max(soliton_densprof_density)], '--', color='grey', label=r'$r_{\rm c}$' ) @@ -170,8 +188,8 @@ def Soliton_fitting_analytical_dens(r, m22, r_c): xy=(0.02,0.02), xycoords='axes fraction' ) # setting for the figure -ax.set_xscale('log') -ax.set_yscale('log') +ax.set_xscale( 'log' ) +ax.set_yscale( 'log' ) ax.set_xlim( 0.5*np.min(soliton_densprof_radius), 2.0*np.max(soliton_densprof_radius) ) ax.set_ylim( 0.0001*np.max(soliton_densprof_density), 3.0*np.max(soliton_densprof_density) ) @@ -183,6 +201,12 @@ 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_SolitonDensityProfile.png' ) + +filename_fig = 'fig_SolitonDensityProfile.png' + +if os.path.exists( filename_fig ): + print( '\nWARNING: file "%s" already exists and will be overwritten !!'%filename_fig ) + +fig.savefig( filename_fig ) plt.close() ################################################################################################################### diff --git a/example/test_problem/ELBDM/HaloMerger/Make_ParCloud_ExtPotTable.py b/example/test_problem/ELBDM/HaloMerger/Make_ParCloud_ExtPotTable.py new file mode 100644 index 0000000000..4d0b84f1d7 --- /dev/null +++ b/example/test_problem/ELBDM/HaloMerger/Make_ParCloud_ExtPotTable.py @@ -0,0 +1,106 @@ +import numpy as np +import matplotlib.pyplot as plt +import os.path + + +################################################################################################################### +# 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 +filename = 'ParCloud_ExtPotTable' + +if os.path.exists( filename ): + print( '\nWARNING: file "%s" already exists and will be overwritten !!'%filename ) + +np.savetxt( filename, + np.column_stack( (ParCloud_ExtPot_table_radius, ParCloud_ExtPot_table_potential) ), + fmt='%23.8e', + header='%21s %23s'%( 'r', 'potential' ) ) +################################################################################################################### + + +################################################################################################################### +# 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 ) + +filename_fig = 'fig_ParCloud_ExtPotTable.png' + +if os.path.exists( filename_fig ): + print( '\nWARNING: file "%s" already exists and will be overwritten !!'%filename_fig ) + +fig.savefig( filename_fig ) +plt.close() +################################################################################################################### diff --git a/example/test_problem/ELBDM/HaloMerger/README b/example/test_problem/ELBDM/HaloMerger/README index ab84f05b10..d0811d5603 100644 --- a/example/test_problem/ELBDM/HaloMerger/README +++ b/example/test_problem/ELBDM/HaloMerger/README @@ -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 "python Make_ParCloud_ExtPotTable.py", + where the default is the potential of a uniform density sphere with 3x mass and 0.3x radius of the default CDM halo + c. Add options `--particle=true --gsl=true` in generate_make.sh to enable compilation options PARTICLE and SUPPORT_GSL + 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 @@ -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 `python plot__density_profile.py -s 0 -e 0` + 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 diff --git a/example/test_problem/ELBDM/HaloMerger/plot_script/plot__density_profile.py b/example/test_problem/ELBDM/HaloMerger/plot_script/plot__density_profile.py index de7e4bc900..5a5ebb6cad 100644 --- a/example/test_problem/ELBDM/HaloMerger/plot_script/plot__density_profile.py +++ b/example/test_problem/ELBDM/HaloMerger/plot_script/plot__density_profile.py @@ -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) ]) @@ -71,8 +71,8 @@ # save the profile to text file np.savetxt( '%s_%s_profile'%(ds,field), np.column_stack( (prof_dens.x.in_units('code_length').d, prof_dens[field].in_units('code_density').d) ), - fmt=' %9.8e', - header=' r (code_length) density (code_density)' ) + fmt='%24.8e', + header='%22s %24s'%( 'r (code_length)', 'density (code_density)' ) ) # load the reference profiles Ref_DensProf_r, Ref_DensProf_dens = np.loadtxt( Ref_DensProf_filename, skiprows=1, unpack=True ) diff --git a/example/test_problem/ELBDM/HaloMerger/plot_script/plot__enclosed_mass.py b/example/test_problem/ELBDM/HaloMerger/plot_script/plot__enclosed_mass.py index dd4b7642b8..489a6df28c 100644 --- a/example/test_problem/ELBDM/HaloMerger/plot_script/plot__enclosed_mass.py +++ b/example/test_problem/ELBDM/HaloMerger/plot_script/plot__enclosed_mass.py @@ -73,8 +73,8 @@ def _TotCellMass(field, data): # save the profile to text file np.savetxt( '%s_EnclosedTotalMass_profile'%(ds), np.column_stack( (prof_mass.x.in_units('code_length').d, prof_mass[field_mass].in_units('code_mass').d) ), - fmt=' %9.8e', - header=' r (code_length) mass (code_mass)' ) + fmt='%24.8e', + header='%22s %24s'%( 'r (code_length)', 'mass (code_mass)' ) ) # decide the units for plotting #UNIT_L_PLOT = 'code_length' diff --git a/src/TestProblem/ELBDM/HaloMerger/Init_TestProb_ELBDM_HaloMerger.cpp b/src/TestProblem/ELBDM/HaloMerger/Init_TestProb_ELBDM_HaloMerger.cpp index 35a8f4a8d3..1055024c73 100644 --- a/src/TestProblem/ELBDM/HaloMerger/Init_TestProb_ELBDM_HaloMerger.cpp +++ b/src/TestProblem/ELBDM/HaloMerger/Init_TestProb_ELBDM_HaloMerger.cpp @@ -54,6 +54,8 @@ static double **HaloMerger_Soliton_DensProf = NULL; // ar double *HaloMerger_ParCloud_DensProf_MaxR = NULL; // maximum radius for particles of each particle cloud int *HaloMerger_ParCloud_RSeed = NULL; // random seed for particles of each particle cloud long *HaloMerger_ParCloud_NPar = NULL; // number of particles of each particle cloud + int *HaloMerger_ParCloud_BuiltWithExtPot = NULL; // whether to consider an external potential for the velocity dispersion of each particle cloud + char (*HaloMerger_ParCloud_ExtPot_Filename)[MAX_STRING] = NULL; // filename of the external potential table of each particle cloud // ======================================================================================= #if ( MODEL == ELBDM && defined GRAVITY ) @@ -431,6 +433,8 @@ void SetParameter() HaloMerger_ParCloud_DensProf_MaxR = new double [HaloMerger_ParCloud_Num]; HaloMerger_ParCloud_RSeed = new int [HaloMerger_ParCloud_Num]; HaloMerger_ParCloud_NPar = new long [HaloMerger_ParCloud_Num]; + HaloMerger_ParCloud_BuiltWithExtPot = new int [HaloMerger_ParCloud_Num]; + HaloMerger_ParCloud_ExtPot_Filename = new char [HaloMerger_ParCloud_Num][MAX_STRING]; } // if ( HaloMerger_ParCloud_InitMode == 1 ) else Aux_Error( ERROR_INFO, "unsupported initialization mode (%s = %d) !!\n", @@ -451,6 +455,8 @@ void SetParameter() char HaloMerger_ParCloud_i_DensProf_MaxR[MAX_STRING]; // maximum radius for particles for the i-th particle cloud (must > 0.0) (HaloMerger_ParCloud_InitMode == 1 only) char HaloMerger_ParCloud_i_RSeed[MAX_STRING]; // random seed for setting particle position and velocity for the i-th particle cloud (must >= 0) (HaloMerger_ParCloud_InitMode == 1 only) char HaloMerger_ParCloud_i_NPar[MAX_STRING]; // number of particles for the i-th particle cloud (must >= 0) (HaloMerger_ParCloud_InitMode == 1 only) + char HaloMerger_ParCloud_i_BuiltWithExtPot[MAX_STRING]; // whether to consider an external potential for the velocity dispersion for the i-th particle cloud (HaloMerger_ParCloud_InitMode == 1 only) + char HaloMerger_ParCloud_i_ExtPot_Filename[MAX_STRING]; // filename of the external potential table for the i-th particle cloud (HaloMerger_ParCloud_InitMode == 1 only) (HaloMerger_ParCloud_i_BuiltWithExtPot == 1 only) for (int index_parcloud=0; index_parcloudAdd( HaloMerger_ParCloud_i_DensProf_MaxR, &HaloMerger_ParCloud_DensProf_MaxR[index_parcloud], 0.5*amr->BoxSize[0], Eps_double, NoMax_double ); ReadPara_ParCloud->Add( HaloMerger_ParCloud_i_RSeed, &HaloMerger_ParCloud_RSeed[index_parcloud], 123, 0, NoMax_int ); ReadPara_ParCloud->Add( HaloMerger_ParCloud_i_NPar, &HaloMerger_ParCloud_NPar[index_parcloud], (long)0, (long)0, NoMax_long ); + ReadPara_ParCloud->Add( HaloMerger_ParCloud_i_BuiltWithExtPot, &HaloMerger_ParCloud_BuiltWithExtPot[index_parcloud], 0, 0, 1 ); + ReadPara_ParCloud->Add( HaloMerger_ParCloud_i_ExtPot_Filename, HaloMerger_ParCloud_ExtPot_Filename[index_parcloud], NoDef_str, Useless_str, Useless_str ); } // if ( HaloMerger_ParCloud_InitMode == 1 ) else Aux_Error( ERROR_INFO, "unsupported initialization mode (%s = %d) !!\n", @@ -878,6 +888,15 @@ void SetParameter() Aux_Error( ERROR_INFO, "ParCloud_%d density profile file \"%s\" does not exist !!\n", index_parcloud_input, HaloMerger_ParCloud_DensProf_Filename[index_parcloud] ); + // check the external potential file exists + if ( HaloMerger_ParCloud_BuiltWithExtPot[index_parcloud] && !Aux_CheckFileExist(HaloMerger_ParCloud_ExtPot_Filename[index_parcloud]) ) + Aux_Error( ERROR_INFO, "ParCloud_%d \"BuiltWithExtPot\" is on, but external potential file \"%s\" does not exist !!\n", + index_parcloud_input, HaloMerger_ParCloud_ExtPot_Filename[index_parcloud] ); + + // rename the external potential filename if it is not needed + if ( !HaloMerger_ParCloud_BuiltWithExtPot[index_parcloud] ) + strcpy( HaloMerger_ParCloud_ExtPot_Filename[index_parcloud], "NONE" ); + // check whether the particle cloud touches the boundary of box for (int d=0; d<3; d++) { @@ -1140,13 +1159,14 @@ void SetParameter() if ( HaloMerger_ParCloud_InitMode == 1 ) { Aux_Message( stdout, "\n particle cloud density profile information:\n" ); - Aux_Message( stdout, " %7s %34s %16s %16s %16s\n", - "ID", "DensProf_Filename", "DensProf_MaxR", "RSeed", "NPar" ); + Aux_Message( stdout, " %7s %34s %16s %16s %16s %16s %34s\n", + "ID", "DensProf_Filename", "DensProf_MaxR", "RSeed", "NPar", "BuiltWithExtPot", "ExtPot_Filename" ); for (int index_parcloud=0; index_parcloud Total number of particles in all particle clouds = %16ld\n", "", HaloMerger_ParCloud_NPar_Total ); @@ -1437,6 +1457,8 @@ void End_HaloMerger() delete [] HaloMerger_ParCloud_DensProf_MaxR; delete [] HaloMerger_ParCloud_RSeed; delete [] HaloMerger_ParCloud_NPar; + delete [] HaloMerger_ParCloud_BuiltWithExtPot; + delete [] HaloMerger_ParCloud_ExtPot_Filename; } // if ( HaloMerger_ParCloud_InitMode == 1 ) else Aux_Error( ERROR_INFO, "unsupported initialization mode (%s = %d) !!\n", diff --git a/src/TestProblem/ELBDM/HaloMerger/Par_Init_ByFunction_HaloMerger.cpp b/src/TestProblem/ELBDM/HaloMerger/Par_Init_ByFunction_HaloMerger.cpp index b4af6f0c45..80629eacb3 100644 --- a/src/TestProblem/ELBDM/HaloMerger/Par_Init_ByFunction_HaloMerger.cpp +++ b/src/TestProblem/ELBDM/HaloMerger/Par_Init_ByFunction_HaloMerger.cpp @@ -11,6 +11,8 @@ extern char (*HaloMerger_ParCloud_DensProf_Filename)[MAX_STRING]; extern double *HaloMerger_ParCloud_DensProf_MaxR; extern int *HaloMerger_ParCloud_RSeed; extern long *HaloMerger_ParCloud_NPar; +extern int *HaloMerger_ParCloud_BuiltWithExtPot; +extern char (*HaloMerger_ParCloud_ExtPot_Filename)[MAX_STRING]; //------------------------------------------------------------------------------------------------------- @@ -103,7 +105,8 @@ void Par_Init_ByFunction_HaloMerger( const long NPar_ThisRank, const long NPar_A Cloud_Constructor.params.Cloud_RSeed = HaloMerger_ParCloud_RSeed[index_parcloud]; Cloud_Constructor.params.Cloud_Par_Num = HaloMerger_ParCloud_NPar[index_parcloud]; Cloud_Constructor.params.Cloud_R0 = 1.0; // will have no effect as long as the value is positive - Cloud_Constructor.params.AddExtPot = 0; // no external potential + Cloud_Constructor.params.AddExtPot = HaloMerger_ParCloud_BuiltWithExtPot[index_parcloud]; + strcpy( Cloud_Constructor.params.ExtPot_Table_Name, HaloMerger_ParCloud_ExtPot_Filename[index_parcloud] ); // initialize the particle cloud Cloud_Constructor.Init();