ThinCurr: Improvements to mesh/model handling and HODLR compression (OpenFUSIONToolkit#332)

- Support resistivity and coil definitions within the mesh file itself, instead of relying on a separate XML file
   * Separate holes and jumpers in mesh definition for improved clarity
   * Holes are now checked when a model is loaded, using the automated script, to avoid accidentally loading malformed meshes
   * Define coils generically in file to enable runtime selection of Icoil vs Vcoil
   * Add a new `OFT_ThinCurr_mesh.py` script that can be used set the above properties and modify/merge ThinCurr mesh files
- Improve HODLR performance and robustness
   * HODLR tolerances can now be set from within the Python API using `set_hodlr_options()`
   * HODLR tolerances now support mesh-relative setting using the average cell size (default)
   * Improve data-locality in parallelization of HODLR operators
- Improve the performance of the hole-finding script `OFT_ThinCurr_holes.py`
- Add support for variable timestep size in `ThinCurr.run_td()` (requires `direct=False`)
- Increase magnetic field operator finite-difference step size and make user-configurable
- Update Icoil source handling in `ThinCurr.run_td()` to be more consistent across Vcoil/Icoil usage
- Reported iterations/time in `ThinCurr.run_td()` are now the average across the printing frequency
- Standardize mesh reading and writing across Python API
- Add documentation pages for all OFT python scripts (eg. `OFT_convert_Cubit.py`)
- Remove legacy ThinCurr fortran executables

This PR makes the following **changes to default settings/behavior**:

- The default tolerances for HODLR operators are now mesh size based (10E-5; relative and scaling with RMS edge size)
- The default step size for mesh B-field calculations is larger and mesh based (1E-2; relative and scaling with RMS edge size)

This PR makes the following **backwards-compatible changes to the Python API**:

- A new `no_verify_holes` argument to `ThinCurr.setup_model()` to optionally skip the newly added hole verification on model load
- A new `B_dx` argument to `ThinCurr.compute_Bmat()` to set the finite difference step size
- A new `times` argument to `ThinCurr.run_td()` to enable specifying variable timesteps in a simulation
- A new `ThinCurr.set_hodlr_options()` method to enable adjusting HODLR tolerances and other parameters
- A new `ThinCurr.update_coil_types()` method to enable setting coil types (`Icoil` or `Vcoil`) at runtime

--------------
Co-authored by: Claude Opus 4.8

Author: hansec

src/CMakeLists.txt

@@ -29,7 +29,7 @@ option(OFT_PACKAGE_PYTHON "Setup Python packaging?" ON)
 option(OFT_COVERAGE "Build with code coverage information?" OFF)
 option(OFT_EXECSTACK "Allow executable stack in shared libraries/executables?" OFF)
 option(OFT_TOKAMAKER_LEGACY "UNSUPPORTED: Build legacy TokaMaker executables?" OFF)
-option(OFT_THINCURR_LEGACY "DEPRECATED: Build legacy ThinCurr executables?" OFF)
+option(OFT_THINCURR_LEGACY "UNSUPPORTED: Build legacy ThinCurr executables?" OFF)
 option(OFT_PY_KERNEL "Python kernel name to use for testing examples" "python3")
 set(OFT_MPI_PLEN "4" CACHE STRING "Size of strings for MPI processor index")
 set(OFT_PATH_SLEN "200" CACHE STRING "Size of path strings")
@@ -59,15 +59,13 @@ if(OFT_PACKAGE_BUILD AND OFT_PACKAGE_NIGHTLY)
   set( OFT_VER_SUFFIX "-${GIT_REV_HASH}" CACHE INTERNAL "Version suffix" )
 endif()
 
-if(OFT_BUILD_DOCS OR OFT_PACKAGE_BUILD OR OFT_DEBUG_STACK)
-  if(OFT_PACKAGE_BUILD AND OFT_PACKAGE_PYTHON)
-    if(APPLE AND (NOT DEFINED CMAKE_OSX_DEPLOYMENT_TARGET))
-      message( FATAL_ERROR "Python package build on macOS requires the MACOSX_DEPLOYMENT_TARGET environment variable to be set" )
-    endif()
-    set(Python3_FIND_VIRTUALENV ONLY)
+if(OFT_PACKAGE_BUILD AND OFT_PACKAGE_PYTHON)
+  if(APPLE AND (NOT DEFINED CMAKE_OSX_DEPLOYMENT_TARGET))
+    message( FATAL_ERROR "Python package build on macOS requires the MACOSX_DEPLOYMENT_TARGET environment variable to be set" )
   endif()
-  find_package(Python3 COMPONENTS Interpreter REQUIRED)
+  set(Python3_FIND_VIRTUALENV ONLY)
 endif()
+find_package(Python3 COMPONENTS Interpreter REQUIRED)
 
 ########################
 # Build documentation
@@ -195,7 +193,7 @@ elseif("${CMAKE_Fortran_COMPILER_ID}" MATCHES "GNU")
   set( CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_C_FLAGS_RELEASE} -Werror=line-truncation -Wconversion" )
   if(NOT OFT_PACKAGE_BUILD)
     if(OFT_DEBUG_CHECK)
-      set( CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} -fcheck=all" )
+      set( CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} -fcheck=all,no-array-temps" )
     elseif(OFT_DEBUG_SANITIZER)
       set( CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} -fsanitize=address" )
     endif()
@@ -430,6 +428,9 @@ endif()
 if( OFT_TOKAMAKER_LEGACY )
   message( FATAL_ERROR "Legacy TokaMaker functionality is no longer available." )
 endif()
+if( OFT_THINCURR_LEGACY )
+  message( FATAL_ERROR "Legacy ThinCurr functionality is no longer available." )
+endif()
 # Include files
 add_subdirectory( include )
 include_directories( ${CMAKE_BINARY_DIR}/include )
@@ -533,9 +534,7 @@ if( OFT_BUILD_PYTHON )
     set(OFT_SHARED_LIBS ${ARGV0} ${OFT_SHARED_LIBS} CACHE INTERNAL "package_shared_libs")
   endfunction(oft_add_shared)
   # Library directories
-  add_subdirectory( python/OpenFUSIONToolkit )
-  add_subdirectory( python/wrappers )
-  add_subdirectory( python/scripts )
+  add_subdirectory( python )
 endif()
 
 ########################

src/base/oft_io.F90

@@ -111,6 +111,7 @@ MODULE oft_io
 !------------------------------------------------------------------------------
 INTERFACE hdf5_read
   MODULE PROCEDURE hdf5_read_string
+  MODULE PROCEDURE hdf5_read_strings
   MODULE PROCEDURE hdf5_read_scalar_r8
   MODULE PROCEDURE hdf5_read_scalar_i4
   MODULE PROCEDURE hdf5_read_scalar_l
@@ -589,7 +590,8 @@ function hdf5_field_exist(filepath,path) result(exists)
 character(LEN=*), intent(in) :: path !< Path of field in file
 integer :: access_flag,error,subpath
 integer(HID_T) :: file_id,dset_id
-logical :: exists
+logical :: exists,exists_lc
+character(LEN=:), allocatable :: path_lc
 exists=oft_file_exist(filepath)
 IF(.NOT.exists)RETURN
 DEBUG_STACK_PUSH
@@ -604,15 +606,27 @@ function hdf5_field_exist(filepath,path) result(exists)
   IF(subpath<LEN_TRIM(path).AND.path(subpath:subpath)/="/")CYCLE
   CALL h5lexists_f(file_id, "/"//path(1:subpath), exists, error)
 END DO
+!---Check for lowercase variant of path (HDF5 is case sensitive)
+IF(.NOT.exists)THEN
+  path_lc=TRIM(path)
+  CALL string_to_lower(path_lc)
+  exists_lc=.TRUE.
+  DO subpath=1,LEN_TRIM(path)
+    IF(.NOT.exists_lc)EXIT
+    IF(subpath<LEN_TRIM(path).AND.path(subpath:subpath)/="/")CYCLE
+    CALL h5lexists_f(file_id, "/"//path_lc(1:subpath), exists_lc, error)
+  END DO
+  IF(exists_lc)THEN
+    CALL oft_warn('Case variant of path exists in HDF5 file: '//TRIM(path))
+  END IF
+END IF
 !---Close file and finalize HDF5
 call h5fclose_f(file_id, error)
 call h5close_f(error)
 DEBUG_STACK_POP
 end function hdf5_field_exist
 !------------------------------------------------------------------------------
 !> Test for exitence of a field in a HDF5 file
-!!
-!! @result Logical flag indicating existence of field and file
 !------------------------------------------------------------------------------
 subroutine hdf5_field_get_sizes(filepath,path,ndims,dim_sizes)
 character(LEN=*), intent(in) :: filepath !< Path to file
@@ -634,15 +648,32 @@ subroutine hdf5_field_get_sizes(filepath,path,ndims,dim_sizes)
   CALL h5dopen_f(file_id, "/"//TRIM(path), dset_id, error)
   IF(error==0)THEN
     CALL h5dget_space_f(dset_id, dspace_id, error)
-    CALL h5sget_simple_extent_ndims_f(dspace_id, ndims, error)
-    ALLOCATE(dim_sizes(ndims),tmp_sizes(ndims),maxdims(ndims))
-    CALL h5sget_simple_extent_dims_f(dspace_id, tmp_sizes, maxdims, error)
-    dim_sizes=INT(tmp_sizes,4)
-    DEALLOCATE(tmp_sizes,maxdims)
-    !---Close dataspace/set
-    call h5sclose_f(dspace_id, error)
-    call h5dclose_f(dset_id, error)
+    IF(error==0)THEN
+      CALL h5sget_simple_extent_ndims_f(dspace_id, ndims, error)
+      IF(error==0)THEN
+        ALLOCATE(dim_sizes(ndims),tmp_sizes(ndims),maxdims(ndims))
+        CALL h5sget_simple_extent_dims_f(dspace_id, tmp_sizes, maxdims, error)
+        IF(error>=0)THEN
+          dim_sizes=INT(tmp_sizes,4)
+        ELSE
+          ndims=-6
+          DEALLOCATE(dim_sizes)
+        END IF
+        DEALLOCATE(tmp_sizes,maxdims)
+        !---Close dataspace/set
+        call h5sclose_f(dspace_id, error)
+      ELSE
+        ndims=-5
+      END IF
+      call h5dclose_f(dset_id, error)
+    ELSE
+      ndims=-4
+    END IF
+  ELSE
+    ndims=-3
   END IF
+ELSE
+  ndims=-2
 END IF
 !---Close file and finalize HDF5
 call h5fclose_f(file_id, error)
@@ -1204,10 +1235,55 @@ subroutine hdf5_read_string(string,filename,path,success)
 character(LEN=*), intent(in) :: filename !< Path to file
 character(LEN=*), intent(in) :: path !< Variable path in file
 logical, optional, intent(out) :: success !< Successful read?
+integer(i4) :: i
+character(LEN=:), allocatable :: strings(:)
+IF(ALLOCATED(string))DEALLOCATE(string)
+CALL hdf5_read_strings(strings,filename,path,success)
+IF(PRESENT(success))THEN
+  IF(.NOT.success)RETURN
+END IF
+IF(.NOT.ALLOCATED(strings))THEN
+  IF(PRESENT(success))THEN
+    success=.FALSE.
+    RETURN
+  ELSE
+    CALL oft_abort('Error reading string from file: '//TRIM(filename)//' path: '//TRIM(path), &
+      'hdf5_read_string',__FILE__)
+  END IF
+END IF
+IF(SIZE(strings)/=1)THEN
+  !---Handle strings stored as a character array
+  IF(LEN(strings(1))==1)THEN
+    ALLOCATE(character(LEN=SIZE(strings)) :: string)
+    DO i=1,SIZE(strings)
+      string(i:i)=strings(i) ! Concatenate single-characters into one string
+    END DO
+    DEALLOCATE(strings)
+    RETURN
+  END IF
+  IF(PRESENT(success))THEN
+    success=.FALSE.
+    RETURN
+  ELSE
+    CALL oft_abort('Error reading string from file: '//TRIM(filename)//' path: '//TRIM(path), &
+      'hdf5_read_string',__FILE__)
+  END IF
+END IF
+string = strings(1)
+DEALLOCATE(strings)
+end subroutine hdf5_read_string
+!------------------------------------------------------------------------------
+!> String implementation of \ref oft_io::hdf5_read
+!------------------------------------------------------------------------------
+subroutine hdf5_read_strings(strings,filename,path,success)
+character(LEN=:), allocatable, intent(inout) :: strings(:) !< String to read from file
+character(LEN=*), intent(in) :: filename !< Path to file
+character(LEN=*), intent(in) :: path !< Variable path in file
+logical, optional, intent(out) :: success !< Successful read?
 integer(i4) :: error
 integer(i4), parameter :: one=1,zero=0
-INTEGER(HSIZE_T) :: space_count
-integer(HID_T) :: file_id,dset_id,dspace_id
+INTEGER(HSIZE_T) :: space_count,string_size
+integer(HID_T) :: file_id,dset_id,dspace_id,type_id
 INTEGER(HSIZE_T), DIMENSION(1) :: dims
 ! DEBUG_STACK_PUSH
 IF(PRESENT(success))THEN
@@ -1217,23 +1293,30 @@ subroutine hdf5_read_string(string,filename,path,success)
 !---Initialize HDF5 and open file
 call h5open_f(error)
 call h5fopen_f(TRIM(filename), H5F_ACC_RDONLY_F, file_id, error)
-IF(error/=0)GOTO 103
+IF(error/=0)GOTO 104
 CALL h5dopen_f(file_id, "/"//TRIM(path), dset_id, error)
+IF(error/=0)GOTO 103
+!---Get string size
+CALL H5Dget_type_f(dset_id, type_id, error)
 IF(error/=0)GOTO 102
-!---
+CALL H5Tget_size_f(type_id, string_size, error)
+IF(error/=0)GOTO 101
+!---Get dataset size (should be 1)
 CALL h5dget_space_f(dset_id, dspace_id, error)
 IF(error/=0)GOTO 101
 CALL h5sget_simple_extent_npoints_f(dspace_id, space_count, error)
 IF(error/=0)GOTO 100
 dims = space_count
-ALLOCATE(CHARACTER(LEN=dims(1)) :: string)
-call h5dread_f(dset_id, H5T_NATIVE_CHARACTER, string, dims, error)
+IF(ALLOCATED(strings))DEALLOCATE(strings)
+ALLOCATE(CHARACTER(LEN=string_size) :: strings(dims(1)))
+call h5dread_f(dset_id, type_id, strings, dims, error)
 IF(error/=0)THEN
-  DEALLOCATE(string)
+  DEALLOCATE(strings)
   GOTO 100
 END IF
 !---Close and finalize HDF5
 CALL h5sclose_f(dspace_id, error)
+CALL h5tclose_f(type_id, error)
 call h5dclose_f(dset_id, error)
 call h5fclose_f(file_id, error)
 call h5close_f(error)
@@ -1244,11 +1327,12 @@ subroutine hdf5_read_string(string,filename,path,success)
 END IF
 RETURN
 100 CALL h5sclose_f(dspace_id, error)
-101 CALL h5dclose_f(dset_id, error)
-102 CALL h5fclose_f(file_id, error)
-103 CALL h5close_f(error)
+101 CALL h5tclose_f(type_id, error)
+102 CALL h5dclose_f(dset_id, error)
+103 CALL h5fclose_f(file_id, error)
+104 CALL h5close_f(error)
 IF(PRESENT(success))CALL h5eset_auto_f(one, error)
-end subroutine hdf5_read_string
+end subroutine hdf5_read_strings
 !------------------------------------------------------------------------------
 !> real(r8) scalar implementation of \ref oft_io::hdf5_read
 !------------------------------------------------------------------------------

src/base/oft_local.F90

@@ -56,12 +56,13 @@ END FUNCTION oft_sleep
 !------------------------------------------------------------------------------
 !> Simple in-memory hashing function for dataset checksumming
 !------------------------------------------------------------------------------
-  FUNCTION oft_simple_hash(key,length)  BIND(C)
+  FUNCTION oft_simple_hash_c(key,length,starting_value)  BIND(C)
     IMPORT c_int, c_long, c_ptr
-    INTEGER(c_int) :: oft_simple_hash !< Hash of data
+    INTEGER(c_int) :: oft_simple_hash_c !< Hash of data
     TYPE(c_ptr), VALUE, INTENT(in) :: key !< Location of data
     INTEGER(c_long), VALUE, INTENT(in) :: length !< Length of data to hash in bytes
-  END FUNCTION oft_simple_hash
+    INTEGER(c_int), VALUE, INTENT(in) :: starting_value !< Starting value for the hash
+  END FUNCTION oft_simple_hash_c
 END INTERFACE
 !------------------------------------------------------------------------------
 !> One dimensional integer set
@@ -116,6 +117,7 @@ END FUNCTION oft_2d_func
 !---------------------------------------------------------------------------------
 TYPE :: oft_timer
   INTEGER(i8) :: count = 0 !< Integer value of system clock at last call
+  INTEGER(i8) :: crate = 0 !< Integer value of system clock rate
 CONTAINS
   !> Start or reset timer
   procedure :: tick => oft_timer_start
@@ -128,6 +130,20 @@ END FUNCTION oft_2d_func
 END TYPE oft_timer
 PRIVATE oft_timer_start, oft_timer_elapsed, oft_timer_intelapsed, oft_timer_timeout
 CONTAINS
+!------------------------------------------------------------------------------
+!> Simple in-memory hashing function for dataset checksumming
+!------------------------------------------------------------------------------
+FUNCTION oft_simple_hash(key,length,starting_value) RESULT(data_hash)
+TYPE(c_ptr), INTENT(in) :: key !< Location of data
+INTEGER(c_long), INTENT(in) :: length !< Length of data to hash in bytes
+INTEGER(c_int), OPTIONAL, INTENT(in) :: starting_value !< Starting value for the hash
+INTEGER(c_int) :: data_hash !< Hash of data
+IF(PRESENT(starting_value))THEN
+  data_hash=oft_simple_hash_c(key,length,starting_value)
+ELSE
+  data_hash=oft_simple_hash_c(key,length,0)
+END IF
+END FUNCTION oft_simple_hash
 !---------------------------------------------------------------------------------
 !> Returns the corresponding lowercase letter, if `c` is an uppercase
 !! ASCII character, otherwise `c` itself.
@@ -141,8 +157,8 @@ ELEMENTAL FUNCTION char_to_lower(c) result(t)
 INTEGER(i4), PARAMETER :: BA=iachar('A')
 INTEGER(i4), PARAMETER :: BZ=iachar('Z')
 INTEGER(i4) :: k
-k = ichar(c) 
-IF(k>=BA.and.k<=BZ)k = k + wp 
+k = ichar(c)
+IF(k>=BA.and.k<=BZ)k = k + wp
 t = char(k)
 END FUNCTION char_to_lower
 !---------------------------------------------------------------------------------
@@ -168,8 +184,8 @@ ELEMENTAL FUNCTION char_to_upper(c) result(t)
 INTEGER(i4), PARAMETER :: BA=iachar('a')
 INTEGER(i4), PARAMETER :: BZ=iachar('z')
 INTEGER(i4) :: k
-k = ichar(c) 
-IF(k>=BA.and.k<=BZ)k = k + wp 
+k = ichar(c)
+IF(k>=BA.and.k<=BZ)k = k + wp
 t = char(k)
 END FUNCTION char_to_upper
 !---------------------------------------------------------------------------------
@@ -187,7 +203,8 @@ END SUBROUTINE string_to_upper
 !---------------------------------------------------------------------------------
 SUBROUTINE oft_timer_start(self)
 CLASS(oft_timer), INTENT(inout) :: self !< Calling timer class
-self%count=oft_time_i8()
+INTEGER(i8) :: cmax
+CALL system_clock(self%count,self%crate,cmax)
 END SUBROUTINE oft_timer_start
 !---------------------------------------------------------------------------------
 !> Set elapsed time since last tick/tock
@@ -204,9 +221,6 @@ FUNCTION oft_timer_elapsed(self) RESULT(time)
 END FUNCTION oft_timer_elapsed
 !---------------------------------------------------------------------------------
 !> Get elapsed time since last tick/tock in integer counts
-!!
-!! @param[in,out] self Calling timer class
-!! @return Number of integer counts since last tick/tock
 !---------------------------------------------------------------------------------
 FUNCTION oft_timer_intelapsed(self) result(dt)
 CLASS(oft_timer), INTENT(inout) :: self !< Calling timer class

src/base/oft_local_c.c

@@ -83,9 +83,9 @@ void oft_set_signal_handlers()
 }
 
 // Jenkins's one_at_a_time hash
-uint32_t oft_simple_hash(const uint8_t* key, long length) {
+uint32_t oft_simple_hash_c(const uint8_t* key, long length, uint32_t starting_value) {
   long i = 0;
-  uint32_t hash = 0;
+  uint32_t hash = starting_value;
   while (i != length) {
     hash += key[i++];
     hash += hash << 10;

src/bin/CMakeLists.txt

@@ -17,9 +17,6 @@ set( THINCURR_EXES
 )
 
 set( BIN_SRCS ${OFT_EXES} ${MARKLIN_EXES} )
-if( OFT_THINCURR_LEGACY )
-  set( BIN_SRCS ${BIN_SRCS} ${THINCURR_EXES} )
-endif()
 
 foreach( file ${BIN_SRCS} )
   oft_add_exe( ${file} )