testEOS.c 10.6 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76
/*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
 *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
 *               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/>.
 *
 ******************************************************************************/

/* Config parameters. */
#include "../config.h"

/* Some standard headers. */
#include <fenv.h>
#include <float.h>
#include <limits.h>
#include <math.h>
#include <pthread.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>

/* Conditional headers. */
#ifdef HAVE_LIBZ
#include <zlib.h>
#endif

/* Local headers. */
#include "equation_of_state.h"
#include "swift.h"

/* Engine policy flags. */
#ifndef ENGINE_POLICY
#define ENGINE_POLICY engine_policy_none
#endif

/**
 * @brief Write a list of densities, energies, and resulting pressures to file
 *  from an equation of state.
 *
 *                      WORK IN PROGRESS
 *
 * So far only does the Hubbard & MacFarlane (1980) equations of state.
 *
 * Usage:
 *      $  ./testEOS  (mat_id)  (do_output)
 *
 * Sys args (optional):
 *      mat_id | int | Material ID, see equation_of_state.h for the options.
 *          Default: 201 (= id_HM80_ice).
 *
 *      do_output | int | Set 1 to write the output file of rho, u, P values,
 *          set 0 for no output. Default: 0.
 *
 * Output text file contains:
 *  header
 *  num_rho num_u   mat_id                      # Header values
 *  rho_0   rho_1   rho_2   ...   rho_num_rho   # Array of densities, rho
 *  u_0     u_1     u_2     ...   u_num_u       # Array of energies, u
 *  P_0_0   P_0_1   ...     P_0_num_u           # Array of pressures, P(rho, u)
 *  P_1_0   ...     ...     P_1_num_u
 *  ...     ...     ...     ...
 *  P_num_rho_0     ...     P_num_rho_num_u
77 78 79 80 81
 *  c_0_0   c_0_1   ...     c_0_num_u           # Array of sound speeds, c(rho,
 * u)
 *  c_1_0   ...     ...     c_1_num_u
 *  ...     ...     ...     ...
 *  c_num_rho_0     ...     c_num_rho_num_u
82 83 84 85 86 87 88 89 90
 *
 * Note that the values tested extend beyond the range that most EOS are
 * designed for (e.g. outside table limits), to help test the EOS in case of
 * unexpected particle behaviour.
 *
 */

