Merge pull request #94 from leonardromano/master

Added options for H2 and dust physics

doc/source/Integration.rst

@@ -480,6 +480,12 @@ and pointers to all field arrays.
    :c:data:`H2_self_shielding` is set to 2.  Field data
    should be in :c:data:`length_units`.
 
+.. c:var:: gr_float *H2_custom_shielding_factor
+
+   Pointer to a field containing attenuation factors to 
+   be multiplied with the H\ :sub:`2`\ dissociation rate.
+   Used when the :c:data:`H2_custom_shielding` flag is set.
+
 .. c:var:: gr_float *isrf_habing
 
    Pointer to a field containing values of the strength

doc/source/Parameters.rst

@@ -210,6 +210,14 @@ For all on/off integer flags, 0 is off and 1 is on.
       directly from equation 2 of `Wolfire et al. (1995)
       <https://ui.adsabs.harvard.edu/abs/1995ApJ...443..152W/abstract>`__.
 
+.. c:var:: int dust_recombination_cooling
+
+   Flag to enable recombination cooling onto dust grains using 
+   equation 9 of `Wolfire et al. (1995) 
+   <https://ui.adsabs.harvard.edu/abs/1995ApJ...443..152W/abstract>`__
+   rescaled by the local dust-to-gas ratio. This option is automatically 
+   set by :c:data:`h2_on_dust` > 0 or :c:data:`dust_chemistry` > 0.
+   Default: 0.
 
 .. note:: With :c:data:`primordial_chemistry` > 0, the electron density
    used to calculate epsilon for :c:data:`photoelectric_heating` = 3
@@ -367,6 +375,18 @@ For all on/off integer flags, 0 is off and 1 is on.
       field.
     - 3: Use the local Jeans length.
 
+.. c:var:: int H2_custom_shielding
+
+   Flag to enable the user to provide an additional field which acts as 
+   an additional attenuation factor for both the UV background dissociation 
+   rate and the H\ :sub:`2`\  dissociation rate given by 
+   :c:data:`RT_H2_dissociation_rate` (if present), that is separate from the 
+   :c:data:`H2_self_shielding` attenuation factor. 
+   The factor, which is intended to be unspecific can e.g. be used in order 
+   to include grain size dependent dust extinction or any other user-specific 
+   source of attenuation.
+   Default: 0.
+
 .. c:var:: int self_shielding_method
 
    Switch to enable approximate self-shielding from the UV background.

src/clib/calculate_cooling_time.c

