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

Improve algorithm and typing of ActorNeutronics #167

Merged
merged 1 commit into from
Oct 12, 2022
Merged

Improve algorithm and typing of ActorNeutronics #167

Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
235 changes: 186 additions & 49 deletions src/actors/neutronics_actor.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,16 +19,18 @@ function Zcoord(n::neutron_particle)
n.z
end

mutable struct ActorNeutronics <: PlasmaAbstractActor
dd::IMAS.dd
mutable struct ActorNeutronics{T} <: PlasmaAbstractActor
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@bclyons12 you raise a good point here about the type of IMAS.dd. This is something I had not originally caught!

In addition of having to define the actors constructors as you have done, we should do the same for all functions that use type annotations. That is ~500 hits across 58 files. Doable but still annoying.

I'll merge this, and actually remove the {T} from the constructor (I just don't like to have the source of an actor being different than the others. I'll open a separate issue specifically to discuss this and decide a course of action.

Copy link
Member Author

@bclyons12 bclyons12 Oct 12, 2022 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@orso82 You need to do this for structs, since IMAS.dd is abstract. You don't need to do it for function arguments though. Since IMAS.dd{T} <: IMAS.dd, having ::IMAS.dd will allow the compiler will dispatch on each separately.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah, ok! you are right! Phew 😅

What threw me off is that you added this constructor.

ActorNeutronics(dd::IMAS.dd{T}, par::ParametersActor; kw...) where {T} = ActorNeutronics{T}(dd, par; kw...)

Is this really necessary? shouldn't ActorNeutronics pick up the right {T} already based on the concrete IMAS.dd type that is passed?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Only if you don't override the default inner constructor. You have here (something I advocate against 😃 ), so it doesn't make the default outer constructor either.

Here's the relevant bit from https://docs.julialang.org/en/v1/manual/constructors/ :
Screen Shot 2022-10-11 at 9 25 45 PM

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You actually might need the inner constructor here, since you use it to accept keywords and Julia doesn't multiple dispatch on keywords. If you made that as an outer constructor, it couldn't distinguish between that and the default inner constructor...

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah, that's also something I had not caught.

As you noted, all actors that only have dd and par as their structure need to use a inner constructor to accept keywords arguments. For all other actors I have defined an outer constructor.

dd::IMAS.dd{T}
par::ParametersActor
function ActorNeutronics(dd::IMAS.dd, par::ParametersActor; kw...)
logging_actor_init(ActorNeutronics)
function ActorNeutronics{T}(dd::IMAS.dd{T}, par::ParametersActor; kw...) where {T}
logging_actor_init(ActorNeutronics{T})
par = par(kw...)
return new(dd, par)
return new{T}(dd, par)
end
end

ActorNeutronics(dd::IMAS.dd{T}, par::ParametersActor; kw...) where {T} = ActorNeutronics{T}(dd, par; kw...)

function ParametersActor(::Type{Val{:ActorNeutronics}})
par = ParametersActor(nothing)
par.N = Entry(Integer, "", "Number of particles"; default=100000)
Expand All @@ -53,10 +55,13 @@ function ActorNeutronics(dd::IMAS.dd, act::ParametersAllActors; kw...)
return actor
end

function _step(actor::ActorNeutronics; N::Integer=actor.par.N, step::Float64=actor.par.step, do_plot::Bool=actor.par.do_plot)
dd = actor.dd
cp1d = dd.core_profiles.profiles_1d[]
eqt = dd.equilibrium.time_slice[]
function define_neutrons(actor::ActorNeutronics; p=nothing)
N::Int = actor.par.N
step::Float64 = actor.par.step
do_plot::Bool = actor.par.do_plot

cp1d = actor.dd.core_profiles.profiles_1d[]
eqt = actor.dd.equilibrium.time_slice[]

