Navigation :

Time propagation

Work in progress!

In the new multisystem framework, the time propagation is handled by the routine time_dependent_run_multisystem()

The main loop

Implementation of time_dependent_run_multisystem()

  subroutine time_dependent_run_multisystem(systems, from_scratch)
    class(multisystem_basic_t), intent(inout) :: systems
    logical,              intent(in)    :: from_scratch

    integer :: internal_loop
    integer, parameter :: MAX_PROPAGATOR_STEPS = 1000
    FLOAT :: final_time

    PUSH_SUB(time_dependent_run_multisystem)

    call multisystem_debug_init("debug/multisystem_propagation.log", global_namespace, systems%grp)

    call messages_write('Info: Running Multi-System time evolution')
    call messages_new_line()
    call messages_new_line()
    call messages_info()

    ! Get final propagation time from input
    ! This variable is also defined (and properly documented) in td/td.F90.
    ! This is temporary, until all the propagators are moved to the new framework.
    call parse_variable(systems%namespace, 'TDPropagationTime', CNST(-1.0), final_time, unit = units_inp%time)
    if (final_time <= M_ZERO) then
      call messages_input_error(systems%namespace, 'TDPropagationTime', 'must be greater than zero')
    end if
    call messages_print_var_value(stdout, 'TDPropagationTime', final_time)

    ! Initialize all propagators
    call systems%init_propagator()

    ! Set initial conditions
    call systems%initial_conditions(from_scratch)

    call messages_print_stress(stdout, "Multi-system propagation", namespace=systems%namespace)

    call systems%propagation_start()

    call multisystem_debug_write_marker(systems%namespace, event_marker_t("propagation_start"))

    ! The full TD loop
    do while (.not. systems%has_reached_final_propagation_time(final_time))
      internal_loop = 1
      do while (internal_loop < MAX_PROPAGATOR_STEPS)
        ! Do one algorithmic operation
        call systems%dt_operation()

        ! Exit loop if a propagation step is done
        if (systems%prop%step_is_done()) exit

        internal_loop = internal_loop + 1
      end do

      write (message(1), '(a)') repeat ('-', 71)
      call messages_info(1)

    end do

    call multisystem_debug_write_marker(systems%namespace, event_marker_t("propagation_finish"))

    call systems%propagation_finish()

    call multisystem_debug_end()

    POP_SUB(time_dependent_run_multisystem)
  end subroutine time_dependent_run_multisystem

The code is quite self-explanatory. The central part of the time propagation is the do ... while loop (see "! The full TD loop"). Here the function system%dt_operation() is called successively, triggering one algorithmic step of the multisystem propagator.

The propagators in general will need several calls of this to finish one time step of the propagator!

At the top level, there is only one system (in general called .), which is of type multisystem_basic_t, which then contains other systems as subsystems. These subsystems can be multisystems themselves.

Updating the system

In the multisystem framework, the highest level system is of type multisystem_basic_t, which inherits the dt_operation() routine from the abstract multisystem_t type.

Implementation of multisystem_dt_operation

This routine first calls the general system_dt_operation() subroutine from the parent class system_t and then loops over all subsystems, which are part of the multisystem, and calls their specific system%dt_operations() routine.

The general system_dt_operation() subroutine is handling all general algorithmic operations, such as starting and ending the SCF loop, or updating interaction.

