Skip to content
Snippets Groups Projects
Commit 77abccd6 authored by lhausamm's avatar lhausamm
Browse files

Deal correctly with input

parent 86869b5a
No related branches found
No related tags found
1 merge request!469Update cooling grackle
Hide whitespace changes
Inline Side-by-side
src/cooling/grackle/cooling.h
+ 67
− 87
View file @ 77abccd6
......@@ -39,35 +39,26 @@
/* include the grackle wrapper */
#include "grackle_wrapper.h"
/*! This function computes the new entropy due to the cooling,
* between step t0 and t1.
*/
static INLINE double do_cooling_grackle(double energy, double density, double dtime, double *ne, double Z, double a_now)
{
static INLINE double do_cooling_grackle(double energy, double density,
double dtime, double* ne, double Z,
double a_now) {
/*********************************************************************
call to the main chemistry solver
*********************************************************************/
if (wrap_do_cooling(density, &energy, dtime,Z, a_now) == 0) {
if (wrap_do_cooling(density, &energy, dtime, Z, a_now) == 0) {
error("Error in do_cooling.\n");
return 0;
}
return energy;
return energy;
}
/**
* @brief Apply the cooling function to a particle.
*
......@@ -84,49 +75,41 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
const struct unit_system* restrict us,
const struct cooling_function_data* restrict cooling,
struct part* restrict p, struct xpart* restrict xp, float dt) {
/* Get current internal energy (dt=0) */
const float u_old = hydro_get_internal_energy(p);
/* Get current density */
const float rho = hydro_get_density(p);
/* Actual scaling fractor */
const float a_now = 1. / (1. + cooling->GrackleRedshift); ; /* must be chaged !!! */
if (cooling->GrackleRedshift == -1) error("TODO time dependant redshift");
const float a_now = 1. / (1. + cooling->GrackleRedshift);
; /* must be chaged !!! */
double ne, Z;
double ne,Z;
Z = 0.02041; /* 0.02041 (= 1 Zsun in Grackle v2.0, but = 1.5761 Zsun in
Grackle v2.1) */
ne = 0.0;
/* mass fraction of eletron */ /* useless for GRACKLE_CHEMISTRY = 0 */
Z = 0.02041; /* 0.02041 (= 1 Zsun in Grackle v2.0, but = 1.5761 Zsun in Grackle v2.1) */
ne= 0.0; /* mass fraction of eletron */ /* useless for GRACKLE_CHEMISTRY = 0 */
float u_new;
float delta_u;
u_new = do_cooling_grackle(u_old, rho, dt, &ne, Z, a_now);
//u_new = u_old * 0.99;
u_new = do_cooling_grackle(u_old, rho, dt, &ne, Z, a_now);
// u_new = u_old * 0.99;
//if (u_new < 0)
//if (p->id==50356)
// printf("WARNING !!! ID=%llu u_old=%g u_new=%g rho=%g dt=%g ne=%g Z=%g a_now=%g\n",p->id,u_old,u_new,rho,dt,ne,Z,a_now);
// if (u_new < 0)
// if (p->id==50356)
// printf("WARNING !!! ID=%llu u_old=%g u_new=%g rho=%g dt=%g ne=%g Z=%g
// a_now=%g\n",p->id,u_old,u_new,rho,dt,ne,Z,a_now);
delta_u = u_new - u_old;
/* record energy lost */
xp->cooling_data.radiated_energy += -delta_u * hydro_get_mass(p);
xp->cooling_data.radiated_energy += -delta_u * hydro_get_mass(p);
/* Update the internal energy */
hydro_set_internal_energy_dt(p, delta_u / dt);
hydro_set_internal_energy_dt(p, delta_u / dt);
}
/**
......@@ -157,7 +140,7 @@ __attribute__((always_inline)) INLINE static float cooling_timestep(
__attribute__((always_inline)) INLINE static void cooling_init_part(
const struct part* restrict p, struct xpart* restrict xp) {
xp->cooling_data.radiated_energy = 0.f;
xp->cooling_data.radiated_energy = 0.f;
}
/**
......@@ -180,59 +163,56 @@ __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 swift_params* parameter_file, const struct unit_system* us,
const struct phys_const* phys_const,
struct cooling_function_data* cooling) {
char cloudytable[200];
double units_density, units_length, units_time;
int grackle_chemistry;
int UVbackground;
parser_get_param_string(parameter_file, "GrackleCooling:GrackleCloudyTable",
cooling->GrackleCloudyTable);
cooling->UVbackground =
parser_get_param_int(parameter_file, "GrackleCooling:UVbackground");
cooling->GrackleRedshift =
parser_get_param_double(parameter_file, "GrackleCooling:GrackleRedshift");
cooling->GrackleHSShieldingDensityThreshold = parser_get_param_double(
parameter_file, "GrackleCooling:GrackleHSShieldingDensityThreshold");
UVbackground = cooling->UVbackground;
grackle_chemistry = 0; /* forced to be zero : read table */
units_density = us->UnitMass_in_cgs / pow(us->UnitLength_in_cgs, 3);
units_length = us->UnitLength_in_cgs;
units_time = us->UnitTime_in_cgs;
#ifdef SWIFT_DEBUG_CHECKS
float threshold = cooling->GrackleHSShieldingDensityThreshold;
threshold /= phys_const->const_proton_mass;
threshold /= pow(us->UnitLength_in_cgs, 3);
message("***************************************");
message("initializing grackle cooling function");
message("");
message("CloudyTable = %s",
cooling->GrackleCloudyTable);
message("UVbackground = %d", UVbackground);
message("GrackleRedshift = %g", cooling->GrackleRedshift);
message("GrackleHSShieldingDensityThreshold = %g atom/cc", threshold);
#endif
if (wrap_init_cooling(cooling->GrackleCloudyTable, UVbackground,
units_density, units_length, units_time,
grackle_chemistry) != 1) {
error("Error in initialize_chemistry_data.");
}
parser_get_param_string(parameter_file,"GrackleCooling:GrackleCloudyTable",cooling->GrackleCloudyTable);
cooling->UVbackground = parser_get_param_int(parameter_file, "GrackleCooling:UVbackground");
cooling->GrackleRedshift = parser_get_param_double(parameter_file, "GrackleCooling:GrackleRedshift");
cooling->GrackleHSShieldingDensityThreshold = parser_get_param_double(parameter_file, "GrackleCooling:GrackleHSShieldingDensityThreshold");
// FIXME : Why a strcpy ?
strcpy(cloudytable,cooling->GrackleCloudyTable);
UVbackground = cooling->UVbackground;
grackle_chemistry = 0; /* forced to be zero : read table */
units_density = us->UnitMass_in_cgs/pow(us->UnitLength_in_cgs,3);
units_length = us->UnitLength_in_cgs;
units_time = us->UnitTime_in_cgs;
printf(" ***************************************\n");
printf(" initializing grackle cooling function\n");
printf(" \n");
printf(" CloudyTable = %s\n",cloudytable);
printf(" UVbackground = %d\n",UVbackground);
printf(" GrackleRedshift = %g\n",cooling->GrackleRedshift);
printf(" GrackleHSShieldingDensityThreshold = %g atom/cc\n",cooling->GrackleHSShieldingDensityThreshold);
if(wrap_init_cooling(cloudytable,UVbackground,units_density, units_length, units_time, grackle_chemistry) != 1)
{
fprintf(stderr,"Error in initialize_chemistry_data.");
exit(-1);
}
printf(" \n");
printf(" ***************************************\n");
#ifdef SWIFT_DEBUG_CHECKS
message("");
message("***************************************");
#endif
}
/**
......
src/cooling/grackle/cooling_struct.h
+ 0
− 3
View file @ 77abccd6
......@@ -33,7 +33,6 @@ struct cooling_function_data {
int UVbackground;
double GrackleRedshift;
double GrackleHSShieldingDensityThreshold;
};
/**
......@@ -43,8 +42,6 @@ struct cooling_xpart_data {
/*! Energy radiated away by this particle since the start of the run */
float radiated_energy;
};
#endif /* SWIFT_COOLING_STRUCT_NONE_H */
src/cooling/grackle/grackle_wrapper.c
+ 76
− 115
View file @ 77abccd6
......@@ -7,26 +7,23 @@
/
/ Distributed under the terms of the Enzo Public Licence.
/
/ The full license is in the file LICENSE, distributed with this
/ The full license is in the file LICENSE, distributed with this
/ software.
************************************************************************/
#include "grackle_wrapper.h"
#define GRACKLE_DEBUG
#ifdef GRACKLE_DEBUG
#ifdef SWIFT_DEBUG_CHECKS
#include <assert.h>
#define GRACKLE_ASSERT(X) assert( (X) )
#define GRACKLE_ASSERT(X) assert((X))
#else
#define GRACKLE_ASSERT(X)
#define GRACKLE_ASSERT(X)
#endif
code_units my_units;
// arrays passed to grackle as input and to be filled
#define FIELD_SIZE 1
#define FIELD_SIZE 1
gr_float HDI_density[FIELD_SIZE];
......@@ -34,26 +31,28 @@ gr_float HDI_density[FIELD_SIZE];
// grid_start and grid_end are used to ignore ghost zones.
const int grid_rank = 1;
int wrap_init_cooling(char* CloudyTable, int UVbackground, double udensity, double ulength, double utime, int grackle_chemistry){
int wrap_init_cooling(char *CloudyTable, int UVbackground, double udensity,
double ulength, double utime, int grackle_chemistry) {
#ifdef GRACKLE_DEBUG
#ifdef SWIFT_DEBUG_CHECKS
grackle_verbose = 1;
#endif
message("cooling");
double velocity_units;
// First, set up the units system.
// These are conversions from code units to cgs.
my_units.a_units = 1.0; // units for the expansion factor (1/1+zi)
my_units.a_units = 1.0; // units for the expansion factor (1/1+zi)
my_units.comoving_coordinates = 0; /* so, according to the doc, we assume here all physical quantities to be in proper coordiname (not comobile) */
my_units.comoving_coordinates =
0; /* so, according to the doc, we assume here all physical quantities to
be in proper coordinate (not comobile) */
my_units.density_units = udensity;
my_units.length_units = ulength;
my_units.time_units = utime;
velocity_units = my_units.a_units * my_units.length_units /my_units.time_units;
my_units.velocity_units = velocity_units;
my_units.time_units = utime;
velocity_units =
my_units.a_units * my_units.length_units / my_units.time_units;
my_units.velocity_units = velocity_units;
// Second, create a chemistry object for parameters and rate data.
if (set_default_chemistry_parameters() == 0) {
......@@ -61,77 +60,65 @@ int wrap_init_cooling(char* CloudyTable, int UVbackground, double udensity, doub
return 0;
}
// Set parameter values for chemistry.
grackle_data.use_grackle = 1;
grackle_data.use_grackle = 1;
grackle_data.with_radiative_cooling = 1;
grackle_data.primordial_chemistry = grackle_chemistry; // molecular network with H, He, D
grackle_data.metal_cooling = 1; // metal cooling on
grackle_data.UVbackground = UVbackground;
grackle_data.primordial_chemistry =
grackle_chemistry; // molecular network with H, He, D
grackle_data.metal_cooling = 1; // metal cooling on
grackle_data.UVbackground = UVbackground;
grackle_data.grackle_data_file = CloudyTable;
// Finally, initialize the chemistry object.
// snl: a_value is not the true initial ae!! This should get set during update_UVbackground
// snl: a_value is not the true initial ae!! This should get set during
// update_UVbackground
double initial_redshift = 0.;
double a_value = 1. / (1. + initial_redshift);
// Finally, initialize the chemistry object.
if (initialize_chemistry_data(&my_units, a_value) == 0) {
fprintf(stderr, "Error in initialize_chemistry_data.\n");
return 0;
}
return 1;
}
//int wrap_update_UVbackground_rates(double auni) {
// // The UV background rates must be updated before
// // calling the other functions.
// /* a_value = auni / my_units.a_units; */
// double a_value = auni / my_units.a_units;
// if (update_UVbackground_rates(my_units, a_value) == 0) {
// return 0;
// }
// return 1;
//}
int wrap_set_UVbackground_On() {
// The UV background rates is enabled
grackle_data.UVbackground = 1;
return 1;
int wrap_set_UVbackground_on() {
// The UV background rates is enabled
grackle_data.UVbackground = 1;
return 1;
}
int wrap_set_UVbackground_Off() {
// The UV background rates is disabled
grackle_data.UVbackground = 0;
return 1;
int wrap_set_UVbackground_off() {
// The UV background rates is disabled
grackle_data.UVbackground = 0;
return 1;
}
int wrap_get_cooling_time(double rho, double u,double Z, double a_now, double *coolingtime)
{
int wrap_get_cooling_time(double rho, double u, double Z, double a_now,
double *coolingtime) {
gr_float den_factor = 1.0;
gr_float u_factor = 1.0;
gr_float x_velocity[FIELD_SIZE] = { 0.0 };
gr_float y_velocity[FIELD_SIZE] = { 0.0 };
gr_float z_velocity[FIELD_SIZE] = { 0.0 };
gr_float Density[FIELD_SIZE] = { rho*den_factor };
gr_float metal_density[FIELD_SIZE] = { Z * Density[0] };
gr_float energy[FIELD_SIZE] = { u * u_factor };
gr_float cooling_time[FIELD_SIZE] = { 0.0 };
gr_float x_velocity[FIELD_SIZE] = {0.0};
gr_float y_velocity[FIELD_SIZE] = {0.0};
gr_float z_velocity[FIELD_SIZE] = {0.0};
gr_float density[FIELD_SIZE] = {rho * den_factor};
gr_float metal_density[FIELD_SIZE] = {Z * density[0]};
gr_float energy[FIELD_SIZE] = {u * u_factor};
gr_float cooling_time[FIELD_SIZE] = {0.0};
int grid_dimension[3] = {1, 0, 0};
int grid_start[3] = {0, 0, 0};
int grid_end[3] = {0, 0, 0};
if(FIELD_SIZE != 1){
fprintf(stderr,"field_size must currently be set to 1.\n");
return 0;
}
if (FIELD_SIZE != 1) {
fprintf(stderr, "field_size must currently be set to 1.\n");
return 0;
}
// passed density and energy are proper
/*
......@@ -143,19 +130,15 @@ int wrap_get_cooling_time(double rho, double u,double Z, double a_now, double *c
u_factor = 1.0;
}
*/
if (calculate_cooling_time_table(&my_units,
a_now,
grid_rank, grid_dimension,
grid_start, grid_end,
Density, energy,
x_velocity, y_velocity, z_velocity,
metal_density,
cooling_time) == 0) {
fprintf(stderr,"Error in calculate_cooling_time.\n");
if (calculate_cooling_time_table(&my_units, a_now, grid_rank, grid_dimension,
grid_start, grid_end, density, energy,
x_velocity, y_velocity, z_velocity,
metal_density, cooling_time) == 0) {
fprintf(stderr, "Error in calculate_cooling_time.\n");
return 0;
}
// return updated chemistry and energy
for (int i = 0; i < FIELD_SIZE; i++) {
*coolingtime = cooling_time[i];
......@@ -164,62 +147,40 @@ int wrap_get_cooling_time(double rho, double u,double Z, double a_now, double *c
return 1;
}
int wrap_do_cooling(double rho, double *u, double dt,double Z, double a_now)
{
int wrap_do_cooling(double rho, double *u, double dt, double Z, double a_now) {
gr_float den_factor = 1.0;
gr_float u_factor = 1.0;
gr_float x_velocity[FIELD_SIZE] = { 0.0 };
gr_float y_velocity[FIELD_SIZE] = { 0.0 };
gr_float z_velocity[FIELD_SIZE] = { 0.0 };
gr_float Density[FIELD_SIZE] = { rho*den_factor };
gr_float metal_density[FIELD_SIZE] = { Z * Density[0] };
gr_float energy[FIELD_SIZE] = { (*u) * u_factor };
gr_float x_velocity[FIELD_SIZE] = {0.0};
gr_float y_velocity[FIELD_SIZE] = {0.0};
gr_float z_velocity[FIELD_SIZE] = {0.0};
gr_float density[FIELD_SIZE] = {rho * den_factor};
gr_float metal_density[FIELD_SIZE] = {Z * density[0]};
gr_float energy[FIELD_SIZE] = {(*u) * u_factor};
int grid_dimension[3] = {1, 0, 0};
int grid_start[3] = {0, 0, 0};
int grid_end[3] = {0, 0, 0};
GRACKLE_ASSERT(FIELD_SIZE == 1);
// passed density and energy are proper
/*
if(my_units.comoving_coordinates){
den_factor = pow(a_now, 3);
u_factor = pow(a_now, 0);
} else {
den_factor = 1.0;
u_factor = 1.0;
}
*/
#ifdef GRACKLE_DEBUG
#ifdef SWIFT_DEBUG_CHECKS
double old_value = energy[0];
#endif
if (solve_chemistry_table(&my_units,
a_now, dt,
grid_rank, grid_dimension,
grid_start, grid_end,
Density, energy,
x_velocity, y_velocity, z_velocity,
metal_density) == 0) {
fprintf(stderr,"Error in solve_chemistry.\n");
if (solve_chemistry_table(&my_units, a_now, dt, grid_rank, grid_dimension,
grid_start, grid_end, density, energy, x_velocity,
y_velocity, z_velocity, metal_density) == 0) {
error("Error in solve_chemistry.");
return 0;
}
}
GRACKLE_ASSERT(old_value != energy[0]);
// return updated chemistry and energy
for (int i = 0; i < FIELD_SIZE; i++) {
*u = energy[i]/u_factor;
*u = energy[i] / u_factor;
}
return 1;
}
src/cooling/grackle/grackle_wrapper.h
+ 11
− 17
View file @ 77abccd6
......@@ -7,37 +7,31 @@
/
/ Distributed under the terms of the Enzo Public Licence.
/
/ The full license is in the file LICENSE, distributed with this
/ The full license is in the file LICENSE, distributed with this
/ software.
************************************************************************/
#ifndef SWIFT_COOLING_GRACKLE_WRAPPER_H
#define SWIFT_COOLING_GRACKLE_WRAPPER_H
#include "config.h"
#include "error.h"
#include <grackle.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "config.h"
#include "error.h"
int wrap_init_cooling(char *CloudyTable, int UVbackground, double udensity,
double ulength, double utime, int grackle_chemistry);
int wrap_set_UVbackground_on();
int wrap_init_cooling(char* CloudyTable,int UVbackground, double udensity, double ulength, double utime,
int grackle_chemistry);
//int wrap_update_UVbackground_rates(double auni);
int wrap_set_UVbackground_On();
int wrap_set_UVbackground_Off();
int wrap_get_cooling_time(double rho, double u, double Z, double a_now, double *coolingtime);
int wrap_do_cooling(double density, double *energy, double dtime, double Z, double a_now);
int wrap_set_UVbackground_off();
int wrap_get_cooling_time(double rho, double u, double Z, double a_now,
double *coolingtime);
int wrap_do_cooling(double density, double *energy, double dtime, double Z,
double a_now);
#endif /* SWIFT_COOLING_GRACKLE_WRAPPER_H */
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment