From 3d2a73ee581b6949408221222347d6534c100ccc Mon Sep 17 00:00:00 2001 From: "gvlipunova@gmail.com" Date: Fri, 19 Apr 2024 23:16:22 +0200 Subject: [PATCH] =?UTF-8?q?1)=20changed=20default=5Fopacity=20from=20"Kram?= =?UTF-8?q?ers"=20to=20"OPAL"=202)=20renamed=20input=20parameter=20Rhot=5F?= =?UTF-8?q?Mdotzero=5Ffactor=20->=20Rfront=5FMdotzero=5Ffactor=20(is=20not?= =?UTF-8?q?=20important=20variable=20now:=20its=20recommended=20value=201?= =?UTF-8?q?=20is=20default)=203)=20added=20input=20parameter=20DIM=5Ffront?= =?UTF-8?q?=5FMdot=5Ffactor=20(default=202.3).=20At=20the=20border=20betwe?= =?UTF-8?q?en=20hod=20and=20cold=20zones,=20there=20is=20a=20negative=20ac?= =?UTF-8?q?cretion=20rate=20Mdot=20=3D=20-=20DIM=5Ffront=5FMdot=5Ffactor?= =?UTF-8?q?=20*=20Mdot(Rin)=204)=20added=20input=20parameter=20DIM=5Ffront?= =?UTF-8?q?=5Fapproach=20(default=20"maxFvis")=20-=20the=20disc=20equation?= =?UTF-8?q?=20is=20solved=20to=20the=20point=20where=20Mdot=3D0=20(default?= =?UTF-8?q?)=20or=20where=20Mdot=20=3D=20-=20DIM=5Ffront=5FMdot=5Ffactor?= =?UTF-8?q?=20*=20Mdot(Rin)=205)=20added=20input=20parameter=20scatter=5Fb?= =?UTF-8?q?y=5Fcorona=20(yes/no)=206)=20added=20DiscEqFailException:=20pub?= =?UTF-8?q?lic=20std::exception=20to=20exceptions.hpp=20to=20report=20abou?= =?UTF-8?q?t=20non-converging=20solution=207)=20added=20household=20variab?= =?UTF-8?q?les=20R=5FdotM0=5Fbefore=5Fshift,=20maxR=5FQirr=5Fno=5Frole=20t?= =?UTF-8?q?o=20freddi=5Fstate=208)=20added=20boost::optional=20Shado?= =?UTF-8?q?w=20:=20ring=20is=20shadowed=20(1)=20or=20not=20(0,=20default)?= =?UTF-8?q?=209)=20added=20inline=20void=20set=5FMdot=5Fouter=5Fboundary()?= =?UTF-8?q?=20to=20change=20current=5F.Mdot=5Fouter=5Fboundary=2010)=20add?= =?UTF-8?q?ed=20inline=20void=20set=5FR=5FdotM0=5Fbefore=5Fshift()=20to=20?= =?UTF-8?q?change=20current=5F.R=5FdotM0=5Fbefore=5Fshift=2011)=20added=20?= =?UTF-8?q?inline=20void=20set=5FmaxR=5FQirr=5Fno=5Frole()=20to=20change?= =?UTF-8?q?=20current=5F.maxR=5FQirr=5Fno=5Frole=2012)=20added=20functions?= =?UTF-8?q?:=20Rfront=5FRhot(double=20r,=20double=20z=5Fr),=20obtain=5FMdo?= =?UTF-8?q?t=5Fouter=5Fboundary(),=20Tirr=5Fcritical(double=20r,=20int=20i?= =?UTF-8?q?i),=20and=20verify=5Fdisc=5Fmass(double=20tau)=2013)=20added=20?= =?UTF-8?q?2nd=20argument=20to=20=20R=5Fcooling=5Ffront(double=20r,=20doub?= =?UTF-8?q?le=20sigma=5Fat=5Fr)=20and=20v=5Fcooling=5Ffront(double=20r,=20?= =?UTF-8?q?double=20sigma=5Fat=5Fr);=2014)=20added=20functions:=20double?= =?UTF-8?q?=20R=5Fvis=5Fstruct=20(double=20r,=20double=20z2r=5Fat=5Fr)=20:?= =?UTF-8?q?=20estimate=20the=20radius,=20to=20which=20viscous=20process=20?= =?UTF-8?q?can=20reach=20in=20time=20tau=2015)=20changed=20name=20of=20fun?= =?UTF-8?q?ction:=20ring=5Fstate=5Fvertical()=20->=20check=5Fring=5Fis=5Fc?= =?UTF-8?q?old()=2016)=20in=20freddi=5Fevolution.cpp:=20changed=20name=20o?= =?UTF-8?q?f=20function:=20Mdot=5Fout()=20->=20Mdot=5Fouter=5Fboundary()?= =?UTF-8?q?=2017)=20in=20freddi=5Fevolution.cpp:=20moved=20all=20condition?= =?UTF-8?q?s=20to=20find=20boundary=20of=20the=20hot=20disc=20to=20functio?= =?UTF-8?q?n=20check=5Fring=5Fis=5Fcold();=201=20:=20ring=20is=20cold=2018?= =?UTF-8?q?)=20removed=20Tirr=5Fexceed=5Fcritical()=2019)=20in=20freddi=5F?= =?UTF-8?q?state.cpp:=20added=20vecd&=20FreddiState::Shadow()=20;=20double?= =?UTF-8?q?=20v=5Fvisc(double=20r,=20double=20z2r)=2020)=20in=20freddi=5Fs?= =?UTF-8?q?tate.cpp:=20in=20const=20vecd&=20FreddiState::Qx()=20:=20added?= =?UTF-8?q?=20(1.0=20-=20Shad[i])=20to=20Qx=2021)=20in=20v=5Fcooling=5Ffro?= =?UTF-8?q?nt()=20added=20check=20that=20if=20(sigma<0.)=20{sigma=3D0.0;}?= =?UTF-8?q?=20in=20the=20Ludwig's=20scheme=2022)=20in=20double=20FreddiSta?= =?UTF-8?q?te::Teff=5Fplus(double=20r)=20const=20minimum=20temperature=20o?= =?UTF-8?q?n=20the=20hot=20branch:=20changed=20from=20Lasota=20et=20al.,?= =?UTF-8?q?=20A&A=20486,=20523=E2=80=93528=20(2008)=20to=20=20Tavleev+23?= =?UTF-8?q?=2023)=20in=20function=20=20R=5Fcooling=5Ffront=20:=20changed?= =?UTF-8?q?=20result:=20R()[last()]=20-=20v=5Fcooling=5Ffront(r)=20*=20arg?= =?UTF-8?q?s().calc->tau;=20->=20r=20-=20v=5Fcooling=5Ffront(r,=20sigma=5F?= =?UTF-8?q?at=5Fr)=20*=20args().calc->tau;=2024)=20in=20freddi=5Fstate.cpp?= =?UTF-8?q?:=20added=20function=20double=20FreddiState::Rfront=5FRhot=20-?= =?UTF-8?q?=20but=20it=20is=20not=20used=20now=20because=20Rfront=5FMdotze?= =?UTF-8?q?ro=5Ffactor=20=3D=3D=201.0=20is=20recommended=2025)=20added=20d?= =?UTF-8?q?ouble=20FreddiState::Tirr=5Fcritical=20(double=20r,=20int=20ii)?= =?UTF-8?q?=20to=20determine=20critical=20level=20Tirr=20which=20keeps=20t?= =?UTF-8?q?he=20ring=20hot.=2026)=20in=20freddi=5Fstate.cpp:=20changed=20n?= =?UTF-8?q?ame=20of=20function=20ring=5Fstate=5Fvertical()=20to=20=20=20Fr?= =?UTF-8?q?eddiState::check=5Fring=5Fis=5Fcold()=20:=20=20returns=201=20fo?= =?UTF-8?q?r=20hot,=200=20for=20cold=2027)=20in=20nonlinear=5Fdiffusion.cp?= =?UTF-8?q?p:=20added=20include=20"exceptions.hpp"=2028)=20in=20nonlinear?= =?UTF-8?q?=5Fdiffusion.cpp:=20added=20vatriables=20iter=5Fsol=20and=20max?= =?UTF-8?q?iter=20to=20keep=20track=20on=20the=20number=20of=20iterations?= =?UTF-8?q?=20:=20if=20DIM=5Ffront=5Fapproach=20is=20"outflow"=20(accretio?= =?UTF-8?q?n=20rate=20is=20negative=20at=20Rhot),=20disk=20equation=20some?= =?UTF-8?q?time=20does=20not=20converge=2029)=20in=20ns=5Fevolution.cpp:?= =?UTF-8?q?=20added=20std=20output=20to=20exit=20from=20error:=20"R=5Fdead?= =?UTF-8?q?=20is=20positive=20and=20less=20than=20R=5Fcor"=2030)=20in=20ou?= =?UTF-8?q?tput.cpp:=20added=20to=20Derived=20values=20freddi->args().basi?= =?UTF-8?q?c->rin=2031)=20in=20output.cpp:=20added=20to=20BasicFreddiFileO?= =?UTF-8?q?utput::diskStructureDump:=20###=20Mdot=20=3D=20freddi->Mdot=5Fi?= =?UTF-8?q?n()=2032)=20input=20parameter=20DIM=5Ffront=5Fapproach=20is=20e?= =?UTF-8?q?ither=20"outflow"=20of=20"maxFvis"=20(default)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- cpp/include/arguments.hpp | 64 +++- cpp/include/exceptions.hpp | 7 + cpp/include/freddi_state.hpp | 27 +- cpp/include/nonlinear_diffusion.hpp | 2 +- cpp/include/ns/ns_arguments.hpp | 3 +- cpp/include/orbit.hpp | 2 +- cpp/src/arguments.cpp | 12 +- cpp/src/freddi_evolution.cpp | 181 ++--------- cpp/src/freddi_state.cpp | 465 +++++++++++++++++++++++++++- cpp/src/nonlinear_diffusion.cpp | 11 +- cpp/src/ns/ns_arguments.cpp | 10 +- cpp/src/ns/ns_evolution.cpp | 6 +- cpp/src/ns/ns_options.cpp | 7 +- cpp/src/opacity_related.cpp | 4 + cpp/src/options.cpp | 32 +- cpp/src/output.cpp | 21 +- 16 files changed, 632 insertions(+), 222 deletions(-) diff --git a/cpp/include/arguments.hpp b/cpp/include/arguments.hpp index f47185a..98f5fb3 100644 --- a/cpp/include/arguments.hpp +++ b/cpp/include/arguments.hpp @@ -175,15 +175,18 @@ class DiskStructureArguments { vecd operator()(const vecd& h) const override; }; public: - constexpr static const char default_opacity[] = "Kramers"; + constexpr static const char default_opacity[] = "OPAL"; constexpr static const double default_Mdotout = 0.; constexpr static const char default_boundcond[] = "Teff"; constexpr static const double default_Thot = 0.; constexpr static const double default_Tirr2Tvishot = 0.; - constexpr static const double default_Rhot_Mdotzero_factor = 1.; + constexpr static const double default_Rfront_Mdotzero_factor = 1.; + constexpr static const double default_DIM_front_Mdot_factor = 2.3; constexpr static const char default_check_state_approach[] = "before2024"; constexpr static const char default_check_Sigma_approach[] = "simple"; constexpr static const char default_check_Temp_approach[] = "const"; + constexpr static const char default_DIM_front_approach[] = "maxFvis"; + constexpr static const char default_scatter_by_corona[] = "yes"; constexpr static const char default_initialcond[] = "powerF"; constexpr static const char default_wind[] = "no"; public: @@ -195,10 +198,13 @@ class DiskStructureArguments { std::string boundcond; double Thot; double Tirr2Tvishot; - double Rhot_Mdotzero_factor; + double Rfront_Mdotzero_factor; + double DIM_front_Mdot_factor; std::string check_state_approach; std::string check_Sigma_approach; std::string check_Temp_approach; + std::string DIM_front_approach; + std::string scatter_by_corona; double F0; double Mdisk0; double Mdot0; @@ -219,14 +225,25 @@ class DiskStructureArguments { const BasicDiskBinaryArguments &bdb_args, const std::string& opacity, double Mdotout, - const std::string& boundcond, double Thot, double Tirr2Tvishot, - double Rhot_Mdotzero_factor, const std::string& check_state_approach, const std::string& check_Sigma_approach, const std::string& check_Temp_approach, + const std::string& boundcond, + double Thot, + double Tirr2Tvishot, + double Rfront_Mdotzero_factor, + double DIM_front_Mdot_factor, + const std::string& check_state_approach, + const std::string& check_Sigma_approach, + const std::string& check_Temp_approach, + const std::string& DIM_front_approach, + const std::string& scatter_by_corona, const std::string& initialcond, std::optional F0, - std::optional Mdisk0, std::optional Mdot0, + std::optional Mdisk0, + std::optional Mdot0, std::optional powerorder, - std::optional gaussmu, std::optional gausssigma, - const std::string& wind, const pard& windparams); + std::optional gaussmu, + std::optional gausssigma, + const std::string& wind, + const pard& windparams); DiskStructureArguments( const std::string &opacity, const OpacityRelated &oprel, @@ -234,23 +251,34 @@ class DiskStructureArguments { const std::string &boundcond, double Thot, double Tirr2Tvishot, - double Rhot_Mdotzero_factor, + double Rfront_Mdotzero_factor, + double DIM_front_Mdot_factor, const std::string &check_state_approach, const std::string &check_Sigma_approach, const std::string &check_Temp_approach, + const std::string &DIM_front_approach, + const std::string &scatter_by_corona, const std::string &initialcond, const std::shared_ptr initial_F_function, - const std::string &wind, const pard &windparams): - opacity(opacity), oprel(oprel), + const std::string &wind, + const pard &windparams): + opacity(opacity), + oprel(oprel), Mdotout(Mdotout), - boundcond(boundcond), Thot(Thot), Tirr2Tvishot(Tirr2Tvishot), - Rhot_Mdotzero_factor(Rhot_Mdotzero_factor), + boundcond(boundcond), + Thot(Thot), + Tirr2Tvishot(Tirr2Tvishot), + Rfront_Mdotzero_factor(Rfront_Mdotzero_factor), + DIM_front_Mdot_factor(DIM_front_Mdot_factor), check_state_approach(check_state_approach), check_Sigma_approach(check_Sigma_approach), check_Temp_approach(check_Temp_approach), + DIM_front_approach(DIM_front_approach), + scatter_by_corona(scatter_by_corona), initialcond(initialcond), initial_F_function(initial_F_function), - wind(wind), windparams(windparams) {} + wind(wind), + windparams(windparams) {} inline vecd initial_F(const vecd& h) const { return (*initial_F_function)(h); } inline size_t initial_first(const vecd& h) const { return initial_F_function->first(h); } }; @@ -328,6 +356,7 @@ class CalculationArguments { constexpr static const unsigned int default_Nt_for_tau = 200; constexpr static const char default_gridscale[] = "log"; constexpr static const unsigned short default_starlod = 3; + constexpr static const int default_verb_level = 0; public: double init_time; double time; @@ -335,18 +364,19 @@ class CalculationArguments { unsigned int Nx; std::string gridscale; unsigned short starlod = 3; + int verb_level; double eps; public: CalculationArguments( double inittime, double time, std::optional tau, - unsigned int Nx, const std::string& gridscale, const unsigned short starlod, - double eps=1e-6): + unsigned int Nx, const std::string& gridscale, const unsigned short starlod, int verb_level, + double eps=1e-6): init_time(inittime), time(time), tau(tau ? *tau : time / default_Nt_for_tau), Nx(Nx), gridscale(gridscale), starlod(starlod), - eps(eps) {} + verb_level(verb_level), eps(eps) {} }; diff --git a/cpp/include/exceptions.hpp b/cpp/include/exceptions.hpp index 9479cce..1799f2f 100644 --- a/cpp/include/exceptions.hpp +++ b/cpp/include/exceptions.hpp @@ -10,4 +10,11 @@ class RadiusCollapseException: public std::exception { } }; +class DiscEqFailException: public std::exception { +public: + virtual const char* what() const noexcept override { + return "Iterations >=2000"; + } +}; + #endif //FREDDI_EXCEPTIONS_HPP diff --git a/cpp/include/freddi_state.hpp b/cpp/include/freddi_state.hpp index 23dac60..41e3555 100644 --- a/cpp/include/freddi_state.hpp +++ b/cpp/include/freddi_state.hpp @@ -243,12 +243,16 @@ class FreddiState { public: double Mdot_out; double Mdot_in_prev = -INFINITY; + double Mdot_outer_boundary; double t; size_t i_t; size_t first; size_t last; vecd F; double F_in; + double R_dotM0_before_shift = 0.0; + double maxR_Qirr_no_role = 0.0; + //double DIM_front_Mdot_factor = 0.0; explicit CurrentState(const DiskStructure& str); CurrentState(const CurrentState&) = default; CurrentState& operator=(const CurrentState&) = default; @@ -261,7 +265,7 @@ class FreddiState { boost::optional Mdisk; boost::optional Lx; boost::optional Mdot_wind; - boost::optional W, Tph, Qx, Tph_vis, Tph_X, Tirr, Kirr, Sigma, Height; + boost::optional W, Tph, Qx, Tph_vis, Tph_X, Tirr, Kirr, Sigma, Height, Shadow; }; protected: @@ -304,6 +308,7 @@ class FreddiState { // current_ public: inline double Mdot_out() const { return current_.Mdot_out; } + inline double Mdot_outer_boundary() const { return current_.Mdot_outer_boundary; } inline double F_in() const { return current_.F_in; } inline const vecd& F() const { return current_.F; } inline double t() const { return current_.t; } @@ -311,15 +316,26 @@ class FreddiState { inline size_t first() const { return current_.first; } inline size_t last() const { return current_.last; } inline double Mdot_in_prev() const { return current_.Mdot_in_prev; } + inline double R_dotM0_before_shift() const { return current_.R_dotM0_before_shift; } + inline double maxR_Qirr_no_role() const { return current_.maxR_Qirr_no_role; } + //inline double DIM_front_Mdot_factor() const { return current_.DIM_front_Mdot_factor; } protected: inline void set_Mdot_in_prev(double Mdot_in) { current_.Mdot_in_prev = Mdot_in; } inline void set_Mdot_in_prev() { set_Mdot_in_prev(Mdot_in()); } + inline void set_Mdot_outer_boundary(double Mdot_boundary) { current_.Mdot_outer_boundary = Mdot_boundary; } + inline void set_R_dotM0_before_shift(double r) { current_.R_dotM0_before_shift = r; } + inline void set_maxR_Qirr_no_role(double r) {current_.maxR_Qirr_no_role = r;} + //inline void set_DIM_front_Mdot_factor(double r) { current_.DIM_front_Mdot_factor = r; } virtual IrradiatedStar::sources_t star_irr_sources(); public: inline double omega_R(double r) const { return std::sqrt(GM() / (r*r*r)); } inline double omega_i(size_t i) const { return omega_R(R()[i]); } virtual double Mdot_in() const; virtual double Lbol_disk() const; + virtual double Rfront_Rhot(double r, double z_r) const; + virtual double obtain_Mdot_outer_boundary(); + virtual double Tirr_critical(double r, int ii); + virtual void verify_disc_mass(double tau); double phase_opt() const; // wind_ public: @@ -397,6 +413,7 @@ class FreddiState { const vecd& Tirr(); const vecd& Kirr(); const vecd& Height(); + const vecd& Shadow(); double Luminosity(const vecd& T, double nu1, double nu2) const; inline double magnitude(const double lambda, const double F0) { return -2.5 * std::log10(I_lambda(lambda) * cosiOverD2() / F0); @@ -425,9 +442,11 @@ class FreddiState { double Sigma_minus(double r) const; double Sigma_plus(double r) const; double Teff_plus(double r) const; - double R_cooling_front(double r); - double v_cooling_front(double r); - int ring_state_vertical(const int i); // check ring hot (1) or cold (0) + double R_cooling_front(double r, double sigma_at_r); + double R_vis_struct (double r, double z2r_at_r); + double v_cooling_front(double r, double sigma_at_r); + double v_visc(double r, double z2r); + int check_ring_is_cold(const int i); // check ring hot (1) or cold (0) int Tirr_exceed_critical(const int i); // check Tirr > Tcrit (1) or not (0) }; diff --git a/cpp/include/nonlinear_diffusion.hpp b/cpp/include/nonlinear_diffusion.hpp index d40a5a0..d9be2b5 100644 --- a/cpp/include/nonlinear_diffusion.hpp +++ b/cpp/include/nonlinear_diffusion.hpp @@ -7,7 +7,7 @@ #include // std::exception #include // std::function #include - +#include typedef std::vector vecd; diff --git a/cpp/include/ns/ns_arguments.hpp b/cpp/include/ns/ns_arguments.hpp index b12a4e2..e3c9228 100644 --- a/cpp/include/ns/ns_arguments.hpp +++ b/cpp/include/ns/ns_arguments.hpp @@ -114,7 +114,8 @@ class NeutronStarDiskStructureArguments: public DiskStructureArguments { const std::string& opacity, double Mdotout, const std::string& boundcond, double Thot, double Tirr2Tvishot, - double Rhot_Mdotzero_factor, const std::string& check_state_approach, const std::string& check_Sigma_approach, const std::string& check_Temp_approach, + double Rfront_Mdotzero_factor, double DIM_front_Mdot_factor, const std::string& check_state_approach, const std::string& check_Sigma_approach, const std::string& check_Temp_approach, + const std::string& DIM_front_approach, const std::string& scatter_by_corona, const std::string& initialcond, std::optional F0, std::optional Mdisk0, std::optional Mdot0, diff --git a/cpp/include/orbit.hpp b/cpp/include/orbit.hpp index 05089d3..1897988 100644 --- a/cpp/include/orbit.hpp +++ b/cpp/include/orbit.hpp @@ -14,7 +14,7 @@ double efficiency_of_accretion(double kerr); class BlackHoleFunctions { public: - static inline double r_kerrISCO(const double Mx, const double kerr) { return rgToCm(r_kerrISCORg(kerr), Mx); } + static inline double r_kerrISCO(const double Mx, const double kerr) { return rgToCm(r_kerrISCORg(kerr), Mx); } }; diff --git a/cpp/src/arguments.cpp b/cpp/src/arguments.cpp index fdab2d1..bc9b265 100644 --- a/cpp/src/arguments.cpp +++ b/cpp/src/arguments.cpp @@ -24,6 +24,8 @@ constexpr const double DiskStructureArguments::default_Mdotout; constexpr const char DiskStructureArguments::default_initialcond[]; constexpr const double DiskStructureArguments::default_Thot; constexpr const char DiskStructureArguments::default_boundcond[]; +constexpr const char DiskStructureArguments::default_scatter_by_corona[]; +constexpr const char DiskStructureArguments::default_DIM_front_approach[]; constexpr const double DiskStructureArguments::mu; constexpr const char DiskStructureArguments::default_wind[]; @@ -32,8 +34,11 @@ DiskStructureArguments::DiskStructureArguments( const std::string& opacity, double Mdotout, const std::string& boundcond, double Thot, double Tirr2Tvishot, - double Rhot_Mdotzero_factor, const std::string& check_state_approach, const std::string& check_Sigma_approach, + double Rfront_Mdotzero_factor, double DIM_front_Mdot_factor, const std::string& check_state_approach, + const std::string& check_Sigma_approach, const std::string& check_Temp_approach, + const std::string& DIM_front_approach, + const std::string& scatter_by_corona, const std::string& initialcond, std::optional F0, std::optional Mdisk0, std::optional Mdot0, @@ -47,8 +52,8 @@ DiskStructureArguments::DiskStructureArguments( boundcond(boundcond), Thot(Thot), Tirr2Tvishot(Tirr2Tvishot), - Rhot_Mdotzero_factor(Rhot_Mdotzero_factor), - check_state_approach(check_state_approach),check_Sigma_approach(check_Sigma_approach), check_Temp_approach(check_Temp_approach), + Rfront_Mdotzero_factor(Rfront_Mdotzero_factor), DIM_front_Mdot_factor(DIM_front_Mdot_factor), check_state_approach(check_state_approach),check_Sigma_approach(check_Sigma_approach), check_Temp_approach(check_Temp_approach), + DIM_front_approach(DIM_front_approach),scatter_by_corona(scatter_by_corona), initialcond(initialcond), initial_F_function(initializeInitialFFunction(oprel, bdb_args, initialcond, @@ -376,3 +381,4 @@ constexpr const unsigned int CalculationArguments::default_Nx; constexpr const unsigned int CalculationArguments::default_Nt_for_tau; constexpr const char CalculationArguments::default_gridscale[]; constexpr const unsigned short CalculationArguments::default_starlod; +constexpr const int CalculationArguments::default_verb_level; diff --git a/cpp/src/freddi_evolution.cpp b/cpp/src/freddi_evolution.cpp index 7b9f930..684c7d6 100644 --- a/cpp/src/freddi_evolution.cpp +++ b/cpp/src/freddi_evolution.cpp @@ -7,6 +7,8 @@ #include "exceptions.hpp" #include "nonlinear_diffusion.hpp" +#define VERB_LEVEL_MESSAGES 30 + using namespace std::placeholders; @@ -15,104 +17,27 @@ FreddiEvolution::FreddiEvolution(const FreddiArguments &args): void FreddiEvolution::step(const double tau) { + //if (args().calc->verb_level > VERB_LEVEL_MESSAGES) {std::cout << "cA_ t="<< [freddi]() {return sToDay(current_.t);} <<"\n" << std::endl;} + if (args().calc->verb_level > VERB_LEVEL_MESSAGES) {std::cout << "c_A_ t="<< sToDay(current_.t) <<"\n" << std::endl;} truncateInnerRadius(); FreddiState::step(tau); + if (args().calc->verb_level > VERB_LEVEL_MESSAGES) {std::cout << "c_A__ t="<< sToDay(current_.t) <<"\n" << std::endl;} nonlinear_diffusion_nonuniform_wind_1_2( args().calc->tau, args().calc->eps, - F_in(), Mdot_out(), + F_in(), Mdot_outer_boundary(), windA(), windB(), windC(), wunc(), h(), current_.F, first(), last()); + if (args().calc->verb_level > VERB_LEVEL_MESSAGES) {std::cout << "c_A___ t="<< sToDay(current_.t) <<"\n" << std::endl;} truncateOuterRadius(); star_.set_sources(star_irr_sources()); } -int FreddiState::Tirr_exceed_critical(const int ii) { - // returns 1 for Tirr > Tcrit, 0 - otherwise - - double radius_popravka = args().disk->Rhot_Mdotzero_factor; - double Rcheck = std::min( R().at(ii-1) * radius_popravka, args().basic->rout); - double radius_popravka_variable = Rcheck/R().at(ii-1); - - radius_popravka = radius_popravka_variable; - - double Tcrit; - if (args().disk->check_Temp_approach == "const") { - - Tcrit = args().disk->Thot * pow(radius_popravka,0.5); - - } else if (args().disk->check_Temp_approach == "Tavleev") { - double Qirr_Qvis = pow(Tirr().at(ii) / Tph_vis().at(ii),4.)*radius_popravka ; - if (Qirr_Qvis > 1.) { - Tcrit = (9040. - 2216.* 1./Qirr_Qvis ) * pow(radius_popravka,0.5) ; - } else { - Tcrit = (9040. - 2216.) * pow(radius_popravka,0.5) ; - } - // DEBUG: - //if (radius_popravka > 1.01) { - // Tcrit = 10000. * pow(radius_popravka,0.5); - //} - } else { - throw std::invalid_argument("Wrong check_Temp_approach [const/Tavleev]"); - } - if ( Tirr().at(ii) >= Tcrit ) { - return 1; - } else { - return 0; - } - throw std::invalid_argument("Logical error at Tirr_exceed_critical"); - return 0; -} -int FreddiState::ring_state_vertical(const int ii) { - // returns 1 for hot, 0 for cold - - double radius_popravka = args().disk->Rhot_Mdotzero_factor; - - // R().at(ii) is the radius where Mdot = 0 - // multiplied by radius_popravka, gives the hot zone radius - // Sigma_minus is the maximum density on the cold branch - // Sigma_plus is the minimum density on the hot branch = Sigma_min(Menou+1999) - - if ( Tirr_exceed_critical(ii) ) { - // if irradiation temperature is greater than critical, disc cannot be cold - return 1; - } else { - // check if surface density at front is larger than critical cold Sigma_minus - // - if (Sigma().at(ii)/pow(radius_popravka, -0.75) > Sigma_minus(R().at(ii) * radius_popravka)) { - // disc cannot be cold if density is greater than critical for cold state - return 1; - } else { - // check if surface density at front is less than critical hot Sigma_plus - double sigma_factor; - if (args().disk->check_Sigma_approach == "Menou99a") { - sigma_factor = 4.3 ; // See fig.8 of Menou+1999; this correctly work only for constant radius_popravka - } else if (args().disk->check_Sigma_approach == "simple") { - sigma_factor = pow(radius_popravka, -0.75); - } else { - throw std::invalid_argument("Wrong check_Sigma_approach [Menou99a/simple]"); - } - if ( Sigma().at(ii)/sigma_factor < Sigma_plus(R().at(ii)*radius_popravka) ) { - //disc cannot be hot if density is lower than critical for hot state - return 0; - } else { - // check cooling front position - if (radius_popravka * R().at(ii) > R_cooling_front ( radius_popravka*R().at(ii)) ) { - //radius is beyond front - return 0; - } else { - return 1; - } - } - } - } - throw std::invalid_argument("ring_state_vertical: logic mistake"); - return 1; -} void FreddiEvolution::truncateOuterRadius() { + if (args().disk->Thot <= 0. ){ return; } @@ -124,90 +49,17 @@ void FreddiEvolution::truncateOuterRadius() { if (Mdot_in() > Mdot_in_prev()) { return; } - + auto ii = last() + 1; + if (args().calc->verb_level > VERB_LEVEL_MESSAGES) {std::cout <<"c_A last ii = "<< last() <<" last R=" <Rhot_Mdotzero_factor; //1.8; 2.1; 1.7; - if (Tirr().at(last()) / Tph_vis().at(last()) < args().disk->Tirr2Tvishot/pow(radius_popravka,0.25)) { - // when irradiation is not important - // (A) hot disc extends as far as Sigma>Sigma_max_cold(alpha_cold) and not farther than R_cooling_front and Tirr >= Thot and Teff_vis >=Teff_plus (condition below is the opposite) - // or (B) hot disk cannot exist for Sigma< Sigma_crit_hot (=Sigma+) while Tirr > Tcrit: this condition ensures fast - // convergence to small disc radius in the case without irradiation - //DEBUG std::cout << Tirr().at(last()) << "\n" << std::endl; - - - do { - ii--; - if (ii <= first()) throw RadiusCollapseException(); - } while( - // conditions for cold zone: - //(A) (1) hot disc extends not farther than R_cooling_front - // sent to DIM team - ((args().disk->check_state_approach == "before2024") && - ( - ( radius_popravka * R().at(ii) > R_cooling_front ( radius_popravka*R().at(ii)) ) - - // Sigma_minus is the maximum density on the cold branch - // Sigma_plus is the minimum density on the hot branch = Sigma_min(Menou+1999) - // (2) hot disc extends as far as Sigma>Sigma_max_cold(alpha_cold) - // && ( Sigma().at(ii)/pow(radius_popravka,0.75) < Sigma_minus(R().at(ii)) ) - // && ( Sigma().at(ii) < Sigma_minus(R().at(ii)) ) - // && ( Sigma().at(ii)/4.3 < Sigma_minus(R().at(ii)) * pow(radius_popravka,1.18) ) // error? - && ((args().disk->check_Sigma_approach == "Menou99a") && ( Sigma().at(ii)/4.3 < Sigma_plus(R().at(ii)) * pow(radius_popravka,1.11) )) // see Fig. 8 of Menou+99 - - // && ( Tirr().at(ii)/pow(radius_popravka,0.5) < args().disk->Thot ) - // && ( Tirr().at(ii)/pow(radius_popravka,0.5) < 9040. - 2216.* 1./(pow(Tirr().at(ii) / Tph_vis().at(ii),4.)*radius_popravka) ) // this is valid only if Tirr>Tvis - - // (3) in the cold disc temperature < critical - && ( Tirr().at(ii)/pow(radius_popravka,0.5) < args().disk->Thot) //6800 was - )) - || - ( (args().disk->check_state_approach == "logic") && ( ring_state_vertical(ii) == 0) ) - - // (B) - // || ( - // ( Sigma().at(ii) < Sigma_plus(R().at(ii)) ) - // && ( Tirr().at(ii) < args().disk->Thot ) - //) - - // (3) next line is added according to Tavleev+23 results: it overrides the line above if Thot=10000. In fact, the line above should be deleted. - // && ( Tirr().at(ii) < 9040. - 2216. ) - // (4 - possibly wrong) next line is added to prevent cooling wave when effective viscous temperature is too high: - //&& ( Tph_vis().at(ii) < Teff_plus(R().at(ii)) ) - ); - //} while( ( R().at(ii) > R_cooling_front ( R().at(ii)) ) && ( Sigma().at(ii) < Sigma_minus(R().at(ii)) ) ); - } else if (args().disk->boundcond == "Teff") { - // irradiation is important, the boundary is at fixed Teff - do { - ii--; - if (ii <= first()) throw RadiusCollapseException(); - } while( Tph().at(ii) < args().disk->Thot ); - - - } else if (args().disk->boundcond == "Tirr") { - -// double Rcheck; -// Rcheck = std::min( R().at(ii-1)*radius_popravka, args().basic->rout); -// double radius_popravka_variable = Rcheck/R().at(ii-1); - // DEBUG std::cout << radius_popravka_variable << " " << 9040. - 2216.* 1./(pow(Tirr().at(ii-1) / Tph_vis().at(ii-1),4.)*radius_popravka_variable) << "\n" << std::endl; - // irradiation is important, the boundary is at fixed Tir - do { - ii--; - if (ii <= first()) throw RadiusCollapseException(); - //} while( Tirr().at(ii) < 9040. - 2216.* Tph_vis().at(ii) / Tirr().at(ii) ); - // according to Tavleev+23 results; it is applicable in quasi-stationary disc (without cooling wave) -// std::cout << Tirr().at(ii) << "\n" << std::endl; -// std::cout << args().disk->Thot << "\n" << std::endl; - // old: -// } while( Tirr().at(ii)/pow(radius_popravka,0.5) < args().disk->Thot); - // } while( Tirr().at(ii)/pow(radius_popravka_variable,0.5) < args().disk->Thot); - } while (Tirr_exceed_critical(ii) == 0); -// } while( Tirr().at(ii)/pow(radius_popravka_variable,0.5) < 9040. - 2216.* 1./(pow(Tirr().at(ii) / Tph_vis().at(ii),4.)*radius_popravka_variable) ); //variable Tcrit does not correspond to SIM model - } else{ - throw std::invalid_argument("Wrong boundcond"); - } + do { + ii--; + if (ii <= first()) throw RadiusCollapseException(); + } while ( check_ring_is_cold(ii) ); + if ( ii <= last() - 1 ){ current_.last = ii; @@ -223,3 +75,6 @@ vecd FreddiEvolution::wunction(const vecd &h, const vecd &F, size_t _first, size } return W; }; + + +#undef VERB_LEVEL_MESSAGES diff --git a/cpp/src/freddi_state.cpp b/cpp/src/freddi_state.cpp index 37b44ff..194896b 100644 --- a/cpp/src/freddi_state.cpp +++ b/cpp/src/freddi_state.cpp @@ -8,6 +8,7 @@ #include "arguments.hpp" #include "nonlinear_diffusion.hpp" #include "orbit.hpp" +#define VERB_LEVEL_MESSAGES 30 FreddiState::DiskStructure::DiskStructure(const FreddiArguments &args, const wunc_t& wunc): @@ -146,13 +147,62 @@ void FreddiState::replaceArgs(const FreddiArguments &args) { void FreddiState::step(double tau) { set_Mdot_in_prev(); + set_Mdot_outer_boundary(obtain_Mdot_outer_boundary()); + verify_disc_mass(tau); invalidate_optional_structure(); + current_.i_t ++; current_.t += tau; wind_->update(*this); } +void FreddiState::verify_disc_mass (double tau) { + if (args().calc->verb_level > VERB_LEVEL_MESSAGES-5) {std::cout << "+++ Mdisk("<< current_.t<<") = " << Mdisk() << " Rlast+1=" << R()[last()+1] << " delta_M_shift= " <<2.0 * M_PI * R()[last()] * Sigma()[last()] *(R()[last()+1]-R()[last()]) << "\n" << std::endl; } + if (args().calc->verb_level > VERB_LEVEL_MESSAGES-5) {std::cout << "+++2 Mdot_in=" << Mdot_in() << " Mdot_out=" << obtain_Mdot_outer_boundary()<<" Mdisk+Mdot*tau=" << Mdisk()+(Mdot_in() - obtain_Mdot_outer_boundary()) * tau << "\n" << std::endl; } + if ( (args().calc->verb_level > VERB_LEVEL_MESSAGES-5) && (R()[last()+1]>0)) {std::cout << "+++ Mdisk+Mdot*tau+Sigma*2*Pi*R*dR =" << Mdisk()+(Mdot_in() + Mdot_in_prev()*2.3) * tau + 2.0 * M_PI * R()[last()] * Sigma()[last()] *(1.-3./4. * (R()[last()+1]-R()[last()])/(R()[last()+1]+R()[last()]) ) * (R()[last()+1]-R()[last()]) << "\n+++" << (1.-3./4. * (R()[last()+1]-R()[last()])/(R()[last()+1]+R()[last()]) ) * (R()[last()+1]-R()[last()]) << std::endl; } +} + + + +double FreddiState::obtain_Mdot_outer_boundary() { + + if (args().disk->DIM_front_approach == "outflow") { + //if (args().disk->boundcond == "Hameury") { + if (args().disk->scatter_by_corona == "yes") { + // assuming there is scatter above the disc + if (last() == Nx()-1) { + return Mdot_out(); + } else { + return -1.0*args().disk->DIM_front_Mdot_factor*Mdot_in(); + //-1.5 make Rfront/RMdot0 = 2 + // -2.5; -2.3 makes very close to Hameury with =Tcrit = (8840. - 2216.* 1./Qirr_Qvis ) * pow(radius_popravka,0.5) ; and Tcrit = 10300.; + } + } + //if (args().disk->boundcond == "Hameury_no_scatter") { + if (args().disk->scatter_by_corona == "no") { + // if there is no scattering, the disc zone beyond maximum of Fvis (where dotM=0) is in the shadow + // This means its temperature + // if irradiation temperature is greater than critical, disc cannot be cold + // Tirr is checked at Thot or Tfront, depending on option boundcond (scatter/no scatter),see Tirr_critical. + // If there is no scattering, the disc beyond r, where dotM=0, is in the shadow from the direct central radiation + std::cout << "T irr_last=" << Tirr().at(last()) << std::endl; + if ((last() == Nx()-1) or ( Tirr().at(last()) >= Tirr_critical (R().at(last()), last()) ) ) { + return Mdot_out(); + } else { + // std::cout << "CF!" << std::endl; + return -1.0*args().disk->DIM_front_Mdot_factor*Mdot_in(); + //-1.5 make Rfront/RMdot0 = 2 + // -2.5; -2.3 makes very close to Hameury with =Tcrit = (8840. - 2216.* 1./Qirr_Qvis ) * pow(radius_popravka,0.5) ; and Tcrit = 10300.; + } + } + + } + // DIM_front_approach == "maxFvis" : + return Mdot_out(); + +} + double FreddiState::Mdot_in() const { return (F()[first() + 1] - F()[first()]) / (h()[first() + 1] - h()[first()]); } @@ -163,6 +213,7 @@ double FreddiState::Lbol_disk() const { } + double FreddiState::phase_opt() const { return 2.0 * M_PI * (t() - args().flux->ephemeris_t0) / args().basic->period; } @@ -228,13 +279,15 @@ const vecd& FreddiState::Tirr() { const vecd& FreddiState::Qx() { + if (!opt_str_.Qx) { vecd x(Nx()); const vecd& K = Kirr(); const vecd& H = Height(); + const vecd& Shad = Shadow(); const double Lbol = Lbol_disk(); for (size_t i = first(); i < Nx(); i++) { - x[i] = K[i] * Lbol * angular_dist_disk(H[i] / R()[i]) / (4. * M_PI * m::pow<2>(R()[i])); + x[i] = (1.0 - Shad[i]) * K[i] * Lbol * angular_dist_disk(H[i] / R()[i]) / (4. * M_PI * m::pow<2>(R()[i])); } opt_str_.Qx = std::move(x); } @@ -273,6 +326,28 @@ const vecd& FreddiState::Height() { } +const vecd& FreddiState::Shadow() { + vecd x(Nx()); + double max_H2R = 0.0; + for (size_t i = first(); i < Nx(); i++) { + const double Hrelative= Height()[i]/R()[i]; + max_H2R = std::max(Hrelative, max_H2R); + if (Hrelative >= max_H2R) { + x[i] = 0.0; + } else { + if (args().disk->scatter_by_corona == "no") { + x[i] = 1.0; + } + } + } + std::cout << "max= " << max_H2R << "H2R last = " << Height()[last()]/R()[last()] << "last = " << last() << std::endl; + opt_str_.Shadow = std::move(x); + return *opt_str_.Shadow; +} + + + + const vecd& FreddiState::Tph_vis() { if (!opt_str_.Tph_vis) { vecd x(Nx(), 0.0); @@ -701,30 +776,402 @@ double FreddiState::Sigma_plus(double r) const { //@@@@@ 28/09/2023 double FreddiState::Teff_plus(double r) const { - // Lasota et al., A&A 486, 523–528 (2008), Eq A.1, DOI: 10.1051/0004-6361:200809658 - return 6890 * std::pow(r / 1e10, -0.09) - * std::pow(args().basic->Mx / GSL_CONST_CGSM_SOLAR_MASS, 0.03); + // minimum temperature on the hot branch; Tavleev+23: + double A = 7341.; + double beta = 0.0290; + double gamma = -0.00484; + double delta = -0.08426; + return A * std::pow(args().basic->Mx / GSL_CONST_CGSM_SOLAR_MASS, beta) * std::pow(args().basic->alpha , gamma) * std::pow(r / 1e10, delta); + + // Lasota et al., A&A 486, 523–528 (2008), Eq A.1, DOI: 10.1051/0004-6361:200809658 + return 6890 * std::pow(r / 1e10, -0.09) + * std::pow(args().basic->Mx / GSL_CONST_CGSM_SOLAR_MASS, 0.03); } +// double FreddiState::Teff_minus(double r) const { +// // maximum temperature on the cold branch Tavleev+23: +// double A = 6152.; +// double beta = 0.0315; +// double gamma = 0.00165; +// double delta = -0.08977; +// return A * std::pow(args().basic->Mx / GSL_CONST_CGSM_SOLAR_MASS, beta) * std::pow(args().basic->alpha , gamma) * std::pow(r / 1e10, delta); +// +// } + +double FreddiState::v_visc(double r, double z2r_at_r) { + //oprel().Height(R()[i], F()[i]) + return args().basic->alpha/oprel().Pi3 * m::pow<2>(z2r_at_r) * std::sqrt(GM()/r); + //Height2R()[] +} -double FreddiState::v_cooling_front(double r) { +double FreddiState::R_vis_struct (double r, double z2r_at_r) { + // estimate the radius, to which viscous process can reach in time tau + return r - v_visc(r, z2r_at_r) * args().calc->tau; +} + + +double FreddiState::v_cooling_front(double r, double sigma_at_r) { // The cooling-front velocity depends on the ratio between the current Sigma and critical Sigmas // Ludwig et al., A&A 290, 473-486 (1994), section 3 + // https://ui.adsabs.harvard.edu/abs/1994A%26A...290..473L // units: cm/s const double Sigma_plus_ = Sigma_plus(r); - const double sigma = std::log( Sigma()[last()] / Sigma_plus_ ) / std::log( Sigma_minus(r)/Sigma_plus_ ) ; - return 1e5 * (1.439-5.305*sigma+10.440*m::pow<2>(sigma)-10.55*m::pow<3>(sigma)+4.142*m::pow<4>(sigma)) + + //double sigma = std::log( Sigma()[last()] / Sigma_plus_ ) / std::log( Sigma_minus(r)/Sigma_plus_ ) ; + //sigma = 0; + //Sigma_plus is the minimum density on the hot branch + + double sigma = std::log( sigma_at_r / Sigma_plus_ ) / std::log( Sigma_minus(r)/Sigma_plus_ ) ; + if (sigma<0.) {sigma=0.0;} + + double v = 1e5 * (1.439-5.305*sigma+10.440*m::pow<2>(sigma)-10.55*m::pow<3>(sigma)+4.142*m::pow<4>(sigma)) * std::pow(args().basic->alpha / 0.2, 0.85-0.69*sigma) * std::pow(args().basic->alphacold / 0.05, 0.05+0.69*sigma) * std::pow(r / 1e10, 0.035) * std::pow(args().basic->Mx / GSL_CONST_CGSM_SOLAR_MASS, -0.012); + if (args().calc->verb_level > VERB_LEVEL_MESSAGES) {std::cout << "c_R v = " << v << " r=" << r << " Sigma="<< sigma_at_r << " sigma=" << sigma<< "\n" << std::endl; } + + return v; } -double FreddiState::R_cooling_front(double r) { +double FreddiState::R_cooling_front (double r, double sigma_at_r) { + // estimate the radius, where cooling front could reach, starting from r, if it moved with corresponding velocity // previous location of Rhot moves with the cooling-front velocity: - return R()[last()] - v_cooling_front(r) * args().calc->tau; + return r - v_cooling_front(r, sigma_at_r) * args().calc->tau; + + + return R()[last()] - v_cooling_front(r,sigma_at_r) * args().calc->tau; //return R()[last()] - v_cooling_front(R()[last()]) * args().calc->tau ; // this variant leads to more abrupt evolution, since the front velocity is larger } + +double FreddiState::Rfront_Rhot (double r, double z_r) const { + + if (args().disk->Rfront_Mdotzero_factor == 1.0) {return 1.0;} + // call: Rfront_Rhot(R().at(ii-1)) or Rfront_Rhot( R().at(ii) ) + double Rfront_Rhot=1.0; + double Rfront_Rhot_variable = std::min(args().disk->Rfront_Mdotzero_factor, args().basic->rout/r); + + /* if (Rfront_Rhot_variable < args().disk->Rfront_Mdotzero_factor) { + if (current_.R_dotM0_before_shift == 0.) { + set_R_dotM0_before_shift(r); + } + } */ + /* + // at start: (1) R_dotM0 = Rout or (2) Rhot_prev < Rout, (2C) if 2 Rhot_prev > Rout + // (1-) no scatter -> R_dotM0(t) moves with v_visc; +// (a) Qirr Yes : Rfront = R_dotM0(t) +// (b) Qirr No : * Rfront increases from R_fotM0(t) until factor = 2 + // (1+) scatter -> +// (a) Qirr Yes : * Rfront increases from R_fotM0(t) until factor = 2 +// (b) Qirr No : * Rfront increases from R_fotM0(t) until factor = 2 +// (2C) At the beginning disc should collapse: +// (2C-) no scatter : +// (a) Qirr Yes : Rhot = R_dotM0 +// (b) Qirr No : Rhot <= 2 R_dotM0 at t=0 and optionally still moves +// (2C+) scatter: +// Qirr Yes : * Rhot <= 2 R_dotM0 at t=0 and optionally still moves +// Qirr No : * Rhot <= 2 R_dotM0 at t=0 and optionally still moves +// (2-) no scatter +// Qirr Yes : +// Qirr No : +// (2+) scatter +// Qirr Yes : +// Qirr No : + // Fmaximum cannot move faster than v_visc + // notice initial position of Fmaximum + return Rfront_Rhot_variable; + + } + if (current_.R_dotM0_before_shift > 0.) && (r > current_.R_dotM0_before_shift/args().disk->Rfront_Mdotzero_factor) { + // Fmaximum moves until it shifts so that Rfront_Rhot_variable = 2: + double R_dotM0 = R_vis_struct(current_.R_dotM0_before_shift,z_r) + + Rfront_Rhot_variable = std::max (args().disk->Rfront_Mdotzero_factor, current_.R_dotM0_before_shift/R_dotM0); + } + } + */ + //if (args().disk->boundcond == "scatter_by_corona") { + if (args().disk->scatter_by_corona == "yes") { + Rfront_Rhot = Rfront_Rhot_variable; + } + //if (args().disk->boundcond == "no_scatter_by_corona") { + if (args().disk->scatter_by_corona == "no") { + if (current_.maxR_Qirr_no_role == 0.) { + // no scattering and irradiation's role is significant + Rfront_Rhot = 1.; // the disk beyond dot M = 0 is shadowed + } else { + // no scattering and irradiation's role is NOT significant + // then conditions of the hot zone are checked at radius > r (dotM=0) + Rfront_Rhot_variable = std::min( args().disk->Rfront_Mdotzero_factor, current_.maxR_Qirr_no_role / r); + Rfront_Rhot = Rfront_Rhot_variable; + //Rfront_Rhot=1.; //?????????????????? + + // здесь надо написать так, чтобы + /* if (Rfront_Rhot_variable < args().disk->Rfront_Mdotzero_factor) { + if (current_.R_dotM0_before_shift == 0.) { + set_R_dotM0_before_shift(r); + // Fmaximum cannot move faster than v_visc + // notice initial position of Fmaximum + return Rfront_Rhot_variable; + + } + if ((current_.R_dotM0_before_shift > 0.) && (r > current_.R_dotM0_before_shift/args().disk->Rfront_Mdotzero_factor)) { + // Fmaximum moves until it shifts so that Rfront_Rhot_variable = 2: + double R_dotM0 = R_vis_struct(current_.R_dotM0_before_shift,z_r) + return std::max (args().disk->Rfront_Mdotzero_factor, current_.R_dotM0_before_shift/R_dotM0); + } + } + */ + } + } + return Rfront_Rhot; +} + + +double FreddiState::Tirr_critical (double r, int ii) { + // determine critical level Tirr which keeps the ring hot. + // if we analyse r=Rfront, we should take into account that Rfront + // is farther than r == R().at(ii) + // + // it might be possible that function can take no arguments, if it calculates only at R()[last()] and use Tirr().at(last()) / Tph_vis().at(last()) + // newest approach: radius_popravka = 1 + + if (args().disk->boundcond == "Teff") { + // never gets here + throw std::invalid_argument("Why do you ask critical irradiation temperature if you set boundcond == Teff? Choose boundcond among Tirr/no_scatter_by_corona/scatter_by_corona\n"); + } + + // henceforth, only two possibilities remain: boundcond = "Tirr" or "no_scatter_by_corona" + + double radius_popravka = Rfront_Rhot( r, oprel().Height(R()[ii], F()[ii])/r); + + if (args().calc->verb_level > VERB_LEVEL_MESSAGES) {std::cout << " c_E R_last=" << R()[last()] << " R_ii="<< R()[ii] << " R_dotM0_before_shift=" << current_.R_dotM0_before_shift << " pop ="<< radius_popravka <<"\n" << std::endl;} + + + // Qirr/Qvis increases linearly with radius, if Qvis ~ r^(-3/4) which is not exactly true beyond radius where dotM=0, because Fvis is not proportional to h there + //thus, generally, Qirr_Qvis is the lower estimate + double Qirr_Qvis = pow(Tirr().at(ii) / Tph_vis().at(ii),4.)*radius_popravka ; + + // determine radius, inside which irradiation is not significant: + if ( (Qirr_Qvis < pow(args().disk->Tirr2Tvishot,4.) ) && (current_.maxR_Qirr_no_role == 0.) ) { + set_maxR_Qirr_no_role (r); + } + + //if ((obtain_Mdot_outer_boundary() < 0.0 ) && (args().disk->scatter_by_corona == "no")){ + // Shadow()[ii] = 1.0; + //} + + + double Tcrit; // value to return + + if (args().disk->check_Temp_approach == "const") { + // since we calculate critical level for irradiation temperature, we use the fact that Tirr ~ R(-1/2) (This is abs true only if Cirr = const) + Tcrit = args().disk->Thot * pow(radius_popravka,0.5); + +// if (args().disk->boundcond == "no_scatter_by_corona") { + if (args().disk->scatter_by_corona == "n_o") { + // do not take into account Tirr if Qirr/Qvis< critical_value + // and override previous assignment; instant return + if (Qirr_Qvis < pow(args().disk->Tirr2Tvishot,4.)) { + if (args().calc->verb_level > VERB_LEVEL_MESSAGES) {std::cout << "c_Y Qirr/Qvisverb_level > VERB_LEVEL_MESSAGES) {std::cout << "c_W Qirr/Qvis=\n" <check_Temp_approach == "Tavleev") || (args().disk->check_Temp_approach == "Hameury") ) { + + //if (args().disk->boundcond == "no_scatter_by_corona") { + if (args().disk->scatter_by_corona == "n_o") { + // if there is no scattering, the Rfront is in the shadow + // do not take into account Tirr if Qirr/Qvis < critical_value: + if (Qirr_Qvis < pow(args().disk->Tirr2Tvishot,4.)) { + if (args().calc->verb_level > VERB_LEVEL_MESSAGES) {std::cout << "c_V Qirr/Qvis=\n" < 1.) { + // since we calculate critical level for irradiation temperature, we use the fact that Tirr ~ R(-1/2) (This is completely true only if Cirr = const) + Tcrit = (9040. - 2216.* 1./Qirr_Qvis ) * pow(radius_popravka,0.5) ; + if (args().disk->check_Temp_approach == "Hameury") { + Tcrit = (8640. - 2216.* 1./Qirr_Qvis ) * pow(radius_popravka,0.5) ; + } + if (args().calc->verb_level > VERB_LEVEL_MESSAGES) {std::cout << "c_F Qirr/Qvis>1 with Tcrit(Rhot)=" << Tcrit << " ii="<verb_level > VERB_LEVEL_MESSAGES) {std::cout << "c_G Qirr/Qvis<1 with Tcrit(Rhot)=" << Tcrit << " ii="<check_Temp_approach == "Hameury") { + if (ii < Nx()-1) { + // + // if Rhot is less than Rtid, then critical temperature is fixed: + // this behaviour is established for solution to agree with Hameury'solution + Tcrit = 10300.; // Hameury email + } + } + } else { + throw std::invalid_argument("Wrong check_Temp_approac: Choose const or Tavleev or Hameury\n"); + } + if (args().calc->verb_level > VERB_LEVEL_MESSAGES) {std::cout << "c_H pop ="<< radius_popravka << " Tirr(Rhot) = " << Tirr().at(ii) << " and Tcrit=" << Tcrit<< "\n" << std::endl;} + + return Tcrit; + +} + +//int FreddiState::ring_state_vertical(const int ii) { +int FreddiState::check_ring_is_cold(const int ii) { + // returns 1 for hot, 0 for cold + + // R().at(ii), the radius where Mdot = 0, or ? + // multiplied by radius_popravka, gives the hot zone radius + // Sigma_minus is the maximum density on the cold branch + // Sigma_plus is the minimum density on the hot branch = Sigma_min(Menou+1999) + + + if (args().disk->boundcond == "Teff") { + if (Tph().at(ii) >= args().disk->Thot) { + return 0; // HOT + } else { + return 1; // COLD + } + } + + if ((args().disk->boundcond == "Tirr") || (args().disk->boundcond == "Tirr_Ham??eury")) { + if (args().calc->verb_level > VERB_LEVEL_MESSAGES) {std::cout << "boundcond=Tirr="<= Tirr_critical (R().at(ii), ii)) { + set_R_dotM0_before_shift( R().at(ii) ); + return 0; // HOT + } else { + return 1; // COLD + + } + } + + // only no_scatter_by_corona is possible ; everything else should be removed here: + if (!(args().disk->boundcond == "DIM")) { + throw std::invalid_argument("Wrong boundcond at check_ring_is_cold() in freddi_state"); + } + + + if ( Tirr().at(ii) >= Tirr_critical (R().at(ii), ii) ) { + // if irradiation temperature is greater than critical, disc cannot be cold + // Tirr is checked at Thot or Tfront, depending on option boundcond (scatter/no scatter),see Tirr_critical. + // If there is no scattering, the disc beyond r, where dotM=0, is in the shadow from the direct central radiation + // In practice, temperature of the disc at R, where dotM=0, is checked against the critical temperature multiplied by some factor + if (args().calc->verb_level > VERB_LEVEL_MESSAGES) {std::cout << "c_K Tirr high - HOT \n" << std::endl;} + // memorise radius before shift of the hot boundary + set_R_dotM0_before_shift( R().at(ii) ); + return 0; // HOT + } else { + // check_Sigma_approach -> "cold-front-approach" ? + // Menou99a -> Sigma_Menou99 ? + // if (args().disk->check_Sigma_approach == "Menou99a") { ? + // when irradiation stops preventing hot-zone shrinking, + // then conditions at Rfront matter: + // newest approach: radius_popravka = 1 + + double radius_popravka = Rfront_Rhot( R().at(ii), oprel().Height(R()[ii], F()[ii])/R().at(ii)); + + + // calculate sigma_popravka: + double sigma_factor; + if (radius_popravka == 1. ) { + // check the conditions at r where Mdot=0 + sigma_factor = 1.; + } else { + if (args().disk->check_Sigma_approach == "Menou99a") { + sigma_factor = 1./4.3 ; // See fig.8 of Menou+1999; this correctly work only for constant radius_popravka + } else if (args().disk->check_Sigma_approach == "simple") { + sigma_factor = pow(radius_popravka, -0.75); + } else if (args().disk->check_Sigma_approach == "critical_Teff") { + + } else { + throw std::invalid_argument("Wrong check_Sigma_approach [Menou99a/simple/critical_Teff]"); + } + } + + // check if surface density at front is larger than critical cold Sigma_minus + // + if (Sigma().at(ii)/pow(radius_popravka, -0.75) > Sigma_minus(R().at(ii) * radius_popravka)) { + // disc cannot be cold if density is greater than critical for cold state + if (args().calc->verb_level > VERB_LEVEL_MESSAGES) {std::cout << "c_L pop ="<< radius_popravka << " Sigma is high - HOT \n" << std::endl;} + set_R_dotM0_before_shift( R().at(ii) ); + return 0; // HOT + } else { + + if (args().calc->verb_level > VERB_LEVEL_MESSAGES) {std::cout << "c_S Sigma is LOW -> cooling wave OR collapse ii="<Nx-1 <<"\n" << std::endl;} + // last = Nx-1 means that code collapses to a starting Rhot : all outer disc is cast COLD + if ( (last() == args().calc->Nx-1) || (current_.R_dotM0_before_shift == 0.) ) { + if (args().calc->verb_level > VERB_LEVEL_MESSAGES) {std::cout << "c_O pop ="<< radius_popravka << " Sigma is low and R_dotM0_before_shift=" << current_.R_dotM0_before_shift<< "- COLD \n" << std::endl;} + return 1; // COLD + } else { + //throw std::invalid_argument("check_ring_is_cold: logic mistake 1"); + } + +// if ( last() < args().calc->Nx-1) { +// if (Sigma_plus(R().at(ii+1))== 0.) { Sigma().at(ii) +// // Sigma is low BUT ring is HOT +// // that means that the front is propagating and it is propagating with +// // meaningful velocity. +// // otherwise we search for a starting Rhot which position can be anything +// if (args().calc->verb_level > VERB_LEVEL_MESSAGES) {std::cout << "c_N pop ="<< radius_popravka << " Sigma is low BUT HOT " << Sigma_plus(R().at(ii+1)) << "\n"<< std::endl;} +// return 0; +// } +// } +// r????eturn 1; + } +// else { + + // R_cooling_front = r - v_cooling_front(r, sigma_at_r) * args().calc->tau; +// if (radius_popravka <= args().disk->Rfront_Mdotzero_factor) { +// if (args().calc->verb_level > VERB_LEVEL_MESSAGES) {std::cout << "c_PP pop ="<< radius_popravka << " radius_popravka < Rfront_Mdotzero_factor->COLD \n" << std::endl;} +// return 1; +// } + + // check cooling front position; cooling front is moving from last Rout + + if (current_.R_dotM0_before_shift == 0 ) { + // COLD ; just collapse to initial position: + return 1; + } + if ((radius_popravka * R().at(ii) > R_cooling_front ( Rfront_Rhot(R()[last()], oprel().Height(R()[last()], F()[last()])/R()[last()] ) * R()[last()] , Sigma().at(last())*sigma_factor) ) && ( last() < args().calc->Nx-1) ){ + //radius is beyond front + if (args().calc->verb_level > VERB_LEVEL_MESSAGES) {std::cout << "++c_P pop ="<< radius_popravka << " beyond front - COLD \n" << std::endl;} + // COLD + return 1; + } else { + if (args().calc->verb_level > VERB_LEVEL_MESSAGES) {std::cout << "++c_Q pop ="<< radius_popravka << " inwards front - HOT \n" << std::endl;} + // in fact front starts from a greater distance? + set_R_dotM0_before_shift( R().at(ii) ); + // HOT + return 0; + } +// } + } + } + throw std::invalid_argument("check_ring_is_cold: logic mistake 2"); + return 0; +} + +#undef VERB_LEVEL_MESSAGES + + diff --git a/cpp/src/nonlinear_diffusion.cpp b/cpp/src/nonlinear_diffusion.cpp index d3419c5..0b8b8be 100644 --- a/cpp/src/nonlinear_diffusion.cpp +++ b/cpp/src/nonlinear_diffusion.cpp @@ -1,5 +1,5 @@ #include "nonlinear_diffusion.hpp" - +#include "exceptions.hpp" // gvlipunova added for DiscEqFailException, did not know where else.... double mean_square_rel(const vecd &A, const vecd &B, size_t first, size_t last){ double rv = 0; @@ -66,7 +66,10 @@ void nonlinear_diffusion_nonuniform_wind_1_2 ( vecd alpha(last + 1), beta(last + 1); double c; + int iter_sol = 0; + int maxiter = 100; do { + iter_sol++; K_0 = K_1; alpha[first + 1] = 0.; beta[first + 1] = left_bounder_cond; @@ -85,5 +88,9 @@ void nonlinear_diffusion_nonuniform_wind_1_2 ( for (size_t i = 1; i <= last - 1; ++i) { K_1[i] = frac[i] * W[i] / y[i]; } - } while (max_dif_rel(K_1, K_0, 1, last - 1) > eps); + } while ((max_dif_rel(K_1, K_0, 1, last - 1) > eps) && (iter_sol <=maxiter)); + if (iter_sol >= maxiter) { + throw std::invalid_argument("Disc equation failed to converge."); + throw DiscEqFailException(); + } } diff --git a/cpp/src/ns/ns_arguments.cpp b/cpp/src/ns/ns_arguments.cpp index 24707e9..91ff0cb 100644 --- a/cpp/src/ns/ns_arguments.cpp +++ b/cpp/src/ns/ns_arguments.cpp @@ -115,10 +115,13 @@ NeutronStarDiskStructureArguments::NeutronStarDiskStructureArguments( const std::string& opacity, double Mdotout, const std::string& boundcond, double Thot, double Tirr2Tvishot, - double Rhot_Mdotzero_factor, + double Rfront_Mdotzero_factor, + double DIM_front_Mdot_factor, const std::string& check_state_approach, const std::string& check_Sigma_approach, const std::string& check_Temp_approach, + const std::string& DIM_front_approach, + const std::string& scatter_by_corona, const std::string& initialcond, std::optional F0, std::optional Mdisk0, std::optional Mdot0, @@ -128,9 +131,10 @@ NeutronStarDiskStructureArguments::NeutronStarDiskStructureArguments( ): DiskStructureArguments(opacity, OpacityRelated(opacity, bdb_args.Mx, bdb_args.alpha, mu), Mdotout, boundcond, Thot, Tirr2Tvishot, - Rhot_Mdotzero_factor, + Rfront_Mdotzero_factor, + DIM_front_Mdot_factor, check_state_approach, check_Sigma_approach, - check_Temp_approach, + check_Temp_approach, DIM_front_approach, scatter_by_corona, initialcond, initializeInitialFFunctionNS(OpacityRelated(opacity, bdb_args.Mx, bdb_args.alpha, mu), ns_args, bdb_args, diff --git a/cpp/src/ns/ns_evolution.cpp b/cpp/src/ns/ns_evolution.cpp index 0a50107..089d8dc 100644 --- a/cpp/src/ns/ns_evolution.cpp +++ b/cpp/src/ns/ns_evolution.cpp @@ -60,7 +60,8 @@ FreddiNeutronStarEvolution::NeutronStarStructure::NeutronStarStructure( dFmagn_dh(initialize_dFmagn_dh(evolution)), d2Fmagn_dh2(initialize_d2Fmagn_dh2(evolution)) { if (args_ns.Rdead > 0. && args_ns.Rdead < R_cor) { - throw std::logic_error("R_dead is positive and less than R_cor, it is unacceptably"); + throw std::logic_error("R_dead is positive and less than R_cor, it is unacceptable"); + std::cout << "R_dead="<< args_ns.Rdead << " R_cor=" << R_cor << "\n" << std::endl; } } @@ -492,11 +493,12 @@ const vecd& FreddiNeutronStarEvolution::Qx() { vecd x(Nx()); const vecd& K = Kirr(); const vecd& H = Height(); + const vecd& Shad = Shadow(); const double L_disk = Lbol_disk(); const double L_ns = Lbol_ns(); for (size_t i = first(); i < Nx(); i++) { const double mu = H[i] / R()[i]; - x[i] = K[i] * (L_disk * angular_dist_disk(mu) + L_ns * angular_dist_ns(mu)) / (4. * M_PI * m::pow<2>(R()[i])); + x[i] = (1.0 - Shad[i]) * K[i] * (L_disk * angular_dist_disk(mu) + L_ns * angular_dist_ns(mu)) / (4. * M_PI * m::pow<2>(R()[i])); } opt_str_.Qx = std::move(x); } diff --git a/cpp/src/ns/ns_options.cpp b/cpp/src/ns/ns_options.cpp index 78f93df..69e3f56 100644 --- a/cpp/src/ns/ns_options.cpp +++ b/cpp/src/ns/ns_options.cpp @@ -186,9 +186,12 @@ NeutronStarDiskStructureOptions::NeutronStarDiskStructureOptions(const po::varia vm["boundcond"].as(), vm["Thot"].as(), std::pow(vm["Qirr2Qvishot"].as(), 0.25), - vm["Rhot_Mdotzero_factor"].as(), + vm["Rfront_Mdotzero_factor"].as(), + vm["DIM_front_Mdot_factor"].as(), vm["check_state_approach"].as(), - vm["check_Sigma_approach"].as(), vm["check_Temp_approach"].as(), + vm["check_Sigma_approach"].as(), vm["check_Temp_approach"].as(), + vm["DIM_front_approach"].as(), + vm["scatter_by_corona"].as(), vm["initialcond"].as(), varToOpt(vm, "F0"), varToOpt(vm, "Mdisk0"), diff --git a/cpp/src/opacity_related.cpp b/cpp/src/opacity_related.cpp index 6caf0bc..6e14e5d 100644 --- a/cpp/src/opacity_related.cpp +++ b/cpp/src/opacity_related.cpp @@ -87,6 +87,10 @@ double OpacityRelated::Height(double R, double F) const { } +// double OpacityRelated::f_F(double xi) const { +// return a0 * xi + a1 * pow(xi, k) + a2 * pow(xi, l); +// } + double OpacityRelated::f_F(double xi) const { return a0 * xi + a1 * pow(xi, k) + a2 * pow(xi, l); } diff --git a/cpp/src/options.cpp b/cpp/src/options.cpp index b5489cc..2f6bc13 100644 --- a/cpp/src/options.cpp +++ b/cpp/src/options.cpp @@ -101,10 +101,13 @@ DiskStructureOptions::DiskStructureOptions(const po::variables_map &vm, const Ba vm["boundcond"].as(), vm["Thot"].as(), std::pow(vm["Qirr2Qvishot"].as(), 0.25), - vm["Rhot_Mdotzero_factor"].as(), + vm["Rfront_Mdotzero_factor"].as(), + vm["DIM_front_Mdot_factor"].as(), vm["check_state_approach"].as(), vm["check_Sigma_approach"].as(), vm["check_Temp_approach"].as(), + vm["DIM_front_approach"].as(), + vm["scatter_by_corona"].as(), vm["initialcond"].as(), varToOpt(vm, "F0"), varToOpt(vm, "Mdisk0"), @@ -217,14 +220,16 @@ po::options_description DiskStructureOptions::description() { "Outer-boundary movement condition\n\n" "Values:\n" " Teff: outer radius of the disk moves inwards to keep photosphere temperature of the disk larger than some value. This value is specified by --Thot option\n" - " Tirr: outer radius of the disk moves inwards to keep irradiation flux of the disk larger than some value. The value of this minimal irradiation flux is [Stefan-Boltzmann constant] * Tirr^4, where Tirr is specified by --Thot option\n" ) // fourSigmaCrit, MdotOut + " Tirr: outer radius of the disk moves inwards to keep irradiation flux of the disk larger than some value. The value of this minimal irradiation flux is [Stefan-Boltzmann constant] * Tirr^4, where Tirr is specified by --Thot option\n" + " DIM: Tirr is checked when Qirr/Qvis>Qirr2Qvishot at Rfront or at R(dotM=0),see option scatter_by_corona (yes or no, respectively); cooling front moves afterwards\n") ( "Thot", po::value()->default_value(default_Thot), "Minimum photosphere or irradiation temperature at the outer edge of the hot disk, Kelvin. For details see --boundcond description\n" ) ( "Qirr2Qvishot", po::value()->default_value(m::pow<4>(default_Tirr2Tvishot)), "Minimum Qirr / Qvis ratio at the outer edge of the hot disk to switch the control over the evolution of the hot disk radius: from temperature-based regime to Sigma-based cooling-front regime (see Lipunova et al. (2021, Section 2.4) and Eq. A.1 in Lasota et al. 2008; --alpha value is used for Sigma_plus and --alphacold value is used for Sigma_minus)\n" ) - ("Rhot_Mdotzero_factor", po::value()->default_value(default_Rhot_Mdotzero_factor), "We check conditions for cooling front at current radius mpltiplied by Rhot_Mdotzero_factor\n" ) + ("Rfront_Mdotzero_factor", po::value()->default_value(default_Rfront_Mdotzero_factor), "We check conditions for cooling front at current radius multiplied by Rfront_Mdotzero_factor\n" ) + ("DIM_front_Mdot_factor", po::value()->default_value(default_DIM_front_Mdot_factor), " = -Mdot(Rfront)/Mdot_in, see DIM_front_approach\n" ) ("check_state_approach", po::value()->default_value(default_check_state_approach), "Type of checking whether the ring is hot or cold\n\n" "Values:\n" " before2024: original version, as published in Lipunova&Malanchev (2017); Lipunova et al (2022); Avakyan et al (2024)\n" - " logic: included option for checking conditions at radius different from the radius where accretion rate is zero. See Rhot_Mdotzero_factor check_Sigma_approach, and check_Temp_approach\n") + " logic: included option for checking conditions at radius different from the radius where accretion rate is zero. See boundcond, DIM_front_approach, scatter_by_corona, check_Sigma_approach, and check_Temp_approach\n") ("check_Sigma_approach", po::value()->default_value(default_check_Sigma_approach), "Type of checking Sigma for hot or cold state\n\n" "Values:\n" " simple: assume that Sigma is proportional to R^(-3/4) between radius where Mdot = 0 and the cooling fron radius\n" @@ -232,7 +237,16 @@ po::options_description DiskStructureOptions::description() { ("check_Temp_approach", po::value()->default_value(default_check_Temp_approach), "Type of checking irradiation temperature for hot or cold state\n\n" "Values:\n" " const: assume that critical Tirr is constant, see --Thot, and --boundcond\n" - " Tavleed: assume that critical Tirr depends on Qvis/Qirr, see Tavleev et al (2023)\n" ) + " Tavleev: assume that critical Tirr depends on Qvis/Qirr, see Tavleev et al (2023)\n" + " Hameury: synthetic approach to study agreement with DIM\n") + ("DIM_front_approach", po::value()->default_value(default_DIM_front_approach), "Type of condition on accretion rate at hot disc radius\n\n" + "Values:\n" + " maxFvis: Rhot corresponds to dot M = 0 \n" + " outflow: there is an outflow from the hot zone at Rhot, see Rfront_Mdotzero_factor\n") + ("scatter_by_corona", po::value()->default_value(default_scatter_by_corona), "Presence of scattering material above the disc allowing irradiation beyond highest H/R. This flag determines 1) when DIM_front_approach=outflow starts to work; 2) Tirr=0 when scatter_by_corona=\"no\" and the ring is not seen directly from the centre.\n\n" + "Values:\n" + " yes - condition 'outflow' for DIM_front_approach takes effect always when Rhot()->default_value(default_initialcond), "Type of the initial condition for viscous torque F or surface density Sigma\n\n" "Values:\n" @@ -285,14 +299,14 @@ SelfIrradiationOptions::SelfIrradiationOptions(const po::variables_map &vm, cons } po::options_description SelfIrradiationOptions::description() { - po::options_description od("Parameters of self-irradiation:\nQirr = Cirr * (H/r / 0.05)^irrindex * L * psi / (4 pi R^2), where psi is the angular distribution of X-ray radiation\n"); + po::options_description od("Parameters of self-irradiation:\nQirr = (1-shadow) * Cirr * (H/r / 0.05)^irrindex * L * psi / (4 pi R^2), where psi is the angular distribution of X-ray radiation; shadow can be 0 (default) or 1. When scatter_by_corona=\"no\", if h/r1 < max (h/r, r()->default_value(default_Cirr), "Irradiation factor for the hot disk\n" ) ( "irrindex", po::value()->default_value(default_irrindex), "Irradiation index for the hot disk\n" ) ( "Cirrcold", po::value()->default_value(default_Cirr_cold), "Irradiation factor for the cold disk\n" ) ( "irrindexcold", po::value()->default_value(default_irrindex_cold), "Irradiation index for the cold disk\n" ) ( "h2rcold", po::value()->default_value(default_height_to_radius_cold), "Semi-height to radius ratio for the cold disk\n" ) - ( "angulardistdisk", po::value()->default_value(default_angular_dist_disk), "Angular distribution of the disk X-ray radiation. Values: isotropic, plane\n" ) + ( "angulardistdisk", po::value()->default_value(default_angular_dist_disk), "Angular distribution of the disk X-ray radiation. Values: isotropic (Psi=1), plane (Psi=2z/R)\n" ) ; return od; } @@ -362,7 +376,8 @@ CalculationOptions::CalculationOptions(const po::variables_map &vm): tauInitializer(vm), vm["Nx"].as(), vm["gridscale"].as(), - vm["starlod"].as()) { + vm["starlod"].as(), + vm["verb_level"].as() ) { if (gridscale != "log" && gridscale != "linear") { throw po::invalid_option_value("Invalid --gridscale value"); } @@ -384,6 +399,7 @@ po::options_description CalculationOptions::description() { ( "Nx", po::value()->default_value(default_Nx), "Size of calculation grid\n" ) ( "gridscale", po::value()->default_value(default_gridscale), "Type of grid for angular momentum h: log or linear\n" ) ( "starlod", po::value()->default_value(default_starlod), "Level of detail of the optical star 3-D model. The optical star is represented by a triangular tile, the number of tiles is 20 * 4^starlod\n" ) + ( "verb_level", po::value()->default_value(default_verb_level), "Verbose level of the code, 0 is muted and default\n" ) ; return od; } diff --git a/cpp/src/output.cpp b/cpp/src/output.cpp index 2cecb7f..c694913 100644 --- a/cpp/src/output.cpp +++ b/cpp/src/output.cpp @@ -62,7 +62,7 @@ BasicFreddiFileOutput::BasicFreddiFileOutput(const std::shared_ptr(&value)) { + out << "# " + << it.first.c_str() + << "=" + << *v + << "\n"; } else if (auto v = boost::any_cast >(&value)) { for (int i = 0; i < v->size(); ++i) { out << "# " @@ -118,9 +124,9 @@ BasicFreddiFileOutput::BasicFreddiFileOutput(const std::shared_ptrargs().general->dir + "/" + freddi->args().general->prefix + "_" + std::to_string(freddi->i_t()) + ".dat"); std::ofstream full_output(filename); full_output.precision(precision); - + full_output << disk_structure_header << "### t = " << sToDay(freddi->t()) << " days" + << "\n### Mdot = " << freddi->Mdot_in() << " g/s" << std::endl; const size_t last = freddi->args().flux->cold_disk ? freddi->Nx() - 1 : freddi->last(); @@ -200,7 +208,7 @@ std::vector FreddiFileOutput::initializeShortFields(const {"Mdisk", "g", "Mass of the hot disk", [freddi]() {return freddi->Mdisk();}}, {"Rhot", "Rsun", "Radius of the hot disk", [freddi]() {return cmToSun(freddi->R()[freddi->last()]);}}, {"Sigmaout", "g/cm^2", "Surface density at the outer radius of the hot disk", [freddi]() {return freddi->Sigma()[freddi->last()];}}, - {"Kirrout", "float", "Irradiation coefficient Kirr at the outer radius of the hot disk", [freddi]() {return freddi->Kirr()[freddi->last()];}}, + {"Kirrout", "float", "Irradiation coefficient Kirr at the outer radius of the hot disk w/o angular distribution: = Cirr * (z/[R*0.05])^irrindex", [freddi]() {return freddi->Kirr()[freddi->last()];}}, {"H2R", "float", "Relative semiheight at the outer radius of the hot disk", [freddi]() {return freddi->Height()[freddi->last()] / freddi->R()[freddi->last()];}}, {"Teffout", "K", "Effective tempreture at the outer radius of the hot disk", [freddi]() {return freddi->Tph()[freddi->last()];}}, {"Tirrout", "K", "Irradiation temperature (Qirr / sigma_SB)^1/4 at the outer radius of the hot disk", [freddi]() {return freddi->Tirr()[freddi->last()];}}, @@ -210,6 +218,7 @@ std::vector FreddiFileOutput::initializeShortFields(const {"Lbol", "erg/s", "Bolometric luminosity of the disk", [freddi]() {return freddi->Lbol_disk();}}, {"Fx", "erg/s/cm^2", "X-ray flux of the disk in the given energy range [emin, emax]", [freddi]() {return freddi->Lx() * freddi->angular_dist_disk(freddi->cosi()) / (FOUR_M_PI * m::pow<2>(freddi->distance()));}}, {"Fbol", "erg/s/cm^2", "Bolometric flux of the disk", [freddi]() {return freddi->Lbol_disk() * freddi->angular_dist_disk(freddi->cosi()) / (FOUR_M_PI * m::pow<2>(freddi->distance()));}}, + {"Rfront_Rhot", "float", "Ratio of cooling front radius to radius with dotM=0", [freddi]() {return freddi->Rfront_Rhot(freddi->R()[freddi->last()], freddi->Height()[freddi->last()] / freddi->R()[freddi->last()]);}}, }; const bool cold_disk = freddi->args().flux->cold_disk; const bool star = freddi->args().flux->star;