From d0e35722682a5cb9c03becbd571d3ad12a9c3a4a Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Mon, 13 Jul 2026 02:03:29 +0200 Subject: [PATCH] Recompute Krook argument after interpolation --- .../background_equilibrium/species_mod.f90 | 9 ++- KIM/tests/test_periodized_background_feed.f90 | 73 +++++++++++++++++-- 2 files changed, 76 insertions(+), 6 deletions(-) diff --git a/KIM/src/background_equilibrium/species_mod.f90 b/KIM/src/background_equilibrium/species_mod.f90 index 53f59122..f3bd5011 100644 --- a/KIM/src/background_equilibrium/species_mod.f90 +++ b/KIM/src/background_equilibrium/species_mod.f90 @@ -739,7 +739,9 @@ subroutine write_species_cc_quantities(spec, r_grid_cc) subroutine interpolate_plasma_backs(plasma_in, grid) use KIM_kinds_m, only: dp + use constants_m, only: com_unit use IO_collection_m, only: plot_profile + use setup_m, only: omega implicit none @@ -784,7 +786,6 @@ subroutine interpolate_plasma_backs(plasma_in, grid) plasma_temp%spec(sp)%omega_c(i) = sum(coef(0,:) * plasma_in%spec(sp)%omega_c(ibeg:iend)) plasma_temp%spec(sp)%lambda_D(i) = sum(coef(0,:) * plasma_in%spec(sp)%lambda_D(ibeg:iend)) plasma_temp%spec(sp)%rho_L(i) = sum(coef(0,:) * plasma_in%spec(sp)%rho_L(ibeg:iend)) - plasma_temp%spec(sp)%z0(i) = sum(coef(0,:) * plasma_in%spec(sp)%z0(ibeg:iend)) end do plasma_temp%B0(i) = sum(coef(0,:) * plasma_in%B0(ibeg:iend)) @@ -797,6 +798,12 @@ subroutine interpolate_plasma_backs(plasma_in, grid) end do + do sp = 0, plasma_temp%n_species-1 + plasma_temp%spec(sp)%z0 = -(plasma_temp%om_E - omega & + - com_unit*plasma_temp%spec(sp)%nu) & + /(abs(plasma_temp%kp)*sqrt(2.0d0)*plasma_temp%spec(sp)%vT) + end do + plasma_in = plasma_temp deallocate(plasma_temp) diff --git a/KIM/tests/test_periodized_background_feed.f90 b/KIM/tests/test_periodized_background_feed.f90 index 6aa88812..131913b5 100644 --- a/KIM/tests/test_periodized_background_feed.f90 +++ b/KIM/tests/test_periodized_background_feed.f90 @@ -5,11 +5,13 @@ program test_periodized_background_feed ! profile rather than periodizing the physical gradient. use KIM_kinds_m, only: dp - use constants_m, only: e_charge, ev - use species_m, only: plasma + use constants_m, only: e_charge, ev, com_unit + use species_m, only: plasma, plasma_t, init_electron_species, & + interpolate_plasma_backs use kernels_m, only: kernel_rho_phi_of_kr_krp_rg use rt_fourier_periodic_m, only: build_periodized_background - use kernel_test_background_m, only: setup_uniform_background, n0, T0, b0 + use kernel_test_background_m, only: setup_uniform_background, n0, T0, b0, & + ks0, nu0 use config_m, only: number_of_ion_species, rescale_density, & ion_flr_scale_factor, collision_frequency_scale use setup_m, only: collisions_off, omega @@ -27,7 +29,7 @@ program test_periodized_background_feed real(dp) :: r, physical, got, centered, physical_slope, deviation real(dp) :: lee, lei, expected_nue - complex(dp) :: kernel_a, kernel_b + complex(dp) :: kernel_a, kernel_b, expected_z0 integer :: i, j, sp number_of_ion_species = 1 @@ -37,6 +39,8 @@ program test_periodized_background_feed collisions_off = .false. omega = 0.0d0 + call check_recomputed_z0 + call setup_uniform_background() allocate (plasma%Er(plasma%grid_size), plasma%B0(plasma%grid_size), & plasma%q(plasma%grid_size)) @@ -61,7 +65,15 @@ program test_periodized_background_feed + 7.7d-6*ni*lei/te**1.5d0) end associate call require_close("collision-frequency scale reaches the background", & - plasma%spec(0)%nu(j), expected_nue, 1.0d-12) + plasma%spec(0)%nu(j), expected_nue, 1.0d-12) + + j = point_index(r_res - 2.0d0) + expected_z0 = -(plasma%om_E(j) - omega & + - com_unit*plasma%spec(0)%nu(j)) & + /(abs(plasma%kp(j))*sqrt(2.0d0)*plasma%spec(0)%vT(j)) + if (abs(plasma%spec(0)%z0(j) - expected_z0) & + > 1.0d-12*max(abs(expected_z0), 1.0d0)) & + error stop 'interpolated z0 violates its defining formula' ! Inside the resonant layer the linear physical profile survives ! exactly (four-point Lagrange interpolation is exact on it). @@ -120,6 +132,57 @@ program test_periodized_background_feed contains + subroutine check_recomputed_z0 + + type(plasma_t) :: coarse + real(dp) :: target_grid(3) + complex(dp) :: want + integer :: k + + coarse%n_species = 1 + coarse%grid_size = 4 + allocate (coarse%spec(0:0)) + call init_electron_species(coarse%spec(0)) + allocate (coarse%r_grid(4), coarse%ks(4), coarse%kp(4), & + coarse%om_E(4), coarse%B0(4), coarse%q(4), & + coarse%dqdr(4), coarse%Er(4)) + allocate (coarse%spec(0)%n(4), coarse%spec(0)%dndr(4), & + coarse%spec(0)%T(4), coarse%spec(0)%dTdr(4), & + coarse%spec(0)%nu(4), coarse%spec(0)%vT(4), & + coarse%spec(0)%omega_c(4), coarse%spec(0)%lambda_D(4), & + coarse%spec(0)%rho_L(4), coarse%spec(0)%z0(4)) + coarse%r_grid = [10.0d0, 20.0d0, 30.0d0, 40.0d0] + coarse%kp = [3.0d-3, 1.0d-3, -1.0d-3, -3.0d-3] + coarse%om_E = [9.3d4, 8.0d4, 7.0d4, 6.0d4] + coarse%ks = ks0 + coarse%B0 = b0 + coarse%q = 1.5d0 + coarse%dqdr = 0.0d0 + coarse%Er = 0.0d0 + coarse%spec(0)%n = n0 + coarse%spec(0)%dndr = 0.0d0 + coarse%spec(0)%T = T0 + coarse%spec(0)%dTdr = 0.0d0 + coarse%spec(0)%nu = nu0 + coarse%spec(0)%vT = 2.0d9 + coarse%spec(0)%omega_c = 1.0d8 + coarse%spec(0)%lambda_D = 1.0d-2 + coarse%spec(0)%rho_L = 1.0d-1 + coarse%spec(0)%z0 = -(coarse%om_E - omega - com_unit*nu0) & + /(abs(coarse%kp)*sqrt(2.0d0)*coarse%spec(0)%vT) + + target_grid = [15.0d0, 25.0d0, 35.0d0] + call interpolate_plasma_backs(coarse, target_grid) + do k = 1, size(target_grid) + want = -(coarse%om_E(k) - omega - com_unit*coarse%spec(0)%nu(k)) & + /(abs(coarse%kp(k))*sqrt(2.0d0)*coarse%spec(0)%vT(k)) + if (abs(coarse%spec(0)%z0(k) - want) & + > 1.0d-12*max(abs(want), 1.0d0)) & + error stop 'interpolated z0 violates its defining formula' + end do + + end subroutine check_recomputed_z0 + integer function point_index(r_target) implicit none