ddeabm_with_event_wrapper

private subroutine ddeabm_with_event_wrapper(me, t, y, tmax, tstop, idid, gval, integration_mode, tstep, continue)

Wrapper routine for ddeabm, with event finding. It will integrate until g(t,x)=0 or t=tmax (whichever comes first). Note that a root at the initial time is ignored (user should check for this manually)

If starting a new problem, must first call me%first_call().

Type Bound

ddeabm_with_event_class

Arguments

Type	Intent	Optional	Attributes		Name
class(ddeabm_with_event_class),	intent(inout)			::	me
real(kind=wp),	intent(inout)			::	t
real(kind=wp),	intent(inout),		dimension(me%neq)	::	y
real(kind=wp),	intent(in)			::	tmax	max time at which a solution is desired. (if root not located, it will integrate to `tmax`)
real(kind=wp),	intent(in),	optional		::	tstop	for some problems it may not be permissible to integrate past a point `tstop` because a discontinuity occurs there or the solution or its derivative is not defined beyond `tstop`. when you have declared a `tstop` point (see `info(4)`), you have told the code not to integrate past `tstop`. in this case any `tmax` beyond `tstop` is invalid input. [not used if not present]
integer,	intent(out)			::	idid	indicator reporting what the code did. you must monitor this integer variable to decide what action to take next. `idid=1000` means a root was found. See ddeabm for other values.
real(kind=wp),	intent(out)			::	gval	value of the event function `g(t,x)` at the final time `t` (not set if the inputs are invalid: `idid=-33`)
integer,	intent(in),	optional		::	integration_mode	Step mode: 1 - normal integration from `t` to `tout`, no reporting [default]. 2 - normal integration from `t` to `tout`, report each step.
real(kind=wp),	intent(in),	optional		::	tstep	Fixed time step to use for reporting and evaluation of the event function (the points are interpolated within each integrator step). If not present, then default integrator steps are used. The event function is also evaluated at the end of each integrator step.
logical,	intent(in),	optional		::	continue	to continue after a previous event location. This option can be used after a previous call to this routine has found an event, to continue integration to the next event without having to restart the integration. It will not report the initial point (which would have been reported as the last point of the previous call). The rest of the integration step in which the previous event was found is checked for events (using interpolation) before integration continues. Limitation: on that remaining interval, the event function is ignored (it is zero at the start), so a second root of it within the same integration step will not be found.

Calls

Help

