diff --git a/src/equilibrium.f90 b/src/equilibrium.f90 index b74bf69..97ffdcb 100644 --- a/src/equilibrium.f90 +++ b/src/equilibrium.f90 @@ -766,7 +766,7 @@ contains implicit none - real(wp_) :: dq + real(wp_) :: dq, gamma btaxis = model%B0 rmaxis = model%R0 @@ -782,9 +782,14 @@ contains zmxm = zbsup zmnm = zbinf - ! FIXME: it should be that Φ(r) = B₀πr² 1/√(1-ε²) - ! where ε=r/R₀ is the tokamak aspect ratio - phitedge = model%B0 * pi * model%a**2 + ! Toroidal flux at r=a: + ! + ! Φ(a) = B₀πa² 2γ/(γ + 1) + ! + ! where γ=1/√(1-ε²), + ! ε=a/R₀ is the tokamak aspect ratio + gamma = 1/sqrt(1 - (model%a/model%R0)**2) + phitedge = model%B0 * pi * model%a**2 * 2*gamma/(gamma + 1) ! In cocos=3 the safety factor is ! @@ -828,7 +833,8 @@ contains ! derivatives wrt (R, z) up to the second order. ! ! Note: all output arguments are optional. - use gray_params, only : iequil + use gray_params, only : iequil + use const_and_precisions, only : one, pi implicit none @@ -838,21 +844,75 @@ contains ddpsidrr, ddpsidzz, ddpsidrz ! local variables - real(wp_) :: rho_t, rho_p ! √Φ_n, √ψ_n + real(wp_) :: r_g, rho_t, rho_p ! geometric radius, √Φ_n, √ψ_n + real(wp_) :: gamma ! γ = 1/√(1 - r²/R₀²) real(wp_) :: dpsidphi ! (∂ψ_n/∂Φ_n) real(wp_) :: ddpsidphidr, ddpsidphidz ! ∇(∂ψ_n/∂Φ_n) - real(wp_) :: dphidr, dphidz ! ∇Φ_n + real(wp_) :: phi_n ! Φ_n + real(wp_) :: dphidr, dphidz ! ∇Φ_n + real(wp_) :: ddphidrdr, ddphidzdz ! ∇∇Φ_n + real(wp_) :: ddphidrdz ! real(wp_) :: q, dq ! q(ρ_p), Δq=(q₁-q₀)/(α/2 + 1) 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) - ! ρ_t is mapped to the normalised geometrical radius, - ! so ρ_t = √[(R-R₀)² + (z-z₀)²]/a and ρ_p(ρ_t) - rho_t = hypot(R - model%R0, z - model%z0)/model%a + ! 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 + + 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) + ! 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²) + ! where ∇(r²) = 2[(R-R₀), (z-z₀)] + dphidr = dphidr2 * 2*(R - model%R0) + dphidz = dphidr2 * 2*(z - model%z0) + + ! ∇∇Φ_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 = ρ_p(ρ_t)² if (present(psi_n)) psi_n = rho_p**2 @@ -871,10 +931,6 @@ contains dq = (model%q1 - model%q0) / (model%alpha/2 + 1) dpsidphi = (model%q0 + dq) / q - ! Since Φ_n = ρ_t²(R,z): ∇Φ_n = 2/a² [R-R₀, z-z₀] - dphidr = 2*(R - model%R0)/model%a**2 - dphidz = 2*(z - model%z0)/model%a**2 - ! Using the above, ∇ψ_n = ∂ψ_n/∂Φ_n ∇Φ_n if (present(dpsidr)) dpsidr = dpsidphi * dphidr if (present(dpsidz)) dpsidz = dpsidphi * dphidz @@ -890,18 +946,18 @@ contains ! ∇q = α/2 (q-q₀) ∇ψ_n/ψ_n ! = α/2 (q-q₀)/ψ_n (∂ψ_n/∂Φ_n) ∇Φ_n. ! - dqdr = model%alpha/2 * (q - model%q0)/rho_p**2 * dpsidphi * dphidr - dqdz = model%alpha/2 * (q - model%q0)/rho_p**2 * dpsidphi * dphidz + 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 - ! Finally, ∇∇Φ_n = ∇∇ ρ_t²(R,z) = 2/a² I, so: + ! Combining all of the above: ! - ! ∇∇ψ_n = ∇(∂ψ_n/∂Φ_n) ∇Φ_n + (∂ψ_n/∂Φ_n) 2/a² I + ! ∇∇ψ_n = ∇(∂ψ_n/∂Φ_n) ∇Φ_n + (∂ψ_n/∂Φ_n) ∇∇Φ_n ! - if (present(ddpsidrr)) ddpsidrr = ddpsidphidr * dphidr + dpsidphi * 2/model%a**2 - if (present(ddpsidzz)) ddpsidzz = ddpsidphidz * dphidz + dpsidphi * 2/model%a**2 - if (present(ddpsidrz)) ddpsidrz = ddpsidphidr * dphidz + 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 (R <= rmxm .and. R >= rmnm .and. &