Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Zi temperature dependent and normalizations #6

Merged
merged 2 commits into from
Oct 16, 2024
Merged

Zi temperature dependent and normalizations #6

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
1a77f90
Zi temperature dep and normalizations
orso82 Oct 16, 2024
30d2a0f
Minor bug fixes
adrianaghiozzi Oct 16, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 1 addition & 2 deletions src/NEO.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,9 +65,8 @@ include("chang_hinton.jl")
"""
run_neo(input_neo::InputNEO)

Saves input.neo file to a temporary directory, runs NEO on that directory and parses output
Saves input.neo file to a temporary directory, runs NEO on that directory and parses output
"""

function run_neo(input_neo::InputNEO)
folder = mktempdir()
save_inputneo(input_neo, joinpath(folder, "input.neo"))
Expand Down
20 changes: 12 additions & 8 deletions src/chang_hinton.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,17 +16,23 @@ function changhinton(
iion::Integer)

eq1d = eqt.profiles_1d
ions = cp1d.ion

e = IMAS.gacode_units.e # statcoul
k = IMAS.gacode_units.k # erg/eV
mp = IMAS.gacode_units.mp # g
m_to_cm = IMAS.gacode_units.m_to_cm
m³_to_cm³ = IMAS.gacode_units.m³_to_cm³

m_to_cm = 1e2
Rmaj0 = eq1d.geometric_axis.r[1] * m_to_cm

rho_eq = eq1d.rho_tor_norm
rho_cp = cp1d.grid.rho_tor_norm
gridpoint_cp = argmin(abs.(rho_cp .- rho_fluxmatch))

ne = cp1d.electrons.density_thermal * 1e-6
ne = cp1d.electrons.density_thermal / m³_to_cm³
Te = cp1d.electrons.temperature[gridpoint_cp]
Ti = cp1d.ion[iion].temperature
Ti = ions[iion].temperature

