Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
32 changes: 17 additions & 15 deletions alpine/AlpineManager.h
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@

#include <memory>

#include "FEM/FEMInterpolate.hpp"
#include "FieldContainer.hpp"
#include "FieldSolver.hpp"
#include "LoadBalancer.hpp"
Expand All @@ -13,7 +14,6 @@
#include "Random/InverseTransformSampling.h"
#include "Random/NormalDistribution.h"
#include "Random/Randn.h"
#include "FEM/FEMInterpolate.hpp"

using view_type = typename ippl::detail::ViewType<ippl::Vector<double, Dim>, 1>::view_type;

Expand All @@ -35,6 +35,7 @@ class AlpineManager : public ippl::PicManager<T, Dim, ParticleContainer<T, Dim>,
std::string solver_m;
std::string stepMethod_m;
std::vector<std::string> preconditioner_params_m;

public:
AlpineManager(size_type totalP_, int nt_, Vector_t<int, Dim>& nr_, double lbt_,
std::string& solver_, std::string& stepMethod_,
Expand Down Expand Up @@ -98,7 +99,7 @@ class AlpineManager : public ippl::PicManager<T, Dim, ParticleContainer<T, Dim>,

std::vector<std::string> getPreconditionerParams() const { return preconditioner_params_m; };

virtual void dump(){/* default does nothing */};
virtual void dump() { /* default does nothing */ };

void pre_step() override {
Inform m("Pre-step");
Expand All @@ -116,7 +117,7 @@ class AlpineManager : public ippl::PicManager<T, Dim, ParticleContainer<T, Dim>,
m << "Finished time step: " << this->it_m << " time: " << this->time_m << endl;
}

void grid2par() override {
void grid2par() override {
if ((getSolver() == "FEM") || (getSolver() == "FEM_PRECON")) {
gatherFEM();
} else {
Expand All @@ -129,15 +130,15 @@ class AlpineManager : public ippl::PicManager<T, Dim, ParticleContainer<T, Dim>,
}

void gatherFEM() {
using exec_space = typename Kokkos::View<const size_t*>::execution_space;
using policy_type = Kokkos::RangePolicy<exec_space>;
using exec_space = typename Kokkos::View<const size_t*>::execution_space;
using policy_type = Kokkos::RangePolicy<exec_space>;
size_type localParticles = this->pcontainer_m->getLocalNum();
policy_type iteration_policy(0, localParticles);

// get the Finite Element space from the solver,
// get the Finite Element space from the solver,
// since the interpolation depends on the space
auto* solver = dynamic_cast<FieldSolver_t*>(this->fsolver_m.get());
auto& space = solver->getSpace();
auto& space = solver->getSpace();

interpolate_grad_to_diracs(this->pcontainer_m->E, this->fcontainer_m->getPhi(),
this->pcontainer_m->R, space, iteration_policy);
Expand Down Expand Up @@ -182,14 +183,14 @@ class AlpineManager : public ippl::PicManager<T, Dim, ParticleContainer<T, Dim>,
double Q = Q_m;
size_type localParticles = this->pcontainer_m->getLocalNum();

using exec_space = typename Kokkos::View<const size_t*>::execution_space;
using exec_space = typename Kokkos::View<const size_t*>::execution_space;
using policy_type = Kokkos::RangePolicy<exec_space>;
policy_type iteration_policy(0, localParticles);

// get the Finite Element space from the solver,
// get the Finite Element space from the solver,
// since the interpolation depends on the space
auto* solver = dynamic_cast<FieldSolver_t*>(this->fsolver_m.get());
auto& space = solver->getSpace();
auto& space = solver->getSpace();

assemble_rhs_from_particles(*q, *rho, *R, space, iteration_policy);

Expand Down Expand Up @@ -220,12 +221,13 @@ class AlpineManager : public ippl::PicManager<T, Dim, ParticleContainer<T, Dim>,
}

void getDensity(Field_t<Dim>* rho) {
Vector_t<double, Dim> rmin = rmin_m;
Vector_t<double, Dim> rmax = rmax_m;
Vector_t<double, Dim> hr = this->hr_m;
double Q = Q_m;
Vector_t<double, Dim> rmin = rmin_m;
Vector_t<double, Dim> rmax = rmax_m;
Vector_t<double, Dim> hr = this->hr_m;
double Q = Q_m;

if ((this->fsolver_m->getStype() != "FEM") && (this->fsolver_m->getStype() != "FEM_PRECON")) {
if ((this->fsolver_m->getStype() != "FEM")
&& (this->fsolver_m->getStype() != "FEM_PRECON")) {
double cellVolume = std::reduce(hr.begin(), hr.end(), 1., std::multiplies<double>());
(*rho) = (*rho) / cellVolume;
}
Expand Down
13 changes: 6 additions & 7 deletions alpine/BumponTailInstabilityManager.h
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
#ifndef IPPL_BUMPON_TAIL_INSTABILITY_MANAGER_H
#define IPPL_BUMPON_TAIL_INSTABILITY_MANAGER_H

#include <memory>
#include <filesystem>
#include <memory>

#include "AlpineManager.h"
#include "FieldContainer.hpp"
Expand Down Expand Up @@ -83,9 +83,9 @@ class BumponTailInstabilityManager : public AlpineManager<T, Dim> {
public:
void pre_run() override {
Inform m("Pre Run");
const double pi = Kokkos::numbers::pi_v<T>;

const double pi = Kokkos::numbers::pi_v<T>;

if (this->solver_m == "OPEN") {
throw IpplException("BumpOnTailInstability",
"Open boundaries solver incompatible with this simulation!");
Expand Down Expand Up @@ -393,9 +393,8 @@ class BumponTailInstabilityManager : public AlpineManager<T, Dim> {
auto& Pi = pc->P;
for (unsigned d = allDims ? 0 : Dim - 1; d < Dim; d++) {
Kokkos::parallel_for(
"Copy phase space", pcount, KOKKOS_CLASS_LAMBDA(const size_t i) {
phase(i) = {Ri(i)[d], Pi(i)[d]};
});
"Copy phase space", pcount,
KOKKOS_CLASS_LAMBDA(const size_t i) { phase(i) = {Ri(i)[d], Pi(i)[d]}; });
phaseSpace = 0;
Kokkos::fence();
scatter(pc->q, phaseSpace, phase);
Expand Down
3 changes: 2 additions & 1 deletion alpine/FieldContainer.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,8 @@ class FieldContainer {
void initializeFields(std::string stype_m = "") {
E_m.initialize(mesh_m, fl_m);
rho_m.initialize(mesh_m, fl_m);
if ((stype_m == "CG") || (stype_m == "PCG") || (stype_m == "FEM") || (stype_m == "FEM_PRECON")) {
if ((stype_m == "CG") || (stype_m == "PCG") || (stype_m == "FEM")
|| (stype_m == "FEM_PRECON")) {
phi_m.initialize(mesh_m, fl_m);
}
}
Expand Down
31 changes: 16 additions & 15 deletions alpine/FieldSolver.hpp
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
#ifndef IPPL_FIELD_SOLVER_H
#define IPPL_FIELD_SOLVER_H

#include <memory>
#include <filesystem>
#include <memory>

#include "Manager/BaseManager.h"
#include "Manager/FieldSolverBase.h"
Expand Down Expand Up @@ -63,11 +63,11 @@ class FieldSolver : public ippl::FieldSolverBase<T, Dim> {
// CG requires explicit periodic boundary conditions while the periodic Poisson solver
// simply assumes them
typedef ippl::BConds<Field<T, Dim>, Dim> bc_type;
if ((this->getStype() == "CG") || (this->getStype() == "PCG") || (this->getStype() == "FEM") ||
(this->getStype() == "FEM_PRECON")) {
if ((this->getStype() == "CG") || (this->getStype() == "PCG") || (this->getStype() == "FEM")
|| (this->getStype() == "FEM_PRECON")) {
bc_type allPeriodic;
for (unsigned int i = 0; i < 2 * Dim; ++i) {
allPeriodic[i] = std::make_shared<ippl::PeriodicFace<Field<T, Dim>>>(i);
allPeriodic[i] = std::make_shared<ippl::PeriodicFace<Field<T, Dim>>>(i);
}
phi_m->setFieldBC(allPeriodic);
if ((this->getStype() == "FEM") || (this->getStype() == "FEM_PRECON")) {
Expand All @@ -77,10 +77,10 @@ class FieldSolver : public ippl::FieldSolverBase<T, Dim> {
}

void runSolver() override {
if ((this->getStype() == "CG") || (this->getStype() == "PCG") || (this->getStype() == "FEM") ||
(this->getStype() == "FEM_PRECON")) {
if ((this->getStype() == "CG") || (this->getStype() == "PCG") || (this->getStype() == "FEM")
|| (this->getStype() == "FEM_PRECON")) {
int iterations = 0;
int residue = 0;
int residue = 0;

if (this->getStype() == "FEM") {
FEMSolver_t<T, Dim>& solver = std::get<FEMSolver_t<T, Dim>>(this->getSolver());
Expand All @@ -89,7 +89,8 @@ class FieldSolver : public ippl::FieldSolverBase<T, Dim> {
iterations = solver.getIterationCount();
residue = solver.getResidue();
} else if (this->getStype() == "FEM_PRECON") {
FEMPreconSolver_t<T, Dim>& solver = std::get<FEMPreconSolver_t<T, Dim>>(this->getSolver());
FEMPreconSolver_t<T, Dim>& solver =
std::get<FEMPreconSolver_t<T, Dim>>(this->getSolver());
solver.solve();

iterations = solver.getIterationCount();
Expand All @@ -105,8 +106,8 @@ class FieldSolver : public ippl::FieldSolverBase<T, Dim> {
if (ippl::Comm->rank() == 0) {
std::filesystem::create_directory("data_CG");
std::stringstream fname;
if ((this->getStype() == "CG") || (this->getStype() == "FEM") ||
(this->getStype() == "FEM_PRECON")) {
if ((this->getStype() == "CG") || (this->getStype() == "FEM")
|| (this->getStype() == "FEM_PRECON")) {
fname << "data_CG/CG_";
} else {
fname << "data_CG/";
Expand Down Expand Up @@ -153,10 +154,10 @@ class FieldSolver : public ippl::FieldSolverBase<T, Dim> {

solver.setRhs(*rho_m);

if constexpr ((std::is_same_v<Solver, CGSolver_t<T, Dim>>) ||
(std::is_same_v<Solver, FEMSolver_t<T, Dim>>) ||
(std::is_same_v<Solver, FEMPreconSolver_t<T, Dim>>)) {
// The CG solver and FEMPoissonSolver compute the potential
if constexpr ((std::is_same_v<Solver, CGSolver_t<T, Dim>>)
|| (std::is_same_v<Solver, FEMSolver_t<T, Dim>>)
|| (std::is_same_v<Solver, FEMPreconSolver_t<T, Dim>>)) {
// The CG solver and FEMPoissonSolver compute the potential
// directly and use this to get the electric field
solver.setLhs(*phi_m);
solver.setGradient(*E_m);
Expand Down Expand Up @@ -292,7 +293,7 @@ class FieldSolver : public ippl::FieldSolverBase<T, Dim> {
sp.add("richardson_iterations", richardson_iterations);
sp.add("communication", communication);
sp.add("ssor_omega", ssor_omega);

initSolverWithParams<FEMPreconSolver_t<T, Dim>>(sp);
}

Expand Down
4 changes: 2 additions & 2 deletions alpine/LandauDamping.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@ int main(int argc, char* argv[]) {
Inform msg(TestName);
Inform msg2all(TestName, INFORM_ALL_NODES);

static IpplTimings::TimerRef mainTimer = IpplTimings::getTimer("total");
static IpplTimings::TimerRef mainTimer = IpplTimings::getTimer("total");
static IpplTimings::TimerRef initializeTimer = IpplTimings::getTimer("initialize");
IpplTimings::startTimer(mainTimer);
IpplTimings::startTimer(initializeTimer);
Expand Down Expand Up @@ -82,7 +82,7 @@ int main(int argc, char* argv[]) {
manager.pre_run();

IpplTimings::stopTimer(initializeTimer);

manager.setTime(0.0);

msg << "Starting iterations ..." << endl;
Expand Down
29 changes: 14 additions & 15 deletions alpine/LandauDampingManager.h
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
#ifndef IPPL_LANDAU_DAMPING_MANAGER_H
#define IPPL_LANDAU_DAMPING_MANAGER_H

#include <memory>
#include <filesystem>
#include <memory>

#include "AlpineManager.h"
#include "FieldContainer.hpp"
Expand Down Expand Up @@ -67,7 +67,7 @@ class LandauDampingManager : public AlpineManager<T, Dim> {
Inform m("Pre Run");

const double pi = Kokkos::numbers::pi_v<T>;

if (this->solver_m == "OPEN") {
throw IpplException("LandauDamping",
"Open boundaries solver incompatible with this simulation!");
Expand All @@ -84,7 +84,7 @@ class LandauDampingManager : public AlpineManager<T, Dim> {
this->rmax_m = 2 * pi / this->kw_m;

bool isFEM = ((this->getSolver() == "FEM") || (this->getSolver() == "FEM_PRECON"));

Vector<int, Dim> nElements = this->nr_m - 1;
if (isFEM) {
this->hr_m = this->rmax_m / nElements;
Expand All @@ -110,8 +110,7 @@ class LandauDampingManager : public AlpineManager<T, Dim> {
this->isAllPeriodic_m));

this->setParticleContainer(std::make_shared<ParticleContainer_t>(
this->fcontainer_m->getMesh(), this->fcontainer_m->getFL(),
isFEM));
this->fcontainer_m->getMesh(), this->fcontainer_m->getFL(), isFEM));

this->fcontainer_m->initializeFields(this->solver_m);

Expand Down Expand Up @@ -191,7 +190,8 @@ class LandauDampingManager : public AlpineManager<T, Dim> {
this->fcontainer_m->getRho().getFieldRangePolicy(),
KOKKOS_LAMBDA(const index_array_type& args) {
// local to global index conversion
Vector_t<double, Dim> xvec = (args + lDom.first() - nghost + 0.5*(!isFEM)) * hr + origin;
Vector_t<double, Dim> xvec =
(args + lDom.first() - nghost + 0.5 * (!isFEM)) * hr + origin;

// ippl::apply accesses the view at the given indices and obtains a
// reference; see src/Expression/IpplOperations.h
Expand Down Expand Up @@ -246,7 +246,7 @@ class LandauDampingManager : public AlpineManager<T, Dim> {
this->pcontainer_m->q = this->Q_m / totalP;

// For FEM need an update due to node-centering, as periodic BCs mean
// that a particle at R=0 is equivalent to R=1 so it could be on the
// that a particle at R=0 is equivalent to R=1 so it could be on the
// wrong rank and needs to be sent over.
if (isFEM) {
this->pcontainer_m->update();
Expand Down Expand Up @@ -325,7 +325,7 @@ class LandauDampingManager : public AlpineManager<T, Dim> {

if ((this->getSolver() == "FEM") || (this->getSolver() == "FEM_PRECON")) {
// When using FEM, we only have E on particles
// so we use the dump function which computes the
// so we use the dump function which computes the
// energy using the particles instead of the field.
dumpLandau();
} else {
Expand Down Expand Up @@ -385,15 +385,15 @@ class LandauDampingManager : public AlpineManager<T, Dim> {
ippl::Comm->barrier();
}

// Overloaded dumpLandau which computes the E-field energy using the particles
// instead of using the E-field on the grid (as above). Since we have E for
// Overloaded dumpLandau which computes the E-field energy using the particles
// instead of using the E-field on the grid (as above). Since we have E for
// each particle, we treat the particles as Monte-Carlo samples to compute
// the energy integral.
void dumpLandau() {
auto Eview = this->pcontainer_m->E.getView();
auto Eview = this->pcontainer_m->E.getView();
size_type localParticles = this->pcontainer_m->getLocalNum();

using exec_space = typename Kokkos::View<const size_t*>::execution_space;
using exec_space = typename Kokkos::View<const size_t*>::execution_space;
using policy_type = Kokkos::RangePolicy<exec_space>;
policy_type iteration_policy(0, localParticles);

Expand All @@ -412,9 +412,8 @@ class LandauDampingManager : public AlpineManager<T, Dim> {

// MC integration: divide by no. of particles N and multiply by volume
ippl::Vector<T, Dim> domain_size = this->rmax_m - this->rmin_m;
double fieldEnergy =
std::reduce(domain_size.begin(), domain_size.end(),
globaltemp, std::multiplies<double>());
double fieldEnergy = std::reduce(domain_size.begin(), domain_size.end(), globaltemp,
std::multiplies<double>());

fieldEnergy = fieldEnergy / this->totalP_m;

Expand Down
Loading
Loading