system.F90

Go to the documentation of this file.

!! Copyright (C) 2019 N. Tancogne-Dejean
!! Copyright (C) 2020 M. Oliveira, Heiko Appel
!! Copyright (C) 2021 S. Ohlmann
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module system_oct_m
  use algorithm_oct_m
  use algorithm_factory_oct_m
  use debug_oct_m
  use ghost_interaction_oct_m
  use global_oct_m
  use interactions_factory_abst_oct_m
  use interaction_partner_oct_m
  use interaction_oct_m
  use iteration_counter_oct_m
  use messages_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use multisystem_debug_oct_m
  use linked_list_oct_m
  use parser_oct_m
  use profiling_oct_m
  use quantity_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use varinfo_oct_m
  implicit none
 
  private
  public ::                                    &
    system_t,                                  &
    system_execute_algorithm,                  &
    system_init_parallelization,               &
    system_new_algorithm,                      &
    system_init_iteration_counters,            &
    system_reset_iteration_counters,           &
    system_create_interactions,                &
    system_algorithm_start,                    &
    system_algorithm_finish,                   &
    system_restart_read,                       &
    system_restart_write,                      &
    system_update_potential_energy,            &
    system_update_total_energy,                &
    system_end,                                &
    system_list_t,                             &
    system_iterator_t
 
  type :: barrier_t
    logical :: active
    real(real64) :: target_time
  end type barrier_t
 
  integer, parameter, public :: &
    NUMBER_BARRIERS = 1,        &
    barrier_restart = 1
 
  type, extends(interaction_partner_t), abstract :: system_t
    private
    type(iteration_counter_t),    public :: iteration
    class(algorithm_t),  pointer, public :: algo => null()
 
    integer, allocatable, public :: supported_interactions(:)
    type(interaction_list_t), public :: interactions
 
    type(mpi_grp_t), public :: grp
 
    type(barrier_t) :: barrier(NUMBER_BARRIERS)
    real(real64), public :: kinetic_energy
    real(real64), public :: potential_energy
    real(real64), public :: internal_energy
    !                                          !! the kinetic energy of its constituents
    real(real64), public :: total_energy
 
  contains
    procedure :: execute_algorithm =>  system_execute_algorithm 
    procedure :: reset_iteration_counters => system_reset_iteration_counters 
    procedure :: new_algorithm => system_new_algorithm 
    procedure :: algorithm_finished => system_algorithm_finished 
    procedure :: init_iteration_counters => system_init_iteration_counters 
    procedure :: create_interactions => system_create_interactions 
    procedure :: init_parallelization => system_init_parallelization 
    procedure :: update_couplings => system_update_couplings  
    procedure :: update_interactions => system_update_interactions 
    procedure :: update_interactions_start => system_update_interactions_start 
    procedure :: update_interactions_finish => system_update_interactions_finish 
    procedure :: algorithm_start => system_algorithm_start 
    procedure :: algorithm_finish => system_algorithm_finish 
    procedure :: iteration_info => system_iteration_info 
    procedure :: restart_write => system_restart_write 
    procedure :: restart_read => system_restart_read 
    procedure :: output_start => system_output_start 
    procedure :: output_write => system_output_write 
    procedure :: output_finish => system_output_finish 
    procedure :: process_is_slave => system_process_is_slave 
    procedure :: start_barrier => system_start_barrier 
    procedure :: end_barrier => system_end_barrier 
    procedure :: arrived_at_barrier => system_arrived_at_barrier 
    procedure :: arrived_at_any_barrier => system_arrived_at_any_barrier 
    procedure :: update_potential_energy => system_update_potential_energy 
    procedure :: update_internal_energy => system_update_internal_energy 
    procedure :: update_total_energy => system_update_total_energy 
    procedure(system_init_interaction),          deferred :: init_interaction
    procedure(system_initialize),                deferred :: initialize
    procedure(system_do_algorithmic_operation),  deferred :: do_algorithmic_operation
    procedure(system_is_tolerance_reached),      deferred :: is_tolerance_reached
    procedure(system_restart_write_data),        deferred :: restart_write_data
    procedure(system_restart_read_data),         deferred :: restart_read_data
    procedure(system_update_kinetic_energy),     deferred :: update_kinetic_energy
  end type system_t
 
  abstract interface
 
    ! ---------------------------------------------------------
    subroutine system_init_interaction(this, interaction)
      import system_t
      import interaction_t
      class(system_t), target, intent(inout) :: this
      class(interaction_t), intent(inout) :: interaction
    end subroutine system_init_interaction
 
    ! ---------------------------------------------------------
    subroutine system_initialize(this)
      import system_t
      class(system_t), intent(inout) :: this
    end subroutine system_initialize
 
    ! ---------------------------------------------------------
    logical function system_do_algorithmic_operation(this, operation, updated_quantities) result(done)
      import system_t
      import algorithmic_operation_t
      class(system_t),                 intent(inout) :: this
      class(algorithmic_operation_t),  intent(in)    :: operation
      character(len=:),   allocatable, intent(out)   :: updated_quantities(:)
    end function system_do_algorithmic_operation
 
    ! ---------------------------------------------------------
    logical function system_is_tolerance_reached(this, tol)
      use, intrinsic :: iso_fortran_env
      import system_t
      class(system_t), intent(in) :: this
      real(real64),    intent(in) :: tol
    end function system_is_tolerance_reached
 
    ! ---------------------------------------------------------
    subroutine system_store_current_status(this)
      import system_t
      class(system_t), intent(inout) :: this
    end subroutine system_store_current_status
 
    ! ---------------------------------------------------------
    subroutine system_restart_write_data(this)
      import system_t
      class(system_t), intent(inout) :: this
    end subroutine system_restart_write_data
 
    ! ---------------------------------------------------------
    ! this function returns true if restart data could be read
    logical function system_restart_read_data(this)
      import system_t
      class(system_t), intent(inout) :: this
    end function system_restart_read_data
    subroutine system_update_kinetic_energy(this)
      import system_t
      class(system_t),      intent(inout) :: this
    end subroutine system_update_kinetic_energy
 
  end interface
 
  type, extends(partner_list_t) :: system_list_t
    private
  contains
    procedure :: add => system_list_add_node 
    procedure :: contains => system_list_contains 
  end type system_list_t
 
  type, extends(linked_list_iterator_t) :: system_iterator_t
    private
  contains
    procedure :: get_next => system_iterator_get_next 
  end type system_iterator_t
 
