update solver test script

fix timesteps to correct for the wrong size of the arrays

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

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",

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)

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

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)]                  \
          -                                      \

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

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