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

make the sign of q consistent with cocos=3

Browse Source
This commit is contained in:
Michele Guerini Rocco 2024-10-18 13:40:03 +02:00 committed by rnhmjoj
parent fdf5ef72fe
commit 8c144c3892
Signed by: rnhmjoj
GPG Key ID: BFBAF4C975F76450
3 changed files with 19 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
src/gray_equil.f90
View File

@ -81,7 +81,7 @@ module gray_equil
R_hi, z_hi, R_lo, z_lo) R_hi, z_hi, R_lo, z_lo)
! Computes a contour of the ψ(R,z)=ψ₀ flux surface. ! Computes a contour of the ψ(R,z)=ψ₀ flux surface.
! Notes: ! Notes:
! - R,z are the contour points ! - R,z are the contour points ordered clockwise
! - R_min is a value such that R>R_min for any contour point ! - R_min is a value such that R>R_min for any contour point
! - (R,z)_hi and (R,z)_lo are a guess for the higher and lower ! - (R,z)_hi and (R,z)_lo are a guess for the higher and lower
! horizontal-tangent points of the contour. These variables ! horizontal-tangent points of the contour. These variables
@ -435,7 +435,7 @@ contains
! q(ρ_p) = q₀ + (q₁-q₀) ρ_p^α ! q(ρ_p) = q₀ + (q₁-q₀) ρ_p^α
! !
rho_p = sqrt(psi_n) rho_p = sqrt(psi_n)
q = abs(self%q0 + (self%q1 - self%q0) * rho_p**self%alpha) q = self%q0 + (self%q1 - self%q0) * rho_p**self%alpha
end function analytic_safety end function analytic_safety
@ -541,8 +541,8 @@ contains
theta = 2*pi / (n - 1) theta = 2*pi / (n - 1)
do concurrent (i=1:n) do concurrent (i=1:n)
R(i) = self%R0 + r_g * cos(theta*(i-1)) R(i) = self%R0 + r_g * cos(-theta*(i-1))
z(i) = self%z0 + r_g * sin(theta*(i-1)) z(i) = self%z0 + r_g * sin(-theta*(i-1))
end do end do
end subroutine analytic_flux_contour end subroutine analytic_flux_contour
@ -955,9 +955,10 @@ contains
real(wp_) :: dx real(wp_) :: dx
integer :: k integer :: k
call self%q_spline%init(data%psi, abs(data%q), npsi) call self%q_spline%init(data%psi, data%q, npsi)
! Toroidal flux as Φ(ψ) = 2π ∫q(ψ)dψ ! Toroidal flux as Φ(ψ) = 2π ∫ q(ψ)dψ
! Specifically, Φ(ψ_n) = 2π|ψ_a| ∫₀ψ_n q(ψ_n)dψ_n, with ψ_n ∈ [0,1]
phi(1) = 0 phi(1) = 0
do k = 1, npsi-1 do k = 1, npsi-1
dx = self%q_spline%data(k+1) - self%q_spline%data(k) dx = self%q_spline%data(k+1) - self%q_spline%data(k)

4
src/gray_tables.f90
View File

@ -366,10 +366,10 @@ contains
if (flag == 1) then if (flag == 1) then
! return f(x) ! return f(x)
f(1) = equil%safety(x(1)) - qvals(i) f(1) = abs(equil%safety(x(1))) - qvals(i)
else else
! return f'(x), computed numerically ! return f'(x), computed numerically
df(1,1) = (equil%safety(x(1) + e) - equil%safety(x(1) - e)) / (2*e) df(1,1) = (abs(equil%safety(x(1) + e)) - abs(equil%safety(x(1) - e))) / (2*e)
end if end if
end subroutine end subroutine

16
src/magsurf_data.f90
View File

@ -108,7 +108,7 @@ contains
block block
real(wp_) :: a, b real(wp_) :: a, b
call equil%pol_flux(equil%axis(1), equil%axis(2), ddpsidrr=a, ddpsidzz=b) call equil%pol_flux(equil%axis(1), equil%axis(2), ddpsidrr=a, ddpsidzz=b)
q(1) = abs(equil%b_axis) / sqrt(a*b)/equil%psi_a q(1) = equil%b_axis / sqrt(a*b)/equil%psi_a
end block end block
! Compute the uniform λ grid and H(0, λ) ! Compute the uniform λ grid and H(0, λ)
@ -140,8 +140,8 @@ contains
block block
! Variables for the previous and current contour point ! Variables for the previous and current contour point
real(wp_) :: B_tot0, B_pol0, B_tot1, B_pol1, J_phi0, J_phi1 real(wp_) :: B_tot0, B_tot1, B_pol0, B_pol1, J_phi0, J_phi1
real(wp_) :: dl ! line element real(wp_) :: dl, dir ! line element, sgn(dl⋅B_p)
! Initialise all integrals ! Initialise all integrals
area(i) = 0 area(i) = 0
@ -184,6 +184,10 @@ contains
! The line differential is difference of the two vertices ! The line differential is difference of the two vertices
dl = norm2(P - Q) dl = norm2(P - Q)
! We always assume dl ∥ B_p (numerically they're not),
! but compute the sign of dl⋅B_p for the case I_p < 0.
dir = sign(one, dot_product([B_R, B_z], P - Q))
end block end block
! Evaluate integrands at the current point for line integrals ∮ dl ! Evaluate integrands at the current point for line integrals ∮ dl
@ -197,18 +201,18 @@ contains
! Do one step of the trapezoid method ! Do one step of the trapezoid method
! N = ∮ dl/B_p ! N = ∮ dl/B_p
! dA/dψ = ∮ (1/R) dl/B_p ! dA/dψ = ∮ (1/R) dl/B_p
! I_p = 1/μ₀ ∫ B_p⋅dS_φ = 1/μ₀ ∮ B_p dl
! ⟨B⟩ = 1/N ∮ B dl/B_p ! ⟨B⟩ = 1/N ∮ B dl/B_p
! ⟨B²⟩ = 1/N ∮ B² dl/B_p ! ⟨B²⟩ = 1/N ∮ B² dl/B_p
! ⟨R⁻²⟩ = 1/N ∮ R⁻² dl/B_p ! ⟨R⁻²⟩ = 1/N ∮ R⁻² dl/B_p
! ⟨J_φ⟩ = 1/N ∮ J_φ dl/B_p ! ⟨J_φ⟩ = 1/N ∮ J_φ dl/B_p
! I_p = 1/μ₀ ∮ B_p⋅dl
norm = norm + dl*(1/B_pol1 + 1/B_pol0)/2 norm = norm + dl*(1/B_pol1 + 1/B_pol0)/2
dAdpsi = dAdpsi + dl*(1/(B_pol1*R(j+1)) + 1/(B_pol0*R(j)))/2 dAdpsi = dAdpsi + dl*(1/(B_pol1*R(j+1)) + 1/(B_pol0*R(j)))/2
I_p(i) = I_p(i) + dl*(B_pol1 + B_pol0)/2 * mu0inv
B_avg(i) = B_avg(i) + dl*(B_tot1/B_pol1 + B_tot0/B_pol0)/2 B_avg(i) = B_avg(i) + dl*(B_tot1/B_pol1 + B_tot0/B_pol0)/2
B2_avg = B2_avg + dl*(B_tot1**2/B_pol1 + B_tot0**2/B_pol0)/2 B2_avg = B2_avg + dl*(B_tot1**2/B_pol1 + B_tot0**2/B_pol0)/2
R2_inv_avg = R2_inv_avg + dl*(1/(B_pol1*R(j+1)**2) + 1/(B_pol0*R(j)**2))/2 R2_inv_avg = R2_inv_avg + dl*(1/(B_pol1*R(j+1)**2) + 1/(B_pol0*R(j)**2))/2
J_phi_avg(i) = J_phi_avg(i) + dl*(J_phi0/(B_pol0*R(j)) + J_phi1/(B_pol1*R(j+1)))/2 J_phi_avg(i) = J_phi_avg(i) + dl*(J_phi0/(B_pol0*R(j)) + J_phi1/(B_pol1*R(j+1)))/2
I_p(i) = I_p(i) + dl*(B_pol1 + B_pol0)/2 * dir * mu0inv
! Update the previous values ! Update the previous values
J_phi0 = J_phi1 J_phi0 = J_phi1
@ -255,7 +259,7 @@ contains
! !
! q(ρ) = F/2π ∮ 1/R² dl/B_p = F/2π N ⟨R⁻²⟩ ! q(ρ) = F/2π ∮ 1/R² dl/B_p = F/2π N ⟨R⁻²⟩
! !
q(i) = abs(fpol/(2*pi) * norm * R2_inv_avg) q(i) = fpol/(2*pi) * norm * R2_inv_avg
! Using Φ=Φ_a⋅ρ_t² and q=1/2π⋅∂Φ/∂ψ (see gray_equil.f90), we have ! Using Φ=Φ_a⋅ρ_t² and q=1/2π⋅∂Φ/∂ψ (see gray_equil.f90), we have
! !