dispersive_medium.F90

Go to the documentation of this file.

!! Copyright (C) 2021 F. Bonafé
!!
!! 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 dispersive_medium_oct_m
  use algorithm_oct_m
  use algorithm_factory_oct_m
  use calc_mode_par_oct_m
  use current_to_mxll_field_oct_m
  use debug_oct_m
  use field_transfer_oct_m
  use global_oct_m
  use interaction_oct_m
  use interaction_surrogate_oct_m
  use interaction_enum_oct_m
  use interaction_oct_m
  use interaction_partner_oct_m
  use io_oct_m
  use messages_oct_m
  use grid_oct_m
  use, intrinsic :: iso_fortran_env
  use lalg_basic_oct_m
  use lattice_vectors_oct_m
  use linear_medium_to_em_field_oct_m
  use linear_medium_oct_m
  use mesh_oct_m
  use mpi_oct_m
  use multicomm_oct_m
  use mxll_e_field_to_matter_oct_m
  use namespace_oct_m
  use output_low_oct_m
  use parser_oct_m
  use propagator_rk4_oct_m
  use propagator_oct_m
  use profiling_oct_m
  use quantity_oct_m
  use regridding_oct_m
  use restart_oct_m
  use space_oct_m
  use system_oct_m
  use unit_system_oct_m
  use unit_oct_m
  use utils_oct_m
  use varinfo_oct_m
  use write_iter_oct_m
 
  implicit none
 
  private
  public ::           &
    dispersive_medium_t,    &
    dispersive_medium_init
 
  type, extends(system_t) :: dispersive_medium_t
    private
    type(space_t)         :: space
    real(real64)          :: omega_p
    real(real64)          :: gamma_p
    real(real64)          :: strength_p
    real(real64), allocatable    :: current_p(:,:)
    real(real64), allocatable    :: e_field(:,:)
    real(real64), allocatable    :: e_field_dt_half(:,:)
    real(real64), allocatable    :: e_field_dt_full(:,:)
    real(real64), allocatable    :: current_at_point(:,:)
    real(real64), allocatable    :: selected_points_coordinate(:,:)
    integer               :: n_output_points
    integer               :: medium_type
    type(grid_t)          :: gr
    type(multicomm_t)     :: mc
    type(c_ptr)           :: write_handle
    type(output_t)        :: outp
    logical               :: from_scratch = .true.
    type(restart_t)       :: restart_load
    type(restart_t)       :: restart_dump
 
  contains
    procedure :: init_interaction => dispersive_medium_init_interaction
    procedure :: init_interaction_as_partner => dispersive_medium_init_interaction_as_partner
    procedure :: initialize => dispersive_medium_initialize
    procedure :: do_algorithmic_operation => dispersive_medium_do_algorithmic_operation
    procedure :: is_tolerance_reached => dispersive_medium_is_tolerance_reached
    procedure :: copy_quantities_to_interaction => dispersive_medium_copy_quantities_to_interaction
    procedure :: restart_write_data => dispersive_medium_restart_write_data
    procedure :: restart_read_data => dispersive_medium_restart_read_data
    procedure :: update_kinetic_energy => dispersive_medium_update_kinetic_energy
    procedure :: output_start => dispersive_medium_output_start
    procedure :: output_write => dispersive_medium_output_write
    procedure :: output_finish => dispersive_medium_output_finish
    procedure :: init_parallelization => dispersive_medium_init_parallelization
    procedure :: get_efield => dispersive_medium_get_efield
    final :: dispersive_medium_finalize
  end type dispersive_medium_t
 
  interface dispersive_medium_t
    procedure dispersive_medium_constructor
  end interface dispersive_medium_t
 
  integer, public, parameter ::      &
    DRUDE_MEDIUM                = 0
 
