Add ForceFromSelfFields

Add function ForceFromSelfFields to calculate space charge self
fields and their force.

src/ImpactX.cpp

@@ -13,6 +13,7 @@
 #include "particles/Push.H"
 #include "particles/diagnostics/DiagnosticOutput.H"
 #include "particles/spacecharge/PoissonSolve.H"
+#include "particles/spacecharge/ForceFromSelfFields.H"
 #include "particles/transformation/CoordinateTransformation.H"
 
 #include <AMReX.H>
@@ -148,9 +149,10 @@ namespace impactx
                     // poisson solve in x,y,z
                     spacecharge::PoissonSolve(*m_particle_container, m_phi, m_rho);
 
-                    // calculate force
-                    // TODO: FDTD stencil m_phi -> m_space_charge
-                    // nodal: (-1 and +1) / 2dx
+                    // calculate force in x,y,z
+                    spacecharge::ForceFromSelfFields(m_space_charge_force,
+                                                     m_phi,
+                                                     this->geom);
 
                     // gather and space-charge push in x,y,z , assuming the space-charge
                     // field is the same before/after transformation

src/particles/spacecharge/CMakeLists.txt

@@ -1,4 +1,5 @@
 target_sources(ImpactX
   PRIVATE
     PoissonSolve.cpp
+    ForceFromSelfFields.cpp
 )

src/particles/spacecharge/ForceFromSelfFields.H

@@ -0,0 +1,41 @@
+/* Copyright 2022 The Regents of the University of California, through Lawrence
+ *           Berkeley National Laboratory (subject to receipt of any required
+ *           approvals from the U.S. Dept. of Energy). All rights reserved.
+ *
+ * This file is part of ImpactX.
+ *
+ * Authors: Marco Garten, Axel Huebl
+ * License: BSD-3-Clause-LBNL
+ */
+#ifndef IMPACTX_FORCEFROMSELFFIELDS_H
+#define IMPACTX_FORCEFROMSELFFIELDS_H
+
+#include "particles/ImpactXParticleContainer.H"
+
+#include <AMReX_Geometry.H>
+#include <AMReX_MultiFab.H>
+#include <AMReX_Vector.H>
+
+#include <unordered_map>
+
+
+namespace impactx::spacecharge
+{
+    /** Calculate the space charge force field from the electric potential
+     *
+     * This resets the values in scf_<component> to zero and then calculates the space
+     * charge force field.
+     *
+     * @param[inout] space_charge_force space charge force component in x,y,z per level
+     * @param[in] phi scalar potential per level
+     * @param[in] geom geometry object
+     */
+    void ForceFromSelfFields (
+        std::unordered_map<int, std::unordered_map<std::string, amrex::MultiFab> > & space_charge_force,
+        std::unordered_map<int, amrex::MultiFab> const & phi,
+        const amrex::Vector<amrex::Geometry>& geom
+    );
+
+} // namespace impactx
+
+#endif // IMPACTX_FORCEFROMSELFFIELDS_H

src/particles/spacecharge/ForceFromSelfFields.cpp

@@ -0,0 +1,62 @@
+/* Copyright 2022 The Regents of the University of California, through Lawrence
+ *           Berkeley National Laboratory (subject to receipt of any required
+ *           approvals from the U.S. Dept. of Energy). All rights reserved.
+ *
+ * This file is part of ImpactX.
+ *
+ * Authors: Marco Garten, Axel Huebl
+ * License: BSD-3-Clause-LBNL
+ */
+#include "ForceFromSelfFields.H"
+
+#include <ablastr/fields/PoissonSolver.H>
+
+#include <AMReX_BLProfiler.H>
+#include <AMReX_REAL.H>       // for Real
+#include <AMReX_SPACE.H>      // for AMREX_D_DECL
+
+
+namespace impactx::spacecharge
+{
+    void ForceFromSelfFields (
+        std::unordered_map<int, std::unordered_map<std::string, amrex::MultiFab> > & space_charge_force,
+        std::unordered_map<int, amrex::MultiFab> const & phi,
+        amrex::Vector<amrex::Geometry> const & geom
+    )
+    {
+        BL_PROFILE("impactx::spacecharge::ForceFromSelfFields");
+
+        using namespace amrex::literals;
+
+        // loop over refinement levels
+        int const finest_level = phi.size() - 1u;
+        for (int lev = 0; lev <= finest_level; ++lev) {
+
+            // stencil coefficients: 0.5 * inverse cell size
+            auto const &gm = geom[lev];
+            auto const dr = gm.CellSizeArray();
+            amrex::GpuArray<amrex::Real, 3> const inv2dr{AMREX_D_DECL(0.5_rt/dr[0], 0.5_rt/dr[1], 0.5_rt/dr[2])};
+
+            // reset the values in space_charge_force to zero
+            space_charge_force.at(lev).at("x").setVal(0.);
+            space_charge_force.at(lev).at("x").setVal(0.);
+            space_charge_force.at(lev).at("x").setVal(0.);
+
+            for (amrex::MFIter mfi(phi.at(lev)); mfi.isValid(); ++mfi) {
+
+                amrex::Box bx = mfi.validbox();
+                auto const phi_arr = (phi.at(lev))[mfi].const_array();
+
+                auto scf_arr_x = space_charge_force[lev]["x"][mfi].array();
+                auto scf_arr_y = space_charge_force[lev]["y"][mfi].array();
+                auto scf_arr_z = space_charge_force[lev]["z"][mfi].array();
+
+                amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
+                    scf_arr_x(i, j, k) = inv2dr[0] * (phi_arr(i-1, j, k) - phi_arr(i+1, j, k));
+                    scf_arr_y(i, j, k) = inv2dr[1] * (phi_arr(i, j-1, k) - phi_arr(i, j+1, k));
+                    scf_arr_z(i, j, k) = inv2dr[2] * (phi_arr(i, j, k-1) - phi_arr(i, j, k+1));
+                });
+            }
+        }
+    }
+} // namespace impactx::spacecharge