hightower8083 / synchrad

Synchrotron Radiation calculator via openCL

Synchrotron Radiation calculator via openCL

Overview

Tool to compute energy density of synchrotron radiation in the spectral volume using the Fourier transformed Lienard-Wiechert potentials (see eq. (14.65) of J.D.Jackson's Classical electrodynamics).

This software processes the 3D (x, y, z, px, py, pz) trajectories of the charged point-like particles with weights, and maps the emittied energy to the 3D spectral domain (omega, theta, phi). The package also contains a number of convenience methods for pre- and post-proccessing.

Language and hardware

SynchRad is written in openCL and is interfaced with Python via PyOpenCL.

Code also uses Mako template manager for tweaking with the data types and native functions. It was tested on GPU and CPU devices using NVIDIA, AMD and Apple platforms, and it demonstrates a reasonable performance on GPUs, while on CPU the simillar OpenMP implementation is signifincantly faster.

Installation

Given that openCL and PyOpenCL are installed (e.g. via conda or pip) and configured on your machine, you can install SynchRad in a standard way, i.e. by cloning the source and running the setup.py script:

git clone https://github.com/hightower8083/synchrad.git
cd synchrad/
python setup.py install

To be able to use the software with multiple GPU or CPU devices via MPI one should also install mpi4py. To output result in VTK format via exportToVTK method, the tvtk.api should be installed.

Usage

A minimal example of SynchRad usage can be found in example/ folder of this repository.

Another common example would be to calculate radiation produced by the particles from, for example, a PIC simulation. In case if PIC software supports output in OpenPMD standard, this can be done with help of openPMD-viewer, using the conversion function:

from openpmd_viewer import OpenPMDTimeSeries, ParticleTracker
from synchrad.utils import tracksFromOPMD

ts = LpaDiagnostics('./run/diags_track/hdf5/', check_all_files=False)
ref_iteration = ts.iterations[-1]
pt = ParticleTracker( ts, iteration=ref_iteration, 
                      preserve_particle_index=True)

tracksFromOPMD( ts, pt, ref_iteration=ref_iteration, fname='tracks.h5')

where one can see doc-lines of tracksFromOPMD for more optional arguments.

The prepared track can be processed and saved to file using following commands:

import numpy as np
from synchrad.calc import SynchRad

calc_input = {'grid':[ (1e3/1.24e-6, 75e3/1.24e-6),
                       (0, 0.04),
                       (0.,2*np.pi),
                       (256, 32, 32) ],
              'dtype':'double',
             }

calc = SynchRad(calc_input)
calc.calculate_spectrum(file_tracks='./tracks.h5',file_spectrum='./spectrum.h5')

where radiation within 40 urad angle is calculated for the energies range [1 keV, 75 keV].

For more details see the example notebooks in tutorials/

Author and Contributions

This software is developed by Igor A Andriyash ([email protected]), and it is on the early stage of development.

Everyone is welcome to contribute either by testing and benchmarking, or by introducing further optimizations and adding utility methods.

synchrad's Issues

GitHub flow

Is there a particular reason for using a more complex dev/master model instead of the simpler GitHub flow based on a single branch?

`numba` is also a dependency

It seems numba should be added to the requirements file.
Also, openpmd_viewer 1.X uses
from openpmd_viewer import .. instead of
from opmd_viewer import .., as stated in the README

Problem trying to run synchrad with tutorials

Hello @hightower8083 and all,

I have been studying Synchrad recently,Thank you for the code!

There was a problem when I was running the Betatron_Example in tutorials and I am not sure why. I ran it on my laptop without a GPU.

    "grid": [
        (1e-3*omega_ch_crit/omega_1m, omega_ch_crit/omega_1m),  
        (0, 2*K0/gamma_0),
        (0.0, 2 * np.pi),
        (512, 32, 32),
    ],  
}
calc = SynchRad(calc_input)
calc.calculate_spectrum( particleTracks=tracks,
                         #comp='cartesian_complex',
                         file_spectrum = './spect_incoh.h5',
                         timeStep=c*dt )
calc.calculate_spectrum( particleTracks=tracks,
                         comp='cartesian_complex',
                         file_spectrum = './spect_coh.h5',
                         timeStep=c*dt )

