Merge branch '333-Add-Lorentz-force-for-charged-particles' into 'develop'

Add Lorentz force for charged classical particles

See merge request !896

src/Makefile.am

@@ -493,6 +493,7 @@ system_f_srcs =                                 \
 	system/forces.F90                       \
 	system/ghost_interaction.F90            \
 	system/interaction_abst.F90             \
+	system/interaction_lorentz_force.F90    \
 	system/interaction_partner.F90          \
 	system/interaction_gravity.F90          \
 	system/interaction_with_partner.F90     \

src/system/charged_particle.F90

@@ -23,6 +23,7 @@ module charged_particle_oct_m
   use clock_oct_m
   use global_oct_m
   use interaction_abst_oct_m
+  use interaction_lorentz_force_oct_m
   use io_oct_m
   use iso_c_binding
   use messages_oct_m
@@ -61,6 +62,7 @@ module charged_particle_oct_m
     procedure :: store_current_status => charged_particle_store_current_status
     procedure :: update_quantity => charged_particle_update_quantity
     procedure :: update_exposed_quantity => charged_particle_update_exposed_quantity
+    procedure :: copy_quantities_to_interaction => charged_particle_copy_quantities_to_interaction
     procedure :: update_interactions_start => charged_particle_update_interactions_start
     procedure :: update_interactions_finish => charged_particle_update_interactions_finish
     final :: charged_particle_finalize
@@ -112,6 +114,9 @@ contains
     call parse_variable(namespace, 'ClassicalParticleCharge', M_ONE, this%charge)
     call messages_print_var_value(stdout, 'ClassicalParticleCharge', this%charge)
 
+    this%quantities(CHARGE)%required = .true.
+    this%quantities(CHARGE)%protected = .true.
+
     POP_SUB(charged_particle_init)
   end subroutine charged_particle_init
 
@@ -120,8 +125,21 @@ contains
     class(charged_particle_t), target, intent(inout) :: this
     class(system_abst_t),              intent(inout) :: partner
 
+    class(interaction_lorentz_force_t), pointer :: lorentz_force
+    type(interaction_lorentz_force_t) :: lorentz_force_t
+
     PUSH_SUB(charged_particle_add_interaction_partner)
 
+    if (partner%has_interaction(lorentz_force_t)) then
+      lorentz_force => interaction_lorentz_force_t(this%space%dim, partner)
+      this%quantities(POSITION)%required = .true.
+      this%quantities(VELOCITY)%required = .true.
+      this%quantities(CHARGE)%required = .true.
+      lorentz_force%system_pos => this%pos
+      lorentz_force%system_vel => this%vel
+      lorentz_force%system_charge => this%charge
+      call this%interactions%add(lorentz_force)
+    end if
     call this%classical_particle_t%add_interaction_partner(partner)
 
     POP_SUB(charged_particle_add_interaction_partner)
@@ -134,7 +152,12 @@ contains
 
     PUSH_SUB(charged_particle_has_interaction)
 
-    charged_particle_has_interaction = this%classical_particle_t%has_interaction(interaction)
+    select type (interaction)
+    type is (interaction_lorentz_force_t)
+      charged_particle_has_interaction = .true.
+    class default
+      charged_particle_has_interaction = this%classical_particle_t%has_interaction(interaction)
+    end select
 
     POP_SUB(charged_particle_has_interaction)
   end function charged_particle_has_interaction
@@ -243,10 +266,10 @@ contains
     case (CHARGE)
       ! The charged particle has a charge, but it is not necessary to update it, as it does not change with time.
       call this%quantities(iq)%clock%set_time(requested_time)
+    case default
+      ! Other quantities should be handled by the parent class
+      call this%classical_particle_t%update_quantity(iq, requested_time)
     end select
-    ! Note: the assert for protected quantities and the default case are taken care of by
-    ! the update of quantities of the classical particle, which follows next:
-    call this%classical_particle_t%update_quantity(iq, requested_time)
 
     POP_SUB(charged_particle_update_quantity)
   end subroutine charged_particle_update_quantity
@@ -263,14 +286,31 @@ contains
     case (CHARGE)
       ! The charged particle has a charge, but it is not necessary to update it, as it does not change with time.
       call this%quantities(iq)%clock%set_time(requested_time)
+    case default
+      ! Other quantities should be handled by the parent class
+      call this%classical_particle_t%update_exposed_quantity(iq, requested_time)
     end select
-    ! Note: the assert for protected quantities and the default case are taken care of by
-    ! the update of exposed quantities of the classical particle, which follows next:
-    call this%classical_particle_t%update_exposed_quantity(iq, requested_time)
 
     POP_SUB(charged_particle_update_exposed_quantity)
   end subroutine charged_particle_update_exposed_quantity
 
+  ! ---------------------------------------------------------
+  subroutine charged_particle_copy_quantities_to_interaction(this, interaction)
+    class(charged_particle_t), intent(inout) :: this
+    class(interaction_abst_t), intent(inout) :: interaction
+
+    PUSH_SUB(classical_particle_copy_quantities_to_interaction)
+
+    select type (interaction)
+    type is (interaction_lorentz_force_t)
+      ! Nothing to copy
+    class default
+      call this%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
+
   ! ---------------------------------------------------------
   subroutine charged_particle_update_interactions_start(this)
     class(charged_particle_t), intent(inout) :: this

src/system/classical_particle.F90

@@ -24,6 +24,7 @@ module classical_particle_oct_m
   use global_oct_m
   use interaction_abst_oct_m
   use interaction_gravity_oct_m
+  use interaction_lorentz_force_oct_m
   use ghost_interaction_oct_m
   use io_oct_m
   use iso_c_binding
@@ -41,13 +42,11 @@ module classical_particle_oct_m
   implicit none
 
   private
-  public ::                 &
-    classical_particle_t,   &
+  public ::                  &
+    classical_particle_t,    &
     classical_particle_init
 
-  type, extends(system_abst_t) :: classical_particle_t
-    private
-
+   type, extends(system_abst_t) :: classical_particle_t
     FLOAT :: mass
     FLOAT :: pos(1:MAX_DIM)
     FLOAT :: vel(1:MAX_DIM)
@@ -238,7 +237,7 @@ contains
   ! ---------------------------------------------------------
   subroutine classical_particle_do_td(this, operation)
     class(classical_particle_t), intent(inout) :: this
-    integer,                 intent(in)    :: operation
+    integer,                     intent(in)    :: operation
 
     integer :: ii
 
@@ -456,7 +455,7 @@ contains
     call write_iter_double(this%output_handle, this%vel, this%space%dim)
     call write_iter_double(this%output_handle, this%tot_force, this%space%dim)
     call write_iter_nl(this%output_handle)
-    
+
     POP_SUB(classical_particle_output_write)
   end subroutine classical_particle_output_write
 
@@ -517,6 +516,8 @@ contains
     type is (interaction_gravity_t)
       interaction%partner_mass = this%mass
       interaction%partner_pos = this%pos
+    type is (interaction_lorentz_force_t)
+      ! Nothing to copy
     type is (ghost_interaction_t)
       ! Nothing to copy
     class default
@@ -552,6 +553,8 @@ contains
     call iter%start(this%interactions)
     do while (iter%has_next())
       select type (interaction => iter%get_next_interaction())
+      type is (interaction_lorentz_force_t)
+        this%tot_force(1:this%space%dim) = this%tot_force(1:this%space%dim) + interaction%force(1:this%space%dim)
       type is (interaction_gravity_t)
         this%tot_force(1:this%space%dim) = this%tot_force(1:this%space%dim) + interaction%force(1:this%space%dim)
       type is (ghost_interaction_t)

src/system/interaction_lorentz_force.F90