# 2D neutron source
neutron_source_1d = IMAS.alpha_heating(cp1d) * 4 # W/m^3
Expand All @@ -76,67 +81,207 @@ function _step(actor::ActorNeutronics; N::Integer=actor.par.N, step::Float64=act
CDF = cumsum(neutron_source_2d[:])
W_per_trace = CDF[end] / N
CDF .= (CDF .- CDF[1]) ./ (CDF[end] - CDF[1])
draw = Int.(ceil.(IMAS.interp1d(CDF, 1:length(CDF), :linear).(rand(N))))
#draw = Int.(ceil.(IMAS.interp1d(CDF, 1:length(CDF), :linear).(rand(N))))
ICDF = IMAS.interp1d(CDF, 1:length(CDF), :linear)

# neutron structures
ϕ = rand(N) .* 2pi
θv = rand(N) .* 2pi
ϕv = acos.(rand(N) .* 2.0 .- 1.0)
neutrons = neutron_particle.(R[draw] .* cos.(ϕ), R[draw] .* sin.(ϕ), Z[draw], step * sin.(ϕv) .* cos.(θv), step * sin.(ϕv) .* sin.(θv), step * cos.(ϕv))
neutrons = Vector{neutron_particle{Float64}}(undef, N)
for k in eachindex(neutrons)
dk = Int(ceil(ICDF(rand())))
ϕ = rand() * 2pi
θv = rand() * 2pi
ϕv = acos(rand() * 2.0 - 1.0)
#dk = draw[k]
Rk = R[dk]
Zk = Z[dk]
yk, xk = Rk .* sincos(ϕ)
zk = Zk
δvyk, δvxk = (step * sin(ϕv)) .* sincos(θv)
δvzk = step * cos(ϕv)
neutrons[k] = neutron_particle(xk, yk, zk, δvxk, δvyk, δvzk)
end

# plot neutrons
if do_plot
histogram2d!(p,
Rcoord.(neutrons),
Zcoord.(neutrons),
nbins=(LinRange(minimum(r), maximum(r), length(r) - 1), LinRange(minimum(z), maximum(z), length(z) - 1)),
aspect_ratio=:equal,
weights=zeros(N) .+ 1 / 40,
)
end

return neutrons, W_per_trace
end

function define_wall(actor::ActorNeutronics)
# resample wall and make sure it's clockwise (for COCOS = 11)
wall = IMAS.first_wall(dd.wall)
wall_r, wall_z = IMAS.resample_2d_line(wall.r, wall.z; step=0.1)
eqt = actor.dd.equilibrium.time_slice[]
wall = IMAS.first_wall(actor.dd.wall)
rwall, zwall = IMAS.resample_2d_line(wall.r, wall.z; step=0.1)
R0 = eqt.global_quantities.magnetic_axis.r
Z0 = eqt.global_quantities.magnetic_axis.z
IMAS.reorder_flux_surface!(wall_r, wall_z, R0, Z0)
IMAS.reorder_flux_surface!(rwall, zwall, R0, Z0)

# advance neutrons until they hit the wall
rz_wall = collect(zip(wall_r, wall_z))
Threads.@threads for n in neutrons
while PolygonOps.inpolygon((Rcoord(n), Zcoord(n)), rz_wall) == 1
n.x += n.δvx
n.y += n.δvy
n.z += n.δvz
return rwall, zwall
end

function intersection(r1::Real, z1::Real, r2::Real, z2::Real, x::Real, y::Real, z::Real, vx::Real, vy::Real, vz::Real, v2::Nothing=nothing, vz2::Nothing=nothing)
return intersection(r1, z1, r2, z2, x, y, z, vx, vy, vz, vx^2 + vy^2, vz^2)
end

function intersection(r1::Real, z1::Real, r2::Real, z2::Real, x::Real, y::Real, z::Real, vx::Real, vy::Real, vz::Real, v2::Real, vz2::Real)
m = (z2 - z1)/(r2 - r1)
z0 = z1 - m * r1
m2 = m^2

a = 0.0
b = 0.0
c = 0.0
if isfinite(m)
z_z0 = z - z0
a = m2 * v2 - vz2
b = 2 * (m2 * (vx * x + vy * y) - z_z0 * vz)
c = m2 * (x^2 + y^2) - z_z0^2
else
a = v2
b = 2 * (vx * x + vy * y)
c = x^2 + y^2 - r1^2
end

t1 = -Inf
t2 = -Inf
if a == 0
b != 0 && (t2 = -c / b)
else
nb_2a = -b / (2a)
b2a_ca = nb_2a^2 - c / a
if b2a_ca >= 0
t_sq = sqrt(b2a_ca)
t1 = nb_2a - t_sq
t2 = nb_2a + t_sq
end
end

t = Inf
if t1 > 0
zi = vz * t1 + z
if zi > z1 && zi < z2
t = t1
end
end
if !isfinite(t) && t2 > 0
zi = vz * t2 + z
if zi > z1 && zi < z2
t = t2
end
end

return t
end

function _step(actor::ActorNeutronics)

do_plot::Bool = actor.par.do_plot
p = do_plot ? plot() : nothing

neutrons, W_per_trace = define_neutrons(actor; p)

rwall, zwall = define_wall(actor)
rz_wall = collect(zip(rwall, zwall))

Nw = length(rz_wall)

for n in neutrons
#println(n)
ti = Inf

xn = n.x
yn = n.y
zn = n.z
vx = n.δvx
vy = n.δvy
vz = n.δvz

v2 = vx^2 + vy^2
vz2 = vz^2

a = v2
b = 2 * (vx * xn + vy * yn)
c = xn^2 + yn^2
rmin = sqrt(c - b^2 / (4a))

@inbounds for k in eachindex(rz_wall)
t = Inf
rw1, zw1 = rz_wall[k]
rw2, zw2 = k == Nw ? rz_wall[1] : rz_wall[k+1]

rw1 == rw2 && zw1 == zw2 && continue

if zw1 > zw2
rw1, rw2 = rw2, rw1
zw1, zw2 = zw2, zw1
end

(vz > 0 && zn > zw2) && continue
(vz < 0 && zn < zw1) && continue
(rw1 < rmin && rw2 < rmin) && continue

t = intersection(rw1, zw1, rw2, zw2, xn, yn, zn, vx, vy, vz, v2, vz2)
t < ti && (ti = t)
end
n.x += vx * ti
n.y += vy * ti
n.z += vz * ti
end

# find neutron flux [counts/s/m²]
# smooth the load of each neutron withing a window
wall_r, wall_z = wall_r[1:end-1], wall_z[1:end-1]
d = sqrt.(IMAS.gradient(vcat(wall_r, wall_r[1])) .^ 2.0 .+ IMAS.gradient(vcat(wall_z, wall_z[1])) .^ 2.0)
d = (d[1:end-1] .+ d[2:end]) / 2.0
wall_r, wall_z = rwall[1:end-1], zwall[1:end-1]
rwall[end] = wall_r[1]
zwall[end] = wall_z[1]
d = sqrt.(IMAS.gradient(rwall) .^ 2.0 .+ IMAS.gradient(zwall) .^ 2.0)
@views d = (d[1:end-1] .+ d[2:end]) ./ 2.0
l = cumsum(d)
wall_s = d .* wall_r .* 2π
ns = 10
stencil = collect(-ns:ns)

nflux_r = zeros(size(wall_r))
nflux_z = zeros(size(wall_z))
ll = stencil .* 0.0
nflux_r = zero(wall_r)
nflux_z = zero(wall_z)
dwall = zero(wall_r)
index = zero(stencil)
ll = zeros(2ns + 1)
window = zeros(2ns + 1)
for n in neutrons
old_r = Rcoord(n)
old_z = Zcoord(n)
index0 = argmin((wall_r .- old_r) .^ 2.0 + (wall_z .- old_z) .^ 2.0)
index = mod.(stencil .+ index0 .- 1, length(wall_r)) .+ 1
@. dwall = (wall_r - old_r) ^ 2 + (wall_z - old_z) ^ 2
index0 = argmin(dwall)
index .= mod.(stencil .+ index0 .- 1, length(wall_r)) .+ 1

n.x += n.δvx
n.y += n.δvy
n.z += n.δvz
new_r = Rcoord(n)
new_z = Zcoord(n)

ll .= cumsum(d[index])
@views cumsum!(ll, d[index])
ll .-= ll[ns+1]
window = exp.(-(ll ./ (l[end] / length(l)) / (2 * ns + 1) * 5) .^ 2)
window = window ./ sum(window)
window .= exp.(.-(ll ./ (l[end] / length(l)) ./ (2ns + 1) .* 5) .^ 2)
window ./= sum(window)
unit_vector = sqrt((new_r - old_r)^2 + (new_z - old_z)^2)

nflux_r[index] += (new_r - old_r) ./ unit_vector .* window .* W_per_trace ./ wall_s[index]
nflux_z[index] += (new_z - old_z) ./ unit_vector .* window .* W_per_trace ./ wall_s[index]
for (k, i) in enumerate(index)
nflux_r[i] += (new_r - old_r) ./ unit_vector .* window[k] .* W_per_trace ./ wall_s[i]
nflux_z[i] += (new_z - old_z) ./ unit_vector .* window[k] .* W_per_trace ./ wall_s[i]
end
end

# IMAS assignements
# IMAS assignments
dd = actor.dd
dd.neutronics.first_wall.r = wall_r
dd.neutronics.first_wall.z = wall_z
ntt = resize!(dd.neutronics.time_slice)
Expand All @@ -146,18 +291,10 @@ function _step(actor::ActorNeutronics; N::Integer=actor.par.N, step::Float64=act

# plot neutron wall loading cx
if do_plot
neutrons = neutron_particle.(R[draw] .* cos.(ϕ), R[draw] .* sin.(ϕ), Z[draw], step * sin.(ϕv) .* cos.(θv), step * sin.(ϕv) .* sin.(θv), step * cos.(ϕv))

histogram2d(
Rcoord.(neutrons),
Zcoord.(neutrons),
nbins=(LinRange(minimum(r), maximum(r), length(r) - 1), LinRange(minimum(z), maximum(z), length(z) - 1)),
aspect_ratio=:equal,
weights=zeros(N) .+ 1 / 40,
)

plot!(ntt.wall_loading, xlim=[minimum(wall.r) * 0.9, maximum(wall.r) * 1.1], title="First wall neutron loading")
display(plot!(wall.r, wall.z, color=:black, label="", xlim=[minimum(wall.r) * 0.9, maximum(wall.r) * 1.1], title=""))
plot!(p, ntt.wall_loading, xlim=[minimum(rwall) * 0.9, maximum(rwall) * 1.1], title="First wall neutron loading")
plot!(p, rwall, zwall, color=:black, label="", xlim=[minimum(rwall) * 0.9, maximum(rwall) * 1.1], title="")
display(p)
end

return
end