diff --git a/.gitignore b/.gitignore index 4f85ed2da..cf8f3579b 100755 --- a/.gitignore +++ b/.gitignore @@ -33,3 +33,4 @@ exec_script.sh.* .DS_Store tables/* *.stats* +.venv* diff --git a/benchmarks/collisions/ionization_multiple.py b/benchmarks/collisions/ionization_multiple.py index 66b2f75bf..6313a72dd 100755 --- a/benchmarks/collisions/ionization_multiple.py +++ b/benchmarks/collisions/ionization_multiple.py @@ -22,7 +22,7 @@ timestep2 = int(np.double(S2.namelist.Main.timestep)) D += [ - S2.ParticleBinning(0,sum={"ekin":[0,1]}, linestyle="", marker=".", color=color ) + S2.ParticleBinning(0,sum={"ekin":[0,1]}, linestyle="", marker=".", color=color, label=elm ) ] diff --git a/benchmarks/tst_collisions3_AM_uniformity.py b/benchmarks/tst_collisions3_AM_uniformity.py new file mode 100644 index 000000000..5c6b63f13 --- /dev/null +++ b/benchmarks/tst_collisions3_AM_uniformity.py @@ -0,0 +1,87 @@ +# --------------------------------------------- +# SIMULATION PARAMETERS FOR THE PIC-CODE SMILEI +# --------------------------------------------- + +import math +L0 = 2.*math.pi # conversion from normalization length to wavelength + + +Main( + geometry = "AMcylindrical", + + number_of_patches = [ 4, 4 ], + + interpolation_order = 2, + number_of_AM = 1, + + timestep = 1 * L0, + simulation_time = 10 * L0, + + time_fields_frozen = 100000000000., + + cell_length = [L0, L0], + grid_length = [64.*L0, 64*L0], + + EM_boundary_conditions = [ ["periodic", "periodic"], ["silver-muller", "buneman"] ], + + solve_poisson = False, + + reference_angular_frequency_SI = L0 * 3e8 /1.e-6, +) + + +ion = "ion" +eon = "eon" + +Species( + name = eon, + position_initialization = "random", + momentum_initialization = "cold", + particles_per_cell= 100, + mass = 1.0, + charge = -1.0, + charge_density = 1., + mean_velocity = [0.8, 0., 0.], + temperature = [0.]*3, + time_frozen = 100000000.0, + boundary_conditions = [ + ["periodic", "periodic"], + ["remove", "remove"], + ], +) + +Species( + name = ion, + position_initialization = "random", + momentum_initialization = "cold", + particles_per_cell= 100, + mass = 1836.0*13., + charge = 3.0, + charge_density = 1., + mean_velocity = [0., 0., 0.], + temperature = [0.]*3, + time_frozen = 100000000.0, + boundary_conditions = [ + ["periodic", "periodic"], + ["remove", "reflective"], + ], + atomic_number = 13 +) + +Collisions( + species1 = [eon], + species2 = [ion], + coulomb_log = 3, +) + +def radius(particles): return np.sqrt(particles.y**2 + particles.z**2) + +DiagParticleBinning( + deposited_quantity = "weight_vy_py", + every = 4, + species = [ion], + axes = [ + ["x", 0, Main.grid_length[0], 64], + [radius, 0, Main.grid_length[1], 64] + ] +) diff --git a/doc/Sphinx/Overview/material.rst b/doc/Sphinx/Overview/material.rst index 5712bfeae..07993d322 100644 --- a/doc/Sphinx/Overview/material.rst +++ b/doc/Sphinx/Overview/material.rst @@ -30,7 +30,7 @@ Papers involving Smilei ^^^^^^^^^^^^^^^^^^^^^^^^ Only papers published in peer-reviewed journals are listed (for the complete list of citing papers see `Google Scholar `_). -As of May 2024, at least 192 papers have been published covering a broad range of topics: +As of October 2024, at least 211 papers have been published covering a broad range of topics: * laser-plasma interaction (LPI) / inertial fusion (FCI) * ultra-high intensity (UHI) applications @@ -51,6 +51,127 @@ Following is the distribution of these topics in the listed publications up to N You can count the number of papers in the list with the vim command :%s/.. \[//gn. +.. [Jirka2024] + + M. Jirka and S. V. Bulanov, + `Effects of Colliding Laser Pulses Polarization on e- e+ Cascade Development in Extreme Focusing`, + `Physical Review Letters 133, 125001 (2024) `_ + +.. [Plotnikov2024] + + I. Plotnikov, A. J. van Marle, C. Guépin, A. Marcowith and Pierrick Martin, + `Kinetic simulations of electron–positron induced streaming instability in the context of gamma-ray halos around pulsars`, + `Astronomy & Astrophysics 688, A134 (2024) `_ + +.. [Mukherjee2024] + + A. Mukherjee and Daniel Seipt, + `Laser polarization control of ionization-injected electron beams and x-ray radiation in laser wakefield accelerators`, + `Plasma Physics and Controlled Fusion 66, 085001 (2024) `_ + +.. [Yao2024] + + W. Yao, R. Lelièvre, T. Waltenspiel, I. Cohen, A. Allaoua, P. Antici, A. Beck, E. Cohen, X. Davoine, E. d’Humières, Q. Ducasse, E. Filippov, C. Gautier, L. Gremillet, P. Koseoglou, D. Michaeli, D. Papadopoulos, S. Pikuz, I. Pomerantz, F. Trompier, Y. Yuan, F. Mathieu and Julien Fuchs, + `Enhanced Energy, Conversion Efficiency and Collimation of Protons Driven by High-Contrast and Ultrashort Laser Pulses`, + `Applied Sciences 14, 6101 (2024) `_ + +.. [Dmitriev2024] + + E. Dmitriev and Ph. Korneev, + `Angular momentum gain by electrons under the action of intense structured light`, + `Physical Review A 110, 013514 (2024) `_ + +.. [Yu2024] + + H. Yu, Q. Xia and Jun Fang, + `Nonthermal Acceleration of Electrons, Positrons, and Protons at a Nonrelativistic Quasi-parallel Collisionless Shock`, + `The Astrophysical Journal 969, 13 (2024) `_ + + +.. [Liu2024] + + B. Liu, B. Lei, Y. Gao, M. Wen, and K. Zhu, + `Plasma opacity induced by laser-driven movement of background ions`, + `Plasma Physics and Controlled Fusion 66, 115004 (2024) `_ + +.. [Martin2024] + + H. Martin, R. W. Paddock, M. W. von der Leyen, V. Eliseev, R. T. Ruskov, R. Timmis, J. J. Lee, A. James, and P. A. Norreys, + `Electrothermal filamentation of igniting plasmas`, + `Physical Review E 110, 035205 (2024) `_ + +.. [Horny2024] + + V. Horný, P. G. Bleotu, D. Ursescu, V. Malka, and P. Tomassini, + `Efficient laser wakefield accelerator in pump depletion dominated bubble regime`, + `Physical Review E 110, 035202 (2024) `_ + +.. [DeAndres2024] + + A. De Andres, S. Bhadoria, J. T. Marmolejo, A. Muschet, P. Fischer, H. R. Barzegar, T. Blackburn, A. Gonoskov, D. Hanstorp, M. Marklund and L. Veisz, + `Unforeseen advantage of looser focusing in vacuum laser acceleration`, + `Nature Communications Physics 7, 293 (2024) `_ + +.. [Yoon2024] + + Y. D. Yoon, M. Laishram, T. E. Moore, and G. S. Yun, + `Non-equilibrium formation and relaxation of magnetic flux ropes at kinetic scales`, + `Nature Communications Physics 7, 297 (2024) `_ + +.. [Laishram2024] + + M. Laishram, G. S. Yun, and Y. D. Yoon, + `Magnetogenesis via the canonical battery effect`, + `Physical Review Research 6, L032052 (2024) `_ + +.. [Kang2024] + + H. L. Kang, Y. D. Yoon, M.-H. Cho and G. S. Yun, + `Fast nonlinear scattering of runaway electron beams through resonant interactions with plasma waves`, + `Nuclear Fusion 64, 10 (2024) `_ + +.. [Gonzalez-Izquierdo2024] + + B. Gonzalez-Izquierdo, P. Fischer, M. Touati, J. Hartmann, M. Speicher, V. Scutelnic, G. Bodini, A. Fazzini, M. M. Guenther, A. K. Haerle, K. Kenney, E. Schork, S. Bruce, M. Spinks, H. J. Quevedo, A. Helal, M. Medina, E. Gaul, H. Ruhl, M. Schollmeier, S. Steinke, and G. Korn, + `Efficient laser-driven proton acceleration from a petawatt contrast-enhanced second harmonic mixed-glass laser system`, + `Physics of Plasmas 31, 083105 (2024) `_ + +.. [Ma2024] + + Z. Ma, Y. Wang, Y. Yang, Y. Wang, K. Zhao, Y. Li, C. Fu, W. He, Y. Ma, + `Simulation of nuclear isomer production in laser-induced plasma`, + `Matter and Radiation at Extremes 9, 055201 (2024) `_ + +.. [Kleij2024] + + P. S. Kleij, S. Marini, M. Caetano de Sousa, M. Grech, C. Riconda, M. Raynaud, + `Photon emission and radiation reaction effects in surface plasma waves in ultra-high intensities`, + `Physics of Plasmas 31, 072111 (2024) `_ + +.. [Ghizzo2024] + + A. Ghizzo, D. Del Sarto, H. Betar, + `Collisionless heating in Vlasov plasma and turbulence-driven filamentation aspects`, + `Physics of Plasmas 31, 072109 (2024) `_ + +.. [Gelfer2024] + + E. G. Gelfer, A. M. Fedotov, O. Klimo, and S. Weber, + `Coherent radiation of an electron bunch colliding with an intense laser pulse`, + `Physical Review Research 6, L032013 (2024) `_ + +.. [Zagidullin2024] + + R. Zagidullin, V. Zorina, J. W. Wang, S. G. Rykovanov, + `Polarization control of attosecond pulses from laser-nanofoil interactions using an external magnetic field`, + `Physics of Plasmas 31, 073303 (2024) `_ + +.. [Marini2024] + + S. Marini, D. F. G. Minenna, F. Massimo, L. Batista, V. Bencini, A. Chancé, N. Chauvin, S. Doebert, J. Farmer, E. Gschwendtner, I. Moulanier, P. Muggli, D. Uriot, B. Cros, and P. A. P. Nghiem, + `Beam physics studies for a high charge and high beam quality laser-plasma accelerator`, + `Physical Review Accelerators and Beams 27, 063401 (2024) `_ + .. [Sikorski2024] P. Sikorski, A. G. R. Thomas, S. S. Bulanov, M. Zepf and D. Seipt, @@ -385,7 +506,7 @@ Following is the distribution of these topics in the listed publications up to N V. Istokskaia, M. Tosca, L. Giuffrida, J. Psikal, F. Grepl, V. Kantarelou, S. Stancek, S. Di Siena, A. Hadjikyriacou, A. McIlvenny, Y. Levy, J. Huynh, M. Cimrman, P. Pleskunov, D. Nikitin, A. Choukourov, F. Belloni, A. Picciotto, S. Kar, M. Borghesi, A. Lucianetti, T. Mocek and D. Margarone, `A multi-MeV alpha particle source via proton-boron fusion driven by a 10-GW tabletop laser`, - `Communications Physics 6, 27 (2023) `_ + `Nature Communications Physics 6, 27 (2023) `_ .. [Yoon2023] @@ -428,6 +549,12 @@ Following is the distribution of these topics in the listed publications up to N S. Marini, M. Grech, P. S. Kleij, M. Raynaud and C. Riconda, `Electron acceleration by laser plasma wedge interaction`, `Physical Review Research 5, 013115 (2023) `_ + +.. [Miloshevsky2023] + + G. Miloshevsky, + `Particle-in-Cell Modeling of Omega Experiments on Ablation of Plasmas`, + `IEEE Transactions on Plasma Science 51, 4 (2023) `_ .. [Blackman2022] @@ -507,7 +634,7 @@ Following is the distribution of these topics in the listed publications up to N `Injection of electron beams into two laser wakefields and generation of electron rings`, `Physical Review E 106, 055202 (2022) `_ -.. [Kumar2022b] +.. [Ku2022] S. Ku., R. Dhawan, D.K. Singh and H. K. Malik, `Diagnostic of laser wakefield acceleration with ultra – Short laser pulse by using SMILEI PIC code`, @@ -519,12 +646,6 @@ Following is the distribution of these topics in the listed publications up to N `Comparative study of ultrashort single-pulse and multi-pulse driven laser wakefield acceleration`, `Laser Physics Letters 20, 026001 (2022) `_ -.. [Miloshevsky2022] - - G. Miloshevsky, - `Pic Modeling of Omega Experiments on Ablation of Plasmas`, - `2022 IEEE International Conference on Plasma Science (ICOPS), ICOPS45751.2022.9813047 (2022) `_ - .. [Zhang2022b] Y. Zhang, F. Wang, J. Liu and J. Sun, @@ -628,12 +749,6 @@ Following is the distribution of these topics in the listed publications up to N `Subcycle terahertz field waveforms clocked by attosecond high-harmonic pulses from relativistic laser plasmas`, `Journal of Applied Physics 131, 103104 (2022) `_ -.. [Umstadter2022] - - D. Umstadter - `Controlled Injection of Electrons for Improved Performance of Laser-Wakefield Acceleration`, - `United States Department of Energy Technical Report (2022) `_ - .. [Massimo2022] F. Massimo, M. Lobet, J. Derouillat, A. Beck, G. Bouchard, M. Grech, F. Pérez, T. Vinci, @@ -880,7 +995,7 @@ Following is the distribution of these topics in the listed publications up to N W. Yao, A. Fazzini, S. N. Chen, K. Burdonov, P. Antici, J. Béard, S. Bolaños, A. Ciardi, R. Diab, E. D. Filippov, S. Kisyov, V. Lelasseux, M. Miceli, Q. Moreno, V. Nastasa, S. Orlando, S. Pikuz, D. C. Popescu, G. Revet, X. Ribeyre, E. d’Humières and J. Fuchs, `Laboratory evidence for proton energization by collisionless shock surfing`, - `Nat. Phys. 17, 1177-1182 (2021) `_ + `Nature Physics 17, 1177-1182 (2021) `_ .. [Gelfer2021] diff --git a/doc/Sphinx/Overview/releases.rst b/doc/Sphinx/Overview/releases.rst index 9e2b1445b..85afde2aa 100755 --- a/doc/Sphinx/Overview/releases.rst +++ b/doc/Sphinx/Overview/releases.rst @@ -28,12 +28,14 @@ Changes made in the repository (not released) * Prescribed fields in AM geometry. * Particle reflective boundary conditions at Rmax in AM geometry. * 1st order Ruyten shape function in AM geometry. + * Support for collisions in single mode AM geometry. * **Bug fixes**: * Tunnel ionization was wrong in some cases for high atomic numbers. * Custom functions in ``ParticleBinning`` crashed with python 3.12. * Species-specific diagnostics in AM geometry with vectorization. + * Frozen particles in AM geometry with adaptive vectorization. * Happi's ``average`` argument would sometimes be missing the last bin. * 1D projector on GPU without diagnostics diff --git a/doc/Sphinx/Understand/collisions.rst b/doc/Sphinx/Understand/collisions.rst index 281ff0764..dc89d3e8c 100755 --- a/doc/Sphinx/Understand/collisions.rst +++ b/doc/Sphinx/Understand/collisions.rst @@ -309,21 +309,38 @@ the ion: \sum\limits_{p=0}^{k-1} R^{i+k}_{i+p} \left(\bar{P}^{i+k} - \bar{P}^{i+p}\right) \prod\limits_{j=0,j\ne p}^{k-1} R^{i+p}_{i+j} & - \quad\mathrm{if}\quad 0U`. + + ---- Test cases for ionization diff --git a/doc/Sphinx/Use/namelist.rst b/doc/Sphinx/Use/namelist.rst index 163aa8654..484512179 100755 --- a/doc/Sphinx/Use/namelist.rst +++ b/doc/Sphinx/Use/namelist.rst @@ -1390,8 +1390,8 @@ Lasers ^^^^^^ A laser consists in applying oscillating boundary conditions for the magnetic -field on one of the box sides. The only boundary condition that supports lasers -is ``"silver-muller"`` (see :py:data:`EM_boundary_conditions`). +field on one of the box sides. The only boundary conditions that support lasers +are ``"silver-muller"`` and ``"PML"`` (see :py:data:`EM_boundary_conditions`). There are several syntaxes to introduce a laser in :program:`Smilei`: .. note:: @@ -1440,10 +1440,10 @@ There are several syntaxes to introduce a laser in :program:`Smilei`: .. py:data:: space_time_profile_AM - :type: A list of maximum 2 x ``number_of_AM`` *python* functions. + :type: A list of maximum 2 x ``number_of_AM`` complex valued *python* functions. - These profiles define the first modes of :math:`B_r` and :math:`B_\theta` in the - order shown in the above example. Undefined modes are considered zero. + These profiles define the first modes of :math:`B_r` and :math:`B_\theta` of the laser in the + order shown in the above example. Higher undefined modes are considered zero. This can be used only in ``AMcylindrical`` geometry. In this geometry a two-dimensional :math:`(x,r)` grid is used and the laser is injected from a :math:`x` boundary, thus the provided profiles must be a function of :math:`(r,t)`. diff --git a/doc/Sphinx/_static/workshopLogo.svg b/doc/Sphinx/_static/workshopLogo.svg new file mode 100644 index 000000000..9cb595c1d --- /dev/null +++ b/doc/Sphinx/_static/workshopLogo.svg @@ -0,0 +1,21 @@ + + + + + + + + + + + + + + + + + + diff --git a/doc/Sphinx/index.rst b/doc/Sphinx/index.rst index b3ca0b413..045a35879 100755 --- a/doc/Sphinx/index.rst +++ b/doc/Sphinx/index.rst @@ -12,16 +12,15 @@ Open-source, collaborative, user-friendly and designed for high performances on it is applied to a wide range of physics studies: from relativistic laser-plasma interaction to astrophysics. -.. - .. raw:: html +.. raw:: html -.. - -
- 4th user & training workshop: 8-10 Nov 2023 in Prague (Czechia) + +
+ workshop +
+ 5th user & training workshop:
19-21st March 2025 in Madrid (Spain)
-..
diff --git a/doc/Sphinx/smilei_theme/layout.html b/doc/Sphinx/smilei_theme/layout.html index 1bed82e81..9a19da78a 100755 --- a/doc/Sphinx/smilei_theme/layout.html +++ b/doc/Sphinx/smilei_theme/layout.html @@ -99,6 +99,12 @@ +
+ + workshop + +
{%- for menu in menus %} diff --git a/src/Collisions/CollisionalIonization.cpp b/src/Collisions/CollisionalIonization.cpp index 96d0489b1..3f32af64d 100755 --- a/src/Collisions/CollisionalIonization.cpp +++ b/src/Collisions/CollisionalIonization.cpp @@ -23,7 +23,7 @@ CollisionalIonization::CollisionalIonization( int Z, Params *params, int ionizat { atomic_number = Z; rate .resize( Z ); - irate.resize( Z ); + rate_product.resize( Z ); prob .resize( Z ); ionization_electrons_ = ionization_electrons; if( params ) { @@ -39,7 +39,7 @@ CollisionalIonization::CollisionalIonization( CollisionalIonization *CI ) { atomic_number = CI->atomic_number; rate .resize( atomic_number ); - irate.resize( atomic_number ); + rate_product.resize( atomic_number ); prob .resize( atomic_number ); ionization_electrons_ = CI->ionization_electrons_; new_electrons.initialize( 0, CI->new_electrons ); @@ -177,7 +177,7 @@ void CollisionalIonization::calculate( double gamma_s, double gammae, double gam // Loop for multiple ionization // k+1 is the number of ionizations const int kmax = atomic_number-Zstar-1; - double cs, w, e, cum_prob = 0; + double cs, w, e, cum_prob = 0, A = 1.; for( int k = 0; k <= kmax; k++ ) { // Calculate the location x (~log of energy) in the databases const double x = a2*log( a1*( gamma_s-1. ) ); @@ -203,24 +203,24 @@ void CollisionalIonization::calculate( double gamma_s, double gammae, double gam } rate[k] = K*cs/gammae ; // k-th ionization rate - irate[k] = 1./rate[k] ; // k-th ionization inverse rate prob[k] = exp( -rate[k] ); // k-th ionization probability + rate_product[k] = 1.; // Calculate the cumulative probability for k-th ionization (Nuter et al, 2011) if( k==0 ) { - cum_prob = prob[k]; + cum_prob = prob[k]; // cumulative probability } else { + double sum = 0.; for( int p=0; p rate; - std::vector irate; + std::vector rate_product; //! Current ionization probability array (one cell per number of ionization events) std::vector prob; diff --git a/src/Diagnostic/Histogram.h b/src/Diagnostic/Histogram.h index 49e970869..9ab9cbab3 100755 --- a/src/Diagnostic/Histogram.h +++ b/src/Diagnostic/Histogram.h @@ -567,24 +567,26 @@ class HistogramAxis_user_function : public HistogramAxis public: HistogramAxis_user_function( PyObject *type_object ) : HistogramAxis(), - function( type_object ), - particleData( 0 ) - { - }; + function( type_object ) + {}; ~HistogramAxis_user_function() { + SMILEI_PY_ACQUIRE_GIL Py_DECREF( function ); + SMILEI_PY_RELEASE_GIL }; private: void calculate_locations( Species *s, double *array, int *index, unsigned int npart, SimWindow * ) { + PyArrayObject *ret; SMILEI_PY_ACQUIRE_GIL - // Expose particle data as numpy arrays - particleData.resize( npart ); - particleData.set( s->particles ); - // run the function - PyArrayObject *ret = ( PyArrayObject * )PyObject_CallFunctionObjArgs( function, particleData.get(), NULL ); - particleData.clear(); + { + // Expose particle data as numpy arrays + ParticleData particleData( npart ); + particleData.set( s->particles ); + // run the function + ret = ( PyArrayObject * )PyObject_CallFunctionObjArgs( function, particleData.get(), NULL ); + } // Copy the result to "array" double *arr = ( double * ) PyArray_GETPTR1( ret, 0 ); for( unsigned int ipart = 0 ; ipart < npart ; ipart++ ) @@ -599,7 +601,6 @@ class HistogramAxis_user_function : public HistogramAxis }; PyObject *function; - ParticleData particleData; }; #endif @@ -1155,24 +1156,27 @@ class Histogram_user_function : public Histogram public: Histogram_user_function( PyObject *deposited_quantity_object ) : Histogram(), - function( deposited_quantity_object ), - particleData( 0 ) + function( deposited_quantity_object ) {}; ~Histogram_user_function() { + SMILEI_PY_ACQUIRE_GIL Py_DECREF( function ); + SMILEI_PY_RELEASE_GIL }; private: void valuate( Species *s, double *array, int *index ) { unsigned int npart = s->getNbrOfParticles(); + PyArrayObject *ret; SMILEI_PY_ACQUIRE_GIL // Expose particle data as numpy arrays - particleData.resize( npart ); - particleData.set( s->particles ); - // run the function - PyArrayObject *ret = ( PyArrayObject * )PyObject_CallFunctionObjArgs( function, particleData.get(), NULL ); - particleData.clear(); + { + ParticleData particleData( npart ); + particleData.set( s->particles ); + // run the function + ret = ( PyArrayObject * )PyObject_CallFunctionObjArgs( function, particleData.get(), NULL ); + } // Copy the result to "array" double *arr = ( double * ) PyArray_GETPTR1( ret, 0 ); for( unsigned int ipart = 0 ; ipart < npart ; ipart++ ) { @@ -1186,7 +1190,6 @@ class Histogram_user_function : public Histogram }; PyObject *function; - ParticleData particleData; }; #endif diff --git a/src/Diagnostic/HistogramFactory.h b/src/Diagnostic/HistogramFactory.h index fe3b53651..bb0f8256c 100755 --- a/src/Diagnostic/HistogramFactory.h +++ b/src/Diagnostic/HistogramFactory.h @@ -86,9 +86,13 @@ class HistogramFactory #ifdef SMILEI_USE_NUMPY // Test the function with temporary, "fake" particles double *dummy = NULL; - ParticleData test( params.nDim_particle, deposited_quantity_object, deposited_quantityPrefix, dummy ); - histogram = new Histogram_user_function( deposited_quantity_object ); - histogram->deposited_quantity = "user_function"; + SMILEI_PY_ACQUIRE_GIL + { + ParticleData test( params.nDim_particle, deposited_quantity_object, deposited_quantityPrefix, dummy ); + histogram = new Histogram_user_function( deposited_quantity_object ); + histogram->deposited_quantity = "user_function"; + } + SMILEI_PY_RELEASE_GIL #else ERROR( deposited_quantityPrefix << " should be a string" ); #endif diff --git a/src/ElectroMagn/Laser.cpp b/src/ElectroMagn/Laser.cpp index eb2efee00..299e39c48 100755 --- a/src/ElectroMagn/Laser.cpp +++ b/src/ElectroMagn/Laser.cpp @@ -453,12 +453,9 @@ void LaserProfileSeparable::initFields( Params ¶ms, Patch *patch, ElectroMag double LaserProfileSeparable::getAmplitude( std::vector, double t, int j, int k ) { double amp; - #pragma omp critical - { - double omega = omega_ * chirpProfile_->valueAt( t ); - double phi = ( *phase )( j, k ); - amp = timeProfile_->valueAt( t-( phi+delay_phase_ )/omega ) * ( *space_envelope )( j, k ) * sin( omega*t - phi ); - } + double omega = omega_ * chirpProfile_->valueAt( t ); + double phi = ( *phase )( j, k ); + amp = timeProfile_->valueAt( t-( phi+delay_phase_ )/omega ) * ( *space_envelope )( j, k ) * sin( omega*t - phi ); return amp; } @@ -562,10 +559,7 @@ double LaserProfileFile::getAmplitude( std::vector pos, double t, int j, for( unsigned int i=0; ivalueAt( pos, t ); - } + amp *= extraProfile->valueAt( pos, t ); return amp; } diff --git a/src/ElectroMagn/Laser.h b/src/ElectroMagn/Laser.h index 1a848b6a9..134feff89 100755 --- a/src/ElectroMagn/Laser.h +++ b/src/ElectroMagn/Laser.h @@ -139,7 +139,6 @@ class LaserProfileNonSeparable : public LaserProfile inline double getAmplitude( std::vector pos, double t, int, int ) override { double amp; - #pragma omp critical amp = spaceAndTimeProfile_->valueAt( pos, t ); return amp; } @@ -147,7 +146,6 @@ class LaserProfileNonSeparable : public LaserProfile inline std::complex getAmplitudecomplex( std::vector pos, double t, int, int ) override { std::complex amp; - #pragma omp critical amp = spaceAndTimeProfile_->complexValueAt( pos, t ); return amp; } diff --git a/src/ElectroMagnBC/ElectroMagnBC1D_refl.cpp b/src/ElectroMagnBC/ElectroMagnBC1D_refl.cpp index 9d7518d78..47a5152b1 100755 --- a/src/ElectroMagnBC/ElectroMagnBC1D_refl.cpp +++ b/src/ElectroMagnBC/ElectroMagnBC1D_refl.cpp @@ -35,20 +35,31 @@ void ElectroMagnBC1D_refl::apply( ElectroMagn *EMfields, double, Patch *patch ) { if( i_boundary_ == 0 ) { if( patch->isXmin() ) { + const Field *B[3]{ EMfields->Bx_, EMfields->By_, EMfields->Bz_ }; + double *const __restrict__ By1D = B[1]->data_; + double *const __restrict__ Bz1D = B[2]->data_; // Application over the full-ghost cell - //Field1D* Ex1D = static_cast(EMfields->Ex_); - //Field1D* Ey1D = static_cast(EMfields->Ey_); - //Field1D* Ez1D = static_cast(EMfields->Ez_); - //Field1D* Bx1D = static_cast(EMfields->Bx_); - Field1D *By1D = static_cast( EMfields->By_ ); - Field1D *Bz1D = static_cast( EMfields->Bz_ ); + //Field1D *By1D = static_cast( EMfields->By_ ); + //Field1D *Bz1D = static_cast( EMfields->Bz_ ); +#ifdef SMILEI_ACCELERATOR_GPU_OACC + const int sizeofB1 = B[1]->number_of_points_; + const int sizeofB2 = B[2]->number_of_points_; +#endif // force constant magnetic fields in the ghost cells +#ifdef SMILEI_ACCELERATOR_GPU_OACC + #pragma acc parallel present(By1D[0:sizeofB1],Bz1D[0:sizeofB2]) + #pragma acc loop gang worker vector +#elif defined( SMILEI_ACCELERATOR_GPU_OMP ) + #pragma omp target + #pragma omp teams distribute parallel for +#endif for( unsigned int i=oversize_; i>0; i-- ) { - //(*Bx1D)(i-1) = (*Bx1D)(i); - ( *By1D )( i-1 ) = ( *By1D )( i ); - ( *Bz1D )( i-1 ) = ( *Bz1D )( i ); + //( *By1D )( i-1 ) = ( *By1D )( i ); + //( *Bz1D )( i-1 ) = ( *Bz1D )( i ); + By1D[i-1] = By1D[i]; + Bz1D[i-1] = Bz1D[i]; } // // force 0 electric fields in the ghost cells @@ -60,7 +71,6 @@ void ElectroMagnBC1D_refl::apply( ElectroMagn *EMfields, double, Patch *patch ) // (*Ez1D)(i) = 0.0; // } - /* DEFINITION BY NICO // The other fields are already defined everywhere, so no need for boundary conditions. @@ -76,21 +86,34 @@ void ElectroMagnBC1D_refl::apply( ElectroMagn *EMfields, double, Patch *patch ) }//if Xmin } else { if( patch->isXmax() ) { + const Field *B[3]{ EMfields->Bx_, EMfields->By_, EMfields->Bz_ }; + double *const __restrict__ By1D = B[1]->data_; + double *const __restrict__ Bz1D = B[2]->data_; +#ifdef SMILEI_ACCELERATOR_GPU_OACC + const int sizeofB1 = B[1]->number_of_points_; + const int sizeofB2 = B[2]->number_of_points_; +#endif + const unsigned int nxd = n_d[0]; // application of Bcs over the full ghost cells - //Field1D* Ex1D = static_cast(EMfields->Ex_); - //Field1D* Ey1D = static_cast(EMfields->Ey_); - //Field1D* Ez1D = static_cast(EMfields->Ez_); - //Field1D* Bx1D = static_cast(EMfields->Bx_); - Field1D *By1D = static_cast( EMfields->By_ ); - Field1D *Bz1D = static_cast( EMfields->Bz_ ); + //Field1D *By1D = static_cast( EMfields->By_ ); + //Field1D *Bz1D = static_cast( EMfields->Bz_ ); // force constant magnetic fields in the ghost cells // for (unsigned int i=n_p[0]-oversize_; i b1( n_p[axis1_], 0. ); double *const __restrict__ db1 = b1.data(); - const unsigned int n1p = n_p[axis1_]; + const unsigned int n1p = n_p[axis1_]; const unsigned int n1d = n_d[axis1_]; const unsigned int nyp = n_p[1]; const unsigned int nyd = n_d[1]; const unsigned int iB0 = iB_[0]; - const unsigned int p0 = iB_[0] - sign_; + const unsigned int p0 = iB_[0] - sign_; const unsigned int p1 = iB_[1] - sign_; const unsigned int iB1 = iB_[1]; - const unsigned int iB2 = iB_[2]; + const unsigned int iB2 = iB_[2]; const unsigned int p2 = iB_[2] - sign_; const int b1_size = n1p ; const int b2_size = n1d ; std::vector pos( 1 ); - if( ! vecLaser.empty() ) { + if( ! vecLaser.empty() ) { for( unsigned int j=isBoundary1min ; jgetDomainLocalMin( axis1_ ) + ( ( int )j - ( int )EMfields->oversize[axis1_] )*d[axis1_]; for( unsigned int ilaser=0; ilaser< vecLaser.size(); ilaser++ ) { diff --git a/src/ElectroMagnBC/ElectroMagnBC2D_refl.cpp b/src/ElectroMagnBC/ElectroMagnBC2D_refl.cpp index 6a6e4fdc7..a3d01be23 100755 --- a/src/ElectroMagnBC/ElectroMagnBC2D_refl.cpp +++ b/src/ElectroMagnBC/ElectroMagnBC2D_refl.cpp @@ -31,35 +31,56 @@ void ElectroMagnBC2D_refl::apply( ElectroMagn *EMfields, double, Patch *patch ) { if( i_boundary_ == 0 && patch->isXmin() ) { + const Field *B[3]{ EMfields->Bx_, EMfields->By_, EMfields->Bz_ }; + double *const __restrict__ By2D = B[1]->data_; + double *const __restrict__ Bz2D = B[2]->data_; +#ifdef SMILEI_ACCELERATOR_GPU_OACC + const int sizeofB1 = B[1]->number_of_points_; + const int sizeofB2 = B[2]->number_of_points_; +#endif + const unsigned int nyp = n_p[1]; + const unsigned int nyd = n_d[1]; // APPLICATION OF BCs OVER THE FULL GHOST CELL REGION // Static cast of the fields - //Field2D* Ex2D = static_cast(EMfields->Ex_); - //Field2D* Ey2D = static_cast(EMfields->Ey_); - //Field2D* Ez2D = static_cast(EMfields->Ez_); - //Field2D* Bx2D = static_cast(EMfields->Bx_); - Field2D *By2D = static_cast( EMfields->By_ ); - Field2D *Bz2D = static_cast( EMfields->Bz_ ); + //Field2D *By2D = static_cast( EMfields->By_ ); + //Field2D *Bz2D = static_cast( EMfields->Bz_ ); // FORCE CONSTANT MAGNETIC FIELDS - // // for Bx^(p,d) - // for (unsigned int i=oversize_; i>0; i--) { - // for (unsigned int j=0 ; j0; i-- ) { - for( unsigned int j=0 ; j0; i-- ) { - for( unsigned int j=0 ; jisXmax() ) { + const Field *B[3]{ EMfields->Bx_, EMfields->By_, EMfields->Bz_ }; + double *const __restrict__ By2D = B[1]->data_; + double *const __restrict__ Bz2D = B[2]->data_; +#ifdef SMILEI_ACCELERATOR_GPU_OACC + const int sizeofB1 = B[1]->number_of_points_; + const int sizeofB2 = B[2]->number_of_points_; +#endif + const unsigned int nxp = n_p[0]; + const unsigned int nxd = n_d[0]; + const unsigned int nyp = n_p[1]; + const unsigned int nyd = n_d[1]; // Static cast of the fields - //Field2D* Ex2D = static_cast(EMfields->Ex_); - //Field2D* Ey2D = static_cast(EMfields->Ey_); - //Field2D* Ez2D = static_cast(EMfields->Ez_); - //Field2D* Bx2D = static_cast(EMfields->Bx_); - Field2D *By2D = static_cast( EMfields->By_ ); - Field2D *Bz2D = static_cast( EMfields->Bz_ ); + //Field2D *By2D = static_cast( EMfields->By_ ); + //Field2D *Bz2D = static_cast( EMfields->Bz_ ); // FORCE CONSTANT MAGNETIC FIELDS - // // for Bx^(p,d) - // for (unsigned int i=n_p[0]-oversize_; iisYmin() ) { + const Field *B[3]{ EMfields->Bx_, EMfields->By_, EMfields->Bz_ }; + double *const __restrict__ Bx2D = B[0]->data_; + double *const __restrict__ Bz2D = B[2]->data_; +#ifdef SMILEI_ACCELERATOR_GPU_OACC + const int sizeofB0 = B[0]->number_of_points_; + const int sizeofB2 = B[2]->number_of_points_; +#endif + const unsigned int nxp = n_p[0]; + const unsigned int nxd = n_d[0]; + const unsigned int nyp = n_p[1]; + const unsigned int nyd = n_d[1]; // APPLICATION OF BCs OVER THE FULL GHOST CELL REGION // Static cast of the fields - //Field2D* Ex2D = static_cast(EMfields->Ex_); - //Field2D* Ey2D = static_cast(EMfields->Ey_); - //Field2D* Ez2D = static_cast(EMfields->Ez_); - Field2D *Bx2D = static_cast( EMfields->Bx_ ); - //Field2D* By2D = static_cast(EMfields->By_); - Field2D *Bz2D = static_cast( EMfields->Bz_ ); + //Field2D *Bx2D = static_cast( EMfields->Bx_ ); + //Field2D *Bz2D = static_cast( EMfields->Bz_ ); // FORCE CONSTANT MAGNETIC FIELDS // for Bx^(p,d) - for( unsigned int i=0; i0 ; j-- ) { - ( *Bx2D )( i, j-1 ) = ( *Bx2D )( i, j ); - }//j - }//i - - // // for By^(d,p) - // for (unsigned int i=0; i0 ; j--) { - // (*By2D)(i,j-1) = (*By2D)(i,j); - // }//j - // }//i + //( *Bx2D )( i, j-1 ) = ( *Bx2D )( i, j ); + Bx2D[i*nyd + j-1] = Bx2D[i*nyd + j]; + } + } // for Bz^(d,d) - for( unsigned int i=0; i0 ; j-- ) { - ( *Bz2D )( i, j-1 ) = ( *Bz2D )( i, j ); +#ifdef SMILEI_ACCELERATOR_GPU_OACC + #pragma acc parallel present(Bz2D[0:sizeofB2],) + #pragma acc loop gang +#elif defined( SMILEI_ACCELERATOR_GPU_OMP ) + #pragma omp target + #pragma omp teams distribute parallel for +#endif + for( unsigned int i=0; i0 ; j-- ) { + //( *Bz2D )( i, j-1 ) = ( *Bz2D )( i, j ); + Bz2D[i*nyd + j-1] = Bz2D[i*nyd + j]; } } @@ -268,34 +333,56 @@ void ElectroMagnBC2D_refl::apply( ElectroMagn *EMfields, double, Patch *patch ) } else if( i_boundary_ == 3 && patch->isYmax() ) { + const Field *B[3]{ EMfields->Bx_, EMfields->By_, EMfields->Bz_ }; + double *const __restrict__ Bx2D = B[0]->data_; + double *const __restrict__ Bz2D = B[2]->data_; +#ifdef SMILEI_ACCELERATOR_GPU_OACC + const int sizeofB0 = B[0]->number_of_points_; + const int sizeofB2 = B[2]->number_of_points_; +#endif + const unsigned int nxp = n_p[0]; + const unsigned int nxd = n_d[0]; + const unsigned int nyp = n_p[1]; + const unsigned int nyd = n_d[1]; // Static cast of the fields - //Field2D* Ex2D = static_cast(EMfields->Ex_); - //Field2D* Ey2D = static_cast(EMfields->Ey_); - //Field2D* Ez2D = static_cast(EMfields->Ez_); - Field2D *Bx2D = static_cast( EMfields->Bx_ ); - //Field2D* By2D = static_cast(EMfields->By_); - Field2D *Bz2D = static_cast( EMfields->Bz_ ); + //Field2D *Bx2D = static_cast( EMfields->Bx_ ); + //Field2D *Bz2D = static_cast( EMfields->Bz_ ); // FORCE CONSTANT MAGNETIC FIELDS // for Bx^(p,d) - for( unsigned int i=0; iBx_, EMfields->By_, EMfields->Bz_ }; + double *const __restrict__ Bx3D = B[0]->data_; + double *const __restrict__ By3D = B[1]->data_; + double *const __restrict__ Bz3D = B[2]->data_; + const unsigned int nzp = n_p[2]; + const unsigned int nzd = n_d[2]; + const unsigned int nxp = n_p[0]; + const unsigned int nxd = n_d[0]; + const unsigned int nyp = n_p[1]; + const unsigned int nyd = n_d[1]; + + const unsigned int nyz_pd = n_p[1] * n_d[2]; + const unsigned int nyz_dp = n_d[1] * n_p[2]; + const unsigned int nyz_dd = n_d[1] * n_d[2]; +#ifdef SMILEI_ACCELERATOR_GPU_OACC + const int sizeofB0 = B[0]->number_of_points_; + const int sizeofB1 = B[1]->number_of_points_; + const int sizeofB2 = B[2]->number_of_points_; +#endif if( i_boundary_ == 0 && patch->isXmin() ) { // APPLICATION OF BCs OVER THE FULL GHOST CELL REGION // Static cast of the fields - Field3D *By3D = static_cast( EMfields->By_ ); - Field3D *Bz3D = static_cast( EMfields->Bz_ ); + //Field3D *By3D = static_cast( EMfields->By_ ); + //Field3D *Bz3D = static_cast( EMfields->Bz_ ); // FORCE CONSTANT MAGNETIC FIELDS // for By^(d,p,d) +#ifdef SMILEI_ACCELERATOR_GPU_OACC + #pragma acc parallel present(By3D[0:sizeofB1]) + #pragma acc loop gang +#elif defined( SMILEI_ACCELERATOR_GPU_OMP ) + #pragma omp target + #pragma omp teams distribute parallel for +#endif for( unsigned int i=oversize_x; i>0; i-- ) { - for( unsigned int j=0 ; j0; i-- ) { - for( unsigned int j=0 ; jisXmax() ) { // Static cast of the fields - Field3D *By3D = static_cast( EMfields->By_ ); - Field3D *Bz3D = static_cast( EMfields->Bz_ ); + //Field3D *By3D = static_cast( EMfields->By_ ); + //Field3D *Bz3D = static_cast( EMfields->Bz_ ); // FORCE CONSTANT MAGNETIC FIELDS // for By^(d,p,d) - for( unsigned int i=n_d[0]-oversize_x; iisYmin() ) { // Static cast of the fields - Field3D *Bx3D = static_cast( EMfields->Bx_ ); - Field3D *Bz3D = static_cast( EMfields->Bz_ ); + //Field3D *Bx3D = static_cast( EMfields->Bx_ ); + //Field3D *Bz3D = static_cast( EMfields->Bz_ ); // FORCE CONSTANT MAGNETIC FIELDS // for Bx^(p,d,d) - for( unsigned int i=0; i0 ; j-- ) { - for( unsigned int k=0; k0 ; j-- ) { - for( unsigned int k=0; kisYmax() ) { // Static cast of the fields - Field3D *Bx3D = static_cast( EMfields->Bx_ ); - Field3D *Bz3D = static_cast( EMfields->Bz_ ); + //Field3D *Bx3D = static_cast( EMfields->Bx_ ); + //Field3D *Bz3D = static_cast( EMfields->Bz_ ); // FORCE CONSTANT MAGNETIC FIELDS // for Bx^(p,d,d) - for( unsigned int i=0; iisZmin() ) { // Static cast of the fields - Field3D *Bx3D = static_cast( EMfields->Bx_ ); - Field3D *By3D = static_cast( EMfields->By_ ); + //Field3D *Bx3D = static_cast( EMfields->Bx_ ); + //Field3D *By3D = static_cast( EMfields->By_ ); // FORCE CONSTANT MAGNETIC FIELDS // for Bx^(p,d,d) - for( unsigned int i=0; i0 ; k-- ) { - ( *Bx3D )( i, j, k-1 ) = ( *Bx3D )( i, j, k ); + //( *Bx3D )( i, j, k-1 ) = ( *Bx3D )( i, j, k ); + Bx3D[i*nyz_dd + j*nzd + k-1] = Bx3D[i*nyz_dd + j*nzd + k]; } } } // for By^(d,p,d) - for( unsigned int i=0; i0 ; k-- ) { - ( *By3D )( i, j, k-1 ) = ( *By3D )( i, j, k ); + //( *By3D )( i, j, k-1 ) = ( *By3D )( i, j, k ); + By3D[i*nyz_pd + j*nzd + k-1] = By3D[i*nyz_pd + j*nzd + k]; } } } @@ -169,25 +295,47 @@ void ElectroMagnBC3D_refl::apply( ElectroMagn *EMfields, double, Patch *patch ) } else if( i_boundary_==5 && patch->isZmax() ) { // Static cast of the fields - Field3D *Bx3D = static_cast( EMfields->Bx_ ); - Field3D *By3D = static_cast( EMfields->By_ ); + //Field3D *Bx3D = static_cast( EMfields->Bx_ ); + //Field3D *By3D = static_cast( EMfields->By_ ); // FORCE CONSTANT MAGNETIC FIELDS // for Bx^(p,d,d) - for( unsigned int i=0; i namelistsFiles ) : string seterr( "seterr" ); string sChar( "s" ); Py_DECREF( PyObject_CallMethod( numpy, &seterr[0], &sChar[0], "ignore" ) ); + string numpy_version = ""; + PyTools::getAttr( numpy, "__version__", numpy_version ); + MESSAGE( "Numpy version " << numpy_version ); Py_DECREF( numpy ); #else WARNING("Numpy not found. Some options will not be available"); @@ -868,7 +871,8 @@ Params::Params( SmileiMPI *smpi, std::vector namelistsFiles ) : if( geometry!="1Dcartesian" && geometry!="2Dcartesian" && geometry!="3Dcartesian" ) { - ERROR_NAMELIST( "Collisions only valid for cartesian geometries for the moment", LINK_NAMELIST + std::string("#collisions-reactions") ); + //ERROR_NAMELIST( "Collisions only valid for cartesian geometries for the moment", LINK_NAMELIST + std::string("#collisions-reactions") ); + WARNING( "Collisions in AM geometry is experimental and valid only with a single mode" ); } } diff --git a/src/Patch/PatchAM.h b/src/Patch/PatchAM.h index 9d8f2c17c..df58a8016 100755 --- a/src/Patch/PatchAM.h +++ b/src/Patch/PatchAM.h @@ -24,10 +24,24 @@ class PatchAM final : public Patch //! Return the volume (or surface or length depending on simulation dimension) //! of one cell at the position of a given particle - double getPrimalCellVolume( Particles *, unsigned int, Params & ) override final + double getPrimalCellVolume( Particles *p, unsigned int ipart, Params ¶ms ) override final { - ERROR( "getPrimalCellVolume not implemented in geometry AM" ); - return cell_volume; + double factor = 1.; + + double halfcell = 0.5 * params.cell_length[0]; + if( p->position(0,ipart) - getDomainLocalMin(0) < halfcell + || getDomainLocalMax(0) - p->position(0,ipart) < halfcell ) { + factor *= 0.5; + } + + halfcell = 0.5 * params.cell_length[1]; + double radius = sqrt(p->position(1,ipart)*p->position(1,ipart) + p->position(2,ipart)*p->position(2,ipart)); + if( radius - getDomainLocalMin(1) < halfcell + || getDomainLocalMax(1) - radius < halfcell ) { + factor *= 0.5; + } + + return factor * cell_volume * radius; }; //! Given several arrays (x,y,z), return indices of points in patch diff --git a/src/Patch/VectorPatch.cpp b/src/Patch/VectorPatch.cpp index bbce82778..b3d2190a2 100755 --- a/src/Patch/VectorPatch.cpp +++ b/src/Patch/VectorPatch.cpp @@ -540,7 +540,7 @@ void VectorPatch::injectParticlesFromBoundaries(Params ¶ms, Timers &timers, timers.particleInjection.restart(); //#pragma omp for schedule(runtime) - #pragma omp single + #pragma omp master for( unsigned int ipatch=0 ; ipatchsize() ; ipatch++ ) { Patch * patch = ( *this )( ipatch ); @@ -779,6 +779,7 @@ void VectorPatch::injectParticlesFromBoundaries(Params ¶ms, Timers &timers, } } // end for ipatch + #pragma omp barrier timers.particleInjection.update( params.printNow( itime ) ); } @@ -1169,6 +1170,7 @@ void VectorPatch::finalizeSyncAndBCFields( Params ¶ms, SmileiMPI *smpi, SimW } timers.maxwellBC.restart(); + SMILEI_PY_SAVE_MASTER_THREAD #pragma omp for schedule(static) for( unsigned int ipatch=0 ; ipatchsize() ; ipatch++ ) { // Applies boundary conditions on B @@ -1176,6 +1178,7 @@ void VectorPatch::finalizeSyncAndBCFields( Params ¶ms, SmileiMPI *smpi, SimW ( *this )( ipatch )->EMfields->boundaryConditions( time_dual, ( *this )( ipatch ), simWindow ); } + SMILEI_PY_RESTORE_MASTER_THREAD SyncVectorPatch::exchangeForPML( params, (*this), smpi ); #pragma omp for schedule(static) @@ -1284,9 +1287,6 @@ void VectorPatch::runAllDiags( Params &/*params*/, SmileiMPI *smpi, unsigned int for( unsigned int idiag = 0 ; idiag < globalDiags.size() ; idiag++ ) { if( globalDiags[idiag]->timeSelection->theTimeIsNow( itime ) ) { - //std::cout << " " << dynamic_cast( globalDiags[idiag] ) - // << std::endl; - if (dynamic_cast( globalDiags[idiag])) { //need_particles = true; //need_fields = true; @@ -1314,7 +1314,7 @@ void VectorPatch::runAllDiags( Params &/*params*/, SmileiMPI *smpi, unsigned int } else if (dynamic_cast(localDiags[idiag])) { need_fields = true; } else if (dynamic_cast(localDiags[idiag])) { - // Nothing to be done + // Nothing to be done } else { need_particles = true; need_fields = true; @@ -1508,12 +1508,6 @@ void VectorPatch::runAllDiags( Params &/*params*/, SmileiMPI *smpi, unsigned int } timers.diags.update(); - if( itime==0 ) { - for( unsigned int idiag = 0 ; idiag < diag_timers_.size() ; idiag++ ) { - diag_timers_[idiag]->reboot(); - } - } - } // END runAllDiags // --------------------------------------------------------------------------------------------------------------------- @@ -1636,13 +1630,14 @@ void VectorPatch::runAllDiagsTasks( Params &, SmileiMPI *smpi, unsigned int itim } timers.diags.update(); - if (itime==0) { - for( unsigned int idiag = 0 ; idiag < diag_timers_.size() ; idiag++ ) - diag_timers_[idiag]->reboot(); - } - } // END runAllDiags +void VectorPatch::rebootDiagTimers() +{ + for( unsigned int idiag = 0 ; idiag < diag_timers_.size() ; idiag++ ) { + diag_timers_[idiag]->reboot(); + } +} // --------------------------------------------------------------------------------------------------------------------- // Check if rho is null (MPI & patch sync) @@ -4033,9 +4028,10 @@ void VectorPatch::applyAntennas( double time ) } else { // Get intensity from antenna of the first patch - #pragma omp single + #pragma omp master antenna_intensity_ = patches_[0]->EMfields->antennas[iAntenna].time_profile->valueAt( time ); - + #pragma omp barrier + // Loop patches to apply #pragma omp for schedule(static) for( unsigned int ipatch=0 ; ipatchdiagPartEventTracing( time_dual, params.timestep); #endif + SMILEI_PY_SAVE_MASTER_THREAD #pragma omp for schedule(runtime) for( unsigned int ipatch=0 ; ipatchsize() ; ipatch++ ) { ( *this )( ipatch )->EMfields->restartRhoJ(); @@ -5017,6 +5014,7 @@ void VectorPatch::dynamicsWithoutTasks( Params ¶ms, } // end loop on species //MESSAGE("species dynamics"); } // end loop on patches + SMILEI_PY_RESTORE_MASTER_THREAD } void VectorPatch::ponderomotiveUpdateSusceptibilityAndMomentumWithoutTasks( Params ¶ms, diff --git a/src/Patch/VectorPatch.h b/src/Patch/VectorPatch.h index d2e12e545..dab9d0ce2 100755 --- a/src/Patch/VectorPatch.h +++ b/src/Patch/VectorPatch.h @@ -238,6 +238,7 @@ public : SimWindow *simWindow ); void runAllDiagsTasks( Params ¶ms, SmileiMPI *smpi, unsigned int itime, Timers &timers, SimWindow *simWindow ); + void rebootDiagTimers(); void initAllDiags( Params ¶ms, SmileiMPI *smpi ); void closeAllDiags( SmileiMPI *smpi ); diff --git a/src/Profiles/Function.cpp b/src/Profiles/Function.cpp index 6a0cf2ac0..d39994722 100755 --- a/src/Profiles/Function.cpp +++ b/src/Profiles/Function.cpp @@ -8,81 +8,144 @@ using namespace std; // 1D double Function_Python1D::valueAt( double time ) { - return PyTools::runPyFunction( py_profile, time ); + double v; + SMILEI_PY_ACQUIRE_GIL + v = PyTools::runPyFunction( py_profile, time ); + SMILEI_PY_RELEASE_GIL + return v; } double Function_Python1D::valueAt( vector, double time ) { - return PyTools::runPyFunction( py_profile, time ); + double v; + SMILEI_PY_ACQUIRE_GIL + v = PyTools::runPyFunction( py_profile, time ); + SMILEI_PY_RELEASE_GIL + return v; } double Function_Python1D::valueAt( vector x_cell ) { - return PyTools::runPyFunction( py_profile, x_cell[0] ); + double v; + SMILEI_PY_ACQUIRE_GIL + v = PyTools::runPyFunction( py_profile, x_cell[0] ); + SMILEI_PY_RELEASE_GIL + return v; } // 2D double Function_Python2D::valueAt( vector x_cell, double time ) { - return PyTools::runPyFunction( py_profile, x_cell[0], time ); + double v; + SMILEI_PY_ACQUIRE_GIL + v = PyTools::runPyFunction( py_profile, x_cell[0], time ); + SMILEI_PY_RELEASE_GIL + return v; } double Function_Python2D::valueAt( vector x_cell ) { - return PyTools::runPyFunction( py_profile, x_cell[0], x_cell[1] ); + double v; + SMILEI_PY_ACQUIRE_GIL + v = PyTools::runPyFunction( py_profile, x_cell[0], x_cell[1] ); + SMILEI_PY_RELEASE_GIL + return v; } // 2D complex std::complex Function_Python2D::complexValueAt( vector x_cell, double time ) { - return PyTools::runPyFunction>( py_profile, x_cell[0], time ); + std::complex v; + SMILEI_PY_ACQUIRE_GIL + v = PyTools::runPyFunction>( py_profile, x_cell[0], time ); + SMILEI_PY_RELEASE_GIL + return v; } std::complex Function_Python2D::complexValueAt( vector x_cell ) { - return PyTools::runPyFunction>( py_profile, x_cell[0], x_cell[1] ); + std::complex v; + SMILEI_PY_ACQUIRE_GIL + v = PyTools::runPyFunction>( py_profile, x_cell[0], x_cell[1] ); + SMILEI_PY_RELEASE_GIL + return v; } // 3D double Function_Python3D::valueAt( vector x_cell, double time ) { - return PyTools::runPyFunction( py_profile, x_cell[0], x_cell[1], time ); + double v; + SMILEI_PY_ACQUIRE_GIL + v = PyTools::runPyFunction( py_profile, x_cell[0], x_cell[1], time ); + SMILEI_PY_RELEASE_GIL + return v; } double Function_Python3D::valueAt( vector x_cell ) { - return PyTools::runPyFunction( py_profile, x_cell[0], x_cell[1], x_cell[2] ); + double v; + SMILEI_PY_ACQUIRE_GIL + v = PyTools::runPyFunction( py_profile, x_cell[0], x_cell[1], x_cell[2] ); + SMILEI_PY_RELEASE_GIL + return v; } // 3D complex std::complex Function_Python3D::complexValueAt( vector x_cell, double time ) { - return PyTools::runPyFunction>( py_profile, x_cell[0], x_cell[1], time ); + std::complex v; + SMILEI_PY_ACQUIRE_GIL + v = PyTools::runPyFunction>( py_profile, x_cell[0], x_cell[1], time ); + SMILEI_PY_RELEASE_GIL + return v; } // 4D double Function_Python4D::valueAt( vector x_cell, double time ) { - return PyTools::runPyFunction( py_profile, x_cell[0], x_cell[1], x_cell[2], time ); + double v; + SMILEI_PY_ACQUIRE_GIL + v = PyTools::runPyFunction( py_profile, x_cell[0], x_cell[1], x_cell[2], time ); + SMILEI_PY_RELEASE_GIL + return v; } // 4D complex std::complex Function_Python4D::complexValueAt( vector x_cell, double time ) { - return PyTools::runPyFunction>( py_profile, x_cell[0], x_cell[1], x_cell[2], time ); + std::complex v; + SMILEI_PY_ACQUIRE_GIL + v = PyTools::runPyFunction>( py_profile, x_cell[0], x_cell[1], x_cell[2], time ); + SMILEI_PY_RELEASE_GIL + return v; } // Special cases for locations specified in numpy arrays #ifdef SMILEI_USE_NUMPY PyArrayObject *Function_Python1D::valueAt( std::vector x ) { - return ( PyArrayObject * )PyObject_CallFunctionObjArgs( py_profile, x[0], NULL ); + PyArrayObject * v; + SMILEI_PY_ACQUIRE_GIL + v = ( PyArrayObject * )PyObject_CallFunctionObjArgs( py_profile, x[0], NULL ); + SMILEI_PY_RELEASE_GIL + return v; } PyArrayObject *Function_Python2D::valueAt( std::vector x ) { - return ( PyArrayObject * )PyObject_CallFunctionObjArgs( py_profile, x[0], x[1], NULL ); + PyArrayObject * v; + SMILEI_PY_ACQUIRE_GIL + v = ( PyArrayObject * )PyObject_CallFunctionObjArgs( py_profile, x[0], x[1], NULL ); + SMILEI_PY_RELEASE_GIL + return v; } PyArrayObject *Function_Python3D::valueAt( std::vector x ) { - return ( PyArrayObject * )PyObject_CallFunctionObjArgs( py_profile, x[0], x[1], x[2], NULL ); + PyArrayObject * v; + SMILEI_PY_ACQUIRE_GIL + v = ( PyArrayObject * )PyObject_CallFunctionObjArgs( py_profile, x[0], x[1], x[2], NULL ); + SMILEI_PY_RELEASE_GIL + return v; } PyArrayObject *Function_Python2D::complexValueAt( std::vector x ) { + PyArrayObject *cvalues; + SMILEI_PY_ACQUIRE_GIL PyObject *values = PyObject_CallFunctionObjArgs( py_profile, x[0], x[1], NULL ); - PyArrayObject *cvalues = ( PyArrayObject * )PyObject_CallMethod( values, const_cast("astype"), const_cast("s"), const_cast("complex"), NULL ); + cvalues = ( PyArrayObject * )PyObject_CallMethod( values, const_cast("astype"), const_cast("s"), const_cast("complex"), NULL ); Py_DECREF( values ); + SMILEI_PY_RELEASE_GIL return cvalues; } @@ -90,45 +153,63 @@ PyArrayObject *Function_Python2D::complexValueAt( std::vector x PyArrayObject *Function_Python1D::valueAt( std::vector, double time ) { + PyArrayObject * ret; + SMILEI_PY_ACQUIRE_GIL PyObject *t = PyFloat_FromDouble( time ); - PyArrayObject * ret = ( PyArrayObject * )PyObject_CallFunctionObjArgs( py_profile, t, NULL ); + ret = ( PyArrayObject * )PyObject_CallFunctionObjArgs( py_profile, t, NULL ); Py_DECREF( t ); + SMILEI_PY_RELEASE_GIL return ret; } PyArrayObject *Function_Python2D::valueAt( std::vector x, double time ) { + PyArrayObject * ret; + SMILEI_PY_ACQUIRE_GIL PyObject *t = PyFloat_FromDouble( time ); - PyArrayObject * ret = ( PyArrayObject * )PyObject_CallFunctionObjArgs( py_profile, x[0], t, NULL ); + ret = ( PyArrayObject * )PyObject_CallFunctionObjArgs( py_profile, x[0], t, NULL ); Py_DECREF( t ); + SMILEI_PY_RELEASE_GIL return ret; } PyArrayObject *Function_Python3D::valueAt( std::vector x, double time ) { + PyArrayObject * ret; + SMILEI_PY_ACQUIRE_GIL PyObject *t = PyFloat_FromDouble( time ); - PyArrayObject * ret = ( PyArrayObject * )PyObject_CallFunctionObjArgs( py_profile, x[0], x[1], t, NULL ); + ret = ( PyArrayObject * )PyObject_CallFunctionObjArgs( py_profile, x[0], x[1], t, NULL ); Py_DECREF( t ); + SMILEI_PY_RELEASE_GIL return ret; } PyArrayObject *Function_Python4D::valueAt( std::vector x, double time ) { + PyArrayObject * ret; + SMILEI_PY_ACQUIRE_GIL PyObject *t = PyFloat_FromDouble( time ); - PyArrayObject * ret = ( PyArrayObject * )PyObject_CallFunctionObjArgs( py_profile, x[0], x[1], x[2], t, NULL ); + ret = ( PyArrayObject * )PyObject_CallFunctionObjArgs( py_profile, x[0], x[1], x[2], t, NULL ); Py_DECREF( t ); + SMILEI_PY_RELEASE_GIL return ret; }PyArrayObject *Function_Python4D::complexValueAt( std::vector x, double time ) { + PyArrayObject *cvalues; + SMILEI_PY_ACQUIRE_GIL PyObject *t = PyFloat_FromDouble( time ); PyObject * values = PyObject_CallFunctionObjArgs( py_profile, x[0], x[1], x[2], t, NULL ); Py_DECREF( t ); - PyArrayObject *cvalues = ( PyArrayObject * )PyObject_CallMethod( values, const_cast("astype"), const_cast("s"), const_cast("complex"), NULL ); + cvalues = ( PyArrayObject * )PyObject_CallMethod( values, const_cast("astype"), const_cast("s"), const_cast("complex"), NULL ); Py_DECREF( values ); + SMILEI_PY_RELEASE_GIL return cvalues; } PyArrayObject *Function_Python4D::complexValueAt( std::vector x, PyArrayObject *t ) { + PyArrayObject *cvalues; + SMILEI_PY_ACQUIRE_GIL PyObject *values = PyObject_CallFunctionObjArgs( py_profile, x[0], x[1], x[2], t, NULL ); - PyArrayObject *cvalues = ( PyArrayObject * )PyObject_CallMethod( values, const_cast("astype"), const_cast("s"), const_cast("complex"), NULL ); + cvalues = ( PyArrayObject * )PyObject_CallMethod( values, const_cast("astype"), const_cast("s"), const_cast("complex"), NULL ); Py_DECREF( values ); + SMILEI_PY_RELEASE_GIL return cvalues; } diff --git a/src/Smilei.cpp b/src/Smilei.cpp index 11e43976e..1c30dfddc 100755 --- a/src/Smilei.cpp +++ b/src/Smilei.cpp @@ -423,12 +423,15 @@ int main( int argc, char *argv[] ) if( !params.restart ) { TITLE( "Running diags at time t = 0" ); + #pragma omp parallel shared( smpi, params, vecPatches, simWindow ) + { #ifdef _OMPTASKS - vecPatches.runAllDiagsTasks( params, &smpi, 0, timers, simWindow ); + vecPatches.runAllDiagsTasks( params, &smpi, 0, timers, simWindow ); #else - vecPatches.runAllDiags( params, &smpi, 0, timers, simWindow ); + vecPatches.runAllDiags( params, &smpi, 0, timers, simWindow ); #endif - + } + vecPatches.rebootDiagTimers(); } TITLE( "Species creation summary" ); @@ -650,7 +653,7 @@ int main( int argc, char *argv[] ) // Standard fields operations (maxwell + comms + boundary conditions) are completed // apply prescribed fields can be considered if requested if( vecPatches(0)->EMfields->prescribedFields.size() ) { - #pragma omp single + #pragma omp master vecPatches.applyPrescribedFields( time_prim ); #pragma omp barrier } diff --git a/src/Species/SpeciesFactory.h b/src/Species/SpeciesFactory.h index edef5e334..0dcb39fcf 100755 --- a/src/Species/SpeciesFactory.h +++ b/src/Species/SpeciesFactory.h @@ -671,7 +671,7 @@ class SpeciesFactory } if( (params.geometry=="AMcylindrical") && ( this_species->boundary_conditions_[1][1] != "remove" ) && ( this_species->boundary_conditions_[1][1] != "stop" ) && ( this_species->boundary_conditions_[1][1] != "reflective" ) ) { ERROR_NAMELIST( - " In AM geometry particle boundary conditions supported in Rmax are 'remove' and 'stop' ", + " In AM geometry particle boundary conditions supported in Rmax are 'remove', 'reflective' and 'stop' ", LINK_NAMELIST + std::string("#species") ); } diff --git a/src/Species/SpeciesVAdaptive.cpp b/src/Species/SpeciesVAdaptive.cpp index 273362561..4152a5e7c 100755 --- a/src/Species/SpeciesVAdaptive.cpp +++ b/src/Species/SpeciesVAdaptive.cpp @@ -346,18 +346,33 @@ void SpeciesVAdaptive::scalarDynamics( double time_dual, unsigned int ispec, if (time_dual <= time_frozen_ && diag_flag &&( !particles->is_test ) ) { //immobile particle (at the moment only project density) + if( params.geometry != "AMcylindrical" ) { - smpi->traceEventIfDiagTracing(diag_PartEventTracing, ithread,0,3); - double *b_rho=nullptr; - for( unsigned int ibin = 0 ; ibin < particles->first_index.size() ; ibin ++ ) { //Loop for projection on buffer_proj - - b_rho = EMfields->rho_s[ispec] ? &( *EMfields->rho_s[ispec] )( 0 ) : &( *EMfields->rho_ )( 0 ) ; - for( iPart=particles->first_index[ibin] ; ( int )iPartlast_index[ibin]; iPart++ ) { - Proj->basic( b_rho, ( *particles ), iPart, 0 ); + smpi->traceEventIfDiagTracing(diag_PartEventTracing, ithread,0,3); + double *b_rho = EMfields->rho_s[ispec] ? &( *EMfields->rho_s[ispec] )( 0 ) : &( *EMfields->rho_ )( 0 ) ; + for( unsigned int scell = 0 ; scell < particles->first_index.size() ; scell ++ ) { //Loop for projection on buffer_proj + for( iPart=particles->first_index[scell] ; ( int )iPartlast_index[scell]; iPart++ ) { + Proj->basic( b_rho, ( *particles ), iPart, 0 ); + } } + smpi->traceEventIfDiagTracing(diag_PartEventTracing, ithread,1,3); + + } else { // AM case + ElectroMagnAM *emAM = static_cast( EMfields ); + int n_species = patch->vecSpecies.size(); + + smpi->traceEventIfDiagTracing(diag_PartEventTracing, ithread,0,3); + for( unsigned int imode = 0; imode *b_rho = emAM->rho_AM_s[ifield] ? &( *emAM->rho_AM_s[ifield] )( 0 ) : &( *emAM->rho_AM_[imode] )( 0 ) ; + for( unsigned int scell = 0 ; scell < particles->first_index.size() ; scell ++ ) { //Loop for projection on buffer_proj + for( int iPart=particles->first_index[scell] ; iPartlast_index[scell]; iPart++ ) { + Proj->basicForComplex( b_rho, ( *particles ), iPart, 0, imode ); + } + } + } + smpi->traceEventIfDiagTracing(diag_PartEventTracing, ithread,1,3); } - smpi->traceEventIfDiagTracing(diag_PartEventTracing, ithread,1,3); - } }//END scalarDynamics diff --git a/src/Tools/PyTools.h b/src/Tools/PyTools.h index a756f3db8..0a5ee96f4 100755 --- a/src/Tools/PyTools.h +++ b/src/Tools/PyTools.h @@ -44,8 +44,8 @@ // } // SMILEI_PY_RESTORE_MASTER_THREAD #if PY_MAJOR_VERSION > 2 && PY_MINOR_VERSION > 11 - #define SMILEI_PY_SAVE_MASTER_THREAD PyThreadState *_save = NULL; _Pragma("omp master") if( Py_IsInitialized() ) _save = PyEval_SaveThread(); - #define SMILEI_PY_RESTORE_MASTER_THREAD _Pragma("omp master") if( Py_IsInitialized() ) PyEval_RestoreThread( _save ); + #define SMILEI_PY_SAVE_MASTER_THREAD PyThreadState *_save = NULL; _Pragma("omp master") if( Py_IsInitialized() ) _save = PyEval_SaveThread(); _Pragma("omp barrier") + #define SMILEI_PY_RESTORE_MASTER_THREAD _Pragma("omp master") if( Py_IsInitialized() ) PyEval_RestoreThread( _save ); _Pragma("omp barrier") #define SMILEI_PY_ACQUIRE_GIL PyGILState_STATE _state = PyGILState_Ensure(); #define SMILEI_PY_RELEASE_GIL PyGILState_Release( _state ); #else diff --git a/validation/analyses/gpu_validate_tst2d_01_plasma_mirror.py b/validation/analyses/gpu_validate_tst2d_01_plasma_mirror.py new file mode 100755 index 000000000..e958d40ee --- /dev/null +++ b/validation/analyses/gpu_validate_tst2d_01_plasma_mirror.py @@ -0,0 +1,17 @@ +import os, re, numpy as np, math +from scipy.ndimage import gaussian_filter as gfilt +import happi + +#S = happi.Open(["./restart*"], verbose=False) + +S = happi.Open(verbose=False) + +# COMPARE THE TOTAL NUMBER OF CREATED IONS +nion = S.Scalar.Ntot_ion(timesteps=0).getData()[0] +Validate("Number of ions at iteration 0", nion) + +# COMPARE THE Ey FIELD +Ey = S.Field.Field0.Ey(timesteps=800).getData()[0] +Ey = gfilt(Ey, 6) # smoothing +Ey = Ey[50:200:4, 300::4] # zoom on the reflected laser +Validate("Ey field at iteration 800", Ey, 0.02) diff --git a/validation/analyses/validate_tst1d_05_tunnel_ionisation.py b/validation/analyses/validate_tst1d_05_tunnel_ionisation.py index 19148bd74..fe26ad1d5 100755 --- a/validation/analyses/validate_tst1d_05_tunnel_ionisation.py +++ b/validation/analyses/validate_tst1d_05_tunnel_ionisation.py @@ -93,7 +93,7 @@ def calculate_ionization(Ip, l): d = S.TrackParticles("electron", axes=["Id","x","Wy"], timesteps=1000).getData() keep = d["Id"] > 0 order = np.argsort(d["x"][keep]) -Validate("Track electron x", d["x"][keep][order][::200], 1e-4) +Validate("Track electron x", d["x"][keep][order][::200], 2e-4) Validate("Track electron Wy", gaussian_filter(maximum_filter1d(d["Wy"][keep][order],20),200)[::200], 1e-5) # NEW PARTICLES DIAGNOSTIC diff --git a/validation/analyses/validate_tst_collisions3_AM_uniformity.py b/validation/analyses/validate_tst_collisions3_AM_uniformity.py new file mode 100644 index 000000000..edb52201a --- /dev/null +++ b/validation/analyses/validate_tst_collisions3_AM_uniformity.py @@ -0,0 +1,7 @@ +import happi + +S = happi.Open(["./restart*"], verbose=False) + +Validate("ion Pyy vs x", S.ParticleBinning(0, sum={"y":"all"}).getData()[-1], 0.4e-10) +Validate("ion Pyy vs y", S.ParticleBinning(0, sum={"x":"all"}).getData()[-1], 0.4e-10) + diff --git a/validation/references/tst_collisions3_AM_uniformity.py.txt b/validation/references/tst_collisions3_AM_uniformity.py.txt new file mode 100644 index 000000000..146b36250 Binary files /dev/null and b/validation/references/tst_collisions3_AM_uniformity.py.txt differ