Code indexing in gitaly is broken and leads to code not being visible to the user. We work on the issue with highest priority.

Skip to content
Snippets Groups Projects
Commit 7e9af2f2 authored by adelmann's avatar adelmann :reminder_ribbon:
Browse files

add some diag. Field solver is wrong: FFT instead of FFTOpen however: FFTOpen crashes

parent a64d5d2a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/PartBunch/CMakeLists.txt
+ 1
− 0
View file @ 7e9af2f2
set (_SRCS set (_SRCS
PartBunch.cpp PartBunch.cpp
FieldSolver.cpp
) )
include_directories ( include_directories (
......
src/PartBunch/FieldSolver.cpp 0 → 100644
+ 145
− 0
View file @ 7e9af2f2
#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
src/PartBunch/FieldSolver.hpp
+ 19
− 133
View file @ 7e9af2f2
...@@ -70,77 +70,14 @@ public: ...@@ -70,77 +70,14 @@ public:
phi_m = phi; phi_m = phi;
} }
void initSolver() override { void initOpenSolver();
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 setPotentialBCs() { void initSolver() override ;
// 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 runSolver() override { void setPotentialBCs();
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 runSolver() override;
template <typename Solver> template <typename Solver>
void initSolverWithParams(const ippl::ParameterList& sp) { void initSolverWithParams(const ippl::ParameterList& sp) {
this->getSolver().template emplace<Solver>(); this->getSolver().template emplace<Solver>();
...@@ -162,74 +99,23 @@ public: ...@@ -162,74 +99,23 @@ public:
} }
} }
void initNullSolver() { 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 initFFTSolver() { void initFFTSolver() {
if constexpr (Dim == 2 || Dim == 3) { ippl::ParameterList sp;
ippl::ParameterList sp; sp.add("output_type", FFTSolver_t<double, 3>::GRAD);
sp.add("output_type", FFTSolver_t<T, Dim>::GRAD); sp.add("use_heffte_defaults", false);
sp.add("use_heffte_defaults", false); sp.add("use_pencils", true);
sp.add("use_pencils", true); sp.add("use_reorder", false);
sp.add("use_reorder", false); sp.add("use_gpu_aware", true);
sp.add("use_gpu_aware", true); sp.add("comm", ippl::p2p_pl);
sp.add("comm", ippl::p2p_pl); sp.add("r2c_direction", 0);
sp.add("r2c_direction", 0); initSolverWithParams<FFTSolver_t<double, 3>>(sp);
initSolverWithParams<FFTSolver_t<T, Dim>>(sp);
} else {
throw std::runtime_error("Unsupported dimensionality for FFT solver");
}
} }
void initCGSolver() { }
void initCGSolver() { void initP3MSolver() { }
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() {
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 #endif
src/PartBunch/PartBunch.cpp
+ 105
− 7
View file @ 7e9af2f2
...@@ -2,6 +2,89 @@ ...@@ -2,6 +2,89 @@
#include <boost/numeric/ublas/io.hpp> #include <boost/numeric/ublas/io.hpp>
#include "Utilities/Util.h" #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<> template<>
void PartBunch<double,3>::calcBeamParameters() { void PartBunch<double,3>::calcBeamParameters() {
// Usefull constants // Usefull constants
...@@ -97,11 +180,24 @@ void PartBunch<double,3>::calcBeamParameters() { ...@@ -97,11 +180,24 @@ void PartBunch<double,3>::calcBeamParameters() {
rmax_m(i) = rmax[i]; rmax_m(i) = rmax[i];
rmin_m(i) = rmin[i]; rmin_m(i) = rmin[i];
} }
ippl::Comm->barrier(); 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<> template<>
Inform& PartBunch<double,3>::print(Inform& os) { Inform& PartBunch<double,3>::print(Inform& os) {
...@@ -285,9 +381,7 @@ void PartBunch<double,3>::computeSelfFields() { ...@@ -285,9 +381,7 @@ void PartBunch<double,3>::computeSelfFields() {
//divide charge by a 'grid-cube' volume to get [C/m^3] //divide charge by a 'grid-cube' volume to get [C/m^3]
rho = rho *tmp2; rho = rho *tmp2;
double Npoints = nr_m[0] * nr_m[1] * nr_m[2];
// auto rhoView = rho.getView();
localDensity2 = 0.; localDensity2 = 0.;
ippl::parallel_reduce( ippl::parallel_reduce(
"Density stats", ippl::getRangePolicy(rhoView), "Density stats", ippl::getRangePolicy(rhoView),
...@@ -297,18 +391,22 @@ void PartBunch<double,3>::computeSelfFields() { ...@@ -297,18 +391,22 @@ void PartBunch<double,3>::computeSelfFields() {
}, },
Kokkos::Sum<double>(localDensity2) ); Kokkos::Sum<double>(localDensity2) );
ippl::Comm->reduce(localDensity2, Density2, 1, std::plus<double>()); 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 ); rmsDensity_m = std::sqrt( (1.0 /Npoints) * Density2 / Physics::q_e / Physics::q_e );
this->calcDebyeLength(); this->calcDebyeLength();
m << "gammaz = " << gammaz << endl; m << "gammaz = " << gammaz << endl;
m << "hr_scaled = " << hr_scaled << endl; m << "hr_scaled = " << hr_scaled << endl;
m << "coupling constant= " << getCouplingConstant() << endl;
this->fsolver_m->runSolver(); this->fsolver_m->runSolver();
gather(this->pcontainer_m->E, this->fcontainer_m->getE(), this->pcontainer_m->R); gather(this->pcontainer_m->E, this->fcontainer_m->getE(), this->pcontainer_m->R);
spaceChargeEFieldCheck();
IpplTimings::stopTimer(SolveTimer); IpplTimings::stopTimer(SolveTimer);
} }
......
src/PartBunch/PartBunch.hpp
+ 18
− 39
View file @ 7e9af2f2
...@@ -62,6 +62,14 @@ diagnostics(): calculate statistics and maybe write tp h5 and stat files ...@@ -62,6 +62,14 @@ diagnostics(): calculate statistics and maybe write tp h5 and stat files
#include "Algorithms/PartData.h" #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; using view_type = typename ippl::detail::ViewType<ippl::Vector<double, 3>, 1>::view_type;
...@@ -276,11 +284,6 @@ public: ...@@ -276,11 +284,6 @@ public:
IpplTimings::startTimer(prerun); IpplTimings::startTimer(prerun);
pre_run(); pre_run();
IpplTimings::stopTimer(prerun); IpplTimings::stopTimer(prerun);
/*
end wrmup
*/
} }
void bunchUpdate(); void bunchUpdate();
...@@ -294,39 +297,10 @@ public: ...@@ -294,39 +297,10 @@ public:
return this->pcontainer_m; return this->pcontainer_m;
} }
void setSolver(std::string solver) { 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 pre_run() override ;
public: public:
void updateMoments(){ void updateMoments(){
this->pcontainer_m->updateMoments(); this->pcontainer_m->updateMoments();
...@@ -383,10 +357,11 @@ public: ...@@ -383,10 +357,11 @@ public:
*/ */
double getCouplingConstant() const { double getCouplingConstant() const {
return 1.0; return couplingConstant_m;
} }
void setCouplingConstant(double c) { void setCouplingConstant(double c) {
couplingConstant_m = c;
} }
void calcBeamParameters(); void calcBeamParameters();
...@@ -850,6 +825,10 @@ public: ...@@ -850,6 +825,10 @@ public:
double calculateAngle(double x, double y) { double calculateAngle(double x, double y) {
return 0.0; return 0.0;
} }
// Sanity check functions
void spaceChargeEFieldCheck();
}; };
/* \todo /* \todo
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment