Fix electrostatic solver (#730)

* Create inputs

* fixed stuff

* Make smaller test

* Re-use EvolveEM for ES calculation

* Use correct variable const_dt

* Avoid out of bound and modify computation of electrostatic field

* Update test file

* Remove tab

* Minor refactoring of space-charge calculation

* Fix electrostatic test

* Rename WarpXEvolveEM to WarpXEvolve

* Use if conditions instead of ifdef

* Correct initialization of `do_electrostatic`

* Fix static variable

* Add documentation for `do_electrostatic`

Co-authored-by: dbizzozero <dbizzozero@users.noreply.github.com>

1 files changed, 34 insertions, 31 deletions

diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
index 159cb7b87..9ba783173 100644
--- a/Source/Particles/PhysicalParticleContainer.cpp
+++ b/Source/Particles/PhysicalParticleContainer.cpp
@@ -1316,41 +1316,44 @@ PhysicalParticleContainer::Evolve (int lev,
                 WARPX_PROFILE_VAR_STOP(blp_ppc_pp);
 
                 //
-                // Current Deposition
+                // Current Deposition (only needed for electromagnetic solver)
                 //
-
-                int* AMREX_RESTRICT ion_lev;
-                if (do_field_ionization){
-                    ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr();
-                } else {
-                    ion_lev = nullptr;
-                }
-                // Deposit inside domains
-                DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, &jx, &jy, &jz,
-                               0, np_current, thread_num,
-                               lev, lev, dt);
-
-                if (has_buffer){
-                    // Deposit in buffers
-                    DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, cjx, cjy, cjz,
-                                   np_current, np-np_current, thread_num,
-                                   lev, lev-1, dt);
-                }
-            }
+                if (!WarpX::do_electrostatic) {
+                    int* AMREX_RESTRICT ion_lev;
+                    if (do_field_ionization){
+                        ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr();
+                    } else {
+                        ion_lev = nullptr;
+                    }
+                    // Deposit inside domains
+                    DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, &jx, &jy, &jz,
+                                   0, np_current, thread_num,
+                                   lev, lev, dt);
+                    if (has_buffer){
+                        // Deposit in buffers
+                        DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, cjx, cjy, cjz,
+                                       np_current, np-np_current, thread_num,
+                                       lev, lev-1, dt);
+                    }
+                } // end of "if !do_electrostatic"
+            } // end of "if do_not_push"
 
             if (rho) {
                 // Deposit charge after particle push, in component 1 of MultiFab rho.
-                int* AMREX_RESTRICT ion_lev;
-                if (do_field_ionization){
-                    ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr();
-                } else {
-                    ion_lev = nullptr;
-                }
-                DepositCharge(pti, wp, ion_lev, rho, 1, 0,
-                              np_current, thread_num, lev, lev);
-                if (has_buffer){
-                    DepositCharge(pti, wp, ion_lev, crho, 1, np_current,
-                                  np-np_current, thread_num, lev, lev-1);
+                // (Skipped for electrostatic solver, as this may lead to out-of-bounds)
+                if (!WarpX::do_electrostatic) {
+                    int* AMREX_RESTRICT ion_lev;
+                    if (do_field_ionization){
+                        ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr();
+                    } else {
+                        ion_lev = nullptr;
+                    }
+                    DepositCharge(pti, wp, ion_lev, rho, 1, 0,
+                                  np_current, thread_num, lev, lev);
+                    if (has_buffer){
+                        DepositCharge(pti, wp, ion_lev, crho, 1, np_current,
+                                      np-np_current, thread_num, lev, lev-1);
+                    }
                 }
             }