Merge branch 'utilities' into 'master'

Add generic utility function to find a value in an array

See merge request !547

.gitignore

@@ -116,6 +116,7 @@ tests/testPotentialPair
 tests/testEOS
 tests/testEOS*.txt
 tests/testEOS*.png
+tests/testUtilities
 
 theory/latex/swift.pdf
 theory/SPH/Kernels/kernels.pdf

examples/MoonFormingImpact/moon_forming_impact.yml

@@ -42,3 +42,7 @@ Gravity:
 InitialConditions:
                                         # The initial conditions file to read
     file_name:  moon_forming_impact.hdf5
+
+# Parameters related to the equation of state
+EoS:
+    planetary_use_Til:    1                       # Whether to prepare the Tillotson EOS

examples/UranusImpact/uranus_impact.yml

@@ -41,3 +41,11 @@ Gravity:
 # Parameters related to the initial conditions
 InitialConditions:
     file_name:      uranus_impact.hdf5  # The initial conditions file to read
+
+# Parameters related to the equation of state
+EoS:
+    planetary_use_HM80:   1                       # Whether to prepare the Hubbard & MacFarlane (1980) EOS
+                                        # Table file paths
+    planetary_HM80_HHe_table_file:    /gpfs/data/dc-kege1/gihr_data/P_rho_u_HHe.txt
+    planetary_HM80_ice_table_file:    /gpfs/data/dc-kege1/gihr_data/P_rho_u_ice.txt
+    planetary_HM80_rock_table_file:   /gpfs/data/dc-kege1/gihr_data/P_rho_u_roc.txt

examples/parameter_example.yml

@@ -131,7 +131,16 @@ DomainDecomposition:
 # Parameters related to the equation of state ------------------------------------------
 
 EoS:
-  isothermal_internal_energy: 20.26784 # Thermal energy per unit mass for the case of isothermal equation of state (in internal units).
+  isothermal_internal_energy: 20.26784  # Thermal energy per unit mass for the case of isothermal equation of state (in internal units).
+
+  planetary_use_Til:    1   # (Optional) Whether to prepare the Tillotson EOS
+  planetary_use_HM80:   0   # (Optional) Whether to prepare the Hubbard & MacFarlane (1980) EOS
+  planetary_use_ANEOS:  0   # (Optional) Whether to prepare the ANEOS EOS
+  planetary_use_SESAME: 0   # (Optional) Whether to prepare the SESAME EOS
+                            # (Optional) Table file paths
+  planetary_HM80_HHe_table_file:    HM80_HHe.txt
+  planetary_HM80_ice_table_file:    HM80_ice.txt
+  planetary_HM80_rock_table_file:   HM80_rock.txt
 
 # Parameters related to external potentials --------------------------------------------
 

src/equation_of_state/planetary/equation_of_state.h

@@ -25,7 +25,7 @@
  *
  * For any/all of the planetary EOS. Each EOS type's functions are set in its
  * own header file: `equation_of_state/planetary/<eos_type>.h`.
- * See `material_id` for the available choices.
+ * See `eos_planetary_material_id` for the available choices.
  *
  * Not all functions are implemented for all EOS types, so not all can be used
  * with all hydro formulations yet.
