charged_particle.F90

Go to the documentation of this file.

!! Copyright (C) 2020 H. Appel, S. Ohlmann, M. Oliveira, N. Tancogne-Dejean
!!
!! 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 charged_particle_oct_m
  use algorithm_oct_m
  use classical_particle_oct_m
  use coulomb_force_oct_m
  use current_to_mxll_field_oct_m
  use debug_oct_m
  use global_oct_m
  use interaction_oct_m
  use interaction_surrogate_oct_m
  use interaction_enum_oct_m
  use lorentz_force_oct_m
  use messages_oct_m
  use multisystem_debug_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use quantity_oct_m
  use space_oct_m
  use system_oct_m
 
  implicit none
 
  private
  public ::               &
    charged_particle_t,   &
    charged_particle_init
 
  type, extends(classical_particle_t) :: charged_particle_t
    ! Note that we extend classical_particle_t and not charged_particles_t. Both
    ! choices have pros and cons. In the end extending classical_particle_t was
    ! chosen because charged_particles_t was only introduced after
    ! charged_particle_t.
    private
 
    real(real64) :: charge(1)
 
  contains
    procedure :: init_interaction => charged_particle_init_interaction
    procedure :: initialize => charged_particle_initialize
    procedure :: is_tolerance_reached => charged_particle_is_tolerance_reached
    procedure :: update_quantity => charged_particle_update_quantity
    procedure :: init_interaction_as_partner => charged_particle_init_interaction_as_partner
    procedure :: copy_quantities_to_interaction => charged_particle_copy_quantities_to_interaction
  end type charged_particle_t
 
  interface charged_particle_t
    procedure charged_particle_constructor
  end interface charged_particle_t
 
contains
 
  ! ---------------------------------------------------------
  function charged_particle_constructor(namespace, grp) result(sys)
    class(charged_particle_t), pointer  :: sys
    type(namespace_t),       intent(in) :: namespace
    type(mpi_grp_t),         intent(in) :: grp
 
    push_sub(charged_particle_constructor)
 
    allocate(sys)
 
    call charged_particle_init(sys, namespace, grp)
 
    pop_sub(charged_particle_constructor)
  end function charged_particle_constructor
 
  ! ---------------------------------------------------------
  subroutine charged_particle_init(this, namespace, grp)
    class(charged_particle_t), intent(inout) :: this
    type(namespace_t),         intent(in)    :: namespace
    type(mpi_grp_t),           intent(in)    :: grp
 
    push_sub(charged_particle_init)
 
    call classical_particle_init(this%classical_particle_t, namespace, grp)
 
    !%Variable ParticleCharge
    !%Type float
    !%Section ClassicalParticles
    !%Description
    !% Charge of classical particle
    !%End
    call parse_variable(namespace, 'ParticleCharge', m_one, this%charge(1))
    call messages_print_var_value('ParticleCharge', this%charge(1), namespace=namespace)
 
    this%supported_interactions = [this%supported_interactions, lorentz_force, coulomb_force]
    this%supported_interactions_as_partner = [this%supported_interactions_as_partner, coulomb_force, current_to_mxll_field]
 
    call this%quantities%add(quantity_t("charge", updated_on_demand = .false., always_available = .true.))
    call this%quantities%add(quantity_t("current", parents=[character(8) :: "charge", "velocity"]))
 
 
    pop_sub(charged_particle_init)
  end subroutine charged_particle_init
 
  ! ---------------------------------------------------------
  subroutine charged_particle_init_interaction(this, interaction)
    class(charged_particle_t), target, intent(inout) :: this
    class(interaction_t),              intent(inout) :: interaction
 
    push_sub(charged_particle_init_interaction)
 
    select type (interaction)
    type is (coulomb_force_t)
      call interaction%init(this%space%dim, 1, this%charge, this%pos)
    type is (lorentz_force_t)
      call interaction%init(this%space%dim, 1, this%charge, this%pos, this%vel, this%namespace)
    class default
      call this%classical_particle_t%init_interaction(interaction)
    end select
 
    pop_sub(charged_particle_init_interaction)
  end subroutine charged_particle_init_interaction
 
  ! ---------------------------------------------------------
  subroutine charged_particle_initialize(this)
    class(charged_particle_t), intent(inout) :: this
 
    push_sub(charged_particle_initialize)
 
    call this%classical_particle_t%initialize()
 
    pop_sub(charged_particle_initialize)
  end subroutine charged_particle_initialize
 
  ! ---------------------------------------------------------
  logical function charged_particle_is_tolerance_reached(this, tol) result(converged)
    class(charged_particle_t), intent(in) :: this
    real(real64),              intent(in) :: tol
 
    push_sub(charged_particle_is_tolerance_reached)
 
    converged = this%classical_particle_t%is_tolerance_reached(tol)
 
    pop_sub(charged_particle_is_tolerance_reached)
  end function charged_particle_is_tolerance_reached
 
  ! ---------------------------------------------------------
  subroutine charged_particle_update_quantity(this, label)
    class(charged_particle_t), intent(inout) :: this
    character(len=*),          intent(in)    :: label
 
    push_sub(charged_particle_update_quantity)
 
    select case (label)
    case ("current")
      ! we do not explicitly compute the current as it is only needed
      ! by the current_to_mxll_field_t interaction where it is mapped
      ! to a grid using charge, position and velocity
    case default
      ! Other quantities should be handled by the parent class
      call this%classical_particle_t%update_quantity(label)
    end select
 
    pop_sub(charged_particle_update_quantity)
  end subroutine charged_particle_update_quantity
 
  ! ---------------------------------------------------------
  subroutine charged_particle_init_interaction_as_partner(partner, interaction)
    class(charged_particle_t), intent(in)    :: partner
    class(interaction_surrogate_t), intent(inout) :: interaction
 
    push_sub(charged_particle_init_interaction_as_partner)
 
    select type (interaction)
    type is (coulomb_force_t)
      interaction%partner_np = 1
      safe_allocate(interaction%partner_charge(1))
      safe_allocate(interaction%partner_pos(1:partner%space%dim, 1))
    type is (current_to_mxll_field_t)
      interaction%partner_np = 1
      interaction%grid_based_partner = .false.
      safe_allocate(interaction%partner_charge(1))
      safe_allocate(interaction%partner_pos(1:partner%space%dim, 1))
      safe_allocate(interaction%partner_vel(1:partner%space%dim, 1))
    class default
      ! Other interactions should be handled by the parent class
      call partner%classical_particle_t%init_interaction_as_partner(interaction)
    end select
 
    pop_sub(charged_particle_init_interaction_as_partner)
  end subroutine charged_particle_init_interaction_as_partner
 
  ! ---------------------------------------------------------
  subroutine charged_particle_copy_quantities_to_interaction(partner, interaction)
    class(charged_particle_t), intent(inout) :: partner
    class(interaction_surrogate_t), intent(inout) :: interaction
 
    push_sub(charged_particle_copy_quantities_to_interaction)
 
    select type (interaction)
    type is (coulomb_force_t)
      interaction%partner_charge(1) = partner%charge(1)
      interaction%partner_pos(:,1) = partner%pos(:, 1)
 
    type is (current_to_mxll_field_t)
      interaction%partner_charge(1) = partner%charge(1)
      interaction%partner_pos(:,1) = partner%pos(:, 1)
      interaction%partner_vel(:,1) = partner%vel(:, 1)
      call interaction%do_mapping()
 
    class default
      call partner%classical_particle_t%copy_quantities_to_interaction(interaction)
    end select
 
    pop_sub(charged_particle_copy_quantities_to_interaction)
  end subroutine charged_particle_copy_quantities_to_interaction
 
end module charged_particle_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: