add albedo option to discs

add pthread warnings

src/binary_c_function_macros.h

@@ -630,4 +630,19 @@
  */
 #define Loop_forever for(;;)
 
+
+
+#ifdef __HAVE_LIBPTHREAD__
+/*
+ * libpthread errors as strings
+ */
+#define Mutex_error_string(X)                   \
+    (X) == EINVAL ? "EINVAL" :                  \
+        (X) == EBUSY ? "EBUSY" :                \
+        (X) == EAGAIN ? "EAGAIN" :              \
+        (X) == EPERM ? "EPERM" :                \
+        (X) == EDEADLK ? "EDEADLK"  :           \
+        "unknown"
+#endif // __HAVE_LIBPTHREAD__
+
 #endif // BINARY_C_FUNCTION_MACROS_H

src/binary_c_structures.h

@@ -930,7 +930,7 @@ struct preferences_t {
     double cbdisc_torquef;
     double cbdisc_init_dM;
     double cbdisc_init_dJdM;
-
+    double cbdisc_albedo;
 #ifdef DISCS_CIRCUMBINARY_FROM_COMENV
     double comenv_disc_angmom_fraction;
     double comenv_disc_mass_fraction;
@@ -1436,6 +1436,7 @@ struct disc_t {
     double M; // disc total mass
     double alpha; // disc viscosity constant (dimensionless)
     double gamma; // adiabatic exponent (dimensionless)
+    double albedo; // albedo
     double torqueF; // binary torque multiplication factor (dimensionless)
     double torqueF_for_Rin; // binary torque multiplication factor used to fix Rin (dimensionless)
     double fuv; // Far UV mass flux in g/cm^2/s

src/binary_star_functions/apply_other_stellar_derivatives.c

@@ -48,7 +48,6 @@ Boolean apply_other_stellar_derivatives(struct stardata_t * RESTRICT const stard
     {
         return FALSE;
     }
-
 #endif //MINT
 
     return TRUE;

src/disc/cbdisc_eccentricity_pumping_rate.c

@@ -156,23 +156,16 @@ void cbdisc_eccentricity_pumping_rate(struct stardata_t * stardata,
              * Note that the l=1, m=2 resonance which is the same
              * as the 1:3 EOLR (eccentric outer Lindblad resonance)
              */
-            if(disc->Rin < Rres_out)
-            {
-                deda = 2.0 * ee / (e + alpha / (100.0 * e)) *
-                    (l/m - 1.0/sqrt(ee)) / a;
-                using = 9;
-                Discdebug(2,"Rin = %g < Rres_Out = %g : strong coupling e = %g -> de/da = %g\n",
-                          disc->Rin,
-                          Rres_out,
-                          e,
-                          deda);
-            }
-            else
-            {
-                Discdebug(2,"Rin = %g >= Rres_Out = %g : no pumping\n",
-                          disc->Rin,
-                          Rres_out);
-            }
+
+            deda = 2.0 * ee / (e + alpha / (100.0 * e)) *
+                (l/m - 1.0/sqrt(ee)) / a;
+            using = 9;
+            Discdebug(2,"Rin = %g < Rres_Out = %g : strong coupling e = %g -> de/da = %g\n",
+                      disc->Rin,
+                      Rres_out,
+                      e,
+                      deda);
+
         }
 
         if(stardata->preferences->cbdisc_resonance_damping == TRUE)

src/disc/disc_build_disc_zones.c

@@ -9,7 +9,9 @@ int disc_build_disc_zones(struct disc_t * const disc,
     /*
      * Make a new zone list in a disc and zone it.
      */
-    const Boolean status = disc_new_zone_list(disc,binary,disc->thermal_zones);
+    const Boolean status = disc_new_zone_list(disc,
+                                              binary,
+                                              disc->thermal_zones);
 
     Discdebug(2,"new_zone_list Rin = %g, Rout = %g : ok ? %d\n",disc->Rin,disc->Rout,status);
 

src/disc/disc_edge_stripping.c

@@ -68,8 +68,6 @@ void disc_edge_stripping(struct stardata_t * const stardata,
             }
             else
             {
-                printf("Disc moves inwards\n");
-
                 /*
                  * Revap_in moves inwards
                  *

src/disc/disc_evolve_disc_structure.c

@@ -498,8 +498,11 @@ static void show_test_expressions(struct binary_system_t * binary,
                                   struct disc_t * disc)
 {
     double B =
-        pow(Pow2(binary->Rstar) * (binary->flux/(4.0*STEFAN_BOLTZMANN_CONSTANT))                                         *sqrt(disc->gamma*BOLTZMANN_CONSTANT/
-                                                                                                                               (GRAVITATIONAL_CONSTANT*binary->mtot*M_PROTON)),
+        pow(Pow2(binary->Rstar) *
+            (binary->flux*(1.0 - disc->albedo)
+             /(4.0*STEFAN_BOLTZMANN_CONSTANT))
+            *sqrt(disc->gamma*BOLTZMANN_CONSTANT/
+                  (GRAVITATIONAL_CONSTANT*binary->mtot*M_PROTON)),
             2.0/7.0);
     double Tin = disc_temperature(disc->Rin,disc);
     Discdebug(2,

src/disc/disc_init_binary_structure.c

@@ -77,6 +77,7 @@ void disc_init_binary_structure(const struct stardata_t * const stardata,
      */
     binary->Rstar = (FR3(0) + FR3(1))/(FR2(0) + FR2(1));
     binary->flux = binary->L/(4.0 * PI * Pow2(binary->Rstar));
+
     /*
      * RLOFing check
      */

src/disc/disc_set_disc_initial_guesses.c

@@ -93,7 +93,8 @@ void disc_set_disc_initial_guesses(int n,/* VLA initial_guesses REQUIRES this to
         disc->Rout = 1e3 * disc->Rin;
         disc->Tvisc0 = 1e15;
     }
-    disc_initial_radiative_guesses(binary,disc);
+    disc_initial_radiative_guesses(binary,
+                                   disc);
     Set_guess(3,
               disc->Tvisc0,
               disc->Rin,

src/disc/disc_structure_functions.c

@@ -841,7 +841,7 @@ void disc_equation_checker(const double r,
     }
     const double B = GRAVITATIONAL_CONSTANT * binary->mtot * mdot /
         (2.0 * PI * binary->separation * Pow2(r) * STEFAN_BOLTZMANN_CONSTANT);
-    const double C = 2.0/(3.0*PI) * binary->flux * Pow3(binary->Rstar / r) / STEFAN_BOLTZMANN_CONSTANT;
+    const double C = 2.0/(3.0*PI) * binary->flux * (1.0 - disc->albedo) * Pow3(binary->Rstar / r) / STEFAN_BOLTZMANN_CONSTANT;
     const double D = 3.0 * PI * H / (4.0 * binary->Rstar) * (dlnHdlnR - 1.0);
     const double T4_RHS = A + 0.0*B + C * (1.0 + D);
 #ifdef DISC_EQUATION_CHECKS_LOG

src/disc/disc_thermal_zones.h

@@ -116,7 +116,9 @@
     (                                                                   \
         pow(                                                            \
             1.0/7.0 * Pow2(binary->Rstar) *                             \
-            (binary->flux/STEFAN_BOLTZMANN_CONSTANT)                    \
+            (binary->flux *                                             \
+             (1.0 - disc->albedo)/                                      \
+             STEFAN_BOLTZMANN_CONSTANT)                                 \
             *sqrt(disc->gamma*BOLTZMANN_CONSTANT/                       \
                   (GRAVITATIONAL_CONSTANT*binary->mtot*disc->mu)),      \
             2.0/7.0)                                                    \
@@ -134,12 +136,13 @@
 #define Exponent_massloss (-1.0/2.0)
 
 // d
-#define Prefactor_radiative_in                          \
-    (                                                   \
-        Pow1d4(                                         \
-            2.0 * Pow3(binary->Rstar) * binary->flux /  \
-            (3.0*PI*STEFAN_BOLTZMANN_CONSTANT)          \
-            )                                           \
+#define Prefactor_radiative_in                              \
+    (                                                       \
+        Pow1d4(                                             \
+            2.0 * Pow3(binary->Rstar) * binary->flux *      \
+            (1.0 - disc->albedo)/   \
+            (3.0*PI*STEFAN_BOLTZMANN_CONSTANT)              \
+            )                                               \
         )
 #define Exponent_radiative_in (-3.0/4.0)
 

src/disc/new_disc.c

@@ -126,6 +126,7 @@ struct disc_t * new_disc(
             cbdisc->gamma = stardata->preferences->cbdisc_gamma;
             cbdisc->torqueF = stardata->preferences->cbdisc_torquef;
             cbdisc->kappa = stardata->preferences->cbdisc_kappa;
+            cbdisc->albedo = stardata->preferences->cbdisc_albedo;
         }
     }
     else

src/evolution/apply_derivatives.c

@@ -4,7 +4,6 @@ No_empty_translation_unit_warning;
 /*
  * Apply derivatives
  */
-
 enum derivatives {
     ORBITAL,
     STELLAR,

src/evolution/evolution_split.c

@@ -62,7 +62,7 @@ static void modulate_stardata(struct stardata_t * RESTRICT const stardata,
 
 // debugging output if required
 #define SDEBUG Dprint
-//#define SDEBUG(...) if(stardata->model.time  > 5.95) { printf("split : ");printf(__VA_ARGS__); }
+//#define SDEBUG(...) if(stardata->model.time  > 0) { printf("split : ");printf(__VA_ARGS__); }
 
 #ifdef NUCSYN
 #if defined NUCSYN_YIELD_COMPRESSION || defined NUCSYN_LOG_BINARY_DX_YIELDS
@@ -81,11 +81,16 @@ int evolution_split(struct stardata_t * RESTRICT const stardata,
 #ifdef EVOLUTION_SPLITTING
     if(stardata->preferences->evolution_splitting)
     {
-        SDEBUG("splitting is on at t=%g : status = %d cf %d (break? %s): stellar types %d %d : splittest = %d, doingsplit = %d, In_splitdepth = %d (current depth %d vs Min(%d,%d))\n",
+#define Break_evolution_Q                               \
+        (abs(*status) == EVOLUTION_BREAK ||             \
+         abs(*status) == EVOLUTION_STOP)
+
+        SDEBUG("splitting is on at t=%g : status = %d cf %d %d (break? %s): stellar types %d %d : splittest = %d, doingsplit = %d, In_splitdepth = %d (current depth %d vs Min(%d,%d))\n",
                stardata->model.time,
                *status,
                EVOLUTION_BREAK,
-               Yesno(abs(*status)==EVOLUTION_BREAK),
+               EVOLUTION_STOP,
+               Yesno(Break_evolution_Q),
                stardata->star[0].stellar_type,
                stardata->star[1].stellar_type,
                splittest(stardata),
@@ -96,7 +101,7 @@ int evolution_split(struct stardata_t * RESTRICT const stardata,
                EVOLUTION_SPLITTING_HARD_MAX_DEPTH
             );
 
-        if(abs(*status) == EVOLUTION_BREAK)
+        if(Break_evolution_Q)
         {
             /*
              * flush yield logs of previous evolution, except
@@ -120,7 +125,7 @@ int evolution_split(struct stardata_t * RESTRICT const stardata,
                  * Check if we need to rerun from a previous
                  * splitpoint
                  */
-                SDEBUG("SPLIT status is evolution_break : currently working on splitinfo struct %d at %p, of which count = %d\n",
+                SDEBUG("SPLIT status is evolution_break : currently working on split-info struct %d at %p, of which count = %d\n",
                        stardata->preferences->current_splitdepth-1,
                        (void*)p,
                        p->count);
@@ -212,7 +217,7 @@ int evolution_split(struct stardata_t * RESTRICT const stardata,
                     = stardata->preferences->splitinfo[stardata->preferences->current_splitdepth]
                     = Malloc(sizeof(struct splitinfo_t));
 
-                SDEBUG("malloced new splitinfo struct at %p\n",(void*)p);
+                SDEBUG("malloced new split-info struct at %p\n",(void*)p);
 
                 p->depth = stardata->preferences->current_splitdepth;
 

src/memory/build_store_contents.c

@@ -28,14 +28,33 @@ void build_store_contents(struct store_t * RESTRICT const store)
      */
     for(size_t i=0; i<BINARY_C_MUTEX_NUMBER; i++)
     {
-        if(pthread_mutex_init(&store->pthread_mutexes[i], NULL) != 0)
+        const int init_status =
+            pthread_mutex_init(&store->pthread_mutexes[i],
+                               NULL);
+        if(init_status != 0)
         {
             fprintf(stderr,
-                    "\n failed to initialize pthread mutex %zu : this will cause problems!\n",
-                    i);
+                    "\n failed to initialize pthread mutex %zu, error %d %s : this will cause problems!\n",
+                    i,
+                    init_status,
+                    Mutex_error_string(init_status)
+                );
         }
     }
-    pthread_mutex_lock(&store->pthread_mutexes[BINARY_C_MUTEX_STORE_LOCK]);
+
+    {
+        const int lock_status =
+            pthread_mutex_lock(&store->pthread_mutexes[BINARY_C_MUTEX_STORE_LOCK]);
+        if(lock_status != 0)
+        {
+            fprintf(stderr,
+                    "\n failed to lock mutex BINARY_C_MUTEX_STORE_LOCK == %d, error %d %s.\n",
+                    BINARY_C_MUTEX_STORE_LOCK,
+                    lock_status,
+                    Mutex_error_string(lock_status));
+        }
+    }
+
 #endif//__HAVE_LIBPTHREAD__
 
     if(store->built==FALSE)
@@ -118,6 +137,17 @@ void build_store_contents(struct store_t * RESTRICT const store)
         store->debug_stopping = 0;
     }
 #ifdef __HAVE_LIBPTHREAD__
-    pthread_mutex_unlock(&store->pthread_mutexes[BINARY_C_MUTEX_STORE_LOCK]);
+    {
+        const int unlock_status =
+            pthread_mutex_unlock(&store->pthread_mutexes[BINARY_C_MUTEX_STORE_LOCK]);
+        if(unlock_status != 0)
+        {
+            fprintf(stderr,
+                    "\n failed to unlock mutex BINARY_C_MUTEX_STORE_LOCK == %d, error %d %s.\n",
+                    BINARY_C_MUTEX_STORE_LOCK,
+                    unlock_status,
+                    Mutex_error_string(unlock_status));
+        }
+    }
 #endif//__HAVE_PTHREAD__
 }

src/memory/free_store_contents.c

@@ -14,11 +14,8 @@ void free_store_contents(struct store_t * RESTRICT const store)
      * NOTHING in here can depend on the current stardata.
      * You will NOT have access to stardata here.
      */
-
     if(store)
     {
-
-
 #ifdef NUCSYN
         nucsyn_free_store_contents(store);
 #endif
@@ -90,17 +87,35 @@ void free_store_contents(struct store_t * RESTRICT const store)
 
 
 #ifdef __HAVE_LIBPTHREAD__
+#define Mutex_error_string(X)                   \
+        (X) == EINVAL ? "EINVAL" :              \
+            (X) == EBUSY ? "EBUSY" :            \
+            (X) == EAGAIN ? "EAGAIN" :          \
+            (X) == EPERM ? "EPERM" :            \
+            (X) == EDEADLK ? "EDEADLK"  :       \
+            "unknown"
+
         /*
          * Make mutexes for later locking
          */
         for(size_t i=0; i<BINARY_C_MUTEX_NUMBER; i++)
         {
-            pthread_mutex_unlock(&store->pthread_mutexes[i]);
-            if(pthread_mutex_destroy(&store->pthread_mutexes[i]) != 0)
+            const int unlock_status = pthread_mutex_unlock(&store->pthread_mutexes[i]);
+            if(unlock_status != 0)
+            {
+                fprintf(stderr,"\n mutex %zu unlock failed with error %d %s\n",
+                        i,
+                        unlock_status,
+                        Mutex_error_string(unlock_status));
+            }
+            const int destroy_status = pthread_mutex_destroy(&store->pthread_mutexes[i]);
+            if(destroy_status != 0)
             {
                 fprintf(stderr,
-                        "\n failed to destroy pthread mutex %zu : may cause problems\n",
-                        i);
+                        "\n failed to destroy pthread mutex %zu error %d %s: may cause problems\n",
+                        i,
+                        destroy_status,
+                        Mutex_error_string(destroy_status));
             }
         }
 #endif//__HAVE_LIBPTHREAD__

src/setup/cmd_line_args_list.h

@@ -2945,6 +2945,15 @@ BATCHMODE_ARGS
         Var(stardata->preferences->lambda_enthalpy),
         1.0
         CMD_LINE_ARG_T_OPTIONAL_ARGS },
+{
+    ARG_SECTION_BINARY,
+        "cbdisc_albedo",
+        "Circumbinary-disc albedo. Requires DISCS.",
+        ARG_FLOAT ,
+        WTTS_USE_DEFAULT,
+        CBDISC_ALBEDO_VAR,
+        1.0
+        CMD_LINE_ARG_T_OPTIONAL_ARGS },
 {
     ARG_SECTION_BINARY,
         "cbdisc_gamma",

src/setup/cmd_line_macros.h

@@ -326,6 +326,7 @@
 #define DISC_TIMESTEP_FACTOR_VAR Var(stardata->preferences->disc_timestep_factor)
 #define DISC_N_MONTE_CARLO_GUESSES_VAR Var(stardata->preferences->disc_n_monte_carlo_guesses)
 #define CBDISC_GAMMA_VAR Var(stardata->preferences->cbdisc_gamma)
+#define CBDISC_ALBEDO_VAR Var(stardata->preferences->cbdisc_albedo)
 #define CBDISC_MAX_LIFETIME_VAR Var(stardata->preferences->cbdisc_max_lifetime)
 #define CBDISC_ALPHA_VAR Var(stardata->preferences->cbdisc_alpha)
 #define CBDISC_KAPPA_VAR Var(stardata->preferences->cbdisc_kappa)
@@ -382,6 +383,7 @@
 #define DISC_TIMESTEP_FACTOR_VAR Not_used_arg
 #define DISC_N_MONTE_CARLO_GUESSES_VAR Not_used_arg
 #define CBDISC_GAMMA_VAR Not_used_arg
+#define CBDISC_ALBEDO_VAR Not_used_arg
 #define CBDISC_MAX_LIFETIME_VAR Not_used_arg
 #define CBDISC_ALPHA_VAR Not_used_arg
 #define CBDISC_KAPPA_VAR Not_used_arg

src/setup/set_default_preferences.c

@@ -437,6 +437,7 @@ preferences->WD_accretion_rate_novae_upper_limit_other_donor = ACCRETION_RATE_NO
     preferences->cbdisc_max_lifetime = 0.0;
     preferences->cbdisc_alpha = 1e-3;
     preferences->cbdisc_kappa = 1e-2;
+    preferences->cbdisc_albedo = 0.0;
     preferences->cbdisc_torquef = 1e-3;
     preferences->cbdisc_init_dM = 0.1;
     preferences->cbdisc_init_dJdM = 0.5;

tbse

@@ -730,6 +730,7 @@ CBDISC_GAMMA=1.4 # 1.4 (molecular hydrogen) or 1.6666666 (monatomic gas)
 CBDISC_ALPHA=1e-3 # 1e-2 or 1e-3
 CBDISC_KAPPA=1e-2 # cgs (cm^2 g^-1) : 1e-2 for cool gas, 1e5 for dust in the visible cf https://ned.ipac.caltech.edu/level5/Sept05/Li2/Li4.html
 CBDISC_TORQUEF=1e-3 # 1e-3 tends to keep discs outside L2
+CBDISC_ALBEDO=0.0
 
 ##############################
 # Initial structure options
@@ -2376,6 +2377,7 @@ $WARMUP \
 --cbdisc_alpha $CBDISC_ALPHA \
 --cbdisc_kappa $CBDISC_KAPPA \
 --cbdisc_torquef $CBDISC_TORQUEF \
+--cbdisc_albedo $CBDISC_ALBEDO \
 --cbdisc_init_dM $CBDISC_INIT_DM \
 --cbdisc_init_dJdM $CBDISC_INIT_DJDM  \
 --cbdisc_minimum_evaporation_timescale $CBDISC_MINIMUM_EVAPORATION_TIMESCALE \