contains
 
  ! ---------------------------------------------------------
  subroutine system_execute_algorithm(this)
    class(system_t),     intent(inout) :: this
 
    type(algorithmic_operation_t) :: operation
    logical :: all_updated, at_barrier, operation_done
    type(event_handle_t) :: debug_handle
    integer :: i
    type(quantity_t), pointer :: quantity
    character(len=:), allocatable :: updated_quantities(:)
 
    push_sub(system_execute_algorithm)
 
    at_barrier = .false.
 
    do while (.not. at_barrier)
 
      operation = this%algo%get_current_operation()
 
      debug_handle = multisystem_debug_write_event_in(this%namespace, event_function_call_t("dt_operation", operation),    &
        system_iteration=this%iteration, algo_iteration=this%algo%iteration)
 
      ! First try to execute the operation as a system specific operation
      operation_done = this%do_algorithmic_operation(operation, updated_quantities)
      if (allocated(updated_quantities)) then
        ! Update the quantities iteration counters
        do i = 1, size(updated_quantities)
          quantity => this%quantities%get(updated_quantities(i))
          call quantity%iteration%set(this%algo%iteration + 1)
          call multisystem_debug_write_marker(this%namespace, event_iteration_update_t("quantity", quantity%label, &
            quantity%iteration, "set"))
        end do
      end if
 
      ! If not done, we try to execute it as an algorithm-specific operation.
      if (.not. operation_done) then
        operation_done = this%algo%do_operation(operation)
      else
        call this%algo%next()
      end if
 
      ! If still not done, the operation must be a generic operation
      if (.not. operation_done) then
 
        select case (operation%id)
        case (skip)
          ! Do nothing
          call this%algo%next()
 
        case (iteration_done)
          ! Increment the system iteration by one
          this%iteration = this%iteration + 1
          call multisystem_debug_write_marker(this%namespace, event_iteration_update_t("system",  "", this%iteration, "tick"))
 
          ! Recompute the total energy
          call this%update_total_energy()
 
          ! Write output
          call this%output_write()
 
          ! Update elapsed time
          call this%algo%update_elapsed_time()
 
          ! Print information about the current iteration
          ! (NB: needs to be done after marking the execution step as finished,
          ! so that the timings are correct)
          call this%iteration_info()
 
          call this%algo%next()
 
        case (rewind_algorithm)
          if (.not. this%algo%finished()) then
            if (.not. this%arrived_at_any_barrier()) then
              ! Reset propagator for next step if not waiting at barrier
              call this%algo%rewind()
            else
              at_barrier = .true.
            end if
          else
            if (this%algo%continues_after_finished()) then
              call this%algo%rewind()
            else
              at_barrier = .true.
            end if
          end if
 
        case (update_couplings)
          ! We increment by one algorithmic step
          this%algo%iteration = this%algo%iteration + 1
          call multisystem_debug_write_marker(this%namespace, event_iteration_update_t("algorithm", "", &
            this%algo%iteration, "tick"))
 
          ! Try to update all needed quantities from the interaction partners
          all_updated = this%update_couplings()
 
          ! Move to next algorithm step if all couplings have been
          ! updated. Otherwise roll back the iteration counter and try again later.
          if (all_updated) then
            call this%algo%next()
          else
            this%algo%iteration = this%algo%iteration - 1
            call multisystem_debug_write_marker(this%namespace, event_iteration_update_t("algorithm", "", &
              this%algo%iteration, "reverse"))
          end if
 
          ! Couplings are implicit barriers
          at_barrier = .true.
 
        case (update_interactions)
          ! Update all the interactions
          call this%update_interactions()
          call this%algo%next()
 
        case default
          message(1) = "Unsupported algorithmic operation."
          write(message(2), '(A,A,A)') trim(operation%id), ": ", trim(operation%label)
          call messages_fatal(2, namespace=this%namespace)
        end select
      end if
 
      call multisystem_debug_write_event_out(debug_handle, system_iteration=this%iteration, algo_iteration=this%algo%iteration)
    end do
 
    pop_sub(system_execute_algorithm)
  end subroutine system_execute_algorithm
 
  ! ---------------------------------------------------------
  subroutine system_reset_iteration_counters(this, accumulated_iterations)
    class(system_t),      intent(inout) :: this
    integer,              intent(in)    :: accumulated_iterations
 
    type(interaction_iterator_t) :: iter
    class(interaction_t), pointer :: interaction
    type(quantity_iterator_t) :: qiter
    class(quantity_t), pointer :: quantity
 
    character(len=MAX_INFO_LEN) :: extended_label
 
    push_sub(system_reset_iteration_counters)
 
    ! Propagator counter
    this%algo%iteration = this%algo%iteration - accumulated_iterations
    call multisystem_debug_write_marker(this%namespace, event_iteration_update_t("algorithm", "", this%algo%iteration, "reset"))
 
    ! Interaction counters
    call iter%start(this%interactions)
    do while (iter%has_next())
      interaction => iter%get_next()
      interaction%iteration = interaction%iteration - accumulated_iterations
 
      extended_label = trim(interaction%label)//"-"//trim(interaction%partner%namespace%get())
      call multisystem_debug_write_marker(this%namespace, event_iteration_update_t( "interaction", extended_label, &
        interaction%iteration, "reset"))
    end do
 
    ! Internal quantities counters
    call qiter%start(this%interactions)
    do while(qiter%has_next())
      quantity => qiter%get_next()
      quantity%iteration = quantity%iteration - accumulated_iterations
      call multisystem_debug_write_marker(this%namespace, event_iteration_update_t("quantity", quantity%label, &
        quantity%iteration, "reset"))
    end do
 
    pop_sub(system_reset_iteration_counters)
  end subroutine system_reset_iteration_counters
 
 
  ! ---------------------------------------------------------
  recursive subroutine system_create_interactions(this, interaction_factory, available_partners)
    class(system_t),                       intent(inout) :: this
    class(interactions_factory_abst_t),    intent(in)    :: interaction_factory
    class(partner_list_t),         target, intent(in)    :: available_partners
    !                                                                           !! for the given system.
 
    logical :: in_list
    integer :: i, ip, iq, interaction_type
    type(partner_list_t) :: partners
    type(partner_iterator_t) :: iter
    class(interaction_partner_t), pointer :: partner
    class(interaction_t), pointer :: interaction
    type(interactions_factory_options_t), allocatable :: options(:)
 
    push_sub(system_create_interactions)
 
    ! Sanity checks
    assert(allocated(this%supported_interactions))
    assert(allocated(this%supported_interactions_as_partner))
 
 
    ! All systems need to be connected to make sure they remain synchronized.
    ! We enforce that be adding a ghost interaction between all systems.
    ! Note that this recursively adds subsystems of containers.
    call iter%start(available_partners)
    do while (iter%has_next())
      partner => iter%get_next()
      call partner%add_partners_to_list(partners)
    end do
 
    call iter%start(partners)
    do while (iter%has_next())
      partner => iter%get_next()
 
      ! No self-interaction
      if (partner%namespace%get() == this%namespace%get()) cycle
 
      call this%interactions%add(ghost_interaction_t(partner))
    end do
 
 
    ! Now create the non-ghost interactions:
    ! Here we only add systems to the parters list, which support the given interaction as a partner
    ! Note that ensemble containers do not add their subsystems to the list.
    ! (see ensemble_oct_m::ensemble_add_parters_to_list())
 
    call partners%empty()
 
    ! Get options to use to create all the interactions supported by this system
    allocate(options(0)) ! Workaround to avoid a gfortran warning
    options = interaction_factory%options(this%namespace, this%supported_interactions)
 
    ! Loop over all interactions supported by this system and create them according to the provided options
    do i = 1, size(this%supported_interactions)
      interaction_type = this%supported_interactions(i)
 
      ! Supported interactions should only appear once in the corresponding
      ! list, otherwise more than one interaction of the same type will be
      ! created
      assert(count(this%supported_interactions == interaction_type) == 1)
 
      ! Get the list of partners that support this interaction type
      call iter%start(available_partners)
      do while (iter%has_next())
        partner => iter%get_next()
        call partner%add_partners_to_list(partners, interaction_type)
      end do
 
      ! Create list of partners for this interaction taking into account the selected mode
      select case (options(i)%mode)
      case (all_partners)
        ! Use all available partners, so no need to modify the list
 
      case (no_partners)
        ! No partners for this interaction, so empty the list
        call partners%empty()
 
      case (only_partners)
        ! Start with full list and remove the partners not listed in the options
        call iter%start(partners)
        do while (iter%has_next())
          partner => iter%get_next()
          in_list = .false.
          do ip = 1, size(options(i)%partners)
            if (partner%namespace%is_contained_in(options(i)%partners(ip))) then
              in_list = .true.
              exit
            end if
          end do
          if (.not. in_list) then
            call partners%delete(partner)
          end if
        end do
 
      case (all_except)
        ! Start with full list and remove the partners listed in the options
        do ip = 1, size(options(i)%partners)
          call iter%start(partners)
          do while (iter%has_next())
            partner => iter%get_next()
            if (partner%namespace%is_contained_in(options(i)%partners(ip))) then
              call partners%delete(partner)
            end if
          end do
        end do
 
      end select
 
      ! Now actually create the interactions for the selected partners
      call iter%start(partners)
      do while (iter%has_next())
        partner => iter%get_next()
 
        interaction => interaction_factory%create(interaction_type, partner)
 
 
        ! Set flag for self-interaction
        interaction%intra_interaction = partner%namespace%get() == this%namespace%get()
 
        ! Set timing and perform sanity checks
        interaction%timing = options(i)%timing
        if (interaction%timing == timing_exact) then
          select type (partner => interaction%partner)
          class is (system_t)
            if (this%algo%iteration%global_step() /= partner%algo%iteration%global_step() .and. &
              .not. all(partner%quantities%always_available(interaction%couplings_from_partner))) then
              write(message(1), '(2a)') "InteractionTiming was set to exact timing, but systems ", &
                trim(this%namespace%get())
              write(message(2), '(3a)') "and ", trim(partner%namespace%get()), " have incompatible steps."
              call messages_fatal(2, namespace=this%namespace)
            end if
          end select
        end if
 
        ! Initialize interaction from system and from partner
        call this%init_interaction(interaction)
        call interaction%partner%init_interaction_as_partner(interaction)
 
        ! Make sure all quantities required by the interaction are available in the system and partner
        if (allocated(interaction%system_quantities)) then
          do iq = 1, size(interaction%system_quantities)
            if (.not. associated(this%quantities%get(interaction%system_quantities(iq)))) then
              write(message(1), '(5a)') "Interaction '", trim(interaction%label), "' requires quantity '", &
                trim(interaction%system_quantities(iq)), "'"
              write(message(2), '(3a)') "from system '", trim(this%namespace%get()),  "' but it is not available."
              call messages_fatal(2, namespace=this%namespace)
            end if
          end do
        end if
        if (allocated(interaction%couplings_from_partner)) then
          do iq = 1, size(interaction%couplings_from_partner)
            if (.not. associated(partner%quantities%get(interaction%couplings_from_partner(iq)))) then
              write(message(1), '(5a)') "Interaction '", trim(interaction%label), "' requires coupling '", &
                trim(interaction%couplings_from_partner(iq)), "'"
              write(message(2), '(3a)') "from partner '", trim(partner%namespace%get()),  "' but it is not available."
              call messages_fatal(2, namespace=this%namespace)
            end if
          end do
        end if
 
        ! Add interaction to list
        call this%interactions%add(interaction)
      end do
 
      ! Empty partner list for next interaction
      call partners%empty()
    end do
 
    pop_sub(system_create_interactions)
  end subroutine system_create_interactions
 
  ! ---------------------------------------------------------
  logical function system_update_couplings(this) result(all_updated)
    class(system_t),      intent(inout) :: this
 
    class(interaction_t),              pointer :: interaction
    type(interaction_iterator_t) :: iter
 
    push_sub(system_update_couplings)
 
    all_updated = .true.
    call iter%start(this%interactions)
    do while (iter%has_next())
      interaction => iter%get_next()
 
      select type (partner => interaction%partner)
      class is (system_t)
        ! If the partner is a system, we can only update its couplings if it
        ! is not too much behind the requested iteration
        if (partner%algo%iteration + 1 >= this%algo%iteration) then
          call interaction%update_partner_couplings(this%algo%iteration)
        end if
 
      class default
        ! Partner that are not systems can have their couplings updated at any iteration
        call interaction%update_partner_couplings(this%algo%iteration)
      end select
 
      all_updated = all_updated .and. interaction%partner_couplings_up_to_date
    end do
 
    pop_sub(system_update_couplings)
  end function system_update_couplings
 
  ! ---------------------------------------------------------
  subroutine system_update_interactions(this)
    class(system_t),      intent(inout) :: this
 
    integer :: iq, n_quantities
    class(interaction_t), pointer :: interaction
    type(interaction_iterator_t)  :: iter
    class(quantity_t), pointer :: quantity
 
    push_sub(system_update_interactions)
 
    ! Some systems might need to perform some specific operations before the
    ! update.
    call this%update_interactions_start()
 
    !Loop over all interactions
    call iter%start(this%interactions)
    do while (iter%has_next())
      interaction => iter%get_next()
 
      ! Update the system quantities that will be needed for computing the interaction
      if (allocated(interaction%system_quantities)) then
        n_quantities = size(interaction%system_quantities)
      else
        n_quantities = 0
      end if
      do iq = 1, n_quantities
        ! Get requested quantity
        quantity => this%quantities%get(interaction%system_quantities(iq))
 
        if (.not. quantity%iteration == this%algo%iteration) then
          ! The requested quantity is not at the requested iteration
 
          ! Sanity check: it should never happen that the quantity is in advance
          ! with respect to the requested iteration.
          if (quantity%iteration > this%algo%iteration) then
            message(1) = "The quantity "//trim(quantity%label)//" is in advance compared to the requested iteration."
            message(2) = "The interaction is "//trim(interaction%label)//"."
            call messages_fatal(2, namespace=this%namespace)
          end if
 
          ! If the quantity can be updated on demand, we try to update it
          if (quantity%updated_on_demand) then
            call this%update_on_demand_quantity(quantity, this%algo%iteration)
          end if
 
          ! Quantity must be at the correct iteration, otherwise the algorithm and
          ! the interaction are incompatible
          if (.not. quantity%iteration == this%algo%iteration .and. .not. quantity%always_available) then
            write(message(1), '(5a)') "Interaction ", trim(interaction%label), " is incompatible with the selected algorithm."
            write(message(2), '(3a)') "The interaction requires the ", trim(quantity%label), &
              " at an iteration it is not available."
            call messages_fatal(2, namespace=this%namespace)
          end if
 
        end if
 
      end do
 
      call interaction%update(this%algo%iteration)
    end do
 
    ! Some systems might need to perform some specific operations after all the
    ! interactions have been updated
    call this%update_interactions_finish()
 
    pop_sub(system_update_interactions)
  end subroutine system_update_interactions
 
  ! ---------------------------------------------------------
  subroutine system_update_interactions_start(this)
    class(system_t), intent(inout) :: this
 
    push_sub(system_update_interactions_start)
 
    ! By default nothing is done just before updating the interactions. Child
    ! classes that wish to change this behaviour should override this method.
 
    pop_sub(system_update_interactions_start)
  end subroutine system_update_interactions_start
 
  ! ---------------------------------------------------------
  subroutine system_update_interactions_finish(this)
    class(system_t), intent(inout) :: this
 
    push_sub(system_update_interactions_finish)
 
    ! By default nothing is done just after updating the interactions. Child
    ! classes that wish to change this behaviour should override this method.
 
    pop_sub(system_update_interactions_finish)
  end subroutine system_update_interactions_finish
 
  ! ---------------------------------------------------------
  subroutine system_restart_write(this)
    class(system_t), intent(inout) :: this
 
    logical :: restart_write
    type(interaction_iterator_t) :: iter
    class(interaction_t),              pointer :: interaction
    type(quantity_iterator_t) :: qiter
    class(quantity_t), pointer :: quantity
 
 
    push_sub(system_restart_write)
 
    call parse_variable(this%namespace, 'RestartWrite', .true., restart_write)
 
    if (restart_write) then
      ! do some generic restart steps here
      ! write iteration counter restart data
      call this%iteration%restart_write('restart_iteration_system', this%namespace)
      call this%algo%iteration%restart_write('restart_iteration_algorithm', this%namespace)
      call iter%start(this%interactions)
      do while (iter%has_next())
        interaction => iter%get_next()
        call interaction%restart_write(this%namespace)
      end do
      call qiter%start(this%quantities)
      do while (qiter%has_next())
        quantity => qiter%get_next()
        call quantity%iteration%restart_write('restart_iteration_quantity_'//trim(quantity%label), &
          this%namespace)
      end do
      ! the following call is delegated to the corresponding system
      call this%restart_write_data()
      message(1) = "Wrote restart data for system "//trim(this%namespace%get())
      call messages_info(1, namespace=this%namespace)
    end if
 
    pop_sub(system_restart_write)
  end subroutine system_restart_write
 
  ! ---------------------------------------------------------
  ! this function returns true if restart data could be read
  logical function system_restart_read(this)
    class(system_t), intent(inout) :: this
 
    type(interaction_iterator_t) :: iter
    class(interaction_t),              pointer :: interaction
    type(quantity_iterator_t) :: qiter
    class(quantity_t), pointer :: quantity
 
    push_sub(system_restart_read)
 
    ! do some generic restart steps here
    ! read iteration data
    system_restart_read = this%iteration%restart_read('restart_iteration_system', this%namespace)
    system_restart_read = system_restart_read .and. &
      this%algo%iteration%restart_read('restart_iteration_algorithm', this%namespace)
    call iter%start(this%interactions)
    do while (iter%has_next())
      interaction => iter%get_next()
      system_restart_read = system_restart_read .and. interaction%restart_read(this%namespace)
      ! reduce by one because of the first UPDATE_INTERACTIONS
      interaction%iteration = interaction%iteration - 1
    end do
    call qiter%start(this%quantities)
    do while (qiter%has_next())
      quantity => qiter%get_next()
      system_restart_read = system_restart_read .and. &
        quantity%iteration%restart_read('restart_iteration_quantity_'//trim(quantity%label), &
        this%namespace)
      if (quantity%updated_on_demand) then
        ! decrease iteration by one for on-demand quantities to account for initial update_interactions step
        quantity%iteration = quantity%iteration - 1
      end if
    end do
    ! the following call is delegated to the corresponding system
    system_restart_read = system_restart_read .and. this%restart_read_data()
 
    if (system_restart_read) then
      message(1) = "Successfully read restart data for system "//trim(this%namespace%get())
      call messages_info(1, namespace=this%namespace)
    end if
 
    pop_sub(system_restart_read)
  end function system_restart_read
 
  ! ---------------------------------------------------------
  subroutine system_output_start(this)
    class(system_t), intent(inout) :: this
 
    push_sub(system_output_start)
 
    ! By default nothing is done to regarding output. Child classes that wish
    ! to change this behaviour should override this method.
 
    pop_sub(system_output_start)
  end subroutine system_output_start
 
  ! ---------------------------------------------------------
  subroutine system_output_write(this)
    class(system_t), intent(inout) :: this
 
    push_sub(system_output_write)
 
    ! By default nothing is done to regarding output. Child classes that wish
    ! to change this behaviour should override this method.
 
    pop_sub(system_output_write)
  end subroutine system_output_write
 
  ! ---------------------------------------------------------
  subroutine system_output_finish(this)
    class(system_t), intent(inout) :: this
 
    push_sub(system_output_finish)
 
    ! By default nothing is done to regarding output. Child classes that wish
    ! to change this behaviour should override this method.
 
    pop_sub(system_output_finish)
  end subroutine system_output_finish
 
  ! ---------------------------------------------------------
  subroutine system_new_algorithm(this, factory)
    class(system_t),            intent(inout) :: this
    class(algorithm_factory_t), intent(in)    :: factory
 
    integer :: ii
 
    push_sub(system_new_algorithm)
 
    call messages_experimental('Multi-system framework')
 
    this%algo => factory%create(this)
 
    call this%init_iteration_counters()
 
    do ii = 1, number_barriers
      this%barrier(ii)%active = .false.
      this%barrier(ii)%target_time = m_zero
    end do
 
    pop_sub(system_new_algorithm)
  end subroutine system_new_algorithm
 
  ! ---------------------------------------------------------------------------------------
  recursive function system_algorithm_finished(this) result(finished)
    class(system_t),       intent(in) :: this
    logical :: finished
 
    finished = this%algo%finished()
 
  end function system_algorithm_finished
 
  ! ---------------------------------------------------------
  !
  subroutine system_init_iteration_counters(this)
    class(system_t),            intent(inout) :: this
 
    type(interaction_iterator_t) :: iter
    class(interaction_t),              pointer :: interaction
    type(quantity_iterator_t) :: qiter
    class(quantity_t), pointer :: quantity
 
    push_sub(system_init_iteration_counters)
 
    ! Initialize algorithm and system counters
    call this%algo%init_iteration_counters()
 
    ! Interactions counters
    call iter%start(this%interactions)
    do while (iter%has_next())
      interaction => iter%get_next()
      interaction%iteration = this%algo%iteration - 1
    end do
 
    ! Quantities counters
    call qiter%start(this%quantities)
    do while (qiter%has_next())
      quantity => qiter%get_next()
      if (quantity%updated_on_demand) then
        ! On demand quantities will be updated before first use,
        ! hence we set the iteration counter one iteration behind the algorithms iteration counter
        quantity%iteration = this%algo%iteration - 1
      else
        quantity%iteration = this%algo%iteration
      end if
    end do
 
    pop_sub(system_init_iteration_counters)
  end subroutine system_init_iteration_counters
 
  ! ---------------------------------------------------------
  subroutine system_algorithm_start(this)
    class(system_t),      intent(inout) :: this
 
    logical :: all_updated
    type(event_handle_t) :: debug_handle
    character(len=:), allocatable :: updated_quantities(:)
 
    push_sub(system_algorithm_start)
 
    debug_handle = multisystem_debug_write_event_in(this%namespace, event_function_call_t("system_algorithm_start"), &
      system_iteration = this%iteration, algo_iteration = this%algo%iteration)
 
    ! Update interactions at initial iteration
    all_updated = this%update_couplings()
    if (.not. all_updated) then
      message(1) = "Unable to update interactions when initializing the algorithm."
      call messages_fatal(1, namespace=this%namespace)
    end if
    call this%update_interactions()
 
    ! System-specific and algorithm-specific initialization step
    if (this%algo%start_operation%id /= skip) then
      if (.not. this%do_algorithmic_operation(this%algo%start_operation, updated_quantities)) then
        message(1) = "Unsupported algorithmic operation."
        write(message(2), '(A,A,A)') trim(this%algo%start_operation%id), ": ", trim(this%algo%start_operation%label)
        call messages_fatal(2, namespace=this%namespace)
      end if
      if (allocated(updated_quantities)) then
        message(1) = "Update of quantities not allowed in algorithmic operation."
        write(message(2), '(A,A,A)') trim(this%algo%start_operation%id), ": ", trim(this%algo%start_operation%label)
        call messages_fatal(2, namespace=this%namespace)
      end if
    end if
 
    ! Compute the total energy at the beginning of the simulation
    call this%update_total_energy()
 
    ! Start output
    call this%output_start()
 
    ! Write header for execution log
    call this%algo%write_output_header()
 
    ! Rewind algorithm (will also set the initial time to compute the elapsed time)
    call this%algo%rewind()
 
    call multisystem_debug_write_event_out(debug_handle, system_iteration = this%iteration, algo_iteration = this%algo%iteration)
 
    pop_sub(system_algorithm_start)
  end subroutine system_algorithm_start
 
  ! ---------------------------------------------------------
  subroutine system_algorithm_finish(this)
    class(system_t),      intent(inout) :: this
    type(event_handle_t) :: debug_handle
 
    character(len=:), allocatable :: updated_quantities(:)
 
    push_sub(system_algorithm_finish)
 
    debug_handle = multisystem_debug_write_event_in(this%namespace, event_function_call_t("system_algorithm_finish"), &
      system_iteration = this%iteration, algo_iteration = this%algo%iteration)
 
    ! Finish output
    call this%output_finish()
 
    ! System-specific and algorithm-specific finalization step
    if (this%algo%final_operation%id /= skip) then
      if (.not.  this%do_algorithmic_operation(this%algo%final_operation, updated_quantities)) then
        message(1) = "Unsupported algorithmic operation."
        write(message(2), '(A,A,A)') trim(this%algo%final_operation%id), ": ", trim(this%algo%final_operation%label)
        call messages_fatal(2, namespace=this%namespace)
      end if
      if (allocated(updated_quantities)) then
        message(1) = "Update of quantities not allowed in algorithmic operation."
        write(message(2), '(A,A,A)') trim(this%algo%final_operation%id), ": ", trim(this%algo%final_operation%label)
        call messages_fatal(2, namespace=this%namespace)
      end if
    end if
 
    call multisystem_debug_write_event_out(debug_handle, system_iteration = this%iteration, algo_iteration = this%algo%iteration)
 
    pop_sub(system_algorithm_finish)
  end subroutine system_algorithm_finish
 
  ! ---------------------------------------------------------
  subroutine system_iteration_info(this)
    class(system_t), intent(in) :: this
 
    real(real64) :: energy
    character(len=40) :: fmt
 
    push_sub(system_iteration_info)
 
    energy = units_from_atomic(units_out%energy, this%total_energy)
    if (abs(energy) >= 1e5) then
      fmt = '(i7,1x,f14.6,1X,es13.6,1X,i9,1X,'
    else
      fmt = '(i7,1x,f14.6,1X,f13.6,1X,i9,1X,'
    end if
    if (this%algo%elapsed_time < 1e-3) then
      fmt = trim(fmt)//'es13.3)'
    else
      fmt = trim(fmt)//'f13.3)'
    end if
 
    write(message(1), fmt) this%iteration%counter(), &
      units_from_atomic(units_out%time, this%iteration%value()), energy, &
      0, this%algo%elapsed_time
    call messages_info(1, namespace=this%namespace)
 
    pop_sub(system_iteration_info)
  end subroutine system_iteration_info
 
  ! ---------------------------------------------------------
  logical function system_process_is_slave(this)
    class(system_t), intent(in) :: this
 
    push_sub(system_process_is_slave)
 
    ! By default an MPI process is not a slave
    system_process_is_slave = .false.
 
    pop_sub(system_process_is_slave)
  end function system_process_is_slave
 
  ! ---------------------------------------------------------
  subroutine system_end(this)
    class(system_t), intent(inout) :: this
 
    type(interaction_iterator_t) :: iter
    class(interaction_t),              pointer :: interaction
 
    push_sub(system_end)
 
    ! No call to safe_deallocate macro here, as it gives an ICE with gfortran
    if (associated(this%algo)) then
      deallocate(this%algo)
    end if
 
    call iter%start(this%interactions)
    do while (iter%has_next())
      interaction => iter%get_next()
      if (associated(interaction)) then
        deallocate(interaction)
      end if
    end do
 
    pop_sub(system_end)
  end subroutine system_end
 
  ! ---------------------------------------------------------
  subroutine system_list_add_node(this, partner)
    class(system_list_t)         :: this
    class(interaction_partner_t), target :: partner
 
    push_sub(system_list_add_node)
 
    select type (partner)
    class is (system_t)
      call this%add_ptr(partner)
    class default
      assert(.false.)
    end select
 
    pop_sub(system_list_add_node)
  end subroutine system_list_add_node
 
  ! ---------------------------------------------------------
  recursive logical function system_list_contains(this, partner) result(contains)
    class(system_list_t)         :: this
    class(interaction_partner_t), target :: partner
 
    type(partner_iterator_t)  :: iterator
    class(interaction_partner_t),  pointer :: system
 
    push_sub(system_list_contains)
 
    contains = .false.
 
    select type (partner)
    class is (system_t)
 
      call iterator%start(this)
      do while (iterator%has_next() .and. .not. contains)
        system => iterator%get_next()
        contains = associated(system, partner)
      end do
 
    class default
      contains = .false.
    end select
 
    pop_sub(system_list_contains)
  end function system_list_contains
 
  ! ---------------------------------------------------------
  function system_iterator_get_next(this) result(system)
    class(system_iterator_t), intent(inout) :: this
    class(system_t),          pointer       :: system
 
    push_sub(system_iterator_get_next)
 
    select type (ptr => this%get_next_ptr())
    class is (system_t)
      system => ptr
    class default
      assert(.false.)
    end select
 
    pop_sub(system_iterator_get_next)
  end function system_iterator_get_next
 
  ! ---------------------------------------------------------
  subroutine system_init_parallelization(this, grp)
    class(system_t), intent(inout) :: this
    type(mpi_grp_t), intent(in)    :: grp
 
    push_sub(system_init_parallelization)
 
    call mpi_grp_copy(this%grp, grp)
    call messages_update_mpi_grp(this%namespace, grp)
 
    pop_sub(system_init_parallelization)
  end subroutine system_init_parallelization
 
 
 
  ! ---------------------------------------------------------
  subroutine system_start_barrier(this, target_time, barrier_index)
    class(system_t), intent(inout) :: this
    real(real64),    intent(in)    :: target_time
    integer,         intent(in)    :: barrier_index
 
    push_sub(system_start_barrier)
 
    this%barrier(barrier_index)%active = .true.
    this%barrier(barrier_index)%target_time = target_time
 
    pop_sub(system_start_barrier)
  end subroutine system_start_barrier
 
  ! ---------------------------------------------------------
  subroutine system_end_barrier(this, barrier_index)
    class(system_t), intent(inout) :: this
    integer,         intent(in)    :: barrier_index
 
    push_sub(system_end_barrier)
 
    this%barrier(barrier_index)%active = .false.
    this%barrier(barrier_index)%target_time = m_zero
 
    pop_sub(system_end_barrier)
  end subroutine system_end_barrier
 
  ! ---------------------------------------------------------
  logical function system_arrived_at_barrier(this, barrier_index)
    class(system_t), intent(inout) :: this
    integer,         intent(in)    :: barrier_index
 
    type(iteration_counter_t) :: iteration
 
    push_sub(system_arrived_at_barrier)
 
    system_arrived_at_barrier = .false.
    if (this%barrier(barrier_index)%active) then
      iteration = this%iteration + 1
      if (iteration%value() > this%barrier(barrier_index)%target_time) then
        system_arrived_at_barrier = .true.
      end if
    end if
 
    pop_sub(system_arrived_at_barrier)
  end function system_arrived_at_barrier
 
  ! ---------------------------------------------------------
  logical function system_arrived_at_any_barrier(this)
    class(system_t), intent(inout) :: this
 
    integer :: ii
 
    push_sub(system_arrived_at_any_barrier)
 
    system_arrived_at_any_barrier = .false.
    do ii = 1, number_barriers
      system_arrived_at_any_barrier = system_arrived_at_any_barrier &
        .or. this%arrived_at_barrier(ii)
    end do
 
    pop_sub(system_arrived_at_any_barrier)
  end function system_arrived_at_any_barrier
 
 
  ! ---------------------------------------------------------
  subroutine system_update_potential_energy(this)
    class(system_t), intent(inout) :: this
 
    type(interaction_iterator_t) :: iter
    class(interaction_t),              pointer :: interaction
 
    push_sub(system_update_potential_energy)
 
    this%potential_energy = m_zero
 
    call iter%start(this%interactions)
    do while (iter%has_next())
      interaction => iter%get_next()
      if(.not. interaction%intra_interaction) then
        call interaction%calculate_energy()
        this%potential_energy = this%potential_energy + interaction%energy
      end if
    end do
 
    pop_sub(system_update_potential_energy)
  end subroutine system_update_potential_energy
 
  ! ---------------------------------------------------------
  subroutine system_update_internal_energy(this)
    class(system_t), intent(inout) :: this
 
    type(interaction_iterator_t) :: iter
    class(interaction_t),              pointer :: interaction
 
    push_sub(system_update_internal_energy)
 
    this%internal_energy = m_zero
    call iter%start(this%interactions)
    do while (iter%has_next())
      interaction => iter%get_next()
      if(interaction%intra_interaction) then
        call interaction%calculate_energy()
        this%internal_energy = this%internal_energy + interaction%energy
      end if
    end do
 
    pop_sub(system_update_internal_energy)
  end subroutine system_update_internal_energy
 
  ! ---------------------------------------------------------
  subroutine system_update_total_energy(this)
    class(system_t), intent(inout) :: this
 
    push_sub(system_update_total_energy)
 
    call this%update_kinetic_energy()
    this%total_energy = this%kinetic_energy
 
    call this%update_potential_energy()
    this%total_energy = this%total_energy + this%potential_energy
 
    call this%update_internal_energy()
    this%total_energy = this%total_energy + this%internal_energy
 
    pop_sub(system_update_total_energy)
  end subroutine system_update_total_energy
 
end module system_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: