diff --git a/.github/ISSUE_TEMPLATE/general-question.md b/.github/ISSUE_TEMPLATE/general-question.md deleted file mode 100644 index 6519650b9..000000000 --- a/.github/ISSUE_TEMPLATE/general-question.md +++ /dev/null @@ -1,11 +0,0 @@ ---- -name: General question -about: Is there something unclear, undocumented, or that you cannot figure how to - achieve? -title: '' -labels: '' -assignees: '' - ---- - - diff --git a/.github/ISSUE_TEMPLATE/installation-errors.md b/.github/ISSUE_TEMPLATE/installation-errors.md deleted file mode 100644 index e724a9a84..000000000 --- a/.github/ISSUE_TEMPLATE/installation-errors.md +++ /dev/null @@ -1,24 +0,0 @@ ---- -name: Installation errors -about: You did not manage to install Smilei or have it run normally -title: '' -labels: installation -assignees: '' - ---- - - - * Need help regarding installation? - Read this doc first: https://smileipic.github.io/Smilei/installation.html - - * Issues installing C++, MPI, HDF5 or python? - Please refer to your system administrator first. - - * If you are still in trouble, please provide the result of these commands in the issue: - git describe --all --long - make config=verbose - make env - echo $LD_LIBRARY_PATH # on mac: echo $DYLD_LIBRARY_PATH - python -m sysconfig - ldd smilei - diff --git a/doc/Sphinx/Overview/material.rst b/doc/Sphinx/Overview/material.rst index aa7c9f06f..45d06b92b 100644 --- a/doc/Sphinx/Overview/material.rst +++ b/doc/Sphinx/Overview/material.rst @@ -48,6 +48,84 @@ As of November 2021, 90 papers have been published covering a broad range of top Use the python script doc/doi2publications.py to generate entries from a DOI number, and paste them here +.. [Gorlova2024] + + D. A. Gorlova, I. N. Tsymbalov, I. P. Tsygvintsev and A. B. Savelev, + `THz transition radiation of electron bunches laser-accelerated in long-scale near-critical-density plasmas`, + `Laser Physics Letters 21, 035001 (2024) `_ + +.. [Seidel2024] + + A. Seidel, B. Lei, C. Zepter, M. C. Kaluza, A. Sävert, M. Zepf, and D. Seipt, + `Polarization and CEP dependence of the transverse phase space in laser driven accelerators`, + `Physical Review Research 6, 013056 (2024) `_ + +.. [Gao2023b] + + X. Gao, + `Anisotropic field ionization in nanoclusters mediated by a Brunel-electron-driven plasma wave`, + `Physical Review A 108, 033109 (2023) `_ + +.. [Yoon2023b] + + Y. D. Yoon, P. M. Bellan and G. S. Yun, + `Phase-space Analysis of Ordered and Disordered Nonthermal Ion Energization during Magnetic Reconnection`, + `The Astrophysical Journal, 956:105 (2023) `_ + +.. [Bhadoria2023] + + S. Bhadoria, M. Marklund and C. H. Keitel, + `Energy enhancement of laser-driven ions by radiation reaction and Breit-Wheeler pair production in the ultra-relativistic transparency regime`, + `High Power Laser Science and Engineering (2023) `_ + +.. [Diab2023] + + R. Diab, S.-G. Baek, P. Bonoli, T. G. Jenkins, M. Ono and D. Smithe, + `Particle-in-cell simulations of parasitic electrostatic wave excitation in the ion cyclotron range of frequencies and high harmonic fast wave regimes`, + `AIP Conference Proceedings 2984, 080001 (2023) `_ + +.. [Sladkov2023] + + A. D. Sladkov and A. V. Korzhimanov, + `Cherenkov Radiation of an Ultrashort Laser Pulse Propagating in a Strongly Magnetized Plasma at Various Intensities and Directions of the Magnetic Field`, + `Radiophysics and Quantum Electronics 65, 888–896 (2023) `_ + +.. [Montefiori2023] + + S. Montefiori and M. Tamburini + `SFQEDtoolkit: A high-performance library for the accurate modeling of strong-field QED processes in PIC and Monte Carlo codes`, + `Computer Physics Communications 292, 108855 (2023) `_ + +.. [Shekhanov2023] + + S. Shekhanov, A. Gintrand, L. Hudec, R. Liska, J. Limpouch, S. Weber and V. Tikhonchuk + `Kinetic modeling of laser absorption in foams`, + `Physics of Plasmas 30, 012708 (2023) `_ + +.. [Yu2023] + + J. Yu, J. Zhong, Y. Ping and W. An + `Electron acceleration in a coil target-driven low-β magnetic reconnection simulation`, + `Matter and Radiation at Extremes 8, 064003 (2023) `_ + +.. [Zagidullin2023] + + R. Zagidullin, S. Tietze, M. Zepf, J. Wang and S. Rykovanov + `Density-dependent carrier-envelope phase shift in attosecond pulse generation from relativistically oscillating mirrors`, + `Matter and Radiation at Extremes 8, 064004 (2023) `_ + +.. [Cai2023] + + J. Cai, Y. Shou, Y. Geng, L. Han, X. Xu, S. Wen, B. Shen, J. Yu and X. Yan + `Extremely powerful and frequency-tunable terahertz pulses from a table-top laser-plasma wiggler`, + `High Power Laser Science and Engineering (2023) `_ + +.. [Jirka2023] + + M. Jirka and H. Kladecová + `Pair production in an electron collision with a radially polarized laser pulse`, + `Physics of Plasmas 30, 113102 (2023) `_ + .. [Guo2023] A. Guo, Q. Lu, S. Lu, S. Wang and R. Wang, @@ -128,7 +206,7 @@ As of November 2021, 90 papers have been published covering a broad range of top .. [Gao2023] - Xiaohui Gao, + X. Gao, `Ionization dynamics of sub-micrometer-sized clusters in intense ultrafast laser pulses`, `Phys. Plasmas 30, 052102 (2023) `_ @@ -146,13 +224,13 @@ As of November 2021, 90 papers have been published covering a broad range of top .. [Ghizzo2023] - Alain Ghizzo, Daniele Del Sarto, and Homam Betar, + A. Ghizzo, D. Del Sarto, and H. Betar, `Collisionless Heating Driven by Vlasov Filamentation in a Counterstreaming Beams Configuration`, `Phys. Rev. Lett. 131, 035101 (2023) `_ .. [Yang2023] - Tong Yang, Zhen Guo, Yang Yan, Minjian Wu, Yadong Xia, Qiangyou He, Hao Cheng, Yuze Li, Yanlv Fang, Yanying Zhao, Xueqing Yan and Chen Lin, + T. Yang, Z. Guo, Y. Yan, M. Wu, Y. Xia, Q. He, H. Cheng, Y. Li, Y. Fang, Y. Zhao, X. Yan and C. Lin, `Measurements of Plasma Density Profile Evolutions with Channel-guided Laser`, `High Power Laser Science and Engineering pp. 1-15 (2023) `_ @@ -164,25 +242,25 @@ As of November 2021, 90 papers have been published covering a broad range of top .. [Pak2023] - Taegyu Pak, Mohammad Rezaei-Pandari, Sang Beom Kim, Geonwoo Lee, Dae Hee Wi, Calin Ioan Hojbota, Mohammad Mirzaie, Hyeongmun Kim, Jae Hee Sung, Seong Ku Lee, Chul Kang and Ki-Yong Kim, + T. Pak, M. Rezaei-Pandari, S. B. Kim, G. Lee, D. H. Wi, C. I. Hojbota, M. Mirzaie, H. Kim, J. H. Sung, S. K. Lee, C. Kang and K.-Y. Kim, `Multi-millijoule terahertz emission from laser-wakefield-accelerated electrons`, `Light Sci Appl 12, 37 (2023) `_ .. [Istokskaia2023] - Valeriia Istokskaia, Marco Tosca, Lorenzo Giuffrida, Jan Psikal, Filip Grepl, Vasiliki Kantarelou, Stanislav Stancek, Sabrina Di Siena, Arsenios Hadjikyriacou, Aodhan McIlvenny, Yoann Levy, Jaroslav Huynh, Martin Cimrman, Pavel Pleskunov, Daniil Nikitin, Andrei Choukourov, Fabio Belloni, Antonino Picciotto, Satyabrata Kar, Marco Borghesi, Antonio Lucianetti, Tomas Mocek and Daniele Margarone, + 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`, `Commun Phys 6, 27 (2023) `_ .. [Yoon2023] - Young Dae Yoon, Deirdre E. Wendel and Gunsu S. Yun, + Y. D. Yoon, D. E. Wendel and G. S. Yun, `Equilibrium selection via current sheet relaxation and guide field amplification`, `Nat Commun 14, 139 (2023) `_ .. [Galbiati2023] - Marta Galbiati, Arianna Formenti, Mickael Grech and Matteo Passoni, + M. Galbiati, A. Formenti, M. Grech and M. Passoni, `Numerical investigation of non-linear inverse Compton scattering in double-layer targets`, `Front. Phys. 11, fphy.2023.1117543 (2023) `_ @@ -200,7 +278,7 @@ As of November 2021, 90 papers have been published covering a broad range of top .. [Miethlinger2023] - Thomas Miethlinger, Nico Hoffmann and Thomas Kluge, + T. Miethlinger, N. Hoffmann and T. Kluge, `Acceptance Rates of Invertible Neural Networks on Electron Spectra from Near-Critical Laser-Plasmas: A Comparison`, `Parallel Processing and Applied Mathematics, 273-284 (2023) `_ @@ -216,13 +294,30 @@ As of November 2021, 90 papers have been published covering a broad range of top `Electron acceleration by laser plasma wedge interaction`, `Phys. Rev. Research 5, 013115 (2023) `_ +.. [Blackman2022] + + D. R. Blackman, Y. Shi, S. R. Klein, M. Cernaianu, D. Doria, P. Ghenuche and A. Arefiev + `Electron acceleration from transparent targets irradiated by ultra-intense helical laser beams`, + `Communications Physics 5, 116 (2022) `_ + +.. [Siminos2022] + + E. Siminos and I. Thiele + `Parametric study of laser wakefield driven generation of intense sub-cycle pulses`, + `Plasma Physics and Controlled Fusion 64, 034006 (2022) `_ + +.. [PChen2022] + + P. Chen, G. Mourou, M. Besancon, Y. Fukuda, J.-F. Glicenstein, J. Nam, C.-E. Lin, K.-N. Lin, S.-X. Liu, Y.-K. Liu, M. Kando, K. Kondo, S. Paganis, A. Pirozhkov, H. Takabe, B. Tuchming, W.-P. Wang, N. Watamura, J. Wheeler and H.-Y. Wu on behalf of the AnaBHEL Collaboration, + `AnaBHEL (Analog Black Hole Evaporation via Lasers) Experiment: Concept, Design, and Status`, + `Photonics 9(12), 1003 (2022) `_ + .. [Bukharskii2022] N. Bukharskii, Iu. Kochetkov and Ph. Korneev, `Terahertz annular antenna driven with a short intense laser pulse`, `Applied Physics Letters 120, 014102 (2022) `_ - .. [Jirka2022] M. Jirka, P. Sasorov and S. V. Bulanov, @@ -273,19 +368,19 @@ As of November 2021, 90 papers have been published covering a broad range of top .. [Chen2022] - Q. Chen, Dominika Maslarova, J. Wang, S. Li, and D. Umstadter, + Q. Chen, D. Maslarova, J. Wang, S. Li, and D. Umstadter, `Injection of electron beams into two laser wakefields and generation of electron rings`, `Phys. Rev. E 106, 055202 (2022) `_ .. [Kumar2022b] - Sonu Kumar, Rajat Dhawan, D.K. Singh and Hitendra K. Malik, + 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`, `Materials Today: Proceedings 62, 3203-3207 (2022) `_ .. [Kumar2022a] - Sonu Kumar, Dhananjay K Singh and Hitendra K Malik, + S. Kumar, D. K. Singh and H. K. Malik, `Comparative study of ultrashort single-pulse and multi-pulse driven laser wakefield acceleration`, `Laser Phys. Lett. 20, 026001 (2022) `_ @@ -297,13 +392,13 @@ As of November 2021, 90 papers have been published covering a broad range of top .. [Zhang2022b] - Yue Zhang, Feng Wang, Jianyong Liu and Jizhong Sun, + Y. Zhang, F. Wang, J. Liu and J. Sun, `Simulation of the inverse bremsstrahlung absorption by plasma plume in laser penetration welding`, `Chemical Physics Letters 793, 139434 (2022) `_ .. [Vladisavlevici2022] - Iuliana-Mariana Vladisavlevici, Daniel Vizman and Emmanuel d’Humières, + I.-M. Vladisavlevici, D. Vizman and E. d’Humières, `Laser Driven Electron Acceleration from Near-Critical Density Targets towards the Generation of High Energy γ-Photons`, `Photonics 9, 953 (2022) `_ @@ -321,32 +416,32 @@ As of November 2021, 90 papers have been published covering a broad range of top .. [Guo2022] - Yinlong Guo, Xuesong Geng, Liangliang Ji, Baifei Shen and Ruxin Li, + Y. Guo, X. Geng, L. Ji, B. Shen and R. Li, `Improving the accuracy of hard photon emission by sigmoid sampling of the quantum-electrodynamic table in particle-in-cell Monte Carlo simulations`, `Phys. Rev. E 105, 025309 (2022) `_ .. [Pae2022] - Ki Hong Pae, Chul Min Kim, Vishwa Bandhu Pathak, Chang-Mo Ryu and Chang Hee Nam, + K. . Pae, C. M. Kim, V. B. Pathak, C.-M. Ryu and C. H. Nam, `Direct laser acceleration of electrons from a plasma mirror by an intense few-cycle Laguerre–Gaussian laser and its dependence on the carrier-envelope phase`, `Plasma Phys. Control. Fusion 64, 055013 (2022) `_ .. [Zhang2022a] - Cui-Wen Zhang, Yi-Xuan Zhu, Jian-Feng Lv and Bai-Song Xie, + C.-W. Zhang, Y.-X. Zhu, J.-F. Lu and B.-S. Xie, `Simulation Study of a Bright Attosecond γ-ray Source Generation by Irradiating an Intense Laser on a Cone Target`, `Applied Sciences 12, 4361 (2022) `_ .. [Han2022] - Qianqian Han, Xuesong Geng, Baifei Shen, Zhizhan Xu and Liangliang Ji, + Q. Han, X. Geng, B. Shen, Z. Xu and L. Ji, `Ultra-fast polarization of a thin electron layer in the rotational standing-wave field driven by double ultra-intense laser pulses`, `New J. Phys. 24, 063013 (2022) `_ .. [Gothel2022] - Ilja Göthel, Constantin Bernert, Michael Bussmann, Marco Garten, Thomas Miethlinger, Martin Rehwald, Karl Zeil, Tim Ziegler, Thomas E Cowan, Ulrich Schramm and Thomas Kluge, + I. Göthel, C. Bernert, M. Bussmann, M. Garten, T. Miethlinger, M. Rehwald, K. Zeil, T. Ziegler, T. E. Cowan, U. Schramm and T. Kluge, `Optimized laser ion acceleration at the relativistic critical density surface`, `Plasma Phys. Control. Fusion 64, 044010 (2022) `_ @@ -364,7 +459,7 @@ As of November 2021, 90 papers have been published covering a broad range of top .. [Sundstrom2022] - Andréas Sundström, Mickael Grech, István Pusztai and Caterina Riconda, + A. Sundström, M. Grech, I. Pusztai and C. Riconda, `Stimulated-Raman-scattering amplification of attosecond XUV pulses with pulse-train pumps and application to local in-depth plasma-density measurement`, `Phys. Rev. E 106, 045208 (2022) `_ @@ -382,13 +477,13 @@ As of November 2021, 90 papers have been published covering a broad range of top .. [Kong2022] - Defeng Kong, Guoqiang Zhang, Yinren Shou, Shirui Xu, Zhusong Mei, Zhengxuan Cao, Zhuo Pan, Pengjie Wang, Guijun Qi, Yao Lou, Zhiguo Ma, Haoyang Lan, Wenzhao Wang, Yunhui Li, Peter Rubovic, Martin Veselsky, Aldo Bonasera, Jiarui Zhao, Yixing Geng, Yanying Zhao, Changbo Fu, Wen Luo, Yugang Ma, Xueqing Yan and Wenjun Ma, + D. Kong, G. Zhang, Y. Shou, S. Xu, Z. Mei, Z. Cao, Z. Pan, P. Wang, G. Qi, Y. Lou, Z. Ma, H. Lan, W. Wang, Y. Li, P. Rubovic, M. Veselsky, A. Bonasera, J. Zhao, Y. Geng, Y. Zhao, C. Fu, W. Luo, Y. Ma, X. Yan and W. Ma, `High-energy-density plasma in femtosecond-laser-irradiated nanowire-array targets for nuclear reactions`, `Matter and Radiation at Extremes 7, 064403 (2022) `_ .. [Davidson2022] - Conor Davidson, Zheng-Ming Sheng, Thomas Wilson and Paul McKenna, + C. Davidson, Z.-M. Sheng, T. Wilson and P. McKenna, `Theoretical and computational studies of the Weibel instability in several beam–plasma interaction configurations`, `J. Plasma Phys. 88, 905880206 (2022) `_ @@ -399,13 +494,14 @@ As of November 2021, 90 papers have been published covering a broad range of top `Journal of Applied Physics 131, 103104 (2022) `_ .. [Umstadter2022] + D. Umstadter `Controlled Injection of Electrons for Improved Performance of Laser-Wakefield Acceleration`, `United States: N. p., (2022) `_ .. [Massimo2022] - Francesco Massimo, Mathieu Lobet, Julien Derouillat, Arnaud Beck, Guillaume Bouchard, Mickael Grech, Frédéric Pérez, Tommaso Vinci, + F. Massimo, M. Lobet, J. Derouillat, A. Beck, G. Bouchard, M. Grech, F. Pérez, T. Vinci, `A Task Programming Implementation for the Particle in Cell Code Smilei`, `PASC '22: Proceedings of the Platform for Advanced Scientific Computing Conference 5, 1 (2022) `_, `arXiv:2204.12837 `_ @@ -429,6 +525,12 @@ As of November 2021, 90 papers have been published covering a broad range of top `In International Conference on High Performance Computing in Asia-Pacific Region Workshops (HPCAsia 2022 Workshop). Association for Computing Machinery, New York, NY, USA, 40–48. (2022) `_ +.. [Romansky2021] + + V. I. Romansky, A. M. Bykov and S. M. Osipov, + `On electron acceleration by mildly-relativistic shocks: PIC simulations`, + `Journal of Physics: Conference Series 2103 012009 (2021) `_ + .. [Tiwary2021] S. Tiwary and N. Kumar, @@ -443,13 +545,13 @@ As of November 2021, 90 papers have been published covering a broad range of top .. [Tomassini2021] - Paolo Tomassini, Francesco Massimo, Luca Labate and Leonida A. Gizzi, + P. Tomassini, F. Massimo, L. Labate and L. A. Gizzi, `Accurate electron beam phase-space theory for ionization-injection schemes driven by laser pulses`, `High Pow Laser Sci Eng 10, e15 (2021) `_ .. [Meinhold2021] - Tim Arniko Meinhold and Naveen Kumar, + T. A. Meinhold and N. Kumar, `Radiation pressure acceleration of protons from structured thin-foil targets`, `J. Plasma Phys. 87, 905870607 (2021) `_ @@ -461,13 +563,13 @@ As of November 2021, 90 papers have been published covering a broad range of top .. [Shi2021b] - Yin Shi, David R Blackman and Alexey Arefiev, + Y. Shi, D. R. Blackman and A. Arefiev, `Electron acceleration using twisted laser wavefronts`, `Plasma Phys. Control. Fusion 63, 125032 (2021) `_ .. [Kumar2021] - Naveen Kumar and Brian Reville, + N. Kumar and B. Reville, `Nonthermal Particle Acceleration at Highly Oblique Nonrelativistic Shocks`, `ApJL 921, L14 (2021) `_ @@ -479,7 +581,7 @@ As of November 2021, 90 papers have been published covering a broad range of top .. [Horny2021] - Vojtěch Horný and László Veisz, + V. Horný and L. Veisz, `Generation of single attosecond relativistic electron bunch from intense laser interaction with a nanosphere`, `Plasma Phys. Control. Fusion 63, 125025 (2021) `_ @@ -509,7 +611,7 @@ As of November 2021, 90 papers have been published covering a broad range of top .. [Shou2021] - Yinren Shou, Dahui Wang, Pengjie Wang, Jianbo Liu, Zhengxuan Cao, Zhusong Mei, Shirui Xu, Zhuo Pan, Defeng Kong, Guijun Qi, Zhipeng Liu, Yulan Liang, Ziyang Peng, Ying Gao, Shiyou Chen, Jiarui Zhao, Yanying Zhao, Han Xu, Jun Zhao, Yanqing Wu, Xueqing Yan and Wenjun Ma, + Y. Shou, D. Wang, P. Wang, J. Liu, Z. Cao, Z. Mei, S. Xu, Z. Pan, D. Kong, G. Qi, Z. Liu, Y. Liang, Z. Peng, Y. Gao, S. Chen, J. Zhao, Y. Zhao, H. Xu, J. Zhao, Y. Wu, X. Yan and W. Ma, `High-efficiency generation of narrowband soft x rays from carbon nanotube foams irradiated by relativistic femtosecond lasers`, `Opt. Lett. 46, 3969 (2021) `_ @@ -521,13 +623,13 @@ As of November 2021, 90 papers have been published covering a broad range of top .. [Hosseinkhani2021] - H. Hosseinkhani, M. Pishdast, J. Yazdanpanah and S.A. Ghasemi, + H. Hosseinkhani, M. Pishdast, J. Yazdanpanah and S. A. Ghasemi, `Investigation of the classical and quantum radiation reaction effect on interaction of ultra high power laser with near critical plasma`, `J. Nuclear Sci. Technol. 42, 27-35 (2021) `_ .. [MercuriBaron2021] - A Mercuri-Baron, M Grech, F Niel, A Grassi, M Lobet, A Di Piazza and C Riconda, + A. Mercuri-Baron, M. Grech, F. Niel, A. Grassi, M. Lobet, A. Di Piazza and C. Riconda, `Impact of the laser spatio-temporal shape on Breit–Wheeler pair production`, `New J. Phys. 23, 085006 (2021) `_ @@ -539,7 +641,7 @@ As of November 2021, 90 papers have been published covering a broad range of top .. [Shi2021a] - Yin Shi, David Blackman, Dan Stutman and Alexey Arefiev, + Y. Shi, D. Blackman, D. Stutman and A. Arefiev, `Generation of Ultrarelativistic Monoenergetic Electron Bunches via a Synergistic Interaction of Longitudinal Electric and Magnetic Fields of a Twisted Laser`, `Phys. Rev. Lett. 126, 234801 (2021) `_ @@ -551,7 +653,7 @@ As of November 2021, 90 papers have been published covering a broad range of top .. [Shekhanov2021] - S A Shekhanov and V T Tikhonchuk, + S. A. Shekhanov and V. T. Tikhonchuk, `SRS-SBS competition and nonlinear laser energy absorption in a high temperature plasma`, `Plasma Phys. Control. Fusion 63, 115016 (2021) `_ @@ -563,7 +665,7 @@ As of November 2021, 90 papers have been published covering a broad range of top .. [Yoon2021b] - Young Dae Yoon, Gunsu S. Yun, Deirdre E. Wendel and James L. Burch, + Y. D. Yoon, G. S. Yun, D. E. Wendel and J. L. Burch, `Collisionless relaxation of a disequilibrated current sheet and implications for bifurcated structures`, `Nat Commun 12, 3774 (2021) `_ @@ -593,7 +695,7 @@ As of November 2021, 90 papers have been published covering a broad range of top .. [Do2021] - Hue Thi Bich Do, Ding Wen Jun, Zackaria Mahfoud, Wu Lin and Michel Bosman, + H. T. B. Do, D. W. Jun, Z. Mahfoud, W. Lin and M. Bosman, `Electron dynamics in plasmons`, `Nanoscale 13, 2801-2810 (2021) `_ @@ -611,7 +713,7 @@ As of November 2021, 90 papers have been published covering a broad range of top .. [Cantono2021] - Giada Cantono, Alexander Permogorov, Julien Ferri, Evgeniya Smetanina, Alexandre Dmitriev, Anders Persson, Tünde Fülöp and Claes-Göran Wahlström, + G. Cantono, A. Permogorov, J. Ferri, E. Smetanina, A. Dmitriev, A. Persson, T. Fülöp and C.-G. Wahlström, `Laser-driven proton acceleration from ultrathin foils with nanoholes`, `Sci Rep 11, 5006 (2021) `_ @@ -623,13 +725,13 @@ As of November 2021, 90 papers have been published covering a broad range of top .. [Yoon2021a] - Young Dae Yoon and Paul M. Bellan, + Y. D. Yoon and P. M. Bellan, `How Hall electric fields intrinsically chaotize and heat ions during collisionless magnetic reconnection`, `Physics of Plasmas 28, 022113 (2021) `_ .. [Sampath2021] - Archana Sampath, Xavier Davoine, Sébastien Corde, Laurent Gremillet, Max Gilljohann, Maitreyi Sangal, Christoph H. Keitel, Robert Ariniello, John Cary, Henrik Ekerfelt, Claudio Emma, Frederico Fiuza, Hiroki Fujii, Mark Hogan, Chan Joshi, Alexander Knetsch, Olena Kononenko, Valentina Lee, Mike Litos, Kenneth Marsh, Zan Nie, Brendan O’Shea, J. Ryan Peterson, Pablo San Miguel Claveria, Doug Storey, Yipeng Wu, Xinlu Xu, Chaojie Zhang and Matteo Tamburini, + A. Sampath, X. Davoine, S. Corde, L. Gremillet, M. Gilljohann, M. Sangal, C. H. Keitel, R. Ariniello, J. Cary, H. Ekerfelt, C. Emma, F. Fiuza, H. Fujii, M. Hogan, C. Joshi, A. Knetsch, O. Kononenko, V. Lee, M. Litos, K. Marsh, Z. Nie, B. O’Shea, J. R. Peterson, P. San Miguel Claveria, D. Storey, Y. Wu, X. Xu, C. Zhang and M. Tamburini, `Extremely Dense Gamma-Ray Pulses in Electron Beam-Multifoil Collisions`, `Phys. Rev. Lett. 126, 064801 (2021) `_ @@ -647,7 +749,7 @@ As of November 2021, 90 papers have been published covering a broad range of top .. [Gelfer2021] - E G Gelfer, A M Fedotov and S Weber, + E. G. Gelfer, A. M, Fedotov and S. Weber, `Radiation induced acceleration of ions in a laser irradiated transparent foil`, `New J. Phys. 23, 095002 (2021) `_ `arXiv:1907.02621 `_ @@ -754,7 +856,7 @@ As of November 2021, 90 papers have been published covering a broad range of top .. [Marcowith2020] - Alexandre Marcowith, Gilles Ferrand, Mickael Grech, Zakaria Meliani, Illya Plotnikov and Rolf Walder, + A. Marcowith, G. Ferrand, M. Grech, Z. Meliani, I. Plotnikov and R. Walder, `Multi-scale simulations of particle acceleration in astrophysical systems`, `Living Rev Comput Astrophys 6, 1 (2020) `_ `arXiv:2002.09411 `_ @@ -843,26 +945,26 @@ As of November 2021, 90 papers have been published covering a broad range of top .. [Fang2019] - Jun Fang, Chun-Yan Lu, Jing-Wen Yan and Huan Yu, + J. Fang, C.-Y. Lu, J.-W. Yan and H. Yu, `Early acceleration of electrons and protons at the nonrelativistic quasiparallel shocks with different obliquity angles`, `Res. Astron. Astrophys. 19, 182 (2019) `_ `arXiv:1908.08170 `_ .. [Yoon2019b] - Young Dae Yoon and Paul M. Bellan, + Y. Yoon and P. M. Bellan, `Kinetic Verification of the Stochastic Ion Heating Mechanism in Collisionless Magnetic Reconnection`, `ApJ 887, L29 (2019) `_ .. [Yoon2019a] - Young Dae Yoon and Paul M. Bellan, + Y. D. Yoon and P. M. Bellan, `The electron canonical battery effect in magnetic reconnection: Completion of the electron canonical vorticity framework`, `Physics of Plasmas 26, 100702 (2019) `_ .. [Massimo2019] - F Massimo, A Beck, J Derouillat, M Grech, M Lobet, F Pérez, I Zemzemi and A Specka, + F. Massimo, A. Beck, J. Derouillat, M. Grech, M. Lobet, F. Pérez, I. Zemzemi and A Specka, `Efficient start-to-end 3D envelope modeling for two-stage laser wakefield acceleration experiments`, `Plasma Phys. Control. Fusion 61, 124001 (2019) `_ `arXiv:1912.04127 `_ @@ -888,6 +990,12 @@ As of November 2021, 90 papers have been published covering a broad range of top `Phys. Rev. Lett. 122, 104803 (2019) `_ `arXiv:1806.04976 `_ +.. [Golovanov2018] + + A. A. Golovanov and I. Yu. Kostyukov, + `Bubble regime of plasma wakefield in 2D and 3D geometries`, + `Physics of Plasmas 25, 103107 (2018) `_ + .. [Massimo2018] F. Massimo, A. Beck, A. Specka, I. Zemzemi, J. Derouillat, M. Grech and F. Pérez, @@ -908,14 +1016,14 @@ As of November 2021, 90 papers have been published covering a broad range of top .. [Niel2018b] - F Niel, C Riconda, F Amiranoff, M Lobet, J Derouillat, F Pérez, T Vinci and M Grech, + F. Niel, C. Riconda, F. Amiranoff, M. Lobet, J. Derouillat, F. Pérez, T. Vinci and M. Grech, `From quantum to classical modeling of radiation reaction: a focus on the radiation spectrum`, `Plasma Phys. Control. Fusion 60, 094002 (2018) `_ `arXiv:1802.02927 `_ .. [Plotnikov2018] - Illya Plotnikov, Anna Grassi and Mickael Grech, + I. Plotnikov, A. Grassi and M. Grech, `Perpendicular relativistic shocks in magnetized pair plasma`, `Monthly Notices of the Royal Astronomical Society 477, 5238-5260 (2018) `_ `arXiv:1712.02883 `_ @@ -936,7 +1044,7 @@ As of November 2021, 90 papers have been published covering a broad range of top .. [Fedeli2017] - L Fedeli, A Formenti, L Cialfi, A Sgattoni, G Cantono and M Passoni, + L. Fedeli, A. Formenti, L. Cialfi, A. Sgattoni, G. Cantono and M. Passoni, `Structured targets for advanced laser-driven sources`, `Plasma Phys. Control. Fusion 60, 014013 (2017) `_ diff --git a/doc/Sphinx/Overview/partners.rst b/doc/Sphinx/Overview/partners.rst index 2170321ba..69b87e746 100755 --- a/doc/Sphinx/Overview/partners.rst +++ b/doc/Sphinx/Overview/partners.rst @@ -56,7 +56,6 @@ Partners | | * `Julien Dérouillat `_ | | | * `Haithem Kallala `_ | | | * `Mathieu Lobet `_ | -| | * `Francesco Massimo `_ | | | * `Charles Prouveur `_ | | | | +------------+---------------------------------------------------------------------------------------------------------+ diff --git a/doc/Sphinx/Overview/releases.rst b/doc/Sphinx/Overview/releases.rst index f9dfbc8cc..1b0bd9ba7 100755 --- a/doc/Sphinx/Overview/releases.rst +++ b/doc/Sphinx/Overview/releases.rst @@ -25,13 +25,21 @@ Changes made in the repository (not released) * Happi: + * In ``Scalar``, it is now possible to make an operation on scalars such as ``"Uelm+Ukin"``. + The list of available scalars can be obtained from ``getScalars()``. * New arguments ``xoffset`` and ``yoffset`` to shift plot coordinates. + * New argument ``timestep_indices`` as an alternative to ``timesteps``. * Changed coordinate reference for 2D probe in 3D or AM geometry (zero is the box origin projected orthogonally on the probe plane). * Documentation: - * Dark theme (switch on the bottom left, or set browser preferences). + * Dark theme (click the switch on the bottom left, or set browser preferences). + +* Bug fixes: + + * ``dump_minutes`` often failed to write some checkpoint files. + * ``"auto"`` limits in ``ParticleBinning`` could fail with only one side on ``"auto"``. ---- diff --git a/doc/Sphinx/Understand/azimuthal_modes_decomposition.rst b/doc/Sphinx/Understand/azimuthal_modes_decomposition.rst index 8cfdab1e6..deb171eca 100644 --- a/doc/Sphinx/Understand/azimuthal_modes_decomposition.rst +++ b/doc/Sphinx/Understand/azimuthal_modes_decomposition.rst @@ -320,7 +320,7 @@ Cancellation on axis """"""""""""""""""""""" The first basic principle is that a mode 0 field defined on axis can only be longitudinal otherwise it would be ill defined. -On the opposite, longitudinal fields on axis can only be of mode 0 since they do not depend on :math:`theta`. +On the opposite, longitudinal fields on axis can only be of mode 0 since they do not depend on :math:`\theta`. From this we can already state that :math:`E_r^{m=0},\ E_t^{m=0},\ B_r^{m=0},\ B_t^{m=0},\ E_l^{m>0},\ B_l^{m>0}` are zero on axis. @@ -415,7 +415,12 @@ To ensure this derivative cancels on axis we simply pick: And equation :eq:`transverse_on_axis` then gives .. math:: - \tilde{E_{\theta}}^{m=1}[2] = -i\tilde{E_r}^{m=1}[3] + \tilde{E_{\theta}}^{m=1}(r=0) = -i\tilde{E_r}^{m=1}(r=0) + +With a finite difference scheme, this is implemented as + +.. math:: + \tilde{E_{\theta}}^{m=1}[2] = -\frac{i}{8}(9\tilde{E_r}^{m=1}[3]-\tilde{E_r}^{m=1}[4]) All the equation derived here are also valid for the magnetic field. But because of a different duality, it is more convenient to use a different approach. @@ -450,7 +455,7 @@ With a similar interpolation we obtain the boundary condition on axis for :math: .. math:: B_{\theta}^{m=1}[2]=-2iB_{r}^{m=1}[2]-B_{\theta}^{m=1}[3] -Longitudinal fields on axis +Longitudinal field on axis """"""""""""""""""""""""""""" We have alreayd established that only modes :math:`m=0` of longitudinal fields are non zero on axis. @@ -468,11 +473,6 @@ Introducing this result in the standard FDTD expression of :math:`E_l` we get: Again, the :math:`n` indice indicates the time step here. -:math:`B_l^{m=0}` is independant of :math:`\theta`. If we assume it is differentiable at :math:`r=0` then its derivative along :math:`r` is zero -on axis (derivative of a pair function is zero at :math:`x=0` ). From this we get: - -.. math:: - B_{l}^{m=0}[2]=B_{l}^{m=0}[3] Below axis """""""""""""""""""""""""""" @@ -481,7 +481,8 @@ Fields "below" axis are primal fields data with indice :math:`j<2` and dual fiel These fields are not physical in the sense that they do not contribute to the reconstruction of any physical field in real space and are not obtained by solving Maxwell's equations. Nevertheless, it is numerically convenient to give them a value in order to facilitate field interpolation for macro-particles near axis. This is already what is done for dual fields in :math:`r` which cancel on axis for instance. -We extend this logic to primal fields in :math:`r`: +We extend this logic to primal fields in :math:`r` as well. +For any given field :math:`F`, the symetric of :math:`F` with respect to the axis is :math:`F` if :math:`F` is non zero on axis and :math:`-F` if :math:`F` is zero on axis: .. math:: @@ -489,19 +490,35 @@ We extend this logic to primal fields in :math:`r`: E_{l}^{m>0}[1] = -E_{l}^{m>0}[3] + E_{r}^{m\neq1}[2] = -E_{r}^{m\neq1}[3] + + E_{r}^{m=1}[2] = E_{r}^{m=1}[3] + E_{\theta}^{m\neq1}[1] = -E_{\theta}^{m\neq1}[3] E_{\theta}^{m=1}[1] = E_{\theta}^{m=1}[3] + B_{l}^{m=0}[2]=B_{l}^{m=0}[3] + + B_{l}^{m>0}[2]=-B_{l}^{m>0}[3] + B_{r}^{m\neq1}[1] = -B_{r}^{m\neq1}[3] B_{r}^{m=1}[1] = B_{r}^{m=1}[3] + B_{t}^{m\neq1}[2] = -B_{t}^{m\neq1}[3] + + B_{t}^{m=1}[2] = B_{t}^{m=1}[3] + Currents near axis """"""""""""""""""""" A specific treatment must be applied to charge and current densities near axis because the projector deposits charge and current "below" axis. -Quantities below axis must be brought back in the "physical" terms on and above axis. +Quantities below axis must be folded back onto their symetric "above" axis. +Mode 0 contributions below axis are added to their "above axis" counterparts. +On the opposite, strictly positive modes contributions are deduced. +This ensures that a particle sitting on axis will have a net contribution to the domain only in mode 0 as expected because in that case :math:`\theta` is not defined and therefore +the deposition can not be a function of :math:`\theta`. Using the continuity equation instead of Gauss law for transverse current of mode :math:`m=1` on axis, we can derive the exact same boundary conditions on axis for current density as for the electric field. diff --git a/doc/Sphinx/Use/installation.rst b/doc/Sphinx/Use/installation.rst index 2ce4fec64..a2f872e36 100755 --- a/doc/Sphinx/Use/installation.rst +++ b/doc/Sphinx/Use/installation.rst @@ -105,9 +105,6 @@ Advanced compilation options make config=vtune # For Intel Vtune make config=inspector # For Intel Inspector make config=detailed_timers # More detailed timers, but somewhat slower execution - make config=omptasks # use OpenMP task parallelization, not supported by old compilers - make config=part_event_tracing_tasks_off # trace the use particle operators, without task parallelization - make config=part_event_tracing_tasks_on # trace the use particle operators, with OpenMP task parallelization make config="gpu_nvidia noopenmp" # For Nvidia GPU acceleration make config="gpu_amd" # For AMD GPU acceleration @@ -143,6 +140,47 @@ executed before compilation. If you successfully write such a file for a common supercomputer, please share it with developpers so that it can be included in the next release of :program:`Smilei`. + +.. rubric:: Compilation for GPU accelerated nodes: + +As each supercomputer has a different environnment to compile for GPUs and since the nvhpc + CUDA/ cray + HIP modules evolve quickly, a machine file is required for the compilation. +Several machine files are already available as an example in smilei/scripts/compile_tools/machine/ ; such as: jean_zay_gpu_V100, jean_zay_gpu_A100, adastra, ruche_gpu2. + +Typically we need it to specify ACCELERATOR_GPU_FLAGS += -ta=tesla:cc80 for nvhpc <23.4 and ACCELERATOR_GPU_FLAGS += -gpu=cc80 -acc for the more recent versions of nvhpc. + +.. code-block:: bash + + make -j 12 machine="jean_zay_gpu_A100" config="gpu_nvidia noopenmp verbose" # for Nvidia GPU + make -j 12 machine="adastra" config="gpu_amd" # for AMD GPU + + +Furthermore, here are 2 examples of known working ennvironments, first for AMD GPUs, second for Nvidia GPUs: + +.. code-block:: bash + + module purge + module load craype-accel-amd-gfx90a craype-x86-trento + module load PrgEnv-cray/8.3.3 + module load cpe/23.02 + module load cray-mpich/8.1.24 cray-hdf5-parallel/1.12.2.1 cray-python/3.9.13.1 + module load amd-mixed/5.2.3 + +.. code-block:: bash + + module purge + module load anaconda-py3/2020.11 # python is fine as well if you can pip install the required modules + module load nvidia-compilers/23.1 + module load cuda/11.2 + module load openmpi/4.1.1-cuda + module load hdf5/1.12.0-mpi-cuda + # For HDF5, note that module show can give you the right path + export HDF5_ROOT_DIR=/DIRECTORY_NAME/hdf5/1.12.0/pgi-20.4-HASH/ + +Note: + +* we are aware of issues with CUDA >12.0, fixes are being tested but are not deployed yet. We recommend CUDA 11.x at the moment. +* The hdf5 module should be compiled with the nvidia/cray compiler ; openmpi as well, but depending on the nvhpc module it might not be needed as it can be included in the nvhpc module + ---- .. _vectorization_flags: diff --git a/doc/Sphinx/Use/post-processing.rst b/doc/Sphinx/Use/post-processing.rst index 78f5083c7..8f90f56f6 100755 --- a/doc/Sphinx/Use/post-processing.rst +++ b/doc/Sphinx/Use/post-processing.rst @@ -101,12 +101,16 @@ about the corresponding diagnostics in the simulation. Returns a list of available diagnostics of the given type * ``diagType``: The diagnostic type (``"Field"``, ``"Probe"``, etc.) + +.. rubric:: Information on specific diagnostics + +.. py:method:: getScalars() + + Returns a list of available scalars. .. py:method:: getTrackSpecies() Returns a list of available tracked species. - -.. rubric:: Information on specific diagnostics .. py:method:: fieldInfo(diag) @@ -145,12 +149,16 @@ Open a Scalar diagnostic .. py:method:: Scalar(scalar=None, timesteps=None, units=[""], data_log=False, data_transform=None, **kwargs) - * ``scalar``: The name of the scalar. - | If not given, then a list of available scalars is printed. - * ``timesteps``: The requested timestep(s). - | If omitted, all timesteps are used. - | If one number given, the nearest timestep available is used. - | If two numbers given, all the timesteps in between are used. + * ``scalar``: The name of the scalar, or an operation on scalars, such as ``"Uelm+Ukin"``. + * ``timesteps`` or ``timestep_indices``: The requested range of timesteps. + + * If omitted, all timesteps are used. + * If one number given, the nearest timestep available is used. + * If two numbers given, all the timesteps in between are used. + + When using ``timesteps``, provide the timesteps themselves, but + when using ``timestep_indices``, provide their indices in the list + of the available timesteps. * ``units``: A unit specification (see :ref:`units`) * ``data_log``: | If ``True``, then :math:`\log_{10}` is applied to the output. @@ -170,7 +178,7 @@ Open a Field diagnostic .. py:method:: Field(diagNumber=None, field=None, timesteps=None, subset=None, average=None, units=[""], data_log=False, data_transform=None, moving=False, export_dir=None, **kwargs) - * ``timesteps``, ``units``, ``data_log``, ``data_transform``: same as before. + * ``timesteps`` (or ``timestep_indices``), ``units``, ``data_log``, ``data_transform``: same as before. * ``diagNumber``: number or ``name`` of the fields diagnostic | If not given, then a list of available diagnostic numbers is printed. * ``field``: The name of a field (``"Ex"``, ``"Ey"``, etc.) @@ -228,7 +236,7 @@ Open a Probe diagnostic .. py:method:: Probe(probeNumber=None, field=None, timesteps=None, subset=None, average=None, units=[""], data_log=False, data_transform=None, **kwargs) - * ``timesteps``, ``units``, ``data_log``, ``data_transform``, ``export_dir``: same as before. + * ``timesteps`` (or ``timestep_indices``), ``units``, ``data_log``, ``data_transform``, ``export_dir``: same as before. * ``probeNumber``: number or ``name`` of the probe (the first one has number 0). | If not given, a list of available probes is printed. * ``field``: name of the field (``"Bx"``, ``"By"``, ``"Bz"``, ``"Ex"``, ``"Ey"``, ``"Ez"``, ``"Jx"``, ``"Jy"``, ``"Jz"`` or ``"Rho"``). @@ -256,7 +264,7 @@ Open a ParticleBinning diagnostic .. py:method:: ParticleBinning(diagNumber=None, timesteps=None, subset=None, average=None, units=[""], data_log=False, data_transform=None, **kwargs) - * ``timesteps``, ``units``, ``data_log``, ``data_transform``, ``export_dir``: same as before. + * ``timesteps`` (or ``timestep_indices``), ``units``, ``data_log``, ``data_transform``, ``export_dir``: same as before. * ``diagNumber``: number or ``name`` of the particle binning diagnostic (starts at 0). | If not given, a list of available diagnostics is printed. | It can also be an operation between several diagnostics. @@ -306,7 +314,7 @@ Open a Screen diagnostic .. py:method:: Screen(diagNumber=None, timesteps=None, subset=None, average=None, units=[""], data_log=False, data_transform=None, **kwargs) - * ``timesteps``, ``units``, ``data_log``, ``data_transform``, ``export_dir``: same as before. + * ``timesteps`` (or ``timestep_indices``), ``units``, ``data_log``, ``data_transform``, ``export_dir``: same as before. * ``diagNumber``, ``subset`` and ``average``: identical to that of ParticleBinning diagnostics. * See also :ref:`otherkwargs` @@ -323,7 +331,7 @@ Open a RadiationSpectrum diagnostic .. py:method:: ParticleBinning(diagNumber=None, timesteps=None, subset=None, average=None, units=[""], data_log=False, data_transform=None, **kwargs) - * ``timesteps``, ``units``, ``data_log``, ``data_transform``, ``export_dir``: same as before. + * ``timesteps`` (or ``timestep_indices``), ``units``, ``data_log``, ``data_transform``, ``export_dir``: same as before. * ``diagNumber``, ``subset`` and ``average``: identical to that of ParticleBinning diagnostics. * See also :ref:`otherkwargs` @@ -344,7 +352,7 @@ Open a TrackParticles diagnostic .. py:method:: TrackParticles(species=None, select="", axes=[], timesteps=None, sort=True, length=None, units=[""], **kwargs) - * ``timesteps``, ``units``, ``export_dir``: same as before. + * ``timesteps`` (or ``timestep_indices``), ``units``, ``export_dir``: same as before. * ``species``: the name of a tracked-particle species. If omitted, a list of available tracked-particle species is printed. * ``select``: Instructions for selecting particles among those available. diff --git a/doc/Sphinx/Use/run.rst b/doc/Sphinx/Use/run.rst index 65a515da1..f5eb48198 100755 --- a/doc/Sphinx/Use/run.rst +++ b/doc/Sphinx/Use/run.rst @@ -141,10 +141,14 @@ With AMD GPUs using cray on Adastra: For the binding scripts themselves, as it depends completely on the node architecture, please contact your admin support team. + +Binding script for adastra can be found here: https://github.com/SmileiPIC/Smilei/issues/672#issuecomment-1820677606 with the example of a slurm script. +it can be used as a template for other AMD GPUs based supercomputers/clusters. Be aware that GPU support is in development and not all features are currently available. Please refer to the list of current supported features. + ---- Debugging diff --git a/happi/_Diagnostics/Diagnostic.py b/happi/_Diagnostics/Diagnostic.py index e588189dd..b17e4c449 100755 --- a/happi/_Diagnostics/Diagnostic.py +++ b/happi/_Diagnostics/Diagnostic.py @@ -1,5 +1,7 @@ from .._Utils import * +from matplotlib import ticker + class Diagnostic(object): """Mother class for all Diagnostics. To create a diagnostic, refer to the doc of the SmileiSimulation class. @@ -24,7 +26,6 @@ def __init__(self, simulation, *args, **kwargs): self._log = [] self._data_log = False self._data_transform = None - self._error = [] self._xoffset = 0. # The 'simulation' is a SmileiSimulation object. It is passed as an instance attribute @@ -41,9 +42,7 @@ def __init__(self, simulation, *args, **kwargs): # Reload the simulation, in case it has been updated self.simulation.reload() - if not self.simulation.valid: - self._error += ["Invalid Smilei simulation"] - return + assert self.simulation.valid, "Invalid Smilei simulation" # Copy some parameters from the simulation self._results_path = self.simulation._results_path @@ -63,21 +62,20 @@ def __init__(self, simulation, *args, **kwargs): self.units = kwargs.pop("units", [""]) if type(self.units) in [list, tuple]: self.units = Units(*self.units , verbose = self._verbose) if type(self.units) is dict : self.units = Units(verbose = self._verbose, **self.units) - if type(self.units) is not Units: - self._error += ["Could not understand the 'units' argument"] - return + assert type(self.units) is Units, "Could not understand the 'units' argument" self.units._initRegistry(self._ureg) + # Manage the 'timesteps' or 'timestep_indices' + assert not ( "timesteps" in kwargs and "timestep_indices" in kwargs ),\ + "Cannot define both `timesteps` and `timestep_indices`" + # Call the '_init' function of the child class remaining_kwargs = self._init(*args, **kwargs) - if remaining_kwargs is not None and len(remaining_kwargs) > 0: - self.valid = False - self._error += ["The following keyword-arguments are unknown: "+", ".join(remaining_kwargs.keys())] - return + assert remaining_kwargs is None or len(remaining_kwargs) == 0,\ + "The following keyword-arguments are unknown: "+", ".join(remaining_kwargs.keys()) self.dim = len(self._shape) - if self.valid: - self._prepareUnits() + self._prepareUnits() # Prepare data_log output self._dataAtTime = self._dataLogAtTime if self._data_log else self._dataLinAtTime @@ -86,19 +84,6 @@ def __init__(self, simulation, *args, **kwargs): def __repr__(self): self.info() return "" - - # Method to verify everything was ok during initialization - def _validate(self): - try: - self.simulation.valid - except Exception as e: - print("No valid Smilei simulation selected") - return False - if not self.simulation.valid or not self.valid: - print("***ERROR*** - Diagnostic is invalid") - print("\n".join(self._error)) - return False - return True # Prepare units for axes def _prepareUnits(self): @@ -141,7 +126,7 @@ def limits(self): # Method to print info on this diag def info(self): - if self._validate() and self._verbose: + if self._verbose: print(self._info()) # Method to get only the arrays of data @@ -157,7 +142,6 @@ def getData(self, timestep=None): A list of arrays: each array corresponding to the diagnostic data at a given timestep. """ - if not self._validate(): return self._prepare1() # prepare the vfactor data = [] @@ -173,7 +157,6 @@ def getData(self, timestep=None): def getTimesteps(self): """Obtains the list of timesteps selected in this diagnostic""" - if not self._validate(): return [] return self._timesteps def getTimes(self): @@ -182,7 +165,6 @@ def getTimes(self): By default, times are in the code's units, but are converted to the diagnostic's units defined by the `units` argument, if provided. """ - if not self._validate(): return [] return self.units.tcoeff * self.timestep * self._np.array(self._timesteps) def _getCenters(self, axis_index, timestep): @@ -232,7 +214,6 @@ def get(self, timestep=None): !!! Deprecated !!! Use functions `getData`, `getTimesteps`, `getTimes` and `getAxis` instead. """ - if not self._validate(): return # obtain the data arrays data = self.getData(timestep=timestep) # format the results into a dictionary @@ -284,7 +265,6 @@ def plot(self, timestep=None, saveAs=None, axes=None, dpi=200, **kwargs): S = happi.Open("path/to/my/results") S.ParticleBinning(1).plot(vmin=0, vmax=1e14) """ - if not self._validate(): return if not self._prepare(): return if not self._setAndCheck(**kwargs): return self.info() @@ -332,7 +312,6 @@ def streak(self, saveAs=None, axes=None, dpi=200, **kwargs): S = happi.Open("path/to/my/results") S.ParticleBinning(1).streak(vmin=0, vmax=1e14) """ - if not self._validate(): return if not self._prepare(): return if not self._setAndCheck(**kwargs): return self.info() @@ -420,7 +399,6 @@ def animate(self, movie="", fps=15, dpi=200, saveAs=None, axes=None, **kwargs): This takes the particle binning diagnostic #1 and plots the resulting array in figure 1 from 0 to 3e14. """ - if not self._validate(): return if not self._prepare(): return if not self._setAndCheck(**kwargs): return self.info() @@ -480,7 +458,6 @@ def slide(self, axes=None, **kwargs): S = happi.Open("path/to/my/results") S.ParticleBinning(1).slide(vmin=0, vmax=1e14) """ - if not self._validate(): return if not self._prepare(): return if not self._setAndCheck(**kwargs): return ax = self._make_axes(axes) @@ -514,20 +491,26 @@ def update(t): self.info() + # Method to select specific timesteps (or indices) among those available + def _selectTimesteps(self, timestep_selection, timestep_index_selection, timesteps): + if timestep_selection is not None: + ts = self._np.array(timestep_selection, ndmin=1, dtype=float) + if ts.size==2: # get all times in between bounds + return timesteps[ (timesteps>=ts[0]) * (timesteps<=ts[1]) ] + elif ts.size==1: # get nearest time + return self._np.array([timesteps[(self._np.abs(timesteps-ts)).argmin()]]) + else: + raise Exception("Argument `timesteps` must be one or two non-negative integers") + elif timestep_index_selection is not None: + ts = self._np.array(timestep_index_selection, ndmin=1, dtype=int) + if ts.size==2: + return timesteps[ ts[0]:ts[1] ] + elif ts.size==1: + return self._np.array(timesteps[ts[0]]) + else: + raise Exception("Argument `timestep_indices` must be one or two non-negative integers") + return timesteps - # Method to select specific timesteps among those available in times - def _selectTimesteps(self, timesteps, times): - ts = self._np.array(self._np.double(timesteps),ndmin=1) - if ts.size==2: - # get all times in between bounds - times = times[ (times>=ts[0]) * (times<=ts[1]) ] - elif ts.size==1: - # get nearest time - times = self._np.array([times[(self._np.abs(times-ts)).argmin()]]) - else: - raise - return times - # Method to select portion of a mesh based on a slice def _selectSubset(self, portion, meshpoints, axisname, axisunits, operation): try: @@ -832,23 +815,18 @@ def _setAxesOptions(self, ax): for option, value in self.options.labels.items(): getattr(ax, "set_"+option)( value, self.options.labels_font[option] ) # Ticklabels + fonts - for option, value in self.options.ticklabels_font.items(): - if option in self.options.ticklabels: - getattr(ax, "set_"+option)( value, self.options.ticklabels_font[option] ) - else: # manage tick label fonts even when not setting tick labels first - ticklabels = getattr(ax, "get_"+option)() - self._plt.setp(ticklabels, **self.options.ticklabels_font[option]) + for tl in ["xticklabels", "yticklabels"]: + font = self.options.ticklabels_font[tl] if tl in self.options.ticklabels_font else {} + if tl in self.options.ticklabels: + getattr(ax, "set_"+tl)( self.options.ticklabels[tl], fontdict=font ) + elif font: # manage tick label fonts even when not setting tick labels first + ticklabels = getattr(ax, "get_"+tl)() + self._plt.setp(ticklabels, **font) # Tick formatting - try: - if self.options.xtick: ax.ticklabel_format(axis="x",**self.options.xtick) - except Exception as e: - if self._verbose: print("Cannot format x ticks (typically happens with log-scale)") - self.options.xtick = [] - try: - if self.options.ytick: ax.ticklabel_format(axis="y",**self.options.ytick) - except Exception as e: - if self._verbose: print("Cannot format y ticks (typically happens with log-scale)") - self.options.ytick = [] + if type(ax.xaxis.get_major_formatter()) == ticker.ScalarFormatter: + ax.ticklabel_format(axis="x",**self.options.xtick) + if type(ax.yaxis.get_major_formatter()) == ticker.ScalarFormatter: + ax.ticklabel_format(axis="y",**self.options.ytick) def _setColorbarOptions(self, ax): # Colorbar tick font if self.options.colorbar_font: @@ -880,7 +858,6 @@ def _dataLogAtTime(self, t): # Convert data to VTK format def toVTK(self, numberOfPieces=1): - if not self._validate(): return # prepare vfactor self._prepare1() diff --git a/happi/_Diagnostics/Field.py b/happi/_Diagnostics/Field.py index cf65a5faa..4b9a678c8 100755 --- a/happi/_Diagnostics/Field.py +++ b/happi/_Diagnostics/Field.py @@ -176,11 +176,8 @@ def fieldTranslator(f): for i,t in enumerate(self._timesteps): self._data.update({ t : i }) # If timesteps is None, then keep all timesteps otherwise, select timesteps - if timesteps is not None: - try: - self._timesteps = self._selectTimesteps(timesteps, self._timesteps) - except Exception as e: - raise Exception("Argument `timesteps` must be one or two non-negative integers") + timestep_indices = kwargs.pop("timestep_indices", None) + self._timesteps = self._selectTimesteps(timesteps, timestep_indices, self._timesteps) # Need at least one timestep if self._timesteps.size < 1: @@ -316,7 +313,6 @@ def _getCenters(self, axis_index, timestep, h5item = None): # get the value of x_moved for a requested timestep def getXmoved(self, t): - if not self._validate(): return # Verify that the timestep is valid if t not in self._timesteps: print("Timestep "+str(t)+" not found in this diagnostic") @@ -330,7 +326,6 @@ def getXmoved(self, t): # Method to obtain the data only def _getDataAtTime(self, t): - if not self._validate(): return # Verify that the timestep is valid if t not in self._timesteps: print("Timestep "+str(t)+" not found in this diagnostic") @@ -372,7 +367,6 @@ def _getDataAtTime(self, t): # Method to obtain the data only # Specific to cylindrical geometry, to reconstruct the plane at some theta def _theta_getDataAtTime(self, t): - if not self._validate(): return # Verify that the timestep is valid if t not in self._timesteps: print("Timestep "+str(t)+" not found in this diagnostic") @@ -429,7 +423,6 @@ def _theta_getDataAtTime(self, t): # Specific to cylindrical geometry, to reconstruct a 3d box def _build3d_getDataAtTime(self, t): from scipy.interpolate import RegularGridInterpolator - if not self._validate(): return # Verify that the timestep is valid if t not in self._timesteps: print("Timestep "+str(t)+" not found in this diagnostic") diff --git a/happi/_Diagnostics/NewParticles.py b/happi/_Diagnostics/NewParticles.py index 8bc6285bc..6b45abc32 100755 --- a/happi/_Diagnostics/NewParticles.py +++ b/happi/_Diagnostics/NewParticles.py @@ -133,13 +133,14 @@ def _init(self, species=None, select="", axes=[], chunksize=20000000, **kwargs): print("Kept "+str(self.nselectedParticles)+" particles") # Finish constructor + assert "timesteps" not in kwargs and "timestep_indices" not in kwargs,\ + "The NewParticles diagnostic does not accept 'timesteps' or 'timestep_indices' as arguments." self._timesteps = [None] self.valid = True return kwargs # We override the get and getData methods def getData(self): - if not self._validate(): return for data in self.iterParticles(chunksize=self.nParticles): return data diff --git a/happi/_Diagnostics/ParticleBinning.py b/happi/_Diagnostics/ParticleBinning.py index 851442ada..bb4cfbf2d 100755 --- a/happi/_Diagnostics/ParticleBinning.py +++ b/happi/_Diagnostics/ParticleBinning.py @@ -95,13 +95,16 @@ def _init(self, diagNumber=None, timesteps=None, subset=None, average=None, sum= self._timesteps = {} self._alltimesteps = {} self._indexOfTime = {} + self._h5files = [] self._h5items = {} + timestep_indices = kwargs.pop("timestep_indices", None) for d in self._diags: # Gather data from all timesteps, and the list of timesteps items = {} for path in self._results_path: try: f = self._h5py.File(path+self._os.sep+self._diagType+str(d)+'.h5', 'r') + self._h5files += [f] except: continue items.update( dict(f) ) @@ -113,12 +116,8 @@ def _init(self, diagNumber=None, timesteps=None, subset=None, average=None, sum= self._indexOfTime.update({ d:{} }) for i,t in enumerate(self._timesteps[d]): self._indexOfTime[d].update({ t : i }) - # If timesteps is None, then keep all timesteps, otherwise, select timesteps - if timesteps is not None: - try: - self._timesteps[d] = self._selectTimesteps(timesteps, self._timesteps[d]) - except Exception as e: - raise Exception("Argument 'timesteps' must be one or two non-negative integers") + # Select timesteps if requested + self._timesteps[d] = self._selectTimesteps(timesteps, timestep_indices, self._timesteps[d]) # Verify that timesteps are the same for all diagnostics if (self._timesteps[d] != self._timesteps[self._diags[0]]).any() : raise Exception( @@ -493,7 +492,6 @@ def getAvailableTimesteps(self, diagNumber=None): # Method to obtain the data only def _getDataAtTime(self, t): - if not self._validate(): return # Auto axes require recalculation of bin size and centers if self.auto_axes: self._updateAxes(t) diff --git a/happi/_Diagnostics/Performances.py b/happi/_Diagnostics/Performances.py index 460fe44e9..c7a8df254 100755 --- a/happi/_Diagnostics/Performances.py +++ b/happi/_Diagnostics/Performances.py @@ -191,16 +191,10 @@ def _init(self, raw=None, map=None, histogram=None, timesteps=None, data_log=Fal self._data = {} for i,t in enumerate(self._timesteps): self._data.update({ t : i }) - # If timesteps is None, then keep all timesteps otherwise, select timesteps - if timesteps is not None: - try: - self._timesteps = self._selectTimesteps(timesteps, self._timesteps) - except Exception as e: - raise Exception("Argument `timesteps` must be one or two non-negative integers") - - # Need at least one timestep - if self._timesteps.size < 1: - raise Exception("Timesteps not found") + # Select timesteps if requested + timestep_indices = kwargs.pop("timestep_indices", None) + self._timesteps = self._selectTimesteps(timesteps, timestep_indices, self._timesteps) + assert self._timesteps.size > 0, "Timesteps not found" # 3 - Manage axes # ------------------------------------------------------------------- @@ -270,7 +264,6 @@ def getAvailableQuantities(self): # Method to obtain the data only def _getDataAtTime(self, t): - if not self._validate(): return # Verify that the timestep is valid if t not in self._timesteps: print("Timestep "+str(t)+" not found in this diagnostic") diff --git a/happi/_Diagnostics/Probe.py b/happi/_Diagnostics/Probe.py index 1fd4958eb..3c6b2c73d 100755 --- a/happi/_Diagnostics/Probe.py +++ b/happi/_Diagnostics/Probe.py @@ -92,16 +92,10 @@ def _init(self, requestedProbe=None, field=None, timesteps=None, subset=None, av # ------------------------------------------------------------------- # If timesteps is None, then keep all timesteps otherwise, select timesteps self._timesteps = self._alltimesteps - if timesteps is not None: - try: - self._timesteps = self._selectTimesteps(timesteps, self._timesteps) - except Exception as e: - raise Exception("Argument `timesteps` must be one or two non-negative integers") - - # Need at least one timestep - if self._timesteps.size < 1: - raise Exception("Timesteps not found") - + timestep_indices = kwargs.pop("timestep_indices", None) + self._timesteps = self._selectTimesteps(timesteps, timestep_indices, self._timesteps) + assert self._timesteps.size > 0, "Timesteps not found" + # 3 - Manage axes # ------------------------------------------------------------------- # Fabricate all axes values @@ -362,7 +356,6 @@ def getFields(self): # get the value of x_moved for a requested timestep def getXmoved(self, t): - if not self._validate(): return # Verify that the timestep is valid if t not in self._timesteps: print("Timestep "+str(t)+" not found in this diagnostic") @@ -379,7 +372,6 @@ def getAvailableTimesteps(self): # Method to obtain the data only def _getDataAtTime(self, t): - if not self._validate(): return # Verify that the timestep is valid if t not in self._timesteps: print("Timestep "+t+" not found in this diagnostic") diff --git a/happi/_Diagnostics/Scalar.py b/happi/_Diagnostics/Scalar.py index 550420e9d..bbacb6484 100755 --- a/happi/_Diagnostics/Scalar.py +++ b/happi/_Diagnostics/Scalar.py @@ -6,39 +6,55 @@ class Scalar(Diagnostic): def _init(self, scalar=None, timesteps=None, data_log=False, data_transform=None, **kwargs): # Get available scalars - scalars = self.getScalars() + scalars = self.simulation.getScalars() # If no scalar chosen, only print the available scalars if scalar is None: + error = [] if len(scalars)>0: - self._error += ["Error: no scalar chosen"] - self._error += ["Printing available scalars:"] - self._error += ["---------------------------"] + error += ["Error: no scalar chosen"] + error += ["Printing available scalars:"] + error += ["---------------------------"] l = [""] for s in scalars: if len(s)>4 and s[:2]!=l[-1][:2] and s[-2:]!=l[-1][-2:]: - if l!=[""]: self._error += ["\t".join(l)] + if l!=[""]: error += ["\t".join(l)] l = [] l.append(s) - if l!=[""]: self._error += ["\t".join(l)] + if l!=[""]: error += ["\t".join(l)] else: - self._error += ["No scalars found"] - return + error += ["No scalars found"] + raise Exception("\n".join(error)) # 1 - verifications, initialization # ------------------------------------------------------------------- - # Find which scalar is requested - if scalar not in scalars: - fs = list(filter(lambda x:scalar in x, scalars)) - if len(fs)==0: - self._error += ["No scalar `"+scalar+"` found"] - return - if len(fs)>1: - self._error += ["Several scalars match: "+(' '.join(fs))] - self._error += ["Please be more specific and retry."] - return - scalar = fs[0] - self._scalarname = scalar + # Find which scalar is requested & associated units + self.operation = scalar + def scalarTranslator(s): + units = "??" + if s == "time": + units = "T_r" + elif s == "Ubal_norm": + units = "" + else: + units = { + "U":"K_r * N_r * L_r^%i" % self._ndim_particles, + "P":"K_r * N_r * L_r^%i" % self._ndim_particles, + "D":"N_r * L_r^%i" % self._ndim_particles, + "E":"E_r", + "B":"B_r", + "J":"J_r", + "R":"Q_r * N_r", + "Z":"Q_r", + "N":"", + }[s[0]] + return units, "S['%s']"%s, s + self._operation = Operation(scalar, scalarTranslator, self._ureg) + self._scalarname = self._operation.variables + self._vunits = self._operation.translated_units + self._title = self._operation.title + if not self._scalarname: + raise Exception("String "+self.operation+" does not seem to include any scalar") # Put data_log as object's variable self._data_log = data_log @@ -47,26 +63,31 @@ def _init(self, scalar=None, timesteps=None, data_log=False, data_transform=None # Already get the data from the file # Loop file line by line self._alltimesteps = [] - self._values = [] - times_values = {} + S = { s:[] for s in self._scalarname } for path in self._results_path: try: - with open(path+'/scalars.txt') as f: - for line in f: - line = line.strip() - if line[0]!="#": break - prevline = line - scalars = prevline[1:].strip().split() # list of scalars - scalarindex = scalars.index(scalar) # index of the requested scalar - line = str(line.strip()).split() - times_values[ int( self._np.round(float(line[0]) / float(self.timestep)) ) ] = float(line[scalarindex]) - for line in f: - line = str(line.strip()).split() - times_values[ int( self._np.round(float(line[0]) / float(self.timestep)) ) ] = float(line[scalarindex]) + f = open(path+'/scalars.txt') except: continue - self._alltimesteps = self._np.array(sorted(times_values.keys())) - self._values = self._np.array([times_values[k] for k in self._alltimesteps]) + for line in f: + if line[0]!="#": break + prevline = line + scalars = prevline[1:].strip().split() # list of scalars + scalarindexes = {s:scalars.index(s) for s in self._scalarname} # indexes of the requested scalars + line = str(line.strip()).split() + self._alltimesteps += [ int( self._np.round(float(line[0]) / float(self.timestep)) ) ] + for s,i in scalarindexes.items(): + S[s] += [float(line[i])] + for line in f: + line = str(line.strip()).split() + self._alltimesteps += [ int( self._np.round(float(line[0]) / float(self.timestep)) ) ] + for s,i in scalarindexes.items(): + S[s] += [float(line[i])] + f.close() + self._alltimesteps = self._np.array(self._alltimesteps) + for s in S: + S[s] = self._np.array(S[s]) + self._values = self._operation.eval(locals()) self._timesteps = self._np.copy(self._alltimesteps) # 2 - Manage timesteps @@ -76,39 +97,11 @@ def _init(self, scalar=None, timesteps=None, data_log=False, data_transform=None for i,t in enumerate(self._timesteps): self._data.update({ t : i }) # If timesteps is None, then keep all timesteps otherwise, select timesteps - if timesteps is not None: - try: - self._timesteps = self._selectTimesteps(timesteps, self._timesteps) - except Exception as e: - self._error += ["Argument `timesteps` must be one or two non-negative integers"] - return + timestep_indices = kwargs.pop("timestep_indices", None) + self._timesteps = self._selectTimesteps(timesteps, timestep_indices, self._timesteps) # Need at least one timestep - if self._timesteps.size < 1: - self._error += ["Timesteps not found"] - return - - - # 3 - Build units - # ------------------------------------------------------------------- - self._vunits = "??" - if self._scalarname == "time": - self._vunits = "T_r" - elif self._scalarname == "Ubal_norm": - self._vunits = "" - else: - self._vunits = { - "U":"K_r * N_r * L_r^%i" % self._ndim_particles, - "P":"K_r * N_r * L_r^%i" % self._ndim_particles, - "D":"N_r * L_r^%i" % self._ndim_particles, - "E":"E_r", - "B":"B_r", - "J":"J_r", - "R":"Q_r * N_r", - "Z":"Q_r", - "N":"", - }[self._scalarname[0]] - self._title =self._scalarname + assert self._timesteps.size > 0, "Timesteps not found" # Finish constructor self.valid = True @@ -116,37 +109,7 @@ def _init(self, scalar=None, timesteps=None, data_log=False, data_transform=None # Method to print info on included scalars def _info(self): - return "Scalar "+self._scalarname - - # get all available scalars - def getScalars(self): - allScalars = None - for path in self._results_path: - try: - file = path+'/scalars.txt' - f = open(file, 'r') - except Exception as e: - continue - try: - # Find last commented line - prevline = "" - for line in f: - line = line.strip() - if line[0]!="#": break - prevline = line[1:].strip() - scalars = str(prevline).split() # list of scalars - scalars = scalars[1:] # remove first, which is "time" - except Exception as e: - scalars = [] - f.close() - if allScalars is None: - allScalars = scalars - else: - allScalars = self._np.intersect1d(allScalars, scalars) - if allScalars is None: - self._error += ["Cannot open 'scalars.txt'"] - return [] - return allScalars + return "Scalar "+self.operation # get all available timesteps def getAvailableTimesteps(self): @@ -154,7 +117,6 @@ def getAvailableTimesteps(self): # Method to obtain the data only def _getDataAtTime(self, t): - if not self._validate(): return # Verify that the timestep is valid if t not in self._timesteps: print("Timestep "+str(t)+" not found in this diagnostic") diff --git a/happi/_Diagnostics/TrackParticles.py b/happi/_Diagnostics/TrackParticles.py index 2d4e8bfb4..de5a63c55 100755 --- a/happi/_Diagnostics/TrackParticles.py +++ b/happi/_Diagnostics/TrackParticles.py @@ -27,6 +27,8 @@ class TrackParticles(ParticleList): def _init(self, species=None, select="", axes=[], timesteps=None, sort=True, sorted_as="", length=None, chunksize=20000000, **kwargs): + timestep_indices = kwargs.pop("timestep_indices", None) + # If argument 'species' not provided, then print available species and leave if species is None: species = self.simulation.getTrackSpecies() @@ -89,13 +91,13 @@ def _init(self, species=None, select="", axes=[], timesteps=None, sort=True, sor self._timesteps = self._np.array(sorted(self._locationForTime)) # If specific timesteps requested, narrow the selection - if sorted_as and timesteps is not None: - self._timesteps = self._filterTimesteps( self._timesteps, timesteps ) - + if sorted_as: + self._timesteps = self._selectTimesteps(timesteps, timestep_indices, self._timesteps) + assert self._timesteps.size > 0, "Timesteps not found" + self._alltimesteps = self._np.copy(self._timesteps) - if not self._locationForTime: - raise Exception("No data found for this diagnostic") + assert self._locationForTime, "No data found for this diagnostic" # List available properties try: # python 2 @@ -132,11 +134,8 @@ def _init(self, species=None, select="", axes=[], timesteps=None, sort=True, sor if self._timesteps.size == 0: raise Exception("No tracked particles found") # If specific timesteps requested, narrow the selection - if timesteps is not None: - self._timesteps = self._filterTimesteps( self._timesteps, timesteps ) - # Need at least one timestep - if self._timesteps.size < 1: - raise Exception("Timesteps not found") + self._timesteps = self._selectTimesteps(timesteps, timestep_indices, self._timesteps) + assert self._timesteps.size > 0, "Timesteps not found" # Select particles @@ -522,7 +521,6 @@ def _orderFiles( self, fileOrdered, chunksize, sort ): # Method to generate the raw data (only done once) def _generateRawData(self, times=None): - if not self._validate(): return self._prepare1() # prepare the vfactor if self._sort: @@ -594,7 +592,6 @@ def _generateRawData(self, times=None): # We override the get and getData methods def getData(self, timestep=None): - if not self._validate(): return self._prepare1() # prepare the vfactor if timestep is None: @@ -629,7 +626,6 @@ def getData(self, timestep=None): # Iterator on UNSORTED particles for a given timestep def iterParticles(self, timestep, chunksize=1): - if not self._validate(): return self._prepare1() # prepare the vfactor if timestep not in self._timesteps: @@ -765,7 +761,6 @@ def _animateOnAxes_2D(self, ax, t, cax_id=0): # Convert data to VTK format def toVTK(self, rendering="trajectory", data_format="xml"): - if not self._validate(): return if self._ndim_particles != 3: print ("Cannot export tracked particles of a "+str(self._ndim_particles)+"D simulation to VTK") diff --git a/happi/_Factories.py b/happi/_Factories.py index 8f47325ec..b862c48b8 100755 --- a/happi/_Factories.py +++ b/happi/_Factories.py @@ -39,10 +39,8 @@ def __init__(self, simulation, scalar=None): # If not a specific scalar (root level), build a list of scalar shortcuts if scalar is None: if simulation._verbose: print("Scanning for Scalar diagnostics") - # Create a temporary, empty scalar diagnostic - tmpDiag = Scalar(simulation) # Get a list of scalars - scalars = tmpDiag.getScalars() + scalars = simulation.getScalars() # Create scalars shortcuts for scalar in scalars: setattr(self, scalar, ScalarFactory(simulation, scalar)) @@ -53,7 +51,25 @@ def __init__(self, simulation, scalar=None): def __call__(self, *args, **kwargs): return Scalar(self._simulation, *(self._additionalArgs+args), **kwargs) - + + def __repr__(self): + msg = object.__repr__(self) + if len(self._additionalArgs) == 0 and self._simulation._scan: + scalars = self._simulation.getScalars() + if scalars: + msg += "\nAvailable Scalar diagnostics:\n" + l = [""] + for s in scalars: + if len(s)>4 and s[:2]!=l[-1][:2] and s[-2:]!=l[-1][-2:]: + if l!=[""]: + msg += "\t".join(l) + "\n" + l = [] + l.append(s) + if l!=[""]: + msg += "\t".join(l) + "\n" + else: + msg += "\nNo Scalar diagnostics available" + return msg class FieldFactory(object): """Import and analyze a Field diagnostic from a Smilei simulation diff --git a/happi/_SmileiSimulation.py b/happi/_SmileiSimulation.py index b8c3a3fbe..279f9f512 100644 --- a/happi/_SmileiSimulation.py +++ b/happi/_SmileiSimulation.py @@ -281,6 +281,35 @@ def _getParticleListSpecies(self, filePrefix): if s: species += [ s.groups()[0] ] return list(set(species)) # unique species + # get all available scalars + def getScalars(self): + allScalars = None + for path in self._results_path: + try: + file = path+'/scalars.txt' + f = open(file, 'r') + except Exception as e: + continue + try: + # Find last commented line + prevline = "" + for line in f: + line = line.strip() + if line[0]!="#": break + prevline = line[1:].strip() + scalars = str(prevline).split() # list of scalars + scalars = scalars[1:] # remove first, which is "time" + except Exception as e: + scalars = [] + f.close() + if allScalars is None: + allScalars = scalars + else: + allScalars = self._np.intersect1d(allScalars, scalars) + if allScalars is None: + return [] + return allScalars + def getTrackSpecies(self): """ List the available tracked species """ return self._getParticleListSpecies("TrackParticlesDisordered") diff --git a/happi/_Utils.py b/happi/_Utils.py index 12bee2315..9fd35a757 100755 --- a/happi/_Utils.py +++ b/happi/_Utils.py @@ -195,7 +195,7 @@ def set(self, **kwargs): kw = kwa[:-5] self.labels_font[kw] = val elif kwa in ["xticklabels", "yticklabels"]: - self.ticklabels[kw] = val + self.ticklabels[kwa] = val elif kwa in ["xticklabels_font", "yticklabels_font"]: kw = kwa[:-5] self.ticklabels_font[kw] = val @@ -331,8 +331,7 @@ class Operation(object): Parameters: ----------- - operation: the user's requested operation - pattern: the regexp pattern to find the variables in the operation + operation: the user's requested operation (string containing variables and operators) QuantityTranslator: function that takes a string as argument (a quantity name) and outputs its units + its replacement string + its displayed name ureg: Pint's unit registry or None for no unit awareness diff --git a/makefile b/makefile index 619f25967..e71133f09 100755 --- a/makefile +++ b/makefile @@ -532,9 +532,9 @@ help: @echo ' inspector : to compile for Intel Inspector analysis' @echo ' gpu_nvidia : to compile for GPU (uses OpenACC)' @echo ' gpu_amd : to compile for GPU (uses OpenMP)' - @echo ' omptasks : to compile with OpenMP tasks' - @echo ' part_event_tracing_tasks_on : to compile particle event tracing and OpenMP tasks' - @echo ' part_event_tracing_tasks_off : to compile particle event tracing without OpenMP tasks' +# @echo ' omptasks : to compile with OpenMP tasks' +# @echo ' part_event_tracing_tasks_on : to compile particle event tracing and OpenMP tasks' +# @echo ' part_event_tracing_tasks_off : to compile particle event tracing without OpenMP tasks' # @echo ' omptasks : to compile with OpenMP tasks' # @echo ' part_event_tracing_tasks_on : to compile particle event tracing and OpenMP tasks' # @echo ' part_event_tracing_tasks_off : to compile particle event tracing without OpenMP tasks' diff --git a/scripts/compile_tools/machine/joliot_curie_rome b/scripts/compile_tools/machine/joliot_curie_rome index d77047de0..c2974e380 100755 --- a/scripts/compile_tools/machine/joliot_curie_rome +++ b/scripts/compile_tools/machine/joliot_curie_rome @@ -4,12 +4,12 @@ # # Load the correct modules: # - +#module load mpi #module load hdf5 #module switch flavor/hdf5/serial flavor/hdf5/parallel #export HDF5_ROOT_DIR=${HDF5_ROOT} #export OMPI_MCA_btl_portals4_use_rdma=0 -#module load python3/3.7.5 +#module load python3/3.8.10 #export PYTHONEXE=python3 #export PYTHONHOME=$PYTHON3_ROOT #export PATH=$PYTHONHOME/bin:$PATH diff --git a/src/Diagnostic/DiagnosticFields.cpp b/src/Diagnostic/DiagnosticFields.cpp index 94ec4edaa..2f6b20444 100755 --- a/src/Diagnostic/DiagnosticFields.cpp +++ b/src/Diagnostic/DiagnosticFields.cpp @@ -154,10 +154,11 @@ DiagnosticFields::DiagnosticFields( Params ¶ms, SmileiMPI *smpi, VectorPatch for( unsigned int ipatch=0 ; ipatchEMfields->allFields_avg.resize( diag_n+1 ); if( time_average > 1 ) { - for( unsigned int ifield=0; ifieldEMfields->allFields_avg[diag_n].push_back( vecPatches( ipatch )->EMfields->createField( fields_names[ifield], params ) ); + } } } } @@ -333,9 +334,16 @@ void DiagnosticFields::run( SmileiMPI *smpi, VectorPatch &vecPatches, int itime, // Write H5Write dset = writeField( iteration_group_, fields_names[ifield] ); // Attributes for openPMD + Field *f = vecPatches( 0 )->EMfields->allFields[fields_indexes[ifield]]; + vector stagger( f->dims().size() ); + for( unsigned int i = 0; i < stagger.size(); i++ ) { + stagger[i] = 0.5 * (double) f->isDual(i); + } + bool ends_with_m = 0 == fields_names[ifield].compare( fields_names[ifield].length()-2, 2, "_m" ); + double stagger_t = ends_with_m ? vecPatches( 0 )->EMfields->timestep*0.5 : 0.; openPMD_->writeFieldAttributes( dset, subgrid_start_, subgrid_step_ ); - openPMD_->writeRecordAttributes( dset, field_type[ifield] ); - openPMD_->writeFieldRecordAttributes( dset ); + openPMD_->writeRecordAttributes( dset, field_type[ifield], stagger_t ); + openPMD_->writeFieldRecordAttributes( dset, stagger ); openPMD_->writeComponentAttributes( dset, field_type[ifield] ); } #pragma omp barrier diff --git a/src/Diagnostic/DiagnosticParticleBinningBase.cpp b/src/Diagnostic/DiagnosticParticleBinningBase.cpp index 33288631d..421956ccf 100755 --- a/src/Diagnostic/DiagnosticParticleBinningBase.cpp +++ b/src/Diagnostic/DiagnosticParticleBinningBase.cpp @@ -246,8 +246,8 @@ void DiagnosticParticleBinningBase::calculate_auto_limits( Patch *patch, SimWind if( !std::isnan( axis->min ) && !std::isnan( axis->max ) ) { continue; } - double axis_min = numeric_limits::max(); - double axis_max = numeric_limits::lowest(); + double axis_min = std::isnan( axis->min ) ? numeric_limits::max() : numeric_limits::lowest(); + double axis_max = std::isnan( axis->max ) ? numeric_limits::lowest() : numeric_limits::max(); for( unsigned int i_s=0; i_svecSpecies[species_indices[i_s]]; unsigned int n = s->getNbrOfParticles(); diff --git a/src/Diagnostic/Histogram.cpp b/src/Diagnostic/Histogram.cpp index 0da6251b4..6090e171d 100755 --- a/src/Diagnostic/Histogram.cpp +++ b/src/Diagnostic/Histogram.cpp @@ -71,7 +71,7 @@ void Histogram::digitize( vector species, } } else { // if the particles out of the "box" must be included - + for( unsigned int ipart = 0 ; ipart < npart ; ipart++ ) { // skip already discarded particles if( int_buffer[ipart] < 0 ) { @@ -82,13 +82,12 @@ void Histogram::digitize( vector species, // move out-of-range indexes back into range if( ind < 0 ) { ind = 0; - } - if( ind >= axis->nbins ) { + } else if( ind >= axis->nbins ) { ind = axis->nbins-1; } int_buffer[ipart] += ind; } - + } } // loop axes diff --git a/src/ElectroMagn/ElectroMagn.h b/src/ElectroMagn/ElectroMagn.h index f21bf53d6..b5bd3e5bc 100755 --- a/src/ElectroMagn/ElectroMagn.h +++ b/src/ElectroMagn/ElectroMagn.h @@ -305,7 +305,6 @@ class ElectroMagn //! Number of elements in arrays without duplicated borders unsigned int bufsize[3][2]; - //!\todo should this be just an integer??? //! Oversize domain to exchange less particles (from params) std::vector oversize; @@ -380,9 +379,6 @@ class ElectroMagn cField *p_AM_; cField *Ap_AM_; - //! \todo check time_dual or time_prim (MG) -// //! method used to solve Maxwell's equation (takes current time and time-step as input parameter) -// virtual void solveMaxwellAmpere() = 0; //! Maxwell Ampere Solver Solver *MaxwellAmpereSolver_; //! Maxwell Faraday Solver diff --git a/src/ElectroMagnBC/ElectroMagnBCAM_BM.cpp b/src/ElectroMagnBC/ElectroMagnBCAM_BM.cpp index 29252c4fe..af72f9906 100755 --- a/src/ElectroMagnBC/ElectroMagnBCAM_BM.cpp +++ b/src/ElectroMagnBC/ElectroMagnBCAM_BM.cpp @@ -188,10 +188,6 @@ void ElectroMagnBCAM_BM::apply( ElectroMagn *EMfields, double, Patch *patch ) unsigned int j= n_d[1]-2; - - // MESSAGE("JGLOB "<< patch->getCellStartingGlobalIndex(1)+j); - //std::cout<<"come heree "<getCellStartingGlobalIndex(1)<<" "<1.){ - //MESSAGE("BlBM"); - //MESSAGE(i); - //MESSAGE(j+1); - //MESSAGE((*Bl)(i,j+1)); - //} }//i ---end Bl // for Bt^(d,d) @@ -215,12 +205,6 @@ void ElectroMagnBCAM_BM::apply( ElectroMagn *EMfields, double, Patch *patch ) + Gamma_Bt_Rmax * ( *Bt_old )( i, j ) - Icpx * ( double )imode * CB_BM * Epsilon_Bt_Rmax * ( ( *Br )( i, j ) + ( *Br_old )( i, j ) ) - CE_BM * Delta_Bt_Rmax * ( ( *Er )( i, j+1 )+( *Er )( i, j )-( *Er )( i-1, j+1 ) -( *Er )( i-1, j ) ) ; - //if (std::abs((*Bt)(i,j+1))>1.){ - // MESSAGE("BtMF"); - // MESSAGE(i); - // MESSAGE(j+1); - // MESSAGE((*Bt)(i,j+1)); - //} }//i ---end Bt } } diff --git a/src/ElectroMagnBC/ElectroMagnBCAM_SM.cpp b/src/ElectroMagnBC/ElectroMagnBCAM_SM.cpp index 28f4b42b7..dd4a790f8 100755 --- a/src/ElectroMagnBC/ElectroMagnBCAM_SM.cpp +++ b/src/ElectroMagnBC/ElectroMagnBCAM_SM.cpp @@ -203,8 +203,8 @@ void ElectroMagnBCAM_SM::apply( ElectroMagn *EMfields, double time_dual, Patch * unsigned int i=0; ( *Br )( i, j ) = Alpha_Xmin * ( *Et )( i, j ) + Beta_Xmin * ( *Br )( i+1, j ) - + Gamma_Xmin * byW; - + Delta_Xmin *( ( *Bl )( i, j+1 )- ( *Bl )( i, j ) ); + + Gamma_Xmin * byW + + Delta_Xmin *( ( *Bl )( i, j+1 )- ( *Bl )( i, j ) ); }//j ---end compute Br @@ -258,7 +258,7 @@ void ElectroMagnBCAM_SM::apply( ElectroMagn *EMfields, double time_dual, Patch * ( *Br )( i, j ) = - Alpha_Xmax * ( *Et )( i-1, j ) + Beta_Xmax * ( *Br )( i-1, j ) + Gamma_Xmax * byE - + Delta_Xmax * ( ( *Bl )( i-1, j+1 )- ( *Bl )( i-1, j ) ); // Check x-index + - Delta_Xmax * ( ( *Bl )( i-1, j+1 )- ( *Bl )( i-1, j ) ); // Check x-index }//j ---end compute Br diff --git a/src/ElectroMagnSolver/MA_SolverAM_norm.cpp b/src/ElectroMagnSolver/MA_SolverAM_norm.cpp index 11d239a53..6c5633fcd 100755 --- a/src/ElectroMagnSolver/MA_SolverAM_norm.cpp +++ b/src/ElectroMagnSolver/MA_SolverAM_norm.cpp @@ -79,7 +79,7 @@ void MA_SolverAM_norm::operator()( ElectroMagn *fields ) ( *El )( i, j-1 )=-( *El )( i, j+1 ); } for( unsigned int i=0 ; i 1 diff --git a/src/Field/Field.h b/src/Field/Field.h index aecd63a15..669106245 100755 --- a/src/Field/Field.h +++ b/src/Field/Field.h @@ -86,7 +86,6 @@ class Field virtual void shift_x( unsigned int delta ) = 0; //! vector containing the dimensions of the Field - //! \todo private/friend/modify (JD) std::vector dims_; //! keep track ofwich direction of the Field is dual diff --git a/src/Params/OpenPMDparams.cpp b/src/Params/OpenPMDparams.cpp index b7f4405e8..001836e9e 100755 --- a/src/Params/OpenPMDparams.cpp +++ b/src/Params/OpenPMDparams.cpp @@ -18,16 +18,12 @@ OpenPMDparams::OpenPMDparams( Params &p ): // Grid parameters string xyz = "xyz"; - gridGlobalOffset.resize( params->nDim_field ); - gridOffset .resize( params->nDim_field ); - position .resize( params->nDim_field ); + gridGlobalOffset.resize( params->nDim_field, 0. ); + gridOffset .resize( params->nDim_field, 0. ); gridSpacing .resize( params->nDim_field ); for( unsigned int idim=0; idimnDim_field; idim++ ) { axisLabels.addString( xyz.substr( idim, 1 ) ); - gridGlobalOffset[idim] = 0.; - gridOffset [idim] = 0.; - position [idim] = 0.; - gridSpacing [idim] = params->cell_length[idim]; + gridSpacing[idim] = params->cell_length[idim]; } fieldSolverParameters = ""; if( params->maxwell_sol == "Yee" ) { @@ -192,15 +188,15 @@ void OpenPMDparams::writeSpeciesAttributes( H5Write & ) { } -void OpenPMDparams::writeRecordAttributes( H5Write &location, unsigned int unit_type ) +void OpenPMDparams::writeRecordAttributes( H5Write &location, unsigned int unit_type, double timeOffset ) { location.attr( "unitDimension", unitDimension[unit_type] ); - location.attr( "timeOffset", 0. ); + location.attr( "timeOffset", timeOffset ); } -void OpenPMDparams::writeFieldRecordAttributes( H5Write &location ) +void OpenPMDparams::writeFieldRecordAttributes( H5Write &location, vector &stagger ) { - location.attr( "position", position ); + location.attr( "position", stagger ); } void OpenPMDparams::writeComponentAttributes( H5Write &location, unsigned int unit_type ) diff --git a/src/Params/OpenPMDparams.h b/src/Params/OpenPMDparams.h index 8971360a6..3a3216ea1 100755 --- a/src/Params/OpenPMDparams.h +++ b/src/Params/OpenPMDparams.h @@ -49,8 +49,6 @@ class OpenPMDparams DividedString fieldBoundary, fieldBoundaryParameters, particleBoundary, particleBoundaryParameters; //! current smoothing description std::string currentSmoothing, currentSmoothingParameters; - //! position of the fields within a cell - std::vector position; //! Returns a time string in the openPMD format std::string getLocalTime(); @@ -74,10 +72,10 @@ class OpenPMDparams void writeSpeciesAttributes( H5Write& ); //! Write the attributes for a record - void writeRecordAttributes( H5Write&, unsigned int ); + void writeRecordAttributes( H5Write&, unsigned int, double timeOffset = 0 ); //! Write the attributes for a field record - void writeFieldRecordAttributes( H5Write& ); + void writeFieldRecordAttributes( H5Write&, std::vector &stagger ); //! Write the attributes for a component void writeComponentAttributes( H5Write&, unsigned int ); diff --git a/src/PartCompTime/PartCompTimeFactory.h b/src/PartCompTime/PartCompTimeFactory.h index fccb45038..e8fd5478f 100644 --- a/src/PartCompTime/PartCompTimeFactory.h +++ b/src/PartCompTime/PartCompTimeFactory.h @@ -70,7 +70,7 @@ class PartCompTimeFactory } } else { - ERROR_NAMELIST( "Unknwon parameters setting the particle computing time evaluation operator: " + ERROR_NAMELIST( "Unknown parameters setting the particle computing time evaluation operator: " << params.geometry << ", Order : " << params.interpolation_order, LINK_NAMELIST + std::string("#main-variables")); } diff --git a/src/Particles/ParticleCreator.cpp b/src/Particles/ParticleCreator.cpp index 05a26fb57..270d0fee1 100644 --- a/src/Particles/ParticleCreator.cpp +++ b/src/Particles/ParticleCreator.cpp @@ -285,13 +285,13 @@ int ParticleCreator::create( struct SubSpace sub_space, } // In AM, normalization of weights might be required - bool renormalize = true; + bool regular_weight = false; if( position_initialization_ == "regular" ) { - renormalize = false; + regular_weight = true; } else if( position_initialization_on_species_ ) { unsigned int ispec = species_->position_initialization_on_species_index_; if( patch->vecSpecies[ispec]->position_initialization_ == "regular" ) { - renormalize = false; + regular_weight = true; } } @@ -328,7 +328,7 @@ int ParticleCreator::create( struct SubSpace sub_space, ParticleCreator::createPosition( position_initialization_, regular_number_array_, particles_, species_, nPart, iPart, indexes, params, patch->rand_ ); } ParticleCreator::createMomentum( momentum_initialization_, particles_, species_, nPart, iPart, &temp[0], &vel[0], patch->rand_ ); - ParticleCreator::createWeight( particles_, nPart, iPart, density( i, j, k ), params, renormalize ); + ParticleCreator::createWeight( particles_, nPart, iPart, density( i, j, k ), params, regular_weight ); ParticleCreator::createCharge( particles_, species_, nPart, iPart, charge( i, j, k ) ); iPart += nPart; @@ -931,33 +931,28 @@ void ParticleCreator::createMomentum( std::string momentum_initialization, //! For all (nPart) particles in a mesh initialize its numerical weight (equivalent to a number density) // --------------------------------------------------------------------------------------------------------------------- void ParticleCreator::createWeight( Particles * particles, unsigned int nPart, unsigned int iPart, double n_real_particles, - Params ¶ms, bool renormalize ) + Params ¶ms, bool regular_weight ) { double w = n_real_particles / nPart; for( unsigned int p=iPart; pweight( p ) = w ; } - // In AM, we have a correction to make : multiply by radius + // In AM, we have a correction to make : multiply by cell radius // because n_real_particles was computed with the cell section, not cell volume // See above : density( i, j, k ) *= params.cell_volume; - // where params.cell_volume is 2*pi*dR*dL (in AM only) + // where params.cell_volume is 2*pi*dR*dL (in AM only). Multiplication by R is missing. if( params.geometry == "AMcylindrical" ) { - double radius; - for( unsigned int p=iPart; pposition(1,p)*particles->position(1,p) + particles->position(2,p)*particles->position(2,p)); - particles->weight( p ) *= radius; - } - // We also need to renormalize in case total weight is not exaclty the same anymore - if( renormalize ) { - double total_weight = 0.; + if (regular_weight){ // Multiply each weight by the radius of the particle. Particles towards the axis are lighter. for( unsigned int p=iPart; pweight( p ); + double radius = sqrt(particles->position(1,p)*particles->position(1,p) + particles->position(2,p)*particles->position(2,p)); + particles->weight( p ) *= radius; } + } else { // Multiply each weight by the radius of the cell. All particles of the cell have identical weights. + double radius = sqrt(particles->position(1,iPart)*particles->position(1,iPart) + particles->position(2,iPart)*particles->position(2,iPart)); double cell_radius = params.cell_length[1] * (floor(radius/params.cell_length[1]) + 0.5); - double coeff = n_real_particles * cell_radius / total_weight; for( unsigned int p=iPart; pweight( p ) *= coeff; + particles->weight( p ) *= cell_radius; } } } diff --git a/src/Particles/ParticleCreator.h b/src/Particles/ParticleCreator.h index 336e98ad6..544acda65 100644 --- a/src/Particles/ParticleCreator.h +++ b/src/Particles/ParticleCreator.h @@ -81,7 +81,7 @@ class ParticleCreator unsigned int iPart, double n_real_particles, Params ¶ms, - bool renormalize ); + bool regular_weight ); // For all particles in a mesh initialize its charge state static void createCharge( Particles * particles, Species * species, diff --git a/src/Patch/VectorPatch.cpp b/src/Patch/VectorPatch.cpp index e7d26bedf..0c2fbb036 100755 --- a/src/Patch/VectorPatch.cpp +++ b/src/Patch/VectorPatch.cpp @@ -339,10 +339,6 @@ void VectorPatch::dynamics( Params ¶ms, timers.particles.restart(); ostringstream t; -#ifdef _PARTEVENTTRACING - bool diag_PartEventTracing {false}; - double reference_time; -#endif #ifdef _OMPTASKS #pragma omp single @@ -362,7 +358,7 @@ void VectorPatch::dynamics( Params ¶ms, bool diag_PartEventTracing {false}; if( !params.Laser_Envelope_model ) { diag_PartEventTracing = smpi->diagPartEventTracing( time_dual, params.timestep); - if (diag_PartEventTracing) smpi->reference_time = MPI_Wtime(); + if (diag_PartEventTracing) smpi->reference_time_ = MPI_Wtime(); } #endif @@ -1406,6 +1402,10 @@ void VectorPatch::runAllDiags( Params &/*params*/, SmileiMPI *smpi, unsigned int for( unsigned int i=0; ipatches_mins[ipatch][i] ); + } + } + for( unsigned int i=0; ipatches_maxs[ipatch][i] ); } } @@ -4559,7 +4559,7 @@ void VectorPatch::ponderomotiveUpdateSusceptibilityAndMomentum( Params ¶ms, # ifdef _PARTEVENTTRACING bool diag_PartEventTracing = smpi->diagPartEventTracing( time_dual, params.timestep); - if (diag_PartEventTracing) smpi->reference_time = MPI_Wtime(); + if (diag_PartEventTracing) smpi->reference_time_ = MPI_Wtime(); # endif @@ -4601,7 +4601,7 @@ void VectorPatch::ponderomotiveUpdatePositionAndCurrents( Params ¶ms, #ifdef _PARTEVENTTRACING bool diag_PartEventTracing {false}; diag_PartEventTracing = smpi->diagPartEventTracing( time_dual, params.timestep); - // if (diag_PartEventTracing) smpi->reference_time = MPI_Wtime(); + // if (diag_PartEventTracing) smpi->reference_time_ = MPI_Wtime(); #endif #ifdef _OMPTASKS diff --git a/src/Projector/ProjectorAM2Order.cpp b/src/Projector/ProjectorAM2Order.cpp index caaeff63a..5487995df 100755 --- a/src/Projector/ProjectorAM2Order.cpp +++ b/src/Projector/ProjectorAM2Order.cpp @@ -77,7 +77,7 @@ void ProjectorAM2Order::currents( ElectroMagnAM *emAM, // arrays used for the Esirkepov projection method double Sl0[5], Sl1[5], Sr0[5], Sr1[5], DSl[5], DSr[5]; complex Jl_p[5], Jr_p[5]; - complex e_delta, e_delta_m1, e_delta_inv, e_bar, e_bar_m1, C_m = 1.; //, C_m_old; + complex e_delta, e_delta_m1, e_bar, e_bar_m1, C_m = 1.; //, C_m_old; complex *Jl, *Jr, *Jt, *rho; for( unsigned int i=0; i<5; i++ ) { @@ -107,7 +107,7 @@ void ProjectorAM2Order::currents( ElectroMagnAM *emAM, //calculate exponential coefficients double rp = sqrt( particles.position( 1, ipart )*particles.position( 1, ipart )+particles.position( 2, ipart )*particles.position( 2, ipart ) ); - std::complex theta_old = array_eitheta_old[0]; + std::complex eitheta_old = array_eitheta_old[0]; std::complex eitheta = ( particles.position( 1, ipart ) + Icpx * particles.position( 2, ipart ) ) / rp ; //exp(i theta) e_bar = 1.; // locate the particle on the primal grid at current time-step & calculate coeff. S1 @@ -138,8 +138,9 @@ void ProjectorAM2Order::currents( ElectroMagnAM *emAM, double r_bar = ((jpo + j_domain_begin_ + deltaold[1*nparts])*dr + rp) * 0.5; // r at t = t0 - dt/2 - e_delta_m1 = std::sqrt(eitheta * (2.*std::real(theta_old) - theta_old)); // std::sqrt keeps the root with positive real part which is what we need here. - e_bar_m1 = theta_old * e_delta_m1; + //e_delta_m1 = std::sqrt(eitheta * (2.*std::real(theta_old) - theta_old)); // std::sqrt keeps the root with positive real part which is what we need here. + e_delta_m1 = std::sqrt(eitheta * std::conj(eitheta_old)); //ei^(theta-theta_old)/2. std::sqrt keeps the root with positive real part which is what we need here. + e_bar_m1 = eitheta_old * e_delta_m1; //ei^(theta+theta_old)/2. ipo -= 2; //This minus 2 come from the order 2 scheme, based on a 5 points stencil from -2 to +2. // i/j/kpo stored with - i/j/k_domain_begin in Interpolator @@ -183,9 +184,9 @@ void ProjectorAM2Order::currents( ElectroMagnAM *emAM, } e_delta = 0.5; - e_delta_inv = 0.5; + //e_delta_inv = 0.5; - //Compute division by R in advance for Jt and rho evaluation. + //Compute division by R in advance for Jt and rho evaluation. for( unsigned int j=0 ; j<5 ; j++ ) { Sr0[j] *= invR__local[j]; Sr1[j] *= invR__local[j]; @@ -197,8 +198,9 @@ void ProjectorAM2Order::currents( ElectroMagnAM *emAM, e_delta *= e_delta_m1; e_bar *= e_bar_m1; C_m = 2. * e_bar ; //multiply modes > 0 by 2 and C_m = 1 otherwise. - e_delta_inv =1./e_delta - 1.; + //e_delta_inv =1./e_delta - 1.; crt_p = charge_weight*Icpx*e_bar / ( dt*( double )imode )*2.*r_bar; + //crt_p = charge_weight*Icpx*e_bar / ( dt*( double )imode )*2.; } // Add contribution J_p to global array @@ -247,7 +249,8 @@ void ProjectorAM2Order::currents( ElectroMagnAM *emAM, iloc = ( i+ipo )*nprimr_ + jpo; for( unsigned int j=0 ; j<5 ; j++ ) { linindex = iloc+j; - Jt [linindex] += crt_p*(Sr1[j]*Sl1[i]*e_delta_inv - Sr0[j]*Sl0[i]*( e_delta-1. )); + //Jt [linindex] += crt_p*(Sr1[j]*Sl1[i]*e_delta_inv - Sr0[j]*Sl0[i]*( e_delta-1. )); + Jt [linindex] -= crt_p*(Sr1[j]*Sl1[i]*( e_delta-1. ) - Sr0[j]*Sl0[i]*(std::conj(e_delta) - 1.)); } } @@ -354,8 +357,6 @@ void ProjectorAM2Order::basicForComplexOnBuffer( complex *rhoj, Particle // This function also assumes that particles position is evaluated at the same time as currents which is usually not true (half time-step difference). // It will therefore fail to evaluate the current accurately at t=0 if a plasma is already in the box. - - // ------------------------------------- // Variable declaration & initialization // ------------------------------------- @@ -366,17 +367,6 @@ void ProjectorAM2Order::basicForComplexOnBuffer( complex *rhoj, Particle if( type > 0 ) { //if current density ERROR("This projector can be used only for charge density at the moment."); - // charge_weight *= 1./sqrt( 1.0 + particles.momentum( 0, ipart )*particles.momentum( 0, ipart ) - // + particles.momentum( 1, ipart )*particles.momentum( 1, ipart ) - // + particles.momentum( 2, ipart )*particles.momentum( 2, ipart ) ); - // if( type == 1 ) { //if Jl - // charge_weight *= particles.momentum( 0, ipart ); - // } else if( type == 2 ) { //if Jr - // charge_weight *= ( particles.momentum( 1, ipart )*particles.position( 1, ipart ) + particles.momentum( 2, ipart )*particles.position( 2, ipart ) )/ r ; - // nr++; - // } else { //if Jt - // charge_weight *= ( -particles.momentum( 1, ipart )*particles.position( 2, ipart ) + particles.momentum( 2, ipart )*particles.position( 1, ipart ) ) / r ; - // } } complex e_theta = ( particles.position( 1, ipart ) + Icpx*particles.position( 2, ipart ) )/r; @@ -460,8 +450,9 @@ void ProjectorAM2Order::axisBC(ElectroMagnAM *emAM, bool diag_flag ) void ProjectorAM2Order::apply_axisBC(std::complex *rhoj,std::complex *Jl, std::complex *Jr, std::complex *Jt, unsigned int imode, bool diag_flag ) { - double sign = 1.; - for (unsigned int i=0; i< imode; i++) sign *= -1; + // Mode 0 contribution "below axis" is added. + // Non zero modes are substracted because a particle sitting exactly on axis has a non defined theta and can not contribute to a theta dependent mode. + double sign = (imode == 0) ? 1 : -1 ; if (diag_flag && rhoj) { for( unsigned int i=2 ; i *rhoj,std::complex 0){ rhoj[i] = 0.; - rhoj[i-1] = - rhoj[i+1]; + rhoj[i-1] = - rhoj[i+1]; // Zero Jl mode > 0 on axis. } else { - //This is just for cosmetics on the picture, rho has no influence on the results - rhoj[i] = (4.*rhoj[i+1] - rhoj[i+2])/3.; - rhoj[i-1] = rhoj[i+1]; + rhoj[i] = rhoj[i+1]; //This smoothing is just for cosmetics on the picture, rho has no influence on the results. + rhoj[i-1] = rhoj[i+1]; // Non zero Jl mode > 0 on axis. } } } @@ -489,10 +479,11 @@ void ProjectorAM2Order::apply_axisBC(std::complex *rhoj,std::complex 0){ Jl [i] = 0. ; - Jl[i-1] = -Jl[i+1]; + Jl[i-1] = -Jl[i+1]; // Zero Jl mode > 0 on axis. } else { - //Jl mode 1 on axis is left as is. It looks over estimated but it is necessary to conserve a correct divergence and a proper evaluation on the field on axis. - Jl [i-1] = Jl [i+1] ; + //Jl mode 0 on axis should be left as is. It looks over estimated but it might be necessary to conserve a correct divergence and a proper evaluation on the field on axis. + //Jl [i] = (4.*Jl [i+1] - Jl [i+2])/3. ; //Here we use smoothig on axis for results appearing more physical. + Jl [i-1] = Jl [i+1] ; // Non zero Jl mode 0 on axis. } } } @@ -503,15 +494,15 @@ void ProjectorAM2Order::apply_axisBC(std::complex *rhoj,std::complex