Merge branch '182-open-solver-should-work-with-origin-0' into 'master'

Resolve "Open solver should work with origin != 0" Closes #182 See merge request OPAL/Libraries/ippl!198
IPPL-framework · Jul 21, 2023 · 092500b · 092500b
2 parents 8e44641 + f152885
commit 092500b
Show file tree

Hide file tree

Showing 2 changed files with 14 additions and 23 deletions.
diff --git a/alpine/PenningTrap.cpp b/alpine/PenningTrap.cpp
@@ -199,11 +199,12 @@ int main(int argc, char* argv[]) {
         // create mesh and layout objects for this problem domain
         Vector_t<double, Dim> rmin = 0;
         Vector_t<double, Dim> rmax = 20;
+        Vector_t<double, Dim> length = rmax - rmin;
 
-        Vector_t<double, Dim> hr     = rmax / nr;
+        Vector_t<double, Dim> hr     = length / nr;
         Vector_t<double, Dim> origin = rmin;
         unsigned int nrMax           = 2048;  // Max grid size in our studies
-        double dxFinest              = rmax[0] / nrMax;
+        double dxFinest              = length[0] / nrMax;
         const double dt              = 0.5 * dxFinest;  // size of timestep
 
         const bool isAllPeriodic = true;
@@ -218,12 +219,10 @@ int main(int argc, char* argv[]) {
 
         P->nr_m = nr;
 
-        Vector_t<double, Dim> length = rmax - rmin;
-
         Vector_t<double, Dim> mu, sd;
 
         for (unsigned d = 0; d < Dim; d++) {
-            mu[d] = 0.5 * length[d];
+            mu[d] = 0.5 * length[d] + origin[d];
         }
         sd[0] = 0.15 * length[0];
         sd[1] = 0.05 * length[1];
@@ -350,15 +349,15 @@ int main(int argc, char* argv[]) {
             auto Rview = P->R.getView();
             auto Pview = P->P.getView();
             auto Eview = P->E.getView();
-            double V0  = 30 * rmax[2];
+            double V0  = 30 * length[2];
             Kokkos::parallel_for(
                 "Kick1", P->getLocalNum(), KOKKOS_LAMBDA(const size_t j) {
                     double Eext_x =
-                        -(Rview(j)[0] - 0.5 * rmax[0]) * (V0 / (2 * Kokkos::pow(rmax[2], 2)));
+                        -(Rview(j)[0] - origin[0] - 0.5 * length[0]) * (V0 / (2 * Kokkos::pow(length[2], 2)));
                     double Eext_y =
-                        -(Rview(j)[1] - 0.5 * rmax[1]) * (V0 / (2 * Kokkos::pow(rmax[2], 2)));
+                        -(Rview(j)[1] - origin[1] - 0.5 * length[1]) * (V0 / (2 * Kokkos::pow(length[2], 2)));
                     double Eext_z =
-                        (Rview(j)[2] - 0.5 * rmax[2]) * (V0 / (Kokkos::pow(rmax[2], 2)));
+                        (Rview(j)[2] - origin[2] - 0.5 * length[2]) * (V0 / (Kokkos::pow(length[2], 2)));
 
                     Eview(j)[0] += Eext_x;
                     Eview(j)[1] += Eext_y;
@@ -410,11 +409,11 @@ int main(int argc, char* argv[]) {
             Kokkos::parallel_for(
                 "Kick2", P->getLocalNum(), KOKKOS_LAMBDA(const size_t j) {
                     double Eext_x =
-                        -(R2view(j)[0] - 0.5 * rmax[0]) * (V0 / (2 * Kokkos::pow(rmax[2], 2)));
+                        -(R2view(j)[0] - origin[0] - 0.5 * length[0]) * (V0 / (2 * Kokkos::pow(length[2], 2)));
                     double Eext_y =
-                        -(R2view(j)[1] - 0.5 * rmax[1]) * (V0 / (2 * Kokkos::pow(rmax[2], 2)));
+                        -(R2view(j)[1] - origin[1] - 0.5 * length[1]) * (V0 / (2 * Kokkos::pow(length[2], 2)));
                     double Eext_z =
-                        (R2view(j)[2] - 0.5 * rmax[2]) * (V0 / (Kokkos::pow(rmax[2], 2)));
+                        (R2view(j)[2] - origin[2] - 0.5 * length[2]) * (V0 / (Kokkos::pow(length[2], 2)));
 
                     E2view(j)[0] += Eext_x;
                     E2view(j)[1] += Eext_y;

diff --git a/src/Solver/FFTPoissonSolver.hpp b/src/Solver/FFTPoissonSolver.hpp
@@ -235,17 +235,9 @@ namespace ippl {
         layout_mp = &(this->rhs_mp->getLayout());
         mesh_mp   = &(this->rhs_mp->get_mesh());
 
-        // get mesh spacing
-        hr_m = mesh_mp->getMeshSpacing();
-
-        // get origin
-        vector_type origin    = mesh_mp->getOrigin();
-        const scalar_type sum = std::abs(origin[0]) + std::abs(origin[1]) + std::abs(origin[2]);
-
-        // origin should always be 0 for Green's function computation to work...
-        if (sum != 0.0) {
-            throw IpplException("FFTPoissonSolver::initializeFields", "Origin is not 0");
-        }
+        // get mesh spacing and origin
+        hr_m               = mesh_mp->getMeshSpacing();
+        vector_type origin = mesh_mp->getOrigin();
 
         // create domain for the real fields
         domain_m = layout_mp->getDomain();