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

Have fewer if conditions for straight guide propagation. #45

Open
wants to merge 3 commits into
base: main
Choose a base branch
from
Open

Have fewer if conditions for straight guide propagation. #45

Changes from 1 commit
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
1cf2c4e
Have fewer if conditions for straight guide propagation.
mtbc Feb 14, 2022
5499b7b
Merge remote-tracking branch 'origin/main' into reduce-ifs
mtbc Feb 23, 2022
ec9a95c
Remove another "if" in inner propagation loop.
mtbc Feb 23, 2022
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
80 changes: 31 additions & 49 deletions acc/components/optics/guide.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
#
# Copyright (c) 2021 by UT-Battelle, LLC.

from math import ceil, sqrt, tanh
from math import ceil, inf, sqrt, tanh
from numba import cuda
from time import time

Expand All @@ -22,11 +22,7 @@ def calc_reflectivity(Q, R0, Qc, alpha, m, W):
"""
R = R0
if Q > Qc:
tmp = (Q - m * Qc) / W
if tmp < 10:
R *= (1 - tanh(tmp)) * (1 - alpha * (Q - Qc)) / 2
else:
R = 0
R *= (1 - tanh((Q - m * Qc) / W)) * (1 - alpha * (Q - Qc)) / 2
return R


Expand Down Expand Up @@ -62,57 +58,43 @@ def propagate(
ch1 = -hh1*l - z*hh; ch2 = y*l
# Compute intersection times.
t1 = (l - z) / vz # for guide exit
i = 0
if vdotn_v1 < 0:
t2 = (cv1 - cv2)/vdotn_v1
if t2 < t1:
t1 = t2
i = 1
if vdotn_v2 < 0:
t2 = (cv1 + cv2)/vdotn_v2
if t2<t1:
t1 = t2
i = 2
if vdotn_h1 < 0:
t2 = (ch1 - ch2)/vdotn_h1
if t2<t1:
t1 = t2
i = 3
if vdotn_h2 < 0:
t2 = (ch1 + ch2)/vdotn_h2
if t2 < t1:
t1 = t2
i = 4
if i == 0:
break # Neutron left guide.
vdots = (-1., vdotn_v1, vdotn_v2, vdotn_h1, vdotn_h2)
times = (t1,
(cv1 - cv2) / vdotn_v1,
(cv1 + cv2) / vdotn_v2,
(ch1 - ch2) / vdotn_h1,
(ch1 + ch2) / vdotn_h2)
best_time = inf
for i in range(0, 5):
if vdots[i] < 0 and best_time > times[i]:
best_time = times[i]
side = i
Copy link
Collaborator Author

Choose a reason for hiding this comment

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

So, here, if I try,

is_better = vdots[i] < 0 and best_time > times[i]
if is_better:
    best_time = times[i]
    side = i

then it's still (relatively) fast but with,

is_better = vdots[i] < 0 and best_time > times[i]
best_time = is_better and times[i] or best_time
side = is_better and i or side

it becomes extremely slow. Have you any tips, @ckendrick?

Copy link
Collaborator

Choose a reason for hiding this comment

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

There might be an issue with how side is being set. I ran into an issue where I was setting the side to an incorrect value which caused neutrons to get stuck in the loop for all max_bounces number of iterations. If you set max_bounces = 100 and it is significantly faster, then this is probably the issue. Also, running test_guide.py in interactive mode should show differences if the side is incorrect.

Sorry I haven't opened a PR yet for the optimized tapered guide. The equivalent loop from mine is

i = 0
for k in range(4):
    side_hit = vdots[k] < 0.0
    t2 = times[k] / vdots[k]
    best_time = side_hit and t2 < t1
    i += best_time * int(k + 1)
    t1 = ((not best_time) * t1) + best_time * t2
if i == 0:
    break

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Aha, thank you. Nice trick with using the numerical aspect of the truth value more directly than I.

And, yes, reducing max_bounces correspondingly speeds the run. Maybe my two alternatives aren't logically equivalent and I've yet to spot why!

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I ran into an issue where I was setting the side to an incorrect value ...

and so was I,

... my two alternatives aren't logically equivalent ...

so they weren't. 😃

Will push a followup commit but still doesn't improve performance significantly.

if side == 0:
# Neutron left guide.
break
t1 = best_time

# propagate time t1 to move to reflection point
x += vx*t1; y += vy*t1; z += vz*t1; t += t1

# reflection
if i == 1: # Left vertical mirror
nlen2 = l*l + ww*ww
q = V2K*(-2)*vdotn_v1/sqrt(nlen2)
d = 2*vdotn_v1/nlen2
vx = vx - d*l
vz = vz - d*ww
elif i == 2: # Right vertical mirror
if side < 3:
# left or right vertical mirror
vdot = (vdotn_v1, vdotn_v2)[side - 1]
nlen2 = l*l + ww*ww
q = V2K*(-2)*vdotn_v2/sqrt(nlen2)
d = 2*vdotn_v2/nlen2
vx = vx + d*l
q = V2K*(-2)*vdot/sqrt(nlen2)
d = 2*vdot/nlen2
vxd = (d*l, -d*l)[side - 1]
vx = vx - vxd
vz = vz - d*ww
elif i == 3: # Lower horizontal mirror
nlen2 = l*l + hh*hh
q = V2K*(-2)*vdotn_h1/sqrt(nlen2)
d = 2*vdotn_h1/nlen2
vy = vy - d*l
vz = vz - d*hh
elif i == 4: # Upper horizontal mirror
else:
# lower or upper horizontal mirror
vdot = (vdotn_h1, vdotn_h2)[side - 3]
nlen2 = l*l + hh*hh
q = V2K*(-2)*vdotn_h2/sqrt(nlen2)
d = 2*vdotn_h2/nlen2
vy = vy + d*l
q = V2K*(-2)*vdot/sqrt(nlen2)
d = 2*vdot/nlen2
vyd = (d*l, -d*l)[side - 3]
vy = vy - vyd
vz = vz - d*hh
R = calc_reflectivity(q, R0, Qc, alpha, m, W)
prob *= R
Expand Down