ComputationalRadiationPhysics / clara2

Clara2 - a parallel classical radiation calculator based on Liénard-Wiechert potentials
GNU General Public License v3.0
13 stars 9 forks source link

Setting up a first Clara2 simulation #96

Open PrometheusPi opened 6 years ago

PrometheusPi commented 6 years ago

This is a continuation of issue #93 Since the discussion changed from how to set up the environment on tianhe2 to how to set up Clara2 in general, I think a new issue should be used (to avid scrolling down through all the issues related to the module system on tinahe2)

This question was asked by @QJohn2017 in #93


In the beginning of settings.hpp you find the basic parameters:

PrometheusPi commented 6 years ago

In the input section of settings.hpp the data in the trajectory file is defined. By default a GPT output trajectory is assumed that looks like:

x_position[meter]    y_position[meter]    z_position[meter]   beta_x[unitless]  beta_y[unitless]   beta_z[unitless]   time[seconds] 

with beta being the velocity normalized to the speed of light and beta_{x,y,z} being its vector components.

Each column contains the value specified above. Each row contains a new point of the trajectory with an increased time.

This results in the following input settings:

namespace input
{
  const unsigned int index_time = 6; /* column id of time */

  const unsigned int index_pos_x = 0; /* column id of x */
  const unsigned int index_pos_y = 1; /* column id of y */
  const unsigned int index_pos_z = 2; /* column id of z */

  const unsigned int index_beta_x = 3; /* column id of beta_x */
  const unsigned int index_beta_y = 4; /* column id of beta_y */
  const unsigned int index_beta_z = 5; /* column id of beta_z */

  /* unit conversion to SI units */
  const double convert_time = 1.0; /* s -> s */
  const double convert_pos = 1.0; /* m -> m */
  const double convert_beta = 1.0; /* None -> None */
}

If you for example prefere cgs units and the trajectory you have is in units of

PrometheusPi commented 6 years ago

If you have a more complex input set, for example if you use momentum values instead of velocities, you should have a look at run_trough_data.hpp. There you could add more complex (non-linear) conversions for your specific trajectory data format.

PrometheusPi commented 6 years ago

The last part of settings.hpp is for the post- processing via process defined here. It collects the radiation of all trajectories and adds them up, to generate a total (incoherent) spectra.

Nothing needs to be changed there.

TheresaBruemmer commented 6 years ago

This last part is what does not work in my case. The single trace spectra are not added up. But I would expect a "could not write output file" ... but I get no error at all...

PrometheusPi commented 6 years ago

Regarding the observation points:

The number of total directions is the product of N_theta and N_phi defined in settings.hpp.

Their individual values in degree are defined here. Be default phi has the values 0 degree and 90 degree.

These values are handed to the single_direction.cpp code for each combination individually here.

In the single_direction.cpp code, the actual unit vector pointing towards the observer is defined here. First of all, all angles are converted to radiant. And then the vector is computed. If you want to change the observation direction, please adjust the code of the vector. By default, the vector is

looking_vector = R_vec(std::cos(angle_theta) ,
                         std::sin(angle_theta) * std::cos(angle_phi)  ,
std::sin(angle_theta) * std::sin(angle_phi) );

so for phi = 0 degree this is (cos(theta), sin(theta), 0) while for phi = 90 degreethis is (cos(theta), 0, sin(theta)). This observer points for theta=0 toward the +x direction (1,0,0) regardless of phi.

PrometheusPi commented 6 years ago

@TheresaBruemmer In order to perform the post-processing you need to run

./process_data

after the parallel job has finished.

The parallel job should have created 3 types of files:

Did you run ./process_data?

Only this last step creates the my_spectrum_all***.txt files.

TheresaBruemmer commented 6 years ago

Nice, thank you. It works!

PrometheusPi commented 6 years ago

In a future version, where no array jobs are supported any more, this step can be executed by clara2 directly. But with array job support, we need a separate executable for this step, because the individual array jobs are not aware of each other. With MPI only, this would not be required.

PrometheusPi commented 6 years ago

@QJohn2017 Did setting up your simulation work for you?

QJohn2017 commented 6 years ago

