MRKISS 2025-09-10
A tiny library with zero dependencies that aims to make it easy to use & experiment with explicit Runge-Kutta methods.
Loading...
Searching...
No Matches
mrkiss_solvers_wt.f90
Go to the documentation of this file.
1! -*- Mode:F90; Coding:us-ascii-unix; fill-column:129 -*-
2!.!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!.H.S.!!
3!>
4!! @file mrkiss_solvers_wt.f90
5!! @author Mitch Richling http://www.mitchr.me/
6!! @brief RK Solvers.@EOL
7!! @keywords runge kutta embedded butcher tableau
8!! @std F2023
9!! @see https://github.com/richmit/MRKISS
10!! @copyright
11!! @parblock
12!! Copyright (c) 2025, Mitchell Jay Richling <http://www.mitchr.me/> All rights reserved.
13!!
14!! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following
15!! conditions are met:
16!!
17!! 1. Redistributions of source code must retain the above copyright notice, this list of conditions, and the following
18!! disclaimer.
19!!
20!! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions, and the following
21!! disclaimer in the documentation and/or other materials provided with the distribution.
22!!
23!! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products
24!! derived from this software without specific prior written permission.
25!!
26!! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
27!! INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
28!! DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
29!! EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
30!! USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31!! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
32!! OF THE POSSIBILITY OF SUCH DAMAGE.
33!! @endparblock
34!.!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!.H.E.!!
35
36
37!----------------------------------------------------------------------------------------------------------------------------------
38!> Solvers.
39!!
40module mrkiss_solvers_wt
41 implicit none
42 private
43
44 public :: step_rk4, step_rkf45, step_dp54 ! Test one step solvers
45 public :: step_all, step_one, step_richardson ! One step solvers
46 public :: fixed_t_steps, fixed_y_steps, sloppy_fixed_y_steps, adaptive_steps ! Many step solvers
47 public :: fixed_t_steps_between, interpolate_solution ! Meta-many step solvers
48
49 !--------------------------------------------------------------------------------------------------------------------------------
50 !> Type for ODE dydt functions.
51 !!
52 !! @param status Exit status
53 !! | Value | Description
54 !! |-----------|------------
55 !! | -inf-0 | Everything worked
56 !! | 1-255 | Error in this routine
57 !! @param t Value for @f$t@f$.
58 !! @param y Value for @f$\mathbf{y}@f$.
59 !! @param param Data payload usually used for constants.
60 !!
61 abstract interface
62 subroutine deq_iface(status, dydt, t, y, param)
63 use mrkiss_config, only: rk
64 implicit none
65 integer, intent(out) :: status
66 real(kind=rk), intent(out) :: dydt(:)
67 real(kind=rk), intent(in) :: t
68 real(kind=rk), intent(in) :: y(:)
69 real(kind=rk), intent(in) :: param(:)
70 end subroutine deq_iface
71 end interface
72
73 !--------------------------------------------------------------------------------------------------------------------------------
74 !> Type step processing subroutine.
75 !!
76 !! @param status Exit status
77 !! | Value | Description
78 !! |-----------|------------
79 !! | -inf-0 | Everything worked
80 !! | 256-511 | Error in this routine
81 !! @param end_run Return used to trigger the calling solver to stop the run when `end_run` is positive.
82 !! @param sdf_flags Return used to trigger SDF run in the calling solver, and pass information to `sdf_iface` function
83 !! @param new_t_delta Returns a new value for @f$\Delta{t}@f$. Used by calling solver if positive.
84 !! @param pnt_idx Index of current point in `solution`
85 !! @param solution Solution
86 !! @param t_delta Value from this step's @f$\Delta{t}@f$
87 !! @param y_delta Value from for this step's @f$\mathbf{\Delta{y}}@f$
88 !!
89 abstract interface
90 subroutine stepp_iface(status, end_run, sdf_flags, new_t_delta, pnt_idx, solution, t_delta, y_delta)
91 use mrkiss_config, only: rk
92 implicit none
93 integer, intent(out) :: status
94 integer, intent(out) :: end_run
95 real(kind=rk), intent(out) :: new_t_delta ! If >0, then redo step with new t_delta
96 integer, intent(out) :: sdf_flags ! If >0, then use bisection to solve SDF and use the result to redo step...
97 integer, intent(in) :: pnt_idx
98 real(kind=rk), intent(in) :: solution(:,:), t_delta, y_delta(:)
99 end subroutine stepp_iface
100 end interface
101
102 !--------------------------------------------------------------------------------------------------------------------------------
103 !> Type SDF function on a solution point.
104 !!
105 !! @param status Exit status
106 !! | Value | Description
107 !! |-----------|------------
108 !! | -inf-0 | Everything worked
109 !! | 512-767 | Error in this routine
110 !! @param dist The distance value of the SDF funciton
111 !! @param sdf_flags Flags passed from an `stepp_iface` routine
112 !! @param t Value for @f$t@f$.
113 !! @param y Value for @f$\mathbf{y}@f$.
114 !! @param param Data payload usually used for constants.
115 !!
116 abstract interface
117 subroutine sdf_iface(status, dist, sdf_flags, t, y)
118 use mrkiss_config, only: rk
119 implicit none
120 integer, intent(out) :: status
121 real(kind=rk), intent(out) :: dist
122 integer, intent(in) :: sdf_flags
123 real(kind=rk), intent(in) :: t, y(:)
124 end subroutine sdf_iface
125 end interface
126
127contains
128
129 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
130 !> @name One Step Solvers
131
132 !--------------------------------------------------------------------------------------------------------------------------------
133 !> Compute one step of a embedded Runge-Kutta method expressed as a Butcher Tableau.
134 !!
135 !! @param status Exit status
136 !! | Value | Description
137 !! |-----------|------------
138 !! | -inf-0 | Everything worked
139 !! | 0-255 | Evaluation of @p deq failed
140 !! | 1232-1247 | Unknown error in this routine
141 !! @param y_deltas Returned @f$\mathbf{\Delta{y}}@f$ values(s) for the method
142 !! @param dy Returned @f$\mathbf{y}'@f$ value at @f$t@f$
143 !! @param deq The equation subroutine returning values for @f$\mathbf{y}'@f$ -- i.e. @f$\mathbf{f}(t, \mathbf{y})@f$
144 !! @param t Initial condition for @f$t@f$ -- i.e. @f$t_0@f$.
145 !! @param y Initial condition for @f$\mathbf{y}@f$ -- i.e. @f$\mathbf{y_0}@f$.
146 !! @param param Data payload passed to @p deq
147 !! @param a The butcher tableau @f$\mathbf{a}@f$ matrix
148 !! @param b The butcher tableau @f$\mathbf{b}@f$ vector
149 !! @param c The butcher tableau @f$\mathbf{c}@f$ vector
150 !! @param t_delta The @f$\Delta{t}@f$ value for this step.
151 !!
152 subroutine step_all(status, y_deltas, dy, deq, t, y, param, a, b, c, t_delta)
153 use mrkiss_config, only: rk, zero_epsilon
154 implicit none
155 ! Arguments
156 integer, intent(out) :: status
157 real(kind=rk), intent(out) :: y_deltas(:,:), dy(:)
158 procedure(deq_iface) :: deq
159 real(kind=rk), intent(in) :: t
160 real(kind=rk), intent(in) :: y(:), param(:), a(:,:), b(:,:), c(:), t_delta
161 ! Variables
162 integer :: i, stage
163 real(kind=rk) :: k(size(y, 1),size(b, 1)), stage_t_delta, y_tmp(size(y, 1)), stage_y_delta(size(y, 1))
164 ! Compute k vectors
165 call deq(status, dy, t, y, param)
166 if (status > 0) return
167 k(:,1) = dy * t_delta
168 do stage=2,size(b, 1)
169 ! Same as 'stage_y_delta = matmul(k(:,1:stage-1), a(:,stage))', but faster...
170 stage_y_delta = 0.0_rk
171 do i=1,(stage-1)
172 if (abs(a(i,stage)) > zero_epsilon) then
173 stage_y_delta = stage_y_delta + a(i,stage) * k(:,i)
174 end if
175 end do
176 stage_t_delta = t_delta*c(stage)
177 call deq(status, y_tmp, t+stage_t_delta, y+stage_y_delta, param)
178 if (status > 0) return
179 k(:,stage) = y_tmp * t_delta
180 end do
181 ! Use k vectors to compute deltas
182 y_deltas = matmul(k, b)
183 status = 0
184 end subroutine step_all
185
186 !--------------------------------------------------------------------------------------------------------------------------------
187 !> Compute one step of a Runge-Kutta method expressed as a Butcher Tableau.
188 !!
189 !! @param status Exit status
190 !! | Value | Description
191 !! |-----------|------------
192 !! | -inf-0 | Everything worked
193 !! | 0-255 | Evaluation of @p deq failed
194 !! | 1232-1247 | Unknown error in this routine
195 !! @param y_delta Returned @f$\mathbf{\Delta{y}}@f$ value for the method
196 !! @param dy Returned @f$\mathbf{y}'@f$ value at @f$t@f$
197 !! @param deq The equation subroutine returning values for @f$\mathbf{y}'@f$ -- i.e. @f$\mathbf{f}(t, \mathbf{y})@f$
198 !! @param t Initial condition for @f$t@f$ -- i.e. @f$t_0@f$.
199 !! @param y Initial condition for @f$\mathbf{y}@f$ -- i.e. @f$\mathbf{y_0}@f$.
200 !! @param param Data payload passed to @p deq
201 !! @param a The butcher tableau @f$\mathbf{a}@f$ matrix
202 !! @param b The butcher tableau @f$\mathbf{b}@f$ vector
203 !! @param c The butcher tableau @f$\mathbf{c}@f$ vector
204 !! @param t_delta The @f$\Delta{t}@f$ value for this step.
205 !!
206 ! SHELLO: sed -n '/^ *subroutine step_all(/,/end subroutine step_all *$/p' mrkiss_solvers_wt.f90 | sed 's/k, b/k, b(:,1)/; s/y_deltas(:,:)/y_delta(:)/; s/y_deltas/y_delta/; s/_all/_one/;'
207 subroutine step_one(status, y_delta, dy, deq, t, y, param, a, b, c, t_delta)
208 use mrkiss_config, only: rk, zero_epsilon
209 implicit none
210 ! Arguments
211 integer, intent(out) :: status
212 real(kind=rk), intent(out) :: y_delta(:), dy(:)
213 procedure(deq_iface) :: deq
214 real(kind=rk), intent(in) :: t
215 real(kind=rk), intent(in) :: y(:), param(:), a(:,:), b(:,:), c(:), t_delta
216 ! Variables
217 integer :: i, stage
218 real(kind=rk) :: k(size(y, 1),size(b, 1)), stage_t_delta, y_tmp(size(y, 1)), stage_y_delta(size(y, 1))
219 ! Compute k vectors
220 call deq(status, dy, t, y, param)
221 if (status > 0) return
222 k(:,1) = dy * t_delta
223 do stage=2,size(b, 1)
224 ! Same as 'stage_y_delta = matmul(k(:,1:stage-1), a(:,stage))', but faster...
225 stage_y_delta = 0.0_rk
226 do i=1,(stage-1)
227 if (abs(a(i,stage)) > zero_epsilon) then
228 stage_y_delta = stage_y_delta + a(i,stage) * k(:,i)
229 end if
230 end do
231 stage_t_delta = t_delta*c(stage)
232 call deq(status, y_tmp, t+stage_t_delta, y+stage_y_delta, param)
233 if (status > 0) return
234 k(:,stage) = y_tmp * t_delta
235 end do
236 ! Use k vectors to compute deltas
237 y_delta = matmul(k, b(:,1))
238 status = 0
239 end subroutine step_one
240
241
242 !--------------------------------------------------------------------------------------------------------------------------------
243 !> Compute one Richardson Extrapolation Step.
244 !!
245 !! This routine produces an improved value for @f$\mathbf{\Delta{y}}@f$, whch we denote by @f$\mathbf{\Delta{\tilde{y}}}@f$, by
246 !! combining two approximations of @f$\mathbf{\Delta{y}}@f$. If the original Runge-Kutta method was @f$\mathcal{O}(p)@f$, then
247 !! this new approximation is @f$\mathcal{O}(p+1)@f$. The first approximation of @f$\mathbf{\Delta{y}}@f$, which we denote by
248 !! @f$\mathbf{\Delta{y_B}}@f$, is computed in the standard way with a single Runge-Kutta step of size @f$\Delta{t}@f$. The
249 !! second approximation of @f$\mathbf{\Delta{y}}@f$, which we denote by @f$\mathbf{\Delta{y_S}}@f$, is computed by taking two
250 !! Runge-Kutta steps of size @f$\frac{\Delta{t}}{2}@f$. We combine these approximations to compute our improved
251 !! @f$\mathbf{\Delta{y}}@f$ as follows:
252 !!
253 !! @f[ \mathbf{\Delta{\tilde{y}}} = \mathbf{\Delta{y_S}} + \frac{\mathbf{\Delta{y_S}} - \mathbf{\Delta{y_B}}}{2^p - 1} @f]
254 !!
255 !! @param status Exit status
256 !! | Value | Description
257 !! |-----------|------------
258 !! | -inf-0 | Everything worked
259 !! | 0-255 | Evaluation of @p deq failed
260 !! | 1216-1231 | Unknown error in this routine
261 !! | 1248-1263 | Error from mrkiss_solvers_wt::step_one()
262 !! @param y_delta Returned @f$\mathbf{\Delta{y}}@f$ for the method
263 !! @param dy Returned @f$\mathbf{y}'@f$ value at @f$t@f$
264 !! @param deq The equation subroutine returning values for @f$\mathbf{y}'@f$ -- i.e. @f$\mathbf{f}(t, \mathbf{y})@f$
265 !! @param t Initial condition for @f$t@f$ -- i.e. @f$t_0@f$.
266 !! @param y Initial condition for @f$\mathbf{y}@f$ -- i.e. @f$\mathbf{y_0}@f$.
267 !! @param param Data payload passed to @p deq
268 !! @param a The butcher tableau @f$\mathbf{a}@f$ matrix
269 !! @param b The butcher tableau @f$\mathbf{{b}}@f$ vector
270 !! @param c The butcher tableau @f$\mathbf{c}@f$ vector
271 !! @param p The order(s) of the Runge-Kutta methods in the butcher tableau
272 !! @param t_delta The @f$\Delta{t}@f$ value for this step.
273 !!
274 subroutine step_richardson(status, y_delta, dy, deq, t, y, param, a, b, c, p, t_delta)
275 use mrkiss_config, only: rk
276 implicit none
277 ! Arguments
278 integer, intent(out) :: status
279 real(kind=rk), intent(out) :: y_delta(:), dy(:)
280 procedure(deq_iface) :: deq
281 real(kind=rk), intent(in) :: t
282 real(kind=rk), intent(in) :: y(:), param(:), a(:,:), b(:,:), c(:), t_delta
283 integer, intent(in) :: p(:)
284 ! Variables
285 real(kind=rk) :: t_delta_tmp, t_tmp, y_tmp(size(y, 1))
286 real(kind=rk) :: y_delta_small_1(size(y, 1)), y_delta_small_2(size(y, 1)), y_delta_big(size(y, 1))
287 ! Compute y_delta_small
288 t_delta_tmp = t_delta / 2.0_rk
289 t_tmp = t
290 y_tmp = y
291 call step_one(status, y_delta_small_1, dy, deq, t_tmp, y_tmp, param, a, b, c, t_delta_tmp)
292 if (status > 0) return
293 t_tmp = t_tmp + t_delta_tmp
294 y_tmp = y_tmp + y_delta_small_1
295 call step_one(status, y_delta_small_2, dy, deq, t_tmp, y_tmp, param, a, b, c, t_delta_tmp)
296 if (status > 0) return
297 ! Compute y_delta_big
298 call step_one(status, y_delta_big, dy, deq, t, y, param, a, b, c, t_delta)
299 if (status > 0) return
300 ! Compute y_delta
301 y_delta = y_delta_small_1 + y_delta_small_2 + (y_delta_small_1 + y_delta_small_2 - y_delta_big) / (2**p(1) - 1)
302 status = 0
303 end subroutine step_richardson
304
305 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
306 !> @name One Step Test Solvers
307
308 !--------------------------------------------------------------------------------------------------------------------------------
309 !> Compute one step of RK4 -- the classical 4th order Runge-Kutta method (mrkiss_erk_kutta_4).
310 !!
311 !! @param status Exit status
312 !! | Value | Description
313 !! |-----------|------------
314 !! | -inf-0 | Everything worked
315 !! | 0-255 | Evaluation of @p deq failed
316 !! | 1200-1215 | Unknown error in this routine
317 !! @param y_delta Returned @f$\mathbf{\Delta{y}}@f$ for the method
318 !! @param dy Returned @f$\mathbf{y}'@f$ value at @f$t@f$
319 !! @param deq The equation subroutine returning values for @f$\mathbf{y}'@f$ -- i.e. @f$\mathbf{f}(t, \mathbf{y})@f$
320 !! @param t Initial condition for @f$t@f$ -- i.e. @f$t_0@f$.
321 !! @param y Initial condition for @f$\mathbf{y}@f$ -- i.e. @f$\mathbf{y_0}@f$.
322 !! @param param Data payload passed to @p deq
323 !! @param t_delta The @f$\Delta{t}@f$ value for this step.
324 !!
325 subroutine step_rk4(status, y_delta, dy, deq, t, y, param, t_delta)
326 use mrkiss_config, only: rk
327 implicit none
328 ! Arguments
329 integer, intent(out) :: status
330 real(kind=rk), intent(out) :: y_delta(:), dy(:)
331 procedure(deq_iface) :: deq
332 real(kind=rk), intent(in) :: t
333 real(kind=rk), intent(in) :: y(:), param(:), t_delta
334 ! Variables
335 real(kind=rk) :: k1(size(y, 1)), k2(size(y, 1)), k3(size(y, 1)), k4(size(y, 1))
336 ! Compute Step
337 call deq(status, k1, t, y, param)
338 if (status > 0) return
339 dy = k1
340 call deq(status, k2, t + 1.0_rk/2.0_rk * t_delta, y + 1.0_rk/2.0_rk * t_delta * k1, param)
341 if (status > 0) return
342 call deq(status, k3, t + 1.0_rk/2.0_rk * t_delta, y + 1.0_rk/2.0_rk * t_delta * k2, param)
343 if (status > 0) return
344 call deq(status, k4, t + t_delta, y + t_delta * k3, param)
345 y_delta = t_delta * (k1 + 2.0_rk * k2 + 2.0_rk * k3 + k4) / 6.0_rk
346 status = 0
347 end subroutine step_rk4
348
349 !--------------------------------------------------------------------------------------------------------------------------------
350 !> Compute one step of RKF45 (mrkiss_eerk_fehlberg_4_5).
351 !!
352 !! @param status Exit status
353 !! | Value | Description
354 !! |-----------|------------
355 !! | -inf-0 | Everything worked
356 !! | 0-255 | Evaluation of @p deq failed
357 !! | 1184-1199 | Unknown error in this routine
358 !! @param y_deltas Returned @f$\mathbf{\Delta{y}}@f$ values
359 !! @param dy Returned @f$\mathbf{y}'@f$ value at @f$t@f$
360 !! @param deq The equation subroutine returning values for @f$\mathbf{y}'@f$ -- i.e. @f$\mathbf{f}(t, \mathbf{y})@f$
361 !! @param t Initial condition for @f$t@f$ -- i.e. @f$t_0@f$.
362 !! @param y Initial condition for @f$\mathbf{y}@f$ -- i.e. @f$\mathbf{y_0}@f$.
363 !! @param param Data payload passed to @p deq
364 !! @param t_delta The @f$\Delta{t}@f$ value for this step.
365 !!
366 subroutine step_rkf45(status, y_deltas, dy, deq, t, y, param, t_delta)
367 use mrkiss_config, only: rk
368 implicit none
369 ! Arguments
370 integer, intent(out) :: status
371 real(kind=rk), intent(out) :: y_deltas(:,:), dy(:)
372 procedure(deq_iface) :: deq
373 real(kind=rk), intent(in) :: t
374 real(kind=rk), intent(in) :: y(:), param(:), t_delta
375 ! Variables
376 real(kind=rk) :: k1(size(y, 1)), k2(size(y, 1)), k3(size(y, 1))
377 real(kind=rk) :: k4(size(y, 1)), k5(size(y, 1)), k6(size(y, 1))
378 ! Compute Step
379 call deq(status, k1, t, y, param)
380 if (status > 0) return
381 dy = k1
382 call deq(status, k2, t + t_delta*1.0_rk/4.0_rk, y + t_delta * ( k1*1.0_rk/4.0_rk), param)
383 if (status > 0) return
384 call deq(status, k3, t + t_delta*3.0_rk/8.0_rk, y + t_delta * ( k1*3.0_rk/32.0_rk + k2*9.0_rk/32.0_rk), param)
385 if (status > 0) return
386 call deq(status, k4, t + t_delta*12.0_rk/13.0_rk, y + t_delta * ( k1*1932.0_rk/2197.0_rk - k2*7200.0_rk/2197.0_rk + k3*7296.0_rk/2197.0_rk), param)
387 if (status > 0) return
388 call deq(status, k5, t + t_delta, y + t_delta * ( k1*439.0_rk/216.0_rk - k2*8.0_rk + k3*3680.0_rk/513.0_rk - k4*845.0_rk/4104.0_rk), param)
389 if (status > 0) return
390 call deq(status, k6, t + t_delta*1.0_rk/2.0_rk, y + t_delta * (-k1*8.0_rk/27.0_rk + k2*2.0_rk - k3*3544.0_rk/2565.0_rk + k4*1859.0_rk/4104.0_rk - k5*11.0_rk/40.0_rk), param)
391 if (status > 0) return
392 y_deltas(:,1) = t_delta * (k1*25.0_rk/216.0_rk + k3*1408.0_rk/2565.0_rk + k4*2197.0_rk/4104.0_rk - k5*1.0_rk/5.0_rk)
393 y_deltas(:,2) = t_delta * (k1*16.0_rk/135.0_rk + k3*6656.0_rk/12825.0_rk + k4*28561.0_rk/56430.0_rk - k5*9.0_rk/50.0_rk + k6*2.0_rk/55.0_rk)
394 status = 0
395 end subroutine step_rkf45
396
397 !--------------------------------------------------------------------------------------------------------------------------------
398 !> Compute one step of DP45 (mrkiss_eerk_dormand_prince_5_4)
399 !!
400 !! @param status Exit status
401 !! | Value | Description
402 !! |-----------|------------
403 !! | -inf-0 | Everything worked
404 !! | 0-255 | Evaluation of @p deq failed
405 !! | 1263-1279 | Unknown error in this routine
406 !! @param y_deltas Returned @f$\mathbf{\Delta{y}}@f$ values
407 !! @param dy Returned @f$\mathbf{y}'@f$ value at @f$t@f$
408 !! @param deq The equation subroutine returning values for @f$\mathbf{y}'@f$ -- i.e. @f$\mathbf{f}(t, \mathbf{y})@f$
409 !! @param t Initial condition for @f$t@f$ -- i.e. @f$t_0@f$.
410 !! @param y Initial condition for @f$\mathbf{y}@f$ -- i.e. @f$\mathbf{y_0}@f$.
411 !! @param param Data payload passed to @p deq
412 !! @param t_delta The @f$\Delta{t}@f$ value for this step.
413 !!
414 subroutine step_dp54(status, y_deltas, dy, deq, t, y, param, t_delta)
415 use mrkiss_config, only: rk
416 implicit none
417 ! Arguments
418 integer, intent(out) :: status
419 real(kind=rk), intent(out) :: y_deltas(:,:), dy(:)
420 procedure(deq_iface) :: deq
421 real(kind=rk), intent(in) :: t
422 real(kind=rk), intent(in) :: y(:), param(:), t_delta
423 ! Variables
424 real(kind=rk) :: k1(size(y, 1)), k2(size(y, 1)), k3(size(y, 1)), k4(size(y, 1))
425 real(kind=rk) :: k5(size(y, 1)), k6(size(y, 1)), k7(size(y, 1))
426 ! Compute Step
427 call deq(status, k1, t, y, param)
428 if (status > 0) return
429 dy = k1;
430 call deq(status, k2, t+t_delta*1.0_rk/5.0_rk, y + t_delta * (k1*1.0_rk/5.0_rk), param)
431 if (status > 0) return
432 call deq(status, k3, t+t_delta*3.0_rk/10.0_rk, y + t_delta * (k1*3.0_rk/40.0_rk + k2*9.0_rk/40.0_rk), param)
433 if (status > 0) return
434 call deq(status, k4, t+t_delta*4.0_rk/5.0_rk, y + t_delta * (k1*44.0_rk/45.0_rk - k2*56.0_rk/15.0_rk + k3*32.0_rk/9.0_rk), param)
435 if (status > 0) return
436 call deq(status, k5, t+t_delta*8.0_rk/9.0_rk, y + t_delta * (k1*19372.0_rk/6561.0_rk - k2*25360.0_rk/2187.0_rk + k3*64448.0_rk/6561.0_rk - k4*212.0_rk/729.0_rk), param)
437 if (status > 0) return
438 call deq(status, k6, t+t_delta, y + t_delta * (k1*9017.0_rk/3168.0_rk - k2*355.0_rk/33.0_rk + k3*46732.0_rk/5247.0_rk + k4*49.0_rk/176.0_rk - k5*5103.0_rk/18656.0_rk), param)
439 if (status > 0) return
440 call deq(status, k7, t+t_delta, y + t_delta * (k1*35.0_rk/384.0_rk + k3*500.0_rk/1113.0_rk + k4*125.0_rk/192.0_rk - k5*2187.0_rk/6784.0_rk + k6*11.0_rk/84.0_rk), param)
441 if (status > 0) return
442 y_deltas(:,1) = t_delta * (k1*35.0_rk/384.0_rk + k3*500.0_rk/1113.0_rk + k4*125.0_rk/192.0_rk - k5*2187.0_rk/6784.0_rk + k6*11.0_rk/84.0_rk)
443 y_deltas(:,2) = t_delta * (k1*5179.0_rk/57600.0_rk + k3*7571.0_rk/16695.0_rk + k4*393.0_rk/640.0_rk - k5*92097.0_rk/339200.0_rk + k6*187.0_rk/2100.0_rk + k7*1.0_rk/40.0_rk)
444 status = 0
445 end subroutine step_dp54
446
447 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
448 !> @name Multistep Solvers
449
450 !--------------------------------------------------------------------------------------------------------------------------------
451 !> Take multiple fixed time steps with a simple Runge-Kutta method.
452 !!
453 !! @param status Exit status
454 !! | Value | Description
455 !! |-----------|------------
456 !! | -inf-0 | Everything worked
457 !! | 0-255 | Evaluation of @p deq failed
458 !! | 1120-1151 | Unnecessary error in this routine
459 !! | 1216-1231 | Error from mrkiss_solvers_wt::step_richardson()
460 !! | 1248-1263 | Error from mrkiss_solvers_wt::step_one()
461 !! @param istats Integer statistics for run
462 !! - See: mrkiss_config::istats_strs for a description of elements.
463 !! - Elements this routine updates
464 !! - mrkiss_config::isi_num_pts
465 !! - mrkiss_config::isi_one_reg
466 !! @param solution Array for solution.
467 !! Each COLUMN is a solution:
468 !! - First element is the @f$t@f$ variable
469 !! - `size(y, 1)` elements starting with 2 have @f$\mathbf{y}@f$ values
470 !! - The next `size(y, 1)` elements have @f$\mathbf{y}'@f$ values
471 !! @param deq The equation subroutine returning values for @f$\mathbf{y}'@f$ -- i.e. @f$\mathbf{f}(t, \mathbf{y})@f$
472 !! @param t Initial condition for @f$t@f$ -- i.e. @f$t_0@f$.
473 !! @param y Initial condition for @f$\mathbf{y}@f$ -- i.e. @f$\mathbf{y_0}@f$.
474 !! @param param Data payload passed to @p deq
475 !! @param a The butcher tableau @f$\mathbf{a}@f$ matrix
476 !! @param b The butcher tableau @f$\mathbf{{b}}@f$ vector
477 !! @param c The butcher tableau @f$\mathbf{c}@f$ vector
478 !! @param p_o The order for the Runge-Kutta method in the butcher tableau to enable Richardson extrapolation
479 !! @param max_pts_o Maximum number of points to put in @p solution.
480 !! If `max_pts_o>1` & `size(solution, 2)==1`, then `max_pts_o-1` steps are taken and
481 !! only the last solution is stored in solution.
482 !! @param t_delta_o Step size to use.
483 !! - Default when @p t_end_o provided: `(t_end - t) / (size(solution, 2) - 1)`
484 !! - Default othewise: mrkiss_config::t_delta_ai
485 !! @param t_end_o End point for last step. Silently ignored if @p t_delta_o is provided.
486 !! @param t_max_o Maximum value for @f$t@f$
487 !!
488 subroutine fixed_t_steps(status, istats, solution, deq, t, y, param, a, b, c, p_o, max_pts_o, t_delta_o, t_end_o, t_max_o)
489 use mrkiss_config, only: rk, t_delta_ai, istats_size, isi_num_pts, isi_one_reg
490 implicit none
491 ! Arguments
492 integer, intent(out) :: status, istats(istats_size)
493 real(kind=rk), intent(out) :: solution(:,:)
494 procedure(deq_iface) :: deq
495 real(kind=rk), intent(in) :: t
496 real(kind=rk), intent(in) :: y(:), param(:), a(:,:), b(:,:), c(:)
497 integer, optional, intent(in) :: p_o(:), max_pts_o
498 real(kind=rk), optional, intent(in) :: t_delta_o, t_end_o, t_max_o
499 ! Vars
500 integer :: cur_pnt_idx, y_dim, cur_step, max_steps, p(size(b, 2))
501 real(kind=rk) :: t_cv, t_delta, y_cv(size(y, 1)), y_delta(size(y, 1)), dy(size(y, 1))
502 logical :: lotsopnts
503 ! Process arguments
504 max_steps = size(solution, 2) - 1
505 lotsopnts = .true.
506 if (present(max_pts_o)) then
507 if ((max_steps == 0) .and. (max_pts_o > 1)) then
508 lotsopnts = .false.
509 max_steps = max_pts_o - 1
510 else
511 max_steps = min(max_pts_o - 1, max_steps)
512 end if
513 end if
514 if (present(t_delta_o)) then
515 t_delta = t_delta_o
516 else
517 if (present(t_end_o)) then
518 t_delta = (t_end_o - t) / max_steps
519 else
520 t_delta = t_delta_ai
521 end if
522 end if
523 p = 0
524 if (present(p_o)) p = p_o
525 ! Compute Solution
526 y_dim = size(y, 1)
527 istats = 0
528 t_cv = t
529 y_cv = y
530 cur_step = 0
531 cur_pnt_idx = 1
532 solution(1, cur_pnt_idx) = t_cv
533 solution(2:(2+y_dim-1), cur_pnt_idx) = y_cv
534 do
535 cur_step = cur_step + 1
536 if (p(1) > 0) then
537 call step_richardson(status, y_delta, dy, deq, t_cv, y_cv, param, a, b, c, p, t_delta)
538 istats(isi_one_reg) = istats(isi_one_reg) + 3
539 else
540 call step_one(status, y_delta, dy, deq, t_cv, y_cv, param, a, b, c, t_delta)
541 istats(isi_one_reg) = istats(isi_one_reg) + 1
542 end if
543 if (status > 0) return
544 y_cv = y_cv + y_delta
545 t_cv = t_cv + t_delta
546 if (lotsopnts) then
547 cur_pnt_idx = cur_pnt_idx + 1
548 solution((2+y_dim):(2+2*y_dim-1), cur_pnt_idx-1) = dy
549 end if
550 solution(1, cur_pnt_idx) = t_cv
551 solution(2:(2+y_dim-1), cur_pnt_idx) = y_cv
552 istats(isi_num_pts) = istats(isi_num_pts) + 1
553 status = 0
554 if (present(t_max_o)) then
555 if (t_cv > t_max_o) exit
556 end if
557 if (cur_step >= max_steps) exit
558 end do
559 ! Compute derivative for final solution point
560 call deq(status, dy, t_cv, y_cv, param) ! This sets return status
561 if (status > 0) return
562 solution((2+y_dim):(2+2*y_dim-1), cur_pnt_idx) = dy
563 istats(isi_num_pts) = istats(isi_num_pts) + 1
564 end subroutine fixed_t_steps
565
566 !--------------------------------------------------------------------------------------------------------------------------------
567 !> Take multiple adaptive steps with constant @f$\mathbf{\Delta{y}}@f$ length using a simple Runge-Kutta method.
568 !!
569 !! This method attempts to precisely control the length of @f$\mathbf{\Delta{y}}@f$ which we denote by
570 !! @f$\vert\vert\mathbf{\Delta{y}}\vert\vert_L@f$. This is done by finding a value for @f$\Delta{t}@f$ for which
571 !! @f$\vert\vert\mathbf{\Delta{y}}\vert\vert_L@f$ is equal to @p y_delta_len_targ with an accuracy of @p y_delta_len_tol_o.
572 !!
573 !! This value for @f$\Delta{t}@f$ is found via bisection. For each solution point we first compute two Runge-Kutta steps with
574 !! @f$\Delta{t}@f$ set to @p t_delta_min_o and @p t_delta_max. We hope that the lengths of the two @f$\mathbf{\Delta{y}}@f$
575 !! values thus obtained will bracket @p y_delta_len_targ. If they don't then we can't isolate an appropriate value for
576 !! @f$\Delta{t}@f$. In this case we ignore the issue and continue to the next step if @p no_bisect_error_o is .TRUE.,
577 !! otherwise we error out with a status code of 1024. If the lengths of our @f$\mathbf{\Delta{y}}@f$ values bracket the
578 !! target, then we use bisection to localize a value for @f$\Delta{t}@f$.
579 !!
580 !! The length, which we have denoted by @f$\vert\vert\mathbf{\Delta{y}}\vert\vert_L@f$, is defined to be the Euclidean
581 !! distance of the vector defined by the elements of @f$\mathbf{\Delta{y}}@f$ that are indexed by @p y_delta_len_idxs_o.
582 !!
583 !! Note there is no mathematical guarantee that a Runge-Kutta step of size @p t_delta_min_o and @p t_delta_max will produce
584 !! solutions that bracket @p y_delta_len_targ. That said, for well behaved functions a @p t_delta_min_o and @p t_delta_max may
585 !! always be found that work. Usually finding good values are @p t_delta_min_o and @p t_delta_max isn't difficult. For
586 !! challenging cases, a good approach is to use fixed_t_steps() over the interval in question, and then examine the values
587 !! of @f$\vert\vert\mathbf{\Delta{y}}\vert\vert_L@f$ in the solution via analyze_solution().
588 !!
589 !! My primary use case for this function is to create uniform sphere sweeps for constructive solid geometry applications.
590 !!
591 !! @param status Exit status
592 !! | Value | Description
593 !! |-----------|------------
594 !! | 0 | Everything worked
595 !! | 0-255 | Evaluation of @p deq failed (see: deq_iface)
596 !! | 1024 | bisection containment anomaly
597 !! | 1025 | `no_bisect_error_o==0` not present and `max_bisect_o` violated
598 !! | 1026-1055 | Unknown error in this routine
599 !! | 1248-1263 | Error from mrkiss_solvers_wt::step_one()
600 !! @param istats Integer statistics for run
601 !! - See: mrkiss_config::istats_strs for a description of elements.
602 !! - Elements this routine updates
603 !! - mrkiss_config::isi_num_pts
604 !! - mrkiss_config::isi_one_reg
605 !! - mrkiss_config::isi_one_y_len
606 !! - mrkiss_config::isi_bic_max
607 !! - mrkiss_config::isi_bic_bnd
608 !! @param solution Array for solution.
609 !! Each COLUMN is a solution:
610 !! - First element is the @f$t@f$ variable
611 !! - `size(y, 1)` elements starting with 2 have @f$\mathbf{y}@f$ values
612 !! - The next `size(y, 1)` elements have dy values if
613 !! @param deq The equation subroutine returning values for @f$\mathbf{y}'@f$ -- i.e. @f$\mathbf{f}(t, \mathbf{y})@f$
614 !! @param t Initial condition for @f$t@f$ -- i.e. @f$t_0@f$.
615 !! @param y Initial condition for @f$\mathbf{y}@f$ -- i.e. @f$\mathbf{y_0}@f$.
616 !! @param param Data payload passed to @p deq
617 !! @param a The butcher tableau @f$\mathbf{a}@f$ matrix
618 !! @param b The butcher tableau @f$\mathbf{{b}}@f$ vector
619 !! @param c The butcher tableau @f$\mathbf{c}@f$ vector
620 !! @param y_delta_len_targ Target length for @f$\mathbf{\Delta{y}}@f$ -- i.e. @f$\vert\vert\mathbf{\Delta{y}}\vert\vert_L@f$
621 !! @param t_delta_max Maximum @f$\Delta{t}@f$
622 !! @param t_delta_min_o Minimum @f$\Delta{t}@f$
623 !! @param y_delta_len_tol_o How close we have to get to @p y_delta_len_targ. Default is `y_delta_len_targ/100`
624 !! @param y_delta_len_idxs_o Components of @f$\mathbf{\Delta{y}}@f$ to use for length computation
625 !! @param max_pts_o Maximum number of points to put in @p solution.
626 !! @param max_bisect_o Maximum number of bisection iterations per each step. Default: mrkiss_config::max_bisect_ai
627 !! @param no_bisect_error_o If `.TRUE.`, then do not exit on bisection errors
628 !! @param y_sol_len_max_o Maximum length of the solution curve
629 !! @param t_max_o Maximum value for @f$t@f$
630 !!
631 subroutine fixed_y_steps(status, istats, solution, deq, t, y, param, a, b, c, y_delta_len_targ, t_delta_max, t_delta_min_o, &
632 & y_delta_len_tol_o, max_bisect_o, no_bisect_error_o, y_delta_len_idxs_o, max_pts_o, y_sol_len_max_o, t_max_o)
633 use mrkiss_config, only: rk, t_delta_min_ai, max_bisect_ai, istats_size, isi_bic_bnd, isi_bic_max, isi_num_pts, isi_one_reg, isi_one_y_len
634 implicit none
635 ! Arguments
636 integer, intent(out) :: status, istats(istats_size)
637 real(kind=rk), intent(out) :: solution(:,:)
638 procedure(deq_iface) :: deq
639 real(kind=rk), intent(in) :: t
640 real(kind=rk), intent(in) :: y(:), param(:), a(:,:), b(:,:), c(:), y_delta_len_targ, t_delta_max
641 real(kind=rk), optional, intent(in) :: t_delta_min_o, y_delta_len_tol_o, y_sol_len_max_o, t_max_o
642 integer, optional, intent(in) :: max_pts_o, max_bisect_o, y_delta_len_idxs_o(:)
643 logical, optional, intent(in) :: no_bisect_error_o
644 ! Variables
645 integer :: max_bisect, max_pts, cur_pnt_idx, biter, y_dim
646 logical :: no_bisect_error
647 real(kind=rk) :: y_delta_len_tol, t_delta_min, y_sol_len, bs_tmp1_y_delta_len
648 real(kind=rk) :: bs_tmp1_t_delta, bs_tmp2_t_delta, t_cv, dy(size(y, 1))
649 real(kind=rk) :: bs_tmpc_dy(size(y, 1)), bs_tmp1_dy(size(y, 1)), bs_tmp2_dy(size(y, 1))
650 real(kind=rk) :: bs_tmp2_y_delta_len, bs_tmpc_y_delta_len, bs_tmpc_t_delta
651 real(kind=rk) :: y_cv(size(y, 1)), bs_tmp1_y_delta(size(y, 1))
652 real(kind=rk) :: bs_tmp2_y_delta(size(y, 1)), bs_tmpc_y_delta(size(y, 1))
653 ! Process arguments
654 max_pts = size(solution, 2)
655 if (present(max_pts_o)) max_pts = min(max_pts, max_pts_o);
656 t_delta_min = t_delta_min_ai
657 if (present(t_delta_min_o)) t_delta_min = t_delta_min_o
658 max_bisect = max_bisect_ai
659 if (present(max_bisect_o)) max_bisect = max_bisect_o
660 y_delta_len_tol = y_delta_len_targ / 100.0_rk
661 if (present(y_delta_len_tol_o)) y_delta_len_tol = y_delta_len_tol_o
662 max_pts = size(solution, 2)
663 if (present(max_pts_o)) max_pts = min(max_pts, max_pts_o);
664 no_bisect_error = .false.
665 if (present(no_bisect_error_o)) no_bisect_error = no_bisect_error_o
666 ! Compute Solution
667 y_dim = size(y, 1)
668 y_sol_len = 0.0_rk
669 istats = 0
670 t_cv = t
671 y_cv = y
672 cur_pnt_idx = 1
673 solution(1, cur_pnt_idx) = t_cv
674 solution(2:(2+y_dim-1), cur_pnt_idx) = y_cv
675 do
676 cur_pnt_idx = cur_pnt_idx + 1
677 ! Compute t_delta 1
678 bs_tmp1_t_delta = t_delta_min
679 call step_one(status, bs_tmp1_y_delta, bs_tmp1_dy, deq, t_cv, y_cv, param, a, b, c, bs_tmp1_t_delta)
680 istats(isi_one_reg) = istats(isi_one_reg) + 1
681 if (status > 0) return
682 if (present(y_delta_len_idxs_o)) then
683 bs_tmp1_y_delta_len = norm2(bs_tmp1_y_delta(y_delta_len_idxs_o))
684 else
685 bs_tmp1_y_delta_len = norm2(bs_tmp1_y_delta)
686 end if
687 ! Compute upper t_delta
688 bs_tmp2_t_delta = t_delta_max
689 call step_one(status, bs_tmp2_y_delta, bs_tmp2_dy, deq, t_cv, y_cv, param, a, b, c, bs_tmp2_t_delta)
690 istats(isi_one_reg) = istats(isi_one_reg) + 1
691 if (status > 0) return
692 if (present(y_delta_len_idxs_o)) then
693 bs_tmp2_y_delta_len = norm2(bs_tmp2_y_delta(y_delta_len_idxs_o))
694 else
695 bs_tmp2_y_delta_len = norm2(bs_tmp2_y_delta)
696 end if
697 ! Swap if required
698 if (bs_tmp2_y_delta_len < bs_tmp1_y_delta_len) then
699 bs_tmpc_t_delta = bs_tmp1_t_delta
700 bs_tmp1_t_delta = bs_tmp2_t_delta
701 bs_tmp2_t_delta = bs_tmpc_t_delta
702 bs_tmpc_y_delta_len = bs_tmp1_y_delta_len
703 bs_tmp1_y_delta_len = bs_tmp2_y_delta_len
704 bs_tmp2_y_delta_len = bs_tmpc_y_delta_len
705 bs_tmpc_y_delta = bs_tmp1_y_delta
706 bs_tmp1_y_delta = bs_tmp2_y_delta
707 bs_tmp2_y_delta = bs_tmpc_y_delta
708 bs_tmpc_dy = bs_tmp1_dy
709 bs_tmp1_dy = bs_tmp2_dy
710 bs_tmp2_dy = bs_tmpc_dy
711 end if
712 ! Initial "current" values for bisection and no-bisection
713 bs_tmpc_t_delta = bs_tmp2_t_delta
714 bs_tmpc_y_delta_len = bs_tmp2_y_delta_len
715 bs_tmpc_y_delta = bs_tmp2_y_delta
716 ! Bisect if required
717 if ((bs_tmp1_y_delta_len < y_delta_len_targ) .and. (bs_tmp2_y_delta_len > y_delta_len_targ)) then
718 biter = 1
719 do while (abs(bs_tmpc_y_delta_len - y_delta_len_targ) > y_delta_len_tol)
720 if (biter > max_bisect) then
721 istats(isi_bic_max) = istats(isi_bic_max) + 1
722 if (no_bisect_error) then
723 exit
724 else
725 status = 1025
726 return
727 end if
728 end if
729 bs_tmpc_t_delta = (bs_tmp1_t_delta + bs_tmp2_t_delta) / 2.0_rk
730 call step_one(status, bs_tmpc_y_delta, bs_tmpc_dy, deq, t_cv, y_cv, param, a, b, c, bs_tmpc_t_delta)
731 istats(isi_one_y_len) = istats(isi_one_y_len) + 1
732 if (status > 0) return
733 if (present(y_delta_len_idxs_o)) then
734 bs_tmpc_y_delta_len = norm2(bs_tmpc_y_delta(y_delta_len_idxs_o))
735 else
736 bs_tmpc_y_delta_len = norm2(bs_tmpc_y_delta)
737 end if
738 if (bs_tmpc_y_delta_len < y_delta_len_targ) then
739 bs_tmp1_y_delta_len = bs_tmpc_y_delta_len
740 bs_tmp1_t_delta = bs_tmpc_t_delta
741 else
742 bs_tmp2_y_delta_len = bs_tmpc_y_delta_len
743 bs_tmp2_t_delta = bs_tmpc_t_delta
744 end if
745 biter = biter + 1;
746 end do
747 else
748 istats(isi_bic_bnd) = istats(isi_bic_bnd) + 1
749 if (no_bisect_error) then
750 status = 1024
751 return
752 end if
753 end if
754 ! Update solution
755 y_cv = y_cv + bs_tmpc_y_delta
756 t_cv = t_cv + bs_tmpc_t_delta
757 solution(1, cur_pnt_idx) = t_cv
758 solution((2+y_dim):(2+2*y_dim-1), cur_pnt_idx-1) = bs_tmpc_dy
759 solution(2:(2+y_dim-1), cur_pnt_idx) = y_cv
760 istats(isi_num_pts) = istats(isi_num_pts) + 1
761 status = 0
762 if (present(y_sol_len_max_o)) then
763 y_sol_len = y_sol_len + bs_tmpc_y_delta_len
764 if (y_sol_len > y_sol_len_max_o) exit
765 end if
766 if (present(t_max_o)) then
767 if (t_cv > t_max_o) exit
768 end if
769 if (cur_pnt_idx >= max_pts) exit
770 end do
771 ! Compute derivative for final solution point
772 call deq(status, dy, t_cv, y_cv, param) ! This sets return status
773 if (status > 0) return
774 solution((2+y_dim):(2+2*y_dim-1), cur_pnt_idx) = dy
775 istats(isi_num_pts) = istats(isi_num_pts) + 1
776 end subroutine fixed_y_steps
777
778 !--------------------------------------------------------------------------------------------------------------------------------
779 !> Take multiple adaptive steps adjusting for the length of @f$\mathbf{\Delta{y}}@f$ using a simple Runge-Kutta method.
780 !!
781 !! This method attempts to control the length of @f$\mathbf{\Delta{y}}@f$ which we denote by
782 !! @f$\vert\vert\mathbf{\Delta{y}}\vert\vert_L@f$. At each solution step it takes a probing step with @f$\Delta{t}@f$ equal to
783 !! @p t_delta_ini and then takes another step with @f$\Delta{t}@f$ equal to: `t_delta_ini * y_delta_len_targ / y_delta_len`
784 !!
785 !! If which will result in a @f$\vert\vert\mathbf{\Delta{y}}\vert\vert_L@f$ approximately @p y_delta_len_targ when
786 !! @f$\Delta{t}@f$ is proportional to @f$\vert\vert\mathbf{\Delta{y}}\vert\vert_L@f$. By default this second step is only
787 !! taken when @f$\vert\vert\mathbf{\Delta{y}}\vert\vert_L@f$ of the probe step is greater than @p y_delta_len_targ; however, it
788 !! will be performed on shorter steps when @p adj_short_o is present -- in this mode it approximates the behavior
789 !! @p fixed_y_steps() but is *much* faster.
790 !!
791 !! Note the assumption that @f$\Delta{t}@f$ is proportional to @f$\vert\vert\mathbf{\Delta{y}}\vert\vert_L@f$. We have no
792 !! mathematical guarantee for this assumption; however, it is a fair approximation with well behaved functions when
793 !! @p t_delta_ini is small enough.
794 !!
795 !! The length, which we have denoted by @f$\vert\vert\mathbf{\Delta{y}}\vert\vert_L@f$, is defined to be the Euclidean
796 !! distance of the vector defined by the elements of @f$\mathbf{\Delta{y}}@f$ that are indexed by @p y_delta_len_idxs_o.
797 !!
798 !! My primary use case for this function is to shrink down long steps for smoother curves/tubes in visualizations.
799 !!
800 !! @param status Exit status
801 !! | Value | Description
802 !! |-----------|------------
803 !! | -inf-0 | Everything worked
804 !! | 0-255 | Evaluation of @p deq failed
805 !! | 1280-1296 | Unknown error in this routine
806 !! | 1248-1263 | Error from mrkiss_solvers_wt::step_one()
807 !! @param istats Integer statistics for run
808 !! - See: mrkiss_config::istats_strs for a description of elements.
809 !! - Elements this routine updates:
810 !! - mrkiss_config::isi_num_pts
811 !! - mrkiss_config::isi_one_reg
812 !! - mrkiss_config::isi_one_y_len
813 !! @param solution Array for solution.
814 !! Each COLUMN is a solution:
815 !! - First element is the @f$t@f$ variable
816 !! - `size(y, 1)` elements starting with 2 have @f$\mathbf{y}@f$ values
817 !! - The next `size(y, 1)` elements have @f$\mathbf{y}'@f$ values
818 !! @param deq The equation subroutine returning values for @f$\mathbf{y}'@f$ -- i.e. @f$\mathbf{f}(t, \mathbf{y})@f$
819 !! @param t Initial condition for @f$t@f$ -- i.e. @f$t_0@f$.
820 !! @param y Initial condition for @f$\mathbf{y}@f$ -- i.e. @f$\mathbf{y_0}@f$.
821 !! @param param Data payload passed to @p deq
822 !! @param a The butcher tableau @f$\mathbf{a}@f$ matrix
823 !! @param b The butcher tableau @f$\mathbf{{b}}@f$ vector
824 !! @param c The butcher tableau @f$\mathbf{c}@f$ vector
825 !! @param y_delta_len_targ Target length for @f$\mathbf{\Delta{y}}@f$ -- i.e. @f$\vert\vert\mathbf{\Delta{y}}\vert\vert_L@f$
826 !! @param t_delta_ini Test step @f$\Delta{t}@f$
827 !! @param t_delta_max_o Maximum @f$\Delta{t}@f$
828 !! @param t_delta_min_o Minimum @f$\Delta{t}@f$
829 !! @param y_delta_len_idxs_o Components of @f$\mathbf{\Delta{y}}@f$ to use for length computation
830 !! @param adj_short_o Adjust when @f$\Delta{t}@f$ is too short as well as when it's too long.
831 !! @param max_pts_o Maximum number of points to put in @p solution.
832 !! @param y_sol_len_max_o Maximum length of the solution curve
833 !! @param t_max_o Maximum value for @f$t@f$
834 !!
835 subroutine sloppy_fixed_y_steps(status, istats, solution, deq, t, y, param, a, b, c, y_delta_len_targ, t_delta_ini, t_max_o, &
836 & t_delta_min_o, t_delta_max_o, y_delta_len_idxs_o, adj_short_o, max_pts_o, y_sol_len_max_o)
837 use mrkiss_config, only: rk, t_delta_min_ai, istats_size, isi_num_pts, isi_one_reg, isi_one_y_len
838 implicit none
839 ! Arguments
840 integer, intent(out) :: status, istats(istats_size)
841 real(kind=rk), intent(out) :: solution(:,:)
842 procedure(deq_iface) :: deq
843 real(kind=rk), intent(in) :: t
844 real(kind=rk), intent(in) :: y(:), param(:), a(:,:), b(:,:), c(:), y_delta_len_targ, t_delta_ini
845 real(kind=rk), optional, intent(in) :: t_delta_min_o, t_delta_max_o, y_sol_len_max_o, t_max_o
846 integer, optional, intent(in) :: max_pts_o, y_delta_len_idxs_o(:), adj_short_o
847 ! Variables
848 integer :: max_pts, cur_pnt_idx, y_dim
849 real(kind=rk) :: t_delta_min, y_sol_len, t_cv, t_delta, y_delta_len
850 real(kind=rk) :: y_cv(size(y, 1)), y_delta(size(y, 1)), dy(size(y, 1))
851 ! Process arguments
852 max_pts = size(solution, 2)
853 if (present(max_pts_o)) max_pts = min(max_pts, max_pts_o);
854 t_delta_min = t_delta_min_ai
855 if (present(t_delta_min_o)) t_delta_min = t_delta_min_o
856 ! Compute Solution
857 y_dim = size(y, 1)
858 y_sol_len = 0.0_rk
859 istats = 0
860 t_cv = t
861 y_cv = y
862 cur_pnt_idx = 1
863 solution(1, cur_pnt_idx) = t_cv
864 solution(2:(2+y_dim-1), cur_pnt_idx) = y_cv
865 do
866 cur_pnt_idx = cur_pnt_idx + 1
867 ! Compute Initial step
868 t_delta = t_delta_ini
869 call step_one(status, y_delta, dy, deq, t_cv, y_cv, param, a, b, c, t_delta)
870 istats(isi_one_reg) = istats(isi_one_reg) + 1
871 if (status > 0) return
872 if (present(y_delta_len_idxs_o)) then
873 y_delta_len = norm2(y_delta(y_delta_len_idxs_o))
874 else
875 y_delta_len = norm2(y_delta)
876 end if
877 if ((y_delta_len > y_delta_len_targ) .or. present(adj_short_o)) then
878 ! Adjust step
879 t_delta = t_delta_ini * y_delta_len_targ / y_delta_len
880 t_delta = max(t_delta_min, t_delta)
881 if (present(t_delta_max_o)) then
882 t_delta = min(t_delta_max_o, t_delta)
883 end if
884 ! Compute adjusted step
885 call step_one(status, y_delta, dy, deq, t_cv, y_cv, param, a, b, c, t_delta)
886 istats(isi_one_y_len) = istats(isi_one_y_len) + 1
887 if (status > 0) return
888 end if
889 ! Update state
890 y_cv = y_cv + y_delta
891 t_cv = t_cv + t_delta
892 solution(1, cur_pnt_idx) = t_cv
893 solution((2+y_dim):(2+2*y_dim-1), cur_pnt_idx-1) = dy
894 solution(2:(2+y_dim-1), cur_pnt_idx) = y_cv
895 istats(isi_num_pts) = istats(isi_num_pts) + 1
896 status = 0
897 ! Process solution length limit
898 if (present(y_sol_len_max_o)) then
899 if (present(y_delta_len_idxs_o)) then
900 y_sol_len = y_sol_len + norm2(y_delta(y_delta_len_idxs_o))
901 else
902 y_sol_len = y_sol_len + norm2(y_delta)
903 end if
904 if (y_sol_len > y_sol_len_max_o) exit
905 end if
906 ! Process t size limit
907 if (present(t_max_o)) then
908 if (t_cv > t_max_o) exit
909 end if
910 ! Max points limit
911 if (cur_pnt_idx >= max_pts) exit
912 end do
913 ! Compute derivative for final solution point
914 call deq(status, dy, t_cv, y_cv, param) ! This sets return status
915 if (status > 0) return
916 solution((2+y_dim):(2+2*y_dim-1), cur_pnt_idx) = dy
917 istats(isi_num_pts) = istats(isi_num_pts) + 1
918 end subroutine sloppy_fixed_y_steps
919
920 !--------------------------------------------------------------------------------------------------------------------------------
921 !> Take multiple adaptive steps with an embedded Runge-Kutta method using relatively traditional step size controls.
922 !!
923 !! Method of controlling step size:
924 !!
925 !! First we compute a combined tolerance vector:
926 !! @f[ \mathbf{E} = [ A_i + R_i \max(\vert y_i\vert, \vert y_i+\Delta t\vert) ] @f]
927 !! Now define a set containing the indexes of the non-zero elements of @f$\mathbf{E}@f$:
928 !! @f[ \mathbf{E_+} = \{i\vert\,E_i\ne0\} @f]
929 !! When @f$\vert E_+\vert=0@f$, i.e. when @f$E_+=\emptyset@f$:
930 !! - We accept the current step
931 !! - Expand the next step by a factor of @p t_delta_fac_max_o
932 !!
933 !! Otherwise, we compute a composite error estimate:
934 !! @f[ \epsilon = \sqrt{\frac{1}{\vert E_+\vert} \sum_{E_+}\left(\frac{\vert\Delta{y}_i-\Delta{y}_i\vert}{E_i}\right)^2} @f]
935 !! From this we compute the ideal step size adjustment ratio:
936 !! @f[ m = \left(\frac{1}{\epsilon}\right)^\frac{1}{1+\min(\p, \check{p})} @f]
937 !! When @f$ \vert\Delta{y}_i-\Delta{y}_i\vert < E_i\,\,\,\,\forall i@f$, we:
938 !! - We accept the current step
939 !! - Expand the next step by a factor of @f$m@f$
940 !!
941 !! Otherwise
942 !! - We recompute the current step with a @f$\Delta{t}@f$ reduced by @f$m@f$
943 !! - Recompute the error conditions
944 !! - Set the next step size based on the error conditions
945 !!
946 !! Notes:
947 !! - When @f$m<1@f$, the factor is adjusted by @p t_delta_fac_fdg_o
948 !! - The final value for @f$m@f$ is constrained by @p t_delta_fac_max_o & @p t_delta_fac_min_o.
949 !! - The final value for @f$\Delta{t}@f$ is always constrained by @p t_delta_min_o & @p t_delta_max_o.
950 !! - This routine is sensitive to compiler optimizations. For example, ifx requires -fp-model=precise in order to insure
951 !! consistent results with other compilers -- the results without this option are not incorrect, just different.
952 !!
953 !! @param status Exit status
954 !! | Value | Description
955 !! |-----------|------------
956 !! | -inf-0 | Everything worked
957 !! | 0-255 | Evaluation of @p deq failed
958 !! | 256-511 | Error in @p stepp_o
959 !! | 512-767 | Error in @p sdf_o
960 !! | 1056 | @p no_bisect_error_o is `.FALSE`. and @p max_bisect_o violated.
961 !! | 1057-1119 | Unknown error in this routine
962 !! | 1232-1247 | Error from mrkiss_solvers_wt::step_all()
963 !! | 1248-1263 | Error from mrkiss_solvers_wt::step_one()
964 !! @param istats Integer statistics for run
965 !! - See: mrkiss_config::istats_strs for a description of elements.
966 !! - Elements this routine updates:
967 !! - mrkiss_config::isi_num_pts
968 !! - mrkiss_config::isi_all_norm
969 !! - mrkiss_config::isi_all_y_err
970 !! - mrkiss_config::isi_one_spp_td
971 !! - mrkiss_config::isi_sdf_step
972 !! - mrkiss_config::isi_bic_max
973 !! - mrkiss_config::isi_bic_bnd
974 !! @param solution Array for solution.
975 !! Each COLUMN is a solution:
976 !! - First element is the @f$t@f$ variable
977 !! - `size(y, 1)` elements starting with 2 have @f$\mathbf{y}@f$ values
978 !! - The next `size(y, 1)` elements have @f$\mathbf{y}'@f$ values
979 !! @param deq The equation subroutine returning values for @f$\mathbf{y}'@f$ -- i.e. @f$\mathbf{f}(t, \mathbf{y})@f$
980 !! @param t Initial condition for @f$t@f$ -- i.e. @f$t_0@f$.
981 !! @param y Initial condition for @f$\mathbf{y}@f$ -- i.e. @f$\mathbf{y_0}@f$.
982 !! @param param Data payload passed to @p deq
983 !! @param a The butcher tableau @f$\mathbf{a}@f$ matrix
984 !! @param b The butcher tableau @f$\mathbf{b}@f$ vector
985 !! @param c The butcher tableau @f$\mathbf{c}@f$ vector
986 !! @param p Orders for the methods in the butcher tableau
987 !! @param t_max_o Stop if @f$t@f$ becomes greater than @p t_max_o. Different from @p t_end_o! Default: NONE
988 !! @param t_end_o Try to stop integration with @f$t@f$ equal to @p t_end. Default: NONE
989 !! @param t_delta_ini_o Initial @f$\Delta{t}@f$.
990 !! - Default when @p t_delta_max provided: `(t_delta_max_o+t_delta_min_o)/2`
991 !! - Default otherwise: mrkiss_config::t_delta_ai
992 !! @param t_delta_min_o Minimum allowed @f$\Delta{t}@f$. Default: mrkiss_config::t_delta_min_ai
993 !! @param t_delta_max_o Maximum allowed @f$\Delta{t}@f$. Default: NONE
994 !! @param t_delta_fac_min_o Minimum @f$\Delta{t}@f$ adaption factor for all but the first step. Default: mrkiss_config::t_delta_fac_min_ai
995 !! @param t_delta_fac_max_o Maximum @f$\Delta{t}@f$ adaption factor. Default: 2.0_rk
996 !! @param t_delta_fac_fdg_o Extra @f$\Delta{t}@f$ adaption factor when shrinking interval. Default: mrkiss_config::t_delta_fac_fdg_ai
997 !! @param error_tol_abs_o Absolute error tolerance. Default: mrkiss_config::error_tol_abs_ai
998 !! @param error_tol_rel_o Relative error tolerance. Default: mrkiss_config::error_tol_rel_ai
999 !! @param max_pts_o Maximum number of points to put in @p solution.
1000 !! @param max_bisect_o Maximum number of bisection iterations to perform for each step. Default: mrkiss_config::max_bisect_ai
1001 !! @param no_bisect_error_o If `.TRUE.`, then not exit on bisection errors
1002 !! @param sdf_o SDF function. Used to set new t_delta for a step. This subroutine may trigger one or
1003 !! more of the following actions:
1004 !! | Return State | Action
1005 !! |-----------------|-------
1006 !! | `status>0` | Immediately return doing nothing else and propagating status to caller.
1007 !! | `new_t_delta>0` | Recompute step using @p new_t_delta for @f$\Delta{t}@f$
1008 !! | `sdf_flags>0` | Redo step with a @f$\Delta{t}@f$ derived from the @p sdf_o via bisection.
1009 !! | `end_run>0` | Routine returns after this step is complete.
1010 !! @param sdf_tol_o How close we have to get to accept an sdf solution. Default: mrkiss_config::sdf_tol_ai
1011 !! @param stepp_o Step processing subroutine. Called after each step.
1012 !!
1013 subroutine adaptive_steps(status, istats, solution, deq, t, y, param, a, b, c, p, t_max_o, t_end_o, t_delta_ini_o, &
1014 & t_delta_min_o, t_delta_max_o, t_delta_fac_min_o, t_delta_fac_max_o, t_delta_fac_fdg_o, error_tol_abs_o, &
1015 & error_tol_rel_o, max_pts_o, max_bisect_o, no_bisect_error_o, sdf_o, sdf_tol_o, stepp_o)
1016 use mrkiss_config
1017 implicit none
1018 ! Arguments
1019 integer, intent(out) :: status, istats(istats_size)
1020 real(kind=rk), intent(out) :: solution(:,:)
1021 procedure(deq_iface) :: deq
1022 real(kind=rk), intent(in) :: t
1023 real(kind=rk), intent(in) :: y(:), param(:), a(:,:), b(:,:), c(:)
1024 integer, intent(in) :: p(:)
1025 real(kind=rk), optional, intent(in) :: t_max_o, t_end_o, t_delta_ini_o, t_delta_min_o, t_delta_max_o
1026 real(kind=rk), optional, intent(in) :: t_delta_fac_min_o, t_delta_fac_max_o, t_delta_fac_fdg_o
1027 real(kind=rk), optional, intent(in) :: error_tol_abs_o(:), error_tol_rel_o(:), sdf_tol_o
1028 integer, optional, intent(in) :: max_pts_o, max_bisect_o
1029 logical, optional, intent(in) :: no_bisect_error_o
1030 procedure(sdf_iface), optional :: sdf_o
1031 procedure(stepp_iface), optional :: stepp_o
1032 ! Variables
1033 integer :: max_pts, cur_pnt_idx, adj_cnt, y_dim, comb_err_nzc
1034 integer :: max_bisect, sp_end_run, sp_sdf_flags, bs_itr
1035 logical :: no_bisect_error, t_delta_end_p, comb_err_msk(size(y, 1))
1036 real(kind=rk) :: t_delta_fac, y_cv(size(y, 1)), y1_delta(size(y, 1)), dy(size(y, 1))
1037 real(kind=rk) :: y2_delta(size(y, 1)), t_delta_ini, t_delta_min, y_deltas(size(y, 1), size(b, 2))
1038 real(kind=rk) :: y_delta_delta(size(y, 1)), t_delta_fac_max, t_delta_fac_min
1039 real(kind=rk) :: t_delta_fac_fdg, t_delta_nxt, sdf_tol, error_tol_abs(size(y, 1))
1040 real(kind=rk) :: error_tol_rel(size(y, 1)), comb_err_tol(size(y, 1)), t_cv, t_delta
1041 real(kind=rk) :: sp_new_t_delta, bs_tmp1_t_delta, bs_tmp2_t_delta, bs_tmp_t_delta
1042 real(kind=rk) :: bs_tmp_y_delta(size(y, 1)), bs_tmp1_dist, bs_tmp2_dist, bs_tmp_dist
1043 ! Process arguments
1044 error_tol_abs = error_tol_abs_ai
1045 if (present(error_tol_abs_o)) then
1046 if (size(error_tol_abs_o, 1) < size(y, 1)) then
1047 error_tol_abs = error_tol_abs_o(1)
1048 else
1049 error_tol_abs = error_tol_abs_o
1050 end if
1051 end if
1052 error_tol_rel = error_tol_rel_ai
1053 if (present(error_tol_rel_o)) then
1054 if (size(error_tol_rel_o, 1) < size(y, 1)) then
1055 error_tol_rel = error_tol_rel_o(1)
1056 else
1057 error_tol_rel = error_tol_rel_o
1058 end if
1059 end if
1060 sdf_tol = sdf_tol_ai
1061 if (present(sdf_tol_o)) sdf_tol = sdf_tol_o
1062 t_delta_fac_fdg = t_delta_fac_fdg_ai
1063 if (present(t_delta_fac_fdg_o)) t_delta_fac_fdg = t_delta_fac_fdg_o
1064 max_bisect = max_bisect_ai
1065 if (present(max_bisect_o)) max_bisect = max_bisect_o
1066 t_delta_fac_max = t_delta_fac_max_ai
1067 if (present(t_delta_fac_max_o)) t_delta_fac_max = t_delta_fac_max_o
1068 t_delta_fac_min = t_delta_fac_min_ai
1069 if (present(t_delta_fac_min_o)) t_delta_fac_min = t_delta_fac_min_o
1070 t_delta_min = t_delta_min_ai
1071 if (present(t_delta_min_o)) t_delta_min = t_delta_min_o
1072 if (present(t_delta_ini_o)) then
1073 t_delta_ini = t_delta_ini_o
1074 else
1075 t_delta_ini = t_delta_ai
1076 if (present(t_delta_max_o)) t_delta_ini = (t_delta_max_o - t_delta_min) / 2
1077 end if
1078 max_pts = size(solution, 2)
1079 if (present(max_pts_o)) max_pts = min(max_pts, max_pts_o);
1080 no_bisect_error = .false.
1081 if (present(no_bisect_error_o)) no_bisect_error = no_bisect_error_o
1082 ! Compute solution
1083 y_dim = size(y, 1)
1084 istats = 0
1085 t_delta = t_delta_ini
1086 t_cv = t
1087 y_cv = y
1088 cur_pnt_idx = 1
1089 solution(1, cur_pnt_idx) = t_cv
1090 solution(2:(2+y_dim-1), cur_pnt_idx) = y_cv
1091 do
1092 ! Finish if we are out of space for points
1093 if (cur_pnt_idx >= max_pts) then
1094 exit
1095 end if
1096 cur_pnt_idx = cur_pnt_idx + 1
1097 ! If close to the end, adjust t_delta to hit t_end_o
1098 t_delta_end_p = .false.
1099 if (present(t_end_o)) then
1100 if (t_cv+t_delta*1.10_rk >= t_end_o) then
1101 t_delta = t_end_o - t_cv
1102 t_delta_end_p = .true.
1103 end if
1104 end if
1105 ! Do step and adaptive step
1106 do adj_cnt=1,2
1107 call step_all(status, y_deltas, dy, deq, t_cv, y_cv, param, a, b(:,1:2), c, t_delta)
1108 y1_delta = y_deltas(:,1)
1109 y2_delta = y_deltas(:,2)
1110 if (adj_cnt > 1) then
1111 istats(isi_all_y_err) = istats(isi_all_y_err) + 1
1112 else
1113 istats(isi_all_norm) = istats(isi_all_norm) + 1
1114 end if
1115 if (status > 0) return
1116 ! Compute new t_delta_nxt based on error estimate.
1117 y_delta_delta = abs(y1_delta-y2_delta)
1118 comb_err_tol = (error_tol_abs + max(abs(solution(2:(y_dim+1), cur_pnt_idx-1)), abs(y_cv + y1_delta)) * error_tol_rel)
1119 comb_err_msk = abs(comb_err_tol) > zero_epsilon
1120 comb_err_nzc = count(comb_err_msk)
1121 if (comb_err_nzc == 0) then
1122 t_delta_fac = t_delta_fac_max
1123 else
1124 t_delta_fac = (1/sqrt(sum((y_delta_delta / comb_err_tol) ** 2, comb_err_msk) / (comb_err_nzc))) ** (1 + 1/min(p(1), p(2)))
1125 if (t_delta_fac < 1.0_rk) t_delta_fac = t_delta_fac * t_delta_fac_fdg
1126 if (cur_pnt_idx > 2) t_delta_fac = max(t_delta_fac_min, t_delta_fac) ! For all but first step, limit fac
1127 t_delta_fac = min(t_delta_fac_max, t_delta_fac)
1128 end if
1129 t_delta_nxt = t_delta * t_delta_fac
1130 if (present(t_delta_min_o)) then
1131 t_delta_nxt = max(t_delta_min_o, t_delta_nxt)
1132 end if
1133 if (present(t_delta_max_o)) then
1134 t_delta_nxt = min(t_delta_max_o, t_delta_nxt)
1135 end if
1136 ! Break out of loop if all is good.
1137 if (all(y_delta_delta <= comb_err_tol)) exit
1138 ! Redoing step
1139 t_delta = t_delta_nxt
1140 t_delta_end_p = .false.
1141 ! When redoing a step close to t_end_o, we don't want t_delta to get us more than half way to t_end_o
1142 if (present(t_end_o)) then
1143 if (t_cv + t_delta * 2 >= t_end_o) then
1144 t_delta = (t_end_o - t_cv) / 2.0_rk
1145 end if
1146 end if
1147 end do
1148 if (present(stepp_o)) then
1149 call stepp_o(status, sp_end_run, sp_sdf_flags, sp_new_t_delta, cur_pnt_idx, solution, t_delta, y1_delta)
1150 if (sp_new_t_delta > 0) then
1151 if (present(t_delta_min_o)) then
1152 sp_new_t_delta = max(t_delta_min_o, sp_new_t_delta)
1153 end if
1154 if (present(t_delta_max_o)) then
1155 sp_new_t_delta = min(t_delta_max_o, sp_new_t_delta)
1156 end if
1157 t_delta = sp_new_t_delta ! Leave t_delta_nxt unchanged...
1158 call step_one(status, y1_delta, dy, deq, t_cv, y_cv, param, a, b(:,1:1), c, t_delta)
1159 istats(isi_one_spp_td) = istats(isi_one_spp_td) + 1
1160 end if
1161 if (sp_sdf_flags > 0) then
1162 bs_tmp1_t_delta = t_delta_min
1163 call step_one(status, bs_tmp_y_delta, dy, deq, t_cv, y_cv, param, a, b(:,1:1), c, bs_tmp1_t_delta)
1164 istats(isi_sdf_step) = istats(isi_sdf_step) + 1
1165 if (status > 0) return
1166 call sdf_o(status, bs_tmp1_dist, sp_sdf_flags, t_cv+bs_tmp1_t_delta, y_cv+bs_tmp_y_delta)
1167 if (status > 0) return
1168 bs_tmp2_t_delta = t_delta
1169 call step_one(status, bs_tmp_y_delta, dy, deq, t_cv, y_cv, param, a, b(:,1:1), c, bs_tmp2_t_delta)
1170 istats(isi_sdf_step) = istats(isi_sdf_step) + 1
1171 if (status > 0) return
1172 call sdf_o(status, bs_tmp2_dist, sp_sdf_flags, t_cv+bs_tmp2_t_delta, y_cv+bs_tmp_y_delta)
1173 if (status > 0) return
1174 if (bs_tmp2_dist < bs_tmp1_dist) then ! Swap if req.
1175 bs_tmp_t_delta = bs_tmp1_t_delta
1176 bs_tmp1_t_delta = bs_tmp2_t_delta
1177 bs_tmp2_t_delta = bs_tmp_t_delta
1178 bs_tmp_dist = bs_tmp1_dist
1179 bs_tmp1_dist = bs_tmp2_dist
1180 bs_tmp2_dist = bs_tmp_dist
1181 end if
1182 if (bs_tmp2_dist >= 0) then ! We don't exit if we can't bisect, we just don't do it.
1183 bs_itr = 1
1184 do
1185 bs_tmp_t_delta = (bs_tmp1_t_delta + bs_tmp2_t_delta) / 2.0_rk
1186 call step_one(status, bs_tmp_y_delta, dy, deq, t_cv, y_cv, param, a, b(:,1:1), c, bs_tmp_t_delta)
1187 istats(isi_sdf_step) = istats(isi_sdf_step) + 1
1188 if (status > 0) return
1189 call sdf_o(status, bs_tmp_dist, sp_sdf_flags, t_cv+bs_tmp_t_delta, y_cv+bs_tmp_y_delta)
1190 if (status > 0) return
1191 if (abs(bs_tmp_dist) <= sdf_tol) exit
1192 if (bs_tmp_dist < 0.0_rk) then
1193 bs_tmp1_t_delta = bs_tmp_t_delta
1194 bs_tmp1_dist = bs_tmp_dist
1195 else
1196 bs_tmp2_t_delta = bs_tmp_t_delta
1197 bs_tmp2_dist = bs_tmp_dist
1198 end if
1199 bs_itr = bs_itr + 1;
1200 if (bs_itr > max_bisect) then
1201 istats(isi_bic_max) = istats(isi_bic_max) + 1
1202 if (no_bisect_error) then
1203 exit
1204 else
1205 status = 1025
1206 return
1207 end if
1208 end if
1209 end do
1210 t_delta = bs_tmp_t_delta
1211 y1_delta = bs_tmp_y_delta
1212 else
1213 istats(isi_bic_bnd) = istats(isi_bic_bnd) + 1
1214 end if
1215 end if
1216 end if
1217 y_cv = y_cv + y1_delta
1218 t_cv = t_cv + t_delta
1219 t_delta = t_delta_nxt
1220 solution(1, cur_pnt_idx) = t_cv
1221 solution((2+y_dim):(2+2*y_dim-1), cur_pnt_idx-1) = dy
1222 solution(2:(2+y_dim-1), cur_pnt_idx) = y_cv
1223 istats(isi_num_pts) = istats(isi_num_pts) + 1
1224 status = 0;
1225 if (present(stepp_o)) then
1226 if (sp_end_run > 0) exit ! If we get here, then status was set by stepp_o
1227 end if
1228 if (present(t_max_o)) then
1229 if (t_cv > t_max_o) exit
1230 end if
1231 if (t_delta_end_p) exit
1232 end do
1233 ! Compute derivative for final solution point
1234 call deq(status, dy, t_cv, y_cv, param) ! This sets return status
1235 if (status > 0) return
1236 solution((2+y_dim):(2+2*y_dim-1), cur_pnt_idx) = dy
1237 istats(isi_num_pts) = istats(isi_num_pts) + 1
1238 end subroutine adaptive_steps
1239
1240 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1241 !> @name Multistep Meta Solvers
1242
1243 !--------------------------------------------------------------------------------------------------------------------------------
1244 !> Take a solution with @f$t@f$ values pre-populated, and take @p steps_per_pnt Runge-Kutta steps between each @f$t@f$ value.
1245 !!
1246 !! @param status Exit status
1247 !! | Value | Description
1248 !! |-----------|------------
1249 !! | -inf-0 | Everything worked
1250 !! | 0-255 | Evaluation of @p deq failed
1251 !! | 1348-1364 | Unknown error in this routine
1252 !! | 1248-1263 | Error from mrkiss_solvers_wt::step_one()
1253 !! @param istats Integer statistics for run
1254 !! - See: mrkiss_config::istats_strs for a description of elements.
1255 !! - Elements this routine updates (via fixed_t_steps() calls):
1256 !! - mrkiss_config::isi_num_pts
1257 !! - mrkiss_config::isi_one_reg
1258 !! @param solution Array for solution.
1259 !! - This array *must* have a populated @f$t@f$ sequence in new_`solution(1,:)`.
1260 !! - Each COLUMN is a solution:
1261 !! - First element is the @f$t@f$ variable if
1262 !! - `size(y, 1)` elements starting with 2 have @f$\mathbf{y}@f$ values
1263 !! - The next `size(y, 1)` elements have @f$\mathbf{y}'@f$ values
1264 !! @param deq The equation subroutine returning values for @f$\mathbf{y}'@f$ -- i.e. @f$\mathbf{f}(t, \mathbf{y})@f$
1265 !! @param y Initial condition for @f$\mathbf{y}@f$ -- i.e. @f$\mathbf{y_0}@f$.
1266 !! @param param Data payload passed to @p deq
1267 !! @param a The butcher tableau @f$\mathbf{a}@f$ matrix
1268 !! @param b The butcher tableau @f$\mathbf{{b}}@f$ vector
1269 !! @param c The butcher tableau @f$\mathbf{c}@f$ vector
1270 !! @param steps_per_pnt Number of Runge-Kutta steps to reach each solution point
1271 !! @param p_o The order for the Runge-Kutta method in the butcher tableau to enable Richardson extrapolation
1272 !!
1273 subroutine fixed_t_steps_between(status, istats, solution, deq, y, param, a, b, c, steps_per_pnt, p_o)
1274 use mrkiss_config, only: rk, istats_size
1275 implicit none
1276 ! Arguments
1277 integer, intent(out) :: status, istats(istats_size)
1278 real(kind=rk), intent(out) :: solution(:,:)
1279 procedure(deq_iface) :: deq
1280 real(kind=rk), intent(in) :: y(:), param(:), a(:,:), b(:,:), c(:)
1281 integer, intent(in) :: steps_per_pnt
1282 integer, optional, intent(in) :: p_o(:)
1283 ! Vars
1284 integer :: cur_pnt_idx, y_dim, jstats(istats_size), p(size(b, 2))
1285 real(kind=rk) :: dy(size(y, 1))
1286 ! Process arguments
1287 p = 0
1288 if (present(p_o)) p = p_o
1289 ! Compute Solution
1290 y_dim = size(y, 1)
1291 istats = 0
1292 solution(2:(2+y_dim-1), 1) = y
1293 call deq(status, dy, & ! wt2nt:IGNORE
1294 & solution(1, 1), & ! wt2nt:DELETE
1295 & y, param)
1296 if (status > 0) return
1297 solution((2+y_dim):(2+2*y_dim-1), 1) = dy
1298 do cur_pnt_idx=2,size(solution, 2)
1299 call fixed_t_steps(status, jstats, solution(:, cur_pnt_idx:cur_pnt_idx), deq, & ! wt2nt:IGNORE
1300 & solution(1, cur_pnt_idx-1), & ! wt2nt:DELETE
1301 & solution(2:(2+y_dim-1), cur_pnt_idx-1), param, a, b, c, p, &
1302 & t_end_o=solution(1,cur_pnt_idx), max_pts_o=steps_per_pnt+1)
1303 istats = istats + jstats
1304 if (status > 0) return
1305 end do
1306 end subroutine fixed_t_steps_between
1307
1308 !--------------------------------------------------------------------------------------------------------------------------------
1309 !> Create an interpolated solution from a source solution.
1310 !!
1311 !! Take a series of @f$t@f$ values and a source solution, and produces a new solution with points at the given @f$t@f$ values.
1312 !!
1313 !! In order to produce a new solution point at a given @f$t@f$ value, we first find two solution points in the source solution
1314 !! that bracket @f$t@f$. That is to say, we find two source solutions with times @f$t_0@f$ and @f$t_1@f$ such that
1315 !! @f[ t_0 \le t \le t_1 @f]
1316 !!
1317 !! If Hermite interpolation is being used (`linear_interp_o==.FALSE.`), we then compute the new @f$\mathbf{\tilde{y}}@f$ value like so:
1318 !! @f[ \mathbf{\tilde{y}}(t) = \sum_{i=0}^3 \mathbf{c_i} \left(\frac{t - t_0}{\Delta{t}}\right)^i @f]
1319 !! Where @f$\Delta{t}=t_1-t_0@f$ and the @f$\mathbf{c_i}@f$ defined as:
1320 !! @f[ \begin{array}{lcl}
1321 !! \mathbf{c_3} & = & 2\cdot\mathbf{y}(t_0) + \Delta{t}\cdot\mathbf{y}'(t_0) - 2\cdot\mathbf{y}(t_1) + \Delta{t}\cdot\mathbf{y}'(t_1) \\
1322 !! \mathbf{c_2} & = & -2\cdot\Delta{t}\mathbf{y}'(t_0) - 3\cdot\mathbf{y}(t_0) + 3\cdot\mathbf{y}(t_1) - \Delta{t}\cdot\mathbf{y}'(t_1) \\
1323 !! \mathbf{c_1} & = & \Delta{t}\cdot\mathbf{y}'(t_0) \\
1324 !! \mathbf{c_0} & = & \mathbf{y}(t_0)
1325 !! \end{array} @f]
1326 !! Note that
1327 !! @f[ \begin{array}{lcl}
1328 !! \mathbf{\tilde{y}}(t_0) & = & \mathbf{y}(t_0) \\
1329 !! \mathbf{\tilde{y}}(t_1) & = & \mathbf{y}(t_1) \\
1330 !! \mathbf{\tilde{y}}'(t_0) & = & \mathbf{y}'(t_0) \\
1331 !! \mathbf{\tilde{y}}'(t_1) & = & \mathbf{y}'(t_1)
1332 !! \end{array} @f]
1333 !! And thus the collection of Hermite approximation polynomials produce a is @f$ \mathcal{C}^\infty@f$ approximation over the
1334 !! entire solution interval. Also note that the Hermite solution provides a solution of @f$\mathcal{O}(3)@f$.
1335 !!
1336 !! If linear interpolation is being used (`linear_interp_o==.TRUE.`), we then we compute the new @f$\mathbf{\tilde{y}}@f$ value
1337 !! like so:
1338 !! @f[ \mathbf{\tilde{y}}(t) = \frac{\mathbf{y}(t_0) (t_1 - t) + \mathbf{y}(t_1) (t - t_0)}{\Delta{t}} @f]
1339 !! Where @f$\Delta{t}=t_1-t_0@f$.
1340 !!
1341 !! While the collection of linear approximation functions provide the same equality at source solution @f$t@f$ values that we
1342 !! had with Hermite interpolation, that equality doesn't extend to the values of the derivatives. So while this collection of
1343 !! linear interpolation functions provide a continuous approximation over the integration interval, this approximation is not
1344 !! generally smooth.
1345 !!
1346 !! Regardless of the interpolation method, the new values for @f$\mathbf{\tilde{y}}'(t)@f$ are directly computed as
1347 !! @f[ \mathbf{f}(t, \mathbf{\tilde{y}}) @f]
1348 !! using the @p deq and @p param arguments.
1349 !!
1350 !! @param status Exit status
1351 !! | Value | Description
1352 !! |-----------|------------
1353 !! | -inf-0 | Everything worked
1354 !! | 0-255 | Evaluation of @p deq failed
1355 !! | 1331 | Solution @f$t@f$ value out of bounds
1356 !! | 1332:1347 | Unknown error in this routine
1357 !! @param istats Integer statistics for run
1358 !! - See: mrkiss_config::istats_strs for a description of elements.
1359 !! - Elements this routine updates:
1360 !! - mrkiss_config::isi_num_pts
1361 !! @param solution Array for new solution.
1362 !! - This array *must* have a populated @f$t@f$ sequence in solution(1,:)
1363 !! - New @f$\mathbf{y}@f$ values are interpolated. New @f$\mathbf{y}'@f$ values are computed from @p deq.
1364 !! - Each COLUMN is a solution:
1365 !! - First element is the @f$t@f$ variable if
1366 !! - `size(y, 1)` elements starting with 2 have @f$\mathbf{y}@f$ values
1367 !! - The next `size(y, 1)` elements have @f$\mathbf{y}'@f$ values
1368 !! @param src_solution Array for source solution.
1369 !! - This array *must* have at least two solution points.
1370 !! @param deq The equation subroutine returning values for @f$\mathbf{y}'@f$ -- i.e. @f$\mathbf{f}(t, \mathbf{y})@f$
1371 !! @param param Data payload passed to @p deq
1372 !! @param num_src_pts_o The number of solutions in src_solution. Default inferred from size of @p src_solution.
1373 !! @param linear_interp_o If `.TRUE.` do linear interpolation, and Hermite otherwise. Default: `.FALSE.`
1374 !!
1375 subroutine interpolate_solution(status, istats, solution, src_solution, deq, param, num_src_pts_o, linear_interp_o)
1376 use :: mrkiss_config, only: rk, istats_size, isi_num_pts
1377 implicit none
1378 ! Arguments
1379 integer, intent(out) :: status, istats(istats_size)
1380 real(kind=rk), intent(inout) :: solution(:,:)
1381 real(kind=rk), intent(in) :: src_solution(:,:)
1382 integer, optional, intent(in) :: num_src_pts_o
1383 procedure(deq_iface) :: deq
1384 real(kind=rk), intent(in) :: param(:)
1385 logical, optional, intent(in) :: linear_interp_o
1386 ! Variables
1387 integer :: new_sol_idx, old_sol_idx, max_idx, num_src_pts, y_dim
1388 logical :: linear_interp
1389 real(kind=rk) :: t, t0, t1, tu, td
1390 real(kind=rk), allocatable :: y0(:), y1(:), dy0(:), dy1(:), yat(:)
1391 ! Process Arguments
1392 linear_interp = .false.
1393 if (present(linear_interp_o)) linear_interp = linear_interp_o
1394 num_src_pts = size(src_solution, 2)
1395 if (present(num_src_pts_o)) num_src_pts = min(num_src_pts_o, size(src_solution, 2))
1396 ! Compute value
1397 istats = 0
1398 y_dim = (size(src_solution, 1) - 1) /2
1399 max_idx = size(solution, 2)
1400 old_sol_idx = 2
1401 do new_sol_idx=1, max_idx
1402 t = solution(1, new_sol_idx)
1403 do while (t > src_solution(1, old_sol_idx))
1404 old_sol_idx = old_sol_idx + 1
1405 if (old_sol_idx > num_src_pts) then
1406 status = 1331
1407 return
1408 end if
1409 end do
1410 ! If we get here, we are in the interval we want
1411 t0 = src_solution(1, old_sol_idx-1)
1412 t1 = src_solution(1, old_sol_idx)
1413 y0 = src_solution(2:(2+y_dim-1), old_sol_idx-1)
1414 y1 = src_solution(2:(2+y_dim-1), old_sol_idx)
1415 td = (t1 - t0)
1416 if (linear_interp) then
1417 yat = (y0 * (t1 - t) + y1 * (t - t0)) / td
1418 else
1419 dy0 = td * src_solution((2+y_dim):(2+2*y_dim-1), old_sol_idx-1)
1420 dy1 = td * src_solution((2+y_dim):(2+2*y_dim-1), old_sol_idx)
1421 tu = (t - t0) / td
1422 yat = tu * (tu * (tu * (2 * y0 + dy0 - 2 * y1 + dy1) - 2 * dy0 - 3 * y0 + 3 * y1 - dy1) + dy0) + y0
1423 end if
1424 solution(2:(2+y_dim-1), new_sol_idx) = yat
1425 call deq(status, solution((2+y_dim):(2+2*y_dim-1), new_sol_idx), & ! wt2nt:IGNORE
1426 & t, & ! wt2nt:DELETE
1427 & yat, param)
1428 if (status > 0) return
1429 istats(isi_num_pts) = istats(isi_num_pts) + 1
1430 end do
1431 status = 0;
1432 end subroutine interpolate_solution
1433
1434end module mrkiss_solvers_wt
mrkiss_solvers_wt::deq_iface
Type for ODE dydt functions.
Definition mrkiss_solvers_wt.f90:62
mrkiss_solvers_wt::sdf_iface
Type SDF function on a solution point.
Definition mrkiss_solvers_wt.f90:117
mrkiss_solvers_wt::stepp_iface
Type step processing subroutine.
Definition mrkiss_solvers_wt.f90:90
mrkiss_config
Configuration for MRKISS == MR RK KISS == Mitch Richling's Runge-Kutta Keep It Simple Stupid.
Definition mrkiss_config.f90:38
mrkiss_config::isi_all_y_err
integer, parameter, public isi_all_y_err
istats index one_stab_y_delta_err
Definition mrkiss_config.f90:65
mrkiss_config::t_delta_fac_fdg_ai
real(kind=rk), parameter, public t_delta_fac_fdg_ai
t_delta_fac_fdg ai default
Definition mrkiss_config.f90:113
mrkiss_config::isi_num_pts
integer, parameter, public isi_num_pts
istats index num_pts
Definition mrkiss_config.f90:61
mrkiss_config::zero_epsilon
real(kind=rk), parameter, public zero_epsilon
Used to test for zero.
Definition mrkiss_config.f90:101
mrkiss_config::isi_bic_max
integer, parameter, public isi_bic_max
istats index bisect_fail_max
Definition mrkiss_config.f90:69
mrkiss_config::t_delta_ai
real(kind=rk), parameter, public t_delta_ai
t_delta ai default
Definition mrkiss_config.f90:107
mrkiss_config::isi_one_y_len
integer, parameter, public isi_one_y_len
istats index one_step_y_delta_len
Definition mrkiss_config.f90:66
mrkiss_config::isi_one_spp_td
integer, parameter, public isi_one_spp_td
istats index one_step_stepp_t_delta
Definition mrkiss_config.f90:67
mrkiss_config::rk
integer, parameter, public rk
Real kind used across the library.
Definition mrkiss_config.f90:46
mrkiss_config::t_delta_min_ai
real(kind=rk), parameter, public t_delta_min_ai
t_delta_min ai default
Definition mrkiss_config.f90:115
mrkiss_config::isi_one_reg
integer, parameter, public isi_one_reg
istats index one_step_norm
Definition mrkiss_config.f90:62
mrkiss_config::sdf_tol_ai
real(kind=rk), parameter, public sdf_tol_ai
sdf_tol ai default
Definition mrkiss_config.f90:114
mrkiss_config::isi_sdf_step
integer, parameter, public isi_sdf_step
istats index one_step_sdf_bisection
Definition mrkiss_config.f90:68
mrkiss_config::error_tol_abs_ai
real(kind=rk), parameter, public error_tol_abs_ai
error_tol_abs ai default
Definition mrkiss_config.f90:108
mrkiss_config::istats_size
integer, parameter, public istats_size
Number of elements in istats(:)
Definition mrkiss_config.f90:58
mrkiss_config::isi_all_norm
integer, parameter, public isi_all_norm
istats index one_etab_norm
Definition mrkiss_config.f90:63
mrkiss_config::isi_bic_bnd
integer, parameter, public isi_bic_bnd
istats index bisect_fail_containment
Definition mrkiss_config.f90:70
mrkiss_config::error_tol_rel_ai
real(kind=rk), parameter, public error_tol_rel_ai
error_tol_rel ai default
Definition mrkiss_config.f90:109
mrkiss_config::t_delta_fac_min_ai
real(kind=rk), parameter, public t_delta_fac_min_ai
t_delta_fac_min ai default
Definition mrkiss_config.f90:111
mrkiss_config::t_delta_fac_max_ai
real(kind=rk), parameter, public t_delta_fac_max_ai
t_delta_fac_max ai default
Definition mrkiss_config.f90:112
mrkiss_config::max_bisect_ai
integer, parameter, public max_bisect_ai
max_bisect ai default
Definition mrkiss_config.f90:110
mrkiss_solvers_wt
Solvers.
Definition mrkiss_solvers_wt.f90:40
mrkiss_solvers_wt::fixed_t_steps_between
subroutine, public fixed_t_steps_between(status, istats, solution, deq, y, param, a, b, c, steps_per_pnt, p_o)
Take a solution with values pre-populated, and take steps_per_pnt Runge-Kutta steps between each va...
Definition mrkiss_solvers_wt.f90:1274
mrkiss_solvers_wt::step_one
subroutine, public step_one(status, y_delta, dy, deq, t, y, param, a, b, c, t_delta)
Compute one step of a Runge-Kutta method expressed as a Butcher Tableau.
Definition mrkiss_solvers_wt.f90:208
mrkiss_solvers_wt::step_all
subroutine, public step_all(status, y_deltas, dy, deq, t, y, param, a, b, c, t_delta)
Compute one step of a embedded Runge-Kutta method expressed as a Butcher Tableau.
Definition mrkiss_solvers_wt.f90:153
mrkiss_solvers_wt::step_rk4
subroutine, public step_rk4(status, y_delta, dy, deq, t, y, param, t_delta)
Compute one step of RK4 – the classical 4th order Runge-Kutta method (mrkiss_erk_kutta_4).
Definition mrkiss_solvers_wt.f90:326
mrkiss_solvers_wt::interpolate_solution
subroutine, public interpolate_solution(status, istats, solution, src_solution, deq, param, num_src_pts_o, linear_interp_o)
Create an interpolated solution from a source solution.
Definition mrkiss_solvers_wt.f90:1376
mrkiss_solvers_wt::adaptive_steps
subroutine, public adaptive_steps(status, istats, solution, deq, t, y, param, a, b, c, p, t_max_o, t_end_o, t_delta_ini_o, t_delta_min_o, t_delta_max_o, t_delta_fac_min_o, t_delta_fac_max_o, t_delta_fac_fdg_o, error_tol_abs_o, error_tol_rel_o, max_pts_o, max_bisect_o, no_bisect_error_o, sdf_o, sdf_tol_o, stepp_o)
Take multiple adaptive steps with an embedded Runge-Kutta method using relatively traditional step si...
Definition mrkiss_solvers_wt.f90:1016
mrkiss_solvers_wt::fixed_y_steps
subroutine, public fixed_y_steps(status, istats, solution, deq, t, y, param, a, b, c, y_delta_len_targ, t_delta_max, t_delta_min_o, y_delta_len_tol_o, max_bisect_o, no_bisect_error_o, y_delta_len_idxs_o, max_pts_o, y_sol_len_max_o, t_max_o)
Take multiple adaptive steps with constant length using a simple Runge-Kutta method.
Definition mrkiss_solvers_wt.f90:633
mrkiss_solvers_wt::sloppy_fixed_y_steps
subroutine, public sloppy_fixed_y_steps(status, istats, solution, deq, t, y, param, a, b, c, y_delta_len_targ, t_delta_ini, t_max_o, t_delta_min_o, t_delta_max_o, y_delta_len_idxs_o, adj_short_o, max_pts_o, y_sol_len_max_o)
Take multiple adaptive steps adjusting for the length of using a simple Runge-Kutta method.
Definition mrkiss_solvers_wt.f90:837
mrkiss_solvers_wt::step_richardson
subroutine, public step_richardson(status, y_delta, dy, deq, t, y, param, a, b, c, p, t_delta)
Compute one Richardson Extrapolation Step.
Definition mrkiss_solvers_wt.f90:275
mrkiss_solvers_wt::fixed_t_steps
subroutine, public fixed_t_steps(status, istats, solution, deq, t, y, param, a, b, c, p_o, max_pts_o, t_delta_o, t_end_o, t_max_o)
Take multiple fixed time steps with a simple Runge-Kutta method.
Definition mrkiss_solvers_wt.f90:489
mrkiss_solvers_wt::step_dp54
subroutine, public step_dp54(status, y_deltas, dy, deq, t, y, param, t_delta)
Compute one step of DP45 (mrkiss_eerk_dormand_prince_5_4)
Definition mrkiss_solvers_wt.f90:415
mrkiss_solvers_wt::step_rkf45
subroutine, public step_rkf45(status, y_deltas, dy, deq, t, y, param, t_delta)
Compute one step of RKF45 (mrkiss_eerk_fehlberg_4_5).
Definition mrkiss_solvers_wt.f90:367