From 1cf2c4e363c736f391e1d54bd13232c1b54e5b8e Mon Sep 17 00:00:00 2001 From: Mark Carroll Date: Mon, 14 Feb 2022 11:54:23 -0500 Subject: [PATCH 1/2] Have fewer if conditions for straight guide propagation. --- acc/components/optics/guide.py | 80 +++++++++++++--------------------- 1 file changed, 31 insertions(+), 49 deletions(-) diff --git a/acc/components/optics/guide.py b/acc/components/optics/guide.py index b6bd2bb8..0e2328bf 100644 --- a/acc/components/optics/guide.py +++ b/acc/components/optics/guide.py @@ -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 @@ -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 @@ -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 times[i]: + best_time = times[i] + side = i + 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 From ec9a95c09867ae050a29916d619b6551daec0555 Mon Sep 17 00:00:00 2001 From: Mark Carroll Date: Wed, 23 Feb 2022 13:50:37 -0500 Subject: [PATCH 2/2] Remove another "if" in inner propagation loop. --- acc/components/optics/guide.py | 36 +++++++++++++++++----------------- 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/acc/components/optics/guide.py b/acc/components/optics/guide.py index a8f61b65..014eff77 100644 --- a/acc/components/optics/guide.py +++ b/acc/components/optics/guide.py @@ -5,7 +5,7 @@ import numpy as np from math import ceil, inf, sqrt, tanh -from numba import cuda, float32, void +from numba import cuda, float32, float64, int64, void from time import time from mcni.AbstractComponent import AbstractComponent @@ -14,7 +14,6 @@ category = 'optics' - @cuda.jit(float32(float32, float32, float32, float32, float32, float32), device=True, inline=True) def calc_reflectivity(Q, R0, Qc, alpha, m, W): @@ -65,18 +64,19 @@ def propagate( ch1 = -hh1*l - z*hh; ch2 = y*l # Compute intersection times. t1 = (l - z) / vz # for guide exit - 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) + vdots = (float64(0), float64(-1), + vdotn_v1, vdotn_v2, + vdotn_h1, vdotn_h2) + times = (float64(0), 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 - if side == 0: + side = -1 + for i in range(1, 6): + is_better = vdots[i] < 0 and best_time > times[i] + best_time = float64(is_better) and times[i] or best_time + side = int64(is_better) and i or side + if side == 1: # Neutron left guide. break t1 = best_time @@ -85,22 +85,22 @@ def propagate( x += vx*t1; y += vy*t1; z += vz*t1; t += t1 # reflection - if side < 3: + if side < 4: # left or right vertical mirror - vdot = (vdotn_v1, vdotn_v2)[side - 1] + vdot = (vdotn_v1, vdotn_v2)[side - 2] nlen2 = l*l + ww*ww q = V2K*(-2)*vdot/sqrt(nlen2) d = 2*vdot/nlen2 - vxd = (d*l, -d*l)[side - 1] + vxd = (d*l, -d*l)[side - 2] vx = vx - vxd vz = vz - d*ww else: # lower or upper horizontal mirror - vdot = (vdotn_h1, vdotn_h2)[side - 3] + vdot = (vdotn_h1, vdotn_h2)[side - 4] nlen2 = l*l + hh*hh q = V2K*(-2)*vdot/sqrt(nlen2) d = 2*vdot/nlen2 - vyd = (d*l, -d*l)[side - 3] + vyd = (d*l, -d*l)[side - 4] vy = vy - vyd vz = vz - d*hh R = calc_reflectivity(q, R0, Qc, alpha, m, W)