From bacbef747f9320a6288ef461bd8ace1c669ae665 Mon Sep 17 00:00:00 2001 From: Edoardo Zoni <59625522+EZoni@users.noreply.github.com> Date: Mon, 7 Oct 2024 16:43:13 -0700 Subject: [PATCH 01/16] AMReX/pyAMReX/PICSAR: weekly update (#5369) Weekly update to latest AMReX/pyAMReX/PICSAR. ```console ./Tools/Release/updateAMReX.py ./Tools/Release/updatepyAMReX.py ./Tools/Release/updatePICSAR.py ``` Co-authored-by: Axel Huebl --- .github/workflows/cuda.yml | 2 +- cmake/dependencies/AMReX.cmake | 2 +- cmake/dependencies/pyAMReX.cmake | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) 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/cmake/dependencies/AMReX.cmake b/cmake/dependencies/AMReX.cmake index 91340066803..51ad361276a 100644 --- a/cmake/dependencies/AMReX.cmake +++ b/cmake/dependencies/AMReX.cmake @@ -279,7 +279,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 "24.10" +set(WarpX_amrex_branch "e1222803739ed2342b9ff6fc2d57316ff0d6cb0c" CACHE STRING "Repository branch for WarpX_amrex_repo if(WarpX_amrex_internal)") diff --git a/cmake/dependencies/pyAMReX.cmake b/cmake/dependencies/pyAMReX.cmake index f7b905c32c3..9543dac2ee2 100644 --- a/cmake/dependencies/pyAMReX.cmake +++ b/cmake/dependencies/pyAMReX.cmake @@ -74,7 +74,7 @@ option(WarpX_pyamrex_internal "Download & build pyAMReX" ON) set(WarpX_pyamrex_repo "https://github.com/AMReX-Codes/pyamrex.git" CACHE STRING "Repository URI to pull and build pyamrex from if(WarpX_pyamrex_internal)") -set(WarpX_pyamrex_branch "24.10" +set(WarpX_pyamrex_branch "3699781e4284921f9ccdbbbbc57169ff79c0de20" CACHE STRING "Repository branch for WarpX_pyamrex_repo if(WarpX_pyamrex_internal)") From 46d88b8311c750efc8f2348a5a25a9fbc12558f2 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 8 Oct 2024 02:21:05 +0000 Subject: [PATCH 02/16] [pre-commit.ci] pre-commit autoupdate (#5373) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit updates: - [github.com/pre-commit/pre-commit-hooks: v4.6.0 → v5.0.0](https://github.com/pre-commit/pre-commit-hooks/compare/v4.6.0...v5.0.0) - [github.com/astral-sh/ruff-pre-commit: v0.6.8 → v0.6.9](https://github.com/astral-sh/ruff-pre-commit/compare/v0.6.8...v0.6.9) Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- .pre-commit-config.yaml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 From 894a699c11d165fdaff53fcc3cd051eeadbdf367 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Mon, 7 Oct 2024 19:46:16 -0700 Subject: [PATCH 03/16] Doc: Governance GitHub Team Links (#5374) Fix the link in the governance doc to the steering committee (formerly: admin) and technical committee (formerly: maintainers). This was a rename when we adopted the governance doc #4743 and fixes two broken links. cc @ECP-WarpX/warpx-steering-committee (for approval) cc @ECP-WarpX/warpx-technical-committee (FYI) --- GOVERNANCE.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 ^^^^ From b7108967b5068c5d942d0ba9e0285efd0d07aa05 Mon Sep 17 00:00:00 2001 From: Edoardo Zoni <59625522+EZoni@users.noreply.github.com> Date: Tue, 8 Oct 2024 14:21:48 -0700 Subject: [PATCH 04/16] CI: update docs and reset tool for checksums (#5372) The major part of this PR is about updating the docs so that it is a bit easier for developers to connect the section on testing with the section on checksums. Here's a couple of screenshots showing the new content organization for the testing and checksums sections:

As part of this PR, I also update the tool that we have to reset checksums locally based on the Azure output. The update is necessary due to a change in #5297, as noted in https://github.com/ECP-WarpX/WarpX/pull/5297#discussion_r1787006065, given that new checksum files are now displayed as follows: ``` New checksums file test_2d_langmuir_multi.json: { "lev=0": { "Bx": 0.0, "By": 5.726296856755232, "Bz": 0.0, "Ex": 3751589134191.326, "Ey": 0.0, "Ez": 3751589134191.332, "jx": 1.0100623329922576e+16, "jy": 0.0, "jz": 1.0100623329922578e+16 }, "electrons": { "particle_momentum_x": 5.668407513430198e-20, "particle_momentum_y": 0.0, "particle_momentum_z": 5.668407513430198e-20, "particle_position_x": 0.6553599999999999, "particle_position_y": 0.65536, "particle_weight": 3200000000000000.5 }, "positrons": { "particle_momentum_x": 5.668407513430198e-20, "particle_momentum_y": 0.0, "particle_momentum_z": 5.668407513430198e-20, "particle_position_x": 0.6553599999999999, "particle_position_y": 0.65536, "particle_weight": 3200000000000000.5 } } ``` as opposed to the old way ``` ---------------- New file for test_2d_langmuir_multi: { "lev=0": { "Bx": 0.0, "By": 5.726296856755232, "Bz": 0.0, "Ex": 3751589134191.326, "Ey": 0.0, "Ez": 3751589134191.332, "jx": 1.0100623329922576e+16, "jy": 0.0, "jz": 1.0100623329922578e+16 }, "electrons": { "particle_momentum_x": 5.668407513430198e-20, "particle_momentum_y": 0.0, "particle_momentum_z": 5.668407513430198e-20, "particle_position_x": 0.6553599999999999, "particle_position_y": 0.65536, "particle_weight": 3200000000000000.5 }, "positrons": { "particle_momentum_x": 5.668407513430198e-20, "particle_momentum_y": 0.0, "particle_momentum_z": 5.668407513430198e-20, "particle_position_x": 0.6553599999999999, "particle_position_y": 0.65536, "particle_weight": 3200000000000000.5 } } ---------------- ``` To-do: - [x] Update docs - [x] Update tool --------- Co-authored-by: Axel Huebl --- Docs/source/developers/checksum.rst | 45 +++++----- Docs/source/developers/testing.rst | 47 +++++++---- Docs/source/developers/workflows.rst | 4 +- .../update_benchmarks_from_azure_output.py | 83 ++++++++----------- 4 files changed, 94 insertions(+), 85 deletions(-) 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..5fd4b498b07 100644 --- a/Docs/source/developers/testing.rst +++ b/Docs/source/developers/testing.rst @@ -3,33 +3,43 @@ 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:`Using checksums `. + +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/Tools/DevUtils/update_benchmarks_from_azure_output.py b/Tools/DevUtils/update_benchmarks_from_azure_output.py index b2be4d17a7b..bcff995b21a 100644 --- a/Tools/DevUtils/update_benchmarks_from_azure_output.py +++ b/Tools/DevUtils/update_benchmarks_from_azure_output.py @@ -1,4 +1,4 @@ -# Copyright 2023 Neil Zaim +# Copyright 2023 Neil Zaim, Edoardo Zoni # # This file is part of WarpX. # @@ -9,56 +9,45 @@ import sys """ -This Python script updates the Azure benchmarks automatically using a raw Azure output textfile -that is given as the first and only argument of the script. - -In the Azure output, we read the lines contained between -"New file for Test_Name:" -and the next occurrence of -"'----------------'" -And use these lines to update the benchmarks +This Python script updates the Azure benchmarks automatically using a raw +Azure output text file that is passed as command line argument of the script. """ -azure_output_filename = sys.argv[1] +# read path to Azure output text file +azure_output = sys.argv[1] -pattern_test_name = "New file for (?P[\w\-]*)" -closing_string = "----------------" -benchmark_path = "../../Regression/Checksum/benchmarks_json/" -benchmark_suffix = ".json" +# string to identify failing tests that require a checksums reset +new_checksums = "New checksums" +failing_test = "" -first_line_read = False -current_test = "" +# path of all checksums benchmark files +benchmark_path = "../../Regression/Checksum/benchmarks_json/" -with open(azure_output_filename, "r") as f: +with open(azure_output, "r") as f: + # find length of Azure prefix to be removed from each line, + # first line of Azure output starts with "##[section]Starting:" + first_line = f.readline() + prefix_length = first_line.find("#") + # loop over lines for line in f: - if current_test == "": - # Here we search lines that read, for example, - # "New file for LaserAcceleration_BTD" - # and we set current_test = "LaserAcceleration_BTD" - match_test_name = re.search(pattern_test_name, line) - if match_test_name: - current_test = match_test_name.group("testname") - new_file_string = "" - + # remove Azure prefix from line + line = line[prefix_length:] + if failing_test == "": + # no failing test found yet + if re.search(new_checksums, line): + # failing test found, set failing test name + failing_test = line[line.find("test_") : line.find(".json")] + json_file_string = "" else: - # We add each line to the new file string until we find the line containing - # "----------------" - # which indicates that we have read the new file entirely - - if closing_string not in line: - if not first_line_read: - # Raw Azure output comes with a prefix at the beginning of each line that we do - # not need here. The first line that we will read is the prefix followed by the - # "{" character, so we determine how long the prefix is by finding the last - # occurrence of the "{" character in this line. - azure_indent = line.rfind("{") - first_line_read = True - new_file_string += line[azure_indent:] - - else: - # We have read the new file entirely. Dump it in the json file. - new_file_json = json.loads(new_file_string) - json_filepath = benchmark_path + current_test + benchmark_suffix - with open(json_filepath, "w") as f_json: - json.dump(new_file_json, f_json, indent=2) - current_test = "" + # extract and dump new checksums of failing test + json_file_string += line + if line.startswith("}"): # end of new checksums + json_file = json.loads(json_file_string) + json_filename = failing_test + ".json" + json_filepath = benchmark_path + json_filename + print(f"\nDumping new checksums file {json_filename}:") + print(json_file_string) + with open(json_filepath, "w") as json_f: + json.dump(json_file, json_f, indent=2) + # reset to empty string to continue search of failing tests + failing_test = "" From 27181aa06b73e21b8c94c8650c753bfa309bd137 Mon Sep 17 00:00:00 2001 From: Edoardo Zoni <59625522+EZoni@users.noreply.github.com> Date: Wed, 9 Oct 2024 09:36:42 -0700 Subject: [PATCH 05/16] Docs: fix checksums section cross-reference (#5376) The checksums section title was changed to "Checksums on Tests" in the latest version of #5372, but the cross-reference in the testing section wasn't updated and still had the old name "Using checksums". --------- Co-authored-by: Axel Huebl --- Docs/source/developers/testing.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Docs/source/developers/testing.rst b/Docs/source/developers/testing.rst index 5fd4b498b07..111e3e7d7cb 100644 --- a/Docs/source/developers/testing.rst +++ b/Docs/source/developers/testing.rst @@ -1,6 +1,6 @@ .. _developers-testing: -Testing the code +Testing the Code ================ When proposing a code change, you want to make sure that @@ -123,7 +123,7 @@ An automated test typically consists of the following components: * analysis script; * checksum file. -To learn more about how to use checksums in automated tests, please see the corresponding section :ref:`Using checksums `. +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 `__. From 1d2910e276b02e4f1c4c1486b710a97a97776809 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Wed, 9 Oct 2024 11:49:34 -0700 Subject: [PATCH 06/16] CMake: Fix List of Pip Options (#5378) We were not yet able to pass lists of options to `pip` commands in our `pip` CMake helper targets. This fixes it. Follow-up to #2822 X-ref: https://github.com/spack/spack/pull/46765 --- .github/workflows/macos.yml | 17 ++++++++--------- CMakeLists.txt | 10 +++++++--- Docs/source/install/cmake.rst | 10 +++++----- 3 files changed, 20 insertions(+), 17 deletions(-) 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/CMakeLists.txt b/CMakeLists.txt index 5065baa0b6d..980b23183fd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -714,9 +714,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) @@ -739,7 +739,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 @@ -754,6 +755,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} ) @@ -771,6 +773,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 @@ -784,6 +787,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/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: From 65e82617839ac70b3f0951eab18fa04a06dfb93e Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Wed, 9 Oct 2024 17:31:50 -0500 Subject: [PATCH 07/16] Disable AMReX_LINEAR_SOLVER_INCFLO by default (#5364) We do not need to compile these linear solvers that are intended for incompressible flow solvers. This should speed up the build process a little bit. Introduced in https://github.com/AMReX-Codes/amrex/pull/4181 --- GNUmakefile | 3 +++ Source/ablastr/Make.package | 1 - Source/ablastr/fields/Make.package | 4 +++- cmake/dependencies/AMReX.cmake | 4 ++++ 4 files changed, 10 insertions(+), 2 deletions(-) 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/Source/ablastr/Make.package b/Source/ablastr/Make.package index b9ff3c72560..edbf43b7802 100644 --- a/Source/ablastr/Make.package +++ b/Source/ablastr/Make.package @@ -1,4 +1,3 @@ -#CEXE_sources += ParticleBoundaries.cpp include $(WARPX_HOME)/Source/ablastr/coarsen/Make.package include $(WARPX_HOME)/Source/ablastr/math/Make.package diff --git a/Source/ablastr/fields/Make.package b/Source/ablastr/fields/Make.package index 727a17b6de8..7441a6a1238 100644 --- a/Source/ablastr/fields/Make.package +++ b/Source/ablastr/fields/Make.package @@ -1,5 +1,7 @@ + +CEXE_sources += MultiFabRegister.cpp + ifeq ($(USE_FFT),TRUE) - CEXE_sources += MultiFabRegister.cpp ifeq ($(DIM),3) CEXE_sources += IntegratedGreenFunctionSolver.cpp endif diff --git a/cmake/dependencies/AMReX.cmake b/cmake/dependencies/AMReX.cmake index 51ad361276a..6513841f327 100644 --- a/cmake/dependencies/AMReX.cmake +++ b/cmake/dependencies/AMReX.cmake @@ -92,6 +92,8 @@ macro(find_amrex) set(AMReX_PARTICLES ON CACHE INTERNAL "") set(AMReX_PROBINIT OFF CACHE INTERNAL "") set(AMReX_TINY_PROFILE ON CACHE BOOL "") + set(AMReX_LINEAR_SOLVERS_EM ON CACHE INTERNAL "") + set(AMReX_LINEAR_SOLVER_INCFLO OFF CACHE INTERNAL "") if(WarpX_ASCENT OR WarpX_SENSEI) set(AMReX_GPU_RDC ON CACHE BOOL "") @@ -200,6 +202,8 @@ macro(find_amrex) mark_as_advanced(AMReX_HYPRE) mark_as_advanced(AMReX_IPO) mark_as_advanced(AMReX_LINEAR_SOLVERS) + mark_as_advanced(AMReX_LINEAR_SOLVERS_INCFLO) + mark_as_advanced(AMReX_LINEAR_SOLVERS_EM) mark_as_advanced(AMReX_MEM_PROFILE) mark_as_advanced(AMReX_MPI) mark_as_advanced(AMReX_MPI_THREAD_MULTIPLE) From 89dc850e6c735b21221f61955204371bf367b773 Mon Sep 17 00:00:00 2001 From: Roelof Groenewald <40245517+roelof-groenewald@users.noreply.github.com> Date: Wed, 9 Oct 2024 19:25:42 -0700 Subject: [PATCH 08/16] Expose `MultiParticleContainer.GetChargeDensity` to Python (#5382) Signed-off-by: roelof-groenewald --- Source/Python/Particles/MultiParticleContainer.cpp | 7 +++++++ 1 file changed, 7 insertions(+) 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") + ) ; } From c045eaf4a525605e69336feb40aeed7bf50bd63a Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Wed, 9 Oct 2024 19:26:28 -0700 Subject: [PATCH 09/16] CMake: No FFTW Needed for SYCL anymore (#5380) We do not need FFTW3 anymore to do FFTs on SYCL GPUs. Follow-up to #5127 X-ref: https://github.com/spack/spack/pull/46765 --- .github/workflows/intel.yml | 1 + cmake/dependencies/FFT.cmake | 16 ++++++++++++---- 2 files changed, 13 insertions(+), 4 deletions(-) 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/cmake/dependencies/FFT.cmake b/cmake/dependencies/FFT.cmake index 571006e8530..df0ef11ae53 100644 --- a/cmake/dependencies/FFT.cmake +++ b/cmake/dependencies/FFT.cmake @@ -48,14 +48,20 @@ if(ABLASTR_FFT) # # cuFFT (CUDA) - # TODO: check if `find_package` search works + if(WarpX_COMPUTE STREQUAL CUDA) + # nothing to do (cuFFT is part of the CUDA SDK) + # TODO: check if `find_package` search works for cuFFT # rocFFT (HIP) - if(WarpX_COMPUTE STREQUAL HIP) + elseif(WarpX_COMPUTE STREQUAL HIP) find_package(rocfft REQUIRED) - # FFTW (NOACC, OMP, SYCL) - elseif(NOT WarpX_COMPUTE STREQUAL CUDA) + elseif(WarpX_COMPUTE STREQUAL SYCL) + # nothing to do (oneMKL is part of oneAPI) + # TODO: check if `find_package` search works for oneMKL + + # FFTW (NOACC, OMP) + else() # On Windows, try searching for FFTW3(f)Config.cmake files first # Installed .pc files wrongly and unconditionally add -lm # https://github.com/FFTW/fftw3/issues/236 @@ -106,6 +112,8 @@ if(ABLASTR_FFT) warpx_make_third_party_includes_system(cufft FFT) elseif(WarpX_COMPUTE STREQUAL HIP) warpx_make_third_party_includes_system(roc::rocfft FFT) + elseif(WarpX_COMPUTE STREQUAL SYCL) + warpx_make_third_party_includes_system(AMReX::SYCL FFT) else() if(WarpX_FFTW_SEARCH STREQUAL CMAKE) warpx_make_third_party_includes_system(FFTW3::fftw3${HFFTWp} FFT) From b2840be3ccd9cd886a9d30eace7d6a01d7b6d1ff Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Wed, 9 Oct 2024 23:58:26 -0700 Subject: [PATCH 10/16] Fix CI: CodeQL Setup (#5385) Fix broken Python setup in CodeQL CI. --- .github/workflows/codeql.yml | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) 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 From a716670ba97e738456241791f933bbccfa3bdbce Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Thu, 10 Oct 2024 09:52:42 -0700 Subject: [PATCH 11/16] Generalize differential luminosity for photons (#5222) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The differential luminosity diagnostic was not valid for particles with mass 0. This PR generalizes the code for photons by expressing the center-of-mass energy with the 4-momentum: Screenshot 2024-09-06 at 9 27 33 PM which is valid for photons as well. I also slightly simplified the code that computes the term Screenshot 2024-09-06 at 9 30 45 PM I also added a test using photons. This is done by extracting an base input script from the existing luminosity test involving **electrons and positrons**, and creating an input script involving photons that leverages the base input script. --- Docs/source/usage/parameters.rst | 9 +-- Examples/Tests/diff_lumi_diag/CMakeLists.txt | 14 +++- Examples/Tests/diff_lumi_diag/analysis.py | 16 ++++- ..._test_3d_diff_lumi_diag => inputs_base_3d} | 29 +-------- .../inputs_test_3d_diff_lumi_diag_leptons | 31 +++++++++ .../inputs_test_3d_diff_lumi_diag_photons | 28 ++++++++ ...on => test_3d_diff_lumi_diag_leptons.json} | 0 .../test_3d_diff_lumi_diag_photons.json | 24 +++++++ .../ReducedDiags/DifferentialLuminosity.cpp | 64 +++++++++++++------ 9 files changed, 162 insertions(+), 53 deletions(-) rename Examples/Tests/diff_lumi_diag/{inputs_test_3d_diff_lumi_diag => inputs_base_3d} (81%) create mode 100644 Examples/Tests/diff_lumi_diag/inputs_test_3d_diff_lumi_diag_leptons create mode 100644 Examples/Tests/diff_lumi_diag/inputs_test_3d_diff_lumi_diag_photons rename Regression/Checksum/benchmarks_json/{test_3d_diff_lumi_diag.json => test_3d_diff_lumi_diag_leptons.json} (100%) create mode 100644 Regression/Checksum/benchmarks_json/test_3d_diff_lumi_diag_photons.json diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index b9d82d5014a..910ce448c0f 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -3466,14 +3466,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/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/Source/Diagnostics/ReducedDiags/DifferentialLuminosity.cpp b/Source/Diagnostics/ReducedDiags/DifferentialLuminosity.cpp index 59a32cf0545..ef5e0da6014 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); From 005ef775f2a35feefe69d4e1f7069c99416b088f Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Thu, 10 Oct 2024 12:52:30 -0700 Subject: [PATCH 12/16] SYCL: 1D EB Compile (#5384) Attempt to fix 1D SYCL EB compile errors (throw not allowed on device). X-ref: https://github.com/spack/spack/pull/46765#issuecomment-2403937237 --- Source/EmbeddedBoundary/DistanceToEB.H | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) 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; From 4434e87ff9ed3dca833f4afec9c610d57b006364 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Fri, 11 Oct 2024 00:07:37 +0200 Subject: [PATCH 13/16] Move `isAnyBoundaryPML` to Warpx.cpp (#5353) `isAnyBoundaryPML` is used only inside `WarpX.cpp`. It does not need to be a member function of the WarpX class and it can be moved into an anonymous namespace inside `WarpX.cpp`. --- Source/WarpX.H | 2 +- Source/WarpX.cpp | 32 +++++++++++++++++++------------- 2 files changed, 20 insertions(+), 14 deletions(-) diff --git a/Source/WarpX.H b/Source/WarpX.H index c61fb92315f..bad63cd44d9 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 250bab273d0..d1e3108e32a 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(); @@ -878,7 +894,7 @@ WarpX::ReadParameters () } #ifdef WARPX_DIM_RZ - 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 not implemented in RZ geometry with FDTD; please set a different boundary condition using boundary.field_lo and boundary.field_hi."); 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."); @@ -2014,7 +2030,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(), @@ -2742,7 +2758,7 @@ void WarpX::AllocLevelSpectralSolverRZ (amrex::Vector Date: Thu, 10 Oct 2024 16:31:33 -0700 Subject: [PATCH 14/16] Implement injection of particles from the embedded boundary (#5208) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit # Overview This PR implements flux injection of particles from the embedded boundary. It also adds a test that emits particles from a sphere in 3D as represented here: ![movie](https://github.com/user-attachments/assets/1e76cf87-fd7d-4fa3-8c83-363956226a42) as well as RZ and 2D versions of this test. (In 2D, the particles are emitted from a cylinder.) As can be seen in the above movie, particles are emitted from a single point within each cell (the centroid of the EB), instead of being emitted uniformly on the surface of the EB within the cell. This could be improved in a future PR. The implementation as well as the user interface largely re-use the infrastructure for the flux injection from a plane. However, as a result, the user interface is perhaps not very intuitive. In particular, when specify the velocity distribution, `uz` represents the direction normal to the EB while `ux`, `uy` represent the tangential directions. This again will be improved in follow-up PR. # Follow-up PRs - [ ] Change the interface of `gaussianflux` so as to specify the tangential and normal distribution. In other words, instead of: ``` 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 ``` we would do: ``` electron.momentum_distribution_type = gaussianflux electron.u_tangential_th = 0.01 # Tangential to the emitting surface electron.u_normal_th = 0.1 # Normal to the emitting surface electron.u_normal_m = 0.07 ``` - [ ] Change the interface so that the user does not need to specify the number of macroparticles per cell (which is problematic for EB, since difference cell contain different EB surface, and should in general emit different numbers of macroparticles). Instead, we would specify the weight of macroparticles, i.e. instead of ``` electron.injection_style = NFluxPerCell electron.num_particles_per_cell = 100 electron.flux_function(x,y,z,t) = "1.” ``` we would do ``` electron.injection_style = NFluxPerCell electron.flux_macroweight = 200 # Number of physical particles per macroparticle electron.flux_function(x,y,z,t) = "4e12” # Number of physical particles emitted per unit time and surface ``` - [ ] Add a way for the user to specify the total flux across the whole emitting surface Example: ``` electron.flux_function(x,y,z,t) = "(x>-1)*(x<1)" electron.total_flux = 4e12 # physical particle / second (not per unit area) ``` (In that case, `flux_function` would be rescaled internally by WarpX so as to emit the right number of particles.) - [ ] Add PICMI interface - [ ] Emit the particles uniformly from the surface of the EB within one cell --- Docs/source/usage/parameters.rst | 13 +- Examples/Tests/flux_injection/CMakeLists.txt | 30 ++++ .../analysis_flux_injection_from_eb.py | 161 +++++++++++++++++ .../Tests/flux_injection/inputs_base_from_eb | 42 +++++ .../inputs_test_2d_flux_injection_from_eb | 13 ++ .../inputs_test_3d_flux_injection_from_eb | 13 ++ .../inputs_test_rz_flux_injection_from_eb | 15 ++ .../test_2d_flux_injection_from_eb.json | 11 ++ .../test_3d_flux_injection_from_eb.json | 12 ++ .../test_rz_flux_injection_from_eb.json | 12 ++ Source/Initialization/PlasmaInjector.H | 2 + Source/Initialization/PlasmaInjector.cpp | 78 ++++---- Source/Particles/AddPlasmaUtilities.H | 139 +++++++++++++- .../Particles/PhysicalParticleContainer.cpp | 169 +++++++++++++----- 14 files changed, 627 insertions(+), 83 deletions(-) create mode 100755 Examples/Tests/flux_injection/analysis_flux_injection_from_eb.py create mode 100644 Examples/Tests/flux_injection/inputs_base_from_eb create mode 100644 Examples/Tests/flux_injection/inputs_test_2d_flux_injection_from_eb create mode 100644 Examples/Tests/flux_injection/inputs_test_3d_flux_injection_from_eb create mode 100644 Examples/Tests/flux_injection/inputs_test_rz_flux_injection_from_eb create mode 100644 Regression/Checksum/benchmarks_json/test_2d_flux_injection_from_eb.json create mode 100644 Regression/Checksum/benchmarks_json/test_3d_flux_injection_from_eb.json create mode 100644 Regression/Checksum/benchmarks_json/test_rz_flux_injection_from_eb.json diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index 910ce448c0f..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`) 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/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_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/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 3d846375a99..76bb7a5be42 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" @@ -303,50 +304,65 @@ 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; #ifdef WARPX_DIM_RZ - 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 (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; + } #ifdef WARPX_DIM_3D - const std::string flux_normal_axis_help = "'x', 'y', or 'z'."; + const std::string flux_normal_axis_help = "'x', 'y', or 'z'."; #else # ifdef WARPX_DIM_RZ - const std::string flux_normal_axis_help = "'r' or 'z'."; + const std::string flux_normal_axis_help = "'r' or 'z'."; # elif WARPX_DIM_XZ - const std::string flux_normal_axis_help = "'x' or 'z'."; + const std::string flux_normal_axis_help = "'x' or 'z'."; # else - const std::string flux_normal_axis_help = "'z'."; + const std::string flux_normal_axis_help = "'z'."; # endif -#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."); + #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."); + } + // 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 8f0489e3921..bb05d7be3c8 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, @@ -71,12 +94,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 26f9fee38d3..7c70c9a35c4 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, @@ -1371,6 +1348,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; @@ -1428,9 +1421,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(); @@ -1450,24 +1454,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, @@ -1537,7 +1574,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); @@ -1548,13 +1593,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 @@ -1595,6 +1659,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 + #ifdef WARPX_DIM_RZ // Conversion from cylindrical to Cartesian coordinates // Replace the x and y, setting an angle theta. @@ -1610,7 +1683,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 From 3cda2c11e1d92ccd4b90c57debc9b399fa1978e6 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Thu, 10 Oct 2024 17:30:11 -0700 Subject: [PATCH 15/16] Avoid interpolating from guard cells in BTD (#5342) BTD diagnostics sometimes show artifacts at the edge of the range of collected data. (See for instance the red curve below.) My understanding is that this happens because the BTD collection planes use data from the guard cells outside of the simulation domain, when interpolating fields that are then used for the Lorentz back-transform. The guard cell data may not be physically correct (e.g. it may not have the right cancellation between `E` and `B`), and could thus cause this artifact. This PR avoids this issue by prevents the collection planes to collect data when it is half a cell from the edge of the simulation domain. See the example below, taken from https://github.com/ECP-WarpX/WarpX/pull/5337 (plot of the laser field, from the BTD diagnostic) ![Figure 40](https://github.com/user-attachments/assets/e4549856-4182-4a87-aa26-2d3bc6ac8e2c) The BTD diagnostics values are identical with this PR, except for the problematic point appearing at the edge of the domain. --- .../test_2d_rigid_injection_btd.json | 30 +++++++++---------- Source/Diagnostics/BTDiagnostics.cpp | 7 +++-- 2 files changed, 20 insertions(+), 17 deletions(-) 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/Source/Diagnostics/BTDiagnostics.cpp b/Source/Diagnostics/BTDiagnostics.cpp index 631de298861..f10f337a1f1 100644 --- a/Source/Diagnostics/BTDiagnostics.cpp +++ b/Source/Diagnostics/BTDiagnostics.cpp @@ -999,12 +999,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 ); From dc68659c1b2d6a689f4b602bddd42f9099fb5928 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Thu, 10 Oct 2024 17:30:47 -0700 Subject: [PATCH 16/16] Update BackTransformed diagnostics to take into account arbitrary moving window velocity (#5341) In the `development` branch, the `BackTransformed` diagnostics assume that the moving window moves exactly at the speed of light. This PR generalizes the code for arbitrary moving window velocity. This PR does not add an automated test, but the upcoming PR #5337 will add a test which features a moving window with a speed different than `c`. This is a follow-up of #5226, which modified the transformation of the simulation box coordinates for arbitrary moving window velocity, but did not yet update the `BackTransformed` diagnostic code. --- Docs/source/usage/faq.rst | 4 ++-- Source/Diagnostics/BTDiagnostics.H | 1 + Source/Diagnostics/BTDiagnostics.cpp | 36 +++++++++++++++------------- 3 files changed, 22 insertions(+), 19 deletions(-) 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/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 f10f337a1f1..312bbc7ec45 100644 --- a/Source/Diagnostics/BTDiagnostics.cpp +++ b/Source/Diagnostics/BTDiagnostics.cpp @@ -69,6 +69,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. @@ -138,7 +139,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 @@ -149,22 +150,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) { @@ -256,7 +256,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, @@ -338,13 +338,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 @@ -403,9 +405,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