when I choose platform Intel(R) OpenCL Graphics the error message:

Choose platform:
[0] <pyopencl.Platform 'Intel(R) OpenCL Graphics' at 0x1df89a9c6d0>
Choice [0]:Intel(R) OpenCL Graphics
Set the environment variable PYOPENCL_CTX='Intel(R) OpenCL Graphics' to avoid being asked again.
Running on 1 devices
  ALL | GPU device: Intel(R) Iris(R) Xe Graphics
Platform: Intel(R) Corporation
Compiler: OpenCL C 1.2 
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[20], line 10
      1 calc_input = {
      2     "grid": [
      3         (1e-3*omega_ch_crit/omega_1m, omega_ch_crit/omega_1m),  
   (...)
      7     ],  
      8 }
---> 10 calc = SynchRad(calc_input)
     11 calc.calculate_spectrum( particleTracks=tracks,
     12                          #comp='cartesian_complex',
     13                          file_spectrum = './spect_incoh.h5',
     14                          timeStep=c*dt )
     16 calc.calculate_spectrum( particleTracks=tracks,
     17                          comp='cartesian_complex',
     18                          file_spectrum = './spect_coh.h5',
     19                          timeStep=c*dt )

File D:\anaconda\Lib\site-packages\synchrad\calc.py:95, in SynchRad.__init__(self, Args, file_spectrum)
     93     self._init_comm()
     94     self._init_data()
---> 95     self._compile_kernels()
     96 else:
     97     self._read_args(file_spectrum)

File D:\anaconda\Lib\site-packages\synchrad\calc.py:596, in SynchRad._compile_kernels(self)
    593     fname += "kernel_nearfield.cl"
    595 src = Template( filename=fname ).render(**agrs)
--> 596 self._mapper = cl.Program(self.ctx, src).build()

File D:\anaconda\Lib\site-packages\pyopencl\__init__.py:535, in Program.build(self, options, devices, cache_dir)
    531 else:
    532     # cached
    534     from pyopencl.cache import create_built_program_from_source_cached
--> 535     self._prg, was_cached = self._build_and_catch_errors(
    536             lambda: create_built_program_from_source_cached(
    537                 self._context, self._source, options_bytes, devices,
    538                 cache_dir=cache_dir, include_path=include_path),
    539             options_bytes=options_bytes, source=self._source)
    541     if was_cached:
    542         build_descr = "cache retrieval"

File D:\anaconda\Lib\site-packages\pyopencl\__init__.py:583, in Program._build_and_catch_errors(self, build_func, options_bytes, source)
    575     err = _cl.RuntimeError(
    576             _cl._ErrorRecord(
    577                 msg=msg,
    578                 code=code,
    579                 routine=routine))
    581 # Python 3.2 outputs the whole list of currently active exceptions
    582 # This serves to remove one (redundant) level from that nesting.
--> 583 raise err

RuntimeError: clBuildProgram failed: BUILD_PROGRAM_FAILURE - clBuildProgram failed: BUILD_PROGRAM_FAILURE - clBuildProgram failed: BUILD_PROGRAM_FAILURE

Build on <pyopencl.Device 'Intel(R) Iris(R) Xe Graphics' on 'Intel(R) OpenCL Graphics' at 0x1df8c15ba40>:

5:3:26: warning: unsupported OpenCL extension 'cl_khr_fp64' - ignoring
#pragma OPENCL EXTENSION cl_khr_fp64 : enable
                         ^
5:7:12: error: use of type 'double' requires cl_khr_fp64 extension to be enabled
  __global double *spectrum,
           ^
5:8:12: error: use of type 'double' requires cl_khr_fp64 extension to be enabled
  __global double *x,
           ^
5:9:12: error: use of type 'double' requires cl_khr_fp64 extension to be enabled
  __global double *y,
           ^
5:10:12: error: use of type 'double' requires cl_khr_fp64 extension to be enabled
  __global double *z,
           ^
5:11:12: error: use of type 'double' requires cl_khr_fp64 extension to be enabled
  __global double *ux,
           ^
5:12:12: error: use of type 'double' requires cl_khr_fp64 extension to be enabled
  __global double *uy,
           ^
5:13:12: error: use of type 'double' requires cl_khr_fp64 extension to be enabled
  __global double *uz,
           ^
5:14:12: error: use of type 'double' requires cl_khr_fp64 extension to be enabled
           double wp,
           ^
5:18:12: error: use of type 'double' requires cl_khr_fp64 extension to be enabled
  __global double *omega,
           ^
5:19:12: error: use of type 'double' requires cl_khr_fp64 extension to be enabled
  __global double *sinTheta,
           ^
5:20:12: error: use of type 'double' requires cl_khr_fp64 extension to be enabled
  __global double *cosTheta,
           ^
5:21:12: error: use of type 'double' requires cl_khr_fp64 extension to be enabled
  __global double *sinPhi,
           ^
5:22:12: error: use of type 'double' requires cl_khr_fp64 extension to be enabled
  __global double *cosPhi,
           ^
5:26:12: error: use of type 'double' requires cl_khr_fp64 extension to be enabled
           double dt,
           ^
5:39:5: error: use of type 'double' requires cl_khr_fp64 extension to be enabled
    double omegaLocal = omega[iOmega];
    ^
5:40:5: error: use of undeclared identifier 'double3'; did you mean 'double'?
    double3 nVec = (double3) { sinTheta[iTheta]*cosPhi[iPhi],
    ^~~~~~~
    double
5:40:5: error: use of type 'double' requires cl_khr_fp64 extension to be enabled
5:40:21: error: use of undeclared identifier 'double3'; did you mean 'double'?
    double3 nVec = (double3) { sinTheta[iTheta]*cosPhi[iPhi],
                    ^~~~~~~
                    double
5:40:29: error: expected ';' at end of declaration
    double3 nVec = (double3) { sinTheta[iTheta]*cosPhi[iPhi],
                            ^
                            ;
5:44:5: error: use of undeclared identifier 'double3'; did you mean 'double'?
    double3 xLocal, uLocal, uNextLocal, aLocal, amplitude;
    ^~~~~~~
    double
5:44:5: error: use of type 'double' requires cl_khr_fp64 extension to be enabled
5:45:5: error: use of type 'double' requires cl_khr_fp64 extension to be enabled
    double time, phase, dPhase, sinPhase, cosPhase, c1, c2, gammaInv;
    ^
5:47:5: error: use of type 'double' requires cl_khr_fp64 extension to be enabled
    double dtInv = (double)1. / dt;
    ^
5:47:28: warning: double precision constant requires cl_khr_fp64, casting to single precision
    double dtInv = (double)1. / dt;
                           ^
5:47:21: error: use of type 'double' requires cl_khr_fp64 extension to be enabled
    double dtInv = (double)1. / dt;
                    ^
5:48:5: error: use of type 'double' requires cl_khr_fp64 extension to be enabled
    double wpdt2 =  wp * dt * dt;
    ^
5:49:5: error: use of type 'double' requires cl_khr_fp64 extension to be enabled
    double phasePrev = (double) 0.;
    ^
5:49:33: warning: double precision constant requires cl_khr_fp64, casting to single precision
    double phasePrev = (double) 0.;
                                ^
5:49:25: error: use of type 'double' requires cl_khr_fp64 extension to be enabled
    double phasePrev = (double) 0.;
                        ^
5:50:5: error: use of undeclared identifier 'double3'; did you mean 'double'?
    double3 spectrLocalRe = (double3) {0., 0., 0.};
    ^~~~~~~
    double
5:50:5: error: use of type 'double' requires cl_khr_fp64 extension to be enabled
5:50:30: error: use of undeclared identifier 'double3'; did you mean 'double'?
    double3 spectrLocalRe = (double3) {0., 0., 0.};
                             ^~~~~~~
                             double
5:50:38: error: expected ';' at end of declaration
    double3 spectrLocalRe = (double3) {0., 0., 0.};
                                     ^
                                     ;
5:51:5: error: use of undeclared identifier 'double3'; did you mean 'double'?
    double3 spectrLocalIm = (double3) {0., 0., 0.};
    ^~~~~~~
    double
5:51:5: error: use of type 'double' requires cl_khr_fp64 extension to be enabled
5:51:30: error: use of undeclared identifier 'double3'; did you mean 'double'?
    double3 spectrLocalIm = (double3) {0., 0., 0.};
                             ^~~~~~~
                             double
5:51:38: error: expected ';' at end of declaration
    double3 spectrLocalIm = (double3) {0., 0., 0.};
                                     ^
                                     ;
5:65:17: error: use of type 'double' requires cl_khr_fp64 extension to be enabled
        time = (double)it * dt;
                ^
5:66:19: error: use of undeclared identifier 'double3'; did you mean 'double'?
        xLocal = (double3) {x[it], y[it], z[it]};
                  ^~~~~~~
                  double

(options: -I D:\anaconda\Lib\site-packages\pyopencl\cl)

Am I getting something wrong?

Thank you for your help!

useful resources for testing/further development

I am compiling a list of useful papers for betatron emission calculation.

• Richard's Thesis
• How to Test and Verify Radiation Diagnostics Simulations
• Quantitatively Consistent Computation of Coherent and Incoherent Radiation
• Computing Angularly-Resolved Far-Field Emission Spectra

Increased clarity in input/output

I have been playing around with Synchrad for a little while now trying to apply it to a non-standard use-case, where there is no obvious "fundamental frequency" like in an undulator or betatron radiation. The tutorials are based around these schemes, and while it is possible for me to extract some methodology from these tutorials, it would be of much help to be able to see what outputs are available from the different functions.

Maybe it is because I am new to Python, but it is not at all obvious to me how to see what is output by the different functions in the code. In Spyder, I can check the main Synchrad module, but it doesn't tell me anything about what information is available after calc = SynchRad( calc_input ) has been run. Moreover, if I then try to look into calc in the variable browser, it says that "the variable is not picklable". I would therefore very much appreciate that something is added at least to the documentation stating what output is available, as well as what units it is given in.

Something which I am currently trying to do is to simply calculate the total photon number / emitted energy, but in for example the undulator tutorial there are some unclear inputs that are (at least to me) not obvious, such as the lambda0_um and phot_num variables. phot_num evidently changes the unit and/or scaling of the output to something I don't understand. Or rather: I guess that calc.get_energy( phot_num=True, normalize_to_weights=True ) should give me the emitted average number of photons per electron, but if I set phot_num=False, I get a very large value which seems to be in neither J nor eV.

I am also interested in getting the radiated power as a function of propagation distance/time (which I guess I can set with the variable nSnaps in the calculate_spectrum?) as well as the physical near-field source size (to be able to calculate the peak brilliance), but again it is not obvious to me how to handle the outputs.

Any clarification here would be much appreciated!
Best regards,
Jonas

no locator available when extracting particle tracks

$ git checkout dev
$ cd example/PIC
$ python fbpic_lwfa_script.py  # takes a while, but works and produces diags/ and diags_track/
$ python extract_tracks.py

Traceback (most recent call last):
  File "extract_tracks.py", line 1, in <module>
    from openpmd_viewer import ParticleTracker
  File "<frozen importlib._bootstrap>", line 991, in _find_and_load
  File "<frozen importlib._bootstrap>", line 975, in _find_and_load_unlocked
  File "<frozen importlib._bootstrap>", line 655, in _load_unlocked
  File "<frozen importlib._bootstrap>", line 618, in _load_backward_compatible
  File "<frozen zipimport>", line 259, in load_module
  File "/data/storage/berceanu/software/pkg/minicondapic/lib/python3.8/site-packages/openPMD_viewer-1.0.1-py3.8.egg/openpmd_viewer/__init__.py", line 9, in <module>
  File "<frozen importlib._bootstrap>", line 991, in _find_and_load
  File "<frozen importlib._bootstrap>", line 975, in _find_and_load_unlocked
  File "<frozen importlib._bootstrap>", line 655, in _load_unlocked
  File "<frozen importlib._bootstrap>", line 618, in _load_backward_compatible
  File "<frozen zipimport>", line 259, in load_module
  File "/data/storage/berceanu/software/pkg/minicondapic/lib/python3.8/site-packages/openPMD_viewer-1.0.1-py3.8.egg/openpmd_viewer/openpmd_timeseries/__init__.py", line 2, in <module>
  File "<frozen importlib._bootstrap>", line 991, in _find_and_load
  File "<frozen importlib._bootstrap>", line 975, in _find_and_load_unlocked
  File "<frozen importlib._bootstrap>", line 655, in _load_unlocked
  File "<frozen importlib._bootstrap>", line 618, in _load_backward_compatible
  File "<frozen zipimport>", line 259, in load_module
  File "/data/storage/berceanu/software/pkg/minicondapic/lib/python3.8/site-packages/openPMD_viewer-1.0.1-py3.8.egg/openpmd_viewer/openpmd_timeseries/main.py", line 13, in <module>
  File "<frozen importlib._bootstrap>", line 991, in _find_and_load
  File "<frozen importlib._bootstrap>", line 975, in _find_and_load_unlocked
  File "<frozen importlib._bootstrap>", line 655, in _load_unlocked
  File "<frozen importlib._bootstrap>", line 618, in _load_backward_compatible
  File "<frozen zipimport>", line 259, in load_module
  File "/data/storage/berceanu/software/pkg/minicondapic/lib/python3.8/site-packages/openPMD_viewer-1.0.1-py3.8.egg/openpmd_viewer/openpmd_timeseries/utilities.py", line 223, in <module>
  File "/data/storage/berceanu/software/pkg/minicondapic/lib/python3.8/site-packages/numba/core/decorators.py", line 211, in wrapper
    disp.enable_caching()
  File "/data/storage/berceanu/software/pkg/minicondapic/lib/python3.8/site-packages/numba/core/dispatcher.py", line 743, in enable_caching
    self._cache = FunctionCache(self.py_func)
  File "/data/storage/berceanu/software/pkg/minicondapic/lib/python3.8/site-packages/numba/core/caching.py", line 620, in __init__
    self._impl = self._impl_class(py_func)
  File "/data/storage/berceanu/software/pkg/minicondapic/lib/python3.8/site-packages/numba/core/caching.py", line 355, in __init__
    raise RuntimeError("cannot cache function %r: no locator available "
RuntimeError: cannot cache function 'histogram_cic_1d': no locator available for file '/data/storage/berceanu/software/pkg/minicondapic/lib/python3.8/site-packages/openPMD_viewer-1.0.1-py3.8.egg/openpmd_viewer/openpmd_timeseries/utilities.py'

I am using openPMD_viewer-1.0.1.

Single trace evaluation

Hi @hightower8083,

I recently started using synchrad for ICS calculations and enjoy it very much, especially the ease/elegance with which it can be combined with other python-based code.
Before, I used clara2 (https://github.com/ComputationalRadiationPhysics/clara2.git). A very useful tool in clara2, I found, is the option to obtain the classical radiation emitted by each macroparticle. This allows for great diagnostics, extracting the dependence of the contribution to the full spectrum on certain particle parameters. It would be of high interest to me, if this could be implemented in your code, as well.

calculation of energy radiated

Apart from computing the spectrum, is the code capable of calculating the net energy radiated?

Radiation calculation

Dear Developers,

I would like to use "synchrad" for calculating the radiation emitted by the oscillating electron beam in LWFA.

Could you please provide me some initial help on how to do it if I have position and momentum information of the electrons.

I am using FBPIC code for studying radiation generated by electron beams in LWFA.

Thanks a lot.
Naka

PIC example broken?

Running on 1 devices
  Starting without device: 
Platform: None
Compiler: None
Traceback (most recent call last):
  File "compute_spectrum.py", line 18, in <module>
    file_spectrum="spectrum.h5",)
  File "/home/berceanu/Development/synchrad/synchrad/calc.py", line 44, in calculate_spectrum
    self._init_raditaion(comp, nSnaps)
  File "/home/berceanu/Development/synchrad/synchrad/calc.py", line 293, in _init_raditaion
    arrcl.zeros(self.queue, radiation_shape, dtype=self.dtype)
AttributeError: 'SynchRad' object has no attribute 'queue'

Steps to reproduce:

$ git checkout dev
$ cd example/PIC
$ python fbpic_lwfa_script.py # produces diags/ and diags_track/
$ python extract_tracks.py # produces tracks.h5
$ python compute_spectrum.py # produces above error, instead of spectrum.h5

Problem trying to run synchrad across multiple nodes

Hello @hightower8083 and all,

I have been using Synchrad recently to calculate coherent radiation of a beam through an undulator.
Thank you for the code!

For my study, it has become clear that I need more macro particles to find convergence of the results, but the simulation already takes about 25 hours in a 4 x GPU (A100) node.
It'd be great to be able to run across multiple nodes to use more GPUs and save some time.
However, I failed on my first try and I am not sure why.

In the submission script, I simply increased the number of requested nodes and adjusted the number of mpi process to use.
This is an example with 2 nodes:

#!/bin/bash -x
#SBATCH --job-name=synchrad
#SBATCH --partition=mpa
#SBATCH --nodes=2
#SBATCH --constraint="GPUx4&A100"
#SBATCH --time=48:00:00
#SBATCH --output=stdout
#SBATCH --error=stderr
#SBATCH --mail-type=END

# Activate environment
source $HOME/synchrad_env.sh
export PYOPENCL_CTX=':'

mpirun -n 8 python undulator_beam.py

The error message is not really helpful to me:

===================================================================================
=   BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
=   PID 140773 RUNNING AT max-mpag009
=   EXIT CODE: 9
=   CLEANING UP REMAINING PROCESSES
=   YOU CAN IGNORE THE BELOW CLEANUP MESSAGES
===================================================================================
YOUR APPLICATION TERMINATED WITH THE EXIT STRING: Killed (signal 9)
This typically refers to a problem with your application.
Please see the FAQ page for debugging suggestions

Am I forgetting anything?

Thank you for your help!

coherent synchrotron radiation

Hi,
We're looking to check the code corectness for the coherent part of the synchrotron emission spectrum by comparing with the PIConGPU radiation plugin output for a simple LWFA setup.

While we get relatively good agreement for the incoherent part of the spectrum, the coherent part seems to be a global upward shift of the whole spectrum in the case of synchrad, while for PIConGPU, the coherent low-energy part actually has a different shape, as can be seen in the attached figure.

We used this code in the case of synchrad:

def main():
    """Main entry point."""
    args = cmdline_args()

    calc_input = {
        "grid": [
            (args.params["omega_min"], args.params["omega_max"]),
            (args.params["theta_min"], args.params["theta_max"]),
            (PHI_MIN, PHI_MAX),
            (N_OMEGA, N_THETA, N_PHI),
        ],
        "dtype": "double",
        "native": False,
        "ctx": "mpi",  # set your context to avoid being asked
    }

    calc = SynchRad(calc_input)
    calc.calculate_spectrum(
        file_tracks=args.infile,
        comp='cartesian_complex',
        file_spectrum=args.outfile,
    )


if __name__ == "__main__":
    main()

adapted from the Betatron_Example.ipynb.

Issue when running example/example_tests.ipynb

Hi! Is this repo still maintained? I could successfully install SynchRad, but I got the following error in cell 7 of the notebook (all previous cells executed successfully):

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
/var/folders/2b/8zpydy_x75s8h46pfd2x5_jc0000gp/T/ipykernel_37319/3885437814.py in <module>
      1 t0 = time.time()
----> 2 calc.calculate_spectrum(particleTracks, timeStep=dt, 
      3                         comp='total', Np_max=Np, nSnaps=24,
      4                         it_range=[0, particleTracks[0][0].size])
      5 

~/miniconda3/envs/synchrad/lib/python3.9/site-packages/synchrad-0.0.1-py3.9.egg/synchrad/calc.py in calculate_spectrum(self, particleTracks, file_tracks, timeStep, comp, Np_max, nSnaps, it_range, file_spectrum, verbose)
     42 
     43         nSnaps = np.uint32(nSnaps)
---> 44         self._init_raditaion(comp, nSnaps)
     45 
     46         if timeStep is not None:

~/miniconda3/envs/synchrad/lib/python3.9/site-packages/synchrad-0.0.1-py3.9.egg/synchrad/calc.py in _init_raditaion(self, comp, nSnaps)
    291         for vec_comp in vec_comps[comp]:
    292             self.Data['radiation'][vec_comp] = \
--> 293                 arrcl.zeros(self.queue, radiation_shape, dtype=self.dtype)
    294 
    295     def _init_data(self):

AttributeError: 'SynchRad' object has no attribute 'queue'

Please let me know if there is another tool that you use preferably. Thanks!