WIP: Automatic Unit Handling

src/clib/CMakeLists.txt

@@ -129,6 +129,7 @@ add_library(Grackle_Grackle
   set_default_chemistry_parameters.c
   solve_chemistry.c
   update_UVbackground_rates.c
+  unit_handling.h
   utils.c
 
   # auto-generated C source files

src/clib/calculate_cooling_time.c

@@ -14,19 +14,17 @@
 #include <stdlib.h>
 #include <stdio.h>
 #include <math.h>
+#include "grackle.h"
 #include "grackle_macros.h"
-#include "grackle_types.h"
-#include "grackle_chemistry_data.h"
 #include "phys_constants.h"
+#include "unit_handling.h"
 #include "utils.h"
 
 extern chemistry_data *grackle_data;
 extern chemistry_data_storage grackle_rates;
 
 /* function prototypes */
 
-double get_temperature_units(code_units *my_units);
-
 int update_UVbackground_rates(chemistry_data *my_chemistry,
                               chemistry_data_storage *my_rates,
                               photo_rate_storage *my_uvb_rates,
@@ -90,6 +88,17 @@ int local_calculate_cooling_time(chemistry_data *my_chemistry,
   if (!my_chemistry->use_grackle)
     return SUCCESS;
 
+  /* do unit-handling */
+  code_units units = determine_code_units(my_units, my_rates,
+                                          my_fields->current_a_value,
+                                          my_chemistry->unit_handling,
+                                          "local_calculate_cooling_time");
+  if (units.a_units < 0) {
+    return FAIL;
+  } else {
+    my_units = &units;
+  }
+
   /* Update UV background rates. */
   photo_rate_storage my_uvb_rates;
 

src/clib/calculate_dust_temperature.c

@@ -14,10 +14,10 @@
 #include <stdlib.h>
 #include <stdio.h>
 #include <math.h>
+#include "grackle.h"
 #include "grackle_macros.h"
-#include "grackle_types.h"
-#include "grackle_chemistry_data.h"
 #include "phys_constants.h"
+#include "unit_handling.h"
 #ifdef _OPENMP
 #include <omp.h>
 #endif
@@ -27,8 +27,6 @@ extern chemistry_data_storage grackle_rates;
 
 /* function prototypes */
 
-double get_temperature_units(code_units *my_units);
-
 int local_calculate_temperature(chemistry_data *my_chemistry,
                                 chemistry_data_storage *my_rates,
                                 code_units *my_units,
@@ -65,6 +63,17 @@ int local_calculate_dust_temperature(chemistry_data *my_chemistry,
   if (my_chemistry->dust_chemistry < 1 && my_chemistry->h2_on_dust < 1)
     return SUCCESS;
 
+  /* do unit-handling */
+  code_units units = determine_code_units(my_units, my_rates,
+                                          my_fields->current_a_value,
+                                          my_chemistry->unit_handling,
+                                          "local_calculate_dust_temperature");
+  if (units.a_units < 0) {
+    return FAIL;
+  } else {
+    my_units = &units;
+  }
+
   double co_length_units, co_density_units;
   if (my_units->comoving_coordinates == TRUE) {
     co_length_units = my_units->length_units;

src/clib/calculate_gamma.c

@@ -14,11 +14,11 @@
 #include <stdlib.h>
 #include <stdio.h>
 #include <math.h>
+#include "grackle.h"
 #include "grackle_macros.h"
-#include "grackle_types.h"
-#include "grackle_chemistry_data.h"
 #include "phys_constants.h"
 #include "index_helper.h"
+#include "unit_handling.h"
 #ifdef _OPENMP
 #include <omp.h>
 #endif
@@ -43,6 +43,18 @@ int local_calculate_gamma(chemistry_data *my_chemistry,
 
   if (!my_chemistry->use_grackle)
     return SUCCESS;
+
+  /* do unit-handling */
+  code_units units = determine_code_units(my_units, my_rates,
+                                          my_fields->current_a_value,
+                                          my_chemistry->unit_handling,
+                                          "local_calculate_gamma");
+  if (units.a_units < 0) {
+    return FAIL;
+  } else {
+    my_units = &units;
+  }
+
 
   const grackle_index_helper ind_helper = _build_index_helper(my_fields);
   int outer_ind, index;

src/clib/calculate_pressure.c

@@ -14,20 +14,18 @@
 #include <stdlib.h>
 #include <stdio.h>
 #include <math.h>
+#include "grackle.h"
 #include "grackle_macros.h"
-#include "grackle_types.h"
-#include "grackle_chemistry_data.h"
 #include "phys_constants.h"
 #include "index_helper.h"
+#include "unit_handling.h"
 #ifdef _OPENMP
 #include <omp.h>
 #endif
 
 extern chemistry_data *grackle_data;
 extern chemistry_data_storage grackle_rates;
 
-double get_temperature_units(code_units *my_units);
-
 int local_calculate_pressure(chemistry_data *my_chemistry,
                              chemistry_data_storage *my_rates,
                              code_units *my_units,
@@ -38,6 +36,18 @@ int local_calculate_pressure(chemistry_data *my_chemistry,
   if (!my_chemistry->use_grackle)
     return SUCCESS;
 
+  /* do unit-handling */
+  code_units units = determine_code_units(my_units, my_rates,
+                                          my_fields->current_a_value,
+                                          my_chemistry->unit_handling,
+                                          "calculate_pressure");
+  if (units.a_units < 0) {
+    return FAIL;
+  } else {
+    my_units = &units;
+  }
+
+
   double tiny_number = 1.e-20;
   const grackle_index_helper ind_helper = _build_index_helper(my_fields);
   int outer_ind, index;

src/clib/calculate_temperature.c

@@ -14,11 +14,11 @@
 #include <stdlib.h>
 #include <stdio.h>
 #include <math.h>
+#include "grackle.h"
 #include "grackle_macros.h"
-#include "grackle_types.h"
-#include "grackle_chemistry_data.h"
 #include "phys_constants.h"
 #include "index_helper.h"
+#include "unit_handling.h"
 #ifdef _OPENMP
 #include <omp.h>
 #endif
@@ -47,8 +47,6 @@ extern void FORTRAN_NAME(calc_temp_cloudy_g)(
         double *priPar1, double *priPar2, double *priPar3, 
  	long long *priDataSize, double *priMMW);
 
-double get_temperature_units(code_units *my_units);
-
 int local_calculate_pressure(chemistry_data *my_chemistry,
                              chemistry_data_storage *my_rates,
                              code_units *my_units,
@@ -71,6 +69,17 @@ int local_calculate_temperature(chemistry_data *my_chemistry,
   if (!my_chemistry->use_grackle)
     return SUCCESS;
 
+  /* do unit-handling */
+  code_units units = determine_code_units(my_units, my_rates,
+                                          my_fields->current_a_value,
+                                          my_chemistry->unit_handling,
+                                          "local_solve_chemistry");
+  if (units.a_units < 0) {
+    return FAIL;
+  } else {
+    my_units = &units;
+  }
+
   /* Compute the pressure first. */
 
   if (my_chemistry->primordial_chemistry > 0) {

src/clib/grackle_chemistry_data_fields.def

@@ -208,6 +208,10 @@ ENTRY(collisional_ionisation_rates, INT, 1) //Collisional ionisation
 ENTRY(recombination_cooling_rates, INT, 1) //Recombination cooling
 ENTRY(bremsstrahlung_cooling_rates, INT, 1) //Bremsstrahlung cooling
 
+/* Flag to specify unit-handling. Use the GR_UNIT_HANDLING_LEGACY and
+   GR_UNIT_HANDLING_AUTOMATIC to set the flag values */
+ENTRY(unit_handling, INT, -3)
+
 /* maximum number of subcycle iterations for solve_chemistry */
 ENTRY(max_iterations, INT, 10000)
 

src/clib/initialize_chemistry_data.c

@@ -236,6 +236,21 @@ int local_initialize_chemistry_data(chemistry_data *my_chemistry,
 //omp_set_schedule( omp_sched_auto,    chunk_size );
 # endif
 
+  if (my_chemistry->unit_handling == -3) {
+    if (grackle_verbose) {
+      fprintf(stderr, ("WARNING: unit_handling is unset. Defaulting to legacy "
+                       "handling.\n"));
+    }
+    my_chemistry->unit_handling = GR_UNIT_HANDLING_LEGACY;
+  } else if ((my_chemistry->unit_handling != GR_UNIT_HANDLING_LEGACY) &&
+             (my_chemistry->unit_handling != GR_UNIT_HANDLING_AUTOMATIC)) {
+    fprintf(stderr, "unit_handling has an invalid value\n");
+    return GR_FAIL;
+  }
+
+  /* store a copy of the initial units */
+  my_rates->initial_units = *my_units;
+
   /* Only allow a units to be one with proper coordinates. */
   if (my_units->comoving_coordinates == FALSE &&
       my_units->a_units != 1.0) {
@@ -342,8 +357,6 @@ int local_initialize_chemistry_data(chemistry_data *my_chemistry,
     return GR_FAIL;
   }
 
-  /* store a copy of the initial units */
-  my_rates->initial_units = *my_units;
 
   if (grackle_verbose) {
     time_t timer;

src/clib/set_default_chemistry_parameters.c

@@ -16,8 +16,7 @@
 #include <math.h>
 
 #include "grackle_macros.h"
-#include "grackle_types.h"
-#include "grackle_chemistry_data.h"
+#include "grackle.h"
 
 int grackle_verbose = FALSE;
 
@@ -53,6 +52,7 @@ int gr_initialize_field_data(grackle_field_data *my_fields)
   my_fields->grid_dimension = NULL;
   my_fields->grid_start = NULL;
   my_fields->grid_end = NULL;
+  my_fields->current_a_value = GR_SPECIFY_INITIAL_A_VALUE;
   my_fields->grid_dx = -1.0;
 
   // now, modify all members holding datafields to have values of NULL

src/clib/solve_chemistry.c

@@ -13,19 +13,17 @@
 
 #include <stdio.h>
 #include <math.h>
+#include "grackle.h"
 #include "grackle_macros.h"
-#include "grackle_types.h"
-#include "grackle_chemistry_data.h"
 #include "phys_constants.h"
+#include "unit_handling.h"
 #include "utils.h"
 
 extern chemistry_data *grackle_data;
 extern chemistry_data_storage grackle_rates;
 
 /* function prototypes */
 
-double get_temperature_units(code_units *my_units);
-
 int update_UVbackground_rates(chemistry_data *my_chemistry,
                               chemistry_data_storage *my_rates,
                               photo_rate_storage *my_uvb_rates,
@@ -102,6 +100,17 @@ int local_solve_chemistry(chemistry_data *my_chemistry,
   if (!my_chemistry->use_grackle)
     return SUCCESS;
 
+  /* do unit-handling */
+  code_units units = determine_code_units(my_units, my_rates,
+                                          my_fields->current_a_value,
+                                          my_chemistry->unit_handling,
+                                          "local_solve_chemistry");
+  if (units.a_units < 0) {
+    return FAIL;
+  } else {
+    my_units = &units;
+  }
+
   /* Update UV background rates. */
   photo_rate_storage my_uvb_rates;