fronts.solve

fronts.solve(D, Si, Sb, dS_dob_bracket=(-1.0, 1.0), radial=False, ob=0.0, Si_tol=1e-06, maxiter=100, verbose=0)

Solve an instance of the general problem.

Given a positive function D, scalars \(S_i\), \(S_b\) and \(o_b\), and coordinate unit vector \(\mathbf{\hat{r}}\), finds a function S of r and t such that:

\[\begin{split}\begin{cases} \dfrac{\partial S}{\partial t} = \nabla\cdot\left[D\left(S\right)\dfrac{\partial S}{\partial r} \mathbf{\hat{r}}\right ] & r>r_b(t),t>0\\ S(r, 0) = S_i & r>0 \\ S(r_b(t), t) = S_b & t>0 \\ r_b(t) = o_b\sqrt t \end{cases}\end{split}\]
Parameters

  • D (callable) – Twice-differentiable function that maps the range of S to positive values. It can be called as D(S) to evaluate it at S. It can also be called as D(S, n) with n equal to 1 or 2, in which case the first n derivatives of the function evaluated at the same S are included (in order) as additional return values. While mathematically a scalar function, D operates in a vectorized fashion with the same semantics when S is a numpy.ndarray.

  • Si (float) – \(S_i\), the initial value of S in the domain.

  • Sb (float) – \(S_b\), the value of S imposed at the boundary.

  • dS_dob_bracket ((float, float), optional) – Search interval that contains the value of the derivative of S with respect to the Boltzmann variable o (i.e., \(dS/do\)) at the boundary in the solution. The interval can be made as wide as desired, at the cost of additional iterations required to obtain the solution. To refine a solution obtained previously with this same function, pass in that solution’s final dS_dob_bracket. This parameter is always checked and a ValueError is raised if a dS_dob_bracket is found not to be valid for the problem.

  • radial ({False, 'cylindrical', 'spherical'}, optional) –

    Choice of coordinate unit vector \(\mathbf{\hat{r}}\). Must be one of the following:

    • False (default)

      \(\mathbf{\hat{r}}\) is any coordinate unit vector in rectangular (Cartesian) coordinates, or an axial unit vector in a cylindrical coordinate system

    • 'cylindrical'

      \(\mathbf{\hat{r}}\) is the radial unit vector in a cylindrical coordinate system

    • 'spherical'

      \(\mathbf{\hat{r}}\) is the radial unit vector in a spherical coordinate system

  • ob (float, optional) – \(o_b\), which determines the behavior of the boundary. The default is zero, which implies that the boundary always exists at \(r=0\). It must be strictly positive if radial is not False. Be aware that a non-zero value implies a moving boundary.

  • Si_tol (float, optional) – Absolute tolerance for \(S_i\).

  • maxiter (int, optional) – Maximum number of iterations. A RuntimeError will be raised if the specified tolerance is not achieved within this number of iterations. Must be >= 2.

  • verbose ({0, 1, 2}, optional) –

    Level of algorithm’s verbosity. Must be one of the following:

    • 0 (default) : work silently.

    • 1 : display a termination report.

    • 2 : display progress during iterations.

Returns

solution – See Solution for a description of the solution object. Additional fields specific to this solver are included in the object:

  • onumpy.ndarray, shape (n,)

    Final solver mesh, in terms of the Boltzmann variable o.

  • niterint

    Number of iterations required to find the solution.

  • dS_dob_bracket(float, float)

    Subinterval of dS_dob_bracket that contains the value of \(dS/do\) at the boundary in the solution. May be used in a subsequent call with a smaller Si_tol to avoid reduntant iterations if wanting to refine a previously obtained solution.

Return type

Solution

See also

solve_from_guess()

Notes

Given the expression of \(r_b\) which specifies the location of the boundary, a fixed boundary can be had only if \(o_b=0\). Any other \(o_b\) implies a moving boundary. This restriction affects radial problems in particular.

This function works by transforming the partial differential equation with the Boltzmann transformation using ode and then solving the resulting ODE repeatedly using the ‘Radau’ method as implemented in scipy.integrate.solve_ivp and a custom shooting algorithm. The boundary condition is satisfied exactly as the starting point, and the shooting algorithm iterates with different values of \(dS/do\) at the boundary (chosen from within dS_dob_bracket using bisection) until it finds the solution that also satisfies the initial condition within the specified tolerance. This scheme assumes that \(dS/do\) at the boundary varies continuously with \(S_i\).