diff --git a/src/ElectromagneticFields.jl b/src/ElectromagneticFields.jl index 05bc86d..559191d 100644 --- a/src/ElectromagneticFields.jl +++ b/src/ElectromagneticFields.jl @@ -9,13 +9,13 @@ module ElectromagneticFields export AnalyticField, AnalyticEquilibrium, AnalyticPerturbation, ZeroPerturbation - export CartesianField, CartesianEquilibrium, CartesianPerturbation - export load_equilibrium, periodicity include("analytic/analytic_equilibrium.jl") include("analytic/cartesian_equilibrium.jl") + using .AnalyticCartesianField + export ABC, AxisymmetricTokamakCartesian, AxisymmetricTokamakCylindrical, AxisymmetricTokamakToroidal, EzCosZ, Solovev, SolovevXpoint, SolovevQuadratic, SymmetricQuadratic, ThetaPinch @@ -36,8 +36,7 @@ module ElectromagneticFields include("analytic/axisymmetric_tokamak_cylindrical.jl") include("analytic/axisymmetric_tokamak_toroidal.jl") include("analytic/ezcosz.jl") - include("analytic/solovev_common.jl") - include("analytic/solovev_psi.jl") + include("analytic/solovev_abstract.jl") include("analytic/solovev.jl") include("analytic/solovev_quadratic.jl") include("analytic/solovev_xpoint.jl") diff --git a/src/analytic/abc.jl b/src/analytic/abc.jl index da7f2e5..009e06f 100644 --- a/src/analytic/abc.jl +++ b/src/analytic/abc.jl @@ -1,5 +1,5 @@ @doc raw""" -ABC equilibrium in (x,y,z) coordinates with covariant components of the vector +ABCEquilibrium equilibrium in (x,y,z) coordinates with covariant components of the vector potential given by ```math A (x,y,z) = \big( a \, \sin(z) + c \, \cos(y) , \, b \, \sin(x) + a \, \cos(z) , \, c \, \sin(y) + b \, \cos(x) \big)^T @@ -8,32 +8,49 @@ resulting in the magnetic field ``B(x,y,z) = A(x,y,z)``. Parameters: `a`, `b`, `c` """ -struct ABC{T <: Number} <: CartesianEquilibrium - name::String - a::T - b::T - c::T +module ABC - ABC{T}(a::T, b::T, c::T) where T <: Number = new("ABCEquilibrium", a, b, c) -end + import ..ElectromagneticFields + import ..ElectromagneticFields: CartesianEquilibrium, ZeroPerturbation + import ..ElectromagneticFields: load_equilibrium, generate_equilibrium_code + import ..AnalyticCartesianField: X, Y, Z -ABC(a::T, b::T, c::T) where T <: Number = ABC{T}(a, b, c) + export ABCEquilibrium + struct ABCEquilibrium{T <: Number} <: CartesianEquilibrium + name::String + a::T + b::T + c::T -function Base.show(io::IO, equ::ABC) - print(io, "ABC Equilibrium with\n") - print(io, " A = ", equ.a, "\n") - print(io, " B = ", equ.b, "\n") - print(io, " C = ", equ.c) -end + ABCEquilibrium{T}(a::T, b::T, c::T) where T <: Number = new("ABCEquilibrium", a, b, c) + end + + ABCEquilibrium(a::T, b::T, c::T) where T <: Number = ABCEquilibrium{T}(a, b, c) + + + function Base.show(io::IO, equ::ABCEquilibrium) + print(io, "ABC Equilibrium with\n") + print(io, " A = ", equ.a, "\n") + print(io, " B = ", equ.b, "\n") + print(io, " C = ", equ.c) + end + + + ElectromagneticFields.A₁(x::AbstractVector, equ::ABCEquilibrium) = equ.a * sin(x[3]) + equ.c * cos(x[2]) + ElectromagneticFields.A₂(x::AbstractVector, equ::ABCEquilibrium) = equ.b * sin(x[1]) + equ.a * cos(x[3]) + ElectromagneticFields.A₃(x::AbstractVector, equ::ABCEquilibrium) = equ.c * sin(x[2]) + equ.b * cos(x[1]) + ElectromagneticFields.get_functions(::ABCEquilibrium) = (X=X, Y=Y, Z=Z) -A₁(x::AbstractVector, equ::ABC) = equ.a * sin(x[3]) + equ.c * cos(x[2]) -A₂(x::AbstractVector, equ::ABC) = equ.b * sin(x[1]) + equ.a * cos(x[3]) -A₃(x::AbstractVector, equ::ABC) = equ.c * sin(x[2]) + equ.b * cos(x[1]) + macro abc_equilibrium(a, b, c) + generate_equilibrium_code(ABCEquilibrium(a, b, c); output=false) + end -get_functions(::ABC) = (X=X, Y=Y, Z=Z) + function init(a, b, c; perturbation=ZeroPerturbation()) + equilibrium = ABCEquilibrium(a, b, c) + load_equilibrium(equilibrium, perturbation; target_module=ABC) + return equilibrium + end -macro abc_equilibrium(a, b, c) - generate_equilibrium_code(ABC(a, b, c); output=false) end diff --git a/src/analytic/analytic_equilibrium.jl b/src/analytic/analytic_equilibrium.jl index 490a897..0aad64a 100644 --- a/src/analytic/analytic_equilibrium.jl +++ b/src/analytic/analytic_equilibrium.jl @@ -8,15 +8,15 @@ abstract type AnalyticPerturbation <: AnalyticField end get_functions(::AnalyticField) = (;) -x¹(::AbstractVector{T}, ::AnalyticField) where {T} = error("x¹() not implemented for ", ET) -x²(::AbstractVector{T}, ::AnalyticField) where {T} = error("x²() not implemented for ", ET) -x³(::AbstractVector{T}, ::AnalyticField) where {T} = error("x³() not implemented for ", ET) +x¹(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("x¹() not implemented for ", ET) +x²(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("x²() not implemented for ", ET) +x³(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("x³() not implemented for ", ET) -ξ¹(::AbstractVector{T}, ::AnalyticField) where {T} = error("ξ¹() not implemented for ", ET) -ξ²(::AbstractVector{T}, ::AnalyticField) where {T} = error("ξ²() not implemented for ", ET) -ξ³(::AbstractVector{T}, ::AnalyticField) where {T} = error("ξ³() not implemented for ", ET) +ξ¹(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("ξ¹() not implemented for ", ET) +ξ²(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("ξ²() not implemented for ", ET) +ξ³(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("ξ³() not implemented for ", ET) -J(::AbstractVector, ::AnalyticField) = error("J() not implemented for ", ET) +J(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("J() not implemented for ", ET) g₁₁(::AbstractVector{T}, ::AnalyticField) where {T} = one(T) g₁₂(::AbstractVector{T}, ::AnalyticField) where {T} = zero(T) @@ -33,9 +33,9 @@ periodicity(x::AbstractVector, ::AnalyticField) = zero(x) from_cartesian(x::AbstractVector, equ::AnalyticField) = [ξ¹(x,equ), ξ²(x,equ), ξ³(x,equ)] to_cartesian(ξ::AbstractVector, equ::AnalyticField) = [x¹(ξ,equ), x²(ξ,equ), x³(ξ,equ)] -A₁(::AbstractVector{T}, ::AnalyticField) where {T} = zero(T) -A₂(::AbstractVector{T}, ::AnalyticField) where {T} = zero(T) -A₃(::AbstractVector{T}, ::AnalyticField) where {T} = zero(T) +A₁(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("A₁() not implemented for ", ET) +A₂(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("A₂() not implemented for ", ET) +A₃(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("A₃() not implemented for ", ET) φ(::AbstractVector{T}, ::AnalyticField) where {T} = zero(T) @@ -55,9 +55,9 @@ struct ZeroPerturbation <: AnalyticPerturbation ZeroPerturbation() = new("ZeroPerturbation") end -A₁(::AbstractVector{T}, ::ZeroPerturbation) where {T} = zero(T) -A₂(::AbstractVector{T}, ::ZeroPerturbation) where {T} = zero(T) -A₃(::AbstractVector{T}, ::ZeroPerturbation) where {T} = zero(T) +A₁(::AbstractVector{T}, ::AnalyticPerturbation) where {T} = zero(T) +A₂(::AbstractVector{T}, ::AnalyticPerturbation) where {T} = zero(T) +A₃(::AbstractVector{T}, ::AnalyticPerturbation) where {T} = zero(T) "Returns the i-th component of the vector corresponding to the one-form α" diff --git a/src/analytic/axisymmetric_tokamak_cartesian.jl b/src/analytic/axisymmetric_tokamak_cartesian.jl index 1bf3473..e16ba7a 100644 --- a/src/analytic/axisymmetric_tokamak_cartesian.jl +++ b/src/analytic/axisymmetric_tokamak_cartesian.jl @@ -11,45 +11,65 @@ B (x,y,z) = \frac{B_0}{q_0} \, \bigg( - \frac{q_0 R_0 y + x z}{R^2} , \, \frac{q where $R = \sqrt{ x^2 + y^2 }$ and $r = \sqrt{ (R - R_0)^2 + z^2 }$. Parameters: - * `R₀`: position of magnetic axis - * `B₀`: B-field at magnetic axis - * `q₀`: safety factor at magnetic axis +* `R₀`: position of magnetic axis +* `B₀`: B-field at magnetic axis +* `q₀`: safety factor at magnetic axis """ -struct AxisymmetricTokamakCartesian{T <: Number} <: CartesianEquilibrium - name::String - R₀::T - B₀::T - q₀::T - - function AxisymmetricTokamakCartesian{T}(R₀::T, B₀::T, q₀::T) where T <: Number - new("AxisymmetricTokamakCartesianEquilibrium", R₀, B₀, q₀) +module AxisymmetricTokamakCartesian + + import ..ElectromagneticFields + import ..ElectromagneticFields: CartesianEquilibrium, ZeroPerturbation + import ..ElectromagneticFields: load_equilibrium, generate_equilibrium_code + import ..AnalyticCartesianField: X, Y, Z + + export AxisymmetricTokamakCartesianEquilibrium + + const DEFAULT_R₀ = 1.0 + const DEFAULT_B₀ = 1.0 + const DEFAULT_q₀ = 2.0 + + struct AxisymmetricTokamakCartesianEquilibrium{T <: Number} <: CartesianEquilibrium + name::String + R₀::T + B₀::T + q₀::T + + function AxisymmetricTokamakCartesianEquilibrium{T}(R₀::T, B₀::T, q₀::T) where T <: Number + new("AxisymmetricTokamakCartesianEquilibrium", R₀, B₀, q₀) + end end -end -AxisymmetricTokamakCartesian(R₀::T=1.0, B₀::T=1.0, q₀::T=2.0) where T <: Number = AxisymmetricTokamakCartesian{T}(R₀, B₀, q₀) + AxisymmetricTokamakCartesianEquilibrium(R₀::T=DEFAULT_R₀, B₀::T=DEFAULT_B₀, q₀::T=DEFAULT_q₀) where T <: Number = AxisymmetricTokamakCartesianEquilibrium{T}(R₀, B₀, q₀) + function Base.show(io::IO, equ::AxisymmetricTokamakCartesianEquilibrium) + print(io, "Axisymmetric Tokamak Equilibrium in (x,y,z) Coordinates with\n") + print(io, " R₀ = ", equ.R₀, "\n") + print(io, " B₀ = ", equ.B₀, "\n") + print(io, " q₀ = ", equ.q₀) + end -function Base.show(io::IO, equ::AxisymmetricTokamakCartesian) - print(io, "Axisymmetric Tokamak Equilibrium in (x,y,z) 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::AxisymmetricTokamakCartesianEquilibrium) = X(x,equ)^2 + Y(x,equ)^2 + r²(x::AbstractVector, equ::AxisymmetricTokamakCartesianEquilibrium) = (R(x,equ) - equ.R₀)^2 + Z(x,equ)^2 + R(x::AbstractVector, equ::AxisymmetricTokamakCartesianEquilibrium) = sqrt(R²(x,equ)) + r(x::AbstractVector, equ::AxisymmetricTokamakCartesianEquilibrium) = sqrt(r²(x,equ)) + θ(x::AbstractVector, equ::AxisymmetricTokamakCartesianEquilibrium) = atan(Z(x,equ), R(x,equ) - equ.R₀) + ϕ(x::AbstractVector, equ::AxisymmetricTokamakCartesianEquilibrium) = atan(Y(x,equ), X(x,equ)) -R²(x::AbstractVector, equ::AxisymmetricTokamakCartesian) = X(x,equ)^2 + Y(x,equ)^2 -r²(x::AbstractVector, equ::AxisymmetricTokamakCartesian) = (R(x,equ) - equ.R₀)^2 + Z(x,equ)^2 -R(x::AbstractVector, equ::AxisymmetricTokamakCartesian) = sqrt(R²(x,equ)) -r(x::AbstractVector, equ::AxisymmetricTokamakCartesian) = sqrt(r²(x,equ)) -θ(x::AbstractVector, equ::AxisymmetricTokamakCartesian) = atan(Z(x,equ), R(x,equ) - equ.R₀) -ϕ(x::AbstractVector, equ::AxisymmetricTokamakCartesian) = atan(Y(x,equ), X(x,equ)) + ElectromagneticFields.A₁(x::AbstractVector, equ::AxisymmetricTokamakCartesianEquilibrium) = + equ.B₀ * (equ.R₀ * X(x,equ) * Z(x,equ) - r²(x,equ) * Y(x,equ) / equ.q₀ ) / R²(x,equ) / 2 + ElectromagneticFields.A₂(x::AbstractVector, equ::AxisymmetricTokamakCartesianEquilibrium) = + equ.B₀ * (equ.R₀ * Y(x,equ) * Z(x,equ) + r²(x,equ) * X(x,equ) / equ.q₀ ) / R²(x,equ) / 2 + ElectromagneticFields.A₃(x::AbstractVector, equ::AxisymmetricTokamakCartesianEquilibrium) = - equ.B₀ * equ.R₀ * log(R(x,equ) / equ.R₀) / 2 -A₁(x::AbstractVector, equ::AxisymmetricTokamakCartesian) = + equ.B₀ * (equ.R₀ * X(x,equ) * Z(x,equ) - r²(x,equ) * Y(x,equ) / equ.q₀ ) / R²(x,equ) / 2 -A₂(x::AbstractVector, equ::AxisymmetricTokamakCartesian) = + equ.B₀ * (equ.R₀ * Y(x,equ) * Z(x,equ) + r²(x,equ) * X(x,equ) / equ.q₀ ) / R²(x,equ) / 2 -A₃(x::AbstractVector, equ::AxisymmetricTokamakCartesian) = - equ.B₀ * equ.R₀ * log(R(x,equ) / equ.R₀) / 2 + ElectromagneticFields.get_functions(::AxisymmetricTokamakCartesianEquilibrium) = (X=X, Y=Y, Z=Z, R=R, r=r, θ=θ, ϕ=ϕ, R²=R², r²=r²) -get_functions(::AxisymmetricTokamakCartesian) = (X=X, Y=Y, Z=Z, R=R, r=r, θ=θ, ϕ=ϕ, R²=R², r²=r²) + macro axisymmetric_tokamak_equilibrium_cartesian(R₀, B₀, q₀) + generate_equilibrium_code(AxisymmetricTokamakCartesianEquilibrium(R₀, B₀, q₀); output=false) + end + + function init(R₀=DEFAULT_R₀, B₀=DEFAULT_B₀, q₀=DEFAULT_q₀; perturbation=ZeroPerturbation()) + equilibrium = AxisymmetricTokamakCartesianEquilibrium(R₀, B₀, q₀) + load_equilibrium(equilibrium, perturbation; target_module=AxisymmetricTokamakCartesian) + return equilibrium + end -macro axisymmetric_tokamak_equilibrium_cartesian(R₀, B₀, q₀) - generate_equilibrium_code(AxisymmetricTokamakCartesian(R₀, B₀, q₀); output=false) end diff --git a/src/analytic/axisymmetric_tokamak_cylindrical.jl b/src/analytic/axisymmetric_tokamak_cylindrical.jl index 9e6c23f..b02692f 100644 --- a/src/analytic/axisymmetric_tokamak_cylindrical.jl +++ b/src/analytic/axisymmetric_tokamak_cylindrical.jl @@ -1,5 +1,3 @@ -using SymEngine - @doc raw""" Axisymmetric tokamak equilibrium in (R,Z,ϕ) coordinates with covariant components of the vector potential given by @@ -17,65 +15,82 @@ Parameters: * `B₀`: B-field at magnetic axis * `q₀`: safety factor at magnetic axis """ -struct AxisymmetricTokamakCylindrical{T <: Number} <: AnalyticEquilibrium - name::String - R₀::T - B₀::T - q₀::T - - function AxisymmetricTokamakCylindrical{T}(R₀::T, B₀::T, q₀::T) where T <: Number - new("AxisymmetricTokamakCylindricalEquilibrium", R₀, B₀, q₀) - end -end +module AxisymmetricTokamakCylindrical -AxisymmetricTokamakCylindrical(R₀::T=1.0, B₀::T=1.0, q₀::T=2.0) where T <: Number = AxisymmetricTokamakCylindrical{T}(R₀, B₀, q₀) + import ..ElectromagneticFields + import ..ElectromagneticFields: AnalyticEquilibrium, ZeroPerturbation + import ..ElectromagneticFields: load_equilibrium, generate_equilibrium_code + export AxisymmetricTokamakCylindricalEquilibrium -function Base.show(io::IO, equ::AxisymmetricTokamakCylindrical) - print(io, "Axisymmetric Tokamak Equilibrium in (R,Z,ϕ) Coordinates with\n") - print(io, " R₀ = ", equ.R₀, "\n") - print(io, " B₀ = ", equ.B₀, "\n") - print(io, " q₀ = ", equ.q₀) -end + const DEFAULT_R₀ = 1.0 + const DEFAULT_B₀ = 1.0 + const DEFAULT_q₀ = 2.0 + struct AxisymmetricTokamakCylindricalEquilibrium{T <: Number} <: AnalyticEquilibrium + name::String + R₀::T + B₀::T + q₀::T -R(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = x[1] -Z(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = x[2] -ϕ(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = x[3] -r²(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = (R(x,equ) - equ.R₀)^2 + Z(x,equ)^2 -r(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = sqrt(r²(x, equ)) -X(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = R(x,equ) * cos(ϕ(x,equ)) -Y(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = R(x,equ) * sin(ϕ(x,equ)) -θ(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = atan(Z(x,equ), R(x,equ) - equ.R₀) + function AxisymmetricTokamakCylindricalEquilibrium{T}(R₀::T, B₀::T, q₀::T) where T <: Number + new("AxisymmetricTokamakCylindricalEquilibrium", R₀, B₀, q₀) + end + end -J(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = R(x,equ) + AxisymmetricTokamakCylindricalEquilibrium(R₀::T=DEFAULT_R₀, B₀::T=DEFAULT_B₀, q₀::T=DEFAULT_q₀) where T <: Number = AxisymmetricTokamakCylindricalEquilibrium{T}(R₀, B₀, q₀) -A₁(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = + equ.B₀ * equ.R₀ * Z(x,equ) / R(x,equ) / 2 -A₂(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = - equ.B₀ * equ.R₀ * log(R(x,equ) / equ.R₀) / 2 -A₃(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = - equ.B₀ * r²(x,equ) / equ.q₀ / 2 + function Base.show(io::IO, equ::AxisymmetricTokamakCylindricalEquilibrium) + print(io, "Axisymmetric Tokamak Equilibrium in (R,Z,ϕ) Coordinates with\n") + print(io, " R₀ = ", equ.R₀, "\n") + print(io, " B₀ = ", equ.B₀, "\n") + print(io, " q₀ = ", equ.q₀) + end -x¹(ξ::AbstractVector, equ::AxisymmetricTokamakCylindrical) = X(ξ,equ) -x²(ξ::AbstractVector, equ::AxisymmetricTokamakCylindrical) = Y(ξ,equ) -x³(ξ::AbstractVector, equ::AxisymmetricTokamakCylindrical) = Z(ξ,equ) -ξ¹(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = sqrt(x[1]^2 + x[2]^2) -ξ²(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = x[3] -ξ³(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = atan(x[2], x[1]) + R(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = x[1] + Z(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = x[2] + ϕ(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = x[3] + r²(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = (R(x,equ) - equ.R₀)^2 + Z(x,equ)^2 + r(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = sqrt(r²(x, equ)) + X(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = R(x,equ) * cos(ϕ(x,equ)) + Y(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = R(x,equ) * sin(ϕ(x,equ)) + θ(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = atan(Z(x,equ), R(x,equ) - equ.R₀) -g₁₁(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = one(eltype(x)) -g₂₂(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = one(eltype(x)) -g₃₃(x::AbstractVector, equ::AxisymmetricTokamakCylindrical) = R(x, equ)^2 + ElectromagneticFields.J(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = R(x,equ) + ElectromagneticFields.A₁(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = + equ.B₀ * equ.R₀ * Z(x,equ) / R(x,equ) / 2 + ElectromagneticFields.A₂(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = - equ.B₀ * equ.R₀ * log(R(x,equ) / equ.R₀) / 2 + ElectromagneticFields.A₃(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = - equ.B₀ * r²(x,equ) / equ.q₀ / 2 -get_functions(::AxisymmetricTokamakCylindrical) = (X=X, Y=Y, Z=Z, R=R, r=r, θ=θ, ϕ=ϕ, r²=r²) + ElectromagneticFields.x¹(ξ::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = X(ξ,equ) + ElectromagneticFields.x²(ξ::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = Y(ξ,equ) + ElectromagneticFields.x³(ξ::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = Z(ξ,equ) -function periodicity(x::AbstractVector, ::AxisymmetricTokamakCylindrical) - p = zero(x) - p[3] = 2π - return p -end + ElectromagneticFields.ξ¹(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = sqrt(x[1]^2 + x[2]^2) + ElectromagneticFields.ξ²(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = x[3] + ElectromagneticFields.ξ³(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = atan(x[2], x[1]) + ElectromagneticFields.g₁₁(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = one(eltype(x)) + ElectromagneticFields.g₂₂(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = one(eltype(x)) + ElectromagneticFields.g₃₃(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = R(x, equ)^2 + + ElectromagneticFields.get_functions(::AxisymmetricTokamakCylindricalEquilibrium) = (X=X, Y=Y, Z=Z, R=R, r=r, θ=θ, ϕ=ϕ, r²=r²) + + function ElectromagneticFields.periodicity(x::AbstractVector, ::AxisymmetricTokamakCylindricalEquilibrium) + p = zero(x) + p[3] = 2π + return p + end + + macro axisymmetric_tokamak_equilibrium_cylindrical(R₀, B₀, q₀) + generate_equilibrium_code(AxisymmetricTokamakCylindricalEquilibrium(R₀, B₀, q₀); output=false) + end + + function init(R₀=DEFAULT_R₀, B₀=DEFAULT_B₀, q₀=DEFAULT_q₀; perturbation=ZeroPerturbation()) + equilibrium = AxisymmetricTokamakCylindricalEquilibrium(R₀, B₀, q₀) + load_equilibrium(equilibrium, perturbation; target_module=AxisymmetricTokamakCylindrical) + return equilibrium + end -macro axisymmetric_tokamak_equilibrium_cylindrical(R₀, B₀, q₀) - generate_equilibrium_code(AxisymmetricTokamakCylindrical(R₀, B₀, q₀); output=false) end diff --git a/src/analytic/axisymmetric_tokamak_toroidal.jl b/src/analytic/axisymmetric_tokamak_toroidal.jl index e21fab0..2037990 100644 --- a/src/analytic/axisymmetric_tokamak_toroidal.jl +++ b/src/analytic/axisymmetric_tokamak_toroidal.jl @@ -1,5 +1,3 @@ -using SymEngine - @doc raw""" Axisymmetric tokamak equilibrium in (r,θ,ϕ) coordinates with covariant components of the vector potential given by @@ -17,65 +15,82 @@ Parameters: * `B₀`: B-field at magnetic axis * `q₀`: safety factor at magnetic axis """ -struct AxisymmetricTokamakToroidal{T <: Number} <: AnalyticEquilibrium - name::String - R₀::T - B₀::T - q₀::T - - function AxisymmetricTokamakToroidal{T}(R₀::T, B₀::T, q₀::T) where T <: Number - new("AxisymmetricTokamakEquilibriumToroidal", R₀, B₀, q₀) - end -end +module AxisymmetricTokamakToroidal -AxisymmetricTokamakToroidal(R₀::T=1.0, B₀::T=1.0, q₀::T=2.0) where T <: Number = AxisymmetricTokamakToroidal{T}(R₀, B₀, q₀) + import ..ElectromagneticFields + import ..ElectromagneticFields: AnalyticEquilibrium, ZeroPerturbation + import ..ElectromagneticFields: load_equilibrium, generate_equilibrium_code + export AxisymmetricTokamakToroidalEquilibrium -function Base.show(io::IO, equ::AxisymmetricTokamakToroidal) - print(io, "Axisymmetric Tokamak Equilibrium in Toroidal Coordinates with\n") - print(io, " R₀ = ", equ.R₀, "\n") - print(io, " B₀ = ", equ.B₀, "\n") - print(io, " q₀ = ", equ.q₀) -end + const DEFAULT_R₀ = 1.0 + const DEFAULT_B₀ = 1.0 + const DEFAULT_q₀ = 2.0 + struct AxisymmetricTokamakToroidalEquilibrium{T <: Number} <: AnalyticEquilibrium + name::String + R₀::T + B₀::T + q₀::T -r(x::AbstractVector, equ::AxisymmetricTokamakToroidal) = x[1] -θ(x::AbstractVector, equ::AxisymmetricTokamakToroidal) = x[2] -ϕ(x::AbstractVector, equ::AxisymmetricTokamakToroidal) = x[3] -R(x::AbstractVector, equ::AxisymmetricTokamakToroidal) = equ.R₀ + r(x,equ) * cos(θ(x,equ)) -X(x::AbstractVector, equ::AxisymmetricTokamakToroidal) = R(x,equ) * cos(ϕ(x,equ)) -Y(x::AbstractVector, equ::AxisymmetricTokamakToroidal) = R(x,equ) * sin(ϕ(x,equ)) -Z(x::AbstractVector, equ::AxisymmetricTokamakToroidal) = r(x,equ) * sin(θ(x,equ)) + function AxisymmetricTokamakToroidalEquilibrium{T}(R₀::T, B₀::T, q₀::T) where T <: Number + new("AxisymmetricTokamakEquilibriumToroidal", R₀, B₀, q₀) + end + end -J(x::AbstractVector, equ::AxisymmetricTokamakToroidal) = r(x,equ) * R(x,equ) + AxisymmetricTokamakToroidalEquilibrium(R₀::T=DEFAULT_R₀, B₀::T=DEFAULT_B₀, q₀::T=DEFAULT_q₀) where T <: Number = AxisymmetricTokamakToroidalEquilibrium{T}(R₀, B₀, q₀) -A₁(x::AbstractVector, equ::AxisymmetricTokamakToroidal) = zero(eltype(x)) -A₂(x::AbstractVector, equ::AxisymmetricTokamakToroidal) = + equ.B₀ * equ.R₀ / cos(θ(x,equ))^2 * ( r(x,equ) * cos(θ(x,equ)) - equ.R₀ * log(R(x,equ) / equ.R₀) ) -A₃(x::AbstractVector, equ::AxisymmetricTokamakToroidal) = - equ.B₀ * r(x,equ)^2 / equ.q₀ / 2 + function Base.show(io::IO, equ::AxisymmetricTokamakToroidalEquilibrium) + print(io, "Axisymmetric Tokamak Equilibrium in Toroidal Coordinates with\n") + print(io, " R₀ = ", equ.R₀, "\n") + print(io, " B₀ = ", equ.B₀, "\n") + print(io, " q₀ = ", equ.q₀) + end -x¹(ξ::AbstractVector, equ::AxisymmetricTokamakToroidal) = X(ξ,equ) -x²(ξ::AbstractVector, equ::AxisymmetricTokamakToroidal) = Y(ξ,equ) -x³(ξ::AbstractVector, equ::AxisymmetricTokamakToroidal) = Z(ξ,equ) -ξ¹(x::AbstractVector, equ::AxisymmetricTokamakToroidal) = sqrt((sqrt(x[1]^2 + x[2]^2)-equ.R₀)^2 + x[3]^2) -ξ²(x::AbstractVector, equ::AxisymmetricTokamakToroidal) = atan(x[3], sqrt(x[1]^2 + x[2]^2)-equ.R₀) -ξ³(x::AbstractVector, equ::AxisymmetricTokamakToroidal) = atan(x[2], x[1]) + 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)) -g₁₁(x::AbstractVector, equ::AxisymmetricTokamakToroidal) = one(eltype(x)) -g₂₂(x::AbstractVector, equ::AxisymmetricTokamakToroidal) = r(x, equ)^2 -g₃₃(x::AbstractVector, equ::AxisymmetricTokamakToroidal) = R(x, equ)^2 + ElectromagneticFields.J(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = r(x,equ) * R(x,equ) + ElectromagneticFields.A₁(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = zero(eltype(x)) + ElectromagneticFields.A₂(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = + 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::AxisymmetricTokamakToroidalEquilibrium) = - equ.B₀ * r(x,equ)^2 / equ.q₀ / 2 -get_functions(::AxisymmetricTokamakToroidal) = (X=X, Y=Y, Z=Z, R=R, r=r, θ=θ, ϕ=ϕ) + 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::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]) -function periodicity(x::AbstractVector, equ::AxisymmetricTokamakToroidal) - p = zero(x) - p[2] = 2π - p[3] = 2π - return p -end + 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.get_functions(::AxisymmetricTokamakToroidalEquilibrium) = (X=X, Y=Y, Z=Z, R=R, r=r, θ=θ, ϕ=ϕ) + + function ElectromagneticFields.periodicity(x::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) + p = zero(x) + p[2] = 2π + p[3] = 2π + return p + end + + macro axisymmetric_tokamak_equilibrium_toroidal(R₀, B₀, q₀) + generate_equilibrium_code(AxisymmetricTokamakToroidalEquilibrium(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) + return equilibrium + end -macro axisymmetric_tokamak_equilibrium_toroidal(R₀, B₀, q₀) - generate_equilibrium_code(AxisymmetricTokamakToroidal(R₀, B₀, q₀); output=false) end diff --git a/src/analytic/cartesian_equilibrium.jl b/src/analytic/cartesian_equilibrium.jl index 67f60dc..e80f65d 100644 --- a/src/analytic/cartesian_equilibrium.jl +++ b/src/analytic/cartesian_equilibrium.jl @@ -1,24 +1,31 @@ +module AnalyticCartesianField -abstract type CartesianEquilibrium <: AnalyticEquilibrium end -abstract type CartesianPerturbation <: AnalyticPerturbation end + import ..ElectromagneticFields + import ..ElectromagneticFields: AnalyticEquilibrium, AnalyticPerturbation -CartesianField = Union{CartesianEquilibrium, CartesianPerturbation} + export CartesianField, CartesianEquilibrium, CartesianPerturbation + abstract type CartesianEquilibrium <: AnalyticEquilibrium end + abstract type CartesianPerturbation <: AnalyticPerturbation end -X(x::AbstractVector, ::CartesianField) = x[1] -Y(x::AbstractVector, ::CartesianField) = x[2] -Z(x::AbstractVector, ::CartesianField) = x[3] + CartesianField = Union{CartesianEquilibrium, CartesianPerturbation} -J(x::AbstractVector, ::CartesianField) = one(eltype(x)) + X(x::AbstractVector, ::CartesianField) = x[1] + Y(x::AbstractVector, ::CartesianField) = x[2] + Z(x::AbstractVector, ::CartesianField) = x[3] -x¹(ξ::AbstractVector, ::CartesianField) = ξ[1] -x²(ξ::AbstractVector, ::CartesianField) = ξ[2] -x³(ξ::AbstractVector, ::CartesianField) = ξ[3] + ElectromagneticFields.J(x::AbstractVector, ::CartesianField) = one(eltype(x)) -ξ¹(x::AbstractVector, ::CartesianField) = x[1] -ξ²(x::AbstractVector, ::CartesianField) = x[2] -ξ³(x::AbstractVector, ::CartesianField) = x[3] + ElectromagneticFields.x¹(ξ::AbstractVector, ::CartesianField) = ξ[1] + ElectromagneticFields.x²(ξ::AbstractVector, ::CartesianField) = ξ[2] + ElectromagneticFields.x³(ξ::AbstractVector, ::CartesianField) = ξ[3] -g₁₁(x::AbstractVector, ::CartesianField) = one(eltype(x)) -g₂₂(x::AbstractVector, ::CartesianField) = one(eltype(x)) -g₃₃(x::AbstractVector, ::CartesianField) = one(eltype(x)) + ElectromagneticFields.ξ¹(x::AbstractVector, ::CartesianField) = x[1] + ElectromagneticFields.ξ²(x::AbstractVector, ::CartesianField) = x[2] + ElectromagneticFields.ξ³(x::AbstractVector, ::CartesianField) = x[3] + + ElectromagneticFields.g₁₁(x::AbstractVector, ::CartesianField) = one(eltype(x)) + ElectromagneticFields.g₂₂(x::AbstractVector, ::CartesianField) = one(eltype(x)) + ElectromagneticFields.g₃₃(x::AbstractVector, ::CartesianField) = one(eltype(x)) + +end diff --git a/src/analytic/ezcosz.jl b/src/analytic/ezcosz.jl index 70a1cda..f2cec18 100644 --- a/src/analytic/ezcosz.jl +++ b/src/analytic/ezcosz.jl @@ -10,28 +10,42 @@ E(x,y,z) = E_0 \, \begin{pmatrix} Parameters: `E₀` """ -struct EzCosZ{T <: Number} <: CartesianPerturbation - name::String - E₀::T - EzCosZ{T}(E₀::T) where T <: Number = new("EzCosZ", E₀) -end +module EzCosZ -EzCosZ(E₀::T) where T <: Number = EzCosZ{T}(E₀) + import ..ElectromagneticFields + import ..ElectromagneticFields: CartesianPerturbation + import ..ElectromagneticFields: load_equilibrium, generate_equilibrium_code + import ..AnalyticCartesianField: X, Y, Z + export EzCosZPerturbation -function Base.show(io::IO, equ::EzCosZ) - print(io, "Simple perturbation in electric field") -end + const DEFAULT_E₀ = 1.0 + struct EzCosZPerturbation{T <: Number} <: CartesianPerturbation + name::String + E₀::T + EzCosZPerturbation{T}(E₀::T) where T <: Number = new("EzCosZ", E₀) + end -function φ(x::AbstractArray{T,1}, equ::EzCosZ) where {T <: Number} - equ.E₀ / (2π) * sin(2π * Z(x,equ)) -end + EzCosZPerturbation(E₀::T=DEFAULT_E₀) where T <: Number = EzCosZPerturbation{T}(E₀) + + + function Base.show(io::IO, equ::EzCosZPerturbation) + print(io, "Simple perturbation in electric field") + end + + + ElectromagneticFields.φ(x::AbstractVector, equ::EzCosZPerturbation) = equ.E₀ / (2π) * sin(2π * Z(x,equ)) + + ElectromagneticFields.get_functions(::EzCosZPerturbation) = (X=X, Y=Y, Z=Z) -get_functions(::EzCosZ) = (X=X, Y=Y, Z=Z) + macro ezcosz_perturbation(E₀=1.) + generate_equilibrium_code(EzCosZPerturbation(E₀); output=false) + end + function init(E₀=DEFAULT_E₀) + EzCosZPerturbation(E₀) + end -macro ezcosz_perturbation(E₀=1.) - generate_equilibrium_code(EzCosZ(E₀); output=false) end diff --git a/src/analytic/solovev.jl b/src/analytic/solovev.jl index dea1570..f759f19 100644 --- a/src/analytic/solovev.jl +++ b/src/analytic/solovev.jl @@ -1,7 +1,3 @@ - -# using SymEngine: N, symbols, diff, expand, subs -using SymPy: N, Sym, diff, expand, solve, subs - """ Axisymmetric Solov'ev equilibra in (R/R₀,Z/R₀,phi) coordinates. Based on Cerfon & Freidberg, Physics of Plasmas 17, 032502, 2010. @@ -14,100 +10,100 @@ Parameters: * `δ`: triangularity * `a`: free constant, determined to match a given beta value """ -struct Solovev{T <: Number} <: AbstractSolovevEquilibrium - name::String - R₀::T - B₀::T - ϵ::T - κ::T - δ::T - a::T - c::Vector{T} - - function Solovev{T}(R₀::T, B₀::T, ϵ::T, κ::T, δ::T, a::T, c::Vector{T}) where T <: Number - new("SolovevEquilibrium", R₀, B₀, ϵ, κ, δ, a, c) +module Solovev + + using SymPy: N, Sym, diff, expand, solve, subs + + import ..ElectromagneticFields + import ..ElectromagneticFields: ZeroPerturbation + import ..ElectromagneticFields: load_equilibrium, generate_equilibrium_code + import ..SolovevAbstract: AbstractSolovevEquilibrium, X, Y, Z, R, r, θ, ϕ, r² + + export SolovevEquilibrium + + include("solovev_psi.jl") + + struct SolovevEquilibrium{T <: Number} <: AbstractSolovevEquilibrium + name::String + R₀::T + B₀::T + ϵ::T + κ::T + δ::T + a::T + c::Vector{T} + + function SolovevEquilibrium{T}(R₀::T, B₀::T, ϵ::T, κ::T, δ::T, a::T, c::Vector{T}) where T <: Number + new("SolovevEquilibrium", R₀, B₀, ϵ, κ, δ, a, c) + end end -end -function Solovev(R₀::T, B₀::T, ϵ::T, κ::T, δ::T, a::T) where T <: Number - - n = 7 - - x = [Sym("x" * string(i)) for i in 1:3] - c = [Sym("c" * string(i)) for i in 1:n] - - ψ = (ψ₀(x,a) + c[1] * ψ₁(x) - + c[2] * ψ₂(x) - + c[3] * ψ₃(x) - + c[4] * ψ₄(x) - + c[5] * ψ₅(x) - + c[6] * ψ₆(x) - + c[7] * ψ₇(x) ) - - eqs = [ - expand(subs(subs(ψ, x[1], 1+ϵ), x[2], 0)), - expand(subs(subs(ψ, x[1], 1-ϵ), x[2], 0)), - expand(subs(subs(ψ, x[1], 1-δ*ϵ), x[2], κ*ϵ)), - expand(subs(subs(diff(ψ, x[1]), x[1], 1-δ*ϵ), x[2], κ*ϵ)), - expand(subs(subs(diff(ψ, x[2], 2), x[1], 1+ϵ), x[2], 0) - (1 + asin(δ))^2 / (ϵ * κ^2) * subs(subs(diff(ψ, x[1]), x[1], 1+ϵ), x[2], 0)), - expand(subs(subs(diff(ψ, x[2], 2), x[1], 1-ϵ), x[2], 0) + (1 - asin(δ))^2 / (ϵ * κ^2) * subs(subs(diff(ψ, x[1]), x[1], 1-ϵ), x[2], 0)), - expand(subs(subs(diff(ψ, x[1], 2), x[1], 1-δ*ϵ), x[2], κ*ϵ) - κ / (ϵ * (1 - δ^2)) * subs(subs(diff(ψ, x[2]), x[1], 1-δ*ϵ), x[2], κ*ϵ)) - ] - - # x = [symbols("x" * string(i)) for i in 1:3] - # c = [symbols("c" * string(i)) for i in 1:n] - # - # ψ = (ψ₀(x,a) + c[1] * ψ₁(x) - # + c[2] * ψ₂(x) - # + c[3] * ψ₃(x) - # + c[4] * ψ₄(x) - # + c[5] * ψ₅(x) - # + c[6] * ψ₆(x) - # + c[7] * ψ₇(x) ) - # - # eqs = [ - # expand(subs(subs(ψ, x[1]=>1+ϵ), x[2]=>0)), - # expand(subs(subs(ψ, x[1]=>1-ϵ), x[2]=>0)), - # expand(subs(subs(ψ, x[1]=>1-δ*ϵ), x[2]=>κ*ϵ)), - # expand(subs(subs(diff(ψ, x[1]), x[1]=>1-δ*ϵ), x[2]=>κ*ϵ)), - # expand(subs(subs(diff(ψ, x[2], 2), x[1]=>1+ϵ), x[2]=>0) - (1 + asin(δ))^2 / (ϵ * κ^2) * subs(subs(diff(ψ, x[1]), x[1]=>1+ϵ), x[2]=>0)), - # expand(subs(subs(diff(ψ, x[2], 2), x[1]=>1-ϵ), x[2]=>0) + (1 - asin(δ))^2 / (ϵ * κ^2) * subs(subs(diff(ψ, x[1]), x[1]=>1-ϵ), x[2]=>0)), - # expand(subs(subs(diff(ψ, x[1], 2), x[1]=>1-δ*ϵ), x[2]=>κ*ϵ) - κ / (ϵ * (1 - δ^2)) * subs(subs(diff(ψ, x[2]), x[1]=>1-δ*ϵ), x[2]=>κ*ϵ)) - # ] - - csym = solve(eqs, c) - cnum = [N(csym[c[i]]) for i in 1:n] - - Solovev{T}(R₀, B₀, ϵ, κ, δ, a, cnum) -end + function SolovevEquilibrium(R₀::T, B₀::T, ϵ::T, κ::T, δ::T, a::T) where T <: Number + n = 7 -SolovevITER() = Solovev(6.2, 5.3, 0.32, 1.7, 0.33, -0.155) -SolovevNSTX() = Solovev(0.85, 0.3, 0.78, 2.0, 0.35, 1.0) + x = [Sym("x" * string(i)) for i in 1:3] + c = [Sym("c" * string(i)) for i in 1:n] + ψ = (ψ₀(x,a) + c[1] * ψ₁(x) + + c[2] * ψ₂(x) + + c[3] * ψ₃(x) + + c[4] * ψ₄(x) + + c[5] * ψ₅(x) + + c[6] * ψ₆(x) + + c[7] * ψ₇(x) ) -function Base.show(io::IO, equ::Solovev) - print(io, "Solovev Equilibrium with\n") - print(io, " R₀ = ", equ.R₀, "\n") - print(io, " B₀ = ", equ.B₀, "\n") - print(io, " ϵ = ", equ.ϵ, "\n") - print(io, " κ = ", equ.κ, "\n") - print(io, " δ = ", equ.δ, "\n") - print(io, " a = ", equ.a) -end + eqs = [ + expand(subs(subs(ψ, x[1], 1+ϵ), x[2], 0)), + expand(subs(subs(ψ, x[1], 1-ϵ), x[2], 0)), + expand(subs(subs(ψ, x[1], 1-δ*ϵ), x[2], κ*ϵ)), + expand(subs(subs(diff(ψ, x[1]), x[1], 1-δ*ϵ), x[2], κ*ϵ)), + expand(subs(subs(diff(ψ, x[2], 2), x[1], 1+ϵ), x[2], 0) - (1 + asin(δ))^2 / (ϵ * κ^2) * subs(subs(diff(ψ, x[1]), x[1], 1+ϵ), x[2], 0)), + expand(subs(subs(diff(ψ, x[2], 2), x[1], 1-ϵ), x[2], 0) + (1 - asin(δ))^2 / (ϵ * κ^2) * subs(subs(diff(ψ, x[1]), x[1], 1-ϵ), x[2], 0)), + expand(subs(subs(diff(ψ, x[1], 2), x[1], 1-δ*ϵ), x[2], κ*ϵ) - κ / (ϵ * (1 - δ^2)) * subs(subs(diff(ψ, x[2]), x[1], 1-δ*ϵ), x[2], κ*ϵ)) + ] + csym = solve(eqs, c) + cnum = [N(csym[c[i]]) for i in 1:n] -function A₃(x::AbstractArray{T,1}, equ::Solovev) where {T <: Number} - (ψ₀(x, equ.a) + equ.c[1] * ψ₁(x) - + equ.c[2] * ψ₂(x) - + equ.c[3] * ψ₃(x) - + equ.c[4] * ψ₄(x) - + equ.c[5] * ψ₅(x) - + equ.c[6] * ψ₆(x) - + equ.c[7] * ψ₇(x) ) -end + SolovevEquilibrium{T}(R₀, B₀, ϵ, κ, δ, a, cnum) + end -macro solovev_equilibrium(R₀, B₀, ϵ, κ, δ, a) - generate_equilibrium_code(Solovev(R₀, B₀, ϵ, κ, δ, a); output=false) + SolovevEquilibriumITER() = SolovevEquilibrium(6.2, 5.3, 0.32, 1.7, 0.33, -0.155) + SolovevEquilibriumNSTX() = SolovevEquilibrium(0.85, 0.3, 0.78, 2.0, 0.35, 1.0) + + + function Base.show(io::IO, equ::SolovevEquilibrium) + print(io, "SolovevEquilibrium Equilibrium with\n") + print(io, " R₀ = ", equ.R₀, "\n") + print(io, " B₀ = ", equ.B₀, "\n") + print(io, " ϵ = ", equ.ϵ, "\n") + print(io, " κ = ", equ.κ, "\n") + print(io, " δ = ", equ.δ, "\n") + print(io, " a = ", equ.a) + end + + + function ElectromagneticFields.A₃(x::AbstractArray{T,1}, equ::SolovevEquilibrium) where {T <: Number} + ψ₀(x, equ.a) + equ.c[1] * ψ₁(x) + + equ.c[2] * ψ₂(x) + + equ.c[3] * ψ₃(x) + + equ.c[4] * ψ₄(x) + + equ.c[5] * ψ₅(x) + + equ.c[6] * ψ₆(x) + + equ.c[7] * ψ₇(x) + end + + + macro solovev_equilibrium(R₀, B₀, ϵ, κ, δ, a) + generate_equilibrium_code(SolovevEquilibrium(R₀, B₀, ϵ, κ, δ, a); output=false) + end + + function init(R₀, B₀, ϵ, κ, δ, a; perturbation=ZeroPerturbation()) + equilibrium = SolovevEquilibrium(R₀, B₀, ϵ, κ, δ, a) + load_equilibrium(equilibrium, perturbation; target_module=Solovev) + return equilibrium + end + end diff --git a/src/analytic/solovev_abstract.jl b/src/analytic/solovev_abstract.jl new file mode 100644 index 0000000..5c4eb0c --- /dev/null +++ b/src/analytic/solovev_abstract.jl @@ -0,0 +1,47 @@ +module SolovevAbstract + + import ..ElectromagneticFields + import ..ElectromagneticFields: AnalyticEquilibrium, AnalyticPerturbation + + export AbstractSolovevEquilibrium + + abstract type AbstractSolovevEquilibrium <: AnalyticEquilibrium end + + + Z(x::AbstractVector, equ::AbstractSolovevEquilibrium) = x[2] * equ.R₀ + R(x::AbstractVector, equ::AbstractSolovevEquilibrium) = x[1] * equ.R₀ + ϕ(x::AbstractVector, equ::AbstractSolovevEquilibrium) = x[3] + + X(x::AbstractVector, equ::AbstractSolovevEquilibrium) = R(x,equ) * cos(ϕ(x,equ)) + Y(x::AbstractVector, equ::AbstractSolovevEquilibrium) = R(x,equ) * sin(ϕ(x,equ)) + θ(x::AbstractVector, equ::AbstractSolovevEquilibrium) = atan(Z(x,equ), R(x,equ) - equ.R₀) + + r²(x::AbstractVector, equ::AbstractSolovevEquilibrium) = (R(x,equ) - equ.R₀)^2 + Z(x,equ)^2 + r(x::AbstractVector, equ::AbstractSolovevEquilibrium) = sqrt(r²(x, equ)) + + ElectromagneticFields.J(x::AbstractVector, equ::AbstractSolovevEquilibrium) = R(x,equ) + + ElectromagneticFields.A₁(x::AbstractVector, equ::AbstractSolovevEquilibrium) = + equ.B₀ * equ.R₀ * x[2] / x[1] / 2 + ElectromagneticFields.A₂(x::AbstractVector, equ::AbstractSolovevEquilibrium) = - equ.B₀ * equ.R₀ * log(x[1]) / 2 + + ElectromagneticFields.x¹(ξ::AbstractVector, equ::AbstractSolovevEquilibrium) = X(ξ,equ) + ElectromagneticFields.x²(ξ::AbstractVector, equ::AbstractSolovevEquilibrium) = Y(ξ,equ) + ElectromagneticFields.x³(ξ::AbstractVector, equ::AbstractSolovevEquilibrium) = Z(ξ,equ) + + ElectromagneticFields.ξ¹(x::AbstractVector, equ::AbstractSolovevEquilibrium) = sqrt(x[1]^2 + x[2]^2) / equ.R₀ + ElectromagneticFields.ξ²(x::AbstractVector, equ::AbstractSolovevEquilibrium) = x[3] / equ.R₀ + ElectromagneticFields.ξ³(x::AbstractVector, equ::AbstractSolovevEquilibrium) = atan(x[2], x[1]) + + ElectromagneticFields.g₁₁(x::AbstractVector, equ::AbstractSolovevEquilibrium) = one(eltype(x)) + ElectromagneticFields.g₂₂(x::AbstractVector, equ::AbstractSolovevEquilibrium) = one(eltype(x)) + ElectromagneticFields.g₃₃(x::AbstractVector, equ::AbstractSolovevEquilibrium) = R(x, equ)^2 + + ElectromagneticFields.get_functions(::AbstractSolovevEquilibrium) = (X=X, Y=Y, Z=Z, R=R, r=r, θ=θ, ϕ=ϕ, r²=r²) + + function ElectromagneticFields.periodicity(x::AbstractArray{T,1}, equ::AbstractSolovevEquilibrium) where {T <: Number} + p = zero(x) + p[3] = 2π + return p + end + +end diff --git a/src/analytic/solovev_common.jl b/src/analytic/solovev_common.jl deleted file mode 100644 index 3d4167e..0000000 --- a/src/analytic/solovev_common.jl +++ /dev/null @@ -1,39 +0,0 @@ - -abstract type AbstractSolovevEquilibrium <: AnalyticEquilibrium end - - -Z(x::AbstractVector, equ::AbstractSolovevEquilibrium) = x[2] * equ.R₀ -R(x::AbstractVector, equ::AbstractSolovevEquilibrium) = x[1] * equ.R₀ -ϕ(x::AbstractVector, equ::AbstractSolovevEquilibrium) = x[3] - -X(x::AbstractVector, equ::AbstractSolovevEquilibrium) = R(x,equ) * cos(ϕ(x,equ)) -Y(x::AbstractVector, equ::AbstractSolovevEquilibrium) = R(x,equ) * sin(ϕ(x,equ)) -θ(x::AbstractVector, equ::AbstractSolovevEquilibrium) = atan(Z(x,equ), R(x,equ) - equ.R₀) - -r²(x::AbstractVector, equ::AbstractSolovevEquilibrium) = (R(x,equ) - equ.R₀)^2 + Z(x,equ)^2 -r(x::AbstractVector, equ::AbstractSolovevEquilibrium) = sqrt(r²(x, equ)) - -J(x::AbstractVector, equ::AbstractSolovevEquilibrium) = R(x,equ) - -A₁(x::AbstractVector, equ::AbstractSolovevEquilibrium) = + equ.B₀ * equ.R₀ * x[2] / x[1] / 2 -A₂(x::AbstractVector, equ::AbstractSolovevEquilibrium) = - equ.B₀ * equ.R₀ * log(x[1]) / 2 - -x¹(ξ::AbstractVector, equ::AbstractSolovevEquilibrium) = X(ξ,equ) -x²(ξ::AbstractVector, equ::AbstractSolovevEquilibrium) = Y(ξ,equ) -x³(ξ::AbstractVector, equ::AbstractSolovevEquilibrium) = Z(ξ,equ) - -ξ¹(x::AbstractVector, equ::AbstractSolovevEquilibrium) = sqrt(x[1]^2 + x[2]^2) / equ.R₀ -ξ²(x::AbstractVector, equ::AbstractSolovevEquilibrium) = x[3] / equ.R₀ -ξ³(x::AbstractVector, equ::AbstractSolovevEquilibrium) = atan(x[2], x[1]) - -g₁₁(x::AbstractVector, equ::AbstractSolovevEquilibrium) = one(eltype(x)) -g₂₂(x::AbstractVector, equ::AbstractSolovevEquilibrium) = one(eltype(x)) -g₃₃(x::AbstractVector, equ::AbstractSolovevEquilibrium) = R(x, equ)^2 - -get_functions(::AbstractSolovevEquilibrium) = (X=X, Y=Y, Z=Z, R=R, r=r, θ=θ, ϕ=ϕ, r²=r²) - -function periodicity(x::AbstractArray{T,1}, equ::AbstractSolovevEquilibrium) where {T <: Number} - p = zero(x) - p[3] = 2π - return p -end diff --git a/src/analytic/solovev_quadratic.jl b/src/analytic/solovev_quadratic.jl index 82c6cb8..930cca9 100644 --- a/src/analytic/solovev_quadratic.jl +++ b/src/analytic/solovev_quadratic.jl @@ -7,37 +7,54 @@ Parameters: * `B₀`: B-field at magnetic axis * `a`, b`: free constants """ -struct SolovevQuadratic{T <: Number} <: AbstractSolovevEquilibrium - name::String - R₀::T - B₀::T - a::T - b::T - - function SolovevQuadratic{T}(R₀::T, B₀::T, a::T, b::T) where T <: Number - new("QuadraticSolovevEquilibrium", R₀, B₀, a, b) +module SolovevQuadratic + + import ..ElectromagneticFields + import ..ElectromagneticFields: ZeroPerturbation + import ..ElectromagneticFields: load_equilibrium, generate_equilibrium_code + import ..SolovevAbstract: AbstractSolovevEquilibrium, X, Y, Z, R, r, θ, ϕ, r² + + export SolovevQuadraticEquilibrium + + struct SolovevQuadraticEquilibrium{T <: Number} <: AbstractSolovevEquilibrium + name::String + R₀::T + B₀::T + a::T + b::T + + function SolovevQuadraticEquilibrium{T}(R₀::T, B₀::T, a::T, b::T) where T <: Number + new("QuadraticSolovevEquilibrium", R₀, B₀, a, b) + end end -end -function SolovevQuadratic(R₀::T, B₀::T, a::T, b::T) where T <: Number - SolovevQuadratic{T}(R₀, B₀, a, b) -end + function SolovevQuadraticEquilibrium(R₀::T, B₀::T, a::T, b::T) where T <: Number + SolovevQuadraticEquilibrium{T}(R₀, B₀, a, b) + end -function Base.show(io::IO, equ::SolovevQuadratic) - print(io, "Quadratic Solovev Equilibrium with\n") - print(io, " R₀ = ", equ.R₀, "\n") - print(io, " B₀ = ", equ.B₀, "\n") - print(io, " a = ", equ.a, "\n") - print(io, " b = ", equ.b) -end + function Base.show(io::IO, equ::SolovevQuadraticEquilibrium) + print(io, "Quadratic Solovev Equilibrium with\n") + print(io, " R₀ = ", equ.R₀, "\n") + print(io, " B₀ = ", equ.B₀, "\n") + print(io, " a = ", equ.a, "\n") + print(io, " b = ", equ.b) + end + + ElectromagneticFields.A₁(x::AbstractVector, equ::SolovevQuadraticEquilibrium) = zero(eltype(x)) + ElectromagneticFields.A₂(x::AbstractVector, equ::SolovevQuadraticEquilibrium) = zero(eltype(x)) + ElectromagneticFields.A₃(x::AbstractVector, equ::SolovevQuadraticEquilibrium) = - equ.a * R(x,equ)^4 / 8 - equ.b * Z(x,equ)^2 / 2 -A₁(x::AbstractVector, equ::SolovevQuadratic) = zero(eltype(x)) -A₂(x::AbstractVector, equ::SolovevQuadratic) = zero(eltype(x)) -A₃(x::AbstractVector, equ::SolovevQuadratic) = - equ.a * R(x,equ)^4 / 8 - equ.b * Z(x,equ)^2 / 2 + macro solovev_equilibrium_quadratic(R₀, B₀, a, b) + generate_equilibrium_code(SolovevQuadraticEquilibrium(R₀, B₀, a, b); output=false) + end + + function init(R₀, B₀, a, b; perturbation=ZeroPerturbation()) + equilibrium = SolovevQuadraticEquilibrium(R₀, B₀, a, b) + load_equilibrium(equilibrium, perturbation; target_module=SolovevQuadratic) + return equilibrium + end -macro solovev_equilibrium_quadratic(R₀, B₀, a, b) - generate_equilibrium_code(SolovevQuadratic(R₀, B₀, a, b); output=false) end diff --git a/src/analytic/solovev_xpoint.jl b/src/analytic/solovev_xpoint.jl index 33c59fd..92c3f92 100644 --- a/src/analytic/solovev_xpoint.jl +++ b/src/analytic/solovev_xpoint.jl @@ -1,7 +1,3 @@ - -# using SymEngine: N, symbols, diff, expand, subs -using SymPy: N, Sym, diff, expand, solve, subs - """ Axisymmetric Solov'ev equilibra with X-point in (R/R₀,Z/R₀,phi) coordinates. Based on Cerfon & Freidberg, Physics of Plasmas 17, 032502, 2010. @@ -16,129 +12,119 @@ Parameters: * `xₛₑₚ`: x position of the X point * `yₛₑₚ`: y position of the X point """ -struct SolovevXpoint{T <: Number} <: AbstractSolovevEquilibrium - name::String - R₀::T - B₀::T - ϵ::T - κ::T - δ::T - a::T - xₛₑₚ::T - yₛₑₚ::T - c::Vector{T} - - function SolovevXpoint{T}(R₀::T, B₀::T, ϵ::T, κ::T, δ::T, a::T, xₛₑₚ::T, yₛₑₚ::T, c::Vector{T}) where T <: Number - new("SolovevXpointEquilibrium", R₀, B₀, ϵ, κ, δ, a, xₛₑₚ, yₛₑₚ, c) +module SolovevXpoint + + using SymPy: N, Sym, diff, expand, solve, subs + + import ..ElectromagneticFields + import ..ElectromagneticFields: ZeroPerturbation + import ..ElectromagneticFields: load_equilibrium, generate_equilibrium_code + import ..SolovevAbstract: AbstractSolovevEquilibrium, X, Y, Z, R, r, θ, ϕ, r² + + export SolovevXpointEquilibrium + + include("solovev_psi.jl") + + struct SolovevXpointEquilibrium{T <: Number} <: AbstractSolovevEquilibrium + name::String + R₀::T + B₀::T + ϵ::T + κ::T + δ::T + a::T + xₛₑₚ::T + yₛₑₚ::T + c::Vector{T} + + function SolovevXpointEquilibrium{T}(R₀::T, B₀::T, ϵ::T, κ::T, δ::T, a::T, xₛₑₚ::T, yₛₑₚ::T, c::Vector{T}) where T <: Number + new("SolovevXpointEquilibrium", R₀, B₀, ϵ, κ, δ, a, xₛₑₚ, yₛₑₚ, c) + end end -end -function SolovevXpoint(R₀::T, B₀::T, ϵ::T, κ::T, δ::T, a::T, xₛₑₚ::T, yₛₑₚ::T) where T <: Number - - n = 12 - - x = [Sym("x" * string(i)) for i in 1:3] - c = [Sym("c" * string(i)) for i in 1:n] - - ψ = (ψ₀(x,a) + c[1] * ψ₁(x) - + c[2] * ψ₂(x) - + c[3] * ψ₃(x) - + c[4] * ψ₄(x) - + c[5] * ψ₅(x) - + c[6] * ψ₆(x) - + c[7] * ψ₇(x) - + c[8] * ψ₈(x) - + c[9] * ψ₉(x) - + c[10] * ψ₁₀(x) - + c[11] * ψ₁₁(x) - + c[12] * ψ₁₂(x) ) - - eqs = [ - expand(subs(subs(ψ, x[1], 1+ϵ), x[2], 0)), - expand(subs(subs(ψ, x[1], 1-ϵ), x[2], 0)), - expand(subs(subs(ψ, x[1], 1-δ*ϵ), x[2], κ*ϵ)), - expand(subs(subs(ψ, x[1], xₛₑₚ), x[2], yₛₑₚ)), - expand(subs(subs(diff(ψ, x[2]), x[1], 1+ϵ), x[2], 0)), - expand(subs(subs(diff(ψ, x[2]), x[1], 1-ϵ), x[2], 0)), - expand(subs(subs(diff(ψ, x[1]), x[1], 1-δ*ϵ), x[2], κ*ϵ)), - expand(subs(subs(diff(ψ, x[1]), x[1], xₛₑₚ), x[2], yₛₑₚ)), - expand(subs(subs(diff(ψ, x[2]), x[1], xₛₑₚ), x[2], yₛₑₚ)), - expand(subs(subs(diff(ψ, x[2], 2), x[1], 1+ϵ), x[2], 0) - (1 + asin(δ))^2 / (ϵ * κ^2) * subs(subs(diff(ψ, x[1]), x[1], 1+ϵ), x[2], 0)), - expand(subs(subs(diff(ψ, x[2], 2), x[1], 1-ϵ), x[2], 0) + (1 - asin(δ))^2 / (ϵ * κ^2) * subs(subs(diff(ψ, x[1]), x[1], 1-ϵ), x[2], 0)), - expand(subs(subs(diff(ψ, x[1], 2), x[1], 1-δ*ϵ), x[2], κ*ϵ) - κ / (ϵ * (1 - δ^2)) * subs(subs(diff(ψ, x[2]), x[1], 1-δ*ϵ), x[2], κ*ϵ)) - ] - - # x = [symbols("x" * string(i)) for i in 1:3] - # c = [symbols("c" * string(i)) for i in 1:n] - # - # ψ = (ψ₀(x,a) + c[1] * ψ₁(x) - # + c[2] * ψ₂(x) - # + c[3] * ψ₃(x) - # + c[4] * ψ₄(x) - # + c[5] * ψ₅(x) - # + c[6] * ψ₆(x) - # + c[7] * ψ₇(x) - # + c[8] * ψ₈(x) - # + c[9] * ψ₉(x) - # + c[10] * ψ₁₀(x) - # + c[11] * ψ₁₁(x) - # + c[12] * ψ₁₂(x) ) - # - # eqs = [ - # expand(subs(subs(ψ, x[1]=>1+ϵ), x[2]=>0)), - # expand(subs(subs(ψ, x[1]=>1-ϵ), x[2]=>0)), - # expand(subs(subs(ψ, x[1]=>1-δ*ϵ), x[2]=>κ*ϵ)), - # expand(subs(subs(ψ, x[1]=>xₛₑₚ), x[2]=>yₛₑₚ)), - # expand(subs(subs(diff(ψ, x[2]), x[1]=>1+ϵ), x[2]=>0)), - # expand(subs(subs(diff(ψ, x[2]), x[1]=>1-ϵ), x[2]=>0)), - # expand(subs(subs(diff(ψ, x[1]), x[1]=>1-δ*ϵ), x[2]=>κ*ϵ)), - # expand(subs(subs(diff(ψ, x[1]), x[1]=>xₛₑₚ), x[2]=>yₛₑₚ)), - # expand(subs(subs(diff(ψ, x[2]), x[1]=>xₛₑₚ), x[2]=>yₛₑₚ)), - # expand(subs(subs(diff(ψ, x[2], 2), x[1]=>1+ϵ), x[2]=>0) - (1 + asin(δ))^2 / (ϵ * κ^2) * subs(subs(diff(ψ, x[1]), x[1]=>1+ϵ), x[2]=>0)), - # expand(subs(subs(diff(ψ, x[2], 2), x[1]=>1-ϵ), x[2]=>0) + (1 - asin(δ))^2 / (ϵ * κ^2) * subs(subs(diff(ψ, x[1]), x[1]=>1-ϵ), x[2]=>0)), - # expand(subs(subs(diff(ψ, x[1], 2), x[1]=>1-δ*ϵ), x[2]=>κ*ϵ) - κ / (ϵ * (1 - δ^2)) * subs(subs(diff(ψ, x[2]), x[1]=>1-δ*ϵ), x[2]=>κ*ϵ)) - # ] - - csym = solve(eqs, c) - cnum = [N(csym[c[i]]) for i in 1:n] - - SolovevXpoint{T}(R₀, B₀, ϵ, κ, δ, a, xₛₑₚ, yₛₑₚ, cnum) -end + function SolovevXpointEquilibrium(R₀::T, B₀::T, ϵ::T, κ::T, δ::T, a::T, xₛₑₚ::T, yₛₑₚ::T) where T <: Number + + n = 12 + + x = [Sym("x" * string(i)) for i in 1:3] + c = [Sym("c" * string(i)) for i in 1:n] + + ψ = (ψ₀(x,a) + c[1] * ψ₁(x) + + c[2] * ψ₂(x) + + c[3] * ψ₃(x) + + c[4] * ψ₄(x) + + c[5] * ψ₅(x) + + c[6] * ψ₆(x) + + c[7] * ψ₇(x) + + c[8] * ψ₈(x) + + c[9] * ψ₉(x) + + c[10] * ψ₁₀(x) + + c[11] * ψ₁₁(x) + + c[12] * ψ₁₂(x) ) + + eqs = [ + expand(subs(subs(ψ, x[1], 1+ϵ), x[2], 0)), + expand(subs(subs(ψ, x[1], 1-ϵ), x[2], 0)), + expand(subs(subs(ψ, x[1], 1-δ*ϵ), x[2], κ*ϵ)), + expand(subs(subs(ψ, x[1], xₛₑₚ), x[2], yₛₑₚ)), + expand(subs(subs(diff(ψ, x[2]), x[1], 1+ϵ), x[2], 0)), + expand(subs(subs(diff(ψ, x[2]), x[1], 1-ϵ), x[2], 0)), + expand(subs(subs(diff(ψ, x[1]), x[1], 1-δ*ϵ), x[2], κ*ϵ)), + expand(subs(subs(diff(ψ, x[1]), x[1], xₛₑₚ), x[2], yₛₑₚ)), + expand(subs(subs(diff(ψ, x[2]), x[1], xₛₑₚ), x[2], yₛₑₚ)), + expand(subs(subs(diff(ψ, x[2], 2), x[1], 1+ϵ), x[2], 0) - (1 + asin(δ))^2 / (ϵ * κ^2) * subs(subs(diff(ψ, x[1]), x[1], 1+ϵ), x[2], 0)), + expand(subs(subs(diff(ψ, x[2], 2), x[1], 1-ϵ), x[2], 0) + (1 - asin(δ))^2 / (ϵ * κ^2) * subs(subs(diff(ψ, x[1]), x[1], 1-ϵ), x[2], 0)), + expand(subs(subs(diff(ψ, x[1], 2), x[1], 1-δ*ϵ), x[2], κ*ϵ) - κ / (ϵ * (1 - δ^2)) * subs(subs(diff(ψ, x[2]), x[1], 1-δ*ϵ), x[2], κ*ϵ)) + ] + + csym = solve(eqs, c) + cnum = [N(csym[c[i]]) for i in 1:n] + + SolovevXpointEquilibrium{T}(R₀, B₀, ϵ, κ, δ, a, xₛₑₚ, yₛₑₚ, cnum) + end -SolovevXpointITER() = SolovevXpoint(6.2, 5.3, 0.32, 1.7, 0.33, -0.155, 0.88, -0.60) -SolovevXpointNSTX() = SolovevXpoint(0.85, 0.3, 0.78, 2.0, 0.35, -0.05, 0.70, -1.71) + SolovevXpointEquilibriumITER() = SolovevXpointEquilibrium(6.2, 5.3, 0.32, 1.7, 0.33, -0.155, 0.88, -0.60) + SolovevXpointEquilibriumNSTX() = SolovevXpointEquilibrium(0.85, 0.3, 0.78, 2.0, 0.35, -0.05, 0.70, -1.71) -function Base.show(io::IO, equ::SolovevXpoint) - print(io, "Solovev Xpoint Equilibrium with\n") - print(io, " R₀ = ", equ.R₀, "\n") - print(io, " B₀ = ", equ.B₀, "\n") - print(io, " ϵ = ", equ.ϵ, "\n") - print(io, " κ = ", equ.κ, "\n") - print(io, " δ = ", equ.δ, "\n") - print(io, " a = ", equ.a, "\n") - print(io, " xₛₑₚ = ", equ.xₛₑₚ, "\n") - print(io, " yₛₑₚ = ", equ.yₛₑₚ) -end + function Base.show(io::IO, equ::SolovevXpointEquilibrium) + print(io, "Solovev Xpoint Equilibrium with\n") + print(io, " R₀ = ", equ.R₀, "\n") + print(io, " B₀ = ", equ.B₀, "\n") + print(io, " ϵ = ", equ.ϵ, "\n") + print(io, " κ = ", equ.κ, "\n") + print(io, " δ = ", equ.δ, "\n") + print(io, " a = ", equ.a, "\n") + print(io, " xₛₑₚ = ", equ.xₛₑₚ, "\n") + print(io, " yₛₑₚ = ", equ.yₛₑₚ) + end -function A₃(x::AbstractArray{T,1}, equ::SolovevXpoint) where {T <: Number} - (ψ₀(x, equ.a) + equ.c[1] * ψ₁(x) - + equ.c[2] * ψ₂(x) - + equ.c[3] * ψ₃(x) - + equ.c[4] * ψ₄(x) - + equ.c[5] * ψ₅(x) - + equ.c[6] * ψ₆(x) - + equ.c[7] * ψ₇(x) - + equ.c[8] * ψ₈(x) - + equ.c[9] * ψ₉(x) - + equ.c[10] * ψ₁₀(x) - + equ.c[11] * ψ₁₁(x) - + equ.c[12] * ψ₁₂(x) ) -end + function ElectromagneticFields.A₃(x::AbstractArray{T,1}, equ::SolovevXpointEquilibrium) where {T <: Number} + ψ₀(x, equ.a) + equ.c[1] * ψ₁(x) + + equ.c[2] * ψ₂(x) + + equ.c[3] * ψ₃(x) + + equ.c[4] * ψ₄(x) + + equ.c[5] * ψ₅(x) + + equ.c[6] * ψ₆(x) + + equ.c[7] * ψ₇(x) + + equ.c[8] * ψ₈(x) + + equ.c[9] * ψ₉(x) + + equ.c[10] * ψ₁₀(x) + + equ.c[11] * ψ₁₁(x) + + equ.c[12] * ψ₁₂(x) + end + + macro solovev_xpoint_equilibrium(R₀, B₀, ϵ, κ, δ, a, xₛₑₚ, yₛₑₚ) + generate_equilibrium_code(SolovevXpointEquilibrium(R₀, B₀, ϵ, κ, δ, a, xₛₑₚ, yₛₑₚ); output=false) + end + + function init(R₀, B₀, ϵ, κ, δ, a, xₛₑₚ, yₛₑₚ; perturbation=ZeroPerturbation()) + equilibrium = SolovevXpointEquilibrium(R₀, B₀, ϵ, κ, δ, a, xₛₑₚ, yₛₑₚ) + load_equilibrium(equilibrium, perturbation; target_module=SolovevXpoint) + return equilibrium + end -macro solovev_xpoint_equilibrium(R₀, B₀, ϵ, κ, δ, a, xₛₑₚ, yₛₑₚ) - generate_equilibrium_code(SolovevXpoint(R₀, B₀, ϵ, κ, δ, a, xₛₑₚ, yₛₑₚ); output=false) end diff --git a/src/analytic/symmetric_quadratic.jl b/src/analytic/symmetric_quadratic.jl index 5d8ed37..970fd7a 100644 --- a/src/analytic/symmetric_quadratic.jl +++ b/src/analytic/symmetric_quadratic.jl @@ -15,32 +15,51 @@ B(x,y,z) = B_0 \, \begin{pmatrix} Parameters: `B₀` """ -struct SymmetricQuadratic{T <: Number} <: CartesianEquilibrium - name::String - B₀::T - SymmetricQuadratic{T}(B₀::T) where T <: Number = new("SymmetricQuadraticEquilibrium", B₀) -end +module SymmetricQuadratic -SymmetricQuadratic(B₀::T) where T <: Number = SymmetricQuadratic{T}(B₀) + import ..ElectromagneticFields + import ..ElectromagneticFields: CartesianEquilibrium, ZeroPerturbation + import ..ElectromagneticFields: load_equilibrium, generate_equilibrium_code + import ..AnalyticCartesianField: X, Y, Z + export SymmetricQuadraticEquilibrium -function Base.show(io::IO, equ::SymmetricQuadratic) - print(io, "Symmetric Quadratic Equilibrium") -end + const DEFAULT_B₀ = 1.0 + + struct SymmetricQuadraticEquilibrium{T <: Number} <: CartesianEquilibrium + name::String + B₀::T + SymmetricQuadraticEquilibrium{T}(B₀::T) where T <: Number = new("SymmetricQuadraticEquilibrium", B₀) + end + + SymmetricQuadraticEquilibrium(B₀::T=DEFAULT_B₀) where T <: Number = SymmetricQuadraticEquilibrium{T}(B₀) + + + function Base.show(io::IO, equ::SymmetricQuadraticEquilibrium) + print(io, "Symmetric Quadratic Equilibrium") + end + + + r²(x::AbstractVector, equ::SymmetricQuadraticEquilibrium) = X(x,equ)^2 + Y(x,equ)^2 + r(x::AbstractVector, equ::SymmetricQuadraticEquilibrium) = sqrt(r²(x,equ)) + R(x::AbstractVector, equ::SymmetricQuadraticEquilibrium) = r(x,equ) + θ(x::AbstractVector, equ::SymmetricQuadraticEquilibrium) = atan(Y(x,equ), X(x,equ)) + ϕ(x::AbstractVector, equ::SymmetricQuadraticEquilibrium) = θ(x,equ) + ElectromagneticFields.A₁(x::AbstractVector, equ::SymmetricQuadraticEquilibrium) = - equ.B₀ * x[2] * (2 + x[1]^2 + x[2]^2) / 4 + ElectromagneticFields.A₂(x::AbstractVector, equ::SymmetricQuadraticEquilibrium) = + equ.B₀ * x[1] * (2 + x[1]^2 + x[2]^2) / 4 + ElectromagneticFields.A₃(x::AbstractVector, equ::SymmetricQuadraticEquilibrium) = zero(eltype(x)) -r²(x::AbstractVector, equ::SymmetricQuadratic) = X(x,equ)^2 + Y(x,equ)^2 -r(x::AbstractVector, equ::SymmetricQuadratic) = sqrt(r²(x,equ)) -R(x::AbstractVector, equ::SymmetricQuadratic) = r(x,equ) -θ(x::AbstractVector, equ::SymmetricQuadratic) = atan(Y(x,equ), X(x,equ)) -ϕ(x::AbstractVector, equ::SymmetricQuadratic) = θ(x,equ) + ElectromagneticFields.get_functions(::SymmetricQuadraticEquilibrium) = (X=X, Y=Y, Z=Z, R=R, r=r, θ=θ, ϕ=ϕ, r²=r²) -A₁(x::AbstractVector, equ::SymmetricQuadratic) = - equ.B₀ * x[2] * (2 + x[1]^2 + x[2]^2) / 4 -A₂(x::AbstractVector, equ::SymmetricQuadratic) = + equ.B₀ * x[1] * (2 + x[1]^2 + x[2]^2) / 4 -A₃(x::AbstractVector, equ::SymmetricQuadratic) = zero(eltype(x)) + macro symmetric_quadratic_equilibrium(B₀=DEFAULT_B₀) + generate_equilibrium_code(SymmetricQuadraticEquilibrium(B₀); output=false) + end -get_functions(::SymmetricQuadratic) = (X=X, Y=Y, Z=Z, R=R, r=r, θ=θ, ϕ=ϕ, r²=r²) + function init(B₀=DEFAULT_B₀; perturbation=ZeroPerturbation()) + equilibrium = SymmetricQuadraticEquilibrium(B₀) + load_equilibrium(equilibrium, perturbation; target_module=SymmetricQuadratic) + return equilibrium + end -macro symmetric_quadratic_equilibrium() - generate_equilibrium_code(SymmetricQuadratic(1.); output=false) end diff --git a/src/analytic/theta_pinch.jl b/src/analytic/theta_pinch.jl index 2f28b08..8f085a1 100644 --- a/src/analytic/theta_pinch.jl +++ b/src/analytic/theta_pinch.jl @@ -12,38 +12,56 @@ B (x,y,z) = \big( 0 , \, 0 , \, B_0 \big)^T . Parameters: B₀: B-field at magnetic axis """ -struct ThetaPinch{T <: Number} <: CartesianEquilibrium - name::String - B₀::T +module ThetaPinch - function ThetaPinch{T}(B₀::T) where T <: Number - new("ThetaPinchEquilibrium", B₀) + import ..ElectromagneticFields + import ..ElectromagneticFields: CartesianEquilibrium, ZeroPerturbation + import ..ElectromagneticFields: load_equilibrium, generate_equilibrium_code + import ..AnalyticCartesianField: X, Y, Z + + export ThetaPinchEquilibrium + + const DEFAULT_B₀ = 1.0 + + struct ThetaPinchEquilibrium{T <: Number} <: CartesianEquilibrium + name::String + B₀::T + + function ThetaPinchEquilibrium{T}(B₀::T) where T <: Number + new("ThetaPinchEquilibrium", B₀) + end end -end -ThetaPinch(B₀::T=1.0) where T <: Number = ThetaPinch{T}(B₀) + ThetaPinchEquilibrium(B₀::T=DEFAULT_B₀) where T <: Number = ThetaPinchEquilibrium{T}(B₀) -function Base.show(io::IO, equ::ThetaPinch) - print(io, "θ-Pinch Equilibrium in (x,y,z) Coordinates with\n") - print(io, " B₀ = ", equ.B₀) -end + function Base.show(io::IO, equ::ThetaPinchEquilibrium) + print(io, "θ-Pinch Equilibrium in (x,y,z) Coordinates with\n") + print(io, " B₀ = ", equ.B₀) + end -R(x::AbstractArray{T,1}, equ::ThetaPinch) where {T <: Number} = r(x,equ) -r(x::AbstractArray{T,1}, equ::ThetaPinch) where {T <: Number} = sqrt(r²(x,equ)) -r²(x::AbstractArray{T,1}, equ::ThetaPinch) where {T <: Number} = X(x,equ)^2 + Y(x,equ)^2 -θ(x::AbstractArray{T,1}, equ::ThetaPinch) where {T <: Number} = atan(Y(x,equ), X(x,equ)) -ϕ(x::AbstractArray{T,1}, equ::ThetaPinch) where {T <: Number} = θ(x,equ) + r²(x::AbstractVector, equ::ThetaPinchEquilibrium) = X(x,equ)^2 + Y(x,equ)^2 + r(x::AbstractVector, equ::ThetaPinchEquilibrium) = sqrt(r²(x,equ)) + R(x::AbstractVector, equ::ThetaPinchEquilibrium) = r(x,equ) + θ(x::AbstractVector, equ::ThetaPinchEquilibrium) = atan(Y(x,equ), X(x,equ)) + ϕ(x::AbstractVector, equ::ThetaPinchEquilibrium) = θ(x,equ) -A₁(x::AbstractArray{T,1}, equ::ThetaPinch) where {T <: Number} = - equ.B₀ * Y(x,equ) / 2 -A₂(x::AbstractArray{T,1}, equ::ThetaPinch) where {T <: Number} = + equ.B₀ * X(x,equ) / 2 -A₃(x::AbstractArray{T,1}, equ::ThetaPinch) where {T <: Number} = zero(eltype(x)) + ElectromagneticFields.A₁(x::AbstractVector, equ::ThetaPinchEquilibrium) = - equ.B₀ * Y(x,equ) / 2 + ElectromagneticFields.A₂(x::AbstractVector, equ::ThetaPinchEquilibrium) = + equ.B₀ * X(x,equ) / 2 + ElectromagneticFields.A₃(x::AbstractVector, equ::ThetaPinchEquilibrium) = zero(eltype(x)) + ElectromagneticFields.get_functions(::ThetaPinchEquilibrium) = (X=X, Y=Y, Z=Z, R=R, r=r, θ=θ, ϕ=ϕ, r²=r²) -get_functions(::ThetaPinch) = (X=X, Y=Y, Z=Z, R=R, r=r, θ=θ, ϕ=ϕ, r²=r²) + macro zpinch_equilibrium(B₀=DEFAULT_B₀) + generate_equilibrium_code(ThetaPinchEquilibrium(B₀); output=false) + end + + function init(B₀=DEFAULT_B₀; perturbation=ZeroPerturbation()) + equilibrium = ThetaPinchEquilibrium(B₀) + load_equilibrium(equilibrium, perturbation; target_module=ThetaPinch) + return equilibrium + end -macro zpinch_equilibrium(R₀, B₀) - generate_equilibrium_code(ThetaPinch(R₀, B₀); output=false) end diff --git a/test/test_analytic.jl b/test/test_analytic.jl index 3bcd6fe..beee3c8 100644 --- a/test/test_analytic.jl +++ b/test/test_analytic.jl @@ -77,39 +77,24 @@ function test_equilibrium(equ_mod, t, x) end -# modules that will hold the generated equilibrium code -module AxisymmetricTokamakCartesianEquilibrium end -module AxisymmetricTokamakCylindricalEquilibrium end -module AxisymmetricTokamakToroidalEquilibrium end -module SymmetricQuadraticEquilibrium end -module ThetaPinchEquilibrium end -module ABCEquilibrium end -module SolovevEquilibrium end -module SolovevXpointEquilibrium end -module SolovevQuadraticEquilibrium end - -module SymmetricQuadraticEquilibriumEzCosZPerturbation end -module ThetaPinchEquilibriumEzCosZPerturbation end - - # equilibrium list (equilibrium, parameters, periodicity, module) eqs = ( - (AxisymmetricTokamakCartesian, (2., 3., 2.), [0., 0., 0.], AxisymmetricTokamakCartesianEquilibrium), - (AxisymmetricTokamakCylindrical, (2., 3., 2.), [0., 0., 2π], AxisymmetricTokamakCylindricalEquilibrium), - (AxisymmetricTokamakToroidal, (2., 3., 2.), [0., 2π, 2π], AxisymmetricTokamakToroidalEquilibrium), - (SymmetricQuadratic, (1.), [0., 0., 0.], SymmetricQuadraticEquilibrium), - (ThetaPinch, (1.), [0., 0., 0.], ThetaPinchEquilibrium), - (ABC, (1., 0.5, 0.5), [0., 0., 0.], ABCEquilibrium), - (Solovev, (6.2, 5.3, 0.32, 1.8, 0.45, -0.155), [0., 0., 2π], SolovevEquilibrium), - (SolovevXpoint, (6.2, 5.3, 0.32, 1.8, 0.45, -0.155, 0.88, -0.60), [0., 0., 2π], SolovevXpointEquilibrium), - (SolovevQuadratic, (6.2, 5.3, 1., 1.), [0., 0., 2π], SolovevQuadraticEquilibrium), + (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.]), + (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π]), ) # perturbation list (equilibrium, parameters, perturbation, parameters, module) perts = ( - (SymmetricQuadratic, (1.), EzCosZ, (2.), SymmetricQuadraticEquilibriumEzCosZPerturbation), - (ThetaPinch, (1.), EzCosZ, (2.), ThetaPinchEquilibriumEzCosZPerturbation), + (SymmetricQuadratic, (1.), EzCosZ, (2.)), + (ThetaPinch, (1.), EzCosZ, (2.)), ) @@ -119,9 +104,9 @@ x = [1., 1., 1.] # test equilibria @testset "$(rpad(teststring(equ[1]),60))" for equ in eqs begin - equ_obj = equ[1](equ[2]...) - load_equilibrium(equ_obj, target_module=equ[4]) - test_equilibrium(equ[4], t, x) + # @eval import .$(Symbol(equ[1])) + equ_obj = equ[1].init(equ[2]...) + test_equilibrium(equ[1], t, x) @test periodicity(x, equ_obj) == equ[3] end end @@ -130,10 +115,10 @@ println() # test perturbations @testset "$(rpad(teststring(equ[1], equ[3]),60))" for equ in perts begin - equ_obj = equ[1](equ[2]...) - prt_obj = equ[3](equ[4]...) - load_equilibrium(equ_obj, prt_obj, target_module=equ[5]) - test_equilibrium(equ[5], t, x) + # @eval import .$(Symbol(equ[1])) + # @eval import .$(Symbol(equ[3])) + equ_obj = equ[1].init(equ[2]..., perturbation=equ[3].init(equ[4]...)) + test_equilibrium(equ[1], t, x) end end println() @@ -227,11 +212,11 @@ function test_abc_equilibrium(equ_mod, t=0., x=[1.0, 0.5, 0.5]) end @testset "$(rpad("Magnetic Fields",60))" begin - test_axisymmetric_tokamak_cartesian_equilibrium(AxisymmetricTokamakCartesianEquilibrium) - test_axisymmetric_tokamak_cylindrical_equilibrium(AxisymmetricTokamakCylindricalEquilibrium) - test_axisymmetric_tokamak_toroidal_equilibrium(AxisymmetricTokamakToroidalEquilibrium) - test_symmetric_quadratic_equilibrium(SymmetricQuadraticEquilibrium) - test_theta_pinch_equilibrium(ThetaPinchEquilibrium) - test_abc_equilibrium(ABCEquilibrium) + test_axisymmetric_tokamak_cartesian_equilibrium(AxisymmetricTokamakCartesian) + 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 println()