From 521f0f2a80a92eb1d98922e6bffb3f98c2012077 Mon Sep 17 00:00:00 2001
From: Robert Izzard <r.izzard@surrey.ac.uk>
Date: Fri, 25 Oct 2019 17:33:19 +0100
Subject: [PATCH] more searching for the bug: narrowing it down
---
src/binary_c_maths.h | 8 +
src/disc/disc_bisection_rooter.c | 60 ++++++-
src/disc/disc_calc_disc_structure.c | 38 +++-
src/disc/disc_constraints.h | 22 ++-
src/disc/disc_logging.c | 2 +-
src/disc/disc_monotonic_check.c | 270 ++++++++++++++++++++++++++++
src/disc/disc_new_zone_list.c | 2 +-
src/disc/disc_new_zone_radii.c | 13 +-
src/disc/disc_nsector.c | 18 +-
src/disc/disc_parameters.h | 6 +-
src/disc/disc_prototypes.h | 7 +
src/disc/disc_show_thermal_zones.c | 130 +++++++-------
src/maths/generic_bisect.c | 2 +-
tbse | 4 +-
14 files changed, 495 insertions(+), 87 deletions(-)
create mode 100644 src/disc/disc_monotonic_check.c
diff --git a/src/binary_c_maths.h b/src/binary_c_maths.h
index 73db25a55..af7ab3bfb 100644
--- a/src/binary_c_maths.h
+++ b/src/binary_c_maths.h
@@ -324,6 +324,14 @@ static inline int func_is_really_zero(const double a){return(Is_really_zero_macr
#define Sign(X) (Is_zero(X) ? +1.0 : (X)/fabs(X))
#define Sign_is_really(X) (Is_really_zero(X) ? +1.0 : (X)/fabs(X))
+
+/*
+ * Sign macros, but give 0 when X is zero
+ */
+#define Sign_int(X) ((int)(Is_zero(X) ? 0 : ((X)/fabs(X) > 0.0 ? +1 : -1)))
+#define Sign_is_really_int(X) ((int)(Is_really_zero(X) ? 0 : ((X)/fabs(X) > 0.0 ? +1 : -1)))
+
+
/* range macros */
#define In_range(A,B,C) ((More_or_equal((A),(B)))&&(Less_or_equal((A),(C))))
#define Limit_range(A,B,C) (Max((B),Min((C),(A))))
diff --git a/src/disc/disc_bisection_rooter.c b/src/disc/disc_bisection_rooter.c
index bbe00d891..406190ad7 100644
--- a/src/disc/disc_bisection_rooter.c
+++ b/src/disc/disc_bisection_rooter.c
@@ -32,7 +32,7 @@ int disc_bisection_rooter(const int n,
int error;
unsigned int i,j;
- printf("BISECTOR %g\n",disc->lifetime);
+ printf("BISECTOR at t = %g y\n",disc->lifetime / YEAR_LENGTH_IN_SECONDS);
while(converged == FALSE && failed == FALSE)
{
@@ -42,14 +42,16 @@ int disc_bisection_rooter(const int n,
i = bisection_count % n;
Discdebug(1,
- "BISECT %lu param %u/%d %u %s : constraint %u %s\n",
+ "BISECT %lu param %u/%d %u %s : constraint %u %s = %g\n",
bisection_count,
i,
n,
parameters[i],
Disc_parameter_string(parameters[i]),
constraints[i],
- Disc_constraint_string(constraints[i]));
+ Disc_constraint_string(constraints[i]),
+ Disc_constraint_value(disc,binary,constraints[i])/Disc_constraint_unit(constraints[i])
+ );
double unknown MAYBE_UNUSED = 0.0;
/*
@@ -95,7 +97,7 @@ int disc_bisection_rooter(const int n,
local_parameters[0]),
DISC_BISECTION_TOLERANCE,
DISC_BISECTION_MAX_ITERATIONS,
- DISC_BISECT_BISECTION_ROOTER_USELOG,
+ FALSE,//DISC_BISECT_BISECTION_ROOTER_USELOG,
1.0,
/* va_args */
@@ -115,6 +117,42 @@ int disc_bisection_rooter(const int n,
"Bisection gave error %d \"%s\"\n",
error,
Bisect_error_string(error));
+
+ if(disc->lifetime > 1000.0 * YEAR_LENGTH_IN_SECONDS)
+ {
+ printf("Want disc->M = %g have %g\n",
+ disc->M,
+ disc_total_mass(disc));
+
+ double Tvwas = disc->Tvisc0;
+ double dTv = 0.001 * disc->Tvisc0;
+ struct disc_t * discwas = New_disc_from(disc);
+ double Tv;
+ double residual[n];
+ int failure;
+ for(Tv = 0.9 * Tvwas;
+ Tv <= 1.1 * Tvwas;
+ Tv += dTv)
+ {
+ Copy_disc(discwas,disc);
+ disc->Tvisc0 = Tv;
+ disc_nsector2(stardata,
+ disc,
+ binary,
+ n,
+ constraints,
+ residual,
+ &failure);
+ printf("CHECKNEAR %30.15e %30.15g eps %g\n",
+ disc->Tvisc0,
+ disc_total_mass(disc),
+ (disc_total_mass(disc) - disc->M)/disc->M
+ );
+
+ }
+ disc_monotonic_check(stardata,disc,binary,n);
+ Flexit;
+ }
if(error != BINARY_C_BISECT_ERROR_NONE)
{
failed = i + 1;
@@ -216,14 +254,14 @@ int disc_bisection_rooter(const int n,
);
- if(bisection_count > 1000)
+ if(bisection_count > 1200)
{
/*
* We're in trouble
*/
- printf("Exit on long bisection\n");
+ //printf("Exit on long bisection\n");
//disc_parameter_space(stardata,disc);
- Flexit;
+ //Flexit;
}
bisection_count++;
@@ -306,6 +344,14 @@ static double disc_nsector_wrapper(const double x,
failed
);
+ printf("disc mass %g = %g Msun : want %g = %g Msun -> residual %g\n",
+ disc_total_mass(disc),
+ disc_total_mass(disc)/M_SUN,
+ disc->M,
+ disc->M/M_SUN,
+ residual[0]
+ );
+
((struct GSL_args_t *)p)->error = *failed;
Discdebug(2,
diff --git a/src/disc/disc_calc_disc_structure.c b/src/disc/disc_calc_disc_structure.c
index 00d7f9d36..c54a20edd 100644
--- a/src/disc/disc_calc_disc_structure.c
+++ b/src/disc/disc_calc_disc_structure.c
@@ -364,6 +364,8 @@ static double converge_with_torque_wrapper(const double torquef,
disc->torqueF = torquef;
discwas->torqueF = torquef;
+
+
double x = converge_with_torque(n,
parameter_types,
constraints,
@@ -996,6 +998,11 @@ static int disc_GSL_multiroot(const int solver,
"Failed to allocated memory for GSL multiroot fsolver\n");
}
#endif//ALLOC_CHECKS
+
+ /*
+ * use logs?
+ */
+#define USELOG
/*
* Set up the test function
@@ -1022,10 +1029,15 @@ static int disc_GSL_multiroot(const int solver,
double unknown = 0.0;
for(i=0;i<n;i++)
{
+ const double x = (double)Disc_parameter(disc,parameter_types[i]);
gsl_vector_set(
guess_vector,
(size_t)i,
- (double)Disc_parameter(disc,parameter_types[i])
+#ifdef USELOG
+ log(x)
+#else
+ x
+#endif//USELOG
);
}
@@ -1104,6 +1116,13 @@ static int disc_GSL_multiroot(const int solver,
status = gsl_multiroot_fsolver_iterate(s);
+ /*
+ disc_monotonic_check(stardata,
+ disc,
+ binary,
+ n);
+ */
+
/*
* If the guess is NaN then something
* has gone very wrong and we should bail as
@@ -1130,6 +1149,17 @@ static int disc_GSL_multiroot(const int solver,
status,
gsl_strerror(status),
iter);
+
+ /* debugging */
+ if(status==27 && 0)
+ {
+ printf("no progress error\n");
+ disc_monotonic_check(stardata,
+ disc,
+ binary,
+ n);
+ Flexit;
+ }
break;
}
@@ -1214,6 +1244,12 @@ static int GSL_multiroot_testfunc(const gsl_vector * x,
{
parameter_values[i] =
gsl_vector_get(x,i);
+#ifdef USELOG
+ parameter_values[i] = exp(parameter_values[i]);
+#endif
+ printf("parameter %d -> %g\n",
+ i,
+ parameter_values[i]);
}
double residuals[n];
diff --git a/src/disc/disc_constraints.h b/src/disc/disc_constraints.h
index 94a02648d..5c52d3be2 100644
--- a/src/disc/disc_constraints.h
+++ b/src/disc/disc_constraints.h
@@ -70,6 +70,16 @@
0.0 \
)
+/*
+ * Unit value for constraint value
+ */
+#define Disc_constraint_unit(TYPE) \
+ ( \
+ (TYPE) == DISC_CONSTRAINT_M_TO_DISC_M ? M_SUN : \
+ ((TYPE) == DISC_CONSTRAINT_RIN_MIN_TO_DISC_RIN || \
+ (TYPE) == DISC_CONSTRAINT_ROUT_MAX_TO_DISC_ROUT) ? R_SUN : \
+ 1.0 \
+ )
/*
* Macro to construct a residual and
@@ -190,7 +200,7 @@
(TYPE) == DISC_PARAMETER_TVISC0 ? "Tvisc0" : \
(TYPE) == DISC_PARAMETER_RIN ? "Rin" : \
(TYPE) == DISC_PARAMETER_ROUT ? "Rout" : \
- (TYPE) == DISC_PARAMETER_ROUT ? "Torquef" : \
+ (TYPE) == DISC_PARAMETER_TORQUEF ? "Torquef" : \
"Unknown")
@@ -240,7 +250,7 @@
#define _EPS 1e-6
#define _JUST_MORE (1.0 + _EPS)
#define _JUST_LESS (1.0 - _EPS)
-
+/*
#define Disc_parameter_min(D,TYPE) \
( \
(TYPE) == DISC_PARAMETER_TVISC0 ? DISC_TVISC0_MIN : \
@@ -249,6 +259,14 @@
(D)->Rin*_JUST_MORE)): \
1e-100 \
)
+*/
+#define Disc_parameter_min(D,TYPE) \
+ ( \
+ (TYPE) == DISC_PARAMETER_TVISC0 ? DISC_TVISC0_MIN : \
+ (TYPE) == DISC_PARAMETER_RIN ? DISC_RIN_MIN : \
+ (TYPE) == DISC_PARAMETER_ROUT ? DISC_RIN_MIN : \
+ 1e-100 \
+ )
/*
diff --git a/src/disc/disc_logging.c b/src/disc/disc_logging.c
index 5be91e1c0..005b26aa1 100644
--- a/src/disc/disc_logging.c
+++ b/src/disc/disc_logging.c
@@ -303,7 +303,7 @@ void disc_logging(struct stardata_t * const stardata,
/*
* File log for splots
*/
- if(Disc_can_log2d &&
+ if(Disc_can_log2d == TRUE &&
Disc_is_disc(disc) &&
Disc_log_filepointer2d==NULL)
{
diff --git a/src/disc/disc_monotonic_check.c b/src/disc/disc_monotonic_check.c
new file mode 100644
index 000000000..a6a855efa
--- /dev/null
+++ b/src/disc/disc_monotonic_check.c
@@ -0,0 +1,270 @@
+#include "../binary_c.h"
+
+#if defined DISCS
+
+Boolean disc_monotonic_check(struct stardata_t * const stardata,
+ struct disc_t * const disc,
+ struct binary_system_t * const binary,
+ const int n
+ )
+{
+ /*
+ * Check that the disc M,J,F change monotonically in
+ * the disc parameters
+ */
+ double unknown MAYBE_UNUSED;
+ unsigned int i;
+
+ /*
+ * Allocate memory
+ */
+ double *ref = Malloc(sizeof(double)*n);
+ disc_parameter * parameters = Malloc(sizeof(disc_parameter)*n);
+ disc_constraint * constraints = Malloc(sizeof(disc_constraint)*n);
+
+ /*
+ * Set the parameters that will be changed
+ */
+ parameters[0] = DISC_PARAMETER_TVISC0;
+ parameters[1] = DISC_PARAMETER_RIN;
+ parameters[2] = DISC_PARAMETER_ROUT;
+
+ /*
+ * Calculate pre-pertubation reference parameters
+ */
+ for(i=0;i<n;i++)
+ {
+ ref[i] = Disc_parameter(disc,parameters[i]);
+ }
+
+ /*
+ * Set up constraints;
+ */
+ constraints[0] = DISC_CONSTRAINT_M_TO_DISC_M;
+ constraints[1] = DISC_CONSTRAINT_J_TO_DISC_J;
+ constraints[2] = DISC_CONSTRAINT_F_TO_DISC_F;
+
+ /*
+ * Backup disc
+ */
+ struct disc_t * discwas = New_disc_from(disc);
+
+ /*
+ * Vary each in turn
+ */
+ const double eps = 1e-2;
+ int j;
+ for(i=0;i<n;i++)
+ {
+ /*
+ * Perturb parameter i by -eps and +eps
+ */
+ const double d = eps * ref[i];
+ double derivative = 0.0;
+ for(j=-1; j<=+1; j+=2)
+ {
+ /*
+ * Restore disc
+ */
+ Copy_disc(discwas,disc);
+
+ /*
+ * Set parameters in the disc
+ */
+ Disc_parameter(disc,parameters[i]) = ref[i] + (double)j * d;
+
+ /*
+ * Calculate the new disc structure
+ */
+ double residual[n];
+ int failure;
+ disc_nsector2(stardata,
+ disc,
+ binary,
+ n,
+ constraints,
+ residual,
+ &failure);
+
+ derivative =
+ j == -1 ? Disc_constraint_value(disc,binary,constraints[i]) :
+ j == +1 ? ((Disc_constraint_value(disc,binary,constraints[i]) - derivative)/(2.0*d)) :
+ 0.0;
+
+ printf("PERTURB %u : %d : %30.15e : ",
+ i,j,Disc_parameter(disc,parameters[i]));
+ printf("Residuals %g %g %g",
+ residual[0],
+ residual[1],
+ residual[2]
+ );
+ if(j==+1)
+ {
+ printf(", derivative %u = %g",i,derivative);
+ }
+ printf("\n");
+ }
+ }
+
+
+#define __FULLRANGE
+#define __LOG
+
+ Boolean debug = FALSE;
+
+#ifdef __FULLRANGE
+ /*
+ * check parameters over the full range
+ */
+ const int N = 100.0;
+ double d[n],prev[n];
+ for(i=0;i<n;i++)
+ {
+ d[i] = 0.0;
+ }
+
+ for(i=0;i<n;i++)
+ {
+ unsigned int j;
+ int sign[n];
+ for(j=0;j<n;j++)
+ {
+ sign[j] = 0;
+ }
+
+ const double logxmin = log10(Disc_parameter_min(disc,parameters[i]));
+ const double logxmax = log10(Disc_parameter_max(disc,parameters[i]));
+ const double dlogx = (logxmax - logxmin)/(1.0*(double)N);
+ double logx = logxmin;
+ int count = 0;
+ while(logx < logxmax)
+ {
+ double x = Pow10(logx);
+ if(debug)printf("Set var %u = %g\n",i,x);
+ /*
+ * Restore disc
+ */
+ Copy_disc(discwas,disc);
+
+ /*
+ * Set parameter i in the disc
+ */
+ Disc_parameter(disc,parameters[i]) = x;
+
+ /*
+ * Calculate the new disc structure
+ */
+ double residual[n];
+ int failure;
+ disc_nsector2(stardata,
+ disc,
+ binary,
+ n,
+ constraints,
+ residual,
+ &failure);
+
+ double c[n];
+ for(j=0;j<n;j++)
+ {
+ c[j] = Disc_constraint_value(disc,binary,constraints[j]);
+ if(debug)printf("Set constraint %u = %30.15e\n",j,c[j]);
+ }
+
+ /*
+ * If we have prev, use it
+ */
+ if(count>0)
+ {
+ for(j=0;j<n;j++)
+ {
+ /* derivative (to within a positive constant) */
+ d[j] = Abs_diff(c[j],prev[j]);
+
+ /* numerical error handling */
+ if(fabs(d[j]) < 1e-8)
+ {
+ d[j]=0.0;
+ }
+ if(debug)printf("Set derivative %u = %g (zero? %s)\n",
+ j,
+ d[j],
+ Yesno(Is_really_zero(d[j])));
+ }
+
+ if(count == 1)
+ {
+ /* first time : save the sign of the derivative */
+ for(j=0;j<n;j++)
+ {
+ sign[j] = Sign_is_really_int(d[j]);
+ if(debug)printf("Set sign (first time) %u = %d\n",
+ j,
+ sign[j]);
+ }
+ }
+ else if(count>=2)
+ {
+ /* subsequent times : check the sign is the same */
+ int newsign[n];
+ for(j=0;j<n;j++)
+ {
+ newsign[j] = Sign_is_really_int(d[j]);
+ if(debug)printf("Set newsign %u = %d\n",j,newsign[j]);
+
+ if(sign[j] != 0 &&
+ newsign[j] != 0 &&
+ newsign[j] != sign[j])
+ {
+ Exit_binary_c(BINARY_C_BISECT_ERROR,
+ "%u is Non-monotonic (sign now %d, was %d)\n",
+ j,
+ Sign_is_really_int(d[j]),
+ sign[j]
+ );
+ }
+ }
+ if(debug)
+ {
+ printf("SIGNS ");
+ for(j=0;j<n;j++)
+ {
+ printf("%d ",Sign_is_really_int(d[j]));
+ }
+ printf("\n");
+ }
+ }
+ }
+
+#ifdef __LOG
+ printf("FULLRANGE %u %30.15e %30.15e %30.15e %30.15e %g %g %g\n",
+ i,
+ x,
+ c[0],
+ c[1],
+ c[2],
+ fabs(residual[0]),
+ fabs(residual[1]),
+ fabs(residual[2])
+ );
+#endif
+ logx += dlogx;
+
+ count++;
+ for(j=0;j<n;j++)
+ {
+ prev[j] = c[j];
+ }
+ }
+ }
+#endif // __FULLRANGE
+
+ /*
+ * Restore the disc
+ */
+ Copy_disc(discwas,disc);
+
+ return TRUE;
+}
+
+#endif // DISCS
diff --git a/src/disc/disc_new_zone_list.c b/src/disc/disc_new_zone_list.c
index 59f209dd0..129b083dd 100644
--- a/src/disc/disc_new_zone_list.c
+++ b/src/disc/disc_new_zone_list.c
@@ -24,7 +24,7 @@ Boolean disc_new_zone_list(struct disc_t * const disc,
/*
* Fail if input parameters are unphysical, e.g. Tvisc0 < 0
*/
- if(unlikely(disc->Tvisc0 < TINY))
+ if(unlikely(disc->Tvisc0 < REALLY_TINY))
{
return FALSE;
}
diff --git a/src/disc/disc_new_zone_radii.c b/src/disc/disc_new_zone_radii.c
index b0e2a1093..21398d80a 100644
--- a/src/disc/disc_new_zone_radii.c
+++ b/src/disc/disc_new_zone_radii.c
@@ -100,6 +100,17 @@ Boolean disc_new_zone_radii(struct disc_t * const disc)
Zprintf("Hottest %u\n",ihottest);
+ /*
+ * If ihottest is -1 we're in trouble
+ */
+ if(ihottest == -1)
+ {
+ Safe_free(P);
+ Safe_free(Q);
+ printf("ihottest -1 error\n");Flexit;
+ return FALSE;
+ }
+
/*
* Zone ihottest, in P, is the hottest at r
*
@@ -249,7 +260,7 @@ Boolean disc_new_zone_radii(struct disc_t * const disc)
disc_set_temperature_power_laws(disc);
Zprintf("\n\n\n");
- disc_show_thermal_zones(disc,disc->thermal_zones,DISCS_MAX_N_ZONES);
+ //disc_show_thermal_zones(disc,disc->thermal_zones,DISCS_MAX_N_ZONES);
Safe_free(P);
Safe_free(Q);
diff --git a/src/disc/disc_nsector.c b/src/disc/disc_nsector.c
index 2bd917208..c036b1b63 100644
--- a/src/disc/disc_nsector.c
+++ b/src/disc/disc_nsector.c
@@ -29,7 +29,8 @@
*
* The constraints and parameter_types are defined in disc_constraints.h
*/
-
+//#undef DISC_DEBUG
+//#define DISC_DEBUG 2
void disc_nsector(const unsigned int n,
const double * parameter_value,
double * residual,
@@ -50,6 +51,11 @@ void disc_nsector(const unsigned int n,
double unknown MAYBE_UNUSED;
struct disc_t * discwas;
unsigned int i,j;
+
+ printf("conv in disc_nsector Tvisc0 = %g Rin = %g Rout = %g\n",
+ disc->Tvisc0,
+ disc->Rin,
+ disc->Rout);
/*
* Back up disc
@@ -62,7 +68,7 @@ void disc_nsector(const unsigned int n,
for(i=0;i<n;i++)
{
Discdebug(2,
- "Nsector parameter %3u/%3u type %3u %40s value = %30.20g\n",
+ "Nsector parameter %3u/%3u type %3u %40s value = %30.15e\n",
i,
n,
parameter_types[i],
@@ -199,7 +205,7 @@ void disc_nsector2(struct stardata_t * stardata,
1.0 - disc->Rin/disc->Rout
);
Discdebug(2,
- "Nsector: disc->M = %g Msun cf M = %g Msun, disc->J = %g cf J = %g, disc->F = %g cf F = %g, residual = %g\n",
+ "Nsector: disc->M = %g Msun cf M = %g Msun, disc->J = %g cf J = %g, disc->F = %g cf F = %g, residual = %30.15e\n",
disc->M/M_SUN,
disc_total_mass(disc)/M_SUN,
disc->J,
@@ -209,7 +215,7 @@ void disc_nsector2(struct stardata_t * stardata,
Disc_residual(disc,binary,DISC_CONSTRAINT_M_TO_DISC_M)
);
Discdebug(2,
- "Nsector: (disc->J = %g) / (J = %g) - 1.0 = %g, residual = %g\n",
+ "Nsector: (disc->J = %g) / (J = %g) - 1.0 = %g, residual = %30.15e\n",
disc->J,
disc_total_angular_momentum(disc,binary),
disc->J/disc_total_angular_momentum(disc,binary)-1.0,
@@ -218,7 +224,7 @@ void disc_nsector2(struct stardata_t * stardata,
double F = disc_total_angular_momentum_flux_from_binary(disc,
binary);
Discdebug(2,
- "Nsector: (disc->F = %g) / (F = %g) - 1.0 = %g\n",
+ "Nsector: (disc->F = %g) / (F = %g) - 1.0 = %30.15e\n",
disc->F,
F,
disc->F / F - 1.0
@@ -235,7 +241,7 @@ void disc_nsector2(struct stardata_t * stardata,
const int constraint = constraints[i];
residual[i] = Disc_residual(disc,binary,constraint);
Discdebug(2,
- "Nsector: constraint %3d type %3d %40s residual = %30.20g\n",
+ "Nsector: constraint %3d type %3d %40s residual = %30.15g\n",
i,
constraint,
Disc_constraint_string(constraint),
diff --git a/src/disc/disc_parameters.h b/src/disc/disc_parameters.h
index 1db0b4577..6a4f07f04 100644
--- a/src/disc/disc_parameters.h
+++ b/src/disc/disc_parameters.h
@@ -184,7 +184,7 @@
* and these work fine). Defaults: 1e-6 and 1000.
*/
-#define DISC_TOLERANCE 1e-6
+#define DISC_TOLERANCE 1e-4
#define DISC_BISECTION_MAX_ITERATIONS 1000
/*
@@ -278,7 +278,7 @@
/*
* Maximum number of iterations through the GSL root finder
*/
-#define DISC_ROOT_FINDER_MAX_ATTEMPTS 100
+#define DISC_ROOT_FINDER_MAX_ATTEMPTS 1000
#define DISC_SMALL_RADIUS (1e2)
#define DISC_LARGE_RADIUS (1e50)
@@ -390,7 +390,7 @@
* Usually this is very slightly slower, but may help explore
* larger parameter spaces faster.
*/
-#define DISC_BISECT_USELOG FALSE
+#define DISC_BISECT_USELOG TRUE
#define DISC_BISECT_OWEN_RADIUS_USELOG DISC_BISECT_USELOG
#define DISC_BISECT_RJ_USELOG DISC_BISECT_USELOG
diff --git a/src/disc/disc_prototypes.h b/src/disc/disc_prototypes.h
index 580e72a8d..15193660c 100644
--- a/src/disc/disc_prototypes.h
+++ b/src/disc/disc_prototypes.h
@@ -487,6 +487,13 @@ void disc_zone_crossing_radii(struct disc_t * const disc,
Boolean disc_new_zone_radii(struct disc_t * const disc);
+
+Boolean disc_monotonic_check(struct stardata_t * const stardata,
+ struct disc_t * const disc,
+ struct binary_system_t * const binary,
+ const int n
+ );
+
#endif // DISCS
#endif // DISC_PROTOTYPES_H
diff --git a/src/disc/disc_show_thermal_zones.c b/src/disc/disc_show_thermal_zones.c
index a4a2062b6..0480de6cb 100644
--- a/src/disc/disc_show_thermal_zones.c
+++ b/src/disc/disc_show_thermal_zones.c
@@ -9,15 +9,17 @@ void disc_show_thermal_zones(struct disc_t * disc,
const int n)
{
const int max = n==-1 ? DISCS_MAX_N_ZONES : n;
- int i,show=0;
+ int i;
+ const Boolean show_invalid = FALSE;
+ Boolean show = FALSE;
for(i=0;i<max;i++)
{
if((zones+i)->valid==TRUE)
{
- show++;
+ show = TRUE;
+ break;
}
}
- show = (show>0) ? 1 : 0;
#define NORMAL_COLOUR(z) (z->valid==TRUE ? BRIGHT_GREEN : RED)
#define ERROR_COLOUR (YELLOW)
@@ -25,73 +27,77 @@ void disc_show_thermal_zones(struct disc_t * disc,
for(i=0;i<max;i++)
{
- const struct disc_thermal_zone_t * z = zones + i;
- printf("Zone %d at %p :",i,z);
- fflush(stdout);
- printf("type %s : %s %s: ",
- DISC_ZONE_TYPE_STRING(z),
- DISC_ZONE_VALID_STRING(z),
- NORMAL_COLOUR(z)
- );
-
- fflush(stdout);
- if(show)
+ if(show_invalid == TRUE || (zones+i)->valid == TRUE)
{
- if(0){
- Boolean xxx=(i!=max-1 && (zones+i+1)->valid==TRUE);
- if(xxx==TRUE)
- {
- printf("R transition=%5.3e ",i!=max-1 ? (disc_Rcross(z+i,z+i+1)/R_SUN) : 0);
- fflush(stdout);
- }
- else
- {
- printf(" ");
+
+ const struct disc_thermal_zone_t * z = zones + i;
+ printf("Zone %d at %p :",i,z);
+ fflush(stdout);
+ printf("type %s : %s %s: ",
+ DISC_ZONE_TYPE_STRING(z),
+ DISC_ZONE_VALID_STRING(z),
+ NORMAL_COLOUR(z)
+ );
+
+ fflush(stdout);
+ if(show == TRUE)
+ {
+ if(0){
+ Boolean xxx=(i!=max-1 && (zones+i+1)->valid==TRUE);
+ if(xxx==TRUE)
+ {
+ printf("R transition=%5.3e ",i!=max-1 ? (disc_Rcross(z+i,z+i+1)/R_SUN) : 0);
+ fflush(stdout);
+ }
+ else
+ {
+ printf(" ");
+ }
}
+ printf(" %sstart=%g%s%s, %send=%g%s%s, %swidth=%g%s%s ",
+ ASSERT_COLOUR(!(More_or_equal(z->rstart/disc->Rin,1.0) &&
+ Less_or_equal(z->rstart/disc->Rout,1.0))),
+ z->rstart/R_SUN,
+ R_SOLAR,
+ NORMAL_COLOUR(z),
+ ASSERT_COLOUR(!(More_or_equal(z->rend/disc->Rin,1.0) &&
+ Less_or_equal(z->rend/disc->Rout,1.0))),
+ z->rend/R_SUN,
+ R_SOLAR,
+ NORMAL_COLOUR(z),
+ ASSERT_COLOUR(z->rstart >= z->rend),
+ (z->rend - z->rstart)/R_SUN,
+ R_SOLAR,
+ NORMAL_COLOUR(z)
+ );
+ fflush(stdout);
+
}
- printf(" %sstart=%g%s%s, %send=%g%s%s, %swidth=%g%s%s ",
- ASSERT_COLOUR(!(More_or_equal(z->rstart/disc->Rin,1.0) &&
- Less_or_equal(z->rstart/disc->Rout,1.0))),
+
+ double dmzone = disc_partial_mass(disc,
+ z->rstart,
+ z->rend);
+ printf(": prefactor=%g exponent=%g%s Tstart=%g (from A0=%g R0=%g %s exponent=%g) T(rstart)=%g Sigma(rstart)=%g [zone M=%20.10g Msun] %s\n",
+ z->Tlaw.A0,
+ z->Tlaw.exponent,
+ CYAN,
+ (
+ Is_really_zero(z->Tlaw.A1) ?
+ // power law not yet complete
+ (z->Tlaw.A0*pow(z->rstart,z->Tlaw.exponent)) :
+ // power law is set up properly
+ power_law(&(z->Tlaw),z->rstart)),
+ z->Tlaw.A0,
z->rstart/R_SUN,
R_SOLAR,
- NORMAL_COLOUR(z),
- ASSERT_COLOUR(!(More_or_equal(z->rend/disc->Rin,1.0) &&
- Less_or_equal(z->rend/disc->Rout,1.0))),
- z->rend/R_SUN,
- R_SOLAR,
- NORMAL_COLOUR(z),
- ASSERT_COLOUR(z->rstart >= z->rend),
- (z->rend - z->rstart)/R_SUN,
- R_SOLAR,
- NORMAL_COLOUR(z)
+ z->Tlaw.exponent,
+ power_law(&z->Tlaw,z->rstart),
+ power_law(&z->power_laws[POWER_LAW_SIGMA],z->rstart),
+ dmzone/M_SUN,
+ COLOUR_RESET
);
fflush(stdout);
-
}
-
- double dmzone = disc_partial_mass(disc,
- z->rstart,
- z->rend);
- printf(": prefactor=%g exponent=%g%s Tstart=%g (from A0=%g R0=%g %s exponent=%g) T(rstart)=%g Sigma(rstart)=%g [zone M=%20.10g Msun] %s\n",
- z->Tlaw.A0,
- z->Tlaw.exponent,
- CYAN,
- (
- Is_really_zero(z->Tlaw.A1) ?
- // power law not yet complete
- (z->Tlaw.A0*pow(z->rstart,z->Tlaw.exponent)) :
- // power law is set up properly
- power_law(&(z->Tlaw),z->rstart)),
- z->Tlaw.A0,
- z->rstart/R_SUN,
- R_SOLAR,
- z->Tlaw.exponent,
- power_law(&z->Tlaw,z->rstart),
- power_law(&z->power_laws[POWER_LAW_SIGMA],z->rstart),
- dmzone/M_SUN,
- COLOUR_RESET
- );
- fflush(stdout);
}
printf("\n");
}
diff --git a/src/maths/generic_bisect.c b/src/maths/generic_bisect.c
index e8c8d2f1c..764301a69 100644
--- a/src/maths/generic_bisect.c
+++ b/src/maths/generic_bisect.c
@@ -90,7 +90,7 @@ double generic_bisect(int * RESTRICT error,
max = log10(Max(REALLY_TINY,max));
guess = log10(guess);
}
- int vb = 0;
+ int vb = 1;
*error = BINARY_C_BISECT_ERROR_NONE;
va_list args_master,args;
double top,bottom;
diff --git a/tbse b/tbse
index 8fe241df4..51f36d68a 100755
--- a/tbse
+++ b/tbse
@@ -578,7 +578,7 @@ DISC_TIMESTEP_FACTOR=1.0
DISC_LOG=DISC_LOG_LEVEL_SUBTIMESTEP_FIRST_DISC_ONLY
# soft 2D disc log switch, same options as DISC_LOG
-DISC_LOG2D=DISC_LOG_LEVEL_SUBTIMESTEP_FIRST_DISC_ONLY
+DISC_LOG2D=DISC_LOG_LEVEL_NONE
# disc log every DISC_LOG_DT Myr only. Ignored if zero.
# If set to -1, is automatically calculated (1Myr except on the GB/AGB)
@@ -2117,7 +2117,7 @@ $WARMUP \
--angelou_lithium_EAGB_massfrac $angelou_lithium_EAGB_massfrac \
--angelou_lithium_TPAGB_massfrac $angelou_lithium_TPAGB_massfrac \
--disc_log $DISC_LOG \
---disc_log2d $DISC_LOG \
+--disc_log2d $DISC_LOG2D \
--disc_log_dt $DISC_LOG_DT \
--disc_n_monte_carlo_guesses $DISC_N_MONTE_CARLO_GUESSES \
--type_Ia_MCh_supernova_algorithm $TYPE_IA_MCH_SUPERNOVA_ALGORITHM \
--
GitLab