Feat: FFT-based self-gravity solver

.gitmodules

@@ -4,3 +4,6 @@
 [submodule "reference"]
 	path = reference
 	url = https://github.com/idefix-code/reference
+[submodule "src/kokkos-fft"]
+	path = src/kokkos-fft
+	url = https://github.com/kokkos/kokkos-fft.git

CMakeLists.txt

@@ -16,6 +16,7 @@ option(Idefix_DEBUG "Enable Idefix debug features (makes the code very slow)" OF
 option(Idefix_RUNTIME_CHECKS "Enable runtime sanity checks" OFF)
 option(Idefix_WERROR "Treat compiler warnings as errors" OFF)
 option(Idefix_PYTHON "Enable python bindings (requires pybind11)" OFF)
+option(Idefix_FFT "Enable FFT (requires dedicated FFT library)" OFF)
 set(Idefix_PROBLEM_DIR "${CMAKE_BINARY_DIR}" CACHE STRING "Problem directory to build for.")
 set(Idefix_CXX_FLAGS "" CACHE STRING "Additional compiler/linker flag")
 set(Idefix_DEFS "definitions.hpp" CACHE FILEPATH "Problem definition header file")
@@ -54,6 +55,12 @@ endif()
 add_subdirectory(src/kokkos ${CMAKE_BINARY_DIR}/build/kokkos)
 include_directories(${Kokkos_INCLUDE_DIRS_RET})
 
+# Add FFT library if requested
+if(Idefix_FFT)
+  add_compile_definitions("WITH_FFT")
+  add_subdirectory(src/kokkos-fft)
+endif()
+
 # Add Idefix CXX Flags
 add_compile_options(${Idefix_CXX_FLAGS})
 
@@ -248,12 +255,17 @@ target_include_directories(idefix PUBLIC
                            src/rkl
                            src/gravity
                            src/utils
+                           src/utils/fft
                            src/utils/iterativesolver
                            src/mpi
                            src
                            )
 
 target_link_libraries(idefix Kokkos::kokkos)
+if(Idefix_FFT)
+  target_link_libraries(idefix KokkosFFT::fft)
+  add_subdirectory(src/utils/fft)
+endif()
 
 message(STATUS "Idefix final configuration")
 if(Idefix_EVOLVE_VECTOR_POTENTIAL)
@@ -263,6 +275,7 @@ else()
 endif()
 message(STATUS "    MPI:  ${Idefix_MPI}")
 message(STATUS "    HDF5: ${Idefix_HDF5}")
+message(STATUS "    FFT:  ${Idefix_FFT}")
 message(STATUS "    Python: ${Idefix_PYTHON}")
 message(STATUS "    Reconstruction: ${Idefix_RECONSTRUCTION}")
 message(STATUS "    Precision: ${Idefix_PRECISION}")

CPPLINT.cfg

@@ -1,7 +1,7 @@
 # Don't search for additional CPPLINT.cfg in parent directories.
 set noparent
 headers=hpp
-linelength=100
+linelength=140
 # Don't use 'SRC_' as the cpp header guard prefix
 root=./src/
 extensions=hpp,cpp

src/gravity/CMakeLists.txt

@@ -5,4 +5,8 @@ target_sources(idefix
   PUBLIC ${CMAKE_CURRENT_LIST_DIR}/laplacian.hpp
   PUBLIC ${CMAKE_CURRENT_LIST_DIR}/selfGravity.hpp
   PUBLIC ${CMAKE_CURRENT_LIST_DIR}/selfGravity.cpp
+  PUBLIC ${CMAKE_CURRENT_LIST_DIR}/selfGravityIterative.hpp
+  PUBLIC ${CMAKE_CURRENT_LIST_DIR}/selfGravityIterative.cpp
+  PUBLIC ${CMAKE_CURRENT_LIST_DIR}/selfGravityFFT.hpp
+  PUBLIC ${CMAKE_CURRENT_LIST_DIR}/selfGravityFFT.cpp
   )

src/gravity/gravity.cpp

@@ -98,7 +98,10 @@ Gravity::Gravity(Input &input, DataBlock *datain) {
 
   // Check SelfGravity object
   if(haveSelfGravityPotential) {
-    selfGravity.Init(input, this->data);
+    if(!haveInitialisedSelfGravity) {
+      selfGravity = SelfGravity::Create(input, this->data);
+    }
+    selfGravity->Init(input, this->data);
     haveInitialisedSelfGravity = true;
   }
 
@@ -125,7 +128,7 @@ void Gravity::ShowConfig() {
     }
     if(haveSelfGravityPotential) {
       idfx::cout << "Gravity: self-gravity ENABLED." << std::endl;
-      selfGravity.ShowConfig();
+      selfGravity->ShowConfig();
     }
     if(havePlanetsPotential) {
       idfx::cout << "Gravity: planet(s) potential ENABLED." << std::endl;
@@ -165,10 +168,10 @@ void Gravity::ComputeGravity(int stepNumber) {
     }
     if(haveSelfGravityPotential) {
       // Solving Poisson for the current gas density distribution
-      if(stepNumber % selfGravity.skipSelfGravity == 0) selfGravity.SolvePoisson();
+      if(stepNumber % selfGravity->skipSelfGravity == 0) selfGravity->SolvePoisson();
 
       // Adding gas self-gravity contribution to global gravity potential
-      selfGravity.AddSelfGravityPotential(phiP);
+      selfGravity->AddSelfGravityPotential(phiP);
     }
   }
   if(haveBodyForce) {

src/gravity/gravity.hpp

@@ -8,6 +8,7 @@
 #ifndef GRAVITY_GRAVITY_HPP_
 #define GRAVITY_GRAVITY_HPP_
 
+#include <memory>
 #include "idefix.hpp"
 #include "input.hpp"
 #include "selfGravity.hpp"
@@ -52,7 +53,7 @@ class Gravity {
   IdefixArray4D<real> bodyForceVector;
 
   // Self gravity
-  SelfGravity selfGravity;
+  std::unique_ptr<SelfGravity> selfGravity;
 
   // JM : moved in public class to handle changing centralMass during computation
   real centralMass{1.0};                    ///< central mass parameter when central mass potential