diff --git a/src/disc_solver/file_format/_containers.py b/src/disc_solver/file_format/_containers.py index f69d2408..0cd2f681 100644 --- a/src/disc_solver/file_format/_containers.py +++ b/src/disc_solver/file_format/_containers.py @@ -63,6 +63,7 @@ class ConfigInput: jump_before_sonic = attr.ib(default=None) use_taylor_jump = attr.ib(default="True") mcmc_vars = attr.ib(default=None) + v_θ_sonic_crit = attr.ib(default=None) @attr.s @@ -93,6 +94,7 @@ class SolutionInput: jump_before_sonic = attr.ib(default=None) use_taylor_jump = attr.ib(default=True) mcmc_vars = attr.ib(default=None) + v_θ_sonic_crit = attr.ib(default=None) class Problems(MutableMapping): diff --git a/src/disc_solver/file_format/_dumpers.py b/src/disc_solver/file_format/_dumpers.py index 13f43ba3..77e82798 100644 --- a/src/disc_solver/file_format/_dumpers.py +++ b/src/disc_solver/file_format/_dumpers.py @@ -69,7 +69,7 @@ def _solution_dumper(solution): ) -@ds_registry.dumper(ConfigInput, "ConfigInput", version=9) +@ds_registry.dumper(ConfigInput, "ConfigInput", version=10) def _config_dumper(config_input): return GroupContainer( attrs={ @@ -98,10 +98,11 @@ def _config_dumper(config_input): }, jump_before_sonic=config_input.jump_before_sonic, mcmc_vars=config_input.mcmc_vars, + v_θ_sonic_crit=config_input.v_θ_sonic_crit, ) -@ds_registry.dumper(SolutionInput, "SolutionInput", version=9) +@ds_registry.dumper(SolutionInput, "SolutionInput", version=10) def _input_dumper(solution_input): return GroupContainer( attrs={ @@ -129,6 +130,7 @@ def _input_dumper(solution_input): }, jump_before_sonic=solution_input.jump_before_sonic, mcmc_vars=solution_input.mcmc_vars, + v_θ_sonic_crit=solution_input.v_θ_sonic_crit, ) diff --git a/src/disc_solver/file_format/_loaders.py b/src/disc_solver/file_format/_loaders.py index bf08c95f..3d6b21ba 100644 --- a/src/disc_solver/file_format/_loaders.py +++ b/src/disc_solver/file_format/_loaders.py @@ -432,6 +432,37 @@ def _config_loader9(group): ) +@ds_registry.loader("ConfigInput", version=10) +def _config_loader10(group): + return ConfigInput( + start=group.attrs["start"], + stop=group.attrs["stop"], + taylor_stop_angle=group.attrs["taylor_stop_angle"], + max_steps=group.attrs["max_steps"], + num_angles=group.attrs["num_angles"], + label=group.attrs["label"], + relative_tolerance=group.attrs["relative_tolerance"], + absolute_tolerance=group.attrs["absolute_tolerance"], + γ=group.attrs["γ"], + nwalkers=group.attrs["nwalkers"], + iterations=group.attrs["iterations"], + threads=group.attrs["threads"], + target_velocity=group.attrs["target_velocity"], + split_method=group.attrs["split_method"], + v_rin_on_c_s=group.attrs["v_rin_on_c_s"], + v_a_on_c_s=group.attrs["v_a_on_c_s"], + c_s_on_v_k=group.attrs["c_s_on_v_k"], + η_O=group.attrs["η_O"], + η_H=group.attrs["η_H"], + η_A=group.attrs["η_A"], + jump_before_sonic=group["jump_before_sonic"], + η_derivs=group.attrs["η_derivs"], + use_taylor_jump=group.attrs["use_taylor_jump"], + mcmc_vars=group["mcmc_vars"], + v_θ_sonic_crit=group["v_θ_sonic_crit"], + ) + + @ds_registry.loader("SolutionInput", version=1) def _input_loader(group): return SolutionInput( @@ -694,6 +725,44 @@ def _input_loader9(group): ) +@ds_registry.loader("SolutionInput", version=10) +def _input_loader10(group): + if group["jump_before_sonic"] is None: + jump_before_sonic = None + else: + jump_before_sonic = group["jump_before_sonic"] + if group["v_θ_sonic_crit"] is None: + v_θ_sonic_crit = None + else: + v_θ_sonic_crit = group["v_θ_sonic_crit"] + return SolutionInput( + start=group.attrs["start"], + stop=group.attrs["stop"], + taylor_stop_angle=group.attrs["taylor_stop_angle"], + max_steps=group.attrs["max_steps"], + num_angles=group.attrs["num_angles"], + relative_tolerance=group.attrs["relative_tolerance"], + absolute_tolerance=group.attrs["absolute_tolerance"], + γ=group.attrs["γ"], + nwalkers=group.attrs["nwalkers"], + iterations=group.attrs["iterations"], + threads=group.attrs["threads"], + target_velocity=group.attrs["target_velocity"], + split_method=group.attrs["split_method"], + v_rin_on_c_s=group.attrs["v_rin_on_c_s"], + v_a_on_c_s=group.attrs["v_a_on_c_s"], + c_s_on_v_k=group.attrs["c_s_on_v_k"], + η_O=group.attrs["η_O"], + η_H=group.attrs["η_H"], + η_A=group.attrs["η_A"], + jump_before_sonic=jump_before_sonic, + η_derivs=group.attrs["η_derivs"], + use_taylor_jump=group.attrs["use_taylor_jump"], + mcmc_vars=group["mcmc_vars"], + v_θ_sonic_crit=v_θ_sonic_crit, + ) + + @ds_registry.loader("Run", version=1) def _run_loader(group): return Run( diff --git a/src/disc_solver/solve/solution.py b/src/disc_solver/solve/solution.py index f705aba3..0ffdb4e0 100644 --- a/src/disc_solver/solve/solution.py +++ b/src/disc_solver/solve/solution.py @@ -35,7 +35,6 @@ INTEGRATOR = "cvode" LINSOLVER = "dense" COORDS = "spherical midplane 0" -V_θ_SONIC_CRIT = 0.05 log = logbook.Logger(__name__) @@ -49,7 +48,8 @@ def ode_system( *, γ, a_0, norm_kepler_sq, init_con, θ_scale=float_type(1), - with_taylor=False, η_derivs=True, store_internal=True, use_E_r=False + with_taylor=False, η_derivs=True, store_internal=True, use_E_r=False, + v_θ_sonic_crit=None ): """ Set up the system we are solving for. @@ -155,7 +155,7 @@ def rhs_equation(x, params, derivs): deriv_v_φ_taylor if with_taylor else deriv_v_φ_normal ) - if np_abs(1 - v_θ) < V_θ_SONIC_CRIT: + if v_θ_sonic_crit is not None and np_abs(1 - v_θ) < v_θ_sonic_crit: deriv_v_θ = deriv_v_θ_sonic( a_0=a_0, ρ=ρ, B_r=B_r, B_φ=B_φ, B_θ=B_θ, η_O=η_O, η_H=η_H, η_A=η_A, θ=θ, v_r=v_r, v_φ=v_φ, deriv_v_r=deriv_v_r, @@ -415,7 +415,7 @@ def main_solution( absolute_tolerance=float_type(1e-10), max_steps=500, onroot_func=None, jump_before_sonic=None, tstop=None, ontstop_func=None, η_derivs=True, store_internal=True, root_func=None, root_func_args=None, - θ_scale=float_type(1), use_E_r=False + θ_scale=float_type(1), use_E_r=False, v_θ_sonic_crit=None ): """ Find solution @@ -439,7 +439,7 @@ def main_solution( γ=γ, a_0=a_0, norm_kepler_sq=norm_kepler_sq, init_con=system_initial_conditions, η_derivs=η_derivs, store_internal=store_internal, with_taylor=False, θ_scale=θ_scale, - use_E_r=use_E_r, + use_E_r=use_E_r, v_θ_sonic_crit=v_θ_sonic_crit, ) solver = ode( @@ -505,6 +505,7 @@ def solution( η_derivs = soln_input.η_derivs use_taylor_jump = soln_input.use_taylor_jump jump_before_sonic = soln_input.jump_before_sonic + v_θ_sonic_crit = soln_input.v_θ_sonic_crit if with_taylor and use_E_r: raise SolverError( @@ -551,6 +552,7 @@ def solution( η_derivs=η_derivs, store_internal=store_internal, jump_before_sonic=jump_before_sonic, root_func=root_func, root_func_args=root_func_args, θ_scale=θ_scale, use_E_r=use_E_r, + v_θ_sonic_crit=v_θ_sonic_crit, ) if jump_before_sonic is not None: @@ -574,7 +576,7 @@ def solution( onroot_func=onroot_func, tstop=tstop, ontstop_func=ontstop_func, η_derivs=η_derivs, store_internal=store_internal, root_func=root_func, root_func_args=root_func_args, - θ_scale=θ_scale, use_E_r=use_E_r, + θ_scale=θ_scale, use_E_r=use_E_r, v_θ_sonic_crit=v_θ_sonic_crit, ) if store_internal and taylor_stop_angle is not None: