diff --git a/examples/main.c b/examples/main.c
index 9a8a6dbe14005cf23153b2530a736143595fd7d7..c5b00e1c5bf8a751615dd743b34d51d61acf8726 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -130,7 +130,7 @@ int main(int argc, char *argv[]) {
struct clocks_time tic, toc;
struct engine e;
-
+
/* Structs used by the engine. Declare now to make sure these are always in
* scope. */
struct chemistry_global_data chemistry;
diff --git a/src/chemistry/EAGLE/chemistry.h b/src/chemistry/EAGLE/chemistry.h
index 535a8a4d649b5d1e8342f72a6bcef8faf3c19b2d..5d7e0e8883d81996c1f8e10288f7d63bdeb9c365 100644
--- a/src/chemistry/EAGLE/chemistry.h
+++ b/src/chemistry/EAGLE/chemistry.h
@@ -41,10 +41,10 @@
* @brief Return a string containing the name of a given #chemistry_element.
*/
//__attribute__((always_inline)) INLINE static const char*
-__attribute__((always_inline)) INLINE const char*
-chemistry_get_element_name(enum chemistry_element elem) {
+__attribute__((always_inline)) INLINE const char* chemistry_get_element_name(
+ enum chemistry_element elem) {
- //static const char* chemistry_element_names[chemistry_element_count] = {
+ // static const char* chemistry_element_names[chemistry_element_count] = {
const char* chemistry_element_names[chemistry_element_count] = {
"Hydrogen", "Helium", "Carbon", "Nitrogen", "Oxygen",
"Neon", "Magnesium", "Silicon", "Iron"};
@@ -139,7 +139,8 @@ static INLINE void chemistry_init_backend(struct swift_params* parameter_file,
parser_get_param_float(parameter_file, "EAGLEChemistry:InitMetallicity");
/* Read the individual mass fractions */
- for (enum chemistry_element elem = chemistry_element_H; elem < chemistry_element_count; ++elem) {
+ for (enum chemistry_element elem = chemistry_element_H;
+ elem < chemistry_element_count; ++elem) {
char buffer[50];
sprintf(buffer, "EAGLEChemistry:InitAbundance_%s",
chemistry_get_element_name((enum chemistry_element)elem));
@@ -174,23 +175,47 @@ static INLINE void chemistry_print_backend(
* @param p particle struct
* @param elem enum value of element
*/
-__attribute__((always_inline)) INLINE static double chemistry_get_number_density(const struct part* restrict p, const struct cosmology* restrict cosmo, enum chemistry_element elem, const struct phys_const* restrict internal_const) {
+__attribute__((always_inline)) INLINE static double
+chemistry_get_number_density(const struct part* restrict p,
+ const struct cosmology* restrict cosmo,
+ enum chemistry_element elem,
+ const struct phys_const* restrict internal_const) {
double number_density;
int atomic_number;
- switch(elem){
- case chemistry_element_H : atomic_number = 1; break;
- case chemistry_element_He: atomic_number = 4; break;
- case chemistry_element_C : atomic_number = 12; break;
- case chemistry_element_N : atomic_number = 14; break;
- case chemistry_element_O : atomic_number = 16; break;
- case chemistry_element_Ne: atomic_number = 20; break;
- case chemistry_element_Mg: atomic_number = 24; break;
- case chemistry_element_Si: atomic_number = 28; break;
- case chemistry_element_Fe: atomic_number = 56; break;
+ switch (elem) {
+ case chemistry_element_H:
+ atomic_number = 1;
+ break;
+ case chemistry_element_He:
+ atomic_number = 4;
+ break;
+ case chemistry_element_C:
+ atomic_number = 12;
+ break;
+ case chemistry_element_N:
+ atomic_number = 14;
+ break;
+ case chemistry_element_O:
+ atomic_number = 16;
+ break;
+ case chemistry_element_Ne:
+ atomic_number = 20;
+ break;
+ case chemistry_element_Mg:
+ atomic_number = 24;
+ break;
+ case chemistry_element_Si:
+ atomic_number = 28;
+ break;
+ case chemistry_element_Fe:
+ atomic_number = 56;
+ break;
}
- double element_mass = internal_const->const_proton_mass*atomic_number;
- number_density = p->chemistry_data.metal_mass_fraction[elem]*hydro_get_physical_density(p,cosmo)/element_mass;
- //number_density = p->chemistry_data.metal_mass_fraction[elem]*hydro_get_comoving_density(p)/element_mass;
+ double element_mass = internal_const->const_proton_mass * atomic_number;
+ number_density = p->chemistry_data.metal_mass_fraction[elem] *
+ hydro_get_physical_density(p, cosmo) / element_mass;
+ // number_density =
+ // p->chemistry_data.metal_mass_fraction[elem]*hydro_get_comoving_density(p)/element_mass;
return number_density;
}
diff --git a/src/chemistry/EAGLE/chemistry_io.h b/src/chemistry/EAGLE/chemistry_io.h
index b55af73190cfc007b70133bbb988f6962e7a6209..d7139141978e81b55e1e9de70fae7809cface0cb 100644
--- a/src/chemistry/EAGLE/chemistry_io.h
+++ b/src/chemistry/EAGLE/chemistry_io.h
@@ -106,7 +106,8 @@ INLINE static int chemistry_write_particles(const struct part* parts,
INLINE static void chemistry_write_flavour(hid_t h_grp) {
io_write_attribute_s(h_grp, "Chemistry Model", "EAGLE");
- for (enum chemistry_element elem = chemistry_element_H; elem < chemistry_element_count; ++elem) {
+ for (enum chemistry_element elem = chemistry_element_H;
+ elem < chemistry_element_count; ++elem) {
char buffer[20];
sprintf(buffer, "Element %d", elem);
io_write_attribute_s(
diff --git a/src/cooling.c b/src/cooling.c
index ddf5e9e74cfd13746ef48e6d2cf4e07f2f42f7be..93d1799c8e5b7b7be774f5cbd0bba766e3b29584 100644
--- a/src/cooling.c
+++ b/src/cooling.c
@@ -29,7 +29,6 @@
#include "cooling.h"
#include "restart.h"
-
/**
* @brief Initialises the cooling properties.
*
diff --git a/src/cooling/EAGLE/cooling.h b/src/cooling/EAGLE/cooling.h
index 8e6ba1ef6ba39ea2a842099111fbaba418e413a2..d7c44b846cad6a44bf9136d336828a7b14a7fa05 100644
--- a/src/cooling/EAGLE/cooling.h
+++ b/src/cooling/EAGLE/cooling.h
@@ -29,21 +29,21 @@
/* Some standard headers. */
#include <float.h>
-#include <math.h>
#include <hdf5.h>
+#include <math.h>
#include <time.h>
/* Local includes. */
+#include "chemistry.h"
#include "cooling_struct.h"
+#include "eagle_cool_tables.h"
#include "error.h"
#include "hydro.h"
-#include "chemistry.h"
#include "io_properties.h"
#include "parser.h"
#include "part.h"
#include "physical_constants.h"
#include "units.h"
-#include "eagle_cool_tables.h"
/* number of calls to eagle cooling rate */
extern int n_eagle_cooling_rate_calls_1;
@@ -64,14 +64,14 @@ enum hdf5_allowed_types {
};
/**
- * @brief Returns the 1d index of element with 2d indices i,j
+ * @brief Returns the 1d index of element with 2d indices i,j
* from a flattened 2d array in row major order
*
* @param i, j Indices of element of interest
* @param nx, ny Sizes of array dimensions
*/
__attribute__((always_inline)) INLINE int row_major_index_2d(int i, int j,
- int nx, int ny) {
+ int nx, int ny) {
int index = i * ny + j;
#ifdef SWIFT_DEBUG_CHECKS
assert(i < nx);
@@ -81,15 +81,15 @@ __attribute__((always_inline)) INLINE int row_major_index_2d(int i, int j,
}
/**
- * @brief Returns the 1d index of element with 3d indices i,j,k
+ * @brief Returns the 1d index of element with 3d indices i,j,k
* from a flattened 3d array in row major order
*
* @param i, j, k Indices of element of interest
* @param nx, ny, nz Sizes of array dimensions
*/
__attribute__((always_inline)) INLINE int row_major_index_3d(int i, int j,
- int k, int nx,
- int ny, int nz) {
+ int k, int nx,
+ int ny, int nz) {
int index = i * ny * nz + j * nz + k;
#ifdef SWIFT_DEBUG_CHECKS
assert(i < nx);
@@ -107,12 +107,12 @@ __attribute__((always_inline)) INLINE int row_major_index_3d(int i, int j,
* @param nx, ny, nz, nw Sizes of array dimensions
*/
__attribute__((always_inline)) INLINE int row_major_index_4d(int i, int j,
- int k, int l,
- int nx, int ny,
- int nz, int nw) {
+ int k, int l,
+ int nx, int ny,
+ int nz, int nw) {
int index = i * ny * nz * nw + j * nz * nw + k * nw + l;
#ifdef SWIFT_DEBUG_CHECKS
- //printf("Eagle cooling.h j, ny %d %d\n",j,ny);
+ // printf("Eagle cooling.h j, ny %d %d\n",j,ny);
assert(i < nx);
assert(j < ny);
assert(k < nz);
@@ -121,7 +121,6 @@ __attribute__((always_inline)) INLINE int row_major_index_4d(int i, int j,
return index;
}
-
/*
* @brief This routine returns the position i of a value x in a 1D table and the
* displacement dx needed for the interpolation. The table is assumed to
@@ -134,7 +133,9 @@ __attribute__((always_inline)) INLINE int row_major_index_4d(int i, int j,
* is the greatest value in the table less than x
* @param dx Pointer to offset of x within table cell
*/
-__attribute__((always_inline)) INLINE void get_index_1d(float *table, int ntable, double x, int *i, float *dx) {
+__attribute__((always_inline)) INLINE void get_index_1d(float *table,
+ int ntable, double x,
+ int *i, float *dx) {
float dxm1;
const float EPS = 1.e-4;
@@ -148,26 +149,32 @@ __attribute__((always_inline)) INLINE void get_index_1d(float *table, int ntable
*dx = 1;
} else {
*i = (int)floor(((float)x - table[0]) * dxm1);
- if (*i >= ntable || *i < 0){
- printf("Eagle cooling.h i, ntable, x, table[0], dxm1 %d %d %.5e %.5e %.5e \n", *i, ntable, x, table[0], dxm1);
+ if (*i >= ntable || *i < 0) {
+ printf(
+ "Eagle cooling.h i, ntable, x, table[0], dxm1 %d %d %.5e %.5e %.5e "
+ "\n",
+ *i, ntable, x, table[0], dxm1);
fflush(stdout);
}
*dx = ((float)x - table[*i]) * dxm1;
}
}
-
/*
- * @brief This routine returns the position i of a value z in the redshift table and the
+ * @brief This routine returns the position i of a value z in the redshift table
+ * and the
* displacement dx needed for the interpolation.
*
- * @param z Redshift whose position within the redshift array we are interested in
+ * @param z Redshift whose position within the redshift array we are interested
+ * in
* @param z_index i Pointer to the index whose corresponding redshift
* is the greatest value in the redshift table less than x
* @param dx Pointer to offset of z within redshift cell
* @param cooling Pointer to cooling structure (contains redshift table)
*/
-__attribute__((always_inline)) INLINE void get_redshift_index(float z, int *z_index, float *dz, const struct cooling_function_data* restrict cooling) {
+__attribute__((always_inline)) INLINE void get_redshift_index(
+ float z, int *z_index, float *dz,
+ const struct cooling_function_data *restrict cooling) {
int i, iz;
if (get_redshift_index_first_call == 0) {
@@ -213,7 +220,6 @@ __attribute__((always_inline)) INLINE void get_redshift_index(float z, int *z_in
}
}
-
/*
* @brief Performs 1d interpolation
*
@@ -221,7 +227,8 @@ __attribute__((always_inline)) INLINE void get_redshift_index(float z, int *z_in
* @param i Index of cell we are interpolating
* @param dx Offset within cell
*/
-__attribute__((always_inline)) INLINE float interpol_1d(float *table, int i, float dx) {
+__attribute__((always_inline)) INLINE float interpol_1d(float *table, int i,
+ float dx) {
float result;
result = (1 - dx) * table[i] + dx * table[i + 1];
@@ -236,7 +243,8 @@ __attribute__((always_inline)) INLINE float interpol_1d(float *table, int i, flo
* @param i Index of cell we are interpolating
* @param dx Offset within cell
*/
-__attribute__((always_inline)) INLINE double interpol_1d_dbl(double *table, int i, float dx) {
+__attribute__((always_inline)) INLINE double interpol_1d_dbl(double *table,
+ int i, float dx) {
double result;
result = (1 - dx) * table[i] + dx * table[i + 1];
@@ -252,23 +260,35 @@ __attribute__((always_inline)) INLINE double interpol_1d_dbl(double *table, int
* @param dx,dy Offset within cell
* @param nx,ny Table dimensions
*/
-__attribute__((always_inline)) INLINE float interpol_2d(float *table, int i, int j, float dx, float dy, int nx, int ny) {
+__attribute__((always_inline)) INLINE float interpol_2d(float *table, int i,
+ int j, float dx,
+ float dy, int nx,
+ int ny) {
float result;
int index[4];
- index[0] = row_major_index_2d(i,j,nx,ny);
- index[1] = row_major_index_2d(i,j+1,nx,ny);
- index[2] = row_major_index_2d(i+1,j,nx,ny);
- index[3] = row_major_index_2d(i+1,j+1,nx,ny);
+ index[0] = row_major_index_2d(i, j, nx, ny);
+ index[1] = row_major_index_2d(i, j + 1, nx, ny);
+ index[2] = row_major_index_2d(i + 1, j, nx, ny);
+ index[3] = row_major_index_2d(i + 1, j + 1, nx, ny);
#ifdef SWIFT_DEBUG_CHECKS
- if(index[0] >= nx*ny || index[0] < 0) fprintf(stderr,"index 0 out of bounds %d, i,j %d, %d, table size %d\n", index[0],i,j,nx*ny);
- if(index[1] >= nx*ny || index[1] < 0) fprintf(stderr,"index 1 out of bounds %d, i,j %d, %d, table size %d\n", index[1],i,j+1,nx*ny);
- if(index[2] >= nx*ny || index[2] < 0) fprintf(stderr,"index 2 out of bounds %d, i,j %d, %d, table size %d\n", index[2],i+1,j,nx*ny);
- if(index[3] >= nx*ny || index[3] < 0) fprintf(stderr,"index 3 out of bounds %d, i,j %d, %d, table size %d\n", index[3],i+1,j+1,nx*ny);
+ if (index[0] >= nx * ny || index[0] < 0)
+ fprintf(stderr, "index 0 out of bounds %d, i,j %d, %d, table size %d\n",
+ index[0], i, j, nx * ny);
+ if (index[1] >= nx * ny || index[1] < 0)
+ fprintf(stderr, "index 1 out of bounds %d, i,j %d, %d, table size %d\n",
+ index[1], i, j + 1, nx * ny);
+ if (index[2] >= nx * ny || index[2] < 0)
+ fprintf(stderr, "index 2 out of bounds %d, i,j %d, %d, table size %d\n",
+ index[2], i + 1, j, nx * ny);
+ if (index[3] >= nx * ny || index[3] < 0)
+ fprintf(stderr, "index 3 out of bounds %d, i,j %d, %d, table size %d\n",
+ index[3], i + 1, j + 1, nx * ny);
#endif
- result = (1 - dx) * (1 - dy) * table[index[0]] + (1 - dx) * dy * table[index[1]] +
- dx * (1 - dy) * table[index[2]] + dx * dy * table[index[3]];
+ result = (1 - dx) * (1 - dy) * table[index[0]] +
+ (1 - dx) * dy * table[index[1]] + dx * (1 - dy) * table[index[2]] +
+ dx * dy * table[index[3]];
return result;
}
@@ -281,23 +301,33 @@ __attribute__((always_inline)) INLINE float interpol_2d(float *table, int i, int
* @param dx,dy Offset within cell
* @param nx,ny Table dimensions
*/
-__attribute__((always_inline)) INLINE double interpol_2d_dbl(double *table, int i, int j, double dx, double dy, int nx, int ny) {
+__attribute__((always_inline)) INLINE double interpol_2d_dbl(
+ double *table, int i, int j, double dx, double dy, int nx, int ny) {
double result;
int index[4];
- index[0] = row_major_index_2d(i,j,nx,ny);
- index[1] = row_major_index_2d(i,j+1,nx,ny);
- index[2] = row_major_index_2d(i+1,j,nx,ny);
- index[3] = row_major_index_2d(i+1,j+1,nx,ny);
+ index[0] = row_major_index_2d(i, j, nx, ny);
+ index[1] = row_major_index_2d(i, j + 1, nx, ny);
+ index[2] = row_major_index_2d(i + 1, j, nx, ny);
+ index[3] = row_major_index_2d(i + 1, j + 1, nx, ny);
#ifdef SWIFT_DEBUG_CHECKS
- if(index[0] >= nx*ny || index[0] < 0) fprintf(stderr,"index 0 out of bounds %d, i,j %d, %d, table size %d\n", index[0],i,j,nx*ny);
- if(index[1] >= nx*ny || index[1] < 0) fprintf(stderr,"index 1 out of bounds %d, i,j %d, %d, table size %d\n", index[1],i,j+1,nx*ny);
- if(index[2] >= nx*ny || index[2] < 0) fprintf(stderr,"index 2 out of bounds %d, i,j %d, %d, table size %d\n", index[2],i+1,j,nx*ny);
- if(index[3] >= nx*ny || index[3] < 0) fprintf(stderr,"index 3 out of bounds %d, i,j %d, %d, table size %d\n", index[3],i+1,j+1,nx*ny);
+ if (index[0] >= nx * ny || index[0] < 0)
+ fprintf(stderr, "index 0 out of bounds %d, i,j %d, %d, table size %d\n",
+ index[0], i, j, nx * ny);
+ if (index[1] >= nx * ny || index[1] < 0)
+ fprintf(stderr, "index 1 out of bounds %d, i,j %d, %d, table size %d\n",
+ index[1], i, j + 1, nx * ny);
+ if (index[2] >= nx * ny || index[2] < 0)
+ fprintf(stderr, "index 2 out of bounds %d, i,j %d, %d, table size %d\n",
+ index[2], i + 1, j, nx * ny);
+ if (index[3] >= nx * ny || index[3] < 0)
+ fprintf(stderr, "index 3 out of bounds %d, i,j %d, %d, table size %d\n",
+ index[3], i + 1, j + 1, nx * ny);
#endif
- result = (1 - dx) * (1 - dy) * table[index[0]] + (1 - dx) * dy * table[index[1]] +
- dx * (1 - dy) * table[index[2]] + dx * dy * table[index[3]];
+ result = (1 - dx) * (1 - dy) * table[index[0]] +
+ (1 - dx) * dy * table[index[1]] + dx * (1 - dy) * table[index[2]] +
+ dx * dy * table[index[3]];
return result;
}
@@ -310,28 +340,55 @@ __attribute__((always_inline)) INLINE double interpol_2d_dbl(double *table, int
* @param dx,dy,dz Offset within cell
* @param nx,ny,nz Table dimensions
*/
-__attribute__((always_inline)) INLINE float interpol_3d(float *table, int i, int j, int k, float dx, float dy,
- float dz, int nx, int ny, int nz) {
+__attribute__((always_inline)) INLINE float interpol_3d(float *table, int i,
+ int j, int k, float dx,
+ float dy, float dz,
+ int nx, int ny,
+ int nz) {
float result;
int index[8];
- index[0] = row_major_index_3d(i,j,k,nx,ny,nz);
- index[1] = row_major_index_3d(i,j,k+1,nx,ny,nz);
- index[2] = row_major_index_3d(i,j+1,k,nx,ny,nz);
- index[3] = row_major_index_3d(i,j+1,k+1,nx,ny,nz);
- index[4] = row_major_index_3d(i+1,j,k,nx,ny,nz);
- index[5] = row_major_index_3d(i+1,j,k+1,nx,ny,nz);
- index[6] = row_major_index_3d(i+1,j+1,k,nx,ny,nz);
- index[7] = row_major_index_3d(i+1,j+1,k+1,nx,ny,nz);
+ index[0] = row_major_index_3d(i, j, k, nx, ny, nz);
+ index[1] = row_major_index_3d(i, j, k + 1, nx, ny, nz);
+ index[2] = row_major_index_3d(i, j + 1, k, nx, ny, nz);
+ index[3] = row_major_index_3d(i, j + 1, k + 1, nx, ny, nz);
+ index[4] = row_major_index_3d(i + 1, j, k, nx, ny, nz);
+ index[5] = row_major_index_3d(i + 1, j, k + 1, nx, ny, nz);
+ index[6] = row_major_index_3d(i + 1, j + 1, k, nx, ny, nz);
+ index[7] = row_major_index_3d(i + 1, j + 1, k + 1, nx, ny, nz);
#ifdef SWIFT_DEBUG_CHECKS
- if(index[0] >= nx*ny*nz || index[0] < 0) fprintf(stderr,"index 0 out of bounds %d, i,j,k %d, %d, %d, table size %d\n", index[0],i,j,k,nx*ny*nz);
- if(index[1] >= nx*ny*nz || index[1] < 0) fprintf(stderr,"index 1 out of bounds %d, i,j,k %d, %d, %d, table size %d\n", index[1],i,j,k+1,nx*ny*nz);
- if(index[2] >= nx*ny*nz || index[2] < 0) fprintf(stderr,"index 2 out of bounds %d, i,j,k %d, %d, %d, table size %d\n", index[2],i,j+1,k,nx*ny*nz);
- if(index[3] >= nx*ny*nz || index[3] < 0) fprintf(stderr,"index 3 out of bounds %d, i,j,k %d, %d, %d, table size %d\n", index[3],i,j+1,k+1,nx*ny*nz);
- if(index[4] >= nx*ny*nz || index[4] < 0) fprintf(stderr,"index 4 out of bounds %d, i,j,k %d, %d, %d, table size %d\n", index[4],i+1,j,k,nx*ny*nz);
- if(index[5] >= nx*ny*nz || index[5] < 0) fprintf(stderr,"index 5 out of bounds %d, i,j,k %d, %d, %d, table size %d\n", index[5],i+1,j,k+1,nx*ny*nz);
- if(index[6] >= nx*ny*nz || index[6] < 0) fprintf(stderr,"index 6 out of bounds %d, i,j,k %d, %d, %d, table size %d\n", index[6],i+1,j+1,k,nx*ny*nz);
- if(index[7] >= nx*ny*nz || index[7] < 0) fprintf(stderr,"index 7 out of bounds %d, i,j,k %d, %d, %d, table size %d\n", index[7],i+1,j+1,k+1,nx*ny*nz);
+ if (index[0] >= nx * ny * nz || index[0] < 0)
+ fprintf(stderr,
+ "index 0 out of bounds %d, i,j,k %d, %d, %d, table size %d\n",
+ index[0], i, j, k, nx * ny * nz);
+ if (index[1] >= nx * ny * nz || index[1] < 0)
+ fprintf(stderr,
+ "index 1 out of bounds %d, i,j,k %d, %d, %d, table size %d\n",
+ index[1], i, j, k + 1, nx * ny * nz);
+ if (index[2] >= nx * ny * nz || index[2] < 0)
+ fprintf(stderr,
+ "index 2 out of bounds %d, i,j,k %d, %d, %d, table size %d\n",
+ index[2], i, j + 1, k, nx * ny * nz);
+ if (index[3] >= nx * ny * nz || index[3] < 0)
+ fprintf(stderr,
+ "index 3 out of bounds %d, i,j,k %d, %d, %d, table size %d\n",
+ index[3], i, j + 1, k + 1, nx * ny * nz);
+ if (index[4] >= nx * ny * nz || index[4] < 0)
+ fprintf(stderr,
+ "index 4 out of bounds %d, i,j,k %d, %d, %d, table size %d\n",
+ index[4], i + 1, j, k, nx * ny * nz);
+ if (index[5] >= nx * ny * nz || index[5] < 0)
+ fprintf(stderr,
+ "index 5 out of bounds %d, i,j,k %d, %d, %d, table size %d\n",
+ index[5], i + 1, j, k + 1, nx * ny * nz);
+ if (index[6] >= nx * ny * nz || index[6] < 0)
+ fprintf(stderr,
+ "index 6 out of bounds %d, i,j,k %d, %d, %d, table size %d\n",
+ index[6], i + 1, j + 1, k, nx * ny * nz);
+ if (index[7] >= nx * ny * nz || index[7] < 0)
+ fprintf(stderr,
+ "index 7 out of bounds %d, i,j,k %d, %d, %d, table size %d\n",
+ index[7], i + 1, j + 1, k + 1, nx * ny * nz);
#endif
result = (1 - dx) * (1 - dy) * (1 - dz) * table[index[0]] +
@@ -346,7 +403,6 @@ __attribute__((always_inline)) INLINE float interpol_3d(float *table, int i, int
return result;
}
-
/*
* @brief Performs 4d interpolation
*
@@ -355,44 +411,109 @@ __attribute__((always_inline)) INLINE float interpol_3d(float *table, int i, int
* @param dx,dy,dz,dw Offset within cell
* @param nx,ny,nz,nw Table dimensions
*/
-__attribute__((always_inline)) INLINE float interpol_4d(float *table, int i, int j, int k, int l, float dx,
- float dy, float dz, float dw, int nx, int ny, int nz, int nw) {
+__attribute__((always_inline)) INLINE float interpol_4d(
+ float *table, int i, int j, int k, int l, float dx, float dy, float dz,
+ float dw, int nx, int ny, int nz, int nw) {
float result;
int index[16];
- index[0] = row_major_index_4d(i,j,k,l,nx,ny,nz,nw);
- index[1] = row_major_index_4d(i,j,k,l+1,nx,ny,nz,nw);
- index[2] = row_major_index_4d(i,j,k+1,l,nx,ny,nz,nw);
- index[3] = row_major_index_4d(i,j,k+1,l+1,nx,ny,nz,nw);
- index[4] = row_major_index_4d(i,j+1,k,l,nx,ny,nz,nw);
- index[5] = row_major_index_4d(i,j+1,k,l+1,nx,ny,nz,nw);
- index[6] = row_major_index_4d(i,j+1,k+1,l,nx,ny,nz,nw);
- index[7] = row_major_index_4d(i,j+1,k+1,l+1,nx,ny,nz,nw);
- index[8] = row_major_index_4d(i+1,j,k,l,nx,ny,nz,nw);
- index[9] = row_major_index_4d(i+1,j,k,l+1,nx,ny,nz,nw);
- index[10] = row_major_index_4d(i+1,j,k+1,l,nx,ny,nz,nw);
- index[11] = row_major_index_4d(i+1,j,k+1,l+1,nx,ny,nz,nw);
- index[12] = row_major_index_4d(i+1,j+1,k,l,nx,ny,nz,nw);
- index[13] = row_major_index_4d(i+1,j+1,k,l+1,nx,ny,nz,nw);
- index[14] = row_major_index_4d(i+1,j+1,k+1,l,nx,ny,nz,nw);
- index[15] = row_major_index_4d(i+1,j+1,k+1,l+1,nx,ny,nz,nw);
+ index[0] = row_major_index_4d(i, j, k, l, nx, ny, nz, nw);
+ index[1] = row_major_index_4d(i, j, k, l + 1, nx, ny, nz, nw);
+ index[2] = row_major_index_4d(i, j, k + 1, l, nx, ny, nz, nw);
+ index[3] = row_major_index_4d(i, j, k + 1, l + 1, nx, ny, nz, nw);
+ index[4] = row_major_index_4d(i, j + 1, k, l, nx, ny, nz, nw);
+ index[5] = row_major_index_4d(i, j + 1, k, l + 1, nx, ny, nz, nw);
+ index[6] = row_major_index_4d(i, j + 1, k + 1, l, nx, ny, nz, nw);
+ index[7] = row_major_index_4d(i, j + 1, k + 1, l + 1, nx, ny, nz, nw);
+ index[8] = row_major_index_4d(i + 1, j, k, l, nx, ny, nz, nw);
+ index[9] = row_major_index_4d(i + 1, j, k, l + 1, nx, ny, nz, nw);
+ index[10] = row_major_index_4d(i + 1, j, k + 1, l, nx, ny, nz, nw);
+ index[11] = row_major_index_4d(i + 1, j, k + 1, l + 1, nx, ny, nz, nw);
+ index[12] = row_major_index_4d(i + 1, j + 1, k, l, nx, ny, nz, nw);
+ index[13] = row_major_index_4d(i + 1, j + 1, k, l + 1, nx, ny, nz, nw);
+ index[14] = row_major_index_4d(i + 1, j + 1, k + 1, l, nx, ny, nz, nw);
+ index[15] = row_major_index_4d(i + 1, j + 1, k + 1, l + 1, nx, ny, nz, nw);
#ifdef SWIFT_DEBUG_CHECKS
- if(index[0] >= nx*ny*nz*nw || index[0] < 0) fprintf(stderr,"index 0 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n", index[0],i,j,k,l,nx*ny*nz*nw);
- if(index[1] >= nx*ny*nz*nw || index[1] < 0) fprintf(stderr,"index 1 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n", index[1],i,j,k,l+1,nx*ny*nz*nw);
- if(index[2] >= nx*ny*nz*nw || index[2] < 0) fprintf(stderr,"index 2 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n", index[2],i,j,k+1,l,nx*ny*nz*nw);
- if(index[3] >= nx*ny*nz*nw || index[3] < 0) fprintf(stderr,"index 3 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n", index[3],i,j,k+1,l+1,nx*ny*nz*nw);
- if(index[4] >= nx*ny*nz*nw || index[4] < 0) fprintf(stderr,"index 4 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n", index[4],i,j+1,k,l,nx*ny*nz*nw);
- if(index[5] >= nx*ny*nz*nw || index[5] < 0) fprintf(stderr,"index 5 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n", index[5],i,j+1,k,l+1,nx*ny*nz*nw);
- if(index[6] >= nx*ny*nz*nw || index[6] < 0) fprintf(stderr,"index 6 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n", index[6],i,j+1,k+1,l,nx*ny*nz*nw);
- if(index[7] >= nx*ny*nz*nw || index[7] < 0) fprintf(stderr,"index 7 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n", index[7],i,j+1,k+1,l+1,nx*ny*nz*nw);
- if(index[8] >= nx*ny*nz*nw || index[8] < 0) fprintf(stderr,"index 8 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n", index[8],i+1,j,k,l,nx*ny*nz*nw);
- if(index[9] >= nx*ny*nz*nw || index[9] < 0) fprintf(stderr,"index 9 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n", index[9],i+1,j,k,l+1,nx*ny*nz*nw);
- if(index[10] >= nx*ny*nz*nw || index[10] < 0) fprintf(stderr,"index 10 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n", index[10],i+1,j,k+1,l,nx*ny*nz*nw);
- if(index[11] >= nx*ny*nz*nw || index[11] < 0) fprintf(stderr,"index 11 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n", index[11],i+1,j,k+1,l+1,nx*ny*nz*nw);
- if(index[12] >= nx*ny*nz*nw || index[12] < 0) fprintf(stderr,"index 12 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n", index[12],i+1,j+1,k,l,nx*ny*nz*nw);
- if(index[13] >= nx*ny*nz*nw || index[13] < 0) fprintf(stderr,"index 13 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n", index[13],i+1,j+1,k,l+1,nx*ny*nz*nw);
- if(index[14] >= nx*ny*nz*nw || index[14] < 0) fprintf(stderr,"index 14 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n", index[14],i+1,j+1,k+1,l,nx*ny*nz*nw);
- if(index[15] >= nx*ny*nz*nw || index[15] < 0) fprintf(stderr,"index 15 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n", index[15],i+1,j+1,k+1,l+1,nx*ny*nz*nw);
+ if (index[0] >= nx * ny * nz * nw || index[0] < 0)
+ fprintf(
+ stderr,
+ "index 0 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
+ index[0], i, j, k, l, nx * ny * nz * nw);
+ if (index[1] >= nx * ny * nz * nw || index[1] < 0)
+ fprintf(
+ stderr,
+ "index 1 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
+ index[1], i, j, k, l + 1, nx * ny * nz * nw);
+ if (index[2] >= nx * ny * nz * nw || index[2] < 0)
+ fprintf(
+ stderr,
+ "index 2 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
+ index[2], i, j, k + 1, l, nx * ny * nz * nw);
+ if (index[3] >= nx * ny * nz * nw || index[3] < 0)
+ fprintf(
+ stderr,
+ "index 3 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
+ index[3], i, j, k + 1, l + 1, nx * ny * nz * nw);
+ if (index[4] >= nx * ny * nz * nw || index[4] < 0)
+ fprintf(
+ stderr,
+ "index 4 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
+ index[4], i, j + 1, k, l, nx * ny * nz * nw);
+ if (index[5] >= nx * ny * nz * nw || index[5] < 0)
+ fprintf(
+ stderr,
+ "index 5 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
+ index[5], i, j + 1, k, l + 1, nx * ny * nz * nw);
+ if (index[6] >= nx * ny * nz * nw || index[6] < 0)
+ fprintf(
+ stderr,
+ "index 6 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
+ index[6], i, j + 1, k + 1, l, nx * ny * nz * nw);
+ if (index[7] >= nx * ny * nz * nw || index[7] < 0)
+ fprintf(
+ stderr,
+ "index 7 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
+ index[7], i, j + 1, k + 1, l + 1, nx * ny * nz * nw);
+ if (index[8] >= nx * ny * nz * nw || index[8] < 0)
+ fprintf(
+ stderr,
+ "index 8 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
+ index[8], i + 1, j, k, l, nx * ny * nz * nw);
+ if (index[9] >= nx * ny * nz * nw || index[9] < 0)
+ fprintf(
+ stderr,
+ "index 9 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
+ index[9], i + 1, j, k, l + 1, nx * ny * nz * nw);
+ if (index[10] >= nx * ny * nz * nw || index[10] < 0)
+ fprintf(
+ stderr,
+ "index 10 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
+ index[10], i + 1, j, k + 1, l, nx * ny * nz * nw);
+ if (index[11] >= nx * ny * nz * nw || index[11] < 0)
+ fprintf(
+ stderr,
+ "index 11 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
+ index[11], i + 1, j, k + 1, l + 1, nx * ny * nz * nw);
+ if (index[12] >= nx * ny * nz * nw || index[12] < 0)
+ fprintf(
+ stderr,
+ "index 12 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
+ index[12], i + 1, j + 1, k, l, nx * ny * nz * nw);
+ if (index[13] >= nx * ny * nz * nw || index[13] < 0)
+ fprintf(
+ stderr,
+ "index 13 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
+ index[13], i + 1, j + 1, k, l + 1, nx * ny * nz * nw);
+ if (index[14] >= nx * ny * nz * nw || index[14] < 0)
+ fprintf(
+ stderr,
+ "index 14 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
+ index[14], i + 1, j + 1, k + 1, l, nx * ny * nz * nw);
+ if (index[15] >= nx * ny * nz * nw || index[15] < 0)
+ fprintf(
+ stderr,
+ "index 15 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
+ index[15], i + 1, j + 1, k + 1, l + 1, nx * ny * nz * nw);
#endif
result = (1 - dx) * (1 - dy) * (1 - dz) * (1 - dw) * table[index[0]] +
@@ -416,7 +537,7 @@ __attribute__((always_inline)) INLINE float interpol_4d(float *table, int i, int
}
/*
- * @brief Interpolates 3d EAGLE table in redshift and hydrogen
+ * @brief Interpolates 3d EAGLE table in redshift and hydrogen
* number density. Produces 1d table depending only
* on log(temperature)
*
@@ -431,37 +552,37 @@ __attribute__((always_inline)) INLINE float interpol_4d(float *table, int i, int
* @param d_n_h Hydrogen number density offset
* @param result_table Pointer to 1d interpolated table
*/
-__attribute__((always_inline)) INLINE static void construct_1d_table_from_3d(const struct part* restrict p,
- const struct cooling_function_data* restrict cooling,
- const struct cosmology* restrict cosmo,
- const struct phys_const *internal_const,
- float *table,
- int z_i,
- float d_z,
- int n_h_i,
- float d_n_h,
- double *result_table){
+__attribute__((always_inline)) INLINE static void construct_1d_table_from_3d(
+ const struct part *restrict p,
+ const struct cooling_function_data *restrict cooling,
+ const struct cosmology *restrict cosmo,
+ const struct phys_const *internal_const, float *table, int z_i, float d_z,
+ int n_h_i, float d_n_h, double *result_table) {
int index[4];
- for(int i = 0; i < cooling->N_Temp; i++){
- index[0] = row_major_index_3d(z_i, n_h_i, i, cooling->N_Redshifts, cooling->N_nH, cooling->N_Temp);
- index[1] = row_major_index_3d(z_i, n_h_i+1, i, cooling->N_Redshifts, cooling->N_nH, cooling->N_Temp);
- index[2] = row_major_index_3d(z_i+1, n_h_i, i, cooling->N_Redshifts, cooling->N_nH, cooling->N_Temp);
- index[3] = row_major_index_3d(z_i+1, n_h_i+1, i, cooling->N_Redshifts, cooling->N_nH, cooling->N_Temp);
-
+ for (int i = 0; i < cooling->N_Temp; i++) {
+ index[0] = row_major_index_3d(z_i, n_h_i, i, cooling->N_Redshifts,
+ cooling->N_nH, cooling->N_Temp);
+ index[1] = row_major_index_3d(z_i, n_h_i + 1, i, cooling->N_Redshifts,
+ cooling->N_nH, cooling->N_Temp);
+ index[2] = row_major_index_3d(z_i + 1, n_h_i, i, cooling->N_Redshifts,
+ cooling->N_nH, cooling->N_Temp);
+ index[3] = row_major_index_3d(z_i + 1, n_h_i + 1, i, cooling->N_Redshifts,
+ cooling->N_nH, cooling->N_Temp);
+
result_table[i] = (1 - d_z) * (1 - d_n_h) * table[index[0]] +
- (1 - d_z) * d_n_h * table[index[1]] +
- d_z * (1 - d_n_h) * table[index[2]] +
- d_z * d_n_h * table[index[3]];
+ (1 - d_z) * d_n_h * table[index[1]] +
+ d_z * (1 - d_n_h) * table[index[2]] +
+ d_z * d_n_h * table[index[3]];
}
}
/*
* @brief Interpolates 4d EAGLE hydrogen and helium cooling
- * table in redshift, helium fraction and hydrogen number
+ * table in redshift, helium fraction and hydrogen number
* density. Produces 1d table depending only on log(temperature)
- * and location of maximum cooling gradient with respect to
- * internal energy. NOTE: calculation of location of maximum
+ * and location of maximum cooling gradient with respect to
+ * internal energy. NOTE: calculation of location of maximum
* cooling gradient is still work in progress.
*
* @param p Particle structure
@@ -480,72 +601,92 @@ __attribute__((always_inline)) INLINE static void construct_1d_table_from_3d(con
* @param u_ini Internal energy at start of hydro timestep
* @param dt Hydro timestep
*/
-__attribute__((always_inline)) INLINE static void construct_1d_table_from_4d_H_He(const struct part* restrict p,
- const struct cooling_function_data* restrict cooling,
- const struct cosmology* restrict cosmo,
- const struct phys_const *internal_const,
- float *table,
- int z_i,
- float d_z,
- int He_i,
- float d_He,
- int n_h_i,
- float d_n_h,
- double *result_table,
- double *x_max_grad,
- double u_ini,
- float dt){
+__attribute__((always_inline)) INLINE static void
+construct_1d_table_from_4d_H_He(
+ const struct part *restrict p,
+ const struct cooling_function_data *restrict cooling,
+ const struct cosmology *restrict cosmo,
+ const struct phys_const *internal_const, float *table, int z_i, float d_z,
+ int He_i, float d_He, int n_h_i, float d_n_h, double *result_table,
+ double *x_max_grad, double u_ini, float dt) {
int index[8];
float x0, x1, y0, y1, old_grad, new_grad;
- float inn_h = chemistry_get_number_density(p,cosmo,chemistry_element_H,internal_const)*cooling->number_density_scale;
+ float inn_h = chemistry_get_number_density(p, cosmo, chemistry_element_H,
+ internal_const) *
+ cooling->number_density_scale;
float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
float ratefact = inn_h * (XH / eagle_proton_mass_cgs);
old_grad = 0.0;
- for(int i = 0; i < cooling->N_Temp; i++){
- index[0] = row_major_index_4d(z_i, n_h_i, He_i, i, cooling->N_Redshifts, cooling->N_nH, cooling->N_He, cooling->N_Temp);
- index[1] = row_major_index_4d(z_i, n_h_i, He_i+1, i, cooling->N_Redshifts, cooling->N_nH, cooling->N_He, cooling->N_Temp);
- index[2] = row_major_index_4d(z_i, n_h_i+1, He_i, i, cooling->N_Redshifts, cooling->N_nH, cooling->N_He, cooling->N_Temp);
- index[3] = row_major_index_4d(z_i, n_h_i+1, He_i+1, i, cooling->N_Redshifts, cooling->N_nH, cooling->N_He, cooling->N_Temp);
- index[4] = row_major_index_4d(z_i+1, n_h_i, He_i, i, cooling->N_Redshifts, cooling->N_nH, cooling->N_He, cooling->N_Temp);
- index[5] = row_major_index_4d(z_i+1, n_h_i, He_i+1, i, cooling->N_Redshifts, cooling->N_nH, cooling->N_He, cooling->N_Temp);
- index[6] = row_major_index_4d(z_i+1, n_h_i+1, He_i, i, cooling->N_Redshifts, cooling->N_nH, cooling->N_He, cooling->N_Temp);
- index[7] = row_major_index_4d(z_i+1, n_h_i+1, He_i+1, i, cooling->N_Redshifts, cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ for (int i = 0; i < cooling->N_Temp; i++) {
+ index[0] =
+ row_major_index_4d(z_i, n_h_i, He_i, i, cooling->N_Redshifts,
+ cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ index[1] =
+ row_major_index_4d(z_i, n_h_i, He_i + 1, i, cooling->N_Redshifts,
+ cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ index[2] =
+ row_major_index_4d(z_i, n_h_i + 1, He_i, i, cooling->N_Redshifts,
+ cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ index[3] =
+ row_major_index_4d(z_i, n_h_i + 1, He_i + 1, i, cooling->N_Redshifts,
+ cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ index[4] =
+ row_major_index_4d(z_i + 1, n_h_i, He_i, i, cooling->N_Redshifts,
+ cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ index[5] =
+ row_major_index_4d(z_i + 1, n_h_i, He_i + 1, i, cooling->N_Redshifts,
+ cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ index[6] =
+ row_major_index_4d(z_i + 1, n_h_i + 1, He_i, i, cooling->N_Redshifts,
+ cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ index[7] = row_major_index_4d(z_i + 1, n_h_i + 1, He_i + 1, i,
+ cooling->N_Redshifts, cooling->N_nH,
+ cooling->N_He, cooling->N_Temp);
result_table[i] = (1 - d_z) * (1 - d_n_h) * (1 - d_He) * table[index[0]] +
- (1 - d_z) * (1 - d_n_h) * d_He * table[index[1]] +
- (1 - d_z) * d_n_h * (1 - d_He) * table[index[2]] +
- (1 - d_z) * d_n_h * d_He * table[index[3]] +
- d_z * (1 - d_n_h) * (1 - d_He) * table[index[4]] +
- d_z * (1 - d_n_h) * d_He * table[index[5]] +
- d_z * d_n_h * (1 - d_He) * table[index[6]] +
- d_z * d_n_h * d_He * table[index[7]];
+ (1 - d_z) * (1 - d_n_h) * d_He * table[index[1]] +
+ (1 - d_z) * d_n_h * (1 - d_He) * table[index[2]] +
+ (1 - d_z) * d_n_h * d_He * table[index[3]] +
+ d_z * (1 - d_n_h) * (1 - d_He) * table[index[4]] +
+ d_z * (1 - d_n_h) * d_He * table[index[5]] +
+ d_z * d_n_h * (1 - d_He) * table[index[6]] +
+ d_z * d_n_h * d_He * table[index[7]];
#if SWIFT_DEBUG_CHECKS
- if (isnan(result_table[i])) printf("Eagle cooling.h 1 i dz dnh dHe table values %d %.5e %.5e %.5e %.5e %.5e %.5e %.5e %.5e %.5e %.5e %.5e \n",i,d_z,d_n_h,d_He,table[index[0]],table[index[1]],table[index[2]],table[index[3]],table[index[4]],table[index[5]],table[index[6]],table[index[7]]);
+ if (isnan(result_table[i]))
+ printf(
+ "Eagle cooling.h 1 i dz dnh dHe table values %d %.5e %.5e %.5e %.5e "
+ "%.5e %.5e %.5e %.5e %.5e %.5e %.5e \n",
+ i, d_z, d_n_h, d_He, table[index[0]], table[index[1]],
+ table[index[2]], table[index[3]], table[index[4]], table[index[5]],
+ table[index[6]], table[index[7]]);
#endif
- if (i > 0 && dt > 0 && x_max_grad != NULL){
- x0 = pow(10.0,cooling->Therm[i-1]);
- x1 = pow(10.0,cooling->Therm[i]);
- y0 = x0 - u_ini - result_table[i-1]*ratefact*dt;
- y1 = x1 - u_ini - result_table[i]*ratefact*dt;
- //new_grad = (y1 - y0)/(x1 - x0);
- new_grad = log10(y1/y0)/log10(x1/x0);
- //printf("Eagle cooling.h gradient max gradient y1 y0 x1 x0 u_ini lambda %.5e %.5e %.5e %.5e %.5e %.5e %.5e %.5e \n", new_grad, old_grad, y1, y0, x1, x0, u_ini, result_table[i-1]*ratefact*dt);
+ if (i > 0 && dt > 0 && x_max_grad != NULL) {
+ x0 = pow(10.0, cooling->Therm[i - 1]);
+ x1 = pow(10.0, cooling->Therm[i]);
+ y0 = x0 - u_ini - result_table[i - 1] * ratefact * dt;
+ y1 = x1 - u_ini - result_table[i] * ratefact * dt;
+ // new_grad = (y1 - y0)/(x1 - x0);
+ new_grad = log10(y1 / y0) / log10(x1 / x0);
+ // printf("Eagle cooling.h gradient max gradient y1 y0 x1 x0 u_ini lambda
+ // %.5e %.5e %.5e %.5e %.5e %.5e %.5e %.5e \n", new_grad, old_grad, y1,
+ // y0, x1, x0, u_ini, result_table[i-1]*ratefact*dt);
if (new_grad > old_grad) {
- *x_max_grad = 0.5*(x0 + x1);
+ *x_max_grad = 0.5 * (x0 + x1);
old_grad = new_grad;
}
}
}
- //printf("Eagle cooling.h max gradient, guess, u_ini %.5e %.5e %.5e\n", old_grad, *x_max_grad, u_ini);
+ // printf("Eagle cooling.h max gradient, guess, u_ini %.5e %.5e %.5e\n",
+ // old_grad, *x_max_grad, u_ini);
}
/*
- * @brief Interpolates 4d EAGLE table in redshift,
- * helium fraction and hydrogen number density.
+ * @brief Interpolates 4d EAGLE table in redshift,
+ * helium fraction and hydrogen number density.
* Produces 1d table depending only on log(temperature)
- *
+ *
* @param p Particle structure
* @param cooling Cooling data structure
* @param cosmo Cosmology structure
@@ -559,50 +700,66 @@ __attribute__((always_inline)) INLINE static void construct_1d_table_from_4d_H_H
* @param d_n_h Hydrogen number density offset
* @param result_table Pointer to 1d interpolated table
*/
-__attribute__((always_inline)) INLINE static void construct_1d_table_from_4d(const struct part* restrict p,
- const struct cooling_function_data* restrict cooling,
- const struct cosmology* restrict cosmo,
- const struct phys_const *internal_const,
- float *table,
- int z_i,
- float d_z,
- int He_i,
- float d_He,
- int n_h_i,
- float d_n_h,
- double *result_table){
+__attribute__((always_inline)) INLINE static void construct_1d_table_from_4d(
+ const struct part *restrict p,
+ const struct cooling_function_data *restrict cooling,
+ const struct cosmology *restrict cosmo,
+ const struct phys_const *internal_const, float *table, int z_i, float d_z,
+ int He_i, float d_He, int n_h_i, float d_n_h, double *result_table) {
int index[8];
- for(int i = 0; i < cooling->N_Temp; i++){
- index[0] = row_major_index_4d(z_i, n_h_i, He_i, i, cooling->N_Redshifts, cooling->N_nH, cooling->N_He, cooling->N_Temp);
- index[1] = row_major_index_4d(z_i, n_h_i, He_i+1, i, cooling->N_Redshifts, cooling->N_nH, cooling->N_He, cooling->N_Temp);
- index[2] = row_major_index_4d(z_i, n_h_i+1, He_i, i, cooling->N_Redshifts, cooling->N_nH, cooling->N_He, cooling->N_Temp);
- index[3] = row_major_index_4d(z_i, n_h_i+1, He_i+1, i, cooling->N_Redshifts, cooling->N_nH, cooling->N_He, cooling->N_Temp);
- index[4] = row_major_index_4d(z_i+1, n_h_i, He_i, i, cooling->N_Redshifts, cooling->N_nH, cooling->N_He, cooling->N_Temp);
- index[5] = row_major_index_4d(z_i+1, n_h_i, He_i+1, i, cooling->N_Redshifts, cooling->N_nH, cooling->N_He, cooling->N_Temp);
- index[6] = row_major_index_4d(z_i+1, n_h_i+1, He_i, i, cooling->N_Redshifts, cooling->N_nH, cooling->N_He, cooling->N_Temp);
- index[7] = row_major_index_4d(z_i+1, n_h_i+1, He_i+1, i, cooling->N_Redshifts, cooling->N_nH, cooling->N_He, cooling->N_Temp);
-
+ for (int i = 0; i < cooling->N_Temp; i++) {
+ index[0] =
+ row_major_index_4d(z_i, n_h_i, He_i, i, cooling->N_Redshifts,
+ cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ index[1] =
+ row_major_index_4d(z_i, n_h_i, He_i + 1, i, cooling->N_Redshifts,
+ cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ index[2] =
+ row_major_index_4d(z_i, n_h_i + 1, He_i, i, cooling->N_Redshifts,
+ cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ index[3] =
+ row_major_index_4d(z_i, n_h_i + 1, He_i + 1, i, cooling->N_Redshifts,
+ cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ index[4] =
+ row_major_index_4d(z_i + 1, n_h_i, He_i, i, cooling->N_Redshifts,
+ cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ index[5] =
+ row_major_index_4d(z_i + 1, n_h_i, He_i + 1, i, cooling->N_Redshifts,
+ cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ index[6] =
+ row_major_index_4d(z_i + 1, n_h_i + 1, He_i, i, cooling->N_Redshifts,
+ cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ index[7] = row_major_index_4d(z_i + 1, n_h_i + 1, He_i + 1, i,
+ cooling->N_Redshifts, cooling->N_nH,
+ cooling->N_He, cooling->N_Temp);
+
result_table[i] = (1 - d_z) * (1 - d_n_h) * (1 - d_He) * table[index[0]] +
- (1 - d_z) * (1 - d_n_h) * d_He * table[index[1]] +
- (1 - d_z) * d_n_h * (1 - d_He) * table[index[2]] +
- (1 - d_z) * d_n_h * d_He * table[index[3]] +
- d_z * (1 - d_n_h) * (1 - d_He) * table[index[4]] +
- d_z * (1 - d_n_h) * d_He * table[index[5]] +
- d_z * d_n_h * (1 - d_He) * table[index[6]] +
- d_z * d_n_h * d_He * table[index[7]];
+ (1 - d_z) * (1 - d_n_h) * d_He * table[index[1]] +
+ (1 - d_z) * d_n_h * (1 - d_He) * table[index[2]] +
+ (1 - d_z) * d_n_h * d_He * table[index[3]] +
+ d_z * (1 - d_n_h) * (1 - d_He) * table[index[4]] +
+ d_z * (1 - d_n_h) * d_He * table[index[5]] +
+ d_z * d_n_h * (1 - d_He) * table[index[6]] +
+ d_z * d_n_h * d_He * table[index[7]];
#if SWIFT_DEBUG_CHECKS
- if (isnan(result_table[i])) printf("Eagle cooling.h 2 i dz dnh dHe table values %d %.5e %.5e %.5e %.5e %.5e %.5e %.5e %.5e %.5e %.5e %.5e \n",i,d_z,d_n_h,d_He,table[index[0]],table[index[1]],table[index[2]],table[index[3]],table[index[4]],table[index[5]],table[index[6]],table[index[7]]);
+ if (isnan(result_table[i]))
+ printf(
+ "Eagle cooling.h 2 i dz dnh dHe table values %d %.5e %.5e %.5e %.5e "
+ "%.5e %.5e %.5e %.5e %.5e %.5e %.5e \n",
+ i, d_z, d_n_h, d_He, table[index[0]], table[index[1]],
+ table[index[2]], table[index[3]], table[index[4]], table[index[5]],
+ table[index[6]], table[index[7]]);
#endif
}
}
/*
- * @brief Interpolates 4d EAGLE metal cooling tables in redshift,
- * helium fraction and hydrogen number density.
+ * @brief Interpolates 4d EAGLE metal cooling tables in redshift,
+ * helium fraction and hydrogen number density.
* Produces 1d table depending only on log(temperature)
* NOTE: Will need to be modified to incorporate particle to solar
* electron abundance ratio.
- *
+ *
* @param p Particle structure
* @param cooling Cooling data structure
* @param cosmo Cosmology structure
@@ -616,35 +773,49 @@ __attribute__((always_inline)) INLINE static void construct_1d_table_from_4d(con
* @param d_n_h Hydrogen number density offset
* @param result_table Pointer to 1d interpolated table
*/
-__attribute__((always_inline)) INLINE static void construct_1d_table_from_4d_elements(const struct part* restrict p,
- const struct cooling_function_data* restrict cooling,
- const struct cosmology* restrict cosmo,
- const struct phys_const *internal_const,
- float *table,
- int z_i,
- float d_z,
- int n_h_i,
- float d_n_h,
- double *result_table){
+__attribute__((always_inline)) INLINE static void
+construct_1d_table_from_4d_elements(
+ const struct part *restrict p,
+ const struct cooling_function_data *restrict cooling,
+ const struct cosmology *restrict cosmo,
+ const struct phys_const *internal_const, float *table, int z_i, float d_z,
+ int n_h_i, float d_n_h, double *result_table) {
int index[4];
- for(int j = 0; j < cooling->N_Elements; j++){
- for(int i = 0; i < cooling->N_Temp; i++){
+ for (int j = 0; j < cooling->N_Elements; j++) {
+ for (int i = 0; i < cooling->N_Temp; i++) {
if (j == 0) result_table[i] = 0.0;
- index[0] = row_major_index_4d(z_i, j, n_h_i, i, cooling->N_Redshifts, cooling->N_Elements, cooling->N_nH, cooling->N_Temp);
- index[1] = row_major_index_4d(z_i, j, n_h_i+1, i, cooling->N_Redshifts, cooling->N_Elements, cooling->N_nH, cooling->N_Temp);
- index[2] = row_major_index_4d(z_i+1, j, n_h_i, i, cooling->N_Redshifts, cooling->N_Elements, cooling->N_nH, cooling->N_Temp);
- index[3] = row_major_index_4d(z_i+1, j, n_h_i+1, i, cooling->N_Redshifts, cooling->N_Elements, cooling->N_nH, cooling->N_Temp);
+ index[0] = row_major_index_4d(z_i, j, n_h_i, i, cooling->N_Redshifts,
+ cooling->N_Elements, cooling->N_nH,
+ cooling->N_Temp);
+ index[1] = row_major_index_4d(z_i, j, n_h_i + 1, i, cooling->N_Redshifts,
+ cooling->N_Elements, cooling->N_nH,
+ cooling->N_Temp);
+ index[2] = row_major_index_4d(z_i + 1, j, n_h_i, i, cooling->N_Redshifts,
+ cooling->N_Elements, cooling->N_nH,
+ cooling->N_Temp);
+ index[3] = row_major_index_4d(z_i + 1, j, n_h_i + 1, i,
+ cooling->N_Redshifts, cooling->N_Elements,
+ cooling->N_nH, cooling->N_Temp);
result_table[i] += (1 - d_z) * (1 - d_n_h) * table[index[0]] +
- (1 - d_z) * d_n_h * table[index[1]] +
- d_z * (1 - d_n_h) * table[index[2]] +
- d_z * d_n_h * table[index[3]];
+ (1 - d_z) * d_n_h * table[index[1]] +
+ d_z * (1 - d_n_h) * table[index[2]] +
+ d_z * d_n_h * table[index[3]];
#if SWIFT_DEBUG_CHECKS
- if (isnan(result_table[i])) printf("Eagle cooling.h 3 i partial sums %d %.5e %.5e %.5e %.5e %.5e \n",i, result_table[i],(1 - d_z) * (1 - d_n_h) * table[index[0]],
- (1 - d_z) * (1 - d_n_h) * table[index[0]] + (1 - d_z) * d_n_h * table[index[1]],
- (1 - d_z) * (1 - d_n_h) * table[index[0]] + (1 - d_z) * d_n_h * table[index[1]] + d_z * (1 - d_n_h) * table[index[2]],
- (1 - d_z) * (1 - d_n_h) * table[index[0]] + (1 - d_z) * d_n_h * table[index[1]] + d_z * (1 - d_n_h) * table[index[2]] + d_z * d_n_h * table[index[3]]);
+ if (isnan(result_table[i]))
+ printf(
+ "Eagle cooling.h 3 i partial sums %d %.5e %.5e %.5e %.5e %.5e \n",
+ i, result_table[i], (1 - d_z) * (1 - d_n_h) * table[index[0]],
+ (1 - d_z) * (1 - d_n_h) * table[index[0]] +
+ (1 - d_z) * d_n_h * table[index[1]],
+ (1 - d_z) * (1 - d_n_h) * table[index[0]] +
+ (1 - d_z) * d_n_h * table[index[1]] +
+ d_z * (1 - d_n_h) * table[index[2]],
+ (1 - d_z) * (1 - d_n_h) * table[index[0]] +
+ (1 - d_z) * d_n_h * table[index[1]] +
+ d_z * (1 - d_n_h) * table[index[2]] +
+ d_z * d_n_h * table[index[3]]);
#endif
}
}
@@ -659,65 +830,71 @@ __attribute__((always_inline)) INLINE static void construct_1d_table_from_4d_ele
* @param cooling Cooling data structure
* @param p Particle structure
*/
-inline int set_cooling_SolarAbundances(const float *element_abundance,
- double *cooling_element_abundance,
- const struct cooling_function_data* restrict cooling,
- const struct part* restrict p) {
+inline int set_cooling_SolarAbundances(
+ const float *element_abundance, double *cooling_element_abundance,
+ const struct cooling_function_data *restrict cooling,
+ const struct part *restrict p) {
int i, index;
int Silicon_SPH_Index = -1;
int Calcium_SPH_Index = -1;
int Sulphur_SPH_Index = -1;
-
+
int Silicon_CoolHeat_Index = -1;
int Calcium_CoolHeat_Index = -1;
int Sulphur_CoolHeat_Index = -1;
- //float *cooling->ElementAbundance_SOLARM1 = malloc(cooling->N_SolarAbundances*sizeof(float));
+ // float *cooling->ElementAbundance_SOLARM1 =
+ // malloc(cooling->N_SolarAbundances*sizeof(float));
- /* determine (inverse of) solar abundance of these elements */
- for (i = 0; i < cooling->N_Elements; i++) {
- index =
- get_element_index(cooling->SolarAbundanceNames,
- cooling->N_SolarAbundances, cooling->ElementNames[i]);
-
- if (index < 0) error("Eagle cooling.h index out of bounds");
+ /* determine (inverse of) solar abundance of these elements */
+ for (i = 0; i < cooling->N_Elements; i++) {
+ index =
+ get_element_index(cooling->SolarAbundanceNames,
+ cooling->N_SolarAbundances, cooling->ElementNames[i]);
- index = cooling->SolarAbundanceNamePointers[i];
+ if (index < 0) error("Eagle cooling.h index out of bounds");
- if(cooling->SolarAbundances[index] != 0) cooling->ElementAbundance_SOLARM1[i] = 1. / cooling->SolarAbundances[index];
- else cooling->ElementAbundance_SOLARM1[i] = 0.0;
+ index = cooling->SolarAbundanceNamePointers[i];
- }
+ if (cooling->SolarAbundances[index] != 0)
+ cooling->ElementAbundance_SOLARM1[i] =
+ 1. / cooling->SolarAbundances[index];
+ else
+ cooling->ElementAbundance_SOLARM1[i] = 0.0;
+ }
- /* Sulphur tracks Silicon: may choose not to follow Sulphur as SPH element
- */
- /* Same is true for Calcium */
- /* We will assume the code tracks Silicon, and may need to scale Calcium and
- * Sulphur accordingly */
-
- Silicon_SPH_Index = element_index("Silicon",cooling);
- Calcium_SPH_Index = element_index("Calcium",cooling);
- Sulphur_SPH_Index = element_index("Sulphur",cooling);
-
- Silicon_CoolHeat_Index =
- get_element_index(cooling->ElementNames, cooling->N_Elements, "Silicon");
- Calcium_CoolHeat_Index =
- get_element_index(cooling->ElementNames, cooling->N_Elements, "Calcium");
- Sulphur_CoolHeat_Index =
- get_element_index(cooling->ElementNames, cooling->N_Elements, "Sulphur");
-
- if (Silicon_CoolHeat_Index == -1 || Calcium_CoolHeat_Index == -1 ||
- Sulphur_CoolHeat_Index == -1) {
- error("Si, Ca, or S index out of bound\n");
- }
+ /* Sulphur tracks Silicon: may choose not to follow Sulphur as SPH element
+ */
+ /* Same is true for Calcium */
+ /* We will assume the code tracks Silicon, and may need to scale Calcium and
+ * Sulphur accordingly */
+
+ Silicon_SPH_Index = element_index("Silicon", cooling);
+ Calcium_SPH_Index = element_index("Calcium", cooling);
+ Sulphur_SPH_Index = element_index("Sulphur", cooling);
+
+ Silicon_CoolHeat_Index =
+ get_element_index(cooling->ElementNames, cooling->N_Elements, "Silicon");
+ Calcium_CoolHeat_Index =
+ get_element_index(cooling->ElementNames, cooling->N_Elements, "Calcium");
+ Sulphur_CoolHeat_Index =
+ get_element_index(cooling->ElementNames, cooling->N_Elements, "Sulphur");
+
+ if (Silicon_CoolHeat_Index == -1 || Calcium_CoolHeat_Index == -1 ||
+ Sulphur_CoolHeat_Index == -1) {
+ error("Si, Ca, or S index out of bound\n");
+ }
int sili_index;
for (i = 0; i < cooling->N_Elements; i++) {
- if (strcmp(chemistry_get_element_name((enum chemistry_element) i), "Silicon") == 0) sili_index = i;
+ if (strcmp(chemistry_get_element_name((enum chemistry_element)i),
+ "Silicon") == 0)
+ sili_index = i;
}
- // Eagle way of identifying and assigning element abundance with strange workaround for calcium and sulphur
- //for (i = 0; i < cooling->N_Elements; i++) {
+ // Eagle way of identifying and assigning element abundance with strange
+ // workaround for calcium and sulphur
+ // for (i = 0; i < cooling->N_Elements; i++) {
// if (i == Calcium_CoolHeat_Index && Calcium_SPH_Index == -1)
// /* SPH does not track Calcium: use Si abundance */
// if (Silicon_SPH_Index == -1)
@@ -737,35 +914,47 @@ inline int set_cooling_SolarAbundances(const float *element_abundance,
// cooling_ElementAbundance_SOLARM1[Silicon_CoolHeat_Index];
// }
// else{
- // cooling_element_abundance[i] = element_abundance[cooling->ElementNamePointers[i]] *
+ // cooling_element_abundance[i] =
+ // element_abundance[cooling->ElementNamePointers[i]] *
// cooling_ElementAbundance_SOLARM1[i];
- // //printf ("Eagle cooling.h element, name, abundance, solar abundance, solarm1, cooling abundance %d %s %.5e %.5e %.5e %.5e\n",cooling->ElementNamePointers[i],cooling->ElementNames[i], element_abundance[cooling->ElementNamePointers[i]],cooling->SolarAbundances[i], cooling_ElementAbundance_SOLARM1[i], cooling_element_abundance[i]);
+ // //printf ("Eagle cooling.h element, name, abundance, solar abundance,
+ // solarm1, cooling abundance %d %s %.5e %.5e %.5e
+ // %.5e\n",cooling->ElementNamePointers[i],cooling->ElementNames[i],
+ // element_abundance[cooling->ElementNamePointers[i]],cooling->SolarAbundances[i],
+ // cooling_ElementAbundance_SOLARM1[i], cooling_element_abundance[i]);
// }
//}
-
+
for (i = 0; i < cooling->N_Elements; i++) {
if (i == Calcium_CoolHeat_Index && Calcium_SPH_Index != -1)
if (Silicon_SPH_Index == -1)
cooling_element_abundance[i] = 0.0;
- else{
+ else {
cooling_element_abundance[i] =
- element_abundance[sili_index] * cooling->calcium_over_silicon_ratio *
+ element_abundance[sili_index] *
+ cooling->calcium_over_silicon_ratio *
cooling->ElementAbundance_SOLARM1[Calcium_CoolHeat_Index];
}
else if (i == Sulphur_CoolHeat_Index && Sulphur_SPH_Index != -1)
/* SPH does not track Sulphur: use Si abundance */
if (Silicon_SPH_Index == -1)
cooling_element_abundance[i] = 0.0;
- else{
+ else {
cooling_element_abundance[i] =
- element_abundance[sili_index] * cooling->sulphur_over_silicon_ratio *
+ element_abundance[sili_index] *
+ cooling->sulphur_over_silicon_ratio *
cooling->ElementAbundance_SOLARM1[Sulphur_CoolHeat_Index];
}
- else{
- cooling_element_abundance[i] = element_abundance[cooling->ElementNamePointers[i]] *
- cooling->ElementAbundance_SOLARM1[i];
+ else {
+ cooling_element_abundance[i] =
+ element_abundance[cooling->ElementNamePointers[i]] *
+ cooling->ElementAbundance_SOLARM1[i];
}
- //printf ("Eagle cooling.h i, solar abundance name pointer, element, name, abundance, solarm1, cooling abundance %d %d %d %s %.5e %.5e %.5e\n", i, cooling->SolarAbundanceNamePointers[i],cooling->ElementNamePointers[i],cooling->ElementNames[i], element_abundance[cooling->ElementNamePointers[i]], cooling->ElementAbundance_SOLARM1[i], cooling_element_abundance[i]);
+ // printf ("Eagle cooling.h i, solar abundance name pointer, element, name,
+ // abundance, solarm1, cooling abundance %d %d %d %s %.5e %.5e %.5e\n", i,
+ // cooling->SolarAbundanceNamePointers[i],cooling->ElementNamePointers[i],cooling->ElementNames[i],
+ // element_abundance[cooling->ElementNamePointers[i]],
+ // cooling->ElementAbundance_SOLARM1[i], cooling_element_abundance[i]);
}
return 0;
@@ -778,48 +967,52 @@ inline int set_cooling_SolarAbundances(const float *element_abundance,
* @param temperature_table Table of EAGLE temperature values
* @param cooling Cooling data structure
*/
-__attribute__((always_inline)) INLINE static double eagle_convert_temp_to_u_1d_table(double temp,
- float *temperature_table,
- const struct cooling_function_data* restrict cooling) {
-
+__attribute__((always_inline)) INLINE static double
+eagle_convert_temp_to_u_1d_table(double temp, float *temperature_table,
+ const struct cooling_function_data *restrict
+ cooling) {
+
int temp_i;
float d_temp, u;
- get_index_1d(temperature_table, cooling->N_Temp, log10(temp), &temp_i, &d_temp);
+ get_index_1d(temperature_table, cooling->N_Temp, log10(temp), &temp_i,
+ &d_temp);
- u = pow(10.0,interpol_1d(cooling->Therm,temp_i,d_temp));
+ u = pow(10.0, interpol_1d(cooling->Therm, temp_i, d_temp));
return u;
}
/*
* @brief Interpolates temperature from internal energy based on table and
- * calculates the size of the internal energy cell for the specified temperature.
+ * calculates the size of the internal energy cell for the specified
+ * temperature.
*
- * @param u Internal energy
+ * @param u Internal energy
* @param delta_u Pointer to size of internal energy cell
* @param temperature_table Table of EAGLE temperature values
* @param cooling Cooling data structure
*/
-__attribute__((always_inline)) INLINE static double eagle_convert_u_to_temp_1d_table(double u,
- float *delta_u,
- double *temperature_table,
- const struct part* restrict p,
- const struct cooling_function_data* restrict cooling,
- const struct cosmology* restrict cosmo,
- const struct phys_const *internal_const) {
-
+__attribute__((always_inline)) INLINE static double
+eagle_convert_u_to_temp_1d_table(
+ double u, float *delta_u, double *temperature_table,
+ const struct part *restrict p,
+ const struct cooling_function_data *restrict cooling,
+ const struct cosmology *restrict cosmo,
+ const struct phys_const *internal_const) {
+
int u_i;
float d_u, logT, T;
get_index_1d(cooling->Therm, cooling->N_Temp, log10(u), &u_i, &d_u);
- logT = interpol_1d_dbl(temperature_table,u_i,d_u);
+ logT = interpol_1d_dbl(temperature_table, u_i, d_u);
T = pow(10.0, logT);
if (u_i == 0 && d_u == 0) T *= u / pow(10.0, cooling->Therm[0]);
- *delta_u = pow(10.0,cooling->Therm[u_i + 1]) - pow(10.0,cooling->Therm[u_i]);
+ *delta_u =
+ pow(10.0, cooling->Therm[u_i + 1]) - pow(10.0, cooling->Therm[u_i]);
return T;
}
@@ -839,24 +1032,19 @@ __attribute__((always_inline)) INLINE static double eagle_convert_u_to_temp_1d_t
* @param cooling Cooling data structure
* @param cosmo Cosmology structure
* @param internal_const Physical constants structure
- * @param element_lambda Pointer for printing contribution to cooling rate
+ * @param element_lambda Pointer for printing contribution to cooling rate
* from each of the metals.
*/
-__attribute__((always_inline)) INLINE static double eagle_metal_cooling_rate(double u,
- double *temp_table,
- int z_index,
- float dz,
- int n_h_i,
- float d_n_h,
- int He_i,
- float d_He,
- const struct part* restrict p,
- const struct cooling_function_data* restrict cooling,
- const struct cosmology* restrict cosmo,
- const struct phys_const *internal_const,
- double* element_lambda) {
+__attribute__((always_inline)) INLINE static double eagle_metal_cooling_rate(
+ double u, double *temp_table, int z_index, float dz, int n_h_i, float d_n_h,
+ int He_i, float d_He, const struct part *restrict p,
+ const struct cooling_function_data *restrict cooling,
+ const struct cosmology *restrict cosmo,
+ const struct phys_const *internal_const, double *element_lambda) {
double T_gam, solar_electron_abundance;
- double n_h = chemistry_get_number_density(p,cosmo,chemistry_element_H,internal_const)*cooling->number_density_scale; // chemistry_data
+ double n_h = chemistry_get_number_density(p, cosmo, chemistry_element_H,
+ internal_const) *
+ cooling->number_density_scale; // chemistry_data
double z = cosmo->z;
double cooling_rate = 0.0, temp_lambda;
float du;
@@ -870,7 +1058,8 @@ __attribute__((always_inline)) INLINE static double eagle_metal_cooling_rate(dou
#if SWIFT_DEBUG_CHECKS
if (isnan(u)) printf("u is nan id %llu\n", p->id);
#endif
- temp = eagle_convert_u_to_temp_1d_table(u,&du,temp_table,p,cooling,cosmo,internal_const);
+ temp = eagle_convert_u_to_temp_1d_table(u, &du, temp_table, p, cooling, cosmo,
+ internal_const);
get_index_1d(cooling->Temp, cooling->N_Temp, log10(temp), &temp_i, &d_temp);
@@ -878,27 +1067,38 @@ __attribute__((always_inline)) INLINE static double eagle_metal_cooling_rate(dou
/* Metal-free cooling */
/* ------------------ */
- /* Collisional cooling */
- //element_lambda[0] =
- // interpol_2d(cooling->table.collisional_cooling.H_plus_He_heating, He_i,
- // temp_i, d_He, d_temp,cooling->N_He,cooling->N_Temp);
- //h_plus_he_electron_abundance =
- // interpol_2d(cooling->table.collisional_cooling.H_plus_He_electron_abundance, He_i,
- // temp_i, d_He, d_temp,cooling->N_He,cooling->N_Temp);
-
- /* Photodissociation */
- //element_lambda[0] =
- // interpol_3d(cooling->table.photodissociation_cooling.H_plus_He_heating, He_i,
- // temp_i, n_h_i, d_He, d_temp, d_n_h,cooling->N_He,cooling->N_Temp,cooling->N_nH);
- //h_plus_he_electron_abundance =
- // interpol_3d(cooling->table.photodissociation_cooling.H_plus_He_electron_abundance, He_i,
- // temp_i, n_h_i, d_He, d_temp, d_n_h,cooling->N_He,cooling->N_Temp,cooling->N_nH);
-
- /* redshift tables */
- temp_lambda = interpol_4d(cooling->table.element_cooling.H_plus_He_heating, z_index, n_h_i, He_i, temp_i, dz, d_n_h, d_He, d_temp,cooling->N_Redshifts,cooling->N_nH,cooling->N_He,cooling->N_Temp);
- h_plus_he_electron_abundance = interpol_4d(cooling->table.element_cooling.H_plus_He_electron_abundance, z_index, n_h_i, He_i, temp_i, dz, d_n_h, d_He, d_temp,cooling->N_Redshifts,cooling->N_nH,cooling->N_He,cooling->N_Temp);
- cooling_rate += temp_lambda;
- if (element_lambda != NULL) element_lambda[0] = temp_lambda;
+ /* Collisional cooling */
+ // element_lambda[0] =
+ // interpol_2d(cooling->table.collisional_cooling.H_plus_He_heating, He_i,
+ // temp_i, d_He, d_temp,cooling->N_He,cooling->N_Temp);
+ // h_plus_he_electron_abundance =
+ // interpol_2d(cooling->table.collisional_cooling.H_plus_He_electron_abundance,
+ // He_i,
+ // temp_i, d_He, d_temp,cooling->N_He,cooling->N_Temp);
+
+ /* Photodissociation */
+ // element_lambda[0] =
+ // interpol_3d(cooling->table.photodissociation_cooling.H_plus_He_heating,
+ // He_i,
+ // temp_i, n_h_i, d_He, d_temp,
+ // d_n_h,cooling->N_He,cooling->N_Temp,cooling->N_nH);
+ // h_plus_he_electron_abundance =
+ // interpol_3d(cooling->table.photodissociation_cooling.H_plus_He_electron_abundance,
+ // He_i,
+ // temp_i, n_h_i, d_He, d_temp,
+ // d_n_h,cooling->N_He,cooling->N_Temp,cooling->N_nH);
+
+ /* redshift tables */
+ temp_lambda = interpol_4d(cooling->table.element_cooling.H_plus_He_heating,
+ z_index, n_h_i, He_i, temp_i, dz, d_n_h, d_He,
+ d_temp, cooling->N_Redshifts, cooling->N_nH,
+ cooling->N_He, cooling->N_Temp);
+ h_plus_he_electron_abundance = interpol_4d(
+ cooling->table.element_cooling.H_plus_He_electron_abundance, z_index,
+ n_h_i, He_i, temp_i, dz, d_n_h, d_He, d_temp, cooling->N_Redshifts,
+ cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ cooling_rate += temp_lambda;
+ if (element_lambda != NULL) element_lambda[0] = temp_lambda;
/* ------------------ */
/* Compton cooling */
@@ -911,8 +1111,9 @@ __attribute__((always_inline)) INLINE static double eagle_metal_cooling_rate(dou
T_gam = eagle_cmb_temperature * (1 + z);
- temp_lambda = -eagle_compton_rate * (temp - eagle_cmb_temperature * (1 + z)) * pow((1 + z), 4) *
- h_plus_he_electron_abundance / n_h;
+ temp_lambda = -eagle_compton_rate *
+ (temp - eagle_cmb_temperature * (1 + z)) * pow((1 + z), 4) *
+ h_plus_he_electron_abundance / n_h;
cooling_rate += temp_lambda;
if (element_lambda != NULL) element_lambda[1] = temp_lambda;
}
@@ -921,48 +1122,62 @@ __attribute__((always_inline)) INLINE static double eagle_metal_cooling_rate(dou
/* Metal cooling */
/* ------------- */
- /* for each element, find the abundance and multiply it
- by the interpolated heating-cooling */
+ /* for each element, find the abundance and multiply it
+ by the interpolated heating-cooling */
- //set_cooling_SolarAbundances(p->chemistry_data.metal_mass_fraction, cooling->solar_abundances, cooling, p);
+ // set_cooling_SolarAbundances(p->chemistry_data.metal_mass_fraction,
+ // cooling->solar_abundances, cooling, p);
- /* Collisional cooling */
- //solar_electron_abundance =
- // interpol_1d(cooling->table.collisional_cooling.electron_abundance, temp_i, d_temp); /* ne/n_h */
+ /* Collisional cooling */
+ // solar_electron_abundance =
+ // interpol_1d(cooling->table.collisional_cooling.electron_abundance,
+ // temp_i, d_temp); /* ne/n_h */
- //for (i = 0; i < cooling->N_Elements; i++){
- // element_lambda[i+2] = interpol_2d(cooling->table.collisional_cooling.metal_heating, i,
- // temp_i, 0.0, d_temp,cooling->N_Elements,cooling->N_Temp) *
- // (h_plus_he_electron_abundance / solar_electron_abundance) *
- // cooling->solar_abundances[i];
- //}
-
- /* Photodissociation */
- //solar_electron_abundance =
- // interpol_2d(cooling->table.photodissociation_cooling.electron_abundance, temp_i, n_h_i, d_temp, d_n_h, cooling->N_Temp, cooling->N_nH); /* ne/n_h */
-
- //for (i = 0; i < cooling->N_Elements; i++){
- // element_lambda[i+2] = interpol_3d(cooling->table.photodissociation_cooling.metal_heating, i,
- // temp_i, n_h_i, 0.0, d_temp, d_n_h,cooling->N_Elements,cooling->N_Temp,cooling->N_nH) *
- // (h_plus_he_electron_abundance / solar_electron_abundance) *
- // cooling->solar_abundances[i];
- //}
-
- /* redshift tables */
- solar_electron_abundance = interpol_3d(cooling->table.element_cooling.electron_abundance, z_index,n_h_i,dz,temp_i,d_n_h,d_temp,cooling->N_Redshifts,cooling->N_nH,cooling->N_Temp);
-
- for (i = 0; i < cooling->N_Elements; i++){
- // WARNING: before doing proper runs need to
- // multiply by ratio of particle abundances to solar
- temp_lambda = interpol_4d(cooling->table.element_cooling.metal_heating,z_index,i,n_h_i,temp_i,dz,0.0,d_n_h,d_temp,cooling->N_Redshifts,cooling->N_Elements,cooling->N_nH,cooling->N_Temp) *
- (h_plus_he_electron_abundance / solar_electron_abundance);
- cooling_rate += temp_lambda;
- if (element_lambda != NULL) element_lambda[i+2] = temp_lambda;
- }
+ // for (i = 0; i < cooling->N_Elements; i++){
+ // element_lambda[i+2] =
+ // interpol_2d(cooling->table.collisional_cooling.metal_heating, i,
+ // temp_i, 0.0, d_temp,cooling->N_Elements,cooling->N_Temp)
+ // *
+ // (h_plus_he_electron_abundance / solar_electron_abundance) *
+ // cooling->solar_abundances[i];
+ //}
- return cooling_rate;
-}
+ /* Photodissociation */
+ // solar_electron_abundance =
+ // interpol_2d(cooling->table.photodissociation_cooling.electron_abundance,
+ // temp_i, n_h_i, d_temp, d_n_h, cooling->N_Temp, cooling->N_nH); /* ne/n_h
+ // */
+
+ // for (i = 0; i < cooling->N_Elements; i++){
+ // element_lambda[i+2] =
+ // interpol_3d(cooling->table.photodissociation_cooling.metal_heating, i,
+ // temp_i, n_h_i, 0.0, d_temp,
+ // d_n_h,cooling->N_Elements,cooling->N_Temp,cooling->N_nH)
+ // *
+ // (h_plus_he_electron_abundance / solar_electron_abundance) *
+ // cooling->solar_abundances[i];
+ //}
+ /* redshift tables */
+ solar_electron_abundance =
+ interpol_3d(cooling->table.element_cooling.electron_abundance, z_index,
+ n_h_i, dz, temp_i, d_n_h, d_temp, cooling->N_Redshifts,
+ cooling->N_nH, cooling->N_Temp);
+
+ for (i = 0; i < cooling->N_Elements; i++) {
+ // WARNING: before doing proper runs need to
+ // multiply by ratio of particle abundances to solar
+ temp_lambda =
+ interpol_4d(cooling->table.element_cooling.metal_heating, z_index, i,
+ n_h_i, temp_i, dz, 0.0, d_n_h, d_temp, cooling->N_Redshifts,
+ cooling->N_Elements, cooling->N_nH, cooling->N_Temp) *
+ (h_plus_he_electron_abundance / solar_electron_abundance);
+ cooling_rate += temp_lambda;
+ if (element_lambda != NULL) element_lambda[i + 2] = temp_lambda;
+ }
+
+ return cooling_rate;
+}
/*
* @brief Calculates cooling rate from 1d tables
@@ -970,9 +1185,11 @@ __attribute__((always_inline)) INLINE static double eagle_metal_cooling_rate(dou
* @param u Internal energy
* @param dlambda_du Pointer to gradient of cooling rate
* @param H_plus_He_heat_table 1D table of heating rates due to H, He
- * @param H_plus_He_electron_abundance_table 1D table of electron abundances due to H,He
+ * @param H_plus_He_electron_abundance_table 1D table of electron abundances due
+ * to H,He
* @param element_cooling_table 1D table of heating rates due to metals
- * @param element_electron_abundance_table 1D table of electron abundances due to metals
+ * @param element_electron_abundance_table 1D table of electron abundances due
+ * to metals
* @param temp_table 1D table of temperatures
* @param z_index Redshift index of particle
* @param dz Redshift offset of particle
@@ -984,29 +1201,23 @@ __attribute__((always_inline)) INLINE static double eagle_metal_cooling_rate(dou
* @param cooling Cooling data structure
* @param cosmo Cosmology structure
* @param internal_const Physical constants structure
- * @param element_lambda Pointer for printing contribution to cooling rate
+ * @param element_lambda Pointer for printing contribution to cooling rate
* from each of the metals.
*/
-__attribute__((always_inline)) INLINE static double eagle_metal_cooling_rate_1d_table(double u,
- double *dlambda_du,
- double *H_plus_He_heat_table,
- double *H_plus_He_electron_abundance_table,
- double *element_cooling_table,
- double *element_electron_abundance_table,
- double *temp_table,
- int z_index,
- float dz,
- int n_h_i,
- float d_n_h,
- int He_i,
- float d_He,
- const struct part* restrict p,
- const struct cooling_function_data* restrict cooling,
- const struct cosmology* restrict cosmo,
- const struct phys_const *internal_const,
- double* element_lambda) {
+__attribute__((always_inline)) INLINE static double
+eagle_metal_cooling_rate_1d_table(
+ double u, double *dlambda_du, double *H_plus_He_heat_table,
+ double *H_plus_He_electron_abundance_table, double *element_cooling_table,
+ double *element_electron_abundance_table, double *temp_table, int z_index,
+ float dz, int n_h_i, float d_n_h, int He_i, float d_He,
+ const struct part *restrict p,
+ const struct cooling_function_data *restrict cooling,
+ const struct cosmology *restrict cosmo,
+ const struct phys_const *internal_const, double *element_lambda) {
double T_gam, solar_electron_abundance;
- double n_h = chemistry_get_number_density(p,cosmo,chemistry_element_H,internal_const)*cooling->number_density_scale; // chemistry_data
+ double n_h = chemistry_get_number_density(p, cosmo, chemistry_element_H,
+ internal_const) *
+ cooling->number_density_scale; // chemistry_data
double z = cosmo->z;
double cooling_rate = 0.0, temp_lambda;
float du;
@@ -1018,11 +1229,12 @@ __attribute__((always_inline)) INLINE static double eagle_metal_cooling_rate_1d_
float d_temp;
*dlambda_du = 0.0;
-
+
#if SWIFT_DEBUG_CHECKS
if (isnan(u)) printf("u is nan id %llu\n", p->id);
#endif
- temp = eagle_convert_u_to_temp_1d_table(u,&du,temp_table,p,cooling,cosmo,internal_const);
+ temp = eagle_convert_u_to_temp_1d_table(u, &du, temp_table, p, cooling, cosmo,
+ internal_const);
get_index_1d(cooling->Temp, cooling->N_Temp, log10(temp), &temp_i, &d_temp);
@@ -1030,28 +1242,36 @@ __attribute__((always_inline)) INLINE static double eagle_metal_cooling_rate_1d_
/* Metal-free cooling */
/* ------------------ */
- /* Collisional cooling */
- //element_lambda[0] =
- // interpol_2d(cooling->table.collisional_cooling.H_plus_He_heating, He_i,
- // temp_i, d_He, d_temp,cooling->N_He,cooling->N_Temp);
- //h_plus_he_electron_abundance =
- // interpol_2d(cooling->table.collisional_cooling.H_plus_He_electron_abundance, He_i,
- // temp_i, d_He, d_temp,cooling->N_He,cooling->N_Temp);
-
- /* Photodissociation */
- //element_lambda[0] =
- // interpol_3d(cooling->table.photodissociation_cooling.H_plus_He_heating, He_i,
- // temp_i, n_h_i, d_He, d_temp, d_n_h,cooling->N_He,cooling->N_Temp,cooling->N_nH);
- //h_plus_he_electron_abundance =
- // interpol_3d(cooling->table.photodissociation_cooling.H_plus_He_electron_abundance, He_i,
- // temp_i, n_h_i, d_He, d_temp, d_n_h,cooling->N_He,cooling->N_Temp,cooling->N_nH);
-
- /* redshift tables */
- temp_lambda = interpol_1d_dbl(H_plus_He_heat_table, temp_i, d_temp);
- h_plus_he_electron_abundance = interpol_1d_dbl(H_plus_He_electron_abundance_table, temp_i, d_temp);
- cooling_rate += temp_lambda;
- *dlambda_du += (((double) H_plus_He_heat_table[temp_i+1]) - ((double) H_plus_He_heat_table[temp_i]))/((double) du);
- if (element_lambda != NULL) element_lambda[0] = temp_lambda;
+ /* Collisional cooling */
+ // element_lambda[0] =
+ // interpol_2d(cooling->table.collisional_cooling.H_plus_He_heating, He_i,
+ // temp_i, d_He, d_temp,cooling->N_He,cooling->N_Temp);
+ // h_plus_he_electron_abundance =
+ // interpol_2d(cooling->table.collisional_cooling.H_plus_He_electron_abundance,
+ // He_i,
+ // temp_i, d_He, d_temp,cooling->N_He,cooling->N_Temp);
+
+ /* Photodissociation */
+ // element_lambda[0] =
+ // interpol_3d(cooling->table.photodissociation_cooling.H_plus_He_heating,
+ // He_i,
+ // temp_i, n_h_i, d_He, d_temp,
+ // d_n_h,cooling->N_He,cooling->N_Temp,cooling->N_nH);
+ // h_plus_he_electron_abundance =
+ // interpol_3d(cooling->table.photodissociation_cooling.H_plus_He_electron_abundance,
+ // He_i,
+ // temp_i, n_h_i, d_He, d_temp,
+ // d_n_h,cooling->N_He,cooling->N_Temp,cooling->N_nH);
+
+ /* redshift tables */
+ temp_lambda = interpol_1d_dbl(H_plus_He_heat_table, temp_i, d_temp);
+ h_plus_he_electron_abundance =
+ interpol_1d_dbl(H_plus_He_electron_abundance_table, temp_i, d_temp);
+ cooling_rate += temp_lambda;
+ *dlambda_du += (((double)H_plus_He_heat_table[temp_i + 1]) -
+ ((double)H_plus_He_heat_table[temp_i])) /
+ ((double)du);
+ if (element_lambda != NULL) element_lambda[0] = temp_lambda;
/* ------------------ */
/* Compton cooling */
@@ -1064,8 +1284,9 @@ __attribute__((always_inline)) INLINE static double eagle_metal_cooling_rate_1d_
T_gam = eagle_cmb_temperature * (1 + z);
- temp_lambda = -eagle_compton_rate * (temp - eagle_cmb_temperature * (1 + z)) * pow((1 + z), 4) *
- h_plus_he_electron_abundance / n_h;
+ temp_lambda = -eagle_compton_rate *
+ (temp - eagle_cmb_temperature * (1 + z)) * pow((1 + z), 4) *
+ h_plus_he_electron_abundance / n_h;
cooling_rate += temp_lambda;
if (element_lambda != NULL) element_lambda[1] = temp_lambda;
}
@@ -1074,61 +1295,80 @@ __attribute__((always_inline)) INLINE static double eagle_metal_cooling_rate_1d_
/* Metal cooling */
/* ------------- */
- /* for each element, find the abundance and multiply it
- by the interpolated heating-cooling */
+ /* for each element, find the abundance and multiply it
+ by the interpolated heating-cooling */
- //set_cooling_SolarAbundances(p->chemistry_data.metal_mass_fraction, cooling->solar_abundances, cooling, p);
+ // set_cooling_SolarAbundances(p->chemistry_data.metal_mass_fraction,
+ // cooling->solar_abundances, cooling, p);
- /* Collisional cooling */
- //solar_electron_abundance =
- // interpol_1d(cooling->table.collisional_cooling.electron_abundance, temp_i, d_temp); /* ne/n_h */
+ /* Collisional cooling */
+ // solar_electron_abundance =
+ // interpol_1d(cooling->table.collisional_cooling.electron_abundance,
+ // temp_i, d_temp); /* ne/n_h */
- //for (i = 0; i < cooling->N_Elements; i++){
- // element_lambda[i+2] = interpol_2d(cooling->table.collisional_cooling.metal_heating, i,
- // temp_i, 0.0, d_temp,cooling->N_Elements,cooling->N_Temp) *
- // (h_plus_he_electron_abundance / solar_electron_abundance) *
- // cooling->solar_abundances[i];
- //}
-
- /* Photodissociation */
- //solar_electron_abundance =
- // interpol_2d(cooling->table.photodissociation_cooling.electron_abundance, temp_i, n_h_i, d_temp, d_n_h, cooling->N_Temp, cooling->N_nH); /* ne/n_h */
-
- //for (i = 0; i < cooling->N_Elements; i++){
- // element_lambda[i+2] = interpol_3d(cooling->table.photodissociation_cooling.metal_heating, i,
- // temp_i, n_h_i, 0.0, d_temp, d_n_h,cooling->N_Elements,cooling->N_Temp,cooling->N_nH) *
- // (h_plus_he_electron_abundance / solar_electron_abundance) *
- // cooling->solar_abundances[i];
- //}
-
- /* redshift tables */
- solar_electron_abundance = interpol_1d_dbl(element_electron_abundance_table, temp_i, d_temp);
-
- for (i = 0; i < cooling->N_Elements; i++){
- // WARNING: before doing proper runs need to
- // multiply by ratio of particle abundances to solar
- temp_lambda = interpol_1d_dbl(element_cooling_table, temp_i, d_temp) *
- (h_plus_he_electron_abundance / solar_electron_abundance);
- cooling_rate += temp_lambda;
- *dlambda_du += (((double) element_cooling_table[temp_i+1])*((double) H_plus_He_electron_abundance_table[temp_i+1])/((double) element_electron_abundance_table[temp_i+1]) -
- ((double) element_cooling_table[temp_i])*((double) H_plus_He_electron_abundance_table[temp_i])/((double) element_electron_abundance_table[temp_i]))/((double) du);
- if (element_lambda != NULL) element_lambda[i+2] = temp_lambda;
- }
+ // for (i = 0; i < cooling->N_Elements; i++){
+ // element_lambda[i+2] =
+ // interpol_2d(cooling->table.collisional_cooling.metal_heating, i,
+ // temp_i, 0.0, d_temp,cooling->N_Elements,cooling->N_Temp)
+ // *
+ // (h_plus_he_electron_abundance / solar_electron_abundance) *
+ // cooling->solar_abundances[i];
+ //}
+
+ /* Photodissociation */
+ // solar_electron_abundance =
+ // interpol_2d(cooling->table.photodissociation_cooling.electron_abundance,
+ // temp_i, n_h_i, d_temp, d_n_h, cooling->N_Temp, cooling->N_nH); /* ne/n_h
+ // */
+
+ // for (i = 0; i < cooling->N_Elements; i++){
+ // element_lambda[i+2] =
+ // interpol_3d(cooling->table.photodissociation_cooling.metal_heating, i,
+ // temp_i, n_h_i, 0.0, d_temp,
+ // d_n_h,cooling->N_Elements,cooling->N_Temp,cooling->N_nH)
+ // *
+ // (h_plus_he_electron_abundance / solar_electron_abundance) *
+ // cooling->solar_abundances[i];
+ //}
- return cooling_rate;
+ /* redshift tables */
+ solar_electron_abundance =
+ interpol_1d_dbl(element_electron_abundance_table, temp_i, d_temp);
+
+ for (i = 0; i < cooling->N_Elements; i++) {
+ // WARNING: before doing proper runs need to
+ // multiply by ratio of particle abundances to solar
+ temp_lambda = interpol_1d_dbl(element_cooling_table, temp_i, d_temp) *
+ (h_plus_he_electron_abundance / solar_electron_abundance);
+ cooling_rate += temp_lambda;
+ *dlambda_du +=
+ (((double)element_cooling_table[temp_i + 1]) *
+ ((double)H_plus_He_electron_abundance_table[temp_i + 1]) /
+ ((double)element_electron_abundance_table[temp_i + 1]) -
+ ((double)element_cooling_table[temp_i]) *
+ ((double)H_plus_He_electron_abundance_table[temp_i]) /
+ ((double)element_electron_abundance_table[temp_i])) /
+ ((double)du);
+ if (element_lambda != NULL) element_lambda[i + 2] = temp_lambda;
+ }
+
+ return cooling_rate;
}
/*
- * @brief Wrapper function previously used to calculate cooling rate and dLambda_du
+ * @brief Wrapper function previously used to calculate cooling rate and
+ * dLambda_du
* Can be used to check whether calculation of dLambda_du is done correctly.
* Should be removed when finished
*
* @param u Internal energy
* @param dlambda_du Pointer to gradient of cooling rate
* @param H_plus_He_heat_table 1D table of heating rates due to H, He
- * @param H_plus_He_electron_abundance_table 1D table of electron abundances due to H,He
+ * @param H_plus_He_electron_abundance_table 1D table of electron abundances due
+ * to H,He
* @param element_cooling_table 1D table of heating rates due to metals
- * @param element_electron_abundance_table 1D table of electron abundances due to metals
+ * @param element_electron_abundance_table 1D table of electron abundances due
+ * to metals
* @param temp_table 1D table of temperatures
* @param z_index Redshift index of particle
* @param dz Redshift offset of particle
@@ -1141,44 +1381,50 @@ __attribute__((always_inline)) INLINE static double eagle_metal_cooling_rate_1d_
* @param cosmo Cosmology structure
* @param internal_const Physical constants structure
*/
-__attribute__((always_inline)) INLINE static double eagle_cooling_rate_1d_table(double u,
- double *dLambdaNet_du,
- double *H_plus_He_heat_table,
- double *H_plus_He_electron_abundance_table,
- double *element_cooling_table,
- double *element_electron_abundance_table,
- double *temp_table,
- int z_index,
- float dz,
- int n_h_i,
- float d_n_h,
- int He_i,
- float d_He,
- const struct part* restrict p,
- const struct cooling_function_data* restrict cooling,
- const struct cosmology* restrict cosmo,
- const struct phys_const *internal_const) {
- double *element_lambda = NULL, lambda_net1 = 0.0;//, lambda_net2 = 0.0, delta, dLambdaNet_du_calc;
- //float d_u;
- //int u_i;
- //get_index_1d(cooling->Therm, cooling->N_Temp, log10(u), &u_i, &d_u);
- //if (u_i >= cooling->N_Temp-2){
+__attribute__((always_inline)) INLINE static double eagle_cooling_rate_1d_table(
+ double u, double *dLambdaNet_du, double *H_plus_He_heat_table,
+ double *H_plus_He_electron_abundance_table, double *element_cooling_table,
+ double *element_electron_abundance_table, double *temp_table, int z_index,
+ float dz, int n_h_i, float d_n_h, int He_i, float d_He,
+ const struct part *restrict p,
+ const struct cooling_function_data *restrict cooling,
+ const struct cosmology *restrict cosmo,
+ const struct phys_const *internal_const) {
+ double *element_lambda = NULL,
+ lambda_net1 = 0.0; //, lambda_net2 = 0.0, delta, dLambdaNet_du_calc;
+ // float d_u;
+ // int u_i;
+ // get_index_1d(cooling->Therm, cooling->N_Temp, log10(u), &u_i, &d_u);
+ // if (u_i >= cooling->N_Temp-2){
// delta = pow(10.0,cooling->Therm[cooling->N_Temp-2]) - u;
//} else {
// delta = pow(10.0,cooling->Therm[u_i+1]) - pow(10.0,cooling->Therm[u_i]);
//}
- //lambda_net1 = eagle_metal_cooling_rate_1d_table(pow(10.0,cooling->Therm[u_i]), dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table, element_cooling_table, element_electron_abundance_table, temp_table, p, cooling, cosmo, internal_const, element_lambda);
- //lambda_net2 = eagle_metal_cooling_rate_1d_table(pow(10.0,cooling->Therm[u_i+1]), dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table, element_cooling_table, element_electron_abundance_table, temp_table, p, cooling, cosmo, internal_const, element_lambda);
+ // lambda_net1 =
+ // eagle_metal_cooling_rate_1d_table(pow(10.0,cooling->Therm[u_i]),
+ // dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table,
+ // element_cooling_table, element_electron_abundance_table, temp_table, p,
+ // cooling, cosmo, internal_const, element_lambda);
+ // lambda_net2 =
+ // eagle_metal_cooling_rate_1d_table(pow(10.0,cooling->Therm[u_i+1]),
+ // dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table,
+ // element_cooling_table, element_electron_abundance_table, temp_table, p,
+ // cooling, cosmo, internal_const, element_lambda);
//
- //dLambdaNet_du_calc = (lambda_net2 - lambda_net1)/delta;
+ // dLambdaNet_du_calc = (lambda_net2 - lambda_net1)/delta;
- //printf( "Eagle cooling.h 2 dlambda_du, dlambda, du, index %0.5e %0.5e %.5e %d\n", dLambdaNet_du_calc, lambda_net2 - lambda_net1, delta, u_i);
-
- lambda_net1 = eagle_metal_cooling_rate_1d_table(u, dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table, element_cooling_table, element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h, He_i, d_He, p, cooling, cosmo, internal_const, element_lambda);
+ // printf( "Eagle cooling.h 2 dlambda_du, dlambda, du, index %0.5e %0.5e %.5e
+ // %d\n", dLambdaNet_du_calc, lambda_net2 - lambda_net1, delta, u_i);
- //printf("Eagle cooling.h double check dlambda_du endpoints, lambda_nets %.5e %.5e\n", *dLambdaNet_du, dLambdaNet_du_calc);
+ lambda_net1 = eagle_metal_cooling_rate_1d_table(
+ u, dLambdaNet_du, H_plus_He_heat_table,
+ H_plus_He_electron_abundance_table, element_cooling_table,
+ element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h,
+ He_i, d_He, p, cooling, cosmo, internal_const, element_lambda);
+ // printf("Eagle cooling.h double check dlambda_du endpoints, lambda_nets %.5e
+ // %.5e\n", *dLambdaNet_du, dLambdaNet_du_calc);
return lambda_net1;
}
@@ -1188,9 +1434,11 @@ __attribute__((always_inline)) INLINE static double eagle_cooling_rate_1d_table(
* to files for checking against Wiersma et al 2009 idl scripts.
*
* @param H_plus_He_heat_table 1D table of heating rates due to H, He
- * @param H_plus_He_electron_abundance_table 1D table of electron abundances due to H,He
+ * @param H_plus_He_electron_abundance_table 1D table of electron abundances due
+ * to H,He
* @param element_cooling_table 1D table of heating rates due to metals
- * @param element_electron_abundance_table 1D table of electron abundances due to metals
+ * @param element_electron_abundance_table 1D table of electron abundances due
+ * to metals
* @param temp_table 1D table of temperatures
* @param z_index Redshift index of particle
* @param dz Redshift offset of particle
@@ -1203,58 +1451,59 @@ __attribute__((always_inline)) INLINE static double eagle_cooling_rate_1d_table(
* @param cosmo Cosmology structure
* @param internal_const Physical constants structure
*/
-__attribute__((always_inline)) INLINE static double eagle_print_metal_cooling_rate_1d_table(double *H_plus_He_heat_table,
- double *H_plus_He_electron_abundance_table,
- double *element_cooling_table,
- double *element_electron_abundance_table,
- double *temp_table,
- int z_index,
- float dz,
- int n_h_i,
- float d_n_h,
- int He_i,
- float d_He,
- const struct part* restrict p,
- const struct cooling_function_data* restrict cooling,
- const struct cosmology* restrict cosmo,
- const struct phys_const *internal_const) {
+__attribute__((always_inline)) INLINE static double
+eagle_print_metal_cooling_rate_1d_table(
+ double *H_plus_He_heat_table, double *H_plus_He_electron_abundance_table,
+ double *element_cooling_table, double *element_electron_abundance_table,
+ double *temp_table, int z_index, float dz, int n_h_i, float d_n_h, int He_i,
+ float d_He, const struct part *restrict p,
+ const struct cooling_function_data *restrict cooling,
+ const struct cosmology *restrict cosmo,
+ const struct phys_const *internal_const) {
double *element_lambda, lambda_net = 0.0;
- element_lambda = malloc((cooling->N_Elements+2)*sizeof(double));
- double u = hydro_get_physical_internal_energy(p,cosmo)*cooling->internal_energy_scale;
+ element_lambda = malloc((cooling->N_Elements + 2) * sizeof(double));
+ double u = hydro_get_physical_internal_energy(p, cosmo) *
+ cooling->internal_energy_scale;
double dLambdaNet_du;
-
+
char output_filename[21];
- FILE** output_file = malloc((cooling->N_Elements+2)*sizeof(FILE*));
- for (int element = 0; element < cooling->N_Elements+2; element++){
- sprintf(output_filename, "%s%d%s", "cooling_output_",element,".dat");
+ FILE **output_file = malloc((cooling->N_Elements + 2) * sizeof(FILE *));
+ for (int element = 0; element < cooling->N_Elements + 2; element++) {
+ sprintf(output_filename, "%s%d%s", "cooling_output_", element, ".dat");
output_file[element] = fopen(output_filename, "a");
- if (output_file == NULL)
- {
- printf("Error opening file!\n");
- exit(1);
+ if (output_file == NULL) {
+ printf("Error opening file!\n");
+ exit(1);
}
}
- for (int j = 0; j < cooling->N_Elements+2; j++) element_lambda[j] = 0.0;
- lambda_net = eagle_metal_cooling_rate_1d_table(u,&dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table, element_cooling_table, element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h, He_i, d_He, p, cooling, cosmo, internal_const, element_lambda);
- for (int j = 0; j < cooling->N_Elements+2; j++) {
- fprintf(output_file[j],"%.5e\n",element_lambda[j]);
+ for (int j = 0; j < cooling->N_Elements + 2; j++) element_lambda[j] = 0.0;
+ lambda_net = eagle_metal_cooling_rate_1d_table(
+ u, &dLambdaNet_du, H_plus_He_heat_table,
+ H_plus_He_electron_abundance_table, element_cooling_table,
+ element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h,
+ He_i, d_He, p, cooling, cosmo, internal_const, element_lambda);
+ for (int j = 0; j < cooling->N_Elements + 2; j++) {
+ fprintf(output_file[j], "%.5e\n", element_lambda[j]);
}
-
- for (int i = 0; i < cooling->N_Elements+2; i++) fclose(output_file[i]);
+
+ for (int i = 0; i < cooling->N_Elements + 2; i++) fclose(output_file[i]);
return lambda_net;
}
/*
- * @brief Bisection integration scheme for when Newton Raphson method fails to converge
+ * @brief Bisection integration scheme for when Newton Raphson method fails to
+ * converge
*
* @param x_init Initial guess for log(internal energy)
* @param u_ini Internal energy at beginning of hydro step
* @param H_plus_He_heat_table 1D table of heating rates due to H, He
- * @param H_plus_He_electron_abundance_table 1D table of electron abundances due to H,He
+ * @param H_plus_He_electron_abundance_table 1D table of electron abundances due
+ * to H,He
* @param element_cooling_table 1D table of heating rates due to metals
- * @param element_electron_abundance_table 1D table of electron abundances due to metals
+ * @param element_electron_abundance_table 1D table of electron abundances due
+ * to metals
* @param temp_table 1D table of temperatures
* @param z_index Redshift index of particle
* @param dz Redshift offset of particle
@@ -1268,54 +1517,63 @@ __attribute__((always_inline)) INLINE static double eagle_print_metal_cooling_ra
* @param phys_const Physical constants structure
* @param dt Hydro timestep
*/
-__attribute__((always_inline)) INLINE static float bisection_iter(float x_init,
- double u_ini,
- double *H_plus_He_heat_table,
- double *H_plus_He_electron_abundance_table,
- double *element_cooling_table,
- double *element_electron_abundance_table,
- double *temp_table,
- int z_index,
- float dz,
- int n_h_i,
- float d_n_h,
- int He_i,
- float d_He,
- struct part* restrict p,
- const struct cosmology* restrict cosmo,
- const struct cooling_function_data* restrict cooling,
- const struct phys_const* restrict phys_const,
- float dt){
+__attribute__((always_inline)) INLINE static float bisection_iter(
+ float x_init, double u_ini, double *H_plus_He_heat_table,
+ double *H_plus_He_electron_abundance_table, double *element_cooling_table,
+ double *element_electron_abundance_table, double *temp_table, int z_index,
+ float dz, int n_h_i, float d_n_h, int He_i, float d_He,
+ struct part *restrict p, const struct cosmology *restrict cosmo,
+ const struct cooling_function_data *restrict cooling,
+ const struct phys_const *restrict phys_const, float dt) {
double x_a, x_b, x_next, f_a, f_b, f_next, dLambdaNet_du;
int i = 0;
float shift = 0.3;
x_a = x_init;
- x_b = cooling->Therm[0]/log10(exp(1.0));
- f_a = eagle_cooling_rate_1d_table(exp(x_a), &dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table, element_cooling_table, element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h, He_i, d_He, p, cooling, cosmo, phys_const);
- f_b = eagle_cooling_rate_1d_table(exp(x_b), &dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table, element_cooling_table, element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h, He_i, d_He, p, cooling, cosmo, phys_const);
-
- while (f_a*f_b >= 0 & i < 10*eagle_max_iterations) {
- if ((x_a+shift)/log(10) < cooling->Therm[cooling->N_Temp-1]) {
+ x_b = cooling->Therm[0] / log10(exp(1.0));
+ f_a = eagle_cooling_rate_1d_table(
+ exp(x_a), &dLambdaNet_du, H_plus_He_heat_table,
+ H_plus_He_electron_abundance_table, element_cooling_table,
+ element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h,
+ He_i, d_He, p, cooling, cosmo, phys_const);
+ f_b = eagle_cooling_rate_1d_table(
+ exp(x_b), &dLambdaNet_du, H_plus_He_heat_table,
+ H_plus_He_electron_abundance_table, element_cooling_table,
+ element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h,
+ He_i, d_He, p, cooling, cosmo, phys_const);
+
+ while (f_a * f_b >= 0 & i < 10 * eagle_max_iterations) {
+ if ((x_a + shift) / log(10) < cooling->Therm[cooling->N_Temp - 1]) {
x_a += shift;
- f_a = eagle_cooling_rate_1d_table(exp(x_a), &dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table, element_cooling_table, element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h, He_i, d_He, p, cooling, cosmo, phys_const);
+ f_a = eagle_cooling_rate_1d_table(
+ exp(x_a), &dLambdaNet_du, H_plus_He_heat_table,
+ H_plus_He_electron_abundance_table, element_cooling_table,
+ element_electron_abundance_table, temp_table, z_index, dz, n_h_i,
+ d_n_h, He_i, d_He, p, cooling, cosmo, phys_const);
}
- //if ((x_b-shift)/log(10) > cooling->Therm[0]) {
+ // if ((x_b-shift)/log(10) > cooling->Therm[0]) {
// x_b -= shift;
// if (isnan(exp(x_b))) printf("4 u is nan id %llu\n",p->id);
- // f_b = eagle_cooling_rate_1d_table(exp(x_b), &dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table, element_cooling_table, element_electron_abundance_table, temp_table, p, cooling, cosmo, phys_const);
+ // f_b = eagle_cooling_rate_1d_table(exp(x_b), &dLambdaNet_du,
+ // H_plus_He_heat_table, H_plus_He_electron_abundance_table,
+ // element_cooling_table, element_electron_abundance_table, temp_table, p,
+ // cooling, cosmo, phys_const);
//}
i++;
}
#if SWIFT_DEBUG_CHECKS
- assert(f_a*f_b < 0);
+ assert(f_a * f_b < 0);
#endif
i = 0;
- do{
- x_next = 0.5*(x_a + x_b);
- f_next = eagle_cooling_rate_1d_table(exp(x_next), &dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table, element_cooling_table, element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h, He_i, d_He, p, cooling, cosmo, phys_const);
- if (f_a*f_next < 0) {
+ do {
+ x_next = 0.5 * (x_a + x_b);
+ f_next = eagle_cooling_rate_1d_table(
+ exp(x_next), &dLambdaNet_du, H_plus_He_heat_table,
+ H_plus_He_electron_abundance_table, element_cooling_table,
+ element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h,
+ He_i, d_He, p, cooling, cosmo, phys_const);
+ if (f_a * f_next < 0) {
x_b = x_next;
f_b = f_next;
} else {
@@ -1323,9 +1581,10 @@ __attribute__((always_inline)) INLINE static float bisection_iter(float x_init,
f_a = f_next;
}
i++;
- //printf("Eagle cooling.h id x_a x_b f_a f_b i %llu %.5e %.5e %.5e %.5e %d\n", p->id, x_a, x_b, f_a, f_b, i );
+ // printf("Eagle cooling.h id x_a x_b f_a f_b i %llu %.5e %.5e %.5e %.5e
+ // %d\n", p->id, x_a, x_b, f_a, f_b, i );
} while (fabs(x_a - x_b) > 1.0e-2 && i < eagle_max_iterations);
-
+
return x_b;
}
@@ -1337,9 +1596,11 @@ __attribute__((always_inline)) INLINE static float bisection_iter(float x_init,
* @param x_init Initial guess for log(internal energy)
* @param u_ini Internal energy at beginning of hydro step
* @param H_plus_He_heat_table 1D table of heating rates due to H, He
- * @param H_plus_He_electron_abundance_table 1D table of electron abundances due to H,He
+ * @param H_plus_He_electron_abundance_table 1D table of electron abundances due
+ * to H,He
* @param element_cooling_table 1D table of heating rates due to metals
- * @param element_electron_abundance_table 1D table of electron abundances due to metals
+ * @param element_electron_abundance_table 1D table of electron abundances due
+ * to metals
* @param temp_table 1D table of temperatures
* @param z_index Redshift index of particle
* @param dz Redshift offset of particle
@@ -1353,131 +1614,169 @@ __attribute__((always_inline)) INLINE static float bisection_iter(float x_init,
* @param phys_const Physical constants structure
* @param dt Hydro timestep
*/
-__attribute__((always_inline)) INLINE static float brent_iter(float x_init,
- double u_ini,
- double *H_plus_He_heat_table,
- double *H_plus_He_electron_abundance_table,
- double *element_cooling_table,
- double *element_electron_abundance_table,
- double *temp_table,
- int z_index,
- float dz,
- int n_h_i,
- float d_n_h,
- int He_i,
- float d_He,
- struct part* restrict p,
- const struct cosmology* restrict cosmo,
- const struct cooling_function_data* restrict cooling,
- const struct phys_const* restrict phys_const,
- float dt){
- /* this routine does the iteration scheme, call one and it iterates to convergence */
-
- double x_a,x_b,x_c,x_d,x_s, x_temp, delta = 1.0e-2;
+__attribute__((always_inline)) INLINE static float brent_iter(
+ float x_init, double u_ini, double *H_plus_He_heat_table,
+ double *H_plus_He_electron_abundance_table, double *element_cooling_table,
+ double *element_electron_abundance_table, double *temp_table, int z_index,
+ float dz, int n_h_i, float d_n_h, int He_i, float d_He,
+ struct part *restrict p, const struct cosmology *restrict cosmo,
+ const struct cooling_function_data *restrict cooling,
+ const struct phys_const *restrict phys_const, float dt) {
+ /* this routine does the iteration scheme, call one and it iterates to
+ * convergence */
+
+ double x_a, x_b, x_c, x_d, x_s, x_temp, delta = 1.0e-2;
double dLambdaNet_du, Lambda_a, Lambda_b, Lambda_c, Lambda_s;
- double f_a,f_b,f_c,f_s, f_temp;
+ double f_a, f_b, f_c, f_s, f_temp;
float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
/* convert Hydrogen mass fraction in Hydrogen number density */
- double inn_h = chemistry_get_number_density(p,cosmo,chemistry_element_H,phys_const)*cooling->number_density_scale;
+ double inn_h =
+ chemistry_get_number_density(p, cosmo, chemistry_element_H, phys_const) *
+ cooling->number_density_scale;
/* ratefact = inn_h * inn_h / rho; Might lead to round-off error: replaced by
* equivalent expression below */
double ratefact = inn_h * (XH / eagle_proton_mass_cgs);
/* set helium and hydrogen reheating term */
- //LambdaTune = eagle_helium_reionization_extraheat(z_index, dz); // INCLUDE HELIUM REIONIZATION????
- //if (zold > z_index) {
+ // LambdaTune = eagle_helium_reionization_extraheat(z_index, dz); // INCLUDE
+ // HELIUM REIONIZATION????
+ // if (zold > z_index) {
// LambdaCumul += LambdaTune;
- // printf(" EagleDoCooling %d %g %g %g\n", z_index, dz, LambdaTune, LambdaCumul);
+ // printf(" EagleDoCooling %d %g %g %g\n", z_index, dz, LambdaTune,
+ // LambdaCumul);
// zold = z_index;
//}
int i = 0;
x_a = x_init;
- //x_b = cooling->Therm[0]*log(10);
- x_b = 0.5*x_init;
+ // x_b = cooling->Therm[0]*log(10);
+ x_b = 0.5 * x_init;
x_c = x_a;
x_d = 0.0;
- Lambda_a = eagle_cooling_rate_1d_table(exp(x_a), &dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table, element_cooling_table, element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h, He_i, d_He, p, cooling, cosmo, phys_const);
- Lambda_b = eagle_cooling_rate_1d_table(exp(x_b), &dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table, element_cooling_table, element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h, He_i, d_He, p, cooling, cosmo, phys_const);
- f_a = exp(x_a) - exp(x_init) - Lambda_a*ratefact*dt;
- f_b = exp(x_b) - exp(x_init) - Lambda_b*ratefact*dt;
- assert(f_a*f_b < 0);
+ Lambda_a = eagle_cooling_rate_1d_table(
+ exp(x_a), &dLambdaNet_du, H_plus_He_heat_table,
+ H_plus_He_electron_abundance_table, element_cooling_table,
+ element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h,
+ He_i, d_He, p, cooling, cosmo, phys_const);
+ Lambda_b = eagle_cooling_rate_1d_table(
+ exp(x_b), &dLambdaNet_du, H_plus_He_heat_table,
+ H_plus_He_electron_abundance_table, element_cooling_table,
+ element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h,
+ He_i, d_He, p, cooling, cosmo, phys_const);
+ f_a = exp(x_a) - exp(x_init) - Lambda_a * ratefact * dt;
+ f_b = exp(x_b) - exp(x_init) - Lambda_b * ratefact * dt;
+ assert(f_a * f_b < 0);
if (f_a < f_b) {
- x_temp = x_a; x_a = x_b; x_b = x_temp;
- f_temp = f_a; f_a = f_b; f_b = f_temp;
+ x_temp = x_a;
+ x_a = x_b;
+ x_b = x_temp;
+ f_temp = f_a;
+ f_a = f_b;
+ f_b = f_temp;
}
int mflag = 1;
do /* iterate to convergence */
- {
- Lambda_c = eagle_cooling_rate_1d_table(exp(x_c), &dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table, element_cooling_table, element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h, He_i, d_He, p, cooling, cosmo, phys_const);
- f_c = exp(x_c) - exp(x_init) - Lambda_c*ratefact*dt;
-
- if ((f_a != f_c) && (f_b != f_c)) {
- x_s = x_a*f_b*f_c/((f_a - f_b)*(f_a - f_c)) +
- x_b*f_a*f_c/((f_b - f_a)*(f_b - f_c)) +
- x_c*f_b*f_a/((f_c - f_b)*(f_c - f_a));
- printf("Eagle cooling.h 1 id x_s, u_s, error, %llu %.5e %.5e %.5e %d\n", p->id, x_s, exp(x_s), x_a - x_b, i);
- printf("Eagle cooling.h 1 id u_a, u_b, u_c, f_a, f_b, f_c, %llu %.5e %.5e %.5e %.5e %.5e %.5e %.5e %d\n", p->id, exp(x_a), exp(x_b), exp(x_c), f_a, f_b, f_c, x_a - x_b, i);
- } else {
- x_s = x_b - f_b*(x_b - x_a)/(f_b - f_a);
- printf("Eagle cooling.h 2 id x_s, u_s, error, %llu %.5e %.5e %.5e %d\n", p->id, x_s, exp(x_s), x_a - x_b, i);
- }
+ {
+ Lambda_c = eagle_cooling_rate_1d_table(
+ exp(x_c), &dLambdaNet_du, H_plus_He_heat_table,
+ H_plus_He_electron_abundance_table, element_cooling_table,
+ element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h,
+ He_i, d_He, p, cooling, cosmo, phys_const);
+ f_c = exp(x_c) - exp(x_init) - Lambda_c * ratefact * dt;
+
+ if ((f_a != f_c) && (f_b != f_c)) {
+ x_s = x_a * f_b * f_c / ((f_a - f_b) * (f_a - f_c)) +
+ x_b * f_a * f_c / ((f_b - f_a) * (f_b - f_c)) +
+ x_c * f_b * f_a / ((f_c - f_b) * (f_c - f_a));
+ printf("Eagle cooling.h 1 id x_s, u_s, error, %llu %.5e %.5e %.5e %d\n",
+ p->id, x_s, exp(x_s), x_a - x_b, i);
+ printf(
+ "Eagle cooling.h 1 id u_a, u_b, u_c, f_a, f_b, f_c, %llu %.5e %.5e "
+ "%.5e %.5e %.5e %.5e %.5e %d\n",
+ p->id, exp(x_a), exp(x_b), exp(x_c), f_a, f_b, f_c, x_a - x_b, i);
+ } else {
+ x_s = x_b - f_b * (x_b - x_a) / (f_b - f_a);
+ printf("Eagle cooling.h 2 id x_s, u_s, error, %llu %.5e %.5e %.5e %d\n",
+ p->id, x_s, exp(x_s), x_a - x_b, i);
+ }
- if ( (x_s < (3.0*x_a + x_b)/4.0 || x_s > x_b) ||
- (mflag == 1 && fabs(x_s - x_b) >= 0.5*fabs(x_b - x_c)) ||
- (mflag != 1 && fabs(x_s - x_b) >= 0.5*fabs(x_c - x_d)) ||
- (mflag == 1 && fabs(x_b - x_c) < delta) ||
- (mflag != 1 && fabs(x_c - x_d) < delta) ) {
- x_s = 0.5*(x_a + x_b);
- printf("Eagle cooling.h 3 id x_s, u_s, error, %llu %.5e %.5e %.5e %d\n", p->id, x_s, exp(x_s), x_a - x_b, i);
- mflag = 1;
- } else {
- mflag = 0;
- }
- x_d = x_c;
- x_c = x_b;
- printf("Eagle cooling.h 4 id x_s, u_s, error, %llu %.5e %.5e %.5e %d\n", p->id, x_s, exp(x_s), x_a - x_b, i);
- Lambda_s = eagle_cooling_rate_1d_table(exp(x_s), &dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table, element_cooling_table, element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h, He_i, d_He, p, cooling, cosmo, phys_const);
- f_s = exp(x_s) - exp(x_init) - Lambda_s*ratefact*dt;
- if (f_s*f_a < 0) {
- x_b = x_s;
- } else {
- x_a = x_s;
- }
- Lambda_a = eagle_cooling_rate_1d_table(exp(x_a), &dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table, element_cooling_table, element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h, He_i, d_He, p, cooling, cosmo, phys_const);
- Lambda_b = eagle_cooling_rate_1d_table(exp(x_b), &dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table, element_cooling_table, element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h, He_i, d_He, p, cooling, cosmo, phys_const);
- f_a = exp(x_a) - exp(x_init) - Lambda_a*ratefact*dt;
- f_b = exp(x_b) - exp(x_init) - Lambda_b*ratefact*dt;
- if (f_a < f_b) {
- x_temp = x_a; x_a = x_b; x_b = x_temp;
- f_temp = f_a; f_a = f_b; f_b = f_temp;
- }
- i++;
+ if ((x_s < (3.0 * x_a + x_b) / 4.0 || x_s > x_b) ||
+ (mflag == 1 && fabs(x_s - x_b) >= 0.5 * fabs(x_b - x_c)) ||
+ (mflag != 1 && fabs(x_s - x_b) >= 0.5 * fabs(x_c - x_d)) ||
+ (mflag == 1 && fabs(x_b - x_c) < delta) ||
+ (mflag != 1 && fabs(x_c - x_d) < delta)) {
+ x_s = 0.5 * (x_a + x_b);
+ printf("Eagle cooling.h 3 id x_s, u_s, error, %llu %.5e %.5e %.5e %d\n",
+ p->id, x_s, exp(x_s), x_a - x_b, i);
+ mflag = 1;
+ } else {
+ mflag = 0;
+ }
+ x_d = x_c;
+ x_c = x_b;
+ printf("Eagle cooling.h 4 id x_s, u_s, error, %llu %.5e %.5e %.5e %d\n",
+ p->id, x_s, exp(x_s), x_a - x_b, i);
+ Lambda_s = eagle_cooling_rate_1d_table(
+ exp(x_s), &dLambdaNet_du, H_plus_He_heat_table,
+ H_plus_He_electron_abundance_table, element_cooling_table,
+ element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h,
+ He_i, d_He, p, cooling, cosmo, phys_const);
+ f_s = exp(x_s) - exp(x_init) - Lambda_s * ratefact * dt;
+ if (f_s * f_a < 0) {
+ x_b = x_s;
+ } else {
+ x_a = x_s;
+ }
+ Lambda_a = eagle_cooling_rate_1d_table(
+ exp(x_a), &dLambdaNet_du, H_plus_He_heat_table,
+ H_plus_He_electron_abundance_table, element_cooling_table,
+ element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h,
+ He_i, d_He, p, cooling, cosmo, phys_const);
+ Lambda_b = eagle_cooling_rate_1d_table(
+ exp(x_b), &dLambdaNet_du, H_plus_He_heat_table,
+ H_plus_He_electron_abundance_table, element_cooling_table,
+ element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h,
+ He_i, d_He, p, cooling, cosmo, phys_const);
+ f_a = exp(x_a) - exp(x_init) - Lambda_a * ratefact * dt;
+ f_b = exp(x_b) - exp(x_init) - Lambda_b * ratefact * dt;
+ if (f_a < f_b) {
+ x_temp = x_a;
+ x_a = x_b;
+ x_b = x_temp;
+ f_temp = f_a;
+ f_a = f_b;
+ f_b = f_temp;
+ }
+ i++;
- } while (fabs(x_b - x_a) > delta && i < eagle_max_iterations);
+ } while (fabs(x_b - x_a) > delta && i < eagle_max_iterations);
- if (i >= eagle_max_iterations) {
- n_eagle_cooling_rate_calls_3++;
- printf("Eagle cooling.h particle %llu with density %.5e, initial energy, %.5e and redshift %.5e not converged \n", p->id, inn_h, exp(x_init), cosmo->z);
- }
+ if (i >= eagle_max_iterations) {
+ n_eagle_cooling_rate_calls_3++;
+ printf(
+ "Eagle cooling.h particle %llu with density %.5e, initial energy, %.5e "
+ "and redshift %.5e not converged \n",
+ p->id, inn_h, exp(x_init), cosmo->z);
+ }
return x_s;
-
}
-
/*
- * @brief Newton Raphson integration scheme to calculate particle cooling over hydro timestep
+ * @brief Newton Raphson integration scheme to calculate particle cooling over
+ * hydro timestep
*
* @param x_init Initial guess for log(internal energy)
* @param u_ini Internal energy at beginning of hydro step
* @param H_plus_He_heat_table 1D table of heating rates due to H, He
- * @param H_plus_He_electron_abundance_table 1D table of electron abundances due to H,He
+ * @param H_plus_He_electron_abundance_table 1D table of electron abundances due
+ * to H,He
* @param element_cooling_table 1D table of heating rates due to metals
- * @param element_electron_abundance_table 1D table of electron abundances due to metals
+ * @param element_electron_abundance_table 1D table of electron abundances due
+ * to metals
* @param temp_table 1D table of temperatures
* @param z_index Redshift index of particle
* @param dz Redshift offset of particle
@@ -1493,94 +1792,102 @@ __attribute__((always_inline)) INLINE static float brent_iter(float x_init,
* @param printflag Flag to identify if scheme failed to converge
* @param relax_factor Factor to increase stability
*/
-__attribute__((always_inline)) INLINE static float newton_iter(float x_init,
- double u_ini,
- double *H_plus_He_heat_table,
- double *H_plus_He_electron_abundance_table,
- double *element_cooling_table,
- double *element_electron_abundance_table,
- double *temp_table,
- int z_index,
- float dz,
- int n_h_i,
- float d_n_h,
- int He_i,
- float d_He,
- struct part* restrict p,
- const struct cosmology* restrict cosmo,
- const struct cooling_function_data* restrict cooling,
- const struct phys_const* restrict phys_const,
- float dt,
- int *printflag,
- float relax_factor){
- /* this routine does the iteration scheme, call one and it iterates to convergence */
+__attribute__((always_inline)) INLINE static float newton_iter(
+ float x_init, double u_ini, double *H_plus_He_heat_table,
+ double *H_plus_He_electron_abundance_table, double *element_cooling_table,
+ double *element_electron_abundance_table, double *temp_table, int z_index,
+ float dz, int n_h_i, float d_n_h, int He_i, float d_He,
+ struct part *restrict p, const struct cosmology *restrict cosmo,
+ const struct cooling_function_data *restrict cooling,
+ const struct phys_const *restrict phys_const, float dt, int *printflag,
+ float relax_factor) {
+ /* this routine does the iteration scheme, call one and it iterates to
+ * convergence */
double x, x_old;
double dLambdaNet_du, LambdaNet;
float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
/* convert Hydrogen mass fraction in Hydrogen number density */
- double inn_h = chemistry_get_number_density(p,cosmo,chemistry_element_H,phys_const)*cooling->number_density_scale;
-
+ double inn_h =
+ chemistry_get_number_density(p, cosmo, chemistry_element_H, phys_const) *
+ cooling->number_density_scale;
+
/* ratefact = inn_h * inn_h / rho; Might lead to round-off error: replaced by
* equivalent expression below */
double ratefact = inn_h * (XH / eagle_proton_mass_cgs);
- float residual ; // added by RGB
-
+ float residual; // added by RGB
+
/* set helium and hydrogen reheating term */
- //LambdaTune = eagle_helium_reionization_extraheat(z_index, dz); // INCLUDE HELIUM REIONIZATION????
- //if (zold > z_index) {
+ // LambdaTune = eagle_helium_reionization_extraheat(z_index, dz); // INCLUDE
+ // HELIUM REIONIZATION????
+ // if (zold > z_index) {
// LambdaCumul += LambdaTune;
- // printf(" EagleDoCooling %d %g %g %g\n", z_index, dz, LambdaTune, LambdaCumul);
+ // printf(" EagleDoCooling %d %g %g %g\n", z_index, dz, LambdaTune,
+ // LambdaCumul);
// zold = z_index;
//}
-
- x_old = x_init ;
+
+ x_old = x_init;
x = x_old;
int i = 0;
do /* iterate to convergence */
- {
- x_old = x ;
- // this version is needed when reionization is included...
- // LambdaNet = (LambdaTune / (dt * ratefact)) + eagle_cooling_rate_1d_table(exp(x_old), &dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table, element_cooling_table, element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h, He_i, d_He, p, cooling, cosmo, phys_const);
- LambdaNet = eagle_cooling_rate_1d_table(exp(x_old), &dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table, element_cooling_table, element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h, He_i, d_He, p, cooling, cosmo, phys_const);
- n_eagle_cooling_rate_calls_1++;
-
- // Newton iterate the variable.
- x = x_old - relax_factor*(1.0 - u_ini*exp(-x_old) - LambdaNet*ratefact*dt*exp(-x_old))/(1.0 - dLambdaNet_du*ratefact*dt);
- // try limiting how far iterations can jump by
- if (*printflag == 1 && x < 0.3*x_old) {
- x = 0.3*x_old;
- } else if (*printflag == 1 && x > 1.3*x_old) {
- x = 1.3*x_old;
- }
+ {
+ x_old = x;
+ // this version is needed when reionization is included...
+ // LambdaNet = (LambdaTune / (dt * ratefact)) +
+ // eagle_cooling_rate_1d_table(exp(x_old), &dLambdaNet_du,
+ // H_plus_He_heat_table, H_plus_He_electron_abundance_table,
+ // element_cooling_table, element_electron_abundance_table, temp_table,
+ // z_index, dz, n_h_i, d_n_h, He_i, d_He, p, cooling, cosmo, phys_const);
+ LambdaNet = eagle_cooling_rate_1d_table(
+ exp(x_old), &dLambdaNet_du, H_plus_He_heat_table,
+ H_plus_He_electron_abundance_table, element_cooling_table,
+ element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h,
+ He_i, d_He, p, cooling, cosmo, phys_const);
+ n_eagle_cooling_rate_calls_1++;
+
+ // Newton iterate the variable.
+ x = x_old -
+ relax_factor * (1.0 - u_ini * exp(-x_old) -
+ LambdaNet * ratefact * dt * exp(-x_old)) /
+ (1.0 - dLambdaNet_du * ratefact * dt);
+ // try limiting how far iterations can jump by
+ if (*printflag == 1 && x < 0.3 * x_old) {
+ x = 0.3 * x_old;
+ } else if (*printflag == 1 && x > 1.3 * x_old) {
+ x = 1.3 * x_old;
+ }
+ // check whether iterations go out of bounds of table,
+ // can be used to try to prevent runaway
+ int table_out_of_bounds = 0;
+ // if (x > cooling->Therm[cooling->N_Temp - 1]*log(10) + 1) {
+ if (exp(x) > 1.0e19) {
+ // x = log(u_ini);
+ table_out_of_bounds = 1;
+ }
+ // if (x < cooling->Therm[0]*log(10) - 1) {
+ if (exp(x) < 1.0e9) {
+ // x = log(u_ini);
+ table_out_of_bounds = 1;
+ }
- // check whether iterations go out of bounds of table,
- // can be used to try to prevent runaway
- int table_out_of_bounds = 0;
- //if (x > cooling->Therm[cooling->N_Temp - 1]*log(10) + 1) {
- if (exp(x) > 1.0e19) {
- //x = log(u_ini);
- table_out_of_bounds = 1;
- }
- //if (x < cooling->Therm[0]*log(10) - 1) {
- if (exp(x) < 1.0e9) {
- //x = log(u_ini);
- table_out_of_bounds = 1;
- }
-
- residual = exp(x_old) - u_ini - LambdaNet*ratefact*dt ;
- if(dt > 0 && *printflag == 1) {
- //printf("z, n_h, x, x_old, u_ini, num, denom, LNet, dL_du %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %d %d\n",
- // cosmo->z, inn_h, exp(x), exp(x_old), u_ini,(1.0 - u_ini*exp(-x_old) - LambdaNet*ratefact*dt*exp(-x_old)),(1.0 - dLambdaNet_du*ratefact*dt), LambdaNet, dLambdaNet_du, x - x_old, i, table_out_of_bounds);
- table_out_of_bounds = 0;
- }
+ residual = exp(x_old) - u_ini - LambdaNet * ratefact * dt;
+ if (dt > 0 && *printflag == 1) {
+ // printf("z, n_h, x, x_old, u_ini, num, denom, LNet, dL_du %.3e %.3e
+ // %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %d %d\n",
+ // cosmo->z, inn_h, exp(x), exp(x_old), u_ini,(1.0 -
+ // u_ini*exp(-x_old) - LambdaNet*ratefact*dt*exp(-x_old)),(1.0 -
+ // dLambdaNet_du*ratefact*dt), LambdaNet, dLambdaNet_du, x - x_old, i,
+ // table_out_of_bounds);
+ table_out_of_bounds = 0;
+ }
- i++;
- if (table_out_of_bounds == 1) i = eagle_max_iterations;
- } while (fabs(x - x_old) > 1.0e-2 && i < eagle_max_iterations);
+ i++;
+ if (table_out_of_bounds == 1) i = eagle_max_iterations;
+ } while (fabs(x - x_old) > 1.0e-2 && i < eagle_max_iterations);
if (i >= eagle_max_iterations) {
n_eagle_cooling_rate_calls_3++;
// failed to converge the first time
@@ -1588,7 +1895,6 @@ __attribute__((always_inline)) INLINE static float newton_iter(float x_init,
}
return x;
-
}
/**
@@ -1603,16 +1909,17 @@ __attribute__((always_inline)) INLINE static float newton_iter(float x_init,
* @param dt The time-step of this particle.
*/
__attribute__((always_inline)) INLINE static void cooling_cool_part(
- const struct phys_const* restrict phys_const,
- const struct unit_system* restrict us,
- const struct cosmology* restrict cosmo,
- const struct cooling_function_data* restrict cooling,
- struct part* restrict p, struct xpart* restrict xp, float dt) {
-
- double u_old = hydro_get_physical_internal_energy(p,cosmo)*cooling->internal_energy_scale;
- float dz;
+ const struct phys_const *restrict phys_const,
+ const struct unit_system *restrict us,
+ const struct cosmology *restrict cosmo,
+ const struct cooling_function_data *restrict cooling,
+ struct part *restrict p, struct xpart *restrict xp, float dt) {
+
+ double u_old = hydro_get_physical_internal_energy(p, cosmo) *
+ cooling->internal_energy_scale;
+ float dz;
int z_index;
- get_redshift_index(cosmo->z,&z_index,&dz,cooling);
+ get_redshift_index(cosmo->z, &z_index, &dz, cooling);
float XH, HeFrac;
double inn_h;
@@ -1620,83 +1927,128 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
double ratefact, u, LambdaNet, LambdaTune = 0, dLambdaNet_du;
static double zold = 100, LambdaCumul = 0;
- dt *= units_cgs_conversion_factor(us,UNIT_CONV_TIME);
+ dt *= units_cgs_conversion_factor(us, UNIT_CONV_TIME);
u = u_old;
double u_ini = u_old, u_temp;
XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
- HeFrac = p->chemistry_data.metal_mass_fraction[chemistry_element_He] / (XH + p->chemistry_data.metal_mass_fraction[chemistry_element_He]);
- //printf("Eagle cooling.h density %.5e\n", hydro_get_physical_density(p,cosmo)*units_cgs_conversion_factor(us,UNIT_CONV_DENSITY));
+ HeFrac = p->chemistry_data.metal_mass_fraction[chemistry_element_He] /
+ (XH + p->chemistry_data.metal_mass_fraction[chemistry_element_He]);
+ // printf("Eagle cooling.h density %.5e\n",
+ // hydro_get_physical_density(p,cosmo)*units_cgs_conversion_factor(us,UNIT_CONV_DENSITY));
/* convert Hydrogen mass fraction in Hydrogen number density */
- inn_h = chemistry_get_number_density(p,cosmo,chemistry_element_H,phys_const)*cooling->number_density_scale;
- //float inn_h_eagle = hydro_get_physical_density(p,cosmo)*units_cgs_conversion_factor(us,UNIT_CONV_DENSITY) * XH / eagle_proton_mass_cgs;
+ inn_h =
+ chemistry_get_number_density(p, cosmo, chemistry_element_H, phys_const) *
+ cooling->number_density_scale;
+ // float inn_h_eagle =
+ // hydro_get_physical_density(p,cosmo)*units_cgs_conversion_factor(us,UNIT_CONV_DENSITY)
+ // * XH / eagle_proton_mass_cgs;
/* ratefact = inn_h * inn_h / rho; Might lead to round-off error: replaced by
* equivalent expression below */
ratefact = inn_h * (XH / eagle_proton_mass_cgs);
/* set helium and hydrogen reheating term */
- //LambdaTune = eagle_helium_reionization_extraheat(z_index, dz); // INCLUDE HELIUM REIONIZATION????
+ // LambdaTune = eagle_helium_reionization_extraheat(z_index, dz); // INCLUDE
+ // HELIUM REIONIZATION????
if (zold > z_index) {
LambdaCumul += LambdaTune;
- printf(" EagleDoCooling %d %g %g %g\n", z_index, dz, LambdaTune, LambdaCumul);
+ printf(" EagleDoCooling %d %g %g %g\n", z_index, dz, LambdaTune,
+ LambdaCumul);
zold = z_index;
}
-
+
int He_i, n_h_i;
float d_He, d_n_h;
get_index_1d(cooling->HeFrac, cooling->N_He, HeFrac, &He_i, &d_He);
get_index_1d(cooling->nH, cooling->N_nH, log10(inn_h), &n_h_i, &d_n_h);
- double temp_table[176]; // WARNING sort out how it is declared/allocated
- construct_1d_table_from_4d(p,cooling,cosmo,phys_const,cooling->table.element_cooling.temperature,z_index,dz,He_i,d_He,n_h_i,d_n_h,temp_table);
+ double temp_table[176]; // WARNING sort out how it is declared/allocated
+ construct_1d_table_from_4d(p, cooling, cosmo, phys_const,
+ cooling->table.element_cooling.temperature,
+ z_index, dz, He_i, d_He, n_h_i, d_n_h, temp_table);
// To be used when amount of cooling is small, put back in when ready
- LambdaNet = eagle_metal_cooling_rate(u_ini, temp_table, z_index, dz, n_h_i, d_n_h, He_i, d_He, p, cooling, cosmo, phys_const,NULL);
+ LambdaNet =
+ eagle_metal_cooling_rate(u_ini, temp_table, z_index, dz, n_h_i, d_n_h,
+ He_i, d_He, p, cooling, cosmo, phys_const, NULL);
if (fabs(ratefact * LambdaNet * dt) < 0.05 * u_old) {
/* cooling rate is small, take the explicit solution */
u = u_old + ratefact * LambdaNet * dt;
} else {
// construct 1d table of cooling rates wrt temperature
- double H_plus_He_heat_table[176]; // WARNING sort out how it is declared/allocated
- double H_plus_He_electron_abundance_table[176]; // WARNING sort out how it is declared/allocated
- double element_cooling_table[176]; // WARNING sort out how it is declared/allocated
- double element_electron_abundance_table[176]; // WARNING sort out how it is declared/allocated
- construct_1d_table_from_4d_H_He(p,cooling,cosmo,phys_const,cooling->table.element_cooling.H_plus_He_heating,z_index,dz,He_i,d_He,n_h_i,d_n_h,H_plus_He_heat_table, &u_temp, u_ini, dt);
- construct_1d_table_from_4d(p,cooling,cosmo,phys_const,cooling->table.element_cooling.H_plus_He_electron_abundance,z_index,dz,He_i,d_He,n_h_i,d_n_h,H_plus_He_electron_abundance_table);
- construct_1d_table_from_4d_elements(p,cooling,cosmo,phys_const,cooling->table.element_cooling.metal_heating,z_index,dz,n_h_i,d_n_h,element_cooling_table);
- construct_1d_table_from_3d(p,cooling,cosmo,phys_const,cooling->table.element_cooling.electron_abundance,z_index,dz,n_h_i,d_n_h,element_electron_abundance_table);
+ double H_plus_He_heat_table[176]; // WARNING sort out how it is
+ // declared/allocated
+ double H_plus_He_electron_abundance_table[176]; // WARNING sort out how it
+ // is declared/allocated
+ double element_cooling_table[176]; // WARNING sort out how it is
+ // declared/allocated
+ double element_electron_abundance_table[176]; // WARNING sort out how it is
+ // declared/allocated
+ construct_1d_table_from_4d_H_He(
+ p, cooling, cosmo, phys_const,
+ cooling->table.element_cooling.H_plus_He_heating, z_index, dz, He_i,
+ d_He, n_h_i, d_n_h, H_plus_He_heat_table, &u_temp, u_ini, dt);
+ construct_1d_table_from_4d(
+ p, cooling, cosmo, phys_const,
+ cooling->table.element_cooling.H_plus_He_electron_abundance, z_index,
+ dz, He_i, d_He, n_h_i, d_n_h, H_plus_He_electron_abundance_table);
+ construct_1d_table_from_4d_elements(
+ p, cooling, cosmo, phys_const,
+ cooling->table.element_cooling.metal_heating, z_index, dz, n_h_i, d_n_h,
+ element_cooling_table);
+ construct_1d_table_from_3d(
+ p, cooling, cosmo, phys_const,
+ cooling->table.element_cooling.electron_abundance, z_index, dz, n_h_i,
+ d_n_h, element_electron_abundance_table);
n_eagle_cooling_rate_calls_2++;
- LambdaNet = eagle_cooling_rate_1d_table(u_ini, &dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table, element_cooling_table, element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h, He_i, d_He, p, cooling, cosmo, phys_const);
+ LambdaNet = eagle_cooling_rate_1d_table(
+ u_ini, &dLambdaNet_du, H_plus_He_heat_table,
+ H_plus_He_electron_abundance_table, element_cooling_table,
+ element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h,
+ He_i, d_He, p, cooling, cosmo, phys_const);
double u_eq = 5.0e12;
u_temp = u_ini;
-
if (dt > 0) {
- float u_temp_guess = u_ini + LambdaNet*ratefact*dt;
- //printf("Eagle cooling.h id u_temp_guess, u_ini, LambdaNet, ratefact, dt %llu %.5e %.5e %.5e %.5e %.5e \n", p->id,u_temp_guess, u_ini, LambdaNet, ratefact, dt);
- //if ((LambdaNet < 0 && u_temp < u_temp_guess) ||
+ float u_temp_guess = u_ini + LambdaNet * ratefact * dt;
+ // printf("Eagle cooling.h id u_temp_guess, u_ini, LambdaNet, ratefact, dt
+ // %llu %.5e %.5e %.5e %.5e %.5e \n", p->id,u_temp_guess, u_ini,
+ // LambdaNet, ratefact, dt);
+ // if ((LambdaNet < 0 && u_temp < u_temp_guess) ||
// (LambdaNet >= 0 && u_temp > u_temp_guess))
// u_temp = u_temp_guess;
u_temp = u_temp_guess;
- if ((LambdaNet < 0 && u_temp < u_eq) ||
- (LambdaNet >= 0 && u_temp > u_eq))
+ if ((LambdaNet < 0 && u_temp < u_eq) || (LambdaNet >= 0 && u_temp > u_eq))
u_temp = u_eq;
if (isnan(u_temp)) u_temp = u_eq;
-
- //printf("Eagle cooling.h id, u_temp, n_h, z %llu %.5e %.5e %.5e\n", p->id, u_temp, inn_h, cosmo->z);
+
+ // printf("Eagle cooling.h id, u_temp, n_h, z %llu %.5e %.5e %.5e\n",
+ // p->id, u_temp, inn_h, cosmo->z);
int printflag = 0;
float relax_factor = 1.0;
- float x = newton_iter(log(u_temp),u_ini,H_plus_He_heat_table,H_plus_He_electron_abundance_table,element_cooling_table,element_electron_abundance_table,temp_table,z_index,dz,n_h_i,d_n_h,He_i,d_He,p,cosmo,cooling,phys_const,dt,&printflag,relax_factor);
+ float x =
+ newton_iter(log(u_temp), u_ini, H_plus_He_heat_table,
+ H_plus_He_electron_abundance_table, element_cooling_table,
+ element_electron_abundance_table, temp_table, z_index, dz,
+ n_h_i, d_n_h, He_i, d_He, p, cosmo, cooling, phys_const,
+ dt, &printflag, relax_factor);
if (printflag == 1) {
- // choose either a relax_factor < 1 to have better stability or bisection method
- //relax_factor = 0.75;
- x = newton_iter(log(u_temp),u_ini,H_plus_He_heat_table,H_plus_He_electron_abundance_table,element_cooling_table,element_electron_abundance_table,temp_table,z_index,dz,n_h_i,d_n_h,He_i,d_He,p,cosmo,cooling,phys_const,dt,&printflag,relax_factor);
- //x = bisection_iter(log(u_temp),u_ini,H_plus_He_heat_table,H_plus_He_electron_abundance_table,element_cooling_table,element_electron_abundance_table,temp_table,z_index,dz,n_h_i,d_n_h,He_i,d_He,p,cosmo,cooling,phys_const,dt);
+ // choose either a relax_factor < 1 to have better stability or
+ // bisection method
+ // relax_factor = 0.75;
+ x = newton_iter(log(u_temp), u_ini, H_plus_He_heat_table,
+ H_plus_He_electron_abundance_table,
+ element_cooling_table, element_electron_abundance_table,
+ temp_table, z_index, dz, n_h_i, d_n_h, He_i, d_He, p,
+ cosmo, cooling, phys_const, dt, &printflag,
+ relax_factor);
+ // x =
+ // bisection_iter(log(u_temp),u_ini,H_plus_He_heat_table,H_plus_He_electron_abundance_table,element_cooling_table,element_electron_abundance_table,temp_table,z_index,dz,n_h_i,d_n_h,He_i,d_He,p,cosmo,cooling,phys_const,dt);
n_eagle_cooling_rate_calls_4++;
printflag = 2;
}
@@ -1704,10 +2056,10 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
}
}
- // calculate du/dt
+ // calculate du/dt
float cooling_du_dt = 0.0;
- if (dt > 0){
- cooling_du_dt = (u - u_ini)/dt/cooling->power_scale;
+ if (dt > 0) {
+ cooling_du_dt = (u - u_ini) / dt / cooling->power_scale;
}
const float hydro_du_dt = hydro_get_internal_energy_dt(p);
@@ -1716,11 +2068,11 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
hydro_set_internal_energy_dt(p, hydro_du_dt + cooling_du_dt);
/* Store the radiated energy */
- xp->cooling_data.radiated_energy += -hydro_get_mass(p) * cooling_du_dt * (dt/units_cgs_conversion_factor(us,UNIT_CONV_TIME));
-
+ xp->cooling_data.radiated_energy +=
+ -hydro_get_mass(p) * cooling_du_dt *
+ (dt / units_cgs_conversion_factor(us, UNIT_CONV_TIME));
}
-
/**
* @brief Writes the current model of SPH to the file
* @param h_grpsph The HDF5 group in which to write
@@ -1741,15 +2093,15 @@ __attribute__((always_inline)) INLINE static void cooling_write_flavour(
* @param p Pointer to the particle data.
*/
__attribute__((always_inline)) INLINE static float cooling_timestep(
- const struct cooling_function_data* restrict cooling,
- const struct phys_const* restrict phys_const,
- const struct cosmology* restrict cosmo,
- const struct unit_system* restrict us, const struct part* restrict p) {
+ const struct cooling_function_data *restrict cooling,
+ const struct phys_const *restrict phys_const,
+ const struct cosmology *restrict cosmo,
+ const struct unit_system *restrict us, const struct part *restrict p) {
/* Remember to update when using an implicit integrator!!!*/
- //const float cooling_rate = cooling->cooling_rate;
- //const float internal_energy = hydro_get_comoving_internal_energy(p);
- //return cooling->cooling_tstep_mult * internal_energy / fabsf(cooling_rate);
+ // const float cooling_rate = cooling->cooling_rate;
+ // const float internal_energy = hydro_get_comoving_internal_energy(p);
+ // return cooling->cooling_tstep_mult * internal_energy / fabsf(cooling_rate);
return FLT_MAX;
}
@@ -1766,10 +2118,10 @@ __attribute__((always_inline)) INLINE static float cooling_timestep(
* @param xp Pointer to the extended particle data.
*/
__attribute__((always_inline)) INLINE static void cooling_first_init_part(
- const struct part* restrict p, struct xpart* restrict xp,
- const struct cooling_function_data* cooling) {
+ const struct part *restrict p, struct xpart *restrict xp,
+ const struct cooling_function_data *cooling) {
- xp->cooling_data.radiated_energy = 0.f;
+ xp->cooling_data.radiated_energy = 0.f;
}
/**
@@ -1778,7 +2130,7 @@ __attribute__((always_inline)) INLINE static void cooling_first_init_part(
* @param xp The extended particle data
*/
__attribute__((always_inline)) INLINE static float cooling_get_radiated_energy(
- const struct xpart* restrict xp) {
+ const struct xpart *restrict xp) {
return xp->cooling_data.radiated_energy;
}
@@ -1792,32 +2144,43 @@ __attribute__((always_inline)) INLINE static float cooling_get_radiated_energy(
* @param cooling The cooling properties to initialize
*/
static INLINE void cooling_init_backend(
- const struct swift_params* parameter_file, const struct unit_system* us,
- const struct phys_const* phys_const,
- struct cooling_function_data* cooling) {
-
+ const struct swift_params *parameter_file, const struct unit_system *us,
+ const struct phys_const *phys_const,
+ struct cooling_function_data *cooling) {
+
char fname[200];
- parser_get_param_string(parameter_file, "EagleCooling:filename",cooling->cooling_table_path);
- cooling->reionisation_redshift = parser_get_param_float(parameter_file, "EagleCooling:reionisation_redshift");
- cooling->calcium_over_silicon_ratio = parser_get_param_float(parameter_file, "EAGLEChemistry:CalciumOverSilicon");
- cooling->sulphur_over_silicon_ratio = parser_get_param_float(parameter_file, "EAGLEChemistry:SulphurOverSilicon");
+ parser_get_param_string(parameter_file, "EagleCooling:filename",
+ cooling->cooling_table_path);
+ cooling->reionisation_redshift = parser_get_param_float(
+ parameter_file, "EagleCooling:reionisation_redshift");
+ cooling->calcium_over_silicon_ratio = parser_get_param_float(
+ parameter_file, "EAGLEChemistry:CalciumOverSilicon");
+ cooling->sulphur_over_silicon_ratio = parser_get_param_float(
+ parameter_file, "EAGLEChemistry:SulphurOverSilicon");
GetCoolingRedshifts(cooling);
sprintf(fname, "%sz_0.000.hdf5", cooling->cooling_table_path);
- ReadCoolingHeader(fname,cooling);
+ ReadCoolingHeader(fname, cooling);
MakeNamePointers(cooling);
- cooling->table = eagle_readtable(cooling->cooling_table_path,cooling);
+ cooling->table = eagle_readtable(cooling->cooling_table_path, cooling);
printf("Eagle cooling.h read table \n");
- cooling->ElementAbundance_SOLARM1 = malloc(cooling->N_SolarAbundances*sizeof(float));
- cooling->solar_abundances = malloc(cooling->N_Elements*sizeof(double));
+ cooling->ElementAbundance_SOLARM1 =
+ malloc(cooling->N_SolarAbundances * sizeof(float));
+ cooling->solar_abundances = malloc(cooling->N_Elements * sizeof(double));
cooling->delta_u = cooling->Therm[1] - cooling->Therm[0];
- cooling->internal_energy_scale = units_cgs_conversion_factor(us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(us,UNIT_CONV_MASS);
- cooling->number_density_scale = units_cgs_conversion_factor(us,UNIT_CONV_DENSITY)/units_cgs_conversion_factor(us,UNIT_CONV_MASS);
- cooling->power_scale = units_cgs_conversion_factor(us,UNIT_CONV_POWER)/units_cgs_conversion_factor(us,UNIT_CONV_MASS);
- cooling->temperature_scale = units_cgs_conversion_factor(us,UNIT_CONV_TEMPERATURE);
+ cooling->internal_energy_scale =
+ units_cgs_conversion_factor(us, UNIT_CONV_ENERGY) /
+ units_cgs_conversion_factor(us, UNIT_CONV_MASS);
+ cooling->number_density_scale =
+ units_cgs_conversion_factor(us, UNIT_CONV_DENSITY) /
+ units_cgs_conversion_factor(us, UNIT_CONV_MASS);
+ cooling->power_scale = units_cgs_conversion_factor(us, UNIT_CONV_POWER) /
+ units_cgs_conversion_factor(us, UNIT_CONV_MASS);
+ cooling->temperature_scale =
+ units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
}
/**
@@ -1826,7 +2189,7 @@ static INLINE void cooling_init_backend(
* @param cooling The properties of the cooling function.
*/
static INLINE void cooling_print_backend(
- const struct cooling_function_data* cooling) {
+ const struct cooling_function_data *cooling) {
message("Cooling function is 'EAGLE'.");
}
diff --git a/src/cooling/EAGLE/cooling_read_table.h b/src/cooling/EAGLE/cooling_read_table.h
index 27c927dbf0d03ef2632677e472de13f4aa3d1e66..82bf92bdd54020bce12c9ffc8a42b19b0ea400fa 100644
--- a/src/cooling/EAGLE/cooling_read_table.h
+++ b/src/cooling/EAGLE/cooling_read_table.h
@@ -7,7 +7,8 @@
* @param fname Filepath
* @param cooling Cooling data structure
*/
-INLINE void ReadCoolingHeader(char *fname, struct cooling_function_data *cooling) {
+INLINE void ReadCoolingHeader(char *fname,
+ struct cooling_function_data *cooling) {
int i;
hid_t tempfile_id, dataset, datatype;
@@ -21,7 +22,8 @@ INLINE void ReadCoolingHeader(char *fname, struct cooling_function_data *cooling
error("[ReadCoolingHeader()]: unable to open file %s\n", fname);
}
- dataset = H5Dopen(tempfile_id, "/Header/Number_of_temperature_bins", H5P_DEFAULT);
+ dataset =
+ H5Dopen(tempfile_id, "/Header/Number_of_temperature_bins", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
&cooling->N_Temp);
status = H5Dclose(dataset);
@@ -31,12 +33,14 @@ INLINE void ReadCoolingHeader(char *fname, struct cooling_function_data *cooling
&cooling->N_nH);
status = H5Dclose(dataset);
- dataset = H5Dopen(tempfile_id, "/Header/Number_of_helium_fractions", H5P_DEFAULT);
+ dataset =
+ H5Dopen(tempfile_id, "/Header/Number_of_helium_fractions", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
&cooling->N_He);
status = H5Dclose(dataset);
- dataset = H5Dopen(tempfile_id, "/Header/Abundances/Number_of_abundances", H5P_DEFAULT);
+ dataset = H5Dopen(tempfile_id, "/Header/Abundances/Number_of_abundances",
+ H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
&cooling->N_SolarAbundances);
status = H5Dclose(dataset);
@@ -47,12 +51,12 @@ INLINE void ReadCoolingHeader(char *fname, struct cooling_function_data *cooling
status = H5Dclose(dataset);
/* allocate arrays for cooling table header */
- //allocate_header_arrays();
+ // allocate_header_arrays();
- cooling->Temp = malloc(cooling->N_Temp*sizeof(float));
- cooling->nH = malloc(cooling->N_Temp*sizeof(float));
- cooling->HeFrac = malloc(cooling->N_Temp*sizeof(float));
- cooling->SolarAbundances = malloc(cooling->N_Temp*sizeof(float));
+ cooling->Temp = malloc(cooling->N_Temp * sizeof(float));
+ cooling->nH = malloc(cooling->N_Temp * sizeof(float));
+ cooling->HeFrac = malloc(cooling->N_Temp * sizeof(float));
+ cooling->SolarAbundances = malloc(cooling->N_Temp * sizeof(float));
/* fill the arrays */
dataset = H5Dopen(tempfile_id, "/Solar/Temperature_bins", H5P_DEFAULT);
@@ -65,17 +69,20 @@ INLINE void ReadCoolingHeader(char *fname, struct cooling_function_data *cooling
cooling->nH);
status = H5Dclose(dataset);
- dataset = H5Dopen(tempfile_id, "/Metal_free/Helium_mass_fraction_bins", H5P_DEFAULT);
+ dataset = H5Dopen(tempfile_id, "/Metal_free/Helium_mass_fraction_bins",
+ H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
cooling->HeFrac);
status = H5Dclose(dataset);
- dataset = H5Dopen(tempfile_id, "/Header/Abundances/Solar_mass_fractions", H5P_DEFAULT);
+ dataset = H5Dopen(tempfile_id, "/Header/Abundances/Solar_mass_fractions",
+ H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
cooling->SolarAbundances);
status = H5Dclose(dataset);
- dataset = H5Dopen(tempfile_id, "/Metal_free/Temperature/Energy_density_bins", H5P_DEFAULT);
+ dataset = H5Dopen(tempfile_id, "/Metal_free/Temperature/Energy_density_bins",
+ H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
cooling->Therm);
status = H5Dclose(dataset);
@@ -94,20 +101,21 @@ INLINE void ReadCoolingHeader(char *fname, struct cooling_function_data *cooling
for (i = 0; i < cooling->N_Elements; i++)
cooling->ElementNames[i] = mystrdup(element_names[i]);
- //char solar_abund_names[cooling_N_SolarAbundances][EL_NAME_LENGTH];
+ // char solar_abund_names[cooling_N_SolarAbundances][EL_NAME_LENGTH];
/* assumed solar abundances used in constructing the tables, and corresponding
* names */
- //dataset = H5Dopen(tempfile_id, "/Header/Abundances/Abund_names", H5P_DEFAULT);
- //status = H5Dread(dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+ // dataset = H5Dopen(tempfile_id, "/Header/Abundances/Abund_names",
+ // H5P_DEFAULT);
+ // status = H5Dread(dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT,
// solar_abund_names);
- //status = H5Dclose(dataset);
- //H5Tclose(datatype);
+ // status = H5Dclose(dataset);
+ // H5Tclose(datatype);
- //for (i = 0; i < cooling_N_SolarAbundances; i++)
+ // for (i = 0; i < cooling_N_SolarAbundances; i++)
// cooling_SolarAbundanceNames[i] = mystrdup(solar_abund_names[i]);
- //status = H5Fclose(tempfile_id);
+ // status = H5Fclose(tempfile_id);
/* Convert to temperature, density and internal energy arrays to log10 */
for (i = 0; i < cooling->N_Temp; i++) {
@@ -119,7 +127,6 @@ INLINE void ReadCoolingHeader(char *fname, struct cooling_function_data *cooling
printf("Done with cooling table header.\n");
fflush(stdout);
-
}
#endif /* SWIFT_COOLING_EAGLE_READ_TABLES_H */
diff --git a/src/cooling/EAGLE/cooling_struct.h b/src/cooling/EAGLE/cooling_struct.h
index bde1cc0c27c944ac95fbe07abd680c4f4717c4ff..3dab0d08f32a71b1e0263cbb864b43123096320c 100644
--- a/src/cooling/EAGLE/cooling_struct.h
+++ b/src/cooling/EAGLE/cooling_struct.h
@@ -24,17 +24,15 @@
// EAGLE defined constants
#define eagle_element_name_length 20
#define eagle_cmb_temperature 2.728
-#define eagle_compton_rate 1.0178e-37*2.728*2.728*2.728*2.728
+#define eagle_compton_rate 1.0178e-37 * 2.728 * 2.728 * 2.728 * 2.728
#define eagle_metal_cooling_on 1
#define eagle_max_iterations 15
#define eagle_proton_mass_cgs 1.6726e-24
-
-
-/*
+/*
* @brief struct containing radiative and heating rates
*/
-//struct radiative_rates{
+// struct radiative_rates{
// float Cooling[NUMBER_OF_CTYPES];
// float Heating[NUMBER_OF_HTYPES];
// float TotalCoolingRate;
@@ -45,14 +43,17 @@
/*
* @brief struct containing cooling tables independent of redshift
- *
+ *
*/
struct cooling_tables_redshift_invariant {
- float *metal_heating; // size: [1][cooling_N_Temp][cooling_N_nH];
- float *H_plus_He_heating; // size: [1][cooling_N_He][cooling_N_Temp][cooling_N_nH];
- float *H_plus_He_electron_abundance; // size: [1][cooling_N_He][cooling_N_Temp][cooling_N_nH];
- float *temperature; // size: [1][cooling_N_He][cooling_N_Temp][cooling_N_nH];
- float *electron_abundance; // size: [1][cooling_N_Temp][cooling_N_nH];
+ float *metal_heating; // size: [1][cooling_N_Temp][cooling_N_nH];
+ float *H_plus_He_heating; // size:
+ // [1][cooling_N_He][cooling_N_Temp][cooling_N_nH];
+ float *
+ H_plus_He_electron_abundance; // size:
+ // [1][cooling_N_He][cooling_N_Temp][cooling_N_nH];
+ float *temperature; // size: [1][cooling_N_He][cooling_N_Temp][cooling_N_nH];
+ float *electron_abundance; // size: [1][cooling_N_Temp][cooling_N_nH];
};
/*
@@ -60,11 +61,20 @@ struct cooling_tables_redshift_invariant {
*
*/
struct cooling_tables {
- float *metal_heating; // size: [cooling_N_Redshifts][cooling_N_Temp][cooling_N_nH];
- float *H_plus_He_heating; // size: [cooling_N_Redshifts][cooling_N_He][cooling_N_Temp][cooling_N_nH];
- float *H_plus_He_electron_abundance; // size: [cooling_N_Redshifts][cooling_N_He][cooling_N_Temp][cooling_N_nH];
- float *temperature; // size: [cooling_N_Redshifts][cooling_N_He][cooling_N_Temp][cooling_N_nH];
- float *electron_abundance; // size: [cooling_N_Redshifts][cooling_N_Temp][cooling_N_nH];
+ float *metal_heating; // size:
+ // [cooling_N_Redshifts][cooling_N_Temp][cooling_N_nH];
+ float *
+ H_plus_He_heating; // size:
+ // [cooling_N_Redshifts][cooling_N_He][cooling_N_Temp][cooling_N_nH];
+ float *
+ H_plus_He_electron_abundance; // size:
+ // [cooling_N_Redshifts][cooling_N_He][cooling_N_Temp][cooling_N_nH];
+ float *
+ temperature; // size:
+ // [cooling_N_Redshifts][cooling_N_He][cooling_N_Temp][cooling_N_nH];
+ float *
+ electron_abundance; // size:
+ // [cooling_N_Redshifts][cooling_N_Temp][cooling_N_nH];
};
/*
@@ -122,7 +132,6 @@ struct cooling_function_data {
float sulphur_over_silicon_ratio;
double delta_u;
-
};
/**
diff --git a/src/cooling/EAGLE/eagle_cool_tables.h b/src/cooling/EAGLE/eagle_cool_tables.h
index 2d152732df29cbc2dd57805b6819a4dad718ecf1..aa1019b940e93461c929e146d91f9a916ba66df8 100644
--- a/src/cooling/EAGLE/eagle_cool_tables.h
+++ b/src/cooling/EAGLE/eagle_cool_tables.h
@@ -1,11 +1,11 @@
#ifndef SWIFT_EAGLE_COOL_TABLES_H
#define SWIFT_EAGLE_COOL_TABLES_H
-#include "cooling_struct.h"
-#include <stdlib.h>
-#include <string.h>
#include <hdf5.h>
#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+#include "cooling_struct.h"
#include "error.h"
/*
@@ -16,7 +16,7 @@
inline char *mystrdup(const char *s) {
char *p;
- p = (char *)malloc((strlen(s) + 1)*sizeof(char));
+ p = (char *)malloc((strlen(s) + 1) * sizeof(char));
strcpy(p, s);
return p;
}
@@ -43,8 +43,7 @@ inline void GetCoolingRedshifts(struct cooling_function_data *cooling) {
if (fscanf(infile, "%s", buffer) != EOF) {
cooling->N_Redshifts = atoi(buffer);
- cooling->Redshifts =
- (float *)malloc(cooling->N_Redshifts * sizeof(float));
+ cooling->Redshifts = (float *)malloc(cooling->N_Redshifts * sizeof(float));
while (fscanf(infile, "%s", buffer) != EOF) {
cooling->Redshifts[i] = atof(buffer);
@@ -52,7 +51,6 @@ inline void GetCoolingRedshifts(struct cooling_function_data *cooling) {
}
}
fclose(infile);
-
}
/*
@@ -61,7 +59,8 @@ inline void GetCoolingRedshifts(struct cooling_function_data *cooling) {
* @param fname Filepath
* @param cooling Cooling data structure
*/
-inline void ReadCoolingHeader(char *fname, struct cooling_function_data *cooling) {
+inline void ReadCoolingHeader(char *fname,
+ struct cooling_function_data *cooling) {
int i;
hid_t tempfile_id, dataset, datatype;
@@ -75,7 +74,8 @@ inline void ReadCoolingHeader(char *fname, struct cooling_function_data *cooling
error("[ReadCoolingHeader()]: unable to open file %s\n", fname);
}
- dataset = H5Dopen(tempfile_id, "/Header/Number_of_temperature_bins", H5P_DEFAULT);
+ dataset =
+ H5Dopen(tempfile_id, "/Header/Number_of_temperature_bins", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
&cooling->N_Temp);
status = H5Dclose(dataset);
@@ -85,12 +85,14 @@ inline void ReadCoolingHeader(char *fname, struct cooling_function_data *cooling
&cooling->N_nH);
status = H5Dclose(dataset);
- dataset = H5Dopen(tempfile_id, "/Header/Number_of_helium_fractions", H5P_DEFAULT);
+ dataset =
+ H5Dopen(tempfile_id, "/Header/Number_of_helium_fractions", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
&cooling->N_He);
status = H5Dclose(dataset);
- dataset = H5Dopen(tempfile_id, "/Header/Abundances/Number_of_abundances", H5P_DEFAULT);
+ dataset = H5Dopen(tempfile_id, "/Header/Abundances/Number_of_abundances",
+ H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
&cooling->N_SolarAbundances);
status = H5Dclose(dataset);
@@ -101,16 +103,18 @@ inline void ReadCoolingHeader(char *fname, struct cooling_function_data *cooling
status = H5Dclose(dataset);
/* allocate arrays for cooling table header */
- //allocate_header_arrays();
-
- cooling->Temp = malloc(cooling->N_Temp*sizeof(float));
- cooling->nH = malloc(cooling->N_nH*sizeof(float));
- cooling->HeFrac = malloc(cooling->N_He*sizeof(float));
- cooling->SolarAbundances = malloc(cooling->N_SolarAbundances*sizeof(float));
- cooling->Therm = malloc(cooling->N_Temp*sizeof(float));
- cooling->ElementNames = malloc(cooling->N_Elements*eagle_element_name_length*sizeof(char));
- cooling->SolarAbundanceNames = malloc(cooling->N_SolarAbundances*eagle_element_name_length*sizeof(char));
-
+ // allocate_header_arrays();
+
+ cooling->Temp = malloc(cooling->N_Temp * sizeof(float));
+ cooling->nH = malloc(cooling->N_nH * sizeof(float));
+ cooling->HeFrac = malloc(cooling->N_He * sizeof(float));
+ cooling->SolarAbundances = malloc(cooling->N_SolarAbundances * sizeof(float));
+ cooling->Therm = malloc(cooling->N_Temp * sizeof(float));
+ cooling->ElementNames =
+ malloc(cooling->N_Elements * eagle_element_name_length * sizeof(char));
+ cooling->SolarAbundanceNames = malloc(
+ cooling->N_SolarAbundances * eagle_element_name_length * sizeof(char));
+
/* fill the arrays */
dataset = H5Dopen(tempfile_id, "/Solar/Temperature_bins", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
@@ -122,17 +126,20 @@ inline void ReadCoolingHeader(char *fname, struct cooling_function_data *cooling
cooling->nH);
status = H5Dclose(dataset);
- dataset = H5Dopen(tempfile_id, "/Metal_free/Helium_mass_fraction_bins", H5P_DEFAULT);
+ dataset = H5Dopen(tempfile_id, "/Metal_free/Helium_mass_fraction_bins",
+ H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
cooling->HeFrac);
status = H5Dclose(dataset);
- dataset = H5Dopen(tempfile_id, "/Header/Abundances/Solar_mass_fractions", H5P_DEFAULT);
+ dataset = H5Dopen(tempfile_id, "/Header/Abundances/Solar_mass_fractions",
+ H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
cooling->SolarAbundances);
status = H5Dclose(dataset);
- dataset = H5Dopen(tempfile_id, "/Metal_free/Temperature/Energy_density_bins", H5P_DEFAULT);
+ dataset = H5Dopen(tempfile_id, "/Metal_free/Temperature/Energy_density_bins",
+ H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
cooling->Therm);
status = H5Dclose(dataset);
@@ -148,10 +155,10 @@ inline void ReadCoolingHeader(char *fname, struct cooling_function_data *cooling
H5Dread(dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, element_names);
status = H5Dclose(dataset);
H5Tclose(datatype);
-
+
for (i = 0; i < cooling->N_Elements; i++)
cooling->ElementNames[i] = mystrdup(element_names[i]);
-
+
char solar_abund_names[cooling->N_SolarAbundances][eagle_element_name_length];
/* assumed solar abundances used in constructing the tables, and corresponding
@@ -176,21 +183,26 @@ inline void ReadCoolingHeader(char *fname, struct cooling_function_data *cooling
}
for (i = 0; i < cooling->N_nH; i++) cooling->nH[i] = log10(cooling->nH[i]);
-
- //printf("eagle_cooling_tables.h temp max, min, N_Temp: %.5e, %.5e, %d\n",cooling->Temp[cooling->N_Temp-1], cooling->Temp[0],cooling->N_Temp);
- //printf("eagle_cooling_tables.h internal energy max, min, N_Temp: %.5e, %.5e, %d\n",cooling->Therm[cooling->N_Temp-1], cooling->Therm[0],cooling->N_Temp);
- //printf("eagle_cooling_tables.h H max, min, N_nH: %.5e, %.5e, %d\n",cooling->nH[cooling->N_nH-1], cooling->nH[0],cooling->N_nH);
- //printf("eagle_cooling_tables.h He max, min, N_He: %.5e, %.5e, %d\n",cooling->HeFrac[cooling->N_He-1], cooling->HeFrac[0],cooling->N_He);
- //printf("eagle_cooling_tables.h Solar abundances max, min, N_SolarAbundances: %.5e, %.5e, %d\n",cooling->SolarAbundances[cooling->N_SolarAbundances-1], cooling->SolarAbundances[0],cooling->N_SolarAbundances);
- //printf("eagle_cooling_tables.h N_Elements: %d\n",cooling->N_Elements);
+ // printf("eagle_cooling_tables.h temp max, min, N_Temp: %.5e, %.5e,
+ // %d\n",cooling->Temp[cooling->N_Temp-1], cooling->Temp[0],cooling->N_Temp);
+ // printf("eagle_cooling_tables.h internal energy max, min, N_Temp: %.5e,
+ // %.5e, %d\n",cooling->Therm[cooling->N_Temp-1],
+ // cooling->Therm[0],cooling->N_Temp);
+ // printf("eagle_cooling_tables.h H max, min, N_nH: %.5e, %.5e,
+ // %d\n",cooling->nH[cooling->N_nH-1], cooling->nH[0],cooling->N_nH);
+ // printf("eagle_cooling_tables.h He max, min, N_He: %.5e, %.5e,
+ // %d\n",cooling->HeFrac[cooling->N_He-1], cooling->HeFrac[0],cooling->N_He);
+ // printf("eagle_cooling_tables.h Solar abundances max, min,
+ // N_SolarAbundances: %.5e, %.5e,
+ // %d\n",cooling->SolarAbundances[cooling->N_SolarAbundances-1],
+ // cooling->SolarAbundances[0],cooling->N_SolarAbundances);
+ // printf("eagle_cooling_tables.h N_Elements: %d\n",cooling->N_Elements);
printf("Done with cooling table header.\n");
fflush(stdout);
-
}
-
/*
* @brief Get the cooling table for photoionized cooling (before redshift ~9)
*
@@ -198,7 +210,9 @@ inline void ReadCoolingHeader(char *fname, struct cooling_function_data *cooling
* @param cooling Cooling data structure
*/
-inline struct cooling_tables_redshift_invariant get_no_compt_table(char *cooling_table_path, const struct cooling_function_data *restrict cooling) {
+inline struct cooling_tables_redshift_invariant get_no_compt_table(
+ char *cooling_table_path,
+ const struct cooling_function_data *restrict cooling) {
struct cooling_tables_redshift_invariant cooling_table;
hid_t file_id, dataset;
@@ -215,24 +229,34 @@ inline struct cooling_tables_redshift_invariant get_no_compt_table(char *cooling
float *he_net_cooling_rate;
float *he_electron_abundance;
- net_cooling_rate = (float *)malloc(cooling->N_Temp*cooling->N_nH*sizeof(float));
- electron_abundance = (float *)malloc(cooling->N_Temp*cooling->N_nH*sizeof(float));
- temperature = (float *)malloc(cooling->N_He*cooling->N_Temp*cooling->N_nH*sizeof(float));
- he_net_cooling_rate = (float *)malloc(cooling->N_He*cooling->N_Temp*cooling->N_nH*sizeof(float));
- he_electron_abundance = (float *)malloc(cooling->N_He*cooling->N_Temp*cooling->N_nH*sizeof(float));
-
- cooling_table.metal_heating = (float *)malloc(cooling->N_Elements*cooling->N_Temp*cooling->N_nH*sizeof(float));
- cooling_table.electron_abundance = (float *)malloc(cooling->N_Temp*cooling->N_nH*sizeof(float));
- cooling_table.temperature = (float *)malloc(cooling->N_He*cooling->N_Temp*cooling->N_nH*sizeof(float));
- cooling_table.H_plus_He_heating = (float *)malloc(cooling->N_He*cooling->N_Temp*cooling->N_nH*sizeof(float));
- cooling_table.H_plus_He_electron_abundance = (float *)malloc(cooling->N_He*cooling->N_Temp*cooling->N_nH*sizeof(float));
+ net_cooling_rate =
+ (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
+ electron_abundance =
+ (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
+ temperature = (float *)malloc(cooling->N_He * cooling->N_Temp *
+ cooling->N_nH * sizeof(float));
+ he_net_cooling_rate = (float *)malloc(cooling->N_He * cooling->N_Temp *
+ cooling->N_nH * sizeof(float));
+ he_electron_abundance = (float *)malloc(cooling->N_He * cooling->N_Temp *
+ cooling->N_nH * sizeof(float));
+
+ cooling_table.metal_heating = (float *)malloc(
+ cooling->N_Elements * cooling->N_Temp * cooling->N_nH * sizeof(float));
+ cooling_table.electron_abundance =
+ (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
+ cooling_table.temperature = (float *)malloc(cooling->N_He * cooling->N_Temp *
+ cooling->N_nH * sizeof(float));
+ cooling_table.H_plus_He_heating = (float *)malloc(
+ cooling->N_He * cooling->N_Temp * cooling->N_nH * sizeof(float));
+ cooling_table.H_plus_He_electron_abundance = (float *)malloc(
+ cooling->N_He * cooling->N_Temp * cooling->N_nH * sizeof(float));
sprintf(fname, "%sz_8.989nocompton.hdf5", cooling_table_path);
file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
- //printf("GetNoCompTable Redshift 1 %ld %s\n", (long int)file_id, fname);
- //fflush(stdout);
+ // printf("GetNoCompTable Redshift 1 %ld %s\n", (long int)file_id, fname);
+ // fflush(stdout);
/* For normal elements */
for (specs = 0; specs < cooling->N_Elements; specs++) {
@@ -242,16 +266,18 @@ inline struct cooling_tables_redshift_invariant get_no_compt_table(char *cooling
net_cooling_rate);
status = H5Dclose(dataset);
- for (j = 0; j < cooling->N_Temp; j++){
- for (k = 0; k < cooling->N_nH; k++){
- table_index = row_major_index_2d(j,k,cooling->N_Temp,cooling->N_nH);
- cooling_index = row_major_index_4d(0,specs,k,j,1,cooling->N_Elements,cooling->N_nH,cooling->N_Temp); //Redshift invariant table!!!
- cooling_table.metal_heating[cooling_index] = -net_cooling_rate[table_index];
+ for (j = 0; j < cooling->N_Temp; j++) {
+ for (k = 0; k < cooling->N_nH; k++) {
+ table_index = row_major_index_2d(j, k, cooling->N_Temp, cooling->N_nH);
+ cooling_index = row_major_index_4d(
+ 0, specs, k, j, 1, cooling->N_Elements, cooling->N_nH,
+ cooling->N_Temp); // Redshift invariant table!!!
+ cooling_table.metal_heating[cooling_index] =
+ -net_cooling_rate[table_index];
}
}
}
-
/* Helium */
sprintf(set_name, "/Metal_free/Net_Cooling");
dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
@@ -271,15 +297,20 @@ inline struct cooling_tables_redshift_invariant get_no_compt_table(char *cooling
he_electron_abundance);
status = H5Dclose(dataset);
- for (i = 0; i < cooling->N_He; i++){
- for (j = 0; j < cooling->N_Temp; j++){
+ for (i = 0; i < cooling->N_He; i++) {
+ for (j = 0; j < cooling->N_Temp; j++) {
for (k = 0; k < cooling->N_nH; k++) {
- table_index = row_major_index_3d(i,j,k,cooling->N_He,cooling->N_Temp,cooling->N_nH);
- cooling_index = row_major_index_4d(0,i,k,j,1,cooling->N_He,cooling->N_nH,cooling->N_Temp); //Redshift invariant table!!!
- cooling_table.H_plus_He_heating[cooling_index] = -he_net_cooling_rate[table_index];
+ table_index = row_major_index_3d(i, j, k, cooling->N_He,
+ cooling->N_Temp, cooling->N_nH);
+ cooling_index =
+ row_major_index_4d(0, i, k, j, 1, cooling->N_He, cooling->N_nH,
+ cooling->N_Temp); // Redshift invariant table!!!
+ cooling_table.H_plus_He_heating[cooling_index] =
+ -he_net_cooling_rate[table_index];
cooling_table.H_plus_He_electron_abundance[cooling_index] =
he_electron_abundance[table_index];
- cooling_table.temperature[cooling_index] = log10(temperature[table_index]);
+ cooling_table.temperature[cooling_index] =
+ log10(temperature[table_index]);
}
}
}
@@ -290,21 +321,25 @@ inline struct cooling_tables_redshift_invariant get_no_compt_table(char *cooling
electron_abundance);
status = H5Dclose(dataset);
- for (i = 0; i < cooling->N_Temp; i++){
- for (j = 0; j < cooling->N_nH; j++){
- table_index = row_major_index_2d(i,j,cooling->N_Temp,cooling->N_nH);
- cooling_index = row_major_index_3d(0,j,i,1,cooling->N_nH,cooling->N_Temp); //Redshift invariant table!!!
- cooling_table.electron_abundance[cooling_index] = electron_abundance[table_index];
+ for (i = 0; i < cooling->N_Temp; i++) {
+ for (j = 0; j < cooling->N_nH; j++) {
+ table_index = row_major_index_2d(i, j, cooling->N_Temp, cooling->N_nH);
+ cooling_index =
+ row_major_index_3d(0, j, i, 1, cooling->N_nH,
+ cooling->N_Temp); // Redshift invariant table!!!
+ cooling_table.electron_abundance[cooling_index] =
+ electron_abundance[table_index];
}
}
status = H5Fclose(file_id);
- //cooling_table.metals_heating = cooling_MetalsNetHeating;
- //cooling_table.H_plus_He_heating = cooling_HplusHeNetHeating;
- //cooling_table.H_plus_He_electron_abundance = cooling_HplusHeElectronAbundance;
- //cooling_table.temperature = cooling_ThermalToTemp;
- //cooling_table.electron_abundance = cooling_SolarElectronAbundance;
+ // cooling_table.metals_heating = cooling_MetalsNetHeating;
+ // cooling_table.H_plus_He_heating = cooling_HplusHeNetHeating;
+ // cooling_table.H_plus_He_electron_abundance =
+ // cooling_HplusHeElectronAbundance;
+ // cooling_table.temperature = cooling_ThermalToTemp;
+ // cooling_table.electron_abundance = cooling_SolarElectronAbundance;
free(net_cooling_rate);
free(electron_abundance);
@@ -324,7 +359,9 @@ inline struct cooling_tables_redshift_invariant get_no_compt_table(char *cooling
* @param cooling Cooling data structure
*/
-inline struct cooling_tables_redshift_invariant get_collisional_table(char *cooling_table_path, const struct cooling_function_data* restrict cooling) {
+inline struct cooling_tables_redshift_invariant get_collisional_table(
+ char *cooling_table_path,
+ const struct cooling_function_data *restrict cooling) {
struct cooling_tables_redshift_invariant cooling_table;
hid_t file_id, dataset;
@@ -341,17 +378,25 @@ inline struct cooling_tables_redshift_invariant get_collisional_table(char *cool
float *he_net_cooling_rate;
float *he_electron_abundance;
- net_cooling_rate = (float *)malloc(cooling->N_Temp*sizeof(float));
- electron_abundance = (float *)malloc(cooling->N_Temp*sizeof(float));
- temperature = (float *)malloc(cooling->N_He*cooling->N_Temp*sizeof(float));
- he_net_cooling_rate = (float *)malloc(cooling->N_He*cooling->N_Temp*sizeof(float));
- he_electron_abundance = (float *)malloc(cooling->N_He*cooling->N_Temp*sizeof(float));
-
- cooling_table.metal_heating = (float *)malloc(cooling->N_Elements*cooling->N_Temp*sizeof(float));
- cooling_table.electron_abundance = (float *)malloc(cooling->N_Temp*sizeof(float));
- cooling_table.temperature = (float *)malloc(cooling->N_He*cooling->N_Temp*sizeof(float));
- cooling_table.H_plus_He_heating = (float *)malloc(cooling->N_He*cooling->N_Temp*sizeof(float));
- cooling_table.H_plus_He_electron_abundance = (float *)malloc(cooling->N_He*cooling->N_Temp*sizeof(float));
+ net_cooling_rate = (float *)malloc(cooling->N_Temp * sizeof(float));
+ electron_abundance = (float *)malloc(cooling->N_Temp * sizeof(float));
+ temperature =
+ (float *)malloc(cooling->N_He * cooling->N_Temp * sizeof(float));
+ he_net_cooling_rate =
+ (float *)malloc(cooling->N_He * cooling->N_Temp * sizeof(float));
+ he_electron_abundance =
+ (float *)malloc(cooling->N_He * cooling->N_Temp * sizeof(float));
+
+ cooling_table.metal_heating =
+ (float *)malloc(cooling->N_Elements * cooling->N_Temp * sizeof(float));
+ cooling_table.electron_abundance =
+ (float *)malloc(cooling->N_Temp * sizeof(float));
+ cooling_table.temperature =
+ (float *)malloc(cooling->N_He * cooling->N_Temp * sizeof(float));
+ cooling_table.H_plus_He_heating =
+ (float *)malloc(cooling->N_He * cooling->N_Temp * sizeof(float));
+ cooling_table.H_plus_He_electron_abundance =
+ (float *)malloc(cooling->N_He * cooling->N_Temp * sizeof(float));
sprintf(fname, "%sz_collis.hdf5", cooling_table_path);
@@ -365,9 +410,11 @@ inline struct cooling_tables_redshift_invariant get_collisional_table(char *cool
net_cooling_rate);
status = H5Dclose(dataset);
- for (j = 0; j < cooling->N_Temp; j++){
- cooling_index = row_major_index_3d(0,specs,j,1,cooling->N_Elements,cooling->N_Temp); //Redshift invariant table!!!
- cooling_table.metal_heating[cooling_index] = -net_cooling_rate[j];
+ for (j = 0; j < cooling->N_Temp; j++) {
+ cooling_index =
+ row_major_index_3d(0, specs, j, 1, cooling->N_Elements,
+ cooling->N_Temp); // Redshift invariant table!!!
+ cooling_table.metal_heating[cooling_index] = -net_cooling_rate[j];
}
}
@@ -390,14 +437,18 @@ inline struct cooling_tables_redshift_invariant get_collisional_table(char *cool
he_electron_abundance);
status = H5Dclose(dataset);
- for (i = 0; i < cooling->N_He; i++){
- for (j = 0; j < cooling->N_Temp; j++){
- table_index = row_major_index_2d(i,j,cooling->N_He,cooling->N_Temp);
- cooling_index = row_major_index_3d(0,i,j,1,cooling->N_He,cooling->N_Temp); //Redshift invariant table!!!
- cooling_table.H_plus_He_heating[cooling_index] = -he_net_cooling_rate[table_index];
+ for (i = 0; i < cooling->N_He; i++) {
+ for (j = 0; j < cooling->N_Temp; j++) {
+ table_index = row_major_index_2d(i, j, cooling->N_He, cooling->N_Temp);
+ cooling_index =
+ row_major_index_3d(0, i, j, 1, cooling->N_He,
+ cooling->N_Temp); // Redshift invariant table!!!
+ cooling_table.H_plus_He_heating[cooling_index] =
+ -he_net_cooling_rate[table_index];
cooling_table.H_plus_He_electron_abundance[cooling_index] =
he_electron_abundance[table_index];
- cooling_table.temperature[cooling_index] = log10(temperature[table_index]);
+ cooling_table.temperature[cooling_index] =
+ log10(temperature[table_index]);
}
}
@@ -407,8 +458,8 @@ inline struct cooling_tables_redshift_invariant get_collisional_table(char *cool
electron_abundance);
status = H5Dclose(dataset);
- for (i = 0; i < cooling->N_Temp; i++){
- cooling_table.electron_abundance[i] = electron_abundance[i];
+ for (i = 0; i < cooling->N_Temp; i++) {
+ cooling_table.electron_abundance[i] = electron_abundance[i];
}
status = H5Fclose(file_id);
@@ -424,13 +475,15 @@ inline struct cooling_tables_redshift_invariant get_collisional_table(char *cool
}
/*
- * @brief Get the cooling table for photodissociation
+ * @brief Get the cooling table for photodissociation
*
* @param cooling_table_path Filepath
* @param cooling Cooling data structure
*/
-inline struct cooling_tables_redshift_invariant get_photodis_table(char *cooling_table_path, const struct cooling_function_data* restrict cooling) {
+inline struct cooling_tables_redshift_invariant get_photodis_table(
+ char *cooling_table_path,
+ const struct cooling_function_data *restrict cooling) {
struct cooling_tables_redshift_invariant cooling_table;
hid_t file_id, dataset;
@@ -447,17 +500,27 @@ inline struct cooling_tables_redshift_invariant get_photodis_table(char *cooling
float *he_net_cooling_rate;
float *he_electron_abundance;
- net_cooling_rate = (float *)malloc(cooling->N_Temp*cooling->N_nH*sizeof(float));
- electron_abundance = (float *)malloc(cooling->N_Temp*cooling->N_nH*sizeof(float));
- temperature = (float *)malloc(cooling->N_He*cooling->N_Temp*cooling->N_nH*sizeof(float));
- he_net_cooling_rate = (float *)malloc(cooling->N_He*cooling->N_Temp*cooling->N_nH*sizeof(float));
- he_electron_abundance = (float *)malloc(cooling->N_He*cooling->N_Temp*cooling->N_nH*sizeof(float));
-
- cooling_table.metal_heating = (float *)malloc(cooling->N_Elements*cooling->N_Temp*cooling->N_nH*sizeof(float));
- cooling_table.electron_abundance = (float *)malloc(cooling->N_Temp*cooling->N_nH*sizeof(float));
- cooling_table.temperature = (float *)malloc(cooling->N_He*cooling->N_Temp*cooling->N_nH*sizeof(float));
- cooling_table.H_plus_He_heating = (float *)malloc(cooling->N_He*cooling->N_Temp*cooling->N_nH*sizeof(float));
- cooling_table.H_plus_He_electron_abundance = (float *)malloc(cooling->N_He*cooling->N_Temp*cooling->N_nH*sizeof(float));
+ net_cooling_rate =
+ (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
+ electron_abundance =
+ (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
+ temperature = (float *)malloc(cooling->N_He * cooling->N_Temp *
+ cooling->N_nH * sizeof(float));
+ he_net_cooling_rate = (float *)malloc(cooling->N_He * cooling->N_Temp *
+ cooling->N_nH * sizeof(float));
+ he_electron_abundance = (float *)malloc(cooling->N_He * cooling->N_Temp *
+ cooling->N_nH * sizeof(float));
+
+ cooling_table.metal_heating = (float *)malloc(
+ cooling->N_Elements * cooling->N_Temp * cooling->N_nH * sizeof(float));
+ cooling_table.electron_abundance =
+ (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
+ cooling_table.temperature = (float *)malloc(cooling->N_He * cooling->N_Temp *
+ cooling->N_nH * sizeof(float));
+ cooling_table.H_plus_He_heating = (float *)malloc(
+ cooling->N_He * cooling->N_Temp * cooling->N_nH * sizeof(float));
+ cooling_table.H_plus_He_electron_abundance = (float *)malloc(
+ cooling->N_He * cooling->N_Temp * cooling->N_nH * sizeof(float));
sprintf(fname, "%sz_photodis.hdf5", cooling_table_path);
@@ -471,11 +534,14 @@ inline struct cooling_tables_redshift_invariant get_photodis_table(char *cooling
net_cooling_rate);
status = H5Dclose(dataset);
- for (j = 0; j < cooling->N_Temp; j++){
- for (k = 0; k < cooling->N_nH; k++){
- table_index = row_major_index_2d(j,k,cooling->N_Temp,cooling->N_nH);
- cooling_index = row_major_index_3d(specs,j,k,cooling->N_Elements,cooling->N_Temp,cooling->N_nH); //Redshift invariant table!!!
- cooling_table.metal_heating[cooling_index] = -net_cooling_rate[table_index];
+ for (j = 0; j < cooling->N_Temp; j++) {
+ for (k = 0; k < cooling->N_nH; k++) {
+ table_index = row_major_index_2d(j, k, cooling->N_Temp, cooling->N_nH);
+ cooling_index = row_major_index_3d(
+ specs, j, k, cooling->N_Elements, cooling->N_Temp,
+ cooling->N_nH); // Redshift invariant table!!!
+ cooling_table.metal_heating[cooling_index] =
+ -net_cooling_rate[table_index];
}
}
}
@@ -499,15 +565,20 @@ inline struct cooling_tables_redshift_invariant get_photodis_table(char *cooling
he_electron_abundance);
status = H5Dclose(dataset);
- for (i = 0; i < cooling->N_He; i++){
- for (j = 0; j < cooling->N_Temp; j++){
+ for (i = 0; i < cooling->N_He; i++) {
+ for (j = 0; j < cooling->N_Temp; j++) {
for (k = 0; k < cooling->N_nH; k++) {
- table_index = row_major_index_3d(i,j,k,cooling->N_He,cooling->N_Temp,cooling->N_nH);
- cooling_index = row_major_index_3d(i,j,k,cooling->N_He,cooling->N_Temp,cooling->N_nH); //Redshift invariant table!!!
- cooling_table.H_plus_He_heating[cooling_index] = -he_net_cooling_rate[table_index];
+ table_index = row_major_index_3d(i, j, k, cooling->N_He,
+ cooling->N_Temp, cooling->N_nH);
+ cooling_index =
+ row_major_index_3d(i, j, k, cooling->N_He, cooling->N_Temp,
+ cooling->N_nH); // Redshift invariant table!!!
+ cooling_table.H_plus_He_heating[cooling_index] =
+ -he_net_cooling_rate[table_index];
cooling_table.H_plus_He_electron_abundance[cooling_index] =
he_electron_abundance[table_index];
- cooling_table.temperature[cooling_index] = log10(temperature[table_index]);
+ cooling_table.temperature[cooling_index] =
+ log10(temperature[table_index]);
}
}
}
@@ -518,11 +589,13 @@ inline struct cooling_tables_redshift_invariant get_photodis_table(char *cooling
electron_abundance);
status = H5Dclose(dataset);
- for (i = 0; i < cooling->N_Temp; i++){
- for (j = 0; j < cooling->N_nH; j++){
- table_index = row_major_index_2d(i,j,cooling->N_Temp,cooling->N_nH);
- cooling_index = row_major_index_2d(i,j,cooling->N_Temp,cooling->N_nH); //Redshift invariant table!!!
- cooling_table.electron_abundance[cooling_index] = electron_abundance[table_index];
+ for (i = 0; i < cooling->N_Temp; i++) {
+ for (j = 0; j < cooling->N_nH; j++) {
+ table_index = row_major_index_2d(i, j, cooling->N_Temp, cooling->N_nH);
+ cooling_index = row_major_index_2d(
+ i, j, cooling->N_Temp, cooling->N_nH); // Redshift invariant table!!!
+ cooling_table.electron_abundance[cooling_index] =
+ electron_abundance[table_index];
}
}
@@ -545,7 +618,9 @@ inline struct cooling_tables_redshift_invariant get_photodis_table(char *cooling
* @param cooling Cooling data structure
*/
-inline struct cooling_tables get_cooling_table(char *cooling_table_path, const struct cooling_function_data* restrict cooling) {
+inline struct cooling_tables get_cooling_table(
+ char *cooling_table_path,
+ const struct cooling_function_data *restrict cooling) {
struct cooling_tables cooling_table;
hid_t file_id, dataset;
@@ -562,21 +637,34 @@ inline struct cooling_tables get_cooling_table(char *cooling_table_path, const s
float *he_net_cooling_rate;
float *he_electron_abundance;
- net_cooling_rate = (float *)malloc(cooling->N_Temp*cooling->N_nH*sizeof(float));
- electron_abundance = (float *)malloc(cooling->N_Temp*cooling->N_nH*sizeof(float));
- temperature = (float *)malloc(cooling->N_He*cooling->N_Temp*cooling->N_nH*sizeof(float));
- he_net_cooling_rate = (float *)malloc(cooling->N_He*cooling->N_Temp*cooling->N_nH*sizeof(float));
- he_electron_abundance = (float *)malloc(cooling->N_He*cooling->N_Temp*cooling->N_nH*sizeof(float));
-
- cooling_table.metal_heating = (float *)malloc(cooling->N_Redshifts*cooling->N_Elements*cooling->N_Temp*cooling->N_nH*sizeof(float));
- cooling_table.electron_abundance = (float *)malloc(cooling->N_Redshifts*cooling->N_Temp*cooling->N_nH*sizeof(float));
- cooling_table.temperature = (float *)malloc(cooling->N_Redshifts*cooling->N_He*cooling->N_Temp*cooling->N_nH*sizeof(float));
- cooling_table.H_plus_He_heating = (float *)malloc(cooling->N_Redshifts*cooling->N_He*cooling->N_Temp*cooling->N_nH*sizeof(float));
- cooling_table.H_plus_He_electron_abundance = (float *)malloc(cooling->N_Redshifts*cooling->N_He*cooling->N_Temp*cooling->N_nH*sizeof(float));
-
+ net_cooling_rate =
+ (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
+ electron_abundance =
+ (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
+ temperature = (float *)malloc(cooling->N_He * cooling->N_Temp *
+ cooling->N_nH * sizeof(float));
+ he_net_cooling_rate = (float *)malloc(cooling->N_He * cooling->N_Temp *
+ cooling->N_nH * sizeof(float));
+ he_electron_abundance = (float *)malloc(cooling->N_He * cooling->N_Temp *
+ cooling->N_nH * sizeof(float));
+
+ cooling_table.metal_heating =
+ (float *)malloc(cooling->N_Redshifts * cooling->N_Elements *
+ cooling->N_Temp * cooling->N_nH * sizeof(float));
+ cooling_table.electron_abundance = (float *)malloc(
+ cooling->N_Redshifts * cooling->N_Temp * cooling->N_nH * sizeof(float));
+ cooling_table.temperature =
+ (float *)malloc(cooling->N_Redshifts * cooling->N_He * cooling->N_Temp *
+ cooling->N_nH * sizeof(float));
+ cooling_table.H_plus_He_heating =
+ (float *)malloc(cooling->N_Redshifts * cooling->N_He * cooling->N_Temp *
+ cooling->N_nH * sizeof(float));
+ cooling_table.H_plus_He_electron_abundance =
+ (float *)malloc(cooling->N_Redshifts * cooling->N_He * cooling->N_Temp *
+ cooling->N_nH * sizeof(float));
/* For normal elements */
- for (int z_index = 0; z_index < cooling->N_Redshifts; z_index++){
+ for (int z_index = 0; z_index < cooling->N_Redshifts; z_index++) {
sprintf(fname, "%sz_%1.3f.hdf5", cooling_table_path,
cooling->Redshifts[z_index]);
file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
@@ -592,12 +680,17 @@ inline struct cooling_tables get_cooling_table(char *cooling_table_path, const s
net_cooling_rate);
status = H5Dclose(dataset);
- for (i = 0; i < cooling->N_nH; i++){
- for (j = 0; j < cooling->N_Temp; j++){
- table_index = row_major_index_2d(j,i,cooling->N_Temp,cooling->N_nH);
- cooling_index = row_major_index_4d(z_index,specs,i,j,cooling->N_Redshifts,cooling->N_Elements,cooling->N_nH,cooling->N_Temp);
- //cooling_index = row_major_index_3d(specs,j,i,cooling->N_Elements,cooling->N_Temp,cooling->N_nH);
- cooling_table.metal_heating[cooling_index] = -net_cooling_rate[table_index];
+ for (i = 0; i < cooling->N_nH; i++) {
+ for (j = 0; j < cooling->N_Temp; j++) {
+ table_index =
+ row_major_index_2d(j, i, cooling->N_Temp, cooling->N_nH);
+ cooling_index = row_major_index_4d(
+ z_index, specs, i, j, cooling->N_Redshifts, cooling->N_Elements,
+ cooling->N_nH, cooling->N_Temp);
+ // cooling_index =
+ // row_major_index_3d(specs,j,i,cooling->N_Elements,cooling->N_Temp,cooling->N_nH);
+ cooling_table.metal_heating[cooling_index] =
+ -net_cooling_rate[table_index];
}
}
}
@@ -621,16 +714,22 @@ inline struct cooling_tables get_cooling_table(char *cooling_table_path, const s
he_electron_abundance);
status = H5Dclose(dataset);
- for (i = 0; i < cooling->N_He; i++){
- for (j = 0; j < cooling->N_Temp; j++){
+ for (i = 0; i < cooling->N_He; i++) {
+ for (j = 0; j < cooling->N_Temp; j++) {
for (k = 0; k < cooling->N_nH; k++) {
- table_index = row_major_index_3d(i,j,k,cooling->N_He,cooling->N_Temp,cooling->N_nH);
- cooling_index = row_major_index_4d(z_index,k,i,j,cooling->N_Redshifts,cooling->N_nH,cooling->N_He,cooling->N_Temp);
- //cooling_index = row_major_index_3d(i,j,k,cooling->N_He,cooling->N_Temp,cooling->N_nH);
- cooling_table.H_plus_He_heating[cooling_index] = -he_net_cooling_rate[table_index];
+ table_index = row_major_index_3d(i, j, k, cooling->N_He,
+ cooling->N_Temp, cooling->N_nH);
+ cooling_index =
+ row_major_index_4d(z_index, k, i, j, cooling->N_Redshifts,
+ cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ // cooling_index =
+ // row_major_index_3d(i,j,k,cooling->N_He,cooling->N_Temp,cooling->N_nH);
+ cooling_table.H_plus_He_heating[cooling_index] =
+ -he_net_cooling_rate[table_index];
cooling_table.H_plus_He_electron_abundance[cooling_index] =
he_electron_abundance[table_index];
- cooling_table.temperature[cooling_index] = log10(temperature[table_index]);
+ cooling_table.temperature[cooling_index] =
+ log10(temperature[table_index]);
}
}
}
@@ -641,12 +740,15 @@ inline struct cooling_tables get_cooling_table(char *cooling_table_path, const s
electron_abundance);
status = H5Dclose(dataset);
- for (i = 0; i < cooling->N_Temp; i++){
- for (j = 0; j < cooling->N_nH; j++){
- table_index = row_major_index_2d(i,j,cooling->N_Temp,cooling->N_nH);
- cooling_index = row_major_index_3d(z_index,j,i,cooling->N_Redshifts,cooling->N_nH,cooling->N_Temp);
- //cooling_index = row_major_index_2d(i,j,cooling->N_Temp,cooling->N_nH);
- cooling_table.electron_abundance[cooling_index] = electron_abundance[table_index];
+ for (i = 0; i < cooling->N_Temp; i++) {
+ for (j = 0; j < cooling->N_nH; j++) {
+ table_index = row_major_index_2d(i, j, cooling->N_Temp, cooling->N_nH);
+ cooling_index = row_major_index_3d(z_index, j, i, cooling->N_Redshifts,
+ cooling->N_nH, cooling->N_Temp);
+ // cooling_index =
+ // row_major_index_2d(i,j,cooling->N_Temp,cooling->N_nH);
+ cooling_table.electron_abundance[cooling_index] =
+ electron_abundance[table_index];
}
}
@@ -670,13 +772,18 @@ inline struct cooling_tables get_cooling_table(char *cooling_table_path, const s
* @param cooling_table_path Filepath
* @param cooling Cooling data structure
*/
-inline struct eagle_cooling_table eagle_readtable(char *cooling_table_path, const struct cooling_function_data* restrict cooling){
+inline struct eagle_cooling_table eagle_readtable(
+ char *cooling_table_path,
+ const struct cooling_function_data *restrict cooling) {
struct eagle_cooling_table table;
- table.photoionisation_cooling = get_no_compt_table(cooling_table_path, cooling);
- table.photodissociation_cooling = get_photodis_table(cooling_table_path, cooling);
- table.collisional_cooling = get_collisional_table(cooling_table_path, cooling);
+ table.photoionisation_cooling =
+ get_no_compt_table(cooling_table_path, cooling);
+ table.photodissociation_cooling =
+ get_photodis_table(cooling_table_path, cooling);
+ table.collisional_cooling =
+ get_collisional_table(cooling_table_path, cooling);
table.element_cooling = get_cooling_table(cooling_table_path, cooling);
return table;
@@ -688,7 +795,8 @@ inline struct eagle_cooling_table eagle_readtable(char *cooling_table_path, cons
* @param element_name Element string we want to match to element index
* @param cooling Cooling data structure
*/
-inline int element_index(char *element_name, const struct cooling_function_data* restrict cooling) {
+inline int element_index(char *element_name,
+ const struct cooling_function_data *restrict cooling) {
int i;
for (i = 0; i < cooling->N_Elements; i++)
@@ -720,12 +828,14 @@ inline int get_element_index(char *table[20], int size, char *element_name) {
*
* @param cooling Cooling data structure
*/
-inline void MakeNamePointers(struct cooling_function_data* cooling) {
+inline void MakeNamePointers(struct cooling_function_data *cooling) {
int i, j, sili_index = 0;
char ElementNames[cooling->N_Elements][eagle_element_name_length];
- /* This is ridiculous, way too many element name arrays. Needs to be changed */
- //ElementNames = malloc(cooling->N_Elements*eagle_element_name_length*sizeof(char));
+ /* This is ridiculous, way too many element name arrays. Needs to be changed
+ */
+ // ElementNames =
+ // malloc(cooling->N_Elements*eagle_element_name_length*sizeof(char));
strcpy(ElementNames[0], "Hydrogen");
strcpy(ElementNames[1], "Helium");
strcpy(ElementNames[2], "Carbon");
@@ -737,7 +847,8 @@ inline void MakeNamePointers(struct cooling_function_data* cooling) {
strcpy(ElementNames[8], "Iron");
cooling->ElementNamePointers = malloc(cooling->N_Elements * sizeof(int));
- cooling->SolarAbundanceNamePointers = malloc(cooling->N_Elements * sizeof(int));
+ cooling->SolarAbundanceNamePointers =
+ malloc(cooling->N_Elements * sizeof(int));
for (i = 0; i < cooling->N_Elements; i++) {
if (strcmp(ElementNames[i], "Silicon") == 0) sili_index = i;
@@ -748,7 +859,8 @@ inline void MakeNamePointers(struct cooling_function_data* cooling) {
cooling->ElementNamePointers[i] = -999;
for (j = 0; j < cooling->N_SolarAbundances; j++) {
- if (strcmp(cooling->ElementNames[i], cooling->SolarAbundanceNames[j]) == 0)
+ if (strcmp(cooling->ElementNames[i], cooling->SolarAbundanceNames[j]) ==
+ 0)
cooling->SolarAbundanceNamePointers[i] = j;
}
diff --git a/src/engine.c b/src/engine.c
index 8f8c17ed588abd23190e6ccf7f9cc559af3b0fc8..d9312ae65e4b18fab1db80c711bc61e645e744f3 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -4743,17 +4743,37 @@ void engine_step(struct engine *e) {
/* Print some information to the screen */
printf(
- " %6d %14e %14e %10.5f %14e %4d %4d %12lld %12lld %12lld %21.3f %6d\n",// %.5e %6d %6d %.5e %6d\n",
+ " %6d %14e %14e %10.5f %14e %4d %4d %12lld %12lld %12lld %21.3f %6d\n", // %.5e %6d %6d %.5e %6d\n",
e->step, e->time, e->cosmology->a, e->cosmology->z, e->time_step,
e->min_active_bin, e->max_active_bin, e->updates, e->g_updates,
- e->s_updates, e->wallclock_time, e->step_props);//, ((float) n_eagle_cooling_rate_calls_1)/((float) n_eagle_cooling_rate_calls_2), n_eagle_cooling_rate_calls_2, n_eagle_cooling_rate_calls_3, ((float) n_eagle_cooling_rate_calls_3)/((float) n_eagle_cooling_rate_calls_2), n_eagle_cooling_rate_calls_4);
+ e->s_updates, e->wallclock_time,
+ e->step_props); //, ((float) n_eagle_cooling_rate_calls_1)/((float)
+ //n_eagle_cooling_rate_calls_2),
+ //n_eagle_cooling_rate_calls_2,
+ //n_eagle_cooling_rate_calls_3, ((float)
+ //n_eagle_cooling_rate_calls_3)/((float)
+ //n_eagle_cooling_rate_calls_2),
+ //n_eagle_cooling_rate_calls_4);
fflush(stdout);
- fprintf(e->file_timesteps,
- " %6d %14e %14e %14e %4d %4d %12lld %12lld %12lld %21.3f %6d\n",// %.5e %6d %6d %.5e %6d\n",
- e->step, e->time, e->cosmology->a, e->time_step, e->min_active_bin,
- e->max_active_bin, e->updates, e->g_updates, e->s_updates,
- e->wallclock_time, e->step_props, );//, ((float) (n_eagle_cooling_rate_calls_1 - 20*n_eagle_cooling_rate_calls_4))/((float) n_eagle_cooling_rate_calls_2), n_eagle_cooling_rate_calls_2, n_eagle_cooling_rate_calls_3, ((float) n_eagle_cooling_rate_calls_3)/((float) n_eagle_cooling_rate_calls_2), n_eagle_cooling_rate_calls_4);
+ fprintf(
+ e->file_timesteps,
+ " %6d %14e %14e %14e %4d %4d %12lld %12lld %12lld %21.3f %6d\n", // %.5e
+ // %6d
+ // %6d
+ // %.5e
+ // %6d\n",
+ e->step, e->time, e->cosmology->a, e->time_step, e->min_active_bin,
+ e->max_active_bin, e->updates, e->g_updates, e->s_updates,
+ e->wallclock_time,
+ e->step_props, ); //, ((float) (n_eagle_cooling_rate_calls_1 -
+ //20*n_eagle_cooling_rate_calls_4))/((float)
+ //n_eagle_cooling_rate_calls_2),
+ //n_eagle_cooling_rate_calls_2,
+ //n_eagle_cooling_rate_calls_3, ((float)
+ //n_eagle_cooling_rate_calls_3)/((float)
+ //n_eagle_cooling_rate_calls_2),
+ //n_eagle_cooling_rate_calls_4);
fflush(e->file_timesteps);
}
n_eagle_cooling_rate_calls_1 = 0;
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index b7a70e9bd2d5bdae8e97fc770f14bf9df90a37de..caf5cdddaa6ae7dfd107fe35e2d2d3453cc3460e 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -62,9 +62,11 @@ __attribute__((always_inline)) INLINE static float
hydro_get_physical_internal_energy(const struct part *restrict p,
const struct cosmology *cosmo) {
- float u = gas_internal_energy_from_entropy(p->rho * cosmo->a3_inv, p->entropy);
- if (u < 0){
- printf("Gadget 2 hydro.h u, rho, a3_inv, entropy %.5e %.5e %.5e %.5e \n", u,p->rho,cosmo->a3_inv,p->entropy);
+ float u =
+ gas_internal_energy_from_entropy(p->rho * cosmo->a3_inv, p->entropy);
+ if (u < 0) {
+ printf("Gadget 2 hydro.h u, rho, a3_inv, entropy %.5e %.5e %.5e %.5e \n", u,
+ p->rho, cosmo->a3_inv, p->entropy);
fflush(stdout);
}
return gas_internal_energy_from_entropy(p->rho * cosmo->a3_inv, p->entropy);
@@ -253,7 +255,12 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
/* CFL condition */
const float dt_cfl = 2.f * kernel_gamma * CFL * cosmo->a * p->h /
(cosmo->a_factor_sound_speed * p->force.v_sig);
- if (dt_cfl == 0.0) printf ("Gadget2 hydro.h id, kernel, cfl, scale factor, h, sound speed, v_sig, u %llu %.5e %.5e %.5e %.5e %.5e %.5e %.5e\n", p->id, kernel_gamma, CFL, cosmo->a, p->h, cosmo->a_factor_sound_speed, p->force.v_sig,hydro_get_physical_internal_energy(p,cosmo));
+ if (dt_cfl == 0.0)
+ printf(
+ "Gadget2 hydro.h id, kernel, cfl, scale factor, h, sound speed, v_sig, "
+ "u %llu %.5e %.5e %.5e %.5e %.5e %.5e %.5e\n",
+ p->id, kernel_gamma, CFL, cosmo->a, p->h, cosmo->a_factor_sound_speed,
+ p->force.v_sig, hydro_get_physical_internal_energy(p, cosmo));
return dt_cfl;
}
@@ -493,8 +500,11 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Predict the entropy */
float entropy_old = p->entropy;
p->entropy += p->entropy_dt * dt_therm;
- if (p->entropy < 0){
- printf("Gadget2 hydro.h entropy less than zero old entropy, entropy, entropy_dt, dt_therm %.5e %.5e %.5e %.5e\n", entropy_old, p->entropy, p->entropy_dt,dt_therm);
+ if (p->entropy < 0) {
+ printf(
+ "Gadget2 hydro.h entropy less than zero old entropy, entropy, "
+ "entropy_dt, dt_therm %.5e %.5e %.5e %.5e\n",
+ entropy_old, p->entropy, p->entropy_dt, dt_therm);
fflush(stdout);
}
diff --git a/src/runner.h b/src/runner.h
index 3d96e7718f9594cee78228c64efaa9654ed0d29e..e33a3e380e6097a67258d116d617483caca35086 100644
--- a/src/runner.h
+++ b/src/runner.h
@@ -86,6 +86,4 @@ void runner_do_unskip_mapper(void *map_data, int num_elements,
void runner_do_drift_all_mapper(void *map_data, int num_elements,
void *extra_data);
-
-
#endif /* SWIFT_RUNNER_H */
diff --git a/tests/testCooling.c b/tests/testCooling.c
index b0ba3a511cd8268b0ca25275fdfa76f6ae850d2c..488e0db0d67dbd28f2138c77fbb4b3854398f825 100644
--- a/tests/testCooling.c
+++ b/tests/testCooling.c
@@ -17,11 +17,11 @@
*
******************************************************************************/
-#include "swift.h"
#include "cooling.h"
+#include "hydro.h"
#include "physical_constants.h"
+#include "swift.h"
#include "units.h"
-#include "hydro.h"
int main() {
@@ -35,9 +35,17 @@ int main() {
struct cooling_function_data cooling;
struct cosmology cosmo;
char *parametersFileName = "../examples/CoolingBox/coolingBox.yml";
- enum table_index {EAGLE_Hydrogen=0,EAGLE_Helium,EAGLE_Carbon,EAGLE_Nitrogen,
- EAGLE_Oxygen,EAGLE_Neon,EAGLE_Magnesium,EAGLE_Silicon,
- EAGLE_Iron};
+ enum table_index {
+ EAGLE_Hydrogen = 0,
+ EAGLE_Helium,
+ EAGLE_Carbon,
+ EAGLE_Nitrogen,
+ EAGLE_Oxygen,
+ EAGLE_Neon,
+ EAGLE_Magnesium,
+ EAGLE_Silicon,
+ EAGLE_Iron
+ };
if (params == NULL) error("Error allocating memory for the parameter file.");
message("Reading runtime parameters from file '%s'", parametersFileName);
@@ -50,64 +58,85 @@ int main() {
phys_const_init(&us, params, &internal_const);
double hydrogen_number_density_cgs = 1e-1;
- double u,u_cgs, hydrogen_number_density, pressure, gamma, cooling_du_dt, temperature_cgs, newton_func, u_ini_cgs, ratefact, dt_cgs;
- u_cgs = 0; cooling_du_dt = 0; temperature_cgs = 0; newton_func = 0;
- //int n_t_i = 2000;
- dt_cgs = 4.0e-8*units_cgs_conversion_factor(&us,UNIT_CONV_TIME);
+ double u, u_cgs, hydrogen_number_density, pressure, gamma, cooling_du_dt,
+ temperature_cgs, newton_func, u_ini_cgs, ratefact, dt_cgs;
+ u_cgs = 0;
+ cooling_du_dt = 0;
+ temperature_cgs = 0;
+ newton_func = 0;
+ // int n_t_i = 2000;
+ dt_cgs = 4.0e-8 * units_cgs_conversion_factor(&us, UNIT_CONV_TIME);
char *output_filename = "cooling_output.dat";
FILE *output_file = fopen(output_filename, "w");
- if (output_file == NULL)
- {
- printf("Error opening file!\n");
- exit(1);
+ if (output_file == NULL) {
+ printf("Error opening file!\n");
+ exit(1);
}
char *output_filename2 = "temperature_output.dat";
FILE *output_file2 = fopen(output_filename2, "w");
- if (output_file2 == NULL)
- {
- printf("Error opening file!\n");
- exit(1);
+ if (output_file2 == NULL) {
+ printf("Error opening file!\n");
+ exit(1);
}
char *output_filename3 = "newton_output.dat";
FILE *output_file3 = fopen(output_filename3, "w");
- if (output_file3 == NULL)
- {
- printf("Error opening file!\n");
- exit(1);
+ if (output_file3 == NULL) {
+ printf("Error opening file!\n");
+ exit(1);
}
- gamma = 5.0/3.0;
+ gamma = 5.0 / 3.0;
chemistry_init(params, &us, &internal_const, &chemistry_data);
- chemistry_first_init_part(&p,&xp,&chemistry_data);
+ chemistry_first_init_part(&p, &xp, &chemistry_data);
chemistry_print(&chemistry_data);
-
+
cosmology_init(params, &us, &internal_const, &cosmo);
cosmology_print(&cosmo);
-
- u = 1.0*pow(10.0,11)/(units_cgs_conversion_factor(&us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS));
- pressure = u*p.rho*(gamma -1.0);
- hydrogen_number_density = hydrogen_number_density_cgs*pow(units_cgs_conversion_factor(&us,UNIT_CONV_LENGTH),3);
- p.rho = hydrogen_number_density*internal_const.const_proton_mass*(1.0+p.chemistry_data.metal_mass_fraction[EAGLE_Helium]/p.chemistry_data.metal_mass_fraction[EAGLE_Hydrogen]);
- p.entropy = pressure/(pow(p.rho,gamma));
-
+ u = 1.0 * pow(10.0, 11) /
+ (units_cgs_conversion_factor(&us, UNIT_CONV_ENERGY) /
+ units_cgs_conversion_factor(&us, UNIT_CONV_MASS));
+ pressure = u * p.rho * (gamma - 1.0);
+ hydrogen_number_density =
+ hydrogen_number_density_cgs *
+ pow(units_cgs_conversion_factor(&us, UNIT_CONV_LENGTH), 3);
+ p.rho = hydrogen_number_density * internal_const.const_proton_mass *
+ (1.0 +
+ p.chemistry_data.metal_mass_fraction[EAGLE_Helium] /
+ p.chemistry_data.metal_mass_fraction[EAGLE_Hydrogen]);
+ p.entropy = pressure / (pow(p.rho, gamma));
cooling_init(params, &us, &internal_const, &cooling);
cooling_print(&cooling);
// construct 1d table of cooling rates wrt temperature
- float H_plus_He_heat_table[176]; // WARNING sort out how it is declared/allocated
- float H_plus_He_electron_abundance_table[176]; // WARNING sort out how it is declared/allocated
- float temp_table[176]; // WARNING sort out how it is declared/allocated
- float element_cooling_table[9*176]; // WARNING sort out how it is declared/allocated
- float element_electron_abundance_table[176]; // WARNING sort out how it is declared/allocated
- construct_1d_table_from_4d(&p,&cooling,&cosmo,&internal_const,cooling.table.element_cooling.H_plus_He_heating,H_plus_He_heat_table);
- construct_1d_table_from_4d(&p,&cooling,&cosmo,&internal_const,cooling.table.element_cooling.H_plus_He_electron_abundance,H_plus_He_electron_abundance_table);
- construct_1d_table_from_4d(&p,&cooling,&cosmo,&internal_const,cooling.table.element_cooling.temperature,temp_table);
- construct_1d_table_from_4d_elements(&p,&cooling,&cosmo,&internal_const,cooling.table.element_cooling.metal_heating,element_cooling_table);
- construct_1d_table_from_3d(&p,&cooling,&cosmo,&internal_const,cooling.table.element_cooling.electron_abundance,element_electron_abundance_table);
+ float H_plus_He_heat_table[176]; // WARNING sort out how it is
+ // declared/allocated
+ float H_plus_He_electron_abundance_table[176]; // WARNING sort out how it is
+ // declared/allocated
+ float temp_table[176]; // WARNING sort out how it is declared/allocated
+ float element_cooling_table[9 * 176]; // WARNING sort out how it is
+ // declared/allocated
+ float element_electron_abundance_table[176]; // WARNING sort out how it is
+ // declared/allocated
+ construct_1d_table_from_4d(&p, &cooling, &cosmo, &internal_const,
+ cooling.table.element_cooling.H_plus_He_heating,
+ H_plus_He_heat_table);
+ construct_1d_table_from_4d(
+ &p, &cooling, &cosmo, &internal_const,
+ cooling.table.element_cooling.H_plus_He_electron_abundance,
+ H_plus_He_electron_abundance_table);
+ construct_1d_table_from_4d(&p, &cooling, &cosmo, &internal_const,
+ cooling.table.element_cooling.temperature,
+ temp_table);
+ construct_1d_table_from_4d_elements(
+ &p, &cooling, &cosmo, &internal_const,
+ cooling.table.element_cooling.metal_heating, element_cooling_table);
+ construct_1d_table_from_3d(&p, &cooling, &cosmo, &internal_const,
+ cooling.table.element_cooling.electron_abundance,
+ element_electron_abundance_table);
//
// printf("cooling table values \n");
@@ -115,55 +144,87 @@ int main() {
// printf(" %.5e", H_plus_He_heat_table[j]+element_cooling_table[j]
float XH = p.chemistry_data.metal_mass_fraction[chemistry_element_H];
- float inn_h = chemistry_get_number_density(&p,&cosmo,chemistry_element_H,&internal_const)*cooling.number_density_scale;
+ float inn_h = chemistry_get_number_density(&p, &cosmo, chemistry_element_H,
+ &internal_const) *
+ cooling.number_density_scale;
ratefact = inn_h * (XH / eagle_proton_mass_cgs);
- u_ini_cgs = pow(10.0,14);
+ u_ini_cgs = pow(10.0, 14);
// Compute contributions to cooling rate from different metals
- //for(int t_i = 0; t_i < n_t_i; t_i++){
- //
+ // for(int t_i = 0; t_i < n_t_i; t_i++){
+ //
// u_cgs = pow(10.0,11 + t_i*6.0/n_t_i);
- // u = u_cgs/(units_cgs_conversion_factor(&us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS));
+ // u =
+ // u_cgs/(units_cgs_conversion_factor(&us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS));
// pressure = u*p.rho*(gamma -1.0);
- // hydrogen_number_density = hydrogen_number_density_cgs*pow(units_cgs_conversion_factor(&us,UNIT_CONV_LENGTH),3);
- // p.rho = hydrogen_number_density*internal_const.const_proton_mass*(1.0+p.chemistry_data.metal_mass_fraction[EAGLE_Helium]/p.chemistry_data.metal_mass_fraction[EAGLE_Hydrogen]);
+ // hydrogen_number_density =
+ // hydrogen_number_density_cgs*pow(units_cgs_conversion_factor(&us,UNIT_CONV_LENGTH),3);
+ // p.rho =
+ // hydrogen_number_density*internal_const.const_proton_mass*(1.0+p.chemistry_data.metal_mass_fraction[EAGLE_Helium]/p.chemistry_data.metal_mass_fraction[EAGLE_Hydrogen]);
// p.entropy = pressure/(pow(p.rho,gamma));
- // //cooling_du_dt = eagle_print_metal_cooling_rate(&p,&cooling,&cosmo,&internal_const);
- // cooling_du_dt = eagle_print_metal_cooling_rate_1d_table(H_plus_He_heat_table,H_plus_He_electron_abundance_table,element_cooling_table,element_electron_abundance_table,temp_table,&p,&cooling,&cosmo,&internal_const);
- // temperature_cgs = eagle_convert_u_to_temp(u_cgs,&p,&cooling,&cosmo,&internal_const);
+ // //cooling_du_dt =
+ // eagle_print_metal_cooling_rate(&p,&cooling,&cosmo,&internal_const);
+ // cooling_du_dt =
+ // eagle_print_metal_cooling_rate_1d_table(H_plus_He_heat_table,H_plus_He_electron_abundance_table,element_cooling_table,element_electron_abundance_table,temp_table,&p,&cooling,&cosmo,&internal_const);
+ // temperature_cgs =
+ // eagle_convert_u_to_temp(u_cgs,&p,&cooling,&cosmo,&internal_const);
// fprintf(output_file,"%.5e %.5e\n",temperature_cgs,cooling_du_dt);
- // fprintf(output_file2,"%.5e %.5e\n",u*(units_cgs_conversion_factor(&us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS)), temperature_cgs);
+ // fprintf(output_file2,"%.5e
+ // %.5e\n",u*(units_cgs_conversion_factor(&us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS)),
+ // temperature_cgs);
// newton_func = u_cgs - u_ini_cgs - cooling_du_dt*ratefact*dt_cgs;
// fprintf(output_file3,"%.5e %.5e\n",u_cgs,newton_func);
- // //printf("u,u_ini,cooling_du_dt,ratefact,dt %.5e %.5e %.5e %.5e %.5e \n",u_cgs,u_ini_cgs,cooling_du_dt,ratefact,dt_cgs);
+ // //printf("u,u_ini,cooling_du_dt,ratefact,dt %.5e %.5e %.5e %.5e %.5e
+ // \n",u_cgs,u_ini_cgs,cooling_du_dt,ratefact,dt_cgs);
//}
- //fclose(output_file);
- //fclose(output_file2);
- //fclose(output_file3);
+ // fclose(output_file);
+ // fclose(output_file2);
+ // fclose(output_file3);
double dLambdaNet_du, LambdaNet; //, LambdaNext;
- float x_init, u_eq = 2.0e12 ;
- for(int j = 5; j < 6; j++){
- float u_ini = eagle_convert_temp_to_u_1d_table(pow(10.0,0.5*(j+5)),temp_table,&p,&cooling,&cosmo,&internal_const),x,du;
- float dt = 2.0e-4*units_cgs_conversion_factor(&us,UNIT_CONV_TIME);
- LambdaNet = eagle_cooling_rate_1d_table(u_ini, &dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table, element_cooling_table, element_electron_abundance_table, temp_table, &p, &cooling, &cosmo, &internal_const);
- float u_temp = u_ini + LambdaNet*ratefact*dt;
+ float x_init, u_eq = 2.0e12;
+ for (int j = 5; j < 6; j++) {
+ float u_ini = eagle_convert_temp_to_u_1d_table(pow(10.0, 0.5 * (j + 5)),
+ temp_table, &p, &cooling,
+ &cosmo, &internal_const),
+ x, du;
+ float dt = 2.0e-4 * units_cgs_conversion_factor(&us, UNIT_CONV_TIME);
+ LambdaNet = eagle_cooling_rate_1d_table(
+ u_ini, &dLambdaNet_du, H_plus_He_heat_table,
+ H_plus_He_electron_abundance_table, element_cooling_table,
+ element_electron_abundance_table, temp_table, &p, &cooling, &cosmo,
+ &internal_const);
+ float u_temp = u_ini + LambdaNet * ratefact * dt;
/* RGB removed this **
- if (u_temp > 0) LambdaNext = eagle_cooling_rate_1d_table(u_temp, &dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table, element_cooling_table, element_electron_abundance_table, temp_table, &p, &cooling, &cosmo, &internal_const);
+ if (u_temp > 0) LambdaNext = eagle_cooling_rate_1d_table(u_temp,
+ &dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table,
+ element_cooling_table, element_electron_abundance_table, temp_table, &p,
+ &cooling, &cosmo, &internal_const);
if (fabs(LambdaNet - LambdaNext)/LambdaNet < 0.5) {
u_temp = u_ini;
} else {
- u_temp = eagle_convert_temp_to_u_1d_table(1.0e4,temp_table,&p,&cooling,&cosmo,&internal_const);
- } */
- if (u_temp > u_eq) {
- x_init = log(u_temp);
- } else {
- x_init = log(u_eq);
+ u_temp =
+ eagle_convert_temp_to_u_1d_table(1.0e4,temp_table,&p,&cooling,&cosmo,&internal_const);
+ } */
+ if (u_temp > u_eq) {
+ x_init = log(u_temp);
+ } else {
+ x_init = log(u_eq);
}
- x = newton_iter(x_init,u_ini,H_plus_He_heat_table,H_plus_He_electron_abundance_table,element_cooling_table,element_electron_abundance_table,temp_table,&p,&cosmo,&cooling,&internal_const,dt);
- printf("testing newton integration, u_ini, u %.5e %.5e, temperature initial, final %.5e %.5e\n", u_ini, exp(x), eagle_convert_u_to_temp_1d_table(u_ini,&du,temp_table,&p,&cooling,&cosmo,&internal_const), eagle_convert_u_to_temp_1d_table(exp(x),&du,temp_table,&p,&cooling,&cosmo,&internal_const));
+ x = newton_iter(x_init, u_ini, H_plus_He_heat_table,
+ H_plus_He_electron_abundance_table, element_cooling_table,
+ element_electron_abundance_table, temp_table, &p, &cosmo,
+ &cooling, &internal_const, dt);
+ printf(
+ "testing newton integration, u_ini, u %.5e %.5e, temperature initial, "
+ "final %.5e %.5e\n",
+ u_ini, exp(x),
+ eagle_convert_u_to_temp_1d_table(u_ini, &du, temp_table, &p, &cooling,
+ &cosmo, &internal_const),
+ eagle_convert_u_to_temp_1d_table(exp(x), &du, temp_table, &p, &cooling,
+ &cosmo, &internal_const));
}
free(params);
diff --git a/tests/testCoolingIntegration.c b/tests/testCoolingIntegration.c
index a027cd5a81a61b02e7c317b60899bd0378d69b62..0fa31e295b4fc57930c6957651a82899142b0808 100644
--- a/tests/testCoolingIntegration.c
+++ b/tests/testCoolingIntegration.c
@@ -17,11 +17,11 @@
*
******************************************************************************/
-#include "swift.h"
#include "cooling.h"
+#include "hydro.h"
#include "physical_constants.h"
+#include "swift.h"
#include "units.h"
-#include "hydro.h"
int main() {
@@ -34,18 +34,26 @@ int main() {
struct phys_const internal_const;
struct cooling_function_data cooling;
struct cosmology cosmo;
- //char *parametersFileName = "../examples/CoolingBox/coolingBox.yml";
+ // char *parametersFileName = "../examples/CoolingBox/coolingBox.yml";
char *parametersFileName = "../examples/EAGLE_12/eagle_12.yml";
- enum table_index {EAGLE_Hydrogen=0,EAGLE_Helium,EAGLE_Carbon,EAGLE_Nitrogen,
- EAGLE_Oxygen,EAGLE_Neon,EAGLE_Magnesium,EAGLE_Silicon,
- EAGLE_Iron};
+ enum table_index {
+ EAGLE_Hydrogen = 0,
+ EAGLE_Helium,
+ EAGLE_Carbon,
+ EAGLE_Nitrogen,
+ EAGLE_Oxygen,
+ EAGLE_Neon,
+ EAGLE_Magnesium,
+ EAGLE_Silicon,
+ EAGLE_Iron
+ };
if (params == NULL) error("Error allocating memory for the parameter file.");
message("Reading runtime parameters from file '%s'", parametersFileName);
parser_read_file(parametersFileName, params);
char output_filename[40];
- FILE** output_file = malloc(10*sizeof(FILE*));
+ FILE **output_file = malloc(10 * sizeof(FILE *));
/* And dump the parameters as used. */
parser_write_params_to_file(params, "used_parameters.yml");
@@ -54,123 +62,203 @@ int main() {
phys_const_init(&us, params, &internal_const);
double hydrogen_number_density_cgs = 1.582e-3;
- //double hydrogen_number_density_cgs = 1.0e-1;
- double u,u_cgs, hydrogen_number_density, pressure, gamma, cooling_du_dt, temperature_cgs, newton_func, ratefact;
- u_cgs = 0; cooling_du_dt = 0; temperature_cgs = 0; newton_func = 0;
- //int n_t_i = 2000;
-
+ // double hydrogen_number_density_cgs = 1.0e-1;
+ double u, u_cgs, hydrogen_number_density, pressure, gamma, cooling_du_dt,
+ temperature_cgs, newton_func, ratefact;
+ u_cgs = 0;
+ cooling_du_dt = 0;
+ temperature_cgs = 0;
+ newton_func = 0;
+ // int n_t_i = 2000;
- gamma = 5.0/3.0;
+ gamma = 5.0 / 3.0;
chemistry_init(params, &us, &internal_const, &chemistry_data);
- chemistry_first_init_part(&p,&xp,&chemistry_data);
+ chemistry_first_init_part(&p, &xp, &chemistry_data);
chemistry_print(&chemistry_data);
-
+
cosmology_init(params, &us, &internal_const, &cosmo);
cosmology_print(&cosmo);
-
- float u_ini = 3.357e15;
- u = u_ini/(units_cgs_conversion_factor(&us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS));
- pressure = u*p.rho*(gamma -1.0);
- hydrogen_number_density = hydrogen_number_density_cgs*pow(units_cgs_conversion_factor(&us,UNIT_CONV_LENGTH),3);
- p.rho = hydrogen_number_density*internal_const.const_proton_mass*(1.0+p.chemistry_data.metal_mass_fraction[EAGLE_Helium]/p.chemistry_data.metal_mass_fraction[EAGLE_Hydrogen]);
- p.entropy = pressure/(pow(p.rho,gamma));
+ float u_ini = 3.357e15;
+ u = u_ini / (units_cgs_conversion_factor(&us, UNIT_CONV_ENERGY) /
+ units_cgs_conversion_factor(&us, UNIT_CONV_MASS));
+ pressure = u * p.rho * (gamma - 1.0);
+ hydrogen_number_density =
+ hydrogen_number_density_cgs *
+ pow(units_cgs_conversion_factor(&us, UNIT_CONV_LENGTH), 3);
+ p.rho = hydrogen_number_density * internal_const.const_proton_mass *
+ (1.0 +
+ p.chemistry_data.metal_mass_fraction[EAGLE_Helium] /
+ p.chemistry_data.metal_mass_fraction[EAGLE_Hydrogen]);
+ p.entropy = pressure / (pow(p.rho, gamma));
cosmo.z = 0.0999744;
cooling_init(params, &us, &internal_const, &cooling);
cooling_print(&cooling);
// construct 1d table of cooling rates wrt temperature
- float H_plus_He_heat_table[176]; // WARNING sort out how it is declared/allocated
- float H_plus_He_electron_abundance_table[176]; // WARNING sort out how it is declared/allocated
- float temp_table[176]; // WARNING sort out how it is declared/allocated
- float element_cooling_table[176]; // WARNING sort out how it is declared/allocated
- float element_electron_abundance_table[176]; // WARNING sort out how it is declared/allocated
- for(int k = 0; k < cooling.N_Temp; k++) H_plus_He_heat_table[k] = 0.0;
- for(int k = 0; k < cooling.N_Temp; k++) H_plus_He_electron_abundance_table[k] = 0.0;
- for(int k = 0; k < cooling.N_Temp; k++) temp_table[k] = 0.0;
- for(int k = 0; k < cooling.N_Temp; k++) element_cooling_table[k] = 0.0;
- for(int k = 0; k < cooling.N_Temp; k++) element_electron_abundance_table[k] = 0.0;
- //construct_1d_table_from_4d(&p,&cooling,&cosmo,&internal_const,cooling.table.element_cooling.H_plus_He_heating,H_plus_He_heat_table);
+ float H_plus_He_heat_table[176]; // WARNING sort out how it is
+ // declared/allocated
+ float H_plus_He_electron_abundance_table[176]; // WARNING sort out how it is
+ // declared/allocated
+ float temp_table[176]; // WARNING sort out how it is declared/allocated
+ float element_cooling_table[176]; // WARNING sort out how it is
+ // declared/allocated
+ float element_electron_abundance_table[176]; // WARNING sort out how it is
+ // declared/allocated
+ for (int k = 0; k < cooling.N_Temp; k++) H_plus_He_heat_table[k] = 0.0;
+ for (int k = 0; k < cooling.N_Temp; k++)
+ H_plus_He_electron_abundance_table[k] = 0.0;
+ for (int k = 0; k < cooling.N_Temp; k++) temp_table[k] = 0.0;
+ for (int k = 0; k < cooling.N_Temp; k++) element_cooling_table[k] = 0.0;
+ for (int k = 0; k < cooling.N_Temp; k++)
+ element_electron_abundance_table[k] = 0.0;
+ // construct_1d_table_from_4d(&p,&cooling,&cosmo,&internal_const,cooling.table.element_cooling.H_plus_He_heating,H_plus_He_heat_table);
float XH = p.chemistry_data.metal_mass_fraction[chemistry_element_H];
- //float inn_h = chemistry_get_number_density(&p,&cosmo,chemistry_element_H,&internal_const)*cooling.number_density_scale;
- float inn_h = hydro_get_physical_density(&p,&cosmo)*units_cgs_conversion_factor(&us,UNIT_CONV_DENSITY) * XH / eagle_proton_mass_cgs;
+ // float inn_h =
+ // chemistry_get_number_density(&p,&cosmo,chemistry_element_H,&internal_const)*cooling.number_density_scale;
+ float inn_h = hydro_get_physical_density(&p, &cosmo) *
+ units_cgs_conversion_factor(&us, UNIT_CONV_DENSITY) * XH /
+ eagle_proton_mass_cgs;
ratefact = inn_h * (XH / eagle_proton_mass_cgs);
printf("XH inn_h ratefact %.5e %.5e %.5e\n", XH, inn_h, ratefact);
double dLambdaNet_du, LambdaNet;
float x_init;
- for(int j = 0; j < 1; j++){
- //float u_ini = eagle_convert_temp_to_u_1d_table(pow(10.0,0.5*(j+5)),temp_table,&p,&cooling,&cosmo,&internal_const);
- //float u_ini = 3.357e15;
- u = u_ini/(units_cgs_conversion_factor(&us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS));
- pressure = u*p.rho*(gamma -1.0);
- hydrogen_number_density = hydrogen_number_density_cgs*pow(units_cgs_conversion_factor(&us,UNIT_CONV_LENGTH),3);
- p.rho = hydrogen_number_density*internal_const.const_proton_mass*(1.0+p.chemistry_data.metal_mass_fraction[EAGLE_Helium]/p.chemistry_data.metal_mass_fraction[EAGLE_Hydrogen]);
- p.entropy = pressure/(pow(p.rho,gamma));
-
- float x,du;
+ for (int j = 0; j < 1; j++) {
+ // float u_ini =
+ // eagle_convert_temp_to_u_1d_table(pow(10.0,0.5*(j+5)),temp_table,&p,&cooling,&cosmo,&internal_const);
+ // float u_ini = 3.357e15;
+ u = u_ini / (units_cgs_conversion_factor(&us, UNIT_CONV_ENERGY) /
+ units_cgs_conversion_factor(&us, UNIT_CONV_MASS));
+ pressure = u * p.rho * (gamma - 1.0);
+ hydrogen_number_density =
+ hydrogen_number_density_cgs *
+ pow(units_cgs_conversion_factor(&us, UNIT_CONV_LENGTH), 3);
+ p.rho = hydrogen_number_density * internal_const.const_proton_mass *
+ (1.0 +
+ p.chemistry_data.metal_mass_fraction[EAGLE_Helium] /
+ p.chemistry_data.metal_mass_fraction[EAGLE_Hydrogen]);
+ p.entropy = pressure / (pow(p.rho, gamma));
+
+ float x, du;
double u_temp;
- //float dt = 2.0e-4*units_cgs_conversion_factor(&us,UNIT_CONV_TIME);
- float dt = 1.73798e-3*units_cgs_conversion_factor(&us,UNIT_CONV_TIME);
- construct_1d_table_from_4d_H_He(&p,&cooling,&cosmo,&internal_const,cooling.table.element_cooling.H_plus_He_heating,H_plus_He_heat_table, &u_temp, u_ini, dt);
- construct_1d_table_from_4d(&p,&cooling,&cosmo,&internal_const,cooling.table.element_cooling.H_plus_He_electron_abundance,H_plus_He_electron_abundance_table);
- construct_1d_table_from_4d(&p,&cooling,&cosmo,&internal_const,cooling.table.element_cooling.temperature,temp_table);
- construct_1d_table_from_4d_elements(&p,&cooling,&cosmo,&internal_const,cooling.table.element_cooling.metal_heating,element_cooling_table);
- construct_1d_table_from_3d(&p,&cooling,&cosmo,&internal_const,cooling.table.element_cooling.electron_abundance,element_electron_abundance_table);
-
- sprintf(output_filename, "%s%d%s", "cooling_integration_output_",j,".dat");
+ // float dt = 2.0e-4*units_cgs_conversion_factor(&us,UNIT_CONV_TIME);
+ float dt = 1.73798e-3 * units_cgs_conversion_factor(&us, UNIT_CONV_TIME);
+ construct_1d_table_from_4d_H_He(
+ &p, &cooling, &cosmo, &internal_const,
+ cooling.table.element_cooling.H_plus_He_heating, H_plus_He_heat_table,
+ &u_temp, u_ini, dt);
+ construct_1d_table_from_4d(
+ &p, &cooling, &cosmo, &internal_const,
+ cooling.table.element_cooling.H_plus_He_electron_abundance,
+ H_plus_He_electron_abundance_table);
+ construct_1d_table_from_4d(&p, &cooling, &cosmo, &internal_const,
+ cooling.table.element_cooling.temperature,
+ temp_table);
+ construct_1d_table_from_4d_elements(
+ &p, &cooling, &cosmo, &internal_const,
+ cooling.table.element_cooling.metal_heating, element_cooling_table);
+ construct_1d_table_from_3d(&p, &cooling, &cosmo, &internal_const,
+ cooling.table.element_cooling.electron_abundance,
+ element_electron_abundance_table);
+
+ sprintf(output_filename, "%s%d%s", "cooling_integration_output_", j,
+ ".dat");
output_file[j] = fopen(output_filename, "w");
- if (output_file[j] == NULL)
- {
- printf("Error opening file!\n");
- exit(1);
+ if (output_file[j] == NULL) {
+ printf("Error opening file!\n");
+ exit(1);
}
- for(int k = 0; k < cooling.N_Temp-1; k++){
- float lambda1 = H_plus_He_heat_table[k] + element_cooling_table[k]*H_plus_He_electron_abundance_table[k]/element_electron_abundance_table[k];
- float lambda2 = H_plus_He_heat_table[k+1] + element_cooling_table[k+1]*H_plus_He_electron_abundance_table[k+1]/element_electron_abundance_table[k+1];
- //printf("temperature %.5e, internal energy %.5e\n",pow(10.0,temp_table[k]), eagle_convert_temp_to_u_1d_table(pow(10.0,temp_table[k]),temp_table,&p,&cooling,&cosmo,&internal_const));
- float u2 = eagle_convert_temp_to_u_1d_table(pow(10.0,temp_table[k+1]),temp_table,&p,&cooling,&cosmo,&internal_const);
- float u1 = eagle_convert_temp_to_u_1d_table(pow(10.0,temp_table[k]),temp_table,&p,&cooling,&cosmo,&internal_const);
+ for (int k = 0; k < cooling.N_Temp - 1; k++) {
+ float lambda1 = H_plus_He_heat_table[k] +
+ element_cooling_table[k] *
+ H_plus_He_electron_abundance_table[k] /
+ element_electron_abundance_table[k];
+ float lambda2 = H_plus_He_heat_table[k + 1] +
+ element_cooling_table[k + 1] *
+ H_plus_He_electron_abundance_table[k + 1] /
+ element_electron_abundance_table[k + 1];
+ // printf("temperature %.5e, internal energy
+ // %.5e\n",pow(10.0,temp_table[k]),
+ // eagle_convert_temp_to_u_1d_table(pow(10.0,temp_table[k]),temp_table,&p,&cooling,&cosmo,&internal_const));
+ float u2 = eagle_convert_temp_to_u_1d_table(pow(10.0, temp_table[k + 1]),
+ temp_table, &p, &cooling,
+ &cosmo, &internal_const);
+ float u1 = eagle_convert_temp_to_u_1d_table(pow(10.0, temp_table[k]),
+ temp_table, &p, &cooling,
+ &cosmo, &internal_const);
float delta_u = u2 - u1;
- //fprintf(output_file[j], "%.5e %.5e %.5e\n", pow(10.0,temp_table[k]), u1 - u_ini - lambda1*ratefact*dt, (lambda2 - lambda1)/delta_u);
- fprintf(output_file[j], "%.5e %.5e %.5e\n", u1, u1 - u_ini - lambda1*ratefact*dt, (lambda2 - lambda1)/delta_u);
+ // fprintf(output_file[j], "%.5e %.5e %.5e\n", pow(10.0,temp_table[k]), u1
+ // - u_ini - lambda1*ratefact*dt, (lambda2 - lambda1)/delta_u);
+ fprintf(output_file[j], "%.5e %.5e %.5e\n", u1,
+ u1 - u_ini - lambda1 * ratefact * dt,
+ (lambda2 - lambda1) / delta_u);
}
fclose(output_file[j]);
- LambdaNet = eagle_cooling_rate_1d_table(u_ini, &dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table, element_cooling_table, element_electron_abundance_table, temp_table, &p, &cooling, &cosmo, &internal_const);
+ LambdaNet = eagle_cooling_rate_1d_table(
+ u_ini, &dLambdaNet_du, H_plus_He_heat_table,
+ H_plus_He_electron_abundance_table, element_cooling_table,
+ element_electron_abundance_table, temp_table, &p, &cooling, &cosmo,
+ &internal_const);
double u_eq = 5.0e12;
- double u_temp_guess = u_ini + LambdaNet*ratefact*dt;
- printf("u_guess, u_temp_guess, u_ini, LambdaNet, dLambdaNet_du %.5e %.5e %.5e %.5e %.5e %.5e %.5e \n", u_temp, u_temp_guess, u_ini, LambdaNet, dLambdaNet_du, ratefact, dt);
+ double u_temp_guess = u_ini + LambdaNet * ratefact * dt;
+ printf(
+ "u_guess, u_temp_guess, u_ini, LambdaNet, dLambdaNet_du %.5e %.5e %.5e "
+ "%.5e %.5e %.5e %.5e \n",
+ u_temp, u_temp_guess, u_ini, LambdaNet, dLambdaNet_du, ratefact, dt);
- //if ((LambdaNet < 0 && u_temp < u_temp_guess) ||
+ // if ((LambdaNet < 0 && u_temp < u_temp_guess) ||
// (LambdaNet >= 0 && u_temp > u_temp_guess))
// u_temp = u_temp_guess;
u_temp = u_temp_guess;
- if ((LambdaNet < 0 && u_temp < u_eq) ||
- (LambdaNet >= 0 && u_temp > u_eq))
+ if ((LambdaNet < 0 && u_temp < u_eq) || (LambdaNet >= 0 && u_temp > u_eq))
u_temp = u_eq;
x_init = log(u_temp);
- //int printflag = 1;
- //x = newton_iter(x_init,u_ini,H_plus_He_heat_table,H_plus_He_electron_abundance_table,element_cooling_table,element_electron_abundance_table,temp_table,&p,&cosmo,&cooling,&internal_const,dt,&printflag);
- x = bisection_iter(x_init,u_ini,H_plus_He_heat_table,H_plus_He_electron_abundance_table,element_cooling_table,element_electron_abundance_table,temp_table,&p,&cosmo,&cooling,&internal_const,dt);
- //printf("testing newton integration, u_ini, u %.5e %.5e, temperature initial, final %.5e %.5e\n", u_ini, exp(x), eagle_convert_u_to_temp_1d_table(u_ini,&du,temp_table,&p,&cooling,&cosmo,&internal_const), eagle_convert_u_to_temp_1d_table(exp(x),&du,temp_table,&p,&cooling,&cosmo,&internal_const));
+ // int printflag = 1;
+ // x =
+ // newton_iter(x_init,u_ini,H_plus_He_heat_table,H_plus_He_electron_abundance_table,element_cooling_table,element_electron_abundance_table,temp_table,&p,&cosmo,&cooling,&internal_const,dt,&printflag);
+ x = bisection_iter(x_init, u_ini, H_plus_He_heat_table,
+ H_plus_He_electron_abundance_table,
+ element_cooling_table, element_electron_abundance_table,
+ temp_table, &p, &cosmo, &cooling, &internal_const, dt);
+ // printf("testing newton integration, u_ini, u %.5e %.5e, temperature
+ // initial, final %.5e %.5e\n", u_ini, exp(x),
+ // eagle_convert_u_to_temp_1d_table(u_ini,&du,temp_table,&p,&cooling,&cosmo,&internal_const),
+ // eagle_convert_u_to_temp_1d_table(exp(x),&du,temp_table,&p,&cooling,&cosmo,&internal_const));
u = u_ini;
double u_next;
int nt = 10000;
- float dt_sub = dt/nt;
- for(int t = 0; t < nt; t++){
- LambdaNet = eagle_cooling_rate_1d_table(u, &dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table, element_cooling_table, element_electron_abundance_table, temp_table, &p, &cooling, &cosmo, &internal_const);
- u_next = u + LambdaNet*ratefact*dt_sub;
- printf("here u_next u lambda_net ratefact dt_sub, du t %.5e %.5e %.5e %.5e %.5e %.5e %d\n",u_next, u, LambdaNet, ratefact, dt_sub, LambdaNet*ratefact*dt_sub,t);
+ float dt_sub = dt / nt;
+ for (int t = 0; t < nt; t++) {
+ LambdaNet = eagle_cooling_rate_1d_table(
+ u, &dLambdaNet_du, H_plus_He_heat_table,
+ H_plus_He_electron_abundance_table, element_cooling_table,
+ element_electron_abundance_table, temp_table, &p, &cooling, &cosmo,
+ &internal_const);
+ u_next = u + LambdaNet * ratefact * dt_sub;
+ printf(
+ "here u_next u lambda_net ratefact dt_sub, du t %.5e %.5e %.5e %.5e "
+ "%.5e %.5e %d\n",
+ u_next, u, LambdaNet, ratefact, dt_sub, LambdaNet * ratefact * dt_sub,
+ t);
u = u_next;
}
- printf("testing newton integration, u_ini, u, u subcycle %.5e %.5e %.5e, temperature initial, final, subcycled %.5e %.5e %.5e\n", u_ini, exp(x), u, eagle_convert_u_to_temp_1d_table(u_ini,&du,temp_table,&p,&cooling,&cosmo,&internal_const), eagle_convert_u_to_temp_1d_table(exp(x),&du,temp_table,&p,&cooling,&cosmo,&internal_const), eagle_convert_u_to_temp_1d_table(u,&du,temp_table,&p,&cooling,&cosmo,&internal_const));
-
+ printf(
+ "testing newton integration, u_ini, u, u subcycle %.5e %.5e %.5e, "
+ "temperature initial, final, subcycled %.5e %.5e %.5e\n",
+ u_ini, exp(x), u,
+ eagle_convert_u_to_temp_1d_table(u_ini, &du, temp_table, &p, &cooling,
+ &cosmo, &internal_const),
+ eagle_convert_u_to_temp_1d_table(exp(x), &du, temp_table, &p, &cooling,
+ &cosmo, &internal_const),
+ eagle_convert_u_to_temp_1d_table(u, &du, temp_table, &p, &cooling,
+ &cosmo, &internal_const));
}
free(params);