@@ -35,7 +35,7 @@ extern void FORTRAN_NAME(cool_multi_time_g)(
 	gr_float *cooltime,
 	int *in, int *jn, int *kn, int *nratec, int *iexpand,
         int *ispecies, int *imetal, int *imcool, int *idust, int *idustall,
-        int *idustfield, int *idim,
+        int *idustfield, int *idustrec, int *idim,
 	int *is, int *js, int *ks, int *ie, int *je, int *ke, 
         int *ih2co, int *ipiht, int *igammah,
 	double *aye, double *temstart, double *temend,
@@ -170,6 +170,7 @@ int local_calculate_cooling_time(chemistry_data *my_chemistry,
        &my_chemistry->h2_on_dust,
        &my_chemistry->dust_chemistry,
        &my_chemistry->use_dust_density_field,
+       &my_chemistry->dust_recombination_cooling,
        &(my_fields->grid_rank),
        my_fields->grid_start,
        my_fields->grid_start+1,

src/clib/cool1d_multi_g.F

@@ -7,7 +7,7 @@ subroutine cool1d_multi_g(
      &                d, e, u, v, w, de, HI, HII, HeI, HeII, HeIII,
      &                in, jn, kn, nratec, 
      &                iexpand, ispecies, imetal, imcool,
-     &                idust, idustall, idustfield,
+     &                idust, idustall, idustfield, idustrec,
      &                idim, is, ie, j, k, ih2co, ipiht, iter, igammah,
      &                aye, temstart, temend, z_solar, fgr,
      &                utem, uxyz, uaye, urho, utim,
@@ -72,7 +72,7 @@ subroutine cool1d_multi_g(
 
       integer in, jn, kn, is, ie, j, k, nratec, idim,
      &        iexpand, ih2co, ipiht, ispecies, imcool,
-     &        idust, idustall, idustfield,
+     &        idust, idustall, idustfield, idustrec,
      &        imetal, igammah, ih2optical, iciecool, clnew,
      &        iVheat, iMheat, iradtrans, iradshield,
      &        iisrffield
@@ -167,7 +167,7 @@ subroutine cool1d_multi_g(
 
 !     Set flag for dust-related options
 
-      anydust = (idust .gt. 0) .or. (idustall .gt. 0)
+      anydust = (idust .gt. 0) .or. (idustall .gt. 0) .or. (idustrec .gt. 0)
 
 !     Set flag for needing interpolation variables
 
@@ -1001,7 +1001,7 @@ subroutine cool1d_multi_g(
 
 !     Electron recombination onto dust grains (eqn. 9 of Wolfire 1995)
 
-      if (idustall .gt. 0) then
+      if ((idustall .gt. 0) .or. (idustrec .gt. 0)) then
 
           do i = is + 1, ie + 1
              if (itmask(i)) then
@@ -1015,7 +1015,7 @@ subroutine cool1d_multi_g(
                 grbeta = 0.74_DKIND / tgas(i)**0.068_DKIND
                 edot(i) = edot(i) -
      &               regr(i) * (myisrf(i)*dom_inv / myde(i))**grbeta *
-     &               myde(i) * rhoH(i)
+     &               myde(i) * rhoH(i) * dust2gas(i) / fgr
              endif
           enddo
 

src/clib/cool_multi_time_g.F

@@ -6,7 +6,7 @@ subroutine cool_multi_time_g(
      &                  cooltime,
      &                in, jn, kn, nratec, iexpand, 
      &                ispecies, imetal, imcool, idust,
-     &                idustall, idustfield, idim,
+     &                idustall, idustfield, idustrec, idim,
      &                is, js, ks, ie, je, ke, ih2co, ipiht, igammah,
      &                aye, temstart, temend,
      &                utem, uxyz, uaye, urho, utim,
@@ -63,7 +63,7 @@ subroutine cool_multi_time_g(
       integer in, jn, kn, is, js, ks, ie, je, ke, nratec, 
      &        iexpand, ih2co, ipiht, ispecies, imetal, idim,
      &        imcool, idust, idustall,
-     &        idustfield, igammah, ih2optical, iciecool,
+     &        idustfield, idustrec, igammah, ih2optical, iciecool,
      &        clnew, iVheat, iMheat, iradtrans, iradshield,
      &        iisrffield
 
@@ -253,7 +253,7 @@ subroutine cool_multi_time_g(
      &                d, e, u, v, w, de, HI, HII, HeI, HeII, HeIII,
      &                in, jn, kn, nratec, 
      &                iexpand, ispecies, imetal, imcool,
-     &                idust, idustall, idustfield,
+     &                idust, idustall, idustfield, idustrec,
      &                idim, is, ie, j, k, ih2co, ipiht, 1, igammah,
      &                aye, temstart, temend, z_solar, fgr,
      &                utem, uxyz, uaye, urho, utim,

src/clib/grackle_chemistry_data.h

@@ -70,6 +70,9 @@ typedef struct
   /* Flag to supply a dust density field */
   int use_dust_density_field;
 
+  /* Flag for enabling recombination cooling on grains */
+  int dust_recombination_cooling;
+
   /* photo-electric heating from irradiated dust */
   int photoelectric_heating;
   double photoelectric_heating_rate;
@@ -138,6 +141,10 @@ typedef struct
      Please refer to the grackle documentation for specifics. */
   int H2_self_shielding;
 
+  /* flag for custom H2-shielding factor. The factor is provided as an additional field 
+     by the user and is multiplied to the rate for radiative H2 dissocitation */
+  int H2_custom_shielding;
+
   /* flag to select which formula for calculating k11 you want to use. 
      Setting to 1 will use Savin 2004, 2 will use Abel et al. 1996  */
   int h2_charge_exchange_rate;

src/clib/grackle_fortran_interface.def

@@ -91,6 +91,7 @@ c     This is the fortran definition of grackle_chemistry_data
          REAL(C_DOUBLE) :: Gamma
          INTEGER(C_INT) :: h2_on_dust
          INTEGER(C_INT) :: use_dust_density_field
+         INTEGER(C_INT) :: dust_recombination_cooling
          INTEGER(C_INT) :: photoelectric_heating
          REAL(C_DOUBLE) :: photoelectric_heating_rate
          INTEGER(C_INT) :: use_isrf_field
@@ -127,6 +128,7 @@ c     This is the fortran definition of grackle_chemistry_data
          INTEGER(C_INT) :: radiative_transfer_hydrogen_only
          INTEGER(C_INT) :: self_shielding_method
          INTEGER(C_INT) :: H2_self_shielding
+         INTEGER(C_INT) :: H2_custom_shielding
          INTEGER(C_INT) :: h2_charge_exchange_rate
          INTEGER(C_INT) :: h2_dust_rate
          INTEGER(C_INT) :: h2_h_cooling_rate

src/clib/grackle_types.h

@@ -68,6 +68,7 @@ typedef struct
   gr_float *RT_H2_dissociation_rate;
 
   gr_float *H2_self_shielding_length;
+  gr_float *H2_custom_shielding_factor;
 
   gr_float *isrf_habing;
 

src/clib/initialize_chemistry_data.c

@@ -73,6 +73,13 @@ int _initialize_chemistry_data(chemistry_data *my_chemistry,
       }
     }
 
+    if (my_chemistry->dust_recombination_cooling < 0) {
+      my_chemistry->dust_recombination_cooling = 1;
+      if (grackle_verbose) {
+        fprintf(stdout, "Dust chemistry enabled, setting dust_recombination_cooling to 1.\n");
+      }
+    }
+
     if (my_chemistry->primordial_chemistry > 1 &&
         my_chemistry->h2_on_dust == 0) {
       my_chemistry->h2_on_dust = 1;
@@ -119,7 +126,7 @@ int _initialize_chemistry_data(chemistry_data *my_chemistry,
     my_chemistry->HydrogenFractionByMass = 1. / (1. + 0.1 * 3.971);
   }
 
-  if (my_chemistry->h2_on_dust > 0 || my_chemistry->dust_chemistry > 0) {
+  if (my_chemistry->h2_on_dust > 0 || my_chemistry->dust_chemistry > 0 || my_chemistry->dust_recombination_cooling > 0) {
     my_rates->gas_grain = malloc(my_chemistry->NumberOfTemperatureBins * sizeof(double));
     my_rates->regr      = malloc(my_chemistry->NumberOfTemperatureBins * sizeof(double));
   }
@@ -265,6 +272,8 @@ void show_parameters(FILE *fp, chemistry_data *my_chemistry)
           my_chemistry->h2_on_dust);
   fprintf(fp, "use_dust_density_field            = %d\n",
           my_chemistry->use_dust_density_field);
+  fprintf(fp, "dust_recombination_cooling        = %d\n",
+          my_chemistry->dust_recombination_cooling);
   fprintf(fp, "photoelectric_heating             = %d\n",
           my_chemistry->photoelectric_heating);
   fprintf(fp, "photoelectric_heating_rate        = %g\n",
@@ -335,6 +344,8 @@ void show_parameters(FILE *fp, chemistry_data *my_chemistry)
           my_chemistry->radiative_transfer_hydrogen_only);
   fprintf(fp, "self_shielding_method             = %d\n",
           my_chemistry->self_shielding_method);
+  fprintf(fp, "H2_custom_shielding               = %d\n",
+          my_chemistry->H2_custom_shielding);
   fprintf(fp, "H2_self_shielding                 = %d\n",
           my_chemistry->H2_self_shielding);
 # ifdef _OPENMP

src/clib/initialize_rates.c

@@ -215,7 +215,7 @@ int initialize_rates(chemistry_data *my_chemistry, chemistry_data_storage *my_ra
 { 
     //* Set the flag for dust calculations.
     int anyDust;
-    if ( my_chemistry->h2_on_dust > 0 || my_chemistry->dust_chemistry > 0) {
+    if ( my_chemistry->h2_on_dust > 0 || my_chemistry->dust_chemistry > 0 || my_chemistry->dust_recombination_cooling > 0) {
         anyDust = TRUE;
     } else {
         anyDust = FALSE;

src/clib/set_default_chemistry_parameters.c

@@ -69,6 +69,7 @@ chemistry_data _set_default_chemistry_parameters(void)
   my_chemistry.recombination_cooling_rates    = 1;
   my_chemistry.bremsstrahlung_cooling_rates   = 1;
 
+  my_chemistry.dust_recombination_cooling     = -1; // unset
   my_chemistry.photoelectric_heating          = -1; // unset
   // epsilon=0.05, G_0=1.7 (in erg s^-1 cm^-3)
   my_chemistry.photoelectric_heating_rate     = 8.5e-26;
@@ -131,6 +132,7 @@ chemistry_data _set_default_chemistry_parameters(void)
   /* approximate self-shielding */
   my_chemistry.self_shielding_method                  = 0;
   my_chemistry.H2_self_shielding                      = 0;
+  my_chemistry.H2_custom_shielding                    = 0;
 
 //number of OpenMP threads
 # ifdef _OPENMP

src/clib/solve_chemistry.c

@@ -36,7 +36,7 @@ extern void FORTRAN_NAME(solve_rate_cool_g)(
         int *ispecies, int *imetal, int *imcool, int *idust, int *idustall,
         int *idustfield, int *idim,
 	int *is, int *js, int *ks, int *ie, int *je, int *ke,
-        int *ih2co, int *ipiht, int *igammah,
+        int *ih2co, int *ipiht, int *idustrec, int *igammah,
 	double *dx, double *dt, double *aye, double *temstart, double *temend,
 	double *utem, double *uxyz, double *uaye, double *urho, double *utim,
 	double *gamma, double *fh, double *dtoh, double *z_solar, double *fgr,
@@ -82,7 +82,8 @@ extern void FORTRAN_NAME(solve_rate_cool_g)(
  	long long *metDataSize, double *metCooling,
         double *metHeating, int *clnew,
         int *iVheat, int *iMheat, gr_float *Vheat, gr_float *Mheat,
-        int *iisrffield, gr_float* isrf_habing);
+        int *iisrffield, gr_float* isrf_habing, 
+        int *iH2shieldcustom, gr_float* f_shield_custom);
 
 int local_solve_chemistry(chemistry_data *my_chemistry,
                           chemistry_data_storage *my_rates,
@@ -201,6 +202,7 @@ int local_solve_chemistry(chemistry_data *my_chemistry,
     my_fields->grid_end+2,
     &my_chemistry->ih2co,
     &my_chemistry->ipiht,
+    &my_chemistry->dust_recombination_cooling,
     &my_chemistry->photoelectric_heating,
     &(my_fields->grid_dx),
     &dt_value,
@@ -358,7 +360,9 @@ int local_solve_chemistry(chemistry_data *my_chemistry,
     my_fields->volumetric_heating_rate,
     my_fields->specific_heating_rate,
     &my_chemistry->use_isrf_field,
-    my_fields->isrf_habing);
+    my_fields->isrf_habing, 
+    &my_chemistry->H2_custom_shielding,
+    my_fields->H2_custom_shielding_factor);
 
   return SUCCESS;