diff --git a/.github/workflows/cuda.yml b/.github/workflows/cuda.yml index 55d5ab171ad..3a14aa81b44 100644 --- a/.github/workflows/cuda.yml +++ b/.github/workflows/cuda.yml @@ -110,7 +110,7 @@ jobs: which nvcc || echo "nvcc not in PATH!" git clone https://github.com/AMReX-Codes/amrex.git ../amrex - cd amrex && git checkout --detach 22.09 && cd - + cd amrex && git checkout --detach 13aa4df0f5a4af40270963ad5b42ac7ce662e045 && cd - make COMP=gcc QED=FALSE USE_MPI=TRUE USE_GPU=TRUE USE_OMP=FALSE USE_PSATD=TRUE USE_CCACHE=TRUE -j 2 build_nvhpc21-11-nvcc: diff --git a/.zenodo.json b/.zenodo.json index bd831642e16..6aa723ea0e3 100644 --- a/.zenodo.json +++ b/.zenodo.json @@ -38,6 +38,11 @@ "name": "Fedeli, Luca", "orcid": "0000-0002-7215-4178" }, + { + "affiliation": "Lawrence Berkeley National Laboratory", + "name": "Garten, Marco", + "orcid": "0000-0001-6994-2475" + }, { "affiliation": "SLAC National Accelerator Laboratory", "name": "Ge, Lixin", @@ -130,7 +135,7 @@ }, { "affiliation": "Lawrence Berkeley National Laboratory", - "name": "Rheaume, Tiberius", + "name": "Rheaume, Elisa", "orcid": "0000-0002-6710-0650" }, { diff --git a/CMakeLists.txt b/CMakeLists.txt index f34e41e539c..b6d6f5a714b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,7 +1,7 @@ # Preamble #################################################################### # cmake_minimum_required(VERSION 3.20.0) -project(WarpX VERSION 22.09) +project(WarpX VERSION 22.10) include(${WarpX_SOURCE_DIR}/cmake/WarpXFunctions.cmake) @@ -262,6 +262,12 @@ if(WarpX_PSATD) if(WarpX_DIMS STREQUAL RZ) target_link_libraries(ablastr PUBLIC blaspp) target_link_libraries(ablastr PUBLIC lapackpp) + + # BLAS++ forgets to declare cuBLAS and cudaRT dependencies + if(WarpX_COMPUTE STREQUAL CUDA) + find_package(CUDAToolkit REQUIRED) + target_link_libraries(ablastr PUBLIC CUDA::cudart CUDA::cublas) + endif() endif() endif() diff --git a/Docs/requirements.txt b/Docs/requirements.txt index 1eedc4e08f2..3da427f88b4 100644 --- a/Docs/requirements.txt +++ b/Docs/requirements.txt @@ -4,6 +4,8 @@ # # License: BSD-3-Clause-LBNL +# WarpX PICMI bindings w/o C++ component (used for autoclass docs) +-e ../Python breathe # docutils 0.17 breaks HTML tags & RTD theme # https://github.com/sphinx-doc/sphinx/issues/9001 @@ -11,10 +13,13 @@ docutils<=0.16 # PICMI API docs # note: keep in sync with version in ../requirements.txt -picmistandard==0.0.19 +picmistandard==0.0.20 # for development against an unreleased PICMI version, use: -#picmistandard @ git+https://github.com/picmi-standard/picmi.git#subdirectory=PICMI_Python +# picmistandard @ git+https://github.com/picmi-standard/picmi.git#subdirectory=PICMI_Python + pygments recommonmark sphinx>=2.0 +sphinx-design sphinx_rtd_theme>=0.3.1 +sphinxcontrib-napoleon diff --git a/Docs/source/conf.py b/Docs/source/conf.py index c8944627df4..ac933290595 100644 --- a/Docs/source/conf.py +++ b/Docs/source/conf.py @@ -45,7 +45,9 @@ # ones. extensions = ['sphinx.ext.autodoc', 'sphinx.ext.mathjax', + 'sphinx.ext.napoleon', 'sphinx.ext.viewcode', + 'sphinx_design', 'breathe' ] @@ -71,16 +73,16 @@ # built documents. # # The short X.Y version. -version = u'22.09' +version = u'22.10' # The full version, including alpha/beta/rc tags. -release = u'22.09' +release = u'22.10' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. # # This is also used if you do content translation via gettext catalogs. # Usually you set "language" from the command line for these cases. -language = None +language = 'en' # List of patterns, relative to source directory, that match files and # directories to ignore when looking for source files. diff --git a/Docs/source/developers/documentation.rst b/Docs/source/developers/documentation.rst index e8cad8e1b7d..f7929eb9dd1 100644 --- a/Docs/source/developers/documentation.rst +++ b/Docs/source/developers/documentation.rst @@ -41,9 +41,11 @@ Breathe documentation --------------------- Your Doxygen documentation is not only useful for people looking into the code, it is also part of the `WarpX online documentation `_ based on `Sphinx `_! -This is done using the Python module `Breathe `_, that allows you to read Doxygen documentation dorectly in the source and include it in your Sphinx documentation, by calling Breathe functions. +This is done using the Python module `Breathe `_, that allows you to write Doxygen documentation directly in the source and have it included it in your Sphinx documentation, by calling Breathe functions. For instance, the following line will get the Doxygen documentation for ``WarpXParticleContainer`` in ``Source/Particles/WarpXParticleContainer.H`` and include it to the html page generated by Sphinx: +.. code-block:: rst + .. doxygenclass:: WarpXParticleContainer Building the documentation diff --git a/Docs/source/developers/faq.rst b/Docs/source/developers/faq.rst index 273deed4c3a..d539c99da28 100644 --- a/Docs/source/developers/faq.rst +++ b/Docs/source/developers/faq.rst @@ -70,3 +70,20 @@ What does const int ``/*i_buffer*/`` mean in argument list? This is often seen in a derived class, overwriting an interface method. It means we do not name the parameter because we do not use it when we overwrite the interface. But we add the name as a comment ``/* ... */`` so that we know what we ignored when looking at the definition of the overwritten method. + + +What is Pinned Memory? +---------------------- + +We need pinned aka "page locked" host memory when we: + +- do asynchronous copies between the host and device +- want to write to CPU memory from a GPU kernel + +A typical use case is initialization of our (filtered/processed) output routines. +AMReX provides pinned memory via the ``amrex::PinnedArenaAllocator`` , which is the last argument passed to constructors of ``ParticleContainer`` and ``MultiFab``. + +Read more on this here: `How to Optimize Data Transfers in CUDA C/C++ `__ (note that pinned memory is a host memory feature and works with all GPU vendors we support) + +Bonus: underneath the hood, asynchronous MPI communications also pin and unpin memory. +One of the benefits of GPU-aware MPI implementations is, besides the possibility to use direct device-device transfers, that MPI and GPU API calls `are aware of each others' pinning ambitions `__ and do not create `data races to unpin the same memory `__. diff --git a/Docs/source/developers/particles.rst b/Docs/source/developers/particles.rst index ef5ae3f13f4..988313079d4 100644 --- a/Docs/source/developers/particles.rst +++ b/Docs/source/developers/particles.rst @@ -85,11 +85,9 @@ On a loop over particles it can be useful to access the fields on the box we are Main functions -------------- -.. doxygenfunction:: PhysicalParticleContainer::FieldGather - .. doxygenfunction:: PhysicalParticleContainer::PushPX -.. doxygenfunction:: WarpXParticleContainer::DepositCurrent +.. doxygenfunction:: WarpXParticleContainer::DepositCurrent(amrex::Vector, 3>> &J, const amrex::Real dt, const amrex::Real relative_time) .. note:: The current deposition is used both by ``PhysicalParticleContainer`` and ``LaserParticleContainer``, so it is in the parent class ``WarpXParticleContainer``. diff --git a/Docs/source/developers/profiling.rst b/Docs/source/developers/profiling.rst index cfc0d2ad4ea..69933a0b423 100644 --- a/Docs/source/developers/profiling.rst +++ b/Docs/source/developers/profiling.rst @@ -66,7 +66,7 @@ behavior of *each* individual MPI rank. The workflow for doing so is the followi cmake -S . -B build -DAMReX_BASE_PROFILE=OFF -DAMReX_TINY_PROFILE=ON - Run the simulation to be profiled. Note that the WarpX executable will create -and new folder `bl_prof`, which contains the profiling data. + a new folder `bl_prof`, which contains the profiling data. .. note:: diff --git a/Docs/source/developers/testing.rst b/Docs/source/developers/testing.rst index 76259389165..06a4cf1ccb9 100644 --- a/Docs/source/developers/testing.rst +++ b/Docs/source/developers/testing.rst @@ -25,7 +25,7 @@ For example, if you like to change the compiler to compilation to build on Nvidi branch = development cmakeSetupOpts = -DAMReX_ASSERTIONS=ON -DAMReX_TESTING=ON -DWarpX_COMPUTE=CUDA -We also support changing compilation options :ref:`via the usual build enviroment variables `__. +We also support changing compilation options via the usual :ref:`build enviroment variables `. For instance, compiling with ``clang++ -Werror`` would be: .. code-block:: sh diff --git a/Docs/source/highlights.rst b/Docs/source/highlights.rst index cc14e297177..49a7b1122ca 100644 --- a/Docs/source/highlights.rst +++ b/Docs/source/highlights.rst @@ -5,7 +5,7 @@ Science Highlights WarpX can be used in many domains of laser-plasma science, plasma physics, accelerator physics and beyond. Below, we collect a series of scientific publications that used WarpX. -Please :ref:`acknowledge WarpX in your works`, so we can find your works. +Please :ref:`acknowledge WarpX in your works `, so we can find your works. Is your publication missing? :ref:`Contact us ` or edit this page via a pull request. @@ -21,12 +21,12 @@ Scientific works in laser-plasma and beam-plasma acceleration. #. Miao B, Shrock JE, Feder L, Hollinger RC, Morrison J, Nedbailo R, Picksley A, Song H, Wang S, Rocca JJ, Milchberg HM. **Multi-GeV electron bunches from an all-optical laser wakefield accelerator**. - *preprint*. under review, 2021. - `arXiv:2112.03489 `__ + Physical Review X **12**, 031038, 2022. + `DOI:10.1103/PhysRevX.12.031038 `__ #. Mirani F, Calzolari D, Formenti A, Passoni M. **Superintense laser-driven photon activation analysis**. - Nature Communications Physics volume **4**.185, 2021 + Nature Communications Physics volume **4**.185, 2021. `DOI:10.1038/s42005-021-00685-2 `__ @@ -37,12 +37,13 @@ Scientific works in laser-ion acceleration and laser-matter interaction. #. Hakimi S, Obst-Huebl L, Huebl A, Nakamura K, Bulanov SS, Steinke S, Leemans WP, Kober Z, Ostermayr TM, Schenkel T, Gonsalves AJ, Vay J-L, Tilborg Jv, Toth C, Schroeder CB, Esarey E, Geddes CGR. **Laser-solid interaction studies enabled by the new capabilities of the iP2 BELLA PW beamline**. - under review, 2022 + Physics of Plasmas **29**, 083102, 2022. + `DOI:10.1063/5.0089331 `__ -#. Levy D, Andriyash IA, Haessler S, Ouille M, Kaur J, Flacco A, Kroupp E, Malka V, Lopez-Martens R. +#. Levy D, Andriyash IA, Haessler S, Kaur J, Ouille M, Flacco A, Kroupp E, Malka V, Lopez-Martens R. **Low-divergence MeV-class proton beams from kHz-driven laser-solid interactions**. - *preprint*. under review, 2021. - `arXiv:2112.12581 `__ + Phys. Rev. Accel. Beams **25**, 093402, 2022. + `DOI:10.1103/PhysRevAccelBeams.25.093402 `__ Particle Accelerator & Beam Physics diff --git a/Docs/source/install/hpc/crusher.rst b/Docs/source/install/hpc/crusher.rst index 34a67fa3374..a33ae0acf11 100644 --- a/Docs/source/install/hpc/crusher.rst +++ b/Docs/source/install/hpc/crusher.rst @@ -18,7 +18,7 @@ If you are new to this system, **please see the following resources**: * `Production directories `_: * ``$PROJWORK/$proj/``: shared with all members of a project, purged every 90 days (recommended) - * ``$MEMBERWORK/$proj/``: single user, purged every 90 days(usually smaller quota) + * ``$MEMBERWORK/$proj/``: single user, purged every 90 days (usually smaller quota) * ``$WORLDWORK/$proj/``: shared with all users, purged every 90 days * Note that the ``$HOME`` directory is mounted as read-only on compute nodes. That means you cannot run in your ``$HOME``. @@ -45,6 +45,21 @@ We recommend to store the above lines in a file, such as ``$HOME/crusher_warpx.p source $HOME/crusher_warpx.profile +And since Crusher does not yet provide a module for them, install BLAS++ and LAPACK++: + +.. code-block:: bash + + # BLAS++ (for PSATD+RZ) + git clone https://bitbucket.org/icl/blaspp.git src/blaspp + rm -rf src/blaspp-crusher-build + cmake -S src/blaspp -B src/blaspp-crusher-build -Duse_openmp=OFF -Dgpu_backend=hip -DCMAKE_CXX_STANDARD=17 -DCMAKE_INSTALL_PREFIX=$HOME/sw/crusher/blaspp-master + cmake --build src/blaspp-crusher-build --target install --parallel 10 + + # LAPACK++ (for PSATD+RZ) + git clone https://bitbucket.org/icl/lapackpp.git src/lapackpp + rm -rf src/lapackpp-crusher-build + cmake -S src/lapackpp -B src/lapackpp-crusher-build -DCMAKE_CXX_STANDARD=17 -Dbuild_tests=OFF -DCMAKE_INSTALL_RPATH_USE_LINK_PATH=ON -DCMAKE_INSTALL_PREFIX=$HOME/sw/crusher/lapackpp-master + cmake --build src/lapackpp-crusher-build --target install --parallel 10 Then, ``cd`` into the directory ``$HOME/src/warpx`` and use the following commands to compile: @@ -69,6 +84,8 @@ Running MI250X GPUs (2x64 GB) ^^^^^^^^^^^^^^^^^^^^^ +ECP WarpX project members, use the ``aph114`` project ID. + After requesting an interactive node with the ``getNode`` alias above, run a simulation like this, here using 8 MPI ranks and a single node: .. code-block:: bash @@ -105,4 +122,16 @@ Known System Issues .. code-block:: bash - export FI_MR_CACHE_MAX_COUNT=0 # libfabric disable caching + #export FI_MR_CACHE_MAX_COUNT=0 # libfabric disable caching + # or, less invasive: + export FI_MR_CACHE_MONITOR=memhooks # alternative cache monitor + +.. warning:: + + Sep 2nd, 2022 (OLCFDEV-1079): + rocFFT in ROCm 5.1+ tries to `write to a cache `__ in the home area by default. + This does not scale, disable it via: + + .. code-block:: bash + + export ROCFFT_RTC_CACHE_PATH=/dev/null diff --git a/Docs/source/install/hpc/frontier.rst b/Docs/source/install/hpc/frontier.rst index bf7cfbe1f4d..1263d681151 100644 --- a/Docs/source/install/hpc/frontier.rst +++ b/Docs/source/install/hpc/frontier.rst @@ -116,4 +116,16 @@ Known System Issues .. code-block:: bash - export FI_MR_CACHE_MAX_COUNT=0 # libfabric disable caching + #export FI_MR_CACHE_MAX_COUNT=0 # libfabric disable caching + # or, less invasive: + export FI_MR_CACHE_MONITOR=memhooks # alternative cache monitor + +.. warning:: + + Sep 2nd, 2022 (OLCFDEV-1079): + rocFFT in ROCm 5.1+ tries to `write to a cache `__ in the home area by default. + This does not scale, disable it via: + + .. code-block:: bash + + export ROCFFT_RTC_CACHE_PATH=/dev/null diff --git a/Docs/source/install/hpc/perlmutter.rst b/Docs/source/install/hpc/perlmutter.rst index dee1909b577..43c047c62ad 100644 --- a/Docs/source/install/hpc/perlmutter.rst +++ b/Docs/source/install/hpc/perlmutter.rst @@ -66,13 +66,13 @@ And since Perlmutter does not yet provide a module for them, install ADIOS2, BLA # BLAS++ (for PSATD+RZ) git clone https://bitbucket.org/icl/blaspp.git src/blaspp rm -rf src/blaspp-pm-build - CXX=$(which CC) cmake -S src/blaspp -B src/blaspp-pm-build -Duse_openmp=ON -Dgpu_backend=CUDA -Duse_cmake_find_blas=ON -DBLAS_LIBRARIES=${CRAY_LIBSCI_PREFIX_DIR}/lib/libsci_gnu.a -DCMAKE_CXX_STANDARD=17 -DCMAKE_INSTALL_PREFIX=$HOME/sw/perlmutter/blaspp-master + CXX=$(which CC) cmake -S src/blaspp -B src/blaspp-pm-build -Duse_openmp=OFF -Dgpu_backend=cuda -DCMAKE_CXX_STANDARD=17 -DCMAKE_INSTALL_PREFIX=$HOME/sw/perlmutter/blaspp-master cmake --build src/blaspp-pm-build --target install --parallel 16 # LAPACK++ (for PSATD+RZ) git clone https://bitbucket.org/icl/lapackpp.git src/lapackpp rm -rf src/lapackpp-pm-build - CXX=$(which CC) CXXFLAGS="-DLAPACK_FORTRAN_ADD_" cmake -S src/lapackpp -B src/lapackpp-pm-build -Duse_cmake_find_lapack=ON -DBLAS_LIBRARIES=${CRAY_LIBSCI_PREFIX_DIR}/lib/libsci_gnu.a -DLAPACK_LIBRARIES=${CRAY_LIBSCI_PREFIX_DIR}/lib/libsci_gnu.a -DCMAKE_CXX_STANDARD=17 -Dbuild_tests=OFF -DCMAKE_INSTALL_RPATH_USE_LINK_PATH=ON -DCMAKE_INSTALL_PREFIX=$HOME/sw/perlmutter/lapackpp-master + CXX=$(which CC) CXXFLAGS="-DLAPACK_FORTRAN_ADD_" cmake -S src/lapackpp -B src/lapackpp-pm-build -DCMAKE_CXX_STANDARD=17 -Dbuild_tests=OFF -DCMAKE_INSTALL_RPATH_USE_LINK_PATH=ON -DCMAKE_INSTALL_PREFIX=$HOME/sw/perlmutter/lapackpp-master cmake --build src/lapackpp-pm-build --target install --parallel 16 Optionally, download and install Python packages for :ref:`PICMI ` or dynamic ensemble optimizations (:ref:`libEnsemble `): diff --git a/Docs/source/install/hpc/summit.rst b/Docs/source/install/hpc/summit.rst index e1212804e7d..0b2cff56b6e 100644 --- a/Docs/source/install/hpc/summit.rst +++ b/Docs/source/install/hpc/summit.rst @@ -45,6 +45,22 @@ We recommend to store the above lines in a file, such as ``$HOME/summit_warpx.pr source $HOME/summit_warpx.profile +For PSATD+RZ simulations, you will need to build BLAS++ and LAPACK++: + +.. code-block:: bash + + # BLAS++ (for PSATD+RZ) + git clone https://bitbucket.org/icl/blaspp.git src/blaspp + rm -rf src/blaspp-summit-build + cmake -S src/blaspp -B src/blaspp-summit-build -Duse_openmp=OFF -Dgpu_backend=cuda -DCMAKE_CXX_STANDARD=17 -DCMAKE_INSTALL_PREFIX=$HOME/sw/summit/blaspp-master + cmake --build src/blaspp-summit-build --target install --parallel 10 + + # LAPACK++ (for PSATD+RZ) + git clone https://bitbucket.org/icl/lapackpp.git src/lapackpp + rm -rf src/lapackpp-summit-build + cmake -S src/lapackpp -B src/lapackpp-summit-build -DCMAKE_CXX_STANDARD=17 -Dbuild_tests=OFF -DCMAKE_INSTALL_RPATH_USE_LINK_PATH=ON -DCMAKE_INSTALL_PREFIX=$HOME/sw/summit/lapackpp-master + cmake --build src/lapackpp-summit-build --target install --parallel 10 + Optionally, download and install Python packages for :ref:`PICMI ` or dynamic ensemble optimizations (:ref:`libEnsemble `): .. code-block:: bash @@ -319,6 +335,10 @@ For post-processing, most users use Python via OLCFs's `Jupyter service ` for more information). @@ -2148,6 +2159,18 @@ BackTransformed Diagnostics (with support for Plotfile/openPMD output) * ``.num_snapshots_lab`` (`integer`) Only used when ``.diag_type`` is ``BackTransformed``. The number of lab-frame snapshots that will be written. + Only this option or ``intervals`` should be specified; + a run-time error occurs if the user attempts to set both ``num_snapshots_lab`` and ``intervals``. + +* ``.intervals`` (`string`) + Only used when ``.diag_type`` is ``BackTransformed``. + Using the `Intervals parser`_ syntax, this string defines the lab frame times at which data is dumped, + given as multiples of the step size ``dt_snapshots_lab`` or ``dz_snapshots_lab`` described below. + Example: ``btdiag1.intervals = 10:11,20:24:2`` and ``btdiag1.dt_snapshots_lab = 1.e-12`` + indicate to dump at lab times ``1e-11``, ``1.1e-11``, ``2e-11``, ``2.2e-11``, and ``2.4e-11`` seconds. + Note that the stop interval, the second number in the slice, must always be specified. + Only this option or ``num_snapshots_lab`` should be specified; + a run-time error occurs if the user attempts to set both ``num_snapshots_lab`` and ``intervals``. * ``.dt_snapshots_lab`` (`float`, in seconds) Only used when ``.diag_type`` is ``BackTransformed``. @@ -2454,39 +2477,42 @@ Reduced Diagnostics sum of the particles' weight of each species. * ``BeamRelevant`` - This type computes properties of a particle beam relevant for particle accelerators, - like position, momentum, emittance, etc. + This type computes properties of a particle beam relevant for particle accelerators, like position, momentum, emittance, etc. - ``.species`` must be provided, - such that the diagnostics are done for this (beam-like) species only. + ``.species`` must be provided, such that the diagnostics are done for this (beam-like) species only. + + The output columns (for 3D-XYZ) are the following, where the average is done over the whole species (typical usage: the particle beam is in a separate species): + + [0]: simulation step (iteration). - The output columns (for 3D-XYZ) are the following, where the average is done over - the whole species (typical usage: the particle beam is in a separate species): + [1]: time (s). - [1], [2], [3]: The mean values of beam positions (m) - :math:`\langle x \rangle`, :math:`\langle y \rangle`, + [2], [3], [4]: The mean values of beam positions (m) + :math:`\langle x \rangle`, + :math:`\langle y \rangle`, :math:`\langle z \rangle`. - [4], [5], [6]: The mean values of beam relativistic momenta (kg m/s) - :math:`\langle p_x \rangle`, :math:`\langle p_y \rangle`, + [5], [6], [7]: The mean values of beam relativistic momenta (kg m/s) + :math:`\langle p_x \rangle`, + :math:`\langle p_y \rangle`, :math:`\langle p_z \rangle`. - [7]: The mean Lorentz factor :math:`\langle \gamma \rangle`. + [8]: The mean Lorentz factor :math:`\langle \gamma \rangle`. - [8], [9], [10]: The RMS values of beam positions (m) + [9], [10], [11]: The RMS values of beam positions (m) :math:`\delta_x = \sqrt{ \langle (x - \langle x \rangle)^2 \rangle }`, :math:`\delta_y = \sqrt{ \langle (y - \langle y \rangle)^2 \rangle }`, :math:`\delta_z = \sqrt{ \langle (z - \langle z \rangle)^2 \rangle }`. - [11], [12], [13]: The RMS values of beam relativistic momenta (kg m/s) + [12], [13], [14]: The RMS values of beam relativistic momenta (kg m/s) :math:`\delta_{px} = \sqrt{ \langle (p_x - \langle p_x \rangle)^2 \rangle }`, :math:`\delta_{py} = \sqrt{ \langle (p_y - \langle p_y \rangle)^2 \rangle }`, :math:`\delta_{pz} = \sqrt{ \langle (p_z - \langle p_z \rangle)^2 \rangle }`. - [14]: The RMS value of the Lorentz factor + [15]: The RMS value of the Lorentz factor :math:`\sqrt{ \langle (\gamma - \langle \gamma \rangle)^2 \rangle }`. - [15], [16], [17]: beam projected transverse RMS normalized emittance (m) + [16], [17], [18]: beam projected transverse RMS normalized emittance (m) :math:`\epsilon_x = \dfrac{1}{mc} \sqrt{\delta_x^2 \delta_{px}^2 - \Big\langle (x-\langle x \rangle) (p_x-\langle p_x \rangle) \Big\rangle^2}`, :math:`\epsilon_y = \dfrac{1}{mc} \sqrt{\delta_y^2 \delta_{py}^2 - @@ -2494,12 +2520,16 @@ Reduced Diagnostics :math:`\epsilon_z = \dfrac{1}{mc} \sqrt{\delta_z^2 \delta_{pz}^2 - \Big\langle (z-\langle z \rangle) (p_z-\langle p_z \rangle) \Big\rangle^2}`. - [18]: The charge of the beam (C). + [19], [20]: beta function for the transverse directions (m) + :math:`\beta_x = \dfrac{{\delta_x}^2}{\epsilon_x}`, + :math:`\beta_y = \dfrac{{\delta_y}^2}{\epsilon_y}`. + + [21]: The charge of the beam (C). For 2D-XZ, :math:`\langle y \rangle`, :math:`\delta_y`, and - :math:`\epsilon_y` will not be outputed. + :math:`\epsilon_y` will not be outputted. * ``LoadBalanceCosts`` This type computes the cost, used in load balancing, for each box on the domain. diff --git a/Docs/source/usage/python.rst b/Docs/source/usage/python.rst index 16836c75418..afd8c5c0cac 100644 --- a/Docs/source/usage/python.rst +++ b/Docs/source/usage/python.rst @@ -7,38 +7,25 @@ WarpX uses the `PICMI standard `__ for Python version 3.7 or newer is required. Example input files can be found in :ref:`the examples section `. -The examples support running in both modes by commenting and uncommenting the appropriate lines. +In the input file, instances of classes are created defining the various aspects of the simulation. +The `Simulation` object is the central object, where the instances are passed, +defining the simulation time, field solver, registered species, etc. .. _usage-picmi-parameters: -Parameters +Classes ---------- Simulation and grid setup ^^^^^^^^^^^^^^^^^^^^^^^^^ -The `Simulation` object is the central object in a PICMI script. -It defines the simulation time, field solver, registered species, etc. - -.. autoclass:: picmistandard.PICMI_Simulation +Simulation +"""""""""" +.. autoclass:: pywarpx.picmi.Simulation :members: step, add_species, add_laser, write_input_file -Field solvers define the updates of electric and magnetic fields. - -.. autoclass:: picmistandard.PICMI_ElectromagneticSolver - -.. autoclass:: picmistandard.PICMI_ElectrostaticSolver - -Grid define the geometry and discretization. - -.. autoclass:: picmistandard.PICMI_Cartesian3DGrid - -.. autoclass:: picmistandard.PICMI_Cartesian2DGrid - -.. autoclass:: picmistandard.PICMI_Cartesian1DGrid - -.. autoclass:: picmistandard.PICMI_CylindricalGrid - +Constants +""""""""" For convenience, the PICMI interface defines the following constants, which can be used directly inside any PICMI script. The values are in SI units. @@ -49,27 +36,83 @@ which can be used directly inside any PICMI script. The values are in SI units. - ``picmi.constants.m_e``: The electron mass - ``picmi.constants.m_p``: The proton mass -Additionally to self-consistent fields from the field solver, external fields can be applied. +Field solvers define the updates of electric and magnetic fields. + +ElectromagneticSolver +""""""""""""""""""""" +.. autoclass:: pywarpx.picmi.ElectromagneticSolver -.. autoclass:: picmistandard.PICMI_ConstantAppliedField +ElectrostaticSolver +""""""""""""""""""" +.. autoclass:: pywarpx.picmi.ElectrostaticSolver -.. autoclass:: picmistandard.PICMI_AnalyticAppliedField +Cartesian3DGrid +""""""""""""""" +.. autoclass:: pywarpx.picmi.Cartesian3DGrid -.. autoclass:: picmistandard.PICMI_Mirror +Cartesian2DGrid +""""""""""""""" +.. autoclass:: pywarpx.picmi.Cartesian2DGrid -Diagnostics can be used to output data. +Cartesian1DGrid +""""""""""""""" +.. autoclass:: pywarpx.picmi.Cartesian1DGrid -.. autoclass:: picmistandard.PICMI_ParticleDiagnostic +CylindricalGrid +""""""""""""""" +.. autoclass:: pywarpx.picmi.CylindricalGrid -.. autoclass:: picmistandard.PICMI_FieldDiagnostic +EmbeddedBoundary +"""""""""""""""" +.. autoclass:: pywarpx.picmi.EmbeddedBoundary -.. autoclass:: picmistandard.PICMI_ElectrostaticFieldDiagnostic +Applied fields +^^^^^^^^^^^^^^ + +ConstantAppliedField +"""""""""""""""""""" +.. autoclass:: pywarpx.picmi.ConstantAppliedField + +AnalyticAppliedField +"""""""""""""""""""" +.. autoclass:: pywarpx.picmi.AnalyticAppliedField + +PlasmaLens +"""""""""" +.. autoclass:: pywarpx.picmi.PlasmaLens + +Mirror +"""""" +.. autoclass:: pywarpx.picmi.Mirror + +Diagnostics +^^^^^^^^^^^ + +ParticleDiagnostic +"""""""""""""""""" +.. autoclass:: pywarpx.picmi.ParticleDiagnostic + +FieldDiagnostic +""""""""""""""" +.. autoclass:: pywarpx.picmi.FieldDiagnostic + +ElectrostaticFieldDiagnostic +"""""""""""""""""""""""""""" +.. autoclass:: pywarpx.picmi.ElectrostaticFieldDiagnostic Lab-frame diagnostics diagnostics are used when running boosted-frame simulations. -.. autoclass:: picmistandard.PICMI_LabFrameParticleDiagnostic +LabFrameParticleDiagnostic +"""""""""""""""""""""""""" +.. autoclass:: pywarpx.picmi.LabFrameParticleDiagnostic + +LabFrameFieldDiagnostic +""""""""""""""""""""""" +.. autoclass:: pywarpx.picmi.LabFrameFieldDiagnostic -.. autoclass:: picmistandard.PICMI_LabFrameFieldDiagnostic +Checkpoint +"""""""""" +.. autoclass:: pywarpx.picmi.Checkpoint Particles ^^^^^^^^^ @@ -77,38 +120,70 @@ Particles Species objects are a collection of particles with similar properties. For instance, background plasma electrons, background plasma ions and an externally injected beam could each be their own particle species. -.. autoclass:: picmistandard.PICMI_Species +Species +""""""" +.. autoclass:: pywarpx.picmi.Species -.. autoclass:: picmistandard.PICMI_MultiSpecies +MultiSpecies +"""""""""""" +.. autoclass:: pywarpx.picmi.MultiSpecies Particle distributions can be used for to initialize particles in a particle species. -.. autoclass:: picmistandard.PICMI_GaussianBunchDistribution +GaussianBunchDistribution +""""""""""""""""""""""""" +.. autoclass:: pywarpx.picmi.GaussianBunchDistribution -.. autoclass:: picmistandard.PICMI_UniformDistribution +UniformDistribution +""""""""""""""""""" +.. autoclass:: pywarpx.picmi.UniformDistribution -.. autoclass:: picmistandard.PICMI_AnalyticDistribution +AnalyticDistribution +"""""""""""""""""""" +.. autoclass:: pywarpx.picmi.AnalyticDistribution -.. autoclass:: picmistandard.PICMI_ParticleListDistribution +ParticleListDistribution +"""""""""""""""""""""""" +.. autoclass:: pywarpx.picmi.ParticleListDistribution Particle layouts determine how to microscopically place macro particles in a grid cell. -.. autoclass:: picmistandard.PICMI_GriddedLayout +GriddedLayout +""""""""""""" +.. autoclass:: pywarpx.picmi.GriddedLayout + +PseudoRandomLayout +"""""""""""""""""" +.. autoclass:: pywarpx.picmi.PseudoRandomLayout -.. autoclass:: picmistandard.PICMI_PseudoRandomLayout +Other operations related to particles + +CoulombCollisions +""""""""""""""""" +.. autoclass:: pywarpx.picmi.CoulombCollisions + +MCCCollisions +""""""""""""" +.. autoclass:: pywarpx.picmi.MCCCollisions Lasers ^^^^^^ Laser profiles can be used to initialize laser pulses in the simulation. -.. autoclass:: picmistandard.PICMI_GaussianLaser +GaussianLaser +""""""""""""" +.. autoclass:: pywarpx.picmi.GaussianLaser -.. autoclass:: picmistandard.PICMI_AnalyticLaser +AnalyticLaser +""""""""""""" +.. autoclass:: pywarpx.picmi.AnalyticLaser Laser injectors control where to initialize laser pulses on the simulation grid. -.. autoclass:: picmistandard.PICMI_LaserAntenna +LaserAntenna +"""""""""""" +.. autoclass:: pywarpx.picmi.LaserAntenna .. _usage-picmi-run: @@ -119,6 +194,7 @@ Running WarpX can be run in one of two modes. It can run as a preprocessor, using the Python input file to generate an input file to be used by the C++ version, or it can be run directly from Python. +The examples support running in both modes by commenting and uncommenting the appropriate lines. In either mode, if using a `virtual environment `__, be sure to activate it before compiling and running WarpX. diff --git a/Docs/spack.yaml b/Docs/spack.yaml index 434d0820929..1e953943610 100644 --- a/Docs/spack.yaml +++ b/Docs/spack.yaml @@ -23,4 +23,5 @@ spack: - py-breathe - py-recommonmark - py-pygments + - py-sphinx-design - py-sphinx-rtd-theme diff --git a/Examples/Modules/RigidInjection/inputs_2d_BoostedFrame b/Examples/Modules/RigidInjection/inputs_2d_BoostedFrame index 73301448ef8..bd8ce220acb 100644 --- a/Examples/Modules/RigidInjection/inputs_2d_BoostedFrame +++ b/Examples/Modules/RigidInjection/inputs_2d_BoostedFrame @@ -53,7 +53,7 @@ diagnostics.diags_names = diag1 diag2 diag1.diag_type = BackTransformed diag1.do_back_transformed_fields = 1 -diag1.num_snapshots_lab = 2 +diag1.intervals = :1 diag1.dt_snapshots_lab = 1.8679589331096515e-13 diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho diag1.format = plotfile diff --git a/Examples/Modules/boosted_diags/analysis_3Dbacktransformed_diag.py b/Examples/Modules/boosted_diags/analysis_3Dbacktransformed_diag.py index 2dac6d245b2..6fa4e9c936f 100755 --- a/Examples/Modules/boosted_diags/analysis_3Dbacktransformed_diag.py +++ b/Examples/Modules/boosted_diags/analysis_3Dbacktransformed_diag.py @@ -16,30 +16,65 @@ between the full back-transformed diagnostic and the reduced diagnostic (i.e., x-z slice) . ''' +import os +import sys + import numpy as np +import openpmd_api as io import read_raw_data +import yt + +yt.funcs.mylog.setLevel(0) + +sys.path.insert(1, '../../../../warpx/Regression/Checksum/') +import checksumAPI + +filename = sys.argv[1] -# Read data from back-transformed diagnostics of entire domain -snapshot = './lab_frame_data/snapshots/snapshot00002' +# Tolerances to check consistency between legacy BTD and new BTD +rtol = 1e-16 +atol = 1e-16 + +# Read data from legacy back-transformed diagnostics (entire domain) +snapshot = './lab_frame_data/snapshots/snapshot00003' header = './lab_frame_data/snapshots/Header' allrd, info = read_raw_data.read_lab_snapshot(snapshot, header) -F = allrd['Ez'] -print("F.shape ", F.shape) -F_1D = np.squeeze(F[F.shape[0]//2,F.shape[1]//2,:]) - +Ez_legacy = allrd['Ez'] +print(f'Ez_legacy.shape = {Ez_legacy.shape}') +Ez_legacy_1D = np.squeeze(Ez_legacy[Ez_legacy.shape[0]//2,Ez_legacy.shape[1]//2,:]) -# Read data from reduced back-transformed diagnostics (i.e. slice) -snapshot_slice = './lab_frame_data/slices/slice00002' +# Read data from reduced back-transformed diagnostics (slice) +snapshot_slice = './lab_frame_data/slices/slice00003' header_slice = './lab_frame_data/slices/Header' allrd, info = read_raw_data.read_lab_snapshot(snapshot_slice, header_slice) -Fs = allrd['Ez'] -print("Fs.shape", Fs.shape) -Fs_1D = np.squeeze(Fs[Fs.shape[0]//2,1,:]) +Ez_legacy_slice = allrd['Ez'] +print(f'Ez_legacy_slice.shape = {Ez_legacy_slice.shape}') +Ez_legacy_slice_1D = np.squeeze(Ez_legacy_slice[Ez_legacy_slice.shape[0]//2,1,:]) + +# Read data from new back-transformed diagnostics (plotfile) +ds_plotfile = yt.load(filename) +data = ds_plotfile.covering_grid( + level=0, + left_edge=ds_plotfile.domain_left_edge, + dims=ds_plotfile.domain_dimensions) +Ez_plotfile = data[('mesh', 'Ez')].to_ndarray() + +# Read data from new back-transformed diagnostics (openPMD) +series = io.Series("./diags/diag2/openpmd_%T.h5", io.Access.read_only) +ds_openpmd = series.iterations[3] +Ez_openpmd = ds_openpmd.meshes['E']['z'].load_chunk() +Ez_openpmd = Ez_openpmd.transpose() +series.flush() -error_rel = np.max(np.abs(Fs_1D - F_1D)) / np.max(np.abs(F_1D)) -tolerance_rel = 1E-15 +# Compare arrays to check consistency between new BTD formats (plotfile and openPMD) +assert(np.allclose(Ez_plotfile, Ez_openpmd, rtol=rtol, atol=atol)) -print("error_rel : " + str(error_rel)) -print("tolerance_rel: " + str(tolerance_rel)) +# Check slicing +err = np.max(np.abs(Ez_legacy_slice_1D-Ez_legacy_1D)) / np.max(np.abs(Ez_legacy_1D)) +tol = 1e-16 +print(f'error = {err}') +print(f'tolerance = {tol}') +assert(err < tol) -assert( error_rel < tolerance_rel ) +test_name = os.path.split(os.getcwd())[1] +checksumAPI.evaluate_checksum(test_name, filename) diff --git a/Examples/Modules/boosted_diags/inputs_3d_slice b/Examples/Modules/boosted_diags/inputs_3d_slice index 5c774fc893e..c7938991935 100644 --- a/Examples/Modules/boosted_diags/inputs_3d_slice +++ b/Examples/Modules/boosted_diags/inputs_3d_slice @@ -113,27 +113,25 @@ slice.dt_slice_snapshots_lab = 3.3356409519815207e-12 slice.particle_slice_width_lab = 2.e-6 # Diagnostics -diagnostics.diags_names = diag1 btd_openpmd btd_pltfile -diag1.intervals = 10000 -diag1.diag_type = Full - -btd_openpmd.diag_type = BackTransformed -btd_openpmd.do_back_transformed_fields = 1 -btd_openpmd.num_snapshots_lab = 4 -btd_openpmd.dz_snapshots_lab = 0.001 -btd_openpmd.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho -btd_openpmd.format = openpmd -btd_openpmd.buffer_size = 32 -btd_openpmd.openpmd_backend = h5 - -btd_pltfile.diag_type = BackTransformed -btd_pltfile.do_back_transformed_fields = 1 -btd_pltfile.num_snapshots_lab = 4 -btd_pltfile.dz_snapshots_lab = 0.001 -btd_pltfile.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho -btd_pltfile.format = plotfile -btd_pltfile.buffer_size = 32 -btd_pltfile.write_species = 1 +diagnostics.diags_names = diag1 diag2 + +diag1.diag_type = BackTransformed +diag1.do_back_transformed_fields = 1 +diag1.num_snapshots_lab = 4 +diag1.dz_snapshots_lab = 0.001 +diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho +diag1.format = plotfile +diag1.buffer_size = 32 +diag1.write_species = 1 + +diag2.diag_type = BackTransformed +diag2.do_back_transformed_fields = 1 +diag2.intervals = 0:4:2, 1:3:2 +diag2.dz_snapshots_lab = 0.001 +diag2.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho +diag2.format = openpmd +diag2.buffer_size = 32 +diag2.openpmd_backend = h5 # old BTD diagnostics warpx.do_back_transformed_diagnostics = 1 diff --git a/Examples/Modules/nuclear_fusion/analysis_deuterium_tritium_fusion.py b/Examples/Modules/nuclear_fusion/analysis_deuterium_tritium_fusion.py index f218df54d70..b5f5da683b2 100755 --- a/Examples/Modules/nuclear_fusion/analysis_deuterium_tritium_fusion.py +++ b/Examples/Modules/nuclear_fusion/analysis_deuterium_tritium_fusion.py @@ -6,6 +6,7 @@ # License: BSD-3-Clause-LBNL import os +import re import sys import yt @@ -66,15 +67,30 @@ E_fusion = 17.5893*MeV_to_Joule # Energy released during the fusion reaction +## Checks whether this is the 2D or the 3D test +warpx_used_inputs = open('./warpx_used_inputs', 'r').read() +if re.search('geometry.dims = RZ', warpx_used_inputs): + is_RZ = True +else: + is_RZ = False + ## Some numerical parameters for this test size_x = 8 size_y = 8 size_z = 16 -dV_total = size_x*size_y*size_z # Total simulation volume +if is_RZ: + dV_slice = np.pi * size_x**2 + yt_z_string = "particle_position_y" + nppcell_1 = 10000*8 + nppcell_2 = 900*8 +else: + dV_slice = size_x*size_y + yt_z_string = "particle_position_z" + nppcell_1 = 10000 + nppcell_2 = 900 # Volume of a slice corresponding to a single cell in the z direction. In tests 1 and 2, all the # particles of a given species in the same slice have the exact same momentum -dV_slice = size_x*size_y -dt = 1./(scc.c*np.sqrt(3.)) + # In test 1 and 2, the energy in cells number i (in z direction) is typically Energy_step * i**2 Energy_step = 22.*keV_to_Joule @@ -89,7 +105,7 @@ def add_existing_species_to_dict(yt_ad, data_dict, species_name, prefix, suffix) data_dict[prefix+"_w_"+suffix] = yt_ad[species_name, "particle_weight"].v data_dict[prefix+"_id_"+suffix] = yt_ad[species_name, "particle_id"].v data_dict[prefix+"_cpu_"+suffix] = yt_ad[species_name, "particle_cpu"].v - data_dict[prefix+"_z_"+suffix] = yt_ad[species_name, "particle_position_z"].v + data_dict[prefix+"_z_"+suffix] = yt_ad[species_name, yt_z_string].v def add_empty_species_to_dict(data_dict, species_name, prefix, suffix): data_dict[prefix+"_px_"+suffix] = np.empty(0) @@ -263,8 +279,11 @@ def expected_weight_com(E_com, reactant0_density, reactant1_density, dV, dt): def check_macroparticle_number(data, fusion_probability_target_value, num_pair_per_cell): ## Checks that the number of macroparticles is as expected for the first and second tests - ## The first slice 0 < z < 1 does not contribute to product species creation - numcells = dV_total - dV_slice + ## The first slice 0 < z < 1 does not contribute to alpha creation + if is_RZ: + numcells = size_x*(size_z-1) + else: + numcells = size_x*size_y*(size_z-1) ## In these tests, the fusion_multiplier is so high that the fusion probability per pair is ## equal to the parameter fusion_probability_target_value fusion_probability_per_pair = fusion_probability_target_value @@ -315,7 +334,7 @@ def compute_E_com2(data): p_reactant0_sq = 2.*mass[reactant_species[0]]*(Energy_step*np.arange(size_z)**2) return p_sq_reactant1_frame_to_E_COM_frame(p_reactant0_sq) -def check_fusion_yield(data, expected_fusion_number, E_com, reactant0_density, reactant1_density): +def check_fusion_yield(data, expected_fusion_number, E_com, reactant0_density, reactant1_density, dt): ## Checks that the fusion yield is as expected for the first and second tests. product_weight_theory = expected_weight_com(E_com/keV_to_Joule, reactant0_density, reactant1_density, dV_slice, dt) @@ -330,24 +349,25 @@ def check_fusion_yield(data, expected_fusion_number, E_com, reactant0_density, r assert(np.all(is_close(product_weight_theory, product_weight_simulation, rtol = 5.*relative_std_weight))) -def specific_check1(data): - check_isotropy(data, relative_tolerance = 3.e-2) +def specific_check1(data, dt): + if not is_RZ: + check_isotropy(data, relative_tolerance = 3.e-2) expected_fusion_number = check_macroparticle_number(data, fusion_probability_target_value = 0.002, - num_pair_per_cell = 10000) + num_pair_per_cell = nppcell_1) E_com = compute_E_com1(data) check_fusion_yield(data, expected_fusion_number, E_com, reactant0_density = 1., - reactant1_density = 1.) + reactant1_density = 1., dt=dt) -def specific_check2(data): +def specific_check2(data, dt): check_xy_isotropy(data) ## Only 900 particles pairs per cell here because we ignore the 10% of reactants that are at rest expected_fusion_number = check_macroparticle_number(data, fusion_probability_target_value = 0.02, - num_pair_per_cell = 900) + num_pair_per_cell = nppcell_2) E_com = compute_E_com2(data) check_fusion_yield(data, expected_fusion_number, E_com, reactant0_density = 1.e20, - reactant1_density = 1.e26) + reactant1_density = 1.e26, dt=dt) def check_charge_conservation(rho_start, rho_end): assert(np.all(is_close(rho_start, rho_end, rtol=2.e-11))) @@ -359,6 +379,7 @@ def main(): ds_start = yt.load(filename_start) ad_end = ds_end.all_data() ad_start = ds_start.all_data() + dt = float(ds_end.current_time - ds_start.current_time) field_data_end = ds_end.covering_grid(level=0, left_edge=ds_end.domain_left_edge, dims=ds_end.domain_dimensions) field_data_start = ds_start.covering_grid(level=0, left_edge=ds_start.domain_left_edge, @@ -379,7 +400,7 @@ def main(): generic_check(data) # Checks that are specific to test number i - eval("specific_check"+str(i)+"(data)") + eval("specific_check"+str(i)+"(data, dt)") rho_start = field_data_start["rho"].to_ndarray() rho_end = field_data_end["rho"].to_ndarray() diff --git a/Examples/Modules/nuclear_fusion/inputs_deuterium_tritium_rz b/Examples/Modules/nuclear_fusion/inputs_deuterium_tritium_rz new file mode 100644 index 00000000000..fb581c82535 --- /dev/null +++ b/Examples/Modules/nuclear_fusion/inputs_deuterium_tritium_rz @@ -0,0 +1,129 @@ +################################# +####### GENERAL PARAMETERS ###### +################################# +## With these parameters, each cell has a size of exactly 1 by 1 by 1 +max_step = 1 +amr.n_cell = 8 16 +amr.max_grid_size = 8 +amr.blocking_factor = 8 +amr.max_level = 0 +geometry.dims = RZ +geometry.prob_lo = 0. 0. +geometry.prob_hi = 8. 16. + +################################# +###### Boundary Condition ####### +################################# +boundary.field_lo = none periodic +boundary.field_hi = pec periodic + +################################# +############ NUMERICS ########### +################################# +warpx.verbose = 1 +warpx.cfl = 1.0 + +# Order of particle shape factors +algo.particle_shape = 1 + +################################# +############ PLASMA ############# +################################# +particles.species_names = deuterium1 tritium1 helium1 neutron1 deuterium2 tritium2 helium2 neutron2 + +my_constants.m_deuterium = 2.01410177812*m_u +my_constants.m_tritium = 3.0160492779*m_u +my_constants.m_reduced = m_deuterium*m_tritium/(m_deuterium+m_tritium) +my_constants.keV_to_J = 1.e3*q_e +my_constants.Energy_step = 22. * keV_to_J + +deuterium1.species_type = deuterium +deuterium1.injection_style = "NRandomPerCell" +deuterium1.num_particles_per_cell = 80000 +deuterium1.profile = constant +deuterium1.density = 1. +deuterium1.momentum_distribution_type = parse_momentum_function +## Thanks to the floor, all particles in the same cell have the exact same momentum +deuterium1.momentum_function_ux(x,y,z) = "u = sqrt(2*m_reduced*Energy_step*(floor(z)**2))/(m_deuterium*clight); if(x*x+y*y>0.0, -u*y/sqrt(x*x+y*y), 0.0)" # azimuthal velocity +deuterium1.momentum_function_uy(x,y,z) = "u = sqrt(2*m_reduced*Energy_step*(floor(z)**2))/(m_deuterium*clight); if(x*x+y*y>0.0, u*x/sqrt(x*x+y*y), 0.0)" # azimuthal velocity +deuterium1.momentum_function_uz(x,y,z) = "0" +deuterium1.do_not_push = 1 +deuterium1.do_not_deposit = 1 + +tritium1.species_type = tritium +tritium1.injection_style = "NRandomPerCell" +tritium1.num_particles_per_cell = 80000 +tritium1.profile = constant +tritium1.density = 1. +tritium1.momentum_distribution_type = "parse_momentum_function" +## Thanks to the floor, all particles in the same cell have the exact same momentum +tritium1.momentum_function_ux(x,y,z) = "u = sqrt(2*m_reduced*Energy_step*(floor(z)**2))/(m_tritium*clight); if(x*x+y*y>0.0, u*y/sqrt(x*x+y*y), 0.0)" # counter-streaming azimuthal velocity +tritium1.momentum_function_uy(x,y,z) = "u = sqrt(2*m_reduced*Energy_step*(floor(z)**2))/(m_tritium*clight); if(x*x+y*y>0.0, -u*x/sqrt(x*x+y*y), 0.0)" # counter-streaming azimuthal velocity +tritium1.momentum_function_uz(x,y,z) = 0 +tritium1.do_not_push = 1 +tritium1.do_not_deposit = 1 + +helium1.species_type = helium4 +helium1.do_not_push = 1 +helium1.do_not_deposit = 1 + +neutron1.species_type = neutron +neutron1.do_not_push = 1 +neutron1.do_not_deposit = 1 + +my_constants.background_dens = 1.e26 +my_constants.beam_dens = 1.e20 + +deuterium2.species_type = deuterium +deuterium2.injection_style = "NRandomPerCell" +deuterium2.num_particles_per_cell = 8000 +deuterium2.profile = "parse_density_function" +## A tenth of the macroparticles in each cell is made of immobile high-density background deuteriums. +## The other nine tenths are made of fast low-density beam deuteriums. +deuterium2.density_function(x,y,z) = if(y - floor(y) < 0.1, 10.*background_dens, 10./9.*beam_dens) +deuterium2.momentum_distribution_type = "parse_momentum_function" +deuterium2.momentum_function_ux(x,y,z) = 0. +deuterium2.momentum_function_uy(x,y,z) = 0. +deuterium2.momentum_function_uz(x,y,z) = "if(y - floor(y) < 0.1, + 0., sqrt(2*m_deuterium*Energy_step*(floor(z)**2))/(m_deuterium*clight))" +deuterium2.do_not_push = 1 +deuterium2.do_not_deposit = 1 + +tritium2.species_type = tritium +tritium2.injection_style = "NRandomPerCell" +tritium2.num_particles_per_cell = 800 +tritium2.profile = constant +tritium2.density = background_dens +tritium2.momentum_distribution_type = "constant" +tritium2.do_not_push = 1 +tritium2.do_not_deposit = 1 + +helium2.species_type = helium4 +helium2.do_not_push = 1 +helium2.do_not_deposit = 1 + +neutron2.species_type = neutron +neutron2.do_not_push = 1 +neutron2.do_not_deposit = 1 + +################################# +############ COLLISION ########## +################################# +collisions.collision_names = DTF1 DTF2 + +DTF1.species = deuterium1 tritium1 +DTF1.product_species = helium1 neutron1 +DTF1.type = nuclearfusion +DTF1.fusion_multiplier = 1.e50 + +DTF2.species = deuterium2 tritium2 +DTF2.product_species = helium2 neutron2 +DTF2.type = nuclearfusion +DTF2.fusion_multiplier = 1.e15 +DTF2.fusion_probability_target_value = 0.02 + +# Diagnostics +diagnostics.diags_names = diag1 +diag1.intervals = 1 +diag1.diag_type = Full +diag1.fields_to_plot = rho diff --git a/Examples/Physics_applications/capacitive_discharge/PICMI_inputs_1d.py b/Examples/Physics_applications/capacitive_discharge/PICMI_inputs_1d.py index 5eac1e172ff..21a631505dd 100644 --- a/Examples/Physics_applications/capacitive_discharge/PICMI_inputs_1d.py +++ b/Examples/Physics_applications/capacitive_discharge/PICMI_inputs_1d.py @@ -17,8 +17,9 @@ class PoissonSolver1D(picmi.ElectrostaticSolver): - """This can be removed and the MLMG solver used once - https://github.com/ECP-WarpX/WarpX/issues/3123 is addressed.""" + """This solver is maintained as an example of the use of Python callbacks. + However, it is not necessarily needed since the 1D code has the direct tridiagonal + solver implemented.""" def __init__(self, grid, **kwargs): """Direct solver for the Poisson equation using superLU. This solver is @@ -60,8 +61,8 @@ def initialize_inputs(self): super(PoissonSolver1D, self).initialize_inputs() - self.nz = self.grid.nx - self.dz = (self.grid.xmax - self.grid.xmin) / self.nz + self.nz = self.grid.number_of_cells[0] + self.dz = (self.grid.upper_bound[0] - self.grid.lower_bound[0]) / self.nz self.nxguardphi = 1 self.nzguardphi = 1 @@ -77,18 +78,6 @@ def decompose_matrix(self): system.""" self.nsolve = self.nz + 1 - # Set up the tridiagonal computation matrix in order to solve A*phi = - # rho for phi. - self.A_ldiag = np.ones(self.nsolve-1) / self.dz**2 - self.A_mdiag = -2.*np.ones(self.nsolve) / self.dz**2 - self.A_udiag = np.ones(self.nsolve-1) / self.dz**2 - - self.A_mdiag[0] = 1. - self.A_udiag[0] = 0.0 - - self.A_mdiag[-1] = 1. - self.A_ldiag[-1] = 0.0 - # Set up the computation matrix in order to solve A*phi = rho A = np.zeros((self.nsolve, self.nsolve)) idx = np.arange(self.nsolve) @@ -101,7 +90,7 @@ def decompose_matrix(self): A[0, 0] = 1.0 A[-1, -1] = 1.0 - A = csc_matrix(A, dtype=np.float32) + A = csc_matrix(A, dtype=np.float64) self.lu = sla.splu(A) def _run_solve(self): @@ -124,7 +113,7 @@ def solve(self): # Construct b vector rho = -self.rho_data / constants.ep0 - b = np.zeros(rho.shape[0], dtype=np.float32) + b = np.zeros(rho.shape[0], dtype=np.float64) b[:] = rho * self.dz**2 b[0] = left_voltage @@ -169,10 +158,11 @@ class CapacitiveDischargeExample(object): # Time (in seconds) between diagnostic evaluations diag_interval = 32 / freq - def __init__(self, n=0, test=False): + def __init__(self, n=0, test=False, pythonsolver=False): """Get input parameters for the specific case (n) desired.""" self.n = n self.test = test + self.pythonsolver = pythonsolver # Case specific input parameters self.voltage = f"{self.voltage[n]}*sin(2*pi*{self.freq:.5e}*t)" @@ -221,11 +211,11 @@ def setup_run(self): # Field solver # ####################################################################### - # self.solver = picmi.ElectrostaticSolver( - # grid=self.grid, method='Multigrid', required_precision=1e-6, - # warpx_self_fields_verbosity=2 - # ) - self.solver = PoissonSolver1D(grid=self.grid) + if self.pythonsolver: + self.solver = PoissonSolver1D(grid=self.grid) + else: + # This will use the tridiagonal solver + self.solver = picmi.ElectrostaticSolver(grid=self.grid) ####################################################################### # Particle types setup # @@ -329,13 +319,18 @@ def setup_run(self): # Add diagnostics for the CI test to be happy # ####################################################################### + if self.pythonsolver: + file_prefix = 'Python_background_mcc_1d_plt' + else: + file_prefix = 'Python_background_mcc_1d_tridiag_plt' + field_diag = picmi.FieldDiagnostic( name='diag1', grid=self.grid, period=0, data_list=['rho_electrons', 'rho_he_ions'], write_dir='.', - warpx_file_prefix='Python_background_mcc_1d_plt' + warpx_file_prefix=file_prefix ) self.sim.add_diagnostic(field_diag) @@ -358,6 +353,19 @@ def run_sim(self): if self.sim.extension.getMyProc() == 0: np.save(f'ion_density_case_{self.n+1}.npy', self.ion_density_array) + # query the particle z-coordinates if this is run during CI testing + # to cover that functionality + if self.test: + nparts = self.sim.extension.get_particle_count( + 'he_ions', local=True + ) + z_coords = np.concatenate( + self.sim.extension.get_particle_z('he_ions') + ) + assert len(z_coords) == nparts + assert np.all(z_coords >= 0.0) and np.all(z_coords <= self.gap) + + ########################## # parse input parameters ########################## @@ -371,11 +379,15 @@ def run_sim(self): '-n', help='Test number to run (1 to 4)', required=False, type=int, default=1 ) +parser.add_argument( + '--pythonsolver', help='toggle whether to use the Python level solver', + action='store_true' +) args, left = parser.parse_known_args() sys.argv = sys.argv[:1]+left if args.n < 1 or args.n > 4: raise AttributeError('Test number must be an integer from 1 to 4.') -run = CapacitiveDischargeExample(n=args.n-1, test=args.test) +run = CapacitiveDischargeExample(n=args.n-1, test=args.test, pythonsolver=args.pythonsolver) run.run_sim() diff --git a/Examples/Physics_applications/capacitive_discharge/PICMI_inputs_2d.py b/Examples/Physics_applications/capacitive_discharge/PICMI_inputs_2d.py index d5a18071ad4..85cc04148d3 100755 --- a/Examples/Physics_applications/capacitive_discharge/PICMI_inputs_2d.py +++ b/Examples/Physics_applications/capacitive_discharge/PICMI_inputs_2d.py @@ -91,10 +91,10 @@ def initialize_inputs(self): super(PoissonSolverPseudo1D, self).initialize_inputs() - self.nx = self.grid.nx - self.nz = self.grid.ny - self.dx = (self.grid.xmax - self.grid.xmin) / self.nx - self.dz = (self.grid.ymax - self.grid.ymin) / self.nz + self.nx = self.grid.number_of_cells[0] + self.nz = self.grid.number_of_cells[1] + self.dx = (self.grid.upper_bound[0] - self.grid.lower_bound[0]) / self.nx + self.dz = (self.grid.upper_bound[1] - self.grid.lower_bound[1]) / self.nz if not np.isclose(self.dx, self.dz): raise RuntimeError('Direct solver requires dx = dz.') diff --git a/Examples/Physics_applications/capacitive_discharge/analysis_1d.py b/Examples/Physics_applications/capacitive_discharge/analysis_1d.py index f9e8be0844e..73b1ce6d600 100755 --- a/Examples/Physics_applications/capacitive_discharge/analysis_1d.py +++ b/Examples/Physics_applications/capacitive_discharge/analysis_1d.py @@ -1,43 +1,43 @@ #!/usr/bin/env python3 -# Copyright 2021 Modern Electron +# Copyright 2022 Modern Electron, David Grote import numpy as np ref_density = np.array([ - 1.29556695e+14, 2.24358819e+14, 2.55381744e+14, 2.55655005e+14, - 2.55796267e+14, 2.55819109e+14, 2.55819687e+14, 2.55751184e+14, - 2.55920806e+14, 2.56072344e+14, 2.55937266e+14, 2.55849080e+14, - 2.55918981e+14, 2.55980835e+14, 2.56054153e+14, 2.56074693e+14, - 2.56036953e+14, 2.56181152e+14, 2.56322617e+14, 2.56253541e+14, - 2.56196222e+14, 2.56353087e+14, 2.56256025e+14, 2.55928999e+14, - 2.56110985e+14, 2.56658916e+14, 2.56832588e+14, 2.56551873e+14, - 2.56491185e+14, 2.56469926e+14, 2.56418623e+14, 2.56541073e+14, - 2.56513772e+14, 2.56424505e+14, 2.56302757e+14, 2.56242393e+14, - 2.56270399e+14, 2.56178952e+14, 2.56071408e+14, 2.56141949e+14, - 2.56419807e+14, 2.56606936e+14, 2.56437775e+14, 2.56252446e+14, - 2.56309518e+14, 2.56383487e+14, 2.56265139e+14, 2.56167672e+14, - 2.56466917e+14, 2.56924869e+14, 2.56901785e+14, 2.56631493e+14, - 2.56643456e+14, 2.56523464e+14, 2.56378270e+14, 2.56571296e+14, - 2.56794304e+14, 2.56788544e+14, 2.56549715e+14, 2.56303160e+14, - 2.56210813e+14, 2.56418356e+14, 2.57314522e+14, 2.58471390e+14, - 2.58169771e+14, 2.56946438e+14, 2.56726546e+14, 2.56853122e+14, - 2.56613699e+14, 2.56509534e+14, 2.56692972e+14, 2.56705133e+14, - 2.56372142e+14, 2.56167556e+14, 2.56296946e+14, 2.56498752e+14, - 2.56523102e+14, 2.56404334e+14, 2.56227096e+14, 2.56398997e+14, - 2.56614907e+14, 2.56436657e+14, 2.56388606e+14, 2.56553679e+14, - 2.56637914e+14, 2.56407785e+14, 2.56104131e+14, 2.56082340e+14, - 2.56095275e+14, 2.56278448e+14, 2.56808134e+14, 2.57127897e+14, - 2.56858174e+14, 2.56326990e+14, 2.56296032e+14, 2.56563349e+14, - 2.56482273e+14, 2.56667481e+14, 2.57072448e+14, 2.56767530e+14, - 2.56433245e+14, 2.56586570e+14, 2.56636412e+14, 2.56765628e+14, - 2.56868130e+14, 2.56783441e+14, 2.56714518e+14, 2.56651014e+14, - 2.56528394e+14, 2.56227520e+14, 2.56163301e+14, 2.56408207e+14, - 2.56433120e+14, 2.56374745e+14, 2.56542028e+14, 2.56748796e+14, - 2.56715201e+14, 2.56298164e+14, 2.56042658e+14, 2.56292455e+14, - 2.56352282e+14, 2.56370562e+14, 2.56487458e+14, 2.56483667e+14, - 2.56741201e+14, 2.56665100e+14, 2.56523784e+14, 2.24741564e+14, - 1.28486944e+14 + 1.29556694e+14, 2.24358818e+14, 2.55381745e+14, 2.55655005e+14, + 2.55796268e+14, 2.55819108e+14, 2.55819686e+14, 2.55751184e+14, + 2.55920806e+14, 2.56072344e+14, 2.55937266e+14, 2.55849080e+14, + 2.55918981e+14, 2.55980835e+14, 2.56054153e+14, 2.56074694e+14, + 2.56036953e+14, 2.56181153e+14, 2.56322618e+14, 2.56253541e+14, + 2.56196224e+14, 2.56353090e+14, 2.56256022e+14, 2.55928997e+14, + 2.56110988e+14, 2.56658917e+14, 2.56832584e+14, 2.56551871e+14, + 2.56491186e+14, 2.56469928e+14, 2.56418625e+14, 2.56541071e+14, + 2.56513773e+14, 2.56424507e+14, 2.56302757e+14, 2.56242392e+14, + 2.56270399e+14, 2.56178952e+14, 2.56071407e+14, 2.56141949e+14, + 2.56419808e+14, 2.56606936e+14, 2.56437774e+14, 2.56252443e+14, + 2.56309513e+14, 2.56383484e+14, 2.56265140e+14, 2.56167674e+14, + 2.56466922e+14, 2.56924871e+14, 2.56901781e+14, 2.56631494e+14, + 2.56643458e+14, 2.56523464e+14, 2.56378273e+14, 2.56571301e+14, + 2.56794308e+14, 2.56788543e+14, 2.56549712e+14, 2.56303156e+14, + 2.56210811e+14, 2.56418363e+14, 2.57314552e+14, 2.58471405e+14, + 2.58169740e+14, 2.56946418e+14, 2.56726550e+14, 2.56853119e+14, + 2.56613698e+14, 2.56509538e+14, 2.56692976e+14, 2.56705132e+14, + 2.56372135e+14, 2.56167561e+14, 2.56296953e+14, 2.56498746e+14, + 2.56523099e+14, 2.56404333e+14, 2.56227098e+14, 2.56399004e+14, + 2.56614905e+14, 2.56436650e+14, 2.56388608e+14, 2.56553683e+14, + 2.56637912e+14, 2.56407782e+14, 2.56104130e+14, 2.56082338e+14, + 2.56095272e+14, 2.56278448e+14, 2.56808134e+14, 2.57127896e+14, + 2.56858173e+14, 2.56326991e+14, 2.56296032e+14, 2.56563348e+14, + 2.56482274e+14, 2.56667483e+14, 2.57072448e+14, 2.56767529e+14, + 2.56433245e+14, 2.56586564e+14, 2.56636403e+14, 2.56765624e+14, + 2.56868122e+14, 2.56783435e+14, 2.56714527e+14, 2.56651030e+14, + 2.56528399e+14, 2.56227514e+14, 2.56163300e+14, 2.56408217e+14, + 2.56433124e+14, 2.56374737e+14, 2.56542023e+14, 2.56748800e+14, + 2.56715205e+14, 2.56298166e+14, 2.56042658e+14, 2.56292458e+14, + 2.56352283e+14, 2.56370559e+14, 2.56487462e+14, 2.56483655e+14, + 2.56741185e+14, 2.56665111e+14, 2.56523794e+14, 2.24741566e+14, + 1.28486948e+14 ]) density_data = np.load( 'ion_density_case_1.npy' ) diff --git a/Examples/Physics_applications/laser_acceleration/inputs_1d b/Examples/Physics_applications/laser_acceleration/inputs_1d index 8d92bfa356c..95e54c7d43e 100644 --- a/Examples/Physics_applications/laser_acceleration/inputs_1d +++ b/Examples/Physics_applications/laser_acceleration/inputs_1d @@ -51,7 +51,7 @@ electrons.addIntegerAttributes = regionofinterest electrons.attribute.regionofinterest(x,y,z,ux,uy,uz,t) = " (z>12.0e-6) * (z<13.0e-6)" ################################# -############ PLASMA ############# +############ LASER ############## ################################# lasers.names = laser1 laser1.profile = Gaussian diff --git a/Examples/Physics_applications/laser_acceleration/inputs_2d b/Examples/Physics_applications/laser_acceleration/inputs_2d index 422e33ddf57..edcb8aba279 100644 --- a/Examples/Physics_applications/laser_acceleration/inputs_2d +++ b/Examples/Physics_applications/laser_acceleration/inputs_2d @@ -68,7 +68,7 @@ beam.uy_th = 2. beam.uz_th = 50. ################################# -############ PLASMA ############# +############ LASER ############## ################################# lasers.names = laser1 laser1.profile = Gaussian diff --git a/Examples/Physics_applications/laser_acceleration/inputs_rz b/Examples/Physics_applications/laser_acceleration/inputs_rz index 971f1b538cd..44350882193 100644 --- a/Examples/Physics_applications/laser_acceleration/inputs_rz +++ b/Examples/Physics_applications/laser_acceleration/inputs_rz @@ -70,7 +70,7 @@ beam.uy_th = 2. beam.uz_th = 50. ################################# -############ PLASMA ############# +############ LASER ############## ################################# lasers.names = laser1 laser1.profile = Gaussian diff --git a/Examples/Tests/FieldProbe/analysis_field_probe.py b/Examples/Tests/FieldProbe/analysis_field_probe.py index 7d7eaa0eb3b..e038bfd0f48 100755 --- a/Examples/Tests/FieldProbe/analysis_field_probe.py +++ b/Examples/Tests/FieldProbe/analysis_field_probe.py @@ -1,6 +1,6 @@ #!/usr/bin/env python3 # -# Copyright 2021-2022 Tiberius Rheaume +# Copyright 2021-2022 Elisa Rheaume # # This file is part of WarpX. # diff --git a/Examples/Tests/FieldProbe/inputs_2d b/Examples/Tests/FieldProbe/inputs_2d index 848dd5adff8..eb0f427e86c 100644 --- a/Examples/Tests/FieldProbe/inputs_2d +++ b/Examples/Tests/FieldProbe/inputs_2d @@ -79,4 +79,4 @@ FP_line.x1_probe = 1.5e-6 FP_line.z1_probe = 1.7e-6 FP_line.resolution = 201 -authors = "Tiberius Rheaume , Axel Huebl " +authors = "Elisa Rheaume , Axel Huebl " diff --git a/LICENSE.txt b/LICENSE.txt index 6b92a5631be..3ccd4f4875a 100644 --- a/LICENSE.txt +++ b/LICENSE.txt @@ -1,4 +1,4 @@ -WarpX v22.09 Copyright (c) 2018-2022, The Regents of the University of California, through Lawrence Berkeley National Laboratory, and Lawrence Livermore National Security, LLC, for the operation of Lawrence Livermore National Laboratory (subject to receipt of any required approvals from the U.S. Dept. of Energy). All rights reserved. +WarpX v22.10 Copyright (c) 2018-2022, The Regents of the University of California, through Lawrence Berkeley National Laboratory, and Lawrence Livermore National Security, LLC, for the operation of Lawrence Livermore National Laboratory (subject to receipt of any required approvals from the U.S. Dept. of Energy). All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: diff --git a/Python/pywarpx/__init__.py b/Python/pywarpx/__init__.py index 5590c7d6fc5..8f49604563c 100644 --- a/Python/pywarpx/__init__.py +++ b/Python/pywarpx/__init__.py @@ -19,3 +19,6 @@ from .Particles import electrons, newspecies, particles, positrons, protons from .WarpX import warpx from ._libwarpx import libwarpx + +# This is a circulor import and must happen after the import of libwarpx +from . import picmi # isort:skip diff --git a/Python/pywarpx/_libwarpx.py b/Python/pywarpx/_libwarpx.py index 691288d0787..d8f7f077f7c 100755 --- a/Python/pywarpx/_libwarpx.py +++ b/Python/pywarpx/_libwarpx.py @@ -831,6 +831,8 @@ def get_particle_x(self, species_name, level=0): return [struct['x'] for struct in structs] elif self.geometry_dim == 'rz': return [struct['x']*np.cos(theta) for struct, theta in zip(structs, self.get_particle_theta(species_name))] + elif self.geometry_dim == '1d': + raise Exception('get_particle_x: There is no x coordinate with 1D Cartesian') def get_particle_y(self, species_name, level=0): ''' @@ -840,10 +842,12 @@ def get_particle_y(self, species_name, level=0): ''' structs = self.get_particle_structs(species_name, level) - if self.geometry_dim == '3d' or self.geometry_dim == '2d': + if self.geometry_dim == '3d': return [struct['y'] for struct in structs] elif self.geometry_dim == 'rz': return [struct['x']*np.sin(theta) for struct, theta in zip(structs, self.get_particle_theta(species_name))] + elif self.geometry_dim == '1d' or self.geometry_dim == '2d': + raise Exception('get_particle_y: There is no y coordinate with 1D or 2D Cartesian') def get_particle_r(self, species_name, level=0): ''' @@ -857,8 +861,8 @@ def get_particle_r(self, species_name, level=0): return [struct['x'] for struct in structs] elif self.geometry_dim == '3d': return [np.sqrt(struct['x']**2 + struct['y']**2) for struct in structs] - elif self.geometry_dim == '2d': - raise Exception('get_particle_r: There is no r coordinate with 2D Cartesian') + elif self.geometry_dim == '2d' or self.geometry_dim == '1d': + raise Exception('get_particle_r: There is no r coordinate with 1D or 2D Cartesian') def get_particle_z(self, species_name, level=0): ''' @@ -872,6 +876,8 @@ def get_particle_z(self, species_name, level=0): return [struct['z'] for struct in structs] elif self.geometry_dim == 'rz' or self.geometry_dim == '2d': return [struct['y'] for struct in structs] + elif self.geometry_dim == '1d': + return [struct['x'] for struct in structs] def get_particle_id(self, species_name, level=0): ''' @@ -946,8 +952,8 @@ def get_particle_theta(self, species_name, level=0): elif self.geometry_dim == '3d': structs = self.get_particle_structs(species_name, level) return [np.arctan2(struct['y'], struct['x']) for struct in structs] - elif self.geometry_dim == '2d': - raise Exception('get_particle_r: There is no theta coordinate with 2D Cartesian') + elif self.geometry_dim == '2d' or self.geometry_dim == '1d': + raise Exception('get_particle_theta: There is no theta coordinate with 1D or 2D Cartesian') def get_particle_comp_index(self, species_name, pid_name): ''' diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index a9a8f80af04..4640a2e2e4d 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -39,7 +39,84 @@ class constants: picmistandard.register_constants(constants) + class Species(picmistandard.PICMI_Species): + """ + See `Input Parameters `_ for more information. + + Parameters + ---------- + warpx_boost_adjust_transverse_positions: bool, default=False + Whether to adjust transverse positions when apply the boost + to the simulation frame + + warpx_self_fields_required_precision: float, default=1.e-11 + Relative precision on the electrostatic solver + (when using the relativistic solver) + + warpx_self_fields_absolute_tolerance: float, default=0. + Absolute precision on the electrostatic solver + (when using the relativistic solver) + + warpx_self_fields_max_iters: integer, default=200 + Maximum number of iterations for the electrostatic + solver for the species + + warpx_self_fields_verbosity: integer, default=2 + Level of verbosity for the electrostatic solver + + warpx_save_previous_position: bool, default=False + Whether to save the old particle positions + + warpx_do_not_deposit: bool, default=False + Whether or not to deposit the charge and current density for + for this species + + warpx_reflection_model_xlo: string, default='0.' + Expression (in terms of the velocity "v") specifying the probability + that the particle will reflect on the lower x boundary + + warpx_reflection_model_xhi: string, default='0.' + Expression (in terms of the velocity "v") specifying the probability + that the particle will reflect on the upper x boundary + + warpx_reflection_model_ylo: string, default='0.' + Expression (in terms of the velocity "v") specifying the probability + that the particle will reflect on the lower y boundary + + warpx_reflection_model_yhi: string, default='0.' + Expression (in terms of the velocity "v") specifying the probability + that the particle will reflect on the upper y boundary + + warpx_reflection_model_zlo: string, default='0.' + Expression (in terms of the velocity "v") specifying the probability + that the particle will reflect on the lower z boundary + + warpx_reflection_model_zhi: string, default='0.' + Expression (in terms of the velocity "v") specifying the probability + that the particle will reflect on the upper z boundary + + warpx_save_particles_at_xlo: bool, default=False + Whether to save particles lost at the lower x boundary + + warpx_save_particles_at_xhi: bool, default=False + Whether to save particles lost at the upper x boundary + + warpx_save_particles_at_ylo: bool, default=False + Whether to save particles lost at the lower y boundary + + warpx_save_particles_at_yhi: bool, default=False + Whether to save particles lost at the upper y boundary + + warpx_save_particles_at_zlo: bool, default=False + Whether to save particles lost at the lower z boundary + + warpx_save_particles_at_zhi: bool, default=False + Whether to save particles lost at the upper z boundary + + warpx_save_particles_at_eb: bool, default=False + Whether to save particles lost at the embedded boundary + """ def init(self, kw): if self.particle_type == 'electron': @@ -404,7 +481,42 @@ def initialize_inputs(self, solver): class CylindricalGrid(picmistandard.PICMI_CylindricalGrid): - """This assumes that WarpX was compiled with USE_RZ = TRUE + """ + This assumes that WarpX was compiled with USE_RZ = TRUE + + See `Input Parameters `_ for more information. + + Parameters + --------- + warpx_max_grid_size: integer, default=32 + Maximum block size in either direction + + warpx_max_grid_size_x: integer, optional + Maximum block size in radial direction + + warpx_max_grid_size_y: integer, optional + Maximum block size in longitudinal direction + + warpx_blocking_factor: integer, optional + Blocking factor (which controls the block size) + + warpx_blocking_factor_x: integer, optional + Blocking factor (which controls the block size) in the radial direction + + warpx_blocking_factor_y: integer, optional + Blocking factor (which controls the block size) in the longitudinal direction + + warpx_potential_lo_r: float, default=0. + Electrostatic potential on the lower radial boundary + + warpx_potential_hi_r: float, default=0. + Electrostatic potential on the upper radial boundary + + warpx_potential_lo_z: float, default=0. + Electrostatic potential on the lower longitudinal boundary + + warpx_potential_hi_z: float, default=0. + Electrostatic potential on the upper longitudinal boundary """ def init(self, kw): self.max_grid_size = kw.pop('warpx_max_grid_size', 32) @@ -472,6 +584,29 @@ def initialize_inputs(self): class Cartesian1DGrid(picmistandard.PICMI_Cartesian1DGrid): + """ + See `Input Parameters `_ for more information. + + Parameters + --------- + warpx_max_grid_size: integer, default=32 + Maximum block size in either direction + + warpx_max_grid_size_x: integer, optional + Maximum block size in longitudinal direction + + warpx_blocking_factor: integer, optional + Blocking factor (which controls the block size) + + warpx_blocking_factor_x: integer, optional + Blocking factor (which controls the block size) in the longitudinal direction + + warpx_potential_lo_z: float, default=0. + Electrostatic potential on the lower longitudinal boundary + + warpx_potential_hi_z: float, default=0. + Electrostatic potential on the upper longitudinal boundary + """ def init(self, kw): self.max_grid_size = kw.pop('warpx_max_grid_size', 32) self.max_grid_size_x = kw.pop('warpx_max_grid_size_x', None) @@ -525,6 +660,41 @@ def initialize_inputs(self): pywarpx.amr.max_level = 0 class Cartesian2DGrid(picmistandard.PICMI_Cartesian2DGrid): + """ + See `Input Parameters `_ for more information. + + Parameters + --------- + warpx_max_grid_size: integer, default=32 + Maximum block size in either direction + + warpx_max_grid_size_x: integer, optional + Maximum block size in x direction + + warpx_max_grid_size_y: integer, optional + Maximum block size in z direction + + warpx_blocking_factor: integer, optional + Blocking factor (which controls the block size) + + warpx_blocking_factor_x: integer, optional + Blocking factor (which controls the block size) in the x direction + + warpx_blocking_factor_y: integer, optional + Blocking factor (which controls the block size) in the z direction + + warpx_potential_lo_x: float, default=0. + Electrostatic potential on the lower x boundary + + warpx_potential_hi_x: float, default=0. + Electrostatic potential on the upper x boundary + + warpx_potential_lo_z: float, default=0. + Electrostatic potential on the lower z boundary + + warpx_potential_hi_z: float, default=0. + Electrostatic potential on the upper z boundary + """ def init(self, kw): self.max_grid_size = kw.pop('warpx_max_grid_size', 32) self.max_grid_size_x = kw.pop('warpx_max_grid_size_x', None) @@ -586,6 +756,53 @@ def initialize_inputs(self): class Cartesian3DGrid(picmistandard.PICMI_Cartesian3DGrid): + """ + See `Input Parameters `_ for more information. + + Parameters + --------- + warpx_max_grid_size: integer, default=32 + Maximum block size in either direction + + warpx_max_grid_size_x: integer, optional + Maximum block size in x direction + + warpx_max_grid_size_y: integer, optional + Maximum block size in z direction + + warpx_max_grid_size_z: integer, optional + Maximum block size in z direction + + warpx_blocking_factor: integer, optional + Blocking factor (which controls the block size) + + warpx_blocking_factor_x: integer, optional + Blocking factor (which controls the block size) in the x direction + + warpx_blocking_factor_y: integer, optional + Blocking factor (which controls the block size) in the z direction + + warpx_blocking_factor_z: integer, optional + Blocking factor (which controls the block size) in the z direction + + warpx_potential_lo_x: float, default=0. + Electrostatic potential on the lower x boundary + + warpx_potential_hi_x: float, default=0. + Electrostatic potential on the upper x boundary + + warpx_potential_lo_y: float, default=0. + Electrostatic potential on the lower z boundary + + warpx_potential_hi_y: float, default=0. + Electrostatic potential on the upper z boundary + + warpx_potential_lo_z: float, default=0. + Electrostatic potential on the lower z boundary + + warpx_potential_hi_z: float, default=0. + Electrostatic potential on the upper z boundary + """ def init(self, kw): self.max_grid_size = kw.pop('warpx_max_grid_size', 32) self.max_grid_size_x = kw.pop('warpx_max_grid_size_x', None) @@ -653,6 +870,38 @@ def initialize_inputs(self): pywarpx.amr.max_level = 0 class ElectromagneticSolver(picmistandard.PICMI_ElectromagneticSolver): + """ + See `Input Parameters `_ for more information. + + Parameters + ---------- + warpx_pml_ncell: integer, optional + The depth of the PML, in number of cells + + warpx_periodic_single_box_fft: bool, default=False + Whether to do the spectral solver FFTs assuming a single + simulation block + + warpx_current_correction: bool, default=True + Whether to do the current correction for the spectral solver. + See documentation for exceptions to the default value. + + warpx_psatd_update_with_rho: bool, optional + Whether to update with the actual rho for the spectral solver + + warpx_psatd_do_time_averaging: bool, optional + Whether to do the time averaging for the spectral solver + + warpx_do_pml_in_domain: bool, default=False + Whether to do the PML boundaries within the domain (versus + in the guard cells) + + warpx_pml_has_particles: bool, default=False + Whether to allow particles in the PML region + + warpx_do_pml_j_damping: bool, default=False + Whether to do damping of J in the PML + """ def init(self, kw): assert self.method is None or self.method in ['Yee', 'CKC', 'PSATD', 'ECT'], Exception("Only 'Yee', 'CKC', 'PSATD', and 'ECT' are supported") @@ -718,6 +967,20 @@ def initialize_inputs(self): pywarpx.warpx.do_pml_j_damping = self.do_pml_j_damping class ElectrostaticSolver(picmistandard.PICMI_ElectrostaticSolver): + """ + See `Input Parameters `_ for more information. + + Parameters + ---------- + warpx_relativistic: bool, default=False + Whether to use the relativistic solver or lab frame solver + + warpx_absolute_tolerance: float, default=0. + Absolute tolerance on the lab fram solver + + warpx_self_fields_verbosity: integer, default=2 + Level of verbosity for the lab frame solver + """ def init(self, kw): self.relativistic = kw.pop('warpx_relativistic', False) self.absolute_tolerance = kw.pop('warpx_absolute_tolerance', None) @@ -862,10 +1125,26 @@ def initialize_inputs(self): class CoulombCollisions(picmistandard.base._ClassWithInit): - """Custom class to handle setup of binary Coulmb collisions in WarpX. If + """ + Custom class to handle setup of binary Coulmb collisions in WarpX. If collision initialization is added to picmistandard this can be changed to - inherit that functionality.""" + inherit that functionality. + + Parameters + ---------- + name: string + Name of instance (used in the inputs file) + species: list of species instances + The species involved in the collision. Must be of length 2. + + CoulombLog: float, optional + Value of the Coulomb log to use in the collision cross section. + If not supplied, it is calculated from the local conditions. + + ndt: integer, optional + The collisions will be applied every "ndt" steps. Must be 1 or larger. + """ def __init__(self, name, species, CoulombLog=None, ndt=None, **kw): self.name = name self.species = species @@ -883,9 +1162,35 @@ def initialize_inputs(self): class MCCCollisions(picmistandard.base._ClassWithInit): - """Custom class to handle setup of MCC collisions in WarpX. If collision + """ + Custom class to handle setup of MCC collisions in WarpX. If collision initialization is added to picmistandard this can be changed to inherit - that functionality.""" + that functionality. + + Parameters + ---------- + name: string + Name of instance (used in the inputs file) + + species: species instance + The species involved in the collision + + background_density: float + The density of the background + + background_temperature: float + The temperature of the background + + scattering_processes: dictionary + The scattering process to use and any needed information + + background_mass: float, optional + The mass of the background particle. If not supplied, the default depends + on the type of scattering process. + + ndt: integer, optional + The collisions will be applied every "ndt" steps. Must be 1 or larger. + """ def __init__(self, name, species, background_density, background_temperature, scattering_processes, @@ -924,14 +1229,31 @@ class EmbeddedBoundary(picmistandard.base._ClassWithInit): changed to inherit that functionality. The geometry can be specified either as an implicit function or as an STL file (ASCII or binary). In the latter case the geometry specified in the STL file can be scaled, translated and inverted. - - implicit_function: Analytic expression describing the embedded boundary - - stl_file: STL file path (string), file contains the embedded boundary geometry - - stl_scale: factor by which the STL geometry is scaled (pure number) - - stl_center: vector by which the STL geometry is translated (in meters) - - stl_reverse_normal: if True inverts the orientation of the STL geometry - - potential: Analytic expression defining the potential. Can only be specified - when the solver is electrostatic. Optional, defaults to 0. - Parameters used in the expressions should be given as additional keyword arguments. + + Parameters + ---------- + implicit_function: string + Analytic expression describing the embedded boundary + + stl_file: string + STL file path (string), file contains the embedded boundary geometry + + stl_scale: float + Factor by which the STL geometry is scaled + + stl_center: vector of floats + Vector by which the STL geometry is translated (in meters) + + stl_reverse_normal: bool + If True inverts the orientation of the STL geometry + + potential: string, default=0. + Analytic expression defining the potential. Can only be specified + when the solver is electrostatic. + + + Parameters used in the analytic expressions should be given as additional keyword arguments. + """ def __init__(self, implicit_function=None, stl_file=None, stl_scale=None, stl_center=None, stl_reverse_normal=False, potential=None, **kw): @@ -990,11 +1312,36 @@ def initialize_inputs(self, solver): class PlasmaLens(picmistandard.base._ClassWithInit): """ Custom class to setup a plasma lens lattice. - The applied fields are dependent on the transverse position - - Ex = x*stengths_E - - Ey = y*stengths_E - - Bx = +y*stengths_B - - By = -x*stengths_B + The applied fields are dependent only on the transverse position. + + Parameters + ---------- + period: float + Periodicity of the lattice (in lab frame, in meters) + + starts: list of floats + The start of each lens relative to the periodic repeat + + lengths: list of floats + The length of each lens + + strengths_E=None: list of floats, default = 0. + The electric field strength of each lens + + strengths_B=None: list of floats, default = 0. + The magnetic field strength of each lens + + + The field that is applied depends on the transverse position of the particle, (x,y) + + - Ex = x*stengths_E + + - Ey = y*stengths_E + + - Bx = +y*stengths_B + + - By = -x*stengths_B + """ def __init__(self, period, starts, lengths, strengths_E=None, strengths_B=None, **kw): self.period = period @@ -1020,6 +1367,78 @@ def initialize_inputs(self): class Simulation(picmistandard.PICMI_Simulation): + """ + See `Input Parameters `_ for more information. + + Parameters + ---------- + warpx_current_deposition_algo: {'direct', 'esirkepov', and 'vay'}, optional + Current deposition algorithm. The default depends on conditions. + + warpx_charge_deposition_algo: {'standard'}, optional + Charge deposition algorithm. + + warpx_field_gathering_algo: {'energy-conserving', 'momentum-conserving'}, optional + Field gathering algorithm. The default depends on conditions. + + warpx_particle_pusher_algo: {'boris', 'vay', 'higuera'}, default='boris' + Particle pushing algorithm. + + warpx_use_filter: bool, optional + Whether to use filtering. The default depends on the conditions. + + warpx_serialize_initial_conditions: bool, default=False + Controls the random numbers used for initialization. + This parameter should only be used for testing and continuous integration. + + warpx_do_dynamic_scheduling: bool, default=True + Whether to do dynamic scheduling with OpenMP + + warpx_load_balance_intervals: string, default='0' + The intervals for doing load balancing + + warpx_load_balance_efficiency_ratio_threshold: float, default=1.1 + (See documentation) + + warpx_load_balance_with_sfc: bool, default=0 + (See documentation) + + warpx_load_balance_knapsack_factor: float, default=1.24 + (See documentation) + + warpx_load_balance_costs_update: {'heuristic' or 'timers' or 'gpuclock'}, optional + (See documentation) + + warpx_costs_heuristic_particles_wt: float, optional + (See documentation) + + warpx_costs_heuristic_cells_wt: float, optional + (See documentation) + + warpx_use_fdtd_nci_corr: bool, optional + Whether to use the NCI correction when using the FDTD solver + + warpx_amr_check_input: bool, optional + Whether AMReX should perform checks on the input + (primarily related to the max grid size and blocking factors) + + warpx_amr_restart: string, optional + The name of the restart to use + + warpx_zmax_plasma_to_compute_max_step: float, optional + Sets the simulation run time based on the maximum z value + + warpx_collisions: collision instance, optional + The collision instance specifying the particle collisions + + warpx_embedded_boundary: embedded boundary instance, optional + + warpx_break_signals: list of strings + Signals on which to break + + warpx_checkpoint_signals: list of strings + Signals on which to write out a checkpoint + """ # Set the C++ WarpX interface (see _libwarpx.LibWarpX) as an extension to # Simulation objects. In the future, LibWarpX objects may actually be owned @@ -1178,6 +1597,32 @@ def finalize(self): class FieldDiagnostic(picmistandard.PICMI_FieldDiagnostic): + """ + See `Input Parameters `_ for more information. + + Parameters + ---------- + warpx_plot_raw_fields: bool, optional + Flag whether to dump the raw fields + + warpx_plot_raw_fields_guards: bool, optional + Flag whether the raw fields should include the guard cells + + warpx_format: {plotfile, checkpoint, openpmd, ascent, sensei}, optional + Diagnostic file format + + warpx_openpmd_backend: {bp, h5, json}, optional + Openpmd backend file format + + warpx_file_prefix: string, optional + Prefix on the diagnostic file name + + warpx_file_min_digits: integer, optional + Minimum number of digits for the time step number in the file name + + warpx_dump_rz_modes: bool, optional + Flag whether to dump the data for all RZ modes + """ def init(self, kw): self.plot_raw_fields = kw.pop('warpx_plot_raw_fields', None) @@ -1274,6 +1719,20 @@ def initialize_inputs(self): class Checkpoint(picmistandard.base._ClassWithInit): + """ + Sets up checkpointing of the simulation, allowing for later restarts + + See `Input Parameters `_ for more information. + + Parameters + ---------- + warpx_file_prefix: string + The prefix to the checkpoint directory names + + warpx_file_min_digits: integer + Minimum number of digits for the time step number in the checkpoint + directory name. + """ def __init__(self, period = 1, write_dir = None, name = None, **kw): @@ -1307,6 +1766,32 @@ def initialize_inputs(self): self.diagnostic.file_prefix = os.path.join(write_dir, file_prefix) class ParticleDiagnostic(picmistandard.PICMI_ParticleDiagnostic): + """ + See `Input Parameters `_ for more information. + + Parameters + ---------- + warpx_format: {plotfile, checkpoint, openpmd, ascent, sensei}, optional + Diagnostic file format + + warpx_openpmd_backend: {bp, h5, json}, optional + Openpmd backend file format + + warpx_file_prefix: string, optional + Prefix on the diagnostic file name + + warpx_file_min_digits: integer, optional + Minimum number of digits for the time step number in the file name + + warpx_random_fraction: float, optional + Random fraction of particles to include in the diagnostic + + warpx_uniform_stride: integer, optional + Stride to down select to the particles to include in the diagnostic + + warpx_plot_filter_function: string, optional + Analytic expression to down select the particles to in the diagnostic + """ def init(self, kw): self.format = kw.pop('warpx_format', 'plotfile') @@ -1407,17 +1892,34 @@ def initialize_inputs(self): class LabFrameFieldDiagnostic(picmistandard.PICMI_LabFrameFieldDiagnostic): """ - Warp specific arguments: - - warpx_new_BTD: Use the new BTD diagnostics - - warpx_format: Passed to .format - - warpx_openpmd_backend: Passed to .openpmd_backend - - warpx_file_prefix: Passed to .file_prefix - - warpx_file_min_digits: Passed to .file_min_digits - - warpx_buffer_size: Passed to .buffer_size - - warpx_lower_bound: Passed to .lower_bound - - warpx_upper_bound: Passed to .upper_bound + See `Input Parameters `_ for more information. + + Parameters + ---------- + warpx_new_BTD: bool, optional + Use the new BTD diagnostics + + warpx_format: string, optional + Passed to .format + + warpx_openpmd_backend: string, optional + Passed to .openpmd_backend + + warpx_file_prefix: string, optional + Passed to .file_prefix + + warpx_file_min_digits: integer, optional + Passed to .file_min_digits + + warpx_buffer_size: integer, optional + Passed to .buffer_size + + warpx_lower_bound: vector of floats, optional + Passed to .lower_bound + + warpx_upper_bound: vector of floats, optional + Passed to .upper_bound """ - __doc__ = picmistandard.PICMI_LabFrameFieldDiagnostic.__doc__ + __doc__ def init(self, kw): self.use_new_BTD = kw.pop('warpx_new_BTD', False) if self.use_new_BTD: diff --git a/Python/setup.py b/Python/setup.py index f662a436777..83abf888e5c 100644 --- a/Python/setup.py +++ b/Python/setup.py @@ -54,12 +54,12 @@ package_data = {} setup(name = 'pywarpx', - version = '22.09', + version = '22.10', packages = ['pywarpx'], package_dir = {'pywarpx': 'pywarpx'}, description = """Wrapper of WarpX""", package_data = package_data, - install_requires = ['numpy', 'picmistandard==0.0.19', 'periodictable'], + install_requires = ['numpy', 'picmistandard==0.0.20', 'periodictable'], python_requires = '>=3.7', zip_safe=False ) diff --git a/Regression/Checksum/benchmark.py b/Regression/Checksum/benchmark.py index 771aa06e1d2..fbbe44f98b0 100644 --- a/Regression/Checksum/benchmark.py +++ b/Regression/Checksum/benchmark.py @@ -13,20 +13,25 @@ class Benchmark: - '''Holds data and functions for referenc benchmark of one checksum test. - ''' + """ + Holds data and functions for referenc benchmark of one checksum test. + """ def __init__(self, test_name, data=None): - '''Constructor - + """ + Benchmark constructor. Store test name and reference checksum value, either from benchmark (used for comparison) or from a plotfile (used to reset a benchmark). - @param self The object pointer. - @param test_name Name of test, as found between [] in .ini file. - @param data checksum value (dictionary). - If None, it is read from benchmark. - ''' + Parameters + ---------- + test_name: string + Name of test, as found between [] in .ini file. + + data: dictionary, optional + Checksum value. + If None, it is read from benchmark. + """ self.test_name = test_name self.json_file = os.path.join(config.benchmark_location, @@ -37,18 +42,16 @@ def __init__(self, test_name, data=None): self.data = data def reset(self): - '''Update the benchmark (overwrites reference json file). - - @param self The object pointer. - ''' + """ + Update the benchmark (overwrites reference json file). + """ with open(self.json_file, 'w') as outfile: json.dump(self.data, outfile, sort_keys=True, indent=2) def get(self): - '''Read benchmark from reference json file. - - @param self The object pointer. - ''' + """ + Read benchmark from reference json file. + """ with open(self.json_file) as infile: data = json.load(infile) diff --git a/Regression/Checksum/benchmarks_json/BTD_ReducedSliceDiag.json b/Regression/Checksum/benchmarks_json/BTD_ReducedSliceDiag.json new file mode 100644 index 00000000000..43b609e30d8 --- /dev/null +++ b/Regression/Checksum/benchmarks_json/BTD_ReducedSliceDiag.json @@ -0,0 +1,32 @@ +{ + "beam": { + "particle_momentum_x": 1.1691227280755267e-17, + "particle_momentum_y": 2.3913921202517482e-17, + "particle_momentum_z": 5.159955875755814e-14, + "particle_position_x": 0.0005606444019625324, + "particle_position_y": 0.0008259894145730366, + "particle_position_z": 2.980004880402436, + "particle_weight": 62415.090744607616 + }, + "electrons": { + "particle_momentum_x": 9.678379481090198e-20, + "particle_momentum_y": 2.5764224165967887e-19, + "particle_momentum_z": 3.0774186320322946e-19, + "particle_position_x": 0.0025853613792120563, + "particle_position_y": 0.0037728464704368673, + "particle_position_z": 0.1901073014926564, + "particle_weight": 1787508743647.5366 + }, + "lev=0": { + "Bx": 498944441.0100795, + "By": 243006013.74301538, + "Bz": 18073698.419804543, + "Ex": 7.004940326170131e+16, + "Ey": 1.4626998046065405e+17, + "Ez": 2.144506972885858e+16, + "jx": 8.583168660352884e+17, + "jy": 2.1430893961423186e+18, + "jz": 2.3459512897800802e+19, + "rho": 74967092858.04996 + } +} \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/Deuterium_Tritium_Fusion_3D.json b/Regression/Checksum/benchmarks_json/Deuterium_Tritium_Fusion_3D.json index eec472698da..2c640846857 100644 --- a/Regression/Checksum/benchmarks_json/Deuterium_Tritium_Fusion_3D.json +++ b/Regression/Checksum/benchmarks_json/Deuterium_Tritium_Fusion_3D.json @@ -74,4 +74,4 @@ "particle_position_z": 819255.8152412223, "particle_weight": 1.0239999999424347e+29 } -} \ No newline at end of file +} diff --git a/Regression/Checksum/benchmarks_json/Deuterium_Tritium_Fusion_RZ.json b/Regression/Checksum/benchmarks_json/Deuterium_Tritium_Fusion_RZ.json new file mode 100644 index 00000000000..a5136dda644 --- /dev/null +++ b/Regression/Checksum/benchmarks_json/Deuterium_Tritium_Fusion_RZ.json @@ -0,0 +1,77 @@ +{ + "deuterium1": { + "particle_momentum_x": 1.8388106511899905e-13, + "particle_momentum_y": 1.837868790009435e-13, + "particle_momentum_z": 0.0, + "particle_position_x": 40959919.499819286, + "particle_position_y": 81919224.48541151, + "particle_theta": 32166860.23003994, + "particle_weight": 3216.984554806547 + }, + "deuterium2": { + "particle_momentum_x": 0.0, + "particle_momentum_y": 0.0, + "particle_momentum_z": 3.336364094249911e-14, + "particle_position_x": 4095908.9083809257, + "particle_position_y": 8192069.080030457, + "particle_theta": 3216444.348910214, + "particle_weight": 3.1898417901971444e+29 + }, + "helium1": { + "particle_momentum_x": 1.858124399143442e-15, + "particle_momentum_y": 1.876715110797694e-15, + "particle_momentum_z": 1.7098432207359157e-15, + "particle_position_x": 152920.23233108618, + "particle_position_y": 323733.9138644398, + "particle_theta": 120064.13771707338, + "particle_weight": 1.603083276067953e-27 + }, + "helium2": { + "particle_momentum_x": 1.5195006688950936e-15, + "particle_momentum_y": 1.52430083815551e-15, + "particle_momentum_z": 1.7654865863613367e-15, + "particle_position_x": 136867.63803188328, + "particle_position_y": 286903.30393175944, + "particle_theta": 107912.20520382549, + "particle_weight": 2.0862696876352987e+19 + }, + "lev=0": { + "rho": 0.0 + }, + "neutron1": { + "particle_momentum_x": 1.7160671487712845e-15, + "particle_momentum_y": 1.7154753069055672e-15, + "particle_momentum_z": 1.7098432207359157e-15, + "particle_position_x": 152920.23233108618, + "particle_position_y": 323733.9138644398, + "particle_theta": 120064.13771707338, + "particle_weight": 1.603083276067953e-27 + }, + "neutron2": { + "particle_momentum_x": 1.5195006688950936e-15, + "particle_momentum_y": 1.52430083815551e-15, + "particle_momentum_z": 1.5463311225724366e-15, + "particle_position_x": 136867.63803188328, + "particle_position_y": 286903.30393175944, + "particle_theta": 107912.20520382549, + "particle_weight": 2.0862696876352987e+19 + }, + "tritium1": { + "particle_momentum_x": 1.8384658063720362e-13, + "particle_momentum_y": 1.8381593257898129e-13, + "particle_momentum_z": 0.0, + "particle_position_x": 40961278.052658774, + "particle_position_y": 81919046.8061561, + "particle_theta": 32163925.891884565, + "particle_weight": 3217.0912552970394 + }, + "tritium2": { + "particle_momentum_x": 0.0, + "particle_momentum_y": 0.0, + "particle_momentum_z": 0.0, + "particle_position_x": 409793.9651940968, + "particle_position_y": 819237.3558155322, + "particle_theta": 321974.4557387621, + "particle_weight": 3.218514276139388e+29 + } +} diff --git a/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_MR.json b/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_MR.json index a067edf6c7f..1fbf0e29f36 100644 --- a/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_MR.json +++ b/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_MR.json @@ -1,40 +1,40 @@ { "electrons": { - "particle_momentum_x": 4.245553180462133e-20, + "particle_momentum_x": 4.244331469283422e-20, "particle_momentum_y": 0.0, - "particle_momentum_z": 4.245553180462118e-20, - "particle_position_x": 0.6553607116882387, - "particle_position_y": 0.6553607116882386, + "particle_momentum_z": 4.244331469283411e-20, + "particle_position_x": 0.655360445407094, + "particle_position_y": 0.6553604454070939, "particle_weight": 3200000000000000.5 }, "lev=0": { "Bx": 0.0, - "By": 35.447497788711175, + "By": 34.880237093101385, "Bz": 0.0, - "Ex": 7569456532381.726, + "Ex": 7573738411730.345, "Ey": 0.0, - "Ez": 7569456532381.752, - "jx": 7305055873426885.0, + "Ez": 7573738411730.373, + "jx": 7301940679621768.0, "jy": 0.0, - "jz": 7305055873426942.0 + "jz": 7301940679621819.0 }, "lev=1": { "Bx": 0.0, - "By": 640.5385076158955, + "By": 663.9253192035617, "Bz": 0.0, - "Ex": 7591603030109.285, + "Ex": 7599937806859.2, "Ey": 0.0, - "Ez": 7591603030109.295, - "jx": 7225962509586529.0, + "Ez": 7599937806859.211, + "jx": 7111706209199148.0, "jy": 0.0, - "jz": 7225962509586578.0 + "jz": 7111706209199176.0 }, "positrons": { - "particle_momentum_x": 4.2453138721853093e-20, + "particle_momentum_x": 4.24387817954511e-20, "particle_momentum_y": 0.0, - "particle_momentum_z": 4.245313872185295e-20, - "particle_position_x": 0.6553599057578553, - "particle_position_y": 0.6553599057578553, + "particle_momentum_z": 4.243878179545099e-20, + "particle_position_x": 0.6553597738022557, + "particle_position_y": 0.6553597738022557, "particle_weight": 3200000000000000.5 } } \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_MR_anisotropic.json b/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_MR_anisotropic.json index 751568a0168..445f514723c 100644 --- a/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_MR_anisotropic.json +++ b/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_MR_anisotropic.json @@ -1,40 +1,40 @@ { "electrons": { - "particle_momentum_x": 4.241333654356747e-20, + "particle_momentum_x": 4.2405417073382216e-20, "particle_momentum_y": 0.0, - "particle_momentum_z": 4.243147160194918e-20, - "particle_position_x": 0.6553604355411431, - "particle_position_y": 0.6553604252756429, + "particle_momentum_z": 4.24319208902951e-20, + "particle_position_x": 0.6553602236954971, + "particle_position_y": 0.6553603460513587, "particle_weight": 3200000000000000.5 }, "lev=0": { "Bx": 0.0, - "By": 31.85950148151796, + "By": 29.016398075047086, "Bz": 0.0, - "Ex": 7575288819909.899, + "Ex": 7578443960945.405, "Ey": 0.0, - "Ez": 7572601899574.373, - "jx": 7297234592450874.0, + "Ez": 7572873581953.129, + "jx": 7295166572198360.0, "jy": 0.0, - "jz": 7303113293544192.0 + "jz": 7303381904136724.0 }, "lev=1": { "Bx": 0.0, - "By": 75.34189668301649, + "By": 72.99870260220474, "Bz": 0.0, - "Ex": 4601932334532.539, + "Ex": 4606426110840.904, "Ey": 0.0, - "Ez": 7011970862051.691, - "jx": 4494391278293099.5, + "Ez": 7017198453908.449, + "jx": 4489353143132240.0, "jy": 0.0, - "jz": 6848792132928508.0 + "jz": 6838171379723720.0 }, "positrons": { - "particle_momentum_x": 4.2408343046906606e-20, + "particle_momentum_x": 4.2397727290339454e-20, "particle_momentum_y": 0.0, - "particle_momentum_z": 4.243183197977201e-20, - "particle_position_x": 0.6553601964895763, - "particle_position_y": 0.6553595823177, + "particle_momentum_z": 4.2431934836877495e-20, + "particle_position_x": 0.6553599859378852, + "particle_position_y": 0.655359657073043, "particle_weight": 3200000000000000.5 } } \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/LaserAcceleration_single_precision_comms.json b/Regression/Checksum/benchmarks_json/LaserAcceleration_single_precision_comms.json index fce66bcd9bf..21b6da776b3 100644 --- a/Regression/Checksum/benchmarks_json/LaserAcceleration_single_precision_comms.json +++ b/Regression/Checksum/benchmarks_json/LaserAcceleration_single_precision_comms.json @@ -12,7 +12,7 @@ }, "lev=0": { "Bx": 5863879.02030842, - "By": 2411.495350685753, + "By": 2411.501904737579, "Bz": 116025.42462998998, "Ex": 6267728094111.663, "Ey": 1670763233105822.0, @@ -22,4 +22,4 @@ "jz": 1045267552192496.5, "rho": 2211742630.7074776 } -} \ No newline at end of file +} diff --git a/Regression/Checksum/benchmarks_json/background_mcc_dp_psp.json b/Regression/Checksum/benchmarks_json/background_mcc_dp_psp.json index 6029c19495a..3c54db63b95 100644 --- a/Regression/Checksum/benchmarks_json/background_mcc_dp_psp.json +++ b/Regression/Checksum/benchmarks_json/background_mcc_dp_psp.json @@ -1,22 +1,22 @@ { "electrons": { - "particle_momentum_x": 1.011437782317433e-18, - "particle_momentum_y": 2.818489864433302e-19, - "particle_momentum_z": 2.810005819942864e-19, - "particle_position_x": 17134.12344905036, - "particle_position_y": 935.6698553080840, + "particle_momentum_x": 1.0114377818220718e-18, + "particle_momentum_y": 2.8184898646509405e-19, + "particle_momentum_z": 2.810005821001351e-19, + "particle_position_x": 17134.123449403327, + "particle_position_y": 935.6698625905589, "particle_weight": 61112621534.71875 }, "he_ions": { - "particle_momentum_x": 2.882592880031919e-18, - "particle_momentum_y": 2.1960207748059422e-18, - "particle_momentum_z": 2.1982341415215343e-18, - "particle_position_x": 17605.83295166959, - "particle_position_y": 1099.9805173814, + "particle_momentum_x": 2.8825928799866215e-18, + "particle_momentum_y": 2.1960207749440838e-18, + "particle_momentum_z": 2.198234141518127e-18, + "particle_position_x": 17605.832952557976, + "particle_position_y": 1099.9805177808296, "particle_weight": 71973184795.40625 }, "lev=0": { - "rho_electrons": 0.03566096818084508, - "rho_he_ions": 0.04192385953761138 + "rho_electrons": 0.03566096817628402, + "rho_he_ions": 0.04192385953761396 } -} +} \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/particles_in_pml_2d_MR.json b/Regression/Checksum/benchmarks_json/particles_in_pml_2d_MR.json index 2042f0e2bd5..8f05b63e7f2 100644 --- a/Regression/Checksum/benchmarks_json/particles_in_pml_2d_MR.json +++ b/Regression/Checksum/benchmarks_json/particles_in_pml_2d_MR.json @@ -1,22 +1,22 @@ { "lev=0": { "Bx": 0.0, - "By": 3.578051588629886e-09, + "By": 3.5780515886298844e-09, "Bz": 0.0, "Ex": 1.9699822913484977, "Ey": 0.0, - "Ez": 0.5356004212513488, + "Ez": 0.5356004212513481, "jx": 0.0, "jy": 0.0, "jz": 0.0 }, "lev=1": { "Bx": 0.0, - "By": 3.2059933548393633e-09, + "By": 3.2059888629493914e-09, "Bz": 0.0, - "Ex": 2.459735202099343, + "Ex": 2.459734923669491, "Ey": 0.0, - "Ez": 0.4991264331012938, + "Ez": 0.4991257213214868, "jx": 0.0, "jy": 0.0, "jz": 0.0 diff --git a/Regression/Checksum/benchmarks_json/particles_in_pml_3d_MR.json b/Regression/Checksum/benchmarks_json/particles_in_pml_3d_MR.json index 43aee659baf..d665b04110f 100644 --- a/Regression/Checksum/benchmarks_json/particles_in_pml_3d_MR.json +++ b/Regression/Checksum/benchmarks_json/particles_in_pml_3d_MR.json @@ -1,24 +1,24 @@ { "lev=0": { - "Bx": 8.432558629746657e-05, - "By": 0.026892609338185998, - "Bz": 0.026892609345544216, - "Ex": 9112221.410831584, - "Ey": 4342060.779985666, - "Ez": 4342060.7789522, + "Bx": 8.432558629750022e-05, + "By": 0.02689260933254305, + "Bz": 0.026892609332543057, + "Ex": 9112221.407998262, + "Ey": 4342060.778180814, + "Ez": 4342060.778180817, "jx": 0.0, "jy": 0.0, "jz": 0.0 }, "lev=1": { - "Bx": 0.00013616524695631051, - "By": 0.04038118398020402, - "Bz": 0.040381183972079154, - "Ex": 13672614.170638269, - "Ey": 6988497.830608286, - "Ez": 6988497.829006851, + "Bx": 0.00013616513675437372, + "By": 0.040381141101323965, + "Bz": 0.040381141101323924, + "Ex": 13672597.083734166, + "Ey": 6988489.814361327, + "Ez": 6988489.814361326, "jx": 0.0, "jy": 0.0, "jz": 0.0 } -} +} \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/pml_psatd_dive_divb_cleaning.json b/Regression/Checksum/benchmarks_json/pml_psatd_dive_divb_cleaning.json index cfae22ad7d6..12f56218081 100644 --- a/Regression/Checksum/benchmarks_json/pml_psatd_dive_divb_cleaning.json +++ b/Regression/Checksum/benchmarks_json/pml_psatd_dive_divb_cleaning.json @@ -1,11 +1,11 @@ { "lev=0": { - "Bx": 1.2401714071102193e-06, - "By": 1.21931343936974e-06, - "Bz": 1.2293954248570514e-06, - "Ex": 614.2618449422512, - "Ey": 622.4788672892859, - "Ez": 613.5910967525642, - "rho": 4.90369639123331e-05 + "Bx": 1.482006352778953e-07, + "By": 1.48205157883426e-07, + "Bz": 1.4954704195856524e-07, + "Ex": 11.789793626679334, + "Ey": 11.78688532983594, + "Ez": 11.770112090435557, + "rho": 4.903696256562049e-05 } } \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/reduced_diags_single_precision.json b/Regression/Checksum/benchmarks_json/reduced_diags_single_precision.json index 670fb8b44e2..b2162221aba 100644 --- a/Regression/Checksum/benchmarks_json/reduced_diags_single_precision.json +++ b/Regression/Checksum/benchmarks_json/reduced_diags_single_precision.json @@ -1,43 +1,43 @@ { "electrons": { - "particle_momentum_x": 2.4417212503278005e-19, - "particle_momentum_y": 2.4564848534619294e-19, - "particle_momentum_z": 2.434157491249898e-19, - "particle_position_x": 16379.247699920481, - "particle_position_y": 16382.836164549284, - "particle_position_z": 16382.833915042822, + "particle_momentum_x": 2.441720869016668e-19, + "particle_momentum_y": 2.456481313629533e-19, + "particle_momentum_z": 2.434157770814882e-19, + "particle_position_x": 16379.247624470314, + "particle_position_y": 16382.836430045034, + "particle_position_z": 16382.83386090216, "particle_weight": 800000003014656.0 }, "lev=0": { - "Bx": 0.08424062892077444, - "By": 0.08475186723996586, - "Bz": 0.08796440505066272, - "Ex": 106872540.64410019, - "Ey": 107805501.93848419, - "Ez": 107783078.61432838, - "jx": 728220.4602749646, - "jy": 736403.7514596283, - "jz": 734813.9692403525, - "rho": 0.027750860941553768, - "rho_electrons": 0.525001227832945, - "rho_protons": 0.5250012272590538 + "Bx": 0.08424077006868202, + "By": 0.08475162318143958, + "Bz": 0.0879638695257583, + "Ex": 106872406.23206902, + "Ey": 107805796.05268478, + "Ez": 107783300.6896553, + "jx": 728221.4094335437, + "jy": 736401.7841227949, + "jz": 734813.4235675633, + "rho": 0.02775090560003779, + "rho_electrons": 0.5250012280648662, + "rho_protons": 0.5250012270062143 }, "photons": { "particle_momentum_x": 1.428291590249666e-18, "particle_momentum_y": 1.4222174024686332e-18, "particle_momentum_z": 1.4246585206989104e-18, - "particle_position_x": 16376.414799568243, - "particle_position_y": 16409.71378429825, - "particle_position_z": 16378.407592666219, + "particle_position_x": 16376.414787647314, + "particle_position_y": 16409.71380140478, + "particle_position_z": 16378.407592070173, "particle_weight": 800000003014656.0 }, "protons": { - "particle_momentum_x": 1.4104802628011656e-19, - "particle_momentum_y": 1.4120349291716893e-19, - "particle_momentum_z": 1.3903171427087353e-19, - "particle_position_x": 16383.951009619981, - "particle_position_y": 16383.99123002775, - "particle_position_z": 16384.033095447347, + "particle_momentum_x": 1.4104804792534928e-19, + "particle_momentum_y": 1.4120362404591128e-19, + "particle_momentum_z": 1.3903172753601875e-19, + "particle_position_x": 16383.951009282842, + "particle_position_y": 16383.991230854765, + "particle_position_z": 16384.03309576772, "particle_weight": 800000003014656.0 } } \ No newline at end of file diff --git a/Regression/Checksum/checksum.py b/Regression/Checksum/checksum.py index 472be941f8d..3e7ac78aa00 100644 --- a/Regression/Checksum/checksum.py +++ b/Regression/Checksum/checksum.py @@ -16,21 +16,29 @@ class Checksum: - '''Class for checksum comparison of one test. - ''' + """Class for checksum comparison of one test. + """ def __init__(self, test_name, plotfile, do_fields=True, do_particles=True): - '''Constructor - + """ + Checksum constructor. Store test_name and plotfile name, and compute checksum from plotfile and store it in self.data. - @param self The object pointer. - @param test_name Name of test, as found between [] in .ini file. - @param plotfile Plotfile from which the checksum is computed. - @param do_fields Whether to compare fields in the checksum. - @param do_particles Whether to compare particles in the checksum. - ''' + Parameters + ---------- + test_name: string + Name of test, as found between [] in .ini file. + + plotfile: string + Plotfile from which the checksum is computed. + + do_fields: bool, default=True + Whether to compare fields in the checksum. + + do_particles: bool, default=True + Whether to compare particles in the checksum. + """ self.test_name = test_name self.plotfile = plotfile @@ -38,16 +46,20 @@ def __init__(self, test_name, plotfile, do_fields=True, do_particles=True): do_particles=do_particles) def read_plotfile(self, do_fields=True, do_particles=True): - '''Get checksum from plotfile. - + """ + Get checksum from plotfile. Read an AMReX plotfile with yt, compute 1 checksum per field and return all checksums in a dictionary. The checksum of quantity Q is max(abs(Q)). - @param self The object pointer. - @param do_fields Whether to read fields from the plotfile. - @param do_particles Whether to read particles from the plotfile. - ''' + Parameters + ---------- + do_fields: bool, default=True + Whether to read fields from the plotfile. + + do_particles: bool, default=True + Whether to read particles from the plotfile. + """ ds = yt.load(self.plotfile) # yt 4.0+ has rounding issues with our domain data: @@ -103,17 +115,21 @@ def read_plotfile(self, do_fields=True, do_particles=True): return data def evaluate(self, rtol=1.e-9, atol=1.e-40): - '''Compare plotfile checksum with benchmark. - + """ + Compare plotfile checksum with benchmark. Read checksum from input plotfile, read benchmark corresponding to test_name, and assert that they are equal. Almost all the body of this functions is for user-readable print statements. - @param self The object pointer. - @param test_name Name of test, as found between [] in .ini file. - @param plotfile Plotfile from which the checksum is computed. - ''' + Parameters + ---------- + rtol: float, default=1.e-9 + Relative tolerance on the benchmark + + atol: float, default=1.e-40 + Absolute tolerance on the benchmark + """ ref_benchmark = Benchmark(self.test_name) diff --git a/Regression/Checksum/checksumAPI.py b/Regression/Checksum/checksumAPI.py index 3d25413bd41..c45b6e233f4 100755 --- a/Regression/Checksum/checksumAPI.py +++ b/Regression/Checksum/checksumAPI.py @@ -16,7 +16,7 @@ from benchmark import Benchmark from checksum import Checksum -''' +""" API for WarpX checksum tests. It can be used in two ways: - Directly use functions below to make a checksum test from a python script. @@ -32,39 +32,61 @@ * Reset a benchmark. From a bash terminal: $ ./checksumAPI.py --reset-benchmark --plotfile \ --test-name -''' +""" def evaluate_checksum(test_name, plotfile, rtol=1.e-9, atol=1.e-40, do_fields=True, do_particles=True): - '''Compare plotfile checksum with benchmark. - + """ + Compare plotfile checksum with benchmark. Read checksum from input plotfile, read benchmark corresponding to test_name, and assert their equality. - @param test_name Name of test, as found between [] in .ini file. - @param plotfile Plotfile from which the checksum is computed. - @param rtol Relative tolerance for the comparison. - @param atol Absolute tolerance for the comparison. - @param do_fields Whether to compare fields in the checksum. - @param do_particles Whether to compare particles in the checksum. - ''' + Parameters + ---------- + test_name: string + Name of test, as found between [] in .ini file. + + plotfile : string + Plotfile from which the checksum is computed. + + rtol: float, default=1.e-9 + Relative tolerance for the comparison. + + atol: float, default=1.e-40 + Absolute tolerance for the comparison. + + do_fields: bool, default=True + Whether to compare fields in the checksum. + + do_particles: bool, default=True + Whether to compare particles in the checksum. + """ test_checksum = Checksum(test_name, plotfile, do_fields=do_fields, do_particles=do_particles) test_checksum.evaluate(rtol=rtol, atol=atol) def reset_benchmark(test_name, plotfile, do_fields=True, do_particles=True): - '''Update the benchmark (overwrites reference json file). - + """ + Update the benchmark (overwrites reference json file). Overwrite value of benchmark corresponding to test_name with checksum read from input plotfile. - @param test_name Name of test, as found between [] in .ini file. - @param plotfile Plotfile from which the checksum is computed. - @param do_fields Whether to write field checksums in the benchmark. - @param do_particles Whether to write particles checksums in the benchmark. - ''' + Parameters + ---------- + test_name: string + Name of test, as found between [] in .ini file. + + plotfile: string + Plotfile from which the checksum is computed. + + do_fields: bool, default=True + Whether to write field checksums in the benchmark. + + do_particles: bool, default=True + Whether to write particles checksums in the benchmark. + """ ref_checksum = Checksum(test_name, plotfile, do_fields=do_fields, do_particles=do_particles) ref_benchmark = Benchmark(test_name, ref_checksum.data) @@ -72,13 +94,17 @@ def reset_benchmark(test_name, plotfile, do_fields=True, do_particles=True): def reset_all_benchmarks(path_to_all_plotfiles): - '''Update all benchmarks (overwrites reference json files) + """ + Update all benchmarks (overwrites reference json files) found in path_to_all_plotfiles - @param path_to_all_plotfiles Path to all plotfiles for which the benchmarks - are to be reset. The plotfiles should be named _plt, which is - what regression_testing.regtests.py does, provided we're careful enough. - ''' + Parameters + ---------- + path_to_all_plotfiles: string + Path to all plotfiles for which the benchmarks + are to be reset. The plotfiles should be named _plt, which is + what regression_testing.regtests.py does, provided we're careful enough. + """ # Get list of plotfiles in path_to_all_plotfiles plotfile_list = glob.glob(path_to_all_plotfiles + '*_plt*[0-9]', diff --git a/Regression/WarpX-GPU-tests.ini b/Regression/WarpX-GPU-tests.ini index 31fe16e1247..63b9f09dcbd 100644 --- a/Regression/WarpX-GPU-tests.ini +++ b/Regression/WarpX-GPU-tests.ini @@ -60,7 +60,7 @@ emailBody = Check https://ccse.lbl.gov/pub/GpuRegressionTesting/WarpX/ for more [AMReX] dir = /home/regtester/git/amrex/ -branch = 22.09 +branch = 13aa4df0f5a4af40270963ad5b42ac7ce662e045 [source] dir = /home/regtester/git/WarpX diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 8479b51ffd2..ba31ef94f5d 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -59,7 +59,7 @@ emailBody = Check https://ccse.lbl.gov/pub/RegressionTesting/WarpX/ for more det [AMReX] dir = /home/regtester/AMReX_RegTesting/amrex/ -branch = 22.09 +branch = 13aa4df0f5a4af40270963ad5b42ac7ce662e045 [source] dir = /home/regtester/AMReX_RegTesting/warpx @@ -414,7 +414,7 @@ analysisOutputImage = langmuir_multi_analysis.png [Langmuir_multi_psatd_multiJ] buildDir = . inputFile = Examples/Tests/Langmuir/inputs_3d_multi_rt -runtime_params = algo.maxwell_solver=psatd warpx.cfl=0.5773502691896258 algo.current_deposition=direct psatd.update_with_rho=1 warpx.do_multi_J=1 warpx.do_multi_J_n_depositions=2 warpx.abort_on_warning_threshold=medium +runtime_params = algo.maxwell_solver=psatd warpx.cfl=0.5773502691896258 algo.current_deposition=direct psatd.update_with_rho=1 warpx.do_multi_J=1 warpx.do_multi_J_n_depositions=2 psatd.J_in_time=linear warpx.abort_on_warning_threshold=medium dim = 3 addToCompileString = USE_PSATD=TRUE cmakeSetupOpts = -DWarpX_DIMS=3 -DWarpX_PSATD=ON @@ -433,7 +433,7 @@ analysisOutputImage = Langmuir_multi_psatd_multiJ.png [Langmuir_multi_psatd_multiJ_nodal] buildDir = . inputFile = Examples/Tests/Langmuir/inputs_3d_multi_rt -runtime_params = algo.maxwell_solver=psatd warpx.cfl=0.5773502691896258 algo.current_deposition=direct psatd.update_with_rho=1 warpx.do_multi_J=1 warpx.do_multi_J_n_depositions=2 warpx.abort_on_warning_threshold=medium warpx.do_nodal=1 +runtime_params = algo.maxwell_solver=psatd warpx.cfl=0.5773502691896258 algo.current_deposition=direct psatd.update_with_rho=1 warpx.do_multi_J=1 warpx.do_multi_J_n_depositions=2 psatd.J_in_time=linear warpx.abort_on_warning_threshold=medium warpx.do_nodal=1 dim = 3 addToCompileString = USE_PSATD=TRUE cmakeSetupOpts = -DWarpX_DIMS=3 -DWarpX_PSATD=ON @@ -699,7 +699,7 @@ analysisOutputImage = langmuir_multi_2d_analysis.png [Langmuir_multi_2d_psatd_multiJ] buildDir = . inputFile = Examples/Tests/Langmuir/inputs_2d_multi_rt -runtime_params = algo.maxwell_solver=psatd warpx.cfl=0.7071067811865475 psatd.update_with_rho=1 warpx.do_multi_J=1 warpx.do_multi_J_n_depositions=2 warpx.abort_on_warning_threshold=medium +runtime_params = algo.maxwell_solver=psatd warpx.cfl=0.7071067811865475 psatd.update_with_rho=1 warpx.do_multi_J=1 warpx.do_multi_J_n_depositions=2 psatd.J_in_time=linear warpx.abort_on_warning_threshold=medium dim = 2 addToCompileString = USE_PSATD=TRUE cmakeSetupOpts = -DWarpX_DIMS=2 -DWarpX_PSATD=ON @@ -718,7 +718,7 @@ analysisOutputImage = Langmuir_multi_2d_psatd_multiJ.png [Langmuir_multi_2d_psatd_multiJ_nodal] buildDir = . inputFile = Examples/Tests/Langmuir/inputs_2d_multi_rt -runtime_params = algo.maxwell_solver=psatd warpx.cfl=0.7071067811865475 psatd.update_with_rho=1 warpx.do_multi_J=1 warpx.do_multi_J_n_depositions=2 warpx.abort_on_warning_threshold=medium warpx.do_nodal=1 +runtime_params = algo.maxwell_solver=psatd warpx.cfl=0.7071067811865475 psatd.update_with_rho=1 warpx.do_multi_J=1 warpx.do_multi_J_n_depositions=2 psatd.J_in_time=linear warpx.abort_on_warning_threshold=medium warpx.do_nodal=1 dim = 2 addToCompileString = USE_PSATD=TRUE cmakeSetupOpts = -DWarpX_DIMS=2 -DWarpX_PSATD=ON @@ -2291,6 +2291,22 @@ compileTest = 0 doVis = 0 analysisRoutine = Examples/Modules/nuclear_fusion/analysis_deuterium_tritium_fusion.py +[Deuterium_Tritium_Fusion_RZ] +buildDir = . +inputFile = Examples/Modules/nuclear_fusion/inputs_deuterium_tritium_rz +runtime_params = warpx.do_dynamic_scheduling=0 warpx.serialize_initial_conditions=1 +dim = 2 +addToCompileString = USE_RZ=TRUE +cmakeSetupOpts = -DWarpX_DIMS=RZ +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 1 +compileTest = 0 +doVis = 0 +analysisRoutine = Examples/Modules/nuclear_fusion/analysis_deuterium_tritium_fusion.py + [Maxwell_Hybrid_QED_solver] buildDir = . inputFile = Examples/Tests/Maxwell_Hybrid_QED/inputs_2d @@ -2703,7 +2719,7 @@ analysisRoutine = Examples/Tests/galilean/analysis.py [multi_J_rz_psatd] buildDir = . inputFile = Examples/Tests/multi_J/inputs_rz -runtime_params = warpx.do_dynamic_scheduling=0 warpx.serialize_initial_conditions=1 warpx.abort_on_warning_threshold=medium +runtime_params = warpx.do_dynamic_scheduling=0 warpx.serialize_initial_conditions=1 warpx.abort_on_warning_threshold=medium psatd.J_in_time=linear dim = 2 addToCompileString = USE_RZ=TRUE USE_PSATD=TRUE cmakeSetupOpts = -DWarpX_DIMS=RZ -DWarpX_PSATD=ON @@ -3269,6 +3285,25 @@ analysisRoutine = Examples/Physics_applications/capacitive_discharge/analysis_2d buildDir = . inputFile = Examples/Physics_applications/capacitive_discharge/PICMI_inputs_1d.py runtime_params = +customRunCmd = python3 PICMI_inputs_1d.py --test --pythonsolver +dim = 1 +addToCompileString = USE_PYTHON_MAIN=TRUE USE_OPENPMD=TRUE QED=FALSE +cmakeSetupOpts = -DWarpX_DIMS=1 -DWarpX_APP=OFF -DWarpX_OPENPMD=ON -DWarpX_QED=OFF +target = pip_install +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 1 +compileTest = 0 +doVis = 0 +compareParticles = 0 +analysisRoutine = Examples/Physics_applications/capacitive_discharge/analysis_1d.py + +[Python_background_mcc_1d_tridiag] +buildDir = . +inputFile = Examples/Physics_applications/capacitive_discharge/PICMI_inputs_1d.py +runtime_params = customRunCmd = python3 PICMI_inputs_1d.py --test dim = 1 addToCompileString = USE_PYTHON_MAIN=TRUE USE_OPENPMD=TRUE QED=FALSE diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index 840f9d825bd..4688c19e508 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -131,7 +131,7 @@ public: int ncell, int delta, amrex::IntVect ref_ratio, amrex::Real dt, int nox_fft, int noy_fft, int noz_fft, bool do_nodal, int do_moving_window, int pml_has_particles, int do_pml_in_domain, - const bool do_multi_J, + const int J_in_time, const int rho_in_time, const bool do_pml_dive_cleaning, const bool do_pml_divb_cleaning, const amrex::IntVect& fill_guards_fields, const amrex::IntVect& fill_guards_current, diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 9c3c17b0bf7..712e287f0a7 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -548,7 +548,7 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri int ncell, int delta, amrex::IntVect ref_ratio, Real dt, int nox_fft, int noy_fft, int noz_fft, bool do_nodal, int do_moving_window, int /*pml_has_particles*/, int do_pml_in_domain, - const bool do_multi_J, + const int J_in_time, const int rho_in_time, const bool do_pml_dive_cleaning, const bool do_pml_divb_cleaning, const amrex::IntVect& fill_guards_fields, const amrex::IntVect& fill_guards_current, @@ -738,7 +738,7 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) { #ifndef WARPX_USE_PSATD - amrex::ignore_unused(lev, dt, do_multi_J); + amrex::ignore_unused(lev, dt, J_in_time, rho_in_time); # if(AMREX_SPACEDIM!=3) amrex::ignore_unused(noy_fft); # endif @@ -759,7 +759,7 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri spectral_solver_fp = std::make_unique(lev, realspace_ba, dm, nox_fft, noy_fft, noz_fft, do_nodal, v_galilean_zero, v_comoving_zero, dx, dt, in_pml, periodic_single_box, update_with_rho, - fft_do_time_averaging, do_multi_J, m_dive_cleaning, m_divb_cleaning); + fft_do_time_averaging, J_in_time, rho_in_time, m_dive_cleaning, m_divb_cleaning); #endif } @@ -878,7 +878,7 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri spectral_solver_cp = std::make_unique(lev, realspace_cba, cdm, nox_fft, noy_fft, noz_fft, do_nodal, v_galilean_zero, v_comoving_zero, cdx, dt, in_pml, periodic_single_box, update_with_rho, - fft_do_time_averaging, do_multi_J, m_dive_cleaning, m_divb_cleaning); + fft_do_time_averaging, J_in_time, rho_in_time, m_dive_cleaning, m_divb_cleaning); #endif } } diff --git a/Source/Diagnostics/BTDiagnostics.H b/Source/Diagnostics/BTDiagnostics.H index beb6c17ec21..fcfc8f13b17 100644 --- a/Source/Diagnostics/BTDiagnostics.H +++ b/Source/Diagnostics/BTDiagnostics.H @@ -10,6 +10,7 @@ #include "Diagnostics.H" #include "Diagnostics/ComputeDiagFunctors/ComputeDiagFunctor.H" #include "Utils/WarpXConst.H" +#include "Utils/IntervalsParser.H" #include #include @@ -37,6 +38,8 @@ private: bool m_plot_raw_fields_guards = false; /** Read relevant parameters for BTD */ void ReadParameters (); + /** Determines timesteps at which BTD diagnostics are written to file */ + BTDIntervalsParser m_intervals; /** \brief Flush m_mf_output and particles to file. * Currently, a temporary customized output format for the buffer * data is implemented and called in this function. diff --git a/Source/Diagnostics/BTDiagnostics.cpp b/Source/Diagnostics/BTDiagnostics.cpp index 72b6244e8ae..beb367309b6 100644 --- a/Source/Diagnostics/BTDiagnostics.cpp +++ b/Source/Diagnostics/BTDiagnostics.cpp @@ -112,9 +112,11 @@ void BTDiagnostics::DerivedInitData () // Turn on do_back_transformed_particles in the particle containers so that // the tmp_particle_data is allocated and the data of the corresponding species is // copied and stored in tmp_particle_data before particles are pushed. - for (auto const& species : m_output_species_names){ + if (m_do_back_transformed_particles) { mpc.SetDoBackTransformedParticles(m_do_back_transformed_particles); - mpc.SetDoBackTransformedParticles(species, m_do_back_transformed_particles); + for (auto const& species : m_output_species_names){ + mpc.SetDoBackTransformedParticles(species, m_do_back_transformed_particles); + } } m_particles_buffer.resize(m_num_buffers); m_totalParticles_flushed_already.resize(m_num_buffers); @@ -159,8 +161,18 @@ BTDiagnostics::ReadParameters () WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_do_back_transformed_fields, " fields must be turned on for the new back-transformed diagnostics"); if (m_do_back_transformed_fields == false) m_varnames.clear(); - getWithParser(pp_diag_name, "num_snapshots_lab", m_num_snapshots_lab); - m_num_buffers = m_num_snapshots_lab; + + std::vector intervals_string_vec = {"0"}; + bool const num_snapshots_specified = queryWithParser(pp_diag_name, "num_snapshots_lab", m_num_snapshots_lab); + bool const intervals_specified = pp_diag_name.queryarr("intervals", intervals_string_vec); + if (num_snapshots_specified) + { + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(!intervals_specified, + "For back-transformed diagnostics, user should specify either num_snapshots_lab or intervals, not both"); + intervals_string_vec = {":" + std::to_string(m_num_snapshots_lab-1)}; + } + m_intervals = BTDIntervalsParser(intervals_string_vec); + m_num_buffers = m_intervals.NumSnapshots(); // Read either dz_snapshots_lab or dt_snapshots_lab bool snapshot_interval_is_specified = false; @@ -239,7 +251,7 @@ BTDiagnostics::InitializeBufferData ( int i_buffer , int lev) auto & warpx = WarpX::GetInstance(); // Lab-frame time for the i^th snapshot amrex::Real zmax_0 = warpx.Geom(lev).ProbHi(m_moving_window_dir); - m_t_lab.at(i_buffer) = i_buffer * m_dt_snapshots_lab + m_t_lab.at(i_buffer) = m_intervals.GetBTDIteration(i_buffer) * m_dt_snapshots_lab + m_gamma_boost*m_beta_boost*zmax_0/PhysConst::c; // Define buffer domain in boosted frame at level, lev, with user-defined lo and hi @@ -1190,7 +1202,7 @@ void BTDiagnostics::UpdateTotalParticlesFlushed(int i_buffer) { for (int isp = 0; isp < m_totalParticles_flushed_already[i_buffer].size(); ++isp) { - m_totalParticles_flushed_already[i_buffer][isp] += m_totalParticles_in_buffer[i_buffer][isp]; + m_totalParticles_flushed_already[i_buffer][isp] += m_particles_buffer[i_buffer][isp]->TotalNumberOfParticles(); } } diff --git a/Source/Diagnostics/ReducedDiags/BeamRelevant.cpp b/Source/Diagnostics/ReducedDiags/BeamRelevant.cpp index 03ce90ac357..4b4d404ba8c 100644 --- a/Source/Diagnostics/ReducedDiags/BeamRelevant.cpp +++ b/Source/Diagnostics/ReducedDiags/BeamRelevant.cpp @@ -48,8 +48,9 @@ BeamRelevant::BeamRelevant (std::string rd_name) // 10,11,12: rms px,py,pz // 13: rms gamma // 14,15,16: emittance x,y,z - // 17: charge - m_data.resize(18, 0.0_rt); + // 17,18: beta-function x,y + // 19: charge + m_data.resize(20, 0.0_rt); #elif (defined WARPX_DIM_XZ) // 0, 1: mean x,z // 2, 3, 4: mean px,py,pz @@ -58,8 +59,9 @@ BeamRelevant::BeamRelevant (std::string rd_name) // 8, 9,10: rms px,py,pz // 11: rms gamma // 12,13: emittance x,z - // 14: charge - m_data.resize(15, 0.0_rt); + // 14: beta-function x + // 15: charge + m_data.resize(16, 0.0_rt); #elif (defined WARPX_DIM_1D_Z) // 0 : mean z // 1,2,3 : mean px,py,pz @@ -101,6 +103,8 @@ BeamRelevant::BeamRelevant (std::string rd_name) ofs << "[" << c++ << "]emittance_x(m)"; ofs << m_sep; ofs << "[" << c++ << "]emittance_y(m)"; ofs << m_sep; ofs << "[" << c++ << "]emittance_z(m)"; ofs << m_sep; + ofs << "[" << c++ << "]beta_x(m)"; ofs << m_sep; + ofs << "[" << c++ << "]beta_y(m)"; ofs << m_sep; ofs << "[" << c++ << "]charge(C)"; ofs << std::endl; #elif (defined WARPX_DIM_XZ) int c = 0; @@ -121,6 +125,7 @@ BeamRelevant::BeamRelevant (std::string rd_name) ofs << "[" << c++ << "]gamma_rms()"; ofs << m_sep; ofs << "[" << c++ << "]emittance_x(m)"; ofs << m_sep; ofs << "[" << c++ << "]emittance_z(m)"; ofs << m_sep; + ofs << "[" << c++ << "]beta_x(m)"; ofs << m_sep; ofs << "[" << c++ << "]charge(C)"; ofs << std::endl; #elif (defined WARPX_DIM_1D_Z) int c = 0; @@ -369,7 +374,9 @@ void BeamRelevant::ComputeDiags (int step) m_data[14] = std::sqrt(x_ms*ux_ms-xux*xux) / PhysConst::c; m_data[15] = std::sqrt(y_ms*uy_ms-yuy*yuy) / PhysConst::c; m_data[16] = std::sqrt(z_ms*uz_ms-zuz*zuz) / PhysConst::c; - m_data[17] = charge; + m_data[17] = (PhysConst::c * x_ms) / std::sqrt(x_ms*ux_ms-xux*xux); + m_data[18] = (PhysConst::c * y_ms) / std::sqrt(y_ms*uy_ms-yuy*yuy); + m_data[19] = charge; #elif (defined WARPX_DIM_XZ) m_data[0] = x_mean; m_data[1] = z_mean; @@ -385,7 +392,8 @@ void BeamRelevant::ComputeDiags (int step) m_data[11] = std::sqrt(gm_ms); m_data[12] = std::sqrt(x_ms*ux_ms-xux*xux) / PhysConst::c; m_data[13] = std::sqrt(z_ms*uz_ms-zuz*zuz) / PhysConst::c; - m_data[14] = charge; + m_data[14] = (PhysConst::c * x_ms) / std::sqrt(x_ms*ux_ms-xux*xux); + m_data[15] = charge; amrex::ignore_unused(y_mean, y_ms, yuy); #elif (defined WARPX_DIM_1D_Z) m_data[0] = z_mean; diff --git a/Source/Diagnostics/ReducedDiags/FieldProbe.H b/Source/Diagnostics/ReducedDiags/FieldProbe.H index 27c52618731..655ebc2a068 100644 --- a/Source/Diagnostics/ReducedDiags/FieldProbe.H +++ b/Source/Diagnostics/ReducedDiags/FieldProbe.H @@ -1,5 +1,5 @@ /* Copyright 2021 Lorenzo Giacomel, Neil Zaim, Yinjian Zhao - * Tiberius Rheaume, Axel Huebl + * Elisa Rheaume, Axel Huebl * * This file is part of WarpX. * diff --git a/Source/Diagnostics/ReducedDiags/FieldProbe.cpp b/Source/Diagnostics/ReducedDiags/FieldProbe.cpp index d5541b7ac19..bd5d5361106 100644 --- a/Source/Diagnostics/ReducedDiags/FieldProbe.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldProbe.cpp @@ -1,4 +1,4 @@ -/* Copyright 2021 Lorenzo Giacomel, Tiberius Rheaume, Axel Huebl +/* Copyright 2021 Lorenzo Giacomel, Elisa Rheaume, Axel Huebl * * This file is part of WarpX. * @@ -151,7 +151,7 @@ FieldProbe::FieldProbe (std::string rd_name) } pp_rd_name.query("integrate", m_field_probe_integrate); pp_rd_name.query("raw_fields", raw_fields); - pp_rd_name.query("interp_order", interp_order); + queryWithParser(pp_rd_name, "interp_order", interp_order); pp_rd_name.query("do_moving_window_FP", do_moving_window_FP); if (WarpX::gamma_boost > 1.0_rt) diff --git a/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.H b/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.H index 0b9ff481e05..dbb55fc0aff 100644 --- a/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.H +++ b/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.H @@ -1,4 +1,4 @@ -/* Copyright 2021 Tiberius Rheaume, Axel Huebl +/* Copyright 2021 Elisa Rheaume, Axel Huebl * * This file is part of WarpX. * diff --git a/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.cpp b/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.cpp index b1f20c4e76d..a49e1e08eaa 100644 --- a/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.cpp @@ -1,4 +1,4 @@ -/* Copyright 2021 Tiberius Rheaume, Axel Huebl +/* Copyright 2021 Elisa Rheaume, Axel Huebl * * This file is part of WarpX. * diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index cc71e78d5c3..202e91643bc 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -564,16 +564,19 @@ WarpX::OneStep_multiJ (const amrex::Real cur_time) // 4) Deposit J at relative time -dt with time step dt // (dt[0] denotes the time step on mesh refinement level 0) - auto& current = (WarpX::do_current_centering) ? current_fp_nodal : current_fp; - mypc->DepositCurrent(current, dt[0], -dt[0]); - // Synchronize J: filter, exchange boundary, and interpolate across levels. - // With current centering, the nodal current is deposited in 'current', - // namely 'current_fp_nodal': SyncCurrent stores the result of its centering - // into 'current_fp' and then performs both filtering, if used, and exchange - // of guard cells. - SyncCurrent(current_fp, current_cp); - // Forward FFT of J - PSATDForwardTransformJ(current_fp, current_cp); + if (J_in_time == JInTime::Linear) + { + auto& current = (WarpX::do_current_centering) ? current_fp_nodal : current_fp; + mypc->DepositCurrent(current, dt[0], -dt[0]); + // Synchronize J: filter, exchange boundary, and interpolate across levels. + // With current centering, the nodal current is deposited in 'current', + // namely 'current_fp_nodal': SyncCurrent stores the result of its centering + // into 'current_fp' and then performs both filtering, if used, and exchange + // of guard cells. + SyncCurrent(current_fp, current_cp); + // Forward FFT of J + PSATDForwardTransformJ(current_fp, current_cp); + } // Number of depositions for multi-J scheme const int n_depose = WarpX::do_multi_J_n_depositions; @@ -587,13 +590,21 @@ WarpX::OneStep_multiJ (const amrex::Real cur_time) for (int i_depose = 0; i_depose < n_loop; i_depose++) { // Move J deposited previously, from new to old - PSATDMoveJNewToJOld(); + if (J_in_time == JInTime::Linear) + { + PSATDMoveJNewToJOld(); + } + + const amrex::Real t_depose_current = (J_in_time == JInTime::Linear) ? + (i_depose-n_depose+1)*sub_dt : (i_depose-n_depose+0.5_rt)*sub_dt; - const amrex::Real t_depose = (i_depose-n_depose+1)*sub_dt; + // TODO Update this when rho quadratic in time is implemented + const amrex::Real t_depose_charge = (i_depose-n_depose+1)*sub_dt; - // Deposit new J at relative time t_depose with time step dt + // Deposit new J at relative time t_depose_current with time step dt // (dt[0] denotes the time step on mesh refinement level 0) - mypc->DepositCurrent(current, dt[0], t_depose); + auto& current = (WarpX::do_current_centering) ? current_fp_nodal : current_fp; + mypc->DepositCurrent(current, dt[0], t_depose_current); // Synchronize J: filter, exchange boundary, and interpolate across levels. // With current centering, the nodal current is deposited in 'current', // namely 'current_fp_nodal': SyncCurrent stores the result of its centering @@ -609,8 +620,8 @@ WarpX::OneStep_multiJ (const amrex::Real cur_time) // Move rho deposited previously, from new to old PSATDMoveRhoNewToRhoOld(); - // Deposit rho at relative time t_depose - mypc->DepositCharge(rho_fp, t_depose); + // Deposit rho at relative time t_depose_charge + mypc->DepositCharge(rho_fp, t_depose_charge); // Filter, exchange boundary, and interpolate across levels SyncRho(); // Forward FFT of rho_new @@ -716,7 +727,8 @@ WarpX::OneStep_sub1 (Real curtime) PushParticlesandDepose(fine_lev, curtime, DtType::FirstHalf); RestrictCurrentFromFineToCoarsePatch(current_fp, current_cp, fine_lev); RestrictRhoFromFineToCoarsePatch(rho_fp, rho_cp, fine_lev); - ApplyFilterandSumBoundaryJ(current_fp, current_cp, fine_lev, PatchType::fine); + if (use_filter) ApplyFilterJ(current_fp, fine_lev); + SumBoundaryJ(current_fp, fine_lev, Geom(fine_lev).periodicity()); ApplyFilterandSumBoundaryRho(rho_fp, rho_cp, fine_lev, PatchType::fine, 0, 2*ncomps); EvolveB(fine_lev, PatchType::fine, 0.5_rt*dt[fine_lev], DtType::FirstHalf); @@ -772,7 +784,8 @@ WarpX::OneStep_sub1 (Real curtime) PushParticlesandDepose(fine_lev, curtime+dt[fine_lev], DtType::SecondHalf); RestrictCurrentFromFineToCoarsePatch(current_fp, current_cp, fine_lev); RestrictRhoFromFineToCoarsePatch(rho_fp, rho_cp, fine_lev); - ApplyFilterandSumBoundaryJ(current_fp, current_cp, fine_lev, PatchType::fine); + if (use_filter) ApplyFilterJ(current_fp, fine_lev); + SumBoundaryJ(current_fp, fine_lev, Geom(fine_lev).periodicity()); ApplyFilterandSumBoundaryRho(rho_fp, rho_cp, fine_lev, PatchType::fine, 0, ncomps); EvolveB(fine_lev, PatchType::fine, 0.5_rt*dt[fine_lev], DtType::FirstHalf); diff --git a/Source/FieldSolver/ElectrostaticSolver.cpp b/Source/FieldSolver/ElectrostaticSolver.cpp index e50098e7d22..3a16b8b70f8 100644 --- a/Source/FieldSolver/ElectrostaticSolver.cpp +++ b/Source/FieldSolver/ElectrostaticSolver.cpp @@ -242,10 +242,24 @@ WarpX::AddSpaceChargeFieldLabFrame () setPhiBC(phi_fp); // Compute the potential phi, by solving the Poisson equation - if ( IsPythonCallBackInstalled("poissonsolver") ) ExecutePythonCallback("poissonsolver"); - else computePhi( rho_fp, phi_fp, beta, self_fields_required_precision, - self_fields_absolute_tolerance, self_fields_max_iters, - self_fields_verbosity ); + if (IsPythonCallBackInstalled("poissonsolver")) { + + // Use the Python level solver (user specified) + ExecutePythonCallback("poissonsolver"); + + } else { + +#if defined(WARPX_DIM_1D_Z) + // Use the tridiag solver with 1D + computePhiTriDiagonal(rho_fp, phi_fp); +#else + // Use the AMREX MLMG solver otherwise + computePhi(rho_fp, phi_fp, beta, self_fields_required_precision, + self_fields_absolute_tolerance, self_fields_max_iters, + self_fields_verbosity); +#endif + + } // Compute the electric field. Note that if an EB is used the electric // field will be calculated in the computePhi call. @@ -671,6 +685,190 @@ WarpX::computeB (amrex::Vector, 3> > } } +/* \brief Compute the potential by solving Poisson's equation with + a 1D tridiagonal solve. + + \param[in] rho The charge density a given species + \param[out] phi The potential to be computed by this function +*/ +void +WarpX::computePhiTriDiagonal (const amrex::Vector >& rho, + amrex::Vector >& phi) const +{ + + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(max_level == 0, + "The tridiagonal solver cannot be used with mesh refinement"); + + const int lev = 0; + + const amrex::Real* dx = Geom(lev).CellSize(); + const amrex::Real xmin = Geom(lev).ProbLo(0); + const amrex::Real xmax = Geom(lev).ProbHi(0); + const int nx_full_domain = static_cast( (xmax - xmin)/dx[0] + 0.5_rt ); + + int nx_solve_min = 1; + int nx_solve_max = nx_full_domain - 1; + + auto field_boundary_lo0 = WarpX::field_boundary_lo[0]; + auto field_boundary_hi0 = WarpX::field_boundary_hi[0]; + if (field_boundary_lo0 == FieldBoundaryType::None || field_boundary_lo0 == FieldBoundaryType::Periodic) { + // Neumann or periodic boundary condition + // Solve for the point on the lower boundary + nx_solve_min = 0; + } + if (field_boundary_hi0 == FieldBoundaryType::None || field_boundary_hi0 == FieldBoundaryType::Periodic) { + // Neumann or periodic boundary condition + // Solve for the point on the upper boundary + nx_solve_max = nx_full_domain; + } + + // Create a 1-D MultiFab that covers all of x, including guard cells on each end. + // The tridiag solve will be done in this MultiFab and then copied out afterwards. + const amrex::IntVect lo_full_domain(AMREX_D_DECL(-1,0,0)); + const amrex::IntVect hi_full_domain(AMREX_D_DECL(nx_full_domain+1,0,0)); + const amrex::Box box_full_domain(lo_full_domain, hi_full_domain); + const BoxArray ba_full_domain(box_full_domain); + amrex::DistributionMapping dm_full_domain; + amrex::Vector pmap = {0}; // The data will only be on processor 0 + dm_full_domain.define(pmap); + const int ncomps1d = 1; + const amrex::IntVect nguard1d(AMREX_D_DECL(1,0,0)); + const BoxArray ba_full_domain_node = amrex::convert(ba_full_domain, amrex::IntVect::TheNodeVector()); + + // Put the data in the pinned arena since the tridiag solver will be done on the CPU, but have + // the data readily accessible from the GPU. + auto phi1d_mf = MultiFab(ba_full_domain_node, dm_full_domain, ncomps1d, nguard1d, MFInfo().SetArena(The_Pinned_Arena())); + auto zwork1d_mf = MultiFab(ba_full_domain_node, dm_full_domain, ncomps1d, nguard1d, MFInfo().SetArena(The_Pinned_Arena())); + auto rho1d_mf = MultiFab(ba_full_domain_node, dm_full_domain, ncomps1d, nguard1d, MFInfo().SetArena(The_Pinned_Arena())); + + // Copy previous phi to get the boundary values + phi1d_mf.ParallelCopy(*phi[lev], 0, 0, 1, Geom(lev).periodicity()); + rho1d_mf.ParallelCopy(*rho[lev], 0, 0, 1, Geom(lev).periodicity()); + + // Multiplier on the charge density + const amrex::Real norm = dx[0]*dx[0]/PhysConst::ep0; + rho1d_mf.mult(norm); + + // Use the MFIter loop since when parallel, only process zero has a FAB. + // This skips the loop on all other processors. + for (MFIter mfi(phi1d_mf); mfi.isValid(); ++mfi) { + + const auto& phi1d_arr = phi1d_mf[mfi].array(); + const auto& zwork1d_arr = zwork1d_mf[mfi].array(); + const auto& rho1d_arr = rho1d_mf[mfi].array(); + + // The loops are always performed on the CPU + + amrex::Real diag = 2._rt; + + // The initial values depend on the boundary condition + if (field_boundary_lo0 == FieldBoundaryType::PEC) { + + phi1d_arr(1,0,0) = (phi1d_arr(0,0,0) + rho1d_arr(1,0,0))/diag; + + } else if (field_boundary_lo0 == FieldBoundaryType::None) { + + // Neumann boundary condition + phi1d_arr(0,0,0) = rho1d_arr(0,0,0)/diag; + + zwork1d_arr(1,0,0) = 2._rt/diag; + diag = 2._rt - zwork1d_arr(1,0,0); + phi1d_arr(1,0,0) = (rho1d_arr(1,0,0) - (-1._rt)*phi1d_arr(1-1,0,0))/diag; + + } else if (field_boundary_lo0 == FieldBoundaryType::Periodic) { + + phi1d_arr(0,0,0) = rho1d_arr(0,0,0)/diag; + + zwork1d_arr(1,0,0) = 1._rt/diag; + diag = 2._rt - zwork1d_arr(1,0,0); + phi1d_arr(1,0,0) = (rho1d_arr(1,0,0) - (-1._rt)*phi1d_arr(1-1,0,0))/diag; + + } + + // Loop upward, calculating the Gaussian elimination multipliers and right hand sides + for (int i_up = 2 ; i_up < nx_solve_max ; i_up++) { + + zwork1d_arr(i_up,0,0) = 1._rt/diag; + diag = 2._rt - zwork1d_arr(i_up,0,0); + phi1d_arr(i_up,0,0) = (rho1d_arr(i_up,0,0) - (-1._rt)*phi1d_arr(i_up-1,0,0))/diag; + + } + + // The last value depend on the boundary condition + int const imax = nx_solve_max; + amrex::Real zwork_product = 1.; // Needed for parallel boundaries + if (field_boundary_hi0 == FieldBoundaryType::PEC) { + + zwork1d_arr(imax,0,0) = 1._rt/diag; + diag = 2._rt - zwork1d_arr(imax,0,0); + phi1d_arr(imax,0,0) = (phi1d_arr(imax+1,0,0) + rho1d_arr(imax,0,0) - (-1._rt)*phi1d_arr(imax-1,0,0))/diag; + + } else if (field_boundary_hi0 == FieldBoundaryType::None) { + + // Neumann boundary condition + zwork1d_arr(imax,0,0) = 1._rt/diag; + diag = 2._rt - 2._rt*zwork1d_arr(imax,0,0); + if (diag == 0._rt) { + // This happens if the lower boundary is also Neumann. + // It this case, the potential is relative to an arbitrary constant, + // so set the upper boundary to zero to force a value. + phi1d_arr(imax,0,0) = 0.; + } else { + phi1d_arr(imax,0,0) = (rho1d_arr(imax,0,0) - (-1._rt)*phi1d_arr(imax-1,0,0))/diag; + } + + } else if (field_boundary_hi0 == FieldBoundaryType::Periodic) { + + zwork1d_arr(imax,0,0) = 1._rt/diag; + + for (int i = 1 ; i <= nx_solve_max ; i++) { + zwork_product *= zwork1d_arr(i,0,0); + } + + diag = 2._rt - zwork1d_arr(imax,0,0) - zwork_product; + phi1d_arr(imax,0,0) = (rho1d_arr(imax,0,0) - (-1._rt)*phi1d_arr(imax-1,0,0))/diag; + + } + + // Loop downward to calculate the phi + if (field_boundary_lo0 == FieldBoundaryType::Periodic) { + + // With periodic, the right hand column adds an extra term for all rows + for (int i_down = nx_solve_max-1 ; i_down >= nx_solve_min ; i_down--) { + zwork_product /= zwork1d_arr(i_down+1,0,0); + phi1d_arr(i_down,0,0) = phi1d_arr(i_down,0,0) + zwork1d_arr(i_down+1,0,0)*phi1d_arr(i_down+1,0,0) + zwork_product*phi1d_arr(imax,0,0); + } + + } else { + + for (int i_down = nx_solve_max-1 ; i_down >= nx_solve_min ; i_down--) { + phi1d_arr(i_down,0,0) = phi1d_arr(i_down,0,0) + zwork1d_arr(i_down+1,0,0)*phi1d_arr(i_down+1,0,0); + } + + } + + // Set the value in the guard cells + // The periodic case is handled in the ParallelCopy below + if (field_boundary_lo0 == FieldBoundaryType::PEC) { + phi1d_arr(-1,0,0) = phi1d_arr(0,0,0); + } else if (field_boundary_lo0 == FieldBoundaryType::None) { + phi1d_arr(-1,0,0) = phi1d_arr(1,0,0); + } + + if (field_boundary_hi0 == FieldBoundaryType::PEC) { + phi1d_arr(nx_full_domain+1,0,0) = phi1d_arr(nx_full_domain,0,0); + } else if (field_boundary_hi0 == FieldBoundaryType::None) { + phi1d_arr(nx_full_domain+1,0,0) = phi1d_arr(nx_full_domain-1,0,0); + } + + } + + // Copy phi1d to phi, including the x guard cell + const IntVect xghost(AMREX_D_DECL(1,0,0)); + phi[lev]->ParallelCopy(phi1d_mf, 0, 0, 1, xghost, xghost, Geom(lev).periodicity()); + +} + void ElectrostaticSolver::PoissonBoundaryHandler::definePhiBCs ( ) { int dim_start = 0; diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt index eb88b1f4e3c..a370d4b2d8d 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt @@ -1,6 +1,6 @@ target_sources(WarpX PRIVATE - PsatdAlgorithm.cpp + PsatdAlgorithmJConstantInTime.cpp PsatdAlgorithmJLinearInTime.cpp PsatdAlgorithmPml.cpp SpectralBaseAlgorithm.cpp diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package index 7c2bf428109..c798ffb01f5 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package @@ -1,5 +1,5 @@ CEXE_sources += SpectralBaseAlgorithm.cpp -CEXE_sources += PsatdAlgorithm.cpp +CEXE_sources += PsatdAlgorithmJConstantInTime.cpp CEXE_sources += PsatdAlgorithmJLinearInTime.cpp CEXE_sources += PsatdAlgorithmPml.cpp CEXE_sources += PsatdAlgorithmComoving.cpp diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.H similarity index 94% rename from Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H rename to Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.H index dd9c6a7fd37..0d2d67434b6 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.H @@ -4,8 +4,8 @@ * * License: BSD-3-Clause-LBNL */ -#ifndef WARPX_PSATD_ALGORITHM_H_ -#define WARPX_PSATD_ALGORITHM_H_ +#ifndef WARPX_PSATD_ALGORITHM_J_CONSTANT_IN_TIME_H_ +#define WARPX_PSATD_ALGORITHM_J_CONSTANT_IN_TIME_H_ #include "FieldSolver/SpectralSolver/SpectralFieldData.H" #include "FieldSolver/SpectralSolver/SpectralKSpace.H" @@ -24,12 +24,12 @@ /* \brief Class that updates the field in spectral space * and stores the coefficients of the corresponding update equation. */ -class PsatdAlgorithm : public SpectralBaseAlgorithm +class PsatdAlgorithmJConstantInTime : public SpectralBaseAlgorithm { public: /** - * \brief Constructor of the class PsatdAlgorithm + * \brief Constructor of the class PsatdAlgorithmJConstantInTime * * \param[in] spectral_kspace spectral space * \param[in] dm distribution mapping @@ -45,7 +45,7 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm * \param[in] dive_cleaning Update F as part of the field update, so that errors in divE=rho propagate away at the speed of light * \param[in] divb_cleaning Update G as part of the field update, so that errors in divB=0 propagate away at the speed of light */ - PsatdAlgorithm ( + PsatdAlgorithmJConstantInTime ( const SpectralKSpace& spectral_kspace, const amrex::DistributionMapping& dm, const SpectralFieldIndex& spectral_index, @@ -142,4 +142,4 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm bool m_is_galilean; }; #endif // WARPX_USE_PSATD -#endif // WARPX_PSATD_ALGORITHM_H_ +#endif // WARPX_PSATD_ALGORITHM_J_CONSTANT_IN_TIME_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.cpp similarity index 98% rename from Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp rename to Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.cpp index 1cbc27f0b1e..8971061f6ce 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.cpp @@ -4,7 +4,7 @@ * * License: BSD-3-Clause-LBNL */ -#include "PsatdAlgorithm.H" +#include "PsatdAlgorithmJConstantInTime.H" #include "Utils/TextMsg.H" #include "Utils/WarpXConst.H" @@ -27,7 +27,7 @@ using namespace amrex; -PsatdAlgorithm::PsatdAlgorithm( +PsatdAlgorithmJConstantInTime::PsatdAlgorithmJConstantInTime( const SpectralKSpace& spectral_kspace, const DistributionMapping& dm, const SpectralFieldIndex& spectral_index, @@ -110,7 +110,7 @@ PsatdAlgorithm::PsatdAlgorithm( } void -PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const +PsatdAlgorithmJConstantInTime::pushSpectralFields (SpectralFieldData& f) const { const bool update_with_rho = m_update_with_rho; const bool time_averaging = m_time_averaging; @@ -340,7 +340,7 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const } } -void PsatdAlgorithm::InitializeSpectralCoefficients ( +void PsatdAlgorithmJConstantInTime::InitializeSpectralCoefficients ( const SpectralKSpace& spectral_kspace, const amrex::DistributionMapping& dm, const amrex::Real dt) @@ -542,7 +542,7 @@ void PsatdAlgorithm::InitializeSpectralCoefficients ( } } -void PsatdAlgorithm::InitializeSpectralCoefficientsAveraging ( +void PsatdAlgorithmJConstantInTime::InitializeSpectralCoefficientsAveraging ( const SpectralKSpace& spectral_kspace, const amrex::DistributionMapping& dm, const amrex::Real dt) @@ -733,10 +733,10 @@ void PsatdAlgorithm::InitializeSpectralCoefficientsAveraging ( } } -void PsatdAlgorithm::CurrentCorrection (SpectralFieldData& field_data) +void PsatdAlgorithmJConstantInTime::CurrentCorrection (SpectralFieldData& field_data) { // Profiling - BL_PROFILE("PsatdAlgorithm::CurrentCorrection"); + BL_PROFILE("PsatdAlgorithmJConstantInTime::CurrentCorrection"); const SpectralFieldIndex& Idx = m_spectral_index; @@ -833,10 +833,10 @@ void PsatdAlgorithm::CurrentCorrection (SpectralFieldData& field_data) } void -PsatdAlgorithm::VayDeposition (SpectralFieldData& field_data) +PsatdAlgorithmJConstantInTime::VayDeposition (SpectralFieldData& field_data) { // Profiling - BL_PROFILE("PsatdAlgorithm::VayDeposition()"); + BL_PROFILE("PsatdAlgorithmJConstantInTime::VayDeposition()"); const SpectralFieldIndex& Idx = m_spectral_index; diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H index 608da5fd500..097a1a9d69b 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H @@ -23,7 +23,8 @@ class PsatdAlgorithmRZ : public SpectralBaseAlgorithmRZ bool const nodal, amrex::Real const dt_step, bool const update_with_rho, const bool time_averaging, - const bool do_multi_J, + const int J_in_time, + const int rho_in_time, const bool dive_cleaning, const bool divb_cleaning); // Redefine functions from base class @@ -62,7 +63,7 @@ class PsatdAlgorithmRZ : public SpectralBaseAlgorithmRZ amrex::Real m_dt; bool m_update_with_rho; bool m_time_averaging; - bool m_do_multi_J; + int m_J_in_time; bool m_dive_cleaning; bool m_divb_cleaning; SpectralRealCoefficients C_coef, S_ck_coef, X1_coef, X2_coef, X3_coef; diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp index fec0ebc8f62..55b58821ce9 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp @@ -6,6 +6,7 @@ */ #include "PsatdAlgorithmRZ.H" #include "Utils/TextMsg.H" +#include "Utils/WarpXAlgorithmSelection.H" #include "Utils/WarpXConst.H" #include "Utils/WarpXProfilerWrapper.H" #include "WarpX.H" @@ -23,7 +24,8 @@ PsatdAlgorithmRZ::PsatdAlgorithmRZ (SpectralKSpaceRZ const & spectral_kspace, bool const nodal, amrex::Real const dt, bool const update_with_rho, const bool time_averaging, - const bool do_multi_J, + const int J_in_time, + const int rho_in_time, const bool dive_cleaning, const bool divb_cleaning) // Initialize members of base class @@ -32,10 +34,11 @@ PsatdAlgorithmRZ::PsatdAlgorithmRZ (SpectralKSpaceRZ const & spectral_kspace, m_dt(dt), m_update_with_rho(update_with_rho), m_time_averaging(time_averaging), - m_do_multi_J(do_multi_J), + m_J_in_time(J_in_time), m_dive_cleaning(dive_cleaning), m_divb_cleaning(divb_cleaning) { + amrex::ignore_unused(rho_in_time); // Allocate the arrays of coefficients amrex::BoxArray const & ba = spectral_kspace.spectralspace_ba; @@ -47,28 +50,28 @@ PsatdAlgorithmRZ::PsatdAlgorithmRZ (SpectralKSpaceRZ const & spectral_kspace, coefficients_initialized = false; - if (time_averaging && do_multi_J) + if (time_averaging && J_in_time == JInTime::Linear) { X5_coef = SpectralRealCoefficients(ba, dm, n_rz_azimuthal_modes, 0); X6_coef = SpectralRealCoefficients(ba, dm, n_rz_azimuthal_modes, 0); } - if (time_averaging && !do_multi_J) + if (time_averaging && J_in_time != JInTime::Linear) { amrex::Abort(Utils::TextMsg::Err( - "RZ PSATD: psatd.do_time_averaging = 1 implemented only with warpx.do_multi_J = 1")); + "RZ PSATD: psatd.do_time_averaging=1 implemented only with psatd.J_in_time=linear")); } - if (dive_cleaning && !do_multi_J) + if (dive_cleaning && J_in_time != JInTime::Linear) { amrex::Abort(Utils::TextMsg::Err( - "RZ PSATD: warpx.do_dive_cleaning = 1 implemented only with warpx.do_multi_J = 1")); + "RZ PSATD: warpx.do_dive_cleaning=1 implemented only with psatd.J_in_time=linear")); } - if (divb_cleaning && !do_multi_J) + if (divb_cleaning && J_in_time != JInTime::Linear) { amrex::Abort(Utils::TextMsg::Err( - "RZ PSATD: warpx.do_divb_cleaning = 1 implemented only with warpx.do_multi_J = 1")); + "RZ PSATD: warpx.do_divb_cleaning=1 implemented only with psatd.J_in_time=linear")); } } @@ -80,7 +83,7 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f) const bool update_with_rho = m_update_with_rho; const bool time_averaging = m_time_averaging; - const bool do_multi_J = m_do_multi_J; + const bool J_in_time_linear = (m_J_in_time == JInTime::Linear) ? true : false; const bool dive_cleaning = m_dive_cleaning; const bool divb_cleaning = m_divb_cleaning; @@ -109,7 +112,7 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f) amrex::Array4 X5_arr; amrex::Array4 X6_arr; - if (time_averaging && do_multi_J) + if (time_averaging && J_in_time_linear) { X5_arr = X5_coef[mfi].array(); X6_arr = X6_coef[mfi].array(); @@ -235,7 +238,7 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f) G_old = fields(i,j,k,G_m); } - if (do_multi_J) + if (J_in_time_linear) { const int Jp_m_new = Idx.Jx_new + Idx.n_fields*mode; const int Jm_m_new = Idx.Jy_new + Idx.n_fields*mode; @@ -332,7 +335,7 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f) void PsatdAlgorithmRZ::InitializeSpectralCoefficients (SpectralFieldDataRZ const & f) { const bool time_averaging = m_time_averaging; - const bool do_multi_J = m_do_multi_J; + const bool J_in_time_linear = (m_J_in_time == JInTime::Linear) ? true : false; // Fill them with the right values: // Loop over boxes and allocate the corresponding coefficients @@ -353,7 +356,7 @@ void PsatdAlgorithmRZ::InitializeSpectralCoefficients (SpectralFieldDataRZ const amrex::Array4 X5; amrex::Array4 X6; - if (time_averaging && do_multi_J) + if (time_averaging && J_in_time_linear) { X5 = X5_coef[mfi].array(); X6 = X6_coef[mfi].array(); @@ -392,7 +395,7 @@ void PsatdAlgorithmRZ::InitializeSpectralCoefficients (SpectralFieldDataRZ const X3(i,j,k,mode) = - c*c * dt*dt / (3._rt*ep0); } - if (time_averaging && do_multi_J) + if (time_averaging && J_in_time_linear) { constexpr amrex::Real c2 = PhysConst::c; const amrex::Real dt3 = dt * dt * dt; diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H index 811fdd4d73c..4ab88f1a378 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -55,7 +55,8 @@ class SpectralFieldIndex */ SpectralFieldIndex (const bool update_with_rho, const bool time_averaging, - const bool do_multi_J, + const int J_in_time, + const int rho_in_time, const bool dive_cleaning, const bool divb_cleaning, const bool pml, @@ -89,8 +90,13 @@ class SpectralFieldIndex int Ex_avg = -1, Ey_avg = -1, Ez_avg = -1; int Bx_avg = -1, By_avg = -1, Bz_avg = -1; - // Multi-J, div(E) and div(B) cleaning + // J linear in time int Jx_new = -1, Jy_new = -1, Jz_new = -1; + + // rho quadratic in time + int rho_mid = -1; + + // div(E) and div(B) cleaning int F = -1, G = -1; // PML diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 00346a2c772..56189e0ff06 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -35,7 +35,8 @@ using namespace amrex; SpectralFieldIndex::SpectralFieldIndex (const bool update_with_rho, const bool time_averaging, - const bool do_multi_J, + const int J_in_time, + const int rho_in_time, const bool dive_cleaning, const bool divb_cleaning, const bool pml, @@ -68,13 +69,18 @@ SpectralFieldIndex::SpectralFieldIndex (const bool update_with_rho, if (divb_cleaning) G = c++; - if (do_multi_J) + if (J_in_time == JInTime::Linear) { Jx_new = c++; Jy_new = c++; Jz_new = c++; } + if (rho_in_time == RhoInTime::Quadratic) + { + rho_mid = c++; + } + if (pml_rz) { Er_pml = c++; Et_pml = c++; diff --git a/Source/FieldSolver/SpectralSolver/SpectralHankelTransform/HankelTransform.H b/Source/FieldSolver/SpectralSolver/SpectralHankelTransform/HankelTransform.H index 5a87031c40c..3f74f5d82ab 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralHankelTransform/HankelTransform.H +++ b/Source/FieldSolver/SpectralSolver/SpectralHankelTransform/HankelTransform.H @@ -8,6 +8,14 @@ #define WARPX_HANKEL_TRANSFORM_H_ #include +#include +#include + +#ifdef AMREX_USE_GPU +# include +#endif + +#include /* \brief This defines the class that performs the Hankel transform. * Original authors: Remi Lehe, Manuel Kirchen @@ -45,6 +53,10 @@ class HankelTransform RealVector m_invM; RealVector m_M; + +#ifdef AMREX_USE_GPU + std::unique_ptr m_queue; +#endif }; #endif diff --git a/Source/FieldSolver/SpectralSolver/SpectralHankelTransform/HankelTransform.cpp b/Source/FieldSolver/SpectralSolver/SpectralHankelTransform/HankelTransform.cpp index 43b26f2ee33..ddd07acad29 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralHankelTransform/HankelTransform.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralHankelTransform/HankelTransform.cpp @@ -11,6 +11,8 @@ #include "Utils/WarpXConst.H" #include "WarpX.H" +#include "Utils/WarpXProfilerWrapper.H" + #include #include @@ -23,10 +25,20 @@ HankelTransform::HankelTransform (int const hankel_order, : m_nr(nr), m_nk(nr) { + WARPX_PROFILE("HankelTransform::HankelTransform"); + // Check that azimuthal_mode has a valid value WARPX_ALWAYS_ASSERT_WITH_MESSAGE(hankel_order-1 <= azimuthal_mode && azimuthal_mode <= hankel_order+1, "azimuthal_mode must be either hankel_order-1, hankel_order or hankel_order+1"); +#ifdef AMREX_USE_GPU + // BLAS setup + // SYCL note: we need to double check AMReX device ID conventions and + // BLAS++ device ID conventions are the same + int const device_id = amrex::Gpu::Device::deviceId(); + m_queue = std::make_unique( device_id, 0 ); +#endif + amrex::Vector alphas; amrex::Vector alpha_errors; @@ -186,6 +198,8 @@ void HankelTransform::HankelForwardTransform (amrex::FArrayBox const& F, int const F_icomp, amrex::FArrayBox & G, int const G_icomp) { + WARPX_PROFILE("HankelTransform::HankelForwardTransform"); + amrex::Box const& F_box = F.box(); amrex::Box const& G_box = G.box(); @@ -198,37 +212,24 @@ HankelTransform::HankelForwardTransform (amrex::FArrayBox const& F, int const F_ AMREX_ALWAYS_ASSERT(ngr >= 0); AMREX_ALWAYS_ASSERT(F_box.bigEnd(0)+1 >= m_nr); -#ifndef AMREX_USE_GPU - // On CPU, the blas::gemm is significantly faster + // We perform stream synchronization since `gemm` may be running + // on a different stream. + amrex::Gpu::streamSynchronize(); // Note that M is flagged to be transposed since it has dimensions (m_nr, m_nk) blas::gemm(blas::Layout::ColMajor, blas::Op::Trans, blas::Op::NoTrans, m_nk, nz, m_nr, 1._rt, m_M.dataPtr(), m_nk, F.dataPtr(F_icomp)+ngr, nrF, 0._rt, - G.dataPtr(G_icomp), m_nk); - -#else - // On GPU, the explicit loop is significantly faster - // It is not clear if the GPU gemm wasn't build properly, it is cycling data out and back - // in to the device, or if it is because gemm is launching its own threads. - - amrex::Real const * M_arr = m_M.dataPtr(); - amrex::Array4 const & F_arr = F.array(); - amrex::Array4< amrex::Real> const & G_arr = G.array(); - - int const nr = m_nr; - - amrex::ParallelFor(G_box, - [=] AMREX_GPU_DEVICE(int ik, int iz, int k3d) noexcept { - G_arr(ik,iz,k3d,G_icomp) = 0.; - for (int ir=0 ; ir < nr ; ir++) { - int const ii = ir + ik*nr; - G_arr(ik,iz,k3d,G_icomp) += M_arr[ii]*F_arr(ir,iz,k3d,F_icomp); - } - }); - + G.dataPtr(G_icomp), m_nk +#ifdef AMREX_USE_GPU + , *m_queue // Calls the GPU version of blas::gemm #endif + ); + + // We perform stream synchronization since `gemm` may be running + // on a different stream. + amrex::Gpu::streamSynchronize(); } @@ -236,6 +237,8 @@ void HankelTransform::HankelInverseTransform (amrex::FArrayBox const& G, int const G_icomp, amrex::FArrayBox & F, int const F_icomp) { + WARPX_PROFILE("HankelTransform::HankelInverseTransform"); + amrex::Box const& G_box = G.box(); amrex::Box const& F_box = F.box(); @@ -248,36 +251,22 @@ HankelTransform::HankelInverseTransform (amrex::FArrayBox const& G, int const G_ AMREX_ALWAYS_ASSERT(ngr >= 0); AMREX_ALWAYS_ASSERT(F_box.bigEnd(0)+1 >= m_nr); -#ifndef AMREX_USE_GPU - // On CPU, the blas::gemm is significantly faster + // We perform stream synchronization since `gemm` may be running + // on a different stream. + amrex::Gpu::streamSynchronize(); // Note that m_invM is flagged to be transposed since it has dimensions (m_nk, m_nr) blas::gemm(blas::Layout::ColMajor, blas::Op::Trans, blas::Op::NoTrans, m_nr, nz, m_nk, 1._rt, m_invM.dataPtr(), m_nr, G.dataPtr(G_icomp), m_nk, 0._rt, - F.dataPtr(F_icomp)+ngr, nrF); - -#else - // On GPU, the explicit loop is significantly faster - // It is not clear if the GPU gemm wasn't build properly, it is cycling data out and back - // in to the device, or if it is because gemm is launching its own threads. - - amrex::Real const * invM_arr = m_invM.dataPtr(); - amrex::Array4 const & G_arr = G.array(); - amrex::Array4< amrex::Real> const & F_arr = F.array(); - - int const nk = m_nk; - - amrex::ParallelFor(G_box, - [=] AMREX_GPU_DEVICE(int ir, int iz, int k3d) noexcept { - F_arr(ir,iz,k3d,F_icomp) = 0.; - for (int ik=0 ; ik < nk ; ik++) { - int const ii = ik + ir*nk; - F_arr(ir,iz,k3d,F_icomp) += invM_arr[ii]*G_arr(ik,iz,k3d,G_icomp); - } - }); - + F.dataPtr(F_icomp)+ngr, nrF +#ifdef AMREX_USE_GPU + , *m_queue // Calls the GPU version of blas::gemm #endif + ); + // We perform stream synchronization since `gemm` may be running + // on a different stream. + amrex::Gpu::streamSynchronize(); } diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H index 684cf9586b8..38b2420105a 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H @@ -58,9 +58,10 @@ class SpectralSolver * (no domain decomposition) * \param[in] update_with_rho whether rho is used in the field update equations * \param[in] fft_do_time_averaging whether the time averaging algorithm is used - * \param[in] do_multi_J whether the multi-J algorithm is used (hence two currents - * computed at the beginning and the end of the time interval - * instead of one current computed at half time) + * \param[in] J_in_time integer that corresponds to the time dependency of J + * (constant, linear) for the PSATD algorithm + * \param[in] rho_in_time integer that corresponds to the time dependency of rho + * (linear, quadratic) for the PSATD algorithm * \param[in] dive_cleaning whether to use div(E) cleaning to account for errors in * Gauss law (new field F in the update equations) * \param[in] divb_cleaning whether to use div(B) cleaning to account for errors in @@ -79,7 +80,8 @@ class SpectralSolver const bool periodic_single_box, const bool update_with_rho, const bool fft_do_time_averaging, - const bool do_multi_J, + const int J_in_time, + const int rho_in_time, const bool dive_cleaning, const bool divb_cleaning); diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp index 75c82319c11..c0d7b89412b 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp @@ -8,7 +8,7 @@ #include "FieldSolver/SpectralSolver/SpectralFieldData.H" #include "SpectralAlgorithms/PsatdAlgorithmComoving.H" #include "SpectralAlgorithms/PsatdAlgorithmPml.H" -#include "SpectralAlgorithms/PsatdAlgorithm.H" +#include "SpectralAlgorithms/PsatdAlgorithmJConstantInTime.H" #include "SpectralAlgorithms/PsatdAlgorithmJLinearInTime.H" #include "SpectralKSpace.H" #include "SpectralSolver.H" @@ -30,7 +30,8 @@ SpectralSolver::SpectralSolver( const bool pml, const bool periodic_single_box, const bool update_with_rho, const bool fft_do_time_averaging, - const bool do_multi_J, + const int J_in_time, + const int rho_in_time, const bool dive_cleaning, const bool divb_cleaning) { @@ -41,8 +42,9 @@ SpectralSolver::SpectralSolver( // as well as the value of the corresponding k coordinates) const SpectralKSpace k_space= SpectralKSpace(realspace_ba, dm, dx); - m_spectral_index = SpectralFieldIndex(update_with_rho, fft_do_time_averaging, - do_multi_J, dive_cleaning, divb_cleaning, pml); + m_spectral_index = SpectralFieldIndex( + update_with_rho, fft_do_time_averaging, J_in_time, rho_in_time, + dive_cleaning, divb_cleaning, pml); // - Select the algorithm depending on the input parameters // Initialize the corresponding coefficients over k space @@ -64,18 +66,18 @@ SpectralSolver::SpectralSolver( } else // PSATD algorithms: standard, Galilean, averaged Galilean, multi-J { - if (do_multi_J) + if (J_in_time == JInTime::Constant) { - algorithm = std::make_unique( + algorithm = std::make_unique( k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal, - dt, fft_do_time_averaging, dive_cleaning, divb_cleaning); + v_galilean, dt, update_with_rho, fft_do_time_averaging, + dive_cleaning, divb_cleaning); } - else // standard, Galilean, averaged Galilean + else // J linear in time { - algorithm = std::make_unique( + algorithm = std::make_unique( k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal, - v_galilean, dt, update_with_rho, fft_do_time_averaging, - dive_cleaning, divb_cleaning); + dt, fft_do_time_averaging, dive_cleaning, divb_cleaning); } } } diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H index 30c92251253..45b55d9d4da 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H @@ -37,7 +37,8 @@ class SpectralSolverRZ bool const with_pml, bool const update_with_rho, const bool fft_do_time_averaging, - const bool do_multi_J, + const int J_in_time, + const int rho_in_time, const bool dive_cleaning, const bool divb_cleaning); diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp index 3f4bddd8b16..23d93fe045b 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp @@ -35,7 +35,8 @@ SpectralSolverRZ::SpectralSolverRZ (const int lev, bool const with_pml, bool const update_with_rho, const bool fft_do_time_averaging, - const bool do_multi_J, + const int J_in_time, + const int rho_in_time, const bool dive_cleaning, const bool divb_cleaning) : k_space(realspace_ba, dm, dx) @@ -47,8 +48,9 @@ SpectralSolverRZ::SpectralSolverRZ (const int lev, // as well as the value of the corresponding k coordinates. const bool is_pml = false; - m_spectral_index = SpectralFieldIndex(update_with_rho, fft_do_time_averaging, - do_multi_J, dive_cleaning, divb_cleaning, is_pml, with_pml); + m_spectral_index = SpectralFieldIndex( + update_with_rho, fft_do_time_averaging, J_in_time, rho_in_time, + dive_cleaning, divb_cleaning, is_pml, with_pml); // - Select the algorithm depending on the input parameters // Initialize the corresponding coefficients over k space @@ -60,7 +62,7 @@ SpectralSolverRZ::SpectralSolverRZ (const int lev, // v_galilean is 0: use standard PSATD algorithm algorithm = std::make_unique( k_space, dm, m_spectral_index, n_rz_azimuthal_modes, norder_z, nodal, dt, - update_with_rho, fft_do_time_averaging, do_multi_J, dive_cleaning, divb_cleaning); + update_with_rho, fft_do_time_averaging, J_in_time, rho_in_time, dive_cleaning, divb_cleaning); } else { // Otherwise: use the Galilean algorithm algorithm = std::make_unique( diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 7499ee140a2..e8f20661a2d 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -280,9 +280,9 @@ void WarpX::PSATDForwardTransformJ ( { Idx = spectral_solver_fp[lev]->m_spectral_index; - idx_jx = (WarpX::do_multi_J) ? static_cast(Idx.Jx_new) : static_cast(Idx.Jx); - idx_jy = (WarpX::do_multi_J) ? static_cast(Idx.Jy_new) : static_cast(Idx.Jy); - idx_jz = (WarpX::do_multi_J) ? static_cast(Idx.Jz_new) : static_cast(Idx.Jz); + idx_jx = (J_in_time == JInTime::Linear) ? static_cast(Idx.Jx_new) : static_cast(Idx.Jx); + idx_jy = (J_in_time == JInTime::Linear) ? static_cast(Idx.Jy_new) : static_cast(Idx.Jy); + idx_jz = (J_in_time == JInTime::Linear) ? static_cast(Idx.Jz_new) : static_cast(Idx.Jz); ForwardTransformVect(lev, *spectral_solver_fp[lev], J_fp[lev], idx_jx, idx_jy, idx_jz); @@ -290,9 +290,9 @@ void WarpX::PSATDForwardTransformJ ( { Idx = spectral_solver_cp[lev]->m_spectral_index; - idx_jx = (WarpX::do_multi_J) ? static_cast(Idx.Jx_new) : static_cast(Idx.Jx); - idx_jy = (WarpX::do_multi_J) ? static_cast(Idx.Jy_new) : static_cast(Idx.Jy); - idx_jz = (WarpX::do_multi_J) ? static_cast(Idx.Jz_new) : static_cast(Idx.Jz); + idx_jx = (J_in_time == JInTime::Linear) ? static_cast(Idx.Jx_new) : static_cast(Idx.Jx); + idx_jy = (J_in_time == JInTime::Linear) ? static_cast(Idx.Jy_new) : static_cast(Idx.Jy); + idx_jz = (J_in_time == JInTime::Linear) ? static_cast(Idx.Jz_new) : static_cast(Idx.Jz); ForwardTransformVect(lev, *spectral_solver_cp[lev], J_cp[lev], idx_jx, idx_jy, idx_jz); } diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp index 3eced8f5d96..07526340f2e 100644 --- a/Source/Initialization/PlasmaInjector.cpp +++ b/Source/Initialization/PlasmaInjector.cpp @@ -156,7 +156,8 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name) charge_is_specified || species_is_specified || (injection_style == "external_file"), - "Need to specify at least one of species_type or charge" + "Need to specify at least one of species_type or charge for species '" + + species_name + "'." ); if ( mass_is_specified && species_is_specified ){ @@ -170,7 +171,8 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name) mass_is_specified || species_is_specified || (injection_style == "external_file"), - "Need to specify at least one of species_type or mass" + "Need to specify at least one of species_type or mass for species '" + + species_name + "'." ); num_particles_per_cell_each_dim.assign(3, 0); @@ -218,7 +220,7 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name) queryWithParser(pp_species_name, "y_cut", y_cut); queryWithParser(pp_species_name, "z_cut", z_cut); getWithParser(pp_species_name, "q_tot", q_tot); - pp_species_name.get("npart", npart); + getWithParser(pp_species_name, "npart", npart); pp_species_name.query("do_symmetrize", do_symmetrize); gaussian_beam = true; parseMomentum(pp_species_name); @@ -240,7 +242,7 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name) // so that inj_pos->getPositionUnitBox calls // InjectorPosition[Random or Regular].getPositionUnitBox. else if (injection_style == "nrandompercell") { - queryWithParser(pp_species_name, "num_particles_per_cell", num_particles_per_cell); + getWithParser(pp_species_name, "num_particles_per_cell", num_particles_per_cell); #if WARPX_DIM_RZ if (WarpX::n_rz_azimuthal_modes > 1) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE( @@ -258,7 +260,7 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name) parseMomentum(pp_species_name); } else if (injection_style == "nfluxpercell") { surface_flux = true; - queryWithParser(pp_species_name, "num_particles_per_cell", num_particles_per_cell_real); + getWithParser(pp_species_name, "num_particles_per_cell", num_particles_per_cell_real); #ifdef WARPX_DIM_RZ if (WarpX::n_rz_azimuthal_modes > 1) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE( diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 24017d32190..46a98b31321 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -21,6 +21,7 @@ #include "Filter/BilinearFilter.H" #include "Filter/NCIGodfreyFilter.H" #include "Particles/MultiParticleContainer.H" +#include "Utils/Logo/GetLogo.H" #include "Utils/MPIInitHelpers.H" #include "Utils/TextMsg.H" #include "Utils/WarpXAlgorithmSelection.H" @@ -30,6 +31,7 @@ #include "Initialization/ReconnectionPerturbation.H" #include +#include #include #include @@ -65,7 +67,6 @@ #include #include #include -#include #include #include #include @@ -354,14 +355,7 @@ WarpX::PrintMainPICparameters () void WarpX::WriteUsedInputsFile (std::string const & filename) const { - amrex::Print() << "For full input parameters, see the file: " << filename << "\n\n"; - - if (ParallelDescriptor::IOProcessor()) { - std::ofstream jobInfoFile; - jobInfoFile.open(filename.c_str(), std::ios::out); - ParmParse::dumpTable(jobInfoFile, true); - jobInfoFile.close(); - } + ablastr::utils::write_used_inputs_file(filename); } void @@ -370,11 +364,14 @@ WarpX::InitData () WARPX_PROFILE("WarpX::InitData()"); utils::warpx_check_mpi_thread_level(); - Print() << "WarpX (" << WarpX::Version() << ")\n"; #ifdef WARPX_QED Print() << "PICSAR (" << WarpX::PicsarVersion() << ")\n"; #endif + Print() << "WarpX (" << WarpX::Version() << ")\n"; + + Print() << utils::logo::get_logo(); + if (restart_chkfile.empty()) { ComputeDt(); @@ -396,10 +393,6 @@ WarpX::InitData () WarpX::InitNCICorrector(); } - if (WarpX::use_filter) { - WarpX::InitFilter(); - } - BuildBufferMasks(); if (WarpX::em_solver_medium==1) { @@ -508,7 +501,7 @@ WarpX::InitPML () pml_ncell, pml_delta, amrex::IntVect::TheZeroVector(), dt[0], nox_fft, noy_fft, noz_fft, do_nodal, do_moving_window, pml_has_particles, do_pml_in_domain, - do_multi_J, + J_in_time, rho_in_time, do_pml_dive_cleaning, do_pml_divb_cleaning, amrex::IntVect(0), amrex::IntVect(0), guard_cells.ng_FieldSolver.max(), @@ -545,7 +538,7 @@ WarpX::InitPML () pml_ncell, pml_delta, refRatio(lev-1), dt[lev], nox_fft, noy_fft, noz_fft, do_nodal, do_moving_window, pml_has_particles, do_pml_in_domain, - do_multi_J, do_pml_dive_cleaning, do_pml_divb_cleaning, + J_in_time, rho_in_time, do_pml_dive_cleaning, do_pml_divb_cleaning, amrex::IntVect(0), amrex::IntVect(0), guard_cells.ng_FieldSolver.max(), v_particle_pml, diff --git a/Source/Parallelization/GuardCellManager.H b/Source/Parallelization/GuardCellManager.H index 05d99de4df6..c122aaa96c9 100644 --- a/Source/Parallelization/GuardCellManager.H +++ b/Source/Parallelization/GuardCellManager.H @@ -71,7 +71,9 @@ public: const bool do_pml, const int do_pml_in_domain, const int pml_ncell, - const amrex::Vector& ref_ratios); + const amrex::Vector& ref_ratios, + const bool use_filter, + const amrex::IntVect& bilinear_filter_stencil_length); // Guard cells allocated for MultiFabs E and B amrex::IntVect ng_alloc_EB = amrex::IntVect::TheZeroVector(); diff --git a/Source/Parallelization/GuardCellManager.cpp b/Source/Parallelization/GuardCellManager.cpp index cfe2c0e4dd2..d32d696190e 100644 --- a/Source/Parallelization/GuardCellManager.cpp +++ b/Source/Parallelization/GuardCellManager.cpp @@ -53,7 +53,9 @@ guardCellManager::Init ( const bool do_pml, const int do_pml_in_domain, const int pml_ncell, - const amrex::Vector& ref_ratios) + const amrex::Vector& ref_ratios, + const bool use_filter, + const amrex::IntVect& bilinear_filter_stencil_length) { // When using subcycling, the particles on the fine level perform two pushes // before being redistributed ; therefore, we need one extra guard cell @@ -160,6 +162,11 @@ guardCellManager::Init ( ng_depos_J = ng_alloc_J; ng_depos_rho = ng_alloc_Rho; + if (use_filter) + { + ng_alloc_J += bilinear_filter_stencil_length - amrex::IntVect(1); + } + // After pushing particle int ng_alloc_F_int = (do_moving_window) ? 2 : 0; // CKC solver requires one additional guard cell diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 27c01b6b3f1..373cb281aa5 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -953,43 +953,79 @@ void WarpX::RestrictCurrentFromFineToCoarsePatch ( CoarsenMR::Coarsen( *crse[2], *fine[2], refinement_ratio ); } -void WarpX::ApplyFilterandSumBoundaryJ ( - const amrex::Vector,3>>& J_fp, - const amrex::Vector,3>>& J_cp, +void WarpX::ApplyFilterJ ( + const amrex::Vector,3>>& current, const int lev, - PatchType patch_type) + const int idim) { - const int glev = (patch_type == PatchType::fine) ? lev : lev-1; - const amrex::Periodicity& period = Geom(glev).periodicity(); - const std::array,3>& j = (patch_type == PatchType::fine) ? - J_fp[lev] : J_cp[lev]; - for (int idim = 0; idim < 3; ++idim) { - IntVect ng = j[idim]->nGrowVect(); - IntVect ng_depos_J = get_ng_depos_J(); - if (WarpX::do_current_centering) - { + amrex::MultiFab& J = *current[lev][idim]; + + const int ncomp = J.nComp(); + const amrex::IntVect ngrow = J.nGrowVect(); + amrex::MultiFab Jf(J.boxArray(), J.DistributionMap(), ncomp, ngrow); + bilinear_filter.ApplyStencil(Jf, J, lev); + + const int srccomp = 0; + const int dstcomp = 0; + amrex::MultiFab::Copy(J, Jf, srccomp, dstcomp, ncomp, ngrow); +} + +void WarpX::ApplyFilterJ ( + const amrex::Vector,3>>& current, + const int lev) +{ + for (int idim=0; idim<3; ++idim) + { + ApplyFilterJ(current, lev, idim); + } +} + +void WarpX::SumBoundaryJ ( + const amrex::Vector,3>>& current, + const int lev, + const int idim, + const amrex::Periodicity& period) +{ + amrex::MultiFab& J = *current[lev][idim]; + + amrex::IntVect ng = J.nGrowVect(); + amrex::IntVect ng_depos_J = get_ng_depos_J(); + + if (WarpX::do_current_centering) + { #if defined(WARPX_DIM_1D_Z) - ng_depos_J[0] += WarpX::current_centering_noz / 2; + ng_depos_J[0] += WarpX::current_centering_noz / 2; #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - ng_depos_J[0] += WarpX::current_centering_nox / 2; - ng_depos_J[1] += WarpX::current_centering_noz / 2; + ng_depos_J[0] += WarpX::current_centering_nox / 2; + ng_depos_J[1] += WarpX::current_centering_noz / 2; #elif defined(WARPX_DIM_3D) - ng_depos_J[0] += WarpX::current_centering_nox / 2; - ng_depos_J[1] += WarpX::current_centering_noy / 2; - ng_depos_J[2] += WarpX::current_centering_noz / 2; + ng_depos_J[0] += WarpX::current_centering_nox / 2; + ng_depos_J[1] += WarpX::current_centering_noy / 2; + ng_depos_J[2] += WarpX::current_centering_noz / 2; #endif - } - if (use_filter) { - ng += bilinear_filter.stencil_length_each_dir-1; - ng_depos_J += bilinear_filter.stencil_length_each_dir-1; - ng_depos_J.min(ng); - MultiFab jf(j[idim]->boxArray(), j[idim]->DistributionMap(), j[idim]->nComp(), ng); - bilinear_filter.ApplyStencil(jf, *j[idim], lev); - WarpXSumGuardCells(*(j[idim]), jf, period, ng_depos_J, 0, (j[idim])->nComp()); - } else { - ng_depos_J.min(ng); - WarpXSumGuardCells(*(j[idim]), period, ng_depos_J, 0, (j[idim])->nComp()); - } + } + + if (use_filter) + { + ng_depos_J += bilinear_filter.stencil_length_each_dir - amrex::IntVect(1); + } + + ng_depos_J.min(ng); + + const amrex::IntVect src_ngrow = ng_depos_J; + const int icomp = 0; + const int ncomp = J.nComp(); + WarpXSumGuardCells(J, period, src_ngrow, icomp, ncomp); +} + +void WarpX::SumBoundaryJ ( + const amrex::Vector,3>>& current, + const int lev, + const amrex::Periodicity& period) +{ + for (int idim=0; idim<3; ++idim) + { + SumBoundaryJ(current, lev, idim, period); } } @@ -1011,88 +1047,73 @@ void WarpX::AddCurrentFromFineLevelandSumBoundary ( const amrex::Vector,3>>& J_cp, const int lev) { - ApplyFilterandSumBoundaryJ(J_fp, J_cp, lev, PatchType::fine); + const amrex::Periodicity& period = Geom(lev).periodicity(); - if (lev < finest_level) { + if (use_filter) + { + ApplyFilterJ(J_fp, lev); + } + SumBoundaryJ(J_fp, lev, period); + + if (lev < finest_level) + { // When there are current buffers, unlike coarse patch, // we don't care about the final state of them. - const amrex::Periodicity& period = Geom(lev).periodicity(); - for (int idim = 0; idim < 3; ++idim) { + for (int idim=0; idim<3; ++idim) + { MultiFab mf(J_fp[lev][idim]->boxArray(), J_fp[lev][idim]->DistributionMap(), J_fp[lev][idim]->nComp(), 0); mf.setVal(0.0); + IntVect ng = J_cp[lev+1][idim]->nGrowVect(); - IntVect ng_depos_J = get_ng_depos_J(); - if (WarpX::do_current_centering) - { -#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - ng_depos_J[0] += WarpX::current_centering_nox / 2; - ng_depos_J[1] += WarpX::current_centering_noz / 2; -#elif defined(WARPX_DIM_3D) - ng_depos_J[0] += WarpX::current_centering_nox / 2; - ng_depos_J[1] += WarpX::current_centering_noy / 2; - ng_depos_J[2] += WarpX::current_centering_noz / 2; -#endif - } + if (use_filter && current_buf[lev+1][idim]) { - // coarse patch of fine level - ng += bilinear_filter.stencil_length_each_dir-1; - ng_depos_J += bilinear_filter.stencil_length_each_dir-1; - ng_depos_J.min(ng); - MultiFab jfc(J_cp[lev+1][idim]->boxArray(), - J_cp[lev+1][idim]->DistributionMap(), J_cp[lev+1][idim]->nComp(), ng); - bilinear_filter.ApplyStencil(jfc, *J_cp[lev+1][idim], lev+1); - - // buffer patch of fine level - MultiFab jfb(current_buf[lev+1][idim]->boxArray(), - current_buf[lev+1][idim]->DistributionMap(), current_buf[lev+1][idim]->nComp(), ng); - bilinear_filter.ApplyStencil(jfb, *current_buf[lev+1][idim], lev+1); - - MultiFab::Add(jfb, jfc, 0, 0, current_buf[lev+1][idim]->nComp(), ng); - ablastr::utils::communication::ParallelAdd(mf, jfb, 0, 0, current_buf[lev + 1][idim]->nComp(), - ng, IntVect::TheZeroVector(), WarpX::do_single_precision_comms, period); - - WarpXSumGuardCells(*J_cp[lev+1][idim], jfc, period, ng_depos_J, 0, J_cp[lev+1][idim]->nComp()); + ApplyFilterJ(J_cp, lev+1, idim); + ApplyFilterJ(current_buf, lev+1, idim); + + MultiFab::Add( + *current_buf[lev+1][idim], *J_cp[lev+1][idim], + 0, 0, current_buf[lev+1][idim]->nComp(), ng); + + ablastr::utils::communication::ParallelAdd( + mf, *current_buf[lev+1][idim], 0, 0, + current_buf[lev+1][idim]->nComp(), + ng, amrex::IntVect(0), + do_single_precision_comms, period); } else if (use_filter) // but no buffer { - // coarse patch of fine level - ng += bilinear_filter.stencil_length_each_dir-1; - ng_depos_J += bilinear_filter.stencil_length_each_dir-1; - ng_depos_J.min(ng); - MultiFab jf(J_cp[lev+1][idim]->boxArray(), - J_cp[lev+1][idim]->DistributionMap(), J_cp[lev+1][idim]->nComp(), ng); - bilinear_filter.ApplyStencil(jf, *J_cp[lev+1][idim], lev+1); - - ablastr::utils::communication::ParallelAdd(mf, jf, 0, 0, J_cp[lev + 1][idim]->nComp(), ng, - IntVect::TheZeroVector(), WarpX::do_single_precision_comms, period); - WarpXSumGuardCells(*J_cp[lev+1][idim], jf, period, ng_depos_J, 0, J_cp[lev+1][idim]->nComp()); + ApplyFilterJ(J_cp, lev+1, idim); + + ablastr::utils::communication::ParallelAdd( + mf, *J_cp[lev+1][idim], 0, 0, + J_cp[lev+1][idim]->nComp(), + ng, amrex::IntVect(0), + do_single_precision_comms, period); } else if (current_buf[lev+1][idim]) // but no filter { - ng_depos_J.min(ng); - MultiFab::Add(*current_buf[lev+1][idim], - *J_cp [lev+1][idim], 0, 0, current_buf[lev+1][idim]->nComp(), - J_cp[lev+1][idim]->nGrowVect()); - ablastr::utils::communication::ParallelAdd(mf, *current_buf[lev + 1][idim], 0, 0, - current_buf[lev + 1][idim]->nComp(), - current_buf[lev + 1][idim]->nGrowVect(), - IntVect::TheZeroVector(), WarpX::do_single_precision_comms, - period); - WarpXSumGuardCells(*(J_cp[lev+1][idim]), period, ng_depos_J, 0, J_cp[lev+1][idim]->nComp()); + MultiFab::Add( + *current_buf[lev+1][idim], *J_cp[lev+1][idim], + 0, 0, current_buf[lev+1][idim]->nComp(), ng); + + ablastr::utils::communication::ParallelAdd( + mf, *current_buf[lev+1][idim], 0, 0, + current_buf[lev+1][idim]->nComp(), + ng, amrex::IntVect(0), + do_single_precision_comms, period); } else // no filter, no buffer { - ng_depos_J.min(ng); - ablastr::utils::communication::ParallelAdd(mf, *J_cp[lev + 1][idim], 0, 0, - J_cp[lev + 1][idim]->nComp(), - J_cp[lev + 1][idim]->nGrowVect(), - IntVect::TheZeroVector(), WarpX::do_single_precision_comms, - period); - WarpXSumGuardCells(*(J_cp[lev+1][idim]), period, ng_depos_J, 0, J_cp[lev+1][idim]->nComp()); + ablastr::utils::communication::ParallelAdd( + mf, *J_cp[lev+1][idim], 0, 0, + J_cp[lev+1][idim]->nComp(), + ng, amrex::IntVect(0), + do_single_precision_comms, period); } + SumBoundaryJ(J_cp, lev+1, idim, period); MultiFab::Add(*J_fp[lev][idim], mf, 0, 0, J_fp[lev+1][idim]->nComp(), 0); } } diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H index a807a426296..391682b469e 100644 --- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H +++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H @@ -196,12 +196,33 @@ public: multiplier_ratio = max_N; } +#if (defined WARPX_DIM_RZ) + amrex::ParticleReal * const AMREX_RESTRICT theta1 = soa_1.m_rdata[PIdx::theta]; + amrex::ParticleReal * const AMREX_RESTRICT theta2 = soa_2.m_rdata[PIdx::theta]; +#endif + for (int k = 0; k < max_N; ++k) { // c1k : how many times the current particle of species 1 is paired with a particle // of species 2. Same for c2k. const int c1k = (k%NI1 < max_N%NI1) ? c1 + 1: c1; const int c2k = (k%NI2 < max_N%NI2) ? c2 + 1: c2; + +#if (defined WARPX_DIM_RZ) + /* In RZ geometry, macroparticles can collide with other macroparticles + * in the same *cylindrical* cell. For this reason, collisions between macroparticles + * are actually not local in space. In this case, the underlying assumption is that + * particles within the same cylindrical cell represent a cylindrically-symmetry + * momentum distribution function. Therefore, here, we temporarily rotate the + * momentum of one of the macroparticles in agreement with this cylindrical symmetry. + * (This is technically only valid if we use only the m=0 azimuthal mode in the simulation; + * there is a corresponding assert statement at initialization.) */ + amrex::ParticleReal const theta = theta2[I2[i2]]-theta1[I1[i1]]; + amrex::ParticleReal const u1xbuf = u1x[I1[i1]]; + u1x[I1[i1]] = u1xbuf*std::cos(theta) - u1y[I1[i1]]*std::sin(theta); + u1y[I1[i1]] = u1xbuf*std::sin(theta) + u1y[I1[i1]]*std::cos(theta); +#endif + SingleNuclearFusionEvent( u1x[ I1[i1] ], u1y[ I1[i1] ], u1z[ I1[i1] ], u2x[ I2[i2] ], u2y[ I2[i2] ], u2z[ I2[i2] ], @@ -211,6 +232,13 @@ public: m_probability_threshold, m_probability_target_value, m_fusion_type, engine); + +#if (defined WARPX_DIM_RZ) + amrex::ParticleReal const u1xbuf_new = u1x[I1[i1]]; + u1x[I1[i1]] = u1xbuf_new*std::cos(-theta) - u1y[I1[i1]]*std::sin(-theta); + u1y[I1[i1]] = u1xbuf_new*std::sin(-theta) + u1y[I1[i1]]*std::cos(-theta); +#endif + p_pair_indices_1[pair_index] = I1[i1]; p_pair_indices_2[pair_index] = I2[i2]; ++i1; if ( i1 == static_cast(I1e) ) { i1 = I1s; } diff --git a/Source/Particles/LaserParticleContainer.cpp b/Source/Particles/LaserParticleContainer.cpp index 565ac3ee1b8..b56129dc0b7 100644 --- a/Source/Particles/LaserParticleContainer.cpp +++ b/Source/Particles/LaserParticleContainer.cpp @@ -99,6 +99,13 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, getArrWithParser(pp_laser_name, "direction", m_nvec); getArrWithParser(pp_laser_name, "polarization", m_p_X); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_position.size() == 3, + m_laser_name + ".position must have three components."); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_nvec.size() == 3, + m_laser_name + ".direction must have three components."); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_p_X.size() == 3, + m_laser_name + ".polarization must have three components."); + getWithParser(pp_laser_name, "wavelength", m_wavelength); AMREX_ALWAYS_ASSERT_WITH_MESSAGE( m_wavelength > 0, "The laser wavelength must be >0."); @@ -115,7 +122,7 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, ); pp_laser_name.query("do_continuous_injection", do_continuous_injection); - pp_laser_name.query("min_particles_per_mode", m_min_particles_per_mode); + queryWithParser(pp_laser_name, "min_particles_per_mode", m_min_particles_per_mode); if (m_e_max == amrex::Real(0.)){ ablastr::warn_manager::WMRecordWarning("Laser", @@ -550,6 +557,8 @@ LaserParticleContainer::Evolve (int lev, amrex::LayoutData* cost = WarpX::getCosts(lev); + const bool has_buffer = cjx; + #ifdef AMREX_USE_OMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif @@ -578,18 +587,21 @@ LaserParticleContainer::Evolve (int lev, auto& uzp = attribs[PIdx::uz]; const long np = pti.numParticles(); - // For now, laser particles do not take the current buffers into account - const long np_current = np; - plane_Xp.resize(np); plane_Yp.resize(np); amplitude_E.resize(np); + // Determine whether particles will deposit on the fine or coarse level + long np_current = np; + if (lev > 0 && m_deposit_on_main_grid && has_buffer) { + np_current = 0; + } + if (rho && ! skip_deposition && ! do_not_deposit) { int* AMREX_RESTRICT ion_lev = nullptr; DepositCharge(pti, wp, ion_lev, rho, 0, 0, np_current, thread_num, lev, lev); - if (crho) { + if (has_buffer) { DepositCharge(pti, wp, ion_lev, crho, 0, np_current, np-np_current, thread_num, lev, lev-1); } @@ -628,7 +640,6 @@ LaserParticleContainer::Evolve (int lev, 0, np_current, thread_num, lev, lev, dt, relative_time); - const bool has_buffer = cjx; if (has_buffer) { // Deposit in buffers @@ -643,7 +654,7 @@ LaserParticleContainer::Evolve (int lev, int* AMREX_RESTRICT ion_lev = nullptr; DepositCharge(pti, wp, ion_lev, rho, 1, 0, np_current, thread_num, lev, lev); - if (crho) { + if (has_buffer) { DepositCharge(pti, wp, ion_lev, crho, 1, np_current, np-np_current, thread_num, lev, lev-1); } diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index bf68dc0681f..e2edc4dd71e 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -371,6 +371,7 @@ protected: //! instead of depositing (current, charge) on the finest patch level, deposit to the coarsest grid std::vector m_deposit_on_main_grid; + std::vector m_laser_deposit_on_main_grid; //! instead of gathering fields from the finest patch level, gather from the coarsest std::vector m_gather_from_main_grid; diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index c6cb9032e1d..662e0b45f6d 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -112,6 +112,7 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) for (int i = nspecies; i < nspecies+nlasers; ++i) { allcontainers[i] = std::make_unique(amr_core, i, lasers_names[i-nspecies]); + allcontainers[i]->m_deposit_on_main_grid = m_laser_deposit_on_main_grid[i-nspecies]; } pc_tmp = std::make_unique(amr_core); @@ -347,6 +348,21 @@ MultiParticleContainer::ReadParameters () ParmParse pp_lasers("lasers"); pp_lasers.queryarr("names", lasers_names); + auto const nlasers = lasers_names.size(); + // Get lasers to deposit on main grid + m_laser_deposit_on_main_grid.resize(nlasers, false); + std::vector tmp; + pp_lasers.queryarr("deposit_on_main_grid", tmp); + for (auto const& name : tmp) { + auto it = std::find(lasers_names.begin(), lasers_names.end(), name); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + it != lasers_names.end(), + "laser '" + name + + "' in lasers.deposit_on_main_grid must be part of lasers.lasers_names"); + int i = std::distance(lasers_names.begin(), it); + m_laser_deposit_on_main_grid[i] = true; + } + #ifdef WARPX_QED ParmParse pp_warpx("warpx"); diff --git a/Source/Utils/CMakeLists.txt b/Source/Utils/CMakeLists.txt index b2754cdaf99..9253f3d0d7f 100644 --- a/Source/Utils/CMakeLists.txt +++ b/Source/Utils/CMakeLists.txt @@ -13,3 +13,5 @@ target_sources(WarpX WarpXUtil.cpp WarpXVersion.cpp ) + +add_subdirectory(Logo) diff --git a/Source/Utils/IntervalsParser.H b/Source/Utils/IntervalsParser.H index c91fe271abc..06258b10999 100644 --- a/Source/Utils/IntervalsParser.H +++ b/Source/Utils/IntervalsParser.H @@ -21,7 +21,7 @@ public: * (0 for the starting point, std::numeric_limits::max() for the stopping point and 1 for * the period). For example SliceParser(":1000:") is equivalent to SliceParser("0:1000:1"). */ - SliceParser (const std::string& instr); + SliceParser (const std::string& instr, bool isBTD=false); /** * \brief A method that returns true if the input integer is contained in the slice. (e.g. if @@ -66,7 +66,14 @@ public: */ int getStop () const; + /** + * @brief A method that returns the number of integers contained by the slice. + * + */ + int numContained() const; + private: + bool m_isBTD = false; int m_start = 0; int m_stop = std::numeric_limits::max(); int m_period = 1; @@ -148,4 +155,53 @@ private: bool m_activated = false; }; +/** + * \brief This class is a parser for multiple slices of the form x,y,z,... where x, y and z are + * slices of the form i:j:k, as defined in the SliceParser class. This class contains a vector of + * SliceParsers. The supported function set differs from the IntervalsParser + */ +class BTDIntervalsParser +{ +public: + /** + * \brief Default constructor of the BTDIntervalsParser class. + */ + BTDIntervalsParser () = default; + + /** + * \brief Constructor of the BTDIntervalsParser class. + * + * @param[in] instr_vec an input vector string, which when concatenated is of the form + * "x,y,z,...". This will call the constructor of SliceParser using x, y and z as input + * arguments. + */ + BTDIntervalsParser (const std::vector& instr_vec); + + /** + * @brief Return the total number of unique labframe snapshots + */ + int NumSnapshots (); + + /** + * @brief Return the iteration number stored at index i_buffer + * + * @param i_buffer buffer or iteration index, between 0 and NumSnapshots + */ + int GetBTDIteration(int i_buffer); + + /** + * \brief A method that returns true if any of the slices contained by the IntervalsParser + * has a strictly positive period. + */ + bool isActivated () const; + +private: + std::vector m_btd_iterations; + std::vector m_slices; + std::vector m_slice_starting_i_buffer; + int m_n_snapshots; + static constexpr char m_separator = ','; + bool m_activated = false; +}; + #endif // WARPX_INTERVALSPARSER_H_ diff --git a/Source/Utils/IntervalsParser.cpp b/Source/Utils/IntervalsParser.cpp index e8425421f29..4da0142a03e 100644 --- a/Source/Utils/IntervalsParser.cpp +++ b/Source/Utils/IntervalsParser.cpp @@ -7,15 +7,18 @@ #include #include -SliceParser::SliceParser (const std::string& instr) +SliceParser::SliceParser (const std::string& instr, const bool isBTD) { + m_isBTD = isBTD; // split string and trim whitespaces auto insplit = WarpXUtilStr::split>(instr, m_separator, true); if(insplit.size() == 1){ // no colon in input string. The input is the period. + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(!m_isBTD, "must specify interval stop for BTD"); m_period = parseStringtoInt(insplit[0], "interval period");} else if(insplit.size() == 2) // 1 colon in input string. The input is start:stop { + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(!m_isBTD || !insplit[1].empty(), "must specify interval stop for BTD"); if (!insplit[0].empty()){ m_start = parseStringtoInt(insplit[0], "interval start");} if (!insplit[1].empty()){ @@ -23,6 +26,7 @@ SliceParser::SliceParser (const std::string& instr) } else // 2 colons in input string. The input is start:stop:period { + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(!m_isBTD || !insplit[1].empty(), "must specify interval stop for BTD"); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( insplit.size() == 3, instr + "' is not a valid syntax for a slice."); @@ -64,6 +68,8 @@ int SliceParser::getStart () const {return m_start;} int SliceParser::getStop () const {return m_stop;} +int SliceParser::numContained () const {return (m_stop - m_start) / m_period + 1;} + IntervalsParser::IntervalsParser (const std::vector& instr_vec) { std::string inconcatenated; @@ -116,3 +122,94 @@ int IntervalsParser::localPeriod (const int n) const } bool IntervalsParser::isActivated () const {return m_activated;} + +BTDIntervalsParser::BTDIntervalsParser (const std::vector& instr_vec) +{ + std::string inconcatenated; + for (const auto& instr_element : instr_vec) inconcatenated +=instr_element; + + auto const insplit = WarpXUtilStr::split>(inconcatenated, std::string(1,m_separator)); + + // parse the Intervals string into Slices and store each slice in m_slices, + // in order of increasing Slice start value + for(const auto& inslc : insplit) + { + bool isBTD = true; + SliceParser temp_slice(inslc, isBTD); + if (m_slices.size() > 0) + { + // find the last index i_slice where + // the start value of m_slices[i_slice] is greater than temp_slices' start_value + int i_slice = 0; + while (temp_slice.getStart() > m_slices[i_slice].getStart() && i_slice < static_cast(m_slices.size())) + { + i_slice++; + } + m_slices.insert(m_slices.begin() + i_slice, temp_slice); + } + else + { + m_slices.push_back(temp_slice); + } + } + // from the vector of slices, m_slices, + // create a vector of integers, m_btd_iterations, containing + // the iteration of every back-transformed snapshot that will be saved + // the iteration values in m_btd_iterations are + // 1. saved in increasing order + // 2. unique, i.e. no duplicate iterations are saved + for (const auto& temp_slice : m_slices) + { + const int start = temp_slice.getStart(); + const int period = temp_slice.getPeriod(); + int btd_iter_ind; + // for Slice temp_slice in m_slices, + // determine the index in m_btd_iterations where temp_slice's starting value goes + // + // Implementation note: + // assuming the user mostly lists slices in ascending order, + // start at the end of m_btd_iterations and search backward + if (m_btd_iterations.size() == 0) + { + btd_iter_ind = 0; + } + else + { + btd_iter_ind = m_btd_iterations.size() - 1; + while (start < m_btd_iterations[btd_iter_ind] and btd_iter_ind>0) + { + btd_iter_ind--; + } + } + // insert each iteration contained in temp_slice into m_btd_iterations + // adding them in increasing sorted order and not adding any iterations + // already contained in m_btd_iterations + for (int ii = start; ii <= temp_slice.getStop(); ii += period) + { + if (m_btd_iterations.size() > 0) + { + // find where iteration ii should go in m_btd_iterations + while (ii > m_btd_iterations[btd_iter_ind] && btd_iter_ind < static_cast(m_btd_iterations.size())) + { + btd_iter_ind++; + } + if (ii != m_btd_iterations[btd_iter_ind]) + { + m_btd_iterations.insert(m_btd_iterations.begin() + btd_iter_ind, ii); + } + } else + { + m_btd_iterations.push_back(ii); + } + } + if ((temp_slice.getPeriod() > 0) && + (temp_slice.getStop() >= start)) m_activated = true; + } +} + +int BTDIntervalsParser::NumSnapshots () { return m_btd_iterations.size(); } + +int BTDIntervalsParser::GetBTDIteration(int i_buffer) +{ + return m_btd_iterations[i_buffer]; +} diff --git a/Source/Utils/Logo/CMakeLists.txt b/Source/Utils/Logo/CMakeLists.txt new file mode 100644 index 00000000000..62b5371ee98 --- /dev/null +++ b/Source/Utils/Logo/CMakeLists.txt @@ -0,0 +1,4 @@ +target_sources(WarpX + PRIVATE + GetLogo.cpp +) diff --git a/Source/Utils/Logo/GetLogo.H b/Source/Utils/Logo/GetLogo.H new file mode 100644 index 00000000000..70a3e38ff29 --- /dev/null +++ b/Source/Utils/Logo/GetLogo.H @@ -0,0 +1,19 @@ +#ifndef WARPX_GET_LOGO_H_ +#define WARPX_GET_LOGO_H_ + +#include + +namespace utils +{ + namespace logo + { + /** + * \brief provides an ASCII art logo for WarpX + * + * \return a string containing an ASCII art logo + */ + std::string get_logo () noexcept; + } +} + +#endif diff --git a/Source/Utils/Logo/GetLogo.cpp b/Source/Utils/Logo/GetLogo.cpp new file mode 100644 index 00000000000..c515ceac1a7 --- /dev/null +++ b/Source/Utils/Logo/GetLogo.cpp @@ -0,0 +1,15 @@ +#include "GetLogo.H" + +std::string utils::logo::get_logo () noexcept +{ + return +R"( + __ __ __ __ + \ \ / /_ _ _ __ _ __\ \/ / + \ \ /\ / / _` | '__| '_ \\ / + \ V V / (_| | | | |_) / \ + \_/\_/ \__,_|_| | .__/_/\_\ + |_| + +)"; +} diff --git a/Source/Utils/Logo/Make.package b/Source/Utils/Logo/Make.package new file mode 100644 index 00000000000..2aa69ae46da --- /dev/null +++ b/Source/Utils/Logo/Make.package @@ -0,0 +1,3 @@ +CEXE_sources += GetLogo.cpp + +VPATH_LOCATIONS += $(WARPX_HOME)/Source/Utils/Logo diff --git a/Source/Utils/Make.package b/Source/Utils/Make.package index c3771508390..b20a1abe291 100644 --- a/Source/Utils/Make.package +++ b/Source/Utils/Make.package @@ -11,4 +11,6 @@ CEXE_sources += MPIInitHelpers.cpp CEXE_sources += RelativeCellPosition.cpp CEXE_sources += ParticleUtils.cpp +include $(WARPX_HOME)/Source/Utils/Logo/Make.package + VPATH_LOCATIONS += $(WARPX_HOME)/Source/Utils diff --git a/Source/Utils/WarpXAlgorithmSelection.H b/Source/Utils/WarpXAlgorithmSelection.H index ba752e012b7..97dd3b89fe3 100644 --- a/Source/Utils/WarpXAlgorithmSelection.H +++ b/Source/Utils/WarpXAlgorithmSelection.H @@ -83,6 +83,20 @@ struct GatheringAlgo { }; }; +struct JInTime { + enum { + Constant = 0, + Linear = 1 + }; +}; + +struct RhoInTime { + enum { + Linear = 1, + Quadratic = 2 + }; +}; + /** Strategy to compute weights for use in load balance. */ struct LoadBalanceCostsUpdateAlgo { diff --git a/Source/Utils/WarpXAlgorithmSelection.cpp b/Source/Utils/WarpXAlgorithmSelection.cpp index d0db13e12e7..088d5322f77 100644 --- a/Source/Utils/WarpXAlgorithmSelection.cpp +++ b/Source/Utils/WarpXAlgorithmSelection.cpp @@ -62,6 +62,18 @@ const std::map gathering_algo_to_int = { {"default", GatheringAlgo::EnergyConserving } }; +const std::map J_in_time_to_int = { + {"constant", JInTime::Constant}, + {"linear", JInTime::Linear}, + {"default", JInTime::Constant} +}; + +const std::map rho_in_time_to_int = { + {"linear", RhoInTime::Linear}, + {"quadratic", RhoInTime::Quadratic}, + {"default", RhoInTime::Linear} +}; + const std::map load_balance_costs_update_algo_to_int = { {"timers", LoadBalanceCostsUpdateAlgo::Timers }, {"gpuclock", LoadBalanceCostsUpdateAlgo::GpuClock }, @@ -131,6 +143,10 @@ GetAlgorithmInteger( amrex::ParmParse& pp, const char* pp_search_key ){ algo_to_int = charge_deposition_algo_to_int; } else if (0 == std::strcmp(pp_search_key, "field_gathering")) { algo_to_int = gathering_algo_to_int; + } else if (0 == std::strcmp(pp_search_key, "J_in_time")) { + algo_to_int = J_in_time_to_int; + } else if (0 == std::strcmp(pp_search_key, "rho_in_time")) { + algo_to_int = rho_in_time_to_int; } else if (0 == std::strcmp(pp_search_key, "load_balance_costs_update")) { algo_to_int = load_balance_costs_update_algo_to_int; } else if (0 == std::strcmp(pp_search_key, "em_solver_medium")) { diff --git a/Source/Utils/WarpXConst.H b/Source/Utils/WarpXConst.H index fc40be2908f..d6f016ede4f 100644 --- a/Source/Utils/WarpXConst.H +++ b/Source/Utils/WarpXConst.H @@ -1,5 +1,5 @@ -/* Copyright 2019 Andrew Myers, Luca Fedeli, Maxence Thevenet - * Weiqun Zhang +/* Copyright 2019-2022 Andrew Myers, Luca Fedeli, Maxence Thevenet, + * Weiqun Zhang, Axel Huebl * * This file is part of WarpX. * @@ -8,43 +8,17 @@ #ifndef WARPX_CONST_H_ #define WARPX_CONST_H_ -#include +#include -#include -#include -#include - -// Math constants namespace MathConst { - static constexpr amrex::Real pi = static_cast(3.14159265358979323846); + using namespace ablastr::constant::math; } -// Physical constants. Values are the 2018 CODATA recommended values -// https://physics.nist.gov/cuu/Constants/index.html -// -// New additions here should also be considered for addition to -// `warpx_constants` in WarpXUtil.cpp's `makeParser`, so that they're -// available in parsing and evaluation of PICMI expressions, as well -// as the corresponding Python definitions namespace PhysConst { - static constexpr auto c = static_cast( 299'792'458. ); - static constexpr auto ep0 = static_cast( 8.8541878128e-12 ); - static constexpr auto mu0 = static_cast( 1.25663706212e-06 ); - static constexpr auto q_e = static_cast( 1.602176634e-19 ); - static constexpr auto m_e = static_cast( 9.1093837015e-31 ); - static constexpr auto m_p = static_cast( 1.67262192369e-27 ); - static constexpr auto m_u = static_cast( 1.66053906660e-27 ); - static constexpr auto hbar = static_cast( 1.054571817e-34 ); - static constexpr auto alpha = static_cast( 0.007297352573748943 );//mu0/(4*MathConst::pi)*q_e*q_e*c/hbar; - static constexpr auto r_e = static_cast( 2.817940326204929e-15 );//1./(4*MathConst::pi*ep0) * q_e*q_e/(m_e*c*c); - static constexpr double xi = 1.3050122447005176e-52; //(2.*alpha*alpha*ep0*ep0*hbar*hbar*hbar)/(45.*m_e*m_e*m_e*m_e*c*c*c*c*c); - static constexpr auto xi_c2 = static_cast( 1.1728865132395492e-35 ); // This should be usable for single precision, though - // very close to smallest number possible (1.2e-38) - - static constexpr auto kb = static_cast( 1.380649e-23 ); // Boltzmann's constant, J/K (exact) + using namespace ablastr::constant::SI; } #endif diff --git a/Source/Utils/WarpXUtil.H b/Source/Utils/WarpXUtil.H index 9360c14d459..275ce7e3143 100644 --- a/Source/Utils/WarpXUtil.H +++ b/Source/Utils/WarpXUtil.H @@ -197,6 +197,17 @@ void getCellCoordinates (int i, int j, int k, int safeCastToInt(amrex::Real x, const std::string& real_name); +/** +* \brief Do a safe cast of a real to a long +* This ensures that the float value is within the range of longs and if not, +* raises an exception. +* +* \param x Real value to cast +* \param real_name String, the name of the variable being casted to use in the error message +*/ +long +safeCastToLong(amrex::Real x, const std::string& real_name); + /** * \brief Initialize an amrex::Parser object from a string containing a math expression * @@ -265,6 +276,10 @@ int queryWithParser (const amrex::ParmParse& a_pp, char const * const str, T& va val = safeCastToInt(std::round(parser.compileHost<0>()()), str); } + else if (std::is_same::value) { + + val = safeCastToLong(std::round(parser.compileHost<0>()()), str); + } else { val = static_cast(parser.compileHost<0>()()); } @@ -289,6 +304,9 @@ int queryArrWithParser (const amrex::ParmParse& a_pp, char const * const str, st if (std::is_same::value) { val[i] = safeCastToInt(std::round(parser.compileHost<0>()()), str); } + else if (std::is_same::value) { + val[i] = safeCastToLong(std::round(parser.compileHost<0>()()), str); + } else { val[i] = static_cast(parser.compileHost<0>()()); } @@ -330,6 +348,9 @@ int queryArrWithParser (const amrex::ParmParse& a_pp, char const * const str, st if (std::is_same::value) { val[i] = safeCastToInt(std::round(parser.compileHost<0>()()), str); } + else if (std::is_same::value) { + val[i] = safeCastToLong(std::round(parser.compileHost<0>()()), str); + } else { val[i] = static_cast(parser.compileHost<0>()()); } @@ -361,6 +382,9 @@ void getWithParser (const amrex::ParmParse& a_pp, char const * const str, T& val if (std::is_same::value) { val = safeCastToInt(std::round(parser.compileHost<0>()()), str); } + else if (std::is_same::value) { + val = safeCastToLong(std::round(parser.compileHost<0>()()), str); + } else { val = static_cast(parser.compileHost<0>()()); } @@ -380,6 +404,9 @@ void getArrWithParser (const amrex::ParmParse& a_pp, char const * const str, std if (std::is_same::value) { val[i] = safeCastToInt(std::round(parser.compileHost<0>()()), str); } + else if (std::is_same::value) { + val[i] = safeCastToLong(std::round(parser.compileHost<0>()()), str); + } else { val[i] = static_cast(parser.compileHost<0>()()); } @@ -416,6 +443,9 @@ void getArrWithParser (const amrex::ParmParse& a_pp, char const * const str, std if (std::is_same::value) { val[i] = safeCastToInt(std::round(parser.compileHost<0>()()), str); } + else if (std::is_same::value) { + val[i] = safeCastToLong(std::round(parser.compileHost<0>()()), str); + } else { val[i] = static_cast(parser.compileHost<0>()()); } diff --git a/Source/Utils/WarpXUtil.cpp b/Source/Utils/WarpXUtil.cpp index afe8b7daa07..eecc04b7802 100644 --- a/Source/Utils/WarpXUtil.cpp +++ b/Source/Utils/WarpXUtil.cpp @@ -276,30 +276,43 @@ void Store_parserString(const amrex::ParmParse& pp, std::string query_string, f.clear(); } -int safeCastToInt(const amrex::Real x, const std::string& real_name) { - int result = 0; - bool error_detected = false; - std::string assert_msg; - // (2.0*(numeric_limits::max()/2+1)) converts numeric_limits::max()+1 to a real ensuring accuracy to all digits - // This accepts x = 2**31-1 but rejects 2**31. - using namespace amrex::literals; - constexpr amrex::Real max_range = (2.0_rt*static_cast(std::numeric_limits::max()/2+1)); - if (x < max_range) { - if (std::ceil(x) >= std::numeric_limits::min()) { - result = static_cast(x); +namespace WarpXUtilSafeCast { + template< typename int_type > + AMREX_FORCE_INLINE + int_type safeCastTo(const amrex::Real x, const std::string& real_name) { + int_type result = int_type(0); + bool error_detected = false; + std::string assert_msg; + // (2.0*(numeric_limits::max()/2+1)) converts numeric_limits::max()+1 to a real ensuring accuracy to all digits + // This accepts x = 2**31-1 but rejects 2**31. + using namespace amrex::literals; + constexpr amrex::Real max_range = (2.0_rt*static_cast(std::numeric_limits::max()/2+1)); + if (x < max_range) { + if (std::ceil(x) >= std::numeric_limits::min()) { + result = static_cast(x); + } else { + error_detected = true; + assert_msg = "Negative overflow detected when casting " + real_name + " = " + + std::to_string(x) + " to integer type"; + } + } else if (x > 0) { + error_detected = true; + assert_msg = "Overflow detected when casting " + real_name + " = " + std::to_string(x) + " to integer type"; } else { error_detected = true; - assert_msg = "Negative overflow detected when casting " + real_name + " = " + std::to_string(x) + " to int"; + assert_msg = "NaN detected when casting " + real_name + " to integer type"; } - } else if (x > 0) { - error_detected = true; - assert_msg = "Overflow detected when casting " + real_name + " = " + std::to_string(x) + " to int"; - } else { - error_detected = true; - assert_msg = "NaN detected when casting " + real_name + " to int"; - } - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(!error_detected, assert_msg); - return result; + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(!error_detected, assert_msg); + return result; + } +} + +int safeCastToInt(const amrex::Real x, const std::string& real_name) { + return WarpXUtilSafeCast::safeCastTo (x, real_name); +} + +long safeCastToLong(const amrex::Real x, const std::string& real_name) { + return WarpXUtilSafeCast::safeCastTo (x, real_name); } Parser makeParser (std::string const& parse_function, amrex::Vector const& varnames) diff --git a/Source/WarpX.H b/Source/WarpX.H index 7a051e30db8..d4acf916e39 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -186,6 +186,10 @@ public: */ static amrex::Vector particle_boundary_hi; + //! Integers that correspond to the time dependency of J (constant, linear) + //! and rho (linear, quadratic) for the PSATD algorithm + static short J_in_time; + static short rho_in_time; //! If true, the current is deposited on a nodal grid and then centered onto a staggered grid //! using finite centering of order given by #current_centering_nox, #current_centering_noy, @@ -804,6 +808,8 @@ public: void computeB (amrex::Vector, 3> >& B, const amrex::Vector >& phi, std::array const beta = {{0,0,0}} ) const; + void computePhiTriDiagonal (const amrex::Vector >& rho, + amrex::Vector >& phi) const; /** * \brief @@ -974,11 +980,22 @@ private: const int lev); void StoreCurrent (const int lev); void RestoreCurrent (const int lev); - void ApplyFilterandSumBoundaryJ ( - const amrex::Vector,3>>& J_fp, - const amrex::Vector,3>>& J_cp, + void ApplyFilterJ ( + const amrex::Vector,3>>& current, const int lev, - PatchType patch_type); + const int idim); + void ApplyFilterJ ( + const amrex::Vector,3>>& current, + const int lev); + void SumBoundaryJ ( + const amrex::Vector,3>>& current, + const int lev, + const int idim, + const amrex::Periodicity& period); + void SumBoundaryJ ( + const amrex::Vector,3>>& current, + const int lev, + const amrex::Periodicity& period); void NodalSyncJ ( const amrex::Vector,3>>& J_fp, const amrex::Vector,3>>& J_cp, @@ -1094,7 +1111,7 @@ private: const int centering_noz); void AllocLevelMFs (int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm, - const amrex::IntVect& ngEB, const amrex::IntVect& ngJ, + const amrex::IntVect& ngEB, amrex::IntVect& ngJ, const amrex::IntVect& ngRho, const amrex::IntVect& ngF, const amrex::IntVect& ngG, const bool aux_is_nodal); #ifdef WARPX_USE_PSATD diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 86961dcaf31..aa6301a6e4b 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -121,6 +121,8 @@ short WarpX::charge_deposition_algo; short WarpX::field_gathering_algo; short WarpX::particle_pusher_algo; short WarpX::maxwell_solver_id; +short WarpX::J_in_time; +short WarpX::rho_in_time; short WarpX::load_balance_costs_update_algo; bool WarpX::do_dive_cleaning = false; bool WarpX::do_divb_cleaning = false; @@ -1154,6 +1156,11 @@ WarpX::ReadParameters () WARPX_ALWAYS_ASSERT_WITH_MESSAGE(noz_fft > 0, "PSATD order must be finite unless psatd.periodic_single_box_fft is used"); } + // Integers that correspond to the time dependency of J (constant, linear) + // and rho (linear, quadratic) for the PSATD algorithm + J_in_time = GetAlgorithmInteger(pp_psatd, "J_in_time"); + rho_in_time = GetAlgorithmInteger(pp_psatd, "rho_in_time"); + // Current correction activated by default, unless a charge-conserving // current deposition (Esirkepov, Vay) or the div(E) cleaning scheme // are used @@ -1308,10 +1315,28 @@ WarpX::ReadParameters () v_galilean_is_zero, "Multi-J algorithm not implemented with Galilean PSATD" ); + } - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(update_with_rho, - "psatd.update_with_rho must be set to 1 when warpx.do_multi_J = 1" - ); + if (J_in_time == JInTime::Constant) + { + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + rho_in_time == RhoInTime::Linear, + "psatd.J_in_time=constant supports only psatd.rho_in_time=linear"); + } + + if (J_in_time == JInTime::Linear) + { + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + update_with_rho, + "psatd.update_with_rho must be set to 1 when psatd.J_in_time=linear"); + + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + v_galilean_is_zero, + "psatd.J_in_time=linear not implemented with Galilean PSATD"); + + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + v_comoving_is_zero, + "psatd.J_in_time=linear not implemented with comoving PSATD"); } for (int dir = 0; dir < AMREX_SPACEDIM; dir++) @@ -1626,6 +1651,13 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d amrex::RealVect dx = {WarpX::CellSize(lev)[0], WarpX::CellSize(lev)[1], WarpX::CellSize(lev)[2]}; #endif + // Initialize filter before guard cells manager + // (needs info on length of filter's stencil) + if (use_filter) + { + InitFilter(); + } + guard_cells.Init( dt[lev], dx, @@ -1648,7 +1680,9 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d WarpX::isAnyBoundaryPML(), WarpX::do_pml_in_domain, WarpX::pml_ncell, - this->refRatio()); + this->refRatio(), + use_filter, + bilinear_filter.stencil_length_each_dir); #ifdef AMREX_USE_EB @@ -1683,7 +1717,7 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d void WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm, - const IntVect& ngEB, const IntVect& ngJ, const IntVect& ngRho, + const IntVect& ngEB, IntVect& ngJ, const IntVect& ngRho, const IntVect& ngF, const IntVect& ngG, const bool aux_is_nodal) { // Declare nodal flags @@ -2207,7 +2241,8 @@ void WarpX::AllocLevelSpectralSolverRZ (amrex::Vector + + +/** Numerical compile-time constants */ +namespace ablastr::constant +{ + /** Mathematical constants */ + namespace math + { + using namespace amrex::literals; + + //! ratio of a circle's circumference to its diameter + static constexpr amrex::Real pi = 3.14159265358979323846_rt; + + //! https://tauday.com/tau-manifesto + static constexpr amrex::Real tau = 2.0_rt * pi; + } // namespace math + + /** Physical constants + * + * Values are the 2018 CODATA recommended values + * https://physics.nist.gov/cuu/Constants/index.html + * + * New additions here should also be considered for addition to + * `warpx_constants` in WarpXUtil.cpp's `makeParser`, so that they're + * available in parsing and evaluation of PICMI expressions, as well + * as the corresponding Python definitions + */ + namespace SI + { + using namespace amrex::literals; + + //! vacuum speed of light [m/s] + static constexpr auto c = 299'792'458._rt; + //! vacuum permittivity: dielectric permittivity of vacuum [F/m] + static constexpr auto ep0 = 8.8541878128e-12_rt; + //! vacuum permeability: magnetic permeability of vacuum = 4.0e-7 * pi [H/m] + static constexpr auto mu0 = 1.25663706212e-06_rt; + //! elementary charge [C] + static constexpr auto q_e = 1.602176634e-19_rt; + //! electron mass [kg] + static constexpr auto m_e = 9.1093837015e-31_rt; + //! proton mass [kg] + static constexpr auto m_p = 1.67262192369e-27_rt; + //! dalton: unified atomic mass unit [kg] + static constexpr auto m_u = 1.66053906660e-27_rt; + + //! reduced Planck Constant = h / tau [J*s] + static constexpr auto hbar = 1.054571817e-34_rt; + //! fine-structure constant = mu0/(4*MathConst::pi)*q_e*q_e*c/hbar [dimensionless] + static constexpr auto alpha = 0.007297352573748943_rt; + //! classical electron radius = 1./(4*MathConst::pi*ep0) * q_e*q_e/(m_e*c*c) [m] + static constexpr auto r_e = 2.817940326204929e-15_rt; + //! xi: nonlinearity parameter of Heisenberg-Euler effective theory = (2.*alpha*alpha*ep0*ep0*hbar*hbar*hbar)/(45.*m_e*m_e*m_e*m_e*c*c*c*c*c) + static constexpr double xi = 1.3050122447005176e-52; + //! xi times c2 = xi*c*c. This should be usable for single precision instead of xi; very close to smallest float32 number possible (1.2e-38) + static constexpr auto xi_c2 = 1.1728865132395492e-35_rt; + + //! Boltzmann constant (exact) [J/K] + static constexpr auto kb = 1.380649e-23_rt; + + //! 1 eV in [J] + static constexpr auto eV = q_e; + //! 1 MeV in [J] + static constexpr auto MeV = q_e * 1e6_rt; + //! 1 eV/c in [kg*m/s] + static constexpr auto eV_invc = eV / c; + //! 1 MeV/c in [kg*m/s] + static constexpr auto MeV_invc = MeV / c; + //! 1 eV/c^2 in [kg] + static constexpr auto eV_invc2 = eV / (c * c); + //! 1 MeV/c^2 in [kg] + static constexpr auto MeV_invc2 = MeV / (c * c); + } // namespace SI +} // namespace ablastr::constant + +#endif // ABLASTR_CONSTANT_H_ diff --git a/Source/ablastr/utils/CMakeLists.txt b/Source/ablastr/utils/CMakeLists.txt index 9ac122caea4..29e7cfb6d20 100644 --- a/Source/ablastr/utils/CMakeLists.txt +++ b/Source/ablastr/utils/CMakeLists.txt @@ -1,8 +1,9 @@ target_sources(ablastr PRIVATE Communication.cpp - TextMsg.cpp SignalHandling.cpp + TextMsg.cpp + UsedInputsFile.cpp ) add_subdirectory(msg_logger) diff --git a/Source/ablastr/utils/Make.package b/Source/ablastr/utils/Make.package index c9be0153acc..6db5f150b94 100644 --- a/Source/ablastr/utils/Make.package +++ b/Source/ablastr/utils/Make.package @@ -1,6 +1,7 @@ -CEXE_sources += TextMsg.cpp -CEXE_sources += SignalHandling.cpp CEXE_sources += Communication.cpp +CEXE_sources += SignalHandling.cpp +CEXE_sources += TextMsg.cpp +CEXE_sources += UsedInputsFile.cpp VPATH_LOCATIONS += $(WARPX_HOME)/Source/ablastr/utils diff --git a/Source/ablastr/utils/UsedInputsFile.H b/Source/ablastr/utils/UsedInputsFile.H new file mode 100644 index 00000000000..543ae98a256 --- /dev/null +++ b/Source/ablastr/utils/UsedInputsFile.H @@ -0,0 +1,27 @@ +/* Copyright 2022 Axel Huebl + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#ifndef ABLASTR_USED_INPUTS_FILE_H +#define ABLASTR_USED_INPUTS_FILE_H + +#include + + +namespace ablastr::utils +{ + /** Write a file that record all inputs: inputs file + command line options + * + * This uses the same syntax as amrex::ParmParse inputs files. + * Only the AMReX IOProcessor writes. + * + * @param filename the name of the text file to write + */ + void + write_used_inputs_file (std::string const & filename); +} + +#endif // ABLASTR_USED_INPUTS_FILE_H diff --git a/Source/ablastr/utils/UsedInputsFile.cpp b/Source/ablastr/utils/UsedInputsFile.cpp new file mode 100644 index 00000000000..175c67619e7 --- /dev/null +++ b/Source/ablastr/utils/UsedInputsFile.cpp @@ -0,0 +1,30 @@ +/* Copyright 2022 Axel Huebl + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#include "UsedInputsFile.H" + +#include +#include +#include + +#include +#include +#include + + +void +ablastr::utils::write_used_inputs_file (std::string const & filename) +{ + amrex::Print() << "For full input parameters, see the file: " << filename << "\n\n"; + + if (amrex::ParallelDescriptor::IOProcessor()) { + std::ofstream jobInfoFile; + jobInfoFile.open(filename.c_str(), std::ios::out); + amrex::ParmParse::dumpTable(jobInfoFile, true); + jobInfoFile.close(); + } +} diff --git a/Tools/machines/crusher-olcf/crusher_warpx.profile.example b/Tools/machines/crusher-olcf/crusher_warpx.profile.example index eeac15cdeec..197f360c2e8 100644 --- a/Tools/machines/crusher-olcf/crusher_warpx.profile.example +++ b/Tools/machines/crusher-olcf/crusher_warpx.profile.example @@ -1,12 +1,14 @@ # please set your project account +# note: WarpX ECP members use aph114_crusher #export proj= # required dependencies -module load cmake/3.22.1 +module load cpe/22.08 +module load cmake/3.23.2 module load craype-accel-amd-gfx90a -module load rocm/5.0.0 +module load rocm/5.2.0 module load cray-mpich -#module load cce/14.0.0 # must be loaded after rocm +module load cce/14.0.2 # must be loaded after rocm # optional: faster builds module load ccache @@ -18,16 +20,20 @@ module load nano # optional: for PSATD in RZ geometry support (not yet available) #module load blaspp #module load lapackpp +export CMAKE_PREFIX_PATH=$HOME/sw/crusher/blaspp-master:$CMAKE_PREFIX_PATH +export CMAKE_PREFIX_PATH=$HOME/sw/crusher/lapackpp-master:$CMAKE_PREFIX_PATH +export LD_LIBRARY_PATH=$HOME/sw/crusher/blaspp-master/lib64:$LD_LIBRARY_PATH +export LD_LIBRARY_PATH=$HOME/sw/crusher/lapackpp-master/lib64:$LD_LIBRARY_PATH # optional: for QED lookup table generation support -module load boost/1.78.0-cxx17 +#module load boost/1.78.0-cxx17 # optional: for openPMD support -module load adios2/2.7.1 -module load cray-hdf5-parallel/1.12.1.1 +module load adios2/2.8.1 +module load cray-hdf5-parallel/1.12.1.5 # optional: for Python bindings or libEnsemble -module load cray-python/3.9.4.2 +module load cray-python/3.9.12.1 # fix system defaults: do not escape $ with a \ on tab completion shopt -s direxpand @@ -51,4 +57,4 @@ export CXX=$(which CC) export FC=$(which ftn) export CFLAGS="-I${ROCM_PATH}/include" export CXXFLAGS="-I${ROCM_PATH}/include" -export LDFLAGS="-L${ROCM_PATH}/lib -lamdhip64" +#export LDFLAGS="-L${ROCM_PATH}/lib -lamdhip64" diff --git a/Tools/machines/crusher-olcf/submit.sh b/Tools/machines/crusher-olcf/submit.sh index 9fbb7748d0a..cefd25c776e 100644 --- a/Tools/machines/crusher-olcf/submit.sh +++ b/Tools/machines/crusher-olcf/submit.sh @@ -1,6 +1,7 @@ #!/usr/bin/env bash #SBATCH -A +# note: WarpX ECP members use aph114 #SBATCH -J warpx #SBATCH -o %x-%j.out #SBATCH -t 00:10:00 @@ -19,10 +20,17 @@ # (GCDs) for a total of 8 GCDs per node. The programmer can think of the 8 GCDs # as 8 separate GPUs, each having 64 GB of high-bandwidth memory (HBM2E). -# note (5-16-22) +# note (5-16-22, OLCFHELP-6888) # this environment setting is currently needed on Crusher to work-around a # known issue with Libfabric -export FI_MR_CACHE_MAX_COUNT=0 +#export FI_MR_CACHE_MAX_COUNT=0 # libfabric disable caching +# or, less invasive: +export FI_MR_CACHE_MONITOR=memhooks # alternative cache monitor + +# note (9-2-22, OLCFDEV-1079) +# this environment setting is needed to avoid that rocFFT writes a cache in +# the home directory, which does not scale. +export ROCFFT_RTC_CACHE_PATH=/dev/null export OMP_NUM_THREADS=8 srun ./warpx inputs > output.txt diff --git a/Tools/machines/frontier-olcf/submit.sh b/Tools/machines/frontier-olcf/submit.sh index d33b4658b83..f5b52a712d3 100644 --- a/Tools/machines/frontier-olcf/submit.sh +++ b/Tools/machines/frontier-olcf/submit.sh @@ -26,7 +26,14 @@ # note (5-16-22 and 7-12-22) # this environment setting is currently needed on Frontier to work-around a # known issue with Libfabric (both in the May and June PE) -export FI_MR_CACHE_MAX_COUNT=0 +#export FI_MR_CACHE_MAX_COUNT=0 # libfabric disable caching +# or, less invasive: +export FI_MR_CACHE_MONITOR=memhooks # alternative cache monitor + +# note (9-2-22, OLCFDEV-1079) +# this environment setting is needed to avoid that rocFFT writes a cache in +# the home directory, which does not scale. +export ROCFFT_RTC_CACHE_PATH=/dev/null export OMP_NUM_THREADS=1 export WARPX_NMPI_PER_NODE=8 diff --git a/Tools/machines/lassen-llnl/lassen_warpx.profile.example b/Tools/machines/lassen-llnl/lassen_warpx.profile.example index 293b60258d6..61e8f971689 100644 --- a/Tools/machines/lassen-llnl/lassen_warpx.profile.example +++ b/Tools/machines/lassen-llnl/lassen_warpx.profile.example @@ -13,7 +13,7 @@ module load fftw/3.3.8 module load boost/1.70.0 # optional: for openPMD support -module load hdf5-parallel/1.10.4 +module load hdf5-parallel/1.12.2 export CMAKE_PREFIX_PATH=$HOME/sw/lassen/c-blosc-1.21.1:$CMAKE_PREFIX_PATH export CMAKE_PREFIX_PATH=$HOME/sw/lassen/adios2-2.7.1:$CMAKE_PREFIX_PATH export LD_LIBRARY_PATH=$HOME/sw/lassen/c-blosc-1.21.1/lib64:$LD_LIBRARY_PATH diff --git a/Tools/machines/summit-olcf/summit_warpx.profile.example b/Tools/machines/summit-olcf/summit_warpx.profile.example index faef6b2145a..216fe9b6f10 100644 --- a/Tools/machines/summit-olcf/summit_warpx.profile.example +++ b/Tools/machines/summit-olcf/summit_warpx.profile.example @@ -13,11 +13,10 @@ module load cuda/11.3.1 module load ccache # optional: for PSATD in RZ geometry support -module load blaspp/2021.04.01-cpu -module load lapackpp/2021.04.00-cpu - -# optional: for PSATD support (CPU only) -#module load fftw/3.3.9 +export CMAKE_PREFIX_PATH=$HOME/sw/summit/blaspp-master:$CMAKE_PREFIX_PATH +export CMAKE_PREFIX_PATH=$HOME/sw/summit/lapackpp-master:$CMAKE_PREFIX_PATH +export LD_LIBRARY_PATH=$HOME/sw/summit/blaspp-master/lib64:$LD_LIBRARY_PATH +export LD_LIBRARY_PATH=$HOME/sw/summit/lapackpp-master/lib64:$LD_LIBRARY_PATH # optional: for QED lookup table generation support module load boost/1.76.0 diff --git a/cmake/dependencies/AMReX.cmake b/cmake/dependencies/AMReX.cmake index 95cfcb6a50e..55cce8ce15e 100644 --- a/cmake/dependencies/AMReX.cmake +++ b/cmake/dependencies/AMReX.cmake @@ -221,7 +221,7 @@ macro(find_amrex) endif() set(COMPONENT_PRECISION ${WarpX_PRECISION} P${WarpX_PARTICLE_PRECISION}) - find_package(AMReX 22.09 CONFIG REQUIRED COMPONENTS ${COMPONENT_ASCENT} ${COMPONENT_DIM} ${COMPONENT_EB} PARTICLES ${COMPONENT_PIC} ${COMPONENT_PRECISION} ${COMPONENT_SENSEI} TINYP LSOLVERS) + find_package(AMReX 22.10 CONFIG REQUIRED COMPONENTS ${COMPONENT_ASCENT} ${COMPONENT_DIM} ${COMPONENT_EB} PARTICLES ${COMPONENT_PIC} ${COMPONENT_PRECISION} ${COMPONENT_SENSEI} TINYP LSOLVERS) message(STATUS "AMReX: Found version '${AMReX_VERSION}'") endif() endmacro() @@ -235,7 +235,7 @@ set(WarpX_amrex_src "" set(WarpX_amrex_repo "https://github.com/AMReX-Codes/amrex.git" CACHE STRING "Repository URI to pull and build AMReX from if(WarpX_amrex_internal)") -set(WarpX_amrex_branch "22.09" +set(WarpX_amrex_branch "13aa4df0f5a4af40270963ad5b42ac7ce662e045" CACHE STRING "Repository branch for WarpX_amrex_repo if(WarpX_amrex_internal)") diff --git a/cmake/dependencies/PICSAR.cmake b/cmake/dependencies/PICSAR.cmake index bd1ab53aa56..3ecee014139 100644 --- a/cmake/dependencies/PICSAR.cmake +++ b/cmake/dependencies/PICSAR.cmake @@ -82,7 +82,7 @@ function(find_picsar) #message(STATUS "PICSAR: Using version '${PICSAR_VERSION}'") else() # not supported by PICSAR (yet) - #find_package(PICSAR 22.09 CONFIG REQUIRED QED) + #find_package(PICSAR 22.10 CONFIG REQUIRED QED) #message(STATUS "PICSAR: Found version '${PICSAR_VERSION}'") message(FATAL_ERROR "PICSAR: Cannot be used as externally installed " "library yet. " diff --git a/requirements.txt b/requirements.txt index d3500be184f..9d6f5bc8840 100644 --- a/requirements.txt +++ b/requirements.txt @@ -4,7 +4,7 @@ periodictable~=1.5 # PICMI # note: don't forget to update the version in Docs/requirements.txt, too -picmistandard==0.0.19 +picmistandard==0.0.20 # for development against an unreleased PICMI version, use: #picmistandard @ git+https://github.com/picmi-standard/picmi.git#subdirectory=PICMI_Python diff --git a/run_test.sh b/run_test.sh index a2f78829ae3..b4e76b8a414 100755 --- a/run_test.sh +++ b/run_test.sh @@ -71,7 +71,7 @@ python3 -m pip install --upgrade -r warpx/Regression/requirements.txt # Clone AMReX and warpx-data git clone https://github.com/AMReX-Codes/amrex.git -cd amrex && git checkout --detach 22.09 && cd - +cd amrex && git checkout --detach 13aa4df0f5a4af40270963ad5b42ac7ce662e045 && cd - # warpx-data contains various required data sets git clone --depth 1 https://github.com/ECP-WarpX/warpx-data.git diff --git a/setup.py b/setup.py index 2ddbf286dd4..695e7328d6a 100644 --- a/setup.py +++ b/setup.py @@ -272,7 +272,7 @@ def build_extension(self, ext): setup( name='pywarpx', # note PEP-440 syntax: x.y.zaN but x.y.z.devN - version = '22.09', + version = '22.10', packages = ['pywarpx'], package_dir = {'pywarpx': 'Python/pywarpx'}, author='Jean-Luc Vay, David P. Grote, Maxence Thévenet, Rémi Lehe, Andrew Myers, Weiqun Zhang, Axel Huebl, et al.',