lorentz_force.F90

Go to the documentation of this file.

!! Copyright (C) 2020 Heiko Appel
!!
!! 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 lorentz_force_oct_m
  use debug_oct_m
  use force_interaction_oct_m
  use global_oct_m
  use interaction_oct_m
  use interaction_partner_oct_m
  use messages_oct_m
  use namespace_oct_m
  use profiling_oct_m
  use quantity_oct_m
 
  implicit none
 
  private
  public ::                &
    lorentz_force_t
 
  type, extends(force_interaction_t) :: lorentz_force_t
    private
    real(real64), pointer :: system_charge(:) => null() 
    real(real64), pointer, public :: system_pos(:,:) => null() 
    real(real64), pointer :: system_vel(:,:) => null() 
 
    real(real64), allocatable, public :: partner_e_field(:,:)
    real(real64), allocatable, public :: partner_b_field(:,:)
 
  contains
    procedure :: init => lorentz_force_init
    procedure :: calculate => lorentz_force_calculate
    procedure :: calculate_energy => lorentz_force_calculate_energy
    final :: lorentz_force_finalize
  end type lorentz_force_t
 
  interface lorentz_force_t
    module procedure lorentz_force_constructor
  end interface lorentz_force_t
 
contains
 
  ! ---------------------------------------------------------
  function lorentz_force_constructor(partner) result(this)
    class(interaction_partner_t), target, intent(inout) :: partner
    class(lorentz_force_t),               pointer       :: this
 
    push_sub(lorentz_force_constructor)
 
    allocate(this)
 
    this%label = "lorenz_force"
 
    this%partner => partner
 
    ! The Lorentz force needs the position, velocity and charge of the system
    this%system_quantities = [position, velocity, charge]
 
    ! The Lorenz force needs the E and B field of the interaction partner at the particle position
    this%couplings_from_partner = [e_field, b_field]
    this%intra_interaction = .false. ! This interaction does not support intra-interactions
 
 
    pop_sub(lorentz_force_constructor)
  end function lorentz_force_constructor
 
  ! ---------------------------------------------------------
  subroutine lorentz_force_init(this, dim, system_np, system_charge, system_pos, system_vel, namespace)
    class(lorentz_force_t),               intent(inout) :: this
    integer,                              intent(in)    :: dim
    integer,                              intent(in)    :: system_np
    real(real64),                 target, intent(in)    :: system_charge(:)
    real(real64),                 target, intent(in)    :: system_pos(:,:)
    real(real64),                 target, intent(in)    :: system_vel(:,:)
    type(namespace_t),                    intent(in)    :: namespace
 
    push_sub(lorentz_force_init)
 
    message(1) = "Energies for Lorentz forces are not yet implemented, and currently set to 0."
    call messages_warning(1, namespace=namespace)
 
    this%dim = dim
    this%system_np = system_np
    this%energy = m_zero
    safe_allocate(this%force(dim, system_np))
    this%force = m_zero
 
    safe_allocate(this%partner_e_field(1:dim, 1:system_np))
    this%partner_e_field = m_zero
    safe_allocate(this%partner_b_field(1:dim, 1:system_np))
    this%partner_b_field = m_zero
 
    this%system_charge => system_charge
    this%system_pos => system_pos
    this%system_vel => system_vel
 
    pop_sub(lorentz_force_init)
  end subroutine lorentz_force_init
 
  ! ---------------------------------------------------------
  subroutine lorentz_force_calculate(this)
    class(lorentz_force_t),             intent(inout) :: this
 
    integer :: ip
 
    push_sub(lorentz_force_calculate)
 
    ! Curl is defined only in 3D
    assert(this%dim == 3)
    assert(.not. this%intra_interaction)
 
    do ip = 1, this%system_np
      this%force(1, ip) = this%partner_e_field(1, ip) + &
        this%system_vel(2, ip)*this%partner_b_field(3, ip) - this%system_vel(3, ip)*this%partner_b_field(2, ip)
 
      this%force(2, ip) = this%partner_e_field(2, ip) + &
        this%system_vel(3, ip)*this%partner_b_field(1, ip) - this%system_vel(1, ip)*this%partner_b_field(3, ip)
 
      this%force(3, ip) = this%partner_e_field(3, ip) + &
        this%system_vel(1, ip)*this%partner_b_field(2, ip) - this%system_vel(2, ip)*this%partner_b_field(1, ip)
 
      this%force(1:3, ip) = this%force(1:3, ip)*this%system_charge(ip)
    end do
 
    pop_sub(lorentz_force_calculate)
  end subroutine lorentz_force_calculate
 
  ! ---------------------------------------------------------
  subroutine lorentz_force_calculate_energy(this)
    class(lorentz_force_t),             intent(inout) :: this
 
    integer :: ip
    real(real64)   :: power
 
    push_sub(lorentz_force_calculate_energy)
 
    ! the rate at which the energy is transferred from the EM field to the particle
    ! the B field does not contribute any energy to the particles
    power = m_zero
    do ip = 1, this%system_np
      power = power - dot_product(this%system_vel(1:3,ip), this%force(1:3,ip))
    end do
 
    this%energy = m_zero
    !TODO: We need to implement the proper time integration of the power.
    !      However, to that end, we first need to implement another clock at which the energies are
    !      calculated.
    ! this%energy = this%energy + power*this%partner%iteration%get_time_step()
 
    pop_sub(lorentz_force_calculate_energy)
  end subroutine lorentz_force_calculate_energy
 
 
  ! ---------------------------------------------------------
  subroutine lorentz_force_finalize(this)
    type(lorentz_force_t), intent(inout) :: this
 
    push_sub(lorentz_force_finalize)
 
    this%force = m_zero
    nullify(this%system_charge)
    nullify(this%system_pos)
    nullify(this%system_vel)
    safe_deallocate_a(this%force)
    safe_deallocate_a(this%partner_e_field)
    safe_deallocate_a(this%partner_b_field)
 
    call interaction_end(this)
 
    pop_sub(lorentz_force_finalize)
  end subroutine lorentz_force_finalize
 
end module lorentz_force_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: