Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Transcribe the scale_fields_table_g function from Fortran to C. #228

Open
wants to merge 6 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
97 changes: 9 additions & 88 deletions src/clib/calc_temp_cloudy_g.F
Original file line number Diff line number Diff line change
Expand Up @@ -6,9 +6,8 @@
subroutine calc_temp_cloudy_g(d, e, metal, temperature,
& in, jn, kn, iexpand, imetal,
& is, js, ks, ie, je, ke,
& aye, temstart, temend,
& utem, uxyz, uaye, urho, utim,
& gamma, fh,
& dom, zr, temstart, temend,
& utem, gamma, fh,
& priGridRank, priGridDim,
& priPar1, priPar2, priPar3,
& priDataSize, priMMW)
Expand All @@ -34,13 +33,10 @@ 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)
!
! 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
Expand Down Expand Up @@ -69,8 +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 temstart, temend, gamma, utem, fh,
& dom, zr

! Density, energy and velocity fields fields

Expand All @@ -84,17 +80,10 @@ 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
integer t, dj, dk
real*8 dom, zr
real*8 dbase1, tbase1, xbase1

! row temporaries

Expand All @@ -107,23 +96,10 @@ 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

! 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

Expand Down Expand Up @@ -173,60 +149,5 @@ 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

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
56 changes: 48 additions & 8 deletions src/clib/calculate_temperature.c
Original file line number Diff line number Diff line change
Expand Up @@ -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 *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);
Expand Down Expand Up @@ -147,6 +146,34 @@ 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)
{
/* 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,
chemistry_data_storage *my_rates,
code_units *my_units,
Expand Down Expand Up @@ -184,6 +211,15 @@ 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));
}

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,
Expand All @@ -200,14 +236,11 @@ int local_calculate_temperature_table(chemistry_data *my_chemistry,
my_fields->grid_end,
my_fields->grid_end+1,
my_fields->grid_end+2,
&my_units->a_value,
&dom,
&zr,
&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,
Expand All @@ -218,6 +251,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;
}

Expand Down