ddeabm_class – ddeabm

type, public :: ddeabm_class

The main DDEABM integration class.

Inherited by

Help

Components

Type	Visibility	Attributes		Name		Initial
integer,	private		::	neq	=	0	the number of (first order) differential equations to be integrated (>=0)
integer,	private		::	maxnum	=	500	the expense of solving the problem is monitored by counting the number of steps attempted. when this exceeds `maxnum`, the counter is reset to zero and the user is informed about possible excessive work.
procedure(deriv_func),	private,	pointer	::	df	=>	null()	to define the system of first order differential equations which is to be solved. for the given values of `x` and the vector `u()=(u(1),u(2),...,u(neq))`, the subroutine must evaluate the `neq` components of the system of differential equations `du/dx=df(x,u)` and store the derivatives in the array `uprime()`, that is, `uprime(i) = du(i)/dx` for equations `i=1,...,neq`.
procedure(report_func),	private,	pointer	::	report	=>	null()	a function for reporting the integration steps (optional). this can be associated in ddeabm_initialize or ddeabm_with_event_initialize, and is enabled using the `mode` argument in ddeabm_wrapper or ddeabm_with_event_wrapper
logical,	private		::	scalar_tols	=	.true.
real(kind=wp),	private,	allocatable, dimension(:)	::	rtol			the user input rel tol
real(kind=wp),	private,	allocatable, dimension(:)	::	atol			the user input abs tol
real(kind=wp),	private,	allocatable, dimension(:)	::	rtol_tmp			rel tol used internally
real(kind=wp),	private,	allocatable, dimension(:)	::	atol_tmp			abs tol used internally
integer,	private,	dimension(4)	::	info	=	0	info array
integer,	private		::	initial_step_mode	=	1	how to choose the initial step `h`: Use dhstrt. Use the older (quicker) algorithm. Use a user-specified value.
real(kind=wp),	private		::	initial_step_size	=	0.0_wp	when `initial_step_mode=3`, the `h` value to use. (for the other modes, this will also contain the value used)
integer,	private		::	icount	=	0	replaced `iwork(liw)` in original code
real(kind=wp),	private		::	tprev	=	0.0_wp	replaced `rwork(itstar)` in original code
real(kind=wp),	private,	dimension(:), allocatable	::	two			$2^i$ array
real(kind=wp),	private,	dimension(:), allocatable	::	gstr			$\| \gamma^{*}_i \|$ array
real(kind=wp),	private,	allocatable, dimension(:)	::	ypout			`(neq)` contains the approximate derivative of the solution component `y(i)`. in ddeabm, it is obtained by calling subroutine `df` to evaluate the differential equation using `t` and `y(*)` when `idid=1` or `2`, and by interpolation when `idid=3`.
real(kind=wp),	private,	allocatable, dimension(:)	::	yp			`(neq)`
real(kind=wp),	private,	allocatable, dimension(:)	::	yy			`(neq)`
real(kind=wp),	private,	allocatable, dimension(:)	::	wt			`(neq)`
real(kind=wp),	private,	allocatable, dimension(:)	::	p			`(neq)`
real(kind=wp),	private,	allocatable, dimension(:, :)	::	phi			`(neq,16)`
real(kind=wp),	private,	dimension(12)	::	alpha	=	0.0_wp
real(kind=wp),	private,	dimension(12)	::	beta	=	0.0_wp
real(kind=wp),	private,	dimension(12)	::	psi	=	0.0_wp
real(kind=wp),	private,	dimension(12)	::	v	=	0.0_wp
real(kind=wp),	private,	dimension(12)	::	w	=	0.0_wp
real(kind=wp),	private,	dimension(13)	::	sig	=	0.0_wp
real(kind=wp),	private,	dimension(13)	::	g	=	0.0_wp
real(kind=wp),	private,	dimension(11)	::	gi	=	0.0_wp
real(kind=wp),	private		::	h	=	0.0_wp	the step size `h` to be attempted on the next step.
real(kind=wp),	private		::	eps	=	0.0_wp	if the tolerances have been increased by the code (`idid=-2`), they were multiplied by `eps`.
real(kind=wp),	private		::	x	=	0.0_wp	contains the current value of the independent variable, i.e. the farthest point integration has reached. this will be different from `t` only when interpolation has been performed (`idid=3`).
real(kind=wp),	private		::	xold	=	0.0_wp
real(kind=wp),	private		::	hold	=	0.0_wp
real(kind=wp),	private		::	told	=	0.0_wp
real(kind=wp),	private		::	delsgn	=	0.0_wp
real(kind=wp),	private		::	tstop	=	0.0_wp	(see `info(4)`)
real(kind=wp),	private		::	twou	=	0.0_wp	machine dependent parameter
real(kind=wp),	private		::	fouru	=	0.0_wp	machine dependent parameter
logical,	private		::	start	=	.false.	indicator for dsteps code
logical,	private		::	phase1	=	.false.	indicator for dsteps code
logical,	private		::	nornd	=	.false.	indicator for dsteps code
logical,	private		::	stiff	=	.false.	indicator for stiffness detection
logical,	private		::	intout	=	.false.	indicator for intermediate-output
integer,	private		::	ns	=	0
integer,	private		::	kord	=	0
integer,	private		::	kold	=	0
integer,	private		::	init	=	0	initialization indicator
integer,	private		::	ksteps	=	0	counter for attempted steps
integer,	private		::	kle4	=	0	step counter for stiffness detection
integer,	private		::	iquit	=	0	termination flag
integer,	private		::	kprev	=	0
integer,	private		::	ivc	=	0
integer,	private,	dimension(10)	::	iv	=	0
integer,	private		::	kgi	=	0
logical,	private		::	error	=	.false.	will be set in `stop_integration`. indicates an error and to abort the integration

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)

