Langmuir Waves

These are examples of Plasma oscillations (Langmuir waves) in a uniform plasma in 1D, 2D, 3D, and RZ.

In each case, a uniform plasma is setup with a sinusoidal perturbation in the electron momentum along each axis. The plasma is followed for a short period of time, long enough so that E fields develop. The resulting fields can be compared to the analytic solutions.

Analyze

We run the following script to analyze correctness:

3D

Script analysis_3d.py

Listing 58 You can copy this file from Examples/Tests/langmuir/analysis_3d.py.

#!/usr/bin/env python3

# Copyright 2019-2022 Jean-Luc Vay, Maxence Thevenet, Remi Lehe, Axel Huebl
#
#
# This file is part of WarpX.
#
# License: BSD-3-Clause-LBNL
#
# This is a script that analyses the simulation results from
# the script `inputs.multi.rt`. This simulates a 3D periodic plasma wave.
# The electric field in the simulation is given (in theory) by:
# $$ E_x = \epsilon \,\frac{m_e c^2 k_x}{q_e}\sin(k_x x)\cos(k_y y)\cos(k_z z)\sin( \omega_p t)$$
# $$ E_y = \epsilon \,\frac{m_e c^2 k_y}{q_e}\cos(k_x x)\sin(k_y y)\cos(k_z z)\sin( \omega_p t)$$
# $$ E_z = \epsilon \,\frac{m_e c^2 k_z}{q_e}\cos(k_x x)\cos(k_y y)\sin(k_z z)\sin( \omega_p t)$$
import os
import re
import sys

import matplotlib.pyplot as plt
import yt
from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable

yt.funcs.mylog.setLevel(50)

import numpy as np
from scipy.constants import c, e, epsilon_0, m_e

# test name
test_name = os.path.split(os.getcwd())[1]

# this will be the name of the plot file
fn = sys.argv[1]

# Parse test name and check if current correction (psatd.current_correction=1) is applied
current_correction = True if re.search("current_correction", test_name) else False

# Parse test name and check if Vay current deposition (algo.current_deposition=vay) is used
vay_deposition = True if re.search("vay_deposition", test_name) else False

# Parse test name and check if div(E)/div(B) cleaning (warpx.do_div<e,b>_cleaning=1) is used
div_cleaning = True if re.search("div_cleaning", test_name) else False

# Parameters (these parameters must match the parameters in `inputs.multi.rt`)
epsilon = 0.01
n = 4.0e24
n_osc_x = 2
n_osc_y = 2
n_osc_z = 2
lo = [-20.0e-6, -20.0e-6, -20.0e-6]
hi = [20.0e-6, 20.0e-6, 20.0e-6]
Ncell = [64, 64, 64]

# Wave vector of the wave
kx = 2.0 * np.pi * n_osc_x / (hi[0] - lo[0])
ky = 2.0 * np.pi * n_osc_y / (hi[1] - lo[1])
kz = 2.0 * np.pi * n_osc_z / (hi[2] - lo[2])
# Plasma frequency
wp = np.sqrt((n * e**2) / (m_e * epsilon_0))

k = {"Ex": kx, "Ey": ky, "Ez": kz}
cos = {"Ex": (0, 1, 1), "Ey": (1, 0, 1), "Ez": (1, 1, 0)}


def get_contribution(is_cos, k, idim):
    du = (hi[idim] - lo[idim]) / Ncell[idim]
    u = lo[idim] + du * (0.5 + np.arange(Ncell[idim]))
    if is_cos[idim] == 1:
        return np.cos(k * u)
    else:
        return np.sin(k * u)


def get_theoretical_field(field, t):
    amplitude = epsilon * (m_e * c**2 * k[field]) / e * np.sin(wp * t)
    cos_flag = cos[field]
    x_contribution = get_contribution(cos_flag, kx, 0)
    y_contribution = get_contribution(cos_flag, ky, 1)
    z_contribution = get_contribution(cos_flag, kz, 2)

    E = (
        amplitude
        * x_contribution[:, np.newaxis, np.newaxis]
        * y_contribution[np.newaxis, :, np.newaxis]
        * z_contribution[np.newaxis, np.newaxis, :]
    )

    return E


# Read the file
ds = yt.load(fn)

# Check that the particle selective output worked:
species = "electrons"
print("ds.field_list", ds.field_list)
for field in ["particle_weight", "particle_momentum_x"]:
    print("assert that this is in ds.field_list", (species, field))
    assert (species, field) in ds.field_list
for field in ["particle_momentum_y", "particle_momentum_z"]:
    print("assert that this is NOT in ds.field_list", (species, field))
    assert (species, field) not in ds.field_list
species = "positrons"
for field in ["particle_momentum_x", "particle_momentum_y"]:
    print("assert that this is NOT in ds.field_list", (species, field))
    assert (species, field) not in ds.field_list

t0 = ds.current_time.to_value()
data = ds.covering_grid(
    level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions
)
edge = np.array(
    [
        (ds.domain_left_edge[2]).item(),
        (ds.domain_right_edge[2]).item(),
        (ds.domain_left_edge[0]).item(),
        (ds.domain_right_edge[0]).item(),
    ]
)

# Check the validity of the fields
error_rel = 0
for field in ["Ex", "Ey", "Ez"]:
    E_sim = data[("mesh", field)].to_ndarray()
    E_th = get_theoretical_field(field, t0)
    max_error = abs(E_sim - E_th).max() / abs(E_th).max()
    print("%s: Max error: %.2e" % (field, max_error))
    error_rel = max(error_rel, max_error)

# Plot the last field from the loop (Ez at iteration 40)
fig, (ax1, ax2) = plt.subplots(1, 2, dpi=100)
# First plot (slice at y=0)
E_plot = E_sim[:, Ncell[1] // 2 + 1, :]
vmin = E_plot.min()
vmax = E_plot.max()
cax1 = make_axes_locatable(ax1).append_axes("right", size="5%", pad="5%")
im1 = ax1.imshow(E_plot, origin="lower", extent=edge, vmin=vmin, vmax=vmax)
cb1 = fig.colorbar(im1, cax=cax1)
ax1.set_xlabel(r"$z$")
ax1.set_ylabel(r"$x$")
ax1.set_title(r"$E_z$ (sim)")
# Second plot (slice at y=0)
E_plot = E_th[:, Ncell[1] // 2 + 1, :]
vmin = E_plot.min()
vmax = E_plot.max()
cax2 = make_axes_locatable(ax2).append_axes("right", size="5%", pad="5%")
im2 = ax2.imshow(E_plot, origin="lower", extent=edge, vmin=vmin, vmax=vmax)
cb2 = fig.colorbar(im2, cax=cax2)
ax2.set_xlabel(r"$z$")
ax2.set_ylabel(r"$x$")
ax2.set_title(r"$E_z$ (theory)")
# Save figure
fig.tight_layout()
fig.savefig("Langmuir_multi_analysis.png", dpi=200)

tolerance_rel = 5e-2

print("error_rel    : " + str(error_rel))
print("tolerance_rel: " + str(tolerance_rel))

assert error_rel < tolerance_rel

# Check relative L-infinity spatial norm of rho/epsilon_0 - div(E)
# with current correction (and periodic single box option) or with Vay current deposition
if current_correction:
    tolerance = 1e-9
elif vay_deposition:
    tolerance = 1e-3
if current_correction or vay_deposition:
    rho = data[("boxlib", "rho")].to_ndarray()
    divE = data[("boxlib", "divE")].to_ndarray()
    error_rel = np.amax(np.abs(divE - rho / epsilon_0)) / np.amax(
        np.abs(rho / epsilon_0)
    )
    print("Check charge conservation:")
    print("error_rel = {}".format(error_rel))
    print("tolerance = {}".format(tolerance))
    assert error_rel < tolerance

if div_cleaning:
    ds_old = yt.load("diags/diag1000038")
    ds_mid = yt.load("diags/diag1000039")
    ds_new = yt.load(fn)  # this is the last plotfile

    ad_old = ds_old.covering_grid(
        level=0, left_edge=ds_old.domain_left_edge, dims=ds_old.domain_dimensions
    )
    ad_mid = ds_mid.covering_grid(
        level=0, left_edge=ds_mid.domain_left_edge, dims=ds_mid.domain_dimensions
    )
    ad_new = ds_new.covering_grid(
        level=0, left_edge=ds_new.domain_left_edge, dims=ds_new.domain_dimensions
    )

    rho = ad_mid["rho"].v.squeeze()
    divE = ad_mid["divE"].v.squeeze()
    F_old = ad_old["F"].v.squeeze()
    F_new = ad_new["F"].v.squeeze()

    # Check max norm of error on dF/dt = div(E) - rho/epsilon_0
    # (the time interval between the old and new data is 2*dt)
    dt = 1.203645751e-15
    x = F_new - F_old
    y = (divE - rho / epsilon_0) * 2 * dt
    error_rel = np.amax(np.abs(x - y)) / np.amax(np.abs(y))
    tolerance = 1e-2
    print("Check div(E) cleaning:")
    print("error_rel = {}".format(error_rel))
    print("tolerance = {}".format(tolerance))
    assert error_rel < tolerance

2D

Script analysis_2d.py

Listing 59 You can copy this file from Examples/Tests/langmuir/analysis_2d.py.

#!/usr/bin/env python3

# Copyright 2019 Jean-Luc Vay, Maxence Thevenet, Remi Lehe
#
#
# This file is part of WarpX.
#
# License: BSD-3-Clause-LBNL
#
# This is a script that analyses the simulation results from
# the script `inputs.multi.rt`. This simulates a 3D periodic plasma wave.
# The electric field in the simulation is given (in theory) by:
# $$ E_x = \epsilon \,\frac{m_e c^2 k_x}{q_e}\sin(k_x x)\cos(k_y y)\cos(k_z z)\sin( \omega_p t)$$
# $$ E_y = \epsilon \,\frac{m_e c^2 k_y}{q_e}\cos(k_x x)\sin(k_y y)\cos(k_z z)\sin( \omega_p t)$$
# $$ E_z = \epsilon \,\frac{m_e c^2 k_z}{q_e}\cos(k_x x)\cos(k_y y)\sin(k_z z)\sin( \omega_p t)$$
import os
import re
import sys

import matplotlib.pyplot as plt
import yt
from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable

yt.funcs.mylog.setLevel(50)

import numpy as np
from scipy.constants import c, e, epsilon_0, m_e

# test name
test_name = os.path.split(os.getcwd())[1]

# this will be the name of the plot file
fn = sys.argv[1]

# Parse test name and check if current correction (psatd.current_correction=1) is applied
current_correction = True if re.search("current_correction", test_name) else False

# Parse test name and check if Vay current deposition (algo.current_deposition=vay) is used
vay_deposition = True if re.search("vay_deposition", test_name) else False

# Parse test name and check if particle_shape = 4 is used
particle_shape_4 = True if re.search("particle_shape_4", test_name) else False

# Parameters (must match the parameters in the inputs)
# FIXME read these parameters from warpx_used_inputs
epsilon = 0.01
n = 4.0e24
n_osc_x = 2
n_osc_z = 2
xmin = -20e-6
xmax = 20.0e-6
Nx = 128
zmin = -20e-6
zmax = 20.0e-6
Nz = 128

# Wave vector of the wave
kx = 2.0 * np.pi * n_osc_x / (xmax - xmin)
kz = 2.0 * np.pi * n_osc_z / (zmax - zmin)
# Plasma frequency
wp = np.sqrt((n * e**2) / (m_e * epsilon_0))

k = {"Ex": kx, "Ez": kz}
cos = {"Ex": (0, 1, 1), "Ez": (1, 1, 0)}


def get_contribution(is_cos, k):
    du = (xmax - xmin) / Nx
    u = xmin + du * (0.5 + np.arange(Nx))
    if is_cos == 1:
        return np.cos(k * u)
    else:
        return np.sin(k * u)


def get_theoretical_field(field, t):
    amplitude = epsilon * (m_e * c**2 * k[field]) / e * np.sin(wp * t)
    cos_flag = cos[field]
    x_contribution = get_contribution(cos_flag[0], kx)
    z_contribution = get_contribution(cos_flag[2], kz)

    E = amplitude * x_contribution[:, np.newaxis] * z_contribution[np.newaxis, :]

    return E


# Read the file
ds = yt.load(fn)
t0 = ds.current_time.to_value()
data = ds.covering_grid(
    level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions
)
edge = np.array(
    [
        (ds.domain_left_edge[1]).item(),
        (ds.domain_right_edge[1]).item(),
        (ds.domain_left_edge[0]).item(),
        (ds.domain_right_edge[0]).item(),
    ]
)

# Check the validity of the fields
error_rel = 0
for field in ["Ex", "Ez"]:
    E_sim = data[("mesh", field)].to_ndarray()[:, :, 0]
    E_th = get_theoretical_field(field, t0)
    max_error = abs(E_sim - E_th).max() / abs(E_th).max()
    print("%s: Max error: %.2e" % (field, max_error))
    error_rel = max(error_rel, max_error)

# Plot the last field from the loop (Ez at iteration 40)
fig, (ax1, ax2) = plt.subplots(1, 2, dpi=100)
# First plot
vmin = E_sim.min()
vmax = E_sim.max()
cax1 = make_axes_locatable(ax1).append_axes("right", size="5%", pad="5%")
im1 = ax1.imshow(E_sim, origin="lower", extent=edge, vmin=vmin, vmax=vmax)
cb1 = fig.colorbar(im1, cax=cax1)
ax1.set_xlabel(r"$z$")
ax1.set_ylabel(r"$x$")
ax1.set_title(r"$E_z$ (sim)")
# Second plot
vmin = E_th.min()
vmax = E_th.max()
cax2 = make_axes_locatable(ax2).append_axes("right", size="5%", pad="5%")
im2 = ax2.imshow(E_th, origin="lower", extent=edge, vmin=vmin, vmax=vmax)
cb2 = fig.colorbar(im2, cax=cax2)
ax2.set_xlabel(r"$z$")
ax2.set_ylabel(r"$x$")
ax2.set_title(r"$E_z$ (theory)")
# Save figure
fig.tight_layout()
fig.savefig("Langmuir_multi_2d_analysis.png", dpi=200)

if particle_shape_4:
    # lower fidelity, due to smoothing
    tolerance_rel = 0.07
else:
    tolerance_rel = 0.05

print("error_rel    : " + str(error_rel))
print("tolerance_rel: " + str(tolerance_rel))

assert error_rel < tolerance_rel

# Check relative L-infinity spatial norm of rho/epsilon_0 - div(E)
# with current correction (and periodic single box option) or with Vay current deposition
if current_correction:
    tolerance = 1e-9
elif vay_deposition:
    tolerance = 1e-3
if current_correction or vay_deposition:
    rho = data[("boxlib", "rho")].to_ndarray()
    divE = data[("boxlib", "divE")].to_ndarray()
    error_rel = np.amax(np.abs(divE - rho / epsilon_0)) / np.amax(
        np.abs(rho / epsilon_0)
    )
    print("Check charge conservation:")
    print("error_rel = {}".format(error_rel))
    print("tolerance = {}".format(tolerance))
    assert error_rel < tolerance

RZ

Script analysis_rz.py

Listing 60 You can copy this file from Examples/Tests/langmuir/analysis_rz.py.

#!/usr/bin/env python3

# Copyright 2019 David Grote, Maxence Thevenet
#
# This file is part of WarpX.
#
# License: BSD-3-Clause-LBNL
#
# This is a script that analyses the simulation results from
# the script `inputs.multi.rz.rt`. This simulates a RZ periodic plasma wave.
# The electric field in the simulation is given (in theory) by:
# $$ E_r = -\partial_r \phi = \epsilon \,\frac{mc^2}{e}\frac{2\,r}{w_0^2} \exp\left(-\frac{r^2}{w_0^2}\right) \sin(k_0 z) \sin(\omega_p t)
# $$ E_z = -\partial_z \phi = - \epsilon \,\frac{mc^2}{e} k_0 \exp\left(-\frac{r^2}{w_0^2}\right) \cos(k_0 z) \sin(\omega_p t)
# Unrelated to the Langmuir waves, we also test the plotfile particle filter function in this
# analysis script.
import os
import re
import sys

import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt
import yt

yt.funcs.mylog.setLevel(50)

import numpy as np
import post_processing_utils
from scipy.constants import c, e, epsilon_0, m_e

# this will be the name of the plot file
fn = sys.argv[1]

# test name
test_name = os.path.split(os.getcwd())[1]

# Parse test name and check if current correction (psatd.current_correction) is applied
current_correction = True if re.search("current_correction", test_name) else False

# Parameters (these parameters must match the parameters in `inputs.multi.rz.rt`)
epsilon = 0.01
n = 2.0e24
w0 = 5.0e-6
n_osc_z = 2
rmin = 0e-6
rmax = 20.0e-6
Nr = 64
zmin = -20e-6
zmax = 20.0e-6
Nz = 128

# Wave vector of the wave
k0 = 2.0 * np.pi * n_osc_z / (zmax - zmin)
# Plasma frequency
wp = np.sqrt((n * e**2) / (m_e * epsilon_0))
kp = wp / c


def Er(z, r, epsilon, k0, w0, wp, t):
    """
    Return the radial electric field as an array
    of the same length as z and r, in the half-plane theta=0
    """
    Er_array = (
        epsilon
        * m_e
        * c**2
        / e
        * 2
        * r
        / w0**2
        * np.exp(-(r**2) / w0**2)
        * np.sin(k0 * z)
        * np.sin(wp * t)
    )
    return Er_array


def Ez(z, r, epsilon, k0, w0, wp, t):
    """
    Return the longitudinal electric field as an array
    of the same length as z and r, in the half-plane theta=0
    """
    Ez_array = (
        -epsilon
        * m_e
        * c**2
        / e
        * k0
        * np.exp(-(r**2) / w0**2)
        * np.cos(k0 * z)
        * np.sin(wp * t)
    )
    return Ez_array


# Read the file
ds = yt.load(fn)
t0 = ds.current_time.to_value()
data = ds.covering_grid(
    level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions
)

# Get cell centered coordinates
dr = (rmax - rmin) / Nr
dz = (zmax - zmin) / Nz
coords = np.indices([Nr, Nz], "d")
rr = rmin + (coords[0] + 0.5) * dr
zz = zmin + (coords[1] + 0.5) * dz

# Check the validity of the fields
overall_max_error = 0
Er_sim = data[("boxlib", "Er")].to_ndarray()[:, :, 0]
Er_th = Er(zz, rr, epsilon, k0, w0, wp, t0)
max_error = abs(Er_sim - Er_th).max() / abs(Er_th).max()
print("Er: Max error: %.2e" % (max_error))
overall_max_error = max(overall_max_error, max_error)

Ez_sim = data[("boxlib", "Ez")].to_ndarray()[:, :, 0]
Ez_th = Ez(zz, rr, epsilon, k0, w0, wp, t0)
max_error = abs(Ez_sim - Ez_th).max() / abs(Ez_th).max()
print("Ez: Max error: %.2e" % (max_error))
overall_max_error = max(overall_max_error, max_error)

# Plot the last field from the loop (Ez at iteration 40)
plt.subplot2grid((1, 2), (0, 0))
plt.imshow(Ez_sim)
plt.colorbar()
plt.title("Ez, last iteration\n(simulation)")
plt.subplot2grid((1, 2), (0, 1))
plt.imshow(Ez_th)
plt.colorbar()
plt.title("Ez, last iteration\n(theory)")
plt.tight_layout()
plt.savefig(test_name + "_analysis.png")

error_rel = overall_max_error

tolerance_rel = 0.12

print("error_rel    : " + str(error_rel))
print("tolerance_rel: " + str(tolerance_rel))

assert error_rel < tolerance_rel

# Check charge conservation (relative L-infinity norm of error) with current correction
if current_correction:
    divE = data[("boxlib", "divE")].to_ndarray()
    rho = data[("boxlib", "rho")].to_ndarray() / epsilon_0
    error_rel = np.amax(np.abs(divE - rho)) / max(np.amax(divE), np.amax(rho))
    tolerance = 1.0e-9
    print("Check charge conservation:")
    print("error_rel = {}".format(error_rel))
    print("tolerance = {}".format(tolerance))
    assert error_rel < tolerance


## In the final past of the test, we verify that the diagnostic particle filter function works as
## expected in RZ geometry. For this, we only use the last simulation timestep.

dim = "rz"
species_name = "electrons"

parser_filter_fn = "diags/diag_parser_filter000080"
parser_filter_expression = "(py-pz < 0) * (r<10e-6) * (z > 0)"
post_processing_utils.check_particle_filter(
    fn, parser_filter_fn, parser_filter_expression, dim, species_name
)

uniform_filter_fn = "diags/diag_uniform_filter000080"
uniform_filter_expression = "ids%3 == 0"
post_processing_utils.check_particle_filter(
    fn, uniform_filter_fn, uniform_filter_expression, dim, species_name
)

random_filter_fn = "diags/diag_random_filter000080"
random_fraction = 0.66
post_processing_utils.check_random_filter(
    fn, random_filter_fn, random_fraction, dim, species_name
)

1D

Script analysis_1d.py

Listing 61 You can copy this file from Examples/Tests/langmuir/analysis_1d.py.

#!/usr/bin/env python3

# Copyright 2019-2022 Jean-Luc Vay, Maxence Thevenet, Remi Lehe, Prabhat Kumar, Axel Huebl
#
#
# This file is part of WarpX.
#
# License: BSD-3-Clause-LBNL
#
# This is a script that analyses the simulation results from
# the script `inputs.multi.rt`. This simulates a 1D periodic plasma wave.
# The electric field in the simulation is given (in theory) by:
# $$ E_z = \epsilon \,\frac{m_e c^2 k_z}{q_e}\sin(k_z z)\sin( \omega_p t)$$
import os
import re
import sys

import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt
import yt

yt.funcs.mylog.setLevel(50)

import numpy as np
from scipy.constants import c, e, epsilon_0, m_e

# test name
test_name = os.path.split(os.getcwd())[1]

# this will be the name of the plot file
fn = sys.argv[1]

# Parse test name and check if current correction (psatd.current_correction=1) is applied
current_correction = True if re.search("current_correction", test_name) else False

# Parse test name and check if Vay current deposition (algo.current_deposition=vay) is used
vay_deposition = True if re.search("vay_deposition", test_name) else False

# Parameters (these parameters must match the parameters in `inputs.multi.rt`)
epsilon = 0.01
n = 4.0e24
n_osc_z = 2
zmin = -20e-6
zmax = 20.0e-6
Nz = 128

# Wave vector of the wave
kz = 2.0 * np.pi * n_osc_z / (zmax - zmin)
# Plasma frequency
wp = np.sqrt((n * e**2) / (m_e * epsilon_0))

k = {"Ez": kz}
cos = {"Ez": (1, 1, 0)}


def get_contribution(is_cos, k):
    du = (zmax - zmin) / Nz
    u = zmin + du * (0.5 + np.arange(Nz))
    if is_cos == 1:
        return np.cos(k * u)
    else:
        return np.sin(k * u)


def get_theoretical_field(field, t):
    amplitude = epsilon * (m_e * c**2 * k[field]) / e * np.sin(wp * t)
    cos_flag = cos[field]
    z_contribution = get_contribution(cos_flag[2], kz)

    E = amplitude * z_contribution

    return E


# Read the file
ds = yt.load(fn)
t0 = ds.current_time.to_value()
data = ds.covering_grid(
    level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions
)
# Check the validity of the fields
error_rel = 0
for field in ["Ez"]:
    E_sim = data[("mesh", field)].to_ndarray()[:, 0, 0]
    E_th = get_theoretical_field(field, t0)
    max_error = abs(E_sim - E_th).max() / abs(E_th).max()
    print("%s: Max error: %.2e" % (field, max_error))
    error_rel = max(error_rel, max_error)

# Plot the last field from the loop (Ez at iteration 80)
plt.subplot2grid((1, 2), (0, 0))
plt.plot(E_sim)
# plt.colorbar()
plt.title("Ez, last iteration\n(simulation)")
plt.subplot2grid((1, 2), (0, 1))
plt.plot(E_th)
# plt.colorbar()
plt.title("Ez, last iteration\n(theory)")
plt.tight_layout()
plt.savefig("langmuir_multi_1d_analysis.png")

tolerance_rel = 0.05

print("error_rel    : " + str(error_rel))
print("tolerance_rel: " + str(tolerance_rel))

assert error_rel < tolerance_rel

# Check relative L-infinity spatial norm of rho/epsilon_0 - div(E) when
# current correction (psatd.do_current_correction=1) is applied or when
# Vay current deposition (algo.current_deposition=vay) is used
if current_correction or vay_deposition:
    rho = data[("boxlib", "rho")].to_ndarray()
    divE = data[("boxlib", "divE")].to_ndarray()
    error_rel = np.amax(np.abs(divE - rho / epsilon_0)) / np.amax(
        np.abs(rho / epsilon_0)
    )
    tolerance = 1.0e-9
    print("Check charge conservation:")
    print("error_rel = {}".format(error_rel))
    print("tolerance = {}".format(tolerance))
    assert error_rel < tolerance

Visualize

Note

This section is TODO.

Langmuir Waves

Run

Analyze

Visualize