Source Code

   subroutine ddeabm_with_event_wrapper(me, t, y, tmax, tstop, idid, gval, integration_mode, tstep, continue)

      implicit none

      class(ddeabm_with_event_class), intent(inout)   :: me
      real(wp), intent(inout)                         :: t
      real(wp), dimension(me%neq), intent(inout)       :: y
      real(wp), intent(in)                            :: tmax    !! max time at which a solution is desired.
                                                              !! (if root not located, it will integrate to `tmax`)
      real(wp), intent(in), optional                   :: tstop   !! for some problems it may not be permissible to integrate
                                                              !! past a point `tstop` because a discontinuity occurs there
                                                              !! or the solution or its derivative is not defined beyond
                                                              !! `tstop`.  when you have declared a `tstop` point (see `info(4)`),
                                                              !! you have told the code not to integrate past `tstop`.
                                                              !! in this case any `tmax` beyond `tstop` is invalid input.
                                                              !! [not used if not present]
      integer, intent(out)                            :: idid    !! indicator reporting what the code did.
                                                              !! you must monitor this integer variable to
                                                              !! decide what action to take next.
                                                              !! `idid=1000` means a root was found.
                                                              !! See [[ddeabm]] for other values.
      real(wp), intent(out)                           :: gval    !! value of the event function `g(t,x)` at the final time `t`
                                                              !! (not set if the inputs are invalid: `idid=-33`)
      integer, intent(in), optional :: integration_mode           !! Step mode:
                                                              !! *1* - normal integration from `t` to `tout`,
                                                              !!  no reporting [default].
                                                              !! *2* - normal integration from `t` to `tout`,
                                                              !!  report each step.
      real(wp), intent(in), optional :: tstep    !! Fixed time step to use for reporting and
                                             !! evaluation of the event function (the points are
                                             !! interpolated within each integrator step). If not
                                             !! present, then default integrator steps are used.
                                             !! The event function is also evaluated at the end of
                                             !! each integrator step.
      logical, intent(in), optional :: continue !! to continue after a previous event location.
                                            !! This option can be used after a previous call to this routine
                                            !! has found an event, to continue integration to the next event
                                            !! without having to restart the integration. It will not report
                                            !! the initial point (which would have been reported as the last
                                            !! point of the previous call). The rest of the integration step
                                            !! in which the previous event was found is checked for events
                                            !! (using interpolation) before integration continues.
                                            !! Limitation: on that remaining interval, the event function
                                            !! is ignored (it is zero at the start), so a second root of it
                                            !! within the same integration step will not be found.

      ! local variables:
      real(wp), dimension(:), allocatable :: y1  !! initial state in an interval
      real(wp), dimension(:), allocatable :: y2  !! final state in an interval
      real(wp), dimension(:), allocatable :: yc  !! interpolated state at `tc`
      real(wp) :: g1          !! value of event function at `t1`
      real(wp) :: g2          !! value of event function at `t2`
      real(wp) :: t1          !! initial time of an interval
      real(wp) :: t2          !! final time of an interval
      integer :: mode         !! local copy of integration_mode
      logical :: first        !! flag for the first interval
      logical :: fixed_step   !! if reporting at the fixed step size `tstep`
      logical :: continuing   !! local copy of optional `continue` argument
      logical :: resuming     !! if the rest of the last step must be checked first
      logical :: found        !! if a root was found (or root finding failed)
      logical :: done         !! if finished after checking the rest of the last step
      real(wp) :: dt          !! for fixed step size: actual signed time step
      real(wp) :: direction   !! direction of integration: `+1: tmax>=t`, `-1: tmax<t`
      real(wp) :: t_anchor    !! for fixed step size: start of the reporting grid
      integer :: n_grid       !! for fixed step size: number of grid points reported after `t_anchor`

      idid = 0

      ! optional inputs:
      if (present(integration_mode)) then
         mode = integration_mode
      else
         mode = 1  !default
      end if
      if (present(continue)) then
         ! if there has not yet been a successful
         ! step yet, then proceed as normal without
         ! continuing. Otherwise, enable continue mode.
         continuing = (continue .and. allocated(me%x_saved))
      else
         continuing = .false.
      end if
      resuming = continuing .and. me%resume_pending
      fixed_step = present(tstep)

      !check for invalid inputs:
      if (mode /= 1 .and. mode /= 2) then
         call report_error('ddeabm_with_event_wrapper', &
                           'integration_mode must be 1 or 2.', 0, 0)
         idid = -33
         return
      end if
      if ((mode == 2) .and. (.not. associated(me%report))) then
         call report_error('ddeabm_with_event_wrapper', &
                           'REPORT procedure must be associated for integration_mode=2.', 0, 0)
         idid = -33
         return
      end if
      if (.not. associated(me%gfunc)) then
         call report_error('ddeabm_with_event_wrapper', &
                           'the event function GFUNC has not been associated.', 0, 0)
         idid = -33
         return
      end if
      ! direction of integration (it cannot be changed when continuing):
      if (continuing) then
         direction = me%delsgn
      else
         direction = sign(1.0_wp, tmax - t)
      end if

      ! the fixed step reporting grid starts at the initial point,
      ! or at the last root found if continuing from one:
      if (resuming) then
         t_anchor = me%t_resume
      else if (continuing) then
         t_anchor = me%t_saved
      else
         t_anchor = t
      end if
      n_grid = 0

      if (resuming) then
         ! (the rest of the last step may be processed using interpolation
         ! only, so these are checked here rather than by ddes)
         if (direction*(tmax - me%t_resume) < 0.0_wp) then
            call report_error('ddeabm_with_event_wrapper', &
                              'cannot change the direction of integration when continuing.', 0, 0)
            idid = -33
            return
         end if
         if (present(tstop)) then
            if (direction*(tstop - me%t_resume) < 0.0_wp .or. direction*(tmax - tstop) > 0.0_wp) then
               call report_error('ddeabm_with_event_wrapper', &
                                 'tmax is beyond tstop. these instructions conflict.', 0, 0)
               idid = -33
               return
            end if
         end if
      end if

      if (fixed_step) then
         ! the step must advance the time everywhere on the interval
         ! (it is not lost in roundoff), and the number of points
         ! must not overflow the counter:
         dt = direction*abs(tstep)
         if (t_anchor + dt == t_anchor .or. tmax - dt == tmax .or. &
             abs(tmax - t_anchor) >= abs(dt)*real(huge(n_grid), wp)) then
            call report_error('ddeabm_with_event_wrapper', &
                              'tstep is too small (or zero) to advance the time.', 0, 0)
            idid = -33
            return
         end if
      end if

      ! work arrays:
      allocate (y1(me%neq), y2(me%neq), yc(me%neq))

      if (continuing) then
         ! if continuing, then we reset the t,y inputs
         ! to the values from the last successful step
         ! (note that an event may have been found in the
         ! last call, so the input values are not at the
         ! end of the step)
         t = me%t_saved
         y = me%x_saved
         me%info(1) = 1   ! continuation call
      end if

      !set info array:

      !info(1) is set when ddeabm_new_problem is called

      !info(2)
      if (me%scalar_tols) then
         me%info(2) = 0
      else
         me%info(2) = 1
      end if

      !info(3): intermediate-output mode: return after each integration step.
      ! the event is checked on each step, since the interpolating
      ! polynomial is only valid on the last step. If using a fixed
      ! step size, the reported points are interpolated within each step.
      me%info(3) = 1

      !info(4)
      if (present(tstop)) then
         me%info(4) = 1
         me%tstop = tstop
      else
         me%info(4) = 0
         me%tstop = 0.0_wp !not used
      end if

      call me%set_working_tols()

      if (resuming) then
         ! the last step was only checked for events up to the
         ! root that was found, so check the rest of it first:
         call check_rest_of_last_step(done)
         if (done) return
      else
         ! evaluate the event function at the initial point:
         t1 = t
         y1 = y
         call me%gfunc(t1, y1, g1)
         if (mode == 2 .and. .not. continuing) call me%report(t, y)  !initial point
      end if
      first = .true.  ! to ignore a root at the start of the first interval

      do

         ! take one step in the direction of tmax:
         call me%ddeabm(neq=me%neq, &
                        t=t, &
                        y=y, &
                        tout=tmax, &
                        info=me%info, &
                        rtol=me%rtol_tmp, &
                        atol=me%atol_tmp, &
                        idid=idid)

         !if there was a user-triggered error, return the last good point:
         if (me%error) then
            idid = -1000
            t = t1
            y = y1
            gval = g1
            me%resume_pending = .false.
            return
         end if

         ! invalid input (see ddeabm or idid codes): no integration was
         ! attempted, so the event function is not evaluated (gval is not set):
         if (idid == -33) then
            me%resume_pending = .false.
            return
         end if

         ! an interrupted task (idid = -1 to -4, see ddeabm or idid codes).
         ! (in intermediate-output mode, the integration has not advanced)
         if (idid < 0) then
            call me%gfunc(t, y, gval)
            me%resume_pending = .false.
            return
         end if

         !save the last successful step in the class:
         me%t_saved = t
         me%x_saved = y

         !evaluate event function at new point:
         t2 = t
         y2 = y
         call me%gfunc(t2, y2, g2)

         ! check for a root on [t1,t2] and report the points. on the first
         ! interval, ignore a root at the initial time (or at the root
         ! where a previous call stopped):
         call process_interval(t1, t2, y1, y2, g1, g2, first .and. abs(g1) <= me%tol, found)
         if (found) return

         ! integration is finished (see ddeabm or idid codes):
         if (idid == 2 .or. idid == 3) then
            gval = g2
            me%resume_pending = .false.
            return
         end if

         ! set up for next step:
         first = .false.
         g1 = g2
         t1 = t2
         y1 = y2

      end do

   contains

      subroutine process_interval(ta, tb, ya, yb, ga, gb, mask, found)

            !! Check for roots of the event function on `[ta,tb]`
            !! (which must be within the last integration step), and
            !! report the points on the interval (if `mode=2`).
            !!
            !! If using a fixed step size, the interval is divided at the points on
            !! the fixed step grid (using interpolation), so that the event function
            !! is evaluated at least every `tstep`, and each subinterval is checked
            !! in order.
            !!
            !! If a root is found, the outputs (`t`, `y`, `idid`, `gval`) are set.

         implicit none

         real(wp), intent(in)               :: ta    !! start of the interval
         real(wp), intent(in)               :: tb    !! end of the interval
         real(wp), dimension(:), intent(in) :: ya    !! state at `ta`
         real(wp), dimension(:), intent(in) :: yb    !! state at `tb`
         real(wp), intent(in)               :: ga    !! event function at `ta`
         real(wp), intent(in)               :: gb    !! event function at `tb`
         logical, intent(in)                :: mask  !! if the function is to be ignored on the first subinterval
         logical, intent(out)               :: found !! if a root was found
                                                     !! (or there was an error locating it)

         real(wp) :: tl                     !! start of the current subinterval
         real(wp) :: tn                     !! next point on the fixed step grid
         real(wp), dimension(me%neq) :: yl  !! state at `tl`
         real(wp), dimension(me%neq) :: yn  !! state at `tn`
         real(wp) :: gl, gn  !! event function at `tl` and `tn`
         logical :: m  !! mask for the current subinterval

         found = .false.
         if (tb == ta) return

         tl = ta
         yl = ya
         gl = ga
         m = mask

         if (fixed_step) then
            ! subintervals ending at the grid points in (ta,tb), before tmax:
            do
               tn = t_anchor + real(n_grid + 1, wp)*dt
               if (direction*(tn - tl) <= 0.0_wp) then
                  ! (should not happen) skip points not after tl
                  n_grid = n_grid + 1
                  cycle
               end if
               if (direction*(tn - tb) >= 0.0_wp .or. direction*(tn - tmax) >= 0.0_wp) exit
               call me%interpolate(tn, yn)
               call me%gfunc(tn, yn, gn)
               call locate_root(tl, tn, yl, yn, gl, gn, m, found)
               if (found) then
                  call root_located()
                  return
               end if
               if (mode == 2) call me%report(tn, yn)
               n_grid = n_grid + 1
               tl = tn
               yl = yn
               gl = gn
               m = .false.
            end do
         end if

         ! the last subinterval (or the whole interval):
         call locate_root(tl, tb, yl, yb, gl, gb, m, found)
         if (found) then
            call root_located()
            return
         end if

         ! (tn is only defined if fixed_step)
         if (fixed_step) then
            if (tb == tmax) then
               if (mode == 2) call me%report(tb, yb) ! the last point
            else if (tn == tb) then
               if (mode == 2) call me%report(tb, yb) ! a grid point at the end of the step
               n_grid = n_grid + 1
            end if
         else
            if (mode == 2) call me%report(tb, yb) ! the end of the integration step (or tmax)
         end if

      end subroutine process_interval

      subroutine root_located()

            !! A root (or an error locating it) was found by `locate_root`:
            !! report it, and set up for continuing from it.

         implicit none

         if (idid > 0) then ! a root was found
            if (mode == 2) call me%report(t, y) ! report the root
            ! a subsequent call with continue=.true. will
            ! check the rest of the step from this point:
            me%t_resume = t
            me%resume_pending = .true.
         end if

      end subroutine root_located

      subroutine locate_root(ta, tb, ya, yb, ga, gb, mask, found)

            !! Check for a root of the event function on `[ta,tb]` and,
            !! if there is one, locate it.
            !! If found, the outputs (`t`, `y`, `idid`, `gval`) are set.

         implicit none

         real(wp), intent(in)               :: ta    !! start of the interval
         real(wp), intent(in)               :: tb    !! end of the interval
         real(wp), dimension(:), intent(in) :: ya    !! state at `ta`
         real(wp), dimension(:), intent(in) :: yb    !! state at `tb`
         real(wp), intent(in)               :: ga    !! event function at `ta`
         real(wp), intent(in)               :: gb    !! event function at `tb`
         logical, intent(in)                :: mask  !! if the function is to be ignored
         logical, intent(out)               :: found !! if a root was found
                                                     !! (or there was an error locating it)

         real(wp) :: tzero   !! time where an event occurs in `[ta,tb]`
         real(wp) :: gzero   !! event function value at `tzero`
         integer :: iflag    !! root finder status flag

         found = .false.
         if (mask) return
         if (ga*gb > 0.0_wp) return  ! no change in sign of the event function

         ! the user bracket function can impose
         ! additional constraints on the root, so
         ! we check that now if it was associated:
         if (associated(me%bracket)) then
            if (.not. me%bracket(ta, tb, ya, yb, ga, gb)) return
         end if

         ! root somewhere on [ta,tb]
         ! call the root finder:
         found = .true.
         call root_scalar(root_method_brent, &
                          fun=zeroin_func, &
                          ax=ta, &
                          bx=tb, &
                          rtol=me%tol, &
                          ftol=me%tol, &
                          xzero=tzero, &
                          fzero=gzero, &
                          iflag=iflag, &
                          fax=ga, &
                          fbx=gb)
         if (iflag == 0) then ! root found at tzero
            idid = 1000 ! root found
            !evaluate again to get the final state for output:
            gval = zeroin_func(tzero)
            t = tzero
            y = yc
         else
            ! unlikely to occur since zeroin is
            ! "guaranteed" to converge, but just in case:
            call report_error('ddeabm_with_event_wrapper', &
                              'Error locating root.', 0, 0)
            idid = -2000 ! no root is found
            t = tb
            y = yb
            gval = gb
            me%resume_pending = .false.
         end if

      end subroutine locate_root

      subroutine check_rest_of_last_step(done)

            !! When continuing after a root was found, check the rest of the
            !! last integration step (`[t_resume, t_saved]`) for events, using
            !! interpolation. Points on this interval are also reported if necessary.
            !! On return, if not `done`, the integration continues from the end of
            !! the last step (`t1`, `y1`, `g1` are set for that).

         implicit none

         logical, intent(out) :: done !! if finished (a root was found or `tmax` was reached)

         logical :: reached_tmax !! if `tmax` is within the last step

         done = .false.

         ! the interval is from the last root to the end of the
         ! last step, or to tmax if it is before that:
         reached_tmax = direction*(me%t_saved - tmax) >= 0.0_wp
         t1 = me%t_resume
         call me%interpolate(t1, y1)
         call me%gfunc(t1, y1, g1)
         if (reached_tmax) then
            t2 = tmax
            call me%interpolate(t2, y2)
         else
            t2 = me%t_saved
            y2 = me%x_saved
         end if
         call me%gfunc(t2, y2, g2)

         ! ignore the function if it is zero at the last root:
         call process_interval(t1, t2, y1, y2, g1, g2, abs(g1) <= me%tol, found)
         if (found) then
            done = .true.
            return
         end if

         if (reached_tmax) then
            ! finished: tmax is within the last step
            t = tmax
            y = y2
            gval = g2
            idid = 3
            me%t_resume = tmax
            me%resume_pending = direction*(me%t_saved - tmax) > 0.0_wp
            done = .true.
            return
         end if

         ! continue the integration from the end of the last step:
         t = me%t_saved
         y = me%x_saved
         t1 = t2
         y1 = y2
         g1 = g2
         me%resume_pending = .false.

      end subroutine check_rest_of_last_step

      function zeroin_func(tc) result(g)

        !! evaluate the g function at `tc`
        !! using interpolation ([[dintp]]).

         implicit none

         real(wp), intent(in)  :: tc  !! current time
         real(wp)             :: g   !! value of event function

         ! interpolate to get the state at tc:
         call me%interpolate(tc, yc)

         ! user defined event function:
         call me%gfunc(tc, yc, g)

      end function zeroin_func

   end subroutine ddeabm_with_event_wrapper

ddeabm_with_event_wrapper Subroutine

Contents

Source Code

private subroutine ddeabm_with_event_wrapper(me, t, y, tmax, tstop, idid, gval, integration_mode, tstep, continue)

See also

Type Bound

Arguments

Calls

Source Code