Oscillation equations
The (inertial) Oscillation equations system is a non linear system of pure ODEs and it can be written as:
The frequency f is constant and it is here selected as f=10^-4.
In the space v1-v2 the path of the oscillation must be a circle, but for the frequency selected some ODE solvers are not stable leading to a wrong path, see [1].
[1] Numerical Methods for Fluid Dynamics With Applications to Geophysics, Dale R. Durran, Springer, 2010.
Before apply a FOODIE solver for solving the above system the system itself must be defined accordingly to FOODIE's type_integrand abstract type (upon which each FOODIE's solver is based).
The proposed FOODIE-compatible Oscillation' system definition is:
type, extends(integrand) :: oscillation
!< Oscillation equations field.
!<
!< It is a FOODIE integrand class.
private
integer(I_P) :: dims=0 !< Space dimensions.
integer(I_P) :: steps=0 !< Number of time steps stored.
real(R_P), dimension(:), allocatable :: state !< Solution vector, [1:dims].
real(R_P), dimension(:,:), allocatable :: previous !< Previous steps solution vector, [1:dims,1:steps].
real(R_P) :: f = 0._R_P !< Oscillation frequency (Hz).
contains
! auxiliary methods
procedure, pass(self), public :: init !< Init field.
procedure, pass(self), public :: output !< Extract Oscillation field.
! type_integrand deferred methods
procedure, pass(self), public :: t => dOscillation_dt !< Time derivative, residuals.
procedure, pass(self), public :: update_previous_steps !< Update previous time steps.
procedure, pass(lhs), public :: integrand_multiply_integrand => oscillation_multiply_oscillation !< Oscillation * oscillation.
procedure, pass(lhs), public :: integrand_multiply_real => oscillation_multiply_real !< Oscillation * real.
procedure, pass(rhs), public :: real_multiply_integrand => real_multiply_oscillation !< Real * Oscillation.
procedure, pass(lhs), public :: add => add_oscillation !< Oscillation + Oscillation.
procedure, pass(lhs), public :: assign_integrand => oscillation_assign_oscillation !< Oscillation = Oscillation.
procedure, pass(lhs), public :: assign_real => oscillation_assign_real !< Oscillation = real.
endtype oscillationThis Oscillation type embeds:
- the 3 parameters values;
- the Oscillation' state variable vector;
- 2 auxiliary (helper) methods, not detailed here they being of minor interest;
- all the
type_integranddeferred methods that are necessary for the definition of a valid concrete extension oftype_integrand.
We chose to define the Oscillation' state vector as a rank 1 array containing the state variables values of the current time step. The eventually allocated previous array is a rank 2 array: in the first subscript there is the state space dimension, namely it has dimension equals to 2 containing the dependent state variables v1, v2; in the second subscript are stored the current and the (eventually) previous time steps solution. The current time step solution is always the upper bound element of the second subscript. For example, let us consider to use a 2 time steps solver, the Oscillation' state vector has the following meaning:
-
oscillation%previous(:, 1)=> solution at timen-1 -
oscillation%previous(:, 2)=> solution at timen=> oscillation%state
The time derivative method oscillation%t => dOscillation_dt, namely the residuals function, depends only on the state vector at time step considered, thus it is very simple to implement:
Using the above oscillation type definition it simply corresponds to:
pure function dOscillation_dt(self, n) result(dState_dt)
!---------------------------------------------------------------------------------------------------------------------------------
!< Time derivative of Oscillation field.
!---------------------------------------------------------------------------------------------------------------------------------
class(oscillation), intent(IN) :: self !< Oscillation field.
integer(I_P), optional, intent(IN) :: n !< Time level.
class(integrand), allocatable :: dState_dt !< Oscillation field time derivative.
type(oscillation), allocatable :: delta !< Delta state used as temporary variables.
integer :: dn !< Time level, dummy variable.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
! preparing temporary delta
delta%dims = self%dims
delta%steps = self%steps
delta%state = self%state
if (allocated(self%previous)) delta%previous = self%previous
delta%f = self%f
! Oscillation equations
if (self%steps>=2) then ! self%previous should be used
dn = self%steps ; if (present(n)) dn = n
delta%state(1) = -self%f * self%previous(2, dn)
delta%state(2) = self%f * self%previous(1, dn)
else ! self%previous should not be used, use directly self%state
delta%state(1) = -self%f * self%state(2)
delta%state(2) = self%f * self%state(1)
endif
call move_alloc(delta, dState_dt)
return
!---------------------------------------------------------------------------------------------------------------------------------
endfunction dOscillation_dtHere delta is our residuals values that is finally copied into the class(integrand) :: dState_dt returning variable. As you can see, the residuals implementation has a very-high level syntax that is easy understand and maintain.
For multi-step (level) time solver (like the Adams-Bashforth class) it is important to provide a method for update the previous time steps solution each time the current solution is integrated over one time step.
Using the above oscillation type definition it simply corresponds to:
pure subroutine update_previous_steps(self)
!---------------------------------------------------------------------------------------------------------------------------------
!< Update previous time steps.
!---------------------------------------------------------------------------------------------------------------------------------
class(oscillation), intent(INOUT) :: self !< Oscillation field.
integer :: s !< Time steps counter.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
if (self%steps>0) then
do s=1, self%steps - 1
self%previous(:, s) = self%previous(:, s + 1)
enddo
self%previous(:, self%steps) = self%state
endif
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine update_previous_stepsIt a simple round-cycle that copy the solution of step n-s+1 into the n-s index, the solution of step n-s+2 into the n-s+1 index and so on until the current step n, s being the time steps used.
In the case you plan to not used multi-step solvers this method is not used, nevertheless it must be always implemented otherwise your Oscillation' definition is not a valid extension of the abstract type type_integrand. A simple implementation could be a do nothing method:
pure subroutine update_previous_steps(self)
!---------------------------------------------------------------------------------------------------------------------------------
!< FAKE update previous time steps for only non multi-step solvers.
!---------------------------------------------------------------------------------------------------------------------------------
class(oscillation), intent(INOUT) :: self !< Oscillation field.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine update_previous_stepsODE solvers typically involve operation between your integrand field (the Oscillation' system here) and numbers (in general reals) and other integrand field instances. This means that you must implement an almost complete set of operators methods for performing summation, multiplication, division, etc... Here we not provide more details on these methods, the interested readers can see the tests suite sources.
Let us now assume to integrate the above Oscillation' system by means of the Adams-Bashforth-Schemes class of solvers. The steps to this aim are very few:
...
real(R_P), parameter :: f=1e-4_R_P !< Frequency.
real(R_P), parameter :: initial_state(1:2)=[0., 1.] !< Initial state.
type(oscillation) :: attractor !< Oscillation field.
...
call attractor%init(initial_state=initial_state, f=f, steps=3)The Adams-Bashforth solvers class must be initialized selecting the time steps used:
...
type(adams_bashforth_integrator) :: ab_integrator !< Adams-Bashforth integrator.
...
call ab_integrator%init(steps=3)do while(.not.finished)
...
call ab_integrator%integrate(field=attractor, dt=0.01)
...
enddoFor a complete example see the Oscillation' test.