Page MenuHomec4science
Files F88762760

compute_lambda_interface.c
No OneTemporary
Actions

File Metadata

Created
Sun, Oct 20, 13:48

compute_lambda_interface.c
View Options

/************************************************************************************************
*
* Compute the cooling function Lambda from CLOUDY hdf5 output tables
* _/\_ = _/\_(z, rho_H, T)
*
* This file includes a Python interface for python2.7
*
* Jeremie Despraz, 2012
*
* Compilation:
* make
*
************************************************************************************************/
/************************************************************
* ATTENTION
*
* hdf5 table files need to be written in the following
* format z_%f.hdf5 (ex: z_0.615.hdf5).
*
* reference files can be downloaded on:
* http://www.strw.leidenuniv.nl/WSS08/
*
* references for the use of the hdf5 library can be
* found on:
* http://www.hdfgroup.org/HDF5/doc/H5.intro.html#Intro-APIs
* http://www.hdfgroup.org/HDF5/doc/UG/UG_frame08TheFile.html
*
************************************************************/
/*********************
*
* TODO:
* efficient error and
* exception handling
*
*********************/
#include <Python.h>
#include "proto.h"
#include <stdlib.h>
#include <math.h>
#include <mpi.h>
//standard C POSIX library for file handling
//(opendir, readdir, etc)
#include <dirent.h>
int withLinInterpolation = 1;
/*******************************************************************************/
int computeLambda(){
int err = 1;
int type = -1;
char* type_name = "";
float lambda = 0.0;
float z_in = 0.0;
float rho_H_in = 0.0;
float T_in = 1.0e4;
float nHe_in = 0.0;
//read user inputs
//determine the type of element for which lambda will be computed
printf("element type:\n "
"[1]Calcium\n [2]Carbon\n [3]Iron\n [4]Magnesium\n [5]Neon\n "
"[6]Nitrogen\n [7]Oxygen\n [8]Silicon\n [9]Sulphur\n "
"[10]Total_Metals\n [11]Solar\n [12]Metal_free\n");
err = scanf("%d", &type);
if(type < 1 || type > 12 || err != 1){
printf("\nerror: element type should be in the interval [1,12]\n\n");
return 1;
}
switch(type){
case 1:
type_name = "/Calcium";
break;
case 2:
type_name = "/Carbon";
break;
case 3:
type_name = "/Iron";
break;
case 4:
type_name = "/Magnesium";
break;
case 5:
type_name = "/Neon";
break;
case 6:
type_name = "/Nitrogen";
break;
case 7:
type_name = "/Oxygen";
break;
case 8:
type_name = "/Silicon";
break;
case 9:
type_name = "/Sulphur";
break;
case 10:
type_name = "/Total_Metals";
break;
case 11:
type_name = "/Solar";
break;
case 12:
type_name = "/Metal_free";
break;
}
//determine the redshift value
printf("redshift z:\n");
err = scanf("%f", &z_in);
if(z_in < 0.0 || z_in > 10.0 || err != 1){
printf("\nerror: redshift value should be in the interval [0,10]\n\n");
return 1;
}
//determine the hydrogen density
printf("H density rho_H:\n");
err = scanf("%f", &rho_H_in);
if(rho_H_in < 0.0 || rho_H_in > 1.0 || err != 1){
printf("\nerror: hydrogen density value should be in the interval [0,1]\n\n");
return 1;
}
//determine the temperature
printf("temperature T:\n");
err = scanf("%f", &T_in);
if(T_in < 0.0 || T_in > 10.0e10 || err != 1){
printf("\nerror: temperature value should be in the interval [0,1e10]\n\n");
return 1;
}
//only in the case of metal free gas
if(strcmp(type_name,"/Metal_free") == 0){
//detetrmine the Helium aboundance
printf("Helium aboundance nHe:\n");
err = scanf("%f", &nHe_in);
if(nHe_in < 0.0 || nHe_in > 1.0 || err != 1){
printf("\nerror: Helium aboundance value should be in the interval [0,1]\n\n");
return 1;
}
}
//compute the value closest possible value of lambda from the tables
lambda = readhdf5CoolingTable(type_name, z_in, rho_H_in, T_in, nHe_in);
printf("\nlambda(%g,%g,%g,%g) = %g\n", z_in, rho_H_in, T_in, nHe_in, lambda);
return 0;
}
/*******************************************************************************/
// ADD DESCRIPTION HERE
//variables ending with "_in" are user input, other variables are closest match
//in the hdf5 tables
float readhdf5CoolingTable(char* type_name, float z_in, float rho_H_in, float T_in, float nHe_in){
float lambda = 0.0;
//hdf5 tables location
char* tables_dir = "../coolingtables/";
printf("tables location:\t%s\n", tables_dir);
//load corresponding hdf5 table
int err = 0;
int hdf5 = 0;
float z_file = 0.0;
float diff = 0.0;
float min = 1.0e10;
DIR *dir;
struct dirent* content;
char* file_name = "";
//open directory containing the tables
dir = opendir(tables_dir);
if(dir == NULL){
printf("an error occured while opening hdf5 directory %s\n", tables_dir);
return 0;
}
/*********************** z ***********************/
float z = 0.0;
//scan directory content
while((content = readdir(dir)) != NULL){
//consider only files having the extension .hdf5
hdf5 = endsWith(content->d_name, ".hdf5");
if(hdf5 == 1){
//scan file name to determine the redshift value
err = sscanf(content->d_name, "z_%f.hdf5", &z_file);
if(err != 1){
printf("an error occured while reading file %s in directory %s\n", content->d_name, tables_dir);
return 0;
}
//determine the file closest to the redshift value given as input
diff = fabs(z_file - z_in);
if(diff < min){
file_name = content->d_name;
z = z_file;
min = diff;
}
}
}
//close directory
err = closedir(dir);
if (err != 0){
printf("an error occured while closing directory %s\n", tables_dir);
return 0;
}
//variables types specific to the hdf5 library
hid_t table;
//store path to table file
char file_path[(int)strlen(tables_dir)+(int)strlen(file_name)];
strcpy(file_path, tables_dir);
strcat(file_path, file_name);
//load hdf5 table
table = H5Fopen(file_path, H5F_ACC_RDONLY, H5P_DEFAULT);
//find the indices for temperature, hydrogen density and
//helium mass fraction that correspond to the closest
//input values
/*********************** T ***********************/
float T [2] = {0.0, 0.0};
int T_index[2] = {0,0};
char* T_key = "/Temperature_bins";
char T_table_key[(int)strlen(type_name)+(int)strlen(T_key)];
strcpy(T_table_key, type_name);
strcat(T_table_key, T_key);
closestMatch1D(table, T_table_key, T_in, T, T_index);
/********************* rho_H *********************/
float rho_H [2] = {0.0, 0.0};
int rho_H_index[2] = {0,0};
char* rho_H_key = "/Hydrogen_density_bins";
char rho_H_table_key[(int)strlen(type_name)+(int)strlen(rho_H_key)];
strcpy(rho_H_table_key, type_name);
strcat(rho_H_table_key, rho_H_key);
closestMatch1D(table, rho_H_table_key, rho_H_in, rho_H, rho_H_index);
/********************** nHe **********************/
float nHe [2] = {0.0, 0.0};
int nHe_index[2] = {-1,-1};
if(strcmp(type_name,"/Metal_free") == 0){
char* nHe_key = "/Helium_mass_fraction_bins";
//or
//char* nHe_key = "/Helium_number_ratio_bins";
char nHe_table_key[(int)strlen(type_name)+(int)strlen(nHe_key)];
strcpy(nHe_table_key, type_name);
strcat(nHe_table_key, nHe_key);
closestMatch1D(table, nHe_table_key, nHe_in, nHe, nHe_index);
}
/********************* lambda *********************/
//load corresponding lambda value
char* lambda_key = "/Net_Cooling"; // /!\ sometimes written "/Net_cooling"
char lambda_table_key[(int)strlen(type_name)+(int)strlen(lambda_key)];
strcpy(lambda_table_key, type_name);
strcat(lambda_table_key, lambda_key);
hid_t dataset;
hid_t dataspace;
hid_t memspace;
herr_t status;
hsize_t dims[3];
int rank;
dataset = H5Dopen2(table, lambda_table_key, H5P_DEFAULT); // open dataset (key)
dataspace = H5Dget_space(dataset); // get dataspace
rank = H5Sget_simple_extent_ndims(dataspace); // compute dataspace rank
status = H5Sget_simple_extent_dims(dataspace, dims, NULL);// save dataspace dimensions in dims
memspace = H5Screate_simple(rank, dims, NULL); // determine memory space required
//read lambda value from indices position
if(nHe_index[0] != -1){
float data[dims[0]][dims[1]][dims[2]];
status = H5Dread(dataset, H5T_NATIVE_FLOAT, memspace, dataspace, H5P_DEFAULT, data);
lambda = data[nHe_index[0]][T_index[0]][rho_H_index[0]];
if(withLinInterpolation == 1){
//linear interpolation of the results
//the derivative is computed by a simple finite difference method.
//f(c) ~ \-/f(a) * (c-a) ~ (f(b)-f(a))/(b-a) * (c-a)
lambda += (data[nHe_index[0]][T_index[0]][rho_H_index[0]]-data[nHe_index[1]][T_index[0]][rho_H_index[0]])
/ (nHe[0]-nHe[1]) * (nHe_in - nHe[0]);
lambda += (data[nHe_index[0]][T_index[0]][rho_H_index[0]]-data[nHe_index[0]][T_index[1]][rho_H_index[0]])
/ (T[0]-T[1]) * (T_in - T[0]);
lambda += (data[nHe_index[0]][T_index[0]][rho_H_index[0]]-data[nHe_index[0]][T_index[0]][rho_H_index[1]])
/ (rho_H[0]-rho_H[1]) * (rho_H_in - rho_H[0]);
}
} else {
float data[dims[0]][dims[1]];
status = H5Dread(dataset, H5T_NATIVE_FLOAT, memspace, dataspace, H5P_DEFAULT, data);
lambda = data[T_index[0]][rho_H_index[0]];
if(withLinInterpolation == 1){
//linear interpolation of the results
//the derivative is computed by a simple finite difference method.
//f(c) ~ \-/f(a) * (c-a) ~ (f(b)-f(a))/(b-a) * (c-a)
lambda += (data[T_index[0]][rho_H_index[0]]-data[T_index[1]][rho_H_index[0]])
/ (T[0]-T[1]) * (T_in - T[0]);
lambda += (data[T_index[0]][rho_H_index[0]]-data[T_index[0]][rho_H_index[1]])
/ (rho_H[0]-rho_H[1]) * (rho_H_in - rho_H[0]);
}
}
//close instances
H5Sclose(memspace);
H5Sclose(dataspace);
H5Dclose(dataset);
//close file
H5Fclose(table);
return lambda;
}
/*******************************************************************************/
//determine if the string str ends by the string suffix
//function used here to determine the extension of file given their name
int endsWith(const char *str, const char *suffix) {
if (!str || !suffix)
return 0;
size_t lenstr = strlen(str);
size_t lensuffix = strlen(suffix);
if (lensuffix > lenstr)
return 0;
return strncmp(str + lenstr - lensuffix, suffix, lensuffix) == 0;
}
/*******************************************************************************/
//return the indices of the 2 closest element of the hdf5 table matching the user
// input. Note that the dimension of the data has to be 1D
void closestMatch1D(hid_t table, char* table_key, float input, float match[2], int index[2]){
hid_t dataset;
hid_t dataspace;
hid_t memspace;
herr_t status;
hsize_t size[1];
int rank;
dataset = H5Dopen2(table, table_key, H5P_DEFAULT); // open dataset (key)
dataspace = H5Dget_space(dataset); // get dataspace
rank = H5Sget_simple_extent_ndims(dataspace); // compute dataspace rank
status = H5Sget_simple_extent_dims(dataspace, size, NULL);// save dataspace dimensions in dims
memspace = H5Screate_simple(rank, size, NULL); // determine memory space required
//read data
float data[size[0]];
status = H5Dread(dataset, H5T_NATIVE_FLOAT, memspace, dataspace, H5P_DEFAULT, data);
printf("%s: rank %d, dimension %d \n", table_key, rank, (int)size[0]);
//determine the index of the closest values (above and below)
int k = 0;
float min1 = 1e30;
float min2 = 1e30;
float diff = 0;
//the value having the smallest difference
//to the input is seeked and the index is stored
for(k=0; k<size[0]; k++){
diff = fabs(input - data[k]);
if(diff < min1){
match[0] = data[k];
index[0] = k;
min1 = diff;
}
else if(diff < min2){
match[1] = data[k];
index[1] = k;
min2 = diff;
}
}
printf("closest match#1 for %f is %f with index %d\n", input, match[0], index[0]);
printf("closest match#2 for %f is %f with index %d\n", input, match[1], index[1]);
//close instances
H5Sclose(memspace);
H5Sclose(dataspace);
H5Dclose(dataset);
}
/************************************* Python Interface ************************************/
/*****************************************************/
/* functions declarations */
/*****************************************************/
static PyObject* computeLambdaInterface(PyObject *self, PyObject *args) {
int error_status;
error_status = computeLambda();
return Py_BuildValue("d",error_status);
}
static PyObject* readhdf5CoolingTableInterface(PyObject *self, PyObject *args) {
const char* type_name;
float z_in;
float rho_H_in;
float T_in;
float nHe_in;
float lambda;
if (!PyArg_ParseTuple(args, "sffff", &type_name, &z_in, &rho_H_in, &T_in, &nHe_in)) {
return NULL;
}
printf("%s\n","computing function with parameters:");
printf("%s: %s\n","type", type_name);
printf("%s: %g\n","redshift", z_in);
printf("%s: %g\n","hydrogen density", rho_H_in);
printf("%s: %g\n","temperature", T_in);
printf("%s: %g\n","helium aboundance", nHe_in);
lambda = readhdf5CoolingTable(type_name, z_in, rho_H_in, T_in, nHe_in);
return Py_BuildValue("f",lambda);
}
/*****************************************************/
/* Method table */
/*****************************************************/
static PyMethodDef computeLambdaMethods[] = {
{"computeLambda", computeLambdaInterface, METH_VARARGS,
"Asks for inputs and compute the cooling accordingly."},
{"readhdf5CoolingTable", readhdf5CoolingTableInterface, METH_VARARGS,
"Compute the cooling from given inputs."},
{NULL, NULL, 0, NULL} /* Sentinel */
};
/*****************************************************/
/* initialization */
/*****************************************************/
//void initcompute_lambda_interface(void)
PyMODINIT_FUNC initcompute_lambda_interface(void) {
(void) Py_InitModule("compute_lambda_interface", computeLambdaMethods);
}

Event Timeline

c4science · Help