Implementation of system_dt_operation()

  subroutine system_dt_operation(this)
    class(system_t),     intent(inout) :: this

    type(algorithmic_operation_t) :: tdop
    logical :: all_updated

    type(event_handle_t) :: debug_handle

    PUSH_SUB(system_dt_operation)

    tdop = this%prop%get_td_operation()

    if (debug%info) then
      write(message(1), '(a,a,1X,a)') "Debug: Start  ", trim(tdop%label), " for '" + trim(this%namespace%get()) + "'"
      call messages_info(1)
    end if

    debug_handle = multisystem_debug_write_event_in(this%namespace, event_function_call_t("dt_operation", tdop),    &
      system_clock=this%clock, prop_clock=this%prop%clock)

    select case (tdop%id)
    case (FINISHED)
      if (.not. this%prop%step_is_done()) then
        ! Increment the system clock by one time-step
        this%clock = this%clock + CLOCK_TICK
        call multisystem_debug_write_marker(this%namespace, event_clock_update_t("system",  "", this%clock, "tick" ) )

        ! Execute additional operations at the end of the time step
        call this%exec_end_of_timestep_tasks()

        ! Print information about the current iteration and write output
        call this%output_write()
        call this%iteration_info()

        ! Mark propagation step as finished
        call this%prop%finished()
      else
        ! Reset propagator for next step
        call this%prop%rewind()
      end if

    case (UPDATE_INTERACTIONS)
      ! We increment by one algorithmic step
      this%prop%clock = this%prop%clock + CLOCK_TICK
      call multisystem_debug_write_marker(this%namespace, event_clock_update_t("propagator", "", this%prop%clock, "tick") )

      ! Try to update all the interactions
      all_updated = this%update_interactions()

      ! Move to next propagator step if all interactions have been
      ! updated. Otherwise try again later.
      if (all_updated) then
        this%accumulated_loop_ticks = this%accumulated_loop_ticks + 1
        call this%prop%next()
      else
        this%prop%clock = this%prop%clock - CLOCK_TICK
        call multisystem_debug_write_marker(this%namespace, event_clock_update_t("propagator", "", this%prop%clock, "reverse") )
      end if

    case (START_SCF_LOOP)
      ASSERT(this%prop%predictor_corrector)

      call this%prop%save_scf_start()
      this%prop%inside_scf = .true.
      this%accumulated_loop_ticks = 0

      if (debug%info) then
        write(message(1), '(a,i3,a)') "Debug: -- SCF iter ", this%prop%scf_count, " for '" + trim(this%namespace%get()) + "'"
        call messages_info(1)
      end if

    case (END_SCF_LOOP)
      ! Here we first check if we did the maximum number of steps.
      ! Otherwise, we need check the tolerance
      if (this%prop%scf_count == this%prop%max_scf_count) then
        if (debug%info) then
          message(1) = "Debug: -- Max SCF Iter reached for '" + trim(this%namespace%get()) + "'"
          call messages_info(1)
        end if
        this%prop%inside_scf = .false.
        call this%prop%next()
      else
        ! We reset the pointer to the begining of the scf loop
        if (this%is_tolerance_reached(this%prop%scf_tol)) then
          if (debug%info) then
            message(1) = "Debug: -- SCF tolerance reached for '" + trim(this%namespace%get()) + "'"
            call messages_info(1)
          end if
          this%prop%inside_scf = .false.
          call this%prop%next()
        else
          ! We rewind the instruction stack
          call this%prop%rewind_scf_loop()

          ! We reset the clocks
          call this%reset_clocks(this%accumulated_loop_ticks)
          this%accumulated_loop_ticks = 0
          if (debug%info) then
            write(message(1), '(a,i3,a,a)') "Debug: -- SCF iter ", this%prop%scf_count, " for '" + trim(this%namespace%get()), "'"
            call messages_info(1)
          end if
        end if
      end if

    case default
      call this%do_td_operation(tdop)
      call this%prop%next()
    end select

    if (debug%info) then
      write(message(1), '(a,a,1X,a, l)') "Debug: Finish ", trim(tdop%label), " for '" + trim(this%namespace%get()) + "' ", &
        this%prop%step_is_done()
      call messages_info(1)
    end if

    call multisystem_debug_write_event_out(debug_handle, system_clock=this%clock, prop_clock=this%prop%clock)

    POP_SUB(system_dt_operation)
  end subroutine system_dt_operation

This routine fetches the next operation from the propagator and executes it. The following operations are implemented here:

UPDATE_INTRERACTIONS
- increment the propagator-clock by one tick (if all interactions can be updated)
- try to update all interactions (this might not always be possible, if some interaction partners are still lagging behind)
- if all interactions are updated succesfully, progress the propagator to the next step.
FINISHED
- increment the system-clock by one tick.
- rewind the propagator
START_SCF_LOOP and END_SCF_LOOP
- handles self-consistency and convergence for predictor-corrector algorithms.
- no clocks are updated.

All remaining operations are passed down to the system specific routine system%do_dt_operation(), which is the place where specific systems need to implement all algorithmic steps, required by the propagators, allowed for that system.