Author: Jacob Williams

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

Arguments

Type	Intent	Optional	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`

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.

Arguments

Type	Intent	Optional	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

private subroutine destroy_ddeabm(me)

Author: Jacob Williams

Destructor for ddeabm_class.

Arguments

Type	Intent	Optional	Attributes		Name
class(ddeabm_class),	intent(out)			::	me

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	Intent	Optional	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)

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

Arguments

Type	Intent	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`

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.

Arguments

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

procedure, private, non_overridable :: ddeabm

private subroutine ddeabm(me, neq, t, y, tout, info, rtol, atol, idid)

solve an initial value problem in ordinary differential equations using an adams-bashforth method.

Arguments

Type	Intent	Attributes		Name
class(ddeabm_class),	intent(inout)		::	me
integer,	intent(in)		::	neq	the number of (first order) differential equations to be integrated. (neq >= 1)
real(kind=wp),	intent(inout)		::	t	value of the independent variable. on input, set it to the initial point of the integration. the code changes its value. on output, the solution was successfully advanced to the output value of t.
real(kind=wp),	intent(inout),	dimension(neq)	::	y
real(kind=wp),	intent(in)		::	tout
integer,	intent(inout),	dimension(4)	::	info
real(kind=wp),	intent(inout),	dimension(:)	::	rtol
real(kind=wp),	intent(inout),	dimension(:)	::	atol
integer,	intent(out)		::	idid

procedure, private, non_overridable :: ddes

private subroutine ddes(me, neq, t, y, tout, info, rtol, atol, idid, ypout, yp, yy, wt, p, phi, alpha, beta, psi, v, w, sig, g, gi, h, eps, x, xold, hold, told, delsgn, tstop, twou, fouru, start, phase1, nornd, stiff, intout, ns, kord, kold, init, ksteps, kle4, iquit, kprev, ivc, iv, kgi)

ddeabm merely allocates storage for ddes to relieve the user of the inconvenience of a long call list. consequently ddes is used as described in the comments for ddeabm.

Arguments

