From 9233c17f1d7d201a2731eb1c70d1528f19175bde Mon Sep 17 00:00:00 2001
From: Robert Izzard <r.izzard@surrey.ac.uk>
Date: Mon, 12 Feb 2018 09:08:25 +0000
Subject: [PATCH] another attempt to fix the bug with discs only taking up
partial timesteps.
I've replaced disc->ftdisc (and removed it) with a local variable,
and the disc->Mdot,Jdot etc. derivatives are now calculated only
for the disc time. disc->dT is the disc time.
The conversion factor ftdisc is now set in the binary_c
derivative calculation function, which is really where it should
be set.
---
failsys | 2 ++
src/binary_c_structures.h | 2 +-
...ar_momentum_and_eccentricity_derivatives.c | 21 +++++++++--------
...ar_mass_and_angular_momentum_derivatives.c | 11 +++++++--
src/disc/disc_apply_derivatives.c | 23 ++++++++++++-------
5 files changed, 39 insertions(+), 20 deletions(-)
mode change 100644 => 100755 failsys
diff --git a/failsys b/failsys
old mode 100644
new mode 100755
index 855c5afae..b356670b3
--- a/failsys
+++ b/failsys
@@ -1 +1,3 @@
+#!/bin/bash
+
./binary_c --M_1 1.16771787551687 --M_2 1.12767845518499 --eccentricity 0 --orbital_period 1.15478198468946 --metallicity 0.02 --max_evolution_time 15000 --probability 8.15193521532146e-07 --AGB_3dup_algorithm 2 --AGB_core_algorithm 2 --AGB_luminosity_algorithm 2 --AGB_radius_algorithm 2 --BH_prescription 0 --Bondi_Hoyle_accretion_factor 1.5 --CRAP_parameter 0 --E2_prescription 0 --NeNaMgAl 1 --RLOF_interpolation_method 0 --RLOF_mdot_factor 1 --RLOF_method 0 --RLOF_timestep_modulator 1 --WRLOF_method 1 --accretion_end_time -1 --accretion_limit_dynamical_multiplier 1 --accretion_limit_eddington_multiplier 1 --accretion_limit_thermal_multiplier 1 --accretion_rate_novae_upper_limit 1.03e-07 --accretion_rate_soft_xray_upper_limit 2.71e-07 --accretion_start_time -1 --agbwindfac 1 --alphaCB 0 --alpha_ce 1 --angular_momentum_accretion_rate1 0 --angular_momentum_accretion_rate2 0 --angular_momentum_accretion_rate_orbit 0 --api_log_filename_prefix /dev/null --beta_reverse_nova 0 --c13_eff 1 --cbdisc_alpha 0.001 --cbdisc_eccentricity_pumping_method 1 --cbdisc_gamma 1.6666666666 --cbdisc_kappa 0.01 --cbdisc_mass_loss_FUV_multiplier 0 --cbdisc_mass_loss_ISM_pressure 3000 --cbdisc_mass_loss_ISM_ram_pressure_multiplier 1 --cbdisc_mass_loss_Xray_multiplier 1 --cbdisc_mass_loss_constant_rate 0 --cbdisc_mass_loss_inner_L2_cross_multiplier 0 --cbdisc_mass_loss_inner_viscous_accretion_method 1 --cbdisc_mass_loss_inner_viscous_multiplier 1 --cbdisc_minimum_fRing 0.2 --cbdisc_minimum_luminosity 0.0001 --cbdisc_minimum_mass 1e-06 --cbdisc_resonance_damping 1 --cbdisc_resonance_multiplier 1 --cbdisc_torquef 0.0001 --cbdisc_viscous_photoevaporation_coupling 1 --chandrasekhar_mass 1.44 --comenv_disc_angmom_fraction -3 --comenv_disc_mass_fraction 0.1 --comenv_ejection_spin_method 0 --comenv_merger_spin_method 0 --comenv_ms_accretion_mass 0 --comenv_ns_accretion_fraction 0 --comenv_ns_accretion_mass 0 --comenv_post_eccentricity 1e-05 --comenv_prescription 0 --comenv_splitmass 1.01 --delta_mcmin 0 --donor_limit_dynamical_multiplier 1 --donor_limit_thermal_multiplier 1 --dtfac 1 --escape_fraction 0 --escape_velocity 1000000000 --evolution_splitting 0 --evolution_splitting_maxdepth 2 --evolution_splitting_sn_n 10 --gb_reimers_eta 0.4 --gbwind 0 --gbwindfac 1 --gce_metallicity 0.02 --gravitational_radiation_model 0 --hachisu_disk_wind 0 --hachisu_qcrit 1.15 --hbbtfac 1 --idum 0 --individual_novae 0 --init_abunds_only 0 --initial_abundance_mix 0 --internal_buffering 2 --internal_buffering_compression 0 --jorb_RLOF_transfer_model 0 --lambda_ce -1 --lambda_enthalpy 0 --lambda_ionisation 0.1 --lambda_min 0 --lambda_multiplier 1 --lithium_GB_post_1DUP 0 --lithium_GB_post_Heflash 0 --lithium_hbb_multiplier 1 --log_filename /tmp/c_log.dat --lw 1 --magnetic_braking_factor 1 --mass_accretion_for_eld 1000000 --mass_accretion_rate1 0 --mass_accretion_rate2 0 --mass_for_Hestar_Ia_lower 0 --mass_for_Hestar_Ia_upper 0 --max_HeWD_mass 0.7 --max_neutron_star_mass 1.8 --max_tpagb_core_mass 1.38 --maximum_nuclear_burning_timestep 1e+20 --maximum_timestep 1e+20 --mc13_pocket_multiplier 1 --merger_angular_momentum_factor 1 --minimum_donor_menv_for_comenv 0.1 --minimum_envelope_mass_for_third_dredgeup 0.5 --minimum_mass_for_carbon_ignition 1.6 --minimum_mass_for_neon_ignition 2.85 --minimum_orbital_period_for_instant_RLOF 0 --minimum_separation_for_instant_RLOF 0 --minimum_timestep 1e-08 --monte_carlo_kicks --nelemans_gamma 1.75 --nelemans_max_frac_j_change 1 --nelemans_minq 0.2 --nelemans_n_comenvs 0 --nelemans_recalc_eccentricity 0 --no_thermohaline_mixing 0 --nonconservative_angmom_gamma -1 --nova_faml_multiplier 0 --nova_irradiation_multiplier 0 --nova_retention_fraction 0.001 --nova_retention_method 0 --nucsyn_metallicity 0.02 --phasevol 0.000307126809344586 --post_SN_orbit_method 0 --pre_main_sequence 0 --pre_main_sequence_fit_lobes 0 --qcrit_BH 3 --qcrit_CHeB 3 --qcrit_COWD 3 --qcrit_EAGB -1 --qcrit_GB -1 --qcrit_HG 4 --qcrit_HeGB 0.78125 --qcrit_HeHG 4 --qcrit_HeMS 3 --qcrit_HeWD 3 --qcrit_LMMS 0.6944 --qcrit_MS 1.6 --qcrit_NS 3 --qcrit_ONeWD 3 --qcrit_TPAGB -1 --qcrit_degenerate_BH 0.625 --qcrit_degenerate_CHeB 3 --qcrit_degenerate_COWD 0.625 --qcrit_degenerate_EAGB 1.15 --qcrit_degenerate_GB 1.15 --qcrit_degenerate_HG 4.7619 --qcrit_degenerate_HeGB 1.15 --qcrit_degenerate_HeHG 4.7619 --qcrit_degenerate_HeMS 3 --qcrit_degenerate_HeWD 0.625 --qcrit_degenerate_LMMS 1 --qcrit_degenerate_MS 1 --qcrit_degenerate_NS 0.625 --qcrit_degenerate_ONeWD 0.625 --qcrit_degenerate_TPAGB 1.15 --random_systems 0 --repeat 1 --rlperi 0 --rotationally_enhanced_exponent 1 --rotationally_enhanced_mass_loss 0 --sn_kick_companion_AIC 0 --sn_kick_companion_AIC_BH 0 --sn_kick_companion_BH_BH 0 --sn_kick_companion_BH_NS 0 --sn_kick_companion_ECAP 0 --sn_kick_companion_GRB_COLLAPSAR 0 --sn_kick_companion_HeStarIa 0 --sn_kick_companion_IA_CHAND 0 --sn_kick_companion_IA_CHAND_Coel 0 --sn_kick_companion_IA_ELD 0 --sn_kick_companion_IA_He 0 --sn_kick_companion_IA_He_Coel 0 --sn_kick_companion_IBC 0 --sn_kick_companion_II 0 --sn_kick_companion_IIa 0 --sn_kick_companion_NS_NS 0 --sn_kick_companion_TZ 0 --sn_kick_companion_WDKICK 0 --sn_kick_dispersion_AIC_BH 0 --sn_kick_dispersion_BH_BH 0 --sn_kick_dispersion_BH_NS 0 --sn_kick_dispersion_ECAP 0 --sn_kick_dispersion_GRB_COLLAPSAR 190 --sn_kick_dispersion_IBC 190 --sn_kick_dispersion_II 190 --sn_kick_dispersion_NS_NS 0 --sn_kick_dispersion_TZ 0 --sn_kick_distribution_AIC_BH 0 --sn_kick_distribution_BH_BH 0 --sn_kick_distribution_BH_NS 0 --sn_kick_distribution_ECAP 0 --sn_kick_distribution_GRB_COLLAPSAR 1 --sn_kick_distribution_IBC 1 --sn_kick_distribution_II 1 --sn_kick_distribution_NS_NS 0 --sn_kick_distribution_TZ 0 --superwind_mira_switchon 500 --third_dup 1 --tidal_strength_factor 1 --timestep_modulator 1 --tpagb_reimers_eta 1 --tpagbwind 0 --vrot1 0 --vrot2 0 --vw93_mira_shift 0 --vw93_multiplier 1 --vwind_multiplier 0.125 --wd_kick_direction 0 --wd_kick_pulse_number 0 --wd_kick_when 0 --wd_sigma 0 --wind_angular_momentum_loss 0 --wind_disc_angmom_fraction 0 --wind_disc_mass_fraction 0 --wind_djorb_fac 1 --wind_mass_loss 1 --wr_wind 0 --wr_wind_fac 1 --wrwindfac 1 --yields_dt 100000 --internal_buffering 0 --internal_buffering_compression 0 --log_filename /tmp/c_log2.dat
diff --git a/src/binary_c_structures.h b/src/binary_c_structures.h
index b18c8c42d..a5a0a6f99 100644
--- a/src/binary_c_structures.h
+++ b/src/binary_c_structures.h
@@ -789,7 +789,7 @@ struct disc_t {
double dM_binary,Mdot_binary;
double dM_ejected,Mdot_ejected;
double dJ_binary,Jdot_binary;
- double dT,ftdisc;
+ double dT;
double de_binary,edot_binary;
double Tvisc,Tradin,Tradout;
double Hin; // inner edge scale height
diff --git a/src/binary_star_functions/angular_momentum_and_eccentricity_derivatives.c b/src/binary_star_functions/angular_momentum_and_eccentricity_derivatives.c
index fe80c347c..80d207abf 100644
--- a/src/binary_star_functions/angular_momentum_and_eccentricity_derivatives.c
+++ b/src/binary_star_functions/angular_momentum_and_eccentricity_derivatives.c
@@ -85,12 +85,12 @@ void angular_momentum_and_eccentricity_derivatives(
*/
struct disc_t * disc = & stardata->common.discs[i];
- double dt = stardata->model.dtm * 1e6;
+ double dt = stardata->model.dtm * 1e6 * YEAR_LENGTH_IN_SECONDS;
+ double ftdisc = disc->dT / dt;
if(disc->dT > TINY && dt > TINY)
{
- double jdot = disc->Jdot_binary * disc->ftdisc *
+ double jdot = disc->Jdot_binary * ftdisc *
YEAR_LENGTH_IN_SECONDS / ANGULAR_MOMENTUM_CGS;
-
stardata->model.derivative[DERIVATIVE_ORBIT_ANGMOM_CBDISC] += jdot;
/*
@@ -101,14 +101,17 @@ void angular_momentum_and_eccentricity_derivatives(
* the disc that falls onto star 0 (the "primary")
*/
double f = disc_inner_edge_accretion_f(stardata);
- double mdot = disc->Mdot_binary * disc->ftdisc *
+ double mdot = disc->Mdot_binary * ftdisc *
YEAR_LENGTH_IN_SECONDS / M_SUN;
-
- if(0)printf("DMBIN deriv mdot %g Msun/y (ftdisc = %g/%g = %g)\n",
+
+ if(0)printf("DMBIN deriv mdot %g Msun/y (ftdisc = %g/%g = %g, disc->M=%g) : DM=%g\n",
mdot,
disc->dT/YEAR_LENGTH_IN_SECONDS,
- stardata->model.dtm*1e6,
- disc->ftdisc);
+ dt,
+ ftdisc,
+ disc->M/M_SUN,
+ mdot * stardata->model.dtm * 1e6
+ );
stardata->star[0].derivative[DERIVATIVE_STELLAR_MASS_CBDISC_GAIN] =
f * mdot;
@@ -121,7 +124,7 @@ void angular_momentum_and_eccentricity_derivatives(
* Resonant interaction with the circumbinary
* disc drives binary eccentricity
*/
- double edot = disc->edot_binary * YEAR_LENGTH_IN_SECONDS * disc->ftdisc;
+ double edot = disc->edot_binary * ftdisc * YEAR_LENGTH_IN_SECONDS;
stardata->model.derivative[DERIVATIVE_ORBIT_ECCENTRICITY_CBDISC] = edot;
diff --git a/src/binary_star_functions/apply_stellar_mass_and_angular_momentum_derivatives.c b/src/binary_star_functions/apply_stellar_mass_and_angular_momentum_derivatives.c
index 857761745..71727cd63 100644
--- a/src/binary_star_functions/apply_stellar_mass_and_angular_momentum_derivatives.c
+++ b/src/binary_star_functions/apply_stellar_mass_and_angular_momentum_derivatives.c
@@ -17,7 +17,7 @@ void apply_stellar_mass_and_angular_momentum_derivatives(struct stardata_t * sta
* Timestep
*/
double dt = RLOF==TRUE ? (stardata->model.dtm*1e6) : stardata->model.dt;
-
+
#if DEBUG==1
if(0&&Debug_expression)
show_derivatives(stardata);
@@ -52,8 +52,15 @@ void apply_stellar_mass_and_angular_momentum_derivatives(struct stardata_t * sta
double dM_in_timestep = star->derivative[DERIVATIVE_STELLAR_MASS]*dt;
star->mass += dM_in_timestep;
- if(star->mass > 90)
+ if(0&&star->mass > 90)
{
+ show_derivatives(stardata);
+ printf("star %d: M=%g dM=%g Mdot=%g\n",
+ star->starnum,
+ star->mass,
+ dM_in_timestep,
+ star->derivative[DERIVATIVE_STELLAR_MASS]
+ );
Exit_binary_c(2,"oops M>90");
}
diff --git a/src/disc/disc_apply_derivatives.c b/src/disc/disc_apply_derivatives.c
index eda0e8e77..1dd59a306 100644
--- a/src/disc/disc_apply_derivatives.c
+++ b/src/disc/disc_apply_derivatives.c
@@ -64,14 +64,15 @@ void disc_apply_derivatives(struct disc_t * disc,
disc->loss[DISC_LOSS_INNER_L2_CROSSING].mdot
));
- if(1)
- printf("DMBIN substep ddt=%g : integrated dM=%g Msun in dT=%g y, rate = %g Msun/y (ddt = %g y, disc->dT = %g y ) disc->M = %g Msun\n",
+ if(0)
+ printf("DMBIN substep ddt=%g : integrated dM=%g Msun in dT=%g y, rate = %g Msun/y (ddt = %g y, disc->dT = %g y, stardata dt = %g y) disc->M = %g Msun\n",
ddt/YEAR_LENGTH_IN_SECONDS,
disc->dM_binary/M_SUN,
disc->dT/YEAR_LENGTH_IN_SECONDS,
(disc->dM_binary * YEAR_LENGTH_IN_SECONDS) / (disc->dT * M_SUN),
ddt/YEAR_LENGTH_IN_SECONDS,
disc->dT/YEAR_LENGTH_IN_SECONDS,
+ stardata->model.dtm * 1e6,
disc->M/M_SUN
);
@@ -171,16 +172,22 @@ void disc_apply_derivatives(struct disc_t * disc,
*/
if(disc->dT > TINY)
{
- disc->ftdisc = disc->dT / (YEAR_LENGTH_IN_SECONDS * 1e6 * stardata->model.dtm);
+
+ /*
+ * The derivatives are only valid while the disc is alive.
+ * This is taken into account during the derivative
+ * application process (with its factor ftdisc).
+ */
disc->Mdot_ejected = disc->dM_ejected / disc->dT;
disc->Mdot_binary = disc->dM_binary / disc->dT;
disc->Jdot_binary = disc->dJ_binary / disc->dT;
disc->edot_binary = disc->de_binary / disc->dT;
- printf("disc dots Mej=%g Mbin=%g : ftdisc = %g\n",
- disc->Mdot_ejected * YEAR_LENGTH_IN_SECONDS / M_SUN,
- disc->Mdot_binary * YEAR_LENGTH_IN_SECONDS / M_SUN,
- disc->ftdisc
- );
+ if(0)
+ printf("disc dMej/dt=%g dMbin/dt=%g Msun/y for dT=%g y\n",
+ disc->Mdot_ejected * YEAR_LENGTH_IN_SECONDS / M_SUN,
+ disc->Mdot_binary * YEAR_LENGTH_IN_SECONDS / M_SUN,
+ disc->dT / YEAR_LENGTH_IN_SECONDS
+ );
}
}
--
GitLab