istp/gray
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases 1 Wiki Activity

src/gray_equil.f90: make Gaussian beam notation consistent

Browse Source 
To keep the notation used in compute_initial_conds consistent with the
source linked (Q = K - iW for the complex curvature tensor), we need to
adds a minus sign in a few places where ∇S_I is computed.
This commit is contained in:
Michele Guerini Rocco 2024-11-08 15:39:00 +01:00
parent 2368a0ba33
commit fd175ad0da
Signed by: rnhmjoj
GPG Key ID: BFBAF4C975F76450
1 changed files with 13 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
src/gray_core.f90
View File

@ -730,7 +730,7 @@ contains
W = matmul(rotate(phi_w), matmul(W, rotate(-phi_w))) W = matmul(rotate(phi_w), matmul(W, rotate(-phi_w)))
! Combine into a single tensor ! Combine into a single tensor
Q = K + im*W Q = K - im*W
end block end block
! Compute the beam → tokamak change-of-basis matrix ! Compute the beam → tokamak change-of-basis matrix
@ -741,12 +741,12 @@ contains
! !
! S_I = Im(z + ½ r̅⋅Q(z)⋅r̅) ! S_I = Im(z + ½ r̅⋅Q(z)⋅r̅)
! = ½ r̅⋅ImQ(z)⋅r̅ ! = ½ r̅⋅ImQ(z)⋅r̅
! = ½ r̅⋅W⋅r̅ ! = -½ r̅⋅W⋅r̅
! = ½ 2/k₀ ξ̅⋅diag(1/w_ξ², 1/w_η²)⋅ξ̅ ! = -½ 2/k₀ ξ̅⋅diag(1/w_ξ², 1/w_η²)⋅ξ̅
! = 1/k₀ ξ̅⋅diag(1/w_ξ², 1/w_η²)⋅ξ̅ ! = -1/k₀ ξ̅⋅diag(1/w_ξ², 1/w_η²)⋅ξ̅
! !
! Using polar coordinates ξ̅ = ρ [w_ξcosα, w_ηsinα] this ! Using polar coordinates ξ̅ = ρ [w_ξcosα, w_ηsinα] this
! simplifies to S_I = ρ²/k₀. To distribue the rays uniformly ! simplifies to S_I = -ρ²/k₀. To distribue the rays uniformly
! on the curves of constant S_I (power flux) we define a grid ! on the curves of constant S_I (power flux) we define a grid
! such that the jk-th ray is traced starting from ! such that the jk-th ray is traced starting from
! !
@ -762,14 +762,14 @@ contains
! To estimate ∇S_I in `gradi_upd` we require the values of ! To estimate ∇S_I in `gradi_upd` we require the values of
! ∇u(ξ̅_jk), where u(x̅) ≡ ρ(x̅)/Δρ and ξ̅_jk is the ray ! ∇u(ξ̅_jk), where u(x̅) ≡ ρ(x̅)/Δρ and ξ̅_jk is the ray
! position on the grid. Since S_I = ρ²/k₀, we have ! position on the grid. Since S_I = -ρ²/k₀, we have
! !
! ∇S_I(ξ̅_jk) = (Δρ²/k₀) ∇u²(ξ̅_jk) ! ∇S_I(ξ̅_jk) = -(Δρ²/k₀) ∇u²(ξ̅_jk)
! = (Δρ²/k₀) 2u(ξ̅_jk) ∇u(ξ̅_jk) ! = -(Δρ²/k₀) 2u(ξ̅_jk) ∇u(ξ̅_jk)
! = (Δρ²/k₀) 2(ρ_j/Δρ) ∇u(ξ̅_jk) ! = -(Δρ²/k₀) 2(ρ_j/Δρ) ∇u(ξ̅_jk)
! = (Δρ²/k₀) 2(j-1) ∇u(ξ̅_jk) ! = -(Δρ²/k₀) 2(j-1) ∇u(ξ̅_jk)
! !
! So, ∇u = k₀/[2Δρ²(j-1)] ∇S_I, for j>1. ! So, ∇u = k₀/[2Δρ²(1-j)] ∇S_I, for j>1.
! For the central ray (j=1) r=0, so ∇S_I=HS_I=0 ! For the central ray (j=1) r=0, so ∇S_I=HS_I=0
jk = 1 jk = 1
@ -787,7 +787,7 @@ contains
xi = beam%w * dr * [cos((k-1)*da), sin((k-1)*da)] xi = beam%w * dr * [cos((k-1)*da), sin((k-1)*da)]
x0 = [matmul(rotate(phi_w), xi), zero] x0 = [matmul(rotate(phi_w), xi), zero]
pos(:,k,1) = x0_c pos(:,k,1) = x0_c
grad_u(:,k,1) = [matmul(Q%im, x0(1:2) * k0/(2*dr**2)), zero] grad_u(:,k,1) = [-matmul(Q%im, x0(1:2) * k0/(2*dr**2)), zero]
end do end do
! Compute x̅, ∇S_I, ∇u for all the other rays ! Compute x̅, ∇S_I, ∇u for all the other rays
@ -803,7 +803,7 @@ contains
S=S, grad_S=grad_S, hess_S=hess_S) S=S, grad_S=grad_S, hess_S=hess_S)
! Compute ∇u using the formula above ! Compute ∇u using the formula above
grad_u(:,k,j) = matmul(M, k0/(2 * dr**2 * (j-1)) * grad_S%im) grad_u(:,k,j) = matmul(M, k0/(2 * dr**2 * (1-j)) * grad_S%im)
! Switch from local beam to tokamak basis ! Switch from local beam to tokamak basis
x0 = x0_c + matmul(M, x0) ! ray position x0 = x0_c + matmul(M, x0) ! ray position