Cropping of particles at boundaries for deposition for charge conservation

Docs/source/usage/parameters.rst

@@ -497,6 +497,7 @@ Domain Boundary Conditions
     * ``damped``: This is the recommended option in the moving direction when using the spectral solver with moving window (currently only supported along z). This boundary condition applies a damping factor to the electric and magnetic fields in the outer half of the guard cells, using a sine squared profile. As the spectral solver is by nature periodic, the damping prevents fields from wrapping around to the other end of the domain when the periodicity is not desired. This boundary condition is only valid when using the spectral solver.
 
     * ``pec``: This option can be used to set a Perfect Electric Conductor at the simulation boundary. Please see the :ref:`PEC theory section <theory-bc-pec>` for more details. Note that PEC boundary is invalid at `r=0` for RZ, RCYLINDER, and RSPHERE. Please use ``none`` option. This boundary condition does not work with the spectral solver.
+      There is the additional input parameter ``particles.crop_on_PEC_boundary`` which sets whethers particle trajectories are cropped when particles cross PEC boundaries, defaulting to false.
 
     * ``pmc``: This option can be used to set a Perfect Magnetic Conductor at the simulation boundary. Please see the :ref:`PEC theory section <theory-bc-pmc>` for more details. This is equivalent to ``Neumann``. This boundary condition does not work with the spectral solver.
 
@@ -506,6 +507,7 @@ Domain Boundary Conditions
       The expressions are given for the low and high boundary on each axis, as listed below. The tangential fields are specified as
       expressions that can depend on the location and time. The tangential fields are in two pairs, the electric fields and the
       magnetic fields. In each pair, if one is specified, the other will be set to zero if not also specified.
+      There is the additional input parameter ``particles.crop_on_PEC_boundary`` which sets whethers particle trajectories are cropped when particles cross pec_insulator boundaries, defaulting to false.
 
       * ``insulator.area_x_lo(y,z)``: For the lower x (or r) boundary, expression specifying the insulator location
 

Examples/Tests/implicit/CMakeLists.txt

@@ -51,6 +51,16 @@ add_warpx_test(
     OFF  # dependency
 )
 
+add_warpx_test(
+    test_2d_theta_implicit_jfnk_vandb_cropping  # name
+    2  # dims
+    2  # nprocs
+    inputs_test_2d_theta_implicit_jfnk_vandb_cropping  # inputs
+    "analysis_vandb_jfnk_2d_cropping.py diags/diag1000010"  # analysis
+    "analysis_default_regression.py --path diags/diag1000010"  # checksum
+    OFF  # dependency
+)
+
 add_warpx_test(
     test_2d_theta_implicit_jfnk_vandb_filtered  # name
     2  # dims

Examples/Tests/implicit/analysis_vandb_jfnk_2d_cropping.py

@@ -0,0 +1,44 @@
+#!/usr/bin/env python3
+
+# Copyright 2024 Justin Angus
+#
+#
+# This file is part of WarpX.
+#
+# License: BSD-3-Clause-LBNL
+#
+# This is a script that analyses the simulation results from the script `inputs_vandb_2d`.
+# This simulates a 2D periodic plasma using the implicit solver
+# with the Villasenor deposition using shape factor 2.
+import sys
+
+import numpy as np
+import yt
+from scipy.constants import e, epsilon_0
+
+# check for machine precision conservation of charge density
+n0 = 1.0e12
+
+pltdir = sys.argv[1]
+ds = yt.load(pltdir)
+data = ds.covering_grid(
+    level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions
+)
+
+divE = data["boxlib", "divE"].value
+rho = data["boxlib", "rho"].value
+
+# compute local error in Gauss's law
+drho = (rho - epsilon_0 * divE) / e / n0
+
+# compute RMS on in error on the grid
+nX = drho.shape[0]
+nZ = drho.shape[1]
+drho_max = np.abs(drho).max()
+
+tolerance_max_charge = 1.0e-13
+
+print(f"max error in charge conservation: {drho_max}")
+print(f"tolerance: {tolerance_max_charge}")
+
+assert drho_max < tolerance_max_charge

Examples/Tests/implicit/inputs_test_2d_theta_implicit_jfnk_vandb_cropping

@@ -0,0 +1,100 @@
+#################################
+########## CONSTANTS ############
+#################################
+
+my_constants.Nppc = 2 # particles per cell along each axis
+
+#################################
+####### GENERAL PARAMETERS ######
+#################################
+max_step = 10
+amr.n_cell =  16 16
+amr.max_level = 0
+geometry.dims = 2
+geometry.prob_lo = 0.0 0.0
+geometry.prob_hi = 1.0e-2 1.0e-2
+
+#################################
+####### Boundary condition ######
+#################################
+boundary.field_lo = pec pec
+boundary.field_hi = pec pec
+boundary.particle_lo = absorbing absorbing
+boundary.particle_hi = absorbing absorbing
+particles.crop_on_PEC_boundary = 1
+
+#################################
+############ NUMERICS ###########
+#################################
+warpx.serialize_initial_conditions = 1
+warpx.verbose = 1
+warpx.const_dt = 1.0e-12
+warpx.use_filter = 0
+
+algo.maxwell_solver = Yee
+algo.evolve_scheme = "theta_implicit_em"
+
+implicit_evolve.theta = 0.5
+implicit_evolve.max_particle_iterations = 21
+implicit_evolve.particle_tolerance = 1.0e-12
+implicit_evolve.particle_suborbits = 1
+
+implicit_evolve.nonlinear_solver = "newton"
+newton.verbose = true
+newton.max_iterations = 100
+newton.relative_tolerance = 1.0e-12
+newton.absolute_tolerance = 0.0
+newton.require_convergence = false
+
+gmres.verbose_int = 2
+gmres.max_iterations = 1000
+gmres.relative_tolerance = 1.0e-4
+gmres.absolute_tolerance = 0.0
+
+algo.particle_pusher = "boris"
+
+algo.particle_shape = 4
+algo.current_deposition = "villasenor"
+
+#################################
+############ PLASMA #############
+#################################
+particles.species_names = electrons positrons
+
+electrons.charge = -q_e
+electrons.mass = m_e
+electrons.injection_style = "NUniformPerCell"
+electrons.num_particles_per_cell_each_dim = Nppc Nppc
+electrons.profile = constant
+electrons.density = 1.0e12
+electrons.xmin = 0.01e-2
+electrons.xmax = 0.99e-2
+electrons.zmin = 0.01e-2
+electrons.zmax = 0.99e-2
+electrons.momentum_distribution_type = "constant"
+electrons.ux = -0.2
+electrons.uy = 0.0
+electrons.uz = -0.2
+
+positrons.charge = q_e
+positrons.mass = m_e
+positrons.injection_style = "NUniformPerCell"
+positrons.num_particles_per_cell_each_dim = Nppc Nppc
+positrons.profile = constant
+positrons.density = 1.0e12
+positrons.xmin = 0.01e-2
+positrons.xmax = 0.99e-2
+positrons.zmin = 0.01e-2
+positrons.zmax = 0.99e-2
+positrons.momentum_distribution_type = "constant"
+positrons.ux = 0.2
+positrons.uy = 0.0
+positrons.uz = 0.2
+
+### Diagnostics
+diagnostics.diags_names = diag1
+diag1.intervals = 1
+diag1.diag_type = Full
+diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho divE
+diag1.electrons.variables = x z w ux uy uz
+diag1.positrons.variables = x z w ux uy uz

Regression/Checksum/benchmarks_json/test_2d_theta_implicit_jfnk_vandb_cropping.json

@@ -0,0 +1,31 @@
+{
+  "lev=0": {
+    "Bx": 0.0,
+    "By": 8.29634805330963e-06,
+    "Bz": 0.0,
+    "Ex": 3564.3506590945044,
+    "Ey": 0.0,
+    "Ez": 3564.350659094638,
+    "divE": 364198.73206478765,
+    "jx": 3997.2735251148706,
+    "jy": 0.0,
+    "jz": 3997.2735251149425,
+    "rho": 3.2246839770712876e-06
+  },
+  "electrons": {
+    "particle_momentum_x": 4.91566293563828e-20,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 4.915662935638171e-20,
+    "particle_position_x": 4.2619942734894956,
+    "particle_position_y": 4.261994273489464,
+    "particle_weight": 87890625.0
+  },
+  "positrons": {
+    "particle_momentum_x": 4.91566293563828e-20,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 4.915662935638171e-20,
+    "particle_position_x": 4.7380057265105044,
+    "particle_position_y": 4.738005726510536,
+    "particle_weight": 87890625.0
+  }
+}