+!! 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 interaction_lorentz_force_oct_m
+  use global_oct_m
+  use interaction_with_partner_oct_m
+  use interaction_partner_oct_m
+  use messages_oct_m
+  use profiling_oct_m
+  use quantity_oct_m
+
+  implicit none
+
+  private
+  public ::                &
+    interaction_lorentz_force_t
+
+  type, extends(interaction_with_partner_t) :: interaction_lorentz_force_t
+    integer :: dim
+
+    FLOAT :: force(MAX_DIM)
+
+    FLOAT, pointer :: system_charge
+    FLOAT, pointer :: system_pos(:)
+    FLOAT, pointer :: system_vel(:)
+
+    FLOAT :: partner_mass
+    FLOAT, allocatable :: partner_E_field(:)
+    FLOAT, allocatable :: partner_B_field(:)
+
+  contains
+    procedure :: calculate => interaction_lorentz_force_calculate
+    procedure :: end => interaction_lorentz_force_end
+    final :: interaction_lorentz_force_finalize
+  end type interaction_lorentz_force_t
+
+  interface interaction_lorentz_force_t
+    module procedure interaction_lorentz_force_init
+  end interface interaction_lorentz_force_t
+
+contains
+
+  ! ---------------------------------------------------------
+
+  function interaction_lorentz_force_init(dim, partner) result(this)
+    integer,                              intent(in)    :: dim
+    class(interaction_partner_t), target, intent(inout) :: partner
+    class(interaction_lorentz_force_t),   pointer       :: this
+
+    PUSH_SUB(interaction_lorentz_force_init)
+
+    SAFE_ALLOCATE(this)
+
+    this%dim = dim
+    this%partner => partner
+
+    ! For the Lorentz force we need the position, velocity and charge of the
+    ! system and the E and B field of the interaction partner at the particle position
+    this%n_system_quantities = 3
+    this%n_partner_quantities = 2
+    SAFE_ALLOCATE(this%system_quantities(this%n_system_quantities))
+    SAFE_ALLOCATE(this%partner_quantities(this%n_partner_quantities))
+    this%system_quantities(1) = POSITION
+    this%system_quantities(2) = VELOCITY
+    this%system_quantities(3) = CHARGE
+    this%partner_quantities(1) = E_FIELD
+    this%partner_quantities(2) = B_FIELD
+    SAFE_ALLOCATE(this%partner_E_field(dim))
+    SAFE_ALLOCATE(this%partner_B_field(dim))
+
+    POP_SUB(interaction_lorentz_force_init)
+  end function interaction_lorentz_force_init
+
+  ! ---------------------------------------------------------
+  subroutine interaction_lorentz_force_calculate(this)
+    class(interaction_lorentz_force_t), intent(inout) :: this
+
+    PUSH_SUB(interaction_lorentz_force_calculate)
+
+    ! Curl is defined only in 3D
+    ASSERT(this%dim == 3)
+
+    this%force(1) = this%system_charge * (this%partner_E_field(1) + &
+      this%system_vel(2) * this%partner_B_field(3) - this%system_vel(3) * this%partner_B_field(2))
+    this%force(2) = this%system_charge * (this%partner_E_field(2) + &
+      this%system_vel(3) * this%partner_B_field(1) - this%system_vel(1) * this%partner_B_field(3))
+    this%force(3) = this%system_charge * (this%partner_E_field(3) + &
+      this%system_vel(1) * this%partner_B_field(2) - this%system_vel(2) * this%partner_B_field(1))
+
+    POP_SUB(interaction_lorentz_force_calculate)
+  end subroutine interaction_lorentz_force_calculate
+
+  ! ---------------------------------------------------------
+  subroutine interaction_lorentz_force_end(this)
+    class(interaction_lorentz_force_t), intent(inout) :: this
+
+    PUSH_SUB(interaction_lorentz_force_end)
+
+    nullify(this%partner)
+    this%force = M_ZERO
+    nullify(this%system_charge)
+    nullify(this%system_pos)
+    nullify(this%system_vel)
+    SAFE_DEALLOCATE_A(this%partner_E_field)
+    SAFE_DEALLOCATE_A(this%partner_B_field)
+
+    call interaction_with_partner_end(this)
+
+    POP_SUB(interaction_lorentz_force_end)
+  end subroutine interaction_lorentz_force_end
+
+  ! ---------------------------------------------------------
+  subroutine interaction_lorentz_force_finalize(this)
+    type(interaction_lorentz_force_t), intent(inout) :: this
+
+    PUSH_SUB(interaction_lorentz_force_finalize)
+
+    call this%end()
+
+    POP_SUB(interaction_lorentz_force_finalize)
+  end subroutine interaction_lorentz_force_finalize
+
+end module interaction_lorentz_force_oct_m
+
+!! Local Variables:
+!! mode: f90
+!! coding: utf-8
+!! End:

