Line data Source code
1 : """
2 : Represent an ODE problem (without initial conditions) as
3 : - [`Solver`](@ref)
4 : - function `dy(t, y)` (see [Evaluation and solver state](@ref))
5 : - [`RK43.AbstractOutput`](@ref), possibly [`pipeline`](@ref)
6 :
7 : See also [`problem`](@ref)
8 : """
9 : struct Problem{N, T, S, Dy, Out}
10 : solver::Solver{N, T, S}
11 : dy::WrapDy{Dy, N, T}
12 : output::Out
13 :
14 5865 : function Problem{N, T, S, Dy, Out}(solver::Solver{N, T, S},
15 : dy::WrapDy{Dy, N, T},
16 : output::Out) where {N, T, S, Dy, Out}
17 5865 : new(solver, dy, output)
18 : end
19 : end
20 :
21 1 : function Base.show(io::IO,
22 : p::Problem{N, T, S, Dy, Out}) where {N, T, S, Dy, Out}
23 1 : print(io, "Problem{")
24 2 : show(io, p.solver)
25 1 : print(io, ", ")
26 1 : show(io, p.output)
27 1 : print(io, "}")
28 : end
29 :
30 45 : function problem(solver::Solver{N, T, S},
31 : dy::Dy,
32 : output::Out=nooutput()) where {N, T, S, Dy, Out}
33 45 : Problem{N, T, S, Dy, Out}(solver, wrap(dy, SVector{N, T}), output)
34 : end
35 :
36 0 : function problem(solver::Solver{N, T, S},
37 : dy::WrapDy{Dy, N, T},
38 : output::Out=nooutput()) where {N, T, S, Dy, Out}
39 0 : Problem{N, T, S, Dy, Out}(solver, dy, output)
40 : end
41 :
42 108 : function problem(p::Problem{N, T, S, Dy, _Out};
43 : replace::Out) where {N, T, S, Dy, _Out, Out}
44 54 : Problem{N, T, S, Dy, Out}(p.solver, p.dy, replace)
45 : end
46 :
47 4889 : function problem(p::Problem{N, T, S, Dy, Out0},
48 : out::Out) where {N, T, S, Dy, Out0, Out}
49 4889 : output = pipeline(p.output, out)
50 4889 : Problem{N, T, S, Dy, typeof(output)}(p.solver, p.dy, output)
51 : end
52 :
53 : """
54 : p = problem(solver, dy[, output])
55 : p = problem(p0::Problem, ouput) # adds output to pipeline
56 : p = problem(p0::Problem; replace=output) # replaces output
57 :
58 : See also [`RK43.Solver`](@ref), [`RK43.AbstractOutput`](@ref),
59 : [`pipeline`](@ref), [`Problem`](@ref)
60 :
61 : !!! note
62 : `replace=output` exchanges the full [`pipeline`](@ref) including
63 : any [`callback`](@ref)s or predicates!
64 :
65 : !!! warning
66 : This method does **not** copy the `solver` instance, i.e.,
67 : problem may share the same solver state. Use [`Base.similar`](@ref)
68 : in case an independent solver (state) is required (e.g., for
69 : [`solve!`](@ref) in *multiple threads*).
70 :
71 : See also [`Base.similar`](@ref)
72 : """ problem
73 :
74 : """
75 : p = similar(p0::Problem)
76 :
77 : Construct a "copy" of `p0` without associated output
78 : ([`nooutput`](@ref)).
79 :
80 : !!! note
81 : This method provides a new and uninitialized solver instance, which
82 : may be required, e.g., for [`solve!`](@ref) in *multiple threads*.
83 :
84 : !!! warning
85 : This method does **not** (and cannot) provide a "copy" of the
86 : function `dy`: The user is responsible for defining `dy`
87 : **without side effects** (or taking into account, e.g.,
88 : `Threads.threadid()` "within" `dy`).
89 :
90 : See also [`problem`](@ref)
91 : """
92 894 : Base.similar(p::Problem{N, T, S, Dy}) where {N, T, S, Dy} =
93 : Problem{N, T, S, Dy, NoOutput}(similar(p.solver), p.dy, NoOutput())
94 :
95 1 : set!(p::Problem{N, T, S, Dy}, opts::SolverOptions{T}) where {N, T, S, Dy} =
96 : set!(p.solver, opts)
97 :
98 : """
99 : opts = options(...)
100 : set!(solver, opts)
101 : set!(problem, opts)
102 :
103 : Change [`SolverOptions`](@ref)
104 :
105 : See also [`options`](@ref)
106 : """ set!
107 :
108 1732 : getsolver(p::Problem{N, T, S, Dy}) where {N, T, S, Dy} = p.solver
109 :
110 : """
111 : solver = getsolver(problem)
112 :
113 : Get solver from `problem`.
114 :
115 : See also [`problem`](@ref)
116 : """ getsolver
117 :
118 0 : getoutput(p::Problem{N, T, S, Dy}) where {N, T, S, Dy} = p.output
119 :
120 : """
121 : output = getoutput(problem)
122 :
123 : Get output from `problem`.
124 :
125 : See also [`problem`](@ref)
126 : """ getoutput
127 :
128 13213 : issuccess(p::Problem) = issuccess(p.solver)
129 :
130 : """
131 : issuccess(solver)
132 : issuccess(problem)
133 :
134 : Does `solver` or solver associated with `problem` have [`RK43.State`](@ref)
135 : `AcceptStep` (equals `RK43.Ok`), i.e., can we continue with the next
136 : step (or [`isdone`](@ref))?
137 :
138 : See also [`Solver`](@ref), [`problem`](@ref)
139 : """ issuccess
140 :
141 2 : isinitialized(p::Problem) = isinitialized(p.solver)
142 :
143 : """
144 : isinitialized(problem)
145 :
146 : Has `solver` or solver associated with `problem` been initialized
147 : by [`initialize!`](ref)?
148 :
149 : See also [`problem`](@ref), [`Solver`](@ref), [`initialize!`](ref)
150 : """
151 : isinitialized
152 :
153 12261 : initialize!(p::Problem{N, T, S, Dy, Out},
154 : t0::T, t1::T, y0::SVector{N, T}) where {N, T, S, Dy, Out} =
155 : initialize!(p.solver, p.dy, t0, t1, y0)
156 :
157 : """
158 : state = initialize!(solver, dy, t0, t1, y0)
159 : state = initialize!(problem, t0, t1, y0)
160 :
161 : Initialize `solver` or solver associated with `problem` for
162 : subsequent call(s) to [`integrate`](@ref) or [`step`](@ref) or
163 : [`solve!`](@ref), respectively.
164 :
165 : See also [`problem`](@ref), [`Solver`](@ref), [`isinitialized`](@ref),
166 : [`RK43.State`](@Ref)
167 : """
168 : initialize!
169 :
170 3 : state(p::Problem{N, T, S, Dy, Out}) where {N, T, S, Dy, Out} = state(p.solver)
171 :
172 : """
173 : s = state(solver)
174 : s = state(problem)
175 :
176 : Query [`RK43.State`](@ref) of `solver` or solver associated with `problem`.
177 :
178 : See also [`problem`](@ref), [`Solver`](@ref), [`RK43.State`](@ref),
179 : [`issuccess`](@ref)
180 : """ state
181 :
182 204 : isforward(p::Problem) = isforward(p.solver)
183 :
184 : """
185 : isforward(solver)
186 : isforward(problem)
187 :
188 : Is integration in *forward direction* (i.e., increasing time) for `solver` or
189 : solver associated with `problem`?
190 :
191 : !!! warning "precondition"
192 : The solver must be initialized (i.e., have `t0` and `t1` defined),
193 : or the result is undefined.
194 :
195 : See also [`problem`](@ref), [`Solver`](@ref)
196 : """ isforward
197 :
198 2051 : isdone(p::Problem) = isdone(p.solver)
199 :
200 : """
201 : isdone(solver)
202 : isdone(problem)
203 :
204 : Is integration finished (i.e., reached destination) for `solver` or
205 : solver associated with `problem`?
206 :
207 : !!! warning "precondition"
208 : The solver must be initialized (i.e., have `t0` and `t1` defined),
209 : or the result is undefined.
210 :
211 : See also [`problem`](@ref), [`Solver`](@ref)
212 : """ isdone
213 :
214 7755 : pos(p::Problem) = pos(p.solver)
215 :
216 : """
217 : pos(solver)
218 : pos(problem)
219 :
220 : Get position tuple `t, y, dy` ,for `solver` or solver associated with
221 : `problem`.
222 :
223 : See also [`problem`](@ref), [`Solver`](@ref)
224 : """ pos
225 :
226 1064 : emptysolution(::Solver{N, T, S}) where {N, T, S} = hermitespline(SVector{N, T})
227 :
228 1050 : emptysolution(p::Problem{N, T, S, Dy, Out}) where {N, T, S, Dy, Out} =
229 : emptysolution(p.solver)
230 :
231 : """
232 : sol = emptysolution(solver)
233 : sol = emptysolution(problem)
234 :
235 : Construct empty solution for `solver` or solvler associated with
236 : `problem`.
237 :
238 : Solutions are represented as splines, i.e., piecewise polynomial curves,
239 : like [`HermiteSpline`](https://visual2.cs.ovgu.de/roessl/HermiteSpline.jl)
240 :
241 : !!! note
242 : `emptysolution` determines the (data) *type* of the solution from
243 : a solver (or problem).
244 :
245 : See also [`problem`](@ref), [`Solver`](@ref)
246 : """ emptysolution
247 :
248 12259 : function solve!(p::Problem{N, T, S, Dy, Out},
249 : t0::T, t1::T, y0::SVector{N, T}) where {N, T, S, Dy, Out}
250 12259 : initialize!(p, t0, t1, y0)
251 12259 : solve!(p)
252 : end
253 :
254 4886 : function solve!(p0::Problem{N, T, S, Dy, Out0},
255 : t0::T, t1::T, y0::SVector{N, T},
256 : out::Out) where {N, T, S, Dy, Out0, Out}
257 4886 : solve!(problem(p0, out), t0, t1, y0)
258 : end
259 :
260 : const MAXSTEPS = 10000 # fallback if opts.maxsteps == 0
261 :
262 12260 : function solve!(p::Problem{N, T, S, Dy, Out}) where {N, T, S, Dy, Out}
263 12260 : @massert isinitialized(p) "missing initialization"
264 :
265 12260 : issuccess(p) || return state(p)
266 :
267 12260 : solver = p.solver
268 :
269 12260 : maxsteps = max(MAXSTEPS, solver.opts.maxsteps)
270 17136 : ostate = _check_output(p.output(solver))
271 :
272 0 : @massert issuccess(ostate) "output() must succeed for initial point ($(solver.t),$(solver.y))"
273 :
274 881961 : while !isdone(solver) && (solver.nsteps < maxsteps)
275 869953 : newstate = step!(solver, p.dy)
276 :
277 870192 : if issuccess(newstate) || (newstate == OutOfDomain)
278 1061549 : ostate = _check_output(p.output(solver))
279 :
280 869953 : if !issuccess(ostate)
281 325 : if ostate == Stopped
282 4 : newstate = Stopped
283 : # TODO [alg:upstream] really commit after Stopped
284 : # TODO [alg:upstream] enable callback to modify step
285 321 : elseif ostate == OutOfDomain
286 : # NOTE: step! handles bisection as well, p.output may "recheck"
287 320 : solver.absh = bisectedstepsize(solver)
288 320 : if solver.absh > hmin(solver)
289 312 : solver.nsteps -= 1
290 312 : solver.rejected = true
291 312 : continue # repeat step
292 : end
293 :
294 8 : newstate = ostate
295 : else
296 1 : @error "invalid state $(ostate): output() may indicate Stop or OutOfDomain!"
297 1 : return (solver.state = Failed)
298 : end
299 : end
300 :
301 869640 : commit!(solver)
302 : end
303 :
304 869891 : issuccess(newstate) || return (solver.state = newstate)
305 869701 : end # while !done
306 :
307 : # Note: special case for output if done !?!
308 :
309 12008 : if solver.nsteps > maxsteps
310 0 : @warn "exceeded maximum number of $(maxsteps) steps"
311 0 : return (solver.state = Stopped)
312 : end
313 :
314 12008 : solver.state
315 : end
316 :
317 : """
318 : state = solve!(problem) # requires initialize!(problem,...)
319 : state = solve!(problem, t0, t1, y0)
320 : state = solve!(problem, t0, t1, y0, output) # solve with additional output
321 :
322 : Solve `problem` and return [`RK43.State`](@ref). The solution is constructed
323 : (as a "side-effect") from the [`RK43.AbstractOutput`](@ref) specified by `problem`.
324 :
325 : !!! note
326 : If the solver stops before reaching `t1` (e.g., due to an `OutOfDomain`),
327 : There will be *no output* after the `tend = first(pos(solver))`.
328 :
329 : See also [`problem`](@ref), [`RK43.state`](@ref), [`initialize!`](@ref)
330 : """ solve!
|
