From 5f15bbad22bcb01c6a3baa9ace1db916b9b17af5 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 28 Feb 2024 17:30:37 -0800 Subject: [PATCH 01/16] Update test that checks absorption of particles in PML (#4733) * Update test that checks absorption of particles in PML * Remove duplicated test file * Update dimensions of data --- .../analysis_particles_in_pml.py | 10 ++- .../analysis_particles_in_pml_2dmr.py | 68 ------------------- Regression/WarpX-tests.ini | 2 +- 3 files changed, 10 insertions(+), 70 deletions(-) delete mode 100755 Examples/Tests/particles_in_pml/analysis_particles_in_pml_2dmr.py diff --git a/Examples/Tests/particles_in_pml/analysis_particles_in_pml.py b/Examples/Tests/particles_in_pml/analysis_particles_in_pml.py index 85878d3add9..1d1a8959edd 100755 --- a/Examples/Tests/particles_in_pml/analysis_particles_in_pml.py +++ b/Examples/Tests/particles_in_pml/analysis_particles_in_pml.py @@ -30,8 +30,16 @@ filename = sys.argv[1] ds = yt.load( filename ) +# When extracting the fields, choose the right dimensions +dimensions = [ n_pts for n_pts in ds.domain_dimensions ] +if ds.max_level == 1: + dimensions[0] *= 2 + dimensions[1] *= 2 + if ds.dimensionality == 3: + dimensions[2] *= 2 + # Check that the field is low enough -ad0 = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions) +ad0 = ds.covering_grid(level=ds.max_level, left_edge=ds.domain_left_edge, dims=dimensions) Ex_array = ad0[('mesh','Ex')].to_ndarray() Ey_array = ad0[('mesh','Ey')].to_ndarray() Ez_array = ad0[('mesh','Ez')].to_ndarray() diff --git a/Examples/Tests/particles_in_pml/analysis_particles_in_pml_2dmr.py b/Examples/Tests/particles_in_pml/analysis_particles_in_pml_2dmr.py deleted file mode 100755 index 81a981e70a6..00000000000 --- a/Examples/Tests/particles_in_pml/analysis_particles_in_pml_2dmr.py +++ /dev/null @@ -1,68 +0,0 @@ -#!/usr/bin/env python3 - -# Copyright 2019-2020 Luca Fedeli, Maxence Thevenet, Remi Lehe -# -# -# This file is part of WarpX. -# -# License: BSD-3-Clause-LBNL - -""" -This script tests the absorption of particles in the PML. - -The input file inputs_2d/inputs is used: it features a positive and a -negative particle, going in opposite direction and eventually -leaving the box. This script tests that the field in the box -is close to 0 once the particles have left. With regular -PML, this test fails, since the particles leave a spurious -charge, with associated fields, behind them. -""" -import os -import sys - -import yt - -yt.funcs.mylog.setLevel(0) -sys.path.insert(1, '../../../../warpx/Regression/Checksum/') -import checksumAPI - -# Open plotfile specified in command line -filename = sys.argv[1] -ds = yt.load( filename ) - -# Check that the field is low enough -ad0 = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions) -Ex_array_lev0 = ad0[('mesh','Ex')].to_ndarray() -Ey_array_lev0 = ad0[('mesh','Ey')].to_ndarray() -Ez_array_lev0 = ad0[('mesh','Ez')].to_ndarray() -max_Efield_lev0 = max(Ex_array_lev0.max(), Ey_array_lev0.max(), Ez_array_lev0.max()) -print( "max_Efield level0 = %s" %max_Efield_lev0 ) - -ad1 = ds.covering_grid(level=1, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions) -Ex_array_lev1 = ad1[('mesh','Ex')].to_ndarray() -Ey_array_lev1 = ad1[('mesh','Ey')].to_ndarray() -Ez_array_lev1 = ad1[('mesh','Ez')].to_ndarray() -max_Efield_lev1 = max(Ex_array_lev1.max(), Ey_array_lev1.max(), Ez_array_lev1.max()) -print( "max_Efield level1 = %s" %max_Efield_lev1 ) - -# The field associated with the particle does not have -# the same amplitude in 2d and 3d -if ds.dimensionality == 2: - if ds.max_level == 0: - tolerance_abs = 0.0003 - elif ds.max_level == 1: - tolerance_abs = 0.0006 -elif ds.dimensionality == 3: - if ds.max_level == 0: - tolerance_abs = 10 - elif ds.max_level == 1: - tolerance_abs = 110 -else: - raise ValueError("Unknown dimensionality") - -print("tolerance_abs: " + str(tolerance_abs)) -assert max_Efield_lev0 < tolerance_abs -assert max_Efield_lev1 < tolerance_abs - -test_name = os.path.split(os.getcwd())[1] -checksumAPI.evaluate_checksum(test_name, filename) diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index b25d41db5d2..128495195da 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -2525,7 +2525,7 @@ numthreads = 1 compileTest = 0 doVis = 0 compareParticles = 0 -analysisRoutine = Examples/Tests/particles_in_pml/analysis_particles_in_pml_2dmr.py +analysisRoutine = Examples/Tests/particles_in_pml/analysis_particles_in_pml.py [particles_in_pml_3d_MR] buildDir = . From 41be7d15d3298f3534ef09433476c720bd9f0408 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Thu, 29 Feb 2024 02:45:48 -0800 Subject: [PATCH 02/16] Fix: Windows add_dll_directory Expand (#4734) Expand paths in `PATH` environment variable before adding them and make them absolute paths. --- Python/pywarpx/__init__.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/Python/pywarpx/__init__.py b/Python/pywarpx/__init__.py index f89926399ca..d5181a0307c 100644 --- a/Python/pywarpx/__init__.py +++ b/Python/pywarpx/__init__.py @@ -19,8 +19,9 @@ # add anything in PATH paths = os.environ.get("PATH", "") for p in paths.split(";"): - if os.path.exists(p): - os.add_dll_directory(p) + p_abs = os.path.abspath(os.path.expanduser(os.path.expandvars(p))) + if os.path.exists(p_abs): + os.add_dll_directory(p_abs) from .Algo import algo from .Amr import amr From 569aa1c4fb24b68d9f0abb1b622c472ad1099677 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Fri, 1 Mar 2024 13:28:51 -0800 Subject: [PATCH 03/16] Doc: Python Bld in Separate Dir (#4739) Python builds require shared AMReX and WarpX libraries. To avoid confusion with users that also build the executables, but do not use the full installer logic, we recommend to use two separate build directories, so the executables can use a static build of their dependendencies to stay relocatable (as in: can be copied without setting environment hints to dependent shared libs). --- Docs/source/install/cmake.rst | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Docs/source/install/cmake.rst b/Docs/source/install/cmake.rst index 0882efd7fe2..4474b3509fe 100644 --- a/Docs/source/install/cmake.rst +++ b/Docs/source/install/cmake.rst @@ -223,11 +223,11 @@ For PICMI Python bindings, configure WarpX to produce a library and call our ``p .. code-block:: bash # find dependencies & configure for all WarpX dimensionalities - cmake -S . -B build -DWarpX_DIMS="1;2;RZ;3" -DWarpX_PYTHON=ON + cmake -S . -B build_py -DWarpX_DIMS="1;2;RZ;3" -DWarpX_PYTHON=ON # build and then call "python3 -m pip install ..." - cmake --build build --target pip_install -j 4 + cmake --build build_py --target pip_install -j 4 **That's it!** You can now :ref:`run a first 3D PICMI script ` from our :ref:`examples `. @@ -324,10 +324,10 @@ This is the workflow most developers will prefer as it allows rapid re-compiles: .. code-block:: bash # build WarpX executables and libraries - cmake -S . -B build -DWarpX_DIMS="1;2;RZ;3" -DWarpX_PYTHON=ON + cmake -S . -B build_py -DWarpX_DIMS="1;2;RZ;3" -DWarpX_PYTHON=ON # build & install Python only - cmake --build build -j 4 --target pip_install + cmake --build build_py -j 4 --target pip_install There is also a ``--target pip_install_nodeps`` option that :ref:`skips pip-based dependency checks `. From b0ae812c1aa09b954b9c2fa5ba8b0f22a652696f Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Sat, 2 Mar 2024 02:15:46 -0800 Subject: [PATCH 04/16] Doc: Pre-Commit Locally (#4741) How to do automatic pre-commit tests locally. --- Docs/source/developers/testing.rst | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/Docs/source/developers/testing.rst b/Docs/source/developers/testing.rst index 2b3fe989f2f..6dcc1886e12 100644 --- a/Docs/source/developers/testing.rst +++ b/Docs/source/developers/testing.rst @@ -34,6 +34,20 @@ For instance, compiling with ``clang++ -Werror`` would be: export CXXFLAGS="-Werror" +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: + +.. code-block:: sh + + python -m pip install -U pre-commit + pre-commit install + +See `pre-commit.com `__ and our ``.pre-commit-config.yaml`` file in the repository for more details. + + Run the test suite locally -------------------------- From df312be5f371fb8b86a59a51a25d5d2d78ae09e6 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Mon, 4 Mar 2024 03:48:35 -0800 Subject: [PATCH 05/16] AMReX/pyAMReX/PICSAR: Weekly Update (#4745) * AMReX: Weekly Update * pyAMReX: Weekly Update --- .github/workflows/cuda.yml | 2 +- Regression/WarpX-GPU-tests.ini | 2 +- Regression/WarpX-tests.ini | 2 +- cmake/dependencies/AMReX.cmake | 2 +- cmake/dependencies/pyAMReX.cmake | 2 +- run_test.sh | 2 +- 6 files changed, 6 insertions(+), 6 deletions(-) diff --git a/.github/workflows/cuda.yml b/.github/workflows/cuda.yml index 368221efa9c..4b23ef931db 100644 --- a/.github/workflows/cuda.yml +++ b/.github/workflows/cuda.yml @@ -115,7 +115,7 @@ jobs: which nvcc || echo "nvcc not in PATH!" git clone https://github.com/AMReX-Codes/amrex.git ../amrex - cd ../amrex && git checkout --detach 2ecafcff40132f56eb2b494e1a374684ff97117a && cd - + cd ../amrex && git checkout --detach 3525b4a3f27eb64f746dd69b6613f71bb02d6e63 && cd - make COMP=gcc QED=FALSE USE_MPI=TRUE USE_GPU=TRUE USE_OMP=FALSE USE_PSATD=TRUE USE_CCACHE=TRUE -j 4 ccache -s diff --git a/Regression/WarpX-GPU-tests.ini b/Regression/WarpX-GPU-tests.ini index ef25ee0dfbe..40338dc1480 100644 --- a/Regression/WarpX-GPU-tests.ini +++ b/Regression/WarpX-GPU-tests.ini @@ -60,7 +60,7 @@ emailBody = Check https://ccse.lbl.gov/pub/GpuRegressionTesting/WarpX/ for more [AMReX] dir = /home/regtester/git/amrex/ -branch = 2ecafcff40132f56eb2b494e1a374684ff97117a +branch = 3525b4a3f27eb64f746dd69b6613f71bb02d6e63 [source] dir = /home/regtester/git/WarpX diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 128495195da..426ecb0717c 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -59,7 +59,7 @@ emailBody = Check https://ccse.lbl.gov/pub/RegressionTesting/WarpX/ for more det [AMReX] dir = /home/regtester/AMReX_RegTesting/amrex/ -branch = 2ecafcff40132f56eb2b494e1a374684ff97117a +branch = 3525b4a3f27eb64f746dd69b6613f71bb02d6e63 [source] dir = /home/regtester/AMReX_RegTesting/warpx diff --git a/cmake/dependencies/AMReX.cmake b/cmake/dependencies/AMReX.cmake index 9657cec146d..5b225de6d50 100644 --- a/cmake/dependencies/AMReX.cmake +++ b/cmake/dependencies/AMReX.cmake @@ -273,7 +273,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 "2ecafcff40132f56eb2b494e1a374684ff97117a" +set(WarpX_amrex_branch "3525b4a3f27eb64f746dd69b6613f71bb02d6e63" 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 71831357e66..e8748f2e920 100644 --- a/cmake/dependencies/pyAMReX.cmake +++ b/cmake/dependencies/pyAMReX.cmake @@ -79,7 +79,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 "0cbf4b08c9045e1845595c836b99f94bb3c1ac9f" +set(WarpX_pyamrex_branch "bfb599fd8361f8ef0765c487fd7bb69409bf78af" CACHE STRING "Repository branch for WarpX_pyamrex_repo if(WarpX_pyamrex_internal)") diff --git a/run_test.sh b/run_test.sh index ff41399a9c8..b5cdc627e35 100755 --- a/run_test.sh +++ b/run_test.sh @@ -68,7 +68,7 @@ python3 -m pip install --upgrade -r warpx/Regression/requirements.txt # Clone AMReX and warpx-data git clone https://github.com/AMReX-Codes/amrex.git -cd amrex && git checkout --detach 2ecafcff40132f56eb2b494e1a374684ff97117a && cd - +cd amrex && git checkout --detach 3525b4a3f27eb64f746dd69b6613f71bb02d6e63 && cd - # warpx-data contains various required data sets git clone --depth 1 https://github.com/ECP-WarpX/warpx-data.git # openPMD-example-datasets contains various required data sets From c4b1266571857b1d4d910848fc7c48e7a8c8dac8 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Mon, 4 Mar 2024 12:50:40 +0100 Subject: [PATCH 06/16] Clang-tidy CI tests: use clang-tidy-15 instead of clang-tidy-14 (#4689) * use clang-tidy-15 instead of clang-tidy-14 for CI tests * fix bug * address issues found with clang-tidy-15 * address more issues found with clang-tidy * address more issues found with clang-tidy * address additional issues found with clang-tidy * fix more issues found with clang-tidy * fixed bugs * fixed bugs * fixed bugs * fixed issue on Windows * fixed bug * address issues found with clang-tidy * make std::string filepath const both for Windows and other OS * make edge_lengths a const both with and without EB * address more issues found with clang-tidy * address issues found with clang-tidy * address the residual issues found with clang-tidy * addressed const-correctness issues * fixed bug * fixed bug * address issue found with clang-tidy * address issue found with clang-tidy --- .github/workflows/clang_tidy.yml | 8 +- .../dependencies/{clang14.sh => clang15.sh} | 8 +- .../AcceleratorLattice/AcceleratorLattice.cpp | 2 +- Source/BoundaryConditions/WarpX_PEC.cpp | 4 +- Source/Diagnostics/BTDiagnostics.cpp | 6 +- .../BackTransformFunctor.cpp | 2 +- .../FlushFormats/FlushFormatPlotfile.cpp | 2 +- Source/Diagnostics/FullDiagnostics.cpp | 4 +- .../ReducedDiags/ColliderRelevant.cpp | 6 +- .../ReducedDiags/FieldReduction.cpp | 2 +- .../ReducedDiags/ParticleHistogram2D.cpp | 6 +- .../Diagnostics/ReducedDiags/ReducedDiags.cpp | 4 +- Source/Diagnostics/WarpXOpenPMD.cpp | 2 +- .../EmbeddedBoundary/WarpXFaceExtensions.cpp | 2 +- Source/Evolve/WarpXOneStepImplicitPicard.cpp | 34 +- .../FiniteDifferenceSolver/ComputeDivE.cpp | 2 +- .../FiniteDifferenceSolver/EvolveE.cpp | 2 +- .../FiniteDifferenceSolver/EvolveF.cpp | 2 +- .../HybridPICModel/HybridPICModel.cpp | 26 +- .../HybridPICSolveE.cpp | 36 +- .../PsatdAlgorithmGalileanRZ.cpp | 10 +- .../SpectralAlgorithms/PsatdAlgorithmPml.cpp | 6 +- .../SpectralAlgorithms/PsatdAlgorithmRZ.cpp | 8 +- .../SpectralBaseAlgorithmRZ.cpp | 2 +- .../SpectralSolver/SpectralFieldDataRZ.cpp | 8 +- .../HankelTransform.cpp | 2 +- .../FieldSolver/WarpXPushFieldsHybridPIC.cpp | 2 +- Source/Fluids/MusclHancockUtils.H | 20 +- Source/Fluids/WarpXFluidContainer.cpp | 411 +++++++++--------- Source/Initialization/ExternalField.cpp | 2 +- Source/Initialization/InjectorPosition.H | 6 +- Source/Initialization/PlasmaInjector.cpp | 6 +- Source/Initialization/WarpXInitData.cpp | 4 +- .../BinaryCollision/BinaryCollision.H | 2 +- .../DSMC/CollisionFilterFunc.H | 2 +- .../Collision/BinaryCollision/DSMC/DSMC.cpp | 18 +- .../BinaryCollision/ParticleCreationFunc.H | 4 +- .../Particles/Deposition/CurrentDeposition.H | 34 +- Source/Particles/Filter/FilterFunctors.H | 6 +- Source/Particles/Gather/FieldGather.H | 80 ++-- Source/Particles/Gather/GetExternalFields.H | 2 +- Source/Particles/LaserParticleContainer.cpp | 8 +- .../ParticleCreation/DefaultInitialization.H | 10 +- .../Particles/PhysicalParticleContainer.cpp | 150 +++---- .../RigidInjectedParticleContainer.cpp | 2 +- Source/Utils/SpeciesUtils.cpp | 4 +- Source/Utils/WarpXMovingWindow.cpp | 2 +- Source/Utils/WarpXTagging.cpp | 4 +- Source/WarpX.cpp | 6 +- Source/ablastr/fields/PoissonSolver.H | 5 + Source/ablastr/utils/SignalHandling.cpp | 2 +- Source/ablastr/utils/msg_logger/MsgLogger.cpp | 3 + 52 files changed, 502 insertions(+), 489 deletions(-) rename .github/workflows/dependencies/{clang14.sh => clang15.sh} (93%) diff --git a/.github/workflows/clang_tidy.yml b/.github/workflows/clang_tidy.yml index 6a1172802a8..c166941f6e2 100644 --- a/.github/workflows/clang_tidy.yml +++ b/.github/workflows/clang_tidy.yml @@ -15,7 +15,7 @@ jobs: - uses: actions/checkout@v4 - name: install dependencies run: | - .github/workflows/dependencies/clang14.sh + .github/workflows/dependencies/clang15.sh - name: set up cache uses: actions/cache@v4 with: @@ -32,8 +32,8 @@ jobs: export CCACHE_LOGFILE=${{ github.workspace }}/ccache.log.txt ccache -z - export CXX=$(which clang++) - export CC=$(which clang) + export CXX=$(which clang++-15) + export CC=$(which clang-15) cmake -S . -B build_clang_tidy \ -DCMAKE_VERBOSE_MAKEFILE=ON \ @@ -51,7 +51,7 @@ jobs: ${{github.workspace}}/.github/workflows/source/makeMakefileForClangTidy.py --input ${{github.workspace}}/ccache.log.txt make -j4 --keep-going -f clang-tidy-ccache-misses.mak \ - CLANG_TIDY=clang-tidy \ + CLANG_TIDY=clang-tidy-15 \ CLANG_TIDY_ARGS="--config-file=${{github.workspace}}/.clang-tidy --warnings-as-errors=*" ccache -s diff --git a/.github/workflows/dependencies/clang14.sh b/.github/workflows/dependencies/clang15.sh similarity index 93% rename from .github/workflows/dependencies/clang14.sh rename to .github/workflows/dependencies/clang15.sh index 27597e06c5d..63d5d70956f 100755 --- a/.github/workflows/dependencies/clang14.sh +++ b/.github/workflows/dependencies/clang15.sh @@ -15,17 +15,17 @@ echo 'Acquire::Retries "3";' | sudo tee /etc/apt/apt.conf.d/80-retries sudo apt-get -qqq update sudo apt-get install -y \ cmake \ - clang-14 \ - clang-tidy-14 \ + clang-15 \ + clang-tidy-15 \ libblas-dev \ - libc++-14-dev \ + libc++-15-dev \ libboost-math-dev \ libfftw3-dev \ libfftw3-mpi-dev \ libhdf5-openmpi-dev \ liblapack-dev \ libopenmpi-dev \ - libomp-dev \ + libomp-15-dev \ ninja-build # ccache diff --git a/Source/AcceleratorLattice/AcceleratorLattice.cpp b/Source/AcceleratorLattice/AcceleratorLattice.cpp index 8dc0983d622..edccae9374a 100644 --- a/Source/AcceleratorLattice/AcceleratorLattice.cpp +++ b/Source/AcceleratorLattice/AcceleratorLattice.cpp @@ -102,6 +102,6 @@ AcceleratorLattice::UpdateElementFinder (int const lev) // NOLINT(readability-ma LatticeElementFinderDevice AcceleratorLattice::GetFinderDeviceInstance (WarpXParIter const& a_pti, int const a_offset) const { - LatticeElementFinder & finder = (*m_element_finder)[a_pti]; + const LatticeElementFinder & finder = (*m_element_finder)[a_pti]; return finder.GetFinderDeviceInstance(a_pti, a_offset, *this); } diff --git a/Source/BoundaryConditions/WarpX_PEC.cpp b/Source/BoundaryConditions/WarpX_PEC.cpp index 4692d8e980f..7555e2aaef2 100644 --- a/Source/BoundaryConditions/WarpX_PEC.cpp +++ b/Source/BoundaryConditions/WarpX_PEC.cpp @@ -500,7 +500,7 @@ PEC::ApplyPECtoElectronPressure (amrex::MultiFab* Pefield, const int lev, amrex::Box domain_box = warpx.Geom(lev).Domain(); if (patch_type == PatchType::coarse) { - amrex::IntVect ref_ratio = ( (lev > 0) ? WarpX::RefRatio(lev-1) : amrex::IntVect(1) ); + const amrex::IntVect ref_ratio = ( (lev > 0) ? WarpX::RefRatio(lev-1) : amrex::IntVect(1) ); domain_box.coarsen(ref_ratio); } domain_box.convert(Pefield->ixType()); @@ -552,7 +552,7 @@ PEC::ApplyPECtoElectronPressure (amrex::MultiFab* Pefield, const int lev, amrex::ignore_unused(j,k); #endif // Store the array index - amrex::IntVect iv(AMREX_D_DECL(i,j,k)); + const amrex::IntVect iv(AMREX_D_DECL(i,j,k)); PEC::SetNeumannOnPEC(n, iv, Pe_array, mirrorfac, is_pec, fabbox); }); diff --git a/Source/Diagnostics/BTDiagnostics.cpp b/Source/Diagnostics/BTDiagnostics.cpp index 41c710805e9..fd9c5c783c1 100644 --- a/Source/Diagnostics/BTDiagnostics.cpp +++ b/Source/Diagnostics/BTDiagnostics.cpp @@ -901,7 +901,7 @@ BTDiagnostics::DefineFieldBufferMultiFab (const int i_buffer, const int lev) const int hi_k_lab = m_buffer_k_index_hi[i_buffer]; m_buffer_box[i_buffer].setSmall( m_moving_window_dir, hi_k_lab - m_buffer_size + 1); m_buffer_box[i_buffer].setBig( m_moving_window_dir, hi_k_lab ); - amrex::BoxArray buffer_ba( m_buffer_box[i_buffer] ); + const amrex::BoxArray buffer_ba( m_buffer_box[i_buffer] ); // Generate a new distribution map for the back-transformed buffer multifab const amrex::DistributionMapping buffer_dmap(buffer_ba); // Number of guard cells for the output buffer is zero. @@ -1040,7 +1040,7 @@ BTDiagnostics::Flush (int i_buffer, bool force_flush) m_buffer_box[i_buffer].setSmall(m_moving_window_dir, (m_buffer_box[i_buffer].smallEnd(m_moving_window_dir) - 1) ); m_buffer_box[i_buffer].setBig(m_moving_window_dir, (m_buffer_box[i_buffer].bigEnd(m_moving_window_dir) + 1) ); const amrex::Box particle_buffer_box = m_buffer_box[i_buffer]; - amrex::BoxArray buffer_ba( particle_buffer_box ); + const amrex::BoxArray buffer_ba( particle_buffer_box ); m_particles_buffer[i_buffer][0]->SetParticleBoxArray(0, buffer_ba); for (int isp = 0; isp < m_particles_buffer.at(i_buffer).size(); ++isp) { // BTD output is single level. Setting particle geometry, dmap, boxarray to level0 @@ -1467,7 +1467,7 @@ BTDiagnostics::PrepareParticleDataForOutput() DefineFieldBufferMultiFab(i_buffer, lev); } const amrex::Box particle_buffer_box = m_buffer_box[i_buffer]; - amrex::BoxArray buffer_ba( particle_buffer_box ); + const amrex::BoxArray buffer_ba( particle_buffer_box ); const amrex::DistributionMapping buffer_dmap(buffer_ba); m_particles_buffer[i_buffer][i]->SetParticleBoxArray(lev, buffer_ba); m_particles_buffer[i_buffer][i]->SetParticleDistributionMap(lev, buffer_dmap); diff --git a/Source/Diagnostics/ComputeDiagFunctors/BackTransformFunctor.cpp b/Source/Diagnostics/ComputeDiagFunctors/BackTransformFunctor.cpp index 22754f51a53..c4809df11a3 100644 --- a/Source/Diagnostics/ComputeDiagFunctors/BackTransformFunctor.cpp +++ b/Source/Diagnostics/ComputeDiagFunctors/BackTransformFunctor.cpp @@ -81,7 +81,7 @@ BackTransformFunctor::operator ()(amrex::MultiFab& mf_dst, int /*dcomp*/, const slice_box.setBig(moving_window_dir, i_boost); // Make it a BoxArray - amrex::BoxArray slice_ba(slice_box); + const amrex::BoxArray slice_ba(slice_box); // Define MultiFab with the distribution map of the destination multifab and // containing all ten components that were in the slice generated from m_mf_src. std::unique_ptr< amrex::MultiFab > tmp_slice_ptr = nullptr; diff --git a/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp b/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp index 880e2df01ff..48361d90fc0 100644 --- a/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp +++ b/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp @@ -303,7 +303,7 @@ FlushFormatPlotfile::WriteWarpXHeader( warpx.GetPartContainer().WriteHeader(HeaderFile); - MultiParticleContainer& mypc = warpx.GetPartContainer(); + const MultiParticleContainer& mypc = warpx.GetPartContainer(); const int n_species = mypc.nSpecies(); for (int i=0; i backward_strings; WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !pp_rd_name.queryarr("reduced_function(x,y,z,Ex,Ey,Ez,Bx,By,Bz)", backward_strings), diff --git a/Source/Diagnostics/ReducedDiags/ParticleHistogram2D.cpp b/Source/Diagnostics/ReducedDiags/ParticleHistogram2D.cpp index 473a53715e6..bcc22eac598 100644 --- a/Source/Diagnostics/ReducedDiags/ParticleHistogram2D.cpp +++ b/Source/Diagnostics/ReducedDiags/ParticleHistogram2D.cpp @@ -267,10 +267,12 @@ void ParticleHistogram2D::WriteToFile (int step) const const std::string fileSuffix = std::string("_%0") + std::to_string(m_file_min_digits) + std::string("T"); filename = filename.append(fileSuffix).append(".").append(m_openpmd_backend); - std::string filepath = m_path + m_rd_name + "/" + filename; // transform paths for Windows #ifdef _WIN32 - filepath = openPMD::auxiliary::replace_all(filepath, "/", "\\"); + const std::string filepath = openPMD::auxiliary::replace_all( + m_path + m_rd_name + "/" + filename, "/", "\\"); + #else + const std::string filepath = m_path + m_rd_name + "/" + filename; #endif // Create the OpenPMD series diff --git a/Source/Diagnostics/ReducedDiags/ReducedDiags.cpp b/Source/Diagnostics/ReducedDiags/ReducedDiags.cpp index ae6dae25a71..483f4c68ece 100644 --- a/Source/Diagnostics/ReducedDiags/ReducedDiags.cpp +++ b/Source/Diagnostics/ReducedDiags/ReducedDiags.cpp @@ -40,7 +40,7 @@ ReducedDiags::ReducedDiags (std::string rd_name) std::string restart_chkfile; const ParmParse pp_amr("amr"); pp_amr.query("restart", restart_chkfile); - bool IsNotRestart = restart_chkfile.empty(); + const bool IsNotRestart = restart_chkfile.empty(); if (ParallelDescriptor::IOProcessor()) { @@ -50,7 +50,7 @@ ReducedDiags::ReducedDiags (std::string rd_name) { amrex::CreateDirectoryFailed(m_path); } // replace / create output file - std::string rd_full_file_name = m_path + m_rd_name + "." + m_extension; + const std::string rd_full_file_name = m_path + m_rd_name + "." + m_extension; m_write_header = IsNotRestart || !amrex::FileExists(rd_full_file_name); // not a restart or file doesn't exist if (m_write_header) { diff --git a/Source/Diagnostics/WarpXOpenPMD.cpp b/Source/Diagnostics/WarpXOpenPMD.cpp index ec8daa52dd9..b55dbe79175 100644 --- a/Source/Diagnostics/WarpXOpenPMD.cpp +++ b/Source/Diagnostics/WarpXOpenPMD.cpp @@ -673,7 +673,7 @@ WarpXOpenPMDPlot::DumpToFile (ParticleContainer* pc, openPMD::ParticleSpecies currSpecies = currIteration.particles[name]; // only BTD writes multiple times into the same step, zero for other methods - unsigned long ParticleFlushOffset = isBTD ? num_already_flushed(currSpecies) : 0; + const unsigned long ParticleFlushOffset = isBTD ? num_already_flushed(currSpecies) : 0; // prepare data structures the first time BTD has non-zero particles // we set some of them to zero extent, so we need to time that well diff --git a/Source/EmbeddedBoundary/WarpXFaceExtensions.cpp b/Source/EmbeddedBoundary/WarpXFaceExtensions.cpp index c885e18b672..21c13f23845 100644 --- a/Source/EmbeddedBoundary/WarpXFaceExtensions.cpp +++ b/Source/EmbeddedBoundary/WarpXFaceExtensions.cpp @@ -333,7 +333,7 @@ ComputeNBorrowOneFaceExtension(const amrex::Dim3 cell, const amrex::Real S_ext, for (int i_n = -1; i_n < 2; i_n++) { for (int j_n = -1; j_n < 2; j_n++) { //This if makes sure that we don't visit the "diagonal neighbours" - if (!(i_n == j_n || i_n == -j_n)) { + if ((i_n != j_n) && (i_n != -j_n)) { // Here a face is available if it doesn't need to be extended itself and if its // area exceeds Sz_ext. Here we need to take into account if the intruded face // has given away already some area, so we use Sz_red rather than Sz. diff --git a/Source/Evolve/WarpXOneStepImplicitPicard.cpp b/Source/Evolve/WarpXOneStepImplicitPicard.cpp index 56ed26db1e2..22d42dfe89d 100644 --- a/Source/Evolve/WarpXOneStepImplicitPicard.cpp +++ b/Source/Evolve/WarpXOneStepImplicitPicard.cpp @@ -150,8 +150,8 @@ WarpX::OneStep_ImplicitPicard(amrex::Real cur_time) // particle velocities by dt, then take average of old and new v, // deposit currents, giving J at n+1/2 // This uses Efield_fp and Bfield_fp, the field at n+1/2 from the previous iteration. - bool skip_current = false; - PushType push_type = PushType::Implicit; + const bool skip_current = false; + const PushType push_type = PushType::Implicit; PushParticlesandDeposit(cur_time, skip_current, push_type); SyncCurrentAndRho(); @@ -209,23 +209,23 @@ WarpX::OneStep_ImplicitPicard(amrex::Real cur_time) Efield_save[0][0]->minus(*Efield_fp[0][0], 0, ncomps, 0); Efield_save[0][1]->minus(*Efield_fp[0][1], 0, ncomps, 0); Efield_save[0][2]->minus(*Efield_fp[0][2], 0, ncomps, 0); - amrex::Real maxE0 = std::max(1._rt, Efield_fp[0][0]->norm0(0, 0)); - amrex::Real maxE1 = std::max(1._rt, Efield_fp[0][1]->norm0(0, 0)); - amrex::Real maxE2 = std::max(1._rt, Efield_fp[0][2]->norm0(0, 0)); - amrex::Real deltaE0 = Efield_save[0][0]->norm0(0, 0)/maxE0; - amrex::Real deltaE1 = Efield_save[0][1]->norm0(0, 0)/maxE1; - amrex::Real deltaE2 = Efield_save[0][2]->norm0(0, 0)/maxE2; + const amrex::Real maxE0 = std::max(1._rt, Efield_fp[0][0]->norm0(0, 0)); + const amrex::Real maxE1 = std::max(1._rt, Efield_fp[0][1]->norm0(0, 0)); + const amrex::Real maxE2 = std::max(1._rt, Efield_fp[0][2]->norm0(0, 0)); + const amrex::Real deltaE0 = Efield_save[0][0]->norm0(0, 0)/maxE0; + const amrex::Real deltaE1 = Efield_save[0][1]->norm0(0, 0)/maxE1; + const amrex::Real deltaE2 = Efield_save[0][2]->norm0(0, 0)/maxE2; deltaE = std::max(std::max(deltaE0, deltaE1), deltaE2); if (evolve_scheme == EvolveScheme::ImplicitPicard) { Bfield_save[0][0]->minus(*Bfield_fp[0][0], 0, ncomps, 0); Bfield_save[0][1]->minus(*Bfield_fp[0][1], 0, ncomps, 0); Bfield_save[0][2]->minus(*Bfield_fp[0][2], 0, ncomps, 0); - amrex::Real maxB0 = std::max(1._rt, Bfield_fp[0][0]->norm0(0, 0)); - amrex::Real maxB1 = std::max(1._rt, Bfield_fp[0][1]->norm0(0, 0)); - amrex::Real maxB2 = std::max(1._rt, Bfield_fp[0][2]->norm0(0, 0)); - amrex::Real deltaB0 = Bfield_save[0][0]->norm0(0, 0)/maxB0; - amrex::Real deltaB1 = Bfield_save[0][1]->norm0(0, 0)/maxB1; - amrex::Real deltaB2 = Bfield_save[0][2]->norm0(0, 0)/maxB2; + const amrex::Real maxB0 = std::max(1._rt, Bfield_fp[0][0]->norm0(0, 0)); + const amrex::Real maxB1 = std::max(1._rt, Bfield_fp[0][1]->norm0(0, 0)); + const amrex::Real maxB2 = std::max(1._rt, Bfield_fp[0][2]->norm0(0, 0)); + const amrex::Real deltaB0 = Bfield_save[0][0]->norm0(0, 0)/maxB0; + const amrex::Real deltaB1 = Bfield_save[0][1]->norm0(0, 0)/maxB1; + const amrex::Real deltaB2 = Bfield_save[0][2]->norm0(0, 0)/maxB2; deltaB = std::max(std::max(deltaB0, deltaB1), deltaB2); } else { deltaB = 0.; @@ -401,9 +401,9 @@ WarpX::FinishImplicitFieldUpdate(amrex::Vector const& Fy_n = Field_n[lev][1]->array(mfi); amrex::Array4 const& Fz_n = Field_n[lev][2]->array(mfi); - amrex::Box tbx = mfi.tilebox(Field_fp[lev][0]->ixType().toIntVect()); - amrex::Box tby = mfi.tilebox(Field_fp[lev][1]->ixType().toIntVect()); - amrex::Box tbz = mfi.tilebox(Field_fp[lev][2]->ixType().toIntVect()); + amrex::Box const tbx = mfi.tilebox(Field_fp[lev][0]->ixType().toIntVect()); + amrex::Box const tby = mfi.tilebox(Field_fp[lev][1]->ixType().toIntVect()); + amrex::Box const tbz = mfi.tilebox(Field_fp[lev][2]->ixType().toIntVect()); amrex::ParallelFor( tbx, ncomps, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/ComputeDivE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/ComputeDivE.cpp index 3268c81ac81..3cc05c58068 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/ComputeDivE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/ComputeDivE.cpp @@ -133,7 +133,7 @@ void FiniteDifferenceSolver::ComputeDivECylindrical ( for ( MFIter mfi(divEfield, TilingIfNotGPU()); mfi.isValid(); ++mfi ) { // Extract field data for this grid/tile - Array4 divE = divEfield.array(mfi); + const Array4 divE = divEfield.array(mfi); Array4 const& Er = Efield[0]->array(mfi); Array4 const& Et = Efield[1]->array(mfi); Array4 const& Ez = Efield[2]->array(mfi); diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp index 74922650a42..e52cca846bb 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp @@ -378,7 +378,7 @@ void FiniteDifferenceSolver::EvolveECylindrical ( if (Ffield) { // Extract field data for this grid/tile - Array4 F = Ffield->array(mfi); + const Array4 F = Ffield->array(mfi); // Loop over the cells and update the fields amrex::ParallelFor(ter, tet, tez, diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveF.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveF.cpp index 83014a08026..8ce578bb52a 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveF.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveF.cpp @@ -148,7 +148,7 @@ void FiniteDifferenceSolver::EvolveFCylindrical ( for ( MFIter mfi(*Ffield, TilingIfNotGPU()); mfi.isValid(); ++mfi ) { // Extract field data for this grid/tile - Array4 F = Ffield->array(mfi); + const Array4 F = Ffield->array(mfi); Array4 const& Er = Efield[0]->array(mfi); Array4 const& Et = Efield[1]->array(mfi); Array4 const& Ez = Efield[2]->array(mfi); diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp index 034bb71efbc..8c2b5a0707a 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp @@ -19,7 +19,7 @@ HybridPICModel::HybridPICModel ( int nlevs_max ) void HybridPICModel::ReadParameters () { - ParmParse pp_hybrid("hybrid_pic_model"); + const ParmParse pp_hybrid("hybrid_pic_model"); // The B-field update is subcycled to improve stability - the number // of sub steps can be specified by the user (defaults to 50). @@ -32,7 +32,7 @@ void HybridPICModel::ReadParameters () if (!utils::parser::queryWithParser(pp_hybrid, "elec_temp", m_elec_temp)) { Abort("hybrid_pic_model.elec_temp must be specified when using the hybrid solver"); } - bool n0_ref_given = utils::parser::queryWithParser(pp_hybrid, "n0_ref", m_n0_ref); + const bool n0_ref_given = utils::parser::queryWithParser(pp_hybrid, "n0_ref", m_n0_ref); if (m_gamma != 1.0 && !n0_ref_given) { Abort("hybrid_pic_model.n0_ref should be specified if hybrid_pic_model.gamma != 1"); } @@ -219,23 +219,23 @@ void HybridPICModel::InitData () // Initialize external current - note that this approach skips the check // if the current is time dependent which is what needs to be done to // write time independent fields on the first step. - std::array< std::unique_ptr, 3 > edge_lengths; - for (int lev = 0; lev <= warpx.finestLevel(); ++lev) { #ifdef AMREX_USE_EB auto& edge_lengths_x = warpx.getedgelengths(lev, 0); - edge_lengths[0] = std::make_unique( - edge_lengths_x, amrex::make_alias, 0, edge_lengths_x.nComp() - ); auto& edge_lengths_y = warpx.getedgelengths(lev, 1); - edge_lengths[1] = std::make_unique( - edge_lengths_y, amrex::make_alias, 0, edge_lengths_y.nComp() - ); auto& edge_lengths_z = warpx.getedgelengths(lev, 2); - edge_lengths[2] = std::make_unique( - edge_lengths_z, amrex::make_alias, 0, edge_lengths_z.nComp() - ); + + const auto edge_lengths = std::array< std::unique_ptr, 3 >{ + std::make_unique( + edge_lengths_x, amrex::make_alias, 0, edge_lengths_x.nComp()), + std::make_unique( + edge_lengths_y, amrex::make_alias, 0, edge_lengths_y.nComp()), + std::make_unique( + edge_lengths_z, amrex::make_alias, 0, edge_lengths_z.nComp()) + }; +#else + const auto edge_lengths = std::array< std::unique_ptr, 3 >(); #endif GetCurrentExternal(edge_lengths, lev); } diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp index 5100eed0df3..066e88b13c0 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp @@ -593,9 +593,9 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical ( // Interpolate current to appropriate staggering to match E field Real jtot_val = 0._rt; if (include_resistivity_term && resistivity_has_J_dependence) { - Real jr_val = Interp(Jr, Jr_stag, Er_stag, coarsen, i, j, 0, 0); - Real jt_val = Interp(Jt, Jt_stag, Er_stag, coarsen, i, j, 0, 0); - Real jz_val = Interp(Jz, Jz_stag, Er_stag, coarsen, i, j, 0, 0); + const Real jr_val = Interp(Jr, Jr_stag, Er_stag, coarsen, i, j, 0, 0); + const Real jt_val = Interp(Jt, Jt_stag, Er_stag, coarsen, i, j, 0, 0); + const Real jz_val = Interp(Jz, Jz_stag, Er_stag, coarsen, i, j, 0, 0); jtot_val = std::sqrt(jr_val*jr_val + jt_val*jt_val + jz_val*jz_val); } @@ -635,9 +635,9 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical ( // Interpolate current to appropriate staggering to match E field Real jtot_val = 0._rt; if (include_resistivity_term && resistivity_has_J_dependence) { - Real jr_val = Interp(Jr, Jr_stag, Et_stag, coarsen, i, j, 0, 0); - Real jt_val = Interp(Jt, Jt_stag, Et_stag, coarsen, i, j, 0, 0); - Real jz_val = Interp(Jz, Jz_stag, Et_stag, coarsen, i, j, 0, 0); + const Real jr_val = Interp(Jr, Jr_stag, Et_stag, coarsen, i, j, 0, 0); + const Real jt_val = Interp(Jt, Jt_stag, Et_stag, coarsen, i, j, 0, 0); + const Real jz_val = Interp(Jz, Jz_stag, Et_stag, coarsen, i, j, 0, 0); jtot_val = std::sqrt(jr_val*jr_val + jt_val*jt_val + jz_val*jz_val); } @@ -669,9 +669,9 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical ( // Interpolate current to appropriate staggering to match E field Real jtot_val = 0._rt; if (include_resistivity_term && resistivity_has_J_dependence) { - Real jr_val = Interp(Jr, Jr_stag, Ez_stag, coarsen, i, j, 0, 0); - Real jt_val = Interp(Jt, Jt_stag, Ez_stag, coarsen, i, j, 0, 0); - Real jz_val = Interp(Jz, Jz_stag, Ez_stag, coarsen, i, j, 0, 0); + const Real jr_val = Interp(Jr, Jr_stag, Ez_stag, coarsen, i, j, 0, 0); + const Real jt_val = Interp(Jt, Jt_stag, Ez_stag, coarsen, i, j, 0, 0); + const Real jz_val = Interp(Jz, Jz_stag, Ez_stag, coarsen, i, j, 0, 0); jtot_val = std::sqrt(jr_val*jr_val + jt_val*jt_val + jz_val*jz_val); } @@ -885,9 +885,9 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian ( // Interpolate current to appropriate staggering to match E field Real jtot_val = 0._rt; if (include_resistivity_term && resistivity_has_J_dependence) { - Real jx_val = Interp(Jx, Jx_stag, Ex_stag, coarsen, i, j, k, 0); - Real jy_val = Interp(Jy, Jy_stag, Ex_stag, coarsen, i, j, k, 0); - Real jz_val = Interp(Jz, Jz_stag, Ex_stag, coarsen, i, j, k, 0); + const Real jx_val = Interp(Jx, Jx_stag, Ex_stag, coarsen, i, j, k, 0); + const Real jy_val = Interp(Jy, Jy_stag, Ex_stag, coarsen, i, j, k, 0); + const Real jz_val = Interp(Jz, Jz_stag, Ex_stag, coarsen, i, j, k, 0); jtot_val = std::sqrt(jx_val*jx_val + jy_val*jy_val + jz_val*jz_val); } @@ -924,9 +924,9 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian ( // Interpolate current to appropriate staggering to match E field Real jtot_val = 0._rt; if (include_resistivity_term && resistivity_has_J_dependence) { - Real jx_val = Interp(Jx, Jx_stag, Ey_stag, coarsen, i, j, k, 0); - Real jy_val = Interp(Jy, Jy_stag, Ey_stag, coarsen, i, j, k, 0); - Real jz_val = Interp(Jz, Jz_stag, Ey_stag, coarsen, i, j, k, 0); + const Real jx_val = Interp(Jx, Jx_stag, Ey_stag, coarsen, i, j, k, 0); + const Real jy_val = Interp(Jy, Jy_stag, Ey_stag, coarsen, i, j, k, 0); + const Real jz_val = Interp(Jz, Jz_stag, Ey_stag, coarsen, i, j, k, 0); jtot_val = std::sqrt(jx_val*jx_val + jy_val*jy_val + jz_val*jz_val); } @@ -957,9 +957,9 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian ( // Interpolate current to appropriate staggering to match E field Real jtot_val = 0._rt; if (include_resistivity_term && resistivity_has_J_dependence) { - Real jx_val = Interp(Jx, Jx_stag, Ez_stag, coarsen, i, j, k, 0); - Real jy_val = Interp(Jy, Jy_stag, Ez_stag, coarsen, i, j, k, 0); - Real jz_val = Interp(Jz, Jz_stag, Ez_stag, coarsen, i, j, k, 0); + const Real jx_val = Interp(Jx, Jx_stag, Ez_stag, coarsen, i, j, k, 0); + const Real jy_val = Interp(Jy, Jy_stag, Ez_stag, coarsen, i, j, k, 0); + const Real jz_val = Interp(Jz, Jz_stag, Ez_stag, coarsen, i, j, k, 0); jtot_val = std::sqrt(jx_val*jx_val + jy_val*jy_val + jz_val*jz_val); } diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmGalileanRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmGalileanRZ.cpp index 71cee95524e..c24648fbf07 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmGalileanRZ.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmGalileanRZ.cpp @@ -196,7 +196,7 @@ void PsatdAlgorithmGalileanRZ::InitializeSpectralCoefficients (SpectralFieldData amrex::Array4 const& T_rho = T_rho_coef[mfi].array(); // Extract real (for portability on GPU) - amrex::Real vz = m_v_galilean[2]; + const amrex::Real vz = m_v_galilean[2]; auto const & kr_modes = f.getKrArray(mfi); amrex::Real const* kr_arr = kr_modes.dataPtr(); @@ -247,7 +247,7 @@ void PsatdAlgorithmGalileanRZ::InitializeSpectralCoefficients (SpectralFieldData // update equation have been modified accordingly so that // the expressions below (with the update equations) // are mathematically equivalent to those of the paper. - Complex x1 = 1._rt/(1._rt-nu*nu) * + const Complex x1 = 1._rt/(1._rt-nu*nu) * (theta_star - C(i,j,k,mode)*theta + I*kv*S_ck(i,j,k,mode)*theta); // x1, above, is identical to the original paper X1(i,j,k,mode) = theta*x1/(ep0*c*c*k_norm*k_norm); @@ -300,7 +300,7 @@ PsatdAlgorithmGalileanRZ::CurrentCorrection (SpectralFieldDataRZ& field_data) amrex::Box const & bx = field_data.fields[mfi].box(); // Extract arrays for the fields to be updated - amrex::Array4 fields = field_data.fields[mfi].array(); + const amrex::Array4 fields = field_data.fields[mfi].array(); // Extract pointers for the k vectors auto const & kr_modes = field_data.getKrArray(mfi); @@ -309,8 +309,8 @@ PsatdAlgorithmGalileanRZ::CurrentCorrection (SpectralFieldDataRZ& field_data) int const nr = bx.length(0); // Local copy of member variables before GPU loop - amrex::Real vz = m_v_galilean[2]; - amrex::Real const dt = m_dt; + const amrex::Real vz = m_v_galilean[2]; + const amrex::Real dt = m_dt; // Loop over indices within one box int const modes = field_data.n_rz_azimuthal_modes; diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPml.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPml.cpp index 92b8b3e900d..792d5126ef7 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPml.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPml.cpp @@ -410,11 +410,11 @@ void PsatdAlgorithmPml::InitializeSpectralCoefficients ( } // Extract Galilean velocity - amrex::Real vg_x = m_v_galilean[0]; + const amrex::Real vg_x = m_v_galilean[0]; #if defined(WARPX_DIM_3D) - amrex::Real vg_y = m_v_galilean[1]; + const amrex::Real vg_y = m_v_galilean[1]; #endif - amrex::Real vg_z = m_v_galilean[2]; + const amrex::Real vg_z = m_v_galilean[2]; // Loop over indices within one box ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp index fcb89c39109..f621d3c7af9 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp @@ -127,9 +127,9 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f) amrex::ParallelFor(bx, modes, [=] AMREX_GPU_DEVICE(int i, int j, int k, int mode) noexcept { - int idx_jx = (J_linear) ? static_cast(Idx.Jx_old) : static_cast(Idx.Jx_mid); - int idx_jy = (J_linear) ? static_cast(Idx.Jy_old) : static_cast(Idx.Jy_mid); - int idx_jz = (J_linear) ? static_cast(Idx.Jz_old) : static_cast(Idx.Jz_mid); + const int idx_jx = (J_linear) ? static_cast(Idx.Jx_old) : static_cast(Idx.Jx_mid); + const int idx_jy = (J_linear) ? static_cast(Idx.Jy_old) : static_cast(Idx.Jy_mid); + const int idx_jz = (J_linear) ? static_cast(Idx.Jz_old) : static_cast(Idx.Jz_mid); // All of the fields of each mode are grouped together int const Ep_m = Idx.Ex + Idx.n_fields*mode; @@ -432,7 +432,7 @@ PsatdAlgorithmRZ::CurrentCorrection (SpectralFieldDataRZ& field_data) amrex::Box const & bx = field_data.fields[mfi].box(); // Extract arrays for the fields to be updated - amrex::Array4 fields = field_data.fields[mfi].array(); + const amrex::Array4 fields = field_data.fields[mfi].array(); // Extract pointers for the k vectors auto const & kr_modes = field_data.getKrArray(mfi); diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.cpp index acbe3815cc1..f8ef0ef4730 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.cpp @@ -38,7 +38,7 @@ SpectralBaseAlgorithmRZ::ComputeSpectralDivE ( Box const & bx = field_data.fields[mfi].box(); // Extract arrays for the fields to be updated - Array4 fields = field_data.fields[mfi].array(); + const Array4 fields = field_data.fields[mfi].array(); // Extract pointers for the k vectors Real const * kr_arr = field_data.getKrArray(mfi).dataPtr(); diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.cpp index c4a693bf25a..b26248fb2d0 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.cpp @@ -662,7 +662,7 @@ SpectralFieldDataRZ::BackwardTransform (const int lev, amrex::MultiFab& field_mf_t, int const field_index_t) { amrex::LayoutData* cost = WarpX::getCosts(lev); - bool do_costs = WarpXUtilLoadBalance::doCosts(cost, field_mf_r.boxArray(), field_mf_r.DistributionMap()); + const bool do_costs = WarpXUtilLoadBalance::doCosts(cost, field_mf_r.boxArray(), field_mf_r.DistributionMap()); // Check field index type, in order to apply proper shift in spectral space. bool const is_nodal_z = field_mf_r.is_nodal(1); @@ -727,7 +727,7 @@ SpectralFieldDataRZ::BackwardTransform (const int lev, sign = -1._rt; } else { // Even modes are anti-symmetric - int imode = (icomp + 1)/2; + const int imode = (icomp + 1)/2; sign = static_cast(std::pow(-1._rt, imode+1)); } } @@ -773,7 +773,7 @@ void SpectralFieldDataRZ::ApplyFilter (const int lev, int const field_index) { amrex::LayoutData* cost = WarpX::getCosts(lev); - bool do_costs = WarpXUtilLoadBalance::doCosts(cost, binomialfilter.boxArray(), binomialfilter.DistributionMap()); + const bool do_costs = WarpXUtilLoadBalance::doCosts(cost, binomialfilter.boxArray(), binomialfilter.DistributionMap()); for (amrex::MFIter mfi(binomialfilter); mfi.isValid(); ++mfi){ @@ -818,7 +818,7 @@ SpectralFieldDataRZ::ApplyFilter (const int lev, int const field_index1, int const field_index2, int const field_index3) { amrex::LayoutData* cost = WarpX::getCosts(lev); - bool do_costs = WarpXUtilLoadBalance::doCosts(cost, binomialfilter.boxArray(), binomialfilter.DistributionMap()); + const bool do_costs = WarpXUtilLoadBalance::doCosts(cost, binomialfilter.boxArray(), binomialfilter.DistributionMap()); for (amrex::MFIter mfi(binomialfilter); mfi.isValid(); ++mfi){ diff --git a/Source/FieldSolver/SpectralSolver/SpectralHankelTransform/HankelTransform.cpp b/Source/FieldSolver/SpectralSolver/SpectralHankelTransform/HankelTransform.cpp index a8930e1d48b..9f4cac71736 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralHankelTransform/HankelTransform.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralHankelTransform/HankelTransform.cpp @@ -54,7 +54,7 @@ HankelTransform::HankelTransform (int const hankel_order, // Calculate the spatial grid (Uniform grid with a half-cell offset) amrex::Vector rmesh(m_nr); - amrex::Real dr = rmax/m_nr; + const amrex::Real dr = rmax/m_nr; for (int ir=0 ; ir < m_nr ; ir++) { rmesh[ir] = dr*(ir + 0.5_rt); } diff --git a/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp b/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp index 50a85a8cdee..ae10ae1e19a 100644 --- a/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp +++ b/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp @@ -54,7 +54,7 @@ void WarpX::HybridPICEvolveFields () } // Get requested number of substeps to use - int sub_steps = m_hybrid_pic_model->m_substeps; + const int sub_steps = m_hybrid_pic_model->m_substeps; // Get the external current m_hybrid_pic_model->GetCurrentExternal(m_edge_lengths); diff --git a/Source/Fluids/MusclHancockUtils.H b/Source/Fluids/MusclHancockUtils.H index 2765f28f3b3..d75e33211b2 100644 --- a/Source/Fluids/MusclHancockUtils.H +++ b/Source/Fluids/MusclHancockUtils.H @@ -36,7 +36,7 @@ amrex::Real V_calc (const amrex::Array4& U, int i, int j, int k, in { using namespace amrex::literals; // comp -> x, y, z -> 0, 1, 2, return Vx, Vy, or Vz: - amrex::Real gamma = std::sqrt(1.0_rt + (U(i,j,k,1)*U(i,j,k,1) + U(i,j,k,2)*U(i,j,k,2) + U(i,j,k,3)*U(i,j,k,3))/(c*c)); + const amrex::Real gamma = std::sqrt(1.0_rt + (U(i,j,k,1)*U(i,j,k,1) + U(i,j,k,2)*U(i,j,k,2) + U(i,j,k,3)*U(i,j,k,3))/(c*c)); return U(i,j,k,comp+1)/gamma; } // mindmod @@ -96,7 +96,7 @@ amrex::Real flux_N (const amrex::Array4& Um, const amrex::Array4& Um, const amrex::Array4 int i, int j, int k, amrex::Real Vm, amrex::Real Vp) { using namespace amrex::literals; - amrex::Real c = std::max( std::abs(Vm) , std::abs(Vp) ); + const amrex::Real c = std::max( std::abs(Vm) , std::abs(Vp) ); return 0.5_rt*(Vm*Um(i,j,k,0)*Um(i,j,k,1) + Vp*Up(i,j,k,0)*Up(i,j,k,1)) - (0.5_rt*c)*(Up(i,j,k,0)*Up(i,j,k,1) - Um(i,j,k,0)*Um(i,j,k,1)); } @@ -115,7 +115,7 @@ amrex::Real flux_NUy (const amrex::Array4& Um, const amrex::Array4 int i, int j, int k, amrex::Real Vm, amrex::Real Vp) { using namespace amrex::literals; - amrex::Real c = std::max( std::abs(Vm) , std::abs(Vp) ); + const amrex::Real c = std::max( std::abs(Vm) , std::abs(Vp) ); return 0.5_rt*(Vm*Um(i,j,k,0)*Um(i,j,k,2) + Vp*Up(i,j,k,0)*Up(i,j,k,2)) - (0.5_rt*c)*(Up(i,j,k,0)*Up(i,j,k,2) - Um(i,j,k,0)*Um(i,j,k,2)); } @@ -125,7 +125,7 @@ amrex::Real flux_NUz (const amrex::Array4& Um, const amrex::Array4 int i, int j, int k, amrex::Real Vm, amrex::Real Vp) { using namespace amrex::literals; - amrex::Real c = std::max( std::abs(Vm) , std::abs(Vp) ); + const amrex::Real c = std::max( std::abs(Vm) , std::abs(Vp) ); return 0.5_rt*(Vm*Um(i,j,k,0)*Um(i,j,k,3) + Vp*Up(i,j,k,0)*Up(i,j,k,3)) - (0.5_rt*c)*(Up(i,j,k,0)*Up(i,j,k,3) - Um(i,j,k,0)*Um(i,j,k,3)); } @@ -170,7 +170,7 @@ amrex::Real ave_stage2 (amrex::Real dQ, amrex::Real a, amrex::Real b, amrex::Rea using namespace amrex::literals; // sigma = sqrt(3) -1 constexpr auto sigma = 0.732050807568877_rt; - amrex::Real dq_min = 2.0_rt*std::min( b - min3(a,b,c), max3(a,b,c) - b); + const amrex::Real dq_min = 2.0_rt*std::min( b - min3(a,b,c), max3(a,b,c) - b); return ( std::abs(dQ)/dQ ) * std::min( std::abs(dQ) , sigma*std::abs(dq_min) ); } // Returns the offset indices for the "plus" grid @@ -499,10 +499,10 @@ const amrex::Array4& U_plus,int i,int j,int k,amrex::Real clight, i int ip, jp, kp; plus_index_offsets(i, j, k, ip, jp, kp, dir); - amrex::Real V_L_minus = V_calc(U_minus,ip,jp,kp,dir,clight); - amrex::Real V_I_minus = V_calc(U_minus,i,j,k,dir,clight); - amrex::Real V_L_plus = V_calc(U_plus,ip,jp,kp,dir,clight); - amrex::Real V_I_plus = V_calc(U_plus,i,j,k,dir,clight); + const amrex::Real V_L_minus = V_calc(U_minus,ip,jp,kp,dir,clight); + const amrex::Real V_I_minus = V_calc(U_minus,i,j,k,dir,clight); + const amrex::Real V_L_plus = V_calc(U_plus,ip,jp,kp,dir,clight); + const amrex::Real V_I_plus = V_calc(U_plus,i,j,k,dir,clight); // Flux differences depending on the component to compute if (comp == 0){ diff --git a/Source/Fluids/WarpXFluidContainer.cpp b/Source/Fluids/WarpXFluidContainer.cpp index 1915e4bf772..11f25678dc0 100644 --- a/Source/Fluids/WarpXFluidContainer.cpp +++ b/Source/Fluids/WarpXFluidContainer.cpp @@ -59,7 +59,7 @@ void WarpXFluidContainer::ReadParameters() { // Extract charge, mass, species type - std::string injection_style = "none"; + const std::string injection_style = "none"; SpeciesUtils::extractSpeciesProperties(species_name, injection_style, charge, mass, physical_species); const ParmParse pp_species_name(species_name); @@ -141,7 +141,7 @@ void WarpXFluidContainer::ReadParameters() void WarpXFluidContainer::AllocateLevelMFs(int lev, const BoxArray &ba, const DistributionMapping &dm) { - int ncomps = 1; + const int ncomps = 1; const amrex::IntVect nguards(AMREX_D_DECL(2, 2, 2)); // set human-readable tag for each MultiFab @@ -189,29 +189,29 @@ void WarpXFluidContainer::InitData(int lev, amrex::Box init_box, amrex::Real cur for (MFIter mfi(*N[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi) { - amrex::Box tile_box = mfi.tilebox(N[lev]->ixType().toIntVect()); + amrex::Box const tile_box = mfi.tilebox(N[lev]->ixType().toIntVect()); amrex::Array4 const &N_arr = N[lev]->array(mfi); amrex::Array4 const &NUx_arr = NU[lev][0]->array(mfi); amrex::Array4 const &NUy_arr = NU[lev][1]->array(mfi); amrex::Array4 const &NUz_arr = NU[lev][2]->array(mfi); // Return the intersection of all cells and the ones we wish to update - amrex::Box init_box_intersection = init_box & tile_box; + amrex::Box const init_box_intersection = init_box & tile_box; amrex::ParallelFor(init_box_intersection, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { #if defined(WARPX_DIM_3D) - amrex::Real x = problo[0] + i * dx[0]; - amrex::Real y = problo[1] + j * dx[1]; + const amrex::Real x = problo[0] + i * dx[0]; + const amrex::Real y = problo[1] + j * dx[1]; amrex::Real z = problo[2] + k * dx[2]; #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - amrex::Real x = problo[0] + i * dx[0]; - amrex::Real y = 0.0_rt; + const amrex::Real x = problo[0] + i * dx[0]; + const amrex::Real y = 0.0_rt; amrex::Real z = problo[1] + j * dx[1]; #else - amrex::Real x = 0.0_rt; - amrex::Real y = 0.0_rt; + const amrex::Real x = 0.0_rt; + const amrex::Real y = 0.0_rt; amrex::Real z = problo[0] + i * dx[0]; #endif @@ -232,9 +232,9 @@ void WarpXFluidContainer::InitData(int lev, amrex::Box init_box, amrex::Real cur // Lorentz transform n, u (from lab to boosted frame) if (n > 0.0){ if (gamma_boost > 1._rt){ - amrex::Real gamma = std::sqrt(1.0_rt + (u.x*u.x + u.y*u.y + u.z*u.z)/(clight*clight)); - amrex::Real n_boosted = gamma_boost*n*( 1.0_rt - beta_boost*u.z/(gamma*clight) ); - amrex::Real uz_boosted = gamma_boost*(u.z - beta_boost*clight*gamma); + const amrex::Real gamma = std::sqrt(1.0_rt + (u.x*u.x + u.y*u.y + u.z*u.z)/(clight*clight)); + const amrex::Real n_boosted = gamma_boost*n*( 1.0_rt - beta_boost*u.z/(gamma*clight) ); + const amrex::Real uz_boosted = gamma_boost*(u.z - beta_boost*clight*gamma); u.z = uz_boosted; n = n_boosted; } @@ -320,10 +320,10 @@ void WarpXFluidContainer::ApplyBcFluidsAndComms (int lev) amrex::Box tile_box = mfi.tilebox(N[lev]->ixType().toIntVect()); - amrex::Array4 N_arr = N[lev]->array(mfi); - amrex::Array4 NUx_arr = NU[lev][0]->array(mfi); - amrex::Array4 NUy_arr = NU[lev][1]->array(mfi); - amrex::Array4 NUz_arr = NU[lev][2]->array(mfi); + const amrex::Array4 N_arr = N[lev]->array(mfi); + const amrex::Array4 NUx_arr = NU[lev][0]->array(mfi); + const amrex::Array4 NUy_arr = NU[lev][1]->array(mfi); + const amrex::Array4 NUz_arr = NU[lev][2]->array(mfi); //Grow the tilebox tile_box.grow(1); @@ -413,28 +413,28 @@ void WarpXFluidContainer::AdvectivePush_Muscl (int lev) const auto dx = geom.CellSizeArray(); const amrex::Real clight = PhysConst::c; #if defined(WARPX_DIM_3D) - amrex::Real dt_over_dx = (dt/dx[0]); - amrex::Real dt_over_dy = (dt/dx[1]); - amrex::Real dt_over_dz = (dt/dx[2]); - amrex::Real dt_over_dx_half = 0.5_rt*(dt/dx[0]); - amrex::Real dt_over_dy_half = 0.5_rt*(dt/dx[1]); - amrex::Real dt_over_dz_half = 0.5_rt*(dt/dx[2]); + const amrex::Real dt_over_dx = (dt/dx[0]); + const amrex::Real dt_over_dy = (dt/dx[1]); + const amrex::Real dt_over_dz = (dt/dx[2]); + const amrex::Real dt_over_dx_half = 0.5_rt*(dt/dx[0]); + const amrex::Real dt_over_dy_half = 0.5_rt*(dt/dx[1]); + const amrex::Real dt_over_dz_half = 0.5_rt*(dt/dx[2]); #elif defined(WARPX_DIM_XZ) - amrex::Real dt_over_dx_half = 0.5_rt*(dt/dx[0]); - amrex::Real dt_over_dz_half = 0.5_rt*(dt/dx[1]); - amrex::Real dt_over_dx = (dt/dx[0]); - amrex::Real dt_over_dz = (dt/dx[1]); + const amrex::Real dt_over_dx_half = 0.5_rt*(dt/dx[0]); + const amrex::Real dt_over_dz_half = 0.5_rt*(dt/dx[1]); + const amrex::Real dt_over_dx = (dt/dx[0]); + const amrex::Real dt_over_dz = (dt/dx[1]); #elif defined(WARPX_DIM_RZ) const auto problo = geom.ProbLoArray(); - amrex::Real dt_over_dx_half = 0.5_rt*(dt/dx[0]); - amrex::Real dt_over_dz_half = 0.5_rt*(dt/dx[1]); - amrex::Box const& domain = geom.Domain(); + const amrex::Real dt_over_dx_half = 0.5_rt*(dt/dx[0]); + const amrex::Real dt_over_dz_half = 0.5_rt*(dt/dx[1]); + const amrex::Box& domain = geom.Domain(); #else - amrex::Real dt_over_dz = (dt/dx[0]); - amrex::Real dt_over_dz_half = 0.5_rt*(dt/dx[0]); + const amrex::Real dt_over_dz = (dt/dx[0]); + const amrex::Real dt_over_dz_half = 0.5_rt*(dt/dx[0]); #endif - amrex::BoxArray ba = N[lev]->boxArray(); + const amrex::BoxArray ba = N[lev]->boxArray(); // Temporary Half-step values #if defined(WARPX_DIM_3D) @@ -464,15 +464,17 @@ void WarpXFluidContainer::AdvectivePush_Muscl (int lev) // Loop over a box with one extra gridpoint in the ghost region to avoid // an extra MPI communication between the edge value computation loop and // the flux calculation loop - amrex::Box tile_box = mfi.growntilebox(1); - - // Limit the grown box for RZ at r = 0, r_max + const amrex::Box tile_box = [&](){ + auto tt = mfi.growntilebox(1); #if defined (WARPX_DIM_RZ) - const int idir = 0; - const int n_cell = -1; - tile_box.growLo(idir, n_cell); - tile_box.growHi(idir, n_cell); + // Limit the grown box for RZ at r = 0, r_max + const int idir = 0; + const int n_cell = -1; + tt.growLo(idir, n_cell); + tt.growHi(idir, n_cell); #endif + return tt; + }(); amrex::Array4 const &N_arr = N[lev]->array(mfi); amrex::Array4 const &NUx_arr = NU[lev][0]->array(mfi); @@ -500,20 +502,20 @@ void WarpXFluidContainer::AdvectivePush_Muscl (int lev) //(i.e. the 4 components correspond to N + the 3 components of U) // Extract the temporary arrays for edge values #if defined(WARPX_DIM_3D) - amrex::Array4 U_minus_x = tmp_U_minus_x.array(mfi); - amrex::Array4 U_plus_x = tmp_U_plus_x.array(mfi); - amrex::Array4 U_minus_y = tmp_U_minus_y.array(mfi); - amrex::Array4 U_plus_y = tmp_U_plus_y.array(mfi); - amrex::Array4 U_minus_z = tmp_U_minus_z.array(mfi); - amrex::Array4 U_plus_z = tmp_U_plus_z.array(mfi); + const amrex::Array4 U_minus_x = tmp_U_minus_x.array(mfi); + const amrex::Array4 U_plus_x = tmp_U_plus_x.array(mfi); + const amrex::Array4 U_minus_y = tmp_U_minus_y.array(mfi); + const amrex::Array4 U_plus_y = tmp_U_plus_y.array(mfi); + const amrex::Array4 U_minus_z = tmp_U_minus_z.array(mfi); + const amrex::Array4 U_plus_z = tmp_U_plus_z.array(mfi); #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - amrex::Array4 U_minus_x = tmp_U_minus_x.array(mfi); - amrex::Array4 U_plus_x = tmp_U_plus_x.array(mfi); - amrex::Array4 U_minus_z = tmp_U_minus_z.array(mfi); - amrex::Array4 U_plus_z = tmp_U_plus_z.array(mfi); + const amrex::Array4 U_minus_x = tmp_U_minus_x.array(mfi); + const amrex::Array4 U_plus_x = tmp_U_plus_x.array(mfi); + const amrex::Array4 U_minus_z = tmp_U_minus_z.array(mfi); + const amrex::Array4 U_plus_z = tmp_U_plus_z.array(mfi); #else - amrex::Array4 U_minus_z = tmp_U_minus_z.array(mfi); - amrex::Array4 U_plus_z = tmp_U_plus_z.array(mfi); + const amrex::Array4 U_minus_z = tmp_U_minus_z.array(mfi); + const amrex::Array4 U_plus_z = tmp_U_plus_z.array(mfi); #endif amrex::ParallelFor(tile_box, @@ -530,87 +532,89 @@ void WarpXFluidContainer::AdvectivePush_Muscl (int lev) amrex::Real Uz = (NUz_arr(i, j, k) / N_arr(i,j,k)); // Compute useful quantities for J - amrex::Real c_sq = clight*clight; - amrex::Real gamma = std::sqrt(1.0_rt + (Ux*Ux + Uy*Uy + Uz*Uz)/(c_sq) ); - amrex::Real inv_c2_gamma3 = 1._rt/(c_sq*gamma*gamma*gamma); + const amrex::Real c_sq = clight*clight; + const amrex::Real gamma = std::sqrt(1.0_rt + (Ux*Ux + Uy*Uy + Uz*Uz)/(c_sq) ); + const amrex::Real inv_c2_gamma3 = 1._rt/(c_sq*gamma*gamma*gamma); // J represents are 4x4 matrices that show up in the advection // equations written as a function of U = {N, Ux, Uy, Uz}: // \partial_t U + Jx \partial_x U + Jy \partial_y U + Jz \partial_z U = 0 #if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ) - amrex::Real Vx = Ux/gamma; + const amrex::Real Vx = Ux/gamma; // Compute the non-zero element of Jx - amrex::Real J00x = Vx; - amrex::Real J01x = N_arr(i,j,k)*(1/gamma)*(1-Vx*Vx/c_sq); - amrex::Real J02x = -N_arr(i,j,k)*Uy*Ux*inv_c2_gamma3; - amrex::Real J03x = -N_arr(i,j,k)*Uz*Ux*inv_c2_gamma3; - amrex::Real J11x = Vx; - amrex::Real J22x = Vx; - amrex::Real J33x = Vx; + const amrex::Real J00x = Vx; + const amrex::Real J01x = N_arr(i,j,k)*(1/gamma)*(1-Vx*Vx/c_sq); + const amrex::Real J02x = -N_arr(i,j,k)*Uy*Ux*inv_c2_gamma3; + const amrex::Real J03x = -N_arr(i,j,k)*Uz*Ux*inv_c2_gamma3; + const amrex::Real J11x = Vx; + const amrex::Real J22x = Vx; + const amrex::Real J33x = Vx; + + amrex::Real dU0x, dU1x, dU2x, dU3x; // Compute the cell slopes x - amrex::Real dU0x = ave( DownDx_N(N_arr,i,j,k), UpDx_N(N_arr,i,j,k) ); - amrex::Real dU1x = ave( DownDx_U(N_arr,NUx_arr,Ux,i,j,k), UpDx_U(N_arr,NUx_arr,Ux,i,j,k) ); - amrex::Real dU2x = ave( DownDx_U(N_arr,NUy_arr,Uy,i,j,k), UpDx_U(N_arr,NUy_arr,Uy,i,j,k) ); - amrex::Real dU3x = ave( DownDx_U(N_arr,NUz_arr,Uz,i,j,k), UpDx_U(N_arr,NUz_arr,Uz,i,j,k) ); + dU0x = ave( DownDx_N(N_arr,i,j,k), UpDx_N(N_arr,i,j,k) ); + dU1x = ave( DownDx_U(N_arr,NUx_arr,Ux,i,j,k), UpDx_U(N_arr,NUx_arr,Ux,i,j,k) ); + dU2x = ave( DownDx_U(N_arr,NUy_arr,Uy,i,j,k), UpDx_U(N_arr,NUy_arr,Uy,i,j,k) ); + dU3x = ave( DownDx_U(N_arr,NUz_arr,Uz,i,j,k), UpDx_U(N_arr,NUz_arr,Uz,i,j,k) ); #endif #if defined(WARPX_DIM_3D) - amrex::Real Vy = Uy/gamma; + const amrex::Real Vy = Uy/gamma; // Compute the non-zero element of Jy - amrex::Real J00y = Vy; - amrex::Real J01y = -N_arr(i,j,k)*Ux*Uy*inv_c2_gamma3; - amrex::Real J02y = N_arr(i,j,k)*(1/gamma)*(1-Vy*Vy/c_sq); - amrex::Real J03y = -N_arr(i,j,k)*Uz*Uy*inv_c2_gamma3; - amrex::Real J11y = Vy; - amrex::Real J22y = Vy; - amrex::Real J33y = Vy; + const amrex::Real J00y = Vy; + const amrex::Real J01y = -N_arr(i,j,k)*Ux*Uy*inv_c2_gamma3; + const amrex::Real J02y = N_arr(i,j,k)*(1/gamma)*(1-Vy*Vy/c_sq); + const amrex::Real J03y = -N_arr(i,j,k)*Uz*Uy*inv_c2_gamma3; + const amrex::Real J11y = Vy; + const amrex::Real J22y = Vy; + const amrex::Real J33y = Vy; // Compute the cell slopes y - amrex::Real dU0y = ave( DownDy_N(N_arr,i,j,k), UpDy_N(N_arr,i,j,k) ); - amrex::Real dU1y = ave( DownDy_U(N_arr,NUx_arr,Ux,i,j,k), UpDy_U(N_arr,NUx_arr,Ux,i,j,k) ); - amrex::Real dU2y = ave( DownDy_U(N_arr,NUy_arr,Uy,i,j,k), UpDy_U(N_arr,NUy_arr,Uy,i,j,k) ); - amrex::Real dU3y = ave( DownDy_U(N_arr,NUz_arr,Uz,i,j,k), UpDy_U(N_arr,NUz_arr,Uz,i,j,k) ); + const amrex::Real dU0y = ave( DownDy_N(N_arr,i,j,k), UpDy_N(N_arr,i,j,k) ); + const amrex::Real dU1y = ave( DownDy_U(N_arr,NUx_arr,Ux,i,j,k), UpDy_U(N_arr,NUx_arr,Ux,i,j,k) ); + const amrex::Real dU2y = ave( DownDy_U(N_arr,NUy_arr,Uy,i,j,k), UpDy_U(N_arr,NUy_arr,Uy,i,j,k) ); + const amrex::Real dU3y = ave( DownDy_U(N_arr,NUz_arr,Uz,i,j,k), UpDy_U(N_arr,NUz_arr,Uz,i,j,k) ); #endif - amrex::Real Vz = Uz/gamma; + const amrex::Real Vz = Uz/gamma; // Compute the non-zero element of Jz - amrex::Real J00z = Vz; - amrex::Real J01z = -N_arr(i,j,k)*Ux*Uz*inv_c2_gamma3; - amrex::Real J02z = -N_arr(i,j,k)*Uy*Uz*inv_c2_gamma3; - amrex::Real J03z = N_arr(i,j,k)*(1/gamma)*(1-Vz*Vz/c_sq); - amrex::Real J11z = Vz; - amrex::Real J22z = Vz; - amrex::Real J33z = Vz; + const amrex::Real J00z = Vz; + const amrex::Real J01z = -N_arr(i,j,k)*Ux*Uz*inv_c2_gamma3; + const amrex::Real J02z = -N_arr(i,j,k)*Uy*Uz*inv_c2_gamma3; + const amrex::Real J03z = N_arr(i,j,k)*(1/gamma)*(1-Vz*Vz/c_sq); + const amrex::Real J11z = Vz; + const amrex::Real J22z = Vz; + const amrex::Real J33z = Vz; // Compute the cell slopes z - amrex::Real dU0z = ave( DownDz_N(N_arr,i,j,k), UpDz_N(N_arr,i,j,k) ); - amrex::Real dU1z = ave( DownDz_U(N_arr,NUx_arr,Ux,i,j,k), UpDz_U(N_arr,NUx_arr,Ux,i,j,k) ); - amrex::Real dU2z = ave( DownDz_U(N_arr,NUy_arr,Uy,i,j,k), UpDz_U(N_arr,NUy_arr,Uy,i,j,k) ); - amrex::Real dU3z = ave( DownDz_U(N_arr,NUz_arr,Uz,i,j,k), UpDz_U(N_arr,NUz_arr,Uz,i,j,k) ); + const amrex::Real dU0z = ave( DownDz_N(N_arr,i,j,k), UpDz_N(N_arr,i,j,k) ); + const amrex::Real dU1z = ave( DownDz_U(N_arr,NUx_arr,Ux,i,j,k), UpDz_U(N_arr,NUx_arr,Ux,i,j,k) ); + const amrex::Real dU2z = ave( DownDz_U(N_arr,NUy_arr,Uy,i,j,k), UpDz_U(N_arr,NUy_arr,Uy,i,j,k) ); + const amrex::Real dU3z = ave( DownDz_U(N_arr,NUz_arr,Uz,i,j,k), UpDz_U(N_arr,NUz_arr,Uz,i,j,k) ); // Select the specific implementation depending on dimensionality #if defined(WARPX_DIM_3D) // Compute U ([ N, U]) at the halfsteps (U_tilde) using the slopes (dU) - amrex::Real JdU0x = J00x*dU0x + J01x*dU1x + J02x*dU2x + J03x*dU3x; - amrex::Real JdU1x = J11x*dU1x ; - amrex::Real JdU2x = J22x*dU2x ; - amrex::Real JdU3x = J33x*dU3x; - amrex::Real JdU0y = J00y*dU0y + J01y*dU1y + J02y*dU2y + J03y*dU3y; - amrex::Real JdU1y = J11y*dU1y; - amrex::Real JdU2y = J22y*dU2y; - amrex::Real JdU3y = J33y*dU3y; - amrex::Real JdU0z = J00z*dU0z + J01z*dU1z + J02z*dU2z + J03z*dU3z; - amrex::Real JdU1z = J11z*dU1z; - amrex::Real JdU2z = J22z*dU2z; - amrex::Real JdU3z = J33z*dU3z; - amrex::Real U_tilde0 = N_arr(i,j,k) - dt_over_dx_half*JdU0x - dt_over_dy_half*JdU0y - dt_over_dz_half*JdU0z; - amrex::Real U_tilde1 = Ux - dt_over_dx_half*JdU1x - dt_over_dy_half*JdU1y - dt_over_dz_half*JdU1z; - amrex::Real U_tilde2 = Uy - dt_over_dx_half*JdU2x - dt_over_dy_half*JdU2y - dt_over_dz_half*JdU2z; - amrex::Real U_tilde3 = Uz - dt_over_dx_half*JdU3x - dt_over_dy_half*JdU3y - dt_over_dz_half*JdU3z; + const amrex::Real JdU0x = J00x*dU0x + J01x*dU1x + J02x*dU2x + J03x*dU3x; + const amrex::Real JdU1x = J11x*dU1x ; + const amrex::Real JdU2x = J22x*dU2x ; + const amrex::Real JdU3x = J33x*dU3x; + const amrex::Real JdU0y = J00y*dU0y + J01y*dU1y + J02y*dU2y + J03y*dU3y; + const amrex::Real JdU1y = J11y*dU1y; + const amrex::Real JdU2y = J22y*dU2y; + const amrex::Real JdU3y = J33y*dU3y; + const amrex::Real JdU0z = J00z*dU0z + J01z*dU1z + J02z*dU2z + J03z*dU3z; + const amrex::Real JdU1z = J11z*dU1z; + const amrex::Real JdU2z = J22z*dU2z; + const amrex::Real JdU3z = J33z*dU3z; + const amrex::Real U_tilde0 = N_arr(i,j,k) - dt_over_dx_half*JdU0x - dt_over_dy_half*JdU0y - dt_over_dz_half*JdU0z; + const amrex::Real U_tilde1 = Ux - dt_over_dx_half*JdU1x - dt_over_dy_half*JdU1y - dt_over_dz_half*JdU1z; + const amrex::Real U_tilde2 = Uy - dt_over_dx_half*JdU2x - dt_over_dy_half*JdU2y - dt_over_dz_half*JdU2z; + const amrex::Real U_tilde3 = Uz - dt_over_dx_half*JdU3x - dt_over_dy_half*JdU3y - dt_over_dz_half*JdU3z; // Predict U at the cell edges (x) @@ -636,12 +640,9 @@ void WarpXFluidContainer::AdvectivePush_Muscl (int lev) #elif defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ) - // Have no RZ-inertial source for primitive vars if in XZ - amrex::Real N_source = 0.0; - #if defined(WARPX_DIM_RZ) - amrex::Real dr = dx[0]; - amrex::Real r = problo[0] + i * dr; + const amrex::Real dr = dx[0]; + const amrex::Real r = problo[0] + i * dr; // Impose "none" boundaries // Condition: dUx = 0 at r = 0 if (i == domain.smallEnd(0)) { @@ -664,24 +665,26 @@ void WarpXFluidContainer::AdvectivePush_Muscl (int lev) } // RZ sources: - if (i != domain.smallEnd(0)) { - N_source = N_arr(i,j,k)*Vx/r; - } + const amrex::Real N_source = + (i != domain.smallEnd(0)) ? N_arr(i,j,k)*Vx/r : 0.0_rt; +#else + // Have no RZ-inertial source for primitive vars if in XZ + const amrex::Real N_source = 0.0; #endif // Compute U ([ N, U]) at the halfsteps (U_tilde) using the slopes (dU) - amrex::Real JdU0x = J00x*dU0x + J01x*dU1x + J02x*dU2x + J03x*dU3x; - amrex::Real JdU1x = J11x*dU1x; - amrex::Real JdU2x = J22x*dU2x; - amrex::Real JdU3x = J33x*dU3x; - amrex::Real JdU0z = J00z*dU0z + J01z*dU1z + J02z*dU2z + J03z*dU3z; - amrex::Real JdU1z = J11z*dU1z; - amrex::Real JdU2z = J22z*dU2z; - amrex::Real JdU3z = J33z*dU3z; - amrex::Real U_tilde0 = N_arr(i,j,k) - dt_over_dx_half*JdU0x - dt_over_dz_half*JdU0z - (dt/2.0_rt)*N_source; - amrex::Real U_tilde1 = Ux - dt_over_dx_half*JdU1x - dt_over_dz_half*JdU1z; - amrex::Real U_tilde2 = Uy - dt_over_dx_half*JdU2x - dt_over_dz_half*JdU2z; - amrex::Real U_tilde3 = Uz - dt_over_dx_half*JdU3x - dt_over_dz_half*JdU3z; + const amrex::Real JdU0x = J00x*dU0x + J01x*dU1x + J02x*dU2x + J03x*dU3x; + const amrex::Real JdU1x = J11x*dU1x; + const amrex::Real JdU2x = J22x*dU2x; + const amrex::Real JdU3x = J33x*dU3x; + const amrex::Real JdU0z = J00z*dU0z + J01z*dU1z + J02z*dU2z + J03z*dU3z; + const amrex::Real JdU1z = J11z*dU1z; + const amrex::Real JdU2z = J22z*dU2z; + const amrex::Real JdU3z = J33z*dU3z; + const amrex::Real U_tilde0 = N_arr(i,j,k) - dt_over_dx_half*JdU0x - dt_over_dz_half*JdU0z - (dt/2.0_rt)*N_source; + const amrex::Real U_tilde1 = Ux - dt_over_dx_half*JdU1x - dt_over_dz_half*JdU1z; + const amrex::Real U_tilde2 = Uy - dt_over_dx_half*JdU2x - dt_over_dz_half*JdU2z; + const amrex::Real U_tilde3 = Uz - dt_over_dx_half*JdU3x - dt_over_dz_half*JdU3z; // Predict U at the cell edges (x) compute_U_edges(U_minus_x, U_plus_x, i, j, k, box_x, U_tilde0, U_tilde1, U_tilde2, U_tilde3, dU0x, dU1x, dU2x, dU3x,0); @@ -700,14 +703,14 @@ void WarpXFluidContainer::AdvectivePush_Muscl (int lev) #else // Compute U ([ N, U]) at the halfsteps (U_tilde) using the slopes (dU) - amrex::Real JdU0z = J00z*dU0z + J01z*dU1z + J02z*dU2z + J03z*dU3z; - amrex::Real JdU1z = J11z*dU1z; - amrex::Real JdU2z = J22z*dU2z; - amrex::Real JdU3z = J33z*dU3z; - amrex::Real U_tilde0 = N_arr(i,j,k) - dt_over_dz_half*JdU0z; - amrex::Real U_tilde1 = Ux - dt_over_dz_half*JdU1z; - amrex::Real U_tilde2 = Uy - dt_over_dz_half*JdU2z; - amrex::Real U_tilde3 = Uz - dt_over_dz_half*JdU3z; + const amrex::Real JdU0z = J00z*dU0z + J01z*dU1z + J02z*dU2z + J03z*dU3z; + const amrex::Real JdU1z = J11z*dU1z; + const amrex::Real JdU2z = J22z*dU2z; + const amrex::Real JdU3z = J33z*dU3z; + const amrex::Real U_tilde0 = N_arr(i,j,k) - dt_over_dz_half*JdU0z; + const amrex::Real U_tilde1 = Ux - dt_over_dz_half*JdU1z; + const amrex::Real U_tilde2 = Uy - dt_over_dz_half*JdU2z; + const amrex::Real U_tilde3 = Uz - dt_over_dz_half*JdU3z; // Predict U at the cell edges (z) compute_U_edges(U_minus_z, U_plus_z, i, j, k, box_z, U_tilde0, U_tilde1, U_tilde2, U_tilde3, dU0z, dU1z, dU2z, dU3z,2); @@ -740,11 +743,11 @@ void WarpXFluidContainer::AdvectivePush_Muscl (int lev) #endif for (MFIter mfi(*N[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi) { - amrex::Box tile_box = mfi.tilebox(N[lev]->ixType().toIntVect()); - amrex::Array4 N_arr = N[lev]->array(mfi); - amrex::Array4 NUx_arr = NU[lev][0]->array(mfi); - amrex::Array4 NUy_arr = NU[lev][1]->array(mfi); - amrex::Array4 NUz_arr = NU[lev][2]->array(mfi); + const amrex::Box tile_box = mfi.tilebox(N[lev]->ixType().toIntVect()); + const amrex::Array4 N_arr = N[lev]->array(mfi); + const amrex::Array4 NUx_arr = NU[lev][0]->array(mfi); + const amrex::Array4 NUy_arr = NU[lev][1]->array(mfi); + const amrex::Array4 NUz_arr = NU[lev][2]->array(mfi); #if defined(WARPX_DIM_3D) amrex::Array4 const &U_minus_x = tmp_U_minus_x.array(mfi); @@ -800,9 +803,9 @@ void WarpXFluidContainer::AdvectivePush_Muscl (int lev) // Compute the flux areas for RZ // Cell-centered radius - amrex::Real dr = dx[0]; - amrex::Real dz = dx[1]; - amrex::Real r = problo[0] + i * dr; + const amrex::Real dr = dx[0]; + const amrex::Real dz = dx[1]; + const amrex::Real r = problo[0] + i * dr; amrex::Real Vij = 0.0_rt; amrex::Real S_Az = 0.0_rt; @@ -830,8 +833,8 @@ void WarpXFluidContainer::AdvectivePush_Muscl (int lev) // Impose "none" boundaries // Condition: Vx(r) = 0 at boundaries - amrex::Real Vx_I_minus = V_calc(U_minus_x,i,j,k,0,clight); - amrex::Real Vx_L_plus = V_calc(U_plus_x,i-1,j,k,0,clight); + const amrex::Real Vx_I_minus = V_calc(U_minus_x,i,j,k,0,clight); + const amrex::Real Vx_L_plus = V_calc(U_plus_x,i-1,j,k,0,clight); // compute the fluxes: // (note that _plus is shifted due to grid location) @@ -897,8 +900,8 @@ void WarpXFluidContainer::centrifugal_source_rz (int lev) amrex::Box const &tile_box = mfi.tilebox(N[lev]->ixType().toIntVect()); amrex::Array4 const &N_arr = N[lev]->array(mfi); - amrex::Array4 NUx_arr = NU[lev][0]->array(mfi); - amrex::Array4 NUy_arr = NU[lev][1]->array(mfi); + const amrex::Array4 NUx_arr = NU[lev][0]->array(mfi); + const amrex::Array4 NUy_arr = NU[lev][1]->array(mfi); amrex::Array4 const &NUz_arr = NU[lev][2]->array(mfi); amrex::ParallelFor(tile_box, @@ -909,20 +912,20 @@ void WarpXFluidContainer::centrifugal_source_rz (int lev) if (N_arr(i,j,k)>0.0_rt) { // Compute r - amrex::Real r = problo[0] + i * dx[0]; + const amrex::Real r = problo[0] + i * dx[0]; // Isolate U from NU amrex::Real u_r = (NUx_arr(i, j, k) / (N_arr(i,j,k) * clight )); amrex::Real u_theta = (NUy_arr(i, j, k) / (N_arr(i,j,k) * clight )); - amrex::Real u_z = (NUz_arr(i, j, k) / (N_arr(i,j,k) * clight )); + const amrex::Real u_z = (NUz_arr(i, j, k) / (N_arr(i,j,k) * clight )); // (SSP-RK3) Push the fluid momentum (R and Theta) // F_r, F_theta are first order euler pushes of our rhs operator if (i != domain.smallEnd(0)) { - amrex::Real u_r_1 = F_r(r,u_r,u_theta,u_z,dt); - amrex::Real u_theta_1 = F_theta(r,u_r,u_theta,u_z,dt); - amrex::Real u_r_2 = (0.75_rt)*(u_r) + (0.25_rt)*F_r(r,u_r_1,u_theta_1,u_z,dt); - amrex::Real u_theta_2 = (0.75_rt)*(u_theta) + (0.25_rt)*F_theta(r,u_r_1,u_theta_1,u_z,dt); + const amrex::Real u_r_1 = F_r(r,u_r,u_theta,u_z,dt); + const amrex::Real u_theta_1 = F_theta(r,u_r,u_theta,u_z,dt); + const amrex::Real u_r_2 = (0.75_rt)*(u_r) + (0.25_rt)*F_r(r,u_r_1,u_theta_1,u_z,dt); + const amrex::Real u_theta_2 = (0.75_rt)*(u_theta) + (0.25_rt)*F_theta(r,u_r_1,u_theta_1,u_z,dt); u_r = (1.0_rt/3.0_rt)*(u_r) + (2.0_rt/3.0_rt)*F_r(r,u_r_2,u_theta_2,u_z,dt); u_theta = (1.0_rt/3.0_rt)*(u_theta) + (2.0_rt/3.0_rt)*F_theta(r,u_r_2,u_theta_2,u_z,dt); @@ -1018,9 +1021,9 @@ void WarpXFluidContainer::GatherAndPush ( amrex::Box const &tile_box = mfi.tilebox(N[lev]->ixType().toIntVect()); amrex::Array4 const &N_arr = N[lev]->array(mfi); - amrex::Array4 NUx_arr = NU[lev][0]->array(mfi); - amrex::Array4 NUy_arr = NU[lev][1]->array(mfi); - amrex::Array4 NUz_arr = NU[lev][2]->array(mfi); + const amrex::Array4 NUx_arr = NU[lev][0]->array(mfi); + const amrex::Array4 NUy_arr = NU[lev][1]->array(mfi); + const amrex::Array4 NUz_arr = NU[lev][2]->array(mfi); amrex::Array4 const& Ex_arr = Ex.array(mfi); amrex::Array4 const& Ey_arr = Ey.array(mfi); @@ -1030,7 +1033,7 @@ void WarpXFluidContainer::GatherAndPush ( amrex::Array4 const& Bz_arr = Bz.array(mfi); // Here, we do not perform any coarsening. - amrex::GpuArray coarsening_ratio = {1, 1, 1}; + const amrex::GpuArray coarsening_ratio = {1, 1, 1}; amrex::ParallelFor(tile_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept @@ -1064,23 +1067,23 @@ void WarpXFluidContainer::GatherAndPush ( // Grab the location #if defined(WARPX_DIM_3D) - amrex::Real x = problo[0] + i * dx[0]; - amrex::Real y = problo[1] + j * dx[1]; - amrex::Real z = problo[2] + k * dx[2]; + const amrex::Real x = problo[0] + i * dx[0]; + const amrex::Real y = problo[1] + j * dx[1]; + const amrex::Real z = problo[2] + k * dx[2]; #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - amrex::Real x = problo[0] + i * dx[0]; - amrex::Real y = 0.0_rt; - amrex::Real z = problo[1] + j * dx[1]; + const amrex::Real x = problo[0] + i * dx[0]; + const amrex::Real y = 0.0_rt; + const amrex::Real z = problo[1] + j * dx[1]; #else - amrex::Real x = 0.0_rt; - amrex::Real y = 0.0_rt; - amrex::Real z = problo[0] + i * dx[0]; + const amrex::Real x = 0.0_rt; + const amrex::Real y = 0.0_rt; + const amrex::Real z = problo[0] + i * dx[0]; #endif // Get the lab frame E and B // Transform (boosted to lab) - amrex::Real t_lab = gamma_boost*(t + beta_boost*z/PhysConst::c); - amrex::Real z_lab = gamma_boost*(z + beta_boost*PhysConst::c*t); + const amrex::Real t_lab = gamma_boost*(t + beta_boost*z/PhysConst::c); + const amrex::Real z_lab = gamma_boost*(z + beta_boost*PhysConst::c*t); // Grab the external fields in the lab frame: if ( external_e_fields ) { @@ -1125,17 +1128,17 @@ void WarpXFluidContainer::GatherAndPush ( // Added external e fields: if ( external_e_fields ){ #if defined(WARPX_DIM_3D) - amrex::Real x = problo[0] + i * dx[0]; - amrex::Real y = problo[1] + j * dx[1]; - amrex::Real z = problo[2] + k * dx[2]; + const amrex::Real x = problo[0] + i * dx[0]; + const amrex::Real y = problo[1] + j * dx[1]; + const amrex::Real z = problo[2] + k * dx[2]; #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - amrex::Real x = problo[0] + i * dx[0]; - amrex::Real y = 0.0_rt; - amrex::Real z = problo[1] + j * dx[1]; + const amrex::Real x = problo[0] + i * dx[0]; + const amrex::Real y = 0.0_rt; + const amrex::Real z = problo[1] + j * dx[1]; #else - amrex::Real x = 0.0_rt; - amrex::Real y = 0.0_rt; - amrex::Real z = problo[0] + i * dx[0]; + const amrex::Real x = 0.0_rt; + const amrex::Real y = 0.0_rt; + const amrex::Real z = problo[0] + i * dx[0]; #endif Ex_Nodal += Exfield_parser(x, y, z, t); @@ -1146,17 +1149,17 @@ void WarpXFluidContainer::GatherAndPush ( // Added external b fields: if ( external_b_fields ){ #if defined(WARPX_DIM_3D) - amrex::Real x = problo[0] + i * dx[0]; - amrex::Real y = problo[1] + j * dx[1]; - amrex::Real z = problo[2] + k * dx[2]; + const amrex::Real x = problo[0] + i * dx[0]; + const amrex::Real y = problo[1] + j * dx[1]; + const amrex::Real z = problo[2] + k * dx[2]; #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - amrex::Real x = problo[0] + i * dx[0]; - amrex::Real y = 0.0_rt; - amrex::Real z = problo[1] + j * dx[1]; + const amrex::Real x = problo[0] + i * dx[0]; + const amrex::Real y = 0.0_rt; + const amrex::Real z = problo[1] + j * dx[1]; #else - amrex::Real x = 0.0_rt; - amrex::Real y = 0.0_rt; - amrex::Real z = problo[0] + i * dx[0]; + const amrex::Real x = 0.0_rt; + const amrex::Real y = 0.0_rt; + const amrex::Real z = problo[0] + i * dx[0]; #endif Bx_Nodal += Bxfield_parser(x, y, z, t); @@ -1217,8 +1220,8 @@ void WarpXFluidContainer::DepositCharge (int lev, amrex::MultiFab &rho, int icom amrex::Box const &tile_box = mfi.tilebox(N[lev]->ixType().toIntVect()); amrex::Array4 const &N_arr = N[lev]->array(mfi); - amrex::Array4 rho_arr = rho.array(mfi); - amrex::Array4 owner_mask_rho_arr = owner_mask_rho->array(mfi); + const amrex::Array4 rho_arr = rho.array(mfi); + const amrex::Array4 owner_mask_rho_arr = owner_mask_rho->array(mfi); // Deposit Rho amrex::ParallelFor(tile_box, @@ -1279,9 +1282,9 @@ void WarpXFluidContainer::DepositCurrent( amrex::Array4 const &NUy_arr = NU[lev][1]->array(mfi); amrex::Array4 const &NUz_arr = NU[lev][2]->array(mfi); - amrex::Array4 tmp_jx_fluid_arr = tmp_jx_fluid.array(mfi); - amrex::Array4 tmp_jy_fluid_arr = tmp_jy_fluid.array(mfi); - amrex::Array4 tmp_jz_fluid_arr = tmp_jz_fluid.array(mfi); + const amrex::Array4 tmp_jx_fluid_arr = tmp_jx_fluid.array(mfi); + const amrex::Array4 tmp_jy_fluid_arr = tmp_jy_fluid.array(mfi); + const amrex::Array4 tmp_jz_fluid_arr = tmp_jz_fluid.array(mfi); amrex::ParallelFor(tile_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept @@ -1311,21 +1314,21 @@ void WarpXFluidContainer::DepositCurrent( amrex::Box const &tile_box_y = mfi.tilebox(jy.ixType().toIntVect()); amrex::Box const &tile_box_z = mfi.tilebox(jz.ixType().toIntVect()); - amrex::Array4 jx_arr = jx.array(mfi); - amrex::Array4 jy_arr = jy.array(mfi); - amrex::Array4 jz_arr = jz.array(mfi); + const amrex::Array4 jx_arr = jx.array(mfi); + const amrex::Array4 jy_arr = jy.array(mfi); + const amrex::Array4 jz_arr = jz.array(mfi); - amrex::Array4 tmp_jx_fluid_arr = tmp_jx_fluid.array(mfi); - amrex::Array4 tmp_jy_fluid_arr = tmp_jy_fluid.array(mfi); - amrex::Array4 tmp_jz_fluid_arr = tmp_jz_fluid.array(mfi); + const amrex::Array4 tmp_jx_fluid_arr = tmp_jx_fluid.array(mfi); + const amrex::Array4 tmp_jy_fluid_arr = tmp_jy_fluid.array(mfi); + const amrex::Array4 tmp_jz_fluid_arr = tmp_jz_fluid.array(mfi); - amrex::Array4 owner_mask_x_arr = owner_mask_x->array(mfi); - amrex::Array4 owner_mask_y_arr = owner_mask_y->array(mfi); - amrex::Array4 owner_mask_z_arr = owner_mask_z->array(mfi); + const amrex::Array4 owner_mask_x_arr = owner_mask_x->array(mfi); + const amrex::Array4 owner_mask_y_arr = owner_mask_y->array(mfi); + const amrex::Array4 owner_mask_z_arr = owner_mask_z->array(mfi); // When using the `Interp` function, one needs to specify whether coarsening is desired. // Here, we do not perform any coarsening. - amrex::GpuArray coarsening_ratio = {1, 1, 1}; + const amrex::GpuArray coarsening_ratio = {1, 1, 1}; // Interpolate fluid current and deposit it @@ -1333,19 +1336,19 @@ void WarpXFluidContainer::DepositCurrent( amrex::ParallelFor( tile_box_x, tile_box_y, tile_box_z, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - amrex::Real jx_tmp = ablastr::coarsen::sample::Interp(tmp_jx_fluid_arr, + const amrex::Real jx_tmp = ablastr::coarsen::sample::Interp(tmp_jx_fluid_arr, j_nodal_type, jx_type, coarsening_ratio, i, j, k, 0); if ( owner_mask_x_arr(i,j,k) ) { jx_arr(i, j, k) += jx_tmp; } }, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - amrex::Real jy_tmp = ablastr::coarsen::sample::Interp(tmp_jy_fluid_arr, + const amrex::Real jy_tmp = ablastr::coarsen::sample::Interp(tmp_jy_fluid_arr, j_nodal_type, jy_type, coarsening_ratio, i, j, k, 0); if ( owner_mask_y_arr(i,j,k) ) { jy_arr(i, j, k) += jy_tmp; } }, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - amrex::Real jz_tmp = ablastr::coarsen::sample::Interp(tmp_jz_fluid_arr, + const amrex::Real jz_tmp = ablastr::coarsen::sample::Interp(tmp_jz_fluid_arr, j_nodal_type, jz_type, coarsening_ratio, i, j, k, 0); if ( owner_mask_z_arr(i,j,k) ) { jz_arr(i, j, k) += jz_tmp; } } diff --git a/Source/Initialization/ExternalField.cpp b/Source/Initialization/ExternalField.cpp index 909e5c01977..e7a272e2e95 100644 --- a/Source/Initialization/ExternalField.cpp +++ b/Source/Initialization/ExternalField.cpp @@ -175,7 +175,7 @@ ExternalFieldParams::ExternalFieldParams(const amrex::ParmParse& pp_warpx) if (E_ext_grid_type == ExternalFieldType::read_from_file || B_ext_grid_type == ExternalFieldType::read_from_file){ - std::string read_fields_from_path="./"; + const std::string read_fields_from_path="./"; pp_warpx.query("read_fields_from_path", external_fields_path); } //___________________________________________________________________________ diff --git a/Source/Initialization/InjectorPosition.H b/Source/Initialization/InjectorPosition.H index c86ed432bd1..7168ac93ebb 100644 --- a/Source/Initialization/InjectorPosition.H +++ b/Source/Initialization/InjectorPosition.H @@ -227,9 +227,9 @@ struct InjectorPosition bool overlapsWith (const amrex::XDim3& lo, const amrex::XDim3& hi) const noexcept { - return ! ( (xmin > hi.x) || (xmax < lo.x) - || (ymin > hi.y) || (ymax < lo.y) - || (zmin > hi.z) || (zmax < lo.z) ); + return ( (xmin <= hi.x) && (xmax >= lo.x) + && (ymin <= hi.y) && (ymax >= lo.y) + && (zmin <= hi.z) && (zmax >= lo.z) ); } private: diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp index 7ae5a9db3d3..2e338b20c09 100644 --- a/Source/Initialization/PlasmaInjector.cpp +++ b/Source/Initialization/PlasmaInjector.cpp @@ -577,9 +577,9 @@ bool PlasmaInjector::insideBounds (amrex::Real x, amrex::Real y, amrex::Real z) bool PlasmaInjector::overlapsWith (const amrex::XDim3& lo, const amrex::XDim3& hi) const noexcept { - return ! ( (xmin > hi.x) || (xmax < lo.x) - || (ymin > hi.y) || (ymax < lo.y) - || (zmin > hi.z) || (zmax < lo.z) ); + return ( (xmin <= hi.x) && (xmax >= lo.x) + && (ymin <= hi.y) && (ymax >= lo.y) + && (zmin <= hi.z) && (zmax >= lo.z) ); } bool diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index eee3012ab4d..e3ec624be8e 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -1478,8 +1478,8 @@ WarpX::ReadExternalFieldFromFile ( #endif #if defined(WARPX_DIM_RZ) - amrex::Array4 fc_array(FC_data, {0,0,0}, {extent0, extent2, extent1}, 1); - double + const amrex::Array4 fc_array(FC_data, {0,0,0}, {extent0, extent2, extent1}, 1); + const double f00 = fc_array(0, iz , ir ), f01 = fc_array(0, iz , ir+1), f10 = fc_array(0, iz+1, ir ), diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H index 78381f22258..a9e081d44b4 100644 --- a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H +++ b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H @@ -420,7 +420,7 @@ public: amrex::Box const& cbx = mfi.tilebox(amrex::IntVect::TheZeroVector()); //Cell-centered box const auto lo = lbound(cbx); const auto hi = ubound(cbx); - int nz = hi.y-lo.y+1; + const int nz = hi.y-lo.y+1; auto dr = geom.CellSize(0); auto dz = geom.CellSize(1); #elif defined(WARPX_DIM_3D) diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/CollisionFilterFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/CollisionFilterFunc.H index c714cfdb133..64f05650d1f 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/CollisionFilterFunc.H +++ b/Source/Particles/Collision/BinaryCollision/DSMC/CollisionFilterFunc.H @@ -68,7 +68,7 @@ void CollisionPairFilter (const amrex::ParticleReal& u1x, const amrex::ParticleR } // calculate total collision probability - amrex::ParticleReal exponent = ( + const amrex::ParticleReal exponent = ( lab_to_COM_factor * multiplier_ratio * w_max * sigma_tot * v_coll * dt / dV ); diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMC.cpp b/Source/Particles/Collision/BinaryCollision/DSMC/DSMC.cpp index 70a14712a0c..92ae3ddf257 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMC.cpp +++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMC.cpp @@ -15,7 +15,7 @@ DSMC::DSMC (const std::string collision_name) : CollisionBase(collision_name) { using namespace amrex::literals; - amrex::ParmParse pp_collision_name(collision_name); + const amrex::ParmParse pp_collision_name(collision_name); #if defined WARPX_DIM_RZ amrex::Abort("DSMC collisions are only implemented for Cartesian coordinates."); @@ -34,7 +34,7 @@ DSMC::DSMC (const std::string collision_name) // create a vector of ScatteringProcess objects from each scattering // process name for (const auto& scattering_process : scattering_process_names) { - std::string kw_cross_section = scattering_process + "_cross_section"; + const std::string kw_cross_section = scattering_process + "_cross_section"; std::string cross_section_file; pp_collision_name.query(kw_cross_section.c_str(), cross_section_file); @@ -43,7 +43,7 @@ DSMC::DSMC (const std::string collision_name) amrex::ParticleReal energy = 0._prt; if (scattering_process.find("excitation") != std::string::npos || scattering_process.find("ionization") != std::string::npos) { - std::string kw_energy = scattering_process + "_energy"; + const std::string kw_energy = scattering_process + "_energy"; utils::parser::getWithParser( pp_collision_name, kw_energy.c_str(), energy); } @@ -83,8 +83,8 @@ DSMC::doCollisions (amrex::Real /*cur_time*/, amrex::Real dt, MultiParticleConta // SmartCopy objects are created that will facilitate the particle splitting // operation involved in DSMC collisions between particles with arbitrary // weights. - SmartCopyFactory copy_factory_species1(species1, species1); - SmartCopyFactory copy_factory_species2(species2, species2); + const auto copy_factory_species1 = SmartCopyFactory{species1, species1}; + const auto copy_factory_species2 = SmartCopyFactory{species2, species2}; auto copy_species1 = copy_factory_species1.getSmartCopy(); auto copy_species2 = copy_factory_species2.getSmartCopy(); @@ -158,13 +158,13 @@ DSMC::doCollisionsWithinTile( // - Species 1 index_type* AMREX_RESTRICT indices_1 = bins_1.permutationPtr(); index_type const* AMREX_RESTRICT cell_offsets_1 = bins_1.offsetsPtr(); - amrex::ParticleReal m1 = species_1.getMass(); + const amrex::ParticleReal m1 = species_1.getMass(); const auto ptd_1 = ptile_1.getParticleTileData(); // - Species 2 index_type* AMREX_RESTRICT indices_2 = bins_2.permutationPtr(); index_type const* AMREX_RESTRICT cell_offsets_2 = bins_2.offsetsPtr(); - amrex::ParticleReal m2 = species_2.getMass(); + const amrex::ParticleReal m2 = species_2.getMass(); const auto ptd_2 = ptile_2.getParticleTileData(); amrex::Geometry const& geom = WarpX::GetInstance().Geom(lev); @@ -176,7 +176,7 @@ DSMC::doCollisionsWithinTile( amrex::Box const& cbx = mfi.tilebox(amrex::IntVect::TheZeroVector()); //Cell-centered box const auto lo = lbound(cbx); const auto hi = ubound(cbx); - int nz = hi.y-lo.y+1; + const int nz = hi.y-lo.y+1; auto dr = geom.CellSize(0); auto dz = geom.CellSize(1); #elif defined(WARPX_DIM_3D) @@ -262,7 +262,7 @@ DSMC::doCollisionsWithinTile( ShuffleFisherYates(indices_1, cell_start_1, cell_stop_1, engine); ShuffleFisherYates(indices_2, cell_start_2, cell_stop_2, engine); #if defined WARPX_DIM_RZ - int ri = (i_cell - i_cell%nz) / nz; + const int ri = (i_cell - i_cell%nz) / nz; auto dV = MathConst::pi*(2.0_prt*ri+1.0_prt)*dr*dr*dz; #endif // Call the function in order to perform collisions diff --git a/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H b/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H index 7a2853e3db5..d53ce4348fb 100644 --- a/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H +++ b/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H @@ -267,8 +267,8 @@ public: // Initialize the user runtime components for (int i = 0; i < m_num_product_species; i++) { - int start_index = int(products_np[i]); - int stop_index = int(products_np[i] + num_added_vec[i]); + const auto start_index = int(products_np[i]); + const auto stop_index = int(products_np[i] + num_added_vec[i]); ParticleCreation::DefaultInitializeRuntimeAttributes(*tile_products[i], 0, 0, pc_products[i]->getUserRealAttribs(), pc_products[i]->getUserIntAttribs(), diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index 2252a63fd07..8d4c5b6a5f4 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -1381,8 +1381,8 @@ void doChargeConservingDepositionShapeNImplicit (const amrex::ParticleReal * con // --- Compute shape factors // Compute shape factors for position as they are now and at old positions // [ijk]_new: leftmost grid point that the particle touches - Compute_shape_factor< depos_order > compute_shape_factor; - Compute_shifted_shape_factor< depos_order > compute_shifted_shape_factor; + const Compute_shape_factor< depos_order > compute_shape_factor; + const Compute_shifted_shape_factor< depos_order > compute_shifted_shape_factor; #if !defined(WARPX_DIM_1D_Z) const int i_new = compute_shape_factor(sx_new+1, x_new); @@ -1748,19 +1748,19 @@ void doVillasenorDepositionShapeNImplicit (const amrex::ParticleReal * const xp_ // compute cell crossings in X-direction const auto i_old = static_cast(x_old-shift); const auto i_new = static_cast(x_new-shift); - int cell_crossings_x = std::abs(i_new-i_old); + const int cell_crossings_x = std::abs(i_new-i_old); num_segments += cell_crossings_x; // compute cell crossings in Y-direction const auto j_old = static_cast(y_old-shift); const auto j_new = static_cast(y_new-shift); - int cell_crossings_y = std::abs(j_new-j_old); + const int cell_crossings_y = std::abs(j_new-j_old); num_segments += cell_crossings_y; // compute cell crossings in Z-direction const auto k_old = static_cast(z_old-shift); const auto k_new = static_cast(z_new-shift); - int cell_crossings_z = std::abs(k_new-k_old); + const int cell_crossings_z = std::abs(k_new-k_old); num_segments += cell_crossings_z; // need to assert that the number of cell crossings in each direction @@ -1783,8 +1783,8 @@ void doVillasenorDepositionShapeNImplicit (const amrex::ParticleReal * const xp_ } // loop over the number of segments and deposit - Compute_shape_factor< depos_order-1 > compute_shape_factor_cell; - Compute_shape_factor_pair< depos_order > compute_shape_factors_node; + const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell; + const Compute_shape_factor_pair< depos_order > compute_shape_factors_node; double dxp_seg, dyp_seg, dzp_seg; double x0_new, y0_new, z0_new; double x0_old = x_old; @@ -1854,7 +1854,7 @@ void doVillasenorDepositionShapeNImplicit (const amrex::ParticleReal * const xp_ const int k0_cell = compute_shape_factor_cell( sz_cell, z0_bar-0.5 ); if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights - Compute_shape_factor_pair compute_shape_factors_cell; + const Compute_shape_factor_pair compute_shape_factors_cell; double sx_old_cell[depos_order] = {0.}; double sx_new_cell[depos_order] = {0.}; double sy_old_cell[depos_order] = {0.}; @@ -1939,13 +1939,13 @@ void doVillasenorDepositionShapeNImplicit (const amrex::ParticleReal * const xp_ // compute cell crossings in X-direction const auto i_old = static_cast(x_old-shift); const auto i_new = static_cast(x_new-shift); - int cell_crossings_x = std::abs(i_new-i_old); + const int cell_crossings_x = std::abs(i_new-i_old); num_segments += cell_crossings_x; // compute cell crossings in Z-direction const auto k_old = static_cast(z_old-shift); const auto k_new = static_cast(z_new-shift); - int cell_crossings_z = std::abs(k_new-k_old); + const int cell_crossings_z = std::abs(k_new-k_old); num_segments += cell_crossings_z; // need to assert that the number of cell crossings in each direction @@ -1965,8 +1965,8 @@ void doVillasenorDepositionShapeNImplicit (const amrex::ParticleReal * const xp_ } // loop over the number of segments and deposit - Compute_shape_factor< depos_order-1 > compute_shape_factor_cell; - Compute_shape_factor_pair< depos_order > compute_shape_factors_node; + const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell; + const Compute_shape_factor_pair< depos_order > compute_shape_factors_node; double dxp_seg, dzp_seg; double x0_new, z0_new; double x0_old = x_old; @@ -2015,7 +2015,7 @@ void doVillasenorDepositionShapeNImplicit (const amrex::ParticleReal * const xp_ const int k0_cell = compute_shape_factor_cell( sz_cell, z0_bar-0.5 ); if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights - Compute_shape_factor_pair compute_shape_factors_cell; + const Compute_shape_factor_pair compute_shape_factors_cell; double sx_old_cell[depos_order] = {0.}; double sx_new_cell[depos_order] = {0.}; double sz_old_cell[depos_order] = {0.}; @@ -2112,7 +2112,7 @@ void doVillasenorDepositionShapeNImplicit (const amrex::ParticleReal * const xp_ // compute cell crossings in Z-direction const auto k_old = static_cast(z_old-shift); const auto k_new = static_cast(z_new-shift); - int cell_crossings_z = std::abs(k_new-k_old); + const int cell_crossings_z = std::abs(k_new-k_old); num_segments += cell_crossings_z; // need to assert that the number of cell crossings in each direction @@ -2125,8 +2125,8 @@ void doVillasenorDepositionShapeNImplicit (const amrex::ParticleReal * const xp_ double Zcell = static_cast(k_old) + shift + 0.5*(1.-dirZ_sign); // loop over the number of segments and deposit - Compute_shape_factor< depos_order-1 > compute_shape_factor_cell; - Compute_shape_factor_pair< depos_order > compute_shape_factors_node; + const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell; + const Compute_shape_factor_pair< depos_order > compute_shape_factors_node; double dzp_seg; double z0_new; double z0_old = z_old; @@ -2152,7 +2152,7 @@ void doVillasenorDepositionShapeNImplicit (const amrex::ParticleReal * const xp_ const int k0_cell = compute_shape_factor_cell( sz_cell, z0_bar-0.5 ); if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights - Compute_shape_factor_pair compute_shape_factors_cell; + const Compute_shape_factor_pair compute_shape_factors_cell; double sz_old_cell[depos_order] = {0.}; double sz_new_cell[depos_order] = {0.}; const int k0_cell_2 = compute_shape_factors_cell( sz_old_cell, sz_new_cell, z0_old-0.5, z0_new-0.5 ); diff --git a/Source/Particles/Filter/FilterFunctors.H b/Source/Particles/Filter/FilterFunctors.H index 354d57e2c98..d1206eb022c 100644 --- a/Source/Particles/Filter/FilterFunctors.H +++ b/Source/Particles/Filter/FilterFunctors.H @@ -165,9 +165,9 @@ struct GeometryFilter bool operator () (const SuperParticleType& p, const amrex::RandomEngine&) const noexcept { if ( !m_is_active ) { return true; } - return ! (AMREX_D_TERM( (p.pos(0) < m_domain.lo(0)) || (p.pos(0) > m_domain.hi(0) ), - || (p.pos(1) < m_domain.lo(1)) || (p.pos(1) > m_domain.hi(1) ), - || (p.pos(2) < m_domain.lo(2)) || (p.pos(2) > m_domain.hi(2) ))); + return AMREX_D_TERM( (p.pos(0) >= m_domain.lo(0)) && (p.pos(0) <= m_domain.hi(0) ), + && (p.pos(1) >= m_domain.lo(1)) && (p.pos(1) <= m_domain.hi(1) ), + && (p.pos(2) >= m_domain.lo(2)) && (p.pos(2) <= m_domain.hi(2) )); } private: /** Whether this diagnostics is activated. Select all particles if false */ diff --git a/Source/Particles/Gather/FieldGather.H b/Source/Particles/Gather/FieldGather.H index d43c6c0741d..e51cfa41c58 100644 --- a/Source/Particles/Gather/FieldGather.H +++ b/Source/Particles/Gather/FieldGather.H @@ -525,12 +525,12 @@ void doGatherShapeNEsirkepovStencilImplicit ( #endif #if !defined(WARPX_DIM_1D_Z) - ParticleReal xp_np1 = 2._prt*xp_nph - xp_n; + const ParticleReal xp_np1 = 2._prt*xp_nph - xp_n; #endif #if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) - ParticleReal yp_np1 = 2._prt*yp_nph - yp_n; + const ParticleReal yp_np1 = 2._prt*yp_nph - yp_n; #endif - ParticleReal zp_np1 = 2._prt*zp_nph - zp_n; + const ParticleReal zp_np1 = 2._prt*zp_nph - zp_n; // computes current and old position in grid units #if defined(WARPX_DIM_RZ) @@ -610,8 +610,8 @@ void doGatherShapeNEsirkepovStencilImplicit ( // --- Compute shape factors // Compute shape factors for position as they are now and at old positions // [ijk]_new: leftmost grid point that the particle touches - Compute_shape_factor< depos_order > compute_shape_factor; - Compute_shifted_shape_factor< depos_order > compute_shifted_shape_factor; + const Compute_shape_factor< depos_order > compute_shape_factor; + const Compute_shifted_shape_factor< depos_order > compute_shifted_shape_factor; #if !defined(WARPX_DIM_1D_Z) const int i_E_new = compute_shape_factor(sx_E_new+1, x_new); @@ -661,8 +661,8 @@ void doGatherShapeNEsirkepovStencilImplicit ( #endif #if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - Compute_shape_factor< depos_order-1 > compute_shape_factor_By; - Compute_shifted_shape_factor< depos_order-1 > compute_shifted_shape_factor_By; + const Compute_shape_factor< depos_order-1 > compute_shape_factor_By; + const Compute_shifted_shape_factor< depos_order-1 > compute_shifted_shape_factor_By; const int i_By_new = compute_shape_factor_By(sx_By_new+1, x_new - 0.5_rt); const int i_By_old = compute_shifted_shape_factor_By(sx_By_old, x_old - 0.5_rt, i_By_new); const int k_By_new = compute_shape_factor_By(sz_By_new+1, z_new - 0.5_rt); @@ -679,7 +679,7 @@ void doGatherShapeNEsirkepovStencilImplicit ( for (int k=dkl_E; k<=depos_order+2-dku_E; k++) { for (int j=djl_E; j<=depos_order+2-dju_E; j++) { - amrex::Real sdzjk = one_third*(sy_E_new[j]*sz_E_new[k] + sy_E_old[j]*sz_E_old[k]) + const amrex::Real sdzjk = one_third*(sy_E_new[j]*sz_E_new[k] + sy_E_old[j]*sz_E_old[k]) +one_sixth*(sy_E_new[j]*sz_E_old[k] + sy_E_old[j]*sz_E_new[k]); amrex::Real sdxi = 0._rt; for (int i=dil_E; i<=depos_order+1-diu_E; i++) { @@ -691,7 +691,7 @@ void doGatherShapeNEsirkepovStencilImplicit ( } for (int k=dkl_E; k<=depos_order+2-dku_E; k++) { for (int i=dil_E; i<=depos_order+2-diu_E; i++) { - amrex::Real sdyik = one_third*(sx_E_new[i]*sz_E_new[k] + sx_E_old[i]*sz_E_old[k]) + const amrex::Real sdyik = one_third*(sx_E_new[i]*sz_E_new[k] + sx_E_old[i]*sz_E_old[k]) +one_sixth*(sx_E_new[i]*sz_E_old[k] + sx_E_old[i]*sz_E_new[k]); amrex::Real sdyj = 0._rt; for (int j=djl_E; j<=depos_order+1-dju_E; j++) { @@ -703,7 +703,7 @@ void doGatherShapeNEsirkepovStencilImplicit ( } for (int j=djl_E; j<=depos_order+2-dju_E; j++) { for (int i=dil_E; i<=depos_order+2-diu_E; i++) { - amrex::Real sdzij = one_third*(sx_E_new[i]*sy_E_new[j] + sx_E_old[i]*sy_E_old[j]) + const amrex::Real sdzij = one_third*(sx_E_new[i]*sy_E_new[j] + sx_E_old[i]*sy_E_old[j]) +one_sixth*(sx_E_new[i]*sy_E_old[j] + sx_E_old[i]*sy_E_new[j]); amrex::Real sdzk = 0._rt; for (int k=dkl_E; k<=depos_order+1-dku_E; k++) { @@ -715,7 +715,7 @@ void doGatherShapeNEsirkepovStencilImplicit ( } for (int k=dkl_B; k<=depos_order+2-dku_B; k++) { for (int j=djl_B; j<=depos_order+2-dju_B; j++) { - amrex::Real sdzjk = one_third*(sy_B_new[j]*sz_B_new[k] + sy_B_old[j]*sz_B_old[k]) + const amrex::Real sdzjk = one_third*(sy_B_new[j]*sz_B_new[k] + sy_B_old[j]*sz_B_old[k]) +one_sixth*(sy_B_new[j]*sz_B_old[k] + sy_B_old[j]*sz_B_new[k]); amrex::Real sdxi = 0._rt; for (int i=dil_B; i<=depos_order+1-diu_B; i++) { @@ -727,7 +727,7 @@ void doGatherShapeNEsirkepovStencilImplicit ( } for (int k=dkl_B; k<=depos_order+2-dku_B; k++) { for (int i=dil_B; i<=depos_order+2-diu_B; i++) { - amrex::Real sdyik = one_third*(sx_B_new[i]*sz_B_new[k] + sx_B_old[i]*sz_B_old[k]) + const amrex::Real sdyik = one_third*(sx_B_new[i]*sz_B_new[k] + sx_B_old[i]*sz_B_old[k]) +one_sixth*(sx_B_new[i]*sz_B_old[k] + sx_B_old[i]*sz_B_new[k]); amrex::Real sdyj = 0._rt; for (int j=djl_B; j<=depos_order+1-dju_B; j++) { @@ -739,7 +739,7 @@ void doGatherShapeNEsirkepovStencilImplicit ( } for (int j=djl_B; j<=depos_order+2-dju_B; j++) { for (int i=dil_B; i<=depos_order+2-diu_B; i++) { - amrex::Real sdzij = one_third*(sx_B_new[i]*sy_B_new[j] + sx_B_old[i]*sy_B_old[j]) + const amrex::Real sdzij = one_third*(sx_B_new[i]*sy_B_new[j] + sx_B_old[i]*sy_B_old[j]) +one_sixth*(sx_B_new[i]*sy_B_old[j] + sx_B_old[i]*sy_B_new[j]); amrex::Real sdzk = 0._rt; for (int k=dkl_B; k<=depos_order+1-dku_B; k++) { @@ -753,7 +753,7 @@ void doGatherShapeNEsirkepovStencilImplicit ( #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) for (int k=dkl_E; k<=depos_order+2-dku_E; k++) { - amrex::Real sdzk = 0.5_rt*(sz_E_new[k] + sz_E_old[k]); + const amrex::Real sdzk = 0.5_rt*(sz_E_new[k] + sz_E_old[k]); amrex::Real sdxi = 0._rt; for (int i=dil_E; i<=depos_order+1-diu_E; i++) { sdxi += (sx_E_old[i] - sx_E_new[i]); @@ -771,7 +771,7 @@ void doGatherShapeNEsirkepovStencilImplicit ( } } for (int i=dil_E; i<=depos_order+2-diu_E; i++) { - amrex::Real sdxi = 0.5_rt*(sx_E_new[i] + sx_E_old[i]); + const amrex::Real sdxi = 0.5_rt*(sx_E_new[i] + sx_E_old[i]); amrex::Real sdzk = 0._rt; for (int k=dkl_E; k<=depos_order+1-dku_E; k++) { sdzk += (sz_E_old[k] - sz_E_new[k]); @@ -797,7 +797,7 @@ void doGatherShapeNEsirkepovStencilImplicit ( // Gather field on particle Exp from field on grid ex_arr // Gather field on particle Bzp from field on grid bz_arr for (int k=dkl_E; k<=depos_order+2-dku_E; k++) { - amrex::Real sdzk = 0.5_rt*(sz_E_new[k] + sz_E_old[k]); + const amrex::Real sdzk = 0.5_rt*(sz_E_new[k] + sz_E_old[k]); amrex::Real sdxi = 0._rt; for (int i=dil_E; i<=depos_order+1-diu_E; i++) { sdxi += (sx_E_old[i] - sx_E_new[i]); @@ -824,7 +824,7 @@ void doGatherShapeNEsirkepovStencilImplicit ( // Gather field on particle Ezp from field on grid ez_arr // Gather field on particle Bxp from field on grid bx_arr for (int i=dil_E; i<=depos_order+2-diu_E; i++) { - amrex::Real sdxi = 0.5_rt*(sx_E_new[i] + sx_E_old[i]); + const amrex::Real sdxi = 0.5_rt*(sx_E_new[i] + sx_E_old[i]); amrex::Real sdzk = 0._rt; for (int k=dkl_E; k<=depos_order+1-dku_E; k++) { sdzk += (sz_E_old[k] - sz_E_new[k]); @@ -950,12 +950,12 @@ void doGatherPicnicShapeN ( Real const zmin = xyzmin[2]; #if !defined(WARPX_DIM_1D_Z) - ParticleReal xp_np1 = 2._prt*xp_nph - xp_n; + const ParticleReal xp_np1 = 2._prt*xp_nph - xp_n; #endif #if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) - ParticleReal yp_np1 = 2._prt*yp_nph - yp_n; + const ParticleReal yp_np1 = 2._prt*yp_nph - yp_n; #endif - ParticleReal zp_np1 = 2._prt*zp_nph - zp_n; + const ParticleReal zp_np1 = 2._prt*zp_nph - zp_n; #if !defined(WARPX_DIM_1D_Z) Real constexpr one_third = 1.0_rt / 3.0_rt; @@ -1019,19 +1019,19 @@ void doGatherPicnicShapeN ( // compute cell crossings in X-direction const auto i_old = static_cast(x_old-shift); const auto i_new = static_cast(x_new-shift); - int cell_crossings_x = std::abs(i_new-i_old); + const int cell_crossings_x = std::abs(i_new-i_old); num_segments += cell_crossings_x; // compute cell crossings in Y-direction const auto j_old = static_cast(y_old-shift); const auto j_new = static_cast(y_new-shift); - int cell_crossings_y = std::abs(j_new-j_old); + const int cell_crossings_y = std::abs(j_new-j_old); num_segments += cell_crossings_y; // compute cell crossings in Z-direction const auto k_old = static_cast(z_old-shift); const auto k_new = static_cast(z_new-shift); - int cell_crossings_z = std::abs(k_new-k_old); + const int cell_crossings_z = std::abs(k_new-k_old); num_segments += cell_crossings_z; // need to assert that the number of cell crossings in each direction @@ -1054,8 +1054,8 @@ void doGatherPicnicShapeN ( } // loop over the number of segments and deposit - Compute_shape_factor< depos_order-1 > compute_shape_factor_cell; - Compute_shape_factor_pair< depos_order > compute_shape_factors_node; + const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell; + const Compute_shape_factor_pair< depos_order > compute_shape_factors_node; double dxp_seg, dyp_seg, dzp_seg; double x0_new, y0_new, z0_new; double x0_old = x_old; @@ -1125,7 +1125,7 @@ void doGatherPicnicShapeN ( const int k0_cell = compute_shape_factor_cell( sz_cell, z0_bar-0.5 ); if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights - Compute_shape_factor_pair compute_shape_factors_cell; + const Compute_shape_factor_pair compute_shape_factors_cell; double sx_old_cell[depos_order] = {0.}; double sx_new_cell[depos_order] = {0.}; double sy_old_cell[depos_order] = {0.}; @@ -1205,7 +1205,7 @@ void doGatherPicnicShapeN ( // gather magnetic field const int depos_order_B = 1; // second template parameter? - Compute_shape_factor< depos_order_B > compute_shape_factor_B; + const Compute_shape_factor< depos_order_B > compute_shape_factor_B; double sz_bar_node[depos_order_B+1] = {0.}; double sz_bar_cell[depos_order_B+1] = {0.}; const int k_bar_node = compute_shape_factor_B(sz_bar_node, z_bar); @@ -1240,13 +1240,13 @@ void doGatherPicnicShapeN ( // compute cell crossings in X-direction const auto i_old = static_cast(x_old-shift); const auto i_new = static_cast(x_new-shift); - int cell_crossings_x = std::abs(i_new-i_old); + const int cell_crossings_x = std::abs(i_new-i_old); num_segments += cell_crossings_x; // compute cell crossings in Z-direction const auto k_old = static_cast(z_old-shift); const auto k_new = static_cast(z_new-shift); - int cell_crossings_z = std::abs(k_new-k_old); + const int cell_crossings_z = std::abs(k_new-k_old); num_segments += cell_crossings_z; // need to assert that the number of cell crossings in each direction @@ -1266,8 +1266,8 @@ void doGatherPicnicShapeN ( } // loop over the number of segments and deposit - Compute_shape_factor< depos_order-1 > compute_shape_factor_cell; - Compute_shape_factor_pair< depos_order > compute_shape_factors_node; + const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell; + const Compute_shape_factor_pair< depos_order > compute_shape_factors_node; double dxp_seg, dzp_seg; double x0_new, z0_new; double x0_old = x_old; @@ -1316,7 +1316,7 @@ void doGatherPicnicShapeN ( const int k0_cell = compute_shape_factor_cell(sz_cell, z0_bar-0.5); if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights - Compute_shape_factor_pair compute_shape_factors_cell; + const Compute_shape_factor_pair compute_shape_factors_cell; double sx_old_cell[depos_order] = {0.}; double sx_new_cell[depos_order] = {0.}; double sz_old_cell[depos_order] = {0.}; @@ -1404,7 +1404,7 @@ void doGatherPicnicShapeN ( // gather magnetic field and out-of-plane electric field const int depos_order_B = 1; // second template parameter? - Compute_shape_factor< depos_order_B > compute_shape_factor_B; + const Compute_shape_factor< depos_order_B > compute_shape_factor_B; double sz_bar_node[depos_order_B+1] = {0.}; double sz_bar_cell[depos_order_B+1] = {0.}; const int k_bar_node = compute_shape_factor_B(sz_bar_node, z_bar); @@ -1459,7 +1459,7 @@ void doGatherPicnicShapeN ( // compute cell crossings in Z-direction const auto k_old = static_cast(z_old-shift); const auto k_new = static_cast(z_new-shift); - int cell_crossings_z = std::abs(k_new-k_old); + const int cell_crossings_z = std::abs(k_new-k_old); num_segments += cell_crossings_z; // need to assert that the number of cell crossings in each direction @@ -1472,8 +1472,8 @@ void doGatherPicnicShapeN ( double Zcell = static_cast(k_old) + shift + 0.5*(1.-dirZ_sign); // loop over the number of segments and deposit - Compute_shape_factor< depos_order-1 > compute_shape_factor_cell; - Compute_shape_factor_pair< depos_order > compute_shape_factors_node; + const Compute_shape_factor< depos_order-1 > compute_shape_factor_cell; + const Compute_shape_factor_pair< depos_order > compute_shape_factors_node; double dzp_seg; double z0_new; double z0_old = z_old; @@ -1499,7 +1499,7 @@ void doGatherPicnicShapeN ( const int k0_cell = compute_shape_factor_cell( sz_cell, z0_bar-0.5 ); if constexpr (depos_order >= 3) { // higher-order correction to the cell-based weights - Compute_shape_factor_pair compute_shape_factors_cell; + const Compute_shape_factor_pair compute_shape_factors_cell; double sz_old_cell[depos_order] = {0.}; double sz_new_cell[depos_order] = {0.}; const int k0_cell_2 = compute_shape_factors_cell( sz_old_cell, sz_new_cell, z0_old-0.5, z0_new-0.5 ); @@ -1536,7 +1536,7 @@ void doGatherPicnicShapeN ( // gather magnetic field const int depos_order_B = 1; // second template parameter? - Compute_shape_factor< depos_order_B > compute_shape_factor_B; + const Compute_shape_factor< depos_order_B > compute_shape_factor_B; double sz_bar_node[depos_order_B+1] = {0.}; double sz_bar_cell[depos_order_B+1] = {0.}; const int k_bar_node = compute_shape_factor_B(sz_bar_node, z_bar); @@ -1591,8 +1591,8 @@ void doGatherShapeN(const GetParticlePosition& getPosition, const int n_rz_azimuthal_modes) { - amrex::GpuArray dx_arr = {dx[0], dx[1], dx[2]}; - amrex::GpuArray xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]}; + const amrex::GpuArray dx_arr = {dx[0], dx[1], dx[2]}; + const amrex::GpuArray xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]}; amrex::Array4 const& ex_arr = exfab->array(); amrex::Array4 const& ey_arr = eyfab->array(); diff --git a/Source/Particles/Gather/GetExternalFields.H b/Source/Particles/Gather/GetExternalFields.H index 9fd534e33d9..7000d6d7c26 100644 --- a/Source/Particles/Gather/GetExternalFields.H +++ b/Source/Particles/Gather/GetExternalFields.H @@ -142,7 +142,7 @@ struct GetExternalEBField // the plasma lens periods do not start before z=0 if (zl > 0) { // find which is the next lens - int i_lens = static_cast(std::floor(zl/m_repeated_plasma_lens_period)); + const auto i_lens = static_cast(std::floor(zl/m_repeated_plasma_lens_period)); if (i_lens < m_n_lenses) { amrex::ParticleReal const lens_start = m_repeated_plasma_lens_starts[i_lens] + i_lens*m_repeated_plasma_lens_period; amrex::ParticleReal const lens_end = lens_start + m_repeated_plasma_lens_lengths[i_lens]; diff --git a/Source/Particles/LaserParticleContainer.cpp b/Source/Particles/LaserParticleContainer.cpp index fd4e9397253..bf846c674f0 100644 --- a/Source/Particles/LaserParticleContainer.cpp +++ b/Source/Particles/LaserParticleContainer.cpp @@ -538,14 +538,14 @@ LaserParticleContainer::InitData (int lev) } } const auto np = static_cast(particle_z.size()); - amrex::Vector particle_ux(np, 0.0); - amrex::Vector particle_uy(np, 0.0); - amrex::Vector particle_uz(np, 0.0); + const amrex::Vector particle_ux(np, 0.0); + const amrex::Vector particle_uy(np, 0.0); + const amrex::Vector particle_uz(np, 0.0); if (Verbose()) { amrex::Print() << Utils::TextMsg::Info("Adding laser particles"); } amrex::Vector> attr; attr.push_back(particle_w); - amrex::Vector> attr_int; + const amrex::Vector> attr_int; // Add particles on level 0. They will be redistributed afterwards AddNParticles(0, np, particle_x, particle_y, particle_z, diff --git a/Source/Particles/ParticleCreation/DefaultInitialization.H b/Source/Particles/ParticleCreation/DefaultInitialization.H index 870fc82bd0f..3f0cbab79d5 100644 --- a/Source/Particles/ParticleCreation/DefaultInitialization.H +++ b/Source/Particles/ParticleCreation/DefaultInitialization.H @@ -161,7 +161,7 @@ void DefaultInitializeRuntimeAttributes (PTile& ptile, if constexpr (amrex::RunOnGpu>::value) { amrex::ParallelForRNG(stop - start, [=] AMREX_GPU_DEVICE (int i, amrex::RandomEngine const& engine) noexcept { - int ip = i + start; + const int ip = i + start; attr_ptr[ip] = quantum_sync_get_opt(engine); }); // Otherwise (e.g. particle tile allocated in pinned memory), run on CPU @@ -183,7 +183,7 @@ void DefaultInitializeRuntimeAttributes (PTile& ptile, if constexpr (amrex::RunOnGpu>::value) { amrex::ParallelForRNG(stop - start, [=] AMREX_GPU_DEVICE (int i, amrex::RandomEngine const& engine) noexcept { - int ip = i + start; + const int ip = i + start; attr_ptr[ip] = breit_wheeler_get_opt(engine); }); // Otherwise (e.g. particle tile allocated in pinned memory), run on CPU @@ -207,7 +207,7 @@ void DefaultInitializeRuntimeAttributes (PTile& ptile, if constexpr (amrex::RunOnGpu>::value) { amrex::ParallelFor(stop - start, [=] AMREX_GPU_DEVICE (int i) noexcept { - int ip = i + start; + const int ip = i + start; amrex::ParticleReal xp, yp, zp; get_position(ip, xp, yp, zp); attr_ptr[ip] = user_real_attrib_parserexec(xp, yp, zp, @@ -238,7 +238,7 @@ void DefaultInitializeRuntimeAttributes (PTile& ptile, if constexpr (amrex::RunOnGpu>::value) { amrex::ParallelFor(stop - start, [=] AMREX_GPU_DEVICE (int i) noexcept { - int ip = i + start; + const int ip = i + start; attr_ptr[ip] = ionization_initial_level; }); } else { @@ -259,7 +259,7 @@ void DefaultInitializeRuntimeAttributes (PTile& ptile, if constexpr (amrex::RunOnGpu>::value) { amrex::ParallelFor(stop - start, [=] AMREX_GPU_DEVICE (int i) noexcept { - int ip = i + start; + const int ip = i + start; amrex::ParticleReal xp, yp, zp; get_position(ip, xp, yp, zp); attr_ptr[ip] = static_cast( diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 3d5cbe0007d..1aa534c871b 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -355,16 +355,16 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp //if the species is not a lepton, do_classical_radiation_reaction //should be false WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - !(do_classical_radiation_reaction && - !(AmIA() || - AmIA() )), + (!do_classical_radiation_reaction) || + AmIA() || + AmIA(), "can't enable classical radiation reaction for non lepton species '" + species_name + "'."); //Only Boris pusher is compatible with radiation reaction WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - !(do_classical_radiation_reaction && - WarpX::particle_pusher_algo != ParticlePusherAlgo::Boris), + (!do_classical_radiation_reaction) || + WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris, "Radiation reaction can be enabled only if Boris pusher is used"); //_____________________________ @@ -570,26 +570,26 @@ PhysicalParticleContainer::AddGaussianBeam ( XDim3 u = plasma_injector.getMomentum(x, y, z); if (plasma_injector.do_focusing){ - XDim3 u_bulk = plasma_injector.getInjectorMomentumHost()->getBulkMomentum(x,y,z); - Real u_bulk_norm = std::sqrt( u_bulk.x*u_bulk.x+u_bulk.y*u_bulk.y+u_bulk.z*u_bulk.z ); + const XDim3 u_bulk = plasma_injector.getInjectorMomentumHost()->getBulkMomentum(x,y,z); + const Real u_bulk_norm = std::sqrt( u_bulk.x*u_bulk.x+u_bulk.y*u_bulk.y+u_bulk.z*u_bulk.z ); // Compute the position of the focal plane // (it is located at a distance `focal_distance` from the beam centroid, in the direction of the bulk velocity) - Real n_x = u_bulk.x/u_bulk_norm; - Real n_y = u_bulk.y/u_bulk_norm; - Real n_z = u_bulk.z/u_bulk_norm; - Real x_f = x_m + focal_distance * n_x; - Real y_f = y_m + focal_distance * n_y; - Real z_f = z_m + focal_distance * n_z; - Real gamma = std::sqrt( 1._rt + (u.x*u.x+u.y*u.y+u.z*u.z) ); - - Real v_x = u.x / gamma * PhysConst::c; - Real v_y = u.y / gamma * PhysConst::c; - Real v_z = u.z / gamma * PhysConst::c; + const Real n_x = u_bulk.x/u_bulk_norm; + const Real n_y = u_bulk.y/u_bulk_norm; + const Real n_z = u_bulk.z/u_bulk_norm; + const Real x_f = x_m + focal_distance * n_x; + const Real y_f = y_m + focal_distance * n_y; + const Real z_f = z_m + focal_distance * n_z; + const Real gamma = std::sqrt( 1._rt + (u.x*u.x+u.y*u.y+u.z*u.z) ); + + const Real v_x = u.x / gamma * PhysConst::c; + const Real v_y = u.y / gamma * PhysConst::c; + const Real v_z = u.z / gamma * PhysConst::c; // Compute the time at which the particle will cross the focal plane - Real v_dot_n = v_x * n_x + v_y * n_y + v_z * n_z; - Real t = ((x_f-x)*n_x + (y_f-y)*n_y + (z_f-z)*n_z) / v_dot_n; + const Real v_dot_n = v_x * n_x + v_y * n_y + v_z * n_z; + const Real t = ((x_f-x)*n_x + (y_f-y)*n_y + (z_f-z)*n_z) / v_dot_n; // Displace particles in the direction orthogonal to the beam bulk momentum // i.e. orthogonal to (n_x, n_y, n_z) @@ -672,18 +672,18 @@ PhysicalParticleContainer::AddGaussianBeam ( // Add the temporary CPU vectors to the particle structure auto const np = static_cast(particle_z.size()); - amrex::Vector xp(particle_x.data(), particle_x.data() + np); - amrex::Vector yp(particle_y.data(), particle_y.data() + np); - amrex::Vector zp(particle_z.data(), particle_z.data() + np); - amrex::Vector uxp(particle_ux.data(), particle_ux.data() + np); - amrex::Vector uyp(particle_uy.data(), particle_uy.data() + np); - amrex::Vector uzp(particle_uz.data(), particle_uz.data() + np); + const amrex::Vector xp(particle_x.data(), particle_x.data() + np); + const amrex::Vector yp(particle_y.data(), particle_y.data() + np); + const amrex::Vector zp(particle_z.data(), particle_z.data() + np); + const amrex::Vector uxp(particle_ux.data(), particle_ux.data() + np); + const amrex::Vector uyp(particle_uy.data(), particle_uy.data() + np); + const amrex::Vector uzp(particle_uz.data(), particle_uz.data() + np); amrex::Vector> attr; - amrex::Vector wp(particle_w.data(), particle_w.data() + np); + const amrex::Vector wp(particle_w.data(), particle_w.data() + np); attr.push_back(wp); - amrex::Vector> attr_int; + const amrex::Vector> attr_int; AddNParticles(0, np, xp, yp, zp, uxp, uyp, uzp, 1, attr, 0, attr_int, 1); @@ -796,18 +796,18 @@ PhysicalParticleContainer::AddPlasmaFromFile(PlasmaInjector & plasma_injector, } } // IO Processor auto const np = static_cast(particle_z.size()); - amrex::Vector xp(particle_x.data(), particle_x.data() + np); - amrex::Vector yp(particle_y.data(), particle_y.data() + np); - amrex::Vector zp(particle_z.data(), particle_z.data() + np); - amrex::Vector uxp(particle_ux.data(), particle_ux.data() + np); - amrex::Vector uyp(particle_uy.data(), particle_uy.data() + np); - amrex::Vector uzp(particle_uz.data(), particle_uz.data() + np); + const amrex::Vector xp(particle_x.data(), particle_x.data() + np); + const amrex::Vector yp(particle_y.data(), particle_y.data() + np); + const amrex::Vector zp(particle_z.data(), particle_z.data() + np); + const amrex::Vector uxp(particle_ux.data(), particle_ux.data() + np); + const amrex::Vector uyp(particle_uy.data(), particle_uy.data() + np); + const amrex::Vector uzp(particle_uz.data(), particle_uz.data() + np); amrex::Vector> attr; - amrex::Vector wp(particle_w.data(), particle_w.data() + np); + const amrex::Vector wp(particle_w.data(), particle_w.data() + np); attr.push_back(wp); - amrex::Vector> attr_int; + const amrex::Vector> attr_int; AddNParticles(0, np, xp, yp, zp, uxp, uyp, uzp, 1, attr, 0, attr_int, 1); @@ -880,14 +880,14 @@ PhysicalParticleContainer::AddParticles (int lev) plasma_injector->single_particle_u[1], plasma_injector->single_particle_u[2]); } - amrex::Vector xp = {plasma_injector->single_particle_pos[0]}; - amrex::Vector yp = {plasma_injector->single_particle_pos[1]}; - amrex::Vector zp = {plasma_injector->single_particle_pos[2]}; - amrex::Vector uxp = {plasma_injector->single_particle_u[0]}; - amrex::Vector uyp = {plasma_injector->single_particle_u[1]}; - amrex::Vector uzp = {plasma_injector->single_particle_u[2]}; - amrex::Vector> attr = {{plasma_injector->single_particle_weight}}; - amrex::Vector> attr_int; + const amrex::Vector xp = {plasma_injector->single_particle_pos[0]}; + const amrex::Vector yp = {plasma_injector->single_particle_pos[1]}; + const amrex::Vector zp = {plasma_injector->single_particle_pos[2]}; + const amrex::Vector uxp = {plasma_injector->single_particle_u[0]}; + const amrex::Vector uyp = {plasma_injector->single_particle_u[1]}; + const amrex::Vector uzp = {plasma_injector->single_particle_u[2]}; + const amrex::Vector> attr = {{plasma_injector->single_particle_weight}}; + const amrex::Vector> attr_int; AddNParticles(lev, 1, xp, yp, zp, uxp, uyp, uzp, 1, attr, 0, attr_int, 0); return; @@ -906,7 +906,7 @@ PhysicalParticleContainer::AddParticles (int lev) } amrex::Vector> attr; attr.push_back(plasma_injector->multiple_particles_weight); - amrex::Vector> attr_int; + const amrex::Vector> attr_int; AddNParticles(lev, static_cast(plasma_injector->multiple_particles_pos_x.size()), plasma_injector->multiple_particles_pos_x, plasma_injector->multiple_particles_pos_y, @@ -958,7 +958,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int const int num_ppc = plasma_injector.num_particles_per_cell; #ifdef WARPX_DIM_RZ - Real rmax = std::min(plasma_injector.xmax, part_realbox.hi(0)); + const Real rmax = std::min(plasma_injector.xmax, part_realbox.hi(0)); #endif const auto dx = geom.CellSizeArray(); @@ -996,7 +996,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int #ifdef WARPX_DIM_RZ const int nmodes = WarpX::n_rz_azimuthal_modes; - bool radially_weighted = plasma_injector.radially_weighted; + const bool radially_weighted = plasma_injector.radially_weighted; #endif MFItInfo info; @@ -1488,7 +1488,7 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, const amrex::Real num_ppc_real = plasma_injector.num_particles_per_cell_real; #ifdef WARPX_DIM_RZ - Real rmax = std::min(plasma_injector.xmax, geom.ProbDomain().hi(0)); + const Real rmax = std::min(plasma_injector.xmax, geom.ProbDomain().hi(0)); #endif const auto dx = geom.CellSizeArray(); @@ -1546,7 +1546,7 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, #ifdef WARPX_DIM_RZ const int nmodes = WarpX::n_rz_azimuthal_modes; const bool rz_random_theta = m_rz_random_theta; - bool radially_weighted = plasma_injector.radially_weighted; + const bool radially_weighted = plasma_injector.radially_weighted; #endif MFItInfo info; @@ -1864,7 +1864,7 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, Real const cos_theta = std::cos(theta); Real const sin_theta = std::sin(theta); // Rotate the position - amrex::Real radial_position = ppos.x; + 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) { @@ -1873,13 +1873,13 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, // the `inj_mom` objects generates a v*Gaussian distribution // along the Cartesian "x" direction by default. This // needs to be rotated along "r". - Real ur = pu.x; - Real ut = pu.y; + const Real ur = pu.x; + const Real ut = pu.y; pu.x = cos_theta*ur - sin_theta*ut; pu.y = sin_theta*ur + cos_theta*ut; } #endif - Real flux = inj_flux->getFlux(ppos.x, ppos.y, ppos.z, t); + const Real flux = inj_flux->getFlux(ppos.x, ppos.y, ppos.z, t); // Remove particle if flux is negative or 0 if (flux <= 0) { pa_idcpu[ip] = amrex::ParticleIdCpus::Invalid; @@ -2505,17 +2505,17 @@ PhysicalParticleContainer::SplitParticles (int lev) // they are not re-split when entering a higher level // AddNParticles calls Redistribute, so that particles // in pctmp_split are in the proper grids and tiles - amrex::Vector xp(psplit_x.data(), psplit_x.data() + np_split_to_add); - amrex::Vector yp(psplit_y.data(), psplit_y.data() + np_split_to_add); - amrex::Vector zp(psplit_z.data(), psplit_z.data() + np_split_to_add); - amrex::Vector uxp(psplit_ux.data(), psplit_ux.data() + np_split_to_add); - amrex::Vector uyp(psplit_uy.data(), psplit_uy.data() + np_split_to_add); - amrex::Vector uzp(psplit_uz.data(), psplit_uz.data() + np_split_to_add); - amrex::Vector wp(psplit_w.data(), psplit_w.data() + np_split_to_add); + const amrex::Vector xp(psplit_x.data(), psplit_x.data() + np_split_to_add); + const amrex::Vector yp(psplit_y.data(), psplit_y.data() + np_split_to_add); + const amrex::Vector zp(psplit_z.data(), psplit_z.data() + np_split_to_add); + const amrex::Vector uxp(psplit_ux.data(), psplit_ux.data() + np_split_to_add); + const amrex::Vector uyp(psplit_uy.data(), psplit_uy.data() + np_split_to_add); + const amrex::Vector uzp(psplit_uz.data(), psplit_uz.data() + np_split_to_add); + const amrex::Vector wp(psplit_w.data(), psplit_w.data() + np_split_to_add); amrex::Vector> attr; attr.push_back(wp); - amrex::Vector> attr_int; + const amrex::Vector> attr_int; pctmp_split.AddNParticles(lev, np_split_to_add, xp, @@ -3025,12 +3025,12 @@ PhysicalParticleContainer::ImplicitPushXP (WarpXParIter& pti, const Dim3 lo = lbound(box); - int depos_type = WarpX::current_deposition_algo; - int nox = WarpX::nox; - int n_rz_azimuthal_modes = WarpX::n_rz_azimuthal_modes; + const int depos_type = WarpX::current_deposition_algo; + const int nox = WarpX::nox; + const int n_rz_azimuthal_modes = WarpX::n_rz_azimuthal_modes; - amrex::GpuArray dx_arr = {dx[0], dx[1], dx[2]}; - amrex::GpuArray xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]}; + const amrex::GpuArray dx_arr = {dx[0], dx[1], dx[2]}; + const amrex::GpuArray xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]}; amrex::Array4 const& ex_arr = exfab->array(); amrex::Array4 const& ey_arr = eyfab->array(); @@ -3062,7 +3062,7 @@ PhysicalParticleContainer::ImplicitPushXP (WarpXParIter& pti, ParticleReal* uy_n = pti.GetAttribs(particle_comps["uy_n"]).dataPtr(); ParticleReal* uz_n = pti.GetAttribs(particle_comps["uz_n"]).dataPtr(); - int do_copy = (m_do_back_transformed_particles && (a_dt_type!=DtType::SecondHalf) ); + const int do_copy = (m_do_back_transformed_particles && (a_dt_type!=DtType::SecondHalf) ); CopyParticleAttribs copyAttribs; if (do_copy) { copyAttribs = CopyParticleAttribs(pti, tmp_particle_data, offset); @@ -3098,11 +3098,11 @@ PhysicalParticleContainer::ImplicitPushXP (WarpXParIter& pti, enum exteb_flags : int { no_exteb, has_exteb }; enum qed_flags : int { no_qed, has_qed }; - int exteb_runtime_flag = getExternalEB.isNoOp() ? no_exteb : has_exteb; + const int exteb_runtime_flag = getExternalEB.isNoOp() ? no_exteb : has_exteb; #ifdef WARPX_QED - int qed_runtime_flag = (local_has_quantum_sync || do_sync) ? has_qed : no_qed; + const int qed_runtime_flag = (local_has_quantum_sync || do_sync) ? has_qed : no_qed; #else - int qed_runtime_flag = no_qed; + const int qed_runtime_flag = no_qed; #endif // Using this version of ParallelFor with compile time options @@ -3118,20 +3118,20 @@ PhysicalParticleContainer::ImplicitPushXP (WarpXParIter& pti, // but uses the most recent velocity. #if (AMREX_SPACEDIM >= 2) amrex::ParticleReal xp = x_n[ip]; - amrex::ParticleReal xp_n = x_n[ip]; + const amrex::ParticleReal xp_n = x_n[ip]; #else amrex::ParticleReal xp = 0._rt; - amrex::ParticleReal xp_n = 0._rt; + const amrex::ParticleReal xp_n = 0._rt; #endif #if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) amrex::ParticleReal yp = y_n[ip]; - amrex::ParticleReal yp_n = y_n[ip]; + const amrex::ParticleReal yp_n = y_n[ip]; #else amrex::ParticleReal yp = 0._rt; - amrex::ParticleReal yp_n = 0._rt; + const amrex::ParticleReal yp_n = 0._rt; #endif amrex::ParticleReal zp = z_n[ip]; - amrex::ParticleReal zp_n = z_n[ip]; + const amrex::ParticleReal zp_n = z_n[ip]; UpdatePositionImplicit(xp, yp, zp, ux_n[ip], uy_n[ip], uz_n[ip], ux[ip], uy[ip], uz[ip], 0.5_rt*dt); setPosition(ip, xp, yp, zp); diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index d50851ffb1b..eed556bf998 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -73,7 +73,7 @@ void RigidInjectedParticleContainer::InitData() { // Perform Lorentz transform of `z_inject_plane` const amrex::Real t_boost = WarpX::GetInstance().gett_new(0); - amrex::Real zinject_plane_boost = zinject_plane/WarpX::gamma_boost + const amrex::Real zinject_plane_boost = zinject_plane/WarpX::gamma_boost - WarpX::beta_boost*PhysConst::c*t_boost; zinject_plane_levels.resize(finestLevel()+1, zinject_plane_boost); diff --git a/Source/Utils/SpeciesUtils.cpp b/Source/Utils/SpeciesUtils.cpp index f8515592a07..56c269c2efe 100644 --- a/Source/Utils/SpeciesUtils.cpp +++ b/Source/Utils/SpeciesUtils.cpp @@ -81,7 +81,7 @@ namespace SpeciesUtils { std::unique_ptr& h_inj_rho, std::unique_ptr& density_parser) { - amrex::ParmParse pp_species(species_name); + const amrex::ParmParse pp_species(species_name); // parse density information std::string rho_prof_s; @@ -126,7 +126,7 @@ namespace SpeciesUtils { { using namespace amrex::literals; - amrex::ParmParse pp_species(species_name); + const amrex::ParmParse pp_species(species_name); // parse momentum information std::string mom_dist_s; diff --git a/Source/Utils/WarpXMovingWindow.cpp b/Source/Utils/WarpXMovingWindow.cpp index fdcb7cb8da4..eb7eceaa2fd 100644 --- a/Source/Utils/WarpXMovingWindow.cpp +++ b/Source/Utils/WarpXMovingWindow.cpp @@ -88,7 +88,7 @@ WarpX::UpdateInjectionPosition (const amrex::Real a_dt) amrex::Real v_shift = 0._rt; if (plasma_injector != nullptr) { - amrex::XDim3 u_bulk = plasma_injector->getInjectorMomentumHost()->getBulkMomentum( + const amrex::XDim3 u_bulk = plasma_injector->getInjectorMomentumHost()->getBulkMomentum( current_injection_position[0], current_injection_position[1], current_injection_position[2]); diff --git a/Source/Utils/WarpXTagging.cpp b/Source/Utils/WarpXTagging.cpp index 192fc9e9152..f88380e2589 100644 --- a/Source/Utils/WarpXTagging.cpp +++ b/Source/Utils/WarpXTagging.cpp @@ -52,10 +52,10 @@ WarpX::ErrorEst (int lev, TagBoxArray& tags, Real /*time*/, int /*ngrow*/) #if defined (WARPX_DIM_3D) tag_val = (ref_parser(pos[0], pos[1], pos[2]) == 1); #elif defined (WARPX_DIM_XZ) || defined (WARPX_DIM_RZ) - amrex::Real unused = 0.0; + const auto unused = 0.0_rt; tag_val = (ref_parser(pos[0], unused, pos[1]) == 1); #elif defined (WARPX_DIM_1D_Z) - amrex::Real unused = 0.0; + const auto unused = 0.0_rt; tag_val = (ref_parser(unused, unused, pos[0]) == 1); #endif } else { diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index b13d7aaf63b..195a6e5fe10 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -1008,10 +1008,10 @@ WarpX::ReadParameters () if (maxLevel() > 0) { Vector lo, hi; - bool fine_tag_lo_specified = utils::parser::queryArrWithParser(pp_warpx, "fine_tag_lo", lo); - bool fine_tag_hi_specified = utils::parser::queryArrWithParser(pp_warpx, "fine_tag_hi", hi); + const bool fine_tag_lo_specified = utils::parser::queryArrWithParser(pp_warpx, "fine_tag_lo", lo); + const bool fine_tag_hi_specified = utils::parser::queryArrWithParser(pp_warpx, "fine_tag_hi", hi); std::string ref_patch_function; - bool parser_specified = pp_warpx.query("ref_patch_function(x,y,z)",ref_patch_function); + const bool parser_specified = pp_warpx.query("ref_patch_function(x,y,z)",ref_patch_function); WARPX_ALWAYS_ASSERT_WITH_MESSAGE( ((fine_tag_lo_specified && fine_tag_hi_specified) || parser_specified ), "For max_level > 0, you need to either set\ diff --git a/Source/ablastr/fields/PoissonSolver.H b/Source/ablastr/fields/PoissonSolver.H index eaeadfbb2fb..62dccc6a09d 100644 --- a/Source/ablastr/fields/PoissonSolver.H +++ b/Source/ablastr/fields/PoissonSolver.H @@ -170,7 +170,12 @@ computePhi (amrex::Vector const & rho, ); } +#if !(defined(AMREX_USE_EB) || defined(WARPX_DIM_RZ)) amrex::LPInfo info; +#else + const amrex::LPInfo info; +#endif + for (int lev=0; lev<=finest_level; lev++) { // Set the value of beta amrex::Array beta_solver = diff --git a/Source/ablastr/utils/SignalHandling.cpp b/Source/ablastr/utils/SignalHandling.cpp index 9521e5f9c3d..4095b9207ce 100644 --- a/Source/ablastr/utils/SignalHandling.cpp +++ b/Source/ablastr/utils/SignalHandling.cpp @@ -87,7 +87,7 @@ SignalHandling::parseSignalNameToNumber (const std::string &str) }; for (const auto& sp : signals_to_parse) { - std::string name_upper = sp.abbrev; + const std::string name_upper = sp.abbrev; std::string name_lower = name_upper; for (char &c : name_lower) { c = static_cast(std::tolower(static_cast(c))); diff --git a/Source/ablastr/utils/msg_logger/MsgLogger.cpp b/Source/ablastr/utils/msg_logger/MsgLogger.cpp index 629eab8c25d..dde041652da 100644 --- a/Source/ablastr/utils/msg_logger/MsgLogger.cpp +++ b/Source/ablastr/utils/msg_logger/MsgLogger.cpp @@ -225,6 +225,7 @@ std::vector Logger::get_msgs() const { auto res = std::vector{}; + res.reserve(m_messages.size()); for (const auto& msg_w_counter : m_messages) { res.emplace_back(msg_w_counter.first); } @@ -236,6 +237,7 @@ std::vector Logger::get_msgs_with_counter() const { auto res = std::vector{}; + res.reserve(m_messages.size()); for (const auto& msg : m_messages) { res.emplace_back(MsgWithCounter{msg.first, msg.second}); } @@ -309,6 +311,7 @@ std::vector Logger::one_rank_gather_msgs_with_counter_and_ranks() const { std::vector res; + res.reserve(m_messages.size()); for (const auto& el : m_messages) { res.emplace_back( From 766a521bea1a840b9ff25ef328dd971dfaa0f9f5 Mon Sep 17 00:00:00 2001 From: Roelof Groenewald <40245517+roelof-groenewald@users.noreply.github.com> Date: Tue, 5 Mar 2024 16:24:23 -0800 Subject: [PATCH 07/16] Enable `WarpX::ComputeEdgeLengths` in RZ (#4749) * Enable `WarpX::ComputeEdgeLengths` in RZ * Update exit message when using EBs in RZ with fully electromagnetic solvers. --- Source/EmbeddedBoundary/WarpXInitEB.cpp | 12 ++++-------- Source/WarpX.cpp | 5 +++-- 2 files changed, 7 insertions(+), 10 deletions(-) diff --git a/Source/EmbeddedBoundary/WarpXInitEB.cpp b/Source/EmbeddedBoundary/WarpXInitEB.cpp index 990ff1445ed..c3036b59aa4 100644 --- a/Source/EmbeddedBoundary/WarpXInitEB.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEB.cpp @@ -114,16 +114,15 @@ WarpX::InitEB () void WarpX::ComputeEdgeLengths (std::array< std::unique_ptr, 3 >& edge_lengths, const amrex::EBFArrayBoxFactory& eb_fact) { -#ifndef WARPX_DIM_RZ BL_PROFILE("ComputeEdgeLengths"); auto const &flags = eb_fact.getMultiEBCellFlagFab(); auto const &edge_centroid = eb_fact.getEdgeCent(); -#ifdef WARPX_DIM_XZ +#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) edge_lengths[1]->setVal(0.); #endif for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi){ -#ifdef WARPX_DIM_XZ +#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) for (int idim = 0; idim < 3; ++idim){ if(idim == 1) continue; #elif defined(WARPX_DIM_3D) @@ -148,7 +147,7 @@ WarpX::ComputeEdgeLengths (std::array< std::unique_ptr, 3 >& ed edge_lengths_dim(i, j, k) = 0.; }); } else { -#ifdef WARPX_DIM_XZ +#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) int idim_amrex = idim; if(idim == 2) idim_amrex = 1; auto const &edge_cent = edge_centroid[idim_amrex]->const_array(mfi); @@ -175,7 +174,6 @@ WarpX::ComputeEdgeLengths (std::array< std::unique_ptr, 3 >& ed } } } -#endif } @@ -246,11 +244,10 @@ WarpX::ComputeFaceAreas (std::array< std::unique_ptr, 3 >& face void WarpX::ScaleEdges (std::array< std::unique_ptr, 3 >& edge_lengths, const std::array& cell_size) { -#ifndef WARPX_DIM_RZ BL_PROFILE("ScaleEdges"); for (amrex::MFIter mfi(*edge_lengths[0]); mfi.isValid(); ++mfi) { -#ifdef WARPX_DIM_XZ +#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) for (int idim = 0; idim < 3; ++idim){ if(idim == 1) continue; #elif defined(WARPX_DIM_3D) @@ -267,7 +264,6 @@ WarpX::ScaleEdges (std::array< std::unique_ptr, 3 >& edge_lengt }); } } -#endif } void diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 195a6e5fe10..d453fa6f941 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -752,8 +752,9 @@ WarpX::ReadParameters () #if defined(AMREX_USE_EB) && defined(WARPX_DIM_RZ) WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - electromagnetic_solver_id==ElectromagneticSolverAlgo::None, - "Currently, the embedded boundary in RZ only works for electrostatic solvers (or no solver)."); + electromagnetic_solver_id==ElectromagneticSolverAlgo::None + || electromagnetic_solver_id==ElectromagneticSolverAlgo::HybridPIC, + "Currently, the embedded boundary in RZ only works for electrostatic solvers, the Ohm's law solver or with no solver installed."); #endif if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrame || From 3982b0005e82c030b2476dc5efd0ff65a18f6fa1 Mon Sep 17 00:00:00 2001 From: Eya D <81635404+EyaDammak@users.noreply.github.com> Date: Wed, 6 Mar 2024 16:25:25 -0800 Subject: [PATCH 08/16] New Python test: Particle-Boundary interaction (#4729) * enable the diagnostic of ParticleScraping in Python * Update picmi.py * Update picmi.py * new test * python update * modification of the script * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Update PICMI_inputs_rz.py * json * update * Update PICMI_inputs_rz.py * Update particle_boundary_interaction.json * Update PICMI_inputs_rz.py * Update PICMI_inputs_rz.py * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Update PICMI_inputs_rz.py * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Fix hanging script in parallel * Make the test executable * Update analysis script * Update particle_containers.py * Update PICMI_inputs_rz.py * Update analysis.py * Update analysis.py * Update particle_containers.py * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Remi Lehe Co-authored-by: Axel Huebl --- .../PICMI_inputs_rz.py | 169 ++++++++++++++++++ .../particle_boundary_interaction/analysis.py | 51 ++++++ Python/pywarpx/particle_containers.py | 3 + .../particle_boundary_interaction.json | 18 ++ Regression/WarpX-tests.ini | 21 +++ 5 files changed, 262 insertions(+) create mode 100644 Examples/Tests/particle_boundary_interaction/PICMI_inputs_rz.py create mode 100755 Examples/Tests/particle_boundary_interaction/analysis.py create mode 100644 Regression/Checksum/benchmarks_json/particle_boundary_interaction.json diff --git a/Examples/Tests/particle_boundary_interaction/PICMI_inputs_rz.py b/Examples/Tests/particle_boundary_interaction/PICMI_inputs_rz.py new file mode 100644 index 00000000000..4017a05d413 --- /dev/null +++ b/Examples/Tests/particle_boundary_interaction/PICMI_inputs_rz.py @@ -0,0 +1,169 @@ +#!/usr/bin/env python3 +# @Eya Dammak supervised by @Remi Lehe, 2024 +# --- Input file for particle-boundary interaction testing in RZ. +# --- This input is a simple case of reflection +# --- of one electron on the surface of a sphere. +import numpy as np + +from pywarpx import callbacks, particle_containers, picmi + +########################## +# numerics parameters +########################## + +dt = 1.0e-11 + +# --- Nb time steps + +max_steps = 23 +diagnostic_interval = 1 + +# --- grid + +nr = 64 +nz= 64 + +rmin = 0.0 +rmax = 2 +zmin = -2 +zmax = 2 + +########################## +# numerics components +########################## + +grid = picmi.CylindricalGrid( + number_of_cells = [nr, nz], + n_azimuthal_modes = 1, + lower_bound = [rmin, zmin], + upper_bound = [rmax, zmax], + lower_boundary_conditions = ['none', 'dirichlet'], + upper_boundary_conditions = ['dirichlet', 'dirichlet'], + lower_boundary_conditions_particles = ['absorbing', 'reflecting'], + upper_boundary_conditions_particles = ['absorbing', 'reflecting'] +) + + +solver = picmi.ElectrostaticSolver( + grid=grid, method='Multigrid', + warpx_absolute_tolerance=1e-7 +) + +embedded_boundary = picmi.EmbeddedBoundary( + implicit_function="-(x**2+y**2+z**2-radius**2)", + radius = 0.2 +) + +########################## +# physics components +########################## + +#one particle +e_dist=picmi.ParticleListDistribution(x=0.0, y=0.0, z=-0.25, ux=0.5e10, uy=0.0, uz=1.0e10, weight=1) + +electrons = picmi.Species( + particle_type='electron', name='electrons', initial_distribution=e_dist, warpx_save_particles_at_eb=1 +) + +########################## +# diagnostics +########################## + +field_diag = picmi.FieldDiagnostic( + name = 'diag1', + grid = grid, + period = diagnostic_interval, + data_list = ['Er', 'Ez', 'phi', 'rho','rho_electrons'], + warpx_format = 'openpmd', + write_dir = '.', + warpx_file_prefix = 'particle_boundary_interaction_plt' +) + +part_diag = picmi.ParticleDiagnostic(name = 'diag1', + period = diagnostic_interval, + species = [electrons], + warpx_format = 'openpmd', + write_dir = '.', + warpx_file_prefix = 'particle_boundary_interaction_plt' +) + +########################## +# simulation setup +########################## + +sim = picmi.Simulation( + solver=solver, + time_step_size = dt, + max_steps = max_steps, + warpx_embedded_boundary=embedded_boundary, + warpx_amrex_the_arena_is_managed=1, +) + +sim.add_species( + electrons, + layout = picmi.GriddedLayout( + n_macroparticle_per_cell=[10, 1, 1], grid=grid + ) +) +sim.add_diagnostic(part_diag) +sim.add_diagnostic(field_diag) + +sim.initialize_inputs() +sim.initialize_warpx() + +########################## +# python particle data access +########################## + +def concat( list_of_arrays ): + if len(list_of_arrays) == 0: + # Return a 1d array of size 0 + return np.empty(0) + else: + return np.concatenate( list_of_arrays ) + + +def mirror_reflection(): + buffer = particle_containers.ParticleBoundaryBufferWrapper() #boundary buffer + + #STEP 1: extract the different parameters of the boundary buffer (normal, time, position) + lev = 0 # level 0 (no mesh refinement here) + delta_t = concat(buffer.get_particle_boundary_buffer("electrons", 'eb', 'deltaTimeScraped', lev)) + r = concat(buffer.get_particle_boundary_buffer("electrons", 'eb', 'x', lev)) + theta = concat(buffer.get_particle_boundary_buffer("electrons", 'eb', 'theta', lev)) + z = concat(buffer.get_particle_boundary_buffer("electrons", 'eb', 'z', lev)) + x= r*np.cos(theta) #from RZ coordinates to 3D coordinates + y= r*np.sin(theta) + ux = concat(buffer.get_particle_boundary_buffer("electrons", 'eb', 'ux', lev)) + uy = concat(buffer.get_particle_boundary_buffer("electrons", 'eb', 'uy', lev)) + uz = concat(buffer.get_particle_boundary_buffer("electrons", 'eb', 'uz', lev)) + w = concat(buffer.get_particle_boundary_buffer("electrons", 'eb', 'w', lev)) + nx = concat(buffer.get_particle_boundary_buffer("electrons", 'eb', 'nx', lev)) + ny = concat(buffer.get_particle_boundary_buffer("electrons", 'eb', 'ny', lev)) + nz = concat(buffer.get_particle_boundary_buffer("electrons", 'eb', 'nz', lev)) + + #STEP 2: use these parameters to inject particle from the same position in the plasma + elect_pc = particle_containers.ParticleContainerWrapper('electrons') #general particle container + + ####this part is specific to the case of simple reflection. + un=ux*nx+uy*ny+uz*nz + ux_reflect=-2*un*nx+ux #for a "mirror reflection" u(sym)=-2(u.n)n+u + uy_reflect=-2*un*ny+uy + uz_reflect=-2*un*nz+uz + elect_pc.add_particles( + x=x + (dt-delta_t)*ux_reflect, y=y + (dt-delta_t)*uy_reflect, z=z + (dt-delta_t)*uz_reflect, + ux=ux_reflect, uy=uy_reflect, uz=uz_reflect, + w=w + ) #adds the particle in the general particle container at the next step + #### Can be modified depending on the model of interaction. + + buffer.clear_buffer() #reinitialise the boundary buffer + +callbacks.installafterstep(mirror_reflection) #mirror_reflection is called at the next step + # using the new particle container modified at the last step + +########################## +# simulation run +########################## + +sim.step(max_steps) #the whole process is done "max_steps" times diff --git a/Examples/Tests/particle_boundary_interaction/analysis.py b/Examples/Tests/particle_boundary_interaction/analysis.py new file mode 100755 index 00000000000..9768751ed74 --- /dev/null +++ b/Examples/Tests/particle_boundary_interaction/analysis.py @@ -0,0 +1,51 @@ +#!/usr/bin/env python +# @Eya Dammak supervised by @Remi Lehe, 2024 +""" +This script tests the last coordinate after adding an electron. +The sphere is centered on O and has a radius of 0.2 (EB) +The electron is initially at: (0,0,-0.25) and moves with a velocity: +(0.5e10,0,1.0e10) with a time step of 1e-11. +An input file PICMI_inputs_rz.py is used. +""" +import os +import sys + +import numpy as np +from openpmd_viewer import OpenPMDTimeSeries +import yt + +yt.funcs.mylog.setLevel(0) +sys.path.insert(1, '../../../../warpx/Regression/Checksum/') +import checksumAPI + +# Open plotfile specified in command line +filename = sys.argv[1] +test_name = os.path.split(os.getcwd())[1] +checksumAPI.evaluate_checksum(test_name, filename, output_format='openpmd') + +ts = OpenPMDTimeSeries('./particle_boundary_interaction_plt') + +it=ts.iterations +x,y,z=ts.get_particle(['x','y','z'], species='electrons', iteration=it[-1]) + +# Analytical results calculated +x_analytic=0.03532 +y_analytic=0.00000 +z_analytic=-0.20531 + +print('NUMERICAL coordinates of the point of contact:') +print('x=%5.5f, y=%5.5f, z=%5.5f' % (x[0], y[0], z[0])) +print('\n') +print('ANALYTICAL coordinates of the point of contact:') +print('x=%5.5f, y=%5.5f, z=%5.5f' % (x_analytic, y_analytic, z_analytic)) + +tolerance=1e-5 + +diff_x=np.abs((x[0]-x_analytic)/x_analytic) +diff_z=np.abs((z[0]-z_analytic)/z_analytic) + +print('\n') +print("percentage error for x = %5.4f %%" %(diff_x *100)) +print("percentage error for z = %5.4f %%" %(diff_z *100)) + +assert (diff_x < tolerance) and (y[0] < 1e-8) and (diff_z < tolerance), 'Test particle_boundary_interaction did not pass' diff --git a/Python/pywarpx/particle_containers.py b/Python/pywarpx/particle_containers.py index 19b0a1ab6c4..f8e687f99bd 100644 --- a/Python/pywarpx/particle_containers.py +++ b/Python/pywarpx/particle_containers.py @@ -761,6 +761,9 @@ def get_particle_boundary_buffer(self, species_name, boundary, comp_name, level) The data for the arrays are not copied, but share the underlying memory buffer with WarpX. The arrays are fully writeable. + You can find `here https://github.com/ECP-WarpX/WarpX/blob/319e55b10ad4f7c71b84a4fb21afbafe1f5b65c2/Examples/Tests/particle_boundary_interaction/PICMI_inputs_rz.py` + an example of a simple case of particle-boundary interaction (reflection). + Parameters ---------- diff --git a/Regression/Checksum/benchmarks_json/particle_boundary_interaction.json b/Regression/Checksum/benchmarks_json/particle_boundary_interaction.json new file mode 100644 index 00000000000..e35516758f2 --- /dev/null +++ b/Regression/Checksum/benchmarks_json/particle_boundary_interaction.json @@ -0,0 +1,18 @@ +{ + "electrons": { + "particle_momentum_x": 7.082238783412152e-21, + "particle_momentum_y": 0.0, + "particle_momentum_z": 7.319015172219235e-21, + "particle_position_x": 0.03532003602713794, + "particle_position_y": 0.0, + "particle_position_z": 0.20531131195378569, + "particle_weight": 1.0 + }, + "lev=0": { + "Er": 1.3287679891086577e-06, + "Ez": 1.791737740110229e-06, + "phi": 3.3983929077087634e-07, + "rho": 7.811529217958718e-16, + "rho_electrons": 7.811529217958718e-16 + } +} diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 426ecb0717c..ab422d5631a 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -4681,3 +4681,24 @@ doVis = 0 compareParticles = 0 doComparison = 0 analysisRoutine = Examples/Tests/gaussian_beam/analysis_focusing_beam.py + +[particle_boundary_interaction] +buildDir = . +inputFile = Examples/Tests/particle_boundary_interaction/PICMI_inputs_rz.py +runtime_params = +customRunCmd = python3 PICMI_inputs_rz.py +dim = 2 +addToCompileString = USE_PYTHON_MAIN=TRUE USE_RZ=TRUE +cmakeSetupOpts = -DWarpX_DIMS="RZ" -DWarpX_EB=ON -DWarpX_PYTHON=ON +target = pip_install +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 1 +compileTest = 0 +doVis = 0 +compareParticles = 1 +particleTypes = electrons +outputFile = particle_boundary_interaction_plt +analysisRoutine = Examples/Tests/particle_boundary_interaction/analysis.py From 642aa5480ec6afc036e1e5133a8cb65e7739feca Mon Sep 17 00:00:00 2001 From: Eya D <81635404+EyaDammak@users.noreply.github.com> Date: Wed, 6 Mar 2024 16:32:11 -0800 Subject: [PATCH 09/16] Adding normal components to regular boundary buffer (#4742) * first draft * adding normal only * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Update ParticleBoundaryBuffer.cpp * Update ParticleBoundaryBuffer.cpp --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- Source/Particles/ParticleBoundaryBuffer.cpp | 23 +++++++++++++++++---- 1 file changed, 19 insertions(+), 4 deletions(-) diff --git a/Source/Particles/ParticleBoundaryBuffer.cpp b/Source/Particles/ParticleBoundaryBuffer.cpp index f991f211d28..d1f9f814bc5 100644 --- a/Source/Particles/ParticleBoundaryBuffer.cpp +++ b/Source/Particles/ParticleBoundaryBuffer.cpp @@ -1,5 +1,5 @@ /* Copyright 2021 Andrew Myers - * + * modified by Remi Lehe, Eya Dammak 2023 * This file is part of WarpX. * * License: BSD-3-Clause-LBNL @@ -22,6 +22,7 @@ #include #include #include +#include using namespace amrex::literals; struct IsOutsideDomainBoundary { @@ -164,8 +165,12 @@ struct FindEmbeddedBoundaryIntersection { struct CopyAndTimestamp { int m_step_index; int m_delta_index; + int m_normal_index; int m_step; const amrex::Real m_dt; + int m_idim; + int m_iside; + template AMREX_GPU_HOST_DEVICE @@ -182,8 +187,16 @@ struct CopyAndTimestamp { for (int j = 0; j < src.m_num_runtime_int; ++j) { dst.m_runtime_idata[j][dst_i] = src.m_runtime_idata[j][src_i]; } + dst.m_runtime_idata[m_step_index][dst_i] = m_step; dst.m_runtime_rdata[m_delta_index][dst_i] = 0._rt; //delta_fraction is initialized to zero + + //calculation of the normal to the boundary + std::array n = {0.0, 0.0, 0.0}; + n[m_idim]=1-2*m_iside; + dst.m_runtime_rdata[m_normal_index][dst_i]= n[0]; + dst.m_runtime_rdata[m_normal_index+1][dst_i]= n[1]; + dst.m_runtime_rdata[m_normal_index+2][dst_i]= n[2]; } }; @@ -351,7 +364,6 @@ void ParticleBoundaryBuffer::gatherParticles (MultiParticleContainer& mypc, const amrex::Geometry& geom = warpx_instance.Geom(0); auto plo = geom.ProbLoArray(); auto phi = geom.ProbHiArray(); - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { if (geom.isPeriodic(idim)) { continue; } @@ -367,6 +379,9 @@ void ParticleBoundaryBuffer::gatherParticles (MultiParticleContainer& mypc, buffer[i] = pc.make_alike(); buffer[i].AddIntComp("stepScraped", false); buffer[i].AddRealComp("deltaTimeScraped", false); + buffer[i].AddRealComp("nx", false); + buffer[i].AddRealComp("ny", false); + buffer[i].AddRealComp("nz", false); } auto& species_buffer = buffer[i]; for (int lev = 0; lev < pc.numLevels(); ++lev) @@ -412,11 +427,11 @@ void ParticleBoundaryBuffer::gatherParticles (MultiParticleContainer& mypc, const int step_scraped_index = string_to_index_intcomp.at("stepScraped"); auto string_to_index_realcomp = buffer[i].getParticleRuntimeComps(); const int delta_index = string_to_index_realcomp.at("deltaTimeScraped"); + const int normal_index = string_to_index_realcomp.at("nx"); const int step = warpx_instance.getistep(0); - amrex::filterAndTransformParticles(ptile_buffer, ptile, predicate, - CopyAndTimestamp{step_scraped_index, delta_index, step, dt}, + CopyAndTimestamp{step_scraped_index, delta_index, normal_index, step, dt, idim, iside}, 0, dst_index); } } From a551793464237297f766a694579ec17f5354edec Mon Sep 17 00:00:00 2001 From: Roelof Groenewald <40245517+roelof-groenewald@users.noreply.github.com> Date: Wed, 6 Mar 2024 16:42:11 -0800 Subject: [PATCH 10/16] Add function to set domain boundary potentials from Python (#4740) * add function to set domain boundary potentials from Python * switch function arguments to form `potential_[lo/hi]_[x/y/z]` and add to docs --- Docs/source/usage/workflows/python_extend.rst | 6 +++++ Source/Python/WarpX.cpp | 22 +++++++++++++++++++ 2 files changed, 28 insertions(+) diff --git a/Docs/source/usage/workflows/python_extend.rst b/Docs/source/usage/workflows/python_extend.rst index d3049508da6..47610e0d7ba 100644 --- a/Docs/source/usage/workflows/python_extend.rst +++ b/Docs/source/usage/workflows/python_extend.rst @@ -90,6 +90,12 @@ This object is the Python equivalent to the C++ ``WarpX`` simulation class. .. py:method:: get_particle_boundary_buffer + .. py:method:: set_potential_on_domain_boundary(potential_[lo/hi]_[x/y/z]: str) + + The potential on the domain boundaries can be modified when using the electrostatic solver. + This function updates the strings and function parsers which set the domain + boundary potentials during the Poisson solve. + .. py:method:: set_potential_on_eb(potential: str) The embedded boundary (EB) conditions can be modified when using the electrostatic solver. diff --git a/Source/Python/WarpX.cpp b/Source/Python/WarpX.cpp index e42b8268fe1..4b0ad56a0ce 100644 --- a/Source/Python/WarpX.cpp +++ b/Source/Python/WarpX.cpp @@ -170,6 +170,28 @@ The physical fields in WarpX have the following naming: "Get the current physical time step size on mesh-refinement level ``lev``." ) + .def("set_potential_on_domain_boundary", + [](WarpX& wx, + std::string potential_lo_x, std::string potential_hi_x, + std::string potential_lo_y, std::string potential_hi_y, + std::string potential_lo_z, std::string potential_hi_z) + { + if (potential_lo_x != "") wx.m_poisson_boundary_handler.potential_xlo_str = potential_lo_x; + if (potential_hi_x != "") wx.m_poisson_boundary_handler.potential_xhi_str = potential_hi_x; + if (potential_lo_y != "") wx.m_poisson_boundary_handler.potential_ylo_str = potential_lo_y; + if (potential_hi_y != "") wx.m_poisson_boundary_handler.potential_yhi_str = potential_hi_y; + if (potential_lo_z != "") wx.m_poisson_boundary_handler.potential_zlo_str = potential_lo_z; + if (potential_hi_z != "") wx.m_poisson_boundary_handler.potential_zhi_str = potential_hi_z; + wx.m_poisson_boundary_handler.buildParsers(); + }, + py::arg("potential_lo_x") = "", + py::arg("potential_hi_x") = "", + py::arg("potential_lo_y") = "", + py::arg("potential_hi_y") = "", + py::arg("potential_lo_z") = "", + py::arg("potential_hi_z") = "", + "Sets the domain boundary potential string(s) and updates the function parser." + ) .def("set_potential_on_eb", [](WarpX& wx, std::string potential) { wx.m_poisson_boundary_handler.setPotentialEB(potential); From 0c6c452fc3b23e7e1b4ddc3c8293f53f98146cf4 Mon Sep 17 00:00:00 2001 From: Arianna Formenti Date: Wed, 6 Mar 2024 16:57:50 -0800 Subject: [PATCH 11/16] clean up `ablastr/fields` (#4753) * move PoissonInterpCPtoFP to Interpolate.H * concatenate nested namespaces --- Source/ablastr/fields/Interpolate.H | 55 +++++++++++++++++++++ Source/ablastr/fields/PoissonSolver.H | 36 +------------- Source/ablastr/fields/VectorPoissonSolver.H | 2 +- 3 files changed, 57 insertions(+), 36 deletions(-) create mode 100644 Source/ablastr/fields/Interpolate.H diff --git a/Source/ablastr/fields/Interpolate.H b/Source/ablastr/fields/Interpolate.H new file mode 100644 index 00000000000..a9f0a7fc75f --- /dev/null +++ b/Source/ablastr/fields/Interpolate.H @@ -0,0 +1,55 @@ +/* Copyright 2019-2022 Axel Huebl, Remi Lehe + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef ABLASTR_INTERPOLATE_H +#define ABLASTR_INTERPOLATE_H + +#include +#include +#include + +#include + +#include +#include + + +namespace ablastr::fields::details { + + /** Local interpolation from phi_cp to phi[lev+1] + * + * This is needed to work-around an NVCC limitation in downstream code (ImpactX), + * when nesting lambdas. Otherwise this could be written directly into the + * ParallelFor. + * + * @param[out] phi_fp_arr phi on the fine level + * @param[in] phi_cp_arr phi on the coarse level + * @param[in] refratio refinement ration + */ + struct PoissonInterpCPtoFP + { + PoissonInterpCPtoFP( + amrex::Array4 const phi_fp_arr, + amrex::Array4 const phi_cp_arr, + amrex::IntVect const refratio) + : m_phi_fp_arr(phi_fp_arr), m_phi_cp_arr(phi_cp_arr), m_refratio(refratio) + {} + + AMREX_GPU_DEVICE AMREX_FORCE_INLINE + void + operator() (int i, int j, int k) const noexcept + { + amrex::mf_nodebilin_interp(i, j, k, 0, m_phi_fp_arr, 0, m_phi_cp_arr, + 0, m_refratio); + } + + amrex::Array4 const m_phi_fp_arr; + amrex::Array4 const m_phi_cp_arr; + amrex::IntVect const m_refratio; + }; +} // namespace ablastr::fields::details + +#endif // ABLASTR_INTERPOLATE_H diff --git a/Source/ablastr/fields/PoissonSolver.H b/Source/ablastr/fields/PoissonSolver.H index 62dccc6a09d..9c6053bed07 100644 --- a/Source/ablastr/fields/PoissonSolver.H +++ b/Source/ablastr/fields/PoissonSolver.H @@ -11,6 +11,7 @@ #include #include #include +#include #include @@ -52,41 +53,6 @@ namespace ablastr::fields { -namespace details -{ - /** Local interpolation from phi_cp to phi[lev+1] - * - * This is needed to work-around an NVCC limitation in downstream code (ImpactX), - * when nesting lambdas. Otherwise this could be written directly into the - * ParallelFor. - * - * @param[out] phi_fp_arr phi on the fine level - * @param[in] phi_cp_arr phi on the coarse level - * @param[in] refratio refinement ration - */ - struct PoissonInterpCPtoFP - { - PoissonInterpCPtoFP( - amrex::Array4 const phi_fp_arr, - amrex::Array4 const phi_cp_arr, - amrex::IntVect const refratio) - : m_phi_fp_arr(phi_fp_arr), m_phi_cp_arr(phi_cp_arr), m_refratio(refratio) - {} - - AMREX_GPU_DEVICE AMREX_FORCE_INLINE - void - operator() (int i, int j, int k) const noexcept - { - amrex::mf_nodebilin_interp(i, j, k, 0, m_phi_fp_arr, 0, m_phi_cp_arr, - 0, m_refratio); - } - - amrex::Array4 const m_phi_fp_arr; - amrex::Array4 const m_phi_cp_arr; - amrex::IntVect const m_refratio; - }; -} - /** Compute the potential `phi` by solving the Poisson equation * * Uses `rho` as a source, assuming that the source moves at a diff --git a/Source/ablastr/fields/VectorPoissonSolver.H b/Source/ablastr/fields/VectorPoissonSolver.H index 8f9021456eb..6e68d0fadb5 100644 --- a/Source/ablastr/fields/VectorPoissonSolver.H +++ b/Source/ablastr/fields/VectorPoissonSolver.H @@ -8,7 +8,7 @@ #define ABLASTR_VECTOR_POISSON_SOLVER_H #include -#include +#include #include #include #include From 5e0d0ec716e9f32aa9e88dd49773ff96c73ededb Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Fri, 8 Mar 2024 06:52:01 +0100 Subject: [PATCH 12/16] Split clang-tidy CI test into 4 to improve performances (#4747) * split clang-tidy checks to improve performances * rename folders and tests * fix concurrency * Simplify --------- Co-authored-by: Axel Huebl --- .github/workflows/clang_tidy.yml | 25 ++++++++++++++----------- 1 file changed, 14 insertions(+), 11 deletions(-) diff --git a/.github/workflows/clang_tidy.yml b/.github/workflows/clang_tidy.yml index c166941f6e2..5abb42352f4 100644 --- a/.github/workflows/clang_tidy.yml +++ b/.github/workflows/clang_tidy.yml @@ -8,7 +8,10 @@ concurrency: jobs: run_clang_tidy: - name: clang-tidy + strategy: + matrix: + dim: [1, 2, RZ, 3] + name: clang-tidy-${{ matrix.dim }}D runs-on: ubuntu-22.04 if: github.event.pull_request.draft == false steps: @@ -35,16 +38,16 @@ jobs: export CXX=$(which clang++-15) export CC=$(which clang-15) - cmake -S . -B build_clang_tidy \ - -DCMAKE_VERBOSE_MAKEFILE=ON \ - -DWarpX_DIMS="1;2;RZ;3" \ - -DWarpX_MPI=ON \ - -DWarpX_COMPUTE=OMP \ - -DWarpX_PSATD=ON \ - -DWarpX_QED=ON \ - -DWarpX_QED_TABLE_GEN=ON \ - -DWarpX_OPENPMD=ON \ - -DWarpX_PRECISION=SINGLE \ + cmake -S . -B build_clang_tidy \ + -DCMAKE_VERBOSE_MAKEFILE=ON \ + -DWarpX_DIMS="${{ matrix.dim }}" \ + -DWarpX_MPI=ON \ + -DWarpX_COMPUTE=OMP \ + -DWarpX_PSATD=ON \ + -DWarpX_QED=ON \ + -DWarpX_QED_TABLE_GEN=ON \ + -DWarpX_OPENPMD=ON \ + -DWarpX_PRECISION=SINGLE \ -DCMAKE_CXX_COMPILER_LAUNCHER=ccache cmake --build build_clang_tidy -j 4 From 510f016771fbfdad49d50277dfd5ec57f3ffc230 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 8 Mar 2024 01:16:54 -0800 Subject: [PATCH 13/16] Replace links to learn git (#4758) * Replace links to learn git --- CONTRIBUTING.rst | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/CONTRIBUTING.rst b/CONTRIBUTING.rst index a3a6a3046a9..77b8200b0d5 100644 --- a/CONTRIBUTING.rst +++ b/CONTRIBUTING.rst @@ -10,10 +10,7 @@ Git workflow ------------ The WarpX project uses `git `_ for version control. -If you are new to git, you can follow one of these tutorials: - -- `Learn git with bitbucket `_ -- `git - the simple guide `_ +If you are new to git, you can follow `this tutorial `__. Configure your GitHub Account & Development Machine ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ From 887a16733611f088ebeb8821304f1735a5a05405 Mon Sep 17 00:00:00 2001 From: Roelof Groenewald <40245517+roelof-groenewald@users.noreply.github.com> Date: Fri, 8 Mar 2024 06:00:53 -0800 Subject: [PATCH 14/16] Bugfix in `fields.py` for GPU run without `cupy` (#4750) * Bugfix in `fields.py` for GPU run without `cupy` * apply suggestion from code review --- Python/pywarpx/fields.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/Python/pywarpx/fields.py b/Python/pywarpx/fields.py index dc0de9d9491..0f680595ef4 100644 --- a/Python/pywarpx/fields.py +++ b/Python/pywarpx/fields.py @@ -47,6 +47,7 @@ import cupy as cp except ImportError: cp = None +from .LoadThirdParty import load_cupy try: from mpi4py import MPI as mpi @@ -516,6 +517,12 @@ def __setitem__(self, index, value): if not isinstance(ic , slice) or len(global_shape) < 4: global_shape[3:3] = [1] value3d.shape = global_shape + if libwarpx.libwarpx_so.Config.have_gpu: + # check if cupy is available for use + xp, cupy_status = load_cupy() + if cupy_status is not None: + libwarpx.amr.Print(cupy_status) + starts = [ixstart, iystart, izstart] stops = [ixstop, iystop, izstop] for mfi in self.mf: @@ -526,7 +533,7 @@ def __setitem__(self, index, value): slice_value = value3d[global_slices] if libwarpx.libwarpx_so.Config.have_gpu: # Copy data from host to device - slice_value = cp.asarray(slice_value) + slice_value = xp.asarray(slice_value) mf_arr[block_slices] = slice_value else: mf_arr[block_slices] = value From ebbe63496274f1a65a1dae3c2daf56bf1d6c1d87 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Sat, 9 Mar 2024 01:06:24 +0900 Subject: [PATCH 15/16] Release 24.03 (#4759) * AMReX: 24.03 * pyAMReX: 24.03 * WarpX: 24.03 --- .github/workflows/cuda.yml | 2 +- CMakeLists.txt | 2 +- Docs/source/conf.py | 4 ++-- Python/setup.py | 2 +- Regression/WarpX-GPU-tests.ini | 2 +- Regression/WarpX-tests.ini | 2 +- cmake/dependencies/AMReX.cmake | 4 ++-- cmake/dependencies/pyAMReX.cmake | 4 ++-- run_test.sh | 2 +- setup.py | 2 +- 10 files changed, 13 insertions(+), 13 deletions(-) diff --git a/.github/workflows/cuda.yml b/.github/workflows/cuda.yml index 4b23ef931db..5e3a3077f5d 100644 --- a/.github/workflows/cuda.yml +++ b/.github/workflows/cuda.yml @@ -115,7 +115,7 @@ jobs: which nvcc || echo "nvcc not in PATH!" git clone https://github.com/AMReX-Codes/amrex.git ../amrex - cd ../amrex && git checkout --detach 3525b4a3f27eb64f746dd69b6613f71bb02d6e63 && cd - + cd ../amrex && git checkout --detach 24.03 && cd - make COMP=gcc QED=FALSE USE_MPI=TRUE USE_GPU=TRUE USE_OMP=FALSE USE_PSATD=TRUE USE_CCACHE=TRUE -j 4 ccache -s diff --git a/CMakeLists.txt b/CMakeLists.txt index 3a947b01dcd..4ddc685c733 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,7 +1,7 @@ # Preamble #################################################################### # cmake_minimum_required(VERSION 3.20.0) -project(WarpX VERSION 24.02) +project(WarpX VERSION 24.03) include(${WarpX_SOURCE_DIR}/cmake/WarpXFunctions.cmake) diff --git a/Docs/source/conf.py b/Docs/source/conf.py index b34c437b829..1a439c8b76e 100644 --- a/Docs/source/conf.py +++ b/Docs/source/conf.py @@ -103,9 +103,9 @@ def __init__(self, *args, **kwargs): # built documents. # # The short X.Y version. -version = u'24.02' +version = u'24.03' # The full version, including alpha/beta/rc tags. -release = u'24.02' +release = u'24.03' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/Python/setup.py b/Python/setup.py index 3c6c0f605c0..08e616cc279 100644 --- a/Python/setup.py +++ b/Python/setup.py @@ -54,7 +54,7 @@ package_data = {} setup(name = 'pywarpx', - version = '24.02', + version = '24.03', packages = ['pywarpx'], package_dir = {'pywarpx': 'pywarpx'}, description = """Wrapper of WarpX""", diff --git a/Regression/WarpX-GPU-tests.ini b/Regression/WarpX-GPU-tests.ini index 40338dc1480..2f0c7b607bb 100644 --- a/Regression/WarpX-GPU-tests.ini +++ b/Regression/WarpX-GPU-tests.ini @@ -60,7 +60,7 @@ emailBody = Check https://ccse.lbl.gov/pub/GpuRegressionTesting/WarpX/ for more [AMReX] dir = /home/regtester/git/amrex/ -branch = 3525b4a3f27eb64f746dd69b6613f71bb02d6e63 +branch = 24.03 [source] dir = /home/regtester/git/WarpX diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index ab422d5631a..1f88db26d7a 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -59,7 +59,7 @@ emailBody = Check https://ccse.lbl.gov/pub/RegressionTesting/WarpX/ for more det [AMReX] dir = /home/regtester/AMReX_RegTesting/amrex/ -branch = 3525b4a3f27eb64f746dd69b6613f71bb02d6e63 +branch = 24.03 [source] dir = /home/regtester/AMReX_RegTesting/warpx diff --git a/cmake/dependencies/AMReX.cmake b/cmake/dependencies/AMReX.cmake index 5b225de6d50..a1045457ec0 100644 --- a/cmake/dependencies/AMReX.cmake +++ b/cmake/dependencies/AMReX.cmake @@ -250,7 +250,7 @@ macro(find_amrex) endif() set(COMPONENT_PRECISION ${WarpX_PRECISION} P${WarpX_PARTICLE_PRECISION}) - find_package(AMReX 24.02 CONFIG REQUIRED COMPONENTS ${COMPONENT_ASCENT} ${COMPONENT_DIMS} ${COMPONENT_EB} PARTICLES ${COMPONENT_PIC} ${COMPONENT_PRECISION} ${COMPONENT_SENSEI} LSOLVERS) + find_package(AMReX 24.03 CONFIG REQUIRED COMPONENTS ${COMPONENT_ASCENT} ${COMPONENT_DIMS} ${COMPONENT_EB} PARTICLES ${COMPONENT_PIC} ${COMPONENT_PRECISION} ${COMPONENT_SENSEI} LSOLVERS) # note: TINYP skipped because user-configured and optional # AMReX CMake helper scripts @@ -273,7 +273,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 "3525b4a3f27eb64f746dd69b6613f71bb02d6e63" +set(WarpX_amrex_branch "24.03" 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 e8748f2e920..6360da84276 100644 --- a/cmake/dependencies/pyAMReX.cmake +++ b/cmake/dependencies/pyAMReX.cmake @@ -64,7 +64,7 @@ function(find_pyamrex) endif() elseif(NOT WarpX_pyamrex_internal) # TODO: MPI control - find_package(pyAMReX 24.02 CONFIG REQUIRED) + find_package(pyAMReX 24.03 CONFIG REQUIRED) message(STATUS "pyAMReX: Found version '${pyAMReX_VERSION}'") endif() endfunction() @@ -79,7 +79,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 "bfb599fd8361f8ef0765c487fd7bb69409bf78af" +set(WarpX_pyamrex_branch "24.03" CACHE STRING "Repository branch for WarpX_pyamrex_repo if(WarpX_pyamrex_internal)") diff --git a/run_test.sh b/run_test.sh index b5cdc627e35..f556dd37e76 100755 --- a/run_test.sh +++ b/run_test.sh @@ -68,7 +68,7 @@ python3 -m pip install --upgrade -r warpx/Regression/requirements.txt # Clone AMReX and warpx-data git clone https://github.com/AMReX-Codes/amrex.git -cd amrex && git checkout --detach 3525b4a3f27eb64f746dd69b6613f71bb02d6e63 && cd - +cd amrex && git checkout --detach 24.03 && cd - # warpx-data contains various required data sets git clone --depth 1 https://github.com/ECP-WarpX/warpx-data.git # openPMD-example-datasets contains various required data sets diff --git a/setup.py b/setup.py index 197a39ce23f..b0f919d707a 100644 --- a/setup.py +++ b/setup.py @@ -278,7 +278,7 @@ def build_extension(self, ext): setup( name='pywarpx', # note PEP-440 syntax: x.y.zaN but x.y.z.devN - version = '24.02', + version = '24.03', packages = ['pywarpx'], package_dir = {'pywarpx': 'Python/pywarpx'}, author='Jean-Luc Vay, David P. Grote, Maxence Thévenet, Rémi Lehe, Andrew Myers, Weiqun Zhang, Axel Huebl, et al.', From e38acc4da155087b3cdccb502fc808639c0c5c30 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 8 Mar 2024 18:12:04 -0800 Subject: [PATCH 16/16] Implement stair-case Yee solver with EB in RZ geometry (#2707) * Allow compilation with RZ EB * Do not push cells for RZ Yee solver, when covered with EB * Fix compilation errors * Fix additional compilation errors * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Fix additional compilation errors * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Add automated test * Add automated test * Fix path in tests * Enable parser in RZ * Update example script * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Clean-up PR * Initialize EB quantities * Modified EM field initialization in 2D with EB * Typo fix * Typo fix * Ignoring unused variables correctly * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Correct condition for updating E * Correct update of B * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Update B push * Update input script * Revert "Update input script" This reverts commit 5087485db606f77806c08849c51af9be635ca622. * Update initialization * Updated test * Move test to a different folder * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Add test for WarpX-test.ini * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Fix path for tests * Update test description * Update test metadata * Add checksum file * Revert changes * Revert changes * Change lx to lr * Revert "Change lx to lr" This reverts commit be3039af229d5dee509c8de0a712addf099a2b82. * Change lx to lr --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: lgiacome --- .../analysis_fields.py | 41 ++++++++++++++++++ .../embedded_boundary_diffraction/inputs_rz | 42 +++++++++++++++++++ .../EmbeddedBoundaryDiffraction.json | 10 +++++ Regression/WarpX-tests.ini | 18 ++++++++ Source/EmbeddedBoundary/WarpXInitEB.cpp | 16 +++---- .../FiniteDifferenceSolver/EvolveB.cpp | 4 +- .../FiniteDifferenceSolver/EvolveE.cpp | 27 ++++++++++-- .../FiniteDifferenceSolver.H | 1 + Source/WarpX.cpp | 7 ---- 9 files changed, 144 insertions(+), 22 deletions(-) create mode 100755 Examples/Tests/embedded_boundary_diffraction/analysis_fields.py create mode 100644 Examples/Tests/embedded_boundary_diffraction/inputs_rz create mode 100644 Regression/Checksum/benchmarks_json/EmbeddedBoundaryDiffraction.json diff --git a/Examples/Tests/embedded_boundary_diffraction/analysis_fields.py b/Examples/Tests/embedded_boundary_diffraction/analysis_fields.py new file mode 100755 index 00000000000..da344f332a1 --- /dev/null +++ b/Examples/Tests/embedded_boundary_diffraction/analysis_fields.py @@ -0,0 +1,41 @@ +#!/usr/bin/env python3 +""" +This test checks the implementation of the embedded boundary in cylindrical geometry, +by checking the diffraction of a laser by an embedded boundary here. +We then check that the first minimum of the diffracted intensity pattern +occurs along the angle given by the theoretical Airy pattern, i.e. +theta_diffraction = 1.22 * lambda / d +""" +import os +import sys + +import numpy as np +from openpmd_viewer import OpenPMDTimeSeries +from scipy.ndimage import gaussian_filter1d + +sys.path.insert(1, '../../../../warpx/Regression/Checksum/') +import checksumAPI + +ts = OpenPMDTimeSeries('./EmbeddedBoundaryDiffraction_plt/') + +# Extract the intensity as a function of r and z +Ex, info = ts.get_field('E', 'x', iteration=300) +I = gaussian_filter1d(Ex**2, sigma=5, axis=0) # Extract intensity by averaging E^2 over wavelength +irmax = np.argmax( I, axis=-1) + +# Find the radius of the first minimum, as a function of z +def r_first_minimum(iz): + ir = len(info.r)//2 + while I[iz, ir+1] < I[iz, ir]: + ir += 1 + return info.r[ir] +r = np.array([ r_first_minimum(iz) for iz in range(len(info.z)) ]) + +# Check that this corresponds to the prediction from the Airy pattern +theta_diffraction = np.arcsin(1.22*0.1/0.4)/2 +assert np.all( abs(r[50:] - theta_diffraction*info.z[50:]) < 0.03 ) + +# Open the right plot file +filename = sys.argv[1] +test_name = os.path.split(os.getcwd())[1] +checksumAPI.evaluate_checksum(test_name, filename, output_format='openpmd') diff --git a/Examples/Tests/embedded_boundary_diffraction/inputs_rz b/Examples/Tests/embedded_boundary_diffraction/inputs_rz new file mode 100644 index 00000000000..705c9641dfb --- /dev/null +++ b/Examples/Tests/embedded_boundary_diffraction/inputs_rz @@ -0,0 +1,42 @@ +# This script tests the diffraction of a laser by +# a cylindrical object, represented by an embedded boundary here + +max_step = 300 +amr.n_cell = 128 256 +amr.max_grid_size = 256 +amr.max_level = 0 + +geometry.dims = RZ +geometry.prob_lo = 0. -0.2 +geometry.prob_hi = 2 1.4 +warpx.cfl = 0.99 +algo.particle_shape = 1 + +boundary.field_lo = none absorbing_silver_mueller +boundary.field_hi = pec absorbing_silver_mueller + +# Make the cylindrical object that the laser will diffract on +my_constants.aperture_l = 0.01 +my_constants.aperture_r = 0.2 +warpx.eb_implicit_function = "if( (abs(z)<0.5*aperture_l) and (x, 3 >& ed void WarpX::ComputeFaceAreas (std::array< std::unique_ptr, 3 >& face_areas, const amrex::EBFArrayBoxFactory& eb_fact) { -#ifndef WARPX_DIM_RZ BL_PROFILE("ComputeFaceAreas"); auto const &flags = eb_fact.getMultiEBCellFlagFab(); -#ifdef WARPX_DIM_XZ +#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) //In 2D the volume frac is actually the area frac. auto const &area_frac = eb_fact.getVolFrac(); #elif defined(WARPX_DIM_3D) @@ -194,12 +193,12 @@ WarpX::ComputeFaceAreas (std::array< std::unique_ptr, 3 >& face "ComputeFaceAreas: Only implemented in 2D3V and 3D3V"); #endif -#ifdef WARPX_DIM_XZ +#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) face_areas[0]->setVal(0.); face_areas[2]->setVal(0.); #endif for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi) { -#ifdef WARPX_DIM_XZ +#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) // In 2D we change the extrema of the for loop so that we only have the case idim=1 for (int idim = 1; idim < AMREX_SPACEDIM; ++idim) { #elif defined(WARPX_DIM_3D) @@ -223,7 +222,7 @@ WarpX::ComputeFaceAreas (std::array< std::unique_ptr, 3 >& face face_areas_dim(i, j, k) = amrex::Real(0.); }); } else { -#ifdef WARPX_DIM_XZ +#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) auto const &face = area_frac.const_array(mfi); #elif defined(WARPX_DIM_3D) auto const &face = area_frac[idim]->const_array(mfi); @@ -237,7 +236,6 @@ WarpX::ComputeFaceAreas (std::array< std::unique_ptr, 3 >& face } } } -#endif } @@ -269,13 +267,12 @@ WarpX::ScaleEdges (std::array< std::unique_ptr, 3 >& edge_lengt void WarpX::ScaleAreas(std::array< std::unique_ptr, 3 >& face_areas, const std::array& cell_size) { -#ifndef WARPX_DIM_RZ BL_PROFILE("ScaleAreas"); amrex::Real full_area; for (amrex::MFIter mfi(*face_areas[0]); mfi.isValid(); ++mfi) { -#ifdef WARPX_DIM_XZ +#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) // In 2D we change the extrema of the for loop so that we only have the case idim=1 for (int idim = 1; idim < AMREX_SPACEDIM; ++idim) { #elif defined(WARPX_DIM_3D) @@ -286,7 +283,7 @@ WarpX::ScaleAreas(std::array< std::unique_ptr, 3 >& face_areas, #endif const amrex::Box& box = mfi.tilebox(face_areas[idim]->ixType().toIntVect(), face_areas[idim]->nGrowVect() ); -#ifdef WARPX_DIM_XZ +#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) full_area = cell_size[0]*cell_size[2]; #elif defined(WARPX_DIM_3D) if (idim == 0) { @@ -308,7 +305,6 @@ WarpX::ScaleAreas(std::array< std::unique_ptr, 3 >& face_areas, } } -#endif } diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp index 4425e655917..fbc1397b413 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp @@ -59,8 +59,8 @@ void FiniteDifferenceSolver::EvolveB ( std::array< std::unique_ptr >, 3 >& borrowing, int lev, amrex::Real const dt ) { -#ifndef AMREX_USE_EB - amrex::ignore_unused(area_mod, ECTRhofield, Venl, flag_info_cell, borrowing); +#if defined(WARPX_DIM_RZ) || !defined(AMREX_USE_EB) + amrex::ignore_unused(area_mod, ECTRhofield, Venl, flag_info_cell, borrowing); #endif // Select algorithm (The choice of algorithm is a runtime option, diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp index e52cca846bb..5f707fbc927 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp @@ -67,8 +67,7 @@ void FiniteDifferenceSolver::EvolveE ( // but we compile code for each algorithm, using templates) #ifdef WARPX_DIM_RZ if (m_fdtd_algo == ElectromagneticSolverAlgo::Yee){ - ignore_unused(edge_lengths); - EvolveECylindrical ( Efield, Bfield, Jfield, Ffield, lev, dt ); + EvolveECylindrical ( Efield, Bfield, Jfield, edge_lengths, Ffield, lev, dt ); #else if (m_grid_type == GridType::Collocated) { @@ -181,7 +180,6 @@ void FiniteDifferenceSolver::EvolveECartesian ( }, [=] AMREX_GPU_DEVICE (int i, int j, int k){ - #ifdef AMREX_USE_EB // Skip field push if this cell is fully covered by embedded boundaries if (lz(i,j,k) <= 0) return; @@ -235,9 +233,14 @@ void FiniteDifferenceSolver::EvolveECylindrical ( std::array< std::unique_ptr, 3 >& Efield, std::array< std::unique_ptr, 3 > const& Bfield, std::array< std::unique_ptr, 3 > const& Jfield, + std::array< std::unique_ptr, 3 > const& edge_lengths, std::unique_ptr const& Ffield, int lev, amrex::Real const dt ) { +#ifndef AMREX_USE_EB + amrex::ignore_unused(edge_lengths); +#endif + amrex::LayoutData* cost = WarpX::getCosts(lev); // Loop through the grids, and over the tiles within each grid @@ -262,6 +265,11 @@ void FiniteDifferenceSolver::EvolveECylindrical ( Array4 const& jt = Jfield[1]->array(mfi); Array4 const& jz = Jfield[2]->array(mfi); +#ifdef AMREX_USE_EB + amrex::Array4 const& lr = edge_lengths[0]->array(mfi); + amrex::Array4 const& lz = edge_lengths[2]->array(mfi); +#endif + // Extract stencil coefficients Real const * const AMREX_RESTRICT coefs_r = m_stencil_coefs_r.dataPtr(); auto const n_coefs_r = static_cast(m_stencil_coefs_r.size()); @@ -284,6 +292,10 @@ void FiniteDifferenceSolver::EvolveECylindrical ( amrex::ParallelFor(ter, tet, tez, [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/){ +#ifdef AMREX_USE_EB + // Skip field push if this cell is fully covered by embedded boundaries + if (lr(i, j, 0) <= 0) return; +#endif Real const r = rmin + (i + 0.5_rt)*dr; // r on cell-centered point (Er is cell-centered in r) Er(i, j, 0, 0) += c2 * dt*( - T_Algo::DownwardDz(Bt, coefs_z, n_coefs_z, i, j, 0, 0) @@ -301,6 +313,11 @@ void FiniteDifferenceSolver::EvolveECylindrical ( }, [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/){ +#ifdef AMREX_USE_EB + // Skip field push if this cell is fully covered by embedded boundaries + // The Et field is at a node, so we need to check if the node is covered + if (lr(i, j, 0)<=0 || lr(i-1, j, 0)<=0 || lz(i, j-1, 0)<=0 || lz(i, j, 0)<=0) return; +#endif Real const r = rmin + i*dr; // r on a nodal grid (Et is nodal in r) if (r != 0) { // Off-axis, regular Maxwell equations Et(i, j, 0, 0) += c2 * dt*( @@ -342,6 +359,10 @@ void FiniteDifferenceSolver::EvolveECylindrical ( }, [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/){ +#ifdef AMREX_USE_EB + // Skip field push if this cell is fully covered by embedded boundaries + if (lz(i, j, 0) <= 0) return; +#endif Real const r = rmin + i*dr; // r on a nodal grid (Ez is nodal in r) if (r != 0) { // Off-axis, regular Maxwell equations Ez(i, j, 0, 0) += c2 * dt*( diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H index 861b2648c1e..d045a30dd44 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H @@ -213,6 +213,7 @@ class FiniteDifferenceSolver std::array< std::unique_ptr, 3 >& Efield, std::array< std::unique_ptr, 3 > const& Bfield, std::array< std::unique_ptr, 3 > const& Jfield, + std::array< std::unique_ptr, 3 > const& edge_lengths, std::unique_ptr const& Ffield, int lev, amrex::Real dt ); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index d453fa6f941..6c3186c5c8d 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -750,13 +750,6 @@ WarpX::ReadParameters () electromagnetic_solver_id = ElectromagneticSolverAlgo::None; } -#if defined(AMREX_USE_EB) && defined(WARPX_DIM_RZ) - WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - electromagnetic_solver_id==ElectromagneticSolverAlgo::None - || electromagnetic_solver_id==ElectromagneticSolverAlgo::HybridPIC, - "Currently, the embedded boundary in RZ only works for electrostatic solvers, the Ohm's law solver or with no solver installed."); -#endif - if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrame || electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic) {