diff --git a/src/PartBunch/CMakeLists.txt b/src/PartBunch/CMakeLists.txt
index 9c9ec923fecb2316911022e4bb5e5113617f7322..6b3d5af74e7714d92957765575b3476781265581 100644
--- a/src/PartBunch/CMakeLists.txt
+++ b/src/PartBunch/CMakeLists.txt
@@ -1,5 +1,6 @@
set (_SRCS
PartBunch.cpp
+ FieldSolver.cpp
)
include_directories (
diff --git a/src/PartBunch/FieldSolver.cpp b/src/PartBunch/FieldSolver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ac359971108d4fe78f7ff14f4f208ddcf8d8a52b
--- /dev/null
+++ b/src/PartBunch/FieldSolver.cpp
@@ -0,0 +1,145 @@
+#include "PartBunch/FieldSolver.hpp"
+
+template <>
+void FieldSolver<double,3>::initOpenSolver() {
+ // if constexpr (Dim == 3) {
+ ippl::ParameterList sp;
+ sp.add("output_type", OpenSolver_t<double, 3>::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("algorithm", OpenSolver_t<double, 3>::HOCKNEY);
+ initSolverWithParams<OpenSolver_t<double, 3>>(sp);
+ // } else {
+ //throw std::runtime_error("Unsupported dimensionality for OPEN solver");
+ //}
+}
+
+template <>
+void FieldSolver<double,3>::initSolver() {
+ Inform m;
+ if (this->getStype() == "FFT") {
+ initFFTSolver();
+ } else if (this->getStype() == "CG") {
+ initCGSolver();
+ } else if (this->getStype() == "P3M") {
+ initP3MSolver();
+ } else if (this->getStype() == "FFTOPEN") {
+ initFFTSolver();
+ } else if (this->getStype() == "NONE") {
+ initNullSolver();
+ }
+ else {
+ m << "No solver matches the argument: " << this->getStype() << endl;
+ }
+}
+
+template <>
+void FieldSolver<double,3>::setPotentialBCs() {
+ // CG requires explicit periodic boundary conditions while the periodic Poisson solver
+ // simply assumes them
+ if (this->getStype() == "CG") {
+ typedef ippl::BConds<Field<double, 3>, 3> bc_type;
+ bc_type allPeriodic;
+ for (unsigned int i = 0; i < 2 * 3; ++i) {
+ allPeriodic[i] = std::make_shared<ippl::PeriodicFace<Field<double, 3>>>(i);
+ }
+ phi_m->setFieldBC(allPeriodic);
+ }
+ }
+
+template<>
+void FieldSolver<double,3>::runSolver() {
+ Inform m("runSolver ");
+ constexpr int Dim = 3;
+
+ m << "solver type= " << this->getStype() << endl;
+
+ if (this->getStype() == "CG") {
+ CGSolver_t<double, 3>& solver = std::get<CGSolver_t<double, 3>>(this->getSolver());
+ solver.solve();
+
+ if (ippl::Comm->rank() == 0) {
+ std::stringstream fname;
+ fname << "data/CG_";
+ fname << ippl::Comm->size();
+ fname << ".csv";
+
+ Inform log(NULL, fname.str().c_str(), Inform::APPEND);
+ int iterations = solver.getIterationCount();
+ // Assume the dummy solve is the first call
+ if (iterations == 0) {
+ log << "residue,iterations" << endl;
+ }
+ // Don't print the dummy solve
+ if (iterations > 0) {
+ log << solver.getResidue() << "," << iterations << endl;
+ }
+ }
+ ippl::Comm->barrier();
+ } else if (this->getStype() == "FFT") {
+ if constexpr (Dim == 2 || Dim == 3) {
+ std::get<FFTSolver_t<double, 3>>(this->getSolver()).solve();
+ }
+ } else if (this->getStype() == "P3M") {
+ if constexpr (Dim == 3) {
+ std::get<P3MSolver_t<double, 3>>(this->getSolver()).solve();
+ }
+ } else if (this->getStype() == "FFTOPEN") {
+ if constexpr (Dim == 3) {
+ m << " about to open.solve " << endl;
+ std::get<FFTSolver_t<double, 3>>(this->getSolver()).solve();
+ m << " done open.solve " << endl;
+ }
+ } else {
+ throw std::runtime_error("Unknown solver type");
+ }
+
+ m << " done xx.solve() " << endl;
+}
+
+/*
+template<>
+void FieldSolver<double,3>::initNullSolver() {
+ // if constexpr (Dim == 2 || Dim == 3) {
+ //ippl::ParameterList sp;
+ throw std::runtime_error("Not implemented Null solver");
+ // } else {
+ //throw std::runtime_error("Unsupported dimensionality for Null solver");
+ //}
+}
+*/
+
+/*
+template <>
+void FieldSolver<double,3>::initCGSolver() {
+ ippl::ParameterList sp;
+ sp.add("output_type", CGSolver_t<double, 3>::GRAD);
+ // Increase tolerance in the 1D case
+ sp.add("tolerance", 1e-10);
+
+ initSolverWithParams<CGSolver_t<double, 3>>(sp);
+}
+
+template<>
+void FieldSolver<double,3>::initP3MSolver() {
+ // if constexpr (Dim == 3) {
+ ippl::ParameterList sp;
+ sp.add("output_type", P3MSolver_t<double, 3>::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);
+
+ initSolverWithParams<P3MSolver_t<double, 3>>(sp);
+ // } else {
+ // throw std::runtime_error("Unsupported dimensionality for P3M solver");
+ // }
+}
+
+*/
\ No newline at end of file
diff --git a/src/PartBunch/FieldSolver.hpp b/src/PartBunch/FieldSolver.hpp
index 2ff8dfa4b3d123adc460568a686db569085ccb71..7573fc2f296a188707bd9ec90d9ebe62ee364dc9 100644
--- a/src/PartBunch/FieldSolver.hpp
+++ b/src/PartBunch/FieldSolver.hpp
@@ -70,77 +70,14 @@ public:
phi_m = phi;
}
- void initSolver() override {
- Inform m("solver ");
- if (this->getStype() == "FFT") {
- initFFTSolver();
- } else if (this->getStype() == "CG") {
- initCGSolver();
- } else if (this->getStype() == "P3M") {
- initP3MSolver();
- } else if (this->getStype() == "OPEN") {
- initOpenSolver();
- } else if (this->getStype() == "NONE") {
- initNullSolver();
- }
- else {
- m << "No solver matches the argument: " << this->getStype() << endl;
- }
- }
+ void initOpenSolver();
- void setPotentialBCs() {
- // CG requires explicit periodic boundary conditions while the periodic Poisson solver
- // simply assumes them
- if (this->getStype() == "CG") {
- typedef ippl::BConds<Field<T, Dim>, Dim> bc_type;
- bc_type allPeriodic;
- for (unsigned int i = 0; i < 2 * Dim; ++i) {
- allPeriodic[i] = std::make_shared<ippl::PeriodicFace<Field<T, Dim>>>(i);
- }
- phi_m->setFieldBC(allPeriodic);
- }
- }
+ void initSolver() override ;
- void runSolver() override {
- if (this->getStype() == "CG") {
- CGSolver_t<T, Dim>& solver = std::get<CGSolver_t<T, Dim>>(this->getSolver());
- solver.solve();
-
- if (ippl::Comm->rank() == 0) {
- std::stringstream fname;
- fname << "data/CG_";
- fname << ippl::Comm->size();
- fname << ".csv";
-
- Inform log(NULL, fname.str().c_str(), Inform::APPEND);
- int iterations = solver.getIterationCount();
- // Assume the dummy solve is the first call
- if (iterations == 0) {
- log << "residue,iterations" << endl;
- }
- // Don't print the dummy solve
- if (iterations > 0) {
- log << solver.getResidue() << "," << iterations << endl;
- }
- }
- ippl::Comm->barrier();
- } else if (this->getStype() == "FFT") {
- if constexpr (Dim == 2 || Dim == 3) {
- std::get<FFTSolver_t<T, Dim>>(this->getSolver()).solve();
- }
- } else if (this->getStype() == "P3M") {
- if constexpr (Dim == 3) {
- std::get<P3MSolver_t<T, Dim>>(this->getSolver()).solve();
- }
- } else if (this->getStype() == "OPEN") {
- if constexpr (Dim == 3) {
- std::get<OpenSolver_t<T, Dim>>(this->getSolver()).solve();
- }
- } else {
- throw std::runtime_error("Unknown solver type");
- }
- }
+ void setPotentialBCs();
+ void runSolver() override;
+
template <typename Solver>
void initSolverWithParams(const ippl::ParameterList& sp) {
this->getSolver().template emplace<Solver>();
@@ -162,74 +99,23 @@ public:
}
}
- void initNullSolver() {
- if constexpr (Dim == 2 || Dim == 3) {
- ippl::ParameterList sp;
- throw std::runtime_error("Not implemented Null solver");
- } else {
- throw std::runtime_error("Unsupported dimensionality for Null solver");
- }
- }
-
+ void initNullSolver() { }
+
void initFFTSolver() {
- if constexpr (Dim == 2 || Dim == 3) {
- ippl::ParameterList sp;
- sp.add("output_type", FFTSolver_t<T, Dim>::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);
-
- initSolverWithParams<FFTSolver_t<T, Dim>>(sp);
- } else {
- throw std::runtime_error("Unsupported dimensionality for FFT solver");
- }
+ ippl::ParameterList sp;
+ sp.add("output_type", FFTSolver_t<double, 3>::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);
+ initSolverWithParams<FFTSolver_t<double, 3>>(sp);
}
+
+ void initCGSolver() { }
- void initCGSolver() {
- ippl::ParameterList sp;
- sp.add("output_type", CGSolver_t<T, Dim>::GRAD);
- // Increase tolerance in the 1D case
- sp.add("tolerance", 1e-10);
-
- initSolverWithParams<CGSolver_t<T, Dim>>(sp);
- }
+ void initP3MSolver() { }
- void initP3MSolver() {
- if constexpr (Dim == 3) {
- ippl::ParameterList sp;
- sp.add("output_type", P3MSolver_t<T, Dim>::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);
-
- initSolverWithParams<P3MSolver_t<T, Dim>>(sp);
- } else {
- throw std::runtime_error("Unsupported dimensionality for P3M solver");
- }
- }
-
- void initOpenSolver() {
- if constexpr (Dim == 3) {
- ippl::ParameterList sp;
- sp.add("output_type", OpenSolver_t<T, Dim>::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("algorithm", OpenSolver_t<T, Dim>::HOCKNEY);
-
- initSolverWithParams<OpenSolver_t<T, Dim>>(sp);
- } else {
- throw std::runtime_error("Unsupported dimensionality for OPEN solver");
- }
- }
};
#endif
diff --git a/src/PartBunch/PartBunch.cpp b/src/PartBunch/PartBunch.cpp
index 938fd2d8d7c14c39f4e29078212a482b2f3cb972..3229e3b986760a7c833924af0c02477e9fbebac3 100644
--- a/src/PartBunch/PartBunch.cpp
+++ b/src/PartBunch/PartBunch.cpp
@@ -2,6 +2,89 @@
#include <boost/numeric/ublas/io.hpp>
#include "Utilities/Util.h"
+template<>
+void PartBunch<double,3>::setSolver(std::string solver) {
+ Inform m("setSolver ");
+
+ if (this->solver_m != "")
+ m << "Warning solver already initiated but overwrite ..." << endl;
+
+ this->solver_m = solver;
+
+ this->fcontainer_m->initializeFields(this->solver_m);
+
+ this->setFieldSolver(std::make_shared<FieldSolver_t>(
+ this->solver_m, &this->fcontainer_m->getRho(), &this->fcontainer_m->getE(),
+ &this->fcontainer_m->getPhi()));
+
+ this->fsolver_m->initSolver();
+
+
+ m << "Solver initialized" << endl;
+
+ /// ADA we need to be able to set a load balancer when not having a field solver
+ this->setLoadBalancer(std::make_shared<LoadBalancer_t>(
+ this->lbt_m, this->fcontainer_m, this->pcontainer_m, this->fsolver_m));
+
+ m << "Load Balancer set" << endl;
+}
+
+template<>
+void PartBunch<double,3>::spaceChargeEFieldCheck() {
+
+ Inform msg("EParticleStats");
+
+ auto pE_view = this->pcontainer_m->E.getView();
+
+ double avgE = 0.0;
+ double minEComponent = std::numeric_limits<double>::max();
+ double maxEComponent = std::numeric_limits<double>::min();
+ double minE = std::numeric_limits<double>::max();
+ double maxE = std::numeric_limits<double>::min();
+ double cc = getCouplingConstant();
+
+ int myRank = ippl::Comm->rank();
+ Kokkos::parallel_reduce(
+ "check e-field", this->getLocalNum(),
+ KOKKOS_LAMBDA(const int i, double& loc_avgE, double& loc_minEComponent,
+ double& loc_maxEComponent, double& loc_minE, double& loc_maxE) {
+ double EX = pE_view[i][0]*cc;
+ double EY = pE_view[i][1]*cc;
+ double EZ = pE_view[i][2]*cc;
+
+ double ENorm = Kokkos::sqrt(EX*EX + EY*EY + EZ*EZ);
+
+ loc_avgE += ENorm;
+
+ loc_minEComponent = EX < loc_minEComponent ? EX : loc_minEComponent;
+ loc_minEComponent = EY < loc_minEComponent ? EY : loc_minEComponent;
+ loc_minEComponent = EZ < loc_minEComponent ? EZ : loc_minEComponent;
+
+ loc_maxEComponent = EX > loc_maxEComponent ? EX : loc_maxEComponent;
+ loc_maxEComponent = EY > loc_maxEComponent ? EY : loc_maxEComponent;
+ loc_maxEComponent = EZ > loc_maxEComponent ? EZ : loc_maxEComponent;
+
+ loc_minE = ENorm < loc_minE ? ENorm : loc_minE;
+ loc_maxE = ENorm > loc_maxE ? ENorm : loc_maxE;
+ },
+ Kokkos::Sum<double>(avgE), Kokkos::Min<double>(minEComponent),
+ Kokkos::Max<double>(maxEComponent), Kokkos::Min<double>(minE),
+ Kokkos::Max<double>(maxE));
+
+ MPI_Reduce(myRank == 0 ? MPI_IN_PLACE : &avgE, &avgE, 1, MPI_DOUBLE, MPI_SUM, 0, ippl::Comm->getCommunicator());
+ MPI_Reduce(myRank == 0 ? MPI_IN_PLACE : &minEComponent, &minEComponent, 1, MPI_DOUBLE, MPI_MIN, 0, ippl::Comm->getCommunicator());
+ MPI_Reduce(myRank == 0 ? MPI_IN_PLACE : &maxEComponent, &maxEComponent, 1, MPI_DOUBLE, MPI_MAX, 0, ippl::Comm->getCommunicator());
+ MPI_Reduce(myRank == 0 ? MPI_IN_PLACE : &minE, &minE, 1, MPI_DOUBLE, MPI_MIN, 0, ippl::Comm->getCommunicator());
+ MPI_Reduce(myRank == 0 ? MPI_IN_PLACE : &maxE, &maxE, 1, MPI_DOUBLE, MPI_MAX, 0, ippl::Comm->getCommunicator());
+
+ avgE /= this->getTotalNum();
+ msg << "avgENorm = " << avgE << endl;
+ msg << "minEComponent = " << minEComponent << endl;
+ msg << "maxEComponent = " << maxEComponent << endl;
+ msg << "minE = " << minE << endl;
+ msg << "maxE = " << maxE << endl;
+}
+
template<>
void PartBunch<double,3>::calcBeamParameters() {
// Usefull constants
@@ -97,11 +180,24 @@ void PartBunch<double,3>::calcBeamParameters() {
rmax_m(i) = rmax[i];
rmin_m(i) = rmin[i];
}
-
ippl::Comm->barrier();
+ }
+template<>
+void PartBunch<double,3>::pre_run() {
+ Inform m("pre_run ");
+
+ static IpplTimings::TimerRef DummySolveTimer = IpplTimings::getTimer("solveWarmup");
+ IpplTimings::startTimer(DummySolveTimer);
- }
+ m << " init rho" << endl;
+ this->fcontainer_m->getRho() = 0.0;
+ m << " init rho done " << endl;
+
+ this->fsolver_m->runSolver();
+ m << " solve done " << endl;
+ IpplTimings::stopTimer(DummySolveTimer);
+}
template<>
Inform& PartBunch<double,3>::print(Inform& os) {
@@ -285,9 +381,7 @@ void PartBunch<double,3>::computeSelfFields() {
//divide charge by a 'grid-cube' volume to get [C/m^3]
rho = rho *tmp2;
- double Npoints = nr_m[0] * nr_m[1] * nr_m[2];
- // auto rhoView = rho.getView();
localDensity2 = 0.;
ippl::parallel_reduce(
"Density stats", ippl::getRangePolicy(rhoView),
@@ -297,18 +391,22 @@ void PartBunch<double,3>::computeSelfFields() {
},
Kokkos::Sum<double>(localDensity2) );
ippl::Comm->reduce(localDensity2, Density2, 1, std::plus<double>());
-
+
+ double Npoints = nr_m[0] * nr_m[1] * nr_m[2];
rmsDensity_m = std::sqrt( (1.0 /Npoints) * Density2 / Physics::q_e / Physics::q_e );
this->calcDebyeLength();
m << "gammaz = " << gammaz << endl;
m << "hr_scaled = " << hr_scaled << endl;
-
+ m << "coupling constant= " << getCouplingConstant() << endl;
+
this->fsolver_m->runSolver();
gather(this->pcontainer_m->E, this->fcontainer_m->getE(), this->pcontainer_m->R);
-
+
+ spaceChargeEFieldCheck();
+
IpplTimings::stopTimer(SolveTimer);
}
diff --git a/src/PartBunch/PartBunch.hpp b/src/PartBunch/PartBunch.hpp
index 75a339f6b1e3032d5416bc19de3463d19464c1bb..c5a29e80de3365feec7ed7744826b4ffcd39215c 100644
--- a/src/PartBunch/PartBunch.hpp
+++ b/src/PartBunch/PartBunch.hpp
@@ -62,6 +62,14 @@ diagnostics(): calculate statistics and maybe write tp h5 and stat files
#include "Algorithms/PartData.h"
+
+template <typename T>
+KOKKOS_INLINE_FUNCTION typename T::value_type L2Norm(T& x) {
+ return sqrt(dot(x, x).apply());
+}
+
+
+
using view_type = typename ippl::detail::ViewType<ippl::Vector<double, 3>, 1>::view_type;
@@ -276,11 +284,6 @@ public:
IpplTimings::startTimer(prerun);
pre_run();
IpplTimings::stopTimer(prerun);
-
- /*
- end wrmup
- */
-
}
void bunchUpdate();
@@ -294,39 +297,10 @@ public:
return this->pcontainer_m;
}
- void setSolver(std::string solver) {
- Inform m("setSolver ");
-
- if (this->solver_m != "")
- m << "Warning solver already initiated but overwrite ..." << endl;
-
- this->solver_m = solver;
-
- this->fcontainer_m->initializeFields(this->solver_m);
-
- this->setFieldSolver(std::make_shared<FieldSolver_t>(
- this->solver_m, &this->fcontainer_m->getRho(), &this->fcontainer_m->getE(),
- &this->fcontainer_m->getPhi()));
-
- this->fsolver_m->initSolver();
-
- /// ADA we need to be able to set a load balancer when not having a field solver
- this->setLoadBalancer(std::make_shared<LoadBalancer_t>(
- this->lbt_m, this->fcontainer_m, this->pcontainer_m, this->fsolver_m));
- }
-
- void pre_run() override {
- static IpplTimings::TimerRef DummySolveTimer = IpplTimings::getTimer("solveWarmup");
- IpplTimings::startTimer(DummySolveTimer);
-
- this->fcontainer_m->getRho() = 0.0;
-
- this->fsolver_m->runSolver();
-
- IpplTimings::stopTimer(DummySolveTimer);
-
- }
+ void setSolver(std::string solver);
+ void pre_run() override ;
+
public:
void updateMoments(){
this->pcontainer_m->updateMoments();
@@ -383,10 +357,11 @@ public:
*/
double getCouplingConstant() const {
- return 1.0;
+ return couplingConstant_m;
}
-
+
void setCouplingConstant(double c) {
+ couplingConstant_m = c;
}
void calcBeamParameters();
@@ -850,6 +825,10 @@ public:
double calculateAngle(double x, double y) {
return 0.0;
}
+
+ // Sanity check functions
+ void spaceChargeEFieldCheck();
+
};
/* \todo