From fd175ad0dacf448e3ea5eae7474cea8abb5815c2 Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Fri, 8 Nov 2024 15:39:00 +0100 Subject: [PATCH] src/gray_equil.f90: make Gaussian beam notation consistent MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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. --- src/gray_core.f90 | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/gray_core.f90 b/src/gray_core.f90 index 3badfde..db73a26 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -730,7 +730,7 @@ contains W = matmul(rotate(phi_w), matmul(W, rotate(-phi_w))) ! Combine into a single tensor - Q = K + im*W + Q = K - im*W end block ! Compute the beam → tokamak change-of-basis matrix @@ -741,12 +741,12 @@ contains ! ! S_I = Im(z + ½ r̅⋅Q(z)⋅r̅) ! = ½ r̅⋅ImQ(z)⋅r̅ - ! = ½ r̅⋅W⋅r̅ - ! = ½ 2/k₀ ξ̅⋅diag(1/w_ξ², 1/w_η²)⋅ξ̅ - ! = 1/k₀ ξ̅⋅diag(1/w_ξ², 1/w_η²)⋅ξ̅ + ! = -½ r̅⋅W⋅r̅ + ! = -½ 2/k₀ ξ̅⋅diag(1/w_ξ², 1/w_η²)⋅ξ̅ + ! = -1/k₀ ξ̅⋅diag(1/w_ξ², 1/w_η²)⋅ξ̅ ! ! 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 ! 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 ! ∇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) - ! = (Δρ²/k₀) 2u(ξ̅_jk) ∇u(ξ̅_jk) - ! = (Δρ²/k₀) 2(ρ_j/Δρ) ∇u(ξ̅_jk) - ! = (Δρ²/k₀) 2(j-1) ∇u(ξ̅_jk) + ! ∇S_I(ξ̅_jk) = -(Δρ²/k₀) ∇u²(ξ̅_jk) + ! = -(Δρ²/k₀) 2u(ξ̅_jk) ∇u(ξ̅_jk) + ! = -(Δρ²/k₀) 2(ρ_j/Δρ) ∇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 jk = 1 @@ -787,7 +787,7 @@ contains xi = beam%w * dr * [cos((k-1)*da), sin((k-1)*da)] x0 = [matmul(rotate(phi_w), xi), zero] 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 ! 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) ! 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 x0 = x0_c + matmul(M, x0) ! ray position