Type	Intent	Attributes		Name
class(ddeabm_class),	intent(inout)		::	me
integer,	intent(in)		::	neq
real(kind=wp),	intent(inout)		::	t
real(kind=wp),	intent(inout),	dimension(neq)	::	y
real(kind=wp),	intent(in)		::	tout
integer,	intent(inout),	dimension(4)	::	info
real(kind=wp),	intent(inout),	dimension(:)	::	rtol
real(kind=wp),	intent(inout),	dimension(:)	::	atol
integer,	intent(inout)		::	idid
real(kind=wp),	intent(inout),	dimension(neq)	::	ypout
real(kind=wp),	intent(inout),	dimension(neq)	::	yp
real(kind=wp),	intent(inout),	dimension(neq)	::	yy
real(kind=wp),	intent(inout),	dimension(neq)	::	wt
real(kind=wp),	intent(inout),	dimension(neq)	::	p
real(kind=wp),	intent(inout),	dimension(neq, 16)	::	phi
real(kind=wp),	intent(inout),	dimension(12)	::	alpha
real(kind=wp),	intent(inout),	dimension(12)	::	beta
real(kind=wp),	intent(inout),	dimension(12)	::	psi
real(kind=wp),	intent(inout),	dimension(12)	::	v
real(kind=wp),	intent(inout),	dimension(12)	::	w
real(kind=wp),	intent(inout),	dimension(13)	::	sig
real(kind=wp),	intent(inout),	dimension(13)	::	g
real(kind=wp),	intent(inout),	dimension(11)	::	gi
real(kind=wp),	intent(inout)		::	h
real(kind=wp),	intent(inout)		::	eps
real(kind=wp),	intent(inout)		::	x
real(kind=wp),	intent(inout)		::	xold
real(kind=wp),	intent(inout)		::	hold
real(kind=wp),	intent(inout)		::	told
real(kind=wp),	intent(inout)		::	delsgn
real(kind=wp),	intent(inout)		::	tstop
real(kind=wp),	intent(inout)		::	twou
real(kind=wp),	intent(inout)		::	fouru
logical,	intent(inout)		::	start
logical,	intent(inout)		::	phase1
logical,	intent(inout)		::	nornd
logical,	intent(inout)		::	stiff
logical,	intent(inout)		::	intout
integer,	intent(inout)		::	ns
integer,	intent(inout)		::	kord
integer,	intent(inout)		::	kold
integer,	intent(inout)		::	init
integer,	intent(inout)		::	ksteps
integer,	intent(inout)		::	kle4
integer,	intent(inout)		::	iquit
integer,	intent(inout)		::	kprev
integer,	intent(inout)		::	ivc
integer,	intent(inout),	dimension(10)	::	iv
integer,	intent(inout)		::	kgi

procedure, private, non_overridable :: dhstrt

private subroutine dhstrt(me, neq, a, b, y, yprime, etol, morder, small, big, spy, pv, yp, sf, h)

Computes a starting step size to be used in solving initial value problems in ordinary differential equations.

Arguments

Type	Intent	Attributes		Name
class(ddeabm_class),	intent(inout)		::	me
integer,	intent(in)		::	neq	the number of (first order) differential equations to be integrated.
real(kind=wp),	intent(in)		::	a	the initial point of integration.
real(kind=wp),	intent(in)		::	b	a value of the independent variable used to define the direction of integration. a reasonable choice is to set `b` to the first point at which a solution is desired. you can also use `b`, if necessary, to restrict the length of the first integration step because the algorithm will not compute a starting step length which is bigger than `abs(b-a)`, unless `b` has been chosen too close to `a`. (it is presumed that dhstrt has been called with `b` different from `a` on the machine being used. also see the discussion about the parameter `small`.)
real(kind=wp),	intent(in),	dimension(neq)	::	y	the vector of initial values of the `neq` solution components at the initial point `a`.
real(kind=wp),	intent(in),	dimension(neq)	::	yprime	the vector of derivatives of the `neq` solution components at the initial point `a`. (defined by the differential equations in subroutine `me%df`)
real(kind=wp),	intent(in),	dimension(neq)	::	etol	the vector of error tolerances corresponding to the `neq` solution components. it is assumed that all elements are positive. following the first integration step, the tolerances are expected to be used by the integrator in an error test which roughly requires that `abs(local error) <= etol` for each vector component.
integer,	intent(in)		::	morder	the order of the formula which will be used by the initial value method for taking the first integration step.
real(kind=wp),	intent(in)		::	small	a small positive machine dependent constant which is used for protecting against computations with numbers which are too small relative to the precision of floating point arithmetic. `small` should be set to (approximately) the smallest positive double precision number such that `(1.0_wp+small) > 1.0_wp`. on the machine being used. the quantity `small*(3.0_wp/8.0_wp)` is used in computing increments of variables for approximating derivatives by differences. also the algorithm will not compute a starting step length which is smaller than `100.0_wpsmall*abs(a)`.
real(kind=wp),	intent(in)		::	big	a large positive machine dependent constant which is used for preventing machine overflows. a reasonable choice is to set big to (approximately) the square root of the largest double precision number which can be held in the machine.
real(kind=wp),	intent(inout),	dimension(neq)	::	spy	work array of length `neq` which provide the routine with needed storage space.
real(kind=wp),	intent(inout),	dimension(neq)	::	pv	work array of length `neq` which provide the routine with needed storage space.
real(kind=wp),	intent(inout),	dimension(neq)	::	yp	work array of length `neq` which provide the routine with needed storage space.
real(kind=wp),	intent(inout),	dimension(neq)	::	sf	work array of length `neq` which provide the routine with needed storage space.
real(kind=wp),	intent(out)		::	h	appropriate starting step size to be attempted by the differential equation method.

