testCooling.c 9.22 KB
Newer Older
Alexei Borissov's avatar
Alexei Borissov committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
/*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (C) 2015 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/>.
 *
 ******************************************************************************/
19
#include "../config.h"
Alexei Borissov's avatar
Alexei Borissov committed
20

21
/* Local headers. */
22
#include "swift.h"
Alexei Borissov's avatar
Alexei Borissov committed
23

24
#if defined(CHEMISTRY_EAGLE) && defined(COOLING_EAGLE) && defined(GADGET2_SPH)
Matthieu Schaller's avatar
Matthieu Schaller committed
25

Alexei Borissov's avatar
Alexei Borissov committed
26 27 28 29 30 31 32 33 34 35 36 37 38 39
/*
 * @brief Assign particle density and entropy corresponding to the
 * hydrogen number density and internal energy specified.
 *
 * @param p Particle data structure
 * @param xp extra particle structure
 * @param us unit system struct
 * @param cooling Cooling function data structure
 * @param cosmo Cosmology data structure
 * @param phys_const Physical constants data structure
 * @param nh_cgs Hydrogen number density (cgs units)
 * @param u Internal energy (cgs units)
 * @param ti_current integertime to set cosmo quantities
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
40
void set_quantities(struct part *restrict p, struct xpart *restrict xp,
Alexei Borissov's avatar
Alexei Borissov committed
41 42 43 44
                    const struct unit_system *restrict us,
                    const struct cooling_function_data *restrict cooling,
                    struct cosmology *restrict cosmo,
                    const struct phys_const *restrict phys_const, float nh_cgs,
45
                    double u_cgs, integertime_t ti_current) {
Alexei Borissov's avatar
Alexei Borissov committed
46
  /* calculate density */
Alexei Borissov's avatar
Alexei Borissov committed
47
  double hydrogen_number_density = nh_cgs / cooling->number_density_to_cgs;
Alexei Borissov's avatar
Alexei Borissov committed
48
  p->rho = hydrogen_number_density * phys_const->const_proton_mass /
49
           p->chemistry_data.metal_mass_fraction[chemistry_element_H];
Alexei Borissov's avatar
Alexei Borissov committed
50 51

  /* update entropy based on internal energy */
52 53
  float pressure = (u_cgs)*cooling->internal_energy_from_cgs * p->rho *
                   (hydro_gamma_minus_one);
Alexei Borissov's avatar
Alexei Borissov committed
54 55
  p->entropy = pressure * (pow(p->rho, -hydro_gamma));
  xp->entropy_full = p->entropy;
56 57

  p->entropy_dt = 0.f;
Alexei Borissov's avatar
Alexei Borissov committed
58 59 60
}

/*
Alexei Borissov's avatar
Alexei Borissov committed
61 62
 * @brief Tests cooling integration scheme by comparing EAGLE
 * integration to subcycled explicit equation.
Alexei Borissov's avatar
Alexei Borissov committed
63 64 65 66 67 68 69 70 71 72 73 74 75
 */
