diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index abf44f6..e5bdc79 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.9.4","generation_timestamp":"2024-01-16T19:33:05","documenter_version":"1.2.1"}} \ No newline at end of file +{"documenter":{"julia_version":"1.9.4","generation_timestamp":"2024-01-16T20:01:24","documenter_version":"1.2.1"}} \ No newline at end of file diff --git a/dev/examples/electromagnetic/beam_heating/index.html b/dev/examples/electromagnetic/beam_heating/index.html index 29bde16..4b00684 100644 --- a/dev/examples/electromagnetic/beam_heating/index.html +++ b/dev/examples/electromagnetic/beam_heating/index.html @@ -1,2 +1,2 @@ -Beam heating · ParticleInCell.jl Documentation
GitHub
+Beam heating · ParticleInCell.jl Documentation
GitHub
diff --git a/dev/examples/electrostatic/beam_cyclotron/index.html b/dev/examples/electrostatic/beam_cyclotron/index.html index 51dc3e3..c31e991 100644 --- a/dev/examples/electrostatic/beam_cyclotron/index.html +++ b/dev/examples/electrostatic/beam_cyclotron/index.html @@ -1,2 +1,2 @@ -Beam-cyclotron instability · ParticleInCell.jl Documentation
GitHub
+Beam-cyclotron instability · ParticleInCell.jl Documentation
GitHub
diff --git a/dev/examples/electrostatic/beam_plasma/index.html b/dev/examples/electrostatic/beam_plasma/index.html index b91d377..14da4ca 100644 --- a/dev/examples/electrostatic/beam_plasma/index.html +++ b/dev/examples/electrostatic/beam_plasma/index.html @@ -1,2 +1,2 @@ -Beam-plasma instability · ParticleInCell.jl Documentation
GitHub
+Beam-plasma instability · ParticleInCell.jl Documentation
GitHub
diff --git a/dev/examples/electrostatic/landau_damping/index.html b/dev/examples/electrostatic/landau_damping/index.html index e76b07b..6f09a44 100644 --- a/dev/examples/electrostatic/landau_damping/index.html +++ b/dev/examples/electrostatic/landau_damping/index.html @@ -1,2 +1,2 @@ -Landau damping · ParticleInCell.jl Documentation
GitHub
+Landau damping · ParticleInCell.jl Documentation
GitHub
diff --git a/dev/examples/electrostatic/two_stream/index.html b/dev/examples/electrostatic/two_stream/index.html index b8654b9..e097827 100644 --- a/dev/examples/electrostatic/two_stream/index.html +++ b/dev/examples/electrostatic/two_stream/index.html @@ -1,2 +1,2 @@ -Two-stream instability · ParticleInCell.jl Documentation
GitHub
+Two-stream instability · ParticleInCell.jl Documentation
GitHub
diff --git a/dev/examples/index.html b/dev/examples/index.html index 1e8bae4..936940d 100644 --- a/dev/examples/index.html +++ b/dev/examples/index.html @@ -41,4 +41,4 @@

Coming soon...

