From 86d5b5a6724944b3ce0d5e2ba99e7ada92ed456b Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Thu, 19 Sep 2024 11:51:36 +0200 Subject: [PATCH] src/gray_core: refactor ic_gb --- src/beamdata.f90 | 72 +++++++ src/beams.f90 | 76 +++++++ src/gray_core.f90 | 514 +++++++++++++++++--------------------------- src/gray_tables.f90 | 44 +--- src/utils.f90 | 28 +++ 5 files changed, 377 insertions(+), 357 deletions(-) diff --git a/src/beamdata.f90 b/src/beamdata.f90 index fdcaadf..4696bd6 100644 --- a/src/beamdata.f90 +++ b/src/beamdata.f90 @@ -65,4 +65,76 @@ contains end if end function unfold_index + + pure function fold_indices(params, j, k) + ! Maps the pair (j, k) of radial (j) and azimuthal (k) + ! indices to a single index jk=1,nray. + use gray_params, only : raytracing_parameters + + type(raytracing_parameters), intent(in) :: params + integer, intent(in) :: j, k + integer :: fold_indices + + fold_indices = (j - 2) * params%nrayth + k + 1 ! jk + end function fold_indices + + + pure function tokamak2beam(n, inverse) result(M) + ! Returns the change-of-basis matrix for switching + ! from the standard tokamak and the local beam basis. + use const_and_precisions, only : zero, one + use utils, only : diag + + ! function arguments + real(wp_), intent(in) :: n(3) ! central ray direction + logical, intent(in), optional :: inverse ! do the opposite transformation + real(wp_) :: M(3,3) ! change-of-basis matrix + + ! local variables + real(wp_) :: c + logical :: inverse_ + + inverse_ = .false. + if (present(inverse)) inverse_ = inverse + + ! If n̅ is the normalised direction of the central ray, + ! and the {e̅₁,e̅₂,e̅₃} the standard basis, the unit vectors of + ! the beam basis are given by: + ! + ! x̅ = n̅×e̅₃ / |n̅×e̅₃| = (n₂/c)e̅₁ + (-n₁/c)e̅₂ + ! y̅ = n̅×(n̅×e̅₃) / |n̅×e̅₃| = (n₁n₃/c)e̅₁ + (n₂n₃/c)e̅₂ - ce̅₃ + ! z̅ = n̅ = n₁e̅₁ + n₂e̅₂ + n₃e̅₃ + ! + ! where c = |n̅×e̅₃| = √(n₁² + n₂²) + ! + ! Or in the case where n̅ ∥ e̅₃ by: + ! + ! x̅ = e̅₁ + ! y̅ = (n₃e̅₃)×e̅₁ = n₃e̅₂ + ! z̅ = n₃e̅₃ + ! + ! So, the change of basis matrix is + ! + ! [n₂/c -n₁/c 0] + ! M = [n₁n₃/c n₂n₃/c -c], c > 0 + ! [n₁ n₂ n₃] + ! + ! or M = diag(1, n₃, n₃) if c = 0. + ! + ! By definition then we have: x̅_loc = M x̅. + ! + c = norm2(n(1:2)) + if (c > 0) then + M = reshape([ n(2)/c, n(1)*n(3)/c, n(1), & + -n(1)/c, n(2)*n(3)/c, n(2), & + zero, -c, n(3) ], [3, 3]) + else + M = diag([one, n(3), n(3)]) + end if + + ! The basis is orthonormal, so M⁻¹ is just the transpose + if (inverse_) M = transpose(M) + + end function tokamak2beam + end module beamdata diff --git a/src/beams.f90 b/src/beams.f90 index 290ac0e..f91a212 100644 --- a/src/beams.f90 +++ b/src/beams.f90 @@ -862,4 +862,80 @@ contains xgcn = 4.0e13_wp_ * pi * qe**2/(me * omega**2) ! [10^-19 m^3] end subroutine xgygcoeff + + pure subroutine gaussian_eikonal(Q, r, z, S, grad_S, hess_S) + ! Computes the complex eikonal, gradient and Hessian of a Gaussian beam + + ! subroutine arguments + complex(wp_), intent(in) :: Q(2,2) ! complex curvature tensor Q(z) + real(wp_), intent(in) :: r(2), z ! coordinates + complex(wp_), intent(out) :: S ! complex eikonal + complex(wp_), intent(out) :: grad_S(3) ! gradient of imaginary eikonal + complex(wp_), intent(out) :: hess_S(3,3) ! Hessian of imaginary eikonal + + ! local variables + complex(wp_) :: Q2(2,2), Q3(2,2) + Q2 = matmul(Q, Q) + Q3 = matmul(Q, Q2) + + ! The most general Gaussian beam is given by + ! + ! E̅(r̅,z) = e̅(z) exp[-ik₀z - ik₀r̅⋅Q(z)⋅r̅/2 + iη(z)] + ! + ! where: + ! - e̅(z) is the amplitude + ! - Q(z) = Q₀ (I + zQ₀)⁻¹ is the complex curvature tensor + ! - η(z) is the Gouy phase + ! - k₀ = ω/c + ! - r̅ = [x, y] + ! + ! See https://doi.org/10.1364/AO.52.006030 for more details. + ! + ! Since the eikonal ansatz in GRAY is E̅(x̅,t) = e̅(r̅) exp[-ik₀S(x̅) + iωt], + ! the eikonal for a Gaussian beam is: + ! + ! S(x̅) = z + ½ r̅⋅Q(z)⋅r̅ - η(z)/k₀ + ! + S = z + dot_product(r, matmul(Q, r))/2 + + ! Splitting up x̅ into (r̅, z), the gradient of S can be written as + ! + ! ∇S(r̅,z) = [Q(z)⋅r̅, 1 + ½ r̅⋅dQ/dz⋅r̅] + ! + ! For dQ/dz, we use the matrix identity dA/dz = - A⁻¹ dA/dz A⁻¹: + ! + ! dQ/dz = Q₀ d/dz (I + zQ₀)⁻¹ + ! = - Q₀ (I + zQ₀)⁻¹ Q₀ (I + zQ₀)⁻¹ + ! = - Q(z)² + ! + ! So ∇S(r̅,z) = [Q(z)⋅r̅, 1 - ½ r̅⋅Q(z)²⋅r̅] + ! + grad_S = [matmul(Q, r), 1 - dot_product(r, matmul(Q2, r))/2] + + ! Carrying on, the Hessian splits into 4 blocks (2x2, 2x1, 1x2, 1x1): + ! + ! HS(r̅,z) = k₀ [ Q(z), dQ/dz⋅r̅] = k₀ [ Q(z), -Q(z)²⋅r̅] + ! [-Q(z)²⋅r̅, -½ r̅⋅dQ²/dz⋅r̅] [-Q(z)²⋅r̅, r̅⋅Q³⋅r̅] + ! + ! having used dQ²/dz = 2Q dQ/dz = -2Q³. + ! + hess_S(1:2,:) = reshape([Q, -matmul(Q2, r)], [2,3]) + hess_S(3, :) = [-matmul(Q2, r), dot_product(r, matmul(Q3, r))] + + ! Note: this ignores the Gouy phase, but for the computation of the + ! initial conditions (x̅, N̅_R=∇S_R, N̅_I=∇S_I for each ray) it is + ! inconsequential because: + ! + ! 1. the initial conditions are given at z=const, so η=const as well + ! 2. N̅_I = ∇S_I does not depend on η(z) + ! 3. N̅_R does depend on η(z), but is obtained indirectly as the + ! solution of the vacuum quasioptical dispersion relation: + ! + ! { N̅_R ⋅ N̅_I = 0 + ! { N_R² - N_I² = N₀² + ! + ! The latter implies N̅_I and N̅_R are not exactly consistent with + ! a Gaussian beam, but this is a necessary trade-off. + end subroutine gaussian_eikonal + end module beams diff --git a/src/gray_core.f90 b/src/gray_core.f90 index a13c574..ba85c87 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -220,10 +220,20 @@ contains lgcpl1 = zero p0ray = p0jk ! * initial beam power - call ic_gb(params, equil, anv0, ak0, yw, ypw, stv, xc, & - du1, gri, ggri, index_rt, results%tables) ! * initial conditions + call compute_initial_conds(params%raytracing, params%antenna, & ! * initial conditions + anv0, ak0, yw, ypw, stv, xc, du1, gri, ggri) do jk=1,params%raytracing%nray + ! Save the step "zero" data + call store_ray_data(params, equil, results%tables, & + i=0, jk=jk, s=stv(jk), P0=one, pos=yw(1:3,jk), & + psi_n=-one, B=zero, b_n=[zero,zero,zero], k0=ak0, & + N_pl=zero, N_pr=zero, N=yw(4:6,jk), N_pr_im=zero, & + n_e=zero, T_e=zero, alpha=zero, tau=zero, dIds=zero, & + nhm=0, nhf=0, iokhawa=0, index_rt=index_rt, & + lambda_r=zero, lambda_i=zero, X=zero, Y=zero, & + grad_X=[zero,zero,zero]) + zzm = yw(3,jk)*0.01_wp_ rrm = sqrt(yw(1,jk)*yw(1,jk)+yw(2,jk)*yw(2,jk))*0.01_wp_ @@ -670,339 +680,209 @@ contains end subroutine vectinit - subroutine ic_gb(params, equil, anv0c, ak0, ywrk0, ypwrk0, & - stv, xc0, du10, gri, ggri, index_rt, & - tables) - ! beam tracing initial conditions igrad=1 - ! !!!!!! check ray tracing initial conditions igrad=0 !!!!!! - use const_and_precisions, only : zero,one,pi,half,two,degree,ui=>im - use gray_params, only : gray_parameters, gray_tables - use gray_equil, only : abstract_equil - use types, only : table - use gray_tables, only : store_ray_data + subroutine compute_initial_conds(rtx, beam, N_c, k0, y, yp, dist, & + pos, grad_u, grad, hess) + ! Computes the initial conditions for tracing a beam + + use const_and_precisions, only : zero, one, pi, im, degree + use gray_params, only : raytracing_parameters, antenna_parameters + use gray_params, only : STEP_ARCLEN, STEP_TIME, STEP_PHASE + use utils, only : diag, rotate + use beams, only : gaussian_eikonal + use beamdata, only : fold_indices, tokamak2beam ! subroutine arguments - type(gray_parameters), intent(in) :: params - class(abstract_equil), intent(in) :: equil - real(wp_), dimension(3), intent(in) :: anv0c - real(wp_), intent(in) :: ak0 - real(wp_), dimension(6, params%raytracing%nray), intent(out) :: ywrk0, ypwrk0 - real(wp_), dimension(params%raytracing%nrayth * params%raytracing%nrayr), intent(out) :: stv - real(wp_), dimension(3, params%raytracing%nrayth, params%raytracing%nrayr), intent(out) :: xc0, du10 - real(wp_), dimension(3, params%raytracing%nray), intent(out) :: gri - real(wp_), dimension(3, 3, params%raytracing%nray), intent(out) :: ggri - integer, intent(in) :: index_rt - ! tables for storing initial rays data - type(gray_tables), intent(inout), optional :: tables + ! Inputs + type(raytracing_parameters), intent(in) :: rtx ! simulation parameters + type(antenna_parameters), intent(in) :: beam ! beam parameters + real(wp_), intent(in) :: N_c(3) ! vacuum refractive index + real(wp_), intent(in) :: k0 ! vacuum wavenumber + + ! Outputs + real(wp_), intent(out) :: y(6, rtx%nray) ! (x̅, N̅) for each ray + real(wp_), intent(out) :: yp(6, rtx%nray) ! (dx̅/dσ, dN̅/dσ) for each ray + real(wp_), intent(out) :: dist(rtx%nrayth * rtx%nrayr) ! distance travelled in σ for each ray + real(wp_), intent(out) :: pos(3, rtx%nrayth, rtx%nrayr) ! ray positions (used in gradi_upd) + real(wp_), intent(out) :: grad_u(3, rtx%nrayth, rtx%nrayr) ! variations Δu (used in gradi_upd) + real(wp_), intent(out) :: grad(3, rtx%nray) ! ∇S_I, gradient of imaginary eikonal + real(wp_), intent(out) :: hess(3, 3, rtx%nray) ! HS_I, Hessian of imaginary eikonal ! local variables - real(wp_), dimension(3) :: xv0c - real(wp_) :: wcsi,weta,rcicsi,rcieta,phiw,phir - integer :: j,k,jk - real(wp_) :: csth,snth,csps,snps,phiwrad,phirrad,csphiw,snphiw,alfak - real(wp_) :: wwcsi,wweta,sk,sw,dk,dw,rci1,ww1,rci2,ww2,wwxx,wwyy,wwxy - real(wp_) :: rcixx,rciyy,rcixy,dwwxx,dwwyy,dwwxy,d2wwxx,d2wwyy,d2wwxy - real(wp_) :: drcixx,drciyy,drcixy,dr,da,ddfu,dcsiw,detaw,dx0t,dy0t - real(wp_) :: x0t,y0t,z0t,dx0,dy0,dz0,x0,y0,z0,gxt,gyt,gzt,gr2 - real(wp_) :: gxxt,gyyt,gzzt,gxyt,gxzt,gyzt,dgr2xt,dgr2yt,dgr2zt - real(wp_) :: dgr2x,dgr2y,dgr2z,pppx,pppy,denpp,ppx,ppy - real(wp_) :: anzt,anxt,anyt,anx,any,anz,an20,an0 - real(wp_) :: du1tx,du1ty,du1tz,denom,ddr,ddi - real(wp_), dimension(params%raytracing%nrayr) :: uj - real(wp_), dimension(params%raytracing%nrayth) :: sna,csa - complex(wp_) :: sss,ddd,phic,qi1,qi2,tc,ts,qqxx,qqxy,qqyy,dqi1,dqi2 - complex(wp_) :: dqqxx,dqqyy,dqqxy,d2qi1,d2qi2,d2qqxx,d2qqyy,d2qqxy + real(wp_) :: x0_c(3), x0(3) ! beam centre, current ray position + real(wp_) :: xi(2) ! position in the principal axes basis (ξ, η) + integer :: j, k, jk ! ray indices, j: radial, k: angular + real(wp_) :: dr, da ! ray grid steps, dρ: radial, dα: angular + real(wp_) :: N(3) ! wavevector + complex(wp_) :: Q(2,2) ! complex cuvature tensor + real(wp_) :: M(3,3) ! local beam → tokamak change of basis matrix + complex(wp_) :: S ! complex eikonal + complex(wp_) :: grad_S(3) ! its gradient + complex(wp_) :: hess_S(3,3) ! its Hessian matrix + real(wp_) :: phi_k, phi_w ! ellipses rotation angles - ! initial beam parameters - xv0c = params%antenna%pos - wcsi = params%antenna%w(1) - weta = params%antenna%w(2) - rcicsi = params%antenna%ri(1) - rcieta = params%antenna%ri(2) - phiw = params%antenna%phi(1) - phir = params%antenna%phi(2) + ! Compute the complex curvature tensor Q from the beam parameters + block + real(wp_) :: K(2,2), W(2,2) + K = diag(beam%ri) ! curvature part + W = diag(2/k0 * 1/beam%w**2) ! beam width part - csth=anv0c(3) - snth=sqrt(one-csth**2) - if(snth > zero) then - csps=anv0c(2)/snth - snps=anv0c(1)/snth - else - csps=one - snps=zero - end if + ! Switch from eigenbasis to local beam basis + ! + ! Notes on geometry: + ! + ! 1. The beam parameters are the eigenvalues of the K, W matrices. + ! These are diagonal in the basis of the principal axes of the + ! respective ellipse (curvature, beam width). + ! + ! 2. The φs are nothing more that the angles of the rotations that + ! diagonalises the respective matrix. + ! + phi_k = beam%phi(2) * degree ! curvature rotation angle + phi_w = beam%phi(1) * degree ! beam width rotation angle - ! Gaussian beam: exp[-ik0 zt] exp[-i k0/2 S(xt,yt,zt)] - ! xt,yt,zt, cartesian coordinate system with zt along the beamline and xt in the z = 0 plane - ! S(xt,yt,zt) = S_real +i S_imag = Qxx(zt) xt^2 + Qyy(zt) yt^2 + 2 Qxy(zt) xt yt - ! (csiw, etaw) and (csiR, etaR) intensity and phase ellipse, rotated by angle phiw and phiR - ! S(xt,yt,zt) = csiR^2 / Rccsi +etaR^2 /Rceta - i (csiw^2 Wcsi +etaw^2 Weta) - ! Rccsi,eta curvature radius at the launching point - ! Wcsi,eta =2/(k0 wcsi,eta^2) with wcsi,eta^2 beam size at the launching point + K = matmul(rotate(phi_k), matmul(K, rotate(-phi_k))) + W = matmul(rotate(phi_w), matmul(W, rotate(-phi_w))) - phiwrad = phiw*degree - phirrad = phir*degree - csphiw = cos(phiwrad) - snphiw = sin(phiwrad) - !csphir = cos(phirrad) - !snphir = sin(phirrad) + ! Combine into a single tensor + Q = K + im*W + end block - wwcsi = two/(ak0*wcsi**2) - wweta = two/(ak0*weta**2) + ! Compute the beam → tokamak change-of-basis matrix + M = tokamak2beam(N_c, inverse=.true.) - if(phir/=phiw) then - sk = rcicsi + rcieta - sw = wwcsi + wweta - dk = rcicsi - rcieta - dw = wwcsi - wweta - ts = -(dk*sin(2*phirrad) - ui*dw*sin(2*phiwrad)) - tc = (dk*cos(2*phirrad) - ui*dw*cos(2*phiwrad)) - phic = half*atan(ts/tc) - ddd = dk*cos(2*(phirrad+phic)) - ui*dw*cos(2*(phiwrad+phic)) - sss = sk - ui*sw - qi1 = half*(sss + ddd) - qi2 = half*(sss - ddd) - rci1 = real(qi1) - rci2 = real(qi2) - ww1 = -imag(qi1) - ww2 = -imag(qi2) - else - rci1 = rcicsi - rci2 = rcieta - ww1 = wwcsi - ww2 = wweta - phic = -phiwrad - qi1 = rci1 - ui*ww1 - qi2 = rci2 - ui*ww2 - end if + ! Expanding the definitions (see `gaussian_eikonal`), the imaginary + ! part of the eikonal becomes: + ! + ! 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_η²)⋅ξ̅ + ! + ! Using polar coordinates ξ̅ = ρ [w_ξcosα, w_ηsinα] this + ! 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 + ! + ! ξ̅_jk = ρ_j [w_ξcos(α_k), w_ηcos(α_k)] + ! + ! where ρ_j = (j-1) Δρ, + ! α_k = (k-1) Δα, + ! Δρ = ρ_max/(N_r - 1), + ! Δα = 2π/N_θ. + ! + dr = merge(rtx%rwmax / (rtx%nrayr - 1), one, rtx%nrayr > 1) + da = 2*pi / rtx%nrayth - !w01=sqrt(2.0_wp_/(ak0*ww1)) - !d01=-rci1/(rci1**2+ww1**2) - !w02=sqrt(2.0_wp_/(ak0*ww2)) - !d02=-rci2/(rci2**2+ww2**2) + ! 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 + ! + ! ∇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. - qqxx = qi1*cos(phic)**2 + qi2*sin(phic)**2 - qqyy = qi1*sin(phic)**2 + qi2*cos(phic)**2 - qqxy = -(qi1 - qi2)*sin(phic)*cos(phic) - wwxx = -imag(qqxx) - wwyy = -imag(qqyy) - wwxy = -imag(qqxy) - rcixx = real(qqxx) - rciyy = real(qqyy) - rcixy = real(qqxy) + ! For the central ray (j=1) r=0, so ∇S_I=HS_I=0 + jk = 1 + grad(:,1) = 0 + hess(:,:,1) = 0 + x0_c = beam%pos + y(:,1) = [x0_c, N_c] + yp(:,1) = [N_c, 0*N_c] - dqi1 = -qi1**2 - dqi2 = -qi2**2 - d2qi1 = 2*qi1**3 - d2qi2 = 2*qi2**3 - dqqxx = dqi1*cos(phic)**2 + dqi2*sin(phic)**2 - dqqyy = dqi1*sin(phic)**2 + dqi2*cos(phic)**2 - dqqxy = -(dqi1 - dqi2)*sin(phic)*cos(phic) - d2qqxx = d2qi1*cos(phic)**2 + d2qi2*sin(phic)**2 - d2qqyy = d2qi1*sin(phic)**2 + d2qi2*cos(phic)**2 - d2qqxy = -(d2qi1 - d2qi2)*sin(phic)*cos(phic) - - dwwxx = -imag(dqqxx) - dwwyy = -imag(dqqyy) - dwwxy = -imag(dqqxy) - d2wwxx = -imag(d2qqxx) - d2wwyy = -imag(d2qqyy) - d2wwxy = -imag(d2qqxy) - drcixx = real(dqqxx) - drciyy = real(dqqyy) - drcixy = real(dqqxy) - - if(params%raytracing%nrayr > 1) then - dr = params%raytracing%rwmax / (params%raytracing%nrayr - 1) - else - dr = 1 - end if - ddfu = 2*dr**2/ak0 - do j = 1, params%raytracing%nrayr - uj(j) = real(j-1, wp_) + ! We consider the central ray as k degenerate rays + ! with positions x̅=x₀_c, but ∇u computed at the points + ! ξ̅_1k = Δρ [w_ξcos(α_k), w_ηcos(α_k)]. This choice + ! simplifies the computation of ∇u in `gradi_upd`. + do k = 1, rtx%nrayth + 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] end do - da=2*pi/params%raytracing%nrayth - do k=1,params%raytracing%nrayth - alfak = (k-1)*da - sna(k) = sin(alfak) - csa(k) = cos(alfak) + ! Compute x̅, ∇S_I, ∇u for all the other rays + do concurrent (j=2:rtx%nrayr, k=1:rtx%nrayth) + jk = fold_indices(rtx, j, k) + + ! Compute the ray position + xi = beam%w * (j-1)*dr * [cos((k-1)*da), sin((k-1)*da)] + x0 = [matmul(rotate(phi_w), xi), zero] + + ! Compute the Eikonal and its derivatives + call gaussian_eikonal(Q, r=x0(1:2), z=x0(3), & + 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) + + ! Switch from local beam to tokamak basis + x0 = x0_c + matmul(M, x0) ! ray position + grad(:,jk) = matmul(M, grad_S%im) ! gradient of S_I + hess(:,:,jk) = matmul(M, matmul(hess_S%im, transpose(M))) ! Hessian of S_I + + ! Compute the refractive index N such that: + ! + ! 1. N² = 1 + ∇S_I² + ! 2. N⋅∇S_I = 0 + ! 3. N₁/N₂ = ∇S_R₁/∇S_R₂ + ! + ! Note: 1. and 2. are the quasioptical dispersion relation, + ! 3. is necessary to obtain a unique solution. + block + real(wp_) :: den, r1, r2, gradS2 + den = dot_product(grad_S(1:2)%re, grad_S(1:2)%im) + gradS2 = norm2(grad_S%im)**2 + + if (den /= 0) then + r1 = -grad_S(1)%re * grad_S(3)%im / den + r2 = -grad_S(2)%re * grad_S(3)%im / den + N = [r1, r2, one] * sqrt((1 + gradS2)/(1 + r1**2 + r2**2)) + else if (grad_S(3)%im /= 0) then + ! In this case the solution reduces to + r1 = grad_S(1)%re + r2 = grad_S(2)%re + N = [r1, r2, zero] * sqrt((1 + gradS2)/(1 + r1**2 + r2**2)) + else + ! When even ∇S_I₃ = 0, the system is underdetermined: + ! we pick the solution N₃ = 1 + ∇S_I², so that N₁=N₂=0. + N = [zero, zero, sqrt(1 + gradS2)] + end if + end block + + ! Compute the actual initial conditions + pos(:,k,j) = x0 + y(:,jk) = [x0, matmul(M, N)] + + ! Compute the r.h.s. of the raytracing eqs. + ! dx̅/dσ = + ∂Λ/∂N̅ / denom + ! dN̅/dσ = - ∂Λ/∂x̅ / denom + block + real(wp_) :: denom + + ! Compute travelled distance in σ and r.h.s. normalisation + select case(rtx%idst) + case (STEP_ARCLEN) ! σ=s, arclength + dist(jk) = 0 + denom = norm2(N) + case (STEP_TIME) ! σ=ct, time + dist(jk) = 0 + denom = 1 + case (STEP_PHASE) ! σ=S_R, phase + dist(jk) = S%re + denom = norm2(N)**2 + end select + + yp(1:3,jk) = matmul(M, N) / denom + yp(4:6,jk) = matmul(M, matmul(hess_S%im, grad_S%im)) / denom + end block end do - ! central ray - jk=1 - gri(:,1) = zero - ggri(:,:,1) = zero - - ywrk0(1:3,1) = xv0c - ywrk0(4:6,1) = anv0c - ypwrk0(1:3,1) = anv0c - ypwrk0(4:6,1) = zero - - do k=1,params%raytracing%nrayth - dcsiw = dr*csa(k)*wcsi - detaw = dr*sna(k)*weta - dx0t = dcsiw*csphiw - detaw*snphiw - dy0t = dcsiw*snphiw + detaw*csphiw - du1tx = (dx0t*wwxx + dy0t*wwxy)/ddfu - du1ty = (dx0t*wwxy + dy0t*wwyy)/ddfu - - xc0(:,k,1) = xv0c - du10(1,k,1) = du1tx*csps + snps*du1ty*csth - du10(2,k,1) = -du1tx*snps + csps*du1ty*csth - du10(3,k,1) = -du1ty*snth - end do - ddr = zero - ddi = zero - - ! loop on rays jk>1 - j=2 - k=0 - do jk = 2, params%raytracing%nray - k=k+1 - if(k > params%raytracing%nrayth) then - j=j+1 - k=1 - end if - - !csiw=u*dcsiw - !etaw=u*detaw - !csir=x0t*csphir+y0t*snphir - !etar=-x0t*snphir+y0t*csphir - dcsiw = dr*csa(k)*wcsi - detaw = dr*sna(k)*weta - dx0t = dcsiw*csphiw - detaw*snphiw - dy0t = dcsiw*snphiw + detaw*csphiw - x0t = uj(j)*dx0t - y0t = uj(j)*dy0t - z0t = 0 - - dx0 = x0t*csps + snps*(y0t*csth + z0t*snth) - dy0 = -x0t*snps + csps*(y0t*csth + z0t*snth) - dz0 = z0t*csth - y0t*snth - x0 = xv0c(1) + dx0 - y0 = xv0c(2) + dy0 - z0 = xv0c(3) + dz0 - - gxt = x0t*wwxx + y0t*wwxy - gyt = x0t*wwxy + y0t*wwyy - gzt = half*(x0t**2*dwwxx + y0t**2*dwwyy ) + x0t*y0t*dwwxy - gr2 = gxt*gxt + gyt*gyt + gzt*gzt - gxxt = wwxx - gyyt = wwyy - gzzt = half*(x0t**2*d2wwxx + y0t**2*d2wwyy) + x0t*y0t*d2wwxy - gxyt = wwxy - gxzt = x0t*dwwxx + y0t*dwwxy - gyzt = x0t*dwwxy + y0t*dwwyy - dgr2xt = 2*(gxt*gxxt + gyt*gxyt + gzt*gxzt) - dgr2yt = 2*(gxt*gxyt + gyt*gyyt + gzt*gyzt) - dgr2zt = 2*(gxt*gxzt + gyt*gyzt + gzt*gzzt) - dgr2x = dgr2xt*csps + snps*(dgr2yt*csth + dgr2zt*snth) - dgr2y = -dgr2xt*snps + csps*(dgr2yt*csth + dgr2zt*snth) - dgr2z = dgr2zt*csth - dgr2yt*snth - - gri(1,jk) = gxt*csps + snps*(gyt*csth + gzt*snth) - gri(2,jk) = -gxt*snps + csps*(gyt*csth + gzt*snth) - gri(3,jk) = gzt*csth - gyt*snth - ggri(1,1,jk) = gxxt*csps**2 & - + snps**2 *(gyyt*csth**2 + gzzt*snth**2 + 2*snth*csth*gyzt) & - +2*snps*csps*(gxyt*csth + gxzt*snth) - ggri(2,1,jk) = csps*snps & - *(-gxxt+csth**2*gyyt + snth**2*gzzt + 2*csth*snth*gyzt) & - +(csps**2 - snps**2)*(snth*gxzt + csth*gxyt) - ggri(3,1,jk) = csth*snth*snps*(gzzt - gyyt) + (csth**2 - snth**2) & - *snps*gyzt + csps*(csth*gxzt - snth*gxyt) - ggri(1,2,jk) = ggri(2,1,jk) - ggri(2,2,jk) = gxxt*snps**2 & - + csps**2 *(gyyt*csth**2 + gzzt*snth**2 + 2*snth*csth*gyzt) & - -2*snps*csps*(gxyt*csth + gxzt*snth) - ggri(3,2,jk) = csth*snth*csps*(gzzt - gyyt) + (csth**2-snth**2) & - *csps*gyzt + snps*(snth*gxyt - csth*gxzt) - ggri(1,3,jk) = ggri(3,1,jk) - ggri(2,3,jk) = ggri(3,2,jk) - ggri(3,3,jk) = gzzt*csth**2 + gyyt*snth**2 - 2*csth*snth*gyzt - - du1tx = (dx0t*wwxx + dy0t*wwxy)/ddfu - du1ty = (dx0t*wwxy + dy0t*wwyy)/ddfu - du1tz = half*uj(j)*(dx0t**2*dwwxx + dy0t**2*dwwyy + 2*dx0t*dy0t*dwwxy)/ddfu - - du10(1,k,j) = du1tx*csps + snps*(du1ty*csth + du1tz*snth) - du10(2,k,j) = -du1tx*snps + csps*(du1ty*csth + du1tz*snth) - du10(3,k,j) = du1tz*csth - du1ty*snth - - pppx = x0t*rcixx + y0t*rcixy - pppy = x0t*rcixy + y0t*rciyy - denpp = pppx*gxt + pppy*gyt - if (denpp/=zero) then - ppx = -pppx*gzt/denpp - ppy = -pppy*gzt/denpp - else - ppx = zero - ppy = zero - end if - - anzt = sqrt((one + gr2)/(one + ppx**2 + ppy**2)) - anxt = ppx*anzt - anyt = ppy*anzt - - anx = anxt*csps + snps*(anyt*csth + anzt*snth) - any =-anxt*snps + csps*(anyt*csth + anzt*snth) - anz = anzt*csth - anyt*snth - - an20 = one + gr2 - an0 = sqrt(an20) - - xc0(1,k,j) = x0 - xc0(2,k,j) = y0 - xc0(3,k,j) = z0 - - ywrk0(1,jk) = x0 - ywrk0(2,jk) = y0 - ywrk0(3,jk) = z0 - ywrk0(4,jk) = anx - ywrk0(5,jk) = any - ywrk0(6,jk) = anz - - ! integration variable - select case(params%raytracing%idst) - case(0) - ! optical path: s - stv(jk) = 0 - denom = an0 - case(1) - ! "time": c⋅t - stv(jk) = 0 - denom = one - case(2) - ! real eikonal: S_R - stv(jk) = half*(rcixx*x0t**2 + rciyy*y0t**2) + rcixy*x0t*y0t - denom = an20 - end select - ypwrk0(1,jk) = anx/denom - ypwrk0(2,jk) = any/denom - ypwrk0(3,jk) = anz/denom - ypwrk0(4,jk) = dgr2x/(2*denom) - ypwrk0(5,jk) = dgr2y/(2*denom) - ypwrk0(6,jk) = dgr2z/(2*denom) - - ddr = anx**2 + any**2 + anz**2 - an20 - ddi = 2*(anxt*gxt + anyt*gyt + anzt*gzt) - - ! save step "zero" data - if (present(tables)) & - call store_ray_data(params, equil, tables, & - i=0, jk=jk, s=stv(jk), P0=one, pos=xc0(:,k,j), & - psi_n=-one, B=zero, b_n=[zero,zero,zero], k0=ak0, & - N_pl=zero, N_pr=zero, N=ywrk0(:, jk), N_pr_im=zero, & - n_e=zero, T_e=zero, alpha=zero, tau=zero, dIds=zero, & - nhm=0, nhf=0, iokhawa=0, index_rt=index_rt, & - lambda_r=ddr, lambda_i=ddi, X=zero, Y=zero, & - grad_X=[zero,zero,zero]) - - end do - end subroutine ic_gb - + end subroutine compute_initial_conds subroutine rkstep(params, equil, plasma, & diff --git a/src/gray_tables.f90 b/src/gray_tables.f90 index 11a1faa..11310ff 100644 --- a/src/gray_tables.f90 +++ b/src/gray_tables.f90 @@ -524,10 +524,10 @@ contains subroutine store_beam_shape(params, tables, s, y, full) ! Computes the beam shape tables - use const_and_precisions, only : zero, one, comp_huge + use const_and_precisions, only : comp_huge use types, only : table, wrap use gray_params, only : raytracing_parameters, gray_tables - use beamdata, only : unfold_index + use beamdata, only : unfold_index, tokamak2beam ! subroutine arguments type(raytracing_parameters), intent(in) :: params @@ -541,7 +541,7 @@ contains ! local variables logical :: full_ integer :: jk, jk0, jkv(2) - real(wp_) :: n(3), M(3,3), dx(3), dx_loc(3), c + real(wp_) :: n(3), M(3,3), dx(3), dx_loc(3) real(wp_) :: r, r_min, r_max type(table), pointer :: beam_shape @@ -557,44 +557,8 @@ contains end if ! Convert to the local beam basis - ! - ! If n̅ is the normalised direction of the central ray, - ! and the {e̅₁,e̅₂,e̅₃} the standard basis, the unit vectors of - ! the beam basis are given by: - ! - ! x̅ = n̅×e̅₃ / |n̅×e̅₃| = (n₂/c)e̅₁ + (-n₁/c)e̅₂ - ! y̅ = n̅×(n̅×e̅₃) / |n̅×e̅₃| = (n₁n₃/c)e̅₁ + (n₂n₃/c)e̅₂ - ce̅₃ - ! z̅ = n̅ = n₁e̅₁ + n₂e̅₂ + n₃e̅₃ - ! - ! where c = |n̅×e̅₃| = √(n₁² + n₂²) - ! - ! Or in the case where n̅ ∥ e̅₃ by: - ! - ! x̅ = e̅₁ - ! y̅ = (n₃e̅₃)×e̅₁ = n₃e̅₂ - ! z̅ = n₃e̅₃ - ! - ! So, the change of basis matrix is - ! - ! [n₂/c -n₁/c 0] - ! M = [n₁n₃/c n₂n₃/c -c], c > 0 - ! [n₁ n₂ n₃] - ! - ! or M = diag(1, n₃, n₃) if c = 0. - ! - ! By definition then we have: x̅_loc = M x̅. - ! n = y(4:6, 1) / norm2(y(4:6, 1)) - c = norm2(n(1:2)) - if (c > 0) then - M = reshape([ n(2)/c, n(1)*n(3)/c, n(1), & - -n(1)/c, n(2)*n(3)/c, n(2), & - zero, -c, n(3) ], [3, 3]) - else - M = reshape([ one, zero, zero, & - zero, n(3), zero, & - zero, zero, n(3) ], [3, 3]) - end if + M = tokamak2beam(n) ! start from the central ray or the last ring jk0 = 1 + merge(0, params%nray - params%nrayth, full_) diff --git a/src/utils.f90 b/src/utils.f90 index 7e38e08..29a89e3 100644 --- a/src/utils.f90 +++ b/src/utils.f90 @@ -190,4 +190,32 @@ contains digits = ceiling((bit_size(abs(n)) - leadz(abs(n))) * log10(2.0_wp_)) end function digits + + pure function diag(v) result(D) + ! Returns a matrix D with v as main diagonal + + ! function arguments + real(wp_), intent(in) :: v(:) + real(wp_) :: D(size(v), size(v)) + + integer :: i + + D = 0 + do concurrent (i = 1:size(v)) + D(i,i) = v(i) + end do + end function + + + pure function rotate(phi) result(R) + ! Returns a 2D rotation matrix + + ! function arguments + real(wp_), intent(in) :: phi + real(wp_) :: R(2,2) + + R = reshape([cos(phi), -sin(phi), & + sin(phi), cos(phi)], shape=[2,2], order=[2,1]) + end function rotate + end module utils