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.

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

Add DynamicalFriction Test Problem #252

Open
wants to merge 39 commits into
base: main
Choose a base branch
from

Conversation

ej-1211
Copy link

@ej-1211 ej-1211 commented Oct 17, 2023

No description provided.

@hyschive hyschive requested review from hsinhaoHHuang and removed request for koarakawaii January 12, 2024 02:12
Copy link
Contributor

@hsinhaoHHuang hsinhaoHHuang left a comment

Choose a reason for hiding this comment

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

Hi @ej-1211 ,

I am really sorry for the delayed review. 😢
Thank you for your patience.
I've had a quick look at this PR recently.

I can run the three default cases and reproduce similar results in this page.
However, the figure for the Quick_Burkert is not 100% identical with the reference. Here is my results:
Compare_result_20240201-1040
I am not sure why and I didn't figure out what may went wrong.

I provided some suggestions in the inline comments. Please take a look.
Please let me know if my messages are not clear or I misunderstood anything. 🙏 Thanks!

example/test_problem/Hydro/DynamicalFriction/clean.sh Outdated Show resolved Hide resolved
example/test_problem/Hydro/DynamicalFriction/clean.sh Outdated Show resolved Hide resolved
Comment on lines 14 to 24
# (1) HALO and STELLAR Parameters
HALO_TYPE NFW #* [Burkert]
HALO_RHO_0 1.96311939345708e-03 #* [Msun/pc^3][665.3026853292199e-3]
HALO_Rs 5.48799489144317e+03 #*Scale radius for the profile [pc][0.25e3]
HALO_Rt 35e3 #*Truncate radius for the halo [pc][6e3]

STELLAR_TYPE None #*[None]
STELLAR_RHO_0 0 #*[Msun/pc^3][0]
STELLAR_Rs 0 #*[pc][0]
STELLAR_Rt 0 #*[pc][0]

Copy link
Contributor

Choose a reason for hiding this comment

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

I found there are several parameters (e.g. HALO_*, STELLAR_*, PURE_TABLE) in Input__TestProb that are not used and are even not read by GAMER. This is weird for me.
The */+/! notation for the variables is also quite not consistent with other test problems in GAMER. I think the input and the output for the Python script are put in the same file is confusing.
Also, it would be better if the default run is as simple as possible such that users can run it without doing many steps.

Therefore, I suggest two workflows for this test problem to improve it.
I personally prefer the 1st one because it is simpler for the user.
But I think the 2nd one is easier to implement from your current code.


The 1st suggestion:

  • Put the calculation in A_SET_IC.py (generating density profile and deciding initial parameters) inside the GAMER source code.
    • Maybe take the Plummer test problem as the reference, see src/TestProblem/Hydro/Plummer/Par_Init_ByFunction_Plummer.cpp.
  • Combine C_Chandrasekhar_dynamical_friction.py and D_Compare.py into the same file to handle the simulation analysis.

Steps:

  1. In Input__Testprob
    • GC_SmallGas
    • FIX_CENTER
    • SEARCH_RADIUS
    • Cloud_Center
    • HALO_TYPE, HALO_Rs, HALO_Rt, HALO_Rho0
    • STELLAR_TYPE, STELLAR_Rs, STELLAR_Rt, STELLAR_Rho0
    • GC_MASS, GC_Initial_R
    • PURE_TABLE
  2. In Init_TestProb_Hydro_DynamicalFriction()
    • Goal:
      • to read the parameters to set the fluid
      • to read the parameters to record the center
      • to read the parameters for the GC
    • Input:
      • GC_SmallGas
      • GC_MASS, GC_Initial_R
      • FIX_CENTER
      • SEARCH_RADIUS
  3. Par_Init_ByFunction_GC()
    • Goal:
      • to set the initial condition of the halo
      • to set the initial condition of the GC
      • to set the derived parameters from the input model
    • Input
      • Cloud_Center
      • HALO_TYPE, HALO_Rs, HALO_Rt, HALO_Rho0
      • STELLAR_TYPE, STELLAR_Rs, STELLAR_Rt, STELLAR_Rho0
      • GC_MASS, GC_Initial_R
      • PURE_TABLE
    • Derived Parameters
      • Density profile and mass profile of the halo
      • GC_POS
      • GC_VEL
      • Cloud_BulkVel
      • Cloud_Rho0
      • Cloud_R0
      • Cloud_MaxR
  4. Init_ExtPot_GC()
    • Goal:
      • to set the external potential
    • Input:
      • DynamicialFriction_ExtPot_Rho0
      • DynamicalFriction_ExtPot_R
      • DynamicalFriction_ExtPot_Center
      • DynamicalFriction_ExtPot_MFrac
  5. Aux_Record_User_GC()
    • Goal:
      • to record the position of minimum potential
      • to record the position of GC
    • Input:
      • FIX_CENTER
      • SEARCH_RADIUS
    • Output:
      • table of the position of the minimum potential
      • table of the position of the GC
  6. plot_compare.py
    • Goal:
      • to provide the Chandrasekhar prediction
      • to plot the comparison of Chandrasekhar's prediction and the simulation results.
    • Input:
      • Chandrasekhar_result_r-t
      • table of the position of the minimum potential
      • table of the position of the GC
      • BOX_SIZE, NXO_TOT_X, END_T
      • LN_LAMBDA_BMAX_FACTOR, LN_LAMBDA_BMIN_FACTOR
      • GC_Sortening
    • Output:
      • Chandrasekhar_result_r-t
      • radius-time comparison figures