procedure, private, non_overridable :: dsteps

private subroutine dsteps(me, neqn, y, x, h, eps, wt, start, hold, k, kold, crash, phi, p, yp, psi, alpha, beta, sig, v, w, g, phase1, ns, nornd, ksteps, twou, fouru, xold, kprev, ivc, iv, kgi, gi)

integrate a system of first order ordinary differential equations one step.

Arguments

Type	Intent	Attributes		Name
class(ddeabm_class),	intent(inout)		::	me
integer,	intent(in)		::	neqn	number of equations to be integrated
real(kind=wp),	intent(inout),	dimension(neqn)	::	y	solution vector at x
real(kind=wp),	intent(inout)		::	x	independent variable
real(kind=wp),	intent(inout)		::	h	appropriate step size for next step. normally determined by code
real(kind=wp),	intent(inout)		::	eps	local error tolerance
real(kind=wp),	intent(inout),	dimension(neqn)	::	wt	vector of weights for error criterion
logical,	intent(inout)		::	start	logical variable set .true. for first step, .false. otherwise
real(kind=wp),	intent(inout)		::	hold	step size used for last successful step
integer,	intent(inout)		::	k	appropriate order for next step (determined by code)
integer,	intent(inout)		::	kold	order used for last successful step
logical,	intent(inout)		::	crash	logical variable set .true. when no step can be taken, .false. otherwise.
real(kind=wp),	intent(inout),	dimension(neqn, 16)	::	phi	required for the interpolation subroutine dintp
real(kind=wp),	intent(inout),	dimension(neqn)	::	p	required for the interpolation subroutine dintp
real(kind=wp),	intent(inout),	dimension(neqn)	::	yp	derivative of solution vector at x after successful step
real(kind=wp),	intent(inout),	dimension(12)	::	psi
real(kind=wp),	intent(inout),	dimension(12)	::	alpha	required for the interpolation subroutine dintp
real(kind=wp),	intent(inout),	dimension(12)	::	beta
real(kind=wp),	intent(inout),	dimension(13)	::	sig
real(kind=wp),	intent(inout),	dimension(12)	::	v
real(kind=wp),	intent(inout),	dimension(12)	::	w	required for the interpolation subroutine dintp
real(kind=wp),	intent(inout),	dimension(13)	::	g	required for the interpolation subroutine dintp
logical,	intent(inout)		::	phase1
integer,	intent(inout)		::	ns	number of dsteps taken with size h
logical,	intent(inout)		::	nornd
integer,	intent(inout)		::	ksteps	counter on attempted steps
real(kind=wp),	intent(inout)		::	twou	2.*u where u is machine unit roundoff quantity
real(kind=wp),	intent(inout)		::	fouru	4.*u where u is machine unit roundoff quantity
real(kind=wp),	intent(inout)		::	xold	required for the interpolation subroutine dintp
integer,	intent(inout)		::	kprev
integer,	intent(inout)		::	ivc	required for the interpolation subroutine dintp
integer,	intent(inout),	dimension(10)	::	iv	required for the interpolation subroutine dintp
integer,	intent(inout)		::	kgi	required for the interpolation subroutine dintp
real(kind=wp),	intent(inout),	dimension(11)	::	gi	required for the interpolation subroutine dintp

