AMR test case: Tracking example.

modified:   ippl/test/AMR/AmrPartBunch.h
modified:   ippl/test/AMR/CMakeLists.txt
new file:   ippl/test/AMR/testTracker.cpp

ippl/test/AMR/AmrPartBunch.h

@@ -153,6 +153,12 @@ public:
         }
     }
     
+    int getLevel(int i) {
+        int l, g, dq;
+        std::tie(l,g,dq) = idxMap_m[i];
+        return l;
+    }
+    
     Vector_t interpolate(int i, MultiFab& quantity) {
         
         int lev, grid, dq;
@@ -173,18 +179,18 @@ public:
         }
         
         const FArrayBox& gfab = (*ac_pointer)[grid];
-        Real grav[1] = { 0.0 };
+        Real grav[3] = { 0.0, 0.0, 0.0 };
         
-        int idx[1] = { 0 };
+        int idx[3] = { 0, 1, 2 };
 
     ParticleBase::Interp(m_particles[lev][grid][dq],
                          m_gdb->Geom(lev),
                          gfab,
                          idx,
                          grav,
-                         1);
+                         BL_SPACEDIM);
         
-        return Vector_t(grav[0], 0, 0);
+        return Vector_t(grav[0], grav[1], grav[2]);
     }
     
 private:

ippl/test/AMR/CMakeLists.txt

@@ -87,6 +87,8 @@ IF (ENABLE_AMR_SOLVER)
     
     add_executable(toFile toFile.cpp H5Reader.cpp Distribution.cpp)
     
+    add_executable (testTracker testTracker.cpp AmrOpal.cpp AmrPartBunch.cpp H5Reader.cpp Distribution.cpp ${F90_source_files} ${CMAKE_CURRENT_SOURCE_DIR}/../../../src/Classic/Physics/Physics.cpp Solver.cpp)
+    
     # Boxlib has a circular dependency between -lcboxlib and -lfboxlib --> resolve by: -lcboxlib -lfboxlib -lcboxlib
 #     target_link_libraries (testAmrPartBunch ${LIBS} -lgfortran ${CCSE_LIBRARIES} -lcboxlib ${OTHER_CMAKE_EXE_LINKER_FLAGS} ${CMAKE_DL_LIBS})
 #     target_link_libraries (testSolver ${LIBS} -lgfortran ${CCSE_LIBRARIES} -lcboxlib ${OTHER_CMAKE_EXE_LINKER_FLAGS} ${CMAKE_DL_LIBS})
@@ -107,4 +109,6 @@ IF (ENABLE_AMR_SOLVER)
     
     target_link_libraries (toFile ${LIBS} -lgfortran ${CCSE_LIBRARIES} -lcboxlib ${OTHER_CMAKE_EXE_LINKER_FLAGS} ${CMAKE_DL_LIBS})
     
+    target_link_libraries (testTracker ${LIBS} -lgfortran ${CCSE_LIBRARIES} -lcboxlib ${OTHER_CMAKE_EXE_LINKER_FLAGS} ${CMAKE_DL_LIBS})
+    
 ENDIF(ENABLE_AMR_SOLVER)

ippl/test/AMR/testTracker.cpp

+/*!
+ * @file testTracker.cpp
+ * @author Matthias Frey
+ * @date November 2016
+ * @details Compute \f$\Delta\phi = -\rho\f$ where the charge
+ * density is a Gaussian distribution. Write plot files
+ * that can be visualized by yt (python visualize.py)
+ * or by AmrVis of the CCSE group at LBNL.
+ * 
+ * Domain:  [-0.5, 0.5] x [-0.5, 0.5] x [-0.5, 0.5]\n
+ * BC:      Dirichlet boundary conditions\n
+ * Charge/particle: elementary charge\n
+ * Gaussian particle distribution N(0.0, 0.01)
+ * 
+ * Call:\n
+ *  mpirun -np [#cores] testTracker [#gridpoints x] [#gridpoints y] [#gridpoints z]
+ *                                  [#particles] [#levels] [max. box size]
+ * 
+ * @brief Computes \f$\Delta\phi = -\rho\f$
+ */
+
+#include "Ippl.h"
+#include <string>
+#include <fstream>
+#include <vector>
+#include <iostream>
+#include <set>
+#include <sstream>
+#include <iomanip>
+
+
+#include <ParmParse.H>
+
+#include "PartBunch.h"
+#include "AmrPartBunch.h"
+
+#include "Distribution.h"
+#include "Solver.h"
+#include "AmrOpal.h"
+
+#include "helper_functions.h"
+
+#include "writePlotFile.H"
+
+#include <cmath>
+
+#include "Physics/Physics.h"
+
+void doSolve(AmrOpal& myAmrOpal, PartBunchBase* bunch,
+             container_t& rhs,
+             container_t& phi,
+             container_t& grad_phi,
+             const Array<Geometry>& geom,
+             const Array<int>& rr,
+             int nLevels,
+             Inform& msg)
+{
+    static IpplTimings::TimerRef allocTimer = IpplTimings::getTimer("alloc-memory-grid");
+    
+    static IpplTimings::TimerRef assignTimer = IpplTimings::getTimer("assign-charge");
+    
+    // =======================================================================                                                                                                                                   
+    // 4. prepare for multi-level solve                                                                                                                                                                          
+    // =======================================================================
+    
+    rhs.resize(nLevels);
+    phi.resize(nLevels);
+    grad_phi.resize(nLevels);
+    
+    IpplTimings::startTimer(allocTimer);
+    
+    for (int lev = 0; lev < nLevels; ++lev) {
+        initGridData(rhs, phi, grad_phi, myAmrOpal.boxArray()[lev], lev);
+    }
+    
+    IpplTimings::stopTimer(allocTimer);
+
+    // Define the density on level 0 from all particles at all levels                                                                                                                                            
+    int base_level   = 0;
+    int finest_level = myAmrOpal.finestLevel();
+
+    IpplTimings::startTimer(assignTimer);
+    dynamic_cast<AmrPartBunch*>(bunch)->AssignDensity(0, false, rhs, base_level, 1, finest_level);
+    IpplTimings::stopTimer(assignTimer);
+    
+    double totCharge = totalCharge(rhs, finest_level, geom);
+    msg << "Total Charge: " << totCharge << " C" << endl
+        << "Vacuum permittivity: " << Physics::epsilon_0 << " F/m (= C/(m V)" << endl;
+    Real vol = (*(geom[0].CellSize()) * *(geom[0].CellSize()) * *(geom[0].CellSize()) );
+    msg << "Cell volume: " << vol << " m^3" << endl;
+    
+    // eps in C / (V * m)
+    double constant = -1.0 / Physics::epsilon_0;  // in [V m / C]
+    for (int i = 0; i <=finest_level; ++i) {
+#ifdef UNIQUE_PTR
+        rhs[i]->mult(constant, 0, 1);       // in [V m]
+#else
+        rhs[i].mult(constant, 0, 1);
+#endif
+    }
+    
+    // **************************************************************************                                                                                                                                
+    // Compute the total charge of all particles in order to compute the offset                                                                                                                                  
+    //     to make the Poisson equations solvable                                                                                                                                                                
+    // **************************************************************************                                                                                                                                
+
+    Real offset = 0.;
+
+    // solve                                                                                                                                                                                                     
+    Solver sol;
+    sol.solve_for_accel(rhs,
+                        phi,
+                        grad_phi,
+                        geom,
+                        base_level,
+                        finest_level,
+                        offset);
+    
+    // for plotting unnormalize
+    for (int i = 0; i <=finest_level; ++i) {
+#ifdef UNIQUE_PTR
+        rhs[i]->mult(1.0 / constant, 0, 1);       // in [V m]
+#else
+        rhs[i].mult(1.0 / constant, 0, 1);
+#endif
+    }
+}
+
+void doBoxLib(const Vektor<size_t, 3>& nr, size_t nParticles,
+              int nLevels, size_t maxBoxSize, Inform& msg)
+{
+    static IpplTimings::TimerRef distTimer = IpplTimings::getTimer("init-distribution");
+    static IpplTimings::TimerRef regridTimer = IpplTimings::getTimer("regrid");
+    static IpplTimings::TimerRef redistTimer = IpplTimings::getTimer("particle-redistr");
+    // ========================================================================
+    // 1. initialize physical domain (just single-level)
+    // ========================================================================
+    
+    std::array<double, BL_SPACEDIM> lower = {{-0.5, -0.5, -0.5}}; // m
+    std::array<double, BL_SPACEDIM> upper = {{ 0.5,  0.5,  0.5}}; // m
+    
+    RealBox domain;
+    
+    // in helper_functions.h
+    init(domain, nr, lower, upper);
+    
+    
+    /*
+     * create an Amr object
+     */
+    ParmParse pp("amr");
+    pp.add("max_grid_size", int(maxBoxSize));
+    
+    Array<int> error_buf(nLevels, 0);
+    
+    pp.addarr("n_error_buf", error_buf);
+    pp.add("grid_eff", 0.95);
+    
+    Array<int> nCells(3);
+    for (int i = 0; i < 3; ++i)
+        nCells[i] = nr[i];
+    
+    AmrOpal myAmrOpal(&domain, nLevels - 1, nCells, 0 /* cartesian */);
+    
+    
+    // ========================================================================
+    // 2. initialize all particles (just single-level)
+    // ========================================================================
+    
+    const Array<BoxArray>& ba = myAmrOpal.boxArray();
+    const Array<DistributionMapping>& dmap = myAmrOpal.DistributionMap();
+    const Array<Geometry>& geom = myAmrOpal.Geom();
+    
+    Array<int> rr(nLevels);
+    for (int i = 0; i < nLevels; ++i)
+        rr[i] = 2;
+    
+    PartBunchBase* bunch = new AmrPartBunch(geom, dmap, ba, rr);
+    
+    
+    // initialize a particle distribution
+    unsigned long int nloc = nParticles / ParallelDescriptor::NProcs();
+    Distribution dist;
+    IpplTimings::startTimer(distTimer);
+    double sig = 0.005;  // m
+    dist.gaussian(0.0, sig, nloc, ParallelDescriptor::MyProc());
+    
+    // copy particles to the PartBunchBase object.
+    dist.injectBeam(*bunch);
+    IpplTimings::stopTimer(distTimer);
+    
+    
+    
+    for (std::size_t i = 0; i < bunch->getLocalNum(); ++i) {
+        bunch->setQM(1.0e-14, i);  // in [C]
+        bunch->setP(1.0e-2, i);
+    }
+    
+    // redistribute on single-level
+    IpplTimings::startTimer(redistTimer);
+    bunch->myUpdate();
+    IpplTimings::stopTimer(redistTimer);
+    
+    msg << "Single-level statistics" << endl;
+    bunch->gatherStatistics();
+    
+    msg << "Bunch radius: " << sig << " m" << endl
+        << "#Particles: " << nParticles << endl
+        << "Charge per particle: " << bunch->getQM(0) << " C" << endl
+        << "Total charge: " << nParticles * bunch->getQM(0) << " C" << endl;
+    
+    
+    myAmrOpal.setBunch(dynamic_cast<AmrPartBunch*>(bunch));
+        
+    const Array<Geometry>& geoms = myAmrOpal.Geom();
+    for (int i = 0; i < nLevels; ++i)
+        msg << "#Cells per dim of level " << i << " for bunch : "
+            << 2.0 * sig / *(geoms[i].CellSize()) << endl;
+    
+    // ========================================================================
+    // 2. tagging (i.e. create BoxArray's, DistributionMapping's of all
+    //    other levels)
+    // ========================================================================
+    
+    /*
+     * do tagging
+     */
+//     dynamic_cast<AmrPartBunch*>(bunch)->Define (myAmrOpal.Geom(),
+//                                                 myAmrOpal.DistributionMap(),
+//                                                 myAmrOpal.boxArray(),
+//                                                 rr);
+//     
+//     
+    // ========================================================================
+    // 3. multi-level redistribute
+    // ========================================================================
+    IpplTimings::startTimer(regridTimer);
+    for (int i = 0; i <= myAmrOpal.finestLevel() && i < myAmrOpal.maxLevel(); ++i)
+        myAmrOpal.regrid(i /*lbase*/, 0.0 /*time*/);
+    IpplTimings::stopTimer(regridTimer);
+    
+    msg << "Multi-level statistics" << endl;
+    bunch->gatherStatistics();
+    
+    container_t rhs;
+    container_t phi;
+    container_t grad_phi;
+    
+    std::string plotsolve = BoxLib::Concatenate("plt", 0, 4);
+    
+    double m = 1.0; // mass
+    double gam = 1.5; // relativistic factor
+    double dt = 0.0005;
+    for (int t = 0; t < 100; ++t) {
+        
+        doSolve(myAmrOpal, bunch, rhs, phi, grad_phi, geoms, rr, nLevels, msg);
+        
+        for (std::size_t i = 0; i < bunch->getLocalNum(); ++i) {
+            int level = dynamic_cast<AmrPartBunch*>(bunch)->getLevel(i);
+            
+            // update momentum half step
+            Vector_t force = dynamic_cast<AmrPartBunch*>(bunch)->interpolate(i, grad_phi[level]);
+            bunch->setP(bunch->getP(i) + 0.5 * dt * m * gam * force, i);
+            
+            // update position full step
+            bunch->setR(bunch->getR(i) + dt / ( m * gam ) * bunch->getP(i), i);
+        }
+        
+        doSolve(myAmrOpal, bunch, rhs, phi, grad_phi, geoms, rr, nLevels, msg);
+        
+        for (std::size_t i = 0; i < bunch->getLocalNum(); ++i) {
+            int level = dynamic_cast<AmrPartBunch*>(bunch)->getLevel(i);
+            Vector_t force = dynamic_cast<AmrPartBunch*>(bunch)->interpolate(i, grad_phi[level]);
+            
+            // update momentum half step
+            bunch->setP(bunch->getP(i) + 0.5 * dt * m * gam * force, i);
+        }
+        
+        msg << "Total field energy: " << totalFieldEnergy(grad_phi, rr) << endl;
+    
+        for (int i = 0; i <= myAmrOpal.finestLevel(); ++i) {
+#ifdef UNIQUE_PTR
+            msg << "Max. potential level " << i << ": "<< phi[i]->max(0) << endl
+                << "Min. potential level " << i << ": " << phi[i]->min(0) << endl
+                << "Max. ex-field level " << i << ": " << grad_phi[i]->max(0) << endl
+                << "Min. ex-field level " << i << ": " << grad_phi[i]->min(0) << endl;
+#else
+            msg << "Max. potential level " << i << ": "<< phi[i].max(0) << endl
+                << "Min. potential level " << i << ": " << phi[i].min(0) << endl
+                << "Max. ex-field level " << i << ": " << grad_phi[i].max(0) << endl
+                << "Min. ex-field level " << i << ": " << grad_phi[i].min(0) << endl;
+#endif
+        }
+        
+        for (int i = 0; i <= myAmrOpal.finestLevel() && i < myAmrOpal.maxLevel(); ++i)
+            myAmrOpal.regrid(i /*lbase*/, 0.0 /*time*/);
+        
+//         writePlotFile(plotsolve, rhs, phi, grad_phi, rr, geoms, t);
+        dynamic_cast<AmrPartBunch*>(bunch)->python_format(t);
+    }
+    
+}
+
+
+int main(int argc, char *argv[]) {
+    
+    Ippl ippl(argc, argv);
+    
+    Inform msg("Solver");
+    
+
+    static IpplTimings::TimerRef mainTimer = IpplTimings::getTimer("main");
+    IpplTimings::startTimer(mainTimer);
+
+    std::stringstream call;
+    call << "Call: mpirun -np [#procs] " << argv[0]
+         << " [#gridpoints x] [#gridpoints y] [#gridpoints z] [#particles] "
+         << "[#levels] [max. box size]";
+    
+    if ( argc < 7 ) {
+        msg << call.str() << endl;
+        return -1;
+    }
+    
+    // number of grid points in each direction
+    Vektor<size_t, 3> nr(std::atoi(argv[1]),
+                         std::atoi(argv[2]),
+                         std::atoi(argv[3]));
+    
+    
+    size_t nParticles = std::atoi(argv[4]);
+    
+    
+    msg << "Particle test running with" << endl
+        << "- #particles = " << nParticles << endl
+        << "- grid       = " << nr << endl;
+        
+    BoxLib::Initialize(argc,argv, false);
+    size_t nLevels = std::atoi(argv[5]) + 1; // i.e. nLevels = 0 --> only single level
+    size_t maxBoxSize = std::atoi(argv[6]);
+    doBoxLib(nr, nParticles, nLevels, maxBoxSize, msg);
+    
+    
+    IpplTimings::stopTimer(mainTimer);
+
+    IpplTimings::print();
+    
+    std::stringstream timefile;
+    timefile << std::string(argv[0]) << "-timing-cores-"
+             << std::setfill('0') << std::setw(6) << Ippl::getNodes()
+             << "-threads-1.dat";
+    
+    IpplTimings::print(timefile.str());
+    
+    return 0;
+}