Add PML Support for multi-J Algorithm

Examples/Tests/multi_J/inputs_2d_pml

@@ -0,0 +1,133 @@
+# Iterations
+max_step = 150
+
+# Domain decomposition
+amr.n_cell = 128 256
+warpx.numprocs = 1 2
+
+# Mesh refinement and geometry
+amr.max_level = 0
+geometry.coord_sys = 0
+geometry.prob_lo = -100e-6 -220e-6
+geometry.prob_hi =  100e-6   10e-6
+
+# Boundary condition
+boundary.field_lo = periodic pml
+boundary.field_hi = periodic pml
+
+# Algorithms
+algo.current_deposition = direct
+algo.field_gathering = energy-conserving
+algo.maxwell_solver = psatd
+algo.particle_pusher = vay
+algo.particle_shape = 3
+
+# Numerics
+warpx.cfl = 3.19
+warpx.do_nodal = 1
+warpx.use_filter = 1
+warpx.verbose = 1
+
+# Boosted frame
+warpx.boost_direction = z
+warpx.gamma_boost = 2.870114028490
+
+# Moving window
+warpx.do_moving_window = 1
+warpx.moving_window_dir = z
+warpx.moving_window_v = 1.
+
+# Spectral solver
+psatd.do_time_averaging = 0
+psatd.J_linear_in_time = 1
+psatd.update_with_rho = 1
+
+# Multi-J scheme
+warpx.do_multi_J = 1
+warpx.do_multi_J_n_depositions = 2
+warpx.do_dive_cleaning = 1
+warpx.do_divb_cleaning = 1
+
+# Particles
+particles.species_names = driver driver_back plasma_e plasma_p
+particles.use_fdtd_nci_corr = 0
+particles.rigid_injected_species = driver driver_back
+
+# Driver (electrons)
+driver.species_type = electron
+driver.injection_style = "gaussian_beam"
+driver.x_rms = 5e-6
+driver.y_rms = 5e-6
+driver.z_rms = 20.1e-6
+driver.x_m = 0.
+driver.y_m = 0.
+driver.z_m = -80e-6
+driver.npart = 100000
+driver.q_tot = -1e-10
+driver.momentum_distribution_type = "constant"
+driver.ux = 0.
+driver.uy = 0.
+driver.uz = 2e9
+driver.zinject_plane = 2e-3
+driver.rigid_advance = true
+driver.initialize_self_fields = 0
+driver.do_symmetrize = 1
+
+# Driver (positrons)
+driver_back.species_type = positron
+driver_back.injection_style = "gaussian_beam"
+driver_back.x_rms = 5e-6
+driver_back.y_rms = 5e-6
+driver_back.z_rms = 20.1e-6
+driver_back.x_m = 0.
+driver_back.y_m = 0.
+driver_back.z_m = -80e-6
+driver_back.npart = 100000
+driver_back.q_tot = 1e-10
+driver_back.momentum_distribution_type = "constant"
+driver_back.ux = 0.
+driver_back.uy = 0.
+driver_back.uz = 2e9
+driver_back.zinject_plane = 2e-3
+driver_back.rigid_advance = true
+driver_back.initialize_self_fields = 0
+driver_back.do_symmetrize = 1
+driver_back.do_backward_propagation = true
+
+# Electrons
+plasma_e.species_type = electron
+plasma_e.injection_style = "NUniformPerCell"
+plasma_e.zmin = 0.
+plasma_e.zmax = 0.05
+plasma_e.xmin = -90e-6
+plasma_e.xmax =  90e-6
+plasma_e.ymin = -90e-6
+plasma_e.ymax =  90e-6
+plasma_e.profile = constant
+plasma_e.density = 1e23
+plasma_e.num_particles_per_cell_each_dim = 1 1 1
+plasma_e.momentum_distribution_type = "at_rest"
+plasma_e.do_continuous_injection = 1
+
+# Hydrogen
+plasma_p.species_type = hydrogen
+plasma_p.injection_style = "NUniformPerCell"
+plasma_p.zmin = 0.
+plasma_p.zmax = 0.05
+plasma_p.xmin = -90e-6
+plasma_p.xmax =  90e-6
+plasma_p.ymin = -90e-6
+plasma_p.ymax =  90e-6
+plasma_p.profile = constant
+plasma_p.density = 1e23
+plasma_p.num_particles_per_cell_each_dim = 1 1 1
+plasma_p.momentum_distribution_type = "at_rest"
+plasma_p.do_continuous_injection = 1
+
+# Diagnostics
+diagnostics.diags_names = diag1
+diag1.intervals = 150
+diag1.diag_type = Full
+diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz F G divE rho rho_driver rho_driver_back rho_plasma_e rho_plasma_p
+diag1.write_species = 1
+diag1.species = driver plasma_e plasma_p

Regression/Checksum/benchmarks_json/multi_J_2d_psatd_pml.json

@@ -0,0 +1,51 @@
+{
+  "driver": {
+    "particle_cpu": 0.0,
+    "particle_id": 5000050000.0,
+    "particle_momentum_x": 1.4499934898508315e-16,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 9.822790382300442e-09,
+    "particle_position_x": 0.40000376757269557,
+    "particle_position_y": 30.100865423602578,
+    "particle_weight": 124830181489215.27
+  },
+  "lev=0": {
+    "Bx": 0.0,
+    "By": 921660.3791780534,
+    "Bz": 0.0,
+    "Ex": 270876373466331.25,
+    "Ey": 0.0,
+    "Ez": 92624942849290.58,
+    "F": 6183.089569602618,
+    "G": 0.0,
+    "divE": 2.519019708946362e+19,
+    "jx": 1.0172819575614014e+16,
+    "jy": 0.0,
+    "jz": 6.385273746740672e+16,
+    "rho": 220331681.67046365,
+    "rho_driver": 2562225.119933118,
+    "rho_driver_back": 0.0,
+    "rho_plasma_e": 1359499350.2687225,
+    "rho_plasma_p": 1361225460.3586755
+  },
+  "plasma_e": {
+    "particle_cpu": 29647.0,
+    "particle_id": 1194856401.0,
+    "particle_momentum_x": 7.187004715206369e-19,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 2.346497069529724e-17,
+    "particle_position_x": 1.3919451702129872,
+    "particle_position_y": 10.07681584795493,
+    "particle_weight": 6.641904330784497e+16
+  },
+  "plasma_p": {
+    "particle_cpu": 29696.0,
+    "particle_id": 1202173888.0,
+    "particle_momentum_x": 1.4494684060818139e-18,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 4.005707216626881e-14,
+    "particle_position_x": 1.3456480686756227,
+    "particle_position_y": 10.082154799552143,
+    "particle_weight": 6.652881944445523e+16
+  }
+}

Regression/WarpX-tests.ini

@@ -2218,6 +2218,24 @@ particleTypes = driver driver_back plasma_e plasma_p
 analysisRoutine = Examples/analysis_default_regression.py
 tolerance = 1e-14
 
+[multi_J_2d_psatd_pml]
+buildDir = .
+inputFile = Examples/Tests/multi_J/inputs_2d_pml
+runtime_params =
+dim = 2
+addToCompileString = USE_PSATD=TRUE
+restartTest = 0
+useMPI = 1
+numprocs = 2
+useOMP = 1
+numthreads = 1
+compileTest = 0
+doVis = 0
+compareParticles = 1
+particleTypes =
+analysisRoutine = Examples/analysis_default_regression.py
+tolerance = 1e-14
+
 [multi_J_rz_psatd]
 buildDir = .
 inputFile = Examples/Tests/multi_J/inputs_rz

Source/Evolve/WarpXEvolve.cpp

@@ -10,6 +10,7 @@
  */
 #include "WarpX.H"
 
+#include "BoundaryConditions/PML.H"
 #include "Diagnostics/BackTransformedDiagnostic.H"
 #include "Diagnostics/MultiDiagnostics.H"
 #include "Diagnostics/ReducedDiags/MultiReducedDiags.H"
@@ -610,23 +611,50 @@ WarpX::OneStep_multiJ (const amrex::Real cur_time)
             }
         }
 
-        // Transform fields back to real space and exchange guard cells
+        // Transform fields back to real space
         if (WarpX::fft_do_time_averaging)
         {
             // We summed the integral of the field over 2*dt
             PSATDScaleAverageFields(1._rt / (2._rt*dt[0]));
             PSATDBackwardTransformEBavg();
         }
+
+        // Evolve fields in PML
+        for (int lev = 0; lev <= finest_level; ++lev)
+        {
+            if (do_pml && pml[lev]->ok())
+            {
+                pml[lev]->PushPSATD(lev);
+            }
+            ApplyEfieldBoundary(lev, PatchType::fine);
+            if (lev > 0) ApplyEfieldBoundary(lev, PatchType::coarse);
+            ApplyBfieldBoundary(lev, PatchType::fine, DtType::FirstHalf);
+            if (lev > 0) ApplyBfieldBoundary(lev, PatchType::coarse, DtType::FirstHalf);
+        }
+
+        // Damp fields in PML before exchanging guard cells
+        if (do_pml)
+        {
+            DampPML();
+        }
+
+        // Exchange guard cells
         FillBoundaryE(guard_cells.ng_alloc_EB);
         FillBoundaryB(guard_cells.ng_alloc_EB);
-        if (WarpX::do_dive_cleaning) FillBoundaryF(guard_cells.ng_alloc_F);
-        if (WarpX::do_divb_cleaning) FillBoundaryG(guard_cells.ng_alloc_G);
+        if (WarpX::do_dive_cleaning || WarpX::do_pml_dive_cleaning) FillBoundaryF(guard_cells.ng_alloc_F);
+        if (WarpX::do_divb_cleaning || WarpX::do_pml_divb_cleaning) FillBoundaryG(guard_cells.ng_alloc_G);
 
         // Synchronize E, B, F, G fields on nodal points
         NodalSync(Efield_fp, Efield_cp);
         NodalSync(Bfield_fp, Bfield_cp);
         if (WarpX::do_dive_cleaning) NodalSync(F_fp, F_cp);
         if (WarpX::do_divb_cleaning) NodalSync(G_fp, G_cp);
+
+        // Synchronize fields on nodal points in PML
+        if (do_pml)
+        {
+            NodalSyncPML();
+        }
     }
     else
     {

Source/WarpX.cpp

@@ -745,11 +745,6 @@ WarpX::ReadParameters ()
         pp_warpx.query("do_pml_j_damping", do_pml_j_damping);
         pp_warpx.query("do_pml_in_domain", do_pml_in_domain);
 
-        if (do_multi_J && isAnyBoundaryPML())
-        {
-            amrex::Abort("Multi-J algorithm not implemented with PMLs");
-        }
-
         // Default values of WarpX::do_pml_dive_cleaning and WarpX::do_pml_divb_cleaning:
         // false for FDTD solver, true for PSATD solver.
         if (maxwell_solver_id != MaxwellSolverAlgo::PSATD)