diff --git a/src/ElectromagneticFields.jl b/src/ElectromagneticFields.jl index 37c6cfb..6fea723 100644 --- a/src/ElectromagneticFields.jl +++ b/src/ElectromagneticFields.jl @@ -18,14 +18,22 @@ module ElectromagneticFields using .AnalyticCartesianField - export ABC, AxisymmetricTokamakCartesian, AxisymmetricTokamakCylindrical, - AxisymmetricTokamakToroidal, EzCosZ, - Solovev, SolovevXpoint, SolovevQuadratic, SymmetricQuadratic, ThetaPinch - export SolovevITER, SolovevNSTX, SolovevXpointITER, SolovevXpointNSTX + export ABC, EzCosZ, + AxisymmetricTokamakCartesian, + AxisymmetricTokamakCircular, + AxisymmetricTokamakCylindrical, + AxisymmetricTokamakToroidalRegularization, + Solovev, SolovevXpoint, SolovevQuadratic, + SymmetricQuadratic, ThetaPinch + + export SolovevITER, SolovevXpointITER, + SolovevNSTX, SolovevXpointNSTX + export @abc_equilibrium, @axisymmetric_tokamak_equilibrium_cartesian, + @axisymmetric_tokamak_equilibrium_circular, @axisymmetric_tokamak_equilibrium_cylindrical, - @axisymmetric_tokamak_equilibrium_toroidal, + @axisymmetric_tokamak_equilibrium_toroidal_regularisation, @ezcosz_perturbation, @solovev_equilibrium, @solovev_xpoint_equilibrium, @@ -35,8 +43,9 @@ module ElectromagneticFields include("analytic/abc.jl") include("analytic/axisymmetric_tokamak_cartesian.jl") + include("analytic/axisymmetric_tokamak_circular.jl") include("analytic/axisymmetric_tokamak_cylindrical.jl") - include("analytic/axisymmetric_tokamak_toroidal.jl") + include("analytic/axisymmetric_tokamak_toroidal_regularization.jl") include("analytic/ezcosz.jl") include("analytic/solovev_abstract.jl") include("analytic/solovev.jl") diff --git a/src/analytic/axisymmetric_tokamak_toroidal.jl b/src/analytic/axisymmetric_tokamak_circular.jl similarity index 63% rename from src/analytic/axisymmetric_tokamak_toroidal.jl rename to src/analytic/axisymmetric_tokamak_circular.jl index 2a159cf..1568efe 100644 --- a/src/analytic/axisymmetric_tokamak_toroidal.jl +++ b/src/analytic/axisymmetric_tokamak_circular.jl @@ -15,32 +15,32 @@ Parameters: * `B₀`: B-field at magnetic axis * `q₀`: safety factor at magnetic axis """ -module AxisymmetricTokamakToroidal +module AxisymmetricTokamakCircular import ..ElectromagneticFields import ..ElectromagneticFields: AnalyticEquilibrium, ZeroPerturbation import ..ElectromagneticFields: load_equilibrium, generate_equilibrium_code - export AxisymmetricTokamakToroidalEquilibrium + export AxisymmetricTokamakCircularEquilibrium const DEFAULT_R₀ = 1.0 const DEFAULT_B₀ = 1.0 const DEFAULT_q₀ = 2.0 - struct AxisymmetricTokamakToroidalEquilibrium{T <: Number} <: AnalyticEquilibrium + struct AxisymmetricTokamakCircularEquilibrium{T <: Number} <: AnalyticEquilibrium name::String R₀::T B₀::T q₀::T - function AxisymmetricTokamakToroidalEquilibrium{T}(R₀::T, B₀::T, q₀::T) where T <: Number + function AxisymmetricTokamakCircularEquilibrium{T}(R₀::T, B₀::T, q₀::T) where T <: Number new("AxisymmetricTokamakEquilibriumToroidal", R₀, B₀, q₀) end end - AxisymmetricTokamakToroidalEquilibrium(R₀::T=DEFAULT_R₀, B₀::T=DEFAULT_B₀, q₀::T=DEFAULT_q₀) where T <: Number = AxisymmetricTokamakToroidalEquilibrium{T}(R₀, B₀, q₀) + AxisymmetricTokamakCircularEquilibrium(R₀::T=DEFAULT_R₀, B₀::T=DEFAULT_B₀, q₀::T=DEFAULT_q₀) where T <: Number = AxisymmetricTokamakCircularEquilibrium{T}(R₀, B₀, q₀) - function Base.show(io::IO, equ::AxisymmetricTokamakToroidalEquilibrium) + function Base.show(io::IO, equ::AxisymmetricTokamakCircularEquilibrium) print(io, "Axisymmetric Tokamak Equilibrium in Toroidal Coordinates with\n") print(io, " R₀ = ", equ.R₀, "\n") print(io, " B₀ = ", equ.B₀, "\n") @@ -48,35 +48,35 @@ module AxisymmetricTokamakToroidal end - r(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = x[1] - θ(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = x[2] - ϕ(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = x[3] - R(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = equ.R₀ + r(x,equ) * cos(θ(x,equ)) - X(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = R(x,equ) * cos(ϕ(x,equ)) - Y(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = R(x,equ) * sin(ϕ(x,equ)) - Z(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = r(x,equ) * sin(θ(x,equ)) + r(x::AbstractVector, equ::AxisymmetricTokamakCircularEquilibrium) = x[1] + θ(x::AbstractVector, equ::AxisymmetricTokamakCircularEquilibrium) = x[2] + ϕ(x::AbstractVector, equ::AxisymmetricTokamakCircularEquilibrium) = x[3] + R(x::AbstractVector, equ::AxisymmetricTokamakCircularEquilibrium) = equ.R₀ + r(x,equ) * cos(θ(x,equ)) + X(x::AbstractVector, equ::AxisymmetricTokamakCircularEquilibrium) = R(x,equ) * cos(ϕ(x,equ)) + Y(x::AbstractVector, equ::AxisymmetricTokamakCircularEquilibrium) = R(x,equ) * sin(ϕ(x,equ)) + Z(x::AbstractVector, equ::AxisymmetricTokamakCircularEquilibrium) = r(x,equ) * sin(θ(x,equ)) - ElectromagneticFields.J(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = r(x,equ) * R(x,equ) + ElectromagneticFields.J(x::AbstractVector, equ::AxisymmetricTokamakCircularEquilibrium) = r(x,equ) * R(x,equ) - ElectromagneticFields.A₁(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = + equ.B₀ * equ.R₀ * ( Z(x,equ) / R(x,equ) * cos(θ(x,equ)) - log(R(x,equ) / equ.R₀) * sin(θ(x,equ)) ) / 2 - ElectromagneticFields.A₂(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = - equ.B₀ * equ.R₀ * ( Z(x,equ) / R(x,equ) * sin(θ(x,equ)) + log(R(x,equ) / equ.R₀) * cos(θ(x,equ)) ) * r(x,equ) / 2 - ElectromagneticFields.A₃(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = + equ.B₀ * r(x,equ)^2 / equ.q₀ / 2 + ElectromagneticFields.A₁(x::AbstractVector, equ::AxisymmetricTokamakCircularEquilibrium) = + equ.B₀ * equ.R₀ * ( Z(x,equ) / R(x,equ) * cos(θ(x,equ)) - log(R(x,equ) / equ.R₀) * sin(θ(x,equ)) ) / 2 + ElectromagneticFields.A₂(x::AbstractVector, equ::AxisymmetricTokamakCircularEquilibrium) = - equ.B₀ * equ.R₀ * ( Z(x,equ) / R(x,equ) * sin(θ(x,equ)) + log(R(x,equ) / equ.R₀) * cos(θ(x,equ)) ) * r(x,equ) / 2 + ElectromagneticFields.A₃(x::AbstractVector, equ::AxisymmetricTokamakCircularEquilibrium) = + equ.B₀ * r(x,equ)^2 / equ.q₀ / 2 - ElectromagneticFields.x¹(ξ::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = X(ξ,equ) - ElectromagneticFields.x²(ξ::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = Y(ξ,equ) - ElectromagneticFields.x³(ξ::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = Z(ξ,equ) + ElectromagneticFields.x¹(ξ::AbstractVector, equ::AxisymmetricTokamakCircularEquilibrium) = X(ξ,equ) + ElectromagneticFields.x²(ξ::AbstractVector, equ::AxisymmetricTokamakCircularEquilibrium) = Y(ξ,equ) + ElectromagneticFields.x³(ξ::AbstractVector, equ::AxisymmetricTokamakCircularEquilibrium) = Z(ξ,equ) - ElectromagneticFields.ξ¹(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = sqrt((sqrt(x[1]^2 + x[2]^2)-equ.R₀)^2 + x[3]^2) - ElectromagneticFields.ξ²(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = atan(x[3], sqrt(x[1]^2 + x[2]^2)-equ.R₀) - ElectromagneticFields.ξ³(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = atan(x[2], x[1]) + ElectromagneticFields.ξ¹(x::AbstractVector, equ::AxisymmetricTokamakCircularEquilibrium) = sqrt((sqrt(x[1]^2 + x[2]^2)-equ.R₀)^2 + x[3]^2) + ElectromagneticFields.ξ²(x::AbstractVector, equ::AxisymmetricTokamakCircularEquilibrium) = atan(x[3], sqrt(x[1]^2 + x[2]^2)-equ.R₀) + ElectromagneticFields.ξ³(x::AbstractVector, equ::AxisymmetricTokamakCircularEquilibrium) = atan(x[2], x[1]) - ElectromagneticFields.g₁₁(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = one(eltype(x)) - ElectromagneticFields.g₂₂(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = r(x, equ)^2 - ElectromagneticFields.g₃₃(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = R(x, equ)^2 + ElectromagneticFields.g₁₁(x::AbstractVector, equ::AxisymmetricTokamakCircularEquilibrium) = one(eltype(x)) + ElectromagneticFields.g₂₂(x::AbstractVector, equ::AxisymmetricTokamakCircularEquilibrium) = r(x, equ)^2 + ElectromagneticFields.g₃₃(x::AbstractVector, equ::AxisymmetricTokamakCircularEquilibrium) = R(x, equ)^2 - ElectromagneticFields.get_functions(::AxisymmetricTokamakToroidalEquilibrium) = (X=X, Y=Y, Z=Z, R=R, r=r, θ=θ, ϕ=ϕ) + ElectromagneticFields.get_functions(::AxisymmetricTokamakCircularEquilibrium) = (X=X, Y=Y, Z=Z, R=R, r=r, θ=θ, ϕ=ϕ) - function ElectromagneticFields.periodicity(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) + function ElectromagneticFields.periodicity(x::AbstractVector, equ::AxisymmetricTokamakCircularEquilibrium) p = zero(x) p[2] = 2π p[3] = 2π @@ -84,12 +84,12 @@ module AxisymmetricTokamakToroidal end macro axisymmetric_tokamak_equilibrium_toroidal(R₀, B₀, q₀) - generate_equilibrium_code(AxisymmetricTokamakToroidalEquilibrium(R₀, B₀, q₀); output=false) + generate_equilibrium_code(AxisymmetricTokamakCircularEquilibrium(R₀, B₀, q₀); output=false) end function init(R₀=DEFAULT_R₀, B₀=DEFAULT_B₀, q₀=DEFAULT_q₀; perturbation=ZeroPerturbation()) - equilibrium = AxisymmetricTokamakToroidalEquilibrium(R₀, B₀, q₀) - load_equilibrium(equilibrium, perturbation; target_module=AxisymmetricTokamakToroidal) + equilibrium = AxisymmetricTokamakCircularEquilibrium(R₀, B₀, q₀) + load_equilibrium(equilibrium, perturbation; target_module=AxisymmetricTokamakCircular) return equilibrium end diff --git a/src/analytic/axisymmetric_tokamak_toroidal_regularization.jl b/src/analytic/axisymmetric_tokamak_toroidal_regularization.jl new file mode 100644 index 0000000..1abae25 --- /dev/null +++ b/src/analytic/axisymmetric_tokamak_toroidal_regularization.jl @@ -0,0 +1,96 @@ +@doc raw""" +Axisymmetric tokamak equilibrium in (r,θ,ϕ) coordinates with covariant +components of the vector potential given by +```math +A (r, \theta, \phi) = B_0 \, \bigg( 0 , \, \frac{r R_0}{\cos (\theta)} - \bigg( \frac{R_0}{\cos (\theta)} \bigg)^2 \, \ln \bigg( \frac{R}{R_0} \bigg) , \, - \frac{r^2}{2 q_0} \bigg)^T , +``` +resulting in the magnetic field with covariant components +```math +B (r, \theta, \phi) = \frac{B_0}{q_0} \, \bigg( 0 , \, \frac{r^2}{R}, \, q_0 R_0 \bigg)^T , +``` +where $R = R_0 + r \cos \theta$. + +Parameters: + * `R₀`: position of magnetic axis + * `B₀`: B-field at magnetic axis + * `q₀`: safety factor at magnetic axis +""" +module AxisymmetricTokamakToroidalRegularization + + import ..ElectromagneticFields + import ..ElectromagneticFields: AnalyticEquilibrium, ZeroPerturbation + import ..ElectromagneticFields: load_equilibrium, generate_equilibrium_code + + export AxisymmetricTokamakToroidalRegularizationEquilibrium + + const DEFAULT_R₀ = 1.0 + const DEFAULT_B₀ = 1.0 + const DEFAULT_q₀ = 2.0 + + struct AxisymmetricTokamakToroidalRegularizationEquilibrium{T <: Number} <: AnalyticEquilibrium + name::String + R₀::T + B₀::T + q₀::T + + function AxisymmetricTokamakToroidalRegularizationEquilibrium{T}(R₀::T, B₀::T, q₀::T) where T <: Number + new("AxisymmetricTokamakEquilibriumToroidalRegularization", R₀, B₀, q₀) + end + end + + AxisymmetricTokamakToroidalRegularizationEquilibrium(R₀::T=DEFAULT_R₀, B₀::T=DEFAULT_B₀, q₀::T=DEFAULT_q₀) where T <: Number = AxisymmetricTokamakToroidalRegularizationEquilibrium{T}(R₀, B₀, q₀) + + function Base.show(io::IO, equ::AxisymmetricTokamakToroidalRegularizationEquilibrium) + print(io, "Axisymmetric Tokamak Equilibrium with Toroidal Regularization in Circular Coordinates with\n") + print(io, " R₀ = ", equ.R₀, "\n") + print(io, " B₀ = ", equ.B₀, "\n") + print(io, " q₀ = ", equ.q₀) + end + + + r(x::AbstractVector, equ::AxisymmetricTokamakToroidalRegularizationEquilibrium) = x[1] + θ(x::AbstractVector, equ::AxisymmetricTokamakToroidalRegularizationEquilibrium) = x[2] + ϕ(x::AbstractVector, equ::AxisymmetricTokamakToroidalRegularizationEquilibrium) = x[3] + R(x::AbstractVector, equ::AxisymmetricTokamakToroidalRegularizationEquilibrium) = equ.R₀ + r(x,equ) * cos(θ(x,equ)) + X(x::AbstractVector, equ::AxisymmetricTokamakToroidalRegularizationEquilibrium) = R(x,equ) * cos(ϕ(x,equ)) + Y(x::AbstractVector, equ::AxisymmetricTokamakToroidalRegularizationEquilibrium) = R(x,equ) * sin(ϕ(x,equ)) + Z(x::AbstractVector, equ::AxisymmetricTokamakToroidalRegularizationEquilibrium) = r(x,equ) * sin(θ(x,equ)) + + ElectromagneticFields.J(x::AbstractVector, equ::AxisymmetricTokamakToroidalRegularizationEquilibrium) = r(x,equ) * R(x,equ) + + ElectromagneticFields.A₁(x::AbstractVector, equ::AxisymmetricTokamakToroidalRegularizationEquilibrium) = zero(eltype(x)) + ElectromagneticFields.A₂(x::AbstractVector, equ::AxisymmetricTokamakToroidalRegularizationEquilibrium) = + equ.B₀ * equ.R₀ / cos(θ(x,equ))^2 * ( r(x,equ) * cos(θ(x,equ)) - equ.R₀ * log(R(x,equ) / equ.R₀) ) + ElectromagneticFields.A₃(x::AbstractVector, equ::AxisymmetricTokamakToroidalRegularizationEquilibrium) = + equ.B₀ * r(x,equ)^2 / equ.q₀ / 2 + + ElectromagneticFields.x¹(ξ::AbstractVector, equ::AxisymmetricTokamakToroidalRegularizationEquilibrium) = X(ξ,equ) + ElectromagneticFields.x²(ξ::AbstractVector, equ::AxisymmetricTokamakToroidalRegularizationEquilibrium) = Y(ξ,equ) + ElectromagneticFields.x³(ξ::AbstractVector, equ::AxisymmetricTokamakToroidalRegularizationEquilibrium) = Z(ξ,equ) + + ElectromagneticFields.ξ¹(x::AbstractVector, equ::AxisymmetricTokamakToroidalRegularizationEquilibrium) = sqrt((sqrt(x[1]^2 + x[2]^2)-equ.R₀)^2 + x[3]^2) + ElectromagneticFields.ξ²(x::AbstractVector, equ::AxisymmetricTokamakToroidalRegularizationEquilibrium) = atan(x[3], sqrt(x[1]^2 + x[2]^2)-equ.R₀) + ElectromagneticFields.ξ³(x::AbstractVector, equ::AxisymmetricTokamakToroidalRegularizationEquilibrium) = atan(x[2], x[1]) + + ElectromagneticFields.g₁₁(x::AbstractVector, equ::AxisymmetricTokamakToroidalRegularizationEquilibrium) = one(eltype(x)) + ElectromagneticFields.g₂₂(x::AbstractVector, equ::AxisymmetricTokamakToroidalRegularizationEquilibrium) = r(x, equ)^2 + ElectromagneticFields.g₃₃(x::AbstractVector, equ::AxisymmetricTokamakToroidalRegularizationEquilibrium) = R(x, equ)^2 + + ElectromagneticFields.get_functions(::AxisymmetricTokamakToroidalRegularizationEquilibrium) = (X=X, Y=Y, Z=Z, R=R, r=r, θ=θ, ϕ=ϕ) + + function ElectromagneticFields.periodicity(x::AbstractVector, ::AxisymmetricTokamakToroidalRegularizationEquilibrium) + p = zero(x) + p[2] = 2π + p[3] = 2π + return p + end + + macro axisymmetric_tokamak_equilibrium_toroidal_regularisation(R₀, B₀, q₀) + generate_equilibrium_code(AxisymmetricTokamakToroidalRegularizationEquilibrium(R₀, B₀, q₀); output=false) + end + + function init(R₀=DEFAULT_R₀, B₀=DEFAULT_B₀, q₀=DEFAULT_q₀; perturbation=ZeroPerturbation()) + equilibrium = AxisymmetricTokamakToroidalRegularizationEquilibrium(R₀, B₀, q₀) + load_equilibrium(equilibrium, perturbation; target_module=AxisymmetricTokamakToroidalRegularization) + return equilibrium + end + +end diff --git a/test/test_analytic.jl b/test/test_analytic.jl index 47aa354..e04737e 100644 --- a/test/test_analytic.jl +++ b/test/test_analytic.jl @@ -94,12 +94,13 @@ end # equilibrium list (equilibrium, parameters, periodicity, module) eqs = ( - (AxisymmetricTokamakCartesian, (2., 3., 2.), [0., 0., 0.]), - (AxisymmetricTokamakCylindrical, (2., 3., 2.), [0., 0., 2π]), - (AxisymmetricTokamakToroidal, (2., 3., 2.), [0., 2π, 2π]), - (SymmetricQuadratic, (1.), [0., 0., 0.]), - (ThetaPinch, (1.), [0., 0., 0.]), - (ABC, (1., 0.5, 0.5), [0., 0., 0.]), + (AxisymmetricTokamakCartesian, (2., 3., 2.), [0., 0., 0.]), + (AxisymmetricTokamakCircular, (2., 3., 2.), [0., 2π, 2π]), + (AxisymmetricTokamakCylindrical, (2., 3., 2.), [0., 0., 2π]), + (AxisymmetricTokamakToroidalRegularization, (2., 3., 2.), [0., 2π, 2π]), + (SymmetricQuadratic, (1.), [0., 0., 0.]), + (ThetaPinch, (1.), [0., 0., 0.]), + (ABC, (1., 0.5, 0.5), [0., 0., 0.]), (Solovev, (6.2, 5.3, 0.32, 1.8, 0.45, -0.155), [0., 0., 2π]), (SolovevXpoint, (6.2, 5.3, 0.32, 1.8, 0.45, -0.155, 0.88, -0.60), [0., 0., 2π]), (SolovevQuadratic, (6.2, 5.3, 1., 1.), [0., 0., 2π]), @@ -148,6 +149,16 @@ function test_axisymmetric_tokamak_cartesian_equilibrium(equ_mod, t=0., x=[1.5, @test equ_mod.B₃(t,x) == + equ_mod.B₀ / equ_mod.q₀ * (equ_mod.R(t,x) - equ_mod.R₀) / equ_mod.R(t,x) end +function test_axisymmetric_tokamak_circular_equilibrium(equ_mod, t=0., x=[0.5, π/10, π/5]) + @test equ_mod.B¹(t,x) == 0 + @test equ_mod.B²(t,x) == - equ_mod.B₀ / equ_mod.q₀ / equ_mod.R(t,x) + @test equ_mod.B³(t,x) ≈ - equ_mod.B₀ * equ_mod.R₀ / equ_mod.R(t,x)^2 atol=1E-14 + + @test equ_mod.B₁(t,x) == 0 + @test equ_mod.B₂(t,x) == - equ_mod.B₀ / equ_mod.q₀ * equ_mod.r(t,x)^2 / equ_mod.R(t,x) + @test equ_mod.B₃(t,x) ≈ - equ_mod.B₀ * equ_mod.R₀ atol=1E-14 +end + function test_axisymmetric_tokamak_cylindrical_equilibrium(equ_mod, t=0., x=[1.5, 0.5, π/5]) @test equ_mod.B¹(t,x) == + equ_mod.B₀ / equ_mod.q₀ * equ_mod.Z(t,x) / equ_mod.R(t,x) @test equ_mod.B²(t,x) == - equ_mod.B₀ / equ_mod.q₀ * (equ_mod.R(t,x) - equ_mod.R₀) / equ_mod.R(t,x) @@ -158,14 +169,17 @@ function test_axisymmetric_tokamak_cylindrical_equilibrium(equ_mod, t=0., x=[1.5 @test equ_mod.B₃(t,x) == - equ_mod.B₀ * equ_mod.R₀ end -function test_axisymmetric_tokamak_toroidal_equilibrium(equ_mod, t=0., x=[0.5, π/10, π/5]) - @test equ_mod.B¹(t,x) == 0 - @test equ_mod.B²(t,x) == - equ_mod.B₀ / equ_mod.q₀ / equ_mod.R(t,x) - @test equ_mod.B³(t,x) ≈ - equ_mod.B₀ * equ_mod.R₀ / equ_mod.R(t,x)^2 atol=1E-14 +function test_consistency_axisymmetric_tokamak_circular_equilibrium(equ_tor, equ_car, t=0., ξ=[0.5, π/10, π/5]) + x = equ_tor.to_cartesian(t,ξ) + DF = equ_tor.DF(t,ξ) - @test equ_mod.B₁(t,x) == 0 - @test equ_mod.B₂(t,x) == - equ_mod.B₀ / equ_mod.q₀ * equ_mod.r(t,x)^2 / equ_mod.R(t,x) - @test equ_mod.B₃(t,x) ≈ - equ_mod.B₀ * equ_mod.R₀ atol=1E-14 + B_tor = [equ_tor.B¹(t,ξ), equ_tor.B²(t,ξ), equ_tor.B³(t,ξ)] + B_car = [equ_car.B¹(t,x), equ_car.B²(t,x), equ_car.B³(t,x)] + + B̂_tor = [equ_tor.B₁(t,ξ), equ_tor.B₂(t,ξ), equ_tor.B₃(t,ξ)] + B̂_car = [equ_car.B₁(t,x), equ_car.B₂(t,x), equ_car.B₃(t,x)] + + @test B_tor' * B̂_tor ≈ B_car' * B̂_car atol=1E-12 end function test_consistency_axisymmetric_tokamak_cylindrical_equilibrium(equ_cyl, equ_car, t=0., ξ=[1.5, 0.5, π/5]) @@ -182,19 +196,6 @@ function test_consistency_axisymmetric_tokamak_cylindrical_equilibrium(equ_cyl, @test B_cyl' * B̂_cyl ≈ B_car' * B̂_car atol=1E-12 end -function test_consistency_axisymmetric_tokamak_toroidal_equilibrium(equ_tor, equ_car, t=0., ξ=[0.5, π/10, π/5]) - x = equ_tor.to_cartesian(t,ξ) - DF = equ_tor.DF(t,ξ) - - B_tor = [equ_tor.B¹(t,ξ), equ_tor.B²(t,ξ), equ_tor.B³(t,ξ)] - B_car = [equ_car.B¹(t,x), equ_car.B²(t,x), equ_car.B³(t,x)] - - B̂_tor = [equ_tor.B₁(t,ξ), equ_tor.B₂(t,ξ), equ_tor.B₃(t,ξ)] - B̂_car = [equ_car.B₁(t,x), equ_car.B₂(t,x), equ_car.B₃(t,x)] - - @test B_tor' * B̂_tor ≈ B_car' * B̂_car atol=1E-12 -end - function test_symmetric_quadratic_equilibrium(equ_mod, t=0., x=[1.0, 0.5, 0.5]) @test equ_mod.B¹(t,x) == equ_mod.B₁(t,x) @test equ_mod.B²(t,x) == equ_mod.B₂(t,x) @@ -251,16 +252,16 @@ end @testset "$(rpad("Magnetic Fields",60))" begin test_axisymmetric_tokamak_cartesian_equilibrium(AxisymmetricTokamakCartesian) + test_axisymmetric_tokamak_circular_equilibrium(AxisymmetricTokamakCircular) test_axisymmetric_tokamak_cylindrical_equilibrium(AxisymmetricTokamakCylindrical) - test_axisymmetric_tokamak_toroidal_equilibrium(AxisymmetricTokamakToroidal) test_symmetric_quadratic_equilibrium(SymmetricQuadratic) test_theta_pinch_equilibrium(ThetaPinch) test_abc_equilibrium(ABC) end @testset "$(rpad("Consistency",60))" begin + test_consistency_axisymmetric_tokamak_circular_equilibrium(AxisymmetricTokamakCircular, AxisymmetricTokamakCartesian) test_consistency_axisymmetric_tokamak_cylindrical_equilibrium(AxisymmetricTokamakCylindrical, AxisymmetricTokamakCartesian) - test_consistency_axisymmetric_tokamak_toroidal_equilibrium(AxisymmetricTokamakToroidal, AxisymmetricTokamakCartesian) end println()