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