diff --git a/.github/workflows/codeql.yml b/.github/workflows/codeql.yml index 5c36b9d9f21..e3549ae340a 100644 --- a/.github/workflows/codeql.yml +++ b/.github/workflows/codeql.yml @@ -31,6 +31,11 @@ jobs: - name: Checkout uses: actions/checkout@v4 + - uses: actions/setup-python@v5 + name: Install Python + with: + python-version: '3.x' + - name: Install Packages (C++) if: ${{ matrix.language == 'cpp' }} run: | @@ -38,9 +43,10 @@ jobs: sudo apt-get install --yes cmake openmpi-bin libopenmpi-dev libhdf5-openmpi-dev libadios-openmpi-dev ccache python -m pip install --upgrade pip + python -m pip install --upgrade pipx python -m pip install --upgrade wheel python -m pip install --upgrade cmake - export CMAKE="$HOME/.local/bin/cmake" && echo "CMAKE=$CMAKE" >> $GITHUB_ENV + python -m pipx install cmake - name: Set Up Cache if: ${{ matrix.language == 'cpp' }} @@ -54,7 +60,7 @@ jobs: - name: Configure (C++) if: ${{ matrix.language == 'cpp' }} run: | - $CMAKE -S . -B build -DWarpX_OPENPMD=ON + cmake -S . -B build -DWarpX_OPENPMD=ON - name: Initialize CodeQL uses: github/codeql-action/init@v3 @@ -75,7 +81,7 @@ jobs: export CCACHE_MAXSIZE=100M ccache -z - $CMAKE --build build -j 4 + cmake --build build -j 4 ccache -s du -hs ~/.cache/ccache @@ -83,7 +89,7 @@ jobs: # Make sure CodeQL has something to do touch Source/Utils/WarpXVersion.cpp export CCACHE_DISABLE=1 - $CMAKE --build build -j 4 + cmake --build build -j 4 - name: Perform CodeQL Analysis uses: github/codeql-action/analyze@v3 diff --git a/.github/workflows/cuda.yml b/.github/workflows/cuda.yml index 4a38872a5f4..2209f425d1f 100644 --- a/.github/workflows/cuda.yml +++ b/.github/workflows/cuda.yml @@ -135,7 +135,7 @@ jobs: which nvcc || echo "nvcc not in PATH!" git clone https://github.com/AMReX-Codes/amrex.git ../amrex - cd ../amrex && git checkout --detach 24.10 && cd - + cd ../amrex && git checkout --detach e1222803739ed2342b9ff6fc2d57316ff0d6cb0c && cd - make COMP=gcc QED=FALSE USE_MPI=TRUE USE_GPU=TRUE USE_OMP=FALSE USE_FFT=TRUE USE_CCACHE=TRUE -j 4 ccache -s diff --git a/.github/workflows/intel.yml b/.github/workflows/intel.yml index 1365fa76865..f27181c2e20 100644 --- a/.github/workflows/intel.yml +++ b/.github/workflows/intel.yml @@ -184,6 +184,7 @@ jobs: -DCMAKE_VERBOSE_MAKEFILE=ON \ -DWarpX_COMPUTE=SYCL \ -DWarpX_EB=ON \ + -DWarpX_FFT=ON \ -DWarpX_PYTHON=ON \ -DWarpX_MPI=OFF \ -DWarpX_OPENPMD=ON \ diff --git a/.github/workflows/macos.yml b/.github/workflows/macos.yml index 463b2dc2501..0afaf6ea451 100644 --- a/.github/workflows/macos.yml +++ b/.github/workflows/macos.yml @@ -22,13 +22,16 @@ jobs: #CMAKE_GENERATOR: Ninja steps: - uses: actions/checkout@v4 - - name: install dependencies + - uses: actions/setup-python@v5 + name: Install Python + with: + python-version: '3.x' + - name: install brew dependencies run: | set +e brew unlink gcc brew update brew upgrade || true - brew install --overwrite python brew install ccache brew install fftw brew install libomp @@ -39,12 +42,12 @@ jobs: set -e brew tap openpmd/openpmd brew install openpmd-api - - python3 -m venv py-venv - source py-venv/bin/activate + - name: install pip dependencies + run: | python3 -m pip install --upgrade pip python3 -m pip install --upgrade build packaging setuptools wheel python3 -m pip install --upgrade mpi4py + python3 -m pip install --upgrade -r Regression/requirements.txt - name: CCache Cache uses: actions/cache@v4 with: @@ -60,8 +63,6 @@ jobs: export CCACHE_SLOPPINESS=time_macros ccache -z - source py-venv/bin/activate - cmake -S . -B build_dp \ -DCMAKE_VERBOSE_MAKEFILE=ON \ -DWarpX_EB=OFF \ @@ -71,7 +72,6 @@ jobs: cmake -S . -B build_sp \ -DCMAKE_VERBOSE_MAKEFILE=ON \ - -DPython_EXECUTABLE=$(which python3) \ -DWarpX_EB=OFF \ -DWarpX_PYTHON=ON \ -DWarpX_OPENPMD=ON \ @@ -85,7 +85,6 @@ jobs: - name: run pywarpx run: | - source py-venv/bin/activate export OMP_NUM_THREADS=1 mpirun -n 2 Examples/Physics_applications/laser_acceleration/inputs_test_3d_laser_acceleration_picmi.py diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index d2b15b8af95..8ba600be560 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -18,7 +18,7 @@ exclude: '^share/openPMD/thirdParty' # See https://pre-commit.com/hooks.html for more hooks repos: - repo: https://github.com/pre-commit/pre-commit-hooks - rev: v4.6.0 + rev: v5.0.0 hooks: - id: trailing-whitespace args: [--markdown-linebreak-ext=md] @@ -69,7 +69,7 @@ repos: # Python: Ruff linter & formatter # https://docs.astral.sh/ruff/ - repo: https://github.com/astral-sh/ruff-pre-commit - rev: v0.6.8 + rev: v0.6.9 hooks: # Run the linter - id: ruff diff --git a/CMakeLists.txt b/CMakeLists.txt index 7b1770cfbf7..75999baaa0f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -724,9 +724,9 @@ endforeach() # if(WarpX_PYTHON) set(PY_PIP_OPTIONS "-v" CACHE STRING - "Additional parameters to pass to `pip`") + "Additional parameters to pass to `pip` as ; separated list") set(PY_PIP_INSTALL_OPTIONS "" CACHE STRING - "Additional parameters to pass to `pip install`") + "Additional parameters to pass to `pip install` as ; separated list") # ensure all targets are built before we package them in a wheel set(pyWarpX_INSTALL_TARGET_NAMES) @@ -749,7 +749,8 @@ if(WarpX_PYTHON) ${CMAKE_COMMAND} -E rm -f -r warpx-whl COMMAND ${CMAKE_COMMAND} -E env PYWARPX_LIB_DIR=$ - ${Python_EXECUTABLE} -m pip ${PY_PIP_OPTIONS} wheel --no-build-isolation --no-deps --wheel-dir=warpx-whl ${WarpX_SOURCE_DIR} + ${Python_EXECUTABLE} -m pip ${PY_PIP_OPTIONS} wheel --no-build-isolation --no-deps --wheel-dir=warpx-whl "${WarpX_SOURCE_DIR}" + COMMAND_EXPAND_LISTS VERBATIM WORKING_DIRECTORY ${WarpX_BINARY_DIR} DEPENDS @@ -764,6 +765,7 @@ if(WarpX_PYTHON) endif() add_custom_target(${WarpX_CUSTOM_TARGET_PREFIX}pip_install_requirements ${Python_EXECUTABLE} -m pip ${PY_PIP_OPTIONS} install ${PY_PIP_INSTALL_OPTIONS} -r "${WarpX_SOURCE_DIR}/${pyWarpX_REQUIREMENT_FILE}" + COMMAND_EXPAND_LISTS VERBATIM WORKING_DIRECTORY ${WarpX_BINARY_DIR} ) @@ -781,6 +783,7 @@ if(WarpX_PYTHON) add_custom_target(${WarpX_CUSTOM_TARGET_PREFIX}pip_install ${CMAKE_COMMAND} -E env WARPX_MPI=${WarpX_MPI} ${Python_EXECUTABLE} -m pip ${PY_PIP_OPTIONS} install --force-reinstall --no-index --no-deps ${PY_PIP_INSTALL_OPTIONS} --find-links=warpx-whl pywarpx + COMMAND_EXPAND_LISTS VERBATIM WORKING_DIRECTORY ${WarpX_BINARY_DIR} DEPENDS @@ -794,6 +797,7 @@ if(WarpX_PYTHON) add_custom_target(${WarpX_CUSTOM_TARGET_PREFIX}pip_install_nodeps ${CMAKE_COMMAND} -E env WARPX_MPI=${WarpX_MPI} ${Python_EXECUTABLE} -m pip ${PY_PIP_OPTIONS} install --force-reinstall --no-index --no-deps ${PY_PIP_INSTALL_OPTIONS} --find-links=warpx-whl pywarpx + COMMAND_EXPAND_LISTS VERBATIM WORKING_DIRECTORY ${WarpX_BINARY_DIR} DEPENDS diff --git a/Docs/source/developers/checksum.rst b/Docs/source/developers/checksum.rst index 2452d074ba1..ccbea3408ef 100644 --- a/Docs/source/developers/checksum.rst +++ b/Docs/source/developers/checksum.rst @@ -1,32 +1,36 @@ .. _developers-checksum: -Checksum regression tests -========================= +Checksums on Tests +================== -WarpX has checksum regression tests: as part of CI testing, when running a given test, the checksum module computes one aggregated number per field (``Ex_checksum = np.sum(np.abs(Ex))``) and compares it to a reference (benchmark). This should be sensitive enough to make the test fail if your PR causes a significant difference, print meaningful error messages, and give you a chance to fix a bug or reset the benchmark if needed. +When running an automated test, we often compare the data of final time step of the test with expected values to catch accidental changes. +Instead of relying on reference files that we would have to store in their full size, we calculate an aggregate checksum. -The checksum module is located in ``Regression/Checksum/``, and the benchmarks are stored as human-readable `JSON `__ files in ``Regression/Checksum/benchmarks_json/``, with one file per benchmark (for instance, test ``Langmuir_2d`` has a corresponding benchmark ``Regression/Checksum/benchmarks_json/Langmuir_2d.json``). +For this purpose, the checksum Python module computes one aggregated number per field (e.g., the sum of the absolute values of the array elements) and compares it to a reference value (benchmark). +This should be sensitive enough to make the test fail if your PR causes a significant difference, print meaningful error messages, and give you a chance to fix a bug or reset the benchmark if needed. -For more details on the implementation, the Python files in ``Regression/Checksum/`` should be well documented. +The checksum module is located in ``Regression/Checksum/``, and the benchmarks are stored as human-readable `JSON `__ files in ``Regression/Checksum/benchmarks_json/``, with one file per benchmark (for example, the test ``test_2d_langmuir_multi`` has a corresponding benchmark ``Regression/Checksum/benchmarks_json/test_2d_langmuir_multi.json``). -From a user point of view, you should only need to use ``checksumAPI.py``. It contains Python functions that can be imported and used from an analysis Python script. It can also be executed directly as a Python script. Here are recipes for the main tasks related to checksum regression tests in WarpX CI. +For more details on the implementation, please refer to the Python implementation in ``Regression/Checksum/``. -Include a checksum regression test in an analysis Python script ---------------------------------------------------------------- +From a user point of view, you should only need to use ``checksumAPI.py``, which contains Python functions that can be imported and used from an analysis Python script or can also be executed directly as a Python script. + +How to compare checksums in your analysis script +------------------------------------------------ This relies on the function ``evaluate_checksum``: .. autofunction:: checksumAPI.evaluate_checksum -For an example, see +Here's an example: -.. literalinclude:: ../../../Examples/analysis_default_regression.py +.. literalinclude:: ../../../Examples/Tests/embedded_circle/analysis.py :language: python -This can also be included in an existing analysis script. Note that the plotfile must be ``_plt?????``, as is generated by the CI framework. +This can also be included as part of an existing analysis script. -Evaluate a checksum regression test from a bash terminal --------------------------------------------------------- +How to evaluate checksums from the command line +----------------------------------------------- You can execute ``checksumAPI.py`` as a Python script for that, and pass the plotfile that you want to evaluate, as well as the test name (so the script knows which benchmark to compare it to). @@ -41,11 +45,8 @@ See additional options * ``--rtol`` relative tolerance for the comparison * ``--atol`` absolute tolerance for the comparison (a sum of both is used by ``numpy.isclose()``) -Create/Reset a benchmark with new values that you know are correct ------------------------------------------------------------------- - -Create/Reset a benchmark from a plotfile generated locally -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +How to create or reset checksums with local benchmark values +------------------------------------------------------------ This is using ``checksumAPI.py`` as a Python script. @@ -65,8 +66,8 @@ Since this will automatically change the JSON file stored on the repo, make a se git add .json git commit -m "reset benchmark for because ..." --author="Tools " -Automated reset of a list of test benchmarks --------------------------------------------- +How to reset checksums for a list of tests with local benchmark values +---------------------------------------------------------------------- If you set the environment variable ``export CHECKSUM_RESET=ON`` before running tests that are compared against existing benchmarks, the test analysis will reset the benchmarks to the new values, skipping the comparison. @@ -80,8 +81,8 @@ With `CTest `__ (coming # ... check and commit changes ... -Reset a benchmark from the Azure pipeline output on Github -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +How to reset checksums for a list of tests with benchmark values from the Azure pipeline output +----------------------------------------------------------------------------------------------- Alternatively, the benchmarks can be reset using the output of the Azure continuous intergration (CI) tests on Github. The output can be accessed by following the steps below: diff --git a/Docs/source/developers/testing.rst b/Docs/source/developers/testing.rst index ee5c82aeea9..111e3e7d7cb 100644 --- a/Docs/source/developers/testing.rst +++ b/Docs/source/developers/testing.rst @@ -1,35 +1,45 @@ .. _developers-testing: -Testing the code +Testing the Code ================ -When adding a new feature, you want to make sure that (i) you did not break the existing code and (ii) your contribution gives correct results. While the code is tested regularly remotely (on the cloud when commits are pushed to an open PR, and every night on local clusters), it can also be useful to run tests on your custom input file. This section details how to use both automated and custom tests. +When proposing a code change, you want to make sure that -Continuous Integration in WarpX -------------------------------- +* the code change does not break the existing code; +* the code change gives correct results (numerics, physics, etc.). -Configuration -^^^^^^^^^^^^^ +WarpX follows the continuous integration (CI) software development practice, where automated builds and tests are run after merging code changes into the main branch. -Our regression tests are run with `CTest `__, an executable that comes with CMake. - -The test suite is ready to run once you have configured and built WarpX with CMake, following the instructions that you find in our :ref:`Users ` or :ref:`Developers ` sections. - -A test that requires a build option that was not configured and built will be skipped automatically. For example, if you configure and build WarpX in 1D only, any test of dimensionality other than 1D, which would require WarpX to be configured and built in the corresponding dimensionality, will be skipped automatically. +While the code is tested regularly remotely (on the cloud when commits are pushed to an open PR, and every night on local clusters), it can also be useful to run tests on your custom input file. How to run pre-commit tests locally ----------------------------------- -When proposing code changes to Warpx, we perform a couple of automated stylistic and correctness checks on the code change. -You can run those locally before you push to save some time, install them once like this: +First, when proposing a code change, we perform a couple of automated style and correctness checks. + +If you install the ``pre-commit`` tool on your local machine via .. code-block:: sh python -m pip install -U pre-commit pre-commit install +the style and correctness checks will run automatically on your local machine, after you commit the change and before you push. + +If you do not install the ``pre-commit`` tool on your local machine, these checks will run automatically as part of our CI workflows and a commit containing style and correctness changes might be added automatically to your branch. +In that case, you will need to pull that automated commit before pushing further changes. + See `pre-commit.com `__ and our ``.pre-commit-config.yaml`` file in the repository for more details. +How to configure the automated tests +------------------------------------ + +Our regression tests are run with `CTest `__, an executable that comes with CMake. + +The test suite is ready to run once you have configured and built WarpX with CMake, following the instructions that you find in our :ref:`Users ` or :ref:`Developers ` sections. + +A test that requires a build option that was not configured and built will be skipped automatically. For example, if you configure and build WarpX in 1D only, any test of dimensionality other than 1D, which would require WarpX to be configured and built in the corresponding dimensionality, will be skipped automatically. + How to run automated tests locally ---------------------------------- @@ -107,7 +117,15 @@ If you modify the code base locally and want to assess the effects of your code How to add automated tests -------------------------- -As mentioned above, the input files and scripts used by the automated tests can be found in the `Examples `__ directory, either under `Physics_applications `__ or `Tests `__. +An automated test typically consists of the following components: + +* input file or PICMI input script; +* analysis script; +* checksum file. + +To learn more about how to use checksums in automated tests, please see the corresponding section :ref:`Checksums on Tests `. + +As mentioned above, the input files and scripts used by the automated tests can be found in the `Examples `__ directory, under either `Physics_applications `__ or `Tests `__. Each test directory must contain a file named ``CMakeLists.txt`` where all tests associated with the input files and scripts in that directory must be listed. @@ -173,7 +191,8 @@ A new test can be added by adding a corresponding entry in ``CMakeLists.txt`` as If you need a new Python package dependency for testing, please add it in `Regression/requirements.txt `__. -Sometimes two or more tests share a large number of input parameters. The shared input parameters can be collected in a "base" input file that can be passed as a runtime parameter in the actual test input files through the parameter ``FILE``. +Sometimes two or more tests share a large number of input parameters. +The shared input parameters can be collected in a "base" input file that can be passed as a runtime parameter in the actual test input files through the parameter ``FILE``. If the new test is added in a new directory that did not exist before, please add the name of that directory with the command ``add_subdirectory`` in `Physics_applications/CMakeLists.txt `__ or `Tests/CMakeLists.txt `__, depending on where the new test directory is located. diff --git a/Docs/source/developers/workflows.rst b/Docs/source/developers/workflows.rst index 00279018e9d..f7c81ae70d8 100644 --- a/Docs/source/developers/workflows.rst +++ b/Docs/source/developers/workflows.rst @@ -8,7 +8,7 @@ Workflows profiling testing - documentation checksum - local_compile run_clang_tidy_locally + local_compile + documentation diff --git a/Docs/source/install/cmake.rst b/Docs/source/install/cmake.rst index 60d9eecc2b4..41e4c40bc85 100644 --- a/Docs/source/install/cmake.rst +++ b/Docs/source/install/cmake.rst @@ -77,9 +77,9 @@ For example, this builds WarpX in all geometries, enables Python bindings and Nv Build Options ------------- -============================= ============================================ ========================================================= +============================= ============================================ =========================================================== CMake Option Default & Values Description -============================= ============================================ ========================================================= +============================= ============================================ =========================================================== ``CMAKE_BUILD_TYPE`` RelWithDebInfo/**Release**/Debug `Type of build, symbols & optimizations `__ ``CMAKE_INSTALL_PREFIX`` system-dependent path `Install path prefix `__ ``CMAKE_VERBOSE_MAKEFILE`` ON/**OFF** `Print all compiler commands to the terminal during build `__ @@ -105,9 +105,9 @@ CMake Option Default & Values Descr ``WarpX_QED_TABLES_GEN_OMP`` **AUTO**/ON/OFF Enables OpenMP support for QED lookup tables generation ``WarpX_SENSEI`` ON/**OFF** SENSEI in situ visualization ``Python_EXECUTABLE`` (newest found) Path to Python executable -``PY_PIP_OPTIONS`` ``-v`` Additional options for ``pip``, e.g., ``-vvv`` -``PY_PIP_INSTALL_OPTIONS`` Additional options for ``pip install``, e.g., ``--user`` -============================= ============================================ ========================================================= +``PY_PIP_OPTIONS`` ``-v`` Additional options for ``pip``, e.g., ``-vvv;-q`` +``PY_PIP_INSTALL_OPTIONS`` Additional options for ``pip install``, e.g., ``--user;-q`` +============================= ============================================ =========================================================== WarpX can be configured in further detail with options from AMReX, which are documented in the AMReX manual: diff --git a/Docs/source/usage/faq.rst b/Docs/source/usage/faq.rst index 67cea8d6621..4ed0f8fa6af 100644 --- a/Docs/source/usage/faq.rst +++ b/Docs/source/usage/faq.rst @@ -74,10 +74,10 @@ Several BTD quantities differ slightly from the lab frame domain described in th In the following discussion, we will use a subscript input (e.g. :math:`\Delta z_{\rm input}`) to denote properties of the lab frame domain. -- The first back-transformed diagnostic (BTD) snapshot may not occur at :math:`t=0`. Rather, it occurs at :math:`t_0=\frac{z_{max}}c \beta(1+\beta)\gamma^2`. This is the first time when the boosted frame can complete the snapshot. +- The first back-transformed diagnostic (BTD) snapshot may not occur at :math:`t=0`. Rather, it occurs at :math:`t_0=\frac{z_{max}}c \beta/(1 - \beta \beta_{mw})`, where :math:`\beta_{mw}` represents the speed of the moving window. This is the first time when the boosted frame can complete the snapshot. - The grid spacing of the BTD snapshot is different from the grid spacing indicated in the input script. It is given by :math:`\Delta z_{\rm grid,snapshot}=\frac{c\Delta t_{\rm boost}}{\gamma\beta}`. For a CFL-limited time step, :math:`\Delta z_{\rm grid,snapshot}\approx \frac{1+\beta}{\beta} \Delta z_{\rm input}\approx 2 \Delta z_{\rm input}`. Hence in many common use cases at large boost, it is expected that the BTD snapshot has a grid spacing twice what is expressed in the input script. - The effective length of the BTD snapshot may be longer than anticipated from the input script because the grid spacing is different. Additionally, the number of grid points in the BTD snapshot is a multiple of ``.buffer_size`` whereas the number of grid cells specified in the input deck may not be. -- The code may require longer than anticipated to complete a BTD snapshot. The code starts filling the :math:`i^{th}` snapshot around step :math:`j_{\rm BTD start}={\rm ceil}\left( i\gamma(1-\beta)\frac{\Delta t_{\rm snapshot}}{\Delta t_{\rm boost}}\right)`. The code then saves information for one BTD cell every time step in the boosted frame simulation. The :math:`i^{th}` snapshot is completed and saved :math:`n_{z,{\rm snapshot}}=n_{\rm buffers}\cdot ({\rm buffer\ size})` time steps after it begins, which is when the effective snapshot length is covered by the simulation. +- The code may require longer than anticipated to complete a BTD snapshot. The code starts filling the :math:`i^{th}` snapshot around step :math:`j_{\rm BTD start}={\rm ceil}\left( i\gamma(1-\beta\beta_{mw})\frac{\Delta t_{\rm snapshot}}{\Delta t_{\rm boost}}\right)`. The code then saves information for one BTD cell every time step in the boosted frame simulation. The :math:`i^{th}` snapshot is completed and saved :math:`n_{z,{\rm snapshot}}=n_{\rm buffers}\cdot ({\rm buffer\ size})` time steps after it begins, which is when the effective snapshot length is covered by the simulation. What kinds of RZ output do you support? --------------------------------------- diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index b9d82d5014a..5014d421fb8 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -971,16 +971,21 @@ Particle initialization The ``external_file`` option is currently implemented for 2D, 3D and RZ geometries, with record components in the cartesian coordinates ``(x,y,z)`` for 3D and RZ, and ``(x,z)`` for 2D. For more information on the `openPMD format `__ and how to build WarpX with it, please visit :ref:`the install section `. - * ``NFluxPerCell``: Continuously inject a flux of macroparticles from a planar surface. + * ``NFluxPerCell``: Continuously inject a flux of macroparticles from a surface. The emitting surface can be chosen to be either a plane + defined by the user (using some of the parameters listed below), or the embedded boundary (see :ref:`Embedded Boundary Conditions `). This requires the additional parameters: * ``.flux_profile`` (see the description of this parameter further below) - * ``.surface_flux_pos`` (`double`, location of the injection plane [meter]) + * ``.inject_from_embedded_boundary`` (`0` or `1`, default `0` ; whether to inject from the embedded boundary or from a user-specified plane. + When injecting from the embedded boundary, the momentum distribution specified by the user along ``z`` (see e.g. ``uz_m``, ``uz_th`` below) is interpreted + as the momentum distribution along the local normal to the embedded boundary.) - * ``.flux_normal_axis`` (`x`, `y`, or `z` for 3D, `x` or `z` for 2D, or `r`, `t`, or `z` for RZ. When `flux_normal_axis` is `r` or `t`, the `x` and `y` components of the user-specified momentum distribution are interpreted as the `r` and `t` components respectively) + * ``.surface_flux_pos`` (only used when injecting from a plane, `double`, location of the injection plane [meter]) - * ``.flux_direction`` (`-1` or `+1`, direction of flux relative to the plane) + * ``.flux_normal_axis`` (only used when injecting from a plane, `x`, `y`, or `z` for 3D, `x` or `z` for 2D, or `r`, `t`, or `z` for RZ. When `flux_normal_axis` is `r` or `t`, the `x` and `y` components of the user-specified momentum distribution are interpreted as the `r` and `t` components respectively) + + * ``.flux_direction`` (only used when injecting from a plane, `-1` or `+1`, direction of flux relative to the plane) * ``.num_particles_per_cell`` (`double`) @@ -3466,14 +3471,15 @@ Reduced Diagnostics \frac{d\mathcal{L}}{d\mathcal{E}^*}(\mathcal{E}^*, t) = \int_0^t dt'\int d\boldsymbol{x}\,d\boldsymbol{p}_1 d\boldsymbol{p}_2\; \sqrt{ |\boldsymbol{v}_1 - \boldsymbol{v}_2|^2 - |\boldsymbol{v}_1\times\boldsymbol{v}_2|^2/c^2} \\ f_1(\boldsymbol{x}, \boldsymbol{p}_1, t')f_2(\boldsymbol{x}, \boldsymbol{p}_2, t') \delta(\mathcal{E}^* - \mathcal{E}^*(\boldsymbol{p}_1, \boldsymbol{p}_2)) - where :math:`\mathcal{E}^*(\boldsymbol{p}_1, \boldsymbol{p}_2) = \sqrt{m_1^2c^4 + m_2^2c^4 + 2(m_1 m_2 c^4 - \gamma_1 \gamma_2 - \boldsymbol{p}_1\cdot\boldsymbol{p}_2 c^2)}` is the energy in the center-of-mass frame, - and :math:`f_i` is the distribution function of species :math:`i`. Note that, if :math:`\sigma^*(\mathcal{E}^*)` + where :math:`f_i` is the distribution function of species :math:`i` and + :math:`\mathcal{E}^*(\boldsymbol{p}_1, \boldsymbol{p}_2) = \sqrt{m_1^2c^4 + m_2^2c^4 + 2 c^2{p_1}^\mu {p_2}_\mu}` + is the energy in the center-of-mass frame, where :math:`p^\mu = (\sqrt{m^2 c^2 + \boldsymbol{p}^2}, \boldsymbol{p})` + represents the 4-momentum. Note that, if :math:`\sigma^*(\mathcal{E}^*)` is the center-of-mass cross-section of a given collision process, then :math:`\int d\mathcal{E}^* \frac{d\mathcal{L}}{d\mathcal{E}^*} (\mathcal{E}^*, t)\sigma^*(\mathcal{E}^*)` gives the total number of collisions of that process (from the beginning of the simulation up until time :math:`t`). - The differential luminosity is given in units of :math:`\text{m}^{-2}.\text{eV}^{-1}`. For collider-relevant WarpX simulations + The differential luminosity is given in units of :math:`\text{m}^{-2}.\text{eV}^{-1}`. For collider-relevant WarpX simulations involving two crossing, high-energy beams of particles, the differential luminosity in :math:`\text{s}^{-1}.\text{m}^{-2}.\text{eV}^{-1}` can be obtained by multiplying the above differential luminosity by the expected repetition rate of the beams. diff --git a/Examples/Tests/diff_lumi_diag/CMakeLists.txt b/Examples/Tests/diff_lumi_diag/CMakeLists.txt index 1651d74115e..481847a023d 100644 --- a/Examples/Tests/diff_lumi_diag/CMakeLists.txt +++ b/Examples/Tests/diff_lumi_diag/CMakeLists.txt @@ -2,10 +2,20 @@ # add_warpx_test( - test_3d_diff_lumi_diag # name + test_3d_diff_lumi_diag_leptons # name 3 # dims 2 # nprocs - inputs_test_3d_diff_lumi_diag # inputs + inputs_test_3d_diff_lumi_diag_leptons # inputs + analysis.py # analysis + diags/diag1000080 # output + OFF # dependency +) + +add_warpx_test( + test_3d_diff_lumi_diag_photons # name + 3 # dims + 2 # nprocs + inputs_test_3d_diff_lumi_diag_photons # inputs analysis.py # analysis diags/diag1000080 # output OFF # dependency diff --git a/Examples/Tests/diff_lumi_diag/analysis.py b/Examples/Tests/diff_lumi_diag/analysis.py index 8f2061ff1dc..41501b1915d 100755 --- a/Examples/Tests/diff_lumi_diag/analysis.py +++ b/Examples/Tests/diff_lumi_diag/analysis.py @@ -37,16 +37,28 @@ * np.exp(-((E_bin - 2 * E_beam) ** 2) / (2 * sigma_E**2)) ) +# Extract test name from path +test_name = os.path.split(os.getcwd())[1] +print("test_name", test_name) + +# Pick tolerance +if "leptons" in test_name: + tol = 1e-2 +elif "photons" in test_name: + # In the photons case, the particles are + # initialized from a density distribution ; + # tolerance is larger due to lower particle statistics + tol = 6e-2 + # Check that the simulation result and analytical result match error = abs(dL_dE_sim - dL_dE_th).max() / abs(dL_dE_th).max() -tol = 1e-2 print("Relative error: ", error) print("Tolerance: ", tol) assert error < tol # compare checksums evaluate_checksum( - test_name=os.path.split(os.getcwd())[1], + test_name=test_name, output_file=sys.argv[1], rtol=1e-2, ) diff --git a/Examples/Tests/diff_lumi_diag/inputs_test_3d_diff_lumi_diag b/Examples/Tests/diff_lumi_diag/inputs_base_3d similarity index 81% rename from Examples/Tests/diff_lumi_diag/inputs_test_3d_diff_lumi_diag rename to Examples/Tests/diff_lumi_diag/inputs_base_3d index e8854937b6e..ba3c823b52b 100644 --- a/Examples/Tests/diff_lumi_diag/inputs_test_3d_diff_lumi_diag +++ b/Examples/Tests/diff_lumi_diag/inputs_base_3d @@ -6,12 +6,11 @@ my_constants.mc2_eV = m_e*clight*clight/q_e # BEAMS my_constants.beam_energy_eV = 125.e9 my_constants.beam_gamma = beam_energy_eV/(mc2_eV) -my_constants.beam_charge = 1.2e10*q_e +my_constants.beam_N = 1.2e10 my_constants.sigmax = 500e-9 my_constants.sigmay = 10e-9 my_constants.sigmaz = 300e-3 -my_constants.muz = -4*sigmaz -my_constants.nmacropart = 2e5 +my_constants.muz = 4*sigmaz # BOX my_constants.Lx = 8*sigmax @@ -62,17 +61,6 @@ warpx.poisson_solver = fft ################################# particles.species_names = beam1 beam2 -beam1.species_type = electron -beam1.injection_style = gaussian_beam -beam1.x_rms = sigmax -beam1.y_rms = sigmay -beam1.z_rms = sigmaz -beam1.x_m = 0 -beam1.y_m = 0 -beam1.z_m = muz -beam1.npart = nmacropart -beam1.q_tot = -beam_charge -beam1.z_cut = 4 beam1.momentum_distribution_type = gaussian beam1.uz_m = beam_gamma beam1.uy_m = 0.0 @@ -82,17 +70,6 @@ beam1.uy_th = 0 beam1.uz_th = 0.02*beam_gamma beam1.do_not_deposit = 1 -beam2.species_type = positron -beam2.injection_style = gaussian_beam -beam2.x_rms = sigmax -beam2.y_rms = sigmay -beam2.z_rms = sigmaz -beam2.x_m = 0 -beam2.y_m = 0 -beam2.z_m = -muz -beam2.npart = nmacropart -beam2.q_tot = beam_charge -beam2.z_cut = 4 beam2.momentum_distribution_type = gaussian beam2.uz_m = -beam_gamma beam2.uy_m = 0.0 @@ -108,7 +85,7 @@ beam2.do_not_deposit = 1 # FULL diagnostics.diags_names = diag1 -diag1.intervals = 1 +diag1.intervals = 80 diag1.diag_type = Full diag1.write_species = 1 diag1.fields_to_plot = rho_beam1 rho_beam2 diff --git a/Examples/Tests/diff_lumi_diag/inputs_test_3d_diff_lumi_diag_leptons b/Examples/Tests/diff_lumi_diag/inputs_test_3d_diff_lumi_diag_leptons new file mode 100644 index 00000000000..1cded30d3af --- /dev/null +++ b/Examples/Tests/diff_lumi_diag/inputs_test_3d_diff_lumi_diag_leptons @@ -0,0 +1,31 @@ +# base input parameters +FILE = inputs_base_3d + +# Test with electrons/positrons: use gaussian beam distribution +# by providing the total charge (q_tot) + +my_constants.nmacropart = 2e5 + +beam1.species_type = electron +beam1.injection_style = gaussian_beam +beam1.x_rms = sigmax +beam1.y_rms = sigmay +beam1.z_rms = sigmaz +beam1.x_m = 0 +beam1.y_m = 0 +beam1.z_m = -muz +beam1.npart = nmacropart +beam1.q_tot = -beam_N*q_e +beam1.z_cut = 4 + +beam2.species_type = positron +beam2.injection_style = gaussian_beam +beam2.x_rms = sigmax +beam2.y_rms = sigmay +beam2.z_rms = sigmaz +beam2.x_m = 0 +beam2.y_m = 0 +beam2.z_m = muz +beam2.npart = nmacropart +beam2.q_tot = beam_N*q_e +beam2.z_cut = 4 diff --git a/Examples/Tests/diff_lumi_diag/inputs_test_3d_diff_lumi_diag_photons b/Examples/Tests/diff_lumi_diag/inputs_test_3d_diff_lumi_diag_photons new file mode 100644 index 00000000000..f0ef254d911 --- /dev/null +++ b/Examples/Tests/diff_lumi_diag/inputs_test_3d_diff_lumi_diag_photons @@ -0,0 +1,28 @@ +# base input parameters +FILE = inputs_base_3d + +# Test with electrons/positrons: use parse_density_function + +beam1.species_type = electron +beam1.injection_style = "NUniformPerCell" +beam1.num_particles_per_cell_each_dim = 1 1 1 +beam1.profile = parse_density_function +beam1.density_function(x,y,z) = "beam_N/(sqrt(2*pi)*2*pi*sigmax*sigmay*sigmaz)*exp(-x*x/(2*sigmax*sigmax)-y*y/(2*sigmay*sigmay)-(z+muz)*(z+muz)/(2*sigmaz*sigmaz))" +beam1.xmin = -4*sigmax +beam1.xmax = 4*sigmax +beam1.ymin = -4*sigmay +beam1.ymax = 4*sigmay +beam1.zmin =-muz-4*sigmaz +beam1.zmax =-muz+4*sigmaz + +beam2.species_type = positron +beam2.injection_style = "NUniformPerCell" +beam2.num_particles_per_cell_each_dim = 1 1 1 +beam2.profile = parse_density_function +beam2.xmin = -4*sigmax +beam2.xmax = 4*sigmax +beam2.ymin = -4*sigmay +beam2.ymax = 4*sigmay +beam2.zmin = muz-4*sigmaz +beam2.zmax = muz+4*sigmaz +beam2.density_function(x,y,z) = "beam_N/(sqrt(2*pi)*2*pi*sigmax*sigmay*sigmaz)*exp(-x*x/(2*sigmax*sigmax)-y*y/(2*sigmay*sigmay)-(z-muz)*(z-muz)/(2*sigmaz*sigmaz))" diff --git a/Examples/Tests/flux_injection/CMakeLists.txt b/Examples/Tests/flux_injection/CMakeLists.txt index d09b83d7618..0929fc3d4c4 100644 --- a/Examples/Tests/flux_injection/CMakeLists.txt +++ b/Examples/Tests/flux_injection/CMakeLists.txt @@ -20,3 +20,33 @@ add_warpx_test( diags/diag1000120 # output OFF # dependency ) + +add_warpx_test( + test_3d_flux_injection_from_eb # name + 3 # dims + 2 # nprocs + inputs_test_3d_flux_injection_from_eb # inputs + analysis_flux_injection_from_eb.py # analysis + diags/diag1000010 # output + OFF # dependency +) + +add_warpx_test( + test_rz_flux_injection_from_eb # name + RZ # dims + 2 # nprocs + inputs_test_rz_flux_injection_from_eb # inputs + analysis_flux_injection_from_eb.py # analysis + diags/diag1000010 # output + OFF # dependency +) + +add_warpx_test( + test_2d_flux_injection_from_eb # name + 2 # dims + 2 # nprocs + inputs_test_2d_flux_injection_from_eb # inputs + analysis_flux_injection_from_eb.py # analysis + diags/diag1000010 # output + OFF # dependency +) diff --git a/Examples/Tests/flux_injection/analysis_flux_injection_from_eb.py b/Examples/Tests/flux_injection/analysis_flux_injection_from_eb.py new file mode 100755 index 00000000000..36ff50bea06 --- /dev/null +++ b/Examples/Tests/flux_injection/analysis_flux_injection_from_eb.py @@ -0,0 +1,161 @@ +#!/usr/bin/env python3 +# +# Copyright 2024 Remi Lehe +# +# This file is part of WarpX. +# +# License: BSD-3-Clause-LBNL + +""" +This script tests the emission of particles from the embedded boundary. +(In this case, the embedded boundary is a sphere in 3D and RZ, a cylinder in 2D.) +We check that the embedded boundary emits the correct number of particles, and that +the particle distributions are consistent with the expected distributions. +""" + +import os +import re +import sys + +import matplotlib.pyplot as plt +import numpy as np +import yt +from scipy.constants import c, m_e +from scipy.special import erf + +sys.path.insert(1, "../../../../warpx/Regression/Checksum/") +import checksumAPI + +yt.funcs.mylog.setLevel(0) + +# Open plotfile specified in command line +fn = sys.argv[1] +ds = yt.load(fn) +ad = ds.all_data() +t_max = ds.current_time.item() # time of simulation + +# Extract the dimensionality of the simulation +with open("./warpx_used_inputs", "r") as f: + warpx_used_inputs = f.read() +if re.search("geometry.dims = 2", warpx_used_inputs): + dims = "2D" +elif re.search("geometry.dims = RZ", warpx_used_inputs): + dims = "RZ" +elif re.search("geometry.dims = 3", warpx_used_inputs): + dims = "3D" + +# Total number of electrons expected: +# Simulation parameters determine the total number of particles emitted (Ntot) +flux = 1.0 # in m^-2.s^-1, from the input script +R = 2.0 # in m, radius of the sphere +if dims == "3D" or dims == "RZ": + emission_surface = 4 * np.pi * R**2 # in m^2 +elif dims == "2D": + emission_surface = 2 * np.pi * R # in m +Ntot = flux * emission_surface * t_max + +# Parameters of the histogram +hist_bins = 50 +hist_range = [-0.5, 0.5] + + +# Define function that histograms and checks the data +def gaussian_dist(u, u_th): + return 1.0 / ((2 * np.pi) ** 0.5 * u_th) * np.exp(-(u**2) / (2 * u_th**2)) + + +def gaussian_flux_dist(u, u_th, u_m): + normalization_factor = u_th**2 * np.exp(-(u_m**2) / (2 * u_th**2)) + ( + np.pi / 2 + ) ** 0.5 * u_m * u_th * (1 + erf(u_m / (2**0.5 * u_th))) + result = ( + 1.0 + / normalization_factor + * np.where(u > 0, u * np.exp(-((u - u_m) ** 2) / (2 * u_th**2)), 0) + ) + return result + + +def compare_gaussian(u, w, u_th, label=""): + du = (hist_range[1] - hist_range[0]) / hist_bins + w_hist, u_hist = np.histogram(u, bins=hist_bins, weights=w / du, range=hist_range) + u_hist = 0.5 * (u_hist[1:] + u_hist[:-1]) + w_th = Ntot * gaussian_dist(u_hist, u_th) + plt.plot(u_hist, w_hist, label=label + ": simulation") + plt.plot(u_hist, w_th, "--", label=label + ": theory") + assert np.allclose(w_hist, w_th, atol=0.07 * w_th.max()) + + +def compare_gaussian_flux(u, w, u_th, u_m, label=""): + du = (hist_range[1] - hist_range[0]) / hist_bins + w_hist, u_hist = np.histogram(u, bins=hist_bins, weights=w / du, range=hist_range) + u_hist = 0.5 * (u_hist[1:] + u_hist[:-1]) + w_th = Ntot * gaussian_flux_dist(u_hist, u_th, u_m) + plt.plot(u_hist, w_hist, label=label + ": simulation") + plt.plot(u_hist, w_th, "--", label=label + ": theory") + assert np.allclose(w_hist, w_th, atol=0.05 * w_th.max()) + + +# Load data and perform check + +plt.figure() + +if dims == "3D": + x = ad["electron", "particle_position_x"].to_ndarray() + y = ad["electron", "particle_position_y"].to_ndarray() + z = ad["electron", "particle_position_z"].to_ndarray() +elif dims == "2D": + x = ad["electron", "particle_position_x"].to_ndarray() + y = np.zeros_like(x) + z = ad["electron", "particle_position_y"].to_ndarray() +elif dims == "RZ": + theta = ad["electron", "particle_theta"].to_ndarray() + r = ad["electron", "particle_position_x"].to_ndarray() + x = r * np.cos(theta) + y = r * np.sin(theta) + z = ad["electron", "particle_position_y"].to_ndarray() +ux = ad["electron", "particle_momentum_x"].to_ndarray() / (m_e * c) +uy = ad["electron", "particle_momentum_y"].to_ndarray() / (m_e * c) +uz = ad["electron", "particle_momentum_z"].to_ndarray() / (m_e * c) +w = ad["electron", "particle_weight"].to_ndarray() + +# Check that the total number of particles emitted is correct +Ntot_sim = np.sum(w) +print("Ntot_sim = ", Ntot_sim) +print("Ntot = ", Ntot) +assert np.isclose(Ntot_sim, Ntot, rtol=0.01) + +# Check that none of the particles are inside the EB +# A factor 0.98 is applied to accomodate +# the cut-cell approximation of the sphere +assert np.all(x**2 + y**2 + z**2 > (0.98 * R) ** 2) + +# Check that the normal component of the velocity is consistent with the expected distribution +r = np.sqrt(x**2 + y**2 + z**2) +nx = x / r +ny = y / r +nz = z / r +u_n = ux * nx + uy * ny + uz * nz # normal component +compare_gaussian_flux(u_n, w, u_th=0.1, u_m=0.07, label="u_n") + +# Pick a direction that is orthogonal to the normal direction, and check the distribution +vx = ny / np.sqrt(nx**2 + ny**2) +vy = -nx / np.sqrt(nx**2 + ny**2) +vz = 0 +u_perp = ux * vx + uy * vy + uz * vz +compare_gaussian(u_perp, w, u_th=0.01, label="u_perp") + +# Pick the other perpendicular direction, and check the distribution +# The third direction is obtained by the cross product (n x v) +wx = ny * vz - nz * vy +wy = nz * vx - nx * vz +wz = nx * vy - ny * vx +u_perp2 = ux * wx + uy * wy + uz * wz +compare_gaussian(u_perp2, w, u_th=0.01, label="u_perp") + +plt.tight_layout() +plt.savefig("Distribution.png") + +# Verify checksum +test_name = os.path.split(os.getcwd())[1] +checksumAPI.evaluate_checksum(test_name, fn) diff --git a/Examples/Tests/flux_injection/inputs_base_from_eb b/Examples/Tests/flux_injection/inputs_base_from_eb new file mode 100644 index 00000000000..3e32d8799b6 --- /dev/null +++ b/Examples/Tests/flux_injection/inputs_base_from_eb @@ -0,0 +1,42 @@ +# Maximum number of time steps +max_step = 10 + +# The lo and hi ends of grids are multipliers of blocking factor +amr.blocking_factor = 8 + +# Maximum allowable size of each subdomain in the problem domain; +# this is used to decompose the domain for parallel calculations. +amr.max_grid_size = 8 + +# Maximum level in hierarchy (for now must be 0, i.e., one level in total) +amr.max_level = 0 + +# Deactivate Maxwell solver +algo.maxwell_solver = none +warpx.const_dt = 1e-9 + +# Embedded boundary +warpx.eb_implicit_function = "-(x**2+y**2+z**2-2**2)" + +# particles +particles.species_names = electron +algo.particle_shape = 3 + +electron.charge = -q_e +electron.mass = m_e +electron.injection_style = NFluxPerCell +electron.inject_from_embedded_boundary = 1 +electron.num_particles_per_cell = 100 +electron.flux_profile = parse_flux_function +electron.flux_function(x,y,z,t) = "1." +electron.momentum_distribution_type = gaussianflux +electron.ux_th = 0.01 +electron.uy_th = 0.01 +electron.uz_th = 0.1 +electron.uz_m = 0.07 + +# Diagnostics +diagnostics.diags_names = diag1 +diag1.intervals = 10 +diag1.diag_type = Full +diag1.fields_to_plot = none diff --git a/Examples/Tests/flux_injection/inputs_test_2d_flux_injection_from_eb b/Examples/Tests/flux_injection/inputs_test_2d_flux_injection_from_eb new file mode 100644 index 00000000000..f2e6f177887 --- /dev/null +++ b/Examples/Tests/flux_injection/inputs_test_2d_flux_injection_from_eb @@ -0,0 +1,13 @@ +FILE = inputs_base_from_eb + +# number of grid points +amr.n_cell = 16 16 + +# Geometry +geometry.dims = 2 +geometry.prob_lo = -4 -4 +geometry.prob_hi = 4 4 + +# Boundary condition +boundary.field_lo = periodic periodic +boundary.field_hi = periodic periodic diff --git a/Examples/Tests/flux_injection/inputs_test_3d_flux_injection_from_eb b/Examples/Tests/flux_injection/inputs_test_3d_flux_injection_from_eb new file mode 100644 index 00000000000..81ddc039977 --- /dev/null +++ b/Examples/Tests/flux_injection/inputs_test_3d_flux_injection_from_eb @@ -0,0 +1,13 @@ +FILE = inputs_base_from_eb + +# number of grid points +amr.n_cell = 16 16 16 + +# Geometry +geometry.dims = 3 +geometry.prob_lo = -4 -4 -4 +geometry.prob_hi = 4 4 4 + +# Boundary condition +boundary.field_lo = periodic periodic periodic +boundary.field_hi = periodic periodic periodic diff --git a/Examples/Tests/flux_injection/inputs_test_rz_flux_injection_from_eb b/Examples/Tests/flux_injection/inputs_test_rz_flux_injection_from_eb new file mode 100644 index 00000000000..4c970257f57 --- /dev/null +++ b/Examples/Tests/flux_injection/inputs_test_rz_flux_injection_from_eb @@ -0,0 +1,15 @@ +FILE = inputs_base_from_eb + +# number of grid points +amr.n_cell = 8 16 + +# Geometry +geometry.dims = RZ +geometry.prob_lo = 0 -4 +geometry.prob_hi = 4 4 + +# Boundary condition +boundary.field_lo = none periodic +boundary.field_hi = pec periodic + +electron.num_particles_per_cell = 300 diff --git a/GNUmakefile b/GNUmakefile index fe10983b780..1cc78403c7b 100644 --- a/GNUmakefile +++ b/GNUmakefile @@ -43,5 +43,8 @@ USE_RZ = FALSE USE_EB = FALSE +USE_LINEAR_SOLVERS_EM = TRUE +USE_LINEAR_SOLVERS_INCFLO = FALSE + WARPX_HOME := . include $(WARPX_HOME)/Source/Make.WarpX diff --git a/GOVERNANCE.rst b/GOVERNANCE.rst index b5253b80f9f..588e8b2df6e 100644 --- a/GOVERNANCE.rst +++ b/GOVERNANCE.rst @@ -16,7 +16,7 @@ Current Roster - Remi Lehe - Axel Huebl -See: `GitHub team `__ +See: `GitHub team `__ Role ^^^^ @@ -66,7 +66,7 @@ Current Roster - Weiqun Zhang - Edoardo Zoni -See: `GitHub team `__ +See: `GitHub team `__ Role ^^^^ diff --git a/Regression/Checksum/benchmarks_json/test_2d_flux_injection_from_eb.json b/Regression/Checksum/benchmarks_json/test_2d_flux_injection_from_eb.json new file mode 100644 index 00000000000..dd489f16e05 --- /dev/null +++ b/Regression/Checksum/benchmarks_json/test_2d_flux_injection_from_eb.json @@ -0,0 +1,11 @@ +{ + "lev=0": {}, + "electron": { + "particle_momentum_x": 6.990772711451971e-19, + "particle_momentum_y": 5.4131306169803364e-20, + "particle_momentum_z": 6.997294931789925e-19, + "particle_position_x": 35518.95120597846, + "particle_position_y": 35517.855675902414, + "particle_weight": 1.25355e-07 + } +} diff --git a/Regression/Checksum/benchmarks_json/test_2d_rigid_injection_btd.json b/Regression/Checksum/benchmarks_json/test_2d_rigid_injection_btd.json index 90cf134201f..9e876d5c23e 100644 --- a/Regression/Checksum/benchmarks_json/test_2d_rigid_injection_btd.json +++ b/Regression/Checksum/benchmarks_json/test_2d_rigid_injection_btd.json @@ -1,22 +1,22 @@ { + "lev=0": { + "Bx": 3.719030475087696e-05, + "By": 0.004843257051761486, + "Bz": 5.522765606391185e-06, + "Ex": 1461264.5033270014, + "Ey": 11205.64142004876, + "Ez": 282020.7784731542, + "jx": 16437877.898892798, + "jy": 2492340.3149980744, + "jz": 215102423.57036853, + "rho": 0.7246235591902171 + }, "beam": { - "particle_momentum_x": 2.2080215038948936e-16, + "particle_momentum_x": 2.2080215038948934e-16, "particle_momentum_y": 2.18711072170811e-16, - "particle_momentum_z": 2.730924530737497e-15, - "particle_position_x": 0.0260823588888081, + "particle_momentum_z": 2.730924530737456e-15, + "particle_position_x": 0.026082358888808558, "particle_position_y": 0.5049438607316916, "particle_weight": 62415.090744607645 - }, - "lev=0": { - "Bx": 3.721807007218884e-05, - "By": 0.004860056238272468, - "Bz": 5.5335765596325185e-06, - "Ex": 1466447.517373168, - "Ey": 11214.10223280318, - "Ez": 283216.0961218869, - "jx": 16437877.898892513, - "jy": 2492340.3149980404, - "jz": 215102423.57036877, - "rho": 0.7246235591902177 } } \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/test_3d_diff_lumi_diag.json b/Regression/Checksum/benchmarks_json/test_3d_diff_lumi_diag_leptons.json similarity index 100% rename from Regression/Checksum/benchmarks_json/test_3d_diff_lumi_diag.json rename to Regression/Checksum/benchmarks_json/test_3d_diff_lumi_diag_leptons.json diff --git a/Regression/Checksum/benchmarks_json/test_3d_diff_lumi_diag_photons.json b/Regression/Checksum/benchmarks_json/test_3d_diff_lumi_diag_photons.json new file mode 100644 index 00000000000..09b2031cdd2 --- /dev/null +++ b/Regression/Checksum/benchmarks_json/test_3d_diff_lumi_diag_photons.json @@ -0,0 +1,24 @@ +{ + "lev=0": { + "rho_beam1": 656097367.2335038, + "rho_beam2": 656097367.2335038 + }, + "beam1": { + "particle_momentum_x": 0.0, + "particle_momentum_y": 0.0, + "particle_momentum_z": 1.7512476113279403e-11, + "particle_position_x": 0.2621440000000001, + "particle_position_y": 0.005242880000000001, + "particle_position_z": 314572.79999473685, + "particle_weight": 11997744756.90957 + }, + "beam2": { + "particle_momentum_x": 0.0, + "particle_momentum_y": 0.0, + "particle_momentum_z": 1.7513431895752007e-11, + "particle_position_x": 0.2621440000000001, + "particle_position_y": 0.005242880000000001, + "particle_position_z": 314572.79999472946, + "particle_weight": 11997744756.909573 + } +} \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/test_3d_flux_injection_from_eb.json b/Regression/Checksum/benchmarks_json/test_3d_flux_injection_from_eb.json new file mode 100644 index 00000000000..e947a8af07b --- /dev/null +++ b/Regression/Checksum/benchmarks_json/test_3d_flux_injection_from_eb.json @@ -0,0 +1,12 @@ +{ + "lev=0": {}, + "electron": { + "particle_momentum_x": 4.371688233196277e-18, + "particle_momentum_y": 4.368885079657374e-18, + "particle_momentum_z": 4.367429424105371e-18, + "particle_position_x": 219746.94401890738, + "particle_position_y": 219690.7015248918, + "particle_position_z": 219689.45580938633, + "particle_weight": 4.954974999999999e-07 + } +} \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/test_rz_flux_injection_from_eb.json b/Regression/Checksum/benchmarks_json/test_rz_flux_injection_from_eb.json new file mode 100644 index 00000000000..23884de9725 --- /dev/null +++ b/Regression/Checksum/benchmarks_json/test_rz_flux_injection_from_eb.json @@ -0,0 +1,12 @@ +{ + "lev=0": {}, + "electron": { + "particle_momentum_x": 6.734984863106283e-19, + "particle_momentum_y": 6.786279785869023e-19, + "particle_momentum_z": 1.0527983828124758e-18, + "particle_position_x": 53309.270966506396, + "particle_position_y": 53302.3776094842, + "particle_theta": 58707.74469425615, + "particle_weight": 4.991396867417661e-07 + } +} \ No newline at end of file diff --git a/Source/Diagnostics/BTDiagnostics.H b/Source/Diagnostics/BTDiagnostics.H index d11db98276b..ab04f30ef18 100644 --- a/Source/Diagnostics/BTDiagnostics.H +++ b/Source/Diagnostics/BTDiagnostics.H @@ -161,6 +161,7 @@ private: * in z-direction for both 2D and 3D simulations in the Cartesian frame of reference. */ int m_moving_window_dir; + amrex::Real m_moving_window_beta; /** Number of back-transformed snapshots in the lab-frame requested by the user */ int m_num_snapshots_lab = std::numeric_limits::lowest(); diff --git a/Source/Diagnostics/BTDiagnostics.cpp b/Source/Diagnostics/BTDiagnostics.cpp index b07a9ddfffb..f9708056944 100644 --- a/Source/Diagnostics/BTDiagnostics.cpp +++ b/Source/Diagnostics/BTDiagnostics.cpp @@ -73,6 +73,7 @@ void BTDiagnostics::DerivedInitData () m_gamma_boost = WarpX::gamma_boost; m_beta_boost = std::sqrt( 1._rt - 1._rt/( m_gamma_boost * m_gamma_boost) ); m_moving_window_dir = WarpX::moving_window_dir; + m_moving_window_beta = WarpX::moving_window_v/PhysConst::c; // Currently, for BTD, all the data is averaged+coarsened to coarsest level // and then sliced+back-transformed+filled_to_buffer. // The number of levels to be output is nlev_output. @@ -142,7 +143,7 @@ void BTDiagnostics::DerivedInitData () const int lev = 0; const amrex::Real dt_boosted_frame = warpx.getdt(lev); const int moving_dir = WarpX::moving_window_dir; - const amrex::Real Lz_lab = warpx.Geom(lev).ProbLength(moving_dir) / WarpX::gamma_boost / (1._rt+WarpX::beta_boost); + const amrex::Real Lz_lab = warpx.Geom(lev).ProbLength(moving_dir) * WarpX::gamma_boost * (1._rt - WarpX::beta_boost*m_moving_window_beta); const int ref_ratio = 1; const amrex::Real dz_snapshot_grid = dz_lab(dt_boosted_frame, ref_ratio); // Need enough buffers so the snapshot length is longer than the lab frame length @@ -153,22 +154,21 @@ void BTDiagnostics::DerivedInitData () // the final snapshot starts filling when the // right edge of the moving window intersects the final snapshot // time of final snapshot : t_sn = t0 + i*dt_snapshot - // where t0 is the time of first BTD snapshot, t0 = zmax / c * beta / (1-beta) + // where t0 is the time of first BTD snapshot, t0 = zmax / c * beta / (1-beta*beta_mw) // // the right edge of the moving window at the time of the final snapshot // has space time coordinates - // time t_intersect = t_sn, position z_intersect=zmax + c*t_sn + // time t_intersect = t_sn, position z_intersect=zmax + v_mw*t_sn // the boosted time of this space time pair is // t_intersect_boost = gamma * (t_intersect - beta * z_intersect_boost/c) - // = gamma * (t_sn * (1 - beta) - beta * zmax / c) - // = gamma * (zmax*beta/c + i*dt_snapshot*(1-beta) - beta*zmax/c) - // = gamma * i * dt_snapshot * (1-beta) - // = i * dt_snapshot / gamma / (1+beta) + // = gamma * (t_sn * (1 - beta*beta_mw) - beta * zmax / c) + // = gamma * (zmax*beta/c + i*dt_snapshot*(1-beta*beta_mw) - beta*zmax/c) + // = gamma * (1-beta*beta_mw) * i * dt_snapshot // // if j = final snapshot starting step, then we want to solve - // j dt_boosted_frame >= t_intersect_boost = i * dt_snapshot / gamma / (1+beta) - // j >= i / gamma / (1+beta) * dt_snapshot / dt_boosted_frame - const int final_snapshot_starting_step = static_cast(std::ceil(final_snapshot_iteration / WarpX::gamma_boost / (1._rt+WarpX::beta_boost) * m_dt_snapshots_lab / dt_boosted_frame)); + // j dt_boosted_frame >= t_intersect_boost = i * gamma * (1-beta*beta_mw) * dt_snapshot + // j >= i * gamma * (1-beta*beta_mw) * dt_snapshot / dt_boosted_frame + const int final_snapshot_starting_step = static_cast(std::ceil(final_snapshot_iteration * WarpX::gamma_boost * (1._rt - WarpX::beta_boost*m_moving_window_beta) * m_dt_snapshots_lab / dt_boosted_frame)); const int final_snapshot_fill_iteration = final_snapshot_starting_step + num_buffers * m_buffer_size - 1; const amrex::Real final_snapshot_fill_time = final_snapshot_fill_iteration * dt_boosted_frame; if (WarpX::compute_max_step_from_btd) { @@ -262,7 +262,7 @@ BTDiagnostics::ReadParameters () bool snapshot_interval_is_specified = utils::parser::queryWithParser( pp_diag_name, "dt_snapshots_lab", m_dt_snapshots_lab); if ( utils::parser::queryWithParser(pp_diag_name, "dz_snapshots_lab", m_dz_snapshots_lab) ) { - m_dt_snapshots_lab = m_dz_snapshots_lab/PhysConst::c; + m_dt_snapshots_lab = m_dz_snapshots_lab/WarpX::moving_window_v; snapshot_interval_is_specified = true; } WARPX_ALWAYS_ASSERT_WITH_MESSAGE(snapshot_interval_is_specified, @@ -344,13 +344,15 @@ BTDiagnostics::InitializeBufferData ( int i_buffer , int lev, bool restart) // When restarting boosted simulations, the code below needs to take // into account the fact that the position of the box at the beginning // of the simulation, is not the one that we had at t=0 (because of the moving window) - const amrex::Real boosted_moving_window_v = (WarpX::moving_window_v - m_beta_boost*PhysConst::c) - / (1._rt - m_beta_boost * WarpX::moving_window_v/PhysConst::c); + const amrex::Real boosted_moving_window_v = (m_moving_window_beta - m_beta_boost) + / (1._rt - m_beta_boost*m_moving_window_beta); // Lab-frame time for the i^th snapshot if (!restart) { - const amrex::Real zmax_0 = warpx.Geom(lev).ProbHi(m_moving_window_dir); + const amrex::Real zmax_boost = warpx.Geom(lev).ProbHi(m_moving_window_dir); 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; + + m_gamma_boost*m_beta_boost*zmax_boost/PhysConst::c; + // Note: gamma_boost*beta_boost*zmax_boost is equal to + // beta_boost*zmax_lab/(1-beta_boost*beta_moving_window) } // Define buffer domain in boosted frame at level, lev, with user-defined lo and hi @@ -409,9 +411,9 @@ BTDiagnostics::InitializeBufferData ( int i_buffer , int lev, bool restart) // Define buffer_domain in lab-frame for the i^th snapshot. // Replace z-dimension with lab-frame co-ordinates. const amrex::Real zmin_buffer_lab = ( diag_dom.lo(m_moving_window_dir) - boosted_moving_window_v * warpx.gett_new(0) ) - / ( (1.0_rt + m_beta_boost) * m_gamma_boost); + * (1.0_rt - m_beta_boost*m_moving_window_beta) * m_gamma_boost; const amrex::Real zmax_buffer_lab = ( diag_dom.hi(m_moving_window_dir) - boosted_moving_window_v * warpx.gett_new(0) ) - / ( (1.0_rt + m_beta_boost) * m_gamma_boost); + * (1.0_rt - m_beta_boost*m_moving_window_beta) * m_gamma_boost; // Initialize buffer counter and z-positions of the i^th snapshot in // boosted-frame and lab-frame @@ -1010,12 +1012,15 @@ BTDiagnostics::GetZSliceInDomainFlag (const int i_buffer, const int lev) { auto & warpx = WarpX::GetInstance(); const amrex::RealBox& boost_domain = warpx.Geom(lev).ProbDomain(); + const amrex::Real boost_cellsize = warpx.Geom(lev).CellSize(m_moving_window_dir); const amrex::Real buffer_zmin_lab = m_snapshot_domain_lab[i_buffer].lo( m_moving_window_dir ); const amrex::Real buffer_zmax_lab = m_snapshot_domain_lab[i_buffer].hi( m_moving_window_dir ); + // Exclude 0.5*boost_cellsize from the edge, to avoid that the interpolation to + // cell centers uses data from the guard cells. const bool slice_not_in_domain = - ( m_current_z_boost[i_buffer] <= boost_domain.lo(m_moving_window_dir) ) || - ( m_current_z_boost[i_buffer] >= boost_domain.hi(m_moving_window_dir) ) || + ( m_current_z_boost[i_buffer] <= boost_domain.lo(m_moving_window_dir) + 0.5_rt*boost_cellsize) || + ( m_current_z_boost[i_buffer] >= boost_domain.hi(m_moving_window_dir) - 0.5_rt*boost_cellsize) || ( m_current_z_lab[i_buffer] <= buffer_zmin_lab ) || ( m_current_z_lab[i_buffer] >= buffer_zmax_lab ); diff --git a/Source/Diagnostics/ReducedDiags/DifferentialLuminosity.cpp b/Source/Diagnostics/ReducedDiags/DifferentialLuminosity.cpp index 7f31f1b56ed..8208458b9d7 100644 --- a/Source/Diagnostics/ReducedDiags/DifferentialLuminosity.cpp +++ b/Source/Diagnostics/ReducedDiags/DifferentialLuminosity.cpp @@ -132,9 +132,8 @@ void DifferentialLuminosity::ComputeDiags (int step) // Since this diagnostic *accumulates* the luminosity in the // array d_data, we add contributions at *each timestep*, but // we only write the data to file at intervals specified by the user. - - const Real c2_over_qe = PhysConst::c*PhysConst::c/PhysConst::q_e; - const Real inv_c2 = 1._rt/(PhysConst::c*PhysConst::c); + const Real c_sq = PhysConst::c*PhysConst::c; + const Real c_over_qe = PhysConst::c/PhysConst::q_e; // get a reference to WarpX instance auto& warpx = WarpX::GetInstance(); @@ -187,6 +186,7 @@ void DifferentialLuminosity::ComputeDiags (int step) amrex::ParticleReal * const AMREX_RESTRICT u1x = soa_1.m_rdata[PIdx::ux]; amrex::ParticleReal * const AMREX_RESTRICT u1y = soa_1.m_rdata[PIdx::uy]; // v*gamma=p/m amrex::ParticleReal * const AMREX_RESTRICT u1z = soa_1.m_rdata[PIdx::uz]; + bool const species1_is_photon = species_1.AmIA(); const auto soa_2 = ptile_2.getParticleTileData(); index_type* AMREX_RESTRICT indices_2 = bins_2.permutationPtr(); @@ -196,6 +196,7 @@ void DifferentialLuminosity::ComputeDiags (int step) amrex::ParticleReal * const AMREX_RESTRICT u2x = soa_2.m_rdata[PIdx::ux]; amrex::ParticleReal * const AMREX_RESTRICT u2y = soa_2.m_rdata[PIdx::uy]; amrex::ParticleReal * const AMREX_RESTRICT u2z = soa_2.m_rdata[PIdx::uz]; + bool const species2_is_photon = species_2.AmIA(); // Extract low-level data auto const n_cells = static_cast(bins_1.numBins()); @@ -218,34 +219,59 @@ void DifferentialLuminosity::ComputeDiags (int step) index_type const j_1 = indices_1[i_1]; index_type const j_2 = indices_2[i_2]; - Real const u1_square = u1x[j_1]*u1x[j_1] + u1y[j_1]*u1y[j_1] + u1z[j_1]*u1z[j_1]; - Real const gamma1 = std::sqrt(1._rt + u1_square*inv_c2); - Real const u2_square = u2x[j_2]*u2x[j_2] + u2y[j_2]*u2y[j_2] + u2z[j_2]*u2z[j_2]; - Real const gamma2 = std::sqrt(1._rt + u2_square*inv_c2); - Real const u1_dot_u2 = u1x[j_1]*u2x[j_2] + u1y[j_1]*u2y[j_2] + u1z[j_1]*u2z[j_2]; + Real p1t=0, p1x=0, p1y=0, p1z=0; // components of 4-momentum of particle 1 + Real const u1_sq = u1x[j_1]*u1x[j_1] + u1y[j_1]*u1y[j_1] + u1z[j_1]*u1z[j_1]; + if (species1_is_photon) { + // photon case (momentum is normalized by m_e in WarpX) + p1t = PhysConst::m_e*std::sqrt( u1_sq ); + p1x = PhysConst::m_e*u1x[j_1]; + p1y = PhysConst::m_e*u1y[j_1]; + p1z = PhysConst::m_e*u1z[j_1]; + } else { + p1t = m1*std::sqrt( c_sq + u1_sq ); + p1x = m1*u1x[j_1]; + p1y = m1*u1y[j_1]; + p1z = m1*u1z[j_1]; + } + + Real p2t=0, p2x=0, p2y=0, p2z=0; // components of 4-momentum of particle 2 + Real const u2_sq = u2x[j_2]*u2x[j_2] + u2y[j_2]*u2y[j_2] + u2z[j_2]*u2z[j_2]; + if (species2_is_photon) { + // photon case (momentum is normalized by m_e in WarpX) + p2t = PhysConst::m_e*std::sqrt(u2_sq); + p2x = PhysConst::m_e*u2x[j_2]; + p2y = PhysConst::m_e*u2y[j_2]; + p2z = PhysConst::m_e*u2z[j_2]; + } else { + p2t = m2*std::sqrt( c_sq + u2_sq ); + p2x = m2*u2x[j_2]; + p2y = m2*u2y[j_2]; + p2z = m2*u2z[j_2]; + } // center of mass energy in eV - Real const E_com = c2_over_qe * std::sqrt(m1*m1 + m2*m2 + 2*m1*m2* (gamma1*gamma2 - u1_dot_u2*inv_c2)); + Real const E_com = c_over_qe * std::sqrt(m1*m1*c_sq + m2*m2*c_sq + 2*(p1t*p2t - p1x*p2x - p1y*p2y - p1z*p2z)); // determine particle bin int const bin = int(Math::floor((E_com-bin_min)/bin_size)); if ( bin<0 || bin>=num_bins ) { continue; } // discard if out-of-range - Real const v1_minus_v2_x = u1x[j_1]/gamma1 - u2x[j_2]/gamma2; - Real const v1_minus_v2_y = u1y[j_1]/gamma1 - u2y[j_2]/gamma2; - Real const v1_minus_v2_z = u1z[j_1]/gamma1 - u2z[j_2]/gamma2; - Real const v1_minus_v2_square = v1_minus_v2_x*v1_minus_v2_x + v1_minus_v2_y*v1_minus_v2_y + v1_minus_v2_z*v1_minus_v2_z; + Real const inv_p1t = 1.0_rt/p1t; + Real const inv_p2t = 1.0_rt/p2t; - Real const u1_cross_u2_x = u1y[j_1]*u2z[j_2] - u1z[j_1]*u2y[j_2]; - Real const u1_cross_u2_y = u1z[j_1]*u2x[j_2] - u1x[j_1]*u2z[j_2]; - Real const u1_cross_u2_z = u1x[j_1]*u2y[j_2] - u1y[j_1]*u2x[j_2]; + Real const beta1_sq = (p1x*p1x + p1y*p1y + p1z*p1z) * inv_p1t*inv_p1t; + Real const beta2_sq = (p2x*p2x + p2y*p2y + p2z*p2z) * inv_p2t*inv_p2t; + Real const beta1_dot_beta2 = (p1x*p2x + p1y*p2y + p1z*p2z) * inv_p1t*inv_p2t; - Real const v1_cross_v2_square = (u1_cross_u2_x*u1_cross_u2_x + u1_cross_u2_y*u1_cross_u2_y + u1_cross_u2_z*u1_cross_u2_z) / (gamma1*gamma1*gamma2*gamma2); + // Here we use the fact that: + // (v1 - v2)^2 = v1^2 + v2^2 - 2 v1.v2 + // and (v1 x v2)^2 = v1^2 v2^2 - (v1.v2)^2 + // we also use beta=v/c instead of v - Real const radicand = v1_minus_v2_square - v1_cross_v2_square * inv_c2; + Real const radicand = beta1_sq + beta2_sq - 2*beta1_dot_beta2 - beta1_sq*beta2_sq + beta1_dot_beta2*beta1_dot_beta2; - Real const dL_dEcom = std::sqrt( radicand ) * w1[j_1] * w2[j_2] / dV / bin_size * dt; // m^-2 eV^-1 + Real const dL_dEcom = PhysConst::c * std::sqrt( radicand ) * w1[j_1] * w2[j_2] / dV / bin_size * dt; // m^-2 eV^-1 amrex::HostDevice::Atomic::Add(&dptr_data[bin], dL_dEcom); diff --git a/Source/EmbeddedBoundary/DistanceToEB.H b/Source/EmbeddedBoundary/DistanceToEB.H index 0c13724380c..0b27fd054cd 100644 --- a/Source/EmbeddedBoundary/DistanceToEB.H +++ b/Source/EmbeddedBoundary/DistanceToEB.H @@ -121,7 +121,13 @@ amrex::RealVect interp_normal (int i, int j, int k, const amrex::Real W[AMREX_SP #else amrex::ignore_unused(i, j, k, ic, jc, kc, W, Wc, phi, dxi); amrex::RealVect normal(0.0); - WARPX_ABORT_WITH_MESSAGE("Error: interp_distance not yet implemented in 1D"); + + AMREX_IF_ON_DEVICE(( + AMREX_DEVICE_ASSERT(0); + )) + AMREX_IF_ON_HOST(( + WARPX_ABORT_WITH_MESSAGE("Error: interp_normal not yet implemented in 1D"); + )) #endif return normal; diff --git a/Source/Initialization/PlasmaInjector.H b/Source/Initialization/PlasmaInjector.H index b9fe2323290..f14720d271c 100644 --- a/Source/Initialization/PlasmaInjector.H +++ b/Source/Initialization/PlasmaInjector.H @@ -131,6 +131,8 @@ public: int flux_normal_axis; int flux_direction; // -1 for left, +1 for right + bool m_inject_from_eb = false; // whether to inject from the embedded boundary + bool radially_weighted = true; std::string str_flux_function; diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp index b5500dace6d..f2db5f4974d 100644 --- a/Source/Initialization/PlasmaInjector.cpp +++ b/Source/Initialization/PlasmaInjector.cpp @@ -9,6 +9,7 @@ */ #include "PlasmaInjector.H" +#include "EmbeddedBoundary/Enabled.H" #include "Initialization/GetTemperature.H" #include "Initialization/GetVelocity.H" #include "Initialization/InjectorDensity.H" @@ -307,57 +308,72 @@ void PlasmaInjector::setupNFluxPerCell (amrex::ParmParse const& pp_species) "(Please visit PR#765 for more information.)"); } #endif - utils::parser::getWithParser(pp_species, source_name, "surface_flux_pos", surface_flux_pos); - utils::parser::queryWithParser(pp_species, source_name, "flux_tmin", flux_tmin); - utils::parser::queryWithParser(pp_species, source_name, "flux_tmax", flux_tmax); - std::string flux_normal_axis_string; - utils::parser::get(pp_species, source_name, "flux_normal_axis", flux_normal_axis_string); - flux_normal_axis = -1; + + // Check whether injection from the embedded boundary is requested + utils::parser::queryWithParser(pp_species, source_name, "inject_from_embedded_boundary", m_inject_from_eb); + if (m_inject_from_eb) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( EB::enabled(), + "Error: Embedded boundary injection is only available when " + "embedded boundaries are enabled."); + flux_normal_axis = 2; // Interpret z as the normal direction to the EB + flux_direction = 1; + } else { + // Injection is through a plane in this case. + // Parse the parameters of the plane (position, normal direction, etc.) + + utils::parser::getWithParser(pp_species, source_name, "surface_flux_pos", surface_flux_pos); + utils::parser::queryWithParser(pp_species, source_name, "flux_tmin", flux_tmin); + utils::parser::queryWithParser(pp_species, source_name, "flux_tmax", flux_tmax); + std::string flux_normal_axis_string; + utils::parser::get(pp_species, source_name, "flux_normal_axis", flux_normal_axis_string); + flux_normal_axis = -1; #if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE) - if (flux_normal_axis_string == "r" || flux_normal_axis_string == "R") { - flux_normal_axis = 0; - } - if (flux_normal_axis_string == "t" || flux_normal_axis_string == "T") { - flux_normal_axis = 1; - } + if (flux_normal_axis_string == "r" || flux_normal_axis_string == "R") { + flux_normal_axis = 0; + } + if (flux_normal_axis_string == "t" || flux_normal_axis_string == "T") { + flux_normal_axis = 1; + } #else # ifndef WARPX_DIM_1D_Z - if (flux_normal_axis_string == "x" || flux_normal_axis_string == "X") { - flux_normal_axis = 0; - } + if (flux_normal_axis_string == "x" || flux_normal_axis_string == "X") { + flux_normal_axis = 0; + } # endif #endif #ifdef WARPX_DIM_3D - if (flux_normal_axis_string == "y" || flux_normal_axis_string == "Y") { - flux_normal_axis = 1; - } + if (flux_normal_axis_string == "y" || flux_normal_axis_string == "Y") { + flux_normal_axis = 1; + } #endif #if defined(WARPX_DIM_RSPHERE) - // "p" for phi - if (flux_normal_axis_string == "p" || flux_normal_axis_string == "P") { - flux_normal_axis = 2; - } + // "p" for phi + if (flux_normal_axis_string == "p" || flux_normal_axis_string == "P") { + flux_normal_axis = 2; + } #else - if (flux_normal_axis_string == "z" || flux_normal_axis_string == "Z") { - flux_normal_axis = 2; - } + if (flux_normal_axis_string == "z" || flux_normal_axis_string == "Z") { + flux_normal_axis = 2; + } #endif -#ifdef WARPX_DIM_3D - const std::string flux_normal_axis_help = "'x', 'y', or 'z'."; +#if defined(WARPX_DIM_3D) + const std::string flux_normal_axis_help = "'x', 'y', or 'z'."; #elif defined(WARPX_DIM_RZ) - const std::string flux_normal_axis_help = "'r' or 'z'."; + const std::string flux_normal_axis_help = "'r' or 'z'."; #elif defined(WARPX_DIM_XZ) - const std::string flux_normal_axis_help = "'x' or 'z'."; + const std::string flux_normal_axis_help = "'x' or 'z'."; #elif defined(WARPX_DIM_1D_Z) - const std::string flux_normal_axis_help = "'z'."; + const std::string flux_normal_axis_help = "'z'."; #elif defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE) - const std::string flux_normal_axis_help = "'r'."; + const std::string flux_normal_axis_help = "'r'."; #endif - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(flux_normal_axis >= 0, - "Error: Invalid value for flux_normal_axis. It must be " + flux_normal_axis_help); - utils::parser::getWithParser(pp_species, source_name, "flux_direction", flux_direction); - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(flux_direction == +1 || flux_direction == -1, - "Error: flux_direction must be -1 or +1."); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(flux_normal_axis >= 0, + "Error: Invalid value for flux_normal_axis. It must be " + flux_normal_axis_help); + utils::parser::getWithParser(pp_species, source_name, "flux_direction", flux_direction); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(flux_direction == +1 || flux_direction == -1, + "Error: flux_direction must be -1 or +1."); + } + // Construct InjectorPosition with InjectorPositionRandom. h_flux_pos = std::make_unique( (InjectorPositionRandomPlane*)nullptr, diff --git a/Source/Particles/AddPlasmaUtilities.H b/Source/Particles/AddPlasmaUtilities.H index e188fb671e2..1081bb69004 100644 --- a/Source/Particles/AddPlasmaUtilities.H +++ b/Source/Particles/AddPlasmaUtilities.H @@ -22,6 +22,29 @@ #include #include +struct PDim3 { + amrex::ParticleReal x, y, z; + + AMREX_GPU_HOST_DEVICE + PDim3(const amrex::XDim3& a): + x{static_cast(a.x)}, + y{static_cast(a.y)}, + z{static_cast(a.z)} + {} + + AMREX_GPU_HOST_DEVICE + ~PDim3() = default; + + AMREX_GPU_HOST_DEVICE + PDim3(PDim3 const &) = default; + AMREX_GPU_HOST_DEVICE + PDim3& operator=(PDim3 const &) = default; + AMREX_GPU_HOST_DEVICE + PDim3(PDim3&&) = default; + AMREX_GPU_HOST_DEVICE + PDim3& operator=(PDim3&&) = default; +}; + /* Finds the overlap region between the given tile_realbox and part_realbox, returning true if an overlap exists and false if otherwise. This also sets the parameters overlap_realbox, @@ -73,12 +96,124 @@ int compute_area_weights (const amrex::IntVect& iv, const int normal_axis) { return r; } + +#ifdef AMREX_USE_EB +/* + * \brief This computes the scale_fac (used for setting the particle weights) on a on area basis + * (used for flux injection from the embedded boundary). + * + * \param[in] dx: cell size in each direction + * \param[in] num_ppc_real: number of particles per cell + * \param[in] eb_bnd_normal_arr: array containing the normal to the embedded boundary + * \param[in] i, j, k: indices of the cell + * + * \return scale_fac: the scaling factor to be applied to the weight of the particles + */ +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +amrex::Real compute_scale_fac_area_eb ( + const amrex::GpuArray& dx, + const amrex::Real num_ppc_real, + amrex::Array4 const& eb_bnd_normal_arr, + int i, int j, int k ) { + using namespace amrex::literals; + // Scale particle weight by the area of the emitting surface, within one cell + // By definition, eb_bnd_area_arr is normalized (unitless). + // Here we undo the normalization (i.e. multiply by the surface used for normalization in amrex: + // see https://amrex-codes.github.io/amrex/docs_html/EB.html#embedded-boundary-data-structures) +#if defined(WARPX_DIM_3D) + const amrex::Real nx = eb_bnd_normal_arr(i,j,k,0); + const amrex::Real ny = eb_bnd_normal_arr(i,j,k,1); + const amrex::Real nz = eb_bnd_normal_arr(i,j,k,2); + amrex::Real scale_fac = std::sqrt(amrex::Math::powi<2>(nx*dx[1]*dx[2]) + + amrex::Math::powi<2>(ny*dx[0]*dx[2]) + + amrex::Math::powi<2>(nz*dx[0]*dx[1])); + +#elif defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ) + const amrex::Real nx = eb_bnd_normal_arr(i,j,k,0); + const amrex::Real nz = eb_bnd_normal_arr(i,j,k,1); + amrex::Real scale_fac = std::sqrt(amrex::Math::powi<2>(nx*dx[1]) + + amrex::Math::powi<2>(nz*dx[0])); +#else + amrex::ignore_unused(dx, eb_bnd_normal_arr, i, j, k); + amrex::Real scale_fac = 0.0_rt; +#endif + // Do not multiply by eb_bnd_area_arr(i,j,k) here because this + // already taken into account by emitting a number of macroparticles + // that is proportional to the area of eb_bnd_area_arr(i,j,k). + scale_fac /= num_ppc_real; + return scale_fac; +} + +/* \brief Rotate the momentum of the particle so that the z direction + * transforms to the normal of the embedded boundary. + * + * More specifically, before calling this function, `pu.z` has the + * momentum distribution that is meant for the direction normal to the + * embedded boundary, and `pu.x`/`pu.y` have the momentum distribution that + * is meant for the tangentional direction. After calling this function, + * `pu.x`, `pu.y`, `pu.z` will have the correct momentum distribution, + * consistent with the local normal to the embedded boundary. + * + * \param[inout] pu momentum of the particle + * \param[in] eb_bnd_normal_arr: array containing the normal to the embedded boundary + * \param[in] i, j, k: indices of the cell + * */ +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void rotate_momentum_eb ( + PDim3 & pu, + amrex::Array4 const& eb_bnd_normal_arr, + int i, int j, int k ) +{ + using namespace amrex::literals; + +#if defined(WARPX_DIM_3D) + // The minus sign below takes into account the fact that eb_bnd_normal_arr + // points towards the covered region, while particles are to be emitted + // *away* from the covered region. + amrex::Real const nx = -eb_bnd_normal_arr(i,j,k,0); + amrex::Real const ny = -eb_bnd_normal_arr(i,j,k,1); + amrex::Real const nz = -eb_bnd_normal_arr(i,j,k,2); + + // Rotate the momentum in theta and phi + amrex::Real const cos_theta = nz; + amrex::Real const sin_theta = std::sqrt(1._rt-nz*nz); + amrex::Real const nperp = std::sqrt(nx*nx + ny*ny); + amrex::Real cos_phi = 1; + amrex::Real sin_phi = 0; + if ( nperp > 0.0 ) { + cos_phi = nx/nperp; + sin_phi = ny/nperp; + } + // Apply rotation matrix + amrex::Real const ux = pu.x*cos_theta*cos_phi - pu.y*sin_phi + pu.z*sin_theta*cos_phi; + amrex::Real const uy = pu.x*cos_theta*sin_phi + pu.y*cos_phi + pu.z*sin_theta*sin_phi; + amrex::Real const uz = -pu.x*sin_theta + pu.z*cos_theta; + pu.x = ux; + pu.y = uy; + pu.z = uz; + +#elif defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ) + // The minus sign below takes into account the fact that eb_bnd_normal_arr + // points towards the covered region, while particles are to be emitted + // *away* from the covered region. + amrex::Real const sin_theta = -eb_bnd_normal_arr(i,j,k,0); + amrex::Real const cos_theta = -eb_bnd_normal_arr(i,j,k,1); + amrex::Real const uz = pu.z*cos_theta - pu.x*sin_theta; + amrex::Real const ux = pu.x*cos_theta + pu.z*sin_theta; + pu.x = ux; + pu.z = uz; +#else + amrex::ignore_unused(pu, eb_bnd_normal_arr, i, j, k); +#endif +} +#endif //AMREX_USE_EB + /* This computes the scale_fac (used for setting the particle weights) on a on area basis - (used for flux injection). + (used for flux injection from a plane). */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -amrex::Real compute_scale_fac_area (const amrex::GpuArray& dx, +amrex::Real compute_scale_fac_area_plane (const amrex::GpuArray& dx, const amrex::Real num_ppc_real, const int flux_normal_axis) { using namespace amrex::literals; amrex::Real scale_fac = AMREX_D_TERM(dx[0],*dx[1],*dx[2])/num_ppc_real; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index d90befdd3a7..bd589ced73d 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -147,29 +147,6 @@ namespace return z0; } - struct PDim3 { - ParticleReal x, y, z; - - AMREX_GPU_HOST_DEVICE - PDim3(const amrex::XDim3& a): - x{static_cast(a.x)}, - y{static_cast(a.y)}, - z{static_cast(a.z)} - {} - - AMREX_GPU_HOST_DEVICE - ~PDim3() = default; - - AMREX_GPU_HOST_DEVICE - PDim3(PDim3 const &) = default; - AMREX_GPU_HOST_DEVICE - PDim3& operator=(PDim3 const &) = default; - AMREX_GPU_HOST_DEVICE - PDim3(PDim3&&) = default; - AMREX_GPU_HOST_DEVICE - PDim3& operator=(PDim3&&) = default; - }; - AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE XDim3 getCellCoords (const GpuArray& lo_corner, const GpuArray& dx, @@ -1428,6 +1405,22 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, const auto dx = geom.CellSizeArray(); const auto problo = geom.ProbLoArray(); +#ifdef AMREX_USE_EB + bool const inject_from_eb = plasma_injector.m_inject_from_eb; // whether to inject from EB or from a plane + // Extract data structures for embedded boundaries + amrex::FabArray const* eb_flag = nullptr; + amrex::MultiCutFab const* eb_bnd_area = nullptr; + amrex::MultiCutFab const* eb_bnd_normal = nullptr; + amrex::MultiCutFab const* eb_bnd_cent = nullptr; + if (inject_from_eb) { + amrex::EBFArrayBoxFactory const& eb_box_factory = WarpX::GetInstance().fieldEBFactory(0); + eb_flag = &eb_box_factory.getMultiEBCellFlagFab(); + eb_bnd_area = &eb_box_factory.getBndryArea(); + eb_bnd_normal = &eb_box_factory.getBndryNormal(); + eb_bnd_cent = &eb_box_factory.getBndryCent(); + } +#endif + amrex::LayoutData* cost = WarpX::getCosts(0); // Create temporary particle container to which particles will be added; @@ -1487,9 +1480,20 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, RealBox overlap_realbox; Box overlap_box; IntVect shifted; - const bool no_overlap = find_overlap_flux(tile_realbox, part_realbox, dx, problo, plasma_injector, overlap_realbox, overlap_box, shifted); - if (no_overlap) { - continue; // Go to the next tile +#ifdef AMREX_USE_EB + if (inject_from_eb) { + // Injection from EB + const amrex::FabType fab_type = (*eb_flag)[mfi].getType(tile_box); + if (fab_type == amrex::FabType::regular) { continue; } // Go to the next tile + if (fab_type == amrex::FabType::covered) { continue; } // Go to the next tile + overlap_box = tile_box; + overlap_realbox = part_realbox; + } else +#endif + { + // Injection from a plane + const bool no_overlap = find_overlap_flux(tile_realbox, part_realbox, dx, problo, plasma_injector, overlap_realbox, overlap_box, shifted); + if (no_overlap) { continue; } // Go to the next tile } const int grid_id = mfi.index(); @@ -1509,24 +1513,57 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, if (refine_injection) { fine_overlap_box = overlap_box & amrex::shift(fine_injection_box, -shifted); } + +#ifdef AMREX_USE_EB + // Extract data structures for embedded boundaries + amrex::Array4::value_type> eb_flag_arr; + amrex::Array4 eb_bnd_area_arr; + amrex::Array4 eb_bnd_normal_arr; + amrex::Array4 eb_bnd_cent_arr; + if (inject_from_eb) { + eb_flag_arr = eb_flag->array(mfi); + eb_bnd_area_arr = eb_bnd_area->array(mfi); + eb_bnd_normal_arr = eb_bnd_normal->array(mfi); + eb_bnd_cent_arr = eb_bnd_cent->array(mfi); + } +#endif + amrex::ParallelForRNG(overlap_box, [=] AMREX_GPU_DEVICE (int i, int j, int k, amrex::RandomEngine const& engine) noexcept { const IntVect iv(AMREX_D_DECL(i, j, k)); amrex::ignore_unused(j,k); - auto lo = getCellCoords(overlap_corner, dx, {0._rt, 0._rt, 0._rt}, iv); - auto hi = getCellCoords(overlap_corner, dx, {1._rt, 1._rt, 1._rt}, iv); - - if (flux_pos->overlapsWith(lo, hi)) + // Determine the number of macroparticles to inject in this cell (num_ppc_int) +#ifdef AMREX_USE_EB + amrex::Real num_ppc_real_in_this_cell = num_ppc_real; // user input: number of macroparticles per cell + if (inject_from_eb) { + // Injection from EB + // Skip cells that are not partially covered by the EB + if (eb_flag_arr(i,j,k).isRegular() || eb_flag_arr(i,j,k).isCovered()) { return; } + // Scale by the (normalized) area of the EB surface in this cell + num_ppc_real_in_this_cell *= eb_bnd_area_arr(i,j,k); + } else +#else + amrex::Real const num_ppc_real_in_this_cell = num_ppc_real; // user input: number of macroparticles per cell +#endif { - auto index = overlap_box.index(iv); - int r = 1; - if (fine_overlap_box.ok() && fine_overlap_box.contains(iv)) { - r = compute_area_weights(rrfac, flux_normal_axis); - } - const int num_ppc_int = static_cast(num_ppc_real*r + amrex::Random(engine)); - pcounts[index] = num_ppc_int; + // Injection from a plane + auto lo = getCellCoords(overlap_corner, dx, {0._rt, 0._rt, 0._rt}, iv); + auto hi = getCellCoords(overlap_corner, dx, {1._rt, 1._rt, 1._rt}, iv); + // Skip cells that do not overlap with the plane + if (!flux_pos->overlapsWith(lo, hi)) { return; } } + + auto index = overlap_box.index(iv); + // Take into account refined injection region + int r = 1; + if (fine_overlap_box.ok() && fine_overlap_box.contains(iv)) { + r = compute_area_weights(rrfac, flux_normal_axis); + } + const int num_ppc_int = static_cast(num_ppc_real_in_this_cell*r + amrex::Random(engine)); + pcounts[index] = num_ppc_int; + + amrex::ignore_unused(j,k); }); // Max number of new particles. All of them are created, @@ -1596,7 +1633,15 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, amrex::ignore_unused(j,k); const auto index = overlap_box.index(iv); - Real scale_fac = compute_scale_fac_area(dx, num_ppc_real, flux_normal_axis); + Real scale_fac; +#ifdef AMREX_USE_EB + if (inject_from_eb) { + scale_fac = compute_scale_fac_area_eb(dx, num_ppc_real, eb_bnd_normal_arr, i, j, k ); + } else +#endif + { + scale_fac = compute_scale_fac_area_plane(dx, num_ppc_real, flux_normal_axis); + } if (fine_overlap_box.ok() && fine_overlap_box.contains(iv)) { scale_fac /= compute_area_weights(rrfac, flux_normal_axis); @@ -1607,13 +1652,32 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, const long ip = poffset[index] + i_part; pa_idcpu[ip] = amrex::SetParticleIDandCPU(pid+ip, cpuid); - // This assumes the flux_pos is of type InjectorPositionRandomPlane - const XDim3 r = (fine_overlap_box.ok() && fine_overlap_box.contains(iv)) ? - // In the refined injection region: use refinement ratio `rrfac` - flux_pos->getPositionUnitBox(i_part, rrfac, engine) : - // Otherwise: use 1 as the refinement ratio - flux_pos->getPositionUnitBox(i_part, amrex::IntVect::TheUnitVector(), engine); - auto pos = getCellCoords(overlap_corner, dx, r, iv); + // Determine the position of the particle within the cell + XDim3 pos; + XDim3 r; +#ifdef AMREX_USE_EB + if (inject_from_eb) { +#if defined(WARPX_DIM_3D) + pos.x = overlap_corner[0] + (iv[0] + 0.5_rt + eb_bnd_cent_arr(i,j,k,0))*dx[0]; + pos.y = overlap_corner[1] + (iv[1] + 0.5_rt + eb_bnd_cent_arr(i,j,k,1))*dx[1]; + pos.z = overlap_corner[2] + (iv[2] + 0.5_rt + eb_bnd_cent_arr(i,j,k,2))*dx[2]; +#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + pos.x = overlap_corner[0] + (iv[0] + 0.5_rt + eb_bnd_cent_arr(i,j,k,0))*dx[0]; + pos.y = 0.0_rt; + pos.z = overlap_corner[1] + (iv[1] + 0.5_rt + eb_bnd_cent_arr(i,j,k,1))*dx[1]; +#endif + } else +#endif + { + // Injection from a plane + // This assumes the flux_pos is of type InjectorPositionRandomPlane + r = (fine_overlap_box.ok() && fine_overlap_box.contains(iv)) ? + // In the refined injection region: use refinement ratio `rrfac` + flux_pos->getPositionUnitBox(i_part, rrfac, engine) : + // Otherwise: use 1 as the refinement ratio + flux_pos->getPositionUnitBox(i_part, amrex::IntVect::TheUnitVector(), engine); + pos = getCellCoords(overlap_corner, dx, r, iv); + } auto ppos = PDim3(pos); // inj_mom would typically be InjectorMomentumGaussianFlux @@ -1660,6 +1724,15 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, continue; } +#ifdef AMREX_USE_EB + if (inject_from_eb) { + // Injection from EB: rotate momentum according to the normal of the EB surface + // (The above code initialized the momentum by assuming that z is the direction + // normal to the EB surface. Thus we need to rotate from z to the normal.) + rotate_momentum_eb(pu, eb_bnd_normal_arr, i, j , k); + } +#endif + #if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_RCYLINDER) // Conversion from cylindrical to Cartesian coordinates // Replace the x and y, setting an angle theta. @@ -1682,7 +1755,11 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, const amrex::Real radial_position = ppos.x; ppos.x = radial_position*cos_theta; ppos.y = radial_position*sin_theta; - if (loc_flux_normal_axis != 2) { + if ((loc_flux_normal_axis != 2) +#ifdef AMREX_USE_EB + || (inject_from_eb) +#endif + ) { // Rotate the momentum // This because, when the flux direction is e.g. "r" // the `inj_mom` objects generates a v*Gaussian distribution diff --git a/Source/Python/Particles/MultiParticleContainer.cpp b/Source/Python/Particles/MultiParticleContainer.cpp index e709f0950b4..7b3b114b080 100644 --- a/Source/Python/Particles/MultiParticleContainer.cpp +++ b/Source/Python/Particles/MultiParticleContainer.cpp @@ -42,5 +42,12 @@ i_lens: int strength_E, strength_B: floats The electric and magnetic focusing strength of the lens)pbdoc" ) + + .def("get_charge_density", + [](MultiParticleContainer& mpc, int lev, bool local) { + return mpc.GetChargeDensity(lev, local); + }, + py::arg("lev"), py::arg("local") + ) ; } diff --git a/Source/WarpX.H b/Source/WarpX.H index 14178891d65..b27c8a32009 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -689,7 +689,7 @@ public: void DampJPML (int lev, PatchType patch_type); void CopyJPML (); - bool isAnyBoundaryPML(); + /** True if any of the particle boundary condition type is Thermal */ static bool isAnyParticleBoundaryThermal(); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 26719916853..31e7923c5cf 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -195,6 +195,22 @@ amrex::IntVect m_rho_nodal_flag; WarpX* WarpX::m_instance = nullptr; +namespace +{ + + [[nodiscard]] bool + isAnyBoundaryPML( + const amrex::Array& field_boundary_lo, + const amrex::Array& field_boundary_hi) + { + constexpr auto is_pml = [](const FieldBoundaryType fbt) {return (fbt == FieldBoundaryType::PML);}; + const auto is_any_pml = + std::any_of(field_boundary_lo.begin(), field_boundary_lo.end(), is_pml) || + std::any_of(field_boundary_hi.begin(), field_boundary_hi.end(), is_pml); + return is_any_pml; + } +} + void WarpX::MakeWarpX () { ParseGeometryInput(); @@ -889,16 +905,16 @@ WarpX::ReadParameters () ); } -#if defined(WARPX_DIM_RZ) - WARPX_ALWAYS_ASSERT_WITH_MESSAGE( field_boundary_lo[1] != FieldBoundaryType::PML && field_boundary_hi[1] != FieldBoundaryType::PML, - "PML are not implemented in RZ geometry along z; please set a different boundary condition using boundary.field_lo and boundary.field_hi."); -#endif #if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE) - WARPX_ALWAYS_ASSERT_WITH_MESSAGE( isAnyBoundaryPML() == false || electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD, + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( ::isAnyBoundaryPML(field_boundary_lo, field_boundary_hi) == false || electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD, "PML are are only implemented with Cartesian geometry with FDTD; please set a different boundary condition using boundary.field_lo and boundary.field_hi."); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( (do_pml_dive_cleaning == false && do_pml_divb_cleaning == false), "do_pml_dive_cleaning and do_pml_divb_cleaning are only implemented in Cartesian geometry." ); #endif +#if defined(WARPX_DIM_RZ) + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( field_boundary_lo[1] != FieldBoundaryType::PML && field_boundary_hi[1] != FieldBoundaryType::PML, + "PML are not implemented in RZ geometry along z; please set a different boundary condition using boundary.field_lo and boundary.field_hi."); +#endif WARPX_ALWAYS_ASSERT_WITH_MESSAGE( (do_pml_j_damping==0)||(do_pml_in_domain==1), @@ -2035,7 +2051,7 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d safe_guard_cells, WarpX::do_multi_J, WarpX::fft_do_time_averaging, - WarpX::isAnyBoundaryPML(), + ::isAnyBoundaryPML(field_boundary_lo, field_boundary_hi), WarpX::do_pml_in_domain, WarpX::pml_ncell, this->refRatio(), @@ -2774,7 +2790,7 @@ void WarpX::AllocLevelSpectralSolverRZ (amrex::Vector