• individual elements are accessed by subscripting the array;
• default lower bound is 1 and can be omitted;
• upto 7 dimensions;
• in contrast to C/C++, Fortran arrays are column major.

                      integer, dimension(5, 3) :: array ! avoid magic numbers
integer                  :: c
integer                  :: r
do r=1, 5
  do c=1, 3
    print*, array(r, c) ! wrong memory access
  enddo
enddo

                    

                      integer, parameter                           :: rows_number=5
integer, parameter                           :: cols_number=3
integer, dimension(rows_number, cols_number) :: array
integer                                      :: c
integer                                      :: r
do c=1, cols_number
  do r=1, rows_number
    print*, array(r, c) ! right column-major memory access
  enddo
enddo

                    

Formatted

                      READ (eunit, format [, advance] [, asynchronous]
      [, blank] [, decimal] [, id] [, pad] [, pos]
      [, round] [, size] [, iostat] [, err] [, end]
      [, eor] [, iomsg]) [io-list]

                    

Formatted, List-Directed

                      READ (eunit, *[, asynchronous] [, blank]
      [, decimal] [, id] [, pad] [, pos]
      [, round] [, size] [, iostat] [, err] [, end]
      [, iomsg]) [io-list]

                    

Formatted, Namelist

                      READ (eunit, nml-group[, asynchronous] [, blank]
      [, decimal] [, id] [, pad] [, pos] [, round]
      [, size] [, iostat] [, err] [, end] [, iomsg])

                    

Unformatted

                      READ (eunit [, asynchronous] [, id] [, pos]
      [, iostat] [, err][, end] [, iomsg]) [io-list]

                    

Formatted

                      READ (eunit, format, rec [, asynchronous] [, blank]
      [, decimal] [, id] [, pad] [, pos] [, round]
      [, size] [, iostat] [, err] [, iomsg])
      [io-list]

                    

Unformatted

                      READ (eunit, rec [, asynchronous] [, id] [, pos]
      [, iostat] [, err] [, iomsg]) [io-list]

                    

                      READ (iunit, format [, nml-group] [, iostat] [, err]
      [, end] [, iomsg]) [io-list]

                    

Namelist

                      READ (iunit, nml-group [, iostat] [, err] [, end]
      [, iomsg]) [io-list]

                    