In a multisystem, the multisystem itself as well as all its subsystems are propagating according to their own propagators. In principle, they could be different, which could lead to unexpected behaviour. If propagators of the subsystems are not specified explicitely in the input file, they inherit automatically the propagator of the parent multisystem. Therefore, they should not be specified individually. The multisystem container itself always uses an ‘empty’ propagator, which only implements the operations OP_UPDATE_INTERACTIONS and OP_FINISHED. For multisystem_t objects, the do_td_operations() routine only implements the SKIP operation.

This routine is triggering the update of the interactions, via the call to system%update_interations(), which loops over the interactions associated with the system, and all required exposed quantities.

Updating the interactions

Each of the systems (i.e. the multisystem and its subsystems) are attempting to update their interactions with system%update_interactions().

Implementation of system%update_interations()

  logical function system_update_interactions(this) result(all_updated)
    class(system_t),      intent(inout) :: this

    logical :: none_updated
    integer :: iq, q_id
    class(interaction_t), pointer :: interaction
    type(interaction_iterator_t) :: iter

    PUSH_SUB(system_update_interactions)

    ! Some systems might need to perform some specific operations before the
    ! update. This should only be done if no interaction has been updated yet,
    ! so that it is only done once.
    none_updated = .true.
    call iter%start(this%interactions)
    do while (iter%has_next())
      interaction => iter%get_next()
      if (interaction%clock == this%prop%clock) then
        none_updated = .false.
        exit
      end if
    end do
    if (none_updated) then
      call this%update_interactions_start()
    end if

    !Loop over all interactions
    all_updated = .true.
    call iter%start(this%interactions)
    do while (iter%has_next())
      interaction => iter%get_next()

      if (.not. interaction%clock == this%prop%clock) then
        ! Update the system quantities that will be needed for computing the interaction
        do iq = 1, interaction%n_system_quantities
          ! Get requested quantity ID
          q_id = interaction%system_quantities(iq)

          ! All needed quantities must have been marked as required. If not, then fix your code!
          ASSERT(this%quantities(q_id)%required)

          ! We do not need to update the protected quantities, the propagator takes care of that
          if (this%quantities(q_id)%protected) cycle

          if (.not. this%quantities(q_id)%clock == this%prop%clock) then
            ! The requested quantity is not at the requested time, so we try to update it

            ! Sanity check: it should never happen that the quantity is in advance
            ! with respect to the requested time.
            if (this%quantities(q_id)%clock > this%prop%clock) then
              message(1) = "The quantity clock is in advance compared to the requested time."
              call messages_fatal(1, namespace=this%namespace)
            end if

            call this%update_quantity(q_id)
          end if

        end do

        ! We can now try to update the interaction
        all_updated = interaction%update(this%prop%clock) .and. all_updated
      end if
    end do

    ! Some systems might need to perform some specific operations after all the
    ! interactions have been updated
    if (all_updated) then
      call this%update_interactions_finish()
    end if

    POP_SUB(system_update_interactions)
  end function system_update_interactions

The first part makes sure that the update_interactions_start() routine is only called when no interaction has been updated yet in this time step, iu.e. if their clocks are behind the system clock.

It is assumed that the interaction can never be ahead of time, compared to the propagator. It is therefore sufficient to ask whether the interaction time equals the propagator time to determine whether an interaction is up-to-date.

In the second part, we actually attempt to update the interactions, if needed. If an interaction is not up-to-date, it first needs to be ensured that the quantities, on which the interaction depends, are up-to-date. Once all quantities are updated, the update() routine of the interaction can be called with the current propagator time as argument.

Finally, if all interactions are updated successfully, there is a step update_interactions_finish(), which can be necessary for some systems.

