Grackle with COMOVING

example/test_problem/Hydro/Grackle_Comoving/Input__TestProb

@@ -0,0 +1,3 @@
+# problem-specific runtime parameters
+GrackleComoving_InitialTemperature 1e8        # initial temperature (set to a 1e8K to draw a cooling function of T<1e8K)
+GrackleComoving_InitialMetallicity 1.295e-2   # assuming solar metallicity

example/test_problem/Hydro/Grackle_Comoving/README

@@ -0,0 +1,29 @@
+Compilation flags:
+========================================
+Enable : MODEL=HYDRO, COMOVING, SUPPORT_GRACKLE
+Disable: PARTICLE
+
+
+Default setup:
+========================================
+1. Default resolution = 16^3
+2. Uniform density = 1.0 (cosmic mean density)
+3. Use the cloudy cooling table for H and He (GRACKLE_PRIMORDIAL = 0)
+   --> Hydrogen's Lyman alpha cooling peak at 10^4 K disappears if non-equilibrium chemistry is used (GRACKLE_PRIMORDIAL > 0) due to the short cooling time compared to the equilibrium time scale.
+4. Use the user-defined timestep criterion (OPT__DT_USER = 1)
+   --> limit the timestep by the cooling time
+
+Note:
+========================================
+1. This test problem computes the thermal evolution on comoving coordinate with the Grackle library
+2. At the end of the simulation, the code will generate the cooling rate table computed on proper coordinates for comparison
+3. Run the python script `plot/Compute_CoolingCurve.py` to draw the cooling functions
+   --> The python script computes the cooling rate from the simulation output and compares it with the reference table
+4. The floating point precision of Grackle must be same with GAMER to compile the code successfully
+5. Run the script `download_datatable.sh` to download the cooling rate table `CloudyData_noUVB.h5` from the Grackle GitHub repository
+6. It is recommended to use double-float precision for the Grackle library
+7. The values of the mean molecular weight and the temperature in the test problems can be different from the original output of GAMER
+   --> The mean molecular weight is computed from the Hydrogen and Helium ionization fraction in Grackle
+       --> can be different from the value of MOLECULAR_WEIGHT in Input__Parameter and Record__Note
+   --> The temperature is calculated considering the mean molecular weight and the total energy density in Grackle
+       --> can be different from the temperature field output by OPT__OUTPUT_TEMP

example/test_problem/Hydro/Grackle_Comoving/clean.sh

@@ -0,0 +1,6 @@
+rm -f Record__Note Record__Timing Record__TimeStep Record__PatchCount Record__Dump Record__MemInfo Record__L1Err \
+      Record__Conservation Data* stderr stdout log XYslice* YZslice* XZslice* Xline* Yline* Zline* \
+      Diag* BaseXYslice* BaseYZslice* BaseXZslice* BaseXline* BaseYline* BaseZline* BaseDiag* \
+      PowerSpec_* Particle_* nohup.out Record__Performance Record__TimingMPI_* \
+      Record__ParticleCount Record__User Patch_* Record__NCorrUnphy FailedPatchGroup* *.pyc Record__LoadBalance \
+      GRACKLE_INFO Record__DivB Record__GrackleComoving coolingrate_z*.dat

example/test_problem/Hydro/Grackle_Comoving/download_datatable.sh

@@ -0,0 +1,2 @@
+filename=CloudyData_noUVB.h5
+curl -L https://github.com/grackle-project/grackle_data_files/raw/main/input/CloudyData_noUVB.h5 -o $filename

example/test_problem/Hydro/Grackle_Comoving/generate_make.sh

@@ -0,0 +1,6 @@
+# This script should run in the same directory as configure.py
+
+PYTHON=python3
+
+${PYTHON} configure.py --machine=eureka_intel --hdf5=true --model=HYDRO --comoving=true --gravity=true --fftw=FFTW3 \
+                       --grackle=true --double=true --passive=1 "$@"