Formatted

                      WRITE (eunit, format [, advance] [, asynchronous]
       [, decimal] [, id] [, pos] [, round]
       [, sign] [, iostat] [ , err] [, [iomsg])
       [io-list]

                    

Formatted, List-Directed

                      WRITE (eunit, * [, asynchronous] [, decimal]
       [, delim] [, id] [, pos] [, round] [, sign]
       [, iostat] [, err] [, iomsg])[io-list]

                    

Formatted, Namelist

                      WRITE (eunit, nml-group [, asynchronous] [, decimal]
       [, delim] [, id] [, pos] [, round] [, sign]
       [, iostat] [, err] [, iomsg])

                    

Unformatted

                      WRITE (eunit [, asynchronous] [, id] [, pos]
       [, iostat] [, err] [, iomsg])[io-list]

                    

Formatted

                      WRITE (eunit, format, rec [, asynchronous]
       [, decimal] [, delim] [, id] [, pos]
       [, round] [, sign] [, iostat] [, err]
       [, iomsg])[io-list]

                    

Unformatted

                      WRITE (eunit, rec [, asynchronous] [, id] [, pos]
       [, iostat] [ , err] [, iomsg])[io-list]

                    

                      WRITE (iunit, format [, nml-group] [, iostat]
       [ , err] [, iomsg])[io-list]

                    

Namelist

                      WRITE (iunit, nml-group [, iostat] [ , err]
       [, iomsg])[io-list]

                    

                      OPEN ([UNIT=] io-unit[, FILE= name] [, ERR= label]
      [, IOMSG=msg-var] [, IOSTAT=i-var], slist)

                    

where slist is a list of one or more specifier, e.g. access='stream', form='unformatted' etc…

                      CLOSE ([UNIT=] io-unit[, STATUS = p] [, ERR= label]
       [, IOMSG=msg-var] [, IOSTAT=i-var])

                    

                      REWIND ([UNIT=] io-unit[, ERR= label]
        [, IOMSG=msg-var] [, IOSTAT=i-var])

                    

Inquiring by File

                      INQUIRE (FILE=name[, ERR=label] [, ID=id-var]
         [, IOMSG=msg-var] [, SIZE=sz]
         [, IOSTAT=i-var] slist)

                    

Inquiring by Unit

                      INQUIRE ([UNIT=]io-unit [, ERR=label] [, ID=id-var]
         [, IOMSG=msg-var] [, SIZE=sz]
         [, IOSTAT=i-var] slist)

                    

Inquiring by Output List

                      INQUIRE (IOLENGTH=len) out-item-list

                    

                      do i=1, N
  a(i) = cos(b(i-1)) + sin(c(i+1))
end do
do i=1, N
  d(i) = cos(e(i-1)) + sin(f(i+1))
end do
do i=1, N
  g(i) = cos(h(i-1)) + sin(p(i+1))
end do

                    

                      subroutine do_something(a, b, d)
...
end subroutine do_something

call do_something(a=a, b=b, c=c)
call do_something(a=d, b=e, c=f)
call do_something(a=g, b=h, c=i)

                    

• interface block is a powerful structure that was introduced in Fortran 90;
• give a calling procedure the full knowledge of the types and characteristics of the dummy arguments;
• provide a way to execute some safety checks when compiling the program;

• intent attribute was introduced in Fortran 90 and is recommended as it
  • allow compilers to check for coding errors;
  • facilitate efficient compilation and optimization;
• declare if a dummy argument is
  • Input: intent(in)
  • Output: intent(out)
  • Both: intent(inout)
• a dummy declared as intent(in) in a procedure cannot be changed during the execution of the procedure.

Old-fashioned Fortran library

                      subroutine foo(arg1, arg2, arg3)
real    :: arg1, agr2(*), arg3
integer :: i
! GOD is real...
god = sum(arg2)
i = bar(arg2)
arg3 = god/i
end subroutine foo

integer function bar(arg)
real :: arg(*)
bar = size(arg, dim=1)
end function bar

                    

Modern, robust Fortran module library

                      module library
implicit none
private
public :: foo
contains
  pure subroutine foo(arg1, arg2, arg3)
    real, intent(in)  :: arg1
    real, intent(in)  :: agr2(*)
    real, intent(out) :: arg3
    real              :: god
    integer           :: i
    god = sum(arg2)
    i = bar(arg2)
    arg3 = god/i
  end subroutine foo

  pure function bar(arg) result(bar_size)
    real, intent(in) :: arg(*)
    integer          :: bar_size
    bar_size = size(arg, dim=1)
  end function bar
end module library

                    

• once upon a time there was COMMON BLOCK… the evil!
• now Module can be used as a more robust GLOBALS CONTAINER

                      module globals
implicit none
private
public :: foo, bar
protected :: baz

integer :: foo = 0
real    :: bar = 0.0
logical :: baz = .false.

contains
  subroutine initialize
  foo = 1
  bar = 2.0
  baz = .true.
  end subroutine initialize
end module globals

                    

• all details of the data structure & order is isolated to a single source file. No need to replicate COMMON BLOCKs throughout each subroutine or use the non-Fortran 77 INCLUDE as a work-around.
• avoid problems with unequal COMMON BLOCK definitions.
• can use all data variables by default, or select and rename a subset.
• encapsulate data together methods!

• all the cons of Globals…
  • thread-un-safe: all use-associated users can (concurrently) modify public data;
  • passing globals by use-association make obscure the signature of procedures;

                      integer           :: Ni
integer           :: Nj
integer           :: Nk
real, allocatable :: centers(:,:,:)
real, allocatable :: nodes(:,:,:)
real, allocatable :: normals(:,:,:)
real, allocatable :: faces(:,:,:)
real, allocatable :: volumes(:,:,:)

call compute_metrics(Ni, Nj, Nk,    &
                     centers,nodes, &
                     normals, fases &
                     volumes)

                    

                        type :: mesh
  integer           :: Ni = 0
  integer           :: Nj = 0
  integer           :: Nk = 0
  real, allocatable :: centers(:,:,:)
  real, allocatable :: nodes(:,:,:)
  real, allocatable :: normals(:,:,:)
  real, allocatable :: faces(:,:,:)
  real, allocatable :: volumes(:,:,:)
end type mesh
type(mesh) :: a_mesh
call compute_metrics(grid=a_mesh)

                      

                      module globals
implicit none
integer           :: Ni
integer           :: Nj
integer           :: Nk
real, allocatable :: centers(:,:,:)
real, allocatable :: nodes(:,:,:)
real, allocatable :: normals(:,:,:)
real, allocatable :: faces(:,:,:)
real, allocatable :: volumes(:,:,:)
end module globals

program bad
use globals
implicit none
call compute_metrics(Ni, Nj, Nk,    &
                     centers,nodes, &
                     normals, fases &
                     volumes)
end program bad

                    

                        module mesh_t
implicit none
private
public :: mesh
type :: mesh
  integer           :: Ni = 0
  integer           :: Nj = 0
  integer           :: Nk = 0
  real, allocatable :: centers(:,:,:)
  real, allocatable :: nodes(:,:,:)
  real, allocatable :: normals(:,:,:)
  real, allocatable :: faces(:,:,:)
  real, allocatable :: volumes(:,:,:)
end type mesh
end module mesh_t
program good
use mesh_t, only: mesh
implicit none
type(mesh) :: a_mesh
call compute_metrics(grid=a_mesh)
end program good

                      

Generic Names avoid the plethora of different names exploiting dynamic dispatching at compile-time (no penalty!)

                      module generic
implicit none
private
public :: string
interface string
  module procedure string_integer, string_real
endinterface
contains
  elemental function string_integer(n) result(str)
  integer, intent(in) :: n
  character(11)       :: str
  write(str, ‘(I11)’) n
  end function string_integer
  elemental function string_real(n) result(str)
  real, intent(in) :: n
  character(13)    :: str
  write(str, ‘(E13.6E2)’) n
  end function string_real
end module generic

                    

Very CONCISE, CLEAR and ROBUST

                      use generic
integer :: step_number
real    :: time
real    :: delta
integer :: unit
integer :: t
delta = 0.1
time  = 0.0
do t=1, 100
  time = time + delta
  open(newunit=unit, &
    file='clock-'//string(t)//'.dat')
  write(unit, "(A)") &
    'current time: '//string(time)
  close(unit)
enddo

                    

A KISS library for exploiting codes portability for modern (2003+) Fortran projects

https://github.com/szaghi/PENF

                      module vector_t
implicit none
private
public :: vector, operator (+), operator(.cross.)
type, public :: Vector
  real :: x = 0.0 ! Cartesian component in x direction.
  real :: y = 0.0 ! Cartesian component in y direction.
  real :: z = 0.0 ! Cartesian component in z direction.
end type Vector
interface operator (+)
  module procedure add
end interface
interface operator (.cross.)
  module procedure crossproduct
end interface
contains
  elemental function add(left, rigth) result(summ)
  type(Vector), intent(in) :: left
  type(Vector), intent(in) :: rigth
  type(Vector)             :: summ
  summ%x = left%x + rigth%x
  summ%y = left%y + rigth%y
  summ%z = left%z + rigth%z
  end function add
  elemental function crossproduct(vec1, vec2) result(cross)
  type(Vector), intent(in) :: vec1
  type(Vector), intent(in) :: vec2
  type(Vector)             :: cross
  cross%x = (vec1%y * vec2%z) - (vec1%z * vec2%y)
  cross%y = (vec1%z * vec2%x) - (vec1%x * vec2%z)
  cross%z = (vec1%x * vec2%y) - (vec1%y * vec2%x)
  end function crossproduct
end module vector_t

                    

A KISS pure Fortran OOD class for computing Vectorial (3D) algebra

https://github.com/szaghi/VecFor

• restrict access to data (and data hiding):
  • more control;
  • easy to debug;
    • easy to maintain;
  • avoid side effects;
    • easy to parallel.

• design by-contract
  • modularize;
  • abstraction;
• develop by tests.

Define What is your Object and How it interacts

                      module gas_t
use vecfor, only : vector
implicit none
private
public :: gas ! Only the class GAS is public!
type :: gas
  private ! hide all members!
  ! the suffix '_' means internal implementation
  real         :: Cp_=0.0
  real         :: Cv_=0.0
  real         :: density_=0.0
  real         :: energy_=0.0
  type(vector) :: momentum_
  contains
    private
    ! the setter (constructor)
    procedure, public :: initialize
    ! return the members value (the getters)
    procedure, public :: density
    procedure, public :: momentum
    procedure, public :: energy
    ! interact with the rest of the World...
    procedure, public :: temperature
    procedure, public :: speed_of_sound
    ...
end type gas
contains
! follow the implementation
...
end module gas_t

                    

Define your World using the Contract

                      use gas_t, only : gas
implicit none

type(gas) :: air     ! a 'gas' object (instance)
real      :: density ! an alert dentisyt value

! respect the contract, use only public methods
call air%initialize(Cp=1004.8, &
                    Cv=716.0,  &
                    density=0.125)

! well designed API is auto-explicative...
if (air%speed()>air%speed_of_sound()) then
  ! do your supersonic stuffs...
end if

                    

So, where is encapsulation? The user of GAS class

• is not aware of how to compute the speed of sound..
  • he/she simply use the speed of sound;
• cannot manipulate gas data directly…
  • a bad user cannot make density negative manually;
• can re-use the same names defined into gas_t module:
  • encapsulation allow easy handling of names-collision (safety names space).

• re-use existing classes:
  • reduce code-duplication:
    • easy debug;
    • easy maintain;
  • inherit members/methods;
  • re-use implementation;
• classes hierarchy;
  • abstract is better, but specials happen…
  • specialize only when/where is necessary.

• design from down to UP…
  • from very down! Abstract all!
    • specialize => inherit from abstract to concrete cases
  • avoid to re-invent the wheel

Extend a Previous Contract adding new behaviors

                      module gas_multifluid_t
use gas_t, only : gas
implicit none
private
public :: gas_multifluid
! inherit from gas
type, extends(gas) :: gas_multifluid
  private
  integer           :: species_number=0
  real, allocatable :: concentrations(:)
  contains
    private
    ! overridded methods
    procedure, public :: initialize
    ...
end type gas_multifluid
contains
! follow the implementation
...
end module gas_multifluid_t

                    

Define your new World using the new Contract

                      use gas_t, only : gas
use gas_multifluid_t, only : gas_multifluid
implicit none
type(gas)            :: air
type(gas_multifluid) :: mixture
! respect the contract, use only public methods
call air%initialize(Cp=1004.8, &
                    Cv=716.0,  &
                    density=0.125)
call mixture%initialize(Cp=1004.8,            &
                        Cv=716.0,             &
                        species_number=3,     &
                        concentrations=[1.0,  &
                                        0.0,  &
                                        0.0], &
                        density=0.125)
! well designed API is auto-explicative...
if (air%speed()>air%speed_of_sound()) then
  ! do your supersonic stuffs...
end if

                    

So, where is inheritance? The developer of MULTIFLUID GAS class

• do not need to re-implement all the behaviors of gas class:
  • multifluid is a special gas… avoid-re-invent the wheel;

While the user of MULTIFLUID GAS class

• found the same facility of the previous gas class:
  • avoid to re-learn from scratch;
• could re-use his/her own applications designed for gas objects…

• design for family of classes (objects);
• exploit inheritance;
• exploit encapsulation.

• be ABSTARCT!

Develop your algorithms for family of classes… as you do for your daily-physics

                      module gas_fluxes_t
use gas_t, only : gas
implicit none
private
public :: fluxes
type :: fluxes
  real         :: mass_=0.0
  real         :: energy_=0.0
  type(vector) :: momentum_
  contains
 ...
   procedure, public :: compute
end type fluxes
contains
  subroutine compute(self, gas_cell)
  ! fluxes at interfaces to be computed
  class(fluxes), intent(inout) :: self
  ! gas into the cell from which compute fluxes
  class(gas),    intent(in)    :: gas_cell
  ...
  select type(gas_cell)
  class is (gas_multifluid)
  ...
  class is (gas)
  ...
  end select
  end subroutine compute
end module gas_fluxes_t

                    

                        use gas_t, only : gas
use gas_multifluid_t, only : gas_multifluid
use gas_library, only : compute_fluxes
implicit none
type(gas)            :: air
type(gas_multifluid) :: mixture
type(fluxes)         :: air_fluxes
type(fluxes)         :: mixture_fluxes
! respect the contract, use only public methods
call air%initialize(Cp=1004.8, &
                    Cv=716.0,  &
                    density=0.125)
call mixture%initialize(Cp=1004.8,            &
                        Cv=716.0,             &
                        species_number=3,     &
                        concentrations=[1.0,  &
                                        0.0,  &
                                        0.0], &
                        density=0.125)
! well designed API is auto-explicative…
if (air%speed()>air%speed_of_sound()) then
  call air_fluxes%compute(gas_cell=air)
  call mixture_fluxes%compute(gas_cell=mixture)
end if

                      

Develop an ODE solver for each system of equations:

• Oscillation Equations (ODE);
• Lorenz Equations (ODE);
• Burgers Equations (PDE);
• Euler Equations (PDE);
• Navier-Stokes Equations (PDE);
• …

One Ring to rule them all, One Ring to find them, One Ring to bring them all and in the darkness bind them

FOODIE, Fortran Object-Oriented Differential-equations Integration Environment

https://github.com/Fortran-FOSS-Programmers/FOODIE

Authors:

• Zaghi, S., CNR-INSEAN;
• Curcic, M., Dept. of Ocean Sciences Rosenstiel School of University of Miami;
• Rouson, D., Sourcery Institute and Sourcery, Inc.
• Beekman, I., Princeton/UMD CCROCCO LAB.

The ABSTRACT keyword implies

• you cannot instance a variable of type integrand, but only of its (concrete) extensions:
  • an abstract type is really a CONTRACT;
• you can (and want) develop for the abstract type integrand, you can ignore the implementation:
  • you must be aware (and use) of only its public interface!

The DEFERRED keyword implies

• you can (and want) ignore the implementation details of this method;
• you can (and want) rely on only it public interface;
• each (concrete) extension of the abstract type integrand must provide the implementation of each deferred method.

                      abstract interface
  function real_op_integrand(lhs, rhs) &
    result(operator_result)
  import :: integrand
  real,             intent(in)  :: lhs
  class(integrand), intent(in)  :: rhs
  class(integrand), allocatable :: operator_result
  end function real_op_integrand

  function time_derivative(self, t) &
    result(dState_dt)
  import :: integrand
  class(integrand),  intent(in) :: self
  real,    optional, intent(in) :: t
  class(integrand), allocatable :: dState_dt
  end function time_derivative

  function integrand_op_real(lhs, rhs) &
    result(operator_result)
  import :: integrand
  class(integrand), intent(in)  :: lhs
  real,             intent(in)  :: rhs
  class(integrand), allocatable :: operator_result
  end function integrand_op_real
end interface

                    

                        type, abstract :: integrand
  contains
    procedure(symmetric_operator),   deferred :: &
      integrand_multiply_integrand
    procedure(integrand_op_real),    deferred :: &
      integrand_multiply_real
    procedure(real_op_integrand),    deferred :: &
      real_multiply_integrand
    procedure(symmetric_operator),   deferred :: add
    …
end type integrand

                      

For all DEFERRED methods you must provide an

EXPLICIT ABSTRACT INTERFACE

                      abstract interface
  function symmetric_operator(lhs, rhs) &
    result(operator_result)
  import :: integrand
  class(integrand), intent(in)  :: lhs
  class(integrand), intent(in)  :: rhs
  class(integrand), allocatable :: operator_result
  endfunction symmetric_operator

  subroutine assignment_integrand(lhs, rhs)
  import :: integrand
  class(integrand), intent(inout) :: lhs
  class(integrand), intent(in)    :: rhs
  endsubroutine assignment_integrand
end interface

                    

                      module foodie_adt_integrand
implicit none
private
public :: integrand

type, abstract :: integrand
  contains
    procedure(time_derivative),      deferred :: t
    procedure(symmetric_operator),   deferred :: &
      integrand_multiply_integrand
    procedure(integrand_op_real),    deferred :: &
      integrand_multiply_real
    procedure(real_op_integrand),    deferred :: &
      real_multiply_integrand
    procedure(symmetric_operator),   deferred :: add
    procedure(symmetric_operator),   deferred :: sub
    procedure(assignment_integrand), deferred :: &
      assign_integrand
    generic :: operator(+) => add
    generic :: operator(-) => sub
    generic :: operator(*) => integrand_multiply_integrand,&
                              real_multiply_integrand, &
                              integrand_multiply_real
    generic :: assignment(=) => assign_integrand
endtype integrand

abstract interface
  function time_derivative(self, t) result(dState_dt)
  import :: integrand, R_P
  class(integrand),       intent(IN) :: self
  real(R_P),    optional, intent(IN) :: t
  class(integrand), allocatable      :: dState_dt
  end function time_derivative
  ...

                    

                          …
  function integrand_op_real(lhs, rhs) &
    result(operator_result)
  import :: integrand, R_P
  class(integrand), intent(IN)  :: lhs
  real(R_P),        intent(IN)  :: rhs
  class(integrand), allocatable :: operator_result
  end function integrand_op_real
  !
  function real_op_integrand(lhs, rhs) &
    result(operator_result)
  import :: integrand, R_P
  real(R_P),        intent(IN)  :: lhs
  class(integrand), intent(IN)  :: rhs
  class(integrand), allocatable :: operator_result
  end function real_op_integrand
  !
  function symmetric_operator(lhs, rhs) &
    result(operator_result)
  import :: integrand
  class(integrand), intent(IN)  :: lhs
  class(integrand), intent(IN)  :: rhs
  class(integrand), allocatable :: operator_result
  end function symmetric_operator
  !
  subroutine assignment_integrand(lhs, rhs)
  import :: integrand
  class(integrand), intent(INOUT) :: lhs
  class(integrand), intent(IN)    :: rhs
  end subroutine assignment_integrand
end interface
end module foodie_adt_integrand

                      

module foodie_integrator_euler_explicit
use foodie_adt_integrand, only : integrand
implicit none
private
public :: euler_explicit_integrator
!
type :: euler_explicit_integrator
  private
  contains
    private
    procedure, nopass, public :: integrate
end type euler_explicit_integrator
contains
  subroutine integrate(U, Dt, t)
  class(integrand),    intent(inout) :: U
  real(R_P),           intent(in)    :: Dt
  real(R_P), optional, intent(in)    :: t
  U = U + U%t(t=t) * Dt
  return
  end subroutine integrate
end module foodie_integrator_euler_explicit

Compare with the math $$ U(t+\Delta t) = U(t) + \Delta t U_t $$ It is almost the same…

• no knowldge of the concrete implementation on U is assumed;
• only the public API of U is used, e.g. U%t, U+U, U*Dt...;
• the euler integrator accepts ALL integrand extensions!

use foodie_integrator_euler_explicit
use my_integrand
type(euler_explicit_integrator) :: solver
type(my_integrand)              :: U
...
call solver%integrate(U=U, Dt=0.1)

module foodie_integrator_low_storage_runge_kutta
use foodie_adt_integrand, only : integrand
implicit none
private
public :: ls_runge_kutta_integrator
!
type :: ls_runge_kutta_integrator
  private
  integer           :: stages=0 !< Number of stages.
  real, allocatable :: A(:)     !< Low storage *A* coefficients.
  real, allocatable :: B(:)     !< Low storage *B* coefficients.
  real, allocatable :: C(:)     !< Low storage *C* coefficients.
  contains
    private
    ...
    procedure, pass(self), public :: integrate
endtype ls_runge_kutta_integrator
contains
  ...
  subroutine integrate(self, U, stage, Dt, t)
  class(ls_runge_kutta_integrator), intent(in)    :: self
  class(integrand),                 intent(inout) :: U
  class(integrand),                 intent(inout) :: stage(1:)
  real(R_P),                        intent(in)    :: Dt
  real(R_P),                        intent(in)    :: t
  integer(I_P)                                    :: s
  select type(stage)
  class is(integrand)
    stage(1) = U
    stage(2) = U*0._R_P
    do s=1, self%stages
      stage(2) = stage(2) * self%A(s) + &
        stage(1)%t(t=t + self%C(s) * Dt) * Dt
      stage(1) = stage(1) + stage(2) * self%B(s)
    enddo
    U = stage(1)
  endselect
  end subroutine integrate
end module foodie_integrator_low_storage_runge_kutta

Compare with the math $$ \begin{matrix} K_1 = U^n \\ K_2 = 0 \\ \left.\begin{matrix} K_2 = A_s K_2 + \Delta t R(t^n + C_s \Delta t, K_1) \\ K_1 = K_1 + B_s K_2 \end{matrix}\right\} s=1,2,...N_s\\ U^{n+1} = K_1 \end{matrix} $$ It is almost the same…

• no knowldge of the concrete implementation on U is assumed;
• only the public API of U is used, e.g. U%t, U+U, U*Dt...;
• the euler integrator accepts ALL integrand extensions!

use foodie_integrator_low_storage_runge_kutta
use my_integrand
type(ls_runge_kutta_integrator) :: solver
type(my_integrand)              :: U
...
call solver%integrate(U=U, Dt=0.1)

                      module oscillation_ode
use foodie, only : integrand
implicit none
private
public :: oscillation
type, extends(integrand) :: oscillation
  private
  real :: f=0.0
  real :: U(1:2)=[0.0, 0.0]
  contains
    procedure, pass(self), public :: t => &
      dOscillation_dt
    procedure, pass(lhs),  public :: &
      integrand_multiply_integrand => &
      oscillation_multiply_oscillation
    procedure, pass(lhs),  public :: &
      integrand_multiply_real => &
      oscillation_multiply_real
    procedure, pass(rhs),  public :: &
      real_multiply_integrand => &
      real_multiply_oscillation
    procedure, pass(lhs),  public :: add => &
      add_oscillation
    procedure, pass(lhs),  public :: sub => &
      sub_oscillation
    procedure, pass(lhs),  public :: &
      assign_integrand => &
      oscillation_assign_oscillation
end type oscillation
contains
  ...
end module oscillation_ode

                    

$$ \begin{matrix} U_t = R(U) \\ U = [v_1, v_2] R(U) = [-fv_2, fv_1] \end{matrix} $$

                        ...
  function dOscillation_dt(self, t) result(dState_dt)
  class(oscillation),     intent(in) :: self
  real(R_P),    optional, intent(in) :: t
  class(integrand),  allocatable     :: dState_dt
  allocate(oscillation :: dState_dt)
  select type(dState_dt)
  class is(oscillation)
    dState_dt = self
    dState_dt%U(1) = -self%f * self%U(2)
    dState_dt%U(2) =  self%f * self%U(1)
  endselect
  end function dOscillation_dt

  function oscillation_multiply_oscillation(lhs, rhs) &
    result(opr)
  class(oscillation), intent(in) :: lhs
  class(integrand),   intent(in) :: rhs
  class(integrand), allocatable  :: opr
  allocate(oscillation :: opr)
  select type(opr)
  class is(oscillation)
    opr = lhs
    select type(rhs)
    class is (oscillation)
      opr%U = lhs%U * rhs%U
    endselect
  endselect
  end function oscillation_multiply_oscillation

                    

• Pure Fortran implementation;
• KISS and user-friendly:
  • simple API, presently based on the Rouson’s Abstract Data Type Pattern;
  • easy building and porting on heterogeneous architectures;
• comprehensive solvers set out-of-the-box:
  • explicit schemes:
    • Adams-Bashforth schemes:
      • 1 step, namely the forward explicit Euler scheme, 1st order accurate;
      • 2 to 16 steps, 2nd to 16th accurate, respectively;
    • Euler (forward explicit) scheme, 1st order accurate;
    • Leapfrog, 2nd order accurate:
      • unfiltered leapfrog, 2nd order accurate, mostly unstable;
      • Robert-Asselin filtered leapfrog, 1st order accurate;
      • Robert-Asselin-Williams filtered leapfrog, 3rd order accurate;
    • Runge-Kutta schemes:
      • low-storage schemes:
        • 1 stage, namely the forward explicit Euler scheme, 1st order accurate;
        • 5 stages, 4th order accurate, 2N registers;
        • 6 stages, 4th order accurate, 2N registers;
        • 7 stages, 4th order accurate, 2N registers;
        • 12 stages, 4th order accurate, 2N registers;
        • 14 stages, 4th order accurate, 2N registers;

• comprehensive solvers set out-of-the-box:
  • explicit schemes:
    • Runge-Kutta schemes:
      • TVD/SSP schemes:
        • 1 stage, namely the forward explicit Euler scheme, 1st order accurate;
        • 2 stages, 2nd order accurate;
        • 3 stages, 3rd order accurate;
        • 5 stages, 4th order accurate;
      • embedded (adaptive) schemes:
        • Heun-Euler, 2 stages, 2nd order accurate;
        • Runge-Kutta-Cash-Karp, 6 stages, 5th order accurate;
        • Prince-Dormand, 7 stages, 4th order accurate;
        • Calvo, 9 stages, 6th order accurate;
        • Feagin, 17 stages, 10th order accurate;
  • implicit schemes:
    • Adams-Moulton schemes:
      • 0 step, 1st order accurate;
      • 1 step, 2nd accurate;
      • 2 steps, 3rd accurate;
      • 3 steps, 4th accurate;
    • Backward Differentiation Formula schemes:
      • 1 step, namely the backward implicit Euler scheme, 1st order accurate;
      • 2 to 6 steps, 2nd to 6th accurate, respectively;
  • predictor-corrector schemes:
    • Adams-Bashforth-Moulton schemes:
      • 1 step, AB(1)-AM(0), 1st order accurate;
      • 2 steps, AB(2)-AM(1), 2nd accurate;
      • 3 steps, AB(3)-AM(2), 3rd accurate;
      • 4 steps, AB(4)-AM(3), 4th accurate;