From 885e2fbdfdb82a3f69a752b9ac3dbd13a6e127ac Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 19 Jul 2024 08:35:20 -0400 Subject: [PATCH 1/5] call scale_fields_table_g outside of calculate_temperature --- src/clib/calc_temp_cloudy_g.F | 21 +++-------------- src/clib/calculate_temperature.c | 39 ++++++++++++++++++++++++++++++++ 2 files changed, 42 insertions(+), 18 deletions(-) diff --git a/src/clib/calc_temp_cloudy_g.F b/src/clib/calc_temp_cloudy_g.F index ac80f6f2..5f79c596 100644 --- a/src/clib/calc_temp_cloudy_g.F +++ b/src/clib/calc_temp_cloudy_g.F @@ -115,15 +115,10 @@ subroutine calc_temp_cloudy_g(d, e, metal, temperature, dbase1 = urho*(aye*uaye)**3 ! urho is [dens]/a^3 = [dens]/([a]*a')^3 ' zr = 1._DKIND/(aye*uaye) - 1._DKIND -! Convert densities from comoving to proper +! We explicitly assume that the densities have already been converted +! from comoving to proper (if iexpand .eq. 1), before this function is +! called. (Historically that logic was handled here. - if (iexpand .eq. 1) then - - call scale_fields_table_g(d, metal, - & is, ie, js, je, ks, ke, - & in, jn, kn, imetal, aye**(-3)) - - endif ! Loop over zones, and do an entire i-column in one go @@ -173,16 +168,6 @@ subroutine calc_temp_cloudy_g(d, e, metal, temperature, !$omp end parallel do #endif -! Convert densities back to comoving from proper - - if (iexpand .eq. 1) then - - call scale_fields_table_g(d, metal, - & is, ie, js, je, ks, ke, - & in, jn, kn, imetal, aye**3) - - endif - return end diff --git a/src/clib/calculate_temperature.c b/src/clib/calculate_temperature.c index 4466db60..37f9cbd5 100644 --- a/src/clib/calculate_temperature.c +++ b/src/clib/calculate_temperature.c @@ -47,6 +47,11 @@ extern void FORTRAN_NAME(calc_temp_cloudy_g)( double *priPar1, double *priPar2, double *priPar3, long long *priDataSize, double *priMMW); +extern void FORTRAN_NAME(scale_fields_table_g)( + gr_float* d, gr_float* metal, + int *is, int *ie, int *js, int *je, int *ks, int *ke, + int *in, int *jn, int *kn, int *imetal, double *factor); + double get_temperature_units(code_units *my_units); int local_calculate_pressure(chemistry_data *my_chemistry, @@ -147,6 +152,27 @@ int local_calculate_temperature(chemistry_data *my_chemistry, return SUCCESS; } +// This routine scales the density (and possibly metal) field(s) from +// comoving units to proper units (and back again) +static void scale_fields_table_g_cversion(grackle_field_data* my_fields, + int imetal, double factor) +{ + int in = my_fields->grid_dimension[0]; + int jn = my_fields->grid_dimension[1]; + int kn = my_fields->grid_dimension[2]; + + int is = my_fields->grid_start[0]; + int js = my_fields->grid_start[1]; + int ks = my_fields->grid_start[2]; + + int ie = my_fields->grid_end[0]; + int je = my_fields->grid_end[1]; + int ke = my_fields->grid_end[2]; + FORTRAN_NAME(scale_fields_table_g)(my_fields->density, my_fields->metal_density, + &is, &ie, &js, &je, &ks, &ke, + &in, &jn, &kn, &imetal, &factor); +} + int local_calculate_temperature_table(chemistry_data *my_chemistry, chemistry_data_storage *my_rates, code_units *my_units, @@ -184,6 +210,12 @@ int local_calculate_temperature_table(chemistry_data *my_chemistry, double temperature_units = get_temperature_units(my_units); + if (my_units->comoving_coordinates == 1) { + // convert density (& possibly metal) field(s) from comoving to proper + scale_fields_table_g_cversion(my_fields, metal_field_present, + pow(my_units->a_value, -3.0)); + } + FORTRAN_NAME(calc_temp_cloudy_g)( my_fields->density, my_fields->internal_energy, @@ -218,6 +250,13 @@ int local_calculate_temperature_table(chemistry_data *my_chemistry, &my_rates->cloudy_primordial.data_size, my_rates->cloudy_primordial.mmw_data); + + if (my_units->comoving_coordinates == 1) { + // convert density (& possibly metal) field(s) back to comoving from proper + scale_fields_table_g_cversion(my_fields, metal_field_present, + pow(my_units->a_value, 3.0)); + } + return SUCCESS; } From 1508037f66c5b024f9c5e0022420dc76c580d565 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 19 Jul 2024 08:50:02 -0400 Subject: [PATCH 2/5] completely transcribe calc_temp_cloudy_g from Fortran to C --- src/clib/calc_temp_cloudy_g.F | 45 -------------------------------- src/clib/calculate_temperature.c | 40 ++++++++++++++-------------- 2 files changed, 21 insertions(+), 64 deletions(-) diff --git a/src/clib/calc_temp_cloudy_g.F b/src/clib/calc_temp_cloudy_g.F index 5f79c596..27dccdd2 100644 --- a/src/clib/calc_temp_cloudy_g.F +++ b/src/clib/calc_temp_cloudy_g.F @@ -170,48 +170,3 @@ subroutine calc_temp_cloudy_g(d, e, metal, temperature, return end - -c ----------------------------------------------------------- -! This routine scales the density fields from comoving to -! proper densities (and back again). - - subroutine scale_fields_table_g(d, metal, - & is, ie, js, je, ks, ke, - & in, jn, kn, imetal, factor) -c ------------------------------------------------------------------- - - implicit NONE -#include "grackle_fortran_types.def" - -! Arguments - - integer in, jn, kn, is, ie, js, je, ks, ke, imetal - R_PREC d(in,jn,kn), metal(in,jn,kn) - real*8 factor - -! locals - - integer i, j, k - -! Multiply all fields by factor (1/a^3 or a^3) - - do k = ks+1, ke+1 - do j = js+1, je+1 - do i = is+1, ie+1 - d(i,j,k) = d(i,j,k) * factor - enddo - enddo - enddo - - if (imetal .eq. 1) then - do k = ks+1, ke+1 - do j = js+1, je+1 - do i = is+1, ie+1 - metal(i,j,k) = metal(i,j,k) * factor - enddo - enddo - enddo - endif - - return - end diff --git a/src/clib/calculate_temperature.c b/src/clib/calculate_temperature.c index 37f9cbd5..c2e582bc 100644 --- a/src/clib/calculate_temperature.c +++ b/src/clib/calculate_temperature.c @@ -47,11 +47,6 @@ extern void FORTRAN_NAME(calc_temp_cloudy_g)( double *priPar1, double *priPar2, double *priPar3, long long *priDataSize, double *priMMW); -extern void FORTRAN_NAME(scale_fields_table_g)( - gr_float* d, gr_float* metal, - int *is, int *ie, int *js, int *je, int *ks, int *ke, - int *in, int *jn, int *kn, int *imetal, double *factor); - double get_temperature_units(code_units *my_units); int local_calculate_pressure(chemistry_data *my_chemistry, @@ -157,20 +152,27 @@ int local_calculate_temperature(chemistry_data *my_chemistry, static void scale_fields_table_g_cversion(grackle_field_data* my_fields, int imetal, double factor) { - int in = my_fields->grid_dimension[0]; - int jn = my_fields->grid_dimension[1]; - int kn = my_fields->grid_dimension[2]; - - int is = my_fields->grid_start[0]; - int js = my_fields->grid_start[1]; - int ks = my_fields->grid_start[2]; - - int ie = my_fields->grid_end[0]; - int je = my_fields->grid_end[1]; - int ke = my_fields->grid_end[2]; - FORTRAN_NAME(scale_fields_table_g)(my_fields->density, my_fields->metal_density, - &is, &ie, &js, &je, &ks, &ke, - &in, &jn, &kn, &imetal, &factor); + /* Compute properties used to index the field. */ + const grackle_index_helper ind_helper = _build_index_helper(my_fields); + + for (int outer_i = 0; outer_i < ind_helper.outer_ind_size; outer_i++) { + const grackle_index_range range = _inner_range(outer_i, &ind_helper); + for (int i = range.start; i <= range.end; i++) { + my_fields->density[i] *= factor; + } // end: loop over i + } // end: loop over outer_ind + + if (imetal) { + + for (int outer_i = 0; outer_i < ind_helper.outer_ind_size; outer_i++) { + const grackle_index_range range = _inner_range(outer_i, &ind_helper); + for (int i = range.start; i <= range.end; i++) { + my_fields->metal_density[i] *= factor; + } // end: loop over i + } // end: loop over outer_ind + + } + } int local_calculate_temperature_table(chemistry_data *my_chemistry, From efc791d57bc8ff5de3972be27ba22320ab5697cd Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 19 Jul 2024 08:53:45 -0400 Subject: [PATCH 3/5] remove some unused variables from calc_temp_cloudy_g --- src/clib/calc_temp_cloudy_g.F | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/clib/calc_temp_cloudy_g.F b/src/clib/calc_temp_cloudy_g.F index 27dccdd2..962ead5a 100644 --- a/src/clib/calc_temp_cloudy_g.F +++ b/src/clib/calc_temp_cloudy_g.F @@ -94,7 +94,6 @@ subroutine calc_temp_cloudy_g(d, e, metal, temperature, integer i, j, k integer t, dj, dk real*8 dom, zr - real*8 dbase1, tbase1, xbase1 ! row temporaries @@ -110,9 +109,6 @@ subroutine calc_temp_cloudy_g(d, e, metal, temperature, ! Set units dom = urho*(aye**3)/mh - tbase1 = utim - xbase1 = uxyz/(aye*uaye) ! uxyz is [x]*a = [x]*[a]*a' ' - dbase1 = urho*(aye*uaye)**3 ! urho is [dens]/a^3 = [dens]/([a]*a')^3 ' zr = 1._DKIND/(aye*uaye) - 1._DKIND ! We explicitly assume that the densities have already been converted From 440edfbad7d7d51f4d14093422cd79f9fda036f5 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 19 Jul 2024 23:13:09 -0400 Subject: [PATCH 4/5] shift 2 variables out of calc_temp_cloudy_g --- src/clib/calc_temp_cloudy_g.F | 10 ++++------ src/clib/calculate_temperature.c | 7 ++++++- 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/src/clib/calc_temp_cloudy_g.F b/src/clib/calc_temp_cloudy_g.F index 962ead5a..004de677 100644 --- a/src/clib/calc_temp_cloudy_g.F +++ b/src/clib/calc_temp_cloudy_g.F @@ -6,6 +6,7 @@ subroutine calc_temp_cloudy_g(d, e, metal, temperature, & in, jn, kn, iexpand, imetal, & is, js, ks, ie, je, ke, + & dom, zr, & aye, temstart, temend, & utem, uxyz, uaye, urho, utim, & gamma, fh, @@ -34,6 +35,8 @@ subroutine calc_temp_cloudy_g(d, e, metal, temperature, ! iexpand - comoving coordinates flag (0 = off, 1 = on) ! imetal - flag if metal field is active (0 = no, 1 = yes) ! +! dom - unit conversion to proper number density in code units +! zr - current redshift ! fh - Hydrogen mass fraction (typically 0.76) ! aye - expansion factor (in code units) ! @@ -71,6 +74,7 @@ subroutine calc_temp_cloudy_g(d, e, metal, temperature, & iexpand, imetal real*8 aye, temstart, temend, gamma, & utim, uxyz, uaye, urho, utem, fh + real*8 dom, zr ! Density, energy and velocity fields fields @@ -93,7 +97,6 @@ subroutine calc_temp_cloudy_g(d, e, metal, temperature, integer i, j, k integer t, dj, dk - real*8 dom, zr ! row temporaries @@ -106,11 +109,6 @@ subroutine calc_temp_cloudy_g(d, e, metal, temperature, !\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\///////////////////////////////// !======================================================================= -! Set units - - dom = urho*(aye**3)/mh - zr = 1._DKIND/(aye*uaye) - 1._DKIND - ! We explicitly assume that the densities have already been converted ! from comoving to proper (if iexpand .eq. 1), before this function is ! called. (Historically that logic was handled here. diff --git a/src/clib/calculate_temperature.c b/src/clib/calculate_temperature.c index c2e582bc..673ca73d 100644 --- a/src/clib/calculate_temperature.c +++ b/src/clib/calculate_temperature.c @@ -40,7 +40,7 @@ extern void FORTRAN_NAME(calc_temp_cloudy_g)( gr_float *d, gr_float *e, gr_float *metal, gr_float *temperature, int *in, int *jn, int *kn, int *iexpand, int *imetal, int *is, int *js, int *ks, int *ie, int *je, int *ke, - double *aye, double *temstart, double *temend, + double *dom, double *zr, double *aye, double *temstart, double *temend, double *utem, double *uxyz, double *uaye, double *urho, double *utim, double *gamma, double *fh, long long *priGridRank, long long *priGridDim, @@ -218,6 +218,9 @@ int local_calculate_temperature_table(chemistry_data *my_chemistry, pow(my_units->a_value, -3.0)); } + double dom = co_density_units * pow(my_units->a_value, 3) / mh; + double zr = 1.0/(my_units->a_value * my_units->a_units) - 1.0; + FORTRAN_NAME(calc_temp_cloudy_g)( my_fields->density, my_fields->internal_energy, @@ -234,6 +237,8 @@ int local_calculate_temperature_table(chemistry_data *my_chemistry, my_fields->grid_end, my_fields->grid_end+1, my_fields->grid_end+2, + &dom, + &zr, &my_units->a_value, &my_chemistry->TemperatureStart, &my_chemistry->TemperatureEnd, From 0c33bcb3875cc5c3410797595781da364aaa7ae6 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 19 Jul 2024 23:39:04 -0400 Subject: [PATCH 5/5] remove extra args from calc_temp_cloudy_g --- src/clib/calc_temp_cloudy_g.F | 21 ++++----------------- src/clib/calculate_temperature.c | 10 ++-------- 2 files changed, 6 insertions(+), 25 deletions(-) diff --git a/src/clib/calc_temp_cloudy_g.F b/src/clib/calc_temp_cloudy_g.F index 004de677..f0857661 100644 --- a/src/clib/calc_temp_cloudy_g.F +++ b/src/clib/calc_temp_cloudy_g.F @@ -6,10 +6,8 @@ subroutine calc_temp_cloudy_g(d, e, metal, temperature, & in, jn, kn, iexpand, imetal, & is, js, ks, ie, je, ke, - & dom, zr, - & aye, temstart, temend, - & utem, uxyz, uaye, urho, utim, - & gamma, fh, + & dom, zr, temstart, temend, + & utem, gamma, fh, & priGridRank, priGridDim, & priPar1, priPar2, priPar3, & priDataSize, priMMW) @@ -38,12 +36,7 @@ subroutine calc_temp_cloudy_g(d, e, metal, temperature, ! dom - unit conversion to proper number density in code units ! zr - current redshift ! fh - Hydrogen mass fraction (typically 0.76) -! aye - expansion factor (in code units) ! -! utim - time units (i.e. code units to CGS conversion factor) -! uaye - expansion factor conversion factor (uaye = 1/(1+zinit)) -! urho - density units -! uxyz - length units ! utem - temperature(-like) units ! ! temstart, temend - start and end of temperature range for rate table @@ -72,9 +65,8 @@ subroutine calc_temp_cloudy_g(d, e, metal, temperature, integer in, jn, kn, is, js, ks, ie, je, ke, & iexpand, imetal - real*8 aye, temstart, temend, gamma, - & utim, uxyz, uaye, urho, utem, fh - real*8 dom, zr + real*8 temstart, temend, gamma, utem, fh, + & dom, zr ! Density, energy and velocity fields fields @@ -88,11 +80,6 @@ subroutine calc_temp_cloudy_g(d, e, metal, temperature, real*8 priPar1(priGridDim(1)), priPar2(priGridDim(2)), & priPar3(priGridDim(3)), priMMW(priDataSize) -! Parameters - - real*8 mh - parameter (mh = mass_h) - ! Locals integer i, j, k diff --git a/src/clib/calculate_temperature.c b/src/clib/calculate_temperature.c index 673ca73d..4841288a 100644 --- a/src/clib/calculate_temperature.c +++ b/src/clib/calculate_temperature.c @@ -40,9 +40,8 @@ extern void FORTRAN_NAME(calc_temp_cloudy_g)( gr_float *d, gr_float *e, gr_float *metal, gr_float *temperature, int *in, int *jn, int *kn, int *iexpand, int *imetal, int *is, int *js, int *ks, int *ie, int *je, int *ke, - double *dom, double *zr, double *aye, double *temstart, double *temend, - double *utem, double *uxyz, double *uaye, double *urho, double *utim, - double *gamma, double *fh, + double *dom, double *zr, double *temstart, double *temend, + double *utem, double *gamma, double *fh, long long *priGridRank, long long *priGridDim, double *priPar1, double *priPar2, double *priPar3, long long *priDataSize, double *priMMW); @@ -239,14 +238,9 @@ int local_calculate_temperature_table(chemistry_data *my_chemistry, my_fields->grid_end+2, &dom, &zr, - &my_units->a_value, &my_chemistry->TemperatureStart, &my_chemistry->TemperatureEnd, &temperature_units, - &co_length_units, - &my_units->a_units, - &co_density_units, - &my_units->time_units, &my_chemistry->Gamma, &my_chemistry->HydrogenFractionByMass, &my_rates->cloudy_primordial.grid_rank,