The 2nd suggestion:

  • Combine A_SET_IC.py and C_Chandrasekhar_dynamical_friction.py to calculate the initial parameters and analytical prediction.
  • Set the parameters manually in Input__TestProb and only feed needed parameters to GAMER.
  • Use a script to analyze the results.

Steps

  1. Use a Python script to get the initial parameters
    • Goal:
      • to provide a density profile of the halo as the initial condition
      • to provide the initial parameters for halos and GC
      • to provide the Chandrasekhar prediction
    • Input
      • Cloud_Center
      • HALO_TYPE, HALO_Rs, HALO_Rt, HALO_Rho0
      • STELLAR_TYPE, STELLAR_Rs, STELLAR_Rt, STELLAR_Rho0
      • GC_MASS, GC_Initial_R
      • PURE_TABLE
      • BOX_SIZE, NXO_TOT_X, END_T
      • LN_LAMBDA_BMAX_FACTOR, LN_LAMBDA_BMIN_FACTOR
      • GC_Sortening
    • Output
      • Density profile file
      • Mass profile file
      • GC_POS
      • GC_VEL
      • Cloud_BulkVel
      • Cloud_Rho0
      • Cloud_R0
      • Cloud_MaxR
      • Chandrasekhar_result_r-t
  2. In Input__TestProb
    • GC_SmallGas
    • FIX_CENTER
    • SEARCH_RADIUS
    • Density profile of the halo
      • Here, maybe we can provide the generated density profiles in this test problem directory and the recommended parameters in README for the three default cases to save the user’s time.
    • Cloud_BulkVel
    • Cloud_Rho0
    • Cloud_R0
    • Cloud_MaxR
    • GC_MASS
    • GC_POS
    • GC_VEL
  3. Init_TestProb_Hydro_DynamicalFriction()
    • Goal:
      • to read the parameters to set the fluid
      • to read the parameters to record the center
      • to read the parameters for the GC
    • Input:
      • GC_SmallGas
      • FIX_CENTER
      • SEARCH_RADIUS
      • GC_MASS
      • GC_POS
      • GC_VEL
  4. Par_Init_ByFunction_GC()
    • Goal:
      • to set the initial condition of the halo
      • to set the initial condition of the GC
    • Input
      • Density profile of the halo
      • Cloud_BulkVel
      • Cloud_Rho0
      • Cloud_R0
      • Cloud_MaxR
      • GC_MASS
      • GC_POS
      • GC_VEL
  5. Init_ExtPot_GC()
    • Goal:
      • to set the external potential
    • Input:
      • DynamicialFriction_ExtPot_Rho0
      • DynamicalFriction_ExtPot_R
      • DynamicalFriction_ExtPot_Center
      • DynamicalFriction_ExtPot_MFrac
  6. Aux_Record_User_GC()
    • Goal:
      • to record the position of the minimum potential
      • to record the position of the GC
    • Input:
      • FIX_CENTER
      • SEARCH_RADIUS
    • Output:
      • table of the position of the minimum potential
      • table of the position of the GC
  7. plot_compare.py
    • Goal:
      • to plot the comparison of Chandrasekhar's prediction and the simulation results.
    • Input:
      • Chandrasskhar_result_r-t
      • table of the position of the minimum potential
      • table of the position of the GC
      • BOX_SIZE, NX0_TOT_X0
    • Output:
      • radius-time comparison figures

@hyschive
Copy link
Contributor

hyschive commented Feb 5, 2024

@JeiwHuang Thanks for the extremely careful and thorough review!

@ej-1211: I agree with @JeiwHuang's most comments. For example,

  • Simplify the initial condition construction, if possible, by utilizing the built-in tools in gamer (e.g., Aux_ComputeProfile()).
  • Make variable names more descriptive
  • Add comments whenever appropriate
  • Follow the gamer coding format (e.g., no tabs, three-space indent, code alignment)
  • Minimize redundant code lines

If you have any questions or different opinions, please feel free to speak out 😉

@ej-1211
Copy link
Author

ej-1211 commented Mar 7, 2024

  • Followed the 1st suggestion above (working but not in GAMER coding style yet)
  • Aux_Record_User_DynamicalFriction() in Init_TestProb_Hydro_DynamicalFriction.cpp is done
    • Remove unnesassary Record__RESULT_Previous_Center.txt
    • changed RESULT to Result
    • remove .txt from the output file name
    • GAMER coding style done
      • three space indenting
      • { } in if else statement
    • Aligned the Record__Result_* header and added #
    • changed FixCenter to boolean and turn to if ( FixCenter )

@ej-1211
Copy link
Author