@PrometheusPi Sorry that I'm busy these days and have no time to set my simulation, i will have a try to create the trajectory files and set the simulation again tonight. My previous trajectory files are generated by the Vsim code and only has the position information of electrons that are being tracked.

QJohn2017 commented 6 years ago

@PrometheusPi , I have set up the first simulation successfully, when I input the command./process_data, there are two output files named my_spectrum_all_000.dat and my_spectrum_all_001.dat created. But I don't know how these two files stand for what and how to analyse them

PrometheusPi commented 6 years ago

@QJohn2017 No problem - thanks for the reply. Great to hear that running Clara2 worked with your data :+1:. That are wonderful news.

The two files ./process_data creates are the directionally resolved spectra for both phi values you used. So in case you were using the default setup my_spectrum_all_000.dat is the phi=0° spectrum and my_spectrum_all_001.dat is the phi=90° spectrum. Both files are text files and are a matrix with columns describing each frequency omega values and the rows describing your theta angle values.

You can load these files directly into your favorite data anaysis/visualization framework like Mathematica, Python (numpy, matplotlib, ...), MatLab, etc.

There is also a (python based) tool provided in clara2/tools/plotRadiation. If you run

./clara2/tools/plotRadiation my_spectrum_all_000.dat

it will plot the directionally resolved spectrum for phi=0° to your screen (if logged in with X-Window via ssh -X ...). Without arguments this uses default values for theta and omega. There is an extensive help that allows you to specify these values as well. You can access this help by running ./plotRadiation -h. If you for example computed radiation from 0 to 100 Hz, and your theta angle goes from -10° to +30°, and you want to produce a pdf instead of a window on your screen, you could run:

./clara2/tools/plotRadiation --dataExtend 0 100 -10 30  --labelOmega "$\omega\ [Hz]$" --pdf my_spectrum_all_000.dat 

If you have questions how to load these output files into python, please feel free to ask me. With all other post processing tools I am unfortunately not that familiar.

@QJohn2017 Did the spectra came out as you expected?

QJohn2017 commented 6 years ago

@PrometheusPi , Why only the first column in spectra file my_spectrum_all_000.dat and my_spectrum_all_001.dat is not zero and the rest of the columns are all zero ?

QJohn2017 commented 6 years ago

@PrometheusPi , When I input the following code, it returns following errors: [ac_siom_jsliu_1@ln1%tianhe2-C src]$ ../tools/plotRadiation my_spectrum_all_000.dat Traceback (most recent call last): File "../tools/plotRadiation", line 53, in <module> matplotlib.__version__, tested_matplotlib_version) ValueError: zero length field name in format [ac_siom_jsliu_1@ln1%tianhe2-C src]$

PrometheusPi commented 6 years ago

@QJohn2017 That might have a variety of reasons. Could you post your line 26 to 32 of your settings.hpp file here? That would help me a lot. Furthermore it would be helpful to know:

My first guess would be that your transversal angular resolution is unsuited for the high frequencies you want to calculate. For example if you have highly relativistic particles moving in direction (0,1,0), and you resolve angles larger than 1/gamma the blue shift of the radiation reduces drastically and only low frequency radiation would occur at large angles which might not be resolved by the fft sampling. But in order to determine that, I need a bit more information on what you simulate (see questions above).

PrometheusPi commented 6 years ago

@QJohn2017 The error from the plotRadiation comes from a non-tested matplotlib version.That might not be critical, but just a python (version) issue.

These values can be determined by typing:

python --version
python
> import matplotlib
> print(matplotlib.__version__)
> quit()
QJohn2017 commented 6 years ago

@PrometheusPi , The line 26 to 32 of my settings.hpp file are:

#pragma once

