Allow controlling of the change to when l’Hôpital is used

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):

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

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(

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: