Grackle with COMOVING #349

YuriOku · 2024-09-07T10:10:09Z

In this PR, I have modified the code to use Grackle with COMOVING. The changes I made are

Define Grackle units in Grackle_Init.cpp and Grackle_Prepare.cpp
Convert the specific thermal energy between GAMER and Grackle units in Grackle_Prepare.cpp and Grackle_Close.cpp
Convert the timestep size to the physical time interval in Grackle_AdvanceDt.cpp
Disable normalization of passive scalars for electron and deuterium densities for consistency with Grackle in Init_Field.cpp

I then created a test problem, Grackle_Comoving. The test simulation evolves uniform density gas with an initial temperature of T=1e8 K from z=9 to z=0. The time evolution of density and temperature is recorded in Record__User, and the cooling curve can be computed from the record with the Python script Compute_CoolingCurve.py. There is a C program coolingrate_proper.c, which generates a cooling curve table using Grackle at each redshift, which can be compared with the result of the test simulation.

In the left panel in the figure below, the solid line shows the cooling curve obtained from the test simulation, and the points indicate the values at each redshift. The dotted lines are cooling curves at each redshift. The points lie just on the dotted lines, indicating that the thermal evolution is consistent with the Grackle cooling curve.

hyschive · 2024-09-09T03:59:29Z

@YuriOku Thanks for the contribution! @ChunYen-Chen will help review your PR.

ChunYen-Chen · 2024-09-09T10:21:34Z

Hi @YuriOku, I will review your PR this week and probably finish the review before the weekend.

ChunYen-Chen · 2024-09-10T03:07:40Z

Hi @YuriOku, can you share the grackle data CloudyData_noUVB.h5 with me?

YuriOku · 2024-09-10T05:05:50Z

Hi Chun-Yen, Please find the data file at https://github.com/grackle-project/grackle_data_files/tree/main/input. Thanks!

…

On Tue, Sep 10, 2024 at 11:08 AM Chun-Yen Chen ***@***.***> wrote: Hi @YuriOku <https://github.com/YuriOku>, can you share the grackle data CloudyData_noUVB.h5 with me? — Reply to this email directly, view it on GitHub <#349 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AO55JJFO5D4IT7KZJ4XIKKDZVZPBDAVCNFSM6AAAAABNZ4SOJCVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGMZZGUYTQOJYHE> . You are receiving this because you were mentioned.Message ID: ***@***.***>

ChunYen-Chen

Hi @YuriOku,

The PR looks great and I can reproduce your results. I have two questions about your analysis script/tool:

Is there any specific reason you choose -O2 to compile coolingrate_proper.c?
Can you run coolingrate with single precision grackle? If not, there should be some bugs in coolingrate.c.
- If your tool does not support single precision, please add a check of FLOAT8 in Init_TestProb_Grackle_Comoving.cpp

Before merging, please fix/change the following:

Please fix all bugs with comments starting with [BUG].
Please update your new code to GAMER code style with the comments starting with [Style]. You can resolve most of the code-style comments by copying the code I wrote in the comment.
Please remove all the trailing spaces. You can search them by command:
```
grep -r -n "[[:blank:]]$" ./
```
Please provide a simple script for downloading the CloudyData_noUVB.h5, e.g. download_ic.sh in ClusterMerger test problem.

The rest of the comments without tags are some minor change suggestions or requests.

Please let me know if any of my comments is unclear.

example/test_problem/Hydro/Grackle_Comoving/Input__DumpTable

example/test_problem/Hydro/Grackle_Comoving/README

src/TestProblem/Hydro/Grackle_Comoving/Init_TestProb_Grackle_Comoving.cpp

src/Grackle/Grackle_AdvanceDt.cpp

example/test_problem/Hydro/Grackle_Comoving/Input__Parameter

src/Grackle/Grackle_Init.cpp

src/TestProblem/Hydro/Grackle_Comoving/Init_TestProb_Grackle_Comoving.cpp

