Phonon Incoherent Elastic sample kernel

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

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):

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

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())

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

tests/components/samples/sampleassemblies/incoh-el/fccNi-sphere-scatterer.xml

@@ -0,0 +1,27 @@
+<?xml version="1.0"?>
+
+<!DOCTYPE scatterer>
+
+<!-- weights: absorption, scattering, transmission -->
+<homogeneous_scatterer mcweights="0, 1, 0">
+
+  <!--
+  IncoherentElastic kernel
+  dw_core: could be computed from DOS
+  scattering_xs and absorption_xs: optional. if not set, they will
+    be computed from unitcell info in Ni.xyz
+    here we set scattering_xs to 18.5 barn, the total scattering cross
+    section that is the sum of incoherent and coherent scattering cross
+    sections. This way, we can use the incoherent kernel to approximate
+    the total scattering.
+    Note: be very careful when customizing cross sections.
+      They are the total scattering cross sections from
+      the whole unit cell.
+  -->
+  <Phonon_IncoherentElastic_Kernel 
+      dw_core='0.1*angstrom**2' 
+      scattering_xs="18.5*barn"
+      >
+  </Phonon_IncoherentElastic_Kernel>
+
+</homogeneous_scatterer>

tests/components/samples/sampleassemblies/incoh-el/sampleassembly.xml

@@ -0,0 +1,27 @@
+<?xml version="1.0"?>
+
+<!DOCTYPE SampleAssembly>
+
+<SampleAssembly name="fccNi-sphere-sampleassembly">
+
+  <PowderSample name='fccNi-sphere' type='sample'>
+    <Shape>
+      <sphere radius="1*cm" />
+    </Shape>
+    <Phase type="crystal">
+      <ChemicalFormula>Ni</ChemicalFormula>
+      <xyzfile>Ni.xyz</xyzfile>
+    </Phase>
+  </PowderSample>
+
+  <LocalGeometer registry-coordinate-system="InstrumentScientist">
+    <Register name="fccNi-sphere" position="(0,0,0)" orientation="(0,0,0)"/>
+  </LocalGeometer>
+<!--
+  <environment>
+    <temperature>300*K</temperature>
+    <pressure>1.*atm</pressure>
+  </environment>
+-->
+
+</SampleAssembly>

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()

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__])