Giter VIP home page Giter VIP logo

xsuite / xsuite Goto Github PK

View Code? Open in Web Editor NEW
24.0 5.0 22.0 27.29 MB

Suite of python packages for multiparticle simulations of particle accelerators.

Home Page: https://xsuite.readthedocs.io

License: Apache License 2.0

Python 79.75% Shell 8.28% Dockerfile 11.31% TeX 0.66%
physics gpu gpu-computing particle-accelerator physics-simulation particle-in-cell particle-tracking beam-optics beam-dynamics

xsuite's Issues

Features needed in monitor

Needed soon:

  • Possibility of passing a Monitor object to tracker to allow stopping and restarting the simulation in python updating the same monitor (needs additional tests)
  • Possibility of monitoring a subset of the particles (e.g. test particles for footprint)

Reorganise xtrack.Line

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.

Unequal tracking results on different contexts

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.

Change of reference momentum

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.

Aperture interpolation for collimation

Status:

  • Introduce Polygon
  • Introduce RectEllipse
  • Introduce Octagon (as Polygon)
  • general aperture interpolation based on polygons and convex hull implemented
  • Make quantitative test using a conical aperture
  • need to handle special case of repeated apertures
  • Need to skip interpolation if drift is very short
  • For now we support only CPU, needs to be asserted

Tests to be introduced in the suite

  • characterization of tilted/shifted apertures upstream and downstream
  • Cone example
  • Repeated aperture example (only in examples, not in tests for now)
  • Test for backtracking

To be moved to new issue:

  • Real test of aperture refinement against sixtrack
  • Need to handle first aperture marker in the line
  • Test for backtracking through crab cavities
  • Real loss map using particles from file or geant interface
  • Is the D1 and D2 aperture handled correctly in the MAD-X model?
  • Clarify meaning of apertures at the same s (what is the interpolation supposed to do?)

Error catching to be improved

Here a couple of situation where the error message should be more explicit:

Size of dynamic field not specified

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

Inconsistent context between kernel and arguments

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'>

Losses in collective elements

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

Update of reference energy now changes zeta. Is that correct?

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

}

Geant interface configuration (BDSIM)

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 ../../

Geant4 installation

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 ../../

BDSIM install

cd build/bdsim
cmake ../../sources/bdsim

Inheritance of elements

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.

To do for LossLocation refinement

To be moved to new issue:

  • Realistic test of aperture refinement against sixtrack
  • Need to handle first aperture marker in the line
  • Generate "partial backtracker" when not all elements are backtrackable (e.g. bb6d)
  • Test for backtracking through crab cavities
  • Real loss map using particles from file or geant interface
  • Is the D1 and D2 aperture handled correctly in the MAD-X model?
  • Clarify meaning of apertures at the same s (what is the interpolation supposed to do?)

`MISALIGNED_SUB_BUFFER_OFFSET` on pocl (intel)

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?

Handling of thick external elements in Xline

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.

Some material on readthedocs documentation

Memory problem in turn-by-turn monitor when using tracker multiple times

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

pmass value hardcoded in particles.py

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

Issue with linear tranformation in (zeta, delta)

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.

update_zeta_example-1

Add instructions for generating `test_data`

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.

How to handle different Particles classes

  • for the track method of the elements, which are compiled "just in time", the Particles class can be inferred from the particles object passed to the track method. We need to keep a dictionary of kernels (one per encountered Particles class) instead of a single kernel.
  • in the tracker the compilation of the track kernel could be delayed to the first track as for the elements and implement the same logic
  • The monitor needs to be rewritten in such a way that the storage is not a particle object but a generic array. In this way the monitor does not depend anymore on the particle class (accesses only the local particle) and can be handled by all the other elements

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    ๐Ÿ–– Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. ๐Ÿ“Š๐Ÿ“ˆ๐ŸŽ‰

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google โค๏ธ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.