Resolve "Implement electrostatic particle-in-Fourier"

CMakeLists.txt

@@ -23,8 +23,12 @@ add_compile_options (-Wunused)
 add_compile_options (-Wextra)
 add_compile_options (-Werror)
 
+
 # allow deprecated functions
 add_compile_options (-Wno-deprecated-declarations)
+add_compile_options (-Wno-stringop-overflow)
+add_compile_options (-Wno-array-bounds)
+add_compile_options (-Wno-restrict)
 
 option (USE_STATIC_LIBRARIES "Link with static libraries if available" ON)
 
@@ -62,6 +66,13 @@ if (ENABLE_FFT)
     message (STATUS "Found Heffte_DIR: ${Heffte_DIR}")
 endif ()
 
+option (ENABLE_NUFFT "Enable NUFFT transform" OFF)
+if (ENABLE_NUFFT)
+    add_definitions (-DENABLE_NUFFT)
+    find_package(CUFINUFFT REQUIRED)
+    message (STATUS "Found CUFINUFFT_DIR: ${CUFINUFFT_DIR}")
+endif ()
+
 option (ENABLE_SOLVERS "Enable IPPL solvers" OFF)
 
 add_subdirectory (src)

CMakeModules/FindCUFINUFFT.cmake

@@ -0,0 +1,32 @@
+#
+# Find CUFINUFFT includes and library
+#
+# CUFINUFFT_INCLUDE_DIR - where to find cufinufft.h
+# CUFINUFFT_LIBRARY     - libcufinufft.so path
+# CUFINUFFT_FOUND       - do not attempt to use if "no" or undefined.
+
+FIND_PATH(CUFINUFFT_INCLUDE_DIR cufinufft.h
+    HINTS $ENV{CUFINUFFT_INCLUDE_PATH} $ENV{CUFINUFFT_INCLUDE_DIR} $ENV{CUFINUFFT_PREFIX}/include $ENV{CUFINUFFT_DIR}/include ${PROJECT_SOURCE_DIR}/include
+    PATHS ENV CPP_INCLUDE_PATH
+)
+#Static library has some issues and gives a cuda error at the end of compilation
+FIND_LIBRARY(CUFINUFFT_LIBRARY_DIR libcufinufft.so
+    HINTS $ENV{CUFINUFFT_LIBRARY_PATH} $ENV{CUFINUFFT_LIBRARY_DIR} $ENV{CUFINUFFT_PREFIX}/lib $ENV{CUFINUFFT_DIR}/lib $ENV{CUFINUFFT}/lib ${PROJECT_SOURCE_DIR}/lib
+    PATHS ENV LIBRARY_PATH
+)
+
+IF(CUFINUFFT_INCLUDE_DIR AND CUFINUFFT_LIBRARY_DIR)
+    SET( CUFINUFFT_FOUND "YES" )
+    SET( CUFINUFFT_DIR $ENV{CUFINUFFT_DIR} )
+ENDIF()
+
+IF (CUFINUFFT_FOUND)
+   IF (NOT CUFINUFFT_FIND_QUIETLY)
+       MESSAGE(STATUS "Found cufinufft library dir: ${CUFINUFFT_LIBRARY_DIR}")
+       MESSAGE(STATUS "Found cufinufft include dir: ${CUFINUFFT_INCLUDE_DIR}")
+   ENDIF (NOT CUFINUFFT_FIND_QUIETLY)
+ELSE (CUFINUFFT_FOUND)
+    IF (CUFINUFFT_FIND_REQUIRED)
+      MESSAGE(FATAL_ERROR "Could not find CUFINUFFT!")
+  ENDIF (CUFINUFFT_FIND_REQUIRED)
+ENDIF (CUFINUFFT_FOUND)

alpine/CMakeLists.txt

