From baa80cf2ee4eef11b5b5f1ab0102ae9f0ac8a516 Mon Sep 17 00:00:00 2001 From: Matthias Frey Date: Wed, 15 Mar 2023 13:58:05 +0000 Subject: [PATCH 01/12] remove default template arguments; change M to Mesh and C to Cell to be more clea --- src/FFT/FFT.h | 2 +- src/Field/BConds.h | 3 +- src/Field/BcTypes.h | 15 ++++------ src/Field/Field.h | 7 ++--- src/Field/Field.hpp | 28 +++++++++--------- src/Field/FieldOperations.hpp | 38 ++++++++++++------------- src/Particle/IntNGP.h | 32 ++++++++++----------- src/Solver/Electrostatics.h | 15 +++++----- src/Solver/ElectrostaticsCG.h | 13 ++++----- src/Solver/FFTPeriodicPoissonSolver.h | 11 ++++--- src/Solver/FFTPeriodicPoissonSolver.hpp | 12 ++++---- src/Solver/FFTPoissonSolver.h | 23 +++++++-------- src/Solver/FFTPoissonSolver.hpp | 38 ++++++++++++------------- src/Solver/PCG.h | 4 +-- src/Solver/Solver.h | 7 ++--- src/Solver/SolverAlgorithm.h | 7 ++--- 16 files changed, 121 insertions(+), 134 deletions(-) diff --git a/src/FFT/FFT.h b/src/FFT/FFT.h index 431cae6d1..f7b626ac7 100644 --- a/src/FFT/FFT.h +++ b/src/FFT/FFT.h @@ -118,7 +118,7 @@ namespace ippl { /** Non-specialized FFT class. We specialize based on Transform tag class */ - template > + template class FFT {}; /** diff --git a/src/Field/BConds.h b/src/Field/BConds.h index b33af4d42..623fc1d79 100644 --- a/src/Field/BConds.h +++ b/src/Field/BConds.h @@ -37,8 +37,7 @@ namespace ippl { template std::ostream& operator<<(std::ostream&, const BConds&); - template , - class Cell = typename Mesh::DefaultCentering> + template class BConds { public: using bc_type = detail::BCondBase; diff --git a/src/Field/BcTypes.h b/src/Field/BcTypes.h index fc55c65e1..16d083360 100644 --- a/src/Field/BcTypes.h +++ b/src/Field/BcTypes.h @@ -98,8 +98,7 @@ namespace ippl { } // namespace detail - template , - class Cell = typename Mesh::DefaultCentering> + template class ExtrapolateFace : public detail::BCondBase { public: // Constructor takes zero, one, or two int's specifying components of @@ -132,8 +131,7 @@ namespace ippl { T slope_m; }; - template , - class Cell = typename Mesh::DefaultCentering> + template class NoBcFace : public detail::BCondBase { public: using Field_t = typename detail::BCondBase::Field_t; @@ -146,8 +144,7 @@ namespace ippl { virtual void write(std::ostream& out) const; }; - template , - class Cell = typename Mesh::DefaultCentering> + template class ConstantFace : public ExtrapolateFace { public: ConstantFace(unsigned int face, T constant) @@ -158,8 +155,7 @@ namespace ippl { virtual void write(std::ostream& out) const; }; - template , - class Cell = typename Mesh::DefaultCentering> + template class ZeroFace : public ConstantFace { public: ZeroFace(unsigned face) @@ -170,8 +166,7 @@ namespace ippl { virtual void write(std::ostream& out) const; }; - template , - class Cell = typename Mesh::DefaultCentering> + template class PeriodicFace : public detail::BCondBase { public: using face_neighbor_type = std::array, 2 * Dim>; diff --git a/src/Field/Field.h b/src/Field/Field.h index 584de9da7..2aabb0f41 100644 --- a/src/Field/Field.h +++ b/src/Field/Field.h @@ -24,15 +24,14 @@ namespace ippl { - template , - class C = typename M::DefaultCentering> + template class Field : public BareField { public: typedef T type; static constexpr unsigned dimension = Dim; - using Mesh_t = M; - using Centering_t = C; + using Mesh_t = Mesh; + using Centering_t = Cell; using Layout_t = FieldLayout; using BareField_t = BareField; using view_type = typename BareField_t::view_type; diff --git a/src/Field/Field.hpp b/src/Field/Field.hpp index c95acf998..8ed19b0a6 100644 --- a/src/Field/Field.hpp +++ b/src/Field/Field.hpp @@ -19,8 +19,8 @@ namespace ippl { namespace detail { - template - struct isExpression> : std::true_type {}; + template + struct isExpression> : std::true_type {}; } // namespace detail ////////////////////////////////////////////////////////////////////////// @@ -28,16 +28,16 @@ namespace ippl { // 'initialize' function before doing anything else. There are no special // checks in the rest of the Field methods to check that the Field has // been properly initialized - template - Field::Field() + template + Field::Field() : BareField() , mesh_m(nullptr) , bc_m() {} ////////////////////////////////////////////////////////////////////////// // Constructors which include a Mesh object as argument - template - Field::Field(Mesh_t& m, Layout_t& l, int nghost) + template + Field::Field(Mesh_t& m, Layout_t& l, int nghost) : BareField(l, nghost) , mesh_m(&m) { for (unsigned int face = 0; face < 2 * Dim; ++face) { @@ -47,8 +47,8 @@ namespace ippl { ////////////////////////////////////////////////////////////////////////// // Initialize the Field, also specifying a mesh - template - void Field::initialize(Mesh_t& m, Layout_t& l, int nghost) { + template + void Field::initialize(Mesh_t& m, Layout_t& l, int nghost) { BareField::initialize(l, nghost); mesh_m = &m; for (unsigned int face = 0; face < 2 * Dim; ++face) { @@ -56,19 +56,19 @@ namespace ippl { } } - template - T Field::getVolumeIntegral() const { + template + T Field::getVolumeIntegral() const { typename M::value_type dV = mesh_m->getCellVolume(); return this->sum() * dV; } - template - T Field::getVolumeAverage() const { + template + T Field::getVolumeAverage() const { return getVolumeIntegral() / mesh_m->getMeshVolume(); } - template - void Field::updateLayout(Layout_t& l, int nghost) { + template + void Field::updateLayout(Layout_t& l, int nghost) { BareField::updateLayout(l, nghost); } diff --git a/src/Field/FieldOperations.hpp b/src/Field/FieldOperations.hpp index d001f4b47..bd7dc4cc4 100644 --- a/src/Field/FieldOperations.hpp +++ b/src/Field/FieldOperations.hpp @@ -23,8 +23,8 @@ namespace ippl { * @param f2 second field * @return Result of f1^T f2 */ - template > - T innerProduct(const Field& f1, const Field& f2) { + template + T innerProduct(const Field& f1, const Field& f2) { T sum = 0; auto view1 = f1.getView(); auto view2 = f2.getView(); @@ -46,8 +46,8 @@ namespace ippl { * @param p desired norm (default 2) * @return The desired norm of the field */ - template - T norm(const Field& field, int p = 2) { + template + T norm(const Field& field, int p = 2) { T local = 0; auto view = field.getView(); switch (p) { @@ -86,8 +86,8 @@ namespace ippl { * User interface of gradient in three dimensions. * @param u field */ - template - detail::meta_grad> grad(Field& u) { + template + detail::meta_grad> grad(Field& u) { u.fillHalo(); BConds& bcField = u.getFieldBC(); bcField.apply(u); @@ -98,15 +98,15 @@ namespace ippl { yvector[1] = 0.5 / mesh.getMeshSpacing(1); typename M::vector_type zvector(0); zvector[2] = 0.5 / mesh.getMeshSpacing(2); - return detail::meta_grad>(u, xvector, yvector, zvector); + return detail::meta_grad>(u, xvector, yvector, zvector); } /*! * User interface of divergence in three dimensions. * @param u field */ - template - detail::meta_div> div(Field& u) { + template + detail::meta_div> div(Field& u) { u.fillHalo(); BConds& bcField = u.getFieldBC(); bcField.apply(u); @@ -117,15 +117,15 @@ namespace ippl { yvector[1] = 0.5 / mesh.getMeshSpacing(1); typename M::vector_type zvector(0); zvector[2] = 0.5 / mesh.getMeshSpacing(2); - return detail::meta_div>(u, xvector, yvector, zvector); + return detail::meta_div>(u, xvector, yvector, zvector); } /*! * User interface of Laplacian in three dimensions. * @param u field */ - template - detail::meta_laplace> laplace(Field& u) { + template + detail::meta_laplace> laplace(Field& u) { u.fillHalo(); BConds& bcField = u.getFieldBC(); bcField.apply(u); @@ -134,15 +134,15 @@ namespace ippl { hvector[0] = 1.0 / std::pow(mesh.getMeshSpacing(0), 2); hvector[1] = 1.0 / std::pow(mesh.getMeshSpacing(1), 2); hvector[2] = 1.0 / std::pow(mesh.getMeshSpacing(2), 2); - return detail::meta_laplace>(u, hvector); + return detail::meta_laplace>(u, hvector); } /*! * User interface of curl in three dimensions. * @param u field */ - template - detail::meta_curl> curl(Field& u) { + template + detail::meta_curl> curl(Field& u) { u.fillHalo(); BConds& bcField = u.getFieldBC(); bcField.apply(u); @@ -155,15 +155,15 @@ namespace ippl { zvector[2] = 1.0; typename M::vector_type hvector(0); hvector = mesh.getMeshSpacing(); - return detail::meta_curl>(u, xvector, yvector, zvector, hvector); + return detail::meta_curl>(u, xvector, yvector, zvector, hvector); } /*! * User interface of Hessian in three dimensions. * @param u field */ - template - detail::meta_hess> hess(Field& u) { + template + detail::meta_hess> hess(Field& u) { u.fillHalo(); BConds& bcField = u.getFieldBC(); bcField.apply(u); @@ -176,6 +176,6 @@ namespace ippl { zvector[2] = 1.0; typename M::vector_type hvector(0); hvector = mesh.getMeshSpacing(); - return detail::meta_hess>(u, xvector, yvector, zvector, hvector); + return detail::meta_hess>(u, xvector, yvector, zvector, hvector); } } // namespace ippl diff --git a/src/Particle/IntNGP.h b/src/Particle/IntNGP.h index a2bf81723..8db494b3c 100644 --- a/src/Particle/IntNGP.h +++ b/src/Particle/IntNGP.h @@ -38,9 +38,9 @@ class IntNGP : public Interpolator { // gather/scatter functions // scatter particle data into Field using particle position and mesh - template - static void scatter(const FT& pdata, Field& f, const Vektor& ppos, - const M& mesh) { + template + static void scatter(const FT& pdata, Field& f, const Vektor& ppos, + const Mesh& mesh) { // find nearest-grid-point for particle position, store in NDIndex obj NDIndex ngp = FindNGP(mesh, ppos, CenteringTag()); // scatter data value to Field ... this assumes that the Field @@ -54,9 +54,9 @@ class IntNGP : public Interpolator { // scatter particle data into Field using particle position and mesh // and cache mesh information for reuse - template - static void scatter(const FT& pdata, Field& f, const Vektor& ppos, - const M& mesh, NDIndex& ngp) { + template + static void scatter(const FT& pdata, Field& f, const Vektor& ppos, + const Mesh& mesh, NDIndex& ngp) { // find nearest-grid-point for particle position, store in NDIndex obj ngp = FindNGP(mesh, ppos, CenteringTag()); // scatter data value to Field ... this assumes that the Field @@ -69,8 +69,8 @@ class IntNGP : public Interpolator { } // scatter particle data into Field using cached mesh information - template - static void scatter(const FT& pdata, Field& f, const NDIndex& ngp) { + template + static void scatter(const FT& pdata, Field& f, const NDIndex& ngp) { // scatter data value to Field ... this assumes that the Field // data point is local to this processor, if not an error will be printed. @@ -81,9 +81,9 @@ class IntNGP : public Interpolator { } // gather particle data from Field using particle position and mesh - template - static void gather(FT& pdata, const Field& f, const Vektor& ppos, - const M& mesh) { + template + static void gather(FT& pdata, const Field& f, const Vektor& ppos, + const Mesh& mesh) { // find nearest-grid-point for particle position, store in NDIndex obj NDIndex ngp = FindNGP(mesh, ppos, CenteringTag()); // gather Field value to particle data ... this assumes that the Field @@ -97,9 +97,9 @@ class IntNGP : public Interpolator { // gather particle data from Field using particle position and mesh // and cache mesh information for reuse - template - static void gather(FT& pdata, const Field& f, const Vektor& ppos, - const M& mesh, NDIndex& ngp) { + template + static void gather(FT& pdata, const Field& f, const Vektor& ppos, + const Mesh& mesh, NDIndex& ngp) { // find nearest-grid-point for particle position, store in NDIndex obj ngp = FindNGP(mesh, ppos, CenteringTag()); // gather Field value to particle data ... this assumes that the Field @@ -112,8 +112,8 @@ class IntNGP : public Interpolator { } // gather particle data from Field using cached mesh information - template - static void gather(FT& pdata, const Field& f, const NDIndex& ngp) { + template + static void gather(FT& pdata, const Field& f, const NDIndex& ngp) { // gather Field value to particle data ... this assumes that the Field // data point is local to this processor, if not an error will be printed. diff --git a/src/Solver/Electrostatics.h b/src/Solver/Electrostatics.h index d2ea3aa12..2ed346447 100644 --- a/src/Solver/Electrostatics.h +++ b/src/Solver/Electrostatics.h @@ -23,13 +23,12 @@ namespace ippl { - template , - class C = typename M::DefaultCentering> - class Electrostatics : public Solver { + template + class Electrostatics : public Solver { public: - using grad_type = Field, Dim, M, C>; - using lhs_type = typename Solver::lhs_type; - using rhs_type = typename Solver::rhs_type; + using grad_type = Field, Dim, Mesh, Cell>; + using lhs_type = typename Solver::lhs_type; + using rhs_type = typename Solver::rhs_type; /*! * Represents the types of fields that should @@ -46,13 +45,13 @@ namespace ippl { * desired output type defaults to solution only */ Electrostatics() - : Solver() + : Solver() , grad_mp(nullptr) { setDefaultParameters(); } Electrostatics(lhs_type& lhs, rhs_type& rhs) - : Solver(lhs, rhs) + : Solver(lhs, rhs) , grad_mp(nullptr) { setDefaultParameters(); } diff --git a/src/Solver/ElectrostaticsCG.h b/src/Solver/ElectrostaticsCG.h index cf26ef921..3c6bcaaf2 100644 --- a/src/Solver/ElectrostaticsCG.h +++ b/src/Solver/ElectrostaticsCG.h @@ -32,15 +32,14 @@ namespace ippl { return fun(arg); \ } - template , - class C = typename M::DefaultCentering> - class ElectrostaticsCG : public Electrostatics { + template + class ElectrostaticsCG : public Electrostatics { public: - using lhs_type = typename Solver::lhs_type; - using rhs_type = typename Solver::rhs_type; + using lhs_type = typename Solver::lhs_type; + using rhs_type = typename Solver::rhs_type; using OpRet = UnaryMinus>; - using algo = PCG; - using Base = Electrostatics; + using algo = PCG; + using Base = Electrostatics; ElectrostaticsCG() : Base() { diff --git a/src/Solver/FFTPeriodicPoissonSolver.h b/src/Solver/FFTPeriodicPoissonSolver.h index 748021ce0..ecd05a844 100644 --- a/src/Solver/FFTPeriodicPoissonSolver.h +++ b/src/Solver/FFTPeriodicPoissonSolver.h @@ -28,9 +28,8 @@ namespace ippl { - template , - class C = typename M::DefaultCentering> - class FFTPeriodicPoissonSolver : public Electrostatics { + template + class FFTPeriodicPoissonSolver : public Electrostatics { public: using Field_t = Field; using FFT_t = FFT; @@ -39,9 +38,9 @@ namespace ippl { using Layout_t = FieldLayout; using Vector_t = Vector; - using Base = Electrostatics; - using lhs_type = typename Solver::lhs_type; - using rhs_type = typename Solver::rhs_type; + using Base = Electrostatics; + using lhs_type = typename Solver::lhs_type; + using rhs_type = typename Solver::rhs_type; FFTPeriodicPoissonSolver() : Base() { diff --git a/src/Solver/FFTPeriodicPoissonSolver.hpp b/src/Solver/FFTPeriodicPoissonSolver.hpp index 09df65e92..1748632ff 100644 --- a/src/Solver/FFTPeriodicPoissonSolver.hpp +++ b/src/Solver/FFTPeriodicPoissonSolver.hpp @@ -19,14 +19,14 @@ namespace ippl { - template - void FFTPeriodicPoissonSolver::setRhs(rhs_type& rhs) { + template + void FFTPeriodicPoissonSolver::setRhs(rhs_type& rhs) { Base::setRhs(rhs); initialize(); } - template - void FFTPeriodicPoissonSolver::initialize() { + template + void FFTPeriodicPoissonSolver::initialize() { const Layout_t& layout_r = this->rhs_mp->getLayout(); domain_m = layout_r.getDomain(); @@ -56,8 +56,8 @@ namespace ippl { fft_mp = std::make_shared(layout_r, *layoutComplex_mp, this->params_m); } - template - void FFTPeriodicPoissonSolver::solve() { + template + void FFTPeriodicPoissonSolver::solve() { fft_mp->transform(1, *this->rhs_mp, fieldComplex_m); auto view = fieldComplex_m.getView(); diff --git a/src/Solver/FFTPoissonSolver.h b/src/Solver/FFTPoissonSolver.h index 84068caa3..ba48eab75 100644 --- a/src/Solver/FFTPoissonSolver.h +++ b/src/Solver/FFTPoissonSolver.h @@ -25,23 +25,22 @@ #include "Types/Vector.h" namespace ippl { - template , - class C = typename M::DefaultCentering> - class FFTPoissonSolver : public Electrostatics { + template + class FFTPoissonSolver : public Electrostatics { public: // types for LHS and RHS - using lhs_type = typename Solver::lhs_type; - using rhs_type = typename Solver::rhs_type; + using lhs_type = typename Solver::lhs_type; + using rhs_type = typename Solver::rhs_type; // type of output - using Base = Electrostatics; + using Base = Electrostatics; // define a type for a 3 dimensional field (e.g. charge density field) // define a type of Field with integers to be used for the helper Green's function // also define a type for the Fourier transformed complex valued fields - typedef Field Field_t; - typedef Field IField_t; - typedef Field, Dim, M> CxField_t; + typedef Field Field_t; + typedef Field IField_t; + typedef Field, Dim, Mesh, Cell> CxField_t; typedef Vector Vector_t; // define type for field layout @@ -110,11 +109,11 @@ namespace ippl { FieldLayout_t* layout_mp; // mesh and layout objects for rho2_m - std::unique_ptr mesh2_m; + std::unique_ptr mesh2_m; std::unique_ptr layout2_m; // mesh and layout objects for the Fourier transformed Complex fields - std::unique_ptr meshComplex_m; + std::unique_ptr meshComplex_m; std::unique_ptr layoutComplex_m; // domains for the various fields @@ -134,7 +133,7 @@ namespace ippl { std::unique_ptr> fft4n_m; - std::unique_ptr mesh4_m; + std::unique_ptr mesh4_m; std::unique_ptr layout4_m; NDIndex domain4_m; diff --git a/src/Solver/FFTPoissonSolver.hpp b/src/Solver/FFTPoissonSolver.hpp index 00316e5ef..b23e6e790 100644 --- a/src/Solver/FFTPoissonSolver.hpp +++ b/src/Solver/FFTPoissonSolver.hpp @@ -130,8 +130,8 @@ namespace ippl { ///////////////////////////////////////////////////////////////////////// // constructor and destructor - template - FFTPoissonSolver::FFTPoissonSolver(rhs_type& rhs, + template + FFTPoissonSolver::FFTPoissonSolver(rhs_type& rhs, ParameterList& fftparams, std::string alg) : mesh_mp(nullptr) @@ -160,8 +160,8 @@ namespace ippl { IpplTimings::stopTimer(initialize); } - template - FFTPoissonSolver::FFTPoissonSolver(lhs_type& lhs, rhs_type& rhs, + template + FFTPoissonSolver::FFTPoissonSolver(lhs_type& lhs, rhs_type& rhs, ParameterList& fftparams, std::string alg, int sol) : mesh_mp(nullptr) @@ -191,15 +191,15 @@ namespace ippl { IpplTimings::stopTimer(initialize); } - template - FFTPoissonSolver::~FFTPoissonSolver(){}; + template + FFTPoissonSolver::~FFTPoissonSolver(){}; ///////////////////////////////////////////////////////////////////////// // allows user to set gradient of phi = Efield instead of spectral // calculation of Efield (which uses FFTs) - template - void FFTPoissonSolver::setGradFD() { + template + void FFTPoissonSolver::setGradFD() { // get the output type (sol, grad, or sol & grad) const int out = this->params_m.template get("output_type"); @@ -215,8 +215,8 @@ namespace ippl { ///////////////////////////////////////////////////////////////////////// // initializeFields method, called in constructor - template - void FFTPoissonSolver::initializeFields() { + template + void FFTPoissonSolver::initializeFields() { // first check if valid algorithm choice if ((alg_m != "VICO") && (alg_m != "HOCKNEY") && (alg_m != "BIHARMONIC")) { throw IpplException( @@ -258,7 +258,7 @@ namespace ippl { } // create double sized mesh and layout objects using the previously defined domain2_m - mesh2_m = std::unique_ptr(new M(domain2_m, hr_m, origin)); + mesh2_m = std::unique_ptr(new M(domain2_m, hr_m, origin)); layout2_m = std::unique_ptr(new FieldLayout_t(domain2_m, decomp)); // create the domain for the transformed (complex) fields @@ -274,7 +274,7 @@ namespace ippl { } // create mesh and layout for the real to complex FFT transformed fields - meshComplex_m = std::unique_ptr(new M(domainComplex_m, hr_m, origin)); + meshComplex_m = std::unique_ptr(new M(domainComplex_m, hr_m, origin)); layoutComplex_m = std::unique_ptr(new FieldLayout_t(domainComplex_m, decomp)); @@ -305,7 +305,7 @@ namespace ippl { } // 4N grid - mesh4_m = std::unique_ptr(new M(domain4_m, hr_m, origin)); + mesh4_m = std::unique_ptr(new M(domain4_m, hr_m, origin)); layout4_m = std::unique_ptr(new FieldLayout_t(domain4_m, decomp)); // initialize fields @@ -424,8 +424,8 @@ namespace ippl { ///////////////////////////////////////////////////////////////////////// // compute electric potential by solving Poisson's eq given a field rho and mesh spacings hr - template - void FFTPoissonSolver::solve() { + template + void FFTPoissonSolver::solve() { // start a timer static IpplTimings::TimerRef solve = IpplTimings::getTimer("Solve"); IpplTimings::startTimer(solve); @@ -841,8 +841,8 @@ namespace ippl { //////////////////////////////////////////////////////////////////////// // calculate FFT of the Green's function - template - void FFTPoissonSolver::greensFunction() { + template + void FFTPoissonSolver::greensFunction() { const double pi = std::acos(-1.0); grn_mr = 0.0; @@ -1061,8 +1061,8 @@ namespace ippl { IpplTimings::stopTimer(fftg); }; - template - void FFTPoissonSolver::communicateVico( + template + void FFTPoissonSolver::communicateVico( Vector size, typename CxField_t::view_type view_g, const ippl::NDIndex ldom_g, const int nghost_g, typename Field_t::view_type view, const ippl::NDIndex ldom, const int nghost) { diff --git a/src/Solver/PCG.h b/src/Solver/PCG.h index e3ec6d2be..36a39cc67 100644 --- a/src/Solver/PCG.h +++ b/src/Solver/PCG.h @@ -24,10 +24,10 @@ namespace ippl { template , class C = typename M::DefaultCentering> + class Mesh, class Cell> class PCG : public SolverAlgorithm { public: - using Base = SolverAlgorithm; + using Base = SolverAlgorithm; using typename Base::lhs_type; using typename Base::rhs_type; using operator_type = std::function; diff --git a/src/Solver/Solver.h b/src/Solver/Solver.h index b3218ee0d..a48800c6e 100644 --- a/src/Solver/Solver.h +++ b/src/Solver/Solver.h @@ -24,12 +24,11 @@ namespace ippl { - template , - class C = typename M::DefaultCentering> + template class Solver { public: - using lhs_type = Field; - using rhs_type = Field; + using lhs_type = Field; + using rhs_type = Field; /*! * Default constructor diff --git a/src/Solver/SolverAlgorithm.h b/src/Solver/SolverAlgorithm.h index 36d6e11b0..47dfab46b 100644 --- a/src/Solver/SolverAlgorithm.h +++ b/src/Solver/SolverAlgorithm.h @@ -24,12 +24,11 @@ namespace ippl { - template , - class C = typename M::DefaultCentering> + template class SolverAlgorithm { public: - using lhs_type = Field; - using rhs_type = Field; + using lhs_type = Field; + using rhs_type = Field; /*! * Solve the problem described by Op(lhs) = rhs, where Op is an unspecified From b613c923670810a70d0f5d9d59e3f11448476175 Mon Sep 17 00:00:00 2001 From: Matthias Frey Date: Wed, 15 Mar 2023 15:46:55 +0000 Subject: [PATCH 02/12] fix some compilation issues --- src/Expression/IpplExpressions.h | 2 ++ src/Field/Field.h | 2 +- src/Field/Field.hpp | 6 ++-- src/Field/FieldOperations.hpp | 50 ++++++++++++++++---------------- 4 files changed, 31 insertions(+), 29 deletions(-) diff --git a/src/Expression/IpplExpressions.h b/src/Expression/IpplExpressions.h index 2986531d0..a782c45b1 100644 --- a/src/Expression/IpplExpressions.h +++ b/src/Expression/IpplExpressions.h @@ -20,6 +20,8 @@ #include +#include + namespace ippl { namespace detail { /*! diff --git a/src/Field/Field.h b/src/Field/Field.h index 2aabb0f41..83193be9a 100644 --- a/src/Field/Field.h +++ b/src/Field/Field.h @@ -35,7 +35,7 @@ namespace ippl { using Layout_t = FieldLayout; using BareField_t = BareField; using view_type = typename BareField_t::view_type; - using BConds_t = BConds; + using BConds_t = BConds; // A default constructor, which should be used only if the user calls the // 'initialize' function before doing anything else. There are no special diff --git a/src/Field/Field.hpp b/src/Field/Field.hpp index 8ed19b0a6..543b4c553 100644 --- a/src/Field/Field.hpp +++ b/src/Field/Field.hpp @@ -41,7 +41,7 @@ namespace ippl { : BareField(l, nghost) , mesh_m(&m) { for (unsigned int face = 0; face < 2 * Dim; ++face) { - bc_m[face] = std::make_shared>(face); + bc_m[face] = std::make_shared>(face); } } @@ -52,13 +52,13 @@ namespace ippl { BareField::initialize(l, nghost); mesh_m = &m; for (unsigned int face = 0; face < 2 * Dim; ++face) { - bc_m[face] = std::make_shared>(face); + bc_m[face] = std::make_shared>(face); } } template T Field::getVolumeIntegral() const { - typename M::value_type dV = mesh_m->getCellVolume(); + typename Mesh::value_type dV = mesh_m->getCellVolume(); return this->sum() * dV; } diff --git a/src/Field/FieldOperations.hpp b/src/Field/FieldOperations.hpp index bd7dc4cc4..da5589b0e 100644 --- a/src/Field/FieldOperations.hpp +++ b/src/Field/FieldOperations.hpp @@ -89,14 +89,14 @@ namespace ippl { template detail::meta_grad> grad(Field& u) { u.fillHalo(); - BConds& bcField = u.getFieldBC(); + BConds& bcField = u.getFieldBC(); bcField.apply(u); - M& mesh = u.get_mesh(); - typename M::vector_type xvector(0); + Mesh& mesh = u.get_mesh(); + typename Mesh::vector_type xvector(0); xvector[0] = 0.5 / mesh.getMeshSpacing(0); - typename M::vector_type yvector(0); + typename Mesh::vector_type yvector(0); yvector[1] = 0.5 / mesh.getMeshSpacing(1); - typename M::vector_type zvector(0); + typename Mesh::vector_type zvector(0); zvector[2] = 0.5 / mesh.getMeshSpacing(2); return detail::meta_grad>(u, xvector, yvector, zvector); } @@ -108,14 +108,14 @@ namespace ippl { template detail::meta_div> div(Field& u) { u.fillHalo(); - BConds& bcField = u.getFieldBC(); + BConds& bcField = u.getFieldBC(); bcField.apply(u); - M& mesh = u.get_mesh(); - typename M::vector_type xvector(0); + Mesh& mesh = u.get_mesh(); + typename Mesh::vector_type xvector(0); xvector[0] = 0.5 / mesh.getMeshSpacing(0); - typename M::vector_type yvector(0); + typename Mesh::vector_type yvector(0); yvector[1] = 0.5 / mesh.getMeshSpacing(1); - typename M::vector_type zvector(0); + typename Mesh::vector_type zvector(0); zvector[2] = 0.5 / mesh.getMeshSpacing(2); return detail::meta_div>(u, xvector, yvector, zvector); } @@ -127,10 +127,10 @@ namespace ippl { template detail::meta_laplace> laplace(Field& u) { u.fillHalo(); - BConds& bcField = u.getFieldBC(); + BConds& bcField = u.getFieldBC(); bcField.apply(u); - M& mesh = u.get_mesh(); - typename M::vector_type hvector(0); + Mesh& mesh = u.get_mesh(); + typename Mesh::vector_type hvector(0); hvector[0] = 1.0 / std::pow(mesh.getMeshSpacing(0), 2); hvector[1] = 1.0 / std::pow(mesh.getMeshSpacing(1), 2); hvector[2] = 1.0 / std::pow(mesh.getMeshSpacing(2), 2); @@ -144,16 +144,16 @@ namespace ippl { template detail::meta_curl> curl(Field& u) { u.fillHalo(); - BConds& bcField = u.getFieldBC(); + BConds& bcField = u.getFieldBC(); bcField.apply(u); - M& mesh = u.get_mesh(); - typename M::vector_type xvector(0); + Mesh& mesh = u.get_mesh(); + typename Mesh::vector_type xvector(0); xvector[0] = 1.0; - typename M::vector_type yvector(0); + typename Mesh::vector_type yvector(0); yvector[1] = 1.0; - typename M::vector_type zvector(0); + typename Mesh::vector_type zvector(0); zvector[2] = 1.0; - typename M::vector_type hvector(0); + typename Mesh::vector_type hvector(0); hvector = mesh.getMeshSpacing(); return detail::meta_curl>(u, xvector, yvector, zvector, hvector); } @@ -165,16 +165,16 @@ namespace ippl { template detail::meta_hess> hess(Field& u) { u.fillHalo(); - BConds& bcField = u.getFieldBC(); + BConds& bcField = u.getFieldBC(); bcField.apply(u); - M& mesh = u.get_mesh(); - typename M::vector_type xvector(0); + Mesh& mesh = u.get_mesh(); + typename Mesh::vector_type xvector(0); xvector[0] = 1.0; - typename M::vector_type yvector(0); + typename Mesh::vector_type yvector(0); yvector[1] = 1.0; - typename M::vector_type zvector(0); + typename Mesh::vector_type zvector(0); zvector[2] = 1.0; - typename M::vector_type hvector(0); + typename Mesh::vector_type hvector(0); hvector = mesh.getMeshSpacing(); return detail::meta_hess>(u, xvector, yvector, zvector, hvector); } From c3600505aa1677678ef1df59c40af2ef45b312f2 Mon Sep 17 00:00:00 2001 From: Matthias Frey Date: Wed, 15 Mar 2023 17:37:21 +0000 Subject: [PATCH 03/12] fix templates in FFT class --- src/FFT/FFT.h | 28 ++++++++--------- src/FFT/FFT.hpp | 77 +++++++++++++++++++++++---------------------- src/Field/Field.hpp | 4 +-- 3 files changed, 55 insertions(+), 54 deletions(-) diff --git a/src/FFT/FFT.h b/src/FFT/FFT.h index f7b626ac7..c6a8d2e40 100644 --- a/src/FFT/FFT.h +++ b/src/FFT/FFT.h @@ -118,18 +118,18 @@ namespace ippl { /** Non-specialized FFT class. We specialize based on Transform tag class */ - template + template class FFT {}; /** complex-to-complex FFT class */ - template - class FFT { + template + class FFT { public: typedef FieldLayout Layout_t; typedef Kokkos::complex Complex_t; - typedef Field ComplexField_t; + typedef Field ComplexField_t; using heffteBackend = typename detail::HeffteBackendType::backend; using workspace_t = @@ -164,18 +164,18 @@ namespace ippl { /** real-to-complex FFT class */ - template - class FFT { + template + class FFT { public: typedef FieldLayout Layout_t; - typedef Field RealField_t; + typedef Field RealField_t; using heffteBackend = typename detail::HeffteBackendType::backend; typedef Kokkos::complex Complex_t; using workspace_t = typename heffte::fft3d_r2c::template buffer_container; - typedef Field ComplexField_t; + typedef Field ComplexField_t; /** Create a new FFT object with the layout for the input and output Fields * and parameters for heffte. @@ -208,11 +208,11 @@ namespace ippl { /** Sine transform class */ - template - class FFT { + template + class FFT { public: typedef FieldLayout Layout_t; - typedef Field Field_t; + typedef Field Field_t; using heffteBackend = typename detail::HeffteBackendType::backendSine; using workspace_t = typename heffte::fft3d::template buffer_container; @@ -243,11 +243,11 @@ namespace ippl { /** Cosine transform class */ - template - class FFT { + template + class FFT { public: typedef FieldLayout Layout_t; - typedef Field Field_t; + typedef Field Field_t; using heffteBackend = typename detail::HeffteBackendType::backendCos; using workspace_t = typename heffte::fft3d::template buffer_container; diff --git a/src/FFT/FFT.hpp b/src/FFT/FFT.hpp index 38815f2f6..fe89a167b 100644 --- a/src/FFT/FFT.hpp +++ b/src/FFT/FFT.hpp @@ -44,8 +44,8 @@ namespace ippl { given layout and heffte parameters. */ - template - FFT::FFT(const Layout_t& layout, const ParameterList& params) { + template + FFT::FFT(const Layout_t& layout, const ParameterList& params) { /** * Heffte requires to pass a 3D array even for 2D and * 1D FFTs we just have to make the length in other @@ -74,10 +74,10 @@ namespace ippl { /** setup performs the initialization necessary. */ - template - void FFT::setup(const std::array& low, - const std::array& high, - const ParameterList& params) { + template + void FFT::setup(const std::array& low, + const std::array& high, + const ParameterList& params) { heffte::box3d inbox = {low, high}; heffte::box3d outbox = {low, high}; @@ -116,9 +116,9 @@ namespace ippl { workspace_m = workspace_t(heffte_m->size_workspace()); } - template - void FFT::transform( - int direction, typename FFT::ComplexField_t& f) { + template + void FFT::transform( + int direction, typename FFT::ComplexField_t& f) { auto fview = f.getView(); const int nghost = f.getNghost(); @@ -175,9 +175,10 @@ namespace ippl { *layouts and heffte parameters. */ - template - FFT::FFT(const Layout_t& layoutInput, const Layout_t& layoutOutput, - const ParameterList& params) { + template + FFT::FFT(const Layout_t& layoutInput, + const Layout_t& layoutOutput, + const ParameterList& params) { /** * Heffte requires to pass a 3D array even for 2D and * 1D FFTs we just have to make the length in other @@ -215,12 +216,12 @@ namespace ippl { /** setup performs the initialization. */ - template - void FFT::setup(const std::array& lowInput, - const std::array& highInput, - const std::array& lowOutput, - const std::array& highOutput, - const ParameterList& params) { + template + void FFT::setup(const std::array& lowInput, + const std::array& highInput, + const std::array& lowOutput, + const std::array& highOutput, + const ParameterList& params) { heffte::box3d inbox = {lowInput, highInput}; heffte::box3d outbox = {lowOutput, highOutput}; @@ -259,10 +260,10 @@ namespace ippl { workspace_m = workspace_t(heffte_m->size_workspace()); } - template - void FFT::transform( - int direction, typename FFT::RealField_t& f, - typename FFT::ComplexField_t& g) { + template + void FFT::transform( + int direction, typename FFT::RealField_t& f, + typename FFT::ComplexField_t& g) { auto fview = f.getView(); auto gview = g.getView(); const int nghostf = f.getNghost(); @@ -342,8 +343,8 @@ namespace ippl { given layout and heffte parameters. */ - template - FFT::FFT(const Layout_t& layout, const ParameterList& params) { + template + FFT::FFT(const Layout_t& layout, const ParameterList& params) { /** * Heffte requires to pass a 3D array even for 2D and * 1D FFTs we just have to make the length in other @@ -372,8 +373,8 @@ namespace ippl { /** setup performs the initialization necessary. */ - template - void FFT::setup(const std::array& low, + template + void FFT::setup(const std::array& low, const std::array& high, const ParameterList& params) { heffte::box3d inbox = {low, high}; @@ -413,9 +414,9 @@ namespace ippl { workspace_m = workspace_t(heffte_m->size_workspace()); } - template - void FFT::transform( - int direction, typename FFT::Field_t& f) { + template + void FFT::transform( + int direction, typename FFT::Field_t& f) { auto fview = f.getView(); const int nghost = f.getNghost(); @@ -470,8 +471,8 @@ namespace ippl { given layout and heffte parameters. */ - template - FFT::FFT(const Layout_t& layout, const ParameterList& params) { + template + FFT::FFT(const Layout_t& layout, const ParameterList& params) { /** * Heffte requires to pass a 3D array even for 2D and * 1D FFTs we just have to make the length in other @@ -500,10 +501,10 @@ namespace ippl { /** setup performs the initialization necessary. */ - template - void FFT::setup(const std::array& low, - const std::array& high, - const ParameterList& params) { + template + void FFT::setup(const std::array& low, + const std::array& high, + const ParameterList& params) { heffte::box3d inbox = {low, high}; heffte::box3d outbox = {low, high}; @@ -541,9 +542,9 @@ namespace ippl { workspace_m = workspace_t(heffte_m->size_workspace()); } - template - void FFT::transform( - int direction, typename FFT::Field_t& f) { + template + void FFT::Field_t& f) { auto fview = f.getView(); const int nghost = f.getNghost(); diff --git a/src/Field/Field.hpp b/src/Field/Field.hpp index 543b4c553..dd03c9e61 100644 --- a/src/Field/Field.hpp +++ b/src/Field/Field.hpp @@ -41,7 +41,7 @@ namespace ippl { : BareField(l, nghost) , mesh_m(&m) { for (unsigned int face = 0; face < 2 * Dim; ++face) { - bc_m[face] = std::make_shared>(face); + bc_m[face] = std::make_shared >(face); } } @@ -52,7 +52,7 @@ namespace ippl { BareField::initialize(l, nghost); mesh_m = &m; for (unsigned int face = 0; face < 2 * Dim; ++face) { - bc_m[face] = std::make_shared>(face); + bc_m[face] = std::make_shared >(face); } } From cc953a44a25273cb7b13909a11b310006fe15503 Mon Sep 17 00:00:00 2001 From: Matthias Frey Date: Thu, 16 Mar 2023 09:52:33 +0000 Subject: [PATCH 04/12] fix compilation errors --- .../OrthogonalRecursiveBisection.h | 6 +-- .../OrthogonalRecursiveBisection.hpp | 48 +++++++++---------- src/FFT/FFT.hpp | 20 ++++---- src/Field/Field.h | 1 + src/Field/FieldOperations.hpp | 20 ++++---- src/IpplCore.h | 1 - 6 files changed, 48 insertions(+), 48 deletions(-) diff --git a/src/Decomposition/OrthogonalRecursiveBisection.h b/src/Decomposition/OrthogonalRecursiveBisection.h index ca7cef6c9..9f04a7eb1 100644 --- a/src/Decomposition/OrthogonalRecursiveBisection.h +++ b/src/Decomposition/OrthogonalRecursiveBisection.h @@ -37,13 +37,13 @@ namespace ippl { * @tparam Dim dimension * @tparam M mesh */ - template + template class OrthogonalRecursiveBisection { public: using view_type = typename detail::ViewType::view_type; // Weight for reduction - Field bf_m; + Field bf_m; /*! * Initialize member field with mesh and field layout @@ -52,7 +52,7 @@ namespace ippl { * @param rho Density field */ void initialize(FieldLayout& fl, UniformCartesian& mesh, - const Field& rho); + const Field& rho); /*! * Performs scatter operation of particle positions in field (weights) and diff --git a/src/Decomposition/OrthogonalRecursiveBisection.hpp b/src/Decomposition/OrthogonalRecursiveBisection.hpp index 48ec6193c..9d167a095 100644 --- a/src/Decomposition/OrthogonalRecursiveBisection.hpp +++ b/src/Decomposition/OrthogonalRecursiveBisection.hpp @@ -1,16 +1,16 @@ #include "Utility/IpplTimings.h" namespace ippl { - template - void OrthogonalRecursiveBisection::initialize(FieldLayout& fl, - UniformCartesian& mesh, - const Field& rho) { + template + void OrthogonalRecursiveBisection::initialize(FieldLayout& fl, + UniformCartesian& mesh, + const Field& rho) { bf_m.initialize(mesh, fl); bf_m = rho; } - template - bool OrthogonalRecursiveBisection::binaryRepartition( + template + bool OrthogonalRecursiveBisection::binaryRepartition( const ParticleAttrib>& R, FieldLayout& fl, const bool& isFirstRepartition) { // Timings @@ -119,8 +119,8 @@ namespace ippl { return true; } - template - int OrthogonalRecursiveBisection::findCutAxis(NDIndex& dom) { + template + int OrthogonalRecursiveBisection::findCutAxis(NDIndex& dom) { int cutAxis = 0; unsigned int maxLength = 0; @@ -136,8 +136,8 @@ namespace ippl { return cutAxis; } - template - void OrthogonalRecursiveBisection::perpendicularReduction( + template + void OrthogonalRecursiveBisection::perpendicularReduction( std::vector& rankWeights, unsigned int cutAxis, NDIndex& dom) { // Check if domains overlap, if not no need for reduction NDIndex lDom = bf_m.getOwned(); @@ -214,8 +214,8 @@ namespace ippl { } } - template - int OrthogonalRecursiveBisection::findMedian(std::vector& w) { + template + int OrthogonalRecursiveBisection::findMedian(std::vector& w) { // Special case when array must be cut in half in order to not have planes if (w.size() == 4) return 1; @@ -250,10 +250,10 @@ namespace ippl { return w.size() - 3; } - template - void OrthogonalRecursiveBisection::cutDomain(std::vector>& domains, - std::vector& procs, int it, - int cutAxis, int median) { + template + void OrthogonalRecursiveBisection::cutDomain(std::vector>& domains, + std::vector& procs, int it, + int cutAxis, int median) { // Cut domains[it] in half at median along cutAxis NDIndex leftDom, rightDom; domains[it].split(leftDom, rightDom, cutAxis, median + domains[it][cutAxis].first()); @@ -266,19 +266,19 @@ namespace ippl { procs.insert(procs.begin() + it + 1, 1, temp - procs[it]); } - template - void OrthogonalRecursiveBisection::scatterR( + template + void OrthogonalRecursiveBisection::scatterR( const ParticleAttrib>& r) { - using vector_type = typename M::vector_type; + using vector_type = typename Mesh::vector_type; // Reset local field bf_m = 0.0; // Get local data - typename Field::view_type view = bf_m.getView(); - const M& mesh = bf_m.get_mesh(); - const FieldLayout& layout = bf_m.getLayout(); - const NDIndex& lDom = layout.getLocalNDIndex(); - const int nghost = bf_m.getNghost(); + typename Field::view_type view = bf_m.getView(); + const Mesh& mesh = bf_m.get_mesh(); + const FieldLayout& layout = bf_m.getLayout(); + const NDIndex& lDom = layout.getLocalNDIndex(); + const int nghost = bf_m.getNghost(); // Get spacings const vector_type& dx = mesh.getMeshSpacing(); diff --git a/src/FFT/FFT.hpp b/src/FFT/FFT.hpp index fe89a167b..5aeeb4458 100644 --- a/src/FFT/FFT.hpp +++ b/src/FFT/FFT.hpp @@ -217,11 +217,11 @@ namespace ippl { setup performs the initialization. */ template - void FFT::setup(const std::array& lowInput, - const std::array& highInput, - const std::array& lowOutput, - const std::array& highOutput, - const ParameterList& params) { + void FFT::setup(const std::array& lowInput, + const std::array& highInput, + const std::array& lowOutput, + const std::array& highOutput, + const ParameterList& params) { heffte::box3d inbox = {lowInput, highInput}; heffte::box3d outbox = {lowOutput, highOutput}; @@ -472,7 +472,7 @@ namespace ippl { */ template - FFT::FFT(const Layout_t& layout, const ParameterList& params) { + FFT::FFT(const Layout_t& layout, const ParameterList& params) { /** * Heffte requires to pass a 3D array even for 2D and * 1D FFTs we just have to make the length in other @@ -502,9 +502,9 @@ namespace ippl { setup performs the initialization necessary. */ template - void FFT::setup(const std::array& low, - const std::array& high, - const ParameterList& params) { + void FFT::setup(const std::array& low, + const std::array& high, + const ParameterList& params) { heffte::box3d inbox = {low, high}; heffte::box3d outbox = {low, high}; @@ -543,7 +543,7 @@ namespace ippl { } template - void FFT::transform( int direction, typename FFT::Field_t& f) { auto fview = f.getView(); const int nghost = f.getNghost(); diff --git a/src/Field/Field.h b/src/Field/Field.h index 83193be9a..6ec90bf21 100644 --- a/src/Field/Field.h +++ b/src/Field/Field.h @@ -19,6 +19,7 @@ #define IPPL_FIELD_H #include "Field/BConds.h" +#include "Field/BcTypes.h" #include "Field/BareField.h" #include "Meshes/UniformCartesian.h" diff --git a/src/Field/FieldOperations.hpp b/src/Field/FieldOperations.hpp index da5589b0e..aa168fbb4 100644 --- a/src/Field/FieldOperations.hpp +++ b/src/Field/FieldOperations.hpp @@ -87,7 +87,7 @@ namespace ippl { * @param u field */ template - detail::meta_grad> grad(Field& u) { + detail::meta_grad > grad(Field& u) { u.fillHalo(); BConds& bcField = u.getFieldBC(); bcField.apply(u); @@ -98,7 +98,7 @@ namespace ippl { yvector[1] = 0.5 / mesh.getMeshSpacing(1); typename Mesh::vector_type zvector(0); zvector[2] = 0.5 / mesh.getMeshSpacing(2); - return detail::meta_grad>(u, xvector, yvector, zvector); + return detail::meta_grad >(u, xvector, yvector, zvector); } /*! @@ -106,7 +106,7 @@ namespace ippl { * @param u field */ template - detail::meta_div> div(Field& u) { + detail::meta_div > div(Field& u) { u.fillHalo(); BConds& bcField = u.getFieldBC(); bcField.apply(u); @@ -117,7 +117,7 @@ namespace ippl { yvector[1] = 0.5 / mesh.getMeshSpacing(1); typename Mesh::vector_type zvector(0); zvector[2] = 0.5 / mesh.getMeshSpacing(2); - return detail::meta_div>(u, xvector, yvector, zvector); + return detail::meta_div >(u, xvector, yvector, zvector); } /*! @@ -125,7 +125,7 @@ namespace ippl { * @param u field */ template - detail::meta_laplace> laplace(Field& u) { + detail::meta_laplace > laplace(Field& u) { u.fillHalo(); BConds& bcField = u.getFieldBC(); bcField.apply(u); @@ -134,7 +134,7 @@ namespace ippl { hvector[0] = 1.0 / std::pow(mesh.getMeshSpacing(0), 2); hvector[1] = 1.0 / std::pow(mesh.getMeshSpacing(1), 2); hvector[2] = 1.0 / std::pow(mesh.getMeshSpacing(2), 2); - return detail::meta_laplace>(u, hvector); + return detail::meta_laplace >(u, hvector); } /*! @@ -142,7 +142,7 @@ namespace ippl { * @param u field */ template - detail::meta_curl> curl(Field& u) { + detail::meta_curl > curl(Field& u) { u.fillHalo(); BConds& bcField = u.getFieldBC(); bcField.apply(u); @@ -155,7 +155,7 @@ namespace ippl { zvector[2] = 1.0; typename Mesh::vector_type hvector(0); hvector = mesh.getMeshSpacing(); - return detail::meta_curl>(u, xvector, yvector, zvector, hvector); + return detail::meta_curl >(u, xvector, yvector, zvector, hvector); } /*! @@ -163,7 +163,7 @@ namespace ippl { * @param u field */ template - detail::meta_hess> hess(Field& u) { + detail::meta_hess > hess(Field& u) { u.fillHalo(); BConds& bcField = u.getFieldBC(); bcField.apply(u); @@ -176,6 +176,6 @@ namespace ippl { zvector[2] = 1.0; typename Mesh::vector_type hvector(0); hvector = mesh.getMeshSpacing(); - return detail::meta_hess>(u, xvector, yvector, zvector, hvector); + return detail::meta_hess >(u, xvector, yvector, zvector, hvector); } } // namespace ippl diff --git a/src/IpplCore.h b/src/IpplCore.h index c00203d71..c159e052b 100644 --- a/src/IpplCore.h +++ b/src/IpplCore.h @@ -18,7 +18,6 @@ #ifndef IPPL_CORE_H #define IPPL_CORE_H -#include "Field/BConds.h" #include "Field/BareField.h" #include "Field/Field.h" From bcaa58172a088e291cb72f574bd7339b38d68e92 Mon Sep 17 00:00:00 2001 From: Matthias Frey Date: Thu, 16 Mar 2023 15:42:17 +0000 Subject: [PATCH 05/12] Cell to Centering --- .../OrthogonalRecursiveBisection.h | 6 +- .../OrthogonalRecursiveBisection.hpp | 32 +++++----- src/FFT/FFT.h | 28 ++++----- src/FFT/FFT.hpp | 58 +++++++++---------- src/Field/BConds.h | 20 +++---- src/Field/BConds.hpp | 16 ++--- src/Field/BcTypes.h | 56 +++++++++--------- src/Field/BcTypes.hpp | 40 ++++++------- src/Field/Field.h | 4 +- src/Field/Field.hpp | 32 +++++----- src/Field/FieldOperations.hpp | 48 +++++++-------- src/Particle/IntNGP.h | 24 ++++---- src/Particle/ParticleBalancer.hpp | 2 +- src/Solver/Electrostatics.h | 14 ++--- src/Solver/ElectrostaticsCG.h | 12 ++-- src/Solver/FFTPeriodicPoissonSolver.h | 10 ++-- src/Solver/FFTPeriodicPoissonSolver.hpp | 12 ++-- src/Solver/FFTPoissonSolver.h | 16 ++--- src/Solver/FFTPoissonSolver.hpp | 32 +++++----- src/Solver/PCG.h | 4 +- src/Solver/Solver.h | 6 +- src/Solver/SolverAlgorithm.h | 6 +- 22 files changed, 239 insertions(+), 239 deletions(-) diff --git a/src/Decomposition/OrthogonalRecursiveBisection.h b/src/Decomposition/OrthogonalRecursiveBisection.h index 9f04a7eb1..aa2582ebe 100644 --- a/src/Decomposition/OrthogonalRecursiveBisection.h +++ b/src/Decomposition/OrthogonalRecursiveBisection.h @@ -37,13 +37,13 @@ namespace ippl { * @tparam Dim dimension * @tparam M mesh */ - template + template class OrthogonalRecursiveBisection { public: using view_type = typename detail::ViewType::view_type; // Weight for reduction - Field bf_m; + Field bf_m; /*! * Initialize member field with mesh and field layout @@ -52,7 +52,7 @@ namespace ippl { * @param rho Density field */ void initialize(FieldLayout& fl, UniformCartesian& mesh, - const Field& rho); + const Field& rho); /*! * Performs scatter operation of particle positions in field (weights) and diff --git a/src/Decomposition/OrthogonalRecursiveBisection.hpp b/src/Decomposition/OrthogonalRecursiveBisection.hpp index 9d167a095..2baa8770b 100644 --- a/src/Decomposition/OrthogonalRecursiveBisection.hpp +++ b/src/Decomposition/OrthogonalRecursiveBisection.hpp @@ -1,16 +1,16 @@ #include "Utility/IpplTimings.h" namespace ippl { - template - void OrthogonalRecursiveBisection::initialize(FieldLayout& fl, + template + void OrthogonalRecursiveBisection::initialize(FieldLayout& fl, UniformCartesian& mesh, - const Field& rho) { + const Field& rho) { bf_m.initialize(mesh, fl); bf_m = rho; } - template - bool OrthogonalRecursiveBisection::binaryRepartition( + template + bool OrthogonalRecursiveBisection::binaryRepartition( const ParticleAttrib>& R, FieldLayout& fl, const bool& isFirstRepartition) { // Timings @@ -119,8 +119,8 @@ namespace ippl { return true; } - template - int OrthogonalRecursiveBisection::findCutAxis(NDIndex& dom) { + template + int OrthogonalRecursiveBisection::findCutAxis(NDIndex& dom) { int cutAxis = 0; unsigned int maxLength = 0; @@ -136,8 +136,8 @@ namespace ippl { return cutAxis; } - template - void OrthogonalRecursiveBisection::perpendicularReduction( + template + void OrthogonalRecursiveBisection::perpendicularReduction( std::vector& rankWeights, unsigned int cutAxis, NDIndex& dom) { // Check if domains overlap, if not no need for reduction NDIndex lDom = bf_m.getOwned(); @@ -214,8 +214,8 @@ namespace ippl { } } - template - int OrthogonalRecursiveBisection::findMedian(std::vector& w) { + template + int OrthogonalRecursiveBisection::findMedian(std::vector& w) { // Special case when array must be cut in half in order to not have planes if (w.size() == 4) return 1; @@ -250,8 +250,8 @@ namespace ippl { return w.size() - 3; } - template - void OrthogonalRecursiveBisection::cutDomain(std::vector>& domains, + template + void OrthogonalRecursiveBisection::cutDomain(std::vector>& domains, std::vector& procs, int it, int cutAxis, int median) { // Cut domains[it] in half at median along cutAxis @@ -266,15 +266,15 @@ namespace ippl { procs.insert(procs.begin() + it + 1, 1, temp - procs[it]); } - template - void OrthogonalRecursiveBisection::scatterR( + template + void OrthogonalRecursiveBisection::scatterR( const ParticleAttrib>& r) { using vector_type = typename Mesh::vector_type; // Reset local field bf_m = 0.0; // Get local data - typename Field::view_type view = bf_m.getView(); + typename Field::view_type view = bf_m.getView(); const Mesh& mesh = bf_m.get_mesh(); const FieldLayout& layout = bf_m.getLayout(); const NDIndex& lDom = layout.getLocalNDIndex(); diff --git a/src/FFT/FFT.h b/src/FFT/FFT.h index c6a8d2e40..f38f32c45 100644 --- a/src/FFT/FFT.h +++ b/src/FFT/FFT.h @@ -118,18 +118,18 @@ namespace ippl { /** Non-specialized FFT class. We specialize based on Transform tag class */ - template + template class FFT {}; /** complex-to-complex FFT class */ - template - class FFT { + template + class FFT { public: typedef FieldLayout Layout_t; typedef Kokkos::complex Complex_t; - typedef Field ComplexField_t; + typedef Field ComplexField_t; using heffteBackend = typename detail::HeffteBackendType::backend; using workspace_t = @@ -164,18 +164,18 @@ namespace ippl { /** real-to-complex FFT class */ - template - class FFT { + template + class FFT { public: typedef FieldLayout Layout_t; - typedef Field RealField_t; + typedef Field RealField_t; using heffteBackend = typename detail::HeffteBackendType::backend; typedef Kokkos::complex Complex_t; using workspace_t = typename heffte::fft3d_r2c::template buffer_container; - typedef Field ComplexField_t; + typedef Field ComplexField_t; /** Create a new FFT object with the layout for the input and output Fields * and parameters for heffte. @@ -208,11 +208,11 @@ namespace ippl { /** Sine transform class */ - template - class FFT { + template + class FFT { public: typedef FieldLayout Layout_t; - typedef Field Field_t; + typedef Field Field_t; using heffteBackend = typename detail::HeffteBackendType::backendSine; using workspace_t = typename heffte::fft3d::template buffer_container; @@ -243,11 +243,11 @@ namespace ippl { /** Cosine transform class */ - template - class FFT { + template + class FFT { public: typedef FieldLayout Layout_t; - typedef Field Field_t; + typedef Field Field_t; using heffteBackend = typename detail::HeffteBackendType::backendCos; using workspace_t = typename heffte::fft3d::template buffer_container; diff --git a/src/FFT/FFT.hpp b/src/FFT/FFT.hpp index 5aeeb4458..de71dc7aa 100644 --- a/src/FFT/FFT.hpp +++ b/src/FFT/FFT.hpp @@ -44,8 +44,8 @@ namespace ippl { given layout and heffte parameters. */ - template - FFT::FFT(const Layout_t& layout, const ParameterList& params) { + template + FFT::FFT(const Layout_t& layout, const ParameterList& params) { /** * Heffte requires to pass a 3D array even for 2D and * 1D FFTs we just have to make the length in other @@ -74,8 +74,8 @@ namespace ippl { /** setup performs the initialization necessary. */ - template - void FFT::setup(const std::array& low, + template + void FFT::setup(const std::array& low, const std::array& high, const ParameterList& params) { heffte::box3d inbox = {low, high}; @@ -116,9 +116,9 @@ namespace ippl { workspace_m = workspace_t(heffte_m->size_workspace()); } - template - void FFT::transform( - int direction, typename FFT::ComplexField_t& f) { + template + void FFT::transform( + int direction, typename FFT::ComplexField_t& f) { auto fview = f.getView(); const int nghost = f.getNghost(); @@ -175,8 +175,8 @@ namespace ippl { *layouts and heffte parameters. */ - template - FFT::FFT(const Layout_t& layoutInput, + template + FFT::FFT(const Layout_t& layoutInput, const Layout_t& layoutOutput, const ParameterList& params) { /** @@ -216,8 +216,8 @@ namespace ippl { /** setup performs the initialization. */ - template - void FFT::setup(const std::array& lowInput, + template + void FFT::setup(const std::array& lowInput, const std::array& highInput, const std::array& lowOutput, const std::array& highOutput, @@ -260,10 +260,10 @@ namespace ippl { workspace_m = workspace_t(heffte_m->size_workspace()); } - template - void FFT::transform( - int direction, typename FFT::RealField_t& f, - typename FFT::ComplexField_t& g) { + template + void FFT::transform( + int direction, typename FFT::RealField_t& f, + typename FFT::ComplexField_t& g) { auto fview = f.getView(); auto gview = g.getView(); const int nghostf = f.getNghost(); @@ -343,8 +343,8 @@ namespace ippl { given layout and heffte parameters. */ - template - FFT::FFT(const Layout_t& layout, const ParameterList& params) { + template + FFT::FFT(const Layout_t& layout, const ParameterList& params) { /** * Heffte requires to pass a 3D array even for 2D and * 1D FFTs we just have to make the length in other @@ -373,8 +373,8 @@ namespace ippl { /** setup performs the initialization necessary. */ - template - void FFT::setup(const std::array& low, + template + void FFT::setup(const std::array& low, const std::array& high, const ParameterList& params) { heffte::box3d inbox = {low, high}; @@ -414,9 +414,9 @@ namespace ippl { workspace_m = workspace_t(heffte_m->size_workspace()); } - template - void FFT::transform( - int direction, typename FFT::Field_t& f) { + template + void FFT::transform( + int direction, typename FFT::Field_t& f) { auto fview = f.getView(); const int nghost = f.getNghost(); @@ -471,8 +471,8 @@ namespace ippl { given layout and heffte parameters. */ - template - FFT::FFT(const Layout_t& layout, const ParameterList& params) { + template + FFT::FFT(const Layout_t& layout, const ParameterList& params) { /** * Heffte requires to pass a 3D array even for 2D and * 1D FFTs we just have to make the length in other @@ -501,8 +501,8 @@ namespace ippl { /** setup performs the initialization necessary. */ - template - void FFT::setup(const std::array& low, + template + void FFT::setup(const std::array& low, const std::array& high, const ParameterList& params) { heffte::box3d inbox = {low, high}; @@ -542,9 +542,9 @@ namespace ippl { workspace_m = workspace_t(heffte_m->size_workspace()); } - template - void FFT::transform( - int direction, typename FFT::Field_t& f) { + template + void FFT::transform( + int direction, typename FFT::Field_t& f) { auto fview = f.getView(); const int nghost = f.getNghost(); diff --git a/src/Field/BConds.h b/src/Field/BConds.h index 623fc1d79..35206ef19 100644 --- a/src/Field/BConds.h +++ b/src/Field/BConds.h @@ -28,19 +28,19 @@ #include namespace ippl { - template + template class Field; - template + template class BConds; - template - std::ostream& operator<<(std::ostream&, const BConds&); + template + std::ostream& operator<<(std::ostream&, const BConds&); - template + template class BConds { public: - using bc_type = detail::BCondBase; + using bc_type = detail::BCondBase; using container = std::array, 2 * Dim>; using iterator = typename container::iterator; using const_iterator = typename container::const_iterator; @@ -48,8 +48,8 @@ namespace ippl { BConds() = default; ~BConds() = default; - void findBCNeighbors(Field& field); - void apply(Field& field); + void findBCNeighbors(Field& field); + void apply(Field& field); bool changesPhysicalCells() const; virtual void write(std::ostream&) const; @@ -62,8 +62,8 @@ namespace ippl { container bc_m; }; - template - inline std::ostream& operator<<(std::ostream& os, const BConds& bc) { + template + inline std::ostream& operator<<(std::ostream& os, const BConds& bc) { bc.write(os); return os; } diff --git a/src/Field/BConds.hpp b/src/Field/BConds.hpp index 540bbc366..eff351a19 100644 --- a/src/Field/BConds.hpp +++ b/src/Field/BConds.hpp @@ -19,8 +19,8 @@ // along with IPPL. If not, see . // namespace ippl { - template - void BConds::write(std::ostream& os) const { + template + void BConds::write(std::ostream& os) const { os << "BConds: (" << std::endl; const_iterator it = bc_m.begin(); for (; it != bc_m.end() - 1; ++it) { @@ -31,8 +31,8 @@ namespace ippl { os << std::endl << ")"; } - template - void BConds::findBCNeighbors(Field& field) { + template + void BConds::findBCNeighbors(Field& field) { for (iterator it = bc_m.begin(); it != bc_m.end(); ++it) { (*it)->findBCNeighbors(field); } @@ -40,8 +40,8 @@ namespace ippl { Ippl::Comm->barrier(); } - template - void BConds::apply(Field& field) { + template + void BConds::apply(Field& field) { for (iterator it = bc_m.begin(); it != bc_m.end(); ++it) { (*it)->apply(field); } @@ -49,8 +49,8 @@ namespace ippl { Ippl::Comm->barrier(); } - template - bool BConds::changesPhysicalCells() const { + template + bool BConds::changesPhysicalCells() const { bool doesChange = false; for (const_iterator it = bc_m.begin(); it != bc_m.end(); ++it) { doesChange |= (*it)->changesPhysicalCells(); diff --git a/src/Field/BcTypes.h b/src/Field/BcTypes.h index 16d083360..812c87fff 100644 --- a/src/Field/BcTypes.h +++ b/src/Field/BcTypes.h @@ -37,7 +37,7 @@ #include "Types/ViewTypes.h" namespace ippl { - template + template class Field; /* @@ -56,16 +56,16 @@ namespace ippl { }; namespace detail { - template + template class BCondBase; - template - std::ostream& operator<<(std::ostream&, const BCondBase&); + template + std::ostream& operator<<(std::ostream&, const BCondBase&); - template + template class BCondBase { public: - using Field_t = Field; + using Field_t = Field; using Layout_t = FieldLayout; // Constructor takes: @@ -78,8 +78,8 @@ namespace ippl { virtual FieldBC getBCType() const { return NO_FACE; } - virtual void findBCNeighbors(Field& field) = 0; - virtual void apply(Field& field) = 0; + virtual void findBCNeighbors(Field& field) = 0; + virtual void apply(Field& field) = 0; virtual void write(std::ostream&) const = 0; // Return face on which BC applies @@ -98,16 +98,16 @@ namespace ippl { } // namespace detail - template - class ExtrapolateFace : public detail::BCondBase { + template + class ExtrapolateFace : public detail::BCondBase { public: // Constructor takes zero, one, or two int's specifying components of // multicomponent types like Vector this BC applies to. // Zero int's specified means apply to all components; one means apply to // component (i), and two means apply to component (i,j), - using base_type = detail::BCondBase; - using Field_t = typename detail::BCondBase::Field_t; - using Layout_t = typename detail::BCondBase::Layout_t; + using base_type = detail::BCondBase; + using Field_t = typename detail::BCondBase::Field_t; + using Layout_t = typename detail::BCondBase::Layout_t; ExtrapolateFace(unsigned face, T offset, T slope) : base_type(face) @@ -131,12 +131,12 @@ namespace ippl { T slope_m; }; - template - class NoBcFace : public detail::BCondBase { + template + class NoBcFace : public detail::BCondBase { public: - using Field_t = typename detail::BCondBase::Field_t; + using Field_t = typename detail::BCondBase::Field_t; NoBcFace(int face) - : detail::BCondBase(face) {} + : detail::BCondBase(face) {} virtual void findBCNeighbors(Field_t& /*field*/) {} virtual void apply(Field_t& /*field*/) {} @@ -144,37 +144,37 @@ namespace ippl { virtual void write(std::ostream& out) const; }; - template - class ConstantFace : public ExtrapolateFace { + template + class ConstantFace : public ExtrapolateFace { public: ConstantFace(unsigned int face, T constant) - : ExtrapolateFace(face, constant, 0) {} + : ExtrapolateFace(face, constant, 0) {} virtual FieldBC getBCType() const { return CONSTANT_FACE; } virtual void write(std::ostream& out) const; }; - template - class ZeroFace : public ConstantFace { + template + class ZeroFace : public ConstantFace { public: ZeroFace(unsigned face) - : ConstantFace(face, 0.0) {} + : ConstantFace(face, 0.0) {} virtual FieldBC getBCType() const { return ZERO_FACE; } virtual void write(std::ostream& out) const; }; - template - class PeriodicFace : public detail::BCondBase { + template + class PeriodicFace : public detail::BCondBase { public: using face_neighbor_type = std::array, 2 * Dim>; - using Field_t = typename detail::BCondBase::Field_t; - using Layout_t = typename detail::BCondBase::Layout_t; + using Field_t = typename detail::BCondBase::Field_t; + using Layout_t = typename detail::BCondBase::Layout_t; PeriodicFace(unsigned face) - : detail::BCondBase(face) {} + : detail::BCondBase(face) {} virtual FieldBC getBCType() const { return PERIODIC_FACE; } diff --git a/src/Field/BcTypes.hpp b/src/Field/BcTypes.hpp index 56a142c1d..ef3f49077 100644 --- a/src/Field/BcTypes.hpp +++ b/src/Field/BcTypes.hpp @@ -33,21 +33,21 @@ namespace ippl { namespace detail { - template - BCondBase::BCondBase(unsigned int face) + template + BCondBase::BCondBase(unsigned int face) : face_m(face) , changePhysical_m(false) {} - template - inline std::ostream& operator<<(std::ostream& os, const BCondBase& bc) { + template + inline std::ostream& operator<<(std::ostream& os, const BCondBase& bc) { bc.write(os); return os; } } // namespace detail - template - void ExtrapolateFace::apply(Field_t& field) { + template + void ExtrapolateFace::apply(Field_t& field) { // We only support constant extrapolation for the moment, other // higher order extrapolation stuffs need to be added. @@ -123,38 +123,38 @@ namespace ippl { } } - template - void ExtrapolateFace::write(std::ostream& out) const { + template + void ExtrapolateFace::write(std::ostream& out) const { out << "Constant Extrapolation Face" << ", Face = " << this->face_m; } - template - void NoBcFace::write(std::ostream& out) const { + template + void NoBcFace::write(std::ostream& out) const { out << "NoBcFace" << ", Face = " << this->face_m; } - template - void ConstantFace::write(std::ostream& out) const { + template + void ConstantFace::write(std::ostream& out) const { out << "ConstantFace" << ", Face = " << this->face_m << ", Constant = " << this->offset_m; } - template - void ZeroFace::write(std::ostream& out) const { + template + void ZeroFace::write(std::ostream& out) const { out << "ZeroFace" << ", Face = " << this->face_m; } - template - void PeriodicFace::write(std::ostream& out) const { + template + void PeriodicFace::write(std::ostream& out) const { out << "PeriodicFace" << ", Face = " << this->face_m; } - template - void PeriodicFace::findBCNeighbors(Field_t& field) { + template + void PeriodicFace::findBCNeighbors(Field_t& field) { // For cell centering only face neighbors are needed unsigned int face = this->face_m; unsigned int d = face / 2; @@ -207,8 +207,8 @@ namespace ippl { } } - template - void PeriodicFace::apply(Field_t& field) { + template + void PeriodicFace::apply(Field_t& field) { unsigned int face = this->face_m; unsigned int d = face / 2; typename Field_t::view_type& view = field.getView(); diff --git a/src/Field/Field.h b/src/Field/Field.h index 6ec90bf21..61f0555f5 100644 --- a/src/Field/Field.h +++ b/src/Field/Field.h @@ -25,7 +25,7 @@ namespace ippl { - template + template class Field : public BareField { public: typedef T type; @@ -36,7 +36,7 @@ namespace ippl { using Layout_t = FieldLayout; using BareField_t = BareField; using view_type = typename BareField_t::view_type; - using BConds_t = BConds; + using BConds_t = BConds; // A default constructor, which should be used only if the user calls the // 'initialize' function before doing anything else. There are no special diff --git a/src/Field/Field.hpp b/src/Field/Field.hpp index dd03c9e61..f61f9f869 100644 --- a/src/Field/Field.hpp +++ b/src/Field/Field.hpp @@ -19,8 +19,8 @@ namespace ippl { namespace detail { - template - struct isExpression> : std::true_type {}; + template + struct isExpression> : std::true_type {}; } // namespace detail ////////////////////////////////////////////////////////////////////////// @@ -28,47 +28,47 @@ namespace ippl { // 'initialize' function before doing anything else. There are no special // checks in the rest of the Field methods to check that the Field has // been properly initialized - template - Field::Field() + template + Field::Field() : BareField() , mesh_m(nullptr) , bc_m() {} ////////////////////////////////////////////////////////////////////////// // Constructors which include a Mesh object as argument - template - Field::Field(Mesh_t& m, Layout_t& l, int nghost) + template + Field::Field(Mesh_t& m, Layout_t& l, int nghost) : BareField(l, nghost) , mesh_m(&m) { for (unsigned int face = 0; face < 2 * Dim; ++face) { - bc_m[face] = std::make_shared >(face); + bc_m[face] = std::make_shared >(face); } } ////////////////////////////////////////////////////////////////////////// // Initialize the Field, also specifying a mesh - template - void Field::initialize(Mesh_t& m, Layout_t& l, int nghost) { + template + void Field::initialize(Mesh_t& m, Layout_t& l, int nghost) { BareField::initialize(l, nghost); mesh_m = &m; for (unsigned int face = 0; face < 2 * Dim; ++face) { - bc_m[face] = std::make_shared >(face); + bc_m[face] = std::make_shared >(face); } } - template - T Field::getVolumeIntegral() const { + template + T Field::getVolumeIntegral() const { typename Mesh::value_type dV = mesh_m->getCellVolume(); return this->sum() * dV; } - template - T Field::getVolumeAverage() const { + template + T Field::getVolumeAverage() const { return getVolumeIntegral() / mesh_m->getMeshVolume(); } - template - void Field::updateLayout(Layout_t& l, int nghost) { + template + void Field::updateLayout(Layout_t& l, int nghost) { BareField::updateLayout(l, nghost); } diff --git a/src/Field/FieldOperations.hpp b/src/Field/FieldOperations.hpp index aa168fbb4..a8e23bc37 100644 --- a/src/Field/FieldOperations.hpp +++ b/src/Field/FieldOperations.hpp @@ -23,8 +23,8 @@ namespace ippl { * @param f2 second field * @return Result of f1^T f2 */ - template - T innerProduct(const Field& f1, const Field& f2) { + template + T innerProduct(const Field& f1, const Field& f2) { T sum = 0; auto view1 = f1.getView(); auto view2 = f2.getView(); @@ -46,8 +46,8 @@ namespace ippl { * @param p desired norm (default 2) * @return The desired norm of the field */ - template - T norm(const Field& field, int p = 2) { + template + T norm(const Field& field, int p = 2) { T local = 0; auto view = field.getView(); switch (p) { @@ -86,10 +86,10 @@ namespace ippl { * User interface of gradient in three dimensions. * @param u field */ - template - detail::meta_grad > grad(Field& u) { + template + detail::meta_grad > grad(Field& u) { u.fillHalo(); - BConds& bcField = u.getFieldBC(); + BConds& bcField = u.getFieldBC(); bcField.apply(u); Mesh& mesh = u.get_mesh(); typename Mesh::vector_type xvector(0); @@ -98,17 +98,17 @@ namespace ippl { yvector[1] = 0.5 / mesh.getMeshSpacing(1); typename Mesh::vector_type zvector(0); zvector[2] = 0.5 / mesh.getMeshSpacing(2); - return detail::meta_grad >(u, xvector, yvector, zvector); + return detail::meta_grad >(u, xvector, yvector, zvector); } /*! * User interface of divergence in three dimensions. * @param u field */ - template - detail::meta_div > div(Field& u) { + template + detail::meta_div > div(Field& u) { u.fillHalo(); - BConds& bcField = u.getFieldBC(); + BConds& bcField = u.getFieldBC(); bcField.apply(u); Mesh& mesh = u.get_mesh(); typename Mesh::vector_type xvector(0); @@ -117,34 +117,34 @@ namespace ippl { yvector[1] = 0.5 / mesh.getMeshSpacing(1); typename Mesh::vector_type zvector(0); zvector[2] = 0.5 / mesh.getMeshSpacing(2); - return detail::meta_div >(u, xvector, yvector, zvector); + return detail::meta_div >(u, xvector, yvector, zvector); } /*! * User interface of Laplacian in three dimensions. * @param u field */ - template - detail::meta_laplace > laplace(Field& u) { + template + detail::meta_laplace > laplace(Field& u) { u.fillHalo(); - BConds& bcField = u.getFieldBC(); + BConds& bcField = u.getFieldBC(); bcField.apply(u); Mesh& mesh = u.get_mesh(); typename Mesh::vector_type hvector(0); hvector[0] = 1.0 / std::pow(mesh.getMeshSpacing(0), 2); hvector[1] = 1.0 / std::pow(mesh.getMeshSpacing(1), 2); hvector[2] = 1.0 / std::pow(mesh.getMeshSpacing(2), 2); - return detail::meta_laplace >(u, hvector); + return detail::meta_laplace >(u, hvector); } /*! * User interface of curl in three dimensions. * @param u field */ - template - detail::meta_curl > curl(Field& u) { + template + detail::meta_curl > curl(Field& u) { u.fillHalo(); - BConds& bcField = u.getFieldBC(); + BConds& bcField = u.getFieldBC(); bcField.apply(u); Mesh& mesh = u.get_mesh(); typename Mesh::vector_type xvector(0); @@ -155,17 +155,17 @@ namespace ippl { zvector[2] = 1.0; typename Mesh::vector_type hvector(0); hvector = mesh.getMeshSpacing(); - return detail::meta_curl >(u, xvector, yvector, zvector, hvector); + return detail::meta_curl >(u, xvector, yvector, zvector, hvector); } /*! * User interface of Hessian in three dimensions. * @param u field */ - template - detail::meta_hess > hess(Field& u) { + template + detail::meta_hess > hess(Field& u) { u.fillHalo(); - BConds& bcField = u.getFieldBC(); + BConds& bcField = u.getFieldBC(); bcField.apply(u); Mesh& mesh = u.get_mesh(); typename Mesh::vector_type xvector(0); @@ -176,6 +176,6 @@ namespace ippl { zvector[2] = 1.0; typename Mesh::vector_type hvector(0); hvector = mesh.getMeshSpacing(); - return detail::meta_hess >(u, xvector, yvector, zvector, hvector); + return detail::meta_hess >(u, xvector, yvector, zvector, hvector); } } // namespace ippl diff --git a/src/Particle/IntNGP.h b/src/Particle/IntNGP.h index 8db494b3c..65cb77b14 100644 --- a/src/Particle/IntNGP.h +++ b/src/Particle/IntNGP.h @@ -38,8 +38,8 @@ class IntNGP : public Interpolator { // gather/scatter functions // scatter particle data into Field using particle position and mesh - template - static void scatter(const FT& pdata, Field& f, const Vektor& ppos, + template + static void scatter(const FT& pdata, Field& f, const Vektor& ppos, const Mesh& mesh) { // find nearest-grid-point for particle position, store in NDIndex obj NDIndex ngp = FindNGP(mesh, ppos, CenteringTag()); @@ -54,8 +54,8 @@ class IntNGP : public Interpolator { // scatter particle data into Field using particle position and mesh // and cache mesh information for reuse - template - static void scatter(const FT& pdata, Field& f, const Vektor& ppos, + template + static void scatter(const FT& pdata, Field& f, const Vektor& ppos, const Mesh& mesh, NDIndex& ngp) { // find nearest-grid-point for particle position, store in NDIndex obj ngp = FindNGP(mesh, ppos, CenteringTag()); @@ -69,8 +69,8 @@ class IntNGP : public Interpolator { } // scatter particle data into Field using cached mesh information - template - static void scatter(const FT& pdata, Field& f, const NDIndex& ngp) { + template + static void scatter(const FT& pdata, Field& f, const NDIndex& ngp) { // scatter data value to Field ... this assumes that the Field // data point is local to this processor, if not an error will be printed. @@ -81,8 +81,8 @@ class IntNGP : public Interpolator { } // gather particle data from Field using particle position and mesh - template - static void gather(FT& pdata, const Field& f, const Vektor& ppos, + template + static void gather(FT& pdata, const Field& f, const Vektor& ppos, const Mesh& mesh) { // find nearest-grid-point for particle position, store in NDIndex obj NDIndex ngp = FindNGP(mesh, ppos, CenteringTag()); @@ -97,8 +97,8 @@ class IntNGP : public Interpolator { // gather particle data from Field using particle position and mesh // and cache mesh information for reuse - template - static void gather(FT& pdata, const Field& f, const Vektor& ppos, + template + static void gather(FT& pdata, const Field& f, const Vektor& ppos, const Mesh& mesh, NDIndex& ngp) { // find nearest-grid-point for particle position, store in NDIndex obj ngp = FindNGP(mesh, ppos, CenteringTag()); @@ -112,8 +112,8 @@ class IntNGP : public Interpolator { } // gather particle data from Field using cached mesh information - template - static void gather(FT& pdata, const Field& f, const NDIndex& ngp) { + template + static void gather(FT& pdata, const Field& f, const NDIndex& ngp) { // gather Field value to particle data ... this assumes that the Field // data point is local to this processor, if not an error will be printed. diff --git a/src/Particle/ParticleBalancer.hpp b/src/Particle/ParticleBalancer.hpp index a13378ae3..540733ce1 100644 --- a/src/Particle/ParticleBalancer.hpp +++ b/src/Particle/ParticleBalancer.hpp @@ -78,7 +78,7 @@ bool BinaryRepartition(IpplParticleBase BF(mesh, FL, GuardCellSizes(1)); + Field BF(mesh, FL, GuardCellSizes(1)); // Now do a number density scatter on this Field // Afterwards, the Field will be deleted, and will checkout of the diff --git a/src/Solver/Electrostatics.h b/src/Solver/Electrostatics.h index 2ed346447..65b9465a9 100644 --- a/src/Solver/Electrostatics.h +++ b/src/Solver/Electrostatics.h @@ -23,12 +23,12 @@ namespace ippl { - template - class Electrostatics : public Solver { + template + class Electrostatics : public Solver { public: - using grad_type = Field, Dim, Mesh, Cell>; - using lhs_type = typename Solver::lhs_type; - using rhs_type = typename Solver::rhs_type; + using grad_type = Field, Dim, Mesh, Centering>; + using lhs_type = typename Solver::lhs_type; + using rhs_type = typename Solver::rhs_type; /*! * Represents the types of fields that should @@ -45,13 +45,13 @@ namespace ippl { * desired output type defaults to solution only */ Electrostatics() - : Solver() + : Solver() , grad_mp(nullptr) { setDefaultParameters(); } Electrostatics(lhs_type& lhs, rhs_type& rhs) - : Solver(lhs, rhs) + : Solver(lhs, rhs) , grad_mp(nullptr) { setDefaultParameters(); } diff --git a/src/Solver/ElectrostaticsCG.h b/src/Solver/ElectrostaticsCG.h index 3c6bcaaf2..797ebddf4 100644 --- a/src/Solver/ElectrostaticsCG.h +++ b/src/Solver/ElectrostaticsCG.h @@ -32,14 +32,14 @@ namespace ippl { return fun(arg); \ } - template - class ElectrostaticsCG : public Electrostatics { + template + class ElectrostaticsCG : public Electrostatics { public: - using lhs_type = typename Solver::lhs_type; - using rhs_type = typename Solver::rhs_type; + using lhs_type = typename Solver::lhs_type; + using rhs_type = typename Solver::rhs_type; using OpRet = UnaryMinus>; - using algo = PCG; - using Base = Electrostatics; + using algo = PCG; + using Base = Electrostatics; ElectrostaticsCG() : Base() { diff --git a/src/Solver/FFTPeriodicPoissonSolver.h b/src/Solver/FFTPeriodicPoissonSolver.h index ecd05a844..48a2718b0 100644 --- a/src/Solver/FFTPeriodicPoissonSolver.h +++ b/src/Solver/FFTPeriodicPoissonSolver.h @@ -28,8 +28,8 @@ namespace ippl { - template - class FFTPeriodicPoissonSolver : public Electrostatics { + template + class FFTPeriodicPoissonSolver : public Electrostatics { public: using Field_t = Field; using FFT_t = FFT; @@ -38,9 +38,9 @@ namespace ippl { using Layout_t = FieldLayout; using Vector_t = Vector; - using Base = Electrostatics; - using lhs_type = typename Solver::lhs_type; - using rhs_type = typename Solver::rhs_type; + using Base = Electrostatics; + using lhs_type = typename Solver::lhs_type; + using rhs_type = typename Solver::rhs_type; FFTPeriodicPoissonSolver() : Base() { diff --git a/src/Solver/FFTPeriodicPoissonSolver.hpp b/src/Solver/FFTPeriodicPoissonSolver.hpp index 1748632ff..99b0eb62c 100644 --- a/src/Solver/FFTPeriodicPoissonSolver.hpp +++ b/src/Solver/FFTPeriodicPoissonSolver.hpp @@ -19,14 +19,14 @@ namespace ippl { - template - void FFTPeriodicPoissonSolver::setRhs(rhs_type& rhs) { + template + void FFTPeriodicPoissonSolver::setRhs(rhs_type& rhs) { Base::setRhs(rhs); initialize(); } - template - void FFTPeriodicPoissonSolver::initialize() { + template + void FFTPeriodicPoissonSolver::initialize() { const Layout_t& layout_r = this->rhs_mp->getLayout(); domain_m = layout_r.getDomain(); @@ -56,8 +56,8 @@ namespace ippl { fft_mp = std::make_shared(layout_r, *layoutComplex_mp, this->params_m); } - template - void FFTPeriodicPoissonSolver::solve() { + template + void FFTPeriodicPoissonSolver::solve() { fft_mp->transform(1, *this->rhs_mp, fieldComplex_m); auto view = fieldComplex_m.getView(); diff --git a/src/Solver/FFTPoissonSolver.h b/src/Solver/FFTPoissonSolver.h index ba48eab75..e0105397c 100644 --- a/src/Solver/FFTPoissonSolver.h +++ b/src/Solver/FFTPoissonSolver.h @@ -25,22 +25,22 @@ #include "Types/Vector.h" namespace ippl { - template - class FFTPoissonSolver : public Electrostatics { + template + class FFTPoissonSolver : public Electrostatics { public: // types for LHS and RHS - using lhs_type = typename Solver::lhs_type; - using rhs_type = typename Solver::rhs_type; + using lhs_type = typename Solver::lhs_type; + using rhs_type = typename Solver::rhs_type; // type of output - using Base = Electrostatics; + using Base = Electrostatics; // define a type for a 3 dimensional field (e.g. charge density field) // define a type of Field with integers to be used for the helper Green's function // also define a type for the Fourier transformed complex valued fields - typedef Field Field_t; - typedef Field IField_t; - typedef Field, Dim, Mesh, Cell> CxField_t; + typedef Field Field_t; + typedef Field IField_t; + typedef Field, Dim, Mesh, Centering> CxField_t; typedef Vector Vector_t; // define type for field layout diff --git a/src/Solver/FFTPoissonSolver.hpp b/src/Solver/FFTPoissonSolver.hpp index b23e6e790..57b88c6f3 100644 --- a/src/Solver/FFTPoissonSolver.hpp +++ b/src/Solver/FFTPoissonSolver.hpp @@ -130,8 +130,8 @@ namespace ippl { ///////////////////////////////////////////////////////////////////////// // constructor and destructor - template - FFTPoissonSolver::FFTPoissonSolver(rhs_type& rhs, + template + FFTPoissonSolver::FFTPoissonSolver(rhs_type& rhs, ParameterList& fftparams, std::string alg) : mesh_mp(nullptr) @@ -160,8 +160,8 @@ namespace ippl { IpplTimings::stopTimer(initialize); } - template - FFTPoissonSolver::FFTPoissonSolver(lhs_type& lhs, rhs_type& rhs, + template + FFTPoissonSolver::FFTPoissonSolver(lhs_type& lhs, rhs_type& rhs, ParameterList& fftparams, std::string alg, int sol) : mesh_mp(nullptr) @@ -191,15 +191,15 @@ namespace ippl { IpplTimings::stopTimer(initialize); } - template - FFTPoissonSolver::~FFTPoissonSolver(){}; + template + FFTPoissonSolver::~FFTPoissonSolver(){}; ///////////////////////////////////////////////////////////////////////// // allows user to set gradient of phi = Efield instead of spectral // calculation of Efield (which uses FFTs) - template - void FFTPoissonSolver::setGradFD() { + template + void FFTPoissonSolver::setGradFD() { // get the output type (sol, grad, or sol & grad) const int out = this->params_m.template get("output_type"); @@ -215,8 +215,8 @@ namespace ippl { ///////////////////////////////////////////////////////////////////////// // initializeFields method, called in constructor - template - void FFTPoissonSolver::initializeFields() { + template + void FFTPoissonSolver::initializeFields() { // first check if valid algorithm choice if ((alg_m != "VICO") && (alg_m != "HOCKNEY") && (alg_m != "BIHARMONIC")) { throw IpplException( @@ -424,8 +424,8 @@ namespace ippl { ///////////////////////////////////////////////////////////////////////// // compute electric potential by solving Poisson's eq given a field rho and mesh spacings hr - template - void FFTPoissonSolver::solve() { + template + void FFTPoissonSolver::solve() { // start a timer static IpplTimings::TimerRef solve = IpplTimings::getTimer("Solve"); IpplTimings::startTimer(solve); @@ -841,8 +841,8 @@ namespace ippl { //////////////////////////////////////////////////////////////////////// // calculate FFT of the Green's function - template - void FFTPoissonSolver::greensFunction() { + template + void FFTPoissonSolver::greensFunction() { const double pi = std::acos(-1.0); grn_mr = 0.0; @@ -1061,8 +1061,8 @@ namespace ippl { IpplTimings::stopTimer(fftg); }; - template - void FFTPoissonSolver::communicateVico( + template + void FFTPoissonSolver::communicateVico( Vector size, typename CxField_t::view_type view_g, const ippl::NDIndex ldom_g, const int nghost_g, typename Field_t::view_type view, const ippl::NDIndex ldom, const int nghost) { diff --git a/src/Solver/PCG.h b/src/Solver/PCG.h index 36a39cc67..d592f58bd 100644 --- a/src/Solver/PCG.h +++ b/src/Solver/PCG.h @@ -24,10 +24,10 @@ namespace ippl { template + class Mesh, class Centering> class PCG : public SolverAlgorithm { public: - using Base = SolverAlgorithm; + using Base = SolverAlgorithm; using typename Base::lhs_type; using typename Base::rhs_type; using operator_type = std::function; diff --git a/src/Solver/Solver.h b/src/Solver/Solver.h index a48800c6e..14db2c07f 100644 --- a/src/Solver/Solver.h +++ b/src/Solver/Solver.h @@ -24,11 +24,11 @@ namespace ippl { - template + template class Solver { public: - using lhs_type = Field; - using rhs_type = Field; + using lhs_type = Field; + using rhs_type = Field; /*! * Default constructor diff --git a/src/Solver/SolverAlgorithm.h b/src/Solver/SolverAlgorithm.h index 47dfab46b..79afa5575 100644 --- a/src/Solver/SolverAlgorithm.h +++ b/src/Solver/SolverAlgorithm.h @@ -24,11 +24,11 @@ namespace ippl { - template + template class SolverAlgorithm { public: - using lhs_type = Field; - using rhs_type = Field; + using lhs_type = Field; + using rhs_type = Field; /*! * Solve the problem described by Op(lhs) = rhs, where Op is an unspecified From 66275e8cd29ba2c0a9f39b2d277965a4a4256be5 Mon Sep 17 00:00:00 2001 From: Matthias Frey Date: Fri, 17 Mar 2023 11:33:09 +0000 Subject: [PATCH 06/12] remove include statement --- src/Expression/IpplExpressions.h | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/Expression/IpplExpressions.h b/src/Expression/IpplExpressions.h index a782c45b1..2986531d0 100644 --- a/src/Expression/IpplExpressions.h +++ b/src/Expression/IpplExpressions.h @@ -20,8 +20,6 @@ #include -#include - namespace ippl { namespace detail { /*! From c4f07519a27162bb4b68e66a4aad9f57455f155a Mon Sep 17 00:00:00 2001 From: Matthias Frey Date: Sun, 19 Mar 2023 11:20:26 +0000 Subject: [PATCH 07/12] make tests compile --- alpine/ChargedParticles.hpp | 7 +- src/Solver/FFTPeriodicPoissonSolver.h | 6 +- src/Solver/FFTPeriodicPoissonSolver.hpp | 6 +- src/Solver/FFTPoissonSolver.h | 6 +- src/Solver/FFTPoissonSolver.hpp | 68 ++++++++++--------- src/Solver/PCG.h | 8 +-- src/Solver/test/Budiardja_plot.cpp | 12 ++-- src/Solver/test/TestCGSolver.cpp | 12 ++-- .../test/TestFFTPeriodicPoissonSolver.cpp | 10 +-- src/Solver/test/TestGaussian.cpp | 15 ++-- src/Solver/test/TestGaussian_biharmonic.cpp | 31 +++++---- src/Solver/test/TestGaussian_convergence.cpp | 34 ++++++---- src/Solver/test/TestSolverDesign.cpp | 10 +-- src/Solver/test/TestSphere.cpp | 18 +++-- test/FFT/TestCos.cpp | 8 ++- test/FFT/TestFFTCC.cpp | 8 ++- test/FFT/TestFFTRC.cpp | 10 +-- test/FFT/TestSine.cpp | 8 ++- test/field/TestCurl.cpp | 6 +- test/field/TestFieldBC.cpp | 15 ++-- test/field/TestHalo.cpp | 6 +- test/field/TestHessian.cpp | 8 ++- test/field/TestLaplace.cpp | 23 ++++--- test/particle/PIC3d.cpp | 5 +- test/particle/TestGather.cpp | 6 +- test/particle/TestScatter.cpp | 6 +- 26 files changed, 203 insertions(+), 149 deletions(-) diff --git a/alpine/ChargedParticles.hpp b/alpine/ChargedParticles.hpp index 8f0e204e6..4ac968a18 100644 --- a/alpine/ChargedParticles.hpp +++ b/alpine/ChargedParticles.hpp @@ -26,8 +26,9 @@ constexpr unsigned Dim = 3; // some typedefs typedef ippl::ParticleSpatialLayout PLayout_t; typedef ippl::UniformCartesian Mesh_t; +typedef Mesh_t::DefaultCentering Centering_t; typedef ippl::FieldLayout FieldLayout_t; -typedef ippl::OrthogonalRecursiveBisection ORB; +typedef ippl::OrthogonalRecursiveBisection ORB; using size_type = ippl::detail::size_type; @@ -35,7 +36,7 @@ template using Vector = ippl::Vector; template -using Field = ippl::Field; +using Field = ippl::Field; template using ParticleAttrib = ippl::ParticleAttrib; @@ -43,7 +44,7 @@ using ParticleAttrib = ippl::ParticleAttrib; typedef Vector Vector_t; typedef Field Field_t; typedef Field VField_t; -typedef ippl::FFTPeriodicPoissonSolver Solver_t; +typedef ippl::FFTPeriodicPoissonSolver Solver_t; const double pi = std::acos(-1.0); diff --git a/src/Solver/FFTPeriodicPoissonSolver.h b/src/Solver/FFTPeriodicPoissonSolver.h index 26b0ff92e..c46cd42ee 100644 --- a/src/Solver/FFTPeriodicPoissonSolver.h +++ b/src/Solver/FFTPeriodicPoissonSolver.h @@ -32,10 +32,10 @@ namespace ippl { template class FFTPeriodicPoissonSolver : public Electrostatics { public: - using Field_t = Field; - using FFT_t = FFT; + using Field_t = Field; + using FFT_t = FFT; using Complex_t = Kokkos::complex; - using CxField_t = Field; + using CxField_t = Field; using Layout_t = FieldLayout; using Vector_t = Vector; diff --git a/src/Solver/FFTPeriodicPoissonSolver.hpp b/src/Solver/FFTPeriodicPoissonSolver.hpp index 99b0eb62c..6a4c27b21 100644 --- a/src/Solver/FFTPeriodicPoissonSolver.hpp +++ b/src/Solver/FFTPeriodicPoissonSolver.hpp @@ -46,7 +46,7 @@ namespace ippl { Vector hComplex = {1.0, 1.0, 1.0}; Vector originComplex = {0.0, 0.0, 0.0}; - M meshComplex(domainComplex, hComplex, originComplex); + Mesh meshComplex(domainComplex, hComplex, originComplex); fieldComplex_m.initialize(meshComplex, *layoutComplex_mp); @@ -65,9 +65,9 @@ namespace ippl { using mdrange_type = Kokkos::MDRangePolicy>; double pi = std::acos(-1.0); - const M& mesh = this->rhs_mp->get_mesh(); + const Mesh& mesh = this->rhs_mp->get_mesh(); const auto& lDomComplex = layoutComplex_mp->getLocalNDIndex(); - using vector_type = typename M::vector_type; + using vector_type = typename Mesh::vector_type; const vector_type& origin = mesh.getOrigin(); const vector_type& hx = mesh.getMeshSpacing(); diff --git a/src/Solver/FFTPoissonSolver.h b/src/Solver/FFTPoissonSolver.h index 741b3cf82..e18bcc5b2 100644 --- a/src/Solver/FFTPoissonSolver.h +++ b/src/Solver/FFTPoissonSolver.h @@ -49,7 +49,7 @@ namespace ippl { typedef FieldLayout FieldLayout_t; // define a type for the 3 dimensional real to complex Fourier transform - typedef FFT FFT_t; + typedef FFT FFT_t; // type for communication buffers using buffer_type = Communicate::buffer_type; @@ -107,7 +107,7 @@ namespace ippl { std::unique_ptr fft_m; // mesh and layout objects for rho_m (RHS) - M* mesh_mp; + Mesh* mesh_mp; FieldLayout_t* layout_mp; // mesh and layout objects for rho2_m @@ -133,7 +133,7 @@ namespace ippl { // members for Vico-Greengard CxField_t grnL_m; - std::unique_ptr> fft4n_m; + std::unique_ptr> fft4n_m; std::unique_ptr mesh4_m; std::unique_ptr layout4_m; diff --git a/src/Solver/FFTPoissonSolver.hpp b/src/Solver/FFTPoissonSolver.hpp index 3118b4f9c..54f9fbc6b 100644 --- a/src/Solver/FFTPoissonSolver.hpp +++ b/src/Solver/FFTPoissonSolver.hpp @@ -24,11 +24,11 @@ #include "Field/HaloCells.h" // Communication specific functions (pack and unpack). -template +template void pack(const ippl::NDIndex<3> intersect, Kokkos::View& view, ippl::detail::FieldBufferData& fd, int nghost, const ippl::NDIndex<3> ldom, ippl::Communicate::size_type& nsends) { - ippl::Field::view_type& buffer = fd.buffer; + typename ippl::Field::view_type& buffer = fd.buffer; size_t size = intersect.size(); nsends = size; @@ -66,10 +66,12 @@ void pack(const ippl::NDIndex<3> intersect, Kokkos::View& view, Kokkos::fence(); } -void unpack(const ippl::NDIndex<3> intersect, const ippl::Field::view_type& view, +template +void unpack(const ippl::NDIndex<3> intersect, + const typename ippl::Field::view_type& view, ippl::detail::FieldBufferData& fd, int nghost, const ippl::NDIndex<3> ldom, bool x = false, bool y = false, bool z = false) { - ippl::Field::view_type& buffer = fd.buffer; + typename ippl::Field::view_type& buffer = fd.buffer; const int first0 = intersect[0].first() + nghost - ldom[0].first(); const int first1 = intersect[1].first() + nghost - ldom[1].first(); @@ -99,10 +101,12 @@ void unpack(const ippl::NDIndex<3> intersect, const ippl::Field::view Kokkos::fence(); } +template void unpack(const ippl::NDIndex<3> intersect, - const ippl::Field, 3>::view_type& view, size_t dim, + const typename ippl::Field, 3, Mesh, Centering>::view_type& view, + size_t dim, ippl::detail::FieldBufferData& fd, int nghost, const ippl::NDIndex<3> ldom) { - ippl::Field::view_type& buffer = fd.buffer; + typename ippl::Field::view_type& buffer = fd.buffer; const int first0 = intersect[0].first() + nghost - ldom[0].first(); const int first1 = intersect[1].first() + nghost - ldom[1].first(); @@ -260,7 +264,7 @@ namespace ippl { } // create double sized mesh and layout objects using the previously defined domain2_m - mesh2_m = std::unique_ptr(new M(domain2_m, hr_m, origin)); + mesh2_m = std::unique_ptr(new Mesh(domain2_m, hr_m, origin)); layout2_m = std::unique_ptr(new FieldLayout_t(domain2_m, decomp)); // create the domain for the transformed (complex) fields @@ -276,7 +280,7 @@ namespace ippl { } // create mesh and layout for the real to complex FFT transformed fields - meshComplex_m = std::unique_ptr(new M(domainComplex_m, hr_m, origin)); + meshComplex_m = std::unique_ptr(new Mesh(domainComplex_m, hr_m, origin)); layoutComplex_m = std::unique_ptr(new FieldLayout_t(domainComplex_m, decomp)); @@ -307,14 +311,14 @@ namespace ippl { } // 4N grid - mesh4_m = std::unique_ptr(new M(domain4_m, hr_m, origin)); + mesh4_m = std::unique_ptr(new Mesh(domain4_m, hr_m, origin)); layout4_m = std::unique_ptr(new FieldLayout_t(domain4_m, decomp)); // initialize fields grnL_m.initialize(*mesh4_m, *layout4_m); // create a Complex-to-Complex FFT object to transform for layout4 - fft4n_m = std::make_unique>(*layout4_m, this->params_m); + fft4n_m = std::make_unique>(*layout4_m, this->params_m); IpplTimings::stopTimer(initialize_vico); } @@ -489,7 +493,7 @@ namespace ippl { requests.resize(requests.size() + 1); Communicate::size_type nsends; - pack(intersection, view1, fd_m, nghost1, ldom1, nsends); + pack(intersection, view1, fd_m, nghost1, ldom1, nsends); buffer_type buf = Ippl::Comm->getBuffer(IPPL_SOLVER_SEND + i, nsends); @@ -516,7 +520,7 @@ namespace ippl { nrecvs); buf->resetReadPos(); - unpack(intersection, view2, fd_m, nghost2, ldom2); + unpack(intersection, view2, fd_m, nghost2, ldom2); } } @@ -612,7 +616,7 @@ namespace ippl { requests.resize(requests.size() + 1); Communicate::size_type nsends; - pack(intersection, view2, fd_m, nghost2, ldom2, nsends); + pack(intersection, view2, fd_m, nghost2, ldom2, nsends); buffer_type buf = Ippl::Comm->getBuffer(IPPL_SOLVER_SEND + i, nsends); @@ -640,7 +644,7 @@ namespace ippl { nrecvs); buf->resetReadPos(); - unpack(intersection, view1, fd_m, nghost1, ldom1); + unpack(intersection, view1, fd_m, nghost1, ldom1); } } @@ -774,7 +778,7 @@ namespace ippl { requests.resize(requests.size() + 1); Communicate::size_type nsends; - pack(intersection, view2, fd_m, nghost2, ldom2, nsends); + pack(intersection, view2, fd_m, nghost2, ldom2, nsends); buffer_type buf = Ippl::Comm->getBuffer(IPPL_SOLVER_SEND + i, nsends); @@ -803,7 +807,7 @@ namespace ippl { nrecvs * sizeof(double), nrecvs); buf->resetReadPos(); - unpack(intersection, viewL, gd, fd_m, nghostL, ldom1); + unpack(intersection, viewL, gd, fd_m, nghostL, ldom1); } } @@ -1128,7 +1132,7 @@ namespace ippl { requests.resize(requests.size() + 1); Communicate::size_type nsends; - pack(intersection, view_g, fd_m, nghost_g, ldom_g, nsends); + pack(intersection, view_g, fd_m, nghost_g, ldom_g, nsends); buffer_type buf = Ippl::Comm->getBuffer(IPPL_VICO_SEND + i, nsends); @@ -1155,7 +1159,7 @@ namespace ippl { requests.resize(requests.size() + 1); Communicate::size_type nsends; - pack(intersection, view_g, fd_m, nghost_g, ldom_g, nsends); + pack(intersection, view_g, fd_m, nghost_g, ldom_g, nsends); buffer_type buf = Ippl::Comm->getBuffer(IPPL_VICO_SEND + 8 + i, nsends); @@ -1182,7 +1186,7 @@ namespace ippl { requests.resize(requests.size() + 1); Communicate::size_type nsends; - pack(intersection, view_g, fd_m, nghost_g, ldom_g, nsends); + pack(intersection, view_g, fd_m, nghost_g, ldom_g, nsends); buffer_type buf = Ippl::Comm->getBuffer(IPPL_VICO_SEND + 2 * 8 + i, nsends); @@ -1210,7 +1214,7 @@ namespace ippl { requests.resize(requests.size() + 1); Communicate::size_type nsends; - pack(intersection, view_g, fd_m, nghost_g, ldom_g, nsends); + pack(intersection, view_g, fd_m, nghost_g, ldom_g, nsends); buffer_type buf = Ippl::Comm->getBuffer(IPPL_VICO_SEND + 3 * 8 + i, nsends); @@ -1240,7 +1244,7 @@ namespace ippl { requests.resize(requests.size() + 1); Communicate::size_type nsends; - pack(intersection, view_g, fd_m, nghost_g, ldom_g, nsends); + pack(intersection, view_g, fd_m, nghost_g, ldom_g, nsends); buffer_type buf = Ippl::Comm->getBuffer(IPPL_VICO_SEND + 4 * 8 + i, nsends); @@ -1270,7 +1274,7 @@ namespace ippl { requests.resize(requests.size() + 1); Communicate::size_type nsends; - pack(intersection, view_g, fd_m, nghost_g, ldom_g, nsends); + pack(intersection, view_g, fd_m, nghost_g, ldom_g, nsends); buffer_type buf = Ippl::Comm->getBuffer(IPPL_VICO_SEND + 5 * 8 + i, nsends); @@ -1300,7 +1304,7 @@ namespace ippl { requests.resize(requests.size() + 1); Communicate::size_type nsends; - pack(intersection, view_g, fd_m, nghost_g, ldom_g, nsends); + pack(intersection, view_g, fd_m, nghost_g, ldom_g, nsends); buffer_type buf = Ippl::Comm->getBuffer(IPPL_VICO_SEND + 6 * 8 + i, nsends); @@ -1332,7 +1336,7 @@ namespace ippl { requests.resize(requests.size() + 1); Communicate::size_type nsends; - pack(intersection, view_g, fd_m, nghost_g, ldom_g, nsends); + pack(intersection, view_g, fd_m, nghost_g, ldom_g, nsends); buffer_type buf = Ippl::Comm->getBuffer(IPPL_VICO_SEND + 7 * 8 + i, nsends); @@ -1364,7 +1368,7 @@ namespace ippl { Ippl::Comm->recv(i, tag, fd_m, *buf, nrecvs * sizeof(double), nrecvs); buf->resetReadPos(); - unpack(intersection, view, fd_m, nghost, ldom); + unpack(intersection, view, fd_m, nghost, ldom); } } @@ -1397,7 +1401,7 @@ namespace ippl { Ippl::Comm->recv(i, tag, fd_m, *buf, nrecvs * sizeof(double), nrecvs); buf->resetReadPos(); - unpack(intersection, view, fd_m, nghost, ldom, true, false, false); + unpack(intersection, view, fd_m, nghost, ldom, true, false, false); } } @@ -1430,7 +1434,7 @@ namespace ippl { Ippl::Comm->recv(i, tag, fd_m, *buf, nrecvs * sizeof(double), nrecvs); buf->resetReadPos(); - unpack(intersection, view, fd_m, nghost, ldom, false, true, false); + unpack(intersection, view, fd_m, nghost, ldom, false, true, false); } } @@ -1463,7 +1467,7 @@ namespace ippl { Ippl::Comm->recv(i, tag, fd_m, *buf, nrecvs * sizeof(double), nrecvs); buf->resetReadPos(); - unpack(intersection, view, fd_m, nghost, ldom, false, false, true); + unpack(intersection, view, fd_m, nghost, ldom, false, false, true); } } @@ -1500,7 +1504,7 @@ namespace ippl { Ippl::Comm->recv(i, tag, fd_m, *buf, nrecvs * sizeof(double), nrecvs); buf->resetReadPos(); - unpack(intersection, view, fd_m, nghost, ldom, true, true, false); + unpack(intersection, view, fd_m, nghost, ldom, true, true, false); } } @@ -1537,7 +1541,7 @@ namespace ippl { Ippl::Comm->recv(i, tag, fd_m, *buf, nrecvs * sizeof(double), nrecvs); buf->resetReadPos(); - unpack(intersection, view, fd_m, nghost, ldom, false, true, true); + unpack(intersection, view, fd_m, nghost, ldom, false, true, true); } } @@ -1574,7 +1578,7 @@ namespace ippl { Ippl::Comm->recv(i, tag, fd_m, *buf, nrecvs * sizeof(double), nrecvs); buf->resetReadPos(); - unpack(intersection, view, fd_m, nghost, ldom, true, false, true); + unpack(intersection, view, fd_m, nghost, ldom, true, false, true); } } @@ -1615,7 +1619,7 @@ namespace ippl { Ippl::Comm->recv(i, tag, fd_m, *buf, nrecvs * sizeof(double), nrecvs); buf->resetReadPos(); - unpack(intersection, view, fd_m, nghost, ldom, true, true, true); + unpack(intersection, view, fd_m, nghost, ldom, true, true, true); } } } diff --git a/src/Solver/PCG.h b/src/Solver/PCG.h index d592f58bd..e46bf6e69 100644 --- a/src/Solver/PCG.h +++ b/src/Solver/PCG.h @@ -25,7 +25,7 @@ namespace ippl { template - class PCG : public SolverAlgorithm { + class PCG : public SolverAlgorithm { public: using Base = SolverAlgorithm; using typename Base::lhs_type; @@ -58,7 +58,7 @@ namespace ippl { // https://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf lhs_type r(mesh, layout), d(mesh, layout); - using bc_type = BConds; + using bc_type = BConds; bc_type lhsBCs = lhs.getFieldBC(); bc_type bc; @@ -67,11 +67,11 @@ namespace ippl { FieldBC bcType = lhsBCs[i]->getBCType(); if (bcType == PERIODIC_FACE) { // If the LHS has periodic BCs, so does the residue - bc[i] = std::make_shared>(i); + bc[i] = std::make_shared>(i); } else if (bcType & CONSTANT_FACE) { // If the LHS has constant BCs, the residue is zero on the BCs // Bitwise AND with CONSTANT_FACE will succeed for ZeroFace or ConstantFace - bc[i] = std::make_shared>(i); + bc[i] = std::make_shared>(i); allFacesPeriodic = false; } else { throw IpplException("PCG::operator()", diff --git a/src/Solver/test/Budiardja_plot.cpp b/src/Solver/test/Budiardja_plot.cpp index fdf62f0d4..494565fa5 100644 --- a/src/Solver/test/Budiardja_plot.cpp +++ b/src/Solver/test/Budiardja_plot.cpp @@ -34,6 +34,9 @@ int main(int argc, char* argv[]) { // number of interations const int n = 5; + using Mesh_t = ippl::UniformCartesian; + using Centering_t = Mesh_t::DefaultCentering; + // number of gridpoints to iterate over std::array N = {48, 144, 288, 384, 576}; @@ -54,13 +57,13 @@ int main(int argc, char* argv[]) { double dx = 2.4 / pt; ippl::Vector hx = {dx, dx, dx}; ippl::Vector origin = {0.0, 0.0, 0.0}; - ippl::UniformCartesian mesh(owned, hx, origin); + Mesh_t mesh(owned, hx, origin); // all parallel layout, standard domain, normal axis order ippl::FieldLayout<3> layout(owned, decomp); // define the L (phi) and R (rho) fields - typedef ippl::Field field; + typedef ippl::Field field; field rho; rho.initialize(mesh, layout); @@ -122,8 +125,9 @@ int main(int argc, char* argv[]) { fftParams.add("r2c_direction", 0); // define an FFTPoissonSolver object - ippl::FFTPoissonSolver, double, 3> FFTsolver(rho, fftParams, - "HOCKNEY"); + ippl::FFTPoissonSolver, double, 3, Mesh_t, Centering_t> FFTsolver(rho, + fftParams, + "HOCKNEY"); // solve the Poisson equation -> rho contains the solution (phi) now FFTsolver.solve(); diff --git a/src/Solver/test/TestCGSolver.cpp b/src/Solver/test/TestCGSolver.cpp index 3672dd53a..ebb3e4721 100644 --- a/src/Solver/test/TestCGSolver.cpp +++ b/src/Solver/test/TestCGSolver.cpp @@ -17,6 +17,8 @@ int main(int argc, char* argv[]) { Ippl ippl(argc, argv); constexpr unsigned int dim = 3; + using Mesh_t = ippl::UniformCartesian; + using Centering_t = Mesh_t::DefaultCentering; int pt = 4, ptY = 4; bool isWeak = false; @@ -55,19 +57,19 @@ int main(int argc, char* argv[]) { double dy = 2.0 / double(ptY); ippl::Vector hx = {dx, dy, dx}; ippl::Vector origin = {-1, -1, -1}; - ippl::UniformCartesian mesh(owned, hx, origin); + Mesh_t mesh(owned, hx, origin); double pi = acos(-1.0); - typedef ippl::Field field_type; + typedef ippl::Field field_type; field_type rhs(mesh, layout), lhs(mesh, layout), solution(mesh, layout); - typedef ippl::BConds bc_type; + typedef ippl::BConds bc_type; bc_type bcField; for (unsigned int i = 0; i < 6; ++i) { - bcField[i] = std::make_shared>(i); + bcField[i] = std::make_shared>(i); } lhs.setFieldBC(bcField); @@ -113,7 +115,7 @@ int main(int argc, char* argv[]) { * sin(sin(pi * z))); }); - ippl::ElectrostaticsCG lapsolver; + ippl::ElectrostaticsCG lapsolver; ippl::ParameterList params; params.add("max_iterations", 2000); diff --git a/src/Solver/test/TestFFTPeriodicPoissonSolver.cpp b/src/Solver/test/TestFFTPeriodicPoissonSolver.cpp index b197c809f..3e196d2b1 100644 --- a/src/Solver/test/TestFFTPeriodicPoissonSolver.cpp +++ b/src/Solver/test/TestFFTPeriodicPoissonSolver.cpp @@ -9,6 +9,8 @@ int main(int argc, char* argv[]) { Ippl ippl(argc, argv); constexpr unsigned int dim = 3; + using Mesh_t = ippl::UniformCartesian; + using Centering_t = Mesh_t::DefaultCentering; const int npts = 7; std::array pts = {2, 4, 8, 16, 32, 64, 128}; @@ -35,18 +37,18 @@ int main(int argc, char* argv[]) { double dx = 2.0 / double(pt); ippl::Vector hx = {dx, dx, dx}; ippl::Vector origin = {-1.0, -1.0, -1.0}; - ippl::UniformCartesian mesh(owned, hx, origin); + Mesh_t mesh(owned, hx, origin); double pi = acos(-1.0); - typedef ippl::Field Field_t; + typedef ippl::Field Field_t; typedef ippl::Vector Vector_t; - typedef ippl::Field VField_t; + typedef ippl::Field VField_t; Field_t field; field.initialize(mesh, layout); - typedef ippl::FFTPeriodicPoissonSolver Solver_t; + typedef ippl::FFTPeriodicPoissonSolver Solver_t; ippl::ParameterList params; params.add("output_type", Solver_t::SOL); diff --git a/src/Solver/test/TestGaussian.cpp b/src/Solver/test/TestGaussian.cpp index d9d532c3c..d65c7d423 100644 --- a/src/Solver/test/TestGaussian.cpp +++ b/src/Solver/test/TestGaussian.cpp @@ -56,6 +56,9 @@ int main(int argc, char* argv[]) { const unsigned int Dim = 3; + using Mesh_t = ippl::UniformCartesian; + using Centering_t = Mesh_t::DefaultCentering; + // start a timer static IpplTimings::TimerRef allTimer = IpplTimings::getTimer("allTimer"); IpplTimings::startTimer(allTimer); @@ -94,19 +97,19 @@ int main(int argc, char* argv[]) { double dz = 1.0 / nr[2]; ippl::Vector hr = {dx, dy, dz}; ippl::Vector origin = {0.0, 0.0, 0.0}; - ippl::UniformCartesian mesh(owned, hr, origin); + Mesh_t mesh(owned, hr, origin); // all parallel layout, standard domain, normal axis order ippl::FieldLayout layout(owned, decomp); // define the R (rho) field - typedef ippl::Field field; + typedef ippl::Field field; field exact, rho; exact.initialize(mesh, layout); rho.initialize(mesh, layout); // define the Vector field E (LHS) - typedef ippl::Field, Dim> fieldV; + typedef ippl::Field, Dim, Mesh_t, Centering_t> fieldV; fieldV exactE, fieldE; exactE.initialize(mesh, layout); fieldE.initialize(mesh, layout); @@ -217,8 +220,10 @@ int main(int argc, char* argv[]) { fftParams.add("r2c_direction", 0); // define an FFTPoissonSolver object - ippl::FFTPoissonSolver, double, Dim> FFTsolver(fieldE, rho, fftParams, - algorithm); + ippl::FFTPoissonSolver, double, Dim, Mesh_t, Centering_t> FFTsolver(fieldE, + rho, + fftParams, + algorithm); // iterate over 5 timesteps for (int times = 0; times < 5; ++times) { diff --git a/src/Solver/test/TestGaussian_biharmonic.cpp b/src/Solver/test/TestGaussian_biharmonic.cpp index d417b3014..143e61677 100644 --- a/src/Solver/test/TestGaussian_biharmonic.cpp +++ b/src/Solver/test/TestGaussian_biharmonic.cpp @@ -10,6 +10,11 @@ #include "FFTPoissonSolver.h" +using Mesh_t = ippl::UniformCartesian; +using Centering_t = Mesh_t::DefaultCentering; +using ScalarField_t = ippl::Field; +using VectorField_t = ippl::Field, 3, Mesh_t, Centering_t>; + KOKKOS_INLINE_FUNCTION double gaussian(double x, double y, double z, double sigma = 0.05, double mu = 0.5) { double pi = std::acos(-1.0); @@ -45,9 +50,9 @@ KOKKOS_INLINE_FUNCTION ippl::Vector exact_grad(double x, double y, do } // Define vtk dump function for plotting the fields -void dumpVTK(std::string path, ippl::Field& rho, int nx, int ny, int nz, int iteration, +void dumpVTK(std::string path, ScalarField_t& rho, int nx, int ny, int nz, int iteration, double dx, double dy, double dz) { - typename ippl::Field::view_type::host_mirror_type host_view = rho.getHostMirror(); + typename ScalarField_t::view_type::host_mirror_type host_view = rho.getHostMirror(); Kokkos::deep_copy(host_view, rho.getView()); std::ofstream vtkout; vtkout.precision(10); @@ -128,24 +133,22 @@ int main(int argc, char* argv[]) { ippl::FieldLayout<3> layout(owned, decomp); // define the R (rho) field - typedef ippl::Field field; - field rho; + ScalarField_t rho; rho.initialize(mesh, layout); // define the exact solution field - field exact; + ScalarField_t exact; exact.initialize(mesh, layout); // field for gradient and exact gradient - typedef ippl::Field, 3> fieldV; - fieldV fieldE, exactE; + VectorField_t fieldE, exactE; fieldE.initialize(mesh, layout); exactE.initialize(mesh, layout); // assign the rho field with a gaussian - typename field::view_type view_rho = rho.getView(); - const int nghost = rho.getNghost(); - const auto& ldom = layout.getLocalNDIndex(); + typename ScalarField_t::view_type view_rho = rho.getView(); + const int nghost = rho.getNghost(); + const auto& ldom = layout.getLocalNDIndex(); Kokkos::parallel_for( "Assign rho field", @@ -167,7 +170,7 @@ int main(int argc, char* argv[]) { }); // assign the exact field with its values (erf function) - typename field::view_type view_exact = exact.getView(); + typename ScalarField_t::view_type view_exact = exact.getView(); Kokkos::parallel_for( "Assign exact field", @@ -220,8 +223,10 @@ int main(int argc, char* argv[]) { fftParams.add("r2c_direction", 0); // define an FFTPoissonSolver object - ippl::FFTPoissonSolver, double, 3> FFTsolver(fieldE, rho, fftParams, - algorithm); + ippl::FFTPoissonSolver, double, 3, Mesh_t, Centering_t> FFTsolver(fieldE, + rho, + fftParams, + algorithm); // solve the Poisson equation -> rho contains the solution (phi) now FFTsolver.solve(); diff --git a/src/Solver/test/TestGaussian_convergence.cpp b/src/Solver/test/TestGaussian_convergence.cpp index af4cb0e46..5a8bc0920 100644 --- a/src/Solver/test/TestGaussian_convergence.cpp +++ b/src/Solver/test/TestGaussian_convergence.cpp @@ -10,6 +10,11 @@ #include "FFTPoissonSolver.h" +using Mesh_t = ippl::UniformCartesian; +using Centering_t = Mesh_t::DefaultCentering; +using ScalarField_t = ippl::Field; +using VectorField_t = ippl::Field, 3, Mesh_t, Centering_t>; + KOKKOS_INLINE_FUNCTION double gaussian(double x, double y, double z, double sigma = 0.05, double mu = 0.5) { double pi = std::acos(-1.0); @@ -40,9 +45,9 @@ KOKKOS_INLINE_FUNCTION ippl::Vector exact_E(double x, double y, doubl } // Define vtk dump function for plotting the fields -void dumpVTK(std::string path, ippl::Field& rho, int nx, int ny, int nz, int iteration, +void dumpVTK(std::string path, ScalarField_t& rho, int nx, int ny, int nz, int iteration, double dx, double dy, double dz) { - typename ippl::Field::view_type::host_mirror_type host_view = rho.getHostMirror(); + typename ScalarField_t::view_type::host_mirror_type host_view = rho.getHostMirror(); Kokkos::deep_copy(host_view, rho.getView()); std::ofstream vtkout; vtkout.precision(10); @@ -117,31 +122,28 @@ int main(int argc, char* argv[]) { double dx = 1.0 / pt; ippl::Vector hx = {dx, dx, dx}; ippl::Vector origin = {0.0, 0.0, 0.0}; - ippl::UniformCartesian mesh(owned, hx, origin); + Mesh_t mesh(owned, hx, origin); // all parallel layout, standard domain, normal axis order ippl::FieldLayout<3> layout(owned, decomp); // define the R (rho) field - typedef ippl::Field field; - field rho; + ScalarField_t rho; rho.initialize(mesh, layout); // define the exact solution field - field exact; + ScalarField_t exact; exact.initialize(mesh, layout); // define the Vector field E and the exact E field - typedef ippl::Field, 3> fieldV; - - fieldV exactE, fieldE; + VectorField_t exactE, fieldE; exactE.initialize(mesh, layout); fieldE.initialize(mesh, layout); // assign the rho field with a gaussian - typename field::view_type view_rho = rho.getView(); - const int nghost = rho.getNghost(); - const auto& ldom = layout.getLocalNDIndex(); + typename ScalarField_t::view_type view_rho = rho.getView(); + const int nghost = rho.getNghost(); + const auto& ldom = layout.getLocalNDIndex(); Kokkos::parallel_for( "Assign rho field", @@ -163,7 +165,7 @@ int main(int argc, char* argv[]) { }); // assign the exact field with its values (erf function) - typename field::view_type view_exact = exact.getView(); + typename ScalarField_t::view_type view_exact = exact.getView(); Kokkos::parallel_for( "Assign exact field", @@ -215,8 +217,10 @@ int main(int argc, char* argv[]) { fftParams.add("comm", ippl::a2av); fftParams.add("r2c_direction", 0); // define an FFTPoissonSolver object - ippl::FFTPoissonSolver, double, 3> FFTsolver(fieldE, rho, fftParams, - algorithm); + ippl::FFTPoissonSolver, double, 3, Mesh_t, Centering_t> FFTsolver(fieldE, + rho, + fftParams, + algorithm); // solve the Poisson equation -> rho contains the solution (phi) now FFTsolver.solve(); diff --git a/src/Solver/test/TestSolverDesign.cpp b/src/Solver/test/TestSolverDesign.cpp index b4448cc96..3ff22ccef 100644 --- a/src/Solver/test/TestSolverDesign.cpp +++ b/src/Solver/test/TestSolverDesign.cpp @@ -8,8 +8,10 @@ #include "Electrostatics.h" constexpr unsigned int dim = 3; +using Mesh_t = ippl::UniformCartesian; +using Centering_t = Mesh_t::DefaultCentering; -class TestSolver : public ippl::Electrostatics { +class TestSolver : public ippl::Electrostatics { public: void solve() override { *rhs_mp = *lhs_mp + *rhs_mp; @@ -38,12 +40,12 @@ int main(int argc, char* argv[]) { double dx = 1.0 / double(pt); ippl::Vector hx = {dx, dx, dx}; ippl::Vector origin = {0, 0, 0}; - ippl::UniformCartesian mesh(owned, hx, origin); + Mesh_t mesh(owned, hx, origin); - typedef ippl::Field field_type; + typedef ippl::Field field_type; field_type lhs(mesh, layout), rhs(mesh, layout); - typedef ippl::Field, dim> vfield_type; + typedef ippl::Field, dim, Mesh_t, Centering_t> vfield_type; vfield_type grad(mesh, layout); lhs = 1.0; diff --git a/src/Solver/test/TestSphere.cpp b/src/Solver/test/TestSphere.cpp index 002710a07..00e2ade3d 100644 --- a/src/Solver/test/TestSphere.cpp +++ b/src/Solver/test/TestSphere.cpp @@ -55,17 +55,21 @@ int main(int argc, char* argv[]) { for (unsigned int d = 0; d < 3; d++) decomp[d] = ippl::PARALLEL; + using Mesh_t = ippl::UniformCartesian; + using Centering_t = Mesh_t::DefaultCentering; + using Vector_t = ippl::Vector; + // unit box - double dx = 2.4 / pt; - ippl::Vector hx = {dx, dx, dx}; - ippl::Vector origin = {0.0, 0.0, 0.0}; - ippl::UniformCartesian mesh(owned, hx, origin); + double dx = 2.4 / pt; + Vector_t hx = {dx, dx, dx}; + Vector_t origin = {0.0, 0.0, 0.0}; + Mesh_t mesh(owned, hx, origin); // all parallel layout, standard domain, normal axis order ippl::FieldLayout<3> layout(owned, decomp); // define the L (phi) and R (rho) fields - typedef ippl::Field field; + typedef ippl::Field field; field rho; rho.initialize(mesh, layout); @@ -126,8 +130,8 @@ int main(int argc, char* argv[]) { fftParams.add("comm", ippl::a2av); fftParams.add("r2c_direction", 0); - ippl::FFTPoissonSolver, double, 3> FFTsolver(rho, fftParams, - algorithm); + ippl::FFTPoissonSolver FFTsolver(rho, fftParams, + algorithm); // solve the Poisson equation -> rho contains the solution (phi) now FFTsolver.solve(); diff --git a/test/FFT/TestCos.cpp b/test/FFT/TestCos.cpp index aa6c96779..8d9253d31 100644 --- a/test/FFT/TestCos.cpp +++ b/test/FFT/TestCos.cpp @@ -11,6 +11,8 @@ int main(int argc, char* argv[]) { Ippl ippl(argc, argv); constexpr unsigned int dim = 3; + using Mesh_t = ippl::UniformCartesian; + using Centering_t = Mesh_t::DefaultCentering; std::array pt = {32, 32, 32}; ippl::Index I(pt[0]); @@ -31,9 +33,9 @@ int main(int argc, char* argv[]) { }; ippl::Vector hx = {dx[0], dx[1], dx[2]}; ippl::Vector origin = {0, 0, 0}; - ippl::UniformCartesian mesh(owned, hx, origin); + Mesh_t mesh(owned, hx, origin); - typedef ippl::Field field_type; + typedef ippl::Field field_type; field_type field(mesh, layout); @@ -44,7 +46,7 @@ int main(int argc, char* argv[]) { fftParams.add("use_gpu_aware", true); fftParams.add("comm", ippl::p2p_pl); - typedef ippl::FFT FFT_type; + typedef ippl::FFT FFT_type; std::unique_ptr fft; diff --git a/test/FFT/TestFFTCC.cpp b/test/FFT/TestFFTCC.cpp index 5688f2eba..e06b9ff36 100644 --- a/test/FFT/TestFFTCC.cpp +++ b/test/FFT/TestFFTCC.cpp @@ -11,6 +11,8 @@ int main(int argc, char* argv[]) { Ippl ippl(argc, argv); constexpr unsigned int dim = 3; + using Mesh_t = ippl::UniformCartesian; + using Centering_t = Mesh_t::DefaultCentering; std::array pt = {32, 32, 32}; ippl::Index I(pt[0]); @@ -31,9 +33,9 @@ int main(int argc, char* argv[]) { }; ippl::Vector hx = {dx[0], dx[1], dx[2]}; ippl::Vector origin = {0, 0, 0}; - ippl::UniformCartesian mesh(owned, hx, origin); + Mesh_t mesh(owned, hx, origin); - typedef ippl::Field, dim> field_type; + typedef ippl::Field, dim, Mesh_t, Centering_t> field_type; field_type field(mesh, layout); @@ -41,7 +43,7 @@ int main(int argc, char* argv[]) { fftParams.add("use_heffte_defaults", true); - typedef ippl::FFT FFT_type; + typedef ippl::FFT FFT_type; std::unique_ptr fft; diff --git a/test/FFT/TestFFTRC.cpp b/test/FFT/TestFFTRC.cpp index 57a6c1138..31c608bdd 100644 --- a/test/FFT/TestFFTRC.cpp +++ b/test/FFT/TestFFTRC.cpp @@ -11,6 +11,8 @@ int main(int argc, char* argv[]) { Ippl ippl(argc, argv); constexpr unsigned int dim = 3; + using Mesh_t = ippl::UniformCartesian; + using Centering_t = Mesh_t::DefaultCentering; std::array pt = {64, 64, 64}; ippl::Index Iinput(pt[0]); @@ -33,8 +35,8 @@ int main(int argc, char* argv[]) { ippl::Vector origin = {0, 0, 0}; ippl::UniformCartesian meshInput(ownedInput, hx, origin); - typedef ippl::Field, dim> field_type_complex; - typedef ippl::Field field_type_real; + typedef ippl::Field, dim, Mesh_t, Centering_t> field_type_complex; + typedef ippl::Field field_type_real; field_type_real fieldInput(meshInput, layoutInput); @@ -65,10 +67,10 @@ int main(int argc, char* argv[]) { } ippl::FieldLayout layoutOutput(ownedOutput, allParallel); - ippl::UniformCartesian meshOutput(ownedOutput, hx, origin); + Mesh_t meshOutput(ownedOutput, hx, origin); field_type_complex fieldOutput(meshOutput, layoutOutput); - typedef ippl::FFT FFT_type; + typedef ippl::FFT FFT_type; std::unique_ptr fft; diff --git a/test/FFT/TestSine.cpp b/test/FFT/TestSine.cpp index dac5393d2..5b5e9784e 100644 --- a/test/FFT/TestSine.cpp +++ b/test/FFT/TestSine.cpp @@ -11,6 +11,8 @@ int main(int argc, char* argv[]) { Ippl ippl(argc, argv); constexpr unsigned int dim = 3; + using Mesh_t = ippl::UniformCartesian; + using Centering_t = Mesh_t::DefaultCentering; std::array pt = {32, 32, 32}; ippl::Index I(pt[0]); @@ -31,9 +33,9 @@ int main(int argc, char* argv[]) { }; ippl::Vector hx = {dx[0], dx[1], dx[2]}; ippl::Vector origin = {0, 0, 0}; - ippl::UniformCartesian mesh(owned, hx, origin); + Mesh_t mesh(owned, hx, origin); - typedef ippl::Field field_type; + typedef ippl::Field field_type; field_type field(mesh, layout); @@ -45,7 +47,7 @@ int main(int argc, char* argv[]) { fftParams.add("use_gpu_aware", true); fftParams.add("comm", ippl::p2p_pl); - typedef ippl::FFT FFT_type; + typedef ippl::FFT FFT_type; std::unique_ptr fft; diff --git a/test/field/TestCurl.cpp b/test/field/TestCurl.cpp index 9750a0172..f87090f07 100644 --- a/test/field/TestCurl.cpp +++ b/test/field/TestCurl.cpp @@ -44,6 +44,8 @@ int main(int argc, char* argv[]) { Ippl ippl(argc, argv); constexpr unsigned int dim = 3; + using Mesh_t = ippl::UniformCartesian; + using Centering_t = Mesh_t::DefaultCentering; int pt = std::atoi(argv[1]); bool gauss_fct = std::atoi(argv[2]); @@ -61,10 +63,10 @@ int main(int argc, char* argv[]) { double dx = 1.0 / double(pt); ippl::Vector hx = {dx, dx, dx}; ippl::Vector origin = {0.0, 0.0, 0.0}; - ippl::UniformCartesian mesh(owned, hx, origin); + Mesh_t mesh(owned, hx, origin); typedef ippl::Vector Vector_t; - typedef ippl::Field Vfield_t; + typedef ippl::Field Vfield_t; Vfield_t vfield(mesh, layout); Vfield_t result(mesh, layout); diff --git a/test/field/TestFieldBC.cpp b/test/field/TestFieldBC.cpp index 5bd91561b..e25fddb00 100644 --- a/test/field/TestFieldBC.cpp +++ b/test/field/TestFieldBC.cpp @@ -25,27 +25,28 @@ int main(int argc, char* argv[]) { ippl::Vector origin = {0, 0, 0}; using Mesh_t = ippl::UniformCartesian; + using Centering_t = Mesh_t::DefaultCentering; Mesh_t mesh(owned, hx, origin); - typedef ippl::Field field_type; + typedef ippl::Field field_type; - typedef ippl::BConds bc_type; + typedef ippl::BConds bc_type; bc_type bcField; // X direction periodic BC for (unsigned int i = 0; i < 2; ++i) { - bcField[i] = std::make_shared>(i); + bcField[i] = std::make_shared>(i); } ////Lower Y face - bcField[2] = std::make_shared>(2); + bcField[2] = std::make_shared>(2); ////Higher Y face - bcField[3] = std::make_shared>(3, 7.0); + bcField[3] = std::make_shared>(3, 7.0); ////Lower Z face - bcField[4] = std::make_shared>(4); + bcField[4] = std::make_shared>(4); ////Higher Z face - bcField[5] = std::make_shared>(5, 0.0, 1.0); + bcField[5] = std::make_shared>(5, 0.0, 1.0); // std::cout << bcField << std::endl; std::cout << layout << std::endl; diff --git a/test/field/TestHalo.cpp b/test/field/TestHalo.cpp index dde58cd46..c0437d189 100644 --- a/test/field/TestHalo.cpp +++ b/test/field/TestHalo.cpp @@ -14,6 +14,8 @@ int main(int argc, char* argv[]) { IpplTimings::startTimer(mainTimer); constexpr unsigned int dim = 3; + using Mesh_t = ippl::UniformCartesian; + using Centering_t = Mesh_t::DefaultCentering; // std::array pt = {8, 7, 13}; std::array pt = {4, 4, 4}; @@ -36,9 +38,9 @@ int main(int argc, char* argv[]) { }; ippl::Vector hx = {dx[0], dx[1], dx[2]}; ippl::Vector origin = {0, 0, 0}; - ippl::UniformCartesian mesh(owned, hx, origin); + Mesh_t mesh(owned, hx, origin); - typedef ippl::Field field_type; + typedef ippl::Field field_type; field_type field(mesh, layout); diff --git a/test/field/TestHessian.cpp b/test/field/TestHessian.cpp index 6ab95928a..66d58e728 100644 --- a/test/field/TestHessian.cpp +++ b/test/field/TestHessian.cpp @@ -43,6 +43,8 @@ int main(int argc, char* argv[]) { Ippl ippl(argc, argv); constexpr unsigned int dim = 3; + using Mesh_t = ippl::UniformCartesian; + using Centering_t = Mesh_t::DefaultCentering; int pt = std::atoi(argv[1]); bool gauss_fct = std::atoi(argv[2]); @@ -59,15 +61,15 @@ int main(int argc, char* argv[]) { // type definitions typedef ippl::Vector Vector_t; - typedef ippl::Field Field_t; + typedef ippl::Field Field_t; typedef ippl::Vector Matrix_t; - typedef ippl::Field MField_t; + typedef ippl::Field MField_t; // domain [0,1]^3 double dx = 1.0 / double(pt); Vector_t hx = {dx, dx, dx}; Vector_t origin = {0.0, 0.0, 0.0}; - ippl::UniformCartesian mesh(owned, hx, origin); + Mesh_t mesh(owned, hx, origin); Field_t field(mesh, layout, 1); MField_t result(mesh, layout, 1); diff --git a/test/field/TestLaplace.cpp b/test/field/TestLaplace.cpp index ae94f01f6..911cb8a55 100644 --- a/test/field/TestLaplace.cpp +++ b/test/field/TestLaplace.cpp @@ -9,6 +9,8 @@ int main(int argc, char* argv[]) { Ippl ippl(argc, argv); constexpr unsigned int dim = 3; + using Mesh_t = ippl::UniformCartesian; + using Centering_t = Mesh_t::DefaultCentering; int pt = 4; ippl::Index I(pt); @@ -26,34 +28,33 @@ int main(int argc, char* argv[]) { double dx = 2.0 / double(pt); ippl::Vector hx = {dx, dx, dx}; ippl::Vector origin = {-1.0, -1.0, -1.0}; - ippl::UniformCartesian mesh(owned, hx, origin); + Mesh_t mesh(owned, hx, origin); double pi = acos(-1.0); - typedef ippl::Field field_type; - typedef ippl::Field, dim> vector_field_type; + typedef ippl::Field Field_t; + typedef ippl::Field, dim, Mesh_t, Centering_t> vector_field_type; typedef ippl::Vector Vector_t; - field_type field(mesh, layout); - field_type Lap(mesh, layout); - field_type Lap_exact(mesh, layout); + Field_t field(mesh, layout); + Field_t Lap(mesh, layout); + Field_t Lap_exact(mesh, layout); vector_field_type vfield(mesh, layout); - typedef ippl::Field Field_t; typename Field_t::view_type& view = field.getView(); typename Field_t::view_type& view_exact = Lap_exact.getView(); - typedef ippl::BConds bc_type; - typedef ippl::BConds vbc_type; + typedef ippl::BConds bc_type; + typedef ippl::BConds vbc_type; bc_type bcField; vbc_type vbcField; // X direction periodic BC for (unsigned int i = 0; i < 6; ++i) { - bcField[i] = std::make_shared>(i); - vbcField[i] = std::make_shared>(i); + bcField[i] = std::make_shared>(i); + vbcField[i] = std::make_shared>(i); } ////Lower Y face // bcField[2] = std::make_shared>(2); diff --git a/test/particle/PIC3d.cpp b/test/particle/PIC3d.cpp index 1261ed894..673e668a6 100644 --- a/test/particle/PIC3d.cpp +++ b/test/particle/PIC3d.cpp @@ -41,13 +41,14 @@ constexpr unsigned Dim = 3; typedef ippl::ParticleSpatialLayout PLayout_t; typedef ippl::UniformCartesian Mesh_t; typedef ippl::FieldLayout FieldLayout_t; -typedef ippl::OrthogonalRecursiveBisection ORB; +typedef Mesh_t::DefaultCentering Centering_t; +typedef ippl::OrthogonalRecursiveBisection ORB; template using Vector = ippl::Vector; template -using Field = ippl::Field; +using Field = ippl::Field; template using ParticleAttrib = ippl::ParticleAttrib; diff --git a/test/particle/TestGather.cpp b/test/particle/TestGather.cpp index 5b3067e86..0a4d8b1b0 100644 --- a/test/particle/TestGather.cpp +++ b/test/particle/TestGather.cpp @@ -20,6 +20,8 @@ int main(int argc, char* argv[]) { typedef ippl::ParticleSpatialLayout playout_type; typedef Bunch bunch_type; + using Mesh_t = ippl::UniformCartesian; + using Centering_t = Mesh_t::DefaultCentering; int pt = 512; ippl::Index I(pt); @@ -35,7 +37,7 @@ int main(int argc, char* argv[]) { double dx = 1.0 / double(pt); ippl::Vector hx = {dx, dx, dx}; ippl::Vector origin = {0, 0, 0}; - ippl::UniformCartesian mesh(owned, hx, origin); + Mesh_t mesh(owned, hx, origin); playout_type pl(layout, mesh); @@ -60,7 +62,7 @@ int main(int argc, char* argv[]) { bunch_type buffer(pl); pl.update(bunch, buffer); - typedef ippl::Field field_type; + typedef ippl::Field field_type; field_type field; diff --git a/test/particle/TestScatter.cpp b/test/particle/TestScatter.cpp index 57dea32bd..efb0421aa 100644 --- a/test/particle/TestScatter.cpp +++ b/test/particle/TestScatter.cpp @@ -20,6 +20,8 @@ int main(int argc, char* argv[]) { typedef ippl::ParticleSpatialLayout playout_type; typedef Bunch bunch_type; + using Mesh_t = ippl::UniformCartesian; + using Centering_t = Mesh_t::DefaultCentering; int pt = 512; ippl::Index I(pt); @@ -35,12 +37,12 @@ int main(int argc, char* argv[]) { double dx = 1.0 / double(pt); ippl::Vector hx = {dx, dx, dx}; ippl::Vector origin = {0, 0, 0}; - ippl::UniformCartesian mesh(owned, hx, origin); + Mesh_t mesh(owned, hx, origin); playout_type pl(layout, mesh); bunch_type bunch(pl); - typedef ippl::Field field_type; + typedef ippl::Field field_type; field_type field; From 345b8dd6574730083323689c01c7c13639e13743 Mon Sep 17 00:00:00 2001 From: Matthias Frey Date: Sun, 19 Mar 2023 11:58:16 +0000 Subject: [PATCH 08/12] update unit-tests --- unit_tests/Field/Field.cpp | 24 +++++++++++++----------- unit_tests/Field/FieldBC.cpp | 14 ++++++++------ unit_tests/PIC/ORB.cpp | 33 +++++++++++++++++---------------- unit_tests/PIC/PIC.cpp | 27 ++++++++++++++------------- 4 files changed, 52 insertions(+), 46 deletions(-) diff --git a/unit_tests/Field/Field.cpp b/unit_tests/Field/Field.cpp index 77a880aa7..299a5fe68 100644 --- a/unit_tests/Field/Field.cpp +++ b/unit_tests/Field/Field.cpp @@ -24,9 +24,12 @@ class FieldTest : public ::testing::Test { public: static constexpr size_t dim = 3; - typedef ippl::Field field_type; - typedef ippl::UniformCartesian mesh_type; - typedef ippl::FieldLayout layout_type; + using Mesh_t = ippl::UniformCartesian; + using Centering_t = Mesh_t::DefaultCentering; + using field_type = ippl::Field ; + using bc_type = ippl::BConds; + using vector_field_type = ippl::Field, dim, Mesh_t, Centering_t>; + using layout_type = ippl::FieldLayout ; FieldTest() : nPoints(8) { @@ -46,13 +49,13 @@ class FieldTest : public ::testing::Test { double dx = 1.0 / double(nPoints); ippl::Vector hx = {dx, dx, dx}; ippl::Vector origin = {0, 0, 0}; - mesh = std::make_shared(owned, hx, origin); + mesh = std::make_shared(owned, hx, origin); field = std::make_unique(*mesh, *layout); } std::unique_ptr field; - std::shared_ptr mesh; + std::shared_ptr mesh; std::shared_ptr layout; size_t nPoints; }; @@ -144,7 +147,7 @@ TEST_F(FieldTest, VolumeIntegral2) { TEST_F(FieldTest, Grad) { *field = 1.; - ippl::Field, dim> vfield(*mesh, *layout); + vector_field_type vfield(*mesh, *layout); vfield = grad(*field); const int shift = vfield.getNghost(); @@ -164,7 +167,7 @@ TEST_F(FieldTest, Grad) { } TEST_F(FieldTest, Curl) { - ippl::Field, dim> vfield(*mesh, *layout); + vector_field_type vfield(*mesh, *layout); const int nghost = vfield.getNghost(); auto view_field = vfield.getView(); @@ -200,7 +203,7 @@ TEST_F(FieldTest, Curl) { Kokkos::deep_copy(view_field, mirror); - ippl::Field, dim> result(*mesh, *layout); + vector_field_type result(*mesh, *layout); result = curl(vfield); const int shift = result.getNghost(); @@ -221,9 +224,8 @@ TEST_F(FieldTest, Curl) { TEST_F(FieldTest, Hessian) { typedef ippl::Vector Vector_t; - typedef ippl::Field, dim> MField_t; - ippl::Field field(*mesh, *layout); + field_type field(*mesh, *layout); int nghost = field.getNghost(); auto view_field = field.getView(); @@ -253,7 +255,7 @@ TEST_F(FieldTest, Hessian) { Kokkos::deep_copy(view_field, mirror); - MField_t result(*mesh, *layout); + vector_field_type result(*mesh, *layout); result = hess(field); nghost = result.getNghost(); diff --git a/unit_tests/Field/FieldBC.cpp b/unit_tests/Field/FieldBC.cpp index f268cf2d0..7782101d3 100644 --- a/unit_tests/Field/FieldBC.cpp +++ b/unit_tests/Field/FieldBC.cpp @@ -27,8 +27,10 @@ class FieldBCTest : public ::testing::Test { public: static constexpr size_t dim = 3; - typedef ippl::Field field_type; - typedef ippl::BConds bc_type; + using Mesh_t = ippl::UniformCartesian; + using Centering_t = Mesh_t::DefaultCentering; + using field_type = ippl::Field ; + using bc_type = ippl::BConds; FieldBCTest() : nPoints(8) { @@ -133,7 +135,7 @@ TEST_F(FieldBCTest, PeriodicBC) { TEST_F(FieldBCTest, NoBC) { for (size_t i = 0; i < 2 * dim; ++i) { - bcField[i] = std::make_shared>(i); + bcField[i] = std::make_shared>(i); } bcField.findBCNeighbors(*field); bcField.apply(*field); @@ -143,7 +145,7 @@ TEST_F(FieldBCTest, NoBC) { TEST_F(FieldBCTest, ZeroBC) { for (size_t i = 0; i < 2 * dim; ++i) { - bcField[i] = std::make_shared>(i); + bcField[i] = std::make_shared>(i); } bcField.findBCNeighbors(*field); bcField.apply(*field); @@ -154,7 +156,7 @@ TEST_F(FieldBCTest, ZeroBC) { TEST_F(FieldBCTest, ConstantBC) { double constant = 7.0; for (size_t i = 0; i < 2 * dim; ++i) { - bcField[i] = std::make_shared>(i, constant); + bcField[i] = std::make_shared>(i, constant); } bcField.findBCNeighbors(*field); bcField.apply(*field); @@ -164,7 +166,7 @@ TEST_F(FieldBCTest, ConstantBC) { TEST_F(FieldBCTest, ExtrapolateBC) { for (size_t i = 0; i < 2 * dim; ++i) { - bcField[i] = std::make_shared>(i, 0.0, 1.0); + bcField[i] = std::make_shared>(i, 0.0, 1.0); } bcField.findBCNeighbors(*field); bcField.apply(*field); diff --git a/unit_tests/PIC/ORB.cpp b/unit_tests/PIC/ORB.cpp index dba3357cb..92f563566 100644 --- a/unit_tests/PIC/ORB.cpp +++ b/unit_tests/PIC/ORB.cpp @@ -26,11 +26,12 @@ class ORBTest : public ::testing::Test { public: static constexpr size_t dim = 3; - typedef ippl::Field field_type; - typedef ippl::FieldLayout flayout_type; - typedef ippl::UniformCartesian mesh_type; - typedef ippl::ParticleSpatialLayout playout_type; - typedef ippl::OrthogonalRecursiveBisection ORB; + using Mesh_t = ippl::UniformCartesian; + using Centering_t = Mesh_t::DefaultCentering; + using Field_t = ippl::Field; + using Flayout_t = ippl::FieldLayout; + using Playout_t = ippl::ParticleSpatialLayout; + using Orb_t = ippl::OrthogonalRecursiveBisection; template struct Bunch : public ippl::ParticleBase { @@ -44,13 +45,13 @@ class ORBTest : public ::testing::Test { typedef ippl::ParticleAttrib charge_container_type; charge_container_type Q; - void updateLayout(flayout_type fl, mesh_type mesh) { + void updateLayout(Flayout_t fl, Mesh_t mesh) { PLayout& layout = this->getLayout(); layout.updateLayout(fl, mesh); } }; - typedef Bunch bunch_type; + typedef Bunch bunch_type; ORBTest() // Original configuration 256^3 particles, 512^3 grid. @@ -68,17 +69,17 @@ class ORBTest : public ::testing::Test { allParallel[d] = ippl::PARALLEL; const bool isAllPeriodic = true; - layout_m = flayout_type(owned, allParallel, isAllPeriodic); + layout_m = Flayout_t(owned, allParallel, isAllPeriodic); double dx = 1.0 / double(nPoints); ippl::Vector hx = {dx, dx, dx}; ippl::Vector origin = {0, 0, 0}; - mesh_m = mesh_type(owned, hx, origin); + mesh_m = Mesh_t(owned, hx, origin); - field = std::make_unique(mesh_m, layout_m); + field = std::make_unique(mesh_m, layout_m); - pl_m = playout_type(layout_m, mesh_m); + pl_m = Playout_t(layout_m, mesh_m); bunch = std::make_unique(pl_m); @@ -119,15 +120,15 @@ class ORBTest : public ::testing::Test { ippl::NDIndex getDomain() { return layout_m.getDomain(); } - std::unique_ptr field; + std::unique_ptr field; std::unique_ptr bunch; size_t nParticles; size_t nPoints; - flayout_type layout_m; - mesh_type mesh_m; - playout_type pl_m; - ORB orb; + Flayout_t layout_m; + Mesh_t mesh_m; + Playout_t pl_m; + Orb_t orb; }; TEST_F(ORBTest, Volume) { diff --git a/unit_tests/PIC/PIC.cpp b/unit_tests/PIC/PIC.cpp index d74da07a5..04642cbba 100644 --- a/unit_tests/PIC/PIC.cpp +++ b/unit_tests/PIC/PIC.cpp @@ -25,10 +25,11 @@ class PICTest : public ::testing::Test { public: static constexpr size_t dim = 3; - typedef ippl::Field field_type; - typedef ippl::FieldLayout flayout_type; - typedef ippl::UniformCartesian mesh_type; - typedef ippl::ParticleSpatialLayout playout_type; + using Mesh_t = ippl::UniformCartesian; + using Centering_t = Mesh_t::DefaultCentering; + using Field_t = ippl::Field; + using Flayout_t = ippl::FieldLayout; + using Playout_t = ippl::ParticleSpatialLayout; template struct Bunch : public ippl::ParticleBase { @@ -43,7 +44,7 @@ class PICTest : public ::testing::Test { charge_container_type Q; }; - typedef Bunch bunch_type; + typedef Bunch bunch_type; PICTest() : nParticles(std::pow(256, 3)) @@ -59,17 +60,17 @@ class PICTest : public ::testing::Test { for (unsigned int d = 0; d < dim; d++) domDec[d] = ippl::PARALLEL; - layout_m = flayout_type(owned, domDec); + layout_m = Flayout_t(owned, domDec); double dx = 1.0 / double(nPoints); ippl::Vector hx = {dx, dx, dx}; ippl::Vector origin = {0, 0, 0}; - mesh_m = mesh_type(owned, hx, origin); + mesh_m = Mesh_t(owned, hx, origin); - field = std::make_unique(mesh_m, layout_m); + field = std::make_unique(mesh_m, layout_m); - pl = playout_type(layout_m, mesh_m); + pl = Playout_t(layout_m, mesh_m); bunch = std::make_unique(pl); @@ -99,15 +100,15 @@ class PICTest : public ::testing::Test { Kokkos::deep_copy(bunch->R.getView(), R_host); } - std::unique_ptr field; + std::unique_ptr field; std::unique_ptr bunch; size_t nParticles; size_t nPoints; - playout_type pl; + Playout_t pl; private: - flayout_type layout_m; - mesh_type mesh_m; + Flayout_t layout_m; + Mesh_t mesh_m; }; TEST_F(PICTest, Scatter) { From 8c017c80840fb73027e912747a835596fd0e9dcb Mon Sep 17 00:00:00 2001 From: Matthias Frey Date: Sun, 19 Mar 2023 18:32:37 +0000 Subject: [PATCH 09/12] fix typos --- unit_tests/Field/Field.cpp | 2 +- unit_tests/Field/FieldBC.cpp | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/unit_tests/Field/Field.cpp b/unit_tests/Field/Field.cpp index 299a5fe68..fd2d242d8 100644 --- a/unit_tests/Field/Field.cpp +++ b/unit_tests/Field/Field.cpp @@ -24,7 +24,7 @@ class FieldTest : public ::testing::Test { public: static constexpr size_t dim = 3; - using Mesh_t = ippl::UniformCartesian; + using Mesh_t = ippl::UniformCartesian; using Centering_t = Mesh_t::DefaultCentering; using field_type = ippl::Field ; using bc_type = ippl::BConds; diff --git a/unit_tests/Field/FieldBC.cpp b/unit_tests/Field/FieldBC.cpp index 7782101d3..3e12814f0 100644 --- a/unit_tests/Field/FieldBC.cpp +++ b/unit_tests/Field/FieldBC.cpp @@ -27,7 +27,7 @@ class FieldBCTest : public ::testing::Test { public: static constexpr size_t dim = 3; - using Mesh_t = ippl::UniformCartesian; + using Mesh_t = ippl::UniformCartesian; using Centering_t = Mesh_t::DefaultCentering; using field_type = ippl::Field ; using bc_type = ippl::BConds; @@ -125,7 +125,7 @@ class FieldBCTest : public ::testing::Test { TEST_F(FieldBCTest, PeriodicBC) { for (size_t i = 0; i < 2 * dim; ++i) { - bcField[i] = std::make_shared>(i); + bcField[i] = std::make_shared>(i); } bcField.findBCNeighbors(*field); bcField.apply(*field); From 4a20927f19f4ba83c84200ce720560a1d9723c0b Mon Sep 17 00:00:00 2001 From: Matthias Frey Date: Sun, 19 Mar 2023 18:41:32 +0000 Subject: [PATCH 10/12] matrix typed field --- unit_tests/Field/Field.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/unit_tests/Field/Field.cpp b/unit_tests/Field/Field.cpp index fd2d242d8..041660913 100644 --- a/unit_tests/Field/Field.cpp +++ b/unit_tests/Field/Field.cpp @@ -228,6 +228,8 @@ TEST_F(FieldTest, Hessian) { field_type field(*mesh, *layout); int nghost = field.getNghost(); auto view_field = field.getView(); + typedef ippl::Field, dim> matrix_field_type; + auto lDom = this->layout->getLocalNDIndex(); ippl::Vector hx = this->mesh->getMeshSpacing(); @@ -255,7 +257,7 @@ TEST_F(FieldTest, Hessian) { Kokkos::deep_copy(view_field, mirror); - vector_field_type result(*mesh, *layout); + matrix_field_type result(*mesh, *layout); result = hess(field); nghost = result.getNghost(); From e681bc24fdfbea125dda581fa7bd2a5f2215231f Mon Sep 17 00:00:00 2001 From: Matthias Frey Date: Sun, 19 Mar 2023 18:42:55 +0000 Subject: [PATCH 11/12] move typedef --- unit_tests/Field/Field.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/unit_tests/Field/Field.cpp b/unit_tests/Field/Field.cpp index 041660913..df8769a8d 100644 --- a/unit_tests/Field/Field.cpp +++ b/unit_tests/Field/Field.cpp @@ -224,12 +224,11 @@ TEST_F(FieldTest, Curl) { TEST_F(FieldTest, Hessian) { typedef ippl::Vector Vector_t; + typedef ippl::Field, dim> matrix_field_type; field_type field(*mesh, *layout); int nghost = field.getNghost(); auto view_field = field.getView(); - typedef ippl::Field, dim> matrix_field_type; - auto lDom = this->layout->getLocalNDIndex(); ippl::Vector hx = this->mesh->getMeshSpacing(); From 94bd21e1990b84a215202df703cd2da2b7caa890 Mon Sep 17 00:00:00 2001 From: Matthias Frey Date: Sun, 19 Mar 2023 18:47:38 +0000 Subject: [PATCH 12/12] add mesh and centering --- unit_tests/Field/Field.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/unit_tests/Field/Field.cpp b/unit_tests/Field/Field.cpp index df8769a8d..c63b18ffc 100644 --- a/unit_tests/Field/Field.cpp +++ b/unit_tests/Field/Field.cpp @@ -224,7 +224,7 @@ TEST_F(FieldTest, Curl) { TEST_F(FieldTest, Hessian) { typedef ippl::Vector Vector_t; - typedef ippl::Field, dim> matrix_field_type; + typedef ippl::Field, dim, Mesh_t, Centering_t> matrix_field_type; field_type field(*mesh, *layout); int nghost = field.getNghost();