Implementation of interaction_with_partner_update()

  logical function interaction_with_partner_update(this, requested_time) result(updated)
    class(interaction_with_partner_t), intent(inout) :: this
    class(clock_t),                    intent(in)    :: requested_time

    logical :: allowed_to_update
    type(event_handle_t) :: debug_handle

    PUSH_SUB(interaction_with_partner_update)

    ! We should only try to update the interaction if it is not yet at the requested time
    ASSERT(.not. (this%clock == requested_time))

    debug_handle = multisystem_debug_write_event_in(event = event_function_call_t("interaction_with_partner_update"), &
      extra="target: "//trim(this%label)//"-"//trim(this%partner%namespace%get()), &
      interaction_clock = this%clock,      &
      partner_clock = this%partner%clock,  &
      requested_clock = requested_time)

    allowed_to_update = this%partner%update_exposed_quantities(requested_time, this)

    if (allowed_to_update) then
      ! We can now compute the interaction from the updated quantities
      call this%calculate()

      ! Update was successful, so set new interaction time
      updated = .true.
      call this%clock%set_time(requested_time)
      call multisystem_debug_write_marker(event = event_clock_update_t( "interaction", &
        trim(this%label)//"-"//trim(this%partner%namespace%get()), &
        this%clock, "set") )

      if (debug%info) then
        write(message(1), '(a,a,a,a,a)') "Debug: -- Updated '", trim(this%label), "' interaction with '", &
          trim(this%partner%namespace%get()), "'"
        write(message(2), '(a,f16.6,a,f16.6)') "Debug: ---- Requested time is ", requested_time%time(), &
          " and partner time is ", this%partner%clock%time()
        call messages_info(2)
      end if
    else
      if (debug%info) then
        write(message(1), '(a,a,a,a,a)') "Debug: -- Cannot update yet the '", trim(this%label), "' interaction with '", &
          trim(this%partner%namespace%get()), "'"
        write(message(2), '(a,f16.6,a,f16.6)') "Debug: ---- Requested time is ", requested_time%time(), &
          " and partner time is ", this%partner%clock%time()
        call messages_info(2)
      end if
      updated = .false.
    end if

    call multisystem_debug_write_event_out(debug_handle, update = updated,   &
      interaction_clock = this%clock,  &
      partner_clock = this%partner%clock, &
      requested_clock = requested_time)

    POP_SUB(interaction_with_partner_update)
  end function interaction_with_partner_update

This function is to be called with for a specific requested_time, which is the time at which the interaction is needed.

The function first tries to update the exposed quantities of the partner. If this is successfull, the function calls the calculate() routine of the interaction, sets the interaction clock to the requested time and returns with the result .true., indicating a successful update. In case the function was blocked by one of the exposed quantities, the interaction is not calculated and the function returns .false..

As can be seen, it is, possible that an interaction cannot be updated at a given time. This can be the case when the interaction also depends on quantities of the other interaction partner, and that partner is not yet at the requested time.

