From 347a967fa9f49f70a8623714f3ea36fc3412209a Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Fri, 24 Mar 2023 13:32:16 -0400 Subject: [PATCH 1/6] Add test sample assembly for phonon incoh elastic kernel --- .../samples/HSS_fccNi_PhononIncohEl_plate.py | 11 ++++++++ .../fccNi_PhononIncohEl_plate_instrument.py | 20 ++++++++++++++ .../samples/sampleassemblies/incoh-el/Ni.xyz | 3 +++ .../incoh-el/fccNi-plate-scatterer.xml | 27 +++++++++++++++++++ .../incoh-el/sampleassembly.xml | 27 +++++++++++++++++++ 5 files changed, 88 insertions(+) create mode 100644 tests/components/samples/HSS_fccNi_PhononIncohEl_plate.py create mode 100644 tests/components/samples/fccNi_PhononIncohEl_plate_instrument.py create mode 100644 tests/components/samples/sampleassemblies/incoh-el/Ni.xyz create mode 100644 tests/components/samples/sampleassemblies/incoh-el/fccNi-plate-scatterer.xml create mode 100644 tests/components/samples/sampleassemblies/incoh-el/sampleassembly.xml diff --git a/tests/components/samples/HSS_fccNi_PhononIncohEl_plate.py b/tests/components/samples/HSS_fccNi_PhononIncohEl_plate.py new file mode 100644 index 0000000..c985d5e --- /dev/null +++ b/tests/components/samples/HSS_fccNi_PhononIncohEl_plate.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_plate_instrument.py b/tests/components/samples/fccNi_PhononIncohEl_plate_instrument.py new file mode 100644 index 0000000..39580fa --- /dev/null +++ b/tests/components/samples/fccNi_PhononIncohEl_plate_instrument.py @@ -0,0 +1,20 @@ +#!/usr/bin/env python + +import os +thisdir = os.path.dirname(__file__) + +from test_powderdiffraction_sample_instrument_factory import construct, Builder + +def instrument(is_acc=True): + if is_acc: + from HSS_fccNi_PhononIncohEl_plate 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=['PSD_4PI'], + 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-plate-scatterer.xml b/tests/components/samples/sampleassemblies/incoh-el/fccNi-plate-scatterer.xml new file mode 100644 index 0000000..1150119 --- /dev/null +++ b/tests/components/samples/sampleassemblies/incoh-el/fccNi-plate-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..2caa5c6 --- /dev/null +++ b/tests/components/samples/sampleassemblies/incoh-el/sampleassembly.xml @@ -0,0 +1,27 @@ + + + + + + + + + + + + Ni + Ni.xyz + + + + + + + + + From 8c96131691b1655b1e5a0704179cc2e2b8ec4cbb Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Thu, 10 Aug 2023 15:49:18 -0400 Subject: [PATCH 2/6] Add initial PhononIncoh Elastic sample kernel --- acc/kernels/Phonon_IncoherentElastic.py | 83 +++++++++++++++++++++++++ acc/kernels/__init__.py | 26 ++++++++ 2 files changed, 109 insertions(+) create mode 100644 acc/kernels/Phonon_IncoherentElastic.py 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..78f3464 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.angstrom + dw_core = kernel.dw_core / AA**2 + + # additional kernel parameters + scattering_xs = kernel.scattering_xs/units.barn \ + if kernel.scattering_xs else 0. + absorption_xs = kernel.absorption_xs/units.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): From 4ad5311db82b05b7c34f72419086bba568bb9265 Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Fri, 11 Aug 2023 11:38:04 -0400 Subject: [PATCH 3/6] Add test files for PhononIncoh Elastic kernel --- .../test_fccNi_Phonon_IncoherentElastic.py | 55 ++++++++++++++++ .../kernels/test_Phonon_IncoherentElastic.py | 65 +++++++++++++++++++ 2 files changed, 120 insertions(+) create mode 100755 tests/components/samples/test_fccNi_Phonon_IncoherentElastic.py create mode 100644 tests/kernels/test_Phonon_IncoherentElastic.py 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..7d2018d --- /dev/null +++ b/tests/components/samples/test_fccNi_Phonon_IncoherentElastic.py @@ -0,0 +1,55 @@ +#!/usr/bin/env python + +import os +import shutil +import pytest +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_plate_instrument.py") + outdir = 'out.fccNi_PhononIncohEl_plate' + 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_plate_instrument.py") + from mcvine.acc.test.compare_acc_nonacc import compare_acc_nonacc + compare_acc_nonacc( + "fccNi_PhononIncohEl_plate", + ["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 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(1e3), 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__]) From 61f26675fcc80a94e3dca0ecb7a961de9fa11ee5 Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Wed, 23 Aug 2023 17:39:28 -0400 Subject: [PATCH 4/6] Fix angstrom/barn unit definitions --- acc/kernels/__init__.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/acc/kernels/__init__.py b/acc/kernels/__init__.py index 78f3464..38c18ef 100644 --- a/acc/kernels/__init__.py +++ b/acc/kernels/__init__.py @@ -153,13 +153,13 @@ def onPhonon_IncoherentElastic_Kernel(self, kernel): print(kernel.scatterer_origin.phase.unitcell) # additional kernel parameters - AA= units.angstrom + AA= units.length.angstrom dw_core = kernel.dw_core / AA**2 # additional kernel parameters - scattering_xs = kernel.scattering_xs/units.barn \ + scattering_xs = kernel.scattering_xs/units.area.barn \ if kernel.scattering_xs else 0. - absorption_xs = kernel.absorption_xs/units.barn \ + absorption_xs = kernel.absorption_xs/units.area.barn \ if kernel.absorption_xs else 0. from .Phonon_IncoherentElastic import scatter From 1e13e0120eb810e09dd84637d9ed90249c8ab6e6 Mon Sep 17 00:00:00 2001 From: Cole Kendrick Date: Thu, 7 Sep 2023 09:25:20 -0400 Subject: [PATCH 5/6] Switch test sample from plate to sphere --- ...e.py => HSS_fccNi_PhononIncohEl_sphere.py} | 0 .../fccNi_PhononIncohEl_plate_instrument.py | 20 ---------- .../fccNi_PhononIncohEl_sphere_instrument.py | 38 +++++++++++++++++++ ...atterer.xml => fccNi-sphere-scatterer.xml} | 0 .../incoh-el/sampleassembly.xml | 8 ++-- .../test_fccNi_Phonon_IncoherentElastic.py | 12 +++--- 6 files changed, 48 insertions(+), 30 deletions(-) rename tests/components/samples/{HSS_fccNi_PhononIncohEl_plate.py => HSS_fccNi_PhononIncohEl_sphere.py} (100%) delete mode 100644 tests/components/samples/fccNi_PhononIncohEl_plate_instrument.py create mode 100644 tests/components/samples/fccNi_PhononIncohEl_sphere_instrument.py rename tests/components/samples/sampleassemblies/incoh-el/{fccNi-plate-scatterer.xml => fccNi-sphere-scatterer.xml} (100%) diff --git a/tests/components/samples/HSS_fccNi_PhononIncohEl_plate.py b/tests/components/samples/HSS_fccNi_PhononIncohEl_sphere.py similarity index 100% rename from tests/components/samples/HSS_fccNi_PhononIncohEl_plate.py rename to tests/components/samples/HSS_fccNi_PhononIncohEl_sphere.py diff --git a/tests/components/samples/fccNi_PhononIncohEl_plate_instrument.py b/tests/components/samples/fccNi_PhononIncohEl_plate_instrument.py deleted file mode 100644 index 39580fa..0000000 --- a/tests/components/samples/fccNi_PhononIncohEl_plate_instrument.py +++ /dev/null @@ -1,20 +0,0 @@ -#!/usr/bin/env python - -import os -thisdir = os.path.dirname(__file__) - -from test_powderdiffraction_sample_instrument_factory import construct, Builder - -def instrument(is_acc=True): - if is_acc: - from HSS_fccNi_PhononIncohEl_plate 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=['PSD_4PI'], - builder=Builder()) 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/fccNi-plate-scatterer.xml b/tests/components/samples/sampleassemblies/incoh-el/fccNi-sphere-scatterer.xml similarity index 100% rename from tests/components/samples/sampleassemblies/incoh-el/fccNi-plate-scatterer.xml rename to tests/components/samples/sampleassemblies/incoh-el/fccNi-sphere-scatterer.xml diff --git a/tests/components/samples/sampleassemblies/incoh-el/sampleassembly.xml b/tests/components/samples/sampleassemblies/incoh-el/sampleassembly.xml index 2caa5c6..a95a449 100644 --- a/tests/components/samples/sampleassemblies/incoh-el/sampleassembly.xml +++ b/tests/components/samples/sampleassemblies/incoh-el/sampleassembly.xml @@ -2,11 +2,11 @@ - + - + - + Ni @@ -15,7 +15,7 @@ - +