Skip to content

Commit

Permalink
Merge pull request #106 from mabruzzo/grid_dim_lim
Browse files Browse the repository at this point in the history 
Modify property calculations to respect `grid_start` and `grid_stop`
  • Loading branch information
@brittonsmith
brittonsmith committed May 18, 2022
2 parents 121e980 + 0d5126a commit cf6dd0f
Show file tree
Hide file tree
Showing 9 changed files with 679 additions and 113 deletions.
1 change: 1 addition & 0 deletions src/clib/Make.config.objects
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@ OBJS_CONFIG_LIB = \
cool1d_multi_g.lo \
cool_multi_time_g.lo \
grackle_units.lo \
index_helper.lo \
initialize_chemistry_data.lo \
initialize_cloudy_data.lo \
initialize_UVbackground_data.lo \
Expand Down
76 changes: 46 additions & 30 deletions src/clib/calculate_gamma.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@
#include "grackle_types.h"
#include "grackle_chemistry_data.h"
#include "phys_constants.h"
#include "index_helper.h"
#ifdef _OPENMP
#include <omp.h>
#endif
Expand All @@ -31,6 +32,8 @@ int local_calculate_temperature(chemistry_data *my_chemistry,
grackle_field_data *my_fields,
gr_float *temperature);

grackle_index_helper _build_index_helper(const grackle_field_data *my_fields);

