Adds the functionality to output force and force gradient to the output in .h5 file, triggered using Lcheck=7

.github/workflows/build.yml

@@ -19,7 +19,7 @@ jobs:
     - name: Install dependencies
       run: |
         sudo apt-get update
-        sudo apt-get install gfortran mpi-default-bin mpi-default-dev libhdf5-103 libhdf5-dev libfftw3-bin libfftw3-dev libopenblas0-openmp libopenblas-dev 
+        sudo apt-get install gfortran mpi-default-bin mpi-default-dev libhdf5-dev libfftw3-bin libfftw3-dev libopenblas0-openmp libopenblas-dev 
         pip3 install --user numpy f90nml scikit-build scipy h5py matplotlib
         cd ${{ github.workspace }}/Utilities/pythontools
         pip3 install -v .

.github/workflows/build_cmake.yml

@@ -15,7 +15,7 @@ jobs:
     - name: Install prereqs
       run: |
         sudo apt-get update
-        sudo apt-get install gfortran mpi-default-bin mpi-default-dev libhdf5-103 libhdf5-dev libfftw3-bin libfftw3-dev libopenblas0-openmp libopenblas-dev 
+        sudo apt-get install gfortran mpi-default-bin mpi-default-dev libhdf5-dev libfftw3-bin libfftw3-dev libopenblas0-openmp libopenblas-dev 
         pip3 install --upgrade pip
         pip3 install --user ninja cmake scipy
         pip3 install --user numpy f90nml scikit-build scipy h5py matplotlib

.github/workflows/python_wrapper.yml

@@ -8,6 +8,11 @@ jobs:
     steps:
     - uses: actions/checkout@v4
 
+    - name: Set up Python 3.11
+        uses: actions/setup-python@v4
+        with:
+          python-version: "3.11"
+
     # Python3 should be pre-installed on 'ubuntu-latest'
     - name: Python version info
       run: |

CMakeLists.txt

@@ -80,36 +80,49 @@ find_package (LAPACK REQUIRED)
 # Find FFTW3.
 # Intel oneAPI has FFT3 available as part of MKL. Just linking against MKL gives FFTW functionality.
 if (NOT ${BLA_VENDOR} MATCHES "^Intel")
-    # For standard FFTW library
-    # We are using findfftw module from github/egpbos. 
-    # PPPL provides fftwconfig.cmake. We could try that first in future
-
-    # To guide cmake to the location of FFT, set FFTW_ROOT to the path of FFTW
-    # Ex: For PPPL gcc suite of compilers, 
-    # We are using FFTW find module from github.com/egpbos
-    # 
-    configure_file(downloadFindFFTW.cmake.in findFFTW-download/CMakeLists.txt)
-    execute_process(COMMAND ${CMAKE_COMMAND} -G "${CMAKE_GENERATOR}" .
-        RESULT_VARIABLE result
-        WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/findFFTW-download )
-    if(result)
-        message(FATAL_ERROR "CMake step for findFFTW failed: ${result}")
-    else()
-        message("CMake step for findFFTW completed (${result}).")
-    endif()
-    execute_process(COMMAND ${CMAKE_COMMAND} --build .
-        RESULT_VARIABLE result
-        WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/findFFTW-download )
-    if(result)
-        message(FATAL_ERROR "Build step for findFFTW failed: ${result}")
-    endif()
+  # Hack to find the proper .so looking first in conda, then in other 
+  set(FFTW3_TEST_PATHS) # Initialize empty list
+
+  if(DEFINED CONDA_PREFIX)
+    list(APPEND FFTW3_TEST_PATHS "${CONDA_PREFIX}/lib") # Add Conda first
+  endif()
+
+  list(APPEND FFTW3_TEST_PATHS
+    /usr/lib
+    /usr/local/lib
+    ${FFTW3_ROOT_DIR}/lib
+  )
+
+  # Find the library using the constructed list of paths. Should find the path to the .so
+  find_library(FFTW3_LIBRARY_TRUE
+    NAMES fftw3
+    PATHS ${FFTW3_TEST_PATHS}
+  )
+
+
+  set(FFTW_INC_TEST_PATHS) # Initialize an empty list for include paths
+
+  foreach(path ${FFTW3_TEST_PATHS})
+    string(REGEX REPLACE "/lib$" "/include" INCLUDE_PATH "${path}")
+    list(APPEND FFTW_INC_TEST_PATHS "${INCLUDE_PATH}")
+  endforeach()
+
+  find_path(FFTW3_INCLUDE_DIRS # Use FFTW3_INCLUDE_DIRS as in the normal find_package call.
+    NAMES fftw3.f90 fftw3.h # Check for both Fortran and C headers
+    PATHS ${FFTW_INC_TEST_PATHS}
+    REQUIRED
+  )
+
+  if(NOT FFTW3_LIBRARY_TRUE)
+      message(FATAL_ERROR "Could not find FFTW3 library with find_library.")
+  endif()
+
+  message(STATUS "found fftw3 library with path: ${FFTW3_LIBRARY_TRUE}")
+  message(STATUS "found fftw3 include dir with path: ${FFTW3_INCLUDE_DIRS}")
 
-    set(findFFTW_DIR ${CMAKE_CURRENT_BINARY_DIR}/findFFTW-src)
-    set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${findFFTW_DIR}")
 
-    find_package(FFTW REQUIRED)
 else()   # Not required for OneAPI, but for older  versions
-    set(FFTW_INCLUDE_DIRS $ENV{MKLROOT}/include/fftw)
+    set(FFTW3_INCLUDE_DIRS $ENV{MKLROOT}/include/fftw)
 endif()
 
 # Find HDF5

src/CMakeLists.txt

@@ -262,11 +262,11 @@ target_include_directories(spec PUBLIC ${HDF5_C_INCLUDE_DIRS} ${HDF5_Fortran_INC
 # MKL could be used for all the three
 target_link_libraries(spec 
     PUBLIC 
-        ${FFTW_LIBRARIES} ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES}
+    ${FFTW3_LIBRARY_TRUE} ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES}
 )
 target_include_directories(spec 
     PUBLIC 
-        ${FFTW_INCLUDE_DIRS} ${LAPACK_INCLUDE_DIRS} ${BLAS_INCLUDE_DIRS}
+        ${FFTW3_INCLUDE_DIRS} ${LAPACK_INCLUDE_DIRS} ${BLAS_INCLUDE_DIRS}
 )
 
 #if(SKBUILD)

src/global.f90

@@ -741,6 +741,8 @@ module allglobal
   LOGICAL              :: Lhessian3Dallocated !< flag to indicate that 2D Hessian matrix is allocated (?)
   REAL,    allocatable :: hessian3D(:,:) !< Hessian 3D
   REAL,    allocatable :: dessian3D(:,:) !< derivative Hessian 3D
+
+  REAL,    allocatable :: force_final(:) !< Final force on the interfaces [inface*mode]
 !> @}
 
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!

src/inputlist.f90

@@ -535,6 +535,7 @@ module inputlist
                                                  !<      <li> must set \c Freebound=1, \c Lfindzero>0, \c mfreeits!=0 </li>
                                                  !<      <li> \c xdiagno must be executed manually </li>
                                                  !<      </ul> </li>
+                                                 !< <li> if \c Lcheck = 7, calculates the force and force_grad, outputs them to the .h5 file </li>
                                                  !< </ul>
   LOGICAL      :: Ltiming          =  .false.    !< to check timing
   LOGICAL      :: Ltransform       = .false.     !< to evaluate iota and straight field line coordinates

src/sphdf5.f90

@@ -1038,7 +1038,8 @@ subroutine hdfint
                         lmns, &
                         TT, &
                         beltramierror, &
-                        IPDt, dlambdaout, lmns
+                        IPDt, dlambdaout, lmns, &
+                        force_final, hessian, NGdof
 
   LOCALS
 
