AMR test case: Add two stream instability test.

modified:   ippl/test/AMR/AmrPartBunch.cpp
modified:   ippl/test/AMR/AmrPartBunch.h
modified:   ippl/test/AMR/CMakeLists.txt
modified:   ippl/test/AMR/Distribution.cpp
modified:   ippl/test/AMR/Distribution.h
modified:   ippl/test/AMR/PartBunchBase.h
modified:   ippl/test/AMR/Solver.cpp
new file:   ippl/test/AMR/twostream.cpp

ippl/test/AMR/AmrPartBunch.cpp

@@ -5,7 +5,7 @@
 // ----------------------------------------------------------------------------
 // STATIC MEMBER VARIABLES
 // ----------------------------------------------------------------------------
-size_t AmrPartBunch::nAttributes = 4;
+size_t AmrPartBunch::nAttributes = 5;
 
 // ----------------------------------------------------------------------------
 // PUBLIC MEMBER FUNCTIONS
@@ -209,6 +209,13 @@ double AmrPartBunch::getQM(int i) {
     return m_particles[l][g][dq].m_data[0];
 }
 
+double AmrPartBunch::getMass(int i) {
+    int l, g, dq;
+    std::tie(l,g,dq) = idxMap_m[i];
+    return m_particles[l][g][dq].m_data[4];
+}
+
+
 
 Vector_t AmrPartBunch::getP(int i) {
     int l, g, dq;
@@ -255,6 +262,12 @@ void AmrPartBunch::setQM(double q, int i) {
     m_particles[l][g][dq].m_data[0] = q;
 }
 
+void AmrPartBunch::setMass(double m, int i) {
+    int l, g, dq;
+    std::tie(l,g,dq) = idxMap_m[i];
+    
+    m_particles[l][g][dq].m_data[4] = m;
+}
 
 void AmrPartBunch::setP(Vector_t v, int i) {
     int l, g, dq;

ippl/test/AMR/AmrPartBunch.h

@@ -28,7 +28,7 @@
 
 /// Particle bunch class for BoxLib
 class AmrPartBunch : public PartBunchBase,
-                     public ParticleContainer<4 /*real attributes*/, 0>
+                     public ParticleContainer<5 /*real attributes*/, 0>
 {
 public:
     typedef std::map<int, std::tuple<int, int, int> > map_t;
@@ -89,6 +89,8 @@ public:
     
     inline double getQM(int i);
     
+    inline double getMass(int i);
+    
     inline Vector_t getP(int i);
     
     inline Vector_t getE(int i);
@@ -99,6 +101,8 @@ public:
     
     inline void setQM(double q, int i);
     
+    inline void setMass(double m, int i);
+    
     inline void setP(Vector_t v, int i);
     
     inline void setE(Vector_t Ef, int i);

ippl/test/AMR/CMakeLists.txt

@@ -89,6 +89,8 @@ IF (ENABLE_AMR_SOLVER)
     
     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)
     
+    add_executable (twostream twostream.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})
@@ -111,4 +113,6 @@ IF (ENABLE_AMR_SOLVER)
     
     target_link_libraries (testTracker ${LIBS} -lgfortran ${CCSE_LIBRARIES} -lcboxlib ${OTHER_CMAKE_EXE_LINKER_FLAGS} ${CMAKE_DL_LIBS})
     
+    target_link_libraries (twostream ${LIBS} -lgfortran ${CCSE_LIBRARIES} -lcboxlib ${OTHER_CMAKE_EXE_LINKER_FLAGS} ${CMAKE_DL_LIBS})
+    
 ENDIF(ENABLE_AMR_SOLVER)

ippl/test/AMR/Distribution.cpp

@@ -15,6 +15,8 @@ Distribution::Distribution():
     px_m(0),
     py_m(0),
     pz_m(0),
+    q_m(0),
+    mass_m(0),
     nloc_m(0),
     ntot_m(0)
 { }
@@ -45,6 +47,7 @@ void Distribution::uniform(double lower, double upper, size_t nloc, int seed) {
     pz_m.resize(nloc);
     
     q_m.resize(nloc);
+    mass_m.resize(nloc);
     
     for (size_t i = 0; i < nloc; ++i) {
         x_m[i] = dist(mt);
@@ -56,6 +59,7 @@ void Distribution::uniform(double lower, double upper, size_t nloc, int seed) {
         pz_m[i] = dist(mt);
         
         q_m[i] = 1.0;
+        mass_m[i] = 1.0;
     }
 }
 
@@ -84,6 +88,7 @@ void Distribution::gaussian(double mean, double stddev, size_t nloc, int seed) {
     pz_m.resize(nloc);
     
     q_m.resize(nloc);
+    mass_m.resize(nloc);
     
     for (size_t i = 0; i < nloc; ++i) {
         x_m[i] = dist(mt);
@@ -96,6 +101,73 @@ void Distribution::gaussian(double mean, double stddev, size_t nloc, int seed) {
         
         
         q_m[i] = 1.0;
+        mass_m[i] = 1.0;
+    }
+}
+
+void Distribution::twostream(const Vector_t& lower, const Vector_t& upper,
+                             const Vektor<std::size_t, 3>& nx, const Vektor<std::size_t, 3>& nv,
+                             const Vektor<double, 3>& vmax, double alpha)
+{
+    Vektor<double, 3> hx = (upper - lower) / Vector_t(nx);
+    Vektor<double, 3> hv = 2.0 * vmax / Vector_t(nv);
+    
+    double thr = 1.0e-12;
+    
+    double factor = 1.0 / ( M_PI * 30.0 );
+    
+    nloc_m = 0;
+    
+    if ( !Ippl::myNode() ) {
+        for (std::size_t i = 0; i < nx[0]; ++i) {
+            for (std::size_t j = 0; j < nx[1]; ++j) {
+                for ( std::size_t k = 0; k < nx[2]; ++k) {
+                    Vektor<double, 3> pos = Vektor<double,3>(
+                                                (0.5 + i) * hx[0] + lower[0],
+                                                (0.5 + j) * hx[1] + lower[1],
+                                                (0.5 + k) * hx[2] + lower[2]
+                                            );
+                    
+                    for (std::size_t iv = 0; iv < nv[0]; ++iv) {
+                        for (std::size_t jv = 0; jv < nv[1]; ++jv) {
+                            for (std::size_t kv = 0; kv < nv[2]; ++kv) {
+                                Vektor<double, 3> vel = -vmax + hv *
+                                                        Vektor<double, 3>(iv + 0.5,
+                                                                          jv + 0.5,
+                                                                          kv + 0.5);
+                                
+                                double v2 = vel[0] * vel[0] +
+                                            vel[1] * vel[1] +
+                                            vel[2] * vel[2];
+                                
+                                double f = factor * std::exp(-0.5 * v2) *
+                                           (1.0 + alpha * std::cos(0.5 * pos[2])) *
+                                           (1.0 + 0.5 * v2);
+                                
+                                double m = hx[0] * hv[0] *
+                                           hx[1] * hv[1] *
+                                           hx[2] * hv[2] * f;
+                                           
+                                
+                                if ( m > thr ) {
+                                    ++nloc_m;
+                                    x_m.push_back( pos[0] );
+                                    y_m.push_back( pos[1] );
+                                    z_m.push_back( pos[2] );
+                                    
+                                    px_m.push_back( vel[0] );
+                                    py_m.push_back( vel[1] );
+                                    pz_m.push_back( vel[2] );
+                                    
+                                    q_m.push_back( -m );
+                                    mass_m.push_back( m );
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
     }
 }
 
@@ -127,11 +199,13 @@ void Distribution::readH5(const std::string& filename, int step) {
     pz_m.resize(nloc_m);
     
     q_m.resize(nloc_m);
+    mass_m.resize(nloc_m);
     
     h5.read(x_m, px_m,
             y_m, py_m,
             z_m, pz_m,
             q_m,
+            mass_m,
             firstParticle,
             lastParticle);
     
@@ -155,6 +229,7 @@ void Distribution::injectBeam(PartBunchBase& bunch, bool doDelete, std::array<do
         bunch.setR(Vector_t(x_m[i] + shift[0], y_m[i] + shift[1], z_m[i] + shift[2]), i + prevnum);
         bunch.setP(Vector_t(px_m[i], py_m[i], pz_m[i]), i + prevnum);
         bunch.setQM(q_m[i], i + prevnum);
+        bunch.setMass(mass_m[i], i + prevnum);
         
         x_m.pop_back();
         y_m.pop_back();
@@ -164,6 +239,7 @@ void Distribution::injectBeam(PartBunchBase& bunch, bool doDelete, std::array<do
         pz_m.pop_back();
         
         q_m.pop_back();
+        mass_m.pop_back();
     }
 }
 

ippl/test/AMR/Distribution.h

@@ -49,6 +49,19 @@ public:
      */
     void gaussian(double mean, double stddev, size_t nloc, int seed);
     
+    /// Generate a two stream instability distribution according to B.\ Ulmer
+    /*!
+     * @param lower boundary of domain
+     * @param upper boundary of domain
+     * @param nx is the number of grid points in cooordinate space
+     * @param nv is the number of grid points in velocity space
+     * @param vmax is the max. velocity
+     * @param alpha is the amplitude of the initial disturbance
+     */
+    void twostream(const Vector_t& lower, const Vector_t& upper,
+                   const Vektor<std::size_t, 3>& nx, const Vektor<std::size_t, 3>& nv,
+                   const Vektor<double, 3>& vmax, double alpha = 0.5);
+    
     /// Read in a distribution from an H5 file
     /*!
      * @param filename is the path and name of the H5 file.
@@ -86,6 +99,7 @@ private:
     container_t py_m;   ///< Vertical particle momentum
     container_t pz_m;   ///< Longitudinal particle momentum
     container_t q_m;    ///< Particle charge (always set to 1.0, except for Distribution::readH5)
+    container_t mass_m; ///< Particl mass (always set to 1.0, except for Distribution::readH5 and Distribution::twostream)
     
     size_t nloc_m;      ///< Local number of particles
     size_t ntot_m;      ///< Total number of particles

ippl/test/AMR/PartBunchBase.h

@@ -94,6 +94,11 @@ public:
      */
     virtual double getQM(int i) = 0;
     
+    /// Access the mass of a particle
+    /*! @returns the mass of the i-th particle
+     */
+    virtual double getMass(int i) = 0;
+    
     /// Access the velocity of a particle
     /*! @returns the i-th particle velocity (px, py, pz)
      */
@@ -123,6 +128,13 @@ public:
      */
     virtual void setQM(double q, int i) = 0;
     
+    /// Set the particle charge-to-mass ratio
+    /*!
+     * @param m is the new mass
+     * @param i specifies the i-th particle
+     */
+    virtual void setMass(double q, int i) = 0;
+    
     /// Set the particle velocity
     /*!
      * @param v is the velocity (vx, vy, vz)

ippl/test/AMR/Solver.cpp

@@ -133,6 +133,16 @@ Solver::solve_with_f90(container_t& rhs,
             mg_bc[2*dir + 0] = MGT_BC_PER;
             mg_bc[2*dir + 1] = MGT_BC_PER;
         }
+    } else if ( Geometry::isAnyPeriodic() ) {
+        for (int dir = 0; dir < BL_SPACEDIM; ++dir) {
+            if ( Geometry::isPeriodic(dir) ) {
+                mg_bc[2*dir + 0] = MGT_BC_PER;
+                mg_bc[2*dir + 1] = MGT_BC_PER;
+            } else {
+                mg_bc[2*dir + 0] = MGT_BC_DIR;
+                mg_bc[2*dir + 1] = MGT_BC_DIR;
+            }
+        }
     } else {
 //         if ( ParallelDescriptor::IOProcessor() )
 //             std::cerr << "Dirichlet BC" << std::endl;

ippl/test/AMR/twostream.cpp

+#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)
+{
+    // =======================================================================                                                                                                                                   
+    // 4. prepare for multi-level solve                                                                                                                                                                          
+    // =======================================================================
+    
+    rhs.resize(nLevels);
+    phi.resize(nLevels);
+    grad_phi.resize(nLevels);
+    
+    for (int lev = 0; lev < nLevels; ++lev) {
+        initGridData(rhs, phi, grad_phi, myAmrOpal.boxArray()[lev], lev);
+    }
+
+    // Define the density on level 0 from all particles at all levels                                                                                                                                            
+    int base_level   = 0;
+    int finest_level = myAmrOpal.finestLevel();
+
+    dynamic_cast<AmrPartBunch*>(bunch)->AssignDensity(0, false, rhs, base_level, 1, finest_level);
+    
+    // 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,
+                        false);
+    
+    // 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 doTwoStream(Vektor<std::size_t, 3> nr,
+                 std::size_t nLevels,
+                 std::size_t maxBoxSize,
+                 double alpha,
+                 double dt,
+                 double epsilon,
+                 std::size_t nIter,
+                 Inform& msg)
+{
+    std::array<double, BL_SPACEDIM> lower = {{0.0, 0.0, 0.0}}; // m
+    std::array<double, BL_SPACEDIM> upper = {{4.0 * Physics::pi,
+                                              4.0 * Physics::pi,
+                                              4.0 * Physics::pi}
+    }; // 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);
+    
+    ParmParse pgeom("geometry");
+    Array<int> is_per = { 0, 0, 1};
+    pgeom.addarr("is_periodic", is_per);
+    
+    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 (std::size_t i = 0; i < nLevels; ++i)
+        rr[i] = 2;
+    
+    PartBunchBase* bunch = new AmrPartBunch(geom, dmap, ba, rr);
+    
+    
+    // initialize a particle distribution
+    Distribution dist;
+    dist.twostream(Vektor<double, 3>(lower[0], lower[1], lower[2]),
+                   Vektor<double, 3>(upper[0], upper[1], upper[2]),
+                   nr, Vektor<std::size_t, 3>(8, 8, 128),
+                   Vektor<double, 3>(6.0, 6.0, 6.0), alpha);
+    
+    // copy particles to the PartBunchBase object.
+    dist.injectBeam(*bunch);
+    
+    // redistribute on single-level
+    bunch->myUpdate();
+    
+    myAmrOpal.setBunch(dynamic_cast<AmrPartBunch*>(bunch));
+    
+    for (int i = 0; i <= myAmrOpal.finestLevel() && i < myAmrOpal.maxLevel(); ++i)
+        myAmrOpal.regrid(i /*lbase*/, 0.0 /*time*/);
+    
+    msg << "Multi-level statistics" << endl;
+    bunch->gatherStatistics();
+    
+    container_t rhs;
+    container_t phi;
+    container_t grad_phi;
+    
+    for (int i = 0; i < nIter; ++i) {
+        
+        for (std::size_t j = 0; j < bunch->getLocalNum(); ++j) {
+            int level = dynamic_cast<AmrPartBunch*>(bunch)->getLevel(j);
+            bunch->setR( bunch->getR(j) + dt * bunch->getP(j), j );
+        }
+            
+        for (int j = 0; j <= myAmrOpal.finestLevel() && j < myAmrOpal.maxLevel(); ++j)
+            myAmrOpal.regrid(j /*lbase*/, 0.0 /*time*/);
+        
+        doSolve(myAmrOpal, bunch, rhs, phi, grad_phi, geoms, rr, nLevels);
+        
+        for (std::size_t j = 0; j < bunch->getLocalNum(); ++j) {
+            int level = dynamic_cast<AmrPartBunch*>(bunch)->getLevel(j);
+            Vector_t externalE = dynamic_cast<AmrPartBunch*>(bunch)->interpolate(j, grad_phi[level]);
+            bunch->setP( bunch->getP(j) + dt * bunch->getQM(j) / bunch->getMass(j) * externalE, j);
+        }
+        
+    }
+    
+    
+    
+}
+
+
+int main(int argc, char *argv[]) {
+    
+    Ippl ippl(argc, argv);
+    BoxLib::Initialize(argc,argv, false);
+    
+    Inform msg("TwoStream");
+    
+
+    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] [#levels]"
+         << " [max. box size] [alpha] [timstep] [epsilon] [#iterations]"
+         << " [out: timing file name (optiona)]";
+    
+    if ( argc < 9 ) {
+        msg << call.str() << endl;
+        return -1;
+    }
+    
+    // number of grid points in each direction
+    Vektor<std::size_t, 3> nr(std::atoi(argv[1]),
+                              std::atoi(argv[2]),
+                              std::atoi(argv[3]));
+    
+    
+    std::size_t nLevels    = std::atoi( argv[4] ) + 1; // i.e. nLevels = 0 --> only single level
+    std::size_t maxBoxSize = std::atoi( argv[5] );
+    double alpha           = std::atof( argv[6] );
+    double dt              = std::atof( argv[7] );
+    double epsilon         = std::atof( argv[8] );
+    std::size_t nIter      = std::atof( argv[9] );
+    
+    msg << "Particle test running with" << endl
+        << "- #level     = " << nLevels << endl
+        << "- grid       = " << nr << endl
+        << "- max. size  = " << maxBoxSize << endl
+        << "- amplitude  = " << alpha << endl
+        << "- epsilon    = " << epsilon << endl
+        << "- time step  = " << dt << endl
+        << "- #steps     = " << nIter << endl;
+    
+    doTwoStream(nr, nLevels, maxBoxSize,
+                alpha, dt, epsilon, nIter, 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";
+    
+    if ( argc == 8 ) {
+        timefile.str(std::string());
+        timefile << std::string(argv[7]);
+    }
+    
+    IpplTimings::print(timefile.str());
+    
+    return 0;
+}
\ No newline at end of file