namespace param
{
const double omega_max                 = 3.0e19;      /* maximum of plotted frequency Hz */
const double theta_max                 = 1.14594939;  /* maximum of theta in degree */
const unsigned int N_spectrum          = 2048; /* number of frequencies "omega" */
const unsigned int N_theta             = 120;     /* number of directions in first angle "theta" */
const unsigned int N_phi               = 2;         /* number of directions in second angle "phi" */
const unsigned int N_trace             = 199;    /* maximum number of traces */

const unsigned int fft_length_factor   = 1; /* needs to be a power of two */

const bool ascii_output = false; /* output spectra as ascii text */

const unsigned int N_char_filename=256; // number of characters
/* template for input traces (using C formatting): */
const char traceFileTemplate[] = "outfile/trace_%d.dat";
/* template for output spectra for each trace (using C formatting): */
const char outputFileTemplate[] = "my_spectrum_trace%06d.dat";

namespace input
{

The direction of my trajectories move mostly is z direction. The trajectories are some electrons move in an magnet undulator and the magnetic direction is y direction. The relativistic factor of electrons is 323. I don't know the phi and theta stand for what directions in your code.

QJohn2017 commented 6 years ago

@PrometheusPi , When I typing the command as your follows, it returns follows:

[ac_siom_jsliu_1@ln1%tianhe2-C src]$ python
Python 2.6.6 (r266:84292, Sep  4 2013, 07:46:00)
[GCC 4.4.7 20120313 (Red Hat 4.4.7-3)] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import matplotlib
>>> print (matplotlib._version_)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
AttributeError: 'module' object has no attribute '_version_'
>>>
PrometheusPi commented 6 years ago

@QJohn2017 You mentioned that you are using Vsim. I assume this is Vsim from Tech-X. That is a quite sophisticated particle-in-cell code. If you are performing laser plasma simulations like LWFA, TNSA, etc. you might be interested in switching from the Vsim + Clara2 workflow and use the particle-in-cell code PIConGPU instead. It comes with a so-called radiation plugin that performs the same calculation as clara2 but directly as a diagnostic within the particle-in-cell code. With the recent use of ALPAKA, this particle-in-cell code can also run on CPU only clusters as thianhe2. However, if you depend on some of the more advanced material implementations or boundary conditions, Vsim + clara2 is probably the better choice.

PrometheusPi commented 6 years ago

@QJohn2017 Could you use two underscores _ before and after matplotlib.__version__ please? The variable with only one underscore is not defined.

QJohn2017 commented 6 years ago

@PrometheusPi , My research field is LWFA and my main simulation tool is VSim. I never use the PIConGPU code, is it easy to get started?

In this clara2 simulation work, the trajectories used in the simulation are produced by a C++ program developed by myself.

QJohn2017 commented 6 years ago

@PrometheusPi , After enter into the python environment, I type the code import matplotlib and print (matplotlib._version_). Is it right?

QJohn2017 commented 6 years ago

@PrometheusPi ,Sorry, I just typed only one underscores _ before and after matplotlib.__version__. After I typed two underscores, it retunrs:

>>> print (matplotlib.__version__)
0.99.1.1
>>>
QJohn2017 commented 6 years ago

@PrometheusPi , The version of thematplotlib is 0.99.1.1.

PrometheusPi commented 6 years ago

@QJohn2017 That is a matplotlib version from August 2009. Sadly there is no documentation for such old matplotlib versions easily available online. From the error message I assume that at this pre-1.0-relese there was not yet a __version__ variable available (or it was defined differently).

Are the more modern python versions as modules available? Best would be python 2.7 or 3.0 to 3.4 with a matplotlib version 1.2.0 or higher. Ideally a recent version 2.0.x would be best. You could test this by typing:

module av python

If there is no more recent version, just comment-out the if-case on line 49 to 54 of plotRadiation by adding a # as a first character. This deactivates this safety check, but I expect other errors caused by the older version.

PrometheusPi commented 6 years ago

@QJohn2017 PIConGPU is definitely not easy to program, but the ease-of-use has improved in recent years. It is not as easy to set uo es e.g. EPOCH, but there is just a few parameters files you need to know, when performing LWFA simulations since this is one of the standard applications of this code.

Since I never used Vsim myself, I have no direct comparison. But you could try to run a simple LWFA default example, to get a felling whether using PIConGPU would be an alternative. One remark in advance, it is capable of providing radiation spectra for billions of macro particles but it will not provide trajectories easily.

PrometheusPi commented 6 years ago

@QJohn2017 Could you post your detector setup from line 67 to 73 from single_direction.cpp? That would allow me to determine in which directions the detectors are set up.

QJohn2017 commented 6 years ago

@PrometheusPi , I have downloaded the spectra file my_spectrum_all_000.dat and my_spectrum_all_001.dat to my own local computer and processed them successfully by the python in my own computer.

QJohn2017 commented 6 years ago

@PrometheusPi , Okay, Thank you ! I will try to use the PIConGPU code in tianhe2 once all the problems in CLARA2 have been solved.

QJohn2017 commented 6 years ago

@PrometheusPi , the detector setup from line 67 to 73 from single_direction.cpp as follows:

  double angle_theta;
  double angle_phi;

  angle_theta = theta_offset/180.0 * M_PI;

  angle_phi =  phi_offset/180.0 * M_PI;

  looking_vector = R_vec(std::cos(angle_theta) ,
                         std::sin(angle_theta) * std::cos(angle_phi)  ,
                         std::sin(angle_theta) * std::sin(angle_phi) );
QJohn2017 commented 6 years ago

@PrometheusPi , I found that both of the spectrum files my_spectrum_all_000.datand my_spectrum_all_001.dat have of a size of 120x2048. So, I guess the result in my simulation that only the first column is not zero but all of the rest columns are zeros can be attributed to the frequency. May be the frequency interval my_delta_omega is too large and only the base frequency is valid for simulation.

QJohn2017 commented 6 years ago

@PrometheusPi , But when I changed the maximum frequency to 3.0e16, the result is still the same as above. So, i don't know what the problem is.

PrometheusPi commented 6 years ago

@QJohn2017 Great to hear that the plotRadiation tool works on your computer.

Regarding your detector setup: With the current implementation, your detectors are located in a narrow angle of ~1.14 degree around the x-direction looking_vector = (1, 0, 0), but your particles propagate in z-direction beta ~ (0,0,1). Due to beaming, nearly no radiation reaches your detector.
If you place your detector around the z-direction by:

