Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 8 additions & 1 deletion KIM/src/background_equilibrium/species_mod.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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))
Expand All @@ -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)
Expand Down
73 changes: 68 additions & 5 deletions KIM/tests/test_periodized_background_feed.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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))
Expand All @@ -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).
Expand Down Expand Up @@ -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
Expand Down