Source Code

   type, public :: ddeabm_class

      !! The main DDEABM integration class.

      private

      integer :: neq = 0
            !! the number of (first order) differential equations to be integrated (>=0)

      integer :: maxnum = 500
            !!  the expense of solving the problem is monitored by counting the
            !!  number of steps attempted. when this exceeds `maxnum`, the counter
            !!  is reset to zero and the user is informed about possible excessive
            !!  work.

      procedure(deriv_func), pointer :: df => null()
            !!  to define the system of first order differential equations
            !!  which is to be solved.  for the given values of `x` and the
            !!  vector `u(*)=(u(1),u(2),...,u(neq))`, the subroutine must
            !!  evaluate the `neq` components of the system of differential
            !!  equations `du/dx=df(x,u)` and store the derivatives in the
            !!  array `uprime(*)`, that is, `uprime(i) = du(i)/dx` for
            !!  equations `i=1,...,neq`.

      procedure(report_func), pointer :: report => null()
            !! a function for reporting the integration steps (optional).
            !! this can be associated in [[ddeabm_initialize]] or
            !! [[ddeabm_with_event_initialize]], and is enabled
            !! using the `mode` argument in [[ddeabm_wrapper]]
            !! or [[ddeabm_with_event_wrapper]]

      !tolerances
      logical :: scalar_tols = .true.
      real(wp), allocatable, dimension(:) :: rtol        !! the user input rel tol
      real(wp), allocatable, dimension(:) :: atol        !! the user input abs tol
      real(wp), allocatable, dimension(:) :: rtol_tmp    !! rel tol used internally
      real(wp), allocatable, dimension(:) :: atol_tmp    !! abs tol used internally

      integer, dimension(4) :: info = 0  !! info array

      ! initial step mode:
      integer :: initial_step_mode = 1    !! how to choose the initial step `h`:
                                          !!
                                          !! 1. Use [[dhstrt]].
                                          !! 2. Use the older (quicker) algorithm.
                                          !! 3. Use a user-specified value.
      real(wp) :: initial_step_size = 0.0_wp !! when `initial_step_mode=3`, the `h` value to use.
                                             !! (for the other modes, this will also contain the value used)

      !variables formerly in work arrays:

      integer :: icount = 0         !! replaced `iwork(liw)` in original code
      real(wp) :: tprev = 0.0_wp    !! replaced `rwork(itstar)` in original code

      !these were fixed in the original code:
      real(wp), dimension(:), allocatable :: two     !! \( 2^i \) array
      real(wp), dimension(:), allocatable :: gstr    !! \( | \gamma^{*}_i | \) array

      real(wp), allocatable, dimension(:)    :: ypout    !! `(neq)` contains the approximate derivative
                                                         !! of the solution component `y(i)`.  in [[ddeabm]], it
                                                         !! is obtained by calling subroutine `df` to
                                                         !! evaluate the differential equation using `t` and
                                                         !! `y(*)` when `idid=1` or `2`, and by interpolation
                                                         !! when `idid=3`.
      real(wp), allocatable, dimension(:)    :: yp       !! `(neq)`
      real(wp), allocatable, dimension(:)    :: yy       !! `(neq)`
      real(wp), allocatable, dimension(:)    :: wt       !! `(neq)`
      real(wp), allocatable, dimension(:)    :: p        !! `(neq)`
      real(wp), allocatable, dimension(:, :)  :: phi      !! `(neq,16)`
      real(wp), dimension(12)               :: alpha = 0.0_wp
      real(wp), dimension(12)               :: beta = 0.0_wp
      real(wp), dimension(12)               :: psi = 0.0_wp
      real(wp), dimension(12)               :: v = 0.0_wp
      real(wp), dimension(12)               :: w = 0.0_wp
      real(wp), dimension(13)               :: sig = 0.0_wp
      real(wp), dimension(13)               :: g = 0.0_wp
      real(wp), dimension(11)               :: gi = 0.0_wp
      real(wp)                             :: h = 0.0_wp  !! the step size `h` to be attempted on the next step.
      real(wp)                             :: eps = 0.0_wp  !! if the tolerances have been increased by the
                                                            !! code (`idid=-2`), they were multiplied by `eps`.
      real(wp)                             :: x = 0.0_wp  !! contains the current value of the
                                                         !! independent variable, i.e. the farthest point
                                                         !! integration has reached. this will be different
                                                         !! from `t` only when interpolation has been
                                                         !! performed (`idid=3`).
      real(wp)                             :: xold = 0.0_wp
      real(wp)                             :: hold = 0.0_wp
      real(wp)                             :: told = 0.0_wp
      real(wp)                             :: delsgn = 0.0_wp
      real(wp)                             :: tstop = 0.0_wp  !! (see `info(4)`)
      real(wp)                             :: twou = 0.0_wp  !! machine dependent parameter
      real(wp)                             :: fouru = 0.0_wp  !! machine dependent parameter
      logical                              :: start = .false. !! indicator for [[dsteps]] code
      logical                              :: phase1 = .false. !! indicator for [[dsteps]] code
      logical                              :: nornd = .false. !! indicator for [[dsteps]] code
      logical                              :: stiff = .false. !! indicator for stiffness detection
      logical                              :: intout = .false. !! indicator for intermediate-output
      integer                              :: ns = 0
      integer                              :: kord = 0
      integer                              :: kold = 0
      integer                              :: init = 0  !! initialization indicator
      integer                              :: ksteps = 0  !! counter for attempted steps
      integer                              :: kle4 = 0  !! step counter for stiffness detection
      integer                              :: iquit = 0  !! termination flag
      integer                              :: kprev = 0
      integer                              :: ivc = 0
      integer, dimension(10)                :: iv = 0
      integer                              :: kgi = 0

      logical :: error = .false.  !! will be set in `stop_integration`.
                                  !! indicates an error and to
                                  !! abort the integration

   contains

      private

      procedure, non_overridable, public :: initialize => ddeabm_initialize
      procedure, non_overridable, public :: integrate => ddeabm_wrapper
      procedure, non_overridable, public :: destroy => destroy_ddeabm
      procedure, non_overridable, public :: first_call => ddeabm_new_problem
      procedure, non_overridable, public :: interpolate => ddeabm_interp !! state interpolation function

      procedure, non_overridable, public :: stop_integration => ddeabm_stop_integration
                !! user can call this in `df` routine to stop the integration

      !support routines:
      procedure, non_overridable :: ddeabm
      procedure, non_overridable :: ddes
      procedure, non_overridable :: dhstrt
      procedure, non_overridable :: dsteps

   end type ddeabm_class

