feat(QL-Balance): add KIM as alternative wave code backend

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)
@@ -128,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)
@@ -143,6 +151,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)
@@ -151,8 +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
-                        cerf
                         slatec
-                        ddeabm
 )

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

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))

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
@@ -882,7 +883,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
@@ -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

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,18 +64,18 @@ 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(:,:)
         real(dp), allocatable :: hz_xl(:), hth_xl(:)
         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
@@ -224,12 +223,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))
@@ -246,12 +273,27 @@ 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))
+        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')
 
+        ! 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)
 

KIM/src/electrostatic_poisson/fields_mod.f90

@@ -16,6 +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(:)   ! 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(:)
@@ -425,6 +428,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

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)
@@ -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