From 7c5b4438470403dc30517301715c9976b85bafeb Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Tue, 30 Jan 2024 10:14:24 +0100 Subject: [PATCH] src/equilibrium.f90: use enums --- src/equilibrium.f90 | 621 +++++++++++++------------ src/gray_core.f90 | 74 +-- src/main.f90 | 10 +- tests/11-vacuum/inputs/equilibrium.txt | 4 - 4 files changed, 380 insertions(+), 329 deletions(-) delete mode 100644 tests/11-vacuum/inputs/equilibrium.txt diff --git a/src/equilibrium.f90 b/src/equilibrium.f90 index ee51ae8..a25dfb1 100644 --- a/src/equilibrium.f90 +++ b/src/equilibrium.f90 @@ -8,6 +8,8 @@ module equilibrium use const_and_precisions, only : wp_ use splines, only : spline_simple, spline_1d, spline_2d, linear_1d + use gray_params, only : EQ_VACUUM, EQ_ANALYTICAL, & + EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL implicit none @@ -360,36 +362,34 @@ contains ! 3. q = 1/2π ∂Φ/∂ψ ~ ∂Φ/∂r⋅∂r/∂ψ < 0. ! 4. In general, sgn(q) = -sgn(I_p)⋅sgn(B_φ). - if (iequil < 2) then - ! Analytical model + select case(iequil) + case (EQ_ANALYTICAL) ! Analytical model + ! Apply signs + if (params%sgni /= 0) then + model%q0 = sign(model%q0, -params%sgni*one) + model%q1 = sign(model%q1, -params%sgni*one) + end if + if (params%sgnb /= 0) then + model%B0 = sign(model%B0, +params%sgnb*one) + end if + ! Rescale + model%B0 = model%B0 * params%factb - ! Apply signs - if (params%sgni /= 0) then - model%q0 = sign(model%q0, -params%sgni*one) - model%q1 = sign(model%q1, -params%sgni*one) - end if - if (params%sgnb /= 0) then - model%B0 = sign(model%B0, +params%sgnb*one) - end if - ! Rescale - model%B0 = model%B0 * params%factb - else - ! Numeric data + case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL) ! Numeric data + ! Apply signs + if (params%sgni /= 0) & + data%psia = sign(data%psia, -params%sgni*one) + if (params%sgnb /= 0) & + data%fpol = sign(data%fpol, +params%sgnb*one) + ! Rescale + data%psia = data%psia * params%factb + data%fpol = data%fpol * params%factb - ! Apply signs - if (params%sgni /= 0) & - data%psia = sign(data%psia, -params%sgni*one) - if (params%sgnb /= 0) & - data%fpol = sign(data%fpol, +params%sgnb*one) - ! Rescale - data%psia = data%psia * params%factb - data%fpol = data%fpol * params%factb - - ! Compute the signs to be shown in the outputs header when cocos≠0,10. - ! Note: In these cases the values sgni,sgnb from gray.ini are unused. - params%sgni = int(sign(one, -data%psia)) - params%sgnb = int(sign(one, +data%fpol(size(data%fpol)))) - end if + ! Compute the signs to be shown in the outputs header when cocos≠0,10. + ! Note: In these cases the values sgni,sgnb from gray.ini are unused. + params%sgni = int(sign(one, -data%psia)) + params%sgnb = int(sign(one, +data%fpol(size(data%fpol)))) + end select end subroutine scale_equil @@ -398,7 +398,7 @@ contains ! in their respective global variables, see the top of this file. use const_and_precisions, only : zero, one use gray_params, only : equilibrium_parameters, equilibrium_data - use gray_params, only : iequil + use gray_params, only : iequil, X_AT_TOP, X_AT_BOTTOM, X_IS_MISSING use utils, only : vmaxmin, vmaxmini, inside use logger, only : log_info @@ -434,77 +434,80 @@ contains ! Spline interpolation of ψ(R, z) - if (iequil>2) then - ! data valid only inside boundary (data%psin=0 outside), e.g. source==ESCO - ! presence of boundary anticipated here to filter invalid data - nbnd = min(size(data%rbnd), size(data%zbnd)) + select case (iequil) - ! allocate knots and spline coefficients arrays - if (allocated(psi_spline%knots_x)) deallocate(psi_spline%knots_x) - if (allocated(psi_spline%knots_y)) deallocate(psi_spline%knots_y) - if (allocated(psi_spline%coeffs)) deallocate(psi_spline%coeffs) - allocate(psi_spline%knots_x(nrest), psi_spline%knots_y(nzest)) - allocate(psi_spline%coeffs(nrz)) + case (EQ_EQDSK_PARTIAL) + ! Data valid only inside boundary (data%psin=0 outside), + ! presence of boundary anticipated here to filter invalid data + nbnd = min(size(data%rbnd), size(data%zbnd)) - ! determine number of valid grid points - nrz=0 - do j=1,nz - do i=1,nr - if (nbnd.gt.0) then - if(.not.inside(data%rbnd,data%zbnd,data%rv(i),data%zv(j))) cycle - else - if(data%psin(i,j).le.0.0d0) cycle - end if - nrz=nrz+1 + ! allocate knots and spline coefficients arrays + if (allocated(psi_spline%knots_x)) deallocate(psi_spline%knots_x) + if (allocated(psi_spline%knots_y)) deallocate(psi_spline%knots_y) + if (allocated(psi_spline%coeffs)) deallocate(psi_spline%coeffs) + allocate(psi_spline%knots_x(nrest), psi_spline%knots_y(nzest)) + allocate(psi_spline%coeffs(nrz)) + + ! determine number of valid grid points + nrz=0 + do j=1,nz + do i=1,nr + if (nbnd.gt.0) then + if(.not.inside(data%rbnd,data%zbnd,data%rv(i),data%zv(j))) cycle + else + if(data%psin(i,j).le.0.0d0) cycle + end if + nrz=nrz+1 + end do end do - end do - ! store valid data - allocate(rv1d(nrz),zv1d(nrz),fvpsi(nrz),wf(nrz)) - ij=0 - do j=1,nz - do i=1,nr - if (nbnd.gt.0) then - if(.not.inside(data%rbnd,data%zbnd,data%rv(i),data%zv(j))) cycle - else - if(data%psin(i,j).le.0.0d0) cycle - end if - ij=ij+1 - rv1d(ij)=data%rv(i) - zv1d(ij)=data%zv(j) - fvpsi(ij)=data%psin(i,j) - wf(ij)=1.0d0 + ! store valid data + allocate(rv1d(nrz),zv1d(nrz),fvpsi(nrz),wf(nrz)) + ij=0 + do j=1,nz + do i=1,nr + if (nbnd.gt.0) then + if(.not.inside(data%rbnd,data%zbnd,data%rv(i),data%zv(j))) cycle + else + if(data%psin(i,j).le.0.0d0) cycle + end if + ij=ij+1 + rv1d(ij)=data%rv(i) + zv1d(ij)=data%zv(j) + fvpsi(ij)=data%psin(i,j) + wf(ij)=1.0d0 + end do end do - end do - ! Fit as a scattered set of points - ! use reduced number of knots to limit memory comsumption ? - psi_spline%nknots_x=nr/4+4 - psi_spline%nknots_y=nz/4+4 - tension = params%ssplps - call scatterspl(rv1d, zv1d, fvpsi, wf, nrz, kspl, tension, & - rmnm, rmxm, zmnm, zmxm, & - psi_spline%knots_x, psi_spline%nknots_x, & - psi_spline%knots_y, psi_spline%nknots_y, & - psi_spline%coeffs, err) - - ! if failed, re-fit with an interpolating spline (zero tension) - if(err == -1) then - err = 0 - tension = 0 + ! Fit as a scattered set of points + ! use reduced number of knots to limit memory comsumption ? psi_spline%nknots_x=nr/4+4 psi_spline%nknots_y=nz/4+4 + tension = params%ssplps call scatterspl(rv1d, zv1d, fvpsi, wf, nrz, kspl, tension, & rmnm, rmxm, zmnm, zmxm, & psi_spline%knots_x, psi_spline%nknots_x, & psi_spline%knots_y, psi_spline%nknots_y, & psi_spline%coeffs, err) - end if - deallocate(rv1d, zv1d, wf, fvpsi) - ! reset nrz to the total number of grid points for next allocations - nrz = nr*nz - else - ! iequil==2: data are valid on the full R,z grid + + ! if failed, re-fit with an interpolating spline (zero tension) + if(err == -1) then + err = 0 + tension = 0 + psi_spline%nknots_x=nr/4+4 + psi_spline%nknots_y=nz/4+4 + call scatterspl(rv1d, zv1d, fvpsi, wf, nrz, kspl, tension, & + rmnm, rmxm, zmnm, zmxm, & + psi_spline%knots_x, psi_spline%nknots_x, & + psi_spline%knots_y, psi_spline%nknots_y, & + psi_spline%coeffs, err) + end if + deallocate(rv1d, zv1d, wf, fvpsi) + ! reset nrz to the total number of grid points for next allocations + nrz = nr*nz + + case (EQ_EQDSK_FULL) + ! Data are valid on the full R,z grid ! reshape 2D ψ array to 1D (transposed) allocate(fvpsi(nrz)) @@ -523,7 +526,7 @@ contains err = 0 end if deallocate(fvpsi) - end if + end select if (err /= 0) then err = 2 @@ -592,11 +595,11 @@ contains 'r', rmaxis, 'z', zmaxis, 'ψ', psinoptmp call log_info(msg, mod='equilibrium', proc='set_equil_spline') - ! search for X-point if params%ixp /= 0 - + ! Search for X-point ixploc = params%ixp - if(ixploc/=0) then - if(ixploc<0) then + + select case (ixploc) + case (X_AT_BOTTOM) call points_ox(rbinf,zbinf,r1,z1,psinxptmp,info) if(psinxptmp/=-1.0_wp_) then write (msg, '("X-point found:", 3(x,a,"=",g0.3))') & @@ -611,7 +614,8 @@ contains else ixploc=0 end if - else + + case (X_AT_TOP) call points_ox(rbsup,zbsup,r1,z1,psinxptmp,info) if(psinxptmp.ne.-1.0_wp_) then write (msg, '("X-point found:", 3(x,a,"=",g0.3))') & @@ -626,26 +630,22 @@ contains else ixploc=0 end if - end if - end if - if (ixploc==0) then - psinop=psinoptmp - psiant=one-psinop - - ! Find upper horizontal tangent point - call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbsup),r1,z1,one,info) - zbsup=z1 - rbsup=r1 - - ! Find lower horizontal tangent point - call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbinf),r1,z1,one,info) - zbinf=z1 - rbinf=r1 - write (msg, '("X-point not found in", 2(x,a,"∈[",g0.3,",",g0.3,"]"))') & - 'r', rbinf, rbsup, 'z', zbinf, zbsup - call log_info(msg, mod='equilibrium', proc='set_equil_spline') - end if + case (X_IS_MISSING) + psinop=psinoptmp + psiant=one-psinop + ! Find upper horizontal tangent point + call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbsup),r1,z1,one,info) + zbsup=z1 + rbsup=r1 + ! Find lower horizontal tangent point + call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbinf),r1,z1,one,info) + zbinf=z1 + rbinf=r1 + write (msg, '("X-point not found in", 2(x,a,"∈[",g0.3,",",g0.3,"]"))') & + 'r', rbinf, rbsup, 'z', zbinf, zbsup + call log_info(msg, mod='equilibrium', proc='set_equil_spline') + end select ! Adjust all the B-spline coefficients ! Note: since ψ_n(R,z) = Σ_ij c_ij B_i(R)B_j(z), to correct ψ_n @@ -908,128 +908,131 @@ contains real(wp_) :: dqdr, dqdz ! ∇q real(wp_) :: dphidr2, ddphidr2dr2 ! dΦ_n/d(r²), d²Φ_n/d(r²)² - if (iequil < 2) then - ! Analytical model - ! - ! The normalised poloidal flux ψ_n(R, z) is computed as follows: - ! 1. ψ_n = ρ_p² - ! 2. ρ_p = ρ_p(ρ_t), using `frhopol`, which in turns uses q(ψ) - ! 3. ρ_t = √Φ_n - ! 4. Φ_n = Φ(r)/Φ(a), where Φ(r) is the flux of B_φ=B₀R₀/R - ! through a circular surface - ! 5. r = √[(R-R₀)²+(z-z₀)²] is the geometric minor radius - r_g = hypot(R - model%R0, z - model%z0) + ! Values for vacuum/outside the domain + if (present(psi_n)) psi_n = -1 + if (present(dpsidr)) dpsidr = 0 + if (present(dpsidz)) dpsidz = 0 + if (present(ddpsidrr)) ddpsidrr = 0 + if (present(ddpsidzz)) ddpsidzz = 0 + if (present(ddpsidrz)) ddpsidrz = 0 - ! The exact flux of the toroidal field B_φ = B₀R₀/R is: - ! - ! Φ(r) = B₀πr² 2γ/(γ + 1) where γ=1/√(1 - r²/R₀²). - ! - ! Notes: - ! 1. the function Φ(r) is defined for r≤R₀ only. - ! 2. as r → 0, γ → 1, so Φ ~ B₀πr². - ! 3. as r → 1⁻, Φ → 2B₀πr² but dΦ/dr → -∞. - ! 4. |B_R|, |B_z| → +-∞. - ! - if (r_g > model%R0) then - if (present(psi_n)) psi_n = -1 - if (present(dpsidr)) dpsidr = 0 - if (present(dpsidz)) dpsidz = 0 - if (present(ddpsidrr)) ddpsidrr = 0 - if (present(ddpsidzz)) ddpsidzz = 0 - if (present(ddpsidrz)) ddpsidrz = 0 - return - end if + select case (iequil) + case (EQ_ANALYTICAL) + ! Analytical model + ! + ! The normalised poloidal flux ψ_n(R, z) is computed as follows: + ! 1. ψ_n = ρ_p² + ! 2. ρ_p = ρ_p(ρ_t), using `frhopol`, which in turns uses q(ψ) + ! 3. ρ_t = √Φ_n + ! 4. Φ_n = Φ(r)/Φ(a), where Φ(r) is the flux of B_φ=B₀R₀/R + ! through a circular surface + ! 5. r = √[(R-R₀)²+(z-z₀)²] is the geometric minor radius + r_g = hypot(R - model%R0, z - model%z0) - gamma = 1 / sqrt(1 - (r_g/model%R0)**2) - phi_n = model%B0 * pi*r_g**2 * 2*gamma/(gamma + 1) / phitedge - rho_t = sqrt(phi_n) - rho_p = frhopol(rho_t) + ! The exact flux of the toroidal field B_φ = B₀R₀/R is: + ! + ! Φ(r) = B₀πr² 2γ/(γ + 1) where γ=1/√(1 - r²/R₀²). + ! + ! Notes: + ! 1. the function Φ(r) is defined for r≤R₀ only. + ! 2. as r → 0, γ → 1, so Φ ~ B₀πr². + ! 3. as r → 1⁻, Φ → 2B₀πr² but dΦ/dr → -∞. + ! 4. |B_R|, |B_z| → +-∞. + ! + if (r_g > model%R0) then + if (present(psi_n)) psi_n = -1 + if (present(dpsidr)) dpsidr = 0 + if (present(dpsidz)) dpsidz = 0 + if (present(ddpsidrr)) ddpsidrr = 0 + if (present(ddpsidzz)) ddpsidzz = 0 + if (present(ddpsidrz)) ddpsidrz = 0 + return + end if - ! For ∇Φ_n and ∇∇Φ_n we also need: - ! - ! ∂Φ∂(r²) = B₀π γ(r) - ! ∂²Φ∂(r²)² = B₀π γ³(r) / (2 R₀²) - ! - dphidr2 = model%B0 * pi * gamma / phitedge - ddphidr2dr2 = model%B0 * pi * gamma**3/(2 * model%R0**2) / phitedge + gamma = 1 / sqrt(1 - (r_g/model%R0)**2) + phi_n = model%B0 * pi*r_g**2 * 2*gamma/(gamma + 1) / phitedge + rho_t = sqrt(phi_n) + rho_p = frhopol(rho_t) - ! ∇Φ_n = ∂Φ_n/∂(r²) ∇(r²) - ! where ∇(r²) = 2[(R-R₀), (z-z₀)] - dphidr = dphidr2 * 2*(R - model%R0) - dphidz = dphidr2 * 2*(z - model%z0) + ! For ∇Φ_n and ∇∇Φ_n we also need: + ! + ! ∂Φ∂(r²) = B₀π γ(r) + ! ∂²Φ∂(r²)² = B₀π γ³(r) / (2 R₀²) + ! + dphidr2 = model%B0 * pi * gamma / phitedge + ddphidr2dr2 = model%B0 * pi * gamma**3/(2 * model%R0**2) / phitedge - ! ∇∇Φ_n = ∇[∂Φ_n/∂(r²)] ∇(r²) + ∂Φ_n/∂(r²) ∇∇(r²) - ! = ∂²Φ_n/∂(r²)² ∇(r²)∇(r²) + ∂Φ_n/∂(r²) ∇∇(r²) - ! where ∇∇(r²) = 2I - ddphidrdr = ddphidr2dr2 * 4*(R - model%R0)*(R - model%R0) + dphidr2*2 - ddphidzdz = ddphidr2dr2 * 4*(z - model%z0)*(z - model%z0) + dphidr2*2 - ddphidrdz = ddphidr2dr2 * 4*(R - model%R0)*(z - model%z0) + ! ∇Φ_n = ∂Φ_n/∂(r²) ∇(r²) + ! where ∇(r²) = 2[(R-R₀), (z-z₀)] + dphidr = dphidr2 * 2*(R - model%R0) + dphidz = dphidr2 * 2*(z - model%z0) - ! ψ_n = ρ_p(ρ_t)² - if (present(psi_n)) psi_n = rho_p**2 + ! ∇∇Φ_n = ∇[∂Φ_n/∂(r²)] ∇(r²) + ∂Φ_n/∂(r²) ∇∇(r²) + ! = ∂²Φ_n/∂(r²)² ∇(r²)∇(r²) + ∂Φ_n/∂(r²) ∇∇(r²) + ! where ∇∇(r²) = 2I + ddphidrdr = ddphidr2dr2 * 4*(R - model%R0)*(R - model%R0) + dphidr2*2 + ddphidzdz = ddphidr2dr2 * 4*(z - model%z0)*(z - model%z0) + dphidr2*2 + ddphidrdz = ddphidr2dr2 * 4*(R - model%R0)*(z - model%z0) - ! Using the definitions in `frhotor`: - ! - ! ∇ψ_n = ∂ψ_n/∂Φ_n ∇Φ_n - ! - ! ∂ψ_n/∂Φ_n = Φ_a/ψ_a ∂ψ/∂Φ - ! = Φ_a/ψ_a 1/2πq - ! - ! Using ψ_a = 1/2π Φ_a / (q₀ + Δq), then: - ! - ! ∂ψ_n/∂Φ_n = (q₀ + Δq)/q - ! - q = model%q0 + (model%q1 - model%q0) * rho_p**model%alpha - dq = (model%q1 - model%q0) / (model%alpha/2 + 1) - dpsidphi = (model%q0 + dq) / q + ! ψ_n = ρ_p(ρ_t)² + if (present(psi_n)) psi_n = rho_p**2 - ! Using the above, ∇ψ_n = ∂ψ_n/∂Φ_n ∇Φ_n - if (present(dpsidr)) dpsidr = dpsidphi * dphidr - if (present(dpsidz)) dpsidz = dpsidphi * dphidz + ! Using the definitions in `frhotor`: + ! + ! ∇ψ_n = ∂ψ_n/∂Φ_n ∇Φ_n + ! + ! ∂ψ_n/∂Φ_n = Φ_a/ψ_a ∂ψ/∂Φ + ! = Φ_a/ψ_a 1/2πq + ! + ! Using ψ_a = 1/2π Φ_a / (q₀ + Δq), then: + ! + ! ∂ψ_n/∂Φ_n = (q₀ + Δq)/q + ! + q = model%q0 + (model%q1 - model%q0) * rho_p**model%alpha + dq = (model%q1 - model%q0) / (model%alpha/2 + 1) + dpsidphi = (model%q0 + dq) / q - ! For the second derivatives: - ! - ! ∇∇ψ_n = ∇(∂ψ_n/∂Φ_n) ∇Φ_n + (∂ψ_n/∂Φ_n) ∇∇Φ_n - ! - ! ∇(∂ψ_n/∂Φ_n) = - (∂ψ_n/∂Φ_n) ∇q/q - ! - ! From q(ψ) = q₀ + (q₁-q₀) ψ_n^α/2, we have: - ! - ! ∇q = α/2 (q-q₀) ∇ψ_n/ψ_n - ! = α/2 (q-q₀)/ψ_n (∂ψ_n/∂Φ_n) ∇Φ_n. - ! - dqdr = model%alpha/2 * (model%q1 - model%q0)*rho_p**(model%alpha-2) * dpsidphi * dphidr - dqdz = model%alpha/2 * (model%q1 - model%q0)*rho_p**(model%alpha-2) * dpsidphi * dphidz - ddpsidphidr = - dpsidphi * dqdr/q - ddpsidphidz = - dpsidphi * dqdz/q + ! Using the above, ∇ψ_n = ∂ψ_n/∂Φ_n ∇Φ_n + if (present(dpsidr)) dpsidr = dpsidphi * dphidr + if (present(dpsidz)) dpsidz = dpsidphi * dphidz - ! Combining all of the above: - ! - ! ∇∇ψ_n = ∇(∂ψ_n/∂Φ_n) ∇Φ_n + (∂ψ_n/∂Φ_n) ∇∇Φ_n - ! - if (present(ddpsidrr)) ddpsidrr = ddpsidphidr * dphidr + dpsidphi * ddphidrdr - if (present(ddpsidzz)) ddpsidzz = ddpsidphidz * dphidz + dpsidphi * ddphidzdz - if (present(ddpsidrz)) ddpsidrz = ddpsidphidr * dphidz + dpsidphi * ddphidrdz - else - ! Numerical data - if (inside(psi_domain%R, psi_domain%z, R, z)) then - ! Within the interpolation range - if (present(psi_n)) psi_n = psi_spline%eval(R, z) - if (present(dpsidr)) dpsidr = psi_spline%deriv(R, z, 1, 0) - if (present(dpsidz)) dpsidz = psi_spline%deriv(R, z, 0, 1) - if (present(ddpsidrr)) ddpsidrr = psi_spline%deriv(R, z, 2, 0) - if (present(ddpsidzz)) ddpsidzz = psi_spline%deriv(R, z, 0, 2) - if (present(ddpsidrz)) ddpsidrz = psi_spline%deriv(R, z, 1, 1) - else - ! Outside - if (present(psi_n)) psi_n = -1 - if (present(dpsidr)) dpsidr = 0 - if (present(dpsidz)) dpsidz = 0 - if (present(ddpsidrr)) ddpsidrr = 0 - if (present(ddpsidzz)) ddpsidzz = 0 - if (present(ddpsidrz)) ddpsidrz = 0 - end if - end if + ! For the second derivatives: + ! + ! ∇∇ψ_n = ∇(∂ψ_n/∂Φ_n) ∇Φ_n + (∂ψ_n/∂Φ_n) ∇∇Φ_n + ! + ! ∇(∂ψ_n/∂Φ_n) = - (∂ψ_n/∂Φ_n) ∇q/q + ! + ! From q(ψ) = q₀ + (q₁-q₀) ψ_n^α/2, we have: + ! + ! ∇q = α/2 (q-q₀) ∇ψ_n/ψ_n + ! = α/2 (q-q₀)/ψ_n (∂ψ_n/∂Φ_n) ∇Φ_n. + ! + dqdr = model%alpha/2 * (model%q1 - model%q0)*rho_p**(model%alpha-2) * dpsidphi * dphidr + dqdz = model%alpha/2 * (model%q1 - model%q0)*rho_p**(model%alpha-2) * dpsidphi * dphidz + ddpsidphidr = - dpsidphi * dqdr/q + ddpsidphidz = - dpsidphi * dqdz/q + + ! Combining all of the above: + ! + ! ∇∇ψ_n = ∇(∂ψ_n/∂Φ_n) ∇Φ_n + (∂ψ_n/∂Φ_n) ∇∇Φ_n + ! + if (present(ddpsidrr)) ddpsidrr = ddpsidphidr * dphidr + dpsidphi * ddphidrdr + if (present(ddpsidzz)) ddpsidzz = ddpsidphidz * dphidz + dpsidphi * ddphidzdz + if (present(ddpsidrz)) ddpsidrz = ddpsidphidr * dphidz + dpsidphi * ddphidrdz + + case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL) + ! Numerical data + if (inside(psi_domain%R, psi_domain%z, R, z)) then + ! Within the interpolation range + if (present(psi_n)) psi_n = psi_spline%eval(R, z) + if (present(dpsidr)) dpsidr = psi_spline%deriv(R, z, 1, 0) + if (present(dpsidz)) dpsidz = psi_spline%deriv(R, z, 0, 1) + if (present(ddpsidrr)) ddpsidrr = psi_spline%deriv(R, z, 2, 0) + if (present(ddpsidzz)) ddpsidzz = psi_spline%deriv(R, z, 0, 2) + if (present(ddpsidrz)) ddpsidrz = psi_spline%deriv(R, z, 1, 1) + end if + + end select end subroutine pol_flux @@ -1043,21 +1046,28 @@ contains real(wp_), intent(out) :: fpol ! poloidal current real(wp_), intent(out), optional :: dfpol ! derivative - if (iequil < 2) then - ! Analytical model - ! F(ψ) = B₀⋅R₀, a constant - fpol = model%B0 * model%R0 - if (present(dfpol)) dfpol = 0 - else - ! Numerical data - if(psi_n <= 1 .and. psi_n >= 0) then - fpol = fpol_spline%eval(psi_n) - if (present(dfpol)) dfpol = fpol_spline%deriv(psi_n) - else - fpol = fpolas + select case (iequil) + case (EQ_VACUUM) + ! Vacuum, no plasma + fpol = 0 if (present(dfpol)) dfpol = 0 - end if - end if + + case (EQ_ANALYTICAL) + ! Analytical model + ! F(ψ) = B₀⋅R₀, a constant + fpol = model%B0 * model%R0 + if (present(dfpol)) dfpol = 0 + + case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL) + ! Numerical data + if(psi_n <= 1 .and. psi_n >= 0) then + fpol = fpol_spline%eval(psi_n) + if (present(dfpol)) dfpol = fpol_spline%deriv(psi_n) + else + fpol = fpolas + if (present(dfpol)) dfpol = 0 + end if + end select end subroutine pol_curr @@ -1069,30 +1079,34 @@ contains real(wp_), intent(in) :: rho_p real(wp_) :: frhotor - if (iequil < 2) then - ! Analytical model - block - ! The change of variable is obtained by integrating - ! - ! q(ψ) = 1/2π ∂Φ/∂ψ - ! - ! and defining ψ = ψ_a ρ_p², Φ = Φ_a ρ_t². - ! The result is: - ! - ! - ψ_a = 1/2π Φ_a / [q₀ + Δq] - ! - ! - ρ_t = ρ_p √[(q₀ + Δq ρ_p^α)/(q₀ + Δq)] - ! - ! where Δq = (q₁ - q₀)/(α/2 + 1) - real(wp_) :: dq - dq = (model%q1 - model%q0) / (model%alpha/2 + 1) - frhotor = rho_p * sqrt((model%q0 + dq*rho_p**model%alpha) & - / (model%q0 + dq)) - end block - else - ! Numerical data - frhotor = rhot_spline%eval(rho_p) - end if + select case (iequil) + + case (EQ_ANALYTICAL) + ! Analytical model + block + ! The change of variable is obtained by integrating + ! + ! q(ψ) = 1/2π ∂Φ/∂ψ + ! + ! and defining ψ = ψ_a ρ_p², Φ = Φ_a ρ_t². + ! The result is: + ! + ! - ψ_a = 1/2π Φ_a / [q₀ + Δq] + ! + ! - ρ_t = ρ_p √[(q₀ + Δq ρ_p^α)/(q₀ + Δq)] + ! + ! where Δq = (q₁ - q₀)/(α/2 + 1) + real(wp_) :: dq + dq = (model%q1 - model%q0) / (model%alpha/2 + 1) + frhotor = rho_p * sqrt((model%q0 + dq*rho_p**model%alpha) & + / (model%q0 + dq)) + end block + + case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL) + ! Numerical data + frhotor = rhot_spline%eval(rho_p) + + end select end function frhotor @@ -1105,26 +1119,32 @@ contains real(wp_), intent(in) :: rho_t real(wp_) :: frhopol - if (iequil < 2) then - ! Analytical model - block - ! In general there is no closed form for ρ_p(ρ_t) in the - ! analytical model, we thus solve numerically the equation - ! ρ_t(ρ_p) = ρ_t₀ for ρ_p. - use minpack, only : hybrj1 + select case (iequil) + case (EQ_VACUUM) + ! Vacuum, no plasma + frhopol = 0 - real(wp_) :: rho_p(1), fvec(1), fjac(1,1), wa(7) - integer :: info + case (EQ_ANALYTICAL) + ! Analytical model + block + ! In general there is no closed form for ρ_p(ρ_t) in the + ! analytical model, we thus solve numerically the equation + ! ρ_t(ρ_p) = ρ_t₀ for ρ_p. + use minpack, only : hybrj1 - rho_p = [rho_t] ! first guess, ρ_p ≈ ρ_t - call hybrj1(equation, n=1, x=rho_p, fvec=fvec, fjac=fjac, & - ldfjac=1, tol=comp_eps, info=info, wa=wa, lwa=7) - frhopol = rho_p(1) - end block - else - ! Numerical data - frhopol = rhop_spline%eval(rho_t) - end if + real(wp_) :: rho_p(1), fvec(1), fjac(1,1), wa(7) + integer :: info + + rho_p = [rho_t] ! first guess, ρ_p ≈ ρ_t + call hybrj1(equation, n=1, x=rho_p, fvec=fvec, fjac=fjac, & + ldfjac=1, tol=comp_eps, info=info, wa=wa, lwa=7) + frhopol = rho_p(1) + end block + + case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL) + ! Numerical data + frhopol = rhop_spline%eval(rho_t) + end select contains @@ -1162,19 +1182,25 @@ contains real(wp_), intent(in) :: psin real(wp_) :: fq - if (iequil < 2) then - ! Analytical model - ! The safety factor is a power law in ρ_p: - ! q(ρ_p) = q₀ + (q₁-q₀) ρ_p^α - block - real(wp_) :: rho_p - rho_p = sqrt(psin) - fq = abs(model%q0 + (model%q1 - model%q0) * rho_p**model%alpha) - end block - else - ! Numerical data - fq = q_spline%eval(psin) - end if + select case(iequil) + case (EQ_VACUUM) + ! Vacuum, q is undefined + fq = 0 + + case (EQ_ANALYTICAL) + ! Analytical model + ! The safety factor is a power law in ρ_p: + ! q(ρ_p) = q₀ + (q₁-q₀) ρ_p^α + block + real(wp_) :: rho_p + rho_p = sqrt(psin) + fq = abs(model%q0 + (model%q1 - model%q0) * rho_p**model%alpha) + end block + + case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL) + ! Numerical data + fq = q_spline%eval(psin) + end select end function fq @@ -1183,6 +1209,7 @@ contains ! (R, z) in cylindrical coordinates ! ! Note: all output arguments are optional. + use gray_params, only : iequil ! subroutine arguments real(wp_), intent(in) :: R, z @@ -1191,6 +1218,14 @@ contains ! local variables real(wp_) :: psi_n, fpol, dpsidr, dpsidz + if (iequil == EQ_VACUUM) then + ! Vacuum, no plasma nor field + if (present(B_R)) B_R = 0 + if (present(B_z)) B_z = 0 + if (present(B_phi)) B_phi = 0 + return + end if + call pol_flux(R, z, psi_n, dpsidr, dpsidz) call pol_curr(psi_n, fpol) diff --git a/src/gray_core.f90 b/src/gray_core.f90 index 5e100b0..bb3c2e8 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -12,7 +12,7 @@ contains use coreprofiles, only : temp, fzeff use dispersion, only : expinit use gray_params, only : gray_parameters, gray_data, gray_results, & - print_parameters + print_parameters, EQ_VACUUM use beams, only : xgygcoeff, launchangles2n use beamdata, only : pweight, rayi2jk use gray_errors, only : is_critical, print_err_raytracing, print_err_ecrh_cd @@ -100,12 +100,14 @@ contains ! Initialise the dispersion module if(params%ecrh_cd%iwarm > 1) call expinit - ! Initialise the magsurf_data module - call flux_average ! requires frhotor for dadrhot,dvdrhot + if (params%equilibrium%iequil /= EQ_VACUUM) then + ! Initialise the magsurf_data module + call flux_average ! requires frhotor for dadrhot,dvdrhot - ! Initialise the output profiles - call pec_init(params%output%ipec, rhout) - nnd = size(rhop_tab) ! number of radial profile points + ! Initialise the output profiles + call pec_init(params%output%ipec, rhout) + nnd = size(rhop_tab) ! number of radial profile points + end if call alloc_multipass(nnd, iwait, iroff, iop, iow, yynext, yypnext, yw0, ypw0, stnext, & stv, p0ray, taus, tau1, etau1, cpls, cpl1, lgcpl1, jphi_beam, & @@ -138,11 +140,13 @@ contains ! print Btot=Bres ! print ne, Te, q, Jphi versus psi, rhop, rhot - call print_bres(bres) - call print_prof(params%profiles) - call print_maps(bres, xgcn, & - norm2(params%antenna%pos(1:2)) * 0.01_wp_, & - sin(params%antenna%beta*degree)) + if (params%equilibrium%iequil /= EQ_VACUUM) then + call print_bres(bres) + call print_prof(params%profiles) + call print_maps(bres, xgcn, & + norm2(params%antenna%pos(1:2)) * 0.01_wp_, & + sin(params%antenna%beta*degree)) + end if ! ========= pre-proc prints END ========= ! =========== main loop BEGIN =========== @@ -518,9 +522,11 @@ contains icd_beam = sum(ccci(:,i)) call vmaxmin(tau0,params%raytracing%nray,taumn,taumx) ! taumn,taumx for print - ! compute power and current density profiles for all rays - call spec(psjki,ppabs,ccci,iiv,pabs_beam,icd_beam,dpdv_beam,jphi_beam,jcd_beam, & - pins_beam,currins_beam) + if (params%equilibrium%iequil /= EQ_VACUUM) then + ! compute power and current density profiles for all rays + call spec(psjki,ppabs,ccci,iiv,pabs_beam,icd_beam,dpdv_beam, & + jphi_beam,jcd_beam,pins_beam,currins_beam) + end if pabs_pass(iox) = pabs_pass(iox) + pabs_beam ! 0D results for current pass, sum on O/X mode beams icd_pass(iox) = icd_pass(iox) + icd_beam @@ -563,17 +569,20 @@ contains end if end if - call print_pec(rhop_tab,rhot_tab,jphi_beam,jcd_beam,dpdv_beam,currins_beam, & - pins_beam,ip) ! *print power and current density profiles for current beam + if (params%equilibrium%iequil /= EQ_VACUUM) then + call print_pec(rhop_tab,rhot_tab,jphi_beam,jcd_beam, & + dpdv_beam,currins_beam,pins_beam,ip) ! *print power and current density profiles for current beam - call postproc_profiles(pabs_beam,icd_beam,rhot_tab,dpdv_beam,jphi_beam, & - rhotpav,drhotpav,rhotjava,drhotjava,dpdvp,jphip,rhotp,drhotp,rhotj, & - drhotj,dpdvmx,jphimx,ratjamx,ratjbmx) ! *compute profiles width for current beam + call postproc_profiles(pabs_beam,icd_beam,rhot_tab,dpdv_beam, & + jphi_beam, rhotpav,drhotpav,rhotjava,drhotjava,dpdvp,jphip, & + rhotp,drhotp,rhotj,drhotj,dpdvmx,jphimx,ratjamx,ratjbmx) ! *compute profiles width for current beam + + call print_finals(pabs_beam,icd_beam,dpdvp,jphip,rhotpav, & + rhotjava,drhotpav,drhotjava,dpdvmx,jphimx,rhotp,rhotj, & + drhotp,drhotj,ratjamx,ratjbmx,stv(1),psipv(index_rt), & + chipv(index_rt),index_rt,sum(p0ray),cpl_beam1,cpl_beam2) ! *print 0D results for current beam + end if - call print_finals(pabs_beam,icd_beam,dpdvp,jphip,rhotpav,rhotjava, & - drhotpav,drhotjava,dpdvmx,jphimx,rhotp,rhotj,drhotp,drhotj,ratjamx, & - ratjbmx,stv(1),psipv(index_rt),chipv(index_rt),index_rt,sum(p0ray), & - cpl_beam1,cpl_beam2) ! *print 0D results for current beam ! ============ post-proc END ============ end do beam_loop @@ -1872,13 +1881,20 @@ contains bv(1)=br*csphi-bphi*snphi bv(2)=br*snphi+bphi*csphi bv(3)=bz - call pol_limit(anv,bv,bres,sox,ext0(jk),eyt0(jk)) - if (jk == 1) then - call stokes_ce(ext0(jk),eyt0(jk),qq,uu,vv) - call polellipse(qq,uu,vv,psipol0,chipol0) - psipol0=psipol0/degree ! convert from rad to degree - chipol0=chipol0/degree + if (norm2(bv) > 0) then + call pol_limit(anv,bv,bres,sox,ext0(jk),eyt0(jk)) + + if (jk == 1) then + call stokes_ce(ext0(jk),eyt0(jk),qq,uu,vv) + call polellipse(qq,uu,vv,psipol0,chipol0) + psipol0=psipol0/degree ! convert from rad to degree + chipol0=chipol0/degree + end if + else + ! X/O mode are undefined if B=0 + psipol0 = 0 + chipol0 = 0 end if else call stokes_ell(chipol0*degree,psipol0*degree,qq,uu,vv) diff --git a/src/main.f90 b/src/main.f90 index 8bd982d..3160d85 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -242,7 +242,7 @@ contains ! Set global variables select case (params%equilibrium%iequil) - case (EQ_VACUUM, EQ_ANALYTICAL) + case (EQ_ANALYTICAL) call set_equil_an case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL) @@ -387,7 +387,7 @@ contains EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL use utils, only : range2rect use limiter, only : limiter_set_globals=>set_globals - use const_and_precisions, only : cm + use const_and_precisions, only : cm, comp_huge ! subroutine arguments type(gray_parameters), intent(inout) :: params @@ -415,7 +415,11 @@ contains ! Max radius, either due to the plasma extent or equilibrium grid select case (params%equilibrium%iequil) - case (EQ_VACUUM, EQ_ANALYTICAL) + case (EQ_VACUUM) + ! Use a very large R, ~ unbounded + R_max = comp_huge + + case (EQ_ANALYTICAL) ! use R₀+a block use equilibrium, only : model diff --git a/tests/11-vacuum/inputs/equilibrium.txt b/tests/11-vacuum/inputs/equilibrium.txt deleted file mode 100644 index c02474f..0000000 --- a/tests/11-vacuum/inputs/equilibrium.txt +++ /dev/null @@ -1,4 +0,0 @@ -0.0 0.0 0.0 : rr0m,zr0m,rpam ! rhot[m] = min(sqrt((r-rr0m)**2+(z-zr0m)**2), rpam); tor flux phi = pi*b0*rhot**2 -0.0 : b0 ! Bphi[T] @ rr0m[m] -1 3 2 : q0, qa, alq ! q = q0 + (qa-q0)*sqrt(psin)**alq -0 : nlim ! number of points in first wall (rlim,zlim) polygon