From 392b12e10c5a1be77a06ea721d9583b2bada1926 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Fri, 27 Mar 2026 11:15:41 +0100 Subject: [PATCH 1/2] remove: strip all GVEC library dependencies and related code --- .gitignore | 1 - CMakeLists.txt | 17 - python/CMakeLists.txt | 3 - requirements.txt | 4 - src/CMakeLists.txt | 19 -- src/boozer_converter.F90 | 4 - src/field.F90 | 24 -- src/field/field_gvec.f90 | 213 ------------- src/field/field_newton.F90 | 13 +- src/field/vmec_field_adapter.f90 | 317 ------------------- src/field/vmec_field_eval.f90 | 1 - src/get_canonical_coordinates.F90 | 12 - test/golden_record/golden_record.sh | 2 - test/test_data/GVEC_elliptok_State_final.dat | 52 --- test/tests/CMakeLists.txt | 74 ----- test/tests/convert_vmec_to_gvec.py | 99 ------ test/tests/export_field_2d.f90 | 119 ++----- test/tests/test_adapter_consistency.f90 | 181 ----------- test/tests/test_canonical_gvec.f90 | 176 ---------- test/tests/test_gvec.f90 | 64 ---- test/tests/test_simple_vmec_gvec.f90 | 144 --------- test/tests/test_vmec_gvec.f90 | 254 --------------- test/tests/test_vmec_gvec_adapter.f90 | 284 ----------------- thirdparty/LICENSE.GVEC | 9 - 24 files changed, 29 insertions(+), 2057 deletions(-) delete mode 100644 src/field/field_gvec.f90 delete mode 100644 src/field/vmec_field_adapter.f90 delete mode 100644 test/test_data/GVEC_elliptok_State_final.dat delete mode 100644 test/tests/convert_vmec_to_gvec.py delete mode 100644 test/tests/test_adapter_consistency.f90 delete mode 100644 test/tests/test_canonical_gvec.f90 delete mode 100644 test/tests/test_gvec.f90 delete mode 100644 test/tests/test_simple_vmec_gvec.f90 delete mode 100644 test/tests/test_vmec_gvec.f90 delete mode 100644 test/tests/test_vmec_gvec_adapter.f90 delete mode 100644 thirdparty/LICENSE.GVEC diff --git a/.gitignore b/.gitignore index d80f9e93..5647313e 100644 --- a/.gitignore +++ b/.gitignore @@ -55,7 +55,6 @@ threed1.* threed1.* timings.txt .aider* -thirdparty/gvec thirdparty/pyplot_module.F90 *.html diff --git a/CMakeLists.txt b/CMakeLists.txt index e2304dcf..17ef1d36 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -73,7 +73,6 @@ endif() option(ENABLE_PYTHON_INTERFACE "Enables the Python-Wrapper." ON) option(SIMPLE_TESTING "Enable testing." ON) option(ENABLE_OPENMP "Enable OpenMP compiler flags." ON) -option(ENABLE_GVEC "Enable GVEC field support (experimental)" OFF) option(ENABLE_COVERAGE "Enable code coverage analysis (Debug/Profile builds only)" OFF) option(SIMPLE_DETERMINISTIC_FP "Disable fast-math for reproducible floating-point" OFF) option(SIMPLE_ENABLE_PYTHON_TOOLS "Enable Python helpers (tests/data generation)" ON) @@ -112,21 +111,6 @@ if(CMAKE_Fortran_COMPILER_ID STREQUAL "NVHPC") endif() endif() -# Conditionally fetch GVEC if enabled -if(ENABLE_GVEC) - message(STATUS "GVEC support enabled - fetching GVEC library") - include(FetchContent) - FetchContent_Declare( - gvec - GIT_REPOSITORY https://gitlab.mpcdf.mpg.de/calbert/gvec.git - GIT_TAG develop - ) - set(GVEC_FOR_SIMPLE TRUE CACHE BOOL "Build minimal GVEC for SIMPLE") - FetchContent_MakeAvailable(gvec) -else() - message(STATUS "GVEC support disabled") -endif() - set(CMAKE_MACOSX_RPATH 1) find_program(NC_CONFIG "nf-config") @@ -248,7 +232,6 @@ message(STATUS "CMake build type: " ${CMAKE_BUILD_TYPE}) include_directories ($ENV{NETCDFF_INCLUDE} ${NFINC}) link_directories ($ENV{NETCDF_LIB} $ENV{NETCDFF_LIB} ${NFLIBS}) -# Set SIMPLE-specific compiler flags after GVEC to avoid affecting it # Compiler-specific flags (not applied to subprojects) if(CMAKE_Fortran_COMPILER_ID STREQUAL "GNU") set(SIMPLE_COMPILE_OPTIONS diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index b8ac740f..5b3d26ca 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -135,9 +135,6 @@ set(GENERATED_MODULE_FILE "${F2PY_MODULE_NAME}.${Python_SOABI}") set(F2PY_WRAPPERS "${CMAKE_CURRENT_BINARY_DIR}/${F2PY_MODULE_C}") set(SIMPLE_F90PP_DEFINES) -if(ENABLE_GVEC) - list(APPEND SIMPLE_F90PP_DEFINES -DGVEC_AVAILABLE) -endif() if(CMAKE_Fortran_COMPILER_ID STREQUAL "NVHPC") list(APPEND SIMPLE_F90PP_DEFINES -DSIMPLE_NVHPC) endif() diff --git a/requirements.txt b/requirements.txt index 633d86a9..77fb5fe1 100644 --- a/requirements.txt +++ b/requirements.txt @@ -17,7 +17,3 @@ netCDF4==1.7.2 matplotlib==3.9.4 scipy==1.15.2 shapely==2.0.5 - -# map2disc is required for map2disc-based chartmap generation via libneo -# Pin to a specific commit for reproducibility. -map2disc @ git+https://gitlab.mpcdf.mpg.de/gvec-group/map2disc.git@14da2f094ab2505a639b9cd04af57f5c0dbfaa54 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index e3668990..24c28c40 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -99,14 +99,6 @@ list(APPEND SOURCES diag/diag_orbit.f90 ) -# Add GVEC-specific sources if enabled -if(ENABLE_GVEC) - list(APPEND SOURCES - field/field_gvec.f90 - field/vmec_field_adapter.f90 - ) -endif() - add_library (simple STATIC ${SOURCES}) if(SIMPLE_ENABLE_CGAL) @@ -142,17 +134,6 @@ if(CMAKE_Fortran_COMPILER_ID STREQUAL "NVHPC") ) endif() -# Conditionally link GVEC if enabled -if(ENABLE_GVEC) - target_include_directories(simple PRIVATE - ${gvec_BINARY_DIR}/include - ) - target_link_libraries(simple PUBLIC - gveclib - ) - target_compile_definitions(simple PUBLIC GVEC_AVAILABLE) -endif() - # Link fortplot (when available) if(NOT CMAKE_Fortran_COMPILER_ID STREQUAL "NVHPC") target_link_libraries(simple PUBLIC diff --git a/src/boozer_converter.F90 b/src/boozer_converter.F90 index 209e01f8..1e0b1090 100644 --- a/src/boozer_converter.F90 +++ b/src/boozer_converter.F90 @@ -482,11 +482,7 @@ subroutine compute_boozer_data use binsrc_sub, only: binsrc use plag_coeff_sub, only: plag_coeff use spline_vmec_sub -#ifdef GVEC_AVAILABLE - use vmec_field_adapter -#else use vmec_field_eval -#endif implicit none diff --git a/src/field.F90 b/src/field.F90 index 81b7cb4c..f1dabdbb 100644 --- a/src/field.F90 +++ b/src/field.F90 @@ -8,9 +8,6 @@ module field use field_vmec, only: vmec_field_t, create_vmec_field use field_coils, only: coils_field_t, create_coils_field use field_splined, only: splined_field_t, create_splined_field -#ifdef GVEC_AVAILABLE - use field_gvec, only: gvec_field_t, create_gvec_field -#endif implicit none @@ -45,16 +42,6 @@ subroutine field_clone(source, dest) class default error stop 'field_clone: Allocation failure (splined)' end select -#ifdef GVEC_AVAILABLE - type is (gvec_field_t) - allocate (gvec_field_t :: dest) - select type (dest) - type is (gvec_field_t) - dest = source - class default - error stop 'field_clone: Allocation failure (gvec)' - end select -#endif class default error stop 'field_clone: Unsupported field type' end select @@ -74,9 +61,6 @@ subroutine field_from_file(filename, field) type(vmec_field_t) :: vmec_field integer :: file_type, ierr character(len=2048) :: message -#ifdef GVEC_AVAILABLE - class(gvec_field_t), allocatable :: gvec_temp -#endif stripped_name = strip_directory(filename) @@ -116,14 +100,6 @@ subroutine field_from_file(filename, field) allocate (splined_coils) call create_splined_field(raw_coils, ref_coords, splined_coils) call move_alloc(splined_coils, field) - else if (endswith(filename, '.dat')) then -#ifdef GVEC_AVAILABLE - call create_gvec_field(filename, gvec_temp) - call move_alloc(gvec_temp, field) -#else - print *, 'ERROR: GVEC support not compiled. Rebuild with -DENABLE_GVEC=ON' - error stop -#endif else print *, 'field_from_file: Unknown file name format ', filename error stop diff --git a/src/field/field_gvec.f90 b/src/field/field_gvec.f90 deleted file mode 100644 index 2739055d..00000000 --- a/src/field/field_gvec.f90 +++ /dev/null @@ -1,213 +0,0 @@ -module field_gvec - !> GVEC magnetic field evaluation from .dat files. - - use, intrinsic :: iso_fortran_env, only: dp => real64 - use field_base, only: magnetic_field_t - use libneo_coordinates, only: coordinate_system_t, make_vmec_coordinate_system - use MODgvec_cubic_spline, only: t_cubspl - use MODgvec_rProfile_bspl, only: t_rProfile_bspl - use MODgvec_ReadState, only: ReadState, eval_phi_r, eval_phiPrime_r, & - eval_iota_r, eval_pres_r, Finalize_ReadState - use MODgvec_ReadState_Vars, only: profiles_1d, sbase_prof, X1_base_r, & - X2_base_r, LA_base_r, X1_r, X2_r, LA_r, hmap_r - - implicit none - - integer, parameter :: DERIV_R = 1 - integer, parameter :: DERIV_THET = 2 - integer, parameter :: DERIV_ZETA = 3 - - real(dp), parameter :: TESLA_IN_GAUSS = 10000.0_dp - real(dp), parameter :: METER_IN_CM = 100.0_dp - - private - public :: gvec_field_t, create_gvec_field - - type, extends(magnetic_field_t) :: gvec_field_t - character(len=256) :: filename = '' - integer :: nfp = 1 - logical :: data_loaded = .false. - contains - procedure :: evaluate => gvec_evaluate - procedure :: load_dat_file - final :: gvec_field_cleanup - end type gvec_field_t - -contains - - subroutine create_gvec_field(gvec_file, gvec_field) - character(*), intent(in) :: gvec_file - class(gvec_field_t), allocatable, intent(out) :: gvec_field - - allocate(gvec_field_t :: gvec_field) - gvec_field%filename = gvec_file - - call gvec_field%load_dat_file() - - if (.not. gvec_field%data_loaded) then - print *, 'ERROR: Failed to load GVEC field from: ', gvec_file - end if - - call make_vmec_coordinate_system(gvec_field%coords) - end subroutine create_gvec_field - - - subroutine load_dat_file(self) - class(gvec_field_t), intent(inout) :: self - logical :: file_exists - - inquire(file=trim(self%filename), exist=file_exists) - if (.not. file_exists) then - print *, 'Error: GVEC file does not exist: ', trim(self%filename) - return - end if - - call ReadState(trim(self%filename)) - - if (allocated(profiles_1d) .and. allocated(sbase_prof)) then - if (allocated(hmap_r)) then - self%nfp = hmap_r%nfp - else - self%nfp = 1 - print *, 'Warning: Could not read nfp from GVEC state, using default nfp=1' - end if - self%data_loaded = .true. - else - print *, 'Error: GVEC state initialization failed' - self%data_loaded = .false. - end if - end subroutine load_dat_file - - - subroutine gvec_evaluate(self, x, Acov, hcov, Bmod, sqgBctr) - class(gvec_field_t), intent(in) :: self - real(dp), intent(in) :: x(3) - real(dp), intent(out) :: Acov(3), hcov(3), Bmod - real(dp), intent(out), optional :: sqgBctr(3) - - real(dp) :: r, theta, varphi, zeta - real(dp) :: gvec_coords(3), RZ_coords(3) - integer :: deriv_flags(2) - real(dp) :: X1_val, X2_val, R_pos, Z_pos - real(dp) :: dX1_ds, dX1_dthet, dX1_dzeta - real(dp) :: dX2_ds, dX2_dthet, dX2_dzeta - real(dp) :: dLA_dr, dLA_dthet, dLA_dzeta - real(dp) :: iota_val, phi_val, phiPrime_val, chi_val - real(dp) :: e_thet(3), e_zeta(3), e_s(3) - real(dp) :: dx_dq1(3), dx_dq2(3), dx_dq3(3) - real(dp) :: g_tt, g_tz, g_zz, g_ss, g_st, g_sz - real(dp) :: Jac_h, Jac_l, Jac - real(dp) :: Bthctr, Bzetactr, Bthcov, Bzetacov, Bscov - - r = x(1) - theta = x(2) - varphi = x(3) - zeta = -varphi - - if (.not. allocated(profiles_1d) .or. .not. allocated(sbase_prof) .or. & - .not. allocated(X1_r) .or. .not. allocated(X2_r) .or. & - .not. allocated(LA_r)) then - call ReadState(trim(self%filename)) - if (.not. allocated(profiles_1d) .or. .not. allocated(sbase_prof) .or. & - .not. allocated(X1_r) .or. .not. allocated(X2_r) .or. & - .not. allocated(LA_r)) then - error stop 'Failed to initialize GVEC state for field evaluation' - end if - end if - - gvec_coords = [r, theta, zeta] - - deriv_flags = [0, 0] - X1_val = X1_base_r%evalDOF_x(gvec_coords, deriv_flags, X1_r) - - deriv_flags = [DERIV_R, 0] - dX1_ds = X1_base_r%evalDOF_x(gvec_coords, deriv_flags, X1_r) - - deriv_flags = [0, DERIV_THET] - dX1_dthet = X1_base_r%evalDOF_x(gvec_coords, deriv_flags, X1_r) - - deriv_flags = [0, DERIV_ZETA] - dX1_dzeta = X1_base_r%evalDOF_x(gvec_coords, deriv_flags, X1_r) - - deriv_flags = [0, 0] - X2_val = X2_base_r%evalDOF_x(gvec_coords, deriv_flags, X2_r) - - deriv_flags = [DERIV_R, 0] - dX2_ds = X2_base_r%evalDOF_x(gvec_coords, deriv_flags, X2_r) - - deriv_flags = [0, DERIV_THET] - dX2_dthet = X2_base_r%evalDOF_x(gvec_coords, deriv_flags, X2_r) - - deriv_flags = [0, DERIV_ZETA] - dX2_dzeta = X2_base_r%evalDOF_x(gvec_coords, deriv_flags, X2_r) - - deriv_flags = [DERIV_R, 0] - dLA_dr = LA_base_r%evalDOF_x(gvec_coords, deriv_flags, LA_r) - - deriv_flags = [0, DERIV_THET] - dLA_dthet = LA_base_r%evalDOF_x(gvec_coords, deriv_flags, LA_r) - - deriv_flags = [0, DERIV_ZETA] - dLA_dzeta = LA_base_r%evalDOF_x(gvec_coords, deriv_flags, LA_r) - - iota_val = eval_iota_r(r) - phi_val = eval_phi_r(r) * TESLA_IN_GAUSS * METER_IN_CM**2 - phiPrime_val = eval_phiPrime_r(r) * TESLA_IN_GAUSS * METER_IN_CM**2 - chi_val = sbase_prof%evalDOF_s(r, 0, profiles_1d(:, 2)) * & - TESLA_IN_GAUSS * METER_IN_CM**2 - - R_pos = X1_val - Z_pos = X2_val - RZ_coords = [R_pos, Z_pos, zeta] - call hmap_r%get_dx_dqi(RZ_coords, dx_dq1, dx_dq2, dx_dq3) - - e_s = dx_dq1 * dX1_ds + dx_dq2 * dX2_ds - e_thet = dx_dq1 * dX1_dthet + dx_dq2 * dX2_dthet - e_zeta = dx_dq1 * dX1_dzeta + dx_dq2 * dX2_dzeta + dx_dq3 - - g_ss = dot_product(e_s, e_s) * METER_IN_CM**2 - g_tt = dot_product(e_thet, e_thet) * METER_IN_CM**2 - g_zz = dot_product(e_zeta, e_zeta) * METER_IN_CM**2 - g_st = dot_product(e_s, e_thet) * METER_IN_CM**2 - g_sz = dot_product(e_s, e_zeta) * METER_IN_CM**2 - g_tz = dot_product(e_thet, e_zeta) * METER_IN_CM**2 - - Jac_h = hmap_r%eval_Jh(RZ_coords) - Jac_l = dX1_ds * dX2_dthet - dX1_dthet * dX2_ds - Jac = Jac_h * Jac_l * METER_IN_CM**3 - - Bthctr = (iota_val - dLA_dzeta) * phiPrime_val / Jac - Bzetactr = (1.0_dp + dLA_dthet) * phiPrime_val / Jac - - if (present(sqgBctr)) then - sqgBctr(1) = 0.0_dp - sqgBctr(2) = (-Jac) * Bthctr - sqgBctr(3) = (-Jac) * (-Bzetactr) - end if - - Bthcov = g_tt * Bthctr + g_tz * Bzetactr - Bzetacov = g_tz * Bthctr + g_zz * Bzetactr - Bscov = g_st * Bthctr + g_sz * Bzetactr - - Bmod = sqrt(Bthctr * Bthcov + Bzetactr * Bzetacov) - - hcov(1) = Bscov / Bmod - hcov(2) = Bthcov / Bmod - hcov(3) = -Bzetacov / Bmod - - Acov(1) = phi_val * dLA_dr - Acov(2) = phi_val * (1.0_dp + dLA_dthet) - Acov(3) = -(-chi_val + phi_val * dLA_dzeta) - end subroutine gvec_evaluate - - - subroutine gvec_field_cleanup(self) - type(gvec_field_t), intent(inout) :: self - - if (self%data_loaded) then - call Finalize_ReadState() - self%data_loaded = .false. - end if - end subroutine gvec_field_cleanup - -end module field_gvec diff --git a/src/field/field_newton.F90 b/src/field/field_newton.F90 index b956f45e..d1d45a2b 100644 --- a/src/field/field_newton.F90 +++ b/src/field/field_newton.F90 @@ -68,12 +68,7 @@ end subroutine newton_theta_from_canonical !> This abstracts the field-specific stream function evaluation subroutine get_stream_function(mag_field, s, theta, varphi, alam, dl_dt) use field, only: vmec_field_t -#ifdef GVEC_AVAILABLE - use field, only: gvec_field_t - use vmec_field_adapter, only: vmec_lambda_interpolate_with_field -#else use vmec_field_eval, only: vmec_lambda_interpolate_with_field -#endif class(magnetic_field_t), intent(in) :: mag_field real(dp), intent(in) :: s, theta, varphi @@ -85,13 +80,7 @@ subroutine get_stream_function(mag_field, s, theta, varphi, alam, dl_dt) type is (vmec_field_t) ! For VMEC, use the existing lambda interpolation call vmec_lambda_interpolate_with_field(mag_field, s, theta, varphi, alam, dl_dt) - -#ifdef GVEC_AVAILABLE - type is (gvec_field_t) - ! For GVEC, Lambda is available as LA - call vmec_lambda_interpolate_with_field(mag_field, s, theta, varphi, alam, dl_dt) -#endif - + class default ! For other fields, stream function may not be defined alam = 0.0_dp diff --git a/src/field/vmec_field_adapter.f90 b/src/field/vmec_field_adapter.f90 deleted file mode 100644 index 6989821b..00000000 --- a/src/field/vmec_field_adapter.f90 +++ /dev/null @@ -1,317 +0,0 @@ -module vmec_field_adapter - !> Module providing GVEC field evaluation functions with VMEC-compatible interface - !> This module extends vmec_field_eval to support GVEC fields - - use, intrinsic :: iso_fortran_env, only : dp => real64 - use field_base, only : magnetic_field_t - use field, only : vmec_field_t, gvec_field_t - use spline_vmec_sub - use MODgvec_ReadState, only: eval_iota_r - use MODgvec_ReadState_Vars, only: profiles_1d, sbase_prof, & - X1_base_r, X2_base_r, LA_base_r, & - X1_r, X2_r, LA_r - - implicit none - private - - ! GVEC derivative constants - integer, parameter :: DERIV_S = 1 - integer, parameter :: DERIV_THET = 2 - integer, parameter :: DERIV_ZETA = 3 - - ! Export functions - public :: vmec_field_evaluate, vmec_field_evaluate_with_field - public :: vmec_iota_interpolate, vmec_iota_interpolate_with_field - public :: vmec_lambda_interpolate, vmec_lambda_interpolate_with_field - public :: vmec_data_interpolate, vmec_data_interpolate_with_field - -contains - - ! Compute vector potential derivatives for GVEC field via numerical differentiation - subroutine compute_vector_potential_derivatives_gvec(field, s, theta, varphi, dA_theta_ds, dA_phi_ds) - class(magnetic_field_t), intent(in) :: field - real(dp), intent(in) :: s, theta, varphi - real(dp), intent(out) :: dA_theta_ds, dA_phi_ds - - real(dp) :: x(3), Acov(3), hcov(3), Bmod - real(dp) :: x_pert(3), Acov_pert(3), hcov_pert(3), Bmod_pert - real(dp), parameter :: ds = 1.0e-7_dp - - ! Check if GVEC state is properly initialized before evaluation - ! This follows the same pattern as field_gvec.f90 - if (.not. allocated(profiles_1d) .or. .not. allocated(sbase_prof) .or. & - .not. allocated(X1_r) .or. .not. allocated(X2_r) .or. .not. allocated(LA_r)) then - ! Need to reinitialize GVEC state - error stop 'GVEC state not initialized in vmec_field_adapter - call field%evaluate first' - end if - - ! Original position - x = [sqrt(s), theta, varphi] - call field%evaluate(x, Acov, hcov, Bmod) - - ! Perturb in s-direction to get dA_theta_ds - x_pert = [sqrt(s + ds), theta, varphi] - call field%evaluate(x_pert, Acov_pert, hcov_pert, Bmod_pert) - - dA_theta_ds = (Acov_pert(2) - Acov(2)) / ds - dA_phi_ds = (Acov_pert(3) - Acov(3)) / ds - end subroutine compute_vector_potential_derivatives_gvec - - !> VMEC field evaluation wrapper - subroutine vmec_field_evaluate(s, theta, varphi, & - A_theta, A_phi, dA_theta_ds, dA_phi_ds, aiota, & - sqg, alam, dl_ds, dl_dt, dl_dp, & - Bctrvr_vartheta, Bctrvr_varphi, & - Bcovar_r, Bcovar_vartheta, Bcovar_varphi) - real(dp), intent(in) :: s, theta, varphi - real(dp), intent(out) :: A_theta, A_phi, dA_theta_ds, dA_phi_ds - real(dp), intent(out) :: aiota, sqg, alam - real(dp), intent(out) :: dl_ds, dl_dt, dl_dp - real(dp), intent(out) :: Bctrvr_vartheta, Bctrvr_varphi - real(dp), intent(out) :: Bcovar_r, Bcovar_vartheta, Bcovar_varphi - - ! Call the existing VMEC routine - call vmec_field(s, theta, varphi, & - A_theta, A_phi, dA_theta_ds, dA_phi_ds, aiota, & - sqg, alam, dl_ds, dl_dt, dl_dp, & - Bctrvr_vartheta, Bctrvr_varphi, & - Bcovar_r, Bcovar_vartheta, Bcovar_varphi) - end subroutine vmec_field_evaluate - - !> Override evaluate function to handle GVEC fields (boozer_converter interface) - subroutine vmec_field_evaluate_with_field(field, s, theta, varphi, & - A_theta, A_phi, dA_theta_ds, dA_phi_ds, aiota, & - sqg, alam, dl_ds, dl_dt, dl_dp, & - Bctrvr_vartheta, Bctrvr_varphi, & - Bcovar_r, Bcovar_vartheta, Bcovar_varphi) - class(magnetic_field_t), intent(in) :: field - real(dp), intent(in) :: s, theta, varphi - real(dp), intent(out) :: A_theta, A_phi, dA_theta_ds, dA_phi_ds - real(dp), intent(out) :: aiota, sqg, alam - real(dp), intent(out) :: dl_ds, dl_dt, dl_dp - real(dp), intent(out) :: Bctrvr_vartheta, Bctrvr_varphi - real(dp), intent(out) :: Bcovar_r, Bcovar_vartheta, Bcovar_varphi - - real(dp) :: x(3), Acov(3), hcov(3), Bmod, sqgBctr(3) - real(dp) :: daiota_ds - integer :: deriv - - select type (field) - type is (vmec_field_t) - ! For VMEC fields, use the existing spline-based evaluation - call vmec_field_evaluate(s, theta, varphi, & - A_theta, A_phi, dA_theta_ds, dA_phi_ds, aiota, & - sqg, alam, dl_ds, dl_dt, dl_dp, & - Bctrvr_vartheta, Bctrvr_varphi, & - Bcovar_r, Bcovar_vartheta, Bcovar_varphi) - - type is (gvec_field_t) - ! For GVEC fields, we need to translate to GVEC coordinates and evaluate - x = [sqrt(s), theta, varphi] - - ! Get basic field quantities - call field%evaluate(x, Acov, hcov, Bmod) - - ! Extract vector potential components - A_theta = Acov(2) - A_phi = Acov(3) - - ! Get derivatives via numerical differentiation - call compute_vector_potential_derivatives_gvec(field, s, theta, varphi, dA_theta_ds, dA_phi_ds) - - ! Get iota - call vmec_iota_interpolate_with_field(field, s, aiota, daiota_ds) - - ! Get Lambda and its derivatives - call vmec_lambda_interpolate_with_field(field, s, theta, varphi, alam, dl_dt) - - ! For dl_ds and dl_dp, we need to compute numerically or from GVEC data - dl_ds = 0.0_dp ! Placeholder - needs proper derivative - dl_dp = 0.0_dp ! Placeholder - needs proper derivative - - ! Magnetic field components - ! These are covariant/contravariant components that need proper computation - sqg = sqrt(Bmod) ! Placeholder - needs proper Jacobian - Bctrvr_vartheta = hcov(2) * Bmod ! Placeholder - Bctrvr_varphi = hcov(3) * Bmod ! Placeholder - Bcovar_r = Acov(1) ! Placeholder - Bcovar_vartheta = Acov(2) ! Placeholder - Bcovar_varphi = Acov(3) ! Placeholder - - class default - error stop 'vmec_field_evaluate_with_field: Unsupported field type' - end select - end subroutine vmec_field_evaluate_with_field - - !> Override iota interpolation to handle GVEC fields - subroutine vmec_iota_interpolate_with_field(field, s, aiota, daiota_ds) - class(magnetic_field_t), intent(in) :: field - real(dp), intent(in) :: s - real(dp), intent(out) :: aiota, daiota_ds - - select type (field) - type is (vmec_field_t) - call vmec_iota_interpolate(s, aiota, daiota_ds) - - type is (gvec_field_t) - ! For GVEC, use eval_iota_r - ! GVEC eval_iota_r takes only radial coordinate - aiota = eval_iota_r(sqrt(s)) - ! For derivative, use finite difference - block - real(dp) :: rho, drho, iota_plus, iota_minus - rho = sqrt(s) - drho = 1.0e-7_dp - iota_plus = eval_iota_r(rho + drho) - iota_minus = eval_iota_r(rho - drho) - daiota_ds = (iota_plus - iota_minus) / (2.0_dp * drho) / (2.0_dp * sqrt(s)) ! Chain rule - end block - - class default - error stop 'vmec_iota_interpolate_with_field: Unsupported field type' - end select - end subroutine vmec_iota_interpolate_with_field - - !> VMEC iota interpolation wrapper - subroutine vmec_iota_interpolate(s, aiota, daiota_ds) - real(dp), intent(in) :: s - real(dp), intent(out) :: aiota, daiota_ds - - call splint_iota(s, aiota, daiota_ds) - end subroutine vmec_iota_interpolate - - !> Override lambda interpolation to handle GVEC fields - subroutine vmec_lambda_interpolate_with_field(field, s, theta, varphi, alam, dl_dt) - class(magnetic_field_t), intent(in) :: field - real(dp), intent(in) :: s, theta, varphi - real(dp), intent(out) :: alam, dl_dt - - real(dp) :: rho - - select type (field) - type is (vmec_field_t) - call vmec_lambda_interpolate(s, theta, varphi, alam, dl_dt) - - type is (gvec_field_t) - ! For GVEC, Lambda is available as LA - rho = sqrt(s) - - ! Use GVEC's Lambda evaluation following field_gvec.f90 pattern - ! LA_base_r is the base object, LA_r is the DOF array - block - real(dp) :: gvec_coords(3) - integer :: deriv_flags(2) - - gvec_coords = [rho, theta, varphi] - - ! No derivatives - deriv_flags = [0, 0] - alam = LA_base_r%evalDOF_x(gvec_coords, deriv_flags, LA_r) - - ! Derivative w.r.t. theta (DERIV_THET = 2) - deriv_flags = [0, 2] ! [radial_deriv, angular_deriv] - dl_dt = LA_base_r%evalDOF_x(gvec_coords, deriv_flags, LA_r) - end block - - class default - error stop 'vmec_lambda_interpolate_with_field: Unsupported field type' - end select - end subroutine vmec_lambda_interpolate_with_field - - !> VMEC lambda interpolation wrapper - subroutine vmec_lambda_interpolate(s, theta, varphi, alam, dl_dt) - real(dp), intent(in) :: s, theta, varphi - real(dp), intent(out) :: alam, dl_dt - - call splint_lambda(s, theta, varphi, alam, dl_dt) - end subroutine vmec_lambda_interpolate - - !> Override data interpolation to handle GVEC fields - subroutine vmec_data_interpolate_with_field(field, s, theta, varphi, & - A_phi, A_theta, dA_phi_ds, dA_theta_ds, aiota, & - R, Z, alam, & - dR_ds, dR_dt, dR_dp, & - dZ_ds, dZ_dt, dZ_dp, & - dl_ds, dl_dt, dl_dp) - class(magnetic_field_t), intent(in) :: field - real(dp), intent(in) :: s, theta, varphi - real(dp), intent(out) :: A_phi, A_theta, dA_phi_ds, dA_theta_ds - real(dp), intent(out) :: aiota, R, Z, alam - real(dp), intent(out) :: dR_ds, dR_dt, dR_dp - real(dp), intent(out) :: dZ_ds, dZ_dt, dZ_dp - real(dp), intent(out) :: dl_ds, dl_dt, dl_dp - - real(dp) :: x(3), Acov(3), hcov(3), Bmod - real(dp) :: daiota_ds - - select type (field) - type is (vmec_field_t) - call vmec_data_interpolate(s, theta, varphi, & - A_phi, A_theta, dA_phi_ds, dA_theta_ds, aiota, & - R, Z, alam, & - dR_ds, dR_dt, dR_dp, & - dZ_ds, dZ_dt, dZ_dp, & - dl_ds, dl_dt, dl_dp) - - type is (gvec_field_t) - ! For GVEC fields, we need to extract these from the field evaluation - x = [sqrt(s), theta, varphi] - call field%evaluate(x, Acov, hcov, Bmod) - - ! Vector potential components - A_theta = Acov(2) - A_phi = Acov(3) - - ! Get derivatives via numerical differentiation - call compute_vector_potential_derivatives_gvec(field, s, theta, varphi, dA_theta_ds, dA_phi_ds) - - ! Get iota - call vmec_iota_interpolate_with_field(field, s, aiota, daiota_ds) - - ! Get Lambda and its derivatives - call vmec_lambda_interpolate_with_field(field, s, theta, varphi, alam, dl_dt) - - ! For dl_ds and dl_dp, we need to compute numerically or from GVEC data - dl_ds = 0.0_dp ! Placeholder - needs proper derivative - dl_dp = 0.0_dp ! Placeholder - needs proper derivative - - ! Coordinate transformation values - these would need proper GVEC evaluation - ! For now, using placeholders - R = 1.0_dp - Z = 0.0_dp - dR_ds = 0.0_dp - dR_dt = 0.0_dp - dR_dp = 0.0_dp - dZ_ds = 0.0_dp - dZ_dt = 0.0_dp - dZ_dp = 0.0_dp - - class default - error stop 'vmec_data_interpolate_with_field: Unsupported field type' - end select - end subroutine vmec_data_interpolate_with_field - - !> VMEC data interpolation wrapper - subroutine vmec_data_interpolate(s, theta, varphi, & - A_phi, A_theta, dA_phi_ds, dA_theta_ds, aiota, & - R, Z, alam, & - dR_ds, dR_dt, dR_dp, & - dZ_ds, dZ_dt, dZ_dp, & - dl_ds, dl_dt, dl_dp) - real(dp), intent(in) :: s, theta, varphi - real(dp), intent(out) :: A_phi, A_theta, dA_phi_ds, dA_theta_ds - real(dp), intent(out) :: aiota, R, Z, alam - real(dp), intent(out) :: dR_ds, dR_dt, dR_dp - real(dp), intent(out) :: dZ_ds, dZ_dt, dZ_dp - real(dp), intent(out) :: dl_ds, dl_dt, dl_dp - - ! Call the existing VMEC routine - call splint_vmec_data(s, theta, varphi, & - A_phi, A_theta, dA_phi_ds, dA_theta_ds, aiota, & - R, Z, alam, & - dR_ds, dR_dt, dR_dp, & - dZ_ds, dZ_dt, dZ_dp, & - dl_ds, dl_dt, dl_dp) - end subroutine vmec_data_interpolate - -end module vmec_field_adapter diff --git a/src/field/vmec_field_eval.f90 b/src/field/vmec_field_eval.f90 index 7edea536..b3ec7b69 100644 --- a/src/field/vmec_field_eval.f90 +++ b/src/field/vmec_field_eval.f90 @@ -1,6 +1,5 @@ module vmec_field_eval !> Module providing VMEC field evaluation functions that work with magnetic_field_t classes - !> This module is always available regardless of GVEC support use, intrinsic :: iso_fortran_env, only : dp => real64 use field_base, only : magnetic_field_t diff --git a/src/get_canonical_coordinates.F90 b/src/get_canonical_coordinates.F90 index db3be568..f8d90ecc 100644 --- a/src/get_canonical_coordinates.F90 +++ b/src/get_canonical_coordinates.F90 @@ -412,11 +412,7 @@ subroutine rhs_cancoord(r, y, dy) Bcovar_vartheta, Bcovar_varphi, & theta, onlytheta use spline_vmec_sub -#ifdef GVEC_AVAILABLE - use vmec_field_adapter -#else use vmec_field_eval -#endif implicit none @@ -844,11 +840,7 @@ subroutine vmec_to_can(r, theta, varphi, vartheta_c, varphi_c) use new_vmec_stuff_mod, only: nper use vector_potentail_mod, only: torflux use chamb_mod, only: rnegflag -#ifdef GVEC_AVAILABLE - use vmec_field_adapter -#else use vmec_field_eval -#endif implicit none @@ -920,11 +912,7 @@ end subroutine vmec_to_can subroutine vmec_to_cyl(s, theta, varphi, Rcyl, Zcyl) use spline_vmec_sub -#ifdef GVEC_AVAILABLE - use vmec_field_adapter -#else use vmec_field_eval -#endif real(dp), intent(in) :: s, theta, varphi real(dp), intent(out) :: Rcyl, Zcyl diff --git a/test/golden_record/golden_record.sh b/test/golden_record/golden_record.sh index ae6d05da..a6464ee7 100755 --- a/test/golden_record/golden_record.sh +++ b/test/golden_record/golden_record.sh @@ -175,8 +175,6 @@ build() { if grep -q "fortplot" CMakeLists.txt 2>/dev/null; then echo "Detected fortplot dependency in CMakeLists.txt" - else - CMAKE_OPTS="$CMAKE_OPTS -DENABLE_GVEC=OFF" fi # Enable deterministic floating-point for reproducible golden record tests diff --git a/test/test_data/GVEC_elliptok_State_final.dat b/test/test_data/GVEC_elliptok_State_final.dat deleted file mode 100644 index 9ab634cb..00000000 --- a/test/test_data/GVEC_elliptok_State_final.dat +++ /dev/null @@ -1,52 +0,0 @@ -## MHD3D Solution... outputLevel and fileID: -0000,00000000 -## grid: nElems, gridType ################################################################################# - 11, 4 -## grid: sp(0:nElems) - 0.737806440415453E-18, 0.823498083671742E-01, 0.179250397055607E+00, 0.286422124794339E+00, 0.399585350312414E+00, 0.514460432338871E+00, 0.626767729602754E+00, 0.732227600833103E+00, 0.826560404758961E+00, 0.905486500109369E+00, 0.964726245613368E+00, 0.100000000000000E+01 -## global: nfp,degGP,mn_nyq(2),hmap ####################################################################### - 1, 8, 21, 1, 1 -## X1_base: s%nbase,s%deg,s%continuity,f%modes,f%sin_cos,f%excl_mn_zero ################################### - 16, 5, 4, 6, 2, 0 -## X2_base: s%nbase,s%deg,s%continuity,f%modes,f%sin_cos,f%excl_mn_zero ################################### - 16, 5, 4, 5, 1, 0 -## LA_base: s%nbase,s%deg,s%continuity,f%modes,f%sin_cos,f%excl_mn_zero ################################### - 15, 4, 3, 5, 1, 1 -## X1: m,n,X1(1:nbase,iMode) ############################################################################## - 0, 0, 0.200000000000000E+01, 0.200000000000000E+01, 0.200000000000000E+01, 0.200000000000000E+01, 0.200000000000000E+01, 0.200000000000000E+01, 0.200000000000000E+01, 0.200000000000000E+01, 0.200000000000000E+01, 0.200000000000000E+01, 0.200000000000000E+01, 0.200000000000000E+01, 0.200000000000000E+01, 0.200000000000000E+01, 0.200000000000000E+01, 0.200000000000000E+01 - 1, 0, 0.484617850462962E-18, 0.658798466937393E-02, 0.209280164338225E-01, 0.438417864173696E-01, 0.758086144423629E-01, 0.116965449029472E+00, 0.160518882728319E+00, 0.204757059030519E+00, 0.247968121427688E+00, 0.288440213411445E+00, 0.324461478473404E+00, 0.354320060105184E+00, 0.375741852038535E+00, 0.389617019657819E+00, 0.397178099649069E+00, 0.400000000000000E+00 - 2, 0, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00 - 3, 0, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00 - 4, 0, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, -0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00 - 5, 0, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, -0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00 -## X2: m,n,X2(1:nbase,iMode) ############################################################################## - 1, 0, 0.726926775694443E-18, 0.988197700406089E-02, 0.313920246507337E-01, 0.657626796260544E-01, 0.113712921663544E+00, 0.175448173544209E+00, 0.240778324092478E+00, 0.307135588545778E+00, 0.371952182141533E+00, 0.432660320117167E+00, 0.486692217710106E+00, 0.531480090157777E+00, 0.563612778057803E+00, 0.584425529486729E+00, 0.595767149473604E+00, 0.600000000000000E+00 - 2, 0, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00 - 3, 0, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00 - 4, 0, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, -0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00 - 5, 0, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, -0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00 -## LA: m,n,LA(1:nbase,iMode) ############################################################################## - 1, 0, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00 - 2, 0, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00 - 3, 0, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00 - 4, 0, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, -0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00 - 5, 0, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, -0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00, 0.000000000000000E+00 -## at X1_base IP point positions (size nBase): spos,phi,chi,iota,pressure ################################ - 0.737806440415453E-18, 0.669577468312232E-37, 0.557981223593527E-37, 0.833333333333333E+00, -0.316230250880711E-31 - 0.164699616734348E-01, 0.333657678827715E-04, 0.278038301168862E-04, 0.833274803613700E+00, -0.157575413756109E+02 - 0.523200410845562E-01, 0.336706964739609E-03, 0.280489713179826E-03, 0.832742811969689E+00, -0.158964733924160E+03 - 0.109604466043424E+00, 0.147765296118844E-02, 0.122946364186213E-02, 0.830743861874486E+00, -0.696786297523252E+03 - 0.189521536105907E+00, 0.441807498552379E-02, 0.366464345487230E-02, 0.825606767650835E+00, -0.207690071861911E+04 - 0.292413622573681E+00, 0.105174667880725E-01, 0.866800139925834E-02, 0.815016762806162E+00, -0.491250507463272E+04 - 0.401297206820797E+00, 0.198083463419793E-01, 0.161659095572980E-01, 0.799052610090220E+00, -0.916187116188146E+04 - 0.511892647576296E+00, 0.322309963200606E-01, 0.259612109997786E-01, 0.778010824508138E+00, -0.147132624611337E+05 - 0.619920303569221E+00, 0.472702248614846E-01, 0.374730745070992E-01, 0.752979109871759E+00, -0.212374985299391E+05 - 0.721100533528612E+00, 0.639598712552214E-01, 0.498119143794922E-01, 0.725733326248043E+00, -0.282304153657131E+05 - 0.811153696183511E+00, 0.809323684579932E-01, 0.618980992170514E-01, 0.698561565218547E+00, -0.350801424172721E+05 - 0.885800150262960E+00, 0.965133624467251E-01, 0.725909663270698E-01, 0.674058718670150E+00, -0.411402202900989E+05 - 0.939354630096339E+00, 0.108536323846437E+00, 0.805828750254613E-01, 0.655419822040941E+00, -0.456695563783200E+05 - 0.974042549144547E+00, 0.116700254809052E+00, 0.858824591485009E-01, 0.642888832294244E+00, -0.486730438618674E+05 - 0.992945249122674E+00, 0.121273678680780E+00, 0.888066953799562E-01, 0.635911328833906E+00, -0.503303259610340E+05 - 0.100000000000000E+01, 0.123003068894718E+00, 0.899041588010694E-01, 0.633280569441086E+00, -0.509523028029541E+05 -## a_minor,r_major,volume ################################################################################ - 0.489897948556636E+00, 0.199999999999999E+01, 0.947482022504577E+01 diff --git a/test/tests/CMakeLists.txt b/test/tests/CMakeLists.txt index 12744bb6..89d6a346 100644 --- a/test/tests/CMakeLists.txt +++ b/test/tests/CMakeLists.txt @@ -215,33 +215,6 @@ add_custom_target(ncsx_test_data DEPENDS ${CHARTMAP_NCSX_FILE} ${COILS_SIMPLE_FI # Symlinks for NCSX test data file(CREATE_LINK "${WOUT_NCSX_FILE}" "${CMAKE_CURRENT_BINARY_DIR}/wout_ncsx.nc" SYMBOLIC) -# Convert VMEC to GVEC format for testing field_vmec vs field_gvec -# This requires GVEC Python API or command-line tools -if(ENABLE_GVEC) - find_program(PYTHON_EXECUTABLE python3) - if(PYTHON_EXECUTABLE AND SIMPLE_ENABLE_PYTHON_TOOLS) - # Try to convert VMEC to GVEC using Python - execute_process( - COMMAND ${PYTHON_EXECUTABLE} -c " -try: - import gvec - print('GVEC Python module found, converting VMEC to GVEC format...') - gvec.read_vmec('${CMAKE_CURRENT_BINARY_DIR}/wout.nc') - gvec.save('${CMAKE_CURRENT_BINARY_DIR}/wout_gvec.dat') - print('Successfully converted VMEC to GVEC format: wout_gvec.dat') -except ImportError: - print('GVEC Python module not found, skipping VMEC to GVEC conversion') -except Exception as e: - print(f'Failed to convert VMEC to GVEC: {e}') -" - WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - RESULT_VARIABLE GVEC_CONVERT_RESULT - ) - else() - message(STATUS "Python not found, skipping VMEC to GVEC conversion") - endif() -endif() - ################################################## # INTERNAL UNIT TEST STUFF ################################################## @@ -324,53 +297,6 @@ add_executable (test_coordinates.x test_coordinates.f90) target_link_libraries(test_coordinates.x simple) add_test(NAME test_coordinates COMMAND test_coordinates.x vmec cyl) -# GVEC-related tests (only if GVEC support is enabled) -if(ENABLE_GVEC) - add_executable (test_gvec.x test_gvec.f90) - target_link_libraries(test_gvec.x simple) - # add_dependencies(test_gvec.x gvec_test_data) - add_test(NAME test_gvec COMMAND test_gvec.x) - - add_executable (test_vmec_gvec.x test_vmec_gvec.f90) - target_link_libraries(test_vmec_gvec.x simple) - add_test(NAME test_vmec_gvec COMMAND test_vmec_gvec.x) - - add_executable (test_canonical_gvec.x test_canonical_gvec.f90) - target_link_libraries(test_canonical_gvec.x simple) - add_test(NAME test_canonical_gvec COMMAND test_canonical_gvec.x) - - add_executable (test_simple_vmec_gvec.x test_simple_vmec_gvec.f90) - target_link_libraries(test_simple_vmec_gvec.x simple) - add_test(NAME test_simple_vmec_gvec COMMAND test_simple_vmec_gvec.x) - set_tests_properties(test_simple_vmec_gvec PROPERTIES - LABELS "slow;integration" - WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - WILL_FAIL TRUE # Expected to fail - ) - - add_executable (test_adapter_consistency.x test_adapter_consistency.f90) - target_link_libraries(test_adapter_consistency.x simple) - add_test(NAME test_adapter_consistency COMMAND test_adapter_consistency.x) - set_tests_properties(test_adapter_consistency PROPERTIES - WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - ) - - add_executable (test_vmec_gvec_adapter.x test_vmec_gvec_adapter.f90) - target_link_libraries(test_vmec_gvec_adapter.x simple) - add_test(NAME test_vmec_gvec_adapter COMMAND test_vmec_gvec_adapter.x) - set_tests_properties(test_vmec_gvec_adapter PROPERTIES - WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - LABELS "adapter" - ) -endif() - - -# Add 2D field export utility (requires GVEC support) -if(ENABLE_GVEC) - add_executable (export_field_2d.x export_field_2d.f90) - target_link_libraries(export_field_2d.x simple) -endif() - add_test(NAME test_magfie COMMAND test_magfie.x) # Unit tests for refactored modules diff --git a/test/tests/convert_vmec_to_gvec.py b/test/tests/convert_vmec_to_gvec.py deleted file mode 100644 index 3a4788d2..00000000 --- a/test/tests/convert_vmec_to_gvec.py +++ /dev/null @@ -1,99 +0,0 @@ -#!/usr/bin/env python3 -""" -Convert VMEC wout file to GVEC state file for testing. -This script uses GVEC's Python interface to read VMEC and save as GVEC format. -""" - -import sys -import os - -try: - # Add GVEC Python path if needed - gvec_path = os.path.join(os.path.dirname(__file__), - "../../../build/_deps/gvec-build/python") - if os.path.exists(gvec_path): - sys.path.insert(0, gvec_path) - - import gvec - - # First, let's see what's available in the GVEC module - print("GVEC module attributes:") - attrs = [attr for attr in dir(gvec) if not attr.startswith('_')] - for attr in attrs: - print(f" {attr}") - - # Check if there's a specific VMEC-related module or function - if hasattr(gvec, 'vmec'): - print("\nGVEC.vmec attributes:") - vmec_attrs = [attr for attr in dir(gvec.vmec) if not attr.startswith('_')] - for attr in vmec_attrs: - print(f" {attr}") - - # Check Run class - if hasattr(gvec, 'Run'): - print("\nGVEC.Run attributes:") - run_attrs = [attr for attr in dir(gvec.Run) if not attr.startswith('_')] - for attr in run_attrs: - print(f" {attr}") - - # Check State class - if hasattr(gvec, 'State'): - print("\nGVEC.State attributes:") - state_attrs = [attr for attr in dir(gvec.State) if not attr.startswith('_')] - for attr in state_attrs: - print(f" {attr}") - - # Check if there's a run function that might help - if hasattr(gvec, 'run'): - print("\nChecking gvec.run function...") - print(f" Type: {type(gvec.run)}") - if hasattr(gvec.run, '__doc__'): - print(f" Docstring: {gvec.run.__doc__}") - - # Try to create a GVEC run from VMEC - vmec_file = "wout.nc" - output_file = "wout_gvec.dat" - - print(f"\nAttempting to convert VMEC file: {vmec_file}") - - # Try using the Run class - if hasattr(gvec, 'Run'): - try: - print("Trying to create GVEC Run from VMEC...") - # Try different initialization methods - run = gvec.Run() # Create empty run - if hasattr(run, 'from_vmec'): - run.from_vmec(vmec_file) - elif hasattr(run, 'load_vmec'): - run.load_vmec(vmec_file) - elif hasattr(run, 'read_vmec'): - run.read_vmec(vmec_file) - else: - print(" Run object methods:") - for attr in dir(run): - if not attr.startswith('_'): - print(f" {attr}") - except Exception as e: - print(f" Failed with Run class: {e}") - - # Try using load_state if available - if hasattr(gvec, 'load_state'): - print("\nChecking load_state function...") - print(f" Type: {type(gvec.load_state)}") - if hasattr(gvec.load_state, '__doc__'): - print(f" Docstring: {gvec.load_state.__doc__}") - - print("\nERROR: Could not determine correct GVEC API for VMEC conversion") - print("Please check GVEC documentation for the correct Python API") - sys.exit(1) - - print("Conversion successful!") - sys.exit(0) - -except ImportError as e: - print(f"ERROR: Could not import gvec module: {e}") - print("GVEC Python module not available") - sys.exit(1) -except Exception as e: - print(f"ERROR: {e}") - sys.exit(1) \ No newline at end of file diff --git a/test/tests/export_field_2d.f90 b/test/tests/export_field_2d.f90 index f236614f..3c7bb6be 100644 --- a/test/tests/export_field_2d.f90 +++ b/test/tests/export_field_2d.f90 @@ -1,24 +1,15 @@ program export_field_2d use, intrinsic :: iso_fortran_env, only: dp => real64 use field_vmec, only: vmec_field_t -#ifdef GVEC_AVAILABLE - use field_gvec, only: gvec_field_t, create_gvec_field -#endif use new_vmec_stuff_mod, only: netcdffile, multharm use spline_vmec_sub, only: spline_vmec_data use params, only: pi - + implicit none - + class(vmec_field_t), allocatable :: vmec_field -#ifdef GVEC_AVAILABLE - class(gvec_field_t), allocatable :: gvec_field -#endif character(len=256) :: vmec_file -#ifdef GVEC_AVAILABLE - character(len=256) :: gvec_file -#endif - + ! Grid parameters integer, parameter :: ns = 50 ! Number of s points integer, parameter :: nt = 100 ! Number of theta points @@ -27,56 +18,40 @@ program export_field_2d real(dp) :: theta_min = 0.0_dp real(dp) :: theta_max = 2.0_dp * pi real(dp) :: phi_fixed = 0.0_dp ! Fixed toroidal angle for 2D slice - + ! Field evaluation variables real(dp) :: x(3) real(dp) :: Acov_vmec(3), hcov_vmec(3), Bmod_vmec -#ifdef GVEC_AVAILABLE - real(dp) :: Acov_gvec(3), hcov_gvec(3), Bmod_gvec -#endif - + ! Grid variables real(dp), allocatable :: s_grid(:), theta_grid(:) real(dp), allocatable :: Bmod_vmec_2d(:,:) -#ifdef GVEC_AVAILABLE - real(dp), allocatable :: Bmod_gvec_2d(:,:) -#endif real(dp), allocatable :: R_grid(:,:), Z_grid(:,:) - + integer :: i, j, unit_out real(dp) :: s, theta, r - + ! File names vmec_file = 'wout.nc' -#ifdef GVEC_AVAILABLE - gvec_file = 'test_vmec_gvec_State_0000_00000000.dat' -#endif - + print *, '================================================================' -#ifdef GVEC_AVAILABLE - print *, 'Exporting 2D magnetic field data for VMEC and GVEC' -#else print *, 'Exporting 2D magnetic field data for VMEC' -#endif print *, '================================================================' - + ! Allocate grids allocate(s_grid(ns), theta_grid(nt)) allocate(Bmod_vmec_2d(ns, nt)) -#ifdef GVEC_AVAILABLE - allocate(Bmod_gvec_2d(ns, nt)) -#endif allocate(R_grid(ns, nt), Z_grid(ns, nt)) - + ! Create grids do i = 1, ns s_grid(i) = s_min + (s_max - s_min) * real(i-1, dp) / real(ns-1, dp) end do - + do j = 1, nt theta_grid(j) = theta_min + (theta_max - theta_min) * real(j-1, dp) / real(nt-1, dp) end do - + ! Initialize VMEC field print *, '' print *, 'Loading VMEC field from ', trim(vmec_file) @@ -84,52 +59,40 @@ program export_field_2d multharm = 7 call spline_vmec_data allocate(vmec_field_t :: vmec_field) - -#ifdef GVEC_AVAILABLE - ! Initialize GVEC field - print *, 'Loading GVEC field from ', trim(gvec_file) - call create_gvec_field(gvec_file, gvec_field) -#endif - + ! Evaluate fields on 2D grid print *, '' - print *, 'Evaluating fields on 2D grid (s, theta) at phi = ', phi_fixed/pi, 'π' + print *, 'Evaluating fields on 2D grid (s, theta) at phi = ', phi_fixed/pi, 'pi' print *, 'Grid size: ', ns, ' x ', nt - + do i = 1, ns s = s_grid(i) r = sqrt(s) ! r = sqrt(s) - + do j = 1, nt theta = theta_grid(j) - + ! Set coordinates (r = sqrt(s), theta, phi) x(1) = r x(2) = theta x(3) = phi_fixed - + ! Evaluate VMEC field call vmec_field%evaluate(x, Acov_vmec, hcov_vmec, Bmod_vmec) Bmod_vmec_2d(i, j) = Bmod_vmec - -#ifdef GVEC_AVAILABLE - ! Evaluate GVEC field - call gvec_field%evaluate(x, Acov_gvec, hcov_gvec, Bmod_gvec) - Bmod_gvec_2d(i, j) = Bmod_gvec -#endif - + ! For visualization, we'll also need R,Z coordinates ! Approximate R,Z for plotting (assuming axisymmetric at phi=0) ! This is a rough approximation - in reality we'd compute from VMEC harmonics R_grid(i, j) = 5.0_dp + r * cos(theta) ! Major radius ~ 5m Z_grid(i, j) = r * sin(theta) end do - + if (mod(i, 10) == 0) then print '(A,I3,A)', ' Progress: ', nint(100.0_dp * real(i, dp) / real(ns, dp)), '%' end if end do - + ! Write VMEC data print *, '' print *, 'Writing VMEC field data to Bmod_vmec_2d.dat' @@ -138,7 +101,7 @@ program export_field_2d write(unit_out, '(A,F8.4,A)') '# Fixed phi = ', phi_fixed/pi, ' pi' write(unit_out, '(A,I5,A,I5)') '# Grid size: ns = ', ns, ', nt = ', nt write(unit_out, '(A)') '# Columns: s, theta, R, Z, |B|' - + do i = 1, ns do j = 1, nt write(unit_out, '(5ES16.8)') s_grid(i), theta_grid(j), & @@ -147,26 +110,7 @@ program export_field_2d write(unit_out, *) ! Empty line for gnuplot end do close(unit_out) - -#ifdef GVEC_AVAILABLE - ! Write GVEC data - print *, 'Writing GVEC field data to Bmod_gvec_2d.dat' - open(newunit=unit_out, file='Bmod_gvec_2d.dat', status='replace') - write(unit_out, '(A)') '# 2D magnetic field magnitude from GVEC' - write(unit_out, '(A,F8.4,A)') '# Fixed phi = ', phi_fixed/pi, ' pi' - write(unit_out, '(A,I5,A,I5)') '# Grid size: ns = ', ns, ', nt = ', nt - write(unit_out, '(A)') '# Columns: s, theta, R, Z, |B|' - - do i = 1, ns - do j = 1, nt - write(unit_out, '(5ES16.8)') s_grid(i), theta_grid(j), & - R_grid(i, j), Z_grid(i, j), Bmod_gvec_2d(i, j) - end do - write(unit_out, *) ! Empty line for gnuplot - end do - close(unit_out) -#endif - + ! Write grid info for Python script print *, 'Writing grid information to grid_info.dat' open(newunit=unit_out, file='grid_info.dat', status='replace') @@ -174,28 +118,21 @@ program export_field_2d write(unit_out, '(4ES16.8)') s_min, s_max, theta_min, theta_max write(unit_out, '(ES16.8)') phi_fixed close(unit_out) - + ! Print summary statistics print *, '' print *, 'Field statistics:' print '(A,ES12.5,A,ES12.5)', ' VMEC |B| range: ', & minval(Bmod_vmec_2d), ' to ', maxval(Bmod_vmec_2d) -#ifdef GVEC_AVAILABLE - print '(A,ES12.5,A,ES12.5)', ' GVEC |B| range: ', & - minval(Bmod_gvec_2d), ' to ', maxval(Bmod_gvec_2d) -#endif - + ! Deallocate deallocate(s_grid, theta_grid) deallocate(Bmod_vmec_2d) -#ifdef GVEC_AVAILABLE - deallocate(Bmod_gvec_2d) -#endif deallocate(R_grid, Z_grid) - + print *, '' print *, 'Export completed successfully!' print *, 'Use plot_fields_2d.py to visualize the results' print *, '================================================================' - -end program export_field_2d \ No newline at end of file + +end program export_field_2d diff --git a/test/tests/test_adapter_consistency.f90 b/test/tests/test_adapter_consistency.f90 deleted file mode 100644 index 0646dc15..00000000 --- a/test/tests/test_adapter_consistency.f90 +++ /dev/null @@ -1,181 +0,0 @@ -program test_adapter_consistency - !> Test that vmec_field_adapter gives consistent results - !> between direct VMEC calls and field object calls - - use, intrinsic :: iso_fortran_env, only : dp => real64 - use vmec_field_adapter, only : vmec_field_evaluate, vmec_field_evaluate_with_field - use field_gvec, only : gvec_field_t, create_gvec_field - use field_base, only : magnetic_field_t - use params, only : pi - - implicit none - - class(gvec_field_t), allocatable :: gvec_field - character(len=256) :: gvec_file - logical :: file_exists, test_passed - - ! Test point - real(dp) :: s_test = 0.5_dp - real(dp) :: theta_test = pi / 3.0_dp - real(dp) :: phi_test = pi / 6.0_dp - - ! Field outputs from direct adapter call - real(dp) :: A_theta_1, A_phi_1, dA_theta_ds_1, dA_phi_ds_1 - real(dp) :: aiota_1, sqg_1, alam_1 - real(dp) :: dl_ds_1, dl_dt_1, dl_dp_1 - real(dp) :: Bctrvr_vartheta_1, Bctrvr_varphi_1 - real(dp) :: Bcovar_r_1, Bcovar_vartheta_1, Bcovar_varphi_1 - - ! Field outputs from field object call - real(dp) :: A_theta_2, A_phi_2, dA_theta_ds_2, dA_phi_ds_2 - real(dp) :: aiota_2, sqg_2, alam_2 - real(dp) :: dl_ds_2, dl_dt_2, dl_dp_2 - real(dp) :: Bctrvr_vartheta_2, Bctrvr_varphi_2 - real(dp) :: Bcovar_r_2, Bcovar_vartheta_2, Bcovar_varphi_2 - - real(dp) :: max_diff - - ! Dummy variables for field initialization - real(dp) :: x_test(3), Acov_dummy(3), hcov_dummy(3), Bmod_dummy - - print *, '=======================================================' - print *, 'Testing vmec_field_adapter consistency' - print *, '=======================================================' - print *, '' - - ! Look for GVEC test file - gvec_file = 'test_vmec_gvec_State_0000_00000000.dat' - inquire(file=gvec_file, exist=file_exists) - if (.not. file_exists) then - print *, 'FAIL: No GVEC file found. Run test_vmec_gvec first.' - error stop 1 - end if - - ! Load GVEC field - call create_gvec_field(gvec_file, gvec_field) - - select type(field => gvec_field) - type is (gvec_field_t) - if (.not. field%data_loaded) then - print *, 'ERROR: Failed to load GVEC field data' - error stop 1 - end if - end select - - print '(A,3(F6.3,A))', 'Test point: (s,θ,φ) = (', & - s_test, ', ', theta_test, ', ', phi_test, ')' - print *, '' - - ! Test 1: Direct adapter call would only work with VMEC loaded - print *, '1. Skipping direct adapter call (requires VMEC data loaded)...' - print *, ' The legacy interface without field object only works with VMEC' - print *, '' - - ! Set values to match field call for comparison - A_theta_1 = 0.0_dp - A_phi_1 = 0.0_dp - dA_theta_ds_1 = 0.0_dp - dA_phi_ds_1 = 0.0_dp - aiota_1 = 0.0_dp - sqg_1 = 0.0_dp - alam_1 = 0.0_dp - dl_ds_1 = 0.0_dp - dl_dt_1 = 0.0_dp - dl_dp_1 = 0.0_dp - Bctrvr_vartheta_1 = 0.0_dp - Bctrvr_varphi_1 = 0.0_dp - Bcovar_r_1 = 0.0_dp - Bcovar_vartheta_1 = 0.0_dp - Bcovar_varphi_1 = 0.0_dp - - ! Test 2: Field object call (uses generic field interface) - print *, '2. Testing field object call (new interface)...' - - ! First initialize GVEC state by calling field%evaluate - x_test = [sqrt(s_test), theta_test, phi_test] - call gvec_field%evaluate(x_test, Acov_dummy, hcov_dummy, Bmod_dummy) - - call vmec_field_evaluate_with_field(gvec_field, s_test, theta_test, phi_test, & - A_theta_2, A_phi_2, dA_theta_ds_2, dA_phi_ds_2, aiota_2, & - sqg_2, alam_2, dl_ds_2, dl_dt_2, dl_dp_2, & - Bctrvr_vartheta_2, Bctrvr_varphi_2, & - Bcovar_r_2, Bcovar_vartheta_2, Bcovar_varphi_2) - - print '(A,ES12.5)', ' sqrt(g) = ', sqg_2 - print '(A,ES12.5)', ' iota = ', aiota_2 - print '(A,ES12.5)', ' Lambda = ', alam_2 - print *, '' - - ! Validate GVEC results - test_passed = .true. - - print *, '3. Validating GVEC field results...' - - ! Check that essential quantities are non-zero and reasonable - if (abs(sqg_2) < 1.0e-6_dp) then - print *, ' FAIL: sqrt(g) too small: ', sqg_2 - test_passed = .false. - else - print '(A,ES12.5)', ' OK: sqrt(g) = ', sqg_2 - end if - - if (abs(aiota_2) < 1.0e-3_dp) then - print *, ' FAIL: iota too small: ', aiota_2 - test_passed = .false. - else - print '(A,ES12.5)', ' OK: iota = ', aiota_2 - end if - - if (abs(Bctrvr_vartheta_2) < 1.0e-10_dp .and. abs(Bctrvr_varphi_2) < 1.0e-10_dp) then - print *, ' FAIL: Both B contravariant components are zero' - test_passed = .false. - else - print '(A,ES12.5)', ' OK: B^theta = ', Bctrvr_vartheta_2 - print '(A,ES12.5)', ' OK: B^phi = ', Bctrvr_varphi_2 - end if - - ! Check vector potential derivatives are computed - if (abs(dA_theta_ds_2) < 1.0e-15_dp .and. abs(dA_phi_ds_2) < 1.0e-15_dp) then - print *, ' WARNING: Both vector potential derivatives are zero' - print *, ' This may indicate numerical differentiation issues' - else - print '(A,ES12.5)', ' OK: dA_theta/ds = ', dA_theta_ds_2 - print '(A,ES12.5)', ' OK: dA_phi/ds = ', dA_phi_ds_2 - end if - - print *, '' - print *, '=======================================================' - - if (test_passed) then - print *, 'TEST PASSED: Adapter correctly handles GVEC fields' - print *, '- All essential field quantities computed correctly' - print *, '- GVEC field works with canonical coordinates via adapter' - else - print *, 'TEST FAILED: Issues with GVEC field adapter' - error stop 1 - end if - -contains - - subroutine check_consistency(name, val1, val2, passed, max_diff) - character(*), intent(in) :: name - real(dp), intent(in) :: val1, val2 - logical, intent(inout) :: passed - real(dp), intent(inout) :: max_diff - - real(dp) :: abs_diff - - abs_diff = abs(val1 - val2) - if (abs_diff > max_diff) max_diff = abs_diff - - if (abs_diff > 1.0e-10_dp) then - print '(A,A,A,2(ES12.5,A),ES12.5)', ' FAIL: ', name, & - ' - Direct: ', val1, ', Field: ', val2, & - ', diff: ', abs_diff - passed = .false. - else - print '(A,A,A,ES12.5)', ' OK: ', name, ' - diff: ', abs_diff - end if - end subroutine check_consistency - -end program test_adapter_consistency diff --git a/test/tests/test_canonical_gvec.f90 b/test/tests/test_canonical_gvec.f90 deleted file mode 100644 index 2f55c0ca..00000000 --- a/test/tests/test_canonical_gvec.f90 +++ /dev/null @@ -1,176 +0,0 @@ -program test_canonical_gvec - !> Integration test for canonical coordinates with GVEC fields - !> Tests that the vmec_field_adapter properly handles GVEC fields - !> and that vector potential derivatives are computed correctly - - use, intrinsic :: iso_fortran_env, only : dp => real64 - use field_gvec, only : gvec_field_t, create_gvec_field - use field_base, only : magnetic_field_t - use vmec_field_adapter, only : vmec_field_evaluate_with_field, & - vmec_lambda_interpolate_with_field, & - vmec_iota_interpolate_with_field - use params, only : pi - - implicit none - - class(gvec_field_t), allocatable :: gvec_field - character(len=256) :: gvec_file - logical :: file_exists - - ! Test coordinates - real(dp) :: s_test, theta_test, varphi_test - - ! Field evaluation outputs - real(dp) :: A_theta, A_phi, dA_theta_ds, dA_phi_ds - real(dp) :: aiota, sqg, alam - real(dp) :: dl_ds, dl_dt, dl_dp - real(dp) :: Bctrvr_vartheta, Bctrvr_varphi - real(dp) :: Bcovar_r, Bcovar_vartheta, Bcovar_varphi - - ! Lambda interpolation outputs - real(dp) :: alam_interp, dl_dt_interp - - ! Iota interpolation outputs - real(dp) :: aiota_interp, daiota_ds - - ! Test flags - logical :: test_passed - - ! Dummy variables for field initialization - real(dp) :: x_test(3), Acov_dummy(3), hcov_dummy(3), Bmod_dummy - - print *, 'Testing canonical coordinates support with GVEC fields...' - print *, '' - - ! Look for GVEC test file created by test_vmec_gvec - gvec_file = 'test_vmec_gvec_State_0000_00000000.dat' - inquire(file=gvec_file, exist=file_exists) - - if (.not. file_exists) then - ! Try alternative name - gvec_file = 'gvec_from_vmec_wout.dat' - inquire(file=gvec_file, exist=file_exists) - end if - - if (.not. file_exists) then - print *, 'FAIL: No GVEC test file found. Run test_vmec_gvec first.' - error stop 1 - end if - - ! Load GVEC field - call create_gvec_field(gvec_file, gvec_field) - - if (.not. gvec_field%data_loaded) then - print *, 'ERROR: Failed to load GVEC field data' - error stop 1 - end if - - print *, 'Successfully loaded GVEC field from: ', trim(gvec_file) - print *, '' - - ! Test coordinates - representative interior point - s_test = 0.25_dp ! Mid-radius flux surface - theta_test = pi / 3.0_dp ! 60 degrees poloidal - varphi_test = pi / 6.0_dp ! 30 degrees toroidal - - print *, 'Testing adapter functions at (s,θ,φ) = (0.25, π/3, π/6):' - print *, '================================================================' - - ! Test 1: Full field evaluation with vector potential derivatives - print *, '1. Testing vmec_field_evaluate_with_field...' - - ! First initialize GVEC state by calling field%evaluate - x_test = [sqrt(s_test), theta_test, varphi_test] - call gvec_field%evaluate(x_test, Acov_dummy, hcov_dummy, Bmod_dummy) - - call vmec_field_evaluate_with_field(gvec_field, s_test, theta_test, varphi_test, & - A_theta, A_phi, dA_theta_ds, dA_phi_ds, aiota, & - sqg, alam, dl_ds, dl_dt, dl_dp, & - Bctrvr_vartheta, Bctrvr_varphi, & - Bcovar_r, Bcovar_vartheta, Bcovar_varphi) - - print '(A,ES12.5)', ' A_theta = ', A_theta - print '(A,ES12.5)', ' A_phi = ', A_phi - print '(A,ES12.5)', ' dA_theta/ds = ', dA_theta_ds - print '(A,ES12.5)', ' dA_phi/ds = ', dA_phi_ds - print '(A,ES12.5)', ' aiota = ', aiota - print '(A,ES12.5)', ' sqrt(g) = ', sqg - print '(A,ES12.5)', ' Lambda = ', alam - print '(A,ES12.5)', ' dLambda/ds = ', dl_ds - print '(A,ES12.5)', ' dLambda/dθ = ', dl_dt - print '(A,ES12.5)', ' dLambda/dφ = ', dl_dp - print *, '' - - ! Test 2: Lambda interpolation - print *, '2. Testing vmec_lambda_interpolate_with_field...' - - call vmec_lambda_interpolate_with_field(gvec_field, s_test, theta_test, varphi_test, & - alam_interp, dl_dt_interp) - - print '(A,ES12.5)', ' Lambda = ', alam_interp - print '(A,ES12.5)', ' dLambda/dθ = ', dl_dt_interp - print '(A,ES12.5)', ' Consistency: Lambda error = ', abs(alam_interp - alam) - print '(A,ES12.5)', ' dΛ/dθ error = ', abs(dl_dt_interp - dl_dt) - print *, '' - - ! Test 3: Iota interpolation - print *, '3. Testing vmec_iota_interpolate_with_field...' - - call vmec_iota_interpolate_with_field(gvec_field, s_test, aiota_interp, daiota_ds) - - print '(A,ES12.5)', ' iota = ', aiota_interp - print '(A,ES12.5)', ' diota/ds = ', daiota_ds - print '(A,ES12.5)', ' Consistency: iota error = ', abs(aiota_interp - aiota) - print *, '' - - ! Test validation criteria - test_passed = .true. - - ! Check that essential quantities are non-zero and reasonable - if (abs(aiota) < 1.0e-6_dp) then - print *, 'ERROR: Rotational transform too small: ', aiota - test_passed = .false. - end if - - if (abs(sqg) < 1.0e-6_dp) then - print *, 'ERROR: Jacobian too small: ', sqg - test_passed = .false. - end if - - ! Check that vector potential derivatives are being computed (non-zero) - if (abs(dA_theta_ds) < 1.0e-15_dp .and. abs(dA_phi_ds) < 1.0e-15_dp) then - print *, 'WARNING: Both vector potential derivatives are zero' - print *, ' This suggests numerical differentiation may not be working' - end if - - ! Check consistency between different interfaces - if (abs(alam_interp - alam) > 1.0e-10_dp) then - print *, 'ERROR: Lambda values inconsistent between interfaces' - test_passed = .false. - end if - - if (abs(dl_dt_interp - dl_dt) > 1.0e-10_dp) then - print *, 'ERROR: dLambda/dtheta values inconsistent between interfaces' - test_passed = .false. - end if - - if (abs(aiota_interp - aiota) > 1.0e-10_dp) then - print *, 'ERROR: iota values inconsistent between interfaces' - test_passed = .false. - end if - - print *, '================================================================' - print *, '' - - if (test_passed) then - print *, 'TEST PASSED: GVEC field adapter functions working correctly' - print *, '- Vector potential derivatives computed via numerical differentiation' - print *, '- Stream function Lambda and derivatives available' - print *, '- Rotational transform iota and derivatives available' - print *, '- All adapter interfaces consistent' - else - print *, 'TEST FAILED: Issues found with GVEC field adapter' - error stop 1 - end if - -end program test_canonical_gvec diff --git a/test/tests/test_gvec.f90 b/test/tests/test_gvec.f90 deleted file mode 100644 index 6ab89549..00000000 --- a/test/tests/test_gvec.f90 +++ /dev/null @@ -1,64 +0,0 @@ -program test_gvec - -use, intrinsic :: iso_fortran_env, only: dp => real64 -use field_gvec, only: gvec_field_t, create_gvec_field - -implicit none - -class(gvec_field_t), allocatable :: gvec_field -character(len=256) :: test_file - -! Test data file path - using relative path from build/test/tests/ directory -test_file = '../../../test/test_data/GVEC_elliptok_State_final.dat' - -print *, 'Testing gvec_field_t creation...' - -! Test: Create gvec_field_t using constructor -print *, 'Creating gvec_field_t with create_gvec_field...' -call create_gvec_field(test_file, gvec_field) - -if (allocated(gvec_field)) then - print *, 'SUCCESS: gvec_field_t created successfully' - print *, 'Filename stored: ', trim(gvec_field%filename) -else - print *, 'FAILED: gvec_field_t creation failed' - error stop 1 -end if - -! Test the evaluate method -print *, '' -print *, 'Testing gvec_field_t evaluate method...' - -block - real(dp) :: x(3), Acov(3), hcov(3), Bmod - real(dp) :: s_test, theta_test, phi_test - - ! Test at a few different flux surfaces - s_test = 0.5_dp ! Half-radius - theta_test = 0.0_dp ! Poloidal angle - phi_test = 0.0_dp ! Toroidal angle - - x(1) = sqrt(s_test) ! r = sqrt(s) - x(2) = theta_test ! theta - x(3) = phi_test ! phi - - call gvec_field%evaluate(x, Acov, hcov, Bmod) - - print *, 'Field evaluation at s=0.5:' - print *, ' Acov = ', Acov - print *, ' hcov = ', hcov - print *, ' Bmod = ', Bmod - - if (Bmod > 0.0_dp) then - print *, 'SUCCESS: evaluate() method working' - else - print *, 'ERROR: evaluate() returned invalid Bmod' - error stop 1 - end if -end block - -print *, '' -print *, 'All tests passed!' -print *, 'gvec_field_t constructor and evaluate() working correctly.' - -end program test_gvec diff --git a/test/tests/test_simple_vmec_gvec.f90 b/test/tests/test_simple_vmec_gvec.f90 deleted file mode 100644 index 42228f40..00000000 --- a/test/tests/test_simple_vmec_gvec.f90 +++ /dev/null @@ -1,144 +0,0 @@ -program test_simple_vmec_gvec - !> Comprehensive test running SIMPLE with both VMEC and GVEC fields - !> Verifies particle confinement statistics are consistent - - use, intrinsic :: iso_fortran_env, only: dp => real64, int64 - use params, only: pi - - implicit none - - character(len=256) :: vmec_file, gvec_file - character(len=256) :: namelist_file - logical :: file_exists - real(dp) :: confined_fraction_vmec, confined_fraction_gvec - real(dp) :: relative_difference - integer :: iostat, unit - - print *, '=======================================================' - print *, 'Testing SIMPLE with both VMEC and GVEC fields' - print *, '=======================================================' - print *, '' - - ! Check for required files - vmec_file = 'wout.nc' - inquire(file=vmec_file, exist=file_exists) - if (.not. file_exists) then - print *, 'FAIL: No VMEC file (wout.nc) found in working directory' - error stop 1 - end if - - ! Look for existing GVEC file created by test_vmec_gvec - gvec_file = 'test_vmec_gvec_State_0000_00000000.dat' - inquire(file=gvec_file, exist=file_exists) - if (.not. file_exists) then - print *, 'FAIL: No GVEC test file found. Run test_vmec_gvec first.' - error stop 1 - end if - - ! Test 1: Run with VMEC field (canonical coordinates) - print *, '1. Running SIMPLE with VMEC field (canonical coordinates)...' - - ! Create namelist for VMEC test - namelist_file = 'test_vmec.in' - open(newunit=unit, file=namelist_file, status='replace') - write(unit, '(A)') '&config' - write(unit, '(A)') 'multharm = 3 ! low order splines' - write(unit, '(A)') 'trace_time = 1d-4 ! very short test run' - write(unit, '(A)') 'sbeg = 0.5d0 ! mid-radius' - write(unit, '(A)') 'ntestpart = 16 ! minimal for quick test' - write(unit, '(A)') 'contr_pp = -1.0d10 ! trace also passing' - write(unit, '(A)') '/' - close(unit) - - ! Run SIMPLE with VMEC - call system('../../simple.x test_vmec.in > test_vmec.log') - - ! Read confined fraction - open(newunit=unit, file='confined_fraction.dat', status='old', iostat=iostat) - if (iostat /= 0) then - print *, 'ERROR: Could not read confined_fraction.dat for VMEC test' - error stop 1 - end if - - ! Skip header and read last line - read(unit, *) ! skip header - do - read(unit, *, iostat=iostat) confined_fraction_vmec - if (iostat /= 0) exit - end do - close(unit) - - print '(A,F8.4)', ' Final confined fraction: ', confined_fraction_vmec - call system('mv confined_fraction.dat confined_fraction_vmec.dat') - call system('mv times_lost.dat times_lost_vmec.dat') - - ! Test 2: Run with GVEC field (canonical coordinates) - print *, '' - print *, '2. Running SIMPLE with GVEC field (canonical coordinates)...' - - ! Create namelist for GVEC test - namelist_file = 'test_gvec.in' - open(newunit=unit, file=namelist_file, status='replace') - write(unit, '(A)') '&config' - write(unit, '(A)') 'trace_time = 1d-4 ! very short test run' - write(unit, '(A)') 'sbeg = 0.5d0 ! mid-radius' - write(unit, '(A)') 'ntestpart = 16 ! minimal for quick test' - write(unit, '(A)') 'field_input = "test_vmec_gvec_State_0000_00000000.dat" ! GVEC file' - write(unit, '(A)') 'isw_field_type = 0 ! canonical coordinates' - write(unit, '(A)') 'startmode = 0 ! volume sampling' - write(unit, '(A)') 'integmode = 4 ! symplectic integrator' - write(unit, '(A)') '/' - close(unit) - - ! Run SIMPLE with GVEC - error stop "Requires fully working GVEC field implementation" - call system('../../simple.x test_gvec.in > test_gvec.log') - - ! Read confined fraction - open(newunit=unit, file='confined_fraction.dat', status='old', iostat=iostat) - if (iostat /= 0) then - print *, 'ERROR: Could not read confined_fraction.dat for GVEC test' - error stop 1 - end if - - ! Skip header and read last line - read(unit, *) ! skip header - do - read(unit, *, iostat=iostat) confined_fraction_gvec - if (iostat /= 0) exit - end do - close(unit) - - print '(A,F8.4)', ' Final confined fraction: ', confined_fraction_gvec - call system('mv confined_fraction.dat confined_fraction_gvec.dat') - call system('mv times_lost.dat times_lost_gvec.dat') - - ! Compare results - print *, '' - print *, '=======================================================' - print *, 'Results Comparison:' - print *, '=======================================================' - print '(A,F8.4)', 'VMEC confined fraction: ', confined_fraction_vmec - print '(A,F8.4)', 'GVEC confined fraction: ', confined_fraction_gvec - - relative_difference = abs(confined_fraction_vmec - confined_fraction_gvec) / & - (0.5_dp * (confined_fraction_vmec + confined_fraction_gvec)) - - print '(A,F8.4,A)', 'Relative difference: ', relative_difference * 100.0_dp, '%' - - ! Clean up - call system('rm -f test_vmec.in test_gvec.in') - call system('rm -f test_vmec.log test_gvec.log') - call system('rm -f start.dat fort.6601') - - ! Check if results are consistent (allow 5% difference due to interpolation) - if (relative_difference < 0.05_dp) then - print *, '' - print *, 'TEST PASSED: VMEC and GVEC fields give consistent results' - else - print *, '' - print *, 'TEST FAILED: Results differ by more than 5%' - error stop 1 - end if - -end program test_simple_vmec_gvec diff --git a/test/tests/test_vmec_gvec.f90 b/test/tests/test_vmec_gvec.f90 deleted file mode 100644 index 31a80752..00000000 --- a/test/tests/test_vmec_gvec.f90 +++ /dev/null @@ -1,254 +0,0 @@ -program test_vmec_gvec - use, intrinsic :: iso_fortran_env, only: dp => real64 - use field_vmec, only: vmec_field_t - use field_gvec, only: gvec_field_t, create_gvec_field - use new_vmec_stuff_mod, only: netcdffile, multharm - use spline_vmec_sub, only: spline_vmec_data - use params, only: pi - - - implicit none - - class(vmec_field_t), allocatable :: vmec_field - class(gvec_field_t), allocatable :: gvec_field - - interface - subroutine export_field_2d_data(vmec_field, gvec_field) - use field_vmec, only: vmec_field_t - use field_gvec, only: gvec_field_t - class(vmec_field_t), intent(in) :: vmec_field - class(gvec_field_t), intent(in) :: gvec_field - end subroutine export_field_2d_data - - subroutine export_field_1d_data(vmec_field, gvec_field) - use field_vmec, only: vmec_field_t - use field_gvec, only: gvec_field_t - class(vmec_field_t), intent(in) :: vmec_field - class(gvec_field_t), intent(in) :: gvec_field - end subroutine export_field_1d_data - end interface - character(len=256) :: vmec_file, gvec_temp_file - logical :: file_exists - - ! Variables for GVEC conversion - character(len=256) :: param_file - character(len=4096) :: param_strings(20) - integer :: num_strings - integer :: unit_param - - ! Test parameters - real(dp) :: x(3) - real(dp) :: Acov_vmec(3), hcov_vmec(3), Bmod_vmec - real(dp) :: Acov_gvec(3), hcov_gvec(3), Bmod_gvec - real(dp) :: s_test, theta_test, phi_test - real(dp) :: rel_tol, abs_tol - integer :: ns, nt, np, i, j, k, icomp - real(dp) :: rel_error, abs_error, rel_err_B - real(dp) :: max_rel_error_Bmod, max_abs_error_Bmod - real(dp) :: max_rel_error_Acov(3), max_abs_error_Acov(3) - real(dp) :: max_rel_error_hcov(3), max_abs_error_hcov(3) - logical :: test_passed - - ! VMEC wout file downloaded by CMake (in build/test/tests directory) - vmec_file = 'wout.nc' - gvec_temp_file = 'gvec_from_vmec_wout.dat' ! GVEC state file created from VMEC - - ! Initialize VMEC field using SIMPLE's method - netcdffile = vmec_file - multharm = 3 ! Use low order splines for faster testing - call spline_vmec_data - allocate(vmec_field_t :: vmec_field) - - ! Check if CMake successfully converted the file using Python GVEC API - inquire(file=gvec_temp_file, exist=file_exists) - if (.not. file_exists) then - - ! Create a minimal GVEC parameter file for VMEC reading - param_file = 'gvec_vmec_params.ini' - open(newunit=unit_param, file=param_file, status='replace') - write(unit_param, '(A)') '! Minimal GVEC parameter file for VMEC reading' - write(unit_param, '(A)') 'whichInitEquilibrium = 1 ! Read from VMEC' - write(unit_param, '(A)') 'VMECwoutfile = ' // trim(vmec_file) - write(unit_param, '(A)') 'VMECwoutfile_format = 0 ! NetCDF format' - write(unit_param, '(A)') '' - write(unit_param, '(A)') '! Grid settings optimized for axis treatment' - write(unit_param, '(A)') 'sgrid_nelems = 21' - write(unit_param, '(A)') 'sgrid_grid_type = 4' - write(unit_param, '(A)') 'sgrid_rmin = 1.0e-6' - write(unit_param, '(A)') 'sgrid_rmax = 0.99' - write(unit_param, '(A)') '' - write(unit_param, '(A)') '! Basis settings with higher accuracy' - write(unit_param, '(A)') 'X1X2_deg = 5' - write(unit_param, '(A)') 'LA_deg = 5' - write(unit_param, '(A)') '' - write(unit_param, '(A)') '! Output control' - write(unit_param, '(A)') 'ProjectName = test_vmec_gvec' - write(unit_param, '(A)') 'outputIter = 0 ! Output initial state only' - write(unit_param, '(A)') 'maxIter = 0 ! No iterations' - close(unit_param) - - ! Use execute_command_line to run GVEC - call execute_command_line('../../_deps/gvec-build/bin/gvec ' // trim(param_file) // ' > /dev/null 2>&1', & - exitstat=num_strings) - - ! The output file will be named based on ProjectName - gvec_temp_file = 'test_vmec_gvec_State_0000_00000000.dat' - - ! Check if the state file was created - inquire(file=gvec_temp_file, exist=file_exists) - if (.not. file_exists) then - print *, 'ERROR: GVEC state file not created' - error stop 1 - end if - end if - call create_gvec_field(gvec_temp_file, gvec_field) - - ! Compare field evaluations - ! Set tolerances based on actual field comparison results - ! The excellent agreement after theta* transformations allows tight tolerances - rel_tol = 2.5e-3_dp ! 0.25% relative tolerance for |B| (actual max ~0.225%) - abs_tol = 1.5e2_dp ! 150 Gauss absolute tolerance (actual max ~136 Gauss) - - ! Grid for comparison including small s values to test axis regularization - ns = 5 ! Number of flux surfaces - nt = 4 ! Number of theta points - np = 2 ! Number of phi points - - max_rel_error_Bmod = 0.0_dp - max_abs_error_Bmod = 0.0_dp - max_rel_error_Acov = 0.0_dp - max_abs_error_Acov = 0.0_dp - max_rel_error_hcov = 0.0_dp - max_abs_error_hcov = 0.0_dp - test_passed = .true. - - do i = 1, ns - ! Test range including small s values: 0.01, 0.05, 0.1, 0.3, 0.6 - if (i == 1) then - s_test = 0.01_dp ! Very small s to test axis regularization - else if (i == 2) then - s_test = 0.05_dp ! Regularization boundary - else - s_test = 0.1_dp + 0.25_dp * real(i-3, dp) / real(ns-3, dp) ! s from 0.1 to 0.6 - end if - - do j = 1, nt - theta_test = 2.0_dp * pi * real(j-1, dp) / real(nt, dp) ! theta from 0 to 2*pi - - do k = 1, np - phi_test = pi * real(k-1, dp) / real(np, dp) ! phi from 0 to pi - - ! Set coordinates (r = sqrt(s), theta, phi) - x(1) = sqrt(s_test) - x(2) = theta_test - x(3) = phi_test - - ! Evaluate VMEC field - call vmec_field%evaluate(x, Acov_vmec, hcov_vmec, Bmod_vmec) - - ! Evaluate GVEC field - call gvec_field%evaluate(x, Acov_gvec, hcov_gvec, Bmod_gvec) - - ! Calculate errors for Bmod - abs_error = abs(Bmod_gvec - Bmod_vmec) - if (abs(Bmod_vmec) > abs_tol) then - rel_error = abs_error / abs(Bmod_vmec) - else - rel_error = 0.0_dp - end if - - max_abs_error_Bmod = max(max_abs_error_Bmod, abs_error) - max_rel_error_Bmod = max(max_rel_error_Bmod, rel_error) - - ! Calculate errors for Acov components - do icomp = 1, 3 - abs_error = abs(Acov_gvec(icomp) - Acov_vmec(icomp)) - if (abs(Acov_vmec(icomp)) > 1.0e-10_dp) then - rel_error = abs_error / abs(Acov_vmec(icomp)) - else - rel_error = 0.0_dp - end if - max_abs_error_Acov(icomp) = max(max_abs_error_Acov(icomp), abs_error) - max_rel_error_Acov(icomp) = max(max_rel_error_Acov(icomp), rel_error) - end do - - ! Calculate errors for hcov components - do icomp = 1, 3 - abs_error = abs(hcov_gvec(icomp) - hcov_vmec(icomp)) - if (abs(hcov_vmec(icomp)) > 1.0e-10_dp) then - rel_error = abs_error / abs(hcov_vmec(icomp)) - else - rel_error = 0.0_dp - end if - max_abs_error_hcov(icomp) = max(max_abs_error_hcov(icomp), abs_error) - max_rel_error_hcov(icomp) = max(max_rel_error_hcov(icomp), rel_error) - end do - - - ! Calculate relative error for |B| - rel_err_B = abs(Bmod_gvec - Bmod_vmec) / abs(Bmod_vmec) - - ! Check if errors exceed tolerance - if (max_rel_error_Bmod > rel_tol .and. max_abs_error_Bmod > abs_tol) then - test_passed = .false. - end if - end do - end do - end do - - ! Evaluate fields at a representative point for typical values - x(1) = sqrt(0.5_dp) ! r = sqrt(s) at s=0.5 - x(2) = pi/4.0_dp ! theta = pi/4 - x(3) = 0.0_dp ! phi = 0 - - call vmec_field%evaluate(x, Acov_vmec, hcov_vmec, Bmod_vmec) - call gvec_field%evaluate(x, Acov_gvec, hcov_gvec, Bmod_gvec) - - ! Print comparison results - print *, '' - print *, 'Field comparison at (s,θ,φ) = (0.5, π/4, 0):' - print *, '================================================================' - print *, ' VMEC GVEC Rel. Error' - print *, '----------------------------------------------------------------' - print '(A,ES16.8,ES16.8,ES16.8)', ' |B| ', Bmod_vmec, Bmod_gvec, & - abs(Bmod_gvec - Bmod_vmec) / abs(Bmod_vmec) - print *, '' - print '(A,ES16.8,ES16.8,ES16.8)', ' Acov(1) ', Acov_vmec(1), Acov_gvec(1), & - merge(abs(Acov_gvec(1) - Acov_vmec(1)) / abs(Acov_vmec(1)), & - 0.0_dp, abs(Acov_vmec(1)) > 1.0e-10_dp) - print '(A,ES16.8,ES16.8,ES16.8)', ' Acov(2) ', Acov_vmec(2), Acov_gvec(2), & - merge(abs(Acov_gvec(2) - Acov_vmec(2)) / abs(Acov_vmec(2)), & - 0.0_dp, abs(Acov_vmec(2)) > 1.0e-10_dp) - print '(A,ES16.8,ES16.8,ES16.8)', ' Acov(3) ', Acov_vmec(3), Acov_gvec(3), & - merge(abs(Acov_gvec(3) - Acov_vmec(3)) / abs(Acov_vmec(3)), & - 0.0_dp, abs(Acov_vmec(3)) > 1.0e-10_dp) - print *, '' - print '(A,ES16.8,ES16.8,ES16.8)', ' hcov(1) ', hcov_vmec(1), hcov_gvec(1), & - merge(abs(hcov_gvec(1) - hcov_vmec(1)) / abs(hcov_vmec(1)), & - 0.0_dp, abs(hcov_vmec(1)) > 1.0e-10_dp) - print '(A,ES16.8,ES16.8,ES16.8)', ' hcov(2) ', hcov_vmec(2), hcov_gvec(2), & - merge(abs(hcov_gvec(2) - hcov_vmec(2)) / abs(hcov_vmec(2)), & - 0.0_dp, abs(hcov_vmec(2)) > 1.0e-10_dp) - print '(A,ES16.8,ES16.8,ES16.8)', ' hcov(3) ', hcov_vmec(3), hcov_gvec(3), & - merge(abs(hcov_gvec(3) - hcov_vmec(3)) / abs(hcov_vmec(3)), & - 0.0_dp, abs(hcov_vmec(3)) > 1.0e-10_dp) - print *, '================================================================' - print *, '' - print *, 'Maximum errors over test grid:' - print *, '----------------------------------------------------------------' - print '(A,ES12.5,A,ES12.5,A)', ' |B| max error: relative = ', max_rel_error_Bmod, & - ', absolute = ', max_abs_error_Bmod, ' Gauss' - print *, '' - - ! Final test result - print *, '' - if (test_passed) then - print *, 'TEST PASSED' - else - print *, 'TEST FAILED' - print '(A,ES12.5)', ' Tolerance (relative): ', rel_tol - print '(A,ES12.5,A)', ' Tolerance (absolute): ', abs_tol, ' Gauss' - error stop 'GVEC field comparison tolerances not met' - end if - -end program test_vmec_gvec diff --git a/test/tests/test_vmec_gvec_adapter.f90 b/test/tests/test_vmec_gvec_adapter.f90 deleted file mode 100644 index a232b02e..00000000 --- a/test/tests/test_vmec_gvec_adapter.f90 +++ /dev/null @@ -1,284 +0,0 @@ -program test_vmec_gvec_adapter - !> Test that compares all adapter quantities between VMEC and GVEC fields - !> This ensures that canonical coordinates work identically for both field types - - use, intrinsic :: iso_fortran_env, only: dp => real64 - use field_vmec, only: vmec_field_t - use field_gvec, only: gvec_field_t, create_gvec_field - use vmec_field_adapter, only: vmec_field_evaluate_with_field, & - vmec_lambda_interpolate_with_field, & - vmec_iota_interpolate_with_field - use new_vmec_stuff_mod, only: netcdffile, multharm - use spline_vmec_sub, only: spline_vmec_data - use params, only: pi - - implicit none - - class(vmec_field_t), allocatable :: vmec_field - class(gvec_field_t), allocatable :: gvec_field - character(len=256) :: vmec_file, gvec_file - logical :: file_exists, test_passed - - ! Test grid - integer, parameter :: n_test = 5 - real(dp) :: s_test(n_test), theta_test(n_test), phi_test(n_test) - - ! Adapter outputs for VMEC - real(dp) :: A_theta_v, A_phi_v, dA_theta_ds_v, dA_phi_ds_v - real(dp) :: aiota_v, sqg_v, alam_v - real(dp) :: dl_ds_v, dl_dt_v, dl_dp_v - real(dp) :: Bctrvr_vartheta_v, Bctrvr_varphi_v - real(dp) :: Bcovar_r_v, Bcovar_vartheta_v, Bcovar_varphi_v - - ! Adapter outputs for GVEC - real(dp) :: A_theta_g, A_phi_g, dA_theta_ds_g, dA_phi_ds_g - real(dp) :: aiota_g, sqg_g, alam_g - real(dp) :: dl_ds_g, dl_dt_g, dl_dp_g - real(dp) :: Bctrvr_vartheta_g, Bctrvr_varphi_g - real(dp) :: Bcovar_r_g, Bcovar_vartheta_g, Bcovar_varphi_g - - ! Lambda interpolation test - real(dp) :: alam_interp_v, dl_dt_interp_v - real(dp) :: alam_interp_g, dl_dt_interp_g - - ! Iota interpolation test - real(dp) :: aiota_interp_v, daiota_ds_v - real(dp) :: aiota_interp_g, daiota_ds_g - - ! Error tracking - real(dp) :: max_rel_diff, rel_diff - real(dp) :: tol_strict = 1.0e-3_dp ! 0.1% for well-defined quantities - real(dp) :: tol_relaxed = 1.0e-2_dp ! 1% for derived quantities - real(dp) :: tol_derivatives = 5.0e-2_dp ! 5% for numerical derivatives - - ! Track what's implemented vs dummy values - logical :: test_essential_only = .true. ! Focus on essential quantities - integer :: n_pass = 0, n_fail = 0, n_dummy = 0 - - integer :: i - - ! Dummy variables for field initialization - real(dp) :: x_test(3), Acov_dummy(3), hcov_dummy(3), Bmod_dummy - - print *, '=======================================================' - print *, 'Testing vmec_field_adapter for VMEC vs GVEC' - print *, '=======================================================' - print *, '' - - ! Setup test files - vmec_file = 'wout.nc' - gvec_file = 'test_vmec_gvec_State_0000_00000000.dat' - - ! Check files exist - inquire(file=vmec_file, exist=file_exists) - if (.not. file_exists) then - print *, 'FAIL: No VMEC file found' - error stop 1 - end if - - inquire(file=gvec_file, exist=file_exists) - if (.not. file_exists) then - print *, 'FAIL: No GVEC file found. Run test_vmec_gvec first.' - error stop 1 - end if - - ! Initialize VMEC field - netcdffile = vmec_file - multharm = 5 - call spline_vmec_data - allocate(vmec_field_t :: vmec_field) - - ! Initialize GVEC field - call create_gvec_field(gvec_file, gvec_field) - - ! Define test points covering different regions - s_test = [0.1_dp, 0.3_dp, 0.5_dp, 0.7_dp, 0.9_dp] - theta_test = [0.0_dp, pi/4.0_dp, pi/2.0_dp, pi, 3.0_dp*pi/2.0_dp] - phi_test = [0.0_dp, pi/6.0_dp, pi/4.0_dp, pi/3.0_dp, pi/2.0_dp] - - test_passed = .true. - max_rel_diff = 0.0_dp - - print *, 'Testing adapter field evaluation at multiple points...' - print *, '------------------------------------------------------' - - do i = 1, n_test - print '(A,I1,A,3(F6.3,A))', 'Point ', i, ': (s,θ,φ) = (', & - s_test(i), ', ', theta_test(i)/pi, 'π, ', phi_test(i)/pi, 'π)' - - ! Evaluate with VMEC field through adapter - call vmec_field_evaluate_with_field(vmec_field, s_test(i), theta_test(i), phi_test(i), & - A_theta_v, A_phi_v, dA_theta_ds_v, dA_phi_ds_v, aiota_v, & - sqg_v, alam_v, dl_ds_v, dl_dt_v, dl_dp_v, & - Bctrvr_vartheta_v, Bctrvr_varphi_v, & - Bcovar_r_v, Bcovar_vartheta_v, Bcovar_varphi_v) - - ! Evaluate with GVEC field through adapter - ! First initialize GVEC state by calling field%evaluate - x_test = [sqrt(s_test(i)), theta_test(i), phi_test(i)] - call gvec_field%evaluate(x_test, Acov_dummy, hcov_dummy, Bmod_dummy) - - call vmec_field_evaluate_with_field(gvec_field, s_test(i), theta_test(i), phi_test(i), & - A_theta_g, A_phi_g, dA_theta_ds_g, dA_phi_ds_g, aiota_g, & - sqg_g, alam_g, dl_ds_g, dl_dt_g, dl_dp_g, & - Bctrvr_vartheta_g, Bctrvr_varphi_g, & - Bcovar_r_g, Bcovar_vartheta_g, Bcovar_varphi_g) - - ! Compare key quantities with detailed diagnostics - - ! Check sqrt(g) - note VMEC uses left-handed system (negative sqrt(g)) - call check_quantity_diagnostic(' sqrt(g)', sqg_v, sqg_g, tol_strict, 'ESSENTIAL', & - 'Jacobian - VMEC is left-handed (negative), GVEC may differ') - - ! Check iota - sign convention may differ - call check_quantity_diagnostic(' iota', aiota_v, aiota_g, tol_strict, 'ESSENTIAL', & - 'Rotational transform - sign convention may differ') - - ! Lambda and derivatives - should match if implemented - call check_quantity_diagnostic(' Lambda', alam_v, alam_g, tol_relaxed, 'COMPUTED', & - 'Stream function for canonical coordinates') - call check_quantity_diagnostic(' dLambda/ds', dl_ds_v, dl_ds_g, tol_relaxed, 'COMPUTED', & - 'Lambda radial derivative') - call check_quantity_diagnostic(' dLambda/dθ', dl_dt_v, dl_dt_g, tol_relaxed, 'COMPUTED', & - 'Lambda poloidal derivative') - - ! Magnetic field components - call check_quantity_diagnostic(' B^theta', Bctrvr_vartheta_v, Bctrvr_vartheta_g, tol_strict, 'ESSENTIAL', & - 'Contravariant poloidal B field') - call check_quantity_diagnostic(' B^phi', Bctrvr_varphi_v, Bctrvr_varphi_g, tol_strict, 'ESSENTIAL', & - 'Contravariant toroidal B field') - - ! Vector potential (gauge-dependent) - call check_quantity_diagnostic(' A_theta', A_theta_v, A_theta_g, tol_relaxed, 'COMPUTED', & - 'Poloidal vector potential - gauge dependent') - call check_quantity_diagnostic(' A_phi', A_phi_v, A_phi_g, tol_relaxed, 'COMPUTED', & - 'Toroidal vector potential') - - ! Derivatives (numerical for GVEC) - if (abs(dA_theta_ds_v) > 1.0e-10_dp .or. abs(dA_theta_ds_g) > 1.0e-10_dp) then - call check_quantity_diagnostic(' dA_theta/ds', dA_theta_ds_v, dA_theta_ds_g, tol_derivatives, 'NUMERICAL', & - 'Computed via finite differences for GVEC') - end if - if (abs(dA_phi_ds_v) > 1.0e-10_dp .or. abs(dA_phi_ds_g) > 1.0e-10_dp) then - call check_quantity_diagnostic(' dA_phi/ds', dA_phi_ds_v, dA_phi_ds_g, tol_derivatives, 'NUMERICAL', & - 'Computed via finite differences for GVEC') - end if - - print *, '' - end do - - ! Test lambda interpolation consistency - print *, 'Testing lambda interpolation consistency...' - print *, '------------------------------------------' - do i = 1, n_test - call vmec_lambda_interpolate_with_field(vmec_field, s_test(i), theta_test(i), phi_test(i), & - alam_interp_v, dl_dt_interp_v) - call vmec_lambda_interpolate_with_field(gvec_field, s_test(i), theta_test(i), phi_test(i), & - alam_interp_g, dl_dt_interp_g) - - call check_quantity_diagnostic('Lambda interp', alam_interp_v, alam_interp_g, tol_relaxed, 'INTERFACE', & - 'Lambda interpolation interface') - call check_quantity_diagnostic('dLambda/dθ interp', dl_dt_interp_v, dl_dt_interp_g, tol_relaxed, 'INTERFACE', & - 'Lambda derivative interpolation') - end do - print *, '' - - ! Test iota interpolation consistency - print *, 'Testing iota interpolation consistency...' - print *, '----------------------------------------' - do i = 1, n_test - call vmec_iota_interpolate_with_field(vmec_field, s_test(i), aiota_interp_v, daiota_ds_v) - call vmec_iota_interpolate_with_field(gvec_field, s_test(i), aiota_interp_g, daiota_ds_g) - - call check_quantity_diagnostic('iota interp', aiota_interp_v, aiota_interp_g, tol_strict, 'INTERFACE', & - 'Iota interpolation - sign convention issue') - call check_quantity_diagnostic('diota/ds', daiota_ds_v, daiota_ds_g, tol_relaxed, 'INTERFACE', & - 'Iota derivative - sign convention issue') - end do - - ! Summary - print *, '' - print *, '=======================================================' - print *, 'TEST SUMMARY - VMEC vs GVEC Adapter Comparison' - print *, '=======================================================' - print '(A,I4)', 'Total tests performed: ', n_pass + n_fail + n_dummy - print '(A,I4,A)', 'PASSED: ', n_pass, ' (within tolerance)' - print '(A,I4,A)', 'FAILED: ', n_fail, ' (excessive difference)' - print '(A,I4,A)', 'DUMMY: ', n_dummy, ' (placeholder values)' - print *, '' - print '(A,ES12.5)', 'Maximum relative difference: ', max_rel_diff - print *, '' - print *, 'KNOWN ISSUES:' - print *, '1. Sign conventions: iota and sqrt(g) have opposite signs' - print *, '2. Coordinate systems: VMEC is left-handed, GVEC may differ' - print *, '3. Some quantities not fully implemented in GVEC adapter' - print *, '' - - if (n_fail == 0) then - print *, 'STATUS: All implemented features working correctly' - else if (test_essential_only .and. n_pass > n_fail) then - print *, 'STATUS: Essential features partially working' - print *, ' Further development needed for full compatibility' - else - print *, 'STATUS: Significant issues found' - print *, ' Canonical coordinates may not work correctly with GVEC' - end if - -contains - - subroutine check_quantity_diagnostic(name, val_vmec, val_gvec, tol, category, description) - character(*), intent(in) :: name, category, description - real(dp), intent(in) :: val_vmec, val_gvec, tol - - real(dp) :: rel_diff, abs_diff - character(len=8) :: status - logical :: is_dummy - - abs_diff = abs(val_vmec - val_gvec) - - ! Check for dummy/placeholder values - is_dummy = (abs(val_vmec) < 1.0e-10_dp .and. category == 'COMPUTED') .or. & - (category == 'NUMERICAL' .and. abs(val_gvec) < 1.0e-10_dp) - - ! Handle relative difference calculation - if (abs(val_vmec) > 1.0e-10_dp) then - rel_diff = abs_diff / abs(val_vmec) - else if (abs(val_gvec) > 1.0e-10_dp) then - rel_diff = abs_diff / abs(val_gvec) - else - ! Both values are near zero - rel_diff = 0.0_dp - end if - - if (rel_diff > max_rel_diff) max_rel_diff = rel_diff - - ! Determine status - if (is_dummy) then - status = 'DUMMY' - n_dummy = n_dummy + 1 - else if (rel_diff <= tol) then - status = 'PASS' - n_pass = n_pass + 1 - else - status = 'FAIL' - n_fail = n_fail + 1 - end if - - ! Print diagnostic info - if (status /= 'PASS' .or. .true.) then ! Set to .true. for verbose output - print '(A,A8,A,ES12.5,A,ES12.5,A,ES12.5,A,A)', & - name, status, ': V=', val_vmec, ', G=', val_gvec, & - ', diff=', rel_diff, ' [', trim(category), ']' - if (status == 'FAIL') then - print '(A,A)', ' Note: ', trim(description) - end if - end if - - ! Update test_passed for backward compatibility - if (status == 'FAIL' .and. category == 'ESSENTIAL') then - test_passed = .false. - end if - end subroutine check_quantity_diagnostic - - -end program test_vmec_gvec_adapter \ No newline at end of file diff --git a/thirdparty/LICENSE.GVEC b/thirdparty/LICENSE.GVEC deleted file mode 100644 index 6abbf7ae..00000000 --- a/thirdparty/LICENSE.GVEC +++ /dev/null @@ -1,9 +0,0 @@ -MIT License - -Copyright (c) 2025 GVEC Contributors, Max Planck Institute for Plasma Physics - -Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the " Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice (including the next paragraph) shall be included in all copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. From 7c63e772581d6551bee11eb6719dc933f076c251 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Fri, 27 Mar 2026 17:40:45 +0100 Subject: [PATCH 2/2] test: align optional Python test registration --- python/CMakeLists.txt | 4 ++++ test/CMakeLists.txt | 4 +++- test/tests/CMakeLists.txt | 15 ++++++++++++++- test/tests/magfie/CMakeLists.txt | 2 +- 4 files changed, 22 insertions(+), 3 deletions(-) diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 5b3d26ca..9545f6a5 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -1,5 +1,7 @@ cmake_minimum_required(VERSION 3.22) +set(SIMPLE_PYTHON_INTERFACE_BUILT OFF PARENT_SCOPE) + # Grab Python find_package(Python COMPONENTS Interpreter Development NumPy) @@ -220,3 +222,5 @@ add_custom_target( install(TARGETS ${F2PY_MODULE_NAME} DESTINATION .) install(FILES ${CMAKE_BINARY_DIR}/simple_backend.py DESTINATION .) + +set(SIMPLE_PYTHON_INTERFACE_BUILT ON PARENT_SCOPE) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 1b9c1ced..a2af8e98 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -1,4 +1,6 @@ add_subdirectory(tests) -if(ENABLE_PYTHON_INTERFACE) +if(ENABLE_PYTHON_INTERFACE AND SIMPLE_PYTHON_INTERFACE_BUILT) add_subdirectory(python) +elseif(ENABLE_PYTHON_INTERFACE) + message(STATUS "Skipping Python tests because the SIMPLE Python backend was not built") endif() diff --git a/test/tests/CMakeLists.txt b/test/tests/CMakeLists.txt index 89d6a346..d249240c 100644 --- a/test/tests/CMakeLists.txt +++ b/test/tests/CMakeLists.txt @@ -31,8 +31,21 @@ add_custom_command( add_custom_target(chartmap_test_data DEPENDS ${CHARTMAP_FILE}) # map2disc-based chartmap generation (via libneo Python API) +set(SIMPLE_HAVE_MAP2DISC OFF) if(SIMPLE_ENABLE_PYTHON_TOOLS) find_package(Python3 REQUIRED COMPONENTS Interpreter) + execute_process( + COMMAND ${Python3_EXECUTABLE} -c "import importlib.util, sys; sys.exit(0 if importlib.util.find_spec('map2disc') else 1)" + RESULT_VARIABLE MAP2DISC_IMPORT_RESULT + ) + if(MAP2DISC_IMPORT_RESULT EQUAL 0) + set(SIMPLE_HAVE_MAP2DISC ON) + else() + message(STATUS "map2disc Python package not found; skipping map2disc-based chartmap tests") + endif() +endif() + +if(SIMPLE_ENABLE_PYTHON_TOOLS AND SIMPLE_HAVE_MAP2DISC) set(CHARTMAP_MAP2DISC_FILE "${CMAKE_CURRENT_BINARY_DIR}/wout.chartmap.map2disc.nc") set(CHARTMAP_NCSX_MAP2DISC_FILE "${CMAKE_CURRENT_BINARY_DIR}/wout_ncsx.chartmap.map2disc.nc") set(CHARTMAP_NCSX_MAP2DISC_WALL_FILE "${CMAKE_CURRENT_BINARY_DIR}/wout_ncsx.chartmap.map2disc.wall0p1m.nc") @@ -411,7 +424,7 @@ set_tests_properties(test_field_equivalence_chartmap PROPERTIES LABELS "integration;slow" TIMEOUT 300) -if(SIMPLE_ENABLE_PYTHON_TOOLS) +if(SIMPLE_ENABLE_PYTHON_TOOLS AND SIMPLE_HAVE_MAP2DISC) add_test(NAME test_field_equivalence_chartmap_map2disc COMMAND test_field_equivalence_chartmap.x wout_ncsx.chartmap.map2disc.nc field_equiv_map2disc WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) diff --git a/test/tests/magfie/CMakeLists.txt b/test/tests/magfie/CMakeLists.txt index fe8d560d..66d473ae 100644 --- a/test/tests/magfie/CMakeLists.txt +++ b/test/tests/magfie/CMakeLists.txt @@ -138,7 +138,7 @@ if(SIMPLE_ENABLE_PYTHON_TOOLS) ) endif() -if(SIMPLE_ENABLE_PYTHON_TOOLS) +if(SIMPLE_ENABLE_PYTHON_TOOLS AND SIMPLE_HAVE_MAP2DISC) add_test(NAME test_orbit_chartmap_comparison_map2disc COMMAND $ wout.chartmap.map2disc.nc orbit_chartmap_comparison_map2disc.nc