diff --git a/.clang-format-ignore b/.clang-format-ignore index 36fa238b6..a75d728e6 100644 --- a/.clang-format-ignore +++ b/.clang-format-ignore @@ -11,11 +11,11 @@ src/clib/c_wrappers/wrap_interpolators_g.c src/clib/calc_tdust_3d.cpp -src/clib/calculate_cooling_time.c -src/clib/calculate_dust_temperature.c -src/clib/calculate_gamma.c -src/clib/calculate_pressure.c -src/clib/calculate_temperature.c +src/clib/calculate_cooling_time.cpp +src/clib/calculate_dust_temperature.cpp +src/clib/calculate_gamma.cpp +src/clib/calculate_pressure.cpp +src/clib/calculate_temperature.cpp src/clib/chemistry_solver_funcs.hpp src/clib/cie_thin_cooling_rate_tables.h src/clib/cool1d_multi_g-cpp.C @@ -42,7 +42,7 @@ src/clib/interpolate.hpp src/clib/phys_constants.h src/clib/rate_functions.c src/clib/set_default_chemistry_parameters.c -src/clib/solve_chemistry.c +src/clib/solve_chemistry.cpp src/clib/solve_rate_cool_g-cpp.cpp src/clib/status_reporting.c src/clib/status_reporting.h diff --git a/src/clib/CMakeLists.txt b/src/clib/CMakeLists.txt index 6e3b757c8..54ebd6778 100644 --- a/src/clib/CMakeLists.txt +++ b/src/clib/CMakeLists.txt @@ -82,11 +82,6 @@ configure_file(auto_general.c.in auto_general.c @ONLY) add_library(Grackle_Grackle # C source files - calculate_cooling_time.c - calculate_dust_temperature.c - calculate_gamma.c - calculate_pressure.c - calculate_temperature.c dynamic_api.c grackle_units.c index_helper.c @@ -94,7 +89,6 @@ add_library(Grackle_Grackle initialize_UVbackground_data.c initialize_UVbackground_data.h rate_functions.c set_default_chemistry_parameters.c - solve_chemistry.c status_reporting.c status_reporting.h update_UVbackground_rates.c utils.c @@ -107,6 +101,11 @@ add_library(Grackle_Grackle # C++ Source (and Private Header Files) calc_tdust_3d.cpp calc_tdust_3d.h calc_temp_cloudy_g.cpp calc_temp_cloudy_g.h + calculate_cooling_time.cpp + calculate_dust_temperature.cpp + calculate_gamma.cpp + calculate_pressure.cpp + calculate_temperature.cpp ceiling_species.hpp collisional_rate_props.cpp collisional_rate_props.hpp cool1d_multi_g.cpp cool1d_multi_g.hpp @@ -121,6 +120,7 @@ add_library(Grackle_Grackle make_consistent.cpp make_consistent.hpp opaque_storage.hpp rate_utils.cpp + solve_chemistry.cpp scale_fields.cpp scale_fields.hpp solve_rate_cool_g-cpp.cpp solve_rate_cool_g-cpp.h step_rate_newton_raphson.hpp diff --git a/src/clib/calculate_cooling_time.c b/src/clib/calculate_cooling_time.cpp similarity index 61% rename from src/clib/calculate_cooling_time.c rename to src/clib/calculate_cooling_time.cpp index 819858bd5..c9a2afb5e 100644 --- a/src/clib/calculate_cooling_time.c +++ b/src/clib/calculate_cooling_time.cpp @@ -11,38 +11,30 @@ / software. ************************************************************************/ -#include -#include -#include +#include #include "grackle.h" -#include "grackle_macros.h" #include "cool_multi_time_g.h" -#include "grackle_types.h" #include "internal_units.h" -#include "phys_constants.h" #include "utils.h" -extern chemistry_data *grackle_data; -extern chemistry_data_storage grackle_rates; - /* function prototypes */ -int update_UVbackground_rates(chemistry_data *my_chemistry, - chemistry_data_storage *my_rates, - photo_rate_storage *my_uvb_rates, - code_units *my_units); +extern "C" int update_UVbackground_rates(chemistry_data *my_chemistry, + chemistry_data_storage *my_rates, + photo_rate_storage *my_uvb_rates, + code_units *my_units); -int local_calculate_cooling_time(chemistry_data *my_chemistry, - chemistry_data_storage *my_rates, - code_units *my_units, - grackle_field_data *my_fields, - gr_float *cooling_time) +extern "C" int local_calculate_cooling_time(chemistry_data *my_chemistry, + chemistry_data_storage *my_rates, + code_units *my_units, + grackle_field_data *my_fields, + gr_float *cooling_time) { /* Return if this doesn't concern us. */ if (!my_chemistry->use_grackle) - return SUCCESS; + return GR_SUCCESS; /* Update UV background rates. */ photo_rate_storage my_uvb_rates; @@ -56,9 +48,9 @@ int local_calculate_cooling_time(chemistry_data *my_chemistry, if (my_chemistry->UVbackground == 1) { if (update_UVbackground_rates(my_chemistry, my_rates, - &my_uvb_rates, my_units) == FAIL) { - fprintf(stderr, "Error in update_UVbackground_rates.\n"); - return FAIL; + &my_uvb_rates, my_units) != GR_SUCCESS) { + std::fprintf(stderr, "Error in update_UVbackground_rates.\n"); + return GR_FAIL; } } else { @@ -81,17 +73,14 @@ int local_calculate_cooling_time(chemistry_data *my_chemistry, } /* Check for a metal field. */ - - int metal_field_present = TRUE; - if (my_fields->metal_density == NULL) - metal_field_present = FALSE; + int metal_field_present = (my_fields->metal_density != NULL) ? TRUE : FALSE; InternalGrUnits internalu = new_internalu_(my_units); /* Error checking for H2 shielding approximation */ if (self_shielding_err_check(my_chemistry, my_fields, - "local_calculate_temperature") == FAIL) { - return FAIL; + "local_calculate_temperature") != GR_SUCCESS) { + return GR_FAIL; } /* Solve cooling equations. */ @@ -100,17 +89,17 @@ int local_calculate_cooling_time(chemistry_data *my_chemistry, my_fields, my_uvb_rates ); - return SUCCESS; + return GR_SUCCESS; } -int calculate_cooling_time(code_units *my_units, - grackle_field_data *my_fields, - gr_float *cooling_time) +extern "C" int calculate_cooling_time(code_units *my_units, + grackle_field_data *my_fields, + gr_float *cooling_time) { if (local_calculate_cooling_time(grackle_data, &grackle_rates, my_units, - my_fields, cooling_time) == FAIL) { - fprintf(stderr, "Error in local_calculate_cooling_time.\n"); - return FAIL; + my_fields, cooling_time) != GR_SUCCESS) { + std::fprintf(stderr, "Error in local_calculate_cooling_time.\n"); + return GR_FAIL; } - return SUCCESS; + return GR_SUCCESS; } diff --git a/src/clib/calculate_dust_temperature.c b/src/clib/calculate_dust_temperature.cpp similarity index 50% rename from src/clib/calculate_dust_temperature.c rename to src/clib/calculate_dust_temperature.cpp index 6e4d64d3e..18ca1caea 100644 --- a/src/clib/calculate_dust_temperature.c +++ b/src/clib/calculate_dust_temperature.cpp @@ -11,27 +11,22 @@ / software. ************************************************************************/ -#include -#include -#include +#include #include "calc_tdust_3d.h" #include "grackle.h" -#include "grackle_macros.h" #include "internal_units.h" -#include "phys_constants.h" -int local_calculate_dust_temperature(chemistry_data *my_chemistry, - chemistry_data_storage *my_rates, - code_units *my_units, - grackle_field_data *my_fields, - gr_float *dust_temperature) +extern "C" int local_calculate_dust_temperature( + chemistry_data *my_chemistry, chemistry_data_storage *my_rates, + code_units *my_units, grackle_field_data *my_fields, + gr_float *dust_temperature) { if (!my_chemistry->use_grackle) - return SUCCESS; + return GR_SUCCESS; if (my_chemistry->dust_chemistry < 1 && my_chemistry->h2_on_dust < 1) - return SUCCESS; + return GR_SUCCESS; InternalGrUnits internalu = new_internalu_(my_units); @@ -43,35 +38,36 @@ int local_calculate_dust_temperature(chemistry_data *my_chemistry, /* Compute the size of the fields. */ - int i, dim, size = 1; - for (dim = 0; dim < my_fields->grid_rank; dim++) + int size = 1; + for (int dim = 0; dim < my_fields->grid_rank; dim++) { size *= my_fields->grid_dimension[dim]; + } - gr_float *temperature = malloc(size * sizeof(gr_float)); + gr_float *temperature = new gr_float[size]; if (local_calculate_temperature(my_chemistry, my_rates, my_units, - my_fields, temperature) == FAIL) { - fprintf(stderr, "Error in local_calculate_temperature.\n"); - return FAIL; + my_fields, temperature) != GR_SUCCESS) { + std::fprintf(stderr, "Error in local_calculate_temperature.\n"); + return GR_FAIL; } calc_tdust_3d_g( temperature, dust_temperature, metal_field_present, my_chemistry, my_rates, my_fields, internalu ); - free(temperature); + delete[] temperature; - return SUCCESS; + return GR_SUCCESS; } -int calculate_dust_temperature(code_units *my_units, - grackle_field_data *my_fields, - gr_float *dust_temperature) +extern "C" int calculate_dust_temperature(code_units *my_units, + grackle_field_data *my_fields, + gr_float *dust_temperature) { if (local_calculate_dust_temperature( grackle_data, &grackle_rates, my_units, - my_fields, dust_temperature) == FAIL) { - fprintf(stderr, "Error in local_calculate_dust_temperature.\n"); - return FAIL; + my_fields, dust_temperature) != GR_SUCCESS) { + std::fprintf(stderr, "Error in local_calculate_dust_temperature.\n"); + return GR_FAIL; } - return SUCCESS; + return GR_SUCCESS; } diff --git a/src/clib/calculate_gamma.c b/src/clib/calculate_gamma.cpp similarity index 70% rename from src/clib/calculate_gamma.c rename to src/clib/calculate_gamma.cpp index bde450dc5..f80cd0f1d 100644 --- a/src/clib/calculate_gamma.c +++ b/src/clib/calculate_gamma.cpp @@ -11,36 +11,24 @@ / software. ************************************************************************/ -#include -#include -#include -#include "grackle_macros.h" -#include "grackle_types.h" -#include "grackle_chemistry_data.h" +#include +#include +#include "grackle.h" #include "phys_constants.h" #include "index_helper.h" #ifdef _OPENMP #include #endif -extern chemistry_data *grackle_data; -extern chemistry_data_storage grackle_rates; - -int local_calculate_temperature(chemistry_data *my_chemistry, - chemistry_data_storage *my_rates, - code_units *my_units, - grackle_field_data *my_fields, - gr_float *temperature); - -int local_calculate_gamma(chemistry_data *my_chemistry, - chemistry_data_storage *my_rates, - code_units *my_units, - grackle_field_data *my_fields, - gr_float *my_gamma) +extern "C" int local_calculate_gamma(chemistry_data *my_chemistry, + chemistry_data_storage *my_rates, + code_units *my_units, + grackle_field_data *my_fields, + gr_float *my_gamma) { if (!my_chemistry->use_grackle) - return SUCCESS; + return GR_SUCCESS; const grackle_index_helper ind_helper = build_index_helper_(my_fields); int outer_ind, index; @@ -63,9 +51,9 @@ int local_calculate_gamma(chemistry_data *my_chemistry, /* Compute the temperature first. */ if (local_calculate_temperature(my_chemistry, my_rates, my_units, - my_fields, my_gamma) == FAIL) { - fprintf(stderr, "Error in local_calculate_temperature.\n"); - return FAIL; + my_fields, my_gamma) != GR_SUCCESS) { + std::fprintf(stderr, "Error in local_calculate_temperature.\n"); + return GR_FAIL; } /* Compute Gamma with molecular Hydrogen formula from Omukau \& Nishi @@ -107,7 +95,7 @@ int local_calculate_gamma(chemistry_data *my_chemistry, if (nH2 / number_density > 1e-3) { x = 6100.0 / my_gamma[index]; if (x < 10.0) - GammaH2Inverse = 0.5*(5 + 2.0 * x*x * exp(x)/POW(exp(x)-1.0,2)); + GammaH2Inverse = 0.5*(5 + 2.0 * x*x * std::exp(x)/std::pow(std::exp(x)-1.0,2.)); } /* Add in H2. */ @@ -120,17 +108,17 @@ int local_calculate_gamma(chemistry_data *my_chemistry, } // end: if (my_chemistry->primordial_chemistry > 1) - return SUCCESS; + return GR_SUCCESS; } -int calculate_gamma(code_units *my_units, - grackle_field_data *my_fields, - gr_float *my_gamma) +extern "C" int calculate_gamma(code_units *my_units, + grackle_field_data *my_fields, + gr_float *my_gamma) { if (local_calculate_gamma(grackle_data, &grackle_rates, my_units, - my_fields, my_gamma) == FAIL) { + my_fields, my_gamma) != GR_SUCCESS) { fprintf(stderr, "Error in local_calculate_gamma.\n"); - return FAIL; + return GR_FAIL; } - return SUCCESS; + return GR_SUCCESS; } diff --git a/src/clib/calculate_pressure.c b/src/clib/calculate_pressure.cpp similarity index 78% rename from src/clib/calculate_pressure.c rename to src/clib/calculate_pressure.cpp index 2f2aa719a..8c92bcc0c 100644 --- a/src/clib/calculate_pressure.c +++ b/src/clib/calculate_pressure.cpp @@ -11,32 +11,24 @@ / software. ************************************************************************/ -#include -#include -#include -#include "grackle_macros.h" -#include "grackle_types.h" -#include "grackle_chemistry_data.h" +#include +#include // std::fmax +#include "grackle.h" #include "phys_constants.h" #include "index_helper.h" #ifdef _OPENMP #include #endif -extern chemistry_data *grackle_data; -extern chemistry_data_storage grackle_rates; - -double get_temperature_units(code_units *my_units); - -int local_calculate_pressure(chemistry_data *my_chemistry, - chemistry_data_storage *my_rates, - code_units *my_units, - grackle_field_data *my_fields, - gr_float *pressure) +extern "C" int local_calculate_pressure(chemistry_data *my_chemistry, + chemistry_data_storage *my_rates, + code_units *my_units, + grackle_field_data *my_fields, + gr_float *pressure) { if (!my_chemistry->use_grackle) - return SUCCESS; + return GR_SUCCESS; double tiny_number = 1.e-20; const grackle_index_helper ind_helper = build_index_helper_(my_fields); @@ -100,7 +92,7 @@ int local_calculate_pressure(chemistry_data *my_chemistry, if (number_density == 0) number_density = tiny_number; - temp = max(temperature_units * pressure[index] / (number_density + nH2), + temp = std::fmax(temperature_units * pressure[index] / (number_density + nH2), 1); /* Only do full computation if there is a reasonable amount of H2. @@ -126,17 +118,17 @@ int local_calculate_pressure(chemistry_data *my_chemistry, } // end: if (my_chemistry->primordial_chemistry > 1) - return SUCCESS; + return GR_SUCCESS; } -int calculate_pressure(code_units *my_units, - grackle_field_data *my_fields, - gr_float *pressure) +extern "C" int calculate_pressure(code_units *my_units, + grackle_field_data *my_fields, + gr_float *pressure) { if (local_calculate_pressure(grackle_data, &grackle_rates, my_units, - my_fields, pressure) == FAIL) { - fprintf(stderr, "Error in local_calculate_pressure.\n"); - return FAIL; + my_fields, pressure) != GR_SUCCESS) { + std::fprintf(stderr, "Error in local_calculate_pressure.\n"); + return GR_FAIL; } - return SUCCESS; + return GR_SUCCESS; } diff --git a/src/clib/calculate_temperature.c b/src/clib/calculate_temperature.cpp similarity index 65% rename from src/clib/calculate_temperature.c rename to src/clib/calculate_temperature.cpp index 4577348e8..d881fdae0 100644 --- a/src/clib/calculate_temperature.c +++ b/src/clib/calculate_temperature.cpp @@ -11,14 +11,11 @@ / software. ************************************************************************/ -#include -#include -#include +#include #include "calc_temp_cloudy_g.h" #include "grackle.h" #include "index_helper.h" #include "internal_units.h" -#include "phys_constants.h" #ifdef _OPENMP #include @@ -33,11 +30,11 @@ #define MINIMUM_TEMPERATURE 1.0 -int local_calculate_temperature(chemistry_data *my_chemistry, - chemistry_data_storage *my_rates, - code_units *my_units, - grackle_field_data *my_fields, - gr_float *temperature) +extern "C" int local_calculate_temperature(chemistry_data *my_chemistry, + chemistry_data_storage *my_rates, + code_units *my_units, + grackle_field_data *my_fields, + gr_float *temperature) { if (!my_chemistry->use_grackle) { return GR_SUCCESS; } @@ -55,44 +52,40 @@ int local_calculate_temperature(chemistry_data *my_chemistry, /* Compute the pressure first. */ if (local_calculate_pressure(my_chemistry, my_rates, my_units, my_fields, temperature) != GR_SUCCESS) { - fprintf(stderr, "Error in calculate_pressure.\n"); + std::fprintf(stderr, "Error in calculate_pressure.\n"); return GR_FAIL; } - /* Calculate temperature units. */ + // Calculate temperature units and fetch some constants - double temperature_units = get_temperature_units(my_units); - - double number_density, tiny_number = 1.-20; - double inv_metal_mol = 1.0 / MU_METAL; + const double temperature_units = get_temperature_units(my_units); + const double tiny_number = 1.-20; + const double inv_metal_mol = 1.0 / MU_METAL; /* Compute properties used to index the field. */ const grackle_index_helper ind_helper = build_index_helper_(my_fields); - int outer_ind, index; /* Compute temperature with mu calculated directly. */ /* parallelize the k and j loops with OpenMP * (these loops are flattened them for better parallelism) */ # ifdef _OPENMP -# pragma omp parallel for schedule( runtime ) \ - private( outer_ind, index, number_density ) +# pragma omp parallel for schedule( runtime ) # endif - for (outer_ind = 0; outer_ind < ind_helper.outer_ind_size; outer_ind++){ + for (int outer_ind = 0; outer_ind < ind_helper.outer_ind_size; outer_ind++){ const field_flat_index_range range = inner_flat_range_(outer_ind, &ind_helper); - for (index = range.start; index <= range.end; index++) { - - if (my_chemistry->primordial_chemistry > 0) { - number_density = - 0.25 * (my_fields->HeI_density[index] + - my_fields->HeII_density[index] + - my_fields->HeIII_density[index]) + - my_fields->HI_density[index] + my_fields->HII_density[index] + - my_fields->e_density[index]; - } + for (int index = range.start; index <= range.end; index++) { + + // we will only be in this loop if primordial_chemistry > 0 + double number_density = + 0.25 * (my_fields->HeI_density[index] + + my_fields->HeII_density[index] + + my_fields->HeIII_density[index]) + + my_fields->HI_density[index] + my_fields->HII_density[index] + + my_fields->e_density[index]; /* Add in H2. */ @@ -118,13 +111,13 @@ int local_calculate_temperature(chemistry_data *my_chemistry, } -int calculate_temperature(code_units *my_units, - grackle_field_data *my_fields, - gr_float *temperature) +extern "C" int calculate_temperature(code_units *my_units, + grackle_field_data *my_fields, + gr_float *temperature) { if (local_calculate_temperature(grackle_data, &grackle_rates, my_units, my_fields, temperature) != GR_SUCCESS) { - fprintf(stderr, "Error in local_calculate_temperature.\n"); + std::fprintf(stderr, "Error in local_calculate_temperature.\n"); return GR_FAIL; } return GR_SUCCESS; diff --git a/src/clib/solve_chemistry.c b/src/clib/solve_chemistry.cpp similarity index 64% rename from src/clib/solve_chemistry.c rename to src/clib/solve_chemistry.cpp index 130293bb8..15491dfdd 100644 --- a/src/clib/solve_chemistry.c +++ b/src/clib/solve_chemistry.cpp @@ -11,33 +11,30 @@ / software. ************************************************************************/ -#include -#include +#include #include "grackle.h" -#include "grackle_macros.h" #include "internal_units.h" -#include "phys_constants.h" #include "solve_rate_cool_g-cpp.h" #include "utils.h" /* function prototypes */ -int update_UVbackground_rates(chemistry_data *my_chemistry, - chemistry_data_storage *my_rates, - photo_rate_storage *my_uvb_rates, - code_units *my_units); +extern "C" int update_UVbackground_rates(chemistry_data *my_chemistry, + chemistry_data_storage *my_rates, + photo_rate_storage *my_uvb_rates, + code_units *my_units); -int local_solve_chemistry(chemistry_data *my_chemistry, - chemistry_data_storage *my_rates, - code_units *my_units, - grackle_field_data *my_fields, - double dt_value) +extern "C" int local_solve_chemistry(chemistry_data *my_chemistry, + chemistry_data_storage *my_rates, + code_units *my_units, + grackle_field_data *my_fields, + double dt_value) { /* Return if this doesn't concern us. */ if (!my_chemistry->use_grackle) - return SUCCESS; + return GR_SUCCESS; /* Update UV background rates. */ photo_rate_storage my_uvb_rates; @@ -51,9 +48,9 @@ int local_solve_chemistry(chemistry_data *my_chemistry, if (my_chemistry->UVbackground == 1) { if (update_UVbackground_rates(my_chemistry, my_rates, - &my_uvb_rates, my_units) == FAIL) { - fprintf(stderr, "Error in update_UVbackground_rates.\n"); - return FAIL; + &my_uvb_rates, my_units) != GR_SUCCESS) { + std::fprintf(stderr, "Error in update_UVbackground_rates.\n"); + return GR_FAIL; } } else { @@ -85,8 +82,8 @@ int local_solve_chemistry(chemistry_data *my_chemistry, /* Error checking for H2 shielding approximation */ if (self_shielding_err_check(my_chemistry, my_fields, - "local_solve_chemistry") == FAIL) { - return FAIL; + "local_solve_chemistry") != GR_SUCCESS) { + return GR_SUCCESS; } /* Call the routine to solve cooling equations. */ @@ -96,22 +93,22 @@ int local_solve_chemistry(chemistry_data *my_chemistry, my_chemistry, my_rates, my_fields, &my_uvb_rates ); - if (ierr == GR_FAIL) { - fprintf(stderr, "Error in solve_rate_cool_g.\n"); + if (ierr != GR_SUCCESS) { + std::fprintf(stderr, "Error in solve_rate_cool_g.\n"); } return ierr; } -int solve_chemistry(code_units *my_units, - grackle_field_data *my_fields, - double dt_value) +extern "C" int solve_chemistry(code_units *my_units, + grackle_field_data *my_fields, + double dt_value) { if (local_solve_chemistry(grackle_data, &grackle_rates, - my_units, my_fields, dt_value) == FAIL) { - fprintf(stderr, "Error in local_solve_chemistry.\n"); - return FAIL; + my_units, my_fields, dt_value) != GR_SUCCESS) { + std::fprintf(stderr, "Error in local_solve_chemistry.\n"); + return GR_FAIL; } - return SUCCESS; + return GR_SUCCESS; } diff --git a/src/clib/utils.h b/src/clib/utils.h index 003794bee..dc74e372e 100644 --- a/src/clib/utils.h +++ b/src/clib/utils.h @@ -11,9 +11,16 @@ / software. ************************************************************************/ +#ifndef UTILS_H +#define UTILS_H + #include "grackle_chemistry_data.h" #include "grackle_types.h" +#ifdef __cplusplus +extern "C" { +#endif // __cplusplus + /// Perform an error check related to self-shielding. If the check fails, the /// function returns FAIL and prints an error message to stderr /// @@ -23,3 +30,10 @@ int self_shielding_err_check(const chemistry_data *my_chemistry, const grackle_field_data *fields, const char* func_name); + + +#ifdef __cplusplus +} // extern "C" +#endif // __cplusplus + +#endif // UTILS_H