fronts.solve

fronts.solve(D, i, b, radial=False, ob=0.0, itol=0.001, d_dob_hint=None, d_dob_bracket=None, method='implicit', maxiter=100, verbose=0)

Solve a problem with a Dirichlet boundary condition.

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

\[\begin{split}\begin{cases} \dfrac{\partial\theta}{\partial t} = \nabla\cdot\left[D(\theta)\dfrac{\partial\theta}{\partial r} \mathbf{\hat{r}}\right ] & r>r_b(t),t>0\\ \theta(r, 0) = \theta_i & r>0 \\ \theta(r_b(t), t) = \theta_b & t>0 \\ r_b(t) = o_b\sqrt t \end{cases}\end{split}\]
Parameters:
  • D (callable or sympy.Expression or str or float) –

    Callable that evaluates \(D\) and its derivatives, obtained from the fronts.D module or defined in the same manner—i.e.:

    • D(theta) evaluates and returns \(D\) at theta

    • D(theta, 1) returns both the value of \(D\) and its first derivative at theta

    • D(theta, 2) returns the value of \(D\), its first derivative, and its second derivative at theta

    When called by this function, theta is always a single float. However, calls as D(theta) should also accept a NumPy array argument.

    Alternatively, instead of a callable, the argument can be the expression of \(D\) in the form of a string or sympy.Expression with a single variable. In this case, the solver will differentiate and evaluate the expression as necessary.

  • i (float) – Initial condition, \(\theta_i\).

  • b (float) – Imposed boundary value, \(\theta_b\).

  • radial ({False, 'cylindrical', 'polar', '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' or 'polar': \(\mathbf{\hat{r}}\) is the radial unit vector in a cylindrical or polar coordinate system

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

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

  • itol (float, optional) – Absolute tolerance for the initial condition.

  • d_dob_hint (None or float, optional) – Optional hint to the solver. If given, it should be a number close to the expected value of the derivative of \(\theta\) with respect to the Boltzmann variable at the boundary (i.e., \(d\theta/do|_b\)) in the solution to be found. This parameter is typically not needed.

  • d_dob_bracket (None or sequence of two floats, optional) – Optional search interval that brackets the value of \(d\theta/do|_b\) in the solution. If given, the solver will use bisection to find a solution in which \(d\theta/do|_b\) falls inside that interval (a ValueError will be raised for an incorrect interval). This parameter cannot be passed together with a d_dob_hint. It is also not needed in typical usage.

  • method ({'implicit', 'explicit'}, optional) –

    Selects the integration method used by the solver:

    • 'implicit' (default): uses a Radau IIA implicit method of order 5. A sensible default choice that will work for any problem

    • 'explicit': uses the DOP853 explicit method of order 8. As an explicit method, it trades off general solver robustness and accuracy for faster results in “well-behaved” cases. With this method, the second derivative of \(D\) is not needed

  • 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 nonnegative.

  • 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: also 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:

  • o (numpy.ndarray, shape (n,)) – Final solver mesh, in terms of the Boltzmann variable.

  • niter (int) – Number of iterations required to find the solution.

  • d_dob_bracket (sequence of two floats or None) – If available, an interval that contains the value of \(d\theta/do|_b\). May be used as the input d_dob_bracket in a subsequent call with a smaller itol for the same problem in order to avoid reduntant iterations. Whether this interval is available or not depends on the strategy used internally by the solver; in particular, this field is never None if a d_dob_bracket is passed when calling the function.

Return type:

Solution

Notes

This function works by transforming the partial differential equation with the Boltzmann transformation using ode() and then solving the resulting ODE repeatedly with the chosen integration method as implemented in the scipy.integrate module and a custom shooting algorithm. The boundary condition is satisfied exactly as the starting point, and the algorithm iterates with different values of \(d\theta/do|_b\) until it finds the solution that also verifies the initial condition within the specified tolerance. Trial values of \(d\theta/do|_b\) are selected automatically by default (using heuristics, which can also take into account an optional hint if passed by the user), or by bisecting an optional search interval. This scheme assumes that \(d\theta/do|_b\) varies continuously with \(\theta_i\).

References

[1] GERLERO, G. S.; BERLI, C. L. A.; KLER, P. A. Open-source high-performance software packages for direct and inverse solving of horizontal capillary flow. Capillarity, 2023, vol. 6, no. 2, pp. 31-40.