Wrapper routine for ddeabm, with event finding for multiple functions.
It will integrate until any 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().
g functionNote
After an interrupted task (idid = -1, -2, -3, -4), call
me%resume_integration() to continue. For idid=-2, the
tolerances increased by the code are used from then on.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(ddeabm_with_event_class_vec), | 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 |
||
| real(kind=wp), | intent(in), | optional | :: | tstop |
for some problems it may not be permissible to integrate
past a point |
|
| integer, | intent(out) | :: | idid |
indicator reporting what the code did.
you must monitor this integer variable to
decide what action to take next.
|
||
| real(kind=wp), | intent(out), | dimension(:) | :: | gval |
value of the event functions |
|
| integer, | intent(in), | optional | :: | integration_mode |
Step mode:
1 - normal integration from |
|
| real(kind=wp), | intent(in), | optional | :: | tstep |
Fixed time step to use for reporting and evaluation of the event functions (the points are interpolated within each integrator step). If not present, then default integrator steps are used. The event functions are 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(s) that triggered the previous event are ignored (they are zero at the start), so a second root of the same function within the same integration step will not be found. Roots of the other functions are found. |
subroutine ddeabm_with_event_wrapper_vec(me, t, y, tmax, tstop, idid, gval, integration_mode, tstep, continue) implicit none class(ddeabm_with_event_class_vec), 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 for the !! (`idid-1000`)th `g` function. !! See [[ddeabm]] for other values. real(wp), dimension(:), intent(out) :: gval !! value of the event functions `g(t,x)` at the final time `t`. !! Must be size `ng` (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 functions (the points are !! interpolated within each integrator step). If not !! present, then default integrator steps are used. !! The event functions are 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(s) !! that triggered the previous event are ignored (they are zero !! at the start), so a second root of the same function within !! the same integration step will not be found. Roots of the !! other functions are found. ! local variables: real(wp), dimension(:), allocatable :: g1 !! value of event functions at `t1` real(wp), dimension(:), allocatable :: g2 !! value of event functions at `t2` logical, dimension(:), allocatable :: mask !! event functions to ignore on an interval !! (those that are zero at its start) 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) :: t1 !! initial time of an interval real(wp) :: t2 !! final time of an interval integer :: mode !! local copy of integration_mode integer :: ig !! `gfunc` to compute in `zeroin_func` 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_vec', & '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_vec', & '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_vec', & 'the event function GFUNC has not been associated.', 0, 0) idid = -33 return end if if (size(gval) /= me%n_g_eqns) then call report_error('ddeabm_with_event_wrapper_vec', & 'gval must have the same size as the number of event functions.', 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_vec', & '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_vec', & '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_vec', & 'tstep is too small (or zero) to advance the time.', 0, 0) idid = -33 return end if end if ! work arrays: allocate (g1(me%n_g_eqns), g2(me%n_g_eqns), mask(me%n_g_eqns)) 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 events are 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 functions at the initial point: t1 = t y1 = y call me%gfunc(t1, y1, g1) ! compute all the g functions 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 all event functions 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 the functions that are zero at the initial time ! (or at the root where a previous call stopped): mask = first .and. abs(g1) <= me%tol call process_interval(t1, t2, y1, y2, g1, g2, mask, 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 functions 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 functions !! 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), dimension(:), intent(in) :: ga !! event functions at `ta` real(wp), dimension(:), intent(in) :: gb !! event functions at `tb` logical, dimension(:), intent(in) :: mask !! functions to ignore 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), dimension(me%n_g_eqns) :: gl, gn !! event functions at `tl` and `tn` logical, dimension(me%n_g_eqns) :: 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 roots of the event functions on `[ta,tb]` and, !! if any, locate the earliest one (in the direction of integration). !! 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), dimension(:), intent(in) :: ga !! event functions at `ta` real(wp), dimension(:), intent(in) :: gb !! event functions at `tb` logical, dimension(:), intent(in) :: mask !! functions to ignore logical, intent(out) :: found !! if a root was found !! (or there was an error locating it) logical, dimension(me%n_g_eqns) :: root_found !! if a root was found real(wp) :: tzero !! time where an event occurs in `[ta,tb]` real(wp) :: gzero !! event function value at `tzero` real(wp) :: tbest !! earliest root found so far integer :: i !! counter integer :: ibest !! function with the earliest root integer :: iflag !! `root_scalar` status flag integer :: istat !! for write statement `iostat` character(len=10) :: istr !! for integer to string conversion found = .false. ! find the one(s) with sign change: do i = 1, me%n_g_eqns root_found(i) = (.not. mask(i)) .and. ga(i)*gb(i) <= 0.0_wp end do ! 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 do i = 1, me%n_g_eqns if (root_found(i)) then root_found(i) = me%bracket(i, ta, tb, ya, yb, ga(i), gb(i)) end if end do end if if (.not. any(root_found)) return ! root somewhere on [ta,tb]. have to check all the functions ! where this is true, and select the one with the earlier root. found = .true. ibest = 0 tbest = tb do i = 1, me%n_g_eqns if (.not. root_found(i)) cycle ! Note: it's possible that zeroing in on one root ! will cause another earlier root to be detected? ! This is not currently handled. Only the ones identified ! by the normal stepping process are considered. ig = i ! when calling zeroin, we only need to compute the ith function call root_scalar(root_method_brent, & fun=zeroin_func, & ax=ta, & bx=tb, & rtol=me%tol(i), & ftol=me%tol(i), & xzero=tzero, & fzero=gzero, & iflag=iflag, & fax=ga(i), & fbx=gb(i)) if (iflag == 0) then ! root found at tzero if (ibest == 0 .or. direction*(tzero - tbest) < 0.0_wp) then ! this is an earlier root so use it ibest = i tbest = tzero end if else ! unlikely to occur since zeroin is ! "guaranteed" to converge, but just in case: write (istr, '(I10)', iostat=istat) i call report_error('ddeabm_with_event_wrapper_vec', & 'Error locating root for function '// & trim(adjustl(istr)), 0, 0) end if end do if (ibest > 0) then idid = 1000 + ibest ! root found t = tbest ! state at the root (note: yc is not necessarily at tbest): call me%interpolate(t, y) call me%gfunc(t, y, gval) ! compute all funcs for output else call report_error('ddeabm_with_event_wrapper_vec', & 'Error locating root.', 0, 0) idid = -2000 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(s) that are zero at the last root: mask = abs(g1) <= me%tol call process_interval(t1, t2, y1, y2, g1, g2, mask, 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 `ig`th `g` function at `tc` !! using interpolation ([[dintp]]). implicit none real(wp), intent(in) :: tc !! current time real(wp) :: g !! value of event function real(wp), dimension(me%n_g_eqns) :: gtmp ! interpolate to get the state at tc: call me%interpolate(tc, yc) ! user defined event function: ! [only compute the one we need] call me%gfunc(tc, yc, gtmp, ig) g = gtmp(ig) end function zeroin_func end subroutine ddeabm_with_event_wrapper_vec