@@ -1,29 +1,22 @@
-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}/..
-)
-
-set (IPPL_LIBS ippl ${MPI_CXX_LIBRARIES})
-set (COMPILE_FLAGS ${OPAL_CXX_FLAGS})
-
-add_executable (PenningTrap PenningTrap.cpp)
-target_link_libraries (PenningTrap ${IPPL_LIBS})
-
-add_executable (UniformPlasmaTest UniformPlasmaTest.cpp)
-target_link_libraries (UniformPlasmaTest ${IPPL_LIBS})
-
-add_executable (LandauDamping LandauDamping.cpp)
-target_link_libraries (LandauDamping ${IPPL_LIBS})
-
-add_executable (BumponTailInstability BumponTailInstability.cpp)
-target_link_libraries (BumponTailInstability ${IPPL_LIBS})
+macro(list_subdirectories retval curdir)
+  file(GLOB sub-dir RELATIVE ${curdir} *)
+  set(list_of_dirs "")
+  foreach(dir ${sub-dir})
+    if(IS_DIRECTORY ${curdir}/${dir})
+      set(list_of_dirs ${list_of_dirs} ${dir})
+    endif()
+  endforeach()
+  set(${retval} ${list_of_dirs})
+endmacro()
+
+#list_subdirectories("TESTS" ${CMAKE_CURRENT_SOURCE_DIR})
+#foreach (test ${TESTS})
+#  add_subdirectory (${test})
+#endforeach()
+
+add_subdirectory (ElectrostaticPIC)
+add_subdirectory (ElectrostaticPIF)
+add_subdirectory (PinT)
 
 # vi: set et ts=4 sw=4 sts=4:
 

alpine/BumponTailInstability.cpp → alpine/ElectrostaticPIC/BumponTailInstability.cpp

@@ -252,7 +252,7 @@ int main(int argc, char *argv[]){
 
     Vector_t hr = {dx, dy, dz};
     Vector_t origin = {rmin[0], rmin[1], rmin[2]};
-    const double dt = 0.5*dx;//0.05
+    const double dt = std::atof(argv[9]);;//0.5*dx;
 
     const bool isAllPeriodic=true;
     Mesh_t mesh(domain, hr, origin);
@@ -383,6 +383,7 @@ int main(int argc, char *argv[]){
 
     IpplTimings::startTimer(dumpDataTimer);
     P->dumpBumponTail();
+    P->dumpEnergy(totalP);
     P->gatherStatistics(totalP);
     //P->dumpLocalDomains(FL, 0);
     IpplTimings::stopTimer(dumpDataTimer);
@@ -442,6 +443,7 @@ int main(int argc, char *argv[]){
         P->time_m += dt;
         IpplTimings::startTimer(dumpDataTimer);
         P->dumpBumponTail();
+        P->dumpEnergy(totalP);
         P->gatherStatistics(totalP);
         IpplTimings::stopTimer(dumpDataTimer);
         msg << "Finished time step: " << it+1 << " time: " << P->time_m << endl;

alpine/ElectrostaticPIC/CMakeLists.txt

@@ -0,0 +1,35 @@
+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}/..
+)
+
+set (IPPL_LIBS ippl ${MPI_CXX_LIBRARIES})
+set (COMPILE_FLAGS ${OPAL_CXX_FLAGS})
+
+add_executable (PenningTrap PenningTrap.cpp)
+target_link_libraries (PenningTrap ${IPPL_LIBS})
+
+add_executable (UniformPlasmaTest UniformPlasmaTest.cpp)
+target_link_libraries (UniformPlasmaTest ${IPPL_LIBS})
+
+add_executable (LandauDamping LandauDamping.cpp)
+target_link_libraries (LandauDamping ${IPPL_LIBS})
+
+add_executable (BumponTailInstability BumponTailInstability.cpp)
+target_link_libraries (BumponTailInstability ${IPPL_LIBS})
+
+# 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:

alpine/ChargedParticles.hpp → alpine/ElectrostaticPIC/ChargedParticles.hpp

@@ -508,7 +508,126 @@ class ChargedParticles : public ippl::ParticleBase<PLayout> {
 
         Ippl::Comm->barrier();
      }
-
+
+     void dumpLandauParticle(size_type totalP) {
+
+        auto Eview = E.getView();
+
+        double fieldEnergy, ExAmp;
+        double temp = 0.0;
+
+        Kokkos::parallel_reduce("Ex energy", this->getLocalNum(),
+                                KOKKOS_LAMBDA(const int i, double& valL){
+                                    double myVal = Eview(i)[0] * Eview(i)[0];
+                                    valL += myVal;
+                                }, Kokkos::Sum<double>(temp));
+
+        double globaltemp = 0.0;
+        MPI_Reduce(&temp, &globaltemp, 1, MPI_DOUBLE, MPI_SUM, 0, Ippl::getComm());
+        double volume = (rmax_m[0] - rmin_m[0]) * (rmax_m[1] - rmin_m[1]) * (rmax_m[2] - rmin_m[2]);
+        fieldEnergy = globaltemp * volume / totalP ;
+
+        double tempMax = 0.0;
+        Kokkos::parallel_reduce("Ex max norm", this->getLocalNum(),
+                                KOKKOS_LAMBDA(const size_t i, double& valL)
+                                {
+                                    double myVal = std::fabs(Eview(i)[0]);
+                                    if(myVal > valL) valL = myVal;
+                                }, Kokkos::Max<double>(tempMax));
+        ExAmp = 0.0;
+        MPI_Reduce(&tempMax, &ExAmp, 1, MPI_DOUBLE, MPI_MAX, 0, Ippl::getComm());
+
+
+        if (Ippl::Comm->rank() == 0) {
+            std::stringstream fname;
+            fname << "data/FieldLandau_";
+            fname << Ippl::Comm->size();
+            fname << ".csv";
+
+
+            Inform csvout(NULL, fname.str().c_str(), Inform::APPEND);
+            csvout.precision(10);
+            csvout.setf(std::ios::scientific, std::ios::floatfield);
+
+            if(time_m == 0.0) {
+                csvout << "time, Ex_field_energy, Ex_max_norm" << endl;
+            }
+
+            csvout << time_m << " "
+                   << fieldEnergy << " "
+                   << ExAmp << endl;
+
+        }
+
+        Ippl::Comm->barrier();
+     }
+
+
+     void dumpEnergy(size_type /*totalP*/) {
+
+
+        double potentialEnergy, kineticEnergy;
+        //auto Eview = E.getView();
+        double temp = 0.0;
+
+        //Kokkos::parallel_reduce("Potential energy", this->getLocalNum(),
+        //                        KOKKOS_LAMBDA(const int i, double& valL){
+        //                            double myVal = dot(Eview(i), Eview(i)).apply();
+        //                            valL += myVal;
+        //                        }, Kokkos::Sum<double>(temp));
+
+        double globaltemp = 0.0;
+        //MPI_Reduce(&temp, &globaltemp, 1, MPI_DOUBLE, MPI_SUM, 0, Ippl::getComm());
+        //double volume = (rmax_m[0] - rmin_m[0]) * (rmax_m[1] - rmin_m[1]) * (rmax_m[2] - rmin_m[2]);
+        //potentialEnergy = 0.5 * globaltemp * volume / totalP ;
+
+        rho_m = dot(E_m, E_m);
+        potentialEnergy = 0.5 * hr_m[0] * hr_m[1] * hr_m[2] * rho_m.sum();
+
+        auto Pview = P.getView();
+        auto qView = q.getView();
+
+        temp = 0.0;
+
+        Kokkos::parallel_reduce("Kinetic Energy", this->getLocalNum(),
+                                KOKKOS_LAMBDA(const int i, double& valL){
+                                    double myVal = dot(Pview(i), Pview(i)).apply();
+                                    myVal *= -qView(i);
+                                    valL += myVal;
+                                }, Kokkos::Sum<double>(temp));
+
+        temp *= 0.5;
+        globaltemp = 0.0;
+        MPI_Reduce(&temp, &globaltemp, 1, MPI_DOUBLE, MPI_SUM, 0, Ippl::getComm());
+
+        kineticEnergy = globaltemp;
+
+        if (Ippl::Comm->rank() == 0) {
+            std::stringstream fname;
+            fname << "data/Energy_";
+            fname << Ippl::Comm->size();
+            fname << ".csv";
+
+
+            Inform csvout(NULL, fname.str().c_str(), Inform::APPEND);
+            csvout.precision(10);
+            csvout.setf(std::ios::scientific, std::ios::floatfield);
+
+            if(time_m == 0.0) {
+                csvout << "time, Potential energy, Kinetic energy, Total energy" << endl;
+            }
+
+            csvout << time_m << " "
+                   << potentialEnergy << " "
+                   << kineticEnergy << " "
+                   << potentialEnergy + kineticEnergy << endl;
+
+        }
+
+        Ippl::Comm->barrier();
+     }
+
+
      void dumpBumponTail() {
 
         const int nghostE = E_m.getNghost();

alpine/LandauDamping.cpp → alpine/ElectrostaticPIC/LandauDamping.cpp

@@ -1,6 +1,6 @@
 // Landau Damping Test
 //   Usage:
-//     srun ./LandauDamping <nx> <ny> <nz> <Np> <Nt> <stype> <lbthres> <ovfactor> --info 10
+//     srun ./LandauDamping <nx> <ny> <nz> <Np> <Nt> <stype> <lbthres> <ovfactor> <dt> --info 10
 //     nx       = No. cell-centered points in the x-direction
 //     ny       = No. cell-centered points in the y-direction
 //     nz       = No. cell-centered points in the z-direction
@@ -13,6 +13,7 @@
 //                simulations.
 //     ovfactor = Over-allocation factor for the buffers used in the communication. Typical
 //                values are 1.0, 2.0. Value 1.0 means no over-allocation.
+//     dt       = Time stepize
 //     Example:
 //     srun ./LandauDamping 128 128 128 10000 10 FFT 0.01 2.0 --info 10
 //
@@ -177,6 +178,7 @@ int main(int argc, char *argv[]){
 
     const size_type totalP = std::atoll(argv[4]);
     const unsigned int nt     = std::atoi(argv[5]);
+    const double dt = std::atof(argv[9]);;//0.5*dx;
 
     msg << "Landau damping"
         << endl
@@ -200,7 +202,7 @@ int main(int argc, char *argv[]){
 
     // create mesh and layout objects for this problem domain
     Vector_t kw = {0.5, 0.5, 0.5};
-    double alpha = 0.05;
+    double alpha = 0.5;
     Vector_t rmin(0.0);
     Vector_t rmax = 2 * pi / kw ;
     double dx = rmax[0] / nr[0];
@@ -209,7 +211,6 @@ int main(int argc, char *argv[]){
 
     Vector_t hr = {dx, dy, dz};
     Vector_t origin = {rmin[0], rmin[1], rmin[2]};
-    const double dt = 0.5*dx;
 
     const bool isAllPeriodic=true;
     Mesh_t mesh(domain, hr, origin);
@@ -331,7 +332,8 @@ int main(int argc, char *argv[]){
     P->gatherCIC();
 
     IpplTimings::startTimer(dumpDataTimer);
-    P->dumpLandau();
+    P->dumpLandauParticle(totalP);
+    P->dumpEnergy(totalP);
     P->gatherStatistics(totalP);
     //P->dumpLocalDomains(FL, 0);
     IpplTimings::stopTimer(dumpDataTimer);
@@ -390,7 +392,8 @@ int main(int argc, char *argv[]){
 
         P->time_m += dt;
         IpplTimings::startTimer(dumpDataTimer);
-        P->dumpLandau();
+        P->dumpLandauParticle(totalP);
+        P->dumpEnergy(totalP);
         P->gatherStatistics(totalP);
         IpplTimings::stopTimer(dumpDataTimer);
         msg << "Finished time step: " << it+1 << " time: " << P->time_m << endl;