Line data    Source code

       1             : # TODO: abstract stepper type Stepper<:AbstractStepper
       2             : 
       3             : # TODO: docs
       4             : 
       5             : """
       6             : Solver for vector type `SVector{N, T}` using `Stepper`
       7             : (e.g., [`RK43StepperState`](@ref)).
       8             : 
       9             : See also [`rk43solver`](@ref)
      10             : """
      11             : mutable struct Solver{N, T, Stepper}
      12             :     stepper::Stepper
      13             : 
      14             :     y::SVector{N, T}
      15             :     dy::SVector{N, T}
      16             :     ynew::SVector{N, T}
      17             :     dynew::SVector{N, T}
      18             :     t::T
      19             :     tnew::T
      20             : 
      21             :     tdir::T
      22             :     tfinal::T
      23             :     absh::T
      24             :     relerr::T
      25             : 
      26             :     nsteps::Int
      27             :     state::State
      28             :     rejected::Bool
      29             :     opts::SolverOptions{T}
      30             : 
      31         932 :     function Solver{N, T, Stepper}(opts::SolverOptions{T}=options()) where {N, T, Stepper}
      32         932 :         unknown = undefined(SVector{N, T})
      33         932 :         new(Stepper(), # may be undefined state
      34        3728 :             ntuple(_ -> unknown, 4)...,
      35        5592 :             ntuple(_ -> T(NaN), 6)...,
      36             :             0, Uninitialized, false, opts)
      37             :     end
      38             : end
      39             : 
      40             : """
      41             :     s = copy(s0::Solver{N, T, Stepper})
      42             : 
      43             : Provide a (deep) copy of solver `s0`.
      44             : """
      45           9 : function Base.copy(s0::Solver{N, T, Stepper}) where {N, T, Stepper}
      46           9 :     s = similar(s0)
      47             : 
      48           9 :     s.stepper = copy(s0.stepper)
      49             : 
      50           9 :     s.y = s0.y
      51           9 :     s.dy = s0.dy
      52           9 :     s.ynew = s0.ynew
      53           9 :     s.dynew = s0.dynew
      54           9 :     s.t = s0.t
      55           9 :     s.tnew = s0.tnew
      56             : 
      57           9 :     s.tdir = s0.tdir
      58           9 :     s.tfinal = s0.tfinal
      59           9 :     s.absh = s0.absh
      60           9 :     s.relerr = s0.relerr
      61             : 
      62           9 :     s.nsteps = s0.nsteps
      63           9 :     s.state = s0.state
      64           9 :     s.rejected = s0.rejected
      65           9 :     s.opts = s0.opts
      66             : 
      67           9 :     s
      68             : end
      69             : 
      70             : """
      71             :     s = similar(s0::Solver)
      72             : 
      73             : Provide a new and uninitialzed "copy" of `s`.
      74             : """
      75         903 : Base.similar(solver::Solver{N, T, S}) where {N, T, S} =
      76             :     Solver{N, T, S}(solver.opts)
      77             : 
      78           2 : function Base.show(io::IO, solver::Solver)
      79           2 :     print(io, "Solver{")
      80           2 :     show(io, solver.stepper)
      81           2 :     print(io, "}[", pos(solver), ", ")
      82           4 :     show(io, state(solver))
      83           2 :     print(io, "]")
      84             : end
      85             : 
      86           1 : function set!(solver::Solver{N, T, S}, opts::SolverOptions{T}) where {N, T, S}
      87           1 :     solver.opts = opts
      88           1 :     solver
      89             : end
      90             : 
      91           0 : set!(solver::Solver{N, T, S}, opts::SolverOptions) where {N, T, S} =
      92             :     set!(solver, options(T, opts))
      93             : 
      94             : 
      95           0 : wrap(dy::F, ::Solver{N, T}) where {F, N, T} = wrap(dy, SVector{N, T})
      96           0 : guard(dy, solver::Solver{N, T}) where {N, T} = wrap(dy, SVector{N, T}) |> guard
      97     1308397 : guard(dy::WrapDy{F, N, T}, solver::Solver{N, T}) where {F, N, T} = guard(dy)
      98             : 
      99             : 
     100          30 : rk43solver(::Type{SVector{N, T}}, opts=options(T)) where {N, T} =
     101             :     Solver{N, T, RK43StepperState{N, T}}(options(T, opts))
     102           0 : rk43solver(::SVector{N, T}, opts=options(T)) where {N, T} =
     103             :     rk43solver(SVector{N, T}, options(T, opts))
     104             : 
     105       13214 : issuccess(solver::Solver) = issuccess(solver.state)
     106           4 : isinitialized(solver::Solver) = (solver.state != Uninitialized)
     107             : 
     108           8 : state(solver::Solver) = solver.state
     109             : 
     110        7767 : pos(solver::Solver{N, T}) where {N, T} = (solver.t, solver.y, solver.dy)
     111             : 
     112           0 : tentativepos(solver::Solver{N, T}) where {N, T} =
     113             :     (solver.tnew, solver.ynew, solver.dynew)
     114             : 
     115     2567989 : stepsize(solver::Solver) = solver.tdir * solver.absh
     116      147858 : isforward(solver::Solver) = solver.tdir > 0
     117             : 
     118      885746 : isdone(solver::Solver) = (solver.t == solver.tfinal)
     119             : 
     120      870273 : hmin(solver::Solver{N, T}) where {N, T} = 16eps(max(abs(solver.t), T(1e-10)))
     121             : 
     122             : """
     123             :     isinside(solver, t)
     124             : 
     125             : Test `t ∈ [solver.t, solver.tnew]`.
     126             : """
     127      307519 : function isinside(solver::Solver{N, T}, t::T) where {N, T}
     128      307519 :     if solver.tdir > 0
     129      282995 :         solver.t <= t <= solver.tnew
     130             :     else
     131       24524 :         solver.tnew <= t <= solver.t
     132             :     end
     133             : end
     134             : 
     135             : """
     136             : Commit [`step!`](@ref).
     137             : """
     138      869640 : function commit!(solver::Solver{N, T}) where {N, T}
     139      869640 :     solver.t = solver.tnew
     140      869640 :     solver.y = solver.ynew
     141      869640 :     solver.dy = solver.dynew
     142             : end
     143             : 
     144             : # TODO: ctor symbols ... not rk43solver (odesolver(...))
     145             : 
     146             : # TODO: Solver/ Stepper traits -- order{0,1}, solution, ...