diff --git a/KIM/CMakeLists.txt b/KIM/CMakeLists.txt index d0d1d71c..b45b4584 100644 --- a/KIM/CMakeLists.txt +++ b/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}) diff --git a/KIM/src/CMakeLists.txt b/KIM/src/CMakeLists.txt index 58777f69..ca2c34f3 100644 --- a/KIM/src/CMakeLists.txt +++ b/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}" diff --git a/KIM/src/CMakeSources.in b/KIM/src/CMakeSources.in index 714aab8f..eba286a8 100644 --- a/KIM/src/CMakeSources.in +++ b/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 diff --git a/KIM/src/background_equilibrium/periodize_m.f90 b/KIM/src/background_equilibrium/periodize_m.f90 new file mode 100644 index 00000000..48a91135 --- /dev/null +++ b/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 diff --git a/KIM/src/kernels/kernel_mod.f90 b/KIM/src/kernels/kernel_mod.f90 new file mode 100644 index 00000000..004296ac --- /dev/null +++ b/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 diff --git a/KIM/src/tests/CMakeLists.txt b/KIM/src/tests/CMakeLists.txt deleted file mode 100644 index e69de29b..00000000 diff --git a/KIM/src/tests/kim_init_for_test.f90 b/KIM/src/tests/kim_init_for_test.f90 deleted file mode 100644 index 43d99c7b..00000000 --- a/KIM/src/tests/kim_init_for_test.f90 +++ /dev/null @@ -1,50 +0,0 @@ -subroutine kim_init_for_test - - use KIM_kinds_m, only: dp - use plasma_parameter, only: r_prof, iprof_length, n_prof, Te_prof, Ti_prof, & - Er_prof, q_prof, ni_prof - use equilibrium_m, only: B0z, B0th, B0, hz, hth - use setup_m, only: btor, R0 - - implicit none - - integer :: ierr = 0 - integer :: i - real(dp) :: h - - call kim_read_config - - iprof_length = 100 - - if (.not. allocated(r_prof)) allocate(r_prof(iprof_length), stat=ierr) - if (ierr /= 0) print *, "array: Allocation request denied" - - h = (70.0d0-3.0d0) / (iprof_length-1) - r_prof(1) = 3.0d0 - do i=2, iprof_length - r_prof(i) = r_prof(i-1) + h - end do - - allocate(n_prof(iprof_length), Te_prof(iprof_length), Ti_prof(1, iprof_length), & - Er_prof(iprof_length), q_prof(iprof_length), ni_prof(1, iprof_length)) - allocate(B0z(iprof_length), B0th(iprof_length), B0(iprof_length), hz(iprof_length), hth(iprof_length)) - - n_prof = 2d13 - ni_prof = 2d13 - Te_prof = 1d3 - Ti_prof = 1d3 - Er_prof = 0d0 - q_prof = 1.5d0 - - ! calculate equilibrium B field and J - B0z = btor - B0th = B0z * r_prof /(q_prof * R0) - B0 = sqrt(B0th**2d0 + B0z**2d0) - - hz = B0z / B0 - hth = B0th / B0 - ! calculate quantities used for the kernels, e.g. A1, A2, dndr, omega_c,... - call allocate_backs - call calculate_backs(.false.) - -end subroutine diff --git a/KIM/src/tests/test_kernel_rho_B.f90 b/KIM/src/tests/test_kernel_rho_B.f90 deleted file mode 100644 index a8ca21c7..00000000 --- a/KIM/src/tests/test_kernel_rho_B.f90 +++ /dev/null @@ -1,105 +0,0 @@ -program test_kernel_rho_B - - use kernels_m, only: kernel_rho_B_of_kr_krp_rg - use KIM_kinds_m, only: dp - use plasma_parameter, only: r_prof, n_prof, ni_prof, Te_prof, Ti_prof, iprof_length, Er_prof - use constants_m, only: pi, ev, e_charge - - implicit none - - real(dp) :: kr, krp, rg - complex(dp) :: kernel_value, kernel_test_value - real(dp) :: lambda_De, lambda_Di, lambda_D - real(dp) :: ne_core - real(dp) :: delta_r - integer :: r_ind = 50 - integer :: i - - call kim_init_for_test - - kr = 1.0d0 - krp = 1.0d0 - - rg = r_prof(r_ind) - print *, "r = ", rg - - kernel_value = kernel_rho_B_of_kr_krp_rg(kr, krp, rg) - - if (abs(kernel_value - 0.0d0) < 1.0d-6) then - print *, "Test constant passed" - print *, "" - else - print *, "Test failed, value: ", kernel_value, " should be: ", 0.0d0 - error stop - end if - - ne_core = n_prof(1) - delta_r = r_prof(iprof_length) - r_prof(1) - - do i=2, iprof_length - n_prof(i) = ne_core * (1.0d0 - r_prof(i) / delta_r) - end do - - call calculate_backs(.false.) - - kernel_value = kernel_rho_B_of_kr_krp_rg(kr, krp, rg) - if (abs(abs(kernel_value) - 0.0158746) < 1.0d-3) then - print *, "Test linear n passed" - print *, "" - else - print *, "Test linear n failed, value: ", abs(kernel_value), " should be: 0.0158746" - error stop - end if - - - do i=2, iprof_length - n_prof(i) = ne_core - Te_prof(i) = Te_prof(1) * (1.0d0 - r_prof(i) / delta_r) - end do - - call calculate_backs(.false.) - - kernel_value = kernel_rho_B_of_kr_krp_rg(kr, krp, rg) - if (abs(abs(kernel_value) - 0.0517352) < 1.0d-3) then - print *, "Test linear Te passed" - print *, "" - else - print *, "Test linear Te failed, value: ", abs(kernel_value), " should be: 0.0517352" - error stop - end if - - - do i=2, iprof_length - Te_prof(i) = Te_prof(1) - Er_prof(i) = 0.3d0 - end do - - call calculate_backs(.false.) - - kernel_value = kernel_rho_B_of_kr_krp_rg(kr, krp, rg) - if (abs(abs(kernel_value) - 0.184567) < 1.0d-1) then - print *, "Test constant Er passed" - print *, "" - else - print *, "Test constant Er failed, value: ", abs(kernel_value), " should be: 0.184567" - error stop - end if - - - do i=2, iprof_length - Te_prof(i) = Te_prof(1) * (1.0d0 - r_prof(i) / delta_r) - Er_prof(i) = 0.3d0 - end do - - call calculate_backs(.false.) - - kernel_value = kernel_rho_B_of_kr_krp_rg(kr, krp, rg) - if (abs(abs(kernel_value) - 0.244917) < 1.0d-1) then - print *, "Test constant Er linear Te passed" - print *, "" - else - print *, "Test constant Er linear Te failed, value: ", abs(kernel_value), " should be: 0.244917" - error stop - end if - -end program diff --git a/KIM/src/tests/test_kernel_rho_phi.f90 b/KIM/src/tests/test_kernel_rho_phi.f90 deleted file mode 100644 index dac0bd74..00000000 --- a/KIM/src/tests/test_kernel_rho_phi.f90 +++ /dev/null @@ -1,142 +0,0 @@ -program test_kernel_rho_phi - - use kernels_m, only: kernel_rho_phi_of_kr_krp_rg - use KIM_kinds_m, only: dp - use plasma_parameter, only: r_prof, n_prof, ni_prof, Te_prof, Ti_prof, iprof_length, Er_prof - use constants_m, only: pi, ev, e_charge - - implicit none - - real(dp) :: kr, krp, rg - complex(dp) :: kernel_value, kernel_test_value - real(dp) :: lambda_De, lambda_Di, lambda_D - real(dp) :: ne_core - real(dp) :: delta_r - integer :: r_ind = 50 - integer :: i - - call kim_init_for_test - - kr = 1.0d0 - krp = 1.0d0 - - rg = r_prof(r_ind) - print *, "r = ", rg - - kernel_value = kernel_rho_phi_of_kr_krp_rg(kr, krp, rg) - - print *, "" - print *, "kernel value = ", kernel_value - - lambda_De = sqrt(Te_prof(r_ind) * ev / (4.0d0 * pi * e_charge**2 * n_prof(r_ind))) - lambda_Di = sqrt(Ti_prof(1, r_ind) * ev / (4.0d0 * pi * e_charge**2 * ni_prof(1, r_ind))) - lambda_D = sqrt(1.0d0/(1.0d0/lambda_De**2 + 1.0d0/lambda_Di**2)) - - kernel_test_value = -1.0d0 / (2.0d0**3.0d0 * pi**2.0d0 * lambda_D**2.0d0) - - print *, "lambda De = ", lambda_De - print *, "lambda Di = ", lambda_Di - print *, "lambda D = ", lambda_D - - if (abs(kernel_value - kernel_test_value) < 1.0d-6) then - print *, "Test constant passed" - print *, "" - else - print *, "Test failed, value: ", kernel_value, " should be: ",kernel_test_value - error stop - end if - - !!!!!! - kr = 1.0d0 - krp = 10.0d0 - - rg = r_prof(r_ind) - print *, "kr = ", kr, " krp = ", krp, " r = ", rg - - kernel_value = kernel_rho_phi_of_kr_krp_rg(kr, krp, rg) - - if (abs(abs(kernel_value) - 491.51) < 0.5) then - print *, "Test constant passed" - print *, "" - else - print *, "Test failed, value: ", abs(kernel_value), " should be: 491.5" - print *, "difference is: ", abs(kernel_value - 491.51) - error stop - end if - - kr = 1.0d0 - krp = 1.0d0 - - !!!! - ne_core = n_prof(1) - delta_r = r_prof(iprof_length) - r_prof(1) - - do i=2, iprof_length - n_prof(i) = ne_core * (1.0d0 - r_prof(i) / delta_r) - end do - - call calculate_backs(.false.) - - kernel_value = kernel_rho_phi_of_kr_krp_rg(kr, krp, rg) - if (abs(abs(kernel_value) - 669.313) < 1.0d-2) then - print *, "Test linear n passed" - print *, "" - else - print *, "Test failed, value: ", abs(kernel_value), " should be: 669.313" - error stop - end if - - - do i=2, iprof_length - n_prof(i) = ne_core - Te_prof(i) = Te_prof(1) * (1.0d0 - r_prof(i) / delta_r) - end do - - call calculate_backs(.false.) - - kernel_value = kernel_rho_phi_of_kr_krp_rg(kr, krp, rg) - if (abs(abs(kernel_value) - 1455.25) < 1.455d0) then - print *, "Test linear Te passed" - print *, "" - else - print *, "Test failed, value: ", abs(kernel_value), " should be: 1455.25" - error stop - end if - - - do i=2, iprof_length - Te_prof(i) = Te_prof(1) - Er_prof(i) = 0.3d0 - end do - - call calculate_backs(.false.) - - kernel_value = kernel_rho_phi_of_kr_krp_rg(kr, krp, rg) - if (abs(abs(kernel_value) - 913.4) < 1.0d-1) then - print *, "Test constant Er passed" - print *, "" - else - print *, "Test failed, value: ", abs(kernel_value), " should be: 913.4" - error stop - end if - - - do i=2, iprof_length - Te_prof(i) = Te_prof(1) * (1.0d0 - r_prof(i) / delta_r) - Er_prof(i) = 0.3d0 - end do - - call calculate_backs(.false.) - - kernel_value = kernel_rho_phi_of_kr_krp_rg(kr, krp, rg) - if (abs(abs(kernel_value) - 1450.1) < 1.0d-1) then - print *, "Test constant Er linear Te passed" - print *, "" - else - print *, "Test failed, value: ", abs(kernel_value), " should be: 1450.1" - end if - - - - -end program diff --git a/KIM/tests/CMakeLists.txt b/KIM/tests/CMakeLists.txt index 2793d388..b228e03e 100644 --- a/KIM/tests/CMakeLists.txt +++ b/KIM/tests/CMakeLists.txt @@ -118,3 +118,34 @@ add_test(NAME test_kim_solver_em COMMAND ${CMAKE_BINARY_DIR}/tests/test_kim_solver_em.x) set_tests_properties(test_kim_solver_em PROPERTIES WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/tests/) + +# Restored continuous-Fourier kernels for the enforced-periodicity solver. +add_library(kernel_test_background STATIC + ${CMAKE_SOURCE_DIR}/KIM/tests/kernel_test_background.f90) +target_link_libraries(kernel_test_background KIM_lib) + +add_executable(test_kernel_rho_phi + ${CMAKE_SOURCE_DIR}/KIM/tests/test_kernel_rho_phi.f90) +set_target_properties(test_kernel_rho_phi PROPERTIES + OUTPUT_NAME test_kernel_rho_phi.x + RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/tests/") +target_link_libraries(test_kernel_rho_phi kernel_test_background KIM_lib kilca_lib lapack cerf fortnum_amos_compat) +add_test(NAME test_kernel_rho_phi + COMMAND ${CMAKE_BINARY_DIR}/tests/test_kernel_rho_phi.x) + +add_executable(test_kernel_rho_B + ${CMAKE_SOURCE_DIR}/KIM/tests/test_kernel_rho_B.f90) +set_target_properties(test_kernel_rho_B PROPERTIES + OUTPUT_NAME test_kernel_rho_B.x + RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/tests/") +target_link_libraries(test_kernel_rho_B kernel_test_background KIM_lib kilca_lib lapack cerf fortnum_amos_compat) +add_test(NAME test_kernel_rho_B + COMMAND ${CMAKE_BINARY_DIR}/tests/test_kernel_rho_B.x) + +add_executable(test_periodize ${CMAKE_SOURCE_DIR}/KIM/tests/test_periodize.f90) +set_target_properties(test_periodize PROPERTIES + OUTPUT_NAME test_periodize.x + RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/tests/") +target_link_libraries(test_periodize KIM_lib kilca_lib lapack cerf fortnum_amos_compat) +add_test(NAME test_periodize + COMMAND ${CMAKE_BINARY_DIR}/tests/test_periodize.x) diff --git a/KIM/tests/kernel_test_background.f90 b/KIM/tests/kernel_test_background.f90 new file mode 100644 index 00000000..c12ca2c9 --- /dev/null +++ b/KIM/tests/kernel_test_background.f90 @@ -0,0 +1,105 @@ +module kernel_test_background_m + ! Uniform two-species background written directly into species_m state, + ! so kernel-test expectations derive from the same first-principles + ! inputs without running the profile pipeline. + + use KIM_kinds_m, only: dp + use constants_m, only: pi, ev, e_charge, sol, com_unit + use species_m, only: plasma, init_electron_species, init_hydrogen_species + + implicit none + + integer, parameter :: npts = 100 + real(dp), parameter :: r_inner = 3.0d0 + real(dp), parameter :: r_outer = 70.0d0 + real(dp), parameter :: n0 = 2.0d13 + real(dp), parameter :: T0 = 1.0d3 + real(dp), parameter :: b0 = 2.0d4 + real(dp), parameter :: ks0 = 5.0d-2 + real(dp), parameter :: kp0 = 1.0d-3 + real(dp), parameter :: nu0 = 1.0d4 + +contains + + subroutine setup_uniform_background() + + implicit none + + integer :: sp, i + + if (.not. allocated(plasma%spec)) then + plasma%n_species = 2 + allocate (plasma%spec(0:1)) + call init_electron_species(plasma%spec(0)) + call init_hydrogen_species(plasma%spec(1)) + end if + + plasma%grid_size = npts + if (.not. allocated(plasma%r_grid)) then + allocate (plasma%r_grid(npts), plasma%ks(npts), plasma%kp(npts), & + plasma%om_E(npts)) + end if + do i = 1, npts + plasma%r_grid(i) = r_inner + (r_outer - r_inner)*(i - 1)/(npts - 1) + end do + plasma%ks = ks0 + plasma%kp = kp0 + plasma%om_E = 0.0d0 + + do sp = 0, plasma%n_species - 1 + associate (spec => plasma%spec(sp)) + if (.not. allocated(spec%vT)) then + allocate (spec%vT(npts), spec%omega_c(npts), & + spec%lambda_D(npts), spec%nu(npts), & + spec%z0(npts), spec%A1(npts), spec%A2(npts)) + end if + spec%vT = sqrt(T0*ev/spec%mass) + spec%omega_c = spec%Zspec*e_charge*b0/(spec%mass*sol) + spec%lambda_D = sqrt(T0*ev/(4.0d0*pi*n0 & + *(spec%Zspec*e_charge)**2)) + spec%nu = nu0 + spec%z0 = -(plasma%om_E - com_unit*nu0) & + /(abs(kp0)*sqrt(2.0d0)*spec%vT) + spec%A1 = 0.0d0 + spec%A2 = 0.0d0 + end associate + end do + + end subroutine setup_uniform_background + + subroutine set_forces(sp, a1_value, a2_value) + + implicit none + + integer, intent(in) :: sp + real(dp), intent(in) :: a1_value, a2_value + + plasma%spec(sp)%A1 = a1_value + plasma%spec(sp)%A2 = a2_value + + end subroutine set_forces + + subroutine require_close(name, got, want, rel_tol) + + implicit none + + character(*), intent(in) :: name + complex(dp), intent(in) :: got, want + real(dp), intent(in) :: rel_tol + + real(dp) :: scale + + scale = max(abs(want), 1.0d0) + if (abs(got - want) <= rel_tol*scale) then + print *, "PASS ", name + else + print *, "FAIL ", name + print *, " got = ", got + print *, " want = ", want + print *, " rel = ", abs(got - want)/scale + error stop + end if + + end subroutine require_close + +end module kernel_test_background_m diff --git a/KIM/tests/test_kernel_rho_B.f90 b/KIM/tests/test_kernel_rho_B.f90 new file mode 100644 index 00000000..d7bcf0e9 --- /dev/null +++ b/KIM/tests/test_kernel_rho_B.f90 @@ -0,0 +1,62 @@ +program test_kernel_rho_B + ! Behavioral checks of the restored Fourier-space rho-B kernel: + ! exact zero without thermodynamic forces, linearity in A1 and A2, + ! and the swap phase relation. + + use KIM_kinds_m, only: dp + use constants_m, only: com_unit + use species_m, only: plasma + use kernels_m, only: kernel_rho_B_of_kr_krp_rg + use kernel_test_background_m, only: setup_uniform_background, & + set_forces, require_close + + implicit none + + real(dp) :: kr, krp, rg + complex(dp) :: kernel_value, with_a1, with_2a1, swapped + + call setup_uniform_background() + rg = plasma%r_grid(50) + kr = 1.0d0 + krp = 3.0d0 + + ! Without thermodynamic forces the magnetic drive vanishes exactly. + kernel_value = kernel_rho_B_of_kr_krp_rg(kr, krp, rg) + if (abs(kernel_value) > 0.0d0) then + print *, "FAIL zero-force kernel is not exactly zero: ", kernel_value + error stop + end if + print *, "PASS zero thermodynamic forces give exactly zero" + + ! Linear in A1 at fixed background. + call set_forces(1, 1.0d-2, 0.0d0) + with_a1 = kernel_rho_B_of_kr_krp_rg(kr, krp, rg) + if (abs(with_a1) <= 0.0d0) then + print *, "FAIL A1 force does not drive the kernel" + error stop + end if + call set_forces(1, 2.0d-2, 0.0d0) + with_2a1 = kernel_rho_B_of_kr_krp_rg(kr, krp, rg) + call require_close("kernel is linear in A1", with_2a1, 2.0d0*with_a1, & + 1.0d-12) + + ! Linear in A2 at fixed background. + call set_forces(1, 0.0d0, 1.0d-2) + with_a1 = kernel_rho_B_of_kr_krp_rg(kr, krp, rg) + if (abs(with_a1) <= 0.0d0) then + print *, "FAIL A2 force does not drive the kernel" + error stop + end if + call set_forces(1, 0.0d0, 2.0d-2) + with_2a1 = kernel_rho_B_of_kr_krp_rg(kr, krp, rg) + call require_close("kernel is linear in A2", with_2a1, 2.0d0*with_a1, & + 1.0d-12) + + ! kr <-> krp swap only reverses the Fourier phase. + kernel_value = kernel_rho_B_of_kr_krp_rg(kr, krp, rg) + swapped = kernel_rho_B_of_kr_krp_rg(krp, kr, rg) + call require_close("kr <-> krp swap reverses the phase", swapped, & + kernel_value*exp(-2.0d0*com_unit*(kr - krp)*rg), & + 1.0d-12) + +end program test_kernel_rho_B diff --git a/KIM/tests/test_kernel_rho_phi.f90 b/KIM/tests/test_kernel_rho_phi.f90 new file mode 100644 index 00000000..56e3bbaf --- /dev/null +++ b/KIM/tests/test_kernel_rho_phi.f90 @@ -0,0 +1,90 @@ +program test_kernel_rho_phi + ! Behavioral checks of the restored Fourier-space rho-phi kernel: + ! adiabatic Debye screening, gyroaverage Gaussian, swap conjugacy, + ! phase translation, Debye switch, and linearity in A1 and A2. + + use KIM_kinds_m, only: dp + use constants_m, only: pi, com_unit + use species_m, only: plasma + use kernels_m, only: kernel_rho_phi_of_kr_krp_rg, kernel_debye_case + use kernel_test_background_m, only: setup_uniform_background, & + set_forces, require_close + + implicit none + + real(dp) :: kr, krp, rg, shift, debye_sum, gauss_sum, rho_l + complex(dp) :: kernel_value, swapped, shifted, base, with_a1, with_2a1 + integer :: sp + + call setup_uniform_background() + rg = plasma%r_grid(50) + + ! Adiabatic screening limit: uniform plasma, zero forces, kr = krp. + kr = 1.0d0 + krp = 1.0d0 + debye_sum = 0.0d0 + do sp = 0, plasma%n_species - 1 + debye_sum = debye_sum + 1.0d0/plasma%spec(sp)%lambda_D(1)**2 + end do + kernel_value = kernel_rho_phi_of_kr_krp_rg(kr, krp, rg) + call require_close("adiabatic Debye screening -1/(8 pi^2 lambda_D^2)", & + kernel_value, cmplx(-debye_sum/(8.0d0*pi**2), 0.0d0, dp), & + 1.0d-12) + + ! Gyroaverage Gaussian and Fourier phase at kr /= krp. + krp = 3.0d0 + gauss_sum = 0.0d0 + do sp = 0, plasma%n_species - 1 + rho_l = plasma%spec(sp)%vT(1)/abs(plasma%spec(sp)%omega_c(1)) + gauss_sum = gauss_sum + exp(-rho_l**2/2.0d0*(kr - krp)**2) & + /plasma%spec(sp)%lambda_D(1)**2 + end do + kernel_value = kernel_rho_phi_of_kr_krp_rg(kr, krp, rg) + call require_close("gyroaverage Gaussian with exp(i (kr-krp) rg) phase", & + kernel_value, & + -gauss_sum/(8.0d0*pi**2)*exp(com_unit*(kr - krp)*rg), & + 1.0d-12) + + ! Swapping kr and krp conjugates the uniform kernel. + swapped = kernel_rho_phi_of_kr_krp_rg(krp, kr, rg) + call require_close("kr <-> krp swap conjugates the uniform kernel", & + swapped, conjg(kernel_value), 1.0d-12) + + ! Translating rg only rotates the phase. + shift = 0.7d0 + shifted = kernel_rho_phi_of_kr_krp_rg(kr, krp, rg + shift) + call require_close("rg translation rotates phase by exp(i (kr-krp) dr)", & + shifted, kernel_value*exp(com_unit*(kr - krp)*shift), & + 1.0d-12) + + ! Debye switch: drops the force response, keeps the adiabatic part. + base = kernel_rho_phi_of_kr_krp_rg(kr, krp, rg) + call set_forces(1, 3.0d-2, 2.0d-2) + kernel_debye_case = .true. + kernel_value = kernel_rho_phi_of_kr_krp_rg(kr, krp, rg) + call require_close("Debye-case switch removes the force response", & + kernel_value, base, 1.0d-12) + kernel_debye_case = .false. + call set_forces(1, 0.0d0, 0.0d0) + + ! Kinetic response is linear in each thermodynamic force. + base = kernel_rho_phi_of_kr_krp_rg(kr, krp, rg) + call set_forces(1, 1.0d-2, 0.0d0) + with_a1 = kernel_rho_phi_of_kr_krp_rg(kr, krp, rg) + call set_forces(1, 2.0d-2, 0.0d0) + with_2a1 = kernel_rho_phi_of_kr_krp_rg(kr, krp, rg) + if (abs(with_a1 - base) <= 0.0d0) then + print *, "FAIL A1 force does not change the kernel" + error stop + end if + call require_close("kernel is linear in A1", with_2a1 - base, & + 2.0d0*(with_a1 - base), 1.0d-12) + + call set_forces(1, 0.0d0, 1.0d-2) + with_a1 = kernel_rho_phi_of_kr_krp_rg(kr, krp, rg) + call set_forces(1, 0.0d0, 2.0d-2) + with_2a1 = kernel_rho_phi_of_kr_krp_rg(kr, krp, rg) + call require_close("kernel is linear in A2", with_2a1 - base, & + 2.0d0*(with_a1 - base), 1.0d-12) + +end program test_kernel_rho_phi diff --git a/KIM/tests/test_periodize.f90 b/KIM/tests/test_periodize.f90 new file mode 100644 index 00000000..63394923 --- /dev/null +++ b/KIM/tests/test_periodize.f90 @@ -0,0 +1,114 @@ +program test_periodize + ! Behavioral checks of the radial periodization: untouched resonant + ! layer, exact L-periodicity, pass-through of an already-periodic + ! profile, and seam smoothness. + + use KIM_kinds_m, only: dp + use constants_m, only: pi + use periodize_m, only: periodized_profile_value + + implicit none + + integer, parameter :: npts = 381 + real(dp), parameter :: r_lo = -10.0d0 + real(dp), parameter :: r_hi = 85.0d0 + real(dp), parameter :: r_mid = 36.5d0 + real(dp), parameter :: dr_layer = 5.0d0 + real(dp), parameter :: dr_transition = 10.0d0 + real(dp), parameter :: period = 2.0d0*(dr_layer + dr_transition) + + real(dp) :: r_grid(npts), quadratic(npts), harmonic(npts) + real(dp) :: x, got, want, slope_left, slope_right, h + integer :: i + + do i = 1, npts + r_grid(i) = r_lo + (r_hi - r_lo)*(i - 1)/(npts - 1) + end do + quadratic = r_grid**2 + harmonic = sin(2.0d0*pi*(r_grid - r_mid)/period) + + ! Inside |x - r_mid| <= dr_layer the profile is unmodified; + ! four-point Lagrange interpolation is exact for a quadratic. + do i = 0, 10 + x = r_mid - dr_layer + 2.0d0*dr_layer*i/10.0d0 + got = periodized_profile_value(r_grid, quadratic, r_mid, dr_layer, & + dr_transition, x) + call require_close("resonant layer is unmodified", got, x**2, 1.0d-12) + end do + + ! Exact L-periodicity across several periods. + do i = 0, 12 + x = r_mid - 2.5d0*period + 5.0d0*period*i/12.0d0 + got = periodized_profile_value(r_grid, quadratic, r_mid, dr_layer, & + dr_transition, x + period) + want = periodized_profile_value(r_grid, quadratic, r_mid, dr_layer, & + dr_transition, x) + call require_close("construction is L-periodic", got, want, 1.0d-12) + end do + + ! A profile that is already L-periodic passes through unchanged. + do i = 0, 12 + x = r_mid - 1.4d0*period + 2.8d0*period*i/12.0d0 + got = periodized_profile_value(r_grid, harmonic, r_mid, dr_layer, & + dr_transition, x) + want = sin(2.0d0*pi*(x - r_mid)/period) + call require_close("periodic profile passes through", got, want, & + 1.0d-8) + end do + + ! The seam is smooth: one-sided slopes agree to the finite-difference + ! error, far below the O(1) jump an unblended wrap would give. + h = 1.0d-3 + x = r_mid - 0.5d0*period + slope_left = (periodized_profile_value(r_grid, quadratic, r_mid, & + dr_layer, dr_transition, x - h) & + - periodized_profile_value(r_grid, quadratic, r_mid, & + dr_layer, dr_transition, x - 3.0d0*h)) & + /(2.0d0*h) + slope_right = (periodized_profile_value(r_grid, quadratic, r_mid, & + dr_layer, dr_transition, x + 3.0d0*h) & + - periodized_profile_value(r_grid, quadratic, r_mid, & + dr_layer, dr_transition, x + h)) & + /(2.0d0*h) + call require_close("seam slopes agree", slope_left, slope_right, 1.0d-3) + + ! The layer edge starts the transition zone: a localizer without + ! vanishing edge derivatives would kink the profile here. + x = r_mid + dr_layer + slope_left = (periodized_profile_value(r_grid, quadratic, r_mid, & + dr_layer, dr_transition, x - h) & + - periodized_profile_value(r_grid, quadratic, r_mid, & + dr_layer, dr_transition, x - 3.0d0*h)) & + /(2.0d0*h) + slope_right = (periodized_profile_value(r_grid, quadratic, r_mid, & + dr_layer, dr_transition, x + 3.0d0*h) & + - periodized_profile_value(r_grid, quadratic, r_mid, & + dr_layer, dr_transition, x + h)) & + /(2.0d0*h) + call require_close("transition edge slopes agree", slope_left, & + slope_right, 1.0d-3) + + print *, "PASS all periodize checks" + +contains + + subroutine require_close(name, got, want, rel_tol) + + implicit none + + character(*), intent(in) :: name + real(dp), intent(in) :: got, want, rel_tol + + real(dp) :: scale + + scale = max(abs(want), 1.0d0) + if (abs(got - want) > rel_tol*scale) then + print *, "FAIL ", name + print *, " got = ", got + print *, " want = ", want + error stop + end if + + end subroutine require_close + +end program test_periodize