maxwell.F90

Go to the documentation of this file.

!! Copyright (C) 2019-2020 Franco Bonafe, Heiko Appel, Rene Jestaedt
!!
!! 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 maxwell_oct_m
  use accel_oct_m
  use batch_oct_m
  use batch_ops_oct_m
  use algorithm_oct_m
  use algorithm_factory_oct_m
  use calc_mode_par_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use distributed_oct_m
  use external_densities_oct_m
  use exponential_oct_m
  use field_transfer_oct_m
  use time_interpolation_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_mxll_oct_m
  use helmholtz_decomposition_m
  use index_oct_m
  use interaction_oct_m
  use interaction_surrogate_oct_m
  use interaction_oct_m
  use interaction_enum_oct_m
  use io_oct_m
  use lalg_basic_oct_m
  use lattice_vectors_oct_m
  use linear_medium_to_em_field_oct_m
  use current_to_mxll_field_oct_m
  use loct_oct_m
  use lorentz_force_oct_m
  use math_oct_m
  use maxwell_boundary_op_oct_m
  use mesh_oct_m
  use messages_oct_m
  use mesh_interpolation_oct_m
  use mpi_oct_m
  use multicomm_oct_m
  use mxll_e_field_to_matter_oct_m
  use mxll_b_field_to_matter_oct_m
  use mxll_vec_pot_to_matter_oct_m
  use namespace_oct_m
  use output_oct_m
  use output_low_oct_m
  use output_mxll_oct_m
  use parser_oct_m
  use poisson_oct_m
  use profiling_oct_m
  use propagator_exp_mid_oct_m
  use propagator_exp_gauss1_oct_m
  use propagator_exp_gauss2_oct_m
  use propagator_oct_m
  use propagator_mxll_oct_m
  use propagator_leapfrog_oct_m
  use quantity_oct_m
  use regridding_oct_m
  use restart_oct_m
  use sort_oct_m
  use space_oct_m
  use system_oct_m
  use states_mxll_oct_m
  use states_mxll_restart_oct_m
  use electrons_oct_m
  use td_write_oct_m
  use types_oct_m
  use unit_oct_m
  use unit_system_oct_m
 
 
  implicit none
 
  private
  public ::               &
    maxwell_t,        &
    maxwell_init
 
  integer, parameter, public ::           &
    MULTIGRID_MX_TO_MA_EQUAL   = 1,       &
    multigrid_mx_to_ma_large   = 2
 
  type, extends(system_t) :: maxwell_t
    type(space_t)                :: space
    type(states_mxll_t)          :: st
    type(hamiltonian_mxll_t)     :: hm
    type(grid_t)                 :: gr
    type(output_t)               :: outp
    type(multicomm_t)            :: mc
 
    type(mesh_interpolation_t)      :: mesh_interpolate
 
    type(propagator_mxll_t)         :: tr_mxll
    type(td_write_t)                :: write_handler
    type(c_ptr)                     :: output_handle
 
    complex(real64), allocatable           :: ff_rs_inhom_t1(:,:), ff_rs_inhom_t2(:,:)
    complex(real64), allocatable           :: rs_state_init(:,:)
    type(time_interpolation_t), pointer :: current_interpolation
    real(real64)                 :: bc_bounds(2, 3), dt_bounds(2, 3)
    integer                      :: energy_update_iter
    type(restart_t)              :: restart_dump
    type(restart_t)              :: restart
 
    type(lattice_vectors_t)         :: latt
 
    type(helmholtz_decomposition_t) :: helmholtz
 
    logical                         :: write_previous_state = .false.
 
  contains
    procedure :: init_interaction => maxwell_init_interaction
    procedure :: init_parallelization => maxwell_init_parallelization
    procedure :: initialize => maxwell_initialize
    procedure :: do_algorithmic_operation => maxwell_do_algorithmic_operation
    procedure :: is_tolerance_reached => maxwell_is_tolerance_reached
    procedure :: init_interaction_as_partner => maxwell_init_interaction_as_partner
    procedure :: copy_quantities_to_interaction => maxwell_copy_quantities_to_interaction
    procedure :: output_start => maxwell_output_start
    procedure :: output_write => maxwell_output_write
    procedure :: output_finish => maxwell_output_finish
    procedure :: update_interactions_start => maxwell_update_interactions_start
    procedure :: update_interactions_finish => maxwell_update_interactions_finish
    procedure :: restart_write_data => maxwell_restart_write_data
    procedure :: restart_read_data => maxwell_restart_read_data
    procedure :: update_kinetic_energy => maxwell_update_kinetic_energy
    procedure :: get_current => maxwell_get_current
    final :: maxwell_finalize
  end type maxwell_t
 
  interface maxwell_t
    procedure maxwell_constructor
  end interface maxwell_t
 
contains
 
  ! ---------------------------------------------------------
  function maxwell_constructor(namespace) result(sys)
    class(maxwell_t),   pointer    :: sys
    type(namespace_t),  intent(in) :: namespace
 
    push_sub(maxwell_constructor)
 
    allocate(sys)
 
    call maxwell_init(sys, namespace)
 
    pop_sub(maxwell_constructor)
  end function maxwell_constructor
 
 
  ! ---------------------------------------------------------
  subroutine maxwell_init(this, namespace)
    class(maxwell_t),     intent(inout) :: this
    type(namespace_t),    intent(in)    :: namespace
 
 
    push_sub(maxwell_init)
 
    call profiling_in("MAXWELL_INIT")
 
    this%namespace = namespace
 
    call messages_obsolete_variable(this%namespace, 'SystemName')
    call messages_print_with_emphasis(msg='Maxwell Simulation Start', namespace=this%namespace)
    this%space = space_t(this%namespace)
 
    call grid_init_stage_1(this%gr, this%namespace, this%space, latt=lattice_vectors_t(namespace, this%space))
    call states_mxll_init(this%st, this%namespace, this%space)
    this%latt = lattice_vectors_t(this%namespace, this%space)
 
    call this%quantities%add(quantity_t("E field", updated_on_demand = .false.))
    call this%quantities%add(quantity_t("B field", updated_on_demand = .false.))
    call this%quantities%add(quantity_t("vector potential", updated_on_demand = .false.))
 
    call mesh_interpolation_init(this%mesh_interpolate, this%gr)
 
    this%supported_interactions = [linear_medium_to_em_field, current_to_mxll_field]
    this%supported_interactions_as_partner = [lorentz_force, mxll_e_field_to_matter, mxll_b_field_to_matter, mxll_vec_pot_to_matter]
 
    call profiling_out("MAXWELL_INIT")
 
    pop_sub(maxwell_init)
  end subroutine maxwell_init
 
  ! ---------------------------------------------------------
  subroutine maxwell_init_interaction(this, interaction)
    class(maxwell_t),     target, intent(inout) :: this
    class(interaction_t),         intent(inout) :: interaction
 
    integer :: depth
 
    push_sub(maxwell_init_interaction)
 
    select type (interaction)
    type is (linear_medium_to_em_field_t)
      call interaction%init(this%gr)
    type is (current_to_mxll_field_t)
      call interaction%init(this%gr, this%st%dim)
      call interaction%init_space_latt(this%space, this%latt)
 
      ! set interpolation depth for interaction
      ! interpolation depth depends on the propagator
      select type (prop => this%algo)
      type is (propagator_exp_mid_t)
        depth = 2
      type is (propagator_leapfrog_t)
        depth = 2
      type is (propagator_exp_gauss1_t)
        depth = 2
      type is (propagator_exp_gauss2_t)
        depth = 4
      class default
        message(1) = "The chosen propagator does not yet support interaction interpolation"
        call messages_fatal(1, namespace=this%namespace)
      end select
      call interaction%init_interpolation(depth, interaction%label, cmplx=.true.)
 
      this%hm%current_density_from_medium = .true.
    class default
      message(1) = "Trying to initialize an unsupported interaction by Maxwell."
      call messages_fatal(1, namespace=this%namespace)
    end select
 
    pop_sub(maxwell_init_interaction)
  end subroutine maxwell_init_interaction
 
  ! ---------------------------------------------------------
  subroutine maxwell_init_parallelization(this, grp)
    class(maxwell_t),     intent(inout) :: this
    type(mpi_grp_t),      intent(in)    :: grp
 
    integer(int64) :: index_range(4)
    integer :: ierr, ip, pos_index, rankmin
    real(real64) :: dmin
 
    push_sub(maxwell_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) = this%st%nst        ! 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 mpi_grp_init(this%st%system_grp, this%mc%master_comm)
 
    call grid_init_stage_2(this%gr, this%namespace, this%space, this%mc)
    call output_mxll_init(this%outp, this%namespace, this%space)
    call hamiltonian_mxll_init(this%hm, this%namespace, this%gr, this%st)
 
    this%st%energy_rate = m_zero
    this%st%delta_energy = m_zero
    this%st%energy_via_flux_calc = m_zero
    this%st%trans_energy_rate = m_zero
    this%st%trans_delta_energy = m_zero
    this%st%trans_energy_via_flux_calc = m_zero
    this%st%plane_waves_energy_rate = m_zero
    this%st%plane_waves_delta_energy = m_zero
    this%st%plane_waves_energy_via_flux_calc = m_zero
 
    safe_allocate(this%rs_state_init(1:this%gr%np_part, 1:this%st%dim))
    this%rs_state_init(:,:) = m_z0
 
    this%energy_update_iter = 1
 
    call poisson_init(this%st%poisson, this%namespace, this%space, this%gr%der, this%mc, this%gr%stencil)
 
    call bc_mxll_init(this%hm%bc, this%namespace, this%space, this%gr, this%st)
    this%bc_bounds(:,1:3) = this%hm%bc%bc_bounds(:,1:3)
    call inner_and_outer_points_mapping(this%gr, this%st, this%bc_bounds)
    this%dt_bounds(2, 1:3) = this%bc_bounds(1, 1:3)
    this%dt_bounds(1, 1:3) = this%bc_bounds(1, 1:3) - this%gr%der%order * this%gr%spacing(1:3)
    call surface_grid_points_mapping(this%gr, this%st, this%dt_bounds)
 
    call propagator_mxll_init(this%gr, this%namespace, this%st, this%hm, this%tr_mxll)
    call states_mxll_allocate(this%st, this%gr)
    call external_current_init(this%st, this%space, this%namespace, this%gr)
    this%hm%propagation_apply = .true.
 
    ! set map for selected points
    do ip = 1, this%st%selected_points_number
      pos_index = mesh_nearest_point(this%gr, this%st%selected_points_coordinate(:,ip), dmin, rankmin)
      if (this%gr%mpi_grp%rank == rankmin) then
        this%st%selected_points_map(ip) = pos_index
      else
        this%st%selected_points_map(ip) = -1
      end if
    end do
    if (accel_is_enabled()) then
      call accel_create_buffer(this%st%buff_selected_points_map, accel_mem_read_only, type_integer, &
        this%st%selected_points_number)
      call accel_write_buffer(this%st%buff_selected_points_map, this%st%selected_points_number, &
        this%st%selected_points_map)
    end if
 
    this%hm%plane_waves_apply = .true.
    this%hm%spatial_constant_apply = .true.
 
    call this%restart%init(this%namespace, restart_td, restart_type_load, this%mc, ierr, mesh=this%gr)
    call this%restart_dump%init(this%namespace, restart_td, restart_type_dump, this%mc, ierr, mesh=this%gr)
 
    ! initialize batches
    call zbatch_init(this%st%rs_stateb, 1, 1, this%st%dim, this%gr%np_part)
    if (this%st%pack_states) call this%st%rs_stateb%do_pack()
    call this%st%rs_stateb%copy_to(this%st%rs_state_prevb)
    call this%st%rs_stateb%copy_to(this%st%inhomogeneousb)
    if (this%tr_mxll%bc_plane_waves) then
      call this%st%rs_stateb%copy_to(this%st%rs_state_plane_wavesb)
    end if
 
    ! the order should depend on the propagator
    !this%current_interpolation => time_interpolation_t(this%gr%np, this%st%dim, 2, .true., "current")
    ! Initialize Helmholtz decomposition
    call this%helmholtz%init(this%namespace, this%gr, this%mc, this%space)
 
    pop_sub(maxwell_init_parallelization)
  end subroutine maxwell_init_parallelization
 
  ! ---------------------------------------------------------
  subroutine maxwell_initialize(this)
    class(maxwell_t), intent(inout) :: this
 
    real(real64)   :: courant
 
 
    push_sub(maxwell_initialize)
 
    call profiling_in("MAXWELL_INIT_CONDITIONS")
 
    select type (algo => this%algo)
    class is (propagator_t)
 
      ! The reference Courant–Friedrichs–Lewy condition for stability here
      ! is calculated with respect to R. Courant, K. Friedrichs, and H. Lewy
      ! On the Partial Difference Equations of
      ! Mathematical Physics, IBM Journal of Research and Development, vol. 11
      ! pp. 215-234, Mar. 1967.
      courant = m_one/(p_c * sqrt(m_one/this%gr%spacing(1)**2 + m_one/this%gr%spacing(2)**2 + &
        m_one/this%gr%spacing(3)**2))
 
      if (algo%dt > (m_one + m_fourth) * courant) then
        write(message(1),'(a)') 'The timestep is too large compared to the Courant-Friedrichs-Lewy'
        write(message(2),'(a)') 'stability condition. Time propagation might not be stable.'
        call messages_warning(2, namespace=this%namespace)
      end if
 
      if (parse_is_defined(this%namespace, 'UserDefinedInitialMaxwellStates')) then
        call states_mxll_read_user_def(this%namespace, this%space, this%gr, this%gr%der, this%st, this%hm%bc,&
          this%rs_state_init)
        call messages_print_with_emphasis(msg="Setting initial EM field inside box", namespace=this%namespace)
        ! TODO: add consistency check that initial state fulfills Gauss laws
        this%st%rs_state(:,:) = this%st%rs_state + this%rs_state_init
        if (this%tr_mxll%bc_plane_waves) then
          this%st%rs_state_plane_waves(:,:) = this%rs_state_init
        end if
      end if
 
      ! initialize the spatial constant field according to the conditions set in the
      ! UserDefinedConstantSpatialMaxwellField block
      if (this%tr_mxll%bc_constant) then
        call spatial_constant_calculation(this%tr_mxll%bc_constant, this%st, this%gr, this%hm, m_zero, &
          algo%dt, m_zero, this%st%rs_state, set_initial_state = .true.)
        ! for mesh parallelization, this needs communication!
        this%st%rs_state_const(:) = this%st%rs_state(mesh_global_index_from_coords(this%gr, [0,0,0]),:)
      end if
 
      if (parse_is_defined(this%namespace, 'UserDefinedInitialMaxwellStates')) then
        safe_deallocate_a(this%rs_state_init)
      end if
 
      call hamiltonian_mxll_update(this%hm, time = m_zero)
 
      ! calculate Maxwell energy
      call energy_mxll_calc(this%gr, this%st, this%hm, this%hm%energy, this%st%rs_state, &
        this%st%rs_state_plane_waves)
 
      ! Get RS states values for selected points
      call get_rs_state_at_point(this%st%selected_points_rs_state(:,:), this%st%rs_state, &
        this%st%selected_points_coordinate(:,:), this%st, this%gr)
 
      call mxll_set_batch(this%st%rs_stateb, this%st%rs_state, this%gr%np, this%st%dim)
      call batch_set_zero(this%st%inhomogeneousb)
      if (this%tr_mxll%bc_plane_waves) then
        call mxll_set_batch(this%st%rs_state_plane_wavesb, this%st%rs_state_plane_waves, this%gr%np, this%st%dim)
      end if
    end select
 
    call profiling_out("MAXWELL_INIT_CONDITIONS")
 
    pop_sub(maxwell_initialize)
  end subroutine maxwell_initialize
 
  ! ---------------------------------------------------------
  logical function maxwell_do_algorithmic_operation(this, operation, updated_quantities) result(done)
    class(maxwell_t),               intent(inout) :: this
    class(algorithmic_operation_t), intent(in)    :: operation
    character(len=:), allocatable,  intent(out)   :: updated_quantities(:)
 
    complex(real64), allocatable :: charge_density_ext(:)
    real(real64) :: current_time
 
    push_sub(maxwell_do_algorithmic_operation)
    call profiling_in(trim(this%namespace%get())//":"//trim(operation%id))
 
    current_time = this%iteration%value()
 
    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 (expmid_start)
        safe_allocate(this%ff_rs_inhom_t1(1:this%gr%np_part, 1:this%hm%dim))
        safe_allocate(this%ff_rs_inhom_t2(1:this%gr%np_part, 1:this%hm%dim))
 
      case (expmid_finish)
        safe_deallocate_a(this%ff_rs_inhom_t1)
        safe_deallocate_a(this%ff_rs_inhom_t2)
 
      case (expmid_extrapolate)
        if (this%hm%current_density_ext_flag .or. this%hm%current_density_from_medium) then
          call this%get_current(current_time, this%st%rs_current_density_t1)
          call this%get_current(current_time+algo%dt, this%st%rs_current_density_t2)
        end if
 
      case (expmid_propagate)
        ! Propagation
 
        !We first compute three external charge and current densities and we convert them as RS vectors
        safe_allocate(charge_density_ext(1:this%gr%np))
 
        !No charge density at the moment
        charge_density_ext = m_z0
 
        call transform_rs_densities(this%hm, this%gr, &
          this%st%rs_current_density_t1, this%ff_rs_inhom_t1, rs_trans_forward)
        call transform_rs_densities(this%hm, this%gr, &
          this%st%rs_current_density_t2, this%ff_rs_inhom_t2, rs_trans_forward)
 
        safe_deallocate_a(charge_density_ext)
 
        ! Propagation dt with H_maxwell
        call mxll_propagation_step(this%hm, this%namespace, this%gr, this%space, this%st, this%tr_mxll,&
          this%st%rs_stateb, this%ff_rs_inhom_t1, this%ff_rs_inhom_t2, current_time, algo%dt)
 
        updated_quantities = [character(16) :: "E field", "B field", "vector potential"]
 
      case (leapfrog_start)
        if (any(this%hm%bc%bc_ab_type == option__maxwellabsorbingboundaries__cpml)) then
          call bc_mxll_initialize_pml_simple(this%hm%bc, this%space, this%gr, this%hm%c_factor, algo%dt)
        end if
 
      case (leapfrog_finish)
 
      case (leapfrog_propagate)
        if (this%hm%current_density_ext_flag .or. this%hm%current_density_from_medium) then
          call this%get_current(current_time, this%st%rs_current_density_t1)
          call mxll_set_batch(this%st%inhomogeneousb, this%st%rs_current_density_t1, this%gr%np, this%st%dim)
          call batch_scal(this%gr%np, -m_one, this%st%inhomogeneousb)
        end if
 
        call mxll_propagate_leapfrog(this%hm, this%namespace, this%gr, this%st, this%tr_mxll, &
          current_time, algo%dt, this%iteration%counter())
 
        updated_quantities = [character(16) :: "E field", "B field", "vector potential"]
 
      case (exp_gauss1_start)
        safe_allocate(this%ff_rs_inhom_t1(1:this%gr%np_part, 1:this%hm%dim))
        if (any(this%hm%bc%bc_ab_type == option__maxwellabsorbingboundaries__cpml)) then
          call bc_mxll_initialize_pml_simple(this%hm%bc, this%space, this%gr, this%hm%c_factor, algo%dt)
        end if
 
      case (exp_gauss1_finish)
        safe_deallocate_a(this%ff_rs_inhom_t1)
 
      case (exp_gauss1_extrapolate)
        if (this%hm%current_density_ext_flag .or. this%hm%current_density_from_medium) then
          call this%get_current(current_time+algo%dt*m_half, this%st%rs_current_density_t1)
        end if
 
      case (exp_gauss1_propagate)
        ! the propagation step is
        ! F_{n+1} = F_n + dt * phi_1(-i*dt*H) [-i*H F_n - J_1]
        ! where J_1 = J(t_n + dt/2)
        call mxll_propagate_expgauss1(this%hm, this%namespace, this%gr, this%st, this%tr_mxll, &
          current_time, algo%dt)
 
        updated_quantities = [character(16) :: "E field", "B field", "vector potential"]
 
      case (exp_gauss2_start)
        safe_allocate(this%ff_rs_inhom_t1(1:this%gr%np_part, 1:this%hm%dim))
        safe_allocate(this%ff_rs_inhom_t2(1:this%gr%np_part, 1:this%hm%dim))
        if (any(this%hm%bc%bc_ab_type == option__maxwellabsorbingboundaries__cpml)) then
          call bc_mxll_initialize_pml_simple(this%hm%bc, this%space, this%gr, this%hm%c_factor, algo%dt)
        end if
 
      case (exp_gauss2_finish)
        safe_deallocate_a(this%ff_rs_inhom_t1)
        safe_deallocate_a(this%ff_rs_inhom_t2)
 
      case (exp_gauss2_extrapolate)
        if (this%hm%current_density_ext_flag .or. this%hm%current_density_from_medium) then
          call this%get_current(current_time+algo%dt*(m_half-sqrt(m_three)/6.0_real64), this%st%rs_current_density_t1)
          call this%get_current(current_time+algo%dt*(m_half+sqrt(m_three)/6.0_real64), this%st%rs_current_density_t2)
        end if
 
      case (exp_gauss2_propagate)
        ! the propagation step is
        ! F_{n+1} = F_n + dt * phi_1(-i*dt*H) [-i*H F_n - a_1 J_1 - a_2 J_2]
        !               + dt * phi_2(-i*dt*H) [-b_1 J_1 - b_2 J_2]
        ! where J_1 = J(t_n + (1/2 - sqrt(3)/6)*dt),
        !       J_2 = J(t_n + (1/2 + sqrt(3)/6)*dt),
        !       a_1 = 1/2*(1+sqrt(3)),
        !       a_2 = 1/2*(1-sqrt(3)),
        !       b_1 = -sqrt(3),
        !       b_2 = sqrt(3)
        call mxll_propagate_expgauss2(this%hm, this%namespace, this%gr, this%st, this%tr_mxll, &
          current_time, algo%dt)
 
        updated_quantities = [character(16) :: "E field", "B field", "vector potential"]
 
      case (iteration_done)
        call maxwell_exec_end_of_timestep_tasks(this)
        done = .false.
 
      case default
        done = .false.
      end select
 
    end select
 
    call profiling_out(trim(this%namespace%get())//":"//trim(operation%id))
    pop_sub(maxwell_do_algorithmic_operation)
  end function maxwell_do_algorithmic_operation
 
  ! ---------------------------------------------------------
  logical function maxwell_is_tolerance_reached(this, tol) result(converged)
    class(maxwell_t),       intent(in)    :: this
    real(real64),           intent(in)    :: tol
 
    push_sub(maxwell_is_tolerance_reached)
 
    converged = .false.
 
    pop_sub(maxwell_is_tolerance_reached)
  end function maxwell_is_tolerance_reached
 
  ! ---------------------------------------------------------
  subroutine maxwell_init_interaction_as_partner(partner, interaction)
    class(maxwell_t),           intent(in)    :: partner
    class(interaction_surrogate_t), intent(inout) :: interaction
 
    push_sub(maxwell_init_interaction_as_partner)
 
    select type (interaction)
    type is (lorentz_force_t)
      ! Nothing to be initialized for the Lorentz force.
    type is (mxll_e_field_to_matter_t)
      call interaction%init_from_partner(partner%gr, partner%space, partner%namespace)
    type is (mxll_b_field_to_matter_t)
      call interaction%init_from_partner(partner%gr, partner%space, partner%namespace)
    type is (mxll_vec_pot_to_matter_t)
      call interaction%init_from_partner(partner%gr, partner%space, partner%namespace)
    class default
      message(1) = "Unsupported interaction."
      call messages_fatal(1, namespace=partner%namespace)
    end select
 
    pop_sub(maxwell_init_interaction_as_partner)
  end subroutine maxwell_init_interaction_as_partner
 
  ! ---------------------------------------------------------
  subroutine maxwell_copy_quantities_to_interaction(partner, interaction)
    class(maxwell_t),           intent(inout) :: partner
    class(interaction_surrogate_t), intent(inout) :: interaction
 
    integer :: ip
    complex(real64) :: interpolated_value(3)
    real(real64), allocatable :: b_field(:,:), vec_pot(:,:)
 
    push_sub(maxwell_copy_quantities_to_interaction)
    call profiling_in(trim(partner%namespace%get())//":"//"COPY_QUANTITY_INTER")
 
    select type (interaction)
    type is (lorentz_force_t)
      call mxll_get_batch(partner%st%rs_stateb, partner%st%rs_state, partner%gr%np, partner%st%dim)
 
      do ip = 1, interaction%system_np
        call mesh_interpolation_evaluate(partner%mesh_interpolate, partner%st%rs_state(:,1), &
          interaction%system_pos(:, ip), interpolated_value(1))
        call mesh_interpolation_evaluate(partner%mesh_interpolate, partner%st%rs_state(:,2), &
          interaction%system_pos(:, ip), interpolated_value(2))
        call mesh_interpolation_evaluate(partner%mesh_interpolate, partner%st%rs_state(:,3), &
          interaction%system_pos(:, ip), interpolated_value(3))
        call get_electric_field_vector(interpolated_value, interaction%partner_e_field(:, ip))
        call get_magnetic_field_vector(interpolated_value, 1, interaction%partner_b_field(:, ip))
      end do
 
    type is (mxll_e_field_to_matter_t)
      call mxll_get_batch(partner%st%rs_stateb, partner%st%rs_state, partner%gr%np, partner%st%dim)
      select case (interaction%type)
 
      case (mxll_field_total)
        call get_electric_field_state(partner%st%rs_state, partner%gr, interaction%partner_field)
 
      case (mxll_field_trans)
        call get_transverse_rs_state(partner%helmholtz, partner%st, partner%namespace)
        call get_electric_field_state(partner%st%rs_state_trans, partner%gr, interaction%partner_field)
 
      case (mxll_field_long)
        call partner%helmholtz%get_long_field(partner%namespace, partner%st%rs_state_long, total_field=partner%st%rs_state)
        call get_electric_field_state(partner%st%rs_state_long, partner%gr, interaction%partner_field)
 
      case default
        message(1) = "Unknown type of field requested by interaction."
        call messages_fatal(1, namespace=partner%namespace)
      end select
      call interaction%do_mapping()
 
    type is (mxll_vec_pot_to_matter_t)
      call mxll_get_batch(partner%st%rs_stateb, partner%st%rs_state, partner%gr%np, partner%st%dim)
      safe_allocate(b_field(1:partner%gr%np_part, 1:partner%gr%box%dim))
      safe_allocate(vec_pot(1:partner%gr%np_part, 1:partner%gr%box%dim))
      ! magnetic field is always transverse
      call get_magnetic_field_state(partner%st%rs_state, partner%gr, partner%st%rs_sign, b_field, &
        partner%st%mu(1:partner%gr%np), partner%gr%np)
      ! vector potential stored in partner_field
      call partner%helmholtz%get_vector_potential(partner%namespace, vec_pot, trans_field=b_field)
      ! in the convention used by the electronic system, the -1/c factor is included in the vector potential
      ! but we do not divide by it, because the Maxwell code is in SI units, so the vector potential units
      ! are suitable for a (p + q.A) minimal coupling interaction
      interaction%partner_field(1:partner%gr%np,1:partner%gr%box%dim) = &
        vec_pot(1:partner%gr%np,1:partner%gr%box%dim)
      safe_deallocate_a(b_field)
      safe_deallocate_a(vec_pot)
      call interaction%do_mapping()
 
    type is (mxll_b_field_to_matter_t)
      call mxll_get_batch(partner%st%rs_stateb, partner%st%rs_state, partner%gr%np, partner%st%dim)
      call get_magnetic_field_state(partner%st%rs_state, partner%gr, partner%st%rs_sign, &
        interaction%partner_field, partner%st%mu(1:partner%gr%np), partner%gr%np)
      call interaction%do_mapping()
 
    class default
      message(1) = "Unsupported interaction."
      call messages_fatal(1, namespace=partner%namespace)
    end select
 
    call profiling_out(trim(partner%namespace%get())//":"//"COPY_QUANTITY_INTER")
    pop_sub(maxwell_copy_quantities_to_interaction)
  end subroutine maxwell_copy_quantities_to_interaction
 
  ! ---------------------------------------------------------
  subroutine maxwell_output_start(this)
    class(maxwell_t), intent(inout) :: this
 
    real(real64) :: itr_value
 
    push_sub(maxwell_output_start)
    call profiling_in(trim(this%namespace%get())//":"//"OUTPUT_START")
 
    select type (algo => this%algo)
    class is (propagator_t)
      call mxll_get_batch(this%st%rs_stateb, this%st%rs_state, this%gr%np, this%st%dim)
 
      call td_write_mxll_init(this%write_handler, this%namespace, this%iteration%counter(), algo%dt, &
        this%st%system_grp)
      if (this%st%fromScratch) then
        call td_write_mxll_iter(this%write_handler, this%space, this%gr, this%st, this%hm, this%helmholtz, algo%dt, &
          this%iteration%counter(), this%namespace)
        itr_value = this%iteration%value()
        call td_write_mxll_free_data(this%namespace, this%space, this%gr, this%st, this%hm, this%helmholtz, &
          this%outp, this%iteration%counter(), itr_value)
      end if
 
    end select
 
    call profiling_out(trim(this%namespace%get())//":"//"OUTPUT_START")
    pop_sub(maxwell_output_start)
  end subroutine maxwell_output_start
 
  ! ---------------------------------------------------------
  subroutine maxwell_output_write(this)
    class(maxwell_t), intent(inout) :: this
 
    logical :: stopping, reached_output_interval
    real(real64) :: itr_value
    integer :: iout
 
    push_sub(maxwell_output_write)
    call profiling_in(trim(this%namespace%get())//":"//"OUTPUT_WRITE")
 
    select type (algo => this%algo)
    class is (propagator_t)
      stopping = clean_stop(this%mc%master_comm)
 
      call td_write_mxll_iter(this%write_handler, this%space, this%gr, this%st, this%hm, this%helmholtz, algo%dt, &
        this%iteration%counter(), this%namespace)
 
      reached_output_interval = .false.
      do iout = 1, out_maxwell_max
        if (this%outp%output_interval(iout) > 0) then
          if (mod(this%iteration%counter(), this%outp%output_interval(iout)) == 0) then
            reached_output_interval = .true.
            exit
          end if
        end if
      end do
 
      if (reached_output_interval .or. stopping) then
        call mxll_get_batch(this%st%rs_stateb, this%st%rs_state, this%gr%np, this%st%dim)
 
        itr_value = this%iteration%value()
        call td_write_mxll_free_data(this%namespace, this%space, this%gr, this%st, this%hm, this%helmholtz, &
          this%outp, this%iteration%counter(), itr_value)
      end if
 
    end select
 
    call profiling_out(trim(this%namespace%get())//":"//"OUTPUT_WRITE")
    pop_sub(maxwell_output_write)
  end subroutine maxwell_output_write
 
  ! ---------------------------------------------------------
  subroutine maxwell_output_finish(this)
    class(maxwell_t), intent(inout) :: this
 
 
    push_sub(maxwell_output_finish)
 
    call profiling_in("MAXWELL_OUTPUT_FINISH")
 
    select type (algo => this%algo)
    class is (propagator_t)
      call td_write_mxll_end(this%write_handler)
    end select
 
    call profiling_out("MAXWELL_OUTPUT_FINISH")
 
    pop_sub(maxwell_output_finish)
  end subroutine maxwell_output_finish
 
  ! ---------------------------------------------------------
  subroutine maxwell_update_interactions_start(this)
    class(maxwell_t), intent(inout) :: this
 
    type(interaction_iterator_t) :: iter
    integer :: int_counter
 
    push_sub(maxwell_update_interactions_start)
 
    int_counter = 0
    call iter%start(this%interactions)
    do while (iter%has_next())
      select type (interaction => iter%get_next())
      class is (linear_medium_to_em_field_t)
        int_counter = int_counter + 1
      end select
    end do
 
    if (int_counter /= 0 .and. .not. allocated(this%hm%medium_boxes)) then
      safe_allocate(this%hm%medium_boxes(1:int_counter))
      this%hm%calc_medium_box = .true.
    end if
 
    pop_sub(maxwell_update_interactions_start)
  end subroutine maxwell_update_interactions_start
 
  ! ---------------------------------------------------------
  subroutine maxwell_update_interactions_finish(this)
    class(maxwell_t), intent(inout) :: this
 
    type(interaction_iterator_t) :: iter
    integer :: iint
 
    push_sub(maxwell_update_interactions_finish)
 
    iint = 0
 
    call iter%start(this%interactions)
    do while (iter%has_next())
      select type (interaction => iter%get_next())
      class is (linear_medium_to_em_field_t)
        if (allocated(this%hm%medium_boxes) .and. .not. this%hm%medium_boxes_initialized) then
          iint = iint + 1
          this%hm%medium_boxes(iint) = interaction%medium_box
        end if
      end select
    end do
 
    if (allocated(this%hm%medium_boxes) .and. .not. this%hm%medium_boxes_initialized) then
      call set_medium_rs_state(this%st, this%gr, this%hm)
      this%hm%medium_boxes_initialized = .true.
    end if
 
    if (this%hm%medium_boxes_initialized .and. this%hm%operator == faraday_ampere) then
      message(1) = "A linear medium has been defined in the input file but the Hamiltonian"
      message(2) = "type you specified is not capable of dealing with the medium."
      message(3) = "Please use MaxwellHamiltonianOperator = faraday_ampere_medium or simple to enable"
      message(4) = "the medium propagation."
      call messages_fatal(4, namespace=this%namespace)
    end if
 
    if (.not. this%hm%medium_boxes_initialized .and. &
      .not. any(this%hm%bc%bc_type == mxll_bc_medium) .and. &
      this%hm%operator == faraday_ampere_medium) then
      message(1) = "The variable MaxwellHamiltonianOperator has been defined as faraday_ampere_medium"
      message(2) = "in the input file but neither a medium boundary nor a linear medium system"
      message(3) = "has been defined. Please either use a different option for"
      message(4) = "MaxwellHamiltonianOperator or define a Maxwell medium."
      call messages_fatal(4, namespace=this%namespace)
    end if
 
    pop_sub(maxwell_update_interactions_finish)
  end subroutine maxwell_update_interactions_finish
 
  ! ---------------------------------------------------------
  subroutine maxwell_restart_write_data(this)
    class(maxwell_t), intent(inout) :: this
 
    integer :: ierr, err, zff_dim, id, id1, id2, ip_in, offset, iout
    logical :: pml_check, write_previous_state
    complex(real64), allocatable :: zff(:,:)
    type(interaction_iterator_t) :: iter
    class(interaction_t), pointer :: interaction
 
    push_sub(maxwell_restart_write_data)
    call profiling_in(trim(this%namespace%get())//":"//"RESTART_WRITE")
 
    ierr = 0
 
    if (mpi_world%is_root()) then
      do iout = 1, out_maxwell_max
        if (this%write_handler%out(iout)%write) then
          call write_iter_flush(this%write_handler%out(iout)%handle)
        end if
      end do
    end if
 
    if (.not. this%restart_dump%skip()) then
      pml_check = any(this%hm%bc%bc_ab_type(1:3) == option__maxwellabsorbingboundaries__cpml)
 
      message(1) = "Debug: Writing td_maxwell restart."
      call messages_info(1, namespace=this%namespace, debug_only=.true.)
 
      if (this%tr_mxll%bc_plane_waves) then
        zff_dim = 2 * this%st%dim
      else
        zff_dim = 1 * this%st%dim
      end if
      if (pml_check) then
        zff_dim = zff_dim + 18
      end if
      select type (prop => this%algo)
      type is (propagator_leapfrog_t)
        write_previous_state = .true.
        zff_dim = zff_dim + this%st%dim
      class  default
        write_previous_state = .false.
      end select
      if (pml_check .and. accel_is_enabled()) then
        call accel_read_buffer(this%hm%bc%pml%buff_conv_plus, &
          this%hm%bc%pml%points_number, 3, 3, this%hm%bc%pml%conv_plus)
        call accel_read_buffer(this%hm%bc%pml%buff_conv_minus, &
          this%hm%bc%pml%points_number, 3, 3, this%hm%bc%pml%conv_minus)
      end if
 
      call mxll_get_batch(this%st%rs_stateb, this%st%rs_state, this%gr%np, this%st%dim)
      if (write_previous_state) then
        call mxll_get_batch(this%st%rs_state_prevb, this%st%rs_state_prev, this%gr%np, this%st%dim)
      end if
 
      safe_allocate(zff(1:this%gr%np,1:zff_dim))
      zff = m_z0
 
      zff(1:this%gr%np, 1:this%st%dim) = this%st%rs_state(1:this%gr%np, 1:this%st%dim)
      if (this%tr_mxll%bc_plane_waves) then
        call mxll_get_batch(this%st%rs_state_plane_wavesb, this%st%rs_state_plane_waves, &
          this%gr%np, this%st%dim)
        zff(1:this%gr%np, this%st%dim+1:this%st%dim+this%st%dim) = &
          this%st%rs_state_plane_waves(1:this%gr%np, 1:this%st%dim)
        offset = 2*this%st%dim
      else
        offset = this%st%dim
      end if
      if (pml_check) then
        id = 0
        do id1 = 1, 3
          do id2 = 1, 3
            id = id + 1
            do ip_in = 1, this%hm%bc%pml%points_number
              zff(ip_in, offset+id) = this%hm%bc%pml%conv_plus(ip_in, id1, id2)
              zff(ip_in, offset+9+id) = this%hm%bc%pml%conv_minus(ip_in, id1, id2)
            end do
          end do
        end do
        offset = offset + 18
      end if
      if (write_previous_state) then
        zff(1:this%gr%np, offset+1:offset+this%st%dim) = &
          this%st%rs_state_prev(1:this%gr%np, 1:this%st%dim)
      end if
 
      call states_mxll_dump(this%restart_dump, this%st, this%space, this%gr, zff, zff_dim, err, this%iteration%counter())
      if (err /= 0) ierr = ierr + 1
 
      if (this%hm%current_density_from_medium) then
        !call this%current_interpolation%write_restart(this%gr, this%space, this%restart_dump, err)
        call iter%start(this%interactions)
        do while (iter%has_next())
          interaction => iter%get_next()
          select type (interaction)
          class is (current_to_mxll_field_t)
            call interaction%write_restart(this%gr, this%space, this%restart_dump, err)
          end select
        end do
        if (err /= 0) ierr = ierr + 1
      end if
 
      message(1) = "Debug: Writing td_maxwell restart done."
      call messages_info(1, namespace=this%namespace, debug_only=.true.)
 
      safe_deallocate_a(zff)
 
      if (ierr /=0) then
        message(1) = "Unable to write time-dependent Maxwell restart information."
        call messages_warning(1, namespace=this%namespace)
      end if
    end if
 
    call profiling_out(trim(this%namespace%get())//":"//"RESTART_WRITE")
    pop_sub(maxwell_restart_write_data)
  end subroutine maxwell_restart_write_data
 
  ! ---------------------------------------------------------
  ! this function returns true if restart data could be read
  logical function maxwell_restart_read_data(this)
    class(maxwell_t), intent(inout) :: this
 
    integer :: ierr, err, zff_dim, id, id1, id2, ip_in, offset
    logical :: pml_check, read_previous_state
    complex(real64), allocatable :: zff(:,:)
    type(interaction_iterator_t) :: iter
    class(interaction_t), pointer :: interaction
 
    push_sub(maxwell_restart_read_data)
    call profiling_in(trim(this%namespace%get())//":"//"RESTART_READ")
 
    if (.not. this%restart%skip()) then
      ierr = 0
      pml_check = any(this%hm%bc%bc_ab_type(1:3) == option__maxwellabsorbingboundaries__cpml)
 
      if (this%restart%skip()) then
        ierr = -1
        call profiling_out(trim(this%namespace%get())//":"//"RESTART_READ")
        pop_sub(maxwell_restart_read_data)
        return
      end if
 
      message(1) = "Debug: Reading td_maxwell restart."
      call messages_info(1, namespace=this%namespace, debug_only=.true.)
 
      if (this%tr_mxll%bc_plane_waves) then
        zff_dim = 2 * this%st%dim
      else
        zff_dim = 1 * this%st%dim
      end if
      if (pml_check) then
        zff_dim = zff_dim + 18
      end if
      select type (prop => this%algo)
      type is (propagator_leapfrog_t)
        read_previous_state = .true.
        zff_dim = zff_dim + this%st%dim
      class  default
        read_previous_state = .false.
      end select
 
      safe_allocate(zff(1:this%gr%np,1:zff_dim))
 
      call states_mxll_load(this%restart, this%st, this%gr, this%namespace, this%space, zff, &
        zff_dim, err, label = ": td_maxwell")
      this%st%rs_current_density_restart = .true.
 
      this%st%rs_state(1:this%gr%np,1:this%st%dim) = zff(1:this%gr%np, 1:this%st%dim)
      if (this%tr_mxll%bc_plane_waves) then
        this%st%rs_state_plane_waves(1:this%gr%np,1:this%st%dim) = &
          zff(1:this%gr%np,this%st%dim+1:this%st%dim+3)
        offset = 2*this%st%dim
      else
        offset = this%st%dim
      end if
      if (pml_check) then
        id = 0
        do id1 = 1, 3
          do id2 = 1, 3
            id = id+1
            do ip_in = 1, this%hm%bc%pml%points_number
              this%hm%bc%pml%conv_plus(ip_in,id1,id2)  = zff(ip_in, offset+  id)
              this%hm%bc%pml%conv_minus(ip_in,id1,id2) = zff(ip_in, offset+9+id)
            end do
          end do
        end do
        this%hm%bc%pml%conv_plus_old = this%hm%bc%pml%conv_plus
        this%hm%bc%pml%conv_minus_old = this%hm%bc%pml%conv_minus
        offset = offset + 18
      end if
      if (read_previous_state) then
        this%st%rs_state_prev(1:this%gr%np, 1:this%st%dim) = &
          zff(1:this%gr%np, offset+1:offset+this%st%dim)
      end if
 
      if (err /= 0) then
        ierr = ierr + 1
      end if
 
      message(1) = "Debug: Reading td restart done."
      call messages_info(1, namespace=this%namespace, debug_only=.true.)
 
      safe_deallocate_a(zff)
 
      if (pml_check .and. accel_is_enabled()) then
        call accel_write_buffer(this%hm%bc%pml%buff_conv_plus, this%hm%bc%pml%points_number, 3, 3, &
          this%hm%bc%pml%conv_plus)
        call accel_write_buffer(this%hm%bc%pml%buff_conv_minus, this%hm%bc%pml%points_number, 3, 3, &
          this%hm%bc%pml%conv_minus)
        call accel_write_buffer(this%hm%bc%pml%buff_conv_plus_old, this%hm%bc%pml%points_number, 3, 3, &
          this%hm%bc%pml%conv_plus_old)
      end if
 
      if (this%hm%current_density_from_medium) then
        !call this%current_interpolation%read_restart(this%gr, this%space, this%restart, err)
        call iter%start(this%interactions)
        do while (iter%has_next())
          interaction => iter%get_next()
          select type (interaction)
          class is (current_to_mxll_field_t)
            call interaction%read_restart(this%gr, this%space, this%restart, err)
          end select
        end do
        if (err /= 0) then
          ierr = ierr + 1
        end if
      end if
 
      ! set batches
      call mxll_set_batch(this%st%rs_stateb, this%st%rs_state, this%gr%np, this%st%dim)
      if (read_previous_state) then
        call mxll_set_batch(this%st%rs_state_prevb, this%st%rs_state_prev, this%gr%np, this%st%dim)
      end if
      call batch_set_zero(this%st%inhomogeneousb)
      if (this%tr_mxll%bc_plane_waves) then
        call mxll_set_batch(this%st%rs_state_plane_wavesb, this%st%rs_state_plane_waves, this%gr%np, this%st%dim)
      end if
 
      this%st%fromScratch = .false.
      maxwell_restart_read_data = .true.
    else
      message(1) = "Unable to read time-dependent Maxwell restart information: Starting from scratch"
      call messages_warning(1, namespace=this%namespace)
      maxwell_restart_read_data = .false.
    end if
 
    call profiling_out(trim(this%namespace%get())//":"//"RESTART_READ")
    pop_sub(maxwell_restart_read_data)
  end function maxwell_restart_read_data
 
  subroutine maxwell_update_kinetic_energy(this)
    class(maxwell_t), intent(inout) :: this
 
    push_sub(maxwell_update_kinetic_energy)
 
    ! the energy has already been computed at the end of the timestep
 
    ! the energy of the EM wave is computed and stored in energy_mxll%energy;
    ! energy_mxll%energy = energy_mxll%e_energy + energy_mxll%b_energy
    ! here this%hm%energy is 'energy_mxll'
    this%kinetic_energy = this%hm%energy%energy
 
    pop_sub(maxwell_update_kinetic_energy)
  end subroutine maxwell_update_kinetic_energy
 
  ! ---------------------------------------------------------
  subroutine maxwell_exec_end_of_timestep_tasks(this)
    class(maxwell_t), intent(inout) :: this
 
    push_sub(maxwell_exec_end_of_timestep_tasks)
    call profiling_in(trim(this%namespace%get())//":"//"END_TIMESTEP")
 
    ! calculate Maxwell energy
    call energy_mxll_calc_batch(this%gr, this%st, this%hm, this%hm%energy, this%st%rs_stateb, this%st%rs_state_plane_wavesb)
 
    ! get RS state values for selected points
    call get_rs_state_batch_selected_points(this%st%selected_points_rs_state(:,:), this%st%rs_stateb, &
      this%st, this%gr)
 
    call profiling_out(trim(this%namespace%get())//":"//"END_TIMESTEP")
    pop_sub(maxwell_exec_end_of_timestep_tasks)
  end subroutine maxwell_exec_end_of_timestep_tasks
 
  ! ---------------------------------------------------------
  subroutine maxwell_get_current(this, time, current)
    class(maxwell_t),             intent(inout) :: this
    real(real64),                 intent(in)    :: time
    complex(real64), contiguous,  intent(inout) :: current(:, :)
 
    type(interaction_iterator_t) :: iter
    complex(real64), allocatable :: current_density_ext(:, :)
 
    push_sub(maxwell_get_current)
 
    safe_allocate(current_density_ext(1:this%gr%np, 1:this%st%dim))
    current = m_z0
    if (this%hm%current_density_from_medium) then
      ! interpolate current from interaction
      call iter%start(this%interactions)
      do while (iter%has_next())
        select type (interaction => iter%get_next())
        class is (current_to_mxll_field_t)
          call interaction%interpolate(time, current_density_ext)
          call lalg_axpy(this%gr%np, 3, m_one, current_density_ext * this%hm%current_factor, current)
        end select
      end do
    end if
    ! calculation of external RS density
    if (this%hm%current_density_ext_flag) then
      call get_rs_density_ext(this%st, this%space, this%gr, time, current_density_ext)
      call lalg_axpy(this%gr%np, 3, m_one, current_density_ext, current)
    end if
    safe_deallocate_a(current_density_ext)
 
    pop_sub(maxwell_get_current)
  end subroutine maxwell_get_current
 
  ! ---------------------------------------------------------
  subroutine maxwell_finalize(this)
    type(maxwell_t), intent(inout) :: this
 
 
    push_sub(maxwell_finalize)
 
    call profiling_in("MAXWELL_FINALIZE")
 
    call system_end(this)
 
    ! free memory
    safe_deallocate_a(this%rs_state_init)
    call hamiltonian_mxll_end(this%hm)
 
    call multicomm_end(this%mc)
 
    call this%st%rs_stateb%end()
    call this%st%rs_state_prevb%end()
    call this%st%inhomogeneousb%end()
    if (this%tr_mxll%bc_plane_waves) then
      call this%st%rs_state_plane_wavesb%end()
    end if
 
    call states_mxll_end(this%st)
 
    call grid_end(this%gr)
 
    call this%restart%end()
    call this%restart_dump%end()
 
    call poisson_end(this%st%poisson)
 
    call profiling_out("MAXWELL_FINALIZE")
 
    pop_sub(maxwell_finalize)
  end subroutine maxwell_finalize
 
end module maxwell_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: