Skip to content

Commit

Permalink
Merge branch 'fix-copies' into 'master'
Browse files Browse the repository at this point in the history 
Bug fixes

See merge request OPAL/Libraries/ippl!171
  • Loading branch information
@Arc676
Arc676 committed May 7, 2023
2 parents 6d08224 + 52d0e0f commit 756b0cb
Show file tree
Hide file tree
Showing 14 changed files with 178 additions and 181 deletions.
8 changes: 4 additions & 4 deletions CMakeLists.txt
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
cmake_minimum_required (VERSION 3.1)
project (IPPL CXX)
set (IPPL_VERSION_MAJOR 1)
set (IPPL_VERSION_MINOR 1.4)
set (IPPL_VERSION_MAJOR 3)
set (IPPL_VERSION_MINOR 0.1)

set (IPPL_VERSION_NAME "V${IPPL_VERSION_MAJOR}.${IPPL_VERSION_MINOR}")

Expand Down Expand Up @@ -37,8 +37,8 @@ set (CMAKE_CXX_FLAGS_RELWITHDEBINFO "-O3 -g ")
set (CMAKE_CXX_FLAGS_RELEASE "-O3")
set (CMAKE_CXX_FLAGS_DEBUG "-O0 -g")


if (${CMAKE_CXX_COMPILER} MATCHES "GNU")
# Compile using C++20 for CPU builds
if (${CMAKE_CXX_COMPILER} MATCHES ".*mpicxx")
set (CMAKE_CXX_STANDARD 20)
endif ()

Expand Down
18 changes: 9 additions & 9 deletions alpine/BumponTailInstability.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,15 +7,15 @@
// ny... = No. cell-centered points in the y-, z-, ...-direction
// Np = Total no. of macro-particles in the simulation
// Nt = Number of time steps
// stype = Field solver type e.g., FFT
// stype = Field solver type (FFT and CG supported)
// 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
// 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
// srun ./BumponTailInstability 128 128 128 10000 10 FFT 0.01 --overallocate 2.0 --info 10
// Change the TestName to TwoStreamInstability or BumponTailInstability
// in order to simulate the Two stream instability or bump on tail instability
// cases
Expand Down Expand Up @@ -314,24 +314,24 @@ int main(int argc, char* argv[]) {
hr[d] = rmax[d] / nr[d];
}
Vector_t<Dim> origin = rmin;
const double dt = 0.5 * hr[0]; // 0.05
const double dt = std::min(.05, 0.5 * *std::min_element(hr.begin(), hr.end()));

const bool isAllPeriodic = true;
Mesh_t<Dim> mesh(domain, hr, origin);
FieldLayout_t<Dim> FL(domain, decomp, isAllPeriodic);
PLayout_t<Dim> PL(FL, mesh);

// Q = -\int\int f dx dv
double Q = std::reduce(rmax.begin(), rmax.end(), -1., std::multiplies<double>());
P = std::make_shared<bunch_type>(PL, hr, rmin, rmax, decomp, Q);
double Q = std::reduce(rmax.begin(), rmax.end(), -1., std::multiplies<double>());
std::string solver = argv[arg++];
P = std::make_shared<bunch_type>(PL, hr, rmin, rmax, decomp, Q, solver);

P->nr_m = nr;

P->initializeFields(mesh, FL);

bunch_type bunchBuffer(PL);

P->stype_m = argv[arg++];
P->initSolver();
P->time_m = 0.0;
P->loadbalancethreshold_m = std::atof(argv[arg++]);
Expand Down Expand Up @@ -433,13 +433,13 @@ int main(int argc, char* argv[]) {

IpplTimings::startTimer(DummySolveTimer);
P->rho_m = 0.0;
P->solver_mp->solve();
P->runSolver();
IpplTimings::stopTimer(DummySolveTimer);

P->scatterCIC(totalP, 0, hr);

IpplTimings::startTimer(SolveTimer);
P->solver_mp->solve();
P->runSolver();
IpplTimings::stopTimer(SolveTimer);

P->gatherCIC();
Expand Down Expand Up @@ -489,7 +489,7 @@ int main(int argc, char* argv[]) {

// Field solve
IpplTimings::startTimer(SolveTimer);
P->solver_mp->solve();
P->runSolver();
IpplTimings::stopTimer(SolveTimer);

// gather E field
Expand Down
118 changes: 86 additions & 32 deletions alpine/ChargedParticles.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,9 +59,16 @@ using VField_t = Field<Vector_t<Dim>, Dim>;

// heFFTe does not support 1D FFTs, so we switch to CG in the 1D case
template <unsigned Dim = 3>
using Solver_t = std::conditional_t<
Dim == 1, ippl::ElectrostaticsCG<double, double, Dim, Mesh_t<Dim>, Centering_t<Dim>>,
ippl::FFTPeriodicPoissonSolver<Vector_t<Dim>, double, Dim, Mesh_t<Dim>, Centering_t<Dim>>>;
using CGSolver_t = ippl::ElectrostaticsCG<double, double, Dim, Mesh_t<Dim>, Centering_t<Dim>>;

template <unsigned Dim = 3>
using FFTSolver_t =
ippl::FFTPeriodicPoissonSolver<Vector_t<Dim>, double, Dim, Mesh_t<Dim>, Centering_t<Dim>>;

template <unsigned Dim = 3>
using Solver_t =
std::conditional_t<Dim == 2 || Dim == 3, std::variant<CGSolver_t<Dim>, FFTSolver_t<Dim>>,
std::variant<CGSolver_t<Dim>>>;

const double pi = std::acos(-1.0);

Expand Down Expand Up @@ -161,12 +168,14 @@ class ChargedParticles : public ippl::ParticleBase<PLayout> {
Vector_t<Dim> rmin_m;
Vector_t<Dim> rmax_m;

double Q_m;

std::string stype_m;

std::shared_ptr<Solver_t<Dim>> solver_mp;
double Q_m;

private:
Solver_t<Dim> solver_m;

public:
double time_m;

double rhoNorm_m;
Expand All @@ -192,11 +201,12 @@ class ChargedParticles : public ippl::ParticleBase<PLayout> {
}

ChargedParticles(PLayout& pl, Vector_t<Dim> hr, Vector_t<Dim> rmin, Vector_t<Dim> rmax,
ippl::e_dim_tag decomp[Dim], double Q)
ippl::e_dim_tag decomp[Dim], double Q, std::string solver)
: ippl::ParticleBase<PLayout>(pl)
, hr_m(hr)
, rmin_m(rmin)
, rmax_m(rmax)
, stype_m(solver)
, Q_m(Q) {
registerAttributes();
for (unsigned int i = 0; i < Dim; i++) {
Expand All @@ -209,7 +219,7 @@ class ChargedParticles : public ippl::ParticleBase<PLayout> {
void setPotentialBCs() {
// CG requires explicit periodic boundary conditions while the periodic Poisson solver
// simply assumes them
if constexpr (Dim == 1) {
if (stype_m == "CG") {
for (unsigned int i = 0; i < 2 * Dim; ++i) {
allPeriodic[i] = std::make_shared<
ippl::PeriodicFace<double, Dim, Mesh_t<Dim>, Centering_t<Dim>>>(i);
Expand All @@ -235,7 +245,7 @@ class ChargedParticles : public ippl::ParticleBase<PLayout> {
IpplTimings::startTimer(tupdateLayout);
this->E_m.updateLayout(fl);
this->rho_m.updateLayout(fl);
if constexpr (Dim == 1) {
if (stype_m == "CG") {
this->phi_m.updateLayout(fl);
phi_m.setFieldBC(allPeriodic);
}
Expand All @@ -255,7 +265,7 @@ class ChargedParticles : public ippl::ParticleBase<PLayout> {
void initializeFields(Mesh_t<Dim>& mesh, FieldLayout_t<Dim>& fl) {
E_m.initialize(mesh, fl);
rho_m.initialize(mesh, fl);
if constexpr (Dim == 1) {
if (stype_m == "CG") {
phi_m.initialize(mesh, fl);
phi_m.setFieldBC(allPeriodic);
}
Expand All @@ -276,7 +286,11 @@ class ChargedParticles : public ippl::ParticleBase<PLayout> {
}
// Update
this->updateLayout(fl, mesh, buffer, isFirstRepartition);
this->solver_mp->setRhs(rho_m);
if constexpr (Dim == 2 || Dim == 3) {
if (stype_m == "FFT") {
std::get<FFTSolver_t<Dim>>(solver_m).setRhs(rho_m);
}
}
}

bool balance(size_type totalP, const unsigned int nstep) {
Expand Down Expand Up @@ -388,43 +402,83 @@ class ChargedParticles : public ippl::ParticleBase<PLayout> {
}
}

void initSolver(const ippl::ParameterList& sp) {
solver_mp = std::make_shared<Solver_t<Dim>>();
void runSolver() {
if (stype_m == "CG") {
CGSolver_t<Dim>& solver = std::get<CGSolver_t<Dim>>(solver_m);
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 (time_m == 0 && iterations == 0) {
log << "time,residue,iterations" << endl;
}
// Don't print the dummy solve
if (time_m > 0 || iterations > 0) {
log << time_m << "," << solver.getResidue() << "," << iterations << endl;
}
}
Ippl::Comm->barrier();
} else if (stype_m == "FFT") {
if constexpr (Dim == 2 || Dim == 3) {
std::get<FFTSolver_t<Dim>>(solver_m).solve();
}
} else {
throw std::runtime_error("Unknown solver type");
}
}

template <typename Solver>
void initSolverWithParams(const ippl::ParameterList& sp) {
solver_m = Solver();
Solver& solver = std::get<Solver>(solver_m);

solver_mp->mergeParameters(sp);
solver.mergeParameters(sp);

solver_mp->setRhs(rho_m);
solver.setRhs(rho_m);

if constexpr (Dim == 1) {
if constexpr (std::is_same_v<Solver, CGSolver_t<Dim>>) {
// The CG solver computes the potential directly and
// uses this to get the electric field
solver_mp->setLhs(phi_m);
solver_mp->setGradient(E_m);
} else {
solver.setLhs(phi_m);
solver.setGradient(E_m);
} else if constexpr (std::is_same_v<Solver, FFTSolver_t<Dim>>) {
// The periodic Poisson solver computes the electric
// field directly
solver_mp->setLhs(E_m);
solver.setLhs(E_m);
}
}

void initCGSolver() {
ippl::ParameterList sp;
sp.add("output_type", Solver_t<Dim>::GRAD);
sp.add("output_type", CGSolver_t<Dim>::GRAD);
// Increase tolerance in the 1D case
sp.add("tolerance", 1e-10);

initSolver(sp);
initSolverWithParams<CGSolver_t<Dim>>(sp);
}

void initFFTSolver() {
ippl::ParameterList sp;
sp.add("output_type", Solver_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);

initSolver(sp);
if constexpr (Dim == 2 || Dim == 3) {
ippl::ParameterList sp;
sp.add("output_type", FFTSolver_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<Dim>>(sp);
} else {
throw std::runtime_error("Unsupported dimensionality for FFT solver");
}
}

void dumpData() {
Expand Down
19 changes: 9 additions & 10 deletions alpine/LandauDamping.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,15 +7,15 @@
// ny... = No. cell-centered points in the y-, z-, ...-direction
// Np = Total no. of macro-particles in the simulation
// Nt = Number of time steps
// stype = Field solver type e.g., FFT
// stype = Field solver type (FFT and CG supported)
// 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
// 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 ./LandauDamping 128 128 128 10000 10 FFT 0.01 2.0 --info 10
// srun ./LandauDamping 128 128 128 10000 10 FFT 0.01 --overallocate 2.0 --info 10
//
// Copyright (c) 2021, Sriramkrishnan Muralikrishnan,
// Paul Scherrer Institut, Villigen PSI, Switzerland
Expand Down Expand Up @@ -203,23 +203,22 @@ int main(int argc, char* argv[]) {
// Q = -\int\int f dx dv
double Q = std::reduce(rmax.begin(), rmax.end(), -1., std::multiplies<double>());
Vector_t<Dim> origin = rmin;
const double dt = 0.5 * hr[0];
const double dt = std::min(.05, 0.5 * *std::min_element(hr.begin(), hr.end()));

const bool isAllPeriodic = true;
Mesh_t<Dim> mesh(domain, hr, origin);
FieldLayout_t<Dim> FL(domain, decomp, isAllPeriodic);
PLayout_t<Dim> PL(FL, mesh);

P = std::make_unique<bunch_type>(PL, hr, rmin, rmax, decomp, Q);
std::string solver = argv[arg++];
P = std::make_unique<bunch_type>(PL, hr, rmin, rmax, decomp, Q, solver);

P->nr_m = nr;

P->E_m.initialize(mesh, FL);
P->rho_m.initialize(mesh, FL);
P->initializeFields(mesh, FL);

bunch_type bunchBuffer(PL);

P->stype_m = argv[arg++];
P->initSolver();
P->time_m = 0.0;
P->loadbalancethreshold_m = std::atof(argv[arg++]);
Expand Down Expand Up @@ -303,13 +302,13 @@ int main(int argc, char* argv[]) {

IpplTimings::startTimer(DummySolveTimer);
P->rho_m = 0.0;
P->solver_mp->solve();
P->runSolver();
IpplTimings::stopTimer(DummySolveTimer);

P->scatterCIC(totalP, 0, hr);

IpplTimings::startTimer(SolveTimer);
P->solver_mp->solve();
P->runSolver();
IpplTimings::stopTimer(SolveTimer);

P->gatherCIC();
Expand Down Expand Up @@ -360,7 +359,7 @@ int main(int argc, char* argv[]) {

// Field solve
IpplTimings::startTimer(SolveTimer);
P->solver_mp->solve();
P->runSolver();
IpplTimings::stopTimer(SolveTimer);

// gather E field
Expand Down
Loading

0 comments on commit 756b0cb

Please sign in to comment.