diff --git a/.clang-format b/.clang-format
new file mode 100644
index 000000000..46bb4bb3b
--- /dev/null
+++ b/.clang-format
@@ -0,0 +1,24 @@
+Language: Cpp
+BasedOnStyle: Google
+AccessModifierOffset: -4
+AlignAfterOpenBracket: Align
+AlignConsecutiveAssignments: true
+AlignConsecutiveMacros: AcrossEmptyLinesAndComments
+AlignEscapedNewlines: Left
+AllowShortFunctionsOnASingleLine: Inline
+AllowShortIfStatementsOnASingleLine: Never
+AllowShortLambdasOnASingleLine: Empty
+AllowShortLoopsOnASingleLine: false
+AllowShortEnumsOnASingleLine: false
+AttributeMacros: ['KOKKOS_INLINE_FUNCTION']
+BreakBeforeBinaryOperators: NonAssignment
+BreakConstructorInitializers: BeforeComma
+ColumnLimit: 100
+DerivePointerAlignment: false
+IndentWidth: 4
+IncludeBlocks: Preserve
+IndentGotoLabels: false
+NamespaceIndentation: All
+PackConstructorInitializers: Never
+ReflowComments: true
+StatementMacros: ['DefineUnaryOperation', 'DefineBinaryOperation']
diff --git a/alpine/BumponTailInstability.cpp b/alpine/BumponTailInstability.cpp
index c15cf60aa..f27143eb7 100644
--- a/alpine/BumponTailInstability.cpp
+++ b/alpine/BumponTailInstability.cpp
@@ -9,13 +9,13 @@
// stype = Field solver type e.g., FFT
// lbthres = Load balancing threshold i.e., lbthres*100 is the maximum load imbalance
// percentage which can be tolerated and beyond which
-// particle load balancing occurs. A value of 0.01 is good for many typical
+// particle load balancing occurs. A value of 0.01 is good for many typical
// simulations.
// ovfactor = Over-allocation factor for the buffers used in the communication. Typical
// values are 1.0, 2.0. Value 1.0 means no over-allocation.
// Example:
// srun ./BumponTailInstability 128 128 128 10000 10 FFT 0.01 2.0 --info 10
-// Change the TestName to TwoStreamInstability or BumponTailInstability
+// Change the TestName to TwoStreamInstability or BumponTailInstability
// in order to simulate the Two stream instability or bump on tail instability
// cases
//
@@ -34,176 +34,165 @@
// along with IPPL. If not, see .
//
-#include "ChargedParticles.hpp"
-#include
-#include
-#include
+#include
#include
+#include
#include
-#include
+#include
+#include
+#include "ChargedParticles.hpp"
-#include
+#include
#include
#include "Utility/IpplTimings.h"
template
struct Newton1D {
+ double tol = 1e-12;
+ int max_iter = 20;
+ double pi = std::acos(-1.0);
- double tol = 1e-12;
- int max_iter = 20;
- double pi = std::acos(-1.0);
-
- T k, delta, u;
-
- KOKKOS_INLINE_FUNCTION
- Newton1D() {}
-
- KOKKOS_INLINE_FUNCTION
- Newton1D(const T& k_, const T& delta_,
- const T& u_)
- : k(k_), delta(delta_), u(u_) {}
-
- KOKKOS_INLINE_FUNCTION
- ~Newton1D() {}
-
- KOKKOS_INLINE_FUNCTION
- T f(T& x) {
- T F;
- F = x + (delta * (std::sin(k * x) / k)) - u;
- return F;
- }
-
- KOKKOS_INLINE_FUNCTION
- T fprime(T& x) {
- T Fprime;
- Fprime = 1 + (delta * std::cos(k * x));
- return Fprime;
- }
-
- KOKKOS_FUNCTION
- void solve(T& x) {
- int iterations = 0;
- while (iterations < max_iter && std::fabs(f(x)) > tol) {
- x = x - (f(x)/fprime(x));
- iterations += 1;
- }
- }
-};
+ T k, delta, u;
+
+ KOKKOS_INLINE_FUNCTION Newton1D() {}
+
+ KOKKOS_INLINE_FUNCTION Newton1D(const T& k_, const T& delta_, const T& u_)
+ : k(k_)
+ , delta(delta_)
+ , u(u_) {}
+
+ KOKKOS_INLINE_FUNCTION ~Newton1D() {}
+
+ KOKKOS_INLINE_FUNCTION T f(T& x) {
+ T F;
+ F = x + (delta * (std::sin(k * x) / k)) - u;
+ return F;
+ }
+
+ KOKKOS_INLINE_FUNCTION T fprime(T& x) {
+ T Fprime;
+ Fprime = 1 + (delta * std::cos(k * x));
+ return Fprime;
+ }
+ KOKKOS_FUNCTION
+ void solve(T& x) {
+ int iterations = 0;
+ while (iterations < max_iter && std::fabs(f(x)) > tol) {
+ x = x - (f(x) / fprime(x));
+ iterations += 1;
+ }
+ }
+};
template
struct generate_random {
+ using view_type = typename ippl::detail::ViewType::view_type;
+ using value_type = typename T::value_type;
+ // Output View for the random numbers
+ view_type x, v;
+
+ // The GeneratorPool
+ GeneratorPool rand_pool;
+
+ value_type delta, sigma, muBulk, muBeam;
+ size_type nlocBulk;
+
+ T k, minU, maxU;
+
+ // Initialize all members
+ generate_random(view_type x_, view_type v_, GeneratorPool rand_pool_, value_type& delta_, T& k_,
+ value_type& sigma_, value_type& muBulk_, value_type& muBeam_,
+ size_type& nlocBulk_, T& minU_, T& maxU_)
+ : x(x_)
+ , v(v_)
+ , rand_pool(rand_pool_)
+ , delta(delta_)
+ , sigma(sigma_)
+ , muBulk(muBulk_)
+ , muBeam(muBeam_)
+ , nlocBulk(nlocBulk_)
+ , k(k_)
+ , minU(minU_)
+ , maxU(maxU_) {}
+
+ KOKKOS_INLINE_FUNCTION void operator()(const size_t i) const {
+ // Get a random number state from the pool for the active thread
+ typename GeneratorPool::generator_type rand_gen = rand_pool.get_state();
+
+ bool isBeam = (i >= nlocBulk);
+
+ value_type muZ = (value_type)(((!isBeam) * muBulk) + (isBeam * muBeam));
+
+ for (unsigned d = 0; d < Dim - 1; ++d) {
+ x(i)[d] = rand_gen.drand(minU[d], maxU[d]);
+ v(i)[d] = rand_gen.normal(0.0, sigma);
+ }
+ v(i)[Dim - 1] = rand_gen.normal(muZ, sigma);
+
+ value_type u = rand_gen.drand(minU[Dim - 1], maxU[Dim - 1]);
+ x(i)[Dim - 1] = u / (1 + delta);
+ Newton1D solver(k[Dim - 1], delta, u);
+ solver.solve(x(i)[Dim - 1]);
- using view_type = typename ippl::detail::ViewType::view_type;
- using value_type = typename T::value_type;
- // Output View for the random numbers
- view_type x, v;
-
- // The GeneratorPool
- GeneratorPool rand_pool;
-
- value_type delta, sigma, muBulk, muBeam;
- size_type nlocBulk;
-
- T k, minU, maxU;
-
- // Initialize all members
- generate_random(view_type x_, view_type v_, GeneratorPool rand_pool_,
- value_type& delta_, T& k_, value_type& sigma_,
- value_type& muBulk_, value_type& muBeam_,
- size_type& nlocBulk_, T& minU_, T& maxU_)
- : x(x_), v(v_), rand_pool(rand_pool_),
- delta(delta_), sigma(sigma_), muBulk(muBulk_), muBeam(muBeam_),
- nlocBulk(nlocBulk_), k(k_), minU(minU_), maxU(maxU_) {}
-
- KOKKOS_INLINE_FUNCTION
- void operator()(const size_t i) const {
- // Get a random number state from the pool for the active thread
- typename GeneratorPool::generator_type rand_gen = rand_pool.get_state();
-
- bool isBeam = (i >= nlocBulk);
-
- value_type muZ = (value_type)(((!isBeam) * muBulk) + (isBeam * muBeam));
-
- for (unsigned d = 0; d < Dim-1; ++d) {
-
- x(i)[d] = rand_gen.drand(minU[d], maxU[d]);
- v(i)[d] = rand_gen.normal(0.0, sigma);
+ // Give the state back, which will allow another thread to acquire it
+ rand_pool.free_state(rand_gen);
}
- v(i)[Dim-1] = rand_gen.normal(muZ, sigma);
-
- value_type u = rand_gen.drand(minU[Dim-1], maxU[Dim-1]);
- x(i)[Dim-1] = u / (1 + delta);
- Newton1D solver(k[Dim-1], delta, u);
- solver.solve(x(i)[Dim-1]);
-
-
- // Give the state back, which will allow another thread to acquire it
- rand_pool.free_state(rand_gen);
- }
};
-double CDF(const double& x, const double& delta, const double& k,
- const unsigned& dim) {
-
- bool isDimZ = (dim == (Dim-1));
- double cdf = x + (double)(isDimZ * ((delta / k) * std::sin(k * x)));
- return cdf;
+double CDF(const double& x, const double& delta, const double& k, const unsigned& dim) {
+ bool isDimZ = (dim == (Dim - 1));
+ double cdf = x + (double)(isDimZ * ((delta / k) * std::sin(k * x)));
+ return cdf;
}
KOKKOS_FUNCTION
-double PDF(const Vector_t& xvec, const double& delta,
- const Vector_t& kw) {
-
- double pdf = 1.0 * 1.0 * (1.0 + delta * std::cos(kw[Dim-1] * xvec[Dim-1]));
+double PDF(const Vector_t& xvec, const double& delta, const Vector_t& kw) {
+ double pdf = 1.0 * 1.0 * (1.0 + delta * std::cos(kw[Dim - 1] * xvec[Dim - 1]));
return pdf;
}
-//const char* TestName = "BumponTailInstability";
+// const char* TestName = "BumponTailInstability";
const char* TestName = "TwoStreamInstability";
-int main(int argc, char *argv[]){
+int main(int argc, char* argv[]) {
Ippl ippl(argc, argv);
Inform msg(TestName);
- Inform msg2all(TestName,INFORM_ALL_NODES);
+ Inform msg2all(TestName, INFORM_ALL_NODES);
Ippl::Comm->setDefaultOverallocation(std::atof(argv[8]));
- auto start = std::chrono::high_resolution_clock::now();
- ippl::Vector nr = {
+ auto start = std::chrono::high_resolution_clock::now();
+ ippl::Vector nr = {
std::atoi(argv[1]),
std::atoi(argv[2]),
- std::atoi(argv[3])
+ std::atoi(argv[3]),
};
- static IpplTimings::TimerRef mainTimer = IpplTimings::getTimer("total");
- static IpplTimings::TimerRef particleCreation = IpplTimings::getTimer("particlesCreation");
- static IpplTimings::TimerRef dumpDataTimer = IpplTimings::getTimer("dumpData");
- static IpplTimings::TimerRef PTimer = IpplTimings::getTimer("pushVelocity");
- static IpplTimings::TimerRef RTimer = IpplTimings::getTimer("pushPosition");
- static IpplTimings::TimerRef updateTimer = IpplTimings::getTimer("update");
- static IpplTimings::TimerRef DummySolveTimer = IpplTimings::getTimer("solveWarmup");
- static IpplTimings::TimerRef SolveTimer = IpplTimings::getTimer("solve");
+ static IpplTimings::TimerRef mainTimer = IpplTimings::getTimer("total");
+ static IpplTimings::TimerRef particleCreation = IpplTimings::getTimer("particlesCreation");
+ static IpplTimings::TimerRef dumpDataTimer = IpplTimings::getTimer("dumpData");
+ static IpplTimings::TimerRef PTimer = IpplTimings::getTimer("pushVelocity");
+ static IpplTimings::TimerRef RTimer = IpplTimings::getTimer("pushPosition");
+ static IpplTimings::TimerRef updateTimer = IpplTimings::getTimer("update");
+ static IpplTimings::TimerRef DummySolveTimer = IpplTimings::getTimer("solveWarmup");
+ static IpplTimings::TimerRef SolveTimer = IpplTimings::getTimer("solve");
static IpplTimings::TimerRef domainDecomposition = IpplTimings::getTimer("loadBalance");
IpplTimings::startTimer(mainTimer);
const size_type totalP = std::atoll(argv[4]);
- const unsigned int nt = std::atoi(argv[5]);
+ const unsigned int nt = std::atoi(argv[5]);
- msg << TestName
- << endl
- << "nt " << nt << " Np= "
- << totalP << " grid = " << nr
- << endl;
+ msg << TestName << endl << "nt " << nt << " Np= " << totalP << " grid = " << nr << endl;
using bunch_type = ChargedParticles;
- std::unique_ptr P;
+ std::unique_ptr P;
ippl::NDIndex domain;
- for (unsigned i = 0; i< Dim; i++) {
+ for (unsigned i = 0; i < Dim; i++) {
domain[i] = ippl::Index(nr[i]);
}
@@ -214,54 +203,51 @@ int main(int argc, char *argv[]){
Vector_t kw;
double sigma, muBulk, muBeam, epsilon, delta;
-
-
- if(std::strcmp(TestName,"TwoStreamInstability") == 0) {
- // Parameters for two stream instability as in
+
+ if (std::strcmp(TestName, "TwoStreamInstability") == 0) {
+ // Parameters for two stream instability as in
// https://www.frontiersin.org/articles/10.3389/fphy.2018.00105/full
- kw = {0.5, 0.5, 0.5};
- sigma = 0.1;
+ kw = {0.5, 0.5, 0.5};
+ sigma = 0.1;
epsilon = 0.5;
- muBulk = -pi / 2.0;
- muBeam = pi / 2.0;
- delta = 0.01;
- }
- else if(std::strcmp(TestName,"BumponTailInstability") == 0) {
- kw = {0.21, 0.21, 0.21};
- sigma = 1.0 / std::sqrt(2.0);
+ muBulk = -pi / 2.0;
+ muBeam = pi / 2.0;
+ delta = 0.01;
+ } else if (std::strcmp(TestName, "BumponTailInstability") == 0) {
+ kw = {0.21, 0.21, 0.21};
+ sigma = 1.0 / std::sqrt(2.0);
epsilon = 0.1;
- muBulk = 0.0;
- muBeam = 4.0;
- delta = 0.01;
- }
- else {
- //Default value is two stream instability
- kw = {0.5, 0.5, 0.5};
- sigma = 0.1;
+ muBulk = 0.0;
+ muBeam = 4.0;
+ delta = 0.01;
+ } else {
+ // Default value is two stream instability
+ kw = {0.5, 0.5, 0.5};
+ sigma = 0.1;
epsilon = 0.5;
- muBulk = -pi / 2.0;
- muBeam = pi / 2.0;
- delta = 0.01;
+ muBulk = -pi / 2.0;
+ muBeam = pi / 2.0;
+ delta = 0.01;
}
Vector_t rmin(0.0);
- Vector_t rmax = 2 * pi / kw ;
- double dx = rmax[0] / nr[0];
- double dy = rmax[1] / nr[1];
- double dz = rmax[2] / nr[2];
+ Vector_t rmax = 2 * pi / kw;
+ double dx = rmax[0] / nr[0];
+ double dy = rmax[1] / nr[1];
+ double dz = rmax[2] / nr[2];
- Vector_t hr = {dx, dy, dz};
+ Vector_t hr = {dx, dy, dz};
Vector_t origin = {rmin[0], rmin[1], rmin[2]};
- const double dt = 0.5*dx;//0.05
+ const double dt = 0.5 * dx; // 0.05
- const bool isAllPeriodic=true;
+ const bool isAllPeriodic = true;
Mesh_t mesh(domain, hr, origin);
FieldLayout_t FL(domain, decomp, isAllPeriodic);
PLayout_t PL(FL, mesh);
- //Q = -\int\int f dx dv
+ // Q = -\int\int f dx dv
double Q = -rmax[0] * rmax[1] * rmax[2];
- P = std::make_unique(PL,hr,rmin,rmax,decomp,Q);
+ P = std::make_unique(PL, hr, rmin, rmax, decomp, Q);
P->nr_m = nr;
@@ -272,7 +258,7 @@ int main(int argc, char *argv[]){
P->stype_m = argv[6];
P->initSolver();
- P->time_m = 0.0;
+ P->time_m = 0.0;
P->loadbalancethreshold_m = std::atof(argv[7]);
bool isFirstRepartition;
@@ -280,99 +266,91 @@ int main(int argc, char *argv[]){
if ((P->loadbalancethreshold_m != 1.0) && (Ippl::Comm->size() > 1)) {
msg << "Starting first repartition" << endl;
IpplTimings::startTimer(domainDecomposition);
- isFirstRepartition = true;
+ isFirstRepartition = true;
const ippl::NDIndex& lDom = FL.getLocalNDIndex();
- const int nghost = P->rho_m.getNghost();
- using mdrange_type = Kokkos::MDRangePolicy>;
- auto rhoview = P->rho_m.getView();
-
- Kokkos::parallel_for("Assign initial rho based on PDF",
- mdrange_type({nghost, nghost, nghost},
- {rhoview.extent(0) - nghost,
- rhoview.extent(1) - nghost,
- rhoview.extent(2) - nghost}),
- KOKKOS_LAMBDA(const int i,
- const int j,
- const int k)
- {
- //local to global index conversion
- const size_t ig = i + lDom[0].first() - nghost;
- const size_t jg = j + lDom[1].first() - nghost;
- const size_t kg = k + lDom[2].first() - nghost;
- double x = (ig + 0.5) * hr[0] + origin[0];
- double y = (jg + 0.5) * hr[1] + origin[1];
- double z = (kg + 0.5) * hr[2] + origin[2];
-
- Vector_t xvec = {x, y, z};
-
- rhoview(i, j, k) = PDF(xvec, delta, kw);
-
- });
+ const int nghost = P->rho_m.getNghost();
+ using mdrange_type = Kokkos::MDRangePolicy>;
+ auto rhoview = P->rho_m.getView();
+
+ Kokkos::parallel_for(
+ "Assign initial rho based on PDF",
+ mdrange_type({nghost, nghost, nghost},
+ {rhoview.extent(0) - nghost, rhoview.extent(1) - nghost,
+ rhoview.extent(2) - nghost}),
+ KOKKOS_LAMBDA(const int i, const int j, const int k) {
+ // local to global index conversion
+ const size_t ig = i + lDom[0].first() - nghost;
+ const size_t jg = j + lDom[1].first() - nghost;
+ const size_t kg = k + lDom[2].first() - nghost;
+ double x = (ig + 0.5) * hr[0] + origin[0];
+ double y = (jg + 0.5) * hr[1] + origin[1];
+ double z = (kg + 0.5) * hr[2] + origin[2];
+
+ Vector_t xvec = {x, y, z};
+
+ rhoview(i, j, k) = PDF(xvec, delta, kw);
+ });
Kokkos::fence();
-
+
P->initializeORB(FL, mesh);
P->repartition(FL, mesh, bunchBuffer, isFirstRepartition);
IpplTimings::stopTimer(domainDecomposition);
}
-
+
msg << "First domain decomposition done" << endl;
IpplTimings::startTimer(particleCreation);
-
typedef ippl::detail::RegionLayout RegionLayout_t;
- const RegionLayout_t& RLayout = PL.getRegionLayout();
+ const RegionLayout_t& RLayout = PL.getRegionLayout();
const typename RegionLayout_t::host_mirror_type Regions = RLayout.gethLocalRegions();
Vector_t Nr, Dr, minU, maxU;
int myRank = Ippl::Comm->rank();
- for (unsigned d = 0; d rank() < rest )
+ if (Ippl::Comm->rank() < rest)
++nloc;
P->create(nloc);
- Kokkos::Random_XorShift64_Pool<> rand_pool64((size_type)(42 + 100*Ippl::Comm->rank()));
+ Kokkos::Random_XorShift64_Pool<> rand_pool64((size_type)(42 + 100 * Ippl::Comm->rank()));
- Kokkos::parallel_for(nloc,
- generate_random, Dim>(
- P->R.getView(), P->P.getView(), rand_pool64, delta, kw, sigma, muBulk,
- muBeam, nlocBulk, minU, maxU));
+ Kokkos::parallel_for(nloc, generate_random, Dim>(
+ P->R.getView(), P->P.getView(), rand_pool64, delta, kw, sigma,
+ muBulk, muBeam, nlocBulk, minU, maxU));
Kokkos::fence();
Ippl::Comm->barrier();
- IpplTimings::stopTimer(particleCreation);
-
- P->q = P->Q_m/totalP;
- //P->dumpParticleData();
+ IpplTimings::stopTimer(particleCreation);
+
+ P->q = P->Q_m / totalP;
+ // P->dumpParticleData();
msg << "particles created and initial conditions assigned " << endl;
isFirstRepartition = false;
- //The update after the particle creation is not needed as the
- //particles are generated locally
+ // The update after the particle creation is not needed as the
+ // particles are generated locally
IpplTimings::startTimer(DummySolveTimer);
P->rho_m = 0.0;
P->solver_mp->solve();
IpplTimings::stopTimer(DummySolveTimer);
-
-
+
P->scatterCIC(totalP, 0, hr);
IpplTimings::startTimer(SolveTimer);
@@ -384,13 +362,12 @@ int main(int argc, char *argv[]){
IpplTimings::startTimer(dumpDataTimer);
P->dumpBumponTail();
P->gatherStatistics(totalP);
- //P->dumpLocalDomains(FL, 0);
+ // P->dumpLocalDomains(FL, 0);
IpplTimings::stopTimer(dumpDataTimer);
// begin main timestep loop
msg << "Starting iterations ..." << endl;
- for (unsigned int it=0; itP = P->P - 0.5 * dt * P->E;
IpplTimings::stopTimer(PTimer);
- //drift
+ // drift
IpplTimings::startTimer(RTimer);
P->R = P->R + dt * P->P;
IpplTimings::stopTimer(RTimer);
- //Since the particles have moved spatially update them to correct processors
- IpplTimings::startTimer(updateTimer);
+ // Since the particles have moved spatially update them to correct processors
+ IpplTimings::startTimer(updateTimer);
PL.update(*P, bunchBuffer);
IpplTimings::stopTimer(updateTimer);
// Domain Decomposition
- if (P->balance(totalP, it+1)) {
- msg << "Starting repartition" << endl;
- IpplTimings::startTimer(domainDecomposition);
- P->repartition(FL, mesh, bunchBuffer, isFirstRepartition);
- IpplTimings::stopTimer(domainDecomposition);
- //IpplTimings::startTimer(dumpDataTimer);
- //P->dumpLocalDomains(FL, it+1);
- //IpplTimings::stopTimer(dumpDataTimer);
+ if (P->balance(totalP, it + 1)) {
+ msg << "Starting repartition" << endl;
+ IpplTimings::startTimer(domainDecomposition);
+ P->repartition(FL, mesh, bunchBuffer, isFirstRepartition);
+ IpplTimings::stopTimer(domainDecomposition);
+ // IpplTimings::startTimer(dumpDataTimer);
+ // P->dumpLocalDomains(FL, it+1);
+ // IpplTimings::stopTimer(dumpDataTimer);
}
+ // scatter the charge onto the underlying grid
+ P->scatterCIC(totalP, it + 1, hr);
- //scatter the charge onto the underlying grid
- P->scatterCIC(totalP, it+1, hr);
-
- //Field solve
+ // Field solve
IpplTimings::startTimer(SolveTimer);
P->solver_mp->solve();
IpplTimings::stopTimer(SolveTimer);
@@ -434,7 +410,7 @@ int main(int argc, char *argv[]){
// gather E field
P->gatherCIC();
- //kick
+ // kick
IpplTimings::startTimer(PTimer);
P->P = P->P - 0.5 * dt * P->E;
IpplTimings::stopTimer(PTimer);
@@ -444,7 +420,7 @@ int main(int argc, char *argv[]){
P->dumpBumponTail();
P->gatherStatistics(totalP);
IpplTimings::stopTimer(dumpDataTimer);
- msg << "Finished time step: " << it+1 << " time: " << P->time_m << endl;
+ msg << "Finished time step: " << it + 1 << " time: " << P->time_m << endl;
}
msg << TestName << ": End." << endl;
@@ -453,7 +429,8 @@ int main(int argc, char *argv[]){
IpplTimings::print(std::string("timing.dat"));
auto end = std::chrono::high_resolution_clock::now();
- std::chrono::duration time_chrono = std::chrono::duration_cast>(end - start);
+ std::chrono::duration time_chrono =
+ std::chrono::duration_cast>(end - start);
std::cout << "Elapsed time: " << time_chrono.count() << std::endl;
return 0;
diff --git a/alpine/ChargedParticles.hpp b/alpine/ChargedParticles.hpp
index e64417e19..9885934db 100644
--- a/alpine/ChargedParticles.hpp
+++ b/alpine/ChargedParticles.hpp
@@ -23,24 +23,24 @@
constexpr unsigned Dim = 3;
// some typedefs
-typedef ippl::ParticleSpatialLayout PLayout_t;
-typedef ippl::UniformCartesian Mesh_t;
+typedef ippl::ParticleSpatialLayout PLayout_t;
+typedef ippl::UniformCartesian Mesh_t;
typedef ippl::FieldLayout FieldLayout_t;
typedef ippl::OrthogonalRecursiveBisection ORB;
using size_type = ippl::detail::size_type;
-template
+template
using Vector = ippl::Vector;
-template
+template
using Field = ippl::Field;
-template
+template
using ParticleAttrib = ippl::ParticleAttrib;
-typedef Vector Vector_t;
-typedef Field Field_t;
+typedef Vector Vector_t;
+typedef Field Field_t;
typedef Field VField_t;
typedef ippl::FFTPeriodicPoissonSolver Solver_t;
@@ -49,10 +49,7 @@ const double pi = std::acos(-1.0);
// Test programs have to define this variable for VTK dump purposes
extern const char* TestName;
-void dumpVTK(VField_t& E, int nx, int ny, int nz, int iteration,
- double dx, double dy, double dz) {
-
-
+void dumpVTK(VField_t& E, int nx, int ny, int nz, int iteration, double dx, double dy, double dz) {
typename VField_t::view_type::host_mirror_type host_view = E.getHostMirror();
std::stringstream fname;
@@ -71,27 +68,23 @@ void dumpVTK(VField_t& E, int nx, int ny, int nz, int iteration,
vtkout << TestName << endl;
vtkout << "ASCII" << endl;
vtkout << "DATASET STRUCTURED_POINTS" << endl;
- vtkout << "DIMENSIONS " << nx+3 << " " << ny+3 << " " << nz+3 << endl;
- vtkout << "ORIGIN " << -dx << " " << -dy << " " << -dz << endl;
+ vtkout << "DIMENSIONS " << nx + 3 << " " << ny + 3 << " " << nz + 3 << endl;
+ vtkout << "ORIGIN " << -dx << " " << -dy << " " << -dz << endl;
vtkout << "SPACING " << dx << " " << dy << " " << dz << endl;
- vtkout << "CELL_DATA " << (nx+2)*(ny+2)*(nz+2) << endl;
+ vtkout << "CELL_DATA " << (nx + 2) * (ny + 2) * (nz + 2) << endl;
vtkout << "VECTORS E-Field float" << endl;
- for (int z=0; z
+template
class ChargedParticles : public ippl::ParticleBase {
public:
VField_t E_m;
@@ -155,55 +147,46 @@ class ChargedParticles : public ippl::ParticleBase {
double rhoNorm_m;
unsigned int loadbalancefreq_m;
-
- double loadbalancethreshold_m;
+ double loadbalancethreshold_m;
public:
- ParticleAttrib q; // charge
+ ParticleAttrib q; // charge
typename ippl::ParticleBase::particle_position_type P; // particle velocity
- typename ippl::ParticleBase::particle_position_type E; // electric field at particle position
-
+ typename ippl::ParticleBase::particle_position_type
+ E; // electric field at particle position
/*
This constructor is mandatory for all derived classes from
ParticleBase as the bunch buffer uses this
*/
ChargedParticles(PLayout& pl)
- : ippl::ParticleBase(pl)
- {
+ : ippl::ParticleBase(pl) {
// register the particle attributes
this->addAttribute(q);
this->addAttribute(P);
this->addAttribute(E);
}
- ChargedParticles(PLayout& pl,
- Vector_t hr,
- Vector_t rmin,
- Vector_t rmax,
- ippl::e_dim_tag decomp[Dim],
- double Q)
- : ippl::ParticleBase(pl)
- , hr_m(hr)
- , rmin_m(rmin)
- , rmax_m(rmax)
- , Q_m(Q)
- {
+ ChargedParticles(PLayout& pl, Vector_t hr, Vector_t rmin, Vector_t rmax,
+ ippl::e_dim_tag decomp[Dim], double Q)
+ : ippl::ParticleBase(pl)
+ , hr_m(hr)
+ , rmin_m(rmin)
+ , rmax_m(rmax)
+ , Q_m(Q) {
// register the particle attributes
this->addAttribute(q);
this->addAttribute(P);
this->addAttribute(E);
setupBCs();
for (unsigned int i = 0; i < Dim; i++)
- decomp_m[i]=decomp[i];
+ decomp_m[i] = decomp[i];
}
- ~ChargedParticles(){ }
+ ~ChargedParticles() {}
- void setupBCs() {
- setBCAllPeriodic();
- }
+ void setupBCs() { setBCAllPeriodic(); }
void updateLayout(FieldLayout_t& fl, Mesh_t& mesh, ChargedParticles& buffer,
bool& isFirstRepartition) {
@@ -219,39 +202,36 @@ class ChargedParticles : public ippl::ParticleBase {
IpplTimings::stopTimer(tupdateLayout);
static IpplTimings::TimerRef tupdatePLayout = IpplTimings::getTimer("updatePB");
IpplTimings::startTimer(tupdatePLayout);
- if(!isFirstRepartition) {
+ if (!isFirstRepartition) {
layout.update(*this, buffer);
}
IpplTimings::stopTimer(tupdatePLayout);
}
- void initializeORB(FieldLayout_t& fl, Mesh_t& mesh) {
- orb.initialize(fl, mesh, rho_m);
- }
+ void initializeORB(FieldLayout_t& fl, Mesh_t& mesh) { orb.initialize(fl, mesh, rho_m); }
- void repartition(FieldLayout_t& fl, Mesh_t& mesh, ChargedParticles& buffer,
+ void repartition(FieldLayout_t& fl, Mesh_t& mesh, ChargedParticles& buffer,
bool& isFirstRepartition) {
// Repartition the domains
bool res = orb.binaryRepartition(this->R, fl, isFirstRepartition);
if (res != true) {
- std::cout << "Could not repartition!" << std::endl;
- return;
+ std::cout << "Could not repartition!" << std::endl;
+ return;
}
// Update
this->updateLayout(fl, mesh, buffer, isFirstRepartition);
this->solver_mp->setRhs(rho_m);
}
- bool balance(size_type totalP, const unsigned int nstep){
- if(std::strcmp(TestName,"UniformPlasmaTest") == 0) {
+ bool balance(size_type totalP, const unsigned int nstep) {
+ if (std::strcmp(TestName, "UniformPlasmaTest") == 0) {
return (nstep % loadbalancefreq_m == 0);
- }
- else {
+ } else {
int local = 0;
std::vector res(Ippl::Comm->size());
- double equalPart = (double) totalP / Ippl::Comm->size();
- double dev = std::abs((double)this->getLocalNum() - equalPart) / totalP;
+ double equalPart = (double)totalP / Ippl::Comm->size();
+ double dev = std::abs((double)this->getLocalNum() - equalPart) / totalP;
if (dev > loadbalancethreshold_m)
local = 1;
MPI_Allgather(&local, 1, MPI_INT, res.data(), 1, MPI_INT, Ippl::getComm());
@@ -266,12 +246,10 @@ class ChargedParticles : public ippl::ParticleBase {
void gatherStatistics(size_type totalP) {
std::vector imb(Ippl::Comm->size());
- double equalPart = (double) totalP / Ippl::Comm->size();
- double dev = (std::abs((double)this->getLocalNum() - equalPart)
- / totalP) * 100.0;
- MPI_Gather(&dev, 1, MPI_DOUBLE, imb.data(), 1, MPI_DOUBLE, 0,
- Ippl::getComm());
-
+ double equalPart = (double)totalP / Ippl::Comm->size();
+ double dev = (std::abs((double)this->getLocalNum() - equalPart) / totalP) * 100.0;
+ MPI_Gather(&dev, 1, MPI_DOUBLE, imb.data(), 1, MPI_DOUBLE, 0, Ippl::getComm());
+
if (Ippl::Comm->rank() == 0) {
std::stringstream fname;
fname << "data/LoadBalance_";
@@ -282,90 +260,79 @@ class ChargedParticles : public ippl::ParticleBase {
csvout.precision(5);
csvout.setf(std::ios::scientific, std::ios::floatfield);
- if(time_m == 0.0) {
+ if (time_m == 0.0) {
csvout << "time, rank, imbalance percentage" << endl;
}
- for(int r=0; r < Ippl::Comm->size(); ++r) {
- csvout << time_m << " "
- << r << " "
- << imb[r] << endl;
+ for (int r = 0; r < Ippl::Comm->size(); ++r) {
+ csvout << time_m << " " << r << " " << imb[r] << endl;
}
}
Ippl::Comm->barrier();
-
}
- void gatherCIC() {
-
- gather(this->E, E_m, this->R);
-
- }
+ void gatherCIC() { gather(this->E, E_m, this->R); }
void scatterCIC(size_type totalP, unsigned int iteration, Vector_t& hrField) {
+ Inform m("scatter ");
+ rho_m = 0.0;
+ scatter(q, rho_m, this->R);
- Inform m("scatter ");
-
- rho_m = 0.0;
- scatter(q, rho_m, this->R);
-
- static IpplTimings::TimerRef sumTimer = IpplTimings::getTimer("Check");
- IpplTimings::startTimer(sumTimer);
- double Q_grid = rho_m.sum();
+ static IpplTimings::TimerRef sumTimer = IpplTimings::getTimer("Check");
+ IpplTimings::startTimer(sumTimer);
+ double Q_grid = rho_m.sum();
- size_type Total_particles = 0;
- size_type local_particles = this->getLocalNum();
+ size_type Total_particles = 0;
+ size_type local_particles = this->getLocalNum();
- MPI_Reduce(&local_particles, &Total_particles, 1,
- MPI_UNSIGNED_LONG, MPI_SUM, 0, Ippl::getComm());
+ MPI_Reduce(&local_particles, &Total_particles, 1, MPI_UNSIGNED_LONG, MPI_SUM, 0,
+ Ippl::getComm());
- double rel_error = std::fabs((Q_m-Q_grid)/Q_m);
- m << "Rel. error in charge conservation = " << rel_error << endl;
+ double rel_error = std::fabs((Q_m - Q_grid) / Q_m);
+ m << "Rel. error in charge conservation = " << rel_error << endl;
- if(Ippl::Comm->rank() == 0) {
- if(Total_particles != totalP || rel_error > 1e-10) {
- m << "Time step: " << iteration << endl;
- m << "Total particles in the sim. " << totalP
- << " " << "after update: "
- << Total_particles << endl;
- m << "Rel. error in charge conservation: "
- << rel_error << endl;
- std::abort();
- }
- }
+ if (Ippl::Comm->rank() == 0) {
+ if (Total_particles != totalP || rel_error > 1e-10) {
+ m << "Time step: " << iteration << endl;
+ m << "Total particles in the sim. " << totalP << " "
+ << "after update: " << Total_particles << endl;
+ m << "Rel. error in charge conservation: " << rel_error << endl;
+ std::abort();
+ }
+ }
- rho_m = rho_m / (hrField[0] * hrField[1] * hrField[2]);
+ rho_m = rho_m / (hrField[0] * hrField[1] * hrField[2]);
- rhoNorm_m = norm(rho_m);
- IpplTimings::stopTimer(sumTimer);
+ rhoNorm_m = norm(rho_m);
+ IpplTimings::stopTimer(sumTimer);
- //dumpVTK(rho_m,nr_m[0],nr_m[1],nr_m[2],iteration,hrField[0],hrField[1],hrField[2]);
+ // dumpVTK(rho_m,nr_m[0],nr_m[1],nr_m[2],iteration,hrField[0],hrField[1],hrField[2]);
- //rho = rho_e - rho_i
- rho_m = rho_m - (Q_m/((rmax_m[0] - rmin_m[0]) * (rmax_m[1] - rmin_m[1]) * (rmax_m[2] - rmin_m[2])));
+ // rho = rho_e - rho_i
+ rho_m =
+ rho_m
+ - (Q_m / ((rmax_m[0] - rmin_m[0]) * (rmax_m[1] - rmin_m[1]) * (rmax_m[2] - rmin_m[2])));
}
void initSolver() {
-
Inform m("solver ");
- if(stype_m == "FFT")
+ if (stype_m == "FFT")
initFFTSolver();
else
m << "No solver matches the argument" << endl;
-
}
void initFFTSolver() {
ippl::ParameterList sp;
sp.add("output_type", Solver_t::GRAD);
- sp.add("use_heffte_defaults", false);
- sp.add("use_pencils", true);
- sp.add("use_reorder", false);
- sp.add("use_gpu_aware", true);
- sp.add("comm", ippl::p2p_pl);
- sp.add("r2c_direction", 0);
+ sp.add("use_heffte_defaults", false);
+ sp.add("use_pencils", true);
+ sp.add("use_reorder", false);
+ sp.add("use_gpu_aware", true);
+ sp.add("comm", ippl::p2p_pl);
+ sp.add("r2c_direction", 0);
solver_mp = std::make_shared();
@@ -376,46 +343,41 @@ class ChargedParticles : public ippl::ParticleBase {
solver_mp->setLhs(E_m);
}
-
-
- void dumpData() {
-
+ void dumpData() {
auto Pview = P.getView();
double Energy = 0.0;
- Kokkos::parallel_reduce("Particle Energy", this->getLocalNum(),
- KOKKOS_LAMBDA(const int i, double& valL){
- double myVal = dot(Pview(i), Pview(i)).apply();
- valL += myVal;
- }, Kokkos::Sum(Energy));
+ Kokkos::parallel_reduce(
+ "Particle Energy", this->getLocalNum(),
+ KOKKOS_LAMBDA(const int i, double& valL) {
+ double myVal = dot(Pview(i), Pview(i)).apply();
+ valL += myVal;
+ },
+ Kokkos::Sum(Energy));
Energy *= 0.5;
double gEnergy = 0.0;
- MPI_Reduce(&Energy, &gEnergy, 1,
- MPI_DOUBLE, MPI_SUM, 0, Ippl::getComm());
-
+ MPI_Reduce(&Energy, &gEnergy, 1, MPI_DOUBLE, MPI_SUM, 0, Ippl::getComm());
const int nghostE = E_m.getNghost();
- auto Eview = E_m.getView();
+ auto Eview = E_m.getView();
Vector_t normE;
using mdrange_type = Kokkos::MDRangePolicy>;
- for (unsigned d=0; d(temp));
+ for (unsigned d = 0; d < Dim; ++d) {
+ double temp = 0.0;
+ Kokkos::parallel_reduce(
+ "Vector E reduce",
+ mdrange_type({nghostE, nghostE, nghostE},
+ {Eview.extent(0) - nghostE, Eview.extent(1) - nghostE,
+ Eview.extent(2) - nghostE}),
+ KOKKOS_LAMBDA(const size_t i, const size_t j, const size_t k, double& valL) {
+ double myVal = std::pow(Eview(i, j, k)[d], 2);
+ valL += myVal;
+ },
+ Kokkos::Sum(temp));
double globaltemp = 0.0;
MPI_Reduce(&temp, &globaltemp, 1, MPI_DOUBLE, MPI_SUM, 0, Ippl::getComm());
normE[d] = std::sqrt(globaltemp);
@@ -431,149 +393,130 @@ class ChargedParticles : public ippl::ParticleBase {
csvout.precision(10);
csvout.setf(std::ios::scientific, std::ios::floatfield);
- if(time_m == 0.0) {
+ if (time_m == 0.0) {
csvout << "time, Kinetic energy, Rho_norm2, Ex_norm2, Ey_norm2, Ez_norm2" << endl;
}
- csvout << time_m << " "
- << gEnergy << " "
- << rhoNorm_m << " "
- << normE[0] << " "
- << normE[1] << " "
- << normE[2] << endl;
+ csvout << time_m << " " << gEnergy << " " << rhoNorm_m << " " << normE[0] << " "
+ << normE[1] << " " << normE[2] << endl;
}
Ippl::Comm->barrier();
- }
-
- void dumpLandau() {
+ }
+ void dumpLandau() {
const int nghostE = E_m.getNghost();
- auto Eview = E_m.getView();
+ auto Eview = E_m.getView();
double fieldEnergy, ExAmp;
using mdrange_type = Kokkos::MDRangePolicy>;
double temp = 0.0;
- Kokkos::parallel_reduce("Ex inner product",
- mdrange_type({nghostE, nghostE, nghostE},
- {Eview.extent(0) - nghostE,
- Eview.extent(1) - nghostE,
- Eview.extent(2) - nghostE}),
- KOKKOS_LAMBDA(const size_t i, const size_t j,
- const size_t k, double& valL)
- {
- double myVal = std::pow(Eview(i, j, k)[0], 2);
- valL += myVal;
- }, Kokkos::Sum(temp));
+ Kokkos::parallel_reduce(
+ "Ex inner product",
+ mdrange_type(
+ {nghostE, nghostE, nghostE},
+ {Eview.extent(0) - nghostE, Eview.extent(1) - nghostE, Eview.extent(2) - nghostE}),
+ KOKKOS_LAMBDA(const size_t i, const size_t j, const size_t k, double& valL) {
+ double myVal = std::pow(Eview(i, j, k)[0], 2);
+ valL += myVal;
+ },
+ Kokkos::Sum(temp));
double globaltemp = 0.0;
MPI_Reduce(&temp, &globaltemp, 1, MPI_DOUBLE, MPI_SUM, 0, Ippl::getComm());
fieldEnergy = globaltemp * hr_m[0] * hr_m[1] * hr_m[2];
double tempMax = 0.0;
- Kokkos::parallel_reduce("Ex max norm",
- mdrange_type({nghostE, nghostE, nghostE},
- {Eview.extent(0) - nghostE,
- Eview.extent(1) - nghostE,
- Eview.extent(2) - nghostE}),
- KOKKOS_LAMBDA(const size_t i, const size_t j,
- const size_t k, double& valL)
- {
- double myVal = std::fabs(Eview(i, j, k)[0]);
- if(myVal > valL) valL = myVal;
- }, Kokkos::Max(tempMax));
+ Kokkos::parallel_reduce(
+ "Ex max norm",
+ mdrange_type(
+ {nghostE, nghostE, nghostE},
+ {Eview.extent(0) - nghostE, Eview.extent(1) - nghostE, Eview.extent(2) - nghostE}),
+ KOKKOS_LAMBDA(const size_t i, const size_t j, const size_t k, double& valL) {
+ double myVal = std::fabs(Eview(i, j, k)[0]);
+ if (myVal > valL)
+ valL = myVal;
+ },
+ Kokkos::Max(tempMax));
ExAmp = 0.0;
MPI_Reduce(&tempMax, &ExAmp, 1, MPI_DOUBLE, MPI_MAX, 0, Ippl::getComm());
-
if (Ippl::Comm->rank() == 0) {
std::stringstream fname;
fname << "data/FieldLandau_";
fname << Ippl::Comm->size();
fname << ".csv";
-
Inform csvout(NULL, fname.str().c_str(), Inform::APPEND);
csvout.precision(10);
csvout.setf(std::ios::scientific, std::ios::floatfield);
- if(time_m == 0.0) {
+ if (time_m == 0.0) {
csvout << "time, Ex_field_energy, Ex_max_norm" << endl;
}
- csvout << time_m << " "
- << fieldEnergy << " "
- << ExAmp << endl;
-
+ csvout << time_m << " " << fieldEnergy << " " << ExAmp << endl;
}
-
+
Ippl::Comm->barrier();
- }
-
- void dumpBumponTail() {
+ }
+ void dumpBumponTail() {
const int nghostE = E_m.getNghost();
- auto Eview = E_m.getView();
+ auto Eview = E_m.getView();
double fieldEnergy, EzAmp;
using mdrange_type = Kokkos::MDRangePolicy>;
double temp = 0.0;
- Kokkos::parallel_reduce("Ex inner product",
- mdrange_type({nghostE, nghostE, nghostE},
- {Eview.extent(0) - nghostE,
- Eview.extent(1) - nghostE,
- Eview.extent(2) - nghostE}),
- KOKKOS_LAMBDA(const size_t i, const size_t j,
- const size_t k, double& valL)
- {
- double myVal = std::pow(Eview(i, j, k)[2], 2);
- valL += myVal;
- }, Kokkos::Sum(temp));
+ Kokkos::parallel_reduce(
+ "Ex inner product",
+ mdrange_type(
+ {nghostE, nghostE, nghostE},
+ {Eview.extent(0) - nghostE, Eview.extent(1) - nghostE, Eview.extent(2) - nghostE}),
+ KOKKOS_LAMBDA(const size_t i, const size_t j, const size_t k, double& valL) {
+ double myVal = std::pow(Eview(i, j, k)[2], 2);
+ valL += myVal;
+ },
+ Kokkos::Sum(temp));
double globaltemp = 0.0;
MPI_Reduce(&temp, &globaltemp, 1, MPI_DOUBLE, MPI_SUM, 0, Ippl::getComm());
fieldEnergy = globaltemp * hr_m[0] * hr_m[1] * hr_m[2];
double tempMax = 0.0;
- Kokkos::parallel_reduce("Ex max norm",
- mdrange_type({nghostE, nghostE, nghostE},
- {Eview.extent(0) - nghostE,
- Eview.extent(1) - nghostE,
- Eview.extent(2) - nghostE}),
- KOKKOS_LAMBDA(const size_t i, const size_t j,
- const size_t k, double& valL)
- {
- double myVal = std::fabs(Eview(i, j, k)[2]);
- if(myVal > valL) valL = myVal;
- }, Kokkos::Max(tempMax));
+ Kokkos::parallel_reduce(
+ "Ex max norm",
+ mdrange_type(
+ {nghostE, nghostE, nghostE},
+ {Eview.extent(0) - nghostE, Eview.extent(1) - nghostE, Eview.extent(2) - nghostE}),
+ KOKKOS_LAMBDA(const size_t i, const size_t j, const size_t k, double& valL) {
+ double myVal = std::fabs(Eview(i, j, k)[2]);
+ if (myVal > valL)
+ valL = myVal;
+ },
+ Kokkos::Max(tempMax));
EzAmp = 0.0;
MPI_Reduce(&tempMax, &EzAmp, 1, MPI_DOUBLE, MPI_MAX, 0, Ippl::getComm());
-
if (Ippl::Comm->rank() == 0) {
std::stringstream fname;
fname << "data/FieldBumponTail_";
fname << Ippl::Comm->size();
fname << ".csv";
-
Inform csvout(NULL, fname.str().c_str(), Inform::APPEND);
csvout.precision(10);
csvout.setf(std::ios::scientific, std::ios::floatfield);
- if(time_m == 0.0) {
+ if (time_m == 0.0) {
csvout << "time, Ez_field_energy, Ez_max_norm" << endl;
}
- csvout << time_m << " "
- << fieldEnergy << " "
- << EzAmp << endl;
-
+ csvout << time_m << " " << fieldEnergy << " " << EzAmp << endl;
}
-
- Ippl::Comm->barrier();
- }
- void dumpParticleData() {
+ Ippl::Comm->barrier();
+ }
+ void dumpParticleData() {
typename ParticleAttrib::HostMirror R_host = this->R.getHostMirror();
typename ParticleAttrib::HostMirror P_host = this->P.getHostMirror();
Kokkos::deep_copy(R_host, this->R.getView());
@@ -586,39 +529,31 @@ class ChargedParticles : public ippl::ParticleBase {
pcsvout.precision(10);
pcsvout.setf(std::ios::scientific, std::ios::floatfield);
pcsvout << "R_x, R_y, R_z, V_x, V_y, V_z" << endl;
- for (size_type i = 0; i< this->getLocalNum(); i++) {
- pcsvout << R_host(i)[0] << " "
- << R_host(i)[1] << " "
- << R_host(i)[2] << " "
- << P_host(i)[0] << " "
- << P_host(i)[1] << " "
- << P_host(i)[2] << endl;
+ for (size_type i = 0; i < this->getLocalNum(); i++) {
+ pcsvout << R_host(i)[0] << " " << R_host(i)[1] << " " << R_host(i)[2] << " "
+ << P_host(i)[0] << " " << P_host(i)[1] << " " << P_host(i)[2] << endl;
}
Ippl::Comm->barrier();
- }
-
- void dumpLocalDomains(const FieldLayout_t& fl, const unsigned int step) {
+ }
+ void dumpLocalDomains(const FieldLayout_t& fl, const unsigned int step) {
if (Ippl::Comm->rank() == 0) {
const typename FieldLayout_t::host_mirror_type domains = fl.getHostLocalDomains();
std::ofstream myfile;
myfile.open("data/domains" + std::to_string(step) + ".txt");
for (unsigned int i = 0; i < domains.size(); ++i) {
- myfile << domains[i][0].first() << " " << domains[i][1].first() << " " << domains[i][2].first() << " "
- << domains[i][0].first() << " " << domains[i][1].last() << " " << domains[i][2].first() << " "
- << domains[i][0].last() << " " << domains[i][1].first() << " " << domains[i][2].first() << " "
- << domains[i][0].first() << " " << domains[i][1].first() << " " << domains[i][2].last()
- << "\n";
+ myfile << domains[i][0].first() << " " << domains[i][1].first() << " "
+ << domains[i][2].first() << " " << domains[i][0].first() << " "
+ << domains[i][1].last() << " " << domains[i][2].first() << " "
+ << domains[i][0].last() << " " << domains[i][1].first() << " "
+ << domains[i][2].first() << " " << domains[i][0].first() << " "
+ << domains[i][1].first() << " " << domains[i][2].last() << "\n";
}
myfile.close();
}
Ippl::Comm->barrier();
- }
-
-private:
- void setBCAllPeriodic() {
-
- this->setParticleBC(ippl::BC::PERIODIC);
}
+private:
+ void setBCAllPeriodic() { this->setParticleBC(ippl::BC::PERIODIC); }
};
diff --git a/alpine/LandauDamping.cpp b/alpine/LandauDamping.cpp
index 26f8a3a88..068711ef4 100644
--- a/alpine/LandauDamping.cpp
+++ b/alpine/LandauDamping.cpp
@@ -9,7 +9,7 @@
// stype = Field solver type e.g., FFT
// lbthres = Load balancing threshold i.e., lbthres*100 is the maximum load imbalance
// percentage which can be tolerated and beyond which
-// particle load balancing occurs. A value of 0.01 is good for many typical
+// particle load balancing occurs. A value of 0.01 is good for many typical
// simulations.
// ovfactor = Over-allocation factor for the buffers used in the communication. Typical
// values are 1.0, 2.0. Value 1.0 means no over-allocation.
@@ -31,113 +31,108 @@
// along with IPPL. If not, see .
//
-#include "ChargedParticles.hpp"
-#include
-#include
-#include
+#include
#include
+#include
#include
-#include
+#include
+#include
+#include "ChargedParticles.hpp"
-#include
+#include
#include
#include "Utility/IpplTimings.h"
template
struct Newton1D {
+ double tol = 1e-12;
+ int max_iter = 20;
+ double pi = std::acos(-1.0);
- double tol = 1e-12;
- int max_iter = 20;
- double pi = std::acos(-1.0);
-
- T k, alpha, u;
-
- KOKKOS_INLINE_FUNCTION
- Newton1D() {}
-
- KOKKOS_INLINE_FUNCTION
- Newton1D(const T& k_, const T& alpha_,
- const T& u_)
- : k(k_), alpha(alpha_), u(u_) {}
-
- KOKKOS_INLINE_FUNCTION
- ~Newton1D() {}
-
- KOKKOS_INLINE_FUNCTION
- T f(T& x) {
- T F;
- F = x + (alpha * (std::sin(k * x) / k)) - u;
- return F;
- }
-
- KOKKOS_INLINE_FUNCTION
- T fprime(T& x) {
- T Fprime;
- Fprime = 1 + (alpha * std::cos(k * x));
- return Fprime;
- }
-
- KOKKOS_FUNCTION
- void solve(T& x) {
- int iterations = 0;
- while (iterations < max_iter && std::fabs(f(x)) > tol) {
- x = x - (f(x)/fprime(x));
- iterations += 1;
- }
- }
-};
-
-
-template
-struct generate_random {
+ T k, alpha, u;
- using view_type = typename ippl::detail::ViewType::view_type;
- using value_type = typename T::value_type;
- // Output View for the random numbers
- view_type x, v;
+ KOKKOS_INLINE_FUNCTION Newton1D() {}
- // The GeneratorPool
- GeneratorPool rand_pool;
+ KOKKOS_INLINE_FUNCTION Newton1D(const T& k_, const T& alpha_, const T& u_)
+ : k(k_)
+ , alpha(alpha_)
+ , u(u_) {}
- value_type alpha;
+ KOKKOS_INLINE_FUNCTION ~Newton1D() {}
- T k, minU, maxU;
+ KOKKOS_INLINE_FUNCTION T f(T& x) {
+ T F;
+ F = x + (alpha * (std::sin(k * x) / k)) - u;
+ return F;
+ }
- // Initialize all members
- generate_random(view_type x_, view_type v_, GeneratorPool rand_pool_,
- value_type& alpha_, T& k_, T& minU_, T& maxU_)
- : x(x_), v(v_), rand_pool(rand_pool_),
- alpha(alpha_), k(k_), minU(minU_), maxU(maxU_) {}
+ KOKKOS_INLINE_FUNCTION T fprime(T& x) {
+ T Fprime;
+ Fprime = 1 + (alpha * std::cos(k * x));
+ return Fprime;
+ }
- KOKKOS_INLINE_FUNCTION
- void operator()(const size_t i) const {
- // Get a random number state from the pool for the active thread
- typename GeneratorPool::generator_type rand_gen = rand_pool.get_state();
+ KOKKOS_FUNCTION
+ void solve(T& x) {
+ int iterations = 0;
+ while (iterations < max_iter && std::fabs(f(x)) > tol) {
+ x = x - (f(x) / fprime(x));
+ iterations += 1;
+ }
+ }
+};
- value_type u;
- for (unsigned d = 0; d < Dim; ++d) {
+template
+struct generate_random {
+ using view_type = typename ippl::detail::ViewType::view_type;
+ using value_type = typename T::value_type;
+ // Output View for the random numbers
+ view_type x, v;
+
+ // The GeneratorPool
+ GeneratorPool rand_pool;
+
+ value_type alpha;
+
+ T k, minU, maxU;
+
+ // Initialize all members
+ generate_random(view_type x_, view_type v_, GeneratorPool rand_pool_, value_type& alpha_, T& k_,
+ T& minU_, T& maxU_)
+ : x(x_)
+ , v(v_)
+ , rand_pool(rand_pool_)
+ , alpha(alpha_)
+ , k(k_)
+ , minU(minU_)
+ , maxU(maxU_) {}
+
+ KOKKOS_INLINE_FUNCTION void operator()(const size_t i) const {
+ // Get a random number state from the pool for the active thread
+ typename GeneratorPool::generator_type rand_gen = rand_pool.get_state();
+
+ value_type u;
+ for (unsigned d = 0; d < Dim; ++d) {
+ u = rand_gen.drand(minU[d], maxU[d]);
+ x(i)[d] = u / (1 + alpha);
+ Newton1D solver(k[d], alpha, u);
+ solver.solve(x(i)[d]);
+ v(i)[d] = rand_gen.normal(0.0, 1.0);
+ }
- u = rand_gen.drand(minU[d], maxU[d]);
- x(i)[d] = u / (1 + alpha);
- Newton1D solver(k[d], alpha, u);
- solver.solve(x(i)[d]);
- v(i)[d] = rand_gen.normal(0.0, 1.0);
+ // Give the state back, which will allow another thread to acquire it
+ rand_pool.free_state(rand_gen);
}
-
- // Give the state back, which will allow another thread to acquire it
- rand_pool.free_state(rand_gen);
- }
};
double CDF(const double& x, const double& alpha, const double& k) {
- double cdf = x + (alpha / k) * std::sin(k * x);
- return cdf;
+ double cdf = x + (alpha / k) * std::sin(k * x);
+ return cdf;
}
KOKKOS_FUNCTION
-double PDF(const Vector_t& xvec, const double& alpha,
- const Vector_t& kw, const unsigned Dim) {
+double PDF(const Vector_t& xvec, const double& alpha, const Vector_t& kw, const unsigned Dim) {
double pdf = 1.0;
for (unsigned d = 0; d < Dim; ++d) {
@@ -148,48 +143,40 @@ double PDF(const Vector_t& xvec, const double& alpha,
const char* TestName = "LandauDamping";
-int main(int argc, char *argv[]){
+int main(int argc, char* argv[]) {
Ippl ippl(argc, argv);
-
+
Inform msg("LandauDamping");
- Inform msg2all("LandauDamping",INFORM_ALL_NODES);
+ Inform msg2all("LandauDamping", INFORM_ALL_NODES);
Ippl::Comm->setDefaultOverallocation(std::atof(argv[8]));
- auto start = std::chrono::high_resolution_clock::now();
- ippl::Vector nr = {
- std::atoi(argv[1]),
- std::atoi(argv[2]),
- std::atoi(argv[3])
- };
-
- static IpplTimings::TimerRef mainTimer = IpplTimings::getTimer("total");
- static IpplTimings::TimerRef particleCreation = IpplTimings::getTimer("particlesCreation");
- static IpplTimings::TimerRef dumpDataTimer = IpplTimings::getTimer("dumpData");
- static IpplTimings::TimerRef PTimer = IpplTimings::getTimer("pushVelocity");
- static IpplTimings::TimerRef RTimer = IpplTimings::getTimer("pushPosition");
- static IpplTimings::TimerRef updateTimer = IpplTimings::getTimer("update");
- static IpplTimings::TimerRef DummySolveTimer = IpplTimings::getTimer("solveWarmup");
- static IpplTimings::TimerRef SolveTimer = IpplTimings::getTimer("solve");
+ auto start = std::chrono::high_resolution_clock::now();
+ ippl::Vector nr = {std::atoi(argv[1]), std::atoi(argv[2]), std::atoi(argv[3])};
+
+ static IpplTimings::TimerRef mainTimer = IpplTimings::getTimer("total");
+ static IpplTimings::TimerRef particleCreation = IpplTimings::getTimer("particlesCreation");
+ static IpplTimings::TimerRef dumpDataTimer = IpplTimings::getTimer("dumpData");
+ static IpplTimings::TimerRef PTimer = IpplTimings::getTimer("pushVelocity");
+ static IpplTimings::TimerRef RTimer = IpplTimings::getTimer("pushPosition");
+ static IpplTimings::TimerRef updateTimer = IpplTimings::getTimer("update");
+ static IpplTimings::TimerRef DummySolveTimer = IpplTimings::getTimer("solveWarmup");
+ static IpplTimings::TimerRef SolveTimer = IpplTimings::getTimer("solve");
static IpplTimings::TimerRef domainDecomposition = IpplTimings::getTimer("loadBalance");
IpplTimings::startTimer(mainTimer);
const size_type totalP = std::atoll(argv[4]);
- const unsigned int nt = std::atoi(argv[5]);
+ const unsigned int nt = std::atoi(argv[5]);
- msg << "Landau damping"
- << endl
- << "nt " << nt << " Np= "
- << totalP << " grid = " << nr
- << endl;
+ msg << "Landau damping" << endl << "nt " << nt << " Np= " << totalP << " grid = " << nr << endl;
using bunch_type = ChargedParticles;
- std::unique_ptr P;
+ std::unique_ptr P;
ippl::NDIndex domain;
- for (unsigned i = 0; i< Dim; i++) {
+ for (unsigned i = 0; i < Dim; i++) {
domain[i] = ippl::Index(nr[i]);
}
@@ -199,26 +186,26 @@ int main(int argc, char *argv[]){
}
// create mesh and layout objects for this problem domain
- Vector_t kw = {0.5, 0.5, 0.5};
+ Vector_t kw = {0.5, 0.5, 0.5};
double alpha = 0.05;
Vector_t rmin(0.0);
- Vector_t rmax = 2 * pi / kw ;
- double dx = rmax[0] / nr[0];
- double dy = rmax[1] / nr[1];
- double dz = rmax[2] / nr[2];
+ Vector_t rmax = 2 * pi / kw;
+ double dx = rmax[0] / nr[0];
+ double dy = rmax[1] / nr[1];
+ double dz = rmax[2] / nr[2];
- Vector_t hr = {dx, dy, dz};
+ Vector_t hr = {dx, dy, dz};
Vector_t origin = {rmin[0], rmin[1], rmin[2]};
- const double dt = 0.5*dx;
+ const double dt = 0.5 * dx;
- const bool isAllPeriodic=true;
+ const bool isAllPeriodic = true;
Mesh_t mesh(domain, hr, origin);
FieldLayout_t FL(domain, decomp, isAllPeriodic);
PLayout_t PL(FL, mesh);
- //Q = -\int\int f dx dv
+ // Q = -\int\int f dx dv
double Q = -rmax[0] * rmax[1] * rmax[2];
- P = std::make_unique(PL,hr,rmin,rmax,decomp,Q);
+ P = std::make_unique(PL, hr, rmin, rmax, decomp, Q);
P->nr_m = nr;
@@ -229,7 +216,7 @@ int main(int argc, char *argv[]){
P->stype_m = argv[6];
P->initSolver();
- P->time_m = 0.0;
+ P->time_m = 0.0;
P->loadbalancethreshold_m = std::atof(argv[7]);
bool isFirstRepartition;
@@ -237,85 +224,80 @@ int main(int argc, char *argv[]){
if ((P->loadbalancethreshold_m != 1.0) && (Ippl::Comm->size() > 1)) {
msg << "Starting first repartition" << endl;
IpplTimings::startTimer(domainDecomposition);
- isFirstRepartition = true;
+ isFirstRepartition = true;
const ippl::NDIndex& lDom = FL.getLocalNDIndex();
- const int nghost = P->rho_m.getNghost();
- using mdrange_type = Kokkos::MDRangePolicy>;
- auto rhoview = P->rho_m.getView();
-
- Kokkos::parallel_for("Assign initial rho based on PDF",
- mdrange_type({nghost, nghost, nghost},
- {rhoview.extent(0) - nghost,
- rhoview.extent(1) - nghost,
- rhoview.extent(2) - nghost}),
- KOKKOS_LAMBDA(const int i,
- const int j,
- const int k)
- {
- //local to global index conversion
- const size_t ig = i + lDom[0].first() - nghost;
- const size_t jg = j + lDom[1].first() - nghost;
- const size_t kg = k + lDom[2].first() - nghost;
- double x = (ig + 0.5) * hr[0] + origin[0];
- double y = (jg + 0.5) * hr[1] + origin[1];
- double z = (kg + 0.5) * hr[2] + origin[2];
-
- Vector_t xvec = {x, y, z};
-
- rhoview(i, j, k) = PDF(xvec, alpha, kw, Dim);
-
- });
+ const int nghost = P->rho_m.getNghost();
+ using mdrange_type = Kokkos::MDRangePolicy>;
+ auto rhoview = P->rho_m.getView();
+
+ Kokkos::parallel_for(
+ "Assign initial rho based on PDF",
+ mdrange_type({nghost, nghost, nghost},
+ {rhoview.extent(0) - nghost, rhoview.extent(1) - nghost,
+ rhoview.extent(2) - nghost}),
+ KOKKOS_LAMBDA(const int i, const int j, const int k) {
+ // local to global index conversion
+ const size_t ig = i + lDom[0].first() - nghost;
+ const size_t jg = j + lDom[1].first() - nghost;
+ const size_t kg = k + lDom[2].first() - nghost;
+ double x = (ig + 0.5) * hr[0] + origin[0];
+ double y = (jg + 0.5) * hr[1] + origin[1];
+ double z = (kg + 0.5) * hr[2] + origin[2];
+
+ Vector_t xvec = {x, y, z};
+
+ rhoview(i, j, k) = PDF(xvec, alpha, kw, Dim);
+ });
Kokkos::fence();
-
+
P->initializeORB(FL, mesh);
P->repartition(FL, mesh, bunchBuffer, isFirstRepartition);
IpplTimings::stopTimer(domainDecomposition);
}
-
+
msg << "First domain decomposition done" << endl;
IpplTimings::startTimer(particleCreation);
typedef ippl::detail::RegionLayout RegionLayout_t;
- const RegionLayout_t& RLayout = PL.getRegionLayout();
+ const RegionLayout_t& RLayout = PL.getRegionLayout();
const typename RegionLayout_t::host_mirror_type Regions = RLayout.gethLocalRegions();
Vector_t Nr, Dr, minU, maxU;
int myRank = Ippl::Comm->rank();
- for (unsigned d = 0; d rank() < rest )
+ if (Ippl::Comm->rank() < rest)
++nloc;
P->create(nloc);
- Kokkos::Random_XorShift64_Pool<> rand_pool64((size_type)(42 + 100*Ippl::Comm->rank()));
+ Kokkos::Random_XorShift64_Pool<> rand_pool64((size_type)(42 + 100 * Ippl::Comm->rank()));
Kokkos::parallel_for(nloc,
generate_random, Dim>(
- P->R.getView(), P->P.getView(), rand_pool64, alpha, kw, minU, maxU));
+ P->R.getView(), P->P.getView(), rand_pool64, alpha, kw, minU, maxU));
Kokkos::fence();
Ippl::Comm->barrier();
- IpplTimings::stopTimer(particleCreation);
-
- P->q = P->Q_m/totalP;
+ IpplTimings::stopTimer(particleCreation);
+
+ P->q = P->Q_m / totalP;
msg << "particles created and initial conditions assigned " << endl;
isFirstRepartition = false;
- //The update after the particle creation is not needed as the
- //particles are generated locally
+ // The update after the particle creation is not needed as the
+ // particles are generated locally
IpplTimings::startTimer(DummySolveTimer);
P->rho_m = 0.0;
@@ -333,13 +315,12 @@ int main(int argc, char *argv[]){
IpplTimings::startTimer(dumpDataTimer);
P->dumpLandau();
P->gatherStatistics(totalP);
- //P->dumpLocalDomains(FL, 0);
+ // P->dumpLocalDomains(FL, 0);
IpplTimings::stopTimer(dumpDataTimer);
// begin main timestep loop
msg << "Starting iterations ..." << endl;
- for (unsigned int it=0; itP = P->P - 0.5 * dt * P->E;
IpplTimings::stopTimer(PTimer);
- //drift
+ // drift
IpplTimings::startTimer(RTimer);
P->R = P->R + dt * P->P;
IpplTimings::stopTimer(RTimer);
- //Since the particles have moved spatially update them to correct processors
- IpplTimings::startTimer(updateTimer);
+ // Since the particles have moved spatially update them to correct processors
+ IpplTimings::startTimer(updateTimer);
PL.update(*P, bunchBuffer);
IpplTimings::stopTimer(updateTimer);
// Domain Decomposition
- if (P->balance(totalP, it+1)) {
- msg << "Starting repartition" << endl;
- IpplTimings::startTimer(domainDecomposition);
- P->repartition(FL, mesh, bunchBuffer, isFirstRepartition);
- IpplTimings::stopTimer(domainDecomposition);
- //IpplTimings::startTimer(dumpDataTimer);
- //P->dumpLocalDomains(FL, it+1);
- //IpplTimings::stopTimer(dumpDataTimer);
+ if (P->balance(totalP, it + 1)) {
+ msg << "Starting repartition" << endl;
+ IpplTimings::startTimer(domainDecomposition);
+ P->repartition(FL, mesh, bunchBuffer, isFirstRepartition);
+ IpplTimings::stopTimer(domainDecomposition);
+ // IpplTimings::startTimer(dumpDataTimer);
+ // P->dumpLocalDomains(FL, it+1);
+ // IpplTimings::stopTimer(dumpDataTimer);
}
+ // scatter the charge onto the underlying grid
+ P->scatterCIC(totalP, it + 1, hr);
- //scatter the charge onto the underlying grid
- P->scatterCIC(totalP, it+1, hr);
-
- //Field solve
+ // Field solve
IpplTimings::startTimer(SolveTimer);
P->solver_mp->solve();
IpplTimings::stopTimer(SolveTimer);
@@ -383,7 +363,7 @@ int main(int argc, char *argv[]){
// gather E field
P->gatherCIC();
- //kick
+ // kick
IpplTimings::startTimer(PTimer);
P->P = P->P - 0.5 * dt * P->E;
IpplTimings::stopTimer(PTimer);
@@ -393,7 +373,7 @@ int main(int argc, char *argv[]){
P->dumpLandau();
P->gatherStatistics(totalP);
IpplTimings::stopTimer(dumpDataTimer);
- msg << "Finished time step: " << it+1 << " time: " << P->time_m << endl;
+ msg << "Finished time step: " << it + 1 << " time: " << P->time_m << endl;
}
msg << "LandauDamping: End." << endl;
@@ -402,9 +382,9 @@ int main(int argc, char *argv[]){
IpplTimings::print(std::string("timing.dat"));
auto end = std::chrono::high_resolution_clock::now();
- std::chrono::duration time_chrono = std::chrono::duration_cast>(end - start);
+ std::chrono::duration time_chrono =
+ std::chrono::duration_cast>(end - start);
std::cout << "Elapsed time: " << time_chrono.count() << std::endl;
-
return 0;
}
diff --git a/alpine/PenningTrap.cpp b/alpine/PenningTrap.cpp
index 87a9a6534..492cb0f76 100644
--- a/alpine/PenningTrap.cpp
+++ b/alpine/PenningTrap.cpp
@@ -9,14 +9,14 @@
// stype = Field solver type e.g., FFT
// lbthres = Load balancing threshold i.e., lbthres*100 is the maximum load imbalance
// percentage which can be tolerated and beyond which
-// particle load balancing occurs. A value of 0.01 is good for many typical
+// particle load balancing occurs. A value of 0.01 is good for many typical
// simulations.
// ovfactor = Over-allocation factor for the buffers used in the communication. Typical
// values are 1.0, 2.0. Value 1.0 means no over-allocation.
// Example:
// srun ./PenningTrap 128 128 128 10000 300 FFT 0.01 1.0 --info 10
//
-// Copyright (c) 2021, Sriramkrishnan Muralikrishnan,
+// Copyright (c) 2021, Sriramkrishnan Muralikrishnan,
// Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
@@ -32,170 +32,152 @@
//
#include "ChargedParticles.hpp"
-#include
-#include
-#include
+#include
#include
+#include
#include
-#include
+#include
+#include
-#include
+#include
#include
#include "Utility/IpplTimings.h"
template
struct Newton1D {
+ double tol = 1e-12;
+ int max_iter = 20;
+ double pi = std::acos(-1.0);
- double tol = 1e-12;
- int max_iter = 20;
- double pi = std::acos(-1.0);
-
- T mu, sigma, u;
-
- KOKKOS_INLINE_FUNCTION
- Newton1D() {}
-
- KOKKOS_INLINE_FUNCTION
- Newton1D(const T& mu_, const T& sigma_,
- const T& u_)
- : mu(mu_), sigma(sigma_), u(u_) {}
-
- KOKKOS_INLINE_FUNCTION
- ~Newton1D() {}
-
- KOKKOS_INLINE_FUNCTION
- T f(T& x) {
- T F;
- F = std::erf((x - mu)/(sigma * std::sqrt(2.0)))
- - 2 * u + 1;
- return F;
- }
-
- KOKKOS_INLINE_FUNCTION
- T fprime(T& x) {
- T Fprime;
- Fprime = (1 / sigma) * std::sqrt(2 / pi) *
- std::exp(-0.5 * (std::pow(((x - mu) / sigma),2)));
- return Fprime;
- }
-
- KOKKOS_FUNCTION
- void solve(T& x) {
- int iterations = 0;
- while ((iterations < max_iter) && (std::fabs(f(x)) > tol)) {
- x = x - (f(x)/fprime(x));
- iterations += 1;
- }
- }
-};
+ T mu, sigma, u;
+ KOKKOS_INLINE_FUNCTION Newton1D() {}
-template
-struct generate_random {
+ KOKKOS_INLINE_FUNCTION Newton1D(const T& mu_, const T& sigma_, const T& u_)
+ : mu(mu_)
+ , sigma(sigma_)
+ , u(u_) {}
- using view_type = typename ippl::detail::ViewType::view_type;
- using value_type = typename T::value_type;
- // Output View for the random numbers
- view_type x, v;
+ KOKKOS_INLINE_FUNCTION ~Newton1D() {}
- // The GeneratorPool
- GeneratorPool rand_pool;
+ KOKKOS_INLINE_FUNCTION T f(T& x) {
+ T F;
+ F = std::erf((x - mu) / (sigma * std::sqrt(2.0))) - 2 * u + 1;
+ return F;
+ }
- T mu, sigma, minU, maxU;
+ KOKKOS_INLINE_FUNCTION T fprime(T& x) {
+ T Fprime;
+ Fprime =
+ (1 / sigma) * std::sqrt(2 / pi) * std::exp(-0.5 * (std::pow(((x - mu) / sigma), 2)));
+ return Fprime;
+ }
- double pi = std::acos(-1.0);
+ KOKKOS_FUNCTION
+ void solve(T& x) {
+ int iterations = 0;
+ while ((iterations < max_iter) && (std::fabs(f(x)) > tol)) {
+ x = x - (f(x) / fprime(x));
+ iterations += 1;
+ }
+ }
+};
- // Initialize all members
- generate_random(view_type x_, view_type v_, GeneratorPool rand_pool_,
- T& mu_, T& sigma_, T& minU_, T& maxU_)
- : x(x_), v(v_), rand_pool(rand_pool_),
- mu(mu_), sigma(sigma_), minU(minU_), maxU(maxU_) {}
+template
+struct generate_random {
+ using view_type = typename ippl::detail::ViewType::view_type;
+ using value_type = typename T::value_type;
+ // Output View for the random numbers
+ view_type x, v;
- KOKKOS_INLINE_FUNCTION
- void operator()(const size_t i) const {
- // Get a random number state from the pool for the active thread
- typename GeneratorPool::generator_type rand_gen = rand_pool.get_state();
+ // The GeneratorPool
+ GeneratorPool rand_pool;
- value_type u;
- for (unsigned d = 0; d < Dim; ++d) {
- u = rand_gen.drand(minU[d], maxU[d]);
- x(i)[d] = (std::sqrt(pi / 2) * (2 * u - 1)) *
- sigma[d] + mu[d];
- Newton1D solver(mu[d], sigma[d], u);
- solver.solve(x(i)[d]);
- v(i)[d] = rand_gen.normal(0.0, 1.0);
- }
+ T mu, sigma, minU, maxU;
- // Give the state back, which will allow another thread to acquire it
- rand_pool.free_state(rand_gen);
- }
+ double pi = std::acos(-1.0);
+
+ // Initialize all members
+ generate_random(view_type x_, view_type v_, GeneratorPool rand_pool_, T& mu_, T& sigma_,
+ T& minU_, T& maxU_)
+ : x(x_)
+ , v(v_)
+ , rand_pool(rand_pool_)
+ , mu(mu_)
+ , sigma(sigma_)
+ , minU(minU_)
+ , maxU(maxU_) {}
+
+ KOKKOS_INLINE_FUNCTION void operator()(const size_t i) const {
+ // Get a random number state from the pool for the active thread
+ typename GeneratorPool::generator_type rand_gen = rand_pool.get_state();
+
+ value_type u;
+ for (unsigned d = 0; d < Dim; ++d) {
+ u = rand_gen.drand(minU[d], maxU[d]);
+ x(i)[d] = (std::sqrt(pi / 2) * (2 * u - 1)) * sigma[d] + mu[d];
+ Newton1D solver(mu[d], sigma[d], u);
+ solver.solve(x(i)[d]);
+ v(i)[d] = rand_gen.normal(0.0, 1.0);
+ }
+
+ // Give the state back, which will allow another thread to acquire it
+ rand_pool.free_state(rand_gen);
+ }
};
double CDF(const double& x, const double& mu, const double& sigma) {
- double cdf = 0.5 * (1.0 + std::erf((x - mu)/(sigma * std::sqrt(2))));
- return cdf;
+ double cdf = 0.5 * (1.0 + std::erf((x - mu) / (sigma * std::sqrt(2))));
+ return cdf;
}
-
KOKKOS_FUNCTION
-double PDF(const Vector_t& xvec, const Vector_t&mu,
- const Vector_t& sigma, const unsigned Dim) {
+double PDF(const Vector_t& xvec, const Vector_t& mu, const Vector_t& sigma, const unsigned Dim) {
double pdf = 1.0;
- double pi = std::acos(-1.0);
+ double pi = std::acos(-1.0);
for (unsigned d = 0; d < Dim; ++d) {
- pdf *= (1.0/ (sigma[d] * std::sqrt(2 * pi))) *
- std::exp(-0.5 * std::pow((xvec[d] - mu[d])/sigma[d],2));
+ pdf *= (1.0 / (sigma[d] * std::sqrt(2 * pi)))
+ * std::exp(-0.5 * std::pow((xvec[d] - mu[d]) / sigma[d], 2));
}
return pdf;
}
const char* TestName = "PenningTrap";
-int main(int argc, char *argv[]){
+int main(int argc, char* argv[]) {
Ippl ippl(argc, argv);
Inform msg("PenningTrap");
- Inform msg2all("PenningTrap",INFORM_ALL_NODES);
-
+ Inform msg2all("PenningTrap", INFORM_ALL_NODES);
Ippl::Comm->setDefaultOverallocation(std::atof(argv[8]));
- auto start = std::chrono::high_resolution_clock::now();
- ippl::Vector nr = {
- std::atoi(argv[1]),
- std::atoi(argv[2]),
- std::atoi(argv[3])
- };
-
- static IpplTimings::TimerRef mainTimer = IpplTimings::getTimer("total");
- static IpplTimings::TimerRef particleCreation = IpplTimings::getTimer("particlesCreation");
- static IpplTimings::TimerRef dumpDataTimer = IpplTimings::getTimer("dumpData");
- static IpplTimings::TimerRef PTimer = IpplTimings::getTimer("pushVelocity");
- static IpplTimings::TimerRef RTimer = IpplTimings::getTimer("pushPosition");
- static IpplTimings::TimerRef updateTimer = IpplTimings::getTimer("update");
- static IpplTimings::TimerRef DummySolveTimer = IpplTimings::getTimer("solveWarmup");
- static IpplTimings::TimerRef SolveTimer = IpplTimings::getTimer("Solve");
+ auto start = std::chrono::high_resolution_clock::now();
+ ippl::Vector nr = {std::atoi(argv[1]), std::atoi(argv[2]), std::atoi(argv[3])};
+
+ static IpplTimings::TimerRef mainTimer = IpplTimings::getTimer("total");
+ static IpplTimings::TimerRef particleCreation = IpplTimings::getTimer("particlesCreation");
+ static IpplTimings::TimerRef dumpDataTimer = IpplTimings::getTimer("dumpData");
+ static IpplTimings::TimerRef PTimer = IpplTimings::getTimer("pushVelocity");
+ static IpplTimings::TimerRef RTimer = IpplTimings::getTimer("pushPosition");
+ static IpplTimings::TimerRef updateTimer = IpplTimings::getTimer("update");
+ static IpplTimings::TimerRef DummySolveTimer = IpplTimings::getTimer("solveWarmup");
+ static IpplTimings::TimerRef SolveTimer = IpplTimings::getTimer("Solve");
static IpplTimings::TimerRef domainDecomposition = IpplTimings::getTimer("loadBalance");
-
-
+
IpplTimings::startTimer(mainTimer);
-
- size_type totalP = std::atol(argv[4]);
- const unsigned int nt = std::atoi(argv[5]);
- msg << "Penning Trap "
- << endl
- << "nt " << nt << " Np= "
- << totalP << " grid = " << nr
- << endl;
+ size_type totalP = std::atol(argv[4]);
+ const unsigned int nt = std::atoi(argv[5]);
+ msg << "Penning Trap " << endl << "nt " << nt << " Np= " << totalP << " grid = " << nr << endl;
using bunch_type = ChargedParticles;
- std::unique_ptr P;
+ std::unique_ptr P;
ippl::NDIndex domain;
- for (unsigned i = 0; i< Dim; i++) {
+ for (unsigned i = 0; i < Dim; i++) {
domain[i] = ippl::Index(nr[i]);
}
@@ -211,137 +193,129 @@ int main(int argc, char *argv[]){
double dy = rmax[1] / nr[1];
double dz = rmax[2] / nr[2];
- Vector_t hr = {dx, dy, dz};
- Vector_t origin = {rmin[0], rmin[1], rmin[2]};
- unsigned int nrMax = 2048;// Max grid size in our studies
- double dxFinest = rmax[0] / nrMax;
- const double dt = 0.5 * dxFinest;//size of timestep
+ Vector_t hr = {dx, dy, dz};
+ Vector_t origin = {rmin[0], rmin[1], rmin[2]};
+ unsigned int nrMax = 2048; // Max grid size in our studies
+ double dxFinest = rmax[0] / nrMax;
+ const double dt = 0.5 * dxFinest; // size of timestep
- const bool isAllPeriodic=true;
+ const bool isAllPeriodic = true;
Mesh_t mesh(domain, hr, origin);
FieldLayout_t FL(domain, decomp, isAllPeriodic);
PLayout_t PL(FL, mesh);
- double Q = -1562.5;
+ double Q = -1562.5;
double Bext = 5.0;
- P = std::make_unique(PL,hr,rmin,rmax,decomp,Q);
+ P = std::make_unique(PL, hr, rmin, rmax, decomp, Q);
P->nr_m = nr;
-
Vector_t length = rmax - rmin;
Vector_t mu, sd;
- for (unsigned d = 0; dE_m.initialize(mesh, FL);
P->rho_m.initialize(mesh, FL);
bunch_type bunchBuffer(PL);
-
P->stype_m = argv[6];
P->initSolver();
- P->time_m = 0.0;
+ P->time_m = 0.0;
P->loadbalancethreshold_m = std::atof(argv[7]);
bool isFirstRepartition;
-
+
if ((P->loadbalancethreshold_m != 1.0) && (Ippl::Comm->size() > 1)) {
msg << "Starting first repartition" << endl;
IpplTimings::startTimer(domainDecomposition);
- isFirstRepartition = true;
+ isFirstRepartition = true;
const ippl::NDIndex& lDom = FL.getLocalNDIndex();
- const int nghost = P->rho_m.getNghost();
- using mdrange_type = Kokkos::MDRangePolicy>;
- auto rhoview = P->rho_m.getView();
-
- Kokkos::parallel_for("Assign initial rho based on PDF",
- mdrange_type({nghost, nghost, nghost},
- {rhoview.extent(0) - nghost,
- rhoview.extent(1) - nghost,
- rhoview.extent(2) - nghost}),
- KOKKOS_LAMBDA(const int i,
- const int j,
- const int k)
- {
- //local to global index conversion
- const size_t ig = i + lDom[0].first() - nghost;
- const size_t jg = j + lDom[1].first() - nghost;
- const size_t kg = k + lDom[2].first() - nghost;
- double x = (ig + 0.5) * hr[0] + origin[0];
- double y = (jg + 0.5) * hr[1] + origin[1];
- double z = (kg + 0.5) * hr[2] + origin[2];
-
- Vector_t xvec = {x, y, z};
-
- rhoview(i, j, k) = PDF(xvec, mu, sd, Dim);
-
- });
+ const int nghost = P->rho_m.getNghost();
+ using mdrange_type = Kokkos::MDRangePolicy>;
+ auto rhoview = P->rho_m.getView();
+
+ Kokkos::parallel_for(
+ "Assign initial rho based on PDF",
+ mdrange_type({nghost, nghost, nghost},
+ {rhoview.extent(0) - nghost, rhoview.extent(1) - nghost,
+ rhoview.extent(2) - nghost}),
+ KOKKOS_LAMBDA(const int i, const int j, const int k) {
+ // local to global index conversion
+ const size_t ig = i + lDom[0].first() - nghost;
+ const size_t jg = j + lDom[1].first() - nghost;
+ const size_t kg = k + lDom[2].first() - nghost;
+ double x = (ig + 0.5) * hr[0] + origin[0];
+ double y = (jg + 0.5) * hr[1] + origin[1];
+ double z = (kg + 0.5) * hr[2] + origin[2];
+
+ Vector_t xvec = {x, y, z};
+
+ rhoview(i, j, k) = PDF(xvec, mu, sd, Dim);
+ });
Kokkos::fence();
-
+
P->initializeORB(FL, mesh);
P->repartition(FL, mesh, bunchBuffer, isFirstRepartition);
IpplTimings::stopTimer(domainDecomposition);
}
-
msg << "First domain decomposition done" << endl;
IpplTimings::startTimer(particleCreation);
typedef ippl::detail::RegionLayout RegionLayout_t;
- const RegionLayout_t& RLayout = PL.getRegionLayout();
+ const RegionLayout_t& RLayout = PL.getRegionLayout();
const typename RegionLayout_t::host_mirror_type Regions = RLayout.gethLocalRegions();
Vector_t Nr, Dr, minU, maxU;
int myRank = Ippl::Comm->rank();
- for (unsigned d = 0; d rank() < rest )
+ if (Ippl::Comm->rank() < rest)
++nloc;
P->create(nloc);
- Kokkos::Random_XorShift64_Pool<> rand_pool64((size_type)(42 + 100*Ippl::Comm->rank()));
+ Kokkos::Random_XorShift64_Pool<> rand_pool64((size_type)(42 + 100 * Ippl::Comm->rank()));
Kokkos::parallel_for(nloc,
generate_random, Dim>(
- P->R.getView(), P->P.getView(), rand_pool64, mu, sd, minU, maxU));
+ P->R.getView(), P->P.getView(), rand_pool64, mu, sd, minU, maxU));
Kokkos::fence();
Ippl::Comm->barrier();
- IpplTimings::stopTimer(particleCreation);
-
- P->q = P->Q_m/totalP;
+ IpplTimings::stopTimer(particleCreation);
+
+ P->q = P->Q_m / totalP;
msg << "particles created and initial conditions assigned " << endl;
isFirstRepartition = false;
- //The update after the particle creation is not needed as the
- //particles are generated locally
-
+ // The update after the particle creation is not needed as the
+ // particles are generated locally
+
IpplTimings::startTimer(DummySolveTimer);
P->rho_m = 0.0;
P->solver_mp->solve();
IpplTimings::stopTimer(DummySolveTimer);
-
+
P->scatterCIC(totalP, 0, hr);
IpplTimings::startTimer(SolveTimer);
@@ -353,16 +327,15 @@ int main(int argc, char *argv[]){
IpplTimings::startTimer(dumpDataTimer);
P->dumpData();
P->gatherStatistics(totalP);
- //P->dumpLocalDomains(FL, 0);
+ // P->dumpLocalDomains(FL, 0);
IpplTimings::stopTimer(dumpDataTimer);
double alpha = -0.5 * dt;
- double DrInv = 1.0 / (1 + (std::pow((alpha * Bext), 2)));
+ double DrInv = 1.0 / (1 + (std::pow((alpha * Bext), 2)));
// begin main timestep loop
msg << "Starting iterations ..." << endl;
- for (unsigned int it=0; itR.getView();
auto Pview = P->P.getView();
auto Eview = P->E.getView();
- double V0 = 30*rmax[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*std::pow(rmax[2],2)));
- double Eext_y = -(Rview(j)[1] - 0.5*rmax[1]) * (V0/(2*std::pow(rmax[2],2)));
- double Eext_z = (Rview(j)[2] - 0.5*rmax[2]) * (V0/(std::pow(rmax[2],2)));
-
- Eview(j)[0] += Eext_x;
- Eview(j)[1] += Eext_y;
- Eview(j)[2] += Eext_z;
-
- Pview(j)[0] += alpha * (Eview(j)[0] + Pview(j)[1] * Bext);
- Pview(j)[1] += alpha * (Eview(j)[1] - Pview(j)[0] * Bext);
- Pview(j)[2] += alpha * Eview(j)[2];
- });
+ double V0 = 30 * rmax[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 * std::pow(rmax[2], 2)));
+ double Eext_y = -(Rview(j)[1] - 0.5 * rmax[1]) * (V0 / (2 * std::pow(rmax[2], 2)));
+ double Eext_z = (Rview(j)[2] - 0.5 * rmax[2]) * (V0 / (std::pow(rmax[2], 2)));
+
+ Eview(j)[0] += Eext_x;
+ Eview(j)[1] += Eext_y;
+ Eview(j)[2] += Eext_z;
+
+ Pview(j)[0] += alpha * (Eview(j)[0] + Pview(j)[1] * Bext);
+ Pview(j)[1] += alpha * (Eview(j)[1] - Pview(j)[0] * Bext);
+ Pview(j)[2] += alpha * Eview(j)[2];
+ });
IpplTimings::stopTimer(PTimer);
- //drift
+ // drift
IpplTimings::startTimer(RTimer);
P->R = P->R + dt * P->P;
IpplTimings::stopTimer(RTimer);
- //Since the particles have moved spatially update them to correct processors
- IpplTimings::startTimer(updateTimer);
+ // Since the particles have moved spatially update them to correct processors
+ IpplTimings::startTimer(updateTimer);
PL.update(*P, bunchBuffer);
IpplTimings::stopTimer(updateTimer);
// Domain Decomposition
- if (P->balance(totalP, it+1)) {
- msg << "Starting repartition" << endl;
- IpplTimings::startTimer(domainDecomposition);
- P->repartition(FL, mesh, bunchBuffer, isFirstRepartition);
- IpplTimings::stopTimer(domainDecomposition);
- //IpplTimings::startTimer(dumpDataTimer);
- //P->dumpLocalDomains(FL, it+1);
- //IpplTimings::stopTimer(dumpDataTimer);
+ if (P->balance(totalP, it + 1)) {
+ msg << "Starting repartition" << endl;
+ IpplTimings::startTimer(domainDecomposition);
+ P->repartition(FL, mesh, bunchBuffer, isFirstRepartition);
+ IpplTimings::stopTimer(domainDecomposition);
+ // IpplTimings::startTimer(dumpDataTimer);
+ // P->dumpLocalDomains(FL, it+1);
+ // IpplTimings::stopTimer(dumpDataTimer);
}
-
- //scatter the charge onto the underlying grid
- P->scatterCIC(totalP, it+1, hr);
- //Field solve
+ // scatter the charge onto the underlying grid
+ P->scatterCIC(totalP, it + 1, hr);
+
+ // Field solve
IpplTimings::startTimer(SolveTimer);
P->solver_mp->solve();
IpplTimings::stopTimer(SolveTimer);
@@ -423,26 +396,32 @@ int main(int argc, char *argv[]){
// gather E field
P->gatherCIC();
- //kick
+ // kick
IpplTimings::startTimer(PTimer);
auto R2view = P->R.getView();
auto P2view = P->P.getView();
auto E2view = P->E.getView();
- Kokkos::parallel_for("Kick2", P->getLocalNum(),
- KOKKOS_LAMBDA(const size_t j){
- double Eext_x = -(R2view(j)[0] - 0.5*rmax[0]) * (V0/(2*std::pow(rmax[2],2)));
- double Eext_y = -(R2view(j)[1] - 0.5*rmax[1]) * (V0/(2*std::pow(rmax[2],2)));
- double Eext_z = (R2view(j)[2] - 0.5*rmax[2]) * (V0/(std::pow(rmax[2],2)));
-
- E2view(j)[0] += Eext_x;
- E2view(j)[1] += Eext_y;
- E2view(j)[2] += Eext_z;
- P2view(j)[0] = DrInv * ( P2view(j)[0] + alpha * (E2view(j)[0]
- + P2view(j)[1] * Bext + alpha * Bext * E2view(j)[1]) );
- P2view(j)[1] = DrInv * ( P2view(j)[1] + alpha * (E2view(j)[1]
- - P2view(j)[0] * Bext - alpha * Bext * E2view(j)[0]) );
- P2view(j)[2] += alpha * E2view(j)[2];
- });
+ Kokkos::parallel_for(
+ "Kick2", P->getLocalNum(), KOKKOS_LAMBDA(const size_t j) {
+ double Eext_x = -(R2view(j)[0] - 0.5 * rmax[0]) * (V0 / (2 * std::pow(rmax[2], 2)));
+ double Eext_y = -(R2view(j)[1] - 0.5 * rmax[1]) * (V0 / (2 * std::pow(rmax[2], 2)));
+ double Eext_z = (R2view(j)[2] - 0.5 * rmax[2]) * (V0 / (std::pow(rmax[2], 2)));
+
+ E2view(j)[0] += Eext_x;
+ E2view(j)[1] += Eext_y;
+ E2view(j)[2] += Eext_z;
+ P2view(j)[0] =
+ DrInv
+ * (P2view(j)[0]
+ + alpha
+ * (E2view(j)[0] + P2view(j)[1] * Bext + alpha * Bext * E2view(j)[1]));
+ P2view(j)[1] =
+ DrInv
+ * (P2view(j)[1]
+ + alpha
+ * (E2view(j)[1] - P2view(j)[0] * Bext - alpha * Bext * E2view(j)[0]));
+ P2view(j)[2] += alpha * E2view(j)[2];
+ });
IpplTimings::stopTimer(PTimer);
P->time_m += dt;
@@ -450,7 +429,7 @@ int main(int argc, char *argv[]){
P->dumpData();
P->gatherStatistics(totalP);
IpplTimings::stopTimer(dumpDataTimer);
- msg << "Finished time step: " << it+1 << " time: " << P->time_m << endl;
+ msg << "Finished time step: " << it + 1 << " time: " << P->time_m << endl;
}
msg << "Penning Trap: End." << endl;
@@ -459,7 +438,8 @@ int main(int argc, char *argv[]){
IpplTimings::print(std::string("timing.dat"));
auto end = std::chrono::high_resolution_clock::now();
- std::chrono::duration time_chrono = std::chrono::duration_cast>(end - start);
+ std::chrono::duration time_chrono =
+ std::chrono::duration_cast>(end - start);
std::cout << "Elapsed time: " << time_chrono.count() << std::endl;
return 0;
diff --git a/alpine/UniformPlasmaTest.cpp b/alpine/UniformPlasmaTest.cpp
index 62924052d..5383dbc1d 100644
--- a/alpine/UniformPlasmaTest.cpp
+++ b/alpine/UniformPlasmaTest.cpp
@@ -30,13 +30,13 @@
//
#include "ChargedParticles.hpp"
-#include
-#include
+#include
#include
#include
-#include
+#include
+#include
-#include
+#include
#include
#include "Utility/IpplTimings.h"
@@ -45,79 +45,71 @@ const char* TestName = "UniformPlasmaTest";
template
struct generate_random {
+ using view_type = typename ippl::detail::ViewType::view_type;
+ // Output View for the random numbers
+ view_type vals;
- using view_type = typename ippl::detail::ViewType::view_type;
- // Output View for the random numbers
- view_type vals;
+ // The GeneratorPool
+ GeneratorPool rand_pool;
- // The GeneratorPool
- GeneratorPool rand_pool;
+ T start, end;
- T start, end;
+ // Initialize all members
+ generate_random(view_type vals_, GeneratorPool rand_pool_, T start_, T end_)
+ : vals(vals_)
+ , rand_pool(rand_pool_)
+ , start(start_)
+ , end(end_) {}
- // Initialize all members
- generate_random(view_type vals_, GeneratorPool rand_pool_,
- T start_, T end_)
- : vals(vals_), rand_pool(rand_pool_),
- start(start_), end(end_) {}
+ KOKKOS_INLINE_FUNCTION void operator()(const size_t i) const {
+ // Get a random number state from the pool for the active thread
+ typename GeneratorPool::generator_type rand_gen = rand_pool.get_state();
- KOKKOS_INLINE_FUNCTION
- void operator()(const size_t i) const {
- // Get a random number state from the pool for the active thread
- typename GeneratorPool::generator_type rand_gen = rand_pool.get_state();
+ // Draw samples numbers from the pool as double in the range [start, end)
+ for (unsigned d = 0; d < Dim; ++d) {
+ vals(i)[d] = rand_gen.drand(start[d], end[d]);
+ }
- // Draw samples numbers from the pool as double in the range [start, end)
- for (unsigned d = 0; d < Dim; ++d) {
- vals(i)[d] = rand_gen.drand(start[d], end[d]);
+ // Give the state back, which will allow another thread to acquire it
+ rand_pool.free_state(rand_gen);
}
-
- // Give the state back, which will allow another thread to acquire it
- rand_pool.free_state(rand_gen);
- }
};
-int main(int argc, char *argv[]){
+int main(int argc, char* argv[]) {
Ippl ippl(argc, argv);
Inform msg("UniformPlasmaTest");
- Inform msg2all(argv[0],INFORM_ALL_NODES);
+ Inform msg2all(argv[0], INFORM_ALL_NODES);
Ippl::Comm->setDefaultOverallocation(std::atof(argv[8]));
- auto start = std::chrono::high_resolution_clock::now();
- ippl::Vector nr = {
- std::atoi(argv[1]),
- std::atoi(argv[2]),
- std::atoi(argv[3])
- };
-
- static IpplTimings::TimerRef mainTimer = IpplTimings::getTimer("total");
- static IpplTimings::TimerRef particleCreation = IpplTimings::getTimer("particlesCreation");
- static IpplTimings::TimerRef dumpDataTimer = IpplTimings::getTimer("dumpData");
- static IpplTimings::TimerRef PTimer = IpplTimings::getTimer("pushVelocity");
- static IpplTimings::TimerRef temp = IpplTimings::getTimer("randomMove");
- static IpplTimings::TimerRef RTimer = IpplTimings::getTimer("pushPosition");
- static IpplTimings::TimerRef updateTimer = IpplTimings::getTimer("update");
- static IpplTimings::TimerRef DummySolveTimer = IpplTimings::getTimer("solveWarmup");
- static IpplTimings::TimerRef SolveTimer = IpplTimings::getTimer("solve");
+ auto start = std::chrono::high_resolution_clock::now();
+ ippl::Vector nr = {std::atoi(argv[1]), std::atoi(argv[2]), std::atoi(argv[3])};
+
+ static IpplTimings::TimerRef mainTimer = IpplTimings::getTimer("total");
+ static IpplTimings::TimerRef particleCreation = IpplTimings::getTimer("particlesCreation");
+ static IpplTimings::TimerRef dumpDataTimer = IpplTimings::getTimer("dumpData");
+ static IpplTimings::TimerRef PTimer = IpplTimings::getTimer("pushVelocity");
+ static IpplTimings::TimerRef temp = IpplTimings::getTimer("randomMove");
+ static IpplTimings::TimerRef RTimer = IpplTimings::getTimer("pushPosition");
+ static IpplTimings::TimerRef updateTimer = IpplTimings::getTimer("update");
+ static IpplTimings::TimerRef DummySolveTimer = IpplTimings::getTimer("solveWarmup");
+ static IpplTimings::TimerRef SolveTimer = IpplTimings::getTimer("solve");
static IpplTimings::TimerRef domainDecomposition = IpplTimings::getTimer("loadBalance");
IpplTimings::startTimer(mainTimer);
const size_type totalP = std::atoll(argv[4]);
- const unsigned int nt = std::atoi(argv[5]);
+ const unsigned int nt = std::atoi(argv[5]);
- msg << "Uniform Plasma Test"
- << endl
- << "nt " << nt << " Np= "
- << totalP << " grid = " << nr
- << endl;
+ msg << "Uniform Plasma Test" << endl
+ << "nt " << nt << " Np= " << totalP << " grid = " << nr << endl;
using bunch_type = ChargedParticles;
- std::unique_ptr P;
+ std::unique_ptr P;
ippl::NDIndex domain;
- for (unsigned i = 0; i< Dim; i++) {
+ for (unsigned i = 0; i < Dim; i++) {
domain[i] = ippl::Index(nr[i]);
}
@@ -133,24 +125,24 @@ int main(int argc, char *argv[]){
double dy = rmax[1] / nr[1];
double dz = rmax[2] / nr[2];
- Vector_t hr = {dx, dy, dz};
+ Vector_t hr = {dx, dy, dz};
Vector_t origin = {rmin[0], rmin[1], rmin[2]};
const double dt = 1.0;
- const bool isAllPeriodic=true;
+ const bool isAllPeriodic = true;
Mesh_t mesh(domain, hr, origin);
FieldLayout_t FL(domain, decomp, isAllPeriodic);
PLayout_t PL(FL, mesh);
double Q = -1562.5;
- P = std::make_unique(PL,hr,rmin,rmax,decomp,Q);
+ P = std::make_unique(PL, hr, rmin, rmax, decomp, Q);
- P->nr_m = nr;
+ P->nr_m = nr;
size_type nloc = totalP / Ippl::Comm->size();
- int rest = (int) (totalP - nloc * Ippl::Comm->size());
+ int rest = (int)(totalP - nloc * Ippl::Comm->size());
- if ( Ippl::Comm->rank() < rest )
+ if (Ippl::Comm->rank() < rest)
++nloc;
IpplTimings::startTimer(particleCreation);
@@ -158,17 +150,16 @@ int main(int argc, char *argv[]){
const ippl::NDIndex& lDom = FL.getLocalNDIndex();
Vector_t Rmin, Rmax;
- for (unsigned d = 0; d rand_pool64((size_type)(42 + 100*Ippl::Comm->rank()));
- Kokkos::parallel_for(nloc,
- generate_random, Dim>(
- P->R.getView(), rand_pool64, Rmin, Rmax));
+ Kokkos::Random_XorShift64_Pool<> rand_pool64((size_type)(42 + 100 * Ippl::Comm->rank()));
+ Kokkos::parallel_for(nloc, generate_random, Dim>(
+ P->R.getView(), rand_pool64, Rmin, Rmax));
Kokkos::fence();
- P->q = P->Q_m/totalP;
+ P->q = P->Q_m / totalP;
P->P = 0.0;
IpplTimings::stopTimer(particleCreation);
@@ -177,7 +168,7 @@ int main(int argc, char *argv[]){
bunch_type bunchBuffer(PL);
- IpplTimings::startTimer(updateTimer);
+ IpplTimings::startTimer(updateTimer);
PL.update(*P, bunchBuffer);
IpplTimings::stopTimer(updateTimer);
@@ -185,15 +176,14 @@ int main(int argc, char *argv[]){
P->stype_m = argv[6];
P->initSolver();
- P->time_m = 0.0;
+ P->time_m = 0.0;
P->loadbalancefreq_m = std::atoi(argv[7]);
-
+
IpplTimings::startTimer(DummySolveTimer);
P->rho_m = 0.0;
P->solver_mp->solve();
IpplTimings::stopTimer(DummySolveTimer);
-
P->scatterCIC(totalP, 0, hr);
P->initializeORB(FL, mesh);
bool fromAnalyticDensity = false;
@@ -209,12 +199,10 @@ int main(int argc, char *argv[]){
P->gatherStatistics(totalP);
IpplTimings::stopTimer(dumpDataTimer);
-
// begin main timestep loop
msg << "Starting iterations ..." << endl;
- //P->gatherStatistics(totalP);
- for (unsigned int it=0; itgatherStatistics(totalP);
+ for (unsigned int it = 0; it < nt; it++) {
// LeapFrog time stepping https://en.wikipedia.org/wiki/Leapfrog_integration
// Here, we assume a constant charge-to-mass ratio of -1 for
// all the particles hence eliminating the need to store mass as
@@ -226,35 +214,33 @@ int main(int argc, char *argv[]){
IpplTimings::startTimer(temp);
Kokkos::parallel_for(P->getLocalNum(),
- generate_random, Dim>(
- P->P.getView(), rand_pool64, -hr, hr));
+ generate_random, Dim>(
+ P->P.getView(), rand_pool64, -hr, hr));
Kokkos::fence();
IpplTimings::stopTimer(temp);
- //drift
+ // drift
IpplTimings::startTimer(RTimer);
P->R = P->R + dt * P->P;
IpplTimings::stopTimer(RTimer);
- //Since the particles have moved spatially update them to correct processors
- IpplTimings::startTimer(updateTimer);
+ // Since the particles have moved spatially update them to correct processors
+ IpplTimings::startTimer(updateTimer);
PL.update(*P, bunchBuffer);
IpplTimings::stopTimer(updateTimer);
-
// Domain Decomposition
- if (P->balance(totalP, it+1)) {
- msg << "Starting repartition" << endl;
- IpplTimings::startTimer(domainDecomposition);
- P->repartition(FL, mesh, bunchBuffer, fromAnalyticDensity);
- IpplTimings::stopTimer(domainDecomposition);
+ if (P->balance(totalP, it + 1)) {
+ msg << "Starting repartition" << endl;
+ IpplTimings::startTimer(domainDecomposition);
+ P->repartition(FL, mesh, bunchBuffer, fromAnalyticDensity);
+ IpplTimings::stopTimer(domainDecomposition);
}
- //scatter the charge onto the underlying grid
- P->scatterCIC(totalP, it+1, hr);
+ // scatter the charge onto the underlying grid
+ P->scatterCIC(totalP, it + 1, hr);
- //Field solve
+ // Field solve
IpplTimings::startTimer(SolveTimer);
P->solver_mp->solve();
IpplTimings::stopTimer(SolveTimer);
@@ -262,7 +248,7 @@ int main(int argc, char *argv[]){
// gather E field
P->gatherCIC();
- //kick
+ // kick
IpplTimings::startTimer(PTimer);
P->P = P->P - 0.5 * dt * P->E;
IpplTimings::stopTimer(PTimer);
@@ -272,7 +258,7 @@ int main(int argc, char *argv[]){
P->dumpData();
P->gatherStatistics(totalP);
IpplTimings::stopTimer(dumpDataTimer);
- msg << "Finished time step: " << it+1 << " time: " << P->time_m << endl;
+ msg << "Finished time step: " << it + 1 << " time: " << P->time_m << endl;
}
msg << "Uniform Plasma Test: End." << endl;
@@ -281,7 +267,8 @@ int main(int argc, char *argv[]){
IpplTimings::print(std::string("timing.dat"));
auto end = std::chrono::high_resolution_clock::now();
- std::chrono::duration time_chrono = std::chrono::duration_cast>(end - start);
+ std::chrono::duration time_chrono =
+ std::chrono::duration_cast>(end - start);
std::cout << "Elapsed time: " << time_chrono.count() << std::endl;
return 0;
diff --git a/hooks/create-hook-symlink.sh b/hooks/create-hook-symlink.sh
new file mode 100755
index 000000000..038393034
--- /dev/null
+++ b/hooks/create-hook-symlink.sh
@@ -0,0 +1,18 @@
+#!/bin/bash
+# OPAL utility script for setting up git hooks (with modifications)
+
+HOOK_NAMES="applypatch-msg pre-applypatch post-applypatch pre-commit prepare-commit-msg commit-msg post-commit pre-rebase post-checkout post-merge pre-receive update post-receive post-update pre-auto-gc"
+
+GITBASEDIR=`git rev-parse --show-toplevel`
+HOOK_DIR=$GITBASEDIR/.git/hooks
+
+for hook in $HOOK_NAMES; do
+ # If the hook already exists, is executable, and is not a symlink
+ if [ ! -h $HOOK_DIR/$hook -a -x $HOOK_DIR/$hook ]; then
+ mv $HOOK_DIR/$hook $HOOK_DIR/$hook.local
+ fi
+ # if the hook is defined in hooks/, then symlink it in the git hooks
+ if [ -f $GITBASEDIR/hooks/$hook ]; then
+ ln -sf $GITBASEDIR/hooks/$hook $HOOK_DIR/$hook
+ fi
+done
diff --git a/hooks/pre-commit b/hooks/pre-commit
new file mode 100755
index 000000000..4483e750b
--- /dev/null
+++ b/hooks/pre-commit
@@ -0,0 +1,5 @@
+#!/bin/sh
+git diff --cached --name-only | grep -e '\.h$' -e '\.cpp$' -e '\.hpp$' | while read FILE; do
+ clang-format -i "$FILE"
+ git add "$FILE"
+done
diff --git a/src/AmrParticle/AmrParticleBase.h b/src/AmrParticle/AmrParticleBase.h
index 3908b86c5..79dd447b4 100644
--- a/src/AmrParticle/AmrParticleBase.h
+++ b/src/AmrParticle/AmrParticleBase.h
@@ -17,8 +17,8 @@
// const ParticleAttrib >& pp,
// int lbase = 0, int lfine = -1) const;
//
-// Copyright (c) 2016 - 2020, Matthias Frey, Uldis Locans, Paul Scherrer Institut, Villigen PSI, Switzerland
-// All rights reserved
+// Copyright (c) 2016 - 2020, Matthias Frey, Uldis Locans, Paul Scherrer Institut, Villigen PSI,
+// Switzerland All rights reserved
//
// Implemented as part of the PhD thesis
// "Precise Simulations of Multibunches in High Intensity Cyclotrons"
@@ -42,47 +42,45 @@
#include "AmrParticleLevelCounter.h"
-template
+template
class AmrParticleBase : public IpplParticleBase {
-
public:
- typedef typename PLayout::ParticlePos_t ParticlePos_t;
- typedef typename PLayout::ParticleIndex_t ParticleIndex_t;
- typedef typename PLayout::SingleParticlePos_t SingleParticlePos_t;
- typedef typename PLayout::AmrField_t AmrField_t;
- typedef typename PLayout::AmrVectorField_t AmrVectorField_t;
- typedef typename PLayout::AmrScalarFieldContainer_t AmrScalarFieldContainer_t;
- typedef typename PLayout::AmrVectorFieldContainer_t AmrVectorFieldContainer_t;
-
- typedef long SortListIndex_t;
- typedef std::vector SortList_t;
- typedef std::vector attrib_container_t;
-
- ParticleIndex_t Level; // m_lev
- ParticleIndex_t Grid; // m_grid
-
+ typedef typename PLayout::ParticlePos_t ParticlePos_t;
+ typedef typename PLayout::ParticleIndex_t ParticleIndex_t;
+ typedef typename PLayout::SingleParticlePos_t SingleParticlePos_t;
+ typedef typename PLayout::AmrField_t AmrField_t;
+ typedef typename PLayout::AmrVectorField_t AmrVectorField_t;
+ typedef typename PLayout::AmrScalarFieldContainer_t AmrScalarFieldContainer_t;
+ typedef typename PLayout::AmrVectorFieldContainer_t AmrVectorFieldContainer_t;
+
+ typedef long SortListIndex_t;
+ typedef std::vector SortList_t;
+ typedef std::vector attrib_container_t;
+
+ ParticleIndex_t Level; // m_lev
+ ParticleIndex_t Grid; // m_grid
+
typedef AmrParticleLevelCounter ParticleLevelCounter_t;
-
+
public:
-
AmrParticleBase();
-
+
AmrParticleBase(PLayout* layout);
-
+
~AmrParticleBase() {}
-
- //initialize AmrParticleBase class - add level and grid variables to attribute list
+
+ // initialize AmrParticleBase class - add level and grid variables to attribute list
void initializeAmr() {
this->addAttribute(Level);
this->addAttribute(Grid);
}
-
+
const ParticleLevelCounter_t& getLocalNumPerLevel() const;
-
+
ParticleLevelCounter_t& getLocalNumPerLevel();
-
+
void setLocalNumPerLevel(const ParticleLevelCounter_t& LocalNumPerLevel);
-
+
/* Functions of IpplParticleBase adpated to
* work with AmrParticleLevelCounter:
* - createWithID()
@@ -90,19 +88,19 @@ class AmrParticleBase : public IpplParticleBase {
* - destroy()
* - performDestroy()
*/
-
+
void createWithID(unsigned id);
void create(size_t M);
void destroy(size_t M, size_t I, bool doNow = false);
-
+
void performDestroy(bool updateLocalNum = false);
-
+
// Update the particle object after a timestep. This routine will change
// our local, total, create particle counts properly.
void update();
-
+
/*!
* There's is NO check performed if lev_min <= lev_max and
* lev_min >= 0.
@@ -111,36 +109,36 @@ class AmrParticleBase : public IpplParticleBase {
* @param isRegrid is true if we are updating the grids (default: false)
*/
void update(int lev_min, int lev_max, bool isRegrid = false);
-
+
// Update the particle object after a timestep. This routine will change
// our local, total, create particle counts properly.
void update(const ParticleAttrib& canSwap);
-
+
// sort particles based on the grid and level that they belong to
void sort();
-
+
// sort the particles given a sortlist
- void sort(SortList_t &sortlist);
-
+ void sort(SortList_t& sortlist);
+
PLayout& getAmrLayout() { return this->getLayout(); }
const PLayout& getAmrLayout() const { return this->getLayout(); }
-
+
/*!
* This method is used in the AmrPartBunch::boundp() function
* in order to avoid multpile particle mappings during the
* mesh regridding process.
- *
+ *
* @param forbidTransform true if we don't want to map particles onto
* \f$[-1, 1]^3\f$
*/
inline void setForbidTransform(bool forbidTransform);
-
+
/*!
* @returns true if we are not mapping the particles onto
* \f$[-1, 1]^3\f$ during an update call.
*/
inline bool isForbidTransform() const;
-
+
/*!
* Linear mapping to AMReX computation domain [-1, 1]^3 including the Lorentz
* transform. All dimensions
@@ -151,45 +149,43 @@ class AmrParticleBase : public IpplParticleBase {
* @returns scaling factor
*/
const double& domainMapping(bool inverse = false);
-
+
/*!
* This function is used during the cell tagging routines.
* @returns the scaling factor of the particle domain mapping.
*/
inline const double& getScalingFactor() const;
-
-
+
void setLorentzFactor(const Vector_t& lorentzFactor);
-
-// void lorentzTransform(bool inverse = false);
-
+
+ // void lorentzTransform(bool inverse = false);
+
private:
- void getLocalBounds_m(Vector_t &rmin, Vector_t &rmax);
- void getGlobalBounds_m(Vector_t &rmin, Vector_t &rmax);
-
+ void getLocalBounds_m(Vector_t& rmin, Vector_t& rmax);
+ void getGlobalBounds_m(Vector_t& rmin, Vector_t& rmax);
+
protected:
IpplTimings::TimerRef updateParticlesTimer_m;
IpplTimings::TimerRef sortParticlesTimer_m;
IpplTimings::TimerRef domainMappingTimer_m;
-
- bool forbidTransform_m; ///< To avoid multiple transformations during regrid
-
+
+ bool forbidTransform_m; ///< To avoid multiple transformations during regrid
+
/*!
* Scaling factor for particle coordinate transform
* (used for Poisson solve and particle-to-core distribution)
*/
double scale_m;
-
+
/*!
* Lorentz factor used for the domain mapping.
* Is updated in AmrBoxLib
- *
+ *
*/
Vector_t lorentzFactor_m;
-
-// bool isLorentzTransformed_m;
-
-
+
+ // bool isLorentzTransformed_m;
+
private:
ParticleLevelCounter_t LocalNumPerLevel_m;
};
diff --git a/src/AmrParticle/AmrParticleBase.hpp b/src/AmrParticle/AmrParticleBase.hpp
index 2b7e8c6fb..d964502f9 100644
--- a/src/AmrParticle/AmrParticleBase.hpp
+++ b/src/AmrParticle/AmrParticleBase.hpp
@@ -17,8 +17,8 @@
// const ParticleAttrib >& pp,
// int lbase = 0, int lfine = -1) const;
//
-// Copyright (c) 2016 - 2020, Matthias Frey, Uldis Locans, Paul Scherrer Institut, Villigen PSI, Switzerland
-// All rights reserved
+// Copyright (c) 2016 - 2020, Matthias Frey, Uldis Locans, Paul Scherrer Institut, Villigen PSI,
+// Switzerland All rights reserved
//
// Implemented as part of the PhD thesis
// "Precise Simulations of Multibunches in High Intensity Cyclotrons"
@@ -36,93 +36,81 @@
#ifndef AMR_PARTICLE_BASE_HPP
#define AMR_PARTICLE_BASE_HPP
-#include
#include
+#include
-template
-AmrParticleBase::AmrParticleBase() : forbidTransform_m(false),
- scale_m(1.0),
- lorentzFactor_m(1.0, 1.0, 1.0),
-// isLorentzTransformed_m(false),
- LocalNumPerLevel_m()
-{
+template
+AmrParticleBase::AmrParticleBase()
+ : forbidTransform_m(false)
+ , scale_m(1.0)
+ , lorentzFactor_m(1.0, 1.0, 1.0)
+ ,
+ // isLorentzTransformed_m(false),
+ LocalNumPerLevel_m() {
updateParticlesTimer_m = IpplTimings::getTimer("AMR update particles");
- sortParticlesTimer_m = IpplTimings::getTimer("AMR sort particles");
- domainMappingTimer_m = IpplTimings::getTimer("AMR map particles");
+ sortParticlesTimer_m = IpplTimings::getTimer("AMR sort particles");
+ domainMappingTimer_m = IpplTimings::getTimer("AMR map particles");
}
-
-template
+template
AmrParticleBase::AmrParticleBase(PLayout* layout)
- : IpplParticleBase(layout),
- forbidTransform_m(false),
- scale_m(1.0),
- lorentzFactor_m(1.0, 1.0, 1.0),
-// isLorentzTransformed_m(false),
- LocalNumPerLevel_m()
-{
+ : IpplParticleBase(layout)
+ , forbidTransform_m(false)
+ , scale_m(1.0)
+ , lorentzFactor_m(1.0, 1.0, 1.0)
+ ,
+ // isLorentzTransformed_m(false),
+ LocalNumPerLevel_m() {
updateParticlesTimer_m = IpplTimings::getTimer("AMR update particles");
- sortParticlesTimer_m = IpplTimings::getTimer("AMR sort particles");
- domainMappingTimer_m = IpplTimings::getTimer("AMR map particles");
+ sortParticlesTimer_m = IpplTimings::getTimer("AMR sort particles");
+ domainMappingTimer_m = IpplTimings::getTimer("AMR map particles");
}
-
-template
+template
const typename AmrParticleBase::ParticleLevelCounter_t&
- AmrParticleBase::getLocalNumPerLevel() const
-{
+AmrParticleBase::getLocalNumPerLevel() const {
return LocalNumPerLevel_m;
}
-
-template
+template
typename AmrParticleBase::ParticleLevelCounter_t&
- AmrParticleBase::getLocalNumPerLevel()
-{
+AmrParticleBase::getLocalNumPerLevel() {
return LocalNumPerLevel_m;
}
-
-template
-void AmrParticleBase::setLocalNumPerLevel(
- const ParticleLevelCounter_t& LocalNumPerLevel)
-{
+template
+void AmrParticleBase