From 60371c15cc366f52b1d5753a0e5306f5f48fd826 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Mon, 26 Aug 2024 10:56:43 -0400 Subject: [PATCH 01/32] Improve maxParticleVelocity --- Source/Particles/MultiParticleContainer.H | 2 + Source/Particles/MultiParticleContainer.cpp | 9 +++++ Source/Particles/WarpXParticleContainer.cpp | 44 +++++++++++++++------ 3 files changed, 44 insertions(+), 11 deletions(-) diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index f007d5088e3..97e4e1bc4da 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -89,6 +89,8 @@ public: return allcontainers[index]->meanParticleVelocity(); } + amrex::ParticleReal maxParticleVelocity(); + void AllocData (); void InitData (); diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index fc496217388..5a48232d49a 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -399,6 +399,15 @@ MultiParticleContainer::GetParticleContainerFromName (const std::string& name) c return *allcontainers[i]; } +amrex::ParticleReal +MultiParticleContainer::maxParticleVelocity() { + amrex::ParticleReal max_v = 0.0_prt; + for (const auto &pc : allcontainers) { + max_v = std::max(max_v, pc->maxParticleVelocity()); + } + return max_v; +} + void MultiParticleContainer::AllocData () { diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index bdce18b7b2b..1cab68f1509 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -1422,26 +1422,48 @@ std::array WarpXParticleContainer::meanParticleVelocity(bool lo amrex::ParticleReal WarpXParticleContainer::maxParticleVelocity(bool local) { - amrex::ParticleReal max_v = 0.0; + amrex::ParticleReal max_usq = 0.0; + const amrex::ParticleReal inv_clight_sq = 1.0_prt/(PhysConst::c*PhysConst::c); const int nLevels = finestLevel(); - for (int lev = 0; lev <= nLevels; ++lev) - { -#ifdef AMREX_USE_OMP -#pragma omp parallel reduction(max:max_v) +#ifdef AMREX_USE_GPU + if (Gpu::inLaunchRegion()) + { + using PType = typename WarpXParticleContainer::ParticleType; + max_usq = amrex::ReduceMax(*this, + [=] AMREX_GPU_HOST_DEVICE (const PType& p) -> amrex::ParticleReal + { + auto ux = p.rdata(3); + auto uy = p.rdata(4); + auto uz = p.rdata(5); + return (ux*ux + uy*uy + uz*uz) * inv_clight_sq; + }); + } + else #endif - for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) + { + for (int lev = 0; lev <= nLevels; ++lev) { - auto& ux = pti.GetAttribs(PIdx::ux); - auto& uy = pti.GetAttribs(PIdx::uy); - auto& uz = pti.GetAttribs(PIdx::uz); - for (unsigned long i = 0; i < ux.size(); i++) { - max_v = std::max(max_v, std::sqrt(ux[i]*ux[i] + uy[i]*uy[i] + uz[i]*uz[i])); +#ifdef AMREX_USE_OMP +#pragma omp parallel reduction(max:max_usq) +#endif + for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) + { + auto& ux = pti.GetAttribs(PIdx::ux); + auto& uy = pti.GetAttribs(PIdx::uy); + auto& uz = pti.GetAttribs(PIdx::uz); + for (unsigned long i = 0; i < ux.size(); i++) { + const amrex::ParticleReal usq = (ux[i]*ux[i] + uy[i]*uy[i] + uz[i]*uz[i]) * inv_clight_sq; + max_usq = std::max(max_usq, usq); + } } } } + const amrex::ParticleReal gaminv = 1.0_prt/std::sqrt(1.0_prt + max_usq); + amrex::ParticleReal max_v = gaminv * std::sqrt(max_usq) * PhysConst::c; + if (!local) { ParallelAllReduce::Max(max_v, ParallelDescriptor::Communicator()); } return max_v; } From 3eb3905176d94a41ca6a20245e641e9c65d77517 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Mon, 26 Aug 2024 11:58:05 -0400 Subject: [PATCH 02/32] gpu fix 1 --- Source/Particles/WarpXParticleContainer.cpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 1cab68f1509..d807d41f540 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -1430,13 +1430,14 @@ amrex::ParticleReal WarpXParticleContainer::maxParticleVelocity(bool local) { #ifdef AMREX_USE_GPU if (Gpu::inLaunchRegion()) { - using PType = typename WarpXParticleContainer::ParticleType; + using PTDType = typename WarpXParticleContainer::ParticleTileType::ParticleTileDataType; max_usq = amrex::ReduceMax(*this, - [=] AMREX_GPU_HOST_DEVICE (const PType& p) -> amrex::ParticleReal + [=] AMREX_GPU_HOST_DEVICE (const PTDType& ptd, const int i) -> amrex::ParticleReal { - auto ux = p.rdata(3); - auto uy = p.rdata(4); - auto uz = p.rdata(5); + const auto& aos = ptd.m_aos[i]; + auto ux = aos.rdata(3); + auto uy = aos.rdata(4); + auto uz = aos.rdata(5); return (ux*ux + uy*uy + uz*uz) * inv_clight_sq; }); } From 91841f0180452cefcbcc488c68ffbb8d08a6515e Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Mon, 26 Aug 2024 12:09:05 -0400 Subject: [PATCH 03/32] gpu fix 2 --- Source/Particles/WarpXParticleContainer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index d807d41f540..47ee1df8705 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -1430,7 +1430,7 @@ amrex::ParticleReal WarpXParticleContainer::maxParticleVelocity(bool local) { #ifdef AMREX_USE_GPU if (Gpu::inLaunchRegion()) { - using PTDType = typename WarpXParticleContainer::ParticleTileType::ParticleTileDataType; + using PTDType = typename WarpXParticleContainer::ParticleTileType::ConstParticleTileDataType; max_usq = amrex::ReduceMax(*this, [=] AMREX_GPU_HOST_DEVICE (const PTDType& ptd, const int i) -> amrex::ParticleReal { From 20b4543e4ee5f4452966b02573d29d3bfbde84ee Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Mon, 26 Aug 2024 14:16:52 -0400 Subject: [PATCH 04/32] gpu fix 3 - change to general reduction --- Source/Particles/WarpXParticleContainer.cpp | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 47ee1df8705..d5ccb84db8b 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -1430,16 +1430,18 @@ amrex::ParticleReal WarpXParticleContainer::maxParticleVelocity(bool local) { #ifdef AMREX_USE_GPU if (Gpu::inLaunchRegion()) { - using PTDType = typename WarpXParticleContainer::ParticleTileType::ConstParticleTileDataType; - max_usq = amrex::ReduceMax(*this, - [=] AMREX_GPU_HOST_DEVICE (const PTDType& ptd, const int i) -> amrex::ParticleReal - { - const auto& aos = ptd.m_aos[i]; - auto ux = aos.rdata(3); - auto uy = aos.rdata(4); - auto uz = aos.rdata(5); + amrex::ReduceOps reduce_ops; + using PType = typename WarpXParticleContainer::SuperParticleType; + auto reduce_res = amrex::ParticleReduce>( + *this, + [=] AMREX_GPU_HOST_DEVICE(const PType &p) noexcept -> amrex::GpuTuple { + const auto ux = p.rdata(PIdx::ux); + const auto uy = p.rdata(PIdx::uy); + const auto uz = p.rdata(PIdx::uz); return (ux*ux + uy*uy + uz*uz) * inv_clight_sq; - }); + }, + reduce_ops); + max_usq = amrex::get<0>(reduce_res); } else #endif From 6f4a291abd7c0d530763cd0862d7d812fb8f134c Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Mon, 26 Aug 2024 14:38:43 -0400 Subject: [PATCH 05/32] Remove seperate OMP bit --- Source/Particles/WarpXParticleContainer.cpp | 53 +++++---------------- 1 file changed, 13 insertions(+), 40 deletions(-) diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index d5ccb84db8b..19cf6335b45 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -1422,48 +1422,21 @@ std::array WarpXParticleContainer::meanParticleVelocity(bool lo amrex::ParticleReal WarpXParticleContainer::maxParticleVelocity(bool local) { - amrex::ParticleReal max_usq = 0.0; const amrex::ParticleReal inv_clight_sq = 1.0_prt/(PhysConst::c*PhysConst::c); - const int nLevels = finestLevel(); - -#ifdef AMREX_USE_GPU - if (Gpu::inLaunchRegion()) - { - amrex::ReduceOps reduce_ops; - using PType = typename WarpXParticleContainer::SuperParticleType; - auto reduce_res = amrex::ParticleReduce>( - *this, - [=] AMREX_GPU_HOST_DEVICE(const PType &p) noexcept -> amrex::GpuTuple { - const auto ux = p.rdata(PIdx::ux); - const auto uy = p.rdata(PIdx::uy); - const auto uz = p.rdata(PIdx::uz); - return (ux*ux + uy*uy + uz*uz) * inv_clight_sq; - }, - reduce_ops); - max_usq = amrex::get<0>(reduce_res); - } - else -#endif - { - for (int lev = 0; lev <= nLevels; ++lev) - { -#ifdef AMREX_USE_OMP -#pragma omp parallel reduction(max:max_usq) -#endif - for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) - { - auto& ux = pti.GetAttribs(PIdx::ux); - auto& uy = pti.GetAttribs(PIdx::uy); - auto& uz = pti.GetAttribs(PIdx::uz); - for (unsigned long i = 0; i < ux.size(); i++) { - const amrex::ParticleReal usq = (ux[i]*ux[i] + uy[i]*uy[i] + uz[i]*uz[i]) * inv_clight_sq; - max_usq = std::max(max_usq, usq); - } - } - } - } - + using PType = typename WarpXParticleContainer::SuperParticleType; + amrex::ReduceOps reduce_ops; + const auto reduce_res = amrex::ParticleReduce>( + *this, + [=] AMREX_GPU_HOST_DEVICE(const PType &p) noexcept -> amrex::GpuTuple { + const auto ux = p.rdata(PIdx::ux); + const auto uy = p.rdata(PIdx::uy); + const auto uz = p.rdata(PIdx::uz); + return (ux*ux + uy*uy + uz*uz) * inv_clight_sq; + }, + reduce_ops); + + const amrex::ParticleReal max_usq = amrex::get<0>(reduce_res); const amrex::ParticleReal gaminv = 1.0_prt/std::sqrt(1.0_prt + max_usq); amrex::ParticleReal max_v = gaminv * std::sqrt(max_usq) * PhysConst::c; From 791c4f289c3fef89c0707d515bfd174586ecdaab Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Tue, 27 Aug 2024 11:23:15 -0400 Subject: [PATCH 06/32] Add const qualifier for clang tidy --- Source/Particles/WarpXParticleContainer.H | 2 +- Source/Particles/WarpXParticleContainer.cpp | 7 ++++--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index b44cb33d66b..fc0865af980 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -284,7 +284,7 @@ public: std::array meanParticleVelocity(bool local = false); - amrex::ParticleReal maxParticleVelocity(bool local = false); + amrex::ParticleReal maxParticleVelocity(bool local = false) const; /** * \brief Adds n particles to the simulation diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 19cf6335b45..865452400bf 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -1420,7 +1420,7 @@ std::array WarpXParticleContainer::meanParticleVelocity(bool lo return mean_v; } -amrex::ParticleReal WarpXParticleContainer::maxParticleVelocity(bool local) { +amrex::ParticleReal WarpXParticleContainer::maxParticleVelocity(bool local) const { const amrex::ParticleReal inv_clight_sq = 1.0_prt/(PhysConst::c*PhysConst::c); @@ -1428,11 +1428,12 @@ amrex::ParticleReal WarpXParticleContainer::maxParticleVelocity(bool local) { amrex::ReduceOps reduce_ops; const auto reduce_res = amrex::ParticleReduce>( *this, - [=] AMREX_GPU_HOST_DEVICE(const PType &p) noexcept -> amrex::GpuTuple { + [=] AMREX_GPU_HOST_DEVICE(const PType &p) noexcept -> amrex::GpuTuple + { const auto ux = p.rdata(PIdx::ux); const auto uy = p.rdata(PIdx::uy); const auto uz = p.rdata(PIdx::uz); - return (ux*ux + uy*uy + uz*uz) * inv_clight_sq; + return {(ux*ux + uy*uy + uz*uz) * inv_clight_sq,}; }, reduce_ops); From 6d6df500e48f614be24c495b024415fa552ec0c4 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 4 Sep 2024 12:12:16 -0700 Subject: [PATCH 07/32] Refactor maxVelocity function --- Source/Particles/WarpXParticleContainer.H | 2 +- Source/Particles/WarpXParticleContainer.cpp | 38 +++++++++++++-------- 2 files changed, 24 insertions(+), 16 deletions(-) diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index fc0865af980..b44cb33d66b 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -284,7 +284,7 @@ public: std::array meanParticleVelocity(bool local = false); - amrex::ParticleReal maxParticleVelocity(bool local = false) const; + amrex::ParticleReal maxParticleVelocity(bool local = false); /** * \brief Adds n particles to the simulation diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 865452400bf..6659abb9b1f 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -1420,24 +1420,32 @@ std::array WarpXParticleContainer::meanParticleVelocity(bool lo return mean_v; } -amrex::ParticleReal WarpXParticleContainer::maxParticleVelocity(bool local) const { +amrex::ParticleReal WarpXParticleContainer::maxParticleVelocity(bool local) { const amrex::ParticleReal inv_clight_sq = 1.0_prt/(PhysConst::c*PhysConst::c); + ReduceOps reduce_op; + ReduceData reduce_data(reduce_op); + + const int nLevels = finestLevel(); + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for (int lev = 0; lev <= nLevels; ++lev) { + for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) + { + auto *const ux = pti.GetAttribs(PIdx::ux).data(); + auto *const uy = pti.GetAttribs(PIdx::uy).data(); + auto *const uz = pti.GetAttribs(PIdx::uz).data(); + + reduce_op.eval(pti.numParticles(), reduce_data, + [=] AMREX_GPU_DEVICE (int ip) + { return (ux[ip]*ux[ip] + uy[ip]*uy[ip] + uz[ip]*uz[ip]) * inv_clight_sq; }); + } + } + + const amrex::ParticleReal max_usq = amrex::get<0>(reduce_data.value()); - using PType = typename WarpXParticleContainer::SuperParticleType; - amrex::ReduceOps reduce_ops; - const auto reduce_res = amrex::ParticleReduce>( - *this, - [=] AMREX_GPU_HOST_DEVICE(const PType &p) noexcept -> amrex::GpuTuple - { - const auto ux = p.rdata(PIdx::ux); - const auto uy = p.rdata(PIdx::uy); - const auto uz = p.rdata(PIdx::uz); - return {(ux*ux + uy*uy + uz*uz) * inv_clight_sq,}; - }, - reduce_ops); - - const amrex::ParticleReal max_usq = amrex::get<0>(reduce_res); const amrex::ParticleReal gaminv = 1.0_prt/std::sqrt(1.0_prt + max_usq); amrex::ParticleReal max_v = gaminv * std::sqrt(max_usq) * PhysConst::c; From 644cfd4e2dd2c89084aad76f5b9b0dc6bef1d90d Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Tue, 27 Aug 2024 12:11:13 -0400 Subject: [PATCH 08/32] Add dt_next parameter to all Evolve + PushPX methods --- Source/Particles/LaserParticleContainer.H | 2 +- Source/Particles/LaserParticleContainer.cpp | 3 ++- Source/Particles/MultiParticleContainer.H | 2 +- Source/Particles/MultiParticleContainer.cpp | 4 ++-- Source/Particles/PhotonParticleContainer.H | 4 +++- Source/Particles/PhotonParticleContainer.cpp | 7 ++++--- Source/Particles/PhysicalParticleContainer.H | 5 ++++- Source/Particles/PhysicalParticleContainer.cpp | 14 ++++++++------ Source/Particles/RigidInjectedParticleContainer.H | 5 ++++- .../Particles/RigidInjectedParticleContainer.cpp | 9 +++++---- Source/Particles/WarpXParticleContainer.H | 2 +- 11 files changed, 35 insertions(+), 22 deletions(-) diff --git a/Source/Particles/LaserParticleContainer.H b/Source/Particles/LaserParticleContainer.H index 197cb897602..42b38113ac2 100644 --- a/Source/Particles/LaserParticleContainer.H +++ b/Source/Particles/LaserParticleContainer.H @@ -72,7 +72,7 @@ public: amrex::MultiFab* rho, amrex::MultiFab* crho, const amrex::MultiFab*, const amrex::MultiFab*, const amrex::MultiFab*, const amrex::MultiFab*, const amrex::MultiFab*, const amrex::MultiFab*, - amrex::Real t, amrex::Real dt, DtType a_dt_type=DtType::Full, + amrex::Real t, amrex::Real dt, amrex::Real dt_next, DtType a_dt_type=DtType::Full, bool skip_deposition=false, PushType push_type=PushType::Explicit) final; void PushP (int lev, amrex::Real dt, diff --git a/Source/Particles/LaserParticleContainer.cpp b/Source/Particles/LaserParticleContainer.cpp index e9509a1ef40..73310f72620 100644 --- a/Source/Particles/LaserParticleContainer.cpp +++ b/Source/Particles/LaserParticleContainer.cpp @@ -566,7 +566,8 @@ LaserParticleContainer::Evolve (int lev, MultiFab* rho, MultiFab* crho, const MultiFab*, const MultiFab*, const MultiFab*, const MultiFab*, const MultiFab*, const MultiFab*, - Real t, Real dt, DtType /*a_dt_type*/, bool skip_deposition, PushType push_type) + Real t, Real dt, Real /*dt_next*/, DtType /*a_dt_type*/, + bool skip_deposition, PushType push_type) { WARPX_PROFILE("LaserParticleContainer::Evolve()"); WARPX_PROFILE_VAR_NS("LaserParticleContainer::Evolve::ParticlePush", blp_pp); diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 97e4e1bc4da..ce84dcb7460 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -110,7 +110,7 @@ public: amrex::MultiFab* rho, amrex::MultiFab* crho, const amrex::MultiFab* cEx, const amrex::MultiFab* cEy, const amrex::MultiFab* cEz, const amrex::MultiFab* cBx, const amrex::MultiFab* cBy, const amrex::MultiFab* cBz, - amrex::Real t, amrex::Real dt, DtType a_dt_type=DtType::Full, bool skip_deposition=false, + amrex::Real t, amrex::Real dt, amrex::Real dt_next, DtType a_dt_type=DtType::Full, bool skip_deposition=false, PushType push_type=PushType::Explicit); /** diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 5a48232d49a..bc502496882 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -467,7 +467,7 @@ MultiParticleContainer::Evolve (int lev, MultiFab* rho, MultiFab* crho, const MultiFab* cEx, const MultiFab* cEy, const MultiFab* cEz, const MultiFab* cBx, const MultiFab* cBy, const MultiFab* cBz, - Real t, Real dt, DtType a_dt_type, bool skip_deposition, + Real t, Real dt, Real dt_next, DtType a_dt_type, bool skip_deposition, PushType push_type) { if (! skip_deposition) { @@ -482,7 +482,7 @@ MultiParticleContainer::Evolve (int lev, } for (auto& pc : allcontainers) { pc->Evolve(lev, Ex, Ey, Ez, Bx, By, Bz, jx, jy, jz, cjx, cjy, cjz, - rho, crho, cEx, cEy, cEz, cBx, cBy, cBz, t, dt, a_dt_type, skip_deposition, push_type); + rho, crho, cEx, cEy, cEz, cBx, cBy, cBz, t, dt, dt_next, a_dt_type, skip_deposition, push_type); } } diff --git a/Source/Particles/PhotonParticleContainer.H b/Source/Particles/PhotonParticleContainer.H index 34afac53482..e0be5076909 100644 --- a/Source/Particles/PhotonParticleContainer.H +++ b/Source/Particles/PhotonParticleContainer.H @@ -69,6 +69,7 @@ public: const amrex::MultiFab* cBz, amrex::Real t, amrex::Real dt, + amrex::Real dt_next, DtType a_dt_type=DtType::Full, bool skip_deposition=false, PushType push_type=PushType::Explicit) override; @@ -84,7 +85,8 @@ public: long offset, long np_to_push, int lev, int gather_lev, - amrex::Real dt, ScaleFields scaleFields, + amrex::Real dt, amrex::Real dt_next, + ScaleFields scaleFields, DtType a_dt_type) override; // Do nothing diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index 1f15d5210f5..f71a7d36595 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -91,7 +91,8 @@ PhotonParticleContainer::PushPX (WarpXParIter& pti, const long offset, const long np_to_push, int lev, int gather_lev, - amrex::Real dt, ScaleFields /*scaleFields*/, DtType a_dt_type) + amrex::Real dt, amrex::Real /*dt_next*/, + ScaleFields /*scaleFields*/, DtType a_dt_type) { // Get inverse cell size on gather_lev const amrex::XDim3 dinv = WarpX::InvCellSize(std::max(gather_lev,0)); @@ -237,7 +238,7 @@ PhotonParticleContainer::Evolve (int lev, MultiFab* rho, MultiFab* crho, const MultiFab* cEx, const MultiFab* cEy, const MultiFab* cEz, const MultiFab* cBx, const MultiFab* cBy, const MultiFab* cBz, - Real t, Real dt, DtType a_dt_type, bool skip_deposition, + Real t, Real dt, Real dt_next, DtType a_dt_type, bool skip_deposition, PushType push_type) { // This does gather, push and deposit. @@ -250,6 +251,6 @@ PhotonParticleContainer::Evolve (int lev, rho, crho, cEx, cEy, cEz, cBx, cBy, cBz, - t, dt, a_dt_type, skip_deposition, push_type); + t, dt, dt_next, a_dt_type, skip_deposition, push_type); } diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 5d9b41b8b75..a04e00a8e05 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -135,6 +135,7 @@ public: const amrex::MultiFab* cBz, amrex::Real t, amrex::Real dt, + amrex::Real dt_next, DtType a_dt_type=DtType::Full, bool skip_deposition=false, PushType push_type=PushType::Explicit) override; @@ -150,7 +151,9 @@ public: long offset, long np_to_push, int lev, int gather_lev, - amrex::Real dt, ScaleFields scaleFields, + amrex::Real dt, + amrex::Real dt_next, + ScaleFields scaleFields, DtType a_dt_type=DtType::Full); void ImplicitPushXP (WarpXParIter& pti, diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 3a30467c20d..7d2e8f0fcdb 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -2020,7 +2020,7 @@ PhysicalParticleContainer::Evolve (int lev, MultiFab* rho, MultiFab* crho, const MultiFab* cEx, const MultiFab* cEy, const MultiFab* cEz, const MultiFab* cBx, const MultiFab* cBy, const MultiFab* cBz, - Real /*t*/, Real dt, DtType a_dt_type, bool skip_deposition, + Real /*t*/, Real dt, Real dt_next, DtType a_dt_type, bool skip_deposition, PushType push_type) { @@ -2154,7 +2154,7 @@ PhysicalParticleContainer::Evolve (int lev, PushPX(pti, exfab, eyfab, ezfab, bxfab, byfab, bzfab, Ex.nGrowVect(), e_is_nodal, - 0, np_to_push, lev, gather_lev, dt, ScaleFields(false), a_dt_type); + 0, np_to_push, lev, gather_lev, dt, dt_next, ScaleFields(false), a_dt_type); } else if (push_type == PushType::Implicit) { ImplicitPushXP(pti, exfab, eyfab, ezfab, bxfab, byfab, bzfab, @@ -2196,7 +2196,7 @@ PhysicalParticleContainer::Evolve (int lev, cbxfab, cbyfab, cbzfab, cEx->nGrowVect(), e_is_nodal, nfine_gather, np-nfine_gather, - lev, lev-1, dt, ScaleFields(false), a_dt_type); + lev, lev-1, dt, dt_next, ScaleFields(false), a_dt_type); } else if (push_type == PushType::Implicit) { ImplicitPushXP(pti, cexfab, ceyfab, cezfab, cbxfab, cbyfab, cbzfab, @@ -2742,7 +2742,9 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, const long offset, const long np_to_push, int lev, int gather_lev, - amrex::Real dt, ScaleFields scaleFields, + amrex::Real dt, + amrex::Real dt_next, + ScaleFields scaleFields, DtType a_dt_type) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE((gather_lev==(lev-1)) || @@ -2931,7 +2933,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, #endif dt); - UpdatePosition(xp, yp, zp, ux[ip], uy[ip], uz[ip], dt); + UpdatePosition(xp, yp, zp, ux[ip], uy[ip], uz[ip], 0.5 * (dt + dt_next)); setPosition(ip, xp, yp, zp); } #ifdef WARPX_QED @@ -2949,7 +2951,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, t_chi_max, dt); - UpdatePosition(xp, yp, zp, ux[ip], uy[ip], uz[ip], dt); + UpdatePosition(xp, yp, zp, ux[ip], uy[ip], uz[ip], 0.5 * (dt + dt_next)); setPosition(ip, xp, yp, zp); } } diff --git a/Source/Particles/RigidInjectedParticleContainer.H b/Source/Particles/RigidInjectedParticleContainer.H index bc20420ea6e..e18931f1f01 100644 --- a/Source/Particles/RigidInjectedParticleContainer.H +++ b/Source/Particles/RigidInjectedParticleContainer.H @@ -84,6 +84,7 @@ public: const amrex::MultiFab* cBz, amrex::Real t, amrex::Real dt, + amrex::Real dt_next, DtType a_dt_type=DtType::Full, bool skip_deposition=false, PushType push_type=PushType::Explicit) override; @@ -99,7 +100,9 @@ public: long offset, long np_to_push, int lev, int gather_lev, - amrex::Real dt, ScaleFields scaleFields, + amrex::Real dt, + amrex::Real dt_next, + ScaleFields scaleFields, DtType a_dt_type=DtType::Full) override; void PushP (int lev, amrex::Real dt, diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index c3ec4c41131..c2a50173759 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -166,7 +166,8 @@ RigidInjectedParticleContainer::PushPX (WarpXParIter& pti, const long offset, const long np_to_push, int lev, int gather_lev, - amrex::Real dt, ScaleFields /*scaleFields*/, + amrex::Real dt, amrex::Real /*dt_next*/, + ScaleFields /*scaleFields*/, DtType a_dt_type) { auto& attribs = pti.GetAttribs(); @@ -242,7 +243,7 @@ RigidInjectedParticleContainer::PushPX (WarpXParIter& pti, const bool do_scale = not done_injecting_lev; const Real v_boost = WarpX::beta_boost*PhysConst::c; PhysicalParticleContainer::PushPX(pti, exfab, eyfab, ezfab, bxfab, byfab, bzfab, - ngEB, e_is_nodal, offset, np_to_push, lev, gather_lev, dt, + ngEB, e_is_nodal, offset, np_to_push, lev, gather_lev, dt, dt, ScaleFields(do_scale, dt, zinject_plane_lev_previous, vzbeam_ave_boosted, v_boost), a_dt_type); @@ -299,7 +300,7 @@ RigidInjectedParticleContainer::Evolve (int lev, MultiFab* rho, MultiFab* crho, const MultiFab* cEx, const MultiFab* cEy, const MultiFab* cEz, const MultiFab* cBx, const MultiFab* cBy, const MultiFab* cBz, - Real t, Real dt, DtType a_dt_type, bool skip_deposition, + Real t, Real dt, Real dt_next, DtType a_dt_type, bool skip_deposition, PushType push_type) { @@ -325,7 +326,7 @@ RigidInjectedParticleContainer::Evolve (int lev, rho, crho, cEx, cEy, cEz, cBx, cBy, cBz, - t, dt, a_dt_type, skip_deposition, push_type); + t, dt, dt_next, a_dt_type, skip_deposition, push_type); } void diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 7e882c151e8..98fb3339e9c 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -153,7 +153,7 @@ public: amrex::MultiFab* rho, amrex::MultiFab* crho, const amrex::MultiFab* cEx, const amrex::MultiFab* cEy, const amrex::MultiFab* cEz, const amrex::MultiFab* cBx, const amrex::MultiFab* cBy, const amrex::MultiFab* cBz, - amrex::Real t, amrex::Real dt, DtType a_dt_type=DtType::Full, bool skip_deposition=false, + amrex::Real t, amrex::Real dt, amrex::Real dt_next, DtType a_dt_type=DtType::Full, bool skip_deposition=false, PushType push_type=PushType::Explicit) = 0; virtual void PostRestart () = 0; From e6804213faf354f5d79b89e8b683f56d7fbb5810 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Tue, 27 Aug 2024 13:02:40 -0400 Subject: [PATCH 09/32] Write first pass update dt function --- Source/Evolve/WarpXComputeDt.cpp | 57 ++++++++++++++++++++++++++------ Source/WarpX.H | 4 +++ 2 files changed, 51 insertions(+), 10 deletions(-) diff --git a/Source/Evolve/WarpXComputeDt.cpp b/Source/Evolve/WarpXComputeDt.cpp index c1a87166920..998c3d1d279 100644 --- a/Source/Evolve/WarpXComputeDt.cpp +++ b/Source/Evolve/WarpXComputeDt.cpp @@ -13,6 +13,7 @@ #else # include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H" #endif +#include "Particles/MultiParticleContainer.H" #include "Utils/TextMsg.H" #include "Utils/WarpXAlgorithmSelection.H" #include "Utils/WarpXConst.H" @@ -27,6 +28,18 @@ #include #include +AMREX_FORCE_INLINE amrex::Real +minDim (const amrex::Real* x) +{ +#if defined(WARPX_DIM_1D_Z) + return x[0]; +#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + return std::min(x[0], x[1]); +#else + return std::min(x[0], std::min(x[1], x[2])); +#endif +} + /** * Determine the timestep of the simulation. */ void @@ -46,9 +59,11 @@ WarpX::ComputeDt () } WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_const_dt.has_value(), errorMsg.str()); - for (int lev=0; lev<=max_level; lev++) { - dt[lev] = m_const_dt.value(); - } + dt.resize(0); + dt_next.resize(0); + dt.resize(max_level+1, m_const_dt.value()); + dt_next.resize(max_level+1, m_const_dt.value()); + return; } @@ -61,13 +76,7 @@ WarpX::ComputeDt () } else if (electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) { // Computation of dt for spectral algorithm // (determined by the minimum cell size in all directions) -#if defined(WARPX_DIM_1D_Z) - deltat = cfl * dx[0] / PhysConst::c; -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - deltat = cfl * std::min(dx[0], dx[1]) / PhysConst::c; -#else - deltat = cfl * std::min(dx[0], std::min(dx[1], dx[2])) / PhysConst::c; -#endif + deltat = cfl / PhysConst::c * minDim(dx); } else { // Computation of dt for FDTD algorithm #ifdef WARPX_DIM_RZ @@ -92,9 +101,37 @@ WarpX::ComputeDt () dt.resize(0); dt.resize(max_level+1,deltat); + dt_next.resize(0); + dt_next.resize(max_level+1,deltat); + if (do_subcycling) { for (int lev = max_level-1; lev >= 0; --lev) { dt[lev] = dt[lev+1] * refRatio(lev)[0]; + dt_next[lev] = dt_next[lev+1] * refRatio(lev)[0]; + } + } +} + +void +WarpX::UpdateDtFromParticleSpeeds () +{ + const amrex::Real* dx = geom[max_level].CellSize(); + amrex::Real dx_min = minDim(dx); + + const amrex::ParticleReal max_v = mypc->maxParticleVelocity(); + const amrex::Real deltat_new = cfl * dx_min / max_v; + + // Set present dt to previous next dt + dt[max_level] = dt_next[max_level]; + dt_next[max_level] = deltat_new; + + for (int lev = max_level-1; lev >= 0; --lev) { + if (do_subcycling) { + dt[lev] = dt[lev+1] * refRatio(lev)[0]; + dt_next[lev] = dt_next[lev+1] * refRatio(lev)[0]; + } else { + dt[lev] = dt[lev+1]; + dt_next[lev] = dt_next[lev+1]; } } } diff --git a/Source/WarpX.H b/Source/WarpX.H index 4573808461e..7aaeccc0e81 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -593,6 +593,9 @@ public: /** Determine the timestep of the simulation. */ void ComputeDt (); + /** Determine the timestep of the simulation. */ + void UpdateDtFromParticleSpeeds (); + /** Print main PIC parameters to stdout */ void PrintMainPICparameters (); @@ -1497,6 +1500,7 @@ private: amrex::Vector t_new; amrex::Vector t_old; amrex::Vector dt; + amrex::Vector dt_next; // Particle container std::unique_ptr mypc; From 6cf016d16da88916368645fe9e48f974471d6823 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Tue, 27 Aug 2024 13:57:53 -0400 Subject: [PATCH 10/32] First pass implementation works, need to add logic to only use it when ES simulation enabled --- Source/Evolve/WarpXEvolve.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 5a2dbdf2f30..164559b068b 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -302,6 +302,9 @@ WarpX::Evolve (int numsteps) if (checkStopSimulation(cur_time)) { break; } + + UpdateDtFromParticleSpeeds(); + } // End loop on time steps // This if statement is needed for PICMI, which allows the Evolve routine to be @@ -1030,7 +1033,7 @@ WarpX::PushParticlesandDeposit (int lev, amrex::Real cur_time, DtType a_dt_type, rho_fp[lev].get(), charge_buf[lev].get(), Efield_cax[lev][0].get(), Efield_cax[lev][1].get(), Efield_cax[lev][2].get(), Bfield_cax[lev][0].get(), Bfield_cax[lev][1].get(), Bfield_cax[lev][2].get(), - cur_time, dt[lev], a_dt_type, skip_current, push_type); + cur_time, dt[lev], dt_next[lev], a_dt_type, skip_current, push_type); if (! skip_current) { #ifdef WARPX_DIM_RZ // This is called after all particles have deposited their current and charge. From 9839b0fa8df341eb4fcfe48397df39199112bf93 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Tue, 27 Aug 2024 14:08:46 -0400 Subject: [PATCH 11/32] Only call update dt if no const dt supplied and solver is electrostatic --- Source/Evolve/WarpXEvolve.cpp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 164559b068b..1123ca544c4 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -299,12 +299,15 @@ WarpX::Evolve (int numsteps) << " s; Avg. per step = " << evolve_time/(step-step_begin+1) << " s\n\n"; } + // Update timestep for electrostatic solver if a constant dt is not provided + if (electrostatic_solver_id != ElectrostaticSolverAlgo::None && + !m_const_dt.has_value()) { + UpdateDtFromParticleSpeeds(); + } + if (checkStopSimulation(cur_time)) { break; } - - UpdateDtFromParticleSpeeds(); - } // End loop on time steps // This if statement is needed for PICMI, which allows the Evolve routine to be From 2d0b59925a57e535d47311231d9bc24bd8e7ca77 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Tue, 27 Aug 2024 14:29:36 -0400 Subject: [PATCH 12/32] No longer errors if const_dt not specified for ES solver. still need to add CFL input to picmi --- Source/Evolve/WarpXComputeDt.cpp | 24 +++++++++--------------- Source/WarpX.H | 5 ++++- 2 files changed, 13 insertions(+), 16 deletions(-) diff --git a/Source/Evolve/WarpXComputeDt.cpp b/Source/Evolve/WarpXComputeDt.cpp index 998c3d1d279..6fd4303a8e8 100644 --- a/Source/Evolve/WarpXComputeDt.cpp +++ b/Source/Evolve/WarpXComputeDt.cpp @@ -46,25 +46,15 @@ void WarpX::ComputeDt () { // Handle cases where the timestep is not limited by the speed of light - if (electromagnetic_solver_id == ElectromagneticSolverAlgo::None || - electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC) { - + // and no constant timestep is provided + if (electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC) { std::stringstream errorMsg; - if (electrostatic_solver_id != ElectrostaticSolverAlgo::None) { - errorMsg << "warpx.const_dt must be specified with the electrostatic solver."; - } else if (electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC) { + if (electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC) { errorMsg << "warpx.const_dt must be specified with the hybrid-PIC solver."; } else { errorMsg << "warpx.const_dt must be specified when not using a field solver."; } WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_const_dt.has_value(), errorMsg.str()); - - dt.resize(0); - dt_next.resize(0); - dt.resize(max_level+1, m_const_dt.value()); - dt_next.resize(max_level+1, m_const_dt.value()); - - return; } // Determine the appropriate timestep as limited by the speed of light @@ -73,7 +63,8 @@ WarpX::ComputeDt () if (m_const_dt.has_value()) { deltat = m_const_dt.value(); - } else if (electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) { + } else if (electrostatic_solver_id != ElectrostaticSolverAlgo::None || + electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) { // Computation of dt for spectral algorithm // (determined by the minimum cell size in all directions) deltat = cfl / PhysConst::c * minDim(dx); @@ -100,7 +91,6 @@ WarpX::ComputeDt () dt.resize(0); dt.resize(max_level+1,deltat); - dt_next.resize(0); dt_next.resize(max_level+1,deltat); @@ -112,6 +102,10 @@ WarpX::ComputeDt () } } +/** + * Determine the simulation timestep from the maximum speed of all particles + * Sets timestep so that a particle can only cross cfl*dx cells per timestep. + */ void WarpX::UpdateDtFromParticleSpeeds () { diff --git a/Source/WarpX.H b/Source/WarpX.H index 7aaeccc0e81..4f5a5075524 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -593,7 +593,10 @@ public: /** Determine the timestep of the simulation. */ void ComputeDt (); - /** Determine the timestep of the simulation. */ + /** + * Determine the simulation timestep from the maximum speed of all particles + * Sets timestep so that a particle can only cross cfl*dx cells per timestep. + */ void UpdateDtFromParticleSpeeds (); /** Print main PIC parameters to stdout */ From d65041622a82b8c51298967eebb38f790031a35b Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Tue, 27 Aug 2024 14:44:24 -0400 Subject: [PATCH 13/32] warpx_cfl option added to picmi --- Python/pywarpx/picmi.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index 0d51a8723b4..bd34029c45a 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -1887,6 +1887,7 @@ def init(self, kw): self.absolute_tolerance = kw.pop("warpx_absolute_tolerance", None) self.self_fields_verbosity = kw.pop("warpx_self_fields_verbosity", None) self.magnetostatic = kw.pop("warpx_magnetostatic", False) + self.cfl = kw.pop("warpx_cfl", None) def solver_initialize_inputs(self): # Open BC means FieldBoundaryType::Open for electrostatic sims, rather than perfectly-matched layer @@ -1894,6 +1895,9 @@ def solver_initialize_inputs(self): self.grid.grid_initialize_inputs() + if self.cfl is not None: + pywarpx.warpx.cfl = self.cfl + if self.relativistic: pywarpx.warpx.do_electrostatic = "relativistic" else: From 39c751d84a1ee98dc346a50693aeb09d42c5a4fc Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Fri, 6 Sep 2024 11:02:07 -0700 Subject: [PATCH 14/32] Address most comments from @dpgrote --- Python/pywarpx/picmi.py | 3 --- Source/Evolve/WarpXComputeDt.cpp | 20 ++++++------------- .../Particles/PhysicalParticleContainer.cpp | 4 ++-- 3 files changed, 8 insertions(+), 19 deletions(-) diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index bd34029c45a..91ef83620db 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -1522,9 +1522,6 @@ def solver_initialize_inputs(self): # --- Same method names are used, though mapped to lower case. pywarpx.algo.maxwell_solver = self.method - if self.cfl is not None: - pywarpx.warpx.cfl = self.cfl - if self.source_smoother is not None: self.source_smoother.smoother_initialize_inputs(self) diff --git a/Source/Evolve/WarpXComputeDt.cpp b/Source/Evolve/WarpXComputeDt.cpp index 6fd4303a8e8..cc603b73887 100644 --- a/Source/Evolve/WarpXComputeDt.cpp +++ b/Source/Evolve/WarpXComputeDt.cpp @@ -28,16 +28,13 @@ #include #include +/** + * Compute the minimum of array x, where x has dimension AMREX_SPACEDIM + */ AMREX_FORCE_INLINE amrex::Real minDim (const amrex::Real* x) { -#if defined(WARPX_DIM_1D_Z) - return x[0]; -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - return std::min(x[0], x[1]); -#else - return std::min(x[0], std::min(x[1], x[2])); -#endif + return std::min({AMREX_D_DECL(x[0], x[1], x[2])}); } /** @@ -49,11 +46,7 @@ WarpX::ComputeDt () // and no constant timestep is provided if (electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC) { std::stringstream errorMsg; - if (electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC) { - errorMsg << "warpx.const_dt must be specified with the hybrid-PIC solver."; - } else { - errorMsg << "warpx.const_dt must be specified when not using a field solver."; - } + errorMsg << "warpx.const_dt must be specified with the hybrid-PIC solver."; WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_const_dt.has_value(), errorMsg.str()); } @@ -120,11 +113,10 @@ WarpX::UpdateDtFromParticleSpeeds () dt_next[max_level] = deltat_new; for (int lev = max_level-1; lev >= 0; --lev) { + dt[lev] = dt_next[lev]; if (do_subcycling) { - dt[lev] = dt[lev+1] * refRatio(lev)[0]; dt_next[lev] = dt_next[lev+1] * refRatio(lev)[0]; } else { - dt[lev] = dt[lev+1]; dt_next[lev] = dt_next[lev+1]; } } diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 7d2e8f0fcdb..8e5393b61a9 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -2933,7 +2933,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, #endif dt); - UpdatePosition(xp, yp, zp, ux[ip], uy[ip], uz[ip], 0.5 * (dt + dt_next)); + UpdatePosition(xp, yp, zp, ux[ip], uy[ip], uz[ip], 0.5_rt * (dt + dt_next)); setPosition(ip, xp, yp, zp); } #ifdef WARPX_QED @@ -2951,7 +2951,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, t_chi_max, dt); - UpdatePosition(xp, yp, zp, ux[ip], uy[ip], uz[ip], 0.5 * (dt + dt_next)); + UpdatePosition(xp, yp, zp, ux[ip], uy[ip], uz[ip], 0.5_rt * (dt + dt_next)); setPosition(ip, xp, yp, zp); } } From 5a4f0041cd5eabc76c802005c969a62ad1fe8114 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Mon, 9 Sep 2024 12:24:05 -0700 Subject: [PATCH 15/32] Move timestep update to start of Evolve --- Source/Evolve/WarpXEvolve.cpp | 12 ++++++------ install-deps.sh | 15 +++++++++++++++ install.sh | 27 +++++++++++++++++++++++++++ 3 files changed, 48 insertions(+), 6 deletions(-) create mode 100755 install-deps.sh create mode 100755 install.sh diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 1123ca544c4..b366a6c3a22 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -86,6 +86,12 @@ WarpX::Evolve (int numsteps) multi_diags->NewIteration(); + // Update timestep for electrostatic solver if a constant dt is not provided + if (electrostatic_solver_id != ElectrostaticSolverAlgo::None && + !m_const_dt.has_value()) { + UpdateDtFromParticleSpeeds(); + } + // Start loop on time steps if (verbose) { amrex::Print() << "STEP " << step+1 << " starts ...\n"; @@ -299,12 +305,6 @@ WarpX::Evolve (int numsteps) << " s; Avg. per step = " << evolve_time/(step-step_begin+1) << " s\n\n"; } - // Update timestep for electrostatic solver if a constant dt is not provided - if (electrostatic_solver_id != ElectrostaticSolverAlgo::None && - !m_const_dt.has_value()) { - UpdateDtFromParticleSpeeds(); - } - if (checkStopSimulation(cur_time)) { break; } diff --git a/install-deps.sh b/install-deps.sh new file mode 100755 index 00000000000..b93d4e98774 --- /dev/null +++ b/install-deps.sh @@ -0,0 +1,15 @@ +brew update +brew tap openpmd/openpmd +brew install adios2 # for openPMD +brew install ccache +brew install cmake +brew install fftw # for PSATD +brew install git +brew install hdf5-mpi # for openPMD +brew install libomp +brew unlink gcc +brew link --force libomp +brew install pkg-config # for fftw +brew install open-mpi +brew install openblas # for PSATD in RZ +brew install openpmd-api # for openPMD diff --git a/install.sh b/install.sh new file mode 100755 index 00000000000..d872dd7b808 --- /dev/null +++ b/install.sh @@ -0,0 +1,27 @@ +#!/bin/bash + +# uninstall old versions +#rm -rf build +#rm -rf *.whl + +# Need to use LLVM's clang instead of apple's clang to use openmp +export CC=/opt/homebrew/Cellar/llvm/18.1.8/bin/clang +export CXX=/opt/homebrew/Cellar/llvm/18.1.8/bin/clang++ + +# Activate python virtual environemnt +source /Users/archermarks/venvs/warpx/bin/activate + +# Build warpx +cmake -S . -B build -Wno-dev \ + -DCMAKE_BUILD_TYPE=RelWithDebInfo \ + -DWarpX_LIB=ON \ + -DWarpX_APP=ON \ + -DWarpX_MPI=ON \ + -DWarpX_DIMS="1;2;3" \ + -DWarpX_PRECISION=DOUBLE \ + -DWarpX_PARTICLE_PRECISION=SINGLE \ + -DWarpX_COMPUTE=OMP + + #-DWarpX_PYTHON=ON \ +cmake --build build -j 20 +#cmake --build build --target pip_install -j 20 From abb53dea6c2191704fcecaabff0281c3f6210e60 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Mon, 9 Sep 2024 12:26:37 -0700 Subject: [PATCH 16/32] Untrack accidentally-tracked files --- install-deps.sh | 15 --------------- install.sh | 27 --------------------------- 2 files changed, 42 deletions(-) delete mode 100755 install-deps.sh delete mode 100755 install.sh diff --git a/install-deps.sh b/install-deps.sh deleted file mode 100755 index b93d4e98774..00000000000 --- a/install-deps.sh +++ /dev/null @@ -1,15 +0,0 @@ -brew update -brew tap openpmd/openpmd -brew install adios2 # for openPMD -brew install ccache -brew install cmake -brew install fftw # for PSATD -brew install git -brew install hdf5-mpi # for openPMD -brew install libomp -brew unlink gcc -brew link --force libomp -brew install pkg-config # for fftw -brew install open-mpi -brew install openblas # for PSATD in RZ -brew install openpmd-api # for openPMD diff --git a/install.sh b/install.sh deleted file mode 100755 index d872dd7b808..00000000000 --- a/install.sh +++ /dev/null @@ -1,27 +0,0 @@ -#!/bin/bash - -# uninstall old versions -#rm -rf build -#rm -rf *.whl - -# Need to use LLVM's clang instead of apple's clang to use openmp -export CC=/opt/homebrew/Cellar/llvm/18.1.8/bin/clang -export CXX=/opt/homebrew/Cellar/llvm/18.1.8/bin/clang++ - -# Activate python virtual environemnt -source /Users/archermarks/venvs/warpx/bin/activate - -# Build warpx -cmake -S . -B build -Wno-dev \ - -DCMAKE_BUILD_TYPE=RelWithDebInfo \ - -DWarpX_LIB=ON \ - -DWarpX_APP=ON \ - -DWarpX_MPI=ON \ - -DWarpX_DIMS="1;2;3" \ - -DWarpX_PRECISION=DOUBLE \ - -DWarpX_PARTICLE_PRECISION=SINGLE \ - -DWarpX_COMPUTE=OMP - - #-DWarpX_PYTHON=ON \ -cmake --build build -j 20 -#cmake --build build --target pip_install -j 20 From cce7cb91e7119fa3938c8be5a65b3fb4ce9e7cd5 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Tue, 10 Sep 2024 10:19:42 -0700 Subject: [PATCH 17/32] Add max_dt param and handle zero-velocity cases --- Python/pywarpx/picmi.py | 6 ++++++ Source/Evolve/WarpXComputeDt.cpp | 17 +++++++++++++++-- Source/WarpX.H | 2 ++ Source/WarpX.cpp | 2 ++ 4 files changed, 25 insertions(+), 2 deletions(-) diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index 91ef83620db..dc4171b02ea 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -2740,6 +2740,9 @@ class Simulation(picmistandard.PICMI_Simulation): warpx_used_inputs_file: string, optional The name of the text file that the used input parameters is written to, + + warpx_max_dt: float, optional + The maximum allowable timestep when using adaptive timestepping (only available for non-electromagnetic solvers) """ # Set the C++ WarpX interface (see _libwarpx.LibWarpX) as an extension to @@ -2819,6 +2822,9 @@ def initialize_inputs(self): if self.time_step_size is not None: pywarpx.warpx.const_dt = self.time_step_size + if self.warpx_max_dt is not None: + pywarpx.warpx.max_dt = self.warpx_max_dt + if self.gamma_boost is not None: pywarpx.warpx.gamma_boost = self.gamma_boost pywarpx.warpx.boost_direction = "z" diff --git a/Source/Evolve/WarpXComputeDt.cpp b/Source/Evolve/WarpXComputeDt.cpp index cc603b73887..3d8d03f8e81 100644 --- a/Source/Evolve/WarpXComputeDt.cpp +++ b/Source/Evolve/WarpXComputeDt.cpp @@ -90,7 +90,7 @@ WarpX::ComputeDt () if (do_subcycling) { for (int lev = max_level-1; lev >= 0; --lev) { dt[lev] = dt[lev+1] * refRatio(lev)[0]; - dt_next[lev] = dt_next[lev+1] * refRatio(lev)[0]; + dt_next[lev] = dt[lev]; } } } @@ -106,7 +106,20 @@ WarpX::UpdateDtFromParticleSpeeds () amrex::Real dx_min = minDim(dx); const amrex::ParticleReal max_v = mypc->maxParticleVelocity(); - const amrex::Real deltat_new = cfl * dx_min / max_v; + amrex::Real deltat_new; + + // Protections from overly-large timesteps + if (max_v == 0) { + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_max_dt.has_value(), "Particles at rest and no maximum timestep specified. Aborting."); + deltat_new == m_max_dt.value(); + } + + deltat_new = cfl * dx_min / max_v; + + // Restrict to be less than user-specified maximum timestep, if present + if (m_max_dt.has_value()) { + deltat_new = std::min(deltat_new, m_max_dt.value()); + } // Set present dt to previous next dt dt[max_level] = dt_next[max_level]; diff --git a/Source/WarpX.H b/Source/WarpX.H index 92045aade72..06b3537bfe8 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -1666,7 +1666,9 @@ private: int num_injected_species = -1; amrex::Vector injected_plasma_species; + // Timestepping parameters std::optional m_const_dt; + std::optional m_max_dt; // Macroscopic properties std::unique_ptr m_macroscopic_properties; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index a705735a541..3c63ec02c82 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -831,7 +831,9 @@ WarpX::ReadParameters () pp_boundary.query("verboncoeur_axis_correction", verboncoeur_axis_correction); #endif + // Read timestepping options utils::parser::queryWithParser(pp_warpx, "const_dt", m_const_dt); + utils::parser::queryWithParser(pp_warpx, "max_dt", m_max_dt); // Filter currently not working with FDTD solver in RZ geometry: turn OFF by default // (see https://github.com/ECP-WarpX/WarpX/issues/1943) From c442ef0b1654176956fd83c4eef7c6eb0f339e33 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Tue, 10 Sep 2024 12:49:41 -0700 Subject: [PATCH 18/32] Fix time-centering of leapfrog push --- .../Diagnostics/ReducedDiags/CMakeLists.txt | 17 ++--- Source/Diagnostics/ReducedDiags/Make.package | 24 +++---- .../ReducedDiags/MultiReducedDiags.cpp | 24 ++++--- Source/Diagnostics/ReducedDiags/Timestep.H | 35 +++++++++ Source/Diagnostics/ReducedDiags/Timestep.cpp | 72 +++++++++++++++++++ Source/Evolve/WarpXComputeDt.cpp | 10 ++- .../Particles/PhysicalParticleContainer.cpp | 6 +- 7 files changed, 152 insertions(+), 36 deletions(-) create mode 100644 Source/Diagnostics/ReducedDiags/Timestep.H create mode 100644 Source/Diagnostics/ReducedDiags/Timestep.cpp diff --git a/Source/Diagnostics/ReducedDiags/CMakeLists.txt b/Source/Diagnostics/ReducedDiags/CMakeLists.txt index 4f0b05f6180..bbf1b6b65b0 100644 --- a/Source/Diagnostics/ReducedDiags/CMakeLists.txt +++ b/Source/Diagnostics/ReducedDiags/CMakeLists.txt @@ -3,26 +3,27 @@ foreach(D IN LISTS WarpX_DIMS) target_sources(lib_${SD} PRIVATE BeamRelevant.cpp + ChargeOnEB.cpp ColliderRelevant.cpp DifferentialLuminosity.cpp FieldEnergy.cpp + FieldMaximum.cpp + FieldMomentum.cpp FieldProbe.cpp FieldProbeParticleContainer.cpp - FieldMomentum.cpp + FieldReduction.cpp + FieldProbe.cpp LoadBalanceCosts.cpp LoadBalanceEfficiency.cpp MultiReducedDiags.cpp ParticleEnergy.cpp - ParticleMomentum.cpp + ParticleExtrema.cpp ParticleHistogram.cpp ParticleHistogram2D.cpp + ParticleMomentum.cpp + ParticleNumber.cpp ReducedDiags.cpp - FieldMaximum.cpp - ParticleExtrema.cpp RhoMaximum.cpp - ParticleNumber.cpp - FieldReduction.cpp - FieldProbe.cpp - ChargeOnEB.cpp + Timestep.cpp ) endforeach() diff --git a/Source/Diagnostics/ReducedDiags/Make.package b/Source/Diagnostics/ReducedDiags/Make.package index e840931f8d3..2611831a3dd 100644 --- a/Source/Diagnostics/ReducedDiags/Make.package +++ b/Source/Diagnostics/ReducedDiags/Make.package @@ -1,24 +1,24 @@ CEXE_sources += MultiReducedDiags.cpp CEXE_sources += ReducedDiags.cpp -CEXE_sources += ParticleEnergy.cpp -CEXE_sources += ParticleMomentum.cpp -CEXE_sources += FieldEnergy.cpp -CEXE_sources += FieldProbe.cpp -CEXE_sources += FieldProbeParticleContainer.cpp -CEXE_sources += FieldMomentum.cpp CEXE_sources += BeamRelevant.cpp +CEXE_sources += ChargeOnEB.cpp CEXE_sources += ColliderRelevant.cpp CEXE_sources += DifferentialLuminosity.cpp +CEXE_sources += FieldEnergy.cpp +CEXE_sources += FieldMaximum.cpp +CEXE_sources += FieldMomentum.cpp +CEXE_sources += FieldProbe.cpp +CEXE_sources += FieldProbeParticleContainer.cpp +CEXE_sources += FieldReduction.cpp CEXE_sources += LoadBalanceCosts.cpp CEXE_sources += LoadBalanceEfficiency.cpp +CEXE_sources += ParticleEnergy.cpp +CEXE_sources += ParticleExtrema.cpp CEXE_sources += ParticleHistogram.cpp CEXE_sources += ParticleHistogram2D.cpp -CEXE_sources += FieldMaximum.cpp -CEXE_sources += FieldProbe.cpp -CEXE_sources += ParticleExtrema.cpp -CEXE_sources += RhoMaximum.cpp +CEXE_sources += ParticleMomentum.cpp CEXE_sources += ParticleNumber.cpp -CEXE_sources += FieldReduction.cpp -CEXE_sources += ChargeOnEB.cpp +CEXE_sources += RhoMaximum.cpp +CEXE_sources += Timestep.cpp VPATH_LOCATIONS += $(WARPX_HOME)/Source/Diagnostics/ReducedDiags diff --git a/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp b/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp index 25ea87d9f54..5035eac58a8 100644 --- a/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp +++ b/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp @@ -12,8 +12,8 @@ #include "DifferentialLuminosity.H" #include "FieldEnergy.H" #include "FieldMaximum.H" -#include "FieldProbe.H" #include "FieldMomentum.H" +#include "FieldProbe.H" #include "FieldReduction.H" #include "LoadBalanceCosts.H" #include "LoadBalanceEfficiency.H" @@ -24,6 +24,7 @@ #include "ParticleMomentum.H" #include "ParticleNumber.H" #include "RhoMaximum.H" +#include "Timestep.H" #include "Utils/TextMsg.H" #include "Utils/WarpXProfilerWrapper.H" @@ -52,24 +53,25 @@ MultiReducedDiags::MultiReducedDiags () using CS = const std::string& ; const auto reduced_diags_dictionary = std::map(CS)>>{ + {"BeamRelevant", [](CS s){return std::make_unique(s);}}, + {"ChargeOnEB", [](CS s){return std::make_unique(s);}}, + {"ColliderRelevant", [](CS s){return std::make_unique(s);}}, + {"DifferentialLuminosity",[](CS s){return std::make_unique(s);}}, {"ParticleEnergy", [](CS s){return std::make_unique(s);}}, + {"ParticleExtrema", [](CS s){return std::make_unique(s);}}, + {"ParticleHistogram", [](CS s){return std::make_unique(s);}}, + {"ParticleHistogram2D", [](CS s){return std::make_unique(s);}}, {"ParticleMomentum", [](CS s){return std::make_unique(s);}}, + {"ParticleNumber", [](CS s){return std::make_unique(s);}}, {"FieldEnergy", [](CS s){return std::make_unique(s);}}, - {"FieldMomentum", [](CS s){return std::make_unique(s);}}, {"FieldMaximum", [](CS s){return std::make_unique(s);}}, + {"FieldMomentum", [](CS s){return std::make_unique(s);}}, {"FieldProbe", [](CS s){return std::make_unique(s);}}, {"FieldReduction", [](CS s){return std::make_unique(s);}}, - {"RhoMaximum", [](CS s){return std::make_unique(s);}}, - {"BeamRelevant", [](CS s){return std::make_unique(s);}}, - {"ColliderRelevant", [](CS s){return std::make_unique(s);}}, - {"DifferentialLuminosity",[](CS s){return std::make_unique(s);}}, {"LoadBalanceCosts", [](CS s){return std::make_unique(s);}}, {"LoadBalanceEfficiency", [](CS s){return std::make_unique(s);}}, - {"ParticleHistogram", [](CS s){return std::make_unique(s);}}, - {"ParticleHistogram2D", [](CS s){return std::make_unique(s);}}, - {"ParticleNumber", [](CS s){return std::make_unique(s);}}, - {"ParticleExtrema", [](CS s){return std::make_unique(s);}}, - {"ChargeOnEB", [](CS s){return std::make_unique(s);}} + {"RhoMaximum", [](CS s){return std::make_unique(s);}}, + {"Timestep", [](CS s){return std::make_unique(s);}} }; // loop over all reduced diags and fill m_multi_rd with requested reduced diags std::transform(m_rd_names.begin(), m_rd_names.end(), std::back_inserter(m_multi_rd), diff --git a/Source/Diagnostics/ReducedDiags/Timestep.H b/Source/Diagnostics/ReducedDiags/Timestep.H new file mode 100644 index 00000000000..bcf4fe6452f --- /dev/null +++ b/Source/Diagnostics/ReducedDiags/Timestep.H @@ -0,0 +1,35 @@ +/* Copyright 2024 The WarpX Community + * + * This file is part of WarpX. + * + * Authors: Thomas Marks + * + * License: BSD-3-Clause-LBNL + */ + +#ifndef WARPX_DIAGNOSTICS_REDUCEDDIAGS_TIMESTEP_H_ +#define WARPX_DIAGNOSTICS_REDUCEDDIAGS_TIMESTEP_H_ + +#include "ReducedDiags.H" +#include + +/** + * This class contains a function for retrieving the current simulation timestep as a diagnostic. + * Useful mainly for simulations using adaptive timestepping. + */ +class Timestep: public ReducedDiags { +public: + /** + * constructor + * @param[in] rd_name reduced diags name + */ + Timestep (const std::string& rd_name); + + /** + * This function gets the current physical timestep of the simulation at all refinement levels. + * @param[in] step current time step + */ + void ComputeDiags (int step) final; +}; + +#endif //WARPX_DIAGNOSTICS_REDUCEDDIAGS_TIMESTEP_H_ diff --git a/Source/Diagnostics/ReducedDiags/Timestep.cpp b/Source/Diagnostics/ReducedDiags/Timestep.cpp new file mode 100644 index 00000000000..3474121db91 --- /dev/null +++ b/Source/Diagnostics/ReducedDiags/Timestep.cpp @@ -0,0 +1,72 @@ +/* Copyright 2024 The WarpX Community + * + * This file is part of WarpX. + * + * Authors: Thomas Marks + * + * License: BSD-3-Clause-LBNL + */ + +#include "Timestep.H" + +#include "WarpX.H" + +#include +#include +#include // TODO: remove this +#include + +#include + +using namespace amrex::literals; + +// constructor +Timestep::Timestep (const std::string& rd_name) +:ReducedDiags{rd_name} +{ + const auto& warpx = WarpX::GetInstance(); + const auto max_level = warpx.maxLevel(); + + // data size should be equal to the number of refinement levels + m_data.resize(max_level + 1, 0.0_rt); + + if (amrex::ParallelDescriptor::IOProcessor() && m_write_header) { + // open file + std::ofstream ofs{m_path + m_rd_name + "." + m_extension, std::ofstream::out}; + + // write header row + int c = 0; + ofs << "#"; + ofs << "[" << c++ << "]step()"; + ofs << m_sep; + ofs << "[" << c++ << "]time(s)"; + ofs << m_sep; + + for (int lev = 0; lev <= max_level; lev++) { + ofs << "[" << c++ << "]timestep[" << lev << "](s)"; + if (lev < max_level) { + ofs << m_sep; + } + } + + // close file + ofs << std::endl; + ofs.close(); + } +} +// end constructor + +// function to get current simulation timestep at all refinement levels +void Timestep::ComputeDiags (int step) { + // Check if diagnostic should be done + if (!m_intervals.contains(step+1)) { return; } + + const auto& warpx = WarpX::GetInstance(); + const auto max_level = warpx.maxLevel(); + const auto dt = warpx.getdt(); + + for (int lev = 0; lev <= max_level; lev++) { + m_data[lev] = dt[lev]; + } +} +// end Timestep::ComputeDiags diff --git a/Source/Evolve/WarpXComputeDt.cpp b/Source/Evolve/WarpXComputeDt.cpp index 3d8d03f8e81..81ff83cf530 100644 --- a/Source/Evolve/WarpXComputeDt.cpp +++ b/Source/Evolve/WarpXComputeDt.cpp @@ -60,7 +60,11 @@ WarpX::ComputeDt () electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) { // Computation of dt for spectral algorithm // (determined by the minimum cell size in all directions) - deltat = cfl / PhysConst::c * minDim(dx); + if (m_max_dt.has_value()) { + deltat = m_max_dt.value(); + } else { + deltat = cfl / PhysConst::c * minDim(dx); + } } else { // Computation of dt for FDTD algorithm #ifdef WARPX_DIM_RZ @@ -106,12 +110,12 @@ WarpX::UpdateDtFromParticleSpeeds () amrex::Real dx_min = minDim(dx); const amrex::ParticleReal max_v = mypc->maxParticleVelocity(); - amrex::Real deltat_new; + amrex::Real deltat_new = 0.; // Protections from overly-large timesteps if (max_v == 0) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_max_dt.has_value(), "Particles at rest and no maximum timestep specified. Aborting."); - deltat_new == m_max_dt.value(); + deltat_new = m_max_dt.value(); } deltat_new = cfl * dx_min / max_v; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 460561fd33c..d0c8fb979eb 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -2931,6 +2931,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, copyAttribs(ip); } + // Push momentum from t^{n-1/2} to t^{n+1/2} doParticleMomentumPush<0>(ux[ip], uy[ip], uz[ip], Exp, Eyp, Ezp, Bxp, Byp, Bzp, ion_lev ? ion_lev[ip] : 1, @@ -2938,9 +2939,10 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, #ifdef WARPX_QED t_chi_max, #endif - dt); + 0.5_rt * (dt + dt_next)); - UpdatePosition(xp, yp, zp, ux[ip], uy[ip], uz[ip], 0.5_rt * (dt + dt_next)); + // Push position from t^n to t^{n+1} + UpdatePosition(xp, yp, zp, ux[ip], uy[ip], uz[ip], dt); setPosition(ip, xp, yp, zp); } #ifdef WARPX_QED From 4eb742324259c042d2dd179d827310c5efe9d633 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Thu, 12 Sep 2024 11:18:44 -0700 Subject: [PATCH 19/32] Fix max dt picmi input --- Python/pywarpx/picmi.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index dc4171b02ea..30cef5fd794 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -2809,6 +2809,8 @@ def init(self, kw): self.checkpoint_signals = kw.pop("warpx_checkpoint_signals", None) self.numprocs = kw.pop("warpx_numprocs", None) + self.max_dt = kw.pop("warpx_max_dt", None) + self.inputs_initialized = False self.warpx_initialized = False @@ -2822,8 +2824,8 @@ def initialize_inputs(self): if self.time_step_size is not None: pywarpx.warpx.const_dt = self.time_step_size - if self.warpx_max_dt is not None: - pywarpx.warpx.max_dt = self.warpx_max_dt + if self.max_dt is not None: + pywarpx.warpx.max_dt = self.max_dt if self.gamma_boost is not None: pywarpx.warpx.gamma_boost = self.gamma_boost From 23685c93575365762774b0f6e38ce05f02a7b6af Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Thu, 12 Sep 2024 11:20:27 -0700 Subject: [PATCH 20/32] Change condition for updating dt --- Source/Evolve/WarpXEvolve.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 571aba59080..b79cbf0ee80 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -88,7 +88,7 @@ WarpX::Evolve (int numsteps) multi_diags->NewIteration(); // Update timestep for electrostatic solver if a constant dt is not provided - if (electrostatic_solver_id != ElectrostaticSolverAlgo::None && + if (electromagnetic_solver_id == ElectromagneticSolverAlgo::None && !m_const_dt.has_value()) { UpdateDtFromParticleSpeeds(); } From af26fcd4933211f8980286045ff6e9fc9fb03688 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Thu, 12 Sep 2024 11:46:23 -0700 Subject: [PATCH 21/32] Fix some more picmi stuff --- Python/pywarpx/picmi.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index 30cef5fd794..534869fa467 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -1522,6 +1522,8 @@ def solver_initialize_inputs(self): # --- Same method names are used, though mapped to lower case. pywarpx.algo.maxwell_solver = self.method + pywarpx.warpx.cfl = self.cfl + if self.source_smoother is not None: self.source_smoother.smoother_initialize_inputs(self) @@ -1869,6 +1871,9 @@ class ElectrostaticSolver(picmistandard.PICMI_ElectrostaticSolver): Parameters ---------- + warpx_cfl: float, optional + Fraction of the CFL condition for particle velocity vs grid size, used to set the timestep + warpx_relativistic: bool, default=False Whether to use the relativistic solver or lab frame solver @@ -1892,8 +1897,8 @@ def solver_initialize_inputs(self): self.grid.grid_initialize_inputs() - if self.cfl is not None: - pywarpx.warpx.cfl = self.cfl + # set CFL number + pywarpx.warpx.cfl = self.cfl if self.relativistic: pywarpx.warpx.do_electrostatic = "relativistic" From eb71b0bdbbaad04f4b32ea9ce6474a9863329703 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Thu, 12 Sep 2024 12:09:34 -0700 Subject: [PATCH 22/32] Add timestep diagnostic to docs --- Docs/source/usage/parameters.rst | 51 +++++++++++++++++--------------- 1 file changed, 27 insertions(+), 24 deletions(-) diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index 1d9c0c14bbf..37b12df3551 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -3448,39 +3448,42 @@ Reduced Diagnostics For 1D-Z, :math:`x`-related and :math:`y`-related quantities are not outputted. RZ geometry is not supported yet. -* ``DifferentialLuminosity`` - This type computes the differential luminosity between two species, defined as: + * ``DifferentialLuminosity`` + This type computes the differential luminosity between two species, defined as: - .. math:: + .. math:: - \frac{d\mathcal{L}}{d\mathcal{E}^*}(\mathcal{E}^*, t) = \int_0^t dt'\int d\boldsymbol{x}\,d\boldsymbol{p}_1 d\boldsymbol{p}_2\; - \sqrt{ |\boldsymbol{v}_1 - \boldsymbol{v}_2|^2 - |\boldsymbol{v}_1\times\boldsymbol{v}_2|^2/c^2} \\ f_1(\boldsymbol{x}, \boldsymbol{p}_1, t')f_2(\boldsymbol{x}, \boldsymbol{p}_2, t') \delta(\mathcal{E}^* - \mathcal{E}^*(\boldsymbol{p}_1, \boldsymbol{p}_2)) + \frac{d\mathcal{L}}{d\mathcal{E}^*}(\mathcal{E}^*, t) = \int_0^t dt'\int d\boldsymbol{x}\,d\boldsymbol{p}_1 d\boldsymbol{p}_2\; + \sqrt{ |\boldsymbol{v}_1 - \boldsymbol{v}_2|^2 - |\boldsymbol{v}_1\times\boldsymbol{v}_2|^2/c^2} \\ f_1(\boldsymbol{x}, \boldsymbol{p}_1, t')f_2(\boldsymbol{x}, \boldsymbol{p}_2, t') \delta(\mathcal{E}^* - \mathcal{E}^*(\boldsymbol{p}_1, \boldsymbol{p}_2)) - where :math:`\mathcal{E}^*(\boldsymbol{p}_1, \boldsymbol{p}_2) = \sqrt{m_1^2c^4 + m_2^2c^4 + 2(m_1 m_2 c^4 - \gamma_1 \gamma_2 - \boldsymbol{p}_1\cdot\boldsymbol{p}_2 c^2)}` is the energy in the center-of-mass frame, - and :math:`f_i` is the distribution function of species :math:`i`. Note that, if :math:`\sigma^*(\mathcal{E}^*)` - is the center-of-mass cross-section of a given collision process, then - :math:`\int d\mathcal{E}^* \frac{d\mathcal{L}}{d\mathcal{E}^*} (\mathcal{E}^*, t)\sigma^*(\mathcal{E}^*)` - gives the total number of collisions of that process (from the beginning of the simulation up until time :math:`t`). + where :math:`\mathcal{E}^*(\boldsymbol{p}_1, \boldsymbol{p}_2) = \sqrt{m_1^2c^4 + m_2^2c^4 + 2(m_1 m_2 c^4 + \gamma_1 \gamma_2 - \boldsymbol{p}_1\cdot\boldsymbol{p}_2 c^2)}` is the energy in the center-of-mass frame, + and :math:`f_i` is the distribution function of species :math:`i`. Note that, if :math:`\sigma^*(\mathcal{E}^*)` + is the center-of-mass cross-section of a given collision process, then + :math:`\int d\mathcal{E}^* \frac{d\mathcal{L}}{d\mathcal{E}^*} (\mathcal{E}^*, t)\sigma^*(\mathcal{E}^*)` + gives the total number of collisions of that process (from the beginning of the simulation up until time :math:`t`). - The differential luminosity is given in units of :math:`\text{m}^{-2}.\text{eV}^{-1}`. For collider-relevant WarpX simulations - involving two crossing, high-energy beams of particles, the differential luminosity in :math:`\text{s}^{-1}.\text{m}^{-2}.\text{eV}^{-1}` - can be obtained by multiplying the above differential luminosity by the expected repetition rate of the beams. + The differential luminosity is given in units of :math:`\text{m}^{-2}.\text{eV}^{-1}`. For collider-relevant WarpX simulations + involving two crossing, high-energy beams of particles, the differential luminosity in :math:`\text{s}^{-1}.\text{m}^{-2}.\text{eV}^{-1}` + can be obtained by multiplying the above differential luminosity by the expected repetition rate of the beams. - In practice, the above expression of the differential luminosity is evaluated over discrete bins in energy :math:`\mathcal{E}^*`, - and by summing over macroparticles. + In practice, the above expression of the differential luminosity is evaluated over discrete bins in energy :math:`\mathcal{E}^*`, + and by summing over macroparticles. - * ``.species`` (`list of two strings`) - The names of the two species for which the differential luminosity is computed. + * ``.species`` (`list of two strings`) + The names of the two species for which the differential luminosity is computed. + + * ``.bin_number`` (`int` > 0) + The number of bins in energy :math:`\mathcal{E}^*` - * ``.bin_number`` (`int` > 0) - The number of bins in energy :math:`\mathcal{E}^*` + * ``.bin_max`` (`float`, in eV) + The minimum value of :math:`\mathcal{E}^*` for which the differential luminosity is computed. - * ``.bin_max`` (`float`, in eV) - The minimum value of :math:`\mathcal{E}^*` for which the differential luminosity is computed. + * ``.bin_min`` (`float`, in eV) + The maximum value of :math:`\mathcal{E}^*` for which the differential luminosity is computed. - * ``.bin_min`` (`float`, in eV) - The maximum value of :math:`\mathcal{E}^*` for which the differential luminosity is computed. + * ``Timestep`` + This type outputs the simulation's physical timestep (in seconds) at each mesh refinement level. * ``.intervals`` (`string`) Using the `Intervals Parser`_ syntax, this string defines the timesteps at which reduced From 5f6df7316eef9d64d497dd295074ccad0db903db Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Thu, 12 Sep 2024 12:41:48 -0700 Subject: [PATCH 23/32] Add timestep diagnostic to picmi --- Python/pywarpx/picmi.py | 1 + 1 file changed, 1 insertion(+) diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index 534869fa467..a888b3c1102 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -3904,6 +3904,7 @@ def __init__( "ParticleNumber", "LoadBalanceCosts", "LoadBalanceEfficiency", + "Timestep", ] # The species diagnostics require a species to be provided self._species_reduced_diagnostics = [ From 3bfb57daebead82106f46ca6a0be2433e89bcec8 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Thu, 12 Sep 2024 13:50:28 -0700 Subject: [PATCH 24/32] Initialize dt_next on restart --- Source/Diagnostics/WarpXIO.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/Source/Diagnostics/WarpXIO.cpp b/Source/Diagnostics/WarpXIO.cpp index a3b902386f6..bba7d24c2c4 100644 --- a/Source/Diagnostics/WarpXIO.cpp +++ b/Source/Diagnostics/WarpXIO.cpp @@ -162,6 +162,11 @@ WarpX::InitFromCheckpoint () lis >> word; dt_lev = static_cast(std::stod(word)); } + // set dt_next + dt_next.resize(dt.size(), dt[0]); + for (int lev = 0; lev <= max_level; lev++) { + dt_next[lev] = dt[lev]; + } } amrex::Real moving_window_x_checkpoint; From bcafbab11df10eec8d9eb393c5d71de53c97e350 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Thu, 12 Sep 2024 15:23:11 -0700 Subject: [PATCH 25/32] Add electrostatic sphere test using adaptive timestep --- .../Tests/electrostatic_sphere/CMakeLists.txt | 10 ++++++++++ .../Tests/electrostatic_sphere/inputs_base_3d | 2 -- .../electrostatic_sphere/inputs_base_constdt_3d | 6 ++++++ .../inputs_test_3d_electrostatic_sphere | 2 +- ...inputs_test_3d_electrostatic_sphere_adaptive | 10 ++++++++++ ...nputs_test_3d_electrostatic_sphere_lab_frame | 2 +- ...d_electrostatic_sphere_lab_frame_mr_emass_10 | 2 +- ...nputs_test_3d_electrostatic_sphere_rel_nodal | 2 +- .../test_3d_electrostatic_sphere_adaptive.json | 17 +++++++++++++++++ 9 files changed, 47 insertions(+), 6 deletions(-) create mode 100644 Examples/Tests/electrostatic_sphere/inputs_base_constdt_3d create mode 100644 Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_adaptive create mode 100644 Regression/Checksum/benchmarks_json/test_3d_electrostatic_sphere_adaptive.json diff --git a/Examples/Tests/electrostatic_sphere/CMakeLists.txt b/Examples/Tests/electrostatic_sphere/CMakeLists.txt index 41a151b7884..2d1682b62bb 100644 --- a/Examples/Tests/electrostatic_sphere/CMakeLists.txt +++ b/Examples/Tests/electrostatic_sphere/CMakeLists.txt @@ -41,6 +41,16 @@ add_warpx_test( OFF # dependency ) +add_warpx_test( + test_3d_electrostatic_sphere_adaptive # name + 3 # dims + 2 # nprocs + inputs_test_3d_electrostatic_sphere_adaptive # inputs + analysis_electrostatic_sphere.py # analysis + diags/diag1000018 # output + OFF # dependency +) + add_warpx_test( test_rz_electrostatic_sphere # name RZ # dims diff --git a/Examples/Tests/electrostatic_sphere/inputs_base_3d b/Examples/Tests/electrostatic_sphere/inputs_base_3d index 1f15b2da0de..3bc302f6225 100644 --- a/Examples/Tests/electrostatic_sphere/inputs_base_3d +++ b/Examples/Tests/electrostatic_sphere/inputs_base_3d @@ -1,4 +1,3 @@ -max_step = 30 amr.n_cell = 64 64 64 amr.max_level = 0 amr.blocking_factor = 8 @@ -8,7 +7,6 @@ geometry.prob_lo = -0.5 -0.5 -0.5 geometry.prob_hi = 0.5 0.5 0.5 boundary.field_lo = pec pec pec boundary.field_hi = pec pec pec -warpx.const_dt = 1e-6 warpx.do_electrostatic = relativistic particles.species_names = electron diff --git a/Examples/Tests/electrostatic_sphere/inputs_base_constdt_3d b/Examples/Tests/electrostatic_sphere/inputs_base_constdt_3d new file mode 100644 index 00000000000..bfa4567e21e --- /dev/null +++ b/Examples/Tests/electrostatic_sphere/inputs_base_constdt_3d @@ -0,0 +1,6 @@ +# base inputs +FILE = inputs_base_3d + +# inputs for tests using constant dt +max_step = 30 +warpx_const_dt = 1e-6 diff --git a/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere b/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere index d89395e9d74..84129349c7a 100644 --- a/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere +++ b/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere @@ -1,5 +1,5 @@ # base input parameters -FILE = inputs_base_3d +FILE = inputs_base_constdt_3d # test input parameters warpx.abort_on_warning_threshold = medium diff --git a/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_adaptive b/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_adaptive new file mode 100644 index 00000000000..99245fc616f --- /dev/null +++ b/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_adaptive @@ -0,0 +1,10 @@ +# base input parameters +FILE = inputs_base_3d + +# adaptive timestepping params +warpx.cfl = 0.2 +warpx.max_dt = 2.0e-6 +stop_time = 30e-6 + +# test input parameters +warpx.abort_on_warning_threshold = medium diff --git a/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_lab_frame b/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_lab_frame index da97ae8afe7..0a44b6a21b0 100644 --- a/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_lab_frame +++ b/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_lab_frame @@ -1,5 +1,5 @@ # base input parameters -FILE = inputs_base_3d +FILE = inputs_base_constdt_3d # test input parameters diag2.electron.variables = x y z ux uy uz w phi diff --git a/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_lab_frame_mr_emass_10 b/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_lab_frame_mr_emass_10 index 481cc65f030..420f6eb1e2d 100644 --- a/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_lab_frame_mr_emass_10 +++ b/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_lab_frame_mr_emass_10 @@ -1,5 +1,5 @@ # base input parameters -FILE = inputs_base_3d +FILE = inputs_base_constdt_3d # test input parameters amr.max_level = 1 diff --git a/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_rel_nodal b/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_rel_nodal index 96bff8aa9c7..9c32802f116 100644 --- a/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_rel_nodal +++ b/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_rel_nodal @@ -1,5 +1,5 @@ # base input parameters -FILE = inputs_base_3d +FILE = inputs_base_constdt_3d # test input parameters warpx.abort_on_warning_threshold = medium diff --git a/Regression/Checksum/benchmarks_json/test_3d_electrostatic_sphere_adaptive.json b/Regression/Checksum/benchmarks_json/test_3d_electrostatic_sphere_adaptive.json new file mode 100644 index 00000000000..2444322b12d --- /dev/null +++ b/Regression/Checksum/benchmarks_json/test_3d_electrostatic_sphere_adaptive.json @@ -0,0 +1,17 @@ +{ + "lev=0": { + "Ex": 6.506693947590733, + "Ey": 6.506693947080617, + "Ez": 6.506693947486778, + "rho": 2.6092567413231304e-10 + }, + "electron": { + "particle_momentum_x": 9.686609803282714e-24, + "particle_momentum_y": 9.686609801152154e-24, + "particle_momentum_z": 9.686609801429006e-24, + "particle_position_x": 512.8986106449738, + "particle_position_y": 512.8986108927056, + "particle_position_z": 512.898610747885, + "particle_weight": 6212.501525878906 + } +} From 813bdad45cdb1604f325c7a42abe64c91ebd7dc4 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Thu, 12 Sep 2024 15:31:10 -0700 Subject: [PATCH 26/32] Fix clang-tidy errors --- Source/Evolve/WarpXComputeDt.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Source/Evolve/WarpXComputeDt.cpp b/Source/Evolve/WarpXComputeDt.cpp index 81ff83cf530..d839de02e5f 100644 --- a/Source/Evolve/WarpXComputeDt.cpp +++ b/Source/Evolve/WarpXComputeDt.cpp @@ -107,7 +107,7 @@ void WarpX::UpdateDtFromParticleSpeeds () { const amrex::Real* dx = geom[max_level].CellSize(); - amrex::Real dx_min = minDim(dx); + const amrex::Real dx_min = minDim(dx); const amrex::ParticleReal max_v = mypc->maxParticleVelocity(); amrex::Real deltat_new = 0.; @@ -116,13 +116,13 @@ WarpX::UpdateDtFromParticleSpeeds () if (max_v == 0) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_max_dt.has_value(), "Particles at rest and no maximum timestep specified. Aborting."); deltat_new = m_max_dt.value(); + } else { + deltat_new = cfl * dx_min / max_v; } - deltat_new = cfl * dx_min / max_v; - // Restrict to be less than user-specified maximum timestep, if present if (m_max_dt.has_value()) { - deltat_new = std::min(deltat_new, m_max_dt.value()); + deltat_new = std::min(deltat_new, cfl * dx_min / max_v); } // Set present dt to previous next dt From 8db80fbad0c1fd8a875429a22c299661cf7b08b8 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Fri, 13 Sep 2024 12:47:01 -0700 Subject: [PATCH 27/32] ruff reformatting --- .../electrostatic_sphere/inputs_base_constdt_3d | 2 +- Python/pywarpx/picmi.py | 17 +++++++---------- Source/Evolve/WarpXComputeDt.cpp | 2 +- 3 files changed, 9 insertions(+), 12 deletions(-) diff --git a/Examples/Tests/electrostatic_sphere/inputs_base_constdt_3d b/Examples/Tests/electrostatic_sphere/inputs_base_constdt_3d index bfa4567e21e..da88ccb8602 100644 --- a/Examples/Tests/electrostatic_sphere/inputs_base_constdt_3d +++ b/Examples/Tests/electrostatic_sphere/inputs_base_constdt_3d @@ -3,4 +3,4 @@ FILE = inputs_base_3d # inputs for tests using constant dt max_step = 30 -warpx_const_dt = 1e-6 +warpx.const_dt = 1e-6 diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index a888b3c1102..c39d87638c0 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -1872,7 +1872,10 @@ class ElectrostaticSolver(picmistandard.PICMI_ElectrostaticSolver): Parameters ---------- warpx_cfl: float, optional - Fraction of the CFL condition for particle velocity vs grid size, used to set the timestep + Fraction of the CFL condition for particle velocity vs grid size, used to set the timestep when employing adaptive timestepping + + warpx_max_dt: float, optional + The maximum allowable timestep when using adaptive timestepping (only available for non-electromagnetic solvers) warpx_relativistic: bool, default=False Whether to use the relativistic solver or lab frame solver @@ -1890,6 +1893,7 @@ def init(self, kw): self.self_fields_verbosity = kw.pop("warpx_self_fields_verbosity", None) self.magnetostatic = kw.pop("warpx_magnetostatic", False) self.cfl = kw.pop("warpx_cfl", None) + self.max_dt = kw.pop("warpx_max_dt", None) def solver_initialize_inputs(self): # Open BC means FieldBoundaryType::Open for electrostatic sims, rather than perfectly-matched layer @@ -1897,8 +1901,9 @@ def solver_initialize_inputs(self): self.grid.grid_initialize_inputs() - # set CFL number + # set adaptive timestepping parameters pywarpx.warpx.cfl = self.cfl + pywarpx.warpx.max_dt = self.max_dt if self.relativistic: pywarpx.warpx.do_electrostatic = "relativistic" @@ -2745,9 +2750,6 @@ class Simulation(picmistandard.PICMI_Simulation): warpx_used_inputs_file: string, optional The name of the text file that the used input parameters is written to, - - warpx_max_dt: float, optional - The maximum allowable timestep when using adaptive timestepping (only available for non-electromagnetic solvers) """ # Set the C++ WarpX interface (see _libwarpx.LibWarpX) as an extension to @@ -2814,8 +2816,6 @@ def init(self, kw): self.checkpoint_signals = kw.pop("warpx_checkpoint_signals", None) self.numprocs = kw.pop("warpx_numprocs", None) - self.max_dt = kw.pop("warpx_max_dt", None) - self.inputs_initialized = False self.warpx_initialized = False @@ -2829,9 +2829,6 @@ def initialize_inputs(self): if self.time_step_size is not None: pywarpx.warpx.const_dt = self.time_step_size - if self.max_dt is not None: - pywarpx.warpx.max_dt = self.max_dt - if self.gamma_boost is not None: pywarpx.warpx.gamma_boost = self.gamma_boost pywarpx.warpx.boost_direction = "z" diff --git a/Source/Evolve/WarpXComputeDt.cpp b/Source/Evolve/WarpXComputeDt.cpp index d839de02e5f..10ebedc6f86 100644 --- a/Source/Evolve/WarpXComputeDt.cpp +++ b/Source/Evolve/WarpXComputeDt.cpp @@ -114,7 +114,7 @@ WarpX::UpdateDtFromParticleSpeeds () // Protections from overly-large timesteps if (max_v == 0) { - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_max_dt.has_value(), "Particles at rest and no maximum timestep specified. Aborting."); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_max_dt.has_value(), "Particles at rest and no constant or maximum timestep specified. Aborting."); deltat_new = m_max_dt.value(); } else { deltat_new = cfl * dx_min / max_v; From 0730b0e209485cae7c5bb53dbd3f4dd70607766c Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Fri, 13 Sep 2024 13:33:12 -0700 Subject: [PATCH 28/32] begin refactoring based on feedback --- Docs/source/usage/parameters.rst | 14 +- Source/Diagnostics/WarpXIO.cpp | 5 - Source/Evolve/WarpXComputeDt.cpp | 21 +- Source/Evolve/WarpXEvolve.cpp | 60 ++- Source/Particles/LaserParticleContainer.H | 2 +- Source/Particles/LaserParticleContainer.cpp | 3 +- Source/Particles/MultiParticleContainer.H | 2 +- Source/Particles/MultiParticleContainer.cpp | 4 +- Source/Particles/PhotonParticleContainer.H | 4 +- Source/Particles/PhotonParticleContainer.cpp | 7 +- Source/Particles/PhysicalParticleContainer.H | 13 +- .../Particles/PhysicalParticleContainer.cpp | 454 +++++------------- .../RigidInjectedParticleContainer.H | 5 +- .../RigidInjectedParticleContainer.cpp | 9 +- Source/Particles/WarpXParticleContainer.H | 2 +- Source/WarpX.H | 3 +- 16 files changed, 190 insertions(+), 418 deletions(-) diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index 37b12df3551..451bc4872ec 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -2129,14 +2129,24 @@ Time step The ratio between the actual timestep that is used in the simulation and the Courant-Friedrichs-Lewy (CFL) limit. (e.g. for `warpx.cfl=1`, the timestep will be exactly equal to the CFL limit.) - This parameter will only be used with the electromagnetic solver. + For some speed v and grid spacing dx, this limits the timestep to `warpx.cfl * dx / v`. + When used with the electromagnetic solver, `v` is the speed of light. + For the electrostatic solver, `v` is the maximum speed among all particles in the domain. * ``warpx.const_dt`` (`float`) Allows direct specification of the time step size, in units of seconds. - When the electrostatic solver is being used, this must be supplied. + When the electrostatic solver is being used, this must be supplied if not using adaptive timestepping. This can be used with the electromagnetic solver, overriding ``warpx.cfl``, but it is up to the user to ensure that the CFL condition is met. +* ``warpx.timestep_adaptation_interval`` (`float`) optional (default `0`) + How many iterations pass between timestep adaptations when using the electrostatic solver. + Must be greater than `0` to use adaptive timestepping, or else ``warpx.const_dt`` must be specified. + +* ``warpx.max_dt`` (`float`) optional + The maximum timestep permitted for the electrostatic solver, when using adaptive timestepping. + If supplied, also sets the initial timestep for these simulations, before the first timestep update. + Filtering ^^^^^^^^^ diff --git a/Source/Diagnostics/WarpXIO.cpp b/Source/Diagnostics/WarpXIO.cpp index bba7d24c2c4..a3b902386f6 100644 --- a/Source/Diagnostics/WarpXIO.cpp +++ b/Source/Diagnostics/WarpXIO.cpp @@ -162,11 +162,6 @@ WarpX::InitFromCheckpoint () lis >> word; dt_lev = static_cast(std::stod(word)); } - // set dt_next - dt_next.resize(dt.size(), dt[0]); - for (int lev = 0; lev <= max_level; lev++) { - dt_next[lev] = dt[lev]; - } } amrex::Real moving_window_x_checkpoint; diff --git a/Source/Evolve/WarpXComputeDt.cpp b/Source/Evolve/WarpXComputeDt.cpp index 10ebedc6f86..7687b6e42e1 100644 --- a/Source/Evolve/WarpXComputeDt.cpp +++ b/Source/Evolve/WarpXComputeDt.cpp @@ -48,6 +48,10 @@ WarpX::ComputeDt () std::stringstream errorMsg; errorMsg << "warpx.const_dt must be specified with the hybrid-PIC solver."; WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_const_dt.has_value(), errorMsg.str()); + } else if (electromagnetic_solver_id == ElectromagneticSolverAlgo::None) { + std::stringstream errorMsg; + errorMsg << "warpx.const_dt must be specified with the electrostatic solver, or the warpx.timestep_adaptation_interval must be > 0."; + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_const_dt.has_value() || timestep_adaptation_interval > 0, errorMsg.str()); } // Determine the appropriate timestep as limited by the speed of light @@ -60,7 +64,7 @@ WarpX::ComputeDt () electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) { // Computation of dt for spectral algorithm // (determined by the minimum cell size in all directions) - if (m_max_dt.has_value()) { + if (m_max_dt.has_value() && electrostatic_solver_id == ElectrostaticSolverAlgo::None) { deltat = m_max_dt.value(); } else { deltat = cfl / PhysConst::c * minDim(dx); @@ -88,13 +92,10 @@ WarpX::ComputeDt () dt.resize(0); dt.resize(max_level+1,deltat); - dt_next.resize(0); - dt_next.resize(max_level+1,deltat); if (do_subcycling) { for (int lev = max_level-1; lev >= 0; --lev) { dt[lev] = dt[lev+1] * refRatio(lev)[0]; - dt_next[lev] = dt[lev]; } } } @@ -125,17 +126,11 @@ WarpX::UpdateDtFromParticleSpeeds () deltat_new = std::min(deltat_new, cfl * dx_min / max_v); } - // Set present dt to previous next dt - dt[max_level] = dt_next[max_level]; - dt_next[max_level] = deltat_new; + // Update dt + dt[max_level] = deltat_new; for (int lev = max_level-1; lev >= 0; --lev) { - dt[lev] = dt_next[lev]; - if (do_subcycling) { - dt_next[lev] = dt_next[lev+1] * refRatio(lev)[0]; - } else { - dt_next[lev] = dt_next[lev+1]; - } + dt[lev] = dt[lev+1] * refRatio(lev)[0]; } } diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index b79cbf0ee80..526eff6e26f 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -60,6 +60,27 @@ using namespace amrex; using ablastr::utils::SignalHandling; +void +WarpX::Synchronize () { + FillBoundaryE(guard_cells.ng_FieldGather); + FillBoundaryB(guard_cells.ng_FieldGather); + if (fft_do_time_averaging) + { + FillBoundaryE_avg(guard_cells.ng_FieldGather); + FillBoundaryB_avg(guard_cells.ng_FieldGather); + } + UpdateAuxilaryData(); + FillBoundaryAux(guard_cells.ng_UpdateAux); + for (int lev = 0; lev <= finest_level; ++lev) { + mypc->PushP(lev, 0.5_rt*dt[lev], + *Efield_aux[lev][0],*Efield_aux[lev][1], + *Efield_aux[lev][2], + *Bfield_aux[lev][0],*Bfield_aux[lev][1], + *Bfield_aux[lev][2]); + } + is_synchronized = true; +} + void WarpX::Evolve (int numsteps) { @@ -87,12 +108,6 @@ WarpX::Evolve (int numsteps) multi_diags->NewIteration(); - // Update timestep for electrostatic solver if a constant dt is not provided - if (electromagnetic_solver_id == ElectromagneticSolverAlgo::None && - !m_const_dt.has_value()) { - UpdateDtFromParticleSpeeds(); - } - // Start loop on time steps if (verbose) { amrex::Print() << "STEP " << step+1 << " starts ...\n"; @@ -101,6 +116,15 @@ WarpX::Evolve (int numsteps) CheckLoadBalance(step); + // Update timestep for electrostatic solver if a constant dt is not provided + // This first synchronizes the position and velocity before setting the new timestep + if (electromagnetic_solver_id == ElectromagneticSolverAlgo::None && + !m_const_dt.has_value() && (step % timestep_adaptation_interval == 0)) { + Synchronize(); + UpdateDtFromParticleSpeeds(); + } + + // If position and velocity are synchronized, push velocity backward one half step if (evolve_scheme == EvolveScheme::Explicit) { ExplicitFillBoundaryEBUpdateAux(); @@ -181,25 +205,9 @@ WarpX::Evolve (int numsteps) // TODO: move out if (evolve_scheme == EvolveScheme::Explicit) { + // At the end of last step, push p by 0.5*dt to synchronize if (cur_time + dt[0] >= stop_time - 1.e-3*dt[0] || step == numsteps_max-1) { - // At the end of last step, push p by 0.5*dt to synchronize - FillBoundaryE(guard_cells.ng_FieldGather); - FillBoundaryB(guard_cells.ng_FieldGather); - if (fft_do_time_averaging) - { - FillBoundaryE_avg(guard_cells.ng_FieldGather); - FillBoundaryB_avg(guard_cells.ng_FieldGather); - } - UpdateAuxilaryData(); - FillBoundaryAux(guard_cells.ng_UpdateAux); - for (int lev = 0; lev <= finest_level; ++lev) { - mypc->PushP(lev, 0.5_rt*dt[lev], - *Efield_aux[lev][0],*Efield_aux[lev][1], - *Efield_aux[lev][2], - *Bfield_aux[lev][0],*Bfield_aux[lev][1], - *Bfield_aux[lev][2]); - } - is_synchronized = true; + Synchronize(); } } @@ -451,7 +459,7 @@ void WarpX::checkEarlyUnusedParams () void WarpX::ExplicitFillBoundaryEBUpdateAux () { WARPX_ALWAYS_ASSERT_WITH_MESSAGE(evolve_scheme == EvolveScheme::Explicit, - "Cannot call WarpX::ExplicitFillBoundaryEBUpdateAux wihtout Explicit evolve scheme set!"); + "Cannot call WarpX::ExplicitFillBoundaryEBUpdateAux without Explicit evolve scheme set!"); // At the beginning, we have B^{n} and E^{n}. // Particles have p^{n} and x^{n}. @@ -1038,7 +1046,7 @@ WarpX::PushParticlesandDeposit (int lev, amrex::Real cur_time, DtType a_dt_type, rho_fp[lev].get(), charge_buf[lev].get(), Efield_cax[lev][0].get(), Efield_cax[lev][1].get(), Efield_cax[lev][2].get(), Bfield_cax[lev][0].get(), Bfield_cax[lev][1].get(), Bfield_cax[lev][2].get(), - cur_time, dt[lev], dt_next[lev], a_dt_type, skip_current, push_type); + cur_time, dt[lev], a_dt_type, skip_current, push_type); if (! skip_current) { #ifdef WARPX_DIM_RZ // This is called after all particles have deposited their current and charge. diff --git a/Source/Particles/LaserParticleContainer.H b/Source/Particles/LaserParticleContainer.H index 42b38113ac2..197cb897602 100644 --- a/Source/Particles/LaserParticleContainer.H +++ b/Source/Particles/LaserParticleContainer.H @@ -72,7 +72,7 @@ public: amrex::MultiFab* rho, amrex::MultiFab* crho, const amrex::MultiFab*, const amrex::MultiFab*, const amrex::MultiFab*, const amrex::MultiFab*, const amrex::MultiFab*, const amrex::MultiFab*, - amrex::Real t, amrex::Real dt, amrex::Real dt_next, DtType a_dt_type=DtType::Full, + amrex::Real t, amrex::Real dt, DtType a_dt_type=DtType::Full, bool skip_deposition=false, PushType push_type=PushType::Explicit) final; void PushP (int lev, amrex::Real dt, diff --git a/Source/Particles/LaserParticleContainer.cpp b/Source/Particles/LaserParticleContainer.cpp index bba38aebfe7..d43fb240756 100644 --- a/Source/Particles/LaserParticleContainer.cpp +++ b/Source/Particles/LaserParticleContainer.cpp @@ -564,8 +564,7 @@ LaserParticleContainer::Evolve (int lev, MultiFab* rho, MultiFab* crho, const MultiFab*, const MultiFab*, const MultiFab*, const MultiFab*, const MultiFab*, const MultiFab*, - Real t, Real dt, Real /*dt_next*/, DtType /*a_dt_type*/, - bool skip_deposition, PushType push_type) + Real t, Real dt, DtType /*a_dt_type*/, bool skip_deposition, PushType push_type) { WARPX_PROFILE("LaserParticleContainer::Evolve()"); WARPX_PROFILE_VAR_NS("LaserParticleContainer::Evolve::ParticlePush", blp_pp); diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index ce84dcb7460..97e4e1bc4da 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -110,7 +110,7 @@ public: amrex::MultiFab* rho, amrex::MultiFab* crho, const amrex::MultiFab* cEx, const amrex::MultiFab* cEy, const amrex::MultiFab* cEz, const amrex::MultiFab* cBx, const amrex::MultiFab* cBy, const amrex::MultiFab* cBz, - amrex::Real t, amrex::Real dt, amrex::Real dt_next, DtType a_dt_type=DtType::Full, bool skip_deposition=false, + amrex::Real t, amrex::Real dt, DtType a_dt_type=DtType::Full, bool skip_deposition=false, PushType push_type=PushType::Explicit); /** diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index e33df518115..23af4177228 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -465,7 +465,7 @@ MultiParticleContainer::Evolve (int lev, MultiFab* rho, MultiFab* crho, const MultiFab* cEx, const MultiFab* cEy, const MultiFab* cEz, const MultiFab* cBx, const MultiFab* cBy, const MultiFab* cBz, - Real t, Real dt, Real dt_next, DtType a_dt_type, bool skip_deposition, + Real t, Real dt, DtType a_dt_type, bool skip_deposition, PushType push_type) { if (! skip_deposition) { @@ -480,7 +480,7 @@ MultiParticleContainer::Evolve (int lev, } for (auto& pc : allcontainers) { pc->Evolve(lev, Ex, Ey, Ez, Bx, By, Bz, jx, jy, jz, cjx, cjy, cjz, - rho, crho, cEx, cEy, cEz, cBx, cBy, cBz, t, dt, dt_next, a_dt_type, skip_deposition, push_type); + rho, crho, cEx, cEy, cEz, cBx, cBy, cBz, t, dt, a_dt_type, skip_deposition, push_type); } } diff --git a/Source/Particles/PhotonParticleContainer.H b/Source/Particles/PhotonParticleContainer.H index e0be5076909..34afac53482 100644 --- a/Source/Particles/PhotonParticleContainer.H +++ b/Source/Particles/PhotonParticleContainer.H @@ -69,7 +69,6 @@ public: const amrex::MultiFab* cBz, amrex::Real t, amrex::Real dt, - amrex::Real dt_next, DtType a_dt_type=DtType::Full, bool skip_deposition=false, PushType push_type=PushType::Explicit) override; @@ -85,8 +84,7 @@ public: long offset, long np_to_push, int lev, int gather_lev, - amrex::Real dt, amrex::Real dt_next, - ScaleFields scaleFields, + amrex::Real dt, ScaleFields scaleFields, DtType a_dt_type) override; // Do nothing diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index f71a7d36595..1f15d5210f5 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -91,8 +91,7 @@ PhotonParticleContainer::PushPX (WarpXParIter& pti, const long offset, const long np_to_push, int lev, int gather_lev, - amrex::Real dt, amrex::Real /*dt_next*/, - ScaleFields /*scaleFields*/, DtType a_dt_type) + amrex::Real dt, ScaleFields /*scaleFields*/, DtType a_dt_type) { // Get inverse cell size on gather_lev const amrex::XDim3 dinv = WarpX::InvCellSize(std::max(gather_lev,0)); @@ -238,7 +237,7 @@ PhotonParticleContainer::Evolve (int lev, MultiFab* rho, MultiFab* crho, const MultiFab* cEx, const MultiFab* cEy, const MultiFab* cEz, const MultiFab* cBx, const MultiFab* cBy, const MultiFab* cBz, - Real t, Real dt, Real dt_next, DtType a_dt_type, bool skip_deposition, + Real t, Real dt, DtType a_dt_type, bool skip_deposition, PushType push_type) { // This does gather, push and deposit. @@ -251,6 +250,6 @@ PhotonParticleContainer::Evolve (int lev, rho, crho, cEx, cEy, cEz, cBx, cBy, cBz, - t, dt, dt_next, a_dt_type, skip_deposition, push_type); + t, dt, a_dt_type, skip_deposition, push_type); } diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index a04e00a8e05..8102fc96a91 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -135,7 +135,6 @@ public: const amrex::MultiFab* cBz, amrex::Real t, amrex::Real dt, - amrex::Real dt_next, DtType a_dt_type=DtType::Full, bool skip_deposition=false, PushType push_type=PushType::Explicit) override; @@ -151,9 +150,7 @@ public: long offset, long np_to_push, int lev, int gather_lev, - amrex::Real dt, - amrex::Real dt_next, - ScaleFields scaleFields, + amrex::Real dt, ScaleFields scaleFields, DtType a_dt_type=DtType::Full); void ImplicitPushXP (WarpXParIter& pti, @@ -395,6 +392,14 @@ public: } protected: + + /* + Finds the box defining the region where refine injection should be used, if that + option is enabled. Currently this only works for numLevels() == 2 and static mesh + refinement. + */ + bool findRefinedInjectionBox (amrex::Box& fine_injection_box, amrex::IntVect& rrfac); + std::string species_name; std::vector> plasma_injectors; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 5d3d31ec14f..d1a19f06993 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -15,6 +15,7 @@ #include "Initialization/InjectorMomentum.H" #include "Initialization/InjectorPosition.H" #include "MultiParticleContainer.H" +#include "Particles/AddPlasmaUtilities.H" #ifdef WARPX_QED # include "Particles/ElementaryProcess/QEDInternals/BreitWheelerEngineWrapper.H" # include "Particles/ElementaryProcess/QEDInternals/QuantumSyncEngineWrapper.H" @@ -217,8 +218,7 @@ namespace const GpuArray& pa, long& ip, const bool& do_field_ionization, int* pi #ifdef WARPX_QED - ,const bool& has_quantum_sync, amrex::ParticleReal* AMREX_RESTRICT p_optical_depth_QSR - ,const bool& has_breit_wheeler, amrex::ParticleReal* AMREX_RESTRICT p_optical_depth_BW + ,const QEDHelper& qed_helper #endif ) noexcept { @@ -227,8 +227,8 @@ namespace } if (do_field_ionization) {pi[ip] = 0;} #ifdef WARPX_QED - if (has_quantum_sync) {p_optical_depth_QSR[ip] = 0._rt;} - if (has_breit_wheeler) {p_optical_depth_BW[ip] = 0._rt;} + if (qed_helper.has_quantum_sync) {qed_helper.p_optical_depth_QSR[ip] = 0._rt;} + if (qed_helper.has_breit_wheeler) {qed_helper.p_optical_depth_BW[ip] = 0._rt;} #endif idcpu[ip] = amrex::ParticleIdCpus::Invalid; @@ -964,22 +964,9 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int amrex::LayoutData* cost = WarpX::getCosts(lev); - const int nlevs = numLevels(); - static bool refine_injection = false; - static Box fine_injection_box; - static amrex::IntVect rrfac(AMREX_D_DECL(1,1,1)); - // This does not work if the mesh is dynamic. But in that case, we should - // not use refined injected either. We also assume there is only one fine level. - if (WarpX::moving_window_active(WarpX::GetInstance().getistep(0)+1) and WarpX::refine_plasma - and do_continuous_injection and nlevs == 2) - { - refine_injection = true; - fine_injection_box = ParticleBoxArray(1).minimalBox(); - fine_injection_box.setSmall(WarpX::moving_window_dir, std::numeric_limits::lowest()/2); - fine_injection_box.setBig(WarpX::moving_window_dir, std::numeric_limits::max()/2); - rrfac = m_gdb->refRatio(0); - fine_injection_box.coarsen(rrfac); - } + Box fine_injection_box; + amrex::IntVect rrfac(AMREX_D_DECL(1,1,1)); + const bool refine_injection = findRefinedInjectionBox(fine_injection_box, rrfac); InjectorPosition* inj_pos = plasma_injector.getInjectorPosition(); InjectorDensity* inj_rho = plasma_injector.getInjectorDensity(); @@ -995,18 +982,12 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int const bool radially_weighted = plasma_injector.radially_weighted; #endif - - // User-defined integer and real attributes: prepare parsers - const auto n_user_int_attribs = static_cast(m_user_int_attribs.size()); - const auto n_user_real_attribs = static_cast(m_user_real_attribs.size()); - amrex::Gpu::PinnedVector< amrex::ParserExecutor<7> > user_int_attrib_parserexec_pinned(n_user_int_attribs); - amrex::Gpu::PinnedVector< amrex::ParserExecutor<7> > user_real_attrib_parserexec_pinned(n_user_real_attribs); - for (int ia = 0; ia < n_user_int_attribs; ++ia) { - user_int_attrib_parserexec_pinned[ia] = m_user_int_attrib_parser[ia]->compile<7>(); - } - for (int ia = 0; ia < n_user_real_attribs; ++ia) { - user_real_attrib_parserexec_pinned[ia] = m_user_real_attrib_parser[ia]->compile<7>(); - } + auto n_user_int_attribs = static_cast(m_user_int_attribs.size()); + auto n_user_real_attribs = static_cast(m_user_real_attribs.size()); + const PlasmaParserWrapper plasma_parser_wrapper (m_user_int_attribs.size(), + m_user_real_attribs.size(), + m_user_int_attrib_parser, + m_user_real_attrib_parser); MFItInfo info; if (do_tiling && Gpu::notInLaunchRegion()) { @@ -1032,30 +1013,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int RealBox overlap_realbox; Box overlap_box; IntVect shifted; - bool no_overlap = false; - - for (int dir=0; dir= part_realbox.lo(dir) ) { - const Real ncells_adjust = std::floor( (part_realbox.hi(dir) - tile_realbox.hi(dir))/dx[dir] ); - overlap_realbox.setHi( dir, part_realbox.hi(dir) - std::max(ncells_adjust, 0._rt) * dx[dir]); - } else { - no_overlap = true; break; - } - // Count the number of cells in this direction in overlap_realbox - overlap_box.setSmall( dir, 0 ); - overlap_box.setBig( dir, - int( std::round((overlap_realbox.hi(dir)-overlap_realbox.lo(dir)) - /dx[dir] )) - 1); - shifted[dir] = - static_cast(std::round((overlap_realbox.lo(dir)-problo[dir])/dx[dir])); - // shifted is exact in non-moving-window direction. That's all we care. - } + const bool no_overlap = find_overlap(tile_realbox, part_realbox, dx, problo, overlap_realbox, overlap_box, shifted); if (no_overlap) { continue; // Go to the next tile } @@ -1110,19 +1068,9 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int } return 0; }; - const int flag_pcount = checker(); - if (flag_pcount == 1) { - pcounts[index] = num_ppc*r; - } else { - pcounts[index] = 0; - } + pcounts[index] = checker() ? num_ppc*r : 0; } -#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - amrex::ignore_unused(k); -#endif -#if defined(WARPX_DIM_1D_Z) amrex::ignore_unused(j,k); -#endif }); // Max number of new particles. All of them are created, @@ -1160,40 +1108,12 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int pa[ia] = soa.GetRealData(ia).data() + old_size; } uint64_t * AMREX_RESTRICT pa_idcpu = soa.GetIdCPUData().data() + old_size; - // user-defined integer and real attributes - amrex::Gpu::PinnedVector pa_user_int_pinned(n_user_int_attribs); - amrex::Gpu::PinnedVector pa_user_real_pinned(n_user_real_attribs); - for (int ia = 0; ia < n_user_int_attribs; ++ia) { - pa_user_int_pinned[ia] = soa.GetIntData(particle_icomps[m_user_int_attribs[ia]]).data() + old_size; - } - for (int ia = 0; ia < n_user_real_attribs; ++ia) { - pa_user_real_pinned[ia] = soa.GetRealData(particle_comps[m_user_real_attribs[ia]]).data() + old_size; - } -#ifdef AMREX_USE_GPU - // To avoid using managed memory, we first define pinned memory vector, initialize on cpu, - // and them memcpy to device from host - amrex::Gpu::DeviceVector d_pa_user_int(n_user_int_attribs); - amrex::Gpu::DeviceVector d_pa_user_real(n_user_real_attribs); - amrex::Gpu::DeviceVector< amrex::ParserExecutor<7> > d_user_int_attrib_parserexec(n_user_int_attribs); - amrex::Gpu::DeviceVector< amrex::ParserExecutor<7> > d_user_real_attrib_parserexec(n_user_real_attribs); - amrex::Gpu::copyAsync(Gpu::hostToDevice, pa_user_int_pinned.begin(), - pa_user_int_pinned.end(), d_pa_user_int.begin()); - amrex::Gpu::copyAsync(Gpu::hostToDevice, pa_user_real_pinned.begin(), - pa_user_real_pinned.end(), d_pa_user_real.begin()); - amrex::Gpu::copyAsync(Gpu::hostToDevice, user_int_attrib_parserexec_pinned.begin(), - user_int_attrib_parserexec_pinned.end(), d_user_int_attrib_parserexec.begin()); - amrex::Gpu::copyAsync(Gpu::hostToDevice, user_real_attrib_parserexec_pinned.begin(), - user_real_attrib_parserexec_pinned.end(), d_user_real_attrib_parserexec.begin()); - int** pa_user_int_data = d_pa_user_int.dataPtr(); - ParticleReal** pa_user_real_data = d_pa_user_real.dataPtr(); - amrex::ParserExecutor<7> const* user_int_parserexec_data = d_user_int_attrib_parserexec.dataPtr(); - amrex::ParserExecutor<7> const* user_real_parserexec_data = d_user_real_attrib_parserexec.dataPtr(); -#else - int** pa_user_int_data = pa_user_int_pinned.dataPtr(); - ParticleReal** pa_user_real_data = pa_user_real_pinned.dataPtr(); - amrex::ParserExecutor<7> const* user_int_parserexec_data = user_int_attrib_parserexec_pinned.dataPtr(); - amrex::ParserExecutor<7> const* user_real_parserexec_data = user_real_attrib_parserexec_pinned.dataPtr(); -#endif + + PlasmaParserHelper plasma_parser_helper (soa, old_size, m_user_int_attribs, m_user_real_attribs, particle_icomps, particle_comps, plasma_parser_wrapper); + int** pa_user_int_data = plasma_parser_helper.getUserIntDataPtrs(); + ParticleReal** pa_user_real_data = plasma_parser_helper.getUserRealDataPtrs(); + amrex::ParserExecutor<7> const* user_int_parserexec_data = plasma_parser_helper.getUserIntParserExecData(); + amrex::ParserExecutor<7> const* user_real_parserexec_data = plasma_parser_helper.getUserRealParserExecData(); int* pi = nullptr; if (do_field_ionization) { @@ -1201,34 +1121,9 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int } #ifdef WARPX_QED - //Pointer to the optical depth component - amrex::ParticleReal* p_optical_depth_QSR = nullptr; - amrex::ParticleReal* p_optical_depth_BW = nullptr; - - // If a QED effect is enabled, the corresponding optical depth - // has to be initialized - const bool loc_has_quantum_sync = has_quantum_sync(); - const bool loc_has_breit_wheeler = has_breit_wheeler(); - if (loc_has_quantum_sync) { - p_optical_depth_QSR = soa.GetRealData( - particle_comps["opticalDepthQSR"]).data() + old_size; - } - if(loc_has_breit_wheeler) { - p_optical_depth_BW = soa.GetRealData( - particle_comps["opticalDepthBW"]).data() + old_size; - } - - //If needed, get the appropriate functors from the engines - QuantumSynchrotronGetOpticalDepth quantum_sync_get_opt; - BreitWheelerGetOpticalDepth breit_wheeler_get_opt; - if(loc_has_quantum_sync){ - quantum_sync_get_opt = - m_shr_p_qs_engine->build_optical_depth_functor(); - } - if(loc_has_breit_wheeler){ - breit_wheeler_get_opt = - m_shr_p_bw_engine->build_optical_depth_functor(); - } + const QEDHelper qed_helper(soa, old_size, particle_comps, + has_quantum_sync(), has_breit_wheeler(), + m_shr_p_qs_engine, m_shr_p_bw_engine); #endif const bool loc_do_field_ionization = do_field_ionization; @@ -1246,18 +1141,14 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int [=] AMREX_GPU_DEVICE (int i, int j, int k, amrex::RandomEngine const& engine) noexcept { const IntVect iv = IntVect(AMREX_D_DECL(i, j, k)); + amrex::ignore_unused(j,k); const auto index = overlap_box.index(iv); #ifdef WARPX_DIM_RZ Real theta_offset = 0._rt; if (rz_random_theta) { theta_offset = amrex::Random(engine) * 2._rt * MathConst::pi; } #endif - Real scale_fac = 0.0_rt; - if( pcounts[index] != 0) { - amrex::Real const dV = AMREX_D_TERM(dx[0], *dx[1], *dx[2]); - scale_fac = dV/pcounts[index]; - } - + const Real scale_fac = compute_scale_fac_volume(dx, pcounts[index]); for (int i_part = 0; i_part < pcounts[index]; ++i_part) { long ip = poffset[index] + i_part; @@ -1281,8 +1172,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int if (!box_contains) { ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi #ifdef WARPX_QED - ,loc_has_quantum_sync, p_optical_depth_QSR - ,loc_has_breit_wheeler, p_optical_depth_BW + ,qed_helper #endif ); continue; @@ -1319,8 +1209,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int if (!inj_pos->insideBounds(xb, yb, z0)) { ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi #ifdef WARPX_QED - ,loc_has_quantum_sync, p_optical_depth_QSR - ,loc_has_breit_wheeler, p_optical_depth_BW + ,qed_helper #endif ); continue; @@ -1333,8 +1222,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int if ( dens < density_min ){ ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi #ifdef WARPX_QED - ,loc_has_quantum_sync, p_optical_depth_QSR - ,loc_has_breit_wheeler, p_optical_depth_BW + ,qed_helper #endif ); continue; @@ -1351,8 +1239,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int if (!inj_pos->insideBounds(xb, yb, z0_lab)) { ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi #ifdef WARPX_QED - ,loc_has_quantum_sync, p_optical_depth_QSR - ,loc_has_breit_wheeler, p_optical_depth_BW + ,qed_helper #endif ); continue; @@ -1363,8 +1250,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int if ( dens < density_min ){ ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi #ifdef WARPX_QED - ,loc_has_quantum_sync, p_optical_depth_QSR - ,loc_has_breit_wheeler, p_optical_depth_BW + ,qed_helper #endif ); continue; @@ -1388,12 +1274,12 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int } #ifdef WARPX_QED - if(loc_has_quantum_sync){ - p_optical_depth_QSR[ip] = quantum_sync_get_opt(engine); + if(qed_helper.has_quantum_sync){ + qed_helper.p_optical_depth_QSR[ip] = qed_helper.quantum_sync_get_opt(engine); } - if(loc_has_breit_wheeler){ - p_optical_depth_BW[ip] = breit_wheeler_get_opt(engine); + if(qed_helper.has_breit_wheeler){ + qed_helper.p_optical_depth_BW[ip] = qed_helper.breit_wheeler_get_opt(engine); } #endif // Initialize user-defined integers with user-defined parser @@ -1481,25 +1367,6 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, const auto dx = geom.CellSizeArray(); const auto problo = geom.ProbLoArray(); - Real scale_fac = 0._rt; - // Scale particle weight by the area of the emitting surface, within one cell -#if defined(WARPX_DIM_3D) - scale_fac = dx[0]*dx[1]*dx[2]/dx[plasma_injector.flux_normal_axis]/num_ppc_real; -#elif defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ) - scale_fac = dx[0]*dx[1]/num_ppc_real; - // When emission is in the r direction, the emitting surface is a cylinder. - // The factor 2*pi*r is added later below. - if (plasma_injector.flux_normal_axis == 0) { scale_fac /= dx[0]; } - // When emission is in the z direction, the emitting surface is an annulus - // The factor 2*pi*r is added later below. - if (plasma_injector.flux_normal_axis == 2) { scale_fac /= dx[1]; } - // When emission is in the theta direction (flux_normal_axis == 1), - // the emitting surface is a rectangle, within the plane of the simulation -#elif defined(WARPX_DIM_1D_Z) - scale_fac = dx[0]/num_ppc_real; - if (plasma_injector.flux_normal_axis == 2) { scale_fac /= dx[0]; } -#endif - amrex::LayoutData* cost = WarpX::getCosts(0); // Create temporary particle container to which particles will be added; @@ -1510,19 +1377,9 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, for (int ic = 0; ic < NumRuntimeIntComps(); ++ic) { tmp_pc.AddIntComp(false); } tmp_pc.defineAllParticleTiles(); - const int nlevs = numLevels(); - static bool refine_injection = false; - static Box fine_injection_box; - static amrex::IntVect rrfac(AMREX_D_DECL(1,1,1)); - // This does not work if the mesh is dynamic. But in that case, we should - // not use refined injected either. We also assume there is only one fine level. - if (WarpX::refine_plasma && nlevs == 2) - { - refine_injection = true; - fine_injection_box = ParticleBoxArray(1).minimalBox(); - rrfac = m_gdb->refRatio(0); - fine_injection_box.coarsen(rrfac); - } + Box fine_injection_box; + amrex::IntVect rrfac(AMREX_D_DECL(1,1,1)); + const bool refine_injection = findRefinedInjectionBox(fine_injection_box, rrfac); InjectorPosition* flux_pos = plasma_injector.getInjectorFluxPosition(); InjectorFlux* inj_flux = plasma_injector.getInjectorFlux(); @@ -1536,6 +1393,13 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, const bool radially_weighted = plasma_injector.radially_weighted; #endif + auto n_user_int_attribs = static_cast(m_user_int_attribs.size()); + auto n_user_real_attribs = static_cast(m_user_real_attribs.size()); + const PlasmaParserWrapper plasma_parser_wrapper (m_user_int_attribs.size(), + m_user_real_attribs.size(), + m_user_int_attrib_parser, + m_user_real_attrib_parser); + MFItInfo info; if (do_tiling && Gpu::notInLaunchRegion()) { info.EnableTiling(tile_size); @@ -1560,61 +1424,7 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, RealBox overlap_realbox; Box overlap_box; IntVect shifted; - bool no_overlap = false; - - for (int dir=0; dir 0) { - if (plasma_injector.surface_flux_pos < tile_realbox.lo(dir) || - plasma_injector.surface_flux_pos >= tile_realbox.hi(dir)) { - no_overlap = true; - break; - } - } else { - if (plasma_injector.surface_flux_pos <= tile_realbox.lo(dir) || - plasma_injector.surface_flux_pos > tile_realbox.hi(dir)) { - no_overlap = true; - break; - } - } - overlap_realbox.setLo( dir, plasma_injector.surface_flux_pos ); - overlap_realbox.setHi( dir, plasma_injector.surface_flux_pos ); - overlap_box.setSmall( dir, 0 ); - overlap_box.setBig( dir, 0 ); - shifted[dir] = - static_cast(std::round((overlap_realbox.lo(dir)-problo[dir])/dx[dir])); - } else { - if ( tile_realbox.lo(dir) <= part_realbox.hi(dir) ) { - const Real ncells_adjust = std::floor( (tile_realbox.lo(dir) - part_realbox.lo(dir))/dx[dir] ); - overlap_realbox.setLo( dir, part_realbox.lo(dir) + std::max(ncells_adjust, 0._rt) * dx[dir]); - } else { - no_overlap = true; break; - } - if ( tile_realbox.hi(dir) >= part_realbox.lo(dir) ) { - const Real ncells_adjust = std::floor( (part_realbox.hi(dir) - tile_realbox.hi(dir))/dx[dir] ); - overlap_realbox.setHi( dir, part_realbox.hi(dir) - std::max(ncells_adjust, 0._rt) * dx[dir]); - } else { - no_overlap = true; break; - } - // Count the number of cells in this direction in overlap_realbox - overlap_box.setSmall( dir, 0 ); - overlap_box.setBig( dir, - int( std::round((overlap_realbox.hi(dir)-overlap_realbox.lo(dir)) - /dx[dir] )) - 1); - shifted[dir] = - static_cast(std::round((overlap_realbox.lo(dir)-problo[dir])/dx[dir])); - // shifted is exact in non-moving-window direction. That's all we care. - } - } + const bool no_overlap = find_overlap_flux(tile_realbox, part_realbox, dx, problo, plasma_injector, overlap_realbox, overlap_box, shifted); if (no_overlap) { continue; // Go to the next tile } @@ -1632,6 +1442,7 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, Gpu::DeviceVector offset(overlap_box.numPts()); auto *pcounts = counts.data(); const amrex::IntVect lrrfac = rrfac; + const int flux_normal_axis = plasma_injector.flux_normal_axis; Box fine_overlap_box; // default Box is NOT ok(). if (refine_injection) { fine_overlap_box = overlap_box & amrex::shift(fine_injection_box, -shifted); @@ -1649,17 +1460,13 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, auto index = overlap_box.index(iv); int r; if (fine_overlap_box.ok() && fine_overlap_box.contains(iv)) { - r = AMREX_D_TERM(lrrfac[0],*lrrfac[1],*lrrfac[2]); + r = compute_area_weights(lrrfac, flux_normal_axis); } else { r = 1; } pcounts[index] = num_ppc_int*r; } -#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - amrex::ignore_unused(k); -#elif defined(WARPX_DIM_1D_Z) amrex::ignore_unused(j,k); -#endif }); // Max number of new particles. All of them are created, @@ -1694,46 +1501,11 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, } uint64_t * AMREX_RESTRICT pa_idcpu = soa.GetIdCPUData().data() + old_size; - // user-defined integer and real attributes - const auto n_user_int_attribs = static_cast(m_user_int_attribs.size()); - const auto n_user_real_attribs = static_cast(m_user_real_attribs.size()); - amrex::Gpu::PinnedVector pa_user_int_pinned(n_user_int_attribs); - amrex::Gpu::PinnedVector pa_user_real_pinned(n_user_real_attribs); - amrex::Gpu::PinnedVector< amrex::ParserExecutor<7> > user_int_attrib_parserexec_pinned(n_user_int_attribs); - amrex::Gpu::PinnedVector< amrex::ParserExecutor<7> > user_real_attrib_parserexec_pinned(n_user_real_attribs); - for (int ia = 0; ia < n_user_int_attribs; ++ia) { - pa_user_int_pinned[ia] = soa.GetIntData(particle_icomps[m_user_int_attribs[ia]]).data() + old_size; - user_int_attrib_parserexec_pinned[ia] = m_user_int_attrib_parser[ia]->compile<7>(); - } - for (int ia = 0; ia < n_user_real_attribs; ++ia) { - pa_user_real_pinned[ia] = soa.GetRealData(particle_comps[m_user_real_attribs[ia]]).data() + old_size; - user_real_attrib_parserexec_pinned[ia] = m_user_real_attrib_parser[ia]->compile<7>(); - } -#ifdef AMREX_USE_GPU - // To avoid using managed memory, we first define pinned memory vector, initialize on cpu, - // and them memcpy to device from host - amrex::Gpu::DeviceVector d_pa_user_int(n_user_int_attribs); - amrex::Gpu::DeviceVector d_pa_user_real(n_user_real_attribs); - amrex::Gpu::DeviceVector< amrex::ParserExecutor<7> > d_user_int_attrib_parserexec(n_user_int_attribs); - amrex::Gpu::DeviceVector< amrex::ParserExecutor<7> > d_user_real_attrib_parserexec(n_user_real_attribs); - amrex::Gpu::copyAsync(Gpu::hostToDevice, pa_user_int_pinned.begin(), - pa_user_int_pinned.end(), d_pa_user_int.begin()); - amrex::Gpu::copyAsync(Gpu::hostToDevice, pa_user_real_pinned.begin(), - pa_user_real_pinned.end(), d_pa_user_real.begin()); - amrex::Gpu::copyAsync(Gpu::hostToDevice, user_int_attrib_parserexec_pinned.begin(), - user_int_attrib_parserexec_pinned.end(), d_user_int_attrib_parserexec.begin()); - amrex::Gpu::copyAsync(Gpu::hostToDevice, user_real_attrib_parserexec_pinned.begin(), - user_real_attrib_parserexec_pinned.end(), d_user_real_attrib_parserexec.begin()); - int** pa_user_int_data = d_pa_user_int.dataPtr(); - ParticleReal** pa_user_real_data = d_pa_user_real.dataPtr(); - amrex::ParserExecutor<7> const* user_int_parserexec_data = d_user_int_attrib_parserexec.dataPtr(); - amrex::ParserExecutor<7> const* user_real_parserexec_data = d_user_real_attrib_parserexec.dataPtr(); -#else - int** pa_user_int_data = pa_user_int_pinned.dataPtr(); - ParticleReal** pa_user_real_data = pa_user_real_pinned.dataPtr(); - amrex::ParserExecutor<7> const* user_int_parserexec_data = user_int_attrib_parserexec_pinned.dataPtr(); - amrex::ParserExecutor<7> const* user_real_parserexec_data = user_real_attrib_parserexec_pinned.dataPtr(); -#endif + PlasmaParserHelper plasma_parser_helper (soa, old_size, m_user_int_attribs, m_user_real_attribs, particle_icomps, particle_comps, plasma_parser_wrapper); + int** pa_user_int_data = plasma_parser_helper.getUserIntDataPtrs(); + ParticleReal** pa_user_real_data = plasma_parser_helper.getUserRealDataPtrs(); + amrex::ParserExecutor<7> const* user_int_parserexec_data = plasma_parser_helper.getUserIntParserExecData(); + amrex::ParserExecutor<7> const* user_real_parserexec_data = plasma_parser_helper.getUserRealParserExecData(); int* p_ion_level = nullptr; if (do_field_ionization) { @@ -1741,34 +1513,9 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, } #ifdef WARPX_QED - //Pointer to the optical depth component - amrex::ParticleReal* p_optical_depth_QSR = nullptr; - amrex::ParticleReal* p_optical_depth_BW = nullptr; - - // If a QED effect is enabled, the corresponding optical depth - // has to be initialized - const bool loc_has_quantum_sync = has_quantum_sync(); - const bool loc_has_breit_wheeler = has_breit_wheeler(); - if (loc_has_quantum_sync) { - p_optical_depth_QSR = soa.GetRealData( - particle_comps["opticalDepthQSR"]).data() + old_size; - } - if(loc_has_breit_wheeler) { - p_optical_depth_BW = soa.GetRealData( - particle_comps["opticalDepthBW"]).data() + old_size; - } - - //If needed, get the appropriate functors from the engines - QuantumSynchrotronGetOpticalDepth quantum_sync_get_opt; - BreitWheelerGetOpticalDepth breit_wheeler_get_opt; - if(loc_has_quantum_sync){ - quantum_sync_get_opt = - m_shr_p_qs_engine->build_optical_depth_functor(); - } - if(loc_has_breit_wheeler){ - breit_wheeler_get_opt = - m_shr_p_bw_engine->build_optical_depth_functor(); - } + const QEDHelper qed_helper(soa, old_size, particle_comps, + has_quantum_sync(), has_breit_wheeler(), + m_shr_p_qs_engine, m_shr_p_bw_engine); #endif const bool loc_do_field_ionization = do_field_ionization; @@ -1787,6 +1534,24 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, { const IntVect iv = IntVect(AMREX_D_DECL(i, j, k)); const auto index = overlap_box.index(iv); + + Real scale_fac = compute_scale_fac_area(dx, num_ppc_real, flux_normal_axis); + + auto lo = getCellCoords(overlap_corner, dx, {0._rt, 0._rt, 0._rt}, iv); + auto hi = getCellCoords(overlap_corner, dx, {1._rt, 1._rt, 1._rt}, iv); + + if (flux_pos->overlapsWith(lo, hi)) + { + int r; + if (fine_overlap_box.ok() && fine_overlap_box.contains(iv)) { + r = compute_area_weights(lrrfac, flux_normal_axis); + } else { + r = 1; + } + scale_fac /= r; + } + amrex::ignore_unused(j,k); + for (int i_part = 0; i_part < pcounts[index]; ++i_part) { const long ip = poffset[index] + i_part; @@ -1878,14 +1643,15 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, } #ifdef WARPX_QED - if (loc_has_quantum_sync) { - p_optical_depth_QSR[ip] = quantum_sync_get_opt(engine); + if(qed_helper.has_quantum_sync){ + qed_helper.p_optical_depth_QSR[ip] = qed_helper.quantum_sync_get_opt(engine); } - if(loc_has_breit_wheeler){ - p_optical_depth_BW[ip] = breit_wheeler_get_opt(engine); + if(qed_helper.has_breit_wheeler){ + qed_helper.p_optical_depth_BW[ip] = qed_helper.breit_wheeler_get_opt(engine); } #endif + // Initialize user-defined integers with user-defined parser for (int ia = 0; ia < n_user_int_attribs; ++ia) { pa_user_int_data[ia][ip] = static_cast(user_int_parserexec_data[ia](pos.x, pos.y, pos.z, u.x, u.y, u.z, t)); @@ -1967,26 +1733,8 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, // are in the right tile.) tmp_pc.Redistribute(); - // Add the particles to the current container, tile by tile - for (int lev=0; levaddParticles(tmp_pc, true); } void @@ -1998,7 +1746,7 @@ PhysicalParticleContainer::Evolve (int lev, MultiFab* rho, MultiFab* crho, const MultiFab* cEx, const MultiFab* cEy, const MultiFab* cEz, const MultiFab* cBx, const MultiFab* cBy, const MultiFab* cBz, - Real /*t*/, Real dt, Real dt_next, DtType a_dt_type, bool skip_deposition, + Real /*t*/, Real dt, DtType a_dt_type, bool skip_deposition, PushType push_type) { @@ -2132,7 +1880,7 @@ PhysicalParticleContainer::Evolve (int lev, PushPX(pti, exfab, eyfab, ezfab, bxfab, byfab, bzfab, Ex.nGrowVect(), e_is_nodal, - 0, np_to_push, lev, gather_lev, dt, dt_next, ScaleFields(false), a_dt_type); + 0, np_to_push, lev, gather_lev, dt, ScaleFields(false), a_dt_type); } else if (push_type == PushType::Implicit) { ImplicitPushXP(pti, exfab, eyfab, ezfab, bxfab, byfab, bzfab, @@ -2174,7 +1922,7 @@ PhysicalParticleContainer::Evolve (int lev, cbxfab, cbyfab, cbzfab, cEx->nGrowVect(), e_is_nodal, nfine_gather, np-nfine_gather, - lev, lev-1, dt, dt_next, 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, @@ -2714,9 +2462,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, const long offset, const long np_to_push, int lev, int gather_lev, - amrex::Real dt, - amrex::Real dt_next, - ScaleFields scaleFields, + amrex::Real dt, ScaleFields scaleFields, DtType a_dt_type) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE((gather_lev==(lev-1)) || @@ -2893,7 +2639,6 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, copyAttribs(ip); } - // Push momentum from t^{n-1/2} to t^{n+1/2} doParticleMomentumPush<0>(ux[ip], uy[ip], uz[ip], Exp, Eyp, Ezp, Bxp, Byp, Bzp, ion_lev ? ion_lev[ip] : 1, @@ -2901,9 +2646,8 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, #ifdef WARPX_QED t_chi_max, #endif - 0.5_rt * (dt + dt_next)); + dt); - // Push position from t^n to t^{n+1} UpdatePosition(xp, yp, zp, ux[ip], uy[ip], uz[ip], dt); setPosition(ip, xp, yp, zp); } @@ -2922,7 +2666,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, t_chi_max, dt); - UpdatePosition(xp, yp, zp, ux[ip], uy[ip], uz[ip], 0.5_rt * (dt + dt_next)); + UpdatePosition(xp, yp, zp, ux[ip], uy[ip], uz[ip], dt); setPosition(ip, xp, yp, zp); } } @@ -3006,7 +2750,7 @@ PhysicalParticleContainer::ImplicitPushXP (WarpXParIter& pti, const Dim3 lo = lbound(box); - const int depos_type = WarpX::current_deposition_algo; + const auto depos_type = WarpX::current_deposition_algo; const int nox = WarpX::nox; const int n_rz_azimuthal_modes = WarpX::n_rz_azimuthal_modes; @@ -3416,6 +3160,28 @@ void PhysicalParticleContainer::resample (const int timestep, const bool verbose WARPX_PROFILE_VAR_STOP(blp_resample_actual); } +bool +PhysicalParticleContainer::findRefinedInjectionBox (amrex::Box& a_fine_injection_box, amrex::IntVect& a_rrfac) +{ + WARPX_PROFILE("PhysicalParticleContainer::findRefinedInjectionBox"); + + // This does not work if the mesh is dynamic. But in that case, we should + // not use refined injected either. We also assume there is only one fine level. + static bool refine_injection = false; + static Box fine_injection_box; + static amrex::IntVect rrfac(AMREX_D_DECL(1,1,1)); + if (!refine_injection and WarpX::moving_window_active(WarpX::GetInstance().getistep(0)+1) and WarpX::refine_plasma and do_continuous_injection and numLevels() == 2) { + refine_injection = true; + fine_injection_box = ParticleBoxArray(1).minimalBox(); + fine_injection_box.setSmall(WarpX::moving_window_dir, std::numeric_limits::lowest()/2); + fine_injection_box.setBig(WarpX::moving_window_dir, std::numeric_limits::max()/2); + rrfac = m_gdb->refRatio(0); + fine_injection_box.coarsen(rrfac); + } + a_fine_injection_box = fine_injection_box; + a_rrfac = rrfac; + return refine_injection; +} #ifdef WARPX_QED diff --git a/Source/Particles/RigidInjectedParticleContainer.H b/Source/Particles/RigidInjectedParticleContainer.H index e18931f1f01..bc20420ea6e 100644 --- a/Source/Particles/RigidInjectedParticleContainer.H +++ b/Source/Particles/RigidInjectedParticleContainer.H @@ -84,7 +84,6 @@ public: const amrex::MultiFab* cBz, amrex::Real t, amrex::Real dt, - amrex::Real dt_next, DtType a_dt_type=DtType::Full, bool skip_deposition=false, PushType push_type=PushType::Explicit) override; @@ -100,9 +99,7 @@ public: long offset, long np_to_push, int lev, int gather_lev, - amrex::Real dt, - amrex::Real dt_next, - ScaleFields scaleFields, + amrex::Real dt, ScaleFields scaleFields, DtType a_dt_type=DtType::Full) override; void PushP (int lev, amrex::Real dt, diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index c2a50173759..c3ec4c41131 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -166,8 +166,7 @@ RigidInjectedParticleContainer::PushPX (WarpXParIter& pti, const long offset, const long np_to_push, int lev, int gather_lev, - amrex::Real dt, amrex::Real /*dt_next*/, - ScaleFields /*scaleFields*/, + amrex::Real dt, ScaleFields /*scaleFields*/, DtType a_dt_type) { auto& attribs = pti.GetAttribs(); @@ -243,7 +242,7 @@ RigidInjectedParticleContainer::PushPX (WarpXParIter& pti, const bool do_scale = not done_injecting_lev; const Real v_boost = WarpX::beta_boost*PhysConst::c; PhysicalParticleContainer::PushPX(pti, exfab, eyfab, ezfab, bxfab, byfab, bzfab, - ngEB, e_is_nodal, offset, np_to_push, lev, gather_lev, dt, dt, + ngEB, e_is_nodal, offset, np_to_push, lev, gather_lev, dt, ScaleFields(do_scale, dt, zinject_plane_lev_previous, vzbeam_ave_boosted, v_boost), a_dt_type); @@ -300,7 +299,7 @@ RigidInjectedParticleContainer::Evolve (int lev, MultiFab* rho, MultiFab* crho, const MultiFab* cEx, const MultiFab* cEy, const MultiFab* cEz, const MultiFab* cBx, const MultiFab* cBy, const MultiFab* cBz, - Real t, Real dt, Real dt_next, DtType a_dt_type, bool skip_deposition, + Real t, Real dt, DtType a_dt_type, bool skip_deposition, PushType push_type) { @@ -326,7 +325,7 @@ RigidInjectedParticleContainer::Evolve (int lev, rho, crho, cEx, cEy, cEz, cBx, cBy, cBz, - t, dt, dt_next, a_dt_type, skip_deposition, push_type); + t, dt, a_dt_type, skip_deposition, push_type); } void diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 98fb3339e9c..7e882c151e8 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -153,7 +153,7 @@ public: amrex::MultiFab* rho, amrex::MultiFab* crho, const amrex::MultiFab* cEx, const amrex::MultiFab* cEy, const amrex::MultiFab* cEz, const amrex::MultiFab* cBx, const amrex::MultiFab* cBy, const amrex::MultiFab* cBz, - amrex::Real t, amrex::Real dt, amrex::Real dt_next, DtType a_dt_type=DtType::Full, bool skip_deposition=false, + amrex::Real t, amrex::Real dt, DtType a_dt_type=DtType::Full, bool skip_deposition=false, PushType push_type=PushType::Explicit) = 0; virtual void PostRestart () = 0; diff --git a/Source/WarpX.H b/Source/WarpX.H index 06b3537bfe8..823b652ef37 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -116,6 +116,7 @@ public: void InitData (); void Evolve (int numsteps = -1); + void Synchronize (); // // Functions used by implicit solvers @@ -1504,7 +1505,7 @@ private: amrex::Vector t_new; amrex::Vector t_old; amrex::Vector dt; - amrex::Vector dt_next; + int timestep_adaptation_interval = 0; // How often to update the timestep when using adaptive timestepping // Particle container std::unique_ptr mypc; From 55368d53c1b3376cb911a1d416538a49f3b14743 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Fri, 13 Sep 2024 14:44:25 -0700 Subject: [PATCH 29/32] finish refactor --- Docs/source/usage/parameters.rst | 2 +- .../Tests/electrostatic_sphere/CMakeLists.txt | 2 +- .../Tests/electrostatic_sphere/inputs_base_3d | 2 + .../inputs_base_constdt_3d | 6 --- .../inputs_test_3d_electrostatic_sphere | 2 +- ...puts_test_3d_electrostatic_sphere_adaptive | 53 ++++++++++++++++--- ...uts_test_3d_electrostatic_sphere_lab_frame | 2 +- ...electrostatic_sphere_lab_frame_mr_emass_10 | 2 +- ...uts_test_3d_electrostatic_sphere_rel_nodal | 2 +- Python/pywarpx/picmi.py | 18 ++++--- ...test_3d_electrostatic_sphere_adaptive.json | 20 +++---- Source/Evolve/WarpXComputeDt.cpp | 8 +-- Source/Evolve/WarpXEvolve.cpp | 5 +- Source/WarpX.H | 1 + Source/WarpX.cpp | 5 ++ 15 files changed, 89 insertions(+), 41 deletions(-) delete mode 100644 Examples/Tests/electrostatic_sphere/inputs_base_constdt_3d diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index 3eb9f81f9cd..b9d82d5014a 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -2139,7 +2139,7 @@ Time step This can be used with the electromagnetic solver, overriding ``warpx.cfl``, but it is up to the user to ensure that the CFL condition is met. -* ``warpx.timestep_adaptation_interval`` (`float`) optional (default `0`) +* ``warpx.dt_update_interval`` (`string`) optional (default `-1`) How many iterations pass between timestep adaptations when using the electrostatic solver. Must be greater than `0` to use adaptive timestepping, or else ``warpx.const_dt`` must be specified. diff --git a/Examples/Tests/electrostatic_sphere/CMakeLists.txt b/Examples/Tests/electrostatic_sphere/CMakeLists.txt index 2d1682b62bb..3d17c4462f8 100644 --- a/Examples/Tests/electrostatic_sphere/CMakeLists.txt +++ b/Examples/Tests/electrostatic_sphere/CMakeLists.txt @@ -47,7 +47,7 @@ add_warpx_test( 2 # nprocs inputs_test_3d_electrostatic_sphere_adaptive # inputs analysis_electrostatic_sphere.py # analysis - diags/diag1000018 # output + diags/diag1000054 # output OFF # dependency ) diff --git a/Examples/Tests/electrostatic_sphere/inputs_base_3d b/Examples/Tests/electrostatic_sphere/inputs_base_3d index 3bc302f6225..1f15b2da0de 100644 --- a/Examples/Tests/electrostatic_sphere/inputs_base_3d +++ b/Examples/Tests/electrostatic_sphere/inputs_base_3d @@ -1,3 +1,4 @@ +max_step = 30 amr.n_cell = 64 64 64 amr.max_level = 0 amr.blocking_factor = 8 @@ -7,6 +8,7 @@ geometry.prob_lo = -0.5 -0.5 -0.5 geometry.prob_hi = 0.5 0.5 0.5 boundary.field_lo = pec pec pec boundary.field_hi = pec pec pec +warpx.const_dt = 1e-6 warpx.do_electrostatic = relativistic particles.species_names = electron diff --git a/Examples/Tests/electrostatic_sphere/inputs_base_constdt_3d b/Examples/Tests/electrostatic_sphere/inputs_base_constdt_3d deleted file mode 100644 index da88ccb8602..00000000000 --- a/Examples/Tests/electrostatic_sphere/inputs_base_constdt_3d +++ /dev/null @@ -1,6 +0,0 @@ -# base inputs -FILE = inputs_base_3d - -# inputs for tests using constant dt -max_step = 30 -warpx.const_dt = 1e-6 diff --git a/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere b/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere index 84129349c7a..d89395e9d74 100644 --- a/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere +++ b/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere @@ -1,5 +1,5 @@ # base input parameters -FILE = inputs_base_constdt_3d +FILE = inputs_base_3d # test input parameters warpx.abort_on_warning_threshold = medium diff --git a/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_adaptive b/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_adaptive index 99245fc616f..f64f6de08ee 100644 --- a/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_adaptive +++ b/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_adaptive @@ -1,10 +1,47 @@ -# base input parameters -FILE = inputs_base_3d - -# adaptive timestepping params +stop_time = 60e-6 warpx.cfl = 0.2 -warpx.max_dt = 2.0e-6 -stop_time = 30e-6 +warpx.dt_update_interval = 10 +warpx.max_dt = 1.5e-6 +amr.n_cell = 64 64 64 +amr.max_level = 0 +amr.blocking_factor = 8 +amr.max_grid_size = 128 +geometry.dims = 3 +geometry.prob_lo = -0.5 -0.5 -0.5 +geometry.prob_hi = 0.5 0.5 0.5 +boundary.field_lo = pec pec pec +boundary.field_hi = pec pec pec +warpx.do_electrostatic = relativistic + +particles.species_names = electron + +algo.field_gathering = momentum-conserving + +# Order of particle shape factors +algo.particle_shape = 1 + +my_constants.n0 = 1.49e6 +my_constants.R0 = 0.1 + +electron.charge = -q_e +electron.mass = m_e +electron.injection_style = "NUniformPerCell" +electron.num_particles_per_cell_each_dim = 2 2 2 +electron.profile = parse_density_function +electron.density_function(x,y,z) = "(x*x + y*y + z*z < R0*R0)*n0" +electron.momentum_distribution_type = at_rest + +diagnostics.diags_names = diag1 diag2 + +diag1.intervals = 30 +diag1.diag_type = Full +diag1.fields_to_plot = Ex Ey Ez rho + +diag2.intervals = 30 +diag2.diag_type = Full +diag2.fields_to_plot = none +diag2.format = openpmd -# test input parameters -warpx.abort_on_warning_threshold = medium +warpx.reduced_diags_names = timestep +timestep.intervals = 1 +timestep.type = Timestep diff --git a/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_lab_frame b/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_lab_frame index 0a44b6a21b0..da97ae8afe7 100644 --- a/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_lab_frame +++ b/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_lab_frame @@ -1,5 +1,5 @@ # base input parameters -FILE = inputs_base_constdt_3d +FILE = inputs_base_3d # test input parameters diag2.electron.variables = x y z ux uy uz w phi diff --git a/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_lab_frame_mr_emass_10 b/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_lab_frame_mr_emass_10 index 420f6eb1e2d..481cc65f030 100644 --- a/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_lab_frame_mr_emass_10 +++ b/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_lab_frame_mr_emass_10 @@ -1,5 +1,5 @@ # base input parameters -FILE = inputs_base_constdt_3d +FILE = inputs_base_3d # test input parameters amr.max_level = 1 diff --git a/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_rel_nodal b/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_rel_nodal index 9c32802f116..96bff8aa9c7 100644 --- a/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_rel_nodal +++ b/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_rel_nodal @@ -1,5 +1,5 @@ # base input parameters -FILE = inputs_base_constdt_3d +FILE = inputs_base_3d # test input parameters warpx.abort_on_warning_threshold = medium diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index c39d87638c0..478b4d5802e 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -1871,12 +1871,6 @@ class ElectrostaticSolver(picmistandard.PICMI_ElectrostaticSolver): Parameters ---------- - warpx_cfl: float, optional - Fraction of the CFL condition for particle velocity vs grid size, used to set the timestep when employing adaptive timestepping - - warpx_max_dt: float, optional - The maximum allowable timestep when using adaptive timestepping (only available for non-electromagnetic solvers) - warpx_relativistic: bool, default=False Whether to use the relativistic solver or lab frame solver @@ -1885,6 +1879,16 @@ class ElectrostaticSolver(picmistandard.PICMI_ElectrostaticSolver): warpx_self_fields_verbosity: integer, default=2 Level of verbosity for the lab frame solver + + warpx_dt_update_interval: string, optional (default = -1) + How frequently the timestep is updated. Adaptive timestepping is disabled when this is <= 0. + + warpx_cfl: float, optional + Fraction of the CFL condition for particle velocity vs grid size, used to set the timestep when `dt_update_interval > 0`. + + warpx_max_dt: float, optional + The maximum allowable timestep when `dt_update_interval > 0`. + """ def init(self, kw): @@ -1893,6 +1897,7 @@ def init(self, kw): self.self_fields_verbosity = kw.pop("warpx_self_fields_verbosity", None) self.magnetostatic = kw.pop("warpx_magnetostatic", False) self.cfl = kw.pop("warpx_cfl", None) + self.dt_update_interval = kw.pop("dt_update_interval", None) self.max_dt = kw.pop("warpx_max_dt", None) def solver_initialize_inputs(self): @@ -1903,6 +1908,7 @@ def solver_initialize_inputs(self): # set adaptive timestepping parameters pywarpx.warpx.cfl = self.cfl + pywarpx.warpx.dt_update_interval = self.dt_update_interval pywarpx.warpx.max_dt = self.max_dt if self.relativistic: diff --git a/Regression/Checksum/benchmarks_json/test_3d_electrostatic_sphere_adaptive.json b/Regression/Checksum/benchmarks_json/test_3d_electrostatic_sphere_adaptive.json index 2444322b12d..7d3609042e4 100644 --- a/Regression/Checksum/benchmarks_json/test_3d_electrostatic_sphere_adaptive.json +++ b/Regression/Checksum/benchmarks_json/test_3d_electrostatic_sphere_adaptive.json @@ -1,17 +1,17 @@ { "lev=0": { - "Ex": 6.506693947590733, - "Ey": 6.506693947080617, - "Ez": 6.506693947486778, - "rho": 2.6092567413231304e-10 + "Ex": 5.17744476240019, + "Ey": 5.177444761672448, + "Ez": 5.17744475997011, + "rho": 2.609256741323131e-10 }, "electron": { - "particle_momentum_x": 9.686609803282714e-24, - "particle_momentum_y": 9.686609801152154e-24, - "particle_momentum_z": 9.686609801429006e-24, - "particle_position_x": 512.8986106449738, - "particle_position_y": 512.8986108927056, - "particle_position_z": 512.898610747885, + "particle_momentum_x": 1.3215018862020975e-23, + "particle_momentum_y": 1.3215018871951551e-23, + "particle_momentum_z": 1.3215018890007144e-23, + "particle_position_x": 912.2309662923217, + "particle_position_y": 912.2309668939561, + "particle_position_z": 912.2309677284211, "particle_weight": 6212.501525878906 } } diff --git a/Source/Evolve/WarpXComputeDt.cpp b/Source/Evolve/WarpXComputeDt.cpp index 7687b6e42e1..a721b210892 100644 --- a/Source/Evolve/WarpXComputeDt.cpp +++ b/Source/Evolve/WarpXComputeDt.cpp @@ -50,8 +50,8 @@ WarpX::ComputeDt () WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_const_dt.has_value(), errorMsg.str()); } else if (electromagnetic_solver_id == ElectromagneticSolverAlgo::None) { std::stringstream errorMsg; - errorMsg << "warpx.const_dt must be specified with the electrostatic solver, or the warpx.timestep_adaptation_interval must be > 0."; - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_const_dt.has_value() || timestep_adaptation_interval > 0, errorMsg.str()); + errorMsg << "warpx.const_dt must be specified with the electrostatic solver, or the warpx.dt_update_interval must be > 0."; + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_const_dt.has_value() || dt_update_interval.isActivated(), errorMsg.str()); } // Determine the appropriate timestep as limited by the speed of light @@ -64,7 +64,7 @@ WarpX::ComputeDt () electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) { // Computation of dt for spectral algorithm // (determined by the minimum cell size in all directions) - if (m_max_dt.has_value() && electrostatic_solver_id == ElectrostaticSolverAlgo::None) { + if (m_max_dt.has_value() && electromagnetic_solver_id == ElectromagneticSolverAlgo::None) { deltat = m_max_dt.value(); } else { deltat = cfl / PhysConst::c * minDim(dx); @@ -123,7 +123,7 @@ WarpX::UpdateDtFromParticleSpeeds () // Restrict to be less than user-specified maximum timestep, if present if (m_max_dt.has_value()) { - deltat_new = std::min(deltat_new, cfl * dx_min / max_v); + deltat_new = std::min(deltat_new, m_max_dt.value()); } // Update dt diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 526eff6e26f..9acbe734405 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -119,7 +119,10 @@ WarpX::Evolve (int numsteps) // Update timestep for electrostatic solver if a constant dt is not provided // This first synchronizes the position and velocity before setting the new timestep if (electromagnetic_solver_id == ElectromagneticSolverAlgo::None && - !m_const_dt.has_value() && (step % timestep_adaptation_interval == 0)) { + !m_const_dt.has_value() && dt_update_interval.contains(step+1)) { + if (verbose) { + amrex::Print() << Utils::TextMsg::Info("updating timestep"); + } Synchronize(); UpdateDtFromParticleSpeeds(); } diff --git a/Source/WarpX.H b/Source/WarpX.H index e66c0373cd2..6115111320b 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -1517,6 +1517,7 @@ private: amrex::Vector t_new; amrex::Vector t_old; amrex::Vector dt; + static utils::parser::IntervalsParser dt_update_interval; int timestep_adaptation_interval = 0; // How often to update the timestep when using adaptive timestepping // Particle container diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 9653e2f6e9a..773ef9ecb30 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -184,6 +184,8 @@ bool WarpX::do_multi_J = false; int WarpX::do_multi_J_n_depositions; bool WarpX::safe_guard_cells = false; +utils::parser::IntervalsParser WarpX::dt_update_interval; + std::map WarpX::multifab_map; std::map WarpX::imultifab_map; @@ -813,6 +815,9 @@ WarpX::ReadParameters () // Read timestepping options utils::parser::queryWithParser(pp_warpx, "const_dt", m_const_dt); utils::parser::queryWithParser(pp_warpx, "max_dt", m_max_dt); + std::vector dt_interval_vec = {"-1"}; + pp_warpx.queryarr("dt_update_interval", dt_interval_vec); + dt_update_interval = utils::parser::IntervalsParser(dt_interval_vec); // Filter currently not working with FDTD solver in RZ geometry: turn OFF by default // (see https://github.com/ECP-WarpX/WarpX/issues/1943) From 4e16d5db73d0b21adaa6de808a4494066d122378 Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Fri, 13 Sep 2024 14:50:37 -0700 Subject: [PATCH 30/32] remove unused variable in WarpX.H --- Source/WarpX.H | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/Source/WarpX.H b/Source/WarpX.H index 6115111320b..199d83706c4 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -1517,8 +1517,7 @@ private: amrex::Vector t_new; amrex::Vector t_old; amrex::Vector dt; - static utils::parser::IntervalsParser dt_update_interval; - int timestep_adaptation_interval = 0; // How often to update the timestep when using adaptive timestepping + static utils::parser::IntervalsParser dt_update_interval; // How often to update the timestep when using adaptive timestepping // Particle container std::unique_ptr mypc; From f236cecfb58849ac6a88f686d0ec5ad44d0c2c0f Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Sat, 14 Sep 2024 07:25:12 -0700 Subject: [PATCH 31/32] Update checksum for new test --- ...test_3d_electrostatic_sphere_adaptive.json | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/Regression/Checksum/benchmarks_json/test_3d_electrostatic_sphere_adaptive.json b/Regression/Checksum/benchmarks_json/test_3d_electrostatic_sphere_adaptive.json index 7d3609042e4..561fbf86669 100644 --- a/Regression/Checksum/benchmarks_json/test_3d_electrostatic_sphere_adaptive.json +++ b/Regression/Checksum/benchmarks_json/test_3d_electrostatic_sphere_adaptive.json @@ -1,17 +1,17 @@ { "lev=0": { - "Ex": 5.17744476240019, - "Ey": 5.177444761672448, - "Ez": 5.17744475997011, - "rho": 2.609256741323131e-10 + "Ex": 5.177444767224255, + "Ey": 5.177444767224254, + "Ez": 5.177444767224256, + "rho": 2.6092568008333797e-10 }, "electron": { - "particle_momentum_x": 1.3215018862020975e-23, - "particle_momentum_y": 1.3215018871951551e-23, - "particle_momentum_z": 1.3215018890007144e-23, - "particle_position_x": 912.2309662923217, - "particle_position_y": 912.2309668939561, - "particle_position_z": 912.2309677284211, + "particle_momentum_x": 1.3215019655285216e-23, + "particle_momentum_y": 1.3215019655285214e-23, + "particle_momentum_z": 1.3215019655285217e-23, + "particle_position_x": 912.2310003741203, + "particle_position_y": 912.2310003741203, + "particle_position_z": 912.2310003741202, "particle_weight": 6212.501525878906 } } From bea09410fedf6de6e2fd9cd9880ca61649f3011f Mon Sep 17 00:00:00 2001 From: Thomas Marks Date: Tue, 17 Sep 2024 11:18:00 -0400 Subject: [PATCH 32/32] Address comments --- Source/Evolve/WarpXComputeDt.cpp | 25 +++++++++++++------------ Source/WarpX.H | 4 ++++ 2 files changed, 17 insertions(+), 12 deletions(-) diff --git a/Source/Evolve/WarpXComputeDt.cpp b/Source/Evolve/WarpXComputeDt.cpp index a721b210892..b82cb6aff26 100644 --- a/Source/Evolve/WarpXComputeDt.cpp +++ b/Source/Evolve/WarpXComputeDt.cpp @@ -45,13 +45,12 @@ WarpX::ComputeDt () // Handle cases where the timestep is not limited by the speed of light // and no constant timestep is provided if (electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC) { - std::stringstream errorMsg; - errorMsg << "warpx.const_dt must be specified with the hybrid-PIC solver."; - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_const_dt.has_value(), errorMsg.str()); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_const_dt.has_value(), "warpx.const_dt must be specified with the hybrid-PIC solver."); } else if (electromagnetic_solver_id == ElectromagneticSolverAlgo::None) { - std::stringstream errorMsg; - errorMsg << "warpx.const_dt must be specified with the electrostatic solver, or the warpx.dt_update_interval must be > 0."; - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_const_dt.has_value() || dt_update_interval.isActivated(), errorMsg.str()); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + m_const_dt.has_value() || dt_update_interval.isActivated(), + "warpx.const_dt must be specified with the electrostatic solver, or warpx.dt_update_interval must be > 0." + ); } // Determine the appropriate timestep as limited by the speed of light @@ -60,15 +59,17 @@ WarpX::ComputeDt () if (m_const_dt.has_value()) { deltat = m_const_dt.value(); - } else if (electrostatic_solver_id != ElectrostaticSolverAlgo::None || - electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) { - // Computation of dt for spectral algorithm - // (determined by the minimum cell size in all directions) - if (m_max_dt.has_value() && electromagnetic_solver_id == ElectromagneticSolverAlgo::None) { + } else if (electrostatic_solver_id != ElectrostaticSolverAlgo::None) { + // Set dt for electrostatic algorithm + if (m_max_dt.has_value()) { deltat = m_max_dt.value(); } else { - deltat = cfl / PhysConst::c * minDim(dx); + deltat = cfl * minDim(dx) / PhysConst::c; } + } else if (electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) { + // Computation of dt for spectral algorithm + // (determined by the minimum cell size in all directions) + deltat = cfl * minDim(dx) / PhysConst::c; } else { // Computation of dt for FDTD algorithm #ifdef WARPX_DIM_RZ diff --git a/Source/WarpX.H b/Source/WarpX.H index 4632d83571a..5065fa73ff9 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -116,6 +116,10 @@ public: void InitData (); void Evolve (int numsteps = -1); + + /** Push momentum one half step forward to synchronize with position. + * Also sets is_synchronized to `true`. + */ void Synchronize (); //