Skip to content

Commit

Permalink
Rename toroidal tokamak equilibrium to circular tokamak equilibrium. …
Browse files Browse the repository at this point in the history
…Add toroidally regularised equilibrium.
  • Loading branch information
michakraus committed Sep 24, 2020
1 parent af12b46 commit ba68b01
Show file tree
Hide file tree
Showing 4 changed files with 171 additions and 65 deletions.
21 changes: 15 additions & 6 deletions src/ElectromagneticFields.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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")
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -15,81 +15,81 @@ 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")
print(io, " q₀ = ", equ.q₀)
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.::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = X(ξ,equ)
ElectromagneticFields.::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = Y(ξ,equ)
ElectromagneticFields.::AbstractVector, equ::AxisymmetricTokamakToroidalEquilibrium) = Z(ξ,equ)
ElectromagneticFields.::AbstractVector, equ::AxisymmetricTokamakCircularEquilibrium) = X(ξ,equ)
ElectromagneticFields.::AbstractVector, equ::AxisymmetricTokamakCircularEquilibrium) = Y(ξ,equ)
ElectromagneticFields.::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π
return p
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

Expand Down
96 changes: 96 additions & 0 deletions src/analytic/axisymmetric_tokamak_toroidal_regularization.jl
Original file line number Diff line number Diff line change
@@ -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.::AbstractVector, equ::AxisymmetricTokamakToroidalRegularizationEquilibrium) = X(ξ,equ)
ElectromagneticFields.::AbstractVector, equ::AxisymmetricTokamakToroidalRegularizationEquilibrium) = Y(ξ,equ)
ElectromagneticFields.::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
Loading

0 comments on commit ba68b01

Please sign in to comment.