src/system/system_mxll.F90

@@ -25,6 +25,7 @@ module system_mxll_oct_m
   use distributed_oct_m
   use geometry_oct_m
   use ghost_interaction_oct_m
+  use interaction_lorentz_force_oct_m
   use global_oct_m
   use grid_oct_m
   use hamiltonian_mxll_oct_m
@@ -34,6 +35,7 @@ module system_mxll_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 namespace_oct_m
@@ -77,6 +79,8 @@ module system_mxll_oct_m
     type(output_t)               :: outp  !< the output
     type(multicomm_t)            :: mc    !< index and domain communicators
 
+    type(mesh_interpolation_t)   :: mesh_interpolate
+
     type(propagator_mxll_t)      :: tr_mxll   !< contains the details of the Maxwell time-evolution
     type(td_write_t)             :: write_handler
     type(c_ptr)                  :: output_handle
@@ -177,6 +181,8 @@ contains
     this%quantities(E_FIELD)%required = .true.
     this%quantities(B_FIELD)%required = .true.
 
+    call mesh_interpolation_init(this%mesh_interpolate, this%gr%mesh)
+
     POP_SUB(system_mxll_init)
 
   contains
@@ -223,7 +229,12 @@ contains
 
     PUSH_SUB(system_mxll_has_interaction)
 
-    system_mxll_has_interaction = .false.
+    select type (interaction)
+    type is (interaction_lorentz_force_t)
+      system_mxll_has_interaction = .true.
+    class default
+      system_mxll_has_interaction = .false.
+    end select
 
     POP_SUB(system_mxll_has_interaction)
   end function system_mxll_has_interaction
@@ -539,7 +550,7 @@ contains
     select case (iq)
     case default
       message(1) = "Incompatible quantity."
-      call messages_fatal(1)
+!      call messages_fatal(1)
     end select
 
     POP_SUB(system_mxll_update_exposed_quantity)
@@ -572,14 +583,26 @@ contains
 
     ! ---------------------------------------------------------
   subroutine system_mxll_copy_quantities_to_interaction(this, interaction)
-    class(system_mxll_t),          intent(inout) :: this
-    class(interaction_abst_t),        intent(inout) :: interaction
+    class(system_mxll_t),      intent(inout) :: this
+    class(interaction_abst_t), intent(inout) :: interaction
+
+    CMPLX :: interpolated_value(3)
+    FLOAT :: e_field(3)
+    FLOAT :: b_field(3)
 
     PUSH_SUB(system_mxll_copy_quantities_to_interaction)
 
     select type (interaction)
     type is (ghost_interaction_t)
       ! Nothing to copy
+    type is (interaction_lorentz_force_t)
+      call mesh_interpolation_evaluate(this%mesh_interpolate, this%st%rs_state(:,1), interaction%system_pos, interpolated_value(1))
+      call mesh_interpolation_evaluate(this%mesh_interpolate, this%st%rs_state(:,2), interaction%system_pos, interpolated_value(2))
+      call mesh_interpolation_evaluate(this%mesh_interpolate, this%st%rs_state(:,3), interaction%system_pos, interpolated_value(3))
+      call get_electric_field_vector(interpolated_value, e_field)
+      call get_magnetic_field_vector(interpolated_value, 1, b_field)
+      interaction%partner_E_field = e_field
+      interaction%partner_B_field = b_field
     class default
       message(1) = "Unsupported interaction."
       call messages_fatal(1)

testsuite/multisystem/05-lorentz-force.01-charged_particle_coupled_to_maxwell.inp