  looking_vector = R_vec(std::sin(angle_theta) * std::sin(angle_phi)  ,
                         std::sin(angle_theta) * std::cos(angle_phi)  ,
                         std::cos(angle_theta) );

You observe the dominant radiation in z-direction. For the above case, phi=0 degree (index 000) is in the plane of the magnetic field, while phi=90 degree is orthogonal to the magnetic field and the propagation direction. After recompilation and rerunning the simulation, the new setup should show more radiation.

Regarding your frequency setup: Since you are simulating an undulator (of undulator period lambda_u) with electrons of Lorentz factor gamma=323 the expected radiation frequency is

# radiation wave length
lambda_rad = lambda_u / (2 * gamma^2)
# radiation frequency
omega_rad = 2 * pi * c / lambda_rad = 4 * gamma^2 * pi * c / lambda_u

So if your undulator period lambda_u is much larger than 13 micrometer you should resolve the radiation. (Please be aware that the above equation is only valid for the fundamental at low k (a_0) values assuming an undulator and not a wiggler.)

@QJohn2017 Did the change in observation direction result in a more reasonable spectrum?

Regarding PIConGPU: If you are planning to use PIConGPU, feel free to ask for support via GitHub by opening an issue in the main repo. If you are simulating the LWFA in combination with the undulator with Vsim, this is possible with PIConGPU if the undulator is in vacuum after the LWFA stage. But if you are simulating an undulator that surrounds the plasma cavity, fiddling around with the exact boundary conditions is still not easy. Adding more advanced boundary conditions is a current project and might be possible in a couple of month.

PrometheusPi commented 6 years ago

The problem with the frequencies should be resolved with the changed observation direction.

With the current observation direction in x-direction, the observable radiation is not blue shifted and any frequency in x-direction will be gamma^2 ~ 100000 times smaller than the expected frequency of the undulator radiation. Since you resolve only 2048 frequency points, this low radiation is lost (due to the interpolation after the fft).

@QJohn2017 How large is your undulator period lambda_u?

TheresaBruemmer commented 6 years ago

Hi, I tried running this clara2 version on hypnos, using trace files I already ran with an old clara2 version. I also got a few non-zero values at the beginning the rest all zero. So in my case, the trajectories are not the problem...maybe it is similar to @QJohn2017

PrometheusPi commented 6 years ago

@TheresaBruemmer Might it be that you are still using the default observation direction?

QJohn2017 commented 6 years ago

@PrometheusPi , Thank you very much that according to your suggestions in setting observation direction, I have got the more reasonable spectrum. When I postprocess the spectrum in my own computer by typing the command as following, it returns:

Kevin@HP /cygdrive/f/Study/python/work/clara2
$ ./plotRadiation --dataExtend 0 3.0E16 0 0.5  --labelOmega "$\omega\ [Hz]$" --pdf my_spectrum_all_001.dat
my_spectrum_all_001.dat  Done
./plotRadiation:54: UserWarning: The matplotlib version differs from that used during
                     development. This might cause the output to appear different.
                     Version used: 2.0.2, tested versions ['1.2.0', '1.3.1']
  warnings.warn(warning_msg)

Kevin@HP /cygdrive/f/Study/python/work/clara2
$

But there is a pdf picture file produced, as the attached file showed. Did the radiation result is right? SpectraOutput.pdf

QJohn2017 commented 6 years ago

@PrometheusPi , The undulator period lambda_u in my simulation is 0.05m and the maximum magnetic is 0.425T.

PrometheusPi commented 6 years ago

In principle your spectrum looks correct - containing higher harmonics (as expected for this high magnetic field) and a blue shift on axis.

An offline discussion with @steindev, our FEL (and undulator) expert, however revealed a discrepancy by a factor 2:

This could only be explained with an undulator period of 10cm instead of 5cm. Or with a bit less electron momentum. Or with a much higher magnetic field. Could this be possible in the Vsim simulation?

I will test the current clara2 version against a known trajectory just to be sure the error is not caused by some new changes in clara2 (to be precise: by the use of fftw, which I did not test properly due to the pre-compiled libraries) .

TheresaBruemmer commented 6 years ago

Ok, sorry, I used the wrong looking vector / detector direction.

TheresaBruemmer commented 6 years ago

Why is the default looking vector in x? I would think z is the more prominent propagation direction ...

QJohn2017 commented 6 years ago

@PrometheusPi , Sorry, I forgot to tell you that in addition to the magnet, there is a seeding laser in co moving with the electrons in the undulator. So, the total parameters are: B0 = 0.425T, lambda_u=0.05m, N_period=20, seeding laser wavelength=1047nm, seed power=1.0e9 W .

PrometheusPi commented 6 years ago

@TheresaBruemmer I do not remember why I choose the x-direction as default. Might be related to the original radiation calculation with trajectories from illumination (see PhD thesis) PIC simulations. But you are right, more common in accelerator physics is the z-direction as main propagation direction. But for me, as a PIConGPU developer, y is the default propagation direction 😉 .

TheresaBruemmer commented 6 years ago

LOL, ok. Not hard to change though, so it is fine... ...maybe I am blind, but I cannot find it: can you tell me where I can change the trace_id to start from 1 instead of 0? so that for N_trace=2000 it will read trace_0001.txt to trace_2000.txt

QJohn2017 commented 6 years ago

@PrometheusPi , However, when I removed the seeding laser by setting the seed power to 0, the results looks no difference than above, as showed in attached file. SpectraOutput.pdf

PrometheusPi commented 6 years ago

@QJohn2017 Please let me run a check on the current dev. I will verify whether the factor two comes from clara2. This might need up to 1-2 hours. I will report back as soon as possible. Perhaps you can check your code in the mean time. Or you wait till tomorrow. It must be already quite late at your place.

@steindev also remarked that a seed laser that co-propagates should have no influence.

PrometheusPi commented 6 years ago

@TheresaBruemmer The main clara2 code does not yet support starting from another value than 0. But it also does not care. It just prints a warning and continues.

On the other hand, the process_data code allows adjusting the start index by adjusting this line of code.

Did you encounter an error with a missing 0th trajectory or just a warning?

TheresaBruemmer commented 6 years ago

@PrometheusPi Yes, I did find that part for process_data. Ok so, if I change N_trace to 1 more, should I then change index_files_last to N_trace-1 or does process_data also not care?

PrometheusPi commented 6 years ago

@TheresaBruemmer index_files_last is the real last index, not the number of spectra. see here

However, both executables just produce warnings and continue if they do not fine files. So just choose them large enough. That will work.

PrometheusPi commented 6 years ago

@QJohn2017 Sorry it might take a bit longer