contains
 
  ! ---------------------------------------------------------
  function dispersive_medium_constructor(namespace) result(sys)
    class(dispersive_medium_t), pointer    :: sys
    type(namespace_t),          intent(in) :: namespace
 
    push_sub(dispersive_medium_constructor)
 
    allocate(sys)
 
    call dispersive_medium_init(sys, namespace)
 
    pop_sub(dispersive_medium_constructor)
  end function dispersive_medium_constructor
 
  ! ---------------------------------------------------------
  ! ---------------------------------------------------------
  subroutine dispersive_medium_init(this, namespace)
    class(dispersive_medium_t), target, intent(inout) :: this
    type(namespace_t),            intent(in)    :: namespace
 
    integer :: nlines, ncols, idim, il
    real(real64) :: pos(3)
    type(block_t) :: blk
 
    push_sub(dispersive_medium_init)
 
    this%namespace = namespace
 
    call profiling_in('DISP_MEDIUM_INIT')
 
    this%space = space_t(this%namespace)
    if (this%space%dim /= 3) then
      call messages_not_implemented('Linear medium for dimensions other than 3', namespace=namespace)
    end if
    call grid_init_stage_1(this%gr, this%namespace, this%space, latt=lattice_vectors_t(namespace, this%space))
 
    ! Parse electromagnetic properties of dispersive media...
 
    !%Variable MediumDispersionType
    !%Type integer
    !%Default drude_medium
    !%Section Maxwell::Medium
    !%Description
    !% Dispersion model used for the medium (only Drude model available for the moment).
    !%Option drude_medium 0
    !% Drude type of dispersion.
    !%End
    call parse_variable(namespace, 'MediumDispersionType', drude_medium, this%medium_type)
    if (.not. varinfo_valid_option('MediumDispersionType', this%medium_type)) then
      call messages_input_error(namespace, 'MediumDispersionType')
    end if
 
    ! Parse electromagnetic properties of Dispersive media...
    !%Variable MediumPoleEnergy
    !%Type float
    !%Default 0
    !%Section Maxwell::Medium
    !%Description
    !% Energy of the pole.
    !%End
    call parse_variable(namespace, 'MediumPoleEnergy', m_zero, this%omega_p, unit_one/units_inp%time)
 
    !%Variable MediumPoleDamping
    !%Type float
    !%Default 0
    !%Section Maxwell::Medium
    !%Description
    !% Damping factor (inverse relaxation time) of the medium.
    !%End
    call parse_variable(namespace, 'MediumPoleDamping', m_zero, this%gamma_p, unit_one/units_inp%time)
 
    !%Variable MediumPoleStrength
    !%Type float
    !%Default 1.0
    !%Section Maxwell::Medium
    !%Description
    !% Strength of the pole (unitless).
    !%End
    call parse_variable(namespace, 'MediumPoleStrength', m_one, this%strength_p, unit_one)
 
    !%Variable MediumCurrentCoordinates
    !%Type block
    !%Section Maxwell::Output
    !%Description
    !%  This allows to output phasor current vectors at particular points in space.
    !%
    !% <tt>%MediumCurrentCoordinates
    !% <br>&nbsp;&nbsp;    -1.0 | 2.0 |  4.0
    !% <br>&nbsp;&nbsp;     0.0 | 1.0 | -2.0
    !% <br>%</tt>
    !%
    !%End
 
    if (parse_block(namespace, 'MediumCurrentCoordinates', blk) == 0) then
      nlines = parse_block_n(blk)
      this%n_output_points = nlines
      safe_allocate(this%selected_points_coordinate(1:nlines,1:3))
      safe_allocate(this%current_at_point(1:nlines,1:3))
      do il = 1, nlines
        ncols = parse_block_cols(blk,0)
        do idim = 1, 3
          call parse_block_float(blk, il-1, idim-1, pos(idim), units_inp%length)
        end do
        this%selected_points_coordinate(il,:) = pos(:)
        this%current_at_point(il,:)  = m_zero
      end do
      call parse_block_end(blk)
    else
      this%n_output_points = 1
      safe_allocate(this%selected_points_coordinate(1,3))
      safe_allocate(this%current_at_point(1,3))
      this%selected_points_coordinate = m_zero
      this%current_at_point = m_zero
    end if
 
    this%supported_interactions = [mxll_e_field_to_matter]
    this%supported_interactions_as_partner = [current_to_mxll_field]
 
    call this%quantities%add(quantity_t("current", updated_on_demand = .false.))
 
    call profiling_out('DISP_MEDIUM_INIT')
 
    pop_sub(dispersive_medium_init)
  end subroutine dispersive_medium_init
 
  ! ---------------------------------------------------------
  subroutine dispersive_medium_init_parallelization(this, grp)
    class(dispersive_medium_t),     intent(inout) :: this
    type(mpi_grp_t),      intent(in)    :: grp
 
    integer(int64) :: index_range(4)
    integer :: ierr
 
    push_sub(dispersive_medium_init_parallelization)
 
    call system_init_parallelization(this, grp)
    ! store the ranges for these two indices (serves as initial guess
    ! for parallelization strategy)
    index_range(1) = this%gr%np_global  ! Number of points in mesh
    index_range(2) = 1                      ! Number of states
    index_range(3) = 1                      ! Number of k-points
    index_range(4) = 100000                 ! Some large number
 
    ! create index and domain communicators
    call multicomm_init(this%mc, this%namespace, mpi_world, calc_mode_par, mpi_world%size, &
      index_range, (/ 5000, 1, 1, 1 /))
    call grid_init_stage_2(this%gr, this%namespace, this%space, this%mc)
 
    call this%restart_dump%init(this%namespace, restart_td, restart_type_dump, &
      this%mc, ierr, mesh=this%gr)
    call this%restart_load%init(this%namespace, restart_td, restart_type_load, &
      this%mc, ierr, mesh=this%gr)
 
    safe_allocate(this%current_p(1:this%gr%np, 1:3))
    safe_allocate(this%e_field(1:this%gr%np, 1:3))
    this%e_field(:,:) = m_zero
 
    pop_sub(dispersive_medium_init_parallelization)
  end subroutine dispersive_medium_init_parallelization
 
  ! ---------------------------------------------------------
  subroutine dispersive_medium_init_interaction(this, interaction)
    class(dispersive_medium_t), target, intent(inout) :: this
    class(interaction_t),               intent(inout) :: interaction
 
    integer :: depth
 
    push_sub(dispersive_medium_init_interaction)
 
    select type (interaction)
    type is (mxll_e_field_to_matter_t)
      call interaction%init(this%gr, 3)
      interaction%type = mxll_field_total
 
      ! set interpolation depth for interaction
      ! interpolation depth depends on the propagator
      select type (prop => this%algo)
      type is (propagator_rk4_t)
        ! TODO: should be five, two is only for backwards compatibility
        depth = 2
      class default
        message(1) = "The chosen propagator does not yet support interaction interpolation"
        call messages_fatal(1)
      end select
      call interaction%init_interpolation(depth, interaction%label)
 
    class default
      message(1) = "Trying to initialize an unsupported interaction by a Dispersive medium."
      call messages_fatal(1)
    end select
 
    pop_sub(dispersive_medium_init_interaction)
  end subroutine dispersive_medium_init_interaction
 
  ! ---------------------------------------------------------
  subroutine dispersive_medium_init_interaction_as_partner(partner, interaction)
    class(dispersive_medium_t),       intent(in)    :: partner
    class(interaction_surrogate_t), intent(inout) :: interaction
 
    push_sub(dispersive_medium_init_interaction_as_partner)
 
    select type (interaction)
    type is (current_to_mxll_field_t)
      call interaction%init_from_partner(partner%gr, partner%space, partner%namespace)
    class default
      message(1) = "Trying to initialize an unsupported interaction by a linear medium."
      call messages_fatal(1)
    end select
 
    pop_sub(dispersive_medium_init_interaction_as_partner)
  end subroutine dispersive_medium_init_interaction_as_partner
 
  ! ---------------------------------------------------------
  subroutine dispersive_medium_initialize(this)
    class(dispersive_medium_t), intent(inout) :: this
 
    push_sub(dispersive_medium_initialize)
 
    this%from_scratch = .true.
    this%current_p(:,:) = m_zero
 
    pop_sub(dispersive_medium_initialize)
  end subroutine dispersive_medium_initialize
 
  ! ---------------------------------------------------------
  logical function dispersive_medium_do_algorithmic_operation(this, operation, updated_quantities) result(done)
    class(dispersive_medium_t),     intent(inout) :: this
    class(algorithmic_operation_t), intent(in)    :: operation
    character(len=:), allocatable,  intent(out)   :: updated_quantities(:)
 
    integer :: ip, idim
    real(real64) :: k1(1:3), k2(1:3), k3(1:3), k4(1:3)
 
    push_sub(dispersive_medium_do_algorithmic_operation)
    call profiling_in(trim(this%namespace%get())//":"//trim(operation%id))
 
    ! calculation of the current using ADE
    ! \partial_t J_P(t) = - \gamma_p * J_P(t) + \epsilon_0 \omega_p^2 E(t)
    ! (Computational Electrodynamics, Taflov and Hagness, 3rd Ed., section 9.4.3, eq. 9.56c)
    ! Analysis of Units:
    ! [e/time^2*a0^2] = (1/time) * (e/time*a0^2) + (e^2/hbar*c) * (1/time)^2 * Ha / (e * a0)
    ! [e/time^2*a0^2] = (e/time^2*a0^2) + (e/hbar*c) * 1/time^2 * Ha/a0, and [c]=a0/time, so [hbar*c]=Ha*a0
    ! [e/time^2*a0^2] = (e/time^2*a0^2) + (e/time^2*a0^2)
 
    select type (algo => this%algo)
    class is (propagator_t)
 
      done = .true.
      select case (operation%id)
      case (store_current_status)
        ! For the moment we do nothing
 
      case (rk4_start)
        safe_allocate(this%e_field_dt_half(1:this%gr%np, 1:3))
        safe_allocate(this%e_field_dt_full(1:this%gr%np, 1:3))
 
      case (rk4_finish)
        safe_deallocate_a(this%e_field_dt_half)
        safe_deallocate_a(this%e_field_dt_full)
 
      case (rk4_extrapolate)
        call this%get_efield(this%iteration%value(), this%e_field)
        call this%get_efield(this%iteration%value()+algo%dt/m_two, this%e_field_dt_half)
        call this%get_efield(this%iteration%value()+algo%dt, this%e_field_dt_full)
 
      case (rk4_propagate)
        do idim=1, 3
          !$omp parallel do private(k1, k2, k3, k4)
          do ip = 1, this%gr%np
            k1(idim) = current_derivative_dir(this%current_p(ip, idim), &
              this%e_field(ip,idim), this%gamma_p, this%omega_p, this%strength_p)
 
            k2(idim) = current_derivative_dir(this%current_p(ip, idim) + algo%dt * k1(idim) / m_two, &
              this%e_field_dt_half(ip, idim), this%gamma_p, this%omega_p, this%strength_p)
 
            k3(idim) = current_derivative_dir(this%current_p(ip, idim) + algo%dt * k2(idim) / m_two, &
              this%e_field_dt_half(ip, idim), this%gamma_p, this%omega_p, this%strength_p)
 
            k4(idim) = current_derivative_dir(this%current_p(ip, idim) + algo%dt * k3(idim), &
              this%e_field_dt_full(ip, idim), this%gamma_p, this%omega_p, this%strength_p)
 
            this%current_p(ip, idim) = this%current_p(ip, idim) + algo%dt / 6.0_real64 * &
              (k1(idim) + m_two * (k2(idim) + k3(idim)) + k4(idim))
          end do
        end do
        updated_quantities = ["current"]
 
      case default
        done = .false.
      end select
    end select
 
    call profiling_out(trim(this%namespace%get())//":"//trim(operation%id))
    pop_sub(dispersive_medium_do_algorithmic_operation)
  contains
 
    function current_derivative_dir(current_p, e_field, gamma_p, omega_p, strength_p) result(current_dot)
      real(real64), intent(in)      :: current_p
      real(real64), intent(in)      :: e_field
      real(real64), intent(in)      :: gamma_p
      real(real64), intent(in)      :: omega_p
      real(real64), intent(in)      :: strength_p
      real(real64)           :: current_dot
 
      current_dot = - gamma_p * current_p + strength_p * p_ep * omega_p**2 * e_field
    end function current_derivative_dir
 
  end function dispersive_medium_do_algorithmic_operation
 
  ! ---------------------------------------------------------
  logical function dispersive_medium_is_tolerance_reached(this, tol) result(converged)
    class(dispersive_medium_t),   intent(in)    :: this
    real(real64),              intent(in)    :: tol
 
    push_sub(dispersive_medium_is_tolerance_reached)
 
    ! this routine is never called at present, no reason to be here
    assert(.false.)
    converged = .false.
 
    pop_sub(dispersive_medium_is_tolerance_reached)
  end function dispersive_medium_is_tolerance_reached
 
  ! ---------------------------------------------------------
  subroutine dispersive_medium_copy_quantities_to_interaction(partner, interaction)
    class(dispersive_medium_t),          intent(inout) :: partner
    class(interaction_surrogate_t),       intent(inout) :: interaction
 
    push_sub(dispersive_medium_copy_quantities_to_interaction)
    call profiling_in(trim(partner%namespace%get())//":"//"COPY_QUANTITY_INTER")
 
    select type (interaction)
    type is (current_to_mxll_field_t)
      interaction%partner_field(:,:) = partner%current_p
      call interaction%do_mapping()
    class default
      message(1) = "Unsupported interaction."
      call messages_fatal(1)
    end select
 
    call profiling_out(trim(partner%namespace%get())//":"//"COPY_QUANTITY_INTER")
    pop_sub(dispersive_medium_copy_quantities_to_interaction)
  end subroutine dispersive_medium_copy_quantities_to_interaction
 
  ! ---------------------------------------------------------
  subroutine dispersive_medium_restart_write_data(this)
    class(dispersive_medium_t), intent(inout) :: this
 
    character(len=256) :: filename
    integer :: idir, err
    type(interaction_iterator_t) :: iter
    class(interaction_t), pointer :: interaction
 
    push_sub(dispersive_medium_restart_write_data)
    call profiling_in(trim(this%namespace%get())//":"//"RESTART_WRITE")
 
    if (.not. this%restart_dump%skip()) then
      do idir = 1, this%space%dim
        write(filename, "(A,A)") "current_p-", index2axis(idir)
        call this%restart_dump%write_mesh_function(this%space, trim(filename), this%gr, &
          this%current_p(:, idir), err)
      end do
 
      call iter%start(this%interactions)
      do while (iter%has_next())
        interaction => iter%get_next()
        select type (interaction)
        class is (mxll_e_field_to_matter_t)
          call interaction%write_restart(this%gr, this%space, this%restart_dump, err)
        end select
      end do
 
      if (err == 0) then
        message(1) = "Successfully wrote restart data for system "//trim(this%namespace%get())
        call messages_info(1, namespace=this%namespace)
      end if
 
    end if
 
    call profiling_out(trim(this%namespace%get())//":"//"RESTART_WRITE")
    pop_sub(dispersive_medium_restart_write_data)
  end subroutine dispersive_medium_restart_write_data
 
  ! ---------------------------------------------------------
  ! this function returns true if restart data could be read
  logical function dispersive_medium_restart_read_data(this)
    class(dispersive_medium_t), intent(inout) :: this
 
    character(len=256) :: filename
    integer :: idir, err
    type(interaction_iterator_t) :: iter
    class(interaction_t), pointer :: interaction
 
    push_sub(dispersive_medium_restart_read_data)
    call profiling_in(trim(this%namespace%get())//":"//"RESTART_READ")
 
    dispersive_medium_restart_read_data = .false.
    if (.not. this%restart_load%skip()) then
      do idir = 1, 3
        write(filename, "(A,A)") "current_p-", index2axis(idir)
        call this%restart_load%read_mesh_function(this%space, trim(filename), this%gr, &
          this%current_p(:, idir), err)
      end do
 
      call iter%start(this%interactions)
      do while (iter%has_next())
        interaction => iter%get_next()
        select type (interaction)
        class is (mxll_e_field_to_matter_t)
          call interaction%read_restart(this%gr, this%space, this%restart_load, err)
        end select
      end do
 
      if (err == 0) then
        dispersive_medium_restart_read_data = .true.
        this%from_scratch = .false.
      else
        ! set to 0 again in case this was read incompletely
        this%current_p(:, :) = m_zero
      end if
    end if
 
    call profiling_out(trim(this%namespace%get())//":"//"RESTART_READ")
    pop_sub(dispersive_medium_restart_read_data)
  end function dispersive_medium_restart_read_data
 
  ! ---------------------------------------------------------
  subroutine dispersive_medium_update_kinetic_energy(this)
    class(dispersive_medium_t), intent(inout) :: this
 
    push_sub(dispersive_medium_update_kinetic_energy)
 
    ! TODO: evaluate proper energy associated with the current distribution
    ! For Drude model: check Giuliani/Vignale book, section 4.6.1
    this%kinetic_energy = m_zero
 
    pop_sub(dispersive_medium_update_kinetic_energy)
 
  end subroutine dispersive_medium_update_kinetic_energy
 
  ! ---------------------------------------------------------
  subroutine dispersive_medium_output_start(this)
    class(dispersive_medium_t), intent(inout) :: this
 
    integer :: first, id, idir
    character(len=130) :: aux
 
    push_sub(dispersive_medium_output_start)
 
    call profiling_in("DISP_MEDIUM_OUTPUT_START")
 
    select type (algo => this%algo)
    class is (propagator_t)
 
      if (this%iteration%counter() == 0) then
        first = 0
      else
        first = this%iteration%counter() + 1
      end if
 
      call io_mkdir('td.general', this%namespace)
      call write_iter_init(this%write_handle, first, units_from_atomic(units_out%time, algo%dt), &
        trim(io_workpath("td.general/current_at_points", this%namespace)))
 
      if (mpi_world%is_root()) then
        if (this%iteration%counter() == 0) then
          call write_iter_clear(this%write_handle)
          call write_iter_string(this%write_handle,&
            '################################################################################')
          call write_iter_nl(this%write_handle)
          call write_iter_string(this%write_handle,'# HEADER')
          call write_iter_nl(this%write_handle)
 
          ! first line
          write(aux, '(a7,e20.12,3a)') '# dt = ', units_from_atomic(units_out%time, algo%dt), &
            " [", trim(units_abbrev(units_out%time)), "]"
          call write_iter_string(this%write_handle, aux)
          call write_iter_nl(this%write_handle)
 
          call write_iter_header_start(this%write_handle)
 
          do id = 1, this%n_output_points
            do idir = 1, 3
              write(aux, '(a,i1,a,i1,a)') 'j(', id, ',', idir, ')'
              call write_iter_header(this%write_handle, aux)
            end do
          end do
 
          call write_iter_nl(this%write_handle)
          call write_iter_header(this%write_handle, '#          [' // trim(units_abbrev(units_out%time)) // ']')
 
          !FIXME: this is not printing the proper unit to output yet, for some reason
          aux = '          [' // trim(units_abbrev(unit_one/(units_out%time*units_out%length**2))) // ']'
          do id = 1, this%n_output_points
            do idir = 1, 3
              call write_iter_header(this%write_handle, aux)
            end do
          end do
          call write_iter_nl(this%write_handle)
          call write_iter_string(this%write_handle,&
            '################################################################################')
          call write_iter_nl(this%write_handle)
        end if
      end if
 
      if (first == 0) call this%output_write()
 
    end select
 
    call profiling_out("DISP_MEDIUM_OUTPUT_START")
 
    pop_sub(dispersive_medium_output_start)
  end subroutine dispersive_medium_output_start
 
  ! ---------------------------------------------------------
  subroutine dispersive_medium_output_write(this)
    class(dispersive_medium_t), intent(inout) :: this
 
    real(real64)   :: dmin, dtmp(3)
    integer :: ip, pos_index, rankmin
 
    push_sub(dispersive_medium_output_write)
    call profiling_in(trim(this%namespace%get())//":"//"OUTPUT_WRITE")
 
    select type (algo => this%algo)
    class is (propagator_t)
      do ip = 1, this%n_output_points
        pos_index = mesh_nearest_point(this%gr, this%selected_points_coordinate(ip,:), dmin, rankmin)
        if (this%gr%mpi_grp%rank == rankmin) then
          dtmp(:) = this%current_p(pos_index,:)
        end if
        if (this%gr%parallel_in_domains) then
          call this%gr%mpi_grp%bcast(dtmp(:), 3, mpi_double_precision, rankmin)
        end if
        this%current_at_point(ip,:) = units_from_atomic((unit_one/units_out%time)/(units_out%length**2), dtmp(:))
      end do
 
      if (.not. mpi_world%is_root()) then
        call profiling_out(trim(this%namespace%get())//":"//"OUTPUT_WRITE")
        pop_sub(dispersive_medium_output_write)
        return ! only first node outputs
      end if
 
      call write_iter_start(this%write_handle)
      do ip = 1, this%n_output_points
        dtmp = this%current_at_point(ip,1:3)
        call write_iter_double(this%write_handle, dtmp, 3)
      end do
 
      call write_iter_nl(this%write_handle)
      call write_iter_flush(this%write_handle)
 
    end select
 
    call profiling_out(trim(this%namespace%get())//":"//"OUTPUT_WRITE")
    pop_sub(dispersive_medium_output_write)
  end subroutine dispersive_medium_output_write
 
  ! ---------------------------------------------------------
  subroutine dispersive_medium_output_finish(this)
    class(dispersive_medium_t), intent(inout) :: this
 
 
    push_sub(dispersive_medium_output_finish)
 
    call profiling_in("DISP_MEDIUM_OUTPUT_FINISH")
 
    call write_iter_end(this%write_handle)
 
    call profiling_out("DISP_MEDIUM_OUTPUT_FINISH")
 
    pop_sub(dispersive_medium_output_finish)
  end subroutine dispersive_medium_output_finish
 
  ! ---------------------------------------------------------
  subroutine dispersive_medium_get_efield(this, time, efield)
    class(dispersive_medium_t), intent(inout) :: this
    real(real64),               intent(in)    :: time
    real(real64), contiguous,   intent(inout) :: efield(:, :)
 
    type(interaction_iterator_t) :: iter
    real(real64), allocatable :: efield_tmp(:, :)
 
    push_sub(dispersive_medium_get_efield)
 
    safe_allocate(efield_tmp(1:this%gr%np, 1:3))
    efield = m_z0
    ! interpolate efield from interaction
    call iter%start(this%interactions)
    do while (iter%has_next())
      select type (interaction => iter%get_next())
      class is (mxll_e_field_to_matter_t)
        call interaction%interpolate(time, efield_tmp)
        call lalg_axpy(this%gr%np, 3, m_one, efield_tmp, efield)
      end select
    end do
    safe_deallocate_a(efield_tmp)
 
    pop_sub(dispersive_medium_get_efield)
  end subroutine dispersive_medium_get_efield
 
  ! ---------------------------------------------------------
  subroutine dispersive_medium_finalize(this)
    type(dispersive_medium_t), intent(inout) :: this
 
    push_sub(dispersive_medium_finalize)
    call system_end(this)
    safe_deallocate_a(this%current_p)
    safe_deallocate_a(this%e_field)
    safe_deallocate_a(this%selected_points_coordinate)
    safe_deallocate_a(this%current_at_point)
    call multicomm_end(this%mc)
    call grid_end(this%gr)
    pop_sub(dispersive_medium_finalize)
  end subroutine dispersive_medium_finalize
 
end module dispersive_medium_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: