diff --git a/src/clib/CMakeLists.txt b/src/clib/CMakeLists.txt index 6e3b757c8..1e3f4fc8e 100644 --- a/src/clib/CMakeLists.txt +++ b/src/clib/CMakeLists.txt @@ -109,6 +109,7 @@ add_library(Grackle_Grackle calc_temp_cloudy_g.cpp calc_temp_cloudy_g.h ceiling_species.hpp collisional_rate_props.cpp collisional_rate_props.hpp + cool1d_cloudy_old_tables_g.cpp cool1d_cloudy_old_tables_g.hpp cool1d_multi_g.cpp cool1d_multi_g.hpp cool_multi_time_g.cpp cool_multi_time_g.h dust_props.hpp @@ -138,7 +139,6 @@ add_library(Grackle_Grackle calc_tdust_1d_g.F calc_temp1d_cloudy_g.F cool1d_cloudy_g.F - cool1d_cloudy_old_tables_g.F interpolators_g.F solve_rate_cool_g.F calc_grain_size_increment_1d.F diff --git a/src/clib/cool1d_cloudy_old_tables_g.F b/src/clib/cool1d_cloudy_old_tables_g.F deleted file mode 100644 index e7ccbb114..000000000 --- a/src/clib/cool1d_cloudy_old_tables_g.F +++ /dev/null @@ -1,335 +0,0 @@ -!======================================================================= -!////////////// SUBROUTINE COOL1D_CLOUDY_OLD_TABLES_G \\\\\\\\\\\\\\\\ - - subroutine cool1D_cloudy_old_tables_g(d, de, rhoH, metallicity, - & in, jn, kn, is, ie, j, k, - & logtem, edot, comp2, ispecies, dom, zr, - & icmbTfloor, iClHeat, - & clEleFra, clGridRank, clGridDim, - & clPar1, clPar2, clPar3, clPar4, clPar5, - & clDataSize, clCooling, clHeating, - & itmask) - -! -! SOLVE CLOUDY METAL COOLING FOR OLD CLOUDY TABLES -! -! written by: Britton Smith -! date: September, 2009 -! -! PURPOSE: -! Solve cloudy cooling by interpolating from the data. -! This version uses tables formatted for the original -! Cloudy cooling functionality in Enzo. -! -! INPUTS: -! in,jn,kn - dimensions of 3D fields -! -! d - total density field -! de - electron density field -! -! rhoH - total H mass density -! metallicity - metallicity -! -! is,ie - start and end indices of active region (zero based) -! ispecies - chemistry module (1 - H/He only, 2 - molecular H, 3 - D) -! logtem - natural log of temperature values -! -! dom - unit conversion to proper number density in code units -! zr - current redshift -! -! icmbTfloor - flag to include temperature floor from cmb -! iClHeat - flag to include cloudy heating -! clEleFra - parameter to account for additional electrons from metals -! clGridRank - rank of cloudy cooling data grid -! clGridDim - array containing dimensions of cloudy data -! clPar1, clPar2, clPar3, clPar4, clPar5 - arrays containing cloudy grid parameter values -! clDataSize - total size of flattened 1D cooling data array -! clCooling - cloudy cooling data -! clHeating - cloudy heating data -! -! itmask - iteration mask -! -! OUTPUTS: -! update edot with heating/cooling contributions from metals -! -! PARAMETERS: -! -!----------------------------------------------------------------------- - - implicit NONE -#include "grackle_fortran_types.def" - -! General Arguments - - integer in, jn, kn, is, ie, j, k, ispecies - - real*8 comp2, dom, zr - R_PREC d(in,jn,kn), de(in,jn,kn) - real*8 rhoH(in), metallicity(in), - & logtem(in), edot(in) - -! Cloudy parameters and data - - integer icmbTfloor, iClHeat - integer*8 clGridRank, clDataSize, clGridDim(5) - real*8 clEleFra - real*8 clPar1(clGridDim(1)), clPar2(clGridDim(2)), - & clPar3(clGridDim(3)), clPar4(clGridDim(4)), - & clPar5(clGridDim(5)) - real*8 clCooling(clDataSize), clHeating(clDataSize) - -! Iteration mask - - MASK_TYPE itmask(in) - -! Parameters - -! Locals - - integer i, q - real*8 dclPar(clGridRank), inv_log10, log10_tCMB - -! Slice locals - - real*8 log_Z(in), e_frac(in), log_e_frac(in), - & cl_e_frac(in), fh(in), log_n_h(in), - & log_cool(in), log_cool_cmb(in), log_heat(in), - & edot_met(in), log10tem(in) - -!\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\///////////////////////////////// -!======================================================================= - - inv_log10 = 1._DKIND / log(10._DKIND) - log10_tCMB = log10(comp2) - -! Calculate parameter value slopes - - dclPar(1) = (clPar1(clGridDim(1)) - clPar1(1)) / - & real(clGridDim(1) - 1, DKIND) - if (clGridRank > 1) then - dclPar(2) = (clPar2(clGridDim(2)) - clPar2(1)) / - & real(clGridDim(2) - 1, DKIND) - endif - if (clGridRank > 2) then - dclPar(3) = (clPar3(clGridDim(3)) - clPar3(1)) / - & real(clGridDim(3) - 1, DKIND) - endif - if (clGridRank > 3) then - dclPar(4) = (clPar4(clGridDim(4)) - clPar4(1)) / - & real(clGridDim(4) - 1, DKIND) - endif - if (clGridRank > 4) then - dclPar(5) = (clPar5(clGridDim(5)) - clPar5(1)) / - & real(clGridDim(5) - 1, DKIND) - endif - - do i=is+1, ie+1 - if ( itmask(i) .ne. MASK_FALSE ) then - - log10tem(i) = logtem(i) * inv_log10 - -! Calcualte H mass fraction - - fh(i) = rhoH(i) / d(i,j,k) - -! Calculate proper log(n_H) - - if (clGridRank > 1) then - - log_n_h(i) = log10(rhoH(i) * dom) - - endif - -! Calculate metallicity - - if (clGridRank > 2) then - - log_Z(i) = log10(metallicity(i)) - - endif - -! Calculate electron fraction - - if (clGridRank > 3) then - - e_frac(i) = 2._DKIND * de(i,j,k) / - & (d(i,j,k) * (1._DKIND + fh(i))) -! Make sure electron fraction is never above 1 -! which can give bad cooling/heating values when -! extrapolating in the Cloudy data. - log_e_frac(i) = min(log10(e_frac(i)), 0._DKIND) - -! Get extra electrons contributed by metals - - cl_e_frac(i) = e_frac(i) * - & (1._DKIND + (2._DKIND * clEleFra * metallicity(i) * - & fh(i)) / (1._DKIND + fh(i))) - - endif - -! Call interpolation functions to get heating/cooling - -! Interpolate over temperature. - if (clGridRank == 1) then - call interpolate_1D_g( - & log10tem(i), clGridDim, clPar1, - & dclPar(1), clDataSize, clCooling, log_cool(i)) - edot_met(i) = -10._DKIND**log_cool(i) - -! Ignore CMB term if T >> T_CMB - if ((icmbTfloor == 1) .and. - & ((log10tem(i) - log10_tCMB) < 2._DKIND)) then - call interpolate_1D_g( - & log10_tCMB, clGridDim, clPar1, - & dclPar(1), clDataSize, clCooling, - & log_cool_cmb(i)) - edot_met(i) = edot_met(i) + 10._DKIND**log_cool_cmb(i) - endif - - if (iClHeat == 1) then - call interpolate_1D_g( - & log10tem(i), clGridDim, clPar1, - & dclPar(1), clDataSize, clHeating, - & log_heat(i)) - edot_met(i) = edot_met(i) + 10._DKIND**log_heat(i) - endif - -! Interpolate over density and temperature. - else if (clGridRank == 2) then - call interpolate_2D_g( - & log_n_h(i), log10tem(i), clGridDim, - & clPar1, dclPar(1), clPar2, dclPar(2), - & clDataSize, clCooling, log_cool(i)) - edot_met(i) = -10._DKIND**log_cool(i) - -! Ignore CMB term if T >> T_CMB - if ((icmbTfloor == 1) .and. - & ((log10tem(i) - log10_tCMB) < 2.0)) then - call interpolate_2D_g( - & log_n_h(i), log10_tCMB, clGridDim, - & clPar1, dclPar(1), clPar2, dclPar(2), - & clDataSize, clCooling, log_cool_cmb(i)) - edot_met(i) = edot_met(i) + 10._DKIND**log_cool_cmb(i) - endif - - if (iClHeat == 1) then - call interpolate_2D_g( - & log_n_h(i), log10tem(i), clGridDim, - & clPar1, dclPar(1), clPar2, dclPar(2), - & clDataSize, clHeating, log_heat(i)) - edot_met(i) = edot_met(i) + 10._DKIND**log_heat(i) - endif - -! Interpolate over density, metallicity, and temperature. - else if (clGridRank == 3) then - call interpolate_3D_g( - & log_n_h(i), log_Z(i), log10tem(i), - & clGridDim, - & clPar1, dclPar(1), clPar2, dclPar(2), - & clPar3, dclPar(3), - & clDataSize, clCooling, log_cool(i)) - edot_met(i) = -10._DKIND**log_cool(i) - -! Ignore CMB term if T >> T_CMB - if ((icmbTfloor == 1) .and. - & ((log10tem(i) - log10_tCMB) < 2._DKIND)) then - call interpolate_3D_g( - & log_n_h(i), log_Z(i), log10_tCMB, - & clGridDim, - & clPar1, dclPar(1), clPar2, dclPar(2), - & clPar3, dclPar(3), - & clDataSize, clCooling, log_cool_cmb(i)) - edot_met(i) = edot_met(i) + 10._DKIND**log_cool_cmb(i) - endif - - if (iClHeat == 1) then - call interpolate_3D_g( - & log_n_h(i), log_Z(i), log10tem(i), - & clGridDim, - & clPar1, dclPar(1), clPar2, dclPar(2), - & clPar3, dclPar(3), - & clDataSize, clHeating, log_heat(i)) - edot_met(i) = edot_met(i) + 10._DKIND**log_heat(i) - endif - -! Interpolate over density, metallicity, electron fraction, and temperature. - else if (clGridRank == 4) then - call interpolate_4D_g( - & log_n_h(i), log_Z(i), - & log_e_frac(i), log10tem(i), clGridDim, - & clPar1, dclPar(1), clPar2, dclPar(2), - & clPar3, dclPar(3), clPar4, dclPar(4), - & clDataSize, clCooling, log_cool(i)) - edot_met(i) = -10._DKIND**log_cool(i) - -! Ignore CMB term if T >> T_CMB - if ((icmbTfloor == 1) .and. - & ((log10tem(i) - log10_tCMB) < 2._DKIND)) then - call interpolate_4D_g( - & log_n_h(i), log_Z(i), - & log_e_frac(i), log10_tCMB, clGridDim, - & clPar1, dclPar(1), clPar2, dclPar(2), - & clPar3, dclPar(3), clPar4, dclPar(4), - & clDataSize, clCooling, log_cool_cmb(i)) - edot_met(i) = edot_met(i) + 10._DKIND**log_cool_cmb(i) - endif - - if (iClHeat == 1) then - call interpolate_4D_g( - & log_n_h(i), log_Z(i), - & log_e_frac(i), log10tem(i), clGridDim, - & clPar1, dclPar(1), clPar2, dclPar(2), - & clPar3, dclPar(3), clPar4, dclPar(4), - & clDataSize, clHeating, log_heat(i)) - edot_met(i) = edot_met(i) + 10._DKIND**log_heat(i) - endif - -! Interpolate over density, metallicity, electron fraction, redshift, -! and temperature. - else - call interpolate_5D_g( - & log_n_h(i), log_Z(i), - & log_e_frac(i), zr, log10tem(i), clGridDim, - & clPar1, dclPar(1), clPar2, dclPar(2), - & clPar3, dclPar(3), clPar4, dclPar(4), - & clPar5, dclPar(5), - & clDataSize, clCooling, log_cool(i)) - edot_met(i) = -10._DKIND**log_cool(i) - -! Ignore CMB term if T >> T_CMB - if ((icmbTfloor == 1) .and. - & ((log10tem(i) - log10_tCMB) < 2._DKIND)) then - call interpolate_5D_g( - & log_n_h(i), log_Z(i), - & log_e_frac(i), zr, log10_tCMB, clGridDim, - & clPar1, dclPar(1), clPar2, dclPar(2), - & clPar3, dclPar(3), clPar4, dclPar(4), - & clPar5, dclPar(5), - & clDataSize, clCooling, log_cool_cmb(i)) - edot_met(i) = edot_met(i) + 10._DKIND**log_cool_cmb(i) - endif - - if (iClHeat == 1) then - call interpolate_5D_g( - & log_n_h(i), log_Z(i), - & log_e_frac(i), zr, log10tem(i), clGridDim, - & clPar1, dclPar(1), clPar2, dclPar(2), - & clPar3, dclPar(3), clPar4, dclPar(4), - & clPar5, dclPar(5), - & clDataSize, clHeating, log_heat(i)) - edot_met(i) = edot_met(i) + 10._DKIND**log_heat(i) - endif - - endif - - if (clGridRank > 3) then - edot_met(i) = edot_met(i) * cl_e_frac(i) - endif - - edot(i) = edot(i) + (edot_met(i) * rhoH(i) * d(i,j,k)) - - end if - enddo - - return - end diff --git a/src/clib/cool1d_cloudy_old_tables_g.cpp b/src/clib/cool1d_cloudy_old_tables_g.cpp new file mode 100644 index 000000000..97ba82461 --- /dev/null +++ b/src/clib/cool1d_cloudy_old_tables_g.cpp @@ -0,0 +1,305 @@ +// See LICENSE file for license and copyright information + +/// @file cool1d_cloudy_old_tables_g-cpp.C +/// @brief Declares signature of cool1d_cloudy_old_tables_g + +// This file was initially generated automatically during conversion of the +// cool1d_cloudy_old_tables_g function from FORTRAN to C++ + +#include +#include +#include + +#include "grackle.h" +#include "fortran_func_decls.h" +#include "fortran_func_wrappers.hpp" +#include "utils-cpp.hpp" + +#include "cool1d_cloudy_old_tables_g.hpp" + +void grackle::impl::cool1d_cloudy_old_tables_g( + const double* rhoH, double* metallicity, const double* logtem, double* edot, + double comp2, double dom, double zr, const gr_mask_type* itmask, + chemistry_data* my_chemistry, cloudy_data cloudy_table, gr_float* density, + gr_float* e_density, grackle_field_data* my_fields, IndexRange idx_range) { + // General Arguments + + grackle::impl::View d(density, idx_range.i_stop, + my_fields->grid_dimension[1], + my_fields->grid_dimension[2]); + grackle::impl::View de(e_density, idx_range.i_stop, + my_fields->grid_dimension[1], + my_fields->grid_dimension[2]); + + // Locals + + int i; + double inv_log10, log10_tCMB; + std::vector dclPar(cloudy_table.grid_rank); + + // Slice locals + + std::vector log_Z(idx_range.i_stop); + std::vector e_frac(idx_range.i_stop); + std::vector log_e_frac(idx_range.i_stop); + std::vector cl_e_frac(idx_range.i_stop); + std::vector fh(idx_range.i_stop); + std::vector log_n_h(idx_range.i_stop); + std::vector log_cool(idx_range.i_stop); + std::vector log_cool_cmb(idx_range.i_stop); + std::vector log_heat(idx_range.i_stop); + std::vector edot_met(idx_range.i_stop); + std::vector log10tem(idx_range.i_stop); + + // \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\///////////////////////////////// + // ======================================================================= + + inv_log10 = 1. / std::log(10.); + log10_tCMB = std::log10(comp2); + + // Calculate parameter value slopes + + dclPar[0] = + (cloudy_table.grid_parameters[0][cloudy_table.grid_dimension[0] - 1] - + cloudy_table.grid_parameters[0][0]) / + (double)(cloudy_table.grid_dimension[0] - 1); + if (cloudy_table.grid_rank > 1) { + dclPar[1] = + (cloudy_table.grid_parameters[1][cloudy_table.grid_dimension[1] - 1] - + cloudy_table.grid_parameters[1][0]) / + (double)(cloudy_table.grid_dimension[1] - 1); + } + if (cloudy_table.grid_rank > 2) { + dclPar[2] = + (cloudy_table.grid_parameters[2][cloudy_table.grid_dimension[2] - 1] - + cloudy_table.grid_parameters[2][0]) / + (double)(cloudy_table.grid_dimension[2] - 1); + } + if (cloudy_table.grid_rank > 3) { + dclPar[3] = + (cloudy_table.grid_parameters[3][cloudy_table.grid_dimension[3] - 1] - + cloudy_table.grid_parameters[3][0]) / + (double)(cloudy_table.grid_dimension[3] - 1); + } + if (cloudy_table.grid_rank > 4) { + dclPar[4] = + (cloudy_table.grid_parameters[4][cloudy_table.grid_dimension[4] - 1] - + cloudy_table.grid_parameters[4][0]) / + (double)(cloudy_table.grid_dimension[4] - 1); + } + + for (i = idx_range.i_start; i <= idx_range.i_end; i++) { + if (itmask[i] != MASK_FALSE) { + log10tem[i] = logtem[i] * inv_log10; + + // Calcualte H mass fraction + + fh[i] = rhoH[i] / d(i, idx_range.j, idx_range.k); + + // Calculate proper log(n_H) + + if (cloudy_table.grid_rank > 1) { + log_n_h[i] = std::log10(rhoH[i] * dom); + } + + // Calculate metallicity + + if (cloudy_table.grid_rank > 2) { + log_Z[i] = std::log10(metallicity[i]); + } + + // Calculate electron fraction + + if (cloudy_table.grid_rank > 3) { + e_frac[i] = 2. * de(i, idx_range.j, idx_range.k) / + (d(i, idx_range.j, idx_range.k) * (1. + fh[i])); + // Make sure electron fraction is never above 1 + // which can give bad cooling/heating values when + // extrapolating in the Cloudy data. + log_e_frac[i] = std::min(std::log10(e_frac[i]), 0.0); + + // Get extra electrons contributed by metals + + cl_e_frac[i] = + e_frac[i] * + (1. + (2. * my_chemistry->cloudy_electron_fraction_factor * + metallicity[i] * fh[i]) / + (1. + fh[i])); + } + + // Call interpolation functions to get heating/cooling + + // Interpolate over temperature. + if (cloudy_table.grid_rank == 1) { + log_cool[i] = grackle::impl::fortran_wrapper::interpolate_1d_g( + log10tem[i], cloudy_table.grid_dimension, + cloudy_table.grid_parameters[0], dclPar[0], cloudy_table.data_size, + cloudy_table.cooling_data); + edot_met[i] = -std::pow(10., log_cool[i]); + + // Ignore CMB term if T >> T_CMB + if ((my_chemistry->cmb_temperature_floor == 1) && + ((log10tem[i] - log10_tCMB) < 2.)) { + log_cool_cmb[i] = grackle::impl::fortran_wrapper::interpolate_1d_g( + log10_tCMB, cloudy_table.grid_dimension, + cloudy_table.grid_parameters[0], dclPar[0], + cloudy_table.data_size, cloudy_table.cooling_data); + edot_met[i] = edot_met[i] + std::pow(10., log_cool_cmb[i]); + } + + if (my_chemistry->UVbackground == 1) { + log_heat[i] = grackle::impl::fortran_wrapper::interpolate_1d_g( + log10tem[i], cloudy_table.grid_dimension, + cloudy_table.grid_parameters[0], dclPar[0], + cloudy_table.data_size, cloudy_table.cooling_data); + edot_met[i] = edot_met[i] + std::pow(10., log_heat[i]); + } + + // Interpolate over density and temperature. + } else if (cloudy_table.grid_rank == 2) { + log_cool[i] = grackle::impl::fortran_wrapper::interpolate_2d_g( + log_n_h[i], log10tem[i], cloudy_table.grid_dimension, + cloudy_table.grid_parameters[0], dclPar[0], + cloudy_table.grid_parameters[1], dclPar[1], cloudy_table.data_size, + cloudy_table.cooling_data); + edot_met[i] = -std::pow(10., log_cool[i]); + + // Ignore CMB term if T >> T_CMB + if ((my_chemistry->cmb_temperature_floor == 1) && + ((log10tem[i] - log10_tCMB) < 2.0f)) { + log_cool_cmb[i] = grackle::impl::fortran_wrapper::interpolate_2d_g( + log_n_h[i], log10_tCMB, cloudy_table.grid_dimension, + cloudy_table.grid_parameters[0], dclPar[0], + cloudy_table.grid_parameters[1], dclPar[1], + cloudy_table.data_size, cloudy_table.cooling_data); + edot_met[i] = edot_met[i] + std::pow(10., log_cool_cmb[i]); + } + + if (my_chemistry->UVbackground == 1) { + log_heat[i] = grackle::impl::fortran_wrapper::interpolate_2d_g( + log_n_h[i], log10tem[i], cloudy_table.grid_dimension, + cloudy_table.grid_parameters[0], dclPar[0], + cloudy_table.grid_parameters[1], dclPar[1], + cloudy_table.data_size, cloudy_table.heating_data); + edot_met[i] = edot_met[i] + std::pow(10., log_heat[i]); + } + + // Interpolate over density, metallicity, and temperature. + } else if (cloudy_table.grid_rank == 3) { + log_cool[i] = grackle::impl::fortran_wrapper::interpolate_3d_g( + log_n_h[i], log_Z[i], log10tem[i], cloudy_table.grid_dimension, + cloudy_table.grid_parameters[0], dclPar[0], + cloudy_table.grid_parameters[1], dclPar[1], + cloudy_table.grid_parameters[2], dclPar[2], cloudy_table.data_size, + cloudy_table.cooling_data); + edot_met[i] = -std::pow(10., log_cool[i]); + + // Ignore CMB term if T >> T_CMB + if ((my_chemistry->cmb_temperature_floor == 1) && + ((log10tem[i] - log10_tCMB) < 2.)) { + log_cool_cmb[i] = grackle::impl::fortran_wrapper::interpolate_3d_g( + log_n_h[i], log_Z[i], log10_tCMB, cloudy_table.grid_dimension, + cloudy_table.grid_parameters[0], dclPar[0], + cloudy_table.grid_parameters[1], dclPar[1], + cloudy_table.grid_parameters[2], dclPar[2], + cloudy_table.data_size, cloudy_table.cooling_data); + edot_met[i] = edot_met[i] + std::pow(10., log_cool_cmb[i]); + } + + if (my_chemistry->UVbackground == 1) { + log_heat[i] = grackle::impl::fortran_wrapper::interpolate_3d_g( + log_n_h[i], log_Z[i], log10tem[i], cloudy_table.grid_dimension, + cloudy_table.grid_parameters[0], dclPar[0], + cloudy_table.grid_parameters[1], dclPar[1], + cloudy_table.grid_parameters[2], dclPar[2], + cloudy_table.data_size, cloudy_table.heating_data); + edot_met[i] = edot_met[i] + std::pow(10., log_heat[i]); + } + + // Interpolate over density, metallicity, electron fraction, and + // temperature. + } else if (cloudy_table.grid_rank == 4) { + log_cool[i] = grackle::impl::fortran_wrapper::interpolate_4d_g( + log_n_h[i], log_Z[i], log_e_frac[i], log10tem[i], + cloudy_table.grid_dimension, cloudy_table.grid_parameters[0], + dclPar[0], cloudy_table.grid_parameters[1], dclPar[1], + cloudy_table.grid_parameters[2], dclPar[2], + cloudy_table.grid_parameters[3], dclPar[3], cloudy_table.data_size, + cloudy_table.cooling_data); + edot_met[i] = -std::pow(10., log_cool[i]); + + // Ignore CMB term if T >> T_CMB + if ((my_chemistry->cmb_temperature_floor == 1) && + ((log10tem[i] - log10_tCMB) < 2.)) { + log_cool_cmb[i] = grackle::impl::fortran_wrapper::interpolate_4d_g( + log_n_h[i], log_Z[i], log_e_frac[i], log10_tCMB, + cloudy_table.grid_dimension, cloudy_table.grid_parameters[0], + dclPar[0], cloudy_table.grid_parameters[1], dclPar[1], + cloudy_table.grid_parameters[2], dclPar[2], + cloudy_table.grid_parameters[3], dclPar[3], + cloudy_table.data_size, cloudy_table.cooling_data); + edot_met[i] = edot_met[i] + std::pow(10., log_cool_cmb[i]); + } + + if (my_chemistry->UVbackground == 1) { + log_heat[i] = grackle::impl::fortran_wrapper::interpolate_4d_g( + log_n_h[i], log_Z[i], log_e_frac[i], log10tem[i], + cloudy_table.grid_dimension, cloudy_table.grid_parameters[0], + dclPar[0], cloudy_table.grid_parameters[1], dclPar[1], + cloudy_table.grid_parameters[2], dclPar[2], + cloudy_table.grid_parameters[3], dclPar[3], + cloudy_table.data_size, cloudy_table.heating_data); + edot_met[i] = edot_met[i] + std::pow(10., log_heat[i]); + } + + // Interpolate over density, metallicity, electron fraction, redshift, + // and temperature. + } else { + log_cool[i] = grackle::impl::fortran_wrapper::interpolate_5d_g( + log_n_h[i], log_Z[i], log_e_frac[i], zr, log10tem[i], + cloudy_table.grid_dimension, cloudy_table.grid_parameters[0], + dclPar[0], cloudy_table.grid_parameters[1], dclPar[1], + cloudy_table.grid_parameters[2], dclPar[2], + cloudy_table.grid_parameters[3], dclPar[3], + cloudy_table.grid_parameters[4], dclPar[4], cloudy_table.data_size, + cloudy_table.cooling_data); + edot_met[i] = -std::pow(10., log_cool[i]); + + // Ignore CMB term if T >> T_CMB + if ((my_chemistry->cmb_temperature_floor == 1) && + ((log10tem[i] - log10_tCMB) < 2.)) { + log_cool_cmb[i] = grackle::impl::fortran_wrapper::interpolate_5d_g( + log_n_h[i], log_Z[i], log_e_frac[i], zr, log10_tCMB, + cloudy_table.grid_dimension, cloudy_table.grid_parameters[0], + dclPar[0], cloudy_table.grid_parameters[1], dclPar[1], + cloudy_table.grid_parameters[2], dclPar[2], + cloudy_table.grid_parameters[3], dclPar[3], + cloudy_table.grid_parameters[4], dclPar[4], + cloudy_table.data_size, cloudy_table.cooling_data); + edot_met[i] = edot_met[i] + std::pow(10., log_cool_cmb[i]); + } + + if (my_chemistry->UVbackground == 1) { + log_heat[i] = grackle::impl::fortran_wrapper::interpolate_5d_g( + log_n_h[i], log_Z[i], log_e_frac[i], zr, log10tem[i], + cloudy_table.grid_dimension, cloudy_table.grid_parameters[0], + dclPar[0], cloudy_table.grid_parameters[1], dclPar[1], + cloudy_table.grid_parameters[2], dclPar[2], + cloudy_table.grid_parameters[3], dclPar[3], + cloudy_table.grid_parameters[4], dclPar[4], + cloudy_table.data_size, cloudy_table.heating_data); + edot_met[i] = edot_met[i] + std::pow(10., log_heat[i]); + } + } + + if (cloudy_table.grid_rank > 3) { + edot_met[i] = edot_met[i] * cl_e_frac[i]; + } + + edot[i] = + edot[i] + (edot_met[i] * rhoH[i] * d(i, idx_range.j, idx_range.k)); + } + } + + return; +} diff --git a/src/clib/cool1d_cloudy_old_tables_g.hpp b/src/clib/cool1d_cloudy_old_tables_g.hpp new file mode 100644 index 000000000..be8e09e80 --- /dev/null +++ b/src/clib/cool1d_cloudy_old_tables_g.hpp @@ -0,0 +1,55 @@ +// See LICENSE file for license and copyright information + +/// @file cool1d_cloudy_old_tables_g-cpp.h +/// @brief Declares signature of cool1d_cloudy_old_tables_g + +// This file was initially generated automatically during conversion of the +// cool1d_cloudy_old_tables_g function from FORTRAN to C++ + +#ifndef COOL1D_CLOUDY_OLD_TABLES_G_HPP +#define COOL1D_CLOUDY_OLD_TABLES_G_HPP + +#include "grackle.h" // gr_float +#include "fortran_func_decls.h" // gr_mask_int +#include "index_helper.h" // IndexRange + +namespace grackle::impl { + +/// Solve cloudy cooling by interpolating from the data. This version uses +/// tables formatted for the original Cloudy cooling functionality in Enzo. +/// +/// Solve cloudy metal cooling for old cloudy metal cooling tables +/// +/// @param[in] rhoH 1D array to hold the computed Hydrogen mass density for the +/// @p idx_range +/// @param[in] metallicity 1D array to hold the computed metallicity for the @p +/// idx_range +/// @param[in] logtem Natural log of temperature values +/// @param[out] edot 1D array to hold the computed the time derivative of the +/// internal energy in the @p idx_range +/// @param[in] comp2 Constant factor to convert cooling rates to code units +/// @param[in] dom Unit conversion to proper number density in code units +/// @param[in] zr Current redshift +/// @param[in] itmask Iteration mask +/// @param[in] my_chemistry holds a number of configuration parameters. +/// @param[in] my_rates Holds assorted rate data and other internal +/// configuration info. +/// @param[in] cloudy_table Cloudy cooling table data +/// @param[in] density total density field +/// @param[in] e_density electron density field +/// @param[in] my_fields Specifies the field data. +/// @param[in] idx_range Specifies the current index-range +/// +/// @par History +/// written by: Britton Smith, 2009 +/// modified1: November, 2025 by Christopher Bignamini & Matthew Abruzzo; C++ +/// port +void cool1d_cloudy_old_tables_g( + const double* rhoH, double* metallicity, const double* logtem, double* edot, + double comp2, double dom, double zr, const gr_mask_type* itmask, + chemistry_data* my_chemistry, cloudy_data cloudy_table, gr_float* density, + gr_float* e_density, grackle_field_data* my_fields, IndexRange idx_range); + +} // namespace grackle::impl + +#endif /* COOL1D_CLOUDY_OLD_TABLES_G_HPP */ \ No newline at end of file diff --git a/src/clib/cool1d_multi_g.cpp b/src/clib/cool1d_multi_g.cpp index 02793cc90..11e5b72e6 100644 --- a/src/clib/cool1d_multi_g.cpp +++ b/src/clib/cool1d_multi_g.cpp @@ -17,6 +17,7 @@ #include #include +#include "cool1d_cloudy_old_tables_g.hpp" #include "cool1d_multi_g.hpp" #include "grackle.h" #include "fortran_func_decls.h" @@ -1674,24 +1675,10 @@ void grackle::impl::cool1d_multi_g( my_rates->cloudy_metal.heating_data, itmask_tab.data()); } else { - FORTRAN_NAME(cool1d_cloudy_old_tables_g)( - d.data(), de.data(), rhoH, metallicity, &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, - logTlininterp_buf.logtem, edot, &comp2, - &my_chemistry->primordial_chemistry, &dom, &zr, - &my_chemistry->cmb_temperature_floor, &my_chemistry->UVbackground, - &my_chemistry->cloudy_electron_fraction_factor, - &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, itmask_tab.data()); + grackle::impl::cool1d_cloudy_old_tables_g( + rhoH, metallicity, logTlininterp_buf.logtem, edot, comp2, dom, zr, + itmask_tab.data(), my_chemistry, my_rates->cloudy_metal, + my_fields->density, my_fields->e_density, my_fields, idx_range); } if (my_chemistry->metal_chemistry == 1) { diff --git a/src/clib/fortran_func_decls.h b/src/clib/fortran_func_decls.h index ff6fd5755..21ade65b7 100644 --- a/src/clib/fortran_func_decls.h +++ b/src/clib/fortran_func_decls.h @@ -139,16 +139,6 @@ void FORTRAN_NAME(cool1d_cloudy_g)( double* clHeating, gr_mask_type* itmask ); -void FORTRAN_NAME(cool1d_cloudy_old_tables_g)( - gr_float* d_data_ptr, gr_float* de_data_ptr, double* rhoH, - double* metallicity, int* in, int* jn, int* kn, int* is, int* ie, int* j, - int* k, double* logtem, double* edot, double* comp2, int* ispecies, - double* dom, double* zr, int* icmbTfloor, int* iClHeat, double* clEleFra, - long long* clGridRank, long long* clGridDim, double* clPar1, double* clPar2, - double* clPar3, double* clPar4, double* clPar5, long long* clDataSize, - double* clCooling, double* clHeating, gr_mask_type* itmask -); - void FORTRAN_NAME(gaussj_g)( int* n, double* a_data_ptr, double* b, int* ierr );