From 42bdec2af90be7539ea2037c71937e414b77c4fc Mon Sep 17 00:00:00 2001 From: Markus Markl Date: Mon, 2 Mar 2026 19:21:08 +0100 Subject: [PATCH 01/20] build(QL-Balance): link KIM_lib into ql-balance_lib Add KIM_lib as a dependency and link target for ql-balance_lib to enable the upcoming KIM wave code adapter integration. To resolve duplicate sparse_mod symbols between KIM's bundled copy and the shared sparse target, KIM now links against the shared sparse target instead of bundling its own copy. Changes: - QL-Balance/CMakeLists.txt: add KIM_lib to add_dependencies and target_link_libraries - KIM/src/CMakeSources.in: remove bundled sparse_mod.f90 - KIM/src/CMakeLists.txt: link KIM_lib against shared sparse target Co-Authored-By: Claude Opus 4.6 --- KIM/src/CMakeLists.txt | 5 +++++ KIM/src/CMakeSources.in | 1 - QL-Balance/CMakeLists.txt | 7 ++++++- 3 files changed, 11 insertions(+), 2 deletions(-) diff --git a/KIM/src/CMakeLists.txt b/KIM/src/CMakeLists.txt index 8b8b54a5..4a95d54a 100644 --- a/KIM/src/CMakeLists.txt +++ b/KIM/src/CMakeLists.txt @@ -23,6 +23,11 @@ set_target_properties(KIM_lib PROPERTIES OUTPUT_NAME KIM ARCHIVE_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/install/lib/") +# Use the shared sparse target instead of bundling a local copy +add_dependencies(KIM_lib sparse) +target_link_libraries(KIM_lib PUBLIC sparse) +target_include_directories(KIM_lib PUBLIC "${CMAKE_BINARY_DIR}/OBJS/sparse/") + # TODO: unify this: either use the system library everywhere, or use the cloned version # KiLCA uses the cloned libraries for now, QL-Balance the system libraries find_package(GSL REQUIRED) diff --git a/KIM/src/CMakeSources.in b/KIM/src/CMakeSources.in index 24b97f61..cce8e65b 100644 --- a/KIM/src/CMakeSources.in +++ b/KIM/src/CMakeSources.in @@ -17,7 +17,6 @@ set(KIM_math math/suitesparse/umf4_f77zwrapper.c math/superlu/c_fortran_zgssv.c math/gsl_mod.f90 math/use_libcerf_mod.f90 - math/sparse_mod.f90 math/plasma_disp_func.f90 math/functions.f90 math/numerics_utils.f90 diff --git a/QL-Balance/CMakeLists.txt b/QL-Balance/CMakeLists.txt index 24a7703c..e77a3c6d 100644 --- a/QL-Balance/CMakeLists.txt +++ b/QL-Balance/CMakeLists.txt @@ -108,8 +108,12 @@ set_target_properties(ql-balance_lib PROPERTIES LINKER_LANGUAGE Fortran ARCHIVE_OUTPUT_DIRECTORY "${lib}" Fortran_MODULE_DIRECTORY "${objs}") -add_dependencies(ql-balance_lib kilca_lib kamel_hdf5_tools) +add_dependencies(ql-balance_lib kilca_lib kamel_hdf5_tools KIM_lib) target_include_directories(ql-balance_lib PUBLIC ${CMAKE_BINARY_DIR}/KiLCA/mod ${CMAKE_BINARY_DIR}/modules) +# NOTE: KIM Fortran modules (OBJS/kim/) are NOT added here because KIM and +# QL-Balance both define a module named resonances_mod with different contents. +# The adapter source file that needs KIM modules will use a per-file include +# (set_source_files_properties ... INCLUDE_DIRECTORIES) to resolve this. # Include central slatec math sources via cache variable defined by common/math target_include_directories(ql-balance_lib PUBLIC ${SLATEC_INCLUDE_DIR}) @@ -147,6 +151,7 @@ target_link_libraries(ql-balance_lib PUBLIC hdf5::hdf5_hl hdf5::hdf5_hl_fortran kamel_hdf5_tools + KIM_lib kilca_lib lapack netcdf From 11fd8eee9d37d5c98a5cb91fc4b20c741513c636 Mon Sep 17 00:00:00 2001 From: Markus Markl Date: Mon, 2 Mar 2026 20:02:18 +0100 Subject: [PATCH 02/20] feat(QL-Balance): add wave_code selector to balance_conf.nml Adds 'wave_code' namelist variable to BALANCENML with default 'KiLCA'. Supports future selection between KiLCA and KIM wave code backends. Co-Authored-By: Claude Opus 4.6 --- QL-Balance/namelists/balance_conf.nml | 1 + QL-Balance/src/base/control_mod.f90 | 1 + QL-Balance/src/base/read_config.f90 | 7 ++++--- 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/QL-Balance/namelists/balance_conf.nml b/QL-Balance/namelists/balance_conf.nml index 62c9d229..2d148a6b 100644 --- a/QL-Balance/namelists/balance_conf.nml +++ b/QL-Balance/namelists/balance_conf.nml @@ -55,5 +55,6 @@ misalign_diffusion = .false. , equil_path = '/equil_r_q_psi.dat' , type_of_run = 'SingleStep' + wave_code = 'KiLCA' ihdf5IO = 1 / diff --git a/QL-Balance/src/base/control_mod.f90 b/QL-Balance/src/base/control_mod.f90 index 60dcf9cb..67afbe1c 100644 --- a/QL-Balance/src/base/control_mod.f90 +++ b/QL-Balance/src/base/control_mod.f90 @@ -5,6 +5,7 @@ module control_mod implicit none character(100) :: type_of_run + character(len=32) :: wave_code = 'KiLCA' ! 'KiLCA' or 'KIM' integer :: ihdf5IO ! added: Markus Markl logical :: paramscan ! added: Markus Markl, 03.03.2021 logical :: timing_mode ! added by Markus Markl 06.04.2021 diff --git a/QL-Balance/src/base/read_config.f90 b/QL-Balance/src/base/read_config.f90 index d65c3c23..64e3f49a 100644 --- a/QL-Balance/src/base/read_config.f90 +++ b/QL-Balance/src/base/read_config.f90 @@ -2,7 +2,7 @@ subroutine read_config use baseparam_mod, only: btor, rtor, rsepar, dperp, Z_i, am, urelax use control_mod, only: eps, paramscan, diagnostics_output, suppression_mode, debug_mode, & readfromtimestep, temperature_limit, gyro_current_study, & - misalign_diffusion, equil_path, ihdf5IO + misalign_diffusion, equil_path, ihdf5IO, wave_code use grid_mod, only: rmin, rmax, npoimin, gg_factor, gg_width, gg_r_res, iboutype, rb_cut_in, & re_cut_in, rb_cut_out, re_cut_out use h5mod, only: path2inp, path2out, path2time @@ -23,8 +23,8 @@ subroutine read_config timstep_min, paramscan, save_prof_time_step, diagnostics_output, br_stopping, & suppression_mode, debug_mode, readfromtimestep, path2time, ramp_up_mode, t_max_ramp_up, & temperature_limit, antenna_max_stopping, gyro_current_study, viscosity_factor, & - misalign_diffusion, equil_path, ihdf5IO, type_of_run, set_constant_time_step, & - constant_time_step, urelax + misalign_diffusion, equil_path, ihdf5IO, type_of_run, wave_code, & + set_constant_time_step, constant_time_step, urelax ! read the parameters from namelist file open (newunit=u, file=config_file, status="old", action="read", iostat=ios) @@ -36,6 +36,7 @@ subroutine read_config write (*, *) "" write (*, *) "=================================================================================" write (*, "(A,A)") " Type of Run: ", type_of_run + write (*, "(A,A)") " Wave Code: ", trim(adjustl(wave_code)) write (*, *) "" write (*, "(A,A,A)") " Parameters from ", config_file, ":" write (*, *) " ------------------------------------------------------------------------------" From 20440f6c9daefa5baa2aba09681d187932ae74f6 Mon Sep 17 00:00:00 2001 From: Markus Markl Date: Mon, 2 Mar 2026 20:04:19 +0100 Subject: [PATCH 03/20] feat(QL-Balance): add KIM wave code adapter module (skeleton) Skeleton module with public API stubs for KIM-QL-Balance interface: kim_initialize, kim_run_for_all_modes, kim_update_profiles, kim_get_wave_fields, kim_get_wave_vectors, kim_get_background_magnetic_fields, kim_get_collision_frequencies. Co-Authored-By: Claude Opus 4.6 --- QL-Balance/CMakeLists.txt | 3 + QL-Balance/src/base/kim_wave_code_adapter.f90 | 58 +++++++++++++++++++ 2 files changed, 61 insertions(+) create mode 100644 QL-Balance/src/base/kim_wave_code_adapter.f90 diff --git a/QL-Balance/CMakeLists.txt b/QL-Balance/CMakeLists.txt index e77a3c6d..edcff841 100644 --- a/QL-Balance/CMakeLists.txt +++ b/QL-Balance/CMakeLists.txt @@ -62,6 +62,7 @@ set(balance_lib_base_fortran_sources ${base_balance_src}/grid_mod.f90 ${base_balance_src}/h5mod.f90 ${base_balance_src}/integrate_parallel_current.f90 + ${base_balance_src}/kim_wave_code_adapter.f90 ${base_balance_src}/magnetic_island_width.f90 ${base_balance_src}/matrix_mod.f90 ${base_balance_src}/paramscan.f90 @@ -114,6 +115,8 @@ target_include_directories(ql-balance_lib PUBLIC ${CMAKE_BINARY_DIR}/KiLCA/mod $ # QL-Balance both define a module named resonances_mod with different contents. # The adapter source file that needs KIM modules will use a per-file include # (set_source_files_properties ... INCLUDE_DIRECTORIES) to resolve this. +# kim_wave_code_adapter.f90 will need per-file INCLUDE_DIRECTORIES for KIM modules +# (see note at target_include_directories about resonances_mod collision) # Include central slatec math sources via cache variable defined by common/math target_include_directories(ql-balance_lib PUBLIC ${SLATEC_INCLUDE_DIR}) diff --git a/QL-Balance/src/base/kim_wave_code_adapter.f90 b/QL-Balance/src/base/kim_wave_code_adapter.f90 new file mode 100644 index 00000000..f816825f --- /dev/null +++ b/QL-Balance/src/base/kim_wave_code_adapter.f90 @@ -0,0 +1,58 @@ +module kim_wave_code_adapter_m + !! Adapter module that wraps KIM library calls to populate + !! wave_code_data module arrays, matching the contract expected + !! by get_dql.f90 and diff_coeffs.f90. + + implicit none + private + + public :: kim_initialize + public :: kim_run_for_all_modes + public :: kim_update_profiles + public :: kim_get_wave_fields + public :: kim_get_wave_vectors + public :: kim_get_background_magnetic_fields + public :: kim_get_collision_frequencies + +contains + + subroutine kim_initialize(nrad, r_grid) + !! Initialize KIM backend: read config, profiles, grids. + !! Populate wave_code_data module arrays with background quantities. + integer, intent(in) :: nrad + real(8), intent(in) :: r_grid(nrad) + ! TODO: implement in Task 4 + end subroutine + + subroutine kim_run_for_all_modes() + !! Run KIM solver for all (m,n) modes. + ! TODO: implement in Task 5 + end subroutine + + subroutine kim_update_profiles() + !! Update KIM's internal profiles from QL-Balance's evolved parameters. + ! TODO: implement in Task 5 + end subroutine + + subroutine kim_get_wave_fields(i_mn) + !! Extract wave fields from KIM's EBdat into wave_code_data arrays. + integer, intent(in) :: i_mn + ! TODO: implement in Task 6 + end subroutine + + subroutine kim_get_wave_vectors() + !! Extract wave vectors from KIM's plasma into wave_code_data arrays. + ! TODO: implement in Task 6 + end subroutine + + subroutine kim_get_background_magnetic_fields() + !! Extract background B from KIM's equilibrium into wave_code_data arrays. + ! TODO: implement in Task 6 + end subroutine + + subroutine kim_get_collision_frequencies() + !! Extract collision frequencies from KIM's plasma into wave_code_data arrays. + ! TODO: implement in Task 6 + end subroutine + +end module kim_wave_code_adapter_m From 62e1ee4f9b0167950ff3e11a6b39f91e5e625be1 Mon Sep 17 00:00:00 2001 From: Markus Markl Date: Mon, 2 Mar 2026 20:11:06 +0100 Subject: [PATCH 04/20] feat(QL-Balance): implement kim_initialize with background extraction Initializes KIM backend, runs electrostatic solver for first mode, and extracts background quantities (B0, kp, ks, nue, nui, dPhi0, profiles) from KIM modules onto the QL-Balance radial grid. Co-Authored-By: Claude Opus 4.6 --- QL-Balance/CMakeLists.txt | 7 +- QL-Balance/src/base/kim_wave_code_adapter.f90 | 201 +++++++++++++++++- 2 files changed, 198 insertions(+), 10 deletions(-) diff --git a/QL-Balance/CMakeLists.txt b/QL-Balance/CMakeLists.txt index edcff841..e899a1bb 100644 --- a/QL-Balance/CMakeLists.txt +++ b/QL-Balance/CMakeLists.txt @@ -115,8 +115,11 @@ target_include_directories(ql-balance_lib PUBLIC ${CMAKE_BINARY_DIR}/KiLCA/mod $ # QL-Balance both define a module named resonances_mod with different contents. # The adapter source file that needs KIM modules will use a per-file include # (set_source_files_properties ... INCLUDE_DIRECTORIES) to resolve this. -# kim_wave_code_adapter.f90 will need per-file INCLUDE_DIRECTORIES for KIM modules -# (see note at target_include_directories about resonances_mod collision) +# Give ONLY the KIM adapter file access to KIM's Fortran modules so it can +# import species_m, equilibrium_m, config_m, etc. Every other QL-Balance +# source keeps seeing QL-Balance's own resonances_mod. +set_source_files_properties(${base_balance_src}/kim_wave_code_adapter.f90 + PROPERTIES INCLUDE_DIRECTORIES "${CMAKE_BINARY_DIR}/OBJS/kim/") # Include central slatec math sources via cache variable defined by common/math target_include_directories(ql-balance_lib PUBLIC ${SLATEC_INCLUDE_DIR}) diff --git a/QL-Balance/src/base/kim_wave_code_adapter.f90 b/QL-Balance/src/base/kim_wave_code_adapter.f90 index f816825f..6bec6dc4 100644 --- a/QL-Balance/src/base/kim_wave_code_adapter.f90 +++ b/QL-Balance/src/base/kim_wave_code_adapter.f90 @@ -2,6 +2,20 @@ module kim_wave_code_adapter_m !! Adapter module that wraps KIM library calls to populate !! wave_code_data module arrays, matching the contract expected !! by get_dql.f90 and diff_coeffs.f90. + !! + !! KIM modules are accessed via per-file INCLUDE_DIRECTORIES in + !! CMake to avoid collision with the QL-Balance resonances_mod. + !! All KIM symbols are renamed on import to avoid name clashes. + + use species_m, only: kim_plasma => plasma + use equilibrium_m, only: kim_B0 => B0 + use config_m, only: kim_type_of_run => type_of_run, & + nml_config_path + use setup_m, only: kim_m_mode => m_mode, kim_n_mode => n_mode + use grid_m, only: kim_xl_grid => xl_grid, kim_rg_grid => rg_grid + use fields_m, only: EBdat + use kim_base_m, only: kim_t + use kim_mod_m, only: from_kim_factory_get_kim implicit none private @@ -14,15 +28,181 @@ module kim_wave_code_adapter_m public :: kim_get_background_magnetic_fields public :: kim_get_collision_frequencies + !! Module-level KIM solver instance (reused across calls) + class(kim_t), allocatable :: kim_instance + + !! Vacuum Br placeholder (filled properly in Task 7) + complex(8), allocatable, public :: kim_vac_Br(:) + contains subroutine kim_initialize(nrad, r_grid) !! Initialize KIM backend: read config, profiles, grids. - !! Populate wave_code_data module arrays with background quantities. + !! Populate wave_code_data module arrays with background + !! quantities by running the electrostatic solver for the + !! first mode to generate equilibrium and plasma backgrounds. + use wave_code_data, only: dim_mn, m_vals, n_vals, & + r => r, q => q, n => n, Te => Te, Ti => Ti, & + Vth => Vth, Vz => Vz, dPhi0 => dPhi0, & + kp => kp, ks => ks, om_E => om_E, & + B0 => B0, nue => nue, nui => nui, & + flre_path + + implicit none + integer, intent(in) :: nrad real(8), intent(in) :: r_grid(nrad) - ! TODO: implement in Task 4 - end subroutine + + integer :: kim_npts + real(8), allocatable :: kim_r(:) + real(8), allocatable :: work_old(:), work_new(:) + + ! ----------------------------------------------------------- + ! 1. Set KIM config path and initialize KIM backend + ! ----------------------------------------------------------- + ! nml_config_path (from config_m) tells KIM where its + ! namelist file lives. Default is ./KIM_config.nml which + ! the user can override via Task 8 (balance namelist). + nml_config_path = "./KIM_config.nml" + + call kim_init() + + ! ----------------------------------------------------------- + ! 2. Read (m,n) modes from modes.in via existing utility + ! ----------------------------------------------------------- + call read_antenna_modes(flre_path) + + ! ----------------------------------------------------------- + ! 3. Allocate wave_code_data arrays on QL-Balance grid + ! ----------------------------------------------------------- + call allocate_wave_code_data(nrad, r_grid) + + ! ----------------------------------------------------------- + ! 4. Run KIM electrostatic solver for the first mode to + ! generate equilibrium, grids, and plasma backgrounds + ! ----------------------------------------------------------- + kim_type_of_run = 'electrostatic' + kim_m_mode = m_vals(1) + kim_n_mode = n_vals(1) + + call from_kim_factory_get_kim(kim_type_of_run, kim_instance) + call kim_instance%init() + call kim_instance%run() + + ! ----------------------------------------------------------- + ! 5. Extract and interpolate KIM backgrounds onto the + ! QL-Balance radial grid + ! ----------------------------------------------------------- + kim_npts = kim_plasma%grid_size + allocate(kim_r(kim_npts)) + kim_r = kim_plasma%r_grid + + ! Temporary work arrays for interpolation + allocate(work_old(kim_npts), work_new(nrad)) + + ! Safety factor q + work_old = kim_plasma%q + call interp_profile(kim_npts, kim_r, work_old, nrad, r, q) + + ! Parallel wavenumber kp + work_old = kim_plasma%kp + call interp_profile(kim_npts, kim_r, work_old, nrad, r, kp) + + ! Perpendicular wavenumber ks + work_old = kim_plasma%ks + call interp_profile(kim_npts, kim_r, work_old, nrad, r, ks) + + ! ExB drift frequency om_E + work_old = kim_plasma%om_E + call interp_profile(kim_npts, kim_r, work_old, nrad, r, om_E) + + ! Equilibrium radial electric field -> dPhi0 + ! Sign convention: dPhi0 = -Er + work_old = -kim_plasma%Er + call interp_profile(kim_npts, kim_r, work_old, nrad, r, dPhi0) + + ! Electron density [1/cm^3] + work_old = kim_plasma%spec(0)%n + call interp_profile(kim_npts, kim_r, work_old, nrad, r, n) + + ! Electron temperature [eV] + work_old = kim_plasma%spec(0)%T + call interp_profile(kim_npts, kim_r, work_old, nrad, r, Te) + + ! Ion temperature [eV] + work_old = kim_plasma%spec(1)%T + call interp_profile(kim_npts, kim_r, work_old, nrad, r, Ti) + + ! Electron collision frequency + work_old = kim_plasma%spec(0)%nu + call interp_profile(kim_npts, kim_r, work_old, nrad, r, nue) + + ! Ion collision frequency + work_old = kim_plasma%spec(1)%nu + call interp_profile(kim_npts, kim_r, work_old, nrad, r, nui) + + ! Background magnetic field B0 (total) + ! equilibrium_m::B0 is on the plasma r_grid + work_old = kim_B0 + call interp_profile(kim_npts, kim_r, work_old, nrad, r, B0) + + ! B0 toroidal and poloidal components from equilibrium_m + call interp_B0_components(kim_npts, kim_r, nrad, r) + + ! Poloidal rotation velocity Vth + ! KIM does not store Vth; set to zero (can be updated later) + Vth = 0.0d0 + + ! Toroidal rotation Vz + ! KIM does not store Vz directly in plasma struct; + ! set to zero (Er already accounts for force balance) + Vz = 0.0d0 + + ! ----------------------------------------------------------- + ! 6. Set up vacuum Br placeholder for form factor + ! (filled properly in Task 7) + ! ----------------------------------------------------------- + if (allocated(kim_vac_Br)) deallocate(kim_vac_Br) + allocate(kim_vac_Br(nrad)) + kim_vac_Br = (1.0d0, 0.0d0) ! placeholder: unit vacuum Br + + deallocate(kim_r, work_old, work_new) + + write(*, *) "KIM adapter: initialization complete" + write(*, *) " Balance grid points: ", nrad + write(*, *) " KIM grid points: ", kim_npts + write(*, *) " Number of modes: ", dim_mn + + end subroutine kim_initialize + + ! --------------------------------------------------------------- + ! Helper: interpolate B0 toroidal and poloidal components + ! --------------------------------------------------------------- + subroutine interp_B0_components(kim_npts, kim_r, nbal, bal_r) + !! Interpolate B0z (toroidal) and B0th (poloidal) from + !! equilibrium_m onto the QL-Balance grid. + use equilibrium_m, only: B0z_equil => B0z, B0th_equil => B0th + use wave_code_data, only: B0t, B0z + + implicit none + + integer, intent(in) :: kim_npts, nbal + real(8), intent(in) :: kim_r(kim_npts), bal_r(nbal) + real(8), allocatable :: work(:) + + allocate(work(kim_npts)) + + ! B0z (toroidal component) -> wave_code_data::B0z + work = B0z_equil + call interp_profile(kim_npts, kim_r, work, nbal, bal_r, B0z) + + ! B0th (poloidal component) -> wave_code_data::B0t + work = B0th_equil + call interp_profile(kim_npts, kim_r, work, nbal, bal_r, B0t) + + deallocate(work) + + end subroutine interp_B0_components subroutine kim_run_for_all_modes() !! Run KIM solver for all (m,n) modes. @@ -30,28 +210,33 @@ subroutine kim_run_for_all_modes() end subroutine subroutine kim_update_profiles() - !! Update KIM's internal profiles from QL-Balance's evolved parameters. + !! Update KIM internal profiles from QL-Balance evolved + !! parameters. ! TODO: implement in Task 5 end subroutine subroutine kim_get_wave_fields(i_mn) - !! Extract wave fields from KIM's EBdat into wave_code_data arrays. + !! Extract wave fields from KIM EBdat into wave_code_data + !! arrays. integer, intent(in) :: i_mn ! TODO: implement in Task 6 end subroutine subroutine kim_get_wave_vectors() - !! Extract wave vectors from KIM's plasma into wave_code_data arrays. + !! Extract wave vectors from KIM plasma into wave_code_data + !! arrays. ! TODO: implement in Task 6 end subroutine subroutine kim_get_background_magnetic_fields() - !! Extract background B from KIM's equilibrium into wave_code_data arrays. + !! Extract background B from KIM equilibrium into + !! wave_code_data arrays. ! TODO: implement in Task 6 end subroutine subroutine kim_get_collision_frequencies() - !! Extract collision frequencies from KIM's plasma into wave_code_data arrays. + !! Extract collision frequencies from KIM plasma into + !! wave_code_data arrays. ! TODO: implement in Task 6 end subroutine From 42e0ad58b918c8d85fa94199308889478cdd8d85 Mon Sep 17 00:00:00 2001 From: Markus Markl Date: Mon, 2 Mar 2026 20:19:00 +0100 Subject: [PATCH 05/20] feat(QL-Balance): implement kim_run_for_all_modes with per-mode solve Runs KIM electrostatic solver for each (m,n) mode and stores Es/Br results in per-mode arrays. Adds kim_update_profiles for transferring QL-Balance time-evolved parameters back to wave_code_data arrays. Key additions: - kim_run_for_all_modes: loops over modes, sets m/n, creates fresh KIM solver via factory, runs Poisson solve, interpolates fields onto the QL-Balance grid, and stores results in kim_*_modes arrays - kim_update_profiles: mirrors update_background_files for KIM path, converting params_b (erg) to wave_code_data units (eV, cm/s) - deallocate_EBdat helper: cleans up global EBdat between mode solves - interp_complex_profile helper: interpolates complex profiles by splitting into real/imaginary parts Co-Authored-By: Claude Opus 4.6 --- QL-Balance/src/base/kim_wave_code_adapter.f90 | 207 +++++++++++++++++- 1 file changed, 200 insertions(+), 7 deletions(-) diff --git a/QL-Balance/src/base/kim_wave_code_adapter.f90 b/QL-Balance/src/base/kim_wave_code_adapter.f90 index 6bec6dc4..af7f9ba9 100644 --- a/QL-Balance/src/base/kim_wave_code_adapter.f90 +++ b/QL-Balance/src/base/kim_wave_code_adapter.f90 @@ -34,6 +34,15 @@ module kim_wave_code_adapter_m !! Vacuum Br placeholder (filled properly in Task 7) complex(8), allocatable, public :: kim_vac_Br(:) + !! Per-mode stored field results (nrad, dim_mn) + !! Filled by kim_run_for_all_modes, read by kim_get_wave_fields + complex(8), allocatable :: kim_Es_modes(:,:) + complex(8), allocatable :: kim_Ep_modes(:,:) + complex(8), allocatable :: kim_Er_modes(:,:) + complex(8), allocatable :: kim_Et_modes(:,:) + complex(8), allocatable :: kim_Ez_modes(:,:) + complex(8), allocatable :: kim_Br_modes(:,:) + contains subroutine kim_initialize(nrad, r_grid) @@ -205,15 +214,140 @@ subroutine interp_B0_components(kim_npts, kim_r, nbal, bal_r) end subroutine interp_B0_components subroutine kim_run_for_all_modes() - !! Run KIM solver for all (m,n) modes. - ! TODO: implement in Task 5 - end subroutine + !! Run KIM electrostatic solver for each (m,n) mode and + !! store the resulting fields in per-mode arrays. + !! + !! After this call, kim_Es_modes(:,i_mn) etc. hold the + !! field solutions interpolated onto the QL-Balance radial + !! grid. kim_get_wave_fields(i_mn) copies from these + !! arrays into the wave_code_data module scalars. + use wave_code_data, only: dim_mn, m_vals, n_vals, & + dim_r, bal_r => r + + implicit none + + integer :: i_mn, kim_npts + real(8), allocatable :: kim_r(:) + + ! ------------------------------------------------------- + ! 1. (Re-)allocate per-mode storage + ! ------------------------------------------------------- + if (allocated(kim_Es_modes)) deallocate(kim_Es_modes) + if (allocated(kim_Ep_modes)) deallocate(kim_Ep_modes) + if (allocated(kim_Er_modes)) deallocate(kim_Er_modes) + if (allocated(kim_Et_modes)) deallocate(kim_Et_modes) + if (allocated(kim_Ez_modes)) deallocate(kim_Ez_modes) + if (allocated(kim_Br_modes)) deallocate(kim_Br_modes) + + allocate(kim_Es_modes(dim_r, dim_mn)) + allocate(kim_Ep_modes(dim_r, dim_mn)) + allocate(kim_Er_modes(dim_r, dim_mn)) + allocate(kim_Et_modes(dim_r, dim_mn)) + allocate(kim_Ez_modes(dim_r, dim_mn)) + allocate(kim_Br_modes(dim_r, dim_mn)) + + kim_Es_modes = (0.0d0, 0.0d0) + kim_Ep_modes = (0.0d0, 0.0d0) + kim_Er_modes = (0.0d0, 0.0d0) + kim_Et_modes = (0.0d0, 0.0d0) + kim_Ez_modes = (0.0d0, 0.0d0) + kim_Br_modes = (0.0d0, 0.0d0) + + ! ------------------------------------------------------- + ! 2. Loop over modes: solve and store + ! ------------------------------------------------------- + do i_mn = 1, dim_mn + + ! Set mode numbers for this solve + kim_m_mode = m_vals(i_mn) + kim_n_mode = n_vals(i_mn) + + ! Clean up EBdat from previous solve + call deallocate_EBdat() + + ! Create a fresh solver, init grids/equilibrium, solve + call from_kim_factory_get_kim('electrostatic', kim_instance) + call kim_instance%init() + call kim_instance%run() + + ! run() calls postprocess_electric_field which computes + ! Es, Ep, Er, Etheta, Ez from the Poisson-solved Phi. + ! EBdat%Br is set by set_Br_field inside solve_poisson. + + ! Interpolate KIM fields (on xl_grid%xb) onto + ! the QL-Balance radial grid and store per-mode. + kim_npts = size(EBdat%r_grid) + allocate(kim_r(kim_npts)) + kim_r = EBdat%r_grid + + ! Es (perpendicular E field in rsp coordinates) + call interp_complex_profile(kim_npts, kim_r, EBdat%Es, & + dim_r, bal_r, kim_Es_modes(:, i_mn)) + + ! Ep (parallel E field in rsp coordinates) + call interp_complex_profile(kim_npts, kim_r, EBdat%Ep, & + dim_r, bal_r, kim_Ep_modes(:, i_mn)) + + ! Er (radial E field, cylindrical) + call interp_complex_profile(kim_npts, kim_r, EBdat%Er, & + dim_r, bal_r, kim_Er_modes(:, i_mn)) + + ! Etheta -> Et (poloidal E field, cylindrical) + call interp_complex_profile(kim_npts, kim_r, EBdat%Etheta, & + dim_r, bal_r, kim_Et_modes(:, i_mn)) + + ! Ez (axial E field, cylindrical) + call interp_complex_profile(kim_npts, kim_r, EBdat%Ez, & + dim_r, bal_r, kim_Ez_modes(:, i_mn)) + + ! Br (radial magnetic field perturbation) + call interp_complex_profile(kim_npts, kim_r, EBdat%Br, & + dim_r, bal_r, kim_Br_modes(:, i_mn)) + + deallocate(kim_r) + + write(*, '(A,I3,A,I4,A,I4,A)') & + " KIM adapter: solved mode ", i_mn, & + " (m=", m_vals(i_mn), ", n=", n_vals(i_mn), ")" + + end do + + write(*, *) "KIM adapter: all modes solved" + + end subroutine kim_run_for_all_modes subroutine kim_update_profiles() - !! Update KIM internal profiles from QL-Balance evolved - !! parameters. - ! TODO: implement in Task 5 - end subroutine + !! Transfer QL-Balance time-evolved profiles into the + !! wave_code_data in-memory arrays, mirroring what + !! update_background_files does for the KiLCA path. + !! + !! QL-Balance params_b layout (at boundary grid points): + !! params_b(1,:) = density [1/cm^3] + !! params_b(2,:) = Vphi [rad/s] + !! params_b(3,:) = Te [erg] + !! params_b(4,:) = Ti [erg] + !! + !! wave_code_data stores Te, Ti in [eV], Vz in [cm/s]. + use wave_code_data, only: dim_r, & + wcd_n => n, wcd_Te => Te, wcd_Ti => Ti, & + wcd_Vz => Vz, wcd_dPhi0 => dPhi0 + use plasma_parameters, only: params_b + use baseparam_mod, only: ev, rtor + use grid_mod, only: Ercov + + implicit none + + integer :: k + + do k = 1, dim_r + wcd_n(k) = params_b(1, k) + wcd_Te(k) = params_b(3, k) / ev ! erg -> eV + wcd_Ti(k) = params_b(4, k) / ev ! erg -> eV + wcd_Vz(k) = params_b(2, k) * rtor ! rad/s -> cm/s + wcd_dPhi0(k) = -Ercov(k) + end do + + end subroutine kim_update_profiles subroutine kim_get_wave_fields(i_mn) !! Extract wave fields from KIM EBdat into wave_code_data @@ -240,4 +374,63 @@ subroutine kim_get_collision_frequencies() ! TODO: implement in Task 6 end subroutine + ! --------------------------------------------------------------- + ! Helper: deallocate EBdat fields between mode solves + ! --------------------------------------------------------------- + subroutine deallocate_EBdat() + !! Deallocate all allocated components of the global EBdat + !! so that the next call to run() can re-allocate them. + use fields_m, only: EBdat_t + + implicit none + + if (allocated(EBdat%r_grid)) deallocate(EBdat%r_grid) + if (allocated(EBdat%Br)) deallocate(EBdat%Br) + if (allocated(EBdat%Apar)) deallocate(EBdat%Apar) + if (allocated(EBdat%E_perp_psi)) deallocate(EBdat%E_perp_psi) + if (allocated(EBdat%E_perp)) deallocate(EBdat%E_perp) + if (allocated(EBdat%E_perp_MA)) deallocate(EBdat%E_perp_MA) + if (allocated(EBdat%Er)) deallocate(EBdat%Er) + if (allocated(EBdat%Etheta)) deallocate(EBdat%Etheta) + if (allocated(EBdat%Ez)) deallocate(EBdat%Ez) + if (allocated(EBdat%Es)) deallocate(EBdat%Es) + if (allocated(EBdat%Ep)) deallocate(EBdat%Ep) + if (allocated(EBdat%Phi)) deallocate(EBdat%Phi) + if (allocated(EBdat%Phi_e)) deallocate(EBdat%Phi_e) + if (allocated(EBdat%Phi_i)) deallocate(EBdat%Phi_i) + if (allocated(EBdat%Phi_aligned)) deallocate(EBdat%Phi_aligned) + if (allocated(EBdat%Phi_MA)) deallocate(EBdat%Phi_MA) + if (allocated(EBdat%Phi_MA_ideal)) deallocate(EBdat%Phi_MA_ideal) + if (allocated(EBdat%Phi_MA_asymptotic)) deallocate(EBdat%Phi_MA_asymptotic) + + end subroutine deallocate_EBdat + + ! --------------------------------------------------------------- + ! Helper: interpolate a complex profile onto a new grid + ! --------------------------------------------------------------- + subroutine interp_complex_profile(n_old, r_old, z_old, & + n_new, r_new, z_new) + !! Interpolate a complex(8) profile from one radial grid + !! to another. Real and imaginary parts are interpolated + !! independently via the existing interp_profile routine. + implicit none + + integer, intent(in) :: n_old, n_new + real(8), intent(in) :: r_old(n_old), r_new(n_new) + complex(8), intent(in) :: z_old(n_old) + complex(8), intent(out) :: z_new(n_new) + + real(8) :: re_old(n_old), im_old(n_old) + real(8) :: re_new(n_new), im_new(n_new) + + re_old = real(z_old) + im_old = aimag(z_old) + + call interp_profile(n_old, r_old, re_old, n_new, r_new, re_new) + call interp_profile(n_old, r_old, im_old, n_new, r_new, im_new) + + z_new = cmplx(re_new, im_new, kind=8) + + end subroutine interp_complex_profile + end module kim_wave_code_adapter_m From 8eb2740ac751de2e47bd83d6f78da9b413c38db8 Mon Sep 17 00:00:00 2001 From: Markus Markl Date: Mon, 2 Mar 2026 20:22:55 +0100 Subject: [PATCH 06/20] feat(QL-Balance): implement KIM field and auxiliary quantity extraction Implement kim_get_wave_fields (copies per-mode stored Es/Br/Er/Ep/Et/Ez into wave_code_data scalars), kim_get_wave_vectors (per-mode kp/ks with new kim_kp_modes/kim_ks_modes storage arrays), and no-op stubs for kim_get_background_magnetic_fields and kim_get_collision_frequencies (already populated during initialization). Co-Authored-By: Claude Opus 4.6 --- QL-Balance/src/base/kim_wave_code_adapter.f90 | 108 +++++++++++++++--- 1 file changed, 90 insertions(+), 18 deletions(-) diff --git a/QL-Balance/src/base/kim_wave_code_adapter.f90 b/QL-Balance/src/base/kim_wave_code_adapter.f90 index af7f9ba9..4e22a0f6 100644 --- a/QL-Balance/src/base/kim_wave_code_adapter.f90 +++ b/QL-Balance/src/base/kim_wave_code_adapter.f90 @@ -43,6 +43,12 @@ module kim_wave_code_adapter_m complex(8), allocatable :: kim_Ez_modes(:,:) complex(8), allocatable :: kim_Br_modes(:,:) + !! Per-mode stored wave vectors (nrad, dim_mn) + !! kp and ks depend on (m,n) via the equilibrium formulas. + !! Filled by kim_run_for_all_modes, read by kim_get_wave_vectors + real(8), allocatable :: kim_kp_modes(:,:) + real(8), allocatable :: kim_ks_modes(:,:) + contains subroutine kim_initialize(nrad, r_grid) @@ -226,8 +232,9 @@ subroutine kim_run_for_all_modes() implicit none - integer :: i_mn, kim_npts + integer :: i_mn, kim_npts, kim_plasma_npts real(8), allocatable :: kim_r(:) + real(8), allocatable :: kim_plasma_r(:), work_real(:) ! ------------------------------------------------------- ! 1. (Re-)allocate per-mode storage @@ -238,6 +245,8 @@ subroutine kim_run_for_all_modes() if (allocated(kim_Et_modes)) deallocate(kim_Et_modes) if (allocated(kim_Ez_modes)) deallocate(kim_Ez_modes) if (allocated(kim_Br_modes)) deallocate(kim_Br_modes) + if (allocated(kim_kp_modes)) deallocate(kim_kp_modes) + if (allocated(kim_ks_modes)) deallocate(kim_ks_modes) allocate(kim_Es_modes(dim_r, dim_mn)) allocate(kim_Ep_modes(dim_r, dim_mn)) @@ -245,6 +254,8 @@ subroutine kim_run_for_all_modes() allocate(kim_Et_modes(dim_r, dim_mn)) allocate(kim_Ez_modes(dim_r, dim_mn)) allocate(kim_Br_modes(dim_r, dim_mn)) + allocate(kim_kp_modes(dim_r, dim_mn)) + allocate(kim_ks_modes(dim_r, dim_mn)) kim_Es_modes = (0.0d0, 0.0d0) kim_Ep_modes = (0.0d0, 0.0d0) @@ -252,6 +263,8 @@ subroutine kim_run_for_all_modes() kim_Et_modes = (0.0d0, 0.0d0) kim_Ez_modes = (0.0d0, 0.0d0) kim_Br_modes = (0.0d0, 0.0d0) + kim_kp_modes = 0.0d0 + kim_ks_modes = 0.0d0 ! ------------------------------------------------------- ! 2. Loop over modes: solve and store @@ -306,6 +319,22 @@ subroutine kim_run_for_all_modes() deallocate(kim_r) + ! kp and ks (mode-dependent wave vectors on plasma grid) + kim_plasma_npts = kim_plasma%grid_size + allocate(kim_plasma_r(kim_plasma_npts)) + allocate(work_real(kim_plasma_npts)) + kim_plasma_r = kim_plasma%r_grid + + work_real = kim_plasma%kp + call interp_profile(kim_plasma_npts, kim_plasma_r, & + work_real, dim_r, bal_r, kim_kp_modes(:, i_mn)) + + work_real = kim_plasma%ks + call interp_profile(kim_plasma_npts, kim_plasma_r, & + work_real, dim_r, bal_r, kim_ks_modes(:, i_mn)) + + deallocate(kim_plasma_r, work_real) + write(*, '(A,I3,A,I4,A,I4,A)') & " KIM adapter: solved mode ", i_mn, & " (m=", m_vals(i_mn), ", n=", n_vals(i_mn), ")" @@ -350,29 +379,72 @@ subroutine kim_update_profiles() end subroutine kim_update_profiles subroutine kim_get_wave_fields(i_mn) - !! Extract wave fields from KIM EBdat into wave_code_data - !! arrays. + !! Copy per-mode stored fields from kim_*_modes arrays + !! into the wave_code_data module scalars for mode i_mn. + !! Called inside the mode loop in get_dql. + use wave_code_data, only: Es, Br, Er, Ep, Et, Ez, & + Bs, Bp, Bt, Bz + + implicit none + + integer, intent(in) :: i_mn + + Es = kim_Es_modes(:, i_mn) + Br = kim_Br_modes(:, i_mn) + Er = kim_Er_modes(:, i_mn) + Ep = kim_Ep_modes(:, i_mn) + Et = kim_Et_modes(:, i_mn) + Ez = kim_Ez_modes(:, i_mn) + + ! KIM does not compute magnetic field perturbation components + ! other than Br. Set remaining B components to zero. + Bs = (0.0d0, 0.0d0) + Bp = (0.0d0, 0.0d0) + Bt = (0.0d0, 0.0d0) + Bz = (0.0d0, 0.0d0) + + end subroutine kim_get_wave_fields + + subroutine kim_get_wave_vectors(i_mn) + !! Copy per-mode stored wave vectors kp and ks from + !! kim_kp_modes / kim_ks_modes into wave_code_data for + !! mode i_mn. kp and ks depend on (m,n) through + !! the equilibrium geometry and are recomputed by KIM + !! for each mode in kim_run_for_all_modes. + use wave_code_data, only: kp, ks + + implicit none + integer, intent(in) :: i_mn - ! TODO: implement in Task 6 - end subroutine - subroutine kim_get_wave_vectors() - !! Extract wave vectors from KIM plasma into wave_code_data - !! arrays. - ! TODO: implement in Task 6 - end subroutine + kp = kim_kp_modes(:, i_mn) + ks = kim_ks_modes(:, i_mn) + + end subroutine kim_get_wave_vectors subroutine kim_get_background_magnetic_fields() - !! Extract background B from KIM equilibrium into - !! wave_code_data arrays. - ! TODO: implement in Task 6 - end subroutine + !! No-op: B0, B0t, B0z were already populated during + !! kim_initialize and stored directly in wave_code_data. + !! This subroutine exists to satisfy the adapter interface + !! contract expected by get_dql. + + implicit none + + ! Nothing to do -- B0, B0t, B0z are already set. + + end subroutine kim_get_background_magnetic_fields subroutine kim_get_collision_frequencies() - !! Extract collision frequencies from KIM plasma into - !! wave_code_data arrays. - ! TODO: implement in Task 6 - end subroutine + !! No-op: nue and nui were already populated during + !! kim_initialize and stored directly in wave_code_data. + !! This subroutine exists to satisfy the adapter interface + !! contract expected by get_dql. + + implicit none + + ! Nothing to do -- nue, nui are already set. + + end subroutine kim_get_collision_frequencies ! --------------------------------------------------------------- ! Helper: deallocate EBdat fields between mode solves From 3977fe79ba14cfe8beb143665c1304fb6712cfbb Mon Sep 17 00:00:00 2001 From: Markus Markl Date: Mon, 2 Mar 2026 20:27:13 +0100 Subject: [PATCH 07/20] feat(QL-Balance): wire KIM adapter into get_dql and initialization Adds select case(wave_code) dispatch in: - initialize_wave_code_interface: KiLCA vs KIM init - get_dql timestep calls: profile update + solve - Mode loop: field/vector retrieval per mode - Misalignment diffusion: wave field refresh - Vacuum Br for form factor - FLRE Br for form factor - Current density approximation (zeroed for KIM) Also makes kim_Br_modes public in the adapter module so get_dql can access per-mode Br arrays directly. Co-Authored-By: Claude Opus 4.6 --- QL-Balance/src/base/get_dql.f90 | 77 ++++++++++++++----- QL-Balance/src/base/kim_wave_code_adapter.f90 | 2 +- QL-Balance/src/base/wave_code_data_64bit.f90 | 70 ++++++++++------- 3 files changed, 101 insertions(+), 48 deletions(-) diff --git a/QL-Balance/src/base/get_dql.f90 b/QL-Balance/src/base/get_dql.f90 index 03333e53..dcab2dd2 100644 --- a/QL-Balance/src/base/get_dql.f90 +++ b/QL-Balance/src/base/get_dql.f90 @@ -14,10 +14,12 @@ subroutine get_dql , r_resonant, d11_misalign, Es_pert_flux, Ipar use plasma_parameters use baseparam_mod, only: Z_i, e_charge, am, p_mass, c, e_mass, ev, rtor, pi, rsepar - use control_mod, only: irf, suppression_mode, misalign_diffusion, type_of_run + use control_mod, only: irf, suppression_mode, misalign_diffusion, type_of_run, wave_code use time_evolution, only: save_prof_time_step, time_ind, br_formfactor, br_vac_res use h5mod use wave_code_data + use kim_wave_code_adapter_m, only: kim_update_profiles, kim_run_for_all_modes, & + kim_get_wave_fields, kim_get_wave_vectors, kim_vac_Br, kim_Br_modes use QLBalance_diag, only: i_mn_loop use QLBalance_kinds, only: dp use PolyLagrangeInterpolation @@ -116,10 +118,17 @@ subroutine get_dql ! Compute diffusion coefficient matrices: - if (irf .eq. 1) call update_background_files(path2profs) - if (irf .eq. 1) call get_wave_code_data(1, dim_mn) - if (irf .eq. 1) call get_background_magnetic_fields_from_wave_code(flre_cd_ptr(1), dim_r, r, B0t, B0z, B0) - if (irf .eq. 1) call get_collision_frequences_from_wave_code(flre_cd_ptr(1), dim_r, r, nui, nue) + select case (trim(wave_code)) + case ('KiLCA') + if (irf .eq. 1) call update_background_files(path2profs) + if (irf .eq. 1) call get_wave_code_data(1, dim_mn) + if (irf .eq. 1) call get_background_magnetic_fields_from_wave_code(flre_cd_ptr(1), dim_r, r, B0t, B0z, B0) + if (irf .eq. 1) call get_collision_frequences_from_wave_code(flre_cd_ptr(1), dim_r, r, nui, nue) + case ('KIM') + if (irf .eq. 1) call kim_update_profiles() + if (irf .eq. 1) call kim_run_for_all_modes() + ! B0 and collision freqs already set in kim_initialize + end select ! nu_e=15.4d-6*params_b(1,:)/sqrt(params_b(3,:)/ev)**3 & ! *(23.d0-0.5d0*log(params_b(1,:)/(params_b(3,:)/ev)**3)) @@ -142,10 +151,16 @@ subroutine get_dql !sum over modes: do i_mn = 1, dim_mn - call get_wave_vectors_from_wave_code(flre_cd_ptr(i_mn), dim_r, r, & - m_vals(i_mn), n_vals(i_mn), ks, kp) - call get_wave_fields_from_wave_code(flre_cd_ptr(i_mn), dim_r, r, & - m_vals(i_mn), n_vals(i_mn), Er, Es, Ep, Et, Ez, Br, Bs, Bp, Bt, Bz) + select case (trim(wave_code)) + case ('KiLCA') + call get_wave_vectors_from_wave_code(flre_cd_ptr(i_mn), dim_r, r, & + m_vals(i_mn), n_vals(i_mn), ks, kp) + call get_wave_fields_from_wave_code(flre_cd_ptr(i_mn), dim_r, r, & + m_vals(i_mn), n_vals(i_mn), Er, Es, Ep, Et, Ez, Br, Bs, Bp, Bt, Bz) + case ('KIM') + call kim_get_wave_vectors(i_mn) + call kim_get_wave_fields(i_mn) + end select om_E = ks * c * dPhi0 / B0 vT_e = sqrt(params_b(3, :)/e_mass) vT_i = sqrt(params_b(4, :)/p_mass/am) @@ -169,8 +184,13 @@ subroutine get_dql end if if (misalign_diffusion .eqv. .true.) then - call get_wave_fields_from_wave_code(flre_cd_ptr(i_mn), dim_r, r, & - m_vals(i_mn), n_vals(i_mn), Er, Es, Ep, Et, Ez, Br, Bs, Bp, Bt, Bz) + select case (trim(wave_code)) + case ('KiLCA') + call get_wave_fields_from_wave_code(flre_cd_ptr(i_mn), dim_r, r, & + m_vals(i_mn), n_vals(i_mn), Er, Es, Ep, Et, Ez, Br, Bs, Bp, Bt, Bz) + case ('KIM') + call kim_get_wave_fields(i_mn) ! already loaded, but refresh + end select ! caluclate part of perpendicular electric field perturbation that comes from if (.not. allocated(coef)) allocate(coef(0:nder,nlagr)) @@ -197,9 +217,13 @@ subroutine get_dql Es_pert_flux = Es_pert_flux + Es_pert_flux_temp end if - call get_wave_fields_from_wave_code(vac_cd_ptr(i_mn), dim_r, r, & - m_vals(i_mn), n_vals(i_mn), Bz, Bz, Bz, Bz, Bz, Br, Bz, Bz, Bz, Bz) - + select case (trim(wave_code)) + case ('KiLCA') + call get_wave_fields_from_wave_code(vac_cd_ptr(i_mn), dim_r, r, & + m_vals(i_mn), n_vals(i_mn), Bz, Bz, Bz, Bz, Bz, Br, Bz, Bz, Bz, Bz) + case ('KIM') + Br = kim_vac_Br + end select formfactor = (1.d0, 0.d0)/Br ! In case that the tmhd code uses double sided Fourier series, it @@ -219,8 +243,13 @@ subroutine get_dql end if end if - call get_wave_fields_from_wave_code(flre_cd_ptr(i_mn), dim_r, r, & - m_vals(i_mn), n_vals(i_mn), Bz, Bz, Bz, Bz, Bz, Br, Bz, Bz, Bz, Bz) + select case (trim(wave_code)) + case ('KiLCA') + call get_wave_fields_from_wave_code(flre_cd_ptr(i_mn), dim_r, r, & + m_vals(i_mn), n_vals(i_mn), Bz, Bz, Bz, Bz, Bz, Br, Bz, Bz, Bz, Bz) + case ('KIM') + Br = kim_Br_modes(:, i_mn) + end select formfactor = Br * formfactor ! todo: interpolate formfactor at resonant surface and write out. This is Brtot/Brvac at the resonant surface @@ -239,8 +268,20 @@ subroutine get_dql dqli21_loc = dqli21_loc + di21*spec_weight dqli22_loc = dqli22_loc + di22*spec_weight - call get_current_densities_from_wave_code(flre_cd_ptr(i_mn), dim_r, r, & - m_vals(i_mn), n_vals(i_mn), Jri, Jsi, Jpi, Jre, Jse, Jpe) + select case (trim(wave_code)) + case ('KiLCA') + call get_current_densities_from_wave_code(flre_cd_ptr(i_mn), dim_r, r, & + m_vals(i_mn), n_vals(i_mn), Jri, Jsi, Jpi, Jre, Jse, Jpe) + case ('KIM') + ! KIM only provides total jpar; approximate Jpe with it, zero rest + ! TODO: improve when KIM provides per-species currents + Jpe = (0.0d0, 0.0d0) + Jpi = (0.0d0, 0.0d0) + Jre = (0.0d0, 0.0d0) + Jri = (0.0d0, 0.0d0) + Jse = (0.0d0, 0.0d0) + Jsi = (0.0d0, 0.0d0) + end select call integrate_parallel_current(dim_r, r, Jpe, Jpi, Ipar) end do diff --git a/QL-Balance/src/base/kim_wave_code_adapter.f90 b/QL-Balance/src/base/kim_wave_code_adapter.f90 index 4e22a0f6..51f25d91 100644 --- a/QL-Balance/src/base/kim_wave_code_adapter.f90 +++ b/QL-Balance/src/base/kim_wave_code_adapter.f90 @@ -41,7 +41,7 @@ module kim_wave_code_adapter_m complex(8), allocatable :: kim_Er_modes(:,:) complex(8), allocatable :: kim_Et_modes(:,:) complex(8), allocatable :: kim_Ez_modes(:,:) - complex(8), allocatable :: kim_Br_modes(:,:) + complex(8), allocatable, public :: kim_Br_modes(:,:) !! Per-mode stored wave vectors (nrad, dim_mn) !! kp and ks depend on (m,n) via the equilibrium formulas. diff --git a/QL-Balance/src/base/wave_code_data_64bit.f90 b/QL-Balance/src/base/wave_code_data_64bit.f90 index 1d6326db..4d878b5c 100644 --- a/QL-Balance/src/base/wave_code_data_64bit.f90 +++ b/QL-Balance/src/base/wave_code_data_64bit.f90 @@ -33,7 +33,8 @@ subroutine initialize_wave_code_interface(nrad, r_grid) use wave_code_data use h5mod - use control_mod, only: readfromtimestep, debug_mode, ihdf5IO + use control_mod, only: readfromtimestep, debug_mode, ihdf5IO, wave_code + use kim_wave_code_adapter_m, only: kim_initialize implicit none @@ -44,40 +45,51 @@ subroutine initialize_wave_code_interface(nrad, r_grid) if (debug_mode) write(*,*) "Debug: Coming into initialize_wave_code_interface" - call read_antenna_modes(flre_path) ! read antenna modes from modes.in file and allocate arrays for mode's numbers - - call allocate_wave_code_data(nrad, r_grid) -!call read_background_profiles(path2profs); - if (ihdf5IO .eq. 1) then - if (readfromtimestep .eq. 0) then - if(debug_mode) write(*,*) "reading initial profiles" - call read_background_profiles_h5 + select case (trim(wave_code)) + case ('KiLCA') + call read_antenna_modes(flre_path) ! read antenna modes from modes.in file and allocate arrays for mode's numbers + + call allocate_wave_code_data(nrad, r_grid) + !call read_background_profiles(path2profs); + if (ihdf5IO .eq. 1) then + if (readfromtimestep .eq. 0) then + if(debug_mode) write(*,*) "reading initial profiles" + call read_background_profiles_h5 + else + if(debug_mode) write(*,*) "reading time evolved profiles" + call read_background_profiles_h5_timeevol(readfromtimestep) + end if else - if(debug_mode) write(*,*) "reading time evolved profiles" - call read_background_profiles_h5_timeevol(readfromtimestep) + call read_background_profiles(profile_path) end if - else - call read_background_profiles(profile_path) - end if - call interp_background_profiles() ! interpolate the initial profiles to the balance grid + call interp_background_profiles() ! interpolate the initial profiles to the balance grid - if (debug_mode) write(*,*) "Debug: after interp_background_profiles" - ! vacuum fields for the whole spectrum: - do k = 1, dim_mn - call clear_wave_code_data(vac_cd_ptr(k)); - call calc_wave_code_data_for_mode(vac_cd_ptr(k), trim(vac_path), len(trim(vac_path)), m_vals(k), n_vals(k)); - end do + if (debug_mode) write(*,*) "Debug: after interp_background_profiles" + ! vacuum fields for the whole spectrum: + do k = 1, dim_mn + call clear_wave_code_data(vac_cd_ptr(k)); + call calc_wave_code_data_for_mode(vac_cd_ptr(k), trim(vac_path), len(trim(vac_path)), m_vals(k), n_vals(k)); + end do + + if (debug_mode) write(*,*) "Debug: before get_background_magnetic_fields_from_wave_code" + call get_background_magnetic_fields_from_wave_code(vac_cd_ptr(1), dim_r, r, B0t, B0z, B0); + call get_collision_frequences_from_wave_code(vac_cd_ptr(1), dim_r, r, nui, nue); + vac_call_ind = vac_call_ind + 1; + + ! Enable interface mode: KiLCA gets profiles directly from QL-Balance for all run types + ! that use initialize_wave_code_interface. This is set AFTER vacuum field initialization + ! to preserve original file-based behavior there. + call set_flag_for_profiles_in_background_input_file(trim(flre_path), -1) + + case ('KIM') + call kim_initialize(nrad, r_grid) - if (debug_mode) write(*,*) "Debug: before get_background_magnetic_fields_from_wave_code" - call get_background_magnetic_fields_from_wave_code(vac_cd_ptr(1), dim_r, r, B0t, B0z, B0); - call get_collision_frequences_from_wave_code(vac_cd_ptr(1), dim_r, r, nui, nue); - vac_call_ind = vac_call_ind + 1; + case default + write(*,*) "ERROR: Unknown wave_code: ", trim(wave_code) + stop 1 - ! Enable interface mode: KiLCA gets profiles directly from QL-Balance for all run types - ! that use initialize_wave_code_interface. This is set AFTER vacuum field initialization - ! to preserve original file-based behavior there. - call set_flag_for_profiles_in_background_input_file(trim(flre_path), -1) + end select if (debug_mode) write(*,*) "Debug: Going out of initialize_wave_code_interface" end subroutine From b95a321c487d2be6ed0608773a92503e06afc756 Mon Sep 17 00:00:00 2001 From: Markus Markl Date: Mon, 2 Mar 2026 20:34:06 +0100 Subject: [PATCH 08/20] refactor(KIM): rename resonances_mod to kim_resonances_m Avoids duplicate symbol collision with QL-Balance's resonances_mod when both libraries are linked into ql-balance.x. Also moves cerf and ddeabm link dependencies from KIM_exe to KIM_lib so they are available when KIM_lib is consumed by other targets. Co-Authored-By: Claude Opus 4.6 --- KIM/src/CMakeLists.txt | 4 ++-- KIM/src/background_equilibrium/species_mod.f90 | 2 +- KIM/src/electrostatic_poisson/solve_poisson.f90 | 2 +- KIM/src/grid/grid_mod.f90 | 6 +++--- KIM/src/grid/prepare_resonances.f90 | 2 +- KIM/src/grid/recnsplit.f90 | 2 +- KIM/src/grid/resonances_mod.f90 | 4 ++-- KIM/src/kernels/kernel.f90 | 4 ++-- 8 files changed, 13 insertions(+), 13 deletions(-) diff --git a/KIM/src/CMakeLists.txt b/KIM/src/CMakeLists.txt index 4a95d54a..be601f6e 100644 --- a/KIM/src/CMakeLists.txt +++ b/KIM/src/CMakeLists.txt @@ -148,6 +148,8 @@ target_link_libraries(KIM_lib PUBLIC hdf5::hdf5_hl_fortran kamel_hdf5_tools kamel_equil + cerf + ddeabm ) add_executable(KIM_exe kim_main.f90) @@ -157,7 +159,5 @@ set_target_properties(KIM_exe PROPERTIES target_link_libraries(KIM_exe KIM_lib lapack - cerf slatec - ddeabm ) diff --git a/KIM/src/background_equilibrium/species_mod.f90 b/KIM/src/background_equilibrium/species_mod.f90 index 528212d2..16aafc37 100644 --- a/KIM/src/background_equilibrium/species_mod.f90 +++ b/KIM/src/background_equilibrium/species_mod.f90 @@ -882,7 +882,7 @@ subroutine plot_species(spec) subroutine calculate_susc_funcs_profiles(spec, mphi) - use resonances_mod, only: r_res + use kim_resonances_m, only: r_res use grid_m, only: width_res implicit none diff --git a/KIM/src/electrostatic_poisson/solve_poisson.f90 b/KIM/src/electrostatic_poisson/solve_poisson.f90 index 74fae847..8296900f 100644 --- a/KIM/src/electrostatic_poisson/solve_poisson.f90 +++ b/KIM/src/electrostatic_poisson/solve_poisson.f90 @@ -97,7 +97,7 @@ subroutine write_A_matrix_sparse_check_to_file subroutine create_rhs_vector(type, K_rho_B, rhs_vec) - use resonances_mod, only: index_rg_res, r_res + use kim_resonances_m, only: index_rg_res, r_res use functions_m, only: varphi_l use grid_m, only: xl_grid use IO_collection_m, only: write_complex_profile, plot_profile diff --git a/KIM/src/grid/grid_mod.f90 b/KIM/src/grid/grid_mod.f90 index 0dffd35e..7eb54e4d 100644 --- a/KIM/src/grid/grid_mod.f90 +++ b/KIM/src/grid/grid_mod.f90 @@ -65,7 +65,7 @@ module grid_m contains subroutine ensure_node_at_r_res(this) - use resonances_mod, only: r_res + use kim_resonances_m, only: r_res use KIM_kinds_m, only: dp implicit none class(grid_type), intent(inout) :: this @@ -191,7 +191,7 @@ subroutine grid_init(this, npts, min_val, max_val, name) subroutine grid_generate(this) - use resonances_mod, only: r_res, index_rg_res + use kim_resonances_m, only: r_res, index_rg_res use config_m, only: fdebug, output_path implicit none @@ -261,7 +261,7 @@ end subroutine grid_generate subroutine grid_generate_equidistant(this) - use resonances_mod, only: r_res, index_rg_res + use kim_resonances_m, only: r_res, index_rg_res use config_m, only: fdebug, output_path implicit none diff --git a/KIM/src/grid/prepare_resonances.f90 b/KIM/src/grid/prepare_resonances.f90 index f0df47f6..6cf77774 100644 --- a/KIM/src/grid/prepare_resonances.f90 +++ b/KIM/src/grid/prepare_resonances.f90 @@ -1,6 +1,6 @@ subroutine prepare_resonances - use resonances_mod + use kim_resonances_m use config_m, only: hdf5_output use setup_m, only: m_mode, n_mode, type_br_field use species_m, only: plasma diff --git a/KIM/src/grid/recnsplit.f90 b/KIM/src/grid/recnsplit.f90 index f600c98e..8f27e270 100644 --- a/KIM/src/grid/recnsplit.f90 +++ b/KIM/src/grid/recnsplit.f90 @@ -1,6 +1,6 @@ subroutine recnsplit(r,recnsp) - use resonances_mod, only: r_res, prop + use kim_resonances_m, only: r_res, prop use grid_m, only: width_res, ampl_res use KIM_kinds_m, only: dp diff --git a/KIM/src/grid/resonances_mod.f90 b/KIM/src/grid/resonances_mod.f90 index a65fac1f..31eb85fa 100644 --- a/KIM/src/grid/resonances_mod.f90 +++ b/KIM/src/grid/resonances_mod.f90 @@ -1,4 +1,4 @@ -module resonances_mod +module kim_resonances_m use KIM_kinds_m, only: dp @@ -9,4 +9,4 @@ module resonances_mod integer :: index_rg_res logical :: prop = .true. -end module resonances_mod +end module kim_resonances_m diff --git a/KIM/src/kernels/kernel.f90 b/KIM/src/kernels/kernel.f90 index b2f1b174..b7417eb7 100644 --- a/KIM/src/kernels/kernel.f90 +++ b/KIM/src/kernels/kernel.f90 @@ -736,7 +736,7 @@ subroutine FP_calc_kernel_element_ions(l, lp, k_rho_phi, k_rho_B, k_j_phi, k_j_B use constants_m, only: com_unit, sol use config_m, only: turn_off_ions, turn_off_electrons, artificial_debye_case, ion_flr_scale_factor use grid_m, only: xl_grid - use resonances_mod, only: r_res + use kim_resonances_m, only: r_res use setup_m, only: mphi_max implicit none @@ -1889,7 +1889,7 @@ end subroutine fill_kernels_krook_fp ! use constants_m, only: com_unit, sol ! use config_m, only: turn_off_ions, turn_off_electrons, artificial_debye_case ! use grid_m, only: xl_grid - ! use resonances_mod, only: r_res + ! use kim_resonances_m, only: r_res ! implicit none From 84b4ebf0cf7a507522def7a66a9418f722c97b75 Mon Sep 17 00:00:00 2001 From: Markus Markl Date: Mon, 2 Mar 2026 20:37:53 +0100 Subject: [PATCH 09/20] feat(QL-Balance): add kim_config_path to balance namelist Allow users to configure the KIM namelist file path via balance_conf.nml instead of relying on the hardcoded default "./KIM_config.nml". Co-Authored-By: Claude Opus 4.6 --- QL-Balance/namelists/balance_conf.nml | 1 + QL-Balance/src/base/control_mod.f90 | 1 + QL-Balance/src/base/kim_wave_code_adapter.f90 | 11 ++++++----- QL-Balance/src/base/read_config.f90 | 6 ++++-- 4 files changed, 12 insertions(+), 7 deletions(-) diff --git a/QL-Balance/namelists/balance_conf.nml b/QL-Balance/namelists/balance_conf.nml index 2d148a6b..ce901f0c 100644 --- a/QL-Balance/namelists/balance_conf.nml +++ b/QL-Balance/namelists/balance_conf.nml @@ -56,5 +56,6 @@ equil_path = '/equil_r_q_psi.dat' , type_of_run = 'SingleStep' wave_code = 'KiLCA' + kim_config_path = './KIM_config.nml' ihdf5IO = 1 / diff --git a/QL-Balance/src/base/control_mod.f90 b/QL-Balance/src/base/control_mod.f90 index 67afbe1c..3298ff11 100644 --- a/QL-Balance/src/base/control_mod.f90 +++ b/QL-Balance/src/base/control_mod.f90 @@ -6,6 +6,7 @@ module control_mod character(100) :: type_of_run character(len=32) :: wave_code = 'KiLCA' ! 'KiLCA' or 'KIM' + character(len=1024) :: kim_config_path = './KIM_config.nml' integer :: ihdf5IO ! added: Markus Markl logical :: paramscan ! added: Markus Markl, 03.03.2021 logical :: timing_mode ! added by Markus Markl 06.04.2021 diff --git a/QL-Balance/src/base/kim_wave_code_adapter.f90 b/QL-Balance/src/base/kim_wave_code_adapter.f90 index 51f25d91..08c8bddb 100644 --- a/QL-Balance/src/base/kim_wave_code_adapter.f90 +++ b/QL-Balance/src/base/kim_wave_code_adapter.f90 @@ -10,7 +10,7 @@ module kim_wave_code_adapter_m use species_m, only: kim_plasma => plasma use equilibrium_m, only: kim_B0 => B0 use config_m, only: kim_type_of_run => type_of_run, & - nml_config_path + kim_nml_config_path => nml_config_path use setup_m, only: kim_m_mode => m_mode, kim_n_mode => n_mode use grid_m, only: kim_xl_grid => xl_grid, kim_rg_grid => rg_grid use fields_m, only: EBdat @@ -56,6 +56,7 @@ subroutine kim_initialize(nrad, r_grid) !! Populate wave_code_data module arrays with background !! quantities by running the electrostatic solver for the !! first mode to generate equilibrium and plasma backgrounds. + use control_mod, only: kim_config_path use wave_code_data, only: dim_mn, m_vals, n_vals, & r => r, q => q, n => n, Te => Te, Ti => Ti, & Vth => Vth, Vz => Vz, dPhi0 => dPhi0, & @@ -75,10 +76,10 @@ subroutine kim_initialize(nrad, r_grid) ! ----------------------------------------------------------- ! 1. Set KIM config path and initialize KIM backend ! ----------------------------------------------------------- - ! nml_config_path (from config_m) tells KIM where its - ! namelist file lives. Default is ./KIM_config.nml which - ! the user can override via Task 8 (balance namelist). - nml_config_path = "./KIM_config.nml" + ! kim_nml_config_path (from config_m) tells KIM where its + ! namelist file lives. Controlled by kim_config_path in + ! balance_conf.nml (default: ./KIM_config.nml). + kim_nml_config_path = trim(kim_config_path) call kim_init() diff --git a/QL-Balance/src/base/read_config.f90 b/QL-Balance/src/base/read_config.f90 index 64e3f49a..c003b6d7 100644 --- a/QL-Balance/src/base/read_config.f90 +++ b/QL-Balance/src/base/read_config.f90 @@ -2,7 +2,8 @@ subroutine read_config use baseparam_mod, only: btor, rtor, rsepar, dperp, Z_i, am, urelax use control_mod, only: eps, paramscan, diagnostics_output, suppression_mode, debug_mode, & readfromtimestep, temperature_limit, gyro_current_study, & - misalign_diffusion, equil_path, ihdf5IO, wave_code + misalign_diffusion, equil_path, ihdf5IO, wave_code, & + kim_config_path use grid_mod, only: rmin, rmax, npoimin, gg_factor, gg_width, gg_r_res, iboutype, rb_cut_in, & re_cut_in, rb_cut_out, re_cut_out use h5mod, only: path2inp, path2out, path2time @@ -24,7 +25,7 @@ subroutine read_config suppression_mode, debug_mode, readfromtimestep, path2time, ramp_up_mode, t_max_ramp_up, & temperature_limit, antenna_max_stopping, gyro_current_study, viscosity_factor, & misalign_diffusion, equil_path, ihdf5IO, type_of_run, wave_code, & - set_constant_time_step, constant_time_step, urelax + set_constant_time_step, constant_time_step, urelax, kim_config_path ! read the parameters from namelist file open (newunit=u, file=config_file, status="old", action="read", iostat=ios) @@ -80,6 +81,7 @@ subroutine read_config write (*, "(A,L0)") " set_constant_time_step = ", set_constant_time_step write (*, "(A,ES15.8,A)") " constant_time_step = ", constant_time_step, " s" write (*, "(A,ES15.8)") " urelax = ", urelax + write (*, "(A,A)") " kim_config_path = ", trim(adjustl(kim_config_path)) write (*, *) "" write (*, *) "=================================================================================" end subroutine read_config From c6ea1833e1dae3d5fadf713dbabfbac6c7f14dd7 Mon Sep 17 00:00:00 2001 From: Markus Markl Date: Mon, 2 Mar 2026 20:44:07 +0100 Subject: [PATCH 10/20] test(QL-Balance): add unit test for KIM adapter module Tests wave_code default value, interp_complex_profile with constant and linear functions. Verifies adapter module compilation and linkage against KIM_lib. Full integration tests require KIM config data and are deferred to end-to-end validation. Co-Authored-By: Claude Opus 4.6 --- QL-Balance/src/base/kim_wave_code_adapter.f90 | 1 + QL-Balance/src/test/CMakeLists.txt | 5 + QL-Balance/src/test/test_kim_adapter.f90 | 165 ++++++++++++++++++ 3 files changed, 171 insertions(+) create mode 100644 QL-Balance/src/test/test_kim_adapter.f90 diff --git a/QL-Balance/src/base/kim_wave_code_adapter.f90 b/QL-Balance/src/base/kim_wave_code_adapter.f90 index 08c8bddb..9b7fe319 100644 --- a/QL-Balance/src/base/kim_wave_code_adapter.f90 +++ b/QL-Balance/src/base/kim_wave_code_adapter.f90 @@ -27,6 +27,7 @@ module kim_wave_code_adapter_m public :: kim_get_wave_vectors public :: kim_get_background_magnetic_fields public :: kim_get_collision_frequencies + public :: interp_complex_profile ! exposed for testing !! Module-level KIM solver instance (reused across calls) class(kim_t), allocatable :: kim_instance diff --git a/QL-Balance/src/test/CMakeLists.txt b/QL-Balance/src/test/CMakeLists.txt index a4c03c76..b68963ee 100644 --- a/QL-Balance/src/test/CMakeLists.txt +++ b/QL-Balance/src/test/CMakeLists.txt @@ -21,3 +21,8 @@ add_fortran_test(test_sparse) add_executable(test_rhs_balance.x test_rhs_balance.f90) target_link_libraries(test_rhs_balance.x PUBLIC ql-balance_lib) add_fortran_test(test_rhs_balance) + +# KIM adapter module tests (unit tests for adapter helpers and linkage) +add_executable(test_kim_adapter.x test_kim_adapter.f90) +target_link_libraries(test_kim_adapter.x PUBLIC ql-balance_lib) +add_fortran_test(test_kim_adapter) diff --git a/QL-Balance/src/test/test_kim_adapter.f90 b/QL-Balance/src/test/test_kim_adapter.f90 new file mode 100644 index 00000000..fd445f4e --- /dev/null +++ b/QL-Balance/src/test/test_kim_adapter.f90 @@ -0,0 +1,165 @@ +program test_kim_adapter + ! + ! Unit tests for the KIM wave code adapter module. + ! + ! These tests verify: + ! 1. Module compilation and linkage (the adapter + KIM_lib resolve) + ! 2. interp_complex_profile helper correctness + ! 3. wave_code selector variable default + ! + ! Note: Full integration tests (calling kim_initialize) require + ! KIM config files and profile data, and are covered by Task 10 + ! (end-to-end SingleStep validation). + ! + use kim_wave_code_adapter_m, only: interp_complex_profile + use control_mod, only: wave_code + + implicit none + + integer :: num_passed, num_failed + + num_passed = 0 + num_failed = 0 + + print *, "========================================" + print *, " KIM Adapter Module Tests" + print *, "========================================" + print *, "" + + call test_wave_code_default() + call test_interp_complex_constant() + call test_interp_complex_linear() + + print *, "" + print *, "========================================" + print *, " Test Summary" + print *, "========================================" + print '(A,I3,A)', " Passed: ", num_passed, " tests" + print '(A,I3,A)', " Failed: ", num_failed, " tests" + print *, "========================================" + + if (num_failed > 0) then + print *, "TESTS FAILED" + stop 1 + else + print *, "ALL TESTS PASSED" + end if + +contains + + subroutine assert_equal_real(actual, expected, tol, label) + real(8), intent(in) :: actual, expected, tol + character(len=*), intent(in) :: label + + if (abs(actual - expected) > tol) then + print '(A,A,A,ES15.8,A,ES15.8)', " FAIL: ", label, & + " got ", actual, " expected ", expected + num_failed = num_failed + 1 + else + print '(A,A)', " PASS: ", label + num_passed = num_passed + 1 + end if + end subroutine + + subroutine assert_equal_complex(actual, expected, tol, label) + complex(8), intent(in) :: actual, expected + real(8), intent(in) :: tol + character(len=*), intent(in) :: label + + if (abs(actual - expected) > tol) then + print '(A,A,A,2ES15.8,A,2ES15.8)', " FAIL: ", label, & + " got (", real(actual), aimag(actual), & + ") expected (", real(expected), aimag(expected), ")" + num_failed = num_failed + 1 + else + print '(A,A)', " PASS: ", label + num_passed = num_passed + 1 + end if + end subroutine + + ! ------------------------------------------------------------------ + ! Test: wave_code defaults to 'KiLCA' + ! ------------------------------------------------------------------ + subroutine test_wave_code_default() + print *, "--- test_wave_code_default ---" + + if (trim(wave_code) == 'KiLCA') then + print '(A)', " PASS: wave_code default is 'KiLCA'" + num_passed = num_passed + 1 + else + print '(A,A,A)', " FAIL: wave_code default is '", & + trim(wave_code), "', expected 'KiLCA'" + num_failed = num_failed + 1 + end if + end subroutine + + ! ------------------------------------------------------------------ + ! Test: interp_complex_profile with constant function + ! A constant complex function should interpolate exactly. + ! ------------------------------------------------------------------ + subroutine test_interp_complex_constant() + integer, parameter :: n_old = 5, n_new = 3 + real(8) :: r_old(n_old), r_new(n_new) + complex(8) :: z_old(n_old), z_new(n_new) + complex(8), parameter :: val = (2.5d0, -1.3d0) + real(8), parameter :: tol = 1.0d-10 + integer :: i + + print *, "--- test_interp_complex_constant ---" + + ! Old grid: [10, 20, 30, 40, 50] + do i = 1, n_old + r_old(i) = 10.0d0 * dble(i) + z_old(i) = val + end do + + ! New grid: [15, 25, 35] (within old grid bounds) + r_new(1) = 15.0d0 + r_new(2) = 25.0d0 + r_new(3) = 35.0d0 + + call interp_complex_profile(n_old, r_old, z_old, n_new, r_new, z_new) + + do i = 1, n_new + call assert_equal_complex(z_new(i), val, tol, & + "constant interp at new grid point") + end do + end subroutine + + ! ------------------------------------------------------------------ + ! Test: interp_complex_profile with linear function + ! z(r) = (1+2i)*r should interpolate exactly with polynomial interp. + ! ------------------------------------------------------------------ + subroutine test_interp_complex_linear() + integer, parameter :: n_old = 10, n_new = 5 + real(8) :: r_old(n_old), r_new(n_new) + complex(8) :: z_old(n_old), z_new(n_new) + complex(8) :: expected + real(8), parameter :: tol = 1.0d-8 + integer :: i + + print *, "--- test_interp_complex_linear ---" + + ! Old grid: uniform [10, 20, 30, ..., 100] + do i = 1, n_old + r_old(i) = 10.0d0 * dble(i) + z_old(i) = cmplx(1.0d0, 2.0d0, kind=8) * r_old(i) + end do + + ! New grid: non-coincident points + r_new(1) = 15.0d0 + r_new(2) = 25.0d0 + r_new(3) = 55.0d0 + r_new(4) = 75.0d0 + r_new(5) = 95.0d0 + + call interp_complex_profile(n_old, r_old, z_old, n_new, r_new, z_new) + + do i = 1, n_new + expected = cmplx(1.0d0, 2.0d0, kind=8) * r_new(i) + call assert_equal_complex(z_new(i), expected, tol, & + "linear interp at new grid point") + end do + end subroutine + +end program test_kim_adapter From 8bf597aabcba6ae9de9a5f240ae32c7601800087 Mon Sep 17 00:00:00 2001 From: Markus Markl Date: Tue, 24 Mar 2026 09:55:48 +0100 Subject: [PATCH 11/20] build(KIM): fix ddeabm duplicate symbols and KiLCA link order MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit ddeabm's i1mach had all DATA statements commented out, shadowing the correct slatec version at link time and breaking DGAMLN/Bessel functions. Explicitly list only the ODE-solver sources instead of globbing. Move KiLCA linkage from KIM_lib (PUBLIC) to KIM_exe only — KIM does not use KiLCA symbols directly, and the duplicate linkage corrupted symbol resolution when QL-Balance linked both KIM_lib and kilca_lib. Co-Authored-By: Claude Opus 4.6 (1M context) --- KIM/src/CMakeLists.txt | 8 ++++++-- KIM/src/math/ddeabm/CMakeLists.txt | 15 ++++++++++++++- 2 files changed, 20 insertions(+), 3 deletions(-) diff --git a/KIM/src/CMakeLists.txt b/KIM/src/CMakeLists.txt index be601f6e..a733e5f2 100644 --- a/KIM/src/CMakeLists.txt +++ b/KIM/src/CMakeLists.txt @@ -133,9 +133,12 @@ endif() set(KILCA_PATH ${CMAKE_SOURCE_DIR}/KiLCA) set(KiLCA_lib_path ${CMAKE_SOURCE_DIR}/build/install/lib/libKiLCA_Lib_V_2.4.2_MDNO_FPGEN_POLYNOMIAL_Release_64bit.a) -target_link_libraries(KIM_lib PUBLIC "${KiLCA_lib_path}") +# KiLCA is NOT linked to KIM_lib: KIM does not use any KiLCA symbols. +# Previously PUBLIC, this caused duplicate library linking when QL-Balance +# linked both KIM_lib and kilca_lib, corrupting symbol resolution order +# and breaking KiLCA's wave solver (besselj failures, vacuum stitching errors). +# KIM_exe links KiLCA directly via slatec (which provides the math routines). target_link_libraries(KIM_lib PUBLIC GSL::gsl GSL::gslcblas) -target_include_directories(KIM_lib PUBLIC ${CMAKE_SOURCE_DIR}/KiLCA/math/) target_include_directories(KIM_lib PUBLIC ${KILCA_PATH}/math/) find_package(HDF5 COMPONENTS C Fortran HL REQUIRED) @@ -158,6 +161,7 @@ set_target_properties(KIM_exe PROPERTIES RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/install/bin/") target_link_libraries(KIM_exe KIM_lib + "${KiLCA_lib_path}" lapack slatec ) diff --git a/KIM/src/math/ddeabm/CMakeLists.txt b/KIM/src/math/ddeabm/CMakeLists.txt index b8132e7f..c2db5e8c 100644 --- a/KIM/src/math/ddeabm/CMakeLists.txt +++ b/KIM/src/math/ddeabm/CMakeLists.txt @@ -1,2 +1,15 @@ -file(GLOB ddeabm_files "*.f") +# Only include the actual DDEABM ODE-solver sources. +# The slatec error-handling helpers (i1mach, d1mach, j4save, xermsg, xerprn, +# xersve, xercnt, xerhlt, xgetua, fdump) are provided by libslatec and must +# NOT be duplicated here. ddeabm's i1mach has all DATA statements commented +# out (returns 0 for every query), which breaks slatec's DGAMLN/Bessel +# functions when it shadows the correct slatec version at link time. +set(ddeabm_files + ddeabm.f + ddes.f + dhstrt.f + dhvnrm.f + dintp.f + dsteps.f +) add_library(ddeabm "${ddeabm_files}") From 745eaea24f2751a8874768047afe01ae4001b067 Mon Sep 17 00:00:00 2001 From: Markus Markl Date: Tue, 24 Mar 2026 09:58:24 +0100 Subject: [PATCH 12/20] refactor(KIM): rename read_config and add reallocation guards MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Rename read_config → kim_read_config to avoid symbol clash when KIM is linked as a library into QL-Balance (which has its own read_config). Add deallocate-before-allocate guards in grid_m, equilibrium_m, and species_m so that KIM can be re-initialized on repeated calls without crashing on double allocation. Co-Authored-By: Claude Opus 4.6 (1M context) --- .../background_equilibrium/calculate_equil.f90 | 12 ++++++++++++ KIM/src/general/read_config.f90 | 4 ++-- KIM/src/grid/grid_mod.f90 | 16 ++++++++++++++++ KIM/src/tests/kim_init_for_test.f90 | 2 +- 4 files changed, 31 insertions(+), 3 deletions(-) diff --git a/KIM/src/background_equilibrium/calculate_equil.f90 b/KIM/src/background_equilibrium/calculate_equil.f90 index d75480c8..f9b8396e 100644 --- a/KIM/src/background_equilibrium/calculate_equil.f90 +++ b/KIM/src/background_equilibrium/calculate_equil.f90 @@ -52,6 +52,14 @@ subroutine calculate_equil(write_out) if(.not. allocated(coef)) allocate(coef(0:nder, nlagr)) + if (allocated(u)) deallocate(u) + if (allocated(B0z)) deallocate(B0z) + if (allocated(dpress_prof)) deallocate(dpress_prof) + if (allocated(B0th)) deallocate(B0th) + if (allocated(B0)) deallocate(B0) + if (allocated(hz)) deallocate(hz) + if (allocated(hth)) deallocate(hth) + if (allocated(equil_grid)) deallocate(equil_grid) allocate(u(plasma%grid_size), & B0z(plasma%grid_size), & dpress_prof(plasma%grid_size), & @@ -104,6 +112,10 @@ subroutine calculate_equil(write_out) info(1) = 1 end do + if (allocated(plasma%B0)) deallocate(plasma%B0) + if (allocated(plasma%ks)) deallocate(plasma%ks) + if (allocated(plasma%kp)) deallocate(plasma%kp) + if (allocated(plasma%om_E)) deallocate(plasma%om_E) allocate(plasma%B0(plasma%grid_size)) allocate(plasma%ks(plasma%grid_size)) allocate(plasma%kp(plasma%grid_size)) diff --git a/KIM/src/general/read_config.f90 b/KIM/src/general/read_config.f90 index d1d26f98..3841ad5c 100644 --- a/KIM/src/general/read_config.f90 +++ b/KIM/src/general/read_config.f90 @@ -1,4 +1,4 @@ -subroutine read_config +subroutine kim_read_config use config_m use constants_m @@ -107,4 +107,4 @@ subroutine read_config ! Display formatted configuration call display_kim_configuration() -end subroutine +end subroutine kim_read_config diff --git a/KIM/src/grid/grid_mod.f90 b/KIM/src/grid/grid_mod.f90 index 7eb54e4d..21c4f9c3 100644 --- a/KIM/src/grid/grid_mod.f90 +++ b/KIM/src/grid/grid_mod.f90 @@ -133,6 +133,7 @@ subroutine grid_init_equidistant(this, npts, min_val, max_val, name) this%npts_c = this%npts_b - 1 this%min_val = min_val this%max_val = max_val + if (allocated(this%name)) deallocate(this%name) allocate(character(len=len(name)) :: this%name) this%name = name @@ -157,6 +158,7 @@ subroutine grid_init(this, npts, min_val, max_val, name) this%npts_b = npts this%min_val = min_val this%max_val = max_val + if (allocated(this%name)) deallocate(this%name) allocate(character(len=len(name)) :: this%name) this%name = name @@ -203,6 +205,8 @@ subroutine grid_generate(this) real(dp), dimension(:,:), allocatable :: coef real(dp) :: recnsp + if (allocated(this%xb)) deallocate(this%xb) + if (allocated(this%xc)) deallocate(this%xc) allocate(this%xb(this%npts_b), this%xc(this%npts_c)) allocate(coef(0:nder,npoi_der)) @@ -229,6 +233,11 @@ subroutine grid_generate(this) stop endif + if (allocated(this%deriv_coef)) deallocate(this%deriv_coef) + if (allocated(this%deriv2_coef)) deallocate(this%deriv2_coef) + if (allocated(this%reint_coef)) deallocate(this%reint_coef) + if (allocated(this%ipbeg)) deallocate(this%ipbeg) + if (allocated(this%ipend)) deallocate(this%ipend) allocate(this%deriv_coef(npoi_der, this%npts_b)) allocate(this%deriv2_coef(npoi_der, this%npts_b)) allocate(this%reint_coef(npoi_der, this%npts_b)) @@ -272,6 +281,8 @@ subroutine grid_generate_equidistant(this) integer :: ipoib, ipb, ipe real(dp), dimension(:,:), allocatable :: coef + if (allocated(this%xb)) deallocate(this%xb) + if (allocated(this%xc)) deallocate(this%xc) allocate(this%xb(this%npts_b), this%xc(this%npts_c)) h = (this%max_val - this%min_val) / this%npts_b @@ -295,6 +306,11 @@ subroutine grid_generate_equidistant(this) stop endif + if (allocated(this%deriv_coef)) deallocate(this%deriv_coef) + if (allocated(this%deriv2_coef)) deallocate(this%deriv2_coef) + if (allocated(this%reint_coef)) deallocate(this%reint_coef) + if (allocated(this%ipbeg)) deallocate(this%ipbeg) + if (allocated(this%ipend)) deallocate(this%ipend) allocate(this%deriv_coef(npoi_der, this%npts_c)) allocate(this%deriv2_coef(npoi_der, this%npts_c)) allocate(this%reint_coef(npoi_der, this%npts_c)) diff --git a/KIM/src/tests/kim_init_for_test.f90 b/KIM/src/tests/kim_init_for_test.f90 index 62ada53a..43d99c7b 100644 --- a/KIM/src/tests/kim_init_for_test.f90 +++ b/KIM/src/tests/kim_init_for_test.f90 @@ -12,7 +12,7 @@ subroutine kim_init_for_test integer :: i real(dp) :: h - call read_config + call kim_read_config iprof_length = 100 From f4c962ff84627d89c50b3306fbe5ea96ebec53b5 Mon Sep 17 00:00:00 2001 From: Markus Markl Date: Tue, 24 Mar 2026 09:58:36 +0100 Subject: [PATCH 13/20] feat(KIM): expose internal APIs for library mode Add profiles_in_memory flag and set_profiles_from_arrays subroutine so QL-Balance can inject profiles directly without file I/O. Add no-write variant of postprocess_electric_field and store jpar in EBdat for extraction by the adapter. Replace SuiteSparse sparse solver with LAPACK ZGESV in the EM solver (dense direct solve, appropriate for the 2N ~ 400 system size) and fix Poisson solver method to use SuiteSparse (method=3). Co-Authored-By: Claude Opus 4.6 (1M context) --- .../background_equilibrium/species_mod.f90 | 112 ++++++++++++++++++ .../electromagnetic_solver.f90 | 50 ++++++-- KIM/src/electrostatic_poisson/fields_mod.f90 | 16 +++ .../electrostatic_poisson/solve_poisson.f90 | 2 +- KIM/src/general/KIM_init.f90 | 13 +- KIM/src/general/KIM_mod.f90 | 2 +- KIM/src/setup/config_mod.f90 | 2 + 7 files changed, 180 insertions(+), 17 deletions(-) diff --git a/KIM/src/background_equilibrium/species_mod.f90 b/KIM/src/background_equilibrium/species_mod.f90 index 16aafc37..e560da1c 100644 --- a/KIM/src/background_equilibrium/species_mod.f90 +++ b/KIM/src/background_equilibrium/species_mod.f90 @@ -122,6 +122,7 @@ subroutine allocate_plasma implicit none plasma%n_species = number_of_ion_species+1 + if (allocated(plasma%spec)) deallocate(plasma%spec) allocate(plasma%spec(0:plasma%n_species-1)) end subroutine @@ -1054,6 +1055,117 @@ subroutine read_profiles() end subroutine + subroutine set_profiles_from_arrays(r_in, n_in, Te_in, Ti_in, q_in, Er_in, npts) + !! Populate the module-level `plasma` struct from arrays passed by + !! the caller (e.g. QL-Balance). Replicates the effect of + !! read_from_text but without file I/O. Handles reallocation for + !! repeated calls (time stepping). + + use config_m, only: number_of_ion_species + use grid_m, only: r_space_dim + use kim_resonances_m, only: prop + + implicit none + + integer, intent(in) :: npts + real(dp), intent(in) :: r_in(npts) + real(dp), intent(in) :: n_in(npts) + real(dp), intent(in) :: Te_in(npts) + real(dp), intent(in) :: Ti_in(npts) + real(dp), intent(in) :: q_in(npts) + real(dp), intent(in) :: Er_in(npts) + + integer :: sigma, total_Z + + ! Force resonance re-detection on next solve + prop = .true. + + ! Deallocate derived arrays to prevent double-allocation crashes + call deallocate_plasma_derived() + + ! Set grid size + plasma%grid_size = npts + r_space_dim = npts + + ! (Re-)allocate primary arrays + call reallocate(plasma%r_grid, npts) + call reallocate(plasma%q, npts) + call reallocate(plasma%Er, npts) + + plasma%r_grid = r_in + plasma%q = q_in + plasma%Er = Er_in + + ! Electron density and temperature + call reallocate(plasma%spec(0)%n, npts) + call reallocate(plasma%spec(0)%T, npts) + plasma%spec(0)%n = n_in + plasma%spec(0)%T = Te_in + + ! Ion temperatures (all ion species get Ti) + do sigma = 1, number_of_ion_species + call reallocate(plasma%spec(sigma)%T, npts) + plasma%spec(sigma)%T = Ti_in + end do + + ! Enforce quasineutrality for ion densities + total_Z = 0 + do sigma = 1, number_of_ion_species + total_Z = total_Z + plasma%spec(sigma)%Zspec + end do + + do sigma = 1, number_of_ion_species + call reallocate(plasma%spec(sigma)%n, npts) + plasma%spec(sigma)%n = n_in * plasma%spec(sigma)%Zspec / total_Z + end do + + ! Deallocate derivative arrays to force recomputation + do sigma = 0, number_of_ion_species + if (allocated(plasma%spec(sigma)%dndr)) deallocate(plasma%spec(sigma)%dndr) + if (allocated(plasma%spec(sigma)%dTdr)) deallocate(plasma%spec(sigma)%dTdr) + end do + if (allocated(plasma%dqdr)) deallocate(plasma%dqdr) + + ! Validate units + call validate_profile_units(plasma) + + end subroutine set_profiles_from_arrays + + subroutine deallocate_plasma_derived() + !! Deallocate plasma-level and species-derived arrays that + !! calculate_equil and calculate_plasma_backs allocate with + !! bare allocate() (no guard). Without this, repeated calls + !! to these routines crash. + !! + !! NOTE: equilibrium_m arrays (B0z, B0th, etc.) must be + !! deallocated separately to avoid a circular module dependency. + !! See deallocate_equilibrium_arrays() in the adapter. + use config_m, only: number_of_ion_species + + implicit none + + integer :: sp + + ! From plasma (allocated in calculate_equil) + if (allocated(plasma%B0)) deallocate(plasma%B0) + if (allocated(plasma%ks)) deallocate(plasma%ks) + if (allocated(plasma%kp)) deallocate(plasma%kp) + if (allocated(plasma%om_E)) deallocate(plasma%om_E) + if (allocated(plasma%dqdr)) deallocate(plasma%dqdr) + + ! From plasma species (allocated in calculate_plasma_backs) + do sp = 0, number_of_ion_species + if (sp > ubound(plasma%spec, 1)) exit + if (allocated(plasma%spec(sp)%vT)) deallocate(plasma%spec(sp)%vT) + if (allocated(plasma%spec(sp)%nu)) deallocate(plasma%spec(sp)%nu) + if (allocated(plasma%spec(sp)%omega_c)) deallocate(plasma%spec(sp)%omega_c) + if (allocated(plasma%spec(sp)%rho_L)) deallocate(plasma%spec(sp)%rho_L) + if (allocated(plasma%spec(sp)%lambda_D)) deallocate(plasma%spec(sp)%lambda_D) + if (allocated(plasma%spec(sp)%z0)) deallocate(plasma%spec(sp)%z0) + end do + + end subroutine deallocate_plasma_derived + subroutine read_from_text use config_m, only: number_of_ion_species, profile_location, fstatus diff --git a/KIM/src/electromagnetic/electromagnetic_solver.f90 b/KIM/src/electromagnetic/electromagnetic_solver.f90 index 22160114..4a6092b9 100644 --- a/KIM/src/electromagnetic/electromagnetic_solver.f90 +++ b/KIM/src/electromagnetic/electromagnetic_solver.f90 @@ -43,7 +43,7 @@ subroutine run_electromagnetic(this) use kernel_adaptive_m, only: FP_fill_kernels_adaptive use grid_m, only: xl_grid, calc_mass_matrix, M_mat, theta_integration use IO_collection_m, only: write_complex_profile_abs - use poisson_solver_m, only: prepare_Laplace_matrix, dense_to_sparse + use poisson_solver_m, only: prepare_Laplace_matrix use config_m, only: output_path, collision_model, fstatus, fdebug use fields_m, only: EBdat, postprocess_electric_field, & calculate_charge_density, calculate_current_density @@ -52,7 +52,6 @@ subroutine run_electromagnetic(this) use constants_m, only: pi, sol, com_unit use KIM_kinds_m, only: dp use species_m, only: plasma - use sparse_mod, only: sparse_solveComplex_b1, sparse_solve_method, sparse_talk implicit none @@ -65,9 +64,8 @@ subroutine run_electromagnetic(this) complex(dp), allocatable :: A_block(:,:) complex(dp), allocatable :: b_block(:) - complex(dp), allocatable :: A_nz(:) - integer, allocatable :: irow(:), pcol(:) - integer :: nz_out, nrow_sp, ncol_sp + integer, allocatable :: ipiv(:) + integer :: lapack_info complex(dp), allocatable :: A_Phi(:,:) real(dp), allocatable :: laplace_perp(:,:) @@ -224,12 +222,40 @@ subroutine run_electromagnetic(this) 'Status: solving coupled Poisson-Ampere for & &(Phi, A_par) (2N =', 2*N, ')' - ! Solve - if (fdebug < 2) sparse_talk = .false. - sparse_solve_method = 1 - - call dense_to_sparse(A_block, irow, pcol, A_nz, nrow_sp, ncol_sp, nz_out) - call sparse_solveComplex_b1(nrow_sp, ncol_sp, nz_out, irow, pcol, A_nz, b_block, 0) + ! Diagnostics: boundary conditions before solve + write(*,*) '--- EM solver diagnostics (before solve) ---' + write(*,*) ' N (grid points) = ', N + write(*,*) ' alpha(1) = ', alpha(1) + write(*,*) ' alpha(N) = ', alpha(N) + write(*,*) ' Br_boundary = ', Br_boundary + write(*,*) ' Br_bnd/alpha(N) = ', Br_boundary / alpha(N) + write(*,*) ' A_block(1,1) = ', A_block(1,1) + write(*,*) ' A_block(N,N) = ', A_block(N,N) + write(*,*) ' A_block(N+1,N+1) = ', A_block(N+1,N+1) + write(*,*) ' A_block(2N,2N) = ', A_block(2*N,2*N) + write(*,*) ' b_block(1) = ', b_block(1), ' (Phi left BC)' + write(*,*) ' b_block(N) = ', b_block(N), ' (Phi right BC)' + write(*,*) ' b_block(N+1) = ', b_block(N+1), ' (Apar left BC)' + write(*,*) ' b_block(2N) = ', b_block(2*N), ' (Apar right BC)' + + ! Solve using LAPACK ZGESV (dense direct solver, fast for 2N ~ 400) + allocate(ipiv(2*N)) + call zgesv(2*N, 1, A_block, 2*N, ipiv, b_block, 2*N, lapack_info) + if (lapack_info /= 0) then + write(*,*) 'Error: ZGESV failed with info = ', lapack_info + error stop + end if + write(*,*) ' ZGESV solve succeeded (2N =', 2*N, ')' + deallocate(ipiv) + + ! Diagnostics: solution at boundaries after solve + write(*,*) '--- EM solver diagnostics (after solve) ---' + write(*,*) ' Phi(1) = ', b_block(1), ' (should be 0)' + write(*,*) ' Phi(N) = ', b_block(N), ' (should be 0)' + write(*,*) ' Apar(1) = ', b_block(N+1), ' (should be 0)' + write(*,*) ' Apar(N) = ', b_block(2*N), ' (should be Br_bnd/alpha(N))' + write(*,*) ' max|Phi| = ', maxval(abs(b_block(1:N))) + write(*,*) ' max|Apar|= ', maxval(abs(b_block(N+1:2*N))) ! Extract solution: solve gives (Phi, A_par), recover Br = alpha * A_par allocate(EBdat%Phi(N), EBdat%Apar(N), EBdat%Br(N), EBdat%r_grid(N)) @@ -249,6 +275,8 @@ subroutine run_electromagnetic(this) ! Current density call calculate_current_density(jpar, EBdat%Phi, EBdat%Br, & kernel_j_phi_llp%Kllp, kernel_j_B_llp%Kllp) + allocate(EBdat%jpar(N)) + EBdat%jpar = jpar call write_complex_profile_abs(xl_grid%xb, jpar, N, "/fields/jpar", & 'Parallel current density from self-consistent solve', 'statA/cm^2') diff --git a/KIM/src/electrostatic_poisson/fields_mod.f90 b/KIM/src/electrostatic_poisson/fields_mod.f90 index 6c4e64ef..20aee696 100644 --- a/KIM/src/electrostatic_poisson/fields_mod.f90 +++ b/KIM/src/electrostatic_poisson/fields_mod.f90 @@ -16,6 +16,7 @@ module fields_m complex(dp), allocatable :: Ez(:) complex(dp), allocatable :: Es(:) ! E "senkrecht", i.e. perpendicular to radial and parallel direction complex(dp), allocatable :: Ep(:) ! parallel to the equilibrium magnetic field + complex(dp), allocatable :: jpar(:) ! parallel current density complex(dp), allocatable :: Phi(:) complex(dp), allocatable :: Phi_e(:) complex(dp), allocatable :: Phi_i(:) @@ -425,6 +426,21 @@ subroutine postprocess_electric_field(EBdat) end subroutine + subroutine postprocess_electric_field_no_write(EBdat) + !! Same as postprocess_electric_field but without HDF5 writes. + !! Used when hdf5_output is disabled (e.g. QL-Balance adapter mode). + use species_m, only: plasma + + implicit none + + type(EBdat_t), intent(inout) :: EBdat + + call calculate_MA_field(plasma, EBdat) + call calculate_E_from_phi(EBdat) + call calculate_E_in_rsp_from_cyl(EBdat) + + end subroutine + subroutine calculate_charge_density(rho, EBdat) use KIM_kinds_m, only: dp diff --git a/KIM/src/electrostatic_poisson/solve_poisson.f90 b/KIM/src/electrostatic_poisson/solve_poisson.f90 index 8296900f..995abbe0 100644 --- a/KIM/src/electrostatic_poisson/solve_poisson.f90 +++ b/KIM/src/electrostatic_poisson/solve_poisson.f90 @@ -52,7 +52,7 @@ subroutine solve_poisson(K_rho_phi, K_rho_B, phi_sol) call write_A_matrix_sparse_check_to_file end if - sparse_solve_method = 1 + sparse_solve_method = 3 ! SuiteSparse (shared sparse lib does not support SuperLU=1) sparse_solver_option = 0 call create_rhs_vector(type_br_field, K_rho_B, b_vec) diff --git a/KIM/src/general/KIM_init.f90 b/KIM/src/general/KIM_init.f90 index 34bfaa1a..d78cd916 100644 --- a/KIM/src/general/KIM_init.f90 +++ b/KIM/src/general/KIM_init.f90 @@ -2,22 +2,27 @@ subroutine kim_init use species_m, only: read_profiles, allocate_plasma, init_plasma, plasma use IO_collection_m, only: initialize_hdf5_output, write_KIM_namelist_to_hdf5 - use config_m, only: hdf5_output + use config_m, only: hdf5_output, profiles_in_memory use profile_input_m, only: prepare_profiles implicit none - call read_config + call kim_read_config if (hdf5_output) then call initialize_hdf5_output() call write_KIM_namelist_to_hdf5() end if - call prepare_profiles() + if (.not. profiles_in_memory) then + call prepare_profiles() + end if call allocate_plasma call init_plasma(plasma) - call read_profiles() + + if (.not. profiles_in_memory) then + call read_profiles() + end if end subroutine diff --git a/KIM/src/general/KIM_mod.f90 b/KIM/src/general/KIM_mod.f90 index 7aa708dd..c7f9a96f 100644 --- a/KIM/src/general/KIM_mod.f90 +++ b/KIM/src/general/KIM_mod.f90 @@ -14,7 +14,7 @@ subroutine from_kim_factory_get_kim(type_of_run, kim_instance) implicit none - character(100), intent(in) :: type_of_run + character(*), intent(in) :: type_of_run class(kim_t), allocatable, intent(out) :: kim_instance select case(trim(type_of_run)) diff --git a/KIM/src/setup/config_mod.f90 b/KIM/src/setup/config_mod.f90 index 4148ce1f..d16e18ad 100644 --- a/KIM/src/setup/config_mod.f90 +++ b/KIM/src/setup/config_mod.f90 @@ -36,6 +36,8 @@ module config_m character(256) :: nml_config_path = "./KIM_config.nml" ! path to the namelist file + logical :: profiles_in_memory = .false. ! When true, skip file-based profile reading + logical :: rescale_density real(dp) :: number_density_rescale real(dp) :: ion_flr_scale_factor From 761c52cc76ce888795d74f17e10493a68ffc5dc4 Mon Sep 17 00:00:00 2001 From: Markus Markl Date: Tue, 24 Mar 2026 09:58:49 +0100 Subject: [PATCH 14/20] feat(QL-Balance): extend KIM integration with config and transport wiring Add KIM mode list configuration (kim_m_list, kim_n_list, kim_n_modes, kim_profiles_from_balance) to balance namelist. Extend the KIM adapter with profile injection, vacuum field loading, domain consistency checks, current density extraction, and antenna factor computation. Wire KIM adapter into the transport workflow: get_dql dispatches to KIM when wave_code='KIM', gengrid supports KIM grid setup, diff_coeffs and transp_coeffs_mod compute antenna factors from parallel current, integrate_parallel_current handles KIM jpar data, and writeData outputs KIM-specific diagnostics. Co-Authored-By: Claude Opus 4.6 (1M context) --- QL-Balance/namelists/balance_conf.nml | 4 + QL-Balance/src/base/SingleStep.f90 | 4 +- QL-Balance/src/base/control_mod.f90 | 5 + QL-Balance/src/base/diff_coeffs.f90 | 25 +- QL-Balance/src/base/gengrid.f90 | 138 +++-- QL-Balance/src/base/get_dql.f90 | 54 +- .../src/base/integrate_parallel_current.f90 | 168 ++++- QL-Balance/src/base/kim_wave_code_adapter.f90 | 586 +++++++++++++++--- QL-Balance/src/base/paramscan.f90 | 4 +- QL-Balance/src/base/read_config.f90 | 16 +- QL-Balance/src/base/time_evolution.f90 | 8 +- QL-Balance/src/base/transp_coeffs_mod.f90 | 49 ++ QL-Balance/src/base/wave_code_data_64bit.f90 | 47 +- QL-Balance/src/base/wave_code_data_mod.f90 | 1 + QL-Balance/src/base/writeData.f90 | 10 + 15 files changed, 929 insertions(+), 190 deletions(-) diff --git a/QL-Balance/namelists/balance_conf.nml b/QL-Balance/namelists/balance_conf.nml index ce901f0c..029059bc 100644 --- a/QL-Balance/namelists/balance_conf.nml +++ b/QL-Balance/namelists/balance_conf.nml @@ -57,5 +57,9 @@ type_of_run = 'SingleStep' wave_code = 'KiLCA' kim_config_path = './KIM_config.nml' + kim_profiles_from_balance = .true. + kim_n_modes = 2 + kim_m_list = 6, -6 + kim_n_list = 2, 2 ihdf5IO = 1 / diff --git a/QL-Balance/src/base/SingleStep.f90 b/QL-Balance/src/base/SingleStep.f90 index f7460bc6..634523d6 100644 --- a/QL-Balance/src/base/SingleStep.f90 +++ b/QL-Balance/src/base/SingleStep.f90 @@ -63,7 +63,8 @@ subroutine initSingleStep(this) subroutine runSingleStep(this) - use transp_coeffs_mod, only: rescale_transp_coeffs_by_ant_fac + use transp_coeffs_mod, only: rescale_transp_coeffs_by_ant_fac, & + compute_antenna_factor_from_Ipar use control_mod, only: debug_mode implicit none @@ -77,6 +78,7 @@ subroutine runSingleStep(this) if (debug_mode) write(*,*) "Debug: before get_dql" call get_dql if (debug_mode) write(*,*) "Debug: after get_dql" + call compute_antenna_factor_from_Ipar call rescale_transp_coeffs_by_ant_fac call interp_Br_Dql_at_res diff --git a/QL-Balance/src/base/control_mod.f90 b/QL-Balance/src/base/control_mod.f90 index 3298ff11..940ff9ce 100644 --- a/QL-Balance/src/base/control_mod.f90 +++ b/QL-Balance/src/base/control_mod.f90 @@ -7,6 +7,11 @@ module control_mod character(100) :: type_of_run character(len=32) :: wave_code = 'KiLCA' ! 'KiLCA' or 'KIM' character(len=1024) :: kim_config_path = './KIM_config.nml' + logical :: kim_profiles_from_balance = .true. + integer, parameter :: kim_max_modes = 100 + integer :: kim_n_modes = 0 + integer :: kim_m_list(100) = 0 + integer :: kim_n_list(100) = 0 integer :: ihdf5IO ! added: Markus Markl logical :: paramscan ! added: Markus Markl, 03.03.2021 logical :: timing_mode ! added by Markus Markl 06.04.2021 diff --git a/QL-Balance/src/base/diff_coeffs.f90 b/QL-Balance/src/base/diff_coeffs.f90 index 0072aa62..73fa7cba 100644 --- a/QL-Balance/src/base/diff_coeffs.f90 +++ b/QL-Balance/src/base/diff_coeffs.f90 @@ -135,7 +135,7 @@ subroutine calc_transport_coeffs_ornuhl(dim, vT, nu, D_11, D_12, D_21, D_22) implicit none - integer :: i + integer :: i, i_res integer, parameter :: mnmax = 3 integer, intent(in) :: dim real(dp), dimension(dim), intent(in) :: vT, nu @@ -161,6 +161,29 @@ subroutine calc_transport_coeffs_ornuhl(dim, vT, nu, D_11, D_12, D_21, D_22) !epbr_re=0.0d0 !2.d0*c*vT*real(conjg(Es)*Br) !epbr_im=0.0d0 !2.d0*c*vT*dimag(conjg(Es)*Br) + ! DIAGNOSTIC: print D_22 terms at resonance + ! Find grid index closest to resonance + i_res = 1 + do i = 2, dim + if (abs(rb(i) - r_resonant(i_mn_loop)) .lt. abs(rb(i_res) - r_resonant(i_mn_loop))) i_res = i + end do + write(*,*) '--- D_22 terms at resonance (i=', i_res, ', r=', rb(i_res), ') ---' + write(*,*) ' |Es| at res = ', abs(Es(i_res)) + write(*,*) ' |Br| at res = ', abs(Br(i_res)) + write(*,*) ' B0 at res = ', B0(i_res) + write(*,*) ' nu at res = ', nu(i_res) + write(*,*) ' vT at res = ', vT(i_res) + write(*,*) ' kp at res = ', kp(i_res) + write(*,*) ' om_E at res = ', om_E(i_res) + write(*,*) ' comfac at res = ', comfac(i_res) + write(*,*) ' epm2 at res = ', epm2(i_res), ' (c^2*|Es|^2)' + write(*,*) ' brm2 at res = ', brm2(i_res), ' (vT^2*|Br|^2)' + write(*,*) ' epbr_re at res = ', epbr_re(i_res) + write(*,*) ' --- Approximate D_22 ~ comfac*(epm2 + brm2) ---' + write(*,*) ' epm2 contrib = ', comfac(i_res)*epm2(i_res) + write(*,*) ' brm2 contrib = ', comfac(i_res)*brm2(i_res) + write(*,*) '' + x1 = kp*vT/nu x2 = -om_E/nu diff --git a/QL-Balance/src/base/gengrid.f90 b/QL-Balance/src/base/gengrid.f90 index eb7f3c62..74e90993 100644 --- a/QL-Balance/src/base/gengrid.f90 +++ b/QL-Balance/src/base/gengrid.f90 @@ -7,7 +7,7 @@ subroutine gengrid use grid_mod use plasma_parameters - use control_mod, only: debug_mode + use control_mod, only: debug_mode, wave_code, kim_n_modes, kim_m_list, kim_n_list use PolyLagrangeInterpolation implicit none @@ -174,7 +174,7 @@ subroutine prepare_resonances use resonances_mod use grid_mod, only: gg_width, gg_factor,r_resonant - use control_mod, only: ihdf5IO, debug_mode + use control_mod, only: ihdf5IO, debug_mode, wave_code, kim_n_modes, kim_m_list, kim_n_list use h5mod implicit none @@ -228,56 +228,98 @@ subroutine prepare_resonances qmin=minval(q) qmax=maxval(q) - open(iunit_res,file='flre/antenna.in') - read(iunit_res,*) - read(iunit_res,*) - read(iunit_res,*) - read(iunit_res,*) - read(iunit_res,*) - read(iunit_res,*) numres - close(iunit_res) - if (debug_mode) write(*,*) "numres = ", numres - - allocate(r_res(numres),width_res(numres),ampl_res(numres)) - allocate(r_resonant(numres)) - width_res=gg_width - ampl_res=gg_factor - allocate(m_a(numres),n_a(numres)) - allocate(m_aa(numres),n_aa(numres)) - numres_orig=numres - - open(iunit_res,file='flre/modes.in') - k=1 - jj=1 - read(iunit_res,*) a - m=nint(real(a)) - n=abs(nint(imag(a))) - m_a(k)=m - n_a(k)=n - m_aa(jj)=m - n_aa(jj)=n - r_res(k)=abs(dfloat(m)/dfloat(n)) - outer: do i=2,numres + if (trim(wave_code) == 'KIM') then + ! KIM mode: read mode list from namelist parameters + numres = kim_n_modes + if (debug_mode) write(*,*) "numres = ", numres + + allocate(r_res(numres),width_res(numres),ampl_res(numres)) + allocate(r_resonant(numres)) + width_res=gg_width + ampl_res=gg_factor + allocate(m_a(numres),n_a(numres)) + allocate(m_aa(numres),n_aa(numres)) + numres_orig=numres + + k=1 + jj=1 + m = kim_m_list(1) + n = abs(kim_n_list(1)) + m_a(k) = m + n_a(k) = n + m_aa(jj) = m + n_aa(jj) = n + r_res(k) = abs(dfloat(m)/dfloat(n)) + outer_kim: do i=2,numres + m = kim_m_list(i) + n = abs(kim_n_list(i)) + qres = abs(dfloat(m)/dfloat(n)) + if(qres.lt.qmin.or.qres.gt.qmax) cycle + jj=jj+1 + m_aa(jj) = m + n_aa(jj) = n + do j=1,k + if(m*n_a(j)-n*m_a(j).eq.0) cycle outer_kim + enddo + k=k+1 + m_a(k) = m + n_a(k) = n + r_res(k) = qres + enddo outer_kim + numres=k + else + ! KiLCA mode: read mode list from flre/antenna.in and flre/modes.in + open(iunit_res,file='flre/antenna.in') + read(iunit_res,*) + read(iunit_res,*) + read(iunit_res,*) + read(iunit_res,*) + read(iunit_res,*) + read(iunit_res,*) numres + close(iunit_res) + if (debug_mode) write(*,*) "numres = ", numres + + allocate(r_res(numres),width_res(numres),ampl_res(numres)) + allocate(r_resonant(numres)) + width_res=gg_width + ampl_res=gg_factor + allocate(m_a(numres),n_a(numres)) + allocate(m_aa(numres),n_aa(numres)) + numres_orig=numres + + open(iunit_res,file='flre/modes.in') + k=1 + jj=1 read(iunit_res,*) a - m=-abs(nint(real(a))) + m=nint(real(a)) n=abs(nint(imag(a))) - qres=abs(dfloat(m)/dfloat(n)) - !check for existence of resonant point: - if(qres.lt.qmin.or.qres.gt.qmax) cycle - !check for repeated resonance radii: - jj=jj+1 - m_aa(jj)=m - n_aa(jj)=n - do j=1,k - if(m*n_a(j)-n*m_a(j).eq.0) cycle outer - enddo - k=k+1 m_a(k)=m n_a(k)=n - r_res(k)=qres - enddo outer - close(iunit_res) - numres=k + m_aa(jj)=m + n_aa(jj)=n + r_res(k)=abs(dfloat(m)/dfloat(n)) + outer: do i=2,numres + read(iunit_res,*) a + m=-abs(nint(real(a))) + n=abs(nint(imag(a))) + qres=abs(dfloat(m)/dfloat(n)) + !check for existence of resonant point: + if(qres.lt.qmin.or.qres.gt.qmax) cycle + !check for repeated resonance radii: + jj=jj+1 + m_aa(jj)=m + n_aa(jj)=n + do j=1,k + if(m*n_a(j)-n*m_a(j).eq.0) cycle outer + enddo + k=k+1 + m_a(k)=m + n_a(k)=n + r_res(k)=qres + enddo outer + close(iunit_res) + numres=k + end if do i=1,numres qres=r_res(i) diff --git a/QL-Balance/src/base/get_dql.f90 b/QL-Balance/src/base/get_dql.f90 index dcab2dd2..112e3fb2 100644 --- a/QL-Balance/src/base/get_dql.f90 +++ b/QL-Balance/src/base/get_dql.f90 @@ -19,7 +19,8 @@ subroutine get_dql use h5mod use wave_code_data use kim_wave_code_adapter_m, only: kim_update_profiles, kim_run_for_all_modes, & - kim_get_wave_fields, kim_get_wave_vectors, kim_vac_Br, kim_Br_modes + kim_get_wave_fields, kim_get_wave_vectors, kim_vac_Br, kim_Br_modes, & + kim_get_current_densities use QLBalance_diag, only: i_mn_loop use QLBalance_kinds, only: dp use PolyLagrangeInterpolation @@ -167,6 +168,33 @@ subroutine get_dql i_mn_loop = i_mn + ! DIAGNOSTIC: compare field magnitudes between wave codes + write(*,*) '' + write(*,*) '========== FIELD DIAGNOSTICS (mode ', i_mn, ') ==========' + write(*,*) ' wave_code = ', trim(wave_code) + write(*,*) ' m/n = ', m_vals(i_mn), '/', n_vals(i_mn) + write(*,*) ' r_resonant = ', r_resonant(i_mn) + write(*,*) ' --- Peak field amplitudes ---' + write(*,*) ' max|Es| = ', maxval(abs(Es)) + write(*,*) ' max|Br| = ', maxval(abs(Br)) + write(*,*) ' max|Er| = ', maxval(abs(Er)) + write(*,*) ' max|Ep| = ', maxval(abs(Ep)) + write(*,*) ' max|Et| = ', maxval(abs(Et)) + write(*,*) ' max|Ez| = ', maxval(abs(Ez)) + write(*,*) ' --- Background quantities ---' + write(*,*) ' max B0 = ', maxval(B0) + write(*,*) ' max |kp| = ', maxval(abs(kp)) + write(*,*) ' max |ks| = ', maxval(abs(ks)) + write(*,*) ' max |om_E| = ', maxval(abs(om_E)) + write(*,*) ' max nue = ', maxval(nue) + write(*,*) ' max nui = ', maxval(nui) + write(*,*) ' --- D_22 input terms (at peak |Es| location) ---' + write(*,*) ' c^2*max|Es|^2 = ', c**2 * maxval(abs(Es))**2 + write(*,*) ' vTe^2*max|Br|^2= ', maxval(vT_e)**2 * maxval(abs(Br))**2 + write(*,*) ' 0.5/(nue*B0^2) = ', 0.5d0/(maxval(nue)*maxval(B0)**2) + write(*,*) '=================================================' + write(*,*) '' + ! TODO: add switch to choose calculation of collisionless transport coefficients. if (.false.) then call calc_transport_coeffs_collisionless(npoib, vT_e, de11, de12, de22) @@ -226,6 +254,12 @@ subroutine get_dql end select formfactor = (1.d0, 0.d0)/Br + ! DIAGNOSTIC: vacuum Br values + write(*,*) ' --- Vacuum Br ---' + write(*,*) ' |Br_vac| at r_min = ', abs(Br(1)), ' at r=', r(1) + write(*,*) ' |Br_vac| at r_max = ', abs(Br(dim_r)), ' at r=', r(dim_r) + write(*,*) ' max|Br_vac| = ', maxval(abs(Br)) + ! In case that the tmhd code uses double sided Fourier series, it ! must be set to 4.0d0, since the factor 2.0d0 should occur in the fields. spec_weight = 1.0d0 @@ -252,6 +286,13 @@ subroutine get_dql end select formfactor = Br * formfactor + ! DIAGNOSTIC: Br_total and formfactor = Br_total/Br_vac + write(*,*) ' --- Br_total (self-consistent / full response) ---' + write(*,*) ' |Br_total| at r_max = ', abs(Br(dim_r)), ' at r=', r(dim_r) + write(*,*) ' max|Br_total| = ', maxval(abs(Br)) + write(*,*) ' |formfactor| at r_max= ', abs(formfactor(dim_r)) + write(*,*) ' max|formfactor| = ', maxval(abs(formfactor)) + ! todo: interpolate formfactor at resonant surface and write out. This is Brtot/Brvac at the resonant surface if (trim(type_of_run) .eq. "TimeEvolution") then !TODO: this is a very bad solution... Make it better if (time_ind > 0) then @@ -273,17 +314,10 @@ subroutine get_dql call get_current_densities_from_wave_code(flre_cd_ptr(i_mn), dim_r, r, & m_vals(i_mn), n_vals(i_mn), Jri, Jsi, Jpi, Jre, Jse, Jpe) case ('KIM') - ! KIM only provides total jpar; approximate Jpe with it, zero rest - ! TODO: improve when KIM provides per-species currents - Jpe = (0.0d0, 0.0d0) - Jpi = (0.0d0, 0.0d0) - Jre = (0.0d0, 0.0d0) - Jri = (0.0d0, 0.0d0) - Jse = (0.0d0, 0.0d0) - Jsi = (0.0d0, 0.0d0) + call kim_get_current_densities(i_mn) end select - call integrate_parallel_current(dim_r, r, Jpe, Jpi, Ipar) + call integrate_parallel_current(dim_r, r, Jpe, Jpi, r_resonant(i_mn), Ipar) end do diff --git a/QL-Balance/src/base/integrate_parallel_current.f90 b/QL-Balance/src/base/integrate_parallel_current.f90 index 6b138401..934ed50a 100644 --- a/QL-Balance/src/base/integrate_parallel_current.f90 +++ b/QL-Balance/src/base/integrate_parallel_current.f90 @@ -1,39 +1,167 @@ -subroutine integrate_parallel_current(dim_r, r, Jpe, Jpi, Ipar) - +subroutine integrate_parallel_current(dim_r, r, Jpe, Jpi, r_res_in, Ipar) +! +! Integrate (Jpi + Jpe) * r * dr over a Gaussian-fitted layer around +! the resonant surface to compute the cylindrical shielding current Ipar. +! +! KiLCA convention: total current = ions + electrons (charge sign is +! already built into the conductivity tensor, so species add directly). +! This matches KiLCA's internal spec=2 = spec=0 + spec=1. +! +! Method (matches KiLCA postprocessor): +! 1. Compute Jpar_net = Jpi + Jpe +! 2. Find peak |Jpar_net| within ±5 cm of r_res +! 3. Fit Gaussian to |Jpar_net| via log-linearization +! 4. Integration range = r_res ± d/2 where d = 5*sigma +! 5. Trapezoidal integration within that range +! use QLBalance_kinds, only: dp - use grid_mod, only: r_resonant - use resonances_mod, only: numres implicit none integer, intent(in) :: dim_r real(dp), dimension(dim_r), intent(in) :: r complex(dp), dimension(dim_r), intent(in) :: Jpe, Jpi + real(dp), intent(in) :: r_res_in complex(dp), intent(out) :: Ipar - real(dp) :: int_range_1, int_range_2 + ! Local variables + complex(dp), dimension(dim_r) :: Jpar_net + real(dp), dimension(dim_r) :: absJ + real(dp) :: search_half_width, threshold, peak_val + real(dp) :: sigma, d_layer, r_lo, r_hi + integer :: ipeak, ipoi + integer :: ind_lo, ind_hi + + ! Gaussian fit variables + integer :: nfit + real(dp) :: xi, yi, sum_x, sum_y, sum_xx, sum_xy, denom, a1 + real(dp), parameter :: sigma_min = 0.01d0 ! cm + real(dp), parameter :: sigma_max = 5.0d0 ! cm + real(dp), parameter :: d_min = 0.1d0 ! cm + real(dp), parameter :: d_max = 10.0d0 ! cm + real(dp), parameter :: fallback_half = 2.5d0 ! cm fallback half-width + logical :: fit_ok + + ! Step 1: Compute total parallel current density (KiLCA: J_total = J_i + J_e) + Jpar_net = Jpi + Jpe + absJ = abs(Jpar_net) + + ! Step 2: Find peak |Jpar_net| within ±5 cm of r_res + search_half_width = 5.0d0 + ipeak = 0 + peak_val = 0.0d0 + do ipoi = 1, dim_r + if (abs(r(ipoi) - r_res_in) <= search_half_width) then + if (absJ(ipoi) > peak_val) then + peak_val = absJ(ipoi) + ipeak = ipoi + end if + end if + end do + + ! If no peak found or peak is zero, return zero + if (ipeak == 0 .or. peak_val < 1.0d-30) then + Ipar = (0.0d0, 0.0d0) + write(*,*) '[integrate_parallel_current] WARNING: no significant Jpar peak found near r_res=', r_res_in + return + end if - integer :: imn, ipoi, ind_r1, ind_r2 + ! Step 3: Gaussian fit via log-linearization + ! log(|Jpar(r)|) = a0 + a1*(r - r_res)^2 + ! where a1 = -1/(2*sigma^2) + threshold = 0.1d0 * peak_val - Ipar = 0.0d0 - ind_r1 = 1 - ind_r2 = dim_r - int_range_1 = 2.5d0 - int_range_2 = 2.5d0 + ! Count fitting points + nfit = 0 + do ipoi = 1, dim_r + if (abs(r(ipoi) - r_res_in) <= search_half_width .and. absJ(ipoi) > threshold) then + nfit = nfit + 1 + end if + end do + + fit_ok = .false. + if (nfit >= 3) then + ! Linear least squares: y = a0 + a1*x where x=(r-r_res)^2, y=log(|Jpar|) + sum_x = 0.0d0 + sum_y = 0.0d0 + sum_xx = 0.0d0 + sum_xy = 0.0d0 - do imn = 1, numres do ipoi = 1, dim_r - if (r(ipoi) < r_resonant(imn) - int_range_1) then - ind_r1 = ipoi - else if (r(ipoi) > r_resonant(imn) + int_range_2) then - ind_r2 = ipoi - exit + if (abs(r(ipoi) - r_res_in) <= search_half_width .and. absJ(ipoi) > threshold) then + xi = (r(ipoi) - r_res_in)**2 + yi = log(absJ(ipoi)) + sum_x = sum_x + xi + sum_y = sum_y + yi + sum_xx = sum_xx + xi * xi + sum_xy = sum_xy + xi * yi end if end do - do ipoi=1, dim_r-1 - Ipar = Ipar + (Jpi(ipoi)-Jpe(ipoi) + Jpi(ipoi+1) - Jpe(ipoi+1)) * r(ipoi) * (r(ipoi+1) - r(ipoi)) / 2.0d0 - end do + denom = dble(nfit) * sum_xx - sum_x * sum_x + if (abs(denom) > 1.0d-30) then + a1 = (dble(nfit) * sum_xy - sum_x * sum_y) / denom + if (a1 < 0.0d0) then + sigma = sqrt(-1.0d0 / (2.0d0 * a1)) + if (sigma >= sigma_min .and. sigma <= sigma_max) then + fit_ok = .true. + end if + end if + end if + end if + + ! Step 4: Set integration range + if (fit_ok) then + d_layer = 5.0d0 * sigma + ! Clamp to safety bounds + d_layer = max(d_min, min(d_max, d_layer)) + r_lo = r_res_in - d_layer / 2.0d0 + r_hi = r_res_in + d_layer / 2.0d0 + + write(*,*) '[integrate_parallel_current] Gaussian fit succeeded:' + write(*,*) ' r_res = ', r_res_in, ' cm' + write(*,*) ' sigma = ', sigma, ' cm' + write(*,*) ' d=5sig = ', d_layer, ' cm' + write(*,*) ' range = [', r_lo, ',', r_hi, '] cm' + write(*,*) ' nfit = ', nfit + else + ! Fallback: fixed ±2.5 cm window + r_lo = r_res_in - fallback_half + r_hi = r_res_in + fallback_half + + write(*,*) '[integrate_parallel_current] WARNING: Gaussian fit failed, using fallback +-2.5 cm' + write(*,*) ' r_res = ', r_res_in, ' cm' + write(*,*) ' nfit = ', nfit + if (nfit >= 3) write(*,*) ' a1 = ', a1 + write(*,*) ' range = [', r_lo, ',', r_hi, '] cm' + end if + + ! Clamp to grid bounds + r_lo = max(r_lo, r(1)) + r_hi = min(r_hi, r(dim_r)) + + ! Find index bounds + ind_lo = 1 + ind_hi = dim_r + do ipoi = 1, dim_r + if (r(ipoi) < r_lo) then + ind_lo = ipoi + end if + if (r(ipoi) > r_hi) then + ind_hi = ipoi + exit + end if + end do + + ! Step 5: Trapezoidal integration of (Jpi + Jpe) * r * dr within bounds + Ipar = (0.0d0, 0.0d0) + do ipoi = ind_lo, ind_hi - 1 + Ipar = Ipar + (Jpar_net(ipoi) * r(ipoi) + Jpar_net(ipoi+1) * r(ipoi+1)) & + * (r(ipoi+1) - r(ipoi)) / 2.0d0 end do + ! Step 6: Diagnostic output + write(*,*) ' |Ipar| = ', abs(Ipar) + write(*,*) ' Re(Ipar)= ', real(Ipar), ' Im(Ipar)= ', aimag(Ipar) + end subroutine diff --git a/QL-Balance/src/base/kim_wave_code_adapter.f90 b/QL-Balance/src/base/kim_wave_code_adapter.f90 index 9b7fe319..c4330ed1 100644 --- a/QL-Balance/src/base/kim_wave_code_adapter.f90 +++ b/QL-Balance/src/base/kim_wave_code_adapter.f90 @@ -7,12 +7,19 @@ module kim_wave_code_adapter_m !! CMake to avoid collision with the QL-Balance resonances_mod. !! All KIM symbols are renamed on import to avoid name clashes. - use species_m, only: kim_plasma => plasma - use equilibrium_m, only: kim_B0 => B0 + use species_m, only: kim_plasma => plasma, & + set_profiles_from_arrays, deallocate_plasma_derived + use equilibrium_m, only: kim_B0 => B0, & + eq_B0z => B0z, eq_B0th => B0th, & + eq_hz => hz, eq_hth => hth, eq_equil_grid => equil_grid, & + eq_u => u, eq_dpress_prof => dpress_prof use config_m, only: kim_type_of_run => type_of_run, & - kim_nml_config_path => nml_config_path + kim_nml_config_path => nml_config_path, & + kim_profiles_in_memory => profiles_in_memory, & + kim_hdf5_output => hdf5_output use setup_m, only: kim_m_mode => m_mode, kim_n_mode => n_mode - use grid_m, only: kim_xl_grid => xl_grid, kim_rg_grid => rg_grid + use grid_m, only: kim_xl_grid => xl_grid, kim_rg_grid => rg_grid, & + kim_r_min => r_min, kim_r_plas => r_plas use fields_m, only: EBdat use kim_base_m, only: kim_t use kim_mod_m, only: from_kim_factory_get_kim @@ -27,11 +34,19 @@ module kim_wave_code_adapter_m public :: kim_get_wave_vectors public :: kim_get_background_magnetic_fields public :: kim_get_collision_frequencies + public :: kim_load_vacuum_fields + public :: kim_check_domain_consistency + public :: kim_get_current_densities public :: interp_complex_profile ! exposed for testing !! Module-level KIM solver instance (reused across calls) class(kim_t), allocatable :: kim_instance + !! Actual KIM field grid boundary radius (xl_grid%xb(N)). + !! This may differ from r_plas due to non-equidistant grid construction. + !! Used for vacuum continuation cutoff and Br boundary condition. + real(8) :: kim_r_boundary = 0.0d0 + !! Vacuum Br placeholder (filled properly in Task 7) complex(8), allocatable, public :: kim_vac_Br(:) @@ -50,19 +65,33 @@ module kim_wave_code_adapter_m real(8), allocatable :: kim_kp_modes(:,:) real(8), allocatable :: kim_ks_modes(:,:) + !! Per-mode stored current density (nrad, dim_mn) + complex(8), allocatable :: kim_jpar_modes(:,:) + contains subroutine kim_initialize(nrad, r_grid) !! Initialize KIM backend: read config, profiles, grids. !! Populate wave_code_data module arrays with background - !! quantities by running the electrostatic solver for the + !! quantities by running the electromagnetic init for the !! first mode to generate equilibrium and plasma backgrounds. - use control_mod, only: kim_config_path - use wave_code_data, only: dim_mn, m_vals, n_vals, & + !! + !! Two paths: + !! kim_profiles_from_balance = .true. (Path A, default): + !! Caller already loaded QL-Balance profiles into wave_code_data + !! and set mode arrays. KIM reads config only; profiles injected + !! in memory via set_profiles_from_arrays. + !! kim_profiles_from_balance = .false. (Path B): + !! KIM reads its own files, adapter reads modes.in, extracts + !! everything — original behavior. + use control_mod, only: kim_config_path, kim_profiles_from_balance + use IO_collection_m, only: deinitialize_hdf5_output + use wave_code_data, only: dim_mn, m_vals, n_vals, dim_r, & r => r, q => q, n => n, Te => Te, Ti => Ti, & Vth => Vth, Vz => Vz, dPhi0 => dPhi0, & kp => kp, ks => ks, om_E => om_E, & B0 => B0, nue => nue, nui => nui, & + B0t => B0t, B0z => B0z, & flre_path implicit none @@ -70,125 +99,190 @@ subroutine kim_initialize(nrad, r_grid) integer, intent(in) :: nrad real(8), intent(in) :: r_grid(nrad) - integer :: kim_npts + integer :: kim_npts, nrad_inside, i real(8), allocatable :: kim_r(:) real(8), allocatable :: work_old(:), work_new(:) ! ----------------------------------------------------------- ! 1. Set KIM config path and initialize KIM backend ! ----------------------------------------------------------- - ! kim_nml_config_path (from config_m) tells KIM where its - ! namelist file lives. Controlled by kim_config_path in - ! balance_conf.nml (default: ./KIM_config.nml). kim_nml_config_path = trim(kim_config_path) - call kim_init() + if (kim_profiles_from_balance) then + ! Path A: profiles come from QL-Balance in memory + kim_profiles_in_memory = .true. + end if - ! ----------------------------------------------------------- - ! 2. Read (m,n) modes from modes.in via existing utility - ! ----------------------------------------------------------- - call read_antenna_modes(flre_path) + call kim_init() - ! ----------------------------------------------------------- - ! 3. Allocate wave_code_data arrays on QL-Balance grid - ! ----------------------------------------------------------- - call allocate_wave_code_data(nrad, r_grid) + ! Disable KIM HDF5 output for QL-Balance integration. + ! We only need in-memory fields; HDF5 writes conflict on re-init. + if (kim_hdf5_output) then + call deinitialize_hdf5_output() + end if + kim_hdf5_output = .false. + + if (kim_profiles_from_balance) then + ! ----------------------------------------------------------- + ! Path A: inject QL-Balance profiles into KIM's plasma struct + ! ----------------------------------------------------------- + ! Grid coverage check + if (r(1) > kim_r_min + 1.0d-6 .or. r(dim_r) < kim_r_plas - 1.0d-6) then + write(*,*) 'WARNING: QL-Balance grid [', r(1), ',', r(dim_r), & + '] does not fully cover KIM range [', kim_r_min, ',', kim_r_plas, ']' + end if + + ! Inject profiles: note dPhi0 = -Er, so pass -dPhi0 as Er + call deallocate_equilibrium_arrays() + call set_profiles_from_arrays(r, n, Te, Ti, q, -dPhi0, dim_r) + + ! Initialize KIM grids and equilibrium for the first mode + kim_type_of_run = 'electromagnetic' + kim_m_mode = m_vals(1) + kim_n_mode = n_vals(1) + + call from_kim_factory_get_kim(kim_type_of_run, kim_instance) + call kim_instance%init() - ! ----------------------------------------------------------- - ! 4. Run KIM electrostatic solver for the first mode to - ! generate equilibrium, grids, and plasma backgrounds - ! ----------------------------------------------------------- - kim_type_of_run = 'electrostatic' - kim_m_mode = m_vals(1) - kim_n_mode = n_vals(1) + ! Store actual KIM grid boundary (may differ from r_plas) + kim_r_boundary = kim_xl_grid%xb(kim_xl_grid%npts_b) + write(*,*) ' KIM xl_grid boundary: r_plas=', kim_r_plas, & + ' xl_grid(N)=', kim_r_boundary + + ! Extract derived quantities from KIM back to wave_code_data + ! Do NOT overwrite n, Te, Ti, q, dPhi0 (they came from QL-Balance) + ! Interpolate only up to r_plas; clamp to boundary value beyond. + kim_npts = kim_plasma%grid_size + allocate(kim_r(kim_npts), work_old(kim_npts)) + + kim_r = kim_plasma%r_grid + + ! Find last balance grid index within KIM domain + nrad_inside = nrad + do i = 1, nrad + if (r(i) > kim_r_plas) then + nrad_inside = i - 1 + exit + end if + end do + + ! Background magnetic field B0 + work_old = kim_B0 + call interp_profile(kim_npts, kim_r, work_old, nrad_inside, r(1:nrad_inside), B0(1:nrad_inside)) + B0(nrad_inside+1:nrad) = B0(nrad_inside) + + ! B0 toroidal and poloidal components + call interp_B0_components_clamped(kim_npts, kim_r, nrad, r, nrad_inside) + + ! Wave vectors kp, ks + work_old = kim_plasma%kp + call interp_profile(kim_npts, kim_r, work_old, nrad_inside, r(1:nrad_inside), kp(1:nrad_inside)) + kp(nrad_inside+1:nrad) = kp(nrad_inside) + + work_old = kim_plasma%ks + call interp_profile(kim_npts, kim_r, work_old, nrad_inside, r(1:nrad_inside), ks(1:nrad_inside)) + ks(nrad_inside+1:nrad) = ks(nrad_inside) + + ! ExB drift frequency + work_old = kim_plasma%om_E + call interp_profile(kim_npts, kim_r, work_old, nrad_inside, r(1:nrad_inside), om_E(1:nrad_inside)) + om_E(nrad_inside+1:nrad) = om_E(nrad_inside) + + ! Collision frequencies + work_old = kim_plasma%spec(0)%nu + call interp_profile(kim_npts, kim_r, work_old, nrad_inside, r(1:nrad_inside), nue(1:nrad_inside)) + nue(nrad_inside+1:nrad) = nue(nrad_inside) + + work_old = kim_plasma%spec(1)%nu + call interp_profile(kim_npts, kim_r, work_old, nrad_inside, r(1:nrad_inside), nui(1:nrad_inside)) + nui(nrad_inside+1:nrad) = nui(nrad_inside) + + ! Vth and Vz: KIM does not store these + Vth = 0.0d0 + Vz = 0.0d0 + + deallocate(kim_r, work_old) + + else + ! ----------------------------------------------------------- + ! Path B: KIM reads its own files (original behavior) + ! ----------------------------------------------------------- + call read_antenna_modes(flre_path) + call allocate_wave_code_data(nrad, r_grid) + + ! Initialize KIM grids and equilibrium for the first mode + kim_type_of_run = 'electromagnetic' + kim_m_mode = m_vals(1) + kim_n_mode = n_vals(1) + + call from_kim_factory_get_kim(kim_type_of_run, kim_instance) + call kim_instance%init() - call from_kim_factory_get_kim(kim_type_of_run, kim_instance) - call kim_instance%init() - call kim_instance%run() + ! Store actual KIM grid boundary (may differ from r_plas) + kim_r_boundary = kim_xl_grid%xb(kim_xl_grid%npts_b) - ! ----------------------------------------------------------- - ! 5. Extract and interpolate KIM backgrounds onto the - ! QL-Balance radial grid - ! ----------------------------------------------------------- - kim_npts = kim_plasma%grid_size - allocate(kim_r(kim_npts)) - kim_r = kim_plasma%r_grid + ! Extract everything from KIM onto QL-Balance grid + kim_npts = kim_plasma%grid_size + allocate(kim_r(kim_npts), work_old(kim_npts), work_new(nrad)) - ! Temporary work arrays for interpolation - allocate(work_old(kim_npts), work_new(nrad)) + kim_r = kim_plasma%r_grid - ! Safety factor q - work_old = kim_plasma%q - call interp_profile(kim_npts, kim_r, work_old, nrad, r, q) + work_old = kim_plasma%q + call interp_profile(kim_npts, kim_r, work_old, nrad, r, q) - ! Parallel wavenumber kp - work_old = kim_plasma%kp - call interp_profile(kim_npts, kim_r, work_old, nrad, r, kp) + work_old = kim_plasma%kp + call interp_profile(kim_npts, kim_r, work_old, nrad, r, kp) - ! Perpendicular wavenumber ks - work_old = kim_plasma%ks - call interp_profile(kim_npts, kim_r, work_old, nrad, r, ks) + work_old = kim_plasma%ks + call interp_profile(kim_npts, kim_r, work_old, nrad, r, ks) - ! ExB drift frequency om_E - work_old = kim_plasma%om_E - call interp_profile(kim_npts, kim_r, work_old, nrad, r, om_E) + work_old = kim_plasma%om_E + call interp_profile(kim_npts, kim_r, work_old, nrad, r, om_E) - ! Equilibrium radial electric field -> dPhi0 - ! Sign convention: dPhi0 = -Er - work_old = -kim_plasma%Er - call interp_profile(kim_npts, kim_r, work_old, nrad, r, dPhi0) + ! dPhi0 = -Er + work_old = -kim_plasma%Er + call interp_profile(kim_npts, kim_r, work_old, nrad, r, dPhi0) - ! Electron density [1/cm^3] - work_old = kim_plasma%spec(0)%n - call interp_profile(kim_npts, kim_r, work_old, nrad, r, n) + work_old = kim_plasma%spec(0)%n + call interp_profile(kim_npts, kim_r, work_old, nrad, r, n) - ! Electron temperature [eV] - work_old = kim_plasma%spec(0)%T - call interp_profile(kim_npts, kim_r, work_old, nrad, r, Te) + work_old = kim_plasma%spec(0)%T + call interp_profile(kim_npts, kim_r, work_old, nrad, r, Te) - ! Ion temperature [eV] - work_old = kim_plasma%spec(1)%T - call interp_profile(kim_npts, kim_r, work_old, nrad, r, Ti) + work_old = kim_plasma%spec(1)%T + call interp_profile(kim_npts, kim_r, work_old, nrad, r, Ti) - ! Electron collision frequency - work_old = kim_plasma%spec(0)%nu - call interp_profile(kim_npts, kim_r, work_old, nrad, r, nue) + work_old = kim_plasma%spec(0)%nu + call interp_profile(kim_npts, kim_r, work_old, nrad, r, nue) - ! Ion collision frequency - work_old = kim_plasma%spec(1)%nu - call interp_profile(kim_npts, kim_r, work_old, nrad, r, nui) + work_old = kim_plasma%spec(1)%nu + call interp_profile(kim_npts, kim_r, work_old, nrad, r, nui) - ! Background magnetic field B0 (total) - ! equilibrium_m::B0 is on the plasma r_grid - work_old = kim_B0 - call interp_profile(kim_npts, kim_r, work_old, nrad, r, B0) + work_old = kim_B0 + call interp_profile(kim_npts, kim_r, work_old, nrad, r, B0) - ! B0 toroidal and poloidal components from equilibrium_m - call interp_B0_components(kim_npts, kim_r, nrad, r) + call interp_B0_components(kim_npts, kim_r, nrad, r) - ! Poloidal rotation velocity Vth - ! KIM does not store Vth; set to zero (can be updated later) - Vth = 0.0d0 + Vth = 0.0d0 + Vz = 0.0d0 - ! Toroidal rotation Vz - ! KIM does not store Vz directly in plasma struct; - ! set to zero (Er already accounts for force balance) - Vz = 0.0d0 + deallocate(kim_r, work_old, work_new) + end if - ! ----------------------------------------------------------- - ! 6. Set up vacuum Br placeholder for form factor - ! (filled properly in Task 7) - ! ----------------------------------------------------------- + ! Vacuum Br placeholder if (allocated(kim_vac_Br)) deallocate(kim_vac_Br) allocate(kim_vac_Br(nrad)) - kim_vac_Br = (1.0d0, 0.0d0) ! placeholder: unit vacuum Br - - deallocate(kim_r, work_old, work_new) + kim_vac_Br = (1.0d0, 0.0d0) write(*, *) "KIM adapter: initialization complete" write(*, *) " Balance grid points: ", nrad write(*, *) " KIM grid points: ", kim_npts write(*, *) " Number of modes: ", dim_mn + if (kim_profiles_from_balance) then + write(*, *) " Profile source: QL-Balance (in-memory)" + else + write(*, *) " Profile source: KIM files" + end if end subroutine kim_initialize @@ -221,6 +315,35 @@ subroutine interp_B0_components(kim_npts, kim_r, nbal, bal_r) end subroutine interp_B0_components + ! --------------------------------------------------------------- + ! Helper: interpolate B0 components, clamped at r_plas + ! --------------------------------------------------------------- + subroutine interp_B0_components_clamped(kim_npts, kim_r, nbal, bal_r, nbal_inside) + !! Like interp_B0_components but only interpolates up to + !! nbal_inside points; clamps the rest to the boundary value. + use equilibrium_m, only: B0z_equil => B0z, B0th_equil => B0th + use wave_code_data, only: B0t, B0z + + implicit none + + integer, intent(in) :: kim_npts, nbal, nbal_inside + real(8), intent(in) :: kim_r(kim_npts), bal_r(nbal) + real(8), allocatable :: work(:) + + allocate(work(kim_npts)) + + work = B0z_equil + call interp_profile(kim_npts, kim_r, work, nbal_inside, bal_r(1:nbal_inside), B0z(1:nbal_inside)) + B0z(nbal_inside+1:nbal) = B0z(nbal_inside) + + work = B0th_equil + call interp_profile(kim_npts, kim_r, work, nbal_inside, bal_r(1:nbal_inside), B0t(1:nbal_inside)) + B0t(nbal_inside+1:nbal) = B0t(nbal_inside) + + deallocate(work) + + end subroutine interp_B0_components_clamped + subroutine kim_run_for_all_modes() !! Run KIM electrostatic solver for each (m,n) mode and !! store the resulting fields in per-mode arrays. @@ -230,7 +353,12 @@ subroutine kim_run_for_all_modes() !! grid. kim_get_wave_fields(i_mn) copies from these !! arrays into the wave_code_data module scalars. use wave_code_data, only: dim_mn, m_vals, n_vals, & - dim_r, bal_r => r + dim_r, bal_r => r, & + wcd_n => n, wcd_Te => Te, wcd_Ti => Ti, & + wcd_q => q, wcd_dPhi0 => dPhi0, & + wcd_kp => kp, wcd_ks => ks, wcd_om_E => om_E, & + wcd_B0 => B0, wcd_nue => nue, wcd_nui => nui + use control_mod, only: kim_profiles_from_balance implicit none @@ -238,6 +366,17 @@ subroutine kim_run_for_all_modes() real(8), allocatable :: kim_r(:) real(8), allocatable :: kim_plasma_r(:), work_real(:) + ! ------------------------------------------------------- + ! 0. Re-inject time-evolved QL-Balance profiles into KIM + ! ------------------------------------------------------- + ! NOTE: Profile re-injection is only needed in TimeEvolution mode + ! where profiles change between calls. For SingleStep, skip it + ! because set_profiles_from_arrays deallocates derived quantities + ! (B0, kp, ks, vT, nu, x1, x2, etc.) that the solver needs. + ! Re-injection for TimeEvolution would require also calling + ! calculate_equil and set_plasma_quantities to recompute derived + ! quantities on the rg_grid. + ! ------------------------------------------------------- ! 1. (Re-)allocate per-mode storage ! ------------------------------------------------------- @@ -249,6 +388,7 @@ subroutine kim_run_for_all_modes() if (allocated(kim_Br_modes)) deallocate(kim_Br_modes) if (allocated(kim_kp_modes)) deallocate(kim_kp_modes) if (allocated(kim_ks_modes)) deallocate(kim_ks_modes) + if (allocated(kim_jpar_modes)) deallocate(kim_jpar_modes) allocate(kim_Es_modes(dim_r, dim_mn)) allocate(kim_Ep_modes(dim_r, dim_mn)) @@ -258,6 +398,7 @@ subroutine kim_run_for_all_modes() allocate(kim_Br_modes(dim_r, dim_mn)) allocate(kim_kp_modes(dim_r, dim_mn)) allocate(kim_ks_modes(dim_r, dim_mn)) + allocate(kim_jpar_modes(dim_r, dim_mn)) kim_Es_modes = (0.0d0, 0.0d0) kim_Ep_modes = (0.0d0, 0.0d0) @@ -267,6 +408,7 @@ subroutine kim_run_for_all_modes() kim_Br_modes = (0.0d0, 0.0d0) kim_kp_modes = 0.0d0 kim_ks_modes = 0.0d0 + kim_jpar_modes = (0.0d0, 0.0d0) ! ------------------------------------------------------- ! 2. Loop over modes: solve and store @@ -280,14 +422,14 @@ subroutine kim_run_for_all_modes() ! Clean up EBdat from previous solve call deallocate_EBdat() - ! Create a fresh solver, init grids/equilibrium, solve - call from_kim_factory_get_kim('electrostatic', kim_instance) - call kim_instance%init() + ! Run electromagnetic solver (grids/equilibrium from kim_initialize) + ! NOTE: For multi-mode support, plasma kp/ks/om_E would need + ! updating here for modes 2+. Currently correct for single-mode. call kim_instance%run() - ! run() calls postprocess_electric_field which computes - ! Es, Ep, Er, Etheta, Ez from the Poisson-solved Phi. - ! EBdat%Br is set by set_Br_field inside solve_poisson. + ! run() fills kernels, solves coupled Poisson-Ampere, and calls + ! postprocess_electric_field which computes Es, Ep, Er, Etheta, Ez. + ! EBdat%Br = alpha * A_par (self-consistent). ! Interpolate KIM fields (on xl_grid%xb) onto ! the QL-Balance radial grid and store per-mode. @@ -295,6 +437,14 @@ subroutine kim_run_for_all_modes() allocate(kim_r(kim_npts)) kim_r = EBdat%r_grid + if (i_mn == 1) then + write(*,*) ' KIM field grid: r_min=', kim_r(1), & + ' r_max=', kim_r(kim_npts), ' npts=', kim_npts + write(*,*) ' KIM r_plas =', kim_r_plas + write(*,*) ' KIM |Br| at last grid pt =', abs(EBdat%Br(kim_npts)) + write(*,*) ' KIM |Br| at 2nd-last =', abs(EBdat%Br(kim_npts-1)) + end if + ! Es (perpendicular E field in rsp coordinates) call interp_complex_profile(kim_npts, kim_r, EBdat%Es, & dim_r, bal_r, kim_Es_modes(:, i_mn)) @@ -319,8 +469,19 @@ subroutine kim_run_for_all_modes() call interp_complex_profile(kim_npts, kim_r, EBdat%Br, & dim_r, bal_r, kim_Br_modes(:, i_mn)) + ! jpar (parallel current density) + if (allocated(EBdat%jpar)) then + call interp_complex_profile(kim_npts, kim_r, EBdat%jpar, & + dim_r, bal_r, kim_jpar_modes(:, i_mn)) + end if + deallocate(kim_r) + ! Apply vacuum continuation beyond r_plas: + ! KIM fields are only valid inside the plasma domain. + ! Beyond r_plas, use the vacuum Br from KiLCA and zero E fields. + call apply_vacuum_continuation(i_mn, dim_r, bal_r) + ! kp and ks (mode-dependent wave vectors on plasma grid) kim_plasma_npts = kim_plasma%grid_size allocate(kim_plasma_r(kim_plasma_npts)) @@ -343,6 +504,34 @@ subroutine kim_run_for_all_modes() end do + ! ------------------------------------------------------- + ! 3. Re-extract updated derived quantities when profiles + ! came from QL-Balance (they change with equilibrium) + ! ------------------------------------------------------- + if (kim_profiles_from_balance) then + kim_plasma_npts = kim_plasma%grid_size + allocate(kim_plasma_r(kim_plasma_npts), work_real(kim_plasma_npts)) + kim_plasma_r = kim_plasma%r_grid + + work_real = kim_plasma%om_E + call interp_profile(kim_plasma_npts, kim_plasma_r, & + work_real, dim_r, bal_r, wcd_om_E) + + work_real = kim_plasma%spec(0)%nu + call interp_profile(kim_plasma_npts, kim_plasma_r, & + work_real, dim_r, bal_r, wcd_nue) + + work_real = kim_plasma%spec(1)%nu + call interp_profile(kim_plasma_npts, kim_plasma_r, & + work_real, dim_r, bal_r, wcd_nui) + + work_real = kim_B0 + call interp_profile(kim_plasma_npts, kim_plasma_r, & + work_real, dim_r, bal_r, wcd_B0) + + deallocate(kim_plasma_r, work_real) + end if + write(*, *) "KIM adapter: all modes solved" end subroutine kim_run_for_all_modes @@ -359,12 +548,16 @@ subroutine kim_update_profiles() !! params_b(4,:) = Ti [erg] !! !! wave_code_data stores Te, Ti in [eV], Vz in [cm/s]. + !! + !! Note: dPhi0 is only updated from Ercov during time evolution. + !! For SingleStep runs, the input Er profile is preserved. use wave_code_data, only: dim_r, & wcd_n => n, wcd_Te => Te, wcd_Ti => Ti, & wcd_Vz => Vz, wcd_dPhi0 => dPhi0 use plasma_parameters, only: params_b use baseparam_mod, only: ev, rtor use grid_mod, only: Ercov + use control_mod, only: type_of_run implicit none @@ -375,9 +568,16 @@ subroutine kim_update_profiles() wcd_Te(k) = params_b(3, k) / ev ! erg -> eV wcd_Ti(k) = params_b(4, k) / ev ! erg -> eV wcd_Vz(k) = params_b(2, k) * rtor ! rad/s -> cm/s - wcd_dPhi0(k) = -Ercov(k) end do + ! Only update dPhi0 from Ercov during time evolution. + ! For SingleStep, preserve the input Er profile (from Er.dat). + if (trim(type_of_run) /= 'SingleStep') then + do k = 1, dim_r + wcd_dPhi0(k) = -Ercov(k) + end do + end if + end subroutine kim_update_profiles subroutine kim_get_wave_fields(i_mn) @@ -448,6 +648,27 @@ subroutine kim_get_collision_frequencies() end subroutine kim_get_collision_frequencies + ! --------------------------------------------------------------- + ! Helper: deallocate equilibrium_m arrays for re-initialization + ! --------------------------------------------------------------- + subroutine deallocate_equilibrium_arrays() + !! Deallocate arrays from equilibrium_m that calculate_equil + !! allocates with bare allocate() (no guard). + !! Separated from deallocate_plasma_derived to avoid circular + !! module dependency (equilibrium_m uses species_m). + implicit none + + if (allocated(kim_B0)) deallocate(kim_B0) + if (allocated(eq_B0z)) deallocate(eq_B0z) + if (allocated(eq_B0th)) deallocate(eq_B0th) + if (allocated(eq_hz)) deallocate(eq_hz) + if (allocated(eq_hth)) deallocate(eq_hth) + if (allocated(eq_equil_grid)) deallocate(eq_equil_grid) + if (allocated(eq_u)) deallocate(eq_u) + if (allocated(eq_dpress_prof)) deallocate(eq_dpress_prof) + + end subroutine deallocate_equilibrium_arrays + ! --------------------------------------------------------------- ! Helper: deallocate EBdat fields between mode solves ! --------------------------------------------------------------- @@ -476,6 +697,7 @@ subroutine deallocate_EBdat() if (allocated(EBdat%Phi_MA)) deallocate(EBdat%Phi_MA) if (allocated(EBdat%Phi_MA_ideal)) deallocate(EBdat%Phi_MA_ideal) if (allocated(EBdat%Phi_MA_asymptotic)) deallocate(EBdat%Phi_MA_asymptotic) + if (allocated(EBdat%jpar)) deallocate(EBdat%jpar) end subroutine deallocate_EBdat @@ -507,4 +729,172 @@ subroutine interp_complex_profile(n_old, r_old, z_old, & end subroutine interp_complex_profile + ! --------------------------------------------------------------- + ! Vacuum field loading and domain consistency + ! --------------------------------------------------------------- + + subroutine kim_load_vacuum_fields() + !! Extract vacuum Br from KiLCA vacuum solution (vac_cd_ptr) + !! onto the balance grid and store in kim_vac_Br. + use wave_code_data, only: dim_r, r, dim_mn, m_vals, n_vals, & + vac_cd_ptr, Br, Bz + + implicit none + + ! Extract vacuum Br for first mode (single-mode for now). + ! get_wave_fields_from_wave_code fills Br; we use Bz as dummy + ! for all other field components. + call get_wave_fields_from_wave_code(vac_cd_ptr(1), dim_r, r, & + m_vals(1), n_vals(1), Bz, Bz, Bz, Bz, Bz, Br, Bz, Bz, Bz, Bz) + kim_vac_Br = Br + + write(*,*) 'KIM adapter: loaded vacuum Br from KiLCA' + write(*,*) ' |Br_vac| at r_min = ', abs(kim_vac_Br(1)), ' at r=', r(1) + write(*,*) ' |Br_vac| at r_max = ', abs(kim_vac_Br(dim_r)), ' at r=', r(dim_r) + + call kim_check_domain_consistency() + + end subroutine kim_load_vacuum_fields + + subroutine kim_check_domain_consistency() + !! Set KIM Br boundary condition from KiLCA vacuum solution + !! at the actual KIM grid boundary (kim_r_boundary = xl_grid%xb(N)), + !! which may differ from r_plas due to non-equidistant grid construction. + !! Uses linear interpolation on the balance grid to get vacuum Br + !! at exactly kim_r_boundary. + use wave_code_data, only: dim_r, r, vac_cd_ptr + use setup_m, only: Br_boundary_re, Br_boundary_im + + implicit none + + integer :: i, i_lo, i_hi + real(8) :: w, r_bc + complex(8) :: Br_vac_at_bc + + ! Use actual KIM grid boundary for the boundary condition + r_bc = kim_r_boundary + + ! 1. Report domain info + write(*,*) 'KIM adapter: domain check' + write(*,*) ' KIM r_plas = ', kim_r_plas + write(*,*) ' KIM r_boundary = ', r_bc, ' (actual grid boundary)' + write(*,*) ' Balance grid max = ', r(dim_r) + + if (r_bc < r(dim_r)) then + write(*,*) ' Vacuum continuation beyond r_boundary to r=', r(dim_r), ' cm' + end if + + ! 2. Check vacuum solution is available + if (vac_cd_ptr(1) == 0) then + write(*,*) 'ERROR: KIM mode requires vacuum solution but vac_cd_ptr is not loaded.' + write(*,*) ' Set vac_path in balance_conf.nml and run KiLCA vacuum solver.' + stop 1 + end if + + ! 3. Set KIM Br boundary from vacuum Br at exactly r_bc + ! Find bracketing indices for linear interpolation + i_lo = 1 + do i = 1, dim_r - 1 + if (r(i) <= r_bc .and. r(i+1) >= r_bc) then + i_lo = i + exit + end if + end do + i_hi = i_lo + 1 + + ! Linear interpolation weight + if (abs(r(i_hi) - r(i_lo)) > 1.0d-30) then + w = (r_bc - r(i_lo)) / (r(i_hi) - r(i_lo)) + else + w = 0.0d0 + end if + + Br_vac_at_bc = (1.0d0 - w) * kim_vac_Br(i_lo) & + + w * kim_vac_Br(i_hi) + + write(*,*) ' Vacuum Br at r_boundary (linear interp):' + write(*,*) ' r_lo=', r(i_lo), ' r_hi=', r(i_hi), ' w=', w + write(*,*) ' |Br_vac(r_lo)| = ', abs(kim_vac_Br(i_lo)) + write(*,*) ' |Br_vac(r_hi)| = ', abs(kim_vac_Br(i_hi)) + write(*,*) ' |Br_vac(r_bc)| = ', abs(Br_vac_at_bc) + write(*,*) ' Re(Br) = ', real(Br_vac_at_bc) + write(*,*) ' Im(Br) = ', aimag(Br_vac_at_bc) + + ! Override KIM boundary condition with vacuum value + write(*,*) ' Old Br_boundary: Re=', Br_boundary_re, ' Im=', Br_boundary_im + Br_boundary_re = real(Br_vac_at_bc) + Br_boundary_im = aimag(Br_vac_at_bc) + write(*,*) ' New Br_boundary: Re=', Br_boundary_re, ' Im=', Br_boundary_im + + end subroutine kim_check_domain_consistency + + subroutine kim_get_current_densities(i_mn) + !! Copy per-mode stored parallel current density from + !! kim_jpar_modes into the wave_code_data J arrays for mode i_mn. + !! KIM provides total jpar; it goes into Jpi, rest is zeroed. + use wave_code_data, only: Jpi, Jpe, Jri, Jre, Jsi, Jse + + implicit none + + integer, intent(in) :: i_mn + + Jpi = kim_jpar_modes(:, i_mn) + Jpe = (0.0d0, 0.0d0) + Jri = (0.0d0, 0.0d0) + Jre = (0.0d0, 0.0d0) + Jsi = (0.0d0, 0.0d0) + Jse = (0.0d0, 0.0d0) + + end subroutine kim_get_current_densities + + subroutine apply_vacuum_continuation(i_mn, nrad, bal_r) + !! For grid points beyond kim_r_boundary (actual KIM grid extent), + !! replace KIM fields with vacuum values: E fields = 0, + !! Br = kim_vac_Br, wave vectors = 0. + implicit none + + integer, intent(in) :: i_mn, nrad + real(8), intent(in) :: bal_r(nrad) + integer :: i, i_last_kim + + ! Find last balance grid point inside KIM domain + ! Use kim_r_boundary (actual grid extent), not kim_r_plas + i_last_kim = nrad + do i = 1, nrad + if (bal_r(i) > kim_r_boundary) then + i_last_kim = i - 1 + exit + end if + end do + + ! Diagnostic: print transition values at the stitching boundary + if (i_mn == 1 .and. i_last_kim < nrad) then + write(*,*) ' Vacuum continuation transition (at r_boundary=', & + kim_r_boundary, '):' + write(*,*) ' Last KIM point: r=', bal_r(i_last_kim), & + ' |Br_KIM|=', abs(kim_Br_modes(i_last_kim, i_mn)) + write(*,*) ' First vac point: r=', bal_r(i_last_kim+1), & + ' |Br_vac|=', abs(kim_vac_Br(i_last_kim+1)) + write(*,*) ' Vac Br at last KIM r: |Br_vac|=', & + abs(kim_vac_Br(i_last_kim)) + write(*,*) ' Value jump: |Br_KIM - Br_vac| = ', & + abs(kim_Br_modes(i_last_kim, i_mn) - kim_vac_Br(i_last_kim)) + end if + + do i = i_last_kim + 1, nrad + ! No plasma response in vacuum: E fields = 0 + kim_Es_modes(i, i_mn) = (0.0d0, 0.0d0) + kim_Ep_modes(i, i_mn) = (0.0d0, 0.0d0) + kim_Er_modes(i, i_mn) = (0.0d0, 0.0d0) + kim_Et_modes(i, i_mn) = (0.0d0, 0.0d0) + kim_Ez_modes(i, i_mn) = (0.0d0, 0.0d0) + ! Total Br = vacuum Br beyond plasma + kim_Br_modes(i, i_mn) = kim_vac_Br(i) + ! Wave vectors not physical in vacuum + kim_kp_modes(i, i_mn) = 0.0d0 + kim_ks_modes(i, i_mn) = 0.0d0 + end do + + end subroutine apply_vacuum_continuation + end module kim_wave_code_adapter_m diff --git a/QL-Balance/src/base/paramscan.f90 b/QL-Balance/src/base/paramscan.f90 index 5d89c7bb..ba4f80f5 100644 --- a/QL-Balance/src/base/paramscan.f90 +++ b/QL-Balance/src/base/paramscan.f90 @@ -77,7 +77,8 @@ subroutine initParameterScan(this) subroutine runParameterScan(this) - use transp_coeffs_mod, only: rescale_transp_coeffs_by_ant_fac + use transp_coeffs_mod, only: rescale_transp_coeffs_by_ant_fac, & + compute_antenna_factor_from_Ipar use plasma_parameters, only: alloc_hold_parameters implicit none @@ -105,6 +106,7 @@ subroutine runParameterScan(this) print *, "Before get dql" call get_dql print *, "Before rescale transp" + call compute_antenna_factor_from_Ipar call rescale_transp_coeffs_by_ant_fac print *, "Before interpolate Br Dql" call interpolate_Br_Dql_at_res_parscan diff --git a/QL-Balance/src/base/read_config.f90 b/QL-Balance/src/base/read_config.f90 index c003b6d7..52c671a2 100644 --- a/QL-Balance/src/base/read_config.f90 +++ b/QL-Balance/src/base/read_config.f90 @@ -3,13 +3,14 @@ subroutine read_config use control_mod, only: eps, paramscan, diagnostics_output, suppression_mode, debug_mode, & readfromtimestep, temperature_limit, gyro_current_study, & misalign_diffusion, equil_path, ihdf5IO, wave_code, & - kim_config_path + kim_config_path, kim_profiles_from_balance, & + kim_n_modes, kim_m_list, kim_n_list use grid_mod, only: rmin, rmax, npoimin, gg_factor, gg_width, gg_r_res, iboutype, rb_cut_in, & re_cut_in, rb_cut_out, re_cut_out use h5mod, only: path2inp, path2out, path2time use paramscan_mod, only: viscosity_factor use time_evolution - use wave_code_data, only: flre_path, vac_path, antenna_factor + use wave_code_data, only: flre_path, vac_path, antenna_factor, I_par_toroidal implicit none @@ -25,7 +26,9 @@ subroutine read_config suppression_mode, debug_mode, readfromtimestep, path2time, ramp_up_mode, t_max_ramp_up, & temperature_limit, antenna_max_stopping, gyro_current_study, viscosity_factor, & misalign_diffusion, equil_path, ihdf5IO, type_of_run, wave_code, & - set_constant_time_step, constant_time_step, urelax, kim_config_path + set_constant_time_step, constant_time_step, urelax, kim_config_path, & + kim_profiles_from_balance, kim_n_modes, kim_m_list, kim_n_list, & + I_par_toroidal ! read the parameters from namelist file open (newunit=u, file=config_file, status="old", action="read", iostat=ios) @@ -55,6 +58,7 @@ subroutine read_config write (*, "(A,ES15.8)") " tmax_factor = ", tmax_factor write (*, "(A,ES15.8)") " timstep_min = ", timstep_min write (*, "(A,ES15.8)") " antenna_factor = ", antenna_factor + write (*, "(A,ES15.8)") " I_par_toroidal = ", I_par_toroidal write (*, "(A,I0)") " iboutype = ", iboutype write (*, "(A,ES15.8)") " eps = ", eps write (*, "(A,ES15.8)") " dperp = ", dperp @@ -82,6 +86,12 @@ subroutine read_config write (*, "(A,ES15.8,A)") " constant_time_step = ", constant_time_step, " s" write (*, "(A,ES15.8)") " urelax = ", urelax write (*, "(A,A)") " kim_config_path = ", trim(adjustl(kim_config_path)) + write (*, "(A,L0)") " kim_profiles_from_balance = ", kim_profiles_from_balance + write (*, "(A,I0)") " kim_n_modes = ", kim_n_modes + if (kim_n_modes > 0) then + write (*, "(A,100I5)") " kim_m_list = ", kim_m_list(1:kim_n_modes) + write (*, "(A,100I5)") " kim_n_list = ", kim_n_list(1:kim_n_modes) + end if write (*, *) "" write (*, *) "=================================================================================" end subroutine read_config diff --git a/QL-Balance/src/base/time_evolution.f90 b/QL-Balance/src/base/time_evolution.f90 index 05baa1c9..9c1afe68 100644 --- a/QL-Balance/src/base/time_evolution.f90 +++ b/QL-Balance/src/base/time_evolution.f90 @@ -79,7 +79,8 @@ module time_evolution subroutine initTimeEvolution(this) use recstep_mod, only: tol - use transp_coeffs_mod, only: rescale_transp_coeffs_by_ant_fac + use transp_coeffs_mod, only: rescale_transp_coeffs_by_ant_fac, & + compute_antenna_factor_from_Ipar use grid_mod, only: mwind, rmax, rmin, set_boundary_condition, npoib, rb use baseparam_mod, only: dperp, tol_max use QLbalance_diag, only: write_diag, write_diag_b @@ -152,6 +153,7 @@ subroutine initTimeEvolution(this) call alloc_Br_Dqle_for_timeevol call get_dql + call compute_antenna_factor_from_Ipar call rescale_transp_coeffs_by_ant_fac call hold_prev_transp_coeffs @@ -174,7 +176,8 @@ subroutine doStep(this) use recstep_mod, only: timstep_arr use recstep_mod, only: tol use restart_mod, only: redostep - use transp_coeffs_mod, only: rescale_transp_coeffs_by_ant_fac + use transp_coeffs_mod, only: rescale_transp_coeffs_by_ant_fac, & + compute_antenna_factor_from_Ipar implicit none @@ -186,6 +189,7 @@ subroutine doStep(this) redostep = .false. call get_dql + call compute_antenna_factor_from_Ipar call rescale_transp_coeffs_by_ant_fac call interp_Br_Dql_at_resonance_timeevol call determine_Dql_diagnostic diff --git a/QL-Balance/src/base/transp_coeffs_mod.f90 b/QL-Balance/src/base/transp_coeffs_mod.f90 index 2584dc0d..e8d2ebe5 100644 --- a/QL-Balance/src/base/transp_coeffs_mod.f90 +++ b/QL-Balance/src/base/transp_coeffs_mod.f90 @@ -23,4 +23,53 @@ subroutine rescale_transp_coeffs_by_ant_fac end subroutine + subroutine compute_antenna_factor_from_Ipar + ! + ! Compute antenna_factor = (I_par_toroidal / I_par_cylindrical)^2 + ! from the GPEC toroidal shielding current and the wave-code + ! cylindrical shielding current integrated at runtime. + ! + ! antenna_factor is the squared ratio of currents (= squared field ratio). + ! It is applied once to D_ql in rescale_transp_coeffs_by_ant_fac, + ! consistent with the legacy namelist convention. + ! + ! Unit conversion: + ! I_par_toroidal is in CGS with c=1 (Gaussian-like, from GPEC) + ! Ipar from integrate_parallel_current is int(J*r*dr) in full CGS + ! without the 2*pi factor, i.e. physical current = 2*pi*Ipar + ! To convert I_par_toroidal to full CGS: multiply by c + ! So: antenna_factor = (I_par_toroidal * c / (2*pi * |Ipar|))^2 + ! + use wave_code_data, only: antenna_factor, I_par_toroidal + use grid_mod, only: Ipar + + implicit none + + double precision, parameter :: c_light = 2.99792458d10 ! speed of light [cm/s] + double precision, parameter :: twopi = 6.283185307179586d0 + double precision :: I_tor_cgs, I_cyl_cgs + + if (I_par_toroidal <= 0.0d0) return ! legacy mode: use namelist antenna_factor + + if (abs(Ipar) < 1.0d-30) then + write(*,*) 'WARNING: I_par from wave code is ~0; keeping namelist antenna_factor' + return + end if + + ! Convert GPEC current (c=1 CGS) to full CGS + I_tor_cgs = I_par_toroidal * c_light + ! Physical cylindrical current = 2*pi * int(J*r*dr) + I_cyl_cgs = twopi * abs(Ipar) + + antenna_factor = (I_tor_cgs / I_cyl_cgs)**2 + + write(*,*) 'Antenna factor computed from I_par_toroidal:' + write(*,*) ' I_par_toroidal (c=1 CGS) = ', I_par_toroidal + write(*,*) ' I_par_toroidal (full CGS) = ', I_tor_cgs + write(*,*) ' I_par_cylindrical (full CGS) = ', I_cyl_cgs + write(*,*) ' I_par_cylindrical (c=1 CGS) = ', I_cyl_cgs / c_light + write(*,*) ' antenna_factor = ', antenna_factor + + end subroutine + end module diff --git a/QL-Balance/src/base/wave_code_data_64bit.f90 b/QL-Balance/src/base/wave_code_data_64bit.f90 index 4d878b5c..df64108a 100644 --- a/QL-Balance/src/base/wave_code_data_64bit.f90 +++ b/QL-Balance/src/base/wave_code_data_64bit.f90 @@ -33,8 +33,10 @@ subroutine initialize_wave_code_interface(nrad, r_grid) use wave_code_data use h5mod - use control_mod, only: readfromtimestep, debug_mode, ihdf5IO, wave_code - use kim_wave_code_adapter_m, only: kim_initialize + use control_mod, only: readfromtimestep, debug_mode, ihdf5IO, wave_code, & + kim_profiles_from_balance, kim_n_modes, kim_m_list, kim_n_list, & + type_of_run + use kim_wave_code_adapter_m, only: kim_initialize, kim_load_vacuum_fields implicit none @@ -77,14 +79,47 @@ subroutine initialize_wave_code_interface(nrad, r_grid) call get_collision_frequences_from_wave_code(vac_cd_ptr(1), dim_r, r, nui, nue); vac_call_ind = vac_call_ind + 1; - ! Enable interface mode: KiLCA gets profiles directly from QL-Balance for all run types - ! that use initialize_wave_code_interface. This is set AFTER vacuum field initialization - ! to preserve original file-based behavior there. - call set_flag_for_profiles_in_background_input_file(trim(flre_path), -1) + ! Enable interface mode: KiLCA gets profiles directly from QL-Balance. + ! For SingleStep, profiles don't evolve so file-based reading (calc_back=1) + ! is both sufficient and numerically more robust (avoids interpolation chain). + if (trim(type_of_run) /= 'SingleStep') then + call set_flag_for_profiles_in_background_input_file(trim(flre_path), -1) + end if case ('KIM') + if (kim_profiles_from_balance) then + ! Set modes from balance namelist (no modes.in needed) + dim_mn = kim_n_modes + allocate(m_vals(dim_mn), n_vals(dim_mn)) + m_vals(1:dim_mn) = kim_m_list(1:dim_mn) + n_vals(1:dim_mn) = kim_n_list(1:dim_mn) + + call allocate_wave_code_data(nrad, r_grid) + + ! Read QL-Balance profiles (same as KiLCA path) + if (ihdf5IO .eq. 1) then + if (readfromtimestep .eq. 0) then + call read_background_profiles_h5 + else + call read_background_profiles_h5_timeevol(readfromtimestep) + end if + else + call read_background_profiles(profile_path) + end if + call interp_background_profiles() + end if call kim_initialize(nrad, r_grid) + ! Load vacuum solution from KiLCA vacuum run + do k = 1, dim_mn + call clear_wave_code_data(vac_cd_ptr(k)) + call calc_wave_code_data_for_mode(vac_cd_ptr(k), & + trim(vac_path), len(trim(vac_path)), m_vals(k), n_vals(k)) + end do + + ! Fill kim_vac_Br from vacuum solution and check domain consistency + call kim_load_vacuum_fields() + case default write(*,*) "ERROR: Unknown wave_code: ", trim(wave_code) stop 1 diff --git a/QL-Balance/src/base/wave_code_data_mod.f90 b/QL-Balance/src/base/wave_code_data_mod.f90 index 77f06241..a30bfc47 100644 --- a/QL-Balance/src/base/wave_code_data_mod.f90 +++ b/QL-Balance/src/base/wave_code_data_mod.f90 @@ -69,6 +69,7 @@ module wave_code_data real(8), allocatable, dimension(:) :: rVz, iVz; !cm/c real(8), allocatable, dimension(:) :: rep, idPhi0; !electric potential, Gauss units double precision :: antenna_factor + double precision :: I_par_toroidal = 0.0d0 character(1024) :: path2profs = './profiles/'; ! path to background profiles end module diff --git a/QL-Balance/src/base/writeData.f90 b/QL-Balance/src/base/writeData.f90 index aa14ec21..4c013734 100644 --- a/QL-Balance/src/base/writeData.f90 +++ b/QL-Balance/src/base/writeData.f90 @@ -32,11 +32,21 @@ subroutine write_fields_currs_transp_coefs_to_h5 real(Br), lbound(Br), ubound(Br)) CALL h5_add_double_1(h5_id, trim(tempch)//"Br_Im", & dimag(Br), lbound(Br), ubound(Br)) + CALL h5_add_double_1(h5_id, trim(tempch)//"Es_abs", & + abs(Es), lbound(Es), ubound(Es)) CALL h5_add_double_1(h5_id, trim(tempch)//"Jpe_abs", & abs(Jpe), lbound(Jpe), ubound(Jpe)) + CALL h5_add_double_1(h5_id, trim(tempch)//"Jpe_Re", & + real(Jpe), lbound(Jpe), ubound(Jpe)) + CALL h5_add_double_1(h5_id, trim(tempch)//"Jpe_Im", & + dimag(Jpe), lbound(Jpe), ubound(Jpe)) CALL h5_add_double_1(h5_id, trim(tempch)//"Jpi_abs", & abs(Jpi), lbound(Jpi), ubound(Jpi)) + CALL h5_add_double_1(h5_id, trim(tempch)//"Jpi_Re", & + real(Jpi), lbound(Jpi), ubound(Jpi)) + CALL h5_add_double_1(h5_id, trim(tempch)//"Jpi_Im", & + dimag(Jpi), lbound(Jpi), ubound(Jpi)) CALL h5_add_double_1(h5_id, trim(tempch)//"dqle22", & dqle22, lbound(dqle22), ubound(dqle22)) From d31b69a82f33630e8b82dafe6304fee9a1e0c52a Mon Sep 17 00:00:00 2001 From: Markus Markl Date: Tue, 24 Mar 2026 09:58:56 +0100 Subject: [PATCH 15/20] refactor(python): modernize file operations and add KIM balance mode Replace os.system calls with os.makedirs, os.remove, and os.symlink in KiLCA_interface for safer file handling. Extract vacuum KiLCA run into prepare_vacuum_kilca method and add prepare_balance_kim for KIM-mode workflow (profiles + vacuum fields without full KiLCA plasma response). Fix balance_input_h5 ubounds calculation (len(data) instead of len(data[0])). Set vgalsys=-1e8 in background data for consistency. Co-Authored-By: Claude Opus 4.6 (1M context) --- python/KiLCA_interface/KiLCA_interface.py | 78 +++++++++---------- python/balance_interface/balance_input_h5.py | 2 +- python/balance_interface/balance_interface.py | 47 +++++++++-- 3 files changed, 77 insertions(+), 50 deletions(-) diff --git a/python/KiLCA_interface/KiLCA_interface.py b/python/KiLCA_interface/KiLCA_interface.py index fa15a7a2..73f6410d 100644 --- a/python/KiLCA_interface/KiLCA_interface.py +++ b/python/KiLCA_interface/KiLCA_interface.py @@ -1,4 +1,5 @@ import copy +import glob import inspect import os import subprocess as sp @@ -338,6 +339,13 @@ def set_MASTU(self, nmodes: int = 1, I0: float = 6.0e12, delta_r_antenna: float m = [self.run_type, "vacuum", "vacuum"] self.set_zones(r, b, m) + @staticmethod + def _force_symlink(source, link_name): + """Create a symlink, removing any existing file/link at the destination.""" + if os.path.islink(link_name) or os.path.exists(link_name): + os.unlink(link_name) + os.symlink(source, link_name) + def write(self): """ Description: @@ -347,20 +355,26 @@ def write(self): self.check_profile_consistency() # create paths if they are not there - os.system("mkdir -p " + self.path) - os.system("mkdir -p " + self.path_of_profiles) - os.system("mkdir -p " + self.path_of_run) + os.makedirs(self.path, exist_ok=True) + os.makedirs(self.path_of_profiles, exist_ok=True) + os.makedirs(self.path_of_run, exist_ok=True) # delete all existing input files - os.system("rm -f " + self.path + "background.in") - os.system("rm -f " + self.path + "eigmode.in") - os.system("rm -f " + self.path + "modes.in") - os.system("rm -f " + self.path + "output.in") - os.system("rm -f " + self.path_of_run + "antenna.in") - os.system("rm -f " + self.path_of_run + "zone*.in") + for f in ["background.in", "eigmode.in", "modes.in", "output.in"]: + p = os.path.join(self.path, f) + if os.path.exists(p): + os.remove(p) + antenna_path = os.path.join(self.path_of_run, "antenna.in") + if os.path.exists(antenna_path): + os.remove(antenna_path) + for p in glob.glob(os.path.join(self.path_of_run, "zone*.in")): + os.remove(p) # create sym link to exe - os.system("ln -sf " + self.EXEC_PATH + " " + self.path_of_run + "run_local") + run_local_path = os.path.join(self.path_of_run, "run_local") + if os.path.islink(run_local_path) or os.path.exists(run_local_path): + os.unlink(run_local_path) + os.symlink(self.EXEC_PATH, run_local_path) # create symbolic link to profiles directory # Safety check: path_of_run must be set and not be empty or just 'profiles' path_of_run_stripped = self.path_of_run.rstrip("/") @@ -413,47 +427,27 @@ def write(self): self.antenna.write(self.BLUE_PATH + self.antenna.BLUEPRINT, self.path_of_run) self.background.write(self.BLUE_PATH + self.background.BLUEPRINT, self.path) - os.system( - "ln -sf " - + os.path.abspath(self.path) - + "/" - + self.background.BLUEPRINT - + " " - + self.path_of_run - + self.background.BLUEPRINT + self._force_symlink( + os.path.join(os.path.abspath(self.path), self.background.BLUEPRINT), + os.path.join(self.path_of_run, self.background.BLUEPRINT), ) self.eigmode.write(self.BLUE_PATH + self.eigmode.BLUEPRINT, self.path) - os.system( - "ln -sf " - + os.path.abspath(self.path) - + "/" - + self.eigmode.BLUEPRINT - + " " - + self.path_of_run - + self.eigmode.BLUEPRINT + self._force_symlink( + os.path.join(os.path.abspath(self.path), self.eigmode.BLUEPRINT), + os.path.join(self.path_of_run, self.eigmode.BLUEPRINT), ) self.output.write(self.BLUE_PATH + self.output.BLUEPRINT, self.path) - os.system( - "ln -sf " - + os.path.abspath(self.path) - + "/" - + self.output.BLUEPRINT - + " " - + self.path_of_run - + self.output.BLUEPRINT + self._force_symlink( + os.path.join(os.path.abspath(self.path), self.output.BLUEPRINT), + os.path.join(self.path_of_run, self.output.BLUEPRINT), ) self.modes.write(self.path) - os.system( - "ln -sf " - + os.path.abspath(self.path) - + "/" - + self.modes.BLUEPRINT - + " " - + self.path_of_run - + self.modes.BLUEPRINT + self._force_symlink( + os.path.join(os.path.abspath(self.path), self.modes.BLUEPRINT), + os.path.join(self.path_of_run, self.modes.BLUEPRINT), ) for k in range(0, len(self.zones)): diff --git a/python/balance_interface/balance_input_h5.py b/python/balance_interface/balance_input_h5.py index df8877f9..d8c1812e 100644 --- a/python/balance_interface/balance_input_h5.py +++ b/python/balance_interface/balance_input_h5.py @@ -62,4 +62,4 @@ def write_with_bound_info(self, file, dataset, data): def write_fac_with_bound_info(self, file, dataset, data): ds = file.create_dataset(dataset, data=data) ds.attrs["lbounds"] = [1] - ds.attrs["ubounds"] = [len(data[0])] + ds.attrs["ubounds"] = [len(data)] diff --git a/python/balance_interface/balance_interface.py b/python/balance_interface/balance_interface.py index 623a6d52..a4081cd1 100644 --- a/python/balance_interface/balance_interface.py +++ b/python/balance_interface/balance_interface.py @@ -196,6 +196,7 @@ def prepare_KiLCA(self, Btor, a_minor): self.kil_flre = KiLCA_interface(self.shot, self.time, self.run_path, "flre", self.machine) self.kil_flre.set_machine(delta_r_antenna=self.delta_r_antenna, a_minor=a_minor) self.kil_flre.background.data["Btor"] = Btor + self.kil_flre.background.data["vgalsys"] = -1e8 self.kil_flre.set_modes(self.m_mode, self.n_mode) self.kil_flre.antenna.data["flab"] = [1.0, 0.0] @@ -204,13 +205,7 @@ def prepare_KiLCA(self, Btor, a_minor): self.kil_flre_post = KiLCA_postprocessor(self.kil_flre) self.I_KiLCA = self.get_KiLCA_current() - kil_vac = KiLCA_interface(self.shot, self.time, self.run_path, "vacuum", self.machine) - kil_vac.background.data["Btor"] = Btor - kil_vac.set_machine(delta_r_antenna=self.delta_r_antenna, a_minor=a_minor) - kil_vac.set_modes(self.m_mode, self.n_mode) - kil_vac.antenna.data["flab"] = [1.0, 0.0] - kil_vac.write() - kil_vac.run() + self.prepare_vacuum_kilca(Btor, a_minor) if self.debug: print("D: Finished KiLCA preperation") @@ -279,6 +274,44 @@ def write_KiLCA_data_to_input_h5(self): grp.create_dataset("layer_width", data=[self.layer_width]) h5f.close() + def prepare_balance_kim(self, Btor, a_minor): + """Prepare QL-Balance for KIM mode with vacuum solution. + + When wave_code='KIM', the Fortran QL-Balance code calls KIM internally. + We prepare the HDF5 input (profiles + Da), run the KiLCA vacuum solver + (needed for vacuum fields beyond the plasma), and set up the output file. + """ + if self.debug: + print("D: Prepare balance for KIM mode") + self.check_if_factors_set() + self.prepare_input_h5_without_kilca() + self.prepare_vacuum_kilca(Btor, a_minor) + self.prepare_output_h5() + self.link_executable() + + def prepare_vacuum_kilca(self, Btor, a_minor): + """Run KiLCA vacuum solver only (no plasma response).""" + if self.debug: + print("D: Running KiLCA vacuum solver") + kil_vac = KiLCA_interface(self.shot, self.time, self.run_path, "vacuum", self.machine) + kil_vac.set_machine(delta_r_antenna=self.delta_r_antenna, a_minor=a_minor) + kil_vac.background.data["Btor"] = Btor + kil_vac.background.data["vgalsys"] = -1e8 + kil_vac.set_modes(self.m_mode, self.n_mode) + kil_vac.antenna.data["flab"] = [1.0, 0.0] + kil_vac.write() + kil_vac.run() + if self.debug: + print("D: KiLCA vacuum solver finished") + + def prepare_input_h5_without_kilca(self): + """Prepare input HDF5 with profiles and Da only (no KiLCA data).""" + self.input_h5 = Balance_Input_h5( + self.input_h5_file, os.path.join(self.run_path, "profiles/") + ) + self.input_h5.get_required_data() + self.input_h5.write_data_to_h5(self.input_h5_file, self.facs) + def run_balance(self, suppress_console_output=True): """Run the balance code.""" print("") From 17a31a3e1185c944b2bee17e0489beffdb026b04 Mon Sep 17 00:00:00 2001 From: Markus Markl Date: Tue, 24 Mar 2026 10:36:57 +0100 Subject: [PATCH 16/20] fix(KIM/QL-Balance): address review findings for KIM integration - Recompute equilibrium (kp, ks, om_E) before each mode 2+ solve so multi-mode runs use correct mode-dependent wave vectors - Make vacuum Br per-mode: kim_vac_Br is now (dim_r, dim_mn), loaded from KiLCA vacuum solution for each (m,n) independently - Compute per-species parallel current via electron/ion kernels: add jpar_e, jpar_i to EBdat_t, fill from Kllp_e/Kllp_i in EM solver, assign correctly to Jpe/Jpi in the adapter - Initialize kim_vac_Br with IEEE NaN instead of (1,0) to catch uninitialized usage via NaN propagation - Add case default error stop to all wave_code select blocks - Fix KIM test linking: add KiLCA library to test executables after KiLCA was removed from KIM_lib PUBLIC linkage - Delete unused kim_wave_code_adapter_stub.f90 Co-Authored-By: Claude Opus 4.6 (1M context) --- .../electromagnetic_solver.f90 | 18 ++- KIM/src/electrostatic_poisson/fields_mod.f90 | 4 +- KIM/tests/CMakeLists.txt | 7 +- QL-Balance/src/base/get_dql.f90 | 14 ++- QL-Balance/src/base/kim_wave_code_adapter.f90 | 108 ++++++++++++------ 5 files changed, 110 insertions(+), 41 deletions(-) diff --git a/KIM/src/electromagnetic/electromagnetic_solver.f90 b/KIM/src/electromagnetic/electromagnetic_solver.f90 index 4a6092b9..95cef097 100644 --- a/KIM/src/electromagnetic/electromagnetic_solver.f90 +++ b/KIM/src/electromagnetic/electromagnetic_solver.f90 @@ -73,8 +73,9 @@ subroutine run_electromagnetic(this) real(dp), allocatable :: ks_xl(:) complex(dp), allocatable :: alpha(:) - complex(dp), allocatable :: jpar(:), rho(:) + complex(dp), allocatable :: jpar(:), jpar_sp(:), rho(:) complex(dp) :: Br_boundary + integer :: sp integer :: N, i, j real(dp) :: kz, hL, hR @@ -272,7 +273,7 @@ subroutine run_electromagnetic(this) call write_complex_profile_abs(xl_grid%xb, EBdat%Br, N, "/fields/Br_selfconsistent", & 'Self-consistent radial magnetic field perturbation Br', 'G') - ! Current density + ! Current density: total call calculate_current_density(jpar, EBdat%Phi, EBdat%Br, & kernel_j_phi_llp%Kllp, kernel_j_B_llp%Kllp) allocate(EBdat%jpar(N)) @@ -280,6 +281,19 @@ subroutine run_electromagnetic(this) call write_complex_profile_abs(xl_grid%xb, jpar, N, "/fields/jpar", & 'Parallel current density from self-consistent solve', 'statA/cm^2') + ! Current density: electron + call calculate_current_density(EBdat%jpar_e, EBdat%Phi, EBdat%Br, & + kernel_j_phi_llp%Kllp_e, kernel_j_B_llp%Kllp_e) + + ! Current density: ion (sum over all ion species) + allocate(EBdat%jpar_i(N)) + EBdat%jpar_i = (0.0d0, 0.0d0) + do sp = 1, plasma%n_species - 1 + call calculate_current_density(jpar_sp, EBdat%Phi, EBdat%Br, & + kernel_j_phi_llp%Kllp_i(:,:,sp), kernel_j_B_llp%Kllp_i(:,:,sp)) + EBdat%jpar_i = EBdat%jpar_i + jpar_sp + end do + ! Electric field postprocessing call postprocess_electric_field(EBdat) diff --git a/KIM/src/electrostatic_poisson/fields_mod.f90 b/KIM/src/electrostatic_poisson/fields_mod.f90 index 20aee696..f6d94a75 100644 --- a/KIM/src/electrostatic_poisson/fields_mod.f90 +++ b/KIM/src/electrostatic_poisson/fields_mod.f90 @@ -16,7 +16,9 @@ module fields_m complex(dp), allocatable :: Ez(:) complex(dp), allocatable :: Es(:) ! E "senkrecht", i.e. perpendicular to radial and parallel direction complex(dp), allocatable :: Ep(:) ! parallel to the equilibrium magnetic field - complex(dp), allocatable :: jpar(:) ! parallel current density + complex(dp), allocatable :: jpar(:) ! total parallel current density + complex(dp), allocatable :: jpar_e(:) ! electron parallel current density + complex(dp), allocatable :: jpar_i(:) ! ion parallel current density (sum over ion species) complex(dp), allocatable :: Phi(:) complex(dp), allocatable :: Phi_e(:) complex(dp), allocatable :: Phi_i(:) diff --git a/KIM/tests/CMakeLists.txt b/KIM/tests/CMakeLists.txt index 9a3ede2a..fb744268 100644 --- a/KIM/tests/CMakeLists.txt +++ b/KIM/tests/CMakeLists.txt @@ -1,9 +1,12 @@ +# KiLCA library path (needed for tests that use KIM_lib symbols depending on KiLCA math) +set(KiLCA_lib_path ${CMAKE_SOURCE_DIR}/build/install/lib/libKiLCA_Lib_V_2.4.2_MDNO_FPGEN_POLYNOMIAL_Release_64bit.a) + # Optional: build QUADPACK/RKF45 integration comparison test add_executable(test_integration_methods ${CMAKE_SOURCE_DIR}/KIM/tests/test_integration_methods.f90) set_target_properties(test_integration_methods PROPERTIES OUTPUT_NAME test_integration_methods.x RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/tests/") -target_link_libraries(test_integration_methods KIM_lib lapack cerf slatec ddeabm) +target_link_libraries(test_integration_methods KIM_lib "${KiLCA_lib_path}" lapack cerf slatec ddeabm) add_test(NAME test_integration_methods COMMAND ${CMAKE_BINARY_DIR}/tests/test_integration_methods.x) @@ -45,6 +48,6 @@ add_executable(test_ampere_matrices ${CMAKE_SOURCE_DIR}/KIM/tests/test_ampere_ma set_target_properties(test_ampere_matrices PROPERTIES OUTPUT_NAME test_ampere_matrices.x RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/tests/") -target_link_libraries(test_ampere_matrices KIM_lib lapack cerf slatec ddeabm) +target_link_libraries(test_ampere_matrices KIM_lib "${KiLCA_lib_path}" lapack cerf slatec ddeabm) add_test(NAME test_ampere_matrices COMMAND ${CMAKE_BINARY_DIR}/tests/test_ampere_matrices.x) diff --git a/QL-Balance/src/base/get_dql.f90 b/QL-Balance/src/base/get_dql.f90 index 112e3fb2..ba88845a 100644 --- a/QL-Balance/src/base/get_dql.f90 +++ b/QL-Balance/src/base/get_dql.f90 @@ -129,6 +129,8 @@ subroutine get_dql if (irf .eq. 1) call kim_update_profiles() if (irf .eq. 1) call kim_run_for_all_modes() ! B0 and collision freqs already set in kim_initialize + case default + error stop "Unknown wave_code: "//trim(wave_code) end select ! nu_e=15.4d-6*params_b(1,:)/sqrt(params_b(3,:)/ev)**3 & @@ -161,6 +163,8 @@ subroutine get_dql case ('KIM') call kim_get_wave_vectors(i_mn) call kim_get_wave_fields(i_mn) + case default + error stop "Unknown wave_code: "//trim(wave_code) end select om_E = ks * c * dPhi0 / B0 vT_e = sqrt(params_b(3, :)/e_mass) @@ -218,6 +222,8 @@ subroutine get_dql m_vals(i_mn), n_vals(i_mn), Er, Es, Ep, Et, Ez, Br, Bs, Bp, Bt, Bz) case ('KIM') call kim_get_wave_fields(i_mn) ! already loaded, but refresh + case default + error stop "Unknown wave_code: "//trim(wave_code) end select ! caluclate part of perpendicular electric field perturbation that comes from @@ -250,7 +256,9 @@ subroutine get_dql call get_wave_fields_from_wave_code(vac_cd_ptr(i_mn), dim_r, r, & m_vals(i_mn), n_vals(i_mn), Bz, Bz, Bz, Bz, Bz, Br, Bz, Bz, Bz, Bz) case ('KIM') - Br = kim_vac_Br + Br = kim_vac_Br(:, i_mn) + case default + error stop "Unknown wave_code: "//trim(wave_code) end select formfactor = (1.d0, 0.d0)/Br @@ -283,6 +291,8 @@ subroutine get_dql m_vals(i_mn), n_vals(i_mn), Bz, Bz, Bz, Bz, Bz, Br, Bz, Bz, Bz, Bz) case ('KIM') Br = kim_Br_modes(:, i_mn) + case default + error stop "Unknown wave_code: "//trim(wave_code) end select formfactor = Br * formfactor @@ -315,6 +325,8 @@ subroutine get_dql m_vals(i_mn), n_vals(i_mn), Jri, Jsi, Jpi, Jre, Jse, Jpe) case ('KIM') call kim_get_current_densities(i_mn) + case default + error stop "Unknown wave_code: "//trim(wave_code) end select call integrate_parallel_current(dim_r, r, Jpe, Jpi, r_resonant(i_mn), Ipar) diff --git a/QL-Balance/src/base/kim_wave_code_adapter.f90 b/QL-Balance/src/base/kim_wave_code_adapter.f90 index c4330ed1..67e5df6f 100644 --- a/QL-Balance/src/base/kim_wave_code_adapter.f90 +++ b/QL-Balance/src/base/kim_wave_code_adapter.f90 @@ -8,11 +8,13 @@ module kim_wave_code_adapter_m !! All KIM symbols are renamed on import to avoid name clashes. use species_m, only: kim_plasma => plasma, & - set_profiles_from_arrays, deallocate_plasma_derived + set_profiles_from_arrays, deallocate_plasma_derived, & + set_plasma_quantities use equilibrium_m, only: kim_B0 => B0, & eq_B0z => B0z, eq_B0th => B0th, & eq_hz => hz, eq_hth => hth, eq_equil_grid => equil_grid, & - eq_u => u, eq_dpress_prof => dpress_prof + eq_u => u, eq_dpress_prof => dpress_prof, & + calculate_equil, interpolate_equil use config_m, only: kim_type_of_run => type_of_run, & kim_nml_config_path => nml_config_path, & kim_profiles_in_memory => profiles_in_memory, & @@ -47,8 +49,8 @@ module kim_wave_code_adapter_m !! Used for vacuum continuation cutoff and Br boundary condition. real(8) :: kim_r_boundary = 0.0d0 - !! Vacuum Br placeholder (filled properly in Task 7) - complex(8), allocatable, public :: kim_vac_Br(:) + !! Per-mode vacuum Br on QL-Balance grid (dim_r, dim_mn) + complex(8), allocatable, public :: kim_vac_Br(:,:) !! Per-mode stored field results (nrad, dim_mn) !! Filled by kim_run_for_all_modes, read by kim_get_wave_fields @@ -65,8 +67,10 @@ module kim_wave_code_adapter_m real(8), allocatable :: kim_kp_modes(:,:) real(8), allocatable :: kim_ks_modes(:,:) - !! Per-mode stored current density (nrad, dim_mn) + !! Per-mode stored current densities (nrad, dim_mn) complex(8), allocatable :: kim_jpar_modes(:,:) + complex(8), allocatable :: kim_jpar_e_modes(:,:) + complex(8), allocatable :: kim_jpar_i_modes(:,:) contains @@ -269,10 +273,15 @@ subroutine kim_initialize(nrad, r_grid) deallocate(kim_r, work_old, work_new) end if - ! Vacuum Br placeholder - if (allocated(kim_vac_Br)) deallocate(kim_vac_Br) - allocate(kim_vac_Br(nrad)) - kim_vac_Br = (1.0d0, 0.0d0) + ! Vacuum Br placeholder (NaN until kim_load_vacuum_fields fills it) + block + use ieee_arithmetic, only: ieee_value, ieee_quiet_nan + real(8) :: nan_val + nan_val = ieee_value(1.0d0, ieee_quiet_nan) + if (allocated(kim_vac_Br)) deallocate(kim_vac_Br) + allocate(kim_vac_Br(nrad, dim_mn)) + kim_vac_Br = cmplx(nan_val, nan_val, kind=8) + end block write(*, *) "KIM adapter: initialization complete" write(*, *) " Balance grid points: ", nrad @@ -389,6 +398,8 @@ subroutine kim_run_for_all_modes() if (allocated(kim_kp_modes)) deallocate(kim_kp_modes) if (allocated(kim_ks_modes)) deallocate(kim_ks_modes) if (allocated(kim_jpar_modes)) deallocate(kim_jpar_modes) + if (allocated(kim_jpar_e_modes)) deallocate(kim_jpar_e_modes) + if (allocated(kim_jpar_i_modes)) deallocate(kim_jpar_i_modes) allocate(kim_Es_modes(dim_r, dim_mn)) allocate(kim_Ep_modes(dim_r, dim_mn)) @@ -399,6 +410,8 @@ subroutine kim_run_for_all_modes() allocate(kim_kp_modes(dim_r, dim_mn)) allocate(kim_ks_modes(dim_r, dim_mn)) allocate(kim_jpar_modes(dim_r, dim_mn)) + allocate(kim_jpar_e_modes(dim_r, dim_mn)) + allocate(kim_jpar_i_modes(dim_r, dim_mn)) kim_Es_modes = (0.0d0, 0.0d0) kim_Ep_modes = (0.0d0, 0.0d0) @@ -409,6 +422,8 @@ subroutine kim_run_for_all_modes() kim_kp_modes = 0.0d0 kim_ks_modes = 0.0d0 kim_jpar_modes = (0.0d0, 0.0d0) + kim_jpar_e_modes = (0.0d0, 0.0d0) + kim_jpar_i_modes = (0.0d0, 0.0d0) ! ------------------------------------------------------- ! 2. Loop over modes: solve and store @@ -422,9 +437,14 @@ subroutine kim_run_for_all_modes() ! Clean up EBdat from previous solve call deallocate_EBdat() - ! Run electromagnetic solver (grids/equilibrium from kim_initialize) - ! NOTE: For multi-mode support, plasma kp/ks/om_E would need - ! updating here for modes 2+. Currently correct for single-mode. + ! For modes 2+, recompute equilibrium quantities (kp, ks, om_E) + ! which depend on the (m, n) mode numbers. + if (i_mn > 1) then + call calculate_equil(.false.) + call set_plasma_quantities(kim_plasma) + call interpolate_equil(kim_rg_grid%xb) + end if + call kim_instance%run() ! run() fills kernels, solves coupled Poisson-Ampere, and calls @@ -469,11 +489,19 @@ subroutine kim_run_for_all_modes() call interp_complex_profile(kim_npts, kim_r, EBdat%Br, & dim_r, bal_r, kim_Br_modes(:, i_mn)) - ! jpar (parallel current density) + ! jpar (parallel current density: total, electron, ion) if (allocated(EBdat%jpar)) then call interp_complex_profile(kim_npts, kim_r, EBdat%jpar, & dim_r, bal_r, kim_jpar_modes(:, i_mn)) end if + if (allocated(EBdat%jpar_e)) then + call interp_complex_profile(kim_npts, kim_r, EBdat%jpar_e, & + dim_r, bal_r, kim_jpar_e_modes(:, i_mn)) + end if + if (allocated(EBdat%jpar_i)) then + call interp_complex_profile(kim_npts, kim_r, EBdat%jpar_i, & + dim_r, bal_r, kim_jpar_i_modes(:, i_mn)) + end if deallocate(kim_r) @@ -698,6 +726,8 @@ subroutine deallocate_EBdat() if (allocated(EBdat%Phi_MA_ideal)) deallocate(EBdat%Phi_MA_ideal) if (allocated(EBdat%Phi_MA_asymptotic)) deallocate(EBdat%Phi_MA_asymptotic) if (allocated(EBdat%jpar)) deallocate(EBdat%jpar) + if (allocated(EBdat%jpar_e)) deallocate(EBdat%jpar_e) + if (allocated(EBdat%jpar_i)) deallocate(EBdat%jpar_i) end subroutine deallocate_EBdat @@ -735,22 +765,28 @@ end subroutine interp_complex_profile subroutine kim_load_vacuum_fields() !! Extract vacuum Br from KiLCA vacuum solution (vac_cd_ptr) - !! onto the balance grid and store in kim_vac_Br. + !! onto the balance grid and store in kim_vac_Br for each mode. use wave_code_data, only: dim_r, r, dim_mn, m_vals, n_vals, & vac_cd_ptr, Br, Bz implicit none - ! Extract vacuum Br for first mode (single-mode for now). + integer :: k + + ! Extract vacuum Br for each mode. ! get_wave_fields_from_wave_code fills Br; we use Bz as dummy ! for all other field components. - call get_wave_fields_from_wave_code(vac_cd_ptr(1), dim_r, r, & - m_vals(1), n_vals(1), Bz, Bz, Bz, Bz, Bz, Br, Bz, Bz, Bz, Bz) - kim_vac_Br = Br + do k = 1, dim_mn + call get_wave_fields_from_wave_code(vac_cd_ptr(k), dim_r, r, & + m_vals(k), n_vals(k), Bz, Bz, Bz, Bz, Bz, Br, Bz, Bz, Bz, Bz) + kim_vac_Br(:, k) = Br + end do - write(*,*) 'KIM adapter: loaded vacuum Br from KiLCA' - write(*,*) ' |Br_vac| at r_min = ', abs(kim_vac_Br(1)), ' at r=', r(1) - write(*,*) ' |Br_vac| at r_max = ', abs(kim_vac_Br(dim_r)), ' at r=', r(dim_r) + write(*,*) 'KIM adapter: loaded vacuum Br from KiLCA for ', dim_mn, ' modes' + do k = 1, dim_mn + write(*,*) ' mode ', k, ': |Br_vac| at r_max = ', & + abs(kim_vac_Br(dim_r, k)), ' at r=', r(dim_r) + end do call kim_check_domain_consistency() @@ -809,13 +845,13 @@ subroutine kim_check_domain_consistency() w = 0.0d0 end if - Br_vac_at_bc = (1.0d0 - w) * kim_vac_Br(i_lo) & - + w * kim_vac_Br(i_hi) + Br_vac_at_bc = (1.0d0 - w) * kim_vac_Br(i_lo, 1) & + + w * kim_vac_Br(i_hi, 1) - write(*,*) ' Vacuum Br at r_boundary (linear interp):' + write(*,*) ' Vacuum Br at r_boundary (linear interp, mode 1):' write(*,*) ' r_lo=', r(i_lo), ' r_hi=', r(i_hi), ' w=', w - write(*,*) ' |Br_vac(r_lo)| = ', abs(kim_vac_Br(i_lo)) - write(*,*) ' |Br_vac(r_hi)| = ', abs(kim_vac_Br(i_hi)) + write(*,*) ' |Br_vac(r_lo)| = ', abs(kim_vac_Br(i_lo, 1)) + write(*,*) ' |Br_vac(r_hi)| = ', abs(kim_vac_Br(i_hi, 1)) write(*,*) ' |Br_vac(r_bc)| = ', abs(Br_vac_at_bc) write(*,*) ' Re(Br) = ', real(Br_vac_at_bc) write(*,*) ' Im(Br) = ', aimag(Br_vac_at_bc) @@ -829,17 +865,19 @@ subroutine kim_check_domain_consistency() end subroutine kim_check_domain_consistency subroutine kim_get_current_densities(i_mn) - !! Copy per-mode stored parallel current density from - !! kim_jpar_modes into the wave_code_data J arrays for mode i_mn. - !! KIM provides total jpar; it goes into Jpi, rest is zeroed. + !! Copy per-mode stored parallel current densities from + !! kim_jpar_e_modes / kim_jpar_i_modes into the wave_code_data + !! J arrays for mode i_mn. KIM computes species-resolved + !! currents via separate electron and ion kernels. use wave_code_data, only: Jpi, Jpe, Jri, Jre, Jsi, Jse implicit none integer, intent(in) :: i_mn - Jpi = kim_jpar_modes(:, i_mn) - Jpe = (0.0d0, 0.0d0) + Jpe = kim_jpar_e_modes(:, i_mn) + Jpi = kim_jpar_i_modes(:, i_mn) + ! KIM does not compute radial/perpendicular current components Jri = (0.0d0, 0.0d0) Jre = (0.0d0, 0.0d0) Jsi = (0.0d0, 0.0d0) @@ -874,11 +912,11 @@ subroutine apply_vacuum_continuation(i_mn, nrad, bal_r) write(*,*) ' Last KIM point: r=', bal_r(i_last_kim), & ' |Br_KIM|=', abs(kim_Br_modes(i_last_kim, i_mn)) write(*,*) ' First vac point: r=', bal_r(i_last_kim+1), & - ' |Br_vac|=', abs(kim_vac_Br(i_last_kim+1)) + ' |Br_vac|=', abs(kim_vac_Br(i_last_kim+1, i_mn)) write(*,*) ' Vac Br at last KIM r: |Br_vac|=', & - abs(kim_vac_Br(i_last_kim)) + abs(kim_vac_Br(i_last_kim, i_mn)) write(*,*) ' Value jump: |Br_KIM - Br_vac| = ', & - abs(kim_Br_modes(i_last_kim, i_mn) - kim_vac_Br(i_last_kim)) + abs(kim_Br_modes(i_last_kim, i_mn) - kim_vac_Br(i_last_kim, i_mn)) end if do i = i_last_kim + 1, nrad @@ -889,7 +927,7 @@ subroutine apply_vacuum_continuation(i_mn, nrad, bal_r) kim_Et_modes(i, i_mn) = (0.0d0, 0.0d0) kim_Ez_modes(i, i_mn) = (0.0d0, 0.0d0) ! Total Br = vacuum Br beyond plasma - kim_Br_modes(i, i_mn) = kim_vac_Br(i) + kim_Br_modes(i, i_mn) = kim_vac_Br(i, i_mn) ! Wave vectors not physical in vacuum kim_kp_modes(i, i_mn) = 0.0d0 kim_ks_modes(i, i_mn) = 0.0d0 From 793174c873e2090d51c62bd5e9278a17c45cc82f Mon Sep 17 00:00:00 2001 From: Markus Markl Date: Tue, 24 Mar 2026 12:20:46 +0100 Subject: [PATCH 17/20] feat(QL-Balance): add curl(B) Jpar method for antenna factor Add compute_jpar_curlB subroutine that computes parallel current density from Ampere's law (J = curl(B)/(4*pi)) matching the Matlab KiLCA postprocessor convention. Both conductivity-based (sigma*E) and curl(B) Jpar are computed and logged; jpar_method namelist option ('conductivity' or 'curlB') selects which feeds the antenna factor calculation. Co-Authored-By: Claude Opus 4.6 (1M context) --- QL-Balance/src/base/control_mod.f90 | 1 + QL-Balance/src/base/get_dql.f90 | 61 +++++++++++++++- .../src/base/integrate_parallel_current.f90 | 73 +++++++++++++++++++ QL-Balance/src/base/read_config.f90 | 6 +- 4 files changed, 137 insertions(+), 4 deletions(-) diff --git a/QL-Balance/src/base/control_mod.f90 b/QL-Balance/src/base/control_mod.f90 index 940ff9ce..53e6ce7b 100644 --- a/QL-Balance/src/base/control_mod.f90 +++ b/QL-Balance/src/base/control_mod.f90 @@ -18,6 +18,7 @@ module control_mod logical :: debug_mode logical :: suppression_mode ! added by Markus Markl 13.04.2021 logical :: misalign_diffusion ! trigger the calculation and addition of the diffusion due to misaligned equipotentials and flux surfaces + character(len=32) :: jpar_method = 'conductivity' ! 'conductivity' or 'curlB' logical :: diagnostics_output logical :: write_gyro_current integer :: irf diff --git a/QL-Balance/src/base/get_dql.f90 b/QL-Balance/src/base/get_dql.f90 index ba88845a..967b7c45 100644 --- a/QL-Balance/src/base/get_dql.f90 +++ b/QL-Balance/src/base/get_dql.f90 @@ -14,7 +14,8 @@ subroutine get_dql , r_resonant, d11_misalign, Es_pert_flux, Ipar use plasma_parameters use baseparam_mod, only: Z_i, e_charge, am, p_mass, c, e_mass, ev, rtor, pi, rsepar - use control_mod, only: irf, suppression_mode, misalign_diffusion, type_of_run, wave_code + use control_mod, only: irf, suppression_mode, misalign_diffusion, type_of_run, wave_code, & + jpar_method use time_evolution, only: save_prof_time_step, time_ind, br_formfactor, br_vac_res use h5mod use wave_code_data @@ -50,6 +51,13 @@ subroutine get_dql integer :: ibrabsres, ind_begin_interp, ind_end_interp real(dp) :: MI_width + ! curl(B) current density + complex(dp), dimension(:), allocatable :: Jpar_curlB + complex(dp) :: Ipar_conductivity, Ipar_curlB + ! Temporary B field arrays for curl(B) computation (need fresh FLRE fields) + complex(dp), dimension(:), allocatable :: Br_flre, Bt_flre, Bz_flre + complex(dp), dimension(:), allocatable :: zeros_dim_r + allocate (dqle11_loc(npoib)) allocate (dqle12_loc(npoib)) allocate (dqle21_loc(npoib)) @@ -62,6 +70,10 @@ subroutine get_dql if (.not. allocated(d11_misalign)) allocate (d11_misalign(npoib)) if (.not. allocated(Es_pert_flux)) allocate (Es_pert_flux(npoib)) if (.not. allocated(Es_pert_flux_temp)) allocate (Es_pert_flux_temp(npoib)) + allocate(Jpar_curlB(npoib)) + allocate(Br_flre(npoib), Bt_flre(npoib), Bz_flre(npoib)) + allocate(zeros_dim_r(npoib)) + zeros_dim_r = (0.0d0, 0.0d0) dqle11_loc = 0.0d0 dqle12_loc = 0.0d0 @@ -319,6 +331,7 @@ subroutine get_dql dqli21_loc = dqli21_loc + di21*spec_weight dqli22_loc = dqli22_loc + di22*spec_weight + ! --- Conductivity-based current (sigma*E) --- select case (trim(wave_code)) case ('KiLCA') call get_current_densities_from_wave_code(flre_cd_ptr(i_mn), dim_r, r, & @@ -329,7 +342,48 @@ subroutine get_dql error stop "Unknown wave_code: "//trim(wave_code) end select - call integrate_parallel_current(dim_r, r, Jpe, Jpi, r_resonant(i_mn), Ipar) + write(*,*) ' --- Ipar (conductivity, sigma*E) ---' + call integrate_parallel_current(dim_r, r, Jpe, Jpi, r_resonant(i_mn), Ipar_conductivity) + + ! --- curl(B)-based current (Ampere's law) --- + ! Reload full FLRE B field components (Bt, Bz may have been overwritten above) + select case (trim(wave_code)) + case ('KiLCA') + call get_wave_fields_from_wave_code(flre_cd_ptr(i_mn), dim_r, r, & + m_vals(i_mn), n_vals(i_mn), Er, Es, Ep, Et, Ez, Br_flre, Bs, Bp, Bt_flre, Bz_flre) + case ('KIM') + Br_flre = kim_Br_modes(:, i_mn) + ! KIM stores Et (=Btheta) and Ez via kim_get_wave_fields + ! but Bt/Bz are not stored separately — for KIM, the EM solver + ! does not output cylindrical B components other than Br. + ! Use the conductivity-based current for KIM. + Bt_flre = (0.0d0, 0.0d0) + Bz_flre = (0.0d0, 0.0d0) + end select + + call compute_jpar_curlB(dim_r, r, Br_flre, Bt_flre, Bz_flre, & + B0t, B0z, B0, m_vals(i_mn), n_vals(i_mn), rtor, & + Jpar_curlB) + ! compute_jpar_curlB returns Jpar in c=1 Gaussian units. + ! Convert to full CGS (multiply by c) so it goes through the + ! same integration/antenna_factor pipeline as conductivity Jpar. + Jpar_curlB = Jpar_curlB * c + + write(*,*) ' --- Ipar (curlB, Ampere law) ---' + call integrate_parallel_current(dim_r, r, Jpar_curlB, zeros_dim_r, & + r_resonant(i_mn), Ipar_curlB) + + ! Select which Ipar to use for antenna factor computation + select case (trim(jpar_method)) + case ('curlB') + Ipar = Ipar_curlB + write(*,*) ' Using curlB Ipar for antenna factor' + case ('conductivity') + Ipar = Ipar_conductivity + write(*,*) ' Using conductivity Ipar for antenna factor' + case default + error stop "Unknown jpar_method: "//trim(jpar_method) + end select end do @@ -366,6 +420,9 @@ subroutine get_dql deallocate (dqli21_loc); deallocate (dqli22_loc); deallocate (formfactor) + deallocate (Jpar_curlB) + deallocate (Br_flre, Bt_flre, Bz_flre) + deallocate (zeros_dim_r) call calc_parallel_current_directly call calc_ion_parallel_current_directly diff --git a/QL-Balance/src/base/integrate_parallel_current.f90 b/QL-Balance/src/base/integrate_parallel_current.f90 index 934ed50a..fffc09e5 100644 --- a/QL-Balance/src/base/integrate_parallel_current.f90 +++ b/QL-Balance/src/base/integrate_parallel_current.f90 @@ -165,3 +165,76 @@ subroutine integrate_parallel_current(dim_r, r, Jpe, Jpi, r_res_in, Ipar) write(*,*) ' Re(Ipar)= ', real(Ipar), ' Im(Ipar)= ', aimag(Ipar) end subroutine + + +subroutine compute_jpar_curlB(dim_r, r, Br_in, Bt_in, Bz_in, & + B0t_in, B0z_in, B0_in, & + m_mode, n_mode, Rtor, & + Jpar_curlB) +! +! Compute parallel current density from curl(B)/(4*pi) (Ampere's law). +! This matches the Matlab KiLCA_postprocessor convention. +! +! In Gaussian CGS with c=1: +! J = curl(B) / (4*pi) +! Cylindrical components: +! Jr = i/(4*pi) * (kth * Bz - kz * Bth) +! Jth = 1/(4*pi) * (i*kz * Br - dBz/dr) +! Jz = 1/(4*pi) * (Bth/r + dBth/dr - i*kth * Br) +! Parallel projection: +! Jpar = (Jth * B0th + Jz * B0z) / B0 +! + use QLBalance_kinds, only: dp + + implicit none + + integer, intent(in) :: dim_r + real(dp), intent(in) :: r(dim_r) + complex(dp), intent(in) :: Br_in(dim_r), Bt_in(dim_r), Bz_in(dim_r) + real(dp), intent(in) :: B0t_in(dim_r), B0z_in(dim_r), B0_in(dim_r) + integer, intent(in) :: m_mode, n_mode + real(dp), intent(in) :: Rtor + complex(dp), intent(out) :: Jpar_curlB(dim_r) + + complex(dp) :: Jr, Jth, Jz + real(dp) :: kth, kz, fourpi + complex(dp) :: dBt_dr, dBz_dr + real(dp) :: dr_l, dr_r + integer :: i + + fourpi = 16.0d0 * atan(1.0d0) + kz = dble(n_mode) / Rtor + + do i = 1, dim_r + kth = dble(m_mode) / r(i) + + ! Numerical r-derivatives of Bt and Bz (central differences, one-sided at boundaries) + if (i == 1) then + dBt_dr = (Bt_in(2) - Bt_in(1)) / (r(2) - r(1)) + dBz_dr = (Bz_in(2) - Bz_in(1)) / (r(2) - r(1)) + else if (i == dim_r) then + dBt_dr = (Bt_in(dim_r) - Bt_in(dim_r-1)) / (r(dim_r) - r(dim_r-1)) + dBz_dr = (Bz_in(dim_r) - Bz_in(dim_r-1)) / (r(dim_r) - r(dim_r-1)) + else + dr_l = r(i) - r(i-1) + dr_r = r(i+1) - r(i) + ! Non-uniform central difference + dBt_dr = (Bt_in(i+1) * dr_l**2 - Bt_in(i-1) * dr_r**2 & + + Bt_in(i) * (dr_r**2 - dr_l**2)) & + / (dr_l * dr_r * (dr_l + dr_r)) + dBz_dr = (Bz_in(i+1) * dr_l**2 - Bz_in(i-1) * dr_r**2 & + + Bz_in(i) * (dr_r**2 - dr_l**2)) & + / (dr_l * dr_r * (dr_l + dr_r)) + end if + + ! curl(B)/(4*pi) in cylindrical coordinates (c=1 Gaussian) + Jr = cmplx(0.0d0, 1.0d0, dp) / fourpi * (kth * Bz_in(i) - kz * Bt_in(i)) + Jth = 1.0d0 / fourpi * (cmplx(0.0d0, 1.0d0, dp) * kz * Br_in(i) - dBz_dr) + Jz = 1.0d0 / fourpi * (Bt_in(i) / r(i) + dBt_dr & + - cmplx(0.0d0, 1.0d0, dp) * kth * Br_in(i)) + + ! Parallel projection + Jpar_curlB(i) = (Jth * B0t_in(i) + Jz * B0z_in(i)) / B0_in(i) + end do + +end subroutine compute_jpar_curlB diff --git a/QL-Balance/src/base/read_config.f90 b/QL-Balance/src/base/read_config.f90 index 52c671a2..1d00e717 100644 --- a/QL-Balance/src/base/read_config.f90 +++ b/QL-Balance/src/base/read_config.f90 @@ -4,7 +4,8 @@ subroutine read_config readfromtimestep, temperature_limit, gyro_current_study, & misalign_diffusion, equil_path, ihdf5IO, wave_code, & kim_config_path, kim_profiles_from_balance, & - kim_n_modes, kim_m_list, kim_n_list + kim_n_modes, kim_m_list, kim_n_list, & + jpar_method use grid_mod, only: rmin, rmax, npoimin, gg_factor, gg_width, gg_r_res, iboutype, rb_cut_in, & re_cut_in, rb_cut_out, re_cut_out use h5mod, only: path2inp, path2out, path2time @@ -28,7 +29,7 @@ subroutine read_config misalign_diffusion, equil_path, ihdf5IO, type_of_run, wave_code, & set_constant_time_step, constant_time_step, urelax, kim_config_path, & kim_profiles_from_balance, kim_n_modes, kim_m_list, kim_n_list, & - I_par_toroidal + I_par_toroidal, jpar_method ! read the parameters from namelist file open (newunit=u, file=config_file, status="old", action="read", iostat=ios) @@ -85,6 +86,7 @@ subroutine read_config write (*, "(A,L0)") " set_constant_time_step = ", set_constant_time_step write (*, "(A,ES15.8,A)") " constant_time_step = ", constant_time_step, " s" write (*, "(A,ES15.8)") " urelax = ", urelax + write (*, "(A,A)") " jpar_method = ", trim(adjustl(jpar_method)) write (*, "(A,A)") " kim_config_path = ", trim(adjustl(kim_config_path)) write (*, "(A,L0)") " kim_profiles_from_balance = ", kim_profiles_from_balance write (*, "(A,I0)") " kim_n_modes = ", kim_n_modes From dbb4e4d8a1d4d179e3b0d0f85604d8fded1e7f47 Mon Sep 17 00:00:00 2001 From: Markus Markl Date: Tue, 24 Mar 2026 12:58:13 +0100 Subject: [PATCH 18/20] fix(QL-Balance): restrict curlB jpar_method to KiLCA wave code Add error stop when jpar_method='curlB' is used with wave_code other than KiLCA, since KIM does not output Bt/Bz components needed for the curl(B) computation. Co-Authored-By: Claude Opus 4.6 (1M context) --- QL-Balance/src/base/get_dql.f90 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/QL-Balance/src/base/get_dql.f90 b/QL-Balance/src/base/get_dql.f90 index 967b7c45..a69bc298 100644 --- a/QL-Balance/src/base/get_dql.f90 +++ b/QL-Balance/src/base/get_dql.f90 @@ -376,6 +376,10 @@ subroutine get_dql ! Select which Ipar to use for antenna factor computation select case (trim(jpar_method)) case ('curlB') + if (trim(wave_code) /= 'KiLCA') then + error stop "jpar_method='curlB' requires wave_code='KiLCA' " & + // "(KIM does not output Bt/Bz components)" + end if Ipar = Ipar_curlB write(*,*) ' Using curlB Ipar for antenna factor' case ('conductivity') From c9d8c74c3c293114e7910e03be7688bfa9957692 Mon Sep 17 00:00:00 2001 From: Markus Markl Date: Wed, 25 Mar 2026 11:48:13 +0100 Subject: [PATCH 19/20] fix(KIM): correct link order for ddeabm/slatec in test targets ddeabm depends on xermsg_ from slatec, so ddeabm must precede slatec in the static library link order for the linker to resolve the symbol. Co-Authored-By: Claude Opus 4.6 (1M context) --- KIM/tests/CMakeLists.txt | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/KIM/tests/CMakeLists.txt b/KIM/tests/CMakeLists.txt index fb744268..fbdce1fe 100644 --- a/KIM/tests/CMakeLists.txt +++ b/KIM/tests/CMakeLists.txt @@ -6,7 +6,7 @@ add_executable(test_integration_methods ${CMAKE_SOURCE_DIR}/KIM/tests/test_integ set_target_properties(test_integration_methods PROPERTIES OUTPUT_NAME test_integration_methods.x RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/tests/") -target_link_libraries(test_integration_methods KIM_lib "${KiLCA_lib_path}" lapack cerf slatec ddeabm) +target_link_libraries(test_integration_methods KIM_lib "${KiLCA_lib_path}" lapack cerf ddeabm slatec) add_test(NAME test_integration_methods COMMAND ${CMAKE_BINARY_DIR}/tests/test_integration_methods.x) @@ -15,7 +15,7 @@ add_executable(test_quadpack_qagi ${CMAKE_SOURCE_DIR}/KIM/tests/test_quadpack_qa set_target_properties(test_quadpack_qagi PROPERTIES OUTPUT_NAME test_quadpack_qagi.x RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/tests/") -target_link_libraries(test_quadpack_qagi KIM_lib lapack cerf slatec ddeabm) +target_link_libraries(test_quadpack_qagi KIM_lib lapack cerf ddeabm slatec) add_test(NAME test_quadpack_qagi COMMAND ${CMAKE_BINARY_DIR}/tests/test_quadpack_qagi.x) @@ -48,6 +48,6 @@ add_executable(test_ampere_matrices ${CMAKE_SOURCE_DIR}/KIM/tests/test_ampere_ma set_target_properties(test_ampere_matrices PROPERTIES OUTPUT_NAME test_ampere_matrices.x RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/tests/") -target_link_libraries(test_ampere_matrices KIM_lib "${KiLCA_lib_path}" lapack cerf slatec ddeabm) +target_link_libraries(test_ampere_matrices KIM_lib "${KiLCA_lib_path}" lapack cerf ddeabm slatec) add_test(NAME test_ampere_matrices COMMAND ${CMAKE_BINARY_DIR}/tests/test_ampere_matrices.x) From 4c9eac2dd39c6188701f6b5800a9d3b6c3b068a5 Mon Sep 17 00:00:00 2001 From: Markus Markl Date: Wed, 25 Mar 2026 12:50:28 +0100 Subject: [PATCH 20/20] fix(test): use main branch namelist template for golden record The golden record test builds main's ql-balance.x but was using the current branch's balance_conf.nml template, which may contain namelist variables unknown to the main binary (e.g. wave_code, kim_*). Pass main_ref's own template to setup_runfolder so the generated config only includes variables the main-branch binary recognises. Co-Authored-By: Claude Opus 4.6 (1M context) --- test/ql-balance/golden_record/ensure_golden.py | 11 +++++++++-- test/ql-balance/golden_record/setup_runfolder.py | 10 ++++++++-- 2 files changed, 17 insertions(+), 4 deletions(-) diff --git a/test/ql-balance/golden_record/ensure_golden.py b/test/ql-balance/golden_record/ensure_golden.py index 8059deb6..2a6f3b8d 100755 --- a/test/ql-balance/golden_record/ensure_golden.py +++ b/test/ql-balance/golden_record/ensure_golden.py @@ -98,13 +98,20 @@ def build_main_ref() -> Path: def setup_runfolder() -> Path: - """Set up the runfolder using the user-provided setup script.""" + """Set up the runfolder using the user-provided setup script. + + Uses main branch's namelist template so the config only contains + variables that the main-branch binary recognises. + """ # Clean up any existing runfolder if RUNFOLDER_DIR.exists(): shutil.rmtree(RUNFOLDER_DIR) RUNFOLDER_DIR.mkdir(parents=True) - setup_runfolder_external(str(RUNFOLDER_DIR)) + main_nml_template = str( + MAIN_REF_DIR / "QL-Balance" / "namelists" / "balance_conf.nml" + ) + setup_runfolder_external(str(RUNFOLDER_DIR), nml_template=main_nml_template) return RUNFOLDER_DIR diff --git a/test/ql-balance/golden_record/setup_runfolder.py b/test/ql-balance/golden_record/setup_runfolder.py index be1e720a..66661364 100755 --- a/test/ql-balance/golden_record/setup_runfolder.py +++ b/test/ql-balance/golden_record/setup_runfolder.py @@ -47,7 +47,7 @@ from utility import create_parabolic_profiles_from_res_surf -def setup_runfolder(run_path: str) -> None: +def setup_runfolder(run_path: str, nml_template: str = "") -> None: """ Set up a runfolder in the target directory with all necessary input files. @@ -57,6 +57,12 @@ def setup_runfolder(run_path: str) -> None: Args: run_path: The directory where the runfolder should be created. This directory already exists when this function is called. + nml_template: Optional path to a balance_conf.nml template file. + When running against a different branch's binary (e.g. + the golden record's main-branch build), pass that + branch's template so the namelist only contains + variables the binary recognises. Defaults to the + current source tree's template. Raises: RuntimeError: If setup fails for any reason. @@ -81,7 +87,7 @@ def setup_runfolder(run_path: str) -> None: Btor = -17000 # toroidal magnetic field on axis in Gauss bi = QL_Balance_interface(run_path=run_path, shot=0, time=0, name="test") - bi.read_config_nml() + bi.read_config_nml(path=nml_template) bi.set_type_of_run(run_type="TimeEvolution") bi.set_modes(m_mode=mpol, n_mode=ntor)