Merge pull request #97 from brittonsmith/uvel

fix internal velocity units

.circleci/config.yml

@@ -18,7 +18,7 @@ commands:
       - run:
           name: "Install dependencies."
           command: |
-            git submodule update --init
+            git submodule update --init --remote
             source $BASH_ENV
             sudo apt update
             sudo apt install -y libhdf5-serial-dev gfortran libtool-bin

.gitmodules

@@ -1,3 +1,4 @@
 [submodule "grackle_data_files"]
 	path = grackle_data_files
 	url = https://github.com/grackle-project/grackle_data_files
+	branch = master

doc/source/Integration.rst

@@ -126,8 +126,9 @@ Code Units
 cosmological simulation.**  The :c:data:`code_units` structure contains
 conversions from code units to CGS.  If :c:data:`comoving_coordinates` is set to
 0, it is assumed that the fields passed into the solver are in the
-proper frame.  All of the units (density, length, time, velocity, and
-expansion factor) must be set.  When using the proper frame, :c:data:`a_units`
+proper frame. Units for length, time, and the expansion factor must be set
+manually. Units for velocity are then set by calling
+:c:data:`set_velocity_units`. When using the proper frame, :c:data:`a_units`
 (units for the expansion factor) must be set to 1.0.
 
 .. c:type:: code_units
@@ -158,11 +159,15 @@ expansion factor) must be set.  When using the proper frame, :c:data:`a_units`
 .. c:var:: double velocity_units
 
    Conversion factor to be multiplied by velocities to return proper cm/s.
+   This should be set units the :c:data:`set_velocity_units` function. Note,
+   units of specific energy (i.e., conversion to erg/g) are then defined
+   as :c:data:`velocity_units`\ :sup:`2` (velocity units squared).
 
 .. c:var:: double a_units
 
    Conversion factor to be multiplied by the expansion factor such that
-   a\ :sub:`true`\  = a\ :sub:`code`\ * :c:data:`a_units`.
+   a\ :sub:`true`\  = a\ :sub:`code`\ * :c:data:`a_units`. When using
+   proper coordinates, :c:data:`a_units` must be set to 1.
 
 .. c:var:: double a_value
 
@@ -181,23 +186,16 @@ expansion factor) must be set.  When using the proper frame, :c:data:`a_units`
   my_units.density_units = 1.67e-24; // 1 m_H/cc
   my_units.length_units = 3.086e21;  // 1 kpc
   my_units.time_units = 3.15569e13;  // 1 Myr
-  my_units.velocity_units = my_units.length_units / my_units.time_units;
   my_units.a_units = 1.0;            // units for the expansion factor
   my_units.a_value = 1. / (1. + current_redshift) / my_units.a_units;
+  // set velocity units
+  set_velocity_units(&my_units);
 
 If :c:data:`comoving_coordinates` is set to 1, it is assumed that the fields being 
-passed to the solver are in the comoving frame.  Hence, the units must 
+passed to the solver are in the comoving frame. Hence, the units must
 convert from code units in the **comoving** frame to CGS in the **proper** 
 frame.  
 
-.. note:: With :c:data:`comoving_coordinate` set to 1, velocity units need to be
-   defined in the following way.
-
-.. code-block:: c++
-
-  my_units.velocity_units = my_units.a_units * 
-    (my_units.length_units / a_value) / my_units.time_units; // since u = a * dx/dt
-
 For an example of using comoving units, see the units system in the 
 `Enzo <http://enzo-project.org/>`_ code.  For cosmological simulations, a
 comoving unit system is preferred, though not required, since it allows the 
@@ -409,7 +407,10 @@ and pointers to all field arrays.
 
 .. c:var:: gr_float* internal_energy
 
-   Pointer to the internal energy field array.
+   Pointer to the internal energy field array. Internal energies should be
+   in units of :c:data:`velocity_units`\ :sup:`2` (velocity units squared).
+   This can be converted to and from a temperature by using the
+   :c:data:`get_temperature_units` function.
 
 .. c:var:: gr_float* x_velocity
 