Definition of system_update_exposed_quantities()

  logical function system_update_exposed_quantities(partner, requested_time, interaction) result(allowed_to_update)
    class(system_t),      intent(inout) :: partner
    type(clock_t),        intent(in)    :: requested_time
    class(interaction_t), intent(inout) :: interaction

    logical :: ahead_in_time, right_on_time, need_to_copy
    integer :: iq, q_id

    type(event_handle_t) :: debug_handle

    PUSH_SUB(system_update_exposed_quantities)

    if (debug%info) then
      write(message(1), '(a,a,a)') "Debug: -- Updating exposed quantities for partner '", trim(partner%namespace%get()), "'"
      call messages_info(1)
    end if

    debug_handle = multisystem_debug_write_event_in(system_namespace = partner%namespace, &
      event = event_function_call_t("system_update_exposed_quantities"), &
      partner_clock = partner%clock, &
      requested_clock = requested_time, &
      interaction_clock = interaction%clock )

    select type (interaction)
    class is (interaction_with_partner_t)

      if (partner%prop%inside_scf .or. partner%prop%clock + CLOCK_TICK < requested_time) then
        ! we are inside an SCF cycle and therefore are not allowed to expose any quantities.
        ! or we are too much behind the requested time
        allowed_to_update = .false.
      else
        allowed_to_update = .true.
        need_to_copy = .true.
        do iq = 1, interaction%n_partner_quantities
          ! Get the requested quantity ID
          q_id = interaction%partner_quantities(iq)

          ! All needed quantities must have been marked as required. If not, then fix your code!
          ASSERT(partner%quantities(q_id)%required)

          ! First update the exposed quantities that are not protected
          if (.not.partner%quantities(q_id)%protected) then
            if (partner%quantities(q_id)%clock /= requested_time .and. &
              partner%quantities(q_id)%clock + CLOCK_TICK <= requested_time) then
              ! We can update because the partner will reach this time in the next sub-timestep
              call partner%update_exposed_quantity(q_id)

              call updated_quantity_debug()
            else
              call not_updated_quantity_debug()
            end if
          else
            call protected_quantity_debug()
          end if


          ! Now compare the times
          ahead_in_time = partner%quantities(q_id)%clock > requested_time
          right_on_time = partner%quantities(q_id)%clock == requested_time

          select case (partner%interaction_timing)
          case (OPTION__INTERACTIONTIMING__TIMING_EXACT)
            ! only allow interaction at exactly the same time
            allowed_to_update = allowed_to_update .and. right_on_time
            need_to_copy = allowed_to_update
          case (OPTION__INTERACTIONTIMING__TIMING_RETARDED)
            ! allow retarded interaction
            allowed_to_update = allowed_to_update .and. &
              (right_on_time .or. ahead_in_time)
            need_to_copy = need_to_copy .and. .not. ahead_in_time
          case default
            call messages_not_implemented("Method for interaction quantity timing")
          end select

        end do

        ! If the quantities have been updated, we copy them to the interaction
        if (allowed_to_update .and. need_to_copy) then
          select type (interaction)
          type is (ghost_interaction_t)
            ! Nothing to copy. We still need to check that we are at the right
            ! time for the update though!
          class default
            call partner%copy_quantities_to_interaction(interaction)
          end select
        end if
      end if

    class default
      message(1) = "A system can only expose quantities to an interaction as a partner."
      call messages_fatal(1, namespace=partner%namespace)
    end select

    if (debug%info) then
      write(message(1), '(a,a,a)') "Debug: -- Finished updating exposed quantities for partner '", &
        trim(partner%namespace%get()), "'"
      call messages_info(1)
    end if

    call multisystem_debug_write_event_out(debug_handle, update=allowed_to_update, &
      partner_clock = partner%clock, &
      requested_clock = requested_time, &
      interaction_clock = interaction%clock )

    POP_SUB(system_update_exposed_quantities)
  contains

    subroutine updated_quantity_debug()

      if (debug%info) then
        write(message(1), '(a,a,a)') "Debug: ---- Updated exposed quantity ", trim(QUANTITY_LABEL(q_id)), "'"
        write(message(2), '(a,f16.6,a,f16.6)') "Debug: ------ Requested time is ", requested_time%time(), &
          ", quantity time is ", partner%quantities(q_id)%clock%time()
        call messages_info(2)
      end if

    end subroutine updated_quantity_debug

    subroutine not_updated_quantity_debug()

      if (debug%info) then
        write(message(1), '(a,a,a)') "Debug: ---- Did not update exposed quantity '", trim(QUANTITY_LABEL(q_id)), "'"
        write(message(2), '(a,f16.6,a,f16.6,a,f16.6)') "Debug: ------ Requested time is ", requested_time%time(), &
          ", quantity time is ", partner%quantities(q_id)%clock%time(), &
          " and partner propagator time is ", partner%prop%clock%time()
        call messages_info(2)
      end if

    end subroutine not_updated_quantity_debug

    subroutine protected_quantity_debug()

      if (debug%info) then
        write(message(1), '(a,a,a)') "Debug: ---- Skip update of quantity '", trim(QUANTITY_LABEL(q_id)), &
          "' as it is protected"
        write(message(2), '(a,f16.6,a,f16.6)') "Debug: ------ Requested time is ", requested_time%time(), &
          ", quantity time is ", partner%quantities(q_id)%clock%time()
        call messages_info(2)
      end if

    end subroutine protected_quantity_debug

  end function system_update_exposed_quantities

This routine demands some more comments:

Firstly, the routine is only allowed if the interaction is of type interaction_with_partner_t, as otherwise there is no partner to update exposed quantities. It is important to remember who called the function, and which quantities are supposed to be updated:

system_update_exposed_quantities() is called from Implementation of interaction_with_partner_update() for each interaction of a system for the interaction partner of that system with respect to a given interaction.

In the following situations, the routine is not allowd to update the exposed quantities:

we are in an internal SCF loop of the propagation step
one of the quantities cannot be updated becuase it is either
- too far behind in time (where two far is more than one clock tick)
- it is one clock tick behind, but the user requested OPTION__INTERACTIONTIMING__TIMING_EXACT In these cases, the routine will return a .false.. The propagator will continue, but the system in question will be held at this step. Other systems, however, can continue to progress, and might be able to get to the time, where then this quantity can be updated.

If the timing conditions are fulfilled (and we are not in an internal SCF loop) the system will call the specific routines to update the exposed quantities, and copies their values to the interaction.

Time propagation

The main loop

Updating the system

Updating the interactions

Deadlocks and how to avoid them