From e9d9aa97a699aace7c3faf6302fe69877759fe4c Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Thu, 27 Apr 2023 22:44:27 -0700 Subject: [PATCH 01/53] added method `forward_distance_in_shape`: acc/geometry/propagation.py --- acc/geometry/propagation.py | 37 ++++++++++++++++++++++++++++++ acc/vec3.py | 14 +++++++++++ tests/geometry/test_propagation.py | 26 ++++++++++++++++++++- 3 files changed, 76 insertions(+), 1 deletion(-) diff --git a/acc/geometry/propagation.py b/acc/geometry/propagation.py index 3df424a4..5de417d5 100644 --- a/acc/geometry/propagation.py +++ b/acc/geometry/propagation.py @@ -4,6 +4,7 @@ from .arrow_intersect import max_intersections from .arrow_intersect import inside, outside, onborder from ..neutron import prop_dt_inplace +from ..vec3 import distance as v3_dist, length as v3_length from .. import test def makePropagateMethods(intersect, locate): @@ -190,6 +191,33 @@ def _tof_before_first_exit(neutron, ts): continue return t + def _forward_distance_in_shape(neutron, end, ts): + "distance to travel in shape to get to `end`. `end` must be on the path of the neutron" + start = neutron[:3] + vv = neutron[3:6] + v = v3_length(vv) + nintersect = forward_intersect(start[0], start[1], start[2], vv[0], vv[1], vv[2], ts) + length = v3_dist(start, end) + tofmax = length/v; + prev = 0; length = 0 + for i in range(nintersect): + tof = ts[i] + if (tof > tofmax): + tof = tofmax + middle = (tof+prev)/2. + # middle point p = start + vv*middle; + x = start[0] + vv[0]*middle + y = start[1] + vv[1]*middle + z = start[2] + vv[2]*middle + # if middle point is inside, count this segment + if locate(x,y,z) == inside: + length += (tof-prev) * v; + if tof > tofmax: + break + prev = tof + continue + return length + if test.USE_CUDASIM: @cuda.jit(device=True, inline=True) @@ -212,6 +240,10 @@ def propagate_to_next_incident_surface(neutron): def propagate_to_next_exiting_surface(neutron): ts = np.zeros(max_intersections, dtype=float) return _propagate_to_next_exiting_surface(neutron, ts) + @cuda.jit(device=True, inline=True) + def forward_distance_in_shape(neutron, end): + ts = np.zeros(max_intersections, dtype=float) + return _forward_distance_in_shape(neutron, end, ts) else: @cuda.jit(device=True, inline=True) def is_exiting(neutron): @@ -233,6 +265,10 @@ def propagate_to_next_incident_surface(neutron): def propagate_to_next_exiting_surface(neutron): ts = cuda.local.array(max_intersections, dtype=numba.float64) return _propagate_to_next_exiting_surface(neutron, ts) + @cuda.jit(device=True, inline=True) + def forward_distance_in_shape(neutron, end): + ts = cuda.local.array(max_intersections, dtype=numba.float64) + return _forward_distance_in_shape(neutron, end, ts) return dict( is_exiting = is_exiting, @@ -241,4 +277,5 @@ def propagate_to_next_exiting_surface(neutron): propagate_out = propagate_out, propagate_to_next_incident_surface = propagate_to_next_incident_surface, propagate_to_next_exiting_surface = propagate_to_next_exiting_surface, + forward_distance_in_shape = forward_distance_in_shape, ) diff --git a/acc/vec3.py b/acc/vec3.py index 42068728..6b35acd4 100644 --- a/acc/vec3.py +++ b/acc/vec3.py @@ -19,6 +19,13 @@ def cross(v1, v2, vout): def length(v): return math.sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]) +@cuda.jit(device=True, inline=True) +def distance(v1, v2): + x = v1[0]-v2[0] + y = v1[1]-v2[1] + z = v1[2]-v2[2] + return math.sqrt(x*x+y*y+z*z) + @cuda.jit(device=True, inline=True) def normalize(v): l = length(v) @@ -39,6 +46,13 @@ def copy(v1, v2): v2[2]=v1[2] return +@cuda.jit(device=True, inline=True) +def add_inplace(v1, v2): # v1+=v2 + v1[0]+=v2[0] + v1[1]+=v2[1] + v1[2]+=v2[2] + return + @cuda.jit(device=True, inline=True) def add(v1, v2, v3): v3[0]=v1[0]+v2[0] diff --git a/tests/geometry/test_propagation.py b/tests/geometry/test_propagation.py index 000582dc..ccc6cc7a 100755 --- a/tests/geometry/test_propagation.py +++ b/tests/geometry/test_propagation.py @@ -134,7 +134,7 @@ def test_propagate_to_next_exiting_surface(): [-0.01,0.,0., 1.,0.,0., 0.,0., 0.015, 1.], [-0.01,0.,0., 1.,0.,0., 0.,0., 0.01, 1.], [-0.01,0.,0., 1.,0.,0., 0.,0., 0.005, 1.], - [-0.01,0.,0., 1.,0.,0., 0.,0., 0.01, 1.], + [-0.01,0.,0., 1.,0.,0., 0.,0., 0.0, 1.], [0.02,0.,0., 1.,0.,0., 0.,0., 0.02, 1.], [0.02,0.,0., 1.,0.,0., 0.,0., 0.01, 1.], [0.02,0.,0., 1.,0.,0., 0.,0., 0.005, 1.], @@ -142,6 +142,29 @@ def test_propagate_to_next_exiting_surface(): [0.025,0.,0., 1.,0.,0., 0.,0., 0., 1.], [1.,0.,0., 1.,0.,0., 0.,0., 0., 1.], ]) + for i, (neutron, out) in enumerate(zip(neutrons, expected)): + propagate_to_next_exiting_surface(neutron) + assert np.allclose(neutron, out), f"{i}: {neutron} != {out}" + return + +@pytest.mark.skipif(not test.USE_CUDASIM, reason='no CUDASIM') +def test_forward_distance_in_shape(): + forward_distance_in_shape = getPropagateMethods()[ + 'forward_distance_in_shape'] + neutrons = np.array([ + [-1.,0.,0., 1.,0.,0., 0.,0., 0., 1.], + ]) + ends = np.array([ + [0.,0.,0.], + [1.,0.,0.], + ]) + expected = np.array([ + 0.01, + 0.02, + ]) + for neutron, end, expected_out in zip(neutrons, ends, expected): + out = forward_distance_in_shape(neutron, end) + assert np.allclose(out, expected_out), f"{out} != {expected_out}" return def main(): @@ -149,6 +172,7 @@ def main(): test_tof_before_first_exit() test_propagate_to_next_incident_surface() test_propagate_to_next_exiting_surface() + test_forward_distance_in_shape() return if __name__ == '__main__': main() From e478f4dcc04b91c616e7b305dc55585b970d3902 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Fri, 28 Apr 2023 19:10:10 -0700 Subject: [PATCH 02/53] added acc.scatterers.homogeneous_scatterer module: currently only implemented `interact_path1`; added a test --- acc/scatterers/__init__.py | 29 +++++ acc/scatterers/homogeneous_scatterer.py | 120 ++++++++++++++++++ acc/scatterers/interaction_types.py | 3 + tests/scatterers/sampleassemblies | 1 + .../scatterers/test_homogeneous_scatterer.py | 33 +++++ 5 files changed, 186 insertions(+) create mode 100644 acc/scatterers/__init__.py create mode 100644 acc/scatterers/homogeneous_scatterer.py create mode 100644 acc/scatterers/interaction_types.py create mode 120000 tests/scatterers/sampleassemblies create mode 100644 tests/scatterers/test_homogeneous_scatterer.py diff --git a/acc/scatterers/__init__.py b/acc/scatterers/__init__.py new file mode 100644 index 00000000..3e30290c --- /dev/null +++ b/acc/scatterers/__init__.py @@ -0,0 +1,29 @@ +import numba +from numba import cuda +from mcni import units + +class ScatterFuncFactory: + + def render(self, scatterer): + """returns cuda device function tuple (scatter, interact_path1) + """ + scatter, interact_path1 = scatterer.identify(self) + return scatter, interact_path1 + + def onCompositeKernel(self, composite): + elements = composite.elements() + weights = [] + for element in elements: + weights.append(getattr(element, 'weight', 1.)) + continue + s = sum(weights) + for w, element in zip(weights, elements): + element.weight = w/s + from .composite import makeKernelMethods + return makeKernelMethods(composite) + + def onHomogeneousScatterer(self, hs): + return + + +scatter_func_factory = ScatterFuncFactory() diff --git a/acc/scatterers/homogeneous_scatterer.py b/acc/scatterers/homogeneous_scatterer.py new file mode 100644 index 00000000..d4cd7e79 --- /dev/null +++ b/acc/scatterers/homogeneous_scatterer.py @@ -0,0 +1,120 @@ +# -*- python -*- +# +# Jiao Lin +# + +import numpy as np, numba +from numba import cuda, void, int64 +from mcvine.acc._numba import xoroshiro128p_uniform_float32 +from math import sqrt, exp + +from .interaction_types import absorption, scattering, none +from .. import test +from ..vec3 import ( + subtract as v3_subtract, + length as v3_length, + add as v3_add, + add_inplace as v3_add_inplace, + scale as v3_scale, + copy as v3_copy, +) +from ..neutron import absorb, prop_dt_inplace +from ..geometry.arrow_intersect import max_intersections +from ..geometry.locate import inside, outside, onborder +from ..geometry.propagation import makePropagateMethods + + +from numba.core import config +if not config.ENABLE_CUDASIM: + from numba.cuda.compiler import Dispatcher, DeviceFunction + +from ..config import get_numba_floattype +NB_FLOAT = get_numba_floattype() + + +category = 'samples' + +def factory(shape, kernel, mcweights, packing_factor): + w_absorption, w_scattering, w_transmission = mcweights + from ..geometry import arrow_intersect + intersect = arrow_intersect.arrow_intersect_func_factory.render(shape) + locate = arrow_intersect.locate_func_factory.render(shape) + propagate_methods = makePropagateMethods(intersect, locate) + propagate_to_next_incident_surface = propagate_methods['propagate_to_next_incident_surface'] + propagate_to_next_exiting_surface = propagate_methods['propagate_to_next_exiting_surface'] + tof_before_first_exit = propagate_methods['tof_before_first_exit'] + forward_distance_in_shape = propagate_methods["forward_distance_in_shape"] + from ..kernels import scatter_func_factory + scatter, calc_scattering_coeff, absorb, calc_absorption_coeff = scatter_func_factory.render(kernel) + def calculate_attenuation(neutron, end): + length = forward_distance_in_shape(neutron, end) + sigma = calc_scattering_coeff(neutron) * packing_factor + mu = calc_absorption_coeff(neutron) * packing_factor + return exp( - (mu+sigma) * length ) + + def _interact_path1(threadindex, rng_states, neutron, tmp_neutron): + x, y, z, vx, vy, vz = neutron[:6] + loc = locate(x,y,z) + if loc != inside: + # propagate to the front surface + propagate_to_next_incident_surface(neutron) + # distance to travel across and leave for the first time + tof = tof_before_first_exit(neutron) + velocity = sqrt(vx*vx+vy*vy+vz*vz) + distance = velocity * tof + # scattering and absorption coeffs + sigma = calc_scattering_coeff(neutron) + mu = calc_absorption_coeff(neutron) + # probability of three interaction types happening + transmission_prob = exp( -(mu+sigma)*distance ) + absorption_prob = (1-transmission_prob)*(mu/(mu+sigma)) + # scattering_prob = (1-transmission_prob)*(sigma/(mu+sigma)) + # + # toss a dice and decide whether we should do transmission, absorption, + # or scattering + transmission_mark = w_transmission + absorption_mark = transmission_mark + w_absorption + sum_of_weights = absorption_mark + w_scattering + r = xoroshiro128p_uniform_float32(rng_states, threadindex) * sum_of_weights + if (r < transmission_mark) : + # transmission + propagate_to_next_exiting_surface( neutron ) + neutron[-1] *= transmission_prob * (sum_of_weights/w_transmission) + return none + if (r >= transmission_mark and r < absorption_mark ): + # absorption + x = xoroshiro128p_uniform_float32(rng_states, threadindex) * distance + prob = mu * distance * exp( -(mu+sigma) * x ) + neutron[-1] *= prob * (sum_of_weights/w_absorption) + prop_dt_inplace( neutron, x/velocity ) + absorb( threadindex, rng_states, neutron ) + neutron[-1] = -1 + return absorption + # scattering + x = xoroshiro128p_uniform_float32(rng_states, threadindex) * distance + atten = exp( -(mu+sigma) * x ) + prob = distance * atten + prob *= sum_of_weights/w_scattering + neutron[-1] *= prob + prop_dt_inplace( neutron, x/velocity ) + scatter( threadindex, rng_states, neutron ) + neutron[-1] *= packing_factor + if neutron[-1] <=0: + return absorption + for i in range(10): + tmp_neutron[i] = neutron[i] + propagate_to_next_exiting_surface( neutron ) + atten2 = calculate_attenuation( tmp_neutron, neutron[:3] ); + neutron[-1] *= atten2; + return scattering + if test.USE_CUDASIM: + def interact_path1(threadindex, rng_states, neutron): + tmp_neutron = np.zeros(10, dtype=float) + _interact_path1(threadindex, rng_states, neutron, tmp_neutron) + else: + def interact_path1(threadindex, rng_states, neutron): + tmp_neutron = cuda.local.array(10, dtype=numba.float64) + _interact_path1(threadindex, rng_states, neutron, tmp_neutron) + return dict( + interact_path1 = interact_path1 + ) diff --git a/acc/scatterers/interaction_types.py b/acc/scatterers/interaction_types.py new file mode 100644 index 00000000..0ccc30f9 --- /dev/null +++ b/acc/scatterers/interaction_types.py @@ -0,0 +1,3 @@ +absorption = 0 +scattering = 1 +none = 2 diff --git a/tests/scatterers/sampleassemblies b/tests/scatterers/sampleassemblies new file mode 120000 index 00000000..e466c9c6 --- /dev/null +++ b/tests/scatterers/sampleassemblies @@ -0,0 +1 @@ +../components/samples/sampleassemblies \ No newline at end of file diff --git a/tests/scatterers/test_homogeneous_scatterer.py b/tests/scatterers/test_homogeneous_scatterer.py new file mode 100644 index 00000000..5f12ea2e --- /dev/null +++ b/tests/scatterers/test_homogeneous_scatterer.py @@ -0,0 +1,33 @@ +import os +thisdir = os.path.dirname(__file__) +import numpy as np + +def test(): + path = os.path.join(thisdir, "sampleassemblies", 'isotropic_sphere', 'sampleassembly.xml') + from mcvine.acc.components.samples import loadFirstHomogeneousScatterer + hs = loadFirstHomogeneousScatterer(path) + shape = hs.shape() + kernel = hs.kernel() + mcweights = 1., 1., 1. + packing_factor = 0.6 + from mcvine.acc.scatterers.homogeneous_scatterer import factory + methods = factory(shape, kernel, mcweights, packing_factor) + interact_path1 = methods['interact_path1'] + threadindex = 0 + rng_states = None + N = 10000 + p = 0. + absorbed = 0 + for i in range(N): + neutron = np.array([0.0,0,0, 0,0,1000, 0,0, 0, 1.]) + interact_path1(threadindex, rng_states, neutron) + p1 = neutron[-1] + if p1>0: + p+=p1 + else: + absorbed += 1 + p/=N + absorbed/=N + assert p<.8 and p>.7 + assert absorbed < 0.36 and absorbed > 0.3 + return From 5e97e551df02b0bdfba424729bf8a54270742aae Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Fri, 28 Apr 2023 19:11:16 -0700 Subject: [PATCH 03/53] added acc.scatterers subpkg: setup.py --- setup.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/setup.py b/setup.py index db57262e..a3dd8ae0 100644 --- a/setup.py +++ b/setup.py @@ -22,6 +22,7 @@ 'mcvine.acc.components.optics', 'mcvine.acc.components.monitors', 'mcvine.acc.components.samples', + 'mcvine.acc.scatterers', ], package_dir = { 'mcvine.acc': "acc", @@ -33,6 +34,7 @@ 'mcvine.acc.components.optics' : "acc/components/optics", 'mcvine.acc.components.monitors' : "acc/components/monitors", 'mcvine.acc.components.samples' : "acc/components/samples", + 'mcvine.acc.scatterers' : "acc/scatterers", }, data_files = [], install_requires = [ From 38a6659df580bb864d3a62c613348ac9c1303ebe Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Fri, 28 Apr 2023 19:11:39 -0700 Subject: [PATCH 04/53] ignore test outputs and build dirs: .gitignore --- .gitignore | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index c9f1dbf5..bfca6d3f 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,8 @@ __pycache__/ .ipynb_checkpoints/ -*.~undo-tree~ \ No newline at end of file +*.~undo-tree~ +out.* +log.* +mcvine.acc.egg-info/ +build/ +dist/ \ No newline at end of file From 4f82ed2dc1d67272b1259fff5e6bb452368f95eb Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Sun, 30 Apr 2023 20:18:36 -0700 Subject: [PATCH 05/53] refactored using new object _numba.coder: acc/kernels/composite.py --- acc/_numba.py | 21 +++++ acc/kernels/composite.py | 184 ++++++++++++++++++++++++++++++--------- 2 files changed, 164 insertions(+), 41 deletions(-) diff --git a/acc/_numba.py b/acc/_numba.py index f9313951..1bf7ab45 100644 --- a/acc/_numba.py +++ b/acc/_numba.py @@ -1,3 +1,4 @@ +import os, tempfile import numpy as np from numba.core import config if config.ENABLE_CUDASIM: @@ -6,3 +7,23 @@ def xoroshiro128p_uniform_float32(rng_states, threadindex): from numba.cuda.random import xoroshiro128p_type else: from numba.cuda.random import xoroshiro128p_uniform_float32, xoroshiro128p_type + +class Coder: + + def __init__(self, workdir=None): + self.workdir = workdir = workdir or os.path.abspath(".mcvine.acc.coder") + if not os.path.exists(workdir): + os.makedirs(workdir) + + def createUniqueDir(self, prefix): + return tempfile.mkdtemp(prefix=prefix, dir=self.workdir) + + @classmethod + def unrollLoop(cls, N, indent=4*' ', before_loop=None, in_loop=None, after_loop=None): + lines=[indent+line for line in before_loop or []] + for i in range(N): + lines+=[indent+line.format(i=i) for line in in_loop or []] + lines+=[indent+line for line in after_loop or []] + return lines + +coder = Coder() diff --git a/acc/kernels/composite.py b/acc/kernels/composite.py index 654d555f..d6385ded 100644 --- a/acc/kernels/composite.py +++ b/acc/kernels/composite.py @@ -14,62 +14,109 @@ def makeKernelMethods(composite): _modules = {} def makeKernelModule(composite): "make cuda device methods for the composite kernel" + from .._numba import coder if composite in _modules: return _modules[composite] import pickle as pkl, tempfile # save composite to be loaded by the "compiled" module - tmpdir = tempfile.mkdtemp(prefix='composite_kernel_', dir=os.curdir) + tmpdir = coder.createUniqueDir("composite_kernel_") pklpath = os.path.abspath(os.path.join(tmpdir, 'composite.pkl')) pkl.dump(composite, open(pklpath, 'wb')) # make "compiled" composite module nkernels = len(composite.elements()) + # + indent = 4*' ' # 1. scatter element_scatter_method_defs = [ f'scatter_{ik}, scattering_coeff_{ik}, absorb_{ik}, absorption_coeff_{ik} = kernel_funcs_list[{ik}]' - for ik in range(nkernels)] + for ik in range(nkernels) + ] element_scatter_method_defs = '\n'.join(element_scatter_method_defs) - indent = 4*' ' - lines = _create_select_kernel_func_lines( + scatter_method = '\n'.join(_create_scatter_method(nkernels, indent=indent)) + # 2. scattering_coeff + scattering_coeff_method = '\n'.join(_create_scattering_coeff_method(nkernels, indent=4*' ')) + # 3. absorb + absorb_method = '\n'.join(_create_absorb_method(nkernels, indent=indent)) + # 4. absorption_coeff + abs_coeff_method = '\n'.join(_create_abs_coeff_method(nkernels, indent=4*' ')) + # all together + content = template.format(**locals()) + modulepath = os.path.join(tmpdir, 'compiled_composite.py') + open(modulepath, 'wt').write(content) + _modules[composite] = modulepath + return modulepath + +def _create_scattering_coeff_method(nkernels, indent=4*' '): + header = [ + '@cuda.jit(device=True)', + 'def scattering_coeff(neutron):' + ] + from .._numba import coder + loop = coder.unrollLoop( + N = nkernels, + indent = indent, + before_loop = ['r = 0.0'], + in_loop = ['r += scattering_coeff_{i}(neutron)'], + after_loop = ['return r*scale_scattering_coeff'] + ) + return header + loop + +def _create_scatter_method(nkernels, indent=4*' '): + header = [ + '@cuda.jit(device=True)', + 'def scatter(threadindex, rng_states, neutron):' + ] + loop = _create_select_kernel_func_lines( nkernels, method='scatter', args = 'threadindex, rng_states, neutron', indent = indent ) - if_clause_for_scatter_method = '\n'.join([indent + line for line in lines]) - # 2. scattering_coeff - lines = [f'r += scattering_coeff_{i}(neutron)' for i in range(nkernels)] - add_scattering_coeff = '\n'.join(indent+l for l in lines) - # 3. absorb - lines = _create_select_kernel_func_lines( + return header + loop + +def _create_abs_coeff_method(nkernels, indent=4*' '): + header = [ + '@cuda.jit(device=True)', + 'def absorption_coeff(neutron):' + ] + from .._numba import coder + loop = coder.unrollLoop( + N = nkernels, + indent = indent, + before_loop = ['r = 0.0'], + in_loop = ['r += absorption_coeff_{i}(neutron)'], + after_loop = ['return r/Nkernels'] + ) + return header + loop + +def _create_absorb_method(nkernels, indent=4*' '): + header = [ + '@cuda.jit(device=True)', + 'def absorb(threadindex, rng_states, neutron):' + ] + loop = _create_select_kernel_func_lines( nkernels, method='absorb', args = 'threadindex, rng_states, neutron', indent = indent ) - if_clause_for_absorb_method = '\n'.join([indent + line for line in lines]) - # 4. absorption_coeff - lines = [f'r += absorption_coeff_{i}(neutron)' for i in range(nkernels)] - add_absorption_coeff = '\n'.join(indent+l for l in lines) - # all together - content = template.format(**locals()) - modulepath = os.path.join(tmpdir, 'compiled_composite.py') - open(modulepath, 'wt').write(content) - _modules[composite] = modulepath - return modulepath + return header + loop def _create_select_kernel_func_lines(nkernels, method, args, indent=4*' '): - lines = [] - for i in range(nkernels-1): - lead = 'if' - if i>0: lead = 'elif' - lines.append(f'{lead} r < cumulative_weights[{i}]:') - lines.append(f'{indent}ret = {method}_{i}({args})') - lines.append(f'{indent}ikernel={i}') - lines.append(f'else:') - lines.append(f'{indent}ret = {method}_{nkernels-1}({args})') - lines.append(f'{indent}ikernel={nkernels-1}') - lines.append(f'neutron[-1]/=weights[ikernel]') - lines.append(f'return ret') - return lines + from .._numba import coder + return coder.unrollLoop( + N = nkernels, + indent = indent, + before_loop = [ + 'r = xoroshiro128p_uniform_float32(rng_states, threadindex)' + ], + in_loop = [ + 'if r < cumulative_weights[{i}]:', + f' ret = {method}_'+'{i}'+f'({args})', + ' neutron[-1]/=weights[{i}]', + ' return ret' + ], + after_loop = [] + ) template = """import os, pickle as pkl, numpy as np from numba import cuda @@ -98,25 +145,80 @@ def _create_select_kernel_func_lines(nkernels, method, args, indent=4*' '): # and assign names {element_scatter_method_defs} +{scatter_method} + +{absorb_method} + +{scattering_coeff_method} + +{abs_coeff_method} +""" + +# Example code generated +""" +import os, pickle as pkl, numpy as np +from numba import cuda +from mcvine.acc._numba import xoroshiro128p_uniform_float32 + +pklpath = '/home/username/mcvinesim/.mcvine.acc.coder/composite_kernel_g9xcahde/composite.pkl' +composite = pkl.load(open(pklpath, 'rb')) +elements = composite.elements() +Nkernels = len(elements) +weights = [float(element.weight) for element in elements] +cumulative_weights = [weights[0]] +for w in weights[1:]: + cumulative_weights.append(w+cumulative_weights[-1]) +cumulative_weights = np.array(cumulative_weights) +weights = np.array(weights) +scale_scattering_coeff = 1./Nkernels if composite.average else 1 + +# create cuda functions for kernels +from mcvine.acc.kernels import scatter_func_factory +kernel_funcs_list = [] +for element in elements: + kernel_funcs = scatter_func_factory.render(element) + kernel_funcs_list.append(kernel_funcs) + continue + +# and assign names +scatter_0, scattering_coeff_0, absorb_0, absorption_coeff_0 = kernel_funcs_list[0] +scatter_1, scattering_coeff_1, absorb_1, absorption_coeff_1 = kernel_funcs_list[1] + @cuda.jit(device=True) def scatter(threadindex, rng_states, neutron): r = xoroshiro128p_uniform_float32(rng_states, threadindex) -{if_clause_for_scatter_method} + if r < cumulative_weights[0]: + ret = scatter_0(threadindex, rng_states, neutron) + neutron[-1]/=weights[0] + return ret + if r < cumulative_weights[1]: + ret = scatter_1(threadindex, rng_states, neutron) + neutron[-1]/=weights[1] + return ret + +@cuda.jit(device=True) +def absorb(threadindex, rng_states, neutron): + r = xoroshiro128p_uniform_float32(rng_states, threadindex) + if r < cumulative_weights[0]: + ret = absorb_0(threadindex, rng_states, neutron) + neutron[-1]/=weights[0] + return ret + if r < cumulative_weights[1]: + ret = absorb_1(threadindex, rng_states, neutron) + neutron[-1]/=weights[1] + return ret @cuda.jit(device=True) def scattering_coeff(neutron): r = 0.0 -{add_scattering_coeff} + r += scattering_coeff_0(neutron) + r += scattering_coeff_1(neutron) return r*scale_scattering_coeff -@cuda.jit(device=True) -def absorb(threadindex, rng_states, neutron): - r = xoroshiro128p_uniform_float32(rng_states, threadindex) -{if_clause_for_absorb_method} - @cuda.jit(device=True) def absorption_coeff(neutron): r = 0.0 -{add_absorption_coeff} + r += absorption_coeff_0(neutron) + r += absorption_coeff_1(neutron) return r/Nkernels """ From bb617f0ca51338f5cafaf9be6d2b2e8de5e21d81 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Sun, 30 Apr 2023 20:33:10 -0700 Subject: [PATCH 06/53] trivial cosmetic changes: acc/kernels/composite.py --- acc/kernels/composite.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/acc/kernels/composite.py b/acc/kernels/composite.py index d6385ded..90e13170 100644 --- a/acc/kernels/composite.py +++ b/acc/kernels/composite.py @@ -55,8 +55,8 @@ def _create_scattering_coeff_method(nkernels, indent=4*' '): N = nkernels, indent = indent, before_loop = ['r = 0.0'], - in_loop = ['r += scattering_coeff_{i}(neutron)'], - after_loop = ['return r*scale_scattering_coeff'] + in_loop = ['r += scattering_coeff_{i}(neutron)'], + after_loop = ['return r*scale_scattering_coeff'] ) return header + loop @@ -67,8 +67,8 @@ def _create_scatter_method(nkernels, indent=4*' '): ] loop = _create_select_kernel_func_lines( nkernels, - method='scatter', - args = 'threadindex, rng_states, neutron', + method = 'scatter', + args = 'threadindex, rng_states, neutron', indent = indent ) return header + loop @@ -83,8 +83,8 @@ def _create_abs_coeff_method(nkernels, indent=4*' '): N = nkernels, indent = indent, before_loop = ['r = 0.0'], - in_loop = ['r += absorption_coeff_{i}(neutron)'], - after_loop = ['return r/Nkernels'] + in_loop = ['r += absorption_coeff_{i}(neutron)'], + after_loop = ['return r/Nkernels'] ) return header + loop @@ -95,8 +95,8 @@ def _create_absorb_method(nkernels, indent=4*' '): ] loop = _create_select_kernel_func_lines( nkernels, - method='absorb', - args = 'threadindex, rng_states, neutron', + method = 'absorb', + args = 'threadindex, rng_states, neutron', indent = indent ) return header + loop From 805d780f4167d61e08821983c1f0411a9578104c Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Sun, 30 Apr 2023 20:56:10 -0700 Subject: [PATCH 07/53] use method factory instead of pickle: acc/kernels/composite.py --- acc/kernels/composite.py | 191 ++++++++++++++++++++------------------- 1 file changed, 98 insertions(+), 93 deletions(-) diff --git a/acc/kernels/composite.py b/acc/kernels/composite.py index 90e13170..377b7d47 100644 --- a/acc/kernels/composite.py +++ b/acc/kernels/composite.py @@ -9,35 +9,41 @@ def makeKernelMethods(composite): mod = makeKernelModule(composite) import imp m = imp.load_source('composite', mod) - return m.scatter, m.scattering_coeff, m.absorb, m.absorption_coeff + return m.createKernelMethods(composite) _modules = {} def makeKernelModule(composite): "make cuda device methods for the composite kernel" from .._numba import coder if composite in _modules: return _modules[composite] - import pickle as pkl, tempfile # save composite to be loaded by the "compiled" module tmpdir = coder.createUniqueDir("composite_kernel_") - pklpath = os.path.abspath(os.path.join(tmpdir, 'composite.pkl')) - pkl.dump(composite, open(pklpath, 'wb')) # make "compiled" composite module nkernels = len(composite.elements()) # indent = 4*' ' + add_indent = lambda lines, n: [indent*n+l for l in lines] # 1. scatter element_scatter_method_defs = [ f'scatter_{ik}, scattering_coeff_{ik}, absorb_{ik}, absorption_coeff_{ik} = kernel_funcs_list[{ik}]' for ik in range(nkernels) ] - element_scatter_method_defs = '\n'.join(element_scatter_method_defs) - scatter_method = '\n'.join(_create_scatter_method(nkernels, indent=indent)) + element_scatter_method_defs = '\n'.join(add_indent(element_scatter_method_defs, 1)) + scatter_method = '\n'.join(add_indent( + _create_scatter_method(nkernels, indent=indent), 1 + )) # 2. scattering_coeff - scattering_coeff_method = '\n'.join(_create_scattering_coeff_method(nkernels, indent=4*' ')) + scattering_coeff_method = '\n'.join(add_indent( + _create_scattering_coeff_method(nkernels, indent=4*' '), 1 + )) # 3. absorb - absorb_method = '\n'.join(_create_absorb_method(nkernels, indent=indent)) + absorb_method = '\n'.join(add_indent( + _create_absorb_method(nkernels, indent=indent), 1 + )) # 4. absorption_coeff - abs_coeff_method = '\n'.join(_create_abs_coeff_method(nkernels, indent=4*' ')) + abs_coeff_method = '\n'.join(add_indent( + _create_abs_coeff_method(nkernels, indent=4*' '), 1 + )) # all together content = template.format(**locals()) modulepath = os.path.join(tmpdir, 'compiled_composite.py') @@ -118,31 +124,30 @@ def _create_select_kernel_func_lines(nkernels, method, args, indent=4*' '): after_loop = [] ) -template = """import os, pickle as pkl, numpy as np +template = """ +import os, numpy as np from numba import cuda from mcvine.acc._numba import xoroshiro128p_uniform_float32 -pklpath = {pklpath!r} -composite = pkl.load(open(pklpath, 'rb')) -elements = composite.elements() -Nkernels = len(elements) -weights = [float(element.weight) for element in elements] -cumulative_weights = [weights[0]] -for w in weights[1:]: - cumulative_weights.append(w+cumulative_weights[-1]) -cumulative_weights = np.array(cumulative_weights) -weights = np.array(weights) -scale_scattering_coeff = 1./Nkernels if composite.average else 1 - -# create cuda functions for kernels -from mcvine.acc.kernels import scatter_func_factory -kernel_funcs_list = [] -for element in elements: - kernel_funcs = scatter_func_factory.render(element) - kernel_funcs_list.append(kernel_funcs) - continue - -# and assign names +def createKernelMethods(composite): + elements = composite.elements() + Nkernels = len(elements) + weights = [float(element.weight) for element in elements] + cumulative_weights = [weights[0]] + for w in weights[1:]: + cumulative_weights.append(w+cumulative_weights[-1]) + cumulative_weights = np.array(cumulative_weights) + weights = np.array(weights) + scale_scattering_coeff = 1./Nkernels if composite.average else 1 + + # create cuda functions for kernels + from mcvine.acc.kernels import scatter_func_factory + kernel_funcs_list = [] + for element in elements: + kernel_funcs = scatter_func_factory.render(element) + kernel_funcs_list.append(kernel_funcs) + continue + {element_scatter_method_defs} {scatter_method} @@ -152,73 +157,73 @@ def _create_select_kernel_func_lines(nkernels, method, args, indent=4*' '): {scattering_coeff_method} {abs_coeff_method} + return scatter, scattering_coeff, absorb, absorption_coeff """ # Example code generated """ -import os, pickle as pkl, numpy as np +import os, numpy as np from numba import cuda from mcvine.acc._numba import xoroshiro128p_uniform_float32 -pklpath = '/home/username/mcvinesim/.mcvine.acc.coder/composite_kernel_g9xcahde/composite.pkl' -composite = pkl.load(open(pklpath, 'rb')) -elements = composite.elements() -Nkernels = len(elements) -weights = [float(element.weight) for element in elements] -cumulative_weights = [weights[0]] -for w in weights[1:]: - cumulative_weights.append(w+cumulative_weights[-1]) -cumulative_weights = np.array(cumulative_weights) -weights = np.array(weights) -scale_scattering_coeff = 1./Nkernels if composite.average else 1 - -# create cuda functions for kernels -from mcvine.acc.kernels import scatter_func_factory -kernel_funcs_list = [] -for element in elements: - kernel_funcs = scatter_func_factory.render(element) - kernel_funcs_list.append(kernel_funcs) - continue - -# and assign names -scatter_0, scattering_coeff_0, absorb_0, absorption_coeff_0 = kernel_funcs_list[0] -scatter_1, scattering_coeff_1, absorb_1, absorption_coeff_1 = kernel_funcs_list[1] - -@cuda.jit(device=True) -def scatter(threadindex, rng_states, neutron): - r = xoroshiro128p_uniform_float32(rng_states, threadindex) - if r < cumulative_weights[0]: - ret = scatter_0(threadindex, rng_states, neutron) - neutron[-1]/=weights[0] - return ret - if r < cumulative_weights[1]: - ret = scatter_1(threadindex, rng_states, neutron) - neutron[-1]/=weights[1] - return ret - -@cuda.jit(device=True) -def absorb(threadindex, rng_states, neutron): - r = xoroshiro128p_uniform_float32(rng_states, threadindex) - if r < cumulative_weights[0]: - ret = absorb_0(threadindex, rng_states, neutron) - neutron[-1]/=weights[0] - return ret - if r < cumulative_weights[1]: - ret = absorb_1(threadindex, rng_states, neutron) - neutron[-1]/=weights[1] - return ret - -@cuda.jit(device=True) -def scattering_coeff(neutron): - r = 0.0 - r += scattering_coeff_0(neutron) - r += scattering_coeff_1(neutron) - return r*scale_scattering_coeff - -@cuda.jit(device=True) -def absorption_coeff(neutron): - r = 0.0 - r += absorption_coeff_0(neutron) - r += absorption_coeff_1(neutron) - return r/Nkernels +def createKernelMethods(composite): + elements = composite.elements() + Nkernels = len(elements) + weights = [float(element.weight) for element in elements] + cumulative_weights = [weights[0]] + for w in weights[1:]: + cumulative_weights.append(w+cumulative_weights[-1]) + cumulative_weights = np.array(cumulative_weights) + weights = np.array(weights) + scale_scattering_coeff = 1./Nkernels if composite.average else 1 + + # create cuda functions for kernels + from mcvine.acc.kernels import scatter_func_factory + kernel_funcs_list = [] + for element in elements: + kernel_funcs = scatter_func_factory.render(element) + kernel_funcs_list.append(kernel_funcs) + continue + + scatter_0, scattering_coeff_0, absorb_0, absorption_coeff_0 = kernel_funcs_list[0] + scatter_1, scattering_coeff_1, absorb_1, absorption_coeff_1 = kernel_funcs_list[1] + + @cuda.jit(device=True) + def scatter(threadindex, rng_states, neutron): + r = xoroshiro128p_uniform_float32(rng_states, threadindex) + if r < cumulative_weights[0]: + ret = scatter_0(threadindex, rng_states, neutron) + neutron[-1]/=weights[0] + return ret + if r < cumulative_weights[1]: + ret = scatter_1(threadindex, rng_states, neutron) + neutron[-1]/=weights[1] + return ret + + @cuda.jit(device=True) + def absorb(threadindex, rng_states, neutron): + r = xoroshiro128p_uniform_float32(rng_states, threadindex) + if r < cumulative_weights[0]: + ret = absorb_0(threadindex, rng_states, neutron) + neutron[-1]/=weights[0] + return ret + if r < cumulative_weights[1]: + ret = absorb_1(threadindex, rng_states, neutron) + neutron[-1]/=weights[1] + return ret + + @cuda.jit(device=True) + def scattering_coeff(neutron): + r = 0.0 + r += scattering_coeff_0(neutron) + r += scattering_coeff_1(neutron) + return r*scale_scattering_coeff + + @cuda.jit(device=True) + def absorption_coeff(neutron): + r = 0.0 + r += absorption_coeff_0(neutron) + r += absorption_coeff_1(neutron) + return r/Nkernels + return scatter, scattering_coeff, absorb, absorption_coeff """ From 9b6601621829f25e192274103c2eead6fc72d2a8 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Sun, 30 Apr 2023 22:16:34 -0700 Subject: [PATCH 08/53] refactored so that only one module needed for all composite kernels with same number of subkernels: acc/kernels/composite.py --- acc/_numba.py | 15 +++++++++++++++ acc/kernels/composite.py | 10 ++++------ 2 files changed, 19 insertions(+), 6 deletions(-) diff --git a/acc/_numba.py b/acc/_numba.py index 1bf7ab45..1d4c6631 100644 --- a/acc/_numba.py +++ b/acc/_numba.py @@ -14,10 +14,25 @@ def __init__(self, workdir=None): self.workdir = workdir = workdir or os.path.abspath(".mcvine.acc.coder") if not os.path.exists(workdir): os.makedirs(workdir) + self.modules = dict() + + def getModule(self, type, N): + container = coder.createDir(type) + modulepath = os.path.join(container, f'compiled_{N}.py') + key = type, N + if os.path.exists(modulepath) and key not in self.modules: + self.modules[key] = modulepath + return modulepath def createUniqueDir(self, prefix): return tempfile.mkdtemp(prefix=prefix, dir=self.workdir) + def createDir(self, name): + wd = os.path.join(self.workdir, name) + if not os.path.exists(wd): + os.makedirs(wd) + return wd + @classmethod def unrollLoop(cls, N, indent=4*' ', before_loop=None, in_loop=None, after_loop=None): lines=[indent+line for line in before_loop or []] diff --git a/acc/kernels/composite.py b/acc/kernels/composite.py index 377b7d47..c246af16 100644 --- a/acc/kernels/composite.py +++ b/acc/kernels/composite.py @@ -11,15 +11,15 @@ def makeKernelMethods(composite): m = imp.load_source('composite', mod) return m.createKernelMethods(composite) -_modules = {} def makeKernelModule(composite): "make cuda device methods for the composite kernel" from .._numba import coder - if composite in _modules: return _modules[composite] + nkernels = len(composite.elements()) + modulepath = coder.getModule("composite_kernel", nkernels) + if os.path.exists(modulepath): + return modulepath # save composite to be loaded by the "compiled" module - tmpdir = coder.createUniqueDir("composite_kernel_") # make "compiled" composite module - nkernels = len(composite.elements()) # indent = 4*' ' add_indent = lambda lines, n: [indent*n+l for l in lines] @@ -46,9 +46,7 @@ def makeKernelModule(composite): )) # all together content = template.format(**locals()) - modulepath = os.path.join(tmpdir, 'compiled_composite.py') open(modulepath, 'wt').write(content) - _modules[composite] = modulepath return modulepath def _create_scattering_coeff_method(nkernels, indent=4*' '): From df0a2eec1162f5b00045bbeb762352e28ba4daf4 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Mon, 1 May 2023 08:05:52 -0700 Subject: [PATCH 09/53] handle overflow: acc/geometry/arrow_intersect.py --- acc/geometry/arrow_intersect.py | 9 +++++++-- tests/geometry/test_arrow_intersect.py | 17 +++++++++++++++++ 2 files changed, 24 insertions(+), 2 deletions(-) diff --git a/acc/geometry/arrow_intersect.py b/acc/geometry/arrow_intersect.py index 830a933e..7f918fe7 100644 --- a/acc/geometry/arrow_intersect.py +++ b/acc/geometry/arrow_intersect.py @@ -150,8 +150,13 @@ def remove_item(idx, l, N): @cuda.jit(device=True) def insert_into_sorted_list(d, l, N): - 'insert data "d" into existing sorted array (low to high) of length N' - if N>=len(l): return N + '''insert data "d" into existing sorted array (low to high) of length N + and keep the lower number when there are too many elements + ''' + if N>=len(l): + if l[N-1] <= d: + return N + N = N - 1 if N==0: l[0] = d return 1 diff --git a/tests/geometry/test_arrow_intersect.py b/tests/geometry/test_arrow_intersect.py index 696c30a0..23de45e6 100755 --- a/tests/geometry/test_arrow_intersect.py +++ b/tests/geometry/test_arrow_intersect.py @@ -58,6 +58,23 @@ def test_insert_into_sorted_list(): assert N==6 np.testing.assert_array_equal(l[:N], [0., 1., 2., 3., 4., 8.]) + l = np.arange(4.) + d = 8 + N = arrow_intersect.insert_into_sorted_list(d, l, 4) + assert N==4 + np.testing.assert_array_equal(l[:N], [0., 1., 2., 3.]) + + l = np.arange(4.) + d = -1 + N = arrow_intersect.insert_into_sorted_list(d, l, 4) + assert N==4 + np.testing.assert_array_equal(l[:N], [-1., 0., 1., 2.]) + + l = np.arange(4.) + d = 1.5 + N = arrow_intersect.insert_into_sorted_list(d, l, 4) + assert N==4 + np.testing.assert_array_equal(l[:N], [0., 1., 1.5, 2.]) return def main(): From 2ad8744b1ef656e45619c156a5a32bd0dc5141ff Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Mon, 1 May 2023 20:12:45 -0700 Subject: [PATCH 10/53] started implementing acc.geometry.composite with an example methods factory. tested `forward_intersect_all` --- acc/geometry/composite.py | 97 +++++++++++++++++++++++++ tests/geometry/test_composite_shape.py | 25 +++++++ tests/geometry/union_three_elements.xml | 11 +++ 3 files changed, 133 insertions(+) create mode 100644 acc/geometry/composite.py create mode 100644 tests/geometry/test_composite_shape.py create mode 100644 tests/geometry/union_three_elements.xml diff --git a/acc/geometry/composite.py b/acc/geometry/composite.py new file mode 100644 index 00000000..329a6da3 --- /dev/null +++ b/acc/geometry/composite.py @@ -0,0 +1,97 @@ +def makeMethods(shapes): + mod = makeModule(shapes) + import imp + m = imp.load_source('mod', mod) + return m.createKernelMethods(composite) + +def makeModule(shapes): + "make cuda device methods for shapes" + from .._numba import coder + nshapes = len(shapes) + modulepath = coder.getModule("shapes", nshapes) + if os.path.exists(modulepath): + return modulepath + indent = 4*' ' + add_indent = lambda lines, n: [indent*n+l for l in lines] + element_scatter_method_defs = [ + f'scatter_{ik}, scattering_coeff_{ik}, absorb_{ik}, absorption_coeff_{ik} = kernel_funcs_list[{ik}]' + for ik in range(nkernels) + ] + + +# Example code +import os, numpy as np +from numba import cuda +from mcvine.acc._numba import xoroshiro128p_uniform_float32 +from mcvine.acc import test + +def createMethods_3(shapes): + assert len(shapes)==3 + Nshapes = len(shapes) + from mcvine.acc.geometry import arrow_intersect + from mcvine.acc.geometry.location import inside, onborder + from mcvine.acc.geometry.arrow_intersect import max_intersections, insert_into_sorted_list + funcs_list = [] + for shape in shapes: + intersect = arrow_intersect.arrow_intersect_func_factory.render(shape) + locate = arrow_intersect.locate_func_factory.render(shape) + funcs = locate, intersect + funcs_list.append(funcs) + continue + + locate_0, intersect_0 = funcs_list[0] + locate_1, intersect_1 = funcs_list[1] + locate_2, intersect_2 = funcs_list[2] + + @cuda.jit(device=True) + def _forward_intersect_all(x,y,z, vx,vy,vz, ts, ts0, ts1, ts2): + N0 = intersect_0(x,y,z, vx,vy,vz, ts0) + N1 = intersect_1(x,y,z, vx,vy,vz, ts1) + N2 = intersect_2(x,y,z, vx,vy,vz, ts2) + N = 0 + for i in range(N0): + t = ts0[i] + if t>0: + N = insert_into_sorted_list(t, ts, N) + for i in range(N1): + t = ts1[i] + if t>0: + N = insert_into_sorted_list(t, ts, N) + for i in range(N2): + t = ts2[i] + if t>0: + N = insert_into_sorted_list(t, ts, N) + return N + + @cuda.jit(device=True) + def find_shape_containing_point(x,y,z): + if locate_0(x,y,z) == inside: + return 0 + if locate_1(x,y,z) == inside: + return 1 + if locate_2(x,y,z) == inside: + return 2 + + @cuda.jit(device=True) + def find_1st_hit(start, direction): + raise NotImplementedError + + if test.USE_CUDASIM: + @cuda.jit(device=True) + def forward_intersect_all(x,y,z, vx,vy,vz, ts): + ts0 = np.zeros(max_intersections, dtype=float) + ts1 = np.zeros(max_intersections, dtype=float) + ts2 = np.zeros(max_intersections, dtype=float) + return _forward_intersect_all(x,y,z, vx,vy,vz, ts, ts0, ts1, ts2) + else: + @cuda.jit(device=True) + def forward_intersect_all(x,y,z, vx,vy,vz, ts): + ts0 = cuda.local.array(max_intersections, dtype=numba.float64) + ts1 = cuda.local.array(max_intersections, dtype=numba.float64) + ts2 = cuda.local.array(max_intersections, dtype=numba.float64) + return _forward_intersect_all(x,y,z, vx,vy,vz, ts, ts0, ts1, ts2) + return dict( + forward_intersect_all = forward_intersect_all, + find_shape_containing_point = find_shape_containing_point, + find_1st_hit = find_1st_hit, + ) diff --git a/tests/geometry/test_composite_shape.py b/tests/geometry/test_composite_shape.py new file mode 100644 index 00000000..1f159a5a --- /dev/null +++ b/tests/geometry/test_composite_shape.py @@ -0,0 +1,25 @@ +#!/usr/bin/env python + +import pytest, os, time +import math, numpy as np +from numba import cuda + +from instrument.nixml import parse_file +from mcvine.acc import test +from mcvine.acc.geometry import locate, location, arrow_intersect +from mcvine.acc.geometry.composite import createMethods_3 + +thisdir = os.path.dirname(__file__) + +# device functions can be tested with CUDASIM only +@pytest.mark.skipif(not test.USE_CUDASIM, reason='no CUDASIM') +def test_forward_intersect_all(): + union = parse_file(os.path.join(thisdir, 'union_three_elements.xml'))[0] + shapes = union.shapes + methods = createMethods_3(shapes) + forward_intersect_all = methods['forward_intersect_all'] + ts = np.zeros(5) + assert forward_intersect_all(0.,0.,0., 1.,0.,0., ts) == 3 + expected = np.array([0.004,0.01,0.025, 0.,0.]) + assert np.allclose(ts, expected) + return diff --git a/tests/geometry/union_three_elements.xml b/tests/geometry/union_three_elements.xml new file mode 100644 index 00000000..d6511ff9 --- /dev/null +++ b/tests/geometry/union_three_elements.xml @@ -0,0 +1,11 @@ + + + + + + + + + + + From 4de8a86d4895c6d551d693568ae3ffcfef2faf07 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Wed, 3 May 2023 17:38:41 -0700 Subject: [PATCH 11/53] implemented `find_1st_hit` and added a test: acc/geometry/composite.py --- acc/geometry/composite.py | 83 ++++++++++++++++++++++++- tests/geometry/test_composite_shape.py | 21 +++++-- tests/geometry/union_three_elements.xml | 12 +++- 3 files changed, 107 insertions(+), 9 deletions(-) diff --git a/acc/geometry/composite.py b/acc/geometry/composite.py index 329a6da3..a283c2d0 100644 --- a/acc/geometry/composite.py +++ b/acc/geometry/composite.py @@ -1,3 +1,72 @@ +from mcvine.acc.geometry.arrow_intersect import max_intersections + +def make_find_1st_hit(forward_intersect_all, is_onborder, find_shape_containing_point): + @cuda.jit(device=True) + def _find_1st_hit(x,y,z, vx,vy,vz, ts): + nIntersections = forward_intersect_all(x,y,z, vx,vy,vz, ts) + # we have two cases + # case 1 + # shape1 | vacuum | shape2 + # start-|---------|---> + # case 2 + # vacuum | shape1 + # start-|-----------> + # case 1: there will be odd number of intersections + # case 2: there will be even number of intersections + # we just need to determine which one of the above two cases is true, + # case1: + if nIntersections % 2 == 1: + ret = find_shape_containing_point( + x + ts[0]/2.*vx, + y + ts[0]/2.*vy, + z + ts[0]/2.*vz, + ) + else : + # case2: + # no intersection + if nIntersections==0: + return -1 + # at least two + midt = (ts[0]+ts[1])/2. + ret = find_shape_containing_point( + x + midt*vx, + y + midt*vy, + z + midt*vz, + ) + # let us determine if the start is on border + isonborder = is_onborder(x,y,z) + # If start is not on border of any shape, it would be easier. + if (not isonborder): return ret + + # on border. that is a bit more difficult. + # we need to go over all intersection pairs, and find the first pair + # whose midlle point is insde a shape. That shape containing + # the middle point is the target. + previous = 0.0 + for point_index in range(nIntersections): + now = ts[point_index] + midt = (previous+now)/2 + ret = find_shape_containing_point( + x + midt*vx, + y + midt*vy, + z + midt*vz, + ) + if ret>=0: return ret + previous = now + return -1 + if test.USE_CUDASIM: + @cuda.jit(device=True) + def find_1st_hit(x,y,z, vx,vy,vz): + ts = np.zeros(max_intersections, dtype=float) + return _find_1st_hit(x,y,z, vx,vy,vz, ts) + else: + @cuda.jit(device=True) + def find_1st_hit(x,y,z, vx,vy,vz): + ts = cuda.local.array(max_intersections, dtype=numba.float64) + return _find_1st_hit(x,y,z, vx,vy,vz, ts) + return find_1st_hit + + def makeMethods(shapes): mod = makeModule(shapes) import imp @@ -71,10 +140,18 @@ def find_shape_containing_point(x,y,z): return 1 if locate_2(x,y,z) == inside: return 2 + return -1 @cuda.jit(device=True) - def find_1st_hit(start, direction): - raise NotImplementedError + def is_onborder(x,y,z): + "check if the point is on the border of any of the shapes" + if locate_0(x,y,z) == onborder: + return True + if locate_1(x,y,z) == onborder: + return True + if locate_2(x,y,z) == onborder: + return True + return False if test.USE_CUDASIM: @cuda.jit(device=True) @@ -93,5 +170,5 @@ def forward_intersect_all(x,y,z, vx,vy,vz, ts): return dict( forward_intersect_all = forward_intersect_all, find_shape_containing_point = find_shape_containing_point, - find_1st_hit = find_1st_hit, + is_onborder = is_onborder, ) diff --git a/tests/geometry/test_composite_shape.py b/tests/geometry/test_composite_shape.py index 1f159a5a..e5f9d99d 100644 --- a/tests/geometry/test_composite_shape.py +++ b/tests/geometry/test_composite_shape.py @@ -7,7 +7,7 @@ from instrument.nixml import parse_file from mcvine.acc import test from mcvine.acc.geometry import locate, location, arrow_intersect -from mcvine.acc.geometry.composite import createMethods_3 +from mcvine.acc.geometry.composite import createMethods_3, make_find_1st_hit thisdir = os.path.dirname(__file__) @@ -18,8 +18,21 @@ def test_forward_intersect_all(): shapes = union.shapes methods = createMethods_3(shapes) forward_intersect_all = methods['forward_intersect_all'] - ts = np.zeros(5) - assert forward_intersect_all(0.,0.,0., 1.,0.,0., ts) == 3 - expected = np.array([0.004,0.01,0.025, 0.,0.]) + ts = np.zeros(6) + assert forward_intersect_all(0.,0.,0., 1.,0.,0., ts) == 5 + expected = np.array([0.025,0.099,0.1, 0.199,0.2, 0]) assert np.allclose(ts, expected) return + +@pytest.mark.skipif(not test.USE_CUDASIM, reason='no CUDASIM') +def test_find_1st_hit(): + union = parse_file(os.path.join(thisdir, 'union_three_elements.xml'))[0] + shapes = union.shapes + methods = createMethods_3(shapes) + find_1st_hit = make_find_1st_hit(**methods) + assert find_1st_hit(-0.25,0.,0., 1.,0.,0.) == 0 + assert find_1st_hit(-0.15,0.,0., 1.,0.,0.) == 1 + assert find_1st_hit(-0.05,0.,0., 1.,0.,0.) == 2 + assert find_1st_hit(0.05,0.,0., 1.,0.,0.) == 1 + assert find_1st_hit(0.15,0.,0., 1.,0.,0.) == 0 + assert find_1st_hit(0.25,0.,0., 1.,0.,0.) == -1 diff --git a/tests/geometry/union_three_elements.xml b/tests/geometry/union_three_elements.xml index d6511ff9..b4112562 100644 --- a/tests/geometry/union_three_elements.xml +++ b/tests/geometry/union_three_elements.xml @@ -3,9 +3,17 @@ - + + + + + + + + + + - From 5369c923738e1e8f1fbac3f3a0252e0b8eddd383 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Sun, 7 May 2023 08:35:58 -0700 Subject: [PATCH 12/53] fixed assertion: tests/components/samples/test_homogeneous_multiple_scatterer.py --- tests/components/samples/test_homogeneous_multiple_scatterer.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/components/samples/test_homogeneous_multiple_scatterer.py b/tests/components/samples/test_homogeneous_multiple_scatterer.py index 2a1c2c39..2c569eeb 100755 --- a/tests/components/samples/test_homogeneous_multiple_scatterer.py +++ b/tests/components/samples/test_homogeneous_multiple_scatterer.py @@ -63,7 +63,7 @@ def run_kernel(rng_states, out_neutrons, N, neutron): rng_states = create_xoroshiro128p_states(1, seed=0) N = np.zeros(1, dtype=int) rt = run_kernel[1,1](rng_states, out_neutrons, N, neutron) - assert N[0] == HMS.max_ms_loops_path1 + 1 + assert N[0] == HMS.max_ms_loops_path1 print(out_neutrons[:N[0]]) return From d15818158006bf09fda16e99b4975e51764065b5 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Sun, 7 May 2023 08:39:28 -0700 Subject: [PATCH 13/53] set pytest marker: tests/scatterers/test_homogeneous_scatterer.py --- tests/scatterers/test_homogeneous_scatterer.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/tests/scatterers/test_homogeneous_scatterer.py b/tests/scatterers/test_homogeneous_scatterer.py index 5f12ea2e..b76ad31c 100644 --- a/tests/scatterers/test_homogeneous_scatterer.py +++ b/tests/scatterers/test_homogeneous_scatterer.py @@ -1,7 +1,10 @@ import os thisdir = os.path.dirname(__file__) import numpy as np +import pytest +from mcvine.acc import test +@pytest.mark.skipif(not test.USE_CUDASIM, reason='no CUDASIM') def test(): path = os.path.join(thisdir, "sampleassemblies", 'isotropic_sphere', 'sampleassembly.xml') from mcvine.acc.components.samples import loadFirstHomogeneousScatterer From 0c6ee4d8fc68ff678d78008f2aed9f2b9adc07e4 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Sun, 7 May 2023 09:34:08 -0700 Subject: [PATCH 14/53] geometry.composite: added method "intersect_all" --- acc/geometry/composite.py | 29 +++++++++++++++++++++++++- tests/geometry/test_composite_shape.py | 12 +++++++++++ 2 files changed, 40 insertions(+), 1 deletion(-) diff --git a/acc/geometry/composite.py b/acc/geometry/composite.py index a283c2d0..4e5cb2c2 100644 --- a/acc/geometry/composite.py +++ b/acc/geometry/composite.py @@ -1,6 +1,6 @@ from mcvine.acc.geometry.arrow_intersect import max_intersections -def make_find_1st_hit(forward_intersect_all, is_onborder, find_shape_containing_point): +def make_find_1st_hit(forward_intersect_all, is_onborder, find_shape_containing_point, **kwds): @cuda.jit(device=True) def _find_1st_hit(x,y,z, vx,vy,vz, ts): nIntersections = forward_intersect_all(x,y,z, vx,vy,vz, ts) @@ -112,6 +112,20 @@ def createMethods_3(shapes): locate_1, intersect_1 = funcs_list[1] locate_2, intersect_2 = funcs_list[2] + @cuda.jit(device=True) + def _intersect_all(x,y,z, vx,vy,vz, ts, ts0, ts1, ts2): + N0 = intersect_0(x,y,z, vx,vy,vz, ts0) + N1 = intersect_1(x,y,z, vx,vy,vz, ts1) + N2 = intersect_2(x,y,z, vx,vy,vz, ts2) + N = 0 + for i in range(N0): + N = insert_into_sorted_list(ts0[i], ts, N) + for i in range(N1): + N = insert_into_sorted_list(ts1[i], ts, N) + for i in range(N2): + N = insert_into_sorted_list(ts2[i], ts, N) + return N + @cuda.jit(device=True) def _forward_intersect_all(x,y,z, vx,vy,vz, ts, ts0, ts1, ts2): N0 = intersect_0(x,y,z, vx,vy,vz, ts0) @@ -155,6 +169,12 @@ def is_onborder(x,y,z): if test.USE_CUDASIM: @cuda.jit(device=True) + def intersect_all(x,y,z, vx,vy,vz, ts): + ts0 = np.zeros(max_intersections, dtype=float) + ts1 = np.zeros(max_intersections, dtype=float) + ts2 = np.zeros(max_intersections, dtype=float) + return _intersect_all(x,y,z, vx,vy,vz, ts, ts0, ts1, ts2) + @cuda.jit(device=True) def forward_intersect_all(x,y,z, vx,vy,vz, ts): ts0 = np.zeros(max_intersections, dtype=float) ts1 = np.zeros(max_intersections, dtype=float) @@ -162,12 +182,19 @@ def forward_intersect_all(x,y,z, vx,vy,vz, ts): return _forward_intersect_all(x,y,z, vx,vy,vz, ts, ts0, ts1, ts2) else: @cuda.jit(device=True) + def intersect_all(x,y,z, vx,vy,vz, ts): + ts0 = cuda.local.array(max_intersections, dtype=numba.float64) + ts1 = cuda.local.array(max_intersections, dtype=numba.float64) + ts2 = cuda.local.array(max_intersections, dtype=numba.float64) + return _intersect_all(x,y,z, vx,vy,vz, ts, ts0, ts1, ts2) + @cuda.jit(device=True) def forward_intersect_all(x,y,z, vx,vy,vz, ts): ts0 = cuda.local.array(max_intersections, dtype=numba.float64) ts1 = cuda.local.array(max_intersections, dtype=numba.float64) ts2 = cuda.local.array(max_intersections, dtype=numba.float64) return _forward_intersect_all(x,y,z, vx,vy,vz, ts, ts0, ts1, ts2) return dict( + intersect_all = intersect_all, forward_intersect_all = forward_intersect_all, find_shape_containing_point = find_shape_containing_point, is_onborder = is_onborder, diff --git a/tests/geometry/test_composite_shape.py b/tests/geometry/test_composite_shape.py index e5f9d99d..10593495 100644 --- a/tests/geometry/test_composite_shape.py +++ b/tests/geometry/test_composite_shape.py @@ -12,6 +12,18 @@ thisdir = os.path.dirname(__file__) # device functions can be tested with CUDASIM only +@pytest.mark.skipif(not test.USE_CUDASIM, reason='no CUDASIM') +def test_intersect_all(): + union = parse_file(os.path.join(thisdir, 'union_three_elements.xml'))[0] + shapes = union.shapes + methods = createMethods_3(shapes) + intersect_all = methods['intersect_all'] + ts = np.zeros(11) + assert intersect_all(0.,0.,0., 1.,0.,0., ts) == 10 + expected = np.array([-0.2,-0.199, -0.1,-0.099, -0.025, 0.025, 0.099,0.1, 0.199,0.2, 0]) + assert np.allclose(ts, expected) + return + @pytest.mark.skipif(not test.USE_CUDASIM, reason='no CUDASIM') def test_forward_intersect_all(): union = parse_file(os.path.join(thisdir, 'union_three_elements.xml'))[0] From 2dc4f88611645c44309509c6a15b20ed86cd788c Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Sun, 7 May 2023 09:59:14 -0700 Subject: [PATCH 15/53] geometry.composite: added locate method for union of 3 shapes and a test --- acc/geometry/composite.py | 31 ++++++++++++++++++++++++++ tests/geometry/test_composite_shape.py | 15 +++++++++++++ 2 files changed, 46 insertions(+) diff --git a/acc/geometry/composite.py b/acc/geometry/composite.py index 4e5cb2c2..058caee1 100644 --- a/acc/geometry/composite.py +++ b/acc/geometry/composite.py @@ -95,6 +95,7 @@ def makeModule(shapes): from mcvine.acc import test def createMethods_3(shapes): + "methods for a list of non-overlapping shapes" assert len(shapes)==3 Nshapes = len(shapes) from mcvine.acc.geometry import arrow_intersect @@ -199,3 +200,33 @@ def forward_intersect_all(x,y,z, vx,vy,vz, ts): find_shape_containing_point = find_shape_containing_point, is_onborder = is_onborder, ) + +def createUnionLocateMethod_3(shapes): + assert len(shapes)==3 + Nshapes = len(shapes) + from mcvine.acc.geometry import arrow_intersect + from mcvine.acc.geometry.location import inside, onborder, outside + locates = [] + for shape in shapes: + locate = arrow_intersect.locate_func_factory.render(shape) + locates.append(locate) + continue + + locate_0 = locates[0] + locate_1 = locates[1] + locate_2 = locates[2] + + @cuda.jit(device=True) + def locate(x,y,z): + loc0 = locate_0(x,y,z) + loc1 = locate_1(x,y,z) + loc2 = locate_2(x,y,z) + if loc0 == inside: return inside + if loc1 == inside: return inside + if loc2 == inside: return inside + if loc0 == onborder: return onborder + if loc1 == onborder: return onborder + if loc2 == onborder: return onborder + return outside + + return locate diff --git a/tests/geometry/test_composite_shape.py b/tests/geometry/test_composite_shape.py index 10593495..0dd9213b 100644 --- a/tests/geometry/test_composite_shape.py +++ b/tests/geometry/test_composite_shape.py @@ -12,6 +12,21 @@ thisdir = os.path.dirname(__file__) # device functions can be tested with CUDASIM only +@pytest.mark.skipif(not test.USE_CUDASIM, reason='no CUDASIM') +def test_locate(): + union = parse_file(os.path.join(thisdir, 'union_three_elements.xml'))[0] + shapes = union.shapes + from mcvine.acc.geometry.composite import createUnionLocateMethod_3 + locate_u3 = createUnionLocateMethod_3(shapes) + assert locate_u3(0, 0, 0) == locate.inside + assert locate_u3(0.025, 0, 0) == locate.onborder + assert locate_u3(0.099, 0, 0) == locate.onborder + assert locate_u3(0.1, 0, 0) == locate.onborder + assert locate_u3(0.199, 0, 0) == locate.onborder + assert locate_u3(0.2, 0, 0) == locate.onborder + assert locate_u3(0.3, 0, 0) == locate.outside + return + @pytest.mark.skipif(not test.USE_CUDASIM, reason='no CUDASIM') def test_intersect_all(): union = parse_file(os.path.join(thisdir, 'union_three_elements.xml'))[0] From 958dbf697750e8dad25b0b5dfc1e33e98e24131e Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Sun, 7 May 2023 10:55:34 -0700 Subject: [PATCH 16/53] partial implementation to handle translation and rotation for locate/intersect --- acc/geometry/arrow_intersect.py | 24 ++++++++++++++++++++++++ acc/geometry/locate.py | 26 +++++++++++++++++++++++++- 2 files changed, 49 insertions(+), 1 deletion(-) diff --git a/acc/geometry/arrow_intersect.py b/acc/geometry/arrow_intersect.py index 7f918fe7..82455ad1 100644 --- a/acc/geometry/arrow_intersect.py +++ b/acc/geometry/arrow_intersect.py @@ -26,6 +26,14 @@ def render(self, shape): return shape.identify(self) def onUnion(self, u): + nelements = len(u.shapes) + if nelements == 2: return self.onUnion2(u) + if nelements == 3: + from .composite import createMethods_3 + return createMethods_3(u.shapes)['intersect_all'] + raise NotImplementedError(f"locate for union of {nelements} elements") + + def onUnion2(self, u): locate1 = self.locate_func_factory.onUnion(u) s1, s2 = u.shapes f1 = s1.identify(self) @@ -85,6 +93,22 @@ def intersectBlock(x, y, z, vx, vy, vz, ts): return intersectBlock + def onTranslation(self, t): + x0,y0,z0 = [_/units.length.meter for _ in t.vector] + body = t.body + intersect_body = self.render(body) + @cuda.jit(device=True, inline=True) + def intersectTranslation(x,y,z, vx,vy,vz, ts): + return intersect_body(x-x0,y-y0,z-z0, vx,vy,vz, ts) + return intersectTranslation + + def onRotation(self, r): + euler_angles = [_/units.angle.rad for _ in r.euler_angles] + # hack + if not np.allclose(euler_angles, [0,0,0]): + raise NotImplementedError + return self.render(r.body) + def onDifference(self, s): locate1 = self.locate_func_factory.onDifference(s) f1 = s.op1.identify(self) diff --git a/acc/geometry/locate.py b/acc/geometry/locate.py index 9fd489a9..c67d3098 100644 --- a/acc/geometry/locate.py +++ b/acc/geometry/locate.py @@ -1,3 +1,4 @@ +import numpy as np from numba import cuda from mcni import units from . import epsilon, location @@ -11,6 +12,14 @@ def render(self, shape): return shape.identify(self) def onUnion(self, u): + nelements = len(u.shapes) + if nelements == 2: return self.onUnion2(u) + if nelements == 3: + from .composite import createUnionLocateMethod_3 + return createUnionLocateMethod_3(u.shapes) + raise NotImplementedError(f"locate for union of {nelements} elements") + + def onUnion2(self, u): s1, s2 = u.shapes f1 = s1.identify(self) f2 = s2.identify(self) @@ -45,6 +54,22 @@ def locateWrtBlock(x, y, z): return locateWrtBlock + def onTranslation(self, t): + x0,y0,z0 = [_/units.length.meter for _ in t.vector] + body = t.body + locate_in_body = self.render(body) + @cuda.jit(device=True, inline=True) + def locateWrtTranslation(x,y,z): + return locate_in_body(x-x0, y-y0, z-z0) + return locateWrtTranslation + + def onRotation(self, r): + euler_angles = [_/units.angle.rad for _ in r.euler_angles] + # hack + if not np.allclose(euler_angles, [0,0,0]): + raise NotImplementedError + return self.render(r.body) + def onDifference(self, s): f1 = s.op1.identify(self) f2 = s.op2.identify(self) @@ -52,7 +77,6 @@ def onDifference(self, s): @cuda.jit(device=True, inline=True) def locateWrtDifference(x, y, z): return cu_device_locate_wrt_difference(x, y, z, f1, f2) - return locateWrtDifference def onIntersection(self, u): From a95a79d65e1b1d49d7b84f64186c42366948ed70 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Sun, 7 May 2023 13:24:24 -0700 Subject: [PATCH 17/53] expose `calculate_attenuation`: acc/scatterers/homogeneous_scatterer.py --- acc/scatterers/homogeneous_scatterer.py | 5 ++--- tests/scatterers/test_homogeneous_scatterer.py | 11 ++++++++++- 2 files changed, 12 insertions(+), 4 deletions(-) diff --git a/acc/scatterers/homogeneous_scatterer.py b/acc/scatterers/homogeneous_scatterer.py index d4cd7e79..0988464a 100644 --- a/acc/scatterers/homogeneous_scatterer.py +++ b/acc/scatterers/homogeneous_scatterer.py @@ -32,8 +32,6 @@ NB_FLOAT = get_numba_floattype() -category = 'samples' - def factory(shape, kernel, mcweights, packing_factor): w_absorption, w_scattering, w_transmission = mcweights from ..geometry import arrow_intersect @@ -116,5 +114,6 @@ def interact_path1(threadindex, rng_states, neutron): tmp_neutron = cuda.local.array(10, dtype=numba.float64) _interact_path1(threadindex, rng_states, neutron, tmp_neutron) return dict( - interact_path1 = interact_path1 + interact_path1 = interact_path1, + calculate_attenuation = calculate_attenuation, ) diff --git a/tests/scatterers/test_homogeneous_scatterer.py b/tests/scatterers/test_homogeneous_scatterer.py index b76ad31c..547e7d49 100644 --- a/tests/scatterers/test_homogeneous_scatterer.py +++ b/tests/scatterers/test_homogeneous_scatterer.py @@ -5,7 +5,7 @@ from mcvine.acc import test @pytest.mark.skipif(not test.USE_CUDASIM, reason='no CUDASIM') -def test(): +def test_homogeneous_scatterer(): path = os.path.join(thisdir, "sampleassemblies", 'isotropic_sphere', 'sampleassembly.xml') from mcvine.acc.components.samples import loadFirstHomogeneousScatterer hs = loadFirstHomogeneousScatterer(path) @@ -34,3 +34,12 @@ def test(): assert p<.8 and p>.7 assert absorbed < 0.36 and absorbed > 0.3 return + +@pytest.mark.skipif(not test.USE_CUDASIM, reason='no CUDASIM') +def test_scatter_func_factory(): + path = os.path.join(thisdir, "sampleassemblies", 'isotropic_sphere', 'sampleassembly.xml') + from mcvine.acc.components.samples import loadFirstHomogeneousScatterer + hs = loadFirstHomogeneousScatterer(path) + from mcvine.acc.scatterers import scatter_func_factory + methods = scatter_func_factory.render(hs) + return From 7df43ed36bda3c9d29e6058b97b24c3cfdcfb877 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Sun, 7 May 2023 13:27:16 -0700 Subject: [PATCH 18/53] composite_scatterer.intersect_path1 preliminary implementation and test --- acc/geometry/composite.py | 2 +- acc/scatterers/__init__.py | 11 +- acc/scatterers/composite_scatterer.py | 133 ++++++++++++++++++ .../sampleassemblies/sample+2cylinders/Al.xyz | 6 + .../sample+2cylinders/sample-scatterer.xml | 10 ++ .../sample+2cylinders/sampleassembly.xml | 47 +++++++ .../sampleenv1-scatterer.xml | 10 ++ .../sampleenv2-scatterer.xml | 10 ++ tests/scatterers/test_composite_scatterer.py | 35 +++++ tests/scatterers/test_scatter_func_factory.py | 14 ++ 10 files changed, 274 insertions(+), 4 deletions(-) create mode 100644 acc/scatterers/composite_scatterer.py create mode 100644 tests/components/samples/sampleassemblies/sample+2cylinders/Al.xyz create mode 100644 tests/components/samples/sampleassemblies/sample+2cylinders/sample-scatterer.xml create mode 100644 tests/components/samples/sampleassemblies/sample+2cylinders/sampleassembly.xml create mode 100644 tests/components/samples/sampleassemblies/sample+2cylinders/sampleenv1-scatterer.xml create mode 100644 tests/components/samples/sampleassemblies/sample+2cylinders/sampleenv2-scatterer.xml create mode 100644 tests/scatterers/test_composite_scatterer.py create mode 100644 tests/scatterers/test_scatter_func_factory.py diff --git a/acc/geometry/composite.py b/acc/geometry/composite.py index 058caee1..6618ca0e 100644 --- a/acc/geometry/composite.py +++ b/acc/geometry/composite.py @@ -16,7 +16,7 @@ def _find_1st_hit(x,y,z, vx,vy,vz, ts): # we just need to determine which one of the above two cases is true, # case1: if nIntersections % 2 == 1: - ret = find_shape_containing_point( + return find_shape_containing_point( x + ts[0]/2.*vx, y + ts[0]/2.*vy, z + ts[0]/2.*vz, diff --git a/acc/scatterers/__init__.py b/acc/scatterers/__init__.py index 3e30290c..50bed903 100644 --- a/acc/scatterers/__init__.py +++ b/acc/scatterers/__init__.py @@ -7,8 +7,7 @@ class ScatterFuncFactory: def render(self, scatterer): """returns cuda device function tuple (scatter, interact_path1) """ - scatter, interact_path1 = scatterer.identify(self) - return scatter, interact_path1 + return scatterer.identify(self) def onCompositeKernel(self, composite): elements = composite.elements() @@ -23,7 +22,13 @@ def onCompositeKernel(self, composite): return makeKernelMethods(composite) def onHomogeneousScatterer(self, hs): - return + from .homogeneous_scatterer import factory + shape = hs.shape() + kernel = hs.kernel() + mcweights = hs.mcweights + packing_factor = hs.packing_factor + methods = factory(shape, kernel, mcweights, packing_factor) + return methods scatter_func_factory = ScatterFuncFactory() diff --git a/acc/scatterers/composite_scatterer.py b/acc/scatterers/composite_scatterer.py new file mode 100644 index 00000000..ec5be283 --- /dev/null +++ b/acc/scatterers/composite_scatterer.py @@ -0,0 +1,133 @@ +# -*- python -*- +# +# Jiao Lin +# + +import numpy as np, numba +from numba import cuda, void, int64 +from mcvine.acc._numba import xoroshiro128p_uniform_float32 +from math import sqrt, exp + +from .interaction_types import absorption, scattering, none +from .. import test +from ..neutron import absorb, prop_dt_inplace, clone +from ..geometry.arrow_intersect import max_intersections +from ..geometry.locate import inside, outside, onborder +from ..geometry.propagation import makePropagateMethods + + +from numba.core import config +if not config.ENABLE_CUDASIM: + from numba.cuda.compiler import Dispatcher, DeviceFunction + +from ..config import get_numba_floattype +NB_FLOAT = get_numba_floattype() + + +def factory_3(composite): + # methods to deal with the whole shape of the composite + shape = composite.shape() + from ..geometry import arrow_intersect + intersect = arrow_intersect.arrow_intersect_func_factory.render(shape) + locate = arrow_intersect.locate_func_factory.render(shape) + propagate_methods = makePropagateMethods(intersect, locate) + propagate_to_next_incident_surface = propagate_methods['propagate_to_next_incident_surface'] + propagate_to_next_exiting_surface = propagate_methods['propagate_to_next_exiting_surface'] + tof_before_first_exit = propagate_methods['tof_before_first_exit'] + forward_distance_in_shape = propagate_methods["forward_distance_in_shape"] + # elements + elements = composite.elements() + nelements = len(elements) + assert nelements == 3 + # find_1st_hit + shapes = [e.shape() for e in elements] + from mcvine.acc.geometry.composite import createMethods_3, make_find_1st_hit + methods = createMethods_3(shapes) + find_1st_hit = make_find_1st_hit(**methods) + # methods for element scatterers + from . import scatter_func_factory + #element_scatter_methods = [scatter_func_factory.render(e) for e in elements] + element_scatter_methods = [] + for e in elements: + m = scatter_func_factory.render(e) + element_scatter_methods.append(m) + element0_interact_path1 = element_scatter_methods[0]['interact_path1'] + element1_interact_path1 = element_scatter_methods[1]['interact_path1'] + element2_interact_path1 = element_scatter_methods[2]['interact_path1'] + element0_calculate_attenuation = element_scatter_methods[0]['calculate_attenuation'] + element1_calculate_attenuation = element_scatter_methods[1]['calculate_attenuation'] + element2_calculate_attenuation = element_scatter_methods[2]['calculate_attenuation'] + + def _interact_path1(threadindex, rng_states, neutron, tmp_neutron): + x,y,z = neutron[:3] + vx,vy,vz = neutron[3:6] + scatterer_index = find_1st_hit(x,y,z, vx,vy,vz) + if scatterer_index < 0 or scatterer_index > nelements: + propagate_to_next_exiting_surface(neutron) + return none + clone(neutron, tmp_neutron) + # global2local(tmp_neutron, scatterer_index) + interaction = element_interact_path1(threadindex, rng_states, tmp_neutron, scatterer_index) + # local2global(tmp_neutron, scatterer_index) + clone(tmp_neutron, neutron) + if interaction == absorption: + absorb(neutron) + return interaction + # nothing happened + if interaction == none: + # if neutron is not inside, we are done + if locate(neutron)!=inside: + return interaction + # otherwise, interact again + return _interact_path1(threadindex, rng_states, neutron, tmp_neutron) + # interaction must be scatter + clone(neutron, tmp_neutron) + # propagate to surface if necessary + if locate(*neutron[:3]) == inside: + propagate_to_next_exiting_surface(neutron) + # apply attenuation + att = calculate_attenuation(tmp_neutron, neutron[:3]) + neutron[-1] *= att + return scattering + + def _calculate_attenuation(neutron, end, tmp_neutron): + ret = 1.0 + for i in range(nelements): + clone(neutron, tmp_neutron) + # global2local(tmp_neutron, i) + ret *= element_calculate_attenuation(neutron, end, i) + return ret + + def element_interact_path1(threadindex, rng_states, neutron, element_index): + if element_index == 0: + return element0_interact_path1(threadindex, rng_states, neutron) + if element_index == 1: + return element1_interact_path1(threadindex, rng_states, neutron) + if element_index == 2: + return element2_interact_path1(threadindex, rng_states, neutron) + + def element_calculate_attenuation(neutron, end, element_index): + if element_index == 0: + return element0_calculate_attenuation(neutron, end) + if element_index == 1: + return element1_calculate_attenuation(neutron, end) + if element_index == 2: + return element2_calculate_attenuation(neutron, end) + + if test.USE_CUDASIM: + def interact_path1(threadindex, rng_states, neutron): + tmp_neutron = np.zeros(10, dtype=float) + return _interact_path1(threadindex, rng_states, neutron, tmp_neutron) + def calculate_attenuation(neutron, end): + tmp_neutron = np.zeros(10, dtype=float) + return _calculate_attenuation(neutron, end, tmp_neutron) + else: + def interact_path1(threadindex, rng_states, neutron): + tmp_neutron = cuda.local.array(10, dtype=numba.float64) + return _interact_path1(threadindex, rng_states, neutron, tmp_neutron) + def calculate_attenuation(neutron, end): + tmp_neutron = cuda.local.array(10, dtype=numba.float64) + return _calculate_attenuation(neutron, end, tmp_neutron) + return dict( + interact_path1 = interact_path1 + ) diff --git a/tests/components/samples/sampleassemblies/sample+2cylinders/Al.xyz b/tests/components/samples/sampleassemblies/sample+2cylinders/Al.xyz new file mode 100644 index 00000000..95e6d938 --- /dev/null +++ b/tests/components/samples/sampleassemblies/sample+2cylinders/Al.xyz @@ -0,0 +1,6 @@ +4 +4.04 0 0 0 4.04 0 0 0 4.04 +Al 0 0 0 +Al 0.5 0.5 0 +Al 0.5 0 0.5 +Al 0 0.5 0.5 diff --git a/tests/components/samples/sampleassemblies/sample+2cylinders/sample-scatterer.xml b/tests/components/samples/sampleassemblies/sample+2cylinders/sample-scatterer.xml new file mode 100644 index 00000000..47b12cc2 --- /dev/null +++ b/tests/components/samples/sampleassemblies/sample+2cylinders/sample-scatterer.xml @@ -0,0 +1,10 @@ + + + + + + + + + + diff --git a/tests/components/samples/sampleassemblies/sample+2cylinders/sampleassembly.xml b/tests/components/samples/sampleassemblies/sample+2cylinders/sampleassembly.xml new file mode 100644 index 00000000..b4be3c18 --- /dev/null +++ b/tests/components/samples/sampleassemblies/sample+2cylinders/sampleassembly.xml @@ -0,0 +1,47 @@ + + + + + + + + + + + Al + Al.xyz + + + + + + + + + + + + Al + Al.xyz + + + + + + + + + + + + Al + Al.xyz + + + + + + + + + diff --git a/tests/components/samples/sampleassemblies/sample+2cylinders/sampleenv1-scatterer.xml b/tests/components/samples/sampleassemblies/sample+2cylinders/sampleenv1-scatterer.xml new file mode 100644 index 00000000..47b12cc2 --- /dev/null +++ b/tests/components/samples/sampleassemblies/sample+2cylinders/sampleenv1-scatterer.xml @@ -0,0 +1,10 @@ + + + + + + + + + + diff --git a/tests/components/samples/sampleassemblies/sample+2cylinders/sampleenv2-scatterer.xml b/tests/components/samples/sampleassemblies/sample+2cylinders/sampleenv2-scatterer.xml new file mode 100644 index 00000000..47b12cc2 --- /dev/null +++ b/tests/components/samples/sampleassemblies/sample+2cylinders/sampleenv2-scatterer.xml @@ -0,0 +1,10 @@ + + + + + + + + + + diff --git a/tests/scatterers/test_composite_scatterer.py b/tests/scatterers/test_composite_scatterer.py new file mode 100644 index 00000000..1e2ed1f9 --- /dev/null +++ b/tests/scatterers/test_composite_scatterer.py @@ -0,0 +1,35 @@ +import os +thisdir = os.path.dirname(__file__) +import numpy as np +import pytest +from mcvine.acc import test + +@pytest.mark.skipif(not test.USE_CUDASIM, reason='no CUDASIM') +def test(): + path = os.path.join(thisdir, "sampleassemblies", 'sample+2cylinders', 'sampleassembly.xml') + from mcvine.acc.components.samples import loadScattererComposite + composite = loadScattererComposite(path) + from mcvine.acc.scatterers.composite_scatterer import factory_3 + methods = factory_3(composite) + interact_path1 = methods['interact_path1'] + neutron = np.array([0.0,0,0, 0,0,1000, 0,0, 0, 1.]) + threadindex = 0 + rng_states = None + N = 10 #000 + p = 0. + absorbed = 0 + for i in range(N): + neutron = np.array([0.0,0,0, 0,0,1000, 0,0, 0, 1.]) + interact_path1(threadindex, rng_states, neutron) + print(neutron) + p1 = neutron[-1] + if p1>0: + p+=p1 + else: + absorbed += 1 + p/=N + absorbed/=N + print(p, absorbed) + assert p < 0.1 + assert absorbed == 0 + return diff --git a/tests/scatterers/test_scatter_func_factory.py b/tests/scatterers/test_scatter_func_factory.py new file mode 100644 index 00000000..dee61bd9 --- /dev/null +++ b/tests/scatterers/test_scatter_func_factory.py @@ -0,0 +1,14 @@ +import os +thisdir = os.path.dirname(__file__) +import numpy as np +import pytest +from mcvine.acc import test + +@pytest.mark.skipif(not test.USE_CUDASIM, reason='no CUDASIM') +def test(): + path = os.path.join(thisdir, "sampleassemblies", 'sample+2cylinders', 'sampleassembly.xml') + from mcvine.acc.components.samples import loadScattererComposite + composite = loadScattererComposite(path) + from mcvine.acc.scatterers import scatter_func_factory + scatter_func_factory.render(composite) + return From ab984dd65143b86157663117815fecf47b793222 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Sun, 7 May 2023 16:46:56 -0700 Subject: [PATCH 19/53] preliminary implementation of acc.scatterers.ScatterFuncFactory.onCompositeScatterer --- acc/scatterers/__init__.py | 15 +++++---------- acc/scatterers/composite_scatterer.py | 3 ++- 2 files changed, 7 insertions(+), 11 deletions(-) diff --git a/acc/scatterers/__init__.py b/acc/scatterers/__init__.py index 50bed903..5fa00e29 100644 --- a/acc/scatterers/__init__.py +++ b/acc/scatterers/__init__.py @@ -9,17 +9,12 @@ def render(self, scatterer): """ return scatterer.identify(self) - def onCompositeKernel(self, composite): + def onCompositeScatterer(self, composite): elements = composite.elements() - weights = [] - for element in elements: - weights.append(getattr(element, 'weight', 1.)) - continue - s = sum(weights) - for w, element in zip(weights, elements): - element.weight = w/s - from .composite import makeKernelMethods - return makeKernelMethods(composite) + if len(elements) != 3: + raise NotImplementedError + from .composite_scatterer import factory_3 + return factory_3(composite) def onHomogeneousScatterer(self, hs): from .homogeneous_scatterer import factory diff --git a/acc/scatterers/composite_scatterer.py b/acc/scatterers/composite_scatterer.py index ec5be283..840eeef1 100644 --- a/acc/scatterers/composite_scatterer.py +++ b/acc/scatterers/composite_scatterer.py @@ -129,5 +129,6 @@ def calculate_attenuation(neutron, end): tmp_neutron = cuda.local.array(10, dtype=numba.float64) return _calculate_attenuation(neutron, end, tmp_neutron) return dict( - interact_path1 = interact_path1 + interact_path1 = interact_path1, + calculate_attenuation = calculate_attenuation, ) From 3ae877a9be67ce3ab99949cb153ff7c0fe94dc20 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Sun, 7 May 2023 17:20:00 -0700 Subject: [PATCH 20/53] added a real CUDA test of homogenouse_scatterer interact_path1 and fixed a few bugs --- acc/geometry/propagation.py | 1 + acc/scatterers/homogeneous_scatterer.py | 4 ++ .../scatterers/test_homogeneous_scatterer.py | 45 ++++++++++++++++++- 3 files changed, 49 insertions(+), 1 deletion(-) diff --git a/acc/geometry/propagation.py b/acc/geometry/propagation.py index 5de417d5..acd4a199 100644 --- a/acc/geometry/propagation.py +++ b/acc/geometry/propagation.py @@ -191,6 +191,7 @@ def _tof_before_first_exit(neutron, ts): continue return t + @cuda.jit(device=True) def _forward_distance_in_shape(neutron, end, ts): "distance to travel in shape to get to `end`. `end` must be on the path of the neutron" start = neutron[:3] diff --git a/acc/scatterers/homogeneous_scatterer.py b/acc/scatterers/homogeneous_scatterer.py index 0988464a..10f2a8f6 100644 --- a/acc/scatterers/homogeneous_scatterer.py +++ b/acc/scatterers/homogeneous_scatterer.py @@ -44,12 +44,14 @@ def factory(shape, kernel, mcweights, packing_factor): forward_distance_in_shape = propagate_methods["forward_distance_in_shape"] from ..kernels import scatter_func_factory scatter, calc_scattering_coeff, absorb, calc_absorption_coeff = scatter_func_factory.render(kernel) + @cuda.jit(device=True) def calculate_attenuation(neutron, end): length = forward_distance_in_shape(neutron, end) sigma = calc_scattering_coeff(neutron) * packing_factor mu = calc_absorption_coeff(neutron) * packing_factor return exp( - (mu+sigma) * length ) + @cuda.jit(device=True) def _interact_path1(threadindex, rng_states, neutron, tmp_neutron): x, y, z, vx, vy, vz = neutron[:6] loc = locate(x,y,z) @@ -106,10 +108,12 @@ def _interact_path1(threadindex, rng_states, neutron, tmp_neutron): neutron[-1] *= atten2; return scattering if test.USE_CUDASIM: + @cuda.jit(device=True, inline=True) def interact_path1(threadindex, rng_states, neutron): tmp_neutron = np.zeros(10, dtype=float) _interact_path1(threadindex, rng_states, neutron, tmp_neutron) else: + @cuda.jit(device=True, inline=True) def interact_path1(threadindex, rng_states, neutron): tmp_neutron = cuda.local.array(10, dtype=numba.float64) _interact_path1(threadindex, rng_states, neutron, tmp_neutron) diff --git a/tests/scatterers/test_homogeneous_scatterer.py b/tests/scatterers/test_homogeneous_scatterer.py index 547e7d49..9c4dacc5 100644 --- a/tests/scatterers/test_homogeneous_scatterer.py +++ b/tests/scatterers/test_homogeneous_scatterer.py @@ -1,6 +1,7 @@ import os thisdir = os.path.dirname(__file__) -import numpy as np +import math, numpy as np +from numba import cuda import pytest from mcvine.acc import test @@ -43,3 +44,45 @@ def test_scatter_func_factory(): from mcvine.acc.scatterers import scatter_func_factory methods = scatter_func_factory.render(hs) return + +def get_interact_path1(): + path = os.path.join(thisdir, "sampleassemblies", 'isotropic_sphere', 'sampleassembly.xml') + from mcvine.acc.components.samples import loadFirstHomogeneousScatterer + hs = loadFirstHomogeneousScatterer(path) + shape = hs.shape() + kernel = hs.kernel() + mcweights = 1., 1., 1. + packing_factor = 0.6 + from mcvine.acc.scatterers.homogeneous_scatterer import factory + methods = factory(shape, kernel, mcweights, packing_factor) + return methods['interact_path1'] +interact_path1 = get_interact_path1() + +@pytest.mark.skipif(not test.USE_CUDA, reason='no CUDA') +def test_homogeneous_scatterer_cuda(): + @cuda.jit + def interact_path1_kernel(rng_states, neutrons, n_neutrons_per_thread): + N = len(neutrons) + thread_index = cuda.grid(1) + start_index = thread_index*n_neutrons_per_thread + end_index = min(start_index+n_neutrons_per_thread, N) + for i in range(start_index, end_index): + interact_path1(thread_index, rng_states, neutrons[i]) + def run(neutrons): + N = len(neutrons) + threads_per_block = 512 + ntotthreads = int(1e5) + nblocks = math.ceil(ntotthreads / threads_per_block) + actual_nthreads = threads_per_block * nblocks + n_neutrons_per_thread = math.ceil(N / actual_nthreads) + from numba.cuda.random import create_xoroshiro128p_states + from mcvine.acc.config import rng_seed + rng_states = create_xoroshiro128p_states(actual_nthreads, seed=rng_seed) + interact_path1_kernel[nblocks, threads_per_block](rng_states, neutrons, n_neutrons_per_thread) + N = 100 + neutrons = np.zeros((N, 10)) + neutron = np.array([0.0,0,0, 0,0,1000, 0,0, 0, 1.]) + neutrons[:] = neutron + run(neutrons) + return + From 5fa2fe96eefea0d00021bcb3a1a5886512a84351 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Sun, 7 May 2023 17:22:37 -0700 Subject: [PATCH 21/53] minor refactor: tests/scatterers/test_homogeneous_scatterer.py --- tests/scatterers/test_homogeneous_scatterer.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/tests/scatterers/test_homogeneous_scatterer.py b/tests/scatterers/test_homogeneous_scatterer.py index 9c4dacc5..9b55dc3e 100644 --- a/tests/scatterers/test_homogeneous_scatterer.py +++ b/tests/scatterers/test_homogeneous_scatterer.py @@ -45,7 +45,8 @@ def test_scatter_func_factory(): methods = scatter_func_factory.render(hs) return -def get_interact_path1(): +@pytest.mark.skipif(not test.USE_CUDA, reason='no CUDA') +def test_homogeneous_scatterer_cuda(): path = os.path.join(thisdir, "sampleassemblies", 'isotropic_sphere', 'sampleassembly.xml') from mcvine.acc.components.samples import loadFirstHomogeneousScatterer hs = loadFirstHomogeneousScatterer(path) @@ -55,11 +56,7 @@ def get_interact_path1(): packing_factor = 0.6 from mcvine.acc.scatterers.homogeneous_scatterer import factory methods = factory(shape, kernel, mcweights, packing_factor) - return methods['interact_path1'] -interact_path1 = get_interact_path1() - -@pytest.mark.skipif(not test.USE_CUDA, reason='no CUDA') -def test_homogeneous_scatterer_cuda(): + interact_path1 = methods['interact_path1'] @cuda.jit def interact_path1_kernel(rng_states, neutrons, n_neutrons_per_thread): N = len(neutrons) From 72ed97adc79646194b068d8960f72e03a2a9cb4b Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Sun, 7 May 2023 17:23:26 -0700 Subject: [PATCH 22/53] minor (increase ncount to 1e5): tests/scatterers/test_homogeneous_scatterer.py --- tests/scatterers/test_homogeneous_scatterer.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/scatterers/test_homogeneous_scatterer.py b/tests/scatterers/test_homogeneous_scatterer.py index 9b55dc3e..1f3b1589 100644 --- a/tests/scatterers/test_homogeneous_scatterer.py +++ b/tests/scatterers/test_homogeneous_scatterer.py @@ -76,7 +76,7 @@ def run(neutrons): from mcvine.acc.config import rng_seed rng_states = create_xoroshiro128p_states(actual_nthreads, seed=rng_seed) interact_path1_kernel[nblocks, threads_per_block](rng_states, neutrons, n_neutrons_per_thread) - N = 100 + N = int(1e5) neutrons = np.zeros((N, 10)) neutron = np.array([0.0,0,0, 0,0,1000, 0,0, 0, 1.]) neutrons[:] = neutron From de67bc2596250524330c3317add9e17e487ae5b0 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Sun, 7 May 2023 17:29:56 -0700 Subject: [PATCH 23/53] simplify with pytest fixtures: tests/scatterers/test_homogeneous_scatterer.py --- .../scatterers/test_homogeneous_scatterer.py | 41 ++++++++----------- 1 file changed, 17 insertions(+), 24 deletions(-) diff --git a/tests/scatterers/test_homogeneous_scatterer.py b/tests/scatterers/test_homogeneous_scatterer.py index 1f3b1589..39fa5f33 100644 --- a/tests/scatterers/test_homogeneous_scatterer.py +++ b/tests/scatterers/test_homogeneous_scatterer.py @@ -5,18 +5,24 @@ import pytest from mcvine.acc import test -@pytest.mark.skipif(not test.USE_CUDASIM, reason='no CUDASIM') -def test_homogeneous_scatterer(): +@pytest.fixture +def hs_example(): path = os.path.join(thisdir, "sampleassemblies", 'isotropic_sphere', 'sampleassembly.xml') from mcvine.acc.components.samples import loadFirstHomogeneousScatterer - hs = loadFirstHomogeneousScatterer(path) - shape = hs.shape() - kernel = hs.kernel() + return loadFirstHomogeneousScatterer(path) + +@pytest.fixture +def hs_example_interact_path1(hs_example): + shape = hs_example.shape() + kernel = hs_example.kernel() mcweights = 1., 1., 1. packing_factor = 0.6 from mcvine.acc.scatterers.homogeneous_scatterer import factory methods = factory(shape, kernel, mcweights, packing_factor) - interact_path1 = methods['interact_path1'] + return methods['interact_path1'] + +@pytest.mark.skipif(not test.USE_CUDASIM, reason='no CUDASIM') +def test_homogeneous_scatterer(hs_example_interact_path1): threadindex = 0 rng_states = None N = 10000 @@ -24,7 +30,7 @@ def test_homogeneous_scatterer(): absorbed = 0 for i in range(N): neutron = np.array([0.0,0,0, 0,0,1000, 0,0, 0, 1.]) - interact_path1(threadindex, rng_states, neutron) + hs_example_interact_path1(threadindex, rng_states, neutron) p1 = neutron[-1] if p1>0: p+=p1 @@ -37,26 +43,13 @@ def test_homogeneous_scatterer(): return @pytest.mark.skipif(not test.USE_CUDASIM, reason='no CUDASIM') -def test_scatter_func_factory(): - path = os.path.join(thisdir, "sampleassemblies", 'isotropic_sphere', 'sampleassembly.xml') - from mcvine.acc.components.samples import loadFirstHomogeneousScatterer - hs = loadFirstHomogeneousScatterer(path) +def test_scatter_func_factory(hs_example): from mcvine.acc.scatterers import scatter_func_factory - methods = scatter_func_factory.render(hs) + methods = scatter_func_factory.render(hs_example) return @pytest.mark.skipif(not test.USE_CUDA, reason='no CUDA') -def test_homogeneous_scatterer_cuda(): - path = os.path.join(thisdir, "sampleassemblies", 'isotropic_sphere', 'sampleassembly.xml') - from mcvine.acc.components.samples import loadFirstHomogeneousScatterer - hs = loadFirstHomogeneousScatterer(path) - shape = hs.shape() - kernel = hs.kernel() - mcweights = 1., 1., 1. - packing_factor = 0.6 - from mcvine.acc.scatterers.homogeneous_scatterer import factory - methods = factory(shape, kernel, mcweights, packing_factor) - interact_path1 = methods['interact_path1'] +def test_homogeneous_scatterer_cuda(hs_example_interact_path1): @cuda.jit def interact_path1_kernel(rng_states, neutrons, n_neutrons_per_thread): N = len(neutrons) @@ -64,7 +57,7 @@ def interact_path1_kernel(rng_states, neutrons, n_neutrons_per_thread): start_index = thread_index*n_neutrons_per_thread end_index = min(start_index+n_neutrons_per_thread, N) for i in range(start_index, end_index): - interact_path1(thread_index, rng_states, neutrons[i]) + hs_example_interact_path1(thread_index, rng_states, neutrons[i]) def run(neutrons): N = len(neutrons) threads_per_block = 512 From 0beaf0249a3ca4d81a71590e36117fe8db37a677 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Sun, 7 May 2023 20:24:24 -0700 Subject: [PATCH 24/53] composite_scatterer.interact_path1 for three elements passed the sanity check --- acc/geometry/composite.py | 1 + acc/scatterers/composite_scatterer.py | 28 ++++++++++--- tests/scatterers/test_composite_scatterer.py | 43 +++++++++++++++++--- 3 files changed, 61 insertions(+), 11 deletions(-) diff --git a/acc/geometry/composite.py b/acc/geometry/composite.py index 6618ca0e..fb7ed273 100644 --- a/acc/geometry/composite.py +++ b/acc/geometry/composite.py @@ -1,3 +1,4 @@ +import numba from mcvine.acc.geometry.arrow_intersect import max_intersections def make_find_1st_hit(forward_intersect_all, is_onborder, find_shape_containing_point, **kwds): diff --git a/acc/scatterers/composite_scatterer.py b/acc/scatterers/composite_scatterer.py index 840eeef1..c4d16691 100644 --- a/acc/scatterers/composite_scatterer.py +++ b/acc/scatterers/composite_scatterer.py @@ -30,11 +30,11 @@ def factory_3(composite): from ..geometry import arrow_intersect intersect = arrow_intersect.arrow_intersect_func_factory.render(shape) locate = arrow_intersect.locate_func_factory.render(shape) + del shape propagate_methods = makePropagateMethods(intersect, locate) - propagate_to_next_incident_surface = propagate_methods['propagate_to_next_incident_surface'] + del intersect propagate_to_next_exiting_surface = propagate_methods['propagate_to_next_exiting_surface'] - tof_before_first_exit = propagate_methods['tof_before_first_exit'] - forward_distance_in_shape = propagate_methods["forward_distance_in_shape"] + del propagate_methods # elements elements = composite.elements() nelements = len(elements) @@ -43,6 +43,7 @@ def factory_3(composite): shapes = [e.shape() for e in elements] from mcvine.acc.geometry.composite import createMethods_3, make_find_1st_hit methods = createMethods_3(shapes) + del shapes find_1st_hit = make_find_1st_hit(**methods) # methods for element scatterers from . import scatter_func_factory @@ -51,13 +52,16 @@ def factory_3(composite): for e in elements: m = scatter_func_factory.render(e) element_scatter_methods.append(m) + del scatter_func_factory element0_interact_path1 = element_scatter_methods[0]['interact_path1'] element1_interact_path1 = element_scatter_methods[1]['interact_path1'] element2_interact_path1 = element_scatter_methods[2]['interact_path1'] element0_calculate_attenuation = element_scatter_methods[0]['calculate_attenuation'] element1_calculate_attenuation = element_scatter_methods[1]['calculate_attenuation'] element2_calculate_attenuation = element_scatter_methods[2]['calculate_attenuation'] + del elements, element_scatter_methods + @cuda.jit(device=True) def _interact_path1(threadindex, rng_states, neutron, tmp_neutron): x,y,z = neutron[:3] vx,vy,vz = neutron[3:6] @@ -76,20 +80,24 @@ def _interact_path1(threadindex, rng_states, neutron, tmp_neutron): # nothing happened if interaction == none: # if neutron is not inside, we are done - if locate(neutron)!=inside: + x,y,z = neutron[:3] + if locate(x,y,z)!=inside: return interaction # otherwise, interact again - return _interact_path1(threadindex, rng_states, neutron, tmp_neutron) + # return _interact_path1(threadindex, rng_states, neutron, tmp_neutron) + return interaction # interaction must be scatter clone(neutron, tmp_neutron) # propagate to surface if necessary - if locate(*neutron[:3]) == inside: + x,y,z = neutron[:3] + if locate(x,y,z) == inside: propagate_to_next_exiting_surface(neutron) # apply attenuation att = calculate_attenuation(tmp_neutron, neutron[:3]) neutron[-1] *= att return scattering + @cuda.jit(device=True) def _calculate_attenuation(neutron, end, tmp_neutron): ret = 1.0 for i in range(nelements): @@ -98,6 +106,7 @@ def _calculate_attenuation(neutron, end, tmp_neutron): ret *= element_calculate_attenuation(neutron, end, i) return ret + @cuda.jit(device=True) def element_interact_path1(threadindex, rng_states, neutron, element_index): if element_index == 0: return element0_interact_path1(threadindex, rng_states, neutron) @@ -105,7 +114,9 @@ def element_interact_path1(threadindex, rng_states, neutron, element_index): return element1_interact_path1(threadindex, rng_states, neutron) if element_index == 2: return element2_interact_path1(threadindex, rng_states, neutron) + return none + @cuda.jit(device=True) def element_calculate_attenuation(neutron, end, element_index): if element_index == 0: return element0_calculate_attenuation(neutron, end) @@ -113,18 +124,23 @@ def element_calculate_attenuation(neutron, end, element_index): return element1_calculate_attenuation(neutron, end) if element_index == 2: return element2_calculate_attenuation(neutron, end) + return 1.0 if test.USE_CUDASIM: + @cuda.jit(device=True, inline=True) def interact_path1(threadindex, rng_states, neutron): tmp_neutron = np.zeros(10, dtype=float) return _interact_path1(threadindex, rng_states, neutron, tmp_neutron) + @cuda.jit(device=True, inline=True) def calculate_attenuation(neutron, end): tmp_neutron = np.zeros(10, dtype=float) return _calculate_attenuation(neutron, end, tmp_neutron) else: + @cuda.jit(device=True, inline=True) def interact_path1(threadindex, rng_states, neutron): tmp_neutron = cuda.local.array(10, dtype=numba.float64) return _interact_path1(threadindex, rng_states, neutron, tmp_neutron) + @cuda.jit(device=True, inline=True) def calculate_attenuation(neutron, end): tmp_neutron = cuda.local.array(10, dtype=numba.float64) return _calculate_attenuation(neutron, end, tmp_neutron) diff --git a/tests/scatterers/test_composite_scatterer.py b/tests/scatterers/test_composite_scatterer.py index 1e2ed1f9..2f502f99 100644 --- a/tests/scatterers/test_composite_scatterer.py +++ b/tests/scatterers/test_composite_scatterer.py @@ -1,17 +1,21 @@ import os thisdir = os.path.dirname(__file__) -import numpy as np +import math, numpy as np +from numba import cuda import pytest from mcvine.acc import test -@pytest.mark.skipif(not test.USE_CUDASIM, reason='no CUDASIM') -def test(): +@pytest.fixture +def composite_example_interact_path1(): path = os.path.join(thisdir, "sampleassemblies", 'sample+2cylinders', 'sampleassembly.xml') from mcvine.acc.components.samples import loadScattererComposite composite = loadScattererComposite(path) from mcvine.acc.scatterers.composite_scatterer import factory_3 methods = factory_3(composite) - interact_path1 = methods['interact_path1'] + return methods['interact_path1'] + +@pytest.mark.skipif(not test.USE_CUDASIM, reason='no CUDASIM') +def test_cudasim(composite_example_interact_path1): neutron = np.array([0.0,0,0, 0,0,1000, 0,0, 0, 1.]) threadindex = 0 rng_states = None @@ -20,7 +24,7 @@ def test(): absorbed = 0 for i in range(N): neutron = np.array([0.0,0,0, 0,0,1000, 0,0, 0, 1.]) - interact_path1(threadindex, rng_states, neutron) + composite_example_interact_path1(threadindex, rng_states, neutron) print(neutron) p1 = neutron[-1] if p1>0: @@ -33,3 +37,32 @@ def test(): assert p < 0.1 assert absorbed == 0 return + +@pytest.mark.skipif(not test.USE_CUDA, reason='no CUDA') +def test_cuda(composite_example_interact_path1): + @cuda.jit + def interact_path1_kernel(rng_states, neutrons, n_neutrons_per_thread): + N = len(neutrons) + thread_index = cuda.grid(1) + start_index = thread_index*n_neutrons_per_thread + end_index = min(start_index+n_neutrons_per_thread, N) + for i in range(start_index, end_index): + composite_example_interact_path1(thread_index, rng_states, neutrons[i]) + def run(neutrons): + N = len(neutrons) + threads_per_block = 128 + ntotthreads = int(1e5) + nblocks = math.ceil(ntotthreads / threads_per_block) + actual_nthreads = threads_per_block * nblocks + n_neutrons_per_thread = math.ceil(N / actual_nthreads) + from numba.cuda.random import create_xoroshiro128p_states + from mcvine.acc.config import rng_seed + rng_states = create_xoroshiro128p_states(actual_nthreads, seed=rng_seed) + interact_path1_kernel[nblocks, threads_per_block](rng_states, neutrons, n_neutrons_per_thread) + N = int(10) + neutrons = np.zeros((N, 10)) + neutron = np.array([0.0,0,0, 0,0,1000, 0,0, 0, 1.]) + neutrons[:] = neutron + run(neutrons) + print(neutrons) + return From 10cfe20e958ebcbf44e91c7a712a2beefa871646 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Mon, 8 May 2023 08:59:42 -0700 Subject: [PATCH 25/53] refactor: moved methods for 3-element composite from acc.geometry.composite to composite_3 --- acc/geometry/arrow_intersect.py | 2 +- acc/geometry/composite.py | 146 +------------------------ acc/geometry/composite_3.py | 142 ++++++++++++++++++++++++ acc/geometry/locate.py | 3 +- acc/scatterers/composite_scatterer.py | 3 +- tests/geometry/test_composite_shape.py | 4 +- 6 files changed, 151 insertions(+), 149 deletions(-) create mode 100644 acc/geometry/composite_3.py diff --git a/acc/geometry/arrow_intersect.py b/acc/geometry/arrow_intersect.py index 82455ad1..ad5bbfb7 100644 --- a/acc/geometry/arrow_intersect.py +++ b/acc/geometry/arrow_intersect.py @@ -29,7 +29,7 @@ def onUnion(self, u): nelements = len(u.shapes) if nelements == 2: return self.onUnion2(u) if nelements == 3: - from .composite import createMethods_3 + from .composite_3 import createMethods_3 return createMethods_3(u.shapes)['intersect_all'] raise NotImplementedError(f"locate for union of {nelements} elements") diff --git a/acc/geometry/composite.py b/acc/geometry/composite.py index fb7ed273..81d8c687 100644 --- a/acc/geometry/composite.py +++ b/acc/geometry/composite.py @@ -1,5 +1,7 @@ -import numba +import os, numpy as np, numba +from numba import cuda from mcvine.acc.geometry.arrow_intersect import max_intersections +from mcvine.acc import test def make_find_1st_hit(forward_intersect_all, is_onborder, find_shape_containing_point, **kwds): @cuda.jit(device=True) @@ -89,145 +91,3 @@ def makeModule(shapes): ] -# Example code -import os, numpy as np -from numba import cuda -from mcvine.acc._numba import xoroshiro128p_uniform_float32 -from mcvine.acc import test - -def createMethods_3(shapes): - "methods for a list of non-overlapping shapes" - assert len(shapes)==3 - Nshapes = len(shapes) - from mcvine.acc.geometry import arrow_intersect - from mcvine.acc.geometry.location import inside, onborder - from mcvine.acc.geometry.arrow_intersect import max_intersections, insert_into_sorted_list - funcs_list = [] - for shape in shapes: - intersect = arrow_intersect.arrow_intersect_func_factory.render(shape) - locate = arrow_intersect.locate_func_factory.render(shape) - funcs = locate, intersect - funcs_list.append(funcs) - continue - - locate_0, intersect_0 = funcs_list[0] - locate_1, intersect_1 = funcs_list[1] - locate_2, intersect_2 = funcs_list[2] - - @cuda.jit(device=True) - def _intersect_all(x,y,z, vx,vy,vz, ts, ts0, ts1, ts2): - N0 = intersect_0(x,y,z, vx,vy,vz, ts0) - N1 = intersect_1(x,y,z, vx,vy,vz, ts1) - N2 = intersect_2(x,y,z, vx,vy,vz, ts2) - N = 0 - for i in range(N0): - N = insert_into_sorted_list(ts0[i], ts, N) - for i in range(N1): - N = insert_into_sorted_list(ts1[i], ts, N) - for i in range(N2): - N = insert_into_sorted_list(ts2[i], ts, N) - return N - - @cuda.jit(device=True) - def _forward_intersect_all(x,y,z, vx,vy,vz, ts, ts0, ts1, ts2): - N0 = intersect_0(x,y,z, vx,vy,vz, ts0) - N1 = intersect_1(x,y,z, vx,vy,vz, ts1) - N2 = intersect_2(x,y,z, vx,vy,vz, ts2) - N = 0 - for i in range(N0): - t = ts0[i] - if t>0: - N = insert_into_sorted_list(t, ts, N) - for i in range(N1): - t = ts1[i] - if t>0: - N = insert_into_sorted_list(t, ts, N) - for i in range(N2): - t = ts2[i] - if t>0: - N = insert_into_sorted_list(t, ts, N) - return N - - @cuda.jit(device=True) - def find_shape_containing_point(x,y,z): - if locate_0(x,y,z) == inside: - return 0 - if locate_1(x,y,z) == inside: - return 1 - if locate_2(x,y,z) == inside: - return 2 - return -1 - - @cuda.jit(device=True) - def is_onborder(x,y,z): - "check if the point is on the border of any of the shapes" - if locate_0(x,y,z) == onborder: - return True - if locate_1(x,y,z) == onborder: - return True - if locate_2(x,y,z) == onborder: - return True - return False - - if test.USE_CUDASIM: - @cuda.jit(device=True) - def intersect_all(x,y,z, vx,vy,vz, ts): - ts0 = np.zeros(max_intersections, dtype=float) - ts1 = np.zeros(max_intersections, dtype=float) - ts2 = np.zeros(max_intersections, dtype=float) - return _intersect_all(x,y,z, vx,vy,vz, ts, ts0, ts1, ts2) - @cuda.jit(device=True) - def forward_intersect_all(x,y,z, vx,vy,vz, ts): - ts0 = np.zeros(max_intersections, dtype=float) - ts1 = np.zeros(max_intersections, dtype=float) - ts2 = np.zeros(max_intersections, dtype=float) - return _forward_intersect_all(x,y,z, vx,vy,vz, ts, ts0, ts1, ts2) - else: - @cuda.jit(device=True) - def intersect_all(x,y,z, vx,vy,vz, ts): - ts0 = cuda.local.array(max_intersections, dtype=numba.float64) - ts1 = cuda.local.array(max_intersections, dtype=numba.float64) - ts2 = cuda.local.array(max_intersections, dtype=numba.float64) - return _intersect_all(x,y,z, vx,vy,vz, ts, ts0, ts1, ts2) - @cuda.jit(device=True) - def forward_intersect_all(x,y,z, vx,vy,vz, ts): - ts0 = cuda.local.array(max_intersections, dtype=numba.float64) - ts1 = cuda.local.array(max_intersections, dtype=numba.float64) - ts2 = cuda.local.array(max_intersections, dtype=numba.float64) - return _forward_intersect_all(x,y,z, vx,vy,vz, ts, ts0, ts1, ts2) - return dict( - intersect_all = intersect_all, - forward_intersect_all = forward_intersect_all, - find_shape_containing_point = find_shape_containing_point, - is_onborder = is_onborder, - ) - -def createUnionLocateMethod_3(shapes): - assert len(shapes)==3 - Nshapes = len(shapes) - from mcvine.acc.geometry import arrow_intersect - from mcvine.acc.geometry.location import inside, onborder, outside - locates = [] - for shape in shapes: - locate = arrow_intersect.locate_func_factory.render(shape) - locates.append(locate) - continue - - locate_0 = locates[0] - locate_1 = locates[1] - locate_2 = locates[2] - - @cuda.jit(device=True) - def locate(x,y,z): - loc0 = locate_0(x,y,z) - loc1 = locate_1(x,y,z) - loc2 = locate_2(x,y,z) - if loc0 == inside: return inside - if loc1 == inside: return inside - if loc2 == inside: return inside - if loc0 == onborder: return onborder - if loc1 == onborder: return onborder - if loc2 == onborder: return onborder - return outside - - return locate diff --git a/acc/geometry/composite_3.py b/acc/geometry/composite_3.py new file mode 100644 index 00000000..0dd47574 --- /dev/null +++ b/acc/geometry/composite_3.py @@ -0,0 +1,142 @@ +# Example code +import os, numpy as np, numba +from numba import cuda +from mcvine.acc._numba import xoroshiro128p_uniform_float32 +from mcvine.acc import test + +def createMethods_3(shapes): + "methods for a list of non-overlapping shapes" + assert len(shapes)==3 + Nshapes = len(shapes) + from mcvine.acc.geometry import arrow_intersect + from mcvine.acc.geometry.location import inside, onborder + from mcvine.acc.geometry.arrow_intersect import max_intersections, insert_into_sorted_list + funcs_list = [] + for shape in shapes: + intersect = arrow_intersect.arrow_intersect_func_factory.render(shape) + locate = arrow_intersect.locate_func_factory.render(shape) + funcs = locate, intersect + funcs_list.append(funcs) + continue + + locate_0, intersect_0 = funcs_list[0] + locate_1, intersect_1 = funcs_list[1] + locate_2, intersect_2 = funcs_list[2] + + @cuda.jit(device=True) + def _intersect_all(x,y,z, vx,vy,vz, ts, ts0, ts1, ts2): + N0 = intersect_0(x,y,z, vx,vy,vz, ts0) + N1 = intersect_1(x,y,z, vx,vy,vz, ts1) + N2 = intersect_2(x,y,z, vx,vy,vz, ts2) + N = 0 + for i in range(N0): + N = insert_into_sorted_list(ts0[i], ts, N) + for i in range(N1): + N = insert_into_sorted_list(ts1[i], ts, N) + for i in range(N2): + N = insert_into_sorted_list(ts2[i], ts, N) + return N + + @cuda.jit(device=True) + def _forward_intersect_all(x,y,z, vx,vy,vz, ts, ts0, ts1, ts2): + N0 = intersect_0(x,y,z, vx,vy,vz, ts0) + N1 = intersect_1(x,y,z, vx,vy,vz, ts1) + N2 = intersect_2(x,y,z, vx,vy,vz, ts2) + N = 0 + for i in range(N0): + t = ts0[i] + if t>0: + N = insert_into_sorted_list(t, ts, N) + for i in range(N1): + t = ts1[i] + if t>0: + N = insert_into_sorted_list(t, ts, N) + for i in range(N2): + t = ts2[i] + if t>0: + N = insert_into_sorted_list(t, ts, N) + return N + + @cuda.jit(device=True) + def find_shape_containing_point(x,y,z): + if locate_0(x,y,z) == inside: + return 0 + if locate_1(x,y,z) == inside: + return 1 + if locate_2(x,y,z) == inside: + return 2 + return -1 + + @cuda.jit(device=True) + def is_onborder(x,y,z): + "check if the point is on the border of any of the shapes" + if locate_0(x,y,z) == onborder: + return True + if locate_1(x,y,z) == onborder: + return True + if locate_2(x,y,z) == onborder: + return True + return False + + if test.USE_CUDASIM: + @cuda.jit(device=True) + def intersect_all(x,y,z, vx,vy,vz, ts): + ts0 = np.zeros(max_intersections, dtype=float) + ts1 = np.zeros(max_intersections, dtype=float) + ts2 = np.zeros(max_intersections, dtype=float) + return _intersect_all(x,y,z, vx,vy,vz, ts, ts0, ts1, ts2) + @cuda.jit(device=True) + def forward_intersect_all(x,y,z, vx,vy,vz, ts): + ts0 = np.zeros(max_intersections, dtype=float) + ts1 = np.zeros(max_intersections, dtype=float) + ts2 = np.zeros(max_intersections, dtype=float) + return _forward_intersect_all(x,y,z, vx,vy,vz, ts, ts0, ts1, ts2) + else: + @cuda.jit(device=True) + def intersect_all(x,y,z, vx,vy,vz, ts): + ts0 = cuda.local.array(max_intersections, dtype=numba.float64) + ts1 = cuda.local.array(max_intersections, dtype=numba.float64) + ts2 = cuda.local.array(max_intersections, dtype=numba.float64) + return _intersect_all(x,y,z, vx,vy,vz, ts, ts0, ts1, ts2) + @cuda.jit(device=True) + def forward_intersect_all(x,y,z, vx,vy,vz, ts): + ts0 = cuda.local.array(max_intersections, dtype=numba.float64) + ts1 = cuda.local.array(max_intersections, dtype=numba.float64) + ts2 = cuda.local.array(max_intersections, dtype=numba.float64) + return _forward_intersect_all(x,y,z, vx,vy,vz, ts, ts0, ts1, ts2) + return dict( + intersect_all = intersect_all, + forward_intersect_all = forward_intersect_all, + find_shape_containing_point = find_shape_containing_point, + is_onborder = is_onborder, + ) + +def createUnionLocateMethod_3(shapes): + assert len(shapes)==3 + Nshapes = len(shapes) + from mcvine.acc.geometry import arrow_intersect + from mcvine.acc.geometry.location import inside, onborder, outside + locates = [] + for shape in shapes: + locate = arrow_intersect.locate_func_factory.render(shape) + locates.append(locate) + continue + + locate_0 = locates[0] + locate_1 = locates[1] + locate_2 = locates[2] + + @cuda.jit(device=True) + def locate(x,y,z): + loc0 = locate_0(x,y,z) + loc1 = locate_1(x,y,z) + loc2 = locate_2(x,y,z) + if loc0 == inside: return inside + if loc1 == inside: return inside + if loc2 == inside: return inside + if loc0 == onborder: return onborder + if loc1 == onborder: return onborder + if loc2 == onborder: return onborder + return outside + + return locate diff --git a/acc/geometry/locate.py b/acc/geometry/locate.py index c67d3098..979c1779 100644 --- a/acc/geometry/locate.py +++ b/acc/geometry/locate.py @@ -15,7 +15,7 @@ def onUnion(self, u): nelements = len(u.shapes) if nelements == 2: return self.onUnion2(u) if nelements == 3: - from .composite import createUnionLocateMethod_3 + from .composite_3 import createUnionLocateMethod_3 return createUnionLocateMethod_3(u.shapes) raise NotImplementedError(f"locate for union of {nelements} elements") @@ -51,7 +51,6 @@ def onBlock(self, b): @cuda.jit(device=True, inline=True) def locateWrtBlock(x, y, z): return cu_device_locate_wrt_box(x, y, z, W, H, D) - return locateWrtBlock def onTranslation(self, t): diff --git a/acc/scatterers/composite_scatterer.py b/acc/scatterers/composite_scatterer.py index c4d16691..f02c0a0e 100644 --- a/acc/scatterers/composite_scatterer.py +++ b/acc/scatterers/composite_scatterer.py @@ -41,7 +41,8 @@ def factory_3(composite): assert nelements == 3 # find_1st_hit shapes = [e.shape() for e in elements] - from mcvine.acc.geometry.composite import createMethods_3, make_find_1st_hit + from mcvine.acc.geometry.composite import make_find_1st_hit + from mcvine.acc.geometry.composite_3 import createMethods_3 methods = createMethods_3(shapes) del shapes find_1st_hit = make_find_1st_hit(**methods) diff --git a/tests/geometry/test_composite_shape.py b/tests/geometry/test_composite_shape.py index 0dd9213b..5aac1fa3 100644 --- a/tests/geometry/test_composite_shape.py +++ b/tests/geometry/test_composite_shape.py @@ -7,7 +7,8 @@ from instrument.nixml import parse_file from mcvine.acc import test from mcvine.acc.geometry import locate, location, arrow_intersect -from mcvine.acc.geometry.composite import createMethods_3, make_find_1st_hit +from mcvine.acc.geometry.composite_3 import createMethods_3, createUnionLocateMethod_3 +from mcvine.acc.geometry.composite import make_find_1st_hit thisdir = os.path.dirname(__file__) @@ -16,7 +17,6 @@ def test_locate(): union = parse_file(os.path.join(thisdir, 'union_three_elements.xml'))[0] shapes = union.shapes - from mcvine.acc.geometry.composite import createUnionLocateMethod_3 locate_u3 = createUnionLocateMethod_3(shapes) assert locate_u3(0, 0, 0) == locate.inside assert locate_u3(0.025, 0, 0) == locate.onborder From 52fd9509c3050ac792e4305e553b0b721963adb3 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Mon, 8 May 2023 09:12:14 -0700 Subject: [PATCH 26/53] minor refactor: acc/geometry/composite_3.py --- acc/geometry/composite_3.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/acc/geometry/composite_3.py b/acc/geometry/composite_3.py index 0dd47574..512f7a4c 100644 --- a/acc/geometry/composite_3.py +++ b/acc/geometry/composite_3.py @@ -3,14 +3,14 @@ from numba import cuda from mcvine.acc._numba import xoroshiro128p_uniform_float32 from mcvine.acc import test +from mcvine.acc.geometry import arrow_intersect +from mcvine.acc.geometry.location import inside, onborder, outside +from mcvine.acc.geometry.arrow_intersect import max_intersections, insert_into_sorted_list def createMethods_3(shapes): "methods for a list of non-overlapping shapes" assert len(shapes)==3 Nshapes = len(shapes) - from mcvine.acc.geometry import arrow_intersect - from mcvine.acc.geometry.location import inside, onborder - from mcvine.acc.geometry.arrow_intersect import max_intersections, insert_into_sorted_list funcs_list = [] for shape in shapes: intersect = arrow_intersect.arrow_intersect_func_factory.render(shape) @@ -114,8 +114,6 @@ def forward_intersect_all(x,y,z, vx,vy,vz, ts): def createUnionLocateMethod_3(shapes): assert len(shapes)==3 Nshapes = len(shapes) - from mcvine.acc.geometry import arrow_intersect - from mcvine.acc.geometry.location import inside, onborder, outside locates = [] for shape in shapes: locate = arrow_intersect.locate_func_factory.render(shape) From 3bb47aea1853be9f83828ebb42282ae4243c32de Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Mon, 8 May 2023 09:15:10 -0700 Subject: [PATCH 27/53] minor refactor: acc/geometry/composite_3.py --- acc/geometry/composite_3.py | 24 +++++++++++------------- 1 file changed, 11 insertions(+), 13 deletions(-) diff --git a/acc/geometry/composite_3.py b/acc/geometry/composite_3.py index 512f7a4c..86677285 100644 --- a/acc/geometry/composite_3.py +++ b/acc/geometry/composite_3.py @@ -11,13 +11,13 @@ def createMethods_3(shapes): "methods for a list of non-overlapping shapes" assert len(shapes)==3 Nshapes = len(shapes) - funcs_list = [] - for shape in shapes: - intersect = arrow_intersect.arrow_intersect_func_factory.render(shape) - locate = arrow_intersect.locate_func_factory.render(shape) - funcs = locate, intersect - funcs_list.append(funcs) - continue + funcs_list = [ + ( + arrow_intersect.locate_func_factory.render(shape), + arrow_intersect.arrow_intersect_func_factory.render(shape), + ) + for shape in shapes + ] locate_0, intersect_0 = funcs_list[0] locate_1, intersect_1 = funcs_list[1] @@ -114,12 +114,10 @@ def forward_intersect_all(x,y,z, vx,vy,vz, ts): def createUnionLocateMethod_3(shapes): assert len(shapes)==3 Nshapes = len(shapes) - locates = [] - for shape in shapes: - locate = arrow_intersect.locate_func_factory.render(shape) - locates.append(locate) - continue - + locates = [ + arrow_intersect.locate_func_factory.render(shape) + for shape in shapes + ] locate_0 = locates[0] locate_1 = locates[1] locate_2 = locates[2] From d4d7598e5da3e591cbad9315d5295d0ef26246c8 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Mon, 8 May 2023 09:35:22 -0700 Subject: [PATCH 28/53] renamed: tests/geometry/test_composite_shape.py -> test_composite_3_elements.py --- .../{test_composite_shape.py => test_composite_3_elements.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename tests/geometry/{test_composite_shape.py => test_composite_3_elements.py} (100%) diff --git a/tests/geometry/test_composite_shape.py b/tests/geometry/test_composite_3_elements.py similarity index 100% rename from tests/geometry/test_composite_shape.py rename to tests/geometry/test_composite_3_elements.py From 9eed2b1e897a5de0839cf958a99d209ac94de3bb Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Mon, 8 May 2023 10:17:02 -0700 Subject: [PATCH 29/53] simplified loops for intersect_all methods: acc/geometry/composite_3.py --- acc/geometry/composite_3.py | 64 ++++++++++++++++--------------------- 1 file changed, 28 insertions(+), 36 deletions(-) diff --git a/acc/geometry/composite_3.py b/acc/geometry/composite_3.py index 86677285..8e076eb0 100644 --- a/acc/geometry/composite_3.py +++ b/acc/geometry/composite_3.py @@ -24,35 +24,35 @@ def createMethods_3(shapes): locate_2, intersect_2 = funcs_list[2] @cuda.jit(device=True) - def _intersect_all(x,y,z, vx,vy,vz, ts, ts0, ts1, ts2): - N0 = intersect_0(x,y,z, vx,vy,vz, ts0) - N1 = intersect_1(x,y,z, vx,vy,vz, ts1) - N2 = intersect_2(x,y,z, vx,vy,vz, ts2) + def _intersect_all(x,y,z, vx,vy,vz, ts, ts_): N = 0 - for i in range(N0): - N = insert_into_sorted_list(ts0[i], ts, N) - for i in range(N1): - N = insert_into_sorted_list(ts1[i], ts, N) - for i in range(N2): - N = insert_into_sorted_list(ts2[i], ts, N) + N_ = intersect_0(x,y,z, vx,vy,vz, ts_) + for i in range(N_): + N = insert_into_sorted_list(ts_[i], ts, N) + N_ = intersect_1(x,y,z, vx,vy,vz, ts_) + for i in range(N_): + N = insert_into_sorted_list(ts_[i], ts, N) + N_ = intersect_2(x,y,z, vx,vy,vz, ts_) + for i in range(N_): + N = insert_into_sorted_list(ts_[i], ts, N) return N @cuda.jit(device=True) - def _forward_intersect_all(x,y,z, vx,vy,vz, ts, ts0, ts1, ts2): - N0 = intersect_0(x,y,z, vx,vy,vz, ts0) - N1 = intersect_1(x,y,z, vx,vy,vz, ts1) - N2 = intersect_2(x,y,z, vx,vy,vz, ts2) + def _forward_intersect_all(x,y,z, vx,vy,vz, ts, ts_): N = 0 - for i in range(N0): - t = ts0[i] + N_ = intersect_0(x,y,z, vx,vy,vz, ts_) + for i in range(N_): + t = ts_[i] if t>0: N = insert_into_sorted_list(t, ts, N) - for i in range(N1): - t = ts1[i] + N_ = intersect_1(x,y,z, vx,vy,vz, ts_) + for i in range(N_): + t = ts_[i] if t>0: N = insert_into_sorted_list(t, ts, N) - for i in range(N2): - t = ts2[i] + N_ = intersect_2(x,y,z, vx,vy,vz, ts_) + for i in range(N_): + t = ts_[i] if t>0: N = insert_into_sorted_list(t, ts, N) return N @@ -81,29 +81,21 @@ def is_onborder(x,y,z): if test.USE_CUDASIM: @cuda.jit(device=True) def intersect_all(x,y,z, vx,vy,vz, ts): - ts0 = np.zeros(max_intersections, dtype=float) - ts1 = np.zeros(max_intersections, dtype=float) - ts2 = np.zeros(max_intersections, dtype=float) - return _intersect_all(x,y,z, vx,vy,vz, ts, ts0, ts1, ts2) + ts_ = np.zeros(max_intersections, dtype=float) + return _intersect_all(x,y,z, vx,vy,vz, ts, ts_) @cuda.jit(device=True) def forward_intersect_all(x,y,z, vx,vy,vz, ts): - ts0 = np.zeros(max_intersections, dtype=float) - ts1 = np.zeros(max_intersections, dtype=float) - ts2 = np.zeros(max_intersections, dtype=float) - return _forward_intersect_all(x,y,z, vx,vy,vz, ts, ts0, ts1, ts2) + ts_ = np.zeros(max_intersections, dtype=float) + return _forward_intersect_all(x,y,z, vx,vy,vz, ts, ts_) else: @cuda.jit(device=True) def intersect_all(x,y,z, vx,vy,vz, ts): - ts0 = cuda.local.array(max_intersections, dtype=numba.float64) - ts1 = cuda.local.array(max_intersections, dtype=numba.float64) - ts2 = cuda.local.array(max_intersections, dtype=numba.float64) - return _intersect_all(x,y,z, vx,vy,vz, ts, ts0, ts1, ts2) + ts_ = cuda.local.array(max_intersections, dtype=numba.float64) + return _intersect_all(x,y,z, vx,vy,vz, ts, ts_) @cuda.jit(device=True) def forward_intersect_all(x,y,z, vx,vy,vz, ts): - ts0 = cuda.local.array(max_intersections, dtype=numba.float64) - ts1 = cuda.local.array(max_intersections, dtype=numba.float64) - ts2 = cuda.local.array(max_intersections, dtype=numba.float64) - return _forward_intersect_all(x,y,z, vx,vy,vz, ts, ts0, ts1, ts2) + ts_ = cuda.local.array(max_intersections, dtype=numba.float64) + return _forward_intersect_all(x,y,z, vx,vy,vz, ts, ts_) return dict( intersect_all = intersect_all, forward_intersect_all = forward_intersect_all, From e6149da2f559330daee44a6af54bc01738708878 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Mon, 8 May 2023 10:20:52 -0700 Subject: [PATCH 30/53] started working on coder -- added _intersect_all: acc/geometry/composite.py --- acc/geometry/composite.py | 43 +++++++++++++++++++++----- tests/geometry/test_composite_shape.py | 17 ++++++++++ 2 files changed, 52 insertions(+), 8 deletions(-) create mode 100644 tests/geometry/test_composite_shape.py diff --git a/acc/geometry/composite.py b/acc/geometry/composite.py index 81d8c687..eeb5347b 100644 --- a/acc/geometry/composite.py +++ b/acc/geometry/composite.py @@ -2,6 +2,7 @@ from numba import cuda from mcvine.acc.geometry.arrow_intersect import max_intersections from mcvine.acc import test +from .._numba import coder def make_find_1st_hit(forward_intersect_all, is_onborder, find_shape_containing_point, **kwds): @cuda.jit(device=True) @@ -76,18 +77,44 @@ def makeMethods(shapes): m = imp.load_source('mod', mod) return m.createKernelMethods(composite) -def makeModule(shapes): - "make cuda device methods for shapes" +def makeModule(N, overwrite=False): + "make cuda device methods for composite with N elements" from .._numba import coder - nshapes = len(shapes) - modulepath = coder.getModule("shapes", nshapes) - if os.path.exists(modulepath): + modulepath = coder.getModule("composite_shape", N) + if os.path.exists(modulepath) and not overwrite: return modulepath indent = 4*' ' add_indent = lambda lines, n: [indent*n+l for l in lines] - element_scatter_method_defs = [ - f'scatter_{ik}, scattering_coeff_{ik}, absorb_{ik}, absorption_coeff_{ik} = kernel_funcs_list[{ik}]' - for ik in range(nkernels) + _intersect_all = _create__intersect_all_method(N, indent=indent) + lines = _intersect_all + with open(modulepath, 'wt') as ostream: + ostream.write("\n".join(lines)) + +def _create__intersect_all_method(N, indent=4*' '): + header = [ + "@cuda.jit(device=True)", + "def _intersect_all(x,y,z, vx,vy,vz, ts, ts_):", ] + loop = coder.unrollLoop( + N = N, + indent = indent, + before_loop = ["N=0"], + in_loop = [ + "N_ = intersect_{i}(x,y,z, vx,vy,vz, ts_)", + "for i in range(N_):", + indent + "N = insert_into_sorted_list(ts_[i], ts, N)", + ], + after_loop = ["return N"] + ) + return header + loop +module_code_template = """ +import os, numpy as np, numba +from numba import cuda +from mcvine.acc._numba import xoroshiro128p_uniform_float32 +from mcvine.acc import test +from mcvine.acc.geometry import arrow_intersect +from mcvine.acc.geometry.location import inside, onborder, outside +from mcvine.acc.geometry.arrow_intersect import max_intersections, insert_into_sorted_list +""" diff --git a/tests/geometry/test_composite_shape.py b/tests/geometry/test_composite_shape.py new file mode 100644 index 00000000..2e19233b --- /dev/null +++ b/tests/geometry/test_composite_shape.py @@ -0,0 +1,17 @@ +#!/usr/bin/env python + +import pytest, os, time +import math, numpy as np +from numba import cuda + +from instrument.nixml import parse_file +from mcvine.acc import test +from mcvine.acc.geometry import locate, location, arrow_intersect +from mcvine.acc.geometry.composite import makeModule + +thisdir = os.path.dirname(__file__) + +def test_makeModule(): + makeModule(4, overwrite=True) + return + From 5b825ad5c6e9ac91da7a72d9da70d9b3c8c6051b Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Mon, 8 May 2023 17:53:06 -0700 Subject: [PATCH 31/53] implemented method to create `createMethods` for N-element composite: acc/geometry/composite.py --- acc/geometry/composite.py | 111 ++++++++++++++++++++++++++++++++++-- acc/geometry/composite_3.py | 2 - 2 files changed, 107 insertions(+), 6 deletions(-) diff --git a/acc/geometry/composite.py b/acc/geometry/composite.py index eeb5347b..d005911e 100644 --- a/acc/geometry/composite.py +++ b/acc/geometry/composite.py @@ -84,13 +84,64 @@ def makeModule(N, overwrite=False): if os.path.exists(modulepath) and not overwrite: return modulepath indent = 4*' ' - add_indent = lambda lines, n: [indent*n+l for l in lines] - _intersect_all = _create__intersect_all_method(N, indent=indent) - lines = _intersect_all + lines = _create_createMethods(N, indent) with open(modulepath, 'wt') as ostream: ostream.write("\n".join(lines)) -def _create__intersect_all_method(N, indent=4*' '): +def _create_createMethods(N, indent=4*' '): + header = f"""assert len(shapes)=={N} +funcs_list = [ + ( + arrow_intersect.locate_func_factory.render(shape), + arrow_intersect.arrow_intersect_func_factory.render(shape), + ) + for shape in shapes +] +""".splitlines() + funcs_loop = coder.unrollLoop( + N = N, + indent = '', + in_loop = ["locate_{i}, intersect_{i} = funcs_list[{i}]"], + ) + end = """ +if test.USE_CUDASIM: + @cuda.jit(device=True) + def intersect_all(x,y,z, vx,vy,vz, ts): + ts_ = np.zeros(max_intersections, dtype=float) + return _intersect_all(x,y,z, vx,vy,vz, ts, ts_) + @cuda.jit(device=True) + def forward_intersect_all(x,y,z, vx,vy,vz, ts): + ts_ = np.zeros(max_intersections, dtype=float) + return _forward_intersect_all(x,y,z, vx,vy,vz, ts, ts_) +else: + @cuda.jit(device=True) + def intersect_all(x,y,z, vx,vy,vz, ts): + ts_ = cuda.local.array(max_intersections, dtype=numba.float64) + return _intersect_all(x,y,z, vx,vy,vz, ts, ts_) + @cuda.jit(device=True) + def forward_intersect_all(x,y,z, vx,vy,vz, ts): + ts_ = cuda.local.array(max_intersections, dtype=numba.float64) + return _forward_intersect_all(x,y,z, vx,vy,vz, ts, ts_) +return dict( + intersect_all = intersect_all, + forward_intersect_all = forward_intersect_all, + find_shape_containing_point = find_shape_containing_point, + is_onborder = is_onborder, +) + """.splitlines() + body = ( + header + funcs_loop + + _create__intersect_all(N, indent) + + _create__forward_intersect_all(N, indent) + + _create_find_shape_containing_point(N, indent) + + _create_is_onborder(N, indent) + + end + ) + add_indent = lambda lines, n: [indent*n+l for l in lines] + return ["def createMethods(shapes):"] + add_indent(body, 1) + + +def _create__intersect_all(N, indent=4*' '): header = [ "@cuda.jit(device=True)", "def _intersect_all(x,y,z, vx,vy,vz, ts, ts_):", @@ -108,6 +159,58 @@ def _create__intersect_all_method(N, indent=4*' '): ) return header + loop +def _create__forward_intersect_all(N, indent=4*' '): + header = [ + "@cuda.jit(device=True)", + "def _forward_intersect_all(x,y,z, vx,vy,vz, ts, ts_):", + ] + loop = coder.unrollLoop( + N = N, + indent = indent, + before_loop = ["N=0"], + in_loop = [ + "N_ = intersect_{i}(x,y,z, vx,vy,vz, ts_)", + "for i in range(N_):", + indent + "t = ts_[i]", + indent + "if t>0:", + indent*2 + "N = insert_into_sorted_list(t, ts, N)", + ], + after_loop = ["return N"] + ) + return header + loop + +def _create_find_shape_containing_point(N, indent=4*' '): + header = [ + "@cuda.jit(device=True)", + "def find_shape_containing_point(x,y,z):", + ] + loop = coder.unrollLoop( + N = N, + indent = indent, + in_loop = [ + "if locate_{i}(x,y,z) == inside:", + indent + "return {i}", + ], + after_loop = ["return -1"] + ) + return header + loop + +def _create_is_onborder(N, indent=4*' '): + header = [ + "@cuda.jit(device=True)", + "def is_onborder(x,y,z):", + ] + loop = coder.unrollLoop( + N = N, + indent = indent, + in_loop = [ + "if locate_{i}(x,y,z) == onborder:", + indent + "return True", + ], + after_loop = ["return False"] + ) + return header + loop + module_code_template = """ import os, numpy as np, numba from numba import cuda diff --git a/acc/geometry/composite_3.py b/acc/geometry/composite_3.py index 8e076eb0..06f8fe23 100644 --- a/acc/geometry/composite_3.py +++ b/acc/geometry/composite_3.py @@ -10,7 +10,6 @@ def createMethods_3(shapes): "methods for a list of non-overlapping shapes" assert len(shapes)==3 - Nshapes = len(shapes) funcs_list = [ ( arrow_intersect.locate_func_factory.render(shape), @@ -105,7 +104,6 @@ def forward_intersect_all(x,y,z, vx,vy,vz, ts): def createUnionLocateMethod_3(shapes): assert len(shapes)==3 - Nshapes = len(shapes) locates = [ arrow_intersect.locate_func_factory.render(shape) for shape in shapes From 9880dcb0c93d25d999cb652f11c3064568e52f9c Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Mon, 8 May 2023 20:14:16 -0700 Subject: [PATCH 32/53] trivial cosmetic: acc/geometry/composite.py --- acc/geometry/composite.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/acc/geometry/composite.py b/acc/geometry/composite.py index d005911e..684a4680 100644 --- a/acc/geometry/composite.py +++ b/acc/geometry/composite.py @@ -130,7 +130,7 @@ def forward_intersect_all(x,y,z, vx,vy,vz, ts): ) """.splitlines() body = ( - header + funcs_loop + header + funcs_loop + [''] + _create__intersect_all(N, indent) + _create__forward_intersect_all(N, indent) + _create_find_shape_containing_point(N, indent) From ee8ae4e78abfb04ff1941429e9b3d378ad01a0c2 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Mon, 8 May 2023 21:38:28 -0700 Subject: [PATCH 33/53] method factories are classes now: acc/geometry/composite.py --- acc/geometry/composite.py | 259 ++++++++++++++++++++++++-------------- 1 file changed, 167 insertions(+), 92 deletions(-) diff --git a/acc/geometry/composite.py b/acc/geometry/composite.py index 684a4680..e328bcab 100644 --- a/acc/geometry/composite.py +++ b/acc/geometry/composite.py @@ -84,12 +84,83 @@ def makeModule(N, overwrite=False): if os.path.exists(modulepath) and not overwrite: return modulepath indent = 4*' ' - lines = _create_createMethods(N, indent) + imports = """import os, numpy as np, numba +from numba import cuda +from mcvine.acc._numba import xoroshiro128p_uniform_float32 +from mcvine.acc import test +from mcvine.acc.geometry import arrow_intersect +from mcvine.acc.geometry.location import inside, onborder, outside +from mcvine.acc.geometry.arrow_intersect import max_intersections, insert_into_sorted_list +""".splitlines() + createMethods = Coder_createMethods(N, indent)() + createUnionLocateMethod = Coder_createUnionLocateMethod(N, indent)() + lines = imports + [''] + createMethods + createUnionLocateMethod with open(modulepath, 'wt') as ostream: ostream.write("\n".join(lines)) -def _create_createMethods(N, indent=4*' '): - header = f"""assert len(shapes)=={N} +class Coder_createUnionLocateMethod: + + def __init__(self, N, indent=4*' '): + self.N = N + self.indent = indent + return + + def __call__(self): + N, indent = self.N, self.indent + header = f"""assert len(shapes)=={N} +locates = [ + arrow_intersect.locate_func_factory.render(shape) + for shape in shapes +] +""".splitlines() + locate_loop = coder.unrollLoop( + N = N, + indent = '', + in_loop = ["locate_{i} = locates[{i}]"], + ) + body = header + locate_loop + [''] + self.locate() + add_indent = lambda lines, n: [indent*n+l for l in lines] + return ( + ["def createUnionLocateMethod(shapes):"] + + add_indent(body, 1) + + [indent+'return locate'] + ) + + def locate(self): + N, indent = self.N, self.indent + header = [ + "@cuda.jit(device=True)", + "def locate(x,y,z):", + ] + assign_loop = coder.unrollLoop( + N = N, + indent = indent, + in_loop = ["loc{i} = locate_{i}(x,y,z)"], + ) + inside_loop = coder.unrollLoop( + N = N, + indent = indent, + in_loop = ["if loc{i} == inside: return inside"], + ) + onborder_loop = coder.unrollLoop( + N = N, + indent = indent, + in_loop = ["if loc{i} == onborder: return onborder"], + after_loop = ['return outside'] + ) + return header + assign_loop + inside_loop + onborder_loop + + +class Coder_createMethods: + + def __init__(self, N, indent=4*' '): + self.N = N + self.indent = indent + return + + def __call__(self): + N, indent = self.N, self.indent + header = f"""assert len(shapes)=={N} funcs_list = [ ( arrow_intersect.locate_func_factory.render(shape), @@ -97,13 +168,13 @@ def _create_createMethods(N, indent=4*' '): ) for shape in shapes ] -""".splitlines() - funcs_loop = coder.unrollLoop( - N = N, - indent = '', - in_loop = ["locate_{i}, intersect_{i} = funcs_list[{i}]"], - ) - end = """ + """.splitlines() + funcs_loop = coder.unrollLoop( + N = N, + indent = '', + in_loop = ["locate_{i}, intersect_{i} = funcs_list[{i}]"], + ) + end = """ if test.USE_CUDASIM: @cuda.jit(device=True) def intersect_all(x,y,z, vx,vy,vz, ts): @@ -128,88 +199,92 @@ def forward_intersect_all(x,y,z, vx,vy,vz, ts): find_shape_containing_point = find_shape_containing_point, is_onborder = is_onborder, ) - """.splitlines() - body = ( - header + funcs_loop + [''] - + _create__intersect_all(N, indent) - + _create__forward_intersect_all(N, indent) - + _create_find_shape_containing_point(N, indent) - + _create_is_onborder(N, indent) - + end - ) - add_indent = lambda lines, n: [indent*n+l for l in lines] - return ["def createMethods(shapes):"] + add_indent(body, 1) - - -def _create__intersect_all(N, indent=4*' '): - header = [ - "@cuda.jit(device=True)", - "def _intersect_all(x,y,z, vx,vy,vz, ts, ts_):", - ] - loop = coder.unrollLoop( - N = N, - indent = indent, - before_loop = ["N=0"], - in_loop = [ - "N_ = intersect_{i}(x,y,z, vx,vy,vz, ts_)", - "for i in range(N_):", - indent + "N = insert_into_sorted_list(ts_[i], ts, N)", - ], - after_loop = ["return N"] - ) - return header + loop - -def _create__forward_intersect_all(N, indent=4*' '): - header = [ - "@cuda.jit(device=True)", - "def _forward_intersect_all(x,y,z, vx,vy,vz, ts, ts_):", - ] - loop = coder.unrollLoop( - N = N, - indent = indent, - before_loop = ["N=0"], - in_loop = [ - "N_ = intersect_{i}(x,y,z, vx,vy,vz, ts_)", - "for i in range(N_):", - indent + "t = ts_[i]", - indent + "if t>0:", - indent*2 + "N = insert_into_sorted_list(t, ts, N)", - ], - after_loop = ["return N"] - ) - return header + loop - -def _create_find_shape_containing_point(N, indent=4*' '): - header = [ - "@cuda.jit(device=True)", - "def find_shape_containing_point(x,y,z):", - ] - loop = coder.unrollLoop( - N = N, - indent = indent, - in_loop = [ - "if locate_{i}(x,y,z) == inside:", - indent + "return {i}", - ], - after_loop = ["return -1"] - ) - return header + loop - -def _create_is_onborder(N, indent=4*' '): - header = [ - "@cuda.jit(device=True)", - "def is_onborder(x,y,z):", - ] - loop = coder.unrollLoop( - N = N, - indent = indent, - in_loop = [ - "if locate_{i}(x,y,z) == onborder:", - indent + "return True", - ], - after_loop = ["return False"] - ) - return header + loop + """.splitlines() + body = ( + header + funcs_loop + [''] + + self._intersect_all() + + self._forward_intersect_all() + + self.find_shape_containing_point() + + self.is_onborder() + + end + ) + add_indent = lambda lines, n: [indent*n+l for l in lines] + return ["def createMethods(shapes):"] + add_indent(body, 1) + + + def _intersect_all(self): + N, indent = self.N, self.indent + header = [ + "@cuda.jit(device=True)", + "def _intersect_all(x,y,z, vx,vy,vz, ts, ts_):", + ] + loop = coder.unrollLoop( + N = N, + indent = indent, + before_loop = ["N=0"], + in_loop = [ + "N_ = intersect_{i}(x,y,z, vx,vy,vz, ts_)", + "for i in range(N_):", + indent + "N = insert_into_sorted_list(ts_[i], ts, N)", + ], + after_loop = ["return N"] + ) + return header + loop + + def _forward_intersect_all(self): + N, indent = self.N, self.indent + header = [ + "@cuda.jit(device=True)", + "def _forward_intersect_all(x,y,z, vx,vy,vz, ts, ts_):", + ] + loop = coder.unrollLoop( + N = N, + indent = indent, + before_loop = ["N=0"], + in_loop = [ + "N_ = intersect_{i}(x,y,z, vx,vy,vz, ts_)", + "for i in range(N_):", + indent + "t = ts_[i]", + indent + "if t>0:", + indent*2 + "N = insert_into_sorted_list(t, ts, N)", + ], + after_loop = ["return N"] + ) + return header + loop + + def find_shape_containing_point(self): + N, indent = self.N, self.indent + header = [ + "@cuda.jit(device=True)", + "def find_shape_containing_point(x,y,z):", + ] + loop = coder.unrollLoop( + N = N, + indent = indent, + in_loop = [ + "if locate_{i}(x,y,z) == inside:", + indent + "return {i}", + ], + after_loop = ["return -1"] + ) + return header + loop + + def is_onborder(self): + N, indent = self.N, self.indent + header = [ + "@cuda.jit(device=True)", + "def is_onborder(x,y,z):", + ] + loop = coder.unrollLoop( + N = N, + indent = indent, + in_loop = [ + "if locate_{i}(x,y,z) == onborder:", + indent + "return True", + ], + after_loop = ["return False"] + ) + return header + loop module_code_template = """ import os, numpy as np, numba From 38feaa7d613c7d267670ee7a915209776b347189 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Mon, 8 May 2023 22:04:36 -0700 Subject: [PATCH 34/53] test the module made by the composite-shape ray-tracing method factory --- acc/geometry/composite.py | 8 ++-- tests/geometry/test_composite_shape.py | 53 +++++++++++++++++++++++++- 2 files changed, 55 insertions(+), 6 deletions(-) diff --git a/acc/geometry/composite.py b/acc/geometry/composite.py index e328bcab..39442419 100644 --- a/acc/geometry/composite.py +++ b/acc/geometry/composite.py @@ -71,11 +71,10 @@ def find_1st_hit(x,y,z, vx,vy,vz): return find_1st_hit -def makeMethods(shapes): - mod = makeModule(shapes) +def importModule(N): + mod = makeModule(N) import imp - m = imp.load_source('mod', mod) - return m.createKernelMethods(composite) + return imp.load_source('mod', mod) def makeModule(N, overwrite=False): "make cuda device methods for composite with N elements" @@ -97,6 +96,7 @@ def makeModule(N, overwrite=False): lines = imports + [''] + createMethods + createUnionLocateMethod with open(modulepath, 'wt') as ostream: ostream.write("\n".join(lines)) + return modulepath class Coder_createUnionLocateMethod: diff --git a/tests/geometry/test_composite_shape.py b/tests/geometry/test_composite_shape.py index 2e19233b..8596739a 100644 --- a/tests/geometry/test_composite_shape.py +++ b/tests/geometry/test_composite_shape.py @@ -6,8 +6,8 @@ from instrument.nixml import parse_file from mcvine.acc import test -from mcvine.acc.geometry import locate, location, arrow_intersect -from mcvine.acc.geometry.composite import makeModule +from mcvine.acc.geometry import locate +from mcvine.acc.geometry.composite import makeModule, importModule, make_find_1st_hit thisdir = os.path.dirname(__file__) @@ -15,3 +15,52 @@ def test_makeModule(): makeModule(4, overwrite=True) return +@pytest.fixture +def ray_tracing_methods(): + union = parse_file(os.path.join(thisdir, 'union_three_elements.xml'))[0] + shapes = union.shapes + mod = importModule(len(shapes)) + methods = mod.createMethods(shapes) + methods['union_locate'] = mod.createUnionLocateMethod(shapes) + return methods + +# device functions can be tested with CUDASIM only +@pytest.mark.skipif(not test.USE_CUDASIM, reason='no CUDASIM') +def test_locate(ray_tracing_methods): + locate_u3 = ray_tracing_methods['union_locate'] + assert locate_u3(0, 0, 0) == locate.inside + assert locate_u3(0.025, 0, 0) == locate.onborder + assert locate_u3(0.099, 0, 0) == locate.onborder + assert locate_u3(0.1, 0, 0) == locate.onborder + assert locate_u3(0.199, 0, 0) == locate.onborder + assert locate_u3(0.2, 0, 0) == locate.onborder + assert locate_u3(0.3, 0, 0) == locate.outside + return + +@pytest.mark.skipif(not test.USE_CUDASIM, reason='no CUDASIM') +def test_intersect_all(ray_tracing_methods): + intersect_all = ray_tracing_methods['intersect_all'] + ts = np.zeros(11) + assert intersect_all(0.,0.,0., 1.,0.,0., ts) == 10 + expected = np.array([-0.2,-0.199, -0.1,-0.099, -0.025, 0.025, 0.099,0.1, 0.199,0.2, 0]) + assert np.allclose(ts, expected) + return + +@pytest.mark.skipif(not test.USE_CUDASIM, reason='no CUDASIM') +def test_forward_intersect_all(ray_tracing_methods): + forward_intersect_all = ray_tracing_methods['forward_intersect_all'] + ts = np.zeros(6) + assert forward_intersect_all(0.,0.,0., 1.,0.,0., ts) == 5 + expected = np.array([0.025,0.099,0.1, 0.199,0.2, 0]) + assert np.allclose(ts, expected) + return + +@pytest.mark.skipif(not test.USE_CUDASIM, reason='no CUDASIM') +def test_find_1st_hit(ray_tracing_methods): + find_1st_hit = make_find_1st_hit(**ray_tracing_methods) + assert find_1st_hit(-0.25,0.,0., 1.,0.,0.) == 0 + assert find_1st_hit(-0.15,0.,0., 1.,0.,0.) == 1 + assert find_1st_hit(-0.05,0.,0., 1.,0.,0.) == 2 + assert find_1st_hit(0.05,0.,0., 1.,0.,0.) == 1 + assert find_1st_hit(0.15,0.,0., 1.,0.,0.) == 0 + assert find_1st_hit(0.25,0.,0., 1.,0.,0.) == -1 From 05ad8a0df03d54b6ea54828c811cd9301d5baac6 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Tue, 9 May 2023 07:48:49 -0700 Subject: [PATCH 35/53] added a real CUDA test: tests/geometry/test_composite_shape.py --- tests/geometry/test_composite_shape.py | 30 ++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/tests/geometry/test_composite_shape.py b/tests/geometry/test_composite_shape.py index 8596739a..729d7795 100644 --- a/tests/geometry/test_composite_shape.py +++ b/tests/geometry/test_composite_shape.py @@ -64,3 +64,33 @@ def test_find_1st_hit(ray_tracing_methods): assert find_1st_hit(0.05,0.,0., 1.,0.,0.) == 1 assert find_1st_hit(0.15,0.,0., 1.,0.,0.) == 0 assert find_1st_hit(0.25,0.,0., 1.,0.,0.) == -1 + +@pytest.mark.skipif(not test.USE_CUDA, reason='no CUDA') +def test_cuda(ray_tracing_methods): + find_1st_hit = make_find_1st_hit(**ray_tracing_methods) + @cuda.jit + def find_1st_hit_kernel(neutrons, out, n_neutrons_per_thread): + N = len(neutrons) + thread_index = cuda.grid(1) + start_index = thread_index*n_neutrons_per_thread + end_index = min(start_index+n_neutrons_per_thread, N) + for i in range(start_index, end_index): + neutron = neutrons[i] + x,y,z, vx,vy,vz = neutron[:6] + out[i] = find_1st_hit(x,y,z, vx,vy,vz) + def run(neutrons, first_hit): + N = len(neutrons) + threads_per_block = 128 + ntotthreads = int(1e5) + nblocks = math.ceil(ntotthreads / threads_per_block) + actual_nthreads = threads_per_block * nblocks + n_neutrons_per_thread = math.ceil(N / actual_nthreads) + find_1st_hit_kernel[nblocks, threads_per_block](neutrons, first_hit, n_neutrons_per_thread) + N = int(100) + neutrons = np.zeros((N, 10)) + neutrons[:, :3] = np.random.random((N, 3)) * 0.1 + neutrons[:, 3:6] = np.random.random((N, 3)) * 1000 + first_hit = np.zeros(N, dtype=int) + run(neutrons, first_hit) + print(first_hit) + return From 6e687ccd2d8dffd534a8b5e359208f3f82752640 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Wed, 10 May 2023 22:15:28 -0700 Subject: [PATCH 36/53] renamed acc.scatterers.composite_scatterer to ..composite_3 --- acc/scatterers/__init__.py | 2 +- acc/scatterers/composite_3.py | 151 +++++++++++++++++++ acc/scatterers/composite_scatterer.py | 126 ---------------- tests/scatterers/test_composite_scatterer.py | 2 +- 4 files changed, 153 insertions(+), 128 deletions(-) create mode 100644 acc/scatterers/composite_3.py diff --git a/acc/scatterers/__init__.py b/acc/scatterers/__init__.py index 5fa00e29..208a41af 100644 --- a/acc/scatterers/__init__.py +++ b/acc/scatterers/__init__.py @@ -13,7 +13,7 @@ def onCompositeScatterer(self, composite): elements = composite.elements() if len(elements) != 3: raise NotImplementedError - from .composite_scatterer import factory_3 + from .composite_3 import factory_3 return factory_3(composite) def onHomogeneousScatterer(self, hs): diff --git a/acc/scatterers/composite_3.py b/acc/scatterers/composite_3.py new file mode 100644 index 00000000..f02c0a0e --- /dev/null +++ b/acc/scatterers/composite_3.py @@ -0,0 +1,151 @@ +# -*- python -*- +# +# Jiao Lin +# + +import numpy as np, numba +from numba import cuda, void, int64 +from mcvine.acc._numba import xoroshiro128p_uniform_float32 +from math import sqrt, exp + +from .interaction_types import absorption, scattering, none +from .. import test +from ..neutron import absorb, prop_dt_inplace, clone +from ..geometry.arrow_intersect import max_intersections +from ..geometry.locate import inside, outside, onborder +from ..geometry.propagation import makePropagateMethods + + +from numba.core import config +if not config.ENABLE_CUDASIM: + from numba.cuda.compiler import Dispatcher, DeviceFunction + +from ..config import get_numba_floattype +NB_FLOAT = get_numba_floattype() + + +def factory_3(composite): + # methods to deal with the whole shape of the composite + shape = composite.shape() + from ..geometry import arrow_intersect + intersect = arrow_intersect.arrow_intersect_func_factory.render(shape) + locate = arrow_intersect.locate_func_factory.render(shape) + del shape + propagate_methods = makePropagateMethods(intersect, locate) + del intersect + propagate_to_next_exiting_surface = propagate_methods['propagate_to_next_exiting_surface'] + del propagate_methods + # elements + elements = composite.elements() + nelements = len(elements) + assert nelements == 3 + # find_1st_hit + shapes = [e.shape() for e in elements] + from mcvine.acc.geometry.composite import make_find_1st_hit + from mcvine.acc.geometry.composite_3 import createMethods_3 + methods = createMethods_3(shapes) + del shapes + find_1st_hit = make_find_1st_hit(**methods) + # methods for element scatterers + from . import scatter_func_factory + #element_scatter_methods = [scatter_func_factory.render(e) for e in elements] + element_scatter_methods = [] + for e in elements: + m = scatter_func_factory.render(e) + element_scatter_methods.append(m) + del scatter_func_factory + element0_interact_path1 = element_scatter_methods[0]['interact_path1'] + element1_interact_path1 = element_scatter_methods[1]['interact_path1'] + element2_interact_path1 = element_scatter_methods[2]['interact_path1'] + element0_calculate_attenuation = element_scatter_methods[0]['calculate_attenuation'] + element1_calculate_attenuation = element_scatter_methods[1]['calculate_attenuation'] + element2_calculate_attenuation = element_scatter_methods[2]['calculate_attenuation'] + del elements, element_scatter_methods + + @cuda.jit(device=True) + def _interact_path1(threadindex, rng_states, neutron, tmp_neutron): + x,y,z = neutron[:3] + vx,vy,vz = neutron[3:6] + scatterer_index = find_1st_hit(x,y,z, vx,vy,vz) + if scatterer_index < 0 or scatterer_index > nelements: + propagate_to_next_exiting_surface(neutron) + return none + clone(neutron, tmp_neutron) + # global2local(tmp_neutron, scatterer_index) + interaction = element_interact_path1(threadindex, rng_states, tmp_neutron, scatterer_index) + # local2global(tmp_neutron, scatterer_index) + clone(tmp_neutron, neutron) + if interaction == absorption: + absorb(neutron) + return interaction + # nothing happened + if interaction == none: + # if neutron is not inside, we are done + x,y,z = neutron[:3] + if locate(x,y,z)!=inside: + return interaction + # otherwise, interact again + # return _interact_path1(threadindex, rng_states, neutron, tmp_neutron) + return interaction + # interaction must be scatter + clone(neutron, tmp_neutron) + # propagate to surface if necessary + x,y,z = neutron[:3] + if locate(x,y,z) == inside: + propagate_to_next_exiting_surface(neutron) + # apply attenuation + att = calculate_attenuation(tmp_neutron, neutron[:3]) + neutron[-1] *= att + return scattering + + @cuda.jit(device=True) + def _calculate_attenuation(neutron, end, tmp_neutron): + ret = 1.0 + for i in range(nelements): + clone(neutron, tmp_neutron) + # global2local(tmp_neutron, i) + ret *= element_calculate_attenuation(neutron, end, i) + return ret + + @cuda.jit(device=True) + def element_interact_path1(threadindex, rng_states, neutron, element_index): + if element_index == 0: + return element0_interact_path1(threadindex, rng_states, neutron) + if element_index == 1: + return element1_interact_path1(threadindex, rng_states, neutron) + if element_index == 2: + return element2_interact_path1(threadindex, rng_states, neutron) + return none + + @cuda.jit(device=True) + def element_calculate_attenuation(neutron, end, element_index): + if element_index == 0: + return element0_calculate_attenuation(neutron, end) + if element_index == 1: + return element1_calculate_attenuation(neutron, end) + if element_index == 2: + return element2_calculate_attenuation(neutron, end) + return 1.0 + + if test.USE_CUDASIM: + @cuda.jit(device=True, inline=True) + def interact_path1(threadindex, rng_states, neutron): + tmp_neutron = np.zeros(10, dtype=float) + return _interact_path1(threadindex, rng_states, neutron, tmp_neutron) + @cuda.jit(device=True, inline=True) + def calculate_attenuation(neutron, end): + tmp_neutron = np.zeros(10, dtype=float) + return _calculate_attenuation(neutron, end, tmp_neutron) + else: + @cuda.jit(device=True, inline=True) + def interact_path1(threadindex, rng_states, neutron): + tmp_neutron = cuda.local.array(10, dtype=numba.float64) + return _interact_path1(threadindex, rng_states, neutron, tmp_neutron) + @cuda.jit(device=True, inline=True) + def calculate_attenuation(neutron, end): + tmp_neutron = cuda.local.array(10, dtype=numba.float64) + return _calculate_attenuation(neutron, end, tmp_neutron) + return dict( + interact_path1 = interact_path1, + calculate_attenuation = calculate_attenuation, + ) diff --git a/acc/scatterers/composite_scatterer.py b/acc/scatterers/composite_scatterer.py index f02c0a0e..3c7244cc 100644 --- a/acc/scatterers/composite_scatterer.py +++ b/acc/scatterers/composite_scatterer.py @@ -23,129 +23,3 @@ from ..config import get_numba_floattype NB_FLOAT = get_numba_floattype() - -def factory_3(composite): - # methods to deal with the whole shape of the composite - shape = composite.shape() - from ..geometry import arrow_intersect - intersect = arrow_intersect.arrow_intersect_func_factory.render(shape) - locate = arrow_intersect.locate_func_factory.render(shape) - del shape - propagate_methods = makePropagateMethods(intersect, locate) - del intersect - propagate_to_next_exiting_surface = propagate_methods['propagate_to_next_exiting_surface'] - del propagate_methods - # elements - elements = composite.elements() - nelements = len(elements) - assert nelements == 3 - # find_1st_hit - shapes = [e.shape() for e in elements] - from mcvine.acc.geometry.composite import make_find_1st_hit - from mcvine.acc.geometry.composite_3 import createMethods_3 - methods = createMethods_3(shapes) - del shapes - find_1st_hit = make_find_1st_hit(**methods) - # methods for element scatterers - from . import scatter_func_factory - #element_scatter_methods = [scatter_func_factory.render(e) for e in elements] - element_scatter_methods = [] - for e in elements: - m = scatter_func_factory.render(e) - element_scatter_methods.append(m) - del scatter_func_factory - element0_interact_path1 = element_scatter_methods[0]['interact_path1'] - element1_interact_path1 = element_scatter_methods[1]['interact_path1'] - element2_interact_path1 = element_scatter_methods[2]['interact_path1'] - element0_calculate_attenuation = element_scatter_methods[0]['calculate_attenuation'] - element1_calculate_attenuation = element_scatter_methods[1]['calculate_attenuation'] - element2_calculate_attenuation = element_scatter_methods[2]['calculate_attenuation'] - del elements, element_scatter_methods - - @cuda.jit(device=True) - def _interact_path1(threadindex, rng_states, neutron, tmp_neutron): - x,y,z = neutron[:3] - vx,vy,vz = neutron[3:6] - scatterer_index = find_1st_hit(x,y,z, vx,vy,vz) - if scatterer_index < 0 or scatterer_index > nelements: - propagate_to_next_exiting_surface(neutron) - return none - clone(neutron, tmp_neutron) - # global2local(tmp_neutron, scatterer_index) - interaction = element_interact_path1(threadindex, rng_states, tmp_neutron, scatterer_index) - # local2global(tmp_neutron, scatterer_index) - clone(tmp_neutron, neutron) - if interaction == absorption: - absorb(neutron) - return interaction - # nothing happened - if interaction == none: - # if neutron is not inside, we are done - x,y,z = neutron[:3] - if locate(x,y,z)!=inside: - return interaction - # otherwise, interact again - # return _interact_path1(threadindex, rng_states, neutron, tmp_neutron) - return interaction - # interaction must be scatter - clone(neutron, tmp_neutron) - # propagate to surface if necessary - x,y,z = neutron[:3] - if locate(x,y,z) == inside: - propagate_to_next_exiting_surface(neutron) - # apply attenuation - att = calculate_attenuation(tmp_neutron, neutron[:3]) - neutron[-1] *= att - return scattering - - @cuda.jit(device=True) - def _calculate_attenuation(neutron, end, tmp_neutron): - ret = 1.0 - for i in range(nelements): - clone(neutron, tmp_neutron) - # global2local(tmp_neutron, i) - ret *= element_calculate_attenuation(neutron, end, i) - return ret - - @cuda.jit(device=True) - def element_interact_path1(threadindex, rng_states, neutron, element_index): - if element_index == 0: - return element0_interact_path1(threadindex, rng_states, neutron) - if element_index == 1: - return element1_interact_path1(threadindex, rng_states, neutron) - if element_index == 2: - return element2_interact_path1(threadindex, rng_states, neutron) - return none - - @cuda.jit(device=True) - def element_calculate_attenuation(neutron, end, element_index): - if element_index == 0: - return element0_calculate_attenuation(neutron, end) - if element_index == 1: - return element1_calculate_attenuation(neutron, end) - if element_index == 2: - return element2_calculate_attenuation(neutron, end) - return 1.0 - - if test.USE_CUDASIM: - @cuda.jit(device=True, inline=True) - def interact_path1(threadindex, rng_states, neutron): - tmp_neutron = np.zeros(10, dtype=float) - return _interact_path1(threadindex, rng_states, neutron, tmp_neutron) - @cuda.jit(device=True, inline=True) - def calculate_attenuation(neutron, end): - tmp_neutron = np.zeros(10, dtype=float) - return _calculate_attenuation(neutron, end, tmp_neutron) - else: - @cuda.jit(device=True, inline=True) - def interact_path1(threadindex, rng_states, neutron): - tmp_neutron = cuda.local.array(10, dtype=numba.float64) - return _interact_path1(threadindex, rng_states, neutron, tmp_neutron) - @cuda.jit(device=True, inline=True) - def calculate_attenuation(neutron, end): - tmp_neutron = cuda.local.array(10, dtype=numba.float64) - return _calculate_attenuation(neutron, end, tmp_neutron) - return dict( - interact_path1 = interact_path1, - calculate_attenuation = calculate_attenuation, - ) diff --git a/tests/scatterers/test_composite_scatterer.py b/tests/scatterers/test_composite_scatterer.py index 2f502f99..44b1e2bf 100644 --- a/tests/scatterers/test_composite_scatterer.py +++ b/tests/scatterers/test_composite_scatterer.py @@ -10,7 +10,7 @@ def composite_example_interact_path1(): path = os.path.join(thisdir, "sampleassemblies", 'sample+2cylinders', 'sampleassembly.xml') from mcvine.acc.components.samples import loadScattererComposite composite = loadScattererComposite(path) - from mcvine.acc.scatterers.composite_scatterer import factory_3 + from mcvine.acc.scatterers.composite_3 import factory_3 methods = factory_3(composite) return methods['interact_path1'] From e6dff8b7e7e14c0aacf0f26b0f99b4513da52ff4 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Wed, 10 May 2023 22:15:58 -0700 Subject: [PATCH 37/53] trivial: acc/geometry/composite_3.py --- acc/geometry/composite_3.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/acc/geometry/composite_3.py b/acc/geometry/composite_3.py index 06f8fe23..8424faf1 100644 --- a/acc/geometry/composite_3.py +++ b/acc/geometry/composite_3.py @@ -1,4 +1,5 @@ -# Example code +# Example code that would be generated automatically + import os, numpy as np, numba from numba import cuda from mcvine.acc._numba import xoroshiro128p_uniform_float32 From 946a5513ad792e9ab4a9b897462f001e0d0b3b54 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Sat, 13 May 2023 08:23:35 -0700 Subject: [PATCH 38/53] make most methods protected in geometry.composite, and added two public methods get_find_1st_hit and get_union_locate --- acc/geometry/composite.py | 31 ++++++++++++++------- acc/geometry/composite_3.py | 2 +- tests/geometry/test_composite_3_elements.py | 15 ++++++---- tests/geometry/test_composite_shape.py | 14 ++++++---- 4 files changed, 39 insertions(+), 23 deletions(-) diff --git a/acc/geometry/composite.py b/acc/geometry/composite.py index 39442419..a5128a8d 100644 --- a/acc/geometry/composite.py +++ b/acc/geometry/composite.py @@ -4,7 +4,18 @@ from mcvine.acc import test from .._numba import coder -def make_find_1st_hit(forward_intersect_all, is_onborder, find_shape_containing_point, **kwds): +def get_find_1st_hit(shapes): + N = len(shapes) + mod = _importModule(N) + methods = mod.createRayTracingMethods_NonOverlappingShapes(shapes) + return _make_find_1st_hit(**methods) + +def get_union_locate(shapes): + N = len(shapes) + mod = _importModule(N) + return mod.createUnionLocateMethod(shapes) + +def _make_find_1st_hit(forward_intersect_all, is_onborder, find_shape_containing_point, **kwds): @cuda.jit(device=True) def _find_1st_hit(x,y,z, vx,vy,vz, ts): nIntersections = forward_intersect_all(x,y,z, vx,vy,vz, ts) @@ -71,12 +82,12 @@ def find_1st_hit(x,y,z, vx,vy,vz): return find_1st_hit -def importModule(N): - mod = makeModule(N) +def _importModule(N): + mod = _makeModule(N) import imp return imp.load_source('mod', mod) -def makeModule(N, overwrite=False): +def _makeModule(N, overwrite=False): "make cuda device methods for composite with N elements" from .._numba import coder modulepath = coder.getModule("composite_shape", N) @@ -91,14 +102,14 @@ def makeModule(N, overwrite=False): from mcvine.acc.geometry.location import inside, onborder, outside from mcvine.acc.geometry.arrow_intersect import max_intersections, insert_into_sorted_list """.splitlines() - createMethods = Coder_createMethods(N, indent)() - createUnionLocateMethod = Coder_createUnionLocateMethod(N, indent)() - lines = imports + [''] + createMethods + createUnionLocateMethod + createRTMethods = _Coder_createRTMethods(N, indent)() + createUnionLocateMethod = _Coder_createUnionLocateMethod(N, indent)() + lines = imports + [''] + createRTMethods + createUnionLocateMethod with open(modulepath, 'wt') as ostream: ostream.write("\n".join(lines)) return modulepath -class Coder_createUnionLocateMethod: +class _Coder_createUnionLocateMethod: def __init__(self, N, indent=4*' '): self.N = N @@ -151,7 +162,7 @@ def locate(self): return header + assign_loop + inside_loop + onborder_loop -class Coder_createMethods: +class _Coder_createRTMethods: def __init__(self, N, indent=4*' '): self.N = N @@ -209,7 +220,7 @@ def forward_intersect_all(x,y,z, vx,vy,vz, ts): + end ) add_indent = lambda lines, n: [indent*n+l for l in lines] - return ["def createMethods(shapes):"] + add_indent(body, 1) + return ["def createRayTracingMethods_NonOverlappingShapes(shapes):"] + add_indent(body, 1) def _intersect_all(self): diff --git a/acc/geometry/composite_3.py b/acc/geometry/composite_3.py index 8424faf1..1585a478 100644 --- a/acc/geometry/composite_3.py +++ b/acc/geometry/composite_3.py @@ -8,7 +8,7 @@ from mcvine.acc.geometry.location import inside, onborder, outside from mcvine.acc.geometry.arrow_intersect import max_intersections, insert_into_sorted_list -def createMethods_3(shapes): +def createRayTracingMethods_3_NonOverlappingShapes(shapes): "methods for a list of non-overlapping shapes" assert len(shapes)==3 funcs_list = [ diff --git a/tests/geometry/test_composite_3_elements.py b/tests/geometry/test_composite_3_elements.py index 5aac1fa3..9e743bdf 100644 --- a/tests/geometry/test_composite_3_elements.py +++ b/tests/geometry/test_composite_3_elements.py @@ -7,8 +7,11 @@ from instrument.nixml import parse_file from mcvine.acc import test from mcvine.acc.geometry import locate, location, arrow_intersect -from mcvine.acc.geometry.composite_3 import createMethods_3, createUnionLocateMethod_3 -from mcvine.acc.geometry.composite import make_find_1st_hit +from mcvine.acc.geometry.composite_3 import ( + createRayTracingMethods_3_NonOverlappingShapes as createRTMethods_3, + createUnionLocateMethod_3 +) +from mcvine.acc.geometry.composite import _make_find_1st_hit thisdir = os.path.dirname(__file__) @@ -31,7 +34,7 @@ def test_locate(): def test_intersect_all(): union = parse_file(os.path.join(thisdir, 'union_three_elements.xml'))[0] shapes = union.shapes - methods = createMethods_3(shapes) + methods = createRTMethods_3(shapes) intersect_all = methods['intersect_all'] ts = np.zeros(11) assert intersect_all(0.,0.,0., 1.,0.,0., ts) == 10 @@ -43,7 +46,7 @@ def test_intersect_all(): def test_forward_intersect_all(): union = parse_file(os.path.join(thisdir, 'union_three_elements.xml'))[0] shapes = union.shapes - methods = createMethods_3(shapes) + methods = createRTMethods_3(shapes) forward_intersect_all = methods['forward_intersect_all'] ts = np.zeros(6) assert forward_intersect_all(0.,0.,0., 1.,0.,0., ts) == 5 @@ -55,8 +58,8 @@ def test_forward_intersect_all(): def test_find_1st_hit(): union = parse_file(os.path.join(thisdir, 'union_three_elements.xml'))[0] shapes = union.shapes - methods = createMethods_3(shapes) - find_1st_hit = make_find_1st_hit(**methods) + methods = createRTMethods_3(shapes) + find_1st_hit = _make_find_1st_hit(**methods) assert find_1st_hit(-0.25,0.,0., 1.,0.,0.) == 0 assert find_1st_hit(-0.15,0.,0., 1.,0.,0.) == 1 assert find_1st_hit(-0.05,0.,0., 1.,0.,0.) == 2 diff --git a/tests/geometry/test_composite_shape.py b/tests/geometry/test_composite_shape.py index 729d7795..fd1bc893 100644 --- a/tests/geometry/test_composite_shape.py +++ b/tests/geometry/test_composite_shape.py @@ -7,20 +7,22 @@ from instrument.nixml import parse_file from mcvine.acc import test from mcvine.acc.geometry import locate -from mcvine.acc.geometry.composite import makeModule, importModule, make_find_1st_hit +from mcvine.acc.geometry.composite import _make_find_1st_hit thisdir = os.path.dirname(__file__) def test_makeModule(): - makeModule(4, overwrite=True) + from mcvine.acc.geometry.composite import _makeModule + _makeModule(4, overwrite=True) return @pytest.fixture def ray_tracing_methods(): + from mcvine.acc.geometry.composite import _importModule union = parse_file(os.path.join(thisdir, 'union_three_elements.xml'))[0] shapes = union.shapes - mod = importModule(len(shapes)) - methods = mod.createMethods(shapes) + mod = _importModule(len(shapes)) + methods = mod.createRayTracingMethods_NonOverlappingShapes(shapes) methods['union_locate'] = mod.createUnionLocateMethod(shapes) return methods @@ -57,7 +59,7 @@ def test_forward_intersect_all(ray_tracing_methods): @pytest.mark.skipif(not test.USE_CUDASIM, reason='no CUDASIM') def test_find_1st_hit(ray_tracing_methods): - find_1st_hit = make_find_1st_hit(**ray_tracing_methods) + find_1st_hit = _make_find_1st_hit(**ray_tracing_methods) assert find_1st_hit(-0.25,0.,0., 1.,0.,0.) == 0 assert find_1st_hit(-0.15,0.,0., 1.,0.,0.) == 1 assert find_1st_hit(-0.05,0.,0., 1.,0.,0.) == 2 @@ -67,7 +69,7 @@ def test_find_1st_hit(ray_tracing_methods): @pytest.mark.skipif(not test.USE_CUDA, reason='no CUDA') def test_cuda(ray_tracing_methods): - find_1st_hit = make_find_1st_hit(**ray_tracing_methods) + find_1st_hit = _make_find_1st_hit(**ray_tracing_methods) @cuda.jit def find_1st_hit_kernel(neutrons, out, n_neutrons_per_thread): N = len(neutrons) From 35f2e536937ec4f22c83fddfc519c67514b5266e Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Sat, 13 May 2023 08:43:23 -0700 Subject: [PATCH 39/53] added a sanity check of public methods in geometry.composite: tests/geometry/test_composite_shape.py --- tests/geometry/test_composite_shape.py | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/tests/geometry/test_composite_shape.py b/tests/geometry/test_composite_shape.py index fd1bc893..558aa396 100644 --- a/tests/geometry/test_composite_shape.py +++ b/tests/geometry/test_composite_shape.py @@ -17,10 +17,13 @@ def test_makeModule(): return @pytest.fixture -def ray_tracing_methods(): - from mcvine.acc.geometry.composite import _importModule +def shapes(): union = parse_file(os.path.join(thisdir, 'union_three_elements.xml'))[0] - shapes = union.shapes + return union.shapes + +@pytest.fixture +def ray_tracing_methods(shapes): + from mcvine.acc.geometry.composite import _importModule mod = _importModule(len(shapes)) methods = mod.createRayTracingMethods_NonOverlappingShapes(shapes) methods['union_locate'] = mod.createUnionLocateMethod(shapes) @@ -96,3 +99,11 @@ def run(neutrons, first_hit): run(neutrons, first_hit) print(first_hit) return + +@pytest.mark.skipif(not test.USE_CUDASIM, reason='no CUDASIM') +def test_public_methods(shapes): + from mcvine.acc.geometry.composite import get_find_1st_hit, get_union_locate + get_find_1st_hit(shapes) + get_union_locate(shapes) + return + From 67b87f37bf1baed01bda1fc93fd1362fde3e9049 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Sat, 13 May 2023 09:09:13 -0700 Subject: [PATCH 40/53] minor: made method names in geometry.composite_3 consistent with those generated in geometry.composite --- acc/geometry/composite_3.py | 4 ++-- tests/geometry/test_composite_3_elements.py | 12 ++++++------ 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/acc/geometry/composite_3.py b/acc/geometry/composite_3.py index 1585a478..02a80c63 100644 --- a/acc/geometry/composite_3.py +++ b/acc/geometry/composite_3.py @@ -8,7 +8,7 @@ from mcvine.acc.geometry.location import inside, onborder, outside from mcvine.acc.geometry.arrow_intersect import max_intersections, insert_into_sorted_list -def createRayTracingMethods_3_NonOverlappingShapes(shapes): +def createRayTracingMethods_NonOverlappingShapes(shapes): "methods for a list of non-overlapping shapes" assert len(shapes)==3 funcs_list = [ @@ -103,7 +103,7 @@ def forward_intersect_all(x,y,z, vx,vy,vz, ts): is_onborder = is_onborder, ) -def createUnionLocateMethod_3(shapes): +def createUnionLocateMethod(shapes): assert len(shapes)==3 locates = [ arrow_intersect.locate_func_factory.render(shape) diff --git a/tests/geometry/test_composite_3_elements.py b/tests/geometry/test_composite_3_elements.py index 9e743bdf..f379aca3 100644 --- a/tests/geometry/test_composite_3_elements.py +++ b/tests/geometry/test_composite_3_elements.py @@ -8,8 +8,8 @@ from mcvine.acc import test from mcvine.acc.geometry import locate, location, arrow_intersect from mcvine.acc.geometry.composite_3 import ( - createRayTracingMethods_3_NonOverlappingShapes as createRTMethods_3, - createUnionLocateMethod_3 + createRayTracingMethods_NonOverlappingShapes as createRTMethods, + createUnionLocateMethod ) from mcvine.acc.geometry.composite import _make_find_1st_hit @@ -20,7 +20,7 @@ def test_locate(): union = parse_file(os.path.join(thisdir, 'union_three_elements.xml'))[0] shapes = union.shapes - locate_u3 = createUnionLocateMethod_3(shapes) + locate_u3 = createUnionLocateMethod(shapes) assert locate_u3(0, 0, 0) == locate.inside assert locate_u3(0.025, 0, 0) == locate.onborder assert locate_u3(0.099, 0, 0) == locate.onborder @@ -34,7 +34,7 @@ def test_locate(): def test_intersect_all(): union = parse_file(os.path.join(thisdir, 'union_three_elements.xml'))[0] shapes = union.shapes - methods = createRTMethods_3(shapes) + methods = createRTMethods(shapes) intersect_all = methods['intersect_all'] ts = np.zeros(11) assert intersect_all(0.,0.,0., 1.,0.,0., ts) == 10 @@ -46,7 +46,7 @@ def test_intersect_all(): def test_forward_intersect_all(): union = parse_file(os.path.join(thisdir, 'union_three_elements.xml'))[0] shapes = union.shapes - methods = createRTMethods_3(shapes) + methods = createRTMethods(shapes) forward_intersect_all = methods['forward_intersect_all'] ts = np.zeros(6) assert forward_intersect_all(0.,0.,0., 1.,0.,0., ts) == 5 @@ -58,7 +58,7 @@ def test_forward_intersect_all(): def test_find_1st_hit(): union = parse_file(os.path.join(thisdir, 'union_three_elements.xml'))[0] shapes = union.shapes - methods = createRTMethods_3(shapes) + methods = createRTMethods(shapes) find_1st_hit = _make_find_1st_hit(**methods) assert find_1st_hit(-0.25,0.,0., 1.,0.,0.) == 0 assert find_1st_hit(-0.15,0.,0., 1.,0.,0.) == 1 From 52b93ab9e21e953351566c5b019cebd5373d016e Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Sat, 13 May 2023 09:55:24 -0700 Subject: [PATCH 41/53] fixed bugs introduced in the last couple commits (method renames): acc/geometry/{arrow_intersect|locate}.py --- acc/geometry/arrow_intersect.py | 4 ++-- acc/geometry/locate.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/acc/geometry/arrow_intersect.py b/acc/geometry/arrow_intersect.py index ad5bbfb7..51d9e68d 100644 --- a/acc/geometry/arrow_intersect.py +++ b/acc/geometry/arrow_intersect.py @@ -29,8 +29,8 @@ def onUnion(self, u): nelements = len(u.shapes) if nelements == 2: return self.onUnion2(u) if nelements == 3: - from .composite_3 import createMethods_3 - return createMethods_3(u.shapes)['intersect_all'] + from .composite_3 import createRayTracingMethods_NonOverlappingShapes as createMethods + return createMethods(u.shapes)['intersect_all'] raise NotImplementedError(f"locate for union of {nelements} elements") def onUnion2(self, u): diff --git a/acc/geometry/locate.py b/acc/geometry/locate.py index 979c1779..5d249f0f 100644 --- a/acc/geometry/locate.py +++ b/acc/geometry/locate.py @@ -15,8 +15,8 @@ def onUnion(self, u): nelements = len(u.shapes) if nelements == 2: return self.onUnion2(u) if nelements == 3: - from .composite_3 import createUnionLocateMethod_3 - return createUnionLocateMethod_3(u.shapes) + from .composite_3 import createUnionLocateMethod + return createUnionLocateMethod(u.shapes) raise NotImplementedError(f"locate for union of {nelements} elements") def onUnion2(self, u): From 994272de6be65193c4d5ff7f4afd5346bea8662f Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Sat, 13 May 2023 09:56:49 -0700 Subject: [PATCH 42/53] refactored codes around composite_scatterer so that we can add "coder" for composites of arbitrary number of elements --- acc/scatterers/__init__.py | 4 +- acc/scatterers/composite_3.py | 116 ++---------------- acc/scatterers/composite_scatterer.py | 118 +++++++++++++++++-- tests/scatterers/test_composite_scatterer.py | 4 +- 4 files changed, 122 insertions(+), 120 deletions(-) diff --git a/acc/scatterers/__init__.py b/acc/scatterers/__init__.py index 208a41af..68f3f09d 100644 --- a/acc/scatterers/__init__.py +++ b/acc/scatterers/__init__.py @@ -13,8 +13,8 @@ def onCompositeScatterer(self, composite): elements = composite.elements() if len(elements) != 3: raise NotImplementedError - from .composite_3 import factory_3 - return factory_3(composite) + from .composite_scatterer import factory + return factory(composite) def onHomogeneousScatterer(self, hs): from .homogeneous_scatterer import factory diff --git a/acc/scatterers/composite_3.py b/acc/scatterers/composite_3.py index f02c0a0e..c4ee27d5 100644 --- a/acc/scatterers/composite_3.py +++ b/acc/scatterers/composite_3.py @@ -4,56 +4,23 @@ # import numpy as np, numba -from numba import cuda, void, int64 -from mcvine.acc._numba import xoroshiro128p_uniform_float32 -from math import sqrt, exp - -from .interaction_types import absorption, scattering, none -from .. import test -from ..neutron import absorb, prop_dt_inplace, clone -from ..geometry.arrow_intersect import max_intersections -from ..geometry.locate import inside, outside, onborder -from ..geometry.propagation import makePropagateMethods - - +from numba import cuda from numba.core import config -if not config.ENABLE_CUDASIM: - from numba.cuda.compiler import Dispatcher, DeviceFunction - -from ..config import get_numba_floattype -NB_FLOAT = get_numba_floattype() +from .interaction_types import none +from .. import test -def factory_3(composite): - # methods to deal with the whole shape of the composite - shape = composite.shape() - from ..geometry import arrow_intersect - intersect = arrow_intersect.arrow_intersect_func_factory.render(shape) - locate = arrow_intersect.locate_func_factory.render(shape) - del shape - propagate_methods = makePropagateMethods(intersect, locate) - del intersect - propagate_to_next_exiting_surface = propagate_methods['propagate_to_next_exiting_surface'] - del propagate_methods +def createHelperMethodsForScatter(composite): # elements elements = composite.elements() nelements = len(elements) assert nelements == 3 - # find_1st_hit - shapes = [e.shape() for e in elements] - from mcvine.acc.geometry.composite import make_find_1st_hit - from mcvine.acc.geometry.composite_3 import createMethods_3 - methods = createMethods_3(shapes) - del shapes - find_1st_hit = make_find_1st_hit(**methods) # methods for element scatterers from . import scatter_func_factory - #element_scatter_methods = [scatter_func_factory.render(e) for e in elements] - element_scatter_methods = [] - for e in elements: - m = scatter_func_factory.render(e) - element_scatter_methods.append(m) - del scatter_func_factory + element_scatter_methods = [ + scatter_func_factory.render(e) + for e in elements + ] element0_interact_path1 = element_scatter_methods[0]['interact_path1'] element1_interact_path1 = element_scatter_methods[1]['interact_path1'] element2_interact_path1 = element_scatter_methods[2]['interact_path1'] @@ -62,51 +29,6 @@ def factory_3(composite): element2_calculate_attenuation = element_scatter_methods[2]['calculate_attenuation'] del elements, element_scatter_methods - @cuda.jit(device=True) - def _interact_path1(threadindex, rng_states, neutron, tmp_neutron): - x,y,z = neutron[:3] - vx,vy,vz = neutron[3:6] - scatterer_index = find_1st_hit(x,y,z, vx,vy,vz) - if scatterer_index < 0 or scatterer_index > nelements: - propagate_to_next_exiting_surface(neutron) - return none - clone(neutron, tmp_neutron) - # global2local(tmp_neutron, scatterer_index) - interaction = element_interact_path1(threadindex, rng_states, tmp_neutron, scatterer_index) - # local2global(tmp_neutron, scatterer_index) - clone(tmp_neutron, neutron) - if interaction == absorption: - absorb(neutron) - return interaction - # nothing happened - if interaction == none: - # if neutron is not inside, we are done - x,y,z = neutron[:3] - if locate(x,y,z)!=inside: - return interaction - # otherwise, interact again - # return _interact_path1(threadindex, rng_states, neutron, tmp_neutron) - return interaction - # interaction must be scatter - clone(neutron, tmp_neutron) - # propagate to surface if necessary - x,y,z = neutron[:3] - if locate(x,y,z) == inside: - propagate_to_next_exiting_surface(neutron) - # apply attenuation - att = calculate_attenuation(tmp_neutron, neutron[:3]) - neutron[-1] *= att - return scattering - - @cuda.jit(device=True) - def _calculate_attenuation(neutron, end, tmp_neutron): - ret = 1.0 - for i in range(nelements): - clone(neutron, tmp_neutron) - # global2local(tmp_neutron, i) - ret *= element_calculate_attenuation(neutron, end, i) - return ret - @cuda.jit(device=True) def element_interact_path1(threadindex, rng_states, neutron, element_index): if element_index == 0: @@ -127,25 +49,7 @@ def element_calculate_attenuation(neutron, end, element_index): return element2_calculate_attenuation(neutron, end) return 1.0 - if test.USE_CUDASIM: - @cuda.jit(device=True, inline=True) - def interact_path1(threadindex, rng_states, neutron): - tmp_neutron = np.zeros(10, dtype=float) - return _interact_path1(threadindex, rng_states, neutron, tmp_neutron) - @cuda.jit(device=True, inline=True) - def calculate_attenuation(neutron, end): - tmp_neutron = np.zeros(10, dtype=float) - return _calculate_attenuation(neutron, end, tmp_neutron) - else: - @cuda.jit(device=True, inline=True) - def interact_path1(threadindex, rng_states, neutron): - tmp_neutron = cuda.local.array(10, dtype=numba.float64) - return _interact_path1(threadindex, rng_states, neutron, tmp_neutron) - @cuda.jit(device=True, inline=True) - def calculate_attenuation(neutron, end): - tmp_neutron = cuda.local.array(10, dtype=numba.float64) - return _calculate_attenuation(neutron, end, tmp_neutron) return dict( - interact_path1 = interact_path1, - calculate_attenuation = calculate_attenuation, + element_interact_path1 = element_interact_path1, + element_calculate_attenuation = element_calculate_attenuation, ) diff --git a/acc/scatterers/composite_scatterer.py b/acc/scatterers/composite_scatterer.py index 3c7244cc..66842330 100644 --- a/acc/scatterers/composite_scatterer.py +++ b/acc/scatterers/composite_scatterer.py @@ -3,23 +3,121 @@ # Jiao Lin # -import numpy as np, numba -from numba import cuda, void, int64 -from mcvine.acc._numba import xoroshiro128p_uniform_float32 from math import sqrt, exp +import numpy as np, numba +from numba import cuda +from numba.core import config -from .interaction_types import absorption, scattering, none from .. import test +from .._numba import xoroshiro128p_uniform_float32 +from .interaction_types import absorption, scattering, none from ..neutron import absorb, prop_dt_inplace, clone -from ..geometry.arrow_intersect import max_intersections from ..geometry.locate import inside, outside, onborder -from ..geometry.propagation import makePropagateMethods -from numba.core import config -if not config.ENABLE_CUDASIM: - from numba.cuda.compiler import Dispatcher, DeviceFunction - from ..config import get_numba_floattype NB_FLOAT = get_numba_floattype() +def _createRayTracingMethods(composite): + """create ray-tracing methods to deal with the overall shape of the composite + """ + shape = composite.shape() + from ..geometry import arrow_intersect + intersect = arrow_intersect.arrow_intersect_func_factory.render(shape) + locate = arrow_intersect.locate_func_factory.render(shape) + del shape + from ..geometry.propagation import makePropagateMethods + propagate_methods = makePropagateMethods(intersect, locate) + propagate_to_next_exiting_surface = propagate_methods['propagate_to_next_exiting_surface'] + del propagate_methods + # find_1st_hit + elements = composite.elements() + shapes = [e.shape() for e in elements] + from mcvine.acc.geometry.composite import get_find_1st_hit + find_1st_hit = get_find_1st_hit(shapes) + return dict( + locate = locate, + intersect = intersect, + find_1st_hit = find_1st_hit, + propagate_to_next_exiting_surface = propagate_to_next_exiting_surface, + ) + +def factory(composite): + elements = composite.elements() + N = len(elements) + assert N==3 + from .composite_3 import createHelperMethodsForScatter + element_methods = createHelperMethodsForScatter(composite) + element_interact_path1 = element_methods["element_interact_path1"] + element_calculate_attenuation = element_methods["element_calculate_attenuation"] + rt_methods = _createRayTracingMethods(composite) + find_1st_hit = rt_methods['find_1st_hit'] + locate = rt_methods['locate'] + propagate_to_next_exiting_surface = rt_methods['propagate_to_next_exiting_surface'] + @cuda.jit(device=True) + def _interact_path1(threadindex, rng_states, neutron, tmp_neutron): + x,y,z = neutron[:3] + vx,vy,vz = neutron[3:6] + scatterer_index = find_1st_hit(x,y,z, vx,vy,vz) + if scatterer_index < 0 or scatterer_index > N: + propagate_to_next_exiting_surface(neutron) + return none + clone(neutron, tmp_neutron) + # global2local(tmp_neutron, scatterer_index) + interaction = element_interact_path1(threadindex, rng_states, tmp_neutron, scatterer_index) + # local2global(tmp_neutron, scatterer_index) + clone(tmp_neutron, neutron) + if interaction == absorption: + absorb(neutron) + return interaction + # nothing happened + if interaction == none: + # if neutron is not inside, we are done + x,y,z = neutron[:3] + if locate(x,y,z)!=inside: + return interaction + # otherwise, interact again + # return _interact_path1(threadindex, rng_states, neutron, tmp_neutron) + return interaction + # interaction must be scatter + clone(neutron, tmp_neutron) + # propagate to surface if necessary + x,y,z = neutron[:3] + if locate(x,y,z) == inside: + propagate_to_next_exiting_surface(neutron) + # apply attenuation + att = calculate_attenuation(tmp_neutron, neutron[:3]) + neutron[-1] *= att + return scattering + + @cuda.jit(device=True) + def _calculate_attenuation(neutron, end, tmp_neutron): + ret = 1.0 + for i in range(N): + clone(neutron, tmp_neutron) + # global2local(tmp_neutron, i) + ret *= element_calculate_attenuation(neutron, end, i) + return ret + + if test.USE_CUDASIM: + @cuda.jit(device=True, inline=True) + def interact_path1(threadindex, rng_states, neutron): + tmp_neutron = np.zeros(10, dtype=float) + return _interact_path1(threadindex, rng_states, neutron, tmp_neutron) + @cuda.jit(device=True, inline=True) + def calculate_attenuation(neutron, end): + tmp_neutron = np.zeros(10, dtype=float) + return _calculate_attenuation(neutron, end, tmp_neutron) + else: + @cuda.jit(device=True, inline=True) + def interact_path1(threadindex, rng_states, neutron): + tmp_neutron = cuda.local.array(10, dtype=numba.float64) + return _interact_path1(threadindex, rng_states, neutron, tmp_neutron) + @cuda.jit(device=True, inline=True) + def calculate_attenuation(neutron, end): + tmp_neutron = cuda.local.array(10, dtype=numba.float64) + return _calculate_attenuation(neutron, end, tmp_neutron) + return dict( + interact_path1 = interact_path1, + calculate_attenuation = calculate_attenuation, + ) diff --git a/tests/scatterers/test_composite_scatterer.py b/tests/scatterers/test_composite_scatterer.py index 44b1e2bf..56bed2f5 100644 --- a/tests/scatterers/test_composite_scatterer.py +++ b/tests/scatterers/test_composite_scatterer.py @@ -10,8 +10,8 @@ def composite_example_interact_path1(): path = os.path.join(thisdir, "sampleassemblies", 'sample+2cylinders', 'sampleassembly.xml') from mcvine.acc.components.samples import loadScattererComposite composite = loadScattererComposite(path) - from mcvine.acc.scatterers.composite_3 import factory_3 - methods = factory_3(composite) + from mcvine.acc.scatterers.composite_scatterer import factory + methods = factory(composite) return methods['interact_path1'] @pytest.mark.skipif(not test.USE_CUDASIM, reason='no CUDASIM') From 2ed20ce8df81d6011d8f41dc8dd8b1bf8d8917ea Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Sat, 13 May 2023 10:26:38 -0700 Subject: [PATCH 43/53] added scatterers.composite_scatterer._makeModule --- acc/geometry/composite.py | 2 +- acc/scatterers/composite_scatterer.py | 111 +++++++++++++++++++ tests/scatterers/test_composite_scatterer.py | 4 + 3 files changed, 116 insertions(+), 1 deletion(-) diff --git a/acc/geometry/composite.py b/acc/geometry/composite.py index a5128a8d..0880f2fc 100644 --- a/acc/geometry/composite.py +++ b/acc/geometry/composite.py @@ -85,7 +85,7 @@ def find_1st_hit(x,y,z, vx,vy,vz): def _importModule(N): mod = _makeModule(N) import imp - return imp.load_source('mod', mod) + return imp.load_source(f'geometry_composite_{N}', mod) def _makeModule(N, overwrite=False): "make cuda device methods for composite with N elements" diff --git a/acc/scatterers/composite_scatterer.py b/acc/scatterers/composite_scatterer.py index 66842330..82ee81c9 100644 --- a/acc/scatterers/composite_scatterer.py +++ b/acc/scatterers/composite_scatterer.py @@ -3,6 +3,7 @@ # Jiao Lin # +import os from math import sqrt, exp import numpy as np, numba from numba import cuda @@ -121,3 +122,113 @@ def calculate_attenuation(neutron, end): interact_path1 = interact_path1, calculate_attenuation = calculate_attenuation, ) + + +def _importModule(N): + mod = _makeModule(N) + import imp + return imp.load_source(f"scatterer_composite_{N}", mod) + +def _makeModule(N, overwrite=False): + "make a python module for cuda device methods for composite scatterer with N elements" + modulepath = coder.getModule("composite_scatterer", N) + if os.path.exists(modulepath) and not overwrite: + return modulepath + indent = 4*' ' + imports = """import os, numpy as np, numba +from numba import cuda +from numba.core import config + +from mcvine.acc.scatterers.interaction_types import none +from mcvine.acc import test +""".splitlines() + createHelperMethodsForScatter = _Coder_createHelperMethodsForScatter(N, indent)() + lines = imports + [''] + createHelperMethodsForScatter + with open(modulepath, 'wt') as ostream: + ostream.write("\n".join(lines)) + return modulepath + +from .._numba import coder +class _Coder_createHelperMethodsForScatter: + + def __init__(self, N, indent=4*' '): + self.N = N + self.indent = indent + return + + def __call__(self): + N, indent = self.N, self.indent + header = f"""# elements +elements = composite.elements() +nelements = len(elements) +# methods for element scatterers +from . import scatter_func_factory +element_scatter_methods = [ + scatter_func_factory.render(e) + for e in elements +] +""".splitlines() + interact_path1_loop = coder.unrollLoop( + N = N, + indent = '', + in_loop = ["element{i}_interact_path1 = element_scatter_methods[{i}]['interact_path1']"], + ) + calculate_attenuation_loop = coder.unrollLoop( + N = N, + indent = '', + in_loop = ["element{i}_calcualte_attenuation = element_scatter_methods[{i}]['calculate_attenuation']"], + ) + footer = """ +return dict( + element_interact_path1 = element_interact_path1, + element_calculate_attenuation = element_calculate_attenuation, +) +""".splitlines() + body = ( + header + + interact_path1_loop + + calculate_attenuation_loop + + [''] + + self.element_interact_path1() + + self.element_calculate_attenuation() + + footer + ) + add_indent = lambda lines, n: [indent*n+l for l in lines] + return ( + ["def createHelperMethodsForScatter(composite):"] + + add_indent(body, 1) + ) + + def element_interact_path1(self): + N, indent = self.N, self.indent + header = [ + "@cuda.jit(device=True)", + "def element_interact_path1(threadindex, rng_states, neutron, element_index):" + ] + loop = coder.unrollLoop( + N = N, + indent = indent, + in_loop = [ + "if element_index == {i}:", + indent + "return element{i}_interact_path1(threadindex, rng_states, neutron)", + ] + ) + return header + loop + [indent + 'return none'] + + def element_calculate_attenuation(self): + N, indent = self.N, self.indent + header = [ + "@cuda.jit(device=True)", + "def element_calculate_attenuation(neutron, end, element_index):" + ] + loop = coder.unrollLoop( + N = N, + indent = indent, + in_loop = [ + "if element_index == {i}:", + indent + "return element{i}_calculate_attenuation(neutron, end)", + ] + ) + return header + loop + [indent + 'return 1.0'] + + diff --git a/tests/scatterers/test_composite_scatterer.py b/tests/scatterers/test_composite_scatterer.py index 56bed2f5..d1dc2ffc 100644 --- a/tests/scatterers/test_composite_scatterer.py +++ b/tests/scatterers/test_composite_scatterer.py @@ -66,3 +66,7 @@ def run(neutrons): run(neutrons) print(neutrons) return + +def test_makeModule(): + from mcvine.acc.scatterers.composite_scatterer import _makeModule + _makeModule(4, overwrite=True) From ed64b716d82563bb4f47288e0cbda089902afbde Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Sat, 13 May 2023 10:34:49 -0700 Subject: [PATCH 44/53] use module created by "coder": acc/scatterers/composite_scatterer.py --- acc/scatterers/composite_scatterer.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/acc/scatterers/composite_scatterer.py b/acc/scatterers/composite_scatterer.py index 82ee81c9..3b3125c4 100644 --- a/acc/scatterers/composite_scatterer.py +++ b/acc/scatterers/composite_scatterer.py @@ -46,9 +46,8 @@ def _createRayTracingMethods(composite): def factory(composite): elements = composite.elements() N = len(elements) - assert N==3 - from .composite_3 import createHelperMethodsForScatter - element_methods = createHelperMethodsForScatter(composite) + mod = _importModule(N) + element_methods = mod.createHelperMethodsForScatter(composite) element_interact_path1 = element_methods["element_interact_path1"] element_calculate_attenuation = element_methods["element_calculate_attenuation"] rt_methods = _createRayTracingMethods(composite) @@ -162,7 +161,7 @@ def __call__(self): elements = composite.elements() nelements = len(elements) # methods for element scatterers -from . import scatter_func_factory +from mcvine.acc.scatterers import scatter_func_factory element_scatter_methods = [ scatter_func_factory.render(e) for e in elements @@ -176,7 +175,7 @@ def __call__(self): calculate_attenuation_loop = coder.unrollLoop( N = N, indent = '', - in_loop = ["element{i}_calcualte_attenuation = element_scatter_methods[{i}]['calculate_attenuation']"], + in_loop = ["element{i}_calculate_attenuation = element_scatter_methods[{i}]['calculate_attenuation']"], ) footer = """ return dict( From 5075426f3bde52c1672dc63ec72331da8b386723 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Sun, 14 May 2023 11:14:53 -0700 Subject: [PATCH 45/53] composite: implement and test single-scatter methods --- acc/scatterers/composite_scatterer.py | 111 ++++++++++++++----- tests/scatterers/test_composite_scatterer.py | 70 +++++++++++- 2 files changed, 151 insertions(+), 30 deletions(-) diff --git a/acc/scatterers/composite_scatterer.py b/acc/scatterers/composite_scatterer.py index 3b3125c4..93db9194 100644 --- a/acc/scatterers/composite_scatterer.py +++ b/acc/scatterers/composite_scatterer.py @@ -12,37 +12,13 @@ from .. import test from .._numba import xoroshiro128p_uniform_float32 from .interaction_types import absorption, scattering, none -from ..neutron import absorb, prop_dt_inplace, clone +from ..neutron import absorb, prop_dt_inplace, clone, is_moving from ..geometry.locate import inside, outside, onborder from ..config import get_numba_floattype NB_FLOAT = get_numba_floattype() -def _createRayTracingMethods(composite): - """create ray-tracing methods to deal with the overall shape of the composite - """ - shape = composite.shape() - from ..geometry import arrow_intersect - intersect = arrow_intersect.arrow_intersect_func_factory.render(shape) - locate = arrow_intersect.locate_func_factory.render(shape) - del shape - from ..geometry.propagation import makePropagateMethods - propagate_methods = makePropagateMethods(intersect, locate) - propagate_to_next_exiting_surface = propagate_methods['propagate_to_next_exiting_surface'] - del propagate_methods - # find_1st_hit - elements = composite.elements() - shapes = [e.shape() for e in elements] - from mcvine.acc.geometry.composite import get_find_1st_hit - find_1st_hit = get_find_1st_hit(shapes) - return dict( - locate = locate, - intersect = intersect, - find_1st_hit = find_1st_hit, - propagate_to_next_exiting_surface = propagate_to_next_exiting_surface, - ) - def factory(composite): elements = composite.elements() N = len(elements) @@ -53,7 +29,47 @@ def factory(composite): rt_methods = _createRayTracingMethods(composite) find_1st_hit = rt_methods['find_1st_hit'] locate = rt_methods['locate'] + intersect = rt_methods['intersect'] + propagate_out = rt_methods['propagate_out'] propagate_to_next_exiting_surface = rt_methods['propagate_to_next_exiting_surface'] + is_exiting = rt_methods["is_exiting"] + is_exiting_union = rt_methods["is_exiting_union"] + propagate_to_next_incident_surface_union = rt_methods["propagate_to_next_incident_surface_union"] + + @cuda.jit(device=True) + def _scatter(threadindex, rng_states, neutron, tmp_neutron, tmp_ts): + if not is_moving(neutron): + absorb(neutron) + return + if is_exiting(neutron): + return + if is_exiting_union(neutron): + propagate_to_next_exiting_surface(neutron) + return + while True: + itype = _interact_path1(threadindex, rng_states, neutron, tmp_neutron) + if itype == absorption: + return + if itype == scattering: + x,y,z, vx,vy,vz = neutron[:6] + Nintersect = intersect(x,y,z, vx,vy,vz, tmp_ts) + if Nintersect == 0: + propagate_out(neutron) + return + clone(neutron, tmp_neutron) + propagate_out(neutron) + att = calculate_attenuation(tmp_neutron, neutron[:3]) + neutron[-1] *= att + return + # itype == none + if is_exiting(neutron): + return + if is_exiting_union(neutron): + propagate_to_next_exiting_surface(neutron) + return + propagate_to_next_incident_surface_union(neutron) + return + @cuda.jit(device=True) def _interact_path1(threadindex, rng_states, neutron, tmp_neutron): x,y,z = neutron[:3] @@ -99,6 +115,7 @@ def _calculate_attenuation(neutron, end, tmp_neutron): ret *= element_calculate_attenuation(neutron, end, i) return ret + from ..geometry.arrow_intersect import max_intersections if test.USE_CUDASIM: @cuda.jit(device=True, inline=True) def interact_path1(threadindex, rng_states, neutron): @@ -108,6 +125,11 @@ def interact_path1(threadindex, rng_states, neutron): def calculate_attenuation(neutron, end): tmp_neutron = np.zeros(10, dtype=float) return _calculate_attenuation(neutron, end, tmp_neutron) + @cuda.jit(device=True, inline=True) + def scatter(threadindex, rng_states, neutron): + tmp_neutron = np.zeros(10, dtype=float) + tmp_ts = np.zeros(max_intersections, dtype=float) + return _scatter(threadindex, rng_states, neutron, tmp_neutron, tmp_ts) else: @cuda.jit(device=True, inline=True) def interact_path1(threadindex, rng_states, neutron): @@ -117,12 +139,51 @@ def interact_path1(threadindex, rng_states, neutron): def calculate_attenuation(neutron, end): tmp_neutron = cuda.local.array(10, dtype=numba.float64) return _calculate_attenuation(neutron, end, tmp_neutron) + @cuda.jit(device=True, inline=True) + def scatter(threadindex, rng_states, neutron): + tmp_neutron = cuda.local.array(10, dtype=numba.float64) + tmp_ts = cuda.local.array(max_intersections, dtype=numba.float64) + return _scatter(threadindex, rng_states, neutron, tmp_neutron, tmp_ts) return dict( + scatter = scatter, interact_path1 = interact_path1, calculate_attenuation = calculate_attenuation, ) +class _Union: + def identify(self, visitor): return visitor.onUnion(self) + +def _createRayTracingMethods(composite): + """create ray-tracing methods to deal with the overall shape of the composite + """ + # methods regarding the overall shape + shape = composite.shape() + from ..geometry import arrow_intersect + intersect = arrow_intersect.arrow_intersect_func_factory.render(shape) + locate = arrow_intersect.locate_func_factory.render(shape) + del shape + from ..geometry.propagation import makePropagateMethods + propagate_methods = makePropagateMethods(intersect, locate) + ret = propagate_methods + # methods regarding union of element shapes + elements = composite.elements() + shapes = [e.shape() for e in elements] + union = _Union(); union.shapes = shapes + u_intersect = arrow_intersect.arrow_intersect_func_factory.render(union) + u_locate = arrow_intersect.locate_func_factory.render(union) + u_propagate_methods = makePropagateMethods(u_intersect, u_locate) + # find_1st_hit + from mcvine.acc.geometry.composite import get_find_1st_hit + ret.update( + locate = locate, + intersect = intersect, + find_1st_hit = get_find_1st_hit(shapes), + is_exiting_union = u_propagate_methods['is_exiting'], + propagate_to_next_incident_surface_union = u_propagate_methods['propagate_to_next_incident_surface'], + ) + return ret + def _importModule(N): mod = _makeModule(N) import imp diff --git a/tests/scatterers/test_composite_scatterer.py b/tests/scatterers/test_composite_scatterer.py index d1dc2ffc..dc0a7f7f 100644 --- a/tests/scatterers/test_composite_scatterer.py +++ b/tests/scatterers/test_composite_scatterer.py @@ -6,16 +6,23 @@ from mcvine.acc import test @pytest.fixture -def composite_example_interact_path1(): +def composite_example_methods(): path = os.path.join(thisdir, "sampleassemblies", 'sample+2cylinders', 'sampleassembly.xml') from mcvine.acc.components.samples import loadScattererComposite composite = loadScattererComposite(path) from mcvine.acc.scatterers.composite_scatterer import factory - methods = factory(composite) - return methods['interact_path1'] + return factory(composite) + +@pytest.fixture +def composite_example_interact_path1(composite_example_methods): + return composite_example_methods['interact_path1'] + +@pytest.fixture +def composite_example_scatter(composite_example_methods): + return composite_example_methods['scatter'] @pytest.mark.skipif(not test.USE_CUDASIM, reason='no CUDASIM') -def test_cudasim(composite_example_interact_path1): +def test_cudasim_interact_path1(composite_example_interact_path1): neutron = np.array([0.0,0,0, 0,0,1000, 0,0, 0, 1.]) threadindex = 0 rng_states = None @@ -38,8 +45,32 @@ def test_cudasim(composite_example_interact_path1): assert absorbed == 0 return +@pytest.mark.skipif(not test.USE_CUDASIM, reason='no CUDASIM') +def test_cudasim_scatter(composite_example_scatter): + neutron = np.array([0.0,0,0, 0,0,1000, 0,0, 0, 1.]) + threadindex = 0 + rng_states = None + N = 10 #000 + p = 0. + absorbed = 0 + for i in range(N): + neutron = np.array([0.0,0,0, 0,0,1000, 0,0, 0, 1.]) + composite_example_scatter(threadindex, rng_states, neutron) + print(neutron) + p1 = neutron[-1] + if p1>0: + p+=p1 + else: + absorbed += 1 + p/=N + absorbed/=N + print(p, absorbed) + assert p < 0.1 + assert absorbed == 0 + return + @pytest.mark.skipif(not test.USE_CUDA, reason='no CUDA') -def test_cuda(composite_example_interact_path1): +def test_cuda_interact_path1(composite_example_interact_path1): @cuda.jit def interact_path1_kernel(rng_states, neutrons, n_neutrons_per_thread): N = len(neutrons) @@ -67,6 +98,35 @@ def run(neutrons): print(neutrons) return +@pytest.mark.skipif(not test.USE_CUDA, reason='no CUDA') +def test_cuda_scatter(composite_example_scatter): + @cuda.jit + def scatter_kernel(rng_states, neutrons, n_neutrons_per_thread): + N = len(neutrons) + thread_index = cuda.grid(1) + start_index = thread_index*n_neutrons_per_thread + end_index = min(start_index+n_neutrons_per_thread, N) + for i in range(start_index, end_index): + composite_example_scatter(thread_index, rng_states, neutrons[i]) + def run(neutrons): + N = len(neutrons) + threads_per_block = 128 + ntotthreads = int(1e5) + nblocks = math.ceil(ntotthreads / threads_per_block) + actual_nthreads = threads_per_block * nblocks + n_neutrons_per_thread = math.ceil(N / actual_nthreads) + from numba.cuda.random import create_xoroshiro128p_states + from mcvine.acc.config import rng_seed + rng_states = create_xoroshiro128p_states(actual_nthreads, seed=rng_seed) + scatter_kernel[nblocks, threads_per_block](rng_states, neutrons, n_neutrons_per_thread) + N = int(10) + neutrons = np.zeros((N, 10)) + neutron = np.array([0.0,0,0, 0,0,1000, 0,0, 0, 1.]) + neutrons[:] = neutron + run(neutrons) + print(neutrons) + return + def test_makeModule(): from mcvine.acc.scatterers.composite_scatterer import _makeModule _makeModule(4, overwrite=True) From dafc8bb93ae5b8fdad6cf7fa198cd60b1f97d54a Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Sun, 11 Jun 2023 08:24:59 -0700 Subject: [PATCH 46/53] added composite sample component and started writing a test module --- acc/components/samples/composite.py | 34 +++++++++ acc/geometry/arrow_intersect.py | 10 ++- acc/geometry/locate.py | 10 ++- acc/scatterers/__init__.py | 2 - .../test_composite_isotropic_sphere.py | 73 +++++++++++++++++++ ...est_instrument_source_composite_monitor.py | 17 +++++ 6 files changed, 138 insertions(+), 8 deletions(-) create mode 100644 acc/components/samples/composite.py create mode 100755 tests/components/samples/test_composite_isotropic_sphere.py create mode 100644 tests/components/samples/test_instrument_source_composite_monitor.py diff --git a/acc/components/samples/composite.py b/acc/components/samples/composite.py new file mode 100644 index 00000000..e3baf935 --- /dev/null +++ b/acc/components/samples/composite.py @@ -0,0 +1,34 @@ +import os +import numpy as np +from numba import cuda, void, int64 +from numba.cuda.random import xoroshiro128p_uniform_float32, xoroshiro128p_type +from ...config import get_numba_floattype +NB_FLOAT = get_numba_floattype() + +from .SampleBase import SampleBase + +category = 'samples' + +def factory(composite): + from ...scatterers import scatter_func_factory + methods = scatter_func_factory.render(composite) + scatter = methods['scatter'] + class Composite(SampleBase): + + def __init__(self, name): + self.name = name + self.propagate_params = () + + # Aim neutrons toward the sample to cause JIT compilation. + import mcni + neutrons = mcni.neutron_buffer(1) + neutrons[0] = mcni.neutron(r=(0, 0, -1), v=(0, 0, 1), prob=1, time=0) + self.process(neutrons) + + @cuda.jit( + void(int64, xoroshiro128p_type[:], NB_FLOAT[:], + ) , device=True) + def propagate(threadindex, rng_states, neutron): + return scatter(threadindex, rng_states, neutron) + + return Composite diff --git a/acc/geometry/arrow_intersect.py b/acc/geometry/arrow_intersect.py index 51d9e68d..c6860082 100644 --- a/acc/geometry/arrow_intersect.py +++ b/acc/geometry/arrow_intersect.py @@ -27,11 +27,15 @@ def render(self, shape): def onUnion(self, u): nelements = len(u.shapes) - if nelements == 2: return self.onUnion2(u) - if nelements == 3: + if nelements == 1: + return u.shapes[0].identify(self) + elif nelements == 2: + return self.onUnion2(u) + elif nelements == 3: from .composite_3 import createRayTracingMethods_NonOverlappingShapes as createMethods return createMethods(u.shapes)['intersect_all'] - raise NotImplementedError(f"locate for union of {nelements} elements") + else: + raise NotImplementedError(f"locate for union of {nelements} elements") def onUnion2(self, u): locate1 = self.locate_func_factory.onUnion(u) diff --git a/acc/geometry/locate.py b/acc/geometry/locate.py index 5d249f0f..0d757ab8 100644 --- a/acc/geometry/locate.py +++ b/acc/geometry/locate.py @@ -13,11 +13,15 @@ def render(self, shape): def onUnion(self, u): nelements = len(u.shapes) - if nelements == 2: return self.onUnion2(u) - if nelements == 3: + if nelements == 1: + return u.shapes[0].identify(self) + elif nelements == 2: + return self.onUnion2(u) + elif nelements == 3: from .composite_3 import createUnionLocateMethod return createUnionLocateMethod(u.shapes) - raise NotImplementedError(f"locate for union of {nelements} elements") + else: + raise NotImplementedError(f"locate for union of {nelements} elements") def onUnion2(self, u): s1, s2 = u.shapes diff --git a/acc/scatterers/__init__.py b/acc/scatterers/__init__.py index 68f3f09d..bbdf9d1f 100644 --- a/acc/scatterers/__init__.py +++ b/acc/scatterers/__init__.py @@ -11,8 +11,6 @@ def render(self, scatterer): def onCompositeScatterer(self, composite): elements = composite.elements() - if len(elements) != 3: - raise NotImplementedError from .composite_scatterer import factory return factory(composite) diff --git a/tests/components/samples/test_composite_isotropic_sphere.py b/tests/components/samples/test_composite_isotropic_sphere.py new file mode 100755 index 00000000..cfd1103f --- /dev/null +++ b/tests/components/samples/test_composite_isotropic_sphere.py @@ -0,0 +1,73 @@ +#!/usr/bin/env python + +import os, shutil +import pytest +from mcvine.acc import test +from mcvine import run_script + +thisdir = os.path.dirname(__file__) + + +def test_1(): + instr = os.path.join(thisdir, "test_instrument_source_composite_monitor.py") + samplexml = "" + outdir = 'out.test_instrument_source_composite_monitor_isotropic_sphere' + if os.path.exists(outdir): shutil.rmtree(outdir) + ncount = 1e5 + run_script.run1( + instr, outdir, + ncount=ncount, + samplexml = "sampleassemblies/isotropic_sphere/sampleassembly.xml" + ) + return + +@pytest.mark.skipif(not test.USE_CUDA, reason='No CUDA') +def test_compare_mcvine(num_neutrons=int(1e7), debug=False, interactive=False): + """ + Tests the acc cpu implementation of a straight guide against mcvine + """ + instr = os.path.join(thisdir, "isotropic_sphere_instrument.py") + from mcvine.acc.test.compare_acc_nonacc import compare_acc_nonacc + compare_acc_nonacc( + "isotropic_sphere", + ["psd_4pi"], + {"float32": 4e-10, "float64": 4e-10}, + num_neutrons, debug, + instr = instr, + interactive=interactive, + acc_component_spec = dict(is_acc=True), + nonacc_component_spec = dict(is_acc=False), + ) + +def psd_monitor_4pi(): + from mcvine.acc.components.monitors.psd_monitor_4pi import PSD_monitor_4Pi + return PSD_monitor_4Pi( + "mon", + nphi=30, ntheta=30, radius=3, + filename = "psd_4pi.h5", + ) + +@pytest.mark.skipif(not test.USE_CUDA, reason='No CUDA') +def test_acc_run_script(ncount = 1e6): + instr = os.path.join(thisdir, "acc_ss_instrument.py") + outdir = 'out.isotropic_sphere-acc_run_script' + ncount = int(ncount) + from mcvine.acc import run_script + run_script.run( + instr, outdir, ncount=ncount, + monitor_factory=psd_monitor_4pi, + z_sample = 1., + ) + return + +def main(): + import journal + journal.info("instrument").activate() + test_1() + # test_compare_mcvine(num_neutrons=int(100), interactive=True, debug=True) + # test_compare_mcvine(num_neutrons=int(1e7), interactive=True) + # test_acc_run_script(ncount=1e7) + return + + +if __name__ == '__main__': main() diff --git a/tests/components/samples/test_instrument_source_composite_monitor.py b/tests/components/samples/test_instrument_source_composite_monitor.py new file mode 100644 index 00000000..74193127 --- /dev/null +++ b/tests/components/samples/test_instrument_source_composite_monitor.py @@ -0,0 +1,17 @@ +#!/usr/bin/env python + +import os +thisdir = os.path.dirname(__file__) + +from test_sample_instrument_factory import construct + +def instrument(samplexml, is_acc=True): + if is_acc: + from mcvine.acc.components.samples import loadScattererComposite + composite = loadScattererComposite(samplexml) + from mcvine.acc.components.samples.composite import factory + target = factory(composite)("sample") + else: + import mcvine.components as mc + target = mc.samples.SampleAssemblyFromXml("sample", samplexml) + return construct(target, size=0.) From 883423c5dcb4dc2b72c4ff8cbe064f8d9a6ca13c Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Sun, 11 Jun 2023 09:14:48 -0700 Subject: [PATCH 47/53] part of an effort to simplify test instruments. see #111 --- acc/components/samples/composite.py | 6 ++++++ .../samples/test_composite_isotropic_sphere.py | 3 ++- .../test_instrument_source_composite_monitor.py | 17 ++++++++--------- 3 files changed, 16 insertions(+), 10 deletions(-) diff --git a/acc/components/samples/composite.py b/acc/components/samples/composite.py index e3baf935..310578be 100644 --- a/acc/components/samples/composite.py +++ b/acc/components/samples/composite.py @@ -9,6 +9,12 @@ category = 'samples' +def sampleassembly_from_xml(name, samplexml): + from . import loadScattererComposite + composite = loadScattererComposite(samplexml) + component = factory(composite) + return component(name) + def factory(composite): from ...scatterers import scatter_func_factory methods = scatter_func_factory.render(composite) diff --git a/tests/components/samples/test_composite_isotropic_sphere.py b/tests/components/samples/test_composite_isotropic_sphere.py index cfd1103f..1862c87a 100755 --- a/tests/components/samples/test_composite_isotropic_sphere.py +++ b/tests/components/samples/test_composite_isotropic_sphere.py @@ -17,7 +17,8 @@ def test_1(): run_script.run1( instr, outdir, ncount=ncount, - samplexml = "sampleassemblies/isotropic_sphere/sampleassembly.xml" + samplexml = "sampleassemblies/isotropic_sphere/sampleassembly.xml", + factory = "mcvine.acc.components.samples.composite.sampleassembly_from_xml", ) return diff --git a/tests/components/samples/test_instrument_source_composite_monitor.py b/tests/components/samples/test_instrument_source_composite_monitor.py index 74193127..f32e3037 100644 --- a/tests/components/samples/test_instrument_source_composite_monitor.py +++ b/tests/components/samples/test_instrument_source_composite_monitor.py @@ -5,13 +5,12 @@ from test_sample_instrument_factory import construct -def instrument(samplexml, is_acc=True): - if is_acc: - from mcvine.acc.components.samples import loadScattererComposite - composite = loadScattererComposite(samplexml) - from mcvine.acc.components.samples.composite import factory - target = factory(composite)("sample") - else: - import mcvine.components as mc - target = mc.samples.SampleAssemblyFromXml("sample", samplexml) +def instrument(samplexml, factory): + tokens = factory.split('.') + module = '.'.join(tokens[:-1]) + method = tokens[-1] + import importlib + module = importlib.import_module(module) + factory = getattr(module, method) + target = factory(name="sample", samplexml=samplexml) return construct(target, size=0.) From 3f65d73e0c9fc2fa241561469c6d55b5a3877a32 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Sun, 11 Jun 2023 09:29:55 -0700 Subject: [PATCH 48/53] refactoring: tests/components/samples/test_composite_isotropic_sphere.py --- .../test_composite_isotropic_sphere.py | 19 +++++++++++++------ ... => test_instrument_with_samplefromxml.py} | 2 +- 2 files changed, 14 insertions(+), 7 deletions(-) rename tests/components/samples/{test_instrument_source_composite_monitor.py => test_instrument_with_samplefromxml.py} (82%) diff --git a/tests/components/samples/test_composite_isotropic_sphere.py b/tests/components/samples/test_composite_isotropic_sphere.py index 1862c87a..79ed7416 100755 --- a/tests/components/samples/test_composite_isotropic_sphere.py +++ b/tests/components/samples/test_composite_isotropic_sphere.py @@ -9,15 +9,15 @@ def test_1(): - instr = os.path.join(thisdir, "test_instrument_source_composite_monitor.py") - samplexml = "" + instr = os.path.join(thisdir, "test_instrument_with_samplefromxml.py") + samplexml = "sampleassemblies/isotropic_sphere/sampleassembly.xml" outdir = 'out.test_instrument_source_composite_monitor_isotropic_sphere' if os.path.exists(outdir): shutil.rmtree(outdir) ncount = 1e5 run_script.run1( instr, outdir, ncount=ncount, - samplexml = "sampleassemblies/isotropic_sphere/sampleassembly.xml", + samplexml = samplexml, factory = "mcvine.acc.components.samples.composite.sampleassembly_from_xml", ) return @@ -27,7 +27,8 @@ def test_compare_mcvine(num_neutrons=int(1e7), debug=False, interactive=False): """ Tests the acc cpu implementation of a straight guide against mcvine """ - instr = os.path.join(thisdir, "isotropic_sphere_instrument.py") + instr = os.path.join(thisdir, "test_instrument_with_samplefromxml.py") + samplexml = "sampleassemblies/isotropic_sphere/sampleassembly.xml", from mcvine.acc.test.compare_acc_nonacc import compare_acc_nonacc compare_acc_nonacc( "isotropic_sphere", @@ -36,8 +37,14 @@ def test_compare_mcvine(num_neutrons=int(1e7), debug=False, interactive=False): num_neutrons, debug, instr = instr, interactive=interactive, - acc_component_spec = dict(is_acc=True), - nonacc_component_spec = dict(is_acc=False), + acc_component_spec = dict( + samplexml=samplexml, + factory = "mcvine.acc.components.samples.composite.sampleassembly_from_xml", + ), + nonacc_component_spec = dict( + samplexml=samplexml, + factory="mcvine.components.samples.SampleAssemblyFromXml", + ), ) def psd_monitor_4pi(): diff --git a/tests/components/samples/test_instrument_source_composite_monitor.py b/tests/components/samples/test_instrument_with_samplefromxml.py similarity index 82% rename from tests/components/samples/test_instrument_source_composite_monitor.py rename to tests/components/samples/test_instrument_with_samplefromxml.py index f32e3037..4d5cd341 100644 --- a/tests/components/samples/test_instrument_source_composite_monitor.py +++ b/tests/components/samples/test_instrument_with_samplefromxml.py @@ -5,7 +5,7 @@ from test_sample_instrument_factory import construct -def instrument(samplexml, factory): +def instrument(samplexml, factory="mcvine.components.samples.SampleAssemblyFromXml"): tokens = factory.split('.') module = '.'.join(tokens[:-1]) method = tokens[-1] From 7725e57575c6867883726231e08b72045fc9b2e7 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Wed, 14 Jun 2023 22:57:51 -0700 Subject: [PATCH 49/53] scatter method - match signature to component.propagate: acc/scatterers/composite_scatterer.py --- acc/scatterers/composite_scatterer.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/acc/scatterers/composite_scatterer.py b/acc/scatterers/composite_scatterer.py index 93db9194..b1ed5be2 100644 --- a/acc/scatterers/composite_scatterer.py +++ b/acc/scatterers/composite_scatterer.py @@ -6,8 +6,9 @@ import os from math import sqrt, exp import numpy as np, numba -from numba import cuda +from numba import cuda, void, int64 from numba.core import config +from numba.cuda.random import xoroshiro128p_type from .. import test from .._numba import xoroshiro128p_uniform_float32 @@ -139,7 +140,10 @@ def interact_path1(threadindex, rng_states, neutron): def calculate_attenuation(neutron, end): tmp_neutron = cuda.local.array(10, dtype=numba.float64) return _calculate_attenuation(neutron, end, tmp_neutron) - @cuda.jit(device=True, inline=True) + @cuda.jit( + void(int64, xoroshiro128p_type[:], NB_FLOAT[:]), + device=True, inline=True + ) def scatter(threadindex, rng_states, neutron): tmp_neutron = cuda.local.array(10, dtype=numba.float64) tmp_ts = cuda.local.array(max_intersections, dtype=numba.float64) @@ -290,5 +294,3 @@ def element_calculate_attenuation(self): ] ) return header + loop + [indent + 'return 1.0'] - - From 235529f2402b2fd574be4021493d189ddbb1c967 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Wed, 14 Jun 2023 23:04:14 -0700 Subject: [PATCH 50/53] composite.py: 1.added SampleAssemblyFromXml method to create the component class; 2.directly register scatter method created by mcvine.acc.scatterers.scatter_func_factory as the propagate method --- acc/components/samples/composite.py | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/acc/components/samples/composite.py b/acc/components/samples/composite.py index 310578be..83b4bb90 100644 --- a/acc/components/samples/composite.py +++ b/acc/components/samples/composite.py @@ -1,5 +1,5 @@ import os -import numpy as np +import numpy as np, numba from numba import cuda, void, int64 from numba.cuda.random import xoroshiro128p_uniform_float32, xoroshiro128p_type from ...config import get_numba_floattype @@ -9,11 +9,14 @@ category = 'samples' -def sampleassembly_from_xml(name, samplexml): +def SampleAssemblyFromXml(samplexml): from . import loadScattererComposite composite = loadScattererComposite(samplexml) - component = factory(composite) - return component(name) + return factory(composite) + +def sampleassembly_from_xml(name, samplexml): + klass = SampleAssemblyFromXml(samplexml) + return klass(name) def factory(composite): from ...scatterers import scatter_func_factory @@ -32,9 +35,11 @@ def __init__(self, name): self.process(neutrons) @cuda.jit( - void(int64, xoroshiro128p_type[:], NB_FLOAT[:], - ) , device=True) + void(int64, xoroshiro128p_type[:], NB_FLOAT[:]), + device=True, inline=True, + ) def propagate(threadindex, rng_states, neutron): - return scatter(threadindex, rng_states, neutron) - + return + Composite.propagate = scatter + Composite.register_propagate_method(scatter) return Composite From fc7e27eb842f32fc1acb4e7296365204e930ae7a Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Wed, 14 Jun 2023 23:06:59 -0700 Subject: [PATCH 51/53] make sure mcweights and packing_factor are floats: acc/scatterers/__init__.py --- acc/scatterers/__init__.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/acc/scatterers/__init__.py b/acc/scatterers/__init__.py index bbdf9d1f..ec52fa13 100644 --- a/acc/scatterers/__init__.py +++ b/acc/scatterers/__init__.py @@ -18,8 +18,8 @@ def onHomogeneousScatterer(self, hs): from .homogeneous_scatterer import factory shape = hs.shape() kernel = hs.kernel() - mcweights = hs.mcweights - packing_factor = hs.packing_factor + mcweights = [float(w) for w in hs.mcweights] + packing_factor = float(hs.packing_factor) methods = factory(shape, kernel, mcweights, packing_factor) return methods From 2c289b6b2f02a6c2569f861b238901a350c97435 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Tue, 25 Jul 2023 22:40:44 -0700 Subject: [PATCH 52/53] minor (inline doc): components/samples/composite.py --- acc/components/samples/composite.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/acc/components/samples/composite.py b/acc/components/samples/composite.py index 83b4bb90..5ff2a4be 100644 --- a/acc/components/samples/composite.py +++ b/acc/components/samples/composite.py @@ -9,15 +9,18 @@ category = 'samples' +def sampleassembly_from_xml(name, samplexml): + """factory method to create sample assembly component. + can only be used in nonacc instrument script + """ + klass = SampleAssemblyFromXml(samplexml) + return klass(name) + def SampleAssemblyFromXml(samplexml): from . import loadScattererComposite composite = loadScattererComposite(samplexml) return factory(composite) -def sampleassembly_from_xml(name, samplexml): - klass = SampleAssemblyFromXml(samplexml) - return klass(name) - def factory(composite): from ...scatterers import scatter_func_factory methods = scatter_func_factory.render(composite) From 4e2367e6c6d2ed86a0f09af62fc9d09a085d8531 Mon Sep 17 00:00:00 2001 From: Jiao Lin Date: Tue, 25 Jul 2023 22:42:19 -0700 Subject: [PATCH 53/53] debugging: tests/components/samples/test_composite_isotropic_sphere.py --- .../samples/test_composite_isotropic_sphere.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/tests/components/samples/test_composite_isotropic_sphere.py b/tests/components/samples/test_composite_isotropic_sphere.py index 79ed7416..d1780996 100755 --- a/tests/components/samples/test_composite_isotropic_sphere.py +++ b/tests/components/samples/test_composite_isotropic_sphere.py @@ -57,8 +57,9 @@ def psd_monitor_4pi(): @pytest.mark.skipif(not test.USE_CUDA, reason='No CUDA') def test_acc_run_script(ncount = 1e6): - instr = os.path.join(thisdir, "acc_ss_instrument.py") - outdir = 'out.isotropic_sphere-acc_run_script' + # instr = os.path.join(thisdir, "acc_ss_instrument.py") + instr = os.path.join(thisdir, "acc_composite_isotropic_sphere_instrument.py") + outdir = 'out.composite_isotropic_sphere-acc_run_script' ncount = int(ncount) from mcvine.acc import run_script run_script.run( @@ -71,10 +72,10 @@ def test_acc_run_script(ncount = 1e6): def main(): import journal journal.info("instrument").activate() - test_1() + # test_1() # test_compare_mcvine(num_neutrons=int(100), interactive=True, debug=True) # test_compare_mcvine(num_neutrons=int(1e7), interactive=True) - # test_acc_run_script(ncount=1e7) + test_acc_run_script(ncount=1e6) return