@@ -1165,6 +1166,11 @@ subroutine hdfint
 !> <li> All quantities marked as real should be treated as double precision. </li>
 !> </ul>
 
+  if( Lcheck.eq.7 ) then
+    HWRITERV(grpOutput, NGdof+1, force_final, force_final(0:NGdof))
+    HWRITERA( grpOutput, NGdof, NGdof, force_final_grad, hessian(1:NGdof,1:NGdof) )
+  end if
+
   HCLOSEGRP( grpOutput )
 
  endif ! myid.eq.0

src/xspech.f90

@@ -317,7 +317,9 @@ subroutine spec
                         first_free_bound, &
                         dMA, dMB, dMD, dMG, MBpsi, solution, IPDt, &
                         version, &
-                        MPI_COMM_SPEC
+                        MPI_COMM_SPEC, &
+                        force_final, Lhessianallocated, LocalConstraint, hessian, dBBdmp, dFFdRZ, dmupfdx, &
+                        dRodR, dRodZ, dZodR, dZodZ, dessian, LGdof
 
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
 
@@ -483,18 +485,37 @@ subroutine spec
 
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
 
-  SALLOCATE( gradient, (0:NGdof), zero )
+  if(Lcheck.eq.7) then
+    SALLOCATE( force_final, (0:NGdof), zero )
+    SALLOCATE( hessian, (1:NGdof,1:NGdof), zero )
+    SALLOCATE( dessian, (1:NGdof,1:LGdof), zero )
+    SALLOCATE( dFFdRZ, (1:LGdof,0:1,1:LGdof,0:1,1:Mvol), zero )
+    SALLOCATE( dBBdmp, (1:LGdof,1:Mvol,0:1,1:2), zero )
+    if( LocalConstraint ) then
+      SALLOCATE( dmupfdx, (1:Mvol,    1:1,1:2,1:LGdof,0:1), zero )
+    else
+      SALLOCATE( dmupfdx, (1:Mvol, 1:Mvol-1,1:2,1:LGdof,1), zero ) ! TODO change the format to put vvol in last index position...
+    endif
+    Lhessianallocated = .true.
+    SALLOCATE( dRodR, (1:Ntz,0:3,1:mn), zero ) ! calculated in rzaxis; 19 Sep 16;
+    SALLOCATE( dRodZ, (1:Ntz,0:3,1:mn), zero )
+    SALLOCATE( dZodR, (1:Ntz,0:3,1:mn), zero )
+    SALLOCATE( dZodZ, (1:Ntz,0:3,1:mn), zero )
+
+    LComputeDerivatives = .true.
 
-  lastcpu = GETTIME
+    LComputeAxis = .true.
+    WCALL( xspech, dforce, ( NGdof, position(0:NGdof), force_final(0:NGdof), LComputeDerivatives, LComputeAxis) )
 
-  LComputeDerivatives = .false.
-  LComputeAxis = .true.
-! vvol = Mvol ; ideriv = 0 ; ii = 1
-! write(ounit,'("xspech : ", 10x ," : sum(Ate(",i3,",",i2,",",i2,")%s) =",99es23.15)') vvol, ideriv, ii, sum(Ate(vvol,ideriv,ii)%s(0:Lrad(vvol)))
+  else
+    SALLOCATE( force_final, (0:NGdof), zero )
 
-  WCALL( xspech, dforce, ( NGdof, position(0:NGdof), gradient(0:NGdof), LComputeDerivatives, LComputeAxis) ) ! (re-)calculate Beltrami fields;
+    LComputeDerivatives = .false.
+    LComputeAxis = .true.
 
-  DALLOCATE(gradient)
+    WCALL( xspech, dforce, ( NGdof, position(0:NGdof), force_final(0:NGdof), LComputeDerivatives, LComputeAxis) )
+  end if
+
 
 #ifdef DEBUG
   do vvol = 1, Mvol-1