ej-1211 commented Mar 9, 2024

  • Separate Input__TestProb to Input__TestProb and Input__Profile_Params
    • Input__Profile_Params is for ParticleEquilibriumIC Test problem
    • reduced the amount of unregcognize parameter
  • Update the input file in /example
    • Only save one set of Input__*
  • Add a feature that can automatically update Cloud_BulkVel* in Input__Profile_Params, which is calculated in the GAMER routine
  • Add "Stellar" part back

@ej-1211
Copy link
Author

ej-1211 commented Mar 11, 2024

  • Change from /b_min_factor and /b_max_factor to *b_min_factor and *b_max_factor
  • Throw away the datetime in the result file
  • read the parameter from the log file and the Record__Note File
  • removed A_ C_ X_

@ej-1211
Copy link
Author

ej-1211 commented Mar 26, 2024

@hsinhaoHHuang (@hyschive) Sorry for the late update, I have chosen the 1st suggestion above and have already fixed most of the problems (comments and suggestions) above, besides the fixing, in this version, we

  • Add the function of subtracting GC's potential so that it will be more accurate to find the center
  • Three cases (NFW, Burkert, DoublePowerLaw) are tested again and results are inside README
  • deleted the unnecessary script (python and .cpp) and unsupported feature

Please review my pull request again and please tell me if you have any suggestion!
Thank you very much !!!!

Copy link
Contributor

@hsinhaoHHuang hsinhaoHHuang left a comment

Choose a reason for hiding this comment

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

Hi @ej-1211, Thanks for your update!

I have finished my second review and left some comments.
Please take another look, and please let me know if I didn't make myself clear.
Thanks!

Here is the result I got from the default set-up (DoublePowerLaw)
The vertical lines look a little weird to me but I am not sure whether this is expected.
Compare_result

SEARCH_RADIUS 2 # search radius for the potential minimum search (unit:spatial resolution)
# If FIX_CENTER = 1, this parameter is useless
# (1) HALO and STELLAR Parameters
HALO_TYPE DoublePowerLaw # the type of the halo (NFW,Burkert,King,DoublePowerLaw)[DoublePowerLaw]
Copy link
Contributor

Choose a reason for hiding this comment

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

How about adding space between different halo types (i.e. (NFW, Burkert, King, DoublePowerLaw))?

Copy link
Author

Choose a reason for hiding this comment

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

@hsinhaoHHuang Do you mean in the comment area? If the answer is yes, I have fixed it!

Comment on lines +179 to +192
Aux_Message( stdout, " GC small Gas = %13.5e\n", GC_SmallGas );
Aux_Message( stdout, " GC initial radius = %13.5e\n", GC_R );
Aux_Message( stdout, " GC Mass = %13.5e\n", GC_MASS );
Aux_Message( stdout, " Fix Center = %d\n", FixCenter );
Aux_Message( stdout, " Search Radius = %13.5e\n", SearchRadius );
Aux_Message( stdout, " Pure Table = %d\n", PURE_TABLE );
Aux_Message( stdout, " Density Table Name = %s\n", TableName );
Aux_Message( stdout, " Halo Type = %s\n", HaloType );
Aux_Message( stdout, " HALO_RHO_0 = %13.5e\n", Halo_Rho0 );
Aux_Message( stdout, " HALO_Rs = %13.5e\n", Halo_Rs );
Aux_Message( stdout, " HALO_Rt = %13.5e\n", Halo_Rt );
Aux_Message( stdout, " Halo_Profle_Param_a = %13.5f\n", Halo_Profile_Param_a );
Aux_Message( stdout, " Halo_Profle_Param_b = %13.5f\n", Halo_Profile_Param_b );
Aux_Message( stdout, " Halo_Profle_Param_c = %13.5f\n", Halo_Profile_Param_c );
Copy link
Contributor

Choose a reason for hiding this comment

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

How about using %13.7e for better alignment?

if ( MPI_Rank == 0 )
{
char Filename[MAX_STRING];
sprintf( Filename, "%s", "Record__Result_Center" );
Copy link
Contributor

Choose a reason for hiding this comment

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

How about specifying that it is the minimum potential in the filename?

{
fprintf(File, "%15s\t%15s\t%15s\t%15s\n", "# Time", " CenterX", " CenterY", " CenterZ");
}
fprintf(File, "%15.7f\t%15.7e\t%15.7e\t%15.7e\n",
Copy link
Contributor

Choose a reason for hiding this comment

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

How about also printing the Time as %15.7e?

src/TestProblem/Hydro/DynamicalFriction/Calculate_IC.cpp Outdated Show resolved Hide resolved
Copy link
Contributor

Choose a reason for hiding this comment

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

I got this warning:

IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
  If increasing the limit yields no improvement it is advised to analyze 
  the integrand in order to determine the difficulties.  If the position of a 
  local difficulty can be determined (singularity, discontinuity) one will 
  probably gain from splitting up the interval and calling the integrator 
  on the subranges.  Perhaps a special-purpose integrator should be used.

@hyschive
Copy link
Contributor

@ej-1211 Just a reminder. Let's get it done this summer 💪

@hyschive hyschive requested a review from hsinhaoHHuang January 4, 2025 07:11
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement particle Particles test Test problems
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants