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

src/equilibrium.f90: re-introduce ρ_t(R,z) mapping

Browse Source 
Directly mapping r(R,z) to ρ_p greatly simplifies the implementation of
pol_flux, but also produces a weird current profile and resonance curve.

This keeps the previous changes to the model but reverts the mapping to
the previous one: r → ρ_t.
This commit is contained in:
Michele Guerini Rocco 2023-12-06 17:26:21 +01:00
parent 24edfdc43a
commit 9a301f5799
Signed by: rnhmjoj
GPG Key ID: BFBAF4C975F76450
1 changed files with 77 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

90
src/equilibrium.f90
View File

@ -782,7 +782,27 @@ 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
! In cocos=3 the safety factor is
!
! q(ψ) = 1/2π ∂Φ/∂ψ.
!
! Given the power law of the model
!
! q(ψ) = q₀ + (q₁-q₀) (ψ/ψa)^(α/2),
!
! we can find ψ_a = ψ(r=a) by integrating:
!
! ∫ q(ψ)dψ = 1/2π ∫ dΦ
! ∫₀^ψ_a q(ψ)dψ = 1/2π Φ_a
! ψa [q₀ + (q₁-q₀)/(α/2+1)] = Φa/2π
!
! ⇒ ψ_a = Φ_a 1/2π 1/(q₀ + Δq)
!
! where Δq = (q₁ - q₀)/(α/2 + 1)
dq = (model%q1 - model%q0) / (model%alpha/2 + 1)
psia = 1/(2*pi) * phitedge / (model%q0 + dq)
end subroutine set_equil_an
@ -818,26 +838,70 @@ contains
ddpsidrr, ddpsidzz, ddpsidrz
! local variables
real(wp_) :: r_p ! poloidal radius
real(wp_) :: rho_p ! poloidal radius normalised
real(wp_) :: rho_t, rho_p ! √Φ_n, √ψ_n
real(wp_) :: dpsidphi ! (∂ψ_n/∂Φ_n)
real(wp_) :: ddpsidphidr, ddpsidphidz ! ∇(∂ψ_n/∂Φ_n)
real(wp_) :: dphidr, dphidz ! ∇Φ_n
real(wp_) :: q, dq ! q(ρ_p), Δq=(q₁-q₀)/(α/2 + 1)
real(wp_) :: dqdr, dqdz ! ∇q
if (iequil < 2) then
! Analytical model
! ρ_p is just the geometrical poloidal radius divided by a
r_p = hypot(R - model%R0, z - model%z0)
rho_p = r_p / model%a
! ψ_n = ρ_p²(R,z) = [(R-R₀)² + (z-z₀)²]/a²
! ρ_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
rho_p = frhopol(rho_t)
! ψ_n = ρ_p(ρ_t)²
if (present(psi_n)) psi_n = rho_p**2
! ∂ψ_n/∂R, ∂ψ_n/∂z
if (present(dpsidr)) dpsidr = 2*(R - model%R0)/model%a**2
if (present(dpsidz)) dpsidz = 2*(z - model%z0)/model%a**2
! 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/∂R², ∂²ψ_n/∂z², ∂²ψ_n/∂R∂z
if (present(ddpsidrr)) ddpsidrr = 2/model%a**2
if (present(ddpsidzz)) ddpsidzz = 2/model%a**2
if (present(ddpsidrz)) ddpsidrz = 0
! 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
! 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 * (q - model%q0)/rho_p**2 * dpsidphi * dphidr
dqdz = model%alpha/2 * (q - model%q0)/rho_p**2 * dpsidphi * dphidz
ddpsidphidr = - dpsidphi * dqdr/q
ddpsidphidz = - dpsidphi * dqdz/q
! Finally, ∇∇Φ_n = ∇∇ ρ_t²(R,z) = 2/a² I, so:
!
! ∇∇ψ_n = ∇(∂ψ_n/∂Φ_n) ∇Φ_n + (∂ψ_n/∂Φ_n) 2/a² I
!
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
else
! Numerical data
if (R <= rmxm .and. R >= rmnm .and. &