From a689ae9e5ff46081db05c9609730df00b90c90eb Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Tue, 12 Mar 2024 10:04:56 +0900 Subject: [PATCH 01/14] AMReX/pyAMReX/PICSAR: Weekly Update (#4763) * 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 5e3a3077f5d..1c384183fbd 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 24.03 && cd - + cd ../amrex && git checkout --detach 155f1f4c3f4fa03528d38d2769beb4b35b38042c && cd - make COMP=gcc QED=FALSE USE_MPI=TRUE USE_GPU=TRUE USE_OMP=FALSE USE_PSATD=TRUE USE_CCACHE=TRUE -j 4 ccache -s diff --git a/Regression/WarpX-GPU-tests.ini b/Regression/WarpX-GPU-tests.ini index 2f0c7b607bb..138cf90534f 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 = 24.03 +branch = 155f1f4c3f4fa03528d38d2769beb4b35b38042c [source] dir = /home/regtester/git/WarpX diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index d6ac81c5607..daeef409dc8 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 = 24.03 +branch = 155f1f4c3f4fa03528d38d2769beb4b35b38042c [source] dir = /home/regtester/AMReX_RegTesting/warpx diff --git a/cmake/dependencies/AMReX.cmake b/cmake/dependencies/AMReX.cmake index a1045457ec0..54b1cc45322 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 "24.03" +set(WarpX_amrex_branch "155f1f4c3f4fa03528d38d2769beb4b35b38042c" 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 6360da84276..eaa0fc016cd 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 "24.03" +set(WarpX_pyamrex_branch "7cd4a4220bc3fc523cd0d78f7320909aa2feb47a" CACHE STRING "Repository branch for WarpX_pyamrex_repo if(WarpX_pyamrex_internal)") diff --git a/run_test.sh b/run_test.sh index f556dd37e76..1bf6e96d236 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 24.03 && cd - +cd amrex && git checkout --detach 155f1f4c3f4fa03528d38d2769beb4b35b38042c && 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 7317aafb9d99b74d077c5bcc1bae3af29f31d9fd Mon Sep 17 00:00:00 2001 From: Arianna Formenti Date: Tue, 12 Mar 2024 06:54:21 -0700 Subject: [PATCH 02/14] clean up (#4761) --- Source/Particles/PhysicalParticleContainer.H | 8 +--- .../Particles/PhysicalParticleContainer.cpp | 40 ++++++++----------- 2 files changed, 18 insertions(+), 30 deletions(-) diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index b77c1147a15..286d0675da6 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -226,13 +226,7 @@ public: amrex::ParticleReal& ux, amrex::ParticleReal& uy, amrex::ParticleReal& uz, amrex::Real t_lab = 0.) const; - void AddGaussianBeam ( - PlasmaInjector const& plasma_injector, - amrex::Real x_m, amrex::Real y_m, amrex::Real z_m, - amrex::Real x_rms, amrex::Real y_rms, amrex::Real z_rms, - amrex::Real x_cut, amrex::Real y_cut, amrex::Real z_cut, - amrex::Real q_tot, long npart, int do_symmetrize, int symmetrization_order, - amrex::Real focal_distance); + void AddGaussianBeam (PlasmaInjector const& plasma_injector); /** Load a particle beam from an external file * @param[in] the PlasmaInjector instance holding the input parameters diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 1aa534c871b..160eac0d19c 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -521,14 +521,22 @@ void PhysicalParticleContainer::MapParticletoBoostedFrame ( } void -PhysicalParticleContainer::AddGaussianBeam ( - PlasmaInjector const& plasma_injector, - const Real x_m, const Real y_m, const Real z_m, - const Real x_rms, const Real y_rms, const Real z_rms, - const Real x_cut, const Real y_cut, const Real z_cut, - const Real q_tot, long npart, - const int do_symmetrize, - const int symmetrization_order, const Real focal_distance) { +PhysicalParticleContainer::AddGaussianBeam (PlasmaInjector const& plasma_injector){ + + const Real x_m = plasma_injector.x_m; + const Real y_m = plasma_injector.y_m; + const Real z_m = plasma_injector.z_m; + const Real x_rms = plasma_injector.x_rms; + const Real y_rms = plasma_injector.y_rms; + const Real z_rms = plasma_injector.z_rms; + const Real x_cut = plasma_injector.x_cut; + const Real y_cut = plasma_injector.y_cut; + const Real z_cut = plasma_injector.z_cut; + const Real q_tot = plasma_injector.q_tot; + long npart = plasma_injector.npart; + const int do_symmetrize = plasma_injector.do_symmetrize; + const int symmetrization_order = plasma_injector.symmetrization_order; + const Real focal_distance = plasma_injector.focal_distance; // Declare temporary vectors on the CPU Gpu::HostVector particle_x; @@ -918,21 +926,7 @@ PhysicalParticleContainer::AddParticles (int lev) } if (plasma_injector->gaussian_beam) { - AddGaussianBeam(*plasma_injector, - plasma_injector->x_m, - plasma_injector->y_m, - plasma_injector->z_m, - plasma_injector->x_rms, - plasma_injector->y_rms, - plasma_injector->z_rms, - plasma_injector->x_cut, - plasma_injector->y_cut, - plasma_injector->z_cut, - plasma_injector->q_tot, - plasma_injector->npart, - plasma_injector->do_symmetrize, - plasma_injector->symmetrization_order, - plasma_injector->focal_distance); + AddGaussianBeam(*plasma_injector); } if (plasma_injector->external_file) { From 966efe04ccfab2ef3369af86a315d5b90b43dec3 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Wed, 13 Mar 2024 17:21:41 +0100 Subject: [PATCH 03/14] remove superfluous include directives (#4770) --- Source/Initialization/WarpXAMReXInit.cpp | 2 -- Source/ablastr/utils/text/StreamUtils.cpp | 1 - 2 files changed, 3 deletions(-) diff --git a/Source/Initialization/WarpXAMReXInit.cpp b/Source/Initialization/WarpXAMReXInit.cpp index d1dd0bbe90b..8cfe78c740f 100644 --- a/Source/Initialization/WarpXAMReXInit.cpp +++ b/Source/Initialization/WarpXAMReXInit.cpp @@ -11,8 +11,6 @@ #include #include -#include - namespace { /** Overwrite defaults in AMReX Inputs * diff --git a/Source/ablastr/utils/text/StreamUtils.cpp b/Source/ablastr/utils/text/StreamUtils.cpp index 965beb179ba..541ef300f7c 100644 --- a/Source/ablastr/utils/text/StreamUtils.cpp +++ b/Source/ablastr/utils/text/StreamUtils.cpp @@ -7,7 +7,6 @@ #include "StreamUtils.H" -#include #include void From 41983ae761c0da0dc6220c25b547f35d4ce3bb9f Mon Sep 17 00:00:00 2001 From: "S. Eric Clark" <25495882+clarkse@users.noreply.github.com> Date: Wed, 13 Mar 2024 21:48:53 -0700 Subject: [PATCH 04/14] Adding hyper-resistivity to generalized ohms law hybrid solver. (#4772) * Adding hyper-resistivity to generalized ohms law hybrid solver. * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Apply suggestions from code review Style changes to clean up Laplacian operators. Co-authored-by: Roelof Groenewald <40245517+roelof-groenewald@users.noreply.github.com> * Removed extra arguments for triggering calculation of hyper_resistivity. --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Roelof Groenewald <40245517+roelof-groenewald@users.noreply.github.com> --- Python/pywarpx/picmi.py | 5 +- .../CartesianYeeAlgorithm.H | 58 +++++++++++++++++++ .../CylindricalYeeAlgorithm.H | 33 +++++++++++ .../HybridPICModel/HybridPICModel.H | 3 + .../HybridPICModel/HybridPICModel.cpp | 2 + .../HybridPICSolveE.cpp | 50 +++++++++++++--- .../FieldSolver/WarpXPushFieldsHybridPIC.cpp | 3 +- 7 files changed, 144 insertions(+), 10 deletions(-) diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index 959db6c701d..59ca1ec789f 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -1142,7 +1142,8 @@ class HybridPICSolver(picmistandard.base._ClassWithInit): Function of space and time specifying external (non-plasma) currents. """ def __init__(self, grid, Te=None, n0=None, gamma=None, - n_floor=None, plasma_resistivity=None, substeps=None, + n_floor=None, plasma_resistivity=None, + plasma_hyper_resistivity=None, substeps=None, Jx_external_function=None, Jy_external_function=None, Jz_external_function=None, **kw): self.grid = grid @@ -1153,6 +1154,7 @@ def __init__(self, grid, Te=None, n0=None, gamma=None, self.gamma = gamma self.n_floor = n_floor self.plasma_resistivity = plasma_resistivity + self.plasma_hyper_resistivity = plasma_hyper_resistivity self.substeps = substeps @@ -1187,6 +1189,7 @@ def solver_initialize_inputs(self): 'plasma_resistivity(rho,J)', pywarpx.my_constants.mangle_expression(self.plasma_resistivity, self.mangle_dict) ) + pywarpx.hybridpicmodel.plasma_hyper_resistivity = self.plasma_hyper_resistivity pywarpx.hybridpicmodel.substeps = self.substeps pywarpx.hybridpicmodel.__setattr__( 'Jx_external_grid_function(x,y,z,t)', diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H index b762530e1a3..c4978287aec 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H @@ -100,6 +100,25 @@ struct CartesianYeeAlgorithm { #endif } + /** + * Perform second derivative along x on a cell-centered grid, from a cell-centered field `F`*/ + template< typename T_Field> + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real Dxx ( + T_Field const& F, + amrex::Real const * const coefs_x, int const /*n_coefs_x*/, + int const i, int const j, int const k, int const ncomp=0 ) { + + using namespace amrex; +#if (defined WARPX_DIM_1D_Z) + amrex::ignore_unused(F, coefs_x, i, j, k, ncomp); + return 0._rt; // 1D Cartesian: derivative along x is 0 +#else + amrex::Real const inv_dx2 = coefs_x[0]*coefs_x[0]; + return inv_dx2*( F(i-1,j,k,ncomp) - 2._rt*F(i,j,k,ncomp) + F(i+1,j,k,ncomp) ); +#endif + } + /** * Perform derivative along y on a cell-centered grid, from a nodal field `F`*/ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE @@ -147,6 +166,25 @@ struct CartesianYeeAlgorithm { #endif } + /** + * Perform derivative along y on a nodal grid, from a cell-centered field `F`*/ + template< typename T_Field> + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real Dyy ( + T_Field const& F, + amrex::Real const * const coefs_y, int const /*n_coefs_y*/, + int const i, int const j, int const k, int const ncomp=0 ) { + + using namespace amrex; +#if defined WARPX_DIM_3D + Real const inv_dy2 = coefs_y[0]*coefs_y[0]; + return inv_dy2*( F(i,j-1,k,ncomp) - 2._rt*F(i,j,k,ncomp) + F(i,j+1,k,ncomp) ); +#elif (defined WARPX_DIM_XZ || WARPX_DIM_1D_Z) + amrex::ignore_unused(F, coefs_y, i, j, k, ncomp); + return 0._rt; // 1D and 2D Cartesian: derivative along y is 0 +#endif + } + /** * Perform derivative along z on a cell-centered grid, from a nodal field `F`*/ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE @@ -186,6 +224,26 @@ struct CartesianYeeAlgorithm { #endif } + /** + * Perform second derivative along z on a cell-centered field `F`*/ + template< typename T_Field> + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real Dzz ( + T_Field const& F, + amrex::Real const * const coefs_z, int const /*n_coefs_z*/, + int const i, int const j, int const k, int const ncomp=0 ) { + + using namespace amrex; + Real const inv_dz2 = coefs_z[0]*coefs_z[0]; +#if defined WARPX_DIM_3D + return inv_dz2*( F(i,j,k-1,ncomp) - 2._rt*F(i,j,k,ncomp) + F(i,j,k+1,ncomp) ); +#elif (defined WARPX_DIM_XZ) + return inv_dz2*( F(i,j-1,k,ncomp) - 2._rt*F(i,j,k,ncomp) + F(i,j+1,k,ncomp) ); +#elif (defined WARPX_DIM_1D_Z) + return inv_dz2*( F(i-1,j,k,ncomp) - 2._rt*F(i,j,k,ncomp) + F(i+1,j,k,ncomp) ); +#endif + } + }; #endif // WARPX_FINITE_DIFFERENCE_ALGORITHM_CARTESIAN_YEE_H_ diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H index e49995557a0..436f0a83f31 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H @@ -134,6 +134,24 @@ struct CylindricalYeeAlgorithm { return inv_dr*( F(i,j,k,comp) - F(i-1,j,k,comp) ); } + /** Applies the differential operator `1/r * d(r * dF/dr)/dr`, + * where `F` is on a *cell-centered* or a nodal grid in `r` + * The input parameter `r` is given at the cell-centered position */ + template< typename T_Field> + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real Dr_rDr_over_r ( + T_Field const& F, + amrex::Real const r, amrex::Real const dr, + amrex::Real const * const coefs_r, int const /*n_coefs_r*/, + int const i, int const j, int const k, int const comp ) { + + using namespace amrex; + + Real const inv_dr2 = coefs_r[0]*coefs_r[0]; + return 1._rt/r * inv_dr2*( (r+0.5_rt*dr)*(F(i+1,j,k,comp) - F(i,j,k,comp)) + - (r-0.5_rt*dr)*(F(i,j,k,comp) - F(i-1,j,k,comp)) ); + } + /** * Perform derivative along z on a cell-centered grid, from a nodal field `F` */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE @@ -162,6 +180,21 @@ struct CylindricalYeeAlgorithm { return inv_dz*( F(i,j,k,comp) - F(i,j-1,k,comp) ); } + /** + * Perform second derivative along z on a cell-centered field `F`*/ + template< typename T_Field> + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real Dzz ( + T_Field const& F, + amrex::Real const * const coefs_z, int const /*n_coefs_z*/, + int const i, int const j, int const k, int const ncomp=0 ) { + + using namespace amrex; + Real const inv_dz2 = coefs_z[0]*coefs_z[0]; + + return inv_dz2*( F(i,j-1,k,ncomp) - 2._rt*F(i,j,k,ncomp) + F(i,j+1,k,ncomp) ); + } + }; #endif // WARPX_FINITE_DIFFERENCE_ALGORITHM_CYLINDRICAL_YEE_H_ diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H index 23ef49b58cb..5e4c6382bfc 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H @@ -175,6 +175,9 @@ public: amrex::ParserExecutor<2> m_eta; bool m_resistivity_has_J_dependence = false; + /** Plasma hyper-resisitivity */ + amrex::Real m_eta_h = 0.0; + /** External current */ std::string m_Jx_ext_grid_function = "0.0"; std::string m_Jy_ext_grid_function = "0.0"; diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp index 8c2b5a0707a..e8330242fe8 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp @@ -40,6 +40,8 @@ void HybridPICModel::ReadParameters () pp_hybrid.query("plasma_resistivity(rho,J)", m_eta_expression); utils::parser::queryWithParser(pp_hybrid, "n_floor", m_n_floor); + utils::parser::queryWithParser(pp_hybrid, "plasma_hyper_resistivity", m_eta_h); + // convert electron temperature from eV to J m_elec_temp *= PhysConst::q_e; diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp index 066e88b13c0..456c542a534 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp @@ -375,7 +375,7 @@ void FiniteDifferenceSolver::HybridPICSolveE ( std::unique_ptr const& Pefield, std::array< std::unique_ptr, 3 > const& edge_lengths, int lev, HybridPICModel const* hybrid_model, - const bool include_resistivity_term ) + const bool include_resistivity_term) { // Select algorithm (The choice of algorithm is a runtime option, // but we compile code for each algorithm, using templates) @@ -432,9 +432,12 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical ( // get hybrid model parameters const auto eta = hybrid_model->m_eta; + const auto eta_h = hybrid_model->m_eta_h; const auto rho_floor = hybrid_model->m_n_floor * PhysConst::q_e; const auto resistivity_has_J_dependence = hybrid_model->m_resistivity_has_J_dependence; + const bool include_hyper_resistivity_term = (eta_h > 0.0) && include_resistivity_term; + // Index type required for interpolating fields from their respective // staggering to the Ex, Ey, Ez locations amrex::GpuArray const& Er_stag = hybrid_model->Ex_IndexType; @@ -612,6 +615,14 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical ( // Add resistivity only if E field value is used to update B if (include_resistivity_term) { Er(i, j, 0) += eta(rho_val, jtot_val) * Jr(i, j, 0); } + + if (include_hyper_resistivity_term) { + // r on cell-centered point (Jr is cell-centered in r) + Real const r = rmin + (i + 0.5_rt)*dr; + + auto nabla2Jr = T_Algo::Dr_rDr_over_r(Jr, r, dr, coefs_r, n_coefs_r, i, j, 0, 0); + Er(i, j, 0) -= eta_h * nabla2Jr; + } }, // Et calculation @@ -655,16 +666,18 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical ( // Add resistivity only if E field value is used to update B if (include_resistivity_term) { Et(i, j, 0) += eta(rho_val, jtot_val) * Jt(i, j, 0); } + + // Note: Hyper-resisitivity should be revisited here when modal decomposition is implemented }, // Ez calculation - [=] AMREX_GPU_DEVICE (int i, int j, int k){ + [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/){ #ifdef AMREX_USE_EB // Skip field solve if this cell is fully covered by embedded boundaries if (lz(i,j,0) <= 0) { return; } #endif // Interpolate to get the appropriate charge density in space - Real rho_val = Interp(rho, nodal, Ez_stag, coarsen, i, j, k, 0); + Real rho_val = Interp(rho, nodal, Ez_stag, coarsen, i, j, 0, 0); // Interpolate current to appropriate staggering to match E field Real jtot_val = 0._rt; @@ -679,15 +692,20 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical ( if (rho_val < rho_floor) { rho_val = rho_floor; } // Get the gradient of the electron pressure - auto grad_Pe = T_Algo::UpwardDz(Pe, coefs_z, n_coefs_z, i, j, k, 0); + auto grad_Pe = T_Algo::UpwardDz(Pe, coefs_z, n_coefs_z, i, j, 0, 0); // interpolate the nodal neE values to the Yee grid - auto enE_z = Interp(enE, nodal, Ez_stag, coarsen, i, j, k, 2); + auto enE_z = Interp(enE, nodal, Ez_stag, coarsen, i, j, 0, 2); - Ez(i, j, k) = (enE_z - grad_Pe) / rho_val; + Ez(i, j, 0) = (enE_z - grad_Pe) / rho_val; // Add resistivity only if E field value is used to update B - if (include_resistivity_term) { Ez(i, j, k) += eta(rho_val, jtot_val) * Jz(i, j, k); } + if (include_resistivity_term) { Ez(i, j, 0) += eta(rho_val, jtot_val) * Jz(i, j, 0); } + + if (include_hyper_resistivity_term) { + auto nabla2Jz = T_Algo::Dzz(Jz, coefs_z, n_coefs_z, i, j, 0, 0); + Ez(i, j, 0) -= eta_h * nabla2Jz; + } } ); @@ -726,9 +744,12 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian ( // get hybrid model parameters const auto eta = hybrid_model->m_eta; + const auto eta_h = hybrid_model->m_eta_h; const auto rho_floor = hybrid_model->m_n_floor * PhysConst::q_e; const auto resistivity_has_J_dependence = hybrid_model->m_resistivity_has_J_dependence; + const bool include_hyper_resistivity_term = (eta_h > 0.) && include_resistivity_term; + // Index type required for interpolating fields from their respective // staggering to the Ex, Ey, Ez locations amrex::GpuArray const& Ex_stag = hybrid_model->Ex_IndexType; @@ -904,6 +925,11 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian ( // Add resistivity only if E field value is used to update B if (include_resistivity_term) { Ex(i, j, k) += eta(rho_val, jtot_val) * Jx(i, j, k); } + + if (include_hyper_resistivity_term) { + auto nabla2Jx = T_Algo::Dxx(Jx, coefs_x, n_coefs_x, i, j, k); + Ex(i, j, k) -= eta_h * nabla2Jx; + } }, // Ey calculation @@ -943,6 +969,11 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian ( // Add resistivity only if E field value is used to update B if (include_resistivity_term) { Ey(i, j, k) += eta(rho_val, jtot_val) * Jy(i, j, k); } + + if (include_hyper_resistivity_term) { + auto nabla2Jy = T_Algo::Dyy(Jy, coefs_y, n_coefs_y, i, j, k); + Ey(i, j, k) -= eta_h * nabla2Jy; + } }, // Ez calculation @@ -976,6 +1007,11 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian ( // Add resistivity only if E field value is used to update B if (include_resistivity_term) { Ez(i, j, k) += eta(rho_val, jtot_val) * Jz(i, j, k); } + + if (include_hyper_resistivity_term) { + auto nabla2Jz = T_Algo::Dzz(Jz, coefs_z, n_coefs_z, i, j, k); + Ez(i, j, k) -= eta_h * nabla2Jz; + } } ); diff --git a/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp b/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp index ae10ae1e19a..4dbe10c4e5a 100644 --- a/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp +++ b/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp @@ -154,8 +154,7 @@ void WarpX::HybridPICEvolveFields () // Update the E field to t=n+1 using the extrapolated J_i^n+1 value m_hybrid_pic_model->CalculateCurrentAmpere(Bfield_fp, m_edge_lengths); m_hybrid_pic_model->HybridPICSolveE( - Efield_fp, current_fp_temp, Bfield_fp, rho_fp, m_edge_lengths, - false + Efield_fp, current_fp_temp, Bfield_fp, rho_fp, m_edge_lengths, false ); FillBoundaryE(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points); From 2630be83f210d32fed4a48bbb28ec3db372cf433 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Thu, 14 Mar 2024 21:43:05 +0900 Subject: [PATCH 05/14] Doc: Conda `-y` Install Block (#4776) * Doc: Update libmamba Usage in Conda * Fix: The configuration key changed slightly. * Cosmetic: Adding `-y` ensures one can execute the whole block at once. * `solver` is the new (conda>=23) key! * Fix: Double backticks for verbatim --- Docs/source/install/dependencies.rst | 4 ++-- Docs/source/install/users.rst | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/Docs/source/install/dependencies.rst b/Docs/source/install/dependencies.rst index 94b2be78bec..4f3544d5c0d 100644 --- a/Docs/source/install/dependencies.rst +++ b/Docs/source/install/dependencies.rst @@ -63,8 +63,8 @@ Conda (Linux/macOS/Windows) .. code-block:: bash - conda update -n base conda - conda install -n base conda-libmamba-solver + conda update -y -n base conda + conda install -y -n base conda-libmamba-solver conda config --set solver libmamba We recommend to deactivate that conda self-activates its ``base`` environment. diff --git a/Docs/source/install/users.rst b/Docs/source/install/users.rst index b7893e248c0..e56d2d8ac43 100644 --- a/Docs/source/install/users.rst +++ b/Docs/source/install/users.rst @@ -45,12 +45,12 @@ A package for WarpX is available via the `Conda `_ package man .. tip:: - We recommend to configure your conda to use the faster `libmamba` `dependency solver `__. + We recommend to configure your conda to use the faster ``libmamba`` `dependency solver `__. .. code-block:: bash - conda update -n base conda - conda install -n base conda-libmamba-solver + conda update -y -n base conda + conda install -y -n base conda-libmamba-solver conda config --set solver libmamba We recommend to deactivate that conda self-activates its ``base`` environment. From 467598cb56a1bdc886c8c2134fe9b33e748f9de6 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Fri, 15 Mar 2024 00:50:28 +0900 Subject: [PATCH 06/14] Zenodo: Add Justin & Kale as Co-Authors (#4760) * Zenodo: Add Justin & Kale as Co-Authors Adding Justin and Kale as a co-author of WarpX. Thank you for your contributions! :sparkle: * Update Order in Kale's Affiliation * Update Justin MN --- .zenodo.json | 22 ++++++++++++++++------ 1 file changed, 16 insertions(+), 6 deletions(-) diff --git a/.zenodo.json b/.zenodo.json index bbba267def9..5431e40c8bb 100644 --- a/.zenodo.json +++ b/.zenodo.json @@ -20,6 +20,11 @@ "name": "Andriyash, Igor", "orcid": "0000-0003-0313-4496" }, + { + "affiliation": "Lawrence Livermore National Laboratory", + "name": "Angus, Justin Ray", + "orcid": "0000-0003-1474-0002" + }, { "affiliation": "Lawrence Berkeley National Laboratory", "name": "Belkin, Daniel", @@ -155,7 +160,7 @@ }, { "affiliation": "Lawrence Berkeley National Laboratory", - "name": "Sandberg, Ryan T.", + "name": "Sandberg, Ryan Thor", "orcid": "0000-0001-7680-8733" }, { @@ -163,11 +168,21 @@ "name": "Scherpelz, Peter", "orcid": "0000-0001-8185-3387" }, + { + "affiliation": "Laboratory for Laser Energetics, University of Rochester", + "name": "Weichman, Kale", + "orcid": "0000-0002-3487-7922" + }, { "affiliation": "Lawrence Berkeley National Laboratory", "name": "Yang, Eloise", "orcid": "0000-0002-9319-4216" }, + { + "affiliation": "LIDYL, CEA-Universit\u00e9 Paris-Saclay, CEA Saclay", + "name": "Zaim, Ne\u00efl", + "orcid": "0000-0003-0313-4496" + }, { "affiliation": "Lawrence Berkeley National Laboratory", "name": "Zhang, Weiqun", @@ -187,11 +202,6 @@ "affiliation": "Lawrence Berkeley National Laboratory", "name": "Zoni, Edoardo", "orcid": "0000-0001-5662-4646" - }, - { - "affiliation": "LIDYL, CEA-Universit\u00e9 Paris-Saclay, CEA Saclay", - "name": "Zaim, Ne\u00efl", - "orcid": "0000-0003-0313-4496" } ], "contributors": [ From 5f96976b18f19fe748e021a3f4ecf6999c68c038 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Thu, 14 Mar 2024 11:33:36 -0700 Subject: [PATCH 07/14] Compute electrostatic fields at the beginning of EM simulations, if potential is specified (#4723) * Compute electrostatic fields at the beginning of EM simulations, if phi is specified * Fix tests * Add documentation * Update test to include initial electrostatic field * Avoid issue in CI * Set the potential to be uniform * Update checksum * Update automated test * Update checksum * Add warning message * Update Source/WarpX.cpp * Update Source/WarpX.cpp * Update Source/WarpX.cpp Co-authored-by: Roelof Groenewald <40245517+roelof-groenewald@users.noreply.github.com> * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --------- Co-authored-by: Roelof Groenewald <40245517+roelof-groenewald@users.noreply.github.com> Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- Docs/source/usage/parameters.rst | 17 ++++++--- .../Tests/embedded_boundary_cube/inputs_2d | 6 +++- .../Tests/embedded_boundary_cube/inputs_3d | 5 +++ Examples/Tests/point_of_contact_EB/inputs_3d | 6 ---- Examples/Tests/point_of_contact_EB/inputs_rz | 6 ---- .../embedded_boundary_cube.json | 15 ++++---- .../embedded_boundary_cube_2d.json | 10 +++--- .../embedded_boundary_cube_macroscopic.json | 10 +++--- Regression/WarpX-tests.ini | 6 ++-- Source/FieldSolver/ElectrostaticSolver.cpp | 3 +- Source/WarpX.H | 1 + Source/WarpX.cpp | 35 +++++++++++++++---- 12 files changed, 74 insertions(+), 46 deletions(-) diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index 0c5fe85973f..a615c9c8a87 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -395,12 +395,18 @@ Domain Boundary Conditions * ``damped``: This is the recommended option in the moving direction when using the spectral solver with moving window (currently only supported along z). This boundary condition applies a damping factor to the electric and magnetic fields in the outer half of the guard cells, using a sine squared profile. As the spectral solver is by nature periodic, the damping prevents fields from wrapping around to the other end of the domain when the periodicity is not desired. This boundary condition is only valid when using the spectral solver. * ``pec``: This option can be used to set a Perfect Electric Conductor at the simulation boundary. Please see the :ref:`PEC theory section ` for more details. Note that PEC boundary is invalid at `r=0` for the RZ solver. Please use ``none`` option. This boundary condition does not work with the spectral solver. - If an electrostatic field solve is used the boundary potentials can also be set through ``boundary.potential_lo_x/y/z`` and ``boundary.potential_hi_x/y/z`` (default `0`). * ``none``: No boundary condition is applied to the fields with the electromagnetic solver. This option must be used for the RZ-solver at `r=0`. * ``neumann``: For the electrostatic solver, a Neumann boundary condition (with gradient of the potential equal to 0) will be applied on the specified boundary. +* ``boundary.potential_lo_x/y/z`` and ``boundary.potential_hi_x/y/z`` (default `0`) + Gives the value of the electric potential at the boundaries, for ``pec`` boundaries. With electrostatic solvers + (i.e., with ``warpx.do_electrostatic = ...``), this is used in order to compute the potential + in the simulation volume at each timestep. When using other solvers (e.g. Maxwell solver), + setting these variables will trigger an electrostatic solve at ``t=0``, to compute the initial + electric field produced by the boundaries. + * ``boundary.particle_lo`` and ``boundary.particle_hi`` (`2 strings` for 2D, `3 strings` for 3D, `absorbing` by default) Options are: @@ -473,9 +479,12 @@ Embedded Boundary Conditions the interior of the embeddded boundary is where the function value is positive. * ``warpx.eb_potential(x,y,z,t)`` (`string`) - Only used when ``warpx.do_electrostatic=labframe``. Gives the value of - the electric potential at the surface of the embedded boundary, - as a function of `x`, `y`, `z` and time. This function is also evaluated + Gives the value of the electric potential at the surface of the embedded boundary, + as a function of `x`, `y`, `z` and `t`. With electrostatic solvers (i.e., with + ``warpx.do_electrostatic = ...``), this is used in order to compute the potential + in the simulation volume at each timestep. When using other solvers (e.g. Maxwell solver), + setting this variable will trigger an electrostatic solve at ``t=0``, to compute the initial + electric field produced by the boundaries. Note that this function is also evaluated inside the embedded boundary. For this reason, it is important to define this function in such a way that it is constant inside the embedded boundary. diff --git a/Examples/Tests/embedded_boundary_cube/inputs_2d b/Examples/Tests/embedded_boundary_cube/inputs_2d index 476992b360d..372e0dc0340 100644 --- a/Examples/Tests/embedded_boundary_cube/inputs_2d +++ b/Examples/Tests/embedded_boundary_cube/inputs_2d @@ -15,10 +15,14 @@ my_constants.xmin = -0.5 my_constants.zmin = -0.5 my_constants.xmax = 0.5 my_constants.zmax = 0.5 -# Alternatively one could use parser to build EB # Note that for amrex EB implicit function, >0 is covered, =0 is boundary and <0 is regular. warpx.eb_implicit_function = "max(max(x+xmin,-(x+xmax)), max(z+zmin,-(z+zmax)))" +# To test the initial electrostatic solver, we also add uniform potential +# Since, the potential is uniform here, it should not modify the fields, so this mainly +# tests that calling the electrostatic solver does not raise any error/crash +warpx.eb_potential(x,y,z,t) = "1" + warpx.B_ext_grid_init_style = parse_B_ext_grid_function warpx.E_ext_grid_init_style = parse_E_ext_grid_function diff --git a/Examples/Tests/embedded_boundary_cube/inputs_3d b/Examples/Tests/embedded_boundary_cube/inputs_3d index a98cd0b75ac..61eb1192e04 100644 --- a/Examples/Tests/embedded_boundary_cube/inputs_3d +++ b/Examples/Tests/embedded_boundary_cube/inputs_3d @@ -19,6 +19,11 @@ eb2.box_has_fluid_inside = true # Note that for amrex EB implicit function, >0 is covered, =0 is boundary and <0 is regular. # warpx.eb_implicit_function = "max(max(max(x-0.5,-0.5-x), max(y-0.5,-0.5-y)), max(z-0.5,-0.5-z))" +# To test the initial electrostatic solver, we also add uniform potential +# Since, the potential is uniform here, it should not modify the fields, so this mainly +# tests that calling the electrostatic solver does not raise any error/crash +warpx.eb_potential(x,y,z,t) = "1" + warpx.B_ext_grid_init_style = parse_B_ext_grid_function my_constants.m = 0 my_constants.n = 1 diff --git a/Examples/Tests/point_of_contact_EB/inputs_3d b/Examples/Tests/point_of_contact_EB/inputs_3d index 005733487e5..8b190c056ff 100644 --- a/Examples/Tests/point_of_contact_EB/inputs_3d +++ b/Examples/Tests/point_of_contact_EB/inputs_3d @@ -12,12 +12,6 @@ geometry.prob_hi = 0.26 0.26 0.26 boundary.field_lo = pec pec pec boundary.field_hi = pec pec pec -boundary.potential_lo_x = 0 -boundary.potential_hi_x = 0 -boundary.potential_lo_y = 0 -boundary.potential_hi_y = 0 -boundary.potential_lo_z = 0 -boundary.potential_hi_z = 0 warpx.const_dt = 1e-10 warpx.eb_implicit_function = "-(x**2+y**2+z**2-0.2**2)" diff --git a/Examples/Tests/point_of_contact_EB/inputs_rz b/Examples/Tests/point_of_contact_EB/inputs_rz index a492c1f705a..169a58adb94 100644 --- a/Examples/Tests/point_of_contact_EB/inputs_rz +++ b/Examples/Tests/point_of_contact_EB/inputs_rz @@ -12,12 +12,6 @@ geometry.prob_hi = 0.26 0.26 boundary.field_lo = none periodic boundary.field_hi = pec periodic -boundary.potential_lo_x = 0 -boundary.potential_hi_x = 0 -boundary.potential_lo_y = 0 -boundary.potential_hi_y = 0 -boundary.potential_lo_z = 0 -boundary.potential_hi_z = 0 warpx.const_dt = 1e-10 warpx.eb_implicit_function = "-(x**2 -0.2**2)" diff --git a/Regression/Checksum/benchmarks_json/embedded_boundary_cube.json b/Regression/Checksum/benchmarks_json/embedded_boundary_cube.json index 69fffed2c6a..d0edafb9f0e 100644 --- a/Regression/Checksum/benchmarks_json/embedded_boundary_cube.json +++ b/Regression/Checksum/benchmarks_json/embedded_boundary_cube.json @@ -1,11 +1,10 @@ { "lev=0": { - "Bx": 0.0, - "By": 0.0066283741197868335, - "Bz": 0.006628374119786833, - "Ex": 5102618.47115243, - "Ey": 0.0, - "Ez": 0.0 + "Bx": 3.769898030127477e-18, + "By": 0.006628374119786834, + "Bz": 0.006628374119786834, + "Ex": 5102618.4711524295, + "Ey": 6.323755130400527e-05, + "Ez": 6.323755130400527e-05 } -} - +} \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/embedded_boundary_cube_2d.json b/Regression/Checksum/benchmarks_json/embedded_boundary_cube_2d.json index 37a82967a23..a3e609bd9a9 100644 --- a/Regression/Checksum/benchmarks_json/embedded_boundary_cube_2d.json +++ b/Regression/Checksum/benchmarks_json/embedded_boundary_cube_2d.json @@ -1,10 +1,10 @@ { "lev=0": { - "Bx": 9.263694545408502e-05, - "By": 0.00031905198933489135, - "Bz": 7.328424783762596e-05, - "Ex": 8553.90669805307, + "Bx": 9.263694545408503e-05, + "By": 0.00031905198933489145, + "Bz": 7.328424783762594e-05, + "Ex": 8553.906698053046, "Ey": 60867.04830538045, - "Ez": 0.0 + "Ez": 8.439422682267567e-07 } } \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/embedded_boundary_cube_macroscopic.json b/Regression/Checksum/benchmarks_json/embedded_boundary_cube_macroscopic.json index 453b83f5e31..c759f36ac5d 100644 --- a/Regression/Checksum/benchmarks_json/embedded_boundary_cube_macroscopic.json +++ b/Regression/Checksum/benchmarks_json/embedded_boundary_cube_macroscopic.json @@ -1,10 +1,10 @@ { "lev=0": { - "Bx": 0.0, - "By": 0.005101824310293575, + "Bx": 3.898540288712651e-18, + "By": 0.005101824310293573, "Bz": 0.005101824310293573, - "Ex": 4414725.184731115, - "Ey": 0.0, - "Ez": 0.0 + "Ex": 4414725.184731114, + "Ey": 6.323754812712063e-05, + "Ez": 6.323754812712063e-05 } } \ No newline at end of file diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index daeef409dc8..56a7fbc8f10 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -605,7 +605,7 @@ analysisRoutine = Examples/Tests/electrostatic_sphere/analysis_electrostatic_sph [embedded_boundary_cube] buildDir = . inputFile = Examples/Tests/embedded_boundary_cube/inputs_3d -runtime_params = +runtime_params = warpx.abort_on_warning_threshold = medium dim = 3 addToCompileString = USE_EB=TRUE cmakeSetupOpts = -DWarpX_DIMS=3 -DWarpX_EB=ON @@ -622,7 +622,7 @@ analysisRoutine = Examples/Tests/embedded_boundary_cube/analysis_fields.py [embedded_boundary_cube_2d] buildDir = . inputFile = Examples/Tests/embedded_boundary_cube/inputs_2d -runtime_params = +runtime_params = warpx.abort_on_warning_threshold = medium dim = 2 addToCompileString = USE_EB=TRUE cmakeSetupOpts = -DWarpX_DIMS=2 -DWarpX_EB=ON @@ -639,7 +639,7 @@ analysisRoutine = Examples/Tests/embedded_boundary_cube/analysis_fields_2d.py [embedded_boundary_cube_macroscopic] buildDir = . inputFile = Examples/Tests/embedded_boundary_cube/inputs_3d -runtime_params = algo.em_solver_medium=macroscopic macroscopic.epsilon=1.5*8.8541878128e-12 macroscopic.sigma=0 macroscopic.mu=1.25663706212e-06 +runtime_params = algo.em_solver_medium=macroscopic macroscopic.epsilon=1.5*8.8541878128e-12 macroscopic.sigma=0 macroscopic.mu=1.25663706212e-06 warpx.abort_on_warning_threshold=medium dim = 3 addToCompileString = USE_EB=TRUE cmakeSetupOpts = -DWarpX_DIMS=3 -DWarpX_EB=ON diff --git a/Source/FieldSolver/ElectrostaticSolver.cpp b/Source/FieldSolver/ElectrostaticSolver.cpp index bd3086c0438..3e3a1729736 100644 --- a/Source/FieldSolver/ElectrostaticSolver.cpp +++ b/Source/FieldSolver/ElectrostaticSolver.cpp @@ -88,7 +88,8 @@ WarpX::ComputeSpaceChargeField (bool const reset_fields) } // Add the field due to the boundary potentials - if (electrostatic_solver_id == ElectrostaticSolverAlgo::Relativistic){ + if (m_boundary_potential_specified || + (electrostatic_solver_id == ElectrostaticSolverAlgo::Relativistic)){ AddBoundaryField(); } } diff --git a/Source/WarpX.H b/Source/WarpX.H index d3900c91593..6526171ff67 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -954,6 +954,7 @@ public: */ [[nodiscard]] amrex::IntVect get_numprocs() const {return numprocs;} + bool m_boundary_potential_specified = false; ElectrostaticSolver::PoissonBoundaryHandler m_poisson_boundary_handler; void ComputeSpaceChargeField (bool reset_fields); void AddBoundaryField (); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 6c3186c5c8d..c929718b043 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -765,13 +765,34 @@ WarpX::ReadParameters () } // Parse the input file for domain boundary potentials const ParmParse pp_boundary("boundary"); - pp_boundary.query("potential_lo_x", m_poisson_boundary_handler.potential_xlo_str); - pp_boundary.query("potential_hi_x", m_poisson_boundary_handler.potential_xhi_str); - pp_boundary.query("potential_lo_y", m_poisson_boundary_handler.potential_ylo_str); - pp_boundary.query("potential_hi_y", m_poisson_boundary_handler.potential_yhi_str); - pp_boundary.query("potential_lo_z", m_poisson_boundary_handler.potential_zlo_str); - pp_boundary.query("potential_hi_z", m_poisson_boundary_handler.potential_zhi_str); - pp_warpx.query("eb_potential(x,y,z,t)", m_poisson_boundary_handler.potential_eb_str); + bool potential_specified = false; + // When reading the potential at the boundary from the input file, set this flag to true if any of the potential is specified + potential_specified |= pp_boundary.query("potential_lo_x", m_poisson_boundary_handler.potential_xlo_str); + potential_specified |= pp_boundary.query("potential_hi_x", m_poisson_boundary_handler.potential_xhi_str); + potential_specified |= pp_boundary.query("potential_lo_y", m_poisson_boundary_handler.potential_ylo_str); + potential_specified |= pp_boundary.query("potential_hi_y", m_poisson_boundary_handler.potential_yhi_str); + potential_specified |= pp_boundary.query("potential_lo_z", m_poisson_boundary_handler.potential_zlo_str); + potential_specified |= pp_boundary.query("potential_hi_z", m_poisson_boundary_handler.potential_zhi_str); +#if defined(AMREX_USE_EB) + potential_specified |= pp_warpx.query("eb_potential(x,y,z,t)", m_poisson_boundary_handler.potential_eb_str); +#endif + m_boundary_potential_specified = potential_specified; + if (potential_specified & (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC)) { + ablastr::warn_manager::WMRecordWarning( + "Algorithms", + "The input script specifies the electric potential (phi) at the boundary, but \ + also uses the hybrid PIC solver based on Ohm’s law. When using this solver, the \ + electric potential does not have any impact on the simulation.", + ablastr::warn_manager::WarnPriority::low); + } + else if (potential_specified & (WarpX::electromagnetic_solver_id != ElectromagneticSolverAlgo::None)) { + ablastr::warn_manager::WMRecordWarning( + "Algorithms", + "The input script specifies the electric potential (phi) at the boundary so \ + an initial Poisson solve will be performed.", + ablastr::warn_manager::WarnPriority::low); + } + m_poisson_boundary_handler.buildParsers(); #ifdef WARPX_DIM_RZ pp_boundary.query("verboncoeur_axis_correction", verboncoeur_axis_correction); From e145007861cdb7988fd99a095055ab38484b99ed Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Thu, 14 Mar 2024 12:31:11 -0700 Subject: [PATCH 08/14] Better error message when Poisson solver does not have supported boundary (#4778) * Better boundary error message * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- Source/FieldSolver/ElectrostaticSolver.cpp | 8 ++++---- .../MagnetostaticSolver/MagnetostaticSolver.cpp | 12 ++++++------ Source/Utils/WarpXAlgorithmSelection.H | 5 +++++ Source/Utils/WarpXAlgorithmSelection.cpp | 12 ++++++++++++ 4 files changed, 27 insertions(+), 10 deletions(-) diff --git a/Source/FieldSolver/ElectrostaticSolver.cpp b/Source/FieldSolver/ElectrostaticSolver.cpp index 3e3a1729736..f1bda102c8b 100644 --- a/Source/FieldSolver/ElectrostaticSolver.cpp +++ b/Source/FieldSolver/ElectrostaticSolver.cpp @@ -1019,9 +1019,9 @@ void ElectrostaticSolver::PoissonBoundaryHandler::definePhiBCs (const amrex::Geo dirichlet_flag[idim*2] = false; } else { - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(false, + WARPX_ABORT_WITH_MESSAGE( "Field boundary conditions have to be either periodic, PEC or neumann " - "when using the electrostatic solver" + "when using the electrostatic solver, but they are " + GetFieldBCTypeString(WarpX::field_boundary_lo[idim]) ); } @@ -1034,9 +1034,9 @@ void ElectrostaticSolver::PoissonBoundaryHandler::definePhiBCs (const amrex::Geo dirichlet_flag[idim*2+1] = false; } else { - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(false, + WARPX_ABORT_WITH_MESSAGE( "Field boundary conditions have to be either periodic, PEC or neumann " - "when using the electrostatic solver" + "when using the electrostatic solver, but they are " + GetFieldBCTypeString(WarpX::field_boundary_hi[idim]) ); } } diff --git a/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.cpp b/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.cpp index 192017656ce..26ac1ac96c8 100644 --- a/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.cpp +++ b/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.cpp @@ -318,9 +318,9 @@ void MagnetostaticSolver::VectorPoissonBoundaryHandler::defineVectorPotentialBCs dirichlet_flag[adim][idim*2] = false; } else { - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(false, - "Field boundary conditions have to be either periodic, PEC, or Neumann " - "when using the magnetostatic solver" + WARPX_ABORT_WITH_MESSAGE( + "Field boundary conditions have to be either periodic, PEC or neumann " + "when using the magnetostatic solver, but they are " + GetFieldBCTypeString(WarpX::field_boundary_lo[idim]) ); } @@ -338,9 +338,9 @@ void MagnetostaticSolver::VectorPoissonBoundaryHandler::defineVectorPotentialBCs dirichlet_flag[adim][idim*2+1] = false; } else { - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(false, - "Field boundary conditions have to be either periodic, PEC, or Neumann " - "when using the magnetostatic solver" + WARPX_ABORT_WITH_MESSAGE( + "Field boundary conditions have to be either periodic, PEC or neumann " + "when using the magnetostatic solver, but they are " + GetFieldBCTypeString(WarpX::field_boundary_lo[idim]) ); } } diff --git a/Source/Utils/WarpXAlgorithmSelection.H b/Source/Utils/WarpXAlgorithmSelection.H index 735fc7993f1..e94b7bb7719 100644 --- a/Source/Utils/WarpXAlgorithmSelection.H +++ b/Source/Utils/WarpXAlgorithmSelection.H @@ -188,4 +188,9 @@ GetFieldBCTypeInteger( std::string BCType ); ParticleBoundaryType GetParticleBCTypeInteger( std::string BCType ); +/** Find the name associated with a BC type + */ +std::string +GetFieldBCTypeString( int fb_type ); + #endif // UTILS_WARPXALGORITHMSELECTION_H_ diff --git a/Source/Utils/WarpXAlgorithmSelection.cpp b/Source/Utils/WarpXAlgorithmSelection.cpp index abaf17f0a2c..b874c151f36 100644 --- a/Source/Utils/WarpXAlgorithmSelection.cpp +++ b/Source/Utils/WarpXAlgorithmSelection.cpp @@ -246,3 +246,15 @@ GetParticleBCTypeInteger( std::string BCType ){ // return ParticleBCType_algo_to_enum[BCType]; // This operator cannot be used for a const map return ParticleBCType_algo_to_enum.at(BCType); } + +std::string +GetFieldBCTypeString( int fb_type ) { + std::string boundary_name = ""; + for (const auto &valid_pair : FieldBCType_algo_to_int) { + if ((valid_pair.second == fb_type)&&(valid_pair.first != "default")){ + boundary_name = valid_pair.first; + break; + } + } + return boundary_name; +} From 8c84e792db66ef8e2d99a4068daf837972d0a7aa Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Fri, 15 Mar 2024 09:57:21 -0700 Subject: [PATCH 09/14] CI: `isort` + `black` compatibility (#4769) * CI: `isort` + `black` compatibility Avoid that these tools follow slightly different rules. * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Remove Unused Import Co-authored-by: Luca Fedeli --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Luca Fedeli --- .pre-commit-config.yaml | 1 + Docs/source/conf.py | 2 +- .../workflows/ml_materials/create_dataset.py | 2 +- .../ml_materials/run_warpx_training.py | 1 - .../usage/workflows/ml_materials/visualize.py | 2 +- .../capacitive_discharge/PICMI_inputs_1d.py | 3 +-- .../capacitive_discharge/PICMI_inputs_2d.py | 3 +-- .../Physics_applications/laser_ion/plot_2d.py | 4 ++-- .../spacecraft_charging/PICMI_inputs_rz.py | 3 +-- .../spacecraft_charging/analysis.py | 2 +- Examples/Tests/AcceleratorLattice/analysis.py | 2 +- Examples/Tests/Implicit/analysis_vandb_2d.py | 2 +- Examples/Tests/boosted_diags/analysis.py | 2 +- Examples/Tests/boundaries/analysis.py | 2 +- .../Tests/collision/analysis_collision_rz.py | 2 +- Examples/Tests/dive_cleaning/analysis.py | 2 +- .../analysis_electrostatic_sphere.py | 2 +- .../Tests/electrostatic_sphere_eb/analysis.py | 1 + .../electrostatic_sphere_eb/analysis_rz.py | 2 +- .../embedded_boundary_cube/analysis_fields.py | 2 +- .../analysis_fields_2d.py | 2 +- .../PICMI_inputs_EB_API.py | 1 - .../analysis_fields.py | 2 +- .../analysis_fields_2d.py | 2 +- .../analysis_flux_injection_3d.py | 2 +- .../analysis_distribution.py | 2 +- .../ion_stopping/analysis_ion_stopping.py | 2 +- Examples/Tests/langmuir/PICMI_inputs_rz.py | 1 - Examples/Tests/langmuir/analysis_2d.py | 2 +- Examples/Tests/langmuir/analysis_3d.py | 2 +- Examples/Tests/langmuir_fluids/analysis_2d.py | 2 +- Examples/Tests/langmuir_fluids/analysis_3d.py | 2 +- Examples/Tests/laser_injection/analysis_2d.py | 4 ++-- .../analysis_from_RZ_file.py | 5 +---- .../Tests/magnetostatic_eb/PICMI_inputs_3d.py | 1 - .../Tests/magnetostatic_eb/PICMI_inputs_rz.py | 1 - .../Tests/ohm_solver_EM_modes/PICMI_inputs.py | 3 +-- .../ohm_solver_EM_modes/PICMI_inputs_rz.py | 3 +-- .../Tests/ohm_solver_EM_modes/analysis.py | 1 - .../Tests/ohm_solver_EM_modes/analysis_rz.py | 7 +++---- .../PICMI_inputs.py | 3 +-- .../ohm_solver_ion_Landau_damping/analysis.py | 1 - .../PICMI_inputs.py | 3 +-- .../analysis.py | 1 - .../PICMI_inputs.py | 5 ++--- .../analysis.py | 2 +- .../PICMI_inputs_rz.py | 1 - .../particle_boundary_interaction/analysis.py | 2 +- .../PICMI_inputs_scrape.py | 1 - .../particle_data_python/PICMI_inputs_2d.py | 1 - .../PICMI_inputs_prev_pos_2d.py | 1 - .../analysis_particle_diags_impl.py | 2 +- .../pass_mpi_communicator/PICMI_inputs_2d.py | 1 - Examples/Tests/plasma_lens/analysis.py | 2 +- .../Tests/point_of_contact_EB/analysis.py | 2 +- .../Tests/python_wrappers/PICMI_inputs_2d.py | 3 +-- .../analysis_reduced_diags_impl.py | 2 +- .../repelling_particles/analysis_repelling.py | 2 +- .../resampling/analysis_leveling_thinning.py | 2 +- .../Tests/restart/PICMI_inputs_id_cpu_read.py | 1 - .../PICMI_inputs_runtime_component_analyze.py | 1 - .../analysis_rigid_injection_BoostedFrame.py | 1 - Examples/Tests/scraping/analysis_rz.py | 2 +- Examples/Tests/scraping/analysis_rz_filter.py | 2 +- .../space_charge_initialization/analysis.py | 2 +- Examples/Tests/vay_deposition/analysis.py | 2 +- Python/pywarpx/WarpX.py | 4 ++-- Python/pywarpx/__init__.py | 4 ++-- Python/pywarpx/particle_containers.py | 2 +- Python/pywarpx/picmi.py | 3 +-- Regression/Checksum/checksum.py | 4 ++-- Regression/prepare_file_ci.py | 1 + Tools/Algorithms/stencil.py | 2 +- Tools/LibEnsemble/run_libensemble_on_warpx.py | 19 +++++++++---------- Tools/LibEnsemble/warpx_simf.py | 2 +- Tools/PerformanceTests/run_automated.py | 4 ++-- .../plot_distribution_mapping.py | 2 +- Tools/PostProcessing/yt3d_mpi.py | 2 +- Tools/Release/updateAMReX.py | 2 +- Tools/Release/updatePICSAR.py | 2 +- Tools/Release/updatepyAMReX.py | 2 +- pyproject.toml | 3 +++ 82 files changed, 85 insertions(+), 110 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index f3d14c123c1..812a046b602 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -80,6 +80,7 @@ repos: hooks: - id: isort name: isort (python) + args: ['--profile black'] # Python: Flake8 (checks only, does this support auto-fixes?) #- repo: https://github.com/PyCQA/flake8 diff --git a/Docs/source/conf.py b/Docs/source/conf.py index 1a439c8b76e..10752c21bda 100644 --- a/Docs/source/conf.py +++ b/Docs/source/conf.py @@ -31,8 +31,8 @@ import urllib.request import pybtex.plugin -from pybtex.style.formatting.unsrt import Style as UnsrtStyle import sphinx_rtd_theme +from pybtex.style.formatting.unsrt import Style as UnsrtStyle sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), '../../Regression/Checksum')) diff --git a/Docs/source/usage/workflows/ml_materials/create_dataset.py b/Docs/source/usage/workflows/ml_materials/create_dataset.py index 08000105479..1dd16f2f95e 100644 --- a/Docs/source/usage/workflows/ml_materials/create_dataset.py +++ b/Docs/source/usage/workflows/ml_materials/create_dataset.py @@ -15,8 +15,8 @@ c = 2.998e8 -from openpmd_viewer import OpenPMDTimeSeries, ParticleTracker import torch +from openpmd_viewer import OpenPMDTimeSeries, ParticleTracker ############### diff --git a/Docs/source/usage/workflows/ml_materials/run_warpx_training.py b/Docs/source/usage/workflows/ml_materials/run_warpx_training.py index 8ee098e3e29..772fe574504 100644 --- a/Docs/source/usage/workflows/ml_materials/run_warpx_training.py +++ b/Docs/source/usage/workflows/ml_materials/run_warpx_training.py @@ -9,7 +9,6 @@ import math import numpy as np - from pywarpx import picmi # Physical constants diff --git a/Docs/source/usage/workflows/ml_materials/visualize.py b/Docs/source/usage/workflows/ml_materials/visualize.py index e6da1029b67..920efe29909 100644 --- a/Docs/source/usage/workflows/ml_materials/visualize.py +++ b/Docs/source/usage/workflows/ml_materials/visualize.py @@ -7,11 +7,11 @@ # Authors: Ryan Sandberg # License: BSD-3-Clause-LBNL # -from matplotlib import pyplot as plt import neural_network_classes as mynn import numpy as np import torch import torch.nn.functional as F +from matplotlib import pyplot as plt c = 2.998e8 diff --git a/Examples/Physics_applications/capacitive_discharge/PICMI_inputs_1d.py b/Examples/Physics_applications/capacitive_discharge/PICMI_inputs_1d.py index 01406b7a0fb..d54800debd6 100644 --- a/Examples/Physics_applications/capacitive_discharge/PICMI_inputs_1d.py +++ b/Examples/Physics_applications/capacitive_discharge/PICMI_inputs_1d.py @@ -8,11 +8,10 @@ import sys import numpy as np +from pywarpx import callbacks, fields, libwarpx, particle_containers, picmi from scipy.sparse import csc_matrix from scipy.sparse import linalg as sla -from pywarpx import callbacks, fields, libwarpx, particle_containers, picmi - constants = picmi.constants diff --git a/Examples/Physics_applications/capacitive_discharge/PICMI_inputs_2d.py b/Examples/Physics_applications/capacitive_discharge/PICMI_inputs_2d.py index 65baabba605..61229a17295 100755 --- a/Examples/Physics_applications/capacitive_discharge/PICMI_inputs_2d.py +++ b/Examples/Physics_applications/capacitive_discharge/PICMI_inputs_2d.py @@ -6,11 +6,10 @@ # --- used for the field solve step. import numpy as np +from pywarpx import callbacks, fields, picmi from scipy.sparse import csc_matrix from scipy.sparse import linalg as sla -from pywarpx import callbacks, fields, picmi - constants = picmi.constants ########################## diff --git a/Examples/Physics_applications/laser_ion/plot_2d.py b/Examples/Physics_applications/laser_ion/plot_2d.py index e85782a7d23..736203e85ea 100644 --- a/Examples/Physics_applications/laser_ion/plot_2d.py +++ b/Examples/Physics_applications/laser_ion/plot_2d.py @@ -14,12 +14,12 @@ import os import re -from matplotlib.colors import TwoSlopeNorm import matplotlib.pyplot as plt import numpy as np -from openpmd_viewer import OpenPMDTimeSeries import pandas as pd import scipy.constants as sc +from matplotlib.colors import TwoSlopeNorm +from openpmd_viewer import OpenPMDTimeSeries plt.rcParams.update({'font.size':16}) diff --git a/Examples/Physics_applications/spacecraft_charging/PICMI_inputs_rz.py b/Examples/Physics_applications/spacecraft_charging/PICMI_inputs_rz.py index 08ebe23f828..3de56a93710 100644 --- a/Examples/Physics_applications/spacecraft_charging/PICMI_inputs_rz.py +++ b/Examples/Physics_applications/spacecraft_charging/PICMI_inputs_rz.py @@ -12,10 +12,9 @@ # --- of electrons - leads to a decrease of the potential on the surface over the time # --- until reaching an equilibrium floating potential of ~144.5 V (*). -from mpi4py import MPI as mpi import numpy as np import scipy.constants as scc - +from mpi4py import MPI as mpi from pywarpx import picmi from pywarpx.callbacks import installafterEsolve, installafterInitEsolve from pywarpx.fields import ExWrapper, EzWrapper, PhiFPWrapper, RhoFPWrapper diff --git a/Examples/Physics_applications/spacecraft_charging/analysis.py b/Examples/Physics_applications/spacecraft_charging/analysis.py index d8e8ac86af8..10786baab61 100755 --- a/Examples/Physics_applications/spacecraft_charging/analysis.py +++ b/Examples/Physics_applications/spacecraft_charging/analysis.py @@ -17,9 +17,9 @@ import matplotlib.pyplot as plt import numpy as np +import yt from openpmd_viewer import OpenPMDTimeSeries from scipy.optimize import curve_fit -import yt yt.funcs.mylog.setLevel(0) sys.path.insert(1, '../../../../warpx/Regression/Checksum/') diff --git a/Examples/Tests/AcceleratorLattice/analysis.py b/Examples/Tests/AcceleratorLattice/analysis.py index 55fb1b2ac96..5f3de4543d6 100755 --- a/Examples/Tests/AcceleratorLattice/analysis.py +++ b/Examples/Tests/AcceleratorLattice/analysis.py @@ -19,8 +19,8 @@ import sys import numpy as np -from scipy.constants import c, e, m_e import yt +from scipy.constants import c, e, m_e yt.funcs.mylog.setLevel(0) sys.path.insert(1, '../../../../warpx/Regression/Checksum/') diff --git a/Examples/Tests/Implicit/analysis_vandb_2d.py b/Examples/Tests/Implicit/analysis_vandb_2d.py index fa3299925a8..85faab61fcc 100755 --- a/Examples/Tests/Implicit/analysis_vandb_2d.py +++ b/Examples/Tests/Implicit/analysis_vandb_2d.py @@ -14,8 +14,8 @@ import sys import numpy as np -from scipy.constants import e, epsilon_0 import yt +from scipy.constants import e, epsilon_0 sys.path.insert(1, '../../../../warpx/Regression/Checksum/') import checksumAPI diff --git a/Examples/Tests/boosted_diags/analysis.py b/Examples/Tests/boosted_diags/analysis.py index c0b03f4a20b..2b21184dc3d 100755 --- a/Examples/Tests/boosted_diags/analysis.py +++ b/Examples/Tests/boosted_diags/analysis.py @@ -21,8 +21,8 @@ import numpy as np import openpmd_api as io -from openpmd_viewer import OpenPMDTimeSeries import yt +from openpmd_viewer import OpenPMDTimeSeries yt.funcs.mylog.setLevel(0) diff --git a/Examples/Tests/boundaries/analysis.py b/Examples/Tests/boundaries/analysis.py index 0a879bb83e5..9c108b16196 100755 --- a/Examples/Tests/boundaries/analysis.py +++ b/Examples/Tests/boundaries/analysis.py @@ -18,8 +18,8 @@ import sys import numpy as np -from scipy.constants import c, m_e import yt +from scipy.constants import c, m_e yt.funcs.mylog.setLevel(0) sys.path.insert(1, '../../../../warpx/Regression/Checksum/') diff --git a/Examples/Tests/collision/analysis_collision_rz.py b/Examples/Tests/collision/analysis_collision_rz.py index 9f9b9c17f5e..b206b2eba7b 100755 --- a/Examples/Tests/collision/analysis_collision_rz.py +++ b/Examples/Tests/collision/analysis_collision_rz.py @@ -16,9 +16,9 @@ # tolerance: 1.0e-30 # Possible running time: ~ 1.0 s -from glob import glob import os import sys +from glob import glob import numpy as np import yt diff --git a/Examples/Tests/dive_cleaning/analysis.py b/Examples/Tests/dive_cleaning/analysis.py index 504368fa85e..a9b52455baa 100755 --- a/Examples/Tests/dive_cleaning/analysis.py +++ b/Examples/Tests/dive_cleaning/analysis.py @@ -22,8 +22,8 @@ import matplotlib.pyplot as plt import numpy as np import scipy.constants as scc -from scipy.special import gammainc import yt +from scipy.special import gammainc yt.funcs.mylog.setLevel(0) diff --git a/Examples/Tests/electrostatic_sphere/analysis_electrostatic_sphere.py b/Examples/Tests/electrostatic_sphere/analysis_electrostatic_sphere.py index 097ea0674f7..35281084c01 100755 --- a/Examples/Tests/electrostatic_sphere/analysis_electrostatic_sphere.py +++ b/Examples/Tests/electrostatic_sphere/analysis_electrostatic_sphere.py @@ -21,8 +21,8 @@ import sys import numpy as np -from scipy.optimize import fsolve import yt +from scipy.optimize import fsolve sys.path.insert(1, '../../../../warpx/Regression/Checksum/') import checksumAPI diff --git a/Examples/Tests/electrostatic_sphere_eb/analysis.py b/Examples/Tests/electrostatic_sphere_eb/analysis.py index c1a5ba33a1b..bf2725616b5 100755 --- a/Examples/Tests/electrostatic_sphere_eb/analysis.py +++ b/Examples/Tests/electrostatic_sphere_eb/analysis.py @@ -9,6 +9,7 @@ sys.path.insert(1, '../../../../warpx/Regression/Checksum/') import checksumAPI + # Check reduced diagnostics for charge on EB import numpy as np from scipy.constants import epsilon_0 diff --git a/Examples/Tests/electrostatic_sphere_eb/analysis_rz.py b/Examples/Tests/electrostatic_sphere_eb/analysis_rz.py index 5d807a2dbd3..18aba4322f0 100755 --- a/Examples/Tests/electrostatic_sphere_eb/analysis_rz.py +++ b/Examples/Tests/electrostatic_sphere_eb/analysis_rz.py @@ -20,8 +20,8 @@ import sys import numpy as np -from unyt import m import yt +from unyt import m sys.path.insert(1, '../../../../warpx/Regression/Checksum/') import checksumAPI diff --git a/Examples/Tests/embedded_boundary_cube/analysis_fields.py b/Examples/Tests/embedded_boundary_cube/analysis_fields.py index 13bd32d73e4..dc6af9d57d2 100755 --- a/Examples/Tests/embedded_boundary_cube/analysis_fields.py +++ b/Examples/Tests/embedded_boundary_cube/analysis_fields.py @@ -5,8 +5,8 @@ import sys import numpy as np -from scipy.constants import c, mu_0, pi import yt +from scipy.constants import c, mu_0, pi sys.path.insert(1, '../../../../warpx/Regression/Checksum/') import checksumAPI diff --git a/Examples/Tests/embedded_boundary_cube/analysis_fields_2d.py b/Examples/Tests/embedded_boundary_cube/analysis_fields_2d.py index 67e893b199a..8faa299025e 100755 --- a/Examples/Tests/embedded_boundary_cube/analysis_fields_2d.py +++ b/Examples/Tests/embedded_boundary_cube/analysis_fields_2d.py @@ -4,8 +4,8 @@ import sys import numpy as np -from scipy.constants import c, mu_0, pi import yt +from scipy.constants import c, mu_0, pi sys.path.insert(1, '../../../../warpx/Regression/Checksum/') import checksumAPI diff --git a/Examples/Tests/embedded_boundary_python_api/PICMI_inputs_EB_API.py b/Examples/Tests/embedded_boundary_python_api/PICMI_inputs_EB_API.py index faec3ed4668..6104e08ef27 100755 --- a/Examples/Tests/embedded_boundary_python_api/PICMI_inputs_EB_API.py +++ b/Examples/Tests/embedded_boundary_python_api/PICMI_inputs_EB_API.py @@ -1,6 +1,5 @@ #!/usr/bin/env python3 import numpy as np - from pywarpx import fields, picmi max_steps = 1 diff --git a/Examples/Tests/embedded_boundary_rotated_cube/analysis_fields.py b/Examples/Tests/embedded_boundary_rotated_cube/analysis_fields.py index 1f7b604365c..e849958468f 100755 --- a/Examples/Tests/embedded_boundary_rotated_cube/analysis_fields.py +++ b/Examples/Tests/embedded_boundary_rotated_cube/analysis_fields.py @@ -11,8 +11,8 @@ import sys import numpy as np -from scipy.constants import c, mu_0, pi import yt +from scipy.constants import c, mu_0, pi sys.path.insert(1, '../../../../warpx/Regression/Checksum/') import checksumAPI diff --git a/Examples/Tests/embedded_boundary_rotated_cube/analysis_fields_2d.py b/Examples/Tests/embedded_boundary_rotated_cube/analysis_fields_2d.py index 5c8c544d0c4..dcdbc83a729 100755 --- a/Examples/Tests/embedded_boundary_rotated_cube/analysis_fields_2d.py +++ b/Examples/Tests/embedded_boundary_rotated_cube/analysis_fields_2d.py @@ -4,8 +4,8 @@ import sys import numpy as np -from scipy.constants import c, mu_0, pi import yt +from scipy.constants import c, mu_0, pi sys.path.insert(1, '../../../../warpx/Regression/Checksum/') import checksumAPI diff --git a/Examples/Tests/flux_injection/analysis_flux_injection_3d.py b/Examples/Tests/flux_injection/analysis_flux_injection_3d.py index 470cbfe6065..d0271f6aa94 100755 --- a/Examples/Tests/flux_injection/analysis_flux_injection_3d.py +++ b/Examples/Tests/flux_injection/analysis_flux_injection_3d.py @@ -26,9 +26,9 @@ import matplotlib.pyplot as plt import numpy as np +import yt from scipy.constants import c, m_e, m_p from scipy.special import erf -import yt sys.path.insert(1, '../../../../warpx/Regression/Checksum/') import checksumAPI diff --git a/Examples/Tests/initial_distribution/analysis_distribution.py b/Examples/Tests/initial_distribution/analysis_distribution.py index d5301c077d2..5a3774133db 100755 --- a/Examples/Tests/initial_distribution/analysis_distribution.py +++ b/Examples/Tests/initial_distribution/analysis_distribution.py @@ -22,9 +22,9 @@ import sys import numpy as np -from read_raw_data import read_reduced_diags, read_reduced_diags_histogram import scipy.constants as scc import scipy.special as scs +from read_raw_data import read_reduced_diags, read_reduced_diags_histogram sys.path.insert(1, '../../../../warpx/Regression/Checksum/') import checksumAPI diff --git a/Examples/Tests/ion_stopping/analysis_ion_stopping.py b/Examples/Tests/ion_stopping/analysis_ion_stopping.py index b090014c0cf..d7774c14d6b 100755 --- a/Examples/Tests/ion_stopping/analysis_ion_stopping.py +++ b/Examples/Tests/ion_stopping/analysis_ion_stopping.py @@ -16,8 +16,8 @@ import sys import numpy as np -from scipy.constants import e, epsilon_0, k, m_e, m_p import yt +from scipy.constants import e, epsilon_0, k, m_e, m_p sys.path.insert(1, '../../../../warpx/Regression/Checksum/') import checksumAPI diff --git a/Examples/Tests/langmuir/PICMI_inputs_rz.py b/Examples/Tests/langmuir/PICMI_inputs_rz.py index 8328aba5185..106b8de5a94 100755 --- a/Examples/Tests/langmuir/PICMI_inputs_rz.py +++ b/Examples/Tests/langmuir/PICMI_inputs_rz.py @@ -9,7 +9,6 @@ matplotlib.use('Agg') import matplotlib.pyplot as plt import numpy as np - from pywarpx import fields, picmi constants = picmi.constants diff --git a/Examples/Tests/langmuir/analysis_2d.py b/Examples/Tests/langmuir/analysis_2d.py index b3327703b82..94f97ca6de8 100755 --- a/Examples/Tests/langmuir/analysis_2d.py +++ b/Examples/Tests/langmuir/analysis_2d.py @@ -18,8 +18,8 @@ import sys import matplotlib.pyplot as plt -from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable import yt +from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable yt.funcs.mylog.setLevel(50) diff --git a/Examples/Tests/langmuir/analysis_3d.py b/Examples/Tests/langmuir/analysis_3d.py index 2e56e64d641..68334f506ff 100755 --- a/Examples/Tests/langmuir/analysis_3d.py +++ b/Examples/Tests/langmuir/analysis_3d.py @@ -18,8 +18,8 @@ import sys import matplotlib.pyplot as plt -from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable import yt +from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable yt.funcs.mylog.setLevel(50) diff --git a/Examples/Tests/langmuir_fluids/analysis_2d.py b/Examples/Tests/langmuir_fluids/analysis_2d.py index cf5d2fb44de..f7244f87137 100755 --- a/Examples/Tests/langmuir_fluids/analysis_2d.py +++ b/Examples/Tests/langmuir_fluids/analysis_2d.py @@ -17,8 +17,8 @@ import sys import matplotlib.pyplot as plt -from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable import yt +from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable yt.funcs.mylog.setLevel(50) diff --git a/Examples/Tests/langmuir_fluids/analysis_3d.py b/Examples/Tests/langmuir_fluids/analysis_3d.py index 0211956859e..686907f103a 100755 --- a/Examples/Tests/langmuir_fluids/analysis_3d.py +++ b/Examples/Tests/langmuir_fluids/analysis_3d.py @@ -18,8 +18,8 @@ import sys import matplotlib.pyplot as plt -from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable import yt +from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable yt.funcs.mylog.setLevel(50) diff --git a/Examples/Tests/laser_injection/analysis_2d.py b/Examples/Tests/laser_injection/analysis_2d.py index 37793171883..4424fe134bc 100755 --- a/Examples/Tests/laser_injection/analysis_2d.py +++ b/Examples/Tests/laser_injection/analysis_2d.py @@ -24,10 +24,10 @@ matplotlib.use('Agg') import matplotlib.pyplot as plt -from mpl_toolkits.axes_grid1 import make_axes_locatable import numpy as np -from scipy.signal import hilbert import yt +from mpl_toolkits.axes_grid1 import make_axes_locatable +from scipy.signal import hilbert sys.path.insert(1, '../../../../warpx/Regression/Checksum/') import checksumAPI diff --git a/Examples/Tests/laser_injection_from_file/analysis_from_RZ_file.py b/Examples/Tests/laser_injection_from_file/analysis_from_RZ_file.py index a00d322bd18..44bbd00a7f4 100755 --- a/Examples/Tests/laser_injection_from_file/analysis_from_RZ_file.py +++ b/Examples/Tests/laser_injection_from_file/analysis_from_RZ_file.py @@ -148,10 +148,7 @@ def launch_analysis(executable): def main() : from lasy.laser import Laser - from lasy.profiles import ( - CombinedLongitudinalTransverseProfile, - GaussianProfile, - ) + from lasy.profiles import CombinedLongitudinalTransverseProfile, GaussianProfile from lasy.profiles.longitudinal import GaussianLongitudinalProfile from lasy.profiles.transverse import LaguerreGaussianTransverseProfile diff --git a/Examples/Tests/magnetostatic_eb/PICMI_inputs_3d.py b/Examples/Tests/magnetostatic_eb/PICMI_inputs_3d.py index 8f205724563..f61823bf223 100755 --- a/Examples/Tests/magnetostatic_eb/PICMI_inputs_3d.py +++ b/Examples/Tests/magnetostatic_eb/PICMI_inputs_3d.py @@ -27,7 +27,6 @@ import matplotlib.pyplot as plt import numpy as np import scipy.constants as con - from pywarpx import fields, picmi ########################## diff --git a/Examples/Tests/magnetostatic_eb/PICMI_inputs_rz.py b/Examples/Tests/magnetostatic_eb/PICMI_inputs_rz.py index 1268f7a02b0..6f7d6daa453 100755 --- a/Examples/Tests/magnetostatic_eb/PICMI_inputs_rz.py +++ b/Examples/Tests/magnetostatic_eb/PICMI_inputs_rz.py @@ -27,7 +27,6 @@ import matplotlib.pyplot as plt import numpy as np import scipy.constants as con - from pywarpx import fields, picmi ########################## diff --git a/Examples/Tests/ohm_solver_EM_modes/PICMI_inputs.py b/Examples/Tests/ohm_solver_EM_modes/PICMI_inputs.py index 4f3cd731323..5101e218129 100644 --- a/Examples/Tests/ohm_solver_EM_modes/PICMI_inputs.py +++ b/Examples/Tests/ohm_solver_EM_modes/PICMI_inputs.py @@ -12,9 +12,8 @@ import sys import dill -from mpi4py import MPI as mpi import numpy as np - +from mpi4py import MPI as mpi from pywarpx import callbacks, fields, libwarpx, picmi constants = picmi.constants diff --git a/Examples/Tests/ohm_solver_EM_modes/PICMI_inputs_rz.py b/Examples/Tests/ohm_solver_EM_modes/PICMI_inputs_rz.py index 81ba65e5b28..e6f9ccc8e65 100644 --- a/Examples/Tests/ohm_solver_EM_modes/PICMI_inputs_rz.py +++ b/Examples/Tests/ohm_solver_EM_modes/PICMI_inputs_rz.py @@ -10,9 +10,8 @@ import sys import dill -from mpi4py import MPI as mpi import numpy as np - +from mpi4py import MPI as mpi from pywarpx import picmi constants = picmi.constants diff --git a/Examples/Tests/ohm_solver_EM_modes/analysis.py b/Examples/Tests/ohm_solver_EM_modes/analysis.py index ad832dc2f50..77be546937c 100755 --- a/Examples/Tests/ohm_solver_EM_modes/analysis.py +++ b/Examples/Tests/ohm_solver_EM_modes/analysis.py @@ -6,7 +6,6 @@ import matplotlib import matplotlib.pyplot as plt import numpy as np - from pywarpx import picmi constants = picmi.constants diff --git a/Examples/Tests/ohm_solver_EM_modes/analysis_rz.py b/Examples/Tests/ohm_solver_EM_modes/analysis_rz.py index 9004c24a05c..b96bc3b4629 100755 --- a/Examples/Tests/ohm_solver_EM_modes/analysis_rz.py +++ b/Examples/Tests/ohm_solver_EM_modes/analysis_rz.py @@ -3,16 +3,15 @@ # --- Analysis script for the hybrid-PIC example producing EM modes. import dill -from matplotlib import colors import matplotlib.pyplot as plt import numpy as np -from openpmd_viewer import OpenPMDTimeSeries import scipy.fft as fft +from matplotlib import colors +from openpmd_viewer import OpenPMDTimeSeries +from pywarpx import picmi from scipy.interpolate import RegularGridInterpolator from scipy.special import j1, jn, jn_zeros -from pywarpx import picmi - constants = picmi.constants # load simulation parameters diff --git a/Examples/Tests/ohm_solver_ion_Landau_damping/PICMI_inputs.py b/Examples/Tests/ohm_solver_ion_Landau_damping/PICMI_inputs.py index d545195a9e5..8249bdb504f 100644 --- a/Examples/Tests/ohm_solver_ion_Landau_damping/PICMI_inputs.py +++ b/Examples/Tests/ohm_solver_ion_Landau_damping/PICMI_inputs.py @@ -11,9 +11,8 @@ import time import dill -from mpi4py import MPI as mpi import numpy as np - +from mpi4py import MPI as mpi from pywarpx import callbacks, fields, libwarpx, particle_containers, picmi constants = picmi.constants diff --git a/Examples/Tests/ohm_solver_ion_Landau_damping/analysis.py b/Examples/Tests/ohm_solver_ion_Landau_damping/analysis.py index ade7187ebe3..93a4962d7e5 100755 --- a/Examples/Tests/ohm_solver_ion_Landau_damping/analysis.py +++ b/Examples/Tests/ohm_solver_ion_Landau_damping/analysis.py @@ -6,7 +6,6 @@ import matplotlib import matplotlib.pyplot as plt import numpy as np - from pywarpx import picmi constants = picmi.constants diff --git a/Examples/Tests/ohm_solver_ion_beam_instability/PICMI_inputs.py b/Examples/Tests/ohm_solver_ion_beam_instability/PICMI_inputs.py index d9a71df6ac4..186ef6cede3 100644 --- a/Examples/Tests/ohm_solver_ion_beam_instability/PICMI_inputs.py +++ b/Examples/Tests/ohm_solver_ion_beam_instability/PICMI_inputs.py @@ -12,9 +12,8 @@ import time import dill -from mpi4py import MPI as mpi import numpy as np - +from mpi4py import MPI as mpi from pywarpx import callbacks, fields, libwarpx, particle_containers, picmi constants = picmi.constants diff --git a/Examples/Tests/ohm_solver_ion_beam_instability/analysis.py b/Examples/Tests/ohm_solver_ion_beam_instability/analysis.py index 6a57b1c1046..6388a4c4085 100755 --- a/Examples/Tests/ohm_solver_ion_beam_instability/analysis.py +++ b/Examples/Tests/ohm_solver_ion_beam_instability/analysis.py @@ -7,7 +7,6 @@ import matplotlib import matplotlib.pyplot as plt import numpy as np - from pywarpx import picmi constants = picmi.constants diff --git a/Examples/Tests/ohm_solver_magnetic_reconnection/PICMI_inputs.py b/Examples/Tests/ohm_solver_magnetic_reconnection/PICMI_inputs.py index a7f2e0da12d..5d4631e9c1e 100644 --- a/Examples/Tests/ohm_solver_magnetic_reconnection/PICMI_inputs.py +++ b/Examples/Tests/ohm_solver_magnetic_reconnection/PICMI_inputs.py @@ -8,14 +8,13 @@ # --- https://aip.scitation.org/doi/10.1063/1.4943893. import argparse -from pathlib import Path import shutil import sys +from pathlib import Path import dill -from mpi4py import MPI as mpi import numpy as np - +from mpi4py import MPI as mpi from pywarpx import callbacks, fields, libwarpx, picmi constants = picmi.constants diff --git a/Examples/Tests/ohm_solver_magnetic_reconnection/analysis.py b/Examples/Tests/ohm_solver_magnetic_reconnection/analysis.py index e48d5c03fb2..0fb1c05ae1a 100755 --- a/Examples/Tests/ohm_solver_magnetic_reconnection/analysis.py +++ b/Examples/Tests/ohm_solver_magnetic_reconnection/analysis.py @@ -5,9 +5,9 @@ import glob import dill -from matplotlib import colors import matplotlib.pyplot as plt import numpy as np +from matplotlib import colors plt.rcParams.update({'font.size': 20}) diff --git a/Examples/Tests/particle_boundary_interaction/PICMI_inputs_rz.py b/Examples/Tests/particle_boundary_interaction/PICMI_inputs_rz.py index 4017a05d413..b9ea76e4a09 100644 --- a/Examples/Tests/particle_boundary_interaction/PICMI_inputs_rz.py +++ b/Examples/Tests/particle_boundary_interaction/PICMI_inputs_rz.py @@ -4,7 +4,6 @@ # --- This input is a simple case of reflection # --- of one electron on the surface of a sphere. import numpy as np - from pywarpx import callbacks, particle_containers, picmi ########################## diff --git a/Examples/Tests/particle_boundary_interaction/analysis.py b/Examples/Tests/particle_boundary_interaction/analysis.py index 9768751ed74..ff83cc1fed7 100755 --- a/Examples/Tests/particle_boundary_interaction/analysis.py +++ b/Examples/Tests/particle_boundary_interaction/analysis.py @@ -11,8 +11,8 @@ import sys import numpy as np -from openpmd_viewer import OpenPMDTimeSeries import yt +from openpmd_viewer import OpenPMDTimeSeries yt.funcs.mylog.setLevel(0) sys.path.insert(1, '../../../../warpx/Regression/Checksum/') diff --git a/Examples/Tests/particle_boundary_scrape/PICMI_inputs_scrape.py b/Examples/Tests/particle_boundary_scrape/PICMI_inputs_scrape.py index e5a9a58f597..63a805c2612 100755 --- a/Examples/Tests/particle_boundary_scrape/PICMI_inputs_scrape.py +++ b/Examples/Tests/particle_boundary_scrape/PICMI_inputs_scrape.py @@ -4,7 +4,6 @@ # --- to access the buffer of scraped particles. import numpy as np - from pywarpx import libwarpx, particle_containers, picmi ########################## diff --git a/Examples/Tests/particle_data_python/PICMI_inputs_2d.py b/Examples/Tests/particle_data_python/PICMI_inputs_2d.py index 572871b8ed5..f43c2215b79 100755 --- a/Examples/Tests/particle_data_python/PICMI_inputs_2d.py +++ b/Examples/Tests/particle_data_python/PICMI_inputs_2d.py @@ -3,7 +3,6 @@ import sys import numpy as np - from pywarpx import callbacks, libwarpx, particle_containers, picmi # Create the parser and add the argument diff --git a/Examples/Tests/particle_data_python/PICMI_inputs_prev_pos_2d.py b/Examples/Tests/particle_data_python/PICMI_inputs_prev_pos_2d.py index 5de9879f0f8..6f70b69fc10 100755 --- a/Examples/Tests/particle_data_python/PICMI_inputs_prev_pos_2d.py +++ b/Examples/Tests/particle_data_python/PICMI_inputs_prev_pos_2d.py @@ -3,7 +3,6 @@ # --- Input file to test the saving of old particle positions import numpy as np - from pywarpx import particle_containers, picmi constants = picmi.constants diff --git a/Examples/Tests/particle_fields_diags/analysis_particle_diags_impl.py b/Examples/Tests/particle_fields_diags/analysis_particle_diags_impl.py index ceb71ffb4b7..5e54fc42d87 100755 --- a/Examples/Tests/particle_fields_diags/analysis_particle_diags_impl.py +++ b/Examples/Tests/particle_fields_diags/analysis_particle_diags_impl.py @@ -16,8 +16,8 @@ import numpy as np import openpmd_api as io -from scipy.constants import c, e, m_e, m_p import yt +from scipy.constants import c, e, m_e, m_p sys.path.insert(1, '../../../../warpx/Regression/Checksum/') import checksumAPI diff --git a/Examples/Tests/pass_mpi_communicator/PICMI_inputs_2d.py b/Examples/Tests/pass_mpi_communicator/PICMI_inputs_2d.py index 66f259da2ef..095702a1145 100755 --- a/Examples/Tests/pass_mpi_communicator/PICMI_inputs_2d.py +++ b/Examples/Tests/pass_mpi_communicator/PICMI_inputs_2d.py @@ -6,7 +6,6 @@ # --- if the correct amount of processors are initialized in AMReX. from mpi4py import MPI - from pywarpx import picmi constants = picmi.constants diff --git a/Examples/Tests/plasma_lens/analysis.py b/Examples/Tests/plasma_lens/analysis.py index 9ce82d903f8..212e71087f9 100755 --- a/Examples/Tests/plasma_lens/analysis.py +++ b/Examples/Tests/plasma_lens/analysis.py @@ -19,8 +19,8 @@ import sys import numpy as np -from scipy.constants import c, e, m_e import yt +from scipy.constants import c, e, m_e yt.funcs.mylog.setLevel(0) sys.path.insert(1, '../../../../warpx/Regression/Checksum/') diff --git a/Examples/Tests/point_of_contact_EB/analysis.py b/Examples/Tests/point_of_contact_EB/analysis.py index cb1fc23ee62..042fc811a62 100755 --- a/Examples/Tests/point_of_contact_EB/analysis.py +++ b/Examples/Tests/point_of_contact_EB/analysis.py @@ -11,8 +11,8 @@ import sys import numpy as np -from openpmd_viewer import OpenPMDTimeSeries import yt +from openpmd_viewer import OpenPMDTimeSeries yt.funcs.mylog.setLevel(0) sys.path.insert(1, '../../../../warpx/Regression/Checksum/') diff --git a/Examples/Tests/python_wrappers/PICMI_inputs_2d.py b/Examples/Tests/python_wrappers/PICMI_inputs_2d.py index 7104963c752..d01835f9fdf 100755 --- a/Examples/Tests/python_wrappers/PICMI_inputs_2d.py +++ b/Examples/Tests/python_wrappers/PICMI_inputs_2d.py @@ -1,9 +1,8 @@ #!/usr/bin/env python3 import matplotlib.pyplot as plt -from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable import numpy as np - +from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable from pywarpx import picmi # Number of time steps diff --git a/Examples/Tests/reduced_diags/analysis_reduced_diags_impl.py b/Examples/Tests/reduced_diags/analysis_reduced_diags_impl.py index 93352d29ff3..21359ed8171 100755 --- a/Examples/Tests/reduced_diags/analysis_reduced_diags_impl.py +++ b/Examples/Tests/reduced_diags/analysis_reduced_diags_impl.py @@ -15,11 +15,11 @@ import sys import numpy as np +import yt from scipy.constants import c from scipy.constants import epsilon_0 as eps0 from scipy.constants import m_e, m_p from scipy.constants import mu_0 as mu0 -import yt sys.path.insert(1, '../../../../warpx/Regression/Checksum/') import checksumAPI diff --git a/Examples/Tests/repelling_particles/analysis_repelling.py b/Examples/Tests/repelling_particles/analysis_repelling.py index 0aa3ed53d1b..ebff010fc40 100755 --- a/Examples/Tests/repelling_particles/analysis_repelling.py +++ b/Examples/Tests/repelling_particles/analysis_repelling.py @@ -29,8 +29,8 @@ import matplotlib.pyplot as plt import numpy as np -from scipy.constants import c, m_e, physical_constants import yt +from scipy.constants import c, m_e, physical_constants yt.funcs.mylog.setLevel(0) diff --git a/Examples/Tests/resampling/analysis_leveling_thinning.py b/Examples/Tests/resampling/analysis_leveling_thinning.py index 50fec628200..5f3dc8ecdff 100755 --- a/Examples/Tests/resampling/analysis_leveling_thinning.py +++ b/Examples/Tests/resampling/analysis_leveling_thinning.py @@ -13,8 +13,8 @@ import sys import numpy as np -from scipy.special import erf import yt +from scipy.special import erf sys.path.insert(1, '../../../../warpx/Regression/Checksum/') import checksumAPI diff --git a/Examples/Tests/restart/PICMI_inputs_id_cpu_read.py b/Examples/Tests/restart/PICMI_inputs_id_cpu_read.py index d400924a378..d76cef82fec 100755 --- a/Examples/Tests/restart/PICMI_inputs_id_cpu_read.py +++ b/Examples/Tests/restart/PICMI_inputs_id_cpu_read.py @@ -6,7 +6,6 @@ import sys import numpy as np - from pywarpx import callbacks, particle_containers, picmi ########################## diff --git a/Examples/Tests/restart/PICMI_inputs_runtime_component_analyze.py b/Examples/Tests/restart/PICMI_inputs_runtime_component_analyze.py index 32c9f4e5808..daac6e1b5ac 100755 --- a/Examples/Tests/restart/PICMI_inputs_runtime_component_analyze.py +++ b/Examples/Tests/restart/PICMI_inputs_runtime_component_analyze.py @@ -7,7 +7,6 @@ import sys import numpy as np - from pywarpx import callbacks, particle_containers, picmi ########################## diff --git a/Examples/Tests/rigid_injection/analysis_rigid_injection_BoostedFrame.py b/Examples/Tests/rigid_injection/analysis_rigid_injection_BoostedFrame.py index ccb55183241..fd51298816b 100755 --- a/Examples/Tests/rigid_injection/analysis_rigid_injection_BoostedFrame.py +++ b/Examples/Tests/rigid_injection/analysis_rigid_injection_BoostedFrame.py @@ -25,7 +25,6 @@ import numpy as np import openpmd_api as io -from scipy.constants import m_e import yt yt.funcs.mylog.setLevel(0) diff --git a/Examples/Tests/scraping/analysis_rz.py b/Examples/Tests/scraping/analysis_rz.py index 7d40bd5edf4..193b2575992 100755 --- a/Examples/Tests/scraping/analysis_rz.py +++ b/Examples/Tests/scraping/analysis_rz.py @@ -24,8 +24,8 @@ import sys import numpy as np -from openpmd_viewer import OpenPMDTimeSeries import yt +from openpmd_viewer import OpenPMDTimeSeries sys.path.insert(1, '../../../../warpx/Regression/Checksum/') import checksumAPI diff --git a/Examples/Tests/scraping/analysis_rz_filter.py b/Examples/Tests/scraping/analysis_rz_filter.py index dd13b993dc5..a4e0dafddbc 100755 --- a/Examples/Tests/scraping/analysis_rz_filter.py +++ b/Examples/Tests/scraping/analysis_rz_filter.py @@ -23,8 +23,8 @@ import sys import numpy as np -from openpmd_viewer import OpenPMDTimeSeries import yt +from openpmd_viewer import OpenPMDTimeSeries tolerance = 0 diff --git a/Examples/Tests/space_charge_initialization/analysis.py b/Examples/Tests/space_charge_initialization/analysis.py index 0df21e6d984..48e19ea0b75 100755 --- a/Examples/Tests/space_charge_initialization/analysis.py +++ b/Examples/Tests/space_charge_initialization/analysis.py @@ -20,8 +20,8 @@ import matplotlib.pyplot as plt import numpy as np import scipy.constants as scc -from scipy.special import gammainc import yt +from scipy.special import gammainc yt.funcs.mylog.setLevel(0) sys.path.insert(1, '../../../../warpx/Regression/Checksum/') diff --git a/Examples/Tests/vay_deposition/analysis.py b/Examples/Tests/vay_deposition/analysis.py index f1680742ded..cfd089f2112 100755 --- a/Examples/Tests/vay_deposition/analysis.py +++ b/Examples/Tests/vay_deposition/analysis.py @@ -10,8 +10,8 @@ import sys import numpy as np -from scipy.constants import epsilon_0 import yt +from scipy.constants import epsilon_0 yt.funcs.mylog.setLevel(50) diff --git a/Python/pywarpx/WarpX.py b/Python/pywarpx/WarpX.py index 3fa59285f26..6b24b35e918 100644 --- a/Python/pywarpx/WarpX.py +++ b/Python/pywarpx/WarpX.py @@ -9,6 +9,7 @@ import sys from . import Particles +from ._libwarpx import libwarpx from .Algo import algo from .Amr import amr from .Amrex import amrex @@ -22,9 +23,8 @@ from .HybridPICModel import hybridpicmodel from .Interpolation import interpolation from .Lasers import lasers, lasers_list -from .PSATD import psatd from .Particles import particles, particles_list -from ._libwarpx import libwarpx +from .PSATD import psatd class WarpX(Bucket): diff --git a/Python/pywarpx/__init__.py b/Python/pywarpx/__init__.py index d5181a0307c..037598f4ed4 100644 --- a/Python/pywarpx/__init__.py +++ b/Python/pywarpx/__init__.py @@ -23,6 +23,7 @@ if os.path.exists(p_abs): os.add_dll_directory(p_abs) +from ._libwarpx import libwarpx from .Algo import algo from .Amr import amr from .Amrex import amrex @@ -36,10 +37,9 @@ from .Interpolation import interpolation from .Lasers import lasers from .LoadThirdParty import load_cupy -from .PSATD import psatd from .Particles import newspecies, particles +from .PSATD import psatd from .WarpX import warpx -from ._libwarpx import libwarpx # This is a circular import and must happen after the import of libwarpx from . import picmi # isort:skip diff --git a/Python/pywarpx/particle_containers.py b/Python/pywarpx/particle_containers.py index f8e687f99bd..f5330d05f84 100644 --- a/Python/pywarpx/particle_containers.py +++ b/Python/pywarpx/particle_containers.py @@ -8,8 +8,8 @@ import numpy as np -from .LoadThirdParty import load_cupy from ._libwarpx import libwarpx +from .LoadThirdParty import load_cupy class ParticleContainerWrapper(object): diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index 59ca1ec789f..911bbe13148 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -8,13 +8,12 @@ """Classes following the PICMI standard """ -from dataclasses import dataclass import os import re +from dataclasses import dataclass import numpy as np import periodictable - import picmistandard import pywarpx diff --git a/Regression/Checksum/checksum.py b/Regression/Checksum/checksum.py index 454edfc2606..727c8beb7f7 100644 --- a/Regression/Checksum/checksum.py +++ b/Regression/Checksum/checksum.py @@ -9,11 +9,11 @@ import json import sys -from benchmark import Benchmark import numpy as np +import yt +from benchmark import Benchmark from openpmd_viewer import OpenPMDTimeSeries from scipy.constants import c -import yt yt.funcs.mylog.setLevel(50) diff --git a/Regression/prepare_file_ci.py b/Regression/prepare_file_ci.py index d52d05b1139..ce7a3398b95 100644 --- a/Regression/prepare_file_ci.py +++ b/Regression/prepare_file_ci.py @@ -6,6 +6,7 @@ # License: BSD-3-Clause-LBNL import os + # This script modifies `WarpX-test.ini` (which is used for nightly builds) # and creates the file `ci-test.ini` (which is used for continuous # integration) diff --git a/Tools/Algorithms/stencil.py b/Tools/Algorithms/stencil.py index 63bb7f11c2e..dde7398daaa 100644 --- a/Tools/Algorithms/stencil.py +++ b/Tools/Algorithms/stencil.py @@ -16,9 +16,9 @@ import sys sys.path.append('../Parser/') -from input_file_parser import parse_input_file import matplotlib.pyplot as plt import numpy as np +from input_file_parser import parse_input_file from scipy.constants import c plt.style.use('tableau-colorblind10') diff --git a/Tools/LibEnsemble/run_libensemble_on_warpx.py b/Tools/LibEnsemble/run_libensemble_on_warpx.py index b86c0249d01..d825288dc56 100644 --- a/Tools/LibEnsemble/run_libensemble_on_warpx.py +++ b/Tools/LibEnsemble/run_libensemble_on_warpx.py @@ -21,20 +21,23 @@ import sys +import numpy as np + # Import libEnsemble modules from libensemble.libE import libE -import numpy as np from warpx_simf import run_warpx # Sim function from current directory if generator_type == 'random': - from libensemble.alloc_funcs.give_sim_work_first import \ - give_sim_work_first as alloc_f + from libensemble.alloc_funcs.give_sim_work_first import ( + give_sim_work_first as alloc_f, + ) from libensemble.gen_funcs.sampling import uniform_random_sample as gen_f elif generator_type == 'aposmm': import libensemble.gen_funcs libensemble.gen_funcs.rc.aposmm_optimizers = 'nlopt' - from libensemble.alloc_funcs.persistent_aposmm_alloc import \ - persistent_aposmm_alloc as alloc_f + from libensemble.alloc_funcs.persistent_aposmm_alloc import ( + persistent_aposmm_alloc as alloc_f, + ) from libensemble.gen_funcs.persistent_aposmm import aposmm as gen_f else: print("you shouldn' hit that") @@ -43,11 +46,7 @@ import all_machine_specs from libensemble import libE_logger from libensemble.executors.mpi_executor import MPIExecutor -from libensemble.tools import ( - add_unique_random_streams, - parse_args, - save_libE_output, -) +from libensemble.tools import add_unique_random_streams, parse_args, save_libE_output # Import machine-specific run parameters if machine == 'local': diff --git a/Tools/LibEnsemble/warpx_simf.py b/Tools/LibEnsemble/warpx_simf.py index 4a1f0f90dee..f310dd0f95e 100644 --- a/Tools/LibEnsemble/warpx_simf.py +++ b/Tools/LibEnsemble/warpx_simf.py @@ -1,9 +1,9 @@ import os import time +import numpy as np from libensemble.executors.executor import Executor from libensemble.message_numbers import TASK_FAILED, WORKER_DONE -import numpy as np from read_sim_output import read_sim_output from write_sim_input import write_sim_input diff --git a/Tools/PerformanceTests/run_automated.py b/Tools/PerformanceTests/run_automated.py index 73ab79b00df..f03ead05376 100644 --- a/Tools/PerformanceTests/run_automated.py +++ b/Tools/PerformanceTests/run_automated.py @@ -13,14 +13,14 @@ import sys import time +import git +import pandas as pd from functions_perftest import ( extract_dataframe, get_file_content, run_batch_nnode, store_git_hash, ) -import git -import pandas as pd # Get name of supercomputer and import configuration functions from # machine-specific file diff --git a/Tools/PostProcessing/plot_distribution_mapping.py b/Tools/PostProcessing/plot_distribution_mapping.py index 19628551f26..4b0cdfd532b 100644 --- a/Tools/PostProcessing/plot_distribution_mapping.py +++ b/Tools/PostProcessing/plot_distribution_mapping.py @@ -1,7 +1,7 @@ # Standard imports +import os from collections import defaultdict from itertools import product -import os # High-performance math import numpy as np diff --git a/Tools/PostProcessing/yt3d_mpi.py b/Tools/PostProcessing/yt3d_mpi.py index 20054756d22..655327aff3d 100644 --- a/Tools/PostProcessing/yt3d_mpi.py +++ b/Tools/PostProcessing/yt3d_mpi.py @@ -19,10 +19,10 @@ import glob -from mpi4py import MPI import numpy as np import scipy.constants as scc import yt +from mpi4py import MPI yt.funcs.mylog.setLevel(50) diff --git a/Tools/Release/updateAMReX.py b/Tools/Release/updateAMReX.py index 4f2d84e4eb7..9dfa7fbeb41 100755 --- a/Tools/Release/updateAMReX.py +++ b/Tools/Release/updateAMReX.py @@ -9,9 +9,9 @@ # when building WarpX. # import datetime -from pathlib import Path import re import sys +from pathlib import Path import requests diff --git a/Tools/Release/updatePICSAR.py b/Tools/Release/updatePICSAR.py index ce1efec0a99..7e61679d371 100755 --- a/Tools/Release/updatePICSAR.py +++ b/Tools/Release/updatePICSAR.py @@ -9,9 +9,9 @@ # when building WarpX. # import datetime -from pathlib import Path import re import sys +from pathlib import Path import requests diff --git a/Tools/Release/updatepyAMReX.py b/Tools/Release/updatepyAMReX.py index 8ba10c895e0..500781e0880 100755 --- a/Tools/Release/updatepyAMReX.py +++ b/Tools/Release/updatepyAMReX.py @@ -9,9 +9,9 @@ # when building WarpX. # import datetime -from pathlib import Path import re import sys +from pathlib import Path import requests diff --git a/pyproject.toml b/pyproject.toml index f53522b5622..5173216a0f4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -6,3 +6,6 @@ requires = [ "packaging>=23", ] build-backend = "setuptools.build_meta" + +[tool.isort] +profile = "black" From 01ee64e4a6f9ab673db299859b2935d002e93dd8 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Fri, 15 Mar 2024 21:53:09 +0100 Subject: [PATCH 10/14] fix redundant string init found with clang-tidy (#4781) --- Source/Utils/WarpXAlgorithmSelection.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Utils/WarpXAlgorithmSelection.cpp b/Source/Utils/WarpXAlgorithmSelection.cpp index b874c151f36..c58318bbd53 100644 --- a/Source/Utils/WarpXAlgorithmSelection.cpp +++ b/Source/Utils/WarpXAlgorithmSelection.cpp @@ -249,7 +249,7 @@ GetParticleBCTypeInteger( std::string BCType ){ std::string GetFieldBCTypeString( int fb_type ) { - std::string boundary_name = ""; + std::string boundary_name; for (const auto &valid_pair : FieldBCType_algo_to_int) { if ((valid_pair.second == fb_type)&&(valid_pair.first != "default")){ boundary_name = valid_pair.first; From a8f48eaf76a7d467cb39cb5ae974317545a17192 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Sat, 16 Mar 2024 00:26:18 +0100 Subject: [PATCH 11/14] Remove WarpX:: from ablastr (#4782) --- Source/ablastr/fields/VectorPoissonSolver.H | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/Source/ablastr/fields/VectorPoissonSolver.H b/Source/ablastr/fields/VectorPoissonSolver.H index 6e68d0fadb5..df2287df7f4 100644 --- a/Source/ablastr/fields/VectorPoissonSolver.H +++ b/Source/ablastr/fields/VectorPoissonSolver.H @@ -185,11 +185,10 @@ computeVectorPotential ( amrex::Vector > co relative_tolerance, absolute_tolerance ); // Synchronize the ghost cells, do halo exchange - ablastr::utils::communication::FillBoundary(*A[lev][adim], - A[lev][adim]->nGrowVect(), - WarpX::do_single_precision_comms, - geom[lev].periodicity(), - false); + ablastr::utils::communication::FillBoundary( + *A[lev][adim], A[lev][adim]->nGrowVect(), + do_single_precision_comms, + geom[lev].periodicity(), false); // needed for solving the levels by levels: // - coarser level is initial guess for finer level From 4e568adffc55c6f9b279098499e7db232aa2c540 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 15 Mar 2024 19:31:02 -0700 Subject: [PATCH 12/14] Define class function `defineAllParticleTiles` for `NamedComponentParticleContainer` (#4780) * Define function to define tiles, for a NamedComponentParticleContainer * Call parent class in WarpXParticleContainer --- Source/Particles/NamedComponentParticleContainer.H | 13 +++++++++++++ Source/Particles/WarpXParticleContainer.cpp | 5 ++++- 2 files changed, 17 insertions(+), 1 deletion(-) diff --git a/Source/Particles/NamedComponentParticleContainer.H b/Source/Particles/NamedComponentParticleContainer.H index e7a7a20fad5..983d417d9c0 100644 --- a/Source/Particles/NamedComponentParticleContainer.H +++ b/Source/Particles/NamedComponentParticleContainer.H @@ -190,6 +190,19 @@ public: } } + void defineAllParticleTiles () noexcept + { + for (int lev = 0; lev <= amrex::ParticleContainerPureSoA::finestLevel(); ++lev) + { + for (auto mfi = amrex::ParticleContainerPureSoA::MakeMFIter(lev); mfi.isValid(); ++mfi) + { + const int grid_id = mfi.index(); + const int tile_id = mfi.LocalTileIndex(); + amrex::ParticleContainerPureSoA::DefineAndReturnParticleTile(lev, grid_id, tile_id); + } + } + } + /** Return the name-to-index map for the compile-time and runtime-time real components */ [[nodiscard]] std::map getParticleComps () const noexcept { return particle_comps;} /** Return the name-to-index map for the compile-time and runtime-time integer components */ diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 6e6c0b27c85..912dd3a5a40 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -1517,6 +1517,10 @@ WarpXParticleContainer::PushX (int lev, amrex::Real dt) // without runtime component). void WarpXParticleContainer::defineAllParticleTiles () noexcept { + // Call the parent class's method + NamedComponentParticleContainer::defineAllParticleTiles(); + + // Resize the tmp_particle_data (no present in parent class) tmp_particle_data.resize(finestLevel()+1); for (int lev = 0; lev <= finestLevel(); ++lev) { @@ -1525,7 +1529,6 @@ void WarpXParticleContainer::defineAllParticleTiles () noexcept const int grid_id = mfi.index(); const int tile_id = mfi.LocalTileIndex(); tmp_particle_data[lev][std::make_pair(grid_id,tile_id)]; - DefineAndReturnParticleTile(lev, grid_id, tile_id); } } } From d22e423f6370e124dcaf61b7e3b2a4a36343df07 Mon Sep 17 00:00:00 2001 From: Ryan Sandberg Date: Sat, 16 Mar 2024 00:55:53 -0700 Subject: [PATCH 13/14] update ml workflow for pasc revisions (#4768) * update ml workflow for pasc revisions * update pasc reference: arxiv link+add andrew * update figures * Update Bibtex * refactor * clean duplication in references * add explicit reference to Zenodo archive --- Docs/source/acknowledge_us.rst | 5 ++ Docs/source/highlights.rst | 11 +++ Docs/source/refs.bib | 8 +- .../usage/workflows/ml_dataset_training.rst | 83 +++++++++++-------- .../workflows/ml_materials/create_dataset.py | 40 ++++----- .../ml_materials/run_warpx_training.py | 21 ++--- .../usage/workflows/ml_materials/train.py | 4 +- 7 files changed, 100 insertions(+), 72 deletions(-) diff --git a/Docs/source/acknowledge_us.rst b/Docs/source/acknowledge_us.rst index 79bb9fa8884..4c7a2d8137c 100644 --- a/Docs/source/acknowledge_us.rst +++ b/Docs/source/acknowledge_us.rst @@ -53,6 +53,11 @@ Prior WarpX references If your project uses a specific algorithm or component, please consider citing the respective publications in addition. +- Sandberg R T, Lehe R, Mitchell C E, Garten M, Myers A, Qiang J, Vay J-L and Huebl A. + **Synthesizing Particle-in-Cell Simulations Through Learning and GPU Computing for Hybrid Particle Accelerator Beamlines**. + Proc. of Platform for Advanced Scientific Computing (PASC'24), *submitted*, 2024. + `preprint __` + - Sandberg R T, Lehe R, Mitchell C E, Garten M, Qiang J, Vay J-L and Huebl A. **Hybrid Beamline Element ML-Training for Surrogates in the ImpactX Beam-Dynamics Code**. 14th International Particle Accelerator Conference (IPAC'23), WEPA101, 2023. diff --git a/Docs/source/highlights.rst b/Docs/source/highlights.rst index 98b0d85391e..91c99ee50ca 100644 --- a/Docs/source/highlights.rst +++ b/Docs/source/highlights.rst @@ -24,6 +24,11 @@ Scientific works in laser-plasma and beam-plasma acceleration. Phys. Rev. Research **5**, 033112, 2023 `DOI:10.1103/PhysRevResearch.5.033112 `__ +#. Sandberg R T, Lehe R, Mitchell C E, Garten M, Myers A, Qiang J, Vay J-L and Huebl A. + **Synthesizing Particle-in-Cell Simulations Through Learning and GPU Computing for Hybrid Particle Accelerator Beamlines**. + Proc. of Platform for Advanced Scientific Computing (PASC'24), *submitted*, 2024. + `preprint __` + #. Sandberg R T, Lehe R, Mitchell C E, Garten M, Qiang J, Vay J-L and Huebl A. **Hybrid Beamline Element ML-Training for Surrogates in the ImpactX Beam-Dynamics Code**. 14th International Particle Accelerator Conference (IPAC'23), WEPA101, 2023. @@ -96,6 +101,11 @@ Particle Accelerator & Beam Physics Scientific works in particle and beam modeling. +#. Sandberg R T, Lehe R, Mitchell C E, Garten M, Myers A, Qiang J, Vay J-L and Huebl A. + **Synthesizing Particle-in-Cell Simulations Through Learning and GPU Computing for Hybrid Particle Accelerator Beamlines**. + Proc. of Platform for Advanced Scientific Computing (PASC'24), *submitted*, 2024. + `preprint __` + #. Sandberg R T, Lehe R, Mitchell C E, Garten M, Qiang J, Vay J-L, Huebl A. **Hybrid Beamline Element ML-Training for Surrogates in the ImpactX Beam-Dynamics Code**. 14th International Particle Accelerator Conference (IPAC'23), WEPA101, *in print*, 2023. @@ -128,6 +138,7 @@ Microelectronics **Characterization of Transmission Lines in Microelectronic Circuits Using the ARTEMIS Solver**. IEEE Journal on Multiscale and Multiphysics Computational Techniques, vol. 8, pp. 31-39, 2023. `DOI:10.1109/JMMCT.2022.3228281 `__ + #. Kumar P, Nonaka A, Jambunathan R, Pahwa G and Salahuddin S, Yao Z. **FerroX: A GPU-accelerated, 3D Phase-Field Simulation Framework for Modeling Ferroelectric Devices**. arXiv preprint, 2022. diff --git a/Docs/source/refs.bib b/Docs/source/refs.bib index 895a6c1392b..6c7665dac84 100644 --- a/Docs/source/refs.bib +++ b/Docs/source/refs.bib @@ -207,13 +207,15 @@ @article{Roedel2010 @misc{SandbergPASC24, address = {Zuerich, Switzerland}, -author = {Ryan Sandberg and Remi Lehe and Chad E Mitchell and Marco Garten and Ji Qiang and Jean-Luc Vay and Axel Huebl}, +author = {Ryan Sandberg and Remi Lehe and Chad E Mitchell and Marco Garten and Andrew Myers and Ji Qiang and Jean-Luc Vay and Axel Huebl}, booktitle = {Proc. of PASC24}, -note = {submitted}, +note = {accepted}, series = {PASC'24 - Platform for Advanced Scientific Computing}, title = {{Synthesizing Particle-in-Cell Simulations Through Learning and GPU Computing for Hybrid Particle Accelerator Beamlines}}, venue = {Zuerich, Switzerland}, -year = {2024} +year = {2024}, +doi = {10.48550/arXiv.2402.17248}, +url = {https://arxiv.org/abs/2402.17248} } @inproceedings{SandbergIPAC23, diff --git a/Docs/source/usage/workflows/ml_dataset_training.rst b/Docs/source/usage/workflows/ml_dataset_training.rst index 6e60a318bee..450e0fe2879 100644 --- a/Docs/source/usage/workflows/ml_dataset_training.rst +++ b/Docs/source/usage/workflows/ml_dataset_training.rst @@ -14,42 +14,42 @@ For example, a simulation determined by the following input script .. literalinclude:: ml_materials/run_warpx_training.py :language: python -In this section we walk through a workflow for data processing and model training. +In this section we walk through a workflow for data processing and model training, using data from this input script as an example. +The simulation output is stored in an online `Zenodo archive `__, in the ``lab_particle_diags`` directory. +In the example scripts provided here, the data is downloaded from the Zenodo archive, properly formatted, and used to train a neural network. This workflow was developed and first presented in :cite:t:`ml-SandbergIPAC23,ml-SandbergPASC24`. - -This assumes you have an up-to-date environment with PyTorch and openPMD. +It assumes you have an up-to-date environment with PyTorch and openPMD. Data Cleaning ------------- -It is important to inspect the data for artifacts to +It is important to inspect the data for artifacts, to check that input/output data make sense. -If we plot the final phase space for beams 1-8, -the particle data is distributed in a single blob, -as shown by :numref:`fig_phase_space_beam_1` for beam 1. -This is as we expect and what is optimal for training neural networks. +If we plot the final phase space of the particle beam, +shown in :numref:`fig_unclean_phase_space`. +we see outlying particles. +Looking closer at the z-pz space, we see that some particles were not trapped in the accelerating region of the wake and have much less energy than the rest of the beam. + +.. _fig_unclean_phase_space: -.. _fig_phase_space_beam_1: +.. figure:: https://gist.githubusercontent.com/RTSandberg/649a81cc0e7926684f103729483eff90/raw/095ac2daccbcf197fa4e18a8f8505711b27e807a/unclean_stage_0.png + :alt: Plot showing the final phase space projections of a particle beam through a laser-plasma acceleration element where some beam particles were not accelerated. -.. figure:: https://user-images.githubusercontent.com/10621396/290010209-c55baf1c-dd98-4d56-a675-ad3729481eee.png - :alt: Plot showing the final phase space projections for beam 1 of the training data, for a surrogate to stage 1. + The final phase space projections of a particle beam through a laser-plasma acceleration element where some beam particles were not accelerated. - The final phase space projections for beam 1 of the training data, for a surrogate to stage 1. +To assist our neural network in learning dynamics of interest, we filter out these particles. +It is sufficient for our purposes to select particles that are not too far back, setting +``particle_selection={'z':[0.280025, None]}``. +After filtering, we can see in :numref:`fig_clean_phase_space` that the beam phase space projections are much cleaner -- this is the beam we want to train on. -.. _fig_phase_space_beam_0: +.. _fig_clean_phase_space: -.. figure:: https://user-images.githubusercontent.com/10621396/290010282-40560ac4-8509-4599-82ca-167bb1739cff.png - :alt: Plot showing the final phase space projections for beam 0 of the training data, for a surrogate to stage 0. +.. figure:: https://gist.githubusercontent.com/RTSandberg/649a81cc0e7926684f103729483eff90/raw/095ac2daccbcf197fa4e18a8f8505711b27e807a/clean_stage_0.png + :alt: Plot showing the final phase space projections of a particle beam through a laser-plasma acceleration element after filtering out outlying particles. - The final phase space projections for beam 0 of the training data, for a surrogate to stage 0 + The final phase space projections of a particle beam through a laser-plasma acceleration element after filtering out outlying particles. -On the other hand, the final phase space for beam 0, shown in :numref:`fig_phase_space_beam_1`, -has a halo of outlying particles. -Looking closer at the z-pz space, we see that some particles got caught in a decelerating -region of the wake, have slipped back and are much slower than the rest of the beam. -To assist our neural network in learning dynamics of interest, we filter out these particles. -It is sufficient for our purposes to select particles that are not too far back, setting -``particle_selection={'z':[0.28002, None]}``. Then a particle tracker is set up to make sure +A particle tracker is set up to make sure we consistently filter out these particles from both the initial and final data. .. literalinclude:: ml_materials/create_dataset.py @@ -58,6 +58,9 @@ we consistently filter out these particles from both the initial and final data. :start-after: # Manual: Particle tracking START :end-before: # Manual: Particle tracking END +This data cleaning ensures that the particle data is distributed in a single blob, +as is optimal for training neural networks. + Create Normalized Dataset ------------------------- @@ -119,7 +122,12 @@ This data are converted to an :math:`N\times 6` numpy array and then to a PyTorc Save Normalizations and Normalized Data ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -With the data properly normalized, it and the normalizations are saved to file for +The data is split into training and testing subsets. +We take most of the data (70%) for training, meaning that data is used to update +the neural network parameters. +The testing data is reserved to determine how well the neural network generalizes; +that is, how well the neural network performs on data that wasn't used to update the neural network parameters. +With the data split and properly normalized, it and the normalizations are saved to file for use in training and inference. .. literalinclude:: ml_materials/create_dataset.py @@ -131,13 +139,22 @@ use in training and inference. Neural Network Structure ------------------------ -It was found in :cite:t:`ml-SandbergPASC24` that reasonable surrogate models are obtained with -shallow feedforward neural networks consisting of fewer than 10 hidden layers and -just under 1000 nodes per layer. +It was found in :cite:t:`ml-SandbergPASC24` that a reasonable surrogate model is obtained with +shallow feedforward neural networks consisting of about 5 hidden layers and 700-900 nodes per layer. The example shown here uses 3 hidden layers and 20 nodes per layer and is trained for 10 epochs. +Some utility functions for creating neural networks are provided in the script below. +These are mostly convenience wrappers and utilities for working with `PyTorch `__ neural network objects. +This script is imported in the training scripts shown later. + +.. dropdown:: Python neural network class definitions + :color: light + :icon: info + :animate: fade-in-slide-down + .. literalinclude:: ml_materials/neural_network_classes.py + :language: python3 Train and Save Neural Network ----------------------------- @@ -188,8 +205,8 @@ which is later divided by the size of the dataset in the training loop. :start-after: # Manual: Test function START :end-before: # Manual: Test function END -Train Loop -^^^^^^^^^^ +Training Loop +^^^^^^^^^^^^^ The full training loop performs ``n_epochs`` number of iterations. At each iteration the training and testing functions are called, @@ -228,14 +245,14 @@ When the test-loss starts to trend flat or even upward, the neural network is no .. _fig_train_test_loss: -.. figure:: https://user-images.githubusercontent.com/10621396/290010428-f83725ab-a08f-494c-b075-314b0d26cb9a.png +.. figure:: https://gist.githubusercontent.com/RTSandberg/649a81cc0e7926684f103729483eff90/raw/095ac2daccbcf197fa4e18a8f8505711b27e807a/beam_stage_0_training_testing_error.png :alt: Plot of training and testing loss curves versus number of training epochs. Training (in blue) and testing (in green) loss curves versus number of training epochs. .. _fig_train_evaluation: -.. figure:: https://user-images.githubusercontent.com/10621396/290010486-4a3541e7-e0be-4cf1-b33b-57d5e5985196.png +.. figure:: https://gist.githubusercontent.com/RTSandberg/649a81cc0e7926684f103729483eff90/raw/095ac2daccbcf197fa4e18a8f8505711b27e807a/beam_stage_0_model_evaluation.png :alt: Plot comparing model prediction with simulation output. A comparison of model prediction (yellow-red dots, colored by mean-squared error) with simulation output (black dots). @@ -243,7 +260,7 @@ When the test-loss starts to trend flat or even upward, the neural network is no A visual inspection of the model prediction can be seen in :numref:`fig_train_evaluation`. This plot compares the model prediction, with dots colored by mean-square error, on the testing data with the actual simulation output in black. The model obtained with the hyperparameters chosen here trains quickly but is not very accurate. -A more accurate model is obtained with 5 hidden layers and 800 nodes per layer, +A more accurate model is obtained with 5 hidden layers and 900 nodes per layer, as discussed in :cite:t:`ml-SandbergPASC24`. These figures can be generated with the following Python script. @@ -261,7 +278,7 @@ Surrogate Usage in Accelerator Physics ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ A neural network such as the one we trained here can be incorporated in other BLAST codes. -`Consider the example using neural networks in ImpactX `__. +Consider this `example using neural network surrogates of WarpX simulations in ImpactX `__. .. bibliography:: :keyprefix: ml- diff --git a/Docs/source/usage/workflows/ml_materials/create_dataset.py b/Docs/source/usage/workflows/ml_materials/create_dataset.py index 1dd16f2f95e..ecd182c8802 100644 --- a/Docs/source/usage/workflows/ml_materials/create_dataset.py +++ b/Docs/source/usage/workflows/ml_materials/create_dataset.py @@ -7,17 +7,16 @@ # Authors: Ryan Sandberg # License: BSD-3-Clause-LBNL # + import os -import tarfile +import zipfile from urllib import request import numpy as np - -c = 2.998e8 - import torch from openpmd_viewer import OpenPMDTimeSeries, ParticleTracker +c = 2.998e8 ############### def sanitize_dir_strings(*dir_strings): @@ -31,14 +30,15 @@ def sanitize_dir_strings(*dir_strings): return dir_strings def download_and_unzip(url, data_dir): - request.urlretrieve(url, data_dir) - with tarfile.open(data_dir) as tar_dataset: - tar_dataset.extractall() + request.urlretrieve(url, data_dir) + with zipfile.ZipFile(data_dir, 'r') as zip_dataset: + zip_dataset.extractall() def create_source_target_data(data_dir, species, source_index=0, target_index=-1, + survivor_select_index=-1, particle_selection=None ): """Create dataset from openPMD files @@ -67,7 +67,7 @@ def create_source_target_data(data_dir, relevant_times = [ts.t[source_index], ts.t[target_index]] # Manual: Particle tracking START - iteration = ts.iterations[target_index] + iteration = ts.iterations[survivor_select_index] pt = ParticleTracker( ts, species=species, iteration=iteration, @@ -116,7 +116,6 @@ def create_source_target_data(data_dir, return source_data, source_means, source_stds, target_data, target_means, target_stds, relevant_times - def save_warpx_surrogate_data(dataset_fullpath_filename, diag_dir, species, @@ -128,12 +127,14 @@ def save_warpx_surrogate_data(dataset_fullpath_filename, particle_selection=None ): - source_target_data = create_source_target_data(data_dir=diag_dir, - species=species, - source_index=source_index, - target_index=target_index, - particle_selection=particle_selection - ) + source_target_data = create_source_target_data( + data_dir=diag_dir, + species=species, + source_index=source_index, + target_index=target_index, + survivor_select_index=survivor_select_index, + particle_selection=particle_selection + ) source_data, source_means, source_stds, target_data, target_means, target_stds, times = source_target_data # Manual: Save dataset START @@ -161,11 +162,10 @@ def save_warpx_surrogate_data(dataset_fullpath_filename, ######## end utility functions ############# ######## start dataset creation ############ -data_url = "https://zenodo.org/records/10368972/files/ml_example_training.tar.gz?download=1" -download_and_unzip(data_url, "training_dataset.tar.gz") +data_url = "https://zenodo.org/records/10810754/files/lab_particle_diags.zip?download=1" +download_and_unzip(data_url, "lab_particle_diags.zip") data_dir = "lab_particle_diags/lab_particle_diags/" - # create data set source_index = 0 @@ -178,7 +178,7 @@ def save_warpx_surrogate_data(dataset_fullpath_filename, # improve stage 0 dataset stage_i = 0 -select = {'z':[0.28002, None]} +select = {'z':[0.280025, None]} species = f'beam_stage_{stage_i}' dataset_filename = f'dataset_{species}.pt' dataset_file = 'datasets/' + dataset_filename @@ -193,7 +193,7 @@ def save_warpx_surrogate_data(dataset_fullpath_filename, particle_selection=select ) -for stage_i in range(1,9): +for stage_i in range(1,15): species = f'beam_stage_{stage_i}' dataset_filename = f'dataset_{species}.pt' dataset_file = 'datasets/' + dataset_filename diff --git a/Docs/source/usage/workflows/ml_materials/run_warpx_training.py b/Docs/source/usage/workflows/ml_materials/run_warpx_training.py index 772fe574504..23fc49e7431 100644 --- a/Docs/source/usage/workflows/ml_materials/run_warpx_training.py +++ b/Docs/source/usage/workflows/ml_materials/run_warpx_training.py @@ -1,11 +1,4 @@ #!/usr/bin/env python3 -# -# Copyright 2022-2023 WarpX contributors -# Authors: WarpX team -# License: BSD-3-Clause-LBNL -# -# -*- coding: utf-8 -*- - import math import numpy as np @@ -21,7 +14,7 @@ # Number of cells dim = '3' nx = ny = 128 -nz = 8832 +nz = 35328 #17664 #8832 if dim == 'rz': nr = nx//2 @@ -33,9 +26,9 @@ # Number of processes for static load balancing # Check with your submit script -num_procs = [1, 1, 16*4] +num_procs = [1, 1, 64*4] if dim == 'rz': - num_procs = [1, 16] + num_procs = [1, 64] # Number of time steps gamma_boost = 60. @@ -73,7 +66,7 @@ # plasma region plasma_rlim = 100.e-6 -N_stage = 9 +N_stage = 15 L_plasma_bulk = 0.28 L_ramp = 1.e-9 L_ramp_up = L_ramp @@ -141,12 +134,12 @@ def get_species_of_accelerator_stage(stage_idx, stage_zmin, stage_zmax, N_beam_particles = int(1e6) beam_centroid_z = -107.e-6 beam_rms_z = 2.e-6 -#beam_gammas = [2000 + 13000 * i_stage for i_stage in range(N_stage)] -beam_gammas = [1957, 15188, 28432, 41678, 54926, 68174, 81423,94672, 107922,121171] # From 3D run +beam_gammas = [1960 + 13246 * i_stage for i_stage in range(N_stage)] +#beam_gammas = [1957, 15188, 28432, 41678, 54926, 68174, 81423,94672, 107922,121171] # From 3D run beams = [] for i_stage in range(N_stage): beam_gamma = beam_gammas[i_stage] - sigma_gamma = 0.10 * beam_gamma + sigma_gamma = 0.06 * beam_gamma gaussian_distribution = picmi.GaussianBunchDistribution( n_physical_particles= abs(beam_charge) / q_e, rms_bunch_size=[2.e-6, 2.e-6, beam_rms_z], diff --git a/Docs/source/usage/workflows/ml_materials/train.py b/Docs/source/usage/workflows/ml_materials/train.py index dc62819061c..4de11b9c99e 100644 --- a/Docs/source/usage/workflows/ml_materials/train.py +++ b/Docs/source/usage/workflows/ml_materials/train.py @@ -20,8 +20,8 @@ stage_i = 0 species = f'beam_stage_{stage_i}' source_index = 0 -target_index = 4 -survivor_select_index = 4 +target_index = 1 +survivor_select_index = 1 data_dim = 6 n_in = data_dim From c6675cddc2d698c37ac4af235f8e5a9dbbf5527e Mon Sep 17 00:00:00 2001 From: RevathiJambunathan Date: Mon, 18 Mar 2024 13:20:12 -0700 Subject: [PATCH 14/14] smooth tanh weighting for fp and cp interpolation in buffer region --- Source/Parallelization/WarpXRegrid.cpp | 1 + Source/Particles/PhysicalParticleContainer.H | 22 +++ .../Particles/PhysicalParticleContainer.cpp | 142 ++++++++++++++++-- Source/WarpX.H | 36 +++++ Source/WarpX.cpp | 88 ++++++++++- 5 files changed, 275 insertions(+), 14 deletions(-) diff --git a/Source/Parallelization/WarpXRegrid.cpp b/Source/Parallelization/WarpXRegrid.cpp index dcea3cab58a..1bd05a2d115 100644 --- a/Source/Parallelization/WarpXRegrid.cpp +++ b/Source/Parallelization/WarpXRegrid.cpp @@ -334,6 +334,7 @@ WarpX::RemakeLevel (int lev, Real /*time*/, const BoxArray& ba, const Distributi // we can avoid redistributing these since we immediately re-build the values via BuildBufferMasks() RemakeMultiFab(current_buffer_masks[lev], dm, false ,lev); RemakeMultiFab(gather_buffer_masks[lev], dm, false ,lev); + RemakeMultiFab(interp_weight_gbuffer[lev], dm, false, lev); if (current_buffer_masks[lev] || gather_buffer_masks[lev]) { BuildBufferMasks(); diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 286d0675da6..73e68fa9a9b 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -320,6 +320,28 @@ public: amrex::FArrayBox const * & ez_ptr, amrex::FArrayBox const * & bx_ptr, amrex::FArrayBox const * & by_ptr, amrex::FArrayBox const * & bz_ptr); + void InterpolateFieldsInGatherBuffer (const amrex::Box& box, + amrex::FArrayBox& bufferEx, amrex::FArrayBox& bufferEy, amrex::FArrayBox& bufferEz, + amrex::FArrayBox& bufferBx, amrex::FArrayBox& bufferBy, amrex::FArrayBox& bufferBz, + amrex::Elixir& buf_exeli, amrex::Elixir& buf_eyeli, amrex::Elixir& buf_ezeli, + amrex::Elixir& buf_bxeli, amrex::Elixir& buf_byeli, amrex::Elixir& buf_bzeli, + amrex::FArrayBox const *& ex_ptr, amrex::FArrayBox const *& ey_ptr, amrex::FArrayBox const* & ez_ptr, + amrex::FArrayBox const *& bx_ptr, amrex::FArrayBox const *& by_ptr, amrex::FArrayBox const* & bz_ptr, + amrex::FArrayBox const *& cex_ptr, amrex::FArrayBox const *& cey_ptr, amrex::FArrayBox const* & cez_ptr, + amrex::FArrayBox const *& cbx_ptr, amrex::FArrayBox const *& cby_ptr, amrex::FArrayBox const* & cbz_ptr, + const amrex::IArrayBox& gather_mask, const amrex::FArrayBox& weight_gbuffer, + const amrex::IntVect ref_ratio); + + void WeightedSumInGatherBuffer (const amrex::Box& tmp_xbox, const amrex::Box& tmp_ybox, const amrex::Box& tmp_zbox, + amrex::Array4 const & xbuf_field_arr, amrex::Array4 const& ybuf_field_arr, + amrex::Array4 const & zbuf_field_arr, + amrex::Array4 cfx_arr, amrex::Array4 cfy_arr, + amrex::Array4 cfz_arr, + amrex::Array4 fx_arr, amrex::Array4 fy_arr, + amrex::Array4 fz_arr, + amrex::Array4 gm_arr, amrex::Array4 wt_arr, + const amrex::IntVect ref_ratio); + /** * \brief This function determines if resampling should be done for the current species, and * if so, performs the resampling. diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 160eac0d19c..f55ced70d13 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -2010,12 +2010,11 @@ PhysicalParticleContainer::Evolve (int lev, WARPX_PROFILE_VAR_NS("PhysicalParticleContainer::Evolve::GatherAndPush", blp_fg); BL_ASSERT(OnSameGrids(lev,jx)); - amrex::LayoutData* cost = WarpX::getCosts(lev); const iMultiFab* current_masks = WarpX::CurrentBufferMasks(lev); const iMultiFab* gather_masks = WarpX::GatherBufferMasks(lev); - + const amrex::MultiFab* weight_gbuffer = WarpX::GetInstance().getInterpWeightGBuffer(lev); const bool has_buffer = cEx || cjx; if (m_do_back_transformed_particles) @@ -2044,7 +2043,8 @@ PhysicalParticleContainer::Evolve (int lev, FArrayBox filtered_Ex, filtered_Ey, filtered_Ez; FArrayBox filtered_Bx, filtered_By, filtered_Bz; - + FArrayBox bufferEx, bufferEy, bufferEz; + FArrayBox bufferBx, bufferBy, bufferBz; for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers) @@ -2073,6 +2073,7 @@ PhysicalParticleContainer::Evolve (int lev, FArrayBox const* bzfab = &Bz[pti]; Elixir exeli, eyeli, ezeli, bxeli, byeli, bzeli; + Elixir buf_exeli, buf_eyeli, buf_ezeli, buf_bxeli, buf_byeli, buf_bzeli; if (WarpX::use_fdtd_nci_corr) { @@ -2114,7 +2115,7 @@ PhysicalParticleContainer::Evolve (int lev, DepositCharge(pti, wp, ion_lev, rho, 0, 0, np_current, thread_num, lev, lev); - if (has_buffer){ + if ((np-np_current)> 0 ){ DepositCharge(pti, wp, ion_lev, crho, 0, np_current, np-np_current, thread_num, lev, lev-1); } @@ -2157,8 +2158,19 @@ PhysicalParticleContainer::Evolve (int lev, FArrayBox const* cbyfab = &(*cBy)[pti]; FArrayBox const* cbzfab = &(*cBz)[pti]; + //// buffer box (Both 2D and 3D) + InterpolateFieldsInGatherBuffer(box, bufferEx, bufferEy, bufferEz, + bufferBx, bufferBy, bufferBz, + buf_exeli, buf_eyeli, buf_ezeli, + buf_bxeli, buf_byeli, buf_bzeli, + exfab, eyfab, ezfab, bxfab, byfab, bzfab, + cexfab, ceyfab, cezfab, cbxfab, cbyfab, cbzfab, + (*gather_masks)[pti], (*weight_gbuffer)[pti], + ref_ratio); + if (WarpX::use_fdtd_nci_corr) { + // should this be bufferEFields* // Filter arrays (*cEx)[pti], store the result in // filtered_Ex and update pointer cexfab so that it // points to filtered_Ex (and do the same for all @@ -2170,15 +2182,17 @@ PhysicalParticleContainer::Evolve (int lev, (*cBx)[pti], (*cBy)[pti], (*cBz)[pti], cexfab, ceyfab, cezfab, cbxfab, cbyfab, cbzfab); } - // Field gather and push for particles in gather buffers e_is_nodal = cEx->is_nodal() and cEy->is_nodal() and cEz->is_nodal(); if (push_type == PushType::Explicit) { - PushPX(pti, cexfab, ceyfab, cezfab, - cbxfab, cbyfab, cbzfab, + //PushPX(pti, cexfab, ceyfab, cezfab, + // cbxfab, cbyfab, cbzfab, + PushPX(pti, &bufferEx, &bufferEy, &bufferEz, + &bufferBx, &bufferBy, &bufferBz, cEx->nGrowVect(), e_is_nodal, nfine_gather, np-nfine_gather, - lev, lev-1, dt, ScaleFields(false), a_dt_type); + lev, gather_lev, dt, ScaleFields(false), a_dt_type); + //lev, lev-1, dt, ScaleFields(false), a_dt_type); } else if (push_type == PushType::Implicit) { ImplicitPushXP(pti, cexfab, ceyfab, cezfab, cbxfab, cbyfab, cbzfab, @@ -2187,7 +2201,6 @@ PhysicalParticleContainer::Evolve (int lev, lev, lev-1, dt, ScaleFields(false), a_dt_type); } } - WARPX_PROFILE_VAR_STOP(blp_fg); // Current Deposition @@ -2204,7 +2217,7 @@ PhysicalParticleContainer::Evolve (int lev, 0, np_current, thread_num, lev, lev, dt, relative_time, push_type); - if (has_buffer) + if ((np-np_current)>0) { // Deposit in buffers DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, cjx, cjy, cjz, @@ -2224,7 +2237,7 @@ PhysicalParticleContainer::Evolve (int lev, DepositCharge(pti, wp, ion_lev, rho, 1, 0, np_current, thread_num, lev, lev); - if (has_buffer){ + if ((np-np_current)>0 ){ DepositCharge(pti, wp, ion_lev, crho, 1, np_current, np-np_current, thread_num, lev, lev-1); } @@ -2327,6 +2340,113 @@ PhysicalParticleContainer::applyNCIFilter ( #endif } +void +PhysicalParticleContainer::InterpolateFieldsInGatherBuffer (const Box& box, + amrex::FArrayBox& bufferEx, amrex::FArrayBox& bufferEy, amrex::FArrayBox& bufferEz, + amrex::FArrayBox& bufferBx, amrex::FArrayBox& bufferBy, amrex::FArrayBox& bufferBz, + Elixir& buf_exeli, Elixir& buf_eyeli, Elixir& buf_ezeli, + Elixir& buf_bxeli, Elixir& buf_byeli, Elixir& buf_bzeli, + FArrayBox const *& ex_ptr, FArrayBox const *& ey_ptr, FArrayBox const* & ez_ptr, + FArrayBox const *& bx_ptr, FArrayBox const *& by_ptr, FArrayBox const* & bz_ptr, + FArrayBox const *& cex_ptr, FArrayBox const *& cey_ptr, FArrayBox const* & cez_ptr, + FArrayBox const *& cbx_ptr, FArrayBox const *& cby_ptr, FArrayBox const* & cbz_ptr, + const IArrayBox& gather_mask, const FArrayBox& weight_gbuffer, + const amrex::IntVect ref_ratio) +{ + amrex::Box tmp_exbox = amrex::convert(box,ex_ptr->box().ixType()); + amrex::Box tmp_eybox = amrex::convert(box,ey_ptr->box().ixType()); + amrex::Box tmp_ezbox = amrex::convert(box,ez_ptr->box().ixType()); + + bufferEx.resize(tmp_exbox); + buf_exeli = bufferEx.elixir(); + bufferEy.resize(tmp_eybox); + buf_eyeli = bufferEy.elixir(); + bufferEz.resize(tmp_ezbox); + buf_ezeli = bufferEz.elixir(); + WeightedSumInGatherBuffer(tmp_exbox, tmp_eybox, tmp_ezbox, + bufferEx.array(), bufferEy.array(), bufferEz.array(), + cex_ptr->array(), cey_ptr->array(), cez_ptr->array(), + ex_ptr->array(), ey_ptr->array(), ez_ptr->array(), + gather_mask.array(), weight_gbuffer.array(), ref_ratio); + + amrex::Box tmp_bxbox = amrex::convert(box,bx_ptr->box().ixType()); + amrex::Box tmp_bybox = amrex::convert(box,by_ptr->box().ixType()); + amrex::Box tmp_bzbox = amrex::convert(box,bz_ptr->box().ixType()); + bufferBx.resize(tmp_bxbox); + buf_bxeli = bufferBx.elixir(); + bufferBy.resize(tmp_bybox); + buf_byeli = bufferBy.elixir(); + bufferBz.resize(tmp_bzbox); + buf_bzeli = bufferBz.elixir(); + + WeightedSumInGatherBuffer(tmp_bxbox, tmp_bybox, tmp_bzbox, + bufferBx.array(), bufferBy.array(), bufferBz.array(), + cbx_ptr->array(), cby_ptr->array(), cbz_ptr->array(), + bx_ptr->array(), by_ptr->array(), bz_ptr->array(), + gather_mask.array(), weight_gbuffer.array(), ref_ratio); +} + +void +PhysicalParticleContainer::WeightedSumInGatherBuffer ( + const amrex::Box& tmp_xbox, const amrex::Box& tmp_ybox, const amrex::Box& tmp_zbox, + amrex::Array4 const & xbuf_field_arr, amrex::Array4 const& ybuf_field_arr, + amrex::Array4 const & zbuf_field_arr, + amrex::Array4 cfx_arr, amrex::Array4 cfy_arr, + amrex::Array4 cfz_arr, + amrex::Array4 fx_arr, amrex::Array4 fy_arr, + amrex::Array4 fz_arr, + amrex::Array4 gm_arr, amrex::Array4 wt_arr, + const amrex::IntVect ref_ratio) +{ + amrex::ParallelFor( tmp_xbox, tmp_ybox, tmp_zbox, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + int ii = amrex::coarsen(i,ref_ratio[0]); + int jj = j; + int kk = k; +#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + jj = amrex::coarsen(j,ref_ratio[1]); +#elif defined(WARPX_DIM_3D) + kk = amrex::coarsen(k,ref_ratio[2]); +#endif + if (gm_arr(i,j,k) == 0) { + xbuf_field_arr(i,j,k) = wt_arr(i,j,k)*fx_arr(i,j,k) + (1._rt-wt_arr(i,j,k))*cfx_arr(ii,jj,kk); + } else { + xbuf_field_arr(i,j,k) = fx_arr(i,j,k); + } + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + int ii = amrex::coarsen(i,ref_ratio[0]); + int jj = j; + int kk = k; +#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + jj = amrex::coarsen(j,ref_ratio[1]); +#elif defined(WARPX_DIM_3D) + kk = amrex::coarsen(k,ref_ratio[2]); +#endif + if (gm_arr(i,j,k) == 0) { + ybuf_field_arr(i,j,k) = wt_arr(i,j,k)*fy_arr(i,j,k) + (1._rt-wt_arr(i,j,k))*cfy_arr(ii,jj,kk); + } else { + ybuf_field_arr(i,j,k) = fy_arr(i,j,k); + } + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + int ii = amrex::coarsen(i,ref_ratio[0]); + int jj = j; + int kk = k; +#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + jj = amrex::coarsen(j,ref_ratio[1]); +#elif defined(WARPX_DIM_3D) + kk = amrex::coarsen(k,ref_ratio[2]); +#endif + if (gm_arr(i,j,k) == 0) { + zbuf_field_arr(i,j,k) = wt_arr(i,j,k)*fz_arr(i,j,k) + (1._rt-wt_arr(i,j,k))*cfz_arr(ii,jj,kk); + } else { + zbuf_field_arr(i,j,k) = fz_arr(i,j,k); + } + } + ); +} + // Loop over all particles in the particle container and // split particles tagged with p.id()=DoSplitParticleID void diff --git a/Source/WarpX.H b/Source/WarpX.H index 6526171ff67..faea121f5ee 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -355,6 +355,32 @@ public: //! #n_field_gather_buffer cells of the edge of the patch, will gather the fields //! from the lower refinement level instead of the refinement patch itself static int n_field_gather_buffer; + + /** + * If true, particles located inside refinement patch but within + * #n_field_gather_buffer cells of the patch edge will gather fields + * that are a weighted sum of the fields from lower refinement + * level and the refinement patch. By default, no interpolation will be done, + * and the fields will be gathered only from the lower refinement level. + * The option can be used to allow for a smooth transition from the fine + * and coarse solutions with a tanh profile, such that, the value + * of the gather fields is close to the solution on the refinement patch + * near the gather buffer edge, and close to the solution on the lower + * refinement level close to the edge of the patch. + * #tanh_midpoint_gather_buffer can be used to set the tanh profile + */ + static bool do_fieldinterp_gather_buffer; + + /** + * With mesh-refinement, finite number of gather buffer cells, and + * do_fieldinterp_gather_buffer, the particles on the fine patch + * gather a weighted sum of fields from fine patch and coarse patch. + * A tanh profile is used to set the smoothing kernel + * #tanh_midpoint_gather_buffer must be a fraction of the buffer width + * measured from coarsefine interface, where the tanh profile is 0.5 + * The default value is 0.5 + */ + static amrex::Real tanh_midpoint_gather_buffer; //! With mesh refinement, particles located inside a refinement patch, but within //! #n_current_deposition_buffer cells of the edge of the patch, will deposit their charge //! and current onto the lower refinement level instead of the refinement patch itself @@ -878,6 +904,10 @@ public: static std::array CellSize (int lev); static amrex::RealBox getRealBox(const amrex::Box& bx, int lev); + [[nodiscard]] const amrex::MultiFab* getInterpWeightGBuffer (int lev) const + { + return interp_weight_gbuffer[lev].get(); + } /** * \brief Return the lower corner of the box in real units. * \param bx The box @@ -1392,6 +1422,7 @@ private: return gather_buffer_masks[lev].get(); } + /** * \brief Re-orders the Fornberg coefficients so that they can be used more conveniently for * finite-order centering operations. For example, for finite-order centering of order 6, @@ -1565,6 +1596,11 @@ private: amrex::Vector, 3 > > Bfield_cax; amrex::Vector > current_buffer_masks; amrex::Vector > gather_buffer_masks; + /** Multifab that stores weights for interpolating fine and coarse solutions + * in the buffer gather region, such that, the weight for fine solution is + * 0 near the coarse-fine interface, and 1 as the buffer gather region terminates in the fine patch + */ + amrex::Vector > interp_weight_gbuffer; // If charge/current deposition buffers are used amrex::Vector, 3 > > current_buf; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index c929718b043..57d287e234d 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -207,6 +207,8 @@ IntVect WarpX::filter_npass_each_dir(1); int WarpX::n_field_gather_buffer = -1; int WarpX::n_current_deposition_buffer = -1; +bool WarpX::do_fieldinterp_gather_buffer = false; +amrex::Real WarpX::tanh_midpoint_gather_buffer = 0.5; short WarpX::grid_type; amrex::IntVect m_rho_nodal_flag; @@ -399,6 +401,7 @@ WarpX::WarpX () gather_buffer_masks.resize(nlevs_max); current_buf.resize(nlevs_max); charge_buf.resize(nlevs_max); + interp_weight_gbuffer.resize(nlevs_max); pml.resize(nlevs_max); #if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD) @@ -888,6 +891,10 @@ WarpX::ReadParameters () pp_warpx, "n_field_gather_buffer", n_field_gather_buffer); utils::parser::queryWithParser( pp_warpx, "n_current_deposition_buffer", n_current_deposition_buffer); + utils::parser::queryWithParser( + pp_warpx, "do_fieldinterp_gather_buffer", do_fieldinterp_gather_buffer); + utils::parser::queryWithParser( + pp_warpx, "tanh_midpoint_gather_buffer", tanh_midpoint_gather_buffer); amrex::Real quantum_xi_tmp; const auto quantum_xi_is_specified = @@ -2067,6 +2074,7 @@ WarpX::ClearLevel (int lev) current_buffer_masks[lev].reset(); gather_buffer_masks[lev].reset(); + interp_weight_gbuffer[lev].reset(); F_fp [lev].reset(); G_fp [lev].reset(); @@ -2688,20 +2696,23 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm AllocInitMultiFab(Efield_cax[lev][1], amrex::convert(cba,Ey_nodal_flag),dm,ncomps,ngEB,lev, "Efield_cax[y]"); AllocInitMultiFab(Efield_cax[lev][2], amrex::convert(cba,Ez_nodal_flag),dm,ncomps,ngEB,lev, "Efield_cax[z]"); } - - AllocInitMultiFab(gather_buffer_masks[lev], ba, dm, ncomps, amrex::IntVect(1), lev, "gather_buffer_masks"); + amrex::Print() << " ba for allocating gahter buffer mask " << "\n"; // Gather buffer masks have 1 ghost cell, because of the fact // that particles may move by more than one cell when using subcycling. + AllocInitMultiFab(gather_buffer_masks[lev], ba, dm, ncomps, amrex::IntVect(1), lev, "gather_buffer_masks"); + AllocInitMultiFab(interp_weight_gbuffer[lev], amrex::convert(ba, IntVect::TheNodeVector()), dm, ncomps, ngEB, lev, "interp_weight_gbuffer", 0.0_rt); } if (n_current_deposition_buffer > 0) { + amrex::Print() << " current dep buffer : " << n_current_deposition_buffer << "\n"; AllocInitMultiFab(current_buf[lev][0], amrex::convert(cba,jx_nodal_flag),dm,ncomps,ngJ,lev, "current_buf[x]"); AllocInitMultiFab(current_buf[lev][1], amrex::convert(cba,jy_nodal_flag),dm,ncomps,ngJ,lev, "current_buf[y]"); AllocInitMultiFab(current_buf[lev][2], amrex::convert(cba,jz_nodal_flag),dm,ncomps,ngJ,lev, "current_buf[z]"); if (rho_cp[lev]) { AllocInitMultiFab(charge_buf[lev], amrex::convert(cba,rho_nodal_flag),dm,2*ncomps,ngRho,lev, "charge_buf"); } - AllocInitMultiFab(current_buffer_masks[lev], ba, dm, ncomps, amrex::IntVect(1), lev, "current_buffer_masks"); + amrex::Print() << " allocate current buffer mask \n"; + AllocInitMultiFab(current_buffer_masks[lev], ba, dm, ncomps, ngJ, lev, "current_buffer_masks"); // Current buffer masks have 1 ghost cell, because of the fact // that particles may move by more than one cell when using subcycling. } @@ -3063,6 +3074,8 @@ WarpX::getLoadBalanceEfficiency (const int lev) void WarpX::BuildBufferMasks () { + bool do_interpolate = WarpX::do_fieldinterp_gather_buffer; + amrex::Real tanh_midpoint = WarpX::tanh_midpoint_gather_buffer; for (int lev = 1; lev <= maxLevel(); ++lev) { for (int ipass = 0; ipass < 2; ++ipass) @@ -3088,6 +3101,75 @@ WarpX::BuildBufferMasks () const Box tbx = mfi.growntilebox(); BuildBufferMasksInBox( tbx, (*bmasks)[mfi], tmp[mfi], ngbuffer ); } + if (ipass==0) continue; + amrex::MultiFab* weight_gbuffer = interp_weight_gbuffer[lev].get(); + // Using tmp to also set weights in the interp_weight_gbuffer multifab + for (MFIter mfi(*weight_gbuffer, true); mfi.isValid(); ++mfi) + { + const Box& tbx = mfi.tilebox(IntVect::TheNodeVector(),weight_gbuffer->nGrowVect()); + auto const& gmsk = tmp[mfi].const_array(); + auto const& bmsk = (*bmasks)[mfi].array(); + auto const& wtmsk = (*weight_gbuffer)[mfi].array(); + amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + wtmsk(i,j,k) = 0._rt; + if (bmsk(i,j,k) == 0 && do_interpolate) { + if(gmsk(i,j,k)==0) { + wtmsk(i,j,k) = 0.; + return; + } + //for (int ii = i-1; ii >= i-ngbuffer; --ii) { + // if (gmsk(ii,j,k)==0) { + // amrex::Real arg = (static_cast(i-ii)-ngbuffer*tanh_midpoint) + // / ((1.-tanh_midpoint)*(ngbuffer/3.)); + // wtmsk(i,j,k) = std::tanh(arg)*0.5 + 0.5; + // amrex::Print() << " i edge wt is " << wtmsk(i,j,k) << "\n"; + // return; + // } + //} + //for (int ii = i+1; ii <= i+ngbuffer; ++ii) { + // if (gmsk(ii,j,k)==0) { + // amrex::Real arg = (static_cast(ii-i)-ngbuffer*tanh_midpoint) + // / ((1.-tanh_midpoint)*(ngbuffer/3.)); + // wtmsk(i,j,k) = std::tanh(arg)*0.5+0.5; + // amrex::Print() << " wt is " << wtmsk(i,j,k) << "\n"; + // return; + // } + //} + for (int jj = j-1; jj >= j-ngbuffer; --jj) { + if (gmsk(i,jj,k)==0) { + amrex::Real arg = (static_cast(j-jj)-ngbuffer*tanh_midpoint) + / ((1.-tanh_midpoint)*(ngbuffer/3.)); + wtmsk(i,j,k) = std::tanh(arg)*0.5 + 0.5; + return; + } + } + for (int jj = j+1; jj <= j+ngbuffer; ++jj) { + if (gmsk(i,jj,k)==0) { + amrex::Real arg = (static_cast(jj - j)-ngbuffer*tanh_midpoint) + / ((1.-tanh_midpoint)*(ngbuffer/3.)); + wtmsk(i,j,k) = std::tanh(arg)*0.5+0.5; + return; + } + } + //for (int kk = k-1; kk >= k-ngbuffer; --kk) { + // if (gmsk(i,j,kk)==0) { + // amrex::Real arg = (static_cast(k-kk)-ngbuffer*tanh_midpoint) + // / ((1.-tanh_midpoint)*(ngbuffer/3.)); + // wtmsk(i,j,k) = std::tanh(arg)*0.5+0.5; + // return; + // } + //} + //for (int kk = k+1; kk <= k+ngbuffer; ++kk) { + // if (gmsk(i,j,kk)==0) { + // amrex::Real arg = (static_cast(kk-k)-ngbuffer*tanh_midpoint) + // / ((1.-tanh_midpoint)*(ngbuffer/3.)); + // wtmsk(i,j,k) = std::tanh(arg)*0.5+0.5; + // return; + // } + //} + } + }); + } } } }