When compiling with the iosthermal equation of state, read the constant...

When compiling with the iosthermal equation of state, read the constant thermal energy from the YAML file.

examples/DiscPatch/HydroStatic/disc-patch-icc.yml

@@ -31,6 +31,9 @@ SPH:
   max_ghost_iterations:  30       # Maximal number of iterations allowed to converge towards the smoothing length.
   h_max:                 60.      # Maximal smoothing length allowed (in internal units).
 
+EoS:
+  isothermal_internal_energy: 20.267845
+  
 # Parameters related to the initial conditions
 InitialConditions:
   file_name:  Disc-Patch.hdf5       # The file to read

examples/DiscPatch/HydroStatic/makeIC.py

@@ -202,7 +202,5 @@ print "--- Runtime parameters (YAML file): ---"
 print "DiscPatchPotential:surface_density:    ", surface_density
 print "DiscPatchPotential:scale_height:       ", scale_height
 print "DiscPatchPotential:x_disc:             ", 0.5 * boxSize_x
+print "EoS:isothermal_internal_energy: %ef"%u_therm
 print ""
-
-print "--- Constant parameters: ---"
-print "const_isothermal_internal_energy: %ef"%u_therm

examples/main.c

@@ -463,6 +463,7 @@ int main(int argc, char *argv[]) {
   /* Initialise the hydro properties */
   struct hydro_props hydro_properties;
   if (with_hydro) hydro_props_init(&hydro_properties, params);
+  if (with_hydro) eos_init(&eos, params);
 
   /* Initialise the gravity properties */
   struct gravity_props gravity_properties;

examples/parameter_example.yml

@@ -50,6 +50,9 @@ SPH:
   max_volume_change:     1.4      # (Optional) Maximal allowed change of kernel volume over one time-step.
   max_ghost_iterations:  30       # (Optional) Maximal number of iterations allowed to converge towards the smoothing length.
 
+EoS:
+  isothermal_internal_energy: 20.26784 # Thermal energy per unit mass for the case of isothermal equation of state.
+  
 # Parameters for the self-gravity scheme
 Gravity:
   eta:          0.025    # Constant dimensionless multiplier for time integration.

src/Makefile.am

@@ -57,7 +57,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
     runner_doiact_fft.c threadpool.c cooling.c sourceterms.c \
     statistics.c runner_doiact_vec.c profiler.c dump.c logger.c \
     part_type.c xmf.c gravity_properties.c gravity.c \
-    collectgroup.c hydro_space.c
+    collectgroup.c hydro_space.c equation_of_state.c
 
 # Include files for distribution, not installation.
 nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \

src/const.h

@@ -36,9 +36,6 @@
 /* Time integration constants. */
 #define const_max_u_change 0.1f
 
-/* Thermal energy per unit mass used as a constant for the isothermal EoS */
-#define const_isothermal_internal_energy 20.2678457288f
-
 /* Type of gradients to use (GIZMO_SPH only) */
 /* If no option is chosen, no gradients are used (first order scheme) */
 //#define GRADIENTS_SPH

src/equation_of_state.c

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2017 Matthieu Schaller (matthieu.schaller@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/>.
+ *
+ ******************************************************************************/
+
+/* This object's header. */
+#include "equation_of_state.h"
+
+/* local headers */
+#include "common_io.h"
+
+/* Equation of state for the physics model
+ * (temporary ugly solution as a global variable) */
+struct eos_parameters eos;
+
+void eos_init(struct eos_parameters *e, const struct swift_params *params) {
+
+#if defined(EOS_IDEAL_GAS)
+/* nothing to do here */
+#elif defined(EOS_ISOTHERMAL_GAS)
+  e->isothermal_internal_energy =
+      parser_get_param_float(params, "EoS:isothermal_internal_energy");
+#endif
+}
+
+void eos_print(const struct eos_parameters *e) {
+
+#if defined(EOS_IDEAL_GAS)
+  message("Equation of state: Ideal gas.");
+#elif defined(EOS_ISOTHERMAL_GAS)
+  message(
+      "Equation of state: Isothermal with internal energy "
+      "per unit mass set to %f.",
+      e->isothermal_internal_energy);
+#endif
+
+  message("Adiabatic index gamma: %f.", hydro_gamma);
+}
+
+#if defined(HAVE_HDF5)
+void eos_print_snapshot(hid_t h_grpsph, const struct eos_parameters *e) {
+
+  io_write_attribute_f(h_grpsph, "Adiabatic index", hydro_gamma);
+
+#if defined(EOS_IDEAL_GAS)
+  io_write_attribute_s(h_grpsph, "Equation of state", "Ideal gas");
+#elif defined(EOS_ISOTHERMAL_GAS)
+  io_write_attribute_s(h_grpsph, "Equation of state", "Isothermal gas");
+  io_write_attribute_f(h_grpsph, "Thermal energy per unit mass",
+                       e->isothermal_internal_energy);
+#endif
+}
+#endif

src/equation_of_state.h

@@ -37,9 +37,13 @@
 #include "debug.h"
 #include "inline.h"
 
+extern struct eos_parameters eos;
+
 /* ------------------------------------------------------------------------- */
 #if defined(EOS_IDEAL_GAS)
 
+struct eos_parameters {};
+
 /**
  * @brief Returns the internal energy given density and entropy
  *
@@ -172,6 +176,12 @@ __attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure(
 /* ------------------------------------------------------------------------- */
 #elif defined(EOS_ISOTHERMAL_GAS)
 
+struct eos_parameters {
+
+  /*! Thermal energy per unit mass */
+  float isothermal_internal_energy;
+};
+
 /**
  * @brief Returns the internal energy given density and entropy
  *
@@ -184,7 +194,7 @@ __attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure(
 __attribute__((always_inline)) INLINE static float
 gas_internal_energy_from_entropy(float density, float entropy) {
 
-  return const_isothermal_internal_energy;
+  return eos.isothermal_internal_energy;
 }
 /**
  * @brief Returns the pressure given density and entropy
@@ -198,7 +208,7 @@ gas_internal_energy_from_entropy(float density, float entropy) {
 __attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
     float density, float entropy) {
 
-  return hydro_gamma_minus_one * const_isothermal_internal_energy * density;
+  return hydro_gamma_minus_one * eos.isothermal_internal_energy * density;
 }
 
 /**
@@ -213,7 +223,7 @@ __attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
 __attribute__((always_inline)) INLINE static float gas_entropy_from_pressure(
     float density, float pressure) {
 
-  return hydro_gamma_minus_one * const_isothermal_internal_energy *
+  return hydro_gamma_minus_one * eos.isothermal_internal_energy *
          pow_minus_gamma_minus_one(density);
 }
 
@@ -229,7 +239,7 @@ __attribute__((always_inline)) INLINE static float gas_entropy_from_pressure(
 __attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy(
     float density, float entropy) {
 
-  return sqrtf(const_isothermal_internal_energy * hydro_gamma *
+  return sqrtf(eos.isothermal_internal_energy * hydro_gamma *
                hydro_gamma_minus_one);
 }
 
@@ -245,7 +255,7 @@ __attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy(
 __attribute__((always_inline)) INLINE static float
 gas_entropy_from_internal_energy(float density, float u) {
 
-  return hydro_gamma_minus_one * const_isothermal_internal_energy *
+  return hydro_gamma_minus_one * eos.isothermal_internal_energy *
          pow_minus_gamma_minus_one(density);
 }
 
@@ -261,7 +271,7 @@ gas_entropy_from_internal_energy(float density, float u) {
 __attribute__((always_inline)) INLINE static float
 gas_pressure_from_internal_energy(float density, float u) {
 
-  return hydro_gamma_minus_one * const_isothermal_internal_energy * density;
+  return hydro_gamma_minus_one * eos.isothermal_internal_energy * density;
 }
 
 /**
@@ -275,7 +285,7 @@ gas_pressure_from_internal_energy(float density, float u) {
  */
 __attribute__((always_inline)) INLINE static float
 gas_internal_energy_from_pressure(float density, float pressure) {
-  return const_isothermal_internal_energy;
+  return eos.isothermal_internal_energy;
 }
 
 /**
@@ -290,7 +300,7 @@ gas_internal_energy_from_pressure(float density, float pressure) {
 __attribute__((always_inline)) INLINE static float
 gas_soundspeed_from_internal_energy(float density, float u) {
 
-  return sqrtf(const_isothermal_internal_energy * hydro_gamma *
+  return sqrtf(eos.isothermal_internal_energy * hydro_gamma *
                hydro_gamma_minus_one);
 }
 
@@ -306,7 +316,7 @@ gas_soundspeed_from_internal_energy(float density, float u) {
 __attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure(
     float density, float P) {
 
-  return sqrtf(const_isothermal_internal_energy * hydro_gamma *
+  return sqrtf(eos.isothermal_internal_energy * hydro_gamma *
                hydro_gamma_minus_one);
 }
 
@@ -317,4 +327,11 @@ __attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure(
 
 #endif
 
+void eos_init(struct eos_parameters *e, const struct swift_params *params);
+void eos_print(const struct eos_parameters *e);
+
+#if defined(HAVE_HDF5)
+void eos_print_snapshot(hid_t h_grpsph, const struct eos_parameters *e);
+#endif
+
 #endif /* SWIFT_EQUATION_OF_STATE_H */

src/hydro_properties.c

@@ -77,17 +77,10 @@ void hydro_props_init(struct hydro_props *p,
 
 void hydro_props_print(const struct hydro_props *p) {
 
-#if defined(EOS_IDEAL_GAS)
-  message("Equation of state: Ideal gas.");
-#elif defined(EOS_ISOTHERMAL_GAS)
-  message(
-      "Equation of state: Isothermal with internal energy "
-      "per unit mass set to %f.",
-      const_isothermal_internal_energy);
-#endif
-
-  message("Adiabatic index gamma: %f.", hydro_gamma);
+  /* Print equation of state first */
+  eos_print(&eos);
 
+  /* Now SPH */
   message("Hydrodynamic scheme: %s in %dD.", SPH_IMPLEMENTATION,
           (int)hydro_dimension);
 
@@ -115,7 +108,8 @@ void hydro_props_print(const struct hydro_props *p) {
 #if defined(HAVE_HDF5)
 void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
 
-  io_write_attribute_f(h_grpsph, "Adiabatic index", hydro_gamma);
+  eos_print_snapshot(h_grpsph, &eos);
+
   io_write_attribute_i(h_grpsph, "Dimension", (int)hydro_dimension);
   io_write_attribute_s(h_grpsph, "Scheme", SPH_IMPLEMENTATION);
   io_write_attribute_s(h_grpsph, "Kernel function", kernel_name);