From 63cfdd99e9378232836cd64fefcfd5d1ef3383b8 Mon Sep 17 00:00:00 2001
From: Robert Izzard <r.izzard@surrey.ac.uk>
Date: Mon, 13 May 2019 20:41:01 +0100
Subject: [PATCH] update solver test script
fix timesteps to correct for the wrong size of the arrays
---
src/binary_star_functions/tides.c | 12 +++++------
src/evolution/evolution.c | 6 ++++++
src/timestep/stellar_timestep.c | 2 +-
src/timestep/timestep.h | 6 +++---
src/timestep/timestep_limits.c | 10 +++++----
tbse | 2 +-
test_solvers.sh | 35 +++++++++++++++++++++++--------
7 files changed, 49 insertions(+), 24 deletions(-)
diff --git a/src/binary_star_functions/tides.c b/src/binary_star_functions/tides.c
index fd4169176..4691210b0 100644
--- a/src/binary_star_functions/tides.c
+++ b/src/binary_star_functions/tides.c
@@ -18,10 +18,10 @@ void tides(struct stardata_t *stardata,
double f2,f3,f4,f5,raa2,raa6,tc,f,tcqr,rg2,ttid;
double cq = companion->q;
double ecc2 = ECC_SQUARED;
- double omecc2 = 1.0 - ecc2;
- double sqome2 = sqrt(omecc2);
- double sqome3 = POW3(sqome2);
- double sq6 = POW6(sqome2);
+ double omecc2 = 1.0 - ecc2; // 1-e^2
+ double sqome2 = sqrt(omecc2); // sqrt(1-e^2)
+ double sqome3 = POW3(sqome2); // (1-e^2)^(3/2)
+ double sq6 = POW3(omecc2); // (1-e^2)^3
double isqome2_13 = 1.0/(sq6*sq6*sqome2);
double starI = moment_of_inertia(star,radius);
const int ndonor MAYBE_UNUSED = stardata->model.ndonor;
@@ -101,7 +101,7 @@ void tides(struct stardata_t *stardata,
(3.0*star->luminosity)),1.0/3.0);
ttid = TWOPI/(1.0e-10 + fabs(stardata->common.orbit.angular_frequency - star->omega));
-
+
/*
* This is correct! It is wrong in Jarrod's BSE code.
* DO NOT CHANGE IT (see MNRAS 329 897 eq 32)
@@ -174,7 +174,7 @@ void tides(struct stardata_t *stardata,
const double x = star->omega;
const double y = f2*stardata->common.orbit.angular_frequency/sqome3f5;
-#undef EXPONENTIAL_TIDES
+//#undef EXPONENTIAL_TIDES
#ifdef EXPONENTIAL_TIDES
/*
* Solve the differential equation analytically, and use the derivative that
diff --git a/src/evolution/evolution.c b/src/evolution/evolution.c
index c3e3080a5..a365ac274 100644
--- a/src/evolution/evolution.c
+++ b/src/evolution/evolution.c
@@ -251,7 +251,13 @@ int evolution(struct stardata_t * RESTRICT stardata,
stardata->model.dt
);
+ if(IS_ZERO(stardata->model.dt))
+ {
+ /* reset in case of repeat */
+ stardata->common.test = 0.0;
+ }
+
double absdiff = ABSDIFF(stardata->model.time*1e6,
stardata->common.test);
printf("POSTTEST %30.12e %30.12e : diff %g\n",
diff --git a/src/timestep/stellar_timestep.c b/src/timestep/stellar_timestep.c
index 753f9a538..c8a7a92ad 100644
--- a/src/timestep/stellar_timestep.c
+++ b/src/timestep/stellar_timestep.c
@@ -90,7 +90,7 @@ void stellar_timestep(Timestep_prototype_args)
star->stellar_timestep = *dt;
-//#define ___LOGGING
+#define ___LOGGING
#ifdef ___LOGGING
const char * string[DT_NUMBER_OF_TIMESTEP_LIMITERS] = { DT_LIMIT_STRINGS };
// if(star->starnum==0)
diff --git a/src/timestep/timestep.h b/src/timestep/timestep.h
index 605a4a3dc..3c0cf10b1 100644
--- a/src/timestep/timestep.h
+++ b/src/timestep/timestep.h
@@ -173,7 +173,7 @@ const Stellar_type stellar_type, \
#define DT_LIMIT_MASS_GAIN2 53
#define DT_LIMIT_MASS_LOSS2 54
-#define DT_NUMBER_OF_TIMESTEP_LIMITERS 53
+#define DT_NUMBER_OF_TIMESTEP_LIMITERS 55
#define DT_LIMIT_STRINGS \
"DT_LIMIT_NONE", \
@@ -226,10 +226,10 @@ const Stellar_type stellar_type, \
"DT_LIMIT_BURN_IN", \
"DT_LIMIT_RADIUS_CHANGES", \
"DT_LIMIT_MASSLESS_REMNANT", \
- "DT_LIMIT_ORBITAL_ANGMOM", \
+ "DT_LIMIT_ORBITAL_ANGMOM",/*50*/ \
"DT_LIMIT_STELLAR_ANGMOM2", \
"DT_LIMIT_STELLAR_TIDES2", \
"DT_LIMIT_MASS_GAIN2", \
- "DT_LIMIT_MASS_LOSS2"
+ "DT_LIMIT_MASS_LOSS2"/*54*/
#endif // TIMESTEP_H
diff --git a/src/timestep/timestep_limits.c b/src/timestep/timestep_limits.c
index 8a6dc3f22..4378737fb 100644
--- a/src/timestep/timestep_limits.c
+++ b/src/timestep/timestep_limits.c
@@ -6,7 +6,7 @@
#include "../nucsyn/nucsyn_angelou.h"
#endif
-#define SECOND_DERIVATIVES
+//#define SECOND_DERIVATIVES
#define SECOND_DERIVATIVE_F 0.1
void timestep_limits(Timestep_prototype_args)
@@ -17,13 +17,15 @@ void timestep_limits(Timestep_prototype_args)
const double fsolver =
stardata->preferences->solver == SOLVER_FORWARD_EULER ? 1.0 :
stardata->preferences->solver == SOLVER_RK2 ? 1.0 :
- stardata->preferences->solver == SOLVER_RK4 ? 10.0 :
+ stardata->preferences->solver == SOLVER_RK4 ? 1.0 :
stardata->preferences->solver == SOLVER_PREDICTOR_CORRECTOR ? 1.0 :
1.0;
-
+
+#ifdef SECOND_DERIVATIVES
struct star_t * prevstar =
&stardata->previous_stardata->star[star->starnum];
-
+#endif
+
#define Second_derivative(N) ( \
(star->derivative[(N)] \
- \
diff --git a/tbse b/tbse
index e66d56d60..e334008a0 100755
--- a/tbse
+++ b/tbse
@@ -302,7 +302,7 @@ GRAVITATIONAL_RADIATION_MODULATOR_E=1.0
# 1 = Andronov, Pinnsoneault and Sills 2003
# 2 = Barnes and Kim 2010
# 3 = Rappaport et al. 1983 (see magnetic_braking_gamma)
-MAGNETIC_BRAKING_ALGORITHM=1
+MAGNETIC_BRAKING_ALGORITHM=0
# Magnetic braking multipier
MAGNETIC_BRAKING_FACTOR=0.0
diff --git a/test_solvers.sh b/test_solvers.sh
index cff6d2767..c604dcd54 100755
--- a/test_solvers.sh
+++ b/test_solvers.sh
@@ -13,26 +13,43 @@ SYSTEM="M_1 1 M_2 1 separation 4 eccentricity 0.5"
############################################################
-PRE="time tbse orbital_period 0 $SYSTEM max_evolution_time 1000 $@"
+PRE="time tbse orbital_period 0 $SYSTEM max_evolution_time 10 repeat 1 $@"
+OUT="/tmp"
+OUT="$HOME/tex/notes/solvers"
+OUT="$HOME/tex/notes/solvers/various"
+OUT="$HOME/tex/notes/solvers/various0.1"
+OUT="$HOME/tex/notes/solvers/various0.01"
+OUT="$HOME/tex/notes/solvers/various0.001"
+OUT="$HOME/tex/notes/solvers/2nd_deriv"
+OUT="$HOME/tex/notes/solvers/anal"
+
+
+# make outdir
+mkdir -p $OUT
echo "Forward Euler"
-CMD="$PRE --log_filename /tmp/c_log.euler --solver 0"
+CMD="$PRE --log_filename $OUT/c_log.euler --solver 0"
echo $CMD
-$CMD | filter JJJ > /tmp/jjj.euler
+
+#$CMD maximum_timestep 1000 | filter JJJ > $OUT/jjj.1000.euler
+#$CMD maximum_timestep 1 | filter JJJ > $OUT/jjj.1.euler
+#$CMD maximum_timestep 0.1 | filter JJJ > $OUT/jjj.0.1.euler
+#$CMD maximum_timestep 0.01 | filter JJJ > $OUT/jjj.0.01.euler
+$CMD | filter JJJ > $OUT/jjj.euler
echo "RK2"
-CMD="$PRE --log_filename /tmp/c_log.RK2 --solver 1"
+CMD="$PRE --log_filename $OUT/c_log.RK2 --solver 1"
echo $CMD
-$CMD | filter JJJ > /tmp/jjj.RK2
+$CMD | filter JJJ > $OUT/jjj.RK2
echo "RK4"
-CMD="$PRE --log_filename /tmp/c_log.RK4 --solver 2"
+CMD="$PRE --log_filename $OUT/c_log.RK4 --solver 2"
echo $CMD
-$CMD | filter JJJ > /tmp/jjj.RK4
+$CMD | filter JJJ > $OUT/jjj.RK4
echo "PEC"
-CMD="$PRE --log_filename /tmp/c_log.PEC --solver 3"
+CMD="$PRE --log_filename $OUT/c_log.PEC --solver 3"
echo $CMD
-$CMD | filter JJJ > /tmp/jjj.PEC
+$CMD | filter JJJ > $OUT/jjj.PEC
# at 1Gyr: Euler = 2.75, RK2 = 2.9 : RK2 is more correct
--
GitLab