diff options
Diffstat (limited to 'Source')
| -rw-r--r-- | Source/Evolve/WarpXEvolve.cpp | 158 | ||||
| -rw-r--r-- | Source/Particles/LaserParticleContainer.H | 3 | ||||
| -rw-r--r-- | Source/Particles/LaserParticleContainer.cpp | 8 | ||||
| -rw-r--r-- | Source/Particles/MultiParticleContainer.H | 2 | ||||
| -rw-r--r-- | Source/Particles/MultiParticleContainer.cpp | 20 | ||||
| -rw-r--r-- | Source/Particles/PhotonParticleContainer.H | 3 | ||||
| -rw-r--r-- | Source/Particles/PhotonParticleContainer.cpp | 4 | ||||
| -rw-r--r-- | Source/Particles/PhysicalParticleContainer.H | 3 | ||||
| -rw-r--r-- | Source/Particles/PhysicalParticleContainer.cpp | 10 | ||||
| -rw-r--r-- | Source/Particles/RigidInjectedParticleContainer.H | 3 | ||||
| -rw-r--r-- | Source/Particles/RigidInjectedParticleContainer.cpp | 4 | ||||
| -rw-r--r-- | Source/Particles/WarpXParticleContainer.H | 2 | ||||
| -rw-r--r-- | Source/WarpX.H | 4 |
13 files changed, 116 insertions, 108 deletions
diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 5efd000c6..3d89c15a0 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -135,7 +135,13 @@ WarpX::Evolve (int numsteps) // Main PIC operation: // gather fields, push particles, deposit sources, update fields - if (do_subcycling == 0 || finest_level == 0) { + if ( do_electrostatic != ElectrostaticSolverAlgo::None ) { + // Special case: electrostatic solver. + // In this case, we only gather fields and push particles + // The deposition and calculation of fields is done further below + bool const skip_deposition=true; + PushParticlesandDepose(cur_time, skip_deposition); + } else if (do_subcycling == 0 || finest_level == 0) { OneStep_nosub(cur_time); // E : guard cells are up-to-date // B : guard cells are NOT up-to-date @@ -327,75 +333,70 @@ WarpX::OneStep_nosub (Real cur_time) if (warpx_py_beforeEsolve) warpx_py_beforeEsolve(); - if( do_electrostatic == ElectrostaticSolverAlgo::None ) { - // Electromagnetic solver: - // Push E and B from {n} to {n+1} - // (And update guard cells immediately afterwards) - if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) { - if (use_hybrid_QED) - { - WarpX::Hybrid_QED_Push(dt); - FillBoundaryE(guard_cells.ng_alloc_EB); - } - PushPSATD(dt[0]); + // Push E and B from {n} to {n+1} + // (And update guard cells immediately afterwards) + if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) { + if (use_hybrid_QED) + { + WarpX::Hybrid_QED_Push(dt); FillBoundaryE(guard_cells.ng_alloc_EB); - FillBoundaryB(guard_cells.ng_alloc_EB); + } + PushPSATD(dt[0]); + FillBoundaryE(guard_cells.ng_alloc_EB); + FillBoundaryB(guard_cells.ng_alloc_EB); - if (use_hybrid_QED) - { - WarpX::Hybrid_QED_Push(dt); - FillBoundaryE(guard_cells.ng_alloc_EB); - } + if (use_hybrid_QED) { + WarpX::Hybrid_QED_Push(dt); + FillBoundaryE(guard_cells.ng_alloc_EB); + } - // Synchronize E and B fields on nodal points - NodalSyncE(); - NodalSyncB(); + // Synchronize E and B fields on nodal points + NodalSyncE(); + NodalSyncB(); - if (do_pml) { - DampPML(); - NodalSyncPML(); - } + if (do_pml) { + DampPML(); + NodalSyncPML(); + } + } else { + EvolveF(0.5_rt * dt[0], DtType::FirstHalf); + FillBoundaryF(guard_cells.ng_FieldSolverF); + EvolveB(0.5_rt * dt[0]); // We now have B^{n+1/2} + + if (do_silver_mueller) ApplySilverMuellerBoundary( dt[0] ); + FillBoundaryB(guard_cells.ng_FieldSolver); + + if (WarpX::em_solver_medium == MediumForEM::Vacuum) { + // vacuum medium + EvolveE(dt[0]); // We now have E^{n+1} + } else if (WarpX::em_solver_medium == MediumForEM::Macroscopic) { + // macroscopic medium + MacroscopicEvolveE(dt[0]); // We now have E^{n+1} } else { - EvolveF(0.5_rt * dt[0], DtType::FirstHalf); - FillBoundaryF(guard_cells.ng_FieldSolverF); - EvolveB(0.5_rt * dt[0]); // We now have B^{n+1/2} - - if (do_silver_mueller) ApplySilverMuellerBoundary( dt[0] ); - FillBoundaryB(guard_cells.ng_FieldSolver); - - if (WarpX::em_solver_medium == MediumForEM::Vacuum) { - // vacuum medium - EvolveE(dt[0]); // We now have E^{n+1} - } else if (WarpX::em_solver_medium == MediumForEM::Macroscopic) { - // macroscopic medium - MacroscopicEvolveE(dt[0]); // We now have E^{n+1} - } else { - amrex::Abort(" Medium for EM is unknown \n"); - } - - FillBoundaryE(guard_cells.ng_FieldSolver); - EvolveF(0.5_rt * dt[0], DtType::SecondHalf); - EvolveB(0.5_rt * dt[0]); // We now have B^{n+1} - - // Synchronize E and B fields on nodal points - NodalSyncE(); - NodalSyncB(); - - if (do_pml) { - FillBoundaryF(guard_cells.ng_alloc_F); - DampPML(); - NodalSyncPML(); - FillBoundaryE(guard_cells.ng_MovingWindow); - FillBoundaryF(guard_cells.ng_MovingWindow); - FillBoundaryB(guard_cells.ng_MovingWindow); - } - // E and B are up-to-date in the domain, but all guard cells are - // outdated. - if (safe_guard_cells) - FillBoundaryB(guard_cells.ng_alloc_EB); - } // !PSATD + amrex::Abort(" Medium for EM is unknown \n"); + } - } // !do_electrostatic + FillBoundaryE(guard_cells.ng_FieldSolver); + EvolveF(0.5_rt * dt[0], DtType::SecondHalf); + EvolveB(0.5_rt * dt[0]); // We now have B^{n+1} + + // Synchronize E and B fields on nodal points + NodalSyncE(); + NodalSyncB(); + + if (do_pml) { + FillBoundaryF(guard_cells.ng_alloc_F); + DampPML(); + NodalSyncPML(); + FillBoundaryE(guard_cells.ng_MovingWindow); + FillBoundaryF(guard_cells.ng_MovingWindow); + FillBoundaryB(guard_cells.ng_MovingWindow); + } + // E and B are up-to-date in the domain, but all guard cells are + // outdated. + if (safe_guard_cells) + FillBoundaryB(guard_cells.ng_alloc_EB); + } // !PSATD if (warpx_py_afterEsolve) warpx_py_afterEsolve(); } @@ -598,17 +599,17 @@ WarpX::doQEDEvents (int lev) #endif void -WarpX::PushParticlesandDepose (amrex::Real cur_time) +WarpX::PushParticlesandDepose (amrex::Real cur_time, bool skip_deposition) { // Evolve particles to p^{n+1/2} and x^{n+1} // Depose current, j^{n+1/2} for (int lev = 0; lev <= finest_level; ++lev) { - PushParticlesandDepose(lev, cur_time); + PushParticlesandDepose(lev, cur_time, DtType::Full, skip_deposition); } } void -WarpX::PushParticlesandDepose (int lev, amrex::Real cur_time, DtType a_dt_type) +WarpX::PushParticlesandDepose (int lev, amrex::Real cur_time, DtType a_dt_type, bool skip_deposition) { // If warpx.do_current_centering = 1, the current is deposited on the nodal MultiFab current_fp_nodal // and then centered onto the staggered MultiFab current_fp @@ -627,18 +628,19 @@ WarpX::PushParticlesandDepose (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); - + cur_time, dt[lev], a_dt_type, skip_deposition); #ifdef WARPX_DIM_RZ - // This is called after all particles have deposited their current and charge. - ApplyInverseVolumeScalingToCurrentDensity(current_fp[lev][0].get(), current_fp[lev][1].get(), current_fp[lev][2].get(), lev); - if (current_buf[lev][0].get()) { - ApplyInverseVolumeScalingToCurrentDensity(current_buf[lev][0].get(), current_buf[lev][1].get(), current_buf[lev][2].get(), lev-1); - } - if (rho_fp[lev].get()) { - ApplyInverseVolumeScalingToChargeDensity(rho_fp[lev].get(), lev); - if (charge_buf[lev].get()) { - ApplyInverseVolumeScalingToChargeDensity(charge_buf[lev].get(), lev-1); + if (! skip_deposition) { + // This is called after all particles have deposited their current and charge. + ApplyInverseVolumeScalingToCurrentDensity(current_fp[lev][0].get(), current_fp[lev][1].get(), current_fp[lev][2].get(), lev); + if (current_buf[lev][0].get()) { + ApplyInverseVolumeScalingToCurrentDensity(current_buf[lev][0].get(), current_buf[lev][1].get(), current_buf[lev][2].get(), lev-1); + } + if (rho_fp[lev].get()) { + ApplyInverseVolumeScalingToChargeDensity(rho_fp[lev].get(), lev); + if (charge_buf[lev].get()) { + ApplyInverseVolumeScalingToChargeDensity(charge_buf[lev].get(), lev-1); + } } } #endif diff --git a/Source/Particles/LaserParticleContainer.H b/Source/Particles/LaserParticleContainer.H index 98273b58a..eabd64bbd 100644 --- a/Source/Particles/LaserParticleContainer.H +++ b/Source/Particles/LaserParticleContainer.H @@ -46,7 +46,8 @@ 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) final; + amrex::Real t, amrex::Real dt, DtType a_dt_type=DtType::Full, + bool skip_deposition=false) final; virtual void PushP (int lev, amrex::Real dt, const amrex::MultiFab& , diff --git a/Source/Particles/LaserParticleContainer.cpp b/Source/Particles/LaserParticleContainer.cpp index 03b623a6d..c5fe3e4fb 100644 --- a/Source/Particles/LaserParticleContainer.cpp +++ b/Source/Particles/LaserParticleContainer.cpp @@ -418,7 +418,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, DtType /*a_dt_type*/) + Real t, Real dt, DtType /*a_dt_type*/, bool skip_deposition) { WARPX_PROFILE("LaserParticleContainer::Evolve()"); WARPX_PROFILE_VAR_NS("LaserParticleContainer::Evolve::ParticlePush", blp_pp); @@ -475,7 +475,7 @@ LaserParticleContainer::Evolve (int lev, plane_Yp.resize(np); amplitude_E.resize(np); - if (rho) { + if (rho && ! skip_deposition) { int* AMREX_RESTRICT ion_lev = nullptr; DepositCharge(pti, wp, ion_lev, rho, 0, 0, np_current, thread_num, lev, lev); @@ -510,7 +510,7 @@ LaserParticleContainer::Evolve (int lev, // Current Deposition // // Deposit inside domains - { + if (! skip_deposition ) { int* ion_lev = nullptr; DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, &jx, &jy, &jz, 0, np_current, thread_num, @@ -525,7 +525,7 @@ LaserParticleContainer::Evolve (int lev, } } - if (rho) { + if (rho && ! skip_deposition) { int* AMREX_RESTRICT ion_lev = nullptr; DepositCharge(pti, wp, ion_lev, rho, 1, 0, np_current, thread_num, lev, lev); diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 58bfd0e79..d6a7967b5 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -100,7 +100,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); + amrex::Real t, amrex::Real dt, DtType a_dt_type=DtType::Full, bool skip_deposition=false); /// /// This pushes the particle positions by one half time step for all the species in the diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 4218e9fdf..34d12a2b3 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -328,16 +328,18 @@ 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) + Real t, Real dt, DtType a_dt_type, bool skip_deposition) { - jx.setVal(0.0); - jy.setVal(0.0); - jz.setVal(0.0); - if (cjx) cjx->setVal(0.0); - if (cjy) cjy->setVal(0.0); - if (cjz) cjz->setVal(0.0); - if (rho) rho->setVal(0.0); - if (crho) crho->setVal(0.0); + if (! skip_deposition) { + jx.setVal(0.0); + jy.setVal(0.0); + jz.setVal(0.0); + if (cjx) cjx->setVal(0.0); + if (cjy) cjy->setVal(0.0); + if (cjz) cjz->setVal(0.0); + if (rho) rho->setVal(0.0); + if (crho) crho->setVal(0.0); + } 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); diff --git a/Source/Particles/PhotonParticleContainer.H b/Source/Particles/PhotonParticleContainer.H index 4d5ec590a..0cfb37729 100644 --- a/Source/Particles/PhotonParticleContainer.H +++ b/Source/Particles/PhotonParticleContainer.H @@ -52,7 +52,8 @@ public: const amrex::MultiFab* cBz, amrex::Real t, amrex::Real dt, - DtType a_dt_type=DtType::Full) override; + DtType a_dt_type=DtType::Full, + bool skip_deposition=false) override; virtual void PushPX(WarpXParIter& pti, amrex::FArrayBox const * exfab, diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index e8ba363a3..bda0f984e 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -198,7 +198,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*/) + Real t, Real dt, DtType a_dt_type, bool skip_deposition) { // This does gather, push and depose. // Push and depose have been re-written for photons @@ -210,6 +210,6 @@ PhotonParticleContainer::Evolve (int lev, rho, crho, cEx, cEy, cEz, cBx, cBy, cBz, - t, dt); + t, dt, a_dt_type, skip_deposition); } diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 250f3fb02..000caf790 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -120,7 +120,8 @@ public: const amrex::MultiFab* cBz, amrex::Real t, amrex::Real dt, - DtType a_dt_type=DtType::Full) override; + DtType a_dt_type=DtType::Full, + bool skip_deposition=false ) override; virtual void PushPX (WarpXParIter& pti, amrex::FArrayBox const * exfab, diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 6fbe442e3..681559fbc 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -953,7 +953,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) + Real /*t*/, Real dt, DtType a_dt_type, bool skip_deposition) { WARPX_PROFILE("PhysicalParticleContainer::Evolve()"); @@ -1055,7 +1055,7 @@ PhysicalParticleContainer::Evolve (int lev, const long np_current = (cjx) ? nfine_current : np; - if (rho) { + if (rho && ! skip_deposition) { // Deposit charge before particle push, in component 0 of MultiFab rho. int* AMREX_RESTRICT ion_lev; if (do_field_ionization){ @@ -1125,9 +1125,9 @@ PhysicalParticleContainer::Evolve (int lev, WARPX_PROFILE_VAR_STOP(blp_fg); // - // Current Deposition (only needed for electromagnetic solver) + // Current Deposition // - if (WarpX::do_electrostatic == ElectrostaticSolverAlgo::None) { + if (! skip_deposition) { int* AMREX_RESTRICT ion_lev; if (do_field_ionization){ ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); @@ -1147,7 +1147,7 @@ PhysicalParticleContainer::Evolve (int lev, } // end of "if do_electrostatic == ElectrostaticSolverAlgo::None" } // end of "if do_not_push" - if (rho) { + if (rho && ! skip_deposition) { // Deposit charge after particle push, in component 1 of MultiFab rho. // (Skipped for electrostatic solver, as this may lead to out-of-bounds) if (WarpX::do_electrostatic == ElectrostaticSolverAlgo::None) { diff --git a/Source/Particles/RigidInjectedParticleContainer.H b/Source/Particles/RigidInjectedParticleContainer.H index 70ac7677e..1f52c0861 100644 --- a/Source/Particles/RigidInjectedParticleContainer.H +++ b/Source/Particles/RigidInjectedParticleContainer.H @@ -66,7 +66,8 @@ public: const amrex::MultiFab* cBz, amrex::Real t, amrex::Real dt, - DtType a_dt_type=DtType::Full) override; + DtType a_dt_type=DtType::Full, + bool skip_deposition=false ) override; virtual void PushPX (WarpXParIter& pti, amrex::FArrayBox const * exfab, diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index dc2f88239..5cb49f69c 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -347,7 +347,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) + Real t, Real dt, DtType a_dt_type, bool skip_deposition) { // Update location of injection plane in the boosted frame @@ -374,7 +374,7 @@ RigidInjectedParticleContainer::Evolve (int lev, rho, crho, cEx, cEy, cEz, cBx, cBy, cBz, - t, dt, a_dt_type); + t, dt, a_dt_type, skip_deposition); } void diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 91d5e64ee..07a9ddcd6 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -164,7 +164,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) = 0; + amrex::Real t, amrex::Real dt, DtType a_dt_type=DtType::Full, bool skip_deposition=false) = 0; virtual void PostRestart () = 0; diff --git a/Source/WarpX.H b/Source/WarpX.H index 73c253d2d..3973ed9c6 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -443,8 +443,8 @@ public: void doQEDEvents (int lev); #endif - void PushParticlesandDepose (int lev, amrex::Real cur_time, DtType a_dt_type=DtType::Full); - void PushParticlesandDepose ( amrex::Real cur_time); + void PushParticlesandDepose (int lev, amrex::Real cur_time, DtType a_dt_type=DtType::Full, bool skip_current=false); + void PushParticlesandDepose ( amrex::Real cur_time, bool skip_current=false); // This function does aux(lev) = fp(lev) + I(aux(lev-1)-cp(lev)). // Caller must make sure fp and cp have ghost cells filled. |