int local_calculate_gamma(chemistry_data *my_chemistry,
chemistry_data_storage *my_rates,
code_units *my_units,
Expand All @@ -41,15 +44,19 @@ int local_calculate_gamma(chemistry_data *my_chemistry,
if (!my_chemistry->use_grackle)
return SUCCESS;

int i, dim, size = 1;
for (dim = 0; dim < my_fields->grid_rank; dim++)
size *= my_fields->grid_dimension[dim];
const grackle_index_helper ind_helper = _build_index_helper(my_fields);
int outer_ind, index;

/* If molecular hydrogen is not being used, just use monotonic.
(this should not really be called, but provide it just in case). */

for (i = 0; i < size; i++) {
my_gamma[i] = my_chemistry->Gamma;

for (outer_ind = 0; outer_ind < ind_helper.outer_ind_size; outer_ind++){

const grackle_index_range range = _inner_range(outer_ind, &ind_helper);

for (index = range.start; index <= range.end; index++) {
my_gamma[index] = my_chemistry->Gamma;
}
}

if (my_chemistry->primordial_chemistry > 1) {
Expand All @@ -68,39 +75,48 @@ int local_calculate_gamma(chemistry_data *my_chemistry,
double x, nH2, number_density, GammaH2Inverse,
GammaInverse = 1 / (my_chemistry->Gamma - 1.0);

/* parallelize the k and j loops with OpenMP
* (these loops are flattened them for better parallelism) */
# ifdef _OPENMP
# pragma omp parallel for schedule( runtime ) \
private( x, nH2, number_density, GammaH2Inverse )
private( outer_ind, index, x, nH2, number_density, GammaH2Inverse )
# endif
for (i = 0; i < size; i++) {

/* Compute relative number abundence of molecular hydrogen. */

number_density =
0.25 * (my_fields->HeI_density[i] + my_fields->HeII_density[i] +
my_fields->HeIII_density[i]) +
my_fields->HI_density[i] + my_fields->HII_density[i] +
my_fields->HM_density[i] + my_fields->e_density[i];
for (outer_ind = 0; outer_ind < ind_helper.outer_ind_size; outer_ind++){

const grackle_index_range range = _inner_range(outer_ind, &ind_helper);

for (index = range.start; index <= range.end; index++) {

nH2 = 0.5 * (my_fields->H2I_density[i] + my_fields->H2II_density[i]);
/* Compute relative number abundence of molecular hydrogen. */

/* Only do full computation if there is a reasonable amount of H2.
The second term in GammaH2Inverse accounts for the vibrational
degrees of freedom. */
number_density =
0.25 * (my_fields->HeI_density[index] +
my_fields->HeII_density[index] +
my_fields->HeIII_density[index]) +
my_fields->HI_density[index] + my_fields->HII_density[index] +
my_fields->HM_density[index] + my_fields->e_density[index];

GammaH2Inverse = 0.5*5.0;
if (nH2 / number_density > 1e-3) {
x = 6100.0 / my_gamma[i];
if (x < 10.0)
GammaH2Inverse = 0.5*(5 + 2.0 * x*x * exp(x)/POW(exp(x)-1.0,2));
}
nH2 = 0.5 * (my_fields->H2I_density[index] +
my_fields->H2II_density[index]);

/* Add in H2. */
/* Only do full computation if there is a reasonable amount of H2.
The second term in GammaH2Inverse accounts for the vibrational
degrees of freedom. */

my_gamma[i] = 1.0 + (nH2 + number_density) /
(nH2 * GammaH2Inverse + number_density * GammaInverse);
GammaH2Inverse = 0.5*5.0;
if (nH2 / number_density > 1e-3) {
x = 6100.0 / my_gamma[index];
if (x < 10.0)
GammaH2Inverse = 0.5*(5 + 2.0 * x*x * exp(x)/POW(exp(x)-1.0,2));
}

/* Add in H2. */

my_gamma[index] = 1.0 + (nH2 + number_density) /
(nH2 * GammaH2Inverse + number_density * GammaInverse);

} // end: loop over i
} // end: loop over index
} // end: loop over outer_ind

} // end: if (my_chemistry->primordial_chemistry > 1)

Expand Down
109 changes: 63 additions & 46 deletions src/clib/calculate_pressure.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@
#include "grackle_types.h"
#include "grackle_chemistry_data.h"
#include "phys_constants.h"
#include "index_helper.h"
#ifdef _OPENMP
#include <omp.h>
#endif
Expand All @@ -38,24 +39,31 @@ int local_calculate_pressure(chemistry_data *my_chemistry,
return SUCCESS;

double tiny_number = 1.e-20;
int i, dim, size = 1;
for (dim = 0; dim < my_fields->grid_rank; dim++)
size *= my_fields->grid_dimension[dim];
const grackle_index_helper ind_helper = _build_index_helper(my_fields);
int outer_ind, index;

/* parallelize the k and j loops with OpenMP
* (these loops are flattened them for better parallelism) */
# ifdef _OPENMP
# pragma omp parallel for schedule( runtime ) private( i )
# pragma omp parallel for schedule( runtime ) private( outer_ind, index )
# endif
for (i = 0; i < size; i++) {

pressure[i] = (my_chemistry->Gamma - 1.0) * my_fields->density[i] *
my_fields->internal_energy[i];
for (outer_ind = 0; outer_ind < ind_helper.outer_ind_size; outer_ind++){

const grackle_index_range range = _inner_range(outer_ind, &ind_helper);

for (index = range.start; index <= range.end; index++) {

pressure[index] = ((my_chemistry->Gamma - 1.0) *
my_fields->density[index] *
my_fields->internal_energy[index]);

if (pressure[i] < tiny_number)
pressure[i] = tiny_number;
}

if (pressure[index] < tiny_number)
pressure[index] = tiny_number;
} // end: loop over i
} // end: loop over outer_ind

/* Correct for Gamma from H2. */

if (my_chemistry->primordial_chemistry > 1) {

/* Calculate temperature units. */
Expand All @@ -67,43 +75,52 @@ int local_calculate_pressure(chemistry_data *my_chemistry,

# ifdef _OPENMP
# pragma omp parallel for schedule( runtime ) \
private( i, number_density, nH2, GammaH2Inverse, x, Gamma1, temp )
private( outer_ind, index, \
number_density, nH2, GammaH2Inverse, x, Gamma1, temp )
# endif
for (i = 0; i < size; i++) {

number_density =
0.25 * (my_fields->HeI_density[i] + my_fields->HeII_density[i] +
my_fields->HeIII_density[i]) +
my_fields->HI_density[i] + my_fields->HII_density[i] +
my_fields->HM_density[i] + my_fields->e_density[i];

nH2 = 0.5 * (my_fields->H2I_density[i] + my_fields->H2II_density[i]);

/* First, approximate temperature. */

if (number_density == 0)
number_density = tiny_number;
temp = max(temperature_units * pressure[i] / (number_density + nH2), 1);

/* Only do full computation if there is a reasonable amount of H2.
The second term in GammaH2Inverse accounts for the vibrational
degrees of freedom. */

GammaH2Inverse = 0.5*5.0;
if (nH2 / number_density > 1e-3) {
x = 6100.0 / temp;
if (x < 10.0)
GammaH2Inverse = 0.5*(5 + 2.0 * x*x * exp(x)/POW(exp(x)-1.0,2));
}

Gamma1 = 1.0 + (nH2 + number_density) /
(nH2 * GammaH2Inverse + number_density * GammaInverse);
for (int outer_ind = 0; outer_ind < ind_helper.outer_ind_size; outer_ind++){

const grackle_index_range range = _inner_range(outer_ind, &ind_helper);

for (index = range.start; index <= range.end; index++) {

number_density =
0.25 * (my_fields->HeI_density[index] +
my_fields->HeII_density[index] +
my_fields->HeIII_density[index]) +
my_fields->HI_density[index] + my_fields->HII_density[index] +
my_fields->HM_density[index] + my_fields->e_density[index];

nH2 = 0.5 * (my_fields->H2I_density[index] +
my_fields->H2II_density[index]);

/* First, approximate temperature. */

if (number_density == 0)
number_density = tiny_number;
temp = max(temperature_units * pressure[index] / (number_density + nH2),
1);

/* Only do full computation if there is a reasonable amount of H2.
The second term in GammaH2Inverse accounts for the vibrational
degrees of freedom. */

GammaH2Inverse = 0.5*5.0;
if (nH2 / number_density > 1e-3) {
x = 6100.0 / temp;
if (x < 10.0)
GammaH2Inverse = 0.5*(5 + 2.0 * x*x * exp(x)/POW(exp(x)-1.0,2));
}

Gamma1 = 1.0 + (nH2 + number_density) /
(nH2 * GammaH2Inverse + number_density * GammaInverse);

/* Correct pressure with improved Gamma. */
/* Correct pressure with improved Gamma. */

pressure[i] *= (Gamma1 - 1.0) / (my_chemistry->Gamma - 1.0);
pressure[index] *= (Gamma1 - 1.0) / (my_chemistry->Gamma - 1.0);

} // end: loop over i
} // end: loop over i
} // end: loop over outer_ind

} // end: if (my_chemistry->primordial_chemistry > 1)

Expand Down
76 changes: 40 additions & 36 deletions src/clib/calculate_temperature.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@
#include "grackle_types.h"
#include "grackle_chemistry_data.h"
#include "phys_constants.h"
#include "index_helper.h"
#ifdef _OPENMP
#include <omp.h>
#endif
Expand Down Expand Up @@ -79,12 +80,6 @@ int local_calculate_temperature(chemistry_data *my_chemistry,
return FAIL;
}
}

/* Compute the size of the fields. */

int i, dim, size = 1;
for (dim = 0; dim < my_fields->grid_rank; dim++)
size *= my_fields->grid_dimension[dim];

/* Calculate temperature units. */

Expand All @@ -102,38 +97,53 @@ int local_calculate_temperature(chemistry_data *my_chemistry,
return SUCCESS;
}

/* Compute temperature with mu calculated directly. */

/* Compute properties used to index the field. */
const grackle_index_helper ind_helper = _build_index_helper(my_fields);
int outer_ind, index;

/* Compute temperature with mu calculated directly. */

/* parallelize the k and j loops with OpenMP
* (these loops are flattened them for better parallelism) */
# ifdef _OPENMP
# pragma omp parallel for schedule( runtime ) private( i, number_density )
# pragma omp parallel for schedule( runtime ) \
private( outer_ind, index, number_density )
# endif
for (i = 0; i < size; i++) {

if (my_chemistry->primordial_chemistry > 0) {
number_density =
0.25 * (my_fields->HeI_density[i] + my_fields->HeII_density[i] +
my_fields->HeIII_density[i]) +
my_fields->HI_density[i] + my_fields->HII_density[i] +
my_fields->e_density[i];
}
for (outer_ind = 0; outer_ind < ind_helper.outer_ind_size; outer_ind++){

const grackle_index_range range = _inner_range(outer_ind, &ind_helper);

/* Add in H2. */
for (index = range.start; index <= range.end; index++) {

if (my_chemistry->primordial_chemistry > 1) {
number_density += my_fields->HM_density[i] +
0.5 * (my_fields->H2I_density[i] + my_fields->H2II_density[i]);
}
if (my_chemistry->primordial_chemistry > 0) {
number_density =
0.25 * (my_fields->HeI_density[index] +
my_fields->HeII_density[index] +
my_fields->HeIII_density[index]) +
my_fields->HI_density[index] + my_fields->HII_density[index] +
my_fields->e_density[index];
}

if (my_fields->metal_density != NULL) {
number_density += my_fields->metal_density[i] * inv_metal_mol;
}
/* Add in H2. */

/* Ignore deuterium. */
if (my_chemistry->primordial_chemistry > 1) {
number_density += my_fields->HM_density[index] +
0.5 * (my_fields->H2I_density[index] +
my_fields->H2II_density[index]);
}

if (my_fields->metal_density != NULL) {
number_density += my_fields->metal_density[index] * inv_metal_mol;
}

temperature[i] *= temperature_units / max(number_density, tiny_number);
temperature[i] = max(temperature[i], MINIMUM_TEMPERATURE);
}
/* Ignore deuterium. */

temperature[index] *= temperature_units / max(number_density,
tiny_number);
temperature[index] = max(temperature[index], MINIMUM_TEMPERATURE);
} // end: loop over i
} // end: loop over outer_ind

return SUCCESS;
}

Expand All @@ -151,12 +161,6 @@ int local_calculate_temperature_table(chemistry_data *my_chemistry,
fprintf(stderr, "ERROR: this function requires primordial_chemistry set to 0.\n");
return FAIL;
}

/* Compute the size of the fields. */

int i, dim, size = 1;
for (dim = 0; dim < my_fields->grid_rank; dim++)
size *= my_fields->grid_dimension[dim];

/* Check for a metal field. */

Expand Down

0 comments on commit cf6dd0f

Please sign in to comment.