+# ----- Calculation mode and parallelization ----------------------------------
+
+CalculationMode   = td
+RestartWrite = no
+ExperimentalFeatures = yes
+
+Dimensions = 3
+
+%Systems
+  'Maxwell' | maxwell
+  'ChargedParticle' | charged_particle
+%
+
+# ----- Charged particle ------------------------------------------------------
+
+ChargedParticle.ClassicalParticleMass = 1
+ChargedParticle.ClassicalParticleCharge = 1
+
+%ChargedParticle.ClassicalParticleInitialPosition
+ 1 | 0 | 0
+%
+%ChargedParticle.ClassicalParticleInitialVelocity
+ 0 | 0 | 0.01
+%
+
+ChargedParticle.TDSystemPropagator = verlet
+
+
+# ----- Maxwell variables -----------------------------------------------------
+
+Maxwell.ParDomains = auto
+Maxwell.ParStates  = no
+
+HamiltonianOperator = faraday_ampere
+
+Maxwell.TDSystemPropagator = exp_mid
+
+%MaxwellBoundaryConditions
+ constant | constant | constant
+%
+
+%MaxwellAbsorbingBoundaries
+ not_absorbing | not_absorbing | not_absorbing
+%
+
+# ----- Maxwell box variables -------------------------------------------------
+
+lsize_mx = 12.0
+dx_mx    = 0.5
+
+Maxwell.BoxShape   = parallelepiped
+
+%Maxwell.Lsize
+ lsize_mx | lsize_mx | lsize_mx
+%
+
+%Maxwell.Spacing
+ dx_mx | dx_mx | dx_mx
+%
+
+
+# ----- Output variables ------------------------------------------------------
+
+OutputFormat = plane_x + plane_y + plane_z + vtk + xyz + axis_x
+
+
+# ----- Maxwell output variables ----------------------------------------------
+
+MaxwellOutput = electric_field + magnetic_field + maxwell_energy_density
+
+MaxwellOutputInterval = 50
+MaxwellTDOutput       = maxwell_energy + maxwell_fields
+
+%MaxwellFieldsCoordinate
+  0.00 | 0.00 | 0.00
+%
+
+
+# ----- Time step variables ---------------------------------------------------
+
+TDTimeStep = 0.2
+TDMaxSteps = 10
+
+
+# ----- Spatially constant magnetic field -------------------------------------
+
+Ez           = 0.00000000
+By           = 0.010000000
+pulse_width  = 500.0
+pulse_shift  = 270.0
+pulse_slope  = 100.0
+
+# Column 1: constant electric field component in x-direction
+# Column 2: constant electric field component in y-direction
+# Column 3: constant electric field component in z-direction
+# Column 4: constant magnetic field component in x-direction
+# Column 5: constant magnetic field component in y-direction
+# Column 6: constant magnetic field component in z-direction
+# Column 7: name of the td function
+
+%UserDefinedConstantSpatialMaxwellField
+0 | 0 | Ez | 0 | By | 0 | "time_function"
+%
+
+PropagateSpatialMaxwellField = no
+
+# Column 1: name of the td function
+# Column 2: function
+# Column 3: amplitude
+# Column 4: pulse slope
+# Column 5: pulse width
+# Column 6: pulse shift
+
+%TDFunctions
+"time_function" | tdf_cw | 1.0
+%

testsuite/multisystem/05-lorentz-force.test

+# -*- coding: utf-8 mode: shell-script -*-
+
+Test       : Multisystem Propagation Lorentz force
+Program    : octopus
+TestGroups : short-run, multisystem
+Enabled    : Yes
+
+Processors: 1
+
+Input: 05-lorentz-force.01-charged_particle_coupled_to_maxwell.inp
+
+match ; Charged particle pos x (t=10 steps) ; LINEFIELD(ChargedParticle/td.general/coordinates, -1, 3) ; 9.998560012040e-01
+match ; Charged particle pos y (t=10 steps) ; LINEFIELD(ChargedParticle/td.general/coordinates, -1, 4) ; 0.000000000000e+00
+match ; Charged particle pos z (t=10 steps) ; LINEFIELD(ChargedParticle/td.general/coordinates, -1, 5) ; 1.999944000145e-02
+match ; Charged particle vel x (t=10 steps) ; LINEFIELD(ChargedParticle/td.general/coordinates, -1, 6) ; -1.699962300068e-04
+match ; Charged particle vel y (t=10 steps) ; LINEFIELD(ChargedParticle/td.general/coordinates, -1, 7) ; 0.000000000000e+00
+match ; Charged particle vel z (t=10 steps) ; LINEFIELD(ChargedParticle/td.general/coordinates, -1, 8) ; 9.998870006810e-03