ddeabm_with_event_class_vec Derived Type

type, public, extends(ddeabm_class) :: ddeabm_with_event_class_vec

A version of ddeabm_class with event location (root finding) of a vector function (it will stop on the first root it encounters)

Call initialize_event to set up the integration.

See also: ddeabm_with_event_class


Inherits

type~~ddeabm_with_event_class_vec~~InheritsGraph type~ddeabm_with_event_class_vec ddeabm_with_event_class_vec type~ddeabm_class ddeabm_class type~ddeabm_with_event_class_vec->type~ddeabm_class

Components

Type Visibility Attributes Name Initial
integer, private :: n_g_eqns = 0

number of gfunc event equations
real(kind=wp), private, dimension(:), allocatable :: tol

tolerance for root finding. this is a vector (size n_g_eqns)
real(kind=wp), private :: t_saved = 0.0_wp

time of last successful step (used to continue after a root finding)
real(kind=wp), private, dimension(:), allocatable :: x_saved

state of last successful step. size (neq). (used to continue after a root finding)
real(kind=wp), private :: t_resume = 0.0_wp

time from which to resume event checking within the last step (the location of the last root found)
logical, private :: resume_pending = .false.

if the interval [t_resume, t_saved] of the last step has not yet been checked for events
procedure(event_func_vec), private, pointer :: gfunc => null()

event function: g(t,x)=0 at event
procedure(bracket_func_vec), private, pointer :: bracket => null()

a function for determining if the root is bracketed. the default just checks for a sign change in the event function. the user can specify this to use other conditions.

Type-Bound Procedures

procedure, public, non_overridable :: initialize => ddeabm_initialize

  • private subroutine ddeabm_initialize(me, neq, maxnum, df, rtol, atol, report, initial_step_mode, initial_step_size, max_step_size, max_order)

    Author
    Jacob Williams

    Initialize the ddeabm_class, and set the variables that cannot be changed during a problem.

    Arguments

    Type IntentOptional Attributes Name
    class(ddeabm_class), intent(inout) :: me
    integer, intent(in) :: neq

    number of equations to be integrated
    integer, intent(in) :: maxnum

    maximum number of integration steps allowed
    procedure(deriv_func) :: df

    derivative function dx/dt
    real(kind=wp), intent(in), dimension(:) :: rtol

    relative tolerance for integration (see ddeabm)
    real(kind=wp), intent(in), dimension(:) :: atol

    absolution tolerance for integration (see ddeabm)
    procedure(report_func), optional :: report

    reporting function
    integer, intent(in), optional :: initial_step_mode

    how to choose the initial step h:

    Read more…
    real(kind=wp), intent(in), optional :: initial_step_size

    for initial_step_mode=3
    real(kind=wp), intent(in), optional :: max_step_size

    maximum step size magnitude (>0). If not present, the step size is not limited.
    integer, intent(in), optional :: max_order

    maximum order of the method (1 to ddeabm_max_order_limit). If not present, ddeabm_default_max_order (12) is used. A lower value can improve stability on some problems. A higher value can reduce the number of steps for very tight tolerances when using higher precision (REAL128).

procedure, public, non_overridable :: integrate => ddeabm_wrapper

  • private subroutine ddeabm_wrapper(me, t, y, tout, tstop, idid, integration_mode, tstep)

    Author
    Jacob Williams

    Wrapper routine for ddeabm.

    Read more…

    Arguments

    Type IntentOptional Attributes Name
    class(ddeabm_class), intent(inout) :: me
    real(kind=wp), intent(inout) :: t
    real(kind=wp), intent(inout), dimension(me%neq) :: y
    real(kind=wp), intent(in) :: tout

    time at which a solution is desired.
    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 tout 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.
    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 integration_mode=2. If not present, then default integrator steps are used. If integration_mode=1, then this is ignored.

procedure, public, non_overridable :: destroy => destroy_ddeabm

procedure, public, non_overridable :: first_call => ddeabm_new_problem

  • private subroutine ddeabm_new_problem(me)

    Author
    Jacob Williams

    Call this to indicate that a new problem is being solved. It sets info(1) = 0 (see ddeabm documentation).

    Arguments

    Type IntentOptional Attributes Name
    class(ddeabm_class), intent(inout) :: me

procedure, public, non_overridable :: interpolate => ddeabm_interp

state interpolation function

  • private subroutine ddeabm_interp(me, tc, yc, istat)

    Interpolation function. Can be used for dense output after a step. It calls the low-level routine dintp.

    Arguments

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

    point at which solution is desired
    real(kind=wp), intent(out), dimension(me%neq) :: yc

    interpolated state at tc
    integer, intent(out), optional :: istat

    status flag:

    Read more…

procedure, public, non_overridable :: stop_integration => ddeabm_stop_integration

user can call this in df routine to stop the integration

  • private subroutine ddeabm_stop_integration(me)

    Author
    Jacob Williams

    Call this to abort the integration if there is an error. (Can be called in df by the user.) The step will be considered invalid and not returned.

    Read more…

    Arguments

    Type IntentOptional Attributes Name
    class(ddeabm_class), intent(inout) :: me

