Detailled example: Verlet
Work in progress!
The Verlet algorithm
According to Wikipedia, the Verlet algorithm is defined as:
- Calculate $\vec{x}(t + \Delta t) = \vec{x}(t) + \vec{v}(t) \Delta t + \tfrac12 \vec{a}(t) \Delta t^2$.
- Derive $\vec{a}(t + \Delta t)$ from the interaction potential using $\vec{x}(t + \Delta t)$.
- Calculate $\vec{v}(t + \Delta t) = \vec{v}(t) + \tfrac12 \big(\vec{a}(t) + \vec{a}(t + \Delta t)\big)\Delta t$.
The Verlet operator is represented by the type propagator_verlet_t which extends propagator_t:
Note, that the class definition does not add anything to the propagator_t class. The only differences are the definition of the operations, and the overloaded constructor.
For the Verlet propagator, we need to define the following operations:
! Specific verlet propagation operations identifiers
character(len=30), public, parameter :: &
VERLET_START = 'VERLET_START', &
VERLET_FINISH = 'VERLET_FINISH', &
VERLET_UPDATE_POS = 'VERLET_UPDATE_POS', &
VERLET_COMPUTE_ACC = 'VERLET_COMPUTE_ACC', &
VERLET_COMPUTE_VEL = 'VERLET_COMPUTE_VEL'
! Specific verlet propagation operations
type(algorithmic_operation_t), public, parameter :: &
OP_VERLET_START = algorithmic_operation_t(VERLET_START, 'Starting Verlet propagation'), &
OP_VERLET_FINISH = algorithmic_operation_t(VERLET_FINISH, 'Finishing Verlet propagation'), &
OP_VERLET_UPDATE_POS = algorithmic_operation_t(VERLET_UPDATE_POS, 'Propagation step - Updating positions'), &
OP_VERLET_COMPUTE_ACC = algorithmic_operation_t(VERLET_COMPUTE_ACC, 'Propagation step - Computing acceleration'), &
OP_VERLET_COMPUTE_VEL = algorithmic_operation_t(VERLET_COMPUTE_VEL, 'Propagation step - Computing velocity')
These are used to define the algorithm, which is done in the constructor of the propagator:
function propagator_verlet_constructor(dt) result(this)
FLOAT, intent(in) :: dt
type(propagator_verlet_t), pointer :: this
PUSH_SUB(propagator_verlet_constructor)
SAFE_ALLOCATE(this)
this%start_step = OP_VERLET_START
this%final_step = OP_VERLET_FINISH
call this%add_operation(OP_VERLET_UPDATE_POS)
call this%add_operation(OP_UPDATE_INTERACTIONS)
call this%add_operation(OP_VERLET_COMPUTE_ACC)
call this%add_operation(OP_VERLET_COMPUTE_VEL)
call this%add_operation(OP_FINISHED)
! Verlet has only one algorithmic step
this%algo_steps = 1
this%dt = dt
POP_SUB(propagator_verlet_constructor)
end function propagator_verlet_constructor
The algorithm also uses steps, which are not specific to the Verlet algorithm, and are defined in the propagator_oct_m module.
Note, the difference between OP_VERLET_FINISH and OP_FINISHED. The former denotes a specific step, to be taken at the end of one time step,
while the latter generally denotes the end of a time step.
As can be seen in this example, the definition of the propagator in terms of the algorithmic operations is a direct translation of the algorithm, shown above.
Implementation of the steps
So far, the propagator is only defined in an abstract way. It is down to the individual systems (or better, their programmers) to implement the individual steps of the algorithm, in terms of the dynamic variables for that system. An easy examample to demonstrate this is the classical particle, as implemented in classical_particles_t
This describes an extremely simple system, consisting of a set of classical, neutral particles. The dynamic variables (i.e. the state) of the system are the positions, velocities and accelerations.
One ‘‘tick’’ of the propagator is defined in the function classical_particle_do_td. As can be seen in the code below, this function implements all possible algorithmic steps for all propagators, allowed for that system.
The timeline explained
The Verlet algorithm is a good (because simple) example to illustrate how the systems and the interaction are updated as the code progresses through the main loop.
This graph illustrates how the state machine is stepping through the algorithm. Each system is picking the next algorithmic step from the propagator. For the containers (i.e. ‘‘root’’ and ‘‘earth’'), the only steps are ‘‘Updating interactions’’ and ‘‘Finished’’. The real systems, on the other hand, are progressing orderly through the operations, defined in the propagator.