example/test_problem/Hydro/Grackle_Comoving/plot/Compute_CoolingCurve.py

@@ -0,0 +1,67 @@
+#!/bin/python
+
+# reconstruct the cooling rate from the time evolution of the temperature
+# and compare with the cooling rate table
+import numpy as np
+import matplotlib.pyplot as plt
+
+# Load the data
+data = np.loadtxt('../Record__GrackleComoving', skiprows=10)
+
+time          = data[:,0] # scale factor
+dt            = data[:,2] # sec
+nden          = data[:,3] # cm^-3
+mu            = data[:,4] # mean molecular weight
+temp          = data[:,5] # K
+eden          = data[:,6] # erg/cm^3
+lcool_grackle = data[:,7] # erg/cm^3/s (cooling rate computed by Grackle)
+
+# reconstruct the cooling rate from Record__GrackleComoving
+dedt = np.gradient(eden) / dt
+dadt = np.gradient(time) / dt
+
+dedt_com = eden * (-5 / time * dadt) # loss due to cosmic expansion
+dedt_rad = dedt - dedt_com           # loss due to radiation
+
+L_cool = (-dedt_rad) / nden**2
+
+fig, ax = plt.subplots(1, 2, figsize=(10, 5))
+# ax[0].plot(temp, lcool, 'r-')
+
+# cooling rate table
+for z in range(0, 10):
+    cr = np.loadtxt('../coolingrate_z{:.1f}.dat'.format(z))
+    ax[0].plot(cr[:,0], cr[:,1], '--', lw=0.5, label='z={}'.format(z))
+
+# cooling rate computed by Grackle during the simulation
+ax[0].plot(temp, lcool_grackle, label='Grackle', ms=0, ls='dashed', lw=1, color='gray')
+
+# result from GAMER
+ax[0].scatter(temp[-1], L_cool[-1], s=3, label='GAMER (reconstructed) (z=0)')
+for z in range(1, 9):
+    a = 1 / (1 + z)
+    idx = np.where(time <= a)[0][-1]
+    a0 = time[idx]
+    a1 = time[idx+1]
+    temp_a = (temp[idx] * (a1 - a) + temp[idx+1] * (a - a0)) / (a1 - a0)
+    L_cool_a = (L_cool[idx] * (a1 - a) + L_cool[idx+1] * (a - a0)) / (a1 - a0)
+    ax[0].scatter(temp_a, L_cool_a, s=3, label='GAMER (reconstructed) (z={})'.format(z), zorder=10)
+ax[0].scatter(temp[0], L_cool[0], s=3, label='GAMER (reconstructed) (z=9)', zorder=10)
+
+ax[0].plot(temp, L_cool, label='GAMER (reconstructed)', ms=0, ls='solid', lw=0.5, color='black')
+
+ax[0].set_xlabel('Temperature [K]')
+ax[0].set_ylabel(r'$\Lambda$ [erg cm$^3$ s$^{-1}$]')
+ax[0].set_xscale('log')
+ax[0].set_yscale('log')
+ax[0].set_xlim(1e4, 2e8)
+ax[0].set_ylim(3e-25, 2e-21)
+ax[0].legend(ncol=2, fontsize=6)
+
+ax[1].plot(time, temp, 'r-')
+ax[1].set_xlabel('a')
+ax[1].set_ylabel('Temperature [K]')
+ax[1].set_yscale('log')
+
+
+plt.savefig('CoolingCurve.png', bbox_inches='tight', dpi=300)

example/test_problem/Template/Input__Parameter

@@ -49,6 +49,7 @@ TESTPROB_ID                   0           # test problem ID [0]
                                           #   23: HYDRO MHD Cosmic Ray Diffusion
                                           #  100: HYDRO CDM cosmological simulation (+GRAVITY & COMOVING & PARTICLE)
                                           #  101: HYDRO Zeldovich pancake collapse (+GRAVITY & COMOVING & PARTICLE)
+                                          #  102: HYDRO Grackle comoving (+GRAVITY & COMOVING & GRACKLE )
                                           # 1000: ELBDM external potential (+GRAVITY)
                                           # 1001: ELBDM Jeans instability in the comoving frame (+GRAVITY, +COMOVING)
                                           # 1002: ELBDM Jeans instability in the physical frame (+GRAVITY)

include/GAMER.h

@@ -58,6 +58,8 @@
 #  define CONFIG_BFLOAT_4
 #endif
 
+#define OMIT_LEGACY_INTERNAL_GRACKLE_FUNC
+
 extern "C" {
 #  include <grackle.h>
 }

include/Typedef.h

@@ -80,6 +80,7 @@ const TestProbID_t
    TESTPROB_HYDRO_JET_ICM_WALL                 =   52,
    TESTPROB_HYDRO_CDM_LSS                      =  100,
    TESTPROB_HYDRO_ZELDOVICH                    =  101,
+   TESTPROB_HYDRO_GRACKLE_COMOVING             =  102,
    TESTPROB_ELBDM_EXTPOT                       = 1000,
    TESTPROB_ELBDM_JEANS_INSTABILITY_COMOVING   = 1001,
    TESTPROB_ELBDM_JEANS_INSTABILITY_PHYSICAL   = 1002,

src/Grackle/Grackle_AdvanceDt.cpp

@@ -26,8 +26,15 @@
 void Grackle_AdvanceDt( const int lv, const double TimeNew, const double TimeOld, const double dt, const int SaveSg,
                         const bool OverlapMPI, const bool Overlap_Sync )
 {
-
-   InvokeSolver( GRACKLE_SOLVER, lv, TimeNew, TimeOld, dt, NULL_REAL, SaveSg, NULL_INT, NULL_INT, OverlapMPI, Overlap_Sync );
+#  ifdef COMOVING
+// convert dt from the comoving time interval to the physical time interval
+// --> see Equation (15) of Schive, Tsai, & Chiueh (2010)
+   const double dt_grackle = dt * SQR(TimeOld);
+#  else
+   const double dt_grackle = dt;
+#  endif
+
+   InvokeSolver( GRACKLE_SOLVER, lv, TimeNew, TimeOld, dt_grackle, NULL_REAL, SaveSg, NULL_INT, NULL_INT, OverlapMPI, Overlap_Sync );
 
 } // FUNCTION : Grackle_AdvanceDt
 

src/Grackle/Grackle_Close.cpp

