diff --git a/.gitignore b/.gitignore index 8cbc70e15..7a3725228 100644 --- a/.gitignore +++ b/.gitignore @@ -1,25 +1,27 @@ +# Build directories ippl-IPPL-1.0.0 build/ -doc/html -doc/latex -manual/*.aux -manual/*.log -manual/*.toc -manual/*.idx -manual/*.lof -manual/*.lot -manual/ippl_user_guide.pdf -build*/ -*~ -*.~ -*.orig +build_*/ +notebooks/ +cmake-build*/ +CMakeFiles/ +CMakeCache.txt +*.cmake +*.sh +vic.* +*.csv +# Spack +.spack-env/ -# ignore build directories -build*/ +# Python virtual environments +.venv/ +venv/ -# ignore .vscode directory -.vscode -.cache +# Python cache +__pycache__/ +*.pyc -# mac +# Editor / OS junk +*~ +*.swp .DS_Store diff --git a/CMakeLists.txt b/CMakeLists.txt index 54c104885..869de6d03 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -38,6 +38,7 @@ option(IPPL_ENABLE_UNIT_TESTS "Enable unit tests using GoogleTest" OFF) option(IPPL_ENABLE_FFT "Enable FFT support" OFF) option(IPPL_ENABLE_SOLVERS "Enable IPPL solvers" OFF) option(IPPL_ENABLE_ALPINE "Enable building the Alpine module" OFF) +option(IPPL_ENABLE_ALVINE "Enable building the Alvine module",OFF) option(IPPL_ENABLE_COSMOLOGY "Enable building the Cosmology module" OFF) option(IPPL_ENABLE_EXAMPLES "Enable building the Example module" OFF) option(IPPL_ENABLE_TESTS "Build integration tests in test/ directory" OFF) @@ -103,7 +104,7 @@ endif() # Define sources for project # ------------------------------------------------------------------------------ add_subdirectory(src) - +add_subdirectory(alvine) # ------------------------------------------------------------------------------ # Include testing related material (BUILD_TESTING is defined in CTest module) # ------------------------------------------------------------------------------ @@ -127,6 +128,10 @@ if(IPPL_ENABLE_ALPINE) add_subdirectory(alpine) endif() +if(IPPL_ENABLE_ALVINE) + add_subdirectory(alvine) +endif() + if(IPPL_ENABLE_COSMOLOGY) add_subdirectory(cosmology) endif() diff --git a/DartConfiguration.tcl b/DartConfiguration.tcl new file mode 100644 index 000000000..76997f069 --- /dev/null +++ b/DartConfiguration.tcl @@ -0,0 +1,106 @@ +# This file is configured by CMake automatically as DartConfiguration.tcl +# If you choose not to use CMake, this file may be hand configured, by +# filling in the required variables. + + +# Configuration directories and files +SourceDirectory: /home/raven/Documents/vic/my_implementation/ippl +BuildDirectory: /home/raven/Documents/vic/my_implementation/ippl + +# Where to place the cost data store +CostDataFile: + +# Site is something like machine.domain, i.e. pragmatic.crd +Site: raven-VivoBook-ASUSLaptop-X515JA + +# Build name is osname-revision-compiler, i.e. Linux-2.4.2-2smp-c++ +BuildName: Linux-c++ + +# Subprojects +LabelsForSubprojects: + +# Submission information +SubmitURL: https://my.cdash.org/submit.php?project=IPPL +SubmitInactivityTimeout: + +# Dashboard start time +NightlyStartTime: 00:00:00 + +# Commands for the build/test/submit cycle +ConfigureCommand: "/usr/bin/cmake" "/home/raven/Documents/vic/my_implementation/ippl" +MakeCommand: /usr/bin/cmake --build . --config "${CTEST_CONFIGURATION_TYPE}" +DefaultCTestConfigurationType: Release + +# version control +UpdateVersionOnly: + +# CVS options +# Default is "-d -P -A" +CVSCommand: +CVSUpdateOptions: + +# Subversion options +SVNCommand: +SVNOptions: +SVNUpdateOptions: + +# Git options +GITCommand: /usr/bin/git +GITInitSubmodules: +GITUpdateOptions: +GITUpdateCustom: + +# Perforce options +P4Command: +P4Client: +P4Options: +P4UpdateOptions: +P4UpdateCustom: + +# Generic update command +UpdateCommand: /usr/bin/git +UpdateOptions: +UpdateType: git + +# Compiler info +Compiler: /usr/bin/c++ +CompilerVersion: 13.3.0 + +# Dynamic analysis (MemCheck) +PurifyCommand: +ValgrindCommand: +ValgrindCommandOptions: +DrMemoryCommand: +DrMemoryCommandOptions: +CudaSanitizerCommand: +CudaSanitizerCommandOptions: +MemoryCheckType: +MemoryCheckSanitizerOptions: +MemoryCheckCommand: MEMORYCHECK_COMMAND-NOTFOUND +MemoryCheckCommandOptions: +MemoryCheckSuppressionFile: + +# Coverage +CoverageCommand: /usr/bin/gcov +CoverageExtraFlags: -l + +# Testing options +# TimeOut is the amount of time in seconds to wait for processes +# to complete during testing. After TimeOut seconds, the +# process will be summarily terminated. +# Currently set to 25 minutes +TimeOut: 1500 + +# During parallel testing CTest will not start a new test if doing +# so would cause the system load to exceed this value. +TestLoad: + +UseLaunchers: +CurlOptions: +# warning, if you add new options here that have to do with submit, +# you have to update cmCTestSubmitCommand.cxx + +# For CTest submissions that timeout, these options +# specify behavior for retrying the submission +CTestSubmitRetryDelay: 5 +CTestSubmitRetryCount: 3 diff --git a/alvine/AlvineManager.h b/alvine/AlvineManager.h new file mode 100644 index 000000000..6b05f80fe --- /dev/null +++ b/alvine/AlvineManager.h @@ -0,0 +1,311 @@ +#ifndef IPPL_ALVINE_MANAGER_H +#define IPPL_ALVINE_MANAGER_H + +#include + +#include "FieldContainer.hpp" +#include "FieldSolver.hpp" +#include "LoadBalancer.hpp" +#include "Manager/BaseManager.h" +#include "Manager/PicManager.h" +#include "ParticleContainer.hpp" +#include "Random/Distribution.h" +#include "Random/InverseTransformSampling.h" +#include "Random/NormalDistribution.h" +#include "Random/Randn.h" + +using view_type = typename ippl::detail::ViewType, 1>::view_type; + +template +class AlvineManager + : public ippl::PicManager, FieldContainer, + LoadBalancer> { +public: + using ParticleContainer_t = ParticleContainer; + using FieldContainer_t = FieldContainer; + using FieldSolver_t= FieldSolver; + using LoadBalancer_t= LoadBalancer; + using Base= ippl::ParticleBase>; + +protected: + unsigned nt_m; + unsigned it_m; + Vector_t nr_m; + unsigned np_m; + std::array decomp_m; + bool isAllPeriodic_m; + ippl::NDIndex domain_m; + std::string solver_m; + int dump_freq_m; + +public: + AlvineManager(unsigned nt_, Vector_t& nr_, unsigned np_, std::string& solver_, int dump_freq_) + : ippl::PicManager, FieldContainer, LoadBalancer>() + , nt_m(nt_) + , nr_m(nr_) + , np_m(np_) + , solver_m(solver_) + , dump_freq_m(dump_freq_) {} + + ~AlvineManager(){} + +protected: + double time_m; + double dt_m; + Vector_t rmin_m; + Vector_t rmax_m; + Vector_t origin_m; + Vector_t hr_m; + double energy0_m = 0.0; + bool energy_initialized_m = false; + double enstrophy0_m = 0.0; + bool enstrophy_initialized_m = false; +public: + + double getTime() { return time_m; } + + void setTime(double time_) { time_m = time_; } + + int getNt() const { return nt_m; } + + void setNt(int nt_) { nt_m = nt_; } + + virtual void dump() { /* default does nothing */ }; + + void pre_step() override { + } + + void post_step() override { + Inform m("Step: "); + this->time_m += this->dt_m; + this->it_m++; + + if(this->it_m % dump_freq_m == 0) { + this->dump(); + } + m << this->it_m << " Done" << endl; + } + + void grid2par() override { + gatherCIC(); + } + + void gatherCIC() { + this->pcontainer_m->P = 0.0; + gather(this->pcontainer_m->P, this->fcontainer_m->getUField(), this->pcontainer_m->R); + } + + void par2grid() override { + scatterCIC(); + } + + void computeVelocityField() { + + VField_t u_field = this->fcontainer_m->getUField(); + u_field = 0.0; + + if constexpr (Dim == 2) { + const int nghost = u_field.getNghost(); + auto view = u_field.getView(); + + auto omega_view = this->fcontainer_m->getOmegaField().getView(); + this->fcontainer_m->getOmegaField().fillHalo(); + + Vector_t hr = hr_m; + Kokkos::parallel_for( + "Assign rhs", ippl::getRangePolicy(view, nghost), + KOKKOS_LAMBDA(const int i, const int j) { + view(i, j) = { + (omega_view(i, j + 1) - omega_view(i, j - 1)) / (2 * hr(1)), + -(omega_view(i + 1, j) - omega_view(i - 1, j)) / (2 * hr(0)) + }; + + }); + } else if constexpr (Dim == 3) { + //TODO compute velocity field in 3D, this should be a simple curl operation (one line) + } + } + +double relativeError(double value, double reference) const { + return std::fabs((value - reference) / std::max(std::fabs(reference), 1e-30)); +} + +double computeParticleCirculation() { + double gamma_local = 0.0; + + auto omega_view = this->pcontainer_m->omega.getView(); + auto nlocal = this->pcontainer_m->getLocalNum(); + + Kokkos::parallel_reduce( + "particle_circulation", + nlocal, + KOKKOS_LAMBDA(const int i, double& lsum) { + lsum += omega_view(i); + }, + gamma_local + ); + + double gamma_global = 0.0; + ippl::Comm->reduce(gamma_local, gamma_global, 1, std::plus()); + + return gamma_global; +} + + +double computeGridCirculation() { + double gamma_local = 0.0; + + auto& omegaField = this->fcontainer_m->getOmegaField(); + auto omega_view = omegaField.getView(); + + const double dA = hr_m[0] * hr_m[1]; + const int nghost = omegaField.getNghost(); + + Kokkos::parallel_reduce( + "grid_circulation", + ippl::getRangePolicy(omega_view, nghost), + KOKKOS_LAMBDA(const int i, const int j, double& lsum) { + lsum += omega_view(i, j); + }, + gamma_local + ); + + gamma_local *= dA; + + double gamma_global = 0.0; + ippl::Comm->reduce(gamma_local, gamma_global, 1, std::plus()); + + return gamma_global; +} + +void checkCirculationConservation(double relError, Inform& m) { + size_type TotalParticles = 0; + size_type localParticles = this->pcontainer_m->getLocalNum(); + + ippl::Comm->reduce(localParticles, TotalParticles, 1, std::plus()); + + if (ippl::Comm->rank() == 0) { + if (TotalParticles != np_m || relError > 1e-12) { + m << "Time step: " << it_m << endl; + m << "Total particles expected: " << np_m + << " after update: " << TotalParticles << endl; + m << "Rel. error in circulation conservation: " << relError << endl; + ippl::Comm->abort(); + } + } +} + + +double computeKineticEnergy() { + double energy_local = 0.0; + + auto& uField = this->fcontainer_m->getUField(); + auto u_view = uField.getView(); + + const double dA = hr_m[0] * hr_m[1]; + const int nghost = uField.getNghost(); + + Kokkos::parallel_reduce( + "kinetic_energy", + ippl::getRangePolicy(u_view, nghost), + KOKKOS_LAMBDA(const int i, const int j, double& lsum) { + const double ux = u_view(i, j)[0]; + const double uy = u_view(i, j)[1]; + lsum += 0.5 * (ux * ux + uy * uy); + }, + energy_local + ); + + energy_local *= dA; + + double energy_global = 0.0; + ippl::Comm->reduce(energy_local, energy_global, 1, std::plus()); + + return energy_global; +} + +void checkEnergyConservation(double energy, double relError, Inform& m) { + if (ippl::Comm->rank() == 0) { + m << "kinetic energy = " << energy + << ", relError = " << relError << endl; + } +} + +double computeEnstrophy() { + double enstrophy_local = 0.0; + + auto& omegaField = this->fcontainer_m->getOmegaField(); + auto omega_view = omegaField.getView(); + const double dA = hr_m[0] * hr_m[1]; + + const int nghost = omegaField.getNghost(); + + Kokkos::parallel_reduce( + "enstrophy", + ippl::getRangePolicy(omega_view, nghost), + KOKKOS_LAMBDA(const int i, const int j, double& lsum) { + const double omega = omega_view(i, j); + lsum += 0.5 * omega * omega; + }, + enstrophy_local + ); + + enstrophy_local *= dA; + + double enstrophy_global = 0.0; + ippl::Comm->reduce(enstrophy_local, enstrophy_global, 1, std::plus()); + + return enstrophy_global; +} + + + +double computeDivergenceL2() { + auto& uField = this->fcontainer_m->getUField(); +// uField.fillHalo(); + + auto divField = this->fcontainer_m->getOmegaField().deepCopy(); + + divField = div(uField); + double N = this->nr_m[0]*this->nr_m[1]; + double div_l2 = norm(divField, 2)/std::sqrt(N); + + // restore omega by recomputing par2grid later if needed + return div_l2; +} + + +void scatterCIC() { + Inform m("scatter "); + + this->fcontainer_m->getOmegaField() = 0.0; + + if constexpr (Dim == 2) { + // Scatter particle strengths to grid + scatter(this->pcontainer_m->omega, + this->fcontainer_m->getOmegaField(), + this->pcontainer_m->R); + + // Convert deposited circulation to vorticity density + this->fcontainer_m->getOmegaField() = + this->fcontainer_m->getOmegaField() / (hr_m[0] * hr_m[1]); + + // Conservation check + double gammaParticles = computeParticleCirculation(); + double gammaGrid = computeGridCirculation(); + + double relError = std::fabs((gammaParticles - gammaGrid) / + std::max(std::fabs(gammaParticles), 1e-30)); + + m << "particle circulation = " << gammaParticles + << ", grid circulation = " << gammaGrid + << ", relError = " << relError << endl; + + checkCirculationConservation(relError, m); + + } else if constexpr (Dim == 3) { + // TODO 3D version + } +} +}; +#endif diff --git a/alvine/CMakeLists.txt b/alvine/CMakeLists.txt new file mode 100644 index 000000000..261337ce9 --- /dev/null +++ b/alvine/CMakeLists.txt @@ -0,0 +1,33 @@ +file (RELATIVE_PATH _relPath "${CMAKE_SOURCE_DIR}" "${CMAKE_CURRENT_SOURCE_DIR}") +message (STATUS "Adding index test found in ${_relPath}") + +include_directories ( + ${CMAKE_SOURCE_DIR}/src +) + +link_directories ( + ${CMAKE_CURRENT_SOURCE_DIR} + ${Kokkos_DIR}/.. +) + +find_package(Heffte CONFIG REQUIRED) +set (IPPL_LIBS ippl ${MPI_CXX_LIBRARIES}) +set (COMPILE_FLAGS ${OPAL_CXX_FLAGS}) + +add_executable (VortexInCell VortexInCell.cpp) +add_executable(Fsl Fsl.cpp) +target_include_directories(VortexInCell PRIVATE ${CMAKE_SOURCE_DIR}/src ${CMAKE_SOURCE_DIR}/alvine ) +target_link_libraries(VortexInCell ${IPPL_LIBS} ippl MPI::MPI_CXX Heffte::Heffte) + +target_include_directories(Fsl PRIVATE ${CMAKE_SOURCE_DIR}/src ${CMAKE_SOURCE_DIR}/alvine ) +target_link_libraries(Fsl ${IPPL_LIBS} ippl MPI::MPI_CXX Heffte::Heffte) + + +# vi: set et ts=4 sw=4 sts=4: + +# Local Variables: +# mode: cmake +# cmake-tab-width: 4 +# indent-tabs-mode: nil +# require-final-newline: nil +# End: diff --git a/alvine/FieldContainer.hpp b/alvine/FieldContainer.hpp new file mode 100644 index 000000000..9272b562c --- /dev/null +++ b/alvine/FieldContainer.hpp @@ -0,0 +1,77 @@ +#ifndef IPPL_FIELD_CONTAINER_H +#define IPPL_FIELD_CONTAINER_H + +#include + +#include "Manager/BaseManager.h" + +// Define the FieldsContainer class +template +class FieldContainer{ + using vorticity_field_type = std::conditional, VField_t>::type; + + +public: + FieldContainer(Vector_t& hr, Vector_t& rmin, + Vector_t& rmax, std::array decomp, + ippl::NDIndex domain, Vector_t origin, + bool isAllPeriodic) + : hr_m(hr) + , rmin_m(rmin) + , rmax_m(rmax) + , decomp_m(decomp) + , mesh_m(domain, hr, origin) + , fl_m(MPI_COMM_WORLD, domain, decomp, isAllPeriodic) {} + + ~FieldContainer(){} + +private: + Vector_t hr_m; + Vector_t rmin_m; + Vector_t rmax_m; + std::array decomp_m; + + vorticity_field_type A_field_m; + vorticity_field_type omega_field_m; + VField_t u_field_m; + + Mesh_t mesh_m; + FieldLayout_t fl_m; + +public: + + vorticity_field_type& getA_field() { return A_field_m; } + void setA_field(vorticity_field_type& A_field) { A_field_m = A_field; } + + vorticity_field_type& getOmegaField() { return omega_field_m; } + void setOmega_field(vorticity_field_type& omega_field) { omega_field_m = omega_field; } + + VField_t& getUField() { return u_field_m; } + void setUField(VField_t& u_field) { u_field_m = u_field; } + + Vector_t& getHr() { return hr_m; } + void setHr(const Vector_t& hr) { hr_m = hr; } + + Vector_t& getRMin() { return rmin_m; } + void setRMin(const Vector_t& rmin) { rmin_m = rmin; } + + Vector_t& getRMax() { return rmax_m; } + void setRMax(const Vector_t& rmax) { rmax_m = rmax; } + + std::array getDecomp() { return decomp_m; } + void setDecomp(std::array decomp) { decomp_m = decomp; } + + Mesh_t& getMesh() { return mesh_m; } + void setMesh(Mesh_t & mesh) { mesh_m = mesh; } + + FieldLayout_t& getFL() { return fl_m; } + void setFL(std::shared_ptr>& fl) { fl_m = fl; } + + void initializeFields() { + A_field_m.initialize(mesh_m, fl_m); + omega_field_m.initialize(mesh_m, fl_m); + u_field_m.initialize(mesh_m, fl_m); + } +}; + +#endif diff --git a/alvine/FieldSolver.hpp b/alvine/FieldSolver.hpp new file mode 100644 index 000000000..0a75e5692 --- /dev/null +++ b/alvine/FieldSolver.hpp @@ -0,0 +1,81 @@ +#ifndef IPPL_FIELD_SOLVER_H +#define IPPL_FIELD_SOLVER_H + +#include +#include "Manager/BaseManager.h" +#include "Manager/FieldSolverBase.h" + +// Define the FieldSolver class +template +class FieldSolver : public ippl::FieldSolverBase { + using omega_field_type = std::conditional, VField_t>::type; + private: + omega_field_type *omega_m; + + public: + FieldSolver(std::string solver, omega_field_type *omega) + : ippl::FieldSolverBase(solver) + , omega_m(omega) {} + + ~FieldSolver(){} + + Field_t *getOmega() const { return omega_m; } + void setOmega(Field_t *omega){ omega_m = omega; } + + void initSolver() override { + + Inform m("solver "); + if (this->getStype() == "FFT") { + initFFTSolver(); + } else { + m << "No solver matches the argument" << endl; + } + } + + void runSolver() override { + if constexpr (Dim == 2) { + if (this->getStype() == "FFT") { + std::get>(this->getSolver()).solve(); + } else { + throw std::runtime_error("Unknown solver type"); + } + } else if constexpr (Dim == 3) { + //TODO: probably need to run solver three times in this case and do ugly conversion. + + } + } + + template + void initSolverWithParams(const ippl::ParameterList& sp) { + this->getSolver().template emplace(); + + Solver& solver = std::get(this->getSolver()); + + solver.mergeParameters(sp); + + if constexpr (Dim == 2) { + solver.setRhs(*omega_m); + } else if constexpr (Dim == 3) { + //TODO probably need ugly construct for 3D case here. + } + } + + void initFFTSolver() { + if constexpr (Dim == 2 || Dim == 3) { + ippl::ParameterList sp; + sp.add("output_type", FFTSolver_t::SOL); + 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"); + } + } + +}; +#endif diff --git a/alvine/Fsl.cpp b/alvine/Fsl.cpp new file mode 100644 index 000000000..d79c6c785 --- /dev/null +++ b/alvine/Fsl.cpp @@ -0,0 +1,96 @@ +// Forward Semi-Lagrangian Test +// Usage: +// srun ./FSL --overallocate 1.0 --info 5 + +constexpr unsigned Dim = 2; +using T = double; +const char* TestName = "FSL"; + +#include "Ippl.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "datatypes.h" + +#include "Utility/IpplTimings.h" + +#include "Manager/PicManager.h" +#include "FslManager.h" +#include "VortexDistributions.h" + +int main(int argc, char* argv[]) { + ippl::initialize(argc, argv); + { + Inform msg(TestName); + + static IpplTimings::TimerRef mainTimer = IpplTimings::getTimer("total"); + IpplTimings::startTimer(mainTimer); + + unsigned arg = 1; + + Vector_t nr; + for (unsigned d = 0; d < Dim; d++) { + nr[d] = std::atoi(argv[arg++]); + } + + int np = std::atoi(argv[arg++]); // unused for FSL, but keep argument format + int nt = std::atoi(argv[arg++]); + + std::string solver = argv[arg++]; + int dump_freq = std::atoi(argv[arg++]); + + msg << " Grid size: " << nr + << " No. of virtual particles per step: " << nr[0] * nr[1] + << " No. of time steps: " << nt << endl; + + ippl::NDIndex domain; + for (unsigned i = 0; i < Dim; i++) { + domain[i] = ippl::Index(nr[i]); + } + + Vector_t rmin(0.0); + Vector_t rmax(10.0); + Vector_t origin = rmin; + Vector_t hr = (rmax - rmin) / nr; + + std::array isParallel; + isParallel.fill(true); + + const bool isAllPeriodic = true; + + Mesh_t mesh(domain, hr, origin); + FieldLayout_t FL(MPI_COMM_WORLD, domain, isParallel, isAllPeriodic); + + FSLManager manager( + nt, + nr, + np, + solver, + dump_freq, + rmin, + rmax, + origin, + FL, + mesh + ); + + manager.pre_run(); + manager.run(manager.getNt()); + + IpplTimings::stopTimer(mainTimer); + IpplTimings::print(); + IpplTimings::print(std::string("timing_fsl.dat")); + } + + ippl::finalize(); + + return 0; +} diff --git a/alvine/FslManager.h b/alvine/FslManager.h new file mode 100644 index 000000000..776d1802b --- /dev/null +++ b/alvine/FslManager.h @@ -0,0 +1,556 @@ +#ifndef IPPL_VORTEX_IN_CELL_MANAGER_H +#define IPPL_VORTEX_IN_CELL_MANAGER_H + +#include + +#include "AlvineManager.h" +#include "FieldContainer.hpp" +#include "FieldSolver.hpp" +#include "LoadBalancer.hpp" +#include "Manager/BaseManager.h" +#include "ParticleContainer.hpp" +#include "Random/Distribution.h" +#include "Random/InverseTransformSampling.h" +#include "Random/NormalDistribution.h" +#include "Random/Randu.h" +#include "SinusoidalJitter.hpp" +#include "VortexDistributions.h" +#include "VtkDump.hpp" + +using view_type = typename ippl::detail::ViewType, 1>::view_type; +using host_type = typename ippl::ParticleAttrib::host_mirror_type; /*using host_type = typename ippl::ParticleAttrib::HostMirror;*/ + + +template +class FSLManager : public AlvineManager { +public: + using ParticleContainer_t = ParticleContainer; + using FieldContainer_t = FieldContainer; + using FieldSolver_t = FieldSolver; + using LoadBalancer_t = LoadBalancer; + FieldLayout_t FL_m; // Store the field layout + Mesh_t mesh_m; // Store the mesh + + // Constructor declaration + FSLManager(unsigned nt_, Vector_t& nr_, unsigned np_, + std::string& solver_, int dump_freq_, + Vector_t rmin_ = 0.0, + Vector_t rmax_ = 10.0, + Vector_t origin_ = 0.0, + FieldLayout_t& FL_ = nullptr, + Mesh_t& mesh_ = nullptr) + : AlvineManager(nt_, nr_, np_, solver_, dump_freq_) { + this->rmin_m = rmin_; + this->rmax_m = rmax_; + this->origin_m = origin_; + this->FL_m = FL_; // Store the layout + this->mesh_m = mesh_; // Store the mesh + } + + ~FSLManager() {} + +void pre_run() override { + for (unsigned i = 0; i < Dim; i++) { + this->domain_m[i] = ippl::Index(this->nr_m[i]); + } + + Vector_t dr = this->rmax_m - this->rmin_m; + + this->hr_m = dr / this->nr_m; + + // Courant condition + this->dt_m = 0.05;//std::min(0.05, 0.5 * ( *std::min_element(this->hr_m.begin(), this->hr_m.end()) ) ); + + this->it_m = 0; + this->time_m = 0.0; + + //this->np_m = 10000; //this->nr_m[0] * this->nr_m[0]; + + this->decomp_m.fill(true); + this->isAllPeriodic_m = true; + + this->setFieldContainer(std::make_shared( + this->hr_m, this->rmin_m, this->rmax_m, this->decomp_m, this->domain_m, this->origin_m, + this->isAllPeriodic_m)); + + this->setParticleContainer(std::make_shared( + this->fcontainer_m->getMesh(), this->fcontainer_m->getFL())); + + this->fcontainer_m->initializeFields(); + + this->setFieldSolver( std::make_shared( this->solver_m, &this->fcontainer_m->getOmegaField()) ); + + this->fsolver_m->initSolver(); + + //this->setLoadBalancer( std::make_shared( this->lbt_m, this->fcontainer_m, this->pcontainer_m, this->fsolver_m) ); + + //initializeParticles(); + + //this->par2grid(); + + initializeGridVorticity(); + normalizeInitialGridCirculation(); + + auto omega0 = this->fcontainer_m->getOmegaField().deepCopy(); + this->fsolver_m->runSolver(); + this->computeVelocityField(); + logEnergyDiagnostics(); + Kokkos::deep_copy(this->fcontainer_m->getOmegaField().getView(), omega0.getView()); + logEnstrophyDiagnostics(); + logDivergenceDiagnostics(); +// this->grid2par(); + double omega_init = computeOmegaL2(); + +if (ippl::Comm->rank() == 0) { + Inform m("debug "); + m << "omega L2 after initialization = " << omega_init << endl; +} + + } + +double computeOmegaL2() { + auto& omegaField = this->fcontainer_m->getOmegaField(); + auto omega = omegaField.getView(); + double local = 0.0; + const int nghost = omegaField.getNghost(); + + Kokkos::parallel_reduce( + "omega_l2", + ippl::getRangePolicy(omega, nghost), + KOKKOS_LAMBDA(const int i, const int j, double& sum) { + sum += omega(i, j) * omega(i, j); + }, + local + ); + + double global = 0.0; + ippl::Comm->reduce(local, global, 1, std::plus()); + + return std::sqrt(global); +} + +void logPushDebug() { + auto pc = this->pcontainer_m; + + auto P_view = pc->P.getView(); + + double dt = this->dt_m; + double dx = this->hr_m[0]; + + double localMaxVel = 0.0; + double localMaxDx = 0.0; + double localMaxDy = 0.0; + + Kokkos::parallel_reduce( + "push_debug", + pc->getLocalNum(), + KOKKOS_LAMBDA(const int p, + double& maxVel, + double& maxDx, + double& maxDy) { + + double ux = P_view(p)[0]; + double uy = P_view(p)[1]; + + double vel = sqrt(ux * ux + uy * uy); + + double dispX = fabs(ux * dt); + double dispY = fabs(uy * dt); + + if (vel > maxVel) maxVel = vel; + if (dispX > maxDx) maxDx = dispX; + if (dispY > maxDy) maxDy = dispY; + + }, + Kokkos::Max(localMaxVel), + Kokkos::Max(localMaxDx), + Kokkos::Max(localMaxDy) + ); + + double globalMaxVel = localMaxVel; + double globalMaxDx = localMaxDx; + double globalMaxDy = localMaxDy; + + if (ippl::Comm->rank() == 0) { + Inform m("push_debug "); + + m << "step = " << this->it_m + << ", max velocity = " << globalMaxVel + << ", max |dx| = " << globalMaxDx + << ", max |dy| = " << globalMaxDy + << ", max displacement/dx = " + << std::max(globalMaxDx, globalMaxDy) / dx + << endl; + } +} + +void logVelocityRatioDebug() { + auto pc = this->pcontainer_m; + auto P_view = pc->P.getView(); + + double localMaxRatio = 0.0; + double localSumRatio = 0.0; + size_type localN = pc->getLocalNum(); + + Kokkos::parallel_reduce( + "velocity_ratio_debug", + localN, + KOKKOS_LAMBDA(const int p, double& maxRatio, double& sumRatio) { + double ux = fabs(P_view(p)[0]); + double uy = fabs(P_view(p)[1]); + + double ratio = uy / (ux + 1.0e-30); + + if (ratio > maxRatio) maxRatio = ratio; + sumRatio += ratio; + }, + Kokkos::Max(localMaxRatio), + localSumRatio + ); + + double globalMaxRatio = 0.0; + double globalSumRatio = 0.0; + size_type globalN = 0; + + MPI_Allreduce(&localMaxRatio, &globalMaxRatio, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); + MPI_Allreduce(&localSumRatio, &globalSumRatio, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&localN, &globalN, 1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, MPI_COMM_WORLD); + + if (ippl::Comm->rank() == 0) { + Inform m("velocity_ratio "); + m << "step = " << this->it_m + << ", max |uy|/|ux| = " << globalMaxRatio + << ", avg |uy|/|ux| = " << globalSumRatio / globalN + << endl; + } +} + +void initializeVirtualParticles(){ + clearVirtualParticles(); + auto* FL = &this->fcontainer_m->getFL(); + auto local = FL->getLocalNDIndex(); + + int i0 = local[0].first(); + int i1 = local[0].last(); + int j0 = local[1].first(); + int j1 = local[1].last(); + + size_type nx = i1 - i0 + 1; + size_type ny = j1 - j0 + 1; + size_type nlocal = nx * ny; + const int nghost = this->fcontainer_m->getOmegaField().getNghost(); + + auto pc = this->pcontainer_m; + pc->create(nlocal); + + + auto R_view = pc->R.getView(); + auto omega_p = pc->omega.getView(); + auto omega_g = this->fcontainer_m->getOmegaField().getView(); + auto P_view = pc->P.getView(); + auto u_g = this->fcontainer_m->getUField().getView(); + + + + Vector_t rmin = this->rmin_m; + Vector_t hr = this->hr_m; + double dA = hr[0] * hr[1]; + Kokkos::parallel_for( + "init_virtual_particles_from_grid", + nlocal, + KOKKOS_LAMBDA(const int p) { + int ii = p % nx; + int jj = p / nx; + + int i = i0 + ii; + int j = j0 + jj; + int li = ii + nghost; + int lj = jj + nghost; + + R_view(p)[0] = rmin[0] + (i+0.5) * hr[0]; + R_view(p)[1] = rmin[1] + (j+0.5) * hr[1]; + + omega_p(p) = omega_g(li, lj) * dA; + P_view(p) = u_g(li, lj); + } + ); + + Kokkos::fence(); +} + + +void initializeGridVorticity() { + auto& omegaField = this->fcontainer_m->getOmegaField(); + auto omega_view = omegaField.getView(); + + auto localND = this->fcontainer_m->getFL().getLocalNDIndex(); + + int i0 = localND[0].first(); + int i1 = localND[0].last(); + int j0 = localND[1].first(); + int j1 = localND[1].last(); + const int nghost = omegaField.getNghost(); + + Vector_t rmin = this->rmin_m; + Vector_t rmax = this->rmax_m; + Vector_t hr = this->hr_m; + + double y_mid = 0.5 * (rmin[1] + rmax[1]); + double y_low = y_mid - 1.0; + double y_high = y_mid + 1.0; + + Kokkos::parallel_for( + "initialize_grid_vorticity", + Kokkos::MDRangePolicy>({nghost, nghost}, + {nghost + i1 - i0 + 1, + nghost + j1 - j0 + 1}), + KOKKOS_LAMBDA(const int li, const int lj) { + int i = i0 + li - nghost; + int j = j0 + lj - nghost; + + double y = rmin[1] + (j + 0.5) * hr[1]; + double perturb = alvine::sinusoidalVorticityPerturbation(i, j); + + omega_view(li, lj) = (y >= y_low && y <= y_high) ? 1.0 + perturb : 0.0; + } + ); + + Kokkos::fence(); +} + +void normalizeInitialGridCirculation() { + double gamma = this->computeGridCirculation(); + double targetGamma = (this->rmax_m[0] - this->rmin_m[0]) * 2.0; + + if (std::fabs(gamma) < 1e-30) { + return; + } + + this->fcontainer_m->getOmegaField() = + this->fcontainer_m->getOmegaField() * (targetGamma / gamma); +} + +void pushVirtualParticlesForward() { + auto pc = this->pcontainer_m; + + // 1. Gather velocity from grid to particles + //this->grid2par(); // fills pc->P using uField + + // Debug before push: checks expected displacement = P * dt + logPushDebug(); + logVelocityRatioDebug(); + // 2. Move particles forward + auto R_view = pc->R.getView(); + auto P_view = pc->P.getView(); + + double dt = this->dt_m; + + Kokkos::parallel_for( + "push_virtual_particles_forward", + pc->getLocalNum(), + KOKKOS_LAMBDA(const int p) { + R_view(p) = R_view(p) + P_view(p) * dt; + } + ); + + Kokkos::fence(); + + // 3. Apply particle boundary conditions / update ownership + pc->update(); +} + +void clearVirtualParticles() { + auto pc = this->pcontainer_m; + size_type nlocal = pc->getLocalNum(); + + if (nlocal == 0) { + return; + } + + Kokkos::View invalid("invalid", nlocal); + + Kokkos::parallel_for( + "mark_all_particles_invalid", + nlocal, + KOKKOS_LAMBDA(const int p) { + invalid(p) = true; + } + ); + + Kokkos::fence(); + + pc->destroy(invalid, nlocal); +} + +void logEnergyDiagnostics() { + double energy = this->computeKineticEnergy(); + + if (!this->energy_initialized_m) { + this->energy0_m = energy; + this->energy_initialized_m = true; + + if (ippl::Comm->rank() == 0) { + std::ofstream out("energy.csv", std::ios::out); + out << "step,time,energy,rel_error,normalized_energy\n"; + out.close(); + } + ippl::Comm->barrier(); + } + + double relErr = this->relativeError(energy, this->energy0_m); + double normalizedEnergy = + energy / (std::fabs(this->energy0_m) > 1e-30 ? this->energy0_m : 1e-30); + + if (ippl::Comm->rank() == 0) { + Inform m("energy "); + m << "kinetic energy = " << energy + << ", relError = " << relErr + << ", normalizedEnergy = " << normalizedEnergy << endl; + + std::ofstream out("energy.csv", std::ios::app); + out.precision(16); + out.setf(std::ios::scientific, std::ios::floatfield); + out << this->it_m << "," + << this->time_m << "," + << energy << "," + << relErr << "," + << normalizedEnergy << "\n"; + out.close(); + } +} + +void logEnstrophyDiagnostics() { + double enstrophy = this->computeEnstrophy(); + + if (!this->enstrophy_initialized_m) { + this->enstrophy0_m = enstrophy; + this->enstrophy_initialized_m = true; + + if (ippl::Comm->rank() == 0) { + std::ofstream out("enstrophy.csv", std::ios::out); + out << "step,time,enstrophy,rel_error\n"; + out.close(); + } + ippl::Comm->barrier(); + } + + double relErr = this->relativeError(enstrophy, this->enstrophy0_m); + + if (ippl::Comm->rank() == 0) { + Inform m("enstrophy "); + m << "enstrophy = " << enstrophy + << ", relError = " << relErr << endl; + + std::ofstream out("enstrophy.csv", std::ios::app); + out.precision(16); + out.setf(std::ios::scientific, std::ios::floatfield); + out << this->it_m << "," + << this->time_m << "," + << enstrophy << "," + << relErr << "\n"; + out.close(); + } +} + + +void logDivergenceDiagnostics() { + double divL2 = this->computeDivergenceL2(); + + if (ippl::Comm->rank() == 0) { + Inform m("divergence "); + m << "L2 = " << divL2 << endl; + + std::ofstream out("divergence.csv", std::ios::app); + + if (this->it_m == 0) { + out << "step,time,div_l2\n"; + } + + out.precision(16); + out.setf(std::ios::scientific, std::ios::floatfield); + + out << this->it_m << "," + << this->time_m << "," + << divL2 << "\n"; + out.close(); + } +} + + void dump() override { + static IpplTimings::TimerRef dumpTimer = IpplTimings::getTimer("vtkDump"); + IpplTimings::startTimer(dumpTimer); + + alvine::vtk::writeScalarField2D("data/FSL", "omega", this->fcontainer_m->getOmegaField(), + this->rmin_m, this->hr_m, this->it_m); + + IpplTimings::stopTimer(dumpTimer); + } + + void logOmegaField() { + alvine::vtk::writeScalarFieldCsv2D("data/FSL/omega_csv", "omega", + this->fcontainer_m->getOmegaField(), this->rmin_m, + this->hr_m, this->it_m); + } + + void post_step() override { + Inform m("Step: "); + this->time_m += this->dt_m; + this->it_m++; + + this->logOmegaField(); + if (this->it_m % this->dump_freq_m == 0) { + // this->dump(); + } + + m << this->it_m << " Done" << endl; + } + + void advance() override { + advectForward(); + } + +void advectForward() { + static IpplTimings::TimerRef PTimer = + IpplTimings::getTimer("pushVelocity"); + static IpplTimings::TimerRef RTimer = + IpplTimings::getTimer("pushPosition"); + static IpplTimings::TimerRef SolveTimer = + IpplTimings::getTimer("solve"); + static IpplTimings::TimerRef par2gridTimer = + IpplTimings::getTimer("par2grid"); + + auto omega_n = this->fcontainer_m->getOmegaField().deepCopy(); + + // 1. Compute velocity u^n from omega^n + // The FFT solver writes the Poisson solution into omegaField. Restore the + // saved vorticity before creating/remapping virtual particles. + IpplTimings::startTimer(SolveTimer); + this->fsolver_m->runSolver(); + IpplTimings::stopTimer(SolveTimer); + + IpplTimings::startTimer(PTimer); + this->computeVelocityField(); + logEnergyDiagnostics(); + IpplTimings::stopTimer(PTimer); + Kokkos::deep_copy(this->fcontainer_m->getOmegaField().getView(), omega_n.getView()); + logEnstrophyDiagnostics(); + logDivergenceDiagnostics(); + + // 2. Create virtual particles from omega^n + initializeVirtualParticles(); + + // 3. Push particles using u^n + IpplTimings::startTimer(RTimer); + pushVirtualParticlesForward(); + IpplTimings::stopTimer(RTimer); + + // 4. Scatter particles to form omega^{n+1} + IpplTimings::startTimer(par2gridTimer); + this->par2grid(); + IpplTimings::stopTimer(par2gridTimer); + + // 5. Delete temporary particles + clearVirtualParticles(); +} +}; +#endif diff --git a/alvine/LoadBalancer.hpp b/alvine/LoadBalancer.hpp new file mode 100644 index 000000000..be9389482 --- /dev/null +++ b/alvine/LoadBalancer.hpp @@ -0,0 +1,108 @@ +#ifndef IPPL_LOAD_BALANCER_H +#define IPPL_LOAD_BALANCER_H + +#include "ParticleContainer.hpp" +#include + +template +class LoadBalancer{ + using Base = ippl::ParticleBase>; + using FieldSolver_t= ippl::FieldSolverBase; + using vorticity_field_type = std::conditional, VField_t>::type; + + private: + double loadbalancethreshold_m; + vorticity_field_type* omega_m; + std::shared_ptr> pc_m; + std::shared_ptr fs_m; + unsigned int loadbalancefreq_m; + ORB orb; + public: + LoadBalancer(double lbs, std::shared_ptr> &fc, std::shared_ptr> &pc, std::shared_ptr &fs) + :loadbalancethreshold_m(lbs), omega_m(&fc->getOmegaField()), pc_m(pc), fs_m(fs) {} + + ~LoadBalancer() { } + + double getLoadBalanceThreshold() const { return loadbalancethreshold_m; } + void setLoadBalanceThreshold(double threshold) { loadbalancethreshold_m = threshold; } + + Field_t* getOmega() const { return omega_m; } + void setOmega(Field_t* omega) { omega_m = omega; } + + std::shared_ptr> getParticleContainer() const { return pc_m; } + void setParticleContainer(std::shared_ptr> pc) { pc_m = pc; } + + std::shared_ptr getFieldSolver() const { return fs_m; } + void setFieldSolver(std::shared_ptr fs) { fs_m = fs; } + + void updateLayout(ippl::FieldLayout* fl, ippl::UniformCartesian* mesh, bool& isFirstRepartition) { + // Update local fields + + static IpplTimings::TimerRef tupdateLayout = IpplTimings::getTimer("updateLayout"); + IpplTimings::startTimer(tupdateLayout); + (*omega_m).updateLayout(*fl); + + // Update layout with new FieldLayout + PLayout_t* layout = &pc_m->getLayout(); + (*layout).updateLayout(*fl, *mesh); + IpplTimings::stopTimer(tupdateLayout); + static IpplTimings::TimerRef tupdatePLayout = IpplTimings::getTimer("updatePB"); + IpplTimings::startTimer(tupdatePLayout); + if (!isFirstRepartition) { + pc_m->update(); + } + IpplTimings::stopTimer(tupdatePLayout); + } + + void initializeORB(ippl::FieldLayout* fl, ippl::UniformCartesian* mesh) { + orb.initialize(*fl, *mesh, *omega_m); + } + + void repartition(ippl::FieldLayout* fl, ippl::UniformCartesian* mesh, bool& isFirstRepartition) { + // Repartition the domains + + using Base = ippl::ParticleBase>; + typename Base::particle_position_type *R; + R = &pc_m->R; + bool res = orb.binaryRepartition(*R, *fl, isFirstRepartition); + if (res != true) { + std::cout << "Could not repartition!" << std::endl; + return; + } + // Update + this->updateLayout(fl, mesh, isFirstRepartition); + if constexpr (Dim == 2 || Dim == 3) { + if (fs_m->getStype() == "FFT") { + std::get>(fs_m->getSolver()).setRhs(*omega_m); + } + } + } + + bool balance(size_type totalP, const unsigned int nstep) { + if (ippl::Comm->size() < 2) { + return false; + } + if (std::strcmp(TestName, "UniformPlasmaTest") == 0) { + return (nstep % loadbalancefreq_m == 0); + } else { + int local = 0; + std::vector res(ippl::Comm->size()); + double equalPart = (double)totalP / ippl::Comm->size(); + double dev = std::abs((double)pc_m->getLocalNum() - equalPart) / totalP; + if (dev > loadbalancethreshold_m) { + local = 1; + } + MPI_Allgather(&local, 1, MPI_INT, res.data(), 1, MPI_INT, + ippl::Comm->getCommunicator()); + + for (unsigned int i = 0; i < res.size(); i++) { + if (res[i] == 1) { + return true; + } + } + return false; + } + } +}; + +#endif diff --git a/alvine/ParticleContainer.hpp b/alvine/ParticleContainer.hpp new file mode 100644 index 000000000..1df902a0a --- /dev/null +++ b/alvine/ParticleContainer.hpp @@ -0,0 +1,46 @@ +#ifndef IPPL_PARTICLE_CONTAINER_H +#define IPPL_PARTICLE_CONTAINER_H + +#include +#include "Manager/BaseManager.h" + +// Define the ParticlesContainer class +template +class ParticleContainer : public ippl::ParticleBase>{ + using Base = ippl::ParticleBase>; + using vorticity_type = std::conditional, typename Base::particle_position_type >::type; + + public: + typename Base::particle_position_type P; + vorticity_type omega; + typename Base::particle_position_type R_old; + + private: + PLayout_t pl_m; + public: + ParticleContainer(Mesh_t& mesh, FieldLayout_t& FL) + : pl_m(FL, mesh) { + this->initialize(pl_m); + registerAttributes(); + setupBCs(); + } + + ~ParticleContainer(){} + + std::shared_ptr> getPL() { return pl_m; } + void setPL(std::shared_ptr>& pl) { pl_m = pl; } + + void registerAttributes() { + // register the particle attributes + + this->addAttribute(P); + this->addAttribute(omega); + this->addAttribute(R_old); + } + void setupBCs() { setBCAllPeriodic(); } + + private: + void setBCAllPeriodic() { this->setParticleBC(ippl::BC::PERIODIC); } +}; + +#endif diff --git a/alvine/SinusoidalJitter.hpp b/alvine/SinusoidalJitter.hpp new file mode 100644 index 000000000..693e50b32 --- /dev/null +++ b/alvine/SinusoidalJitter.hpp @@ -0,0 +1,29 @@ +#ifndef IPPL_ALVINE_SINUSOIDAL_JITTER_HPP +#define IPPL_ALVINE_SINUSOIDAL_JITTER_HPP + +#include + +namespace alvine { + +KOKKOS_INLINE_FUNCTION +double sinusoidalPositionJitter(double spacing, double xIndex, double yIndex, int component) { + double seed = Kokkos::sin(12.9898 * (xIndex + 1.0) + + 78.233 * (yIndex + 1.0) + + 37.719 * component) + * 43758.5453; + double random01 = seed - Kokkos::floor(seed); + + return spacing * 0.2 * (random01 - 0.5); +} + +KOKKOS_INLINE_FUNCTION +double sinusoidalVorticityPerturbation(double xIndex, double yIndex) { + double seed = Kokkos::sin(12.9898 * (xIndex + 1.0) + 78.233 * (yIndex + 1.0)) * 43758.5453; + double random01 = seed - Kokkos::floor(seed); + + return 0.1 * (random01 - 0.5); +} + +} // namespace alvine + +#endif diff --git a/alvine/VortexDistributions.h b/alvine/VortexDistributions.h new file mode 100644 index 000000000..f762a7468 --- /dev/null +++ b/alvine/VortexDistributions.h @@ -0,0 +1,89 @@ +#ifndef IPPL_VORTEX_IN_CELL_DITRIBUTIONS_H +#define IPPL_VORTEX_IN_CELL_DITRIBUTIONS_H + +#include + +#include "Manager/BaseManager.h" +#include "ParticleContainer.hpp" + +using view_type = typename ippl::detail::ViewType, 1>::view_type; +using host_type = typename ippl::ParticleAttrib::host_mirror_type;;/*typename ippl::ParticleAttrib::HostMirror;*/ +using vector_type = ippl::Vector; + +class BaseDistribution { +public: + view_type r; + host_type omega; + vector_type rmin, rmax, origin, center; + unsigned np; + + BaseDistribution(view_type r_, host_type omega_, vector_type r_min, vector_type r_max, + vector_type origin, unsigned np) + : r(r_) + , omega(omega_) + , rmin(r_min) + , rmax(r_max) + , origin(origin) + , np(np) { + this->center = rmin + 0.5 * (rmax - rmin); + } + + KOKKOS_INLINE_FUNCTION virtual void operator()(const size_t i) const = 0; +}; + +class UnitDisk : BaseDistribution { +public: + UnitDisk(view_type r_, host_type omega_, vector_type r_min, vector_type r_max, + vector_type origin, unsigned np) + : BaseDistribution(r_, omega_, r_min, r_max, origin, np) {} + + KOKKOS_INLINE_FUNCTION void operator()(const size_t i) const { + vector_type dist = this->r(i) - this->center; + double norm = 0.0; + for (unsigned int d = 0; d < Dim; d++) { + norm += std::pow(dist(d), 2); + } + norm = std::sqrt(norm); + float radius = 3.0; + if (norm > radius) { + this->omega(i) = 0; // 15/radius_core seemed to be too strong + } else { + //this->omega(i) = 1; + this->omega(i) = 1 * ((rmax[1] - rmin[1]) * (rmax[0] - rmin[0])) / np; + } + } +}; + +class HalfPlane : BaseDistribution { +public: + HalfPlane(view_type r_, host_type omega_, vector_type r_min, vector_type r_max, + vector_type origin, unsigned np) + : BaseDistribution(r_, omega_, r_min, r_max, origin, np) {} + + KOKKOS_INLINE_FUNCTION void operator()(const size_t i) const { + if (this->r(i)(1) > this->center(1)) { + //this->omega(i) = 1; + this->omega(i) = 1 * ((rmax[1] - rmin[1]) * (rmax[0] - rmin[0])) / np; + } else { + //this->omega(i) = -1; + this->omega(i) = -1 * ((rmax[1] - rmin[1]) * (rmax[0] - rmin[0])) / np; + } + } +}; + +class Band : BaseDistribution { +public: + Band(view_type r_, host_type omega_, vector_type r_min, vector_type r_max, vector_type origin, unsigned np) + : BaseDistribution(r_, omega_, r_min, r_max, origin, np) {} + + KOKKOS_INLINE_FUNCTION void operator()(const size_t i) const { + //if (this->r(i)(1) > this->center(1) + 1 or this->r(i)(1) < this->center(1) - 1) { + // // Outside of the band + // this->omega(i) = 0; + //} else { + this->omega(i) = (2.0 * (rmax[0] - rmin[0])) / np; + //} + } +}; + +#endif diff --git a/alvine/VortexInCell.cpp b/alvine/VortexInCell.cpp new file mode 100644 index 000000000..ff80097ac --- /dev/null +++ b/alvine/VortexInCell.cpp @@ -0,0 +1,101 @@ +// Vortex In Cell Test +// Usage: +// srun ./VortexInCell +// [...] --overallocate --info 10 +// nx = No. cell-centered points in the x-direction +// ny... = No. cell-centered points in the y-, z-, ...-direction +// Np = No. of vortex particles in the simulation +// Nt = Number of time steps +// stype = Field solver type (FFT and CG supported) +// dump_freq= Dumping frequency of particle output +// 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: +// makdir build_*/alvine/data +// chmod +x data +// srun ./VortexInCell 128 128 10000 100 FFT 100 --overallocate 1.0 --info 5 +// to build, call +// make VortexInCell +// in the build directory to only build this target + +constexpr unsigned Dim = 2; +using T = double; +const char* TestName = "VortexInCell"; + +#include "Ippl.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "datatypes.h" + +#include "Utility/IpplTimings.h" + +#include "Manager/PicManager.h" +#include "VortexInCellManager.h" +#include "VortexDistributions.h" + + +int main(int argc, char* argv[]) { + ippl::initialize(argc, argv); + { + Inform msg(TestName); + + static IpplTimings::TimerRef mainTimer = IpplTimings::getTimer("total"); + IpplTimings::startTimer(mainTimer); + + unsigned arg = 1; + Vector_t nr; + for (unsigned d = 0; d < Dim; d++) { + nr[d] = std::atoi(argv[arg++]); + } + + int np = std::atoi(argv[arg++]); + int nt = std::atoi(argv[arg++]); + std::string solver = argv[arg++]; + int dump_freq = std::atoi(argv[arg++]); + + msg << " Grid size: " << nr << " No. of particles: " << np << " No. of time steps: " << nt << endl; + + // ===== CRITICAL: Create mesh and layout with proper MPI decomposition ===== + ippl::NDIndex domain; + for (unsigned i = 0; i < Dim; i++) { + domain[i] = ippl::Index(nr[i]); + } + + // Domain bounds (adjust as needed for your vortex problem) + Vector_t rmin(0.0); + Vector_t rmax(10.0); + Vector_t origin = rmin; + Vector_t hr = (rmax - rmin) / nr; + + std::array isParallel; + isParallel.fill(true); + const bool isAllPeriodic = true; + + // Create the mesh and layout with MPI communicator + Mesh_t mesh(domain, hr, origin); + FieldLayout_t FL(MPI_COMM_WORLD, domain, isParallel, isAllPeriodic); + + // Now create manager WITH the layout info + VortexInCellManager manager(nt, nr, np, solver, dump_freq, + rmin, rmax, origin, FL, mesh); + + manager.pre_run(); + manager.run(manager.getNt()); + + IpplTimings::stopTimer(mainTimer); + IpplTimings::print(); + IpplTimings::print(std::string("timing.dat")); + } + ippl::finalize(); + + return 0; +} diff --git a/alvine/VortexInCellManager.h b/alvine/VortexInCellManager.h new file mode 100644 index 000000000..1e69c0f69 --- /dev/null +++ b/alvine/VortexInCellManager.h @@ -0,0 +1,470 @@ +#ifndef IPPL_VORTEX_IN_CELL_MANAGER_H +#define IPPL_VORTEX_IN_CELL_MANAGER_H + +#include + +#include "AlvineManager.h" +#include "FieldContainer.hpp" +#include "FieldSolver.hpp" +#include "LoadBalancer.hpp" +#include "Manager/BaseManager.h" +#include "ParticleContainer.hpp" +#include "Random/Distribution.h" +#include "Random/InverseTransformSampling.h" +#include "Random/NormalDistribution.h" +#include "Random/Randu.h" +#include "SinusoidalJitter.hpp" +#include "VortexDistributions.h" +#include "VtkDump.hpp" + +using view_type = typename ippl::detail::ViewType, 1>::view_type; +using host_type = typename ippl::ParticleAttrib::host_mirror_type; /*using host_type = typename ippl::ParticleAttrib::HostMirror;*/ + + +template +class VortexInCellManager : public AlvineManager { +public: + using ParticleContainer_t = ParticleContainer; + using FieldContainer_t = FieldContainer; + using FieldSolver_t = FieldSolver; + using LoadBalancer_t = LoadBalancer; + FieldLayout_t FL_m; // Store the field layout + Mesh_t mesh_m; // Store the mesh + + // Constructor declaration + VortexInCellManager(unsigned nt_, Vector_t& nr_, unsigned np_, + std::string& solver_, int dump_freq_, + Vector_t rmin_ = 0.0, + Vector_t rmax_ = 10.0, + Vector_t origin_ = 0.0, + FieldLayout_t& FL_ = nullptr, + Mesh_t& mesh_ = nullptr) + : AlvineManager(nt_, nr_, np_, solver_, dump_freq_) { + this->rmin_m = rmin_; + this->rmax_m = rmax_; + this->origin_m = origin_; + this->FL_m = FL_; // Store the layout + this->mesh_m = mesh_; // Store the mesh + } + + ~VortexInCellManager() {} + +void pre_run() override { + + Inform csvout(NULL, "particles.csv", Inform::OVERWRITE); + csvout.precision(16); + csvout.setf(std::ios::scientific, std::ios::floatfield); + + if constexpr (Dim == 2) { + csvout << "time,index,pos_x,pos_y,vorticity" << endl; + } else { + csvout << "time,index,pos_x,pos_y,pos_z" << endl; + } + + for (unsigned i = 0; i < Dim; i++) { + this->domain_m[i] = ippl::Index(this->nr_m[i]); + } + + + Vector_t dr = this->rmax_m - this->rmin_m; + + this->hr_m = dr / this->nr_m; + + // Courant condition + this->dt_m = 0.05;//std::min(0.05, 0.5 * ( *std::min_element(this->hr_m.begin(), this->hr_m.end()) ) ); + + this->it_m = 0; + this->time_m = 0.0; + + //this->np_m = 10000; //this->nr_m[0] * this->nr_m[0]; + + this->decomp_m.fill(true); + this->isAllPeriodic_m = true; + + this->setFieldContainer(std::make_shared( + this->hr_m, this->rmin_m, this->rmax_m, this->decomp_m, this->domain_m, this->origin_m, + this->isAllPeriodic_m)); + + this->setParticleContainer(std::make_shared( + this->fcontainer_m->getMesh(), this->fcontainer_m->getFL())); + + this->fcontainer_m->initializeFields(); + + this->setFieldSolver( std::make_shared( this->solver_m, &this->fcontainer_m->getOmegaField()) ); + + this->fsolver_m->initSolver(); + + //this->setLoadBalancer( std::make_shared( this->lbt_m, this->fcontainer_m, this->pcontainer_m, this->fsolver_m) ); + + initializeParticles(); + + this->par2grid(); + auto omega0 = this->fcontainer_m->getOmegaField().deepCopy(); + + this->fsolver_m->runSolver(); + this->computeVelocityField(); + logEnergyDiagnostics(); + Kokkos::deep_copy(this->fcontainer_m->getOmegaField().getView(), omega0.getView()); + logEnstrophyDiagnostics(); + logDivergenceDiagnostics(); + this->grid2par(); + + std::shared_ptr pc = this->pcontainer_m; + + pc->R_old = pc->R; + pc->R = pc->R_old + pc->P * this->dt_m; + pc->update(); + + } + + +void initializeParticles() { + auto* mesh = &this->fcontainer_m->getMesh(); + auto* FL = &this->fcontainer_m->getFL(); + std::shared_ptr pc = this->pcontainer_m; + ippl::detail::RegionLayout> rlayout; + const bool isFEM = (this->solver_m == "FEM") || (this->solver_m == "FEM_PRECON"); + rlayout = ippl::detail::RegionLayout>(*FL, *mesh, isFEM); + const ippl::NDIndex& local = FL->getLocalNDIndex(); + + // 1. Global lattice dimensions based on user‑supplied particle count + unsigned nxp_global = static_cast(std::sqrt(this->np_m)); + unsigned nyp_global = this->np_m / nxp_global; + size_type totalP_global = nxp_global * nyp_global; + + // 2. Physical bounds and spacing + double xmin_global = this->rmin_m[0]; + double xmax_global = this->rmax_m[0]; + double ymin_global = this->rmin_m[1]; + double ymax_global = this->rmax_m[1]; + double ymin_band = (ymin_global + ymax_global) / 2.0 - 1.0; + double ymax_band = (ymin_global + ymax_global) / 2.0 + 1.0; + + double dxp = (xmax_global - xmin_global) / nxp_global; + double dyp = (ymax_band - ymin_band) / nyp_global; + + // 3. Local domain from grid decomposition + int local_start_x = local[0].first(); + int local_end_x = local[0].last(); + int local_start_y = local[1].first(); + int local_end_y = local[1].last(); + + double xmin_local = xmin_global + local_start_x * this->hr_m[0]; + double xmax_local = xmin_global + (local_end_x + 1) * this->hr_m[0]; + double ymin_local = ymin_global + local_start_y * this->hr_m[1]; + double ymax_local = ymin_global + (local_end_y + 1) * this->hr_m[1]; + + // 4. Intersect rank's physical rectangle with the vortex band + double y_low = std::max(ymin_local, ymin_band); + double y_high = std::min(ymax_local, ymax_band); + // If the intersection is empty, the rank gets no particles + if (y_low >= y_high) { + pc->create(0); + return; + } + + // 5. Find the range of lattice indices that fall inside the rank's rectangle + int ix_start = static_cast(std::ceil((xmin_local - xmin_global - 0.5 * dxp) / dxp)); + int ix_end = static_cast(std::floor((xmax_local - xmin_global - 0.5 * dxp) / dxp)); + ix_start = std::max(0, ix_start); + ix_end = std::min(static_cast(nxp_global - 1), ix_end); + + int iy_start = static_cast(std::ceil((y_low - ymin_band - 0.5 * dyp) / dyp)); + int iy_end = static_cast(std::floor((y_high - ymin_band - 0.5 * dyp) / dyp)); + iy_start = std::max(0, iy_start); + iy_end = std::min(static_cast(nyp_global - 1), iy_end); + + unsigned nxp_local = ix_end - ix_start + 1; + unsigned nyp_local = iy_end - iy_start + 1; + size_type nlocal = nxp_local * nyp_local; + + // 6. Create particles + pc->create(nlocal); + + // 7. Fill positions on the lattice (device side) + auto R_view = pc->R.getView(); + + Kokkos::parallel_for( + "init_particle_positions", + nlocal, + KOKKOS_LAMBDA(const int i) { + unsigned ix_local = i % nxp_local; + unsigned iy_local = i / nxp_local; + unsigned ix_global = ix_start + ix_local; + unsigned iy_global = iy_start + iy_local; + + double jitter_x = alvine::sinusoidalPositionJitter(dxp, ix_global, iy_global, 0); + double jitter_y = alvine::sinusoidalPositionJitter(dyp, ix_global, iy_global, 1); + + double x = xmin_global + (ix_global + 0.5) * dxp + jitter_x; + double y = ymin_band + (iy_global + 0.5) * dyp + jitter_y; + + R_view(i)[0] = x; + R_view(i)[1] = y; + } + ); + +// 9. Particle circulation strength (2D VIC) +auto omega_view = pc->omega.getView(); +double omega0 = 1.0; // physical vorticity amplitude +double Ap = dxp * dyp; // particle area + +Kokkos::parallel_for( + "init_particle_vorticity", + nlocal, + KOKKOS_LAMBDA(const int i) { + omega_view(i) = omega0 * Ap; + } +); + + Kokkos::fence(); + +} + +void logEnergyDiagnostics() { + double energy = this->computeKineticEnergy(); + + if (!this->energy_initialized_m) { + this->energy0_m = energy; + this->energy_initialized_m = true; + + if (ippl::Comm->rank() == 0) { + std::ofstream out("energy.csv", std::ios::out); + out << "step,time,energy,rel_error,normalized_energy\n"; + out.close(); + } + ippl::Comm->barrier(); + } + + double relErr = this->relativeError(energy, this->energy0_m); + double normalizedEnergy = + energy / (std::fabs(this->energy0_m) > 1e-30 ? this->energy0_m : 1e-30); + + if (ippl::Comm->rank() == 0) { + Inform m("energy "); + m << "kinetic energy = " << energy + << ", relError = " << relErr + << ", normalizedEnergy = " << normalizedEnergy << endl; + + std::ofstream out("energy.csv", std::ios::app); + out.precision(16); + out.setf(std::ios::scientific, std::ios::floatfield); + out << this->it_m << "," + << this->time_m << "," + << energy << "," + << relErr << "," + << normalizedEnergy << "\n"; + out.close(); + } +} + +void logEnstrophyDiagnostics() { + double enstrophy = this->computeEnstrophy(); + + if (!this->enstrophy_initialized_m) { + this->enstrophy0_m = enstrophy; + this->enstrophy_initialized_m = true; + + if (ippl::Comm->rank() == 0) { + std::ofstream out("enstrophy.csv", std::ios::out); + out << "step,time,enstrophy,rel_error\n"; + out.close(); + } + ippl::Comm->barrier(); + } + + double relErr = this->relativeError(enstrophy, this->enstrophy0_m); + + if (ippl::Comm->rank() == 0) { + Inform m("enstrophy "); + m << "enstrophy = " << enstrophy + << ", relError = " << relErr << endl; + + std::ofstream out("enstrophy.csv", std::ios::app); + out.precision(16); + out.setf(std::ios::scientific, std::ios::floatfield); + out << this->it_m << "," + << this->time_m << "," + << enstrophy << "," + << relErr << "\n"; + out.close(); + } +} + + +void logDivergenceDiagnostics() { + double divL2 = this->computeDivergenceL2(); + + if (ippl::Comm->rank() == 0) { + Inform m("divergence "); + m << "L2 = " << divL2 << endl; + + std::ofstream out("divergence.csv", std::ios::app); + + if (this->it_m == 0) { + out << "step,time,div_l2\n"; + } + + out.precision(16); + out.setf(std::ios::scientific, std::ios::floatfield); + + out << this->it_m << "," + << this->time_m << "," + << divL2 << "\n"; + out.close(); + } +} + void advance() override { + LeapFrogStep(); + } + + void LeapFrogStep() { + + static IpplTimings::TimerRef PTimer = IpplTimings::getTimer("pushVelocity"); + static IpplTimings::TimerRef RTimer = IpplTimings::getTimer("pushPosition"); + static IpplTimings::TimerRef updateTimer = IpplTimings::getTimer("update"); + static IpplTimings::TimerRef SolveTimer = IpplTimings::getTimer("solve"); + static IpplTimings::TimerRef par2gridTimer = IpplTimings::getTimer("par2grid"); + static IpplTimings::TimerRef grid2parTimer = IpplTimings::getTimer("grid2par"); + + std::shared_ptr pc = this->pcontainer_m; + + // scatter the vorticity to the underlying grid + IpplTimings::startTimer(par2gridTimer); + this->par2grid(); + IpplTimings::stopTimer(par2gridTimer); + auto omega_n = this->fcontainer_m->getOmegaField().deepCopy(); + + // claculate stream function + IpplTimings::startTimer(SolveTimer); + this->fsolver_m->runSolver(); + IpplTimings::stopTimer(SolveTimer); + + // calculate velocity from stream function + IpplTimings::startTimer(PTimer); + this->computeVelocityField(); + logEnergyDiagnostics(); + Kokkos::deep_copy(this->fcontainer_m->getOmegaField().getView(), omega_n.getView()); + logEnstrophyDiagnostics(); + logDivergenceDiagnostics(); + IpplTimings::stopTimer(PTimer); + + // gather velocity field + IpplTimings::startTimer(grid2parTimer); + this->grid2par(); + IpplTimings::stopTimer(grid2parTimer); + + //drift + IpplTimings::startTimer(RTimer); + typename ippl::ParticleBase>::particle_position_type R_old_temp = pc->R_old; + + pc->R_old = pc->R; + pc->R = R_old_temp + 2 * pc->P * this->dt_m; + IpplTimings::stopTimer(RTimer); + + + IpplTimings::startTimer(updateTimer); + pc->update(); + IpplTimings::stopTimer(updateTimer); + + } +#include +#include +#include + +void dumpParticleDataPerRank() { + auto pc = this->pcontainer_m; + + auto R_host = pc->R.getHostMirror(); + auto omega_host = pc->omega.getHostMirror(); + + Kokkos::deep_copy(R_host, pc->R.getView()); + Kokkos::deep_copy(omega_host, pc->omega.getView()); + + std::stringstream fname; + fname << "particles_rank_" << ippl::Comm->rank() << ".csv"; + + bool write_header = (this->it_m == 1); + + std::ofstream csvout; + if (write_header) { + csvout.open(fname.str(), std::ios::out); + } else { + csvout.open(fname.str(), std::ios::app); + } + + if constexpr (Dim == 2) { + if (write_header) { + csvout << "time,index,pos_x,pos_y,vorticity\n"; + } + + for (size_type i = 0; i < pc->getLocalNum(); i++) { + csvout << this->it_m << "," << i + << "," << R_host(i)[0] + << "," << R_host(i)[1] + << "," << omega_host(i) << "\n"; + } + } else { + if (write_header) { + csvout << "time,index,pos_x,pos_y,pos_z\n"; + } + + for (size_type i = 0; i < pc->getLocalNum(); i++) { + csvout << this->it_m << "," << i + << "," << R_host(i)[0] + << "," << R_host(i)[1] + << "," << R_host(i)[2] << "\n"; + } + } + + csvout.close(); + ippl::Comm->barrier(); +} + + void dump() override { + static IpplTimings::TimerRef dumpTimer = IpplTimings::getTimer("vtkDump"); + IpplTimings::startTimer(dumpTimer); + + this->par2grid(); + auto omega_current = this->fcontainer_m->getOmegaField().deepCopy(); + this->fsolver_m->runSolver(); + this->computeVelocityField(); + Kokkos::deep_copy(this->fcontainer_m->getOmegaField().getView(), omega_current.getView()); + + alvine::vtk::writeScalarField2D("data/VortexInCell", "omega", + this->fcontainer_m->getOmegaField(), + this->rmin_m, this->hr_m, this->it_m); + + IpplTimings::stopTimer(dumpTimer); + } + + void logOmegaField() { + this->par2grid(); + alvine::vtk::writeScalarFieldCsv2D("data/VortexInCell/omega_csv", "omega", + this->fcontainer_m->getOmegaField(), this->rmin_m, + this->hr_m, this->it_m); + } + + void post_step() override { + Inform m("Step: "); + this->time_m += this->dt_m; + this->it_m++; + + this->logOmegaField(); + if (this->it_m % this->dump_freq_m == 0) { + // this->dump(); + } + + m << this->it_m << " Done" << endl; + } + + /* void dump() override { + static IpplTimings::TimerRef dumpTimer = IpplTimings::getTimer("dump"); + IpplTimings::startTimer(dumpTimer); + dumpParticleDataPerRank(); + IpplTimings::stopTimer(dumpTimer); + + }*/ + +}; +#endif diff --git a/alvine/VtkDump.hpp b/alvine/VtkDump.hpp new file mode 100644 index 000000000..4e5b45d99 --- /dev/null +++ b/alvine/VtkDump.hpp @@ -0,0 +1,223 @@ +#ifndef IPPL_ALVINE_VTK_DUMP_HPP +#define IPPL_ALVINE_VTK_DUMP_HPP + +#include + +#include +#include +#include +#include +#include +#include + +#include "Ippl.h" + +namespace alvine::vtk { + +inline std::string stepString(int step) { + std::ostringstream os; + os << std::setw(6) << std::setfill('0') << step; + return os.str(); +} + +inline std::filesystem::path legacyFileName(const std::string& outputDir, const std::string& name, + int step) { + std::ostringstream os; + os << name << "_" << stepString(step); + if (ippl::Comm->size() > 1) { + os << "_rank" << ippl::Comm->rank(); + } + os << ".vtk"; + return std::filesystem::path(outputDir) / os.str(); +} + +inline std::filesystem::path csvFileName(const std::string& outputDir, const std::string& name, + int step) { + std::ostringstream os; + os << name << "_" << stepString(step); + if (ippl::Comm->size() > 1) { + os << "_rank" << ippl::Comm->rank(); + } + os << ".csv"; + return std::filesystem::path(outputDir) / os.str(); +} + +inline void ensureOutputDirectory(const std::string& outputDir) { + if (ippl::Comm->rank() == 0) { + std::filesystem::create_directories(outputDir); + } + ippl::Comm->barrier(); +} + +template +void writeScalarField2D(const std::string& outputDir, const std::string& name, FieldType& field, + const VectorType& origin, const VectorType& spacing, int step) { + static_assert(FieldType::dim == 2, "Legacy Alvine VTK output expects a 2D field."); + + ensureOutputDirectory(outputDir); + + auto host = field.getHostMirror(); + Kokkos::deep_copy(host, field.getView()); + + const auto& local = field.getLayout().getLocalNDIndex(); + const int nghost = field.getNghost(); + const int nx = local[0].last() - local[0].first() + 1; + const int ny = local[1].last() - local[1].first() + 1; + + const auto file = legacyFileName(outputDir, name, step); + std::ofstream vtkout(file, std::ios::out); + if (!vtkout) { + throw std::runtime_error("Could not open VTK file: " + file.string()); + } + + vtkout.precision(10); + vtkout.setf(std::ios::scientific, std::ios::floatfield); + + vtkout << "# vtk DataFile Version 2.0\n"; + vtkout << name << "\n"; + vtkout << "ASCII\n"; + vtkout << "DATASET STRUCTURED_POINTS\n"; + vtkout << "DIMENSIONS " << nx + 1 << " " << ny + 1 << " 2\n"; + vtkout << "ORIGIN " << origin[0] + local[0].first() * spacing[0] << " " + << origin[1] + local[1].first() * spacing[1] << " 0\n"; + vtkout << "SPACING " << spacing[0] << " " << spacing[1] << " 1\n"; + vtkout << "CELL_DATA " << nx * ny << "\n"; + vtkout << "SCALARS " << name << " float\n"; + vtkout << "LOOKUP_TABLE default\n"; + + for (int j = local[1].first(); j <= local[1].last(); ++j) { + for (int i = local[0].first(); i <= local[0].last(); ++i) { + const int li = i - local[0].first() + nghost; + const int lj = j - local[1].first() + nghost; + vtkout << host(li, lj) << "\n"; + } + } +} + +template +void writeScalarFieldCsv2D(const std::string& outputDir, const std::string& name, FieldType& field, + const VectorType& origin, const VectorType& spacing, int step) { + static_assert(FieldType::dim == 2, "Alvine CSV output expects a 2D field."); + + ensureOutputDirectory(outputDir); + + auto host = field.getHostMirror(); + Kokkos::deep_copy(host, field.getView()); + + const auto& local = field.getLayout().getLocalNDIndex(); + const int nghost = field.getNghost(); + + const auto file = csvFileName(outputDir, name, step); + std::ofstream csvout(file, std::ios::out); + if (!csvout) { + throw std::runtime_error("Could not open CSV file: " + file.string()); + } + + csvout.precision(16); + csvout.setf(std::ios::scientific, std::ios::floatfield); + csvout << "step,i,j,x,y," << name << "\n"; + + for (int j = local[1].first(); j <= local[1].last(); ++j) { + for (int i = local[0].first(); i <= local[0].last(); ++i) { + const int li = i - local[0].first() + nghost; + const int lj = j - local[1].first() + nghost; + const double x = origin[0] + (i + 0.5) * spacing[0]; + const double y = origin[1] + (j + 0.5) * spacing[1]; + csvout << step << "," << i << "," << j << "," << x << "," << y << "," + << host(li, lj) << "\n"; + } + } +} + +template +void writeVectorField2D(const std::string& outputDir, const std::string& name, FieldType& field, + const VectorType& origin, const VectorType& spacing, int step) { + static_assert(FieldType::dim == 2, "Legacy Alvine VTK output expects a 2D field."); + + ensureOutputDirectory(outputDir); + + auto host = field.getHostMirror(); + Kokkos::deep_copy(host, field.getView()); + + const auto& local = field.getLayout().getLocalNDIndex(); + const int nghost = field.getNghost(); + const int nx = local[0].last() - local[0].first() + 1; + const int ny = local[1].last() - local[1].first() + 1; + + const auto file = legacyFileName(outputDir, name, step); + std::ofstream vtkout(file, std::ios::out); + if (!vtkout) { + throw std::runtime_error("Could not open VTK file: " + file.string()); + } + + vtkout.precision(10); + vtkout.setf(std::ios::scientific, std::ios::floatfield); + + vtkout << "# vtk DataFile Version 2.0\n"; + vtkout << name << "\n"; + vtkout << "ASCII\n"; + vtkout << "DATASET STRUCTURED_POINTS\n"; + vtkout << "DIMENSIONS " << nx + 1 << " " << ny + 1 << " 2\n"; + vtkout << "ORIGIN " << origin[0] + local[0].first() * spacing[0] << " " + << origin[1] + local[1].first() * spacing[1] << " 0\n"; + vtkout << "SPACING " << spacing[0] << " " << spacing[1] << " 1\n"; + vtkout << "CELL_DATA " << nx * ny << "\n"; + vtkout << "VECTORS " << name << " float\n"; + + for (int j = local[1].first(); j <= local[1].last(); ++j) { + for (int i = local[0].first(); i <= local[0].last(); ++i) { + const int li = i - local[0].first() + nghost; + const int lj = j - local[1].first() + nghost; + vtkout << host(li, lj)[0] << "\t" << host(li, lj)[1] << "\t0\n"; + } + } +} + +template +void writeParticles2D(const std::string& outputDir, const std::string& name, + ParticleContainerType& particles, int step) { + ensureOutputDirectory(outputDir); + + auto rHost = particles.R.getHostMirror(); + auto omegaHost = particles.omega.getHostMirror(); + + Kokkos::deep_copy(rHost, particles.R.getView()); + Kokkos::deep_copy(omegaHost, particles.omega.getView()); + + const auto nlocal = particles.getLocalNum(); + const auto file = legacyFileName(outputDir, name, step); + + std::ofstream vtkout(file, std::ios::out); + if (!vtkout) { + throw std::runtime_error("Could not open VTK file: " + file.string()); + } + + vtkout.precision(10); + vtkout.setf(std::ios::scientific, std::ios::floatfield); + + vtkout << "# vtk DataFile Version 2.0\n"; + vtkout << name << "\n"; + vtkout << "ASCII\n"; + vtkout << "DATASET POLYDATA\n"; + vtkout << "POINTS " << nlocal << " float\n"; + + for (ippl::detail::size_type i = 0; i < nlocal; ++i) { + vtkout << rHost(i)[0] << "\t" << rHost(i)[1] << "\t0\n"; + } + + vtkout << "VERTICES " << nlocal << " " << 2 * nlocal << "\n"; + for (ippl::detail::size_type i = 0; i < nlocal; ++i) { + vtkout << "1 " << i << "\n"; + } + + vtkout << "POINT_DATA " << nlocal << "\n"; + vtkout << "SCALARS omega float\n"; + vtkout << "LOOKUP_TABLE default\n"; + for (ippl::detail::size_type i = 0; i < nlocal; ++i) { + vtkout << omegaHost(i) << "\n"; + } +} + +} // namespace alvine::vtk + +#endif diff --git a/alvine/datatypes.h b/alvine/datatypes.h new file mode 100644 index 000000000..fe680a076 --- /dev/null +++ b/alvine/datatypes.h @@ -0,0 +1,64 @@ +#include "PoissonSolvers/FFTOpenPoissonSolver.h" +#include "PoissonSolvers/FFTPeriodicPoissonSolver.h" +//#include "PoissonSolvers/P3MSolver.h" +#include "PoissonSolvers/PoissonCG.h" + +// some typedefs +template +using Mesh_t = ippl::UniformCartesian; + +template +using PLayout_t = typename ippl::ParticleSpatialLayout>; + +template +using Centering_t = typename Mesh_t::DefaultCentering; + +template +using FieldLayout_t = ippl::FieldLayout; + +using size_type = ippl::detail::size_type; + +template +using Vector = ippl::Vector; + +template +using Field = ippl::Field, Centering_t, ViewArgs...>; + +template +using ORB = ippl::OrthogonalRecursiveBisection, T>; + +template +using ParticleAttrib = ippl::ParticleAttrib; + +template +using Vector_t = ippl::Vector; + +template +using Field_t = Field; + +template +using VField_t = Field, Dim, ViewArgs...>; + +template +using CGSolver_t = ippl::PoissonCG, Field_t>; + +using ippl::detail::ConditionalType, ippl::detail::VariantFromConditionalTypes; + +template +using FFTSolver_t = ConditionalType, Field_t>>; + +/*template +using P3MSolver_t = ConditionalType, Field_t>>;*/ + +template +using OpenSolver_t = + ConditionalType, Field_t>>; + +/*template +using Solver_t = VariantFromConditionalTypes, FFTSolver_t, + P3MSolver_t, OpenSolver_t>; +*/ +const double pi = Kokkos::numbers::pi_v; + +extern const char* TestName; diff --git a/plotting/particle_positions.py b/plotting/particle_positions.py new file mode 100644 index 000000000..340faf236 --- /dev/null +++ b/plotting/particle_positions.py @@ -0,0 +1,65 @@ +import matplotlib.pyplot as plt +import pandas as pd +import numpy as np +import glob +from matplotlib.animation import FuncAnimation + +# Change this path to the path of the txt file +PATH = "../" + +# Find all rank files and combine them +rank_files = glob.glob(f'{PATH}/particles_rank_*.csv') +print(f"Found {len(rank_files)} rank files: {rank_files}") + +# Load and combine data from all rank files +dfs = [] +for file in rank_files: + df = pd.read_csv(file) + dfs.append(df) + +df = pd.concat(dfs, ignore_index=True) +# ... after loading and concatenating +df = pd.concat(dfs, ignore_index=True) + +# Convert vorticity column to numeric, coercing errors to NaN +df['vorticity'] = pd.to_numeric(df['vorticity'], errors='coerce') + +# Optionally drop rows with NaN vorticity (if any) +df.dropna(subset=['vorticity'], inplace=True) + +# Then continue as before +df = df.sort_values(['time', 'index']).reset_index(drop=True) + +# ... rest of the code +# Sort by time and index to ensure proper ordering +df = df.sort_values(['time', 'index']).reset_index(drop=True) + +print(f"Combined data: {len(df)} total rows") +print(f"Time steps: {df['time'].nunique()}") +print(f"Particles per time step: {df.groupby('time').size().iloc[0] if len(df) > 0 else 0}") + +# Find the unique times and indices for particles +times = df['time'].unique() +particle_indices = df['index'].unique() + +# Create a figure and axis +fig, ax = plt.subplots() +ax.set_xlim((0,10)) +ax.set_ylim((0,10)) + +# Initialize scatter plot with empty data and color map +scat = ax.scatter([], [], c=[], s=2) + +def update(frame): + current_time = times[frame] + frame_data = df[df['time'] == current_time] + scat.set_offsets(np.c_[frame_data['pos_x'], frame_data['pos_y']]) + scat.set_array(frame_data['vorticity']) # Set colors based on vorticity + return scat, + +# Create animation +ani = FuncAnimation(fig, update, frames=len(times), blit=True, interval=50) + +# Show animation +# plt.show() +ani.save(f'particles_combined.gif', fps=30) diff --git a/plotting/plot.error b/plotting/plot.error new file mode 100644 index 000000000..731377e2b --- /dev/null +++ b/plotting/plot.error @@ -0,0 +1,35 @@ ++ '[' -z '' ']' ++ case "$-" in ++ __lmod_vx=x ++ '[' -n x ']' ++ set +x +Shell debugging temporarily silenced: export LMOD_SH_DBG_ON=1 for this output (/p/software/juwelsbooster/lmod/8.7.64/init/bash) +Shell debugging restarted ++ unset __lmod_vx ++ source /p/home/jusers/mostafa3/juwels/ippl/plotting/.venv/bin/activate +++ deactivate nondestructive +++ '[' -n '' ']' +++ '[' -n '' ']' +++ '[' -n /bin/bash -o -n '' ']' +++ hash -r +++ '[' -n '' ']' +++ unset VIRTUAL_ENV +++ '[' '!' nondestructive = nondestructive ']' +++ VIRTUAL_ENV=/p/home/jusers/mostafa3/juwels/ippl/plotting/.venv +++ export VIRTUAL_ENV +++ _OLD_VIRTUAL_PATH=/p/software/default/stages/2026/software/nano/8.5-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/NCCL/default-GCCcore-14.3.0-CUDA-13/bin:/p/software/default/stages/2026/software/CMake/4.0.3-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/libarchive/3.8.1-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/zstd/1.5.7-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/lz4/1.10.0-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/gzip/1.14-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/cURL/8.14.1-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/libpsl/0.21.5-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/libiconv/1.18-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/libidn2/2.3.8-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/bzip2/1.0.8-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/ncurses/6.5-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/OpenMPI/5.0.8-GCC-14.3.0/bin:/p/software/default/stages/2026/software/UCC/default-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/PMIx/5.0.8-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/UCX/default-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/libevent/2.1.12-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/OpenSSL/3/bin:/p/software/default/stages/2026/software/hwloc/2.12.1-GCCcore-14.3.0/sbin:/p/software/default/stages/2026/software/hwloc/2.12.1-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/libxml2/2.14.3-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/XZ/5.8.1-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/numactl/2.0.19-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/CUDA/13/nvvm/bin:/p/software/default/stages/2026/software/CUDA/13/bin:/p/software/default/stages/2026/software/GCCcore/14.3.0/bin:/p/software/default/stages/2026/software/make/4.4.1-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/binutils/2.44-GCCcore-14.3.0/bin:/usr/local/bin:/usr/bin:/usr/local/sbin:/usr/sbin:/opt/jsc/bin:/usr/local/jsc/bin:/opt/parastation/bin:/p/software/default/bin +++ PATH=/p/home/jusers/mostafa3/juwels/ippl/plotting/.venv/bin:/p/software/default/stages/2026/software/nano/8.5-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/NCCL/default-GCCcore-14.3.0-CUDA-13/bin:/p/software/default/stages/2026/software/CMake/4.0.3-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/libarchive/3.8.1-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/zstd/1.5.7-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/lz4/1.10.0-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/gzip/1.14-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/cURL/8.14.1-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/libpsl/0.21.5-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/libiconv/1.18-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/libidn2/2.3.8-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/bzip2/1.0.8-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/ncurses/6.5-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/OpenMPI/5.0.8-GCC-14.3.0/bin:/p/software/default/stages/2026/software/UCC/default-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/PMIx/5.0.8-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/UCX/default-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/libevent/2.1.12-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/OpenSSL/3/bin:/p/software/default/stages/2026/software/hwloc/2.12.1-GCCcore-14.3.0/sbin:/p/software/default/stages/2026/software/hwloc/2.12.1-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/libxml2/2.14.3-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/XZ/5.8.1-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/numactl/2.0.19-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/CUDA/13/nvvm/bin:/p/software/default/stages/2026/software/CUDA/13/bin:/p/software/default/stages/2026/software/GCCcore/14.3.0/bin:/p/software/default/stages/2026/software/make/4.4.1-GCCcore-14.3.0/bin:/p/software/default/stages/2026/software/binutils/2.44-GCCcore-14.3.0/bin:/usr/local/bin:/usr/bin:/usr/local/sbin:/usr/sbin:/opt/jsc/bin:/usr/local/jsc/bin:/opt/parastation/bin:/p/software/default/bin +++ export PATH +++ '[' -n '' ']' +++ '[' -z '' ']' +++ _OLD_VIRTUAL_PS1= +++ PS1='(.venv) ' +++ export PS1 +++ '[' -n /bin/bash -o -n '' ']' +++ hash -r ++ python /p/home/jusers/mostafa3/juwels/ippl/plotting/particle_positions.py +/p/home/jusers/mostafa3/juwels/ippl/plotting/particle_positions.py:20: FutureWarning: The behavior of DataFrame concatenation with empty or all-NA entries is deprecated. In a future version, this will no longer exclude empty or all-NA columns when determining the result dtypes. To retain the old behavior, exclude the relevant entries before the concat operation. + df = pd.concat(dfs, ignore_index=True) +/p/home/jusers/mostafa3/juwels/ippl/plotting/particle_positions.py:22: FutureWarning: The behavior of DataFrame concatenation with empty or all-NA entries is deprecated. In a future version, this will no longer exclude empty or all-NA columns when determining the result dtypes. To retain the old behavior, exclude the relevant entries before the concat operation. + df = pd.concat(dfs, ignore_index=True) +MovieWriter ffmpeg unavailable; using Pillow instead. diff --git a/plotting/plot.out b/plotting/plot.out new file mode 100644 index 000000000..e69de29bb diff --git a/plotting/requirements.txt b/plotting/requirements.txt new file mode 100644 index 000000000..061f0734c --- /dev/null +++ b/plotting/requirements.txt @@ -0,0 +1,4 @@ +Pandas +matplotlib +numpy +imageio diff --git a/plotting/setup_venv.sh b/plotting/setup_venv.sh new file mode 100755 index 000000000..81add56d3 --- /dev/null +++ b/plotting/setup_venv.sh @@ -0,0 +1,7 @@ +#!/usr/bin/env bash +set -evo pipefail +python3 -m venv .venv +source .venv/bin/activate +pip install --upgrade pip +pip install -r requirements.txt +echo "Done. Activate with: source .venv/bin/activate" diff --git a/timing.dat b/timing.dat new file mode 100644 index 000000000..7c34e14f9 --- /dev/null +++ b/timing.dat @@ -0,0 +1,28 @@ + ranks Wall tot +===================================== +total............... 8 135.9 + + ranks Wall max Wall min Wall avg +========================================================= +total............... 8 135.9 135.9 135.9 +findInternal........ 8 8.064e-05 4.933e-05 6.52e-05 +findPeriodic........ 8 8.466e-05 6.226e-05 6.982e-05 +scatter............. 8 0.4892 0.2079 0.323 +accumulateHalo...... 8 17.53 8.394 11.87 +fillHalo............ 8 2.52 2.3 2.353 +gather.............. 8 0.3393 0.1474 0.2261 +particleBC.......... 8 0.2042 0.1215 0.1531 +updateParticle...... 8 77.6 68.32 73.97 +locateParticles..... 8 12.86 6.146 9.2 +neighborSearch...... 8 1.479 0.6279 0.9633 +nonNeighboringParti. 8 0.01173 0.01126 0.01157 +sendPreprocess...... 8 9.363 1.997 6.051 +particleSend........ 8 48.69 42.75 44.8 +particleDestroy..... 8 0.956 0.7294 0.826 +particleRecv........ 8 18.51 7.345 13.01 +pushVelocity........ 8 10.07 8.848 9.312 +pushPosition........ 8 0.2197 0.1073 0.1424 +update.............. 8 78.45 68.96 74.78 +solve............... 8 35.03 33.82 34.66 +par2grid............ 8 19.3 9.769 13.26 +grid2par............ 8 3.006 2.553 2.703 diff --git a/timing_fsl.dat b/timing_fsl.dat new file mode 100644 index 000000000..400ba3b0e --- /dev/null +++ b/timing_fsl.dat @@ -0,0 +1,26 @@ + ranks Wall tot +===================================== +total............... 8 109.2 + + ranks Wall max Wall min Wall avg +========================================================= +total............... 8 109.2 109.1 109.1 +findInternal........ 8 8.199e-05 4.951e-05 6.533e-05 +findPeriodic........ 8 0.0001248 6.02e-05 8.12e-05 +pushVelocity........ 8 15.03 13.45 14.1 +pushPosition........ 8 18.37 18.16 18.26 +solve............... 8 69.03 66.98 68.05 +par2grid............ 8 3.926 3.538 3.68 +fillHalo............ 8 2.439 2.342 2.386 +gather.............. 8 0.1444 0.1356 0.1399 +particleBC.......... 8 0.1565 0.1518 0.154 +updateParticle...... 8 12.46 12.41 12.45 +locateParticles..... 8 11.98 11.16 11.58 +neighborSearch...... 8 1.011 0.8896 0.9457 +nonNeighboringParti. 8 0.01199 0.01136 0.01164 +sendPreprocess...... 8 1.113 0.3346 0.7261 +particleSend........ 8 0.02682 0.0247 0.02552 +particleDestroy..... 8 0.02234 0.02128 0.02178 +particleRecv........ 8 0.01135 0.01079 0.01109 +scatter............. 8 0.1702 0.1592 0.1624 +accumulateHalo...... 8 2.477 2.41 2.443