fix(KIM): apply ion FLR scale consistently

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

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)
 

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 &

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

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