Rmaj = IMAS.interp1d(eq1d.rho_tor_norm, 0.5 * (eq1d.r_outboard .+ eq1d.r_inboard)).(cp1d.grid.rho_tor_norm)
rmin = IMAS.interp1d(rho_eq, 0.5 * (eq1d.r_outboard - eq1d.r_inboard)).(rho_cp)
Expand All @@ -43,10 +49,8 @@ function changhinton(
dlnndr = dlnndr * a / m_to_cm
Ti = Ti[gridpoint_cp]
ne = ne[gridpoint_cp]
k = 1.6e-12
e = 4.8e-10

mi = 1.6726e-24 * cp1d.ion[iion].element[1].a
mi = ions[iion].element[1].a * mp

c_s = sqrt(k * Te / mi)
loglam = 24.0 - log(sqrt(ne / Te))
Expand All @@ -56,7 +60,7 @@ function changhinton(
b0 = 0.31
c0 = 0.74

Zi = IMAS.avgZ(cp1d.ion[iion].element[1].z_n, Ti)
Zi = IMAS.avgZ(ions[iion].element[1].z_n, Ti)
alpha = Zi - 1.0
nui = sqrt(2.0) * pi * ne * (Zi * e)^4 * loglam / sqrt(mi) / (k * Ti)^1.5
nu = nui * a / c_s / sqrt(Ti / Ti) * (Ti / Te)^1.5
Expand Down Expand Up @@ -99,7 +103,7 @@ function changhinton(
efluxi = (CH_I_div_psip^2
*
Ti / Te
* (cp1d.ion[iion].element[1].a / Zi^2 * neo_rho_star_in^2 * sqrt(eps) * nui_HH)
* (ions[iion].element[1].a / Zi^2 * neo_rho_star_in^2 * sqrt(eps) * nui_HH)
* ((K2 + K1) * dlntdr + K1 * dlnndr))

qneo_gb = neo_rho_star_in^2
Expand Down
56 changes: 27 additions & 29 deletions src/hirshman_sigmar.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,6 @@ end

Populates equilibrium_geometry structure with equilibrium quantities from dd
"""

function get_equilibrium_parameters(dd::IMAS.dd)
equilibrium_geometry = NEO.equilibrium_geometry()

Expand Down Expand Up @@ -58,9 +57,8 @@ end
"""
get_ion_electron_parameters(dd::IMAS.dd)

Populates parameter_matrices structure with profile data from dd using NEO normalizations
Populates parameter_matrices structure with profile data from dd using NEO normalizations
"""

function get_ion_electron_parameters(dd::IMAS.dd)
parameter_matrices = NEO.parameter_matrices()

Expand All @@ -75,54 +73,55 @@ function get_ion_electron_parameters(dd::IMAS.dd)

e = IMAS.gacode_units.e # statcoul
k = IMAS.gacode_units.k # erg/eV
mp = IMAS.constants.m_p * 1e3 # g
md = 3.34358e-27 * 1e3
mp = IMAS.gacode_units.mp # g
me = IMAS.gacode_units.me # g
md = 2 * mp # g
m³_to_cm³ = IMAS.gacode_units.m³_to_cm³

loglam = 24.0 .- log.(sqrt.(cp1d.electrons.density ./ 1e6) ./ (cp1d.electrons.temperature))
Te = cp1d.electrons.temperature # ev
ne = cp1d.electrons.density / m³_to_cm³ # cm^-3

n_norm = cp1d.electrons.density ./ 1e6
t_norm = cp1d.electrons.temperature
loglam = 24.0 .- log.(sqrt.(ne) ./ Te)

m_norm = 2.0
nu_norm = sqrt.(k .* cp1d.electrons.temperature ./ md) ./ a
n_norm = ne
m_norm = md
t_norm = Te
nu_norm = sqrt.(k .* Te ./ md) ./ a

Z = Vector{Float64}(undef, num_ions + 1)
mass = Vector{Float64}(undef, num_ions + 1)

dens = zeros(Float64, n, num_ions + 1)
temp = zeros(Float64, n, num_ions + 1)
vth = zeros(Float64, n, num_ions + 1)
nu = zeros(Float64, n, num_ions + 1)
dlnndr = zeros(Float64, n, num_ions + 1)
dlntdr = zeros(Float64, n, num_ions + 1)

for i in 1:num_ions
Z[i] = cp1d.ion[i].element[1].z_n
mass[i] = cp1d.ion[i].element[1].a ./ m_norm
T1 = cp1d.ion[i].temperature[Int(ceil(end / 2))]
Z[i] = IMAS.avgZ(cp1d.ion[i].element[1].z_n, T1)
mass[i] = cp1d.ion[i].element[1].a * mp / m_norm

dens[:, i] = cp1d.ion[i].density ./ 1e6 ./ n_norm
dens[:, i] = cp1d.ion[i].density ./ m³_to_cm³ ./ n_norm
temp[:, i] = cp1d.ion[i].temperature ./ t_norm

nu[:, i] = (@. sqrt(2) * pi * (cp1d.ion[i].density ./ 1e6) * Z[i]^4.0 * e^4.0 * loglam / sqrt(cp1d.ion[i].element[1].a * mp) / (k * cp1d.ion[i].temperature)^1.5) ./ nu_norm
nu[:, i] =
(@. sqrt(2) * pi * (cp1d.ion[i].density ./ m³_to_cm³) * Z[i]^4.0 * e^4.0 * loglam / sqrt(cp1d.ion[i].element[1].a * mp) / (k * cp1d.ion[i].temperature)^1.5) ./ nu_norm

dlnndr[:, i] = -IMAS.calc_z(rmin / a, cp1d.ion[i].density, :backward)
dlntdr[:, i] = -IMAS.calc_z(rmin / a, cp1d.ion[i].temperature, :backward)

vth[:, i] = sqrt.((cp1d.ion[i].temperature ./ t_norm) ./ (cp1d.ion[i].element[1].a ./ m_norm))
vth[:, i] = sqrt.((cp1d.ion[i].temperature ./ t_norm) ./ (cp1d.ion[i].element[1].a * mp / m_norm))
end

# tacking on electron parameters at the end
Z[end] = -1.0
mass[end] = 0.00054858 / m_norm # 0.00054858 is the mass of an electron in AMU
dens[:, end] = cp1d.electrons.density ./ 1e6 ./ n_norm
temp[:, end] = cp1d.electrons.temperature ./ t_norm

nu[:, end] = (@. sqrt(2) * pi * (cp1d.electrons.density ./ 1e6) * (Z[end] * e)^4.0 * loglam / sqrt(0.00054858 * mp) / (k .* cp1d.electrons.temperature)^1.5) ./ nu_norm

dlnndr[:, end] = -IMAS.calc_z(rmin / a, cp1d.electrons.density, :backward)
dlntdr[:, end] = -IMAS.calc_z(rmin / a, cp1d.electrons.temperature, :backward)

vth[:, end] = sqrt.((cp1d.electrons.temperature ./ t_norm) ./ (0.00054858 ./ m_norm))
mass[end] = me / md
dens[:, end] = ne ./ n_norm
temp[:, end] = Te ./ t_norm
nu[:, end] = (@. sqrt(2) * pi * ne * (Z[end] * e)^4.0 * loglam / sqrt(me) / (k .* Te)^1.5) ./ nu_norm
dlnndr[:, end] = -IMAS.calc_z(rmin / a, ne, :backward)
dlntdr[:, end] = -IMAS.calc_z(rmin / a, Te, :backward)
vth[:, end] = sqrt.((Te ./ t_norm) ./ (me ./ md))

parameter_matrices.Z = Z
parameter_matrices.mass = mass
Expand Down Expand Up @@ -199,8 +198,7 @@ function gauss_integ(
ietype::Int,
equilibrium_geometry::NEO.equilibrium_geometry,
is_globalHS::Int,
ir_global::Int
)
ir_global::Int)

x0, w0 = gauss_legendre(0, 1, order)

Expand Down
36 changes: 16 additions & 20 deletions src/input_neo.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -221,9 +221,7 @@ end
InputNEO(dd::IMAS.dd, gridpoint_cp)

Populates InputNEO structure with quantities from dd using NEO normalizations

"""

function InputNEO(dd::IMAS.dd, gridpoint_cp)
input_neo = InputNEO()

Expand All @@ -233,24 +231,22 @@ function InputNEO(dd::IMAS.dd, gridpoint_cp)
cp1d = dd.core_profiles.profiles_1d[]
ions = cp1d.ion

mp = IMAS.constants.m_p * 1e3 # g
me = IMAS.constants.m_e * 1e3

e = IMAS.gacode_units.e # statcoul
k = IMAS.gacode_units.k # erg/eV
mp = IMAS.gacode_units.mp # g
me = IMAS.gacode_units.me # g
md = 2 * mp # g
m_to_cm = IMAS.gacode_units.m_to_cm
m³_to_cm³ = IMAS.gacode_units.m³_to_cm³

rmin = IMAS.r_min_core_profiles(cp1d, eqt)
a = rmin[end]

e = IMAS.gacode_units.e # statcoul
k = IMAS.gacode_units.k # erg/eV

temp_1 = ions[1].temperature
T1 = temp_1[gridpoint_cp]
dens_1 = ions[1].density[gridpoint_cp] ./ m³_to_cm³

dens_1 = ions[1].density ./ 1e6
n1 = dens_1[gridpoint_cp]

dens_e = cp1d.electrons.density ./ 1e6
dens_e = cp1d.electrons.density ./ m³_to_cm³
dlnnedr = -IMAS.calc_z(rmin ./ a, dens_e, :backward)
ne = dens_e[gridpoint_cp]
dlnnedr = dlnnedr[gridpoint_cp]
Expand All @@ -261,7 +257,7 @@ function InputNEO(dd::IMAS.dd, gridpoint_cp)
dlntedr = dlntedr[gridpoint_cp]

n_norm = ne
m_norm = 3.3452e-27 * 1e3 # mass of deuterium in grams
m_norm = md
t_norm = Te
v_norm = sqrt(k .* t_norm ./ m_norm)

Expand All @@ -275,15 +271,15 @@ function InputNEO(dd::IMAS.dd, gridpoint_cp)
loglam = IMAS.lnΛ_ei(cp1d.electrons.density[gridpoint_cp], cp1d.electrons.temperature[gridpoint_cp])
Z1 = IMAS.avgZ(ions[1].element[1].z_n, T1)
m1 = ions[1].element[1].a * mp
nu1 = @. sqrt(2) * pi * dens_1 * Z1^4.0 * e^4.0 * loglam / (sqrt(m1) * (k * temp_1)^1.5)
nu1 = sqrt(2) * pi * dens_1 * Z1^4.0 * e^4.0 * loglam ./ (sqrt(m1) * (k * temp_1).^1.5)

input_neo.NU_1 = (nu1/(v_norm./a))[gridpoint_cp]
input_neo.NU_1 = (nu1 ./ (v_norm ./ a))[gridpoint_cp]

w0 = -1 .* cp1d.rotation_frequency_tor_sonic
w0p = IMAS.gradient(rmin, w0)
w0 = -cp1d.rotation_frequency_tor_sonic[gridpoint_cp]
w0p = -IMAS.gradient(rmin, cp1d.rotation_frequency_tor_sonic)[gridpoint_cp]

input_neo.OMEGA_ROT = w0[gridpoint_cp] / (v_norm / a)
input_neo.OMEGA_ROT_DERIV = w0p[gridpoint_cp] * a^2 / v_norm
input_neo.OMEGA_ROT = w0 / (v_norm / a)
input_neo.OMEGA_ROT_DERIV = w0p * a^2 / v_norm

####################

Expand Down Expand Up @@ -334,7 +330,7 @@ function InputNEO(dd::IMAS.dd, gridpoint_cp)
Ti = Ti[gridpoint_cp]
dlntidr = dlntidr[gridpoint_cp]

ni = ions[iion].density ./ 1e6 / n_norm
ni = ions[iion].density ./ m³_to_cm³ / n_norm
dlnnidr = -IMAS.calc_z(rmin ./ a, ni, :backward)
ni = ni[gridpoint_cp]
dlnnidr = dlnnidr[gridpoint_cp]
Expand Down
Loading