diff --git a/CMakeLists.txt b/CMakeLists.txt index 91e072bf8..c87050f3c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,7 +1,7 @@ cmake_minimum_required (VERSION 3.1) project (IPPL CXX) -set (IPPL_VERSION_MAJOR 1) -set (IPPL_VERSION_MINOR 1.4) +set (IPPL_VERSION_MAJOR 3) +set (IPPL_VERSION_MINOR 0.1) set (IPPL_VERSION_NAME "V${IPPL_VERSION_MAJOR}.${IPPL_VERSION_MINOR}") @@ -37,8 +37,8 @@ set (CMAKE_CXX_FLAGS_RELWITHDEBINFO "-O3 -g ") set (CMAKE_CXX_FLAGS_RELEASE "-O3") set (CMAKE_CXX_FLAGS_DEBUG "-O0 -g") - -if (${CMAKE_CXX_COMPILER} MATCHES "GNU") +# Compile using C++20 for CPU builds +if (${CMAKE_CXX_COMPILER} MATCHES ".*mpicxx") set (CMAKE_CXX_STANDARD 20) endif () diff --git a/alpine/BumponTailInstability.cpp b/alpine/BumponTailInstability.cpp index c990ccfe1..f9e7428d6 100644 --- a/alpine/BumponTailInstability.cpp +++ b/alpine/BumponTailInstability.cpp @@ -7,7 +7,7 @@ // ny... = No. cell-centered points in the y-, z-, ...-direction // Np = Total no. of macro-particles in the simulation // Nt = Number of time steps -// stype = Field solver type e.g., FFT +// stype = Field solver type (FFT and CG supported) // lbthres = Load balancing threshold i.e., lbthres*100 is the maximum load imbalance // percentage which can be tolerated and beyond which // particle load balancing occurs. A value of 0.01 is good for many typical @@ -15,7 +15,7 @@ // ovfactor = Over-allocation factor for the buffers used in the communication. Typical // values are 1.0, 2.0. Value 1.0 means no over-allocation. // Example: -// srun ./BumponTailInstability 128 128 128 10000 10 FFT 0.01 2.0 --info 10 +// srun ./BumponTailInstability 128 128 128 10000 10 FFT 0.01 --overallocate 2.0 --info 10 // Change the TestName to TwoStreamInstability or BumponTailInstability // in order to simulate the Two stream instability or bump on tail instability // cases @@ -314,7 +314,7 @@ int main(int argc, char* argv[]) { hr[d] = rmax[d] / nr[d]; } Vector_t origin = rmin; - const double dt = 0.5 * hr[0]; // 0.05 + const double dt = std::min(.05, 0.5 * *std::min_element(hr.begin(), hr.end())); const bool isAllPeriodic = true; Mesh_t mesh(domain, hr, origin); @@ -322,8 +322,9 @@ int main(int argc, char* argv[]) { PLayout_t PL(FL, mesh); // Q = -\int\int f dx dv - double Q = std::reduce(rmax.begin(), rmax.end(), -1., std::multiplies()); - P = std::make_shared(PL, hr, rmin, rmax, decomp, Q); + double Q = std::reduce(rmax.begin(), rmax.end(), -1., std::multiplies()); + std::string solver = argv[arg++]; + P = std::make_shared(PL, hr, rmin, rmax, decomp, Q, solver); P->nr_m = nr; @@ -331,7 +332,6 @@ int main(int argc, char* argv[]) { bunch_type bunchBuffer(PL); - P->stype_m = argv[arg++]; P->initSolver(); P->time_m = 0.0; P->loadbalancethreshold_m = std::atof(argv[arg++]); @@ -433,13 +433,13 @@ int main(int argc, char* argv[]) { IpplTimings::startTimer(DummySolveTimer); P->rho_m = 0.0; - P->solver_mp->solve(); + P->runSolver(); IpplTimings::stopTimer(DummySolveTimer); P->scatterCIC(totalP, 0, hr); IpplTimings::startTimer(SolveTimer); - P->solver_mp->solve(); + P->runSolver(); IpplTimings::stopTimer(SolveTimer); P->gatherCIC(); @@ -489,7 +489,7 @@ int main(int argc, char* argv[]) { // Field solve IpplTimings::startTimer(SolveTimer); - P->solver_mp->solve(); + P->runSolver(); IpplTimings::stopTimer(SolveTimer); // gather E field diff --git a/alpine/ChargedParticles.hpp b/alpine/ChargedParticles.hpp index 8d02e6258..3ed092def 100644 --- a/alpine/ChargedParticles.hpp +++ b/alpine/ChargedParticles.hpp @@ -59,9 +59,16 @@ using VField_t = Field, Dim>; // heFFTe does not support 1D FFTs, so we switch to CG in the 1D case template -using Solver_t = std::conditional_t< - Dim == 1, ippl::ElectrostaticsCG, Centering_t>, - ippl::FFTPeriodicPoissonSolver, double, Dim, Mesh_t, Centering_t>>; +using CGSolver_t = ippl::ElectrostaticsCG, Centering_t>; + +template +using FFTSolver_t = + ippl::FFTPeriodicPoissonSolver, double, Dim, Mesh_t, Centering_t>; + +template +using Solver_t = + std::conditional_t, FFTSolver_t>, + std::variant>>; const double pi = std::acos(-1.0); @@ -161,12 +168,14 @@ class ChargedParticles : public ippl::ParticleBase { Vector_t rmin_m; Vector_t rmax_m; - double Q_m; - std::string stype_m; - std::shared_ptr> solver_mp; + double Q_m; + +private: + Solver_t solver_m; +public: double time_m; double rhoNorm_m; @@ -192,11 +201,12 @@ class ChargedParticles : public ippl::ParticleBase { } ChargedParticles(PLayout& pl, Vector_t hr, Vector_t rmin, Vector_t rmax, - ippl::e_dim_tag decomp[Dim], double Q) + ippl::e_dim_tag decomp[Dim], double Q, std::string solver) : ippl::ParticleBase(pl) , hr_m(hr) , rmin_m(rmin) , rmax_m(rmax) + , stype_m(solver) , Q_m(Q) { registerAttributes(); for (unsigned int i = 0; i < Dim; i++) { @@ -209,7 +219,7 @@ class ChargedParticles : public ippl::ParticleBase { void setPotentialBCs() { // CG requires explicit periodic boundary conditions while the periodic Poisson solver // simply assumes them - if constexpr (Dim == 1) { + if (stype_m == "CG") { for (unsigned int i = 0; i < 2 * Dim; ++i) { allPeriodic[i] = std::make_shared< ippl::PeriodicFace, Centering_t>>(i); @@ -235,7 +245,7 @@ class ChargedParticles : public ippl::ParticleBase { IpplTimings::startTimer(tupdateLayout); this->E_m.updateLayout(fl); this->rho_m.updateLayout(fl); - if constexpr (Dim == 1) { + if (stype_m == "CG") { this->phi_m.updateLayout(fl); phi_m.setFieldBC(allPeriodic); } @@ -255,7 +265,7 @@ class ChargedParticles : public ippl::ParticleBase { void initializeFields(Mesh_t& mesh, FieldLayout_t& fl) { E_m.initialize(mesh, fl); rho_m.initialize(mesh, fl); - if constexpr (Dim == 1) { + if (stype_m == "CG") { phi_m.initialize(mesh, fl); phi_m.setFieldBC(allPeriodic); } @@ -276,7 +286,11 @@ class ChargedParticles : public ippl::ParticleBase { } // Update this->updateLayout(fl, mesh, buffer, isFirstRepartition); - this->solver_mp->setRhs(rho_m); + if constexpr (Dim == 2 || Dim == 3) { + if (stype_m == "FFT") { + std::get>(solver_m).setRhs(rho_m); + } + } } bool balance(size_type totalP, const unsigned int nstep) { @@ -388,43 +402,83 @@ class ChargedParticles : public ippl::ParticleBase { } } - void initSolver(const ippl::ParameterList& sp) { - solver_mp = std::make_shared>(); + void runSolver() { + if (stype_m == "CG") { + CGSolver_t& solver = std::get>(solver_m); + solver.solve(); + + if (Ippl::Comm->rank() == 0) { + std::stringstream fname; + fname << "data/CG_"; + fname << Ippl::Comm->size(); + fname << ".csv"; + + Inform log(NULL, fname.str().c_str(), Inform::APPEND); + int iterations = solver.getIterationCount(); + // Assume the dummy solve is the first call + if (time_m == 0 && iterations == 0) { + log << "time,residue,iterations" << endl; + } + // Don't print the dummy solve + if (time_m > 0 || iterations > 0) { + log << time_m << "," << solver.getResidue() << "," << iterations << endl; + } + } + Ippl::Comm->barrier(); + } else if (stype_m == "FFT") { + if constexpr (Dim == 2 || Dim == 3) { + std::get>(solver_m).solve(); + } + } else { + throw std::runtime_error("Unknown solver type"); + } + } + + template + void initSolverWithParams(const ippl::ParameterList& sp) { + solver_m = Solver(); + Solver& solver = std::get(solver_m); - solver_mp->mergeParameters(sp); + solver.mergeParameters(sp); - solver_mp->setRhs(rho_m); + solver.setRhs(rho_m); - if constexpr (Dim == 1) { + if constexpr (std::is_same_v>) { // The CG solver computes the potential directly and // uses this to get the electric field - solver_mp->setLhs(phi_m); - solver_mp->setGradient(E_m); - } else { + solver.setLhs(phi_m); + solver.setGradient(E_m); + } else if constexpr (std::is_same_v>) { // The periodic Poisson solver computes the electric // field directly - solver_mp->setLhs(E_m); + solver.setLhs(E_m); } } void initCGSolver() { ippl::ParameterList sp; - sp.add("output_type", Solver_t::GRAD); + sp.add("output_type", CGSolver_t::GRAD); + // Increase tolerance in the 1D case + sp.add("tolerance", 1e-10); - initSolver(sp); + initSolverWithParams>(sp); } void initFFTSolver() { - ippl::ParameterList sp; - sp.add("output_type", Solver_t::GRAD); - sp.add("use_heffte_defaults", false); - sp.add("use_pencils", true); - sp.add("use_reorder", false); - sp.add("use_gpu_aware", true); - sp.add("comm", ippl::p2p_pl); - sp.add("r2c_direction", 0); - - initSolver(sp); + if constexpr (Dim == 2 || Dim == 3) { + ippl::ParameterList sp; + sp.add("output_type", FFTSolver_t::GRAD); + sp.add("use_heffte_defaults", false); + sp.add("use_pencils", true); + sp.add("use_reorder", false); + sp.add("use_gpu_aware", true); + sp.add("comm", ippl::p2p_pl); + sp.add("r2c_direction", 0); + + initSolverWithParams>(sp); + } else { + throw std::runtime_error("Unsupported dimensionality for FFT solver"); + } } void dumpData() { diff --git a/alpine/LandauDamping.cpp b/alpine/LandauDamping.cpp index 143976d10..46d12628e 100644 --- a/alpine/LandauDamping.cpp +++ b/alpine/LandauDamping.cpp @@ -7,7 +7,7 @@ // ny... = No. cell-centered points in the y-, z-, ...-direction // Np = Total no. of macro-particles in the simulation // Nt = Number of time steps -// stype = Field solver type e.g., FFT +// stype = Field solver type (FFT and CG supported) // lbthres = Load balancing threshold i.e., lbthres*100 is the maximum load imbalance // percentage which can be tolerated and beyond which // particle load balancing occurs. A value of 0.01 is good for many typical @@ -15,7 +15,7 @@ // ovfactor = Over-allocation factor for the buffers used in the communication. Typical // values are 1.0, 2.0. Value 1.0 means no over-allocation. // Example: -// srun ./LandauDamping 128 128 128 10000 10 FFT 0.01 2.0 --info 10 +// srun ./LandauDamping 128 128 128 10000 10 FFT 0.01 --overallocate 2.0 --info 10 // // Copyright (c) 2021, Sriramkrishnan Muralikrishnan, // Paul Scherrer Institut, Villigen PSI, Switzerland @@ -203,23 +203,22 @@ int main(int argc, char* argv[]) { // Q = -\int\int f dx dv double Q = std::reduce(rmax.begin(), rmax.end(), -1., std::multiplies()); Vector_t origin = rmin; - const double dt = 0.5 * hr[0]; + const double dt = std::min(.05, 0.5 * *std::min_element(hr.begin(), hr.end())); const bool isAllPeriodic = true; Mesh_t mesh(domain, hr, origin); FieldLayout_t FL(domain, decomp, isAllPeriodic); PLayout_t PL(FL, mesh); - P = std::make_unique(PL, hr, rmin, rmax, decomp, Q); + std::string solver = argv[arg++]; + P = std::make_unique(PL, hr, rmin, rmax, decomp, Q, solver); P->nr_m = nr; - P->E_m.initialize(mesh, FL); - P->rho_m.initialize(mesh, FL); + P->initializeFields(mesh, FL); bunch_type bunchBuffer(PL); - P->stype_m = argv[arg++]; P->initSolver(); P->time_m = 0.0; P->loadbalancethreshold_m = std::atof(argv[arg++]); @@ -303,13 +302,13 @@ int main(int argc, char* argv[]) { IpplTimings::startTimer(DummySolveTimer); P->rho_m = 0.0; - P->solver_mp->solve(); + P->runSolver(); IpplTimings::stopTimer(DummySolveTimer); P->scatterCIC(totalP, 0, hr); IpplTimings::startTimer(SolveTimer); - P->solver_mp->solve(); + P->runSolver(); IpplTimings::stopTimer(SolveTimer); P->gatherCIC(); @@ -360,7 +359,7 @@ int main(int argc, char* argv[]) { // Field solve IpplTimings::startTimer(SolveTimer); - P->solver_mp->solve(); + P->runSolver(); IpplTimings::stopTimer(SolveTimer); // gather E field diff --git a/alpine/PenningTrap.cpp b/alpine/PenningTrap.cpp index 3536abf89..e1e14da08 100644 --- a/alpine/PenningTrap.cpp +++ b/alpine/PenningTrap.cpp @@ -8,7 +8,7 @@ // nz = No. cell-centered points in the z-direction // Np = Total no. of macro-particles in the simulation // Nt = Number of time steps -// stype = Field solver type e.g., FFT +// stype = Field solver type (FFT and CG supported) // lbthres = Load balancing threshold i.e., lbthres*100 is the maximum load imbalance // percentage which can be tolerated and beyond which // particle load balancing occurs. A value of 0.01 is good for many typical @@ -16,7 +16,7 @@ // ovfactor = Over-allocation factor for the buffers used in the communication. Typical // values are 1.0, 2.0. Value 1.0 means no over-allocation. // Example: -// srun ./PenningTrap 128 128 128 10000 300 FFT 0.01 1.0 --info 10 +// srun ./PenningTrap 128 128 128 10000 300 FFT 0.01 --overallocate 1.0 --info 10 // // Copyright (c) 2021, Sriramkrishnan Muralikrishnan, // Paul Scherrer Institut, Villigen PSI, Switzerland @@ -151,6 +151,7 @@ double PDF(const Vector_t& xvec, const Vector_t& mu, const Vector_t FL(domain, decomp, isAllPeriodic); PLayout_t PL(FL, mesh); - double Q = -1562.5; - double Bext = 5.0; - P = std::make_unique(PL, hr, rmin, rmax, decomp, Q); + double Q = -1562.5; + double Bext = 5.0; + std::string solver = argv[6]; + P = std::make_unique(PL, hr, rmin, rmax, decomp, Q, solver); P->nr_m = nr; @@ -221,12 +223,10 @@ int main(int argc, char* argv[]) { sd[1] = 0.05 * length[1]; sd[2] = 0.20 * length[2]; - P->E_m.initialize(mesh, FL); - P->rho_m.initialize(mesh, FL); + P->initializeFields(mesh, FL); bunch_type bunchBuffer(PL); - P->stype_m = argv[6]; P->initSolver(); P->time_m = 0.0; P->loadbalancethreshold_m = std::atof(argv[7]); @@ -313,13 +313,13 @@ int main(int argc, char* argv[]) { IpplTimings::startTimer(DummySolveTimer); P->rho_m = 0.0; - P->solver_mp->solve(); + P->runSolver(); IpplTimings::stopTimer(DummySolveTimer); P->scatterCIC(totalP, 0, hr); IpplTimings::startTimer(SolveTimer); - P->solver_mp->solve(); + P->runSolver(); IpplTimings::stopTimer(SolveTimer); P->gatherCIC(); @@ -390,7 +390,7 @@ int main(int argc, char* argv[]) { // Field solve IpplTimings::startTimer(SolveTimer); - P->solver_mp->solve(); + P->runSolver(); IpplTimings::stopTimer(SolveTimer); // gather E field diff --git a/alpine/UniformPlasmaTest.cpp b/alpine/UniformPlasmaTest.cpp index 9dcf6d7eb..1fbbff961 100644 --- a/alpine/UniformPlasmaTest.cpp +++ b/alpine/UniformPlasmaTest.cpp @@ -7,13 +7,13 @@ // ny... = No. cell-centered points in the y-, z-, ...direction // Np = Total no. of macro-particles in the simulation // Nt = Number of time steps -// stype = Field solver type e.g., FFT +// stype = Field solver type (FFT and CG supported) // lbfreq = Load balancing frequency i.e., Number of time steps after which particle // load balancing should happen // ovfactor = Over-allocation factor for the buffers used in the communication. Typical // values are 1.0, 2.0. Value 1.0 means no over-allocation. // Example: -// srun ./UniformPlasmaTest 128 128 128 10000 10 FFT 10 1.0 --info 10 +// srun ./UniformPlasmaTest 128 128 128 10000 10 FFT 10 --overallocate 1.0 --info 10 // // Copyright (c) 2021, Sriramkrishnan Muralikrishnan, // Paul Scherrer Institut, Villigen PSI, Switzerland @@ -136,8 +136,9 @@ int main(int argc, char* argv[]) { FieldLayout_t FL(domain, decomp, isAllPeriodic); PLayout_t PL(FL, mesh); - double Q = -1562.5; - P = std::make_unique(PL, hr, rmin, rmax, decomp, Q); + double Q = -1562.5; + std::string solver = argv[arg++]; + P = std::make_unique(PL, hr, rmin, rmax, decomp, Q, solver); P->nr_m = nr; size_type nloc = totalP / Ippl::Comm->size(); @@ -166,8 +167,7 @@ int main(int argc, char* argv[]) { P->P = 0.0; IpplTimings::stopTimer(particleCreation); - P->E_m.initialize(mesh, FL); - P->rho_m.initialize(mesh, FL); + P->initializeFields(mesh, FL); bunch_type bunchBuffer(PL); @@ -177,14 +177,13 @@ int main(int argc, char* argv[]) { msg << "particles created and initial conditions assigned " << endl; - P->stype_m = argv[arg++]; P->initSolver(); P->time_m = 0.0; P->loadbalancefreq_m = std::atoi(argv[arg++]); IpplTimings::startTimer(DummySolveTimer); P->rho_m = 0.0; - P->solver_mp->solve(); + P->runSolver(); IpplTimings::stopTimer(DummySolveTimer); P->scatterCIC(totalP, 0, hr); @@ -192,7 +191,7 @@ int main(int argc, char* argv[]) { bool fromAnalyticDensity = false; IpplTimings::startTimer(SolveTimer); - P->solver_mp->solve(); + P->runSolver(); IpplTimings::stopTimer(SolveTimer); P->gatherCIC(); @@ -245,7 +244,7 @@ int main(int argc, char* argv[]) { // Field solve IpplTimings::startTimer(SolveTimer); - P->solver_mp->solve(); + P->runSolver(); IpplTimings::stopTimer(SolveTimer); // gather E field diff --git a/src/Field/BareField.h b/src/Field/BareField.h index c7413ef25..bc6c38eb4 100644 --- a/src/Field/BareField.h +++ b/src/Field/BareField.h @@ -72,15 +72,19 @@ namespace ippl { */ BareField(); + BareField(const BareField&) = default; + /*! Constructor for a BareField. The default constructor is deleted. * @param l of field * @param nghost number of ghost layers */ BareField(Layout_t& l, int nghost = 1); - BareField(const BareField&); - - BareField& operator=(const BareField&); + /*! + * Creates a new BareField with the same properties and contents + * @return A deep copy of the field + */ + BareField deepCopy() const; // Destroy the BareField. ~BareField() = default; @@ -203,9 +207,6 @@ namespace ippl { T min(int nghost = 0) const; T prod(int nghost = 0) const; - protected: - virtual void swap(BareField& other); - private: //! Number of ghost layers on each field boundary int nghost_m; diff --git a/src/Field/BareField.hpp b/src/Field/BareField.hpp index 0200b8f62..813643c12 100644 --- a/src/Field/BareField.hpp +++ b/src/Field/BareField.hpp @@ -38,28 +38,10 @@ namespace ippl { , layout_m(nullptr) {} template - BareField::BareField(const BareField& other) - : nghost_m(other.nghost_m) - , owned_m(other.owned_m) - , layout_m(other.layout_m) { - setup(); - Kokkos::deep_copy(dview_m, other.dview_m); - } - - template - BareField& BareField::operator=(const BareField& other) { - BareField copy(other); - swap(copy); - - return *this; - } - - template - void BareField::swap(BareField& other) { - nghost_m = other.nghost_m; - owned_m = other.owned_m; - layout_m = other.layout_m; - std::swap(dview_m, other.dview_m); + BareField BareField::deepCopy() const { + BareField copy(*layout_m, nghost_m); + Kokkos::deep_copy(copy.dview_m, dview_m); + return copy; } template diff --git a/src/Field/Field.h b/src/Field/Field.h index f1a5e251e..594c818cb 100644 --- a/src/Field/Field.h +++ b/src/Field/Field.h @@ -44,9 +44,13 @@ namespace ippl { // been properly initialized. Field(); - Field(const Field&); + Field(const Field&) = default; - Field& operator=(const Field&); + /*! + * Creates a new Field with the same properties and contents + * @return A deep copy of the field + */ + Field deepCopy() const; virtual ~Field() = default; @@ -83,9 +87,6 @@ namespace ippl { // Assignment from constants and other arrays. using BareField::operator=; - protected: - virtual void swap(Field& other); - private: // The Mesh object, and a flag indicating if we constructed it Mesh_t* mesh_m; diff --git a/src/Field/Field.hpp b/src/Field/Field.hpp index 0247b457c..8ae3815f9 100644 --- a/src/Field/Field.hpp +++ b/src/Field/Field.hpp @@ -35,25 +35,11 @@ namespace ippl { , bc_m() {} template - Field::Field(const Field& other) - : BareField(other) - , mesh_m(other.mesh_m) - , bc_m(other.bc_m) {} + Field Field::deepCopy() const { + Field copy(*mesh_m, this->getLayout(), this->getNghost()); + Kokkos::deep_copy(copy.getView(), this->getView()); - template - Field& Field::operator=( - const Field& other) { - Field copy(other); - swap(copy); - - return *this; - } - - template - void Field::swap(Field& other) { - BareField::swap(other); - mesh_m = other.mesh_m; - std::swap(bc_m, other.bc_m); + return copy; } ////////////////////////////////////////////////////////////////////////// diff --git a/src/Solver/ElectrostaticsCG.h b/src/Solver/ElectrostaticsCG.h index 797ebddf4..9177ab380 100644 --- a/src/Solver/ElectrostaticsCG.h +++ b/src/Solver/ElectrostaticsCG.h @@ -57,7 +57,7 @@ namespace ippl { int output = this->params_m.template get("output_type"); if (output & Base::GRAD) { - *(this->grad_mp) = grad(*(this->lhs_mp)); + *(this->grad_mp) = -grad(*(this->lhs_mp)); } } @@ -68,6 +68,12 @@ namespace ippl { */ int getIterationCount() { return algo_m.getIterationCount(); } + /*! + * Query the residue + * @return Residue norm from last solve + */ + Tlhs getResidue() const { return algo_m.getResidue(); } + protected: algo algo_m = algo(); diff --git a/src/Solver/PCG.h b/src/Solver/PCG.h index abafbd27e..5945408ac 100644 --- a/src/Solver/PCG.h +++ b/src/Solver/PCG.h @@ -26,8 +26,10 @@ namespace ippl { template class PCG : public SolverAlgorithm { - public: using Base = SolverAlgorithm; + typedef typename Base::lhs_type::type T; + + public: using typename Base::lhs_type; using typename Base::rhs_type; using operator_type = std::function; @@ -46,8 +48,6 @@ namespace ippl { int getIterationCount() { return iterations_m; } void operator()(lhs_type& lhs, rhs_type& rhs, const ParameterList& params) override { - typedef typename lhs_type::type T; - typename lhs_type::Mesh_t mesh = lhs.get_mesh(); typename lhs_type::Layout_t layout = lhs.getLayout(); @@ -83,16 +83,16 @@ namespace ippl { r = rhs - op_m(lhs); - lhs_type d(r); + lhs_type d = r.deepCopy(); d.setFieldBC(bc); T delta1 = innerProduct(r, r); - T rNorm = std::sqrt(delta1); + residueNorm = std::sqrt(delta1); const T tolerance = params.get("tolerance") * norm(rhs); lhs_type q(mesh, layout); - while (iterations_m < maxIterations && rNorm > tolerance) { + while (iterations_m < maxIterations && residueNorm > tolerance) { q = op_m(d); T alpha = delta1 / innerProduct(d, q); lhs = lhs + alpha * d; @@ -110,7 +110,7 @@ namespace ippl { delta1 = innerProduct(r, r); T beta = delta1 / delta0; - rNorm = std::sqrt(delta1); + residueNorm = std::sqrt(delta1); d = r + beta * d; @@ -123,8 +123,11 @@ namespace ippl { } } + T getResidue() const { return residueNorm; } + protected: operator_type op_m; + T residueNorm = 0; int iterations_m = 0; }; diff --git a/unit_tests/BareField/BareField.cpp b/unit_tests/BareField/BareField.cpp index 9b5cc4822..b3257915e 100644 --- a/unit_tests/BareField/BareField.cpp +++ b/unit_tests/BareField/BareField.cpp @@ -77,9 +77,11 @@ struct FieldVal { } }; -TEST_F(BareFieldTest, CopyConstruction) { +TEST_F(BareFieldTest, DeepCopy) { auto check = [](std::shared_ptr>& field) { - field_type copy(*field); + *field = 0; + field_type copy = field->deepCopy(); + copy = copy + 1; auto mirrorA = field->getHostMirror(); auto mirrorB = copy.getHostMirror(); @@ -87,27 +89,8 @@ TEST_F(BareFieldTest, CopyConstruction) { Kokkos::deep_copy(mirrorA, field->getView()); Kokkos::deep_copy(mirrorB, copy.getView()); - nestedViewLoop(mirrorA, 0, [&](const Idx... args) { - ASSERT_DOUBLE_EQ(mirrorA(args...), mirrorB(args...)); - }); - }; - - apply(check, fields); -} - -TEST_F(BareFieldTest, CopyAssignment) { - auto check = [](std::shared_ptr>& field) { - field_type copy; - copy = *field; - - auto mirrorA = field->getHostMirror(); - auto mirrorB = copy.getHostMirror(); - - Kokkos::deep_copy(mirrorA, field->getView()); - Kokkos::deep_copy(mirrorB, copy.getView()); - - nestedViewLoop(mirrorA, 0, [&](const Idx... args) { - ASSERT_DOUBLE_EQ(mirrorA(args...), mirrorB(args...)); + nestedViewLoop(mirrorA, field->getNghost(), [&](const Idx... args) { + ASSERT_DOUBLE_EQ(mirrorA(args...) + 1, mirrorB(args...)); }); }; diff --git a/unit_tests/Field/Field.cpp b/unit_tests/Field/Field.cpp index 5dfb6b755..b593dab63 100644 --- a/unit_tests/Field/Field.cpp +++ b/unit_tests/Field/Field.cpp @@ -158,9 +158,11 @@ struct FieldVal { } }; -TEST_F(FieldTest, CopyConstruction) { +TEST_F(FieldTest, DeepCopy) { auto check = [](std::shared_ptr>& field) { - field_type copy(*field); + *field = 0; + field_type copy = field->deepCopy(); + copy = copy + 1.; auto mirrorA = field->getHostMirror(); auto mirrorB = copy.getHostMirror(); @@ -168,27 +170,8 @@ TEST_F(FieldTest, CopyConstruction) { Kokkos::deep_copy(mirrorA, field->getView()); Kokkos::deep_copy(mirrorB, copy.getView()); - nestedViewLoop(mirrorA, 0, [&](const Idx... args) { - ASSERT_DOUBLE_EQ(mirrorA(args...), mirrorB(args...)); - }); - }; - - apply(check, fields); -} - -TEST_F(FieldTest, CopyAssignment) { - auto check = [](std::shared_ptr>& field) { - field_type copy; - copy = *field; - - auto mirrorA = field->getHostMirror(); - auto mirrorB = copy.getHostMirror(); - - Kokkos::deep_copy(mirrorA, field->getView()); - Kokkos::deep_copy(mirrorB, copy.getView()); - - nestedViewLoop(mirrorA, 0, [&](const Idx... args) { - ASSERT_DOUBLE_EQ(mirrorA(args...), mirrorB(args...)); + nestedViewLoop(mirrorA, field->getNghost(), [&](const Idx... args) { + ASSERT_DOUBLE_EQ(mirrorA(args...) + 1, mirrorB(args...)); }); };