From 1f95a1fd1015a090163a3c3693d748ff0a86ebaa Mon Sep 17 00:00:00 2001 From: marjohma <73427010+marjohma@users.noreply.github.com> Date: Tue, 28 Jul 2026 12:27:44 +0200 Subject: [PATCH] fix(KIM): apply ion FLR scale consistently --- .../background_equilibrium/species_mod.f90 | 3 +- KIM/src/kernels/FP_kernel_plasma_prefacs.f90 | 40 +++++-- KIM/src/kernels/kernel.f90 | 20 +++- KIM/tests/CMakeLists.txt | 10 ++ KIM/tests/test_ion_flr_scale_factor.f90 | 105 ++++++++++++++++++ 5 files changed, 164 insertions(+), 14 deletions(-) create mode 100644 KIM/tests/test_ion_flr_scale_factor.f90 diff --git a/KIM/src/background_equilibrium/species_mod.f90 b/KIM/src/background_equilibrium/species_mod.f90 index 98a52173..9d3cef90 100644 --- a/KIM/src/background_equilibrium/species_mod.f90 +++ b/KIM/src/background_equilibrium/species_mod.f90 @@ -404,8 +404,7 @@ subroutine calculate_thermodynamic_forces_and_susc(plasma_in) use setup_m, only: omega, mphi_max use grid_m, only: rg_grid use KIM_kinds_m, only: dp - use config_m, only: ion_flr_scale_factor, ifunc_model_for_species, & - ion_temperature_gradient_model, & + use config_m, only: ifunc_model_for_species, ion_temperature_gradient_model, & temperature_gradient_force_terms implicit none diff --git a/KIM/src/kernels/FP_kernel_plasma_prefacs.f90 b/KIM/src/kernels/FP_kernel_plasma_prefacs.f90 index 43b5927d..41b45867 100644 --- a/KIM/src/kernels/FP_kernel_plasma_prefacs.f90 +++ b/KIM/src/kernels/FP_kernel_plasma_prefacs.f90 @@ -5,6 +5,24 @@ module FP_kernel_plasma_prefacs_m contains + function FP_flr_scale(spec) result(scale) + + use species_m, only: species_t + use KIM_kinds_m, only: dp + use config_m, only: ion_flr_scale_factor + + implicit none + + type(species_t), intent(in) :: spec + real(dp) :: scale + + ! KIM supports negatively charged electrons and positively charged + ! ions; the artificial ion scale must never alter electron FLR terms. + scale = 1.0_dp + if (spec%Zspec > 0) scale = ion_flr_scale_factor + + end function FP_flr_scale + function FP_kappa_rho_phi(j, spec) result(val) use species_m, only: species_t @@ -25,7 +43,6 @@ function FP_kappa_rho_B(j, spec) result(val) use species_m, only: species_t use KIM_kinds_m, only: dp use constants_m, only: sol - use config_m, only: ion_flr_scale_factor implicit none @@ -39,7 +56,8 @@ function FP_kappa_rho_B(j, spec) result(val) nu = spec%nu_cc(j) omega_c = spec%omega_c_cc(j) - val = - vT**3.0d0 / (lambda**2.0d0 * omega_c * nu * sol) * ion_flr_scale_factor + val = - vT**3.0d0 / (lambda**2.0d0 * omega_c * nu * sol) & + * FP_flr_scale(spec) end function FP_kappa_rho_B @@ -61,7 +79,8 @@ function FP_kappa_j_phi(j, spec) result(val) nu = spec%nu_cc(j) omega_c = spec%omega_c_cc(j) - val = com_unit * vT**3.0d0 / (lambda**2.0d0 * omega_c * nu) + val = com_unit * vT**3.0d0 / (lambda**2.0d0 * omega_c * nu) & + * FP_flr_scale(spec) end function FP_kappa_j_phi @@ -83,7 +102,8 @@ function FP_kappa_j_B(j, spec) result(val) nu = spec%nu_cc(j) omega_c = spec%omega_c_cc(j) - val = - vT**4.0d0 / (lambda**2.0d0 * omega_c * nu * sol) + val = - vT**4.0d0 / (lambda**2.0d0 * omega_c * nu * sol) & + * FP_flr_scale(spec) end function FP_kappa_j_B @@ -107,7 +127,6 @@ function FP_G1_rho_phi(j, spec, mphi) result(val) use species_m, only: species_t, plasma use KIM_kinds_m, only: dp use constants_m, only: com_unit - use config_m, only: ion_flr_scale_factor implicit none @@ -124,7 +143,8 @@ function FP_G1_rho_phi(j, spec, mphi) result(val) I20 = spec%I20_cc(j, mphi) prefactor = FP_kappa_rho_phi(j, spec) * com_unit * spec%vT_cc(j)**2.0d0 / & - (spec%omega_c_cc(j) * spec%nu_cc(j)) * plasma%ks_cc(j) * ion_flr_scale_factor + (spec%omega_c_cc(j) * spec%nu_cc(j)) * plasma%ks_cc(j) & + * FP_flr_scale(spec) val = (I00 * (A1 + A2 * (1.0d0 - mphi)) + 0.5d0 * A2 * I20) * prefactor @@ -136,7 +156,6 @@ function FP_G2_rho_phi(j, spec, mphi) result(val) use species_m, only: species_t, plasma use KIM_kinds_m, only: dp use constants_m, only: com_unit - use config_m, only: ion_flr_scale_factor implicit none @@ -151,7 +170,8 @@ function FP_G2_rho_phi(j, spec, mphi) result(val) I00 = spec%I00_cc(j, mphi) prefactor = com_unit * spec%vT_cc(j)**2.0d0 / & - (spec%omega_c_cc(j) * spec%nu_cc(j)) * plasma%ks_cc(j) * ion_flr_scale_factor + (spec%omega_c_cc(j) * spec%nu_cc(j)) * plasma%ks_cc(j) & + * FP_flr_scale(spec) val = - I00 * A2 * prefactor * FP_kappa_rho_phi(j, spec) @@ -162,7 +182,6 @@ function FP_G3_rho_phi(j, spec, mphi) result(val) use species_m, only: species_t, plasma use KIM_kinds_m, only: dp use constants_m, only: com_unit - use config_m, only: ion_flr_scale_factor implicit none @@ -177,7 +196,8 @@ function FP_G3_rho_phi(j, spec, mphi) result(val) I00 = spec%I00_cc(j, mphi) prefactor = com_unit * spec%vT_cc(j)**2.0d0 / & - (spec%omega_c_cc(j) * spec%nu_cc(j)) * plasma%ks_cc(j) * ion_flr_scale_factor + (spec%omega_c_cc(j) * spec%nu_cc(j)) * plasma%ks_cc(j) & + * FP_flr_scale(spec) val = I00 * A2 * prefactor * FP_kappa_rho_phi(j, spec) diff --git a/KIM/src/kernels/kernel.f90 b/KIM/src/kernels/kernel.f90 index 1294296b..abd1b0b4 100644 --- a/KIM/src/kernels/kernel.f90 +++ b/KIM/src/kernels/kernel.f90 @@ -41,6 +41,22 @@ module kernel_m contains + subroutine set_fp_integration_radius(point, spec, j) + + use integrands_gauss_m, only: integration_point_t + use species_m, only: species_t + + implicit none + + type(integration_point_t), intent(inout) :: point + type(species_t), intent(in) :: spec + integer, intent(in) :: j + + ! rho_L_cc already contains the ion-only artificial FLR scale. + point%rhoT = max(spec%rho_L_cc(j), 0.0_dp) + + end subroutine set_fp_integration_radius + ! Thread-safe wrapper to update loading bar from within OpenMP subroutine update_bar(cur, tot, sc, cr) use loading_bar_m, only: updateLoadingBarWithETA @@ -933,7 +949,7 @@ subroutine FP_calc_kernel_element_ions(l, lp, k_rho_phi, k_rho_B, k_j_phi, k_j_B use FP_kernel_plasma_prefacs_m, only: FP_G0_rho_phi use grid_m, only: Larmor_skip_factor, rg_grid 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 config_m, only: turn_off_ions, turn_off_electrons, artificial_debye_case use grid_m, only: xl_grid use kim_resonances_m, only: r_res use setup_m, only: mphi_max @@ -965,7 +981,7 @@ subroutine FP_calc_kernel_element_ions(l, lp, k_rho_phi, k_rho_B, k_j_phi, k_j_B do j = 1, rg_grid%npts_b-1 int_point%j = j - int_point%rhoT = max(plasma%spec(sigma)%rho_L_cc(j), 0.0d0) * ion_flr_scale_factor + call set_fp_integration_radius(int_point, plasma%spec(sigma), j) ! Debye term (F0) integration if (abs(l-lp)<=1 .and. artificial_debye_case /= 2 & diff --git a/KIM/tests/CMakeLists.txt b/KIM/tests/CMakeLists.txt index 68f88716..eb2e6f22 100644 --- a/KIM/tests/CMakeLists.txt +++ b/KIM/tests/CMakeLists.txt @@ -234,6 +234,16 @@ target_link_libraries(test_collision_scale_z0 KIM_lib kilca_lib lapack cerf fort add_test(NAME test_collision_scale_z0 COMMAND ${CMAKE_BINARY_DIR}/tests/test_collision_scale_z0.x) +add_executable(test_ion_flr_scale_factor + ${CMAKE_SOURCE_DIR}/KIM/tests/test_ion_flr_scale_factor.f90) +set_target_properties(test_ion_flr_scale_factor PROPERTIES + OUTPUT_NAME test_ion_flr_scale_factor.x + RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/tests/") +target_link_libraries(test_ion_flr_scale_factor + KIM_lib kilca_lib lapack cerf fortnum_amos_compat) +add_test(NAME test_ion_flr_scale_factor + COMMAND ${CMAKE_BINARY_DIR}/tests/test_ion_flr_scale_factor.x) + # Regression test for the vendored periodization utility. add_executable(test_periodization ${CMAKE_SOURCE_DIR}/KIM/tests/test_periodization.f90) set_target_properties(test_periodization PROPERTIES diff --git a/KIM/tests/test_ion_flr_scale_factor.f90 b/KIM/tests/test_ion_flr_scale_factor.f90 new file mode 100644 index 00000000..fe85ec36 --- /dev/null +++ b/KIM/tests/test_ion_flr_scale_factor.f90 @@ -0,0 +1,105 @@ +program test_ion_flr_scale_factor + use KIM_kinds_m, only: dp + use species_m, only: plasma, species_t, init_electron_species, & + init_deuterium_species + use config_m, only: ion_flr_scale_factor + use FP_kernel_plasma_prefacs_m, only: FP_kappa_rho_B, FP_G1_rho_phi, & + FP_G2_rho_phi, FP_G3_rho_phi, & + FP_kappa_j_phi, FP_kappa_j_B + use integrands_gauss_m, only: integration_point_t + use kernel_m, only: set_fp_integration_radius + + implicit none + + real(dp), parameter :: scale = 4.0_dp + type(species_t) :: electron, ion + complex(dp) :: electron_base(6), electron_scaled(6) + complex(dp) :: ion_base(6), ion_scaled(6) + type(integration_point_t) :: point + + call init_electron_species(electron) + call init_deuterium_species(ion) + call populate_species(electron) + call populate_species(ion) + allocate(plasma%ks_cc(1)) + plasma%ks_cc = 0.7_dp + + ion_flr_scale_factor = 1.0_dp + call evaluate_scaled_prefactors(electron, electron_base) + call evaluate_scaled_prefactors(ion, ion_base) + + ion_flr_scale_factor = scale + call evaluate_scaled_prefactors(electron, electron_scaled) + call evaluate_scaled_prefactors(ion, ion_scaled) + + call require_close(electron_scaled, electron_base, & + 'ion FLR scale must not change electron prefactors') + call require_close(ion_scaled, scale * ion_base, & + 'ion FLR scale must reach ion prefactors exactly once') + call set_fp_integration_radius(point, ion, 1) + call require_real_close(point%rhoT, ion%rho_L_cc(1), & + 'global geometry must consume stored ion rho_L once') + + print *, 'Ion FLR scale-factor tests PASSED' + +contains + + subroutine populate_species(spec) + type(species_t), intent(inout) :: spec + + allocate(spec%lambda_D_cc(1), spec%vT_cc(1), spec%nu_cc(1)) + allocate(spec%omega_c_cc(1), spec%rho_L_cc(1)) + allocate(spec%A1_cc(1), spec%A2_cc(1)) + allocate(spec%I00_cc(1, 0:0), spec%I20_cc(1, 0:0)) + + spec%lambda_D_cc = 2.0_dp + spec%vT_cc = 3.0_dp + spec%nu_cc = 5.0_dp + spec%omega_c_cc = 7.0_dp + spec%rho_L_cc = 0.37_dp + spec%A1_cc = 0.11_dp + spec%A2_cc = -0.04_dp + spec%I00_cc(:, 0) = cmplx(0.8_dp, -0.2_dp, dp) + spec%I20_cc(:, 0) = cmplx(-0.3_dp, 0.1_dp, dp) + end subroutine populate_species + + subroutine evaluate_scaled_prefactors(spec, values) + type(species_t), intent(in) :: spec + complex(dp), intent(out) :: values(6) + + values(1) = FP_kappa_rho_B(1, spec) + values(2) = FP_G1_rho_phi(1, spec, 0) + values(3) = FP_G2_rho_phi(1, spec, 0) + values(4) = FP_G3_rho_phi(1, spec, 0) + values(5) = FP_kappa_j_phi(1, spec) + values(6) = FP_kappa_j_B(1, spec) + end subroutine evaluate_scaled_prefactors + + subroutine require_close(got, want, label) + complex(dp), intent(in) :: got(:), want(:) + character(*), intent(in) :: label + real(dp) :: error + + error = maxval(abs(got - want) / max(abs(want), 1.0_dp)) + if (error > 1.0e-12_dp) then + print *, 'FAIL: ', label + print *, ' got: ', got + print *, ' want: ', want + print *, ' max relative error: ', error + error stop 'ion FLR scale-factor mismatch' + end if + end subroutine require_close + + subroutine require_real_close(got, want, label) + real(dp), intent(in) :: got, want + character(*), intent(in) :: label + + if (abs(got - want) > 1.0e-12_dp * max(abs(want), 1.0_dp)) then + print *, 'FAIL: ', label + print *, ' got: ', got + print *, ' want: ', want + error stop 'ion FLR geometry mismatch' + end if + end subroutine require_real_close + +end program test_ion_flr_scale_factor