ddeabm_class Derived Type

Contents

Variables

Type-Bound Procedures

Source Code

type, public :: ddeabm_class

Inherited by

Components

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)

Arguments

procedure, public, non_overridable :: integrate => ddeabm_wrapper

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

Arguments

procedure, public, non_overridable :: destroy => destroy_ddeabm

private subroutine destroy_ddeabm(me)

Arguments

procedure, public, non_overridable :: first_call => ddeabm_new_problem

private subroutine ddeabm_new_problem(me)

Arguments

procedure, public, non_overridable :: interpolate => ddeabm_interp

private subroutine ddeabm_interp(me, tc, yc)

Arguments

procedure, public, non_overridable :: stop_integration => ddeabm_stop_integration

private subroutine ddeabm_stop_integration(me)

Arguments

procedure, private, non_overridable :: ddeabm

private subroutine ddeabm(me, neq, t, y, tout, info, rtol, atol, idid)

Arguments

procedure, private, non_overridable :: ddes

private subroutine ddes(me, neq, t, y, tout, info, rtol, atol, idid, ypout, yp, yy, wt, p, phi, alpha, beta, psi, v, w, sig, g, gi, h, eps, x, xold, hold, told, delsgn, tstop, twou, fouru, start, phase1, nornd, stiff, intout, ns, kord, kold, init, ksteps, kle4, iquit, kprev, ivc, iv, kgi)

Arguments

procedure, private, non_overridable :: dhstrt

private subroutine dhstrt(me, neq, a, b, y, yprime, etol, morder, small, big, spy, pv, yp, sf, h)

Arguments

procedure, private, non_overridable :: dsteps

private subroutine dsteps(me, neqn, y, x, h, eps, wt, start, hold, k, kold, crash, phi, p, yp, psi, alpha, beta, sig, v, w, g, phase1, ns, nornd, ksteps, twou, fouru, xold, kprev, ivc, iv, kgi, gi)

Arguments

Source Code