integrate_to_event Subroutine

private subroutine integrate_to_event(me, t0, x0, h, tmax, tol, tf, xf, gf)

Uses

  • proc~~integrate_to_event~2~~UsesGraph proc~integrate_to_event~2 rk_module::rk_class%integrate_to_event module~brent_module brent_module proc~integrate_to_event~2->module~brent_module module~kind_module kind_module module~brent_module->module~kind_module module~numbers_module numbers_module module~brent_module->module~numbers_module iso_fortran_env iso_fortran_env module~kind_module->iso_fortran_env module~numbers_module->module~kind_module

Event-finding integration routine for the rk_class. Integrates until g(t,x)=0, or until t=tf (whichever happens first).

Note

There are some efficiency improvements that could be made here. This is a work in progress.

Type Bound

rk_class

Arguments

Type IntentOptional Attributes Name
class(rk_class), intent(inout) :: me
real(kind=wp), intent(in) :: t0

initial time
real(kind=wp), intent(in), dimension(me%n) :: x0

initial state
real(kind=wp), intent(in) :: h

abs(time step)
real(kind=wp), intent(in) :: tmax

max final time if event not located
real(kind=wp), intent(in) :: tol

function tolerance for root finding
real(kind=wp), intent(out) :: tf

actual final time reached
real(kind=wp), intent(out), dimension(me%n) :: xf

final state (at tf)
real(kind=wp), intent(out) :: gf

g value at tf

Calls

proc~~integrate_to_event~2~~CallsGraph proc~integrate_to_event~2 rk_module::rk_class%integrate_to_event proc~set_function brent_module::brent_class%set_function proc~integrate_to_event~2->proc~set_function proc~zeroin brent_module::brent_class%zeroin proc~integrate_to_event~2->proc~zeroin step step proc~integrate_to_event~2->step

Called by

proc~~integrate_to_event~2~~CalledByGraph proc~integrate_to_event~2 rk_module::rk_class%integrate_to_event proc~halo_to_rv_diffcorr halo_orbit_module::halo_to_rv_diffcorr proc~halo_to_rv_diffcorr->proc~integrate_to_event~2 proc~rk_test rk_module::rk_test proc~rk_test->proc~integrate_to_event~2

Source Code

    subroutine integrate_to_event(me,t0,x0,h,tmax,tol,tf,xf,gf)

    use brent_module

    implicit none

    class(rk_class),intent(inout)        :: me
    real(wp),intent(in)                  :: t0      !! initial time
    real(wp),dimension(me%n),intent(in)  :: x0      !! initial state
    real(wp),intent(in)                  :: h       !! abs(time step)
    real(wp),intent(in)                  :: tmax    !! max final time if event not located
    real(wp),intent(in)                  :: tol     !! function tolerance for root finding
    real(wp),intent(out)                 :: tf      !! actual final time reached
    real(wp),dimension(me%n),intent(out) :: xf      !! final state (at tf)
    real(wp),intent(out)                 :: gf      !! g value at tf

    !local variables:
    real(wp) :: t,dt,t2,ga,gb,dt_root,dum
    real(wp),dimension(me%n) :: x,g_xf
    logical :: first,last,export
    procedure(report_func),pointer :: report
    type(brent_class) :: solver
    integer :: iflag

    if (.not. associated(me%f)) error stop 'Error in integrate_to_event: f is not associated.'
    if (.not. associated(me%g)) error stop 'Error in integrate_to_event: g is not associated.'
    if (h==zero) error stop 'Error in integrate_to_event: h must not be zero.'

    !If the points are being exported:
    export = associated(me%report)

    !first point:
    if (export) call me%report(t0,x0)

    if (t0==tmax) then
        xf = x0
        tf = t0
        call me%g(t0,x0,gf)
    else

        first = .true.
        t = t0
        x = x0
        call me%g(t0,x0,ga)     !evaluate event function
        dt = sign(h,tmax-t0)    !time step (correct sign)

        do

            t2 = t + dt
            last = ((dt>=zero .and. t2>=tmax) .or. &  !adjust last time step
                    (dt<zero .and. t2<=tmax))         !
            if (last) then
                dt = tmax-t
                t2 = tmax
            end if
            call me%step(t,x,dt,xf)
            call me%g(t2,xf,gb)     !evaluate event function

            if (first .and. abs(ga)<=tol) then

                !we ignore a root at t0 after the first step
                if (abs(gb)<=tol) then !check this one since it could have landed on a root
                    gf = gb
                    tf = t2
                    exit
                else
                    if (last) then  !exiting without having found a root
                        tf = t2
                        gf = gb
                        exit
                    end if
                    if (export) call me%report(t2,xf)   !intermediate point
                    x = xf
                    t = t2
                    ga = gb
                end if

            elseif (ga*gb<=zero) then !there is a root somewhere on [t,t+dt]

                !find the root:
                call solver%set_function(solver_func)
                call solver%find_zero(zero,dt,tol,dt_root,dum,iflag,ga,gb)
                t2 = t + dt_root
                gf = solver_func(solver,dt_root)
                tf = t2
                xf = g_xf !computed in the solver function
                exit

            else  !no root yet, continue

                if (last) then  !exiting without having found a root
                    tf = t2
                    gf = gb
                    exit
                end if
                if (export) call me%report(t2,xf)   !intermediate point
                x = xf
                t = t2
                ga = gb

            end if

            if (first) first = .false.

        end do

    end if

    if (export) call me%report(t2,xf)   !last point

    contains

        function solver_func(this,delt) result(g)

        !! root solver function. The input is the dt offset from time t.

        implicit none

        class(brent_class),intent(inout) :: this
        real(wp),intent(in) :: delt  !! from [0 to dt]
        real(wp) :: g

        !take a step from t to t+delt and evaluate g function:
        call me%step(t,x,delt,g_xf)
        call me%g(t+delt,g_xf,g)

        end function solver_func

    end subroutine integrate_to_event