#ifdef EOS_PLANETARY
int main(int argc, char *argv[]) {
91
  float rho, u, log_rho, log_u, P, c;
92
  struct unit_system us;
93 94 95
  struct swift_params *params =
      (struct swift_params *)malloc(sizeof(struct swift_params));
  if (params == NULL) error("Error allocating memory for the parameter file.");
96
  const struct phys_const *phys_const = 0;  // Unused placeholder
97
  const float J_kg_to_erg_g = 1e4;          // Convert J/kg to erg/g
98 99 100
  char filename[64];
  // Output table params
  const int num_rho = 100, num_u = 100;
101 102 103 104
  float log_rho_min = logf(1e-4f), log_rho_max = logf(1e3f),  // Densities (cgs)
      log_u_min = logf(1e4f),
        log_u_max = logf(1e13f),  // Sp. int. energies (SI)
      log_rho_step = (log_rho_max - log_rho_min) / (num_rho - 1.f),
105 106 107
        log_u_step = (log_u_max - log_u_min) / (num_u - 1.f);
  float A1_rho[num_rho], A1_u[num_u];
  // Sys args
108
  int mat_id_in, do_output;
109
  // Default sys args
110
  const int mat_id_def = eos_planetary_id_HM80_ice;
111 112 113 114 115 116
  const int do_output_def = 0;

  // Check the number of system arguments and use defaults if not provided
  switch (argc) {
    case 1:
      // Default both
117
      mat_id_in = mat_id_def;
118 119 120 121 122
      do_output = do_output_def;
      break;

    case 2:
      // Read mat_id, default do_output
123
      mat_id_in = atoi(argv[1]);
124 125 126 127 128
      do_output = do_output_def;
      break;

    case 3:
      // Read both
129
      mat_id_in = atoi(argv[1]);
130 131 132 133 134
      do_output = atoi(argv[2]);
      break;

    default:
      error("Invalid number of system arguments!\n");
135
      mat_id_in = mat_id_def;  // Ignored, just here to keep the compiler happy
136 137 138
      do_output = do_output_def;
  };

139 140 141
  enum eos_planetary_material_id mat_id =
      (enum eos_planetary_material_id)mat_id_in;

142 143 144 145
  /* Greeting message */
  printf("This is %s\n", package_description());

  // Check material ID
146 147 148 149 150
  const enum eos_planetary_type_id type =
      (enum eos_planetary_type_id)(mat_id / eos_planetary_type_factor);

  // Select the material base type
  switch (type) {
151
    // Tillotson
152
    case eos_planetary_type_Til:
153
      switch (mat_id) {
154
        case eos_planetary_id_Til_iron:
155 156 157
          printf("  Tillotson iron \n");
          break;

158
        case eos_planetary_id_Til_granite:
159 160 161
          printf("  Tillotson granite \n");
          break;

162
        case eos_planetary_id_Til_water:
163 164 165 166 167 168 169 170 171
          printf("  Tillotson water \n");
          break;

        default:
          error("Unknown material ID! mat_id = %d \n", mat_id);
      };
      break;

    // Hubbard & MacFarlane (1980)
172
    case eos_planetary_type_HM80:
173
      switch (mat_id) {
174
        case eos_planetary_id_HM80_HHe:
175 176 177
          printf("  Hubbard & MacFarlane (1980) hydrogen-helium atmosphere \n");
          break;

178
        case eos_planetary_id_HM80_ice:
179 180 181
          printf("  Hubbard & MacFarlane (1980) ice mix \n");
          break;

182
        case eos_planetary_id_HM80_rock:
183 184 185 186 187 188 189 190
          printf("  Hubbard & MacFarlane (1980) rock mix \n");
          break;

        default:
          error("Unknown material ID! mat_id = %d \n", mat_id);
      };
      break;

191 192
    // SESAME
    case eos_planetary_type_SESAME:
193
      switch (mat_id) {
194 195
        case eos_planetary_id_SESAME_iron:
          printf("  SESAME basalt 7530 \n");
196 197
          break;

198 199
        case eos_planetary_id_SESAME_basalt:
          printf("  SESAME basalt 7530 \n");
200 201
          break;

202 203 204
        case eos_planetary_id_SESAME_water:
          printf("  SESAME water 7154 \n");
          break;
205

206 207
        case eos_planetary_id_SS08_water:
          printf("  Senft & Stewart (2008) SESAME-like water \n");
208 209 210 211 212 213 214 215 216 217 218
          break;

        default:
          error("Unknown material ID! mat_id = %d \n", mat_id);
      };
      break;

    default:
      error("Unknown material type! mat_id = %d \n", mat_id);
  }

219 220 221 222 223
  // Convert to internal units
  // Earth masses and radii
  //  units_init(&us, 5.9724e27, 6.3710e8, 1.f, 1.f, 1.f);
  // SI
  units_init(&us, 1000.f, 100.f, 1.f, 1.f, 1.f);
224 225
  log_rho_min -= logf(units_cgs_conversion_factor(&us, UNIT_CONV_DENSITY));
  log_rho_max -= logf(units_cgs_conversion_factor(&us, UNIT_CONV_DENSITY));
226 227 228 229
  log_u_min += logf(J_kg_to_erg_g / units_cgs_conversion_factor(
                                        &us, UNIT_CONV_ENERGY_PER_UNIT_MASS));
  log_u_max += logf(J_kg_to_erg_g / units_cgs_conversion_factor(
                                        &us, UNIT_CONV_ENERGY_PER_UNIT_MASS));
230

231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258
  // Set the input parameters
  // Which EOS to initialise
  parser_set_param(params, "EoS:planetary_use_Til:1");
  parser_set_param(params, "EoS:planetary_use_HM80:1");
  parser_set_param(params, "EoS:planetary_use_SESAME:1");
  // Table file names
  parser_set_param(params,
                   "EoS:planetary_HM80_HHe_table_file:"
                   "../examples/planetary_HM80_HHe.txt");
  parser_set_param(params,
                   "EoS:planetary_HM80_ice_table_file:"
                   "../examples/planetary_HM80_ice.txt");
  parser_set_param(params,
                   "EoS:planetary_HM80_rock_table_file:"
                   "../examples/planetary_HM80_rock.txt");
  parser_set_param(params,
                   "EoS:planetary_SESAME_iron_table_file:"
                   "../examples/planetary_SESAME_iron_2140.txt");
  parser_set_param(params,
                   "EoS:planetary_SESAME_basalt_table_file:"
                   "../examples/planetary_SESAME_basalt_7530.txt");
  parser_set_param(params,
                   "EoS:planetary_SESAME_water_table_file:"
                   "../examples/planetary_SESAME_water_7154.txt");
  parser_set_param(params,
                   "EoS:planetary_SS08_water_table_file:"
                   "../examples/planetary_SS08_water.txt");

259 260 261
  // Initialise the EOS materials
  eos_init(&eos, phys_const, &us, params);

262 263 264 265 266 267 268 269 270 271 272 273
  // Manual debug testing
  if (1) {
    printf("\n ### MANUAL DEBUG TESTING ### \n");

    rho = 5960;
    u = 1.7e8;
    P = gas_pressure_from_internal_energy(rho, u, eos_planetary_id_HM80_ice);
    printf("u = %.2e,    rho = %.2e,    P = %.2e \n", u, rho, P);

    return 0;
  }

274
  // Output file
275
  sprintf(filename, "testEOS_rho_u_P_c_%d.txt", mat_id);
276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305
  FILE *f = fopen(filename, "w");
  if (f == NULL) {
    printf("Could not open output file!\n");
    exit(EXIT_FAILURE);
  }

  if (do_output == 1) {
    fprintf(f, "Density  Sp.Int.Energy  mat_id \n");
    fprintf(f, "%d      %d            %d \n", num_rho, num_u, mat_id);
  }

  // Densities
  log_rho = log_rho_min;
  for (int i = 0; i < num_rho; i++) {
    A1_rho[i] = exp(log_rho);
    log_rho += log_rho_step;

    if (do_output == 1)
      fprintf(f, "%.6e ",
              A1_rho[i] * units_cgs_conversion_factor(&us, UNIT_CONV_DENSITY));
  }
  if (do_output == 1) fprintf(f, "\n");

  // Sp. int. energies
  log_u = log_u_min;
  for (int i = 0; i < num_u; i++) {
    A1_u[i] = exp(log_u);
    log_u += log_u_step;

    if (do_output == 1)
306 307 308
      fprintf(f, "%.6e ",
              A1_u[i] * units_cgs_conversion_factor(
                            &us, UNIT_CONV_ENERGY_PER_UNIT_MASS));
309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325
  }
  if (do_output == 1) fprintf(f, "\n");

  // Pressures
  for (int i = 0; i < num_rho; i++) {
    rho = A1_rho[i];

    for (int j = 0; j < num_u; j++) {
      P = gas_pressure_from_internal_energy(rho, A1_u[j], mat_id);

      if (do_output == 1)
        fprintf(f, "%.6e ",
                P * units_cgs_conversion_factor(&us, UNIT_CONV_PRESSURE));
    }

    if (do_output == 1) fprintf(f, "\n");
  }
326 327 328 329 330 331 332 333 334 335 336 337 338 339 340

  // Sound speeds
  for (int i = 0; i < num_rho; i++) {
    rho = A1_rho[i];

    for (int j = 0; j < num_u; j++) {
      c = gas_soundspeed_from_internal_energy(rho, A1_u[j], mat_id);

      if (do_output == 1)
        fprintf(f, "%.6e ",
                c * units_cgs_conversion_factor(&us, UNIT_CONV_SPEED));
    }

    if (do_output == 1) fprintf(f, "\n");
  }
341 342 343 344 345
  fclose(f);

  return 0;
}
#else
346
int main(int argc, char *argv[]) { return 0; }
347
#endif