@@ -95,7 +95,7 @@ enum eos_planetary_material_id {
   eos_planetary_id_ANEOS_iron =
       eos_planetary_type_ANEOS * eos_planetary_type_factor,
 
-  /*! ANEOS forsterite */
+  /*! MANEOS forsterite */
   eos_planetary_id_MANEOS_forsterite =
       eos_planetary_type_ANEOS * eos_planetary_type_factor + 1,
 
@@ -1077,46 +1077,62 @@ __attribute__((always_inline)) INLINE static void eos_init(
     struct eos_parameters *e, const struct phys_const *phys_const,
     const struct unit_system *us, const struct swift_params *params) {
 
+  // Table file names
+  char HM80_HHe_table_file[PARSER_MAX_LINE_SIZE];
+  char HM80_ice_table_file[PARSER_MAX_LINE_SIZE];
+  char HM80_rock_table_file[PARSER_MAX_LINE_SIZE];
+
   // Set the parameters and material IDs, load tables, etc. for each material
+  // and convert to internal units
   // Tillotson
-  set_Til_iron(&e->Til_iron, eos_planetary_id_Til_iron);
-  set_Til_granite(&e->Til_granite, eos_planetary_id_Til_granite);
-  set_Til_water(&e->Til_water, eos_planetary_id_Til_water);
+  if (parser_get_opt_param_int(params, "EoS:planetary_use_Til", 0)) {
+      set_Til_iron(&e->Til_iron, eos_planetary_id_Til_iron);
+      set_Til_granite(&e->Til_granite, eos_planetary_id_Til_granite);
+      set_Til_water(&e->Til_water, eos_planetary_id_Til_water);
+
+      convert_units_Til(&e->Til_iron, us);
+      convert_units_Til(&e->Til_granite, us);
+      convert_units_Til(&e->Til_water, us);
+  }
 
   // Hubbard & MacFarlane (1980)
-  set_HM80_HHe(&e->HM80_HHe, eos_planetary_id_HM80_HHe);
-  set_HM80_ice(&e->HM80_ice, eos_planetary_id_HM80_ice);
-  set_HM80_rock(&e->HM80_rock, eos_planetary_id_HM80_rock);
-
-  load_HM80_table(&e->HM80_HHe, HM80_HHe_table_file);
-  load_HM80_table(&e->HM80_ice, HM80_ice_table_file);
-  load_HM80_table(&e->HM80_rock, HM80_rock_table_file);
+  if (parser_get_opt_param_int(params, "EoS:planetary_use_HM80", 0)) {
+      set_HM80_HHe(&e->HM80_HHe, eos_planetary_id_HM80_HHe);
+      set_HM80_ice(&e->HM80_ice, eos_planetary_id_HM80_ice);
+      set_HM80_rock(&e->HM80_rock, eos_planetary_id_HM80_rock);
+
+      parser_get_param_string(params, "EoS:planetary_HM80_HHe_table_file",
+                              HM80_HHe_table_file);
+      parser_get_param_string(params, "EoS:planetary_HM80_ice_table_file",
+                              HM80_ice_table_file);
+      parser_get_param_string(params, "EoS:planetary_HM80_rock_table_file",
+                              HM80_rock_table_file);
+
+      load_HM80_table(&e->HM80_HHe, HM80_HHe_table_file);
+      load_HM80_table(&e->HM80_ice, HM80_ice_table_file);
+      load_HM80_table(&e->HM80_rock, HM80_rock_table_file);
+
+      convert_units_HM80(&e->HM80_HHe, us);
+      convert_units_HM80(&e->HM80_ice, us);
+      convert_units_HM80(&e->HM80_rock, us);
+  }
 
   // ANEOS
-  set_ANEOS_iron(&e->ANEOS_iron, eos_planetary_id_ANEOS_iron);
-  set_MANEOS_forsterite(&e->MANEOS_forsterite,
-                        eos_planetary_id_MANEOS_forsterite);
-
-  // SESAME
-  set_SESAME_iron(&e->SESAME_iron, eos_planetary_id_SESAME_iron);
-
-  // Convert to internal units
-  // Tillotson
-  convert_units_Til(&e->Til_iron, us);
-  convert_units_Til(&e->Til_granite, us);
-  convert_units_Til(&e->Til_water, us);
-
-  // Hubbard & MacFarlane (1980)
-  convert_units_HM80(&e->HM80_HHe, us);
-  convert_units_HM80(&e->HM80_ice, us);
-  convert_units_HM80(&e->HM80_rock, us);
+  if (parser_get_opt_param_int(params, "EoS:planetary_use_ANEOS", 0)) {
+      set_ANEOS_iron(&e->ANEOS_iron, eos_planetary_id_ANEOS_iron);
+      set_MANEOS_forsterite(&e->MANEOS_forsterite,
+                            eos_planetary_id_MANEOS_forsterite);
 
-  // ANEOS
-  convert_units_ANEOS(&e->ANEOS_iron, us);
-  convert_units_ANEOS(&e->MANEOS_forsterite, us);
+      convert_units_ANEOS(&e->ANEOS_iron, us);
+      convert_units_ANEOS(&e->MANEOS_forsterite, us);
+  }
 
   // SESAME
-  convert_units_SESAME(&e->SESAME_iron, us);
+  if (parser_get_opt_param_int(params, "EoS:planetary_use_SESAME", 0)) {
+    set_SESAME_iron(&e->SESAME_iron, eos_planetary_id_SESAME_iron);
+
+    convert_units_SESAME(&e->SESAME_iron, us);
+  }
 }
 
 /**

src/equation_of_state/planetary/hm80.h

@@ -41,19 +41,13 @@
 
 // Hubbard & MacFarlane (1980) parameters
 struct HM80_params {
-  float **table_P_rho_u;
+  float *table_P_rho_u;
   int num_rho, num_u;
   float log_rho_min, log_rho_max, log_rho_step, inv_log_rho_step, log_u_min,
       log_u_max, log_u_step, inv_log_u_step, bulk_mod;
   enum eos_planetary_material_id mat_id;
 };
 
-// Table file names
-/// to be read in from the parameter file instead once finished testing...
-#define HM80_HHe_table_file "/gpfs/data/dc-kege1/gihr_data/P_rho_u_HHe.txt"
-#define HM80_ice_table_file "/gpfs/data/dc-kege1/gihr_data/P_rho_u_ice.txt"
-#define HM80_rock_table_file "/gpfs/data/dc-kege1/gihr_data/P_rho_u_roc.txt"
-
 // Parameter values for each material (cgs units)
 INLINE static void set_HM80_HHe(struct HM80_params *mat,
                                 enum eos_planetary_material_id mat_id) {
@@ -107,17 +101,15 @@ INLINE static void set_HM80_rock(struct HM80_params *mat,
 // Read the table from file
 INLINE static void load_HM80_table(struct HM80_params *mat, char *table_file) {
   // Allocate table memory
-  mat->table_P_rho_u = (float **)malloc(mat->num_rho * sizeof(float *));
-  for (int i = 0; i < mat->num_rho; i++) {
-    mat->table_P_rho_u[i] = (float *)malloc(mat->num_u * sizeof(float));
-  }
+  mat->table_P_rho_u = (float *)malloc(mat->num_rho * mat->num_u *
+                                       sizeof(float *));
 
   // Load table contents from file
   FILE *f = fopen(table_file, "r");
   int c;
   for (int i = 0; i < mat->num_rho; i++) {
     for (int j = 0; j < mat->num_u; j++) {
-      c = fscanf(f, "%f", &mat->table_P_rho_u[i][j]);
+      c = fscanf(f, "%f", &mat->table_P_rho_u[i*mat->num_rho + j]);
       if (c != 1) {
         error("Failed to read EOS table");
       }
@@ -147,7 +139,7 @@ INLINE static void convert_units_HM80(struct HM80_params *mat,
   // Table Pressures in Mbar
   for (int i = 0; i < mat->num_rho; i++) {
     for (int j = 0; j < mat->num_u; j++) {
-      mat->table_P_rho_u[i][j] *=
+      mat->table_P_rho_u[i*mat->num_rho + j] *=
           Mbar_to_Ba / units_cgs_conversion_factor(us, UNIT_CONV_PRESSURE);
     }
   }
@@ -194,7 +186,6 @@ INLINE static float HM80_soundspeed_from_entropy(
 // gas_entropy_from_internal_energy
 INLINE static float HM80_entropy_from_internal_energy(
     float density, float u, const struct HM80_params *mat) {
-  error("This EOS function is not yet implemented!");
 
   return 0;
 }
@@ -205,7 +196,7 @@ INLINE static float HM80_pressure_from_internal_energy(
 
   float P;
 
-  if (u <= 0) {
+  if (u <= 0.f) {
     return 0.f;
   }
 
@@ -226,32 +217,32 @@ INLINE static float HM80_pressure_from_internal_energy(
   // Return zero pressure if below the table minimum/a
   // Extrapolate the pressure for low densities
   if (rho_idx < 0) {  // Too-low rho
-    P = expf(logf((1 - intp_u) * mat->table_P_rho_u[0][u_idx] +
-                  intp_u * mat->table_P_rho_u[0][u_idx + 1]) +
+    P = expf(logf((1 - intp_u) * mat->table_P_rho_u[u_idx] +
+                  intp_u * mat->table_P_rho_u[u_idx + 1]) +
              log_rho - mat->log_rho_min);
     if (u_idx < 0) {  // and too-low u
-      P = 0;
+      P = 0.f;
     }
   } else if (u_idx < 0) {  // Too-low u
-    P = 0;
+    P = 0.f;
   }
   // Return an edge value if above the table maximum/a
   else if (rho_idx >= mat->num_rho - 1) {  // Too-high rho
     if (u_idx >= mat->num_u - 1) {         // and too-high u
-      P = mat->table_P_rho_u[mat->num_rho - 1][mat->num_u - 1];
+      P = mat->table_P_rho_u[(mat->num_rho - 1)*mat->num_u + mat->num_u - 1];
     } else {
-      P = mat->table_P_rho_u[mat->num_rho - 1][u_idx];
+      P = mat->table_P_rho_u[(mat->num_rho - 1)*mat->num_u + u_idx];
     }
   } else if (u_idx >= mat->num_u - 1) {  // Too-high u
-    P = mat->table_P_rho_u[rho_idx][mat->num_u - 1];
+    P = mat->table_P_rho_u[rho_idx*mat->num_u + mat->num_u - 1];
   }
   // Normal interpolation within the table
   else {
     P = (1.f - intp_rho) *
-            ((1.f - intp_u) * mat->table_P_rho_u[rho_idx][u_idx] +
-             intp_u * mat->table_P_rho_u[rho_idx][u_idx + 1]) +
-        intp_rho * ((1 - intp_u) * mat->table_P_rho_u[rho_idx + 1][u_idx] +
-                    intp_u * mat->table_P_rho_u[rho_idx + 1][u_idx + 1]);
+            ((1.f - intp_u) * mat->table_P_rho_u[rho_idx*mat->num_u + u_idx] +
+             intp_u * mat->table_P_rho_u[rho_idx*mat->num_u + u_idx + 1]) +
+        intp_rho * ((1 - intp_u) * mat->table_P_rho_u[(rho_idx + 1)*mat->num_u + u_idx] +
+                    intp_u * mat->table_P_rho_u[(rho_idx + 1)*mat->num_u + u_idx + 1]);
   }
 
   return P;

src/equation_of_state/planetary/tillotson.h

@@ -147,7 +147,6 @@ INLINE static float Til_soundspeed_from_entropy(float density, float entropy,
 // gas_entropy_from_internal_energy
 INLINE static float Til_entropy_from_internal_energy(
     float density, float u, const struct Til_params *mat) {
-  error("This EOS function is not yet implemented!");
 
   return 0;
 }

src/hydro/MinimalMultiMat/hydro_io.h

@@ -72,6 +72,12 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
       io_make_input_field("MaterialID", INT, 1, OPTIONAL, 1, parts, mat_id);
 }
 
+void convert_S(const struct engine* e, const struct part* p,
+               const struct xpart* xp, float* ret) {
+
+  ret[0] = hydro_get_comoving_entropy(p);
+}
+
 void convert_P(const struct engine* e, const struct part* p,
                const struct xpart* xp, float* ret) {
 
@@ -146,7 +152,7 @@ void convert_part_potential(const struct engine* e, const struct part* p,
 void hydro_write_particles(const struct part* parts, const struct xpart* xparts,
                            struct io_props* list, int* num_fields) {
 
-  *num_fields = 10;
+  *num_fields = 11;
 
   /* List what we want to write */
   list[0] = io_make_output_field_convert_part("Coordinates", DOUBLE, 3,
@@ -164,13 +170,16 @@ void hydro_write_particles(const struct part* parts, const struct xpart* xparts,
                                  UNIT_CONV_NO_UNITS, parts, id);
   list[6] =
       io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
-  list[7] = io_make_output_field("MaterialID", INT, 1, UNIT_CONV_NO_UNITS,
+  list[7] = io_make_output_field_convert_part("Entropy", FLOAT, 1,
+                                          UNIT_CONV_ENTROPY_PER_UNIT_MASS,
+                                          parts, xparts, convert_S);
+  list[8] = io_make_output_field("MaterialID", INT, 1, UNIT_CONV_NO_UNITS,
                                  parts, mat_id);
-  list[8] = io_make_output_field_convert_part(
+  list[9] = io_make_output_field_convert_part(
       "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, xparts, convert_P);
-  list[9] = io_make_output_field_convert_part("Potential", FLOAT, 1,
-                                              UNIT_CONV_POTENTIAL, parts,
-                                              xparts, convert_part_potential);
+  list[10] = io_make_output_field_convert_part("Potential", FLOAT, 1,
+                                               UNIT_CONV_POTENTIAL, parts,
+                                               xparts, convert_part_potential);
 }
 
 /**

src/utilities.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2018 Jacob Kegerreis (jacob.kegerreis@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#ifndef SWIFT_UTILITIES_H
+#define SWIFT_UTILITIES_H
+
+/**
+ * @brief Search for a value in a monotonically increasing array to find the
+ *      index such that array[index] < value < array[index + 1]
+ *
+ * @param x The value to find
+ * @param array The array to search
+ * @param n The length of the array
+ *
+ * Return -1 and n for x below and above the array edge values respectively.
+ */
+INLINE static int find_value_in_monot_incr_array(
+    const float x, const float *array, const int n) {
+
+    int index_mid, index_low = 0, index_high = n;
+
+    // Until array[index_low] < x < array[index_high=index_low+1]
+    while (index_high - index_low > 1) {
+        index_mid = (index_high + index_low) / 2;  // Middle index
+
+        // Replace the low or high index with the middle
+        if (array[index_mid] <= x)
+            index_low = index_mid;
+        else
+            index_high = index_mid;
+    }
+
+    // Set index with the found index_low or an error value if outside the array
+    if (x < array[0])
+        return -1;
+    else if (array[n-1] <= x)
+        return n;
+    else
+        return index_low;
+}
+
+#endif /* SWIFT_UTILITIES_H */

tests/Makefile.am

@@ -27,7 +27,7 @@ TESTS = testGreetings testMaths testReading.sh testSingle testKernel testSymmetr
         testMatrixInversion testThreadpool testDump testLogger testInteractions.sh \
         testVoronoi1D testVoronoi2D testVoronoi3D testGravityDerivatives \
 	testPeriodicBC.sh testPeriodicBCPerturbed.sh testPotentialSelf \
-	testPotentialPair testEOS
+	testPotentialPair testEOS testUtilities
 
 # List of test programs to compile
 check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
@@ -37,7 +37,7 @@ check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
                  testAdiabaticIndex testRiemannExact testRiemannTRRS \
                  testRiemannHLLC testMatrixInversion testDump testLogger \
 		 testVoronoi1D testVoronoi2D testVoronoi3D testPeriodicBC \
-		 testGravityDerivatives testPotentialSelf testPotentialPair testEOS
+		 testGravityDerivatives testPotentialSelf testPotentialPair testEOS testUtilities
 
 # Rebuild tests when SWIFT is updated.
 $(check_PROGRAMS): ../src/.libs/libswiftsim.a
@@ -107,6 +107,8 @@ testPotentialPair_SOURCES = testPotentialPair.c
 
 testEOS_SOURCES = testEOS.c
 
+testUtilities_SOURCES = testUtilities.c
+
 # Files necessary for distribution
 EXTRA_DIST = testReading.sh makeInput.py testActivePair.sh \
 	     test27cells.sh test27cellsPerturbed.sh testParser.sh testPeriodicBC.sh \

tests/testUtilities.c

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (C) 2018 Jacob Kegerreis (jacob.kegerreis@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#include "swift.h"
+#include "utilities.h"
+
+/**
+ * @brief Test generic utility functions
+ */
+int main() {
+    /// Test find_value_in_monot_incr_array()
+    int n = 100;
+    float array[n];
+    int index;
+    float x;
+
+    // Initialise test array
+    for (int j = 0; j < n; j++) {
+        array[j] = j;
+    }
+
+    // Typical value
+    x = 42.42f;
+    index = find_value_in_monot_incr_array(x, array, n);
+    if (index != 42) {
+        error("Failed with a typical value ");
+    }
+
+    // Value on array element
+    x = 33.f;
+    index = find_value_in_monot_incr_array(x, array, n);
+    if (index != 33) {
+        error("Failed with an array element ");
+    }
+
+    // Value below array
+    x = -123.f;
+    index = find_value_in_monot_incr_array(x, array, n);
+    if (index != -1) {
+        error("Failed with a value below the array ");
+    }
+
+    // Value above array
+    x = 123.f;
+    index = find_value_in_monot_incr_array(x, array, n);
+    if (index != n) {
+        error("Failed with a value above the array ");
+    }
+
+    // Array slice with typical value
+    x = 9.81f;
+    n = 10;
+    index = find_value_in_monot_incr_array(x, array + 5, n);
+    if (index != 4) {
+        error("Failed with an array slice ");
+    }
+
+    return 0;
+}