/******************************************************************************* * 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 . * ******************************************************************************/ #include /* Some standard headers. */ #include #include /* Local headers. */ #include "swift.h" #if defined(COOLING_COLIBRE) && defined(GADGET2_SPH) #include "cooling/COLIBRE/cooling_rates.h" #include "cooling/COLIBRE/cooling_tables.h" /** * @brief Assign particle density and entropy corresponding to the * hydrogen number density and internal energy specified. * * @param p Particle data structure * @param cooling Cooling function data structure * @param cosmo Cosmology data structure * @param internal_const Physical constants data structure * @param nh Hydrogen number density (cgs units) * @param u Internal energy (cgs units) */ void set_quantities(struct part *restrict p, struct xpart *restrict xp, const struct unit_system *restrict us, const struct cooling_function_data *restrict cooling, const struct cosmology *restrict cosmo, const struct phys_const *restrict internal_const, float nh, double u) { double hydrogen_number_density = nh * pow(units_cgs_conversion_factor(us, UNIT_CONV_LENGTH), 3); p->rho = hydrogen_number_density * internal_const->const_proton_mass / p->chemistry_data.metal_mass_fraction[chemistry_element_H]; float pressure = (u * cosmo->a * cosmo->a) * cooling->internal_energy_from_cgs * p->rho * (hydro_gamma_minus_one); p->entropy = pressure * (pow(p->rho, -hydro_gamma)); xp->entropy_full = p->entropy; } /** * @brief Produces contributions to cooling rates for different * hydrogen number densities, from different metals, * tests 1d and 4d table interpolations produce * same results for cooling rate, dlambda/du and temperature. */ int main(int argc, char **argv) { // Declare relevant structs struct swift_params *params = malloc(sizeof(struct swift_params)); struct unit_system us; struct dustevo_props dustevo_properties; struct chemistry_global_data chem_data; struct part p; struct xpart xp; struct phys_const internal_const; struct hydro_props hydro_properties; struct cooling_function_data cooling; struct cosmology cosmo; const char *parametersFileName = "./cooling_rates.yml"; /* Initialize CPU frequency, this also starts time. */ unsigned long long cpufreq = 0; clocks_set_cpufreq(cpufreq); /* Choke on FP-exceptions */ #ifdef HAVE_FE_ENABLE_EXCEPT feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW); #endif const int npts = 250; // number of values for the internal energy at which // cooling rate is evaluated // Set some default values float redshift = 0.0, log_10_nh = -1; // Read options int param; while ((param = getopt(argc, argv, "z:d:")) != -1) switch (param) { case 'z': // read redshift redshift = atof(optarg); break; case 'd': // read log10 of hydrogen number density log_10_nh = atof(optarg); break; case '?': if (optopt == 'z') printf("Option -%c requires an argument.\n", optopt); else printf("Unknown option character `\\x%x'.\n", optopt); error("invalid option(s) to cooling_rates"); } // 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, &internal_const); // Init properties of hydro hydro_props_init(&hydro_properties, &internal_const, &us, params); // Init chemistry chemistry_init(params, &us, &internal_const, &chem_data); chemistry_first_init_part(&internal_const, &us, &cosmo, &chem_data, &p, &xp); chemistry_print(&chem_data); // Init cosmology cosmology_init(params, &us, &internal_const, &cosmo); /* Initialise the dust evolution structure */ bzero(&dustevo_properties, sizeof(struct dustevo_props)); // Set redshift and associated quantities const float scale_factor = 1.0 / (1.0 + redshift); integertime_t ti_current = log(scale_factor / cosmo.a_begin) / cosmo.time_base; cosmology_update(&cosmo, &internal_const, ti_current); message("Redshift is %f", cosmo.z); message("Density is [log nH]%f", log_10_nh); // Init cooling cooling_init(params, &us, &internal_const, &hydro_properties, &cooling, &dustevo_properties); cooling_print(&cooling); // extract mass fractions, calculate table indices and offsets float XH = p.chemistry_data.metal_mass_fraction[chemistry_element_H]; char cool_outputfile[35]; sprintf(cool_outputfile, "cooling_output_lognH%.2f.dat", log_10_nh); char heat_outputfile[35]; sprintf(heat_outputfile, "heating_output_lognH%.2f.dat", log_10_nh); char net_outputfile[35]; sprintf(net_outputfile, "netcooling_output_lognH%.2f.dat", log_10_nh); // Calculate contributions from metals to cooling rate // open file FILE *output_file = fopen(cool_outputfile, "w"); if (output_file == NULL) { error("Error opening output file!\n"); } FILE *output_file_heat = fopen(heat_outputfile, "w"); if (output_file_heat == NULL) { error("Error opening output file (heat)!\n"); } FILE *output_file_net = fopen(net_outputfile, "w"); if (output_file_net == NULL) { error("Error opening output file (net)!\n"); } // set hydrogen number density const float nh = exp(M_LN10 * log_10_nh); /* Initial internal energy */ double u = 1.0e14; // set internal energy to dummy value, will get reset when looping over // internal energies set_quantities(&p, &xp, &us, &cooling, &cosmo, &internal_const, nh, u); float abundance_ratio[chemistry_element_count + 3]; float logZZsol = abundance_ratio_to_solar(&p, &cooling, abundance_ratio); float inn_h = hydro_get_physical_density(&p, &cosmo) * XH / internal_const.const_proton_mass * cooling.number_density_to_cgs; float d_red, d_met, d_n_H; int red_index, met_index, n_H_index; double lambda_net, temperature; if (cosmo.z < cooling.H_reion_z) { get_index_1d(cooling.Redshifts, colibre_cooling_N_redshifts, cosmo.z, &red_index, &d_red); } else { red_index = colibre_cooling_N_redshifts - 2; d_red = 1.0; } message("cooling.Redshifts[red_index] = %.4f", cooling.Redshifts[red_index]); message("red_index = %i", red_index); message("d_red = %.4f", d_red); get_index_1d(cooling.Metallicity, colibre_cooling_N_metallicity, logZZsol, &met_index, &d_met); get_index_1d(cooling.nH, colibre_cooling_N_density, log10(inn_h), &n_H_index, &d_n_H); // Loop over internal energy for (int j = 0; j < npts; j++) { // Update the particle with the new values set_quantities(&p, &xp, &us, &cooling, &cosmo, &internal_const, nh, pow(10.0, 9.0 + j * 9.5 / npts)); // New internal energy u = hydro_get_physical_internal_energy(&p, &xp, &cosmo) * cooling.internal_energy_to_cgs; // calculate temperature temperature = colibre_convert_u_to_temp(log10(u), cosmo.z, n_H_index, d_n_H, met_index, d_met, red_index, d_red, &cooling); // calculate cooling rates lambda_net = colibre_cooling_rate(log10(u), cosmo.z, nh, abundance_ratio, n_H_index, d_n_H, met_index, d_met, red_index, d_red, &cooling, 0, 1, 0, -1); fprintf(output_file, "%.5e %.5e", exp(M_LN10 * temperature), lambda_net); for (int icool = 0; icool < colibre_cooling_N_cooltypes - 2; icool++) { lambda_net = colibre_cooling_rate( log10(u), cosmo.z, nh, abundance_ratio, n_H_index, d_n_H, met_index, d_met, red_index, d_red, &cooling, 1, 1, icool, -1); fprintf(output_file, " %.5e", lambda_net); } fprintf(output_file, "\n"); // calculate heating rates lambda_net = colibre_cooling_rate(log10(u), cosmo.z, nh, abundance_ratio, n_H_index, d_n_H, met_index, d_met, red_index, d_red, &cooling, 1, 0, -1, 0); fprintf(output_file_heat, "%.5e %.5e", exp(M_LN10 * temperature), lambda_net); for (int iheat = 0; iheat < colibre_cooling_N_heattypes - 2; iheat++) { lambda_net = colibre_cooling_rate( log10(u), cosmo.z, nh, abundance_ratio, n_H_index, d_n_H, met_index, d_met, red_index, d_red, &cooling, 1, 1, -1, iheat); fprintf(output_file_heat, " %.5e", lambda_net); } fprintf(output_file_heat, "\n"); // calculate standard net cooling lambda_net = colibre_cooling_rate(log10(u), cosmo.z, nh, abundance_ratio, n_H_index, d_n_H, met_index, d_met, red_index, d_red, &cooling, 0, 0, 0, 0); fprintf(output_file_net, "%.5e %.5e\n", exp(M_LN10 * temperature), lambda_net); } fclose(output_file); fclose(output_file_heat); fclose(output_file_net); /* Clean everything */ cosmology_clean(&cosmo); cooling_clean(&cooling); free(params); return 0; } #else int main(int argc, char **argv) { /* Initialize CPU frequency, this also starts time. */ unsigned long long cpufreq = 0; clocks_set_cpufreq(cpufreq); message( "This test is only defined for the COLIBRE cooling model and Gadget-2 " "SPH."); return 0; } #endif