adding infrastructure for computing the mass matrices to be used by the implicit EM solvers (#5768)

This PR adds the basic infrastructure for computing the mass matrices
needed by the implicit EM solvers to efficiently take time steps larger
than the plasma period.

The mass matrices are deposited to the grid from the particles in the
same fashion as the current is deposited to the grid. The relationship
between dJ and dE is written as
dJx = MassMatrices_xx * dEx + MassMatrices_xy * dEy + MassMatrices_xz *
dEz
dJy = MassMatrices_yx * dEx + MassMatrices_yy * dEy + MassMatrices_yz *
dEz
dJz = MassMatrices_zx * dEx + MassMatrices_zy * dEy + MassMatrices_zz *
dEz

At each grid point, each element in the expressions above is a sum over
all local grid points depending on the particle shape factor. **For this
PR, we are only computing the diagonal elements of the diagonal mass
matrices**:
dJx(i,j,k) = **MassMatrices_xx(i,j,k)** * dEx(i,j,k)
dJy(i,j,k) = **MassMatrices_yy(i,j,k)** * dEy(i,j,k)
dJz(i,j,k) = **MassMatrices_zz(i,j,k)** * dEz(i,j,k)

The diagonal elements are used in the preconditioner for capturing
plasma response.

The mass matrices are described in
Justin Ray Angus, et. al., An implicit particle code with exact energy
and charge conservation for electromagnetic studies of dense plasmas,
Journal of Computational Physics, Volume 491, 2023, 112383,
https://doi.org/10.1016/j.jcp.2023.112383.

<img width="1029" alt="dane_gmres_planar_pinch_longer"
src="https://github.com/user-attachments/assets/edf66b6a-a9ce-4eb1-a519-a183c14f0ce8"
/>


This figure illustrates the effectiveness of including the plasma
response in the preconditioner. These results are from 1D planar pinch
simulations as described in the above reference. The parameters are such
that \omega_{pe}*\Delta_t ranges from 18 at t = 0 to 33 for t > 15 ns
where the density is compressed. The CFL number is 4.2.

The number of GMRES iterations per Newton iteration increases as the
plasma density increases when not using a preconditioner. Using the curl
curl preconditioner without the plasma response actually result in more
iterations. This is not totally unexpected since the stiffness of the
system is dominated by the plasma rather than light waves. When the
diagonal mass matrices are included in the curl curl preconditioner, the
number of iterations decreases and becomes relatively insensitive to the
change in density. This behavior is similar to that for the same
conditions as reported in the above reference.

The number of GMRES iterations can be reduced further by including the
first layer of off-diagonal elements of the mass matrices in the PC.
That will be done in a future PR.

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

Author: JustinRayAngus

Docs/source/usage/parameters.rst

@@ -126,6 +126,10 @@ Overall simulation parameters
     , this sets the relative tolerance for the iterative method used to obtain a self-consistent update of the particles at
     each iteration in the JFNK process.
 
+* ``implicit_evolve.use_mass_matrices`` (`bool`, default: false)
+    When `algo.evolve_scheme` is either `theta_implicit_em`, `strang_implicit_spectral_em`, or `semi_implicit_em` and `implicit_evolve.nonlinear_solver = newton` and a preconditioner is being used
+    , the diagonal components of the diagonal mass matrices are used to capture the plasma response in the preconditioner.
+
 * ``picard.verbose`` (`bool`, default: 1)
     When `implicit_evolve.nonlinear_solver = picard`, this sets the verbosity of the Picard solver. If true, then information
     on the nonlinear error are printed to screen at each nonlinear iteration.

Examples/Tests/implicit/CMakeLists.txt

@@ -1,6 +1,16 @@
 # Add tests (alphabetical order) ##############################################
 #
 
+add_warpx_test(
+    test_1d_planar_pinch_withPC  # name
+    1  # dims
+    2  # nprocs
+    inputs_test_1d_planar_pinch_withPC  # inputs
+    "analysis_1d_pinch.py"  # analysis
+    "analysis_default_regression.py --path diags/diag1000020"  # checksum
+    OFF  # dependency
+)
+
 add_warpx_test(
     test_1d_semi_implicit_picard  # name
     1  # dims

Examples/Tests/implicit/analysis_1d_pinch.py

@@ -0,0 +1,44 @@
+#!/usr/bin/env python3
+
+# This file is part of WarpX.
+#
+# License: BSD-3-Clause-LBNL
+#
+# This is a script that analyses the simulation results from
+# the script `inputs_test_1d_planar_pinch_withPC`.
+# This simulates a 1D planar pinch using the implicit solver with
+# the curl curl PC including the diagonal response from mass matrices.
+
+import numpy as np
+
+newton_solver = np.loadtxt("diags/reduced_files/newton_solver.txt", skiprows=1)
+gmres_iters = newton_solver[:, 5]
+
+field_energy = np.loadtxt("diags/reduced_files/field_energy.txt", skiprows=1)
+particle_energy = np.loadtxt("diags/reduced_files/particle_energy.txt", skiprows=1)
+poynting_flux = np.loadtxt("diags/reduced_files/poynting_flux.txt", skiprows=1)
+
+E_energy = field_energy[:, 3]
+B_energy = field_energy[:, 4]
+Efields = E_energy + B_energy
+ele_energy = particle_energy[:, 3]
+ion_energy = particle_energy[:, 4]
+Eplasma = ele_energy + ion_energy
+
+Eout_lo_x = poynting_flux[:, 4]
+Eout_hi_x = poynting_flux[:, 5]
+
+# check that violation of energy conservation is below tolerance
+dE = Efields + Eplasma + Eout_hi_x + Eout_lo_x
+rel_net_energy = np.abs(dE - dE[0]) / Eplasma
+max_rel_net_energy = rel_net_energy.max()
+rel_net_energy_tol = 1.0e-9
+print(f"max relative delta energy : {max_rel_net_energy}")
+print(f"relative delta energy tolerance : {rel_net_energy_tol}")
+assert max_rel_net_energy < rel_net_energy_tol
+
+# check that the maximum gmres iterations is below tolerance
+max_gmres_iters_tol = 30
+print(f"max gmres iters: {gmres_iters.max()}")
+print(f"gmres iters tolerance: {max_gmres_iters_tol}")
+assert gmres_iters.max() < max_gmres_iters_tol

Examples/Tests/implicit/inputs_test_1d_planar_pinch_withPC

@@ -0,0 +1,155 @@
+#################################
+########## CONSTANTS ############
+#################################
+
+###   characteristic plasma values
+my_constants.n0 = 1.e23    # plasma density [1/m^3]
+my_constants.T0 = 1.       # plasma temperature [eV]
+my_constants.R0 = 1.5e-2   # plasma radius [m]
+my_constants.I0 = 200.e3   # pinch current [Amps]
+
+my_constants.AD = 2.01410  # atomic number deuterium
+my_constants.mD = AD*m_u - m_e   # deuterium mass [kg]
+
+###   derived characteristic pinch values
+my_constants.B0 = mu0*I0/(2*pi*R0)  # B at r=R0 [Tesla]
+my_constants.PB0 = B0^2/(2.0*mu0)   # B pressure at r=R0
+my_constants.rho0 = AD*m_u*n0       # mass density [kg/m^3]
+my_constants.u0 = sqrt(PB0/rho0)    # implosion speed [m/s]
+my_constants.t0 = R0/u0             # implosion time [m/s]
+
+###   simulation parameters
+my_constants.tsim = 220.e-9  # sim time = 1.35*t0 [s]
+my_constants.rise_time = t0/20.0  # current rise time [s]
+my_constants.dt = 1.0e-3*1.0e-9   # time step [s]
+my_constants.Nppc = 1028      # particles per cell
+
+#################################
+####### GENERAL PARAMETERS ######
+#################################
+stop_time = tsim
+max_step = 20
+amr.n_cell =  216
+amr.max_level = 0
+warpx.numprocs = 2
+geometry.dims = 1
+geometry.prob_lo = 0.0
+geometry.prob_hi = 1.54e-2
+
+#################################
+####### Boundary condition ######
+#################################
+boundary.field_lo = pmc
+boundary.field_hi = pec_insulator
+boundary.particle_lo = reflecting
+boundary.particle_hi = absorbing
+
+insulator.area_z_hi(x,y) = 1.0
+insulator.Bx_z_hi(x,y,t) = min(t/rise_time,1)*B0
+
+#################################
+############ NUMERICS ###########
+#################################
+warpx.serialize_initial_conditions = 1
+warpx.verbose = 1
+warpx.const_dt = dt
+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-10
+
+implicit_evolve.nonlinear_solver = "newton"
+newton.verbose = true
+newton.max_iterations = 100
+newton.relative_tolerance = 1.0e-8
+newton.absolute_tolerance = 0.0
+newton.require_convergence = false
+newton.diagnostic_file = "diags/reduced_files/newton_solver.txt"
+newton.diagnostic_interval = 10
+
+gmres.verbose_int = 2
+gmres.max_iterations = 1000
+gmres.relative_tolerance = 1.0e-4
+gmres.absolute_tolerance = 0.0
+
+implicit_evolve.use_mass_matrices = true
+jacobian.pc_type = "pc_curl_curl_mlmg"
+pc_curl_curl_mlmg.verbose = true
+pc_curl_curl_mlmg.max_iter = 10
+pc_curl_curl_mlmg.relative_tolerance = 1e-4
+
+algo.particle_pusher = "boris"
+
+algo.particle_shape = 2
+algo.current_deposition = "villasenor"
+
+#################################
+############ PLASMA #############
+#################################
+particles.species_names = electrons deuterium
+
+electrons.charge = -q_e
+electrons.mass = m_e
+electrons.injection_style = "NUniformPerCell"
+electrons.num_particles_per_cell_each_dim = Nppc
+electrons.profile = constant
+electrons.density = n0
+electrons.zmax = R0
+electrons.momentum_distribution_type = "gaussian"
+electrons.ux_th = sqrt(T0*q_e/m_e)/clight
+electrons.uy_th = sqrt(T0*q_e/m_e)/clight
+electrons.uz_th = sqrt(T0*q_e/m_e)/clight
+
+deuterium.charge = q_e
+deuterium.mass = mD
+deuterium.injection_style = "NUniformPerCell"
+deuterium.num_particles_per_cell_each_dim = Nppc
+deuterium.profile = constant
+deuterium.density = n0
+deuterium.zmax = R0
+deuterium.momentum_distribution_type = "gaussian"
+deuterium.ux_th = sqrt(T0*q_e/mD)/clight
+deuterium.uy_th = sqrt(T0*q_e/mD)/clight
+deuterium.uz_th = sqrt(T0*q_e/mD)/clight
+
+###   Diagnostics
+diagnostics.diags_names = diag1
+diag1.intervals = 20
+diag1.diag_type = Full
+diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho divE
+diag1.electrons.variables = z w ux uy uz
+diag1.deuterium.variables = z w ux uy uz
+
+###   reduced diagnostics
+warpx.reduced_diags_names = particle_energy field_energy poynting_flux
+
+particle_energy.type = ParticleEnergy
+particle_energy.intervals = 10
+particle_energy.precision = 18
+particle_energy.path = diags/reduced_files/
+
+field_energy.type = FieldEnergy
+field_energy.intervals = 10
+field_energy.precision = 18
+field_energy.path = diags/reduced_files/
+
+poynting_flux.type = FieldPoyntingFlux
+poynting_flux.intervals = 10
+poynting_flux.precision = 18
+poynting_flux.path = diags/reduced_files/
+
+#################################
+############ COLLISION ##########
+#################################
+collisions.collision_names = collision1 collision2 collision3
+
+collision1.species = electrons electrons
+collision2.species = deuterium deuterium
+collision3.species = electrons deuterium
+collision1.CoulombLog = 10.0
+collision2.CoulombLog = 10.0
+collision3.CoulombLog = 10.0

Regression/Checksum/benchmarks_json/test_1d_planar_pinch_withPC.json

@@ -0,0 +1,29 @@
+{
+  "lev=0": {
+    "Bx": 0.6884832387863651,
+    "By": 0.6521976230101406,
+    "Bz": 0.0,
+    "Ex": 146532976.36091387,
+    "Ey": 149594679.4515056,
+    "Ez": 220206653.62773055,
+    "divE": 2463670936036.45,
+    "jx": 5671672121.997759,
+    "jy": 5780685336.896632,
+    "jz": 3889054995.637104,
+    "rho": 21.813805174064612
+  },
+  "deuterium": {
+    "particle_momentum_x": 4.004247985489721e-18,
+    "particle_momentum_y": 3.99862983197008e-18,
+    "particle_momentum_z": 3.985428948884101e-18,
+    "particle_position_x": 1622.111089442873,
+    "particle_weight": 1.5000033326127676e+21
+  },
+  "electrons": {
+    "particle_momentum_x": 6.586489162124888e-20,
+    "particle_momentum_y": 6.583939351411483e-20,
+    "particle_momentum_z": 6.595354090513763e-20,
+    "particle_position_x": 1622.1110484337041,
+    "particle_weight": 1.5000033326127676e+21
+  }
+}

Source/Evolve/WarpXEvolve.cpp

@@ -1139,18 +1139,20 @@ WarpX::doQEDEvents ()
 #endif
 
 void
-WarpX::PushParticlesandDeposit (amrex::Real cur_time, bool skip_current, PushType push_type)
+WarpX::PushParticlesandDeposit (amrex::Real cur_time, bool skip_current,
+                                bool deposit_mass_matrices, PushType push_type)
 {
     // Evolve particles to p^{n+1/2} and x^{n+1}
     // Deposit current, j^{n+1/2}
     for (int lev = 0; lev <= finest_level; ++lev) {
-        PushParticlesandDeposit(lev, cur_time, DtType::Full, skip_current, push_type);
+        PushParticlesandDeposit(lev, cur_time, DtType::Full, skip_current,
+                                deposit_mass_matrices, push_type);
     }
 }
 
 void
 WarpX::PushParticlesandDeposit (int lev, amrex::Real cur_time, DtType a_dt_type, bool skip_current,
-                               PushType push_type)
+                               bool deposit_mass_matrices, PushType push_type)
 {
     using ablastr::fields::Direction;
     using warpx::fields::FieldType;
@@ -1178,6 +1180,7 @@ WarpX::PushParticlesandDeposit (int lev, amrex::Real cur_time, DtType a_dt_type,
         dt[lev],
         a_dt_type,
         skip_current,
+        deposit_mass_matrices,
         push_type
     );
     if (! skip_current) {

Source/FieldSolver/ImplicitSolvers/ImplicitSolver.H

@@ -66,12 +66,6 @@ public:
                            amrex::Real  a_dt,
                            int          a_step ) = 0;
 
-
-    /**
-     * \brief Return pointer to MultiFab array for mass matrix
-     */
-    [[nodiscard]] virtual const amrex::Vector<amrex::Array<amrex::MultiFab*,3>>* GetSigmaCoeff() const { return nullptr; }
-
     //
     // the following routines are called by the linear and nonlinear solvers
     //
@@ -102,6 +96,14 @@ public:
     [[nodiscard]] amrex::Array<amrex::LinOpBCType,AMREX_SPACEDIM> GetLinOpBCLo () const;
     [[nodiscard]] amrex::Array<amrex::LinOpBCType,AMREX_SPACEDIM> GetLinOpBCHi () const;
 
+    /**
+     * \brief Return pointer to MultiFab array for mass matrix
+     */
+    [[nodiscard]] const amrex::Vector<amrex::Array<amrex::MultiFab*,3>>* GetMassMatricesCoeff() const {
+        if (m_use_mass_matrices) { return &m_mmpc_mfarrvec; }
+        else { return nullptr; }
+    }
+
     [[nodiscard]] amrex::Real  GetTheta () const { return m_theta; }
 
 protected:
@@ -146,6 +148,16 @@ protected:
      */
     int m_max_particle_iterations = 21;
 
+    /**
+     * \brief Whether to use mass matrices for the implicit solver
+     */
+    bool m_use_mass_matrices = false;
+
+    /**
+     * \brief Array of multifab pointers to mass matrix preconditioner
+     */
+    amrex::Vector<amrex::Array<amrex::MultiFab*,3>> m_mmpc_mfarrvec;
+
     /**
      * \brief parse nonlinear solver parameters (if one is used)
      */
@@ -166,6 +178,7 @@ protected:
             m_nlsolver = std::make_unique<NewtonSolver<WarpXSolverVec,ImplicitSolver>>();
             pp.query("max_particle_iterations", m_max_particle_iterations);
             pp.query("particle_tolerance", m_particle_tolerance);
+            pp.query("use_mass_matrices", m_use_mass_matrices);
         }
         else {
             WARPX_ABORT_WITH_MESSAGE(
@@ -174,6 +187,11 @@ protected:
 
     }
 
+    /**
+     * \brief Initialize the Mass Matrices used for plasma response in nonlinear Newton solver
+     */
+    void InitializeMassMatrices ();
+
     /**
      * \brief Convert from WarpX FieldBoundaryType to amrex::LinOpBCType
      */