doc/source/Reference.rst

@@ -41,6 +41,43 @@ Primary Functions
    :rtype: int
    :returns: 1 (success) or 0 (failure)
 
+.. c:function:: void set_velocity_units(code_units *my_units);
+
+   Sets the :c:data:`velocity_units` value of the input ``my_units``
+   :c:data:`code_units` struct. For proper coordinates, velocity units are equal
+   to :c:data:`length_units` / :c:data:`time_units`. For comoving coordinates,
+   velocity units are equal to (:c:data:`length_units` / :c:data:`a_value`) /
+   :c:data:`time_units`.
+
+   :param code_units* my_units: code units conversions
+
+.. c:function:: double get_velocity_units(code_units *my_units);
+
+   Returns the appropriate value for velocity units given the values of
+   :c:data:`length_units`, :c:data:`a_value`, and :c:data:`time_units`
+   in the input ``my_units`` :c:data:`code_units` struct. For proper coordinates,
+   velocity units are equal to :c:data:`length_units` / :c:data:`time_units`.
+   For comoving coordinates, velocity units are equal to (:c:data:`length_units`
+   / :c:data:`a_value`) / :c:data:`time_units`. Note, this function only returns
+   a value, but does not set it in the struct. To set the value in the struct, use
+   :c:data:`set_velocity_units`.
+
+   :param code_units* my_units: code units conversions
+   :rtype: double
+   :returns: velocity_units
+
+.. c:function:: double get_temperature_units(code_units *my_units);
+
+   Returns the conversion factor between specific internal energy and temperature
+   assuming gamma (the adiabatic index) = 1, such that temperature in K is equal to
+   :c:data:`internal_energy` * ``temperature_units``. This unit conversion is
+   defined as m\ :sub:`H` * :c:data:`velocity_units`\ :sup:`2` / k\ :sub:`b`,
+   where m\ :sub:`H` is the Hydrogen mass and k\ :sub:`b` is the Boltzmann constant.
+
+   :param code_units* my_units: code units conversions
+   :rtype: double
+   :returns: temperature_units
+
 .. c:function:: int solve_chemistry(code_units *my_units, grackle_field_data *my_fields, double dt_value);
 
    Evolves the species densities and internal energies over a given timestep 

src/clib/Make.config.objects

@@ -28,6 +28,7 @@ OBJS_CONFIG_LIB = \
         cool1d_cloudy_old_tables_g.lo \
         cool1d_multi_g.lo \
         cool_multi_time_g.lo \
+        grackle_units.lo \
         initialize_chemistry_data.lo \
         initialize_cloudy_data.lo \
         initialize_UVbackground_data.lo \

src/clib/calculate_cooling_time.c

