/******************************************************************************* * This file is part of SWIFT. * Copyright (c) 2020 Mladen Ivkovic (mladen.ivkovic@hotmail.com) * * 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 . * ******************************************************************************/ #ifndef SWIFT_RT_IO_GEAR_H #define SWIFT_RT_IO_GEAR_H #define RT_LABELS_SIZE 10 /** * @file src/rt/GEAR/rt_io.h * @brief Main header file for GEAR M1 Closure radiative transfer * scheme IO routines. */ /** * @brief Specifies which particle fields to read from a dataset * * @param parts The particle array. * @param list The list of i/o properties to read. * * @return Returns the number of fields to read. */ INLINE static int rt_read_particles(const struct part* parts, struct io_props* list) { /* List what we want to read */ char fieldname[30]; int count = 0; for (int phg = 0; phg < RT_NGROUPS; phg++) { sprintf(fieldname, "PhotonEnergiesGroup%d", phg + 1); list[count++] = io_make_input_field(fieldname, FLOAT, 1, OPTIONAL, UNIT_CONV_ENERGY, parts, rt_data.radiation[phg].energy_density); sprintf(fieldname, "PhotonFluxesGroup%d", phg + 1); list[count++] = io_make_input_field(fieldname, FLOAT, 3, OPTIONAL, UNIT_CONV_ENERGY_VELOCITY, parts, rt_data.radiation[phg].flux); } list[count++] = io_make_input_field("MassFractionHI", FLOAT, 1, OPTIONAL, UNIT_CONV_NO_UNITS, parts, rt_data.tchem.mass_fraction_HI); list[count++] = io_make_input_field("MassFractionHII", FLOAT, 1, OPTIONAL, UNIT_CONV_NO_UNITS, parts, rt_data.tchem.mass_fraction_HII); list[count++] = io_make_input_field("MassFractionHeI", FLOAT, 1, OPTIONAL, UNIT_CONV_NO_UNITS, parts, rt_data.tchem.mass_fraction_HeI); list[count++] = io_make_input_field("MassFractionHeII", FLOAT, 1, OPTIONAL, UNIT_CONV_NO_UNITS, parts, rt_data.tchem.mass_fraction_HeII); list[count++] = io_make_input_field("MassFractionHeIII", FLOAT, 1, OPTIONAL, UNIT_CONV_NO_UNITS, parts, rt_data.tchem.mass_fraction_HeIII); return count; } /** * @brief Specifies which particle fields to read from a dataset * * @param sparts The star particle array. * @param list The list of i/o properties to read. * * @return Returns the number of fields to read. */ INLINE static int rt_read_stars(const struct spart* sparts, struct io_props* list) { return 0; } /** * @brief Extract radiation energies of radiation struct for all photon groups * Note: "allocation" of `float* ret` happens in io_copy_temp_buffer() * * @param engine the engine * @param part the particle to extract data from * @param xpart the according xpart to extract data from * @param ret (return) the extracted data */ INLINE static void rt_convert_radiation_energies(const struct engine* engine, const struct part* part, const struct xpart* xpart, float* ret) { for (int g = 0; g < RT_NGROUPS; g++) { ret[g] = part->rt_data.radiation[g].energy_density * part->geometry.volume; } } /** * @brief Extract radiation fluxes of radiation struct for all photon groups * Note: "allocation" of `float* ret` happens in io_copy_temp_buffer() * * @param engine the engine * @param part the particle to extract data from * @param xpart the according xpart to extract data from * @param ret (return) the extracted data */ INLINE static void rt_convert_radiation_fluxes(const struct engine* engine, const struct part* part, const struct xpart* xpart, float* ret) { int i = 0; for (int g = 0; g < RT_NGROUPS; g++) { ret[i++] = part->rt_data.radiation[g].flux[0]; ret[i++] = part->rt_data.radiation[g].flux[1]; ret[i++] = part->rt_data.radiation[g].flux[2]; } } /** * @brief Extract mass fractions of constituent species from tchem struct. * Note: "allocation" of `float* ret` happens in io_copy_temp_buffer() * * @param engine the engine * @param part the particle to extract data from * @param xpart the according xpart to extract data from * @param ret (return) the extracted data */ INLINE static void rt_convert_mass_fractions(const struct engine* engine, const struct part* part, const struct xpart* xpart, float* ret) { ret[0] = part->rt_data.tchem.mass_fraction_HI; ret[1] = part->rt_data.tchem.mass_fraction_HII; ret[2] = part->rt_data.tchem.mass_fraction_HeI; ret[3] = part->rt_data.tchem.mass_fraction_HeII; ret[4] = part->rt_data.tchem.mass_fraction_HeIII; } /** * @brief Creates additional output fields for the radiative * transfer data of hydro particles. * * @param parts The particle array. * @param list The list of i/o properties to write. * * @return Returns the number of fields to write. */ INLINE static int rt_write_particles(const struct part* parts, struct io_props* list) { int num_elements = 3; list[0] = io_make_physical_output_field_convert_part( "PhotonEnergies", FLOAT, RT_NGROUPS, UNIT_CONV_ENERGY, 0, parts, /*xparts=*/NULL, /*convertible to comoving=*/1, rt_convert_radiation_energies, "Photon Energies (all groups)"); list[1] = io_make_physical_output_field_convert_part( "PhotonFluxes", FLOAT, 3 * RT_NGROUPS, UNIT_CONV_RADIATION_FLUX, 0, parts, /*xparts=*/NULL, /*convertible to comoving=*/1, rt_convert_radiation_fluxes, "Photon Fluxes (all groups; x, y, and z coordinates)"); list[2] = io_make_output_field_convert_part( "IonMassFractions", FLOAT, 5, UNIT_CONV_NO_UNITS, 0, parts, /*xparts=*/NULL, rt_convert_mass_fractions, "Mass fractions of all constituent species"); #ifdef SWIFT_RT_DEBUG_CHECKS num_elements += 8; list[3] = io_make_output_field("RTDebugInjectionDone", INT, 1, UNIT_CONV_NO_UNITS, 0, parts, rt_data.debug_injection_done, "How many times rt_injection_update_photon_density " "has been called"); list[4] = io_make_output_field( "RTDebugCallsIactGradientInteractions", INT, 1, UNIT_CONV_NO_UNITS, 0, parts, rt_data.debug_calls_iact_gradient_interaction, "number of calls to this particle during the gradient interaction loop " "from the actual interaction function"); list[5] = io_make_output_field("RTDebugCallsIactTransportInteractions", INT, 1, UNIT_CONV_NO_UNITS, 0, parts, rt_data.debug_calls_iact_transport_interaction, "number of calls to this particle during the " "transport interaction loop from the actual " "interaction function"); list[6] = io_make_output_field("RTDebugGradientsDone", INT, 1, UNIT_CONV_NO_UNITS, 0, parts, rt_data.debug_gradients_done, "How many times finalise_gradients was called"); list[7] = io_make_output_field("RTDebugTransportDone", INT, 1, UNIT_CONV_NO_UNITS, 0, parts, rt_data.debug_transport_done, "How many times finalise_transport was called"); list[8] = io_make_output_field( "RTDebugThermochemistryDone", INT, 1, UNIT_CONV_NO_UNITS, 0, parts, rt_data.debug_thermochem_done, "How many times rt_tchem was called"); list[9] = io_make_output_field( "RTDebugRadAbsorbedTot", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0, parts, rt_data.debug_radiation_absorbed_tot, "Radiation absorbed by this part during its lifetime"); list[10] = io_make_output_field( "RTDebugSubcycles", INT, 1, UNIT_CONV_NO_UNITS, 0, parts, rt_data.debug_nsubcycles, "How many times this part was subcycled"); #endif return num_elements; } /** * @brief Creates additional output fields for the radiative * transfer data of star particles. * * @param sparts The star particle array. * @param list The list of i/o properties to write. * * @return Returns the number of fields to write. */ INLINE static int rt_write_stars(const struct spart* sparts, struct io_props* list) { int num_elements = 0; #ifdef SWIFT_RT_DEBUG_CHECKS num_elements += 4; list[0] = io_make_output_field("RTDebugHydroIact", INT, 1, UNIT_CONV_NO_UNITS, 0, sparts, rt_data.debug_iact_hydro_inject, "number of interactions between this star " "particle and any particle during injection"); list[1] = io_make_output_field( "RTDebugEmissionRateSet", INT, 1, UNIT_CONV_NO_UNITS, 0, sparts, rt_data.debug_emission_rate_set, "Stellar photon emission rates set?"); list[2] = io_make_output_field( "RTDebugRadEmittedTot", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0, sparts, rt_data.debug_radiation_emitted_tot, "Total radiation emitted during the lifetime of this star"); list[3] = io_make_output_field("RTDebugInjectedPhotonEnergy", FLOAT, RT_NGROUPS, UNIT_CONV_ENERGY, 0, sparts, rt_data.debug_injected_energy_tot, "Total radiation actually injected into gas"); #endif return num_elements; } /** * @brief Write the RT model properties to the snapshot. * * @param h_grp The HDF5 group in which to write * @param h_grp_columns The HDF5 group containing named columns * @param e The engine * @param internal_units The internal unit system * @param snapshot_units Units used for the snapshot * @param rtp The #rt_props */ INLINE static void rt_write_flavour(hid_t h_grp, hid_t h_grp_columns, const struct engine* e, const struct unit_system* internal_units, const struct unit_system* snapshot_units, const struct rt_props* rtp) { #if defined(HAVE_HDF5) /* Write scheme name */ /* ----------------- */ io_write_attribute_s(h_grp, "RT Scheme", RT_IMPLEMENTATION); io_write_attribute_s(h_grp, "RT Riemann Solver", RT_RIEMANN_SOLVER_NAME); /* Write photon group counts */ /* ------------------------- */ io_write_attribute_i(h_grp, "PhotonGroupNumber", RT_NGROUPS); /* Write photon group bin edges */ /* ---------------------------- */ /* Note: photon frequency bin edges are kept in cgs. Convert them here to * internal units so we're still compatible with swiftsimio. */ const float Hz_internal = units_cgs_conversion_factor(internal_units, UNIT_CONV_INV_TIME); const float Hz_internal_inv = 1.f / Hz_internal; float photon_groups_internal[RT_NGROUPS]; for (int g = 0; g < RT_NGROUPS; g++) photon_groups_internal[g] = rtp->photon_groups[g] * Hz_internal_inv; hid_t type_float = H5Tcopy(io_hdf5_type(FLOAT)); hsize_t dims[1] = {RT_NGROUPS}; hid_t space = H5Screate_simple(1, dims, NULL); hid_t dset = H5Dcreate(h_grp, "PhotonGroupEdges", type_float, space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); H5Dwrite(dset, type_float, H5S_ALL, H5S_ALL, H5P_DEFAULT, photon_groups_internal); /* Write unit conversion factors for this data set */ char buffer[FIELD_BUFFER_SIZE] = {0}; units_cgs_conversion_string(buffer, snapshot_units, UNIT_CONV_INV_TIME, /*scale_factor_exponent=*/0); float baseUnitsExp[5]; units_get_base_unit_exponents_array(baseUnitsExp, UNIT_CONV_INV_TIME); io_write_attribute_f(dset, "U_M exponent", baseUnitsExp[UNIT_MASS]); io_write_attribute_f(dset, "U_L exponent", baseUnitsExp[UNIT_LENGTH]); io_write_attribute_f(dset, "U_t exponent", baseUnitsExp[UNIT_TIME]); io_write_attribute_f(dset, "U_I exponent", baseUnitsExp[UNIT_CURRENT]); io_write_attribute_f(dset, "U_T exponent", baseUnitsExp[UNIT_TEMPERATURE]); io_write_attribute_f(dset, "h-scale exponent", 0.f); io_write_attribute_f(dset, "a-scale exponent", 0.f); io_write_attribute_s(dset, "Expression for physical CGS units", buffer); io_write_attribute_b(dset, "Value stored as physical", 1); io_write_attribute_b(dset, "Property can be converted to comoving", 0); /* Write the actual number this conversion factor corresponds to */ const double factor = units_cgs_conversion_factor(snapshot_units, UNIT_CONV_INV_TIME); io_write_attribute_d( dset, "Conversion factor to CGS (not including cosmological corrections)", factor); io_write_attribute_d( dset, "Conversion factor to physical CGS (including cosmological corrections)", factor * pow(e->cosmology->a, 0.f)); H5Dclose(dset); /* H5Tclose(type_float); [> close this later <] */ /* If without RT, we have nothing more to do. */ const int with_rt = e->policy & engine_policy_rt; if (!with_rt) return; /* Write photon group names */ /* -------------------------*/ /* Generate Energy Group names */ char names_energy[RT_NGROUPS * RT_LABELS_SIZE]; for (int g = 0; g < RT_NGROUPS; g++) { char newEname[RT_LABELS_SIZE]; sprintf(newEname, "Group%d", g + 1); strcpy(names_energy + g * RT_LABELS_SIZE, newEname); } /* Now write them down */ hid_t type_string_label = H5Tcopy(H5T_C_S1); H5Tset_size(type_string_label, RT_LABELS_SIZE); hsize_t dimsE[1] = {RT_NGROUPS}; hid_t spaceE = H5Screate_simple(1, dimsE, NULL); hid_t dsetE = H5Dcreate(h_grp_columns, "PhotonEnergies", type_string_label, spaceE, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); H5Dwrite(dsetE, type_string_label, H5S_ALL, H5S_ALL, H5P_DEFAULT, names_energy); H5Dclose(dsetE); /* Generate Fluxes Group Names */ char names_fluxes[3 * RT_NGROUPS * RT_LABELS_SIZE]; int i = 0; for (int g = 0; g < RT_NGROUPS; g++) { char newFnameX[RT_LABELS_SIZE]; sprintf(newFnameX, "Group%dX", g + 1); strcpy(names_fluxes + i * RT_LABELS_SIZE, newFnameX); i++; char newFnameY[RT_LABELS_SIZE]; sprintf(newFnameY, "Group%dY", g + 1); strcpy(names_fluxes + i * RT_LABELS_SIZE, newFnameY); i++; char newFnameZ[RT_LABELS_SIZE]; sprintf(newFnameZ, "Group%dZ", g + 1); strcpy(names_fluxes + i * RT_LABELS_SIZE, newFnameZ); i++; } /* Now write them down */ hsize_t dimsF[1] = {3 * RT_NGROUPS}; hid_t spaceF = H5Screate_simple(1, dimsF, NULL); hid_t dsetF = H5Dcreate(h_grp_columns, "PhotonFluxes", type_string_label, spaceF, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); H5Dwrite(dsetF, type_string_label, H5S_ALL, H5S_ALL, H5P_DEFAULT, names_fluxes); H5Dclose(dsetF); /* H5Tclose(type_string_label); [> close this later <] */ /* Write reduced speed of light */ /* ---------------------------- */ /* hid_t type2 = H5Tcopy(io_hdf5_type(FLOAT)); */ hsize_t dims_cred[1] = {1}; hid_t space_cred = H5Screate_simple(1, dims_cred, NULL); hid_t dset_cred = H5Dcreate(h_grp, "ReducedLightspeed", type_float, space_cred, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); H5Dwrite(dset_cred, type_float, H5S_ALL, H5S_ALL, H5P_DEFAULT, &rt_params.reduced_speed_of_light); /* Write unit conversion factors for this data set */ char buffer_cred[FIELD_BUFFER_SIZE] = {0}; units_cgs_conversion_string(buffer_cred, snapshot_units, UNIT_CONV_VELOCITY, /*scale_factor_exponent=*/0); float baseUnitsExp_cred[5]; units_get_base_unit_exponents_array(baseUnitsExp_cred, UNIT_CONV_VELOCITY); io_write_attribute_f(dset_cred, "U_M exponent", baseUnitsExp_cred[UNIT_MASS]); io_write_attribute_f(dset_cred, "U_L exponent", baseUnitsExp_cred[UNIT_LENGTH]); io_write_attribute_f(dset_cred, "U_t exponent", baseUnitsExp_cred[UNIT_TIME]); io_write_attribute_f(dset_cred, "U_I exponent", baseUnitsExp_cred[UNIT_CURRENT]); io_write_attribute_f(dset_cred, "U_T exponent", baseUnitsExp_cred[UNIT_TEMPERATURE]); io_write_attribute_f(dset_cred, "h-scale exponent", 0.f); io_write_attribute_f(dset_cred, "a-scale exponent", 0.f); io_write_attribute_s(dset_cred, "Expression for physical CGS units", buffer_cred); io_write_attribute_b(dset_cred, "Value stored as physical", 1); io_write_attribute_b(dset_cred, "Property can be converted to comoving", 0); /* Write the actual number this conversion factor corresponds to */ /* TODO Mladen: check cosmology. reduced_speed_of_light is physical only for * now. */ const double factor_cred = units_cgs_conversion_factor(snapshot_units, UNIT_CONV_VELOCITY); io_write_attribute_d( dset_cred, "Conversion factor to CGS (not including cosmological corrections)", factor_cred); io_write_attribute_d( dset_cred, "Conversion factor to physical CGS (including cosmological corrections)", factor_cred * pow(e->cosmology->a, 0.f)); H5Dclose(dset_cred); /* Write constituent species mass fractions */ /* ---------------------------------------- */ char names_mf[5 * RT_LABELS_SIZE]; strcpy(names_mf + 0 * RT_LABELS_SIZE, "HI\0"); strcpy(names_mf + 1 * RT_LABELS_SIZE, "HII\0"); strcpy(names_mf + 2 * RT_LABELS_SIZE, "HeI\0"); strcpy(names_mf + 3 * RT_LABELS_SIZE, "HeII\0"); strcpy(names_mf + 4 * RT_LABELS_SIZE, "HeIII\0"); hsize_t dims_mf[1] = {5}; hid_t space_mf = H5Screate_simple(1, dims_mf, NULL); hid_t dset_mf = H5Dcreate(h_grp_columns, "IonMassFractions", type_string_label, space_mf, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); H5Dwrite(dset_mf, type_string_label, H5S_ALL, H5S_ALL, H5P_DEFAULT, names_mf); H5Dclose(dset_mf); /* Clean up after yourself */ /* ----------------------- */ /* Close up the types */ H5Tclose(type_float); H5Tclose(type_string_label); #endif /* HAVE_HDF5 */ } #endif /* SWIFT_RT_IO_GEAR_H */