diff --git a/src/clib/CMakeLists.txt b/src/clib/CMakeLists.txt index e717de268..68800c8f9 100644 --- a/src/clib/CMakeLists.txt +++ b/src/clib/CMakeLists.txt @@ -109,6 +109,8 @@ add_library(Grackle_Grackle cool1d_multi_g.cpp cool1d_multi_g.hpp cool_multi_time_g.cpp cool_multi_time_g.h dust_props.hpp + dust/calc_all_tdust_gasgr_1d_g.cpp dust/calc_all_tdust_gasgr_1d_g.hpp + dust/calc_tdust_1d_g.cpp dust/calc_tdust_1d_g.hpp dust/grain_species_info.cpp dust/grain_species_info.hpp dust/lookup_dust_rates1d.hpp dust/multi_grain_species/calc_grain_size_increment_1d.hpp @@ -149,7 +151,6 @@ add_library(Grackle_Grackle visitor/printer.hpp # Fortran Source Files - calc_all_tdust_gasgr_1d_g.F calc_tdust_1d_g.F interpolators_g.F calc_grain_size_increment_1d.F diff --git a/src/clib/Make.config.objects b/src/clib/Make.config.objects index 7b6ece514..a7c587ee8 100644 --- a/src/clib/Make.config.objects +++ b/src/clib/Make.config.objects @@ -25,6 +25,8 @@ OBJS_CONFIG_LIB = \ cool1d_cloudy_old_tables_g.lo \ cool1d_multi_g.lo \ cool_multi_time_g.lo \ + dust/calc_all_tdust_gasgr_1d_g.lo \ + dust/calc_tdust_1d_g.lo \ dust/grain_species_info.lo \ dynamic_api.lo \ grackle_units.lo \ @@ -45,7 +47,6 @@ OBJS_CONFIG_LIB = \ solve_rate_cool.lo \ status_reporting.lo \ update_UVbackground_rates.lo \ - calc_all_tdust_gasgr_1d_g.lo \ calc_tdust_1d_g.lo \ calc_tdust_3d.lo \ calc_grain_size_increment_1d.lo \ diff --git a/src/clib/calc_all_tdust_gasgr_1d_g.F b/src/clib/calc_all_tdust_gasgr_1d_g.F deleted file mode 100644 index 7e5c10c75..000000000 --- a/src/clib/calc_all_tdust_gasgr_1d_g.F +++ /dev/null @@ -1,353 +0,0 @@ - subroutine calc_all_tdust_gasgr_1d_g(in, jn, kn, nratec, - & idustfield, is, ie, j, k, fgr, gamma_isrfa, - & trad, gasgra, indixe, tdef, tgas, tdust, - & metallicity, dust2gas, nh, gasgr_tdust, - & itmask_metal, - & idspecies, itdmulti, gr_N, gr_Size, gr_dT, - & gr_Td, tSiM, tFeM, tMg2SiO4, tMgSiO3, tFe3O4, - & tAC, tSiO2D, tMgO, tFeS, tAl2O3, treforg, - & tvolorg, tH2Oice, gasgr2a, gamma_isrf2a, - & coolunit, gasgr, myisrf, sgSiM, sgFeM, sgMg2SiO4, - & sgMgSiO3, sgFe3O4, sgAC, sgSiO2D, sgMgO, sgFeS, - & sgAl2O3, sgreforg, sgvolorg, sgH2Oice, sgtot, - & alSiM, alFeM, alMg2SiO4, alMgSiO3, alFe3O4, alAC, - & alSiO2D, alMgO, alFeS, alAl2O3, alreforg, - & alvolorg, alH2Oice, altot, kpSiM, kpFeM, - & kpMg2SiO4, kpMgSiO3, kpFe3O4, kpAC, kpSiO2D, - & kpMgO, kpFeS, kpAl2O3, kpreforg, kpvolorg, - & kpH2Oice, kptot, gasSiM, gasFeM, gasMg2SiO4, - & gasMgSiO3, gasFe3O4, gasAC, gasSiO2D, gasMgO, - & gasFeS, gasAl2O3, gasreforg, gasvolorg, - & gasH2Oice - & ) - -! PURPOSE: -! Calculate all dust temperature(s) and the gas to grain heat -! transfer rate(s) (the latter is commonly called gasgr) -! -! An argument could be made for directly using gasgr to compute -! contributions to edot within this routine, rather than returning -! the gasgr value(s). But that should be reconsidered in the future. -! -! We could significantly reduce the amount of buffer space allocated -! within the routine. But, we will hold off on doing that until after -! we have transcribed to C/C++ -! -! INPUTS: -! is,ie - start and end indicies of active region (zero-based!) -! -! PARAMETERS: -! - -!----------------------------------------------------------------------- - - implicit NONE -#include "grackle_fortran_types.def" -#include "phys_const.def" - -! Arguments - - integer in, jn, kn, is, ie, j, k, nratec, idustfield - - real*8 fgr - real*8 gasgra(nratec), gamma_isrfa - -! Parameters - - real*8, parameter :: mh = mass_h !DPC - -! Locals - - integer i - real*8 trad, coolunit - -! Slice locals - - integer*8 indixe(in) - real*8 tdef(in), tgas(in), tdust(in), nh(in), metallicity(in), - & dust2gas(in) - -! Cooling/heating slice locals - - real*8 gasgr(in), gasgr_tdust(in), myisrf(in) - integer idspecies, itdmulti -! opacity table - integer gr_N(2), gr_Size - real*8 gr_dT, gr_Td(gr_N(2)) -! grain growth - real*8 sgSiM(in), sgFeM(in), sgMg2SiO4(in) - & , sgMgSiO3(in), sgFe3O4(in), sgAC(in) - & , sgSiO2D(in), sgMgO(in), sgFeS(in) - & , sgAl2O3(in) - & , sgreforg(in), sgvolorg(in), sgH2Oice(in) - & , sgtot(in) - real*8 alSiM(gr_N(2),in), alFeM(gr_N(2),in) - & , alMg2SiO4(gr_N(2),in), alMgSiO3(gr_N(2),in) - & , alFe3O4(gr_N(2),in), alAC(gr_N(2),in) - & , alSiO2D(gr_N(2),in), alMgO(gr_N(2),in) - & , alFeS(gr_N(2),in), alAl2O3(gr_N(2),in) - & , alreforg(gr_N(2),in) - & , alvolorg(gr_N(2),in), alH2Oice(gr_N(2),in) - & , altot(gr_N(2),in) - real*8 kpSiM(in), kpFeM(in), kpMg2SiO4(in) - & , kpMgSiO3(in), kpFe3O4(in), kpAC(in) - & , kpSiO2D(in), kpMgO(in), kpFeS(in) - & , kpAl2O3(in) - & , kpreforg(in), kpvolorg(in), kpH2Oice(in) - & , kptot(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) - real*8 gasgr2a(nratec), gamma_isrf2a - real*8 gasSiM(in), gasFeM(in), gasMg2SiO4(in) - & , gasMgSiO3(in), gasFe3O4(in), gasAC(in) - & , gasSiO2D(in), gasMgO(in), gasFeS(in) - & , gasAl2O3(in) - & , gasreforg(in), gasvolorg(in), gasH2Oice(in) - real*8 gasgr_tSiM(in), gasgr_tFeM(in), gasgr_tMg2SiO4(in) - & , gasgr_tMgSiO3(in), gasgr_tFe3O4(in), gasgr_tAC(in) - & , gasgr_tSiO2D(in), gasgr_tMgO(in), gasgr_tFeS(in) - & , gasgr_tAl2O3(in) - & , gasgr_treforg(in), gasgr_tvolorg(in), gasgr_tH2Oice(in) - real*8 mygisrf(in), fv2k, fac - real*8 gisrfSiM(in), gisrfFeM(in), gisrfMg2SiO4(in) - & , gisrfMgSiO3(in), gisrfFe3O4(in), gisrfAC(in) - & , gisrfSiO2D(in), gisrfMgO(in), gisrfFeS(in) - & , gisrfAl2O3(in) - & , gisrfreforg(in), gisrfvolorg(in), gisrfH2Oice(in) -! Iteration mask - - MASK_TYPE itmask_metal(in) - - ! Calculate heating from interstellar radiation field - ! -> this is ONLY used when `itmask_metal .eq. MASK_TRUE` - - do i = is+1, ie+1 - if ( itmask_metal(i) .ne. MASK_FALSE ) then - - if (idspecies .eq. 0 ) then - if (idustfield .gt. 0) then - mygisrf(i) = gamma_isrfa - & * fgr / dust2gas(i) * metallicity(i) - !! correct with the depletion or enhancement of condensation rate. - else - mygisrf(i) = gamma_isrfa - endif - - else ! idspecies - - if (itdmulti .eq. 0) then - - mygisrf(i) = gamma_isrf2a * sgtot(i) - !! in UV, absorption coefficient Q ~ 1 (Goldsmith 2001) - !! so we use the geometrical cross-section of grains [cgs] - - else - - if (idspecies .gt. 0) then - gisrfMgSiO3 (i) = gamma_isrf2a * sgMgSiO3 (i) -!! write(*,*) 'sil', d(i,j,k), gamma_isrf2a, sgMgSiO3(i) - gisrfAC (i) = gamma_isrf2a * sgAC (i) -!! write(*,*) 'car', d(i,j,k), gamma_isrf2a, sgMgSiO3(i) - endif - if (idspecies .gt. 1) then - gisrfSiM (i) = gamma_isrf2a * sgSiM (i) - gisrfFeM (i) = gamma_isrf2a * sgFeM (i) - gisrfMg2SiO4 (i) = gamma_isrf2a * sgMg2SiO4 (i) - gisrfFe3O4 (i) = gamma_isrf2a * sgFe3O4 (i) - gisrfSiO2D (i) = gamma_isrf2a * sgSiO2D (i) - gisrfMgO (i) = gamma_isrf2a * sgMgO (i) - gisrfFeS (i) = gamma_isrf2a * sgFeS (i) - gisrfAl2O3 (i) = gamma_isrf2a * sgAl2O3 (i) - endif - if (idspecies .gt. 2) then - gisrfreforg (i) = gamma_isrf2a * sgreforg (i) - gisrfvolorg (i) = gamma_isrf2a * sgvolorg (i) - gisrfH2Oice (i) = gamma_isrf2a * sgH2Oice (i) - endif - - endif - - endif ! idspecies - - endif - enddo - -! --- Gas to grain heat transfer --- - -! Look up gas/grain heat transfer rates - - do i = is+1, ie+1 - if ( itmask_metal(i) .ne. MASK_FALSE ) then - - if(idspecies .eq. 0) then - - gasgr(i) = gasgra(indixe(i)) + tdef(i) - & *(gasgra(indixe(i)+1) -gasgra(indixe(i))) - -!! gasgr_tdust(i) = fgr * gasgr(i) * coolunit / mh - gasgr_tdust(i) = (dust2gas(i) / metallicity(i)) - & * gasgr(i) * coolunit / mh - !! apply to (idustfield .eq. 1) GC20200701 - - else ! idspecies - - fv2k = gasgr2a(indixe(i)) + tdef(i) - & *(gasgr2a(indixe(i)+1) -gasgr2a(indixe(i))) - - fac = coolunit / mh - - if ( itdmulti .eq. 0 ) then - - gasgr(i) = fv2k * sgtot(i) - - gasgr_tdust(i) = gasgr(i) * fac - - else - - if (idspecies .gt. 0) then - gasMgSiO3 (i) = fv2k * sgMgSiO3 (i) - gasAC (i) = fv2k * sgAC (i) - endif - if (idspecies .gt. 1) then - gasSiM (i) = fv2k * sgSiM (i) - gasFeM (i) = fv2k * sgFeM (i) - gasMg2SiO4 (i) = fv2k * sgMg2SiO4 (i) - gasFe3O4 (i) = fv2k * sgFe3O4 (i) - gasSiO2D (i) = fv2k * sgSiO2D (i) - gasMgO (i) = fv2k * sgMgO (i) - gasFeS (i) = fv2k * sgFeS (i) - gasAl2O3 (i) = fv2k * sgAl2O3 (i) - endif - if (idspecies .gt. 2) then - gasreforg (i) = fv2k * sgreforg (i) - gasvolorg (i) = fv2k * sgvolorg (i) - gasH2Oice (i) = fv2k * sgH2Oice (i) - endif - - if (idspecies .gt. 0) then - gasgr_tMgSiO3 (i) = gasMgSiO3 (i) * fac - gasgr_tAC (i) = gasAC (i) * fac - endif - if (idspecies .gt. 1) then - gasgr_tSiM (i) = gasSiM (i) * fac - gasgr_tFeM (i) = gasFeM (i) * fac - gasgr_tMg2SiO4 (i) = gasMg2SiO4 (i) * fac - gasgr_tFe3O4 (i) = gasFe3O4 (i) * fac - gasgr_tSiO2D (i) = gasSiO2D (i) * fac - gasgr_tMgO (i) = gasMgO (i) * fac - gasgr_tFeS (i) = gasFeS (i) * fac - gasgr_tAl2O3 (i) = gasAl2O3 (i) * fac - endif - if (idspecies .gt. 2) then - gasgr_treforg (i) = gasreforg (i) * fac - gasgr_tvolorg (i) = gasvolorg (i) * fac - gasgr_tH2Oice (i) = gasH2Oice (i) * fac - endif - - endif - - endif ! idspecies - - endif - enddo - -! Compute dust temperature - - if (itdmulti .eq. 0) then - - call calc_tdust_1d_g(tdust, tgas, nh, gasgr_tdust, - & mygisrf, myisrf, itmask_metal, trad, in, is, ie, j, k - & , gr_N, gr_Size, gr_dT, gr_Td, altot, kptot, idspecies) - - else - - if (idspecies .gt. 0) then - call calc_tdust_1d_g(tMgSiO3 , tgas, nh, gasgr_tMgSiO3 , - & gisrfMgSiO3 , myisrf, itmask_metal, trad, - & in, is, ie, j, k - & , gr_N, gr_Size, gr_dT, gr_Td, alMgSiO3 , kpMgSiO3 - & , idspecies) - - call calc_tdust_1d_g(tAC , tgas, nh, gasgr_tAC , - & gisrfAC , myisrf, itmask_metal, trad, - & in, is, ie, j, k - & , gr_N, gr_Size, gr_dT, gr_Td, alAC , kpAC - & , idspecies) - endif - - if (idspecies .gt. 1) then - call calc_tdust_1d_g(tSiM , tgas, nh, gasgr_tSiM , - & gisrfSiM , myisrf, itmask_metal, trad, - & in, is, ie, j, k - & , gr_N, gr_Size, gr_dT, gr_Td, alSiM , kpSiM - & , idspecies) - - call calc_tdust_1d_g(tFeM , tgas, nh, gasgr_tFeM , - & gisrfFeM , myisrf, itmask_metal, trad, - & in, is, ie, j, k - & , gr_N, gr_Size, gr_dT, gr_Td, alFeM , kpFeM - & , idspecies) - - call calc_tdust_1d_g(tMg2SiO4 , tgas, nh, gasgr_tMg2SiO4 , - & gisrfMg2SiO4 , myisrf, itmask_metal, trad, - & in, is, ie, j, k - & , gr_N, gr_Size, gr_dT, gr_Td, alMg2SiO4 , kpMg2SiO4 - & , idspecies) - - call calc_tdust_1d_g(tFe3O4 , tgas, nh, gasgr_tFe3O4 , - & gisrfFe3O4 , myisrf, itmask_metal, trad, - & in, is, ie, j, k - & , gr_N, gr_Size, gr_dT, gr_Td, alFe3O4 , kpFe3O4 - & , idspecies) - - call calc_tdust_1d_g(tSiO2D , tgas, nh, gasgr_tSiO2D , - & gisrfSiO2D , myisrf, itmask_metal, trad, - & in, is, ie, j, k - & , gr_N, gr_Size, gr_dT, gr_Td, alSiO2D , kpSiO2D - & , idspecies) - - call calc_tdust_1d_g(tMgO , tgas, nh, gasgr_tMgO , - & gisrfMgO , myisrf, itmask_metal, trad, - & in, is, ie, j, k - & , gr_N, gr_Size, gr_dT, gr_Td, alMgO , kpMgO - & , idspecies) - - call calc_tdust_1d_g(tFeS , tgas, nh, gasgr_tFeS , - & gisrfFeS , myisrf, itmask_metal, trad, - & in, is, ie, j, k - & , gr_N, gr_Size, gr_dT, gr_Td, alFeS , kpFeS - & , idspecies) - - call calc_tdust_1d_g(tAl2O3 , tgas, nh, gasgr_tAl2O3 , - & gisrfAl2O3 , myisrf, itmask_metal, trad, - & in, is, ie, j, k - & , gr_N, gr_Size, gr_dT, gr_Td, alAl2O3 , kpAl2O3 - & , idspecies) - endif - - if (idspecies .gt. 2) then - call calc_tdust_1d_g(treforg , tgas, nh, gasgr_treforg , - & gisrfreforg , myisrf, itmask_metal, trad, - & in, is, ie, j, k - & , gr_N, gr_Size, gr_dT, gr_Td, alreforg , kpreforg - & , idspecies) - - call calc_tdust_1d_g(tvolorg , tgas, nh, gasgr_tvolorg , - & gisrfvolorg , myisrf, itmask_metal, trad, - & in, is, ie, j, k - & , gr_N, gr_Size, gr_dT, gr_Td, alvolorg , kpvolorg - & , idspecies) - - call calc_tdust_1d_g(tH2Oice , tgas, nh, gasgr_tH2Oice , - & gisrfH2Oice , myisrf, itmask_metal, trad, - & in, is, ie, j, k - & , gr_N, gr_Size, gr_dT, gr_Td, alH2Oice , kpH2Oice - & , idspecies) - endif - - endif - end diff --git a/src/clib/calc_tdust_1d_g.F b/src/clib/calc_tdust_1d_g.F index bbeb792ec..6d3b19202 100644 --- a/src/clib/calc_tdust_1d_g.F +++ b/src/clib/calc_tdust_1d_g.F @@ -1,348 +1,5 @@ #include "phys_const.def" -!======================================================================= -!//////////////////// SUBROUTINE CALC_TDUST_1D_G \\\\\\\\\\\\\\\\\\\\\ - - subroutine calc_tdust_1d_g( - & tdust, tgas, nh, gasgr, gamma_isrfa, isrf, itmask, - & trad, in, is, ie, j, k - & , gr_N, gr_Size, gr_dT, gr_Td, alsp, kgr, idspecies - & ) - -! CALCULATE EQUILIBRIUM DUST TEMPERATURE -! -! written by: Britton Smith -! date: February, 2011 -! modified1: -! -! PURPOSE: -! Calculate dust temperature. -! -! TODO: -! this docstring, and the docstrings of all helper functions should -! EXPLICITLY document how the meaning of arguments change based on -! whether we are using the single-field dust model or the -! multi-species dust model. The different meaning of each variable -! gets VERY confusing -! -! INPUTS: -! in - dimension of 1D slice -! -! tdust - dust temperature -! -! tgas - gas temperature -! nh - H number density -! gasgr - gas/grain heat transfer rate -! gamma_isrfa - heating from interstellar radiation field -! isrf - interstellar radiation field in Habing units -! -! trad - CMB temperature -! -! is,ie - start and end indices of active region (zero based) -! j,k - indices of 1D slice -! -! itmask - iteration mask -! -! PARAMETERS: -! -!----------------------------------------------------------------------- - - implicit NONE -#include "grackle_fortran_types.def" - -! Arguments - - integer in, is, ie, j, k - - real*8 tdust(in), tgas(in), nh(in), gasgr(in), isrf(in) - real*8 gamma_isrfa(in), trad - -! opacity table of a grain species - integer gr_N(2), gr_Size - real*8 gr_dT, gr_Td(gr_N(2)) - R_PREC alsp(gr_N(2),in) - R_PREC logalsp(gr_N(2),in) - integer Td_N(1), Td_Size - real*8 logkgr - -! Iteration mask - - MASK_TYPE itmask(in) - -! Parameters - - integer idspecies - ! grain sublimation temperature - real*8, parameter :: t_subl = 1.5e3_DKIND - real*8, parameter :: radf = 4._DKIND * sigma_sb - -! grain opacity from Omukai (2000, equation 17) normalized by -! the local dust-to-gas ratio, which in this work is 0.934e-2. - real*8, parameter :: kgr1 = 4.0e-4_DKIND / 0.00934d0 - - real*8 gamma_isrf(in) - real*8, parameter :: tol = 1.e-5_DKIND - real*8, parameter :: bi_tol = 1.e-3_DKIND - real*8, parameter :: minpert = 1.e-10_DKIND - integer, parameter :: itmax = 50 - integer, parameter :: bi_itmax = 30 - -! Locals - - integer i, iter, c_done, c_total, nm_done - - real*8 pert_i, trad4 - -! Slice Locals - - real*8 kgr(in), kgrplus(in), sol(in), solplus(in), - & slope(in), tdplus(in), tdustnow(in), tdustold(in), - & pert(in), - & bi_t_mid(in), bi_t_high(in) - MASK_TYPE nm_itmask(in), bi_itmask(in) - -!\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\///////////////////////////////// -!======================================================================= -!! write(*,*) 'aa', gasgr, gamma_isrfa, isrf - - pert_i = 1.e-3_DKIND - - trad = max(1._DKIND, trad) - trad4 = trad**4 - -! \sum rho_SN kappa_SN / \sum rho_SN ndust_SN - - Td_N(1) = gr_N(2) - Td_Size = gr_N(2) - do i = is+1, ie+1 - if ( itmask(i) .ne. MASK_FALSE ) then - logalsp(:,i) = log10(alsp(:,i)) - endif - enddo - -! Set total cells for calculation - - c_done = 0 - nm_done = 0 - c_total = ie - is + 1 - -! Set local iteration mask and initial guess - - do i = is+1, ie+1 - if ( itmask(i) .ne. MASK_FALSE ) then - gamma_isrf(i) = isrf(i) * gamma_isrfa(i) - endif - enddo - - do i = is+1, ie+1 - nm_itmask(i) = itmask(i) - bi_itmask(i) = itmask(i) - if ( nm_itmask(i) .ne. MASK_FALSE ) then - - if (trad .ge. tgas(i)) then - tdustnow(i) = trad - nm_itmask(i) = MASK_FALSE - bi_itmask(i) = MASK_FALSE - c_done = c_done + 1 - nm_done = nm_done + 1 - else if (tgas(i) .gt. t_subl) then -! Use bisection if T_gas > grain sublimation temperature. - nm_itmask(i) = MASK_FALSE - nm_done = nm_done + 1 - else - tdustnow(i) = max(trad, - & (gamma_isrf(i) / radf / kgr1)**0.17_DKIND) - pert(i) = pert_i - endif - - else - c_done = c_done + 1 - nm_done = nm_done + 1 - endif - enddo - -! Iterate to convergence with Newton's method - - do iter = 1, itmax - -! Loop over slice - - do i = is+1, ie+1 - if ( nm_itmask(i) .ne. MASK_FALSE ) then - - tdplus(i) = max(1.e-3_DKIND, ((1._DKIND + pert(i)) - & * tdustnow(i))) - - endif - enddo - -! Calculate grain opacities - - call calc_kappa_gr_g(tdustnow, kgr, nm_itmask, - & in, is, ie, t_subl - & , Td_N, Td_Size, gr_dT, gr_Td, logalsp, idspecies) - - call calc_kappa_gr_g(tdplus, kgrplus, nm_itmask, - & in, is, ie, t_subl - & , Td_N, Td_Size, gr_dT, gr_Td, logalsp, idspecies) - -! Calculate heating/cooling balance - - call calc_gr_balance_g(tdustnow, tgas, kgr, trad4, gasgr, - & gamma_isrf, nh, nm_itmask, sol, in, is, ie) - - call calc_gr_balance_g(tdplus, tgas, kgrplus, trad4, gasgr, - & gamma_isrf, nh, nm_itmask, solplus, in, is, ie) - - do i = is+1, ie+1 - if ( nm_itmask(i) .ne. MASK_FALSE ) then - -! Use Newton's method to solve for Tdust - - slope(i) = (solplus(i) - sol(i)) / - & (pert(i) * tdustnow(i)) - - tdustold(i) = tdustnow(i) -! tdustnow(i) = tdustnow(i) - (sol(i) / slope(i)) - tdustnow(i) = min(tdustnow(i) - (sol(i) / slope(i)) - & , 3e3_DKIND) - - pert(i) = max(min(pert(i), - & (0.5_DKIND * abs(tdustnow(i) - tdustold(i)) / - & tdustnow(i))), minpert) - -! If negative solution calculated, give up and wait for bisection step. - if (tdustnow(i) .lt. trad) then - nm_itmask(i) = MASK_FALSE - nm_done = nm_done + 1 -! Check for convergence of solution - else if (abs(sol(i) / solplus(i)) .lt. tol) then - nm_itmask(i) = MASK_FALSE - c_done = c_done + 1 - bi_itmask(i) = MASK_FALSE - nm_done = nm_done + 1 - endif - -! if ( nm_itmask(i) ) - endif - -! End loop over slice - enddo - -! Check for all cells converged - if (c_done .ge. c_total) go to 666 - -! Check for all cells done with Newton method -! This includes attempts where a negative solution was found - if (nm_done .ge. c_total) go to 555 - -! End iteration loop for Newton's method - enddo - - 555 continue - -! If iteration count exceeded, try once more with bisection - if (c_done .lt. c_total) then - do i = is+1, ie+1 - if ( bi_itmask(i) .ne. MASK_FALSE ) then - tdustnow(i) = trad -! bi_t_high(i) = tgas(i) - bi_t_high(i) = 3e3_DKIND - endif - enddo - - do iter = 1, bi_itmax - - do i = is+1, ie+1 - if ( bi_itmask(i) .ne. MASK_FALSE ) then - - bi_t_mid(i) = 0.5_DKIND * (tdustnow(i) + bi_t_high(i)) - if (iter .eq. 1) then - bi_t_mid(i) = min(bi_t_mid(i), t_subl) - endif - - endif - enddo - - call calc_kappa_gr_g(bi_t_mid, kgr, bi_itmask, - & in, is, ie, t_subl - & , Td_N, Td_Size, gr_dT, gr_Td, logalsp, idspecies) - - call calc_gr_balance_g(bi_t_mid, tgas, kgr, trad4, gasgr, - & gamma_isrf, nh, bi_itmask, sol, in, is, ie) - - do i = is+1, ie+1 - if ( bi_itmask(i) .ne. MASK_FALSE ) then - - if (sol(i) .gt. 0._DKIND) then - tdustnow(i) = bi_t_mid(i) - else - bi_t_high(i) = bi_t_mid(i) - endif - - if ((abs(bi_t_high(i) - tdustnow(i)) / tdustnow(i)) - & .le. bi_tol) then - bi_itmask(i) = MASK_FALSE - c_done = c_done + 1 - endif - -! Check for all cells converged - if (c_done .ge. c_total) go to 666 - -! if ( bi_itmask(i) ) - endif - -! End loop over slice - enddo - -! End iteration loop for bisection - enddo - -! If iteration count exceeded with bisection, end of the line. - if (iter .gt. itmax) then -#ifdef _OPENMP -!$omp critical -#endif - write(6,*) 'CALC_TDUST_1D_G failed using both methods for ', - & (c_total - c_done), 'cells.' -#ifdef _OPENMP -!$omp end critical -#endif -c ERROR_MESSAGE - endif - -! if (iter .gt. itmax) then - endif - - 666 continue - -! Copy values back to thrown slice - do i = is+1, ie+1 - if ( itmask(i) .ne. MASK_FALSE ) then - -! Check for bad solutions - if (tdustnow(i) .lt. 0._DKIND) then -#ifdef _OPENMP -!$omp critical -#endif - write(6, *) 'CALC_TDUST_1D_G Newton method - ', - & 'T_dust < 0: i = ', i, 'j = ', j, - & 'k = ', k, 'nh = ', nh(i), - & 't_gas = ', tgas(i), 't_rad = ', trad, - & 't_dust = ', tdustnow(i) -#ifdef _OPENMP -!$omp end critical -#endif -c ERROR_MESSAGE - endif - - tdust(i) = tdustnow(i) - endif - enddo - - return - end - !======================================================================= !//////////////////// SUBROUTINE CALC_KAPPA_GR_G \\\\\\\\\\\\\\\\\\\\\ diff --git a/src/clib/calc_tdust_3d.cpp b/src/clib/calc_tdust_3d.cpp index 2035b424d..24670ae9d 100644 --- a/src/clib/calc_tdust_3d.cpp +++ b/src/clib/calc_tdust_3d.cpp @@ -16,6 +16,7 @@ #include #include +#include "dust/calc_all_tdust_gasgr_1d_g.hpp" #include "calc_tdust_3d.h" #include "dust_props.hpp" #include "dust/multi_grain_species/calc_grain_size_increment_1d.hpp" @@ -37,8 +38,6 @@ void calc_tdust_3d_g( grackle_field_data* my_fields, InternalGrUnits internalu ) { - // shorten `grackle::impl::fortran_wrapper` to `f_wrap` within this function - namespace f_wrap = ::grackle::impl::fortran_wrapper; const double mh_local_var = mh_grflt; @@ -255,13 +254,13 @@ void calc_tdust_3d_g( } // Compute dust temperature(s) in the index-range - f_wrap::calc_all_tdust_gasgr_1d_g( - trad, tgas.data(), tdust.data(), metallicity.data(), dust2gas.data(), - nh.data(), gasgr_tdust.data(), itmask_metal.data(), internalu.coolunit, - gasgr.data(), myisrf.data(), kappa_tot.data(), my_chemistry, my_rates, - my_fields, idx_range, grain_temperatures, gas_grainsp_heatrate, - grain_kappa, logTlininterp_buf, internal_dust_prop_buf - ); + grackle::impl::calc_all_tdust_gasgr_1d_g( + trad, tgas.data(), tdust.data(), metallicity.data(), + dust2gas.data(), nh.data(), gasgr_tdust.data(), itmask_metal.data(), + internalu.coolunit, gasgr.data(), myisrf.data(), kappa_tot.data(), + my_chemistry, my_rates, my_fields, idx_range, grain_temperatures, + gas_grainsp_heatrate, logTlininterp_buf, internal_dust_prop_buf, + grain_kappa); // Copy slice values back to grid diff --git a/src/clib/cool1d_multi_g.cpp b/src/clib/cool1d_multi_g.cpp index 02cc901c5..3d2a4aaab 100644 --- a/src/clib/cool1d_multi_g.cpp +++ b/src/clib/cool1d_multi_g.cpp @@ -25,6 +25,7 @@ #include "fortran_func_decls.h" #include "fortran_func_wrappers.hpp" #include "dust_props.hpp" +#include "dust/calc_all_tdust_gasgr_1d_g.hpp" #include "dust/multi_grain_species/calc_grain_size_increment_1d.hpp" #include "inject_model/grain_metal_inject_pathways.hpp" #include "internal_types.hpp" @@ -1155,12 +1156,12 @@ void grackle::impl::cool1d_multi_g( // compute dust temperature and cooling due to dust if (anydust != MASK_FALSE) { // TODO: trad -> comp2 - grackle::impl::fortran_wrapper::calc_all_tdust_gasgr_1d_g( + grackle::impl::calc_all_tdust_gasgr_1d_g( comp2, tgas, tdust, metallicity, dust2gas, cool1dmulti_buf.mynh, cool1dmulti_buf.gasgr_tdust, itmask_metal, coolunit, gasgr.data(), myisrf.data(), kappa_tot.data(), my_chemistry, my_rates, my_fields, - idx_range, grain_temperatures, gas_grainsp_heatrate, grain_kappa, - logTlininterp_buf, internal_dust_prop_buf); + idx_range, grain_temperatures, gas_grainsp_heatrate, logTlininterp_buf, + internal_dust_prop_buf, grain_kappa); } // Calculate dust cooling rate diff --git a/src/clib/dust/calc_all_tdust_gasgr_1d_g.cpp b/src/clib/dust/calc_all_tdust_gasgr_1d_g.cpp new file mode 100644 index 000000000..b66c26e16 --- /dev/null +++ b/src/clib/dust/calc_all_tdust_gasgr_1d_g.cpp @@ -0,0 +1,417 @@ +//===----------------------------------------------------------------------===// +// +// See the LICENSE file for license and copyright information +// SPDX-License-Identifier: NCSA AND BSD-3-Clause +// +//===----------------------------------------------------------------------===// +/// +/// @file +/// Implements the calc_all_tdust_gasgr_1d_g function +/// +//===----------------------------------------------------------------------===// + +// This file was initially generated automatically during conversion of the +// calc_all_tdust_gasgr_1d_g function from FORTRAN to C++ + +#include +#include + +#include "grackle.h" +#include "utils-cpp.hpp" + +#include "calc_tdust_1d_g.hpp" +#include "calc_all_tdust_gasgr_1d_g.hpp" +#include "inject_model/grain_metal_inject_pathways.hpp" +#include "opaque_storage.hpp" + +void grackle::impl::calc_all_tdust_gasgr_1d_g( + double trad, double* tgas, double* tdust, const double* metallicity, + const double* dust2gas, double* nh, double* gasgr_tdust, + gr_mask_type* itmask_metal, double coolunit, double* gasgr, double* myisrf, + double* kptot, chemistry_data* my_chemistry, + chemistry_data_storage* my_rates, grackle_field_data* my_fields, + IndexRange idx_range, + grackle::impl::GrainSpeciesCollection grain_temperatures, + grackle::impl::GrainSpeciesCollection gas_grainsp_heatrate, + grackle::impl::LogTLinInterpScratchBuf logTlininterp_buf, + grackle::impl::InternalDustPropBuf internal_dust_prop_buf, + grackle::impl::GrainSpeciesCollection grain_kappa) { + const double mh_local_var = mh_grflt; + int i; + + // Cooling/heating slice locals + + std::vector gasgr_tSiM(my_fields->grid_dimension[0]); + std::vector gasgr_tFeM(my_fields->grid_dimension[0]); + std::vector gasgr_tMg2SiO4(my_fields->grid_dimension[0]); + std::vector gasgr_tMgSiO3(my_fields->grid_dimension[0]); + std::vector gasgr_tFe3O4(my_fields->grid_dimension[0]); + std::vector gasgr_tAC(my_fields->grid_dimension[0]); + std::vector gasgr_tSiO2D(my_fields->grid_dimension[0]); + std::vector gasgr_tMgO(my_fields->grid_dimension[0]); + std::vector gasgr_tFeS(my_fields->grid_dimension[0]); + std::vector gasgr_tAl2O3(my_fields->grid_dimension[0]); + std::vector gasgr_treforg(my_fields->grid_dimension[0]); + std::vector gasgr_tvolorg(my_fields->grid_dimension[0]); + std::vector gasgr_tH2Oice(my_fields->grid_dimension[0]); + double fv2k, fac; + std::vector mygisrf(my_fields->grid_dimension[0]); + std::vector gisrfSiM(my_fields->grid_dimension[0]); + std::vector gisrfFeM(my_fields->grid_dimension[0]); + std::vector gisrfMg2SiO4(my_fields->grid_dimension[0]); + std::vector gisrfMgSiO3(my_fields->grid_dimension[0]); + std::vector gisrfFe3O4(my_fields->grid_dimension[0]); + std::vector gisrfAC(my_fields->grid_dimension[0]); + std::vector gisrfSiO2D(my_fields->grid_dimension[0]); + std::vector gisrfMgO(my_fields->grid_dimension[0]); + std::vector gisrfFeS(my_fields->grid_dimension[0]); + std::vector gisrfAl2O3(my_fields->grid_dimension[0]); + std::vector gisrfreforg(my_fields->grid_dimension[0]); + std::vector gisrfvolorg(my_fields->grid_dimension[0]); + std::vector gisrfH2Oice(my_fields->grid_dimension[0]); + + grackle::impl::GrainMetalInjectPathways* inject_pathway_props = + my_rates->opaque_storage->inject_pathway_props; + + double dlog10Tdust = 0.0; + double* log10Tdust_vals = nullptr; + + // NOTE: gr_N is a historical names + // -> it is pretty uninformative and should be changed! + int gr_N[2] = {0, 0}; + if (inject_pathway_props != nullptr) { + dlog10Tdust = + inject_pathway_props->log10Tdust_interp_props.parameter_spacing[0]; + log10Tdust_vals = + inject_pathway_props->log10Tdust_interp_props.parameters[0]; + + gr_N[0] = inject_pathway_props->n_opac_poly_coef; + gr_N[1] = static_cast( + inject_pathway_props->log10Tdust_interp_props.dimension[0]); + }; + + // Calculate heating from interstellar radiation field + // -> this is ONLY used when `itmask_metal .eq. MASK_TRUE` + + for (i = idx_range.i_start; i <= idx_range.i_end; i++) { + if (itmask_metal[i] != MASK_FALSE) { + if (my_chemistry->dust_species == 0) { + if (my_chemistry->use_dust_density_field > 0) { + mygisrf[i] = my_rates->gamma_isrf * + my_chemistry->local_dust_to_gas_ratio / dust2gas[i] * + metallicity[i]; + // ! correct with the depletion or enhancement of condensation rate. + } else { + mygisrf[i] = my_rates->gamma_isrf; + } + + } else { + if (my_chemistry->use_multiple_dust_temperatures == 0) { + mygisrf[i] = my_rates->gamma_isrf2 * + internal_dust_prop_buf.sigma_per_gas_mass_tot[i]; + // ! in UV, absorption coefficient Q ~ 1 (Goldsmith 2001) + // ! so we use the geometrical cross-section of grains [cgs] + + } else { + if (my_chemistry->dust_species > 0) { + gisrfMgSiO3[i] = my_rates->gamma_isrf2 * + internal_dust_prop_buf.grain_sigma_per_gas_mass + .data[OnlyGrainSpLUT::MgSiO3_dust][i]; + gisrfAC[i] = my_rates->gamma_isrf2 * + internal_dust_prop_buf.grain_sigma_per_gas_mass + .data[OnlyGrainSpLUT::AC_dust][i]; + } + if (my_chemistry->dust_species > 1) { + gisrfSiM[i] = my_rates->gamma_isrf2 * + internal_dust_prop_buf.grain_sigma_per_gas_mass + .data[OnlyGrainSpLUT::SiM_dust][i]; + gisrfFeM[i] = my_rates->gamma_isrf2 * + internal_dust_prop_buf.grain_sigma_per_gas_mass + .data[OnlyGrainSpLUT::FeM_dust][i]; + gisrfMg2SiO4[i] = my_rates->gamma_isrf2 * + internal_dust_prop_buf.grain_sigma_per_gas_mass + .data[OnlyGrainSpLUT::Mg2SiO4_dust][i]; + gisrfFe3O4[i] = my_rates->gamma_isrf2 * + internal_dust_prop_buf.grain_sigma_per_gas_mass + .data[OnlyGrainSpLUT::Fe3O4_dust][i]; + gisrfSiO2D[i] = my_rates->gamma_isrf2 * + internal_dust_prop_buf.grain_sigma_per_gas_mass + .data[OnlyGrainSpLUT::SiO2_dust][i]; + gisrfMgO[i] = my_rates->gamma_isrf2 * + internal_dust_prop_buf.grain_sigma_per_gas_mass + .data[OnlyGrainSpLUT::MgO_dust][i]; + gisrfFeS[i] = my_rates->gamma_isrf2 * + internal_dust_prop_buf.grain_sigma_per_gas_mass + .data[OnlyGrainSpLUT::FeS_dust][i]; + gisrfAl2O3[i] = my_rates->gamma_isrf2 * + internal_dust_prop_buf.grain_sigma_per_gas_mass + .data[OnlyGrainSpLUT::Al2O3_dust][i]; + } + if (my_chemistry->dust_species > 2) { + gisrfreforg[i] = my_rates->gamma_isrf2 * + internal_dust_prop_buf.grain_sigma_per_gas_mass + .data[OnlyGrainSpLUT::ref_org_dust][i]; + gisrfvolorg[i] = my_rates->gamma_isrf2 * + internal_dust_prop_buf.grain_sigma_per_gas_mass + .data[OnlyGrainSpLUT::vol_org_dust][i]; + gisrfH2Oice[i] = my_rates->gamma_isrf2 * + internal_dust_prop_buf.grain_sigma_per_gas_mass + .data[OnlyGrainSpLUT::H2O_ice_dust][i]; + } + } + } + } + } + + // --- Gas to grain heat transfer --- + + // Look up gas/grain heat transfer rates + + for (i = idx_range.i_start; i <= idx_range.i_end; i++) { + if (itmask_metal[i] != MASK_FALSE) { + if (my_chemistry->dust_species == 0) { + gasgr[i] = + my_rates->gas_grain[logTlininterp_buf.indixe[i] - 1] + + logTlininterp_buf.tdef[i] * + (my_rates->gas_grain[logTlininterp_buf.indixe[i] + 1 - 1] - + my_rates->gas_grain[logTlininterp_buf.indixe[i] - 1]); + + gasgr_tdust[i] = + (dust2gas[i] / metallicity[i]) * gasgr[i] * coolunit / mh_local_var; + + } else { + fv2k = my_rates->gas_grain2[logTlininterp_buf.indixe[i] - 1] + + logTlininterp_buf.tdef[i] * + (my_rates->gas_grain2[logTlininterp_buf.indixe[i] + 1 - 1] - + my_rates->gas_grain2[logTlininterp_buf.indixe[i] - 1]); + + fac = coolunit / mh_local_var; + + if (my_chemistry->use_multiple_dust_temperatures == 0) { + gasgr[i] = fv2k * internal_dust_prop_buf.sigma_per_gas_mass_tot[i]; + + gasgr_tdust[i] = gasgr[i] * fac; + + } else { + if (my_chemistry->dust_species > 0) { + gas_grainsp_heatrate.data[OnlyGrainSpLUT::MgSiO3_dust][i] = + fv2k * internal_dust_prop_buf.grain_sigma_per_gas_mass + .data[OnlyGrainSpLUT::MgSiO3_dust][i]; + gas_grainsp_heatrate.data[OnlyGrainSpLUT::AC_dust][i] = + fv2k * internal_dust_prop_buf.grain_sigma_per_gas_mass + .data[OnlyGrainSpLUT::AC_dust][i]; + } + if (my_chemistry->dust_species > 1) { + gas_grainsp_heatrate.data[OnlyGrainSpLUT::SiM_dust][i] = + fv2k * internal_dust_prop_buf.grain_sigma_per_gas_mass + .data[OnlyGrainSpLUT::SiM_dust][i]; + gas_grainsp_heatrate.data[OnlyGrainSpLUT::FeM_dust][i] = + fv2k * internal_dust_prop_buf.grain_sigma_per_gas_mass + .data[OnlyGrainSpLUT::FeM_dust][i]; + gas_grainsp_heatrate.data[OnlyGrainSpLUT::Mg2SiO4_dust][i] = + fv2k * internal_dust_prop_buf.grain_sigma_per_gas_mass + .data[OnlyGrainSpLUT::Mg2SiO4_dust][i]; + gas_grainsp_heatrate.data[OnlyGrainSpLUT::Fe3O4_dust][i] = + fv2k * internal_dust_prop_buf.grain_sigma_per_gas_mass + .data[OnlyGrainSpLUT::Fe3O4_dust][i]; + gas_grainsp_heatrate.data[OnlyGrainSpLUT::SiO2_dust][i] = + fv2k * internal_dust_prop_buf.grain_sigma_per_gas_mass + .data[OnlyGrainSpLUT::SiO2_dust][i]; + gas_grainsp_heatrate.data[OnlyGrainSpLUT::MgO_dust][i] = + fv2k * internal_dust_prop_buf.grain_sigma_per_gas_mass + .data[OnlyGrainSpLUT::MgO_dust][i]; + gas_grainsp_heatrate.data[OnlyGrainSpLUT::FeS_dust][i] = + fv2k * internal_dust_prop_buf.grain_sigma_per_gas_mass + .data[OnlyGrainSpLUT::FeS_dust][i]; + gas_grainsp_heatrate.data[OnlyGrainSpLUT::Al2O3_dust][i] = + fv2k * internal_dust_prop_buf.grain_sigma_per_gas_mass + .data[OnlyGrainSpLUT::Al2O3_dust][i]; + } + if (my_chemistry->dust_species > 2) { + gas_grainsp_heatrate.data[OnlyGrainSpLUT::ref_org_dust][i] = + fv2k * internal_dust_prop_buf.grain_sigma_per_gas_mass + .data[OnlyGrainSpLUT::ref_org_dust][i]; + gas_grainsp_heatrate.data[OnlyGrainSpLUT::vol_org_dust][i] = + fv2k * internal_dust_prop_buf.grain_sigma_per_gas_mass + .data[OnlyGrainSpLUT::vol_org_dust][i]; + gas_grainsp_heatrate.data[OnlyGrainSpLUT::H2O_ice_dust][i] = + fv2k * internal_dust_prop_buf.grain_sigma_per_gas_mass + .data[OnlyGrainSpLUT::H2O_ice_dust][i]; + } + + if (my_chemistry->dust_species > 0) { + gasgr_tMgSiO3[i] = + gas_grainsp_heatrate.data[OnlyGrainSpLUT::MgSiO3_dust][i] * fac; + gasgr_tAC[i] = + gas_grainsp_heatrate.data[OnlyGrainSpLUT::AC_dust][i] * fac; + } + if (my_chemistry->dust_species > 1) { + gasgr_tSiM[i] = + gas_grainsp_heatrate.data[OnlyGrainSpLUT::SiM_dust][i] * fac; + gasgr_tFeM[i] = + gas_grainsp_heatrate.data[OnlyGrainSpLUT::FeM_dust][i] * fac; + gasgr_tMg2SiO4[i] = + gas_grainsp_heatrate.data[OnlyGrainSpLUT::Mg2SiO4_dust][i] * + fac; + gasgr_tFe3O4[i] = + gas_grainsp_heatrate.data[OnlyGrainSpLUT::Fe3O4_dust][i] * fac; + gasgr_tSiO2D[i] = + gas_grainsp_heatrate.data[OnlyGrainSpLUT::SiO2_dust][i] * fac; + gasgr_tMgO[i] = + gas_grainsp_heatrate.data[OnlyGrainSpLUT::MgO_dust][i] * fac; + gasgr_tFeS[i] = + gas_grainsp_heatrate.data[OnlyGrainSpLUT::FeS_dust][i] * fac; + gasgr_tAl2O3[i] = + gas_grainsp_heatrate.data[OnlyGrainSpLUT::Al2O3_dust][i] * fac; + } + if (my_chemistry->dust_species > 2) { + gasgr_treforg[i] = + gas_grainsp_heatrate.data[OnlyGrainSpLUT::ref_org_dust][i] * + fac; + gasgr_tvolorg[i] = + gas_grainsp_heatrate.data[OnlyGrainSpLUT::vol_org_dust][i] * + fac; + gasgr_tH2Oice[i] = + gas_grainsp_heatrate.data[OnlyGrainSpLUT::H2O_ice_dust][i] * + fac; + } + } + } + } + } + + // Compute dust temperature + + if (my_chemistry->use_multiple_dust_temperatures == 0) { + grackle::impl::calc_tdust_1d_g( + tdust, tgas, nh, gasgr_tdust, mygisrf.data(), myisrf, itmask_metal, + trad, my_fields->grid_dimension[0], gr_N[1], &dlog10Tdust, + log10Tdust_vals, internal_dust_prop_buf.dyntab_kappa_tot, kptot, + &my_chemistry->dust_species, idx_range); + + } else { + if (my_chemistry->dust_species > 0) { + grackle::impl::calc_tdust_1d_g( + grain_temperatures.data[OnlyGrainSpLUT::MgSiO3_dust], tgas, nh, + gasgr_tMgSiO3.data(), gisrfMgSiO3.data(), myisrf, itmask_metal, trad, + my_fields->grid_dimension[0], gr_N[1], &dlog10Tdust, log10Tdust_vals, + internal_dust_prop_buf.grain_dyntab_kappa + .data[OnlyGrainSpLUT::MgSiO3_dust], + grain_kappa.data[OnlyGrainSpLUT::MgSiO3_dust], + &my_chemistry->dust_species, idx_range); + + grackle::impl::calc_tdust_1d_g( + grain_temperatures.data[OnlyGrainSpLUT::AC_dust], tgas, nh, + gasgr_tAC.data(), gisrfAC.data(), myisrf, itmask_metal, trad, + my_fields->grid_dimension[0], gr_N[1], &dlog10Tdust, log10Tdust_vals, + internal_dust_prop_buf.grain_dyntab_kappa + .data[OnlyGrainSpLUT::AC_dust], + grain_kappa.data[OnlyGrainSpLUT::AC_dust], + &my_chemistry->dust_species, idx_range); + } + + if (my_chemistry->dust_species > 1) { + grackle::impl::calc_tdust_1d_g( + grain_temperatures.data[OnlyGrainSpLUT::SiM_dust], tgas, nh, + gasgr_tSiM.data(), gisrfSiM.data(), myisrf, itmask_metal, trad, + my_fields->grid_dimension[0], gr_N[1], &dlog10Tdust, log10Tdust_vals, + internal_dust_prop_buf.grain_dyntab_kappa + .data[OnlyGrainSpLUT::SiM_dust], + grain_kappa.data[OnlyGrainSpLUT::SiM_dust], + &my_chemistry->dust_species, idx_range); + + grackle::impl::calc_tdust_1d_g( + grain_temperatures.data[OnlyGrainSpLUT::FeM_dust], tgas, nh, + gasgr_tFeM.data(), gisrfFeM.data(), myisrf, itmask_metal, trad, + my_fields->grid_dimension[0], gr_N[1], &dlog10Tdust, log10Tdust_vals, + internal_dust_prop_buf.grain_dyntab_kappa + .data[OnlyGrainSpLUT::FeM_dust], + grain_kappa.data[OnlyGrainSpLUT::FeM_dust], + &my_chemistry->dust_species, idx_range); + + grackle::impl::calc_tdust_1d_g( + grain_temperatures.data[OnlyGrainSpLUT::Mg2SiO4_dust], tgas, nh, + gasgr_tMg2SiO4.data(), gisrfMg2SiO4.data(), myisrf, itmask_metal, + trad, my_fields->grid_dimension[0], gr_N[1], &dlog10Tdust, + log10Tdust_vals, + internal_dust_prop_buf.grain_dyntab_kappa + .data[OnlyGrainSpLUT::Mg2SiO4_dust], + grain_kappa.data[OnlyGrainSpLUT::Mg2SiO4_dust], + &my_chemistry->dust_species, idx_range); + + grackle::impl::calc_tdust_1d_g( + grain_temperatures.data[OnlyGrainSpLUT::Fe3O4_dust], tgas, nh, + gasgr_tFe3O4.data(), gisrfFe3O4.data(), myisrf, itmask_metal, trad, + my_fields->grid_dimension[0], gr_N[1], &dlog10Tdust, log10Tdust_vals, + internal_dust_prop_buf.grain_dyntab_kappa + .data[OnlyGrainSpLUT::Fe3O4_dust], + grain_kappa.data[OnlyGrainSpLUT::Fe3O4_dust], + &my_chemistry->dust_species, idx_range); + + grackle::impl::calc_tdust_1d_g( + grain_temperatures.data[OnlyGrainSpLUT::SiO2_dust], tgas, nh, + gasgr_tSiO2D.data(), gisrfSiO2D.data(), myisrf, itmask_metal, trad, + my_fields->grid_dimension[0], gr_N[1], &dlog10Tdust, log10Tdust_vals, + internal_dust_prop_buf.grain_dyntab_kappa + .data[OnlyGrainSpLUT::SiO2_dust], + grain_kappa.data[OnlyGrainSpLUT::SiO2_dust], + &my_chemistry->dust_species, idx_range); + + grackle::impl::calc_tdust_1d_g( + grain_temperatures.data[OnlyGrainSpLUT::MgO_dust], tgas, nh, + gasgr_tMgO.data(), gisrfMgO.data(), myisrf, itmask_metal, trad, + my_fields->grid_dimension[0], gr_N[1], &dlog10Tdust, log10Tdust_vals, + internal_dust_prop_buf.grain_dyntab_kappa + .data[OnlyGrainSpLUT::MgO_dust], + grain_kappa.data[OnlyGrainSpLUT::MgO_dust], + &my_chemistry->dust_species, idx_range); + + grackle::impl::calc_tdust_1d_g( + grain_temperatures.data[OnlyGrainSpLUT::FeS_dust], tgas, nh, + gasgr_tFeS.data(), gisrfFeS.data(), myisrf, itmask_metal, trad, + my_fields->grid_dimension[0], gr_N[1], &dlog10Tdust, log10Tdust_vals, + internal_dust_prop_buf.grain_dyntab_kappa + .data[OnlyGrainSpLUT::FeS_dust], + grain_kappa.data[OnlyGrainSpLUT::FeS_dust], + &my_chemistry->dust_species, idx_range); + + grackle::impl::calc_tdust_1d_g( + grain_temperatures.data[OnlyGrainSpLUT::Al2O3_dust], tgas, nh, + gasgr_tAl2O3.data(), gisrfAl2O3.data(), myisrf, itmask_metal, trad, + my_fields->grid_dimension[0], gr_N[1], &dlog10Tdust, log10Tdust_vals, + internal_dust_prop_buf.grain_dyntab_kappa + .data[OnlyGrainSpLUT::Al2O3_dust], + grain_kappa.data[OnlyGrainSpLUT::Al2O3_dust], + &my_chemistry->dust_species, idx_range); + } + + if (my_chemistry->dust_species > 2) { + grackle::impl::calc_tdust_1d_g( + grain_temperatures.data[OnlyGrainSpLUT::ref_org_dust], tgas, nh, + gasgr_treforg.data(), gisrfreforg.data(), myisrf, itmask_metal, trad, + my_fields->grid_dimension[0], gr_N[1], &dlog10Tdust, log10Tdust_vals, + internal_dust_prop_buf.grain_dyntab_kappa + .data[OnlyGrainSpLUT::ref_org_dust], + grain_kappa.data[OnlyGrainSpLUT::ref_org_dust], + &my_chemistry->dust_species, idx_range); + + grackle::impl::calc_tdust_1d_g( + grain_temperatures.data[OnlyGrainSpLUT::vol_org_dust], tgas, nh, + gasgr_tvolorg.data(), gisrfvolorg.data(), myisrf, itmask_metal, trad, + my_fields->grid_dimension[0], gr_N[1], &dlog10Tdust, log10Tdust_vals, + internal_dust_prop_buf.grain_dyntab_kappa + .data[OnlyGrainSpLUT::vol_org_dust], + grain_kappa.data[OnlyGrainSpLUT::vol_org_dust], + &my_chemistry->dust_species, idx_range); + + grackle::impl::calc_tdust_1d_g( + grain_temperatures.data[OnlyGrainSpLUT::H2O_ice_dust], tgas, nh, + gasgr_tH2Oice.data(), gisrfH2Oice.data(), myisrf, itmask_metal, trad, + my_fields->grid_dimension[0], gr_N[1], &dlog10Tdust, log10Tdust_vals, + internal_dust_prop_buf.grain_dyntab_kappa + .data[OnlyGrainSpLUT::H2O_ice_dust], + grain_kappa.data[OnlyGrainSpLUT::H2O_ice_dust], + &my_chemistry->dust_species, idx_range); + } + } +} \ No newline at end of file diff --git a/src/clib/dust/calc_all_tdust_gasgr_1d_g.hpp b/src/clib/dust/calc_all_tdust_gasgr_1d_g.hpp new file mode 100644 index 000000000..713e39d71 --- /dev/null +++ b/src/clib/dust/calc_all_tdust_gasgr_1d_g.hpp @@ -0,0 +1,81 @@ +//===----------------------------------------------------------------------===// +// +// See the LICENSE file for license and copyright information +// SPDX-License-Identifier: NCSA AND BSD-3-Clause +// +//===----------------------------------------------------------------------===// +/// +/// @file +/// Declares the calc_all_tdust_gasgr_1d_g function +/// +//===----------------------------------------------------------------------===// + +// This file was initially generated automatically during conversion of the +// calc_all_tdust_gasgr_1d_g function from FORTRAN to C++ + +#ifndef CALC_ALL_TDUST_GASGR_1D_G_HPP +#define CALC_ALL_TDUST_GASGR_1D_G_HPP + +#include "dust_props.hpp" +#include "fortran_func_decls.h" +#include "grackle.h" +#include "index_helper.h" +#include "internal_types.hpp" + +namespace grackle::impl { + +/// Calculate all dust temperature(s) and the gas to grain heat +/// transfer rate(s) (the latter is commonly called gasgr) +/// +/// An argument could be made for directly using gasgr to compute +/// contributions to edot within this routine, rather than returning +/// the gasgr value(s). But that should be reconsidered in the future. +/// +/// We could significantly reduce the amount of buffer space allocated +/// within the routine. +/// +/// @param[in] trad CMB ratiation temperature +/// @param[in] tgas 1D array to hold the computed gas temperatures +/// @param[out] tdust 1D array to hold the computed dust temperatures +/// @param[in] metallicity 1D array to hold the computed metallicity +/// @param[out] dust2gas Holds the computed dust-to-gas ratio at each +/// location in the index range. In other words, this holds the dust mass +/// per unit gas mass (only used in certain configuration) +/// @param[in] nh 1D array of hydrogen number densities +/// @param[in] gasgr_tdust Array of gas-grain dust temperatures +/// @param[in] itmask_metal Mask array indicating presence of metals +/// @param[in] coolunit Cooling unit conversion factor +/// @param[out] gasgr Array of gas-grain heat transfer rates +/// @param[in] myisrf Array of interstellar radiation field strengths +/// @param[in] kptot Array of total grain opacities +/// @param[in] my_chemistry Pointer to chemistry data structure +/// @param[in] my_rates Pointer to chemistry data storage structure +/// @param[in] my_fields Pointer to field data structure +/// @param[in] idx_range Specifies the current index-range +/// @param[in] grain_temperatures buffers to hold individual grain species +/// temperatures. This is only used in certain configurations (i.e. when we +/// aren't using the tdust argument) +/// @param[in] gas_grainsp_heatrate Gas-grain species heat rate collection +/// @param[in] logTlininterp_buf Scratch space used to temporarily hold values +/// for each location in @p idx_range . +/// @param[in] internal_dust_prop_buf Internal dust properties buffer +/// @param[in] grain_kappa Grain species opacity collection +/// +/// @par History +/// modified: January, 2026 by Christopher Bignamini & Matthew Abruzzo; C++ port +void calc_all_tdust_gasgr_1d_g( + double trad, double* tgas, double* tdust, const double* metallicity, + const double* dust2gas, double* nh, double* gasgr_tdust, + gr_mask_type* itmask_metal, double coolunit, double* gasgr, double* myisrf, + double* kptot, chemistry_data* my_chemistry, + chemistry_data_storage* my_rates, grackle_field_data* my_fields, + IndexRange idx_range, + grackle::impl::GrainSpeciesCollection grain_temperatures, + grackle::impl::GrainSpeciesCollection gas_grainsp_heatrate, + grackle::impl::LogTLinInterpScratchBuf logTlininterp_buf, + grackle::impl::InternalDustPropBuf internal_dust_prop_buf, + grackle::impl::GrainSpeciesCollection grain_kappa); + +} // namespace grackle::impl + +#endif /* CALC_ALL_TDUST_GASGR_1D_G_HPP */ diff --git a/src/clib/dust/calc_tdust_1d_g.cpp b/src/clib/dust/calc_tdust_1d_g.cpp new file mode 100644 index 000000000..8a07fce47 --- /dev/null +++ b/src/clib/dust/calc_tdust_1d_g.cpp @@ -0,0 +1,316 @@ +//===----------------------------------------------------------------------===// +// +// See the LICENSE file for license and copyright information +// SPDX-License-Identifier: NCSA AND BSD-3-Clause +// +//===----------------------------------------------------------------------===// +/// +/// @file +/// Implements the calc_tdust_1d_g function +/// +//===----------------------------------------------------------------------===// + +// This file was initially generated automatically during conversion of the +// calc_tdust_1d_g function from FORTRAN to C++ + +#include +#include + +#include "grackle.h" +#include "fortran_func_decls.h" +#include "fortran_func_wrappers.hpp" +#include "utils-cpp.hpp" + +#include "calc_tdust_1d_g.hpp" + +void grackle::impl::calc_tdust_1d_g( + double* tdust, double* tgas, double* nh, double* gasgr, + const double* gamma_isrfa, const double* isrf, const gr_mask_type* itmask, + double trad, int in, int gr_N, double* gr_dT, double* gr_Td, + gr_float* alsp_data_, double* kgr, int* idspecies, IndexRange idx_range) { + // opacity table of a grain species + // + // In some configurations gr_N can be 0 while the backing buffer may still be + // non-null. The View invariant disallows non-null data with a zero leading + // extent, so pass nullptr for the zero-length case. + gr_float* alsp_ptr = (gr_N > 0) ? alsp_data_ : nullptr; + grackle::impl::View alsp(alsp_ptr, gr_N, in); + std::vector logalsp_data_(gr_N * in); + gr_float* logalsp_ptr = (gr_N > 0) ? logalsp_data_.data() : nullptr; + grackle::impl::View logalsp(logalsp_ptr, gr_N, in); + int Td_Size; + int Td_N; + + // Parameters + + // grain sublimation temperature + double t_subl = 1.5e3; // TODO: should be const + const double radf = 4.; + + // grain opacity from Omukai (2000, equation 17) normalized by + // the local dust-to-gas ratio, which in this work is 0.934e-2. + const double kgr1 = 4.0e-4; + + std::vector gamma_isrf(in); + const double tol = 1.e-5; + const double bi_tol = 1.e-3; + const double minpert = 1.e-10; + const int itmax = 50; + const int bi_itmax = 30; + + // Locals + + int i, iter, c_done, c_total, nm_done; + + double pert_i, trad4; + + // Slice Locals + + std::vector kgrplus(in); + std::vector sol(in); + std::vector solplus(in); + std::vector slope(in); + std::vector tdplus(in); + std::vector tdustnow(in); + std::vector tdustold(in); + std::vector pert(in); + std::vector bi_t_mid(in); + std::vector bi_t_high(in); + std::vector nm_itmask(in); + std::vector bi_itmask(in); + + pert_i = 1.e-3; + + trad = std::fmax( + 1., + trad); // TODO: do we really want to write in a passed by copy input par? + trad4 = std::pow(trad, 4); + + // \sum rho_SN kappa_SN / \sum rho_SN ndust_SN + + Td_N = gr_N; + Td_Size = gr_N; + for (i = idx_range.i_start; i <= idx_range.i_end; i++) { + if (itmask[i] != MASK_FALSE) { + for (int j = 0; j < gr_N; j++) { + logalsp(j, i) = std::log10(alsp(j, i)); + } + } + } + + // Set total cells for calculation + + c_done = 0; + nm_done = 0; + c_total = idx_range.i_end - idx_range.i_start + 1; + + // Set local iteration mask and initial guess + + for (i = idx_range.i_start; i <= idx_range.i_end; i++) { + if (itmask[i] != MASK_FALSE) { + gamma_isrf[i] = isrf[i] * gamma_isrfa[i]; + } + } + + for (i = idx_range.i_start; i <= idx_range.i_end; i++) { + nm_itmask[i] = itmask[i]; + bi_itmask[i] = itmask[i]; + if (nm_itmask[i] != MASK_FALSE) { + if (trad >= tgas[i]) { + tdustnow[i] = trad; + nm_itmask[i] = MASK_FALSE; + bi_itmask[i] = MASK_FALSE; + c_done = c_done + 1; + nm_done = nm_done + 1; + } else if (tgas[i] > t_subl) { + // Use bisection if T_gas > grain sublimation temperature. + nm_itmask[i] = MASK_FALSE; + nm_done = nm_done + 1; + } else { + tdustnow[i] = + std::fmax(trad, std::pow((gamma_isrf[i] / radf / kgr1), 0.17)); + pert[i] = pert_i; + } + + } else { + c_done = c_done + 1; + nm_done = nm_done + 1; + } + } + + // Iterate to convergence with Newton's method + + for (iter = 1; iter <= (itmax); iter++) { + // Loop over slice + + for (i = idx_range.i_start; i <= idx_range.i_end; i++) { + if (nm_itmask[i] != MASK_FALSE) { + tdplus[i] = std::fmax(1.e-3, ((1. + pert[i]) * tdustnow[i])); + } + } + + // Calculate grain opacities + + FORTRAN_NAME(calc_kappa_gr_g)(tdustnow.data(), kgr, nm_itmask.data(), &in, + &idx_range.i_start, &idx_range.i_end, &t_subl, + &Td_N, &Td_Size, gr_dT, gr_Td, logalsp.data(), + idspecies); + + FORTRAN_NAME(calc_kappa_gr_g)(tdplus.data(), kgrplus.data(), + nm_itmask.data(), &in, &idx_range.i_start, + &idx_range.i_end, &t_subl, &Td_N, &Td_Size, + gr_dT, gr_Td, logalsp.data(), idspecies); + + // Calculate heating/cooling balance + + FORTRAN_NAME(calc_gr_balance_g)(tdustnow.data(), tgas, kgr, &trad4, gasgr, + gamma_isrf.data(), nh, nm_itmask.data(), + sol.data(), &in, &idx_range.i_start, + &idx_range.i_end); + + FORTRAN_NAME(calc_gr_balance_g)(tdplus.data(), tgas, kgrplus.data(), &trad4, + gasgr, gamma_isrf.data(), nh, + nm_itmask.data(), solplus.data(), &in, + &idx_range.i_start, &idx_range.i_end); + + for (i = idx_range.i_start; i <= idx_range.i_end; i++) { + if (nm_itmask[i] != MASK_FALSE) { + // Use Newton's method to solve for Tdust + + slope[i] = (solplus[i] - sol[i]) / (pert[i] * tdustnow[i]); + + tdustold[i] = tdustnow[i]; + // tdustnow(i) = tdustnow(i) - (sol(i) / slope(i)) + tdustnow[i] = std::fmin(tdustnow[i] - (sol[i] / slope[i]), 3e3); + + pert[i] = std::fmax( + std::fmin(pert[i], (0.5 * std::fabs(tdustnow[i] - tdustold[i]) / + tdustnow[i])), + minpert); + + // If negative solution calculated, give up and wait for bisection step. + if (tdustnow[i] < trad) { + nm_itmask[i] = MASK_FALSE; + nm_done = nm_done + 1; + // Check for convergence of solution + } else if (std::fabs(sol[i] / solplus[i]) < tol) { + nm_itmask[i] = MASK_FALSE; + c_done = c_done + 1; + bi_itmask[i] = MASK_FALSE; + nm_done = nm_done + 1; + } + + // if ( nm_itmask(i) ) + } + + // End loop over slice + } + + // Check for all cells converged + if (c_done >= c_total) { + break; + } + + // Check for all cells done with Newton method + // This includes attempts where a negative solution was found + if (nm_done >= c_total) { + break; + } + + // End iteration loop for Newton's method + } + + // If iteration count exceeded, try once more with bisection + if (c_done < c_total) { + for (i = idx_range.i_start; i <= idx_range.i_end; i++) { + if (bi_itmask[i] != MASK_FALSE) { + tdustnow[i] = trad; + // bi_t_high(i) = tgas(i) + bi_t_high[i] = 3e3; + } + } + + for (iter = 1; iter <= (bi_itmax); iter++) { + for (i = idx_range.i_start; i <= idx_range.i_end; i++) { + if (bi_itmask[i] != MASK_FALSE) { + bi_t_mid[i] = 0.5 * (tdustnow[i] + bi_t_high[i]); + if (iter == 1) { + bi_t_mid[i] = std::fmin(bi_t_mid[i], t_subl); + } + } + } + + FORTRAN_NAME(calc_kappa_gr_g)(bi_t_mid.data(), kgr, bi_itmask.data(), &in, + &idx_range.i_start, &idx_range.i_end, + &t_subl, &Td_N, &Td_Size, gr_dT, gr_Td, + logalsp.data(), idspecies); + + FORTRAN_NAME(calc_gr_balance_g)(bi_t_mid.data(), tgas, kgr, &trad4, gasgr, + gamma_isrf.data(), nh, bi_itmask.data(), + sol.data(), &in, &idx_range.i_start, + &idx_range.i_end); + + for (i = idx_range.i_start; i <= idx_range.i_end; i++) { + if (bi_itmask[i] != MASK_FALSE) { + if (sol[i] > 0.) { + tdustnow[i] = bi_t_mid[i]; + } else { + bi_t_high[i] = bi_t_mid[i]; + } + + if ((std::fabs(bi_t_high[i] - tdustnow[i]) / tdustnow[i]) <= bi_tol) { + bi_itmask[i] = MASK_FALSE; + c_done = c_done + 1; + } + + // Check for all cells converged + if (c_done >= c_total) { + break; + } + + // if ( bi_itmask(i) ) + } + + // End loop over slice + } + + if (c_done >= c_total) { + break; + } + // End iteration loop for bisection + } + + // If iteration count exceeded with bisection, end of the line. + if (iter > itmax) { + OMP_PRAGMA_CRITICAL { + eprintf( + "CALC_TDUST_1D_G failed to converge after %d iterations for %d " + "cells.\n", + iter, (c_total - c_done)); + } + } + + // if (iter .gt. itmax) then + } + + // Copy values back to thrown slice + for (i = idx_range.i_start; i <= idx_range.i_end; i++) { + if (itmask[i] != MASK_FALSE) { + // Check for bad solutions + if (tdustnow[i] < 0.) { + OMP_PRAGMA_CRITICAL { + eprintf( + "CALC_TDUST_1D_G Newton method - T_dust < 0: i = %d j = %d k " + "= %d nh = %g t_gas = %g t_rad = %g t_dust = %g\n", + i, idx_range.jp1, idx_range.kp1, nh[i], tgas[i], trad, + tdustnow[i]); + } + // ERROR_MESSAGE + } + + tdust[i] = tdustnow[i]; + } + } + + return; +} diff --git a/src/clib/dust/calc_tdust_1d_g.hpp b/src/clib/dust/calc_tdust_1d_g.hpp new file mode 100644 index 000000000..79b01941b --- /dev/null +++ b/src/clib/dust/calc_tdust_1d_g.hpp @@ -0,0 +1,62 @@ +//===----------------------------------------------------------------------===// +// +// See the LICENSE file for license and copyright information +// SPDX-License-Identifier: NCSA AND BSD-3-Clause +// +//===----------------------------------------------------------------------===// +/// +/// @file +/// Declares the calc_tdust_1d_g function +/// +//===----------------------------------------------------------------------===// + +// This file was initially generated automatically during conversion of the +// calc_tdust_1d_g function from FORTRAN to C++ + +#ifndef CALC_TDUST_1D_G_HPP +#define CALC_TDUST_1D_G_HPP + +#include "fortran_func_decls.h" // gr_mask_int +#include "grackle.h" // gr_float +#include "index_helper.h" + +namespace grackle::impl { + +/// Calculate equilibrium dust temperature +/// +/// TODO: this docstring, and the docstrings of all helper functions should +/// EXPLICITLY document how the meaning of arguments change based on +/// whether we are using the single-field dust model or the +/// multi-species dust model. The different meaning of each variable +/// gets VERY confusing +/// +/// @param[out] tdust 1D array to hold dust temperatures +/// @param[in] tgas 1D array to hold gas temperatures +/// @param[in] nh 1D array of hydrogen number densities +/// @param[in] gasgr Array of gas-grain heat transfer rates +/// @param[in] gamma_isrfa Heating from interstellar radiation field +/// @param[in] isrf Interstellar radiation field strength in Habing units +/// @param[in] itmask Itaration mask +/// @param[in] trad CMB ratiation temperature +/// @param[in] in Length of the 1D slice +/// @param[in] gr_N Number of temperature points in the grain opacity table +/// @param[in] gr_dT Temperature spacing of the grain opacity table +/// @param[in] gr_Td Temperature values of the grain opacity table +/// @param[in] alsp_data_ Grain opacity table data +/// @param[in] kgr Array to hold computed grain opacities +/// @param[in] idspecies Array of grain species IDs (only used in certain +/// configurations) +/// @param[in] idx_range Index range specifying the portion of the grid to +/// operate on +/// +/// @par History +/// written by: Britton Smith, 2011 +/// modified: January, 2026 by Christopher Bignamini & Matthew Abruzzo; C++ port +void calc_tdust_1d_g(double* tdust, double* tgas, double* nh, double* gasgr, + const double* gamma_isrfa, const double* isrf, + const gr_mask_type* itmask, double trad, int in, int gr_N, + double* gr_dT, double* gr_Td, gr_float* alsp_data_, + double* kgr, int* idspecies, IndexRange idx_range); + +} // namespace grackle::impl +#endif /* CALC_TDUST_1D_G_HPP */ diff --git a/src/clib/fortran_func_decls.h b/src/clib/fortran_func_decls.h index 42422e16d..1b9de9da2 100644 --- a/src/clib/fortran_func_decls.h +++ b/src/clib/fortran_func_decls.h @@ -19,34 +19,6 @@ typedef long long gr_i64; #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_species_1d)( const int* igrgr, const gr_mask_type* itmask, const int* SN0_N, int* in, int* jn, int* kn, int* is, int* ie, int* j, int* k, double* dom, gr_float* d_data_ptr, int* nSN, @@ -60,13 +32,6 @@ void FORTRAN_NAME(solve_cubic_equation)( double* a, double* b, double* c, double* root ); -void FORTRAN_NAME(calc_tdust_1d_g)( - double* tdust, double* tgas, double* nh, double* gasgr, double* gamma_isrfa, - double* isrf, gr_mask_type* itmask, double* trad, int* in, int* is, int* ie, - int* j, int* k, int* gr_N, int* gr_Size, double* gr_dT, double* gr_Td, - gr_float* alsp_data_ptr, double* kgr, int* idspecies -); - void FORTRAN_NAME(calc_kappa_gr_g)( double* tdust, double* kgr, gr_mask_type* itmask, int* in, int* is, int* ie, double* t_subl, int* gr_N, int* gr_Size, double* gr_dT, double* gr_Td, diff --git a/src/clib/fortran_func_wrappers.hpp b/src/clib/fortran_func_wrappers.hpp index 4bbc921e8..df7a509ec 100644 --- a/src/clib/fortran_func_wrappers.hpp +++ b/src/clib/fortran_func_wrappers.hpp @@ -42,71 +42,6 @@ // namespace name when they call these routines namespace grackle::impl::fortran_wrapper { -inline void calc_all_tdust_gasgr_1d_g( - double trad, double* tgas, double* tdust, double* metallicity, - double* dust2gas, double* nh, double* gasgr_tdust, gr_mask_type* itmask_metal, - double coolunit, double* gasgr, double* myisrf, double* kappa_tot, - chemistry_data* my_chemistry, chemistry_data_storage* my_rates, - grackle_field_data* my_fields, IndexRange idx_range, - grackle::impl::GrainSpeciesCollection grain_temperatures, - grackle::impl::GrainSpeciesCollection gas_grainsp_heatrate, - grackle::impl::GrainSpeciesCollection grain_kappa, - grackle::impl::LogTLinInterpScratchBuf logTlininterp_buf, - grackle::impl::InternalDustPropBuf internal_dust_prop_buf -) { - - // after transcription, we should obviously move this logic inside of - // the transcribed function - grackle::impl::GrainMetalInjectPathways* inject_pathway_props = - my_rates->opaque_storage->inject_pathway_props; - - double dlog10Tdust = 0.0; - double* log10Tdust_vals = nullptr; - - // NOTE: gr_N and gr_Size are historical names - // -> they are pretty uninformative and should be changed! - int gr_N[2] = {0, 0}; - int gr_Size = 0; - if (inject_pathway_props != nullptr) { - dlog10Tdust = - inject_pathway_props->log10Tdust_interp_props.parameter_spacing[0]; - log10Tdust_vals = - inject_pathway_props->log10Tdust_interp_props.parameters[0]; - - gr_N[0] = inject_pathway_props->n_opac_poly_coef; - gr_N[1] = static_cast( - inject_pathway_props->log10Tdust_interp_props.dimension[0]); - }; - gr_Size = gr_N[0] * gr_N[1]; - - - FORTRAN_NAME(calc_all_tdust_gasgr_1d_g)( - &my_fields->grid_dimension[0], &my_fields->grid_dimension[1], &my_fields->grid_dimension[2], &my_chemistry->NumberOfTemperatureBins, - &my_chemistry->use_dust_density_field, &idx_range.i_start, &idx_range.i_end, &idx_range.jp1, &idx_range.kp1, &my_chemistry->local_dust_to_gas_ratio, &my_rates->gamma_isrf, - &trad, my_rates->gas_grain, logTlininterp_buf.indixe, logTlininterp_buf.tdef, tgas, tdust, - metallicity, dust2gas, nh, gasgr_tdust, - itmask_metal, - &my_chemistry->dust_species, &my_chemistry->use_multiple_dust_temperatures, - gr_N, &gr_Size, &dlog10Tdust, log10Tdust_vals, - grain_temperatures.data[OnlyGrainSpLUT::SiM_dust], grain_temperatures.data[OnlyGrainSpLUT::FeM_dust], grain_temperatures.data[OnlyGrainSpLUT::Mg2SiO4_dust], grain_temperatures.data[OnlyGrainSpLUT::MgSiO3_dust], grain_temperatures.data[OnlyGrainSpLUT::Fe3O4_dust], - grain_temperatures.data[OnlyGrainSpLUT::AC_dust], grain_temperatures.data[OnlyGrainSpLUT::SiO2_dust], grain_temperatures.data[OnlyGrainSpLUT::MgO_dust], grain_temperatures.data[OnlyGrainSpLUT::FeS_dust], grain_temperatures.data[OnlyGrainSpLUT::Al2O3_dust], grain_temperatures.data[OnlyGrainSpLUT::ref_org_dust], - grain_temperatures.data[OnlyGrainSpLUT::vol_org_dust], grain_temperatures.data[OnlyGrainSpLUT::H2O_ice_dust], my_rates->gas_grain2, &my_rates->gamma_isrf2, - &coolunit, gasgr, myisrf, internal_dust_prop_buf.grain_sigma_per_gas_mass.data[OnlyGrainSpLUT::SiM_dust], internal_dust_prop_buf.grain_sigma_per_gas_mass.data[OnlyGrainSpLUT::FeM_dust], internal_dust_prop_buf.grain_sigma_per_gas_mass.data[OnlyGrainSpLUT::Mg2SiO4_dust], - internal_dust_prop_buf.grain_sigma_per_gas_mass.data[OnlyGrainSpLUT::MgSiO3_dust], internal_dust_prop_buf.grain_sigma_per_gas_mass.data[OnlyGrainSpLUT::Fe3O4_dust], internal_dust_prop_buf.grain_sigma_per_gas_mass.data[OnlyGrainSpLUT::AC_dust], internal_dust_prop_buf.grain_sigma_per_gas_mass.data[OnlyGrainSpLUT::SiO2_dust], internal_dust_prop_buf.grain_sigma_per_gas_mass.data[OnlyGrainSpLUT::MgO_dust], internal_dust_prop_buf.grain_sigma_per_gas_mass.data[OnlyGrainSpLUT::FeS_dust], - internal_dust_prop_buf.grain_sigma_per_gas_mass.data[OnlyGrainSpLUT::Al2O3_dust], internal_dust_prop_buf.grain_sigma_per_gas_mass.data[OnlyGrainSpLUT::ref_org_dust], internal_dust_prop_buf.grain_sigma_per_gas_mass.data[OnlyGrainSpLUT::vol_org_dust], internal_dust_prop_buf.grain_sigma_per_gas_mass.data[OnlyGrainSpLUT::H2O_ice_dust], internal_dust_prop_buf.sigma_per_gas_mass_tot, - internal_dust_prop_buf.grain_dyntab_kappa.data[OnlyGrainSpLUT::SiM_dust], internal_dust_prop_buf.grain_dyntab_kappa.data[OnlyGrainSpLUT::FeM_dust], internal_dust_prop_buf.grain_dyntab_kappa.data[OnlyGrainSpLUT::Mg2SiO4_dust], internal_dust_prop_buf.grain_dyntab_kappa.data[OnlyGrainSpLUT::MgSiO3_dust], internal_dust_prop_buf.grain_dyntab_kappa.data[OnlyGrainSpLUT::Fe3O4_dust], internal_dust_prop_buf.grain_dyntab_kappa.data[OnlyGrainSpLUT::AC_dust], - internal_dust_prop_buf.grain_dyntab_kappa.data[OnlyGrainSpLUT::SiO2_dust], internal_dust_prop_buf.grain_dyntab_kappa.data[OnlyGrainSpLUT::MgO_dust], internal_dust_prop_buf.grain_dyntab_kappa.data[OnlyGrainSpLUT::FeS_dust], internal_dust_prop_buf.grain_dyntab_kappa.data[OnlyGrainSpLUT::Al2O3_dust], internal_dust_prop_buf.grain_dyntab_kappa.data[OnlyGrainSpLUT::ref_org_dust], - internal_dust_prop_buf.grain_dyntab_kappa.data[OnlyGrainSpLUT::vol_org_dust], internal_dust_prop_buf.grain_dyntab_kappa.data[OnlyGrainSpLUT::H2O_ice_dust], internal_dust_prop_buf.dyntab_kappa_tot, grain_kappa.data[OnlyGrainSpLUT::SiM_dust], grain_kappa.data[OnlyGrainSpLUT::FeM_dust], - grain_kappa.data[OnlyGrainSpLUT::Mg2SiO4_dust], grain_kappa.data[OnlyGrainSpLUT::MgSiO3_dust], grain_kappa.data[OnlyGrainSpLUT::Fe3O4_dust], grain_kappa.data[OnlyGrainSpLUT::AC_dust], grain_kappa.data[OnlyGrainSpLUT::SiO2_dust], - grain_kappa.data[OnlyGrainSpLUT::MgO_dust], grain_kappa.data[OnlyGrainSpLUT::FeS_dust], grain_kappa.data[OnlyGrainSpLUT::Al2O3_dust], grain_kappa.data[OnlyGrainSpLUT::ref_org_dust], grain_kappa.data[OnlyGrainSpLUT::vol_org_dust], - grain_kappa.data[OnlyGrainSpLUT::H2O_ice_dust], kappa_tot, gas_grainsp_heatrate.data[OnlyGrainSpLUT::SiM_dust], gas_grainsp_heatrate.data[OnlyGrainSpLUT::FeM_dust], gas_grainsp_heatrate.data[OnlyGrainSpLUT::Mg2SiO4_dust], - gas_grainsp_heatrate.data[OnlyGrainSpLUT::MgSiO3_dust], gas_grainsp_heatrate.data[OnlyGrainSpLUT::Fe3O4_dust], gas_grainsp_heatrate.data[OnlyGrainSpLUT::AC_dust], gas_grainsp_heatrate.data[OnlyGrainSpLUT::SiO2_dust], gas_grainsp_heatrate.data[OnlyGrainSpLUT::MgO_dust], - gas_grainsp_heatrate.data[OnlyGrainSpLUT::FeS_dust], gas_grainsp_heatrate.data[OnlyGrainSpLUT::Al2O3_dust], gas_grainsp_heatrate.data[OnlyGrainSpLUT::ref_org_dust], gas_grainsp_heatrate.data[OnlyGrainSpLUT::vol_org_dust], - gas_grainsp_heatrate.data[OnlyGrainSpLUT::H2O_ice_dust] - ); - -} - /// Performs Gauss-Jordan elimination to solve the specified system of linear /// equations and inverts the coefficient matrix. ///