Particle Steppers

Within the particle-in-cell capabilities in SALAMANDER, two main particle steppers are provided, the Leapfrog stepper and the Boris stepper. Both of these particle steppers are second order accurate in time. In order to verify that these steppers work properly, a convergence study can be utilized to ensure that second order accuracy is observed.

Demonstration Problems

To verify and demonstrate the implementation of the particle steppers, a classic example of single particle motion, cyclotron motion, is utilized. While cyclotron motion requires a magnetic field, we will utilize the analytic solution for particle motion to configure electric fields to reproduce the same particle paths.

Cyclotron Motion

The problem setup of cyclotron motion starts with a static magnetic field perpendicular to the plane in which the particle is traveling. For this example, the particle is traveling within the xy plane, and the magentic field is aligned along the positive z direction:

$var element = document.getElementById("moose-equation-5847a04f-222f-4753-b6ca-e2bea5e4b823");katex.render(" \\vec{B} \\left(\\vec{r}, t \\right) = B_0 \\hat{z},", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-5701a9ad-0f48-4144-9a67-2fe58ceaeed2");katex.render("B", element, {displayMode:false,throwOnError:false});$ represents the magnetic field and $var element = document.getElementById("moose-equation-b8a58333-b781-439c-a92d-b8fd3d9935e2");katex.render("B_0", element, {displayMode:false,throwOnError:false});$ denotes the magnitude of the magnetic field. Starting from Newton's equation of motion

$var element = document.getElementById("moose-equation-ef2cd04f-68f9-44b0-adb7-5ea4582999f0");katex.render(" \\frac{d\\vec{v}}{dt} = \\frac{q}{m} \\vec{v} \\times \\vec{B},", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-0c585a2c-67ad-474b-867a-b9f1e33f9535");katex.render("\\vec{v}", element, {displayMode:false,throwOnError:false});$ is the particle's velocity, $var element = document.getElementById("moose-equation-b2357f10-9db8-4e19-bbee-660aea3ee569");katex.render("q", element, {displayMode:false,throwOnError:false});$ is the particle's charge, and $var element = document.getElementById("moose-equation-ab3aee33-f18b-47ff-8671-6581156bfe22");katex.render("m", element, {displayMode:false,throwOnError:false});$ is the particle's mass.

The analytic solution to these equations are given by

$var element = document.getElementById("moose-equation-143ec12d-3020-425e-844c-f3088427b933");katex.render(" x(t) = r_L \\sin \\left( \\omega_c t \\right),", element, {displayMode:true,throwOnError:false});$ $var element = document.getElementById("moose-equation-d4aad9e1-8c2d-4b1c-8a4b-ee2516a83970");katex.render(" y(t) = r_L \\cos \\left( \\omega_c t \\right),", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-acf514a4-6672-43ed-8de8-54c9c2bec0ac");katex.render("r_L", element, {displayMode:false,throwOnError:false});$ , the Larmor radius, is given by

$var element = document.getElementById("moose-equation-152e3c42-def8-43dc-8d8a-136167717a17");katex.render(" r_L = \\frac{v_\\perp}{\\omega_c}.", element, {displayMode:true,throwOnError:false});$

Here, $var element = document.getElementById("moose-equation-11430a73-8bc0-4fd0-9d98-f29f5a775abf");katex.render("v_\\perp", element, {displayMode:false,throwOnError:false});$ is the initial speed of the particle in the plane perpendicular to the magnetic field, and $var element = document.getElementById("moose-equation-3d628772-b7ea-4358-bacf-93e4717e8fd9");katex.render("\\omega_c", element, {displayMode:false,throwOnError:false});$ is the gyro frequency. The gyrofrequency is defined as

$var element = document.getElementById("moose-equation-f0e45aab-f7da-4917-b64e-23b6fa4330a8");katex.render(" \\omega_c = \\frac{q B_0}{m}.", element, {displayMode:true,throwOnError:false});$

Time-varying Electric Field

The leapfrog stepper is designed to move particles through an electric field. Starting from the equations of motion for a particle in an electric field, the particle velocity is given by

$var element = document.getElementById("moose-equation-3f556790-a499-4813-aeb8-8725f83c6804");katex.render(" \\frac{d\\vec{v}}{dt} = \\frac{q}{m} \\vec{E}.", element, {displayMode:true,throwOnError:false});$

We can then use the solution for the particle's path under cyclotron motion to calculate the electric field which will reproduce the same particle motion. Doing so yields the electric field:

$var element = document.getElementById("moose-equation-77bd2454-df7e-48ab-941b-585b124f43bc");katex.render(" \\vec{E} \\left( \\vec{r}, t \\right) = - \\frac{m}{q} \\omega_c^2 r_L \\sin \\left( \\omega_c t \\right) \\hat{x} + - \\frac{m}{q} \\omega_c^2 r_L \\cos \\left( \\omega_c t \\right) \\hat{y}.", element, {displayMode:true,throwOnError:false});$

Convergence Studies

For the sake of simplicity, all particle properties are set to unity, in S.I. units, and the magnetic field magnitude, $var element = document.getElementById("moose-equation-d0b7c6b5-e1b4-41bb-8427-e94aeacba29d");katex.render("B_0", element, {displayMode:false,throwOnError:false});$ , is set to 1 T. Four time steps are used, each one an order of magnitude smaller than the last, $var element = document.getElementById("moose-equation-d2de47f2-c26a-4b1b-8661-0eb6fed0dfd9");katex.render("\\Delta t \\in [1, 10^{-1}, 10^{-2}, 10^{-3}]", element, {displayMode:false,throwOnError:false});$ seconds. After a simulation at a given time step is completed, the relative $var element = document.getElementById("moose-equation-a6ca82bf-e3b5-4e06-ab76-bb3e99d58c6a");katex.render("l_2", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-daccf5db-36ba-4ad8-94fb-54930e2b5340");katex.render("l_\\infty", element, {displayMode:false,throwOnError:false});$ norms of the error between the analytic and simulated solutions are taken. Plotting these errors on a log-log plot and then calculating the slope reveals the order of accuracy of the time integration; both implemented steppers are second order accurate in time, so the slope of each line should be two.

Boris Stepper: Cyclotron Motion

This test specifically tests the part of the stepper which calculates the force on the particle due to the magnetic field. The relative $var element = document.getElementById("moose-equation-3f318996-9c37-4410-b8f1-554ae2d3cc2d");katex.render("l_2", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-b9fd9034-74d2-4f75-9eee-ad33a809f74b");katex.render("l_\\infty", element, {displayMode:false,throwOnError:false});$ errors seen in Figure 1 and Figure 2, respectively, show the error over the course of the convergence study. Both errors indicate that the implementation of the Boris Stepper achieves the expected second order accuracy in time when moving a particle through a magnetic field. Second order accuracy is observed for both particle position and velocity.

Relative $l_2$ error for particle position and velocity for a particle subject to cyclotron motion being pushed via the Boris Stepper for four different time step sizes.

Figure 1: Relative $var element = document.getElementById("moose-equation-6abae716-a534-467f-ae88-ec429c914e51");katex.render("l_2", element, {displayMode:false,throwOnError:false});$ error for particle position and velocity for a particle subject to cyclotron motion being pushed via the Boris Stepper for four different time step sizes.

Relative $l_\infty$ error for particle position and velocity for a particle subject to cyclotron motion being pushed via the Boris Stepper for four different time step sizes.

Figure 2: Relative $var element = document.getElementById("moose-equation-2781949b-8631-45e1-a236-c09b136cc869");katex.render("l_\\infty", element, {displayMode:false,throwOnError:false});$ error for particle position and velocity for a particle subject to cyclotron motion being pushed via the Boris Stepper for four different time step sizes.

Boris Stepper: Time-varying Electric Field

In this example, the part of the Boris stepper which handles the acceleration on the particle due to the electric field is tested. Both the relative $var element = document.getElementById("moose-equation-358cbf2c-6abb-4219-a0f5-4bec4670bf99");katex.render("l_2", element, {displayMode:false,throwOnError:false});$ , Figure 3, and $var element = document.getElementById("moose-equation-6ba279be-f8bc-4c4e-aa10-ff797b3315ed");katex.render("l_\\infty", element, {displayMode:false,throwOnError:false});$ , Figure 4, errors are examined. The errors show that the implementation of the Boris stepper achieves the expected second order accuracy in time when moving a particle through an electric field.

Relative $l_2$ error for particle position and velocity for a particle being moved in a time varying electric field via the Boris Stepper for four different time step sizes.

Figure 3: Relative $var element = document.getElementById("moose-equation-a866d9e3-13ce-4a16-8b88-176ba06a2e86");katex.render("l_2", element, {displayMode:false,throwOnError:false});$ error for particle position and velocity for a particle being moved in a time varying electric field via the Boris Stepper for four different time step sizes.

Relative $l_\infty$ error for particle position and velocity for a particle being moved in a time varying electric field via the Boris Stepper for four different time step sizes.

Figure 4: Relative $var element = document.getElementById("moose-equation-998180f0-4294-4325-a796-f02a9b0ff554");katex.render("l_\\infty", element, {displayMode:false,throwOnError:false});$ error for particle position and velocity for a particle being moved in a time varying electric field via the Boris Stepper for four different time step sizes.

Leapfrog Stepper: Time-varying Electric Field

Since the Leapfrog stepper can be thought of as a special case of the Boris stepper where no magnetic field is present, the resulting errors should be the same as those produced in the previous test shown in Figure 3 and Figure 4. In fact, the errors and order of accuracy of the stepper implementation should be the same in both tests. In examining the relative $var element = document.getElementById("moose-equation-d95cbe68-080e-4c29-acc9-2e506b78154a");katex.render("l_2", element, {displayMode:false,throwOnError:false});$ , Figure 5, and $var element = document.getElementById("moose-equation-7c86218e-987f-4b8b-941c-8be5d6dce984");katex.render("l_\\infty", element, {displayMode:false,throwOnError:false});$ , Figure 6, errors for this test case, it can be seen that they are, in fact, the same as those observed in Figure 3 for $var element = document.getElementById("moose-equation-922d5362-24df-485e-abfa-55d2132957de");katex.render("l_2", element, {displayMode:false,throwOnError:false});$ and Figure 4 for $var element = document.getElementById("moose-equation-9e8d6b56-3cec-4617-a390-7deed8a21362");katex.render("l_\\infty", element, {displayMode:false,throwOnError:false});$ .

Relative $l_2$ error for particle position and velocity for a particle being moved in a time-varying electric field via the Leapfrog Stepper for four different time step sizes.

Figure 5: Relative $var element = document.getElementById("moose-equation-b189e24f-9d91-422f-b489-bf2f24461643");katex.render("l_2", element, {displayMode:false,throwOnError:false});$ error for particle position and velocity for a particle being moved in a time-varying electric field via the Leapfrog Stepper for four different time step sizes.

Relative $l_\infty$ error for particle position and velocity for a particle being moved in a time-varying electric field via the Leapfrog Stepper for four different time step sizes.

Figure 6: Relative $var element = document.getElementById("moose-equation-395abe5e-36be-4440-a934-1197a38c4e0f");katex.render("l_\\infty", element, {displayMode:false,throwOnError:false});$ error for particle position and velocity for a particle being moved in a time-varying electric field via the Leapfrog Stepper for four different time step sizes.

Through this verification test both the Boris and the Leapfrog steppers have been verified. Both portions of the Boris stepper, that which handles the acceleration due to a magnetic field and that which handles the acceleration due to an electric field, produce results which are, as expected, second order accurate in time. Additionally, the Leapfrog stepper produces these expected results as well.

Input Files

There are several input files used for this verification case. The two for the Boris stepper can be found at (test/tests/userobjects/particle_stepper/boris_stepper/cyclotron_motion.i) and (test/tests/userobjects/particle_stepper/boris_stepper/circular_e_field.i). The one for the Leapfrog stepper can be found at (test/tests/userobjects/particle_stepper/leapfrog_stepper/circular_e_field.i) Since the majority of the input file contents are the same for all cases, the files are broken up into several common files.

(test/tests/userobjects/particle_stepper/stepper_base.i) disables the field solver, so that particle motion tests are only conducted in fields that are known exactly, and sets up the PICStudy.
(test/tests/userobjects/particle_stepper/boris_stepper/boris_base.i) adds magnetic field components to the list of AuxVariables and selects the BorisStepper as the particle stepper.
(test/tests/userobjects/particle_stepper/leapfrog_stepper/leapfrog_base.i) selects the LeapFrogStepper as the particle stepper.
(test/tests/userobjects/particle_stepper/boris_stepper/cyclotron_motion.i), (test/tests/userobjects/particle_stepper/boris_stepper/circular_e_field.i), and (test/tests/userobjects/particle_stepper/leapfrog_stepper/circular_e_field.i) set up the initial conditions for the particles and fields for each case.

To combine them into one input file when running the simulation, the !include feature within MOOSE is utilized.

tip:Input file include syntax information

To learn more about the !include feature, refer to the Input File Syntax page.

(test/tests/userobjects/particle_stepper/boris_stepper/cyclotron_motion.i)

# The boris_base.i file sets up everything that the simulation needs to utilize
# the BorisStepper and sample the proper field variables.
!include boris_base.i
# This file gives the field variables the proper state for the cyclotron motion case,
# and creates a particle with the proper initial conditions.

[Mesh/gmg]
  nx = 5
  ny = 5
  xmin = -2
  xmax = 2
  ymin = -2
  ymax = 2
[]

[Functions]
  [E_x_ic]
    expression = '0'
  []
  [E_y_ic]
    expression = '0'
  []
  [E_z_ic]
    expression = '0'
  []
  [B_x_ic]
    expression = '0'
  []
  [B_y_ic]
    expression = '0'
  []
  [B_z_ic]
    expression = '1'
  []
[]

[UserObjects]
  [velocity_initializer]
    velocities = '1 0 0'
  []

  [particle_initializer]
    start_points = '0 1 0'
    mass = 1
    weight = 1
    charge = 1
  []
[]

[Executioner]
  dt = 1e-1
[]

(test/tests/userobjects/particle_stepper/boris_stepper/circular_e_field.i)

# The boris_base.i file sets up everything that the simulation needs to utilize
# the BorisStepper and sample the proper field variables.
!include boris_base.i
# This file gives the field variables the proper time-dependent values for the
# time varying electric field case, and creates a particle with the proper initial conditions.

[Mesh/gmg]
  nx = 5
  ny = 5
  xmin = -2
  xmax = 2
  ymin = -2
  ymax = 2
[]

[AuxKernels]
  [E_x]
    type = FunctionAux
    variable = Ex
    function = E_x_ic
  []
  [E_y]
    type = FunctionAux
    variable = Ey
    function = E_y_ic
  []
  [E_z]
    type = FunctionAux
    variable = Ez
    function = E_z_ic
  []
[]

[Functions]
  [E_x_ic]
    expression = '-sin(t)'
  []
  [E_y_ic]
    expression = '-cos(t)'
  []
  [E_z_ic]
    expression = '0'
  []
    [B_x_ic]
    expression = '0'
  []
  [B_y_ic]
    expression = '0'
  []
  [B_z_ic]
    expression = '0'
  []
[]

[UserObjects]
  [velocity_initializer]
    velocities = '1 0 0'
  []

  [particle_initializer]
    start_points = '0 1 0'
    mass = 1
    charge = 1
  []
[]


[Executioner]
  dt = 1e-1
[]

(test/tests/userobjects/particle_stepper/leapfrog_stepper/circular_e_field.i)

# The boris_base.i file sets up everything that the simulation needs to utilize
# the LeapFrogStepper and sample the proper field variables.
!include leapfrog_base.i
# This file gives the field variables the proper time-dependent values for the
# time varying electric field case, and creates a particle with the proper initial conditions.

[Mesh/gmg]
  nx = 5
  ny = 5
  xmin = -2
  xmax = 2
  ymin = -2
  ymax = 2
[]

[AuxKernels]
  [E_x]
    type = FunctionAux
    variable = Ex
    function = E_x_ic
  []
  [E_y]
    type = FunctionAux
    variable = Ey
    function = E_y_ic
  []
  [E_z]
    type = FunctionAux
    variable = Ez
    function = E_z_ic
  []
[]

[Functions]
  [E_x_ic]
    expression = '-sin(t)'
  []
  [E_y_ic]
    expression = '-cos(t)'
  []
  [E_z_ic]
    expression = '0'
  []
[]

[UserObjects]
  [velocity_initializer]
    velocities = '1 0 0'
  []

  [particle_initializer]
    start_points = '0 1 0'
    mass = 1
    charge = 1
  []
[]


[Executioner]
  dt = 1e-1
[]

(test/tests/userobjects/particle_stepper/stepper_base.i)

# This file serves as the common base case for all of the particle stepping test cases.
# It provides the bare minimum required to test a particle's motion.
# It also provides the required output capabilities to test a single particle's motion.
# The output method has been designed to properly output all particle data on time steps
# when a single particle is being used and the mesh is divided across several processors.
# Files that build off of this one will need to provide the following
# - The actual expression for the field initial conditions.
# - The type of stepper that will be tested
# - The initial conditions for the particle that will be utilized.
# Since this is only testing the particle stepper, we do not actually need to solve any
# field equations, and we are not going to be using any kernels.
[Problem]
  solve = false
  kernel_coverage_check=false
[]

[Mesh/gmg]
  type = GeneratedMeshGenerator
  dim = 2
[]

[AuxVariables]
  [Ex]
  []
  [Ey]
  []
  [Ez]
  []
[]

[Functions]
  [E_x_ic]
    type = ParsedFunction
  []
  [E_y_ic]
    type = ParsedFunction
  []
  [E_z_ic]
    type = ParsedFunction
  []
[]

[ICs]
  [Ex_ic]
    type = FunctionIC
    variable = Ex
    function = E_x_ic
  []
  [Ey_ic]
    type = FunctionIC
    variable = Ey
    function = E_y_ic
  []
  [Ez_ic]
    type = FunctionIC
    variable = Ez
    function = E_z_ic
  []
[]

[UserObjects]
  [velocity_initializer]
    type = ConstantVelocityInitializer
  []

  [stepper]
  []

  [particle_initializer]
    type = TestPlacedParticleInitializer
    velocity_initializer = 'velocity_initializer'
  []

  [study]
    type = TestInitializedPICStudy
    stepper = stepper
    particle_initializer = particle_initializer
    use_custom_rayids = false
    always_cache_traces = true
    data_on_cache_traces = true
    execute_on = 'TIMESTEP_BEGIN'
    ray_kernel_coverage_check = false
  []
[]


[VectorPostprocessors/particle_data]
  type = TestSingleParticleDataVectorPostprocessor
  study = study
  execute_on = 'TIMESTEP_END'
  additional_ray_data_outputs = 'charge mass'
[]

[Executioner]
  type = Transient
  num_steps = 10
[]

[Outputs]
  exodus = true
  [csv]
    type = CSV
    execute_on = 'FINAL'
  []
[]

(test/tests/userobjects/particle_stepper/boris_stepper/boris_base.i)

# The stepper base file provides the basic objects that are required for both
# the LeapFrogStepper and the BorisSteper tests
!include ../stepper_base.i
# The rest of this file adds the AuxVariables for the magnetic field, sets type of
# the particle stepper to be the BorisStepper and tells the stepper to sample the
# proper variables for the electric and magnetic fields.

[AuxVariables]
  [Bx]
  []
  [By]
  []
  [Bz]
  []
[]

[Functions]
  [B_x_ic]
    type = ParsedFunction
  []
  [B_y_ic]
    type = ParsedFunction
  []
  [B_z_ic]
    type = ParsedFunction
  []
[]

[ICs]
  [Bx_ic]
    type = FunctionIC
    variable = Bx
    function = B_x_ic
  []
  [By_ic]
    type = FunctionIC
    variable = By
    function = B_y_ic
  []
  [Bz_ic]
    type = FunctionIC
    variable = Bz
    function = B_z_ic
  []
[]

[UserObjects]
  [stepper]
    type = BorisStepper
    efield_components = 'Ex Ey Ez'
    bfield_components = 'Bx By Bz'
  []
[]

(test/tests/userobjects/particle_stepper/leapfrog_stepper/leapfrog_base.i)

# The stepper base file provides the basic objects that are required for both
# the LeapFrogStepper and the BorisSteper tests
!include ../stepper_base.i
# The rest of this file simply sets the type of stepper to be the LeapFrogStepper
# and tells the stepper which field variables to sample.

[UserObjects]
  [stepper]
    type = LeapFrogStepper
    field_components = 'Ex Ey Ez'
  []
[]

(test/tests/userobjects/particle_stepper/boris_stepper/cyclotron_motion.i)

# The boris_base.i file sets up everything that the simulation needs to utilize
# the BorisStepper and sample the proper field variables.
!include boris_base.i
# This file gives the field variables the proper state for the cyclotron motion case,
# and creates a particle with the proper initial conditions.

[Mesh/gmg]
  nx = 5
  ny = 5
  xmin = -2
  xmax = 2
  ymin = -2
  ymax = 2
[]

[Functions]
  [E_x_ic]
    expression = '0'
  []
  [E_y_ic]
    expression = '0'
  []
  [E_z_ic]
    expression = '0'
  []
  [B_x_ic]
    expression = '0'
  []
  [B_y_ic]
    expression = '0'
  []
  [B_z_ic]
    expression = '1'
  []
[]

[UserObjects]
  [velocity_initializer]
    velocities = '1 0 0'
  []

  [particle_initializer]
    start_points = '0 1 0'
    mass = 1
    weight = 1
    charge = 1
  []
[]

[Executioner]
  dt = 1e-1
[]

(test/tests/userobjects/particle_stepper/boris_stepper/circular_e_field.i)

# The boris_base.i file sets up everything that the simulation needs to utilize
# the BorisStepper and sample the proper field variables.
!include boris_base.i
# This file gives the field variables the proper time-dependent values for the
# time varying electric field case, and creates a particle with the proper initial conditions.

[Mesh/gmg]
  nx = 5
  ny = 5
  xmin = -2
  xmax = 2
  ymin = -2
  ymax = 2
[]

[AuxKernels]
  [E_x]
    type = FunctionAux
    variable = Ex
    function = E_x_ic
  []
  [E_y]
    type = FunctionAux
    variable = Ey
    function = E_y_ic
  []
  [E_z]
    type = FunctionAux
    variable = Ez
    function = E_z_ic
  []
[]

[Functions]
  [E_x_ic]
    expression = '-sin(t)'
  []
  [E_y_ic]
    expression = '-cos(t)'
  []
  [E_z_ic]
    expression = '0'
  []
    [B_x_ic]
    expression = '0'
  []
  [B_y_ic]
    expression = '0'
  []
  [B_z_ic]
    expression = '0'
  []
[]

[UserObjects]
  [velocity_initializer]
    velocities = '1 0 0'
  []

  [particle_initializer]
    start_points = '0 1 0'
    mass = 1
    charge = 1
  []
[]


[Executioner]
  dt = 1e-1
[]

(test/tests/userobjects/particle_stepper/leapfrog_stepper/circular_e_field.i)

# The boris_base.i file sets up everything that the simulation needs to utilize
# the LeapFrogStepper and sample the proper field variables.
!include leapfrog_base.i
# This file gives the field variables the proper time-dependent values for the
# time varying electric field case, and creates a particle with the proper initial conditions.

[Mesh/gmg]
  nx = 5
  ny = 5
  xmin = -2
  xmax = 2
  ymin = -2
  ymax = 2
[]

[AuxKernels]
  [E_x]
    type = FunctionAux
    variable = Ex
    function = E_x_ic
  []
  [E_y]
    type = FunctionAux
    variable = Ey
    function = E_y_ic
  []
  [E_z]
    type = FunctionAux
    variable = Ez
    function = E_z_ic
  []
[]

[Functions]
  [E_x_ic]
    expression = '-sin(t)'
  []
  [E_y_ic]
    expression = '-cos(t)'
  []
  [E_z_ic]
    expression = '0'
  []
[]

[UserObjects]
  [velocity_initializer]
    velocities = '1 0 0'
  []

  [particle_initializer]
    start_points = '0 1 0'
    mass = 1
    charge = 1
  []
[]


[Executioner]
  dt = 1e-1
[]

Cyclotron Motion
Time - varying Electric Field
Convergence Studies
Input Files