- + diff --git a/dev/examples/tutorial/langmuir_oscillation.ipynb b/dev/examples/tutorial/langmuir_oscillation.ipynb index b3c4353..e9dbd04 100644 --- a/dev/examples/tutorial/langmuir_oscillation.ipynb +++ b/dev/examples/tutorial/langmuir_oscillation.ipynb @@ -3640,7 +3640,7 @@ { "output_type": "execute_result", "data": { - "text/plain": "4-element Vector{Float64}:\n 5.654866776461627e8\n 5.654866776461627e8\n 3.1415926535897934e8\n 1.2566370614359173e8" + "text/plain": "4-element Vector{Float64}:\n 5.654866776461627e8\n 3.7699111843077517e8\n 5.654866776461627e8\n 5.654866776461627e8" }, "metadata": {}, "execution_count": 21 diff --git a/dev/examples/tutorial/langmuir_oscillation/index.html b/dev/examples/tutorial/langmuir_oscillation/index.html index 7f18881..916abcd 100644 --- a/dev/examples/tutorial/langmuir_oscillation/index.html +++ b/dev/examples/tutorial/langmuir_oscillation/index.html @@ -175,8 +175,8 @@ measure_plasma_frequency.(1e14, temperatures, 1 / 2 * pi)
4-element Vector{Float64}:
  5.654866776461627e8
  5.654866776461627e8
- 5.654866776461627e8
- 5.654866776461627e8

We can compare this to the expected result.

function warm_plasma_freq(number_density, temperature, wavenumber)
+ 4.3982297150257105e8
+ 6.2831853071795866e7

We can compare this to the expected result.

function warm_plasma_freq(number_density, temperature, wavenumber)
     epsilon_0 = 8.8e-12
     elec_charge = 1.6e-19
     elec_mass = 9e-31
@@ -189,4 +189,4 @@
  5.685352436149611e8
  5.685352436182908e8
  5.685352436482581e8
- 5.685352439479299e8

Clearly, the restoring force of the pressure is not large enough to make a difference in this case.

This page was generated using DemoCards.jl and Literate.jl.

+ 5.685352439479299e8

Clearly, the restoring force of the pressure is not large enough to make a difference in this case.

This page was generated using DemoCards.jl and Literate.jl.

diff --git a/dev/index.html b/dev/index.html index 2cf25fe..5a7137f 100644 --- a/dev/index.html +++ b/dev/index.html @@ -1,3 +1,3 @@ Introduction · ParticleInCell.jl Documentation
GitHub

ParticleInCell.jl

ParticleInCell — Module

ParticleInCell.jl Logo

Latest Documentation CI Status Code Coverage Statistics

ParticleInCell.jl is a Julia package for kinetic plasma physics simulation. Specifically, it focuses on the simulation of kinetic (non-thermal) plasmas using particle-in-cell (PIC) algorithms. Currently, this package is in in a pre-1.0.0 state, and thus breaking changes should be expected. However, this also means that I am willing to entertain radical suggestions to improve the functionality of the package. If you are interested in using ParticleInCell for your plasma research, and you find that it does not meet you needs, please reach out on either GitHub, or over email, so that we can discuss how the package can be modified to suite your needs.

Getting Started

ParticleInCell is currently not registered in the Julia package registry. Thus, to install this package, you should use Pkg.develop:

using Pkg
-Pkg.develop(url="https://github.com/JuliaPlasma/ParticleInCell.jl")

Documentation

You can view the latest documentation here.

Goals

  • Fast: aim to have core time of less than 1 microsecond per particle per step without collisions.
  • Flexible: it should be possible to implement essentially any kinetic plasma simulation in ParticleInCell.jl. For common types of simulations, this might mean just piecing together components that are already included. More esoteric problems might require writing custom types that implement the desired algorithms. The advantage of writing this package in Julia is that these custom types will be just as performant as native components that are included in the package.
  • Scalable: the eventual goal is to enable scaling across an essentially unlimited number of cores using Julia's native multithreading for parallelization on a single node, and MPI.jl for communication across nodes. The goal is to support two different modes of parallelization:
    • Each core is responsible for a single rectangular subdomain. The domain assigned to an entire node is also rectangular, which imposes constraints on how the node domain can be subdivided into subdomains for each core. Load balancing is achieved by varying the relative sizes of the domains such that each core has a similar amount of work per step.
    • The simulation domain is subdivided into subdomains called 'patches', and every node is assigned a list of patches that it is responsible for updating. The cores on each node work collaboratively on the list, each choosing one patch to work on, and then selecting another when they are finished. Load balancing is achieved by swapping patches between nodes to balance the workload while also seeking to minimize communication time by keeping the surface area of each node's responsibilities minimized. In order for this scheme to effectively load balance, it must be the case that the total number of patches is larger (ideally much larger) that the total number of cores.
source

Documentation sections

To get started, look at the Tutorial, which includes step-by-step instructions for running your first simulation using ParticleInCell. After that, you may want to browse the Example Gallery, to see other problems that this package can be used to solve. The Plasma Simulation Theory section contains information about the art of computational plasma physics, and the numerical constraints that require the specialized tools developed here. Finally, the Manual contains detailed information about the entire public interface of the package.

+Pkg.develop(url="https://github.com/JuliaPlasma/ParticleInCell.jl")

Documentation

You can view the latest documentation here.

Goals

source

Documentation sections

To get started, look at the Tutorial, which includes step-by-step instructions for running your first simulation using ParticleInCell. After that, you may want to browse the Example Gallery, to see other problems that this package can be used to solve. The Plasma Simulation Theory section contains information about the art of computational plasma physics, and the numerical constraints that require the specialized tools developed here. Finally, the Manual contains detailed information about the entire public interface of the package.

diff --git a/dev/manual/index.html b/dev/manual/index.html index 3671882..a53f5f6 100644 --- a/dev/manual/index.html +++ b/dev/manual/index.html @@ -1,5 +1,5 @@ -Manual · ParticleInCell.jl Documentation
GitHub

Manual

Grids

ParticleInCell.AbstractGrid — Type

Parent type for all grid objects, which are used to define the simulation domain, and to convert between coordinate systems. There are three different numbering systems that can refer to a location in the simulation domain:

  1. The 'physical coordinates' of a point are the real (dimensionalful) coordinates associated with that point. This value can range from the lower bounds to the upper bounds of the simulation. This value will typically take the form Vector{T} or NTuple{N, T} where T <: Real.
  2. The 'cell coordinates' of a point is the non-dimensional location of the point in units of cell lengths. This value can range from 0 to num_cells - 1, or outside this range if guard cells are included. The value will typically have the type NTuple{N, Int}.
  3. The 'cell index' of a point is the CartesianIndex that can be used to index into field arrays at that point. This value must strictly be confined to axes(field.values), which, for any given dimension, will typically range from 1 to numcells + 2*numguard_cells + 1.

The first two types of indexing, physcoords and cellcoords, are independent of the number of guard cells in a given field, and depend only on grid quantities. Thus utilities for converting between these systems require only a reference to a grid object. On the other hand, the utilities for cell_index are specific the field being used, and thus those must be provided an AbstractField to do the coordinate conversion.

In general, physical coordinates are useful when considering the location of a particle, while the cell index is used to interpolate to and from the particle locations. The cell coordinates are useful for some interpolation algorithms, especially those that are defined for non-uniform grids.

source
ParticleInCell.UniformCartesianGrid — Type
UniformCartesianGrid(lower_bounds, upper_bounds, num_cells, periodic)

The simplest grid type, which represents a set of equally spaced rectangular cells. The grid can optionally be periodic in one or more dimensions.

source

Grid utility functions

ParticleInCell.cell_lengths — Function
cell_lengths(grid, [cell_coords])

Returns the length of the cell located at cell_coords. For uniform grids, the cell_coords argument is optional.

source
ParticleInCell.cell_coords_to_phys_coords — Function
cell_coords_to_phys_coords(grid, idxs, [offset, component])

Converts the cell coordinates idxs to a physical coordinate using the geometry specified in grid. Optionally, an offset and component can be specified to get the physical coordinates of a specific edge or face of the cell. The component argument is one-indexed.

For more information on the different types of coordinate systems, see AbstractGrid.

source
ParticleInCell.phys_coords_to_cell_coords — Function
phys_coords_to_cell_coords(grid, xs)

Converts the physical coordinate xs to a grid coordinate using the geometry specified in grid.

For more information on the different types of coordinate systems, see AbstractGrid.

source

Species

ParticleInCell.AbstractSpecies — Type

Subtypes of AbstractSpecies represent a group of macroparticles with a single physical charge and mass. Each particle of a species can optionally represent a different number of physical particles.

source
ParticleInCell.VariableWeightSpecies — Type
VariableWeightSpecies(positions::Matrix, momentums::Matrix, weight::Vector charge, mass)

Stores the positions and momentums of particles that share a common charge and mass. Each particle can have a different weight, and the momentum space can have a larger dimension than the position space.

source
ParticleInCell.electrons — Function
electrons(positions::Matrix, momentums::Matrix, weights::Vector)
-electrons(positions::Matrix, momentums::Matrix, weight)

Creates an AbsractSpecies with the given positions, momentums, and weight, and the charge and mass of a single electron.

source

Species utility functions

Base.eachindex — Function
eachindex(species)

Returns an iterator to the index of each of the macroparticles contained in species.

source

Fields

ParticleInCell.Field — Type
Field(grid, offset, num_elements,[ lower_guard_cells = 0,
-    [upper_guard_cells = lower_guard_cells + 1]])

Stores the values of a field.

source

Update steps

ParticleInCell.AbstractUpdateStep
-step!
ParticleInCell.ElectrostaticParticlePush — Type
ElectrostaticParticlePush(species, E, timestep, [interpolation_order=1])

An update step that moves and accelerates particles of species according to the electric field E. The field will be interpolated to the particle positions using b-splines of interpolation_order.

source
ParticleInCell.BorisParticlePush — Type
BorisParticlePush(species, E, B, timestep)

This update step moves the particles of species subject to both an electric field, E, and a magnetic field, B. The method is frequently referred to by the shorthand accelerate-rotate-accelerate because the acceleration from the electric field is split in half, and applied before and after the magnetic field rotation.

An applied magnetic field forces charged particles to travel in circular orbits perpendicular to the magnetic field. Thus, a simulation using the Boris method only makes sense when there are at least two velocity components, and in this case, the applied magnetic field must be strictly perpendicular to the velocity components. So for example, one could have a 1d2v simulation with a spatial grid along the x-axis, and velocity components in the x and y directions. In this case, the magnetic field would be forced to point along the z axis. If all three velocity components are included in the simulation, then the magnetic field can point in any direction.

For more details on the method, see Sections 4.3 and 4.4 of Birdsall and Langdon, or the Boris' original conference proceedings.

source

Misc

ParticleInCell.compute_bspline_coeffs — Function
compute_bspline_coeffs(degree, [T])

Returns a vector of vectors of the coefficients for the b-spline with polynomial degree, and unit-spaced knots, centered at zero.

source
+Manual · ParticleInCell.jl Documentation
GitHub

Manual

Grids

ParticleInCell.AbstractGrid — Type

Parent type for all grid objects, which are used to define the simulation domain, and to convert between coordinate systems. There are three different numbering systems that can refer to a location in the simulation domain:

  1. The 'physical coordinates' of a point are the real (dimensionalful) coordinates associated with that point. This value can range from the lower bounds to the upper bounds of the simulation. This value will typically take the form Vector{T} or NTuple{N, T} where T <: Real.
  2. The 'cell coordinates' of a point is the non-dimensional location of the point in units of cell lengths. This value can range from 0 to num_cells - 1, or outside this range if guard cells are included. The value will typically have the type NTuple{N, Int}.
  3. The 'cell index' of a point is the CartesianIndex that can be used to index into field arrays at that point. This value must strictly be confined to axes(field.values), which, for any given dimension, will typically range from 1 to numcells + 2*numguard_cells + 1.

The first two types of indexing, physcoords and cellcoords, are independent of the number of guard cells in a given field, and depend only on grid quantities. Thus utilities for converting between these systems require only a reference to a grid object. On the other hand, the utilities for cell_index are specific the field being used, and thus those must be provided an AbstractField to do the coordinate conversion.

In general, physical coordinates are useful when considering the location of a particle, while the cell index is used to interpolate to and from the particle locations. The cell coordinates are useful for some interpolation algorithms, especially those that are defined for non-uniform grids.

source
ParticleInCell.UniformCartesianGrid — Type
UniformCartesianGrid(lower_bounds, upper_bounds, num_cells, periodic)

The simplest grid type, which represents a set of equally spaced rectangular cells. The grid can optionally be periodic in one or more dimensions.

source

Grid utility functions

ParticleInCell.cell_lengths — Function
cell_lengths(grid, [cell_coords])

Returns the length of the cell located at cell_coords. For uniform grids, the cell_coords argument is optional.

source
ParticleInCell.cell_coords_to_phys_coords — Function
cell_coords_to_phys_coords(grid, idxs, [offset, component])

Converts the cell coordinates idxs to a physical coordinate using the geometry specified in grid. Optionally, an offset and component can be specified to get the physical coordinates of a specific edge or face of the cell. The component argument is one-indexed.

For more information on the different types of coordinate systems, see AbstractGrid.

source
ParticleInCell.phys_coords_to_cell_coords — Function
phys_coords_to_cell_coords(grid, xs)

Converts the physical coordinate xs to a grid coordinate using the geometry specified in grid.

For more information on the different types of coordinate systems, see AbstractGrid.

source

Species

ParticleInCell.AbstractSpecies — Type

Subtypes of AbstractSpecies represent a group of macroparticles with a single physical charge and mass. Each particle of a species can optionally represent a different number of physical particles.

source
ParticleInCell.VariableWeightSpecies — Type
VariableWeightSpecies(positions::Matrix, momentums::Matrix, weight::Vector charge, mass)

Stores the positions and momentums of particles that share a common charge and mass. Each particle can have a different weight, and the momentum space can have a larger dimension than the position space.

source
ParticleInCell.electrons — Function
electrons(positions::Matrix, momentums::Matrix, weights::Vector)
+electrons(positions::Matrix, momentums::Matrix, weight)

Creates an AbsractSpecies with the given positions, momentums, and weight, and the charge and mass of a single electron.

source

Species utility functions

Base.eachindex — Function
eachindex(species)

Returns an iterator to the index of each of the macroparticles contained in species.

source

Fields

ParticleInCell.Field — Type
Field(grid, offset, num_elements,[ lower_guard_cells = 0,
+    [upper_guard_cells = lower_guard_cells + 1]])

Stores the values of a field.

source

Update steps

ParticleInCell.AbstractUpdateStep
+step!
ParticleInCell.ElectrostaticParticlePush — Type
ElectrostaticParticlePush(species, E, timestep, [interpolation_order=1])

An update step that moves and accelerates particles of species according to the electric field E. The field will be interpolated to the particle positions using b-splines of interpolation_order.

source
ParticleInCell.BorisParticlePush — Type
BorisParticlePush(species, E, B, timestep)

This update step moves the particles of species subject to both an electric field, E, and a magnetic field, B. The method is frequently referred to by the shorthand accelerate-rotate-accelerate because the acceleration from the electric field is split in half, and applied before and after the magnetic field rotation.

An applied magnetic field forces charged particles to travel in circular orbits perpendicular to the magnetic field. Thus, a simulation using the Boris method only makes sense when there are at least two velocity components, and in this case, the applied magnetic field must be strictly perpendicular to the velocity components. So for example, one could have a 1d2v simulation with a spatial grid along the x-axis, and velocity components in the x and y directions. In this case, the magnetic field would be forced to point along the z axis. If all three velocity components are included in the simulation, then the magnetic field can point in any direction.

For more details on the method, see Sections 4.3 and 4.4 of Birdsall and Langdon, or the Boris' original conference proceedings.

source

Misc

ParticleInCell.compute_bspline_coeffs — Function
compute_bspline_coeffs(degree, [T])

Returns a vector of vectors of the coefficients for the b-spline with polynomial degree, and unit-spaced knots, centered at zero.

source
diff --git a/dev/references/index.html b/dev/references/index.html index 6480d53..a30cf40 100644 --- a/dev/references/index.html +++ b/dev/references/index.html @@ -1,2 +1,2 @@ -References · ParticleInCell.jl Documentation
GitHub

References

  • Birdsall, C. K. and Langdon, A. B. (2004). Plasma Physics via Computer Simulation. 1 edition Edition (CRC Press).
  • Boris, J. P. (1970). Relativistic plasma simulation-optimization of a hybrid code. Proceeding of Fourth Conference on Numerical Simulations of Plasmas, 3–8.
  • Hockney, R. W. and Eastwood, J. W. (1989). Computer Simulation Using Particles. First Edition Edition (CRC Press).
+References · ParticleInCell.jl Documentation
GitHub

References

  • Birdsall, C. K. and Langdon, A. B. (2004). Plasma Physics via Computer Simulation. 1 edition Edition (CRC Press).
  • Boris, J. P. (1970). Relativistic plasma simulation-optimization of a hybrid code. Proceeding of Fourth Conference on Numerical Simulations of Plasmas, 3–8.
  • Hockney, R. W. and Eastwood, J. W. (1989). Computer Simulation Using Particles. First Edition Edition (CRC Press).
diff --git a/dev/search_index.js b/dev/search_index.js index 158eef0..1d8fa36 100644 --- a/dev/search_index.js +++ b/dev/search_index.js @@ -1,3 +1,3 @@ var documenterSearchIndex = {"docs": -[{"location":"references/#References","page":"References","title":"References","text":"","category":"section"},{"location":"references/","page":"References","title":"References","text":"Birdsall, C. K. and Langdon, A. B. (2004). Plasma Physics via Computer Simulation. 1 edition Edition (CRC Press).\n\n\n\nBoris, J. P. (1970). Relativistic plasma simulation-optimization of a hybrid code. Proceeding of Fourth Conference on Numerical Simulations of Plasmas, 3–8.\n\n\n\nHockney, R. W. and Eastwood, J. W. (1989). Computer Simulation Using Particles. First Edition Edition (CRC Press).\n\n\n\n","category":"page"},{"location":"examples/electrostatic/landau_damping/","page":"Landau damping","title":"Landau damping","text":"EditURL = \"/home/runner/work/ParticleInCell.jl/ParticleInCell.jl/examples/electrostatic/landau_damping.jl\"","category":"page"},{"location":"examples/electrostatic/landau_damping/#Landau-damping","page":"Landau damping","title":"Landau damping","text":"","category":"section"},{"location":"examples/electrostatic/landau_damping/","page":"Landau damping","title":"Landau damping","text":"(Image: Source code) (Image: notebook)","category":"page"},{"location":"examples/electrostatic/landau_damping/","page":"Landau damping","title":"Landau damping","text":"Coming soon...","category":"page"},{"location":"examples/electrostatic/landau_damping/","page":"Landau damping","title":"Landau damping","text":"","category":"page"},{"location":"examples/electrostatic/landau_damping/","page":"Landau damping","title":"Landau damping","text":"This page was generated using DemoCards.jl and Literate.jl.","category":"page"},{"location":"manual/#Manual","page":"Manual","title":"Manual","text":"","category":"section"},{"location":"manual/#Grids","page":"Manual","title":"Grids","text":"","category":"section"},{"location":"manual/","page":"Manual","title":"Manual","text":"AbstractGrid\nUniformCartesianGrid","category":"page"},{"location":"manual/#ParticleInCell.AbstractGrid","page":"Manual","title":"ParticleInCell.AbstractGrid","text":"Parent type for all grid objects, which are used to define the simulation domain, and to convert between coordinate systems. There are three different numbering systems that can refer to a location in the simulation domain:\n\nThe 'physical coordinates' of a point are the real (dimensionalful) coordinates associated with that point. This value can range from the lower bounds to the upper bounds of the simulation. This value will typically take the form Vector{T} or NTuple{N, T} where T <: Real.\nThe 'cell coordinates' of a point is the non-dimensional location of the point in units of cell lengths. This value can range from 0 to num_cells - 1, or outside this range if guard cells are included. The value will typically have the type NTuple{N, Int}.\nThe 'cell index' of a point is the CartesianIndex that can be used to index into field arrays at that point. This value must strictly be confined to axes(field.values), which, for any given dimension, will typically range from 1 to numcells + 2*numguard_cells + 1.\n\nThe first two types of indexing, physcoords and cellcoords, are independent of the number of guard cells in a given field, and depend only on grid quantities. Thus utilities for converting between these systems require only a reference to a grid object. On the other hand, the utilities for cell_index are specific the field being used, and thus those must be provided an AbstractField to do the coordinate conversion.\n\nIn general, physical coordinates are useful when considering the location of a particle, while the cell index is used to interpolate to and from the particle locations. The cell coordinates are useful for some interpolation algorithms, especially those that are defined for non-uniform grids.\n\n\n\n\n\n","category":"type"},{"location":"manual/#ParticleInCell.UniformCartesianGrid","page":"Manual","title":"ParticleInCell.UniformCartesianGrid","text":"UniformCartesianGrid(lower_bounds, upper_bounds, num_cells, periodic)\n\nThe simplest grid type, which represents a set of equally spaced rectangular cells. The grid can optionally be periodic in one or more dimensions.\n\n\n\n\n\n","category":"type"},{"location":"manual/#Grid-utility-functions","page":"Manual","title":"Grid utility functions","text":"","category":"section"},{"location":"manual/","page":"Manual","title":"Manual","text":"ParticleInCell.cell_lengths\nParticleInCell.sim_lengths\nParticleInCell.cell_coords_to_phys_coords\nParticleInCell.phys_coords_to_cell_coords","category":"page"},{"location":"manual/#ParticleInCell.cell_lengths","page":"Manual","title":"ParticleInCell.cell_lengths","text":"cell_lengths(grid, [cell_coords])\n\nReturns the length of the cell located at cell_coords. For uniform grids, the cell_coords argument is optional.\n\n\n\n\n\n","category":"function"},{"location":"manual/#ParticleInCell.sim_lengths","page":"Manual","title":"ParticleInCell.sim_lengths","text":"sim_lengths(grid)\n\nReturns a tuple of the length of the simulation domain in each dimension.\n\n\n\n\n\n","category":"function"},{"location":"manual/#ParticleInCell.cell_coords_to_phys_coords","page":"Manual","title":"ParticleInCell.cell_coords_to_phys_coords","text":"cell_coords_to_phys_coords(grid, idxs, [offset, component])\n\nConverts the cell coordinates idxs to a physical coordinate using the geometry specified in grid. Optionally, an offset and component can be specified to get the physical coordinates of a specific edge or face of the cell. The component argument is one-indexed.\n\nFor more information on the different types of coordinate systems, see AbstractGrid.\n\n\n\n\n\n","category":"function"},{"location":"manual/#ParticleInCell.phys_coords_to_cell_coords","page":"Manual","title":"ParticleInCell.phys_coords_to_cell_coords","text":"phys_coords_to_cell_coords(grid, xs)\n\nConverts the physical coordinate xs to a grid coordinate using the geometry specified in grid.\n\nFor more information on the different types of coordinate systems, see AbstractGrid.\n\n\n\n\n\n","category":"function"},{"location":"manual/#Species","page":"Manual","title":"Species","text":"","category":"section"},{"location":"manual/","page":"Manual","title":"Manual","text":"AbstractSpecies\nVariableWeightSpecies\nelectrons","category":"page"},{"location":"manual/#ParticleInCell.AbstractSpecies","page":"Manual","title":"ParticleInCell.AbstractSpecies","text":"Subtypes of AbstractSpecies represent a group of macroparticles with a single physical charge and mass. Each particle of a species can optionally represent a different number of physical particles.\n\n\n\n\n\n","category":"type"},{"location":"manual/#ParticleInCell.VariableWeightSpecies","page":"Manual","title":"ParticleInCell.VariableWeightSpecies","text":"VariableWeightSpecies(positions::Matrix, momentums::Matrix, weight::Vector charge, mass)\n\nStores the positions and momentums of particles that share a common charge and mass. Each particle can have a different weight, and the momentum space can have a larger dimension than the position space.\n\n\n\n\n\n","category":"type"},{"location":"manual/#ParticleInCell.electrons","page":"Manual","title":"ParticleInCell.electrons","text":"electrons(positions::Matrix, momentums::Matrix, weights::Vector)\nelectrons(positions::Matrix, momentums::Matrix, weight)\n\nCreates an AbsractSpecies with the given positions, momentums, and weight, and the charge and mass of a single electron.\n\n\n\n\n\n","category":"function"},{"location":"manual/#Species-utility-functions","page":"Manual","title":"Species utility functions","text":"","category":"section"},{"location":"manual/","page":"Manual","title":"Manual","text":"particle_charge\nphysical_charge\nparticle_mass\nphysical_mass\nparticle_weight\nparticle_position\nparticle_position!\nparticle_momentum\nparticle_momentum!\nparticle_velocity\nphysical_momentum\neachindex","category":"page"},{"location":"manual/#ParticleInCell.particle_charge","page":"Manual","title":"ParticleInCell.particle_charge","text":"particle_charge(species, idx)\n\nReturns the charge of the macroparticle of species with index idx. See also physical_charge.\n\n\n\n\n\n","category":"function"},{"location":"manual/#ParticleInCell.physical_charge","page":"Manual","title":"ParticleInCell.physical_charge","text":"physical_charge(species, idx)\n\nReturns the charge of a physical particle represented by the macroparticle of species with index idx. See also particle_charge.\n\n\n\n\n\n","category":"function"},{"location":"manual/#ParticleInCell.particle_mass","page":"Manual","title":"ParticleInCell.particle_mass","text":"particle_mass(species, idx)\n\nReturns the mass of the macroparticle of species with index idx. See also physical_mass.\n\n\n\n\n\n","category":"function"},{"location":"manual/#ParticleInCell.physical_mass","page":"Manual","title":"ParticleInCell.physical_mass","text":"physical_mass(species, idx)\n\nReturns the mass of a physical particle represented by the macroparticle of species with index idx. See also particle_charge.\n\n\n\n\n\n","category":"function"},{"location":"manual/#ParticleInCell.particle_weight","page":"Manual","title":"ParticleInCell.particle_weight","text":"particle_weight(species, idx)\n\nReturns the number of physical particles represented by the macroparticle of species with index idx.\n\n\n\n\n\n","category":"function"},{"location":"manual/#ParticleInCell.particle_position","page":"Manual","title":"ParticleInCell.particle_position","text":"particle_position(species, idx)\n\nReturns the position of the macroparticle of species with index idx. See also particle_momentum and particle_position!.\n\n\n\n\n\n","category":"function"},{"location":"manual/#ParticleInCell.particle_position!","page":"Manual","title":"ParticleInCell.particle_position!","text":"particle_position!(species, idx, value)\n\nSets the position of the macroparticle of species with index idx to value. See also particle_momentum! and particle_position.\n\n\n\n\n\n","category":"function"},{"location":"manual/#ParticleInCell.particle_momentum","page":"Manual","title":"ParticleInCell.particle_momentum","text":"particle_momentum(species, idx)\n\nReturns the momentum of the macroparticle of species with index idx. See also particle_momentum!, particle_velocity, and physical_momentum.\n\n\n\n\n\n","category":"function"},{"location":"manual/#ParticleInCell.particle_momentum!","page":"Manual","title":"ParticleInCell.particle_momentum!","text":"particle_momentum!(species, idx, value)\n\nSets the momentum of the macroparticle of species with index idx to value. See also particle_momentum, particle_velocity, and physical_momentum.\n\n\n\n\n\n","category":"function"},{"location":"manual/#ParticleInCell.particle_velocity","page":"Manual","title":"ParticleInCell.particle_velocity","text":"particle_velocity(species, idx)\n\nReturns the velocity of the macroparticle of species with index idx. See also particle_momentum, particle_momentum!, and physical_momentum.\n\n\n\n\n\n","category":"function"},{"location":"manual/#ParticleInCell.physical_momentum","page":"Manual","title":"ParticleInCell.physical_momentum","text":"physical_momentum(species, idx)\n\nReturns the momentum of a physical particle represented by the macroparticle of species with index idx. See also particle_momentum, particle_momentum!, and particle_velocity.\n\n\n\n\n\n","category":"function"},{"location":"manual/#Base.eachindex","page":"Manual","title":"Base.eachindex","text":"eachindex(species)\n\nReturns an iterator to the index of each of the macroparticles contained in species.\n\n\n\n\n\n","category":"function"},{"location":"manual/#Fields","page":"Manual","title":"Fields","text":"","category":"section"},{"location":"manual/","page":"Manual","title":"Manual","text":"Field","category":"page"},{"location":"manual/#ParticleInCell.Field","page":"Manual","title":"ParticleInCell.Field","text":"Field(grid, offset, num_elements,[ lower_guard_cells = 0,\n [upper_guard_cells = lower_guard_cells + 1]])\n\nStores the values of a field.\n\n\n\n\n\n","category":"type"},{"location":"manual/#Update-steps","page":"Manual","title":"Update steps","text":"","category":"section"},{"location":"manual/","page":"Manual","title":"Manual","text":"ParticleInCell.AbstractUpdateStep\nstep!","category":"page"},{"location":"manual/","page":"Manual","title":"Manual","text":"ElectrostaticParticlePush\nBorisParticlePush","category":"page"},{"location":"manual/#ParticleInCell.ElectrostaticParticlePush","page":"Manual","title":"ParticleInCell.ElectrostaticParticlePush","text":"ElectrostaticParticlePush(species, E, timestep, [interpolation_order=1])\n\nAn update step that moves and accelerates particles of species according to the electric field E. The field will be interpolated to the particle positions using b-splines of interpolation_order.\n\n\n\n\n\n","category":"type"},{"location":"manual/#ParticleInCell.BorisParticlePush","page":"Manual","title":"ParticleInCell.BorisParticlePush","text":"BorisParticlePush(species, E, B, timestep)\n\nThis update step moves the particles of species subject to both an electric field, E, and a magnetic field, B. The method is frequently referred to by the shorthand accelerate-rotate-accelerate because the acceleration from the electric field is split in half, and applied before and after the magnetic field rotation.\n\nAn applied magnetic field forces charged particles to travel in circular orbits perpendicular to the magnetic field. Thus, a simulation using the Boris method only makes sense when there are at least two velocity components, and in this case, the applied magnetic field must be strictly perpendicular to the velocity components. So for example, one could have a 1d2v simulation with a spatial grid along the x-axis, and velocity components in the x and y directions. In this case, the magnetic field would be forced to point along the z axis. If all three velocity components are included in the simulation, then the magnetic field can point in any direction.\n\nFor more details on the method, see Sections 4.3 and 4.4 of Birdsall and Langdon, or the Boris' original conference proceedings.\n\n\n\n\n\n","category":"type"},{"location":"manual/#Misc","page":"Manual","title":"Misc","text":"","category":"section"},{"location":"manual/","page":"Manual","title":"Manual","text":"ParticleInCell.compute_knots\nParticleInCell.compute_bspline_coeffs","category":"page"},{"location":"manual/#ParticleInCell.compute_knots","page":"Manual","title":"ParticleInCell.compute_knots","text":"compute_knots(degree)\n\nReturns a vector of the knot locations for a polynomial b-spline with degree.\n\n\n\n\n\n","category":"function"},{"location":"manual/#ParticleInCell.compute_bspline_coeffs","page":"Manual","title":"ParticleInCell.compute_bspline_coeffs","text":"compute_bspline_coeffs(degree, [T])\n\nReturns a vector of vectors of the coefficients for the b-spline with polynomial degree, and unit-spaced knots, centered at zero.\n\n\n\n\n\n","category":"function"},{"location":"examples/electrostatic/beam_cyclotron/","page":"Beam-cyclotron instability","title":"Beam-cyclotron instability","text":"EditURL = \"/home/runner/work/ParticleInCell.jl/ParticleInCell.jl/examples/electrostatic/beam_cyclotron.jl\"","category":"page"},{"location":"examples/electrostatic/beam_cyclotron/#Beam-cyclotron-instability","page":"Beam-cyclotron instability","title":"Beam-cyclotron instability","text":"","category":"section"},{"location":"examples/electrostatic/beam_cyclotron/","page":"Beam-cyclotron instability","title":"Beam-cyclotron instability","text":"(Image: Source code) (Image: notebook)","category":"page"},{"location":"examples/electrostatic/beam_cyclotron/","page":"Beam-cyclotron instability","title":"Beam-cyclotron instability","text":"Coming soon...","category":"page"},{"location":"examples/electrostatic/beam_cyclotron/","page":"Beam-cyclotron instability","title":"Beam-cyclotron instability","text":"","category":"page"},{"location":"examples/electrostatic/beam_cyclotron/","page":"Beam-cyclotron instability","title":"Beam-cyclotron instability","text":"This page was generated using DemoCards.jl and Literate.jl.","category":"page"},{"location":"examples/electrostatic/two_stream/","page":"Two-stream instability","title":"Two-stream instability","text":"EditURL = \"/home/runner/work/ParticleInCell.jl/ParticleInCell.jl/examples/electrostatic/two_stream.jl\"","category":"page"},{"location":"examples/electrostatic/two_stream/#Two-stream-instability","page":"Two-stream instability","title":"Two-stream instability","text":"","category":"section"},{"location":"examples/electrostatic/two_stream/","page":"Two-stream instability","title":"Two-stream instability","text":"(Image: Source code) (Image: notebook)","category":"page"},{"location":"examples/electrostatic/two_stream/","page":"Two-stream instability","title":"Two-stream instability","text":"Coming soon...","category":"page"},{"location":"examples/electrostatic/two_stream/","page":"Two-stream instability","title":"Two-stream instability","text":"","category":"page"},{"location":"examples/electrostatic/two_stream/","page":"Two-stream instability","title":"Two-stream instability","text":"This page was generated using DemoCards.jl and Literate.jl.","category":"page"},{"location":"theory/intro_to_pic/#Introduction-to-particle-in-cell-simulation","page":"PIC simulation","title":"Introduction to particle-in-cell simulation","text":"","category":"section"},{"location":"theory/intro_to_pic/","page":"PIC simulation","title":"PIC simulation","text":"Classical plasma physics considers the motion of charged particles. The dynamics of these particles will be effected by the presence of externally imposed electric and magnetic fields, which is relatively easy to model since the motion of each particle is independent. However, the particles will themselves source an electric field–-and, if they are moving fast enough, a magnetic field–-due to their charge. The dynamics of other particles will be influenced by these 'self-consistent' fields, which corresponding source fields of their own. Thus, the dynamics of all of the particles is coupled, and so their equations of motion must be solve together. For a typical plasma, the resulting differential equations cannot be solved analytically, and the number of degrees of freedom means that direct computational integration of the equations is also impossible. Thus, plasma physics relies on various simplifications and assumptions to create reduced models that can be solved–- either exactly, or approximately.","category":"page"},{"location":"theory/intro_to_pic/#Approximations-to-the-true-particle-distribution-function","page":"PIC simulation","title":"Approximations to the true particle distribution function","text":"","category":"section"},{"location":"theory/intro_to_pic/","page":"PIC simulation","title":"PIC simulation","text":"It is convenient to represent the locations and momenta of the particles using the distribution function","category":"page"},{"location":"theory/intro_to_pic/","page":"PIC simulation","title":"PIC simulation","text":"f(xpt) = sum_i=1^N delta(x - x_i(t)) delta(p - p_i(t))","category":"page"},{"location":"theory/intro_to_pic/","page":"PIC simulation","title":"PIC simulation","text":"where N is the total number of particles. The evolution of this distribution function can be written using the Maxwell-Boltzmann system","category":"page"},{"location":"theory/intro_to_pic/","page":"PIC simulation","title":"PIC simulation","text":"beginaligned\nTODO\nTODO-Maxwell\nendaligned","category":"page"},{"location":"theory/intro_to_pic/","page":"PIC simulation","title":"PIC simulation","text":"Unfortunately, for almost all plasmas of interest, N is enormous; thus direct simulation of the equations of motion is computationally intractable. The solution is to recognize that if N is large, then any physically small phase-space region will contain many particles, and so we c","category":"page"},{"location":"theory/intro_to_pic/#Todo","page":"PIC simulation","title":"Todo","text":"","category":"section"},{"location":"theory/intro_to_pic/","page":"PIC simulation","title":"PIC simulation","text":"course-graining\ncollisionless plasmas\nbriefly mention thermalization, two-fluid, and MHD","category":"page"},{"location":"theory/intro_to_pic/#Discretized-solutions-of-the-Boltzmann-Poisson-equations","page":"PIC simulation","title":"Discretized solutions of the Boltzmann-Poisson equations","text":"","category":"section"},{"location":"theory/intro_to_pic/","page":"PIC simulation","title":"PIC simulation","text":"We have seen that a kinetic plasma can be approximated using the Vlasov-Boltzmann equation, along with an appropriate equation of motion for the electromagnetic fields. However, the resulting equation is still not easy to analyse for an arbitrary f. We therefore turn to computational simulation of the equations of motion to understand the plasma dynamics.","category":"page"},{"location":"theory/intro_to_pic/","page":"PIC simulation","title":"PIC simulation","text":"In order to simulate the plasma, we must first discretize the equations so they can be represented on a computer. One option, called a Vlasov code, represents the distribution function as","category":"page"},{"location":"theory/intro_to_pic/#The-standard-PIC-cycle","page":"PIC simulation","title":"The standard PIC cycle","text":"","category":"section"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"EditURL = \"/home/runner/work/ParticleInCell.jl/ParticleInCell.jl/examples/index.md\"","category":"page"},{"location":"examples/#Example-Gallery","page":"Example Gallery","title":"Example Gallery","text":"","category":"section"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"These examples are in progress...","category":"page"},{"location":"examples/#Tutorial","page":"Example Gallery","title":"Tutorial","text":"","category":"section"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n\n \n \n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"(Image: list-card-cover-image)","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"Tutorial: Langmuir oscillations","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"In this tutorial, you will use ParticleInCell to model one of the simplest phenomena in plasma physics: an electrostatic (or Langmuir) oscillation. This tutorial is part of a series of examples that uses ParticleInCell to demonstrate the basic plasma physics concepts that are covered in Birdsall and Langdon's classic PIC textbook.","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n
\n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
","category":"page"},{"location":"examples/#Electrostatic","page":"Example Gallery","title":"Electrostatic","text":"","category":"section"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n\n \n \n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"(Image: list-card-cover-image)","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"Two-stream instability","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"Coming soon...","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n
\n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n\n \n \n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"(Image: list-card-cover-image)","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"Beam-plasma instability","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"Coming soon...","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n
\n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n\n \n \n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"(Image: list-card-cover-image)","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"Beam-cyclotron instability","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"Coming soon...","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n
\n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n\n \n \n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"(Image: list-card-cover-image)","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"Landau damping","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"Coming soon...","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n
\n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
","category":"page"},{"location":"examples/#Electromagnetic","page":"Example Gallery","title":"Electromagnetic","text":"","category":"section"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n\n \n \n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"(Image: list-card-cover-image)","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"Beam heating","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"Coming soon...","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n
\n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
","category":"page"},{"location":"examples/electromagnetic/beam_heating/","page":"Beam heating","title":"Beam heating","text":"EditURL = \"/home/runner/work/ParticleInCell.jl/ParticleInCell.jl/examples/electromagnetic/beam_heating.jl\"","category":"page"},{"location":"examples/electromagnetic/beam_heating/#Beam-heating","page":"Beam heating","title":"Beam heating","text":"","category":"section"},{"location":"examples/electromagnetic/beam_heating/","page":"Beam heating","title":"Beam heating","text":"(Image: Source code) (Image: notebook)","category":"page"},{"location":"examples/electromagnetic/beam_heating/","page":"Beam heating","title":"Beam heating","text":"Coming soon...","category":"page"},{"location":"examples/electromagnetic/beam_heating/","page":"Beam heating","title":"Beam heating","text":"","category":"page"},{"location":"examples/electromagnetic/beam_heating/","page":"Beam heating","title":"Beam heating","text":"This page was generated using DemoCards.jl and Literate.jl.","category":"page"},{"location":"theory/#Plasma-Simulation-Theory","page":"Introduction","title":"Plasma Simulation Theory","text":"","category":"section"},{"location":"theory/","page":"Introduction","title":"Introduction","text":"The fundamental physics governing the dynamics of a plasma have been well understood for over a century, and yet plasma physics remains an active area of research. This is because the dynamics of a plasma are highly nonlinear, and it is therefore difficult to make analytic statements about how a given plasma will behave over long time-spans. Instead, theoretical plasma physicists often rely on simulation to understand plasma dynamics, and to make general statements about the behaviors of particular classes of plasmas.","category":"page"},{"location":"theory/","page":"Introduction","title":"Introduction","text":"Unfortunately, the simulation of plasmas is itself a quite challenging task because plasmas are composed of many, many, charged particles. As a result, there are far too many degrees of freedom to exactly solve the full equations of motion for a given plasma. Instead, physicists rely on approximations to derive physically relevant models for the systems in question.","category":"page"},{"location":"theory/","page":"Introduction","title":"Introduction","text":"The most \"realistic\" class of models are called kinetic models. In these models, the individual charged particles of each species are averaged to create distribution functions that depends on both position and velocity. The distribution functions give the likelihood of finding a charged particle of a particular species at any point in phase space. These distribution functions, along with a method for calculating the self-consistent interaction between the particles, yields a set of approximate equations of motion describing the plasma dynamics.","category":"page"},{"location":"theory/","page":"Introduction","title":"Introduction","text":"If the particle species are nonrelativistic, then the particles are not influenced by self-consistent magnetic fields, and so the interactions can be modeled using electrostatics. However, once the particles (typically the lightweight electrons) become relativistic, the sourced magnetic field must be taken into account, and so full electromagnetic interactions are required.","category":"page"},{"location":"theory/","page":"Introduction","title":"Introduction","text":"Over long periods of time, the plasma will begin to thermalize–-the distribution of particle velocities will become closer and closer to Maxwellian. This fact can be used to drastically improve the size and speed of simulations by using the two-fluid and magnetohydrodynamic (MHD) approximations of plasma dynamics. As this package does not implement algorithms for these simulation methods, we will not discuss the details of these methods further.","category":"page"},{"location":"theory/","page":"Introduction","title":"Introduction","text":"In the following sections, we describe the details of particle-in-cell algorithms, a specific type of kinetic simulation algorithm. For a more in depth introduction of PIC simulation theory, see the books by Birdsall and Langdon (2004) and Hockney and Eastwood (1989).","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"EditURL = \"/home/runner/work/ParticleInCell.jl/ParticleInCell.jl/examples/tutorial/langmuir_oscillation.jl\"","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/#tutorial_langmuir","page":"Tutorial","title":"Tutorial: Langmuir oscillations","text":"","category":"section"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"(Image: Source code) (Image: notebook)","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"In this tutorial, you will use ParticleInCell to model one of the simplest phenomena in plasma physics: an electrostatic (or Langmuir) oscillation. This tutorial is part of a series of examples that uses ParticleInCell to demonstrate the basic plasma physics concepts that are covered in Birdsall and Langdon's classic PIC textbook.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"A Langmuir oscillation occurs when a slab of charge in a uniform plasma is displaced. The resulting charge density gradient creates a restoring force that causes the displaced slab of charge to return to its original position. But–-just as in a classical pendulum oscillation–-the momentum of the charge carries it past its equilibrium point, creating an opposite charge gradient, and a restoring force in the opposite direction. As a result, the slab of charge oscillates around its equilibrium forever (at least in this idealized model that ignores possible damping mechanisms). For a plasma composed of a single mobile species s with mass m_s and charge q_s, the frequency of this oscillation is given by","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"omega_ps = sqrtfracn_s q_s^2epsilon_0 m_s","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"where n_s is the number density of the plasma and epsilon_0 is the permitivity of free space. Notice that the plasma frequency has a m_s^-12 dependence, and thus the lightest species (typically electrons) will dominate the dynamics of a plasma oscillation. For this reason, we will only model the dynamics of the electrons in our simulation.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/#Simulating-a-cold-electron-plasma","page":"Tutorial","title":"Simulating a cold electron plasma","text":"","category":"section"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"We begin by loading the ParticleInCell package. Additionally, we load CairoMakie which is a backend for Makie that can generate beautiful, publication-quality graphics.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"using ParticleInCell\nusing CairoMakie","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"We begin by creating some electrons to move in the simulation. For even a tiny simulation volume, there are far too many physical electrons to simulate each one individually. Instead, PIC algorithms group physical particles into 'macroparticles'. The distribution of macroparticles in phase space serves as an approximation for the phase space distribution of physical particles. We arbitrarily choose a simulation domain of length one, and a nominal electrons number density of 10^14. Then, for a given number of macroparticles, we can calculate the number of physical electrons represented by each.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"sim_length = 1.0\nnumber_density = 1e14\nnum_macroparticles = 320\nparticles_per_macro = number_density * sim_length / num_macroparticles","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"3.125e11","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"We distribute the macroparticles evenly across the simulation domain.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"positions = collect(0:num_macroparticles-1) ./ num_macroparticles;","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"In order to seed a Langmuir oscillation, we give the electrons a sinusoidal velocity perturbation. This corresponds to the moment in a Langmuir oscillation when the slab of charge has reached equilibrium, but is being carried past by its momentum. This perturbation is defined by a wavenumber k and an amplitude.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"k = 1 * 2pi / sim_length\namplitude = 1e3\nelec_mass = 9e-31\nmomentums = (particles_per_macro * elec_mass * amplitude) .* sin.(positions .* k);","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"We can visualize the initial condition of the electron macroparticle by plotting the initial phase space.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"scatter(\n positions,\n momentums;\n axis = (;\n title = \"Electron phase space\",\n xlabel = \"Position (m)\",\n ylabel = \"Momentum (kg m / s)\",\n ),\n)","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"(Image: )","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"Finally, we create a VariableWeightSpecies which holds the all of the macroparticles. Additionally, we must pass the value of particles_per_macro, which is used to calculate the charge and mass of the macroparticles.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"electrons = ParticleInCell.electrons(positions, momentums, particles_per_macro);","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"Now we address the 'cell' piece of particle-in-cell by creating a grid. Because Langmuir oscillations are a one-dimensional phenomena, we will choose to perform a 1D simulation.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"The choice of grid resolution is determined by the scale of the smallest relevant dynamics begin simulated. For a Langmuir oscillation, the scale of the dynamics is set by k, and so the simulation could likely accomplished with as few as 4 or 8 cells. However, this is not a computationally demanding simulation, and so we arbitrarily choose to use 32 equally spaced (i.e. uniform) grid points. Additionally, we make the simulation domain periodic.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"num_cells = 32\ndx = sim_length / num_cells\nperiodic = true\ngrid = UniformCartesianGrid((0.0,), (sim_length,), (num_cells,), (periodic,));","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"Next, we set up the required fields for an electrostatic PIC simulation. In a basic PIC cycle, we first compute the charge density, rho, on the grid points. We then compute the corresponding electric potential, phi. The electric field is conventionally determined in a two step process. First the potential, which is located at the nodes of the grid cells, is finite differenced to the edges of the cells, producing an edge electric field, Eedge. The edge electric fields are then averaged to get the electric fields located at the nodes, Enode.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"When creating the fields, we must specify the underlying grid on which the field is based, the location of the values of the field (i.e. are the field values located at the nodes of each cell? The edge?), the dimension of the field, and the number of guard cells surrounding the field.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"field_dimension = 1\nlower_guard_cells = 1\nrho = Field(grid, ParticleInCell.node, field_dimension, lower_guard_cells)\nphi = Field(grid, ParticleInCell.node, field_dimension, lower_guard_cells)\nEedge = Field(grid, ParticleInCell.edge, field_dimension, lower_guard_cells)\nEnode = Field(grid, ParticleInCell.node, field_dimension, lower_guard_cells);","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"At this point, we must choose a timestep for the simulation. We would like to use a large timestep so that more of the systems dynamics can be observed with the same number of steps. However, we must resolve the fastest timescale of the dynamics that we are trying to simulate. In this case, we must resolve the plasma frequency. Additionally, we must choose a timestep that is short enough that particles do not cross more than one cell per timestep to prevent numerical instabilities from arising. For the oscillation amplitude that we have chosen, the particles do not move fast enough for the CFL condition to matter, and so we will choose our timestep based on the expected plasma frequency.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"epsilon_0 = 8.8e-12\nelec_charge = 1.6e-19\nelec_mass = 9e-31\nexpected_plasma_freq = sqrt(number_density * elec_charge^2 / elec_mass / epsilon_0)\nexpected_plasma_period = 2pi / expected_plasma_freq","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"1.1051531770007306e-8","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"Once again, this is not a computationally demanding simulation, and so we will choose a relatively small timestep for improve numerical accuracy. You can play with increasing the timestep, and see when the simulation results begin to deteriorate.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"dt = 5e-11","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"5.0e-11","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"In the final step of the setup, we create all of the simulation steps required to do the electrostatic simulation. In this tutorial, we will not discuss the details of PIC simulation, but you can find more information about the PIC simulation cycle elsewhere in this documentation.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"charge_interp = BSplineChargeInterpolation(electrons, rho, 1)\ncomm_rho = CommunicateGuardCells(rho, true)\nfield_solve = PoissonSolveFFT(rho, phi)\ncomm_phi = CommunicateGuardCells(phi)\nfinite_diff = FiniteDifferenceToEdges(phi, Eedge)\ncomm_Eedge = CommunicateGuardCells(Eedge)\nelec_ave = AverageEdgesToNodes(Eedge, Enode)\ncomm_Enode = CommunicateGuardCells(Enode)\npush = ElectrostaticParticlePush(electrons, Enode, dt)\ncomm_electrons = CommunicateSpecies(electrons, grid);","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"Now we are ready to run the simulation. We will simulate the plasma for 1000 timesteps, and at each step, we will calculate the electric field energy,","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"U_E = int E(x)^2 mathrmdx","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"This field energy will oscillate as the electrons move in and out of equilibrium, and so we can use it to observe the Langmuir oscillation.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"n_steps = 1000\n\nelectric_field_energy = Vector{Float64}(undef, n_steps)\n\nfor n = 1:n_steps\n # Calculate the electric field energy\n electric_field_energy[n] = 0\n for I in eachindex(Enode)\n electric_field_energy[n] += (dx * epsilon_0 / 2) * (Enode.values[I])^2\n end\n\n # TODO\n rho.values .= 0\n\n step!(charge_interp)\n step!(comm_rho)\n step!(field_solve)\n step!(comm_phi)\n step!(finite_diff)\n step!(comm_Eedge)\n step!(elec_ave)\n step!(comm_Enode)\n step!(push)\n step!(comm_electrons)\nend","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"We can now visualize the electric field energy to see the plasma oscillation.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"times = collect(range(1, n_steps)) .* dt\nlines(\n times,\n electric_field_energy;\n axis = (; title = \"Electric Field Energy\", xlabel = \"Time (s)\", ylabel = \"Energy\"),\n)","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"(Image: )","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"Notice that the electric field energy is slowly growing over time, which is unphysical. We will discuss where this numerical instability comes from– and how it can be avoided–later. But for now, we can still use the electric-field-energy time series to calculate the plasma frequency. First, let's plot the Fourier transform of the electric field energy.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"using FFTW\n\nfreqs = fftfreq(n_steps, 1 / dt) .* 2pi\nfreq_amps = abs.(fft(electric_field_energy))\n\nlines(\n freqs,\n freq_amps;\n axis = (;\n title = \"Electric Field Energy Frequency Spectrum\",\n xlabel = \"Frequency (1/s)\",\n ylabel = \"Amplitude\",\n ),\n)","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"(Image: )","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"It is hard to see what is happening at the low frequencies, so let's zoom in on the positive low frequencies.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"cutoff_index = round(Int, n_steps * 0.05)\nlines(\n freqs[1:cutoff_index],\n freq_amps[1:cutoff_index];\n axis = (;\n title = \"Electric Field Energy Frequency Spectrum\",\n xlabel = \"Frequency (1/s)\",\n ylabel = \"Amplitude\",\n ),\n)","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"(Image: )","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"Next, we find the maximum frequency. We don't care about the spike at zero frequency (that is just a consequence of the electric field energy being a strictly positive quantity) so we first set its amplitude to zero, and then find the largest remaining amplitude, and it's corresponding frequency.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"freq_amps[1] = 0\nmax_index = findmax(freq_amps)[2]\nmax_freq = freqs[max_index]\n\n# Divide by 2 because the electric field energy goes through a maximum twice\n# per plasma oscillation, and take the absolute value because we don't care\n# about the phase of the oscillation.\nplasma_freq = abs(max_freq / 2)","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"5.654866776461627e8","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"Finally, we can compare this to the theoretically expected result:","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"relative_error = (plasma_freq - expected_plasma_freq) / expected_plasma_freq","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"-0.005362140699342522","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"Less than 1% error. Not bad!","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/#Adding-temperature-to-the-plasma","page":"Tutorial","title":"Adding temperature to the plasma","text":"","category":"section"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"In the previous simulation, the electric field energy grew unphysically throughout the simulation. This was a result of the \"grid-heating instability\", which occurs when the grid does not resolve the plasma Debye length, which for an electron plasma is given by","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"lambda_De = sqrtfracepsilon_0 k_B Tn_e q_e^2","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"where k_B is the Boltzmann constant and T_e is the electron temperature.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"When the Debye length of a plasma is underresolved, the plasma will unphysically heat, causing the Debye length to grow until it is resolved by the grid. Many strategies have been developed to mitigate this instability, but in this tutorial, we will simply give our plasma sufficient thermal energy to begin so that the simulation will be stable against the grid-heating instability.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"Let's calculate the required electron temperature in the previous simulation so that the 32 cell grid will resolve the Debye length. We set lambda_De = Delta x, and solve for T to find","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"boltzmann_constant = 1.381e-23\ndx^2 * number_density * elec_charge^2 / epsilon_0 / boltzmann_constant","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"2.0571390955170825e7","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"Alternatively, we can express this temperature in terms of electron volts as","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"dx^2 * number_density * elec_charge / epsilon_0","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"1775.5681818181818","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"We will define a function that takes an electron density, electron temperature, and oscillation wavenumber, and returns a measured plasma frequency.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"function measure_plasma_frequency(number_density, temperature, wavenumber)\n sim_length = 1.0\n\n num_cells = 32\n dx = sim_length / num_cells\n\n num_macroparticles = 320\n particles_per_macro = number_density * sim_length / num_macroparticles\n\n perturb_amplitude = 1e3\n elec_mass = 9e-31\n boltzmann_constant = 1.381e-23\n thermal_velocity = sqrt(3 * boltzmann_constant * temperature / elec_mass)\n\n positions = collect(0:num_macroparticles-1) ./ num_macroparticles\n momentums =\n (particles_per_macro * elec_mass) .*\n (perturb_amplitude .* sin.(positions .* wavenumber) .+ thermal_velocity .* randn.())\n\n electrons = ParticleInCell.electrons(positions, momentums, particles_per_macro)\n\n grid = UniformCartesianGrid((0.0,), (sim_length,), (num_cells,), (true,))\n\n field_dimension = 1\n lower_guard_cells = 1\n rho = Field(grid, ParticleInCell.node, field_dimension, lower_guard_cells)\n phi = Field(grid, ParticleInCell.node, field_dimension, lower_guard_cells)\n Eedge = Field(grid, ParticleInCell.edge, field_dimension, lower_guard_cells)\n Enode = Field(grid, ParticleInCell.node, field_dimension, lower_guard_cells)\n\n dt = 5e-11\n charge_interp = BSplineChargeInterpolation(electrons, rho, 1)\n comm_rho = CommunicateGuardCells(rho, true)\n field_solve = PoissonSolveFFT(rho, phi)\n comm_phi = CommunicateGuardCells(phi)\n finite_diff = FiniteDifferenceToEdges(phi, Eedge)\n comm_Eedge = CommunicateGuardCells(Eedge)\n elec_ave = AverageEdgesToNodes(Eedge, Enode)\n comm_Enode = CommunicateGuardCells(Enode)\n push = ElectrostaticParticlePush(electrons, Enode, dt)\n comm_electrons = CommunicateSpecies(electrons, grid)\n\n n_steps = 1000\n electric_field_energy = Vector{Float64}(undef, n_steps)\n\n epsilon_0 = 8.8e-12\n for n = 1:n_steps\n # Calculate the electric field energy\n electric_field_energy[n] = 0\n for I in eachindex(Enode)\n electric_field_energy[n] += (dx * epsilon_0 / 2) * (Enode.values[I])^2\n end\n\n # TODO\n rho.values .= 0\n\n step!(charge_interp)\n step!(comm_rho)\n step!(field_solve)\n step!(comm_phi)\n step!(finite_diff)\n step!(comm_Eedge)\n step!(elec_ave)\n step!(comm_Enode)\n step!(push)\n step!(comm_electrons)\n end\n\n freqs = fftfreq(n_steps, 1 / dt) .* 2pi\n freq_amps = abs.(fft(electric_field_energy))\n\n freq_amps[1] = 0\n max_index = findmax(freq_amps)[2]\n max_freq = freqs[max_index]\n plasma_freq = abs(max_freq / 2)\n\n return plasma_freq\nend","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"measure_plasma_frequency (generic function with 1 method)","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"For a warm plasma, the thermal pressure acts as an additional restoring force on the plasma oscillation. It can be show that the modified dispersion relation (frequency response as a function of wavenumber) is given by","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"omega^2 = omega_pe^2 + frac3 k_B T_em_e k^2","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"Notice that when T_e = 0, the frequency is the plasma frequency, regardless of the wavenumber.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"Let's run a few simulations and verify that this relationship holds.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"temperatures = [0, 0.1, 1, 10]\nmeasure_plasma_frequency.(1e14, temperatures, 1 / 2 * pi)","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"4-element Vector{Float64}:\n 5.654866776461627e8\n 5.654866776461627e8\n 5.654866776461627e8\n 5.654866776461627e8","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"We can compare this to the expected result.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"function warm_plasma_freq(number_density, temperature, wavenumber)\n epsilon_0 = 8.8e-12\n elec_charge = 1.6e-19\n elec_mass = 9e-31\n boltzmann_constant = 1.381e-23\n square_plasma_freq = number_density * elec_charge^2 / elec_mass / epsilon_0\n correction_factor = boltzmann_constant * temperature * wavenumber^2 / elec_mass\n return sqrt(square_plasma_freq + correction_factor)\nend\nwarm_plasma_freq.(1e14, temperatures, 1 / 2 * pi)","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"4-element Vector{Float64}:\n 5.685352436149611e8\n 5.685352436182908e8\n 5.685352436482581e8\n 5.685352439479299e8","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"Clearly, the restoring force of the pressure is not large enough to make a difference in this case.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"This page was generated using DemoCards.jl and Literate.jl.","category":"page"},{"location":"examples/electrostatic/beam_plasma/","page":"Beam-plasma instability","title":"Beam-plasma instability","text":"EditURL = \"/home/runner/work/ParticleInCell.jl/ParticleInCell.jl/examples/electrostatic/beam_plasma.jl\"","category":"page"},{"location":"examples/electrostatic/beam_plasma/#Beam-plasma-instability","page":"Beam-plasma instability","title":"Beam-plasma instability","text":"","category":"section"},{"location":"examples/electrostatic/beam_plasma/","page":"Beam-plasma instability","title":"Beam-plasma instability","text":"(Image: Source code) (Image: notebook)","category":"page"},{"location":"examples/electrostatic/beam_plasma/","page":"Beam-plasma instability","title":"Beam-plasma instability","text":"Coming soon...","category":"page"},{"location":"examples/electrostatic/beam_plasma/","page":"Beam-plasma instability","title":"Beam-plasma instability","text":"","category":"page"},{"location":"examples/electrostatic/beam_plasma/","page":"Beam-plasma instability","title":"Beam-plasma instability","text":"This page was generated using DemoCards.jl and Literate.jl.","category":"page"},{"location":"#ParticleInCell.jl","page":"Introduction","title":"ParticleInCell.jl","text":"","category":"section"},{"location":"","page":"Introduction","title":"Introduction","text":"ParticleInCell","category":"page"},{"location":"#ParticleInCell","page":"Introduction","title":"ParticleInCell","text":"(Image: ParticleInCell.jl Logo)\n\n(Image: Latest Documentation) (Image: CI Status) (Image: Code Coverage Statistics)\n\nParticleInCell.jl is a Julia package for kinetic plasma physics simulation. Specifically, it focuses on the simulation of kinetic (non-thermal) plasmas using particle-in-cell (PIC) algorithms. Currently, this package is in in a pre-1.0.0 state, and thus breaking changes should be expected. However, this also means that I am willing to entertain radical suggestions to improve the functionality of the package. If you are interested in using ParticleInCell for your plasma research, and you find that it does not meet you needs, please reach out on either GitHub, or over email, so that we can discuss how the package can be modified to suite your needs.\n\nGetting Started\n\nParticleInCell is currently not registered in the Julia package registry. Thus, to install this package, you should use Pkg.develop:\n\nusing Pkg\nPkg.develop(url=\"https://github.com/JuliaPlasma/ParticleInCell.jl\")\n\nDocumentation\n\nYou can view the latest documentation here.\n\nGoals\n\nFast: aim to have core time of less than 1 microsecond per particle per step without collisions.\nFlexible: it should be possible to implement essentially any kinetic plasma simulation in ParticleInCell.jl. For common types of simulations, this might mean just piecing together components that are already included. More esoteric problems might require writing custom types that implement the desired algorithms. The advantage of writing this package in Julia is that these custom types will be just as performant as native components that are included in the package.\nScalable: the eventual goal is to enable scaling across an essentially unlimited number of cores using Julia's native multithreading for parallelization on a single node, and MPI.jl for communication across nodes. The goal is to support two different modes of parallelization:\nEach core is responsible for a single rectangular subdomain. The domain assigned to an entire node is also rectangular, which imposes constraints on how the node domain can be subdivided into subdomains for each core. Load balancing is achieved by varying the relative sizes of the domains such that each core has a similar amount of work per step.\nThe simulation domain is subdivided into subdomains called 'patches', and every node is assigned a list of patches that it is responsible for updating. The cores on each node work collaboratively on the list, each choosing one patch to work on, and then selecting another when they are finished. Load balancing is achieved by swapping patches between nodes to balance the workload while also seeking to minimize communication time by keeping the surface area of each node's responsibilities minimized. In order for this scheme to effectively load balance, it must be the case that the total number of patches is larger (ideally much larger) that the total number of cores.\n\n\n\n\n\n","category":"module"},{"location":"#Documentation-sections","page":"Introduction","title":"Documentation sections","text":"","category":"section"},{"location":"","page":"Introduction","title":"Introduction","text":"To get started, look at the Tutorial, which includes step-by-step instructions for running your first simulation using ParticleInCell. After that, you may want to browse the Example Gallery, to see other problems that this package can be used to solve. The Plasma Simulation Theory section contains information about the art of computational plasma physics, and the numerical constraints that require the specialized tools developed here. Finally, the Manual contains detailed information about the entire public interface of the package.","category":"page"}] +[{"location":"references/#References","page":"References","title":"References","text":"","category":"section"},{"location":"references/","page":"References","title":"References","text":"Birdsall, C. K. and Langdon, A. B. (2004). Plasma Physics via Computer Simulation. 1 edition Edition (CRC Press).\n\n\n\nBoris, J. P. (1970). Relativistic plasma simulation-optimization of a hybrid code. Proceeding of Fourth Conference on Numerical Simulations of Plasmas, 3–8.\n\n\n\nHockney, R. W. and Eastwood, J. W. (1989). Computer Simulation Using Particles. First Edition Edition (CRC Press).\n\n\n\n","category":"page"},{"location":"examples/electrostatic/landau_damping/","page":"Landau damping","title":"Landau damping","text":"EditURL = \"/home/runner/work/ParticleInCell.jl/ParticleInCell.jl/examples/electrostatic/landau_damping.jl\"","category":"page"},{"location":"examples/electrostatic/landau_damping/#Landau-damping","page":"Landau damping","title":"Landau damping","text":"","category":"section"},{"location":"examples/electrostatic/landau_damping/","page":"Landau damping","title":"Landau damping","text":"(Image: Source code) (Image: notebook)","category":"page"},{"location":"examples/electrostatic/landau_damping/","page":"Landau damping","title":"Landau damping","text":"Coming soon...","category":"page"},{"location":"examples/electrostatic/landau_damping/","page":"Landau damping","title":"Landau damping","text":"","category":"page"},{"location":"examples/electrostatic/landau_damping/","page":"Landau damping","title":"Landau damping","text":"This page was generated using DemoCards.jl and Literate.jl.","category":"page"},{"location":"manual/#Manual","page":"Manual","title":"Manual","text":"","category":"section"},{"location":"manual/#Grids","page":"Manual","title":"Grids","text":"","category":"section"},{"location":"manual/","page":"Manual","title":"Manual","text":"AbstractGrid\nUniformCartesianGrid","category":"page"},{"location":"manual/#ParticleInCell.AbstractGrid","page":"Manual","title":"ParticleInCell.AbstractGrid","text":"Parent type for all grid objects, which are used to define the simulation domain, and to convert between coordinate systems. There are three different numbering systems that can refer to a location in the simulation domain:\n\nThe 'physical coordinates' of a point are the real (dimensionalful) coordinates associated with that point. This value can range from the lower bounds to the upper bounds of the simulation. This value will typically take the form Vector{T} or NTuple{N, T} where T <: Real.\nThe 'cell coordinates' of a point is the non-dimensional location of the point in units of cell lengths. This value can range from 0 to num_cells - 1, or outside this range if guard cells are included. The value will typically have the type NTuple{N, Int}.\nThe 'cell index' of a point is the CartesianIndex that can be used to index into field arrays at that point. This value must strictly be confined to axes(field.values), which, for any given dimension, will typically range from 1 to numcells + 2*numguard_cells + 1.\n\nThe first two types of indexing, physcoords and cellcoords, are independent of the number of guard cells in a given field, and depend only on grid quantities. Thus utilities for converting between these systems require only a reference to a grid object. On the other hand, the utilities for cell_index are specific the field being used, and thus those must be provided an AbstractField to do the coordinate conversion.\n\nIn general, physical coordinates are useful when considering the location of a particle, while the cell index is used to interpolate to and from the particle locations. The cell coordinates are useful for some interpolation algorithms, especially those that are defined for non-uniform grids.\n\n\n\n\n\n","category":"type"},{"location":"manual/#ParticleInCell.UniformCartesianGrid","page":"Manual","title":"ParticleInCell.UniformCartesianGrid","text":"UniformCartesianGrid(lower_bounds, upper_bounds, num_cells, periodic)\n\nThe simplest grid type, which represents a set of equally spaced rectangular cells. The grid can optionally be periodic in one or more dimensions.\n\n\n\n\n\n","category":"type"},{"location":"manual/#Grid-utility-functions","page":"Manual","title":"Grid utility functions","text":"","category":"section"},{"location":"manual/","page":"Manual","title":"Manual","text":"ParticleInCell.cell_lengths\nParticleInCell.sim_lengths\nParticleInCell.cell_coords_to_phys_coords\nParticleInCell.phys_coords_to_cell_coords","category":"page"},{"location":"manual/#ParticleInCell.cell_lengths","page":"Manual","title":"ParticleInCell.cell_lengths","text":"cell_lengths(grid, [cell_coords])\n\nReturns the length of the cell located at cell_coords. For uniform grids, the cell_coords argument is optional.\n\n\n\n\n\n","category":"function"},{"location":"manual/#ParticleInCell.sim_lengths","page":"Manual","title":"ParticleInCell.sim_lengths","text":"sim_lengths(grid)\n\nReturns a tuple of the length of the simulation domain in each dimension.\n\n\n\n\n\n","category":"function"},{"location":"manual/#ParticleInCell.cell_coords_to_phys_coords","page":"Manual","title":"ParticleInCell.cell_coords_to_phys_coords","text":"cell_coords_to_phys_coords(grid, idxs, [offset, component])\n\nConverts the cell coordinates idxs to a physical coordinate using the geometry specified in grid. Optionally, an offset and component can be specified to get the physical coordinates of a specific edge or face of the cell. The component argument is one-indexed.\n\nFor more information on the different types of coordinate systems, see AbstractGrid.\n\n\n\n\n\n","category":"function"},{"location":"manual/#ParticleInCell.phys_coords_to_cell_coords","page":"Manual","title":"ParticleInCell.phys_coords_to_cell_coords","text":"phys_coords_to_cell_coords(grid, xs)\n\nConverts the physical coordinate xs to a grid coordinate using the geometry specified in grid.\n\nFor more information on the different types of coordinate systems, see AbstractGrid.\n\n\n\n\n\n","category":"function"},{"location":"manual/#Species","page":"Manual","title":"Species","text":"","category":"section"},{"location":"manual/","page":"Manual","title":"Manual","text":"AbstractSpecies\nVariableWeightSpecies\nelectrons","category":"page"},{"location":"manual/#ParticleInCell.AbstractSpecies","page":"Manual","title":"ParticleInCell.AbstractSpecies","text":"Subtypes of AbstractSpecies represent a group of macroparticles with a single physical charge and mass. Each particle of a species can optionally represent a different number of physical particles.\n\n\n\n\n\n","category":"type"},{"location":"manual/#ParticleInCell.VariableWeightSpecies","page":"Manual","title":"ParticleInCell.VariableWeightSpecies","text":"VariableWeightSpecies(positions::Matrix, momentums::Matrix, weight::Vector charge, mass)\n\nStores the positions and momentums of particles that share a common charge and mass. Each particle can have a different weight, and the momentum space can have a larger dimension than the position space.\n\n\n\n\n\n","category":"type"},{"location":"manual/#ParticleInCell.electrons","page":"Manual","title":"ParticleInCell.electrons","text":"electrons(positions::Matrix, momentums::Matrix, weights::Vector)\nelectrons(positions::Matrix, momentums::Matrix, weight)\n\nCreates an AbsractSpecies with the given positions, momentums, and weight, and the charge and mass of a single electron.\n\n\n\n\n\n","category":"function"},{"location":"manual/#Species-utility-functions","page":"Manual","title":"Species utility functions","text":"","category":"section"},{"location":"manual/","page":"Manual","title":"Manual","text":"particle_charge\nphysical_charge\nparticle_mass\nphysical_mass\nparticle_weight\nparticle_position\nparticle_position!\nparticle_momentum\nparticle_momentum!\nparticle_velocity\nphysical_momentum\neachindex","category":"page"},{"location":"manual/#ParticleInCell.particle_charge","page":"Manual","title":"ParticleInCell.particle_charge","text":"particle_charge(species, idx)\n\nReturns the charge of the macroparticle of species with index idx. See also physical_charge.\n\n\n\n\n\n","category":"function"},{"location":"manual/#ParticleInCell.physical_charge","page":"Manual","title":"ParticleInCell.physical_charge","text":"physical_charge(species, idx)\n\nReturns the charge of a physical particle represented by the macroparticle of species with index idx. See also particle_charge.\n\n\n\n\n\n","category":"function"},{"location":"manual/#ParticleInCell.particle_mass","page":"Manual","title":"ParticleInCell.particle_mass","text":"particle_mass(species, idx)\n\nReturns the mass of the macroparticle of species with index idx. See also physical_mass.\n\n\n\n\n\n","category":"function"},{"location":"manual/#ParticleInCell.physical_mass","page":"Manual","title":"ParticleInCell.physical_mass","text":"physical_mass(species, idx)\n\nReturns the mass of a physical particle represented by the macroparticle of species with index idx. See also particle_charge.\n\n\n\n\n\n","category":"function"},{"location":"manual/#ParticleInCell.particle_weight","page":"Manual","title":"ParticleInCell.particle_weight","text":"particle_weight(species, idx)\n\nReturns the number of physical particles represented by the macroparticle of species with index idx.\n\n\n\n\n\n","category":"function"},{"location":"manual/#ParticleInCell.particle_position","page":"Manual","title":"ParticleInCell.particle_position","text":"particle_position(species, idx)\n\nReturns the position of the macroparticle of species with index idx. See also particle_momentum and particle_position!.\n\n\n\n\n\n","category":"function"},{"location":"manual/#ParticleInCell.particle_position!","page":"Manual","title":"ParticleInCell.particle_position!","text":"particle_position!(species, idx, value)\n\nSets the position of the macroparticle of species with index idx to value. See also particle_momentum! and particle_position.\n\n\n\n\n\n","category":"function"},{"location":"manual/#ParticleInCell.particle_momentum","page":"Manual","title":"ParticleInCell.particle_momentum","text":"particle_momentum(species, idx)\n\nReturns the momentum of the macroparticle of species with index idx. See also particle_momentum!, particle_velocity, and physical_momentum.\n\n\n\n\n\n","category":"function"},{"location":"manual/#ParticleInCell.particle_momentum!","page":"Manual","title":"ParticleInCell.particle_momentum!","text":"particle_momentum!(species, idx, value)\n\nSets the momentum of the macroparticle of species with index idx to value. See also particle_momentum, particle_velocity, and physical_momentum.\n\n\n\n\n\n","category":"function"},{"location":"manual/#ParticleInCell.particle_velocity","page":"Manual","title":"ParticleInCell.particle_velocity","text":"particle_velocity(species, idx)\n\nReturns the velocity of the macroparticle of species with index idx. See also particle_momentum, particle_momentum!, and physical_momentum.\n\n\n\n\n\n","category":"function"},{"location":"manual/#ParticleInCell.physical_momentum","page":"Manual","title":"ParticleInCell.physical_momentum","text":"physical_momentum(species, idx)\n\nReturns the momentum of a physical particle represented by the macroparticle of species with index idx. See also particle_momentum, particle_momentum!, and particle_velocity.\n\n\n\n\n\n","category":"function"},{"location":"manual/#Base.eachindex","page":"Manual","title":"Base.eachindex","text":"eachindex(species)\n\nReturns an iterator to the index of each of the macroparticles contained in species.\n\n\n\n\n\n","category":"function"},{"location":"manual/#Fields","page":"Manual","title":"Fields","text":"","category":"section"},{"location":"manual/","page":"Manual","title":"Manual","text":"Field","category":"page"},{"location":"manual/#ParticleInCell.Field","page":"Manual","title":"ParticleInCell.Field","text":"Field(grid, offset, num_elements,[ lower_guard_cells = 0,\n [upper_guard_cells = lower_guard_cells + 1]])\n\nStores the values of a field.\n\n\n\n\n\n","category":"type"},{"location":"manual/#Update-steps","page":"Manual","title":"Update steps","text":"","category":"section"},{"location":"manual/","page":"Manual","title":"Manual","text":"ParticleInCell.AbstractUpdateStep\nstep!","category":"page"},{"location":"manual/","page":"Manual","title":"Manual","text":"ElectrostaticParticlePush\nBorisParticlePush","category":"page"},{"location":"manual/#ParticleInCell.ElectrostaticParticlePush","page":"Manual","title":"ParticleInCell.ElectrostaticParticlePush","text":"ElectrostaticParticlePush(species, E, timestep, [interpolation_order=1])\n\nAn update step that moves and accelerates particles of species according to the electric field E. The field will be interpolated to the particle positions using b-splines of interpolation_order.\n\n\n\n\n\n","category":"type"},{"location":"manual/#ParticleInCell.BorisParticlePush","page":"Manual","title":"ParticleInCell.BorisParticlePush","text":"BorisParticlePush(species, E, B, timestep)\n\nThis update step moves the particles of species subject to both an electric field, E, and a magnetic field, B. The method is frequently referred to by the shorthand accelerate-rotate-accelerate because the acceleration from the electric field is split in half, and applied before and after the magnetic field rotation.\n\nAn applied magnetic field forces charged particles to travel in circular orbits perpendicular to the magnetic field. Thus, a simulation using the Boris method only makes sense when there are at least two velocity components, and in this case, the applied magnetic field must be strictly perpendicular to the velocity components. So for example, one could have a 1d2v simulation with a spatial grid along the x-axis, and velocity components in the x and y directions. In this case, the magnetic field would be forced to point along the z axis. If all three velocity components are included in the simulation, then the magnetic field can point in any direction.\n\nFor more details on the method, see Sections 4.3 and 4.4 of Birdsall and Langdon, or the Boris' original conference proceedings.\n\n\n\n\n\n","category":"type"},{"location":"manual/#Misc","page":"Manual","title":"Misc","text":"","category":"section"},{"location":"manual/","page":"Manual","title":"Manual","text":"ParticleInCell.compute_knots\nParticleInCell.compute_bspline_coeffs","category":"page"},{"location":"manual/#ParticleInCell.compute_knots","page":"Manual","title":"ParticleInCell.compute_knots","text":"compute_knots(degree)\n\nReturns a vector of the knot locations for a polynomial b-spline with degree.\n\n\n\n\n\n","category":"function"},{"location":"manual/#ParticleInCell.compute_bspline_coeffs","page":"Manual","title":"ParticleInCell.compute_bspline_coeffs","text":"compute_bspline_coeffs(degree, [T])\n\nReturns a vector of vectors of the coefficients for the b-spline with polynomial degree, and unit-spaced knots, centered at zero.\n\n\n\n\n\n","category":"function"},{"location":"examples/electrostatic/beam_cyclotron/","page":"Beam-cyclotron instability","title":"Beam-cyclotron instability","text":"EditURL = \"/home/runner/work/ParticleInCell.jl/ParticleInCell.jl/examples/electrostatic/beam_cyclotron.jl\"","category":"page"},{"location":"examples/electrostatic/beam_cyclotron/#Beam-cyclotron-instability","page":"Beam-cyclotron instability","title":"Beam-cyclotron instability","text":"","category":"section"},{"location":"examples/electrostatic/beam_cyclotron/","page":"Beam-cyclotron instability","title":"Beam-cyclotron instability","text":"(Image: Source code) (Image: notebook)","category":"page"},{"location":"examples/electrostatic/beam_cyclotron/","page":"Beam-cyclotron instability","title":"Beam-cyclotron instability","text":"Coming soon...","category":"page"},{"location":"examples/electrostatic/beam_cyclotron/","page":"Beam-cyclotron instability","title":"Beam-cyclotron instability","text":"","category":"page"},{"location":"examples/electrostatic/beam_cyclotron/","page":"Beam-cyclotron instability","title":"Beam-cyclotron instability","text":"This page was generated using DemoCards.jl and Literate.jl.","category":"page"},{"location":"examples/electrostatic/two_stream/","page":"Two-stream instability","title":"Two-stream instability","text":"EditURL = \"/home/runner/work/ParticleInCell.jl/ParticleInCell.jl/examples/electrostatic/two_stream.jl\"","category":"page"},{"location":"examples/electrostatic/two_stream/#Two-stream-instability","page":"Two-stream instability","title":"Two-stream instability","text":"","category":"section"},{"location":"examples/electrostatic/two_stream/","page":"Two-stream instability","title":"Two-stream instability","text":"(Image: Source code) (Image: notebook)","category":"page"},{"location":"examples/electrostatic/two_stream/","page":"Two-stream instability","title":"Two-stream instability","text":"Coming soon...","category":"page"},{"location":"examples/electrostatic/two_stream/","page":"Two-stream instability","title":"Two-stream instability","text":"","category":"page"},{"location":"examples/electrostatic/two_stream/","page":"Two-stream instability","title":"Two-stream instability","text":"This page was generated using DemoCards.jl and Literate.jl.","category":"page"},{"location":"theory/intro_to_pic/#Introduction-to-particle-in-cell-simulation","page":"PIC simulation","title":"Introduction to particle-in-cell simulation","text":"","category":"section"},{"location":"theory/intro_to_pic/","page":"PIC simulation","title":"PIC simulation","text":"Classical plasma physics considers the motion of charged particles. The dynamics of these particles will be effected by the presence of externally imposed electric and magnetic fields, which is relatively easy to model since the motion of each particle is independent. However, the particles will themselves source an electric field–-and, if they are moving fast enough, a magnetic field–-due to their charge. The dynamics of other particles will be influenced by these 'self-consistent' fields, which corresponding source fields of their own. Thus, the dynamics of all of the particles is coupled, and so their equations of motion must be solve together. For a typical plasma, the resulting differential equations cannot be solved analytically, and the number of degrees of freedom means that direct computational integration of the equations is also impossible. Thus, plasma physics relies on various simplifications and assumptions to create reduced models that can be solved–- either exactly, or approximately.","category":"page"},{"location":"theory/intro_to_pic/#Approximations-to-the-true-particle-distribution-function","page":"PIC simulation","title":"Approximations to the true particle distribution function","text":"","category":"section"},{"location":"theory/intro_to_pic/","page":"PIC simulation","title":"PIC simulation","text":"It is convenient to represent the locations and momenta of the particles using the distribution function","category":"page"},{"location":"theory/intro_to_pic/","page":"PIC simulation","title":"PIC simulation","text":"f(xpt) = sum_i=1^N delta(x - x_i(t)) delta(p - p_i(t))","category":"page"},{"location":"theory/intro_to_pic/","page":"PIC simulation","title":"PIC simulation","text":"where N is the total number of particles. The evolution of this distribution function can be written using the Maxwell-Boltzmann system","category":"page"},{"location":"theory/intro_to_pic/","page":"PIC simulation","title":"PIC simulation","text":"beginaligned\nTODO\nTODO-Maxwell\nendaligned","category":"page"},{"location":"theory/intro_to_pic/","page":"PIC simulation","title":"PIC simulation","text":"Unfortunately, for almost all plasmas of interest, N is enormous; thus direct simulation of the equations of motion is computationally intractable. The solution is to recognize that if N is large, then any physically small phase-space region will contain many particles, and so we c","category":"page"},{"location":"theory/intro_to_pic/#Todo","page":"PIC simulation","title":"Todo","text":"","category":"section"},{"location":"theory/intro_to_pic/","page":"PIC simulation","title":"PIC simulation","text":"course-graining\ncollisionless plasmas\nbriefly mention thermalization, two-fluid, and MHD","category":"page"},{"location":"theory/intro_to_pic/#Discretized-solutions-of-the-Boltzmann-Poisson-equations","page":"PIC simulation","title":"Discretized solutions of the Boltzmann-Poisson equations","text":"","category":"section"},{"location":"theory/intro_to_pic/","page":"PIC simulation","title":"PIC simulation","text":"We have seen that a kinetic plasma can be approximated using the Vlasov-Boltzmann equation, along with an appropriate equation of motion for the electromagnetic fields. However, the resulting equation is still not easy to analyse for an arbitrary f. We therefore turn to computational simulation of the equations of motion to understand the plasma dynamics.","category":"page"},{"location":"theory/intro_to_pic/","page":"PIC simulation","title":"PIC simulation","text":"In order to simulate the plasma, we must first discretize the equations so they can be represented on a computer. One option, called a Vlasov code, represents the distribution function as","category":"page"},{"location":"theory/intro_to_pic/#The-standard-PIC-cycle","page":"PIC simulation","title":"The standard PIC cycle","text":"","category":"section"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"EditURL = \"/home/runner/work/ParticleInCell.jl/ParticleInCell.jl/examples/index.md\"","category":"page"},{"location":"examples/#Example-Gallery","page":"Example Gallery","title":"Example Gallery","text":"","category":"section"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"These examples are in progress...","category":"page"},{"location":"examples/#Tutorial","page":"Example Gallery","title":"Tutorial","text":"","category":"section"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n\n \n \n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"(Image: list-card-cover-image)","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"Tutorial: Langmuir oscillations","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"In this tutorial, you will use ParticleInCell to model one of the simplest phenomena in plasma physics: an electrostatic (or Langmuir) oscillation. This tutorial is part of a series of examples that uses ParticleInCell to demonstrate the basic plasma physics concepts that are covered in Birdsall and Langdon's classic PIC textbook.","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n
\n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
","category":"page"},{"location":"examples/#Electrostatic","page":"Example Gallery","title":"Electrostatic","text":"","category":"section"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n\n \n \n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"(Image: list-card-cover-image)","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"Two-stream instability","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"Coming soon...","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n
\n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n\n \n \n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"(Image: list-card-cover-image)","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"Beam-plasma instability","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"Coming soon...","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n
\n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n\n \n \n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"(Image: list-card-cover-image)","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"Beam-cyclotron instability","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"Coming soon...","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n
\n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n\n \n \n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"(Image: list-card-cover-image)","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"Landau damping","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"Coming soon...","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n
\n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
","category":"page"},{"location":"examples/#Electromagnetic","page":"Example Gallery","title":"Electromagnetic","text":"","category":"section"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n\n \n \n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"(Image: list-card-cover-image)","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"Beam heating","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"Coming soon...","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
\n
\n
","category":"page"},{"location":"examples/","page":"Example Gallery","title":"Example Gallery","text":"
","category":"page"},{"location":"examples/electromagnetic/beam_heating/","page":"Beam heating","title":"Beam heating","text":"EditURL = \"/home/runner/work/ParticleInCell.jl/ParticleInCell.jl/examples/electromagnetic/beam_heating.jl\"","category":"page"},{"location":"examples/electromagnetic/beam_heating/#Beam-heating","page":"Beam heating","title":"Beam heating","text":"","category":"section"},{"location":"examples/electromagnetic/beam_heating/","page":"Beam heating","title":"Beam heating","text":"(Image: Source code) (Image: notebook)","category":"page"},{"location":"examples/electromagnetic/beam_heating/","page":"Beam heating","title":"Beam heating","text":"Coming soon...","category":"page"},{"location":"examples/electromagnetic/beam_heating/","page":"Beam heating","title":"Beam heating","text":"","category":"page"},{"location":"examples/electromagnetic/beam_heating/","page":"Beam heating","title":"Beam heating","text":"This page was generated using DemoCards.jl and Literate.jl.","category":"page"},{"location":"theory/#Plasma-Simulation-Theory","page":"Introduction","title":"Plasma Simulation Theory","text":"","category":"section"},{"location":"theory/","page":"Introduction","title":"Introduction","text":"The fundamental physics governing the dynamics of a plasma have been well understood for over a century, and yet plasma physics remains an active area of research. This is because the dynamics of a plasma are highly nonlinear, and it is therefore difficult to make analytic statements about how a given plasma will behave over long time-spans. Instead, theoretical plasma physicists often rely on simulation to understand plasma dynamics, and to make general statements about the behaviors of particular classes of plasmas.","category":"page"},{"location":"theory/","page":"Introduction","title":"Introduction","text":"Unfortunately, the simulation of plasmas is itself a quite challenging task because plasmas are composed of many, many, charged particles. As a result, there are far too many degrees of freedom to exactly solve the full equations of motion for a given plasma. Instead, physicists rely on approximations to derive physically relevant models for the systems in question.","category":"page"},{"location":"theory/","page":"Introduction","title":"Introduction","text":"The most \"realistic\" class of models are called kinetic models. In these models, the individual charged particles of each species are averaged to create distribution functions that depends on both position and velocity. The distribution functions give the likelihood of finding a charged particle of a particular species at any point in phase space. These distribution functions, along with a method for calculating the self-consistent interaction between the particles, yields a set of approximate equations of motion describing the plasma dynamics.","category":"page"},{"location":"theory/","page":"Introduction","title":"Introduction","text":"If the particle species are nonrelativistic, then the particles are not influenced by self-consistent magnetic fields, and so the interactions can be modeled using electrostatics. However, once the particles (typically the lightweight electrons) become relativistic, the sourced magnetic field must be taken into account, and so full electromagnetic interactions are required.","category":"page"},{"location":"theory/","page":"Introduction","title":"Introduction","text":"Over long periods of time, the plasma will begin to thermalize–-the distribution of particle velocities will become closer and closer to Maxwellian. This fact can be used to drastically improve the size and speed of simulations by using the two-fluid and magnetohydrodynamic (MHD) approximations of plasma dynamics. As this package does not implement algorithms for these simulation methods, we will not discuss the details of these methods further.","category":"page"},{"location":"theory/","page":"Introduction","title":"Introduction","text":"In the following sections, we describe the details of particle-in-cell algorithms, a specific type of kinetic simulation algorithm. For a more in depth introduction of PIC simulation theory, see the books by Birdsall and Langdon (2004) and Hockney and Eastwood (1989).","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"EditURL = \"/home/runner/work/ParticleInCell.jl/ParticleInCell.jl/examples/tutorial/langmuir_oscillation.jl\"","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/#tutorial_langmuir","page":"Tutorial","title":"Tutorial: Langmuir oscillations","text":"","category":"section"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"(Image: Source code) (Image: notebook)","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"In this tutorial, you will use ParticleInCell to model one of the simplest phenomena in plasma physics: an electrostatic (or Langmuir) oscillation. This tutorial is part of a series of examples that uses ParticleInCell to demonstrate the basic plasma physics concepts that are covered in Birdsall and Langdon's classic PIC textbook.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"A Langmuir oscillation occurs when a slab of charge in a uniform plasma is displaced. The resulting charge density gradient creates a restoring force that causes the displaced slab of charge to return to its original position. But–-just as in a classical pendulum oscillation–-the momentum of the charge carries it past its equilibrium point, creating an opposite charge gradient, and a restoring force in the opposite direction. As a result, the slab of charge oscillates around its equilibrium forever (at least in this idealized model that ignores possible damping mechanisms). For a plasma composed of a single mobile species s with mass m_s and charge q_s, the frequency of this oscillation is given by","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"omega_ps = sqrtfracn_s q_s^2epsilon_0 m_s","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"where n_s is the number density of the plasma and epsilon_0 is the permitivity of free space. Notice that the plasma frequency has a m_s^-12 dependence, and thus the lightest species (typically electrons) will dominate the dynamics of a plasma oscillation. For this reason, we will only model the dynamics of the electrons in our simulation.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/#Simulating-a-cold-electron-plasma","page":"Tutorial","title":"Simulating a cold electron plasma","text":"","category":"section"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"We begin by loading the ParticleInCell package. Additionally, we load CairoMakie which is a backend for Makie that can generate beautiful, publication-quality graphics.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"using ParticleInCell\nusing CairoMakie","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"We begin by creating some electrons to move in the simulation. For even a tiny simulation volume, there are far too many physical electrons to simulate each one individually. Instead, PIC algorithms group physical particles into 'macroparticles'. The distribution of macroparticles in phase space serves as an approximation for the phase space distribution of physical particles. We arbitrarily choose a simulation domain of length one, and a nominal electrons number density of 10^14. Then, for a given number of macroparticles, we can calculate the number of physical electrons represented by each.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"sim_length = 1.0\nnumber_density = 1e14\nnum_macroparticles = 320\nparticles_per_macro = number_density * sim_length / num_macroparticles","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"3.125e11","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"We distribute the macroparticles evenly across the simulation domain.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"positions = collect(0:num_macroparticles-1) ./ num_macroparticles;","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"In order to seed a Langmuir oscillation, we give the electrons a sinusoidal velocity perturbation. This corresponds to the moment in a Langmuir oscillation when the slab of charge has reached equilibrium, but is being carried past by its momentum. This perturbation is defined by a wavenumber k and an amplitude.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"k = 1 * 2pi / sim_length\namplitude = 1e3\nelec_mass = 9e-31\nmomentums = (particles_per_macro * elec_mass * amplitude) .* sin.(positions .* k);","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"We can visualize the initial condition of the electron macroparticle by plotting the initial phase space.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"scatter(\n positions,\n momentums;\n axis = (;\n title = \"Electron phase space\",\n xlabel = \"Position (m)\",\n ylabel = \"Momentum (kg m / s)\",\n ),\n)","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"(Image: )","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"Finally, we create a VariableWeightSpecies which holds the all of the macroparticles. Additionally, we must pass the value of particles_per_macro, which is used to calculate the charge and mass of the macroparticles.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"electrons = ParticleInCell.electrons(positions, momentums, particles_per_macro);","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"Now we address the 'cell' piece of particle-in-cell by creating a grid. Because Langmuir oscillations are a one-dimensional phenomena, we will choose to perform a 1D simulation.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"The choice of grid resolution is determined by the scale of the smallest relevant dynamics begin simulated. For a Langmuir oscillation, the scale of the dynamics is set by k, and so the simulation could likely accomplished with as few as 4 or 8 cells. However, this is not a computationally demanding simulation, and so we arbitrarily choose to use 32 equally spaced (i.e. uniform) grid points. Additionally, we make the simulation domain periodic.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"num_cells = 32\ndx = sim_length / num_cells\nperiodic = true\ngrid = UniformCartesianGrid((0.0,), (sim_length,), (num_cells,), (periodic,));","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"Next, we set up the required fields for an electrostatic PIC simulation. In a basic PIC cycle, we first compute the charge density, rho, on the grid points. We then compute the corresponding electric potential, phi. The electric field is conventionally determined in a two step process. First the potential, which is located at the nodes of the grid cells, is finite differenced to the edges of the cells, producing an edge electric field, Eedge. The edge electric fields are then averaged to get the electric fields located at the nodes, Enode.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"When creating the fields, we must specify the underlying grid on which the field is based, the location of the values of the field (i.e. are the field values located at the nodes of each cell? The edge?), the dimension of the field, and the number of guard cells surrounding the field.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"field_dimension = 1\nlower_guard_cells = 1\nrho = Field(grid, ParticleInCell.node, field_dimension, lower_guard_cells)\nphi = Field(grid, ParticleInCell.node, field_dimension, lower_guard_cells)\nEedge = Field(grid, ParticleInCell.edge, field_dimension, lower_guard_cells)\nEnode = Field(grid, ParticleInCell.node, field_dimension, lower_guard_cells);","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"At this point, we must choose a timestep for the simulation. We would like to use a large timestep so that more of the systems dynamics can be observed with the same number of steps. However, we must resolve the fastest timescale of the dynamics that we are trying to simulate. In this case, we must resolve the plasma frequency. Additionally, we must choose a timestep that is short enough that particles do not cross more than one cell per timestep to prevent numerical instabilities from arising. For the oscillation amplitude that we have chosen, the particles do not move fast enough for the CFL condition to matter, and so we will choose our timestep based on the expected plasma frequency.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"epsilon_0 = 8.8e-12\nelec_charge = 1.6e-19\nelec_mass = 9e-31\nexpected_plasma_freq = sqrt(number_density * elec_charge^2 / elec_mass / epsilon_0)\nexpected_plasma_period = 2pi / expected_plasma_freq","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"1.1051531770007306e-8","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"Once again, this is not a computationally demanding simulation, and so we will choose a relatively small timestep for improve numerical accuracy. You can play with increasing the timestep, and see when the simulation results begin to deteriorate.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"dt = 5e-11","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"5.0e-11","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"In the final step of the setup, we create all of the simulation steps required to do the electrostatic simulation. In this tutorial, we will not discuss the details of PIC simulation, but you can find more information about the PIC simulation cycle elsewhere in this documentation.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"charge_interp = BSplineChargeInterpolation(electrons, rho, 1)\ncomm_rho = CommunicateGuardCells(rho, true)\nfield_solve = PoissonSolveFFT(rho, phi)\ncomm_phi = CommunicateGuardCells(phi)\nfinite_diff = FiniteDifferenceToEdges(phi, Eedge)\ncomm_Eedge = CommunicateGuardCells(Eedge)\nelec_ave = AverageEdgesToNodes(Eedge, Enode)\ncomm_Enode = CommunicateGuardCells(Enode)\npush = ElectrostaticParticlePush(electrons, Enode, dt)\ncomm_electrons = CommunicateSpecies(electrons, grid);","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"Now we are ready to run the simulation. We will simulate the plasma for 1000 timesteps, and at each step, we will calculate the electric field energy,","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"U_E = int E(x)^2 mathrmdx","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"This field energy will oscillate as the electrons move in and out of equilibrium, and so we can use it to observe the Langmuir oscillation.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"n_steps = 1000\n\nelectric_field_energy = Vector{Float64}(undef, n_steps)\n\nfor n = 1:n_steps\n # Calculate the electric field energy\n electric_field_energy[n] = 0\n for I in eachindex(Enode)\n electric_field_energy[n] += (dx * epsilon_0 / 2) * (Enode.values[I])^2\n end\n\n # TODO\n rho.values .= 0\n\n step!(charge_interp)\n step!(comm_rho)\n step!(field_solve)\n step!(comm_phi)\n step!(finite_diff)\n step!(comm_Eedge)\n step!(elec_ave)\n step!(comm_Enode)\n step!(push)\n step!(comm_electrons)\nend","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"We can now visualize the electric field energy to see the plasma oscillation.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"times = collect(range(1, n_steps)) .* dt\nlines(\n times,\n electric_field_energy;\n axis = (; title = \"Electric Field Energy\", xlabel = \"Time (s)\", ylabel = \"Energy\"),\n)","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"(Image: )","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"Notice that the electric field energy is slowly growing over time, which is unphysical. We will discuss where this numerical instability comes from– and how it can be avoided–later. But for now, we can still use the electric-field-energy time series to calculate the plasma frequency. First, let's plot the Fourier transform of the electric field energy.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"using FFTW\n\nfreqs = fftfreq(n_steps, 1 / dt) .* 2pi\nfreq_amps = abs.(fft(electric_field_energy))\n\nlines(\n freqs,\n freq_amps;\n axis = (;\n title = \"Electric Field Energy Frequency Spectrum\",\n xlabel = \"Frequency (1/s)\",\n ylabel = \"Amplitude\",\n ),\n)","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"(Image: )","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"It is hard to see what is happening at the low frequencies, so let's zoom in on the positive low frequencies.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"cutoff_index = round(Int, n_steps * 0.05)\nlines(\n freqs[1:cutoff_index],\n freq_amps[1:cutoff_index];\n axis = (;\n title = \"Electric Field Energy Frequency Spectrum\",\n xlabel = \"Frequency (1/s)\",\n ylabel = \"Amplitude\",\n ),\n)","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"(Image: )","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"Next, we find the maximum frequency. We don't care about the spike at zero frequency (that is just a consequence of the electric field energy being a strictly positive quantity) so we first set its amplitude to zero, and then find the largest remaining amplitude, and it's corresponding frequency.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"freq_amps[1] = 0\nmax_index = findmax(freq_amps)[2]\nmax_freq = freqs[max_index]\n\n# Divide by 2 because the electric field energy goes through a maximum twice\n# per plasma oscillation, and take the absolute value because we don't care\n# about the phase of the oscillation.\nplasma_freq = abs(max_freq / 2)","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"5.654866776461627e8","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"Finally, we can compare this to the theoretically expected result:","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"relative_error = (plasma_freq - expected_plasma_freq) / expected_plasma_freq","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"-0.005362140699342522","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"Less than 1% error. Not bad!","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/#Adding-temperature-to-the-plasma","page":"Tutorial","title":"Adding temperature to the plasma","text":"","category":"section"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"In the previous simulation, the electric field energy grew unphysically throughout the simulation. This was a result of the \"grid-heating instability\", which occurs when the grid does not resolve the plasma Debye length, which for an electron plasma is given by","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"lambda_De = sqrtfracepsilon_0 k_B Tn_e q_e^2","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"where k_B is the Boltzmann constant and T_e is the electron temperature.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"When the Debye length of a plasma is underresolved, the plasma will unphysically heat, causing the Debye length to grow until it is resolved by the grid. Many strategies have been developed to mitigate this instability, but in this tutorial, we will simply give our plasma sufficient thermal energy to begin so that the simulation will be stable against the grid-heating instability.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"Let's calculate the required electron temperature in the previous simulation so that the 32 cell grid will resolve the Debye length. We set lambda_De = Delta x, and solve for T to find","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"boltzmann_constant = 1.381e-23\ndx^2 * number_density * elec_charge^2 / epsilon_0 / boltzmann_constant","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"2.0571390955170825e7","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"Alternatively, we can express this temperature in terms of electron volts as","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"dx^2 * number_density * elec_charge / epsilon_0","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"1775.5681818181818","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"We will define a function that takes an electron density, electron temperature, and oscillation wavenumber, and returns a measured plasma frequency.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"function measure_plasma_frequency(number_density, temperature, wavenumber)\n sim_length = 1.0\n\n num_cells = 32\n dx = sim_length / num_cells\n\n num_macroparticles = 320\n particles_per_macro = number_density * sim_length / num_macroparticles\n\n perturb_amplitude = 1e3\n elec_mass = 9e-31\n boltzmann_constant = 1.381e-23\n thermal_velocity = sqrt(3 * boltzmann_constant * temperature / elec_mass)\n\n positions = collect(0:num_macroparticles-1) ./ num_macroparticles\n momentums =\n (particles_per_macro * elec_mass) .*\n (perturb_amplitude .* sin.(positions .* wavenumber) .+ thermal_velocity .* randn.())\n\n electrons = ParticleInCell.electrons(positions, momentums, particles_per_macro)\n\n grid = UniformCartesianGrid((0.0,), (sim_length,), (num_cells,), (true,))\n\n field_dimension = 1\n lower_guard_cells = 1\n rho = Field(grid, ParticleInCell.node, field_dimension, lower_guard_cells)\n phi = Field(grid, ParticleInCell.node, field_dimension, lower_guard_cells)\n Eedge = Field(grid, ParticleInCell.edge, field_dimension, lower_guard_cells)\n Enode = Field(grid, ParticleInCell.node, field_dimension, lower_guard_cells)\n\n dt = 5e-11\n charge_interp = BSplineChargeInterpolation(electrons, rho, 1)\n comm_rho = CommunicateGuardCells(rho, true)\n field_solve = PoissonSolveFFT(rho, phi)\n comm_phi = CommunicateGuardCells(phi)\n finite_diff = FiniteDifferenceToEdges(phi, Eedge)\n comm_Eedge = CommunicateGuardCells(Eedge)\n elec_ave = AverageEdgesToNodes(Eedge, Enode)\n comm_Enode = CommunicateGuardCells(Enode)\n push = ElectrostaticParticlePush(electrons, Enode, dt)\n comm_electrons = CommunicateSpecies(electrons, grid)\n\n n_steps = 1000\n electric_field_energy = Vector{Float64}(undef, n_steps)\n\n epsilon_0 = 8.8e-12\n for n = 1:n_steps\n # Calculate the electric field energy\n electric_field_energy[n] = 0\n for I in eachindex(Enode)\n electric_field_energy[n] += (dx * epsilon_0 / 2) * (Enode.values[I])^2\n end\n\n # TODO\n rho.values .= 0\n\n step!(charge_interp)\n step!(comm_rho)\n step!(field_solve)\n step!(comm_phi)\n step!(finite_diff)\n step!(comm_Eedge)\n step!(elec_ave)\n step!(comm_Enode)\n step!(push)\n step!(comm_electrons)\n end\n\n freqs = fftfreq(n_steps, 1 / dt) .* 2pi\n freq_amps = abs.(fft(electric_field_energy))\n\n freq_amps[1] = 0\n max_index = findmax(freq_amps)[2]\n max_freq = freqs[max_index]\n plasma_freq = abs(max_freq / 2)\n\n return plasma_freq\nend","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"measure_plasma_frequency (generic function with 1 method)","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"For a warm plasma, the thermal pressure acts as an additional restoring force on the plasma oscillation. It can be show that the modified dispersion relation (frequency response as a function of wavenumber) is given by","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"omega^2 = omega_pe^2 + frac3 k_B T_em_e k^2","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"Notice that when T_e = 0, the frequency is the plasma frequency, regardless of the wavenumber.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"Let's run a few simulations and verify that this relationship holds.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"temperatures = [0, 0.1, 1, 10]\nmeasure_plasma_frequency.(1e14, temperatures, 1 / 2 * pi)","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"4-element Vector{Float64}:\n 5.654866776461627e8\n 5.654866776461627e8\n 4.3982297150257105e8\n 6.2831853071795866e7","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"We can compare this to the expected result.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"function warm_plasma_freq(number_density, temperature, wavenumber)\n epsilon_0 = 8.8e-12\n elec_charge = 1.6e-19\n elec_mass = 9e-31\n boltzmann_constant = 1.381e-23\n square_plasma_freq = number_density * elec_charge^2 / elec_mass / epsilon_0\n correction_factor = boltzmann_constant * temperature * wavenumber^2 / elec_mass\n return sqrt(square_plasma_freq + correction_factor)\nend\nwarm_plasma_freq.(1e14, temperatures, 1 / 2 * pi)","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"4-element Vector{Float64}:\n 5.685352436149611e8\n 5.685352436182908e8\n 5.685352436482581e8\n 5.685352439479299e8","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"Clearly, the restoring force of the pressure is not large enough to make a difference in this case.","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"","category":"page"},{"location":"examples/tutorial/langmuir_oscillation/","page":"Tutorial","title":"Tutorial","text":"This page was generated using DemoCards.jl and Literate.jl.","category":"page"},{"location":"examples/electrostatic/beam_plasma/","page":"Beam-plasma instability","title":"Beam-plasma instability","text":"EditURL = \"/home/runner/work/ParticleInCell.jl/ParticleInCell.jl/examples/electrostatic/beam_plasma.jl\"","category":"page"},{"location":"examples/electrostatic/beam_plasma/#Beam-plasma-instability","page":"Beam-plasma instability","title":"Beam-plasma instability","text":"","category":"section"},{"location":"examples/electrostatic/beam_plasma/","page":"Beam-plasma instability","title":"Beam-plasma instability","text":"(Image: Source code) (Image: notebook)","category":"page"},{"location":"examples/electrostatic/beam_plasma/","page":"Beam-plasma instability","title":"Beam-plasma instability","text":"Coming soon...","category":"page"},{"location":"examples/electrostatic/beam_plasma/","page":"Beam-plasma instability","title":"Beam-plasma instability","text":"","category":"page"},{"location":"examples/electrostatic/beam_plasma/","page":"Beam-plasma instability","title":"Beam-plasma instability","text":"This page was generated using DemoCards.jl and Literate.jl.","category":"page"},{"location":"#ParticleInCell.jl","page":"Introduction","title":"ParticleInCell.jl","text":"","category":"section"},{"location":"","page":"Introduction","title":"Introduction","text":"ParticleInCell","category":"page"},{"location":"#ParticleInCell","page":"Introduction","title":"ParticleInCell","text":"(Image: ParticleInCell.jl Logo)\n\n(Image: Latest Documentation) (Image: CI Status) (Image: Code Coverage Statistics)\n\nParticleInCell.jl is a Julia package for kinetic plasma physics simulation. Specifically, it focuses on the simulation of kinetic (non-thermal) plasmas using particle-in-cell (PIC) algorithms. Currently, this package is in in a pre-1.0.0 state, and thus breaking changes should be expected. However, this also means that I am willing to entertain radical suggestions to improve the functionality of the package. If you are interested in using ParticleInCell for your plasma research, and you find that it does not meet you needs, please reach out on either GitHub, or over email, so that we can discuss how the package can be modified to suite your needs.\n\nGetting Started\n\nParticleInCell is currently not registered in the Julia package registry. Thus, to install this package, you should use Pkg.develop:\n\nusing Pkg\nPkg.develop(url=\"https://github.com/JuliaPlasma/ParticleInCell.jl\")\n\nDocumentation\n\nYou can view the latest documentation here.\n\nGoals\n\nFast: aim to have core time of less than 1 microsecond per particle per step without collisions.\nFlexible: it should be possible to implement essentially any kinetic plasma simulation in ParticleInCell.jl. For common types of simulations, this might mean just piecing together components that are already included. More esoteric problems might require writing custom types that implement the desired algorithms. The advantage of writing this package in Julia is that these custom types will be just as performant as native components that are included in the package.\nScalable: the eventual goal is to enable scaling across an essentially unlimited number of cores using Julia's native multithreading for parallelization on a single node, and MPI.jl for communication across nodes. The goal is to support two different modes of parallelization:\nEach core is responsible for a single rectangular subdomain. The domain assigned to an entire node is also rectangular, which imposes constraints on how the node domain can be subdivided into subdomains for each core. Load balancing is achieved by varying the relative sizes of the domains such that each core has a similar amount of work per step.\nThe simulation domain is subdivided into subdomains called 'patches', and every node is assigned a list of patches that it is responsible for updating. The cores on each node work collaboratively on the list, each choosing one patch to work on, and then selecting another when they are finished. Load balancing is achieved by swapping patches between nodes to balance the workload while also seeking to minimize communication time by keeping the surface area of each node's responsibilities minimized. In order for this scheme to effectively load balance, it must be the case that the total number of patches is larger (ideally much larger) that the total number of cores.\n\n\n\n\n\n","category":"module"},{"location":"#Documentation-sections","page":"Introduction","title":"Documentation sections","text":"","category":"section"},{"location":"","page":"Introduction","title":"Introduction","text":"To get started, look at the Tutorial, which includes step-by-step instructions for running your first simulation using ParticleInCell. After that, you may want to browse the Example Gallery, to see other problems that this package can be used to solve. The Plasma Simulation Theory section contains information about the art of computational plasma physics, and the numerical constraints that require the specialized tools developed here. Finally, the Manual contains detailed information about the entire public interface of the package.","category":"page"}] } diff --git a/dev/theory/index.html b/dev/theory/index.html index 3668f9b..4bb88d5 100644 --- a/dev/theory/index.html +++ b/dev/theory/index.html @@ -1,2 +1,2 @@ -Introduction · ParticleInCell.jl Documentation
GitHub

Plasma Simulation Theory

The fundamental physics governing the dynamics of a plasma have been well understood for over a century, and yet plasma physics remains an active area of research. This is because the dynamics of a plasma are highly nonlinear, and it is therefore difficult to make analytic statements about how a given plasma will behave over long time-spans. Instead, theoretical plasma physicists often rely on simulation to understand plasma dynamics, and to make general statements about the behaviors of particular classes of plasmas.

Unfortunately, the simulation of plasmas is itself a quite challenging task because plasmas are composed of many, many, charged particles. As a result, there are far too many degrees of freedom to exactly solve the full equations of motion for a given plasma. Instead, physicists rely on approximations to derive physically relevant models for the systems in question.

The most "realistic" class of models are called kinetic models. In these models, the individual charged particles of each species are averaged to create distribution functions that depends on both position and velocity. The distribution functions give the likelihood of finding a charged particle of a particular species at any point in phase space. These distribution functions, along with a method for calculating the self-consistent interaction between the particles, yields a set of approximate equations of motion describing the plasma dynamics.

If the particle species are nonrelativistic, then the particles are not influenced by self-consistent magnetic fields, and so the interactions can be modeled using electrostatics. However, once the particles (typically the lightweight electrons) become relativistic, the sourced magnetic field must be taken into account, and so full electromagnetic interactions are required.

Over long periods of time, the plasma will begin to thermalize–-the distribution of particle velocities will become closer and closer to Maxwellian. This fact can be used to drastically improve the size and speed of simulations by using the two-fluid and magnetohydrodynamic (MHD) approximations of plasma dynamics. As this package does not implement algorithms for these simulation methods, we will not discuss the details of these methods further.

In the following sections, we describe the details of particle-in-cell algorithms, a specific type of kinetic simulation algorithm. For a more in depth introduction of PIC simulation theory, see the books by Birdsall and Langdon (2004) and Hockney and Eastwood (1989).

+Introduction · ParticleInCell.jl Documentation
GitHub

Plasma Simulation Theory

The fundamental physics governing the dynamics of a plasma have been well understood for over a century, and yet plasma physics remains an active area of research. This is because the dynamics of a plasma are highly nonlinear, and it is therefore difficult to make analytic statements about how a given plasma will behave over long time-spans. Instead, theoretical plasma physicists often rely on simulation to understand plasma dynamics, and to make general statements about the behaviors of particular classes of plasmas.

Unfortunately, the simulation of plasmas is itself a quite challenging task because plasmas are composed of many, many, charged particles. As a result, there are far too many degrees of freedom to exactly solve the full equations of motion for a given plasma. Instead, physicists rely on approximations to derive physically relevant models for the systems in question.

The most "realistic" class of models are called kinetic models. In these models, the individual charged particles of each species are averaged to create distribution functions that depends on both position and velocity. The distribution functions give the likelihood of finding a charged particle of a particular species at any point in phase space. These distribution functions, along with a method for calculating the self-consistent interaction between the particles, yields a set of approximate equations of motion describing the plasma dynamics.

If the particle species are nonrelativistic, then the particles are not influenced by self-consistent magnetic fields, and so the interactions can be modeled using electrostatics. However, once the particles (typically the lightweight electrons) become relativistic, the sourced magnetic field must be taken into account, and so full electromagnetic interactions are required.

Over long periods of time, the plasma will begin to thermalize–-the distribution of particle velocities will become closer and closer to Maxwellian. This fact can be used to drastically improve the size and speed of simulations by using the two-fluid and magnetohydrodynamic (MHD) approximations of plasma dynamics. As this package does not implement algorithms for these simulation methods, we will not discuss the details of these methods further.

In the following sections, we describe the details of particle-in-cell algorithms, a specific type of kinetic simulation algorithm. For a more in depth introduction of PIC simulation theory, see the books by Birdsall and Langdon (2004) and Hockney and Eastwood (1989).

diff --git a/dev/theory/intro_to_pic/index.html b/dev/theory/intro_to_pic/index.html index 0d22352..30f3e94 100644 --- a/dev/theory/intro_to_pic/index.html +++ b/dev/theory/intro_to_pic/index.html @@ -2,4 +2,4 @@ PIC simulation · ParticleInCell.jl Documentation
GitHub

Introduction to particle-in-cell simulation

Classical plasma physics considers the motion of charged particles. The dynamics of these particles will be effected by the presence of externally imposed electric and magnetic fields, which is relatively easy to model since the motion of each particle is independent. However, the particles will themselves source an electric field–-and, if they are moving fast enough, a magnetic field–-due to their charge. The dynamics of other particles will be influenced by these 'self-consistent' fields, which corresponding source fields of their own. Thus, the dynamics of all of the particles is coupled, and so their equations of motion must be solve together. For a typical plasma, the resulting differential equations cannot be solved analytically, and the number of degrees of freedom means that direct computational integration of the equations is also impossible. Thus, plasma physics relies on various simplifications and assumptions to create reduced models that can be solved–- either exactly, or approximately.

Approximations to the true particle distribution function

It is convenient to represent the locations and momenta of the particles using the distribution function

\[f(x,p,t) = \sum_{i=1}^N \delta(x - x_i(t)) \delta(p - p_i(t)),\]

where $N$ is the total number of particles. The evolution of this distribution function can be written using the Maxwell-Boltzmann system

\[\begin{aligned} TODO TODO-Maxwell -\end{aligned}\]

Unfortunately, for almost all plasmas of interest, $N$ is enormous; thus direct simulation of the equations of motion is computationally intractable. The solution is to recognize that if $N$ is large, then any physically small phase-space region will contain many particles, and so we c

Todo

  • course-graining
  • collisionless plasmas
  • briefly mention thermalization, two-fluid, and MHD

Discretized solutions of the Boltzmann-Poisson equations

We have seen that a kinetic plasma can be approximated using the Vlasov-Boltzmann equation, along with an appropriate equation of motion for the electromagnetic fields. However, the resulting equation is still not easy to analyse for an arbitrary $f$. We therefore turn to computational simulation of the equations of motion to understand the plasma dynamics.

In order to simulate the plasma, we must first discretize the equations so they can be represented on a computer. One option, called a Vlasov code, represents the distribution function as

The standard PIC cycle

+\end{aligned}\]

Unfortunately, for almost all plasmas of interest, $N$ is enormous; thus direct simulation of the equations of motion is computationally intractable. The solution is to recognize that if $N$ is large, then any physically small phase-space region will contain many particles, and so we c

Todo

Discretized solutions of the Boltzmann-Poisson equations

We have seen that a kinetic plasma can be approximated using the Vlasov-Boltzmann equation, along with an appropriate equation of motion for the electromagnetic fields. However, the resulting equation is still not easy to analyse for an arbitrary $f$. We therefore turn to computational simulation of the equations of motion to understand the plasma dynamics.

In order to simulate the plasma, we must first discretize the equations so they can be represented on a computer. One option, called a Vlasov code, represents the distribution function as

The standard PIC cycle