YuriOku · 2024-09-10T13:06:36Z

Hi @ChunYen-Chen, thank you for your review! I will address your points.

Is there any specific reason you choose -O2 to compile coolingrate_proper.c?

There is no reason for the option.

Can you run coolingrate with single precision grackle? If not, there should be some bugs in coolingrate.c.
If your tool does not support single precision, please add a check of FLOAT8 in Init_TestProb_Grackle_Comoving.cpp

I have error in running coolingrate with single precision grackle. I will add a check of FLOAT8 in Init_TestProb_Grackle_Comoving.cpp

ChunYen-Chen · 2024-09-10T13:48:38Z

Hi @YuriOku,

Is there any specific reason you choose -O2 to compile coolingrate_proper.c?

There is no reason for the option.

How about changing it to -O3?

Can you run coolingrate with single precision grackle? If not, there should be some bugs in coolingrate.c.
If your tool does not support single precision, please add a check of FLOAT8 in Init_TestProb_Grackle_Comoving.cpp

I have error in running coolingrate with single precision grackle. I will add a check of FLOAT8 in Init_TestProb_Grackle_Comoving.cpp

Great, please also update the README in your test problem.

hyschive · 2024-09-10T14:43:11Z

@YuriOku @ChunYen-Chen

It's fine to adopt -O2, which is also the default set-up in GAMER (e.g., configs/spock_intel.config).
GAMER prefers single precision for higher GPU performance. So it would be great if both coolingrate_proper.c and Init_TestProb_Grackle_Comoving.cpp could support both single and double precision. Besides, by doing so, we can run GAMER and coolingrate_proper.c using the same floating-point precision so that one only needs to compile Grackle once.
Related to the above point, we should validate whether Grackle, GAMER (presumably in Grackle_Init.cpp , and coolingrate_proper.c adopt the same floating-point precision. Does Grackle provide an API to inquire this information?

YuriOku · 2024-09-11T08:19:26Z

Hi @hyschive,

I see. I will try to run the code with a single precision.

The floating-point precision of Grackle can be checked by sizeof(gr_float).

hyschive · 2024-09-11T08:30:13Z

So maybe we should re-enable this check. I forget why it was commented out. Please test it.

YuriOku · 2024-09-13T08:38:51Z

@hyschive @ChunYen-Chen
I have tried to my code in single-float precision, but I could not make it. The Grackle's example C code also fails in single-float precision. I could not identify the cause, but I think this is a bug in Grackle.

Given that Grackle is not recommended run in single-float precision (Grackle Installation Doc), I think it would be better to run Grackle in double precision regardless of the precision of fluid solver. I propose to define h_Che_Array in gr_float, and then, the type conversion is automatically made in Grackle_Prepare and Grackle_Close. An option like GRACKLE_FLOAT can be added to check the precision of Grackle. How do you think?

hyschive · 2024-09-13T09:00:28Z

@YuriOku

I agree with your suggestion and have just filed a feature request Grackle floating-point precision #354.
On the other hand, IIUC, @ChunYen-Chen already performed the Grackle_Comoving test in both single and double precision and the results are consistent. So I guess the problem lies in coolingrate_proper.c. @ChunYen-Chen Do you agree?

hyschive · 2024-09-13T09:04:09Z

@ChunYen-Chen After merging this PR, please help file a follow-up PR to address #354. Thanks.

YuriOku · 2024-09-14T13:49:46Z

@ChunYen-Chen
I removed the trailing space and added a script to download CloudyData_noUVB.h5. I also updated the code according to your suggestion, leaving two conversations unresolved that I am not sure I have resolved. Could you review them?

After these fixes, I have confirmed that I can reproduce the same results in the test problem Grackle_Comoving.

ChunYen-Chen · 2024-09-14T14:46:00Z

@YuriOku Thanks for the update. The PR looks great and I think it is ready to merge.
@hyschive Yes, I do believe the problem is in coolingrate_proper.c.

hyschive · 2024-12-31T09:22:45Z

@YuriOku Sorry for the delayed progress. We’re planning to merge this PR soon. Could you please update it to resolve the conflicts first?

YuriOku · 2024-12-31T10:02:09Z

resolved the conflict

hsinhaoHHuang

Hi @YuriOku,
I have finished the review of this PR.
It is exciting to know GAMER can support Grackle with COMOVING. Thank you for implementing it.

I have raised some questions and added some suggestions in the comments. Could you help me take a look when you have time?
Thanks!

hsinhaoHHuang · 2025-01-03T05:58:01Z

src/Grackle/Grackle_Prepare.cpp

+   Che_Units.velocity_units       = UNIT_V;
+   Che_Units.a_units              = 1.0;
+   Che_Units.a_value              = Time[lv];
+#  endif

   const int  Size1pg          = CUBE(PS2);
   const int  Size1v           = NPG*Size1pg;


Please help align "=" in the next line with Line. 107 and 106.

hsinhaoHHuang · 2025-01-03T06:26:43Z

src/Grackle/Grackle_Close.cpp

+//          convert from the proper frame to the comoving specific internal energy
+            Eint = Ptr_sEint[idx_pg] * Dens * SQR(Time[lv]);


To avoid confusion, how about "convert from the proper frame specific internal energy to the comoving internal energy density"?

hsinhaoHHuang · 2025-01-03T06:34:57Z

src/Grackle/Grackle_Init.cpp

+// see https://grackle.readthedocs.io/en/latest/Interaction.html#comoving-coordinates and
+// https://github.com/grackle-project/grackle/issues/192
+// --> the density and length units are conversion factors from code values to cgs in proper coordinates,
+//     while the time and velocity units should remain constant, i.e., velocity_units = length_units / (Time * time_units)


Although the "Time" here refers to Time in GAMER, which is scale_factor, I think it would cause less confusion if we wrote:
"velocity_units = length_units / (a_value * time_units)".

hsinhaoHHuang · 2025-01-03T09:54:25Z

src/Grackle/Grackle_Init.cpp

+// https://github.com/grackle-project/grackle/issues/192
+// --> the density and length units are conversion factors from code values to cgs in proper coordinates,
+//     while the time and velocity units should remain constant, i.e., velocity_units = length_units / (Time * time_units)
+// --> the specific energy passed to Grackle is u = k_B * T / ((gamma - 1) * mu * m_p)


How about adding a reminder here that the specific energy has to be converted to the proper frame as the velocity_units does not change with a?

hsinhaoHHuang · 2025-01-03T10:39:10Z

include/GAMER.h

@@ -58,6 +58,8 @@
 #  define CONFIG_BFLOAT_4
 #endif

+#define OMIT_LEGACY_INTERNAL_GRACKLE_FUNC


Could you provide more information on why this is necessary? I cannot find this in my Grackle.

hsinhaoHHuang · 2025-01-04T13:48:04Z

src/TestProblem/Hydro/Grackle_Comoving/Init_TestProb_Grackle_Comoving.cpp

+      if ( calculate_cooling_time( &Che_Units, &my_fields, my_cooling_time ) == 0 )
+         Aux_Error( ERROR_INFO, "Error in calculate_cooling_time.\n" );
+//    calculate temperature
+      if ( calculate_temperature(&Che_Units, &my_fields, my_temperature ) == 0 )


Suggested change

if ( calculate_temperature(&Che_Units, &my_fields, my_temperature ) == 0 )

if ( calculate_temperature( &Che_Units, &my_fields, my_temperature ) == 0 )

hsinhaoHHuang · 2025-01-04T14:01:30Z

src/TestProblem/Hydro/Grackle_Comoving/Init_TestProb_Grackle_Comoving.cpp

+      const double dt_SubStep        = Mis_dTime2dt( Time[0], dTime_Base ) * SQR(Time[0]) * UNIT_T;
+      const double temperature_units = get_temperature_units(&Che_Units); // set temperature units
+      const double mu                = my_temperature[0] / (my_fields.internal_energy[0] * (my_gamma[0] - 1.) * temperature_units);
+      const double n                 = Dens / CUBE(Time[0]) * UNIT_D / mu / Const_mp;
+      const double Edens             = Eint * UNIT_P / SQR(Time[0]) / CUBE(Time[0]);
+      const double Temp              = my_temperature[0];
+      const double Lcool             = Edens / fabs(my_cooling_time[0] * UNIT_T) / n / n;


I think it would be better to provide some comments to describe these variables.

hsinhaoHHuang · 2025-01-04T14:32:55Z

src/TestProblem/Hydro/Grackle_Comoving/Init_TestProb_Grackle_Comoving.cpp

+      const double Temp              = my_temperature[0];
+      const double Lcool             = Edens / fabs(my_cooling_time[0] * UNIT_T) / n / n;
+
+      fprintf( File_User, "%14.7e%14ld%3s%14.7e%14.7e%14.7e%14.7e%14.7e%14.7e", Time[0], Step, "", dt_SubStep, n, mu, Temp, Edens, Lcool );


I think one thing we should remind the user of is that these values are in the physical frame, which is different from other output data where the values are in the comoving frames when COMOVING is enabled.
Actually, I think it would be better to output the values (e.g., number density, temperature, and internal energy density) both in the comoving and physical frames. It could possibly avoid confusion and also be helpful for analysis.

hsinhaoHHuang · 2025-01-04T14:41:42Z

src/TestProblem/Hydro/Grackle_Comoving/Init_TestProb_Grackle_Comoving.cpp

+         if ( Aux_CheckFileExist(FileName) )    Aux_Message( stderr, "WARNING : file \"%s\" already exists !!\n", FileName );
+
+         FILE *File_User = fopen( FileName, "a" );
+         fprintf( File_User, "#%13s%14s%3s%14s%14s%14s%14s%14s",  "Time", "Step", "", "dt [s]", "n [1/cc]", "mu", "Temp [K]", "Edens [erg/cc], Lcool[erg cm^3 s^-1]" );


It is true that Time refers to the scale factor in COMOVING, and it is used in other Record__* files. But with the presence of dt and also all the other variables converted into the physical frame, would you think it is clearer to just use scale factor a for Time in the header?

hsinhaoHHuang · 2025-01-05T07:30:26Z

src/TestProblem/Hydro/Grackle_Comoving/Init_TestProb_Grackle_Comoving.cpp

+      const double n                 = Dens / CUBE(Time[0]) * UNIT_D / mu / Const_mp;
+      const double Edens             = Eint * UNIT_P / SQR(Time[0]) / CUBE(Time[0]);
+      const double Temp              = my_temperature[0];
+      const double Lcool             = Edens / fabs(my_cooling_time[0] * UNIT_T) / n / n;


If I understand correctly, the cooling rate calculated here is not used to compare the results of the simulation (in Compute_CoolingCurve.py). Instead, we use coolingrate_proper.c to calculate the cooling rate at each redshift to validate the evolution of internal energy in simulations.

It seems that coolingrate_proper.c performs a function similar to the cooling rate calculation here in GAMER, but operates in proper frames and with different temperatures.
Would it be possible to integrate the functionality of coolingrate_proper.c into Aux_Record_GrackleComoving() in GAMER?
We could iterate through various temperatures and get the cooling rate in proper frames as a table of temperatures (possibly output in a separate file) at each step.

This approach would allow us to plot the results after the simulations are complete, eliminating the need to use Grackle for additional calculations.
It would be beneficial if users did not need to compile an additional C program to generate the cooling rate table.
This integration could also prevent potential inconsistencies in parameters used in GAMER and Grackle (e.g., metallicity, primordial_chemistry).

YuriOku added 2 commits September 5, 2024 21:19

grackle with comoving

cb8d46a

added analysis scripts

123bd1d

hyschive requested a review from ChunYen-Chen September 9, 2024 03:45

hyschive assigned ChunYen-Chen Sep 9, 2024

hyschive added enhancement hydro Hydrodynamics labels Sep 9, 2024

ChunYen-Chen requested changes Sep 10, 2024

View reviewed changes

YuriOku added 8 commits September 13, 2024 17:24

Update README

3aad958

Update clean.sh

cc1b3b3

Update Compute_CoolingCurve.py

c23ab25

updated style

9b8ed51

updated style

25da4ce

removed trailing space

39bd356

updated README

275cde5

updated README

a8a83a7

ChunYen-Chen approved these changes Sep 14, 2024

View reviewed changes

hyschive requested a review from hsinhaoHHuang December 31, 2024 09:22

hyschive assigned hsinhaoHHuang Dec 31, 2024

Merge branch 'main' into grackle

a1081f1

hsinhaoHHuang reviewed Jan 5, 2025

View reviewed changes

hyschive mentioned this pull request Jan 11, 2025

Support Grackle in cosmological simulations #242

Open

7 tasks

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Grackle with COMOVING #349

Grackle with COMOVING #349

YuriOku commented Sep 7, 2024

hyschive commented Sep 9, 2024

ChunYen-Chen commented Sep 9, 2024

ChunYen-Chen commented Sep 10, 2024

YuriOku commented Sep 10, 2024 via email

ChunYen-Chen left a comment •

edited

Loading

YuriOku commented Sep 10, 2024

ChunYen-Chen commented Sep 10, 2024

hyschive commented Sep 10, 2024

YuriOku commented Sep 11, 2024

hyschive commented Sep 11, 2024

YuriOku commented Sep 13, 2024

hyschive commented Sep 13, 2024

hyschive commented Sep 13, 2024

YuriOku commented Sep 14, 2024 •

edited

Loading

ChunYen-Chen commented Sep 14, 2024

hyschive commented Dec 31, 2024

YuriOku commented Dec 31, 2024

hsinhaoHHuang left a comment

hsinhaoHHuang Jan 3, 2025

hsinhaoHHuang Jan 3, 2025

hsinhaoHHuang Jan 3, 2025

hsinhaoHHuang Jan 3, 2025

hsinhaoHHuang Jan 3, 2025

hsinhaoHHuang Jan 4, 2025

hsinhaoHHuang Jan 4, 2025

hsinhaoHHuang Jan 4, 2025

hsinhaoHHuang Jan 4, 2025

hsinhaoHHuang Jan 5, 2025

		// convert from the proper frame to the comoving specific internal energy
		Eint = Ptr_sEint[idx_pg] * Dens * SQR(Time[lv]);

	if ( calculate_temperature(&Che_Units, &my_fields, my_temperature ) == 0 )
	if ( calculate_temperature( &Che_Units, &my_fields, my_temperature ) == 0 )

Grackle with COMOVING #349

Are you sure you want to change the base?

Grackle with COMOVING #349

Conversation

YuriOku commented Sep 7, 2024

hyschive commented Sep 9, 2024

ChunYen-Chen commented Sep 9, 2024

ChunYen-Chen commented Sep 10, 2024

YuriOku commented Sep 10, 2024 via email

ChunYen-Chen left a comment • edited Loading

Choose a reason for hiding this comment

YuriOku commented Sep 10, 2024

ChunYen-Chen commented Sep 10, 2024

hyschive commented Sep 10, 2024

YuriOku commented Sep 11, 2024

hyschive commented Sep 11, 2024

YuriOku commented Sep 13, 2024

hyschive commented Sep 13, 2024

hyschive commented Sep 13, 2024

YuriOku commented Sep 14, 2024 • edited Loading

ChunYen-Chen commented Sep 14, 2024

hyschive commented Dec 31, 2024

YuriOku commented Dec 31, 2024

hsinhaoHHuang left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ChunYen-Chen left a comment •

edited

Loading

YuriOku commented Sep 14, 2024 •

edited

Loading