diff --git a/.clang-format-ignore b/.clang-format-ignore index 36fa238b6..bda0eb1bb 100644 --- a/.clang-format-ignore +++ b/.clang-format-ignore @@ -46,6 +46,7 @@ src/clib/solve_chemistry.c src/clib/solve_rate_cool_g-cpp.cpp src/clib/status_reporting.c src/clib/status_reporting.h +src/clib/step_rate_gauss_seidel.hpp src/clib/step_rate_newton_raphson.hpp src/clib/time_deriv_0d.hpp src/clib/update_UVbackground_rates.c diff --git a/.clang-tidy b/.clang-tidy index 34b430c8a..dc5220ca0 100644 --- a/.clang-tidy +++ b/.clang-tidy @@ -40,12 +40,15 @@ Checks: > cert-flp30-c, cert-mem57-cpp, + clang-analyzer-*, + cppcoreguidelines-avoid-goto, cppcoreguidelines-avoid-non-const-global-variables, cppcoreguidelines-virtual-class-destructor google-build-using-namespace, google-explicit-constructor, + google-global-names-in-headers, google-readability-avoid-underscore-in-googletest-name, google-runtime-float, google-upgrade-googletest-case, diff --git a/src/clib/CMakeLists.txt b/src/clib/CMakeLists.txt index 6e3b757c8..e12681ce3 100644 --- a/src/clib/CMakeLists.txt +++ b/src/clib/CMakeLists.txt @@ -123,6 +123,7 @@ add_library(Grackle_Grackle rate_utils.cpp scale_fields.cpp scale_fields.hpp solve_rate_cool_g-cpp.cpp solve_rate_cool_g-cpp.h + step_rate_gauss_seidel.hpp step_rate_newton_raphson.hpp time_deriv_0d.hpp utils-cpp.cpp utils-cpp.hpp diff --git a/src/clib/calc_temp_cloudy_g.cpp b/src/clib/calc_temp_cloudy_g.cpp index 2a8ed9522..fc75dd17e 100644 --- a/src/clib/calc_temp_cloudy_g.cpp +++ b/src/clib/calc_temp_cloudy_g.cpp @@ -41,7 +41,16 @@ void calc_temp_cloudy_g(gr_float* temperature_data_, int imetal, const double dom = internalu_calc_dom_(internalu); const double zr = 1. / (internalu.a_value * internalu.a_units) - 1.; - // Convert densities from comoving to proper + // this assertion is a hint to clang-analyzer about the relationship between + // `imetal` and whether `metal_density` is a nullptr + // -> for context, `scale_fields_table` is inlined into this function. Rather + // than look at `imetal`, it checks if `metal_density` is a nullptr + // -> without this assertion clang-tidy will infer that since there's an + // explicit check for whether `metal_density` is null, regardless of + // `imetal`'s value, then it must be possible for `metal_density` to be + // null when `imetal` is 1. Then, it would report an error + GR_INTERNAL_REQUIRE((imetal != 1) || (my_fields->metal_density != nullptr), + "imetal has an incorrect value"); if (internalu.extfields_in_comoving == 1) { double factor = std::pow(internalu.a_value, -3); diff --git a/src/clib/chemistry_solver_funcs.hpp b/src/clib/chemistry_solver_funcs.hpp index b099352f7..f36db0a7f 100644 --- a/src/clib/chemistry_solver_funcs.hpp +++ b/src/clib/chemistry_solver_funcs.hpp @@ -19,11 +19,1421 @@ #define CHEMISTRY_SOLVER_FUNCS_HPP #include "grackle.h" +#include "fortran_func_decls.h" // gr_mask_type +#include "index_helper.h" #include "internal_types.hpp" #include "LUT.hpp" namespace grackle::impl::chemistry { +/// Perform an implicit Gauss-Seidel sweep of a backward-Euler time integrator +/// to advance the rate equations by one (sub-)cycle (for each index in the +/// index-range that is selected by the given itmask) +/// +/// @param[out] out_spdens Holds buffers that are used to store the updated +/// species densities for the @p idx_range. +/// @param[in] idx_range Specifies the current index-range +/// @param[in] dtit Specifies the timestep of the current sub-cycle for each +/// index in @p idx_range. +/// @param[in] anydust Indicates whether we are modelling dust +/// @param[in] h2dust Specifies the rate of H2 dust-formation on dust grains +/// for eacg k +/// @param[in] rhoH Indicates the mass density of all Hydrogen +/// @param[in] itmask The general iteration mask for @p idx_range. +/// @param[in] itmask_metal The iteration mask @p idx_range that specifies +/// where we should account for metal-species chemistry and grain-species +/// chemistry. +/// @param[in] my_chemistry Provides various runtime parameters (we probably +/// don't need to pass the whole thing) +/// @param[in] my_fields Specifies the current values of the field data +/// @param[in] my_uvb_rates specifies precomputed rxn rates dependent on the +/// UV background, without accounting for self-shield (we probably don't +/// need to pass the whole thing since we also pass kshield_buf) +/// @param[in] grain_growth_rates, kcr_buf, kshield_buf specifies the +/// precomputed rxn rates (depends on local physical conditions) +/// +/// Refactoring Goals +/// ----------------- +/// +/// Right now the implementation of every rate is **very** manual, which makes +/// it very easy to forget to add a rate. +/// +/// The shorter-term goal is to deduplicate as much as possible with between +/// this function & grackle::impl::chemistry::species_density_derivatives_0d. +/// Some important considerations: +/// - it will be easiest to deduplicate the metal-chemistry and grain-species +/// growth rates (we only need to introduce some light usage of templates) +/// - deduplicating other logic will involve additional templates +/// - introducing light-weight templates to deduplicate the metal-chemistry +/// and grain-species growth rate logic is well worth the cost to improve +/// the maintenance burden. Deduplicating the other logic is probably also +/// worth the cost (but it's more debatable) +/// +/// The longer-term goal is to overhaul this function. +/// - In the longer term, it would make more sense to use more of a "table +/// based approach" in the regime where we have lots of chemistry. +/// - We elaborate a little more down below: +/// - we already have a table of 1d rate buffers (i.e. @p kcol_buf ). +/// - we might also track a table (with matching indices) where we track the +/// indexes associated with each reactant/product of a rate as well as +/// stoichiometric coefficients. +/// - Furthermore, we could envision tracking rate-lists for each species of +/// all creational reactions and destructive reactions. +/// - If we have all of this, then we could just iterate over each species' +/// reactions and use the tables to dynamically access reactant densities +/// by index. Essentially, we would dynamically build up the creational +/// & destructive coefficients (``acoef`` & ``scoef``) for a full row and +/// use the values to then perform the update. +/// - in practice things would be a little more complex since we have lots +/// of rates not stored in @p kcol_buf, but it's still **very** doable +/// - once this infrastructure is in place, we could replace +/// grackle::impl::chemistry::species_density_derivatives_0d +/// that iterates through all of the rates in a very different manner +/// and directly compute partial derivatives in a much more efficient +/// manner (currently, it's used with finite derivatives for computing +/// these partial derivatives) +/// - Importantly, the table-based approach is *probably* slower than the +/// manual, hand-coded approach in the limit of a small chemical network. +/// I suspect that we'll preserve the hand-coded approach for +/// primordial_chemistry == 1, probably primordial_chemistry == 2, and +/// maybe even primordial_chemistry == 3 (but, we'll obviously need to test +/// performance). +inline void species_density_updates_gauss_seidel( + grackle::impl::SpeciesCollection out_spdens, IndexRange idx_range, + const double* dtit, gr_mask_type anydust, const double* h2dust, + const double* rhoH, const gr_mask_type* itmask, + const gr_mask_type* itmask_metal, chemistry_data* my_chemistry, + grackle_field_data* my_fields, photo_rate_storage my_uvb_rates, + grackle::impl::GrainSpeciesCollection grain_growth_rates, + grackle::impl::CollisionalRxnRateCollection kcol_buf, + grackle::impl::PhotoRxnRateCollection kshield_buf +) +{ + + // Construct views of various species fields + // ----------------------------------------- + grackle::impl::View de(my_fields->e_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View HI(my_fields->HI_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View HII(my_fields->HII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View HeI(my_fields->HeI_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View HeII(my_fields->HeII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View HeIII(my_fields->HeIII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View HM(my_fields->HM_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View H2I(my_fields->H2I_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View H2II(my_fields->H2II_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View DI(my_fields->DI_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View DII(my_fields->DII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View HDI(my_fields->HDI_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View DM(my_fields->DM_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View HDII(my_fields->HDII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View HeHII(my_fields->HeHII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View CI(my_fields->CI_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View CII(my_fields->CII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View CO(my_fields->CO_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View CO2(my_fields->CO2_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View OI(my_fields->OI_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View OH(my_fields->OH_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View H2O(my_fields->H2O_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View O2(my_fields->O2_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View SiI(my_fields->SiI_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View SiOI(my_fields->SiOI_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View SiO2I(my_fields->SiO2I_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View CH(my_fields->CH_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View CH2(my_fields->CH2_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View COII(my_fields->COII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View OII(my_fields->OII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View OHII(my_fields->OHII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View H2OII(my_fields->H2OII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View H3OII(my_fields->H3OII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View O2II(my_fields->O2II_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View Mg(my_fields->Mg_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View Al(my_fields->Al_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View S(my_fields->S_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View Fe(my_fields->Fe_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View SiM(my_fields->SiM_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View FeM(my_fields->FeM_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View Mg2SiO4(my_fields->Mg2SiO4_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View MgSiO3(my_fields->MgSiO3_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View Fe3O4(my_fields->Fe3O4_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View AC(my_fields->AC_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View SiO2D(my_fields->SiO2_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View MgO(my_fields->MgO_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View FeS(my_fields->FeS_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View Al2O3(my_fields->Al2O3_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View reforg(my_fields->ref_org_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View volorg(my_fields->vol_org_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View H2Oice(my_fields->H2O_ice_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + + // Radiation Fields + grackle::impl::View kphHI(my_fields->RT_HI_ionization_rate, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View kphHeI(my_fields->RT_HeI_ionization_rate, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View kphHeII(my_fields->RT_HeII_ionization_rate, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View kdissHDI(my_fields->RT_HDI_dissociation_rate, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View kphCI(my_fields->RT_CI_ionization_rate, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View kphOI(my_fields->RT_OI_ionization_rate, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View kdissCO(my_fields->RT_CO_dissociation_rate, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View kdissOH(my_fields->RT_OH_dissociation_rate, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View kdissH2O(my_fields->RT_H2O_dissociation_rate, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + + // locals + + int i; + double scoef, acoef; + + const int j = idx_range.j; + const int k = idx_range.k; + + // A) the 6-species integrator + if (my_chemistry->primordial_chemistry == 1) { + + for (i = idx_range.i_start; i < idx_range.i_stop; i++) { + if (itmask[i] != MASK_FALSE) { + + // 1) HI + + scoef = kcol_buf.data[CollisionalRxnLUT::k2][i]*HII(i,j,k)*de(i,j,k); + acoef = kcol_buf.data[CollisionalRxnLUT::k1][i]*de(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k57][i]*HI(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k58][i]*HeI(i,j,k)/4. + + kshield_buf.k24[i]; + if (my_chemistry->use_radiative_transfer == 1) { acoef = acoef + kphHI(i,j,k); } + out_spdens.data[SpLUT::HI][i] = (scoef*dtit[i] + HI(i,j,k))/ + (1. + acoef*dtit[i]); + if (out_spdens.data[SpLUT::HI][i] != out_spdens.data[SpLUT::HI][i]) { + OMP_PRAGMA_CRITICAL + { + std::printf("HUGE HIp! :: %d %d %d %g %g %g %g %g %g %g %g\n", + i, + j, + k, + out_spdens.data[SpLUT::HI] [ i ], + HI ( i, j, k ), + HII ( i, j, k ), + de ( i, j, k ), + kphHI ( i, j, k ), + scoef, + acoef, + dtit [ i ]); + } + // ERROR_MESSAGE + } + + // 2) HII + scoef = kcol_buf.data[CollisionalRxnLUT::k1][i]*out_spdens.data[SpLUT::HI][i]*de(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k57][i]*out_spdens.data[SpLUT::HI][i]*out_spdens.data[SpLUT::HI][i] + + kcol_buf.data[CollisionalRxnLUT::k58][i]*out_spdens.data[SpLUT::HI][i]*HeI(i,j,k)/4. + + kshield_buf.k24[i]*out_spdens.data[SpLUT::HI][i]; + if (my_chemistry->use_radiative_transfer == 1) + { scoef = scoef + kphHI(i,j,k)*out_spdens.data[SpLUT::HI][i]; } + acoef = kcol_buf.data[CollisionalRxnLUT::k2][i]*de (i,j,k); + out_spdens.data[SpLUT::HII][i] = (scoef*dtit[i] + HII(i,j,k))/ + (1. +acoef*dtit[i]); + // + if (out_spdens.data[SpLUT::HII][i] <= 0.) //##### + { + OMP_PRAGMA_CRITICAL + { + std::printf("negative HIIp! :: %d %d %d %g %g %g %g %g %g %g %g %g %g\n", + i, + j, + k, + out_spdens.data[SpLUT::HII] [ i ], + scoef, + dtit [ i ], + HII ( i, j, k ), + acoef, + kcol_buf.data[CollisionalRxnLUT::k2] [ i ], + de ( i, j, k ), + kphHI ( i, j, k ), + out_spdens.data[SpLUT::HI] [ i ], + kshield_buf.k24 [ i ]); + } + } + + // 3) Electron density + + scoef = 0. + + kcol_buf.data[CollisionalRxnLUT::k57][i]*out_spdens.data[SpLUT::HI][i]*out_spdens.data[SpLUT::HI][i] + + kcol_buf.data[CollisionalRxnLUT::k58][i]*out_spdens.data[SpLUT::HI][i]*HeI(i,j,k)/4. + + kshield_buf.k24[i]*HI(i,j,k) + + kshield_buf.k25[i]*HeII(i,j,k)/4. + + kshield_buf.k26[i]*HeI(i,j,k)/4.; + + if ( (my_chemistry->use_radiative_transfer == 1) && ( my_chemistry->radiative_transfer_hydrogen_only == 0) ) + { scoef = scoef + kphHI(i,j,k) * HI(i,j,k) + + kphHeI(i,j,k) * HeI(i,j,k) / 4. + + kphHeII(i,j,k) * HeII(i,j,k) / 4.; } + if ( (my_chemistry->use_radiative_transfer == 1) && ( my_chemistry->radiative_transfer_hydrogen_only == 1) ) + { scoef = scoef + kphHI(i,j,k) * HI(i,j,k); } + + + + acoef = -(kcol_buf.data[CollisionalRxnLUT::k1][i]*HI(i,j,k) - kcol_buf.data[CollisionalRxnLUT::k2][i]*HII(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k3][i]*HeI(i,j,k)/4. - + kcol_buf.data[CollisionalRxnLUT::k6][i]*HeIII(i,j,k)/4. + + kcol_buf.data[CollisionalRxnLUT::k5][i]*HeII(i,j,k)/4. - + kcol_buf.data[CollisionalRxnLUT::k4][i]*HeII(i,j,k)/4.); + out_spdens.data[SpLUT::e][i] = (scoef*dtit[i] + de(i,j,k)) + / (1. + acoef*dtit[i]); + + } + } + + } + + // --- (B) Do helium chemistry in any case: (for all ispecies values) --- + + for (i = idx_range.i_start; i < idx_range.i_stop; i++) { + if (itmask[i] != MASK_FALSE) { + + // 4) HeI + + scoef = kcol_buf.data[CollisionalRxnLUT::k4][i]*HeII(i,j,k)*de(i,j,k); + acoef = kcol_buf.data[CollisionalRxnLUT::k3][i]*de(i,j,k) + + kshield_buf.k26[i]; + + if ( (my_chemistry->use_radiative_transfer == 1) && (my_chemistry->radiative_transfer_hydrogen_only == 0)) + { acoef = acoef + kphHeI(i,j,k); } + if (my_chemistry->primordial_chemistry > 3) { + scoef = scoef + 4. * ( 0. + + kcol_buf.data[CollisionalRxnLUT::k152][i] * HeHII(i,j,k) * HI(i,j,k) / 5. + + kcol_buf.data[CollisionalRxnLUT::k153][i] * HeHII(i,j,k) * de(i,j,k) / 5. + ); + acoef = acoef + + kcol_buf.data[CollisionalRxnLUT::k148][i] * HII(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k149][i] * HII(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k150][i] * H2II(i,j,k) / 2.; + } + out_spdens.data[SpLUT::HeI][i] = ( scoef*dtit[i] + HeI(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + // 5) HeII + + scoef = kcol_buf.data[CollisionalRxnLUT::k3][i]*out_spdens.data[SpLUT::HeI][i]*de(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k6][i]*HeIII(i,j,k)*de(i,j,k) + + kshield_buf.k26[i]*out_spdens.data[SpLUT::HeI][i]; + + if ( (my_chemistry->use_radiative_transfer == 1) && (my_chemistry->radiative_transfer_hydrogen_only == 0)) + { scoef = scoef + kphHeI(i,j,k)*out_spdens.data[SpLUT::HeI][i]; } + + acoef = kcol_buf.data[CollisionalRxnLUT::k4][i]*de(i,j,k) + kcol_buf.data[CollisionalRxnLUT::k5][i]*de(i,j,k) + + kshield_buf.k25[i]; + + if ( (my_chemistry->use_radiative_transfer == 1) && (my_chemistry->radiative_transfer_hydrogen_only == 0)) + { acoef = acoef + kphHeII(i,j,k); } + if (my_chemistry->primordial_chemistry > 3) { + acoef = acoef + + kcol_buf.data[CollisionalRxnLUT::k151][i] * HI(i,j,k); + } + out_spdens.data[SpLUT::HeII][i] = ( scoef*dtit[i] + HeII(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + // 6) HeIII + + scoef = kcol_buf.data[CollisionalRxnLUT::k5][i]*out_spdens.data[SpLUT::HeII][i]*de(i,j,k) + + kshield_buf.k25[i]*out_spdens.data[SpLUT::HeII][i]; + if ((my_chemistry->use_radiative_transfer == 1) && (my_chemistry->radiative_transfer_hydrogen_only == 0)) + { scoef = scoef + kphHeII(i,j,k) * out_spdens.data[SpLUT::HeII][i]; } + acoef = kcol_buf.data[CollisionalRxnLUT::k6][i]*de(i,j,k); + out_spdens.data[SpLUT::HeIII][i] = ( scoef*dtit[i] + HeIII(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + } + } + + // --- (C) Now do extra 3-species for molecular hydrogen --- + + if (my_chemistry->primordial_chemistry > 1) { + + // First, do HI/HII with molecular hydrogen terms + + for (i = idx_range.i_start; i < idx_range.i_stop; i++) { + if (itmask[i] != MASK_FALSE) { + + // 1) HI + scoef = kcol_buf.data[CollisionalRxnLUT::k2][i] * HII(i,j,k) * de(i,j,k) + + 2.*kcol_buf.data[CollisionalRxnLUT::k13][i]* HI(i,j,k) * H2I(i,j,k)/2. + + kcol_buf.data[CollisionalRxnLUT::k11][i]* HII(i,j,k) * H2I(i,j,k)/2. + + 2.*kcol_buf.data[CollisionalRxnLUT::k12][i]* de(i,j,k) * H2I(i,j,k)/2. + + kcol_buf.data[CollisionalRxnLUT::k14][i]* HM(i,j,k) * de(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k15][i]* HM(i,j,k) * HI(i,j,k) + + 2.*kcol_buf.data[CollisionalRxnLUT::k16][i]* HM(i,j,k) * HII(i,j,k) + + 2.*kcol_buf.data[CollisionalRxnLUT::k18][i]* H2II(i,j,k)* de(i,j,k)/2. + + kcol_buf.data[CollisionalRxnLUT::k19][i]* H2II(i,j,k)* HM(i,j,k)/2. + + 2.*kshield_buf.k31[i] * H2I(i,j,k)/2.; + + acoef = kcol_buf.data[CollisionalRxnLUT::k1][i] * de(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k7][i] * de(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k8][i] * HM(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k9][i] * HII(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k10][i]* H2II(i,j,k)/2. + + 2.*kcol_buf.data[CollisionalRxnLUT::k22][i]* std::pow(HI(i,j,k),2) + + kcol_buf.data[CollisionalRxnLUT::k57][i]* HI(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k58][i]* HeI(i,j,k)/4. + + kshield_buf.k24[i]; + + if (my_chemistry->use_radiative_transfer == 1) { acoef = acoef + kphHI(i,j,k); } + if (my_chemistry->use_radiative_transfer == 1) { + if ((my_chemistry->primordial_chemistry > 2) && (my_chemistry->radiative_transfer_HDI_dissociation > 0)) { + scoef = scoef + + kdissHDI(i,j,k) * HDI(i,j,k)/3.0; + } + if ( (my_chemistry->metal_chemistry == 1) && + (itmask_metal[i] != MASK_FALSE) ) { + if (my_chemistry->radiative_transfer_metal_dissociation > 0) { + scoef = scoef + + kdissOH (i,j,k) * OH(i,j,k) /17.0 + + kdissH2O(i,j,k) * H2O(i,j,k)/18.0; + } + } + } + + if (anydust != MASK_FALSE) { + if(itmask_metal[i] != MASK_FALSE) { + acoef = acoef + 2. * h2dust[i] * rhoH[i]; + } + } +#ifdef CONTRIBUTION_OF_MINOR_SPECIES + if (my_chemistry->primordial_chemistry > 2) { + scoef = scoef + + kcol_buf.data[CollisionalRxnLUT::k50][i] * HII(i,j,k) * DI(i,j,k) / 2. + + kcol_buf.data[CollisionalRxnLUT::k54][i] * H2I(i,j,k) * DI(i,j,k) / 4.; + acoef = acoef + + kcol_buf.data[CollisionalRxnLUT::k51][i] * DII(i,j,k) / 2. + + kcol_buf.data[CollisionalRxnLUT::k55][i] * HDI(i,j,k) / 3.; + } +#endif + if (my_chemistry->primordial_chemistry > 3) { + scoef = scoef + + kcol_buf.data[CollisionalRxnLUT::k131][i] * HDII(i,j,k) * de(i,j,k) / 3. + + kcol_buf.data[CollisionalRxnLUT::k134][i] * HII(i,j,k) * DM(i,j,k) / 2. + + kcol_buf.data[CollisionalRxnLUT::k135][i] * HM(i,j,k) * DI(i,j,k) / 2. + + kcol_buf.data[CollisionalRxnLUT::k150][i] * HeI(i,j,k) * H2II(i,j,k) / 8. + + kcol_buf.data[CollisionalRxnLUT::k153][i] * HeHII(i,j,k) * de(i,j,k) / 5.; + acoef = acoef + + kcol_buf.data[CollisionalRxnLUT::k125][i] * HDII(i,j,k) / 3. + + kcol_buf.data[CollisionalRxnLUT::k130][i] * DII(i,j,k) / 2. + + kcol_buf.data[CollisionalRxnLUT::k136][i] * DM(i,j,k) / 2. + + kcol_buf.data[CollisionalRxnLUT::k137][i] * DM(i,j,k) / 2. + + kcol_buf.data[CollisionalRxnLUT::k151][i] * HeII(i,j,k) / 4. + + kcol_buf.data[CollisionalRxnLUT::k152][i] * HeHII(i,j,k) / 5.; + } + + if ( (my_chemistry->metal_chemistry == 1) && + (itmask_metal[i] != MASK_FALSE) ) { + scoef = scoef + + kcol_buf.data[CollisionalRxnLUT::kz20][i] * CI(i,j,k) * H2I(i,j,k) / 24. + + kcol_buf.data[CollisionalRxnLUT::kz21][i] * OI(i,j,k) * H2I(i,j,k) / 32. + + kcol_buf.data[CollisionalRxnLUT::kz22][i] * HII(i,j,k) * OI(i,j,k) / 16. + + kcol_buf.data[CollisionalRxnLUT::kz23][i] * H2I(i,j,k) * CH(i,j,k) / 26. + + kcol_buf.data[CollisionalRxnLUT::kz24][i] * H2I(i,j,k) * OH(i,j,k) / 34. + + kcol_buf.data[CollisionalRxnLUT::kz26][i] * OH(i,j,k) * CO(i,j,k) / 476. + + kcol_buf.data[CollisionalRxnLUT::kz28][i] * CI(i,j,k) * OH(i,j,k) / 204. + + kcol_buf.data[CollisionalRxnLUT::kz32][i] * OI(i,j,k) * CH(i,j,k) / 208. + + kcol_buf.data[CollisionalRxnLUT::kz33][i] * OI(i,j,k) * OH(i,j,k) / 272. + + kcol_buf.data[CollisionalRxnLUT::kz34][i] * HII(i,j,k) * OH(i,j,k) / 17. + + kcol_buf.data[CollisionalRxnLUT::kz35][i] * HII(i,j,k) * H2O(i,j,k) / 18. + + kcol_buf.data[CollisionalRxnLUT::kz36][i] * HII(i,j,k) * O2(i,j,k) / 32. + + kcol_buf.data[CollisionalRxnLUT::kz37][i] * CII(i,j,k) * OH(i,j,k) / 204. + + kcol_buf.data[CollisionalRxnLUT::kz40][i] * OII(i,j,k) * H2I(i,j,k) / 32. + + kcol_buf.data[CollisionalRxnLUT::kz41][i] * OHII(i,j,k) * H2I(i,j,k) / 34. + + kcol_buf.data[CollisionalRxnLUT::kz42][i] * H2OII(i,j,k) * H2I(i,j,k) / 36. + + kcol_buf.data[CollisionalRxnLUT::kz46][i] * H2OII(i,j,k) * de(i,j,k) / 18. + + kcol_buf.data[CollisionalRxnLUT::kz48][i] * H3OII(i,j,k) * de(i,j,k) / 19. + + kcol_buf.data[CollisionalRxnLUT::kz49][i] * H3OII(i,j,k) * de(i,j,k) / 9.5 + + kcol_buf.data[CollisionalRxnLUT::kz52][i] * SiI(i,j,k) * OH(i,j,k) / 476. + + kcol_buf.data[CollisionalRxnLUT::kz54][i] * SiOI(i,j,k) * OH(i,j,k) / 748.; + acoef = acoef + + kcol_buf.data[CollisionalRxnLUT::kz15][i] * CH(i,j,k) / 13. + + kcol_buf.data[CollisionalRxnLUT::kz16][i] * CH2(i,j,k) / 14. + + kcol_buf.data[CollisionalRxnLUT::kz17][i] * OH(i,j,k) / 17. + + kcol_buf.data[CollisionalRxnLUT::kz18][i] * H2O(i,j,k) / 18. + + kcol_buf.data[CollisionalRxnLUT::kz19][i] * O2(i,j,k) / 32. + + kcol_buf.data[CollisionalRxnLUT::kz27][i] * CI(i,j,k) / 12. + + kcol_buf.data[CollisionalRxnLUT::kz30][i] * OI(i,j,k) / 16. + + kcol_buf.data[CollisionalRxnLUT::kz39][i] * OII(i,j,k) / 16. + + kcol_buf.data[CollisionalRxnLUT::kz43][i] * COII(i,j,k) / 28.; + } + out_spdens.data[SpLUT::HI][i] = ( scoef*dtit[i] + HI(i,j,k) ) / + ( 1.f + acoef*dtit[i] ); + if (out_spdens.data[SpLUT::HI][i] != out_spdens.data[SpLUT::HI][i]) { + OMP_PRAGMA_CRITICAL + { + std::printf("HUGE HIp! :: %d %d %d %g %g %g %g %g %g\n", + i, + j, + k, + out_spdens.data[SpLUT::HI] [ i ], + HI ( i, j, k ), + HII ( i, j, k ), + de ( i, j, k ), + H2I ( i, j, k ), + kphHI ( i, j, k )); + } + } + + // 2) HII + + scoef = kcol_buf.data[CollisionalRxnLUT::k1][i] * HI(i,j,k) * de(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k10][i] * H2II(i,j,k)*HI(i,j,k)/2. + + kcol_buf.data[CollisionalRxnLUT::k57][i] * HI(i,j,k) * HI(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k58][i] * HI(i,j,k) * HeI(i,j,k)/4. + + kshield_buf.k24[i]*HI(i,j,k); + + if (my_chemistry->use_radiative_transfer == 1) + { scoef = scoef + kphHI(i,j,k) * HI(i,j,k); } + + acoef = kcol_buf.data[CollisionalRxnLUT::k2][i] * de(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k9][i] * HI(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k11][i] * H2I(i,j,k)/2. + + kcol_buf.data[CollisionalRxnLUT::k16][i] * HM(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k17][i] * HM(i,j,k); +#ifdef CONTRIBUTION_OF_MINOR_SPECIES + if (my_chemistry->primordial_chemistry > 2) { + scoef = scoef + + kcol_buf.data[CollisionalRxnLUT::k51][i] * HI (i,j,k) * DII(i,j,k) / 2. + + kcol_buf.data[CollisionalRxnLUT::k52][i] * H2I(i,j,k) * DII(i,j,k) / 4.; + acoef = acoef + + kcol_buf.data[CollisionalRxnLUT::k50][i] * DI (i,j,k) / 2. + + kcol_buf.data[CollisionalRxnLUT::k53][i] * HDI(i,j,k) / 3.; + } +#endif + if (my_chemistry->primordial_chemistry > 3) { + scoef = scoef + + kcol_buf.data[CollisionalRxnLUT::k125][i] * HDII(i,j,k) * HI(i,j,k) / 3.; + acoef = acoef + + kcol_buf.data[CollisionalRxnLUT::k129][i] * DI(i,j,k) / 2. + + kcol_buf.data[CollisionalRxnLUT::k134][i] * DM(i,j,k) / 2. + + kcol_buf.data[CollisionalRxnLUT::k148][i] * HeI(i,j,k) / 4. + + kcol_buf.data[CollisionalRxnLUT::k149][i] * HeI(i,j,k) / 4.; + } + + if ( (my_chemistry->metal_chemistry == 1) && + (itmask_metal[i] != MASK_FALSE) ) { + scoef = scoef + + kcol_buf.data[CollisionalRxnLUT::kz39][i] * OII(i,j,k) * HI(i,j,k) / 16. + + kcol_buf.data[CollisionalRxnLUT::kz43][i] * COII(i,j,k) * HI(i,j,k) / 28.; + acoef = acoef + + kcol_buf.data[CollisionalRxnLUT::kz22][i] * OI(i,j,k) / 16. + + kcol_buf.data[CollisionalRxnLUT::kz34][i] * OH(i,j,k) / 17. + + kcol_buf.data[CollisionalRxnLUT::kz35][i] * H2O(i,j,k) / 18. + + kcol_buf.data[CollisionalRxnLUT::kz36][i] * O2(i,j,k) / 32.; + } + out_spdens.data[SpLUT::HII][i] = ( scoef*dtit[i] + HII(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + // 3) electrons: + + scoef = kcol_buf.data[CollisionalRxnLUT::k8][i] * HM(i,j,k) * HI(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k15][i]* HM(i,j,k) * HI(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k17][i]* HM(i,j,k) * HII(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k57][i]* HI(i,j,k) * HI(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k58][i]* HI(i,j,k) * HeI(i,j,k)/4. + // + + kshield_buf.k24[i]*out_spdens.data[SpLUT::HI][i] + + kshield_buf.k25[i]*out_spdens.data[SpLUT::HeII][i]/4. + + kshield_buf.k26[i]*out_spdens.data[SpLUT::HeI][i]/4.; + + if ( (my_chemistry->use_radiative_transfer == 1) && (my_chemistry->radiative_transfer_hydrogen_only == 0) ) + { scoef = scoef + kphHI(i,j,k) * out_spdens.data[SpLUT::HI][i] + + kphHeI(i,j,k) * out_spdens.data[SpLUT::HeI][i] / 4. + + kphHeII(i,j,k) * out_spdens.data[SpLUT::HeII][i] / 4.; } + if ( (my_chemistry->use_radiative_transfer == 1) && (my_chemistry->radiative_transfer_hydrogen_only == 1) ) + { scoef = scoef + kphHI(i,j,k) * out_spdens.data[SpLUT::HI][i]; } + if (my_chemistry->use_radiative_transfer == 1) { + if ( (my_chemistry->metal_chemistry == 1) && + (itmask_metal[i] != MASK_FALSE) ) { + if (my_chemistry->radiative_transfer_metal_ionization > 0) { + scoef = scoef + + kphCI(i,j,k) * CI(i,j,k)/12.0 + + kphOI(i,j,k) * OI(i,j,k)/16.0; + } + } + } + + acoef = - (kcol_buf.data[CollisionalRxnLUT::k1][i] *HI(i,j,k) - kcol_buf.data[CollisionalRxnLUT::k2][i]*HII(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k3][i] *HeI(i,j,k)/4. - + kcol_buf.data[CollisionalRxnLUT::k6][i]*HeIII(i,j,k)/4. + + kcol_buf.data[CollisionalRxnLUT::k5][i] *HeII(i,j,k)/4. - + kcol_buf.data[CollisionalRxnLUT::k4][i]*HeII(i,j,k)/4. + + kcol_buf.data[CollisionalRxnLUT::k14][i]*HM(i,j,k) + - kcol_buf.data[CollisionalRxnLUT::k7][i] *HI(i,j,k) + - kcol_buf.data[CollisionalRxnLUT::k18][i]*H2II(i,j,k)/2.); +#ifdef CONTRIBUTION_OF_MINOR_SPECIES + if (my_chemistry->primordial_chemistry > 2) { + scoef = scoef + + kcol_buf.data[CollisionalRxnLUT::k56][i] * DI (i,j,k) * HM(i,j,k) / 2.; + acoef = acoef + - kcol_buf.data[CollisionalRxnLUT::k1] [i] * DI (i,j,k) / 2. + + kcol_buf.data[CollisionalRxnLUT::k2] [i] * DII(i,j,k) / 2.; + } +#endif + if (my_chemistry->primordial_chemistry > 3) { + scoef = scoef + + kcol_buf.data[CollisionalRxnLUT::k137][i] * DM(i,j,k) * HI(i,j,k) / 2.; + acoef = acoef + + kcol_buf.data[CollisionalRxnLUT::k131][i] * HDII(i,j,k) / 3. + + kcol_buf.data[CollisionalRxnLUT::k132][i] * DI(i,j,k) / 2. + + kcol_buf.data[CollisionalRxnLUT::k153][i] * HeHII(i,j,k) / 5.; + } + + if ( (my_chemistry->metal_chemistry == 1) && + (itmask_metal[i] != MASK_FALSE) ) { + scoef = scoef; + acoef = acoef + + kcol_buf.data[CollisionalRxnLUT::kz44][i] * CII(i,j,k) / 12. + + kcol_buf.data[CollisionalRxnLUT::kz45][i] * OII(i,j,k) / 16. + + kcol_buf.data[CollisionalRxnLUT::kz46][i] * H2OII(i,j,k) / 18. + + kcol_buf.data[CollisionalRxnLUT::kz47][i] * H2OII(i,j,k) / 18. + + kcol_buf.data[CollisionalRxnLUT::kz48][i] * H3OII(i,j,k) / 19. + + kcol_buf.data[CollisionalRxnLUT::kz49][i] * H3OII(i,j,k) / 19. + + kcol_buf.data[CollisionalRxnLUT::kz50][i] * O2II(i,j,k) / 32.; + } + out_spdens.data[SpLUT::e][i] = ( scoef*dtit[i] + de(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + // 7) H2 + + scoef = 2.*(kcol_buf.data[CollisionalRxnLUT::k8][i] * HM(i,j,k) * HI(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k10][i] * H2II(i,j,k) * HI(i,j,k)/2. + + kcol_buf.data[CollisionalRxnLUT::k19][i] * H2II(i,j,k) * HM(i,j,k)/2. + + kcol_buf.data[CollisionalRxnLUT::k22][i] * HI(i,j,k) * std::pow((HI(i,j,k)),2.)); + acoef = ( kcol_buf.data[CollisionalRxnLUT::k13][i]*HI(i,j,k) + kcol_buf.data[CollisionalRxnLUT::k11][i]*HII(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k12][i]*de(i,j,k) ) + + kshield_buf.k29[i] + kshield_buf.k31[i]; + + if (anydust != MASK_FALSE) { + if(itmask_metal[i] != MASK_FALSE) { + scoef = scoef + 2. * h2dust[i] * + HI(i,j,k) * rhoH[i]; + } + } +#ifdef CONTRIBUTION_OF_MINOR_SPECIES + if (my_chemistry->primordial_chemistry > 2) { + scoef = scoef + 2. * ( + kcol_buf.data[CollisionalRxnLUT::k53][i] * HDI(i,j,k) * HII(i,j,k) / 3. + + kcol_buf.data[CollisionalRxnLUT::k55][i] * HDI(i,j,k) * HI (i,j,k) / 3. + ); + acoef = acoef + + kcol_buf.data[CollisionalRxnLUT::k52][i] * DII(i,j,k) / 2. + + kcol_buf.data[CollisionalRxnLUT::k54][i] * DI (i,j,k) / 2.; + } +#endif + if ( (my_chemistry->metal_chemistry == 1) && + (itmask_metal[i] != MASK_FALSE) ) { + scoef = scoef + 2. * ( 0. + + kcol_buf.data[CollisionalRxnLUT::kz15][i] * HI(i,j,k) * CH(i,j,k) / 13. + + kcol_buf.data[CollisionalRxnLUT::kz16][i] * HI(i,j,k) * CH2(i,j,k) / 14. + + kcol_buf.data[CollisionalRxnLUT::kz17][i] * HI(i,j,k) * OH(i,j,k) / 17. + + kcol_buf.data[CollisionalRxnLUT::kz18][i] * HI(i,j,k) * H2O(i,j,k) / 18. + + kcol_buf.data[CollisionalRxnLUT::kz47][i] * H2OII(i,j,k) * de(i,j,k) / 18. + ); + acoef = acoef + + kcol_buf.data[CollisionalRxnLUT::kz20][i] * CI(i,j,k) / 12. + + kcol_buf.data[CollisionalRxnLUT::kz21][i] * OI(i,j,k) / 16. + + kcol_buf.data[CollisionalRxnLUT::kz23][i] * CH(i,j,k) / 13. + + kcol_buf.data[CollisionalRxnLUT::kz24][i] * OH(i,j,k) / 17. + + kcol_buf.data[CollisionalRxnLUT::kz40][i] * OII(i,j,k) / 16. + + kcol_buf.data[CollisionalRxnLUT::kz41][i] * OHII(i,j,k) / 17. + + kcol_buf.data[CollisionalRxnLUT::kz42][i] * H2OII(i,j,k) / 18. + + kcol_buf.data[CollisionalRxnLUT::kz51][i] * CI(i,j,k) / 12.; + if ((my_chemistry->grain_growth == 1) || (my_chemistry->dust_sublimation == 1)) { + if (my_chemistry->dust_species > 0) { + scoef = scoef + 2. * + grain_growth_rates.data[OnlyGrainSpLUT::MgSiO3_dust] [i] * 2.; + + } + if (my_chemistry->dust_species > 1) { + scoef = scoef + 2. * ( + grain_growth_rates.data[OnlyGrainSpLUT::Mg2SiO4_dust] [i] * 3. + + grain_growth_rates.data[OnlyGrainSpLUT::Fe3O4_dust] [i] * 4. + + grain_growth_rates.data[OnlyGrainSpLUT::MgO_dust] [i] + + grain_growth_rates.data[OnlyGrainSpLUT::Al2O3_dust] [i] * 3. + ); + } + if (my_chemistry->dust_species > 2) { + acoef = acoef + + grain_growth_rates.data[OnlyGrainSpLUT::vol_org_dust] [i] / H2I(i,j,k) * 2. * 2.; + } + } + } + out_spdens.data[SpLUT::H2I][i] = ( scoef*dtit[i] + H2I(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + // 8) H- + + scoef = kcol_buf.data[CollisionalRxnLUT::k7][i] * HI(i,j,k) * de(i,j,k); + acoef = (kcol_buf.data[CollisionalRxnLUT::k8][i] + kcol_buf.data[CollisionalRxnLUT::k15][i]) * HI(i,j,k) + + (kcol_buf.data[CollisionalRxnLUT::k16][i] + kcol_buf.data[CollisionalRxnLUT::k17][i]) * HII(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k14][i] * de(i,j,k) + kcol_buf.data[CollisionalRxnLUT::k19][i] * H2II(i,j,k)/2.0f + + my_uvb_rates.k27; +#ifdef CONTRIBUTION_OF_MINOR_SPECIES + if (my_chemistry->primordial_chemistry > 2) { + acoef = acoef + + kcol_buf.data[CollisionalRxnLUT::k56][i] * DI (i,j,k) / 2.; + } +#endif + if (my_chemistry->primordial_chemistry > 3) { + scoef = scoef + + kcol_buf.data[CollisionalRxnLUT::k136][i] * DM(i,j,k) * HI(i,j,k) / 2.; + acoef = acoef + + kcol_buf.data[CollisionalRxnLUT::k135][i] * DI(i,j,k) / 2.; + } + out_spdens.data[SpLUT::HM][i] = (scoef*dtit[i] + HM(i,j,k)) + / (1.0f + acoef*dtit[i]); + + + // 9) H2+ + + out_spdens.data[SpLUT::H2II][i] = 2.*( kcol_buf.data[CollisionalRxnLUT::k9] [i]*out_spdens.data[SpLUT::HI][i]*out_spdens.data[SpLUT::HII][i] + + kcol_buf.data[CollisionalRxnLUT::k11][i]*out_spdens.data[SpLUT::H2I][i]/2.*out_spdens.data[SpLUT::HII][i] + + kcol_buf.data[CollisionalRxnLUT::k17][i]*out_spdens.data[SpLUT::HM][i]*out_spdens.data[SpLUT::HII][i] + + kshield_buf.k29[i]*out_spdens.data[SpLUT::H2I][i] + ) + / ( kcol_buf.data[CollisionalRxnLUT::k10][i]*out_spdens.data[SpLUT::HI][i] + kcol_buf.data[CollisionalRxnLUT::k18][i]*out_spdens.data[SpLUT::e][i] + + kcol_buf.data[CollisionalRxnLUT::k19][i]*out_spdens.data[SpLUT::HM][i] + + (kshield_buf.k28[i]+kshield_buf.k30[i]) + ); + if (my_chemistry->primordial_chemistry > 3) { + out_spdens.data[SpLUT::H2II][i] = 2. * ( kcol_buf.data[CollisionalRxnLUT::k9] [i]*out_spdens.data[SpLUT::HI][i]*out_spdens.data[SpLUT::HII][i] + + kcol_buf.data[CollisionalRxnLUT::k11][i]*out_spdens.data[SpLUT::H2I][i]/2.*out_spdens.data[SpLUT::HII][i] + + kcol_buf.data[CollisionalRxnLUT::k17][i]*out_spdens.data[SpLUT::HM][i]*out_spdens.data[SpLUT::HII][i] + + kshield_buf.k29[i]*out_spdens.data[SpLUT::H2I][i] + + kcol_buf.data[CollisionalRxnLUT::k152][i]*HeHII(i,j,k)*out_spdens.data[SpLUT::HI][i]/5. + ) + / ( kcol_buf.data[CollisionalRxnLUT::k10][i]*out_spdens.data[SpLUT::HI][i] + kcol_buf.data[CollisionalRxnLUT::k18][i]*out_spdens.data[SpLUT::e][i] + + kcol_buf.data[CollisionalRxnLUT::k19][i]*out_spdens.data[SpLUT::HM][i] + + (kshield_buf.k28[i]+kshield_buf.k30[i]) + + kcol_buf.data[CollisionalRxnLUT::k150][i]*out_spdens.data[SpLUT::HeI][i]/4. + ); + } + } + } + // + } + + // --- (D) Now do extra 3-species for molecular HD --- + if (my_chemistry->primordial_chemistry > 2) { + for (i = idx_range.i_start; i < idx_range.i_stop; i++) { + if (itmask[i] != MASK_FALSE) { + + // 1) DI + scoef = ( kcol_buf.data[CollisionalRxnLUT::k2][i] * DII(i,j,k) * de(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k51][i]* DII(i,j,k) * HI(i,j,k) + + 2.*kcol_buf.data[CollisionalRxnLUT::k55][i]* HDI(i,j,k) * + HI(i,j,k)/3. + ); + acoef = kcol_buf.data[CollisionalRxnLUT::k1][i] * de(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k50][i] * HII(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k54][i] * H2I(i,j,k)/2. + + kcol_buf.data[CollisionalRxnLUT::k56][i] * HM(i,j,k) + + kshield_buf.k24[i]; + if (my_chemistry->use_radiative_transfer == 1) { acoef = acoef + kphHI(i,j,k); } + if (my_chemistry->primordial_chemistry > 3) { + scoef = scoef + 2. * ( 0. + + kcol_buf.data[CollisionalRxnLUT::k131][i] * HDII(i,j,k) * de(i,j,k) / 3. + + kcol_buf.data[CollisionalRxnLUT::k133][i] * DII(i,j,k) * DM(i,j,k) / 2. + + kcol_buf.data[CollisionalRxnLUT::k134][i] * HII(i,j,k) * DM(i,j,k) / 2. + + kcol_buf.data[CollisionalRxnLUT::k136][i] * DM(i,j,k) * HI(i,j,k) / 2. + ); + acoef = acoef + + kcol_buf.data[CollisionalRxnLUT::k129][i] * HII(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k132][i] * de(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k135][i] * HM(i,j,k); + } + if (my_chemistry->use_radiative_transfer == 1) { + if (my_chemistry->radiative_transfer_HDI_dissociation > 0) { + scoef = scoef + + 2. * kdissHDI(i,j,k) * HDI(i,j,k)/3.0; + } + } + out_spdens.data[SpLUT::DI][i] = ( scoef*dtit[i] + DI(i,j,k) ) / + ( 1. + acoef*dtit[i] ); + + // 2) DII + scoef = ( kcol_buf.data[CollisionalRxnLUT::k1][i] * DI(i,j,k) * de(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k50][i] * HII(i,j,k)* DI(i,j,k) + + 2.*kcol_buf.data[CollisionalRxnLUT::k53][i] * HII(i,j,k)* HDI(i,j,k)/3. + ) + + kshield_buf.k24[i]*DI(i,j,k); + acoef = 0.; + // ! initialize GC202002 + if (my_chemistry->use_radiative_transfer == 1) { scoef = scoef + kphHI(i,j,k)*DI(i,j,k); } + acoef = kcol_buf.data[CollisionalRxnLUT::k2][i] * de(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k51][i] * HI(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k52][i] * H2I(i,j,k)/2.; + if (my_chemistry->primordial_chemistry > 3) { + acoef = acoef + + kcol_buf.data[CollisionalRxnLUT::k130][i] * HI(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k133][i] * DM(i,j,k) / 2.; + } + out_spdens.data[SpLUT::DII][i] = ( scoef*dtit[i] + DII(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + // 3) HDI + scoef = 3.*(kcol_buf.data[CollisionalRxnLUT::k52][i] * DII(i,j,k)* + H2I(i,j,k)/2./2. + + kcol_buf.data[CollisionalRxnLUT::k54][i] * DI(i,j,k) * H2I(i,j,k)/2./2. + // ! & + 2._DKIND*k56(i) * DI(i,j,k) * HM(i,j,k)/2._DKIND + //- ! corrected by GC202005 + + kcol_buf.data[CollisionalRxnLUT::k56][i] * DI(i,j,k) * HM(i,j,k)/2. + ); + acoef = kcol_buf.data[CollisionalRxnLUT::k53][i] * HII(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k55][i] * HI(i,j,k); + if (my_chemistry->use_radiative_transfer == 1) { + if (my_chemistry->radiative_transfer_HDI_dissociation > 0) { + acoef = acoef + + kdissHDI(i,j,k); + } + } + if (my_chemistry->primordial_chemistry > 3) { + scoef = scoef + 3. * ( 0. + + kcol_buf.data[CollisionalRxnLUT::k125][i] * HDII(i,j,k) * HI(i,j,k) / 3. + + kcol_buf.data[CollisionalRxnLUT::k137][i] * DM(i,j,k) * HI(i,j,k) / 2. + ); + } + out_spdens.data[SpLUT::HDI][i] = ( scoef*dtit[i] + HDI(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + } + } + } + + // --- (D2) Now do extra 3-species for minor primordial species --- + if (my_chemistry->primordial_chemistry > 3) { + for (i = idx_range.i_start; i < idx_range.i_stop; i++) { + if (itmask[i] != MASK_FALSE) { + + // 1) DM + scoef = + kcol_buf.data[CollisionalRxnLUT::k132][i] * DI(i,j,k) * de(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k135][i] * HM(i,j,k) * DI(i,j,k); + acoef = + kcol_buf.data[CollisionalRxnLUT::k133][i] * DII(i,j,k) / 2. + + kcol_buf.data[CollisionalRxnLUT::k134][i] * HII(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k136][i] * HI(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k137][i] * HI(i,j,k); + + out_spdens.data[SpLUT::DM][i] = ( scoef*dtit[i] + DM(i,j,k) ) / + ( 1. + acoef*dtit[i] ); + + // 2) HDII + scoef = 3. * ( + kcol_buf.data[CollisionalRxnLUT::k129][i] * DI(i,j,k) * HII(i,j,k) / 2. + + kcol_buf.data[CollisionalRxnLUT::k130][i] * DII(i,j,k) * HI(i,j,k) / 2. + ); + acoef = + kcol_buf.data[CollisionalRxnLUT::k125][i] * HI(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k131][i] * de(i,j,k); + + out_spdens.data[SpLUT::HDII][i] = ( scoef*dtit[i] + HDII(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + // 3) HeHII + scoef = 5. * ( + kcol_buf.data[CollisionalRxnLUT::k148][i] * HeI(i,j,k) * HII(i,j,k) / 4. + + kcol_buf.data[CollisionalRxnLUT::k149][i] * HeI(i,j,k) * HII(i,j,k) / 4. + + kcol_buf.data[CollisionalRxnLUT::k150][i] * HeI(i,j,k) * H2II(i,j,k) / 8. + + kcol_buf.data[CollisionalRxnLUT::k151][i] * HeII(i,j,k) * HI(i,j,k) / 4. + ); + acoef = + kcol_buf.data[CollisionalRxnLUT::k152][i] * HI(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::k153][i] * de(i,j,k); + + out_spdens.data[SpLUT::HeHII][i] = ( scoef*dtit[i] + HeHII(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + } + } + } + + // --- (D3) Now do metal species --- + if (my_chemistry->metal_chemistry == 1) { + for (i = idx_range.i_start; i < idx_range.i_stop; i++) { + if (itmask_metal[i] != MASK_FALSE) { + + // ***** CI ********** + scoef = 0. + 12. * ( 0. + + kcol_buf.data[CollisionalRxnLUT::kz15][i] * HI(i,j,k) * CH(i,j,k) / 13. + + kcol_buf.data[CollisionalRxnLUT::kz44][i] * CII(i,j,k) * de(i,j,k) / 12. + ); + acoef = 0. + + kcol_buf.data[CollisionalRxnLUT::kz20][i] * H2I(i,j,k) / 2. + + kcol_buf.data[CollisionalRxnLUT::kz27][i] * HI(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::kz28][i] * OH(i,j,k) / 17. + + kcol_buf.data[CollisionalRxnLUT::kz29][i] * O2(i,j,k) / 32. + + kcol_buf.data[CollisionalRxnLUT::kz51][i] * H2I(i,j,k) / 2.; + if ( ( my_chemistry->grain_growth == 1 ) || ( my_chemistry->dust_sublimation == 1) ) { + if (my_chemistry->dust_species > 0) { + acoef = acoef + + grain_growth_rates.data[OnlyGrainSpLUT::AC_dust] [i] / CI(i,j,k) * 12.; + } + } + if (my_chemistry->use_radiative_transfer == 1) { + if (my_chemistry->radiative_transfer_metal_ionization > 0) { + acoef = acoef + + kphCI(i,j,k); + } + if (my_chemistry->radiative_transfer_metal_dissociation > 0) { + scoef = scoef + 12. * + kdissCO (i,j,k) * CO(i,j,k) /28.0; + } + } + + out_spdens.data[SpLUT::CI][i] = ( scoef*dtit[i] + CI(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + // ***** CII ********** + scoef = 0. + 12. * ( 0. + ); + acoef = 0. + + kcol_buf.data[CollisionalRxnLUT::kz37][i] * OH(i,j,k) / 17. + + kcol_buf.data[CollisionalRxnLUT::kz38][i] * O2(i,j,k) / 32. + + kcol_buf.data[CollisionalRxnLUT::kz44][i] * de(i,j,k); + if (my_chemistry->use_radiative_transfer == 1) { + if (my_chemistry->radiative_transfer_metal_ionization > 0) { + scoef = scoef + + kphCI(i,j,k) * CI(i,j,k); + } + } + + out_spdens.data[SpLUT::CII][i] = ( scoef*dtit[i] + CII(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + // ***** CO ********** + scoef = 0. + 28. * ( 0. + + kcol_buf.data[CollisionalRxnLUT::kz28][i] * CI(i,j,k) * OH(i,j,k) / 204. + + kcol_buf.data[CollisionalRxnLUT::kz29][i] * CI(i,j,k) * O2(i,j,k) / 384. + + kcol_buf.data[CollisionalRxnLUT::kz32][i] * OI(i,j,k) * CH(i,j,k) / 208. + + kcol_buf.data[CollisionalRxnLUT::kz38][i] * CII(i,j,k) * O2(i,j,k) / 384. + + kcol_buf.data[CollisionalRxnLUT::kz43][i] * COII(i,j,k) * HI(i,j,k) / 28. + ); + acoef = 0. + + kcol_buf.data[CollisionalRxnLUT::kz26][i] * OH(i,j,k) / 17.; + if ( ( my_chemistry->grain_growth == 1 ) || ( my_chemistry->dust_sublimation == 1) ) { + if (my_chemistry->dust_species > 2) { + acoef = acoef + + grain_growth_rates.data[OnlyGrainSpLUT::ref_org_dust] [i] / CO(i,j,k) * 17. * 0.5 + + grain_growth_rates.data[OnlyGrainSpLUT::vol_org_dust] [i] / CO(i,j,k) * 17.; + } + } + if (my_chemistry->use_radiative_transfer == 1) { + if (my_chemistry->radiative_transfer_metal_dissociation > 0) { + acoef = acoef + + kdissCO (i,j,k); + } + } + + out_spdens.data[SpLUT::CO][i] = ( scoef*dtit[i] + CO(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + // ***** CO2 ********** + scoef = 0. + 44. * ( 0. + + kcol_buf.data[CollisionalRxnLUT::kz26][i] * OH(i,j,k) * CO(i,j,k) / 476. + ); + acoef = 0.; + + out_spdens.data[SpLUT::CO2][i] = ( scoef*dtit[i] + CO2(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + // ***** OI ********** + scoef = 0. + 16. * ( 0. + + kcol_buf.data[CollisionalRxnLUT::kz17][i] * HI(i,j,k) * OH(i,j,k) / 17. + + kcol_buf.data[CollisionalRxnLUT::kz19][i] * HI(i,j,k) * O2(i,j,k) / 32. + + kcol_buf.data[CollisionalRxnLUT::kz25][i] * OH(i,j,k) * OH(i,j,k) / 289. + + kcol_buf.data[CollisionalRxnLUT::kz29][i] * CI(i,j,k) * O2(i,j,k) / 384. + + kcol_buf.data[CollisionalRxnLUT::kz39][i] * OII(i,j,k) * HI(i,j,k) / 16. + + kcol_buf.data[CollisionalRxnLUT::kz45][i] * OII(i,j,k) * de(i,j,k) / 16. + + kcol_buf.data[CollisionalRxnLUT::kz47][i] * H2OII(i,j,k) * de(i,j,k) / 18. + + kcol_buf.data[CollisionalRxnLUT::kz50][i] * O2II(i,j,k) * de(i,j,k) / 16. + + kcol_buf.data[CollisionalRxnLUT::kz53][i] * SiI(i,j,k) * O2(i,j,k) / 896. + ); + acoef = 0. + + kcol_buf.data[CollisionalRxnLUT::kz21][i] * H2I(i,j,k) / 2. + + kcol_buf.data[CollisionalRxnLUT::kz22][i] * HII(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::kz30][i] * HI(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::kz31][i] * OI(i,j,k) / 8. + + kcol_buf.data[CollisionalRxnLUT::kz32][i] * CH(i,j,k) / 13. + + kcol_buf.data[CollisionalRxnLUT::kz33][i] * OH(i,j,k) / 17.; + if (my_chemistry->use_radiative_transfer == 1) { + if (my_chemistry->radiative_transfer_metal_ionization > 0) { + acoef = acoef + + kphOI(i,j,k); + } + if (my_chemistry->radiative_transfer_metal_dissociation > 0) { + scoef = scoef + 16. * + ( kdissOH (i,j,k) * OH(i,j,k) /17.0 + + kdissCO (i,j,k) * CO(i,j,k) /28.0); + } + } + + out_spdens.data[SpLUT::OI][i] = ( scoef*dtit[i] + OI(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + // ***** OH ********** + scoef = 0. + 17. * ( 0. + + kcol_buf.data[CollisionalRxnLUT::kz18][i] * HI(i,j,k) * H2O(i,j,k) / 18. + + kcol_buf.data[CollisionalRxnLUT::kz19][i] * HI(i,j,k) * O2(i,j,k) / 32. + + kcol_buf.data[CollisionalRxnLUT::kz21][i] * OI(i,j,k) * H2I(i,j,k) / 32. + + kcol_buf.data[CollisionalRxnLUT::kz30][i] * OI(i,j,k) * HI(i,j,k) / 16. + + kcol_buf.data[CollisionalRxnLUT::kz46][i] * H2OII(i,j,k) * de(i,j,k) / 18. + + kcol_buf.data[CollisionalRxnLUT::kz49][i] * H3OII(i,j,k) * de(i,j,k) / 19. + ); + acoef = 0. + + kcol_buf.data[CollisionalRxnLUT::kz17][i] * HI(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::kz24][i] * H2I(i,j,k) / 2. + + kcol_buf.data[CollisionalRxnLUT::kz25][i] * OH(i,j,k) / 8.5 + + kcol_buf.data[CollisionalRxnLUT::kz26][i] * CO(i,j,k) / 28. + + kcol_buf.data[CollisionalRxnLUT::kz28][i] * CI(i,j,k) / 12. + + kcol_buf.data[CollisionalRxnLUT::kz33][i] * OI(i,j,k) / 16. + + kcol_buf.data[CollisionalRxnLUT::kz34][i] * HII(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::kz37][i] * CII(i,j,k) / 12. + + kcol_buf.data[CollisionalRxnLUT::kz52][i] * SiI(i,j,k) / 28. + + kcol_buf.data[CollisionalRxnLUT::kz54][i] * SiOI(i,j,k) / 44.; + if (my_chemistry->use_radiative_transfer == 1) { + if (my_chemistry->radiative_transfer_metal_dissociation > 0) { + acoef = acoef + + kdissOH (i,j,k); + scoef = scoef + 17. * + kdissH2O(i,j,k) * H2O(i,j,k)/18.0; + } + } + + out_spdens.data[SpLUT::OH][i] = ( scoef*dtit[i] + OH(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + // ***** H2O ********** + scoef = 0. + 18. * ( 0. + + kcol_buf.data[CollisionalRxnLUT::kz24][i] * H2I(i,j,k) * OH(i,j,k) / 34. + + kcol_buf.data[CollisionalRxnLUT::kz25][i] * OH(i,j,k) * OH(i,j,k) / 289. + + kcol_buf.data[CollisionalRxnLUT::kz48][i] * H3OII(i,j,k) * de(i,j,k) / 19. + ); + acoef = 0. + + kcol_buf.data[CollisionalRxnLUT::kz18][i] * HI(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::kz35][i] * HII(i,j,k); + if ( ( my_chemistry->grain_growth == 1 ) || ( my_chemistry->dust_sublimation == 1) ) { + if (my_chemistry->dust_species > 0) { + acoef = acoef + + grain_growth_rates.data[OnlyGrainSpLUT::MgSiO3_dust] [i] / H2O(i,j,k) * 18. * 2.; + } + if (my_chemistry->dust_species > 1) { + acoef = acoef + + grain_growth_rates.data[OnlyGrainSpLUT::Mg2SiO4_dust] [i] / H2O(i,j,k) * 18. * 3. + + grain_growth_rates.data[OnlyGrainSpLUT::Fe3O4_dust] [i] / H2O(i,j,k) * 18. * 4. + + grain_growth_rates.data[OnlyGrainSpLUT::MgO_dust] [i] / H2O(i,j,k) * 18. + + grain_growth_rates.data[OnlyGrainSpLUT::Al2O3_dust] [i] / H2O(i,j,k) * 18. * 3.; + } + if (my_chemistry->dust_species > 2) { + acoef = acoef + + grain_growth_rates.data[OnlyGrainSpLUT::H2O_ice_dust] [i] / H2O(i,j,k) * 18.; + } + } + if (my_chemistry->use_radiative_transfer == 1) { + if (my_chemistry->radiative_transfer_metal_dissociation > 0) { + acoef = acoef + + kdissH2O(i,j,k); + } + } + + out_spdens.data[SpLUT::H2O][i] = ( scoef*dtit[i] + H2O(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + // ***** O2 ********** + scoef = 0. + 32. * ( 0. + + kcol_buf.data[CollisionalRxnLUT::kz31][i] * OI(i,j,k) * OI(i,j,k) / 256. + + kcol_buf.data[CollisionalRxnLUT::kz33][i] * OI(i,j,k) * OH(i,j,k) / 272. + ); + acoef = 0. + + kcol_buf.data[CollisionalRxnLUT::kz19][i] * HI(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::kz29][i] * CI(i,j,k) / 12. + + kcol_buf.data[CollisionalRxnLUT::kz36][i] * HII(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::kz38][i] * CII(i,j,k) / 12. + + kcol_buf.data[CollisionalRxnLUT::kz53][i] * SiI(i,j,k) / 28.; + + out_spdens.data[SpLUT::O2][i] = ( scoef*dtit[i] + O2(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + // ***** SiI ********** + scoef = 0. + 28. * ( 0. + ); + acoef = 0. + + kcol_buf.data[CollisionalRxnLUT::kz52][i] * OH(i,j,k) / 17. + + kcol_buf.data[CollisionalRxnLUT::kz53][i] * O2(i,j,k) / 32.; + if ( ( my_chemistry->grain_growth == 1 ) || ( my_chemistry->dust_sublimation == 1) ) { + if (my_chemistry->dust_species > 1) { + acoef = acoef + + grain_growth_rates.data[OnlyGrainSpLUT::SiM_dust] [i] / SiI(i,j,k) * 28.; + } + } + + out_spdens.data[SpLUT::SiI][i] = ( scoef*dtit[i] + SiI(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + // ***** SiOI ********** + scoef = 0. + 44. * ( 0. + + kcol_buf.data[CollisionalRxnLUT::kz52][i] * SiI(i,j,k) * OH(i,j,k) / 476. + + kcol_buf.data[CollisionalRxnLUT::kz53][i] * SiI(i,j,k) * O2(i,j,k) / 896. + ); + acoef = 0. + + kcol_buf.data[CollisionalRxnLUT::kz54][i] * OH(i,j,k) / 17.; + if ( ( my_chemistry->grain_growth == 1 ) || ( my_chemistry->dust_sublimation == 1) ) { + if (my_chemistry->dust_species > 0) { + acoef = acoef + + grain_growth_rates.data[OnlyGrainSpLUT::MgSiO3_dust] [i] / SiOI(i,j,k) * 44.; + } + if (my_chemistry->dust_species > 1) { + acoef = acoef + + grain_growth_rates.data[OnlyGrainSpLUT::Mg2SiO4_dust] [i] / SiOI(i,j,k) * 44.; + } + } + + out_spdens.data[SpLUT::SiOI][i] = ( scoef*dtit[i] + SiOI(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + // ***** SiO2I ********** + scoef = 0. + 60. * ( 0. + + kcol_buf.data[CollisionalRxnLUT::kz54][i] * SiOI(i,j,k) * OH(i,j,k) / 748. + ); + acoef = 0.; + if ( ( my_chemistry->grain_growth == 1 ) || ( my_chemistry->dust_sublimation == 1) ) { + if (my_chemistry->dust_species > 1) { + acoef = acoef + + grain_growth_rates.data[OnlyGrainSpLUT::SiO2_dust] [i] / SiO2I(i,j,k) * 60.; + } + } + + out_spdens.data[SpLUT::SiO2I][i] = ( scoef*dtit[i] + SiO2I(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + // MINOR BUT IMPORTANT SPECIES FOR MOLECULAR FORMATION + //- ***** CH ********** + scoef = 0. + 13. * ( 0. + + kcol_buf.data[CollisionalRxnLUT::kz16][i] * HI(i,j,k) * CH2(i,j,k) / 14. + + kcol_buf.data[CollisionalRxnLUT::kz20][i] * CI(i,j,k) * H2I(i,j,k) / 24. + + kcol_buf.data[CollisionalRxnLUT::kz27][i] * CI(i,j,k) * HI(i,j,k) / 12. + ); + acoef = 0. + + kcol_buf.data[CollisionalRxnLUT::kz15][i] * HI(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::kz23][i] * H2I(i,j,k) / 2. + + kcol_buf.data[CollisionalRxnLUT::kz32][i] * OI(i,j,k) / 16.; + + out_spdens.data[SpLUT::CH][i] = ( scoef*dtit[i] + CH(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + // ***** CH2 ********** + scoef = 0. + 14. * ( 0. + + kcol_buf.data[CollisionalRxnLUT::kz23][i] * H2I(i,j,k) * CH(i,j,k) / 26. + + kcol_buf.data[CollisionalRxnLUT::kz51][i] * H2I(i,j,k) * CI(i,j,k) / 24. + ); + acoef = 0. + + kcol_buf.data[CollisionalRxnLUT::kz16][i] * HI(i,j,k); + if ( ( my_chemistry->grain_growth == 1 ) || ( my_chemistry->dust_sublimation == 1) ) { + if (my_chemistry->dust_species > 2) { + acoef = acoef + + grain_growth_rates.data[OnlyGrainSpLUT::ref_org_dust] [i] / CH2(i,j,k) * 14. * 0.5; + } + } + + out_spdens.data[SpLUT::CH2][i] = ( scoef*dtit[i] + CH2(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + // ***** COII ********** + scoef = 0. + 28. * ( 0. + + kcol_buf.data[CollisionalRxnLUT::kz37][i] * CII(i,j,k) * OH(i,j,k) / 204. + ); + acoef = 0. + + kcol_buf.data[CollisionalRxnLUT::kz43][i] * HI(i,j,k); + + out_spdens.data[SpLUT::COII][i] = ( scoef*dtit[i] + COII(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + // ***** OII ********** + scoef = 0. + 16. * ( 0. + + kcol_buf.data[CollisionalRxnLUT::kz22][i] * HII(i,j,k) * OI(i,j,k) / 16. + + kcol_buf.data[CollisionalRxnLUT::kz38][i] * CII(i,j,k) * O2(i,j,k) / 384. + ); + acoef = 0. + + kcol_buf.data[CollisionalRxnLUT::kz39][i] * HI(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::kz40][i] * H2I(i,j,k) / 2. + + kcol_buf.data[CollisionalRxnLUT::kz45][i] * de(i,j,k); + if (my_chemistry->use_radiative_transfer == 1) { + if (my_chemistry->radiative_transfer_metal_ionization > 0) { + scoef = scoef + + kphOI(i,j,k) * OI(i,j,k); + } + } + + out_spdens.data[SpLUT::OII][i] = ( scoef*dtit[i] + OII(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + // ***** OHII ********** + scoef = 0. + 17. * ( 0. + + kcol_buf.data[CollisionalRxnLUT::kz34][i] * HII(i,j,k) * OH(i,j,k) / 17. + + kcol_buf.data[CollisionalRxnLUT::kz40][i] * OII(i,j,k) * H2I(i,j,k) / 32. + ); + acoef = 0. + + kcol_buf.data[CollisionalRxnLUT::kz41][i] * H2I(i,j,k) / 2.; + + out_spdens.data[SpLUT::OHII][i] = ( scoef*dtit[i] + OHII(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + // ***** H2OII ********** + scoef = 0. + 18. * ( 0. + + kcol_buf.data[CollisionalRxnLUT::kz35][i] * HII(i,j,k) * H2O(i,j,k) / 18. + + kcol_buf.data[CollisionalRxnLUT::kz41][i] * OHII(i,j,k) * H2I(i,j,k) / 34. + ); + acoef = 0. + + kcol_buf.data[CollisionalRxnLUT::kz42][i] * H2I(i,j,k) / 2. + + kcol_buf.data[CollisionalRxnLUT::kz46][i] * de(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::kz47][i] * de(i,j,k); + + out_spdens.data[SpLUT::H2OII][i] = ( scoef*dtit[i] + H2OII(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + // ***** H3OII ********** + scoef = 0. + 19. * ( 0. + + kcol_buf.data[CollisionalRxnLUT::kz42][i] * H2OII(i,j,k) * H2I(i,j,k) / 36. + ); + acoef = 0. + + kcol_buf.data[CollisionalRxnLUT::kz48][i] * de(i,j,k) + + kcol_buf.data[CollisionalRxnLUT::kz49][i] * de(i,j,k); + + out_spdens.data[SpLUT::H3OII][i] = ( scoef*dtit[i] + H3OII(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + // ***** O2II ********** + scoef = 0. + 32. * ( 0. + + kcol_buf.data[CollisionalRxnLUT::kz36][i] * HII(i,j,k) * O2(i,j,k) / 32. + ); + acoef = 0. + + kcol_buf.data[CollisionalRxnLUT::kz50][i] * de(i,j,k); + + out_spdens.data[SpLUT::O2II][i] = ( scoef*dtit[i] + O2II(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + if ( ( my_chemistry->grain_growth == 1 ) || ( my_chemistry->dust_sublimation == 1) ) { + if (my_chemistry->dust_species > 0) { + // ***** Mg ********** + scoef = 0.; + acoef = 0.; + acoef = acoef + + grain_growth_rates.data[OnlyGrainSpLUT::MgSiO3_dust] [i] / Mg(i,j,k) * 24.; + if (my_chemistry->dust_species > 1) { + acoef = acoef + + grain_growth_rates.data[OnlyGrainSpLUT::Mg2SiO4_dust] [i] / Mg(i,j,k) * 24. * 2. + + grain_growth_rates.data[OnlyGrainSpLUT::MgO_dust] [i] / Mg(i,j,k) * 24.; + } + + out_spdens.data[SpLUT::Mg][i] = ( scoef*dtit[i] + Mg(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + } + + if (my_chemistry->dust_species > 1) { + // ***** Al ********** + scoef = 0.; + acoef = 0.; + acoef = acoef + + grain_growth_rates.data[OnlyGrainSpLUT::Al2O3_dust] [i] / Al(i,j,k) * 27. * 2.; + + out_spdens.data[SpLUT::Al][i] = ( scoef*dtit[i] + Al(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + // ***** S ********** + scoef = 0.; + acoef = 0.; + acoef = acoef + + grain_growth_rates.data[OnlyGrainSpLUT::FeS_dust] [i] / S(i,j,k) * 32.; + + out_spdens.data[SpLUT::S][i] = ( scoef*dtit[i] + S(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + // ***** Fe ********** + scoef = 0.; + acoef = 0.; + acoef = acoef + + grain_growth_rates.data[OnlyGrainSpLUT::FeM_dust] [i] / Fe(i,j,k) * 56. + + grain_growth_rates.data[OnlyGrainSpLUT::Fe3O4_dust] [i] / Fe(i,j,k) * 56. * 3. + + grain_growth_rates.data[OnlyGrainSpLUT::FeS_dust] [i] / Fe(i,j,k) * 56.; + + out_spdens.data[SpLUT::Fe][i] = ( scoef*dtit[i] + Fe(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + } + } + + } + } + } + + // --- (D4) Now do dust species --- + if ( ( my_chemistry->grain_growth == 1 ) || ( my_chemistry->dust_sublimation == 1) ) { + for (i = idx_range.i_start; i < idx_range.i_stop; i++) { + if (itmask_metal[i] != MASK_FALSE) { + + if (my_chemistry->dust_species > 0) { + // ***** MgSiO3 ********** + scoef = 0.; + scoef = scoef + + grain_growth_rates.data[OnlyGrainSpLUT::MgSiO3_dust] [i] * 100.; + acoef = 0.; + + out_spdens.data[SpLUT::MgSiO3_dust][i] = ( scoef*dtit[i] + MgSiO3(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + // ***** AC ********** + scoef = 0.; + scoef = scoef + + grain_growth_rates.data[OnlyGrainSpLUT::AC_dust] [i] * 12.; + acoef = 0.; + + out_spdens.data[SpLUT::AC_dust][i] = ( scoef*dtit[i] + AC(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + } + + if (my_chemistry->dust_species > 1) { + // ***** SiM ********** + scoef = 0.; + scoef = scoef + + grain_growth_rates.data[OnlyGrainSpLUT::SiM_dust] [i] * 28.; + acoef = 0.; + + out_spdens.data[SpLUT::SiM_dust][i] = ( scoef*dtit[i] + SiM(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + // ***** FeM ********** + scoef = 0.; + scoef = scoef + + grain_growth_rates.data[OnlyGrainSpLUT::FeM_dust] [i] * 56.; + acoef = 0.; + + out_spdens.data[SpLUT::FeM_dust][i] = ( scoef*dtit[i] + FeM(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + // ***** Mg2SiO4 ********** + scoef = 0.; + scoef = scoef + + grain_growth_rates.data[OnlyGrainSpLUT::Mg2SiO4_dust] [i] * 140.; + acoef = 0.; + + out_spdens.data[SpLUT::Mg2SiO4_dust][i] = ( scoef*dtit[i] + Mg2SiO4(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + // ***** Fe3O4 ********** + scoef = 0.; + scoef = scoef + + grain_growth_rates.data[OnlyGrainSpLUT::Fe3O4_dust] [i] * 232.; + acoef = 0.; + + out_spdens.data[SpLUT::Fe3O4_dust][i] = ( scoef*dtit[i] + Fe3O4(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + // ***** SiO2D ********** + scoef = 0.; + scoef = scoef + + grain_growth_rates.data[OnlyGrainSpLUT::SiO2_dust] [i] * 60.; + acoef = 0.; + + out_spdens.data[SpLUT::SiO2_dust][i] = ( scoef*dtit[i] + SiO2D(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + // ***** MgO ********** + scoef = 0.; + scoef = scoef + + grain_growth_rates.data[OnlyGrainSpLUT::MgO_dust] [i] * 40.; + acoef = 0.; + + out_spdens.data[SpLUT::MgO_dust][i] = ( scoef*dtit[i] + MgO(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + // ***** FeS ********** + scoef = 0.; + scoef = scoef + + grain_growth_rates.data[OnlyGrainSpLUT::FeS_dust] [i] * 88.; + acoef = 0.; + + out_spdens.data[SpLUT::FeS_dust][i] = ( scoef*dtit[i] + FeS(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + // ***** Al2O3 ********** + scoef = 0.; + scoef = scoef + + grain_growth_rates.data[OnlyGrainSpLUT::Al2O3_dust] [i] * 102.; + acoef = 0.; + + out_spdens.data[SpLUT::Al2O3_dust][i] = ( scoef*dtit[i] + Al2O3(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + } + + if (my_chemistry->dust_species > 2) { + // ***** reforg ********** + scoef = 0.; + scoef = scoef + + grain_growth_rates.data[OnlyGrainSpLUT::ref_org_dust] [i] * 22.68; + acoef = 0.; + + out_spdens.data[SpLUT::ref_org_dust][i] = ( scoef*dtit[i] + reforg(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + // ***** volorg ********** + scoef = 0.; + scoef = scoef + + grain_growth_rates.data[OnlyGrainSpLUT::vol_org_dust] [i] * 32.; + acoef = 0.; + + out_spdens.data[SpLUT::vol_org_dust][i] = ( scoef*dtit[i] + volorg(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + + // ***** H2Oice ********** + scoef = 0.; + scoef = scoef + + grain_growth_rates.data[OnlyGrainSpLUT::H2O_ice_dust] [i] * 18.; + acoef = 0.; + + out_spdens.data[SpLUT::H2O_ice_dust][i] = ( scoef*dtit[i] + H2Oice(i,j,k) ) + / ( 1. + acoef*dtit[i] ); + + } + + } + } + } +} + + + /// Computes the time derivatives of the mass densities of various species /// /// This function is specialized to operate on a single cell at a time. This is @@ -1380,6 +2790,8 @@ inline void species_density_derivatives_0d( } + + } // namespace grackle::impl::chemistry #endif /* CHEMISTRY_SOLVER_FUNCS_HPP */ diff --git a/src/clib/cool1d_multi_g.cpp b/src/clib/cool1d_multi_g.cpp index 02793cc90..301052408 100644 --- a/src/clib/cool1d_multi_g.cpp +++ b/src/clib/cool1d_multi_g.cpp @@ -959,7 +959,6 @@ void grackle::impl::cool1d_multi_g( if (itmask[i] != MASK_FALSE) { // Only calculate if H2I(i) is a substantial fraction if (d(i, idx_range.j, idx_range.k) * dom > 1e10) { - ciefudge = 1.; tau = std::pow(((d(i, idx_range.j, idx_range.k) / 2e16) * dom), 2.8); // 2e16 is in units of cm^-3 tau = std::fmax(tau, 1.e-5); @@ -1969,4 +1968,4 @@ void grackle::impl::cool1d_multi_g( grackle::impl::drop_GrainSpeciesCollection(&gas_grainsp_heatrate); return; -} \ No newline at end of file +} diff --git a/src/clib/fortran_func_decls.h b/src/clib/fortran_func_decls.h index ff6fd5755..3d885a126 100644 --- a/src/clib/fortran_func_decls.h +++ b/src/clib/fortran_func_decls.h @@ -266,70 +266,6 @@ void FORTRAN_NAME(rate_timestep_g)( gr_float* kdissH2O_data_ptr ); -void FORTRAN_NAME(step_rate_g)( - gr_float* de_data_ptr, gr_float* HI_data_ptr, gr_float* HII_data_ptr, - gr_float* HeI_data_ptr, gr_float* HeII_data_ptr, gr_float* HeIII_data_ptr, - gr_float* d_data_ptr, gr_float* HM_data_ptr, gr_float* H2I_data_ptr, - gr_float* H2II_data_ptr, gr_float* DI_data_ptr, gr_float* DII_data_ptr, - gr_float* HDI_data_ptr, double* dtit, int* in, int* jn, int* kn, int* is, - int* ie, int* j, int* k, int* ispecies, gr_mask_type* anydust, double* k1, - double* k2, double* k3, double* k4, double* k5, double* k6, double* k7, - double* k8, double* k9, double* k10, double* k11, double* k12, double* k13, - double* k14, double* k15, double* k16, double* k17, double* k18, double* k19, - double* k22, double* k24, double* k25, double* k26, double* k27, double* k28, - double* k29, double* k30, double* k50, double* k51, double* k52, double* k53, - double* k54, double* k55, double* k56, double* k57, double* k58, - double* h2dust, double* rhoH, double* k24shield, double* k25shield, - double* k26shield, double* k28shield, double* k29shield, double* k30shield, - double* k31shield, double* HIp, double* HIIp, double* HeIp, double* HeIIp, - double* HeIIIp, double* dep, double* HMp, double* H2Ip, double* H2IIp, - double* DIp, double* DIIp, double* HDIp, double* dedot_prev, - double* HIdot_prev, int* iradtrans, int* irt_honly, gr_float* kphHI_data_ptr, - gr_float* kphHeI_data_ptr, gr_float* kphHeII_data_ptr, gr_mask_type* itmask, - gr_mask_type* itmask_metal, gr_float* DM_data_ptr, gr_float* HDII_data_ptr, - gr_float* HeHII_data_ptr, int* imetal, gr_float* metal_data_ptr, int* imchem, - int* idspecies, int* igrgr, int* idsub, gr_float* CI_data_ptr, - gr_float* CII_data_ptr, gr_float* CO_data_ptr, gr_float* CO2_data_ptr, - gr_float* OI_data_ptr, gr_float* OH_data_ptr, gr_float* H2O_data_ptr, - gr_float* O2_data_ptr, gr_float* SiI_data_ptr, gr_float* SiOI_data_ptr, - gr_float* SiO2I_data_ptr, gr_float* CH_data_ptr, gr_float* CH2_data_ptr, - gr_float* COII_data_ptr, gr_float* OII_data_ptr, gr_float* OHII_data_ptr, - gr_float* H2OII_data_ptr, gr_float* H3OII_data_ptr, gr_float* O2II_data_ptr, - gr_float* Mg_data_ptr, gr_float* Al_data_ptr, gr_float* S_data_ptr, - gr_float* Fe_data_ptr, gr_float* SiM_data_ptr, gr_float* FeM_data_ptr, - gr_float* Mg2SiO4_data_ptr, gr_float* MgSiO3_data_ptr, - gr_float* Fe3O4_data_ptr, gr_float* AC_data_ptr, gr_float* SiO2D_data_ptr, - gr_float* MgO_data_ptr, gr_float* FeS_data_ptr, gr_float* Al2O3_data_ptr, - gr_float* reforg_data_ptr, gr_float* volorg_data_ptr, - gr_float* H2Oice_data_ptr, double* k125, double* k129, double* k130, - double* k131, double* k132, double* k133, double* k134, double* k135, - double* k136, double* k137, double* k148, double* k149, double* k150, - double* k151, double* k152, double* k153, double* kz15, double* kz16, - double* kz17, double* kz18, double* kz19, double* kz20, double* kz21, - double* kz22, double* kz23, double* kz24, double* kz25, double* kz26, - double* kz27, double* kz28, double* kz29, double* kz30, double* kz31, - double* kz32, double* kz33, double* kz34, double* kz35, double* kz36, - double* kz37, double* kz38, double* kz39, double* kz40, double* kz41, - double* kz42, double* kz43, double* kz44, double* kz45, double* kz46, - double* kz47, double* kz48, double* kz49, double* kz50, double* kz51, - double* kz52, double* kz53, double* kz54, double* DMp, double* HDIIp, - double* HeHIIp, double* CIp, double* CIIp, double* COp, double* CO2p, - double* OIp, double* OHp, double* H2Op, double* O2p, double* SiIp, - double* SiOIp, double* SiO2Ip, double* CHp, double* CH2p, double* COIIp, - double* OIIp, double* OHIIp, double* H2OIIp, double* H3OIIp, double* O2IIp, - double* Mgp, double* Alp, double* Sp, double* Fep, gr_float* SiMp, - gr_float* FeMp, gr_float* Mg2SiO4p, gr_float* MgSiO3p, gr_float* Fe3O4p, - gr_float* ACp, gr_float* SiO2Dp, gr_float* MgOp, gr_float* FeSp, - gr_float* Al2O3p, gr_float* reforgp, gr_float* volorgp, gr_float* H2Oicep, - double* kdSiM, double* kdFeM, double* kdMg2SiO4, double* kdMgSiO3, - double* kdFe3O4, double* kdAC, double* kdSiO2D, double* kdMgO, double* kdFeS, - double* kdAl2O3, double* kdreforg, double* kdvolorg, double* kdH2Oice, - int* idissHDI, gr_float* kdissHDI_data_ptr, int* iionZ, - gr_float* kphCI_data_ptr, gr_float* kphOI_data_ptr, int* idissZ, - gr_float* kdissCO_data_ptr, gr_float* kdissOH_data_ptr, - gr_float* kdissH2O_data_ptr -); - void FORTRAN_NAME(make_consistent_g)( gr_float* de_data_ptr, gr_float* HI_data_ptr, gr_float* HII_data_ptr, gr_float* HeI_data_ptr, gr_float* HeII_data_ptr, gr_float* HeIII_data_ptr, diff --git a/src/clib/fortran_func_wrappers.hpp b/src/clib/fortran_func_wrappers.hpp index e311723a3..6cf0b6b1b 100644 --- a/src/clib/fortran_func_wrappers.hpp +++ b/src/clib/fortran_func_wrappers.hpp @@ -31,6 +31,8 @@ #include "LUT.hpp" #include "utils-cpp.hpp" +#include "step_rate_gauss_seidel.hpp" + // callers of these functions are generally expected to locally shorten the // namespace name when they call these routines namespace grackle::impl::fortran_wrapper { @@ -370,79 +372,6 @@ inline void rate_timestep_g( ); } -/// Uses one linearly implicit Gauss-Seidel sweep of a backward-Euler time -/// integrator to advance the rate equations by one (sub-)cycle (dtit). -inline void step_rate_g( - double* dtit, IndexRange idx_range, gr_mask_type anydust, double* h2dust, - double* rhoH, double* dedot_prev, double* HIdot_prev, - gr_mask_type* itmask, gr_mask_type* itmask_metal, int imetal, - chemistry_data* my_chemistry, grackle_field_data* my_fields, - photo_rate_storage my_uvb_rates, - grackle::impl::GrainSpeciesCollection grain_growth_rates, - grackle::impl::SpeciesCollection species_tmpdens, - grackle::impl::CollisionalRxnRateCollection kcr_buf, - grackle::impl::PhotoRxnRateCollection kshield_buf -) { - FORTRAN_NAME(step_rate_g)(my_fields->e_density, my_fields->HI_density, my_fields->HII_density, my_fields->HeI_density, my_fields->HeII_density, my_fields->HeIII_density, my_fields->density, - my_fields->HM_density, my_fields->H2I_density, my_fields->H2II_density, my_fields->DI_density, my_fields->DII_density, my_fields->HDI_density, dtit, - &my_fields->grid_dimension[0], &my_fields->grid_dimension[1], &my_fields->grid_dimension[2], &idx_range.i_start, &idx_range.i_end, &idx_range.jp1, &idx_range.kp1, - &my_chemistry->primordial_chemistry, &anydust, - kcr_buf.data[CollisionalRxnLUT::k1], kcr_buf.data[CollisionalRxnLUT::k2], kcr_buf.data[CollisionalRxnLUT::k3], kcr_buf.data[CollisionalRxnLUT::k4], kcr_buf.data[CollisionalRxnLUT::k5], kcr_buf.data[CollisionalRxnLUT::k6], kcr_buf.data[CollisionalRxnLUT::k7], kcr_buf.data[CollisionalRxnLUT::k8], kcr_buf.data[CollisionalRxnLUT::k9], kcr_buf.data[CollisionalRxnLUT::k10], kcr_buf.data[CollisionalRxnLUT::k11], - kcr_buf.data[CollisionalRxnLUT::k12], kcr_buf.data[CollisionalRxnLUT::k13], kcr_buf.data[CollisionalRxnLUT::k14], kcr_buf.data[CollisionalRxnLUT::k15], kcr_buf.data[CollisionalRxnLUT::k16], kcr_buf.data[CollisionalRxnLUT::k17], kcr_buf.data[CollisionalRxnLUT::k18], kcr_buf.data[CollisionalRxnLUT::k19], kcr_buf.data[CollisionalRxnLUT::k22], - &my_uvb_rates.k24, &my_uvb_rates.k25, &my_uvb_rates.k26, &my_uvb_rates.k27, &my_uvb_rates.k28, &my_uvb_rates.k29, &my_uvb_rates.k30, - kcr_buf.data[CollisionalRxnLUT::k50], kcr_buf.data[CollisionalRxnLUT::k51], kcr_buf.data[CollisionalRxnLUT::k52], kcr_buf.data[CollisionalRxnLUT::k53], kcr_buf.data[CollisionalRxnLUT::k54], kcr_buf.data[CollisionalRxnLUT::k55], kcr_buf.data[CollisionalRxnLUT::k56], kcr_buf.data[CollisionalRxnLUT::k57], kcr_buf.data[CollisionalRxnLUT::k58], - h2dust, rhoH, - kshield_buf.k24, kshield_buf.k25, kshield_buf.k26, - kshield_buf.k28, kshield_buf.k29, kshield_buf.k30, kshield_buf.k31, - species_tmpdens.data[SpLUT::HI], species_tmpdens.data[SpLUT::HII], species_tmpdens.data[SpLUT::HeI], species_tmpdens.data[SpLUT::HeII], species_tmpdens.data[SpLUT::HeIII], species_tmpdens.data[SpLUT::e], - species_tmpdens.data[SpLUT::HM], species_tmpdens.data[SpLUT::H2I], species_tmpdens.data[SpLUT::H2II], species_tmpdens.data[SpLUT::DI], species_tmpdens.data[SpLUT::DII], species_tmpdens.data[SpLUT::HDI], - dedot_prev, HIdot_prev, - &my_chemistry->use_radiative_transfer, &my_chemistry->radiative_transfer_hydrogen_only, - my_fields->RT_HI_ionization_rate, my_fields->RT_HeI_ionization_rate, my_fields->RT_HeII_ionization_rate, - itmask, itmask_metal, - my_fields->DM_density, my_fields->HDII_density, my_fields->HeHII_density, &imetal, my_fields->metal_density, - &my_chemistry->metal_chemistry, &my_chemistry->dust_species, &my_chemistry->grain_growth, &my_chemistry->dust_sublimation, - my_fields->CI_density, my_fields->CII_density, my_fields->CO_density, my_fields->CO2_density, - my_fields->OI_density, my_fields->OH_density, my_fields->H2O_density, my_fields->O2_density, - my_fields->SiI_density, my_fields->SiOI_density, my_fields->SiO2I_density, - my_fields->CH_density, my_fields->CH2_density, my_fields->COII_density, my_fields->OII_density, - my_fields->OHII_density, my_fields->H2OII_density, my_fields->H3OII_density, my_fields->O2II_density, - my_fields->Mg_density, my_fields->Al_density, my_fields->S_density, my_fields->Fe_density, - my_fields->SiM_dust_density, my_fields->FeM_dust_density, my_fields->Mg2SiO4_dust_density, my_fields->MgSiO3_dust_density, my_fields->Fe3O4_dust_density, - my_fields->AC_dust_density, my_fields->SiO2_dust_density, my_fields->MgO_dust_density, my_fields->FeS_dust_density, my_fields->Al2O3_dust_density, - my_fields->ref_org_dust_density, my_fields->vol_org_dust_density, my_fields->H2O_ice_dust_density, - kcr_buf.data[CollisionalRxnLUT::k125], kcr_buf.data[CollisionalRxnLUT::k129], kcr_buf.data[CollisionalRxnLUT::k130], kcr_buf.data[CollisionalRxnLUT::k131], kcr_buf.data[CollisionalRxnLUT::k132], - kcr_buf.data[CollisionalRxnLUT::k133], kcr_buf.data[CollisionalRxnLUT::k134], kcr_buf.data[CollisionalRxnLUT::k135], kcr_buf.data[CollisionalRxnLUT::k136], kcr_buf.data[CollisionalRxnLUT::k137], - kcr_buf.data[CollisionalRxnLUT::k148], kcr_buf.data[CollisionalRxnLUT::k149], kcr_buf.data[CollisionalRxnLUT::k150], kcr_buf.data[CollisionalRxnLUT::k151], kcr_buf.data[CollisionalRxnLUT::k152], - kcr_buf.data[CollisionalRxnLUT::k153], - kcr_buf.data[CollisionalRxnLUT::kz15], kcr_buf.data[CollisionalRxnLUT::kz16], kcr_buf.data[CollisionalRxnLUT::kz17], kcr_buf.data[CollisionalRxnLUT::kz18], kcr_buf.data[CollisionalRxnLUT::kz19], - kcr_buf.data[CollisionalRxnLUT::kz20], kcr_buf.data[CollisionalRxnLUT::kz21], kcr_buf.data[CollisionalRxnLUT::kz22], kcr_buf.data[CollisionalRxnLUT::kz23], kcr_buf.data[CollisionalRxnLUT::kz24], - kcr_buf.data[CollisionalRxnLUT::kz25], kcr_buf.data[CollisionalRxnLUT::kz26], kcr_buf.data[CollisionalRxnLUT::kz27], kcr_buf.data[CollisionalRxnLUT::kz28], kcr_buf.data[CollisionalRxnLUT::kz29], - kcr_buf.data[CollisionalRxnLUT::kz30], kcr_buf.data[CollisionalRxnLUT::kz31], kcr_buf.data[CollisionalRxnLUT::kz32], kcr_buf.data[CollisionalRxnLUT::kz33], kcr_buf.data[CollisionalRxnLUT::kz34], - kcr_buf.data[CollisionalRxnLUT::kz35], kcr_buf.data[CollisionalRxnLUT::kz36], kcr_buf.data[CollisionalRxnLUT::kz37], kcr_buf.data[CollisionalRxnLUT::kz38], kcr_buf.data[CollisionalRxnLUT::kz39], - kcr_buf.data[CollisionalRxnLUT::kz40], kcr_buf.data[CollisionalRxnLUT::kz41], kcr_buf.data[CollisionalRxnLUT::kz42], kcr_buf.data[CollisionalRxnLUT::kz43], kcr_buf.data[CollisionalRxnLUT::kz44], - kcr_buf.data[CollisionalRxnLUT::kz45], kcr_buf.data[CollisionalRxnLUT::kz46], kcr_buf.data[CollisionalRxnLUT::kz47], kcr_buf.data[CollisionalRxnLUT::kz48], kcr_buf.data[CollisionalRxnLUT::kz49], - kcr_buf.data[CollisionalRxnLUT::kz50], kcr_buf.data[CollisionalRxnLUT::kz51], kcr_buf.data[CollisionalRxnLUT::kz52], kcr_buf.data[CollisionalRxnLUT::kz53], kcr_buf.data[CollisionalRxnLUT::kz54], - species_tmpdens.data[SpLUT::DM], species_tmpdens.data[SpLUT::HDII], species_tmpdens.data[SpLUT::HeHII], - species_tmpdens.data[SpLUT::CI], species_tmpdens.data[SpLUT::CII], species_tmpdens.data[SpLUT::CO], species_tmpdens.data[SpLUT::CO2], - species_tmpdens.data[SpLUT::OI], species_tmpdens.data[SpLUT::OH], species_tmpdens.data[SpLUT::H2O], species_tmpdens.data[SpLUT::O2], - species_tmpdens.data[SpLUT::SiI], species_tmpdens.data[SpLUT::SiOI], species_tmpdens.data[SpLUT::SiO2I], - species_tmpdens.data[SpLUT::CH], species_tmpdens.data[SpLUT::CH2], species_tmpdens.data[SpLUT::COII], species_tmpdens.data[SpLUT::OII], - species_tmpdens.data[SpLUT::OHII], species_tmpdens.data[SpLUT::H2OII], species_tmpdens.data[SpLUT::H3OII], species_tmpdens.data[SpLUT::O2II], - species_tmpdens.data[SpLUT::Mg], species_tmpdens.data[SpLUT::Al], species_tmpdens.data[SpLUT::S], species_tmpdens.data[SpLUT::Fe], - species_tmpdens.data[SpLUT::SiM_dust], species_tmpdens.data[SpLUT::FeM_dust], species_tmpdens.data[SpLUT::Mg2SiO4_dust], species_tmpdens.data[SpLUT::MgSiO3_dust], species_tmpdens.data[SpLUT::Fe3O4_dust], - species_tmpdens.data[SpLUT::AC_dust], species_tmpdens.data[SpLUT::SiO2_dust], species_tmpdens.data[SpLUT::MgO_dust], species_tmpdens.data[SpLUT::FeS_dust], species_tmpdens.data[SpLUT::Al2O3_dust], - species_tmpdens.data[SpLUT::ref_org_dust], species_tmpdens.data[SpLUT::vol_org_dust], species_tmpdens.data[SpLUT::H2O_ice_dust], - grain_growth_rates.data[OnlyGrainSpLUT::SiM_dust], grain_growth_rates.data[OnlyGrainSpLUT::FeM_dust], grain_growth_rates.data[OnlyGrainSpLUT::Mg2SiO4_dust], grain_growth_rates.data[OnlyGrainSpLUT::MgSiO3_dust], grain_growth_rates.data[OnlyGrainSpLUT::Fe3O4_dust], - grain_growth_rates.data[OnlyGrainSpLUT::AC_dust], grain_growth_rates.data[OnlyGrainSpLUT::SiO2_dust], grain_growth_rates.data[OnlyGrainSpLUT::MgO_dust], grain_growth_rates.data[OnlyGrainSpLUT::FeS_dust], grain_growth_rates.data[OnlyGrainSpLUT::Al2O3_dust], - grain_growth_rates.data[OnlyGrainSpLUT::ref_org_dust], grain_growth_rates.data[OnlyGrainSpLUT::vol_org_dust], grain_growth_rates.data[OnlyGrainSpLUT::H2O_ice_dust], - &my_chemistry->radiative_transfer_HDI_dissociation, my_fields->RT_HDI_dissociation_rate, &my_chemistry->radiative_transfer_metal_ionization, my_fields->RT_CI_ionization_rate, my_fields->RT_OI_ionization_rate, - &my_chemistry->radiative_transfer_metal_dissociation, my_fields->RT_CO_dissociation_rate, my_fields->RT_OH_dissociation_rate, my_fields->RT_H2O_dissociation_rate - ); - -} - - } // namespace grackle::impl::fortran_wrapper #endif /* FORTRAN_FUNC_WRAPPERS_HPP */ diff --git a/src/clib/initialize_chemistry_data.cpp b/src/clib/initialize_chemistry_data.cpp index 294dbf7d6..8ead3cd27 100644 --- a/src/clib/initialize_chemistry_data.cpp +++ b/src/clib/initialize_chemistry_data.cpp @@ -463,16 +463,34 @@ extern "C" int local_initialize_chemistry_data(chemistry_data *my_chemistry, tm_info = localtime(&timer); strftime(tstr, 26, "%Y-%m-%d %H:%M:%S", tm_info); - FILE *fptr = fopen("GRACKLE_INFO", "w"); - fprintf(fptr, "%s\n", tstr); - show_version(fptr); - fprintf(fptr, "Grackle build options:\n"); - auto_show_config(fptr); - fprintf(fptr, "Grackle build flags:\n"); - auto_show_flags(fptr); - fprintf(fptr, "Grackle run-time parameters:\n"); - show_parameters(fptr, my_chemistry); - fclose(fptr); + const char* fname = "GRACKLE_INFO"; + + FILE *fptr = fopen(fname, "w"); + if (fptr == nullptr) { + // on posix platforms we could give a more detailed error message + // - posix (not the C or C++ standard) mandates that `fopen` set errno + // upon failure + // - properly informing the user of the error gets a little "hairy." We + // should use POSIX's strerror_r to get the error (since regular + // strerror may not be threadsafe). But we would need to account for + // the fact that glibc describes an incompatible signature for + // strerror_r (https://www.club.cc.cmu.edu/~cmccabe/blog_strerror.html) + // Since this isn't central to Grackle's functionality, we'll punt on + // this... + fprintf(stderr, "Failed to open \"%s\" file (to record config info)\n", + fname); + // an argument could be made that we should return with an error + } else { + fprintf(fptr, "%s\n", tstr); + show_version(fptr); + fprintf(fptr, "Grackle build options:\n"); + auto_show_config(fptr); + fprintf(fptr, "Grackle build flags:\n"); + auto_show_flags(fptr); + fprintf(fptr, "Grackle run-time parameters:\n"); + show_parameters(fptr, my_chemistry); + fclose(fptr); + } fprintf(stdout, "Grackle run-time parameters:\n"); show_parameters(stdout, my_chemistry); diff --git a/src/clib/initialize_rates.cpp b/src/clib/initialize_rates.cpp index 69091a05b..0fee8cec1 100644 --- a/src/clib/initialize_rates.cpp +++ b/src/clib/initialize_rates.cpp @@ -365,6 +365,7 @@ int init_kcol_rate_tables( if (n_kcol_rate_indices > (unsigned int)CollisionalRxnLUT::NUM_ENTRIES) { GRIMPL_ERROR("something went very wrong when counting up rates"); } else if (n_kcol_rate_indices == 0) { + delete tables; return GrPrintAndReturnErr("this logic shouldn't be executed in a config " "with 0 \"ordinary\" collisional rxn rates"); } diff --git a/src/clib/solve_rate_cool_g-cpp.cpp b/src/clib/solve_rate_cool_g-cpp.cpp index c7305d6aa..bde08114e 100644 --- a/src/clib/solve_rate_cool_g-cpp.cpp +++ b/src/clib/solve_rate_cool_g-cpp.cpp @@ -791,7 +791,7 @@ int solve_rate_cool_g( // declare 2 variables (primarily used for subcycling, but also used in // error reporting) int iter; - double ttmin; + double ttmin = huge8; // ------------------ Loop over subcycles ---------------- @@ -900,12 +900,12 @@ int solve_rate_cool_g( // Solve rate equations with one linearly implicit Gauss-Seidel // sweep of a backward Euler method (for all cells specified by // itmask_gs) - f_wrap::step_rate_g( + grackle::impl::step_rate_gauss_seidel( dtit.data(), idx_range, anydust, spsolvbuf.h2dust, rhoH.data(), spsolvbuf.dedot_prev, spsolvbuf.HIdot_prev, spsolvbuf.itmask_gs, - itmask_metal.data(), imetal, my_chemistry, my_fields, - *my_uvb_rates, spsolvbuf.grain_growth_rates, - spsolvbuf.species_tmpdens, spsolvbuf.kcr_buf, spsolvbuf.kshield_buf + itmask_metal.data(), my_chemistry, my_fields, *my_uvb_rates, + spsolvbuf.grain_growth_rates, spsolvbuf.species_tmpdens, + spsolvbuf.kcr_buf, spsolvbuf.kshield_buf ); // Solve rate equations with one linearly implicit Gauss-Seidel diff --git a/src/clib/solve_rate_cool_g.F b/src/clib/solve_rate_cool_g.F index 4386f2871..d8a6cc630 100644 --- a/src/clib/solve_rate_cool_g.F +++ b/src/clib/solve_rate_cool_g.F @@ -536,1536 +536,3 @@ subroutine rate_timestep_g( return end - -! ----------------------------------------------------------- -! This routine uses one linearly implicit Gauss-Seidel sweep of -! a backward-Euler time integrator to advance the rate equations -! by one (sub-)cycle (dtit). - - subroutine step_rate_g(de, HI, HII, HeI, HeII, HeIII, d, - & HM, H2I, H2II, DI, DII, HDI, dtit, - & in, jn, kn, is, ie, j, k, ispecies, anydust, - & k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, k11, - & k12, k13, k14, k15, k16, k17, k18, k19, k22, - & k24, k25, k26, k27, k28, k29, k30, - & k50, k51, k52, k53, k54, k55, k56, k57, k58, - & h2dust, rhoH, - & k24shield, k25shield, k26shield, - & k28shield, k29shield, k30shield, k31shield, - & HIp, HIIp, HeIp, HeIIp, HeIIIp, dep, - & HMp, H2Ip, H2IIp, DIp, DIIp, HDIp, - & dedot_prev, HIdot_prev, - & iradtrans, irt_honly, - & kphHI, kphHeI, kphHeII, - & itmask, itmask_metal - & , DM, HDII, HeHII, imetal, metal - & , imchem, idspecies, igrgr, idsub - & , CI, CII, CO, CO2 - & , OI, OH, H2O, O2 - & , SiI, SiOI, SiO2I - & , CH, CH2, COII, OII - & , OHII, H2OII, H3OII, O2II - & , Mg, Al, S, Fe - & , SiM, FeM, Mg2SiO4, MgSiO3, Fe3O4 - & , AC, SiO2D, MgO, FeS, Al2O3 - & , reforg, volorg, H2Oice - & , k125, k129, k130, k131, k132 - & , k133, k134, k135, k136, k137 - & , k148, k149, k150, k151, k152 - & , k153 - & , kz15 , kz16 , kz17 , kz18 , kz19 - & , kz20 , kz21 , kz22 , kz23 , kz24 - & , kz25 , kz26 , kz27 , kz28 , kz29 - & , kz30 , kz31 , kz32 , kz33 , kz34 - & , kz35 , kz36 , kz37 , kz38 , kz39 - & , kz40 , kz41 , kz42 , kz43 , kz44 - & , kz45 , kz46 , kz47 , kz48 , kz49 - & , kz50 , kz51 , kz52 , kz53 , kz54 - & , DMp, HDIIp, HeHIIp - & , CIp, CIIp, COp, CO2p - & , OIp, OHp, H2Op, O2p - & , SiIp, SiOIp, SiO2Ip - & , CHp, CH2p, COIIp, OIIp - & , OHIIp, H2OIIp, H3OIIp, O2IIp - & , Mgp, Alp, Sp, Fep - & , SiMp, FeMp, Mg2SiO4p, MgSiO3p, Fe3O4p - & , ACp, SiO2Dp, MgOp, FeSp, Al2O3p - & , reforgp, volorgp, H2Oicep - & , kdSiM, kdFeM, kdMg2SiO4, kdMgSiO3, kdFe3O4 - & , kdAC, kdSiO2D, kdMgO, kdFeS, kdAl2O3 - & , kdreforg, kdvolorg, kdH2Oice - & , idissHDI, kdissHDI, iionZ, kphCI, kphOI - & , idissZ, kdissCO, kdissOH, kdissH2O - & ) -c ------------------------------------------------------------------- - - implicit NONE -#include "grackle_fortran_types.def" - -! arguments - - integer ispecies, in, jn, kn, is, ie, j, k, - & iradtrans, irt_honly - real*8 dtit(in), dedot_prev(in), HIdot_prev(in) - MASK_TYPE itmask(in), itmask_metal(in), anydust - -! Density fields - - R_PREC de(in,jn,kn), HI(in,jn,kn), HII(in,jn,kn), - & HeI(in,jn,kn), HeII(in,jn,kn), HeIII(in,jn,kn), - & d(in,jn,kn), - & HM(in,jn,kn), H2I(in,jn,kn), H2II(in,jn,kn), - & DI(in,jn,kn), DII(in,jn,kn), HDI(in,jn,kn) - integer imetal, imchem, igrgr, idspecies, idsub - R_PREC metal(in,jn,kn) - R_PREC DM(in,jn,kn) , HDII(in,jn,kn) , HeHII(in,jn,kn) - & , CI(in,jn,kn) , CII(in,jn,kn) , CO(in,jn,kn) - & , CO2(in,jn,kn) , OI(in,jn,kn) , OH(in,jn,kn) - & , H2O(in,jn,kn) , O2(in,jn,kn) , SiI(in,jn,kn) - & , SiOI(in,jn,kn) , SiO2I(in,jn,kn) , CH(in,jn,kn) - & , CH2(in,jn,kn) , COII(in,jn,kn) , OII(in,jn,kn) - & , OHII(in,jn,kn) , H2OII(in,jn,kn) , H3OII(in,jn,kn) - & , O2II(in,jn,kn) , Mg(in,jn,kn) , Al(in,jn,kn) - & , S(in,jn,kn) , Fe(in,jn,kn) - R_PREC SiM(in,jn,kn), FeM(in,jn,kn), Mg2SiO4(in,jn,kn) - & , MgSiO3(in,jn,kn), Fe3O4(in,jn,kn), AC(in,jn,kn) - & , SiO2D(in,jn,kn), MgO(in,jn,kn), FeS(in,jn,kn) - & , Al2O3(in,jn,kn) - & , reforg(in,jn,kn), volorg(in,jn,kn), H2Oice(in,jn,kn) - -! Radiation Fields - R_PREC kphHI(in,jn,kn), kphHeI(in,jn,kn), kphHeII(in,jn,kn) - - integer idissHDI, iionZ, idissZ - R_PREC kdissHDI(in,jn,kn), kphCI(in,jn,kn), kphOI(in,jn,kn), - & kdissCO(in,jn,kn), kdissOH(in,jn,kn), kdissH2O(in,jn,kn) - - -! Rate values - - real*8 k1 (in), k2 (in), k3 (in), k4 (in), k5 (in), - & k6 (in), k7 (in), k8 (in), k9 (in), k10(in), - & k11(in), k12(in), k13(in), k14(in), k15(in), - & k16(in), k17(in), k18(in), k19(in), k22(in), - & k50(in), k51(in), k52(in), k53(in), k54(in), - & k55(in), k56(in), k57(in), k58(in), - & h2dust(in), rhoH(in), - & k24shield(in), k25shield(in), k26shield(in), - & k28shield(in), k29shield(in), k30shield(in), - & k31shield(in), - & k24, k25, k26, k27, k28, k29, k30 - real*8 k125(in), k129(in), k130(in), k131(in), k132(in) - & , k133(in), k134(in), k135(in), k136(in), k137(in) - & , k148(in), k149(in), k150(in), k151(in), k152(in) - & , k153(in) - & , kz15(in), kz16(in), kz17(in), kz18(in), kz19(in) - & , kz20(in), kz21(in), kz22(in), kz23(in), kz24(in) - & , kz25(in), kz26(in), kz27(in), kz28(in), kz29(in) - & , kz30(in), kz31(in), kz32(in), kz33(in), kz34(in) - & , kz35(in), kz36(in), kz37(in), kz38(in), kz39(in) - & , kz40(in), kz41(in), kz42(in), kz43(in), kz44(in) - & , kz45(in), kz46(in), kz47(in), kz48(in), kz49(in) - & , kz50(in), kz51(in), kz52(in), kz53(in), kz54(in) - -! temporaries (passed in) - - real*8 HIp(in), HIIp(in), HeIp(in), HeIIp(in), HeIIIp(in), - & HMp(in), H2Ip(in), H2IIp(in), dep(in), - & DIp(in), DIIp(in), HDIp(in) - real*8 DMp(in) , HDIIp(in) , HeHIIp(in) - & , CIp(in) , CIIp(in) , COp(in) - & , CO2p(in) , OIp(in) , OHp(in) - & , H2Op(in) , O2p(in) , SiIp(in) - & , SiOIp(in) , SiO2Ip(in) , CHp(in) - & , CH2p(in) , COIIp(in) , OIIp(in) - & , OHIIp(in) , H2OIIp(in) , H3OIIp(in) - & , O2IIp(in) , Mgp(in) , Alp(in) - & , Sp(in) , Fep(in) - R_PREC SiMp(in), FeMp(in), Mg2SiO4p(in) - & , MgSiO3p(in), Fe3O4p(in), ACp(in) - & , SiO2Dp(in), MgOp(in), FeSp(in) - & , Al2O3p(in) - & , reforgp(in), volorgp(in), H2Oicep(in) - real*8 kdSiM(in), kdFeM(in), kdMg2SiO4(in) - & , kdMgSiO3(in), kdFe3O4(in), kdAC(in) - & , kdSiO2D(in), kdMgO(in), kdFeS(in) - & , kdAl2O3(in) - & , kdreforg(in), kdvolorg(in), kdH2Oice(in) - -! locals - - integer i - real*8 scoef, acoef - -! A) the 6-species integrator -! - if (ispecies .eq. 1) then - - do i = is+1, ie+1 - if (itmask(i) .ne. MASK_FALSE) then - -! 1) HI - - scoef = k2(i)*HII(i,j,k)*de(i,j,k) - acoef = k1(i)*de(i,j,k) - & + k57(i)*HI(i,j,k) - & + k58(i)*HeI(i,j,k)/4._DKIND - & + k24shield(i) - if (iradtrans .eq. 1) acoef = acoef + kphHI(i,j,k) - HIp(i) = (scoef*dtit(i) + HI(i,j,k))/ - & (1._DKIND + acoef*dtit(i)) - if (HIp(i) .ne. HIp(i)) then -#ifdef _OPENMP -!$omp critical -#endif - write(*,*) 'HUGE HIp! :: ', i, j, k, HIp(i), HI(i,j,k), - $ HII(i,j,k), de(i,j,k), kphHI(i,j,k), - $ scoef, acoef, dtit(i) -#ifdef _OPENMP -!$omp end critical -#endif -c ERROR_MESSAGE - endif - -! 2) HII -c - scoef = k1(i)*HIp(i)*de(i,j,k) - & + k57(i)*HIp(i)*HIp(i) - & + k58(i)*HIp(i)*HeI(i,j,k)/4._DKIND - & + k24shield(i)*HIp(i) - if (iradtrans .eq. 1) - & scoef = scoef + kphHI(i,j,k)*HIp(i) - acoef = k2(i)*de (i,j,k) - HIIp(i) = (scoef*dtit(i) + HII(i,j,k))/ - & (1._DKIND +acoef*dtit(i)) -! - if (HIIp(i) .le. 0._DKIND) then !##### -#ifdef _OPENMP -!$omp critical -#endif - write(*,*) 'negative HIIp! :: ', i, j, k, HIIp(i), - $ scoef, dtit(i), HII(i,j,k), acoef, - $ k2(i), de(i,j,k), - $ kphHI(i,j,k), HIp(i), - $ k24shield(i) -#ifdef _OPENMP -!$omp end critical -#endif - endif - -! 3) Electron density - - scoef = 0._DKIND - & + k57(i)*HIp(i)*HIp(i) - & + k58(i)*HIp(i)*HeI(i,j,k)/4._DKIND - & + k24shield(i)*HI(i,j,k) - & + k25shield(i)*HeII(i,j,k)/4._DKIND - & + k26shield(i)*HeI(i,j,k)/4._DKIND - - if ( (iradtrans .eq. 1) .and. ( irt_honly .eq. 0) ) - & scoef = scoef + kphHI(i,j,k) * HI(i,j,k) - & + kphHeI(i,j,k) * HeI(i,j,k) / 4._DKIND - & + kphHeII(i,j,k) * HeII(i,j,k) / 4._DKIND - if ( (iradtrans .eq. 1) .and. ( irt_honly .eq. 1) ) - & scoef = scoef + kphHI(i,j,k) * HI(i,j,k) - - - - acoef = -(k1(i)*HI(i,j,k) - k2(i)*HII(i,j,k) - & + k3(i)*HeI(i,j,k)/4._DKIND - - & k6(i)*HeIII(i,j,k)/4._DKIND - & + k5(i)*HeII(i,j,k)/4._DKIND - - & k4(i)*HeII(i,j,k)/4._DKIND) - dep(i) = (scoef*dtit(i) + de(i,j,k)) - & / (1._DKIND + acoef*dtit(i)) - - endif ! itmask - enddo - - endif ! (ispecies .eq. 1) - -! --- (B) Do helium chemistry in any case: (for all ispecies values) --- - - do i = is+1, ie+1 - if (itmask(i) .ne. MASK_FALSE) then - -! 4) HeI - - scoef = k4(i)*HeII(i,j,k)*de(i,j,k) - acoef = k3(i)*de(i,j,k) - & + k26shield(i) - - if ( (iradtrans .eq. 1) .and. (irt_honly .eq. 0)) - & acoef = acoef + kphHeI(i,j,k) - if (ispecies .gt. 3) then - scoef = scoef + 4._DKIND * ( 0._DKIND - & + k152(i) * HeHII(i,j,k) * HI(i,j,k) / 5._DKIND - & + k153(i) * HeHII(i,j,k) * de(i,j,k) / 5._DKIND - & ) - acoef = acoef - & + k148(i) * HII(i,j,k) - & + k149(i) * HII(i,j,k) - & + k150(i) * H2II(i,j,k) / 2._DKIND - endif - HeIp(i) = ( scoef*dtit(i) + HeI(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - -! 5) HeII - - scoef = k3(i)*HeIp(i)*de(i,j,k) - & + k6(i)*HeIII(i,j,k)*de(i,j,k) - & + k26shield(i)*HeIp(i) - - if ( (iradtrans .eq. 1) .and. (irt_honly .eq. 0)) - & scoef = scoef + kphHeI(i,j,k)*HeIp(i) - - acoef = k4(i)*de(i,j,k) + k5(i)*de(i,j,k) - & + k25shield(i) - - if ( (iradtrans .eq. 1) .and. (irt_honly .eq. 0)) - & acoef = acoef + kphHeII(i,j,k) - if (ispecies .gt. 3) then - acoef = acoef - & + k151(i) * HI(i,j,k) - endif - HeIIp(i) = ( scoef*dtit(i) + HeII(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - -! 6) HeIII - - scoef = k5(i)*HeIIp(i)*de(i,j,k) - & + k25shield(i)*HeIIp(i) - if ((iradtrans .eq. 1) .and. (irt_honly .eq. 0)) - & scoef = scoef + kphHeII(i,j,k) * HeIIp(i) - acoef = k6(i)*de(i,j,k) - HeIIIp(i) = ( scoef*dtit(i) + HeIII(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - endif ! itmask - enddo - -c --- (C) Now do extra 3-species for molecular hydrogen --- - - if (ispecies .gt. 1) then - -! First, do HI/HII with molecular hydrogen terms - - do i = is+1, ie+1 - if (itmask(i) .ne. MASK_FALSE) then - -! 1) HI -! - scoef = k2(i) * HII(i,j,k) * de(i,j,k) - & + 2._DKIND*k13(i)* HI(i,j,k) * H2I(i,j,k)/2._DKIND - & + k11(i)* HII(i,j,k) * H2I(i,j,k)/2._DKIND - & + 2._DKIND*k12(i)* de(i,j,k) * H2I(i,j,k)/2._DKIND - & + k14(i)* HM(i,j,k) * de(i,j,k) - & + k15(i)* HM(i,j,k) * HI(i,j,k) - & + 2._DKIND*k16(i)* HM(i,j,k) * HII(i,j,k) - & + 2._DKIND*k18(i)* H2II(i,j,k)* de(i,j,k)/2._DKIND - & + k19(i)* H2II(i,j,k)* HM(i,j,k)/2._DKIND - & + 2._DKIND*k31shield(i) * H2I(i,j,k)/2._DKIND - - acoef = k1(i) * de(i,j,k) - & + k7(i) * de(i,j,k) - & + k8(i) * HM(i,j,k) - & + k9(i) * HII(i,j,k) - & + k10(i)* H2II(i,j,k)/2._DKIND - & + 2._DKIND*k22(i)* HI(i,j,k)**2 - & + k57(i)* HI(i,j,k) - & + k58(i)* HeI(i,j,k)/4._DKIND - & + k24shield(i) - - if (iradtrans .eq. 1) acoef = acoef + kphHI(i,j,k) - if (iradtrans .eq. 1) then - if ((ispecies .gt. 2).and.(idissHDI .gt. 0)) then - scoef = scoef - & + kdissHDI(i,j,k) * HDI(i,j,k)/3.0_DKIND - endif - if ( (imchem .eq. 1) .and. - & (itmask_metal(i) .ne. MASK_FALSE) ) then - if (idissZ .gt. 0) then - scoef = scoef - & + kdissOH (i,j,k) * OH(i,j,k) /17.0_DKIND - & + kdissH2O(i,j,k) * H2O(i,j,k)/18.0_DKIND - endif - endif - endif - - if (anydust .ne. MASK_FALSE) then - if(itmask_metal(i) .ne. MASK_FALSE) then - acoef = acoef + 2._DKIND * h2dust(i) * rhoH(i) - endif - endif -#ifdef CONTRIBUTION_OF_MINOR_SPECIES - if (ispecies .gt. 2) then - scoef = scoef - & + k50(i) * HII(i,j,k) * DI(i,j,k) / 2._DKIND - & + k54(i) * H2I(i,j,k) * DI(i,j,k) / 4._DKIND - acoef = acoef - & + k51(i) * DII(i,j,k) / 2._DKIND - & + k55(i) * HDI(i,j,k) / 3._DKIND - endif -#endif - if (ispecies .gt. 3) then - scoef = scoef - & + k131(i) * HDII(i,j,k) * de(i,j,k) / 3._DKIND - & + k134(i) * HII(i,j,k) * DM(i,j,k) / 2._DKIND - & + k135(i) * HM(i,j,k) * DI(i,j,k) / 2._DKIND - & + k150(i) * HeI(i,j,k) * H2II(i,j,k) / 8._DKIND - & + k153(i) * HeHII(i,j,k) * de(i,j,k) / 5._DKIND - acoef = acoef - & + k125(i) * HDII(i,j,k) / 3._DKIND - & + k130(i) * DII(i,j,k) / 2._DKIND - & + k136(i) * DM(i,j,k) / 2._DKIND - & + k137(i) * DM(i,j,k) / 2._DKIND - & + k151(i) * HeII(i,j,k) / 4._DKIND - & + k152(i) * HeHII(i,j,k) / 5._DKIND - endif - - if ( (imchem .eq. 1) .and. - & (itmask_metal(i) .ne. MASK_FALSE) ) then - scoef = scoef - & + kz20(i) * CI(i,j,k) * H2I(i,j,k) / 24._DKIND - & + kz21(i) * OI(i,j,k) * H2I(i,j,k) / 32._DKIND - & + kz22(i) * HII(i,j,k) * OI(i,j,k) / 16._DKIND - & + kz23(i) * H2I(i,j,k) * CH(i,j,k) / 26._DKIND - & + kz24(i) * H2I(i,j,k) * OH(i,j,k) / 34._DKIND - & + kz26(i) * OH(i,j,k) * CO(i,j,k) / 476._DKIND - & + kz28(i) * CI(i,j,k) * OH(i,j,k) / 204._DKIND - & + kz32(i) * OI(i,j,k) * CH(i,j,k) / 208._DKIND - & + kz33(i) * OI(i,j,k) * OH(i,j,k) / 272._DKIND - & + kz34(i) * HII(i,j,k) * OH(i,j,k) / 17._DKIND - & + kz35(i) * HII(i,j,k) * H2O(i,j,k) / 18._DKIND - & + kz36(i) * HII(i,j,k) * O2(i,j,k) / 32._DKIND - & + kz37(i) * CII(i,j,k) * OH(i,j,k) / 204._DKIND - & + kz40(i) * OII(i,j,k) * H2I(i,j,k) / 32._DKIND - & + kz41(i) * OHII(i,j,k) * H2I(i,j,k) / 34._DKIND - & + kz42(i) * H2OII(i,j,k) * H2I(i,j,k) / 36._DKIND - & + kz46(i) * H2OII(i,j,k) * de(i,j,k) / 18._DKIND - & + kz48(i) * H3OII(i,j,k) * de(i,j,k) / 19._DKIND - & + kz49(i) * H3OII(i,j,k) * de(i,j,k) / 9.5_DKIND - & + kz52(i) * SiI(i,j,k) * OH(i,j,k) / 476._DKIND - & + kz54(i) * SiOI(i,j,k) * OH(i,j,k) / 748._DKIND - acoef = acoef - & + kz15(i) * CH(i,j,k) / 13._DKIND - & + kz16(i) * CH2(i,j,k) / 14._DKIND - & + kz17(i) * OH(i,j,k) / 17._DKIND - & + kz18(i) * H2O(i,j,k) / 18._DKIND - & + kz19(i) * O2(i,j,k) / 32._DKIND - & + kz27(i) * CI(i,j,k) / 12._DKIND - & + kz30(i) * OI(i,j,k) / 16._DKIND - & + kz39(i) * OII(i,j,k) / 16._DKIND - & + kz43(i) * COII(i,j,k) / 28._DKIND - endif - HIp(i) = ( scoef*dtit(i) + HI(i,j,k) ) / - & ( 1. + acoef*dtit(i) ) - if (HIp(i) .ne. HIp(i)) then -#ifdef _OPENMP -!$omp critical -#endif - write(*,*) 'HUGE HIp! :: ', i, j, k, HIp(i), HI(i,j,k), - $ HII(i,j,k), de(i,j,k), H2I(i,j,k), - $ kphHI(i,j,k) -#ifdef _OPENMP -!$omp end critical -#endif - endif - -! 2) HII - - scoef = k1(i) * HI(i,j,k) * de(i,j,k) - & + k10(i) * H2II(i,j,k)*HI(i,j,k)/2._DKIND - & + k57(i) * HI(i,j,k) * HI(i,j,k) - & + k58(i) * HI(i,j,k) * HeI(i,j,k)/4._DKIND - & + k24shield(i)*HI(i,j,k) - - if (iradtrans .eq. 1) - & scoef = scoef + kphHI(i,j,k) * HI(i,j,k) - - acoef = k2(i) * de(i,j,k) - & + k9(i) * HI(i,j,k) - & + k11(i) * H2I(i,j,k)/2._DKIND - & + k16(i) * HM(i,j,k) - & + k17(i) * HM(i,j,k) -#ifdef CONTRIBUTION_OF_MINOR_SPECIES - if (ispecies .gt. 2) then - scoef = scoef - & + k51(i) * HI (i,j,k) * DII(i,j,k) / 2._DKIND - & + k52(i) * H2I(i,j,k) * DII(i,j,k) / 4._DKIND - acoef = acoef - & + k50(i) * DI (i,j,k) / 2._DKIND - & + k53(i) * HDI(i,j,k) / 3._DKIND - endif -#endif - if (ispecies .gt. 3) then - scoef = scoef - & + k125(i) * HDII(i,j,k) * HI(i,j,k) / 3._DKIND - acoef = acoef - & + k129(i) * DI(i,j,k) / 2._DKIND - & + k134(i) * DM(i,j,k) / 2._DKIND - & + k148(i) * HeI(i,j,k) / 4._DKIND - & + k149(i) * HeI(i,j,k) / 4._DKIND - endif - - if ( (imchem .eq. 1) .and. - & (itmask_metal(i) .ne. MASK_FALSE) ) then - scoef = scoef - & + kz39(i) * OII(i,j,k) * HI(i,j,k) / 16._DKIND - & + kz43(i) * COII(i,j,k) * HI(i,j,k) / 28._DKIND - acoef = acoef - & + kz22(i) * OI(i,j,k) / 16._DKIND - & + kz34(i) * OH(i,j,k) / 17._DKIND - & + kz35(i) * H2O(i,j,k) / 18._DKIND - & + kz36(i) * O2(i,j,k) / 32._DKIND - endif - HIIp(i) = ( scoef*dtit(i) + HII(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) -! -! 3) electrons: - - scoef = k8(i) * HM(i,j,k) * HI(i,j,k) - & + k15(i)* HM(i,j,k) * HI(i,j,k) - & + k17(i)* HM(i,j,k) * HII(i,j,k) - & + k57(i)* HI(i,j,k) * HI(i,j,k) - & + k58(i)* HI(i,j,k) * HeI(i,j,k)/4._DKIND -! - & + k24shield(i)*HIp(i) - & + k25shield(i)*HeIIp(i)/4._DKIND - & + k26shield(i)*HeIp(i)/4._DKIND - - if ( (iradtrans .eq. 1) .and. (irt_honly .eq. 0) ) - & scoef = scoef + kphHI(i,j,k) * HIp(i) - & + kphHeI(i,j,k) * HeIp(i) / 4._DKIND - & + kphHeII(i,j,k) * HeIIp(i) / 4._DKIND - if ( (iradtrans .eq. 1) .and. (irt_honly .eq. 1) ) - & scoef = scoef + kphHI(i,j,k) * HIp(i) - if (iradtrans .eq. 1) then - if ( (imchem .eq. 1) .and. - & (itmask_metal(i) .ne. MASK_FALSE) ) then - if (iionZ .gt. 0) then - scoef = scoef - & + kphCI(i,j,k) * CI(i,j,k)/12.0_DKIND - & + kphOI(i,j,k) * OI(i,j,k)/16.0_DKIND - endif - endif - endif - - acoef = - (k1(i) *HI(i,j,k) - k2(i)*HII(i,j,k) - & + k3(i) *HeI(i,j,k)/4._DKIND - - & k6(i)*HeIII(i,j,k)/4._DKIND - & + k5(i) *HeII(i,j,k)/4._DKIND - - & k4(i)*HeII(i,j,k)/4._DKIND - & + k14(i)*HM(i,j,k) - & - k7(i) *HI(i,j,k) - & - k18(i)*H2II(i,j,k)/2._DKIND) -#ifdef CONTRIBUTION_OF_MINOR_SPECIES - if (ispecies .gt. 2) then - scoef = scoef - & + k56(i) * DI (i,j,k) * HM(i,j,k) / 2._DKIND - acoef = acoef - & - k1 (i) * DI (i,j,k) / 2._DKIND - & + k2 (i) * DII(i,j,k) / 2._DKIND - endif -#endif - if (ispecies .gt. 3) then - scoef = scoef - & + k137(i) * DM(i,j,k) * HI(i,j,k) / 2._DKIND - acoef = acoef - & + k131(i) * HDII(i,j,k) / 3._DKIND - & + k132(i) * DI(i,j,k) / 2._DKIND - & + k153(i) * HeHII(i,j,k) / 5._DKIND - endif - - if ( (imchem .eq. 1) .and. - & (itmask_metal(i) .ne. MASK_FALSE) ) then - scoef = scoef - acoef = acoef - & + kz44(i) * CII(i,j,k) / 12._DKIND - & + kz45(i) * OII(i,j,k) / 16._DKIND - & + kz46(i) * H2OII(i,j,k) / 18._DKIND - & + kz47(i) * H2OII(i,j,k) / 18._DKIND - & + kz48(i) * H3OII(i,j,k) / 19._DKIND - & + kz49(i) * H3OII(i,j,k) / 19._DKIND - & + kz50(i) * O2II(i,j,k) / 32._DKIND - endif - dep(i) = ( scoef*dtit(i) + de(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - -! 7) H2 - - scoef = 2._DKIND*(k8(i) * HM(i,j,k) * HI(i,j,k) - & + k10(i) * H2II(i,j,k) * HI(i,j,k)/2._DKIND - & + k19(i) * H2II(i,j,k) * HM(i,j,k)/2._DKIND - & + k22(i) * HI(i,j,k) * (HI(i,j,k))**2._DKIND) - acoef = ( k13(i)*HI(i,j,k) + k11(i)*HII(i,j,k) - & + k12(i)*de(i,j,k) ) - & + k29shield(i) + k31shield(i) - - if (anydust .ne. MASK_FALSE) then - if(itmask_metal(i) .ne. MASK_FALSE) then - scoef = scoef + 2._DKIND * h2dust(i) * - & HI(i,j,k) * rhoH(i) - endif - endif -#ifdef CONTRIBUTION_OF_MINOR_SPECIES - if (ispecies .gt. 2) then - scoef = scoef + 2._DKIND * ( - & k53(i) * HDI(i,j,k) * HII(i,j,k) / 3._DKIND - & + k55(i) * HDI(i,j,k) * HI (i,j,k) / 3._DKIND - & ) - acoef = acoef - & + k52(i) * DII(i,j,k) / 2._DKIND - & + k54(i) * DI (i,j,k) / 2._DKIND - endif -#endif - if ( (imchem .eq. 1) .and. - & (itmask_metal(i) .ne. MASK_FALSE) ) then - scoef = scoef + 2._DKIND * ( 0._DKIND - & + kz15(i) * HI(i,j,k) * CH(i,j,k) / 13._DKIND - & + kz16(i) * HI(i,j,k) * CH2(i,j,k) / 14._DKIND - & + kz17(i) * HI(i,j,k) * OH(i,j,k) / 17._DKIND - & + kz18(i) * HI(i,j,k) * H2O(i,j,k) / 18._DKIND - & + kz47(i) * H2OII(i,j,k) * de(i,j,k) / 18._DKIND - & ) - acoef = acoef - & + kz20(i) * CI(i,j,k) / 12._DKIND - & + kz21(i) * OI(i,j,k) / 16._DKIND - & + kz23(i) * CH(i,j,k) / 13._DKIND - & + kz24(i) * OH(i,j,k) / 17._DKIND - & + kz40(i) * OII(i,j,k) / 16._DKIND - & + kz41(i) * OHII(i,j,k) / 17._DKIND - & + kz42(i) * H2OII(i,j,k) / 18._DKIND - & + kz51(i) * CI(i,j,k) / 12._DKIND - if ((igrgr .eq. 1) .or. (idsub .eq. 1)) then - if (idspecies .gt. 0) then - scoef = scoef + 2._DKIND * - & kdMgSiO3 (i) * 2._DKIND - - endif - if (idspecies .gt. 1) then - scoef = scoef + 2._DKIND * ( - & kdMg2SiO4 (i) * 3._DKIND - & + kdFe3O4 (i) * 4._DKIND - & + kdMgO (i) - & + kdAl2O3 (i) * 3._DKIND - & ) - endif - if (idspecies .gt. 2) then - acoef = acoef - & + kdvolorg (i) / H2I(i,j,k) * 2._DKIND * 2._DKIND - endif - endif - endif - H2Ip(i) = ( scoef*dtit(i) + H2I(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - -! 8) H- - - scoef = k7(i) * HI(i,j,k) * de(i,j,k) - acoef = (k8(i) + k15(i)) * HI(i,j,k) + - & (k16(i) + k17(i)) * HII(i,j,k) + - & k14(i) * de(i,j,k) + k19(i) * H2II(i,j,k)/2.0 + - & k27 -#ifdef CONTRIBUTION_OF_MINOR_SPECIES - if (ispecies .gt. 2) then - acoef = acoef - & + k56(i) * DI (i,j,k) / 2._DKIND - endif -#endif - if (ispecies .gt. 3) then - scoef = scoef - & + k136(i) * DM(i,j,k) * HI(i,j,k) / 2._DKIND - acoef = acoef - & + k135(i) * DI(i,j,k) / 2._DKIND - endif - HMp(i) = (scoef*dtit(i) + HM(i,j,k)) - & / (1.0 + acoef*dtit(i)) - - -! 9) H2+ - - H2IIp(i) = 2._DKIND*( k9 (i)*HIp(i)*HIIp(i) - & + k11(i)*H2Ip(i)/2._DKIND*HIIp(i) - & + k17(i)*HMp(i)*HIIp(i) - & + k29shield(i)*H2Ip(i) - & ) - & / ( k10(i)*HIp(i) + k18(i)*dep(i) - & + k19(i)*HMp(i) - & + (k28shield(i)+k30shield(i)) - & ) - if (ispecies .gt. 3) then - H2IIp(i) = 2._DKIND * ( k9 (i)*HIp(i)*HIIp(i) - & + k11(i)*H2Ip(i)/2._DKIND*HIIp(i) - & + k17(i)*HMp(i)*HIIp(i) - & + k29shield(i)*H2Ip(i) - & + k152(i)*HeHII(i,j,k)*HIp(i)/5._DKIND - & ) - & / ( k10(i)*HIp(i) + k18(i)*dep(i) - & + k19(i)*HMp(i) - & + (k28shield(i)+k30shield(i)) - & + k150(i)*HeIp(i)/4._DKIND - & ) - endif - endif ! itmask - enddo -! - endif ! H2 - -! --- (D) Now do extra 3-species for molecular HD --- -! - if (ispecies .gt. 2) then - do i = is+1, ie+1 - if (itmask(i) .ne. MASK_FALSE) then -! -! 1) DI -! - scoef = ( k2(i) * DII(i,j,k) * de(i,j,k) - & + k51(i)* DII(i,j,k) * HI(i,j,k) - & + 2._DKIND*k55(i)* HDI(i,j,k) * - & HI(i,j,k)/3._DKIND - & ) - acoef = k1(i) * de(i,j,k) - & + k50(i) * HII(i,j,k) - & + k54(i) * H2I(i,j,k)/2._DKIND - & + k56(i) * HM(i,j,k) - & + k24shield(i) - if (iradtrans .eq. 1) acoef = acoef + kphHI(i,j,k) - if (ispecies .gt. 3) then - scoef = scoef + 2._DKIND * ( 0._DKIND - & + k131(i) * HDII(i,j,k) * de(i,j,k) / 3._DKIND - & + k133(i) * DII(i,j,k) * DM(i,j,k) / 2._DKIND - & + k134(i) * HII(i,j,k) * DM(i,j,k) / 2._DKIND - & + k136(i) * DM(i,j,k) * HI(i,j,k) / 2._DKIND - & ) - acoef = acoef - & + k129(i) * HII(i,j,k) - & + k132(i) * de(i,j,k) - & + k135(i) * HM(i,j,k) - endif - if (iradtrans .eq. 1) then - if (idissHDI .gt. 0) then - scoef = scoef - & + 2._DKIND * kdissHDI(i,j,k) * HDI(i,j,k)/3.0_DKIND - endif - endif - DIp(i) = ( scoef*dtit(i) + DI(i,j,k) ) / - & ( 1._DKIND + acoef*dtit(i) ) - -! 2) DII -c - scoef = ( k1(i) * DI(i,j,k) * de(i,j,k) - & + k50(i) * HII(i,j,k)* DI(i,j,k) - & + 2._DKIND*k53(i) * HII(i,j,k)* HDI(i,j,k)/3._DKIND - & ) - & + k24shield(i)*DI(i,j,k) - acoef = 0._DKIND - !! initialize GC202002 - if (iradtrans .eq. 1) scoef = scoef + kphHI(i,j,k)*DI(i,j,k) - acoef = k2(i) * de(i,j,k) - & + k51(i) * HI(i,j,k) - & + k52(i) * H2I(i,j,k)/2._DKIND - if (ispecies .gt. 3) then - acoef = acoef - & + k130(i) * HI(i,j,k) - & + k133(i) * DM(i,j,k) / 2._DKIND - endif - DIIp(i) = ( scoef*dtit(i) + DII(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - -! 3) HDI -c - scoef = 3._DKIND*(k52(i) * DII(i,j,k)* - & H2I(i,j,k)/2._DKIND/2._DKIND - & + k54(i) * DI(i,j,k) * H2I(i,j,k)/2._DKIND/2._DKIND -!! & + 2._DKIND*k56(i) * DI(i,j,k) * HM(i,j,k)/2._DKIND - !! corrected by GC202005 - & + k56(i) * DI(i,j,k) * HM(i,j,k)/2._DKIND - & ) - acoef = k53(i) * HII(i,j,k) - & + k55(i) * HI(i,j,k) - if (iradtrans .eq. 1) then - if (idissHDI .gt. 0) then - acoef = acoef - & + kdissHDI(i,j,k) - endif - endif - if (ispecies .gt. 3) then - scoef = scoef + 3._DKIND * ( 0._DKIND - & + k125(i) * HDII(i,j,k) * HI(i,j,k) / 3._DKIND - & + k137(i) * DM(i,j,k) * HI(i,j,k) / 2._DKIND - & ) - endif - HDIp(i) = ( scoef*dtit(i) + HDI(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - endif ! itmask - enddo - endif - -! --- (D2) Now do extra 3-species for minor primordial species --- -! - if (ispecies .gt. 3) then - do i = is+1, ie+1 - if (itmask(i) .ne. MASK_FALSE) then -! -! 1) DM -! - scoef = - & k132(i) * DI(i,j,k) * de(i,j,k) - & + k135(i) * HM(i,j,k) * DI(i,j,k) - acoef = - & k133(i) * DII(i,j,k) / 2._DKIND - & + k134(i) * HII(i,j,k) - & + k136(i) * HI(i,j,k) - & + k137(i) * HI(i,j,k) - - DMp(i) = ( scoef*dtit(i) + DM(i,j,k) ) / - & ( 1._DKIND + acoef*dtit(i) ) - -! 2) HDII -c - scoef = 3._DKIND * ( - & k129(i) * DI(i,j,k) * HII(i,j,k) / 2._DKIND - & + k130(i) * DII(i,j,k) * HI(i,j,k) / 2._DKIND - & ) - acoef = - & k125(i) * HI(i,j,k) - & + k131(i) * de(i,j,k) - - HDIIp(i) = ( scoef*dtit(i) + HDII(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - -! 3) HeHII -c - scoef = 5._DKIND * ( - & k148(i) * HeI(i,j,k) * HII(i,j,k) / 4._DKIND - & + k149(i) * HeI(i,j,k) * HII(i,j,k) / 4._DKIND - & + k150(i) * HeI(i,j,k) * H2II(i,j,k) / 8._DKIND - & + k151(i) * HeII(i,j,k) * HI(i,j,k) / 4._DKIND - & ) - acoef = - & k152(i) * HI(i,j,k) - & + k153(i) * de(i,j,k) - - HeHIIp(i) = ( scoef*dtit(i) + HeHII(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - endif ! itmask - enddo - endif - -! --- (D3) Now do metal species --- -! - if (imchem .eq. 1) then - do i = is+1, ie+1 - if (itmask_metal(i) .ne. MASK_FALSE) then - -C***** CI ********** - scoef = 0._DKIND + 12._DKIND * ( 0._DKIND - & + kz15(i) * HI(i,j,k) * CH(i,j,k) / 13._DKIND - & + kz44(i) * CII(i,j,k) * de(i,j,k) / 12._DKIND - & ) - acoef = 0._DKIND - & + kz20(i) * H2I(i,j,k) / 2._DKIND - & + kz27(i) * HI(i,j,k) - & + kz28(i) * OH(i,j,k) / 17._DKIND - & + kz29(i) * O2(i,j,k) / 32._DKIND - & + kz51(i) * H2I(i,j,k) / 2._DKIND - if ( ( igrgr .eq. 1 ) .or. ( idsub .eq. 1) ) then - if (idspecies .gt. 0) then - acoef = acoef - & + kdAC (i) / CI(i,j,k) * 12._DKIND - endif - endif - if (iradtrans .eq. 1) then - if (iionZ .gt. 0) then - acoef = acoef - & + kphCI(i,j,k) - endif - if (idissZ .gt. 0) then - scoef = scoef + 12._DKIND * - & kdissCO (i,j,k) * CO(i,j,k) /28.0_DKIND - endif - endif - - CIp(i) = ( scoef*dtit(i) + CI(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - -C***** CII ********** - scoef = 0._DKIND + 12._DKIND * ( 0._DKIND - & ) - acoef = 0._DKIND - & + kz37(i) * OH(i,j,k) / 17._DKIND - & + kz38(i) * O2(i,j,k) / 32._DKIND - & + kz44(i) * de(i,j,k) - if (iradtrans .eq. 1) then - if (iionZ .gt. 0) then - scoef = scoef - & + kphCI(i,j,k) * CI(i,j,k) - endif - endif - - CIIp(i) = ( scoef*dtit(i) + CII(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - -C***** CO ********** - scoef = 0._DKIND + 28._DKIND * ( 0._DKIND - & + kz28(i) * CI(i,j,k) * OH(i,j,k) / 204._DKIND - & + kz29(i) * CI(i,j,k) * O2(i,j,k) / 384._DKIND - & + kz32(i) * OI(i,j,k) * CH(i,j,k) / 208._DKIND - & + kz38(i) * CII(i,j,k) * O2(i,j,k) / 384._DKIND - & + kz43(i) * COII(i,j,k) * HI(i,j,k) / 28._DKIND - & ) - acoef = 0._DKIND - & + kz26(i) * OH(i,j,k) / 17._DKIND - if ( ( igrgr .eq. 1 ) .or. ( idsub .eq. 1) ) then - if (idspecies .gt. 2) then - acoef = acoef - & + kdreforg (i) / CO(i,j,k) * 17._DKIND * 0.5_DKIND - & + kdvolorg (i) / CO(i,j,k) * 17._DKIND - endif - endif - if (iradtrans .eq. 1) then - if (idissZ .gt. 0) then - acoef = acoef - & + kdissCO (i,j,k) - endif - endif - - COp(i) = ( scoef*dtit(i) + CO(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - -C***** CO2 ********** - scoef = 0._DKIND + 44._DKIND * ( 0._DKIND - & + kz26(i) * OH(i,j,k) * CO(i,j,k) / 476._DKIND - & ) - acoef = 0._DKIND - - CO2p(i) = ( scoef*dtit(i) + CO2(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - -C***** OI ********** - scoef = 0._DKIND + 16._DKIND * ( 0._DKIND - & + kz17(i) * HI(i,j,k) * OH(i,j,k) / 17._DKIND - & + kz19(i) * HI(i,j,k) * O2(i,j,k) / 32._DKIND - & + kz25(i) * OH(i,j,k) * OH(i,j,k) / 289._DKIND - & + kz29(i) * CI(i,j,k) * O2(i,j,k) / 384._DKIND - & + kz39(i) * OII(i,j,k) * HI(i,j,k) / 16._DKIND - & + kz45(i) * OII(i,j,k) * de(i,j,k) / 16._DKIND - & + kz47(i) * H2OII(i,j,k) * de(i,j,k) / 18._DKIND - & + kz50(i) * O2II(i,j,k) * de(i,j,k) / 16._DKIND - & + kz53(i) * SiI(i,j,k) * O2(i,j,k) / 896._DKIND - & ) - acoef = 0._DKIND - & + kz21(i) * H2I(i,j,k) / 2._DKIND - & + kz22(i) * HII(i,j,k) - & + kz30(i) * HI(i,j,k) - & + kz31(i) * OI(i,j,k) / 8._DKIND - & + kz32(i) * CH(i,j,k) / 13._DKIND - & + kz33(i) * OH(i,j,k) / 17._DKIND - if (iradtrans .eq. 1) then - if (iionZ .gt. 0) then - acoef = acoef - & + kphOI(i,j,k) - endif - if (idissZ .gt. 0) then - scoef = scoef + 16._DKIND * - & ( kdissOH (i,j,k) * OH(i,j,k) /17.0_DKIND - & + kdissCO (i,j,k) * CO(i,j,k) /28.0_DKIND) - endif - endif - - OIp(i) = ( scoef*dtit(i) + OI(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - -C***** OH ********** - scoef = 0._DKIND + 17._DKIND * ( 0._DKIND - & + kz18(i) * HI(i,j,k) * H2O(i,j,k) / 18._DKIND - & + kz19(i) * HI(i,j,k) * O2(i,j,k) / 32._DKIND - & + kz21(i) * OI(i,j,k) * H2I(i,j,k) / 32._DKIND - & + kz30(i) * OI(i,j,k) * HI(i,j,k) / 16._DKIND - & + kz46(i) * H2OII(i,j,k) * de(i,j,k) / 18._DKIND - & + kz49(i) * H3OII(i,j,k) * de(i,j,k) / 19._DKIND - & ) - acoef = 0._DKIND - & + kz17(i) * HI(i,j,k) - & + kz24(i) * H2I(i,j,k) / 2._DKIND - & + kz25(i) * OH(i,j,k) / 8.5_DKIND - & + kz26(i) * CO(i,j,k) / 28._DKIND - & + kz28(i) * CI(i,j,k) / 12._DKIND - & + kz33(i) * OI(i,j,k) / 16._DKIND - & + kz34(i) * HII(i,j,k) - & + kz37(i) * CII(i,j,k) / 12._DKIND - & + kz52(i) * SiI(i,j,k) / 28._DKIND - & + kz54(i) * SiOI(i,j,k) / 44._DKIND - if (iradtrans .eq. 1) then - if (idissZ .gt. 0) then - acoef = acoef - & + kdissOH (i,j,k) - scoef = scoef + 17._DKIND * - & kdissH2O(i,j,k) * H2O(i,j,k)/18.0_DKIND - endif - endif - - OHp(i) = ( scoef*dtit(i) + OH(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - -C***** H2O ********** - scoef = 0._DKIND + 18._DKIND * ( 0._DKIND - & + kz24(i) * H2I(i,j,k) * OH(i,j,k) / 34._DKIND - & + kz25(i) * OH(i,j,k) * OH(i,j,k) / 289._DKIND - & + kz48(i) * H3OII(i,j,k) * de(i,j,k) / 19._DKIND - & ) - acoef = 0._DKIND - & + kz18(i) * HI(i,j,k) - & + kz35(i) * HII(i,j,k) - if ( ( igrgr .eq. 1 ) .or. ( idsub .eq. 1) ) then - if (idspecies .gt. 0) then - acoef = acoef - & + kdMgSiO3 (i) / H2O(i,j,k) * 18._DKIND * 2._DKIND - endif - if (idspecies .gt. 1) then - acoef = acoef - & + kdMg2SiO4 (i) / H2O(i,j,k) * 18._DKIND * 3._DKIND - & + kdFe3O4 (i) / H2O(i,j,k) * 18._DKIND * 4._DKIND - & + kdMgO (i) / H2O(i,j,k) * 18._DKIND - & + kdAl2O3 (i) / H2O(i,j,k) * 18._DKIND * 3._DKIND - endif - if (idspecies .gt. 2) then - acoef = acoef - & + kdH2Oice (i) / H2O(i,j,k) * 18._DKIND - endif - endif - if (iradtrans .eq. 1) then - if (idissZ .gt. 0) then - acoef = acoef - & + kdissH2O(i,j,k) - endif - endif - - H2Op(i) = ( scoef*dtit(i) + H2O(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - -C***** O2 ********** - scoef = 0._DKIND + 32._DKIND * ( 0._DKIND - & + kz31(i) * OI(i,j,k) * OI(i,j,k) / 256._DKIND - & + kz33(i) * OI(i,j,k) * OH(i,j,k) / 272._DKIND - & ) - acoef = 0._DKIND - & + kz19(i) * HI(i,j,k) - & + kz29(i) * CI(i,j,k) / 12._DKIND - & + kz36(i) * HII(i,j,k) - & + kz38(i) * CII(i,j,k) / 12._DKIND - & + kz53(i) * SiI(i,j,k) / 28._DKIND - - O2p(i) = ( scoef*dtit(i) + O2(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - -C***** SiI ********** - scoef = 0._DKIND + 28._DKIND * ( 0._DKIND - & ) - acoef = 0._DKIND - & + kz52(i) * OH(i,j,k) / 17._DKIND - & + kz53(i) * O2(i,j,k) / 32._DKIND - if ( ( igrgr .eq. 1 ) .or. ( idsub .eq. 1) ) then - if (idspecies .gt. 1) then - acoef = acoef - & + kdSiM (i) / SiI(i,j,k) * 28._DKIND - endif - endif - - SiIp(i) = ( scoef*dtit(i) + SiI(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - -C***** SiOI ********** - scoef = 0._DKIND + 44._DKIND * ( 0._DKIND - & + kz52(i) * SiI(i,j,k) * OH(i,j,k) / 476._DKIND - & + kz53(i) * SiI(i,j,k) * O2(i,j,k) / 896._DKIND - & ) - acoef = 0._DKIND - & + kz54(i) * OH(i,j,k) / 17._DKIND - if ( ( igrgr .eq. 1 ) .or. ( idsub .eq. 1) ) then - if (idspecies .gt. 0) then - acoef = acoef - & + kdMgSiO3 (i) / SiOI(i,j,k) * 44._DKIND - endif - if (idspecies .gt. 1) then - acoef = acoef - & + kdMg2SiO4 (i) / SiOI(i,j,k) * 44._DKIND - endif - endif - - SiOIp(i) = ( scoef*dtit(i) + SiOI(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - -C***** SiO2I ********** - scoef = 0._DKIND + 60._DKIND * ( 0._DKIND - & + kz54(i) * SiOI(i,j,k) * OH(i,j,k) / 748._DKIND - & ) - acoef = 0._DKIND - if ( ( igrgr .eq. 1 ) .or. ( idsub .eq. 1) ) then - if (idspecies .gt. 1) then - acoef = acoef - & + kdSiO2D (i) / SiO2I(i,j,k) * 60._DKIND - endif - endif - - SiO2Ip(i) = ( scoef*dtit(i) + SiO2I(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - ! MINOR BUT IMPORTANT SPECIES FOR MOLECULAR FORMATION -C***** CH ********** - scoef = 0._DKIND + 13._DKIND * ( 0._DKIND - & + kz16(i) * HI(i,j,k) * CH2(i,j,k) / 14._DKIND - & + kz20(i) * CI(i,j,k) * H2I(i,j,k) / 24._DKIND - & + kz27(i) * CI(i,j,k) * HI(i,j,k) / 12._DKIND - & ) - acoef = 0._DKIND - & + kz15(i) * HI(i,j,k) - & + kz23(i) * H2I(i,j,k) / 2._DKIND - & + kz32(i) * OI(i,j,k) / 16._DKIND - - CHp(i) = ( scoef*dtit(i) + CH(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - -C***** CH2 ********** - scoef = 0._DKIND + 14._DKIND * ( 0._DKIND - & + kz23(i) * H2I(i,j,k) * CH(i,j,k) / 26._DKIND - & + kz51(i) * H2I(i,j,k) * CI(i,j,k) / 24._DKIND - & ) - acoef = 0._DKIND - & + kz16(i) * HI(i,j,k) - if ( ( igrgr .eq. 1 ) .or. ( idsub .eq. 1) ) then - if (idspecies .gt. 2) then - acoef = acoef - & + kdreforg (i) / CH2(i,j,k) * 14._DKIND * 0.5_DKIND - endif - endif - - CH2p(i) = ( scoef*dtit(i) + CH2(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - -C***** COII ********** - scoef = 0._DKIND + 28._DKIND * ( 0._DKIND - & + kz37(i) * CII(i,j,k) * OH(i,j,k) / 204._DKIND - & ) - acoef = 0._DKIND - & + kz43(i) * HI(i,j,k) - - COIIp(i) = ( scoef*dtit(i) + COII(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - -C***** OII ********** - scoef = 0._DKIND + 16._DKIND * ( 0._DKIND - & + kz22(i) * HII(i,j,k) * OI(i,j,k) / 16._DKIND - & + kz38(i) * CII(i,j,k) * O2(i,j,k) / 384._DKIND - & ) - acoef = 0._DKIND - & + kz39(i) * HI(i,j,k) - & + kz40(i) * H2I(i,j,k) / 2._DKIND - & + kz45(i) * de(i,j,k) - if (iradtrans .eq. 1) then - if (iionZ .gt. 0) then - scoef = scoef - & + kphOI(i,j,k) * OI(i,j,k) - endif - endif - - OIIp(i) = ( scoef*dtit(i) + OII(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - -C***** OHII ********** - scoef = 0._DKIND + 17._DKIND * ( 0._DKIND - & + kz34(i) * HII(i,j,k) * OH(i,j,k) / 17._DKIND - & + kz40(i) * OII(i,j,k) * H2I(i,j,k) / 32._DKIND - & ) - acoef = 0._DKIND - & + kz41(i) * H2I(i,j,k) / 2._DKIND - - OHIIp(i) = ( scoef*dtit(i) + OHII(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - -C***** H2OII ********** - scoef = 0._DKIND + 18._DKIND * ( 0._DKIND - & + kz35(i) * HII(i,j,k) * H2O(i,j,k) / 18._DKIND - & + kz41(i) * OHII(i,j,k) * H2I(i,j,k) / 34._DKIND - & ) - acoef = 0._DKIND - & + kz42(i) * H2I(i,j,k) / 2._DKIND - & + kz46(i) * de(i,j,k) - & + kz47(i) * de(i,j,k) - - H2OIIp(i) = ( scoef*dtit(i) + H2OII(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - -C***** H3OII ********** - scoef = 0._DKIND + 19._DKIND * ( 0._DKIND - & + kz42(i) * H2OII(i,j,k) * H2I(i,j,k) / 36._DKIND - & ) - acoef = 0._DKIND - & + kz48(i) * de(i,j,k) - & + kz49(i) * de(i,j,k) - - H3OIIp(i) = ( scoef*dtit(i) + H3OII(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - -C***** O2II ********** - scoef = 0._DKIND + 32._DKIND * ( 0._DKIND - & + kz36(i) * HII(i,j,k) * O2(i,j,k) / 32._DKIND - & ) - acoef = 0._DKIND - & + kz50(i) * de(i,j,k) - - O2IIp(i) = ( scoef*dtit(i) + O2II(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - - if ( ( igrgr .eq. 1 ) .or. ( idsub .eq. 1) ) then - if (idspecies .gt. 0) then -C***** Mg ********** - scoef = 0._DKIND - acoef = 0._DKIND - acoef = acoef - & + kdMgSiO3 (i) / Mg(i,j,k) * 24._DKIND - if (idspecies .gt. 1) then - acoef = acoef - & + kdMg2SiO4 (i) / Mg(i,j,k) * 24._DKIND * 2._DKIND - & + kdMgO (i) / Mg(i,j,k) * 24._DKIND - endif - - Mgp(i) = ( scoef*dtit(i) + Mg(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - endif - - if (idspecies .gt. 1) then -C***** Al ********** - scoef = 0._DKIND - acoef = 0._DKIND - acoef = acoef - & + kdAl2O3 (i) / Al(i,j,k) * 27._DKIND * 2._DKIND - - Alp(i) = ( scoef*dtit(i) + Al(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - -C***** S ********** - scoef = 0._DKIND - acoef = 0._DKIND - acoef = acoef - & + kdFeS (i) / S(i,j,k) * 32._DKIND - - Sp(i) = ( scoef*dtit(i) + S(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - -C***** Fe ********** - scoef = 0._DKIND - acoef = 0._DKIND - acoef = acoef - & + kdFeM (i) / Fe(i,j,k) * 56._DKIND - & + kdFe3O4 (i) / Fe(i,j,k) * 56._DKIND * 3._DKIND - & + kdFeS (i) / Fe(i,j,k) * 56._DKIND - - Fep(i) = ( scoef*dtit(i) + Fe(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - endif - endif - - endif ! itmask_metal - enddo - endif ! imchem - -! --- (D4) Now do dust species --- -! - if ( ( igrgr .eq. 1 ) .or. ( idsub .eq. 1) ) then - do i = is+1, ie+1 - if (itmask_metal(i) .ne. MASK_FALSE) then - - if (idspecies .gt. 0) then -C***** MgSiO3 ********** - scoef = 0._DKIND - scoef = scoef - & + kdMgSiO3 (i) * 100._DKIND - acoef = 0._DKIND - - MgSiO3p(i) = ( scoef*dtit(i) + MgSiO3(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - -C***** AC ********** - scoef = 0._DKIND - scoef = scoef - & + kdAC (i) * 12._DKIND - acoef = 0._DKIND - - ACp(i) = ( scoef*dtit(i) + AC(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - endif - - if (idspecies .gt. 1) then -C***** SiM ********** - scoef = 0._DKIND - scoef = scoef - & + kdSiM (i) * 28._DKIND - acoef = 0._DKIND - - SiMp(i) = ( scoef*dtit(i) + SiM(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - -C***** FeM ********** - scoef = 0._DKIND - scoef = scoef - & + kdFeM (i) * 56._DKIND - acoef = 0._DKIND - - FeMp(i) = ( scoef*dtit(i) + FeM(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - -C***** Mg2SiO4 ********** - scoef = 0._DKIND - scoef = scoef - & + kdMg2SiO4 (i) * 140._DKIND - acoef = 0._DKIND - - Mg2SiO4p(i) = ( scoef*dtit(i) + Mg2SiO4(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - -C***** Fe3O4 ********** - scoef = 0._DKIND - scoef = scoef - & + kdFe3O4 (i) * 232._DKIND - acoef = 0._DKIND - - Fe3O4p(i) = ( scoef*dtit(i) + Fe3O4(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - -C***** SiO2D ********** - scoef = 0._DKIND - scoef = scoef - & + kdSiO2D (i) * 60._DKIND - acoef = 0._DKIND - - SiO2Dp(i) = ( scoef*dtit(i) + SiO2D(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - -C***** MgO ********** - scoef = 0._DKIND - scoef = scoef - & + kdMgO (i) * 40._DKIND - acoef = 0._DKIND - - MgOp(i) = ( scoef*dtit(i) + MgO(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - -C***** FeS ********** - scoef = 0._DKIND - scoef = scoef - & + kdFeS (i) * 88._DKIND - acoef = 0._DKIND - - FeSp(i) = ( scoef*dtit(i) + FeS(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - -C***** Al2O3 ********** - scoef = 0._DKIND - scoef = scoef - & + kdAl2O3 (i) * 102._DKIND - acoef = 0._DKIND - - Al2O3p(i) = ( scoef*dtit(i) + Al2O3(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - endif - - if (idspecies .gt. 2) then -C***** reforg ********** - scoef = 0._DKIND - scoef = scoef - & + kdreforg (i) * 22.68_DKIND - acoef = 0._DKIND - - reforgp(i) = ( scoef*dtit(i) + reforg(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - -C***** volorg ********** - scoef = 0._DKIND - scoef = scoef - & + kdvolorg (i) * 32._DKIND - acoef = 0._DKIND - - volorgp(i) = ( scoef*dtit(i) + volorg(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - -C***** H2Oice ********** - scoef = 0._DKIND - scoef = scoef - & + kdH2Oice (i) * 18._DKIND - acoef = 0._DKIND - - H2Oicep(i) = ( scoef*dtit(i) + H2Oice(i,j,k) ) - & / ( 1._DKIND + acoef*dtit(i) ) - - endif - - endif ! itmask_metal - enddo - endif ! igrgr or idsub - -! --- (E) Set densities from 1D temps to 3D fields --- - - do i = is+1, ie+1 - if (itmask(i) .ne. MASK_FALSE) then - HIdot_prev(i) = abs(HI(i,j,k)-HIp(i)) / - & max(real(dtit(i), DKIND), tiny8) - HI(i,j,k) = max(real(HIp(i), RKIND), tiny) - HII(i,j,k) = max(real(HIIp(i), RKIND), tiny) - HeI(i,j,k) = max(real(HeIp(i), RKIND), tiny) - HeII(i,j,k) = max(real(HeIIp(i), RKIND), tiny) - HeIII(i,j,k) = max(real(HeIIIp(i), RKIND), 1e-5_RKIND*tiny) - -! de(i,j,k) = dep(i) - -! Use charge conservation to determine electron fraction - - dedot_prev(i) = de(i,j,k) - de(i,j,k) = HII(i,j,k) + HeII(i,j,k)/4._RKIND + - & HeIII(i,j,k)/2._RKIND - if (ispecies .gt. 1) - & de(i,j,k) = de(i,j,k) - HM(i,j,k) + H2II(i,j,k)/2._RKIND - - if (ispecies .gt. 2) - & de(i,j,k) = de(i,j,k) + DII(i,j,k)/2._RKIND - if (ispecies .gt. 3) - & de(i,j,k) = de(i,j,k) - DM(i,j,k)/2._RKIND - & + HDII(i,j,k)/3._RKIND + HeHII(i,j,k)/5._RKIND - if ( (imchem .eq. 1) .and. - & (itmask_metal(i) .ne. MASK_FALSE) ) - & de(i,j,k) = de(i,j,k) - & + CII(i,j,k)/12._RKIND + COII(i,j,k)/28._RKIND - & + OII(i,j,k)/16._RKIND + OHII(i,j,k)/17._RKIND - & + H2OII(i,j,k)/18._RKIND + H3OII(i,j,k)/19._RKIND - & + O2II(i,j,k)/32._RKIND - - dedot_prev(i) = abs(de(i,j,k)-dedot_prev(i))/ - & max(dtit(i),tiny8) - - if (ispecies .gt. 1) then - HM(i,j,k) = max(real(HMp(i), RKIND), tiny) - H2I(i,j,k) = max(real(H2Ip(i),RKIND), tiny) - H2II(i,j,k) = max(real(H2IIp(i), RKIND), tiny) - endif - - if (ispecies .gt. 2) then - DI(i,j,k) = max(real(DIp(i), RKIND), tiny) - DII(i,j,k) = max(real(DIIp(i), RKIND), tiny) - HDI(i,j,k) = max(real(HDIp(i), RKIND), tiny) - endif - - if (ispecies .gt. 3) then - DM(i,j,k) = max(real(DMp(i), RKIND), tiny) - HDII(i,j,k) = max(real(HDIIp(i), RKIND), tiny) - HeHII(i,j,k) = max(real(HeHIIp(i), RKIND), tiny) - endif - - if ( (imchem .eq. 1) .and. - & (itmask_metal(i) .ne. MASK_FALSE) ) then - CI(i,j,k) = max(real(CIp(i) , RKIND), tiny) - CII(i,j,k) = max(real(CIIp(i) , RKIND), tiny) - CO(i,j,k) = max(real(COp(i) , RKIND), tiny) - CO2(i,j,k) = max(real(CO2p(i) , RKIND), tiny) - OI(i,j,k) = max(real(OIp(i) , RKIND), tiny) - OH(i,j,k) = max(real(OHp(i) , RKIND), tiny) - H2O(i,j,k) = max(real(H2Op(i) , RKIND), tiny) - O2(i,j,k) = max(real(O2p(i) , RKIND), tiny) - SiI(i,j,k) = max(real(SiIp(i) , RKIND), tiny) - SiOI(i,j,k) = max(real(SiOIp(i) , RKIND), tiny) - SiO2I(i,j,k) = max(real(SiO2Ip(i) , RKIND), tiny) - CH(i,j,k) = max(real(CHp(i) , RKIND), tiny) - CH2(i,j,k) = max(real(CH2p(i) , RKIND), tiny) - COII(i,j,k) = max(real(COIIp(i) , RKIND), tiny) - OII(i,j,k) = max(real(OIIp(i) , RKIND), tiny) - OHII(i,j,k) = max(real(OHIIp(i) , RKIND), tiny) - H2OII(i,j,k) = max(real(H2OIIp(i) , RKIND), tiny) - H3OII(i,j,k) = max(real(H3OIIp(i) , RKIND), tiny) - O2II(i,j,k) = max(real(O2IIp(i) , RKIND), tiny) - if ( ( igrgr .eq. 1 ) .or. ( idsub .eq. 1) ) then - if (idspecies .gt. 0) then - Mg(i,j,k) = max(real(Mgp(i) , RKIND), tiny) - endif - if (idspecies .gt. 1) then - Al(i,j,k) = max(real(Alp(i) , RKIND), tiny) - S(i,j,k) = max(real(Sp(i) , RKIND), tiny) - Fe(i,j,k) = max(real(Fep(i) , RKIND), tiny) - endif - endif - endif - - if ( ( igrgr .eq. 1 ) .or. ( idsub .eq. 1) ) then - if (idspecies .gt. 0) then - MgSiO3(i,j,k) = max(real(MgSiO3p(i) , RKIND), tiny) - AC(i,j,k) = max(real(ACp(i) , RKIND), tiny) - endif - if (idspecies .gt. 1) then - SiM(i,j,k) = max(real(SiMp(i) , RKIND), tiny) - FeM(i,j,k) = max(real(FeMp(i) , RKIND), tiny) - Mg2SiO4(i,j,k) = max(real(Mg2SiO4p(i), RKIND), tiny) - Fe3O4(i,j,k) = max(real(Fe3O4p(i) , RKIND), tiny) - SiO2D(i,j,k) = max(real(SiO2Dp(i) , RKIND), tiny) - MgO(i,j,k) = max(real(MgOp(i) , RKIND), tiny) - FeS(i,j,k) = max(real(FeSp(i) , RKIND), tiny) - Al2O3(i,j,k) = max(real(Al2O3p(i) , RKIND), tiny) - endif - if (idspecies .gt. 2) then - reforg(i,j,k) = max(real(reforgp(i) , RKIND), tiny) - volorg(i,j,k) = max(real(volorgp(i) , RKIND), tiny) - H2Oice(i,j,k) = max(real(H2Oicep(i) , RKIND), tiny) - endif - endif - - endif ! itmask -! - - if (HI(i,j,k) .ne. HI(i,j,k)) then -#ifdef _OPENMP -!$omp critical -#endif - write(*,*) 'HUGE HI! :: ', i, j, k, HI(i,j,k) -#ifdef _OPENMP -!$omp end critical -#endif - endif - - enddo ! end loop over i - - return - end diff --git a/src/clib/step_rate_gauss_seidel.hpp b/src/clib/step_rate_gauss_seidel.hpp new file mode 100644 index 000000000..3fe3928a6 --- /dev/null +++ b/src/clib/step_rate_gauss_seidel.hpp @@ -0,0 +1,260 @@ +//===----------------------------------------------------------------------===// +// +// See the LICENSE file for license and copyright information +// SPDX-License-Identifier: NCSA AND BSD-3-Clause +// +//===----------------------------------------------------------------------===// +/// +/// @file +/// Implements step_rate_gauss_seidel +/// +//===----------------------------------------------------------------------===// + +// This file was initially generated automatically during conversion of the +// step_rate_g function from FORTRAN to C++ + +#ifndef STEP_RATE_GAUSS_SEIDEL_HPP +#define STEP_RATE_GAUSS_SEIDEL_HPP + +#include +#include + +#include "grackle.h" // gr_float +#include "chemistry_solver_funcs.hpp" +#include "fortran_func_decls.h" // gr_mask_type +#include "index_helper.h" +#include "internal_types.hpp" +#include "LUT.hpp" + +namespace grackle::impl { + + +/// This function overwrites the species density field values using the values +/// in @p species_tmpdens, which holds the values from the end of the current +/// (sub-)cycle. +inline void update_fields_from_tmpdens_gauss_seidel( + const double* dtit, IndexRange idx_range, double* dedot_prev, + double* HIdot_prev, const gr_mask_type* itmask, + const gr_mask_type* itmask_metal, chemistry_data* my_chemistry, + grackle_field_data* my_fields, + grackle::impl::SpeciesCollection species_tmpdens +) { + + // Construct views of various species fields + // -> TODO: stop doing this! + + grackle::impl::View de(my_fields->e_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View HI(my_fields->HI_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View HII(my_fields->HII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View HeI(my_fields->HeI_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View HeII(my_fields->HeII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View HeIII(my_fields->HeIII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View HM(my_fields->HM_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View H2I(my_fields->H2I_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View H2II(my_fields->H2II_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View DI(my_fields->DI_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View DII(my_fields->DII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View HDI(my_fields->HDI_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View DM(my_fields->DM_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View HDII(my_fields->HDII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View HeHII(my_fields->HeHII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View CI(my_fields->CI_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View CII(my_fields->CII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View CO(my_fields->CO_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View CO2(my_fields->CO2_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View OI(my_fields->OI_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View OH(my_fields->OH_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View H2O(my_fields->H2O_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View O2(my_fields->O2_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View SiI(my_fields->SiI_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View SiOI(my_fields->SiOI_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View SiO2I(my_fields->SiO2I_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View CH(my_fields->CH_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View CH2(my_fields->CH2_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View COII(my_fields->COII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View OII(my_fields->OII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View OHII(my_fields->OHII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View H2OII(my_fields->H2OII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View H3OII(my_fields->H3OII_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View O2II(my_fields->O2II_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View Mg(my_fields->Mg_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View Al(my_fields->Al_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View S(my_fields->S_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View Fe(my_fields->Fe_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View SiM(my_fields->SiM_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View FeM(my_fields->FeM_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View Mg2SiO4(my_fields->Mg2SiO4_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View MgSiO3(my_fields->MgSiO3_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View Fe3O4(my_fields->Fe3O4_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View AC(my_fields->AC_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View SiO2D(my_fields->SiO2_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View MgO(my_fields->MgO_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View FeS(my_fields->FeS_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View Al2O3(my_fields->Al2O3_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View reforg(my_fields->ref_org_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View volorg(my_fields->vol_org_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View H2Oice(my_fields->H2O_ice_dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + + // --- (E) Set densities from 1D temps to 3D fields --- + + const int j = idx_range.j; + const int k = idx_range.k; + + for (int i = idx_range.i_start; i < idx_range.i_stop; i++) { + if (itmask[i] != MASK_FALSE) { + HIdot_prev[i] = std::fabs(HI(i,j,k)-species_tmpdens.data[SpLUT::HI][i]) / + std::fmax((double)(dtit[i] ), tiny8); + HI(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::HI][i] ), tiny_fortran_val); + HII(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::HII][i] ), tiny_fortran_val); + HeI(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::HeI][i] ), tiny_fortran_val); + HeII(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::HeII][i] ), tiny_fortran_val); + HeIII(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::HeIII][i] ), (gr_float)(1e-5)*tiny_fortran_val); + + // temporarily store the electron density from the start of the current + // subcycle. + dedot_prev[i] = de(i,j,k); + + // Use charge conservation to determine electron fraction + // -> in other words, we ignore species_tmpdens.data[SpLUT::e] (in + // practice, I think that array holds garbage values) + de(i,j,k) = HII(i,j,k) + HeII(i,j,k)/(gr_float)(4.) + + HeIII(i,j,k)/(gr_float)(2.); + if (my_chemistry->primordial_chemistry > 1) + { de(i,j,k) = de(i,j,k) - HM(i,j,k) + H2II(i,j,k)/(gr_float)(2.); } + + if (my_chemistry->primordial_chemistry > 2) + { de(i,j,k) = de(i,j,k) + DII(i,j,k)/(gr_float)(2.); } + if (my_chemistry->primordial_chemistry > 3) + { de(i,j,k) = de(i,j,k) - DM(i,j,k)/(gr_float)(2.) + + HDII(i,j,k)/(gr_float)(3.) + HeHII(i,j,k)/(gr_float)(5.); } + if ( (my_chemistry->metal_chemistry == 1) && + (itmask_metal[i] != MASK_FALSE) ) + { de(i,j,k) = de(i,j,k) + + CII(i,j,k)/(gr_float)(12.) + COII(i,j,k)/(gr_float)(28.) + + OII(i,j,k)/(gr_float)(16.) + OHII(i,j,k)/(gr_float)(17.) + + H2OII(i,j,k)/(gr_float)(18.) + H3OII(i,j,k)/(gr_float)(19.) + + O2II(i,j,k)/(gr_float)(32.); } + + // store the time-derivative of the electron-density in dedot + // (don't forget that we previously stored the value from the start of + // the current cycle within dedot_prev) + dedot_prev[i] = std::fabs(de(i,j,k)-dedot_prev[i])/ + std::fmax(dtit[i],tiny8); + + if (my_chemistry->primordial_chemistry > 1) { + HM(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::HM][i] ), tiny_fortran_val); + H2I(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::H2I][i]), tiny_fortran_val); + H2II(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::H2II][i] ), tiny_fortran_val); + } + + if (my_chemistry->primordial_chemistry > 2) { + DI(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::DI][i] ), tiny_fortran_val); + DII(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::DII][i] ), tiny_fortran_val); + HDI(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::HDI][i] ), tiny_fortran_val); + } + + if (my_chemistry->primordial_chemistry > 3) { + DM(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::DM][i] ), tiny_fortran_val); + HDII(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::HDII][i] ), tiny_fortran_val); + HeHII(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::HeHII][i] ), tiny_fortran_val); + } + + if ( (my_chemistry->metal_chemistry == 1) && + (itmask_metal[i] != MASK_FALSE) ) { + CI(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::CI][i] ), tiny_fortran_val); + CII(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::CII][i] ), tiny_fortran_val); + CO(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::CO][i] ), tiny_fortran_val); + CO2(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::CO2][i] ), tiny_fortran_val); + OI(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::OI][i] ), tiny_fortran_val); + OH(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::OH][i] ), tiny_fortran_val); + H2O(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::H2O][i] ), tiny_fortran_val); + O2(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::O2][i] ), tiny_fortran_val); + SiI(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::SiI][i] ), tiny_fortran_val); + SiOI(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::SiOI][i] ), tiny_fortran_val); + SiO2I(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::SiO2I][i] ), tiny_fortran_val); + CH(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::CH][i] ), tiny_fortran_val); + CH2(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::CH2][i] ), tiny_fortran_val); + COII(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::COII][i] ), tiny_fortran_val); + OII(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::OII][i] ), tiny_fortran_val); + OHII(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::OHII][i] ), tiny_fortran_val); + H2OII(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::H2OII][i] ), tiny_fortran_val); + H3OII(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::H3OII][i] ), tiny_fortran_val); + O2II(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::O2II][i] ), tiny_fortran_val); + if ( ( my_chemistry->grain_growth == 1 ) || ( my_chemistry->dust_sublimation == 1) ) { + if (my_chemistry->dust_species > 0) { + Mg(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::Mg][i] ), tiny_fortran_val); + } + if (my_chemistry->dust_species > 1) { + Al(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::Al][i] ), tiny_fortran_val); + S(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::S][i] ), tiny_fortran_val); + Fe(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::Fe][i] ), tiny_fortran_val); + } + } + } + + if ( ( my_chemistry->grain_growth == 1 ) || ( my_chemistry->dust_sublimation == 1) ) { + if (my_chemistry->dust_species > 0) { + MgSiO3(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::MgSiO3_dust][i] ), tiny_fortran_val); + AC(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::AC_dust][i] ), tiny_fortran_val); + } + if (my_chemistry->dust_species > 1) { + SiM(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::SiM_dust][i] ), tiny_fortran_val); + FeM(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::FeM_dust][i] ), tiny_fortran_val); + Mg2SiO4(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::Mg2SiO4_dust][i] ), tiny_fortran_val); + Fe3O4(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::Fe3O4_dust][i] ), tiny_fortran_val); + SiO2D(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::SiO2_dust][i] ), tiny_fortran_val); + MgO(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::MgO_dust][i] ), tiny_fortran_val); + FeS(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::FeS_dust][i] ), tiny_fortran_val); + Al2O3(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::Al2O3_dust][i] ), tiny_fortran_val); + } + if (my_chemistry->dust_species > 2) { + reforg(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::ref_org_dust][i] ), tiny_fortran_val); + volorg(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::vol_org_dust][i] ), tiny_fortran_val); + H2Oice(i,j,k) = std::fmax((gr_float)(species_tmpdens.data[SpLUT::H2O_ice_dust][i] ), tiny_fortran_val); + } + } + + } + // + + if (HI(i,j,k) != HI(i,j,k)) { + OMP_PRAGMA_CRITICAL + { + std::printf("HUGE HI! :: %d %d %d %g\n", + i, j, k, HI ( i, j, k )); + } + } + } + +} + +/// Uses one linearly implicit Gauss-Seidel sweep of a backward-Euler time +/// integrator to advance the rate equations by one (sub-)cycle (dtit). +inline void step_rate_gauss_seidel( + const double* dtit, IndexRange idx_range, gr_mask_type anydust, + const double* h2dust, const double* rhoH, double* dedot_prev, + double* HIdot_prev, const gr_mask_type* itmask, + const gr_mask_type* itmask_metal, chemistry_data* my_chemistry, + grackle_field_data* my_fields, photo_rate_storage my_uvb_rates, + grackle::impl::GrainSpeciesCollection grain_growth_rates, + grackle::impl::SpeciesCollection species_tmpdens, + grackle::impl::CollisionalRxnRateCollection kcol_buf, + grackle::impl::PhotoRxnRateCollection kshield_buf +) { + + // perform the Gauss-Seidel sweep to compute the species densities at the + // end of the current timestep. The results are saved in species_tmpdens + grackle::impl::chemistry::species_density_updates_gauss_seidel( + species_tmpdens, idx_range, dtit, anydust, h2dust, rhoH, itmask, + itmask_metal, my_chemistry, my_fields, my_uvb_rates, grain_growth_rates, + kcol_buf, kshield_buf); + + // update the entries from my_fields with the values in species_tmpdens + update_fields_from_tmpdens_gauss_seidel( + dtit, idx_range, dedot_prev, HIdot_prev, itmask, itmask_metal, + my_chemistry, my_fields, species_tmpdens); +} + +} // namespace grackle::impl + +#endif /* STEP_RATE_GAUSS_SEIDEL_HPP */ diff --git a/tests/unit/grtest_os.cpp b/tests/unit/grtest_os.cpp index eba65b404..a927d782a 100644 --- a/tests/unit/grtest_os.cpp +++ b/tests/unit/grtest_os.cpp @@ -86,7 +86,9 @@ std::unique_ptr grtest::CaptureSink::create(FILE* stream) { } // 5. prepare the container - std::unique_ptr out{new grtest::CaptureSink}; + + // NOLINTNEXTLINE(clang-analyzer-cplusplus.NewDeleteLeaks) + std::unique_ptr out(new grtest::CaptureSink); out->redirected_stream_ = stream; out->redirected_fd_ = underlying_fd; out->backup_fd_ = backup_fd; diff --git a/tests/unit/test_chemistry_struct_synced.cpp b/tests/unit/test_chemistry_struct_synced.cpp index 6a3bd9642..2c4f764ea 100644 --- a/tests/unit/test_chemistry_struct_synced.cpp +++ b/tests/unit/test_chemistry_struct_synced.cpp @@ -105,7 +105,11 @@ MemberTypeNameMapResult query_struct_members(std::string struct_name) "The build-system failed to create the xml-file. (THIS SHOULDN'T HAPPEN)", {} }; - } else if (std::fgetc(f) == EOF) { + } + bool file_is_empty = (std::fgetc(f) == EOF); + std::fclose(f); + + if (file_is_empty) { std::string msg( "The build-system wrote an empty xml-file, which usually means that " "castxml wasn't installed. If you install castxml, you need to instruct " @@ -114,7 +118,6 @@ MemberTypeNameMapResult query_struct_members(std::string struct_name) ); return MemberTypeNameMapResult{msg, {}}; } - std::fclose(f); // step 2: parse the xml // ===================== diff --git a/tests/unit/test_visitor.cpp b/tests/unit/test_visitor.cpp index 73e8af2ae..39e247469 100644 --- a/tests/unit/test_visitor.cpp +++ b/tests/unit/test_visitor.cpp @@ -42,8 +42,12 @@ TEST(VisitorTest, Simple) { MyDummyStruct tmp{nullptr, nullptr}; visit_member(&tmp, grackle::impl::visitor::AllocateMembers{ctx}); - ASSERT_NE(tmp.attr1, nullptr); - ASSERT_NE(tmp.attr2, nullptr); + EXPECT_NE(tmp.attr1, nullptr); + EXPECT_NE(tmp.attr2, nullptr); + if ((tmp.attr1 == nullptr) || (tmp.attr1 == nullptr)) { + visit_member(&tmp, grackle::impl::visitor::FreeMembers{}); + return; + } tmp.attr1[0] = 1.0; tmp.attr1[1] = 2.0; @@ -54,10 +58,10 @@ TEST(VisitorTest, Simple) { visit_member(&tmp2, grackle::impl::visitor::AllocateMembers{ctx}); visit_member_pair(tmp, tmp2, grackle::impl::visitor::CopyMembers{ctx}); - ASSERT_EQ(tmp2.attr1[0], 1.0); - ASSERT_EQ(tmp2.attr1[1], 2.0); - ASSERT_EQ(tmp2.attr2[0], -2); - ASSERT_EQ(tmp2.attr2[1], -1); + EXPECT_EQ(tmp2.attr1[0], 1.0); + EXPECT_EQ(tmp2.attr1[1], 2.0); + EXPECT_EQ(tmp2.attr2[0], -2); + EXPECT_EQ(tmp2.attr2[1], -1); visit_member(&tmp, grackle::impl::visitor::FreeMembers{}); visit_member(&tmp2, grackle::impl::visitor::FreeMembers{}); @@ -94,9 +98,14 @@ TEST(VisitorTest, Nested) { MyNestedStruct tmp{nullptr, {nullptr, nullptr}}; visit_member(&tmp, grackle::impl::visitor::AllocateMembers{ctx}); - ASSERT_NE(tmp.attr0, nullptr); - ASSERT_NE(tmp.dummy.attr1, nullptr); - ASSERT_NE(tmp.dummy.attr2, nullptr); + EXPECT_NE(tmp.attr0, nullptr); + EXPECT_NE(tmp.dummy.attr1, nullptr); + EXPECT_NE(tmp.dummy.attr2, nullptr); + if ((tmp.attr0 == nullptr) || (tmp.dummy.attr1 == nullptr) || + (tmp.dummy.attr2 == nullptr)) { + visit_member(&tmp, grackle::impl::visitor::FreeMembers{}); + return; + } tmp.attr0[0] = 100; tmp.attr0[1] = -100;