Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F85236840
grackle_wrapper.cxx
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Fri, Sep 27, 17:29
Size
10 KB
Mime Type
text/x-c
Expires
Sun, Sep 29, 17:29 (2 d)
Engine
blob
Format
Raw Data
Handle
21144829
Attached To
rGEAR Gear
grackle_wrapper.cxx
View Options
/***********************************************************************
/
/ Grackle c wrapper
/
/
/ Copyright (c) 2013, Enzo/Grackle Development Team.
/
/ Distributed under the terms of the Enzo Public Licence.
/
/ The full license is in the file LICENSE, distributed with this
/ software.
************************************************************************/
#include <grackle.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "grackle_wrapper.h"
gr_float a_value;
gr_float temperature_units ;
gr_float velocity_units ;
gr_float cool_units ;
#define griddim 1
#define kboltz 1.3806504e-16
#define mass_h 1.67262171e-24
#define mass_e 9.10938215e-28
#define pi_val 3.14159265
#define hplanck 6.6260693e-27
#define ev2erg 1.60217653e-12
#define c_light 2.99792458e10
#define GravConst 6.6726e-8
#define sigma_sb 5.670373e-5
#define SolarMass 1.9891e33
#define Mpc 3.0857e24
#define kpc 3.0857e21
#define pc 3.0857e18
#define yr_to_s 3.15569e7
// arrays passed to grackle as input and to be filled
#define field_size 1
//gr_float density[field_size];
gr_float internal_energy[field_size];
gr_float x_velocity[field_size];
gr_float y_velocity[field_size];
gr_float z_velocity[field_size];
// for primordial_chemistry >= 1
gr_float HI_density[field_size];
gr_float HII_density[field_size];
gr_float HeI_density[field_size];
gr_float HeII_density[field_size];
gr_float HeIII_density[field_size];
gr_float ie_density[field_size];
// for primordial_chemistry >= 2
gr_float HM_density[field_size];
gr_float H2I_density[field_size];
gr_float H2II_density[field_size];
// for primordial_chemistry >= 3
gr_float DI_density[field_size];
gr_float DII_density[field_size];
gr_float HDI_density[field_size];
// for metal_cooling = 1
gr_float metal_density[field_size];
// calculated
gr_float cooling_time[field_size];
gr_float outgamma[field_size];
gr_float Density[field_size];
gr_float energy[field_size];
gr_float e_density[field_size];
// Set grid dimension and size.
// grid_start and grid_end are used to ignore ghost zones.
gr_int grid_rank = griddim;
// If grid rank is less than 3, set the other dimensions,
// start indices, and end indices to 0.
gr_int grid_dimension[griddim], grid_start[griddim], grid_end[griddim];
static chemistry_data my_chemistry;
code_units my_units;
int wrap_init_cooling(char* CloudyTable, int UVbackground, double udensity, double ulength, double utime, int grackle_chemistry){
/*********************************************************************
/ Initial setup of units and chemistry objects.
/ This should be done at simulation start.
*********************************************************************/
int i;
// First, set up the units system.
// These are conversions from code units to cgs.
my_units.a_units = 1.0; // units for the expansion factor (1/1+zi)
my_units.comoving_coordinates = 0; /* so, according to the doc, we assume here all physical quantities to be in proper coordiname (not comobile) */
my_units.density_units = udensity;
my_units.length_units = ulength;
my_units.time_units = utime;
velocity_units = my_units.a_units * my_units.length_units /my_units.time_units;
my_units.velocity_units = velocity_units;
//temperature_units = mass_h * velocity_units*velocity_units / kboltz;
// Second, create a chemistry object for parameters and rate data.
my_chemistry = set_default_chemistry_parameters();
// Set parameter values for chemistry.
my_chemistry.use_grackle = 1;
my_chemistry.with_radiative_cooling = 1;
my_chemistry.primordial_chemistry = grackle_chemistry; // molecular network with H, He, D
my_chemistry.metal_cooling = 1; // metal cooling on
my_chemistry.UVbackground = UVbackground;
my_chemistry.grackle_data_file = CloudyTable;
//my_chemistry.cmb_temperature_floor = 0;
//my_chemistry.LWbackground_sawtooth_suppression = 1;
// Finally, initialize the chemistry object.
// snl: a_value is not the true initial ae!! This should get set during update_UVbackground
gr_float initial_redshift = 0.;
a_value = 1. / (1. + initial_redshift);
if (initialize_chemistry_data(my_chemistry, my_units, a_value) == 0) {
return 0;
}
for (i = 0;i < griddim;i++) {
grid_dimension[i] = 0;
grid_start[i] = 0;
grid_end[i] = 0;
}
grid_dimension[0] = field_size;
grid_end[0] = field_size - 1;
return 1;
}
int wrap_update_UVbackground_rates(double auni) {
// The UV background rates must be updated before
// calling the other functions.
/* a_value = auni / my_units.a_units; */
a_value = auni / my_units.a_units;
if (update_UVbackground_rates(my_chemistry,
my_units, a_value) == 0) {
return 0;
}
return 1;
}
int wrap_get_cooling_time(double rho, double u, double vx, double vy, double vz,
double HI, double HII, double HeI, double HeII, double HeIII,
double HM, double H2I, double H2II, double DI, double DII, double HDI,
double ne, double Z, double a_now, double *coolingtime)
{
int i;
gr_float den_factor, u_factor;
if(field_size != 1){
fprintf(stderr,"field_size must currently be set to 1.\n");
return 0;
}
// passed density and energy are proper
/*
if(my_units.comoving_coordinates){
den_factor = pow(a_now, 3);
u_factor = pow(a_now, 0);
} else {
den_factor = 1.0;
u_factor = 1.0;
}
*/
den_factor = 1.0;
u_factor = 1.0;
for (i = 0;i < field_size;i++) {
Density[i] = rho * den_factor;
HI_density[i] = HI * Density[i];
HII_density[i] = HII * Density[i];
HeI_density[i] = HeI * Density[i];
HeII_density[i] = HeII * Density[i];
HeIII_density[i] = HeIII * Density[i];
HM_density[i] = HM * Density[i];
H2I_density[i] = H2I * Density[i];
H2II_density[i] = H2II * Density[i];
DI_density[i] = DI * Density[i];
DII_density[i] = DII * Density[i];
HDI_density[i] = HDI * Density[i];
e_density[i] = ne * Density[i];
metal_density[i] = Z * Density[i];
x_velocity[i] = vx;
y_velocity[i] = vy;
z_velocity[i] = vz;
energy[i] = u* u_factor;
}
if (my_chemistry.primordial_chemistry)
{
if (calculate_cooling_time(my_chemistry, my_units,
a_now,
grid_rank, grid_dimension,
grid_start, grid_end,
Density, energy,
x_velocity, y_velocity, z_velocity,
HI_density, HII_density, HM_density,
HeI_density, HeII_density, HeIII_density,
H2I_density, H2II_density,
DI_density, DII_density, HDI_density,
e_density, metal_density,
cooling_time) == 0) {
fprintf(stderr,"Error in calculate_cooling_time.\n");
return 0;
}
}
else
{
if (calculate_cooling_time(my_chemistry, my_units,
a_now,
grid_rank, grid_dimension,
grid_start, grid_end,
Density, energy,
x_velocity, y_velocity, z_velocity,
metal_density,
cooling_time) == 0) {
fprintf(stderr,"Error in calculate_cooling_time.\n");
return 0;
}
}
// return updated chemistry and energy
for (i = 0;i < field_size;i++)
*coolingtime = cooling_time[i];
return 1;
}
int wrap_do_cooling(double rho, double *u, double dt, double vx, double vy, double vz,
double *HI, double *HII, double *HeI, double *HeII, double *HeIII,
double *HM, double *H2I, double *H2II, double *DI, double *DII, double *HDI,
double *ne, double Z, double a_now)
{
int i;
gr_float den_factor, u_factor;
if(field_size != 1){
fprintf(stderr,"field_size must currently be set to 1.\n");
return 0;
}
// passed density and energy are proper
/*
if(my_units.comoving_coordinates){
den_factor = pow(a_now, 3);
u_factor = pow(a_now, 0);
} else {
den_factor = 1.0;
u_factor = 1.0;
}
*/
den_factor = 1.0;
u_factor = 1.0;
for (i = 0; i < field_size; i++) {
Density[i] = rho * den_factor;
HI_density[i] = *HI * Density[i];
HII_density[i] = *HII * Density[i];
HeI_density[i] = *HeI * Density[i];
HeII_density[i] = *HeII * Density[i];
HeIII_density[i] = *HeIII * Density[i];
HM_density[i] = *HM * Density[i];
H2I_density[i] = *H2I * Density[i];
H2II_density[i] = *H2II * Density[i];
DI_density[i] = *DI * Density[i];
DII_density[i] = *DII * Density[i];
HDI_density[i] = *HDI * Density[i];
e_density[i] = *ne * Density[i];
metal_density[i] = Z * Density[i];
x_velocity[i] = vx;
y_velocity[i] = vy;
z_velocity[i] = vz;
energy[i] = *u * u_factor;
}
if (my_chemistry.primordial_chemistry)
{
if (solve_chemistry(my_chemistry, my_units,
a_now, dt,
grid_rank, grid_dimension,
grid_start, grid_end,
Density, energy,
x_velocity, y_velocity, z_velocity,
HI_density, HII_density, HM_density,
HeI_density, HeII_density, HeIII_density,
H2I_density, H2II_density,
DI_density, DII_density, HDI_density,
e_density, metal_density) == 0) {
fprintf(stderr,"Error in solve_chemistry.\n");
return 0;
}
}
else
{
if (solve_chemistry(my_chemistry, my_units,
a_now, dt,
grid_rank, grid_dimension,
grid_start, grid_end,
Density, energy,
x_velocity, y_velocity, z_velocity,
metal_density) == 0) {
fprintf(stderr,"Error in solve_chemistry.\n");
return 0;
}
}
// return updated chemistry and energy
for (i = 0;i < field_size;i++) {
*HI = HI_density[i]/Density[i];
*HII = HII_density[i]/Density[i];
*HeI = HeI_density[i]/Density[i];
*HeII = HeII_density[i]/Density[i];
*HeIII = HeIII_density[i]/Density[i];
*HM = HM_density[i]/Density[i];
*H2I = H2I_density[i]/Density[i];
*H2II = H2II_density[i]/Density[i];
*DI = DI_density[i]/Density[i];
*DII = DII_density[i]/Density[i];
*HDI = HDI_density[i]/Density[i];
*ne = e_density[i]/Density[i];
*u = energy[i]/u_factor;
}
return 1;
}
Event Timeline
Log In to Comment