diff --git a/.gitignore b/.gitignore index 4a8a0b583aec90cffe3e450155fa54a4885778f0..b6f80a5ad0a549843b6261c439bcc2a704b43124 100644 --- a/.gitignore +++ b/.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 diff --git a/examples/MoonFormingImpact/moon_forming_impact.yml b/examples/MoonFormingImpact/moon_forming_impact.yml index 7d726bc02ba4fb6c8dd02f74907ffc48c3ed9431..323adf7f3ac73f41b45b50eaa76a95033dca35d7 100644 --- a/examples/MoonFormingImpact/moon_forming_impact.yml +++ b/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 diff --git a/examples/UranusImpact/uranus_impact.yml b/examples/UranusImpact/uranus_impact.yml index aae9f66847a9f9b55984ea5d2ea1c79099c01e95..fabddca00f80fcdd79ff6114ff0544cd251046f4 100644 --- a/examples/UranusImpact/uranus_impact.yml +++ b/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 diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml index c9bc9811803051b274cd09cdcf6b58d1a65a0d7d..791db2758290e400d4ed9ffe5b8e0d4303057874 100644 --- a/examples/parameter_example.yml +++ b/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 -------------------------------------------- diff --git a/src/equation_of_state/planetary/equation_of_state.h b/src/equation_of_state/planetary/equation_of_state.h index fcc1c959152801253d465915d8b28b1552c913b8..b82b26f328b1f0757b0edd8e1e2e5188d961dd4c 100644 --- a/src/equation_of_state/planetary/equation_of_state.h +++ b/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); + } } /** diff --git a/src/equation_of_state/planetary/hm80.h b/src/equation_of_state/planetary/hm80.h index b0669c129b2b4d642ae27e58f2f2bec221e3df9b..69f213f38184bcc15ec3ad82a66b7f197a7bb502 100644 --- a/src/equation_of_state/planetary/hm80.h +++ b/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; diff --git a/src/equation_of_state/planetary/tillotson.h b/src/equation_of_state/planetary/tillotson.h index ed280079b040776dfc20480c50fbc9a3b6ff4470..d5b6d5c35d5edf9e114fe7f010c4f5b1e2327a83 100644 --- a/src/equation_of_state/planetary/tillotson.h +++ b/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; } diff --git a/src/hydro/MinimalMultiMat/hydro_io.h b/src/hydro/MinimalMultiMat/hydro_io.h index 6e55497dc42111de291b8984a7dca048bb774722..5419582f3c5e268eb0dbf59bb5544bf700973838 100644 --- a/src/hydro/MinimalMultiMat/hydro_io.h +++ b/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); } /** diff --git a/src/utilities.h b/src/utilities.h new file mode 100644 index 0000000000000000000000000000000000000000..12038a6fc82df6861430ec56b2e8a702c1d00b27 --- /dev/null +++ b/src/utilities.h @@ -0,0 +1,58 @@ +/******************************************************************************* + * 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 */ diff --git a/tests/Makefile.am b/tests/Makefile.am index fb3384fa56936b3c7c1fe1b415a2201e419ba566..77aad7beb58c8ab80da04c28fc68b64afba3d99a 100644 --- a/tests/Makefile.am +++ b/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 \ diff --git a/tests/testUtilities.c b/tests/testUtilities.c new file mode 100644 index 0000000000000000000000000000000000000000..9234141e98eac63d39bc70bd79ab3204c7ee106d --- /dev/null +++ b/tests/testUtilities.c @@ -0,0 +1,75 @@ +/******************************************************************************* + * 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; +}