@@ -117,7 +117,12 @@ void Grackle_Close( const int lv, const int SaveSg, const real h_Che_Array[], co
          {
 //          apply internal energy floor
             Dens = Ptr_Dens [idx_pg];
-            Eint = Ptr_sEint[idx_pg]*Dens;
+#           ifdef COMOVING
+//          convert from the proper frame specific internal energy to the comoving internal energy density
+            Eint = Ptr_sEint[idx_pg] * Dens * SQR(Time[lv]);
+#           else
+            Eint = Ptr_sEint[idx_pg] * Dens;
+#           endif
             Eint = Hydro_CheckMinEint( Eint, MIN_EINT );
 
 //          update the total energy density

src/Grackle/Grackle_Init.cpp

@@ -11,7 +11,6 @@
 // Description :  Initialize the chemistry and radiative cooling library Grackle
 //
 // Note        :  1. Must be called AFTER Init_Load_Parameter(), Init_Unit(), and Init_OpenMP()
-//                2. COMOVING is not supported yet
 //
 // Parameter   :  None
 //
@@ -35,10 +34,6 @@ void Grackle_Init()
                  sizeof(real), sizeof(gr_float) );
                  */
 
-// comoving frame is not supported yet
-#  ifdef COMOVING
-   Aux_Error( ERROR_INFO, "SUPPORT_GRACKLE does not work with COMOVING yet !!\n" );
-#  endif
 
    if (  ( GRACKLE_PRIMORDIAL == GRACKLE_PRI_CHE_CLOUDY || GRACKLE_METAL || GRACKLE_UV )  &&
          !Aux_CheckFileExist(GRACKLE_CLOUDY_TABLE)  )
@@ -52,11 +47,16 @@ void Grackle_Init()
 // units in cgs
 // --> Che_Units is declared as a global variable since all Grackle solvers require that as well
 #  ifdef COMOVING
+// 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 / (a_scale * time_units)
+// --> the specific energy passed to Grackle is defined on proper coordinates, i.e., u = k_B * T / ((gamma - 1) * mu * m_p)
    Che_Units.comoving_coordinates = 1;
-   Che_Units.density_units        = NULL_REAL;  // not sure how to set the units in the comoving coordinates yet...
-   Che_Units.length_units         = NULL_REAL;  // --> see http://grackle.readthedocs.io/en/latest/Integration.html
-   Che_Units.time_units           = NULL_REAL;
-   Che_Units.velocity_units       = NULL_REAL;
+   Che_Units.density_units        = UNIT_D / CUBE(Time[0]);
+   Che_Units.length_units         = UNIT_L * Time[0];
+   Che_Units.time_units           = UNIT_T;
+   Che_Units.velocity_units       = UNIT_V;
    Che_Units.a_units              = 1.0;
    Che_Units.a_value              = Time[0];
 

src/Grackle/Grackle_Prepare.cpp

@@ -34,6 +34,7 @@ extern int CheIdx_Metal;
 //                   --> Che_NField and the corresponding array indices in h_Che_Array[] (e.g., CheIdx_Dens)
 //                       are declared and set by Init_MemAllocate_Grackle()
 //                2. This function always prepares the latest FluSg data
+//                3. If COMOVING is on, update grackle units
 //
 // Parameter   :  lv          : Target refinement level
 //                h_Che_Array : Host array to store the prepared data
@@ -91,14 +92,24 @@ void Grackle_Prepare( const int lv, real h_Che_Array[], const int NPG, const int
    }
 #  endif // #ifdef GAMER_DEBUG
 
+#  ifdef COMOVING
+// update grackle units
+   Che_Units.comoving_coordinates = 1;
+   Che_Units.density_units        = UNIT_D / CUBE(Time[lv]);
+   Che_Units.length_units         = UNIT_L * Time[lv];
+   Che_Units.time_units           = UNIT_T;
+   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;
-   const real MassRatio_pe    = Const_mp / Const_me;
+   const int  Size1pg             = CUBE(PS2);
+   const int  Size1v              = NPG*Size1pg;
+   const real MassRatio_pe        = Const_mp / Const_me;
 #  ifdef DUAL_ENERGY
-   const bool CheckMinPres_No  = false;
+   const bool CheckMinPres_No     = false;
 #  else
-   const bool CheckMinEint_Yes = true;
+   const bool CheckMinEint_Yes    = true;
 #  endif
 
    real *Ptr_Dens0  = h_Che_Array + CheIdx_Dens *Size1v;
@@ -203,9 +214,15 @@ void Grackle_Prepare( const int lv, real h_Che_Array[], const int NPG, const int
 
 //          mandatory fields
             Ptr_Dens [idx_pg] = Dens;
-            Ptr_sEint[idx_pg] = Eint / Dens;
             Ptr_Ent  [idx_pg] = Etot - Eint; // non-thermal energy density
 
+#           ifdef COMOVING
+//          convert from the comoving specific internal energy to the proper frame
+            Ptr_sEint[idx_pg] = Eint / Dens / SQR(Time[lv]);
+#           else
+            Ptr_sEint[idx_pg] = Eint / Dens;
+#           endif
+
 //          6-species network
             if ( GRACKLE_PRIMORDIAL >= GRACKLE_PRI_CHE_NSPE6 ) {
             Ptr_e    [idx_pg] = *( fluid[Idx_e    ][0][0] + idx_p ) * MassRatio_pe;

src/Init/Init_Field.cpp

@@ -100,8 +100,10 @@ void Init_Field()
 
 // 2. add other predefined fields
 #  ifdef SUPPORT_GRACKLE
+// The masses of electron and deuterium are ignored in normalization in Grackle
+// --> see the subroutine make_consistent_g in solve_rate_cool_g.F in the Grackle source code
    if ( GRACKLE_PRIMORDIAL >= GRACKLE_PRI_CHE_NSPE6 ) {
-   Idx_e       = AddField( "Electron", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_YES, INTERP_FRAC_YES );
+   Idx_e       = AddField( "Electron", FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_NO,  INTERP_FRAC_YES ); // electron mass is neglected in Grackle
    Idx_HI      = AddField( "HI",       FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_YES, INTERP_FRAC_YES );
    Idx_HII     = AddField( "HII",      FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_YES, INTERP_FRAC_YES );
    Idx_HeI     = AddField( "HeI",      FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_YES, INTERP_FRAC_YES );
@@ -116,9 +118,9 @@ void Init_Field()
    }
 
    if ( GRACKLE_PRIMORDIAL >= GRACKLE_PRI_CHE_NSPE12 ) {
-   Idx_DI      = AddField( "DI",       FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_YES, INTERP_FRAC_YES );
-   Idx_DII     = AddField( "DII",      FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_YES, INTERP_FRAC_YES );
-   Idx_HDI     = AddField( "HDI",      FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_YES, INTERP_FRAC_YES );
+   Idx_DI      = AddField( "DI",       FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_NO,  INTERP_FRAC_YES ); // deuterium mass is neglected in Grackle
+   Idx_DII     = AddField( "DII",      FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_NO,  INTERP_FRAC_YES ); // deuterium mass is neglected in Grackle
+   Idx_HDI     = AddField( "HDI",      FIXUP_FLUX_YES, FIXUP_REST_YES, NORMALIZE_NO,  INTERP_FRAC_YES ); // deuterium mass is neglected in Grackle
    }
 
 // normalize the metallicity field only when adopting the non-equilibrium chemistry

src/Init/Init_TestProb.cpp

@@ -26,6 +26,7 @@ void Init_TestProb_Hydro_BarredPot();
 void Init_TestProb_Hydro_ParticleTest();
 void Init_TestProb_Hydro_CDM_LSS();
 void Init_TestProb_Hydro_Zeldovich();
+void Init_TestProb_Hydro_Grackle_Comoving();
 void Init_TestProb_Hydro_EnergyPowerSpectrum();
 void Init_TestProb_Hydro_CR_SoundWave();
 void Init_TestProb_Hydro_CR_ShockTube();
@@ -92,6 +93,7 @@ void Init_TestProb()
       case TESTPROB_HYDRO_PARTICLE_TEST :                Init_TestProb_Hydro_ParticleTest();                break;
       case TESTPROB_HYDRO_CDM_LSS :                      Init_TestProb_Hydro_CDM_LSS();                     break;
       case TESTPROB_HYDRO_ZELDOVICH :                    Init_TestProb_Hydro_Zeldovich();                   break;
+      case TESTPROB_HYDRO_GRACKLE_COMOVING :             Init_TestProb_Hydro_Grackle_Comoving();            break;
       case TESTPROB_HYDRO_ENERGY_POWER_SPECTRUM :        Init_TestProb_Hydro_EnergyPowerSpectrum();         break;
       case TESTPROB_HYDRO_CR_SOUNDWAVE :                 Init_TestProb_Hydro_CR_SoundWave();                break;
       case TESTPROB_HYDRO_CR_SHOCKTUBE :                 Init_TestProb_Hydro_CR_ShockTube();                break;