From f1af6f40636f8ae3085d350de425ce52f06c5479 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Mon, 22 Jul 2024 11:38:32 -0700 Subject: [PATCH 01/12] AMReX/pyAMReX/PICSAR: Weekly Update (#5074) * 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 e6ba115925a..c2ec9aaa763 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 dcb9cc0383dcc71e38dee9070574e325a812f8bf && cd - + cd ../amrex && git checkout --detach 0c3273f5e591815909180f8ffaf5b793cabbf9bc && cd - make COMP=gcc QED=FALSE USE_MPI=TRUE USE_GPU=TRUE USE_OMP=FALSE USE_FFT=TRUE USE_CCACHE=TRUE -j 4 ccache -s diff --git a/Regression/WarpX-GPU-tests.ini b/Regression/WarpX-GPU-tests.ini index b32ca61c709..a2447ac3cf5 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 = dcb9cc0383dcc71e38dee9070574e325a812f8bf +branch = 0c3273f5e591815909180f8ffaf5b793cabbf9bc [source] dir = /home/regtester/git/WarpX diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 6d7c57f0961..17ecfd64f29 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 = dcb9cc0383dcc71e38dee9070574e325a812f8bf +branch = 0c3273f5e591815909180f8ffaf5b793cabbf9bc [source] dir = /home/regtester/AMReX_RegTesting/warpx diff --git a/cmake/dependencies/AMReX.cmake b/cmake/dependencies/AMReX.cmake index 8acbd974779..025b4d3b6e0 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 "dcb9cc0383dcc71e38dee9070574e325a812f8bf" +set(WarpX_amrex_branch "0c3273f5e591815909180f8ffaf5b793cabbf9bc" 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 a2fbd0ddee8..a4598953432 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 "18f0026b1dd9b2aa4869c96f74e4b77262a067d0" +set(WarpX_pyamrex_branch "ff4643869c63d4ee40a87054b901f61eefcb97a3" CACHE STRING "Repository branch for WarpX_pyamrex_repo if(WarpX_pyamrex_internal)") diff --git a/run_test.sh b/run_test.sh index c01c3fb1176..1692f66263c 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 dcb9cc0383dcc71e38dee9070574e325a812f8bf && cd - +cd amrex && git checkout --detach 0c3273f5e591815909180f8ffaf5b793cabbf9bc && 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 396c06106ac7d635c70d46e112fb1b7ff2b720ad Mon Sep 17 00:00:00 2001 From: David Grote Date: Wed, 24 Jul 2024 23:05:44 -0700 Subject: [PATCH 02/12] Fix again the collisionXYZ benchmark (#5083) --- .../Tests/collision/analysis_collision_2d.py | 2 +- .../benchmarks_json/collisionXYZ.json | 28 +++++++++---------- .../Checksum/benchmarks_json/collisionXZ.json | 16 +++++++++++ 3 files changed, 31 insertions(+), 15 deletions(-) diff --git a/Examples/Tests/collision/analysis_collision_2d.py b/Examples/Tests/collision/analysis_collision_2d.py index bd013898a8a..29482909b2e 100755 --- a/Examples/Tests/collision/analysis_collision_2d.py +++ b/Examples/Tests/collision/analysis_collision_2d.py @@ -114,4 +114,4 @@ dim, species_name) test_name = os.path.split(os.getcwd())[1] -checksumAPI.evaluate_checksum(test_name, fn, do_particles=False) +checksumAPI.evaluate_checksum(test_name, fn) diff --git a/Regression/Checksum/benchmarks_json/collisionXYZ.json b/Regression/Checksum/benchmarks_json/collisionXYZ.json index 39df1b9eee1..2aa0e4686eb 100644 --- a/Regression/Checksum/benchmarks_json/collisionXYZ.json +++ b/Regression/Checksum/benchmarks_json/collisionXYZ.json @@ -6,25 +6,25 @@ "Ex": 0.0, "Ey": 0.0, "Ez": 0.0, - "T_electron": 353604.6247926339, - "T_ion": 347976.6168136309 + "T_electron": 381488.6528070591, + "T_ion": 320091.2785835478 }, "electron": { - "particle_momentum_x": 8.370755929299189e-19, - "particle_momentum_y": 8.228112213603589e-19, - "particle_momentum_z": 8.204295817378347e-19, - "particle_position_x": 21284971.94721422, - "particle_position_y": 21212829.42991966, - "particle_position_z": 21214774.536558084, + "particle_momentum_x": 8.667025573698235e-19, + "particle_momentum_y": 8.457499789250831e-19, + "particle_momentum_z": 8.482438182280524e-19, + "particle_position_x": 21262567.138872623, + "particle_position_y": 21245135.070665065, + "particle_position_z": 21232644.283726066, "particle_weight": 7.168263344048695e+28 }, "ion": { - "particle_momentum_x": 2.0074097598289766e-18, - "particle_momentum_y": 1.8203553942782305e-18, - "particle_momentum_z": 1.823420185235695e-18, - "particle_position_x": 21227192.857240494, - "particle_position_y": 21286501.692027714, - "particle_position_z": 21245587.6706009, + "particle_momentum_x": 1.9300495097720012e-18, + "particle_momentum_y": 1.747257416857836e-18, + "particle_momentum_z": 1.7510296287537058e-18, + "particle_position_x": 21217348.883301035, + "particle_position_y": 21300859.0630925, + "particle_position_z": 21237901.246521123, "particle_weight": 7.168263344048695e+28 } } diff --git a/Regression/Checksum/benchmarks_json/collisionXZ.json b/Regression/Checksum/benchmarks_json/collisionXZ.json index 927848745a8..49138f0a1df 100644 --- a/Regression/Checksum/benchmarks_json/collisionXZ.json +++ b/Regression/Checksum/benchmarks_json/collisionXZ.json @@ -6,5 +6,21 @@ "Ex": 0.0, "Ey": 0.0, "Ez": 0.0 + }, + "ion": { + "particle_momentum_x": 2.458306853810186e-19, + "particle_momentum_y": 2.272685285153902e-19, + "particle_momentum_z": 2.281205462681013e-19, + "particle_position_x": 2645436.647039526, + "particle_position_y": 2672571.48688055, + "particle_weight": 1.7256099431746894e+26 + }, + "electron": { + "particle_momentum_x": 1.0454942263455085e-19, + "particle_momentum_y": 1.0323735347957779e-19, + "particle_momentum_z": 1.0199134968670343e-19, + "particle_position_x": 2681776.3648108337, + "particle_position_y": 2663907.8843079703, + "particle_weight": 1.7256099431746894e+26 } } From 49df0ca46930229678c3bebe063f949e24b63013 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Thu, 25 Jul 2024 19:15:13 +0200 Subject: [PATCH 03/12] [WIP] Recompute the macroscopic properties everytime the moving window moves. (#5082) * Recompute the macroscopic properties everytime the moving window moves. * Minor cleanup * Separate allocation and initialization --- .../MacroscopicProperties.H | 14 +++++++++---- .../MacroscopicProperties.cpp | 20 ++++++++++--------- Source/Initialization/WarpXInitData.cpp | 5 +---- Source/Utils/WarpXMovingWindow.cpp | 12 +++++++++++ Source/WarpX.cpp | 7 +++++++ 5 files changed, 41 insertions(+), 17 deletions(-) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H index 129b2e2f782..ca9c77334ea 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H @@ -38,20 +38,26 @@ public: void ReadParameters (); /** - * \brief Initialize multifabs storing macroscopic multifabs + * \brief Allocate multifabs storing macroscopic multifabs * * @param[in] ba the box array associated to the multifabs E and B * @param[in] dmap the distribution mapping * @param[in] ng_EB_alloc guard cells allocated for multifabs E and B + */ + void AllocateLevelMFs ( + const amrex::BoxArray& ba, + const amrex::DistributionMapping& dm, + const amrex::IntVect& ng_EB_alloc ); + + /** + * \brief Initialize multifabs storing macroscopic multifabs + * * @param[in] geom the geometry * @param[in] Ex_stag staggering of the Ex field * @param[in] Ey_stag staggering of the Ey field * @param[in] Ez_stag staggering of the Ez field */ void InitData ( - const amrex::BoxArray& ba, - const amrex::DistributionMapping& dmap, - const amrex::IntVect& ng_EB_alloc, const amrex::Geometry& geom, const amrex::IntVect& Ex_stag, const amrex::IntVect& Ey_stag, diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp index 1467c7a8c0c..a6a389fe056 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp @@ -119,24 +119,26 @@ MacroscopicProperties::ReadParameters () } void -MacroscopicProperties::InitData ( +MacroscopicProperties::AllocateLevelMFs ( const amrex::BoxArray& ba, const amrex::DistributionMapping& dmap, - const amrex::IntVect& ng_EB_alloc, - const amrex::Geometry& geom, - const amrex::IntVect& Ex_stag, - const amrex::IntVect& Ey_stag, - const amrex::IntVect& Ez_stag) + const amrex::IntVect& ng_EB_alloc ) { - amrex::Print() << Utils::TextMsg::Info("we are in init data of macro"); - - // Define material property multifabs using ba and dmap from WarpX instance // sigma is cell-centered MultiFab m_sigma_mf = std::make_unique(ba, dmap, 1, ng_EB_alloc); // epsilon is cell-centered MultiFab m_eps_mf = std::make_unique(ba, dmap, 1, ng_EB_alloc); // mu is cell-centered MultiFab m_mu_mf = std::make_unique(ba, dmap, 1, ng_EB_alloc); +} + +void +MacroscopicProperties::InitData ( + const amrex::Geometry& geom, + const amrex::IntVect& Ex_stag, + const amrex::IntVect& Ey_stag, + const amrex::IntVect& Ez_stag) +{ // Initialize sigma if (m_sigma_s == "constant") { diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 8af0ce23ce8..e9f63e34bb0 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -474,12 +474,9 @@ WarpX::InitData () BuildBufferMasks(); - if (WarpX::em_solver_medium==1) { + if (WarpX::em_solver_medium == MediumForEM::Macroscopic) { const int lev_zero = 0; m_macroscopic_properties->InitData( - boxArray(lev_zero), - DistributionMap(lev_zero), - getngEB(), Geom(lev_zero), getField(warpx::fields::FieldType::Efield_fp, lev_zero,0).ixType().toIntVect(), getField(warpx::fields::FieldType::Efield_fp, lev_zero,1).ixType().toIntVect(), diff --git a/Source/Utils/WarpXMovingWindow.cpp b/Source/Utils/WarpXMovingWindow.cpp index 5fcaccf761a..d64394dafc9 100644 --- a/Source/Utils/WarpXMovingWindow.cpp +++ b/Source/Utils/WarpXMovingWindow.cpp @@ -19,6 +19,7 @@ #include "Utils/TextMsg.H" #include "Utils/WarpXConst.H" #include "Utils/WarpXProfilerWrapper.H" +#include "FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H" #include @@ -452,6 +453,17 @@ WarpX::MoveWindow (const int step, bool move_j) } } + // Recompute macroscopic properties of the medium + if (WarpX::em_solver_medium == MediumForEM::Macroscopic) { + const int lev_zero = 0; + m_macroscopic_properties->InitData( + Geom(lev_zero), + getField(warpx::fields::FieldType::Efield_fp, lev_zero,0).ixType().toIntVect(), + getField(warpx::fields::FieldType::Efield_fp, lev_zero,1).ixType().toIntVect(), + getField(warpx::fields::FieldType::Efield_fp, lev_zero,2).ixType().toIntVect() + ); + } + return num_shift_base; } diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index fcd10418b23..b725431c2f1 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -2373,6 +2373,13 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm myfl->InitData(lev, geom[lev].Domain(),cur_time); } + // Allocate extra multifabs for macroscopic properties of the medium + if (em_solver_medium == MediumForEM::Macroscopic) { + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( lev==0, + "Macroscopic properties are not supported with mesh refinement."); + m_macroscopic_properties->AllocateLevelMFs(ba, dm, ngEB); + } + if (fft_do_time_averaging) { AllocInitMultiFab(Bfield_avg_fp[lev][0], amrex::convert(ba, Bx_nodal_flag), dm, ncomps, ngEB, lev, "Bfield_avg_fp[x]", 0.0_rt); From eda6af5525d89eb251305f340f2f9b6a955b1b48 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 25 Jul 2024 19:26:15 +0200 Subject: [PATCH 04/12] Update script to compile WarpX on Frontier (OLCF) (#5076) * update profile script for Frontier * satisfy dependency for module load * change module versions and manually install adios2 * fix bug * fix bug * update instructions * update scripts * remove line for LibEnsemble --- .../frontier-olcf/frontier_warpx.profile.example | 14 +++++++------- .../machines/frontier-olcf/install_dependencies.sh | 3 --- 2 files changed, 7 insertions(+), 10 deletions(-) diff --git a/Tools/machines/frontier-olcf/frontier_warpx.profile.example b/Tools/machines/frontier-olcf/frontier_warpx.profile.example index 50e12e99a38..f59f2d3d058 100644 --- a/Tools/machines/frontier-olcf/frontier_warpx.profile.example +++ b/Tools/machines/frontier-olcf/frontier_warpx.profile.example @@ -8,9 +8,9 @@ if [ -z ${proj-} ]; then echo "WARNING: The 'proj' variable is not yet set in yo # required dependencies module load cmake/3.23.2 module load craype-accel-amd-gfx90a -module load rocm/5.2.0 # waiting for 5.6 for next bump -module load cray-mpich -module load cce/15.0.0 # must be loaded after rocm +module load rocm/5.7.1 +module load cray-mpich/8.1.28 +module load cce/17.0.0 # must be loaded after rocm # optional: faster builds module load ccache @@ -26,14 +26,14 @@ export LD_LIBRARY_PATH=${HOME}/sw/frontier/gpu/blaspp-2024.05.31/lib64:$LD_LIBRA export LD_LIBRARY_PATH=${HOME}/sw/frontier/gpu/lapackpp-2024.05.31/lib64:$LD_LIBRARY_PATH # optional: for QED lookup table generation support -module load boost/1.79.0-cxx17 +module load boost/1.79.0 # optional: for openPMD support -module load adios2/2.8.3 -module load hdf5/1.14.0 +module load adios2/2.8.3-mpi +module load hdf5/1.12.1-mpi # optional: for Python bindings or libEnsemble -module load cray-python/3.9.13.1 +module load cray-python/3.11.5 if [ -d "${HOME}/sw/frontier/gpu/venvs/warpx-frontier" ] then diff --git a/Tools/machines/frontier-olcf/install_dependencies.sh b/Tools/machines/frontier-olcf/install_dependencies.sh index 503ab4525e9..fd1d28e76b5 100755 --- a/Tools/machines/frontier-olcf/install_dependencies.sh +++ b/Tools/machines/frontier-olcf/install_dependencies.sh @@ -74,7 +74,6 @@ CXX=$(which CC) CXXFLAGS="-DLAPACK_FORTRAN_ADD_" cmake -S $HOME/src/lapackpp -B cmake --build $HOME/src/lapackpp-frontier-gpu-build --target install --parallel 16 rm -rf $HOME/src/lapackpp-frontier-gpu-build - # Python ###################################################################### # python3 -m pip install --upgrade pip @@ -109,8 +108,6 @@ CUPY_INSTALL_USE_HIP=1 \ ROCM_HOME=${ROCM_PATH} \ HCC_AMDGPU_TARGET=${AMREX_AMD_ARCH} \ python3 -m pip install -v cupy -# optional: for libEnsemble -python3 -m pip install -r $HOME/src/warpx/Tools/LibEnsemble/requirements.txt # optional: for optimas (based on libEnsemble & ax->botorch->gpytorch->pytorch) #python3 -m pip install --upgrade torch --index-url https://download.pytorch.org/whl/rocm5.4.2 #python3 -m pip install -r $HOME/src/warpx/Tools/optimas/requirements.txt From 675369d9797aa2bcfda8b64c8e500fff6a2b9f04 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Thu, 25 Jul 2024 10:29:56 -0700 Subject: [PATCH 05/12] PICMI (Bucket): NumPy 2.0 Compatibility (#5075) * PICMI (Bucket): NumPy 2.0 Compatibility `repr` returns `'np.float64(value)'` instead of `'value'` in NumPy 2.0 * Use fstrings * Remove outdated comment * `requirements.txt`: Relax `numpy` version Version 1 and 2 support. * Remove outdated comment --------- Co-authored-by: David Grote --- Python/pywarpx/Bucket.py | 6 ++---- requirements.txt | 2 +- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/Python/pywarpx/Bucket.py b/Python/pywarpx/Bucket.py index 91a34b9d3d6..9dbe4def88e 100644 --- a/Python/pywarpx/Bucket.py +++ b/Python/pywarpx/Bucket.py @@ -59,8 +59,6 @@ def attrlist(self): for attr, value in self.argvattrs.items(): if value is None: continue - # --- repr is applied to value so that for floats, all of the digits are included. - # --- The strip of "'" is then needed when value is a string. if isinstance(value, str): if value.find('=') > -1: # --- Expressions with temporary variables need to be inside quotes @@ -73,11 +71,11 @@ def attrlist(self): continue # --- For lists, tuples, and arrays make a space delimited string of the values. # --- The lambda is needed in case this is a list of strings. - rhs = ' '.join(map(lambda s : repr(s).strip("'"), value)) + rhs = ' '.join(map(lambda s : f'{s}', value)) elif isinstance(value, bool): rhs = 1 if value else 0 else: rhs = value - attrstring = '{0}.{1} = {2}'.format(self.instancename, attr, repr(rhs).strip("'")) + attrstring = f'{self.instancename}.{attr} = {rhs}' result += [attrstring] return result diff --git a/requirements.txt b/requirements.txt index bb0d8ebb704..8e664ae3096 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,5 +1,5 @@ # basic dependencies -numpy~=1.15 +numpy>=1.15 periodictable~=1.5 # PICMI From 5343aa87aa1d98fb68aafd6546db09b0e69f8d8b Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Mon, 29 Jul 2024 11:24:32 +0200 Subject: [PATCH 06/12] Add PRL citation (#5093) --- Docs/source/highlights.rst | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/Docs/source/highlights.rst b/Docs/source/highlights.rst index 7baec74d606..108f685a551 100644 --- a/Docs/source/highlights.rst +++ b/Docs/source/highlights.rst @@ -14,6 +14,11 @@ Plasma-Based Acceleration Scientific works in laser-plasma and beam-plasma acceleration. +#. Shrock JE, Rockafellow E, Miao B, Le M, Hollinger RC, Wang S, Gonsalves AJ, Picksley A, Rocca JJ, and Milchberg HM + **Guided Mode Evolution and Ionization Injection in Meter-Scale Multi-GeV Laser Wakefield Accelerators**. + Phys. Rev. Lett. **133**, 045002, 2024 + `DOI:10.1103/PhysRevLett.133.045002 `__ + #. Ross AJ, Chappell J, van de Wetering JJ, Cowley J, Archer E, Bourgeois N, Corner L, Emerson DR, Feder L, Gu XJ, Jakobsson O, Jones H, Picksley A, Reid L, Wang W, Walczak R, Hooker SM. **Resonant excitation of plasma waves in a plasma channel**. Phys. Rev. Research **6**, L022001, 2024 From b1e7932cc2df0f07ede608f2a2eaf9cbde99d600 Mon Sep 17 00:00:00 2001 From: Justin Ray Angus Date: Mon, 29 Jul 2024 12:05:13 -0700 Subject: [PATCH 07/12] Improving Coulomb collision method for weighted-particles (#5091) * using corrected weighted-particle Coulomb collision method. * adding CI test for weighted Coulomb collisions. * using n12 in UpdateMomentumPerezElastic * updating Checksum * fixing type issue. * updating checksums * checksum --- .../Tests/collision/analysis_collision_1d.py | 126 ++++++++++++++++++ Examples/Tests/collision/inputs_1d | 90 +++++++++++++ .../benchmarks_json/collisionISO.json | 12 +- .../benchmarks_json/collisionXYZ.json | 28 ++-- .../Checksum/benchmarks_json/collisionXZ.json | 20 +-- .../Checksum/benchmarks_json/collisionZ.json | 20 +++ Regression/WarpX-tests.ini | 14 ++ .../BinaryCollision/BinaryCollision.H | 61 +-------- .../Coulomb/ElasticCollisionPerez.H | 31 +++-- .../Coulomb/PairWiseCoulombCollisionFunc.H | 7 +- .../Coulomb/UpdateMomentumPerezElastic.H | 11 +- .../Collision/BinaryCollision/DSMC/DSMCFunc.H | 1 - .../NuclearFusion/NuclearFusionFunc.H | 1 - 13 files changed, 313 insertions(+), 109 deletions(-) create mode 100755 Examples/Tests/collision/analysis_collision_1d.py create mode 100644 Examples/Tests/collision/inputs_1d create mode 100644 Regression/Checksum/benchmarks_json/collisionZ.json diff --git a/Examples/Tests/collision/analysis_collision_1d.py b/Examples/Tests/collision/analysis_collision_1d.py new file mode 100755 index 00000000000..7775a476dae --- /dev/null +++ b/Examples/Tests/collision/analysis_collision_1d.py @@ -0,0 +1,126 @@ +#!/usr/bin/env python3 + +# Copyright 2024 Justin Angus +# +# +# This file is part of WarpX. +# +# License: BSD-3-Clause-LBNL +# +# This is a script that analyses the simulation results from the script `inputs_1d`. +# run locally: python analysis_vandb_1d.py diags/diag1000600/ +# +# This is a 1D intra-species Coulomb scattering relaxation test consisting +# of a low-density population streaming into a higher density population at rest. +# Both populations belong to the same carbon12 ion species. +# See test T1b from JCP 413 (2020) by D. Higginson, et al. +# +import os +import sys + +import numpy as np +import yt +from scipy.constants import e + +sys.path.insert(1, '../../../../warpx/Regression/Checksum/') +import checksumAPI + +# this will be the name of the plot file +fn = sys.argv[1] +ds = yt.load(fn) +data = ds.covering_grid(level = 0, left_edge = ds.domain_left_edge, dims = ds.domain_dimensions) + +# carbon 12 ion (mass = 12*amu - 6*me) +mass = 1.992100316897910e-26 + +# Separate macroparticles from group A (low weight) and group B (high weight) +# by sorting based on weight +sorted_indices = data['ions','particle_weight'].argsort() +sorted_wp = data['ions', 'particle_weight'][sorted_indices].value +sorted_px = data['ions', 'particle_momentum_x'][sorted_indices].value +sorted_py = data['ions', 'particle_momentum_y'][sorted_indices].value +sorted_pz = data['ions', 'particle_momentum_z'][sorted_indices].value + +# Find the index 'Npmin' that separates macroparticles from group A and group B +Np = len(sorted_wp) +wpmin = sorted_wp.min(); +wpmax = sorted_wp.max(); +for i in range(len(sorted_wp)): + if sorted_wp[i] > wpmin: + Npmin = i + break + +NpA = Npmin +wpA = wpmin; +NpB = Np - Npmin +wpB = wpmax; +NpAs = 0 +NpAe = Npmin +NpBs = Npmin +NpBe = Np + +############# + +sorted_px_sum = np.abs(sorted_px).sum(); +sorted_py_sum = np.abs(sorted_py).sum(); +sorted_pz_sum = np.abs(sorted_pz).sum(); +sorted_wp_sum = np.abs(sorted_wp).sum(); + +# compute mean velocities +wAtot = wpA*NpA +wBtot = wpB*NpB + +uBx = uBy = uBz = 0. +for i in range(NpBs,NpBe): + uBx += wpB*sorted_px[i] + uBy += wpB*sorted_py[i] + uBz += wpB*sorted_pz[i] +uBx /= (mass*wBtot) # [m/s] +uBy /= (mass*wBtot) # [m/s] +uBz /= (mass*wBtot) # [m/s] + +uAx = uAy = uAz = 0. +for i in range(NpAs,NpAe): + uAx += wpA*sorted_px[i] + uAy += wpA*sorted_py[i] + uAz += wpA*sorted_pz[i] +uAx /= (mass*wAtot) # [m/s] +uAy /= (mass*wAtot) # [m/s] +uAz /= (mass*wAtot) # [m/s] + +# compute temperatures +TBx = TBy = TBz = 0. +for i in range(NpBs,NpBe): + TBx += wpB*(sorted_px[i]/mass - uBx)**2 + TBy += wpB*(sorted_py[i]/mass - uBy)**2 + TBz += wpB*(sorted_pz[i]/mass - uBz)**2 +TBx *= mass/(e*wBtot) +TBy *= mass/(e*wBtot) +TBz *= mass/(e*wBtot) + +TAx = TAy = TAz = 0. +for i in range(NpAs,NpAe): + TAx += wpA*(sorted_px[i]/mass - uAx)**2 + TAy += wpA*(sorted_py[i]/mass - uAy)**2 + TAz += wpA*(sorted_pz[i]/mass - uAz)**2 +TAx *= mass/(e*wAtot) +TAy *= mass/(e*wAtot) +TAz *= mass/(e*wAtot) + +TApar = TAz +TAperp = (TAx + TAy)/2.0 +TA = (TAx + TAy + TAz)/3.0 + +TBpar = TBz +TBperp = (TBx + TBy)/2.0 +TB = (TBx + TBy + TBz)/3.0 + +TApar_30ps_soln = 6.15e3 # TA parallel solution at t = 30 ps +error = np.abs(TApar-TApar_30ps_soln)/TApar_30ps_soln +tolerance = 0.02 +print('TApar at 30ps error = ', error); +print('tolerance = ', tolerance); +assert error < tolerance + +test_name = os.path.split(os.getcwd())[1] +checksumAPI.evaluate_checksum(test_name, fn) diff --git a/Examples/Tests/collision/inputs_1d b/Examples/Tests/collision/inputs_1d new file mode 100644 index 00000000000..b2de17192ae --- /dev/null +++ b/Examples/Tests/collision/inputs_1d @@ -0,0 +1,90 @@ +################################# +########## CONSTANTS ############ +################################# + +my_constants.nA = 1.e25 # m^-3 +my_constants.NpA = 400 # m^-3 +my_constants.UA = 6.55e5 # m/s +my_constants.TA = 500. # eV +# +my_constants.nB = 1.e26 # m^-3 +my_constants.NpB = 400 # m^-3 +my_constants.UB = 0. # m/s +my_constants.TB = 500. # eV +# +my_constants.q_c12 = 6.*q_e +my_constants.m_c12 = 12.*m_u - 6.*m_e + +################################# +####### GENERAL PARAMETERS ###### +################################# +max_step = 600 +amr.n_cell = 180 +amr.max_level = 0 +amr.blocking_factor = 4 +geometry.dims = 1 +geometry.prob_lo = 0. +geometry.prob_hi = 0.01 + +################################# +###### Boundary Condition ####### +################################# +boundary.field_lo = periodic +boundary.field_hi = periodic + +################################# +############ NUMERICS ########### +################################# +warpx.serialize_initial_conditions = 1 +warpx.verbose = 1 +warpx.const_dt = 0.05e-12 +warpx.use_filter = 0 + +# Do not evolve the E and B fields +algo.maxwell_solver = none + +# Order of particle shape factors +algo.particle_shape = 1 + +################################# +############ PLASMA ############# +################################# +particles.species_names = ions + +ions.charge = q_c12 +ions.mass = m_c12 +ions.do_not_deposit = 1 + +ions.injection_sources = groupA groupB + +ions.groupA.injection_style = "NUniformPerCell" +ions.groupA.num_particles_per_cell_each_dim = NpA +ions.groupA.profile = constant +ions.groupA.density = nA # number per m^3 +ions.groupA.momentum_distribution_type = "gaussian" +ions.groupA.uz_m = UA/clight +ions.groupA.ux_th = sqrt(TA*q_e/m_c12)/clight +ions.groupA.uy_th = sqrt(TA*q_e/m_c12)/clight +ions.groupA.uz_th = sqrt(TA*q_e/m_c12)/clight + +ions.groupB.injection_style = "NUniformPerCell" +ions.groupB.num_particles_per_cell_each_dim = NpB +ions.groupB.profile = constant +ions.groupB.density = nB # number per m^3 +ions.groupB.momentum_distribution_type = "gaussian" +ions.groupB.uz_m = UB/clight +ions.groupB.ux_th = sqrt(TB*q_e/m_c12)/clight +ions.groupB.uy_th = sqrt(TB*q_e/m_c12)/clight +ions.groupB.uz_th = sqrt(TB*q_e/m_c12)/clight + +################################# +############ COLLISION ########## +################################# +collisions.collision_names = collision1 +collision1.species = ions ions +collision1.CoulombLog = 10.0 + +# Diagnostics +diagnostics.diags_names = diag1 +diag1.intervals = 600 +diag1.diag_type = Full diff --git a/Regression/Checksum/benchmarks_json/collisionISO.json b/Regression/Checksum/benchmarks_json/collisionISO.json index a2fd6116cb8..350848d4aee 100644 --- a/Regression/Checksum/benchmarks_json/collisionISO.json +++ b/Regression/Checksum/benchmarks_json/collisionISO.json @@ -11,12 +11,12 @@ "jz": 0.0 }, "electron": { - "particle_momentum_x": 3.579989064013309e-19, - "particle_momentum_y": 3.5822945977746767e-19, - "particle_momentum_z": 3.579753452653627e-19, - "particle_position_x": 1.0241322532163375, - "particle_position_y": 1.0238995904625479, - "particle_position_z": 1.02402135051502, + "particle_momentum_x": 3.5790777034053853e-19, + "particle_momentum_y": 3.5815348106229496e-19, + "particle_momentum_z": 3.577963316718249e-19, + "particle_position_x": 1.024180253191667, + "particle_position_y": 1.023919590453571, + "particle_position_z": 1.0240653505082926, "particle_weight": 714240000000.0 } } diff --git a/Regression/Checksum/benchmarks_json/collisionXYZ.json b/Regression/Checksum/benchmarks_json/collisionXYZ.json index 2aa0e4686eb..c13f404857a 100644 --- a/Regression/Checksum/benchmarks_json/collisionXYZ.json +++ b/Regression/Checksum/benchmarks_json/collisionXYZ.json @@ -6,25 +6,25 @@ "Ex": 0.0, "Ey": 0.0, "Ez": 0.0, - "T_electron": 381488.6528070591, - "T_ion": 320091.2785835478 + "T_electron": 381957.7394223898, + "T_ion": 319565.6269784763 }, "electron": { - "particle_momentum_x": 8.667025573698235e-19, - "particle_momentum_y": 8.457499789250831e-19, - "particle_momentum_z": 8.482438182280524e-19, - "particle_position_x": 21262567.138872623, - "particle_position_y": 21245135.070665065, - "particle_position_z": 21232644.283726066, + "particle_momentum_x": 8.631833485301232e-19, + "particle_momentum_y": 8.476316745254719e-19, + "particle_momentum_z": 8.514139891331418e-19, + "particle_position_x": 21253105.73686561, + "particle_position_y": 21282643.519070115, + "particle_position_z": 21239057.457948968, "particle_weight": 7.168263344048695e+28 }, "ion": { - "particle_momentum_x": 1.9300495097720012e-18, - "particle_momentum_y": 1.747257416857836e-18, - "particle_momentum_z": 1.7510296287537058e-18, - "particle_position_x": 21217348.883301035, - "particle_position_y": 21300859.0630925, - "particle_position_z": 21237901.246521123, + "particle_momentum_x": 1.9215585867122464e-18, + "particle_momentum_y": 1.7471481568315848e-18, + "particle_momentum_z": 1.7510887207292533e-18, + "particle_position_x": 21228900.948879313, + "particle_position_y": 21304564.011731848, + "particle_position_z": 21250585.221808463, "particle_weight": 7.168263344048695e+28 } } diff --git a/Regression/Checksum/benchmarks_json/collisionXZ.json b/Regression/Checksum/benchmarks_json/collisionXZ.json index 49138f0a1df..f90c34bc86d 100644 --- a/Regression/Checksum/benchmarks_json/collisionXZ.json +++ b/Regression/Checksum/benchmarks_json/collisionXZ.json @@ -8,19 +8,19 @@ "Ez": 0.0 }, "ion": { - "particle_momentum_x": 2.458306853810186e-19, - "particle_momentum_y": 2.272685285153902e-19, - "particle_momentum_z": 2.281205462681013e-19, - "particle_position_x": 2645436.647039526, - "particle_position_y": 2672571.48688055, + "particle_momentum_x": 2.4932317055825563e-19, + "particle_momentum_y": 2.274916403334278e-19, + "particle_momentum_z": 2.2528161767665816e-19, + "particle_position_x": 2648815.601036139, + "particle_position_y": 2662836.7581390506, "particle_weight": 1.7256099431746894e+26 }, "electron": { - "particle_momentum_x": 1.0454942263455085e-19, - "particle_momentum_y": 1.0323735347957779e-19, - "particle_momentum_z": 1.0199134968670343e-19, - "particle_position_x": 2681776.3648108337, - "particle_position_y": 2663907.8843079703, + "particle_momentum_x": 1.0489203687862582e-19, + "particle_momentum_y": 1.0209657029567292e-19, + "particle_momentum_z": 1.0248962872393911e-19, + "particle_position_x": 2657004.8285825616, + "particle_position_y": 2670174.272797987, "particle_weight": 1.7256099431746894e+26 } } diff --git a/Regression/Checksum/benchmarks_json/collisionZ.json b/Regression/Checksum/benchmarks_json/collisionZ.json new file mode 100644 index 00000000000..3be8d5ae893 --- /dev/null +++ b/Regression/Checksum/benchmarks_json/collisionZ.json @@ -0,0 +1,20 @@ +{ + "lev=0": { + "Bx": 0.0, + "By": 0.0, + "Bz": 0.0, + "Ex": 0.0, + "Ey": 0.0, + "Ez": 0.0, + "jx": 0.0, + "jy": 0.0, + "jz": 0.0 + }, + "ions": { + "particle_momentum_x": 3.425400072687143e-16, + "particle_momentum_y": 3.421937133999805e-16, + "particle_momentum_z": 5.522701882677923e-16, + "particle_position_x": 720.0011611411148, + "particle_weight": 1.0999999999999999e+24 + } +} diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 17ecfd64f29..c1a6316e7e2 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -195,6 +195,20 @@ useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/collider_relevant_diags/analysis_multiple_particles.py +[collisionZ] +buildDir = . +inputFile = Examples/Tests/collision/inputs_1d +runtime_params = +dim = 1 +addToCompileString = +cmakeSetupOpts = -DWarpX_DIMS=1 +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 1 +analysisRoutine = Examples/Tests/collision/analysis_collision_1d.py + [collisionISO] buildDir = . inputFile = Examples/Tests/collision/inputs_3d_isotropization diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H index 0ba49e9b09f..6f8f61e66b1 100644 --- a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H +++ b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H @@ -363,17 +363,14 @@ public: // create vectors to store density and temperature on cell level amrex::Gpu::DeviceVector n1_vec; - amrex::Gpu::DeviceVector n12_vec; amrex::Gpu::DeviceVector T1_vec; if (binary_collision_functor.m_computeSpeciesDensities) { n1_vec.resize(n_cells); - n12_vec.resize(n_cells); } if (binary_collision_functor.m_computeSpeciesTemperatures) { T1_vec.resize(n_cells); } amrex::ParticleReal* AMREX_RESTRICT n1_in_each_cell = n1_vec.dataPtr(); - amrex::ParticleReal* AMREX_RESTRICT n12_in_each_cell = n12_vec.dataPtr(); amrex::ParticleReal* AMREX_RESTRICT T1_in_each_cell = T1_vec.dataPtr(); // Loop over cells @@ -384,7 +381,6 @@ public: // given by the `indices_1[cell_start_1:cell_stop_1]` index_type const cell_start_1 = cell_offsets_1[i_cell]; index_type const cell_stop_1 = cell_offsets_1[i_cell+1]; - index_type const cell_half_1 = (cell_start_1+cell_stop_1)/2; // Do not collide if there is only one particle in the cell if ( cell_stop_1 - cell_start_1 <= 1 ) { return; } @@ -415,29 +411,6 @@ public: // shuffle ShuffleFisherYates(indices_1, cell_start_1, cell_stop_1, engine); - - // compute n12 for intra-species - if (binary_collision_functor.m_computeSpeciesDensities) { - amrex::ParticleReal n12 = 0.0; - index_type const NI1 = cell_half_1 - cell_start_1; - index_type const NI2 = cell_stop_1 - cell_half_1; - index_type const max_N = amrex::max(NI1,NI2); - index_type i1 = cell_start_1; - index_type i2 = cell_half_1; - amrex::ParticleReal * const AMREX_RESTRICT w1 = soa_1.m_rdata[PIdx::w]; - for (index_type k=0; k n1_vec, n2_vec; - amrex::Gpu::DeviceVector n12_vec; amrex::Gpu::DeviceVector T1_vec, T2_vec; if (binary_collision_functor.m_computeSpeciesDensities) { n1_vec.resize(n_cells); n2_vec.resize(n_cells); - n12_vec.resize(n_cells); } if (binary_collision_functor.m_computeSpeciesTemperatures) { T1_vec.resize(n_cells); @@ -655,7 +624,6 @@ public: } amrex::ParticleReal* AMREX_RESTRICT n1_in_each_cell = n1_vec.dataPtr(); amrex::ParticleReal* AMREX_RESTRICT n2_in_each_cell = n2_vec.dataPtr(); - amrex::ParticleReal* AMREX_RESTRICT n12_in_each_cell = n12_vec.dataPtr(); amrex::ParticleReal* AMREX_RESTRICT T1_in_each_cell = T1_vec.dataPtr(); amrex::ParticleReal* AMREX_RESTRICT T2_in_each_cell = T2_vec.dataPtr(); @@ -719,29 +687,6 @@ public: // shuffle ShuffleFisherYates(indices_1, cell_start_1, cell_stop_1, engine); ShuffleFisherYates(indices_2, cell_start_2, cell_stop_2, engine); - - // compute n12 for inter-species - if (binary_collision_functor.m_computeSpeciesDensities) { - amrex::ParticleReal n12 = 0.0; - index_type const NI1 = cell_stop_1 - cell_start_1; - index_type const NI2 = cell_stop_2 - cell_start_2; - index_type const max_N = amrex::max(NI1,NI2); - index_type i1 = cell_start_1; - index_type i2 = cell_start_2; - amrex::ParticleReal * const AMREX_RESTRICT w1 = soa_1.m_rdata[PIdx::w]; - amrex::ParticleReal * const AMREX_RESTRICT w2 = soa_2.m_rdata[PIdx::w]; - for (index_type k=0; k( 1.0/std::cbrt(4.0*MathConst::pi/3.0*maxn) ); + + // bmax (screening length) cannot be smaller than atomic spacing + const T_PR bmax = amrex::max(lmdD, rmin); #if (defined WARPX_DIM_RZ) T_PR * const AMREX_RESTRICT theta1 = soa_1.m_rdata[PIdx::theta]; @@ -108,12 +110,23 @@ void ElasticCollisionPerez ( u1y[I1[i1]] = u1xbuf*std::sin(theta) + u1y[I1[i1]]*std::cos(theta); #endif + // Compute the effective density n12 used to compute the normalized + // scattering path s12 in UpdateMomentumPerezElastic(). + // s12 is defined such that the expected value of the change in particle + // velocity is equal to that from the full NxN pairing method, as described + // here https://arxiv.org/submit/5758216/view. This method is a direct extension + // of the original method by Takizuka and Abe JCP 25 (1977) to weighted particles. + T_PR n12; + const T_PR wpmax = amrex::max(w1[ I1[i1] ],w2[ I2[i2] ]); + if (isSameSpecies) { n12 = wpmax*static_cast(min_N+max_N-1)/dV; } + else { n12 = wpmax*static_cast(min_N)/dV; } + UpdateMomentumPerezElastic( u1x[ I1[i1] ], u1y[ I1[i1] ], u1z[ I1[i1] ], u2x[ I2[i2] ], u2y[ I2[i2] ], u2z[ I2[i2] ], n1, n2, n12, q1, m1, w1[ I1[i1] ], q2, m2, w2[ I2[i2] ], - dt, L, lmdD, + dt, L, bmax, engine); #if (defined WARPX_DIM_RZ) diff --git a/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H b/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H index cd3cd5a5c83..3322c19ed2d 100644 --- a/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H @@ -91,11 +91,10 @@ public: const SoaData_type& soa_1, const SoaData_type& soa_2, GetParticlePosition /*get_position_1*/, GetParticlePosition /*get_position_2*/, amrex::ParticleReal const n1, amrex::ParticleReal const n2, - amrex::ParticleReal const n12, amrex::ParticleReal const T1, amrex::ParticleReal const T2, amrex::ParticleReal const q1, amrex::ParticleReal const q2, amrex::ParticleReal const m1, amrex::ParticleReal const m2, - amrex::Real const dt, amrex::Real const /*dV*/, index_type coll_idx, + amrex::Real const dt, amrex::Real const dV, index_type coll_idx, index_type const /*cell_start_pair*/, index_type* /*p_mask*/, index_type* /*p_pair_indices_1*/, index_type* /*p_pair_indices_2*/, amrex::ParticleReal* /*p_pair_reaction_weight*/, @@ -105,9 +104,9 @@ public: ElasticCollisionPerez( I1s, I1e, I2s, I2e, I1, I2, - soa_1, soa_2, n1, n2, n12, T1, T2, + soa_1, soa_2, n1, n2, T1, T2, q1, q2, m1, m2, - dt, m_CoulombLog, engine, coll_idx); + dt, m_CoulombLog, dV, engine, m_isSameSpecies, coll_idx); } amrex::ParticleReal m_CoulombLog; diff --git a/Source/Particles/Collision/BinaryCollision/Coulomb/UpdateMomentumPerezElastic.H b/Source/Particles/Collision/BinaryCollision/Coulomb/UpdateMomentumPerezElastic.H index 3101047e211..93acc3d8aef 100644 --- a/Source/Particles/Collision/BinaryCollision/Coulomb/UpdateMomentumPerezElastic.H +++ b/Source/Particles/Collision/BinaryCollision/Coulomb/UpdateMomentumPerezElastic.H @@ -18,9 +18,10 @@ /* \brief Update particle velocities according to * F. Perez et al., Phys.Plasmas.19.083104 (2012), * which is based on Nanbu's method, PhysRevE.55.4642 (1997). - * @param[in] LmdD is max(Debye length, minimal interparticle distance). + * @param[in] bmax is max(Debye length, minimal interparticle distance). * @param[in] L is the Coulomb log. A fixed L will be used if L > 0, * otherwise L will be calculated based on the algorithm. + * @param[in] n12 = max(w1,w2)*min(N1,N2)/dV is the effective density used for s12 * To see if there are nan or inf updated velocities, * compile with USE_ASSERTION=TRUE. * @@ -36,7 +37,7 @@ void UpdateMomentumPerezElastic ( T_PR const n1, T_PR const n2, T_PR const n12, T_PR const q1, T_PR const m1, T_PR const w1, T_PR const q2, T_PR const m2, T_PR const w2, - T_R const dt, T_PR const L, T_PR const lmdD, + T_R const dt, T_PR const L, T_PR const bmax, amrex::RandomEngine const& engine) { @@ -128,13 +129,13 @@ void UpdateMomentumPerezElastic ( // Compute the Coulomb log lnLmd lnLmd = amrex::max( T_PR(2.0), - T_PR(0.5)*std::log(T_PR(1.0)+lmdD*lmdD/(bmin*bmin)) ); + T_PR(0.5)*std::log(T_PR(1.0) + bmax*bmax/(bmin*bmin)) ); } // Compute s const auto tts = m1*g1s*m2*g2s/(inv_c2*p1sm*p1sm) + T_PR(1.0); const auto tts2 = tts*tts; - s = n1*n2/n12 * dt*lnLmd*q1*q1*q2*q2 / + s = n12 * dt*lnLmd*q1*q1*q2*q2 / ( T_PR(4.0) * MathConst::pi * PhysConst::ep0 * PhysConst::ep0 * m1*g1*m2*g2/(inv_c2*inv_c2) ) * gc*p1sm/mass_g * tts2; @@ -144,7 +145,7 @@ void UpdateMomentumPerezElastic ( const auto coeff = static_cast( std::pow(4.0*MathConst::pi/3.0,1.0/3.0)); T_PR const vrel = mass_g*p1sm/(m1*g1s*m2*g2s*gc); - T_PR const sp = coeff * n1*n2/n12 * dt * vrel * (m1+m2) / + T_PR const sp = coeff * n12 * dt * vrel * (m1+m2) / amrex::max( m1*cbrt_n1*cbrt_n1, m2*cbrt_n2*cbrt_n2); diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H index 152dca4e1a1..30b466e2ec2 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H +++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H @@ -97,7 +97,6 @@ public: const SoaData_type& soa_1, const SoaData_type& soa_2, GetParticlePosition /*get_position_1*/, GetParticlePosition /*get_position_2*/, amrex::ParticleReal const /*n1*/, amrex::ParticleReal const /*n2*/, - amrex::ParticleReal const /*n12*/, amrex::ParticleReal const /*T1*/, amrex::ParticleReal const /*T2*/, amrex::ParticleReal const /*q1*/, amrex::ParticleReal const /*q2*/, amrex::ParticleReal const m1, amrex::ParticleReal const m2, diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H index 1f3e5e56582..3113dc69839 100644 --- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H +++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H @@ -132,7 +132,6 @@ public: const SoaData_type& soa_1, const SoaData_type& soa_2, GetParticlePosition /*get_position_1*/, GetParticlePosition /*get_position_2*/, amrex::ParticleReal const /*n1*/, amrex::ParticleReal const /*n2*/, - amrex::ParticleReal const /*n12*/, amrex::ParticleReal const /*T1*/, amrex::ParticleReal const /*T2*/, amrex::ParticleReal const /*q1*/, amrex::ParticleReal const /*q2*/, amrex::ParticleReal const m1, amrex::ParticleReal const m2, From e26f7ef0d42903a96cb67c3bf1fbe56af88712c1 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Tue, 30 Jul 2024 02:34:10 +0300 Subject: [PATCH 08/12] Fix FFT python interface, take 2 (#5081) --- Python/pywarpx/picmi.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index 42c0313a6e9..d06fe90d2d6 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -1351,6 +1351,9 @@ def init(self, kw): def solver_initialize_inputs(self): + # Open BC means FieldBoundaryType::Open for electrostatic sims, rather than perfectly-matched layer + BC_map['open'] = 'open' + self.grid.grid_initialize_inputs() if self.relativistic: @@ -1371,6 +1374,8 @@ def solver_initialize_inputs(self): pywarpx.boundary.potential_hi_y = self.grid.potential_ymax pywarpx.boundary.potential_hi_z = self.grid.potential_zmax + pywarpx.warpx.poisson_solver = self.method + class GaussianLaser(picmistandard.PICMI_GaussianLaser): def laser_initialize_inputs(self): From 24c0711fc92577b3ed8c59af316c78d76894e10f Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Mon, 29 Jul 2024 18:35:13 -0500 Subject: [PATCH 09/12] AMReX/pyAMReX/PICSAR: Weekly Update (#5095) * AMReX: Weekly Update * pyAMReX: Weekly Update * Increase Clang-Tidy `timeout-minutes` * Cast long long to int when using IParser --------- Co-authored-by: Edoardo Zoni --- .github/workflows/clang_tidy.yml | 1 + .github/workflows/cuda.yml | 2 +- Regression/WarpX-GPU-tests.ini | 2 +- Regression/WarpX-tests.ini | 2 +- Source/ablastr/utils/SignalHandling.cpp | 2 +- cmake/dependencies/AMReX.cmake | 2 +- cmake/dependencies/pyAMReX.cmake | 2 +- run_test.sh | 2 +- 8 files changed, 8 insertions(+), 7 deletions(-) diff --git a/.github/workflows/clang_tidy.yml b/.github/workflows/clang_tidy.yml index 28b46b52530..af8632ed698 100644 --- a/.github/workflows/clang_tidy.yml +++ b/.github/workflows/clang_tidy.yml @@ -13,6 +13,7 @@ jobs: dim: [1, 2, RZ, 3] name: clang-tidy-${{ matrix.dim }}D runs-on: ubuntu-22.04 + timeout-minutes: 120 if: github.event.pull_request.draft == false steps: - uses: actions/checkout@v4 diff --git a/.github/workflows/cuda.yml b/.github/workflows/cuda.yml index c2ec9aaa763..8f5e59ccaa4 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 0c3273f5e591815909180f8ffaf5b793cabbf9bc && cd - + cd ../amrex && git checkout --detach 20e6f2eadf0c297517588ba38973ec7c7084fa31 && cd - make COMP=gcc QED=FALSE USE_MPI=TRUE USE_GPU=TRUE USE_OMP=FALSE USE_FFT=TRUE USE_CCACHE=TRUE -j 4 ccache -s diff --git a/Regression/WarpX-GPU-tests.ini b/Regression/WarpX-GPU-tests.ini index a2447ac3cf5..68d59d850a0 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 = 0c3273f5e591815909180f8ffaf5b793cabbf9bc +branch = 20e6f2eadf0c297517588ba38973ec7c7084fa31 [source] dir = /home/regtester/git/WarpX diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index c1a6316e7e2..8e2f723b8c5 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 = 0c3273f5e591815909180f8ffaf5b793cabbf9bc +branch = 20e6f2eadf0c297517588ba38973ec7c7084fa31 [source] dir = /home/regtester/AMReX_RegTesting/warpx diff --git a/Source/ablastr/utils/SignalHandling.cpp b/Source/ablastr/utils/SignalHandling.cpp index 4095b9207ce..5eeaeec259f 100644 --- a/Source/ablastr/utils/SignalHandling.cpp +++ b/Source/ablastr/utils/SignalHandling.cpp @@ -104,7 +104,7 @@ SignalHandling::parseSignalNameToNumber (const std::string &str) auto spf = signals_parser.compileHost<0>(); - const int sig = spf(); + const auto sig = int(spf()); ABLASTR_ALWAYS_ASSERT_WITH_MESSAGE(sig < NUM_SIGNALS, "Parsed signal value is outside the supported range of [1, 31]"); diff --git a/cmake/dependencies/AMReX.cmake b/cmake/dependencies/AMReX.cmake index 025b4d3b6e0..f1e5b3f62e2 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 "0c3273f5e591815909180f8ffaf5b793cabbf9bc" +set(WarpX_amrex_branch "20e6f2eadf0c297517588ba38973ec7c7084fa31" 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 a4598953432..133747aaeae 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 "ff4643869c63d4ee40a87054b901f61eefcb97a3" +set(WarpX_pyamrex_branch "e007e730d48cb5fdbe1e10462d7d0a14e2bc8f8e" CACHE STRING "Repository branch for WarpX_pyamrex_repo if(WarpX_pyamrex_internal)") diff --git a/run_test.sh b/run_test.sh index 1692f66263c..f56a957c17a 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 0c3273f5e591815909180f8ffaf5b793cabbf9bc && cd - +cd amrex && git checkout --detach 20e6f2eadf0c297517588ba38973ec7c7084fa31 && 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 3413568d0d372d508de3c56ffdebd25b884f8d97 Mon Sep 17 00:00:00 2001 From: Edoardo Zoni <59625522+EZoni@users.noreply.github.com> Date: Mon, 29 Jul 2024 16:36:10 -0700 Subject: [PATCH 10/12] CI: remove unused parameter `analysisOutputImage` (#5097) --- Regression/WarpX-tests.ini | 45 -------------------------------------- 1 file changed, 45 deletions(-) diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 8e2f723b8c5..a1dc0a168f7 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -380,7 +380,6 @@ useOMP = 1 numthreads = 1 runtime_params = geometry.dims=2 analysisRoutine = Examples/Tests/dive_cleaning/analysis.py -analysisOutputImage = Comparison.png [dive_cleaning_3d] buildDir = . @@ -395,7 +394,6 @@ useOMP = 1 numthreads = 1 runtime_params = analysisRoutine = Examples/Tests/dive_cleaning/analysis.py -analysisOutputImage = Comparison.png [ElectrostaticSphere] buildDir = . @@ -946,7 +944,6 @@ numprocs = 2 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_3d.py -analysisOutputImage = langmuir_multi_analysis.png [Langmuir_fluid_1D] buildDir = . @@ -961,7 +958,6 @@ numprocs = 2 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir_fluids/analysis_1d.py -analysisOutputImage = langmuir_fluid_multi_1d_analysis.png [Langmuir_fluid_RZ] buildDir = . @@ -976,7 +972,6 @@ numprocs = 2 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir_fluids/analysis_rz.py -analysisOutputImage = langmuir_fluid_rz_analysis.png [Langmuir_fluid_2D] buildDir = . @@ -991,7 +986,6 @@ numprocs = 2 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir_fluids/analysis_2d.py -analysisOutputImage = langmuir_fluid_multi_2d_analysis.png [Langmuir_fluid_multi] buildDir = . @@ -1006,7 +1000,6 @@ numprocs = 2 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir_fluids/analysis_3d.py -analysisOutputImage = langmuir_fluid_multi_analysis.png [Langmuir_multi_1d] buildDir = . @@ -1021,7 +1014,6 @@ numprocs = 2 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_1d.py -analysisOutputImage = langmuir_multi_1d_analysis.png [Langmuir_multi_2d_MR] buildDir = . @@ -1036,7 +1028,6 @@ numprocs = 2 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_2d.py -analysisOutputImage = Langmuir_multi_2d_MR.png [Langmuir_multi_2d_MR_anisotropic] buildDir = . @@ -1051,7 +1042,6 @@ numprocs = 2 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_2d.py -analysisOutputImage = Langmuir_multi_2d_MR.png [Langmuir_multi_2d_MR_momentum_conserving] buildDir = . @@ -1066,7 +1056,6 @@ numprocs = 2 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_2d.py -analysisOutputImage = Langmuir_multi_2d_MR_momentum_conserving.png [Langmuir_multi_2d_MR_psatd] buildDir = . @@ -1081,7 +1070,6 @@ numprocs = 2 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_2d.py -analysisOutputImage = Langmuir_multi_2d_MR_psatd.png [Langmuir_multi_2d_nodal] buildDir = . @@ -1096,7 +1084,6 @@ numprocs = 2 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_2d.py -analysisOutputImage = langmuir_multi_2d_analysis.png [Langmuir_multi_2d_psatd] buildDir = . @@ -1111,7 +1098,6 @@ numprocs = 2 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_2d.py -analysisOutputImage = langmuir_multi_2d_analysis.png [Langmuir_multi_2d_psatd_current_correction] buildDir = . @@ -1126,7 +1112,6 @@ numprocs = 1 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_2d.py -analysisOutputImage = langmuir_multi_2d_analysis.png [Langmuir_multi_2d_psatd_current_correction_nodal] buildDir = . @@ -1141,7 +1126,6 @@ numprocs = 1 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_2d.py -analysisOutputImage = langmuir_multi_2d_analysis.png [Langmuir_multi_2d_psatd_momentum_conserving] buildDir = . @@ -1156,7 +1140,6 @@ numprocs = 2 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_2d.py -analysisOutputImage = langmuir_multi_2d_analysis.png [Langmuir_multi_2d_psatd_multiJ] buildDir = . @@ -1171,7 +1154,6 @@ numprocs = 2 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_2d.py -analysisOutputImage = Langmuir_multi_2d_psatd_multiJ.png [Langmuir_multi_2d_psatd_multiJ_nodal] buildDir = . @@ -1186,7 +1168,6 @@ numprocs = 2 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_2d.py -analysisOutputImage = Langmuir_multi_2d_psatd_multiJ_nodal.png [Langmuir_multi_2d_psatd_nodal] buildDir = . @@ -1201,7 +1182,6 @@ numprocs = 2 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_2d.py -analysisOutputImage = langmuir_multi_2d_analysis.png [Langmuir_multi_2d_psatd_Vay_deposition] buildDir = . @@ -1216,7 +1196,6 @@ numprocs = 1 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_2d.py -analysisOutputImage = langmuir_multi_2d_analysis.png [Langmuir_multi_2d_psatd_Vay_deposition_particle_shape_4] buildDir = . @@ -1231,7 +1210,6 @@ numprocs = 1 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_2d.py -analysisOutputImage = langmuir_multi_2d_analysis.png [Langmuir_multi_2d_psatd_Vay_deposition_nodal] buildDir = . @@ -1246,7 +1224,6 @@ numprocs = 1 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_2d.py -analysisOutputImage = langmuir_multi_2d_analysis.png [Langmuir_multi_nodal] buildDir = . @@ -1261,7 +1238,6 @@ numprocs = 2 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_3d.py -analysisOutputImage = langmuir_multi_analysis.png [Langmuir_multi_psatd] buildDir = . @@ -1276,7 +1252,6 @@ numprocs = 2 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_3d.py -analysisOutputImage = langmuir_multi_analysis.png [Langmuir_multi_psatd_current_correction] buildDir = . @@ -1291,7 +1266,6 @@ numprocs = 1 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_3d.py -analysisOutputImage = langmuir_multi_analysis.png [Langmuir_multi_psatd_current_correction_nodal] buildDir = . @@ -1306,7 +1280,6 @@ numprocs = 1 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_3d.py -analysisOutputImage = langmuir_multi_analysis.png [Langmuir_multi_psatd_div_cleaning] buildDir = . @@ -1321,7 +1294,6 @@ numprocs = 2 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_3d.py -analysisOutputImage = langmuir_multi_analysis.png [Langmuir_multi_psatd_momentum_conserving] buildDir = . @@ -1336,7 +1308,6 @@ numprocs = 2 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_3d.py -analysisOutputImage = langmuir_multi_analysis.png [Langmuir_multi_psatd_multiJ] buildDir = . @@ -1351,7 +1322,6 @@ numprocs = 2 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_3d.py -analysisOutputImage = Langmuir_multi_psatd_multiJ.png [Langmuir_multi_psatd_multiJ_nodal] buildDir = . @@ -1366,7 +1336,6 @@ numprocs = 2 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_3d.py -analysisOutputImage = Langmuir_multi_psatd_multiJ_nodal.png [Langmuir_multi_psatd_nodal] buildDir = . @@ -1381,7 +1350,6 @@ numprocs = 2 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_3d.py -analysisOutputImage = langmuir_multi_analysis.png [Langmuir_multi_psatd_single_precision] buildDir = . @@ -1396,7 +1364,6 @@ numprocs = 2 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_3d.py -analysisOutputImage = langmuir_multi_analysis.png [Langmuir_multi_psatd_Vay_deposition] buildDir = . @@ -1411,7 +1378,6 @@ numprocs = 1 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_3d.py -analysisOutputImage = langmuir_multi_analysis.png [Langmuir_multi_psatd_Vay_deposition_nodal] buildDir = . @@ -1426,7 +1392,6 @@ numprocs = 1 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_3d.py -analysisOutputImage = langmuir_multi_analysis.png [Langmuir_multi_rz] buildDir = . @@ -1441,7 +1406,6 @@ numprocs = 2 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_rz.py -analysisOutputImage = Langmuir_multi_rz_analysis.png aux1File = Regression/PostProcessingUtils/post_processing_utils.py [Langmuir_multi_rz_psatd] @@ -1457,7 +1421,6 @@ numprocs = 2 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_rz.py -analysisOutputImage = Langmuir_multi_rz_psatd_analysis.png aux1File = Regression/PostProcessingUtils/post_processing_utils.py [Langmuir_multi_rz_psatd_current_correction] @@ -1473,7 +1436,6 @@ numprocs = 1 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_rz.py -analysisOutputImage = Langmuir_multi_rz_psatd_analysis.png aux1File = Regression/PostProcessingUtils/post_processing_utils.py [Langmuir_multi_rz_psatd_multiJ] @@ -1489,7 +1451,6 @@ numprocs = 2 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_rz.py -analysisOutputImage = Langmuir_multi_rz_psatd_multiJ_analysis.png aux1File = Regression/PostProcessingUtils/post_processing_utils.py [Langmuir_multi_single_precision] @@ -1505,7 +1466,6 @@ numprocs = 2 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/langmuir/analysis_3d.py -analysisOutputImage = langmuir_multi_analysis.png [Larmor] buildDir = . @@ -1677,7 +1637,6 @@ numprocs = 2 useOMP = 1 numthreads = 1 analysisRoutine = Examples/Tests/laser_injection/analysis_laser.py -analysisOutputImage = laser_analysis.png [LaserInjection_1d] buildDir = . @@ -3487,7 +3446,6 @@ useOMP = 1 numthreads = 1 runtime_params = analysisRoutine = Examples/Tests/relativistic_space_charge_initialization/analysis.py -analysisOutputImage = Comparison.png [RepellingParticles] buildDir = . @@ -3701,7 +3659,6 @@ useOMP = 1 numthreads = 1 runtime_params = analysisRoutine = Examples/Tests/space_charge_initialization/analysis.py -analysisOutputImage = Comparison.png [space_charge_initialization_2d] buildDir = . @@ -3716,7 +3673,6 @@ useOMP = 1 numthreads = 1 runtime_params = geometry.dims=2 analysisRoutine = Examples/Tests/space_charge_initialization/analysis.py -analysisOutputImage = Comparison.png [subcyclingMR] buildDir = . @@ -3848,7 +3804,6 @@ useOMP = 1 numthreads = 1 outputFile = spacecraft_charging_plt analysisRoutine = Examples/Physics_applications/spacecraft_charging/analysis.py -analysisOutputImage = min_phi_analysis.png [Point_of_contact_EB_3d] buildDir = . From 4ac5962ef4fce6e4f812ab64bed2f584ea412e91 Mon Sep 17 00:00:00 2001 From: Olga Shapoval <30510597+oshapoval@users.noreply.github.com> Date: Mon, 29 Jul 2024 17:07:18 -0700 Subject: [PATCH 11/12] Fix bug with ES solver and MR: `E_aux=E_fp` in `UpdateAuxilaryData` (#4922) * Removal of asserttion which prevented from usung the averaged PSATD algorithms with PML BC * Clean-up * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Fixed to arr_aux(j,k,l) = fine when ES solve is used * Removed temporary print statements * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Clean-up * Added CI test ElectrostaticSphereEB_RZ_MR_lev_1 to check the fields on the level=1 * Updated becnmarks for ElectrostaticSphereLabFrame_MR_emass_10 * Imported regular expression (re) in the analysis script. * Fixed typo * United two CI tests for different levels of MR in one test. * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Updated CI test ElectrostaticSphereEB_RZ_MR and the corresponding analysis script with smaller MR patch. * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Clean-up * Do deepcopy for lev>0 and collocated grid * Fix bugs to resolve failure of CI tests * Preserve plotfile output, update benchmark file * Working on CI test * Update benchmark of `ElectrostaticSphereEB_RZ_MR` * Initialize with value all `aux`, `cax` fields * Remove changes related to averaged Galilean PSATD with PML * Remove style changes (e.g., changes to empty lines) * Replace `MFIter`/`ParallelFor` loop with simple copy * Apply suggestions from code review * Revert part of the code to its previous, equivalent state * Removed DeepCopy & no need for ghost cells in temp phi_cp. * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Update Source/ablastr/fields/PoissonSolver.H * Add inline comments * Update analysis script * Use `TilingIfNotGPU`, `growntilebox` in E loop --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Edoardo Zoni Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com> Co-authored-by: Remi Lehe Co-authored-by: Weiqun Zhang --- .../electrostatic_sphere_eb/analysis_rz_mr.py | 99 ++++ .../electrostatic_sphere_eb/inputs_rz_mr | 1 + .../ElectrostaticSphereEB_RZ_MR.json | 8 +- ...ectrostaticSphereLabFrame_MR_emass_10.json | 12 +- Regression/WarpX-tests.ini | 5 +- Source/Parallelization/WarpXComm.cpp | 510 ++++++++++-------- Source/Parallelization/WarpXComm_K.H | 90 +++- Source/ablastr/fields/PoissonSolver.H | 21 +- 8 files changed, 508 insertions(+), 238 deletions(-) create mode 100755 Examples/Tests/electrostatic_sphere_eb/analysis_rz_mr.py diff --git a/Examples/Tests/electrostatic_sphere_eb/analysis_rz_mr.py b/Examples/Tests/electrostatic_sphere_eb/analysis_rz_mr.py new file mode 100755 index 00000000000..0b01b128362 --- /dev/null +++ b/Examples/Tests/electrostatic_sphere_eb/analysis_rz_mr.py @@ -0,0 +1,99 @@ +#!/usr/bin/env python + +# Copyright 2024 Olga Shapoval, Edoardo Zoni +# +# This file is part of WarpX. +# +# License: BSD-3-Clause-LBNL + +# This script tests the embedded boundary in RZ. +# A cylindrical surface (r=0.1) has a fixed potential 1 V. +# The outer surface has 0 V fixed. +# Thus the analytical solution has the form: +# phi(r) = A+B*log(r), Er(r) = -B/r. + +import os +import sys + +import numpy as np +from openpmd_viewer import OpenPMDTimeSeries + +sys.path.insert(1, '../../../../warpx/Regression/Checksum/') +import checksumAPI + +tolerance = 0.004 +print(f'tolerance = {tolerance}') + +fn = sys.argv[1] + +def find_first_non_zero_from_bottom_left(matrix): + for i in range(matrix.shape[0]): + for j in range(matrix.shape[1]): + if (matrix[i][j] != 0) and (matrix[i][j] != np.nan): + return (i, j) + return i, j + +def find_first_non_zero_from_upper_right(matrix): + for i in range(matrix.shape[0]-1, -1, -1): + for j in range(matrix.shape[1]-1, -1, -1): + if (matrix[i][j] != 0) and (matrix[i][j] != np.nan): + return (i, j) + return i,j + +def get_fields(ts, level): + if level == 0: + Er, info = ts.get_field('E', 'r', iteration=0) + phi, info = ts.get_field('phi', iteration=0) + else: + Er, info = ts.get_field(f'E_lvl{level}', 'r', iteration=0) + phi, info = ts.get_field(f'phi_lvl{level}', iteration=0) + return Er, phi, info + +def get_error_per_lev(ts,level): + Er, phi, info = get_fields(ts, level) + + nr_half = info.r.shape[0] // 2 + dr = info.dr + + Er_patch = Er[:,nr_half:] + phi_patch = phi[:,nr_half:] + r1 = info.r[nr_half:] + patch_left_lower_i, patch_left_lower_j = find_first_non_zero_from_bottom_left(Er_patch) + patch_right_upper_i, patch_right_upper_j = find_first_non_zero_from_upper_right(Er_patch) + + # phi and Er field on the MR patch + phi_sim = phi_patch[patch_left_lower_i:patch_right_upper_i+1, patch_left_lower_j:patch_right_upper_j+1] + Er_sim = Er_patch[patch_left_lower_i:patch_right_upper_i+1, patch_left_lower_j:patch_right_upper_j+1] + r = r1[patch_left_lower_j:patch_right_upper_j+1] + + B = 1.0/np.log(0.1/0.5) + A = -B*np.log(0.5) + + # outside EB and last cutcell + rmin = np.min(np.argwhere(r >= (0.1+dr))) + rmax = -1 + r = r[rmin:rmax] + phi_sim = phi_sim[:,rmin:rmax] + Er_sim = Er_sim[:,rmin:rmax] + + phi_theory = A + B*np.log(r) + phi_theory = np.tile(phi_theory, (phi_sim.shape[0],1)) + phi_error = np.max(np.abs(phi_theory-phi_sim) / np.abs(phi_theory)) + + Er_theory = -B/r + Er_theory = np.tile(Er_theory, (Er_sim.shape[0],1)) + Er_error = np.max(np.abs(Er_theory-Er_sim) / np.abs(Er_theory)) + + print(f'max error of phi[lev={level}]: {phi_error}') + print(f'max error of Er[lev={level}]: {Er_error}') + assert(phi_error < tolerance) + assert(Er_error < tolerance) + +ts = OpenPMDTimeSeries(fn) +level_fields = [field for field in ts.avail_fields if 'lvl' in field] +nlevels = 0 if level_fields == [] else int(level_fields[-1][-1]) +for level in range(nlevels+1): + get_error_per_lev(ts,level) + +test_name = os.path.split(os.getcwd())[1] +checksumAPI.evaluate_checksum(test_name, fn, output_format="openpmd") diff --git a/Examples/Tests/electrostatic_sphere_eb/inputs_rz_mr b/Examples/Tests/electrostatic_sphere_eb/inputs_rz_mr index 3bea63d76fb..722fc916416 100644 --- a/Examples/Tests/electrostatic_sphere_eb/inputs_rz_mr +++ b/Examples/Tests/electrostatic_sphere_eb/inputs_rz_mr @@ -30,3 +30,4 @@ diagnostics.diags_names = diag1 diag1.intervals = 1 diag1.diag_type = Full diag1.fields_to_plot = Er phi +diag1.format = openpmd diff --git a/Regression/Checksum/benchmarks_json/ElectrostaticSphereEB_RZ_MR.json b/Regression/Checksum/benchmarks_json/ElectrostaticSphereEB_RZ_MR.json index ffa8d68f9d9..6bbfce0e3b3 100644 --- a/Regression/Checksum/benchmarks_json/ElectrostaticSphereEB_RZ_MR.json +++ b/Regression/Checksum/benchmarks_json/ElectrostaticSphereEB_RZ_MR.json @@ -1,10 +1,10 @@ { "lev=0": { - "Er": 8487.661571739109, - "phi": 2036.0428085225362 + "Er": 16975.32314347822, + "phi": 4072.085617045073 }, "lev=1": { - "Er": 19519.172334977942, - "phi": 3291.0262856782897 + "Er": 26818.189739547757, + "phi_lvl1": 8731.176548788893 } } diff --git a/Regression/Checksum/benchmarks_json/ElectrostaticSphereLabFrame_MR_emass_10.json b/Regression/Checksum/benchmarks_json/ElectrostaticSphereLabFrame_MR_emass_10.json index 024127a1bd2..21d5208c59a 100644 --- a/Regression/Checksum/benchmarks_json/ElectrostaticSphereLabFrame_MR_emass_10.json +++ b/Regression/Checksum/benchmarks_json/ElectrostaticSphereLabFrame_MR_emass_10.json @@ -6,15 +6,15 @@ "rho": 0.0 }, "lev=1": { - "Ex": 14.281015560380963, - "Ey": 14.281015560380965, - "Ez": 14.281015560380965, + "Ex": 7.170105936287823, + "Ey": 7.17010593628782, + "Ez": 7.170105936287821, "rho": 2.6092568008333786e-10 }, "electron": { - "particle_momentum_x": 1.80842228672388e-24, - "particle_momentum_y": 1.8084222867238806e-24, - "particle_momentum_z": 1.7598771525647628e-24, + "particle_momentum_x": 9.257577597262615e-25, + "particle_momentum_y": 9.257577597262618e-25, + "particle_momentum_z": 9.257577597262624e-25, "particle_position_x": 327.46875, "particle_position_y": 327.46875, "particle_position_z": 327.46875, diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index a1dc0a168f7..8048318de7a 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -483,7 +483,7 @@ analysisRoutine = Examples/Tests/electrostatic_sphere_eb/analysis_rz.py [ElectrostaticSphereEB_RZ_MR] buildDir = . inputFile = Examples/Tests/electrostatic_sphere_eb/inputs_rz_mr -runtime_params = warpx.abort_on_warning_threshold = medium +runtime_params = warpx.abort_on_warning_threshold = medium amr.ref_ratio_vect = 2 2 2 dim = 2 addToCompileString = USE_EB=TRUE USE_RZ=TRUE cmakeSetupOpts = -DWarpX_DIMS=RZ -DWarpX_EB=ON @@ -492,7 +492,8 @@ useMPI = 1 numprocs = 2 useOMP = 1 numthreads = 1 -analysisRoutine = Examples/Tests/electrostatic_sphere_eb/analysis_rz.py +outputFile = ElectrostaticSphereEB_RZ_MR_plt +analysisRoutine = Examples/Tests/electrostatic_sphere_eb/analysis_rz_mr.py [ElectrostaticSphereLabFrame] buildDir = . diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index e7df489236e..2887bd4d056 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -171,136 +171,191 @@ WarpX::UpdateAuxilaryDataStagToNodal () // Bfield { - Array,3> Btmp; - if (Bfield_cax[lev][0]) { - for (int i = 0; i < 3; ++i) { - Btmp[i] = std::make_unique( - *Bfield_cax[lev][i], amrex::make_alias, 0, 1); + if (electromagnetic_solver_id != ElectromagneticSolverAlgo::None) { + Array,3> Btmp; + if (Bfield_cax[lev][0]) { + for (int i = 0; i < 3; ++i) { + Btmp[i] = std::make_unique( + *Bfield_cax[lev][i], amrex::make_alias, 0, 1); + } + } else { + const IntVect ngtmp = Bfield_aux[lev-1][0]->nGrowVect(); + for (int i = 0; i < 3; ++i) { + Btmp[i] = std::make_unique(cnba, dm, 1, ngtmp); + } } - } else { - const IntVect ngtmp = Bfield_aux[lev-1][0]->nGrowVect(); + Btmp[0]->setVal(0.0); + Btmp[1]->setVal(0.0); + Btmp[2]->setVal(0.0); + // ParallelCopy from coarse level for (int i = 0; i < 3; ++i) { - Btmp[i] = std::make_unique(cnba, dm, 1, ngtmp); + const IntVect ng = Btmp[i]->nGrowVect(); + // Guard cells may not be up to date beyond ng_FieldGather + const amrex::IntVect& ng_src = guard_cells.ng_FieldGather; + // Copy Bfield_aux to Btmp, using up to ng_src (=ng_FieldGather) guard cells from + // Bfield_aux and filling up to ng (=nGrow) guard cells in Btmp + ablastr::utils::communication::ParallelCopy(*Btmp[i], *Bfield_aux[lev - 1][i], 0, 0, 1, + ng_src, ng, WarpX::do_single_precision_comms, cperiod); } - } - Btmp[0]->setVal(0.0); - Btmp[1]->setVal(0.0); - Btmp[2]->setVal(0.0); - // ParallelCopy from coarse level - for (int i = 0; i < 3; ++i) { - const IntVect ng = Btmp[i]->nGrowVect(); - // Guard cells may not be up to date beyond ng_FieldGather - const amrex::IntVect& ng_src = guard_cells.ng_FieldGather; - // Copy Bfield_aux to Btmp, using up to ng_src (=ng_FieldGather) guard cells from - // Bfield_aux and filling up to ng (=nGrow) guard cells in Btmp - ablastr::utils::communication::ParallelCopy(*Btmp[i], *Bfield_aux[lev - 1][i], 0, 0, 1, - ng_src, ng, WarpX::do_single_precision_comms, cperiod); - } - const amrex::IntVect& refinement_ratio = refRatio(lev-1); + const amrex::IntVect& refinement_ratio = refRatio(lev-1); - const amrex::IntVect& Bx_fp_stag = Bfield_fp[lev][0]->ixType().toIntVect(); - const amrex::IntVect& By_fp_stag = Bfield_fp[lev][1]->ixType().toIntVect(); - const amrex::IntVect& Bz_fp_stag = Bfield_fp[lev][2]->ixType().toIntVect(); + const amrex::IntVect& Bx_fp_stag = Bfield_fp[lev][0]->ixType().toIntVect(); + const amrex::IntVect& By_fp_stag = Bfield_fp[lev][1]->ixType().toIntVect(); + const amrex::IntVect& Bz_fp_stag = Bfield_fp[lev][2]->ixType().toIntVect(); - const amrex::IntVect& Bx_cp_stag = Bfield_cp[lev][0]->ixType().toIntVect(); - const amrex::IntVect& By_cp_stag = Bfield_cp[lev][1]->ixType().toIntVect(); - const amrex::IntVect& Bz_cp_stag = Bfield_cp[lev][2]->ixType().toIntVect(); + const amrex::IntVect& Bx_cp_stag = Bfield_cp[lev][0]->ixType().toIntVect(); + const amrex::IntVect& By_cp_stag = Bfield_cp[lev][1]->ixType().toIntVect(); + const amrex::IntVect& Bz_cp_stag = Bfield_cp[lev][2]->ixType().toIntVect(); #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif - for (MFIter mfi(*Bfield_aux[lev][0], TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - Array4 const& bx_aux = Bfield_aux[lev][0]->array(mfi); - Array4 const& by_aux = Bfield_aux[lev][1]->array(mfi); - Array4 const& bz_aux = Bfield_aux[lev][2]->array(mfi); - Array4 const& bx_fp = Bfield_fp[lev][0]->const_array(mfi); - Array4 const& by_fp = Bfield_fp[lev][1]->const_array(mfi); - Array4 const& bz_fp = Bfield_fp[lev][2]->const_array(mfi); - Array4 const& bx_cp = Bfield_cp[lev][0]->const_array(mfi); - Array4 const& by_cp = Bfield_cp[lev][1]->const_array(mfi); - Array4 const& bz_cp = Bfield_cp[lev][2]->const_array(mfi); - Array4 const& bx_c = Btmp[0]->const_array(mfi); - Array4 const& by_c = Btmp[1]->const_array(mfi); - Array4 const& bz_c = Btmp[2]->const_array(mfi); - - const Box& bx = mfi.growntilebox(); - amrex::ParallelFor(bx, - [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept + for (MFIter mfi(*Bfield_aux[lev][0], TilingIfNotGPU()); mfi.isValid(); ++mfi) { - warpx_interp(j, k, l, bx_aux, bx_fp, bx_cp, bx_c, Bx_fp_stag, Bx_cp_stag, refinement_ratio); - warpx_interp(j, k, l, by_aux, by_fp, by_cp, by_c, By_fp_stag, By_cp_stag, refinement_ratio); - warpx_interp(j, k, l, bz_aux, bz_fp, bz_cp, bz_c, Bz_fp_stag, Bz_cp_stag, refinement_ratio); - }); + Array4 const& bx_aux = Bfield_aux[lev][0]->array(mfi); + Array4 const& by_aux = Bfield_aux[lev][1]->array(mfi); + Array4 const& bz_aux = Bfield_aux[lev][2]->array(mfi); + Array4 const& bx_fp = Bfield_fp[lev][0]->const_array(mfi); + Array4 const& by_fp = Bfield_fp[lev][1]->const_array(mfi); + Array4 const& bz_fp = Bfield_fp[lev][2]->const_array(mfi); + Array4 const& bx_cp = Bfield_cp[lev][0]->const_array(mfi); + Array4 const& by_cp = Bfield_cp[lev][1]->const_array(mfi); + Array4 const& bz_cp = Bfield_cp[lev][2]->const_array(mfi); + Array4 const& bx_c = Btmp[0]->const_array(mfi); + Array4 const& by_c = Btmp[1]->const_array(mfi); + Array4 const& bz_c = Btmp[2]->const_array(mfi); + + const Box& bx = mfi.growntilebox(); + amrex::ParallelFor(bx, + [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept + { + warpx_interp(j, k, l, bx_aux, bx_fp, bx_cp, bx_c, Bx_fp_stag, Bx_cp_stag, refinement_ratio); + warpx_interp(j, k, l, by_aux, by_fp, by_cp, by_c, By_fp_stag, By_cp_stag, refinement_ratio); + warpx_interp(j, k, l, bz_aux, bz_fp, bz_cp, bz_c, Bz_fp_stag, Bz_cp_stag, refinement_ratio); + }); + } + } + else { // electrostatic + const amrex::IntVect& Bx_fp_stag = Bfield_fp[lev][0]->ixType().toIntVect(); + const amrex::IntVect& By_fp_stag = Bfield_fp[lev][1]->ixType().toIntVect(); + const amrex::IntVect& Bz_fp_stag = Bfield_fp[lev][2]->ixType().toIntVect(); +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(*Bfield_aux[lev][0], TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + Array4 const& bx_aux = Bfield_aux[lev][0]->array(mfi); + Array4 const& by_aux = Bfield_aux[lev][1]->array(mfi); + Array4 const& bz_aux = Bfield_aux[lev][2]->array(mfi); + Array4 const& bx_fp = Bfield_fp[lev][0]->const_array(mfi); + Array4 const& by_fp = Bfield_fp[lev][1]->const_array(mfi); + Array4 const& bz_fp = Bfield_fp[lev][2]->const_array(mfi); + + const Box& bx = mfi.growntilebox(); + amrex::ParallelFor(bx, + [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept + { + warpx_interp(j, k, l, bx_aux, bx_fp, Bx_fp_stag); + warpx_interp(j, k, l, by_aux, by_fp, By_fp_stag); + warpx_interp(j, k, l, bz_aux, bz_fp, Bz_fp_stag); + }); + } } } - // Efield { - Array,3> Etmp; - if (Efield_cax[lev][0]) { - for (int i = 0; i < 3; ++i) { - Etmp[i] = std::make_unique( - *Efield_cax[lev][i], amrex::make_alias, 0, 1); + if (electromagnetic_solver_id != ElectromagneticSolverAlgo::None) { + Array,3> Etmp; + if (Efield_cax[lev][0]) { + for (int i = 0; i < 3; ++i) { + Etmp[i] = std::make_unique( + *Efield_cax[lev][i], amrex::make_alias, 0, 1); + } + } else { + const IntVect ngtmp = Efield_aux[lev-1][0]->nGrowVect(); + for (int i = 0; i < 3; ++i) { + Etmp[i] = std::make_unique( + cnba, dm, 1, ngtmp); + } } - } else { - const IntVect ngtmp = Efield_aux[lev-1][0]->nGrowVect(); + Etmp[0]->setVal(0.0); + Etmp[1]->setVal(0.0); + Etmp[2]->setVal(0.0); + // ParallelCopy from coarse level for (int i = 0; i < 3; ++i) { - Etmp[i] = std::make_unique( - cnba, dm, 1, ngtmp); + const IntVect ng = Etmp[i]->nGrowVect(); + // Guard cells may not be up to date beyond ng_FieldGather + const amrex::IntVect& ng_src = guard_cells.ng_FieldGather; + // Copy Efield_aux to Etmp, using up to ng_src (=ng_FieldGather) guard cells from + // Efield_aux and filling up to ng (=nGrow) guard cells in Etmp + ablastr::utils::communication::ParallelCopy(*Etmp[i], *Efield_aux[lev - 1][i], 0, 0, 1, + ng_src, ng, WarpX::do_single_precision_comms, cperiod); } - } - Etmp[0]->setVal(0.0); - Etmp[1]->setVal(0.0); - Etmp[2]->setVal(0.0); - // ParallelCopy from coarse level - for (int i = 0; i < 3; ++i) { - const IntVect ng = Etmp[i]->nGrowVect(); - // Guard cells may not be up to date beyond ng_FieldGather - const amrex::IntVect& ng_src = guard_cells.ng_FieldGather; - // Copy Efield_aux to Etmp, using up to ng_src (=ng_FieldGather) guard cells from - // Efield_aux and filling up to ng (=nGrow) guard cells in Etmp - ablastr::utils::communication::ParallelCopy(*Etmp[i], *Efield_aux[lev - 1][i], 0, 0, 1, - ng_src, ng, WarpX::do_single_precision_comms, cperiod); - } - const amrex::IntVect& refinement_ratio = refRatio(lev-1); + const amrex::IntVect& refinement_ratio = refRatio(lev-1); - const amrex::IntVect& Ex_fp_stag = Efield_fp[lev][0]->ixType().toIntVect(); - const amrex::IntVect& Ey_fp_stag = Efield_fp[lev][1]->ixType().toIntVect(); - const amrex::IntVect& Ez_fp_stag = Efield_fp[lev][2]->ixType().toIntVect(); + const amrex::IntVect& Ex_fp_stag = Efield_fp[lev][0]->ixType().toIntVect(); + const amrex::IntVect& Ey_fp_stag = Efield_fp[lev][1]->ixType().toIntVect(); + const amrex::IntVect& Ez_fp_stag = Efield_fp[lev][2]->ixType().toIntVect(); - const amrex::IntVect& Ex_cp_stag = Efield_cp[lev][0]->ixType().toIntVect(); - const amrex::IntVect& Ey_cp_stag = Efield_cp[lev][1]->ixType().toIntVect(); - const amrex::IntVect& Ez_cp_stag = Efield_cp[lev][2]->ixType().toIntVect(); + const amrex::IntVect& Ex_cp_stag = Efield_cp[lev][0]->ixType().toIntVect(); + const amrex::IntVect& Ey_cp_stag = Efield_cp[lev][1]->ixType().toIntVect(); + const amrex::IntVect& Ez_cp_stag = Efield_cp[lev][2]->ixType().toIntVect(); #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif - for (MFIter mfi(*Efield_aux[lev][0]); mfi.isValid(); ++mfi) - { - Array4 const& ex_aux = Efield_aux[lev][0]->array(mfi); - Array4 const& ey_aux = Efield_aux[lev][1]->array(mfi); - Array4 const& ez_aux = Efield_aux[lev][2]->array(mfi); - Array4 const& ex_fp = Efield_fp[lev][0]->const_array(mfi); - Array4 const& ey_fp = Efield_fp[lev][1]->const_array(mfi); - Array4 const& ez_fp = Efield_fp[lev][2]->const_array(mfi); - Array4 const& ex_cp = Efield_cp[lev][0]->const_array(mfi); - Array4 const& ey_cp = Efield_cp[lev][1]->const_array(mfi); - Array4 const& ez_cp = Efield_cp[lev][2]->const_array(mfi); - Array4 const& ex_c = Etmp[0]->const_array(mfi); - Array4 const& ey_c = Etmp[1]->const_array(mfi); - Array4 const& ez_c = Etmp[2]->const_array(mfi); - - const Box& bx = mfi.fabbox(); - amrex::ParallelFor(bx, - [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept + for (MFIter mfi(*Efield_aux[lev][0]); mfi.isValid(); ++mfi) { - warpx_interp(j, k, l, ex_aux, ex_fp, ex_cp, ex_c, Ex_fp_stag, Ex_cp_stag, refinement_ratio); - warpx_interp(j, k, l, ey_aux, ey_fp, ey_cp, ey_c, Ey_fp_stag, Ey_cp_stag, refinement_ratio); - warpx_interp(j, k, l, ez_aux, ez_fp, ez_cp, ez_c, Ez_fp_stag, Ez_cp_stag, refinement_ratio); - }); + Array4 const& ex_aux = Efield_aux[lev][0]->array(mfi); + Array4 const& ey_aux = Efield_aux[lev][1]->array(mfi); + Array4 const& ez_aux = Efield_aux[lev][2]->array(mfi); + Array4 const& ex_fp = Efield_fp[lev][0]->const_array(mfi); + Array4 const& ey_fp = Efield_fp[lev][1]->const_array(mfi); + Array4 const& ez_fp = Efield_fp[lev][2]->const_array(mfi); + Array4 const& ex_cp = Efield_cp[lev][0]->const_array(mfi); + Array4 const& ey_cp = Efield_cp[lev][1]->const_array(mfi); + Array4 const& ez_cp = Efield_cp[lev][2]->const_array(mfi); + Array4 const& ex_c = Etmp[0]->const_array(mfi); + Array4 const& ey_c = Etmp[1]->const_array(mfi); + Array4 const& ez_c = Etmp[2]->const_array(mfi); + + const Box& bx = mfi.fabbox(); + amrex::ParallelFor(bx, + [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept + { + warpx_interp(j, k, l, ex_aux, ex_fp, ex_cp, ex_c, Ex_fp_stag, Ex_cp_stag, refinement_ratio); + warpx_interp(j, k, l, ey_aux, ey_fp, ey_cp, ey_c, Ey_fp_stag, Ey_cp_stag, refinement_ratio); + warpx_interp(j, k, l, ez_aux, ez_fp, ez_cp, ez_c, Ez_fp_stag, Ez_cp_stag, refinement_ratio); + }); + } + } + else { // electrostatic + const amrex::IntVect& Ex_fp_stag = Efield_fp[lev][0]->ixType().toIntVect(); + const amrex::IntVect& Ey_fp_stag = Efield_fp[lev][1]->ixType().toIntVect(); + const amrex::IntVect& Ez_fp_stag = Efield_fp[lev][2]->ixType().toIntVect(); +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(*Efield_aux[lev][0], TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + Array4 const& ex_aux = Efield_aux[lev][0]->array(mfi); + Array4 const& ey_aux = Efield_aux[lev][1]->array(mfi); + Array4 const& ez_aux = Efield_aux[lev][2]->array(mfi); + Array4 const& ex_fp = Efield_fp[lev][0]->const_array(mfi); + Array4 const& ey_fp = Efield_fp[lev][1]->const_array(mfi); + Array4 const& ez_fp = Efield_fp[lev][2]->const_array(mfi); + + const Box& bx = mfi.growntilebox(); + amrex::ParallelFor(bx, + [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept + { + warpx_interp(j, k, l, ex_aux, ex_fp, Ex_fp_stag); + warpx_interp(j, k, l, ey_aux, ey_fp, Ey_fp_stag); + warpx_interp(j, k, l, ez_aux, ez_fp, Ez_fp_stag); + }); + } } } } @@ -341,141 +396,158 @@ WarpX::UpdateAuxilaryDataSameType () // B field { - MultiFab dBx(Bfield_cp[lev][0]->boxArray(), dm, Bfield_cp[lev][0]->nComp(), ng); - MultiFab dBy(Bfield_cp[lev][1]->boxArray(), dm, Bfield_cp[lev][1]->nComp(), ng); - MultiFab dBz(Bfield_cp[lev][2]->boxArray(), dm, Bfield_cp[lev][2]->nComp(), ng); - dBx.setVal(0.0); - dBy.setVal(0.0); - dBz.setVal(0.0); - - // Copy Bfield_aux to the dB MultiFabs, using up to ng_src (=ng_FieldGather) guard - // cells from Bfield_aux and filling up to ng (=nGrow) guard cells in the dB MultiFabs - - ablastr::utils::communication::ParallelCopy(dBx, *Bfield_aux[lev - 1][0], 0, 0, - Bfield_aux[lev - 1][0]->nComp(), ng_src, ng, WarpX::do_single_precision_comms, - crse_period); - ablastr::utils::communication::ParallelCopy(dBy, *Bfield_aux[lev - 1][1], 0, 0, - Bfield_aux[lev - 1][1]->nComp(), ng_src, ng, WarpX::do_single_precision_comms, - crse_period); - ablastr::utils::communication::ParallelCopy(dBz, *Bfield_aux[lev - 1][2], 0, 0, - Bfield_aux[lev - 1][2]->nComp(), ng_src, ng, WarpX::do_single_precision_comms, - crse_period); - - if (Bfield_cax[lev][0]) + if (electromagnetic_solver_id != ElectromagneticSolverAlgo::None) { - MultiFab::Copy(*Bfield_cax[lev][0], dBx, 0, 0, Bfield_cax[lev][0]->nComp(), ng); - MultiFab::Copy(*Bfield_cax[lev][1], dBy, 0, 0, Bfield_cax[lev][1]->nComp(), ng); - MultiFab::Copy(*Bfield_cax[lev][2], dBz, 0, 0, Bfield_cax[lev][2]->nComp(), ng); - } - MultiFab::Subtract(dBx, *Bfield_cp[lev][0], 0, 0, Bfield_cp[lev][0]->nComp(), ng); - MultiFab::Subtract(dBy, *Bfield_cp[lev][1], 0, 0, Bfield_cp[lev][1]->nComp(), ng); - MultiFab::Subtract(dBz, *Bfield_cp[lev][2], 0, 0, Bfield_cp[lev][2]->nComp(), ng); + MultiFab dBx(Bfield_cp[lev][0]->boxArray(), dm, Bfield_cp[lev][0]->nComp(), ng); + MultiFab dBy(Bfield_cp[lev][1]->boxArray(), dm, Bfield_cp[lev][1]->nComp(), ng); + MultiFab dBz(Bfield_cp[lev][2]->boxArray(), dm, Bfield_cp[lev][2]->nComp(), ng); + dBx.setVal(0.0); + dBy.setVal(0.0); + dBz.setVal(0.0); + + // Copy Bfield_aux to the dB MultiFabs, using up to ng_src (=ng_FieldGather) guard + // cells from Bfield_aux and filling up to ng (=nGrow) guard cells in the dB MultiFabs + + ablastr::utils::communication::ParallelCopy(dBx, *Bfield_aux[lev - 1][0], 0, 0, + Bfield_aux[lev - 1][0]->nComp(), ng_src, ng, WarpX::do_single_precision_comms, + crse_period); + ablastr::utils::communication::ParallelCopy(dBy, *Bfield_aux[lev - 1][1], 0, 0, + Bfield_aux[lev - 1][1]->nComp(), ng_src, ng, WarpX::do_single_precision_comms, + crse_period); + ablastr::utils::communication::ParallelCopy(dBz, *Bfield_aux[lev - 1][2], 0, 0, + Bfield_aux[lev - 1][2]->nComp(), ng_src, ng, WarpX::do_single_precision_comms, + crse_period); + + if (Bfield_cax[lev][0]) + { + MultiFab::Copy(*Bfield_cax[lev][0], dBx, 0, 0, Bfield_cax[lev][0]->nComp(), ng); + MultiFab::Copy(*Bfield_cax[lev][1], dBy, 0, 0, Bfield_cax[lev][1]->nComp(), ng); + MultiFab::Copy(*Bfield_cax[lev][2], dBz, 0, 0, Bfield_cax[lev][2]->nComp(), ng); + } + MultiFab::Subtract(dBx, *Bfield_cp[lev][0], 0, 0, Bfield_cp[lev][0]->nComp(), ng); + MultiFab::Subtract(dBy, *Bfield_cp[lev][1], 0, 0, Bfield_cp[lev][1]->nComp(), ng); + MultiFab::Subtract(dBz, *Bfield_cp[lev][2], 0, 0, Bfield_cp[lev][2]->nComp(), ng); - const amrex::IntVect& refinement_ratio = refRatio(lev-1); + const amrex::IntVect& refinement_ratio = refRatio(lev-1); - const amrex::IntVect& Bx_stag = Bfield_aux[lev-1][0]->ixType().toIntVect(); - const amrex::IntVect& By_stag = Bfield_aux[lev-1][1]->ixType().toIntVect(); - const amrex::IntVect& Bz_stag = Bfield_aux[lev-1][2]->ixType().toIntVect(); + const amrex::IntVect& Bx_stag = Bfield_aux[lev-1][0]->ixType().toIntVect(); + const amrex::IntVect& By_stag = Bfield_aux[lev-1][1]->ixType().toIntVect(); + const amrex::IntVect& Bz_stag = Bfield_aux[lev-1][2]->ixType().toIntVect(); #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif - for (MFIter mfi(*Bfield_aux[lev][0]); mfi.isValid(); ++mfi) - { - Array4 const& bx_aux = Bfield_aux[lev][0]->array(mfi); - Array4 const& by_aux = Bfield_aux[lev][1]->array(mfi); - Array4 const& bz_aux = Bfield_aux[lev][2]->array(mfi); - Array4 const& bx_fp = Bfield_fp[lev][0]->const_array(mfi); - Array4 const& by_fp = Bfield_fp[lev][1]->const_array(mfi); - Array4 const& bz_fp = Bfield_fp[lev][2]->const_array(mfi); - Array4 const& bx_c = dBx.const_array(mfi); - Array4 const& by_c = dBy.const_array(mfi); - Array4 const& bz_c = dBz.const_array(mfi); - - amrex::ParallelFor(Box(bx_aux), Box(by_aux), Box(bz_aux), - [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept - { - warpx_interp(j, k, l, bx_aux, bx_fp, bx_c, Bx_stag, refinement_ratio); - }, - [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept - { - warpx_interp(j, k, l, by_aux, by_fp, by_c, By_stag, refinement_ratio); - }, - [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept + for (MFIter mfi(*Bfield_aux[lev][0]); mfi.isValid(); ++mfi) { - warpx_interp(j, k, l, bz_aux, bz_fp, bz_c, Bz_stag, refinement_ratio); - }); + Array4 const& bx_aux = Bfield_aux[lev][0]->array(mfi); + Array4 const& by_aux = Bfield_aux[lev][1]->array(mfi); + Array4 const& bz_aux = Bfield_aux[lev][2]->array(mfi); + Array4 const& bx_fp = Bfield_fp[lev][0]->const_array(mfi); + Array4 const& by_fp = Bfield_fp[lev][1]->const_array(mfi); + Array4 const& bz_fp = Bfield_fp[lev][2]->const_array(mfi); + Array4 const& bx_c = dBx.const_array(mfi); + Array4 const& by_c = dBy.const_array(mfi); + Array4 const& bz_c = dBz.const_array(mfi); + + amrex::ParallelFor(Box(bx_aux), Box(by_aux), Box(bz_aux), + [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept + { + warpx_interp(j, k, l, bx_aux, bx_fp, bx_c, Bx_stag, refinement_ratio); + }, + [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept + { + warpx_interp(j, k, l, by_aux, by_fp, by_c, By_stag, refinement_ratio); + }, + [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept + { + warpx_interp(j, k, l, bz_aux, bz_fp, bz_c, Bz_stag, refinement_ratio); + }); + } + } + else // electrostatic + { + MultiFab::Copy(*Bfield_aux[lev][0], *Bfield_fp[lev][0], 0, 0, Bfield_aux[lev][0]->nComp(), Bfield_aux[lev][0]->nGrowVect()); + MultiFab::Copy(*Bfield_aux[lev][1], *Bfield_fp[lev][1], 0, 0, Bfield_aux[lev][1]->nComp(), Bfield_aux[lev][1]->nGrowVect()); + MultiFab::Copy(*Bfield_aux[lev][2], *Bfield_fp[lev][2], 0, 0, Bfield_aux[lev][2]->nComp(), Bfield_aux[lev][2]->nGrowVect()); } } - // E field { - MultiFab dEx(Efield_cp[lev][0]->boxArray(), dm, Efield_cp[lev][0]->nComp(), ng); - MultiFab dEy(Efield_cp[lev][1]->boxArray(), dm, Efield_cp[lev][1]->nComp(), ng); - MultiFab dEz(Efield_cp[lev][2]->boxArray(), dm, Efield_cp[lev][2]->nComp(), ng); - dEx.setVal(0.0); - dEy.setVal(0.0); - dEz.setVal(0.0); - - // Copy Efield_aux to the dE MultiFabs, using up to ng_src (=ng_FieldGather) guard - // cells from Efield_aux and filling up to ng (=nGrow) guard cells in the dE MultiFabs - ablastr::utils::communication::ParallelCopy(dEx, *Efield_aux[lev - 1][0], 0, 0, - Efield_aux[lev - 1][0]->nComp(), ng_src, ng, - WarpX::do_single_precision_comms, - crse_period); - ablastr::utils::communication::ParallelCopy(dEy, *Efield_aux[lev - 1][1], 0, 0, - Efield_aux[lev - 1][1]->nComp(), ng_src, ng, - WarpX::do_single_precision_comms, - crse_period); - ablastr::utils::communication::ParallelCopy(dEz, *Efield_aux[lev - 1][2], 0, 0, - Efield_aux[lev - 1][2]->nComp(), ng_src, ng, - WarpX::do_single_precision_comms, - crse_period); - - if (Efield_cax[lev][0]) + if (electromagnetic_solver_id != ElectromagneticSolverAlgo::None) { - MultiFab::Copy(*Efield_cax[lev][0], dEx, 0, 0, Efield_cax[lev][0]->nComp(), ng); - MultiFab::Copy(*Efield_cax[lev][1], dEy, 0, 0, Efield_cax[lev][1]->nComp(), ng); - MultiFab::Copy(*Efield_cax[lev][2], dEz, 0, 0, Efield_cax[lev][2]->nComp(), ng); - } - MultiFab::Subtract(dEx, *Efield_cp[lev][0], 0, 0, Efield_cp[lev][0]->nComp(), ng); - MultiFab::Subtract(dEy, *Efield_cp[lev][1], 0, 0, Efield_cp[lev][1]->nComp(), ng); - MultiFab::Subtract(dEz, *Efield_cp[lev][2], 0, 0, Efield_cp[lev][2]->nComp(), ng); + MultiFab dEx(Efield_cp[lev][0]->boxArray(), dm, Efield_cp[lev][0]->nComp(), ng); + MultiFab dEy(Efield_cp[lev][1]->boxArray(), dm, Efield_cp[lev][1]->nComp(), ng); + MultiFab dEz(Efield_cp[lev][2]->boxArray(), dm, Efield_cp[lev][2]->nComp(), ng); + dEx.setVal(0.0); + dEy.setVal(0.0); + dEz.setVal(0.0); + + // Copy Efield_aux to the dE MultiFabs, using up to ng_src (=ng_FieldGather) guard + // cells from Efield_aux and filling up to ng (=nGrow) guard cells in the dE MultiFabs + ablastr::utils::communication::ParallelCopy(dEx, *Efield_aux[lev - 1][0], 0, 0, + Efield_aux[lev - 1][0]->nComp(), ng_src, ng, + WarpX::do_single_precision_comms, + crse_period); + ablastr::utils::communication::ParallelCopy(dEy, *Efield_aux[lev - 1][1], 0, 0, + Efield_aux[lev - 1][1]->nComp(), ng_src, ng, + WarpX::do_single_precision_comms, + crse_period); + ablastr::utils::communication::ParallelCopy(dEz, *Efield_aux[lev - 1][2], 0, 0, + Efield_aux[lev - 1][2]->nComp(), ng_src, ng, + WarpX::do_single_precision_comms, + crse_period); + + if (Efield_cax[lev][0]) + { + MultiFab::Copy(*Efield_cax[lev][0], dEx, 0, 0, Efield_cax[lev][0]->nComp(), ng); + MultiFab::Copy(*Efield_cax[lev][1], dEy, 0, 0, Efield_cax[lev][1]->nComp(), ng); + MultiFab::Copy(*Efield_cax[lev][2], dEz, 0, 0, Efield_cax[lev][2]->nComp(), ng); + } + MultiFab::Subtract(dEx, *Efield_cp[lev][0], 0, 0, Efield_cp[lev][0]->nComp(), ng); + MultiFab::Subtract(dEy, *Efield_cp[lev][1], 0, 0, Efield_cp[lev][1]->nComp(), ng); + MultiFab::Subtract(dEz, *Efield_cp[lev][2], 0, 0, Efield_cp[lev][2]->nComp(), ng); - const amrex::IntVect& refinement_ratio = refRatio(lev-1); + const amrex::IntVect& refinement_ratio = refRatio(lev-1); - const amrex::IntVect& Ex_stag = Efield_aux[lev-1][0]->ixType().toIntVect(); - const amrex::IntVect& Ey_stag = Efield_aux[lev-1][1]->ixType().toIntVect(); - const amrex::IntVect& Ez_stag = Efield_aux[lev-1][2]->ixType().toIntVect(); + const amrex::IntVect& Ex_stag = Efield_aux[lev-1][0]->ixType().toIntVect(); + const amrex::IntVect& Ey_stag = Efield_aux[lev-1][1]->ixType().toIntVect(); + const amrex::IntVect& Ez_stag = Efield_aux[lev-1][2]->ixType().toIntVect(); #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif - for (MFIter mfi(*Efield_aux[lev][0]); mfi.isValid(); ++mfi) - { - Array4 const& ex_aux = Efield_aux[lev][0]->array(mfi); - Array4 const& ey_aux = Efield_aux[lev][1]->array(mfi); - Array4 const& ez_aux = Efield_aux[lev][2]->array(mfi); - Array4 const& ex_fp = Efield_fp[lev][0]->const_array(mfi); - Array4 const& ey_fp = Efield_fp[lev][1]->const_array(mfi); - Array4 const& ez_fp = Efield_fp[lev][2]->const_array(mfi); - Array4 const& ex_c = dEx.const_array(mfi); - Array4 const& ey_c = dEy.const_array(mfi); - Array4 const& ez_c = dEz.const_array(mfi); - - amrex::ParallelFor(Box(ex_aux), Box(ey_aux), Box(ez_aux), - [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept + for (MFIter mfi(*Efield_aux[lev][0]); mfi.isValid(); ++mfi) { - warpx_interp(j, k, l, ex_aux, ex_fp, ex_c, Ex_stag, refinement_ratio); - }, - [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept - { - warpx_interp(j, k, l, ey_aux, ey_fp, ey_c, Ey_stag, refinement_ratio); - }, - [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept - { - warpx_interp(j, k, l, ez_aux, ez_fp, ez_c, Ez_stag, refinement_ratio); - }); + Array4 const& ex_aux = Efield_aux[lev][0]->array(mfi); + Array4 const& ey_aux = Efield_aux[lev][1]->array(mfi); + Array4 const& ez_aux = Efield_aux[lev][2]->array(mfi); + Array4 const& ex_fp = Efield_fp[lev][0]->const_array(mfi); + Array4 const& ey_fp = Efield_fp[lev][1]->const_array(mfi); + Array4 const& ez_fp = Efield_fp[lev][2]->const_array(mfi); + Array4 const& ex_c = dEx.const_array(mfi); + Array4 const& ey_c = dEy.const_array(mfi); + Array4 const& ez_c = dEz.const_array(mfi); + + amrex::ParallelFor(Box(ex_aux), Box(ey_aux), Box(ez_aux), + [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept + { + warpx_interp(j, k, l, ex_aux, ex_fp, ex_c, Ex_stag, refinement_ratio); + }, + [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept + { + warpx_interp(j, k, l, ey_aux, ey_fp, ey_c, Ey_stag, refinement_ratio); + }, + [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept + { + warpx_interp(j, k, l, ez_aux, ez_fp, ez_c, Ez_stag, refinement_ratio); + }); + } + } + else // electrostatic + { + MultiFab::Copy(*Efield_aux[lev][0], *Efield_fp[lev][0], 0, 0, Efield_aux[lev][0]->nComp(), Efield_aux[lev][0]->nGrowVect()); + MultiFab::Copy(*Efield_aux[lev][1], *Efield_fp[lev][1], 0, 0, Efield_aux[lev][1]->nComp(), Efield_aux[lev][1]->nGrowVect()); + MultiFab::Copy(*Efield_aux[lev][2], *Efield_fp[lev][2], 0, 0, Efield_aux[lev][2]->nComp(), Efield_aux[lev][2]->nGrowVect()); } } } diff --git a/Source/Parallelization/WarpXComm_K.H b/Source/Parallelization/WarpXComm_K.H index a2b8fe38ed4..c3362087ad9 100644 --- a/Source/Parallelization/WarpXComm_K.H +++ b/Source/Parallelization/WarpXComm_K.H @@ -12,7 +12,7 @@ /** * \brief Interpolation function called within WarpX::UpdateAuxilaryDataSameType - * to interpolate data from the coarse and fine grids to the fine aux grid, + * with electromagnetic solver to interpolate data from the coarse and fine grids to the fine aux grid, * assuming that all grids have the same staggering (either collocated or staggered). * * \param[in] j index along x of the output array @@ -285,6 +285,94 @@ void warpx_interp (int j, int k, int l, // Final result arr_aux(j,k,l) = tmp + (fine - coarse); } +/** + * \brief Interpolation function called within WarpX::UpdateAuxilaryDataStagToNodal + * to interpolate data from the coarse and fine grids to the fine aux grid, + * with momentum-conserving field gathering, hence between grids with different staggering, + * and assuming that the aux grid is collocated. + * + * \param[in] j index along x of the output array + * \param[in] k index along y (in 3D) or z (in 2D) of the output array + * \param[in] l index along z (in 3D, l=0 in 2D) of the output array + * \param[in,out] arr_aux output array where interpolated values are stored + * \param[in] arr_fine input fine-patch array storing the values to interpolate + * \param[in] arr_fine_stag IndexType of the fine-patch arrays + */ +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void warpx_interp (int j, int k, int l, + amrex::Array4 const& arr_aux, + amrex::Array4 const& arr_fine, + const amrex::IntVect& arr_fine_stag) +{ + using namespace amrex; + + // Pad input arrays with zeros beyond ghost cells + // for out-of-bound accesses due to large-stencil operations + const auto arr_fine_zeropad = [arr_fine] (const int jj, const int kk, const int ll) noexcept + { + return arr_fine.contains(jj,kk,ll) ? arr_fine(jj,kk,ll) : 0.0_rt; + }; + + // NOTE Indices (j,k,l) in the following refer to: + // - (z,-,-) in 1D + // - (x,z,-) in 2D + // - (r,z,-) in RZ + // - (x,y,z) in 3D + + // Staggering of fine array (0: cell-centered; 1: nodal) + const int sj_fp = arr_fine_stag[0]; +#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + const int sk_fp = arr_fine_stag[1]; +#elif defined(WARPX_DIM_3D) + const int sk_fp = arr_fine_stag[1]; + const int sl_fp = arr_fine_stag[2]; +#endif + + // Number of points used for interpolation from coarse grid to fine grid + int nj; + int nk; + int nl; + + amrex::Real fine = 0.0_rt; + + // 3) Interpolation from fine staggered to fine nodal + + nj = (sj_fp == 0) ? 2 : 1; +#if defined(WARPX_DIM_1D_Z) + nk = 1; + nl = 1; +#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + nk = (sk_fp == 0) ? 2 : 1; + nl = 1; +#else + nk = (sk_fp == 0) ? 2 : 1; + nl = (sl_fp == 0) ? 2 : 1; +#endif + + const int jm = (sj_fp == 0) ? j-1 : j; +#if defined(WARPX_DIM_1D_Z) + const int km = k; + const int lm = l; +#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + const int km = (sk_fp == 0) ? k-1 : k; + const int lm = l; +#else + const int km = (sk_fp == 0) ? k-1 : k; + const int lm = (sl_fp == 0) ? l-1 : l; +#endif + + for (int jj = 0; jj < nj; jj++) { + for (int kk = 0; kk < nk; kk++) { + for (int ll = 0; ll < nl; ll++) { + fine += arr_fine_zeropad(jm+jj,km+kk,lm+ll); + } + } + } + fine = fine/static_cast(nj*nk*nl); + + // Final result + arr_aux(j,k,l) = fine; +} /** * \brief Arbitrary-order interpolation function used to center a given MultiFab between two grids diff --git a/Source/ablastr/fields/PoissonSolver.H b/Source/ablastr/fields/PoissonSolver.H index ca262981010..589c6ec1835 100644 --- a/Source/ablastr/fields/PoissonSolver.H +++ b/Source/ablastr/fields/PoissonSolver.H @@ -263,10 +263,15 @@ computePhi (amrex::Vector const & rho, mlmg.setVerbose(verbosity); mlmg.setMaxIter(max_iters); mlmg.setAlwaysUseBNorm(always_use_bnorm); + if (WarpX::grid_type == GridType::Collocated) { + // In this case, computeE needs to use ghost nodes data. So we + // ask MLMG to fill BC for us after it solves the problem. + mlmg.setFinalFillBC(true); + } // Solve Poisson equation at lev mlmg.solve( {phi[lev]}, {rho[lev]}, - relative_tolerance, absolute_tolerance ); + relative_tolerance, absolute_tolerance ); // needed for solving the levels by levels: // - coarser level is initial guess for finer level @@ -280,10 +285,14 @@ computePhi (amrex::Vector const & rho, const amrex::IntVect& refratio = rel_ref_ratio.value()[lev]; ba.coarsen(refratio); const int ncomp = linop.getNComp(); - amrex::MultiFab phi_cp(ba, phi[lev+1]->DistributionMap(), ncomp, 1); + const int ng = (WarpX::grid_type == GridType::Collocated) ? 1 : 0; + amrex::MultiFab phi_cp(ba, phi[lev+1]->DistributionMap(), ncomp, ng); + if (ng > 0) { + // Set all values outside the domain to zero + phi_cp.setDomainBndry(0.0_rt, geom[lev]); + } // Copy from phi[lev] to phi_cp (in parallel) - const amrex::IntVect& ng = amrex::IntVect::TheUnitVector(); const amrex::Periodicity& crse_period = geom[lev].periodicity(); ablastr::utils::communication::ParallelCopy( @@ -292,8 +301,8 @@ computePhi (amrex::Vector const & rho, 0, 0, 1, - ng, - ng, + amrex::IntVect(0), + amrex::IntVect(ng), do_single_precision_comms, crse_period ); @@ -308,7 +317,7 @@ computePhi (amrex::Vector const & rho, details::PoissonInterpCPtoFP const interp(phi_fp_arr, phi_cp_arr, refratio); - amrex::Box const b = mfi.tilebox(phi[lev + 1]->ixType().toIntVect()); + amrex::Box const& b = mfi.growntilebox(ng); amrex::ParallelFor(b, interp); } From c8a2479a8b4cb94e8be7e6d2c268ec9bd8581a0f Mon Sep 17 00:00:00 2001 From: Edoardo Zoni <59625522+EZoni@users.noreply.github.com> Date: Tue, 30 Jul 2024 06:57:53 -0700 Subject: [PATCH 12/12] Enable tiling in some `MFIter` loops (#5096) --- Source/Parallelization/WarpXComm.cpp | 4 ++-- Source/Utils/WarpXMovingWindow.cpp | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 2887bd4d056..6c44df061fd 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -306,7 +306,7 @@ WarpX::UpdateAuxilaryDataStagToNodal () #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif - for (MFIter mfi(*Efield_aux[lev][0]); mfi.isValid(); ++mfi) + for (MFIter mfi(*Efield_aux[lev][0], TilingIfNotGPU()); mfi.isValid(); ++mfi) { Array4 const& ex_aux = Efield_aux[lev][0]->array(mfi); Array4 const& ey_aux = Efield_aux[lev][1]->array(mfi); @@ -321,7 +321,7 @@ WarpX::UpdateAuxilaryDataStagToNodal () Array4 const& ey_c = Etmp[1]->const_array(mfi); Array4 const& ez_c = Etmp[2]->const_array(mfi); - const Box& bx = mfi.fabbox(); + const Box& bx = mfi.growntilebox(); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept { diff --git a/Source/Utils/WarpXMovingWindow.cpp b/Source/Utils/WarpXMovingWindow.cpp index d64394dafc9..73696838cd4 100644 --- a/Source/Utils/WarpXMovingWindow.cpp +++ b/Source/Utils/WarpXMovingWindow.cpp @@ -534,7 +534,7 @@ WarpX::shiftMF (amrex::MultiFab& mf, const amrex::Geometry& geom, #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif - for (amrex::MFIter mfi(tmpmf); mfi.isValid(); ++mfi ) + for (amrex::MFIter mfi(tmpmf, TilingIfNotGPU()); mfi.isValid(); ++mfi ) { if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers) { @@ -545,7 +545,7 @@ WarpX::shiftMF (amrex::MultiFab& mf, const amrex::Geometry& geom, auto const& dstfab = mf.array(mfi); auto const& srcfab = tmpmf.array(mfi); - const amrex::Box& outbox = mfi.fabbox() & adjBox; + const amrex::Box& outbox = mfi.growntilebox() & adjBox; if (outbox.ok()) { if (!useparser) {