procedure, public, non_overridable :: resume_integration => ddeabm_resume_integration

call this after an interrupted task (idid = -1, -2, -3, -4) to continue the integration (equivalent to setting info(1)=1)

  • private subroutine ddeabm_resume_integration(me)

    Author
    Jacob Williams

    Call this to continue the integration after an interrupted task (idid = -1, -2, -3, or -4). This is equivalent to setting info(1) = 1 as described in the ddeabm documentation. It has no effect if the previous task was not interrupted.

    Read more…

    Arguments

    Type IntentOptional Attributes Name
    class(ddeabm_class), intent(inout) :: me

procedure, public, non_overridable :: initialize_event => ddeabm_with_event_initialize_vec

  • private subroutine ddeabm_with_event_initialize_vec(me, neq, maxnum, df, rtol, atol, ng, g, root_tol, report, bracket, initial_step_mode, initial_step_size, max_step_size, max_order)

    Author
    Jacob Williams

    Initialize ddeabm_with_event_class_vec class, and set the variables that cannot be changed during a problem.

    Read more…

    Arguments

    Type IntentOptional Attributes Name
    class(ddeabm_with_event_class_vec), intent(inout) :: me
    integer, intent(in) :: neq

    number of equations to be integrated
    integer, intent(in) :: maxnum

    maximum number of integration steps allowed
    procedure(deriv_func) :: df

    derivative function

    real(kind=wp), intent(in), dimension(:) :: rtol

    relative tolerance for integration (see ddeabm)
    real(kind=wp), intent(in), dimension(:) :: atol

    absolution tolerance for integration (see ddeabm)
    integer, intent(in) :: ng

    number of event functions in g
    procedure(event_func_vec) :: g

    Event function . This should be a smooth function than can have values and . When (within the tolerance), then a root has been located and the integration will stop.
    real(kind=wp), intent(in), dimension(:) :: root_tol

    tolerance for the root finding. This should be sized ng or 1 (in which case the value is used for all elements)
    procedure(report_func), optional :: report

    reporting function
    procedure(bracket_func_vec), optional :: bracket

    root bracketing function. if not present, the default is used.
    integer, intent(in), optional :: initial_step_mode

    how to choose the initial step h:

    Read more…
    real(kind=wp), intent(in), optional :: initial_step_size

    for initial_step_mode=3
    real(kind=wp), intent(in), optional :: max_step_size

    maximum step size magnitude (>0). If not present, the step size is not limited.
    integer, intent(in), optional :: max_order

    maximum order of the method (1 to ddeabm_max_order_limit). If not present, ddeabm_default_max_order (12) is used. A lower value can improve stability on some problems. A higher value can reduce the number of steps for very tight tolerances when using higher precision (REAL128).

procedure, public, non_overridable :: integrate_to_event => ddeabm_with_event_wrapper_vec

main routine for integration to an event (note that integrate can also still be used from ddeabm_class)

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

    Author
    Jacob Williams

    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)

    Read more…

    Arguments

    Type IntentOptional 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 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 for the (idid-1000)th g function. See ddeabm for other values.
    real(kind=wp), intent(out), dimension(:) :: 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(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.

Source Code

   type, extends(ddeabm_class), public :: ddeabm_with_event_class_vec

        !! A version of [[ddeabm_class]] with event location (root finding)
        !! of a vector function (it will stop on the first root it encounters)
        !!
        !! Call *initialize_event* to set up the integration.
        !!
        !! See also: [[ddeabm_with_event_class]]

      private

      integer :: n_g_eqns = 0 !! number of `gfunc` event equations

      real(wp), dimension(:), allocatable :: tol   !! tolerance for root finding.
                                                   !! this is a vector (size `n_g_eqns`)

      real(wp) :: t_saved = 0.0_wp            !! time of last successful step
                                              !! (used to continue after a root finding)
      real(wp), dimension(:), allocatable :: x_saved !! state of last successful step.  size `(neq)`.
                                                     !! (used to continue after a root finding)
      real(wp) :: t_resume = 0.0_wp           !! time from which to resume event checking
                                              !! within the last step (the location of the
                                              !! last root found)
      logical :: resume_pending = .false.     !! if the interval `[t_resume, t_saved]` of the
                                              !! last step has not yet been checked for events

      procedure(event_func_vec), pointer :: gfunc => null()  !! event function: g(t,x)=0 at event

      procedure(bracket_func_vec), pointer :: bracket => null()
            !! a function for determining if the root is bracketed.
            !! the default just checks for a sign change in the event function.
            !! the user can specify this to use other conditions.

   contains

      private

      procedure, non_overridable, public :: initialize_event => ddeabm_with_event_initialize_vec
      procedure, non_overridable, public :: integrate_to_event => ddeabm_with_event_wrapper_vec
            !! main routine for integration to an event
            !! (note that *integrate* can also still be used from [[ddeabm_class]])

   end type ddeabm_with_event_class_vec