diff --git a/src/MINT/MINT_shellular_burning.c b/src/MINT/MINT_shellular_burning.c
index f15d09d4184be4827e42f181aff8b91cf8b98ac1..3f22969c3e239aa24502a3bffecd6c8e6528f9b5 100644
--- a/src/MINT/MINT_shellular_burning.c
+++ b/src/MINT/MINT_shellular_burning.c
@@ -90,7 +90,6 @@ void MINT_shellular_burning(struct stardata_t * RESTRICT const stardata,
convective_sigmav[j] /= mconv_squared;
}
-
/*
* Finalise the Xconv calculation
*/
@@ -119,6 +118,8 @@ void MINT_shellular_burning(struct stardata_t * RESTRICT const stardata,
nucsyn_burning_cycles(Nconv,
convective_sigmav,
stardata->model.dt, /* dt in years */
+ star->mint->shells[0].T,
+ star->mint->shells[0].rho,
stardata);
Nconv[Xe] = 0.0;
@@ -170,6 +171,8 @@ void MINT_shellular_burning(struct stardata_t * RESTRICT const stardata,
nucsyn_burning_cycles(Nshell,
sigmav,
stardata->model.dt, /* dt in years */
+ shell->T,
+ shell->rho,
stardata);
Nshell[Xe] = 0.0;
nancheck("MINT burn (radiative N) ",Nshell,ISOTOPE_ARRAY_SIZE);
diff --git a/src/binary_c.h b/src/binary_c.h
index 3657e39d816c685b3d65ae7ac2fa4a6a97eea31b..17433903635e7b67ab545ed8908332df41dfae2a 100644
--- a/src/binary_c.h
+++ b/src/binary_c.h
@@ -60,8 +60,6 @@
/* binary_c's headers */
#include "binary_c_version_macros.h"
#include "binary_c_code_options.h"
-
-
#include "binary_c_debug.h"
#include "binary_c_error_codes.h"
#include "misc/misc.h"
diff --git a/src/binary_c_debug.h b/src/binary_c_debug.h
index acd42bea4642a0e90243b3dbca655021c6d1dfc6..0d288a0d4481bcb269c085a1448fde980a325d30 100644
--- a/src/binary_c_debug.h
+++ b/src/binary_c_debug.h
@@ -22,6 +22,7 @@
#define BINARY_C_DEBUG_H
#include "binary_c_string_macros.h"
+#include "logging/ansi_colours.h"
/*
* Do you want debugging output?
diff --git a/src/binary_c_structures.h b/src/binary_c_structures.h
index 8e2f3a58d34c1c132f2c898acaa1a4867f28f710..03947bb09d28ac09f76544188dba65298ce6608d 100644
--- a/src/binary_c_structures.h
+++ b/src/binary_c_structures.h
@@ -259,6 +259,7 @@ struct store_t {
char ** isotope_strings;/* strings describing each element */
struct nuclear_network_range_t nuclear_network_range[NUCSYN_NETWORK_NUMBER];
nucsyn_burnfunc nucsyn_network_function[NUCSYN_SOLVER_NUMBER][NUCSYN_NETWORK_NUMBER];
+ nucsyn_checkfunc nucsyn_network_checkfunction[NUCSYN_NETWORK_NUMBER];
#endif // NUCSYN
@@ -1208,6 +1209,8 @@ struct preferences_t {
Stellar_type MINT_Kippenhahn_companion_stellar_type;
#endif // MINT
+ int float_overflow_checks;
+
Stellar_type max_stellar_type[4];
Stellar_structure_algorithm stellar_structure_algorithm;
@@ -1218,10 +1221,8 @@ struct preferences_t {
/* Boolean parameters */
Boolean tasks[BINARY_C_TASK_NUMBER];
-
Boolean require_drdm;
Boolean save_pre_events_stardata;
-
Boolean PN_fast_wind;
Boolean PN_resolve;
diff --git a/src/debug/debug_colours.h b/src/debug/debug_colours.h
index 8e887697d293ce28d3b258f0580ffffaef364bc8..e7e2b9901f49de4a18c08bdf90cf0cfde897a6f7 100644
--- a/src/debug/debug_colours.h
+++ b/src/debug/debug_colours.h
@@ -8,12 +8,14 @@
#define COLOUR_ENDLINE COLOUR_RESET
#define COLOUR_FILENAME RED
#define COLOUR_CALLER YELLOW
+#define COLOUR_FLOAT_BITSTRING MAGENTA
#else
#define COLOUR_TIME ""
#define COLOUR_LINENUMBERS ""
#define COLOUR_ENDLINE ""
#define COLOUR_FILENAME ""
#define COLOUR_CALLER ""
+#define COLOUR_FLOAT_BITSTRING ""
#endif
#endif // DEBUG_COLOURS_H
diff --git a/src/debug/debug_fprintf.c b/src/debug/debug_fprintf.c
index 21f3f9ae39207bf5c6858db8a56bf5ae8b97f841..90ccbc72086b32246764edd92043325cdb0cf3dd 100644
--- a/src/debug/debug_fprintf.c
+++ b/src/debug/debug_fprintf.c
@@ -147,36 +147,43 @@ void Gnu_format_args(5,6) debug_fprintf(struct stardata_t * const stardata,
caller[0] = '\0';
#endif // DEBUG_CALLER && BACKTRACE
+ char * float_test_bitstring =
+ floating_point_exception_string();
+
/*
* Make output string
*/
char out[MAX_DEBUG_PRINT_SIZE];
#ifdef DEBUG_LINENUMBERS
snprintf_returned = snprintf(out,
- MAX_DEBUG_PRINT_SIZE,
- "%s%s%20.12g %d %s%s% 6d %s :%s : %s\n",
- sid,
- COLOUR_TIME,
- t,
- stardata!=NULL ? stardata->model.intpol : 0,
- f,
- COLOUR_LINENUMBERS,
- fileline,
- caller,
- COLOUR_ENDLINE,
- s);
+ MAX_DEBUG_PRINT_SIZE,
+ "%s%s%s%s%20.12g %d %s%s% 6d %s :%s : %s\n",
+ sid,
+ COLOUR_FLOAT_BITSTRING,
+ float_test_bitstring,
+ COLOUR_TIME,
+ t,
+ stardata!=NULL ? stardata->model.intpol : 0,
+ f,
+ COLOUR_LINENUMBERS,
+ fileline,
+ caller,
+ COLOUR_ENDLINE,
+ s);
#else
snprintf_returned = snprintf(out,
- MAX_DEBUG_PRINT_SIZE,
- "%s%s%20.12g %d %s%s%s : %s\n",
- sid,
- COLOUR_TIME,
- t,
- stardata!=NULL ? stardata->model.intpol : 0,
- f,
- caller,
- COLOUR_ENDLINE,
- s);
+ MAX_DEBUG_PRINT_SIZE,
+ "%s%s%s%s%20.12g %d %s%s%s : %s\n",
+ sid,
+ COLOUR_FLOAT_BITSTRING,
+ float_test_bitstring,
+ COLOUR_TIME,
+ t,
+ stardata!=NULL ? stardata->model.intpol : 0,
+ f,
+ caller,
+ COLOUR_ENDLINE,
+ s);
#endif //DEBUG_LINENUMBERS
if(snprintf_returned < 0)
diff --git a/src/evolution/evolve_system_binary_c.c b/src/evolution/evolve_system_binary_c.c
index 836a42ea74dc73862cb9361560abe83797c2c3bf..0e9891744a7dda7ed1e1ea9a6a1e698cadbd7f09 100644
--- a/src/evolution/evolve_system_binary_c.c
+++ b/src/evolution/evolve_system_binary_c.c
@@ -71,6 +71,11 @@ int evolve_system_binary_c(struct stardata_t * RESTRICT const stardata)
while(evolving == TRUE)
{
+ if(stardata->preferences->float_overflow_checks > 0)
+ {
+ floating_point_clear_flags(stardata);
+ }
+
stardata->model.model_number++;
eprint("Evolution loop top : model %d\n",
stardata->model.model_number);
@@ -262,6 +267,14 @@ int evolve_system_binary_c(struct stardata_t * RESTRICT const stardata)
catch_events(stardata);
}
+ /*
+ * Check for floating-point exceptions
+ */
+ if(stardata->preferences->float_overflow_checks > 0)
+ {
+ floating_point_exception_checks(stardata);
+ }
+
/*
* Check for evolution split before
* we save_to_previous. If evolution_split
diff --git a/src/evolution/initialize_generic_system_every_timestep.c b/src/evolution/initialize_generic_system_every_timestep.c
index 4cab8e997668dd4cc5fad32c5dc1f2406b1820fa..02fbc30603556e4594d5f3ee5d7f2329e414b41f 100644
--- a/src/evolution/initialize_generic_system_every_timestep.c
+++ b/src/evolution/initialize_generic_system_every_timestep.c
@@ -27,6 +27,7 @@ void initialize_generic_system_every_timestep(struct stardata_t * RESTRICT const
initialize_generic_system_every_timestep_star(star);
}
}
+ Dprint("init done\n");
}
void initialize_generic_system_every_timestep2(struct stardata_t * RESTRICT const stardata)
diff --git a/src/logging/ansi_colours.h b/src/logging/ansi_colours.h
index c7938b06ab60d06de3fcd21c676441ec58821c9c..83251ee464dec9a12975dd3cfd9009ad9f386e19 100644
--- a/src/logging/ansi_colours.h
+++ b/src/logging/ansi_colours.h
@@ -33,6 +33,7 @@
#define COLOUR_RESET "\x1B[0;0m"
#define DEBUG_COLOURS
+
#define FILE_LOG_COLOURS
#else//USE_ANSI_COLOURS
#define BLACK ""
diff --git a/src/maths/floating_point_clear_flags.c b/src/maths/floating_point_clear_flags.c
new file mode 100644
index 0000000000000000000000000000000000000000..60f7c0db09e06daf3ce3414739ed3b5062419bfd
--- /dev/null
+++ b/src/maths/floating_point_clear_flags.c
@@ -0,0 +1,12 @@
+#include "../binary_c.h"
+No_empty_translation_unit_warning;
+#include <fenv.h>
+
+/*
+ * Clear floating-point exception flags
+ */
+
+void floating_point_clear_flags(struct stardata_t * const stardata MAYBE_UNUSED)
+{
+ feclearexcept(FE_ALL_EXCEPT);
+}
diff --git a/src/maths/floating_point_exception_checks.c b/src/maths/floating_point_exception_checks.c
new file mode 100644
index 0000000000000000000000000000000000000000..9c5f7981ffb2592f01b5bfe1a3f59a1af52bded0
--- /dev/null
+++ b/src/maths/floating_point_exception_checks.c
@@ -0,0 +1,25 @@
+#include "../binary_c.h"
+No_empty_translation_unit_warning;
+
+/*
+ * Check for floating-point exceptions after each timestep
+ */
+
+void floating_point_exception_checks(struct stardata_t * const stardata)
+{
+ if(fetestexcept(FE_ALL_EXCEPT))
+ {
+ /*
+ * Exception found!
+ */
+ if(stardata->preferences->float_overflow_checks == 1)
+ {
+ fprintf(stderr,"Floating point exception detected\n");
+ }
+ else
+ {
+ Exit_binary_c(BINARY_C_FLOATING_POINT_ERROR,
+ "Floating point exception detected\n");
+ }
+ }
+}
diff --git a/src/maths/floating_point_exception_string.c b/src/maths/floating_point_exception_string.c
new file mode 100644
index 0000000000000000000000000000000000000000..e138b299bd37236e4378f5dc172b8b30154128fd
--- /dev/null
+++ b/src/maths/floating_point_exception_string.c
@@ -0,0 +1,47 @@
+#include "../binary_c.h"
+No_empty_translation_unit_warning;
+
+/*
+ * Make floating point exception bitstring,
+ * useful for debugging
+ */
+char * floating_point_exception_string(void)
+{
+ unsigned int n = 0;
+#ifdef FE_INEXACT
+ n++;
+#endif
+#ifdef FE_DIVBYZERO
+ n++;
+#endif
+#ifdef FE_UNDERFLOW
+ n++;
+#endif
+#ifdef FE_OVERFLOW
+ n++;
+#endif
+#ifdef FE_INVALID
+ n++;
+#endif
+ char * const s = Malloc(sizeof(char)*(n+1));
+ n=0;
+#undef _TEST
+#define _TEST(X) s[n++] = (fetestexcept(X) ? '!' : '.')
+#ifdef FE_INEXACT
+ _TEST(FE_INEXACT);
+#endif
+#ifdef FE_DIVBYZERO
+ _TEST(FE_DIVBYZERO);
+#endif
+#ifdef FE_UNDERFLOW
+ _TEST(FE_UNDERFLOW);
+#endif
+#ifdef FE_OVERFLOW
+ _TEST(FE_OVERFLOW);
+#endif
+#ifdef FE_INVALID
+ _TEST(FE_INVALID);
+#endif
+ s[n] = '\0';
+ return s;
+}
diff --git a/src/maths/kaps_rentrop_GSL.c b/src/maths/kaps_rentrop_GSL.c
index 5d985601bf3e49c7a4be5e12d1e359773ee74763..d8fdf323332af91244dff4e6fe0e0feb2b5929f7 100644
--- a/src/maths/kaps_rentrop_GSL.c
+++ b/src/maths/kaps_rentrop_GSL.c
@@ -115,9 +115,9 @@ double kaps_rentrop_GSL(struct stardata_t * RESTRICT const stardata,
double d; // dummy variable
const double ih = 1.0 / h; // saves multiple calculations of 1/h
- // used in matrix calculations (Numerical Recipes code)
+ // used in matrix calculations
const unsigned int nn = n + 1; // size of arrays
- const size_t nnd = nn * sizeof(double); // speedup
+ const size_t nnd = nn * sizeof(double);
#define m stardata->tmpstore->kaps_rentrop_GSL_m
#define p stardata->tmpstore->kaps_rentrop_GSL_p
diff --git a/src/maths/maths_prototypes.h b/src/maths/maths_prototypes.h
index 620de1012cea3f3172e7cbc085b21165fde91281..90a5d955de3c9841de482a796577bbe2758f7ff9 100644
--- a/src/maths/maths_prototypes.h
+++ b/src/maths/maths_prototypes.h
@@ -295,4 +295,9 @@ long strtol10(char * ptr,
double exp10(double x);
#endif // __HAVE_NATIVE_EXP10
+
+void floating_point_exception_checks(struct stardata_t * const stardata);
+void floating_point_clear_flags(struct stardata_t * const stardata);
+char * floating_point_exception_string(void);
+
#endif // MATHS_PROTOTYPES
diff --git a/src/nucsyn/nucsyn_build_store_contents.c b/src/nucsyn/nucsyn_build_store_contents.c
index c7c4d8a52f376da1678377d5ee96fc95a11c4e5c..5da6e1838bfca4f12656baafccdb0eee183dbc9c 100644
--- a/src/nucsyn/nucsyn_build_store_contents.c
+++ b/src/nucsyn/nucsyn_build_store_contents.c
@@ -83,12 +83,14 @@ void nucsyn_build_store_contents(struct store_t * RESTRICT const store)
STRING, \
KAPS_RENTROP_FUNCTION, \
LSODA_FUNCTION, \
- CVODE_FUNCTION) \
+ CVODE_FUNCTION, \
+ CHECKFUNC) \
if(NUCSYN_NETWORK_##CODE != NUCSYN_NETWORK_NUMBER) \
{ \
store->nucsyn_network_function[NUCSYN_SOLVER_KAPS_RENTROP][NUCSYN_NETWORK_##CODE] = KAPS_RENTROP_FUNCTION; \
store->nucsyn_network_function[NUCSYN_SOLVER_LSODA][NUCSYN_NETWORK_##CODE] = LSODA_FUNCTION; \
store->nucsyn_network_function[NUCSYN_SOLVER_CVODE][NUCSYN_NETWORK_##CODE] = _CVODE(CVODE_FUNCTION); \
+ store->nucsyn_network_checkfunction[NUCSYN_NETWORK_##CODE] = CHECKFUNC; \
}
NUCSYN_NETWORKS_LIST;
#undef X
diff --git a/src/nucsyn/nucsyn_burn.c b/src/nucsyn/nucsyn_burn.c
index 480e6f531b1dcfa6d3e012b6bddb222706a75640..76b0a6781be3ed76b1915473ae3b6bcca70f13d7 100644
--- a/src/nucsyn/nucsyn_burn.c
+++ b/src/nucsyn/nucsyn_burn.c
@@ -1,7 +1,7 @@
#include "../binary_c.h"
No_empty_translation_unit_warning;
-#if (defined NUCSYN && defined NUCSYN_TPAGB_HBB)
+#ifdef NUCSYN
#include <string.h>
#include "nucsyn_isotopes.h"
@@ -17,9 +17,10 @@ No_empty_translation_unit_warning;
/*
* dt is the timestep (years), N is the number density (cgs), passed in is the
* old number density, passed out is the new. oN stores the old values of the
- * number density
+ * number density. logT = log10(temperature/K), density is in g/cm^3.
*/
-void nucsyn_burn(const double T,
+void nucsyn_burn(const double logT,
+ const double density,
const double dt,
Number_density * RESTRICT const N,
struct star_t * const star MAYBE_UNUSED,
@@ -38,7 +39,7 @@ void nucsyn_burn(const double T,
N[Xe] = nucsyn_free_electron_density(N,stardata->store->atomic_number);
nucsyn_set_sigmav(stardata,
stardata->store,
- T,
+ logT,
sigmav,
stardata->preferences->reaction_rate_multipliers
#ifdef RATES_OF_AMANDA_NE22PG_FIX
@@ -63,7 +64,7 @@ void nucsyn_burn(const double T,
#endif
/* call the burning cycles (dt in years) */
- nucsyn_burning_cycles(N,sigmav,dt,stardata);
+ nucsyn_burning_cycles(N,sigmav,dt,exp10(logT),density,stardata);
Dprint("nucsyn_burn: post-burning cycles %g %g\n",N[0],N[1]);
@@ -144,7 +145,7 @@ void nucsyn_burn(const double T,
nucsyn_set_sigmav(stardata,
stardata->store,
- T,
+ logT,
sigmav,
stardata->preferences->reaction_rate_multipliers
#ifdef RATES_OF_AMANDA_NE22PG_FIX
@@ -363,7 +364,7 @@ void nucsyn_burn(const double T,
* Check the temperature here, if it's too low (<1e6K)
* then do NOT activate the NeNa and MgAl cycles
*/
- if(T<6)
+ if(logT<6)
{
return;
}
@@ -376,7 +377,7 @@ void nucsyn_burn(const double T,
#ifdef NUCSYN_TPAGB_HBB_NeNa_CYCLE
/* Analytic form of the NeNa cycle */
/* Set up tau_20,21,22,23 */
- // if(T>T_MIN_NeNa_HBB)
+ // if(logT>T_MIN_NeNa_HBB)
if(1==1)
{
timescales[0]=Nuclear_burning_timescale(SIGMAV_F19A);
@@ -387,7 +388,7 @@ void nucsyn_burn(const double T,
timescales[5]=Nuclear_burning_timescale(SIGMAV_T23P);
Dprint("NenaT T9=%g (T=%g) tscls %g %g %g %g leak %g\n",
- exp10(T)/1e9,T,
+ exp10(T)/1e9,logT,
timescales[1],
timescales[2],
timescales[3],
@@ -469,7 +470,7 @@ void nucsyn_burn(const double T,
{
if(sigmavset==FALSE) nucsyn_set_sigmav(stardata,
stardata->store,
- T,
+ logT,
sigmav,
stardata->preferences->reaction_rate_multipliers
#ifdef RATES_OF_AMANDA_NE22PG_FIX
@@ -597,7 +598,7 @@ void nucsyn_burn(const double T,
-#endif /* NUCSYN_TPAGB_HBB */
+#endif // NUCSYN
diff --git a/src/nucsyn/nucsyn_burning_cycles.c b/src/nucsyn/nucsyn_burning_cycles.c
index 8d14adeea591583fe7e1a5b928d05719aa279620..d3edc5130a5a3ff64c8fa360e12775bf3d7f42ae 100644
--- a/src/nucsyn/nucsyn_burning_cycles.c
+++ b/src/nucsyn/nucsyn_burning_cycles.c
@@ -10,6 +10,8 @@ No_empty_translation_unit_warning;
void nucsyn_burning_cycles(Abundance * RESTRICT const N,
Reaction_rate * RESTRICT const sigmav,
const double dt,
+ const double temperature,
+ const double density,
struct stardata_t * RESTRICT stardata)
{
#include "nucsyn_solver_functions.h"
@@ -60,10 +62,25 @@ void nucsyn_burning_cycles(Abundance * RESTRICT const N,
*/
for(Reaction_network i=0; i<NUCSYN_NETWORK_NUMBER; i++)
{
- Dprint("Network %u \"%s\" : burn? %s : for time %g\n",
+
+ /*
+ * If we have a function which tests T,rho to see
+ * if they are in the appropriate range for burning, use it,
+ * otherwise assume we must burn
+ */
+ const Boolean test_T_rho =
+ stardata->store->nucsyn_network_checkfunction[i] != NULL ?
+ stardata->store->nucsyn_network_checkfunction[i](stardata,
+ N,
+ temperature,
+ density) :
+ TRUE;
+
+ Dprint("Network %u \"%s\" : burn? %s && %s : for time %g\n",
i,
nuclear_network_strings[i],
Yesno(stardata->preferences->nucsyn_network[i]),
+ Yesno(test_T_rho),
dt);
#ifdef TESTING
@@ -83,7 +100,8 @@ void nucsyn_burning_cycles(Abundance * RESTRICT const N,
const double time_step = dt;
#endif
- if(stardata->preferences->nucsyn_network[i] == TRUE)
+ if(test_T_rho == TRUE &&
+ stardata->preferences->nucsyn_network[i] == TRUE)
{
#ifdef TESTING
double time = 0.0; // years
diff --git a/src/nucsyn/nucsyn_functions.h b/src/nucsyn/nucsyn_functions.h
index 48131ac6a0eedb2577ec9d4318644a24921d89ed..2162a1e5b2fdb8337cfeef74a089dfb478a477c6 100644
--- a/src/nucsyn/nucsyn_functions.h
+++ b/src/nucsyn/nucsyn_functions.h
@@ -9,6 +9,12 @@ typedef double(*nucsyn_burnfunc)(struct stardata_t * const,
const double * const,
const double);
+/* nucsyn network check function */
+typedef Boolean(*nucsyn_checkfunc)(struct stardata_t * const,
+ double * const,
+ const double,
+ const double);
+
/* nucsyn network logging function */
typedef void(*nucsyn_logfunc)(struct stardata_t * const stardata,
const double * const Nin,
diff --git a/src/nucsyn/nucsyn_hbb.c b/src/nucsyn/nucsyn_hbb.c
index b15649e2b5ae9ce51231ec7b506ad13756b3561c..ecb6b35ed77c23b027cafcb1dea5510e810ba269 100644
--- a/src/nucsyn/nucsyn_hbb.c
+++ b/src/nucsyn/nucsyn_hbb.c
@@ -26,27 +26,27 @@ void nucsyn_hbb(const double T, /* (log10) temperature */
struct stardata_t * const stardata
)
{
- const Nuclear_mass * mnuc = stardata->store->mnuc;
- const Nuclear_mass *imnuc = stardata->store->imnuc;
+ const Nuclear_mass * const mnuc = stardata->store->mnuc;
+ const Nuclear_mass * const imnuc = stardata->store->imnuc;
//#define __TEST_HE_BURN
#ifdef __TEST_HE_BURN
Abundance * X = New_clear_isotope_array;
Number_density * N = New_isotope_array;
+
+ const double timemax = 100;
+ const double dtime = timemax/100.0;
+ const double logtemp = log10(1e8);
const double testrho = 1e3;
/* start with a lot of helium */
X[XHe4] = 1.0;
X_to_N(imnuc,testrho,N,X,ISOTOPE_ARRAY_SIZE);
- const double timemax = 100;
double time = 0.0;
- double dtime = timemax/100.0;
- double logtemp = log10(3e8);
-
while(time < timemax)
{
- nucsyn_burn(logtemp,dtime,N,star,stardata);
+ nucsyn_burn(logtemp,testrho,dtime,N,star,stardata);
N_to_X(mnuc,testrho,N,X,ISOTOPE_ARRAY_SIZE);
printf("HEBURN %g %g %g %g %g %g %g\n",
time,
@@ -117,7 +117,7 @@ void nucsyn_hbb(const double T, /* (log10) temperature */
#endif
/* do the burn ( N(t) -> N(t+dt) ) */
- nucsyn_burn(T,dt,Nhbb,star,stardata);
+ nucsyn_burn(T,rho,dt,Nhbb,star,stardata);
#ifdef NANCHECKS
Isotope_loop(i)
diff --git a/src/nucsyn/nucsyn_helium_checkfunction.c b/src/nucsyn/nucsyn_helium_checkfunction.c
new file mode 100644
index 0000000000000000000000000000000000000000..d862cf832e3bb08ed1bf41d0dbc1cc0481d1c74b
--- /dev/null
+++ b/src/nucsyn/nucsyn_helium_checkfunction.c
@@ -0,0 +1,15 @@
+#include "../binary_c.h"
+No_empty_translation_unit_warning;
+
+Boolean nucsyn_helium_checkfunction(struct stardata_t * const stardata MAYBE_UNUSED,
+ double * const N,
+ const double temperature,
+ const double density MAYBE_UNUSED)
+{
+ /*
+ * Check temperature, density and input abundance
+ * array (number densities): return TRUE if we should
+ * burn helium, FALSE otherwise
+ */
+ return (N[XHe4]>0.0 && temperature >= 0.8e8) ? TRUE : FALSE;
+}
diff --git a/src/nucsyn/nucsyn_networks.def b/src/nucsyn/nucsyn_networks.def
index afc80e37690f7696db0287a9c3c4fcb043599db8..e9881971885664445317d0cbde83893024d036f1 100644
--- a/src/nucsyn/nucsyn_networks.def
+++ b/src/nucsyn/nucsyn_networks.def
@@ -12,6 +12,7 @@
* 4: function to call to do the burning with the Kaps-Rentrop solver
* 5: function to call to do the burning with the LSODA solver
* 6: function to call to do the burning with the CVODE solver
+ * 7: function to test T,rho to determine whether we should attempt a burn (or NULL)
*/
#ifdef __HAVE_LIBSUNDIALS_CVODE__
#define _CVODE(X) (X)
@@ -20,15 +21,17 @@
#endif // __HAVE_LIBSUNDIALS_CVODE__
#define NUCSYN_NETWORKS_LIST \
- X(PP, TRUE, "pp", &nucsyn_burn_kaps_rentrop_pp, &nucsyn_burn_lsoda_pp, _CVODE(&nucsyn_burn_cvode_pp)) \
- X(PPfast, FALSE, "pp fast", &nucsyn_burn_ppfast,NULL,_CVODE(NULL)) \
- X(COLDCNO, TRUE, "Cold CNO",&nucsyn_burn_kaps_rentrop_coldCNO,&nucsyn_burn_lsoda_coldCNO,_CVODE(&nucsyn_burn_cvode_coldCNO)) \
- X(HOTCNO, FALSE, "Hot CNO",&nucsyn_burn_hotCNO,NULL,_CVODE(NULL)) \
- X(NeNa, FALSE, "NeNa",NULL,NULL,_CVODE(NULL)) \
- X(NeNaMgAlnoleak, FALSE, "NeNaMgAl no leak",NULL,NULL,_CVODE(NULL)) \
- X(NeNaMgAl, TRUE, "NeNaMgAl",&nucsyn_burn_kaps_rentrop_NeNaMgAl,&nucsyn_burn_lsoda_coldCNO,_CVODE(&nucsyn_burn_cvode_NeNaMgAl)) \
- X(helium,TRUE,"helium",&nucsyn_burn_helium,NULL,_CVODE(NULL)) \
- X(NUMBER, FALSE, "",NULL,NULL,NULL)
+ X(PP, TRUE, "pp", &nucsyn_burn_kaps_rentrop_pp, &nucsyn_burn_lsoda_pp, _CVODE(&nucsyn_burn_cvode_pp),NULL) \
+ X(PPfast, FALSE, "pp fast", &nucsyn_burn_ppfast,NULL,_CVODE(NULL),NULL) \
+ X(COLDCNO, TRUE, "Cold CNO",&nucsyn_burn_kaps_rentrop_coldCNO,&nucsyn_burn_lsoda_coldCNO,_CVODE(&nucsyn_burn_cvode_coldCNO),NULL) \
+ X(HOTCNO, FALSE, "Hot CNO",&nucsyn_burn_hotCNO,NULL,_CVODE(NULL),NULL) \
+ X(NeNa, FALSE, "NeNa",NULL,NULL,_CVODE(NULL),NULL) \
+ X(NeNaMgAlnoleak, FALSE, "NeNaMgAl no leak",NULL,NULL,_CVODE(NULL),NULL) \
+ X(NeNaMgAl, TRUE, "NeNaMgAl",&nucsyn_burn_kaps_rentrop_NeNaMgAl,&nucsyn_burn_lsoda_coldCNO,_CVODE(&nucsyn_burn_cvode_NeNaMgAl),NULL) \
+ X(helium,TRUE,"helium",&nucsyn_burn_helium,NULL,_CVODE(NULL),&nucsyn_helium_checkfunction) \
+ \
+ /* this last line should not be changed */ \
+ X(NUMBER, FALSE, "",NULL,NULL,NULL,NULL)
diff --git a/src/nucsyn/nucsyn_networks.h b/src/nucsyn/nucsyn_networks.h
index cdb224dfd2d0d8bb22df4301a9acb522bcf26ddb..93d7fa62d341701036b9db6c1738b1a21a3c9544 100644
--- a/src/nucsyn/nucsyn_networks.h
+++ b/src/nucsyn/nucsyn_networks.h
@@ -16,14 +16,14 @@
* Make enum list for the NUCSYN_NETWORK_XXX values
*/
#undef X
-#define X(CODE,ONOFF,STRING,KAPS_RENTROP,LSODA,CVODE) NUCSYN_NETWORK_##CODE,
- enum { NUCSYN_NETWORKS_LIST };
+#define X(CODE,ONOFF,STRING,KAPS_RENTROP,LSODA,CVODE,CHECKFUNC) NUCSYN_NETWORK_##CODE,
+enum { NUCSYN_NETWORKS_LIST };
#undef X
/*
* Make strings for logging
*/
-#define X(CODE,ONOFF,STRING,KAPS_RENTROP,LSODA,CVODE) STRING,
+#define X(CODE,ONOFF,STRING,KAPS_RENTROP,LSODA,CVODE,CHECKFUNC) STRING,
static const char * const nuclear_network_strings[] = { NUCSYN_NETWORKS_LIST };
#undef X
diff --git a/src/nucsyn/nucsyn_prototypes.h b/src/nucsyn/nucsyn_prototypes.h
index 842295885a13cdca7c75a01f851ecbe5dd3766ad..d66837d5183d8fd7be9d2558104b39439603eee6 100644
--- a/src/nucsyn/nucsyn_prototypes.h
+++ b/src/nucsyn/nucsyn_prototypes.h
@@ -382,7 +382,8 @@ void nucsyn_r_process_Arlandini1999(Abundance * RESTRICT const Xr,
/*
* Nuclear network solvers
*/
-void nucsyn_burn(const double T,
+void nucsyn_burn(const double logT,
+ const double density,
const double dt,
Number_density * RESTRICT const N,
struct star_t * const star MAYBE_UNUSED,
@@ -504,7 +505,10 @@ double nucsyn_burn_helium(struct stardata_t * RESTRICT const stardata,
const double * RESTRICT const sigmav,
const double h);
-
+Boolean nucsyn_helium_checkfunction(struct stardata_t * const stardata,
+ double * const N,
+ const double temperature,
+ const double density);
double Constant_function nucsyn_CN_timescale(const double t9,
@@ -533,6 +537,8 @@ double Pure_function nucsyn_free_electron_density(const double * RESTRICT const
void nucsyn_burning_cycles(Abundance * RESTRICT const N,
Reaction_rate * RESTRICT const sigmav,
const double dt,
+ const double temperature,
+ const double density,
struct stardata_t * RESTRICT const stardata);
void nucsyn_reset_reaction_rates(struct preferences_t * RESTRICT const preferences);
diff --git a/src/nucsyn/nucsyn_solvers.h b/src/nucsyn/nucsyn_solvers.h
index 2c7abd046fa94c1485fa7748d9434671cfbabb5e..23914805900bda546f4302d89b3c41c9a364aa3f 100644
--- a/src/nucsyn/nucsyn_solvers.h
+++ b/src/nucsyn/nucsyn_solvers.h
@@ -5,7 +5,9 @@
#undef X
#define X(CODE,STRING,FUNC) NUCSYN_SOLVER_##CODE,
-enum{ NUCSYN_SOLVERS_LIST };
+enum{
+ NUCSYN_SOLVERS_LIST
+};
#undef X
#define X(CODE,STRING,FUNC) STRING,
@@ -13,4 +15,5 @@ static char * nucsyn_solver_strings[] MAYBE_UNUSED = { NUCSYN_SOLVERS_LIST };
#undef X
#undef X
+
#endif // NUCSYN_SOLVERS_H
diff --git a/src/setup/cmd_line_args_list.h b/src/setup/cmd_line_args_list.h
index 2f959e2d0d7c575a5b8beb796aade3fbace94715..563659d69ec16a4dd715ae748c3b1e1f5072ca1b 100644
--- a/src/setup/cmd_line_args_list.h
+++ b/src/setup/cmd_line_args_list.h
@@ -1943,6 +1943,15 @@ BATCHMODE_ARGS
NUCSYN_CF_AMANDA_VAR,
1.0
CMD_LINE_ARG_T_OPTIONAL_ARGS },
+{
+ ARG_SECTION_OUTPUT,
+ "float_overflow_checks",
+ "Turn on to enable floating-point overflow checks at the end of each timestep, if they are available. 0=off, 1=warn (stderr) on failure, 2=exit on failure (0)",
+ ARG_INTEGER ,
+ WTTS_USE_DEFAULT,
+ Var(stardata->preferences->float_overflow_checks),
+ 1.0
+ CMD_LINE_ARG_T_OPTIONAL_ARGS },
{
ARG_SECTION_OUTPUT,
"save_pre_events_stardata",
diff --git a/src/setup/set_cmd_line_macro_pairs.c b/src/setup/set_cmd_line_macro_pairs.c
index 77ee16c1cb80b43d1a09ad1e9ebbe81820935697..d48794b9435b8339a44ae93c4f3c7cf0892c88fd 100644
--- a/src/setup/set_cmd_line_macro_pairs.c
+++ b/src/setup/set_cmd_line_macro_pairs.c
@@ -383,10 +383,10 @@ void set_cmd_line_macro_pairs(struct stardata_t * RESTRICT const stardata,
#ifdef NUCSYN
-#define X(CODE,ONOFF,STRING,FUNCTION,LSODA,CVODE) Add_argument_Xmacro_pair("nucsyn_network%d",CODE,NUCSYN_NETWORK_##CODE);
+#define X(CODE,ONOFF,STRING,FUNCTION,LSODA,CVODE,CHECKFUNC) Add_argument_Xmacro_pair("nucsyn_network%d",CODE,NUCSYN_NETWORK_##CODE);
NUCSYN_NETWORKS_LIST;
#undef X
-#define X(CODE,ONOFF,STRING,FUNCTION,LSODA,CVODE) Add_argument_Xmacro_pair("nucsyn_network_error%d",CODE,NUCSYN_NETWORK_##CODE);
+#define X(CODE,ONOFF,STRING,FUNCTION,LSODA,CVODE,CHECKFUNC) Add_argument_Xmacro_pair("nucsyn_network_error%d",CODE,NUCSYN_NETWORK_##CODE);
NUCSYN_NETWORKS_LIST;
#undef X
diff --git a/src/setup/set_default_preferences.c b/src/setup/set_default_preferences.c
index 4fd7e5ab3ee41760da4fd46eab745d3540dd1d7e..0d2396cde378621fd500f163bfad4d51e7aab0d2 100644
--- a/src/setup/set_default_preferences.c
+++ b/src/setup/set_default_preferences.c
@@ -85,7 +85,7 @@ void set_default_preferences(struct preferences_t * RESTRICT const preferences)
preferences->nucsyn_network_error[i] = NUCSYN_NETWORK_DEFAULT_ERROR;
}
#undef X
-#define X(CODE,ONOFF,STRING,FUNCTION,LSODA,CVODE) \
+#define X(CODE,ONOFF,STRING,FUNCTION,LSODA,CVODE,CHECKFUNC) \
if(NUCSYN_NETWORK_##CODE!=NUCSYN_NETWORK_NUMBER) \
{ \
preferences->nucsyn_network[NUCSYN_NETWORK_##CODE] = (ONOFF); \
diff --git a/src/setup/version.c b/src/setup/version.c
index 2cabc7309db586704f295fbcde2658a6257e4011..fce970d520e3640df2be5d231e150fad65dc1ade 100644
--- a/src/setup/version.c
+++ b/src/setup/version.c
@@ -1745,14 +1745,16 @@ void version(struct stardata_t * RESTRICT const stardata)
{
int nnet = 0;
#undef X
-#define X(CODE,ONOFF,STRING,KAPS_RENTROP,LSODA,CVODE) \
- Printf("Network %d==%s is %s : implemented in %s%s%s\n", \
+#define X(CODE,ONOFF,STRING,KAPS_RENTROP,LSODA,CVODE,CHECKFUNC) \
+ Printf("Network %d==%s is %s : implemented in %s%s%s (checkfunc %s)\n", \
nnet++, \
#CODE, \
STRING, \
KAPS_RENTROP!=NULL ? "Kaps-Rentrop " : "", \
LSODA!=NULL ? "LSODA " : "", \
- CVODE!=NULL ? "CVODE " : "");
+ CVODE!=NULL ? "CVODE " : "", \
+ CHECKFUNC!=NULL ? #CHECKFUNC : "none" \
+ );
NUCSYN_NETWORKS_LIST;
#undef X
}
diff --git a/src/supernovae/supernova.c b/src/supernovae/supernova.c
index 5d03a44138ed38eb219833a95921e934e826590f..47fd93bb39ac45fc6bd0d98476e695391abe4228 100644
--- a/src/supernovae/supernova.c
+++ b/src/supernovae/supernova.c
@@ -35,9 +35,9 @@ void supernova(struct stardata_t * const stardata,
);
- struct star_t * pre_explosion_star = New_star_from(star);
- struct star_t * companion MAYBE_UNUSED = Other_star_struct(star);
- double mass_ejected MAYBE_UNUSED =
+ struct star_t * const pre_explosion_star = New_star_from(star);
+ struct star_t * const companion MAYBE_UNUSED = Other_star_struct(star);
+ const double mass_ejected MAYBE_UNUSED =
star->mass - star->new_supernova->new_stellar_structure->mass;
#ifdef LOG_SUPERNOVAE
@@ -119,7 +119,7 @@ void supernova(struct stardata_t * const stardata,
* Assume it has the same specific angular momentum
* as the pre-explosion stellar core.
*/
- double l = pre_explosion_star->angular_momentum/
+ const double l = pre_explosion_star->angular_momentum/
pre_explosion_star->mass;
@@ -231,10 +231,9 @@ void supernova(struct stardata_t * const stardata,
#if defined DISCS && defined DISC_LOG_POPSYN
Dprint("Logging (popsyn)\n");
- int i;
- for(i=0; i<stardata->common.ndiscs; i++)
+ for(int i=0; i<stardata->common.ndiscs; i++)
{
- struct disc_t * disc = & stardata->common.discs[i];
+ struct disc_t * const disc = & stardata->common.discs[i];
struct binary_system_t binary;
disc_init_binary_structure(stardata,&binary,disc);
if(disc->n_thermal_zones > 0 &&
diff --git a/src/timestep/timestep_RLOF.c b/src/timestep/timestep_RLOF.c
index 6da1a046c577e9983fc527c3da714bf68fbc8a8b..a433c60d27c43b1f075dea743fd2e06fce8a191a 100644
--- a/src/timestep/timestep_RLOF.c
+++ b/src/timestep/timestep_RLOF.c
@@ -14,7 +14,7 @@ void timestep_RLOF(Timestep_prototype_args)
if(stardata->model.sgl==FALSE)
{
/* binary star conditions */
- const double fr=star->radius/star->roche_radius;
+ const double fr = star->radius / star->roche_radius;
if(fr > RLOF_REDUCE_TIMESTEP_THRESHOLD)
{
diff --git a/src/timestep/timestep_limits.c b/src/timestep/timestep_limits.c
index 637b9ab35b7d9fbc0958debc1e028203f3fe5284..3fa81b4e54ff627d8f8d891238271897487fb717 100644
--- a/src/timestep/timestep_limits.c
+++ b/src/timestep/timestep_limits.c
@@ -192,6 +192,7 @@ void timestep_limits(Timestep_prototype_args)
fabs(star->derivative[DERIVATIVE_STELLAR_ANGULAR_VELOCITY_TIDES] /
Second_derivative(DERIVATIVE_STELLAR_ANGULAR_VELOCITY_TIDES));
Limit_timestep(*dt,dt_tides2,star,DT_LIMIT_TIDES2);
+ Dprint("DT second derivative\n");
}
#endif
diff --git a/tbse b/tbse
index d9afcb9155b29c1aedb7acaca8853f607a7fca47..5c80adb4a8d2981f1da8000a14fbee006babc0b8 100755
--- a/tbse
+++ b/tbse
@@ -1867,6 +1867,13 @@ TIMESTEP_LOGGING=F
# van den Heuvel diagram logging
VANDENHEUVEL_LOGGING=F
+# floating-point overflow checks at the
+# end of each timestep
+# 0 = no checks
+# 1 = warn on check failure
+# 2 = exit on check failure
+FLOAT_OVERFLOW_CHECKS=0
+
############################################################
# Orbiting objects
#
@@ -2586,6 +2593,7 @@ $ORBITING_OBJECTS \
--HeWD_HeWD_ignition_mass $HEWD_HEWD_IGNITION_MASS \
--wind_Nieuwenhuijzen_luminosity_lower_limit $WIND_NIEUWENHUIJZEN_LUMINOSITY_LOWER_LIMT \
--wind_LBV_luminosity_lower_limit $WIND_LBV_LUMINOSITY_LOWER_LIMT \
+--float_overflow_checks $FLOAT_OVERFLOW_CHECKS \
";
# valgrind options