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`
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 event function. If not present, then default integrator steps are used.
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).

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`
      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 event function. If not present,
                                             !! then default integrator steps are used.
      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).

      ! 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
      real(wp) :: tzero       !! time where an event occurs in `[t1,t2]`
      integer :: iflag        !! root finder status flag
      logical :: first        !! flag for the first step
      integer :: mode         !! local copy of integration_mode
      logical :: fixed_step   !! if using the fixed step size `tstep`
      logical :: last         !! for fixed step size: the last step
      real(wp) :: dt          !! for fixed step size: actual signed time step
      real(wp) :: direction   !! direction of integration for
                            !! fixed step size: `+1: dt>=0`, `-1: dt<0`
      logical :: continuing   !! local copy of optional `continue` argument
      real(wp) :: first_dt    !! for fixed step size: the `dt` for the
                            !! first step  (can differ from `dt` if
                            !! continuing from a previous event solve)
      logical :: root_found   !! if a root was found

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

      ! 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

      ! if we are reporting the output at a fixed step size:
      fixed_step = present(tstep)
      if (fixed_step) then
         direction = sign(1.0_wp, tmax - t)
         dt = direction*abs(tstep)
         first_dt = dt
         last = .false.
      end if

      if (continuing) then
         ! if continuing, then we reset the t,y inputs
         ! to the values from the last successful step
         ! (note than an event may have been found in the
         ! last call, so the input values are not correct)
         if (fixed_step) then
            ! we need the first step to be from the
            ! last reported point, not the last
            ! successful step:
            first_dt = (t + dt) - me%t_saved
            ! make sure not a 0 dt or in the wrong direction:
            if (first_dt == 0.0_wp .or. sign(1.0_wp, first_dt) /= sign(1.0_wp, dt)) first_dt = dt
         end if
         t = me%t_saved
         y = me%x_saved
         me%info(1) = 1   ! necessary? (does not seem to matter)
      end if

      !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

      !make sure the event function was set:
      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

      !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)
      if (fixed_step) then
         me%info(3) = 0  !return after integrating to the specified (fixed step) time
      else
         me%info(3) = 1  !intermediate output mode: return after each default step
      end if

      !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

      !make a copy of the tols, since the routine might change them:
      me%rtol_tmp = me%rtol
      me%atol_tmp = me%atol

      !evaluate the event function at the initial point:
      first = .true.  !to avoid catching a root at the initial time
      t1 = t
      y1 = y
      call me%gfunc(t1, y1, g1)
      if (mode == 2 .and. .not. continuing) call me%report(t, y)  !initial point

      do

         if (fixed_step) then
            ! take one step to t2 and return
            if (first) then
               t2 = t + first_dt
            else
               t2 = t + dt
            end if
            last = direction*(tmax - t2) <= 0.0_wp  ! if last point
            if (last) t2 = tmax  ! adjust last time step if necessary
         else
            t2 = tmax  ! take a step in the direction of tmax and return
         end if

         !call the lower-level routine:
         call me%ddeabm(neq=me%neq, &
                        t=t, &
                        y=y, &
                        tout=t2, &
                        info=me%info, &
                        rtol=me%rtol_tmp, &
                        atol=me%atol_tmp, &
                        idid=idid)

         ! if there was an integration error (see ddeabm or idid codes):
         if (idid < 0) exit

         !if there was a user-triggered error:
         if (me%error) then
            idid = -1000
            exit
         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)

         if (first .and. abs(g1) <= me%tol) then

            ! ignore a root at the initial time (first point)
            ! (or if continuing integration from a previous root)

            ! warning: is this check always sufficient? or
            !          is it possible for zeroin to have
            !          converged to a root outside this tol?

         else if (g1*g2 <= 0.0_wp) then ! change in sign of the event function

            root_found = .true.

            ! 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
               root_found = me%bracket(t1, t2, y1, y2, g1, g2)
            end if

            if (root_found) then
               ! root somewhere on [t1,t2]
               ! call the root finder:
               call root_scalar(root_method_brent, &
                                fun=zeroin_func, &
                                ax=t1, &
                                bx=t2, &
                                rtol=me%tol, &
                                ftol=me%tol, &
                                xzero=tzero, &
                                fzero=gval, &
                                iflag=iflag, &
                                fax=g1, &
                                fbx=g2)
               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
                  if (mode == 2) call me%report(t, y) ! report this point
                  return
               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
                  return
               end if
            end if

         end if

         ! report this point:
         if (mode == 2) call me%report(t, y)

         ! integration is finished (see ddeabm or idid codes):
         if (fixed_step) then
            if (last) return
         else
            if (idid == 2 .or. idid == 3) return
         end if

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

      end do

   contains

      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