Homec4science
Diffusion Gear (master) Browse

Gear/src/lambda_depraz.c
aefcf288a995master

Last Change
Raw File

lambda_depraz.c
View Options

#ifdef PY_INTERFACE
#include <Python.h>
#else
#include <stdio.h>
#include <string.h>
#endif
#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>
#ifdef COOLING
#ifdef LAMBDA_DEPRAZ
#include <hdf5.h>
int withLinInterpolation = 1;
int withDebugMessages = 1;
int Z_global_index_solar = 4;
/*******************************************************************************/
// Update gloval variable All.*_TABLE with current values
int updateCoolingTable(){
printf("Updating cooling tables...\n");
//FIXME (make general use of type_name)
char* type_name = "/Total_Metals";
//hdf5 tables location
char* tables_dir = "/home/epfl/revaz/code/gear/PyCool/tables_wiersma/coolingtables/";
if(withDebugMessages == 1){
printf("tables location: %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 1;
}
/*********************** z ***********************/
float z = 0.0;
#ifndef PY_INTERFACE
float a = get_a_from_CosmicTime(All.Time);
float actual_z = get_Redshift_from_a(a);
#else
float actual_z = 4.0;
#endif
if(withDebugMessages == 1){
printf("Redshift z = %g\n", actual_z);
}
//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 1;
}
//determine the file closest to the redshift value given as input
diff = fabs(z_file - actual_z);
if(diff < min){
file_name = content->d_name;
z = z_file;
min = diff;
}
}
}
if(withDebugMessages == 1){
printf("loading data from file %s\n", file_name);
}
All.CURRENT_TABLE_REDSHIFT = z;
//close directory
err = closedir(dir);
if (err != 0){
printf("an error occured while closing directory %s\n", tables_dir);
return 1;
}
//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);
/*********************** T ***********************/
char* T_key = "/Temperature_bins";
char T_table_key[(int)strlen(type_name)+(int)strlen(T_key)+1];
strcpy(T_table_key, type_name);
strcat(T_table_key, T_key);
loadDataInTable1D(table, T_table_key, &All.TEMPERATURE_TABLES, &All.SIZE_TEMPERATURE_TABLES);
/********************* rho_H *********************/
char* rho_H_key = "/Hydrogen_density_bins";
char rho_H_table_key[(int)strlen(type_name)+(int)strlen(rho_H_key)+1];
strcpy(rho_H_table_key, type_name);
strcat(rho_H_table_key, rho_H_key);
loadDataInTable1D(table, rho_H_table_key, &All.HYDROGEN_TABLES, &All.SIZE_HYDROGEN_TABLES);
/********************** nHe **********************/
char* nHe_key = "/Helium_mass_fraction_bins";
//or
//char* nHe_key = "/Helium_number_ratio_bins";
char nHe_table_key[(int)strlen("/Metal_free")+(int)strlen(nHe_key)+1];
strcpy(nHe_table_key, "/Metal_free");
strcat(nHe_table_key, nHe_key);
loadDataInTable1D(table, nHe_table_key, &All.HELIUM_ABOUNDANCE_TABLES, &All.SIZE_HELIUM_ABOUNDANCE_TABLES);
/********************* ne/n_H *********************/
char* ne_over_nH_key = "/Electron_density_over_n_h";
char ne_over_nH_table_key[(int)strlen("/Metal_free")+(int)strlen(ne_over_nH_key)+1];
strcpy(ne_over_nH_table_key, "/Metal_free");
strcat(ne_over_nH_table_key, ne_over_nH_key);
loadDataInTable3D(table, ne_over_nH_table_key, &All.ELECTRON_DENSITY_OVER_N_H_TABLES);
/******************** ne/n_H (Solar) ********************/
char* ne_over_nH_solar_key = "/Electron_density_over_n_h";
char ne_over_nH_solar_table_key[(int)strlen("/Solar")+(int)strlen(ne_over_nH_solar_key)+1];
strcpy(ne_over_nH_solar_table_key, "/Solar");
strcat(ne_over_nH_solar_table_key, ne_over_nH_solar_key);
loadDataInTable2D(table, ne_over_nH_solar_table_key, &All.ELECTRON_DENSITY_OVER_N_H_TABLES_SOLAR);
/********************* lambda *********************/
//load corresponding lambda value for metal
char* lambda_key = "/Net_cooling"; // /!\ sometimes written "/Net_Cooling"
char lambda_table_key[(int)strlen(type_name)+(int)strlen(lambda_key)+1];
strcpy(lambda_table_key, type_name);
strcat(lambda_table_key, lambda_key);
loadDataInTable2D(table, lambda_table_key, &All.COOLING_TABLES_TOTAL_METAL);
/**************** lambda Metal Free ****************/
//load corresponding lambda value for metal free
char* lambda_key_mf = "/Net_Cooling"; // /!\ sometimes written "/Net_cooling"
char lambda_table_key_mf[(int)strlen("/Metal_free")+(int)strlen(lambda_key_mf)+1];
strcpy(lambda_table_key_mf, "/Metal_free");
strcat(lambda_table_key_mf, lambda_key_mf);
loadDataInTable3D(table, lambda_table_key_mf, &All.COOLING_TABLES_METAL_FREE);
//close file
H5Fclose(table);
return 0;
}
/*******************************************************************************/
// Compute the resulting value of the cooling function using the data stored in
// the global variables
//
// variables ending with "_in" are user input, other variables are closest match
// in the hdf5 tables
float computeLambda(float rho_H_in, float T_in, float nHe_in, float metalicity){
float lambda_metal_free = 0.0;
float lambda_metals = 0.0;
float total_lambda = 0.0;
float ne_over_nH = 0.0;
float ne_over_nH_solar = 0.0;
float epsilon_diff = 1.0e-7;
float epsilon = 1.0e-38;
//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};
closestMatch1D(All.TEMPERATURE_TABLES, All.SIZE_TEMPERATURE_TABLES, T_in, T, T_index);
/********************* rho_H *********************/
float rho_H [2] = {0.0, 0.0};
int rho_H_index[2] = {0,0};
closestMatch1D(All.HYDROGEN_TABLES, All.SIZE_HYDROGEN_TABLES, rho_H_in, rho_H, rho_H_index);
/********************** nHe **********************/
float nHe [2] = {0.0, 0.0};
int nHe_index[2] = {0,0};
closestMatch1D(All.HELIUM_ABOUNDANCE_TABLES, All.SIZE_HELIUM_ABOUNDANCE_TABLES, nHe_in, nHe, nHe_index);
/********************* ne/n_H *********************/
//read electron density value from indices position
ne_over_nH = All.ELECTRON_DENSITY_OVER_N_H_TABLES[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)
ne_over_nH += (All.ELECTRON_DENSITY_OVER_N_H_TABLES[nHe_index[0]][T_index[0]][rho_H_index[0]]
- All.ELECTRON_DENSITY_OVER_N_H_TABLES[nHe_index[1]][T_index[0]][rho_H_index[0]])
/ fmax(nHe[0]-nHe[1], epsilon_diff) * (nHe_in - nHe[0]);
ne_over_nH += (All.ELECTRON_DENSITY_OVER_N_H_TABLES[nHe_index[0]][T_index[0]][rho_H_index[0]]
- All.ELECTRON_DENSITY_OVER_N_H_TABLES[nHe_index[0]][T_index[1]][rho_H_index[0]])
/ fmax(T[0]-T[1], epsilon_diff) * (T_in - T[0]);
ne_over_nH += (All.ELECTRON_DENSITY_OVER_N_H_TABLES[nHe_index[0]][T_index[0]][rho_H_index[0]]
- All.ELECTRON_DENSITY_OVER_N_H_TABLES[nHe_index[0]][T_index[0]][rho_H_index[1]])
/ fmax(rho_H[0]-rho_H[1], epsilon_diff) * (rho_H_in - rho_H[0]);
}
//contribution of metals
ne_over_nH_solar = All.ELECTRON_DENSITY_OVER_N_H_TABLES_SOLAR[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)
ne_over_nH_solar += (All.ELECTRON_DENSITY_OVER_N_H_TABLES_SOLAR[T_index[0]][rho_H_index[0]]
- All.ELECTRON_DENSITY_OVER_N_H_TABLES_SOLAR[T_index[1]][rho_H_index[0]])
/ fmax(T[0]-T[1], epsilon_diff) * (T_in - T[0]);
ne_over_nH_solar += (All.ELECTRON_DENSITY_OVER_N_H_TABLES_SOLAR[T_index[0]][rho_H_index[0]]
- All.ELECTRON_DENSITY_OVER_N_H_TABLES_SOLAR[T_index[0]][rho_H_index[1]])
/ fmax(rho_H[0]-rho_H[1], epsilon_diff) * (rho_H_in - rho_H[0]);
}
/********************* lambda *********************/
//read lambda value from indices position
//contribution of hydrogen and helium
lambda_metal_free = All.COOLING_TABLES_METAL_FREE[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_metal_free += (All.COOLING_TABLES_METAL_FREE[nHe_index[0]][T_index[0]][rho_H_index[0]]
- All.COOLING_TABLES_METAL_FREE[nHe_index[1]][T_index[0]][rho_H_index[0]])
/ fmax(nHe[0]-nHe[1], epsilon_diff) * (nHe_in - nHe[0]);
lambda_metal_free += (All.COOLING_TABLES_METAL_FREE[nHe_index[0]][T_index[0]][rho_H_index[0]]
- All.COOLING_TABLES_METAL_FREE[nHe_index[0]][T_index[1]][rho_H_index[0]])
/ fmax(T[0]-T[1], epsilon_diff) * (T_in - T[0]);
lambda_metal_free += (All.COOLING_TABLES_METAL_FREE[nHe_index[0]][T_index[0]][rho_H_index[0]]
- All.COOLING_TABLES_METAL_FREE[nHe_index[0]][T_index[0]][rho_H_index[1]])
/ fmax(rho_H[0]-rho_H[1], epsilon_diff) * (rho_H_in - rho_H[0]);
}
//contribution of metals
lambda_metals = All.COOLING_TABLES_TOTAL_METAL[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_metals += (All.COOLING_TABLES_TOTAL_METAL[T_index[0]][rho_H_index[0]]
- All.COOLING_TABLES_TOTAL_METAL[T_index[1]][rho_H_index[0]])
/ fmax(T[0]-T[1], epsilon_diff) * (T_in - T[0]);
lambda_metals += (All.COOLING_TABLES_TOTAL_METAL[T_index[0]][rho_H_index[0]]
- All.COOLING_TABLES_TOTAL_METAL[T_index[0]][rho_H_index[1]])
/ fmax(rho_H[0]-rho_H[1], epsilon_diff) * (rho_H_in - rho_H[0]);
}
// formula from taken from:
// "The effect of photoionization on the cooling rates of enriched,
// astrophysical plasmas"
// Robert P.C. Wiersma, J. Schaye & B.D. Smith
// p.101, eq. (5)
// http://onlinelibrary.wiley.com/doi/10.1111/j.1365-2966.2008.14191.x/pdf
//
#ifndef PY_INTERFACE
total_lambda = lambda_metal_free + lambda_metals * ne_over_nH/ne_over_nH_solar
* metalicity/fmax(chimie_extraheader.SolarAbundances[Z_global_index_solar], epsilon);
#else
//chimie_extraheader.SolarAbundances[] not defined, use constant value instead
total_lambda = lambda_metal_free + lambda_metals * ne_over_nH/ne_over_nH_solar
* metalicity/0.02;
#endif
//denormalized value (not used for now...)
/*
if(withLinInterpolation == 1){
total_lambda *= pow( ((All.HYDROGEN_TABLES[rho_H_index[0]] - All.HYDROGEN_TABLES[rho_H_index[1]])
/ fmax(rho_H[0]-rho_H[1], epsilon_diff) * (rho_H_in - rho_H[0])) ,2);
} else {
total_lambda *= pow(All.HYDROGEN_TABLES[rho_H_index[0]], 2);
}
*/
if(withDebugMessages == 1){
printf("lambda = %g\n", total_lambda);
}
return total_lambda;
}
/*******************************************************************************/
// write data from hdf5 tables to global variable TABLE (1D) and write the
// corresponding array size in SIZE_TABLE
void loadDataInTable1D(hid_t table, char* table_key, float** TABLE, int* SIZE_TABLE){
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);
if(withDebugMessages == 1){
printf("%s: rank %d, dimension %d \n", table_key, rank, (int)size[0]);
}
*TABLE = (float*)malloc(size[0]*sizeof(float));
if(TABLE == NULL){
printf("\nMemory allocation failed in loadDataInTable1D()\n");
return;
}
//store data in global variables
*SIZE_TABLE = size[0];
int k=0;
for(k=0; k<size[0]; k++){
(*TABLE)[k] = data[k];
}
//close instances
H5Sclose(memspace);
H5Sclose(dataspace);
H5Dclose(dataset);
}
/*******************************************************************************/
// write data from hdf5 tables to global variable TABLE (2D)
void loadDataInTable2D(hid_t table, char* table_key, float*** TABLE){
hid_t dataset;
hid_t dataspace;
hid_t memspace;
herr_t status;
hsize_t dims[2];
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, dims, NULL);// save dataspace dimensions in dims
memspace = H5Screate_simple(rank, dims, NULL); // determine memory space required
float* temp;
int i = 0;
int j = 0;
//format of the data: data[T][nH]
float data[dims[0]][dims[1]];
status = H5Dread(dataset, H5T_NATIVE_FLOAT, memspace, dataspace, H5P_DEFAULT, data);
//allocate contiguous memory to enable the later use of MPI_Bcast()
temp = (float*)malloc(dims[0]*dims[1]*sizeof(float));
//allocate memory for TABLE accordingly
*TABLE = (float**)malloc(dims[0]*sizeof(float*));
for(i=0; i<dims[0]; i++){
(*TABLE)[i] = &(temp[i*dims[1]]);
}
if(TABLE == NULL){
printf("\nMemory allocation failed in loadDataInTable2D()\n");
free(temp);
return;
}
for(i=0; i<dims[0]; i++){
for(j=0; j<dims[1]; j++){
(*TABLE)[i][j] = data[i][j];
}
}
if(withDebugMessages == 1){
printf("%s loaded with dimensions: %d x %d\n", table_key, (int)dims[0], (int)dims[1]);
}
//close instances
H5Sclose(memspace);
H5Sclose(dataspace);
H5Dclose(dataset);
}
/*******************************************************************************/
// write data from hdf5 tables to global variable TABLE (3D)
void loadDataInTable3D(hid_t table, char* table_key, float**** TABLE){
hid_t dataset;
hid_t dataspace;
hid_t memspace;
herr_t status;
hsize_t dims[3];
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, dims, NULL);// save dataspace dimensions in dims
memspace = H5Screate_simple(rank, dims, NULL); // determine memory space required
float** temp1;
float* temp2;
int i = 0;
int j = 0;
int k = 0;
//format of the data: data[nHe][T][nH]
float data[dims[0]][dims[1]][dims[2]];
status = H5Dread(dataset, H5T_NATIVE_FLOAT, memspace, dataspace, H5P_DEFAULT, data);
//allocate contiguous memory to enable the later use of MPI_Bcast()
temp1 = (float**)malloc(dims[0]*dims[1]*sizeof(float*));
temp2 = (float*)malloc(dims[0]*dims[1]*dims[2]*sizeof(float));
for(i=0, j=0; i<dims[0]*dims[1]; i++, j++){
temp1[i] = &(temp2[j*dims[2]]);
}
*TABLE = (float***)malloc(dims[0]*sizeof(float**));
for(i=0, j=0; j<dims[0]; i++, j++){
(*TABLE)[i] = &(temp1[j*dims[1]]);
}
if(TABLE == NULL){
printf("\nMemory allocation failed in loadDataInTable3D()\n");
free(temp1);
free(temp2);
return;
}
for(i=0; i<dims[0]; i++){
for(j=0; j<dims[1]; j++){
for(k=0; k<dims[2]; k++){
(*TABLE)[i][j][k] = data[i][j][k];
}
}
}
if(withDebugMessages == 1){
printf("%s loaded with dimensions: %d x %d x %d\n",
table_key, (int)dims[0], (int)dims[1], (int)dims[2]);
}
//close instances
H5Sclose(memspace);
H5Sclose(dataspace);
H5Dclose(dataset);
}
/*******************************************************************************/
//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(float* TABLE, int SIZE, float input, float match[2], int index[2]){
//determine the index of the closest values (above and below)
int k = 0;
int oldk = 0;
float min1 = 1e30;
float min2 = 1e30;
float oldmin = 1e30;
float oldmatch = 0.0;
float diff = 0.0;
//the values having the smallest difference
//to the input are computed and the indeces are stored
for(k=0; k<SIZE; k++){
diff = fabs(input - TABLE[k]);
if(diff < min1){
oldmatch = match[0];
match[0] = TABLE[k];
oldk = index[0];
index[0] = k;
oldmin = min1;
min1 = diff;
}
if(diff < min2 && diff != min1){
if(min2 < oldmin){
match[1] = TABLE[k];
index[1] = k;
min2 = diff;
} else {
match[1] = oldmatch;
index[1] = oldk;
min2 = oldmin;
}
}
}
if(withDebugMessages == 1){
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]);
}
}
/*******************************************************************************/
// Check if the file that is read currently still has the correct redshift value
void checkRedshiftForUpdate(){
#ifndef PY_INTERFACE
float a = get_a_from_CosmicTime(All.Time);
float z = get_Redshift_from_a(a);
#else
float z = 0.0;
#endif
if(fabs(All.CURRENT_TABLE_REDSHIFT - z) >= 0.04){
freeMemory();
if(ThisTask == 0){
updateCoolingTable();
}
BroadcastTablesToAllFromMaster();
}
}
/*******************************************************************************/
// Broadcast tables to all other procs from proc 0
void BroadcastTablesToAllFromMaster(){
if(withDebugMessages == 1 && ThisTask == 0){
printf("broadcasting cooling tables...\n");
}
//broadcast array sizes
MPI_Bcast(&All.SIZE_HYDROGEN_TABLES, 1, MPI_FLOAT, 0, MPI_COMM_WORLD);
MPI_Bcast(&All.SIZE_TEMPERATURE_TABLES, 1, MPI_FLOAT, 0, MPI_COMM_WORLD);
MPI_Bcast(&All.SIZE_HELIUM_ABOUNDANCE_TABLES, 1, MPI_FLOAT, 0, MPI_COMM_WORLD);
//for safety
MPI_Barrier(MPI_COMM_WORLD);
//allocate memory for all procs (except proc 0)
if(ThisTask != 0){
//temporary arrays used to allocate contiguous memory
float* temp1;
float* temp2;
float* temp3;
float** temp4;
float* temp5;
float** temp6;
int i=0;
int j=0;
// 1D arrays
All.HYDROGEN_TABLES = (float*)malloc(All.SIZE_HYDROGEN_TABLES*sizeof(float));
All.TEMPERATURE_TABLES = (float*)malloc(All.SIZE_TEMPERATURE_TABLES*sizeof(float));
All.HELIUM_ABOUNDANCE_TABLES = (float*)malloc(All.SIZE_HELIUM_ABOUNDANCE_TABLES*sizeof(float));
// 2D arrays
// allocate contiguous memory
temp1 = (float*)malloc(All.SIZE_TEMPERATURE_TABLES*All.SIZE_HYDROGEN_TABLES*sizeof(float));
temp2 = (float*)malloc(All.SIZE_TEMPERATURE_TABLES*All.SIZE_HYDROGEN_TABLES*sizeof(float));
//allocate arrays
All.COOLING_TABLES_TOTAL_METAL = (float**)malloc(All.SIZE_TEMPERATURE_TABLES*sizeof(float*));
for(i=0; i<All.SIZE_TEMPERATURE_TABLES; i++){
All.COOLING_TABLES_TOTAL_METAL[i] = &(temp1[i*All.SIZE_HYDROGEN_TABLES]);
}
All.ELECTRON_DENSITY_OVER_N_H_TABLES_SOLAR = (float**)malloc(All.SIZE_TEMPERATURE_TABLES*sizeof(float*));
for(i=0; i<All.SIZE_TEMPERATURE_TABLES; i++){
All.ELECTRON_DENSITY_OVER_N_H_TABLES_SOLAR[i] = &(temp2[i*All.SIZE_HYDROGEN_TABLES]);
}
// 3D arrays
// allocate contiguous memory
temp3 = (float*)malloc(All.SIZE_HELIUM_ABOUNDANCE_TABLES*All.SIZE_TEMPERATURE_TABLES*All.SIZE_HYDROGEN_TABLES*sizeof(float));
temp4 = (float**)malloc(All.SIZE_HELIUM_ABOUNDANCE_TABLES*All.SIZE_TEMPERATURE_TABLES*sizeof(float*));
for(i=0, j=0; i<All.SIZE_HELIUM_ABOUNDANCE_TABLES*All.SIZE_TEMPERATURE_TABLES; i++, j++){
temp4[i] = &(temp3[j*All.SIZE_HYDROGEN_TABLES]);
}
All.COOLING_TABLES_METAL_FREE = (float***)malloc(All.SIZE_HELIUM_ABOUNDANCE_TABLES*sizeof(float**));
for(i=0, j=0; j<All.SIZE_HELIUM_ABOUNDANCE_TABLES; i++, j++){
All.COOLING_TABLES_METAL_FREE[i] = &(temp4[j*All.SIZE_TEMPERATURE_TABLES]);
}
temp5 = (float*)malloc(All.SIZE_HELIUM_ABOUNDANCE_TABLES*All.SIZE_TEMPERATURE_TABLES*All.SIZE_HYDROGEN_TABLES*sizeof(float));
temp6 = (float**)malloc(All.SIZE_HELIUM_ABOUNDANCE_TABLES*All.SIZE_TEMPERATURE_TABLES*sizeof(float*));
for(i=0, j=0; i<All.SIZE_HELIUM_ABOUNDANCE_TABLES*All.SIZE_TEMPERATURE_TABLES; i++, j++){
temp6[i] = &(temp5[j*All.SIZE_HYDROGEN_TABLES]);
}
All.ELECTRON_DENSITY_OVER_N_H_TABLES = (float***)malloc(All.SIZE_HELIUM_ABOUNDANCE_TABLES*sizeof(float**));
for(i=0, j=0; j<All.SIZE_HELIUM_ABOUNDANCE_TABLES; i++, j++){
All.ELECTRON_DENSITY_OVER_N_H_TABLES[i] = &(temp6[j*All.SIZE_TEMPERATURE_TABLES]);
}
}
//broadcast global arrays to all procs from 0
MPI_Bcast(&(All.COOLING_TABLES_TOTAL_METAL[0][0]),
All.SIZE_TEMPERATURE_TABLES*All.SIZE_HYDROGEN_TABLES, MPI_FLOAT,0, MPI_COMM_WORLD);
MPI_Bcast(&(All.COOLING_TABLES_METAL_FREE[0][0][0]),
All.SIZE_HELIUM_ABOUNDANCE_TABLES*All.SIZE_TEMPERATURE_TABLES*All.SIZE_HYDROGEN_TABLES, MPI_FLOAT, 0, MPI_COMM_WORLD);
MPI_Bcast(&(All.ELECTRON_DENSITY_OVER_N_H_TABLES_SOLAR[0][0]),
All.SIZE_TEMPERATURE_TABLES*All.SIZE_HYDROGEN_TABLES, MPI_FLOAT, 0, MPI_COMM_WORLD);
MPI_Bcast(&(All.ELECTRON_DENSITY_OVER_N_H_TABLES[0][0][0]),
All.SIZE_HELIUM_ABOUNDANCE_TABLES*All.SIZE_TEMPERATURE_TABLES*All.SIZE_HYDROGEN_TABLES, MPI_FLOAT, 0, MPI_COMM_WORLD);
MPI_Bcast(All.HYDROGEN_TABLES, All.SIZE_HYDROGEN_TABLES, MPI_FLOAT, 0, MPI_COMM_WORLD);
MPI_Bcast(All.TEMPERATURE_TABLES, All.SIZE_TEMPERATURE_TABLES, MPI_FLOAT, 0, MPI_COMM_WORLD);
MPI_Bcast(All.HELIUM_ABOUNDANCE_TABLES, All.SIZE_HELIUM_ABOUNDANCE_TABLES, MPI_FLOAT, 0, MPI_COMM_WORLD);
MPI_Bcast(&All.CURRENT_TABLE_REDSHIFT, 1, MPI_FLOAT, 0, MPI_COMM_WORLD);
if(withDebugMessages == 1 && ThisTask == 0){
printf("cooling tables broadcasted to all procs\n");
}
}
/*******************************************************************************/
// Function that free memory to avoid overflows
// copy in allocate.c
int freeMemory(){
int error = 0;
// free 3D arrays
if(All.COOLING_TABLES_METAL_FREE){
free(&((All.COOLING_TABLES_METAL_FREE)[0][0]));
free(All.COOLING_TABLES_METAL_FREE);
All.COOLING_TABLES_METAL_FREE = NULL;
}
if(All.ELECTRON_DENSITY_OVER_N_H_TABLES){
free(&((All.ELECTRON_DENSITY_OVER_N_H_TABLES)[0][0]));
free(All.ELECTRON_DENSITY_OVER_N_H_TABLES);
All.ELECTRON_DENSITY_OVER_N_H_TABLES = NULL;
}
// free 2D arrays
if(All.COOLING_TABLES_TOTAL_METAL){
free(All.COOLING_TABLES_TOTAL_METAL);
All.COOLING_TABLES_TOTAL_METAL = NULL;
}
if(All.ELECTRON_DENSITY_OVER_N_H_TABLES_SOLAR){
free(All.ELECTRON_DENSITY_OVER_N_H_TABLES_SOLAR);
All.ELECTRON_DENSITY_OVER_N_H_TABLES_SOLAR = NULL;
}
// free 1D arrays
if(All.HYDROGEN_TABLES){
free(All.HYDROGEN_TABLES);
All.HYDROGEN_TABLES = NULL;
}
if(All.TEMPERATURE_TABLES){
free(All.TEMPERATURE_TABLES);
All.TEMPERATURE_TABLES = NULL;
}
if(All.HELIUM_ABOUNDANCE_TABLES){
free(All.HELIUM_ABOUNDANCE_TABLES);
All.HELIUM_ABOUNDANCE_TABLES = NULL;
}
if(withDebugMessages == 1 && ThisTask == 0){
printf("cooling table memory freed\n");
}
return error;
}
/************************************* Python Interface ************************************/
#ifdef PY_INTERFACE
/*****************************************************/
/* functions declarations */
/*****************************************************/
static PyObject* updateCoolingTableInterface(PyObject *self, PyObject *args) {
int error_status = 0;
error_status = updateCoolingTable();
return Py_BuildValue("d", error_status);
}
static PyObject* BroadcastTablesToAllFromMasterInterface(PyObject *self, PyObject *args) {
BroadcastTablesToAllFromMaster();
return Py_BuildValue("d", 0);
}
static PyObject* computeLambdaInterface(PyObject *self, PyObject *args) {
float rho_H_in;
float T_in;
float nHe_in;
float metalicity;
float lambda = 0.0;
if (!PyArg_ParseTuple(args, "ffff", &rho_H_in, &T_in, &nHe_in, &metalicity)) {
return NULL;
}
//printf("%s\n","computing function with parameters:");
//printf("%s: %g\n","hydrogen density", rho_H_in);
//printf("%s: %g\n","temperature", T_in);
//printf("%s: %g\n","helium aboundance", nHe_in);
//printf("%s: %g\n","metalicity", metalicity);
lambda = computeLambda(rho_H_in, T_in, nHe_in, metalicity);
return Py_BuildValue("f",lambda);
}
static PyObject* freeMemoryInterface(PyObject *self, PyObject *args) {
int error_status;
error_status = freeMemory();
return Py_BuildValue("d",error_status);
}
/*****************************************************/
/* Method table */
/*****************************************************/
static PyMethodDef computeLambdaMethods[] = {
{"updateCoolingTable", updateCoolingTableInterface, METH_VARARGS,
"Update cooling table data using HDF5 tables."},
{"BroadcastTablesToAllFromMaster", BroadcastTablesToAllFromMasterInterface, METH_VARARGS,
"Broadcast tables to all other procs from proc 0."},
{"computeLambda", computeLambdaInterface, METH_VARARGS,
"Compute the cooling from given inputs."},
{"freeMemory", freeMemoryInterface, METH_VARARGS,
"Free allocated memory."},
{NULL, NULL, 0, NULL} /* Sentinel */
};
/*****************************************************/
/* initialization */
/*****************************************************/
//void initcompute_lambda_interface(void)
PyMODINIT_FUNC initcompute_lambda_interface(void) {
(void) Py_InitModule("compute_lambda_interface", computeLambdaMethods);
}
#endif //PY_INTERFACE
#endif //LAMBDA_DEPRAZ
#endif //COOLING
c4science · Help