Add radial background periodization

KIM/CMakeLists.txt

@@ -28,7 +28,6 @@ enable_testing()
 add_subdirectory(src)
 # fortnum_amos_compat is provided by the top-level cmake/Dependencies.cmake; do not add its subdir here
 add_subdirectory(src/math/libcerf-main)
-add_subdirectory(src/tests)
 add_subdirectory(tests)
 
 set(CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake" ${CMAKE_MODULE_PATH})

KIM/src/CMakeLists.txt

@@ -9,7 +9,6 @@ add_library(KIM_lib STATIC "${KIM_setup}"
                         "${KIM_general}"
                         "${KIM_poisson}"
                         "${KIM_util}"
-                        #"${KIM_kernels}"
                         "${KIM_math}"
                         "${KIM_fokkerplanck}"
                         "${KIM_dispersion}"

KIM/src/CMakeSources.in

@@ -39,15 +39,10 @@ set(KIM_util util/findIndex.f90
             util/gauss_quadrature_m.f90
 )
 
-#set(KIM_kernels kernels/integrands.f90
-                ##kernels/kernel_functions.f90
-                #kernels/cut_off_integration.f90
-                #kernels/kernel_mod.f90
-#)
-
 set(KIM_background_equilibrium background_equilibrium/species_mod.f90
                             background_equilibrium/calculate_equil.f90
                             background_equilibrium/profile_input_m.f90
+                            background_equilibrium/periodize_m.f90
 )
 
 set(KIM_dispersion  dispersion/fun_input.f90
@@ -59,6 +54,7 @@ set(KIM_poisson electrostatic_poisson/poisson.f90
                 electrostatic_poisson/flr2_benchmark.f90
                 electrostatic_poisson/solve_poisson.f90
                 kernels/kernel.f90
+                kernels/kernel_mod.f90
                 kernels/integrands.f90
                 electrostatic_poisson/fields_mod.f90
                 kernels/integrals.f90

KIM/src/background_equilibrium/periodize_m.f90

@@ -0,0 +1,107 @@
+module periodize_m
+    ! Radial periodization for the enforced-periodicity solver (issue #175,
+    ! Kasilov note 2026-07-04): background profiles stay untouched inside
+    ! the resonant layer |r - r_m| <= dr_layer and are blended with their
+    ! period-shifted copies inside the transition zones so that the result
+    ! is periodic with L = 2 (dr_layer + dr_transition) and keeps all
+    ! derivatives continuous. Gradients must be computed from the
+    ! periodized profile afterwards; periodizing a gradient is wrong in
+    ! the transition zone (mathematica/27 in the programme repository).
+
+    use KIM_kinds_m, only: dp
+    use constants_m, only: pi
+
+    implicit none
+
+contains
+
+    ! C-infinity transition weight of the prototype localizer: 1 below x1,
+    ! 0 above x2, essential-singularity decay at both ends.
+    real(dp) function localizer_weight(x1, x2, x)
+
+        implicit none
+
+        real(dp), intent(in) :: x1, x2, x
+
+        real(dp) :: t
+
+        t = (x - x1)/(x2 - x1)
+        if (t <= 0.0d0) then
+            localizer_weight = 1.0d0
+        else if (t >= 1.0d0) then
+            localizer_weight = 0.0d0
+        else
+            localizer_weight = exp(-2.0d0*pi/(1.0d0 - t)*exp(-sqrt(2.0d0)/t))
+        end if
+
+    end function localizer_weight
+
+    ! Periodized profile value at x from a tabulation (r_grid, f) of the
+    ! original profile. The tabulation must cover three half-periods on
+    ! both sides of r_mid because the blend samples the shifted copies.
+    real(dp) function periodized_profile_value(r_grid, f, r_mid, dr_layer, &
+                                               dr_transition, x)
+
+        implicit none
+
+        real(dp), intent(in) :: r_grid(:), f(:)
+        real(dp), intent(in) :: r_mid, dr_layer, dr_transition, x
+
+        real(dp) :: half_period, period, x_leftboundary, x_inperiod
+        real(dp) :: value_in, value_left, value_right, weight
+
+        if (dr_layer <= 0.0d0 .or. dr_transition <= 0.0d0) then
+            error stop "periodize_m: layer and transition widths must be positive"
+        end if
+        half_period = dr_layer + dr_transition
+        period = 2.0d0*half_period
+        if (r_grid(1) > r_mid - 3.0d0*half_period) then
+            error stop "periodize_m: tabulation does not cover the left copies"
+        end if
+        if (r_grid(size(r_grid)) < r_mid + 3.0d0*half_period) then
+            error stop "periodize_m: tabulation does not cover the right copies"
+        end if
+
+        x_leftboundary = r_mid - half_period
+        x_inperiod = x_leftboundary + modulo(x - x_leftboundary, period)
+
+        value_in = tabulated_value(r_grid, f, x_inperiod)
+        value_left = tabulated_value(r_grid, f, x_inperiod - period)
+        value_right = tabulated_value(r_grid, f, x_inperiod + period)
+
+        weight = localizer_weight(r_mid + dr_layer, &
+                                  r_mid + dr_layer + 2.0d0*dr_transition, &
+                                  x_inperiod)
+        periodized_profile_value = value_in*weight + value_left*(1.0d0 - weight)
+
+        weight = localizer_weight(r_mid - dr_layer - 2.0d0*dr_transition, &
+                                  r_mid - dr_layer, x_inperiod)
+        periodized_profile_value = periodized_profile_value*(1.0d0 - weight) &
+                                   + value_right*weight
+
+    end function periodized_profile_value
+
+    real(dp) function tabulated_value(r_grid, f, x)
+
+        implicit none
+
+        real(dp), intent(in) :: r_grid(:), f(:), x
+
+        integer, parameter :: nlagr = 4
+        integer, parameter :: nder = 0
+        real(dp) :: coef(0:nder, nlagr)
+        integer :: ir, ibeg, iend
+
+        call binsrc(r_grid, 1, size(r_grid), x, ir)
+        ibeg = max(1, ir - nlagr/2)
+        iend = ibeg + nlagr - 1
+        if (iend > size(r_grid)) then
+            iend = size(r_grid)
+            ibeg = iend - nlagr + 1
+        end if
+        call plag_coeff(nlagr, nder, x, r_grid(ibeg:iend), coef)
+        tabulated_value = sum(coef(0, :)*f(ibeg:iend))
+
+    end function tabulated_value
+
+end module periodize_m

KIM/src/kernels/kernel_mod.f90

@@ -0,0 +1,209 @@
+module kernels_m
+    ! Continuous-Fourier KIM kernels G(k_r, k_r', r_g) restored from
+    ! kernel_mod.f90 at 1ae0aeeb~1 for the enforced-periodicity solver:
+    ! the periodic matrix element is K_{m,m'} = (2 pi / L) times the
+    ! one-period r_g integral of these functions (Kasilov, "Enforced
+    ! periodicity", 2026-07-04). Fourier phase convention
+    ! exp(i k_r r); CGS-Gaussian units; the rho_phi kernel carries
+    ! 1/(8 pi^2), the rho_B kernel i/(8 pi^2 c).
+
+    use KIM_kinds_m, only: dp
+
+    implicit none
+
+    ! Adiabatic-only response: drop the thermodynamic-force terms while
+    ! keeping the sign, gyroaverage, and Fourier phase of the full
+    ! expression. Matches the artificial_debye_case semantics of the
+    ! hat-basis path; the historical branch flipped the sign and dropped
+    ! the phase and is not restored.
+    logical :: kernel_debye_case = .false.
+
+    real(dp) :: bessel_large_arg_limit = 10d0
+
+    integer, parameter :: nlagr = 4
+    integer, parameter :: nder = 0
+
+contains
+
+    subroutine lagrange_weights(val_rg, ibeg, iend, weights)
+
+        use species_m, only: plasma
+
+        implicit none
+
+        real(dp), intent(in) :: val_rg
+        integer, intent(out) :: ibeg, iend
+        real(dp), intent(out) :: weights(nlagr)
+
+        real(dp) :: coef(0:nder, nlagr)
+        integer :: ir
+
+        call binsrc(plasma%r_grid, 1, plasma%grid_size, val_rg, ir)
+        ibeg = max(1, ir - nlagr/2)
+        iend = ibeg + nlagr - 1
+        if (iend > plasma%grid_size) then
+            iend = plasma%grid_size
+            ibeg = iend - nlagr + 1
+        end if
+        call plag_coeff(nlagr, nder, val_rg, plasma%r_grid(ibeg:iend), coef)
+        weights = coef(0, :)
+
+    end subroutine lagrange_weights
+
+    ! This is without the exp(i k_r(r_g - x_l)) factor
+    complex(dp) function kernel_rho_phi_of_kr_krp_rg(val_kr, val_krp, val_rg)
+
+        use constants_m, only: pi, com_unit
+        use species_m, only: plasma
+        use fortnum_special, only: bessel_in
+
+        implicit none
+
+        real(dp), intent(in) :: val_kr, val_krp, val_rg
+
+        complex(dp), external :: plasma_Z
+
+        real(dp) :: weights(nlagr)
+        integer :: ibeg, iend, sp
+        real(dp) :: eval_bp, eval_bt
+        complex(dp) :: z0_interp, eval_besselI0, eval_besselIm1
+        real(dp) :: vT_interp, omc_interp, ks_interp, kp_interp, A1_interp, &
+                    A2_interp, lambda_D_interp, rhoL_interp, kperp, kperpp
+
+        kernel_rho_phi_of_kr_krp_rg = (0.0d0, 0.0d0)
+
+        call lagrange_weights(val_rg, ibeg, iend, weights)
+        ks_interp = sum(weights*plasma%ks(ibeg:iend))
+        kp_interp = sum(weights*plasma%kp(ibeg:iend))
+        kperp = sqrt(ks_interp**2 + val_kr**2)
+        kperpp = sqrt(ks_interp**2 + val_krp**2)
+
+        do sp = 0, plasma%n_species - 1
+            associate (spec => plasma%spec(sp))
+                vT_interp = sum(weights*spec%vT(ibeg:iend))
+                omc_interp = sum(weights*spec%omega_c(ibeg:iend))
+                A1_interp = sum(weights*spec%A1(ibeg:iend))
+                A2_interp = sum(weights*spec%A2(ibeg:iend))
+                lambda_D_interp = sum(weights*spec%lambda_D(ibeg:iend))
+                z0_interp = sum(weights*spec%z0(ibeg:iend))
+            end associate
+
+            rhoL_interp = vT_interp/abs(omc_interp)
+
+            eval_bp = rhoL_interp**2/2.0d0*(kperp**2 + kperpp**2)
+            eval_bt = rhoL_interp**2*kperp*kperpp
+
+            if (kernel_debye_case) then
+                A1_interp = 0.0d0
+                A2_interp = 0.0d0
+            end if
+
+            if (eval_bt > bessel_large_arg_limit) then
+                ! limit close to magnetic axis (k_s -> infinity) and
+                ! large k_r and k_rp: asymptotics for Bessel I functions
+                eval_besselI0 = exp(eval_bt - eval_bp) &
+                                /sqrt(2.0d0*pi*eval_bt)
+                eval_besselIm1 = exp(-eval_bp + asinh(-1.0d0/eval_bt) &
+                                     + eval_bt*sqrt(1.0d0 + 1.0d0/eval_bt**2)) &
+                                 /sqrt(2.0d0*pi*eval_bt &
+                                       *sqrt(1.0d0 + 1.0d0/eval_bt**2))
+            else
+                eval_besselI0 = bessel_in(0, eval_bt)*exp(-eval_bp)
+                eval_besselIm1 = bessel_in(-1, eval_bt)*exp(-eval_bp)
+            end if
+
+            kernel_rho_phi_of_kr_krp_rg = kernel_rho_phi_of_kr_krp_rg &
+                + 1.0d0/lambda_D_interp**2 &
+                *exp(com_unit*(val_kr - val_krp)*val_rg) &
+                *(-exp(-rhoL_interp**2/2.0d0*(val_kr - val_krp)**2) &
+                  + ks_interp*rhoL_interp/(kp_interp*sqrt(2.0d0)) &
+                  *(A1_interp*eval_besselI0*plasma_Z(z0_interp) &
+                    + A2_interp*(plasma_Z(z0_interp)*eval_besselI0 &
+                                 *(1.0d0 + eval_bp + z0_interp**2) &
+                                 + eval_besselIm1*eval_bt &
+                                 + z0_interp*eval_besselI0)))
+        end do
+
+        kernel_rho_phi_of_kr_krp_rg = kernel_rho_phi_of_kr_krp_rg &
+                                      /(2.0d0**3*pi**2)
+
+    end function kernel_rho_phi_of_kr_krp_rg
+
+    complex(dp) function kernel_rho_B_of_kr_krp_rg(val_kr, val_krp, val_rg)
+
+        use constants_m, only: pi, com_unit, sol
+        use species_m, only: plasma
+        use fortnum_special, only: bessel_in
+
+        implicit none
+
+        real(dp), intent(in) :: val_kr, val_krp, val_rg
+
+        complex(dp), external :: plasma_Z
+
+        real(dp) :: weights(nlagr)
+        integer :: ibeg, iend, sp
+        real(dp) :: eval_bp, eval_bt
+        complex(dp) :: z0_interp, eval_besselI0, eval_besselIm1
+        real(dp) :: vT_interp, omc_interp, ks_interp, kp_interp, A1_interp, &
+                    A2_interp, lambda_D_interp, rhoL_interp, kperp, kperpp
+
+        kernel_rho_B_of_kr_krp_rg = (0.0d0, 0.0d0)
+
+        call lagrange_weights(val_rg, ibeg, iend, weights)
+        ks_interp = sum(weights*plasma%ks(ibeg:iend))
+        kp_interp = sum(weights*plasma%kp(ibeg:iend))
+        kperp = sqrt(ks_interp**2 + val_kr**2)
+        kperpp = sqrt(ks_interp**2 + val_krp**2)
+
+        do sp = 0, plasma%n_species - 1
+            associate (spec => plasma%spec(sp))
+                vT_interp = sum(weights*spec%vT(ibeg:iend))
+                omc_interp = sum(weights*spec%omega_c(ibeg:iend))
+                A1_interp = sum(weights*spec%A1(ibeg:iend))
+                A2_interp = sum(weights*spec%A2(ibeg:iend))
+                lambda_D_interp = sum(weights*spec%lambda_D(ibeg:iend))
+                z0_interp = sum(weights*spec%z0(ibeg:iend))
+            end associate
+
+            rhoL_interp = vT_interp/abs(omc_interp)
+
+            eval_bp = rhoL_interp**2/2.0d0*(kperp**2 + kperpp**2)
+            eval_bt = rhoL_interp**2*kperp*kperpp
+
+            if (kernel_debye_case) then
+                A1_interp = 0.0d0
+                A2_interp = 0.0d0
+            end if
+
+            if (eval_bt > bessel_large_arg_limit) then
+                ! limit close to magnetic axis (k_s -> infinity) and
+                ! large k_r and k_rp: asymptotics for Bessel I functions
+                eval_besselI0 = exp(eval_bt - eval_bp) &
+                                /sqrt(2.0d0*pi*eval_bt)
+                eval_besselIm1 = exp(-eval_bp + asinh(-1.0d0/eval_bt) &
+                                     + eval_bt*sqrt(1.0d0 + 1.0d0/eval_bt**2)) &
+                                 /sqrt(2.0d0*pi*eval_bt &
+                                       *sqrt(1.0d0 + 1.0d0/eval_bt**2))
+            else
+                eval_besselI0 = bessel_in(0, eval_bt)*exp(-eval_bp)
+                eval_besselIm1 = bessel_in(-1, eval_bt)*exp(-eval_bp)
+            end if
+
+            kernel_rho_B_of_kr_krp_rg = kernel_rho_B_of_kr_krp_rg &
+                + exp(com_unit*(val_kr - val_krp)*val_rg) &
+                *vT_interp**2/(lambda_D_interp**2*omc_interp*kp_interp) &
+                *(0.5d0*A1_interp*eval_besselI0 &
+                  *(z0_interp*plasma_Z(z0_interp) + 1.0d0) &
+                  + A2_interp*(0.5d0*eval_besselI0 &
+                               + (z0_interp*plasma_Z(z0_interp) + 1.0d0) &
+                               *((1.0d0 + eval_bp + z0_interp**2)*eval_besselI0 &
+                                 + eval_bp*eval_besselIm1)))
+        end do
+
+        kernel_rho_B_of_kr_krp_rg = kernel_rho_B_of_kr_krp_rg &
+                                    *com_unit/(8.0d0*pi**2*sol)
+
+    end function kernel_rho_B_of_kr_krp_rg
+
+end module kernels_m