ensemble.F90

Go to the documentation of this file.

!! Copyright (C) 2024 S. de la Pena, M. Lueders
!!
!! This Source Code Form is subject to the terms of the Mozilla Public
!! License, v. 2.0. If a copy of the MPL was not distributed with this
!! file, You can obtain one at https://mozilla.org/MPL/2.0/.
!!
#include "global.h"
 
!
module ensemble_oct_m
 
  use algorithm_oct_m
  use debug_oct_m
  use global_oct_m
  use interactions_factory_abst_oct_m
  use interaction_partner_oct_m
  use io_oct_m
  use messages_oct_m
  use multisystem_basic_oct_m
  use multisystem_oct_m
  use namespace_oct_m
  use system_factory_abst_oct_m
  use system_oct_m
  use parser_oct_m
  use profiling_oct_m
 
  implicit none
 
  private
  public ::               &
    ensemble_t
 
  type, extends(multisystem_t) :: ensemble_t
    integer                   :: nreplica
    !                                        !!
    !                                        !! This is the number of replicas in, planned for a calculation.
    !                                        !! The actual number of replicas in one run is determined by
    !                                        !! last-first+1.
    !                                        !!
    !                                        !! The total number of replicas can be different can differ if
    !                                        !! further runs have been added to improve the convergence.
    integer                   :: first
    integer                   :: last
    real(real64), allocatable :: weights(:)
 
  contains
    procedure :: create_interactions => ensemble_create_interactions
    procedure :: add_partners_to_list => ensemble_add_partners_to_list
    final :: ensemble_finalizer
  end type ensemble_t
 
  interface ensemble_t
    procedure ensemble_constructor
  end interface ensemble_t
 
contains
 
  ! ---------------------------------------------------------
  !
  recursive function ensemble_constructor(namespace, n_replicas, first, last, system_factory, names, types, calc_mode_id) &
    result(sys)
    type(namespace_t),                  intent(in) :: namespace
    integer,                            intent(in) :: n_replicas
    class(system_factory_abst_t),       intent(in) :: system_factory
    integer,                            intent(in) :: first
    integer,                            intent(in) :: last
    character(len=128),                 intent(in) :: names(:)
    integer,                            intent(in) :: types(:)
    integer,                            intent(in) :: calc_mode_id
    class(ensemble_t), pointer    :: sys
 
    integer                  :: ireplica
    class(system_t), pointer :: replica
    character(len=128)       :: replica_name
 
    push_sub(ensemble_constructor)
 
    allocate(sys)
    sys%namespace = namespace
    sys%nreplica = n_replicas
    sys%first = max(1, first)
    sys%last = min(n_replicas, last)
 
 
    ! No interaction directly supported by this system (but classes that extend it can add their own)
    allocate(sys%supported_interactions(0))
    allocate(sys%supported_interactions_as_partner(0))
 
    do ireplica = sys%first, sys%last
      write(replica_name, '(I8.8)') ireplica
 
      ! Be aware that we need to use `sys%namespace` here. Passing only `namespace` leads to an error,
      ! as `namespace` seems to be deleted when the routine exits.
      !
      replica => multisystem_basic_t(namespace_t(replica_name, parent=sys%namespace), names, types, system_factory, calc_mode_id)
 
      call sys%list%add(replica)
 
    end do
 
    ! Initialize weights to uniform distribution
 
    safe_allocate(sys%weights(1:sys%nreplica))
    sys%weights = 1.0_real64/sys%nreplica
 
    pop_sub(ensemble_constructor)
  end function ensemble_constructor
 
  recursive subroutine ensemble_finalizer(this)
    type(ensemble_t), intent(inout) :: this
 
    push_sub(ensemble_finalizer)
 
    call ensemble_end(this)
 
    pop_sub(ensemble_finalizer)
  end subroutine ensemble_finalizer
 
  recursive subroutine ensemble_end(this)
    type(ensemble_t), intent(inout) :: this
 
    push_sub(ensemble_end)
 
    safe_deallocate_a(this%weights)
    call multisystem_end(this)
 
    pop_sub(ensemble_end)
  end subroutine ensemble_end
 
  ! ---------------------------------------------------------
  !
  recursive subroutine ensemble_create_interactions(this, interaction_factory, available_partners)
    class(ensemble_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
    type(partner_list_t) :: partners
 
    push_sub(ensemble_create_interactions)
 
    ! Create interactions of the ensemble container
    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()
      ! partners = available_partners
      ! the above line would provide inter-replica ghost interactions, which does not seem
      ! to be necessary, as the replicas can propagate independent of each other.
      ! If we need more synchronization, e.g. for averages, this should be provided by the
      ! real interaction.
      ! In order to restore the inter-replica ghost, uncomment the line above and remove
      ! the line below (partners.empty()).
      call partners%empty()
      call partners%add(subsystem)
      call subsystem%create_interactions(interaction_factory, partners)
    end do
 
    pop_sub(ensemble_create_interactions)
  end subroutine ensemble_create_interactions
 
 
  recursive subroutine ensemble_add_partners_to_list(this, list, interaction_type)
    class(ensemble_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
 
    push_sub(ensemble_add_partners_to_list)
 
    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)
 
      call iter%start(this%list)
      do while (iter%has_next())
        system => iter%get_next()
        call system%add_partners_to_list(list)
      end do
 
    end if
 
    pop_sub(ensemble_add_partners_to_list)
  end subroutine ensemble_add_partners_to_list
 
end module ensemble_oct_m