@@ -24,6 +24,8 @@ 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,
@@ -142,7 +144,7 @@ int local_calculate_cooling_time(chemistry_data *my_chemistry,
 
   /* Calculate temperature units. */
 
-  double temperature_units = mh * POW(my_units->velocity_units, 2) / kboltz;
+  double temperature_units = get_temperature_units(my_units);
 
   /* Call the fortran routine to solve cooling equations. */
 

src/clib/calculate_dust_temperature.c

@@ -27,6 +27,8 @@ 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,
@@ -74,7 +76,7 @@ int local_calculate_dust_temperature(chemistry_data *my_chemistry,
     co_density_units = my_units->density_units /
       POW(my_units->a_value * my_units->a_units, 3);
   }
-  double temperature_units = mh * POW(my_units->velocity_units, 2) / kboltz;
+  double temperature_units = get_temperature_units(my_units);
 
   /* Compute the size of the fields. */
 

src/clib/calculate_pressure.c

@@ -25,6 +25,8 @@
 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,
@@ -58,7 +60,7 @@ int local_calculate_pressure(chemistry_data *my_chemistry,
 
     /* Calculate temperature units. */
 
-    double temperature_units =  mh * POW(my_units->velocity_units, 2) / kboltz;
+    double temperature_units = get_temperature_units(my_units);
 
     double number_density, nH2, GammaH2Inverse,
       GammaInverse = 1.0/(my_chemistry->Gamma-1.0), x, Gamma1, temp;

src/clib/calculate_temperature.c

@@ -46,6 +46,8 @@ 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,
@@ -86,7 +88,7 @@ int local_calculate_temperature(chemistry_data *my_chemistry,
 
   /* Calculate temperature units. */
 
-  double temperature_units =  mh * POW(my_units->velocity_units, 2) / kboltz;
+  double temperature_units = get_temperature_units(my_units);
 
   double number_density, tiny_number = 1.-20;
   double inv_metal_mol = 1.0 / MU_METAL;
@@ -176,7 +178,7 @@ int local_calculate_temperature_table(chemistry_data *my_chemistry,
 
   /* Calculate temperature units. */
 
-  double temperature_units = mh * POW(my_units->velocity_units, 2) / kboltz;
+  double temperature_units = get_temperature_units(my_units);
 
   FORTRAN_NAME(calc_temp_cloudy_g)(
         my_fields->density,

src/clib/grackle.h

@@ -22,6 +22,12 @@ extern int grackle_verbose;
 extern chemistry_data *grackle_data;
 extern chemistry_data_storage grackle_rates;
 
+double get_velocity_units(code_units *my_units);
+
+void set_velocity_units(code_units *my_units);
+
+double get_temperature_units(code_units *my_units);
+
 int set_default_chemistry_parameters(chemistry_data *my_grackle);
 
 chemistry_data _set_default_chemistry_parameters(void);

src/clib/grackle_fortran_interface.def

@@ -141,6 +141,29 @@ cc       INTEGER(C_INT) :: omp_nthreads // not supported in fortran
 
 c     The following define the fortran interfaces to the C routines
 
+      INTERFACE
+         REAL(C_DOUBLE) FUNCTION get_temperature_units(my_units)
+     &        bind(C)
+            IMPORT
+            TYPE(grackle_units), INTENT(IN) :: my_units
+         END FUNCTION get_temperature_units
+      END INTERFACE
+
+      INTERFACE
+         REAL(C_DOUBLE) FUNCTION get_velocity_units(my_units)
+     &        bind(C)
+            IMPORT
+            TYPE(grackle_units), INTENT(IN) :: my_units
+         END FUNCTION get_velocity_units
+      END INTERFACE
+
+      INTERFACE
+         SUBROUTINE set_velocity_units(my_units) bind(C)
+            IMPORT
+            TYPE(grackle_units), INTENT(INOUT) :: my_units
+         END SUBROUTINE set_velocity_units
+      END INTERFACE
+
       INTERFACE
          INTEGER(C_INT) FUNCTION set_default_chemistry_parameters
      &        (my_grackle) bind(C)

src/clib/grackle_units.c

@@ -0,0 +1,42 @@
+/***********************************************************************
+/
+/ Initialize chemistry and cooling rate data
+/
+/
+/ Copyright (c) Enzo/Grackle Development Team. All rights reserved.
+/
+/ Distributed under the terms of the Enzo Public Licence.
+/
+/ The full license is in the file LICENSE, distributed with this 
+/ software.
+************************************************************************/
+
+#include <stdlib.h>
+#include <string.h>
+#include <stdio.h>
+#include <time.h>
+#include <math.h>
+#include "grackle_macros.h"
+#include "grackle_types.h"
+#include "grackle_chemistry_data.h"
+#include "phys_constants.h"
+
+double get_velocity_units(code_units *my_units)
+{
+  double velocity_units = my_units->length_units / my_units->time_units;
+  if (my_units->comoving_coordinates == 1) {
+    velocity_units /= my_units->a_value;
+  }
+  return velocity_units;
+}
+
+void set_velocity_units(code_units *my_units)
+{
+  my_units->velocity_units = get_velocity_units(my_units);
+}
+
+double get_temperature_units(code_units *my_units)
+{
+  double velocity_units = get_velocity_units(my_units);
+  return mh * POW(velocity_units, 2) / kboltz;
+}

src/clib/initialize_chemistry_data.c

@@ -11,7 +11,7 @@
 / software.
 ************************************************************************/
 
-#include <stdlib.h> 
+#include <stdlib.h>
 #include <string.h>
 #include <stdio.h>
 #include <time.h>

src/clib/solve_chemistry.c

@@ -23,6 +23,8 @@ 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,
@@ -163,7 +165,7 @@ int local_solve_chemistry(chemistry_data *my_chemistry,
 
   /* Calculate temperature units. */
 
-  double temperature_units =  mh * POW(my_units->velocity_units, 2) / kboltz;
+  double temperature_units = get_temperature_units(my_units);
 
   /* Call the fortran routine to solve cooling equations. */
 

src/clib/solve_rate_cool_g.F

@@ -333,9 +333,8 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de,
       xbase1   = uxyz/(aye*uaye)    ! uxyz is [x]*a      = [x]*[a]*a'        '
       dbase1   = urho*(aye*uaye)**3 ! urho is [dens]/a^3 = [dens]/([a]*a')^3 '
       coolunit = (uaye**5 * xbase1**2 * mh**2) / (tbase1**3 * dbase1)
-      uvel     = uxyz / utim
-c      chunit = (7.17775e-12_DKIND)/(2._DKIND*uvel*uvel*mh)   ! 4.5 eV per H2 formed
-      chunit = (1.60218e-12_DKIND)/(2._DKIND*uvel*uvel*mh)   ! 1 eV per H2 formed
+      uvel = (uxyz/aye) / utim
+      chunit   = (1.60218e-12_DKIND)/(2._DKIND*uvel*uvel*mh)   ! 1 eV per H2 formed
 
       dx_cgs = dx * xbase1
       c_ljeans = sqrt((gamma * pi * kboltz) /

src/example/c_example.c

@@ -43,10 +43,10 @@ int main(int argc, char *argv[])
   my_units.density_units = 1.67e-24;
   my_units.length_units = 1.0;
   my_units.time_units = 1.0e12;
-  my_units.velocity_units = my_units.length_units / my_units.time_units;
   my_units.a_units = 1.0; // units for the expansion factor
   // Set expansion factor to 1 for non-cosmological simulation.
   my_units.a_value = 1. / (1. + initial_redshift) / my_units.a_units;
+  set_velocity_units(&my_units);
 
   // Second, create a chemistry object for parameters.  This needs to be a pointer.
   chemistry_data *my_grackle_data;
@@ -133,9 +133,7 @@ int main(int argc, char *argv[])
   my_fields.RT_heating_rate = malloc(field_size * sizeof(gr_float));
 
   // set temperature units
-  double temperature_units =  mh * pow(my_units.a_units * 
-                                         my_units.length_units /
-                                         my_units.time_units, 2) / kboltz;
+  double temperature_units = get_temperature_units(&my_units);
 
   for (i = 0;i < field_size;i++) {
     my_fields.density[i] = 1.0;

src/example/cxx_example.C

@@ -45,10 +45,10 @@ int main(int argc, char *argv[])
   my_units.density_units = 1.67e-24;
   my_units.length_units = 1.0;
   my_units.time_units = 1.0e12;
-  my_units.velocity_units = my_units.length_units / my_units.time_units;
   my_units.a_units = 1.0; // units for the expansion factor
   // Set expansion factor to 1 for non-cosmological simulation.
   my_units.a_value = 1. / (1. + initial_redshift) / my_units.a_units;
+  set_velocity_units(&my_units);
 
   // Second, create a chemistry object for parameters.  This needs to be a pointer.
   chemistry_data *my_grackle_data;
@@ -136,9 +136,7 @@ int main(int argc, char *argv[])
   my_fields.isrf_habing = new gr_float[field_size];
 
   // set temperature units
-  double temperature_units = mh * pow(my_units.a_units * 
-                                      my_units.length_units /
-                                      my_units.time_units, 2) / kboltz;
+  double temperature_units = get_temperature_units(&my_units);
 
   int i;
   for (i = 0;i < field_size;i++) {