Page MenuHomec4science
Files F85234371

cooling_grackle.c
No OneTemporary
Actions

File Metadata

Created
Fri, Sep 27, 17:03

cooling_grackle.c
View Options

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>
#include "allvars.h"
#include "proto.h"
#ifdef COOLING
#ifdef COOLING_GRACKLE
#include "grackle_wrapper.h"
#define mass_h 1.67262171e-24
#define Zsol 0.02041;
#define kboltz 1.3806504e-16
#define GRACKLE_CHEMISTRY 0
/*! This function initialize the cooling parameters for Grackle
*/
void init_cooling_GRACKLE()
{
char * cloudytable;
double units_density, units_length, units_time;
int grackle_chemistry;
int UVbackground;
cloudytable = All.GrackleCloudyTable;
UVbackground = All.GrackleUVbackground;
units_density = All.UnitDensity_in_cgs * All.HubbleParam * All.HubbleParam;
units_length = All.UnitLength_in_cm / All.HubbleParam;
units_time = All.UnitTime_in_s / All.HubbleParam;
grackle_chemistry = GRACKLE_CHEMISTRY;
//printf("Initializing cooling...\n");
if(wrap_init_cooling(cloudytable,UVbackground,units_density, units_length, units_time, grackle_chemistry) != 1)
{
fprintf(stderr,"Error in initialize_chemistry_data.");
exit(-1);
}
//printf("Cooling Initialization done.\n\n");
}
/*! This function computes the cooling time.
*/
#if (GRACKLE_CHEMISTRY == 0)
double GetCoolingTime_GRACKLE(double energy, double density, double ne, double Z, double a_now)
#endif
#if (GRACKLE_CHEMISTRY == 1)
double GetCoolingTime_GRACKLE(double energy, double density, double ne, double HI, double HII, double HeI, double HeII, double HeIII, double Z, double a_now)
#endif
#if (GRACKLE_CHEMISTRY == 2)
double GetCoolingTime_GRACKLE(double energy, double density, double ne, double HI, double HII, double HeI, double HeII, double HeIII, double HM, double H2I, double H2II, double Z, double a_now)
#endif
#if (GRACKLE_CHEMISTRY == 3)
double GetCoolingTime_GRACKLE(double energy, double density, double ne, double HI, double HII, double HeI, double HeII, double HeIII, double HM, double H2I, double H2II, double DI, double DII, double HDI, double Z, double a_now)
#endif
{
double x_velocity, y_velocity, z_velocity;
double cooling_time;
#if (GRACKLE_CHEMISTRY < 3)
double DI=0, DII=0, HDI=0;
#endif
#if (GRACKLE_CHEMISTRY < 2)
double HM=0, H2I=0, H2II=0;
#endif
#if (GRACKLE_CHEMISTRY < 1)
double HI=0, HII=0, HeI=0, HeII = 0, HeIII=0;
#endif
x_velocity = 0.0;
y_velocity = 0.0;
z_velocity = 0.0;
/*********************************************************************
call to the main chemistry solver
*********************************************************************/
if (wrap_get_cooling_time(density, energy,
x_velocity, y_velocity, z_velocity,
HI, HII, HeI, HeII, HeIII,
HM, H2I, H2II, DI, DII, HDI,
ne, Z, a_now, &cooling_time) == 0) {
fprintf(stderr, "Error in calculate_cooling_time.\n");
}
return cooling_time;
}
/*! This function computes the new entropy due to the cooling,
* between step t0 and t1.
*/
#if (GRACKLE_CHEMISTRY == 0)
double DoCooling_GRACKLE(double energy, double density, double dtime, double *ne, double Z, double a_now)
#endif
#if (GRACKLE_CHEMISTRY == 1)
double DoCooling_GRACKLE(double energy, double density, double dtime, double *ne, double *HI, double *HII, double *HeI, double *HeII, double *HeIII, double Z, double a_now)
#endif
#if (GRACKLE_CHEMISTRY == 2)
double DoCooling_GRACKLE(double energy, double density, double dtime, double *ne, double *HI, double *HII, double *HeI, double *HeII, double *HeIII, double *HM, double *H2I, double *H2II, double Z, double a_now)
#endif
#if (GRACKLE_CHEMISTRY == 3)
double DoCooling_GRACKLE(double energy, double density, double dtime, double *ne, double *HI, double *HII, double *HeI, double *HeII, double *HeIII, double *HM, double *H2I, double *H2II, double *DI, double *DII, double *HDI, double Z, double a_now)
#endif
{
double x_velocity, y_velocity, z_velocity;
double HI_g, HII_g, HeI_g, HeII_g, HeIII_g, ne_g;
double HM_g, H2I_g, H2II_g, DI_g, DII_g, HDI_g;
x_velocity = 0.0;
y_velocity = 0.0;
z_velocity = 0.0;
#if (GRACKLE_CHEMISTRY < 3)
DI_g = 0;
DII_g = 0;
HDI_g = 0;
#else
DI_g = *DI;
DII_g = *DII;
HDI_g = *HDI;
#endif
#if (GRACKLE_CHEMISTRY < 2)
HM_g = 0;
H2I_g = 0;
H2II_g = 0;
#else
HM_g = *HM;
H2I_g = *H2I;
H2II_g = *H2II;
#endif
#if (GRACKLE_CHEMISTRY < 1)
HI_g = 0;
HII_g = 0;
HeI_g = 0;
HeII_g = 0;
HeIII_g = 0;
ne_g = 0;
#else
HI_g = *HI;
HII_g = *HII;
HeI_g = *HeI;
HeII_g = *HeII;
HeIII_g = *HeIII;
ne_g = *ne;
#endif
/*********************************************************************
call to the main chemistry solver
*********************************************************************/
if (wrap_do_cooling(density, &energy, dtime,
x_velocity, y_velocity, z_velocity,
&HI_g, &HII_g, &HeI_g, &HeII_g, &HeIII_g,
&HM_g, &H2I_g, &H2II_g, &DI_g, &DII_g, &HDI_g,
&ne_g, Z, a_now) == 0) {
fprintf(stderr, "Error in do_cooling.\n");
return 0;
}
return energy;
}
/*! This function computes the new entropy due to the cooling,
* between step t0 and t1.
*/
void CoolingForOne_GRACKLE(int i, double dt_entr, double dt_entr2, double dt_entr3, double a_now, double a3inv)
{
/*
Ti_begstep Ti_endstep Ti_endstep + ti_step
n n+1/2 n+1 n+1+1/2 n+2
| | | | |
x ----------------------------------------------------------------------------------------
| | | | |
| |
v ---------------------------------------------
| |
tstart tend
(Ti_begstep+Ti_endstep)/2 Ti_endstep + ti_step/2
we to the cooling in 4 steps
effective cooling
1) adiabatic cooling : tstart -> tend
2) radiative cooling : tstart -> tend
predicted cooling
3) adiabatic cooling : tend -> Ti_endstep + ti_step
4) radiative cooling : tend -> Ti_endstep + ti_step
dt_entr = tend - tstart n+1/2 -> n+1+1/2
dt_entr2 = tend - Ti_endstep n+1+1/2 -> n+1
dt_entr3 = ti_step / 2 n+1+1/2 -> n+2
*/
double density,egyspec,ne,Z;
double Dg;
double egyspec_new,egyspec_init,egyspec_old;
double minentropy;
double DtEgySpec,DtEgySpec_rad,DtEgySpec_ad,EgySpec;
#ifdef ENTROPYPRED
double EgySpecPred;
#endif
if(All.ComovingIntegrationOn)
{
printf("dtime must be computed correctly...\n");
endrun(7765789);
}
density = SphP[i].Density*a3inv ; /* in proper */
#ifdef CHIMIE
Z = SphP[i].Metal[METALS];
#else
Z = All.GasMetal; //0.02041;
#endif
ne = 0.0; /* mass fraction of eletron */ /* useless for GRACKLE_CHEMISTRY = 0 */
/* compute the initial energy */
#ifdef DENSITY_INDEPENDENT_SPH
Dg = pow(SphP[i].EgyWtDensity * a3inv, GAMMA_MINUS1);
#else
Dg = pow(SphP[i].Density * a3inv, GAMMA_MINUS1);
#endif
#ifdef MULTIPHASE
if (SphP[i].Phase== GAS_SPH)
{
egyspec_init = SphP[i].Entropy / (GAMMA_MINUS1) * Dg;
DtEgySpec_rad = SphP[i].DtEntropy / (GAMMA_MINUS1) * Dg; /* check that... */
}
else
{
egyspec_init = SphP[i].Entropy;
DtEgySpec_rad = SphP[i].DtEntropy;
}
#else
egyspec_init = SphP[i].Entropy / (GAMMA_MINUS1) * Dg;
DtEgySpec_rad = SphP[i].DtEntropy / (GAMMA_MINUS1) * Dg;
#endif
/*******************************************************
1) do the adiabatic cooling/heating n+1/2 -> n+1+1/2
********************************************************/
/* In case of cooling, we prevent that the energy (and hence temperature decreases by more than a factor 0.5 */
if(DtEgySpec_ad * dt_entr > -0.5 * egyspec_init)
egyspec = egyspec_init + DtEgySpec_ad * dt_entr;
else
{
egyspec *= 0.5;
DtEgySpec_ad = -0.5*egyspec/dt_entr;
}
/* keep minimum allowed value */
if(egyspec < All.MinEgySpec)
{
egyspec = All.MinEgySpec;
DtEgySpec_ad = 0; /* this is needed for point 3) */
}
/*******************************************************
2) do the radiative cooling/heating n+1/2 -> n+1+1/2
********************************************************/
egyspec_old = egyspec;
/* now, compute the new energy */
/* here, egyspec, density, dtime are in physical units */
egyspec_new = DoCooling_GRACKLE(egyspec, density, dt_entr, &ne, Z, a_now);
/* compute the rate of change of energy */
DtEgySpec_rad = (egyspec_new-egyspec)/dt_entr;
/* In case of cooling, we prevent that the energy (and hence temperature decreases by more than a factor 0.5 */
if(DtEgySpec_rad * dt_entr > -0.5 * egyspec)
egyspec = egyspec + DtEgySpec_rad * dt_entr;
else
{
egyspec *= 0.5;
DtEgySpec_rad = -0.5*egyspec/dt_entr;
}
/* keep minimum allowed value */
if(egyspec < All.MinEgySpec)
{
egyspec = All.MinEgySpec;
DtEgySpec_rad = 0; /* this is needed for point 4) */
}
/* count energy lost by radiative cooling/heating */
LocalSysState.RadiatedEnergy += (egyspec-egyspec_old) * P[i].Mass;
/* EgySpec is now the energy at n+1+1/2 */
EgySpec = egyspec;
DtEgySpec = DtEgySpec_ad +DtEgySpec_rad;
/*******************************************************
now, we can predic energy at n+1, from n+1+1/2
********************************************************/
#ifdef ENTROPYPRED
EgySpecPred = EgySpec - dt_entr2*DtEgySpec;
#endif
/*******************************************************
3) check the adiabatic cooling/heating validity for n+1+1/2 -> n+2
********************************************************/
/* In case the timestep increases in the new step, we
make sure that we do not 'overcool' when deriving
predicted temperatures. The maximum timespan over
which prediction can occur is ti_step/2, i.e. from
the middle to the end of the current step */
if(egyspec + DtEgySpec_ad * dt_entr3 < 0.5 * egyspec)
{
egyspec *= 0.5;
DtEgySpec_ad = -0.5 * egyspec / dt_entr3;
}
/*******************************************************
4) check the adiabatic cooling/heating validity for n+1+1/2 -> n+2
********************************************************/
/* In case the timestep increases in the new step, we
make sure that we do not 'overcool' when deriving
predicted temperatures. The maximum timespan over
which prediction can occur is ti_step/2, i.e. from
the middle to the end of the current step */
if(egyspec + DtEgySpec_rad * dt_entr3 < 0.5 * egyspec)
{
egyspec *= 0.5;
DtEgySpec_rad = -0.5 * egyspec / dt_entr3;
}
/*******************************************************
transform into entropy
********************************************************/
DtEgySpec = DtEgySpec_ad + DtEgySpec_rad; /* now the value used in predict */
#ifdef MULTIPHASE
if (SphP[i].Phase== GAS_SPH)
{
SphP[i].Entropy = GAMMA_MINUS1 * EgySpec / Dg;
SphP[i].DtEntropy = GAMMA_MINUS1 * DtEgySpec / Dg;
#ifdef ENTROPYPRED
SphP[i].EntropyPred = GAMMA_MINUS1 * EgySpecPred / Dg;
#endif
}
else
{
SphP[i].Entropy = EgySpec;
SphP[i].DtEntropy = DtEgySpec;
SphP[i].EntropyPred = EgySpecPred;
}
#else
SphP[i].Entropy = GAMMA_MINUS1 * EgySpec / Dg;
SphP[i].DtEntropy = GAMMA_MINUS1 * DtEgySpec / Dg;
#ifdef ENTROPYPRED
SphP[i].EntropyPred = GAMMA_MINUS1 * EgySpecPred / Dg;
#endif
#endif /* MULTIPHASE */
}
/****************************************************************************************/
/*
/*
/*
/* PYTHON INTERFACE
/*
/*
/*
/****************************************************************************************/
#ifdef PY_INTERFACE
#include <Python.h>
#include <numpy/arrayobject.h>
static PyObject * grackle_InitDefaultParameters(void)
{
/* grackle parameters */
strcpy(All.GrackleCloudyTable,"table.hdf5");
All.GrackleUVbackground = 1.;
/* System of units */
All.UnitLength_in_cm = 3.085e+21; /* 1.0 kpc */
All.UnitMass_in_g = 1.989e+43; /* 1.0e10 solar masses */
All.UnitVelocity_in_cm_per_s = 20725573.785998672; /* 207 km/sec */
All.HubbleParam = 1;
return Py_BuildValue("i",1);
}
static PyObject * SetParameters(PyObject *dict)
{
PyObject *key;
PyObject *value;
int ivalue;
float fvalue;
double dvalue;
/* check that it is a PyDictObject */
if(!PyDict_Check(dict))
{
PyErr_SetString(PyExc_AttributeError, "argument is not a dictionary.");
return NULL;
}
if (PyDict_Size(dict)==0)
return Py_BuildValue("i",0);
Py_ssize_t pos=0;
while(PyDict_Next(dict,&pos,&key,&value))
{
if(PyString_Check(key))
{
/* grackle parameters */
if(strcmp(PyString_AsString(key), "GrackleCloudyTable")==0)
{
if(PyString_Check(value))
strcpy(All.GrackleCloudyTable,PyString_AsString(value));
}
if(strcmp(PyString_AsString(key), "GrackleUVbackground")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.GrackleUVbackground = PyInt_AsLong(value);
}
/* System of units */
if(strcmp(PyString_AsString(key), "UnitLength_in_cm")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.UnitLength_in_cm = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "UnitMass_in_g")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.UnitMass_in_g = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "UnitVelocity_in_cm_per_s")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.UnitVelocity_in_cm_per_s = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "HubbleParam")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.HubbleParam = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "GravityConstantInternal")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.GravityConstantInternal = PyFloat_AsDouble(value);
}
}
}
return Py_BuildValue("i",1);
}
static PyObject * grackle_SetParameters(PyObject *self, PyObject *args)
{
PyObject *dict;
/* here, we can have either arguments or dict directly */
if(PyDict_Check(args))
{
dict = args;
}
else
{
if (! PyArg_ParseTuple(args, "O",&dict))
return NULL;
}
SetParameters(dict);
return Py_BuildValue("i",1);
}
static PyObject * grackle_GetParameters(void)
{
PyObject *dict;
PyObject *key;
PyObject *value;
dict = PyDict_New();
/* grackle parameters */
key = PyString_FromString("GrackleCloudyTable");
value = PyString_FromString(All.GrackleCloudyTable);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("GrackleUVbackground");
value = PyInt_FromLong(All.GrackleUVbackground);
PyDict_SetItem(dict,key,value);
/* System of units */
key = PyString_FromString("UnitLength_in_cm");
value = PyFloat_FromDouble(All.UnitLength_in_cm);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("UnitMass_in_g");
value = PyFloat_FromDouble(All.UnitMass_in_g);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("UnitVelocity_in_cm_per_s");
value = PyFloat_FromDouble(All.UnitVelocity_in_cm_per_s);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("HubbleParam");
value = PyFloat_FromDouble(All.HubbleParam);
PyDict_SetItem(dict,key,value);
return Py_BuildValue("O",dict);
}
static PyObject *
grackle_init_grackle(PyObject *self, PyObject *args, PyObject *kwds)
{
PyObject *paramsDict=NULL;
static char *kwlist[] = {"params", NULL};
/* this fails with python2.6, I do not know why ??? */
if (! PyArg_ParseTupleAndKeywords(args, kwds, "|O",kwlist,&paramsDict))
{
PyErr_SetString(PyExc_ValueError,"init_grackle, error in parsing arguments.");
return NULL;
}
/* use default parameters */
grackle_InitDefaultParameters();
/* check that it is a PyDictObject */
if(!PyDict_Check(paramsDict))
{
PyErr_SetString(PyExc_AttributeError, "argument is not a dictionary.");
return NULL;
}
else
{
SetParameters(paramsDict);
}
All.UnitTime_in_s = All.UnitLength_in_cm / All.UnitVelocity_in_cm_per_s;
All.UnitDensity_in_cgs = All.UnitMass_in_g / pow(All.UnitLength_in_cm, 3);
init_cooling_GRACKLE();
return Py_BuildValue("O",Py_None);
}
static PyObject *
grackle_GetCoolingTime(self, args)
PyObject *self;
PyObject *args;
{
double energy, density, ne, Z, a_now;
double cooling_time;
if (!PyArg_ParseTuple(args, "ddddd", &energy, &density, &ne, &Z, &a_now))
return NULL;
cooling_time = GetCoolingTime_GRACKLE(energy, density, ne, Z, a_now);
return Py_BuildValue("d",cooling_time);
}
static PyMethodDef grackleMethods[] = {
{"InitDefaultParameters", (PyCFunction)grackle_InitDefaultParameters, METH_VARARGS,
"Init default parameters"},
{"SetParameters", (PyCFunction)grackle_SetParameters, METH_VARARGS,
"Set gadget parameters"},
{"GetParameters", (PyCFunction)grackle_GetParameters, METH_VARARGS,
"get some gadget parameters"},
{"init_grackle", grackle_init_grackle, METH_VARARGS| METH_KEYWORDS,
"Init grackle."},
{"GetCoolingTime", (PyCFunction)grackle_GetCoolingTime, METH_VARARGS,
"get cooling time"},
{NULL, NULL, 0, NULL} /* Sentinel */
};
void initgrackle(void)
{
(void) Py_InitModule("grackle", grackleMethods);
import_array();
}
#endif /* PY_INTERFACE */
#endif /* COOLING_GRACKLE */
#endif /* COOLING */

Event Timeline

c4science · Help