multisystem.F90

Go to the documentation of this file.

!! Copyright (C) 2019-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 multisystem_oct_m
  use algorithm_oct_m
  use algorithm_factory_oct_m
  use debug_oct_m
  use global_oct_m
  use ghost_interaction_oct_m
  use interaction_surrogate_oct_m
  use interaction_partner_oct_m
  use interaction_oct_m
  use interactions_factory_abst_oct_m
  use iteration_counter_oct_m
  use loct_oct_m
  use messages_oct_m
  use mpi_oct_m
  use multisystem_debug_oct_m
  use namespace_oct_m
  use profiling_oct_m
  use quantity_oct_m
  use system_oct_m
  implicit none
 
  private
  public ::               &
    multisystem_t,        &
    multisystem_end
 
  type, extends(system_t), abstract :: multisystem_t
    type(system_list_t) :: list
  contains
    procedure :: execute_algorithm =>  multisystem_execute_algorithm
    procedure :: init_parallelization => multisystem_init_parallelization
    procedure :: next_time_on_largest_dt => multisystem_next_time_on_largest_dt
    procedure :: reset_iteration_counters => multisystem_reset_iteration_counters
    procedure :: new_algorithm => multisystem_new_algorithm
    procedure :: algorithm_finished => multisystem_algorithm_finished
    procedure :: init_iteration_counters => multisystem_init_iteration_counters
    procedure :: add_partners_to_list => multisystem_add_partners_to_list
    procedure :: create_interactions => multisystem_create_interactions
    procedure :: algorithm_start => multisystem_algorithm_start
    procedure :: algorithm_finish => multisystem_algorithm_finish
    procedure :: init_interaction => multisystem_init_interaction
    procedure :: write_interaction_graph => multisystem_write_interaction_graph
    procedure :: initialize => multisystem_initialize
    procedure :: do_algorithmic_operation => multisystem_do_algorithmic_operation
    procedure :: is_tolerance_reached => multisystem_is_tolerance_reached
    procedure :: update_quantity => multisystem_update_quantity
    procedure :: init_interaction_as_partner => multisystem_init_interaction_as_partner
    procedure :: copy_quantities_to_interaction => multisystem_copy_quantities_to_interaction
    procedure :: process_is_slave => multisystem_process_is_slave
    procedure :: start_barrier => multisystem_start_barrier
    procedure :: end_barrier => multisystem_end_barrier
    procedure :: arrived_at_barrier => multisystem_arrived_at_barrier
    procedure :: restart_write => multisystem_restart_write
    procedure :: restart_read => multisystem_restart_read
    procedure :: restart_write_data => multisystem_restart_write_data
    procedure :: restart_read_data => multisystem_restart_read_data
    procedure :: update_kinetic_energy => multisystem_update_kinetic_energy
    procedure :: update_potential_energy => multisystem_update_potential_energy
    procedure :: update_internal_energy => multisystem_update_internal_energy
    procedure :: get_flat_list => multisystem_get_flat_list
  end type multisystem_t
 
