diff --git a/src/clib/CMakeLists.txt b/src/clib/CMakeLists.txt index 1de0d0304..b89f7035e 100644 --- a/src/clib/CMakeLists.txt +++ b/src/clib/CMakeLists.txt @@ -109,7 +109,9 @@ add_library(Grackle_Grackle # C++ Source (and Private Header Files) cool_multi_time_g-cpp.C cool_multi_time_g-cpp.h internal_types.C internal_types.hpp + solve_rate_cool_g-cpp.C solve_rate_cool_g-cpp.h utils-cpp.C utils-cpp.hpp + fortran_func_wrappers.hpp # Fortran Source Files calc_all_tdust_gasgr_1d_g.F diff --git a/src/clib/LUT.hpp b/src/clib/LUT.hpp new file mode 100644 index 000000000..64326ab3d --- /dev/null +++ b/src/clib/LUT.hpp @@ -0,0 +1,168 @@ +// See LICENSE file for license and copyright information + +/// @file LUT.hpp +/// @brief Declares some lookup-tables (LUTs) used internally by Grackle + +#ifndef LUT_HPP +#define LUT_HPP + +#ifndef __cplusplus +#error "This file can only be read by a c++ compiler" +#endif + +// once we have transcribed more code, we should really put this header's +// contents inside of the grackle::impl namespace +// - if a lookup-table is called {name}, then its fully qualified name will be +// `grackle::impl::{name}`. +// - when using the LUT in a function, we can shorten its name, within the +// function to just {name} by inserting `using grackle::impl::{name};` near +// the start of the function +// - we should hold off on doing this until more code is transcribed (since it +// will be hard for the transcription tools to automatically handle the +// shortenning of the fully qualified name) + + + +/// This is collection of enumerators (localized to the `SpLUT::` scope), with +/// an enumerator named for EVERY species (primordial-species, metal-species, +/// grain-species, etc). The enumerator values are intended to be used in a +/// compile-time lookup table (LUT), that incur no overhead at runtime. +/// +/// To understand this entity's purpose, imagine that we want to temporarily +/// track data for all N species known to Grackle (the data for each species is +/// stored in a buffer of common size determined at runtime). Consider the +/// following 3 ways we could accomplish this: +/// 1. define a separate variable holding a pointer for every single species +/// - this is what the fortran routines historically did +/// - this won't scale well with N (e.g. a function that requires all of +/// this info needs to take N arguments) +/// 2. define a single struct with a separate data member holding a pointer +/// for every single species. +/// - still doesn't scale well if for applying a common operation to the +/// data of each species (e.g. allocating/deallocating/zeroing-out) +/// 3. track each data-member in a statically sized array that holds pointers +/// for each species +/// - the array length is `SpLUT::NUM_ENTRIES` (known at compile time) +/// - it is now trivial to apply a common operation to each pointer +/// - we use the enumerator values in this LUT to lookup the value +/// corresponding to a given species (for example, we could use +/// `SpLUT::HI` to look up values related to HI) +/// +/// > [!important] +/// > This is a private implementation detail. We should never unintentionally +/// > expose this information directly in the public API (i.e. we want to +/// > maintain our flexibility). With that said, there are hypothetical +/// > scenarios where we might deliberately expose it (for performance purposes) +/// > -- but we should think long and hard about it before we do +/// +/// Implementation +/// -------------- +/// We implement choose to implement this entity by defining the enumerators +/// within a "stateless struct" (i.e. has a size of 0) that won't be +/// instantiated. The enclosing struct acts exactly like a pseudo-namespace +/// - doing this lets us use short/simple names for each enumerator without +/// fear of collision +/// - there are 2 advantages of using a struct over a namespace: +/// 1. it forces us to localize the definition in 1 location (in comparison, +/// C++ lets us add identifiers to a namespace in other files) +/// 2. it provides the flexibility to experiment with passing the LUT as a +/// template parameter +/// - a function that take accepts a specialized LUT (e.g. picked based +/// on primordial_chemistry) could reduce branching within loops and +/// and plausibly produce faster code. +/// - **WE DON'T CURRENTLY HAVE PLANS TO PURSUE THIS.** There are +/// much more immediately useful things to pursue first +/// +/// Other Thoughts +/// -------------- +/// - maybe we should store the dust species in a separate lookup table? +/// - adopting data-structures that leverage this LUT, is probably an important +/// stepping stone to implementing a system where reaction rates are +/// dynamically specified. After adopting such a system, we may not even need +/// a lookup table any more (this all assumes, of course, that such a system +/// is adequately performant) +struct SpLUT { + + // in the future, we may want to reimplement the following in terms of the + // XMacros provided in grackle_field_data_fdatamembers.def (or we may need to + // slightly revise the system?) + enum { + e, + HI, + HII, + HeI, + HeII, + HeIII, + + HM, + H2I, + H2II, + + DI, + DII, + HDI, + + DM, + HDII, + HeHII, + + 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, + + NUM_ENTRIES // <- always last (so it specifies the number of species) + }; // enum + +}; // SpLUT struct + + +/// Defines the LUT for Standard Collisional reaction rates +struct CollisionalRxnLUT { + + enum { + #define ENTRY(NAME) NAME, + #include "collisional_rxn_rate_members.def" + #undef ENTRY + + NUM_ENTRIES // <- always last (so it specifies the number of species) + }; // enum + +}; // CollisionalRxnLUT struct + + +#endif /* LUT_HPP */ diff --git a/src/clib/Make.config.objects b/src/clib/Make.config.objects index d567ece47..12be26063 100644 --- a/src/clib/Make.config.objects +++ b/src/clib/Make.config.objects @@ -38,6 +38,7 @@ OBJS_CONFIG_LIB = \ set_default_chemistry_parameters.lo \ solve_chemistry.lo \ solve_rate_cool_g.lo \ + solve_rate_cool_g-cpp.lo \ status_reporting.lo \ update_UVbackground_rates.lo \ lookup_cool_rates0d.lo \ diff --git a/src/clib/collisional_rxn_rate_members.def b/src/clib/collisional_rxn_rate_members.def new file mode 100644 index 000000000..141636222 --- /dev/null +++ b/src/clib/collisional_rxn_rate_members.def @@ -0,0 +1,132 @@ +/*********************************************************************** +/ +/ this file lists each member of the collisional species reaction rates +/ that we interpolate in 1d as a function of temperature. This list is +/ intended to be used with X-Macros (to reduce the amount of code +/ required to interact with these data-members) +/ +/ Currently, we just specify each rate name. In the future, we will +/ probably add other metadata to track the initializer function. +/ +/ Copyright (c) 2013, Enzo/Grackle Development Team. +/ +/ Distributed under the terms of the Enzo Public Licence. +/ +/ The full license is in the file LICENSE, distributed with this +/ software. +************************************************************************/ + + +/********************************** + * primordial chemistry rate data * + **********************************/ + +/* 6 species rates */ +ENTRY(k1) +ENTRY(k2) +ENTRY(k3) +ENTRY(k4) +ENTRY(k5) +ENTRY(k6) + +/* 9 species rates (including H2) */ +ENTRY(k7) +ENTRY(k8) +ENTRY(k9) +ENTRY(k10) +ENTRY(k11) +ENTRY(k12) +ENTRY(k13) +ENTRY(k14) +ENTRY(k15) +ENTRY(k16) +ENTRY(k17) +ENTRY(k18) +ENTRY(k19) +ENTRY(k20) /* currently not used */ +ENTRY(k21) /* currently not used */ +ENTRY(k22) /* 3-body H2 formation */ +ENTRY(k23) /* H2-H2 dissociation */ + +// k13dd is currently omitted from this list. It holds the density dependent +// version of k13 (collisional H2 dissociation) + +// radiative rates for 6-species (k24, k25, k26) are omitted +// radiative rates for 9-species (k27, k28, k29, k30, k31) are omitted + +// 12 species rates (with Deuterium) +ENTRY(k50) +ENTRY(k51) +ENTRY(k52) +ENTRY(k53) +ENTRY(k54) +ENTRY(k55) +ENTRY(k56) + +// New H-ionizing reactions, used for 6, 9 & 12 species chemistry +ENTRY(k57) +ENTRY(k58) + +// 15 species rates (with DM, HDII, HeHII) +ENTRY(k125) +ENTRY(k129) +ENTRY(k130) +ENTRY(k131) +ENTRY(k132) +ENTRY(k133) +ENTRY(k134) +ENTRY(k135) +ENTRY(k136) +ENTRY(k137) +ENTRY(k148) +ENTRY(k149) +ENTRY(k150) +ENTRY(k151) +ENTRY(k152) +ENTRY(k153) + +// Metal species chemistry rate data +// --------------------------------- + +ENTRY(kz15) +ENTRY(kz16) +ENTRY(kz17) +ENTRY(kz18) +ENTRY(kz19) +ENTRY(kz20) +ENTRY(kz21) +ENTRY(kz22) +ENTRY(kz23) +ENTRY(kz24) +ENTRY(kz25) +ENTRY(kz26) +ENTRY(kz27) +ENTRY(kz28) +ENTRY(kz29) +ENTRY(kz30) +ENTRY(kz31) +ENTRY(kz32) +ENTRY(kz33) +ENTRY(kz34) +ENTRY(kz35) +ENTRY(kz36) +ENTRY(kz37) +ENTRY(kz38) +ENTRY(kz39) +ENTRY(kz40) +ENTRY(kz41) +ENTRY(kz42) +ENTRY(kz43) +ENTRY(kz44) +ENTRY(kz45) +ENTRY(kz46) +ENTRY(kz47) +ENTRY(kz48) +ENTRY(kz49) +ENTRY(kz50) +ENTRY(kz51) +ENTRY(kz52) +ENTRY(kz53) +ENTRY(kz54) + + diff --git a/src/clib/fortran_func_decls.h b/src/clib/fortran_func_decls.h index d9782ac0f..4da25aa4d 100644 --- a/src/clib/fortran_func_decls.h +++ b/src/clib/fortran_func_decls.h @@ -18,6 +18,34 @@ typedef int32_t gr_mask_type; #define MASK_TRUE 1 #define MASK_FALSE 0 +void FORTRAN_NAME(calc_all_tdust_gasgr_1d_g)( + int* in, int* jn, int* kn, int* nratec, int* idustfield, int* is, int* ie, + int* j, int* k, double* fgr, double* gamma_isrfa, double* trad, + double* gasgra, long long* indixe, double* tdef, double* tgas, double* tdust, + double* metallicity, double* dust2gas, double* nh, double* gasgr_tdust, + gr_mask_type* itmask_metal, int* idspecies, int* itdmulti, int* gr_N, + int* gr_Size, double* gr_dT, double* gr_Td, double* tSiM, double* tFeM, + double* tMg2SiO4, double* tMgSiO3, double* tFe3O4, double* tAC, + double* tSiO2D, double* tMgO, double* tFeS, double* tAl2O3, double* treforg, + double* tvolorg, double* tH2Oice, double* gasgr2a, double* gamma_isrf2a, + double* coolunit, double* gasgr, double* myisrf, double* sgSiM, double* sgFeM, + double* sgMg2SiO4, double* sgMgSiO3, double* sgFe3O4, double* sgAC, + double* sgSiO2D, double* sgMgO, double* sgFeS, double* sgAl2O3, + double* sgreforg, double* sgvolorg, double* sgH2Oice, double* sgtot, + double* alSiM_data_ptr, double* alFeM_data_ptr, double* alMg2SiO4_data_ptr, + double* alMgSiO3_data_ptr, double* alFe3O4_data_ptr, double* alAC_data_ptr, + double* alSiO2D_data_ptr, double* alMgO_data_ptr, double* alFeS_data_ptr, + double* alAl2O3_data_ptr, double* alreforg_data_ptr, + double* alvolorg_data_ptr, double* alH2Oice_data_ptr, double* altot_data_ptr, + double* kpSiM, double* kpFeM, double* kpMg2SiO4, double* kpMgSiO3, + double* kpFe3O4, double* kpAC, double* kpSiO2D, double* kpMgO, double* kpFeS, + double* kpAl2O3, double* kpreforg, double* kpvolorg, double* kpH2Oice, + double* kptot, double* gasSiM, double* gasFeM, double* gasMg2SiO4, + double* gasMgSiO3, double* gasFe3O4, double* gasAC, double* gasSiO2D, + double* gasMgO, double* gasFeS, double* gasAl2O3, double* gasreforg, + double* gasvolorg, double* gasH2Oice +); + void FORTRAN_NAME(calc_grain_size_increment_1d)( int* immulti, int* imabund, int* idspecies, int* igrgr, gr_mask_type* itmask, int* in, int* jn, int* kn, int* is, int* ie, int* j, int* k, double* dom, @@ -488,133 +516,6 @@ void FORTRAN_NAME(scale_fields_g)( gr_float* metal_Y19_data_ptr ); -void FORTRAN_NAME(solve_rate_cool_g)( - int* icool, gr_float* d_data_ptr, gr_float* e_data_ptr, gr_float* u_data_ptr, - gr_float* v_data_ptr, gr_float* w_data_ptr, 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, int* in, int* jn, int* kn, - int* nratec, int* iexpand, int* ispecies, int* imetal, int* imcool, - int* idust, int* idustall, int* idustfield, int* idim, int* is, int* js, - int* ks, int* ie, int* je, int* ke, int* ih2co, int* ipiht, int* idustrec, - int* igammah, double* dx, double* dt, double* aye, double* temstart, - double* temend, double* utem, double* uxyz, double* uaye, double* urho, - double* utim, double* gamma, double* fh, double* dtoh, double* z_solar, - double* fgr, double* k1a, double* k2a, double* k3a, double* k4a, double* k5a, - double* k6a, double* k7a, double* k8a, double* k9a, double* k10a, - double* k11a, double* k12a, double* k13a, double* k13dda_data_ptr, - double* k14a, double* k15a, double* k16a, double* k17a, double* k18a, - double* k19a, double* k22a, double* k24, double* k25, double* k26, - double* k27, double* k28, double* k29, double* k30, double* k31, double* k50a, - double* k51a, double* k52a, double* k53a, double* k54a, double* k55a, - double* k56a, double* k57a, double* k58a, int* ndratec, double* dtemstart, - double* dtemend, double* h2dusta_data_ptr, double* ncrna, double* ncrd1a, - double* ncrd2a, double* ceHIa, double* ceHeIa, double* ceHeIIa, double* ciHIa, - double* ciHeIa, double* ciHeISa, double* ciHeIIa, double* reHIIa, - double* reHeII1a, double* reHeII2a, double* reHeIIIa, double* brema, - double* compa, double* gammaha, double* isrf, double* regra, - double* gamma_isrfa, double* comp_xraya, double* comp_temp, double* piHI, - double* piHeI, double* piHeII, 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, gr_float* metal_data_ptr, gr_float* dust_data_ptr, - double* hyd01ka, double* h2k01a, double* vibha, double* rotha, double* rotla, - double* gpldla, double* gphdla, double* hdltea, double* hdlowa, double* gaHIa, - double* gaH2a, double* gaHea, double* gaHpa, double* gaela, double* h2ltea, - double* gasgra, int* iH2shield, int* iradshield, double* avgsighi, - double* avgsighei, double* avgsigheii, int* iradtrans, int* iradcoupled, - int* iradstep, int* irt_honly, gr_float* kphHI_data_ptr, - gr_float* kphHeI_data_ptr, gr_float* kphHeII_data_ptr, - gr_float* kdissH2I_data_ptr, gr_float* photogamma_data_ptr, - gr_float* xH2shield_data_ptr, int* ierr, int* ih2optical, int* iciecool, - int* ithreebody, int* ih2cr, int* ihdcr, double* ciecoa, int* icmbTfloor, - int* iClHeat, double* clEleFra, long long* priGridRank, long long* priGridDim, - double* priPar1, double* priPar2, double* priPar3, double* priPar4, - double* priPar5, long long* priDataSize, double* priCooling, - double* priHeating, double* priMMW, long long* metGridRank, - long long* metGridDim, double* metPar1, double* metPar2, double* metPar3, - double* metPar4, double* metPar5, long long* metDataSize, double* metCooling, - double* metHeating, int* clnew, int* iVheat, int* iMheat, - gr_float* Vheat_data_ptr, gr_float* Mheat_data_ptr, int* iTfloor, - double* Tfloor_scalar, gr_float* Tfloor_data_ptr, int* imchem, int* igrgr, - int* ipcont, double* tmcool, gr_float* DM_data_ptr, gr_float* HDII_data_ptr, - gr_float* HeHII_data_ptr, 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* k125a, - double* k129a, double* k130a, double* k131a, double* k132a, double* k133a, - double* k134a, double* k135a, double* k136a, double* k137a, double* k148a, - double* k149a, double* k150a, double* k151a, double* k152a, double* k153a, - double* kz15a, double* kz16a, double* kz17a, double* kz18a, double* kz19a, - double* kz20a, double* kz21a, double* kz22a, double* kz23a, double* kz24a, - double* kz25a, double* kz26a, double* kz27a, double* kz28a, double* kz29a, - double* kz30a, double* kz31a, double* kz32a, double* kz33a, double* kz34a, - double* kz35a, double* kz36a, double* kz37a, double* kz38a, double* kz39a, - double* kz40a, double* kz41a, double* kz42a, double* kz43a, double* kz44a, - double* kz45a, double* kz46a, double* kz47a, double* kz48a, double* kz49a, - double* kz50a, double* kz51a, double* kz52a, double* kz53a, double* kz54a, - double* cieY06a, long long* LH2_N, long long* LH2_Size, double* LH2_D, - double* LH2_T, double* LH2_H, double* LH2_dD, double* LH2_dT, double* LH2_dH, - double* LH2_L, long long* LHD_N, long long* LHD_Size, double* LHD_D, - double* LHD_T, double* LHD_H, double* LHD_dD, double* LHD_dT, double* LHD_dH, - double* LHD_L, long long* LCI_N, long long* LCI_Size, double* LCI_D, - double* LCI_T, double* LCI_H, double* LCI_dD, double* LCI_dT, double* LCI_dH, - double* LCI_L, long long* LCII_N, long long* LCII_Size, double* LCII_D, - double* LCII_T, double* LCII_H, double* LCII_dD, double* LCII_dT, - double* LCII_dH, double* LCII_L, long long* LOI_N, long long* LOI_Size, - double* LOI_D, double* LOI_T, double* LOI_H, double* LOI_dD, double* LOI_dT, - double* LOI_dH, double* LOI_L, long long* LCO_N, long long* LCO_Size, - double* LCO_D, double* LCO_T, double* LCO_H, double* LCO_dD, double* LCO_dT, - double* LCO_dH, double* LCO_L, long long* LOH_N, long long* LOH_Size, - double* LOH_D, double* LOH_T, double* LOH_H, double* LOH_dD, double* LOH_dT, - double* LOH_dH, double* LOH_L, long long* LH2O_N, long long* LH2O_Size, - double* LH2O_D, double* LH2O_T, double* LH2O_H, double* LH2O_dD, - double* LH2O_dT, double* LH2O_dH, double* LH2O_L, long long* alphap_N, - long long* alphap_Size, double* alphap_D, double* alphap_T, double* alphap_dD, - double* alphap_dT, double* alphap_Data, int* immulti, int* imabund, - int* idspecies, int* itdmulti, int* idsub, gr_float* metal_loc_data_ptr, - gr_float* metal_C13_data_ptr, gr_float* metal_C20_data_ptr, - gr_float* metal_C25_data_ptr, gr_float* metal_C30_data_ptr, - gr_float* metal_F13_data_ptr, gr_float* metal_F15_data_ptr, - gr_float* metal_F50_data_ptr, gr_float* metal_F80_data_ptr, - gr_float* metal_P170_data_ptr, gr_float* metal_P200_data_ptr, - gr_float* metal_Y19_data_ptr, int* SN0_N, double* SN0_XC, double* SN0_XO, - double* SN0_XMg, double* SN0_XAl, double* SN0_XSi, double* SN0_XS, - double* SN0_XFe, double* SN0_fC, double* SN0_fO, double* SN0_fMg, - double* SN0_fAl, double* SN0_fSi, double* SN0_fS, double* SN0_fFe, - double* SN0_fSiM, double* SN0_fFeM, double* SN0_fMg2SiO4, double* SN0_fMgSiO3, - double* SN0_fFe3O4, double* SN0_fAC, double* SN0_fSiO2D, double* SN0_fMgO, - double* SN0_fFeS, double* SN0_fAl2O3, double* SN0_freforg, - double* SN0_fvolorg, double* SN0_fH2Oice, double* SN0_r0SiM_data_ptr, - double* SN0_r0FeM_data_ptr, double* SN0_r0Mg2SiO4_data_ptr, - double* SN0_r0MgSiO3_data_ptr, double* SN0_r0Fe3O4_data_ptr, - double* SN0_r0AC_data_ptr, double* SN0_r0SiO2D_data_ptr, - double* SN0_r0MgO_data_ptr, double* SN0_r0FeS_data_ptr, - double* SN0_r0Al2O3_data_ptr, double* SN0_r0reforg_data_ptr, - double* SN0_r0volorg_data_ptr, double* SN0_r0H2Oice_data_ptr, int* gr_N, - int* gr_Size, double* gr_dT, double* gr_Td, double* SN0_kpSiM_data_ptr, - double* SN0_kpFeM_data_ptr, double* SN0_kpMg2SiO4_data_ptr, - double* SN0_kpMgSiO3_data_ptr, double* SN0_kpFe3O4_data_ptr, - double* SN0_kpAC_data_ptr, double* SN0_kpSiO2D_data_ptr, - double* SN0_kpMgO_data_ptr, double* SN0_kpFeS_data_ptr, - double* SN0_kpAl2O3_data_ptr, double* SN0_kpreforg_data_ptr, - double* SN0_kpvolorg_data_ptr, double* SN0_kpH2Oice_data_ptr, - double* h2dustSa, double* h2dustCa, double* gasgr2a, double* gamma_isrf2a, - double* grogra, 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, int* iuseH2shield, int* iisrffield, - gr_float* isrf_habing_data_ptr, int* iH2shieldcustom, - gr_float* f_shield_custom_data_ptr, int* itmax, int* exititmax -); - void FORTRAN_NAME(ceiling_species_g)( gr_float* d_data_ptr, 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, @@ -674,7 +575,7 @@ void FORTRAN_NAME(lookup_cool_rates1d_g)( double* k29shield, double* k30shield, double* k31shield, double* h2dust, double* ncrn, double* ncrd1, double* ncrd2, double* t1, double* t2, double* tdef, double* logtem, long long* indixe, double* dom, - double* coolunit, double* tbase1, double* xbase1, double* dx_cgs, + double* coolunit, double* tbase1, double* uxyz, double* xbase1, double* dx_cgs, double* c_ljeans, int* iradtrans, gr_float* kdissH2I_data_ptr, gr_float* xH2shield_data_ptr, gr_mask_type* itmask, gr_mask_type* itmask_metal, double* fh, gr_float* metal_data_ptr, @@ -992,7 +893,7 @@ void FORTRAN_NAME(step_rate_newton_raphson)( gr_float* kphOI_data_ptr, int* idissZ, gr_float* kdissCO_data_ptr, gr_float* kdissOH_data_ptr, gr_float* kdissH2O_data_ptr, int* iuseH2shield, int* iisrffield, gr_float* isrf_habing_data_ptr, int* iH2shieldcustom, - gr_float* f_shield_custom_data_ptr, int* ierror, int* j, int* k, int* iter, + gr_float* f_shield_custom_data_ptr, int* j, int* k, int* iter, double* dom, double* comp1, double* comp2, double* coolunit, double* tbase1, double* xbase1, double* chunit, double* dx_cgs, double* c_ljeans, long long* indixe, double* t1, double* t2, double* logtem, double* tdef, @@ -1008,7 +909,7 @@ void FORTRAN_NAME(step_rate_newton_raphson)( double* reHeIII, double* brem, double* edot, double* hyd01k, double* h2k01, double* vibh, double* roth, double* rotl, double* gpldl, double* gphdl, double* hdlte, double* hdlow, double* cieco, gr_mask_type* anydust, - gr_mask_type* itmask_nr, gr_mask_type* itmask_metal, int* itr, int* imp_eng + gr_mask_type* itmask_nr, gr_mask_type* itmask_metal, int* imp_eng ); diff --git a/src/clib/fortran_func_wrappers.hpp b/src/clib/fortran_func_wrappers.hpp new file mode 100644 index 000000000..3796a04ec --- /dev/null +++ b/src/clib/fortran_func_wrappers.hpp @@ -0,0 +1,591 @@ +// See LICENSE file for license and copyright information + +/// @file fortran_func_wrappers.hpp +/// @brief Declares wrappers around routines that haven't been transcribed yet +/// from Fortran. +/// +/// THIS FILE WILL BE DELETED ONCE WE COMPLETE TRANSCRIPTION +/// +/// This file exists to aid the transcription process +/// - This file holds C++ functions that wrap untranscribed Fortran routines. +/// The idea is that these wrapper functions take reduced argument lists +/// (that may includes structs) before converting extended arg lists used +/// within the fortran subroutine. +/// - ideally the reduced argument list used by the wrapper function should +/// roughly approximate the final argument list that the routines will use +/// after transcription is complete. + +#ifndef FORTRAN_FUNC_WRAPPERS_HPP +#define FORTRAN_FUNC_WRAPPERS_HPP + +#ifndef __cplusplus +#error "This file must be read by a c++ compiler" +#endif + +#include "grackle.h" +#include "fortran_func_decls.h" +#include "index_helper.h" +#include "internal_types.hpp" +#include "internal_units.h" +#include "LUT.hpp" + + +// callers of these functions are generally expected to locally shorten the +// namespace name when they call these routines +namespace grackle::impl::fortran_wrapper { + +inline void ceiling_species_g( + int imetal, chemistry_data* my_chemistry, grackle_field_data* my_fields +) { + + FORTRAN_NAME(ceiling_species_g)(my_fields->density, 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->HM_density, my_fields->H2I_density, my_fields->H2II_density, my_fields->DI_density, my_fields->DII_density, my_fields->HDI_density, my_fields->metal_density, my_fields->dust_density, + &my_fields->grid_start[0], &my_fields->grid_end[0], &my_fields->grid_start[1], &my_fields->grid_end[1], &my_fields->grid_start[2], &my_fields->grid_end[2], + &my_fields->grid_dimension[0], &my_fields->grid_dimension[1], &my_fields->grid_dimension[2], &my_chemistry->primordial_chemistry, &imetal, &my_chemistry->use_dust_density_field, + my_fields->DM_density, my_fields->HDII_density, my_fields->HeHII_density, + &my_chemistry->metal_abundances, &my_chemistry->metal_chemistry, &my_chemistry->dust_species, &my_chemistry->multi_metals, + &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, + my_fields->local_ISM_metal_density, + my_fields->ccsn13_metal_density, my_fields->ccsn20_metal_density, my_fields->ccsn25_metal_density, my_fields->ccsn30_metal_density, + my_fields->fsn13_metal_density, my_fields->fsn15_metal_density, my_fields->fsn50_metal_density, my_fields->fsn80_metal_density, + my_fields->pisn170_metal_density, my_fields->pisn200_metal_density, my_fields->y19_metal_density); + +} + +inline void cool1d_multi_g( + int imetal, IndexRange idx_range, int iter, double* edot, double* tgas, + double* mmw, double* p2d, double* tdust, double* metallicity, + double* dust2gas, double* rhoH, gr_mask_type* itmask, + gr_mask_type* itmask_metal, chemistry_data* my_chemistry, + chemistry_data_storage* my_rates, grackle_field_data* my_fields, + photo_rate_storage my_uvb_rates, InternalGrUnits internalu, + grackle::impl::GrainSpeciesCollection grain_temperatures, + grackle::impl::LogTLinInterpScratchBuf logTlininterp_buf, + grackle::impl::Cool1DMultiScratchBuf cool1dmulti_buf, + grackle::impl::CoolHeatScratchBuf coolingheating_buf +) { + + // the following 2 variables should not be arguments (they are only inside + // of cool1d_multi_g to temporarily store a local value) + // + // We will fix this in the future + double comp1, comp2; + + // TODO: we should really pass in a value for min_metallicity (this is + // computed in a few independent places) + + FORTRAN_NAME(cool1d_multi_g)( + my_fields->density, my_fields->internal_energy, my_fields->x_velocity, my_fields->y_velocity, my_fields->z_velocity, 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->grid_dimension[0], &my_fields->grid_dimension[1], &my_fields->grid_dimension[2], &my_chemistry->NumberOfTemperatureBins, + &internalu.extfields_in_comoving, &my_chemistry->primordial_chemistry, &imetal, &my_chemistry->metal_cooling, + &my_chemistry->h2_on_dust, &my_chemistry->dust_chemistry, &my_chemistry->use_dust_density_field, &my_chemistry->dust_recombination_cooling, + &my_fields->grid_rank, &idx_range.i_start, &idx_range.i_end, &idx_range.jp1, &idx_range.kp1, &my_chemistry->ih2co, &my_chemistry->ipiht, &iter, &my_chemistry->photoelectric_heating, + &internalu.a_value, &my_chemistry->TemperatureStart, &my_chemistry->TemperatureEnd, &my_chemistry->SolarMetalFractionByMass, &my_chemistry->local_dust_to_gas_ratio, + &internalu.utem, &internalu.uxyz, &internalu.a_units, &internalu.urho, &internalu.tbase1, + &my_chemistry->Gamma, &my_chemistry->HydrogenFractionByMass, + my_rates->ceHI, my_rates->ceHeI, my_rates->ceHeII, my_rates->ciHI, my_rates->ciHeI, + my_rates->ciHeIS, my_rates->ciHeII, my_rates->reHII, my_rates->reHeII1, + my_rates->reHeII2, my_rates->reHeIII, my_rates->brem, &my_rates->comp, &my_rates->gammah, + &my_chemistry->interstellar_radiation_field, my_rates->regr, &my_rates->gamma_isrf, &my_uvb_rates.comp_xray, &my_uvb_rates.temp_xray, + &my_uvb_rates.piHI, &my_uvb_rates.piHeI, &my_uvb_rates.piHeII, &comp1, &comp2, + my_fields->HM_density, my_fields->H2I_density, my_fields->H2II_density, my_fields->DI_density, my_fields->DII_density, my_fields->HDI_density, my_fields->metal_density, my_fields->dust_density, + my_rates->hyd01k, my_rates->h2k01, my_rates->vibh, my_rates->roth, my_rates->rotl, + coolingheating_buf.hyd01k, coolingheating_buf.h2k01, coolingheating_buf.vibh, coolingheating_buf.roth, coolingheating_buf.rotl, + my_rates->GP99LowDensityLimit, my_rates->GP99HighDensityLimit, coolingheating_buf.gpldl, coolingheating_buf.gphdl, + my_rates->HDlte, my_rates->HDlow, coolingheating_buf.hdlte, coolingheating_buf.hdlow, + my_rates->GAHI, my_rates->GAH2, my_rates->GAHe, my_rates->GAHp, my_rates->GAel, + my_rates->H2LTE, my_rates->gas_grain, + coolingheating_buf.ceHI, coolingheating_buf.ceHeI, coolingheating_buf.ceHeII, coolingheating_buf.ciHI, coolingheating_buf.ciHeI, coolingheating_buf.ciHeIS, coolingheating_buf.ciHeII, + coolingheating_buf.reHII, coolingheating_buf.reHeII1, coolingheating_buf.reHeII2, coolingheating_buf.reHeIII, coolingheating_buf.brem, + logTlininterp_buf.indixe, logTlininterp_buf.t1, logTlininterp_buf.t2, logTlininterp_buf.logtem, logTlininterp_buf.tdef, edot, + tgas, cool1dmulti_buf.tgasold, mmw, p2d, tdust, metallicity, + dust2gas, rhoH, cool1dmulti_buf.mynh, cool1dmulti_buf.myde, + cool1dmulti_buf.gammaha_eff, cool1dmulti_buf.gasgr_tdust, cool1dmulti_buf.regr, + &my_chemistry->self_shielding_method, &my_uvb_rates.crsHI, &my_uvb_rates.crsHeI, &my_uvb_rates.crsHeII, + &my_uvb_rates.k24, &my_uvb_rates.k26, + &my_chemistry->use_radiative_transfer, my_fields->RT_heating_rate, + &my_chemistry->h2_optical_depth_approximation, &my_chemistry->cie_cooling, &my_chemistry->h2_cooling_rate, &my_chemistry->hd_cooling_rate, my_rates->cieco, coolingheating_buf.cieco, + &my_chemistry->cmb_temperature_floor, &my_chemistry->UVbackground, &my_chemistry->cloudy_electron_fraction_factor, + &my_rates->cloudy_primordial.grid_rank, my_rates->cloudy_primordial.grid_dimension, + my_rates->cloudy_primordial.grid_parameters[0], my_rates->cloudy_primordial.grid_parameters[1], my_rates->cloudy_primordial.grid_parameters[2], my_rates->cloudy_primordial.grid_parameters[3], my_rates->cloudy_primordial.grid_parameters[4], + &my_rates->cloudy_primordial.data_size, my_rates->cloudy_primordial.cooling_data, my_rates->cloudy_primordial.heating_data, my_rates->cloudy_primordial.mmw_data, + &my_rates->cloudy_metal.grid_rank, my_rates->cloudy_metal.grid_dimension, + my_rates->cloudy_metal.grid_parameters[0], my_rates->cloudy_metal.grid_parameters[1], my_rates->cloudy_metal.grid_parameters[2], my_rates->cloudy_metal.grid_parameters[3], my_rates->cloudy_metal.grid_parameters[4], + &my_rates->cloudy_metal.data_size, my_rates->cloudy_metal.cooling_data, my_rates->cloudy_metal.heating_data, &my_rates->cloudy_data_new, + &my_chemistry->use_volumetric_heating_rate, &my_chemistry->use_specific_heating_rate, my_fields->volumetric_heating_rate, my_fields->specific_heating_rate, + &my_chemistry->use_temperature_floor, &my_chemistry->temperature_floor_scalar, my_fields->temperature_floor, + &my_chemistry->use_isrf_field, my_fields->isrf_habing, itmask, itmask_metal, + &my_chemistry->metal_chemistry, &my_chemistry->grain_growth, &my_chemistry->use_primordial_continuum_opacity, &my_chemistry->tabulated_cooling_minimum_temperature, + my_fields->DM_density, my_fields->HDII_density, my_fields->HeHII_density, + 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, + my_rates->cieY06, + my_rates->LH2.props.dimension, &my_rates->LH2.props.data_size, + my_rates->LH2.props.parameters[0], my_rates->LH2.props.parameters[1], my_rates->LH2.props.parameters[2], + &my_rates->LH2.props.parameter_spacing[0], &my_rates->LH2.props.parameter_spacing[1], &my_rates->LH2.props.parameter_spacing[2], my_rates->LH2.data, + my_rates->LHD.props.dimension, &my_rates->LHD.props.data_size, + my_rates->LHD.props.parameters[0], my_rates->LHD.props.parameters[1], my_rates->LHD.props.parameters[2], + &my_rates->LHD.props.parameter_spacing[0], &my_rates->LHD.props.parameter_spacing[1], &my_rates->LHD.props.parameter_spacing[2], my_rates->LHD.data, + my_rates->LCI.props.dimension, &my_rates->LCI.props.data_size, + my_rates->LCI.props.parameters[0], my_rates->LCI.props.parameters[1], my_rates->LCI.props.parameters[2], + &my_rates->LCI.props.parameter_spacing[0], &my_rates->LCI.props.parameter_spacing[1], &my_rates->LCI.props.parameter_spacing[2], my_rates->LCI.data, + my_rates->LCII.props.dimension, &my_rates->LCII.props.data_size, + my_rates->LCII.props.parameters[0], my_rates->LCII.props.parameters[1], my_rates->LCII.props.parameters[2], + &my_rates->LCII.props.parameter_spacing[0], &my_rates->LCII.props.parameter_spacing[1], &my_rates->LCII.props.parameter_spacing[2], my_rates->LCII.data, + my_rates->LOI.props.dimension, &my_rates->LOI.props.data_size, + my_rates->LOI.props.parameters[0], my_rates->LOI.props.parameters[1], my_rates->LOI.props.parameters[2], + &my_rates->LOI.props.parameter_spacing[0], &my_rates->LOI.props.parameter_spacing[1], &my_rates->LOI.props.parameter_spacing[2], my_rates->LOI.data, + my_rates->LCO.props.dimension, &my_rates->LCO.props.data_size, + my_rates->LCO.props.parameters[0], my_rates->LCO.props.parameters[1], my_rates->LCO.props.parameters[2], + &my_rates->LCO.props.parameter_spacing[0], &my_rates->LCO.props.parameter_spacing[1], &my_rates->LCO.props.parameter_spacing[2], my_rates->LCO.data, + my_rates->LOH.props.dimension, &my_rates->LOH.props.data_size, + my_rates->LOH.props.parameters[0], my_rates->LOH.props.parameters[1], my_rates->LOH.props.parameters[2], + &my_rates->LOH.props.parameter_spacing[0], &my_rates->LOH.props.parameter_spacing[1], &my_rates->LOH.props.parameter_spacing[2], my_rates->LOH.data, + my_rates->LH2O.props.dimension, &my_rates->LH2O.props.data_size, + my_rates->LH2O.props.parameters[0], my_rates->LH2O.props.parameters[1], my_rates->LH2O.props.parameters[2], + &my_rates->LH2O.props.parameter_spacing[0], &my_rates->LH2O.props.parameter_spacing[1], &my_rates->LH2O.props.parameter_spacing[2], my_rates->LH2O.data, + my_rates->alphap.props.dimension, &my_rates->alphap.props.data_size, + my_rates->alphap.props.parameters[0], my_rates->alphap.props.parameters[1], &my_rates->alphap.props.parameter_spacing[0], &my_rates->alphap.props.parameter_spacing[1], + my_rates->alphap.data, + &my_chemistry->multi_metals, &my_chemistry->metal_abundances, &my_chemistry->dust_species, &my_chemistry->use_multiple_dust_temperatures, &my_chemistry->dust_sublimation, + my_fields->local_ISM_metal_density, + my_fields->ccsn13_metal_density, my_fields->ccsn20_metal_density, my_fields->ccsn25_metal_density, my_fields->ccsn30_metal_density, + my_fields->fsn13_metal_density, my_fields->fsn15_metal_density, my_fields->fsn50_metal_density, my_fields->fsn80_metal_density, + my_fields->pisn170_metal_density, my_fields->pisn200_metal_density, my_fields->y19_metal_density, + &my_rates->SN0_N, + my_rates->SN0_fSiM, my_rates->SN0_fFeM, my_rates->SN0_fMg2SiO4, my_rates->SN0_fMgSiO3, + my_rates->SN0_fFe3O4, my_rates->SN0_fAC, my_rates->SN0_fSiO2D, my_rates->SN0_fMgO, + my_rates->SN0_fFeS, my_rates->SN0_fAl2O3, + my_rates->SN0_freforg, my_rates->SN0_fvolorg, my_rates->SN0_fH2Oice, + my_rates->SN0_r0SiM, my_rates->SN0_r0FeM, my_rates->SN0_r0Mg2SiO4, my_rates->SN0_r0MgSiO3, + my_rates->SN0_r0Fe3O4, my_rates->SN0_r0AC, my_rates->SN0_r0SiO2D, my_rates->SN0_r0MgO, + my_rates->SN0_r0FeS, my_rates->SN0_r0Al2O3, + my_rates->SN0_r0reforg, my_rates->SN0_r0volorg, my_rates->SN0_r0H2Oice, + my_rates->gr_N, &my_rates->gr_Size, &my_rates->gr_dT, my_rates->gr_Td, + my_rates->SN0_kpSiM, my_rates->SN0_kpFeM, my_rates->SN0_kpMg2SiO4, my_rates->SN0_kpMgSiO3, + my_rates->SN0_kpFe3O4, my_rates->SN0_kpAC, my_rates->SN0_kpSiO2D, my_rates->SN0_kpMgO, + my_rates->SN0_kpFeS, my_rates->SN0_kpAl2O3, + my_rates->SN0_kpreforg, my_rates->SN0_kpvolorg, my_rates->SN0_kpH2Oice, + grain_temperatures.SiM, grain_temperatures.FeM, grain_temperatures.Mg2SiO4, grain_temperatures.MgSiO3, grain_temperatures.Fe3O4, + grain_temperatures.AC, grain_temperatures.SiO2D, grain_temperatures.MgO, grain_temperatures.FeS, grain_temperatures.Al2O3, + grain_temperatures.reforg, grain_temperatures.volorg, grain_temperatures.H2Oice, + my_rates->gas_grain2, &my_rates->gamma_isrf2 + ); + +} + +/// This routine uses the temperature to look up the chemical rates that are +/// tabulated in a log table as a function of temperature +inline void lookup_cool_rates1d_g( + IndexRange idx_range, gr_mask_type anydust, double* tgas1d, double* mmw, + double* tdust, double* dust2gas, double* k13dd, double* h2dust, + double dom, double dx_cgs, double c_ljeans, gr_mask_type* itmask, + gr_mask_type* itmask_metal, int imetal, gr_float* rhoH, double dt, + chemistry_data* my_chemistry, chemistry_data_storage* my_rates, + grackle_field_data* my_fields, photo_rate_storage my_uvb_rates, + InternalGrUnits internalu, + grackle::impl::GrainSpeciesCollection grain_growth_rates, + grackle::impl::GrainSpeciesCollection grain_temperatures, + grackle::impl::LogTLinInterpScratchBuf logTlininterp_buf, + grackle::impl::CollisionalRxnRateCollection kcr_buf, + grackle::impl::PhotoRxnRateCollection kshield_buf, + grackle::impl::ChemHeatingRates chemheatrates_buf +) { + + FORTRAN_NAME(lookup_cool_rates1d_g)(&my_chemistry->TemperatureStart, &my_chemistry->TemperatureEnd, &my_chemistry->NumberOfTemperatureBins, &idx_range.jp1, &idx_range.kp1, + &idx_range.i_start, &idx_range.i_end, &my_chemistry->three_body_rate, + &my_fields->grid_dimension[0], &my_fields->grid_dimension[1], &my_fields->grid_dimension[2], &my_chemistry->primordial_chemistry, &anydust, + &my_chemistry->H2_self_shielding, &my_chemistry->self_shielding_method, + tgas1d, mmw, my_fields->density, my_fields->HI_density, my_fields->HII_density, my_fields->HeI_density, my_fields->HeII_density, my_fields->HeIII_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, + tdust, dust2gas, + my_rates->k1, my_rates->k2, my_rates->k3, my_rates->k4, my_rates->k5, my_rates->k6, my_rates->k7, my_rates->k8, my_rates->k9, my_rates->k10, + my_rates->k11, my_rates->k12, my_rates->k13, my_rates->k13dd, my_rates->k14, my_rates->k15, my_rates->k16, + my_rates->k17, my_rates->k18, my_rates->k19, my_rates->k22, + my_rates->k50, my_rates->k51, my_rates->k52, my_rates->k53, my_rates->k54, my_rates->k55, my_rates->k56, + my_rates->k57, my_rates->k58, &my_chemistry->NumberOfDustTemperatureBins, &my_chemistry->DustTemperatureStart, &my_chemistry->DustTemperatureEnd, my_rates->h2dust, + my_rates->n_cr_n, my_rates->n_cr_d1, my_rates->n_cr_d2, + &my_uvb_rates.crsHI, &my_uvb_rates.crsHeI, &my_uvb_rates.crsHeII, &my_uvb_rates.piHI, &my_uvb_rates.piHeI, + 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.k28, &my_uvb_rates.k29, &my_uvb_rates.k30, &my_uvb_rates.k31, + 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], k13dd, kshield_buf.k24, kshield_buf.k25, kshield_buf.k26, + kshield_buf.k28, kshield_buf.k29, kshield_buf.k30, + kshield_buf.k31, h2dust, chemheatrates_buf.n_cr_n, chemheatrates_buf.n_cr_d1, chemheatrates_buf.n_cr_d2, + logTlininterp_buf.t1, logTlininterp_buf.t2, logTlininterp_buf.tdef, logTlininterp_buf.logtem, logTlininterp_buf.indixe, + &dom, &internalu.coolunit, &internalu.tbase1, &internalu.uxyz, &internalu.xbase1, &dx_cgs, &c_ljeans, + &my_chemistry->use_radiative_transfer, my_fields->RT_H2_dissociation_rate, my_fields->H2_self_shielding_length, itmask, + itmask_metal, + &my_chemistry->HydrogenFractionByMass, my_fields->metal_density, + my_fields->DM_density, my_fields->HDII_density, my_fields->HeHII_density, &imetal, &my_chemistry->metal_chemistry, &my_chemistry->grain_growth, + 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, + my_rates->k125, my_rates->k129, my_rates->k130, my_rates->k131, my_rates->k132, + my_rates->k133, my_rates->k134, my_rates->k135, my_rates->k136, my_rates->k137, + my_rates->k148, my_rates->k149, my_rates->k150, my_rates->k151, my_rates->k152, + my_rates->k153, + my_rates->kz15, my_rates->kz16, my_rates->kz17, my_rates->kz18, my_rates->kz19, + my_rates->kz20, my_rates->kz21, my_rates->kz22, my_rates->kz23, my_rates->kz24, + my_rates->kz25, my_rates->kz26, my_rates->kz27, my_rates->kz28, my_rates->kz29, + my_rates->kz30, my_rates->kz31, my_rates->kz32, my_rates->kz33, my_rates->kz34, + my_rates->kz35, my_rates->kz36, my_rates->kz37, my_rates->kz38, my_rates->kz39, + my_rates->kz40, my_rates->kz41, my_rates->kz42, my_rates->kz43, my_rates->kz44, + my_rates->kz45, my_rates->kz46, my_rates->kz47, my_rates->kz48, my_rates->kz49, + my_rates->kz50, my_rates->kz51, my_rates->kz52, my_rates->kz53, my_rates->kz54, + 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], + &my_chemistry->multi_metals, &my_chemistry->metal_abundances, &my_chemistry->dust_species, &my_chemistry->use_multiple_dust_temperatures, &my_chemistry->dust_sublimation, + my_fields->local_ISM_metal_density, + my_fields->ccsn13_metal_density, my_fields->ccsn20_metal_density, my_fields->ccsn25_metal_density, my_fields->ccsn30_metal_density, + my_fields->fsn13_metal_density, my_fields->fsn15_metal_density, my_fields->fsn50_metal_density, my_fields->fsn80_metal_density, + my_fields->pisn170_metal_density, my_fields->pisn200_metal_density, my_fields->y19_metal_density, + &my_rates->SN0_N, + my_rates->SN0_fSiM, my_rates->SN0_fFeM, my_rates->SN0_fMg2SiO4, my_rates->SN0_fMgSiO3, + my_rates->SN0_fFe3O4, my_rates->SN0_fAC, my_rates->SN0_fSiO2D, my_rates->SN0_fMgO, + my_rates->SN0_fFeS, my_rates->SN0_fAl2O3, + my_rates->SN0_freforg, my_rates->SN0_fvolorg, my_rates->SN0_fH2Oice, + my_rates->SN0_r0SiM, my_rates->SN0_r0FeM, my_rates->SN0_r0Mg2SiO4, my_rates->SN0_r0MgSiO3, + my_rates->SN0_r0Fe3O4, my_rates->SN0_r0AC, my_rates->SN0_r0SiO2D, my_rates->SN0_r0MgO, + my_rates->SN0_r0FeS, my_rates->SN0_r0Al2O3, + my_rates->SN0_r0reforg, my_rates->SN0_r0volorg, my_rates->SN0_r0H2Oice, + my_rates->gr_N, &my_rates->gr_Size, &my_rates->gr_dT, my_rates->gr_Td, + my_rates->SN0_kpSiM, my_rates->SN0_kpFeM, my_rates->SN0_kpMg2SiO4, my_rates->SN0_kpMgSiO3, + my_rates->SN0_kpFe3O4, my_rates->SN0_kpAC, my_rates->SN0_kpSiO2D, my_rates->SN0_kpMgO, + my_rates->SN0_kpFeS, my_rates->SN0_kpAl2O3, + my_rates->SN0_kpreforg, my_rates->SN0_kpvolorg, my_rates->SN0_kpH2Oice, + my_rates->h2dustS, my_rates->h2dustC, rhoH, my_rates->grain_growth_rate, &dt, + grain_growth_rates.SiM, grain_growth_rates.FeM, grain_growth_rates.Mg2SiO4, + grain_growth_rates.MgSiO3, grain_growth_rates.Fe3O4, grain_growth_rates.AC, grain_growth_rates.SiO2D, grain_growth_rates.MgO, grain_growth_rates.FeS, + grain_growth_rates.Al2O3, grain_growth_rates.reforg, grain_growth_rates.volorg, grain_growth_rates.H2Oice, + grain_temperatures.SiM, grain_temperatures.FeM, grain_temperatures.Mg2SiO4, grain_temperatures.MgSiO3, grain_temperatures.Fe3O4, + grain_temperatures.AC, grain_temperatures.SiO2D, grain_temperatures.MgO, grain_temperatures.FeS, grain_temperatures.Al2O3, + grain_temperatures.reforg, grain_temperatures.volorg, grain_temperatures.H2Oice, &my_chemistry->radiative_transfer_use_H2_shielding, + &my_chemistry->H2_custom_shielding, my_fields->H2_custom_shielding_factor + ); + +} + +// the following case was handcoded (so the argument order may shift when we +// actually transcribe the routine) +inline void make_consistent_g( + int imetal, double dom, chemistry_data* my_chemistry, + chemistry_data_storage* my_rates, grackle_field_data* my_fields +){ + FORTRAN_NAME(make_consistent_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->HM_density, my_fields->H2I_density, my_fields->H2II_density, my_fields->DI_density, my_fields->DII_density, my_fields->HDI_density, my_fields->metal_density, my_fields->dust_density, + my_fields->density, &my_fields->grid_start[0], &my_fields->grid_end[0], &my_fields->grid_start[1], &my_fields->grid_end[1], &my_fields->grid_start[2], &my_fields->grid_end[2], + &my_fields->grid_dimension[0], &my_fields->grid_dimension[1], &my_fields->grid_dimension[2], &my_chemistry->primordial_chemistry, &imetal, &my_chemistry->HydrogenFractionByMass, &my_chemistry->DeuteriumToHydrogenRatio, + &my_chemistry->use_dust_density_field, &my_chemistry->metal_chemistry, &my_chemistry->grain_growth, &dom, + my_fields->DM_density, my_fields->HDII_density, my_fields->HeHII_density, + 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, + &my_chemistry->multi_metals, &my_chemistry->metal_abundances, &my_chemistry->dust_species, &my_chemistry->use_multiple_dust_temperatures, &my_chemistry->dust_sublimation, + my_fields->local_ISM_metal_density, + my_fields->ccsn13_metal_density, my_fields->ccsn20_metal_density, my_fields->ccsn25_metal_density, my_fields->ccsn30_metal_density, + my_fields->fsn13_metal_density, my_fields->fsn15_metal_density, my_fields->fsn50_metal_density, my_fields->fsn80_metal_density, + my_fields->pisn170_metal_density, my_fields->pisn200_metal_density, my_fields->y19_metal_density, + &my_rates->SN0_N, + my_rates->SN0_XC, my_rates->SN0_XO, my_rates->SN0_XMg, my_rates->SN0_XAl, my_rates->SN0_XSi, + my_rates->SN0_XS, my_rates->SN0_XFe, + my_rates->SN0_fC, my_rates->SN0_fO, my_rates->SN0_fMg, my_rates->SN0_fAl, my_rates->SN0_fSi, + my_rates->SN0_fS, my_rates->SN0_fFe + ); +} + +/// This routine calculates the electron and HI rates of change in order to +/// determine the maximum permitted timestep +inline void rate_timestep_g( + double* dedot, double* HIdot, gr_mask_type anydust, IndexRange idx_range, + double* h2dust, double* rhoH, gr_mask_type* itmask, double* edot, + double chunit, double dom, chemistry_data* my_chemistry, + grackle_field_data* my_fields, photo_rate_storage my_uvb_rates, + grackle::impl::CollisionalRxnRateCollection kcr_buf, + grackle::impl::PhotoRxnRateCollection kshield_buf, + grackle::impl::ChemHeatingRates chemheatrates_buf +) { + FORTRAN_NAME(rate_timestep_g)( + dedot, HIdot, &my_chemistry->primordial_chemistry, &anydust, + 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->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, + 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, chemheatrates_buf.n_cr_n, chemheatrates_buf.n_cr_d1, chemheatrates_buf.n_cr_d2, rhoH, + kshield_buf.k24, kshield_buf.k25, kshield_buf.k26, + kshield_buf.k28, kshield_buf.k29, kshield_buf.k30, kshield_buf.k31, + &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, edot, &chunit, &dom, my_fields->metal_density, + my_fields->HDI_density, &my_chemistry->metal_chemistry, my_fields->CI_density, my_fields->OI_density, my_fields->OH_density, my_fields->CO_density, my_fields->H2O_density, + &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 + ); +} + + +inline void scale_fields_g( + int imetal, gr_float factor, chemistry_data* my_chemistry, + grackle_field_data* my_fields +) { + + FORTRAN_NAME(scale_fields_g)( + &my_chemistry->primordial_chemistry, &imetal, &my_chemistry->use_dust_density_field, &my_chemistry->metal_chemistry, + &my_chemistry->metal_abundances, &my_chemistry->dust_species, &my_chemistry->multi_metals, &my_chemistry->grain_growth, &my_chemistry->dust_sublimation, &factor, + &my_fields->grid_start[0], &my_fields->grid_end[0], &my_fields->grid_start[1], &my_fields->grid_end[1], &my_fields->grid_start[2], &my_fields->grid_end[2], &my_fields->grid_dimension[0], &my_fields->grid_dimension[1], &my_fields->grid_dimension[2], + my_fields->density, 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->HM_density, my_fields->H2I_density, my_fields->H2II_density, my_fields->DI_density, my_fields->DII_density, my_fields->HDI_density, my_fields->DM_density, my_fields->HDII_density, my_fields->HeHII_density, + my_fields->metal_density, my_fields->dust_density, + 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, + my_fields->local_ISM_metal_density, my_fields->ccsn13_metal_density, my_fields->ccsn20_metal_density, my_fields->ccsn25_metal_density, my_fields->ccsn30_metal_density, + my_fields->fsn13_metal_density, my_fields->fsn15_metal_density, my_fields->fsn50_metal_density, my_fields->fsn80_metal_density, + my_fields->pisn170_metal_density, my_fields->pisn200_metal_density, my_fields->y19_metal_density); + +} + +/// 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], species_tmpdens.data[SpLUT::FeM], species_tmpdens.data[SpLUT::Mg2SiO4], species_tmpdens.data[SpLUT::MgSiO3], species_tmpdens.data[SpLUT::Fe3O4], + species_tmpdens.data[SpLUT::AC], species_tmpdens.data[SpLUT::SiO2D], species_tmpdens.data[SpLUT::MgO], species_tmpdens.data[SpLUT::FeS], species_tmpdens.data[SpLUT::Al2O3], + species_tmpdens.data[SpLUT::reforg], species_tmpdens.data[SpLUT::volorg], species_tmpdens.data[SpLUT::H2Oice], + grain_growth_rates.SiM, grain_growth_rates.FeM, grain_growth_rates.Mg2SiO4, grain_growth_rates.MgSiO3, grain_growth_rates.Fe3O4, + grain_growth_rates.AC, grain_growth_rates.SiO2D, grain_growth_rates.MgO, grain_growth_rates.FeS, grain_growth_rates.Al2O3, + grain_growth_rates.reforg, grain_growth_rates.volorg, grain_growth_rates.H2Oice, + &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 + ); + +} + +inline void step_rate_newton_raphson( + int imetal, IndexRange idx_range, int iter, double dom, double chunit, + double dx_cgs, double c_ljeans, double* dtit, double* p2d, double* tgas, + double* tdust, double* metallicity, double* dust2gas, double* rhoH, + double* mmw, double* h2dust, double* edot, gr_mask_type anydust, + gr_mask_type* itmask_nr, gr_mask_type* itmask_metal, int* imp_eng, + chemistry_data* my_chemistry, chemistry_data_storage* my_rates, + grackle_field_data* my_fields, photo_rate_storage my_uvb_rates, + InternalGrUnits internalu, + grackle::impl::GrainSpeciesCollection grain_temperatures, + grackle::impl::LogTLinInterpScratchBuf logTlininterp_buf, + grackle::impl::Cool1DMultiScratchBuf cool1dmulti_buf, + grackle::impl::CoolHeatScratchBuf coolingheating_buf, + grackle::impl::ChemHeatingRates chemheatrates_buf +) { + + // the following 2 variables should not be arguments (they are only inside + // of cool1d_multi_g to temporarily store a local value) + // + // We will fix this in the future + double comp1, comp2; + + FORTRAN_NAME(step_rate_newton_raphson)(&my_chemistry->with_radiative_cooling, my_fields->density, my_fields->internal_energy, my_fields->x_velocity, my_fields->y_velocity, my_fields->z_velocity, 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->grid_dimension[0], &my_fields->grid_dimension[1], &my_fields->grid_dimension[2], &my_chemistry->NumberOfTemperatureBins, + &internalu.extfields_in_comoving, &my_chemistry->primordial_chemistry, &imetal, &my_chemistry->metal_cooling, &my_chemistry->h2_on_dust, + &my_chemistry->dust_chemistry, &my_chemistry->use_dust_density_field, &idx_range.i_start, &idx_range.i_end, &my_chemistry->ih2co, &my_chemistry->ipiht, + &my_chemistry->dust_recombination_cooling, &my_chemistry->photoelectric_heating, &internalu.a_value, &my_chemistry->TemperatureStart, &my_chemistry->TemperatureEnd, &internalu.utem, + &internalu.uxyz, &internalu.a_units, &internalu.urho, &internalu.tbase1, &my_chemistry->Gamma, &my_chemistry->HydrogenFractionByMass, &my_chemistry->SolarMetalFractionByMass, &my_chemistry->local_dust_to_gas_ratio, + my_rates->k1, my_rates->k2, my_rates->k3, my_rates->k4, my_rates->k5, my_rates->k6, my_rates->k7, my_rates->k8, my_rates->k9, + my_rates->k10, my_rates->k11, my_rates->k12, my_rates->k13, my_rates->k13dd, my_rates->k14, my_rates->k15, my_rates->k16, + my_rates->k17, my_rates->k18, my_rates->k19, my_rates->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, &my_uvb_rates.k31, my_rates->k50, my_rates->k51, my_rates->k52, my_rates->k53, my_rates->k54, + my_rates->k55, my_rates->k56, my_rates->k57, my_rates->k58, &my_chemistry->NumberOfDustTemperatureBins, &my_chemistry->DustTemperatureStart, + &my_chemistry->DustTemperatureEnd, my_rates->h2dust, my_rates->n_cr_n, my_rates->n_cr_d1, my_rates->n_cr_d2, my_rates->ceHI, + my_rates->ceHeI, my_rates->ceHeII, my_rates->ciHI, my_rates->ciHeI, my_rates->ciHeIS, my_rates->ciHeII, + my_rates->reHII, my_rates->reHeII1, my_rates->reHeII2, my_rates->reHeIII, my_rates->brem, + &my_rates->comp, &my_rates->gammah, &my_chemistry->interstellar_radiation_field, my_rates->regr, &my_rates->gamma_isrf, + &my_uvb_rates.comp_xray, &my_uvb_rates.temp_xray, &my_uvb_rates.piHI, &my_uvb_rates.piHeI, &my_uvb_rates.piHeII, my_fields->HM_density, + my_fields->H2I_density, my_fields->H2II_density, my_fields->DI_density, my_fields->DII_density, my_fields->HDI_density, my_fields->metal_density, my_fields->dust_density, my_rates->hyd01k, + my_rates->h2k01, my_rates->vibh, my_rates->roth, my_rates->rotl, my_rates->GP99LowDensityLimit, my_rates->GP99HighDensityLimit, + my_rates->HDlte, my_rates->HDlow, my_rates->GAHI, my_rates->GAH2, my_rates->GAHe, my_rates->GAHp, + my_rates->GAel, my_rates->H2LTE, my_rates->gas_grain, &my_chemistry->H2_self_shielding, &my_chemistry->self_shielding_method, + &my_uvb_rates.crsHI, &my_uvb_rates.crsHeI, &my_uvb_rates.crsHeII, &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, my_fields->RT_H2_dissociation_rate, + my_fields->RT_heating_rate, my_fields->H2_self_shielding_length, &my_chemistry->h2_optical_depth_approximation, &my_chemistry->cie_cooling, + &my_chemistry->three_body_rate, &my_chemistry->h2_cooling_rate, &my_chemistry->hd_cooling_rate, my_rates->cieco, &my_chemistry->cmb_temperature_floor, + &my_chemistry->UVbackground, &my_chemistry->cloudy_electron_fraction_factor, &my_rates->cloudy_primordial.grid_rank, my_rates->cloudy_primordial.grid_dimension, + my_rates->cloudy_primordial.grid_parameters[0], my_rates->cloudy_primordial.grid_parameters[1], my_rates->cloudy_primordial.grid_parameters[2], my_rates->cloudy_primordial.grid_parameters[3], my_rates->cloudy_primordial.grid_parameters[4], + &my_rates->cloudy_primordial.data_size, my_rates->cloudy_primordial.cooling_data, my_rates->cloudy_primordial.heating_data, my_rates->cloudy_primordial.mmw_data, + &my_rates->cloudy_metal.grid_rank, my_rates->cloudy_metal.grid_dimension, my_rates->cloudy_metal.grid_parameters[0], my_rates->cloudy_metal.grid_parameters[1], + my_rates->cloudy_metal.grid_parameters[2], my_rates->cloudy_metal.grid_parameters[3], my_rates->cloudy_metal.grid_parameters[4], &my_rates->cloudy_metal.data_size, + my_rates->cloudy_metal.cooling_data, my_rates->cloudy_metal.heating_data, &my_rates->cloudy_data_new, &my_chemistry->use_volumetric_heating_rate, &my_chemistry->use_specific_heating_rate, + my_fields->volumetric_heating_rate, my_fields->specific_heating_rate, &my_chemistry->use_temperature_floor, &my_chemistry->temperature_floor_scalar, my_fields->temperature_floor, + &my_chemistry->metal_chemistry, &my_chemistry->grain_growth, &my_chemistry->use_primordial_continuum_opacity, &my_chemistry->tabulated_cooling_minimum_temperature, my_fields->DM_density, my_fields->HDII_density, my_fields->HeHII_density, + 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, my_rates->k125, my_rates->k129, my_rates->k130, my_rates->k131, + my_rates->k132, my_rates->k133, my_rates->k134, my_rates->k135, my_rates->k136, my_rates->k137, my_rates->k148, + my_rates->k149, my_rates->k150, my_rates->k151, my_rates->k152, my_rates->k153, my_rates->kz15, my_rates->kz16, + my_rates->kz17, my_rates->kz18, my_rates->kz19, my_rates->kz20, my_rates->kz21, my_rates->kz22, my_rates->kz23, + my_rates->kz24, my_rates->kz25, my_rates->kz26, my_rates->kz27, my_rates->kz28, my_rates->kz29, my_rates->kz30, + my_rates->kz31, my_rates->kz32, my_rates->kz33, my_rates->kz34, my_rates->kz35, my_rates->kz36, my_rates->kz37, + my_rates->kz38, my_rates->kz39, my_rates->kz40, my_rates->kz41, my_rates->kz42, my_rates->kz43, my_rates->kz44, + my_rates->kz45, my_rates->kz46, my_rates->kz47, my_rates->kz48, my_rates->kz49, my_rates->kz50, my_rates->kz51, + my_rates->kz52, my_rates->kz53, my_rates->kz54, my_rates->cieY06, my_rates->LH2.props.dimension, &my_rates->LH2.props.data_size, + my_rates->LH2.props.parameters[0], my_rates->LH2.props.parameters[1], my_rates->LH2.props.parameters[2], &my_rates->LH2.props.parameter_spacing[0], &my_rates->LH2.props.parameter_spacing[1], &my_rates->LH2.props.parameter_spacing[2], + my_rates->LH2.data, my_rates->LHD.props.dimension, &my_rates->LHD.props.data_size, my_rates->LHD.props.parameters[0], my_rates->LHD.props.parameters[1], my_rates->LHD.props.parameters[2], + &my_rates->LHD.props.parameter_spacing[0], &my_rates->LHD.props.parameter_spacing[1], &my_rates->LHD.props.parameter_spacing[2], my_rates->LHD.data, my_rates->LCI.props.dimension, &my_rates->LCI.props.data_size, + my_rates->LCI.props.parameters[0], my_rates->LCI.props.parameters[1], my_rates->LCI.props.parameters[2], &my_rates->LCI.props.parameter_spacing[0], &my_rates->LCI.props.parameter_spacing[1], &my_rates->LCI.props.parameter_spacing[2], + my_rates->LCI.data, my_rates->LCII.props.dimension, &my_rates->LCII.props.data_size, my_rates->LCII.props.parameters[0], my_rates->LCII.props.parameters[1], my_rates->LCII.props.parameters[2], + &my_rates->LCII.props.parameter_spacing[0], &my_rates->LCII.props.parameter_spacing[1], &my_rates->LCII.props.parameter_spacing[2], my_rates->LCII.data, my_rates->LOI.props.dimension, + &my_rates->LOI.props.data_size, my_rates->LOI.props.parameters[0], my_rates->LOI.props.parameters[1], my_rates->LOI.props.parameters[2], &my_rates->LOI.props.parameter_spacing[0], &my_rates->LOI.props.parameter_spacing[1], + &my_rates->LOI.props.parameter_spacing[2], my_rates->LOI.data, my_rates->LCO.props.dimension, &my_rates->LCO.props.data_size, my_rates->LCO.props.parameters[0], my_rates->LCO.props.parameters[1], + my_rates->LCO.props.parameters[2], &my_rates->LCO.props.parameter_spacing[0], &my_rates->LCO.props.parameter_spacing[1], &my_rates->LCO.props.parameter_spacing[2], my_rates->LCO.data, my_rates->LOH.props.dimension, + &my_rates->LOH.props.data_size, my_rates->LOH.props.parameters[0], my_rates->LOH.props.parameters[1], my_rates->LOH.props.parameters[2], &my_rates->LOH.props.parameter_spacing[0], &my_rates->LOH.props.parameter_spacing[1], + &my_rates->LOH.props.parameter_spacing[2], my_rates->LOH.data, my_rates->LH2O.props.dimension, &my_rates->LH2O.props.data_size, my_rates->LH2O.props.parameters[0], my_rates->LH2O.props.parameters[1], + my_rates->LH2O.props.parameters[2], &my_rates->LH2O.props.parameter_spacing[0], &my_rates->LH2O.props.parameter_spacing[1], &my_rates->LH2O.props.parameter_spacing[2], my_rates->LH2O.data, + my_rates->alphap.props.dimension, &my_rates->alphap.props.data_size, my_rates->alphap.props.parameters[0], my_rates->alphap.props.parameters[1], + &my_rates->alphap.props.parameter_spacing[0], &my_rates->alphap.props.parameter_spacing[1], my_rates->alphap.data, &my_chemistry->multi_metals, + &my_chemistry->metal_abundances, &my_chemistry->dust_species, &my_chemistry->use_multiple_dust_temperatures, &my_chemistry->dust_sublimation, my_fields->local_ISM_metal_density, + my_fields->ccsn13_metal_density, my_fields->ccsn20_metal_density, my_fields->ccsn25_metal_density, my_fields->ccsn30_metal_density, + my_fields->fsn13_metal_density, my_fields->fsn15_metal_density, my_fields->fsn50_metal_density, my_fields->fsn80_metal_density, + my_fields->pisn170_metal_density, my_fields->pisn200_metal_density, my_fields->y19_metal_density, &my_rates->SN0_N, + my_rates->SN0_fSiM, my_rates->SN0_fFeM, my_rates->SN0_fMg2SiO4, my_rates->SN0_fMgSiO3, + my_rates->SN0_fFe3O4, my_rates->SN0_fAC, my_rates->SN0_fSiO2D, my_rates->SN0_fMgO, + my_rates->SN0_fFeS, my_rates->SN0_fAl2O3, my_rates->SN0_freforg, my_rates->SN0_fvolorg, + my_rates->SN0_fH2Oice, my_rates->SN0_r0SiM, my_rates->SN0_r0FeM, my_rates->SN0_r0Mg2SiO4, + my_rates->SN0_r0MgSiO3, my_rates->SN0_r0Fe3O4, my_rates->SN0_r0AC, my_rates->SN0_r0SiO2D, + my_rates->SN0_r0MgO, my_rates->SN0_r0FeS, my_rates->SN0_r0Al2O3, my_rates->SN0_r0reforg, + my_rates->SN0_r0volorg, my_rates->SN0_r0H2Oice, my_rates->gr_N, &my_rates->gr_Size, &my_rates->gr_dT, + my_rates->gr_Td, my_rates->SN0_kpSiM, my_rates->SN0_kpFeM, my_rates->SN0_kpMg2SiO4, + my_rates->SN0_kpMgSiO3, my_rates->SN0_kpFe3O4, my_rates->SN0_kpAC, my_rates->SN0_kpSiO2D, + my_rates->SN0_kpMgO, my_rates->SN0_kpFeS, my_rates->SN0_kpAl2O3, my_rates->SN0_kpreforg, + my_rates->SN0_kpvolorg, my_rates->SN0_kpH2Oice, my_rates->h2dustS, my_rates->h2dustC, + my_rates->gas_grain2, &my_rates->gamma_isrf2, my_rates->grain_growth_rate, &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, &my_chemistry->radiative_transfer_use_H2_shielding, &my_chemistry->use_isrf_field, + my_fields->isrf_habing, &my_chemistry->H2_custom_shielding, my_fields->H2_custom_shielding_factor, + &idx_range.jp1, &idx_range.kp1, &iter, &dom, &comp1, + &comp2, &internalu.coolunit, &internalu.tbase1, &internalu.xbase1, &chunit, &dx_cgs, + &c_ljeans, logTlininterp_buf.indixe, logTlininterp_buf.t1, logTlininterp_buf.t2, logTlininterp_buf.logtem, logTlininterp_buf.tdef, dtit, + p2d, tgas, cool1dmulti_buf.tgasold, tdust, metallicity, dust2gas, + rhoH, mmw, cool1dmulti_buf.mynh, cool1dmulti_buf.myde, cool1dmulti_buf.gammaha_eff, cool1dmulti_buf.gasgr_tdust, + cool1dmulti_buf.regr, h2dust, chemheatrates_buf.n_cr_n, chemheatrates_buf.n_cr_d1, chemheatrates_buf.n_cr_d2, grain_temperatures.SiM, grain_temperatures.FeM, + grain_temperatures.Mg2SiO4, grain_temperatures.MgSiO3, grain_temperatures.Fe3O4, grain_temperatures.AC, grain_temperatures.SiO2D, grain_temperatures.MgO, + grain_temperatures.FeS, grain_temperatures.Al2O3, grain_temperatures.reforg, grain_temperatures.volorg, grain_temperatures.H2Oice, coolingheating_buf.ceHI, + coolingheating_buf.ceHeI, coolingheating_buf.ceHeII, coolingheating_buf.ciHI, coolingheating_buf.ciHeI, coolingheating_buf.ciHeIS, coolingheating_buf.ciHeII, + coolingheating_buf.reHII, coolingheating_buf.reHeII1, coolingheating_buf.reHeII2, coolingheating_buf.reHeIII, coolingheating_buf.brem, edot, + coolingheating_buf.hyd01k, coolingheating_buf.h2k01, coolingheating_buf.vibh, coolingheating_buf.roth, coolingheating_buf.rotl, coolingheating_buf.gpldl, coolingheating_buf.gphdl, + coolingheating_buf.hdlte, coolingheating_buf.hdlow, coolingheating_buf.cieco, &anydust, itmask_nr, + itmask_metal, imp_eng + ); + + + + +} + +} // namespace grackle::impl::fortran_wrapper + +#endif /* FORTRAN_FUNC_WRAPPERS_HPP */ + diff --git a/src/clib/index_helper.h b/src/clib/index_helper.h index 143f6d5b8..b08677af3 100644 --- a/src/clib/index_helper.h +++ b/src/clib/index_helper.h @@ -60,6 +60,77 @@ typedef struct int end; } field_flat_index_range; +/// Specifies the range of indices for grackle's 3D fields. +/// +/// Essentially, you use this struct to easily get the 3D indices for elements +/// in an "islice". +/// - Grackle commonly uses buffers that are used to hold temporary variable +/// that you just use the i-indices to access +/// - If you want to use a 3D index to get an element from one of grackle's 3D +/// field, special logic is required to remap the 3 indices to a single index +/// in the underlying flat 1d buffer that holds a field's data. In some cases +/// (e.g. Fortran Arrays) the remapping is done behind the scenes and in +/// cases, it is more explicit (but the logic is always somewhere) +/// - For context, the `field_flat_index_range` can be used to specify the same +/// range for 3D fields, but the remapping logic is pre-applied (this makes +/// it much less useful when you have these custom buffers) +/// +/// @par Role in Transcription +/// This is intended to aide the transcription process of fortran subroutines +/// - prior to transcription, the 1-based `j` and `k` index values are passed +/// throughout the call-stack. Furthermore, the 0-based `i_start` & `i_end` +/// values are passed around. +/// - the automatic transcription process ends up making the indexing look a +/// little messy when transcribing from Fortran to C++ (in order to account +/// for the difference between 1-based & 0-based indexing). Thus gradual, +/// manual intervention is required to make indexing appear idiomatic. +/// - the manual process introduces a lot of opportunities for error if we +/// aren't extremely careful with keeping track of whether a local variable +/// still refers to a 1-based index or been converted to a 0-based index. +/// By using this struct, we can be very explicit and certain about whether +/// an expressions refers to a 1-based index or a 0-based index' +/// +/// @note +/// Naively, this all might sound like it involves a bunch of extra work, but +/// the following it will work really well with the transcription machinery: +/// - Consider a fortran subroutine, `my_subroutine`, where one of the +/// arguments is the `j` index (which is expected to use 1-based indexing) +/// and the signature looks like `void my_subroutine(..., int* j, ...)` +/// - Now suppose we have a C/C++ function that calls `my_subroutine` and that +/// there is an `IndexRange` in the function called `idx_range`. +/// - If the call to `my_subroutine` in the C/C++ function looks like +/// `my_subroutine(..., &idx_range.jp1, ...)`, then when applying the +/// transcription tools to `my_subroutine` we can have the transcription +/// machinery replace every reference to `j` in the transcribed function +/// with `idx_range.jp1` (it will also modify the argument list properly). +/// +/// @note +/// ultimately, it will be more idiomatic to make for-loops use `i_stop` rather +/// than `i_end`. Obviously, this is a **very low** priority (but this is why +/// we put both variables in this struct) +/// +/// @par Future Usage +/// If we continue using this type after we complete transcription, we can have +/// it take on a prominent role in adding GPU-support. +typedef struct IndexRange +{ + // specifies the fixed (0-based) j and k indices to be used with the range + int j; + int k; + + // j & k indices for 1-based indexing (remove when transcription is done) + int jp1; + int kp1; + + // specifies the starting and stopping (0-based) indices specifying the range + // along the i-axis (the stopping index isn't included in the range) + int i_start; + int i_stop; + + // holds the value of `i_stop-1` (remove when transcription is done) + int i_end; +} IndexRange; + /*********************************************************************** / / FUNCTION DECLARATIONS @@ -81,6 +152,23 @@ static inline field_flat_index_range inner_flat_range_( return out; } +/// constructs an IndexRange, which holds the 3D index information for an +/// "islice." +static inline IndexRange make_idx_range_( + int outer_index, const grackle_index_helper* idx_helper +) +{ + IndexRange out; + out.k = (outer_index / idx_helper->num_j_inds) + idx_helper->k_start; + out.j = (outer_index % idx_helper->num_j_inds) + idx_helper->j_start; + out.jp1 = out.j+1; + out.kp1 = out.k+1; + out.i_start = idx_helper->i_start; + out.i_end = idx_helper->i_end; + out.i_stop = idx_helper->i_end+1; + return out; +} + grackle_index_helper build_index_helper_(const grackle_field_data *my_fields); #ifdef __cplusplus diff --git a/src/clib/internal_types.C b/src/clib/internal_types.C index d5f4ca660..5ca6ab01a 100644 --- a/src/clib/internal_types.C +++ b/src/clib/internal_types.C @@ -191,4 +191,109 @@ void grackle::impl::drop_GrainSpeciesCollection( for_each_grainspeciescol_member(ptr, &cleanup_member_, NULL); } +// ----------------------------------------------------------------- + +grackle::impl::SpeciesCollection grackle::impl::new_SpeciesCollection( + int nelem +) { + GRIMPL_REQUIRE(nelem > 0, "nelem must be positive"); + grackle::impl::SpeciesCollection out; + double* ptr = (double*)malloc(sizeof(double) * nelem * SpLUT::NUM_ENTRIES); + for (int i = 0; i < SpLUT::NUM_ENTRIES; i++) { + out.data[i] = ptr + (i * nelem); + } + return out; +} + +void grackle::impl::drop_SpeciesCollection( + grackle::impl::SpeciesCollection *ptr +) { + // since we only allocate a single pointer, we only need to call free once + free(ptr->data[0]); +} + +// ----------------------------------------------------------------- + +grackle::impl::CollisionalRxnRateCollection grackle::impl::new_CollisionalRxnRateCollection( + int nelem +) { + GRIMPL_REQUIRE(nelem > 0, "nelem must be positive"); + grackle::impl::CollisionalRxnRateCollection out; + double* ptr = (double*)malloc( + sizeof(double) * nelem * CollisionalRxnLUT::NUM_ENTRIES + ); + for (int i = 0; i < CollisionalRxnLUT::NUM_ENTRIES; i++) { + out.data[i] = ptr + (i * nelem); + } + return out; +} + +void grackle::impl::drop_CollisionalRxnRateCollection( + grackle::impl::CollisionalRxnRateCollection *ptr +) { + // since we only allocate a single pointer, we only need to call free once + free(ptr->data[0]); +} + +// ----------------------------------------------------------------- + +/// Apply a function to each data member of PhotoRxnRateCollection +/// +/// @note +/// If we are willing to embrace C++, then this should accept a template +/// argument (instead of a function pointer and the callback_ctx pointer) +static void for_each_kphotoCollection_member( + grackle::impl::PhotoRxnRateCollection* ptr, + modify_member_callback* fn, + void* callback_ctx +) { + fn(&ptr->k24, callback_ctx); + fn(&ptr->k25, callback_ctx); + fn(&ptr->k26, callback_ctx); + fn(&ptr->k27, callback_ctx); + fn(&ptr->k28, callback_ctx); + fn(&ptr->k29, callback_ctx); + fn(&ptr->k30, callback_ctx); + fn(&ptr->k31, callback_ctx); +} + +grackle::impl::PhotoRxnRateCollection grackle::impl::new_PhotoRxnRateCollection( + int nelem +) { + GRIMPL_REQUIRE(nelem > 0, "nelem must be positive"); + grackle::impl::PhotoRxnRateCollection out; + MemberAllocCtx_ ctx{nelem}; + for_each_kphotoCollection_member(&out, &allocate_member_, (void*)(&ctx)); + return out; +} + +void grackle::impl::drop_PhotoRxnRateCollection( + grackle::impl::PhotoRxnRateCollection* ptr +) +{ + for_each_kphotoCollection_member(ptr, &cleanup_member_, NULL); +} + + +// ----------------------------------------------------------------- + +grackle::impl::ChemHeatingRates grackle::impl::new_ChemHeatingRates( + int nelem +) +{ + GRIMPL_REQUIRE(nelem > 0, "nelem must be positive"); + ChemHeatingRates out; + out.n_cr_n = (double*)malloc(sizeof(double)*nelem); + out.n_cr_d1 = (double*)malloc(sizeof(double)*nelem); + out.n_cr_d2 = (double*)malloc(sizeof(double)*nelem); + return out; +} +void grackle::impl::drop_ChemHeatingRates( + grackle::impl::ChemHeatingRates* ptr +) +{ + GRACKLE_FREE(ptr->n_cr_n); + GRACKLE_FREE(ptr->n_cr_d1); + GRACKLE_FREE(ptr->n_cr_d2); +} diff --git a/src/clib/internal_types.hpp b/src/clib/internal_types.hpp index f620b9d0d..11f0d2075 100644 --- a/src/clib/internal_types.hpp +++ b/src/clib/internal_types.hpp @@ -98,6 +98,8 @@ #error "This file must be used by a c++ compiler" #endif +#include "LUT.hpp" + namespace grackle::impl { /// Holds 1D arrays used for cooling and heating @@ -244,6 +246,124 @@ GrainSpeciesCollection new_GrainSpeciesCollection(int nelem); /// This effectively invokes a destructor void drop_GrainSpeciesCollection(GrainSpeciesCollection*); + +/// holds properties about each kind of species +/// +/// The basic premise is that we can use `SpLUT::` to lookup values +/// for the desired species +/// +/// @note +/// The way this is currently a lot like a struct of arrays (i.e. there is an +/// extra level of indirection). +/// - For concreteness, to access index `idx` of the data for HeI, you +/// would write `obj.data[SpLUT::HeI][idx]`. +/// - in terms of performance, this is very similar to what came before (and +/// is good enough for now) +/// - in the long term, it would be even better to replace this with a +/// `View`, but I want to wait until after we have transcribed the +/// bulk of grackle before we do that. +/// +/// @note +/// At the time of writing, both this data structure and the +/// GrainSpeciesCollection data structure both reserve space for each of the +/// dust species. Maybe this data structure shouldn't do this? (If we remove +/// the grain species, maybe we rename this so that it is called +/// ChemSpeciesCollection?) +struct SpeciesCollection { + double* data[SpLUT::NUM_ENTRIES]; +}; + + +/// allocates the contents of a new SpeciesCollection +/// +/// @param nelem The number of elements in each buffer +SpeciesCollection new_SpeciesCollection(int nelem); + +/// performs cleanup of the contents of SpeciesCollection +/// +/// This effectively invokes a destructor +void drop_SpeciesCollection(SpeciesCollection*); + + +/// holds properties about Collisional Reaction Rates that behave in the +/// "standard" way (i.e. we interpolate in 1D with respect to log T) +/// +/// This operates in a similar manner to SpeciesCollection (i.e. we use +/// `CollisionalRxnLUT::` to lookup values for the desired rate). +struct CollisionalRxnRateCollection { + double* data[CollisionalRxnLUT::NUM_ENTRIES]; +}; + +/// allocates the contents of a new CollisionalRxnRateCollection +/// +/// @param nelem The number of elements in each buffer +CollisionalRxnRateCollection new_CollisionalRxnRateCollection(int nelem); + +/// performs cleanup of the contents of CollisionalRxnRateCollection +/// +/// This effectively invokes a destructor +void drop_CollisionalRxnRateCollection(CollisionalRxnRateCollection*); + +/// holds radiative reaction rate buffers +/// +/// @note +/// In the future, it would be nice to use this within the general rate +/// storage struct (and maybe also within the photo_rate_storage struct). +/// Since the storage struct only needs to store these rates as scalars we +/// would need to adopt the LUT strategy (or templates) +struct PhotoRxnRateCollection { + /* Radiative rates for 6-species. */ + double* k24; + double* k25; + double* k26; + + /* Radiative rates for 6-species. */ + double* k27; + double* k28; + double* k29; + double* k30; + double* k31; +}; + +/// allocates the contents of a new PhotoRxnRateCollection +/// +/// @param nelem The number of elements in each buffer +PhotoRxnRateCollection new_PhotoRxnRateCollection(int nelem); + +/// performs cleanup of the contents of PhotoRxnRateCollection +/// +/// This effectively invokes a destructor +void drop_PhotoRxnRateCollection(PhotoRxnRateCollection*); + +/// holds reaction rates chemical heating reaction rates +/// +/// @note +/// In the future, it would be nice to use this within the general rate +/// storage struct +/// +/// @note +/// For now, the modelled rates are fairly limited, but the idea is this could +/// be extended in the future +struct ChemHeatingRates{ + // Chemical heating from H2 formation. + // numerator and denominator of Eq 23 of Omukai ea. 2000. + double *n_cr_n; + double *n_cr_d1; + double *n_cr_d2; +}; + +/// allocates the contents of a new ChemHeatingRates +/// +/// @param nelem The number of elements in each buffer +ChemHeatingRates new_ChemHeatingRates(int nelem); + +/// performs cleanup of the contents of ChemHeatingRates +/// +/// This effectively invokes a destructor +void drop_ChemHeatingRates(ChemHeatingRates*); + + + } // namespace grackle::impl #endif /* INTERNAL_TYPES_HPP */ diff --git a/src/clib/internal_units.h b/src/clib/internal_units.h new file mode 100644 index 000000000..6a1cdfb9f --- /dev/null +++ b/src/clib/internal_units.h @@ -0,0 +1,236 @@ +// See LICENSE file for license and copyright information + +/// @file internal_units.h +/// @brief Declares/Implements the InternalGrUnits struct + +/// This file creates the InternalGrUnits type. +/// - it is fully C-compatible, but behaves in a loosely object-oriented way +/// (i.e. it has associated functions that act a little like methods) +/// - the impetus for creating the constructs in this file is to aide with +/// transcription from Fortran to C++ +/// - you should think of this type (and its associated methods) as a mechansim +/// for grouping together and encapsulating some common logic that was +/// previously repeated across the Fortran versions of the subroutines. +/// +/// As a consequence of this type's role in the transcription process: +/// - the logic was preserved almost **EXACTLY** from the original Fortran files +/// - This is done to maximize consistency and minimize the drift in the results +/// in Grackle's answer tests (modifying the logic could have cascading +/// impacts and plausibly make the results drift outside of tolerances, which +/// is obviously something we want to avoid during transcription) +/// - This means that we may be doing some round-about calculations and using +/// some slightly weird versions of some physical constants +/// +/// **AFTER WE COMPLETE TRANSCRIPTION**, we should simplify this logic +#ifndef INTERNAL_UNITS_HPP +#define INTERNAL_UNITS_HPP + +#ifdef __cplusplus +extern "C" { +#endif /* __cplusplus */ + + +#include + +#include "grackle.h" +#include "grackle_macros.h" +#include "phys_constants.h" + +/// Encapsulates Grackle’s Internal Unit System. +/// +/// Grackle internally employs an Enzo-style comoving coordinate system (and +/// internally tracks various quantities/rates in this unit system). +/// +/// This struct is commonly employed by Grackle’s internal functions that do +/// most of the heavy lifting. We coerce the user-specified units (described by +/// the frontend's `code_units` type) near the top of the Grackle’s internal +/// function call-stack, to an equivalent description in terms of Grackle’s +/// comoving unit system (the latter description encapsulated by this object). +/// +/// @note +/// Once we finish with transcription, we may wish to: +/// - better unify with the code units frontend +/// - adopt more meaningful variable names (most are preserved from the names +/// used in the routines from before transcription) +/// - consider whether or not we even need to track uxyz and urho (I suspect +/// the only reason that these quantities are passed around is that we use +/// them to recompute the other units) +typedef struct InternalGrUnits{ + /// gives the comoving scale-factor when multiplied by a_units + double a_value; + double a_units; + + /// specifies whether user specifies fields in Grackle's comoving coordinates + /// - while this info has historically been passed deeply into the callstack + /// I've confirmed that it's only used for converting the field units + int extfields_in_comoving; + + /// proper-frame: cm per 1 comoving length unit (depends on scale factor, a) + double uxyz; + /// comoving-frame: cm per 1 comoving length unit (independent of a) + double xbase1; + + /// proper frame: g/cm^3 per 1 comoving density unit (depends on a) + double urho; + /// comoving frame: g/cm^3 per 1 comoving density unit (independent of a) + double dbase1; + + /// seconds per time unit (this is defined in our comoving unit system + /// such that it's the same in the proper & comoving frames. In other words, + /// it is independent of a) + double tbase1; + + /// erg*cm^3/s per cooling rate unit (this is defined in our comoving unit + /// system such that it's the same in the proper & comoving frames. In other + /// words, it's independent of a). + /// + /// This is needed in the vast majority of cases where internal units is + /// accessed so we precompute it + double coolunit; + + /// constant coefficients (including unit conversions) that must be + /// multiplied by `internal_energy` (in squared velocity code units) when + /// calculating temperature (this term neglects contributions from mmw and + /// adiabatic index). For more information, see the documentation of + /// `get_temperature_units` in the website docs. + /// + /// This isn't really a "unit", but we are including it in this struct + /// because we historically have called it a unit and it used in the vast + /// majority of cases where this struct gets used. + /// + /// Because of the way that relevant units (namely code velocity) in our + /// comoving unit system are defined, this is the same in the proper & + /// comoving frame + double utem; + +} InternalGrUnits; + +/// Return the version of hydrogen mass constant used by the internal units +static inline double internalu_get_mh_(InternalGrUnits internalu) { + // purely for the sake of consistency, this value of the hydrogen mass is + // initialized as a floating point literal (with the same precision as + // gr_float and then it is casted to a double) + return (double)(mh_grflt); +} + +/// Return the cm*a_unit/s per 1 velocity unit +/// +/// The velocity unit is defined in our comoving unit system such that it's the +/// same in the proper & comoving frames. In other words, it's independent of a +static inline double internalu_get_uvel_(InternalGrUnits internalu) { + return (internalu.uxyz/internalu.a_value) / internalu.tbase1; +} + +/// Return the chunit conversion factor +/// +/// This quantity is used internally to help calculate the amount of heating +/// performed contributed by certain chemical reactions. +static inline double internalu_get_chunit_(InternalGrUnits internalu){ + double mh_local_var = internalu_get_mh_(internalu); + double uvel = internalu_get_uvel_(internalu); + + // before it was moved to this function, there were originally 2 versions of + // this calculation: + // -> the older, commented-out version had an extra factor of 2 in the + // denominator and had a comment stating "1 eV per H2 formed" + // -> the newer version, without the factor of 2 (shown below), had a + // comment stating "1 eV per REACTION (Feb 2020, Gen Chiaki)" + return (1.60218e-12)/(uvel*uvel*mh_local_var); +} + +/// calculates a standard quantity used throughout the codebase +static inline double internalu_calc_dom_(InternalGrUnits internalu) { + const double mh_local_var = internalu_get_mh_(internalu); + return internalu.urho*(pow(internalu.a_value,3))/mh_local_var; +} + +/// calculates coefficients used for computing the Jeans' length +/// +/// The main benefit of factoring this out is that we can be sure that we are +/// using a consistent choices in the versions of constants +static inline double internalu_calc_coef_ljeans_(InternalGrUnits internalu, + double gamma) { + const double mh_local_var = internalu_get_mh_(internalu); + return sqrt((gamma * pi_fortran_val * kboltz_grflt) / + (GravConst_grflt * mh_local_var * internalu.dbase1)); +} + +/// Construct an instance of InternalGrUnits from the frontend_units +static inline InternalGrUnits new_internalu_( + const code_units* frontend_units +) { + // rename the frontend_units object for convenience + const code_units* my_units = frontend_units; + + // Part 1: Common logic extracted from solve_chemistry and + // calculate_(cooling_time|temperature|dust_temperature) that was + // historically used to determine unit-related quantities passed into + // the fortran subroutines + + // calculate temperature units (obviously, don't change this until we finish + // transcription -- under the hood, this calculation implicitly chooses a + // definition of mh that may differ from the value used elesewhere within the + // functions in the current file) + double temperature_units = get_temperature_units(my_units); + + // determine the size of comoving length & comoving density units + // (that are equivalent to the frontend units), as measured in the proper + // reference frame + // - this logic (including variable names) has been preserved from the + // standard pattern for handling units pre-transcription + // - in the future, we can rename some of this + double co_length_units, co_density_units; + if (my_units->comoving_coordinates == TRUE) { + // in this case, the frontend unit-system is exactly the same as grackle's + // internal comoving unit system + co_length_units = my_units->length_units; + co_density_units = my_units->density_units; + } else { + // in this case, the frontend unit-system uses proper units. Here we need + // to do some coercion to properly represent the values using Grackle's + // internal unit system + co_length_units = my_units->length_units * + my_units->a_value * my_units->a_units; + co_density_units = my_units->density_units / + POW(my_units->a_value * my_units->a_units, 3); + } + + // Part 2: initialize output units and copy some stuff from frontend units + InternalGrUnits internalu; + internalu.a_value = my_units->a_value; + internalu.a_units = my_units->a_units; + internalu.extfields_in_comoving = my_units->comoving_coordinates; + + // store the temperature_units (this name remapping would historically occur) + internalu.utem = temperature_units; + + // this represents the name remapping that would occur when passing + // co_length_units and co_density_units into a Fortran routine + internalu.uxyz = co_length_units; + internalu.urho = co_density_units; + + // Part 3: Employ logic from fortran subroutines to compute fundamental unit + // basis of the unit-system and store them in internalu. Then compute + // the cooling units. + + // TODO: (AFTER FINISIHING TRANSCRIPTION) Consider refactoring this logic for + // the case with `my_units->comoving_coordinates==0` (in that case, we do an + // unnecessary round-trip) + internalu.tbase1 = my_units->time_units; + internalu.xbase1 = internalu.uxyz/(my_units->a_value*my_units->a_units); // uxyz is [x]*a = [x]*[a]*a' + internalu.dbase1 = internalu.urho*pow((my_units->a_value*my_units->a_units),3); // urho is [dens]/a^3 = [dens]/([a]*a')^3 ' + + // lastly, compute coolunit (make sure we use the correct version of mh) + const double mh_local_var = internalu_get_mh_(internalu); + internalu.coolunit = ( + pow(my_units->a_units,5) * pow(internalu.xbase1,2) * pow(mh_local_var,2) + ) / (pow(internalu.tbase1,3) * internalu.dbase1); + + return internalu; +} + +#ifdef __cplusplus +} // extern "C" +#endif /* __cplusplus */ + +#endif /* INTERNAL_UNITS_HPP */ diff --git a/src/clib/solve_chemistry.c b/src/clib/solve_chemistry.c index 6856d109e..130293bb8 100644 --- a/src/clib/solve_chemistry.c +++ b/src/clib/solve_chemistry.c @@ -13,157 +13,20 @@ #include #include +#include "grackle.h" #include "grackle_macros.h" -#include "grackle_types.h" -#include "grackle_chemistry_data.h" +#include "internal_units.h" #include "phys_constants.h" +#include "solve_rate_cool_g-cpp.h" #include "utils.h" -extern chemistry_data *grackle_data; -extern chemistry_data_storage grackle_rates; - /* function prototypes */ -double get_temperature_units(code_units *my_units); - int update_UVbackground_rates(chemistry_data *my_chemistry, chemistry_data_storage *my_rates, photo_rate_storage *my_uvb_rates, code_units *my_units); -extern void FORTRAN_NAME(solve_rate_cool_g)( - int *icool, - gr_float *d, gr_float *e, gr_float *u, gr_float *v, gr_float *w, gr_float *de, - gr_float *HI, gr_float *HII, gr_float *HeI, gr_float *HeII, gr_float *HeIII, - int *in, int *jn, int *kn, int *nratec, int *iexpand, - int *ispecies, int *imetal, int *imcool, int *idust, int *idustall, - int *idustfield, int *idim, - int *is, int *js, int *ks, int *ie, int *je, int *ke, - int *ih2co, int *ipiht, int *idustrec, int *igammah, - double *dx, double *dt, double *aye, double *temstart, double *temend, - double *utem, double *uxyz, double *uaye, double *urho, double *utim, - double *gamma, double *fh, double *dtoh, double *z_solar, double *fgr, - double *k1a, double *k2a, double *k3a, double *k4a, double *k5a, - double *k6a, double *k7a, double *k8a, double *k9a, double *k10a, - double *k11a, double *k12a, double *k13a, double *k13dda, double *k14a, - double *k15a, double *k16a, double *k17a, double *k18a, double *k19a, - double *k22a, double *k24, double *k25, double *k26, double *k27, - double *k28, double *k29, double *k30, double *k31, - double *k50a, double *k51a, double *k52a, double *k53a, double *k54a, - double *k55a, double *k56a, double *k57a, double *k58a, - int *ndratec, double *dtemstart, double *dtemend, double *h2dusta, - double *ncrna, double *ncrd1a, double *ncrd2a, - double *ceHIa, double *ceHeIa, double *ceHeIIa, double *ciHIa, - double *ciHeIa, double *ciHeISa, double *ciHeIIa, - double *reHIIa, double *reHeII1a, double *reHeII2a, double *reHeIIIa, - double *brema, double *compa, double *gammaha, double *isrf, - double *regra, double *gamma_isrfa, double *comp_xraya, double *comp_temp, - double *piHI, double *piHeI, double *piHeII, - gr_float *HM, gr_float *H2I, gr_float *H2II, - gr_float *DI, gr_float *DII, gr_float *HDI, - gr_float *metal, gr_float *dust, - double *hyd01ka, double *h2k01a, double *vibha, - double *rotha, double *rotla, - double *gpldl, double *gphdl, double *HDltea, double *HDlowa, - double *gaHIa, double *gaH2a, double *gaHea, double *gaHpa, double *gaela, - double *h2ltea, double *gasgra, int *iH2shield, - int *iradshield, double *avgsighi, double *avgsighei, double *avgsigheii, - int *iradtrans, int *iradcoupled, int *iradstep, int *irt_honly, - gr_float *kphHI, gr_float *kphHeI, gr_float *kphHeII, gr_float *kdissH2I, - gr_float *photogamma, gr_float *xH2shield, - int *ierr, - int *ih2optical, int *iciecool, int *ithreebody, int *ih2cr, int *ihdcr, - double *ciecoa, - int *icmbTfloor, int *iClHeat, double *clEleFra, - long long *priGridRank, long long *priGridDim, - double *priPar1, double *priPar2, double *priPar3, - double *priPar4, double *priPar5, - long long *priDataSize, double *priCooling, - double *priHeating, double *priMMW, - long long *metGridRank, long long *metGridDim, - double *metPar1, double *metPar2, double *metPar3, - double *metPar4, double *metPar5, - long long *metDataSize, double *metCooling, - double *metHeating, int *clnew, - int *iVheat, int *iMheat, gr_float *Vheat, gr_float *Mheat, - int *iTfloor, gr_float *Tfloor_scalar, gr_float *Tfloor, - int *imchem, int *igrgr, int *ipcont, double *tmcool, - gr_float *DM, gr_float *HDII, gr_float *HeHII, - gr_float *CI, gr_float *CII, gr_float *CO, gr_float *CO2, - gr_float *OI, gr_float *OH, gr_float *H2O, gr_float *O2, - gr_float *SiI, gr_float *SiOI, gr_float *SiO2I, - gr_float *CH, gr_float *CH2, gr_float *COII, gr_float *OII, - gr_float *OHII, gr_float *H2OII, gr_float *H3OII, gr_float *O2II, - gr_float *Mg, gr_float *Al, gr_float *S, gr_float *Fe, - gr_float *SiM, gr_float *FeM, gr_float *Mg2SiO4, gr_float *MgSiO3, gr_float *Fe3O4, - gr_float *AC, gr_float *SiO2D, gr_float *MgO, gr_float *FeS, gr_float *Al2O3, - gr_float *reforg, gr_float *volorg, gr_float *H2Oice, - double *k125a, double *k129a, double *k130a, double *k131a, double *k132a, - double *k133a, double *k134a, double *k135a, double *k136a, double *k137a, - double *k148a, double *k149a, double *k150a, double *k151a, double *k152a, - double *k153a, - double *kz15a, double *kz16a, double *kz17a, double *kz18a, double *kz19a, - double *kz20a, double *kz21a, double *kz22a, double *kz23a, double *kz24a, - double *kz25a, double *kz26a, double *kz27a, double *kz28a, double *kz29a, - double *kz30a, double *kz31a, double *kz32a, double *kz33a, double *kz34a, - double *kz35a, double *kz36a, double *kz37a, double *kz38a, double *kz39a, - double *kz40a, double *kz41a, double *kz42a, double *kz43a, double *kz44a, - double *kz45a, double *kz46a, double *kz47a, double *kz48a, double *kz49a, - double *kz50a, double *kz51a, double *kz52a, double *kz53a, double *kz54a, - double *cieY06, - long long *LH2_N, long long *LH2_Size, - double *LH2_D, double *LH2_T, double *LH2_H, - double *LH2_dD, double *LH2_dT, double *LH2_dH, double *LH2_L, - long long *LHD_N, long long *LHD_Size, - double *LHD_D, double *LHD_T, double *LHD_H, - double *LHD_dD, double *LHD_dT, double *LHD_dH, double *LHD_L, - long long *LCI_N, long long *LCI_Size, - double *LCI_D, double *LCI_T, double *LCI_H, - double *LCI_dD, double *LCI_dT, double *LCI_dH, double *LCI_L, - long long *LCII_N, long long *LCII_Size, - double *LCII_D, double *LCII_T, double *LCII_H, - double *LCII_dD, double *LCII_dT, double *LCII_dH, double *LCII_L, - long long *LOI_N, long long *LOI_Size, - double *LOI_D, double *LOI_T, double *LOI_H, - double *LOI_dD, double *LOI_dT, double *LOI_dH, double *LOI_L, - long long *LCO_N, long long *LCO_Size, - double *LCO_D, double *LCO_T, double *LCO_H, - double *LCO_dD, double *LCO_dT, double *LCO_dH, double *LCO_L, - long long *LOH_N, long long *LOH_Size, - double *LOH_D, double *LOH_T, double *LOH_H, - double *LOH_dD, double *LOH_dT, double *LOH_dH, double *LOH_L, - long long *LH2O_N, long long *LH2O_Size, - double *LH2O_D, double *LH2O_T, double *LH2O_H, - double *LH2O_dD, double *LH2O_dT, double *LH2O_dH, double *LH2O_L, - long long *alphap_N, long long *alphap_Size, - double *alphap_D, double *alphap_T, double *alphap_dD, double *alphap_dT, - double *alphap_Data, - int *immulti, int *imabund, int *idspecies, int *itdmulti, int *idsub, - gr_float *metal_loc, gr_float *metal_C13, gr_float *metal_C20, gr_float *metal_C25, gr_float *metal_C30, - gr_float *metal_F13, gr_float *metal_F15, gr_float *metal_F50, gr_float *metal_F80, - gr_float *metal_P170, gr_float *metal_P200, gr_float *metal_Y19, - int *SN0_N, - double *SN0_XC , double *SN0_XO, double *SN0_XMg, double *SN0_XAl, - double *SN0_XSi, double *SN0_XS, double *SN0_XFe, - double *SN0_fC , double *SN0_fO, double *SN0_fMg, double *SN0_fAl, - double *SN0_fSi, double *SN0_fS, double *SN0_fFe, - double *SN0_fSiM, double *SN0_fFeM, double *SN0_fMg2SiO4, double *SN0_fMgSiO3, double *SN0_fFe3O4, - double *SN0_fAC, double *SN0_fSiO2D, double *SN0_fMgO, double *SN0_fFeS, double *SN0_fAl2O3, - double *SN0_freforg , double *SN0_fvolorg , double *SN0_fH2Oice, - double *SN0_r0SiM, double *SN0_r0FeM, double *SN0_r0Mg2SiO4, double *SN0_r0MgSiO3, double *SN0_r0Fe3O4, - double *SN0_r0AC, double *SN0_r0SiO2D, double *SN0_r0MgO, double *SN0_r0FeS, double *SN0_r0Al2O3, - double *SN0_r0reforg , double *SN0_r0volorg , double *SN0_r0H2Oice, - int *gr_N, int *gr_Size, double *gr_dT, double *gr_Td, - double *SN0_kpSiM, double *SN0_kpFeM, double *SN0_kpMg2SiO4, double *SN0_kpMgSiO3, double *SN0_kpFe3O4, - double *SN0_kpAC, double *SN0_kpSiO2D, double *SN0_kpMgO, double *SN0_kpFeS, double *SN0_kpAl2O3, - double *SN0_kpreforg , double *SN0_kpvolorg , double *SN0_kpH2Oice, - double *h2dustSa, double *h2dustCa, double *gasgr2a, double *gamma_isrf2a, double *grogra, - int *idissHDI, gr_float *kdissHDI, int *iionZ, gr_float *kphCI, gr_float *kphOI, - int *idissZ, gr_float *kdissCO, gr_float *kdissOH, gr_float *kdissH2O, int *iuseH2shield, - int *iisrffield, gr_float* isrf_habing, - int *iH2shieldcustom, gr_float* f_shield_custom, - int *itmax, int *exititmax); - int local_solve_chemistry(chemistry_data *my_chemistry, chemistry_data_storage *my_rates, code_units *my_units, @@ -218,17 +81,7 @@ int local_solve_chemistry(chemistry_data *my_chemistry, if (my_fields->metal_density == NULL) metal_field_present = FALSE; - double co_length_units, co_density_units; - if (my_units->comoving_coordinates == TRUE) { - co_length_units = my_units->length_units; - co_density_units = my_units->density_units; - } - else { - co_length_units = my_units->length_units * - my_units->a_value * my_units->a_units; - co_density_units = my_units->density_units / - POW(my_units->a_value * my_units->a_units, 3); - } + InternalGrUnits internalu = new_internalu_(my_units); /* Error checking for H2 shielding approximation */ if (self_shielding_err_check(my_chemistry, my_fields, @@ -236,486 +89,14 @@ int local_solve_chemistry(chemistry_data *my_chemistry, return FAIL; } - /* Calculate temperature units. */ - - double temperature_units = get_temperature_units(my_units); - - /* Call the fortran routine to solve cooling equations. */ - - int ierr = 0; + /* Call the routine to solve cooling equations. */ - FORTRAN_NAME(solve_rate_cool_g)( - &my_chemistry->with_radiative_cooling, - my_fields->density, - my_fields->internal_energy, - my_fields->x_velocity, - my_fields->y_velocity, - my_fields->z_velocity, - 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->grid_dimension+0, - my_fields->grid_dimension+1, - my_fields->grid_dimension+2, - &my_chemistry->NumberOfTemperatureBins, - &my_units->comoving_coordinates, - &my_chemistry->primordial_chemistry, - &metal_field_present, - &my_chemistry->metal_cooling, - &my_chemistry->h2_on_dust, - &my_chemistry->dust_chemistry, - &my_chemistry->use_dust_density_field, - &(my_fields->grid_rank), - my_fields->grid_start+0, - my_fields->grid_start+1, - my_fields->grid_start+2, - my_fields->grid_end+0, - my_fields->grid_end+1, - my_fields->grid_end+2, - &my_chemistry->ih2co, - &my_chemistry->ipiht, - &my_chemistry->dust_recombination_cooling, - &my_chemistry->photoelectric_heating, - &(my_fields->grid_dx), - &dt_value, - &my_units->a_value, - &my_chemistry->TemperatureStart, - &my_chemistry->TemperatureEnd, - &temperature_units, - &co_length_units, - &my_units->a_units, - &co_density_units, - &my_units->time_units, - &my_chemistry->Gamma, - &my_chemistry->HydrogenFractionByMass, - &my_chemistry->DeuteriumToHydrogenRatio, - &my_chemistry->SolarMetalFractionByMass, - &my_chemistry->local_dust_to_gas_ratio, - my_rates->k1, - my_rates->k2, - my_rates->k3, - my_rates->k4, - my_rates->k5, - my_rates->k6, - my_rates->k7, - my_rates->k8, - my_rates->k9, - my_rates->k10, - my_rates->k11, - my_rates->k12, - my_rates->k13, - my_rates->k13dd, - my_rates->k14, - my_rates->k15, - my_rates->k16, - my_rates->k17, - my_rates->k18, - my_rates->k19, - my_rates->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, - &my_uvb_rates.k31, - my_rates->k50, - my_rates->k51, - my_rates->k52, - my_rates->k53, - my_rates->k54, - my_rates->k55, - my_rates->k56, - my_rates->k57, - my_rates->k58, - &my_chemistry->NumberOfDustTemperatureBins, - &my_chemistry->DustTemperatureStart, - &my_chemistry->DustTemperatureEnd, - my_rates->h2dust, - my_rates->n_cr_n, - my_rates->n_cr_d1, - my_rates->n_cr_d2, - my_rates->ceHI, - my_rates->ceHeI, - my_rates->ceHeII, - my_rates->ciHI, - my_rates->ciHeI, - my_rates->ciHeIS, - my_rates->ciHeII, - my_rates->reHII, - my_rates->reHeII1, - my_rates->reHeII2, - my_rates->reHeIII, - my_rates->brem, - &my_rates->comp, - &my_rates->gammah, - &my_chemistry->interstellar_radiation_field, - my_rates->regr, - &my_rates->gamma_isrf, - &my_uvb_rates.comp_xray, - &my_uvb_rates.temp_xray, - &my_uvb_rates.piHI, - &my_uvb_rates.piHeI, - &my_uvb_rates.piHeII, - my_fields->HM_density, - my_fields->H2I_density, - my_fields->H2II_density, - my_fields->DI_density, - my_fields->DII_density, - my_fields->HDI_density, - my_fields->metal_density, - my_fields->dust_density, - my_rates->hyd01k, - my_rates->h2k01, - my_rates->vibh, - my_rates->roth, - my_rates->rotl, - my_rates->GP99LowDensityLimit, - my_rates->GP99HighDensityLimit, - my_rates->HDlte, - my_rates->HDlow, - my_rates->GAHI, - my_rates->GAH2, - my_rates->GAHe, - my_rates->GAHp, - my_rates->GAel, - my_rates->H2LTE, - my_rates->gas_grain, - &my_chemistry->H2_self_shielding, - &my_chemistry->self_shielding_method, - &my_uvb_rates.crsHI, - &my_uvb_rates.crsHeI, - &my_uvb_rates.crsHeII, - &my_chemistry->use_radiative_transfer, - &my_chemistry->radiative_transfer_coupled_rate_solver, - &my_chemistry->radiative_transfer_intermediate_step, - &my_chemistry->radiative_transfer_hydrogen_only, - my_fields->RT_HI_ionization_rate, - my_fields->RT_HeI_ionization_rate, - my_fields->RT_HeII_ionization_rate, - my_fields->RT_H2_dissociation_rate, - my_fields->RT_heating_rate, - my_fields-> H2_self_shielding_length, - &ierr, - &my_chemistry->h2_optical_depth_approximation, - &my_chemistry->cie_cooling, - &my_chemistry->three_body_rate, - &my_chemistry->h2_cooling_rate, - &my_chemistry->hd_cooling_rate, - my_rates->cieco, - &my_chemistry->cmb_temperature_floor, - &my_chemistry->UVbackground, - &my_chemistry->cloudy_electron_fraction_factor, - &my_rates->cloudy_primordial.grid_rank, - my_rates->cloudy_primordial.grid_dimension, - my_rates->cloudy_primordial.grid_parameters[0], - my_rates->cloudy_primordial.grid_parameters[1], - my_rates->cloudy_primordial.grid_parameters[2], - my_rates->cloudy_primordial.grid_parameters[3], - my_rates->cloudy_primordial.grid_parameters[4], - &my_rates->cloudy_primordial.data_size, - my_rates->cloudy_primordial.cooling_data, - my_rates->cloudy_primordial.heating_data, - my_rates->cloudy_primordial.mmw_data, - &my_rates->cloudy_metal.grid_rank, - my_rates->cloudy_metal.grid_dimension, - my_rates->cloudy_metal.grid_parameters[0], - my_rates->cloudy_metal.grid_parameters[1], - my_rates->cloudy_metal.grid_parameters[2], - my_rates->cloudy_metal.grid_parameters[3], - my_rates->cloudy_metal.grid_parameters[4], - &my_rates->cloudy_metal.data_size, - my_rates->cloudy_metal.cooling_data, - my_rates->cloudy_metal.heating_data, - &my_rates->cloudy_data_new, - &my_chemistry->use_volumetric_heating_rate, - &my_chemistry->use_specific_heating_rate, - my_fields->volumetric_heating_rate, - my_fields->specific_heating_rate, - &my_chemistry->use_temperature_floor, - &my_chemistry->temperature_floor_scalar, - my_fields->temperature_floor, - &my_chemistry->metal_chemistry, - &my_chemistry->grain_growth, - &my_chemistry->use_primordial_continuum_opacity, - &my_chemistry->tabulated_cooling_minimum_temperature, - my_fields->DM_density, - my_fields->HDII_density, - my_fields->HeHII_density, - 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, - my_rates->k125, - my_rates->k129, - my_rates->k130, - my_rates->k131, - my_rates->k132, - my_rates->k133, - my_rates->k134, - my_rates->k135, - my_rates->k136, - my_rates->k137, - my_rates->k148, - my_rates->k149, - my_rates->k150, - my_rates->k151, - my_rates->k152, - my_rates->k153, - my_rates->kz15, - my_rates->kz16, - my_rates->kz17, - my_rates->kz18, - my_rates->kz19, - my_rates->kz20, - my_rates->kz21, - my_rates->kz22, - my_rates->kz23, - my_rates->kz24, - my_rates->kz25, - my_rates->kz26, - my_rates->kz27, - my_rates->kz28, - my_rates->kz29, - my_rates->kz30, - my_rates->kz31, - my_rates->kz32, - my_rates->kz33, - my_rates->kz34, - my_rates->kz35, - my_rates->kz36, - my_rates->kz37, - my_rates->kz38, - my_rates->kz39, - my_rates->kz40, - my_rates->kz41, - my_rates->kz42, - my_rates->kz43, - my_rates->kz44, - my_rates->kz45, - my_rates->kz46, - my_rates->kz47, - my_rates->kz48, - my_rates->kz49, - my_rates->kz50, - my_rates->kz51, - my_rates->kz52, - my_rates->kz53, - my_rates->kz54, - my_rates->cieY06, - my_rates->LH2.props.dimension, - &my_rates->LH2.props.data_size, - my_rates->LH2.props.parameters[0], - my_rates->LH2.props.parameters[1], - my_rates->LH2.props.parameters[2], - &my_rates->LH2.props.parameter_spacing[0], - &my_rates->LH2.props.parameter_spacing[1], - &my_rates->LH2.props.parameter_spacing[2], - my_rates->LH2.data, - my_rates->LHD.props.dimension, - &my_rates->LHD.props.data_size, - my_rates->LHD.props.parameters[0], - my_rates->LHD.props.parameters[1], - my_rates->LHD.props.parameters[2], - &my_rates->LHD.props.parameter_spacing[0], - &my_rates->LHD.props.parameter_spacing[1], - &my_rates->LHD.props.parameter_spacing[2], - my_rates->LHD.data, - my_rates->LCI.props.dimension, - &my_rates->LCI.props.data_size, - my_rates->LCI.props.parameters[0], - my_rates->LCI.props.parameters[1], - my_rates->LCI.props.parameters[2], - &my_rates->LCI.props.parameter_spacing[0], - &my_rates->LCI.props.parameter_spacing[1], - &my_rates->LCI.props.parameter_spacing[2], - my_rates->LCI.data, - my_rates->LCII.props.dimension, - &my_rates->LCII.props.data_size, - my_rates->LCII.props.parameters[0], - my_rates->LCII.props.parameters[1], - my_rates->LCII.props.parameters[2], - &my_rates->LCII.props.parameter_spacing[0], - &my_rates->LCII.props.parameter_spacing[1], - &my_rates->LCII.props.parameter_spacing[2], - my_rates->LCII.data, - my_rates->LOI.props.dimension, - &my_rates->LOI.props.data_size, - my_rates->LOI.props.parameters[0], - my_rates->LOI.props.parameters[1], - my_rates->LOI.props.parameters[2], - &my_rates->LOI.props.parameter_spacing[0], - &my_rates->LOI.props.parameter_spacing[1], - &my_rates->LOI.props.parameter_spacing[2], - my_rates->LOI.data, - my_rates->LCO.props.dimension, - &my_rates->LCO.props.data_size, - my_rates->LCO.props.parameters[0], - my_rates->LCO.props.parameters[1], - my_rates->LCO.props.parameters[2], - &my_rates->LCO.props.parameter_spacing[0], - &my_rates->LCO.props.parameter_spacing[1], - &my_rates->LCO.props.parameter_spacing[2], - my_rates->LCO.data, - my_rates->LOH.props.dimension, - &my_rates->LOH.props.data_size, - my_rates->LOH.props.parameters[0], - my_rates->LOH.props.parameters[1], - my_rates->LOH.props.parameters[2], - &my_rates->LOH.props.parameter_spacing[0], - &my_rates->LOH.props.parameter_spacing[1], - &my_rates->LOH.props.parameter_spacing[2], - my_rates->LOH.data, - my_rates->LH2O.props.dimension, - &my_rates->LH2O.props.data_size, - my_rates->LH2O.props.parameters[0], - my_rates->LH2O.props.parameters[1], - my_rates->LH2O.props.parameters[2], - &my_rates->LH2O.props.parameter_spacing[0], - &my_rates->LH2O.props.parameter_spacing[1], - &my_rates->LH2O.props.parameter_spacing[2], - my_rates->LH2O.data, - my_rates->alphap.props.dimension, - &my_rates->alphap.props.data_size, - my_rates->alphap.props.parameters[0], - my_rates->alphap.props.parameters[1], - &my_rates->alphap.props.parameter_spacing[0], - &my_rates->alphap.props.parameter_spacing[1], - my_rates->alphap.data, - &my_chemistry->multi_metals, - &my_chemistry->metal_abundances, - &my_chemistry->dust_species, - &my_chemistry->use_multiple_dust_temperatures, - &my_chemistry->dust_sublimation, - my_fields->local_ISM_metal_density, - my_fields->ccsn13_metal_density, - my_fields->ccsn20_metal_density, - my_fields->ccsn25_metal_density, - my_fields->ccsn30_metal_density, - my_fields->fsn13_metal_density, - my_fields->fsn15_metal_density, - my_fields->fsn50_metal_density, - my_fields->fsn80_metal_density, - my_fields->pisn170_metal_density, - my_fields->pisn200_metal_density, - my_fields->y19_metal_density, - &my_rates->SN0_N, - my_rates->SN0_XC, - my_rates->SN0_XO, - my_rates->SN0_XMg, - my_rates->SN0_XAl, - my_rates->SN0_XSi, - my_rates->SN0_XS, - my_rates->SN0_XFe, - my_rates->SN0_fC, - my_rates->SN0_fO, - my_rates->SN0_fMg, - my_rates->SN0_fAl, - my_rates->SN0_fSi, - my_rates->SN0_fS, - my_rates->SN0_fFe, - my_rates->SN0_fSiM, - my_rates->SN0_fFeM, - my_rates->SN0_fMg2SiO4, - my_rates->SN0_fMgSiO3, - my_rates->SN0_fFe3O4, - my_rates->SN0_fAC, - my_rates->SN0_fSiO2D, - my_rates->SN0_fMgO, - my_rates->SN0_fFeS, - my_rates->SN0_fAl2O3, - my_rates->SN0_freforg, - my_rates->SN0_fvolorg, - my_rates->SN0_fH2Oice, - my_rates->SN0_r0SiM, - my_rates->SN0_r0FeM, - my_rates->SN0_r0Mg2SiO4, - my_rates->SN0_r0MgSiO3, - my_rates->SN0_r0Fe3O4, - my_rates->SN0_r0AC, - my_rates->SN0_r0SiO2D, - my_rates->SN0_r0MgO, - my_rates->SN0_r0FeS, - my_rates->SN0_r0Al2O3, - my_rates->SN0_r0reforg, - my_rates->SN0_r0volorg, - my_rates->SN0_r0H2Oice, - my_rates->gr_N, - &my_rates->gr_Size, - &my_rates->gr_dT, - my_rates->gr_Td, - my_rates->SN0_kpSiM, - my_rates->SN0_kpFeM, - my_rates->SN0_kpMg2SiO4, - my_rates->SN0_kpMgSiO3, - my_rates->SN0_kpFe3O4, - my_rates->SN0_kpAC, - my_rates->SN0_kpSiO2D, - my_rates->SN0_kpMgO, - my_rates->SN0_kpFeS, - my_rates->SN0_kpAl2O3, - my_rates->SN0_kpreforg, - my_rates->SN0_kpvolorg, - my_rates->SN0_kpH2Oice, - my_rates->h2dustS, - my_rates->h2dustC, - my_rates->gas_grain2, - &my_rates->gamma_isrf2, - my_rates->grain_growth_rate, - &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, - &my_chemistry->radiative_transfer_use_H2_shielding, - &my_chemistry->use_isrf_field, - my_fields->isrf_habing, - &my_chemistry->H2_custom_shielding, - my_fields->H2_custom_shielding_factor, - &my_chemistry->max_iterations, - &my_chemistry->exit_after_iterations_exceeded); + int ierr = solve_rate_cool_g( + metal_field_present, dt_value, internalu, + my_chemistry, my_rates, my_fields, &my_uvb_rates + ); - if (ierr == FAIL) { + if (ierr == GR_FAIL) { fprintf(stderr, "Error in solve_rate_cool_g.\n"); } diff --git a/src/clib/solve_rate_cool_g-cpp.C b/src/clib/solve_rate_cool_g-cpp.C new file mode 100644 index 000000000..7e02b076e --- /dev/null +++ b/src/clib/solve_rate_cool_g-cpp.C @@ -0,0 +1,1012 @@ +// See LICENSE file for license and copyright information + +/// @file solve_rate_cool_g-cpp.C +/// @brief Declares signature of solve_rate_cool_g + +// This file was initially generated automatically during conversion of the +// solve_rate_cool_g function from FORTRAN to C++ + +#include +#include // std::malloc, std::free +#include // std::memcpy +#include + +#include "grackle.h" +#include "fortran_func_wrappers.hpp" +#include "index_helper.h" +#include "internal_types.hpp" +#include "internal_units.h" +#include "utils-cpp.hpp" + +#include "solve_rate_cool_g-cpp.h" + +/// overrides the subcycle timestep (for each index in the index-range that is +/// selected by the given itmask) with the maximum allowed heating/cooling +/// timestep when the current value is larger. +/// +/// @param[out] dtit buffer tracking the current subcycle timestep for each +/// index in the index-range. If the current value exceeds the maximum +/// allowed heating/cooling timestep, the values will overwritten +/// @param[in] idx_range Specifies the current index-range +/// @param[in] dt tracks the full timestep that all the subcycles will +/// eventually add up to +/// @param[in] ttot tracks the total time that has already elapsed from +/// previous subcycles for each location in `idx_range` +/// @param[in] itmask Specifies the `idx_range`'s iteration-mask for this +/// calculation +/// @param[in] tgas specifies the gas temperatures for the `idx_range` +/// @param[in] p2d specifies the pressures for the `idx_range`. This is +/// computed user-specified nominal adiabatic index value (i.e. no attempts +/// are made to correct for presence of H2) +/// @param[in,out] edot specifies the time derivative of internal energy +/// density for each location in `idx_range`. This may be overwritten to +/// enforce the floor. +static void enforce_max_heatcool_subcycle_dt_( + double* dtit, IndexRange idx_range, double dt, const double* ttot, + const gr_mask_type* itmask, const double* tgas, const double* p2d, + double* edot, const chemistry_data* my_chemistry +) { + + for (int i = idx_range.i_start; i < idx_range.i_stop; i++) { + + if (itmask[i] != MASK_FALSE) { + // Set energy per unit volume of this cell based in the pressure + // (the gamma used here is the right one even for H2 since p2d + // is calculated with this gamma). + double energy = std::fmax(p2d[i]/(my_chemistry->Gamma-1.), tiny8); + + // If the temperature is at the bottom of the temperature look-up + // table and edot < 0, then shut off the cooling. + if (tgas[i] <= 1.01*my_chemistry->TemperatureStart && edot[i] < 0.) { + edot[i] = tiny8; + } + + // enforce the floor + if (std::fabs(edot[i]) < tiny8) { edot[i] = tiny8; } + + // Compute timestep for 10% change + dtit[i] = grackle::impl::fmin( + (double)(std::fabs(0.1 * energy / edot[i])), dt - ttot[i], dtit[i] + ); + + if (dtit[i] != dtit[i]) { + OMP_PRAGMA_CRITICAL + { + eprintf("HUGE dtit :: %g %g %g %g %g %g %g\n", + energy, edot[i], dtit[i], dt, ttot[i], + std::fabs(0.1 * energy / edot[i]), + (double)(std::fabs(0.1 * energy / edot [i])) ); + } + } + + } + } + +} + +// ------------------------------------------------------------- + +/// Selects the scheme that will be used to evolve the chemistry network and +/// then accordingly updates the iteration masks (and `imp_eng`) +/// +/// There are 2 schemes: +/// 1. Gauss-Seidel (for low-density zones) +/// 2. Newton-Raphson (for high-density zones) in 2 modes: +/// - internal energy evolution is operator-split (handled separately from +/// chemistry network) +/// - internal energy is coupled with the rest of the chemistry network +/// +/// @param[in] idx_range Specifies the current index-range +/// @param[in,out] itmask Initially specifies all locations to be evolved +/// during the current subcycle (in `idx_range`). Will be updated to only +/// specify the locations to apply Gauss-Seidel scheme +/// @param[out] itmask_nr Buffer for `idx_range` that is used to specify +/// locations where we will apply Newton-Raphson scheme +/// @param[out] imp_eng Buffer for `idx_range` where the choice of +/// energy-evolution handling is recorded for the Newton-Raphson scheme +/// @param[out] itmask_tmp Buffer where the initial values of itmask are +/// copied into +/// @param[in] mask_len the length of the iteration masks +/// @param[in] imetal specifies whether or not the caller provided a metal +/// density field +/// @param[in] min_metallicity specifies the minimum metallicity where we +/// consider metal chemistry/cooling +/// @param[in] ddom specifies precomputed product of mass density and the +/// `dom` quantity for each location in `idx_range` +/// @param[in] tgas specifies the gas temperatures for the `idx_range` +/// @param[in] metallicity specifies the metallicity for the `idx_range` +/// +/// @todo +/// It might make more sense to create `itmask_tmp` before calling this +/// function and then completely ignore the initial values in `itmask` (this +/// would be far less confusing) +static void select_chem_scheme_update_masks_( + IndexRange idx_range, gr_mask_type* itmask, gr_mask_type* itmask_nr, + int* imp_eng, gr_mask_type* itmask_tmp, int mask_len, int imetal, + double min_metallicity, const double* ddom, const double* tgas, + const double* metallicity, const chemistry_data* my_chemistry +) { + + std::memcpy(itmask_tmp, itmask, sizeof(gr_mask_type)*mask_len); + std::memcpy(itmask_nr, itmask, sizeof(gr_mask_type)*mask_len); + + // would it be more robust to use my_chemistry->metal_cooling than imetal? + + for (int i = idx_range.i_start; i < idx_range.i_stop; i++) { + if ( itmask_tmp[i] != MASK_FALSE ) { + + if ( (imetal == 0) && (ddom[i] < 1.e8) ) { + itmask_nr[i] = MASK_FALSE; + } else if ( + (imetal == 1) && + (((metallicity[i] <= min_metallicity) && (ddom[i] < 1.e8)) || + ((metallicity[i] > min_metallicity) && (ddom[i] < 1.e6))) + ) { + itmask_nr[i] = MASK_FALSE; + } else { + itmask[i] = MASK_FALSE; + } + + } + } + + for (int i = idx_range.i_start; i < idx_range.i_stop; i++) { + if (itmask_nr[i] != MASK_FALSE) { + if ((my_chemistry->with_radiative_cooling == 1) && + (my_chemistry->primordial_chemistry > 1) && + (ddom[i] > 1.e7) && (tgas[i] > 1650.0)) { + imp_eng[i] = 1; + } else { + imp_eng[i] = 0; + } + } + } +} + +// ------------------------------------------------------------- + +/// Computes the timescale given by `ndens_Heq / (d ndens_Heq / d t)` +/// +/// This is used to help compute the subcycle timestep when using a primordial +/// chemistry solver +/// +/// @param my_chemistry holds a number of configuration parameters +/// @param my_rates holds assorted rate data. In this function, this is being +/// used to specify some the interpolation tables of some relevant reaction +/// rates (they are tabulated with respect to logT) +/// @param dlogtem Specifies the constant spacing shared by the relevant rate +/// interpolation tables +/// @param logTlininterp_buf Specifies the information related to the position +/// in the logT interpolations (for a number of chemistry zones) +/// @param k13, k22 1D arrays specifying the previously looked up, local values +/// of the k13 and k22 rates. +/// @param local_rho specifies the local (total) mass density +/// @param tgas 1D array specifying the temperature +/// @param p2d 1D array specifying the pressures values. This is computed from +/// the user-specified nominal adiabatic index value (i.e. no attempts +/// are made to correct for presence of H2) +/// @param edot 1D array specifying the time derivative of the internal energy +/// density +/// @param i Specifies the index of the relevant zone in the 1D array. (**BE +/// AWARE:** this is a 0-based index) +/// +/// @note +/// The `static` annotation indicates that this function is only visible to the +/// current translation unit +static double calc_Heq_div_dHeqdt_( + const chemistry_data* my_chemistry, + const chemistry_data_storage* my_rates, + double dlogtem, + const grackle::impl::LogTLinInterpScratchBuf logTlininterp_buf, + const double* k13, + const double* k22, + double local_rho, + const double* tgas, + const double* p2d, + const double* edot, + int i +) { + + // Equilibrium value for H is: + // Heq = (-1._DKIND / (4*k22)) * (k13 - sqrt(8 k13 k22 rho + k13^2)) + // We want to know dH_eq/dt. + // - We can trivially get dH_eq/dT. + // - We have de/dt. + // - We need dT/de. + // + // T = (g-1)*p2d*utem/N; tgas == (g-1)(p2d*utem/N) + // dH_eq / dt = (dH_eq/dT) * (dT/de) * (de/dt) + // dH_eq / dT (see above; we can calculate the derivative here) + // dT / de = utem * (gamma - 1._DKIND) / N == tgas / p2d + // de / dt = edot + // Now we use our estimate of dT/de to get the estimated + // difference in the equilibrium + double eqt2 = std::fmin(std::log(tgas[i]) + 0.1*dlogtem, logTlininterp_buf.t2[i]); + double eqtdef = (eqt2 - logTlininterp_buf.t1[i])/(logTlininterp_buf.t2[i] - logTlininterp_buf.t1[i]); + double eqk222 = my_rates->k22[logTlininterp_buf.indixe[i]-1] + + (my_rates->k22[logTlininterp_buf.indixe[i]+1-1] -my_rates->k22[logTlininterp_buf.indixe[i]-1])*eqtdef; + double eqk132 = my_rates->k13[logTlininterp_buf.indixe[i]-1] + + (my_rates->k13[logTlininterp_buf.indixe[i]+1-1] -my_rates->k13[logTlininterp_buf.indixe[i]-1])*eqtdef; + double heq2 = (-1. / (4.*eqk222)) * (eqk132- + std::sqrt(8.*eqk132*eqk222* + my_chemistry->HydrogenFractionByMass*local_rho+ + std::pow(eqk132,2.))); + + double eqt1 = std::fmax(std::log(tgas[i]) - 0.1*dlogtem, logTlininterp_buf.t1[i]); + eqtdef = (eqt1 - logTlininterp_buf.t1[i])/(logTlininterp_buf.t2[i] - logTlininterp_buf.t1[i]); + double eqk221 = my_rates->k22[logTlininterp_buf.indixe[i]-1] + + (my_rates->k22[logTlininterp_buf.indixe[i]+1-1] -my_rates->k22[logTlininterp_buf.indixe[i]-1])*eqtdef; + double eqk131 = my_rates->k13[logTlininterp_buf.indixe[i]-1] + + (my_rates->k13[logTlininterp_buf.indixe[i]+1-1] -my_rates->k13[logTlininterp_buf.indixe[i]-1])*eqtdef; + double heq1 = (-1. / (4.*eqk221)) * (eqk131- + std::sqrt(8.*eqk131*eqk221* + my_chemistry->HydrogenFractionByMass*local_rho+std::pow(eqk131,2.))); + + double dheq = (std::fabs(heq2-heq1)/(std::exp(eqt2) - std::exp(eqt1))) + * (tgas[i]/p2d[i]) * edot[i]; + double heq = (-1. / (4.*k22[i])) * (k13[i]- + std::sqrt(8.*k13[i]*k22[i]* + my_chemistry->HydrogenFractionByMass*local_rho+std::pow(k13[i],2.))); + + return heq / dheq; +} + +// ------------------------------------------------------------- + +/// Sets the current subcycle timestep for each index in the index-range +/// if it exceeds maximum the allowed chemistry-rate timestep. +/// +/// @param[out] dtit buffer tracking the current subcycle timestep for each +/// index in the index-range. Values will be modified in place. +/// @param[in] idx_range Specifies the current index-range +/// @param[in] iter current subcycle iteration +/// @param[in] dt tracks the full timestep that all the subcycles will +/// eventually add up to +/// @param[in] ttot tracks the total time that has already elapsed from +/// previous subcycles for each location in `idx_range` +/// @param[in] itmask Specifies the `idx_range`'s iteration-mask for the +/// Gauss-Seidel scheme +/// @param[in] itmask_nr Specifies the `idx_range`'s iteration-mask for the +/// Newton-Raphson scheme +/// @param[in] imp_eng Specifies how Newton-Raphson scheme handles energy +/// evolution at each `idx_range` location +/// @param[in] dedot, HIdot respectively specify the time derivative of the +/// free electrons and HI for the `idx_range` +/// @param[in] dedot_prev, HIdot_prev respectively specify the time derivative +/// of the free electron density and HI density for the `idx_range` from the +/// previous subcycle (they're allowed to hold garbage data in 1st subcycle) +/// @param[in] ddom specifies precomputed product of mass density and the +/// `dom` quantity for each location in `idx_range` +/// @param[in] tgas specifies the gas temperatures for the `idx_range` +/// @param[in] p2d specifies the pressures for the `idx_range`. This is +/// computed user-specified nominal adiabatic index value (i.e. no attempts +/// are made to correct for presence of H2) +/// @param[in] edot specifies the time derivative of the internal energy +/// density for the `idx_range`. +/// @param[in] my_chemistry holds a number of configuration parameters +/// @param[in] my_rates holds assorted rate data. In this function, this is +/// being used to specify the interpolation tables of some relevant reaction +/// rates (they are tabulated with respect to logT) +/// @param[in] dlogtem Specifies the constant spacing shared by the relevant +/// rate interpolation tables +/// @param[in] logTlininterp_buf Specifies the information related to the +/// position in the logT interpolations (for a number of chemistry zones) +/// @param[in] my_fields specifies the field data +/// @param[in] kcr_buf holds various pre-computed chemical reaction rates for +/// each location in `idx_range`. +/// +/// @todo +/// Consider breaking this into 2 functions that separately determines dtit for +/// Gauss-Seidel and Newton-Raphson. (At the time of writing, the included +/// logic for Newton-Raphson doesn't care about the chemistry-rates, instead +/// it sets the timestep based on the energy evolution) +static void set_subcycle_dt_from_chemistry_scheme_( + double* dtit, IndexRange idx_range, int iter, double dt, const double* ttot, + const gr_mask_type* itmask, const gr_mask_type* itmask_nr, const int* imp_eng, + double* dedot, double* HIdot, + const double* dedot_prev, const double* HIdot_prev, + const double* ddom, const double* tgas, const double* p2d, const double* edot, + const chemistry_data* my_chemistry, const chemistry_data_storage* my_rates, + double dlogtem, + const grackle::impl::LogTLinInterpScratchBuf logTlininterp_buf, + grackle_field_data* my_fields, + grackle::impl::CollisionalRxnRateCollection kcr_buf +) { + const int j = idx_range.j; + const int k = idx_range.k; + + 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 d(my_fields->density, + my_fields->grid_dimension[0], + my_fields->grid_dimension[1], + my_fields->grid_dimension[2]); + grackle::impl::View e(my_fields->internal_energy, + my_fields->grid_dimension[0], + my_fields->grid_dimension[1], + my_fields->grid_dimension[2]); + + + for (int i = idx_range.i_start; i < idx_range.i_stop; i++) { + if (itmask[i] != MASK_FALSE) { + // in this case, the chemical network will be evolved with Gauss-Seidel + + // Part 1 of 2: adjust values of dedot and HIdot + // --------------------------------------------- + + // Bound from below to prevent numerical errors + if (std::fabs(dedot[i]) < tiny8) { + dedot[i] = std::fmin(tiny_fortran_val, de(i,j,k)); + } + if (std::fabs(HIdot[i]) < tiny8){ + HIdot[i] = std::fmin(tiny_fortran_val, HI(i,j,k)); + } + + // If the net rate is almost perfectly balanced then set + // it to zero (since it is zero to available precision) + { + double ion_rate = std::fabs(kcr_buf.data[CollisionalRxnLUT::k1][i] * + de(i,j,k) * HI(i,j,k)); + double recomb_rate = std::fabs(kcr_buf.data[CollisionalRxnLUT::k2][i] * + HII(i,j,k) * de(i,j,k)); + double ratio = (std::fmin(ion_rate, recomb_rate) / + std::fmax(std::fabs(dedot[i]), std::fabs(HIdot[i]))); + if (ratio > 1.0e6) { + dedot[i] = tiny8; + HIdot[i] = tiny8; + } + } + + // If the iteration count is high then take the smaller of + // the calculated dedot and last time step's actual dedot. + // This is intended to get around the problem of a low + // electron or HI fraction which is in equilibrium with high + // individual terms (which all nearly cancel). + if (iter > 50) { + dedot[i] = std::fmin(std::fabs(dedot[i]), std::fabs(dedot_prev[i])); + HIdot[i] = std::fmin(std::fabs(HIdot[i]), std::fabs(HIdot_prev[i])); + } + + // Part 2 of 2: compute minimum rate timestep + // ------------------------------------------ + + double olddtit = dtit[i]; + dtit[i] = grackle::impl::fmin(std::fabs(0.1*de(i,j,k)/dedot[i]), + std::fabs(0.1*HI(i,j,k)/HIdot[i]), + dt-ttot[i], + 0.5*dt); + + if (ddom[i] > 1.e8 && edot[i] > 0. && + my_chemistry->primordial_chemistry > 1) { + // here, we ensure that that the equilibrium mass density of + // Hydrogen changes by 10% or less + double Heq_div_dHeqdt = calc_Heq_div_dHeqdt_( + my_chemistry, my_rates, dlogtem, logTlininterp_buf, + kcr_buf.data[CollisionalRxnLUT::k13], + kcr_buf.data[CollisionalRxnLUT::k22], + d(i,j,k), tgas, p2d, edot, i + ); + + dtit[i] = std::fmin(dtit[i], 0.1*Heq_div_dHeqdt); + } + + if (iter > 10) { + dtit[i] = std::fmin(olddtit*1.5, dtit[i]); + } + + } else if ((itmask_nr[i]!=MASK_FALSE) && (imp_eng[i]==0)) { + // we may want to handle this case and the next case in a separate + // function (they determine the timestep using very different logic than + // in the above case) + dtit[i] = grackle::impl::fmin(std::fabs(0.1*e(i,j,k)/edot[i]*d(i,j,k)), + dt-ttot[i], + 0.5*dt); + + } else if ((itmask_nr[i]!=MASK_FALSE) && (imp_eng[i]==1)) { + dtit[i] = dt - ttot[i]; + + } else { + dtit[i] = dt; + } + } +} + +// ------------------------------------------------------------- + +/// Updates the iteration mask in the case where the user has specified that we +/// are using grackle as a part of a coupled radiative transfer calculation +/// +/// @param[out] itmask the mask that will be overriden +/// @param[in] idx_range specifies the index-range +/// @param[in] my_chemistry specifies grackle settings (we probably don't need +/// to pass in everything) +/// @param[in] my_fields used to access HI photo-ionization rate field + +static inline void coupled_rt_modify_itmask_( + gr_mask_type* itmask, + IndexRange idx_range, + const chemistry_data* my_chemistry, + grackle_field_data* my_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]); + + // adjust iteration mask if the caller indicates that they're using + // Grackle in a coupled radiative-transfer/chemistry-energy calculation + // (that has intermediate steps) + if (my_chemistry->use_radiative_transfer == 1 && + my_chemistry->radiative_transfer_coupled_rate_solver == 1) { + // we only define behavior for radiative_transfer_intermediate_step + // values of 0 or 1 + + const int j = idx_range.j; + const int k = idx_range.k; + + if (my_chemistry->radiative_transfer_intermediate_step == 1) { + // the caller has invoked this chemistry-energy solver as an + // intermediate step of a coupled radiative-transfer/chemistry-energy + // calculation and they only want the solver consider cells where + // the radiation is non-zero + for (int i = idx_range.i_start; i < idx_range.i_stop; i++) { + itmask[i] = (kphHI(i,j,k) > 0) ? MASK_TRUE : MASK_FALSE; + } + } else if (my_chemistry->radiative_transfer_intermediate_step == 0) { + // the caller has invoked this chemistry-energy solver outside + // of their coupled radiative-transfer/chemistry-energy calculation. + // They want to apply the solver to cells where radiation is 0 (i.e. + // locations where skipped by the coupled calculation) + for (int i = idx_range.i_start; i < idx_range.i_stop; i++) { + itmask[i] = (kphHI(i,j,k) > 0) ? MASK_FALSE : MASK_TRUE; + } + } + } +} + +// ------------------------------------------------------------- + +namespace grackle::impl { + +/// Aggregates buffers used as scratch space in rate-related calculations +/// +/// This exists to encapsulate the logic for all of the local buffers used in +/// rate calculations within solve_rate_cool_g (this is useful given the size +/// of the function). It is not expected to be used outside of the function. +/// +/// @note +/// The purpose of this type is similar in spirit to the purpose of the +/// ScratchBuf data structrues described in the inernal_types.hpp header and +/// we consequently observe those conventions relating to (con|de)structors. +/// +/// @note +/// The majority of the time, functions only need a subset of data stored by +/// this struct. In these cases, you should strongly prefer to only pass the +/// members that are required as function arguments (if we unnecessarily pass +/// this whole struct, that can make it difficult to visualize the data flow) +struct SpeciesRateSolverScratchBuf { + + /// specifies precomputed product of mass density and the `dom` quantity for + /// each location in the index-range. Used to pick the scheme for solving the + /// rate equations and in setting the max allowed chemistry-timstep + double* ddom; + + /// buffers to hold time derivatives of free electron and HI mass densities + /// for index_range + double *dedot, *HIdot; + /// buffers to hold time derivatives of free electron and HI mass densities + /// for index_range computed during the previous cycle + double *dedot_prev, *HIdot_prev; + + /// buffer used to track the rate of H2 formation on dust grains + double* h2dust; + + /// scratch space used only within lookup_cool_rates1d_g. This is 14 times + /// larger than most of the other buffers. + /// + /// (with minimal refactoring, this buffer could probably be removed) + double *k13dd; + + /// iteration mask denoting where the Newton-Raphson scheme will be used + gr_mask_type* itmask_nr; + + /// buffer specifying how the Newton-Raphson scheme handles energy evolution + int* imp_eng; + + // buffers in the following data structure are used to temporarily hold + // the evolved density of various species as we evolve over a subcycle + // (currently only used by step_rate_g) + grackle::impl::SpeciesCollection species_tmpdens; + + // buffers in the following data structure are used to temporarily hold + // the interpolated Collisional Rxn Rates that have been + // interpolated using the standard 1D log temperature table. + grackle::impl::CollisionalRxnRateCollection kcr_buf; + + // buffers in the following data structure are used to temporarily hold + // the computed radiative reaction rates + grackle::impl::PhotoRxnRateCollection kshield_buf; + + // buffers in the following data structure are used to temporarily hold + // the interpolated chemistry-heating rates at each index-range location + grackle::impl::ChemHeatingRates chemheatrates_buf; + + // holds computed grain growth/destruction rates: + grackle::impl::GrainSpeciesCollection grain_growth_rates; + +}; + +/// allocates the contents of a new SpeciesRateSolverScratchBuf +/// +/// @param nelem The number of elements a buffer is expected to have in order +/// to store values for the standard sized index-range +SpeciesRateSolverScratchBuf new_SpeciesRateSolverScratchBuf(int nelem) { + SpeciesRateSolverScratchBuf out; + + out.ddom = (double*)malloc(sizeof(double)*nelem); + + out.dedot = (double*)malloc(sizeof(double)*nelem); + out.HIdot = (double*)malloc(sizeof(double)*nelem); + out.dedot_prev = (double*)malloc(sizeof(double)*nelem); + out.HIdot_prev = (double*)malloc(sizeof(double)*nelem); + + out.k13dd = (double*)malloc(sizeof(double)*nelem*14); + + out.h2dust = (double*)malloc(sizeof(double)*nelem); + + out.itmask_nr = (gr_mask_type*)malloc(sizeof(gr_mask_type)*nelem); + + out.imp_eng = (int*)malloc(sizeof(int)*nelem); + + out.species_tmpdens = grackle::impl::new_SpeciesCollection(nelem); + out.kcr_buf = grackle::impl::new_CollisionalRxnRateCollection(nelem); + out.kshield_buf = grackle::impl::new_PhotoRxnRateCollection(nelem); + out.chemheatrates_buf = grackle::impl::new_ChemHeatingRates(nelem); + out.grain_growth_rates = grackle::impl::new_GrainSpeciesCollection(nelem); + + return out; +} + +/// performs cleanup of the contents of SpeciesRateSolverScratchBuf +/// +/// This effectively invokes a destructor +void drop_SpeciesRateSolverScratchBuf(SpeciesRateSolverScratchBuf* ptr) { + free(ptr->ddom); + + free(ptr->dedot); + free(ptr->HIdot); + free(ptr->dedot_prev); + free(ptr->HIdot_prev); + + free(ptr->k13dd); + + free(ptr->h2dust); + + free(ptr->itmask_nr); + free(ptr->imp_eng); + + grackle::impl::drop_SpeciesCollection(&ptr->species_tmpdens); + grackle::impl::drop_CollisionalRxnRateCollection(&ptr->kcr_buf); + grackle::impl::drop_PhotoRxnRateCollection(&ptr->kshield_buf); + grackle::impl::drop_ChemHeatingRates(&ptr->chemheatrates_buf); + grackle::impl::drop_GrainSpeciesCollection(&ptr->grain_growth_rates); +} + + +} // namespace grackle::impl + +// ------------------------------------------------------------- + +#ifdef __cplusplus +extern "C" { +#endif /* __cplusplus */ + +int solve_rate_cool_g( + int imetal, double dt, InternalGrUnits internalu, + chemistry_data* my_chemistry, chemistry_data_storage* my_rates, + grackle_field_data* my_fields, photo_rate_storage* my_uvb_rates +) +{ + // shorten `grackle::impl::fortran_wrapper` to `f_wrap` within this function + namespace f_wrap = ::grackle::impl::fortran_wrapper; + +#ifdef GRACKLE_FLOAT_4 + const gr_float tolerance = (gr_float)(1.0e-05); +#else + const gr_float tolerance = (gr_float)(1.0e-10); +#endif + + // Set error indicator (we will return this value) + int ierr = GR_SUCCESS; + + // Set flag for dust-related options + const gr_mask_type anydust = + ((my_chemistry->h2_on_dust > 0) || (my_chemistry->dust_chemistry > 0)) + ? MASK_TRUE + : MASK_FALSE; + + // ignore metal chemistry/cooling below this metallicity + const double min_metallicity = 1.e-9 / my_chemistry->SolarMetalFractionByMass; + + // Set units + const double dom = internalu_calc_dom_(internalu); + const double chunit = internalu_get_chunit_(internalu); + + const double dx_cgs = my_fields->grid_dx * internalu.xbase1; + const double c_ljeans = internalu_calc_coef_ljeans_(internalu, + my_chemistry->Gamma); + + const double dlogtem = ( + (std::log(my_chemistry->TemperatureEnd) - + std::log(my_chemistry->TemperatureStart)) / + (double)(my_chemistry->NumberOfTemperatureBins-1 ) + ); + + // Convert densities from comoving to proper + + if (internalu.extfields_in_comoving == 1) { + gr_float factor = (gr_float)(std::pow(internalu.a_value,(-3)) ); + f_wrap::scale_fields_g(imetal, factor, my_chemistry, my_fields); + } + + f_wrap::ceiling_species_g(imetal, my_chemistry, my_fields); + + const grackle_index_helper idx_helper = build_index_helper_(my_fields); + + OMP_PRAGMA("omp parallel") + { + // each OMP thread separately initializes/allocates variables defined in + // the current scope and then enters the for-loop + + // holds computed grain temperatures: + grackle::impl::GrainSpeciesCollection grain_temperatures = + grackle::impl::new_GrainSpeciesCollection(my_fields->grid_dimension[0]); + + grackle::impl::LogTLinInterpScratchBuf logTlininterp_buf = + grackle::impl::new_LogTLinInterpScratchBuf(my_fields->grid_dimension[0]); + + grackle::impl::Cool1DMultiScratchBuf cool1dmulti_buf = + grackle::impl::new_Cool1DMultiScratchBuf(my_fields->grid_dimension[0]); + + grackle::impl::CoolHeatScratchBuf coolingheating_buf = + grackle::impl::new_CoolHeatScratchBuf(my_fields->grid_dimension[0]); + + // holds buffers exclusively used for solving species rate equations + // (i.e. in the future, we could have the constructor skip allocations of + // all contained data structures when using primordial_chemistry == 0) + grackle::impl::SpeciesRateSolverScratchBuf spsolvbuf = + grackle::impl::new_SpeciesRateSolverScratchBuf( + my_fields->grid_dimension[0] + ); + + // the following variables aren't embedded in structs because they are used + // in a number of different internal routines. Sorting these into + // additional structs (or leaving them free-standing) will become more + // obvious as we transcribe more routines. + std::vector dtit(my_fields->grid_dimension[0]); + std::vector ttot(my_fields->grid_dimension[0]); + std::vector p2d(my_fields->grid_dimension[0]); + std::vector tgas(my_fields->grid_dimension[0]); + std::vector tdust(my_fields->grid_dimension[0]); + std::vector metallicity(my_fields->grid_dimension[0]); + std::vector dust2gas(my_fields->grid_dimension[0]); + std::vector rhoH(my_fields->grid_dimension[0]); + std::vector mmw(my_fields->grid_dimension[0]); + std::vector edot(my_fields->grid_dimension[0]); + + // iteration masks + std::vector itmask(my_fields->grid_dimension[0]); + std::vector itmask_tmp(my_fields->grid_dimension[0]); + std::vector itmask_metal(my_fields->grid_dimension[0]); + + // create views of density and internal energy fields to support 3D access + grackle::impl::View d(my_fields->density, + my_fields->grid_dimension[0], + my_fields->grid_dimension[1], + my_fields->grid_dimension[2]); + grackle::impl::View e(my_fields->internal_energy, + my_fields->grid_dimension[0], + my_fields->grid_dimension[1], + my_fields->grid_dimension[2]); + + // The following for-loop is a flattened loop over every k,j combination. + // OpenMP divides this loop between all threads. Within the loop, we + // complete calculations for the constructed index-range construct + // (an index range corresponds to an "i-slice") + OMP_PRAGMA("omp for schedule(runtime)") + for (int t = 0; t < idx_helper.outer_ind_size; t++) { + // construct an index-range corresponding to "i-slice" + const IndexRange idx_range = make_idx_range_(t, &idx_helper); + const int k = idx_range.k; // use 0-based index + const int j = idx_range.j; // use 0-based index + + // `tolerance = 1.0e-06_DKIND * dt` was some commented logic in the + // original fortran subroutine in this location + + // Initialize iteration mask to true for all cells. + for (int i = idx_range.i_start; i < idx_range.i_stop; i++) { + itmask[i] = MASK_TRUE; + } + + // adjust iteration mask (but only if using Grackle in a coupled + // radiative-transfer calculation) + coupled_rt_modify_itmask_(itmask.data(), idx_range, my_chemistry, + my_fields); + + // Set time elapsed to zero for each cell in 1D section + + for (int i = idx_range.i_start; i < idx_range.i_stop; i++) { + ttot[i] = 0.; + } + + // A useful slice variable since we do this a lot + // -> we don't need it for primordial_chemistry==0 + for (int i = idx_range.i_start; i < idx_range.i_stop; i++) { + spsolvbuf.ddom[i] = d(i,j,k) * dom; + } + + // declare 2 variables (primarily used for subcycling, but also used in + // error reporting) + int iter; + double ttmin; + + // ------------------ Loop over subcycles ---------------- + + for (iter = 1; iter<=(my_chemistry->max_iterations); iter++) { + + for (int i = idx_range.i_start; i < idx_range.i_stop; i++) { + if (itmask[i] != MASK_FALSE) { + dtit[i] = huge8; + } + } + + // Compute the cooling rate, tgas, tdust, and metallicity for this row + f_wrap::cool1d_multi_g( + imetal, idx_range, iter, edot.data(), tgas.data(), + mmw.data(), p2d.data(), tdust.data(), metallicity.data(), + dust2gas.data(), rhoH.data(), itmask.data(), itmask_metal.data(), + my_chemistry, my_rates, my_fields, *my_uvb_rates, internalu, + grain_temperatures, logTlininterp_buf, cool1dmulti_buf, + coolingheating_buf + ); + + if (my_chemistry->primordial_chemistry == 0) { + // This is some basic book-keeping to ensure that itmask_tmp has + // sensible values when ispecies is 0 + // -> see the comment following this if-else statement suggesting how + // we could refactor itmask-handling (eliminating this branch) + std::memcpy(itmask_tmp.data(), itmask.data(), + sizeof(gr_mask_type)*my_fields->grid_dimension[0]); + + } else { + + // Look-up rates as a function of temperature for 1D set of zones + // (maybe should add itmask to this call) + + f_wrap::lookup_cool_rates1d_g( + idx_range, anydust, tgas.data(), mmw.data(), tdust.data(), + dust2gas.data(), spsolvbuf.k13dd, spsolvbuf.h2dust, + dom, dx_cgs, c_ljeans, itmask.data(), itmask_metal.data(), + imetal, rhoH.data(), dt, my_chemistry, my_rates, my_fields, + *my_uvb_rates, internalu, spsolvbuf.grain_growth_rates, + grain_temperatures, logTlininterp_buf, spsolvbuf.kcr_buf, + spsolvbuf.kshield_buf, spsolvbuf.chemheatrates_buf + ); + + // Compute dedot and HIdot, the rates of change of de and HI + // (should add itmask to this call) + + f_wrap::rate_timestep_g( + spsolvbuf.dedot, spsolvbuf.HIdot, anydust, idx_range, + spsolvbuf.h2dust, rhoH.data(), itmask.data(), edot.data(), + chunit, dom, my_chemistry, my_fields, *my_uvb_rates, + spsolvbuf.kcr_buf, spsolvbuf.kshield_buf, + spsolvbuf.chemheatrates_buf + ); + + // First, copy itmask's current values into a temporary array + // (itmask_tmp). Then setup masks to identify which chemistry schemes + // to use. We split cells by density: + // => low-density: Gauss-Seidel scheme, tracked by itmask + // => high-density: Newton-Raphson scheme, tracked by itmask_nr + // + // (the values stored within itmask will change within the function) + select_chem_scheme_update_masks_( + idx_range, itmask.data(), spsolvbuf.itmask_nr, spsolvbuf.imp_eng, + itmask_tmp.data(), my_fields->grid_dimension[0], imetal, + min_metallicity, spsolvbuf.ddom, tgas.data(), + metallicity.data(), my_chemistry + ); + + // Set the max timestep for the current subcycle based on our scheme + // for updating the chemical network: + // - for Gauss-Seidel, pick a timestep that keeps relative chemical + // changes below 10% + // - do something else for Newton-Raphson + + set_subcycle_dt_from_chemistry_scheme_( + dtit.data(), idx_range, iter, dt, ttot.data(), itmask.data(), + spsolvbuf.itmask_nr, spsolvbuf.imp_eng, + spsolvbuf.dedot, spsolvbuf.HIdot, + spsolvbuf.dedot_prev, spsolvbuf.HIdot_prev, + spsolvbuf.ddom, tgas.data(), p2d.data(), edot.data(), + my_chemistry, my_rates, dlogtem, logTlininterp_buf, my_fields, + spsolvbuf.kcr_buf + ); + } + + // TODO: Consider refactoring the iteration mask handling: + // 1. introduce `itmask_gs` as a member of `spsolvbuf` and have the + // `select_chem_scheme_update_masks_` function store locations in + // `spsolvbuf.itmask_gs` where we will apply Gauss-Seidel scheme + // 2. replace `itmask` with `spsolvbuf.itmask_gs` in arg-lists of + // `set_subcycle_dt_from_chemistry_scheme_` & `f_wrap::step_rate_g` + // 3. insert following chunk of logic RIGHT HERE + // > const gr_mask_type* energy_itmask = + // > (my_chemistry->primordial_chemistry == 0) + // > ? itmask.data() : spsolvbuf.itmask_gs; + // 4. replace `itmask` with `energy_itmask` in arg-list of + // `enforce_max_heatcool_subcycle_dt_` & within energy update loop + // 5. stop modifying `itmask` in `select_chem_scheme_update_masks_` + // 6. remove declaration of `itmask_tmp`, logic that initializes, and + // the loop that uses it to override `itmask` + + // Update dtit (the current subcycle timestep) to ensure it doesn't + // exceed the max timestep for cooling/heating + // -> zones that will use Newton-Raphson scheme are ignored + enforce_max_heatcool_subcycle_dt_( + dtit.data(), idx_range, dt, ttot.data(), itmask.data(), tgas.data(), + p2d.data(), edot.data(), my_chemistry + ); + + // Update total and gas energy + // -> zones that will use Newton-Raphson scheme are ignored + if (my_chemistry->with_radiative_cooling == 1) { + for (int i = idx_range.i_start; i < idx_range.i_stop; i++) { + if (itmask[i] != MASK_FALSE) { + e(i,j,k) = e(i,j,k) + (gr_float)(edot[i]/d(i,j,k)*dtit[i]); + } + } + } + + if (my_chemistry->primordial_chemistry > 0) { + + // Solve rate equations with one linearly implicit Gauss-Seidel + // sweep of a backward Euler method (for all cells specified by + // itmask) + + f_wrap::step_rate_g( + dtit.data(), idx_range, anydust, spsolvbuf.h2dust, rhoH.data(), + spsolvbuf.dedot_prev, spsolvbuf.HIdot_prev, itmask.data(), + itmask_metal.data(), imetal, 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 + // sweep of a backward Euler method (for all cells specified by + // itmask_nr) + + f_wrap::step_rate_newton_raphson( + imetal, idx_range, iter, dom, chunit, dx_cgs, c_ljeans, + dtit.data(), p2d.data(), tgas.data(), tdust.data(), + metallicity.data(), dust2gas.data(), rhoH.data(), mmw.data(), + spsolvbuf.h2dust, edot.data(), anydust, spsolvbuf.itmask_nr, + itmask_metal.data(), spsolvbuf.imp_eng, my_chemistry, my_rates, + my_fields, *my_uvb_rates, internalu, grain_temperatures, + logTlininterp_buf, cool1dmulti_buf, coolingheating_buf, + spsolvbuf.chemheatrates_buf + ); + + } + + // restore the values of itmask from the backup that we made earlier + for (int i = idx_range.i_start; i < idx_range.i_stop; i++) { + itmask[i] = itmask_tmp[i]; + } + + // Add the timestep to the elapsed time for each cell and find + // minimum elapsed time step in this row + ttmin = huge8; + for (int i = idx_range.i_start; i < idx_range.i_stop; i++) { + ttot[i] = std::fmin(ttot[i] + dtit[i], dt); + + if (std::fabs(dt-ttot[i]) < tolerance*dt) { itmask[i] = MASK_FALSE; } + + if (ttot[i] my_chemistry->max_iterations) { + OMP_PRAGMA_CRITICAL + { + printf("inside if statement solve rate cool: %d %d\n", + my_fields->grid_start[0], + my_fields->grid_end[0]); + eprintf("MULTI_COOL iter > %d at j_0based,k_0based = %d %d\n", + my_chemistry->max_iterations, idx_range.j, idx_range.k); + printf("FATAL error (2) in MULTI_COOL\n"); + printf("( dt = %.17e ttmin = %.17e )", dt, ttmin); + grackle::impl::print_contiguous_row_( + dtit.data(), my_fields->grid_start[0], my_fields->grid_end[0]+1 + ); + grackle::impl::print_contiguous_row_( + ttot.data(), my_fields->grid_start[0], my_fields->grid_end[0]+1 + ); + grackle::impl::print_contiguous_row_( + edot.data(), my_fields->grid_start[0], my_fields->grid_end[0]+1 + ); + grackle::impl::print_contiguous_row_( + itmask.data(), my_fields->grid_start[0], my_fields->grid_end[0]+1 + ); + + if (my_chemistry->exit_after_iterations_exceeded == 1) { + ierr = GR_FAIL; + } + } // OMP_PRAGMA_CRITICAL + } + + if (iter > my_chemistry->max_iterations/2) { // WARNING_MESSAGE + OMP_PRAGMA_CRITICAL + { + eprintf("MULTI_COOL iter,j_0based,k_0based = %d %d %d\n", + iter, idx_range.j, idx_range.k); + } + } + + } // outer-loop (index t) - each of these correspond to j,k pairs + + // cleanup manually allocated temporaries + grackle::impl::drop_GrainSpeciesCollection(&grain_temperatures); + grackle::impl::drop_LogTLinInterpScratchBuf(&logTlininterp_buf); + grackle::impl::drop_Cool1DMultiScratchBuf(&cool1dmulti_buf); + grackle::impl::drop_CoolHeatScratchBuf(&coolingheating_buf); + + grackle::impl::drop_SpeciesRateSolverScratchBuf(&spsolvbuf); + + } // OMP_PRAGMA("omp parallel") + + // If an error has been produced, return now. + + if (ierr != GR_SUCCESS) { + return ierr; + } + + // Convert densities back to comoving from proper + + if (internalu.extfields_in_comoving == 1) { + gr_float factor = (gr_float)(std::pow(internalu.a_value,3) ); + f_wrap::scale_fields_g(imetal, factor, my_chemistry, my_fields); + } + + if (my_chemistry->primordial_chemistry > 0) { + + // Correct the species to ensure consistency (i.e. type conservation) + + f_wrap::make_consistent_g(imetal, dom, my_chemistry, my_rates, my_fields); + + } + + return ierr; +} + +#ifdef __cplusplus +} // extern "C" +#endif /* __cplusplus */ diff --git a/src/clib/solve_rate_cool_g-cpp.h b/src/clib/solve_rate_cool_g-cpp.h new file mode 100644 index 000000000..76a2868ea --- /dev/null +++ b/src/clib/solve_rate_cool_g-cpp.h @@ -0,0 +1,50 @@ +// See LICENSE file for license and copyright information + +/// @file solve_rate_cool_g-cpp.h +/// @brief Declares signature of solve_rate_cool_g + +// This file was initially generated automatically during conversion of the +// solve_rate_cool_g function from FORTRAN to C++ + +#ifndef MY_FILE_CPP_H +#define MY_FILE_CPP_H + +#include "grackle.h" // gr_float +#include "fortran_func_decls.h" // gr_mask_int +#include "internal_units.h" // InternalGrUnits + +#ifdef __cplusplus +extern "C" { +#endif /* __cplusplus */ +// the following function can be called from C or C++ + +/// Solve the multi-species rate and cooling equations +/// +/// @par History +/// written by: Yu Zhang, Peter Anninos and Tom Abel +/// modified1: January, 1996 by Greg Bryan; converted to KRONOS +/// modified2: October, 1996 by GB; adapted to AMR +/// modified3: May, 1999 by GB and Tom Abel, 3bodyH2, solver, HD +/// modified4: June, 2005 by GB to solve rate & cool at same time +/// modified5: April, 2009 by JHW to include radiative transfer +/// modified6: September, 2009 by BDS to include cloudy cooling +/// modified7: January, 2025 by Matthew Abruzzo; ported to C++ +/// +/// @return Returns GR_SUCCESS or GR_FAIL to indicate whether there was an error +/// +/// @todo +/// Once the file where this routine called is adjusted to be compiled with a +/// C++ compiler, modify this function (prototype & implementation) such that: +/// - it's not enclosed by a `extern "C"` block +/// - it's defined within a `grackle::impl` namespace +int solve_rate_cool_g( + int imetal, double dt, InternalGrUnits internalu, + chemistry_data* my_chemistry, chemistry_data_storage* my_rates, + grackle_field_data* my_fields, photo_rate_storage* my_uvb_rates +); + +#ifdef __cplusplus +} // extern "C" +#endif /* __cplusplus */ + +#endif /* MY_FILE_CPP_H */ diff --git a/src/clib/solve_rate_cool_g.F b/src/clib/solve_rate_cool_g.F index 18d30e54b..cfbc200be 100644 --- a/src/clib/solve_rate_cool_g.F +++ b/src/clib/solve_rate_cool_g.F @@ -1,1460 +1,9 @@ #include "phys_const.def" -!======================================================================= -!/////////////////// SUBROUTINE SOLVE_RATE_COOL_G \\\\\\\\\\\\\\\\\\\\\ - - subroutine solve_rate_cool_g(icool, d, e, u, v, w, de, - & HI, HII, HeI, HeII, HeIII, - & in, jn, kn, nratec, iexpand, - & ispecies, imetal, imcool, idust, - & idustall, idustfield, idim, - & is, js, ks, ie, je, ke, ih2co, ipiht, idustrec, - & igammah, - & dx, dt, aye, temstart, temend, - & utem, uxyz, uaye, urho, utim, - & gamma, fh, dtoh, z_solar, fgr, - & k1a, k2a, k3a, k4a, k5a, k6a, k7a, k8a, k9a, k10a, - & k11a, k12a, k13a, k13dda, k14a, k15a, - & k16a, k17a, k18a, k19a, k22a, - & k24, k25, k26, k27, k28, k29, k30, k31, - & k50a, k51a, k52a, k53a, k54a, k55a, k56a, - & k57a, k58a, - & ndratec, dtemstart, dtemend, h2dusta, - & ncrna, ncrd1a, ncrd2a, - & ceHIa, ceHeIa, ceHeIIa, ciHIa, ciHeIa, - & ciHeISa, ciHeIIa, reHIIa, reHeII1a, - & reHeII2a, reHeIIIa, brema, compa, gammaha, isrf, - & regra, gamma_isrfa, - & comp_xraya, comp_temp, piHI, piHeI, piHeII, - & HM, H2I, H2II, DI, DII, HDI, metal, dust, - & hyd01ka, h2k01a, vibha, rotha, rotla, - & gpldla, gphdla, hdltea, hdlowa, - & gaHIa, gaH2a, gaHea, gaHpa, gaela, - & h2ltea, gasgra, iH2shield, iradshield, - & avgsighi, avgsighei, avgsigheii, - & iradtrans, iradcoupled, iradstep, - & irt_honly, kphHI, kphHeI, kphHeII, kdissH2I, - & photogamma, xH2shield, - & ierr, - & ih2optical, iciecool, ithreebody, ih2cr, ihdcr, - & ciecoa, - & icmbTfloor, iClHeat, clEleFra, - & priGridRank, priGridDim, - & priPar1, priPar2, priPar3, priPar4, priPar5, - & priDataSize, priCooling, priHeating, priMMW, - & metGridRank, metGridDim, - & metPar1, metPar2, metPar3, metPar4, metPar5, - & metDataSize, metCooling, metHeating, clnew, - & iVheat, iMheat, Vheat, Mheat, - & iTfloor, Tfloor_scalar, Tfloor, - & imchem, igrgr, ipcont, tmcool - & , DM, HDII, HeHII - & , 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 - & , k125a, k129a, k130a, k131a, k132a - & , k133a, k134a, k135a, k136a, k137a - & , k148a, k149a, k150a, k151a, k152a - & , k153a - & , kz15a, kz16a, kz17a, kz18a, kz19a - & , kz20a, kz21a, kz22a, kz23a, kz24a - & , kz25a, kz26a, kz27a, kz28a, kz29a - & , kz30a, kz31a, kz32a, kz33a, kz34a - & , kz35a, kz36a, kz37a, kz38a, kz39a - & , kz40a, kz41a, kz42a, kz43a, kz44a - & , kz45a, kz46a, kz47a, kz48a, kz49a - & , kz50a, kz51a, kz52a, kz53a, kz54a - & , cieY06a - & , LH2_N, LH2_Size - & , LH2_D, LH2_T, LH2_H - & , LH2_dD, LH2_dT, LH2_dH, LH2_L - & , LHD_N, LHD_Size - & , LHD_D, LHD_T, LHD_H - & , LHD_dD, LHD_dT, LHD_dH, LHD_L - & , LCI_N, LCI_Size - & , LCI_D, LCI_T, LCI_H - & , LCI_dD, LCI_dT, LCI_dH, LCI_L - & , LCII_N, LCII_Size - & , LCII_D, LCII_T, LCII_H - & , LCII_dD, LCII_dT, LCII_dH, LCII_L - & , LOI_N, LOI_Size - & , LOI_D, LOI_T, LOI_H - & , LOI_dD, LOI_dT, LOI_dH, LOI_L - & , LCO_N, LCO_Size - & , LCO_D, LCO_T, LCO_H - & , LCO_dD, LCO_dT, LCO_dH, LCO_L - & , LOH_N, LOH_Size - & , LOH_D, LOH_T, LOH_H - & , LOH_dD, LOH_dT, LOH_dH, LOH_L - & , LH2O_N, LH2O_Size - & , LH2O_D, LH2O_T, LH2O_H - & , LH2O_dD, LH2O_dT, LH2O_dH, LH2O_L - & , alphap_N, alphap_Size - & , alphap_D, alphap_T, alphap_dD, alphap_dT - & , alphap_Data - & , immulti, imabund, idspecies, itdmulti, idsub - & , metal_loc - & , metal_C13, metal_C20, metal_C25, metal_C30 - & , metal_F13, metal_F15, metal_F50, metal_F80 - & , metal_P170, metal_P200, metal_Y19 - & , SN0_N - & , SN0_XC , SN0_XO, SN0_XMg, SN0_XAl - & , SN0_XSi, SN0_XS, SN0_XFe - & , SN0_fC , SN0_fO, SN0_fMg, SN0_fAl - & , SN0_fSi, SN0_fS, SN0_fFe - & , SN0_fSiM, SN0_fFeM, SN0_fMg2SiO4, SN0_fMgSiO3 - & , SN0_fFe3O4, SN0_fAC, SN0_fSiO2D, SN0_fMgO - & , SN0_fFeS, SN0_fAl2O3 - & , SN0_freforg, SN0_fvolorg, SN0_fH2Oice - & , SN0_r0SiM, SN0_r0FeM, SN0_r0Mg2SiO4, SN0_r0MgSiO3 - & , SN0_r0Fe3O4, SN0_r0AC, SN0_r0SiO2D, SN0_r0MgO - & , SN0_r0FeS, SN0_r0Al2O3 - & , SN0_r0reforg, SN0_r0volorg, SN0_r0H2Oice - & , gr_N, gr_Size, gr_dT, gr_Td - & , SN0_kpSiM, SN0_kpFeM, SN0_kpMg2SiO4, SN0_kpMgSiO3 - & , SN0_kpFe3O4, SN0_kpAC, SN0_kpSiO2D, SN0_kpMgO - & , SN0_kpFeS, SN0_kpAl2O3 - & , SN0_kpreforg, SN0_kpvolorg, SN0_kpH2Oice - & , h2dustSa, h2dustCa, gasgr2a, gamma_isrf2a, grogra - & , idissHDI, kdissHDI, iionZ, kphCI, kphOI - & , idissZ, kdissCO, kdissOH, kdissH2O, iuseH2shield, - & iisrffield, isrf_habing, - & iH2shieldcustom, f_shield_custom, - & itmax, exititmax) - -! -! SOLVE MULTI-SPECIES RATE EQUATIONS AND RADIATIVE COOLING -! -! written by: Yu Zhang, Peter Anninos and Tom Abel -! date: -! modified1: January, 1996 by Greg Bryan; converted to KRONOS -! modified2: October, 1996 by GB; adapted to AMR -! modified3: May, 1999 by GB and Tom Abel, 3bodyH2, solver, HD -! modified4: June, 2005 by GB to solve rate & cool at same time -! modified5: April, 2009 by JHW to include radiative transfer -! modified6: September, 2009 by BDS to include cloudy cooling -! -! PURPOSE: -! Solve the multi-species rate and cool equations. -! -! INPUTS: -! icool - flag to update energy from radiative cooling -! in,jn,kn - dimensions of 3D fields -! -! d - total density field -! de - electron density field -! HI,HII - H density fields (neutral & ionized) -! HeI/II/III - He density fields -! DI/II - D density fields (neutral & ionized) -! HDI - neutral HD molecule density field -! HM - H- density field -! H2I - H_2 (molecular H) density field -! H2II - H_2+ density field -! metal - metal density field -! dust - dust density field -! kph* - photoionization fields -! gamma* - photoheating fields -! f_shield_custom - custom H2 shielding factor -! -! is,ie - start and end indices of active region (zero based) -! iexpand - comoving coordinates flag (0 = off, 1 = on) -! idim - dimensionality (rank) of problem -! ispecies - chemistry module (1 - H/He only, 2 - molecular H, 3 - D) -! imetal - flag if metal field is active (0 = no, 1 = yes) -! imcool - flag if there is metal cooling -! idust - flag for H2 formation on dust grains -! idustall - flag for dust (0 - none, 1 - heating/cooling + H2 form.) -! idustfield - flag if a dust density field is present -! iisrffield - flag if a field for the interstellar radiation field is present -! ih2co - flag to include H2 cooling (1 = on, 0 = off) -! ipiht - flag to include photoionization heating (1 = on, 0 = off) -! idustrec - flag to include dust recombination cooling (1 = on, -1 = off) -! iH2shield - flag for approximate self-shielding of H2 (Wolcott-Green+ 2011) -! iradshield - flag for approximate self-shielding of UV background -! avgsighi - spectrum averaged ionization crs for HI for use with shielding -! avgsighei - spectrum averaged ionization crs for HeI for use with shielding -! avgsigheii - spectrum averaged ionization crs for HeII for use with shielding -! iradtrans - flag to include radiative transfer (1 = on, 0 = off) -! iradcoupled - flag to indicate coupled radiative transfer -! iradstep - flag to indicate intermediate coupled radiative transfer timestep -! irt_honly - flag to indicate applying RT ionization and heating to HI only -! iH2shieldcustom - flag to indicate a custom H2 shielding factor is provided - -! fh - Hydrogen mass fraction (typically 0.76) -! dtoh - Deuterium to H mass ratio -! z_solar - Solar metal mass fraction -! fgr - the local dust to gas ratio (by mass) -! dt - timestep to integrate over -! aye - expansion factor (in code units) -! -! utim - time units (i.e. code units to CGS conversion factor) -! uaye - expansion factor conversion factor (uaye = 1/(1+zinit)) -! urho - density units -! uxyz - length units -! utem - temperature(-like) units -! -! temstart, temend - start and end of temperature range for rate table -! nratec - dimensions of chemical rate arrays (functions of temperature) -! dtemstart, dtemend - start and end of dust temperature range -! ndratec - extra dimension for H2 formation on dust rate (dust temperature) -! -! icmbTfloor - flag to include temperature floor from cmb -! iClHeat - flag to include cloudy heating -! priGridRank - rank of cloudy primordial cooling data grid -! priGridDim - array containing dimensions of cloudy primordial data -! priPar1, priPar2, priPar3 - arrays containing primordial grid parameter values -! priDataSize - total size of flattened 1D primordial cooling data array -! priCooling - primordial cooling data -! priHeating - primordial heating data -! priMMW - primordial mmw data -! metGridRank - rank of cloudy metal cooling data grid -! metGridDim - array containing dimensions of cloudy metal data -! metPar1, metPar2, metPar3 - arrays containing metal grid parameter values -! metDataSize - total size of flattened 1D metal cooling data array -! metCooling - metal cooling data -! metHeating - metal heating data -! iVheat - flag for using volumetric heating rate -! iMheat - flag for using specific heating rate -! Vheat - array of volumetric heating rates -! Mheat - array of specific heating rates -! iTfloor - flag for using temperature floor field -! Tfloor_scalar - scalar temperature floor value -! Tfloor - array of temperature floor values -! itmax - maximum allowed sub-cycle iterations -! exititmax - flag to exit if max iterations exceeded -! -! OUTPUTS: -! update chemical rate densities (HI, HII, etc) -! -! PARAMETERS: -! mh - H mass in cgs units -! -!----------------------------------------------------------------------- - - implicit NONE -#include "grackle_fortran_types.def" -#ifdef _OPENMP -#include "omp_lib.h" -#endif - -! General Arguments - - integer icool, in, jn, kn, is, js, ks, ie, je, ke, nratec, - & iexpand, ih2co, ipiht, ispecies, imetal, idim, - & ierror, imcool, idust, idustall, idustfield, idustrec, - & igammah, ih2optical, iciecool, ithreebody, ih2cr, ihdcr, - & ndratec, clnew, iVheat, iMheat, iTfloor, - & iH2shield, iradshield, - & iradtrans, iradcoupled, iradstep, irt_honly, - & imchem, igrgr, ipcont, - & iisrffield, iH2shieldcustom, ierr, itmax, exititmax - - real*8 dx, dt, aye, temstart, temend, gamma, - & utim, uxyz, uaye, urho, utem, fh, dtoh, z_solar, - & fgr, dtemstart, dtemend, clEleFra, Tfloor_scalar, - & tmcool - -! Density, energy and velocity fields 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), - & HM(in,jn,kn), H2I(in,jn,kn), H2II(in,jn,kn), - & DI(in,jn,kn), DII(in,jn,kn), HDI(in,jn,kn), - & d(in,jn,kn), e(in,jn,kn), - & u(in,jn,kn), v(in,jn,kn), w(in,jn,kn), - & metal(in,jn,kn), dust(in,jn,kn), - & Vheat(in,jn,kn), Mheat(in,jn,kn), Tfloor(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) - R_PREC metal_loc(in,jn,kn) - & , metal_C13(in,jn,kn), metal_C20(in,jn,kn) - & , metal_C25(in,jn,kn), metal_C30(in,jn,kn) - & , metal_F13(in,jn,kn), metal_F15(in,jn,kn) - & , metal_F50(in,jn,kn), metal_F80(in,jn,kn) - & , metal_P170(in,jn,kn), metal_P200(in,jn,kn) - & , metal_Y19(in,jn,kn) - -! Radiative transfer fields - - R_PREC kphHI(in,jn,kn), kphHeI(in,jn,kn), kphHeII(in,jn,kn), - & kdissH2I(in,jn,kn), photogamma(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) - integer iuseH2shield - -! H2 self-shielding length-scale field - - R_PREC xH2shield(in,jn,kn) - -! Interstellar radiation field for dust heating - - R_PREC isrf_habing(in,jn,kn) - -! Custom H2 shielding factor - - R_PREC f_shield_custom(in, jn, kn) - -! Cooling tables (coolings rates as a function of temperature) - - real*8 hyd01ka(nratec), h2k01a(nratec), vibha(nratec), - & rotha(nratec), rotla(nratec), gpldla(nratec), - & gphdla(nratec), hdltea(nratec), hdlowa(nratec), - & gaHIa(nratec), gaH2a(nratec), gaHea(nratec), - & gaHpa(nratec), gaela(nratec), h2ltea(nratec), - & gasgra(nratec), ciecoa(nratec), - & ceHIa(nratec), ceHeIa(nratec), ceHeIIa(nratec), - & ciHIa(nratec), ciHeIa(nratec), ciHeISa(nratec), - & ciHeIIa(nratec), reHIIa(nratec), reHeII1a(nratec), - & reHeII2a(nratec), reHeIIIa(nratec), brema(nratec), - & compa, piHI, piHeI, piHeII, comp_xraya, comp_temp, - & gammaha, isrf, regra(nratec), gamma_isrfa - real*8 cieY06a(nratec) - integer*8 LH2_N(3), LH2_Size - real*8 LH2_D(LH2_N(1)), LH2_T(LH2_N(2)), LH2_H(LH2_N(3)) - & , LH2_dD, LH2_dT, LH2_dH, LH2_L(LH2_Size) - integer*8 LHD_N(3), LHD_Size - real*8 LHD_D(LHD_N(1)), LHD_T(LHD_N(2)), LHD_H(LHD_N(3)) - & , LHD_dD, LHD_dT, LHD_dH, LHD_L(LHD_Size) - integer*8 LCI_N(3), LCI_Size - real*8 LCI_D(LCI_N(1)), LCI_T(LCI_N(2)), LCI_H(LCI_N(3)) - & , LCI_dD, LCI_dT, LCI_dH, LCI_L(LCI_Size) - integer*8 LCII_N(3), LCII_Size - real*8 LCII_D(LCII_N(1)), LCII_T(LCII_N(2)), LCII_H(LCII_N(3)) - & , LCII_dD, LCII_dT, LCII_dH, LCII_L(LCII_Size) - integer*8 LOI_N(3), LOI_Size - real*8 LOI_D(LOI_N(1)), LOI_T(LOI_N(2)), LOI_H(LOI_N(3)) - & , LOI_dD, LOI_dT, LOI_dH, LOI_L(LOI_Size) - integer*8 LCO_N(3), LCO_Size - real*8 LCO_D(LCO_N(1)), LCO_T(LCO_N(2)), LCO_H(LCO_N(3)) - & , LCO_dD, LCO_dT, LCO_dH, LCO_L(LCO_Size) - integer*8 LOH_N(3), LOH_Size - real*8 LOH_D(LOH_N(1)), LOH_T(LOH_N(2)), LOH_H(LOH_N(3)) - & , LOH_dD, LOH_dT, LOH_dH, LOH_L(LOH_Size) - integer*8 LH2O_N(3), LH2O_Size - real*8 LH2O_D(LH2O_N(1)), LH2O_T(LH2O_N(2)), LH2O_H(LH2O_N(3)) - & , LH2O_dD, LH2O_dT, LH2O_dH, LH2O_L(LH2O_Size) - integer*8 alphap_N(2), alphap_Size - real*8 alphap_D(alphap_N(1)), alphap_T(alphap_N(2)) - & , alphap_dD, alphap_dT - & , alphap_Data(alphap_Size) - integer immulti, imabund, idspecies, itdmulti, idsub - integer SN0_N - real*8 SN0_XC (SN0_N), SN0_XO(SN0_N), SN0_XMg(SN0_N) - & , SN0_XAl(SN0_N), SN0_XSi(SN0_N), SN0_XS(SN0_N) - & , SN0_XFe(SN0_N) - real*8 SN0_fC (SN0_N), SN0_fO(SN0_N), SN0_fMg(SN0_N) - & , SN0_fAl(SN0_N), SN0_fSi(SN0_N), SN0_fS(SN0_N) - & , SN0_fFe(SN0_N) - real*8 SN0_fSiM(SN0_N), SN0_fFeM(SN0_N), SN0_fMg2SiO4(SN0_N) - & , SN0_fMgSiO3(SN0_N), SN0_fFe3O4(SN0_N), SN0_fAC(SN0_N) - & , SN0_fSiO2D(SN0_N), SN0_fMgO(SN0_N), SN0_fFeS(SN0_N) - & , SN0_fAl2O3(SN0_N) - & , SN0_freforg(SN0_N), SN0_fvolorg(SN0_N), SN0_fH2Oice(SN0_N) - real*8 SN0_r0SiM(3,SN0_N), SN0_r0FeM(3,SN0_N) - & , SN0_r0Mg2SiO4(3,SN0_N), SN0_r0MgSiO3(3,SN0_N) - & , SN0_r0Fe3O4(3,SN0_N), SN0_r0AC(3,SN0_N) - & , SN0_r0SiO2D(3,SN0_N), SN0_r0MgO(3,SN0_N) - & , SN0_r0FeS(3,SN0_N), SN0_r0Al2O3(3,SN0_N) - & , SN0_r0reforg(3,SN0_N) - & , SN0_r0volorg(3,SN0_N), SN0_r0H2Oice(3,SN0_N) -! opacity table - integer gr_N(2), gr_Size - real*8 gr_dT, gr_Td(gr_N(2)) - real*8 SN0_kpSiM(gr_Size,SN0_N), SN0_kpFeM(gr_Size,SN0_N) - & , SN0_kpMg2SiO4(gr_Size,SN0_N), SN0_kpMgSiO3(gr_Size,SN0_N) - & , SN0_kpFe3O4(gr_Size,SN0_N), SN0_kpAC(gr_Size,SN0_N) - & , SN0_kpSiO2D(gr_Size,SN0_N), SN0_kpMgO(gr_Size,SN0_N) - & , SN0_kpFeS(gr_Size,SN0_N), SN0_kpAl2O3(gr_Size,SN0_N) - & , SN0_kpreforg(gr_Size,SN0_N) - & , SN0_kpvolorg(gr_Size,SN0_N), SN0_kpH2Oice(gr_Size,SN0_N) - real*8 gasgr2a(nratec), gamma_isrf2a - - real*8 avgsighi, avgsighei, avgsigheii - -! Chemistry tables (rates as a function of temperature) - - real*8 k1a (nratec), k2a (nratec), k3a (nratec), k4a (nratec), - & k5a (nratec), k6a (nratec), k7a (nratec), k8a (nratec), - & k9a (nratec), k10a(nratec), k11a(nratec), k12a(nratec), - & k13a(nratec), k14a(nratec), k15a(nratec), k16a(nratec), - & k17a(nratec), k18a(nratec), k19a(nratec), k22a(nratec), - & k50a(nratec), k51a(nratec), k52a(nratec), k53a(nratec), - & k54a(nratec), k55a(nratec), k56a(nratec), - & k57a(nratec), k58a(nratec), - & k13dda(nratec, 14), h2dusta(nratec, ndratec), - & ncrna(nratec), ncrd1a(nratec), ncrd2a(nratec), - & k24, k25, k26, k27, k28, k29, k30, k31 - real*8 k125a(nratec), k129a(nratec), k130a(nratec) - & , k131a(nratec), k132a(nratec), k133a(nratec) - & , k134a(nratec), k135a(nratec), k136a(nratec) - & , k137a(nratec), k148a(nratec), k149a(nratec) - & , k150a(nratec), k151a(nratec), k152a(nratec) - & , k153a(nratec) - & , kz15a(nratec), kz16a(nratec), kz17a(nratec) - & , kz18a(nratec), kz19a(nratec), kz20a(nratec) - & , kz21a(nratec), kz22a(nratec), kz23a(nratec) - & , kz24a(nratec), kz25a(nratec), kz26a(nratec) - & , kz27a(nratec), kz28a(nratec), kz29a(nratec) - & , kz30a(nratec), kz31a(nratec), kz32a(nratec) - & , kz33a(nratec), kz34a(nratec), kz35a(nratec) - & , kz36a(nratec), kz37a(nratec), kz38a(nratec) - & , kz39a(nratec), kz40a(nratec), kz41a(nratec) - & , kz42a(nratec), kz43a(nratec), kz44a(nratec) - & , kz45a(nratec), kz46a(nratec), kz47a(nratec) - & , kz48a(nratec), kz49a(nratec), kz50a(nratec) - & , kz51a(nratec), kz52a(nratec), kz53a(nratec) - & , kz54a(nratec) - -! Cloudy cooling data - - integer icmbTfloor, iClHeat - integer*8 priGridRank, priDataSize, - & metGridRank, metDataSize, - & priGridDim(priGridRank), metGridDim(metGridRank) - real*8 priPar1(priGridDim(1)), priPar2(priGridDim(2)), - & priPar3(priGridDim(3)), priPar4(priGridDim(4)), - & priPar5(priGridDim(5)), - & metPar1(metGridDim(1)), metPar2(metGridDim(2)), - & metPar3(metGridDim(3)), metPar4(metGridDim(4)), - & metPar5(metGridDim(5)), - & priCooling(priDataSize), priHeating(priDataSize), - & priMMW(priDataSize), - & metCooling(metDataSize), metHeating(metDataSize) - -! Parameters - -#ifdef GRACKLE_FLOAT_4 - R_PREC, parameter :: tolerance = 1.0e-05_RKIND -#else - R_PREC, parameter :: tolerance = 1.0e-10_RKIND -#endif - - real*8, parameter :: mh = mass_h - real*8, parameter :: pi = pi_val - -! Locals - - integer i, j, k, iter - integer t, dj, dk - real*8 ttmin, dom, energy, comp1, comp2 - real*8 coolunit, dbase1, tbase1, xbase1, chunit, uvel - real*8 heq1, heq2, eqk221, eqk222, eqk131, eqk132, - & eqt1, eqt2, eqtdef, dheq, heq, dlogtem, dx_cgs, - & c_ljeans, min_metallicity - R_PREC factor - -! row temporaries - - integer*8 indixe(in) - real*8 t1(in), t2(in), logtem(in), tdef(in), - & dtit(in), ttot(in), p2d(in), tgas(in), tgasold(in), - & tdust(in), metallicity(in), dust2gas(in), - & rhoH(in), mmw(in), mynh(in), myde(in), gammaha_eff(in), - & gasgr_tdust(in), regr(in), ddom(in), olddtit - -! Rate equation row temporaries - - real*8 HIp(in), HIIp(in), HeIp(in), HeIIp(in), HeIIIp(in), - & HMp(in), H2Ip(in), H2IIp(in), - & dep(in), dedot(in),HIdot(in), dedot_prev(in), - & DIp(in), DIIp(in), HDIp(in), HIdot_prev(in), - & k24shield(in), k25shield(in), k26shield(in), - & k28shield(in), k29shield(in), k30shield(in), - & k31shield(in), - & 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), - & k13dd(in, 14), h2dust(in), - & ncrn(in), ncrd1(in), ncrd2(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 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) - real*8 h2dustSa(nratec * ndratec), h2dustCa(nratec * ndratec) - real*8 grogra(nratec) - 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) -! grain temperature - real*8 tSiM(in), tFeM(in), tMg2SiO4(in) - & , tMgSiO3(in), tFe3O4(in), tAC(in) - & , tSiO2D(in), tMgO(in), tFeS(in) - & , tAl2O3(in) - & , treforg(in), tvolorg(in), tH2Oice(in) - -! Cooling/heating row locals - - real*8 ceHI(in), ceHeI(in), ceHeII(in), - & ciHI(in), ciHeI(in), ciHeIS(in), ciHeII(in), - & reHII(in), reHeII1(in), reHeII2(in), reHeIII(in), - & brem(in), edot(in) - real*8 hyd01k(in), h2k01(in), vibh(in), roth(in), rotl(in), - & gpldl(in), gphdl(in), hdlte(in), hdlow(in), cieco(in) - -! Iteration mask - - MASK_TYPE itmask(in), anydust - MASK_TYPE itmask_tmp(in), itmask_nr(in) - MASK_TYPE itmask_metal(in) - integer itr, imp_eng(in) - -! -!\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\///////////////////////////////// -!======================================================================= - -! Set error indicator - - ierr = 1 - -! Set flag for dust-related options - - if ((idust .gt. 0) .or. (idustall .gt. 0)) then - anydust = MASK_TRUE - else - anydust = MASK_FALSE - endif - -! ignore metal chemistry/cooling below this metallicity - min_metallicity = 1.d-9 / z_solar - -! Set units - - dom = urho*(aye**3)/mh - tbase1 = utim - xbase1 = uxyz/(aye*uaye) ! uxyz is [x]*a = [x]*[a]*a' ' - dbase1 = urho*(aye*uaye)**3 ! urho is [dens]/a^3 = [dens]/([a]*a')^3 ' - coolunit = (uaye**5 * xbase1**2 * mh**2) / (tbase1**3 * dbase1) - uvel = (uxyz/aye) / utim -c chunit = (1.60218e-12_DKIND)/(2._DKIND*uvel*uvel*mh) ! 1 eV per H2 formed - chunit = (1.60218e-12_DKIND)/(uvel*uvel*mh) ! 1 eV per REACTION (Feb 2020, Gen Chiaki) - - dx_cgs = dx * xbase1 - c_ljeans = sqrt((gamma * pi * kboltz) / - & (GravConst * mh * dbase1)) - - dlogtem = (log(temend) - log(temstart))/real(nratec-1, DKIND) - -! Convert densities from comoving to proper - - if (iexpand .eq. 1) then - - factor = real(aye**(-3), RKIND) - - call scale_fields_g( - & ispecies, imetal, idustfield, imchem, - & imabund, idspecies, immulti, igrgr, idsub, factor, - & is, ie, js, je, ks, ke, in, jn, kn, - & d, de, HI, HII, HeI, HeII, HeIII, - & HM, H2I, H2II, DI, DII, HDI, DM, HDII, HeHII, - & metal, dust, - & 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, - & metal_loc, metal_C13, metal_C20, metal_C25, metal_C30, - & metal_F13, metal_F15, metal_F50, metal_F80, - & metal_P170, metal_P200, metal_Y19) - - endif - - call ceiling_species_g(d, de, HI, HII, HeI, HeII, HeIII, - & HM, H2I, H2II, DI, DII, HDI, metal, dust, - & is, ie, js, je, ks, ke, - & in, jn, kn, ispecies, imetal, idustfield - & , DM, HDII, HeHII - & , imabund, imchem, idspecies, immulti - & , 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 - & , metal_loc - & , metal_C13, metal_C20, metal_C25, metal_C30 - & , metal_F13, metal_F15, metal_F50, metal_F80 - & , metal_P170, metal_P200, metal_Y19) - - -! Loop over zones, and do an entire i-column in one go - dk = ke - ks + 1 - dj = je - js + 1 - -! parallelize the k and j loops with OpenMP -! flat j and k loops for better parallelism -#ifdef _OPENMP -! ierr is declared as shared and should be modified with atomic operation -!$omp parallel do schedule(runtime) private( -!$omp& i, j, k, iter, -!$omp& ttmin, energy, comp1, comp2, -!$omp& heq1, heq2, eqk221, eqk222, eqk131, eqk132, -!$omp& eqt1, eqt2, eqtdef, dheq, heq, -!$omp& indixe, -!$omp& t1, t2, logtem, tdef, -!$omp& dtit, ttot, p2d, tgas, tgasold, -!$omp& tdust, metallicity, dust2gas, rhoH, mmw, -!$omp& mynh, myde, gammaha_eff, gasgr_tdust, regr, ddom, -!$omp& olddtit, -!$omp& HIp, HIIp, HeIp, HeIIp, HeIIIp, -!$omp& HMp, H2Ip, H2IIp, -!$omp& dep, dedot,HIdot, dedot_prev, -!$omp& DIp, DIIp, HDIp, HIdot_prev, -!$omp& k24shield, k25shield, k26shield, -!$omp& k28shield, k29shield, k30shield, -!$omp& k31shield, -!$omp& k1 , k2 , k3 , k4 , k5, -!$omp& k6 , k7 , k8 , k9 , k10, -!$omp& k11, k12, k13, k14, k15, -!$omp& k16, k17, k18, k19, k22, -!$omp& k50, k51, k52, k53, k54, -!$omp& k55, k56, k57, k58, k13dd, h2dust, -!$omp& ncrn, ncrd1, ncrd2, -!$omp& ceHI, ceHeI, ceHeII, -!$omp& ciHI, ciHeI, ciHeIS, ciHeII, -!$omp& reHII, reHeII1, reHeII2, reHeIII, -!$omp& brem, edot, -!$omp& hyd01k, h2k01, vibh, roth, rotl, -!$omp& gpldl, gphdl, hdlte, hdlow, cieco, -!$omp& itmask, itmask_metal ) -#endif - do t = 0, dk*dj-1 - k = t/dj + ks+1 - j = mod(t,dj) + js+1 - -! tolerance = 1.0e-06_DKIND * dt - -! Initialize iteration mask to true for all cells. - - do i = is+1, ie+1 - itmask(i) = MASK_TRUE - enddo - -! If we are using coupled radiation with intermediate stepping, -! set iteration mask to include only cells with radiation in the -! intermediate coupled chemistry / energy step - if (iradtrans .eq. 1) then - if (iradcoupled .eq. 1 .and. iradstep .eq. 1) then - do i = is+1, ie+1 - if (kphHI(i,j,k) .gt. 0) then - itmask(i) = MASK_TRUE - else - itmask(i) = MASK_FALSE - endif - enddo - endif - -! Normal rate solver, but don't double count cells with radiation - if (iradcoupled .eq. 1 .and. iradstep .eq. 0) then - do i = is+1, ie + 1 - if (kphHI(i,j,k) .gt. 0) then - itmask(i) = MASK_FALSE - else - itmask(i) = MASK_TRUE - endif - enddo - endif - endif ! end rad trans check (divergent from original code) - -! Set time elapsed to zero for each cell in 1D section - - do i = is+1, ie+1 - ttot(i) = 0._DKIND - enddo - -! A useful slice variable since we do this a lot - - do i = is+1, ie+1 - ddom(i) = d(i,j,k) * dom - enddo - -! ------------------ Loop over subcycles ---------------- - - do iter = 1, itmax - - do i = is+1, ie+1 - if (itmask(i) .ne. MASK_FALSE) then - dtit(i) = huge8 - endif - enddo - -! Compute the cooling rate, tgas, tdust, and metallicity for this row - - call cool1d_multi_g( - & d, e, u, v, w, de, HI, HII, HeI, HeII, HeIII, - & in, jn, kn, nratec, - & iexpand, ispecies, imetal, imcool, - & idust, idustall, idustfield, idustrec, - & idim, is, ie, j, k, ih2co, ipiht, iter, igammah, - & aye, temstart, temend, z_solar, fgr, - & utem, uxyz, uaye, urho, utim, - & gamma, fh, - & ceHIa, ceHeIa, ceHeIIa, ciHIa, ciHeIa, - & ciHeISa, ciHeIIa, reHIIa, reHeII1a, - & reHeII2a, reHeIIIa, brema, compa, gammaha, - & isrf, regra, gamma_isrfa, comp_xraya, comp_temp, - & piHI, piHeI, piHeII, comp1, comp2, - & HM, H2I, H2II, DI, DII, HDI, metal, dust, - & hyd01ka, h2k01a, vibha, rotha, rotla, - & hyd01k, h2k01, vibh, roth, rotl, - & gpldla, gphdla, gpldl, gphdl, - & hdltea, hdlowa, hdlte, hdlow, - & gaHIa, gaH2a, gaHea, gaHpa, gaela, - & h2ltea, gasgra, - & ceHI, ceHeI, ceHeII, ciHI, ciHeI, ciHeIS, ciHeII, - & reHII, reHeII1, reHeII2, reHeIII, brem, - & indixe, t1, t2, logtem, tdef, edot, - & tgas, tgasold, mmw, p2d, tdust, metallicity, - & dust2gas, rhoH, mynh, myde, - & gammaha_eff, gasgr_tdust, regr, - & iradshield, avgsighi, avgsighei, avgsigheii, - & k24, k26, - & iradtrans, photogamma, - & ih2optical, iciecool, ih2cr, ihdcr, ciecoa, cieco, - & icmbTfloor, iClHeat, clEleFra, - & priGridRank, priGridDim, - & priPar1, priPar2, priPar3, priPar4, priPar5, - & priDataSize, priCooling, priHeating, priMMW, - & metGridRank, metGridDim, - & metPar1, metPar2, metPar3, metPar4, metPar5, - & metDataSize, metCooling, metHeating, clnew, - & iVheat, iMheat, Vheat, Mheat, - & iTfloor, Tfloor_scalar, Tfloor, - & iisrffield, isrf_habing, itmask, itmask_metal - & , imchem, igrgr, ipcont, tmcool - & , DM, HDII, HeHII - & , 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 - & , cieY06a - & , LH2_N, LH2_Size - & , LH2_D, LH2_T, LH2_H - & , LH2_dD, LH2_dT, LH2_dH, LH2_L - & , LHD_N, LHD_Size - & , LHD_D, LHD_T, LHD_H - & , LHD_dD, LHD_dT, LHD_dH, LHD_L - & , LCI_N, LCI_Size - & , LCI_D, LCI_T, LCI_H - & , LCI_dD, LCI_dT, LCI_dH, LCI_L - & , LCII_N, LCII_Size - & , LCII_D, LCII_T, LCII_H - & , LCII_dD, LCII_dT, LCII_dH, LCII_L - & , LOI_N, LOI_Size - & , LOI_D, LOI_T, LOI_H - & , LOI_dD, LOI_dT, LOI_dH, LOI_L - & , LCO_N, LCO_Size - & , LCO_D, LCO_T, LCO_H - & , LCO_dD, LCO_dT, LCO_dH, LCO_L - & , LOH_N, LOH_Size - & , LOH_D, LOH_T, LOH_H - & , LOH_dD, LOH_dT, LOH_dH, LOH_L - & , LH2O_N, LH2O_Size - & , LH2O_D, LH2O_T, LH2O_H - & , LH2O_dD, LH2O_dT, LH2O_dH, LH2O_L - & , alphap_N, alphap_Size - & , alphap_D, alphap_T, alphap_dD, alphap_dT - & , alphap_Data - & , immulti, imabund, idspecies, itdmulti, idsub - & , metal_loc - & , metal_C13, metal_C20, metal_C25, metal_C30 - & , metal_F13, metal_F15, metal_F50, metal_F80 - & , metal_P170, metal_P200, metal_Y19 - & , SN0_N - & , SN0_fSiM, SN0_fFeM, SN0_fMg2SiO4, SN0_fMgSiO3 - & , SN0_fFe3O4, SN0_fAC, SN0_fSiO2D, SN0_fMgO - & , SN0_fFeS, SN0_fAl2O3 - & , SN0_freforg, SN0_fvolorg, SN0_fH2Oice - & , SN0_r0SiM, SN0_r0FeM, SN0_r0Mg2SiO4, SN0_r0MgSiO3 - & , SN0_r0Fe3O4, SN0_r0AC, SN0_r0SiO2D, SN0_r0MgO - & , SN0_r0FeS, SN0_r0Al2O3 - & , SN0_r0reforg, SN0_r0volorg, SN0_r0H2Oice - & , gr_N, gr_Size, gr_dT, gr_Td - & , SN0_kpSiM, SN0_kpFeM, SN0_kpMg2SiO4, SN0_kpMgSiO3 - & , SN0_kpFe3O4, SN0_kpAC, SN0_kpSiO2D, SN0_kpMgO - & , SN0_kpFeS, SN0_kpAl2O3 - & , SN0_kpreforg, SN0_kpvolorg, SN0_kpH2Oice - & , tSiM, tFeM, tMg2SiO4, tMgSiO3, tFe3O4 - & , tAC, tSiO2D, tMgO, tFeS, tAl2O3 - & , treforg, tvolorg, tH2Oice - & , gasgr2a, gamma_isrf2a - & ) - - if (ispecies .eq. 0) then -! This is some basic book-keeping to ensure that itmask_tmp has -! sensible values when ispecies is 0 -! -> There's room for optimization: when ispecies is 0, there is -! need to ever touch this variable. But, for now we focus on -! correct behavior before implementing this optimization - itmask_tmp = itmask - - else - -! Look-up rates as a function of temperature for 1D set of zones -! (maybe should add itmask to this call) - - call lookup_cool_rates1d_g(temstart, temend, nratec, j, k, - & is, ie, ithreebody, - & in, jn, kn, ispecies, anydust, - & iH2shield, iradshield, - & tgas, mmw, d, HI, HII, HeI, HeII, HeIII, - & HM, H2I, H2II, DI, DII, HDI, - & tdust, dust2gas, - & k1a, k2a, k3a, k4a, k5a, k6a, k7a, k8a, k9a, k10a, - & k11a, k12a, k13a, k13dda, k14a, k15a, k16a, - & k17a, k18a, k19a, k22a, - & k50a, k51a, k52a, k53a, k54a, k55a, k56a, - & k57a, k58a, ndratec, dtemstart, dtemend, h2dusta, - & ncrna, ncrd1a, ncrd2a, - & avgsighi, avgsighei, avgsigheii, piHI, piHeI, - & k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, - & k11, k12, k13, k14, k15, k16, k17, k18, - & k19, k22, k24, k25, k26, k28, k29, k30, k31, - & k50, k51, k52, k53, k54, k55, k56, k57, - & k58, k13dd, k24shield, k25shield, k26shield, - & k28shield, k29shield, k30shield, - & k31shield, h2dust, ncrn, ncrd1, ncrd2, - & t1, t2, tdef, logtem, indixe, - & dom, coolunit, tbase1, uxyz, xbase1, dx_cgs, - & c_ljeans, - & iradtrans, kdissH2I, xH2shield, itmask, - & itmask_metal - & , fh, metal - & , DM, HDII, HeHII, imetal, imchem, igrgr - & , 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 - & , k125a, k129a, k130a, k131a, k132a - & , k133a, k134a, k135a, k136a, k137a - & , k148a, k149a, k150a, k151a, k152a - & , k153a - & , kz15a, kz16a, kz17a, kz18a, kz19a - & , kz20a, kz21a, kz22a, kz23a, kz24a - & , kz25a, kz26a, kz27a, kz28a, kz29a - & , kz30a, kz31a, kz32a, kz33a, kz34a - & , kz35a, kz36a, kz37a, kz38a, kz39a - & , kz40a, kz41a, kz42a, kz43a, kz44a - & , kz45a, kz46a, kz47a, kz48a, kz49a - & , kz50a, kz51a, kz52a, kz53a, kz54a - & , 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 - & , immulti, imabund, idspecies, itdmulti, idsub - & , metal_loc - & , metal_C13, metal_C20, metal_C25, metal_C30 - & , metal_F13, metal_F15, metal_F50, metal_F80 - & , metal_P170, metal_P200, metal_Y19 - & , SN0_N - & , SN0_fSiM, SN0_fFeM, SN0_fMg2SiO4, SN0_fMgSiO3 - & , SN0_fFe3O4, SN0_fAC, SN0_fSiO2D, SN0_fMgO - & , SN0_fFeS, SN0_fAl2O3 - & , SN0_freforg, SN0_fvolorg, SN0_fH2Oice - & , SN0_r0SiM, SN0_r0FeM, SN0_r0Mg2SiO4, SN0_r0MgSiO3 - & , SN0_r0Fe3O4, SN0_r0AC, SN0_r0SiO2D, SN0_r0MgO - & , SN0_r0FeS, SN0_r0Al2O3 - & , SN0_r0reforg, SN0_r0volorg, SN0_r0H2Oice - & , gr_N, gr_Size, gr_dT, gr_Td - & , SN0_kpSiM, SN0_kpFeM, SN0_kpMg2SiO4, SN0_kpMgSiO3 - & , SN0_kpFe3O4, SN0_kpAC, SN0_kpSiO2D, SN0_kpMgO - & , SN0_kpFeS, SN0_kpAl2O3 - & , SN0_kpreforg, SN0_kpvolorg, SN0_kpH2Oice - & , h2dustSa, h2dustCa, rhoH, grogra, dt - & , kdSiM, kdFeM, kdMg2SiO4 - & , kdMgSiO3, kdFe3O4, kdAC, kdSiO2D, kdMgO, kdFeS - & , kdAl2O3, kdreforg, kdvolorg, kdH2Oice - & , tSiM, tFeM, tMg2SiO4, tMgSiO3, tFe3O4 - & , tAC, tSiO2D, tMgO, tFeS, tAl2O3 - & , treforg, tvolorg, tH2Oice, iuseH2shield - & , iH2shieldcustom, f_shield_custom - & ) - -! Compute dedot and HIdot, the rates of change of de and HI -! (should add itmask to this call) - - call rate_timestep_g( - & dedot, HIdot, ispecies, anydust, - & de, HI, HII, HeI, HeII, HeIII, d, - & HM, H2I, H2II, - & in, jn, kn, is, ie, j, k, - & 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, ncrn, ncrd1, ncrd2, rhoH, - & k24shield, k25shield, k26shield, - & k28shield, k29shield, k30shield, k31shield, - & iradtrans, irt_honly, - & kphHI, kphHeI, kphHeII, - & itmask, edot, chunit, dom, metal - & , HDI, imchem, CI, OI, OH, CO, H2O - & , idissHDI, kdissHDI, iionZ, kphCI, kphOI - & , idissZ, kdissCO, kdissOH, kdissH2O - & ) - -! move itmask temporary array -! then split cells with low densities -! => Gauss-Seidel scheme -! and with high densities -! => Newton-Raphson scheme - - itmask_tmp = itmask - itmask_nr = itmask - do i = is+1, ie+1 - if ( itmask_tmp(i) .ne. MASK_FALSE ) then - - if ( ( (imetal .eq. 0) - & .and. (ddom(i) .lt. 1.d8) ) - & .or. ( (imetal .eq. 1) - & .and. ( ( (metallicity(i) .le. min_metallicity) - & .and. (ddom(i) .lt. 1.d8) ) - & .or. ( (metallicity(i) .gt. min_metallicity) - & .and. (ddom(i) .lt. 1.d6) ) ) ) ) then - itmask_nr(i) = MASK_FALSE - else - itmask(i) = MASK_FALSE - endif - - endif - enddo - - do i = is+1, ie+1 - if (itmask_nr(i) .ne. MASK_FALSE) then - if ( (icool .eq. 1) .and. (ispecies .gt. 1) .and. - & ((ddom(i) .gt. 1.d7) - & .and.(tgas(i) .gt. 1650.d0)) ) then - imp_eng(i) = 1 - else - imp_eng(i) = 0 - endif - endif - enddo - -! Find timestep that keeps relative chemical changes below 10% - - do i = is+1, ie+1 - if (itmask(i) .ne. MASK_FALSE) then -! Bound from below to prevent numerical errors - - if (abs(dedot(i)) .lt. tiny8) - & dedot(i) = min(tiny,de(i,j,k)) - if (abs(HIdot(i)) .lt. tiny8) - & HIdot(i) = min(tiny,HI(i,j,k)) - -! If the net rate is almost perfectly balanced then set -! it to zero (since it is zero to available precision) - - if (min(abs(k1(i)* de(i,j,k)*HI(i,j,k)), - & abs(k2(i)*HII(i,j,k)*de(i,j,k)))/ - & max(abs(dedot(i)),abs(HIdot(i))) .gt. - & 1.0e6_DKIND) then - dedot(i) = tiny8 - HIdot(i) = tiny8 - endif - -! If the iteration count is high then take the smaller of -! the calculated dedot and last time step's actual dedot. -! This is intended to get around the problem of a low -! electron or HI fraction which is in equilibrium with high -! individual terms (which all nearly cancel). - - if (iter .gt. 50) then - dedot(i) = min(abs(dedot(i)), abs(dedot_prev(i))) - HIdot(i) = min(abs(HIdot(i)), abs(HIdot_prev(i))) - endif - -! compute minimum rate timestep - - olddtit = dtit(i) - dtit(i) = min(abs(0.1_DKIND*de(i,j,k)/dedot(i)), - & abs(0.1_DKIND*HI(i,j,k)/HIdot(i)), - & dt-ttot(i), 0.5_DKIND*dt) - - if (ddom(i) .gt. 1.d8 .and. - & edot(i) .gt. 0._DKIND .and. - & ispecies .gt. 1) then -! Equilibrium value for H is: -! H = (-1._DKIND / (4*k22)) * (k13 - sqrt(8 k13 k22 rho + k13^2)) -! We now want this to change by 10% or less, but we're only -! differentiating by dT. We have de/dt. We need dT/de. -! T = (g-1)*p2d*utem/N; tgas == (g-1)(p2d*utem/N) -! dH_eq / dt = (dH_eq/dT) * (dT/de) * (de/dt) -! dH_eq / dT (see above; we can calculate the derivative here) -! dT / de = utem * (gamma - 1._DKIND) / N == tgas / p2d -! de / dt = edot -! Now we use our estimate of dT/de to get the estimated -! difference in the equilibrium - eqt2 = min(log(tgas(i)) + 0.1_DKIND*dlogtem, t2(i)) - eqtdef = (eqt2 - t1(i))/(t2(i) - t1(i)) - eqk222 = k22a(indixe(i)) + - & (k22a(indixe(i)+1) -k22a(indixe(i)))*eqtdef - eqk132 = k13a(indixe(i)) + - & (k13a(indixe(i)+1) -k13a(indixe(i)))*eqtdef - heq2 = (-1._DKIND / (4._DKIND*eqk222)) * (eqk132- - & sqrt(8._DKIND*eqk132*eqk222* - & fh*d(i,j,k)+eqk132**2._DKIND)) - - eqt1 = max(log(tgas(i)) - 0.1_DKIND*dlogtem, t1(i)) - eqtdef = (eqt1 - t1(i))/(t2(i) - t1(i)) - eqk221 = k22a(indixe(i)) + - & (k22a(indixe(i)+1) -k22a(indixe(i)))*eqtdef - eqk131 = k13a(indixe(i)) + - & (k13a(indixe(i)+1) -k13a(indixe(i)))*eqtdef - heq1 = (-1._DKIND / (4._DKIND*eqk221)) * (eqk131- - & sqrt(8._DKIND*eqk131*eqk221* - & fh*d(i,j,k)+eqk131**2._DKIND)) - - dheq = (abs(heq2-heq1)/(exp(eqt2) - exp(eqt1))) - & * (tgas(i)/p2d(i)) * edot(i) - heq = (-1._DKIND / (4._DKIND*k22(i))) * (k13(i)- - & sqrt(8._DKIND*k13(i)*k22(i)* - & fh*d(i,j,k)+k13(i)**2._DKIND)) - dtit(i) = min(dtit(i), 0.1_DKIND*heq/dheq) - endif - if (iter.gt.10_DKIND) then - dtit(i) = min(olddtit*1.5_DKIND, dtit(i)) - endif - - else if ((itmask_nr(i).ne.MASK_FALSE).and. - & (imp_eng(i).eq.0)) then - dtit(i) = min(abs(0.1_DKIND*e(i,j,k)/edot(i)*d(i,j,k)), - & dt-ttot(i), 0.5_DKIND*dt) - - else if ((itmask_nr(i).ne.MASK_FALSE).and. - & (imp_eng(i).eq.1)) then - dtit(i) = dt - ttot(i) - - else ! itmask - dtit(i) = dt - endif - enddo ! end loop over i - - endif ! end if (ispecies .gt. 0) - -! Compute maximum timestep for cooling/heating - - do i = is+1, ie+1 - if (itmask(i) .ne. MASK_FALSE) then -! Set energy per unit volume of this cell based in the pressure -! (the gamma used here is the right one even for H2 since p2d -! is calculated with this gamma). - - energy = max(p2d(i)/(gamma-1._DKIND), tiny8) - -! If the temperature is at the bottom of the temperature look-up -! table and edot < 0, then shut off the cooling. - - if (tgas(i) .le. 1.01_DKIND*temstart .and. - & edot(i) .lt. 0._DKIND) - & edot(i) = tiny8 - if (abs(edot(i)) .lt. tiny8) edot(i) = tiny8 - -! Compute timestep for 10% change - - dtit(i) = min(real(abs(0.1_DKIND* - & energy/edot(i)), DKIND), - & dt-ttot(i), dtit(i)) - - if (dtit(i) .ne. dtit(i)) then -#ifdef _OPENMP -!$omp critical -#endif - write(6,*) 'HUGE dtit :: ', energy, edot(i), dtit(i), - & dt, ttot(i), abs(0.1_DKIND*energy/edot(i)), - & real(abs(0.1_DKIND*energy/edot(i)), DKIND) -#ifdef _OPENMP -!$omp end critical -#endif - endif - - endif ! itmask - enddo ! end loop over i - -! Update total and gas energy - - if (icool .eq. 1) then - do i = is+1, ie+1 - if (itmask(i) .ne. MASK_FALSE) then - e(i,j,k) = e(i,j,k) + - & real(edot(i)/d(i,j,k)*dtit(i), RKIND) - - endif ! itmask - enddo - endif - - if (ispecies .gt. 0) then - -! Solve rate equations with one linearly implicit Gauss-Seidel -! sweep of a backward Euler method (for all cells specified by -! itmask) - - call 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 - & ) - - -! Solve rate equations with one linearly implicit Gauss-Seidel -! sweep of a backward Euler method (for all cells specified by -! itmask_nr) - - call step_rate_newton_raphson(icool, d, e, u, v, w, de, HI, - & HII, HeI, HeII, HeIII, in, jn, kn, nratec, - & iexpand, ispecies, imetal, imcool, idust, - & idustall, idustfield, is, ie, ih2co, ipiht, - & idustrec, igammah, aye, temstart, temend, utem, - & uxyz, uaye, urho, utim, gamma, fh, z_solar, fgr, - & k1a, k2a, k3a, k4a, k5a, k6a, k7a, k8a, k9a, - & k10a, k11a, k12a, k13a, k13dda, k14a, k15a, k16a, - & k17a, k18a, k19a, k22a, k24, k25, k26, k27, k28, - & k29, k30, k31, k50a, k51a, k52a, k53a, k54a, - & k55a, k56a, k57a, k58a, ndratec, dtemstart, - & dtemend, h2dusta, ncrna, ncrd1a, ncrd2a, ceHIa, - & ceHeIa, ceHeIIa, ciHIa, ciHeIa, ciHeISa, ciHeIIa, - & reHIIa, reHeII1a, reHeII2a, reHeIIIa, brema, - & compa, gammaha, isrf, regra, gamma_isrfa, - & comp_xraya, comp_temp, piHI, piHeI, piHeII, HM, - & H2I, H2II, DI, DII, HDI, metal, dust, hyd01ka, - & h2k01a, vibha, rotha, rotla, gpldla, gphdla, - & hdltea, hdlowa, gaHIa, gaH2a, gaHea, gaHpa, - & gaela, h2ltea, gasgra, iH2shield, iradshield, - & avgsighi, avgsighei, avgsigheii, iradtrans, - & irt_honly, kphHI, kphHeI, kphHeII, kdissH2I, - & photogamma, xH2shield, ih2optical, iciecool, - & ithreebody, ih2cr, ihdcr, ciecoa, icmbTfloor, - & iClHeat, clEleFra, priGridRank, priGridDim, - & priPar1, priPar2, priPar3, priPar4, priPar5, - & priDataSize, priCooling, priHeating, priMMW, - & metGridRank, metGridDim, metPar1, metPar2, - & metPar3, metPar4, metPar5, metDataSize, - & metCooling, metHeating, clnew, iVheat, iMheat, - & Vheat, Mheat, iTfloor, Tfloor_scalar, Tfloor, - & imchem, igrgr, ipcont, tmcool, DM, HDII, HeHII, - & 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, k125a, k129a, k130a, k131a, - & k132a, k133a, k134a, k135a, k136a, k137a, k148a, - & k149a, k150a, k151a, k152a, k153a, kz15a, kz16a, - & kz17a, kz18a, kz19a, kz20a, kz21a, kz22a, kz23a, - & kz24a, kz25a, kz26a, kz27a, kz28a, kz29a, kz30a, - & kz31a, kz32a, kz33a, kz34a, kz35a, kz36a, kz37a, - & kz38a, kz39a, kz40a, kz41a, kz42a, kz43a, kz44a, - & kz45a, kz46a, kz47a, kz48a, kz49a, kz50a, kz51a, - & kz52a, kz53a, kz54a, cieY06a, LH2_N, LH2_Size, - & LH2_D, LH2_T, LH2_H, LH2_dD, LH2_dT, LH2_dH, - & LH2_L, LHD_N, LHD_Size, LHD_D, LHD_T, LHD_H, - & LHD_dD, LHD_dT, LHD_dH, LHD_L, LCI_N, LCI_Size, - & LCI_D, LCI_T, LCI_H, LCI_dD, LCI_dT, LCI_dH, - & LCI_L, LCII_N, LCII_Size, LCII_D, LCII_T, LCII_H, - & LCII_dD, LCII_dT, LCII_dH, LCII_L, LOI_N, - & LOI_Size, LOI_D, LOI_T, LOI_H, LOI_dD, LOI_dT, - & LOI_dH, LOI_L, LCO_N, LCO_Size, LCO_D, LCO_T, - & LCO_H, LCO_dD, LCO_dT, LCO_dH, LCO_L, LOH_N, - & LOH_Size, LOH_D, LOH_T, LOH_H, LOH_dD, LOH_dT, - & LOH_dH, LOH_L, LH2O_N, LH2O_Size, LH2O_D, LH2O_T, - & LH2O_H, LH2O_dD, LH2O_dT, LH2O_dH, LH2O_L, - & alphap_N, alphap_Size, alphap_D, alphap_T, - & alphap_dD, alphap_dT, alphap_Data, immulti, - & imabund, idspecies, itdmulti, idsub, metal_loc, - & metal_C13, metal_C20, metal_C25, metal_C30, - & metal_F13, metal_F15, metal_F50, metal_F80, - & metal_P170, metal_P200, metal_Y19, SN0_N, - & SN0_fSiM, SN0_fFeM, SN0_fMg2SiO4, SN0_fMgSiO3, - & SN0_fFe3O4, SN0_fAC, SN0_fSiO2D, SN0_fMgO, - & SN0_fFeS, SN0_fAl2O3, SN0_freforg, SN0_fvolorg, - & SN0_fH2Oice, SN0_r0SiM, SN0_r0FeM, SN0_r0Mg2SiO4, - & SN0_r0MgSiO3, SN0_r0Fe3O4, SN0_r0AC, SN0_r0SiO2D, - & SN0_r0MgO, SN0_r0FeS, SN0_r0Al2O3, SN0_r0reforg, - & SN0_r0volorg, SN0_r0H2Oice, gr_N, gr_Size, gr_dT, - & gr_Td, SN0_kpSiM, SN0_kpFeM, SN0_kpMg2SiO4, - & SN0_kpMgSiO3, SN0_kpFe3O4, SN0_kpAC, SN0_kpSiO2D, - & SN0_kpMgO, SN0_kpFeS, SN0_kpAl2O3, SN0_kpreforg, - & SN0_kpvolorg, SN0_kpH2Oice, h2dustSa, h2dustCa, - & gasgr2a, gamma_isrf2a, grogra, idissHDI, - & kdissHDI, iionZ, kphCI, kphOI, idissZ, kdissCO, - & kdissOH, kdissH2O, iuseH2shield, iisrffield, - & isrf_habing, iH2shieldcustom, f_shield_custom, - & ierror, j, k, iter, dom, comp1, - & comp2, coolunit, tbase1, xbase1, chunit, dx_cgs, - & c_ljeans, indixe, t1, t2, logtem, tdef, dtit, - & p2d, tgas, tgasold, tdust, metallicity, dust2gas, - & rhoH, mmw, mynh, myde, gammaha_eff, gasgr_tdust, - & regr, h2dust, ncrn, ncrd1, ncrd2, tSiM, tFeM, - & tMg2SiO4, tMgSiO3, tFe3O4, tAC, tSiO2D, tMgO, - & tFeS, tAl2O3, treforg, tvolorg, tH2Oice, ceHI, - & ceHeI, ceHeII, ciHI, ciHeI, ciHeIS, ciHeII, - & reHII, reHeII1, reHeII2, reHeIII, brem, edot, - & hyd01k, h2k01, vibh, roth, rotl, gpldl, gphdl, - & hdlte, hdlow, cieco, anydust, itmask_nr, - & itmask_metal, itr, imp_eng - & ) - - endif ! if (ispecies .gt. 0) then - -! return itmask - do i = is+1, ie+1 - itmask(i) = itmask_tmp(i) - enddo - -! Add the timestep to the elapsed time for each cell and find -! minimum elapsed time step in this row - - ttmin = huge8 - do i = is+1, ie+1 - ttot(i) = min(ttot(i) + dtit(i), dt) - if (abs(dt-ttot(i)) .lt. - & tolerance*dt) itmask(i) = MASK_FALSE - if (ttot(i).lt.ttmin) ttmin = ttot(i) - enddo - -! If all cells are done (on this slice), then exit - - if (abs(dt-ttmin) .lt. tolerance*dt) go to 9999 - -! Next subcycle iteration - - enddo ! end of 'do iter = 1, itmax' - - 9999 continue - -! Abort if iteration count exceeds maximum - - if (iter .gt. itmax) then -#ifdef _OPENMP -!$omp critical -#endif - write(0,*) 'inside if statement solve rate cool:',is,ie - write(6,*) 'MULTI_COOL iter > ',itmax,' at j,k =',j,k - write(0,*) 'FATAL error (2) in MULTI_COOL' - write(0,'(" dt = ",1pe10.3," ttmin = ",1pe10.3)') dt, ttmin - write(0,'((16(1pe8.1)))') (dtit(i),i=is+1,ie+1) - write(0,'((16(1pe8.1)))') (ttot(i),i=is+1,ie+1) - write(0,'((16(1pe8.1)))') (edot(i),i=is+1,ie+1) - write(0,'((16(I3)))') (itmask(i),i=is+1,ie+1) - - if (exititmax .eq. 1) then - ierr = 0 - endif -#ifdef _OPENMP -!$omp end critical -#endif -c WARNING_MESSAGE - endif - - if (iter .gt. itmax/2) then -#ifdef _OPENMP -!$omp critical -#endif - write(6,*) 'MULTI_COOL iter,j,k =',iter,j,k -#ifdef _OPENMP -!$omp end critical -#endif - end if -! -! Next j,k -! - enddo ! end of 'do t = 0, dk*dj-1' -#ifdef _OPENMP -!$omp end parallel do -#endif - -! If an error has been produced, return now. - - if (ierr .eq. 0) then - return - endif - -! Convert densities back to comoving from proper - - if (iexpand .eq. 1) then - - factor = real(aye**3, RKIND) - - call scale_fields_g( - & ispecies, imetal, idustfield, imchem, - & imabund, idspecies, immulti, igrgr, idsub, factor, - & is, ie, js, je, ks, ke, in, jn, kn, - & d, de, HI, HII, HeI, HeII, HeIII, - & HM, H2I, H2II, DI, DII, HDI, DM, HDII, HeHII, - & metal, dust, - & 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, - & metal_loc, metal_C13, metal_C20, metal_C25, metal_C30, - & metal_F13, metal_F15, metal_F50, metal_F80, - & metal_P170, metal_P200, metal_Y19) - - endif - - if (ispecies .gt. 0) then - -! Correct the species to ensure consistency (i.e. type conservation) - - call make_consistent_g(de, HI, HII, HeI, HeII, HeIII, - & HM, H2I, H2II, DI, DII, HDI, metal, dust, - & d, is, ie, js, je, ks, ke, - & in, jn, kn, ispecies, imetal, fh, dtoh - & , idustfield, imchem, igrgr, dom - & , DM, HDII, HeHII - & , 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 - & , immulti, imabund, idspecies, itdmulti, idsub - & , metal_loc - & , metal_C13, metal_C20, metal_C25, metal_C30 - & , metal_F13, metal_F15, metal_F50, metal_F80 - & , metal_P170, metal_P200, metal_Y19 - & , SN0_N - & , SN0_XC, SN0_XO, SN0_XMg, SN0_XAl, SN0_XSi - & , SN0_XS, SN0_XFe - & , SN0_fC, SN0_fO, SN0_fMg, SN0_fAl, SN0_fSi - & , SN0_fS, SN0_fFe - & ) - - endif - - return - end +! This file originally defined the solve_rate_cool_g subroutine and a +! number of helper subroutines +! - solve_rate_cool_g has now been transcribed to C++ +! - the remaining helper subroutines still need to be transcribed c ----------------------------------------------------------- ! This routine ensures that the species aren't below tiny. diff --git a/src/clib/step_rate_newton_raphson.F b/src/clib/step_rate_newton_raphson.F index 9aff620e1..cf0b4273d 100644 --- a/src/clib/step_rate_newton_raphson.F +++ b/src/clib/step_rate_newton_raphson.F @@ -76,7 +76,7 @@ subroutine step_rate_newton_raphson(icool, d, e, u, v, w, de, HI, & kdissHDI, iionZ, kphCI, kphOI, idissZ, kdissCO, & kdissOH, kdissH2O, iuseH2shield, iisrffield, & isrf_habing, iH2shieldcustom, f_shield_custom, - & ierror, j, k, iter, dom, comp1, + & j, k, iter, dom, comp1, & comp2, coolunit, tbase1, xbase1, chunit, dx_cgs, & c_ljeans, indixe, t1, t2, logtem, tdef, dtit, & p2d, tgas, tgasold, tdust, metallicity, dust2gas, @@ -88,7 +88,7 @@ subroutine step_rate_newton_raphson(icool, d, e, u, v, w, de, HI, & reHII, reHeII1, reHeII2, reHeIII, brem, edot, & hyd01k, h2k01, vibh, roth, rotl, gpldl, gphdl, & hdlte, hdlow, cieco, anydust, itmask_nr, - & itmask_metal, itr, imp_eng + & itmask_metal, imp_eng & ) ! PURPOSE: @@ -101,7 +101,7 @@ subroutine step_rate_newton_raphson(icool, d, e, u, v, w, de, HI, ! General Arguments integer icool, in, jn, kn, is, ie, nratec, iexpand, ih2co, ipiht, - & ispecies, imetal, ierror, imcool, idust, idustall, + & ispecies, imetal, imcool, idust, idustall, & idustfield, idustrec, igammah, ih2optical, iciecool, & ithreebody, ih2cr, ihdcr, ndratec, clnew, iVheat, iMheat, & iTfloor, iH2shield, iradshield, iradtrans, irt_honly, @@ -330,9 +330,16 @@ subroutine step_rate_newton_raphson(icool, d, e, u, v, w, de, HI, MASK_TYPE anydust MASK_TYPE itmask_nr(in) MASK_TYPE itmask_metal(in) - integer itr, imp_eng(in) + integer imp_eng(in) + +! ierror local variable +! - this variable is only used internally by this subroutine for +! determining control-flow +! - it does NOT report whether or not this function succeeded (maybe +! we should change the name?) + integer ierror ! Local variable - integer itr_time + integer itr, itr_time integer nsp, isp, jsp, id real*8 dspj, err, err_max integer,parameter :: i_eng = 52