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 (1 d, 22 h)
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