contains
 
  ! ---------------------------------------------------------------------------------------
  recursive subroutine multisystem_init_parallelization(this, grp)
    class(multisystem_t), intent(inout) :: this
    type(mpi_grp_t),      intent(in)    :: grp
 
    type(system_iterator_t) :: iter
    class(system_t), pointer :: sys
    type(mpi_grp_t) :: sys_grp
 
    push_sub(multisystem_init_parallelization)
 
    call system_init_parallelization(this, grp)
 
    ! Now parallelize over systems in this multisystem
    call iter%start(this%list)
    do while (iter%has_next())
      sys => iter%get_next()
      ! for now, duplicate communicator - more complicated parallelization schemes can be implemented here
      call mpi_grp_duplicate(sys_grp, grp)
      call sys%init_parallelization(sys_grp)
    end do
 
    pop_sub(multisystem_init_parallelization)
  end subroutine multisystem_init_parallelization
 
  ! ---------------------------------------------------------------------------------------
  recursive function multisystem_next_time_on_largest_dt(this) result(next_time_on_largest_dt)
    class(multisystem_t), intent(inout) :: this
    real(real64) :: next_time_on_largest_dt
 
    type(system_iterator_t) :: iter
    class(system_t), pointer :: system
    type(iteration_counter_t) :: iteration
 
    push_sub(multisystem_next_time_on_largest_dt)
 
    next_time_on_largest_dt = m_zero
    call iter%start(this%list)
    do while (iter%has_next())
      system => iter%get_next()
      select type (system)
      class is (multisystem_t)
        next_time_on_largest_dt = max(next_time_on_largest_dt, system%next_time_on_largest_dt())
      class default
        iteration = system%iteration + 1
        next_time_on_largest_dt = max(next_time_on_largest_dt, iteration%value())
      end select
    end do
 
    pop_sub(multisystem_next_time_on_largest_dt)
  end function multisystem_next_time_on_largest_dt
 
  ! ---------------------------------------------------------------------------------------
  recursive subroutine multisystem_execute_algorithm(this)
    class(multisystem_t),     intent(inout) :: this
 
    type(system_iterator_t) :: iter
    class(system_t), pointer :: system
 
    type(event_handle_t) :: debug_handle
 
    push_sub(multisystem_execute_algorithm)
 
    debug_handle = multisystem_debug_write_event_in(this%namespace, event_function_call_t("multisystem_dt_operation"), &
      system_iteration = this%iteration, algo_iteration = this%algo%iteration)
 
    ! Multisystem
    call system_execute_algorithm(this)
 
    ! Subsystems
    call iter%start(this%list)
    do while (iter%has_next())
      system => iter%get_next()
      call system%execute_algorithm()
    end do
 
    call multisystem_debug_write_event_out(debug_handle, system_iteration = this%iteration, algo_iteration = this%algo%iteration)
 
    pop_sub(multisystem_execute_algorithm)
  end subroutine multisystem_execute_algorithm
 
  ! ---------------------------------------------------------------------------------------
  recursive subroutine multisystem_reset_iteration_counters(this, accumulated_iterations)
    class(multisystem_t),      intent(inout) :: this
    integer,                   intent(in)    :: accumulated_iterations
 
    type(system_iterator_t) :: iter
    class(system_t), pointer :: system
 
    push_sub(multisystem_reset_iteration_counters)
 
    ! Multisystem iteration counters
    call system_reset_iteration_counters(this, accumulated_iterations)
 
    ! Subsystems iteration counters
    call iter%start(this%list)
    do while (iter%has_next())
      system => iter%get_next()
      call system%reset_iteration_counters(accumulated_iterations)
    end do
 
    pop_sub(multisystem_reset_iteration_counters)
  end subroutine multisystem_reset_iteration_counters
 
  ! ---------------------------------------------------------------------------------------
  recursive subroutine multisystem_new_algorithm(this, factory)
    class(multisystem_t),       intent(inout) :: this
    class(algorithm_factory_t), intent(in)    :: factory
 
    type(system_iterator_t) :: iter
    class(system_t), pointer :: system
 
    push_sub(multisystem_new_algorithm)
 
    ! Now initialized the algorithms of the subsystems
    call iter%start(this%list)
    do while (iter%has_next())
      system => iter%get_next()
      call system%new_algorithm(factory)
    end do
 
    ! Initialize the algorithm of the multisystem. By default the
    ! multisystem itself and its own quantities are kept unchanged
    ! by using the static propagator. However, the subsystems are allowed to have
    ! their own propagators and those do not have to be static.
    ! Needs to be done after initializing the subsystems propagators,
    ! as we use the largest dt of the subsystems.
    this%algo => factory%create_static(this)
    call this%algo%rewind()
 
    call system_init_iteration_counters(this)
 
    pop_sub(multisystem_new_algorithm)
  end subroutine multisystem_new_algorithm
 
  ! ---------------------------------------------------------------------------------------
  recursive function multisystem_algorithm_finished(this) result(finished)
    class(multisystem_t),       intent(in) :: this
    logical :: finished
 
    type(system_iterator_t) :: iter
    class(system_t), pointer :: system
 
    ! Check if multisystem itself is finished
    finished = this%algo%finished()
 
    ! Check subsystems
    call iter%start(this%list)
    do while (iter%has_next())
      system => iter%get_next()
      finished = finished .and. system%algorithm_finished()
    end do
 
  end function multisystem_algorithm_finished
 
  ! ---------------------------------------------------------------------------------------
  recursive subroutine multisystem_init_iteration_counters(this)
    class(multisystem_t),       intent(inout) :: this
 
    type(system_iterator_t) :: iter
    class(system_t), pointer :: system
 
    push_sub(multisystem_init_iteration_counters)
 
    ! initialize multisystem iteration counters
    call system_init_iteration_counters(this)
 
    ! initialize iteration counters of subsystems
    call iter%start(this%list)
    do while (iter%has_next())
      system => iter%get_next()
      call system%init_iteration_counters()
    end do
 
    pop_sub(multisystem_init_iteration_counters)
  end subroutine multisystem_init_iteration_counters
 
  ! ---------------------------------------------------------------------------------------
  recursive subroutine multisystem_algorithm_start(this)
    class(multisystem_t),      intent(inout) :: this
 
    type(system_iterator_t) :: iter
    class(system_t), pointer :: system
 
    type(event_handle_t) :: debug_handle
 
    push_sub(multisystem_algorithm_start)
 
    debug_handle = multisystem_debug_write_event_in(this%namespace, event_function_call_t("multisystem_algorithm_start"), &
      system_iteration = this%iteration, algo_iteration = this%algo%iteration)
 
    ! Now start the execution of the subsystems
    call iter%start(this%list)
    do while (iter%has_next())
      system => iter%get_next()
      call system%algorithm_start()
    end do
 
    ! Once the subsystems are initialized start the propagation of the multisystem
    call system_algorithm_start(this)
 
    call multisystem_debug_write_event_out(debug_handle, system_iteration = this%iteration, algo_iteration = this%algo%iteration)
 
    pop_sub(multisystem_algorithm_start)
  end subroutine multisystem_algorithm_start
 
  ! ---------------------------------------------------------------------------------------
  recursive subroutine multisystem_algorithm_finish(this)
    class(multisystem_t),      intent(inout) :: this
 
    type(system_iterator_t) :: iter
    class(system_t), pointer :: system
 
    type(event_handle_t) :: debug_handle
 
    push_sub(multisystem_algorithm_finish)
 
    debug_handle = multisystem_debug_write_event_in(this%namespace, event_function_call_t("multisystem_algorithm_finish"), &
      system_iteration = this%iteration, algo_iteration = this%algo%iteration)
 
    ! Finish the execution of the multisystem
    call system_algorithm_finish(this)
 
    ! Now finish the execution of the subsystems
    call iter%start(this%list)
    do while (iter%has_next())
      system => iter%get_next()
      call system%algorithm_finish()
    end do
 
    call multisystem_debug_write_event_out(debug_handle, system_iteration = this%iteration, algo_iteration = this%algo%iteration)
 
    pop_sub(multisystem_algorithm_finish)
  end subroutine multisystem_algorithm_finish
 
  ! ---------------------------------------------------------
  recursive subroutine multisystem_add_partners_to_list(this, list, interaction_type)
    class(multisystem_t),             intent(in)    :: this
    class(partner_list_t),            intent(inout) :: list
    integer,                optional, intent(in)    :: interaction_type
 
    type(system_iterator_t) :: iter
    class(system_t), pointer :: system
 
    if (present(interaction_type)) then
      if (any(this%supported_interactions_as_partner == interaction_type)) then
        call list%add(this)
      end if
    else
      call list%add(this)
    end if
 
    call iter%start(this%list)
    do while (iter%has_next())
      system => iter%get_next()
      call system%add_partners_to_list(list, interaction_type)
    end do
 
  end subroutine multisystem_add_partners_to_list
 
  ! ---------------------------------------------------------
 
  recursive subroutine multisystem_create_interactions(this, interaction_factory, available_partners)
    class(multisystem_t),                  intent(inout) :: this
    class(interactions_factory_abst_t),    intent(in)    :: interaction_factory
    class(partner_list_t),         target, intent(in)    :: available_partners
 
    type(system_iterator_t) :: iter
    class(system_t), pointer :: subsystem
 
    ! Create the multisystem interactions
    call system_create_interactions(this, interaction_factory, available_partners)
 
    ! Create the subsystems interactions
    call iter%start(this%list)
    do while (iter%has_next())
      subsystem => iter%get_next()
      call subsystem%create_interactions(interaction_factory, available_partners)
    end do
 
  end subroutine multisystem_create_interactions
 
  ! ---------------------------------------------------------
 
  subroutine multisystem_init_interaction(this, interaction)
    class(multisystem_t), target, intent(inout) :: this
    class(interaction_t),         intent(inout) :: interaction
 
    push_sub(multisystem_init_interaction)
 
    ! The multisystem class should never know about any specific interaction.
    ! Only classes that extend it can know about specific interactions.
    ! Such classes should override this method to add new supported interactions.
    message(1) = "Trying to initialize an interaction in the multi-system container class"
    call messages_fatal(1, namespace=this%namespace)
 
    pop_sub(multisystem_init_interaction)
  end subroutine multisystem_init_interaction
 
  ! ---------------------------------------------------------------------------------------
  recursive subroutine multisystem_write_interaction_graph(this, iunit, include_ghosts)
    class(multisystem_t), intent(in) :: this
    integer,              intent(in) :: iunit
    logical,              intent(in) :: include_ghosts
 
    class(system_t), pointer :: system
    class(interaction_t),              pointer :: interaction
    type(system_iterator_t) :: sys_iter
    type(interaction_iterator_t) :: inter_iter
 
    push_sub(multisystem_write_interaction_graph)
 
    ! Loop over all the subsystems
    call sys_iter%start(this%list)
    do while (sys_iter%has_next())
      system => sys_iter%get_next()
 
      ! Loop over the interactions that this subsystem has
      call inter_iter%start(system%interactions)
      do while (inter_iter%has_next())
        interaction => inter_iter%get_next()
 
        ! Write interaction to DOT graph if this interaction has a partner
        select type (interaction)
        type is (ghost_interaction_t)
          if (include_ghosts) then
            write(iunit, '(2x,a)') '"' + trim(system%namespace%get()) + '" <- "' + trim(interaction%partner%namespace%get()) + &
              '" [label="'+ interaction%label + '"];'
          end if
          ! Do not include systems connected by ghost interactions
        class default
          write(iunit, '(2x,a)') '"' + trim(system%namespace%get()) + '" <- "' + trim(interaction%partner%namespace%get()) + &
            '" [label="'+ interaction%label + '"];'
        end select
      end do
 
      ! If this subsystem is also a multisystem, then we also need to traverse it
      select type (system)
      class is (multisystem_t)
        call system%write_interaction_graph(iunit, include_ghosts)
      end select
    end do
 
    pop_sub(multisystem_write_interaction_graph)
  end subroutine multisystem_write_interaction_graph
 
  ! ---------------------------------------------------------
  recursive subroutine multisystem_initialize(this)
    class(multisystem_t), intent(inout) :: this
 
    type(system_iterator_t) :: iter
    class(system_t), pointer :: system
 
    push_sub(multisystem_initialize)
 
    call iter%start(this%list)
    do while (iter%has_next())
      system => iter%get_next()
      call system%initialize()
    end do
 
    pop_sub(multisystem_initialize)
  end subroutine multisystem_initialize
 
  ! ---------------------------------------------------------
  logical function multisystem_do_algorithmic_operation(this, operation, updated_quantities) result(done)
    class(multisystem_t),           intent(inout) :: this
    class(algorithmic_operation_t), intent(in)    :: operation
    character(len=:),  allocatable, intent(out)   :: updated_quantities(:)
 
    push_sub(multisystem_do_algorithmic_operation)
 
    ! Currently there are no multisystem specific algorithmic operations
    done = .false.
 
    pop_sub(multisystem_do_algorithmic_operation)
  end function multisystem_do_algorithmic_operation
 
  ! ---------------------------------------------------------
  recursive logical function multisystem_is_tolerance_reached(this, tol) result(converged)
    class(multisystem_t), intent(in)    :: this
    real(real64),         intent(in)    :: tol
 
    type(system_iterator_t) :: iter
    class(system_t), pointer :: system
 
    push_sub(multisystem_is_tolerance_reached)
 
    converged = .true.
    call iter%start(this%list)
    do while (iter%has_next())
      system => iter%get_next()
      if (.not. system%is_tolerance_reached(tol)) converged = .false.
    end do
 
    pop_sub(multisystem_is_tolerance_reached)
  end function multisystem_is_tolerance_reached
 
  ! ---------------------------------------------------------
  subroutine multisystem_update_quantity(this, label)
    class(multisystem_t), intent(inout) :: this
    character(len=*), intent(in)    :: label
 
    push_sub(multisystem_update_quantity)
 
    ! The multisystem class should never know about any specific quantities.
    ! Only classes that extend it can know about specific quantities.
    ! Such classes should override this method to add new supported quantities.
    message(1) = "Trying to update a quantity in the multi-system container class"
    call messages_fatal(1, namespace=this%namespace)
 
    pop_sub(multisystem_update_quantity)
  end subroutine multisystem_update_quantity
 
  ! ---------------------------------------------------------
  subroutine multisystem_init_interaction_as_partner(partner, interaction)
    class(multisystem_t),           intent(in)    :: partner
    class(interaction_surrogate_t), intent(inout) :: interaction
 
    push_sub(multisystem_init_interaction_as_partner)
 
    ! The multisystem class should never know about any specific interaction.
    ! Only classes that extend it can know about specific interactions.
    ! Such classes should override this method to add new supported interactions.
    message(1) = "Trying to initialize an interaction as partner in the multi-system container class"
    call messages_fatal(1, namespace=partner%namespace)
 
    pop_sub(multisystem_init_interaction_as_partner)
  end subroutine multisystem_init_interaction_as_partner
 
  ! ---------------------------------------------------------
  subroutine multisystem_copy_quantities_to_interaction(partner, interaction)
    class(multisystem_t),           intent(inout) :: partner
    class(interaction_surrogate_t), intent(inout) :: interaction
 
    push_sub(multisystem_copy_quantities_to_interaction)
 
    ! The multisystem class should never know about any specific quantities.
    ! Only classes that extend it can know about specific quantities.
    ! Such classes should override this method to add new supported quantities.
    message(1) = "Trying to copy quantities to interaction in the multi-system container class"
    call messages_fatal(1, namespace=partner%namespace)
 
    pop_sub(multisystem_copy_quantities_to_interaction)
  end subroutine multisystem_copy_quantities_to_interaction
 
  ! ---------------------------------------------------------
  recursive logical function multisystem_process_is_slave(this) result(is_slave)
    class(multisystem_t), intent(in) :: this
 
    type(system_iterator_t) :: iter
    class(system_t), pointer :: system
 
    push_sub(multisystem_process_is_slave)
 
    is_slave = .false.
    call iter%start(this%list)
    do while (iter%has_next())
      system => iter%get_next()
      if (system%process_is_slave()) is_slave = .true.
    end do
 
    pop_sub(multisystem_process_is_slave)
  end function multisystem_process_is_slave
 
  !--------------------------------------------------------------------
  recursive subroutine multisystem_update_kinetic_energy(this)
    class(multisystem_t), intent(inout) :: this
 
    push_sub(multisystem_update_kinetic_energy)
 
    ! We currently do not have the center of mass coordinates implemented for multisystems,
    ! hence we set the kinetic energy to zero.
    ! The kinetic energies of the constituents are contributing to the internal energy.
 
    this%kinetic_energy = m_zero
 
    pop_sub(multisystem_update_kinetic_energy)
  end subroutine multisystem_update_kinetic_energy
 
  !---------------------------------------------------------
  recursive subroutine multisystem_update_internal_energy(this)
    class(multisystem_t), intent(inout) :: this
 
    class(system_t),              pointer :: system
    class(system_t),              pointer :: system_2
    type(system_iterator_t)      :: system_iter
    type(system_iterator_t)      :: system_iter_2
 
    push_sub(multisystem_update_internal_energy)
 
    ! The internal energy of the multisystem contains the kinetic and internal energies of the consistuents
    !TODO: the kinetic energy wrt the center of mass motion should be subtracted.
 
    this%internal_energy = m_zero
 
    call system_iter%start(this%list)
    do while (system_iter%has_next())
 
      system => system_iter%get_next()
 
      ! First add the kinetic energies of the subsystems
      call system%update_kinetic_energy()
      this%internal_energy = this%internal_energy + system%kinetic_energy
 
      ! First add the internal energies of the subsystems
      call system%update_internal_energy()
      this%internal_energy = this%internal_energy + system%internal_energy
 
      ! Now add the (inter-) interactions between the systems in the container.
      call system_iter_2%start(this%list)
      do while(system_iter_2%has_next())
 
        system_2 => system_iter_2%get_next()
 
        ! exclude self-interactions (intra-interactions) as they are included in the internal energy
        ! of the subsystem, which was already added above.
        if(.not. associated(system, system_2)) then
          this%internal_energy = this%internal_energy + multisystem_pair_energy(system, system_2)
        end if
      end do ! system_iter_2
 
    end do ! system_iter
 
    pop_sub(multisystem_update_internal_energy)
  end subroutine multisystem_update_internal_energy
 
  ! ---------------------------------------------------------
  subroutine multisystem_update_potential_energy(this)
    class(multisystem_t), intent(inout) :: this
 
    type(system_iterator_t)  :: system_iter
    class(system_t), pointer :: system
    type(interaction_iterator_t)  :: interaction_iter
    class(interaction_t),              pointer :: interaction
    type(system_list_t) :: flat_list
 
    push_sub(multisystem_update_potential_energy)
 
    this%potential_energy = m_zero
 
    ! We need to handle interactions of the container itself:
    call system_update_potential_energy(this)
 
    ! generate a list of all systems inside the container and its subcontainers:
    call this%get_flat_list(flat_list)
 
    ! loop over all systems inside the container
    call system_iter%start(flat_list)
    do while (system_iter%has_next())
 
      system => system_iter%get_next()
 
      ! Even though we are not using the potential energy of the subsystems here, we need to trigger their calculation
      call system%update_potential_energy()
 
      ! loop over all interactions and discard those with partners inside the container
      call interaction_iter%start(system%interactions)
      do while (interaction_iter%has_next())
        interaction => interaction_iter%get_next()
        if(.not. flat_list%contains(interaction%partner) .and. .not. interaction%intra_interaction) then
          call interaction%calculate_energy()
          this%potential_energy = this%potential_energy + interaction%energy
        end if
      end do
 
    end do
 
    pop_sub(multisystem_update_potential_energy)
  end subroutine multisystem_update_potential_energy
 
  ! ---------------------------------------------------------
  recursive real(real64) function multisystem_pair_energy(partner_A, partner_B) result(pair_energy)
    class(interaction_partner_t), intent(in) :: partner_a
    class(interaction_partner_t), intent(in) :: partner_b
 
    class(system_t), pointer :: system_a
    class(system_t), pointer :: system_b
    type(system_iterator_t) :: system_iterator_a
    type(system_iterator_t) :: system_iterator_b
 
    push_sub(multisystem_pair_energy)
 
    pair_energy = m_zero
 
    select type(partner_a)
    class is (multisystem_t) ! partner_A is container
 
      call system_iterator_a%start(partner_a%list)
      do while( system_iterator_a%has_next() )
 
        system_a => system_iterator_a%get_next()
 
        select type(partner_b)
        class is (multisystem_t)
 
          call system_iterator_b%start(partner_b%list)
          do while( system_iterator_b%has_next() )
            system_b => system_iterator_b%get_next()
            pair_energy = pair_energy + multisystem_pair_energy(system_a, system_b)
          end do
 
        class is (system_t)
          pair_energy = pair_energy + interaction_energy(partner_a, partner_b)
        class default
          assert(.false.) ! partner_A must be a system_t
        end select
      end do
 
    class is (system_t) ! partner_A is non-container system
 
      select type(partner_b)
      class is (multisystem_t) ! partner_B is container
 
        call system_iterator_b%start(partner_b%list)
        do while( system_iterator_b%has_next() )
          system_b => system_iterator_b%get_next()
          pair_energy = pair_energy + multisystem_pair_energy(partner_a, system_b)
        end do
 
      class default ! both partner_A and partner_B are explicit: we need to calculate
        pair_energy = pair_energy + interaction_energy(partner_a, partner_b)
      end select
 
    class default
      assert(.false.)
    end select
 
    pop_sub(multisystem_pair_energy)
 
  contains
 
    real(real64) function interaction_energy(system, partner) result (energy)
      class(system_t), target, intent(in) :: system
      class(interaction_partner_t), target, intent(in) :: partner
 
      type(interaction_iterator_t) :: interaction_iterator
      class(interaction_t),              pointer :: interaction
 
      energy = m_zero
 
      call interaction_iterator%start(system%interactions)
      do while(interaction_iterator%has_next())
        interaction => interaction_iterator%get_next()
        if( associated(interaction%partner, partner)) then
          call interaction%calculate_energy()
          energy = energy + interaction%energy
        end if
      end do
    end function interaction_energy
 
  end function multisystem_pair_energy
 
 
  ! ---------------------------------------------------------
  recursive subroutine multisystem_get_flat_list(this, flat_list)
    class(multisystem_t), intent(in) :: this
    type(system_list_t), intent(out) :: flat_list
 
    class(interaction_partner_t), pointer :: partner
    type(partner_iterator_t) :: iterator
 
    push_sub(multisystem_get_flat_list)
 
    call iterator%start(this%list)
    do while (iterator%has_next())
      partner => iterator%get_next()
 
      call flat_list%add(partner)
 
      select type (partner)
      class is (multisystem_t)
        ! Also include the subsystems of a multisystem
        call partner%get_flat_list(flat_list)
      end select
 
    end do
 
    pop_sub(multisystem_get_flat_list)
 
  end subroutine multisystem_get_flat_list
 
  ! ---------------------------------------------------------
  recursive subroutine multisystem_end(this)
    class(multisystem_t), intent(inout) :: this
 
    type(system_iterator_t) :: iter
    class(system_t), pointer :: system
 
    push_sub(multisystem_end)
 
    call iter%start(this%list)
    do while (iter%has_next())
      system => iter%get_next()
      if (associated(system)) then
        deallocate(system)
      end if
    end do
 
    call system_end(this)
 
    pop_sub(multisystem_end)
  end subroutine multisystem_end
 
  ! ---------------------------------------------------------
  recursive subroutine multisystem_start_barrier(this, target_time, barrier_index)
    class(multisystem_t), intent(inout) :: this
    real(real64),         intent(in)    :: target_time
    integer,              intent(in)    :: barrier_index
 
    type(system_iterator_t) :: iter
    class(system_t), pointer :: system
 
    push_sub(multisystem_start_barrier)
 
    call iter%start(this%list)
    do while (iter%has_next())
      system => iter%get_next()
      call system%start_barrier(target_time, barrier_index)
    end do
 
    pop_sub(multisystem_start_barrier)
  end subroutine multisystem_start_barrier
 
  ! ---------------------------------------------------------
  recursive subroutine multisystem_end_barrier(this, barrier_index)
    class(multisystem_t), intent(inout) :: this
    integer,              intent(in)    :: barrier_index
 
    type(system_iterator_t) :: iter
    class(system_t), pointer :: system
 
    push_sub(multisystem_end_barrier)
 
    call iter%start(this%list)
    do while (iter%has_next())
      system => iter%get_next()
      call system%end_barrier(barrier_index)
    end do
 
    pop_sub(multisystem_end_barrier)
  end subroutine multisystem_end_barrier
 
  ! ---------------------------------------------------------
  recursive logical function multisystem_arrived_at_barrier(this, barrier_index)
    class(multisystem_t), intent(inout) :: this
    integer,              intent(in)    :: barrier_index
 
    type(system_iterator_t) :: iter
    class(system_t), pointer :: system
 
    push_sub(multisystem_arrived_at_barrier)
 
    multisystem_arrived_at_barrier = .true.
    call iter%start(this%list)
    do while (iter%has_next())
      system => iter%get_next()
      multisystem_arrived_at_barrier = multisystem_arrived_at_barrier .and. &
        system%arrived_at_barrier(barrier_index)
    end do
 
    pop_sub(multisystem_arrived_at_barrier)
  end function multisystem_arrived_at_barrier
 
  ! ---------------------------------------------------------
  recursive subroutine multisystem_restart_write(this)
    class(multisystem_t), intent(inout) :: this
 
    type(system_iterator_t) :: iter
    class(system_t), pointer :: system
 
    push_sub(multisystem_restart_write)
 
    ! do generic restart steps
    call system_restart_write(this)
 
    ! loop over all subsystems
    call iter%start(this%list)
    do while (iter%has_next())
      system => iter%get_next()
      call system%restart_write()
    end do
    message(1) = "Wrote restart data for multisystem "//trim(this%namespace%get())
    call messages_info(1, namespace=this%namespace)
 
    pop_sub(multisystem_restart_write)
  end subroutine multisystem_restart_write
 
  ! ---------------------------------------------------------
  recursive logical function multisystem_restart_read(this)
    class(multisystem_t), intent(inout) :: this
 
    type(system_iterator_t) :: iter
    class(system_t), pointer :: system
 
    push_sub(multisystem_restart_read)
 
    ! read generic restart data
    multisystem_restart_read = system_restart_read(this)
    call iter%start(this%list)
    do while (iter%has_next())
      system => iter%get_next()
      ! TODO: adapt logics here for consistent restarting
      multisystem_restart_read = multisystem_restart_read .and. &
        system%restart_read()
    end do
 
    if (multisystem_restart_read) then
      message(1) = "Successfully read restart data for multisystem "//trim(this%namespace%get())
      call messages_info(1, namespace=this%namespace)
    end if
 
    pop_sub(multisystem_restart_read)
  end function multisystem_restart_read
 
  ! ---------------------------------------------------------
  subroutine multisystem_restart_write_data(this)
    class(multisystem_t), intent(inout) :: this
 
    push_sub(multisystem_restart_write_data)
 
    ! do not write restart data for multisystem_t
 
    pop_sub(multisystem_restart_write_data)
  end subroutine multisystem_restart_write_data
 
  ! ---------------------------------------------------------
  ! this function returns true if restart data could be read
  logical function multisystem_restart_read_data(this)
    class(multisystem_t), intent(inout) :: this
 
    push_sub(multisystem_restart_read_data)
 
    multisystem_restart_read_data = .true.
 
    pop_sub(multisystem_restart_read_data)
  end function multisystem_restart_read_data
 
end module multisystem_oct_m