diff --git a/acc/kernels/Phonon_IncoherentElastic.py b/acc/kernels/Phonon_IncoherentElastic.py
new file mode 100644
index 0000000..9c30787
--- /dev/null
+++ b/acc/kernels/Phonon_IncoherentElastic.py
@@ -0,0 +1,83 @@
+import numpy as np
+from mcni.utils import conversion
+from mcvine.acc.neutron import v2e
+from numba import cuda, float64
+from numba.cuda.random import xoroshiro128p_uniform_float32, xoroshiro128p_type
+import math
+
+from .. import vec3
+from ..config import get_numba_floattype
+
+NB_FLOAT = get_numba_floattype()
+
+epsilon = 1e-7
+
+
+@cuda.jit(device=True)
+def scatter(threadindex, rng_states, neutron, dw_core):
+ """
+ Scatters to a specific target location (within circle of given radius) and at a specific TOF (within +/- delta_tof / 2)
+ """
+
+ v = neutron[3:6]
+
+ # incident neutron velocity
+ vi = vec3.length(v)
+
+ # randomly pick a theta between [0, PI]
+ r = xoroshiro128p_uniform_float32(rng_states, threadindex)
+ theta = r * math.pi
+
+ # randomly pick a phi between [0, 2*PI]
+ r = xoroshiro128p_uniform_float32(rng_states, threadindex)
+ phi = r * 2.0 * math.pi
+
+ Q = conversion.V2K * vi * 2.0 * math.sin(0.5 * theta)
+
+ # Debye Waller factor
+ dw = math.exp(-dw_core * Q * Q)
+
+ e1 = cuda.local.array(3, dtype=float64)
+ e2 = cuda.local.array(3, dtype=float64)
+ norm = cuda.local.array(3, dtype=float64)
+
+ # Scattered neutron velocity vector
+ vec3.copy(v, e1)
+ vec3.normalize(e1)
+
+ # if e1 is not in the z-direction we set e2 to be the cross product of e1 and (0, 0, 1)
+ # if e1 is right on the z-direction, that means e1 = (0, 0, 1) and we set e2 = (1, 0, 0)
+ if cuda.libdevice.fabs(e1[0]) > epsilon or cuda.libdevice.fabs(e1[1]) > epsilon:
+ #norm = [0.0, 0.0, 1.0]
+ norm[0] = 0.0
+ norm[1] = 0.0
+ norm[2] = 1.0
+ vec3.cross(norm, e1, e2)
+ vec3.normalize(e2)
+ else:
+ #e2 = [1.0, 0.0, 0.0]
+ e2[0] = 1.0
+ e2[1] = 0.0
+ e2[2] = 0.0
+
+ # calculate e3 = e1 * e2
+ vec3.cross(e1, e2, norm) # norm = e3
+
+ sint = cuda.libdevice.sin(theta)
+
+ # V_f = sin(theta) * cos(phi) * e2 +
+ # sin(theta) * sin(phi) * e3 +
+ # cos(theta) * e1
+
+ vec3.scale(e2, sint * cuda.libdevice.cos(phi)) # sin(theta) * cos(phi) * e2
+ vec3.scale(norm, sint * cuda.libdevice.sin(phi)) # sin(theta) * sin(phi) * e3
+ vec3.scale(e1, cuda.libdevice.cos(theta)) # cos(theta) * e1
+
+ # elastic scattering
+ neutron[3] = vi * (e2[0] + norm[0] + e1[0])
+ neutron[4] = vi * (e2[1] + norm[1] + e1[1])
+ neutron[5] = vi * (e2[2] + norm[2] + e1[2])
+
+ # adjust probability of neutron event
+ # normalization factor is 2 * PI^2 / 4PI = PI / 2
+ neutron[-1] *= sint * (math.pi * 0.5) * dw
diff --git a/acc/kernels/__init__.py b/acc/kernels/__init__.py
index bb5016d..38c18ef 100644
--- a/acc/kernels/__init__.py
+++ b/acc/kernels/__init__.py
@@ -144,6 +144,32 @@ def dgssxres_scatter(threadindex, rng_states, neutron):
tof_target, dtof)
return dgssxres_scatter, None, None, None
+ def onPhonon_IncoherentElastic_Kernel(self, kernel):
+ from ..components.samples import getAbsScttCoeffs
+ mu, sigma = getAbsScttCoeffs(kernel)
+
+ print(kernel)
+ print(kernel.scatterer_origin)
+ print(kernel.scatterer_origin.phase.unitcell)
+
+ # additional kernel parameters
+ AA= units.length.angstrom
+ dw_core = kernel.dw_core / AA**2
+
+ # additional kernel parameters
+ scattering_xs = kernel.scattering_xs/units.area.barn \
+ if kernel.scattering_xs else 0.
+ absorption_xs = kernel.absorption_xs/units.area.barn \
+ if kernel.absorption_xs else 0.
+
+ from .Phonon_IncoherentElastic import scatter
+ @cuda.jit(device=True)
+ def phonon_incoherentelastic_scatter(threadindex, rng_states, neutron):
+ neutron[-1] *= sigma
+ return scatter(threadindex, rng_states, neutron, dw_core)
+ return phonon_incoherentelastic_scatter, None, None, None
+
+
scatter_func_factory = ScatterFuncFactory()
def make_calc_sctt_coeff_func(sigma):
diff --git a/tests/components/samples/HSS_fccNi_PhononIncohEl_sphere.py b/tests/components/samples/HSS_fccNi_PhononIncohEl_sphere.py
new file mode 100644
index 0000000..c985d5e
--- /dev/null
+++ b/tests/components/samples/HSS_fccNi_PhononIncohEl_sphere.py
@@ -0,0 +1,11 @@
+import os
+from mcvine.acc.components.samples.homogeneous_single_scatterer import factory
+thisdir = os.path.dirname(__file__)
+
+path = os.path.join(thisdir, "sampleassemblies", 'incoh-el', 'sampleassembly.xml')
+from mcvine.acc.components.samples import loadFirstHomogeneousScatterer
+hs = loadFirstHomogeneousScatterer(path)
+shape = hs.shape()
+kernel = hs.kernel()
+HSSbase = factory(shape = shape, kernel = kernel)
+class HSS(HSSbase): pass
diff --git a/tests/components/samples/fccNi_PhononIncohEl_sphere_instrument.py b/tests/components/samples/fccNi_PhononIncohEl_sphere_instrument.py
new file mode 100644
index 0000000..1012260
--- /dev/null
+++ b/tests/components/samples/fccNi_PhononIncohEl_sphere_instrument.py
@@ -0,0 +1,38 @@
+#!/usr/bin/env python
+
+import os
+thisdir = os.path.dirname(__file__)
+
+from test_sample_instrument_factory import construct, Builder as base
+
+class Builder(base):
+
+ def addIQEMonitor(self, **kwds):
+ params = dict(
+ Ei = 70.0,
+ Qmin=0., Qmax=6.0, nQ = 160,
+ Emin=-.15, Emax=0.15, nE = 120,
+ min_angle_in_plane=0., max_angle_in_plane=359.,
+ min_angle_out_of_plane=-90., max_angle_out_of_plane=90.,
+ )
+ params.update(kwds)
+ mon = self.get_monitor(
+ subtype="IQE_monitor", name = "IQE",
+ **params
+ )
+ self.add(mon, gap=0)
+
+
+def instrument(is_acc=True):
+ if is_acc:
+ from HSS_fccNi_PhononIncohEl_sphere import HSS
+ target = HSS(name='sample')
+ else:
+ import mcvine.components as mc
+ xml=os.path.join(thisdir, "sampleassemblies", "incoh-el", "sampleassembly.xml")
+ target = mc.samples.SampleAssemblyFromXml("sample", xml)
+ source_params = dict(E0 = 70.0, dE=0.1, Lambda0=0, dLambda=0.)
+ return construct(
+ target, size=0.,
+ source_params=source_params, monitors=['IQE'],
+ builder=Builder())
diff --git a/tests/components/samples/sampleassemblies/incoh-el/Ni.xyz b/tests/components/samples/sampleassemblies/incoh-el/Ni.xyz
new file mode 100644
index 0000000..fd0e4f2
--- /dev/null
+++ b/tests/components/samples/sampleassemblies/incoh-el/Ni.xyz
@@ -0,0 +1,3 @@
+1
+1.76 1.76 0 1.76 0 1.76 0 -1.76 -1.76
+Ni 0 0 0
diff --git a/tests/components/samples/sampleassemblies/incoh-el/fccNi-sphere-scatterer.xml b/tests/components/samples/sampleassemblies/incoh-el/fccNi-sphere-scatterer.xml
new file mode 100644
index 0000000..1150119
--- /dev/null
+++ b/tests/components/samples/sampleassemblies/incoh-el/fccNi-sphere-scatterer.xml
@@ -0,0 +1,27 @@
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/tests/components/samples/sampleassemblies/incoh-el/sampleassembly.xml b/tests/components/samples/sampleassemblies/incoh-el/sampleassembly.xml
new file mode 100644
index 0000000..a95a449
--- /dev/null
+++ b/tests/components/samples/sampleassemblies/incoh-el/sampleassembly.xml
@@ -0,0 +1,27 @@
+
+
+
+
+
+
+
+
+
+
+
+ Ni
+ Ni.xyz
+
+
+
+
+
+
+
+
+
diff --git a/tests/components/samples/test_fccNi_Phonon_IncoherentElastic.py b/tests/components/samples/test_fccNi_Phonon_IncoherentElastic.py
new file mode 100755
index 0000000..5260638
--- /dev/null
+++ b/tests/components/samples/test_fccNi_Phonon_IncoherentElastic.py
@@ -0,0 +1,81 @@
+#!/usr/bin/env python
+
+import os
+import shutil
+import pytest
+import numpy as np
+from mcvine.acc import test
+from mcvine import run_script
+
+thisdir = os.path.dirname(__file__)
+
+
+def test_cpu():
+ instr = os.path.join(thisdir, "fccNi_PhononIncohEl_sphere_instrument.py")
+ outdir = 'out.fccNi_PhononIncohEl_sphere'
+ if os.path.exists(outdir):
+ shutil.rmtree(outdir)
+ ncount = 1e5
+ run_script.run1(
+ instr, outdir,
+ ncount=ncount, buffer_size=int(ncount),
+ is_acc=False,
+ )
+ return
+
+
+@pytest.mark.skipif(not test.USE_CUDA, reason='No CUDA')
+def test_compare_mcvine(num_neutrons=int(1024), debug=False, interactive=False):
+ """
+ Tests the acc cpu implementation of a straight guide against mcvine
+ """
+ instr = os.path.join(thisdir, "fccNi_PhononIncohEl_sphere_instrument.py")
+ from mcvine.acc.test.compare_acc_nonacc import compare_acc_nonacc
+ compare_acc_nonacc(
+ "fccNi_PhononIncohEl_sphere",
+ ["IQE"],
+ {"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),
+ )
+
+ gpu_file = os.path.join("./out.debug-fccni_phononincohel_sphere_gpu_instrument", "IQE.h5")
+ cpu_file = os.path.join("./out.debug-mcvine_fccni_phononincohel_sphere_cpu_instrument", "IQE.h5")
+
+ import histogram.hdf as hh
+ gpu_hist = hh.load(gpu_file)
+ cpu_hist = hh.load(cpu_file)
+
+ # integrate the energies
+ gpu_E_hist = gpu_hist.sum("energy")
+ cpu_E_hist = cpu_hist.sum("energy")
+
+ # subtract the integrated CPU and GPU intensities
+ diff_hist = gpu_E_hist - cpu_E_hist
+
+ if interactive:
+ from histogram import plot as plotHist
+ plotHist(gpu_hist)
+ plotHist(cpu_hist)
+
+ plotHist(gpu_E_hist)
+ plotHist(cpu_E_hist)
+ plotHist(diff_hist)
+
+ assert np.mean(diff_hist.I) < 1.0e-12
+
+
+def main():
+ import journal
+ journal.info("instrument").activate()
+ # test_cpu()
+ # test_compare_mcvine(num_neutrons=int(100), interactive=True, debug=True)
+ test_compare_mcvine(num_neutrons=int(1e6), interactive=True)
+ return
+
+
+if __name__ == '__main__':
+ main()
diff --git a/tests/kernels/test_Phonon_IncoherentElastic.py b/tests/kernels/test_Phonon_IncoherentElastic.py
new file mode 100644
index 0000000..6bd0c1b
--- /dev/null
+++ b/tests/kernels/test_Phonon_IncoherentElastic.py
@@ -0,0 +1,65 @@
+import numpy as np
+import pytest
+from mcni import neutron_buffer, neutron
+from mcni.neutron_storage import neutrons_as_npyarr, ndblsperneutron
+from mcni.utils import conversion
+from mcvine.acc import test
+from mcvine.acc.config import rng_seed
+from mcvine.acc.kernels import Phonon_IncoherentElastic
+from numba import cuda
+from numba.cuda.random import create_xoroshiro128p_states
+
+
+# Simple wrapper kernel to test the scatter device function
+@cuda.jit()
+def scatter_test_kernel(
+ rng_states, N, n_neutrons_per_thread, neutrons, dw_core
+):
+ thread_index = cuda.grid(1)
+ start_index = thread_index * n_neutrons_per_thread
+ end_index = min(start_index + n_neutrons_per_thread, N)
+ for neutron_index in range(start_index, end_index):
+ Phonon_IncoherentElastic.scatter(
+ neutron_index, rng_states, neutrons[neutron_index], dw_core)
+
+
+@pytest.mark.skipif(not test.USE_CUDA, reason='No CUDA')
+def test_IncoherentElastic_kernel():
+
+ tof_at_sample = 1.0
+ dw_core = 0.1
+
+ vil = 3000.0
+ n = neutron([0., 0., -5.0], [0., 0., vil], prob=1.0, time=tof_at_sample)
+
+ # calculate initial vi and ei
+ vi = np.linalg.norm(n.state.velocity)
+ Ei = conversion.v2e(vi)
+
+ # create neutron buffer
+ ncount = 10
+ buffer = neutron_buffer(ncount)
+ for i in range(ncount):
+ buffer[i] = n
+ tmp = neutrons_as_npyarr(buffer)
+ tmp.shape = -1, ndblsperneutron
+ buffer_d = cuda.to_device(tmp)
+
+ # setup test kernel with 1 neutron
+ nblocks = ncount
+ threads_per_block = 1
+ rng_states = create_xoroshiro128p_states(
+ nblocks * threads_per_block, seed=rng_seed)
+
+ scatter_test_kernel[nblocks, threads_per_block](
+ rng_states, ncount, 1, buffer_d, dw_core)
+
+ cuda.synchronize()
+ buffer = buffer_d.copy_to_host()
+
+ # calculate Q and E and compare against expected
+ print(buffer)
+
+
+if __name__ == '__main__':
+ pytest.main([__file__])