easydop853/easydop853_module.f90 at main · GasinAn/easydop853 · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
!*****************************************************************************************
!> author: Jiachen An
!
!  Simplified Modern Fortran Edition of the DOP853 ODE Solver.
!
!### License
!
!        Simplified Modern Fortran Edition of the DOP853 ODE Solver
!        https://github.com/GasinAn/easydop853
!
!        Copyright (c) 2020-2024, Jiachen An
!        All rights reserved.
!
!        Redistribution and use in source and binary forms, with or without modification,
!        are permitted provided that the following conditions are met:
!
!        * Redistributions of source code must retain the above copyright notice, this
!          list of conditions and the following disclaimer.
!
!        * Redistributions in binary form must reproduce the above copyright notice, this
!          list of conditions and the following disclaimer in the documentation and/or
!          other materials provided with the distribution.
!
!        * The names of its contributors may not be used to endorse or promote products
!          derived from this software without specific prior written permission.
!
!        THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
!        ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
!        WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
!        DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
!        ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
!        (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
!        LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
!        ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
!        (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
!        SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
!*****************************************************************************************

module easydop853_module

    use dop853_module
    use dop853_constants

    implicit none

    contains

    subroutine easydop853(fcn, x, xend, y)

        implicit none

        procedure(deriv_func) :: fcn
        !! subroutine computing the value of \(dy/dx=f(x,y)\)
        real(wp), intent(inout) :: x
        !! `x` value (input is initial value and output is final value)
        real(wp), intent(in) :: xend
        !! endpoint of integration (final value of `x`)
        real(wp), dimension(:), intent(inout) :: y
        !! `y` value (input is initial value and output is final value)

        integer, parameter :: nstiff = 1
        !! nstiff parameter for stiffness detection,
        !! which will occur at step 1*nstiff, 2*nstiff, 3*nstiff ... if nstiff>0
        !! and will not occur if nstiff<=0
        integer, parameter :: nmax = 2250000
        !! maximal number of allowed steps
        real(wp), dimension(1), parameter :: rtol = 1.0e-12_wp
        !! relative tolerance
        real(wp), dimension(1), parameter :: atol = 1.0e-24_wp
        !! absolute tolerance

        type(dop853_class) :: prop
        logical :: status_ok
        integer :: idid

        call prop%initialize(n=size(y), fcn=fcn, nstiff=nstiff, nmax=nmax, &
                             status_ok=status_ok)
        if (.not. status_ok) error stop 'initialization error'

        call prop%integrate(x, y, xend, rtol, atol, iout=0, idid=idid)
        if (idid < 0) error stop 'integration failure'

    end subroutine easydop853

end module easydop853_module