diff --git a/acc/components/__init__.py b/acc/components/__init__.py new file mode 100644 index 0000000..3b4d236 --- /dev/null +++ b/acc/components/__init__.py @@ -0,0 +1 @@ +from . import monitors, optics, samples, sources \ No newline at end of file diff --git a/acc/components/monitors/__init__.py b/acc/components/monitors/__init__.py new file mode 100644 index 0000000..f9b6449 --- /dev/null +++ b/acc/components/monitors/__init__.py @@ -0,0 +1,5 @@ +from .iqe_monitor import IQE_monitor +from .posdiv_monitor import PosDiv_monitor +from .psd_monitor import PSD_monitor +from .psd_monitor_4pi import PSD_monitor_4Pi +from .wavelength_monitor import Wavelength_monitor \ No newline at end of file diff --git a/acc/components/optics/__init__.py b/acc/components/optics/__init__.py new file mode 100644 index 0000000..046151d --- /dev/null +++ b/acc/components/optics/__init__.py @@ -0,0 +1,7 @@ +from .arm import Arm +from .beamstop import Beamstop +from .guide import Guide +# from .guide_anyshape_gravity import Guide_anyshape_gravity +from .guide_anyshape import Guide_anyshape +from .guide_tapering import Guide as Guide_tapering +from .slit import Slit diff --git a/acc/components/optics/geometry2d.py b/acc/components/optics/geometry2d.py index d2b80fc..89484b1 100644 --- a/acc/components/optics/geometry2d.py +++ b/acc/components/optics/geometry2d.py @@ -11,7 +11,7 @@ # modified from # https://stackoverflow.com/questions/1119627/how-to-test-if-a-point-is-inside-of-a-convex-polygon-in-2d-integer-coordinates#:~:text=If%20it%20is%20convex%2C%20a,traversed%20in%20the%20same%20order) @cuda.jit(device=True) -def inside_convex_polygon(point, vertices): +def inside_convex_polygon(point, vertices, l_epsilon): # previous_side = get_side(v_sub(vertices[1], vertices[0]), v_sub(point, vertices[0])) previous_sign = cosine_sign(v_sub(vertices[1], vertices[0]), v_sub(point, vertices[0])) previous_side = previous_sign <= 0 @@ -21,11 +21,11 @@ def inside_convex_polygon(point, vertices): affine_segment = v_sub(b, a) affine_point = v_sub(point, a) sign = cosine_sign(affine_segment, affine_point) - if fabs(sign) < 1E-14: + if fabs(sign) < l_epsilon: continue # current_side = get_side(affine_segment, affine_point) current_side = sign <= 0 - if fabs(previous_sign) < 1E-14: + if fabs(previous_sign) < l_epsilon: previous_side = current_side continue if previous_side != current_side: diff --git a/acc/components/optics/guide_anyshape.py b/acc/components/optics/guide_anyshape.py index d7a00a7..8318d2b 100644 --- a/acc/components/optics/guide_anyshape.py +++ b/acc/components/optics/guide_anyshape.py @@ -1,6 +1,8 @@ # -*- python -*- # +import logging +logger = logging.getLogger(__name__) from math import sqrt, fabs import numpy as np, numba from numba import cuda, void @@ -11,14 +13,16 @@ from ...config import get_numba_floattype NB_FLOAT = get_numba_floattype() from ._guide_utils import calc_reflectivity -from ...neutron import absorb +from ...neutron import absorb, clone, prop_dt_inplace from ...geometry._utils import insert_into_sorted_list_with_indexes from .geometry2d import inside_convex_polygon from ... import vec3 -max_bounces = 1000 -max_numfaces = 100 - +max_bounces = 30 +max_numfaces = 5000 +history_size_for_bounces_in_epsilon_neighborhood = 10 # This normally would not exceed 3 for typical guide design. +t_epsilon = 1E-14 +l_epsilon = 1E-11 @cuda.jit(NB_FLOAT(NB_FLOAT[:], NB_FLOAT[:], NB_FLOAT[:], NB_FLOAT[:, :], NB_FLOAT[:]), device=True, inline=True) @@ -85,60 +89,134 @@ def load_scaled_centered_faces(path, xwidth=0, yheight=0, zdepth=0, center=False faces = np.array([[vertices[i] for i in face] for face in faces]) return faces +@cuda.jit(numba.boolean(NB_FLOAT[:], NB_FLOAT[:], NB_FLOAT[:, :], NB_FLOAT[:, :]), + device=True, inline=True) +def likely_inside_face(postmp, face_center, face_uvecs, face2d_bounds): + vec3.subtract(postmp, face_center, postmp) + ex = face_uvecs[0] + ey = face_uvecs[1] + x = vec3.dot(ex, postmp) + y = vec3.dot(ey, postmp) + return ( + (x>face2d_bounds[0, 0]-l_epsilon) and (xface2d_bounds[0, 1]-l_epsilon) and (y0: - ninter = insert_into_sorted_list_with_indexes(iface, intersection, face_indexes, intersections, ninter) + found_in_history = False + # logger.debug(f" face {iface}, {face_hist_start_ind}, {face_hist_start_ind+face_hist_size}") + for ifh in range(face_hist_start_ind, face_hist_start_ind+face_hist_size): + # logger.debug(f" face {iface}, {ifh}, {tmp_face_hist[ifh]}") + if iface == tmp_face_hist[ifh]: + found_in_history = True + break + # logger.debug(f" face {iface} of {nfaces} faces: found_in_history={found_in_history}") + if found_in_history: continue + x0 = in_neutron[:3]; v0 = in_neutron[3:6] + face_center = centers[iface] + face_uvecs = unitvecs[iface] + intersection = intersect_plane( x0, v0, face_center, face_uvecs, tmpv3) + # logger.debug(f" face {iface}: x0={x0}, v0={v0}, face_center={face_center}, face_uvecs={face_uvecs}, intersection={intersection}") + # if intersection<=-t_epsilon: + # if intersection<0: + if intersection<=-t_epsilon/5: + continue + face2d_bounds = bounds2d[iface] + clone(in_neutron, tmp_neutron) + prop_dt_inplace(tmp_neutron, intersection) + vec3.copy(tmp_neutron[:3], tmpv3) + if not likely_inside_face(tmpv3, face_center, face_uvecs, face2d_bounds): + continue + # logger.debug(f" face {iface}: intersection likely inside face") + # vec3.copy(tmp_neutron[:3], tmpv3) + # if intersection < t_epsilon and on_the_other_side(tmpv3, face_center, face_uvecs[2]): + # continue + # logger.debug(f" face {iface}: intersection on the right side") + ninter = insert_into_sorted_list_with_indexes(iface, intersection, face_indexes, intersections, ninter) + # logger.debug(f" face {iface}: new intersection list {intersections[:ninter]} for faces {face_indexes[:ninter]}") + # logger.debug(f"bounce {nb}: ninter={ninter}") + # import pdb; pdb.set_trace() if not ninter: break # find the smallest intersection that is inside the mirror - found = False + found = -1 for iinter in range(ninter): face_index = face_indexes[iinter] intersection = intersections[iinter] # calc position of intersection - vec3.copy(in_neutron[3:6], tmp1) - vec3.scale(tmp1, intersection) - vec3.add(tmp1, in_neutron[:3], tmp1) + vec3.copy(in_neutron[3:6], tmpv3) + vec3.scale(tmpv3, intersection) + vec3.add(tmpv3, in_neutron[:3], tmpv3) # calc 2d coordinates and use it to check if it is inside the mirror - vec3.subtract(tmp1, centers[face_index], tmp1) + vec3.subtract(tmpv3, centers[face_index], tmpv3) e0 = unitvecs[face_index, 0, :] e1 = unitvecs[face_index, 1, :] e2 = unitvecs[face_index, 2, :] face2d = faces2d[face_index] - if inside_convex_polygon((vec3.dot(tmp1, e0), vec3.dot(tmp1, e1)), face2d): - found = True + # logger.debug(f"check intersection {iinter}: face={face_center}, {face_uvecs[2]}, intersection={intersection}") + if inside_convex_polygon((vec3.dot(tmpv3, e0), vec3.dot(tmpv3, e1)), face2d, l_epsilon): + found = face_index + # logger.debug(f"check intersection {iinter}: found={found}") + to_travel_after_bounce = l_epsilon/vec3.length(in_neutron[3:6]) # move a tiny distance + if iinter0: - ninter = insert_into_sorted_list_with_indexes(iface, tmp1[0], face_indexes, intersections, ninter) - if tmp1[1]>0: - ninter = insert_into_sorted_list_with_indexes(iface, tmp1[1], face_indexes, intersections, ninter) + if tmp1[0]>-t_epsilon: + propagate_with_gravity(in_neutron[:3], in_neutron[3:6], gravity, tmp1[0], tmp_neutron[:3], tmp_neutron[3:6], tmp2) + if likely_inside_face(tmp_neutron[:3], face_center, face_uvecs, face2d_bounds): + ninter = insert_into_sorted_list_with_indexes(iface, tmp1[0], face_indexes, intersections, ninter) + if tmp1[1]>-t_epsilon: + propagate_with_gravity(in_neutron[:3], in_neutron[3:6], gravity, tmp1[1], tmp_neutron[:3], tmp_neutron[3:6], tmp2) + if likely_inside_face(tmp_neutron[:3], face_center, face_uvecs, face2d_bounds): + ninter = insert_into_sorted_list_with_indexes(iface, tmp1[1], face_indexes, intersections, ninter) if not ninter: break # find the smallest intersection that is inside the mirror - found = False + found = -1 for iinter in range(ninter): face_index = face_indexes[iinter] intersection = intersections[iinter] # calc position and velocity at intersection propagate_with_gravity(in_neutron[:3], in_neutron[3:6], gravity, intersection, tmp1, tmp2, tmp3) e2 = unitvecs[face_index, 2, :] - # The next bounce should not be very close in time. - # This handles the case of numeric error where the last intersection - # calculated is a bit larger than exact value and the neutron - # was propagated to just slightly behind the mirror. - # However, this will also exclude the case where neutron hit right next to the corner. - # There could be leakages for those cases. - # Also this assumes that initially (before this guide component) - # the neutron is not right next to a mirror. - if intersection < 1E-10: # and vec3.dot(tmp2, e2)<0: - continue # calc 2d coordinates and use it to check if it is inside the mirror vec3.subtract(tmp1, centers[face_index], tmp3) e0 = unitvecs[face_index, 0, :] e1 = unitvecs[face_index, 1, :] face2d = faces2d[face_index] - if inside_convex_polygon((vec3.dot(tmp3, e0), vec3.dot(tmp3, e1)), face2d): - found = True + if inside_convex_polygon((vec3.dot(tmp3, e0), vec3.dot(tmp3, e1)), face2d, l_epsilon): + found = face_index + # logger.debug(f"check intersection {iinter}: found={found}") + to_travel_after_bounce = l_epsilon/vec3.length(in_neutron[3:6]) # move a tiny distance + if iinter