int main(int argc, char **argv) {
  // Declare relevant structs
  struct swift_params *params = malloc(sizeof(struct swift_params));
  struct unit_system us;
  struct chemistry_global_data chem_data;
  struct part p;
  struct xpart xp;
  struct phys_const phys_const;
  struct cooling_function_data cooling;
  struct cosmology cosmo;
  char *parametersFileName = "./testCooling.yml";

76 77 78 79
  float nh_cgs;  // hydrogen number density
  double u_cgs;  // internal energy

  const float seconds_per_year = 3.154e7;
Alexei Borissov's avatar
Alexei Borissov committed
80

Matthieu Schaller's avatar
Matthieu Schaller committed
81
  /* Number of values to test for in redshift,
Alexei Borissov's avatar
Alexei Borissov committed
82
   * hydrogen number density and internal energy */
83 84 85
  const int n_z = 10;
  const int n_nh = 10;
  const int n_u = 10;
Alexei Borissov's avatar
Alexei Borissov committed
86 87

  /* Number of subcycles and tolerance used to compare
Matthieu Schaller's avatar
Matthieu Schaller committed
88
   * subcycled and implicit solution. Note, high value
Alexei Borissov's avatar
Alexei Borissov committed
89 90
   * of tolerance due to mismatch between explicit and
   * implicit solution for large timesteps */
91
  const int n_subcycle = 1000;
Alexei Borissov's avatar
Alexei Borissov committed
92 93 94 95 96 97 98 99 100 101 102 103 104 105
  const float integration_tolerance = 0.2;

  /* Read the parameter file */
  if (params == NULL) error("Error allocating memory for the parameter file.");
  message("Reading runtime parameters from file '%s'", parametersFileName);
  parser_read_file(parametersFileName, params);

  /* Init units */
  units_init_from_params(&us, params, "InternalUnitSystem");
  phys_const_init(&us, params, &phys_const);

  /* Init chemistry */
  chemistry_init(params, &us, &phys_const, &chem_data);
  chemistry_first_init_part(&phys_const, &us, &cosmo, &chem_data, &p, &xp);
106
  chemistry_part_has_no_neighbours(&p, &xp, &chem_data, &cosmo);
Alexei Borissov's avatar
Alexei Borissov committed
107 108 109 110 111
  chemistry_print(&chem_data);

  /* Init cosmology */
  cosmology_init(params, &us, &phys_const, &cosmo);
  cosmology_print(&cosmo);
112 113 114 115 116

  /* Set dt */
  const int timebin = 38;
  float dt_cool, dt_therm;

Alexei Borissov's avatar
Alexei Borissov committed
117 118 119
  /* Init hydro_props */
  struct hydro_props hydro_properties;
  hydro_props_init(&hydro_properties, &phys_const, &us, params);
Alexei Borissov's avatar
Alexei Borissov committed
120 121

  /* Init cooling */
Alexei Borissov's avatar
Alexei Borissov committed
122
  cooling_init(params, &us, &phys_const, &hydro_properties, &cooling);
Alexei Borissov's avatar
Alexei Borissov committed
123 124 125
  cooling_print(&cooling);
  cooling_update(&cosmo, &cooling, 0);

Alexei Borissov's avatar
Alexei Borissov committed
126 127 128
  /* Init entropy floor */
  struct entropy_floor_properties floor_props;
  entropy_floor_init(&floor_props, &phys_const, &us, &hydro_properties, params);
129

Alexei Borissov's avatar
Alexei Borissov committed
130 131 132 133 134
  /* Cooling function needs to know the minimal energy. Set it to the lowest
   * internal energy in the cooling table. */
  hydro_properties.minimal_internal_energy =
      exp(M_LN10 * cooling.Therm[0]) * cooling.internal_energy_from_cgs;

Alexei Borissov's avatar
Alexei Borissov committed
135 136 137 138 139 140 141 142 143 144 145 146
  /* Calculate abundance ratios */
  float *abundance_ratio;
  abundance_ratio = malloc((chemistry_element_count + 2) * sizeof(float));
  abundance_ratio_to_solar(&p, &cooling, abundance_ratio);

  /* extract mass fractions, calculate table indices and offsets */
  float XH = p.chemistry_data.metal_mass_fraction[chemistry_element_H];
  float HeFrac =
      p.chemistry_data.metal_mass_fraction[chemistry_element_He] /
      (XH + p.chemistry_data.metal_mass_fraction[chemistry_element_He]);
  int He_i;
  float d_He;
Alexei Borissov's avatar
Alexei Borissov committed
147
  get_index_1d(cooling.HeFrac, eagle_cooling_N_He_frac, HeFrac, &He_i, &d_He);
Alexei Borissov's avatar
Alexei Borissov committed
148

Matthieu Schaller's avatar
Matthieu Schaller committed
149
  /* calculate spacing in nh and u */
150 151 152 153
  const float log_u_min_cgs = 11, log_u_max_cgs = 17;
  const float log_nh_min_cgs = -6, log_nh_max_cgs = 3;
  const float delta_log_nh_cgs = (log_nh_max_cgs - log_nh_min_cgs) / n_nh;
  const float delta_log_u_cgs = (log_u_max_cgs - log_u_min_cgs) / n_u;
Alexei Borissov's avatar
Alexei Borissov committed
154 155

  /* Declare variables we will be checking */
156
  double du_dt_implicit, du_dt_check;
Alexei Borissov's avatar
Alexei Borissov committed
157
  integertime_t ti_current;
Matthieu Schaller's avatar
Matthieu Schaller committed
158

159
  /* Loop over values of nh and u */
Alexei Borissov's avatar
Alexei Borissov committed
160
  for (int nh_i = 0; nh_i < n_nh; nh_i++) {
161
    nh_cgs = exp(M_LN10 * log_nh_min_cgs + delta_log_nh_cgs * nh_i);
Alexei Borissov's avatar
Alexei Borissov committed
162
    for (int u_i = 0; u_i < n_u; u_i++) {
163
      u_cgs = exp(M_LN10 * log_u_min_cgs + delta_log_u_cgs * u_i);
Matthieu Schaller's avatar
Matthieu Schaller committed
164

165
      /* Loop over z */
166 167
      for (int z_i = 0; z_i <= n_z; z_i++) {
        ti_current = max_nr_timesteps / n_z * z_i + 1;
Alexei Borissov's avatar
Alexei Borissov committed
168

169 170 171 172
        /* update nh, u, z */
        cosmology_update(&cosmo, &phys_const, ti_current);
        cooling_init(params, &us, &phys_const, &hydro_properties, &cooling);
        cooling_update(&cosmo, &cooling, 0);
173 174
        set_quantities(&p, &xp, &us, &cooling, &cosmo, &phys_const, nh_cgs,
                       u_cgs, ti_current);
175 176 177 178 179 180 181 182 183 184 185 186 187 188

        /* Set dt */
        const integertime_t ti_step = get_integer_timestep(timebin);
        const integertime_t ti_begin =
            get_integer_time_begin(ti_current - 1, timebin);
        dt_cool =
            cosmology_get_delta_time(&cosmo, ti_begin, ti_begin + ti_step);
        dt_therm = cosmology_get_therm_kick_factor(&cosmo, ti_begin,
                                                   ti_begin + ti_step);

        /* calculate subcycled solution */
        for (int t_subcycle = 0; t_subcycle < n_subcycle; t_subcycle++) {
          p.entropy_dt = 0;
          cooling_cool_part(&phys_const, &us, &cosmo, &hydro_properties,
189 190
                            &floor_props, &cooling, &p, &xp,
                            dt_cool / n_subcycle, dt_therm / n_subcycle);
191
          xp.entropy_full += p.entropy_dt * dt_therm / n_subcycle;
Alexei Borissov's avatar
Alexei Borissov committed
192
        }
193 194
        du_dt_check = hydro_get_physical_internal_energy_dt(&p, &cosmo);

195 196
        /* reset quantities to nh, u, and z that we want to test */
        cosmology_update(&cosmo, &phys_const, ti_current);
197 198 199
        set_quantities(&p, &xp, &us, &cooling, &cosmo, &phys_const, nh_cgs,
                       u_cgs, ti_current);

200
        /* compute implicit solution */
201 202 203
        cooling_cool_part(&phys_const, &us, &cosmo, &hydro_properties,
                          &floor_props, &cooling, &p, &xp, dt_cool, dt_therm);
        du_dt_implicit = hydro_get_physical_internal_energy_dt(&p, &cosmo);
Matthieu Schaller's avatar
Matthieu Schaller committed
204 205

        /* check if the two solutions are consistent */
206
        if (fabs((du_dt_implicit - du_dt_check) / du_dt_check) >
207
                integration_tolerance ||
208
            (du_dt_check == 0.0 && du_dt_implicit != 0.0))
209
          error(
210
              "Solutions do not match. scale factor %.5e z %.5e nh_cgs %.5e "
211 212 213 214 215 216 217 218 219 220 221 222 223
              "u_cgs %.5e dt (years) %.5e du cgs implicit %.5e reference %.5e "
              "error %.5e",
              cosmo.a, cosmo.z, nh_cgs, u_cgs,
              dt_cool * units_cgs_conversion_factor(&us, UNIT_CONV_TIME) /
                  seconds_per_year,
              du_dt_implicit *
                  units_cgs_conversion_factor(&us,
                                              UNIT_CONV_ENERGY_PER_UNIT_MASS) *
                  dt_therm,
              du_dt_check *
                  units_cgs_conversion_factor(&us,
                                              UNIT_CONV_ENERGY_PER_UNIT_MASS) *
                  dt_therm,
224
              fabs((du_dt_implicit - du_dt_check) / du_dt_check));
Alexei Borissov's avatar
Alexei Borissov committed
225 226 227
      }
    }
  }
228
  message("done explicit subcycling cooling test");
Alexei Borissov's avatar
Alexei Borissov committed
229 230 231 232

  free(params);
  return 0;
}
Matthieu Schaller's avatar
Matthieu Schaller committed
233 234 235 236 237 238

#else

int main(int argc, char **argv) { return 0; }

#endif