Suite of python packages for multiparticle simulations for particle accelerators.
Documentation is on: https://xsuite.readthedocs.io/en/latest/.
Suite of python packages for multiparticle simulations of particle accelerators.
Home Page: https://xsuite.readthedocs.io
License: Apache License 2.0
Suite of python packages for multiparticle simulations for particle accelerators.
Documentation is on: https://xsuite.readthedocs.io/en/latest/.
Status:
Tests to be introduced in the suite
To be moved to new issue:
Some material
WSL https://docs.microsoft.com/en-us/windows/wsl/about
On gpu and docker
https://www.docker.com/blog/wsl-2-gpu-support-is-here/
On GPU and singularity
https://sylabs.io/guides/3.5/user-guide/gpu.html
Currently only some elements (drift) are taken into account when computing the total length and the S position. A thick BeamInteraction element, for example, will thus result in an incorrect machine length and S positions for donwstream elements.
Needs correct declarations in MANIFEST file
Needed soon:
Currently, xline.Line uses two main lists: elements and element_names. Lattice modification requires synchronous changes to both lists, which can sometimes be tedious and error-prone. Consider using tuples or class instances to tie together the elements and their names for easier lattice processing.
Good sphinx tutorial:
https://sphinx-tutorial.readthedocs.io/start/
Readthedocs tutorial
https://sphinx-rtd-tutorial.readthedocs.io/en/latest/index.html
Instructions to get a local installation (also on Mac):
https://docs.readthedocs.io/en/stable/development/install.html
On the table of contents
https://sphinx-rtd-theme.readthedocs.io/en/stable/configuring.html#how-the-table-of-contents-displays
A nice documentation:
https://docs.cupy.dev/en/stable/reference/generated/cupy.array.html#cupy.array
Using commit 16bea69, the test_collective_tracker.py
test fails on my machine with the following output:
================================================================================================================================== FAILURES ===================================================================================================================================
___________________________________________________________________________________________________________________________ test_collective_tracker ___________________________________________________________________________________________________________________________
def test_collective_tracker():
for CTX in xo.ContextCpu, xo.ContextPyopencl, xo.ContextCupy:
if CTX not in available:
continue
# ....
print('Track a few turns...')
n_turns = 10
> tracker.track(particles, num_turns=n_turns,
turn_by_turn_monitor=True)
tests/test_collective_tracker.py:72:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
xtrack/tracker.py:406: in _track_with_collective
pp.track(particles)
../xfields/xfields/beam_elements/spacecharge.py:231: in track
super().track(particles)
xtrack/base_element.py:119: in track
self._track_kernel(el=self._xobject, particles=particles)
../xobjects/xobjects/context_pyopencl.py:462: in __call__
arg_list.append(self.to_function_arg(arg, vv))
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
self = <xobjects.context_pyopencl.KernelPyopencl object at 0x7fee0d22ce50>, arg = <xobjects.context.Arg object at 0x7fee033dc4c0>, value = SpaceChargeBiGaussianData(...)
def to_function_arg(self, arg, value):
if arg.pointer:
if hasattr(arg.atype, "_dtype"): # it is numerical scalar
if isinstance(value, cl.Buffer):
return value
elif hasattr(value, "dtype"): # nparray
assert isinstance(value, cla.Array)
return value.base_data[value.offset :]
elif hasattr(value, "_shape"): # xobject array
raise NotImplementedError
else:
raise NotImplementedError
else:
raise ValueError(
f"Invalid value {value} for argument {arg.name} "
f"of kernel {self.description.pyname}"
)
else:
if hasattr(arg.atype, "_dtype"): # it is numerical scalar
return arg.atype(value) # try to return a numpy scalar
elif hasattr(arg.atype, "_size"): # it is a compound xobject
> return value._buffer.buffer[value._offset :]
E pyopencl._cl.RuntimeError: clCreateSubBuffer failed: MISALIGNED_SUB_BUFFER_OFFSET
../xobjects/xobjects/context_pyopencl.py:447: RuntimeError
The environment is PoCL 1.6 on *buntu 21.04 x86_64 on an Intel 11600K CPU. Does anybody else get the same issue?
replace_spaceharge_*
should be replace_spacecharge_*
Running pytest
locally on using commit 16bea69 fails two tests because the test_data/sps_w_spacecharge
directory does not contain the .json
files generated by 000_prepare_line.py
.
In order to successfully run the script to generate the missing files, installing lhcopt/lhcmask seems to be required (rather than pymask
from upstream). Maybe this is obvious to others, but would suggest to document this dependency along with the other ones in the install documentation.
Now scattered and in part duplicated across xline and pymask
Users would like:
The rest of the python code takes constants from scipy.constants
The C code has them hardcoded in constants.h
We shovel remote the hardcoded values in particles.py to have the same strategy across the python code and then review this whole constants business
To be moved to new issue:
In this way update on track can change the intensity
Currently inheritance of elements is not supported as the class name is hard coded in the c code. A workaround might be found on the python side, but the use case is not clear at the moment.
Energy update should use psigma instead of calculating everything from delta:
https://github.com/alatina/xtrack/blob/7781789bf246870048f154f0ac76801879290320/xtrack/particles/particles.py#L438
Have a look also at the definitions here:
https://github.com/xsuite/xline/blob/7e001261a273780fbb6c7a8a55069e1caadb707e/xline/particles.py#L20
Update of zeta is present in add_to_energy but not in update delta:
Most likely the same problem is present in xline (to be checked).
For sure it is there in sixtracklib:
See SixTrack/sixtracklib#139
This script triggers the issue:
import numpy as np
import xobjects as xo
import xline as xl
import xtrack as xt
## Generate a simple sequence
sequence = xl.Line(
elements=[xl.Drift(length=2.),
xl.Multipole(knl=[0, 1.], ksl=[0,0]),
xl.Drift(length=1.),
xl.Multipole(knl=[0, -1.], ksl=[0,0])],
element_names=['drift_0', 'quad_0', 'drift_1', 'quad_1'])
## Choose a context
context = xo.ContextCpu() # For CPU
# context = xo.ContextCupy() # For CUDA GPUs
# context = xo.ContextPyopencl() # For OpenCL GPUs
## Transfer lattice on context and compile tracking code
tracker = xt.Tracker(_context=context, sequence=sequence)
## Build particle object on context
n_part = 200
particles = xt.Particles(_context=context,
p0c=6500e9,
x=np.random.uniform(-1e-3, 1e-3, n_part),
px=np.random.uniform(-1e-5, 1e-5, n_part),
y=np.random.uniform(-2e-3, 2e-3, n_part),
py=np.random.uniform(-3e-5, 3e-5, n_part),
zeta=np.random.uniform(-1e-2, 1e-2, n_part),
delta=np.random.uniform(-1e-4, 1e-4, n_part),
)
## Track (saving turn-by-turn data)
n_turns = 100
tracker.track(particles, num_turns=n_turns,
turn_by_turn_monitor=True)
## Track second time (saving turn-by-turn data)
tracker.track(particles, num_turns=n_turns,
turn_by_turn_monitor=True)
## Try to access turn-by-turn data
tracker.record_last_track.x
On CPU I get:
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-3-1fe53bc9bc1a> in <module>
43
44 ## Try to access turn-by-turn data
---> 45 tracker.record_last_track.x
~/Desktop/20210303_xfields_gpudev/xtrack/xtrack/monitors.py in __get__(self, container, ContainerType)
59
60 def __get__(self, container, ContainerType=None):
---> 61 vv = getattr(container.data, self.name)
62 n_cols = container.stop_at_turn - container.start_at_turn
63 n_rows = container.n_records // n_cols
~/Desktop/20210303_xfields_gpudev/xtrack/xtrack/dress.py in __get__(self, container, ContainerType)
14 def __get__(self, container, ContainerType=None):
15 if self.isnplikearray:
---> 16 return getattr(container._xobject, self.name).to_nplike()
17 elif hasattr(container, '_dressed_'+self.name):
18 return getattr(container, '_dressed_'+self.name)
~/Desktop/20210303_xfields_gpudev/xobjects/xobjects/array.py in to_nplike(self)
573 cshape = [shape[ii] for ii in self._order]
574 if hasattr(self._itemtype, "_dtype"):
--> 575 arr = self._buffer.to_nplike(
576 self._offset + self._data_offset, self._itemtype._dtype, cshape
577 ).transpose(self._order)
~/Desktop/20210303_xfields_gpudev/xobjects/xobjects/context_cpu.py in to_nplike(self, offset, dtype, shape)
445 # dtype=np.dtype(dtype)
446 # return self.buffer[offset:].view(dtype)[:count].reshape(*shape)
--> 447 return np.frombuffer(
448 self.buffer, dtype=dtype, count=count, offset=offset
449 ).reshape(*shape)
ValueError: buffer is smaller than requested size
Hi,
I was giving a go at the simple FODO example from the documentation, and tried the different contexts available to me.
However, it seems that doing the same tracking on two different contexts yields different results.
I am on a macbook pro machine (Intel CPU + Radeon Pro Vega 20 GPU), and tried the CPU
and OpenCL
contexts.
Here is a script which reproduces the issue for me:
from typing import Union
import numpy as np
import xline as xl
import xobjects as xo
import xtrack as xt
SEQUENCE = xl.Line(
elements=[
xl.Drift(length=2.0),
xl.Multipole(knl=[0, 1.0], ksl=[0, 0]),
xl.Drift(length=1.0),
xl.Multipole(knl=[0, -1.0], ksl=[0, 0]),
],
element_names=["drift_0", "quad_0", "drift_1", "quad_1"],
)
def track_on_context(
_context: Union[xo.ContextCpu, xo.ContextCupy, xo.ContextPyopencl],
n_part: int = 200,
n_turns: int = 10_000,
) -> np.ndarray:
"""
Convenience function to track as the documentation's simple example on a given context, returning the
resulting horizontal coordinates.
Args:
_context (Union[xo.ContextCpu, xo.ContextCupy, xo.ContextPyopencl]): `xobjects` context to build for.
n_part (int): number of particles to track, defaults to 200.
n_turns (int): number of turns to track for, defaults to 10000.
Returns:
Resulting X coordinates from tracking, as a `numpy.ndarray`.
"""
# Transfer lattice on context and compile tracking code
tracker = xt.Tracker(_context=_context, sequence=SEQUENCE)
# Build particle object on context
particles = xt.Particles(
_context=_context,
p0c=6500e9,
x=np.random.uniform(-1e-3, 1e-3, n_part),
px=np.random.uniform(-1e-5, 1e-5, n_part),
y=np.random.uniform(-2e-3, 2e-3, n_part),
py=np.random.uniform(-3e-5, 3e-5, n_part),
zeta=np.random.uniform(-1e-2, 1e-2, n_part),
delta=np.random.uniform(-1e-4, 1e-4, n_part),
)
# Track and return turn-by-turn
tracker.track(particles, num_turns=n_turns, turn_by_turn_monitor=True)
return tracker.record_last_track.x
# Track on CPU Context
context_cpu = xo.ContextCpu()
x_cpu = track_on_context(_context=context_cpu)
# Track on GPU Context
# macOS machine, using the OpenGL context as Metal is not supported
context_gpu = xo.ContextPyopencl()
x_gpu = track_on_context(_context=context_gpu)
# Comparing Tracking Data
print(np.array_equal(x_cpu, x_gpu))
The output in my case is, disregarding the compilation information, False
.
Having a look at the arrays themselves, it is apparent they're very different.
I am available if you needs help debugging this on my machine.
Not ideal...
Compiled code does not show a cpu dependence. Numpy instead does due to the use of simd intrinsics and usage of optimised blas and lapack libraries:
https://numpy.org/devdocs/reference/simd/simd-optimizations.html
https://numpy.org/doc/stable/user/building.html
To achieve numerical portability (mostly of pymask) we need to compile numpy and scipy with these features disabled and with the reference netlib implementation for blas and lapack.
A change of reference momentum is needed to study ramping machines. It will be needed for the muon collider study (featuring re-circulating LINACs and rapid-cycling synchrotrons). To run on GPUs, this feature requires that the reference momentum (and associated parameters) become attributes of the local particles. For the moment a workaround is used by adjusting the reference momentum on the python side.
Losses are incompatible with pyht interface as x, y, ... are not cut at _num_active_particles
How do we implement it?
(It should not degrade the performance when not uses)
Could be mentioned in installation pages with instructions on how to get them
http://www.pp.rhul.ac.uk/bdsim/manual/installation.html
sudo apt install cmake-curses-gui
sudo apt install libxmu-dev
sudo apt install libmotif-dev
sudo apt-get install libglfw3-dev libgles2-mesa-dev
conda install matplotlib
conda install -c anaconda xerces-c
conda install -c conda-forge flex
conda install -c conda-forge bison
conda install -c conda-forge clhep
conda install -c conda-forge root # be patient!
conda install -c anaconda mesa-libgl-cos6-x86_64
conda install mesa-libgl-devel-cos6-x86_64 mesa-dri-drivers-cos6-x86_64 libselinux-cos6-x86_64 libxdamage-cos6-x86_64 libxxf86vm-cos6-x86_64 clhep=2.4.4.0 expat freetype libglu qt=5.12 xerces-c xorg-libx11 xorg-libxfixes xorg-libxmu zlib xorg-libxfixes libxext-cos6-x86_64
mkdir bdsim_env
cd bdsim_env
mkdir sources
cd sources
git clone --depth 1 --branch v1.6.0 https://bitbucket.org/jairhul/bdsim.git
git clone --depth 1 --branch v10.4.3 https://github.com/Geant4/geant4.git
cd ../
mkdir build
mkdir install
mkdir build/geant4
mkdir build/bdsim
mkdir install/geant4
mkdir install/bdsim
make -j 4 install
cd ../../
cd build/geant4/
cmake ../../sources/geant4/
ccmake .
Set
CMAKE_BUILD_TYPE Release
CMAKE_INSTALL_PREFIX /home/giadarol/Desktop/20210303_xfields_gpudev/bdsim_env/install/geant4
GEANT4_BUILD_MULTITHREADED OFF
GEANT4_INSTALL_DATA ON
GEANT4_INSTALL_DATADIR
GEANT4_USE_G3TOG4 OFF
GEANT4_USE_GDML OFF
GEANT4_USE_INVENTOR OFF
GEANT4_USE_OPENGL_X11 ON
GEANT4_USE_QT ON
GEANT4_USE_RAYTRACER_X11 ON
GEANT4_USE_SYSTEM_CLHEP ON
GEANT4_USE_SYSTEM_EXPAT ON
GEANT4_USE_SYSTEM_ZLIB OFF
GEANT4_USE_XM ON
Hit c
to configure and then g
to generate
need also (when asked):
CLHEP_DIR /home/giadarol/Desktop/20210303_xfields_gpudev/bdsim_env/install/clhep/lib/CLHEP-2.4.5.1
make -j 24 install
cd ../../
cd build/bdsim
cmake ../../sources/bdsim
Here a couple of situation where the error message should be more explicit:
import xobjects as xo
class MyStruct(xoStruct):
a = xo.Float64[:]
s = MyStruct()
gives:
---------------------------------------------------------------------------
IndexError Traceback (most recent call last)
<ipython-input-4-a5a954f0bd97> in <module>
----> 1 s = MyStruct()
~/Desktop/20210303_xfields_gpudev/xobjects/xobjects/struct.py in __init__(self, _context, _buffer, _offset, *args, **kwags)
326 args, kwargs = cls._pre_init(*args, **kwargs)
327 # compute resources
--> 328 info = cls._inspect_args(*args, **kwargs)
329 self._size = info.size
330 # acquire buffer
[...]
~/Desktop/20210303_xfields_gpudev/xobjects/xobjects/array.py in _inspect_args(cls, *args)
310 value = args[0]
311 else: # complete dimensions
--> 312 if not is_integer(args[0]): # init with array
313 value = args[0]
314 shape = get_shape_from_array(value)
IndexError: tuple index out of range
Using xtrack for simplicity:
import xobjects as xo
import xtrack as xt
ctx1 = xo.ContextCpu()
ctx2 = xo.ContextCupy()
part = xt.Particles(_context=ctx1, p0c=6.5e12, x=[1,2,3])
drift = xt.Drift(_context=ctx2, length=1.)
drift.track(part)
gives:
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-5-d4e36c877d76> in <module>
7 drift = xt.Drift(_context=ctx2, length=1.)
8
----> 9 drift.track(part)
~/Desktop/20210705_xsuite_dev/xtrack/xtrack/dress_element.py in track(self, particles)
63
64 self._track_kernel.description.n_threads = particles.num_particles
---> 65 self._track_kernel(el=self._xobject, particles=particles)
66
67 DressedElement.compile_track_kernel = compile_track_kernel
~/Desktop/20210705_xsuite_dev/xobjects/xobjects/context_cupy.py in __call__(self, **kwargs)
376
377 grid_size = int(np.ceil(n_threads / self.block_size))
--> 378 self.function((grid_size,), (self.block_size,), arg_list)
379
380
cupy/_core/raw.pyx in cupy._core.raw.RawKernel.__call__()
cupy/cuda/function.pyx in cupy.cuda.function.Function.__call__()
cupy/cuda/function.pyx in cupy.cuda.function._launch()
cupy/cuda/function.pyx in cupy.cuda.function._pointer()
TypeError: Unsupported type <class 'numpy.ndarray'>
At the moment it is assumed that the collective element takes care of incrementing at_element
To be checked with @rdemaria
Here the parts of the code
/*gpufun*/
void LocalParticle_update_p0c(LocalParticle* part, double new_p0c_value){
double const mass0 = LocalParticle_get_mass0(part);
double const old_p0c = LocalParticle_get_p0c(part);
double const old_delta = LocalParticle_get_delta(part);
double const ppc = old_p0c * old_delta + old_p0c;
double const new_delta = (ppc - new_p0c_value)/new_p0c_value;
double const new_energy0 = sqrt(new_p0c_value*new_p0c_value + mass0 * mass0);
double const new_beta0 = new_p0c_value / new_energy0;
double const new_gamma0 = new_energy0 / mass0;
LocalParticle_set_p0c(part, new_p0c_value);
LocalParticle_set_gamma0(part, new_gamma0);
LocalParticle_set_beta0(part, new_beta0);
LocalParticle_update_delta(part, new_delta);
// TODO: This changes zeta. Is this correct?
}
where:
/*gpufun*/
void LocalParticle_update_delta(LocalParticle* part, double new_delta_value){
double const beta0 = LocalParticle_get_beta0(part);
double const delta_beta0 = new_delta_value * beta0;
double const ptau_beta0 = sqrt( delta_beta0 * delta_beta0 +
2. * delta_beta0 * beta0 + 1. ) - 1.;
double const one_plus_delta = 1. + new_delta_value;
double const rvv = ( one_plus_delta ) / ( 1. + ptau_beta0 );
double const rpp = 1. / one_plus_delta;
double const psigma = ptau_beta0 / ( beta0 * beta0 );
LocalParticle_set_delta(part, new_delta_value);
LocalParticle_scale_zeta(part,
rvv / LocalParticle_get_rvv(part));
LocalParticle_set_rvv(part, rvv );
LocalParticle_set_rpp(part, rpp );
LocalParticle_set_psigma(part, psigma );
}
From @xbuffat:
On my version of xtrack (which is up to date with the main one, but with an additional element that tracks with sigma/psigma instead of zeta/delta, similar to the one in my pull request) I put an example that compares the two options:
https://github.com/xbuffat/xtrack/blob/main/examples/update_zeta_example.py
The problem with the zeta/delta tracking is that the function update_zeta modifies zeta according to beta which varies during linear synchrotron motion. This appears as an additional contribution to the linear motion that adds up over the turns (red dots in the figure attached). On the other hand, with a linear matrix in sigma/psigma, the motion remains stable (blue dots). If one looks carefully, the phase space is not exactly a circle in zeta/delta, it is slightly distorted due to the beta dependence of zeta, whereas in sigma/psigma we do get the perfect circle.
At the moment losses are not correctly handled in the presence of collective elements. In particular check_is_not_lost
is not called before calling the .track
python method of the collective element
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.