Unphysical structure when radius grid is smaller than sonic radius

[ladsantos]

This issue was pointed out by Yassin Jaziri (U. of Geneva). When a p-winds model is run in a radius grid smaller than the sonic radius, the code produces a discontinuity in the densities and velocities profile. For example, if the quickstart notebook is run with r = np.logspace(0, np.log10(3), 100), the resulting is structure is as follows:

Furthermore, for even smaller grids, I get an error when calculating the H ion_fraction:

RuntimeError                              Traceback (most recent call last)
/var/folders/0q/ph_vqx8j58s5t8b4smvzpcf80002vw/T/ipykernel_32187/1288101433.py in <module>
      2 r = np.logspace(0, np.log10(3), 100)  # Radial distance profile in unit of planetary radii
      3 
----> 4 f_r, mu_bar = hydrogen.ion_fraction(r, R_pl, T_0, h_fraction, 
      5                             m_dot, M_pl, mu_0,
      6                             spectrum_at_planet=spectrum, exact_phi=True,

~/Coding/p-winds/p_winds/hydrogen.py in ion_fraction(radius_profile, planet_radius, temperature, h_fraction, mass_loss_rate, planet_mass, mean_molecular_weight_0, spectrum_at_planet, flux_euv, initial_f_ion, relax_solution, convergence, max_n_relax, exact_phi, return_mu, **options_solve_ivp)
    461 
    462             # And solve it again
--> 463             sol = solve_ivp(_fun, (r[0], r[-1],), np.array([initial_f_ion, ]),
    464                             t_eval=r, args=(phi, k2), **options_solve_ivp)
    465             f_r = sol['y'][0]

~/Applications/miniconda3/envs/dev/lib/python3.9/site-packages/scipy/integrate/_ivp/ivp.py in solve_ivp(fun, t_span, y0, method, t_eval, dense_output, events, vectorized, args, **options)
    574     status = None
    575     while status is None:
--> 576         message = solver.step()
    577 
    578         if solver.status == 'finished':

~/Applications/miniconda3/envs/dev/lib/python3.9/site-packages/scipy/integrate/_ivp/base.py in step(self)
    179         else:
    180             t = self.t
--> 181             success, message = self._step_impl()
    182 
    183             if not success:

~/Applications/miniconda3/envs/dev/lib/python3.9/site-packages/scipy/integrate/_ivp/rk.py in _step_impl(self)
    142             h_abs = np.abs(h)
    143 
--> 144             y_new, f_new = rk_step(self.fun, t, y, self.f, h, self.A,
    145                                    self.B, self.C, self.K)
    146             scale = atol + np.maximum(np.abs(y), np.abs(y_new)) * rtol

~/Applications/miniconda3/envs/dev/lib/python3.9/site-packages/scipy/integrate/_ivp/rk.py in rk_step(fun, t, y, f, h, A, B, C, K)
     62     for s, (a, c) in enumerate(zip(A[1:], C[1:]), start=1):
     63         dy = np.dot(K[:s].T, a[:s]) * h
---> 64         K[s] = fun(t + c * h, y + dy)
     65 
     66     y_new = y + h * np.dot(K[:-1].T, B)

~/Applications/miniconda3/envs/dev/lib/python3.9/site-packages/scipy/integrate/_ivp/base.py in fun(t, y)
    136         def fun(t, y):
    137             self.nfev += 1
--> 138             return self.fun_single(t, y)
    139 
    140         self.fun = fun

~/Applications/miniconda3/envs/dev/lib/python3.9/site-packages/scipy/integrate/_ivp/base.py in fun_wrapped(t, y)
     18 
     19     def fun_wrapped(t, y):
---> 20         return np.asarray(fun(t, y), dtype=dtype)
     21 
     22     return fun_wrapped, y0

~/Applications/miniconda3/envs/dev/lib/python3.9/site-packages/scipy/integrate/_ivp/ivp.py in <lambda>(t, x, fun)
    512         # additional parameters.  Pass in the original fun as a keyword
    513         # argument to keep it in the scope of the lambda.
--> 514         fun = lambda t, x, fun=fun: fun(t, x, *args)
    515         jac = options.get('jac')
    516         if callable(jac):

~/Coding/p-winds/p_winds/hydrogen.py in _fun(_r, _f, _phi, _k2)
    400             _phi_prime = np.exp(-_t)*_phi
    401         _v_guess = _v_fun(_r)
--> 402         _v, _rho = parker.structure(_r, _v_guess)
    403         # In terms 1 and 2 we use the values of k2 and phi from above
    404         term1 = (1. - _f) / _v * _phi_prime

~/Coding/p-winds/p_winds/parker.py in structure(r, v_guess)
    194 
    195     if v_guess is not None:
--> 196         velocity_r = so.newton(_eq_to_solve, x0=v_guess, args=(r,))
    197 
    198     # Otherwise, we set it automatically: If the radius is below the sonic

~/Applications/miniconda3/envs/dev/lib/python3.9/site-packages/scipy/optimize/zeros.py in newton(func, x0, fprime, args, tol, maxiter, fprime2, x1, rtol, full_output, disp)
    338                             " Failed to converge after %d iterations, value is %s."
    339                             % (itr + 1, p1))
--> 340                         raise RuntimeError(msg)
    341                     warnings.warn(msg, RuntimeWarning)
    342                 p = (p1 + p0) / 2.0

RuntimeError: Tolerance of -87175409.78639884 reached. Failed to converge after 5 iterations, value is -87193373.31058452.