From d1018ef5b4b263ba5859b334af18223130dd9788 Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Mon, 6 Aug 2018 15:17:27 -0700 Subject: [PATCH 01/43] extend the boosted frame diagnostics to multiple levels of refinement. --- Source/ParticleContainer.cpp | 61 +++++++++++---------- Source/PhysicalParticleContainer.cpp | 18 ++++--- Source/WarpX.cpp | 60 ++++++++++----------- Source/WarpXBoostedFrameDiagnostic.cpp | 1 - Source/WarpXIO.cpp | 74 ++++++++++++++------------ Source/WarpXParticleContainer.H | 2 +- 6 files changed, 114 insertions(+), 102 deletions(-) diff --git a/Source/ParticleContainer.cpp b/Source/ParticleContainer.cpp index e0539e799a3..d87b41132d9 100644 --- a/Source/ParticleContainer.cpp +++ b/Source/ParticleContainer.cpp @@ -297,36 +297,39 @@ ::GetLabFrameData(const std::string& snapshot_name, WarpXParticleContainer* pc = allcontainers[i].get(); WarpXParticleContainer::DiagnosticParticles diagnostic_particles; pc->GetParticleSlice(direction, z_old, z_new, t_boost, t_lab, dt, diagnostic_particles); - - for (auto it = diagnostic_particles.begin(); it != diagnostic_particles.end(); ++it) - { - parts[i].GetRealData(DiagIdx::w).insert( parts[i].GetRealData(DiagIdx::w ).end(), - it->second.GetRealData(DiagIdx::w ).begin(), - it->second.GetRealData(DiagIdx::w ).end()); - - parts[i].GetRealData(DiagIdx::x).insert( parts[i].GetRealData(DiagIdx::x ).end(), - it->second.GetRealData(DiagIdx::x ).begin(), - it->second.GetRealData(DiagIdx::x ).end()); - - parts[i].GetRealData(DiagIdx::y).insert( parts[i].GetRealData(DiagIdx::y ).end(), - it->second.GetRealData(DiagIdx::y ).begin(), - it->second.GetRealData(DiagIdx::y ).end()); - - parts[i].GetRealData(DiagIdx::z).insert( parts[i].GetRealData(DiagIdx::z ).end(), - it->second.GetRealData(DiagIdx::z ).begin(), - it->second.GetRealData(DiagIdx::z ).end()); - parts[i].GetRealData(DiagIdx::ux).insert( parts[i].GetRealData(DiagIdx::ux).end(), - it->second.GetRealData(DiagIdx::ux).begin(), - it->second.GetRealData(DiagIdx::ux).end()); - - parts[i].GetRealData(DiagIdx::uy).insert( parts[i].GetRealData(DiagIdx::uy).end(), - it->second.GetRealData(DiagIdx::uy).begin(), - it->second.GetRealData(DiagIdx::uy).end()); - - parts[i].GetRealData(DiagIdx::uz).insert( parts[i].GetRealData(DiagIdx::uz).end(), - it->second.GetRealData(DiagIdx::uz).begin(), - it->second.GetRealData(DiagIdx::uz).end()); + for (int lev = 0; lev <= pc->finestLevel(); ++lev) + { + for (auto it = diagnostic_particles[lev].begin(); it != diagnostic_particles[lev].end(); ++it) + { + parts[i].GetRealData(DiagIdx::w).insert( parts[i].GetRealData(DiagIdx::w ).end(), + it->second.GetRealData(DiagIdx::w ).begin(), + it->second.GetRealData(DiagIdx::w ).end()); + + parts[i].GetRealData(DiagIdx::x).insert( parts[i].GetRealData(DiagIdx::x ).end(), + it->second.GetRealData(DiagIdx::x ).begin(), + it->second.GetRealData(DiagIdx::x ).end()); + + parts[i].GetRealData(DiagIdx::y).insert( parts[i].GetRealData(DiagIdx::y ).end(), + it->second.GetRealData(DiagIdx::y ).begin(), + it->second.GetRealData(DiagIdx::y ).end()); + + parts[i].GetRealData(DiagIdx::z).insert( parts[i].GetRealData(DiagIdx::z ).end(), + it->second.GetRealData(DiagIdx::z ).begin(), + it->second.GetRealData(DiagIdx::z ).end()); + + parts[i].GetRealData(DiagIdx::ux).insert( parts[i].GetRealData(DiagIdx::ux).end(), + it->second.GetRealData(DiagIdx::ux).begin(), + it->second.GetRealData(DiagIdx::ux).end()); + + parts[i].GetRealData(DiagIdx::uy).insert( parts[i].GetRealData(DiagIdx::uy).end(), + it->second.GetRealData(DiagIdx::uy).begin(), + it->second.GetRealData(DiagIdx::uy).end()); + + parts[i].GetRealData(DiagIdx::uz).insert( parts[i].GetRealData(DiagIdx::uz).end(), + it->second.GetRealData(DiagIdx::uz).begin(), + it->second.GetRealData(DiagIdx::uz).end()); + } } } } diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp index 2e4ff73ba41..035e78b496e 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/PhysicalParticleContainer.cpp @@ -1196,6 +1196,8 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real slice_box.setLo(direction, z_min); slice_box.setHi(direction, z_max); + diagnostic_particles.resize(finestLevel()+1); + for (int lev = 0; lev < nlevs; ++lev) { const Real* dx = Geom(lev).CellSize(); @@ -1205,7 +1207,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { auto index = std::make_pair(pti.index(), pti.LocalTileIndex()); - diagnostic_particles[index]; + diagnostic_particles[lev][index]; } #ifdef _OPENMP @@ -1276,15 +1278,15 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real Real uyp = uyp_old[i]*weight_old + uyp_new[i]*weight_new; Real uzp = uz_old_p *weight_old + uz_new_p *weight_new; - diagnostic_particles[index].GetRealData(DiagIdx::w).push_back(wp[i]); + diagnostic_particles[lev][index].GetRealData(DiagIdx::w).push_back(wp[i]); - diagnostic_particles[index].GetRealData(DiagIdx::x).push_back(xp); - diagnostic_particles[index].GetRealData(DiagIdx::y).push_back(yp); - diagnostic_particles[index].GetRealData(DiagIdx::z).push_back(zp); + diagnostic_particles[lev][index].GetRealData(DiagIdx::x).push_back(xp); + diagnostic_particles[lev][index].GetRealData(DiagIdx::y).push_back(yp); + diagnostic_particles[lev][index].GetRealData(DiagIdx::z).push_back(zp); - diagnostic_particles[index].GetRealData(DiagIdx::ux).push_back(uxp); - diagnostic_particles[index].GetRealData(DiagIdx::uy).push_back(uyp); - diagnostic_particles[index].GetRealData(DiagIdx::uz).push_back(uzp); + diagnostic_particles[lev][index].GetRealData(DiagIdx::ux).push_back(uxp); + diagnostic_particles[lev][index].GetRealData(DiagIdx::uy).push_back(uyp); + diagnostic_particles[lev][index].GetRealData(DiagIdx::uz).push_back(uzp); } } } diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index d550212e54b..5367250a177 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -216,30 +216,30 @@ WarpX::ReadParameters () pp.query("verbose", verbose); pp.query("regrid_int", regrid_int); - // Boosted-frame parameters - pp.query("gamma_boost", gamma_boost); - beta_boost = std::sqrt(1.-1./pow(gamma_boost,2)); - if( gamma_boost > 1. ){ - // Read the boost direction - std::string s; + // Boosted-frame parameters + pp.query("gamma_boost", gamma_boost); + beta_boost = std::sqrt(1.-1./pow(gamma_boost,2)); + if( gamma_boost > 1. ){ + // Read the boost direction + std::string s; pp.get("boost_direction", s); if (s == "x" || s == "X") { - boost_direction[0] = 1.; - } + boost_direction[0] = 1.; + } #if (AMREX_SPACEDIM == 3) else if (s == "y" || s == "Y") { - boost_direction[1] = 1.; + boost_direction[1] = 1.; } #endif else if (s == "z" || s == "Z") { - boost_direction[2] = 1.; + boost_direction[2] = 1.; } - else { + else { const std::string msg = "Unknown boost_dir: "+s; amrex::Abort(msg.c_str()); } - } - + } + pp.queryarr("B_external", B_external); pp.query("do_moving_window", do_moving_window); @@ -262,33 +262,33 @@ WarpX::ReadParameters () const std::string msg = "Unknown moving_window_dir: "+s; amrex::Abort(msg.c_str()); } - + moving_window_x = geom[0].ProbLo(moving_window_dir); - + pp.get("moving_window_v", moving_window_v); moving_window_v *= PhysConst::c; } - + pp.query("do_plasma_injection", do_plasma_injection); if (do_plasma_injection) { - pp.get("num_injected_species", num_injected_species); - injected_plasma_species.resize(num_injected_species); - pp.getarr("injected_plasma_species", injected_plasma_species, - 0, num_injected_species); - if (moving_window_v >= 0){ - // Inject particles continuously from the right end of the box - current_injection_position = geom[0].ProbHi(moving_window_dir); - } else { - // Inject particles continuously from the left end of the box - current_injection_position = geom[0].ProbLo(moving_window_dir); - } + pp.get("num_injected_species", num_injected_species); + injected_plasma_species.resize(num_injected_species); + pp.getarr("injected_plasma_species", injected_plasma_species, + 0, num_injected_species); + if (moving_window_v >= 0){ + // Inject particles continuously from the right end of the box + current_injection_position = geom[0].ProbHi(moving_window_dir); + } else { + // Inject particles continuously from the left end of the box + current_injection_position = geom[0].ProbLo(moving_window_dir); + } } - + pp.query("do_boosted_frame_diagnostic", do_boosted_frame_diagnostic); if (do_boosted_frame_diagnostic) { - + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(gamma_boost > 1.0, - "gamma_boost must be > 1 to use the boosted frame diagnostic."); + "gamma_boost must be > 1 to use the boosted frame diagnostic."); std::string s; pp.get("boost_direction", s); diff --git a/Source/WarpXBoostedFrameDiagnostic.cpp b/Source/WarpXBoostedFrameDiagnostic.cpp index 9f04436bbde..7115f836de2 100644 --- a/Source/WarpXBoostedFrameDiagnostic.cpp +++ b/Source/WarpXBoostedFrameDiagnostic.cpp @@ -219,7 +219,6 @@ BoostedFrameDiagnostic:: writeParticleData(const WarpXParticleContainer::DiagnosticParticleData& pdata, const std::string& name, const int i_lab) { - BL_PROFILE("BoostedFrameDiagnostic::writeParticleData"); std::string field_name; diff --git a/Source/WarpXIO.cpp b/Source/WarpXIO.cpp index 2e6fd264d1a..05d3ccfeae5 100644 --- a/Source/WarpXIO.cpp +++ b/Source/WarpXIO.cpp @@ -401,47 +401,55 @@ WarpX::GetCellCenteredData() { const int ng = 1; const int nc = 10; - const int lev = 0; - auto cc = std::unique_ptr( new MultiFab(boxArray(lev), - DistributionMap(lev), - nc, ng) ); - Vector srcmf(AMREX_SPACEDIM); - int dcomp = 0; - - // first the electric field - PackPlotDataPtrs(srcmf, Efield_aux[lev]); - amrex::average_edge_to_cellcenter(*cc, dcomp, srcmf); + Vector > cc(finest_level+1); + + for (int lev = 0; lev <= finest_level; ++lev) + { + cc[lev].reset( new MultiFab(grids[lev], dmap[lev], nc, ng) ); + + Vector srcmf(AMREX_SPACEDIM); + int dcomp = 0; + + // first the electric field + PackPlotDataPtrs(srcmf, Efield_aux[lev]); + amrex::average_edge_to_cellcenter(*cc[lev], dcomp, srcmf); #if (AMREX_SPACEDIM == 2) - MultiFab::Copy(*cc, *cc, dcomp+1, dcomp+2, 1, ng); - amrex::average_node_to_cellcenter(*cc, dcomp+1, *Efield_aux[lev][1], 0, 1); + MultiFab::Copy(*cc[lev], *cc[lev], dcomp+1, dcomp+2, 1, ng); + amrex::average_node_to_cellcenter(*cc[lev], dcomp+1, *Efield_aux[lev][1], 0, 1); #endif - dcomp += 3; - - // then the magnetic field - PackPlotDataPtrs(srcmf, Bfield_aux[lev]); - amrex::average_face_to_cellcenter(*cc, dcomp, srcmf); + dcomp += 3; + + // then the magnetic field + PackPlotDataPtrs(srcmf, Bfield_aux[lev]); + amrex::average_face_to_cellcenter(*cc[lev], dcomp, srcmf); #if (AMREX_SPACEDIM == 2) - MultiFab::Copy(*cc, *cc, dcomp+1, dcomp+2, 1, ng); - MultiFab::Copy(*cc, *Bfield_aux[lev][1], 0, dcomp+1, 1, ng); + MultiFab::Copy(*cc[lev], *cc[lev], dcomp+1, dcomp+2, 1, ng); + MultiFab::Copy(*cc[lev], *Bfield_aux[lev][1], 0, dcomp+1, 1, ng); #endif - dcomp += 3; - - // then the current density - PackPlotDataPtrs(srcmf, current_fp[lev]); - amrex::average_edge_to_cellcenter(*cc, dcomp, srcmf); + dcomp += 3; + + // then the current density + PackPlotDataPtrs(srcmf, current_fp[lev]); + amrex::average_edge_to_cellcenter(*cc[lev], dcomp, srcmf); #if (AMREX_SPACEDIM == 2) - MultiFab::Copy(*cc, *cc, dcomp+1, dcomp+2, 1, ng); - amrex::average_node_to_cellcenter(*cc, dcomp+1, *current_fp[lev][1], 0, 1); + MultiFab::Copy(*cc[lev], *cc[lev], dcomp+1, dcomp+2, 1, ng); + amrex::average_node_to_cellcenter(*cc[lev], dcomp+1, *current_fp[lev][1], 0, 1); #endif - dcomp += 3; - - const std::unique_ptr& charge_density = mypc->GetChargeDensity(lev); - amrex::average_node_to_cellcenter(*cc, dcomp, *charge_density, 0, 1); - - cc->FillBoundary(geom[lev].periodicity()); + dcomp += 3; + + const std::unique_ptr& charge_density = mypc->GetChargeDensity(lev); + amrex::average_node_to_cellcenter(*cc[lev], dcomp, *charge_density, 0, 1); + + cc[lev]->FillBoundary(geom[lev].periodicity()); + } + + for (int lev = finest_level; lev > 0; --lev) + { + amrex::average_down(*cc[lev], *cc[lev-1], 0, nc, refRatio(lev-1)); + } - return cc; + return std::move(cc[0]); } void diff --git a/Source/WarpXParticleContainer.H b/Source/WarpXParticleContainer.H index f97c5cafed6..ad52b92b9d5 100644 --- a/Source/WarpXParticleContainer.H +++ b/Source/WarpXParticleContainer.H @@ -70,7 +70,7 @@ public: friend MultiParticleContainer; using DiagnosticParticleData = amrex::StructOfArrays; - using DiagnosticParticles = std::map, DiagnosticParticleData>; + using DiagnosticParticles = amrex::Vector, DiagnosticParticleData> >; WarpXParticleContainer (amrex::AmrCore* amr_core, int ispecies); virtual ~WarpXParticleContainer() {} From e9887911d0a97a55e11c19d6d876e0c810ecb275 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Mon, 10 Sep 2018 15:37:06 -0700 Subject: [PATCH 02/43] start reimplementing subcycling method 1 --- Source/WarpX.H | 11 +++++++++++ Source/WarpX.cpp | 14 ++++++++++++++ Source/WarpXEvolve.cpp | 10 ++++++++-- 3 files changed, 33 insertions(+), 2 deletions(-) diff --git a/Source/WarpX.H b/Source/WarpX.H index f8188dfe880..b7f68a569a9 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -34,6 +34,12 @@ enum struct DtType : int SecondHalf }; +enum struct PatchType : int +{ + fine, + coarse +}; + class WarpX : public amrex::AmrCore { @@ -365,6 +371,9 @@ private: amrex::Vector, 3 > > Efield_fp; amrex::Vector, 3 > > Bfield_fp; + // store fine patch + amrex::Vector, 3 > > current_store; + // Coarse patch amrex::Vector< std::unique_ptr > F_cp; amrex::Vector< std::unique_ptr > rho_cp; @@ -410,6 +419,8 @@ private: // Other runtime parameters int verbose = 1; + int do_subcycling = 0; + int max_step = std::numeric_limits::max(); amrex::Real stop_time = std::numeric_limits::max(); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 894ca509be2..e3e1b43d30e 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -150,6 +150,8 @@ WarpX::WarpX () Efield_fp.resize(nlevs_max); Bfield_fp.resize(nlevs_max); + current_store.resize(nlevs_max); + F_cp.resize(nlevs_max); rho_cp.resize(nlevs_max); current_cp.resize(nlevs_max); @@ -231,6 +233,10 @@ WarpX::ReadParameters () pp.query("cfl", cfl); pp.query("verbose", verbose); pp.query("regrid_int", regrid_int); + pp.query("do_subcycling", do_subcycling); + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(do_subcycling != 1 || max_level <= 1, + "Subcycling method 1 only works for 2 levels."); ReadBoostedFrameParameters(gamma_boost, beta_boost, boost_direction); @@ -444,6 +450,8 @@ WarpX::ClearLevel (int lev) Efield_fp [lev][i].reset(); Bfield_fp [lev][i].reset(); + current_store[lev][i].reset(); + current_cp[lev][i].reset(); Efield_cp [lev][i].reset(); Bfield_cp [lev][i].reset(); @@ -537,6 +545,12 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d current_fp[lev][1].reset( new MultiFab(amrex::convert(ba,jy_nodal_flag),dm,1,ngJ)); current_fp[lev][2].reset( new MultiFab(amrex::convert(ba,jz_nodal_flag),dm,1,ngJ)); + if (do_subcycling == 1 && lev == 1) { + current_store[lev][0].reset( new MultiFab(amrex::convert(ba,jx_nodal_flag),dm,1,ngJ)); + current_store[lev][1].reset( new MultiFab(amrex::convert(ba,jy_nodal_flag),dm,1,ngJ)); + current_store[lev][2].reset( new MultiFab(amrex::convert(ba,jz_nodal_flag),dm,1,ngJ)); + } + if (do_dive_cleaning) { F_fp[lev].reset (new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,1, ngF)); diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index 2906deb1d1c..bdf35c620ca 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -675,8 +675,8 @@ WarpX::ComputeDt () if (maxwell_fdtd_solver_id == 0) { // CFL time step Yee solver deltat = cfl * 1./( std::sqrt(AMREX_D_TERM( 1./(dx[0]*dx[0]), - + 1./(dx[1]*dx[1]), - + 1./(dx[2]*dx[2]))) * PhysConst::c ); + + 1./(dx[1]*dx[1]), + + 1./(dx[2]*dx[2]))) * PhysConst::c ); } else { // CFL time step CKC solver #if (BL_SPACEDIM == 3) @@ -689,6 +689,12 @@ WarpX::ComputeDt () dt.resize(0); dt.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]; + } + } + if (do_electrostatic) { dt[0] = const_dt; } From 2de5551d59fe3c2d0bfa10f7c49f1ff149486ea8 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Mon, 10 Sep 2018 17:25:38 -0700 Subject: [PATCH 03/43] WIP --- Source/WarpX.H | 10 +- Source/WarpXComm.cpp | 37 +++ Source/WarpXEvolve.cpp | 554 ++++++++++++++++++++++++----------------- Source/WarpXPML.H | 3 + Source/WarpXPML.cpp | 14 +- 5 files changed, 390 insertions(+), 228 deletions(-) diff --git a/Source/WarpX.H b/Source/WarpX.H index b7f68a569a9..d0fb8d9f7fc 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -235,6 +235,15 @@ private: /// void EvolveEM(int numsteps); + void OneStep_nosub (amrex::Real t); + void OneStep_sub1 (amrex::Real t); + + void RestrictCurrentFromFineToCoarsePatch (int lev); + void SumBoundaryJ (int lev, PatchType patch_type); + + void EvolveB (int lev, PatchType patch_type, amrex::Real dt); + void EvolveE (int lev, PatchType patch_type, amrex::Real dt); + #ifdef WARPX_DO_ELECTROSTATIC /// /// Advance the simulation by numsteps steps, electrostatic case. @@ -332,7 +341,6 @@ private: void LoadBalance (); - static void applyFilter (amrex::MultiFab& dstmf, const amrex::MultiFab& srcmf); void BuildBufferMasks (); diff --git a/Source/WarpXComm.cpp b/Source/WarpXComm.cpp index 156be97a2d0..b5297e04ae3 100644 --- a/Source/WarpXComm.cpp +++ b/Source/WarpXComm.cpp @@ -562,3 +562,40 @@ WarpX::SyncRho (const MultiFab& fine, MultiFab& crse, int ref_ratio) } } +void +WarpX::RestrictCurrentFromFineToCoarsePatch (int lev) +{ + current_cp[lev][0]->setVal(0.0); + current_cp[lev][1]->setVal(0.0); + current_cp[lev][2]->setVal(0.0); + + const IntVect& ref_ratio = refRatio(lev-1); + + std::array fine { current_fp[lev][0].get(), + current_fp[lev][1].get(), + current_fp[lev][2].get() }; + std::array< MultiFab*,3> crse { current_cp[lev][0].get(), + current_cp[lev][1].get(), + current_cp[lev][2].get() }; + SyncCurrent(fine, crse, ref_ratio[0]); +} + +void +WarpX::SumBoundaryJ (int lev, PatchType patch_type) +{ + if (patch_type == PatchType::fine) + { + const auto& period = Geom(lev).periodicity(); + current_fp[lev][0]->SumBoundary(period); + current_fp[lev][1]->SumBoundary(period); + current_fp[lev][2]->SumBoundary(period); + } + else if (patch_type == PatchType::coarse) + { + const auto& cperiod = Geom(lev-1).periodicity(); + current_cp[lev][0]->SumBoundary(cperiod); + current_cp[lev][1]->SumBoundary(cperiod); + current_cp[lev][2]->SumBoundary(cperiod); + } +} + diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index bdf35c620ca..30415cf5c14 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -82,7 +82,7 @@ WarpX::EvolveEM (int numsteps) UpdateAuxilaryData(); // on first step, push p by -0.5*dt for (int lev = 0; lev <= finest_level; ++lev) { - mypc->PushP(lev, -0.5*dt[0], + mypc->PushP(lev, -0.5*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]); } @@ -91,44 +91,22 @@ WarpX::EvolveEM (int numsteps) // Beyond one step, we have E^{n} and B^{n}. // Particles have p^{n-1/2} and x^{n}. UpdateAuxilaryData(); - } - - // Push particle from x^{n} to x^{n+1} - // from p^{n-1/2} to p^{n+1/2} - // Deposit current j^{n+1/2} - // Deposit charge density rho^{n} - PushParticlesandDepose(cur_time); - - SyncCurrent(); - - SyncRho(rho_fp, rho_cp); + } - // Push E and B from {n} to {n+1} - // (And update guard cells immediately afterwards) -#ifdef WARPX_USE_PSATD - PushPSATD(dt[0]); - FillBoundaryE(); - FillBoundaryB(); -#else - EvolveF(0.5*dt[0], DtType::FirstHalf); - EvolveB(0.5*dt[0]); // We now have B^{n+1/2} - FillBoundaryB(); - EvolveE(dt[0]); // We now have E^{n+1} - FillBoundaryE(); - EvolveF(0.5*dt[0], DtType::SecondHalf); - EvolveB(0.5*dt[0]); // We now have B^{n+1} - if (do_pml) { - DampPML(); - FillBoundaryE(); + if (do_subcycling == 0 || finest_level == 0) { + OneStep_nosub(cur_time); + } else if (do_subcycling == 1 && finest_level == 1) { + OneStep_sub1(cur_time); + } else { + amrex::Print() << "Error: do_subcycling = " << do_subcycling << std::endl; + amrex::Abort("Unsupported do_subcycling type"); } - FillBoundaryB(); -#endif 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 UpdateAuxilaryData(); for (int lev = 0; lev <= finest_level; ++lev) { - mypc->PushP(lev, 0.5*dt[0], + mypc->PushP(lev, 0.5*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]); } @@ -227,6 +205,141 @@ WarpX::EvolveEM (int numsteps) } } +void +WarpX::OneStep_nosub (Real cur_time) +{ + // Push particle from x^{n} to x^{n+1} + // from p^{n-1/2} to p^{n+1/2} + // Deposit current j^{n+1/2} + // Deposit charge density rho^{n} + PushParticlesandDepose(cur_time); + + SyncCurrent(); + + SyncRho(rho_fp, rho_cp); + + // Push E and B from {n} to {n+1} + // (And update guard cells immediately afterwards) +#ifdef WARPX_USE_PSATD + PushPSATD(dt[0]); + FillBoundaryE(); + FillBoundaryB(); +#else + EvolveF(0.5*dt[0], DtType::FirstHalf); + EvolveB(0.5*dt[0]); // We now have B^{n+1/2} + FillBoundaryB(); + EvolveE(dt[0]); // We now have E^{n+1} + FillBoundaryE(); + EvolveF(0.5*dt[0], DtType::SecondHalf); + EvolveB(0.5*dt[0]); // We now have B^{n+1} + if (do_pml) { + DampPML(); + FillBoundaryE(); + } + FillBoundaryB(); +#endif +} + +void +WarpX::OneStep_sub1 (Real curtime) +{ + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(finest_level == 1, "Must have exactly two levels"); + const int fine_lev = 1; + const int coarse_lev = 0; + // i) + PushParticlesandDepose(fine_lev, curtime); + RestrictCurrentFromFineToCoarsePatch(fine_lev); + SumBoundaryJ(fine_lev, PatchType::fine); + + amrex::Abort("xxxxx"); +#if 0 + // i) + PushParticlesandDepose( 1, cur_time ); + RestrictFineToCoarse(); + SumBoundaryJFine( 1 ); + EvolvePatchB( 1, 0.5*dt[1], fine_tag ); + EvolvePatchBPML( 1, 0.5*dt[1], fine_tag ); + FillBoundaryB(); + EvolvePatchE( 1, dt[1], fine_tag ); + EvolvePatchEPML( 1, dt[1], fine_tag ); + FillBoundaryE(); + EvolvePatchB( 1, 0.5*dt[1], fine_tag ); + EvolvePatchBPML( 1, 0.5*dt[1], fine_tag ); + // xxxxx + if (do_pml) { + DampPatchPML(1, fine_tag); // level 1 fp + FillBoundaryE(); // xxxxx make it level only + } + FillBoundaryB(); // xxxxx make it level only + // ii) + PushParticlesandDepose( 0, cur_time ); + CurrentFineToStore(); + AddCoarseToPatchBelow(); + SumBoundaryJCoarse(1); // xxxxx we may need to do nodal sync + SumBoundaryJFine(0); + EvolvePatchB( 1, dt[1], coarse_tag ); + EvolvePatchBPML( 1, dt[1], coarse_tag ); + FillBoundaryB(); + EvolvePatchE( 1, dt[1], coarse_tag ); + EvolvePatchEPML( 1, dt[1], coarse_tag ); + FillBoundaryE(); + EvolvePatchB( 0, 0.5*dt[0], fine_tag ); + EvolvePatchBPML( 0, 0.5*dt[0], fine_tag ); + FillBoundaryB(); + EvolvePatchE( 0, 0.5*dt[0], fine_tag ); + EvolvePatchEPML( 0, 0.5*dt[0], fine_tag ); + FillBoundaryE(); + // iii) + UpdateAuxilaryData(); + // iv) + PushParticlesandDepose( 1, cur_time+dt[1] ); + RestrictFineToCoarse(); + SumBoundaryJFine(1); + EvolvePatchB( 1, 0.5*dt[1], fine_tag ); + EvolvePatchBPML( 1, 0.5*dt[1], fine_tag ); + FillBoundaryB(); + EvolvePatchE( 1, dt[1], fine_tag ); + EvolvePatchEPML( 1, dt[1], fine_tag ); + FillBoundaryE(); + EvolvePatchB( 1, 0.5*dt[1], fine_tag ); + EvolvePatchBPML( 1, 0.5*dt[1], fine_tag ); + // xxxxx + if (do_pml) { + DampPatchPML(1, fine_tag); // level 1 fp + FillBoundaryE(); // xxxxx make it level only + } + FillBoundaryB(); + AddCoarseToPatchBelow(); + SumBoundaryJCoarse(1); + SumBoundaryJFine(0); + // v) + EvolvePatchE( 1, dt[1], coarse_tag ); + EvolvePatchEPML( 1, dt[1], coarse_tag ); + FillBoundaryE(); + EvolvePatchB( 1, dt[1], coarse_tag ); + EvolvePatchBPML( 1, dt[1], coarse_tag ); + // xxxxx + if (do_pml) { + DampPatchPML(1, coarse_tag); // level 1 cp + DampPatchPML(1, coarse_tag); // level 1 cp + FillBoundaryE(); // xxxxx make it level only + } + FillBoundaryB(); + + EvolvePatchE( 0, 0.5*dt[0], fine_tag ); + EvolvePatchEPML( 0, 0.5*dt[0], fine_tag ); + FillBoundaryE(); + EvolvePatchB( 0, 0.5*dt[0], fine_tag ); + EvolvePatchBPML( 0, 0.5*dt[0], fine_tag ); + // xxxxx + if (do_pml) { + DampPatchPML(0, fine_tag); // level 0 + FillBoundaryE(); // xxxxx make it level only + } + FillBoundaryB(); +#endif +} + void WarpX::EvolveB (Real dt) { @@ -239,107 +352,102 @@ void WarpX::EvolveB (int lev, Real dt) { BL_PROFILE("WarpX::EvolveB()"); + EvolveB(lev, PatchType::fine, dt); + if (lev > 0) { + EvolveB(lev, PatchType::coarse, dt); + } +} - // Parameters of the solver: order and mesh spacing - const int norder = 2; +void +WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real dt) +{ + const int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; + const std::array& dx = WarpX::CellSize(patch_level); + const std::array dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]}; + + MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz; + if (patch_type == PatchType::fine) + { + Ex = Efield_fp[lev][0].get(); + Ey = Efield_fp[lev][1].get(); + Ez = Efield_fp[lev][2].get(); + Bx = Bfield_fp[lev][0].get(); + By = Bfield_fp[lev][1].get(); + Bz = Bfield_fp[lev][2].get(); + } + else + { + Ex = Efield_cp[lev][0].get(); + Ey = Efield_cp[lev][1].get(); + Ez = Efield_cp[lev][2].get(); + Bx = Bfield_cp[lev][0].get(); + By = Bfield_cp[lev][1].get(); + Bz = Bfield_cp[lev][2].get(); + } - int npatches = (lev == 0) ? 1 : 2; + MultiFab* cost = costs[lev].get(); + const IntVect& rr = (lev < finestLevel()) ? refRatio(lev) : IntVect::TheUnitVector(); - for (int ipatch = 0; ipatch < npatches; ++ipatch) + // Loop through the grids, and over the tiles within each grid +#ifdef _OPENMP +#pragma omp parallel +#endif + for ( MFIter mfi(*Bx,true); mfi.isValid(); ++mfi ) { - int patch_level = (ipatch == 0) ? lev : lev-1; - const std::array& dx = WarpX::CellSize(patch_level); - const std::array dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]}; - - MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz; - if (ipatch == 0) - { - Ex = Efield_fp[lev][0].get(); - Ey = Efield_fp[lev][1].get(); - Ez = Efield_fp[lev][2].get(); - Bx = Bfield_fp[lev][0].get(); - By = Bfield_fp[lev][1].get(); - Bz = Bfield_fp[lev][2].get(); - } - else - { - Ex = Efield_cp[lev][0].get(); - Ey = Efield_cp[lev][1].get(); - Ez = Efield_cp[lev][2].get(); - Bx = Bfield_cp[lev][0].get(); - By = Bfield_cp[lev][1].get(); - Bz = Bfield_cp[lev][2].get(); + Real wt = ParallelDescriptor::second(); + + const Box& tbx = mfi.tilebox(Bx_nodal_flag); + const Box& tby = mfi.tilebox(By_nodal_flag); + const Box& tbz = mfi.tilebox(Bz_nodal_flag); + + // Call picsar routine for each tile + warpx_push_bvec( + tbx.loVect(), tbx.hiVect(), + tby.loVect(), tby.hiVect(), + tbz.loVect(), tbz.hiVect(), + BL_TO_FORTRAN_3D((*Ex)[mfi]), + BL_TO_FORTRAN_3D((*Ey)[mfi]), + BL_TO_FORTRAN_3D((*Ez)[mfi]), + BL_TO_FORTRAN_3D((*Bx)[mfi]), + BL_TO_FORTRAN_3D((*By)[mfi]), + BL_TO_FORTRAN_3D((*Bz)[mfi]), + &dtsdx[0], &dtsdx[1], &dtsdx[2], + &WarpX::maxwell_fdtd_solver_id); + + if (cost) { + Box cbx = mfi.tilebox(IntVect{AMREX_D_DECL(0,0,0)}); + if (patch_type == PatchType::coarse) cbx.refine(rr); + wt = (ParallelDescriptor::second() - wt) / cbx.d_numPts(); + (*cost)[mfi].plus(wt, cbx); } + } - MultiFab* cost = costs[lev].get(); - const IntVect& rr = (lev < finestLevel()) ? refRatio(lev) : IntVect::TheUnitVector(); + if (do_pml && pml[lev]->ok()) + { + const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp(); + const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); - // Loop through the grids, and over the tiles within each grid #ifdef _OPENMP #pragma omp parallel #endif - for ( MFIter mfi(*Bx,true); mfi.isValid(); ++mfi ) + for ( MFIter mfi(*pml_B[0],true); mfi.isValid(); ++mfi ) { - Real wt = ParallelDescriptor::second(); - const Box& tbx = mfi.tilebox(Bx_nodal_flag); const Box& tby = mfi.tilebox(By_nodal_flag); const Box& tbz = mfi.tilebox(Bz_nodal_flag); - - // Call picsar routine for each tile - warpx_push_bvec( + + WRPX_PUSH_PML_BVEC( tbx.loVect(), tbx.hiVect(), tby.loVect(), tby.hiVect(), tbz.loVect(), tbz.hiVect(), - BL_TO_FORTRAN_3D((*Ex)[mfi]), - BL_TO_FORTRAN_3D((*Ey)[mfi]), - BL_TO_FORTRAN_3D((*Ez)[mfi]), - BL_TO_FORTRAN_3D((*Bx)[mfi]), - BL_TO_FORTRAN_3D((*By)[mfi]), - BL_TO_FORTRAN_3D((*Bz)[mfi]), - &dtsdx[0], &dtsdx[1], &dtsdx[2], - &WarpX::maxwell_fdtd_solver_id); - - if (cost) { - Box cbx = mfi.tilebox(IntVect{AMREX_D_DECL(0,0,0)}); - if (ipatch == 1) cbx.refine(rr); - wt = (ParallelDescriptor::second() - wt) / cbx.d_numPts(); - (*cost)[mfi].plus(wt, cbx); - } - } - } - - if (do_pml && pml[lev]->ok()) - { - for (int ipatch = 0; ipatch < npatches; ++ipatch) - { - const auto& pml_B = (ipatch==0) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp(); - const auto& pml_E = (ipatch==0) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); - int patch_level = (ipatch == 0) ? lev : lev-1; - const std::array& dx = WarpX::CellSize(patch_level); - const std::array dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]}; -#ifdef _OPENMP -#pragma omp parallel -#endif - for ( MFIter mfi(*pml_B[0],true); mfi.isValid(); ++mfi ) - { - const Box& tbx = mfi.tilebox(Bx_nodal_flag); - const Box& tby = mfi.tilebox(By_nodal_flag); - const Box& tbz = mfi.tilebox(Bz_nodal_flag); - - WRPX_PUSH_PML_BVEC( - tbx.loVect(), tbx.hiVect(), - tby.loVect(), tby.hiVect(), - tbz.loVect(), tbz.hiVect(), - BL_TO_FORTRAN_3D((*pml_E[0])[mfi]), - BL_TO_FORTRAN_3D((*pml_E[1])[mfi]), - BL_TO_FORTRAN_3D((*pml_E[2])[mfi]), - BL_TO_FORTRAN_3D((*pml_B[0])[mfi]), - BL_TO_FORTRAN_3D((*pml_B[1])[mfi]), - BL_TO_FORTRAN_3D((*pml_B[2])[mfi]), - &dtsdx[0], &dtsdx[1], &dtsdx[2], - &WarpX::maxwell_fdtd_solver_id); - } + BL_TO_FORTRAN_3D((*pml_E[0])[mfi]), + BL_TO_FORTRAN_3D((*pml_E[1])[mfi]), + BL_TO_FORTRAN_3D((*pml_E[2])[mfi]), + BL_TO_FORTRAN_3D((*pml_B[0])[mfi]), + BL_TO_FORTRAN_3D((*pml_B[1])[mfi]), + BL_TO_FORTRAN_3D((*pml_B[2])[mfi]), + &dtsdx[0], &dtsdx[1], &dtsdx[2], + &WarpX::maxwell_fdtd_solver_id); } } } @@ -356,146 +464,142 @@ void WarpX::EvolveE (int lev, Real dt) { BL_PROFILE("WarpX::EvolveE()"); + EvolveE(lev, PatchType::fine, dt); + if (lev > 0) { + EvolveE(lev, PatchType::coarse, dt); + } +} +void +WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real dt) +{ // Parameters of the solver: order and mesh spacing - const int norder = 2; static constexpr Real c2 = PhysConst::c*PhysConst::c; const Real mu_c2_dt = (PhysConst::mu0*PhysConst::c*PhysConst::c) * dt; const Real c2dt = (PhysConst::c*PhysConst::c) * dt; - int npatches = (lev == 0) ? 1 : 2; + int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; + const std::array& dx = WarpX::CellSize(patch_level); + const std::array dtsdx_c2 {c2dt/dx[0], c2dt/dx[1], c2dt/dx[2]}; - for (int ipatch = 0; ipatch < npatches; ++ipatch) + MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz, *jx, *jy, *jz, *F; + if (patch_type == PatchType::fine) { - int patch_level = (ipatch == 0) ? lev : lev-1; - const std::array& dx = WarpX::CellSize(patch_level); - const std::array dtsdx_c2 {c2dt/dx[0], c2dt/dx[1], c2dt/dx[2]}; + Ex = Efield_fp[lev][0].get(); + Ey = Efield_fp[lev][1].get(); + Ez = Efield_fp[lev][2].get(); + Bx = Bfield_fp[lev][0].get(); + By = Bfield_fp[lev][1].get(); + Bz = Bfield_fp[lev][2].get(); + jx = current_fp[lev][0].get(); + jy = current_fp[lev][1].get(); + jz = current_fp[lev][2].get(); + F = F_fp[lev].get(); + } + else if (patch_type == PatchType::coarse) + { + Ex = Efield_cp[lev][0].get(); + Ey = Efield_cp[lev][1].get(); + Ez = Efield_cp[lev][2].get(); + Bx = Bfield_cp[lev][0].get(); + By = Bfield_cp[lev][1].get(); + Bz = Bfield_cp[lev][2].get(); + jx = current_cp[lev][0].get(); + jy = current_cp[lev][1].get(); + jz = current_cp[lev][2].get(); + F = F_cp[lev].get(); + } - MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz, *jx, *jy, *jz, *F; - if (ipatch == 0) - { - Ex = Efield_fp[lev][0].get(); - Ey = Efield_fp[lev][1].get(); - Ez = Efield_fp[lev][2].get(); - Bx = Bfield_fp[lev][0].get(); - By = Bfield_fp[lev][1].get(); - Bz = Bfield_fp[lev][2].get(); - jx = current_fp[lev][0].get(); - jy = current_fp[lev][1].get(); - jz = current_fp[lev][2].get(); - F = F_fp[lev].get(); + MultiFab* cost = costs[lev].get(); + const IntVect& rr = (lev < finestLevel()) ? refRatio(lev) : IntVect::TheUnitVector(); + + // Loop through the grids, and over the tiles within each grid +#ifdef _OPENMP +#pragma omp parallel +#endif + for ( MFIter mfi(*Ex,true); mfi.isValid(); ++mfi ) + { + Real wt = ParallelDescriptor::second(); + + const Box& tex = mfi.tilebox(Ex_nodal_flag); + const Box& tey = mfi.tilebox(Ey_nodal_flag); + const Box& tez = mfi.tilebox(Ez_nodal_flag); + + // Call picsar routine for each tile + warpx_push_evec( + tex.loVect(), tex.hiVect(), + tey.loVect(), tey.hiVect(), + tez.loVect(), tez.hiVect(), + BL_TO_FORTRAN_3D((*Ex)[mfi]), + BL_TO_FORTRAN_3D((*Ey)[mfi]), + BL_TO_FORTRAN_3D((*Ez)[mfi]), + BL_TO_FORTRAN_3D((*Bx)[mfi]), + BL_TO_FORTRAN_3D((*By)[mfi]), + BL_TO_FORTRAN_3D((*Bz)[mfi]), + BL_TO_FORTRAN_3D((*jx)[mfi]), + BL_TO_FORTRAN_3D((*jy)[mfi]), + BL_TO_FORTRAN_3D((*jz)[mfi]), + &mu_c2_dt, + &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2]); + + if (F) { + WRPX_CLEAN_EVEC(tex.loVect(), tex.hiVect(), + tey.loVect(), tey.hiVect(), + tez.loVect(), tez.hiVect(), + BL_TO_FORTRAN_3D((*Ex)[mfi]), + BL_TO_FORTRAN_3D((*Ey)[mfi]), + BL_TO_FORTRAN_3D((*Ez)[mfi]), + BL_TO_FORTRAN_3D((*F)[mfi]), + &dtsdx_c2[0]); } - else - { - Ex = Efield_cp[lev][0].get(); - Ey = Efield_cp[lev][1].get(); - Ez = Efield_cp[lev][2].get(); - Bx = Bfield_cp[lev][0].get(); - By = Bfield_cp[lev][1].get(); - Bz = Bfield_cp[lev][2].get(); - jx = current_cp[lev][0].get(); - jy = current_cp[lev][1].get(); - jz = current_cp[lev][2].get(); - F = F_cp[lev].get(); + + if (cost) { + Box cbx = mfi.tilebox(IntVect{AMREX_D_DECL(0,0,0)}); + if (patch_type == PatchType::coarse) cbx.refine(rr); + wt = (ParallelDescriptor::second() - wt) / cbx.d_numPts(); + (*cost)[mfi].plus(wt, cbx); } + } - MultiFab* cost = costs[lev].get(); - const IntVect& rr = (lev < finestLevel()) ? refRatio(lev) : IntVect::TheUnitVector(); + if (do_pml && pml[lev]->ok()) + { + if (F) pml[lev]->ExchangeF(patch_type, F); - // Loop through the grids, and over the tiles within each grid + const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp(); + const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); + const auto& pml_F = (patch_type == PatchType::fine) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp(); #ifdef _OPENMP #pragma omp parallel #endif - for ( MFIter mfi(*Ex,true); mfi.isValid(); ++mfi ) + for ( MFIter mfi(*pml_E[0],true); mfi.isValid(); ++mfi ) { - Real wt = ParallelDescriptor::second(); - const Box& tex = mfi.tilebox(Ex_nodal_flag); const Box& tey = mfi.tilebox(Ey_nodal_flag); const Box& tez = mfi.tilebox(Ez_nodal_flag); - - // Call picsar routine for each tile - warpx_push_evec( + + WRPX_PUSH_PML_EVEC( tex.loVect(), tex.hiVect(), tey.loVect(), tey.hiVect(), tez.loVect(), tez.hiVect(), - BL_TO_FORTRAN_3D((*Ex)[mfi]), - BL_TO_FORTRAN_3D((*Ey)[mfi]), - BL_TO_FORTRAN_3D((*Ez)[mfi]), - BL_TO_FORTRAN_3D((*Bx)[mfi]), - BL_TO_FORTRAN_3D((*By)[mfi]), - BL_TO_FORTRAN_3D((*Bz)[mfi]), - BL_TO_FORTRAN_3D((*jx)[mfi]), - BL_TO_FORTRAN_3D((*jy)[mfi]), - BL_TO_FORTRAN_3D((*jz)[mfi]), - &mu_c2_dt, + BL_TO_FORTRAN_3D((*pml_E[0])[mfi]), + BL_TO_FORTRAN_3D((*pml_E[1])[mfi]), + BL_TO_FORTRAN_3D((*pml_E[2])[mfi]), + BL_TO_FORTRAN_3D((*pml_B[0])[mfi]), + BL_TO_FORTRAN_3D((*pml_B[1])[mfi]), + BL_TO_FORTRAN_3D((*pml_B[2])[mfi]), &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2]); - - if (F) { - WRPX_CLEAN_EVEC(tex.loVect(), tex.hiVect(), - tey.loVect(), tey.hiVect(), - tez.loVect(), tez.hiVect(), - BL_TO_FORTRAN_3D((*Ex)[mfi]), - BL_TO_FORTRAN_3D((*Ey)[mfi]), - BL_TO_FORTRAN_3D((*Ez)[mfi]), - BL_TO_FORTRAN_3D((*F)[mfi]), - &dtsdx_c2[0]); - } - - if (cost) { - Box cbx = mfi.tilebox(IntVect{AMREX_D_DECL(0,0,0)}); - if (ipatch == 1) cbx.refine(rr); - wt = (ParallelDescriptor::second() - wt) / cbx.d_numPts(); - (*cost)[mfi].plus(wt, cbx); - } - } - } - - if (do_pml && pml[lev]->ok()) - { - pml[lev]->ExchangeF(F_fp[lev].get(), F_cp[lev].get()); - - for (int ipatch = 0; ipatch < npatches; ++ipatch) - { - const auto& pml_B = (ipatch==0) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp(); - const auto& pml_E = (ipatch==0) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); - const auto& pml_F = (ipatch==0) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp(); - int patch_level = (ipatch == 0) ? lev : lev-1; - const std::array& dx = WarpX::CellSize(patch_level); - const std::array dtsdx_c2 {c2dt/dx[0], c2dt/dx[1], c2dt/dx[2]}; -#ifdef _OPENMP -#pragma omp parallel -#endif - for ( MFIter mfi(*pml_E[0],true); mfi.isValid(); ++mfi ) + + if (pml_F) { - const Box& tex = mfi.tilebox(Ex_nodal_flag); - const Box& tey = mfi.tilebox(Ey_nodal_flag); - const Box& tez = mfi.tilebox(Ez_nodal_flag); - - WRPX_PUSH_PML_EVEC( + WRPX_PUSH_PML_EVEC_F( tex.loVect(), tex.hiVect(), tey.loVect(), tey.hiVect(), tez.loVect(), tez.hiVect(), BL_TO_FORTRAN_3D((*pml_E[0])[mfi]), BL_TO_FORTRAN_3D((*pml_E[1])[mfi]), BL_TO_FORTRAN_3D((*pml_E[2])[mfi]), - BL_TO_FORTRAN_3D((*pml_B[0])[mfi]), - BL_TO_FORTRAN_3D((*pml_B[1])[mfi]), - BL_TO_FORTRAN_3D((*pml_B[2])[mfi]), - &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2]); - - if (pml_F) - { - WRPX_PUSH_PML_EVEC_F( - tex.loVect(), tex.hiVect(), - tey.loVect(), tey.hiVect(), - tez.loVect(), tez.hiVect(), - BL_TO_FORTRAN_3D((*pml_E[0])[mfi]), - BL_TO_FORTRAN_3D((*pml_E[1])[mfi]), - BL_TO_FORTRAN_3D((*pml_E[2])[mfi]), - BL_TO_FORTRAN_3D((*pml_F )[mfi]), - &dtsdx_c2[0]); - } + BL_TO_FORTRAN_3D((*pml_F )[mfi]), + &dtsdx_c2[0]); } } } diff --git a/Source/WarpXPML.H b/Source/WarpXPML.H index 44088e7366b..4b0c3f79efe 100644 --- a/Source/WarpXPML.H +++ b/Source/WarpXPML.H @@ -83,6 +83,8 @@ private: amrex::Real dt_E = -1.e10; }; +enum struct PatchType : int; + class PML { public: @@ -112,6 +114,7 @@ public: const std::array& E_cp); void ExchangeF (amrex::MultiFab* F_fp, amrex::MultiFab* F_cp); + void ExchangeF (PatchType patch_type, amrex::MultiFab* F_fp); void FillBoundary (); void FillBoundaryE (); diff --git a/Source/WarpXPML.cpp b/Source/WarpXPML.cpp index 42bb913b3c1..c4aa568e14b 100644 --- a/Source/WarpXPML.cpp +++ b/Source/WarpXPML.cpp @@ -554,8 +554,18 @@ PML::ExchangeE (const std::array& E_fp, void PML::ExchangeF (MultiFab* F_fp, MultiFab* F_cp) { - if (pml_F_fp) Exchange(*pml_F_fp, *F_fp, *m_geom); - if (pml_F_cp) Exchange(*pml_F_cp, *F_cp, *m_cgeom); + ExchangeF(PatchType::fine, F_fp); + ExchangeF(PatchType::coarse, F_cp); +} + +void +PML::ExchangeF (PatchType patch_type, MultiFab* Fp) +{ + if (patch_type == PatchType::fine && pml_F_fp) { + Exchange(*pml_F_fp, *Fp, *m_geom); + } else if (patch_type == PatchType::coarse && pml_F_cp) { + Exchange(*pml_F_cp, *Fp, *m_cgeom); + } } void From 91efd23401e424608e6f4fda8095e177ba6704e7 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Tue, 11 Sep 2018 12:55:46 -0700 Subject: [PATCH 04/43] WIP --- Source/WarpX.H | 2 + Source/WarpXComm.cpp | 87 +++++++++++++++++++++++++++++++----------- Source/WarpXEvolve.cpp | 12 +++--- Source/WarpXPML.H | 8 +++- Source/WarpXPML.cpp | 76 ++++++++++++++++++++++++------------ 5 files changed, 132 insertions(+), 53 deletions(-) diff --git a/Source/WarpX.H b/Source/WarpX.H index d0fb8d9f7fc..535824f8680 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -240,6 +240,8 @@ private: void RestrictCurrentFromFineToCoarsePatch (int lev); void SumBoundaryJ (int lev, PatchType patch_type); + void FillBoundaryB (int lev, PatchType patch_type); + void FillBoundaryE (int lev, PatchType patch_type); void EvolveB (int lev, PatchType patch_type, amrex::Real dt); void EvolveE (int lev, PatchType patch_type, amrex::Real dt); diff --git a/Source/WarpXComm.cpp b/Source/WarpXComm.cpp index b5297e04ae3..b4499797c46 100644 --- a/Source/WarpXComm.cpp +++ b/Source/WarpXComm.cpp @@ -218,45 +218,86 @@ WarpX::FillBoundaryE () } void -WarpX::FillBoundaryE(int lev) +WarpX::FillBoundaryE (int lev) { - const auto& period = Geom(lev).periodicity(); - - if (do_pml && pml[lev]->ok()) { - ExchangeWithPmlE(lev); - pml[lev]->FillBoundaryE(); - } + FillBoundaryE(lev, PatchType::fine); + if (lev > 0) FillBoundaryE(lev, PatchType::coarse); +} - (*Efield_fp[lev][0]).FillBoundary( period ); - (*Efield_fp[lev][1]).FillBoundary( period ); - (*Efield_fp[lev][2]).FillBoundary( period ); +void +WarpX::FillBoundaryE (int lev, PatchType patch_type) +{ + if (patch_type == PatchType::fine) + { + if (do_pml && pml[lev]->ok()) + { + pml[lev]->ExchangeE(patch_type, + { Efield_fp[lev][0].get(), + Efield_fp[lev][1].get(), + Efield_fp[lev][2].get() }); + pml[lev]->FillBoundaryE(patch_type); + } - if (lev > 0) + const auto& period = Geom(lev).periodicity(); + (*Efield_fp[lev][0]).FillBoundary(period); + (*Efield_fp[lev][1]).FillBoundary(period); + (*Efield_fp[lev][2]).FillBoundary(period); + } + else if (patch_type == PatchType::coarse) { + if (do_pml && pml[lev]->ok()) + { + pml[lev]->ExchangeE(patch_type, + { Efield_cp[lev][0].get(), + Efield_cp[lev][1].get(), + Efield_cp[lev][2].get() }); + pml[lev]->FillBoundaryE(patch_type); + } + const auto& cperiod = Geom(lev-1).periodicity(); (*Efield_cp[lev][0]).FillBoundary(cperiod); (*Efield_cp[lev][1]).FillBoundary(cperiod); (*Efield_cp[lev][2]).FillBoundary(cperiod); - } + } } void -WarpX::FillBoundaryB(int lev) +WarpX::FillBoundaryB (int lev) { - const auto& period = Geom(lev).periodicity(); + FillBoundaryB(lev, PatchType::fine); + if (lev > 0) FillBoundaryB(lev, PatchType::coarse); +} - if (do_pml && pml[lev]->ok()) +void +WarpX::FillBoundaryB (int lev, PatchType patch_type) +{ + if (patch_type == PatchType::fine) { - ExchangeWithPmlB(lev); - pml[lev]->FillBoundaryB(); - } - - (*Bfield_fp[lev][0]).FillBoundary(period); - (*Bfield_fp[lev][1]).FillBoundary(period); - (*Bfield_fp[lev][2]).FillBoundary(period); + if (do_pml && pml[lev]->ok()) + { + pml[lev]->ExchangeB(patch_type, + { Bfield_fp[lev][0].get(), + Bfield_fp[lev][1].get(), + Bfield_fp[lev][2].get() }); + pml[lev]->FillBoundaryB(patch_type); + } - if (lev > 0) + const auto& period = Geom(lev).periodicity(); + (*Bfield_fp[lev][0]).FillBoundary(period); + (*Bfield_fp[lev][1]).FillBoundary(period); + (*Bfield_fp[lev][2]).FillBoundary(period); + } + else if (patch_type == PatchType::coarse) { + if (do_pml && pml[lev]->ok()) + { + pml[lev]->ExchangeB(patch_type, + { Bfield_cp[lev][0].get(), + Bfield_cp[lev][1].get(), + Bfield_cp[lev][2].get() }); + pml[lev]->FillBoundaryB(patch_type); + } + const auto& cperiod = Geom(lev-1).periodicity(); (*Bfield_cp[lev][0]).FillBoundary(cperiod); (*Bfield_cp[lev][1]).FillBoundary(cperiod); diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index 30415cf5c14..a8c2599cd04 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -250,15 +250,17 @@ WarpX::OneStep_sub1 (Real curtime) PushParticlesandDepose(fine_lev, curtime); RestrictCurrentFromFineToCoarsePatch(fine_lev); SumBoundaryJ(fine_lev, PatchType::fine); + EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]); + FillBoundaryB(fine_lev, PatchType::fine); amrex::Abort("xxxxx"); #if 0 // i) - PushParticlesandDepose( 1, cur_time ); - RestrictFineToCoarse(); - SumBoundaryJFine( 1 ); - EvolvePatchB( 1, 0.5*dt[1], fine_tag ); - EvolvePatchBPML( 1, 0.5*dt[1], fine_tag ); +// PushParticlesandDepose( 1, cur_time ); +// RestrictFineToCoarse(); +// SumBoundaryJFine( 1 ); +// EvolvePatchB( 1, 0.5*dt[1], fine_tag ); +// EvolvePatchBPML( 1, 0.5*dt[1], fine_tag ); FillBoundaryB(); EvolvePatchE( 1, dt[1], fine_tag ); EvolvePatchEPML( 1, dt[1], fine_tag ); diff --git a/Source/WarpXPML.H b/Source/WarpXPML.H index 4b0c3f79efe..e3f6a8e901b 100644 --- a/Source/WarpXPML.H +++ b/Source/WarpXPML.H @@ -112,13 +112,19 @@ public: const std::array& B_cp); void ExchangeE (const std::array& E_fp, const std::array& E_cp); + void ExchangeB (PatchType patch_type, + const std::array& Bp); + void ExchangeE (PatchType patch_type, + const std::array& Ep); void ExchangeF (amrex::MultiFab* F_fp, amrex::MultiFab* F_cp); - void ExchangeF (PatchType patch_type, amrex::MultiFab* F_fp); + void ExchangeF (PatchType patch_type, amrex::MultiFab* Fp); void FillBoundary (); void FillBoundaryE (); void FillBoundaryB (); + void FillBoundaryE (PatchType patch_type); + void FillBoundaryB (PatchType patch_type); bool ok () const { return m_ok; } diff --git a/Source/WarpXPML.cpp b/Source/WarpXPML.cpp index c4aa568e14b..3146a53d4b3 100644 --- a/Source/WarpXPML.cpp +++ b/Source/WarpXPML.cpp @@ -519,17 +519,25 @@ void PML::ExchangeB (const std::array& B_fp, const std::array& B_cp) { - if (pml_B_fp[0]) + ExchangeB(PatchType::fine, B_fp); + ExchangeB(PatchType::coarse, B_cp); +} + +void +PML::ExchangeB (PatchType patch_type, + const std::array& Bp) +{ + if (patch_type == PatchType::fine && pml_B_fp[0] && Bp[0]) { - Exchange(*pml_B_fp[0], *B_fp[0], *m_geom); - Exchange(*pml_B_fp[1], *B_fp[1], *m_geom); - Exchange(*pml_B_fp[2], *B_fp[2], *m_geom); + Exchange(*pml_B_fp[0], *Bp[0], *m_geom); + Exchange(*pml_B_fp[1], *Bp[1], *m_geom); + Exchange(*pml_B_fp[2], *Bp[2], *m_geom); } - if (B_cp[0] && pml_B_cp[0]) + else if (patch_type == PatchType::coarse && pml_B_cp[0] && Bp[0]) { - Exchange(*pml_B_cp[0], *B_cp[0], *m_cgeom); - Exchange(*pml_B_cp[1], *B_cp[1], *m_cgeom); - Exchange(*pml_B_cp[2], *B_cp[2], *m_cgeom); + Exchange(*pml_B_cp[0], *Bp[0], *m_cgeom); + Exchange(*pml_B_cp[1], *Bp[1], *m_cgeom); + Exchange(*pml_B_cp[2], *Bp[2], *m_cgeom); } } @@ -537,17 +545,25 @@ void PML::ExchangeE (const std::array& E_fp, const std::array& E_cp) { - if (pml_E_fp[0]) + ExchangeB(PatchType::fine, E_fp); + ExchangeB(PatchType::coarse, E_cp); +} + +void +PML::ExchangeE (PatchType patch_type, + const std::array& Ep) +{ + if (patch_type == PatchType::fine && pml_E_fp[0] && Ep[0]) { - Exchange(*pml_E_fp[0], *E_fp[0], *m_geom); - Exchange(*pml_E_fp[1], *E_fp[1], *m_geom); - Exchange(*pml_E_fp[2], *E_fp[2], *m_geom); + Exchange(*pml_E_fp[0], *Ep[0], *m_geom); + Exchange(*pml_E_fp[1], *Ep[1], *m_geom); + Exchange(*pml_E_fp[2], *Ep[2], *m_geom); } - if (E_cp[0] && pml_E_cp[0]) + else if (patch_type == PatchType::coarse && pml_E_cp[0] && Ep[0]) { - Exchange(*pml_E_cp[0], *E_cp[0], *m_cgeom); - Exchange(*pml_E_cp[1], *E_cp[1], *m_cgeom); - Exchange(*pml_E_cp[2], *E_cp[2], *m_cgeom); + Exchange(*pml_E_cp[0], *Ep[0], *m_cgeom); + Exchange(*pml_E_cp[1], *Ep[1], *m_cgeom); + Exchange(*pml_E_cp[2], *Ep[2], *m_cgeom); } } @@ -561,9 +577,9 @@ PML::ExchangeF (MultiFab* F_fp, MultiFab* F_cp) void PML::ExchangeF (PatchType patch_type, MultiFab* Fp) { - if (patch_type == PatchType::fine && pml_F_fp) { + if (patch_type == PatchType::fine && pml_F_fp && Fp) { Exchange(*pml_F_fp, *Fp, *m_geom); - } else if (patch_type == PatchType::coarse && pml_F_cp) { + } else if (patch_type == PatchType::coarse && pml_F_cp && Fp) { Exchange(*pml_F_cp, *Fp, *m_cgeom); } } @@ -621,15 +637,21 @@ PML::FillBoundary () void PML::FillBoundaryE () { - if (pml_E_fp[0] && pml_E_fp[0]->nGrowVect().max() > 0) + FillBoundaryE(PatchType::fine); + FillBoundaryE(PatchType::coarse); +} + +void +PML::FillBoundaryE (PatchType patch_type) +{ + if (patch_type == PatchType::fine && pml_E_fp[0] && pml_E_fp[0]->nGrowVect().max() > 0) { const auto& period = m_geom->periodicity(); pml_E_fp[0]->FillBoundary(period); pml_E_fp[1]->FillBoundary(period); pml_E_fp[2]->FillBoundary(period); } - - if (pml_E_cp[0] && pml_E_cp[0]->nGrowVect().max() > 0) + else if (patch_type == PatchType::coarse && pml_E_cp[0] && pml_E_cp[0]->nGrowVect().max() > 0) { const auto& period = m_cgeom->periodicity(); pml_E_cp[0]->FillBoundary(period); @@ -641,15 +663,21 @@ PML::FillBoundaryE () void PML::FillBoundaryB () { - if (pml_B_fp[0]) + FillBoundaryB(PatchType::fine); + FillBoundaryB(PatchType::coarse); +} + +void +PML::FillBoundaryB (PatchType patch_type) +{ + if (patch_type == PatchType::fine && pml_B_fp[0]) { const auto& period = m_geom->periodicity(); pml_B_fp[0]->FillBoundary(period); pml_B_fp[1]->FillBoundary(period); pml_B_fp[2]->FillBoundary(period); } - - if (pml_B_cp[0]) + else if (patch_type == PatchType::coarse && pml_B_cp[0]) { const auto& period = m_cgeom->periodicity(); pml_B_cp[0]->FillBoundary(period); From 3c0bbe5760499e748ac8a655a1570ea6143df104 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Tue, 11 Sep 2018 14:38:07 -0700 Subject: [PATCH 05/43] first pass --- Source/WarpX.H | 5 + Source/WarpX.cpp | 22 +++- Source/WarpXComm.cpp | 15 +++ Source/WarpXEvolve.cpp | 223 ++++++++++++++++++++--------------------- 4 files changed, 150 insertions(+), 115 deletions(-) diff --git a/Source/WarpX.H b/Source/WarpX.H index 535824f8680..3e4d4a9750b 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -239,12 +239,17 @@ private: void OneStep_sub1 (amrex::Real t); void RestrictCurrentFromFineToCoarsePatch (int lev); + void AddCurrentFromFineLevelandSumBoundary (int lev); + void StoreCurrent (int lev); + void RestoreCurrent (int lev); void SumBoundaryJ (int lev, PatchType patch_type); + void FillBoundaryB (int lev, PatchType patch_type); void FillBoundaryE (int lev, PatchType patch_type); void EvolveB (int lev, PatchType patch_type, amrex::Real dt); void EvolveE (int lev, PatchType patch_type, amrex::Real dt); + void DampPML (int lev, PatchType patch_type); #ifdef WARPX_DO_ELECTROSTATIC /// diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index e3e1b43d30e..491f4648d21 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -545,7 +545,7 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d current_fp[lev][1].reset( new MultiFab(amrex::convert(ba,jy_nodal_flag),dm,1,ngJ)); current_fp[lev][2].reset( new MultiFab(amrex::convert(ba,jz_nodal_flag),dm,1,ngJ)); - if (do_subcycling == 1 && lev == 1) { + if (do_subcycling == 1 && lev == 0) { current_store[lev][0].reset( new MultiFab(amrex::convert(ba,jx_nodal_flag),dm,1,ngJ)); current_store[lev][1].reset( new MultiFab(amrex::convert(ba,jy_nodal_flag),dm,1,ngJ)); current_store[lev][2].reset( new MultiFab(amrex::convert(ba,jz_nodal_flag),dm,1,ngJ)); @@ -822,3 +822,23 @@ WarpX::GatherBufferMasks (int lev) return GetInstance().getGatherBufferMasks(lev); } +void +WarpX::StoreCurrent (int lev) +{ + for (int idim = 0; idim < 3; ++idim) { + if (current_store[lev][idim]) { + MultiFab::Copy(*current_store[lev][idim], *current_fp[lev][idim], + 0, 0, 1, current_store[lev][idim]->nGrowVect()); + } + } +} + +void +WarpX::RestoreCurrent (int lev) +{ + for (int idim = 0; idim < 3; ++idim) { + if (current_store[lev][idim]) { + std::swap(current_fp[lev][idim], current_store[lev][idim]); + } + } +} diff --git a/Source/WarpXComm.cpp b/Source/WarpXComm.cpp index b4499797c46..f07ad9815bf 100644 --- a/Source/WarpXComm.cpp +++ b/Source/WarpXComm.cpp @@ -640,3 +640,18 @@ WarpX::SumBoundaryJ (int lev, PatchType patch_type) } } +void +WarpX::AddCurrentFromFineLevelandSumBoundary (int lev) +{ + const auto& period = Geom(lev).periodicity(); + for (int idim = 0; idim < 3; ++idim) { + MultiFab mf(current_fp[lev][idim]->boxArray(), + current_fp[lev][idim]->DistributionMap(), 1, 0); + mf.setVal(0.0); + mf.ParallelAdd(*current_cp[lev+1][idim], 0, 0, 1, + current_cp[lev+1][idim]->nGrowVect(), IntVect::TheZeroVector(), + period); + current_fp[lev][idim]->SumBoundary(period); + MultiFab::Copy(*current_fp[lev][idim], mf, 0, 0, 1, 0); + } +} diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index a8c2599cd04..61ee5f6b1ca 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -246,100 +246,95 @@ WarpX::OneStep_sub1 (Real curtime) AMREX_ALWAYS_ASSERT_WITH_MESSAGE(finest_level == 1, "Must have exactly two levels"); const int fine_lev = 1; const int coarse_lev = 0; + // i) PushParticlesandDepose(fine_lev, curtime); RestrictCurrentFromFineToCoarsePatch(fine_lev); SumBoundaryJ(fine_lev, PatchType::fine); + EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]); FillBoundaryB(fine_lev, PatchType::fine); - amrex::Abort("xxxxx"); -#if 0 - // i) -// PushParticlesandDepose( 1, cur_time ); -// RestrictFineToCoarse(); -// SumBoundaryJFine( 1 ); -// EvolvePatchB( 1, 0.5*dt[1], fine_tag ); -// EvolvePatchBPML( 1, 0.5*dt[1], fine_tag ); - FillBoundaryB(); - EvolvePatchE( 1, dt[1], fine_tag ); - EvolvePatchEPML( 1, dt[1], fine_tag ); - FillBoundaryE(); - EvolvePatchB( 1, 0.5*dt[1], fine_tag ); - EvolvePatchBPML( 1, 0.5*dt[1], fine_tag ); - // xxxxx + EvolveE(fine_lev, PatchType::fine, dt[fine_lev]); + FillBoundaryE(fine_lev, PatchType::fine); + + EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]); + if (do_pml) { - DampPatchPML(1, fine_tag); // level 1 fp - FillBoundaryE(); // xxxxx make it level only + DampPML(fine_lev, PatchType::fine); + FillBoundaryE(fine_lev, PatchType::fine); } - FillBoundaryB(); // xxxxx make it level only + + FillBoundaryB(fine_lev, PatchType::fine); + // ii) - PushParticlesandDepose( 0, cur_time ); - CurrentFineToStore(); - AddCoarseToPatchBelow(); - SumBoundaryJCoarse(1); // xxxxx we may need to do nodal sync - SumBoundaryJFine(0); - EvolvePatchB( 1, dt[1], coarse_tag ); - EvolvePatchBPML( 1, dt[1], coarse_tag ); - FillBoundaryB(); - EvolvePatchE( 1, dt[1], coarse_tag ); - EvolvePatchEPML( 1, dt[1], coarse_tag ); - FillBoundaryE(); - EvolvePatchB( 0, 0.5*dt[0], fine_tag ); - EvolvePatchBPML( 0, 0.5*dt[0], fine_tag ); - FillBoundaryB(); - EvolvePatchE( 0, 0.5*dt[0], fine_tag ); - EvolvePatchEPML( 0, 0.5*dt[0], fine_tag ); - FillBoundaryE(); + PushParticlesandDepose(coarse_lev, curtime); + StoreCurrent(coarse_lev); + AddCurrentFromFineLevelandSumBoundary(coarse_lev); + + EvolveB(fine_lev, PatchType::coarse, dt[fine_lev]); + FillBoundaryB(fine_lev, PatchType::coarse); + + EvolveE(fine_lev, PatchType::coarse, dt[fine_lev]); + FillBoundaryE(fine_lev, PatchType::coarse); + + EvolveB(coarse_lev, PatchType::fine, 0.5*dt[coarse_lev]); + FillBoundaryB(coarse_lev, PatchType::fine); + + EvolveE(coarse_lev, PatchType::fine, 0.5*dt[coarse_lev]); + FillBoundaryE(coarse_lev, PatchType::fine); + // iii) UpdateAuxilaryData(); + // iv) - PushParticlesandDepose( 1, cur_time+dt[1] ); - RestrictFineToCoarse(); - SumBoundaryJFine(1); - EvolvePatchB( 1, 0.5*dt[1], fine_tag ); - EvolvePatchBPML( 1, 0.5*dt[1], fine_tag ); - FillBoundaryB(); - EvolvePatchE( 1, dt[1], fine_tag ); - EvolvePatchEPML( 1, dt[1], fine_tag ); - FillBoundaryE(); - EvolvePatchB( 1, 0.5*dt[1], fine_tag ); - EvolvePatchBPML( 1, 0.5*dt[1], fine_tag ); - // xxxxx + PushParticlesandDepose(fine_lev, curtime+dt[fine_lev]); + RestrictCurrentFromFineToCoarsePatch(fine_lev); + SumBoundaryJ(fine_lev, PatchType::fine); + + EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]); + FillBoundaryB(fine_lev, PatchType::fine); + + EvolveE(fine_lev, PatchType::fine, dt[fine_lev]); + FillBoundaryE(fine_lev, PatchType::fine); + + EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]); + if (do_pml) { - DampPatchPML(1, fine_tag); // level 1 fp - FillBoundaryE(); // xxxxx make it level only + DampPML(fine_lev, PatchType::fine); + FillBoundaryE(fine_lev, PatchType::fine); } - FillBoundaryB(); - AddCoarseToPatchBelow(); - SumBoundaryJCoarse(1); - SumBoundaryJFine(0); + + FillBoundaryB(fine_lev, PatchType::fine); + // v) - EvolvePatchE( 1, dt[1], coarse_tag ); - EvolvePatchEPML( 1, dt[1], coarse_tag ); - FillBoundaryE(); - EvolvePatchB( 1, dt[1], coarse_tag ); - EvolvePatchBPML( 1, dt[1], coarse_tag ); - // xxxxx + RestoreCurrent(coarse_lev); + AddCurrentFromFineLevelandSumBoundary(coarse_lev); + + EvolveE(fine_lev, PatchType::coarse, dt[fine_lev]); + FillBoundaryE(fine_lev, PatchType::coarse); + + EvolveB(fine_lev, PatchType::coarse, dt[fine_lev]); + if (do_pml) { - DampPatchPML(1, coarse_tag); // level 1 cp - DampPatchPML(1, coarse_tag); // level 1 cp - FillBoundaryE(); // xxxxx make it level only + DampPML(fine_lev, PatchType::coarse); // do it twice + DampPML(fine_lev, PatchType::coarse); + FillBoundaryE(fine_lev, PatchType::coarse); } - FillBoundaryB(); + + FillBoundaryB(fine_lev, PatchType::coarse); + + EvolveE(coarse_lev, PatchType::fine, 0.5*dt[coarse_lev]); + FillBoundaryE(coarse_lev, PatchType::fine); - EvolvePatchE( 0, 0.5*dt[0], fine_tag ); - EvolvePatchEPML( 0, 0.5*dt[0], fine_tag ); - FillBoundaryE(); - EvolvePatchB( 0, 0.5*dt[0], fine_tag ); - EvolvePatchBPML( 0, 0.5*dt[0], fine_tag ); - // xxxxx + EvolveB(coarse_lev, PatchType::fine, 0.5*dt[coarse_lev]); + if (do_pml) { - DampPatchPML(0, fine_tag); // level 0 - FillBoundaryE(); // xxxxx make it level only + DampPML(coarse_lev, PatchType::fine); + FillBoundaryE(coarse_lev, PatchType::fine); } - FillBoundaryB(); -#endif + + FillBoundaryB(coarse_lev, PatchType::fine); } void @@ -685,8 +680,6 @@ WarpX::EvolveF (int lev, Real dt, DtType dt_type) void WarpX::DampPML () { - if (!do_pml) return; - for (int lev = 0; lev <= finest_level; ++lev) { DampPML(lev); } @@ -694,55 +687,57 @@ WarpX::DampPML () void WarpX::DampPML (int lev) +{ + DampPML(lev, PatchType::fine); + if (lev > 0) DampPML(lev, PatchType::coarse); +} + +void +WarpX::DampPML (int lev, PatchType patch_type) { if (!do_pml) return; BL_PROFILE("WarpX::DampPML()"); - int npatches = (lev == 0) ? 1 : 2; - - for (int ipatch = 0; ipatch < npatches; ++ipatch) + if (pml[lev]->ok()) { - if (pml[lev]->ok()) - { - const auto& pml_E = (ipatch==0) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); - const auto& pml_B = (ipatch==0) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp(); - const auto& pml_F = (ipatch==0) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp(); - const auto& sigba = (ipatch==0) ? pml[lev]->GetMultiSigmaBox_fp() - : pml[lev]->GetMultiSigmaBox_cp(); + const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); + const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp(); + const auto& pml_F = (patch_type == PatchType::fine) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp(); + const auto& sigba = (patch_type == PatchType::fine) ? pml[lev]->GetMultiSigmaBox_fp() + : pml[lev]->GetMultiSigmaBox_cp(); #ifdef _OPENMP #pragma omp parallel #endif - for ( MFIter mfi(*pml_E[0],true); mfi.isValid(); ++mfi ) - { - const Box& tex = mfi.tilebox(Ex_nodal_flag); - const Box& tey = mfi.tilebox(Ey_nodal_flag); - const Box& tez = mfi.tilebox(Ez_nodal_flag); - const Box& tbx = mfi.tilebox(Bx_nodal_flag); - const Box& tby = mfi.tilebox(By_nodal_flag); - const Box& tbz = mfi.tilebox(Bz_nodal_flag); - - WRPX_DAMP_PML(tex.loVect(), tex.hiVect(), - tey.loVect(), tey.hiVect(), - tez.loVect(), tez.hiVect(), - tbx.loVect(), tbx.hiVect(), - tby.loVect(), tby.hiVect(), - tbz.loVect(), tbz.hiVect(), - BL_TO_FORTRAN_3D((*pml_E[0])[mfi]), - BL_TO_FORTRAN_3D((*pml_E[1])[mfi]), - BL_TO_FORTRAN_3D((*pml_E[2])[mfi]), - BL_TO_FORTRAN_3D((*pml_B[0])[mfi]), - BL_TO_FORTRAN_3D((*pml_B[1])[mfi]), - BL_TO_FORTRAN_3D((*pml_B[2])[mfi]), - WRPX_PML_TO_FORTRAN(sigba[mfi])); - - if (pml_F) { - const Box& tnd = mfi.nodaltilebox(); - WRPX_DAMP_PML_F(tnd.loVect(), tnd.hiVect(), - BL_TO_FORTRAN_3D((*pml_F)[mfi]), - WRPX_PML_TO_FORTRAN(sigba[mfi])); - } + for ( MFIter mfi(*pml_E[0],true); mfi.isValid(); ++mfi ) + { + const Box& tex = mfi.tilebox(Ex_nodal_flag); + const Box& tey = mfi.tilebox(Ey_nodal_flag); + const Box& tez = mfi.tilebox(Ez_nodal_flag); + const Box& tbx = mfi.tilebox(Bx_nodal_flag); + const Box& tby = mfi.tilebox(By_nodal_flag); + const Box& tbz = mfi.tilebox(Bz_nodal_flag); + + WRPX_DAMP_PML(tex.loVect(), tex.hiVect(), + tey.loVect(), tey.hiVect(), + tez.loVect(), tez.hiVect(), + tbx.loVect(), tbx.hiVect(), + tby.loVect(), tby.hiVect(), + tbz.loVect(), tbz.hiVect(), + BL_TO_FORTRAN_3D((*pml_E[0])[mfi]), + BL_TO_FORTRAN_3D((*pml_E[1])[mfi]), + BL_TO_FORTRAN_3D((*pml_E[2])[mfi]), + BL_TO_FORTRAN_3D((*pml_B[0])[mfi]), + BL_TO_FORTRAN_3D((*pml_B[1])[mfi]), + BL_TO_FORTRAN_3D((*pml_B[2])[mfi]), + WRPX_PML_TO_FORTRAN(sigba[mfi])); + + if (pml_F) { + const Box& tnd = mfi.nodaltilebox(); + WRPX_DAMP_PML_F(tnd.loVect(), tnd.hiVect(), + BL_TO_FORTRAN_3D((*pml_F)[mfi]), + WRPX_PML_TO_FORTRAN(sigba[mfi])); } } } From 23231af9c9558ffa7da9dfefd477a6971970b78e Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Tue, 11 Sep 2018 15:41:18 -0700 Subject: [PATCH 06/43] fix a bug --- Source/WarpXComm.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/WarpXComm.cpp b/Source/WarpXComm.cpp index f07ad9815bf..e069e24415f 100644 --- a/Source/WarpXComm.cpp +++ b/Source/WarpXComm.cpp @@ -652,6 +652,6 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev) current_cp[lev+1][idim]->nGrowVect(), IntVect::TheZeroVector(), period); current_fp[lev][idim]->SumBoundary(period); - MultiFab::Copy(*current_fp[lev][idim], mf, 0, 0, 1, 0); + MultiFab::Add(*current_fp[lev][idim], mf, 0, 0, 1, 0); } } From 93dcb155de3fe4978a578109bc422ac8ed28e7ae Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Tue, 11 Sep 2018 17:31:04 -0700 Subject: [PATCH 07/43] dive cleaning --- Source/WarpX.H | 5 ++ Source/WarpXComm.cpp | 44 +++++++++++++++++ Source/WarpXEvolve.cpp | 108 ++++++++++++++++++++++++----------------- 3 files changed, 113 insertions(+), 44 deletions(-) diff --git a/Source/WarpX.H b/Source/WarpX.H index 3e4d4a9750b..a2eddb07822 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -244,11 +244,16 @@ private: void RestoreCurrent (int lev); void SumBoundaryJ (int lev, PatchType patch_type); + void RestrictRhoFromFineToCoarsePatch (int lev); + void SumBoundaryRho (int lev, PatchType patch_type); + void AddRhoFromFineLevelandSumBoundary (int lev, DtType dt_type); + void FillBoundaryB (int lev, PatchType patch_type); void FillBoundaryE (int lev, PatchType patch_type); void EvolveB (int lev, PatchType patch_type, amrex::Real dt); void EvolveE (int lev, PatchType patch_type, amrex::Real dt); + void EvolveF (int lev, PatchType patch_type, amrex::Real dt, DtType dt_type); void DampPML (int lev, PatchType patch_type); #ifdef WARPX_DO_ELECTROSTATIC diff --git a/Source/WarpXComm.cpp b/Source/WarpXComm.cpp index e069e24415f..b3989d25697 100644 --- a/Source/WarpXComm.cpp +++ b/Source/WarpXComm.cpp @@ -655,3 +655,47 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev) MultiFab::Add(*current_fp[lev][idim], mf, 0, 0, 1, 0); } } + +void +WarpX::RestrictRhoFromFineToCoarsePatch (int lev) +{ + if (rho_fp[lev]) { + rho_cp[lev]->setVal(0.0); + const IntVect& ref_ratio = refRatio(lev-1); + SyncRho(*rho_fp[lev], *rho_cp[lev], ref_ratio[0]); + } +} + +void +WarpX::SumBoundaryRho (int lev, PatchType patch_type) +{ + if (patch_type == PatchType::fine && rho_fp[lev]) + { + const auto& period = Geom(lev).periodicity(); + rho_fp[lev]->SumBoundary(period); + } + else if (patch_type == PatchType::coarse && rho_cp[lev]) + { + const auto& cperiod = Geom(lev-1).periodicity(); + rho_cp[lev]->SumBoundary(cperiod); + } +} + +void +WarpX::AddRhoFromFineLevelandSumBoundary(int lev, DtType dt_type) +{ + if (rho_fp[lev]) { + const auto& period = Geom(lev).periodicity(); + const int icomp = (dt_type == DtType::SecondHalf) ? 1 : 0; + const int ncomp = (dt_type == DtType::Full) ? 2 : 1; + MultiFab mf(rho_fp[lev]->boxArray(), + rho_fp[lev]->DistributionMap(), + ncomp, 0); + mf.setVal(0.0); + mf.ParallelAdd(*rho_cp[lev+1], icomp, 0, ncomp, + rho_cp[lev+1]->nGrowVect(), IntVect::TheZeroVector(), + period); + rho_fp[lev]->SumBoundary(period); + MultiFab::Add(*rho_fp[lev], mf, 0, icomp, ncomp, 0); + } +} diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index 61ee5f6b1ca..b46b8f9e0a3 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -243,6 +243,9 @@ WarpX::OneStep_nosub (Real cur_time) void WarpX::OneStep_sub1 (Real curtime) { + // TODO: we could save some charge depositions + // TODO: We need to redistribute particles after the first sub-step on the fine level + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(finest_level == 1, "Must have exactly two levels"); const int fine_lev = 1; const int coarse_lev = 0; @@ -250,15 +253,19 @@ WarpX::OneStep_sub1 (Real curtime) // i) PushParticlesandDepose(fine_lev, curtime); RestrictCurrentFromFineToCoarsePatch(fine_lev); + RestrictRhoFromFineToCoarsePatch(fine_lev); SumBoundaryJ(fine_lev, PatchType::fine); + SumBoundaryRho(fine_lev, PatchType::fine); EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]); + EvolveF(fine_lev, PatchType::fine, 0.5*dt[fine_lev], DtType::FirstHalf); FillBoundaryB(fine_lev, PatchType::fine); EvolveE(fine_lev, PatchType::fine, dt[fine_lev]); FillBoundaryE(fine_lev, PatchType::fine); EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]); + EvolveF(fine_lev, PatchType::fine, 0.5*dt[fine_lev], DtType::SecondHalf); if (do_pml) { DampPML(fine_lev, PatchType::fine); @@ -271,14 +278,17 @@ WarpX::OneStep_sub1 (Real curtime) PushParticlesandDepose(coarse_lev, curtime); StoreCurrent(coarse_lev); AddCurrentFromFineLevelandSumBoundary(coarse_lev); + AddRhoFromFineLevelandSumBoundary(coarse_lev, DtType::FirstHalf); EvolveB(fine_lev, PatchType::coarse, dt[fine_lev]); + EvolveF(fine_lev, PatchType::coarse, dt[fine_lev], DtType::FirstHalf); FillBoundaryB(fine_lev, PatchType::coarse); EvolveE(fine_lev, PatchType::coarse, dt[fine_lev]); FillBoundaryE(fine_lev, PatchType::coarse); EvolveB(coarse_lev, PatchType::fine, 0.5*dt[coarse_lev]); + EvolveF(coarse_lev, PatchType::fine, 0.5*dt[coarse_lev], DtType::FirstHalf); FillBoundaryB(coarse_lev, PatchType::fine); EvolveE(coarse_lev, PatchType::fine, 0.5*dt[coarse_lev]); @@ -291,14 +301,17 @@ WarpX::OneStep_sub1 (Real curtime) PushParticlesandDepose(fine_lev, curtime+dt[fine_lev]); RestrictCurrentFromFineToCoarsePatch(fine_lev); SumBoundaryJ(fine_lev, PatchType::fine); + SumBoundaryRho(fine_lev, PatchType::fine); EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]); + EvolveF(fine_lev, PatchType::fine, 0.5*dt[fine_lev], DtType::FirstHalf); FillBoundaryB(fine_lev, PatchType::fine); EvolveE(fine_lev, PatchType::fine, dt[fine_lev]); FillBoundaryE(fine_lev, PatchType::fine); EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]); + EvolveF(fine_lev, PatchType::fine, 0.5*dt[fine_lev], DtType::SecondHalf); if (do_pml) { DampPML(fine_lev, PatchType::fine); @@ -310,11 +323,13 @@ WarpX::OneStep_sub1 (Real curtime) // v) RestoreCurrent(coarse_lev); AddCurrentFromFineLevelandSumBoundary(coarse_lev); + AddRhoFromFineLevelandSumBoundary(coarse_lev, DtType::SecondHalf); EvolveE(fine_lev, PatchType::coarse, dt[fine_lev]); FillBoundaryE(fine_lev, PatchType::coarse); EvolveB(fine_lev, PatchType::coarse, dt[fine_lev]); + EvolveF(fine_lev, PatchType::coarse, dt[fine_lev], DtType::SecondHalf); if (do_pml) { DampPML(fine_lev, PatchType::coarse); // do it twice @@ -328,6 +343,7 @@ WarpX::OneStep_sub1 (Real curtime) FillBoundaryE(coarse_lev, PatchType::fine); EvolveB(coarse_lev, PatchType::fine, 0.5*dt[coarse_lev]); + EvolveF(coarse_lev, PatchType::fine, 0.5*dt[coarse_lev], DtType::SecondHalf); if (do_pml) { DampPML(coarse_lev, PatchType::fine); @@ -617,62 +633,66 @@ WarpX::EvolveF (int lev, Real dt, DtType dt_type) { if (!do_dive_cleaning) return; + EvolveF(lev, PatchType::fine, dt, dt_type); + if (lev > 0) EvolveF(lev, PatchType::coarse, dt, dt_type); +} + +void +WarpX::EvolveF (int lev, PatchType patch_type, Real dt, DtType dt_type) +{ + if (!do_dive_cleaning) return; + BL_PROFILE("WarpX::EvolveF()"); static constexpr Real c2inv = 1.0/(PhysConst::c*PhysConst::c); static constexpr Real mu_c2 = PhysConst::mu0*PhysConst::c*PhysConst::c; - int npatches = (lev == 0) ? 1 : 2; + int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; + const auto& dx = WarpX::CellSize(patch_level); + const std::array dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]}; - for (int ipatch = 0; ipatch < npatches; ++ipatch) + MultiFab *Ex, *Ey, *Ez, *rho, *F; + if (patch_type == PatchType::fine) { - int patch_level = (ipatch == 0) ? lev : lev-1; - const auto& dx = WarpX::CellSize(patch_level); - const std::array dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]}; - - MultiFab *Ex, *Ey, *Ez, *rho, *F; - if (ipatch == 0) - { - Ex = Efield_fp[lev][0].get(); - Ey = Efield_fp[lev][1].get(); - Ez = Efield_fp[lev][2].get(); - rho = rho_fp[lev].get(); - F = F_fp[lev].get(); - } - else - { - Ex = Efield_cp[lev][0].get(); - Ey = Efield_cp[lev][1].get(); - Ez = Efield_cp[lev][2].get(); - rho = rho_cp[lev].get(); - F = F_cp[lev].get(); - } - - const int rhocomp = (dt_type == DtType::FirstHalf) ? 0 : 1; - - MultiFab src(rho->boxArray(), rho->DistributionMap(), 1, 0); - ComputeDivE(src, 0, {Ex,Ey,Ez}, dx); - MultiFab::Saxpy(src, -mu_c2, *rho, rhocomp, 0, 1, 0); - MultiFab::Saxpy(*F, dt, src, 0, 0, 1, 0); + Ex = Efield_fp[lev][0].get(); + Ey = Efield_fp[lev][1].get(); + Ez = Efield_fp[lev][2].get(); + rho = rho_fp[lev].get(); + F = F_fp[lev].get(); + } + else + { + Ex = Efield_cp[lev][0].get(); + Ey = Efield_cp[lev][1].get(); + Ez = Efield_cp[lev][2].get(); + rho = rho_cp[lev].get(); + F = F_cp[lev].get(); + } + + const int rhocomp = (dt_type == DtType::FirstHalf) ? 0 : 1; - if (do_pml && pml[lev]->ok()) - { - const auto& pml_F = (ipatch==0) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp(); - const auto& pml_E = (ipatch==0) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); + MultiFab src(rho->boxArray(), rho->DistributionMap(), 1, 0); + ComputeDivE(src, 0, {Ex,Ey,Ez}, dx); + MultiFab::Saxpy(src, -mu_c2, *rho, rhocomp, 0, 1, 0); + MultiFab::Saxpy(*F, dt, src, 0, 0, 1, 0); + + if (do_pml && pml[lev]->ok()) + { + const auto& pml_F = (patch_type == PatchType::fine) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp(); + const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); #ifdef _OPENMP #pragma omp parallel #endif - for ( MFIter mfi(*pml_F,true); mfi.isValid(); ++mfi ) - { - const Box& bx = mfi.tilebox(); - WRPX_PUSH_PML_F(bx.loVect(), bx.hiVect(), - BL_TO_FORTRAN_ANYD((*pml_F )[mfi]), - BL_TO_FORTRAN_ANYD((*pml_E[0])[mfi]), - BL_TO_FORTRAN_ANYD((*pml_E[1])[mfi]), - BL_TO_FORTRAN_ANYD((*pml_E[2])[mfi]), - &dtsdx[0], &dtsdx[1], &dtsdx[2]); - } + for ( MFIter mfi(*pml_F,true); mfi.isValid(); ++mfi ) + { + const Box& bx = mfi.tilebox(); + WRPX_PUSH_PML_F(bx.loVect(), bx.hiVect(), + BL_TO_FORTRAN_ANYD((*pml_F )[mfi]), + BL_TO_FORTRAN_ANYD((*pml_E[0])[mfi]), + BL_TO_FORTRAN_ANYD((*pml_E[1])[mfi]), + BL_TO_FORTRAN_ANYD((*pml_E[2])[mfi]), + &dtsdx[0], &dtsdx[1], &dtsdx[2]); } } } From fe73bd37bec213331a8ebd3a3948c5a4f21583eb Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Wed, 12 Sep 2018 09:32:11 -0700 Subject: [PATCH 08/43] ParallelDescriptor::second -> amrex::second for more accurate clock --- Source/WarpXEvolve.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index 88237624702..b1c5cf213bd 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -407,7 +407,7 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real dt) #endif for ( MFIter mfi(*Bx,true); mfi.isValid(); ++mfi ) { - Real wt = ParallelDescriptor::second(); + Real wt = amrex::second(); const Box& tbx = mfi.tilebox(Bx_nodal_flag); const Box& tby = mfi.tilebox(By_nodal_flag); @@ -430,7 +430,7 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real dt) if (cost) { Box cbx = mfi.tilebox(IntVect{AMREX_D_DECL(0,0,0)}); if (patch_type == PatchType::coarse) cbx.refine(rr); - wt = (ParallelDescriptor::second() - wt) / cbx.d_numPts(); + wt = (amrex::second() - wt) / cbx.d_numPts(); (*cost)[mfi].plus(wt, cbx); } } @@ -532,7 +532,7 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real dt) #endif for ( MFIter mfi(*Ex,true); mfi.isValid(); ++mfi ) { - Real wt = ParallelDescriptor::second(); + Real wt = amrex::second(); const Box& tex = mfi.tilebox(Ex_nodal_flag); const Box& tey = mfi.tilebox(Ey_nodal_flag); @@ -569,7 +569,7 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real dt) if (cost) { Box cbx = mfi.tilebox(IntVect{AMREX_D_DECL(0,0,0)}); if (patch_type == PatchType::coarse) cbx.refine(rr); - wt = (ParallelDescriptor::second() - wt) / cbx.d_numPts(); + wt = (amrex::second() - wt) / cbx.d_numPts(); (*cost)[mfi].plus(wt, cbx); } } From 77f4e1a6ea0e9e34feca6d1b6111f07fe19ee093 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Wed, 12 Sep 2018 11:13:07 -0700 Subject: [PATCH 09/43] use_filter --- Source/WarpX.H | 9 ++++--- Source/WarpX.cpp | 8 +++--- Source/WarpXComm.cpp | 60 ++++++++++++++++++++++-------------------- Source/WarpXEvolve.cpp | 12 ++++----- 4 files changed, 46 insertions(+), 43 deletions(-) diff --git a/Source/WarpX.H b/Source/WarpX.H index a2eddb07822..fc7fd60920f 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -242,11 +242,11 @@ private: void AddCurrentFromFineLevelandSumBoundary (int lev); void StoreCurrent (int lev); void RestoreCurrent (int lev); - void SumBoundaryJ (int lev, PatchType patch_type); + void ApplyFilterandSumBoundaryJ (int lev, PatchType patch_type); void RestrictRhoFromFineToCoarsePatch (int lev); - void SumBoundaryRho (int lev, PatchType patch_type); - void AddRhoFromFineLevelandSumBoundary (int lev, DtType dt_type); + void ApplyFilterandSumBoundaryRho (int lev, PatchType patch_type, int icomp, int ncomp); + void AddRhoFromFineLevelandSumBoundary (int lev, int icomp, int ncomp); void FillBoundaryB (int lev, PatchType patch_type); void FillBoundaryE (int lev, PatchType patch_type); @@ -353,7 +353,8 @@ private: void LoadBalance (); - static void applyFilter (amrex::MultiFab& dstmf, const amrex::MultiFab& srcmf); + static void applyFilter (amrex::MultiFab& dstmf, const amrex::MultiFab& srcmf, + int scomp = 0, int dcomp = 0, int ncomp = 10000); void BuildBufferMasks (); const amrex::iMultiFab* getCurrentBufferMasks (int lev) const { diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 491f4648d21..e81bd1c4f56 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -748,9 +748,9 @@ WarpX::ComputeDivE (MultiFab& divE, int dcomp, } void -WarpX::applyFilter (MultiFab& dstmf, const MultiFab& srcmf) +WarpX::applyFilter (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dcomp, int ncomp) { - const int ncomp = srcmf.nComp(); + ncomp = std::min(ncomp, srcmf.nComp()); #ifdef _OPENMP #pragma omp parallel #endif @@ -765,10 +765,10 @@ WarpX::applyFilter (MultiFab& dstmf, const MultiFab& srcmf) tmpfab.resize(gbx,ncomp); tmpfab.setVal(0.0, gbx, 0, ncomp); const Box& ibx = gbx & srcfab.box(); - tmpfab.copy(srcfab, ibx, 0, ibx, 0, ncomp); + tmpfab.copy(srcfab, ibx, scomp, ibx, 0, ncomp); WRPX_FILTER(BL_TO_FORTRAN_BOX(tbx), BL_TO_FORTRAN_ANYD(tmpfab), - BL_TO_FORTRAN_ANYD(dstfab), + BL_TO_FORTRAN_N_ANYD(dstfab,dcomp), ncomp); } } diff --git a/Source/WarpXComm.cpp b/Source/WarpXComm.cpp index b3989d25697..88774f0d8a0 100644 --- a/Source/WarpXComm.cpp +++ b/Source/WarpXComm.cpp @@ -622,21 +622,22 @@ WarpX::RestrictCurrentFromFineToCoarsePatch (int lev) } void -WarpX::SumBoundaryJ (int lev, PatchType patch_type) +WarpX::ApplyFilterandSumBoundaryJ (int lev, PatchType patch_type) { - if (patch_type == PatchType::fine) - { - const auto& period = Geom(lev).periodicity(); - current_fp[lev][0]->SumBoundary(period); - current_fp[lev][1]->SumBoundary(period); - current_fp[lev][2]->SumBoundary(period); - } - else if (patch_type == PatchType::coarse) - { - const auto& cperiod = Geom(lev-1).periodicity(); - current_cp[lev][0]->SumBoundary(cperiod); - current_cp[lev][1]->SumBoundary(cperiod); - current_cp[lev][2]->SumBoundary(cperiod); + const int glev = (patch_type == PatchType::fine) ? lev : lev-1; + const auto& period = Geom(glev).periodicity(); + auto& j = (patch_type == PatchType::fine) ? current_fp[lev] : current_cp[lev]; + for (int idim = 0; idim < 3; ++idim) { + if (use_filter) { + IntVect ng = j[idim]->nGrowVect(); + ng += 1; + MultiFab jf(j[idim]->boxArray(), j[idim]->DistributionMap(), 1, ng); + applyFilter(jf, *j[idim]); + jf.SumBoundary(period); + MultiFab::Copy(*j[idim], jf, 0, 0, 1, 0); + } else { + j[idim]->SumBoundary(period); + } } } @@ -651,7 +652,7 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev) mf.ParallelAdd(*current_cp[lev+1][idim], 0, 0, 1, current_cp[lev+1][idim]->nGrowVect(), IntVect::TheZeroVector(), period); - current_fp[lev][idim]->SumBoundary(period); + ApplyFilterandSumBoundaryJ(lev, PatchType::fine); MultiFab::Add(*current_fp[lev][idim], mf, 0, 0, 1, 0); } } @@ -667,27 +668,28 @@ WarpX::RestrictRhoFromFineToCoarsePatch (int lev) } void -WarpX::SumBoundaryRho (int lev, PatchType patch_type) +WarpX::ApplyFilterandSumBoundaryRho (int lev, PatchType patch_type, int icomp, int ncomp) { - if (patch_type == PatchType::fine && rho_fp[lev]) - { - const auto& period = Geom(lev).periodicity(); - rho_fp[lev]->SumBoundary(period); - } - else if (patch_type == PatchType::coarse && rho_cp[lev]) - { - const auto& cperiod = Geom(lev-1).periodicity(); - rho_cp[lev]->SumBoundary(cperiod); + const int glev = (patch_type == PatchType::fine) ? lev : lev-1; + const auto& period = Geom(glev).periodicity(); + auto& r = (patch_type == PatchType::fine) ? rho_fp[lev] : rho_cp[lev]; + if (use_filter) { + IntVect ng = r->nGrowVect(); + ng += 1; + MultiFab rf(r->boxArray(), r->DistributionMap(), ncomp, ng); + applyFilter(rf, *r, icomp, 0, ncomp); + rf.SumBoundary(period); + MultiFab::Copy(*r, rf, 0, icomp, ncomp, 0); + } else { + r->SumBoundary(icomp, ncomp, period); } } void -WarpX::AddRhoFromFineLevelandSumBoundary(int lev, DtType dt_type) +WarpX::AddRhoFromFineLevelandSumBoundary(int lev, int icomp, int ncomp) { if (rho_fp[lev]) { const auto& period = Geom(lev).periodicity(); - const int icomp = (dt_type == DtType::SecondHalf) ? 1 : 0; - const int ncomp = (dt_type == DtType::Full) ? 2 : 1; MultiFab mf(rho_fp[lev]->boxArray(), rho_fp[lev]->DistributionMap(), ncomp, 0); @@ -695,7 +697,7 @@ WarpX::AddRhoFromFineLevelandSumBoundary(int lev, DtType dt_type) mf.ParallelAdd(*rho_cp[lev+1], icomp, 0, ncomp, rho_cp[lev+1]->nGrowVect(), IntVect::TheZeroVector(), period); - rho_fp[lev]->SumBoundary(period); + ApplyFilterandSumBoundaryRho(lev, PatchType::fine, icomp, ncomp); MultiFab::Add(*rho_fp[lev], mf, 0, icomp, ncomp, 0); } } diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index b1c5cf213bd..796ed3a2a8b 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -254,8 +254,8 @@ WarpX::OneStep_sub1 (Real curtime) PushParticlesandDepose(fine_lev, curtime); RestrictCurrentFromFineToCoarsePatch(fine_lev); RestrictRhoFromFineToCoarsePatch(fine_lev); - SumBoundaryJ(fine_lev, PatchType::fine); - SumBoundaryRho(fine_lev, PatchType::fine); + ApplyFilterandSumBoundaryJ(fine_lev, PatchType::fine); + ApplyFilterandSumBoundaryRho(fine_lev, PatchType::fine, 0, 2); EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]); EvolveF(fine_lev, PatchType::fine, 0.5*dt[fine_lev], DtType::FirstHalf); @@ -278,7 +278,7 @@ WarpX::OneStep_sub1 (Real curtime) PushParticlesandDepose(coarse_lev, curtime); StoreCurrent(coarse_lev); AddCurrentFromFineLevelandSumBoundary(coarse_lev); - AddRhoFromFineLevelandSumBoundary(coarse_lev, DtType::FirstHalf); + AddRhoFromFineLevelandSumBoundary(coarse_lev, 0, 1); EvolveB(fine_lev, PatchType::coarse, dt[fine_lev]); EvolveF(fine_lev, PatchType::coarse, dt[fine_lev], DtType::FirstHalf); @@ -300,8 +300,8 @@ WarpX::OneStep_sub1 (Real curtime) // iv) PushParticlesandDepose(fine_lev, curtime+dt[fine_lev]); RestrictCurrentFromFineToCoarsePatch(fine_lev); - SumBoundaryJ(fine_lev, PatchType::fine); - SumBoundaryRho(fine_lev, PatchType::fine); + ApplyFilterandSumBoundaryJ(fine_lev, PatchType::fine); + ApplyFilterandSumBoundaryRho(fine_lev, PatchType::fine, 0, 2); EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]); EvolveF(fine_lev, PatchType::fine, 0.5*dt[fine_lev], DtType::FirstHalf); @@ -323,7 +323,7 @@ WarpX::OneStep_sub1 (Real curtime) // v) RestoreCurrent(coarse_lev); AddCurrentFromFineLevelandSumBoundary(coarse_lev); - AddRhoFromFineLevelandSumBoundary(coarse_lev, DtType::SecondHalf); + AddRhoFromFineLevelandSumBoundary(coarse_lev, 1, 1); EvolveE(fine_lev, PatchType::coarse, dt[fine_lev]); FillBoundaryE(fine_lev, PatchType::coarse); From 009c74af02b500cf79f44a9596854d4275314d1d Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Fri, 14 Sep 2018 13:30:26 -0700 Subject: [PATCH 10/43] if test nullptr --- Source/WarpXComm.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/Source/WarpXComm.cpp b/Source/WarpXComm.cpp index 88774f0d8a0..5ff2ddf7632 100644 --- a/Source/WarpXComm.cpp +++ b/Source/WarpXComm.cpp @@ -673,6 +673,7 @@ WarpX::ApplyFilterandSumBoundaryRho (int lev, PatchType patch_type, int icomp, i const int glev = (patch_type == PatchType::fine) ? lev : lev-1; const auto& period = Geom(glev).periodicity(); auto& r = (patch_type == PatchType::fine) ? rho_fp[lev] : rho_cp[lev]; + if (r == nullptr) return; if (use_filter) { IntVect ng = r->nGrowVect(); ng += 1; From c893145f621155d6734579aaaed4db15ccfc3ca7 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Fri, 14 Sep 2018 15:39:16 -0700 Subject: [PATCH 11/43] filter coarse patch too --- Source/WarpXComm.cpp | 32 ++++++++++++++++++++++++++------ 1 file changed, 26 insertions(+), 6 deletions(-) diff --git a/Source/WarpXComm.cpp b/Source/WarpXComm.cpp index 5ff2ddf7632..03f5fa6767a 100644 --- a/Source/WarpXComm.cpp +++ b/Source/WarpXComm.cpp @@ -649,9 +649,20 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev) MultiFab mf(current_fp[lev][idim]->boxArray(), current_fp[lev][idim]->DistributionMap(), 1, 0); mf.setVal(0.0); - mf.ParallelAdd(*current_cp[lev+1][idim], 0, 0, 1, - current_cp[lev+1][idim]->nGrowVect(), IntVect::TheZeroVector(), - period); + if (use_filter) { + IntVect ng = current_cp[lev+1][idim]->nGrowVect(); + ng += 1; + MultiFab jf(current_cp[lev+1][idim]->boxArray(), + current_cp[lev+1][idim]->DistributionMap(), 1, ng); + applyFilter(jf, *current_cp[lev+1][idim]); + mf.ParallelAdd(jf, 0, 0, 1, ng, IntVect::TheZeroVector(), period); + MultiFab::Copy(*current_cp[lev+1][idim], jf, 0, 0, 1, 0); + + } else { + mf.ParallelAdd(*current_cp[lev+1][idim], 0, 0, 1, + current_cp[lev+1][idim]->nGrowVect(), IntVect::TheZeroVector(), + period); + } ApplyFilterandSumBoundaryJ(lev, PatchType::fine); MultiFab::Add(*current_fp[lev][idim], mf, 0, 0, 1, 0); } @@ -695,9 +706,18 @@ WarpX::AddRhoFromFineLevelandSumBoundary(int lev, int icomp, int ncomp) rho_fp[lev]->DistributionMap(), ncomp, 0); mf.setVal(0.0); - mf.ParallelAdd(*rho_cp[lev+1], icomp, 0, ncomp, - rho_cp[lev+1]->nGrowVect(), IntVect::TheZeroVector(), - period); + if (use_filter) { + IntVect ng = rho_cp[lev+1]->nGrowVect(); + ng += 1; + MultiFab rf(rho_cp[lev+1]->boxArray(), rho_cp[lev+1]->DistributionMap(), ncomp, ng); + applyFilter(rf, *rho_cp[lev+1], icomp, 0, ncomp); + mf.ParallelAdd(rf, 0, 0, ncomp, ng, IntVect::TheZeroVector(), period); + MultiFab::Copy(*rho_cp[lev+1], rf, 0, icomp, ncomp, 0); + } else { + mf.ParallelAdd(*rho_cp[lev+1], icomp, 0, ncomp, + rho_cp[lev+1]->nGrowVect(), IntVect::TheZeroVector(), + period); + } ApplyFilterandSumBoundaryRho(lev, PatchType::fine, icomp, ncomp); MultiFab::Add(*rho_fp[lev], mf, 0, icomp, ncomp, 0); } From dd6cb0466cc372215b2b3e7b2a8f56b91bd8929b Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Fri, 14 Sep 2018 17:15:51 -0700 Subject: [PATCH 12/43] missed a restriction of rho --- Source/WarpXEvolve.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index 796ed3a2a8b..8685b07d9a7 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -300,6 +300,7 @@ WarpX::OneStep_sub1 (Real curtime) // iv) PushParticlesandDepose(fine_lev, curtime+dt[fine_lev]); RestrictCurrentFromFineToCoarsePatch(fine_lev); + RestrictRhoFromFineToCoarsePatch(fine_lev); ApplyFilterandSumBoundaryJ(fine_lev, PatchType::fine); ApplyFilterandSumBoundaryRho(fine_lev, PatchType::fine, 0, 2); From 458c7b6c46673184e3f18a54d5d2806690eb9510 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Tue, 18 Sep 2018 14:32:17 -0700 Subject: [PATCH 13/43] Remove ApplyFilterAndSumBoundary from loop --- Source/WarpXComm.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Source/WarpXComm.cpp b/Source/WarpXComm.cpp index bd8572c48fb..d449eadf0c5 100644 --- a/Source/WarpXComm.cpp +++ b/Source/WarpXComm.cpp @@ -643,6 +643,8 @@ WarpX::ApplyFilterandSumBoundaryJ (int lev, PatchType patch_type) void WarpX::AddCurrentFromFineLevelandSumBoundary (int lev) { + ApplyFilterandSumBoundaryJ(lev, PatchType::fine); + const auto& period = Geom(lev).periodicity(); for (int idim = 0; idim < 3; ++idim) { MultiFab mf(current_fp[lev][idim]->boxArray(), @@ -662,7 +664,6 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev) current_cp[lev+1][idim]->nGrowVect(), IntVect::TheZeroVector(), period); } - ApplyFilterandSumBoundaryJ(lev, PatchType::fine); MultiFab::Add(*current_fp[lev][idim], mf, 0, 0, 1, 0); } } From 22c73414e110f77d10a59e723eef0f2c92fa03df Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Thu, 20 Sep 2018 16:02:02 -0700 Subject: [PATCH 14/43] add one extra ghost cell to J --- Source/WarpX.cpp | 14 +++++++++----- Source/WarpXEvolve.cpp | 1 - 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index e81bd1c4f56..aacb6d17407 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -513,17 +513,21 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d ngz = ngz_nonci; } + int ngJx = WarpX::nox + 1; + int ngJy = WarpX::noy + 1; + int ngJz = WarpX::noz + 1; + #if (AMREX_SPACEDIM == 3) IntVect ngE(ngx,ngy,ngz); - IntVect ngJ(ngx,ngy,ngz_nonci); - IntVect ngRho = ngJ + 1; // One extra ghost cell, so that it's safe to deposit charge density - // after pushing particle. + IntVect ngJ(ngJx,ngJy,ngJz); #elif (AMREX_SPACEDIM == 2) IntVect ngE(ngx,ngz); - IntVect ngJ(ngx,ngz_nonci); - IntVect ngRho = ngJ + 1; + IntVect ngJ(ngJx,ngJz); #endif + IntVect ngRho = ngJ+1; //One extra ghost cell, so that it's safe to deposit charge density + // after pushing particle. + if (n_current_deposition_buffer < 0) { n_current_deposition_buffer = ngJ.max(); } diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index 910d49c8a99..2754be9a5d4 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -263,7 +263,6 @@ void WarpX::OneStep_sub1 (Real curtime) { // TODO: we could save some charge depositions - // TODO: We need to redistribute particles after the first sub-step on the fine level AMREX_ALWAYS_ASSERT_WITH_MESSAGE(finest_level == 1, "Must have exactly two levels"); const int fine_lev = 1; From f5cce0a197652de08b74bfd0b7a74ec5c91c4a33 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Fri, 21 Sep 2018 16:28:34 -0700 Subject: [PATCH 15/43] summboundary on coarse patch and do nodal sync --- Source/WarpX.H | 2 ++ Source/WarpXComm.cpp | 49 ++++++++++++++++++++++++++++++++++++++++-- Source/WarpXEvolve.cpp | 4 ++++ 3 files changed, 53 insertions(+), 2 deletions(-) diff --git a/Source/WarpX.H b/Source/WarpX.H index fc7fd60920f..8502d1703a9 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -243,10 +243,12 @@ private: void StoreCurrent (int lev); void RestoreCurrent (int lev); void ApplyFilterandSumBoundaryJ (int lev, PatchType patch_type); + void NodalSyncJ (int lev, PatchType patch_type); void RestrictRhoFromFineToCoarsePatch (int lev); void ApplyFilterandSumBoundaryRho (int lev, PatchType patch_type, int icomp, int ncomp); void AddRhoFromFineLevelandSumBoundary (int lev, int icomp, int ncomp); + void NodalSyncRho (int lev, PatchType patch_type, int icomp, int ncomp); void FillBoundaryB (int lev, PatchType patch_type); void FillBoundaryE (int lev, PatchType patch_type); diff --git a/Source/WarpXComm.cpp b/Source/WarpXComm.cpp index d449eadf0c5..2099920bc7d 100644 --- a/Source/WarpXComm.cpp +++ b/Source/WarpXComm.cpp @@ -460,7 +460,7 @@ WarpX::SyncCurrent (const std::array& fine, int ref_ratio) { BL_ASSERT(ref_ratio == 2); - const IntVect& ng = fine[0]->nGrowVect()/ref_ratio; + const IntVect& ng = (fine[0]->nGrowVect() + 1) /ref_ratio; #ifdef _OPEMP #pragma omp parallel @@ -578,7 +578,7 @@ void WarpX::SyncRho (const MultiFab& fine, MultiFab& crse, int ref_ratio) { BL_ASSERT(ref_ratio == 2); - const IntVect& ng = fine.nGrowVect()/ref_ratio; + const IntVect& ng = (fine.nGrowVect()+1)/ref_ratio; const int nc = fine.nComp(); #ifdef _OPEMP @@ -657,15 +657,19 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev) current_cp[lev+1][idim]->DistributionMap(), 1, ng); applyFilter(jf, *current_cp[lev+1][idim]); mf.ParallelAdd(jf, 0, 0, 1, ng, IntVect::TheZeroVector(), period); + jf.SumBoundary(period); MultiFab::Copy(*current_cp[lev+1][idim], jf, 0, 0, 1, 0); } else { mf.ParallelAdd(*current_cp[lev+1][idim], 0, 0, 1, current_cp[lev+1][idim]->nGrowVect(), IntVect::TheZeroVector(), period); + current_cp[lev+1][idim]->SumBoundary(period); } MultiFab::Add(*current_fp[lev][idim], mf, 0, 0, 1, 0); } + NodalSyncJ(lev, PatchType::fine); + NodalSyncJ(lev+1, PatchType::coarse); } void @@ -712,13 +716,54 @@ WarpX::AddRhoFromFineLevelandSumBoundary(int lev, int icomp, int ncomp) MultiFab rf(rho_cp[lev+1]->boxArray(), rho_cp[lev+1]->DistributionMap(), ncomp, ng); applyFilter(rf, *rho_cp[lev+1], icomp, 0, ncomp); mf.ParallelAdd(rf, 0, 0, ncomp, ng, IntVect::TheZeroVector(), period); + rf.SumBoundary(icomp, ncomp, period); MultiFab::Copy(*rho_cp[lev+1], rf, 0, icomp, ncomp, 0); } else { mf.ParallelAdd(*rho_cp[lev+1], icomp, 0, ncomp, rho_cp[lev+1]->nGrowVect(), IntVect::TheZeroVector(), period); + rho_cp[lev+1]->SumBoundary(icomp, ncomp, period); } ApplyFilterandSumBoundaryRho(lev, PatchType::fine, icomp, ncomp); MultiFab::Add(*rho_fp[lev], mf, 0, icomp, ncomp, 0); + + NodalSyncRho(lev, PatchType::fine, icomp, ncomp); + NodalSyncRho(lev+1, PatchType::coarse, icomp, ncomp); + } +} + +void +WarpX::NodalSyncJ (int lev, PatchType patch_type) +{ + if (patch_type == PatchType::fine) + { + const auto& period = Geom(lev).periodicity(); + current_fp[lev][0]->OverrideSync(period); + current_fp[lev][1]->OverrideSync(period); + current_fp[lev][2]->OverrideSync(period); + } + else if (patch_type == PatchType::coarse) + { + const auto& cperiod = Geom(lev-1).periodicity(); + current_cp[lev][0]->OverrideSync(cperiod); + current_cp[lev][1]->OverrideSync(cperiod); + current_cp[lev][2]->OverrideSync(cperiod); + } +} + +void +WarpX::NodalSyncRho (int lev, PatchType patch_type, int icomp, int ncomp) +{ + if (patch_type == PatchType::fine && rho_fp[lev]) + { + const auto& period = Geom(lev).periodicity(); + MultiFab rhof(*rho_fp[lev], amrex::make_alias, icomp, ncomp); + rhof.OverrideSync(period); + } + else if (patch_type == PatchType::coarse && rho_cp[lev]) + { + const auto& cperiod = Geom(lev-1).periodicity(); + MultiFab rhoc(*rho_cp[lev], amrex::make_alias, icomp, ncomp); + rhoc.OverrideSync(cperiod); } } diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index 2754be9a5d4..d6939db972f 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -273,7 +273,9 @@ WarpX::OneStep_sub1 (Real curtime) RestrictCurrentFromFineToCoarsePatch(fine_lev); RestrictRhoFromFineToCoarsePatch(fine_lev); ApplyFilterandSumBoundaryJ(fine_lev, PatchType::fine); + NodalSyncJ(fine_lev, PatchType::fine); ApplyFilterandSumBoundaryRho(fine_lev, PatchType::fine, 0, 2); + NodalSyncRho(fine_lev, PatchType::fine, 0, 2); EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]); EvolveF(fine_lev, PatchType::fine, 0.5*dt[fine_lev], DtType::FirstHalf); @@ -320,7 +322,9 @@ WarpX::OneStep_sub1 (Real curtime) RestrictCurrentFromFineToCoarsePatch(fine_lev); RestrictRhoFromFineToCoarsePatch(fine_lev); ApplyFilterandSumBoundaryJ(fine_lev, PatchType::fine); + NodalSyncJ(fine_lev, PatchType::fine); ApplyFilterandSumBoundaryRho(fine_lev, PatchType::fine, 0, 2); + NodalSyncRho(fine_lev, PatchType::fine, 0, 2); EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]); EvolveF(fine_lev, PatchType::fine, 0.5*dt[fine_lev], DtType::FirstHalf); From d5dc055e1ea5e80fc17e0f0c2667f4dbdb26766b Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Fri, 21 Sep 2018 18:30:58 -0700 Subject: [PATCH 16/43] fix the component index in sumboundary --- Source/WarpXComm.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/WarpXComm.cpp b/Source/WarpXComm.cpp index 2099920bc7d..16c7108a0f5 100644 --- a/Source/WarpXComm.cpp +++ b/Source/WarpXComm.cpp @@ -716,7 +716,7 @@ WarpX::AddRhoFromFineLevelandSumBoundary(int lev, int icomp, int ncomp) MultiFab rf(rho_cp[lev+1]->boxArray(), rho_cp[lev+1]->DistributionMap(), ncomp, ng); applyFilter(rf, *rho_cp[lev+1], icomp, 0, ncomp); mf.ParallelAdd(rf, 0, 0, ncomp, ng, IntVect::TheZeroVector(), period); - rf.SumBoundary(icomp, ncomp, period); + rf.SumBoundary(0, ncomp, period); MultiFab::Copy(*rho_cp[lev+1], rf, 0, icomp, ncomp, 0); } else { mf.ParallelAdd(*rho_cp[lev+1], icomp, 0, ncomp, From b0a1347d47138155c966e11f2bddc7930aca1249 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Sun, 23 Sep 2018 18:17:18 -0700 Subject: [PATCH 17/43] current deposition buffer for subcycling method 1 --- Source/WarpX.cpp | 9 +++++---- Source/WarpXComm.cpp | 42 +++++++++++++++++++++++++++++++++++++++--- 2 files changed, 44 insertions(+), 7 deletions(-) diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index aacb6d17407..1d9a3c88c05 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -643,14 +643,14 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d Efield_cax[lev][1].reset( new MultiFab(amrex::convert(cba,Ey_nodal_flag),dm,1,ngE)); Efield_cax[lev][2].reset( new MultiFab(amrex::convert(cba,Ez_nodal_flag),dm,1,ngE)); - gather_buffer_masks[lev].reset( new iMultiFab(ba, dm, 1, 0) ); + gather_buffer_masks[lev].reset( new iMultiFab(ba, dm, 1, 1) ); } if (n_current_deposition_buffer > 0) { current_buf[lev][0].reset( new MultiFab(amrex::convert(cba,jx_nodal_flag),dm,1,ngJ)); current_buf[lev][1].reset( new MultiFab(amrex::convert(cba,jy_nodal_flag),dm,1,ngJ)); current_buf[lev][2].reset( new MultiFab(amrex::convert(cba,jz_nodal_flag),dm,1,ngJ)); - current_buffer_masks[lev].reset( new iMultiFab(ba, dm, 1, 0) ); + current_buffer_masks[lev].reset( new iMultiFab(ba, dm, 1, 1) ); } } @@ -790,7 +790,8 @@ WarpX::BuildBufferMasks () iMultiFab* bmasks = (ipass == 0) ? current_buffer_masks[lev].get() : gather_buffer_masks[lev].get(); if (bmasks) { - iMultiFab tmp(bmasks->boxArray(), bmasks->DistributionMap(), 1, ngbuffer); + const int ngtmp = ngbuffer + bmasks->nGrow(); + iMultiFab tmp(bmasks->boxArray(), bmasks->DistributionMap(), 1, ngtmp); const int covered = 1; const int notcovered = 0; const int physbnd = 1; @@ -803,7 +804,7 @@ WarpX::BuildBufferMasks () #endif for (MFIter mfi(*bmasks, true); mfi.isValid(); ++mfi) { - const Box& tbx = mfi.tilebox(); + const Box& tbx = mfi.growntilebox(); warpx_build_buffer_masks (BL_TO_FORTRAN_BOX(tbx), BL_TO_FORTRAN_ANYD((*bmasks)[mfi]), BL_TO_FORTRAN_ANYD(tmp[mfi]), diff --git a/Source/WarpXComm.cpp b/Source/WarpXComm.cpp index 16c7108a0f5..33e26e51a8e 100644 --- a/Source/WarpXComm.cpp +++ b/Source/WarpXComm.cpp @@ -645,12 +645,37 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev) { ApplyFilterandSumBoundaryJ(lev, PatchType::fine); + // When there are current buffers, unlike coarse patch, + // we don't care about the final state of them. + const auto& period = Geom(lev).periodicity(); for (int idim = 0; idim < 3; ++idim) { MultiFab mf(current_fp[lev][idim]->boxArray(), current_fp[lev][idim]->DistributionMap(), 1, 0); mf.setVal(0.0); - if (use_filter) { + if (use_filter && current_buf[lev+1][idim]) + { + // coarse patch of fine level + IntVect ng = current_cp[lev+1][idim]->nGrowVect(); + ng += 1; + MultiFab jfc(current_cp[lev+1][idim]->boxArray(), + current_cp[lev+1][idim]->DistributionMap(), 1, ng); + applyFilter(jfc, *current_cp[lev+1][idim]); + + // buffer patch of fine level + MultiFab jfb(current_buf[lev+1][idim]->boxArray(), + current_buf[lev+1][idim]->DistributionMap(), 1, ng); + applyFilter(jfb, *current_buf[lev+1][idim]); + + MultiFab::Copy(jfb, jfc, 0, 0, 1, ng); + mf.ParallelAdd(jfb, 0, 0, 1, ng, IntVect::TheZeroVector(), period); + + jfc.SumBoundary(period); + MultiFab::Copy(*current_cp[lev+1][idim], jfc, 0, 0, 1, 0); + } + else if (use_filter) // but no buffer + { + // coarse patch of fine level IntVect ng = current_cp[lev+1][idim]->nGrowVect(); ng += 1; MultiFab jf(current_cp[lev+1][idim]->boxArray(), @@ -659,8 +684,19 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev) mf.ParallelAdd(jf, 0, 0, 1, ng, IntVect::TheZeroVector(), period); jf.SumBoundary(period); MultiFab::Copy(*current_cp[lev+1][idim], jf, 0, 0, 1, 0); - - } else { + } + else if (current_buf[lev+1][idim]) // but no filter + { + MultiFab::Copy(*current_buf[lev+1][idim], + *current_cp [lev+1][idim], 0, 0, 1, + current_cp[lev+1][idim]->nGrow()); + mf.ParallelAdd(*current_buf[lev+1][idim], 0, 0, 1, + current_buf[lev+1][idim]->nGrowVect(), IntVect::TheZeroVector(), + period); + current_cp[lev+1][idim]->SumBoundary(period); + } + else // no filter, no buffer + { mf.ParallelAdd(*current_cp[lev+1][idim], 0, 0, 1, current_cp[lev+1][idim]->nGrowVect(), IntVect::TheZeroVector(), period); From c2645488c025dcdab745a8b42fe31847006484fe Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Sun, 23 Sep 2018 18:31:17 -0700 Subject: [PATCH 18/43] make sure E and B have enough ghost cells for subcycling --- Source/WarpX.cpp | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 1d9a3c88c05..bdf65236df7 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -498,24 +498,28 @@ WarpX::ClearLevel (int lev) void WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& dm) { + const int ngx_tmp = (maxLevel() > 0 && do_subcycling == 1) ? WarpX::nox+1 : WarpX::nox; + const int ngy_tmp = (maxLevel() > 0 && do_subcycling == 1) ? WarpX::noy+1 : WarpX::noy; + const int ngz_tmp = (maxLevel() > 0 && do_subcycling == 1) ? WarpX::noz+1 : WarpX::noz; + // Ex, Ey, Ez, Bx, By, and Bz have the same number of ghost cells. // jx, jy, jz and rho have the same number of ghost cells. // E and B have the same number of ghost cells as j and rho if NCI filter is not used, // but different number of ghost cells in z-direction if NCI filter is used. - int ngx = (WarpX::nox % 2) ? WarpX::nox+1 : WarpX::nox; // Always even number - int ngy = (WarpX::noy % 2) ? WarpX::noy+1 : WarpX::noy; // Always even number - int ngz_nonci = (WarpX::noz % 2) ? WarpX::noz+1 : WarpX::noz; // Always even number + int ngx = (ngx_tmp % 2) ? ngx_tmp+1 : ngx_tmp; // Always even number + int ngy = (ngy_tmp % 2) ? ngy_tmp+1 : ngy_tmp; // Always even number + int ngz_nonci = (ngz_tmp % 2) ? ngz_tmp+1 : ngz_tmp; // Always even number int ngz; if (warpx_use_fdtd_nci_corr()) { - int ng = WarpX::noz + (mypc->nstencilz_fdtd_nci_corr-1); + int ng = ngz_tmp + (mypc->nstencilz_fdtd_nci_corr-1); ngz = (ng % 2) ? ng+1 : ng; } else { ngz = ngz_nonci; } - int ngJx = WarpX::nox + 1; - int ngJy = WarpX::noy + 1; - int ngJz = WarpX::noz + 1; + int ngJx = ngx_tmp; + int ngJy = ngy_tmp; + int ngJz = ngz_tmp; #if (AMREX_SPACEDIM == 3) IntVect ngE(ngx,ngy,ngz); From 3425ed34ebc36ab6a96772d8e2fb1473861728d2 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Sun, 23 Sep 2018 21:14:01 -0700 Subject: [PATCH 19/43] add particles.deposit_on_main_grid option to allow deposit some particles to the main grid --- Source/ParticleContainer.H | 10 +++ Source/ParticleContainer.cpp | 11 ++++ Source/PhysicalParticleContainer.cpp | 93 +++++++++++++++------------- Source/WarpX.cpp | 4 ++ Source/WarpXParticleContainer.H | 2 + 5 files changed, 77 insertions(+), 43 deletions(-) diff --git a/Source/ParticleContainer.H b/Source/ParticleContainer.H index 2634ada815e..f643823527c 100644 --- a/Source/ParticleContainer.H +++ b/Source/ParticleContainer.H @@ -152,6 +152,14 @@ public: int nSpecies() const {return nspecies;} + int nSpeciesDepositOnMainGrid () const { + int r = 0; + for (int i : deposit_on_main_grid) { + if (i) ++r; + } + return r; + } + int L_lower_order_in_v() {return l_lower_order_in_v;} bool Use_fdtd_nci_corr() {return use_fdtd_nci_corr;} @@ -179,6 +187,8 @@ protected: std::vector species_names; + std::vector deposit_on_main_grid; + std::vector species_types; private: diff --git a/Source/ParticleContainer.cpp b/Source/ParticleContainer.cpp index c40e0c1bbe0..aa6759496c3 100644 --- a/Source/ParticleContainer.cpp +++ b/Source/ParticleContainer.cpp @@ -23,6 +23,7 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) else if (species_types[i] == PCTypes::RigidInjected) { allcontainers[i].reset(new RigidInjectedParticleContainer(amr_core, i, species_names[i])); } + allcontainers[i]->deposit_on_main_grid = deposit_on_main_grid[i]; } if (WarpX::use_laser) { allcontainers[n-1].reset(new LaserParticleContainer(amr_core,n-1)); @@ -44,6 +45,16 @@ MultiParticleContainer::ReadParameters () pp.getarr("species_names", species_names); BL_ASSERT(species_names.size() == nspecies); + deposit_on_main_grid.resize(nspecies, 0); + std::vector tmp; + pp.queryarr("deposit_on_main_grid", tmp); + for (auto const& name : tmp) { + auto it = std::find(species_names.begin(), species_names.end(), name); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(it != species_names.end(), "ERROR: species in particles.deposit_on_main_grid must be part of particles.species_names"); + int i = std::distance(species_names.begin(), it); + deposit_on_main_grid[i] = 1; + } + species_types.resize(nspecies, PCTypes::Physical); std::vector rigid_injected_species; diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp index 23e1bb19e82..6f8b05a70b3 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/PhysicalParticleContainer.cpp @@ -891,6 +891,10 @@ PhysicalParticleContainer::Evolve (int lev, } } + if (deposit_on_main_grid && lev > 0) { + nfine_current = 0; + } + if (nfine_current != np || nfine_gather != np) { particle_tmp.resize(np); @@ -1065,51 +1069,54 @@ PhysicalParticleContainer::Evolve (int lev, const std::array& xyzmin = xyzmin_tile; - tbx.grow(ngJ); - tby.grow(ngJ); - tbz.grow(ngJ); - - local_jx.resize(tbx); - local_jy.resize(tby); - local_jz.resize(tbz); - - local_jx = 0.0; - local_jy = 0.0; - local_jz = 0.0; - - jx_ptr = local_jx.dataPtr(); - jy_ptr = local_jy.dataPtr(); - jz_ptr = local_jz.dataPtr(); - - jxntot = local_jx.length(); - jyntot = local_jy.length(); - jzntot = local_jz.length(); - - warpx_current_deposition( - jx_ptr, &ngJ, jxntot, - jy_ptr, &ngJ, jyntot, - jz_ptr, &ngJ, jzntot, - &np_current, xp.data(), yp.data(), zp.data(), - uxp.data(), uyp.data(), uzp.data(), - giv.data(), wp.data(), &this->charge, - &xyzmin[0], &xyzmin[1], &xyzmin[2], - &dt, &dx[0], &dx[1], &dx[2], - &WarpX::nox,&WarpX::noy,&WarpX::noz, - &lvect,&WarpX::current_deposition_algo); - - BL_PROFILE_VAR_STOP(blp_pxr_cd); - - BL_PROFILE_VAR_START(blp_accumulate); - const int ncomp = 1; - amrex_atomic_accumulate_fab(BL_TO_FORTRAN_3D(local_jx), - BL_TO_FORTRAN_3D(jxfab), ncomp); + if (np_current > 0) + { + tbx.grow(ngJ); + tby.grow(ngJ); + tbz.grow(ngJ); + + local_jx.resize(tbx); + local_jy.resize(tby); + local_jz.resize(tbz); - amrex_atomic_accumulate_fab(BL_TO_FORTRAN_3D(local_jy), - BL_TO_FORTRAN_3D(jyfab), ncomp); + local_jx = 0.0; + local_jy = 0.0; + local_jz = 0.0; - amrex_atomic_accumulate_fab(BL_TO_FORTRAN_3D(local_jz), - BL_TO_FORTRAN_3D(jzfab), ncomp); - BL_PROFILE_VAR_STOP(blp_accumulate); + jx_ptr = local_jx.dataPtr(); + jy_ptr = local_jy.dataPtr(); + jz_ptr = local_jz.dataPtr(); + + jxntot = local_jx.length(); + jyntot = local_jy.length(); + jzntot = local_jz.length(); + + warpx_current_deposition( + jx_ptr, &ngJ, jxntot, + jy_ptr, &ngJ, jyntot, + jz_ptr, &ngJ, jzntot, + &np_current, xp.data(), yp.data(), zp.data(), + uxp.data(), uyp.data(), uzp.data(), + giv.data(), wp.data(), &this->charge, + &xyzmin[0], &xyzmin[1], &xyzmin[2], + &dt, &dx[0], &dx[1], &dx[2], + &WarpX::nox,&WarpX::noy,&WarpX::noz, + &lvect,&WarpX::current_deposition_algo); + + BL_PROFILE_VAR_STOP(blp_pxr_cd); + + BL_PROFILE_VAR_START(blp_accumulate); + const int ncomp = 1; + amrex_atomic_accumulate_fab(BL_TO_FORTRAN_3D(local_jx), + BL_TO_FORTRAN_3D(jxfab), ncomp); + + amrex_atomic_accumulate_fab(BL_TO_FORTRAN_3D(local_jy), + BL_TO_FORTRAN_3D(jyfab), ncomp); + + amrex_atomic_accumulate_fab(BL_TO_FORTRAN_3D(local_jz), + BL_TO_FORTRAN_3D(jzfab), ncomp); + BL_PROFILE_VAR_STOP(blp_accumulate); + } if (np_current < np) { diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index bdf65236df7..2644061fa41 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -532,6 +532,10 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d IntVect ngRho = ngJ+1; //One extra ghost cell, so that it's safe to deposit charge density // after pushing particle. + if (mypc->nSpeciesDepositOnMainGrid() && n_current_deposition_buffer == 0) { + n_current_deposition_buffer = 1; + } + if (n_current_deposition_buffer < 0) { n_current_deposition_buffer = ngJ.max(); } diff --git a/Source/WarpXParticleContainer.H b/Source/WarpXParticleContainer.H index 8eaff7241ae..dd94b52978f 100644 --- a/Source/WarpXParticleContainer.H +++ b/Source/WarpXParticleContainer.H @@ -179,6 +179,8 @@ protected: amrex::Real charge; amrex::Real mass; + bool deposit_on_main_grid = false; + static int do_not_push; }; From 9d306501e7b20d6d0b7ee976bab81592ca2b4772 Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Mon, 24 Sep 2018 10:07:45 -0700 Subject: [PATCH 20/43] adding option to plot F and rho --- Source/WarpX.H | 2 ++ Source/WarpX.cpp | 2 ++ Source/WarpXIO.cpp | 22 ++++++++++++++++++++++ 3 files changed, 26 insertions(+) diff --git a/Source/WarpX.H b/Source/WarpX.H index 8502d1703a9..20b36fd67c3 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -464,6 +464,8 @@ private: bool plot_proc_number = false; bool plot_dive = false; bool plot_divb = true; + bool plot_rho = false; + bool plot_F = false; bool plot_finepatch = false; bool plot_crsepatch = false; bool plot_raw_fields = false; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 2644061fa41..5ab5dd91175 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -338,6 +338,8 @@ WarpX::ReadParameters () pp.query("plot_proc_number" , plot_proc_number); pp.query("plot_dive" , plot_dive); pp.query("plot_divb" , plot_divb); + pp.query("plot_rho" , plot_rho); + pp.query("plot_F" , plot_F); if (maxLevel() > 0) { pp.query("plot_finepatch", plot_finepatch); diff --git a/Source/WarpXIO.cpp b/Source/WarpXIO.cpp index cf4c2e01b52..0102c3ccd80 100644 --- a/Source/WarpXIO.cpp +++ b/Source/WarpXIO.cpp @@ -468,6 +468,8 @@ WarpX::WritePlotFile () const + static_cast(plot_proc_number) + static_cast(plot_divb) + static_cast(plot_dive) + + static_cast(plot_rho) + + static_cast(plot_F) + static_cast(plot_finepatch)*6 + static_cast(plot_crsepatch)*6 + static_cast(costs[0] != nullptr); @@ -618,6 +620,26 @@ WarpX::WritePlotFile () const dcomp += 1; } + if (plot_rho) + { + amrex::average_node_to_cellcenter(*mf[lev], dcomp, *rho_fp[lev], 0, 1); + if (lev == 0) + { + varnames.push_back("rho"); + } + dcomp += 1; + } + + if (plot_F) + { + amrex::average_node_to_cellcenter(*mf[lev], dcomp, *F_fp[lev], 0, 1); + if (lev == 0) + { + varnames.push_back("F"); + } + dcomp += 1; + } + if (plot_finepatch) { PackPlotDataPtrs(srcmf, Efield_fp[lev]); From a6a10fa7eae751fc324775ae18eb7bfcfee843b1 Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Mon, 24 Sep 2018 12:52:23 -0700 Subject: [PATCH 21/43] also use a buffer region for charge deposition --- Source/LaserParticleContainer.H | 2 +- Source/LaserParticleContainer.cpp | 2 +- Source/ParticleContainer.H | 2 +- Source/ParticleContainer.cpp | 4 +- Source/PhysicalParticleContainer.H | 1 + Source/PhysicalParticleContainer.cpp | 110 ++++++++++++++++------ Source/RigidInjectedParticleContainer.H | 1 + Source/RigidInjectedParticleContainer.cpp | 4 +- Source/WarpX.H | 4 +- Source/WarpX.cpp | 4 + Source/WarpXComm.cpp | 71 +++++++++++++- Source/WarpXEvolve.cpp | 2 +- Source/WarpXParticleContainer.H | 2 +- Source/WarpXRegrid.cpp | 8 ++ 14 files changed, 172 insertions(+), 45 deletions(-) diff --git a/Source/LaserParticleContainer.H b/Source/LaserParticleContainer.H index b5dcf72a6ad..8fe012fc092 100644 --- a/Source/LaserParticleContainer.H +++ b/Source/LaserParticleContainer.H @@ -28,7 +28,7 @@ public: const amrex::MultiFab&, const amrex::MultiFab&, const amrex::MultiFab&, amrex::MultiFab& jx, amrex::MultiFab& jy, amrex::MultiFab& jz, amrex::MultiFab*, amrex::MultiFab*, amrex::MultiFab*, - amrex::MultiFab* rho, + 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) final; diff --git a/Source/LaserParticleContainer.cpp b/Source/LaserParticleContainer.cpp index 10a9117df74..abb0f0fe5b8 100644 --- a/Source/LaserParticleContainer.cpp +++ b/Source/LaserParticleContainer.cpp @@ -280,7 +280,7 @@ LaserParticleContainer::Evolve (int lev, const MultiFab&, const MultiFab&, const MultiFab&, MultiFab& jx, MultiFab& jy, MultiFab& jz, MultiFab*, MultiFab*, MultiFab*, - MultiFab* rho, + MultiFab* rho, MultiFab* crho, const MultiFab*, const MultiFab*, const MultiFab*, const MultiFab*, const MultiFab*, const MultiFab*, Real t, Real dt) diff --git a/Source/ParticleContainer.H b/Source/ParticleContainer.H index f643823527c..868ddc8d557 100644 --- a/Source/ParticleContainer.H +++ b/Source/ParticleContainer.H @@ -98,7 +98,7 @@ public: const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz, amrex::MultiFab& jx, amrex::MultiFab& jy, amrex::MultiFab& jz, amrex::MultiFab* cjx, amrex::MultiFab* cjy, amrex::MultiFab* cjz, - amrex::MultiFab* rho, + 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); diff --git a/Source/ParticleContainer.cpp b/Source/ParticleContainer.cpp index aa6759496c3..ac7b3800635 100644 --- a/Source/ParticleContainer.cpp +++ b/Source/ParticleContainer.cpp @@ -208,7 +208,7 @@ MultiParticleContainer::Evolve (int lev, const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz, MultiFab& jx, MultiFab& jy, MultiFab& jz, MultiFab* cjx, MultiFab* cjy, MultiFab* cjz, - MultiFab* rho, + 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) @@ -222,7 +222,7 @@ MultiParticleContainer::Evolve (int lev, if (rho) rho->setVal(0.0); for (auto& pc : allcontainers) { pc->Evolve(lev, Ex, Ey, Ez, Bx, By, Bz, jx, jy, jz, cjx, cjy, cjz, - rho, cEx, cEy, cEz, cBx, cBy, cBz, t, dt); + rho, crho, cEx, cEy, cEz, cBx, cBy, cBz, t, dt); } } diff --git a/Source/PhysicalParticleContainer.H b/Source/PhysicalParticleContainer.H index 1f3e28c9ce4..52e59713b10 100644 --- a/Source/PhysicalParticleContainer.H +++ b/Source/PhysicalParticleContainer.H @@ -50,6 +50,7 @@ public: amrex::MultiFab* cjy, amrex::MultiFab* cjz, amrex::MultiFab* rho, + amrex::MultiFab* crho, const amrex::MultiFab* cEx, const amrex::MultiFab* cEy, const amrex::MultiFab* cEz, diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp index 6f8b05a70b3..346c18037e6 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/PhysicalParticleContainer.cpp @@ -687,7 +687,7 @@ PhysicalParticleContainer::Evolve (int lev, const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz, MultiFab& jx, MultiFab& jy, MultiFab& jz, MultiFab* cjx, MultiFab* cjy, MultiFab* cjz, - MultiFab* rho, + 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) @@ -927,6 +927,8 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_STOP(blp_partition); } + const long np_current = (cjx) ? nfine_current : np; + // // copy data from particle container to temp arrays // @@ -943,40 +945,88 @@ PhysicalParticleContainer::Evolve (int lev, auto depositCharge = [&] (MultiFab* rhomf, int icomp) { long ngRho = rhomf->nGrow(); - Real* data_ptr; - const int *rholen; - FArrayBox& rhofab = (*rhomf)[pti]; Box tile_box = convert(pti.tilebox(), IntVect::TheUnitVector()); - Box grown_box; - const std::array& xyzmin = xyzmin_tile; - tile_box.grow(ngRho); - local_rho.resize(tile_box); - local_rho = 0.0; - data_ptr = local_rho.dataPtr(); - rholen = local_rho.length(); + const int *rholen; + + if (np_current > 0) + { + FArrayBox& rhofab = (*rhomf)[pti]; + const std::array& xyzmin = xyzmin_tile; + tile_box.grow(ngRho); + local_rho.resize(tile_box); + local_rho = 0.0; + data_ptr = local_rho.dataPtr(); + rholen = local_rho.length(); + +#if (AMREX_SPACEDIM == 3) + const long nx = rholen[0]-1-2*ngRho; + const long ny = rholen[1]-1-2*ngRho; + const long nz = rholen[2]-1-2*ngRho; +#else + const long nx = rholen[0]-1-2*ngRho; + const long ny = 0; + const long nz = rholen[1]-1-2*ngRho; +#endif + warpx_charge_deposition(data_ptr, &np, + xp.data(), yp.data(), zp.data(), wp.data(), + &this->charge, + &xyzmin[0], &xyzmin[1], &xyzmin[2], + &dx[0], &dx[1], &dx[2], &nx, &ny, &nz, + &ngRho, &ngRho, &ngRho, + &WarpX::nox,&WarpX::noy,&WarpX::noz, + &lvect, &WarpX::charge_deposition_algo); + + const int ncomp = 1; + amrex_atomic_accumulate_fab(BL_TO_FORTRAN_3D(local_rho), + BL_TO_FORTRAN_N_3D(rhofab,icomp), ncomp); + } + + if (np_current < np) + { + const IntVect& ref_ratio = WarpX::RefRatio(lev-1); + const Box& ctilebox = amrex::coarsen(pti.tilebox(), ref_ratio); + const std::array& cxyzmin_tile = WarpX::LowerCorner(ctilebox, lev-1); + + tile_box = amrex::convert(ctilebox, IntVect::TheUnitVector()); + tile_box.grow(ngRho); + + local_rho.resize(tile_box); + + local_rho = 0.0; + + data_ptr = local_rho.dataPtr(); + rholen = local_rho.length(); #if (AMREX_SPACEDIM == 3) - const long nx = rholen[0]-1-2*ngRho; - const long ny = rholen[1]-1-2*ngRho; - const long nz = rholen[2]-1-2*ngRho; + const long nx = rholen[0]-1-2*ngRho; + const long ny = rholen[1]-1-2*ngRho; + const long nz = rholen[2]-1-2*ngRho; #else - const long nx = rholen[0]-1-2*ngRho; - const long ny = 0; - const long nz = rholen[1]-1-2*ngRho; + const long nx = rholen[0]-1-2*ngRho; + const long ny = 0; + const long nz = rholen[1]-1-2*ngRho; #endif - warpx_charge_deposition(data_ptr, &np, - xp.data(), yp.data(), zp.data(), wp.data(), - &this->charge, - &xyzmin[0], &xyzmin[1], &xyzmin[2], - &dx[0], &dx[1], &dx[2], &nx, &ny, &nz, - &ngRho, &ngRho, &ngRho, - &WarpX::nox,&WarpX::noy,&WarpX::noz, - &lvect, &WarpX::charge_deposition_algo); - - const int ncomp = 1; - amrex_atomic_accumulate_fab(BL_TO_FORTRAN_3D(local_rho), - BL_TO_FORTRAN_N_3D(rhofab,icomp), ncomp); + + long ncrse = np - nfine_current; + warpx_charge_deposition(data_ptr, &ncrse, + xp.data() + nfine_current, + yp.data() + nfine_current, + zp.data() + nfine_current, + wp.data() + nfine_current, + &this->charge, + &cxyzmin_tile[0], &cxyzmin_tile[1], &cxyzmin_tile[2], + &cdx[0], &cdx[1], &cdx[2], &nx, &ny, &nz, + &ngRho, &ngRho, &ngRho, + &WarpX::nox,&WarpX::noy,&WarpX::noz, + &lvect, &WarpX::charge_deposition_algo); + + FArrayBox& crhofab = (*crho)[pti]; + + const int ncomp = 1; + amrex_atomic_accumulate_fab(BL_TO_FORTRAN_3D(local_rho), + BL_TO_FORTRAN_N_3D(crhofab,icomp), ncomp); + } }; if (rho) depositCharge(rho,0); @@ -1065,8 +1115,6 @@ PhysicalParticleContainer::Evolve (int lev, Box tbz = convert(pti.tilebox(), WarpX::jz_nodal_flag); Box gtbx, gtby, gtbz; - const long np_current = (cjx) ? nfine_current : np; - const std::array& xyzmin = xyzmin_tile; if (np_current > 0) diff --git a/Source/RigidInjectedParticleContainer.H b/Source/RigidInjectedParticleContainer.H index 18e9743ce33..f2440ad7a6c 100644 --- a/Source/RigidInjectedParticleContainer.H +++ b/Source/RigidInjectedParticleContainer.H @@ -32,6 +32,7 @@ public: amrex::MultiFab* cjy, amrex::MultiFab* cjz, amrex::MultiFab* rho, + amrex::MultiFab* crho, const amrex::MultiFab* cEx, const amrex::MultiFab* cEy, const amrex::MultiFab* cEz, diff --git a/Source/RigidInjectedParticleContainer.cpp b/Source/RigidInjectedParticleContainer.cpp index e5d4015801c..b3208dcba4b 100644 --- a/Source/RigidInjectedParticleContainer.cpp +++ b/Source/RigidInjectedParticleContainer.cpp @@ -307,7 +307,7 @@ RigidInjectedParticleContainer::Evolve (int lev, const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz, MultiFab& jx, MultiFab& jy, MultiFab& jz, MultiFab* cjx, MultiFab* cjy, MultiFab* cjz, - MultiFab* rho, + 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) @@ -332,7 +332,7 @@ RigidInjectedParticleContainer::Evolve (int lev, Bx, By, Bz, jx, jy, jz, cjx, cjy, cjz, - rho, + rho, crho, cEx, cEy, cEz, cBx, cBy, cBz, t, dt); diff --git a/Source/WarpX.H b/Source/WarpX.H index 20b36fd67c3..40d3c323d6b 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -409,8 +409,10 @@ private: amrex::Vector, 3 > > Bfield_cax; amrex::Vector > current_buffer_masks; amrex::Vector > gather_buffer_masks; - // If current deposition buffers are used + + // If charge/current deposition buffers are used amrex::Vector, 3 > > current_buf; + amrex::Vector > charge_buf; // div E cleaning int do_dive_cleaning = 0; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 5ab5dd91175..2c7bdecc55e 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -163,6 +163,7 @@ WarpX::WarpX () current_buffer_masks.resize(nlevs_max); gather_buffer_masks.resize(nlevs_max); current_buf.resize(nlevs_max); + charge_buf.resize(nlevs_max); pml.resize(nlevs_max); @@ -463,6 +464,8 @@ WarpX::ClearLevel (int lev) current_buf[lev][i].reset(); } + charge_buf[lev].reset(); + current_buffer_masks[lev].reset(); gather_buffer_masks[lev].reset(); @@ -660,6 +663,7 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d current_buf[lev][0].reset( new MultiFab(amrex::convert(cba,jx_nodal_flag),dm,1,ngJ)); current_buf[lev][1].reset( new MultiFab(amrex::convert(cba,jy_nodal_flag),dm,1,ngJ)); current_buf[lev][2].reset( new MultiFab(amrex::convert(cba,jz_nodal_flag),dm,1,ngJ)); + charge_buf[lev].reset( new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,2,ngRho)); current_buffer_masks[lev].reset( new iMultiFab(ba, dm, 1, 1) ); } } diff --git a/Source/WarpXComm.cpp b/Source/WarpXComm.cpp index 33e26e51a8e..031b4d14263 100644 --- a/Source/WarpXComm.cpp +++ b/Source/WarpXComm.cpp @@ -502,7 +502,8 @@ WarpX::SyncRho (amrex::Vector >& rhof, Vector > rho_f_g(finest_level+1); Vector > rho_c_g(finest_level+1); - + Vector > rho_buf_g(finest_level+1); + if (WarpX::use_filter) { for (int lev = 0; lev <= finest_level; ++lev) { const int ncomp = rhof[lev]->nComp(); @@ -524,6 +525,18 @@ WarpX::SyncRho (amrex::Vector >& rhof, applyFilter(*rho_c_g[lev], *rhoc[lev]); std::swap(rho_c_g[lev], rhoc[lev]); } + for (int lev = 1; lev <= finest_level; ++lev) { + if (charge_buf[lev]) { + const int ncomp = charge_buf[lev]->nComp(); + IntVect ng = charge_buf[lev]->nGrowVect(); + ng += 1; + rho_buf_g[lev].reset(new MultiFab(charge_buf[lev]->boxArray(), + charge_buf[lev]->DistributionMap(), + ncomp, ng)); + applyFilter(*rho_buf_g[lev], *charge_buf[lev]); + std::swap(*rho_buf_g[lev], *charge_buf[lev]); + } + } } // Sum up fine patch @@ -540,7 +553,14 @@ WarpX::SyncRho (amrex::Vector >& rhof, const int ncomp = rhoc[lev+1]->nComp(); const IntVect& ngsrc = rhoc[lev+1]->nGrowVect(); const IntVect ngdst = IntVect::TheZeroVector(); - rhof[lev]->copy(*rhoc[lev+1],0,0,ncomp,ngsrc,ngdst,period,FabArrayBase::ADD); + const MultiFab* crho = rhoc[lev+1].get(); + if (charge_buf[lev+1]) + { + MultiFab::Add(*charge_buf[lev+1], *rhoc[lev+1], 0, 0, ncomp, ngsrc); + crho = charge_buf[lev+1].get(); + } + + rhof[lev]->copy(*crho,0,0,ncomp,ngsrc,ngdst,period,FabArrayBase::ADD); } // Sum up coarse patch @@ -559,6 +579,13 @@ WarpX::SyncRho (amrex::Vector >& rhof, std::swap(rho_c_g[lev], rhoc[lev]); MultiFab::Copy(*rhoc[lev], *rho_c_g[lev], 0, 0, rhoc[lev]->nComp(), 0); } + for (int lev = 1; lev <= finest_level; ++lev) + { + if (rho_buf_g[lev]) { + std::swap(rho_buf_g[lev], charge_buf[lev]); + MultiFab::Copy(*charge_buf[lev], *rho_buf_g[lev], 0, 0, rhoc[lev]->nComp(), 0); + } + } } // sync shared nodal points @@ -746,7 +773,29 @@ WarpX::AddRhoFromFineLevelandSumBoundary(int lev, int icomp, int ncomp) rho_fp[lev]->DistributionMap(), ncomp, 0); mf.setVal(0.0); - if (use_filter) { + if (use_filter && charge_buf[lev]) + { + // coarse patch of fine level + IntVect ng = rho_cp[lev+1]->nGrowVect(); + const int ncomp = rho_cp[lev+1]->nComp(); + ng += 1; + MultiFab rhofc(rho_cp[lev+1]->boxArray(), + rho_cp[lev+1]->DistributionMap(), ncomp, ng); + applyFilter(rhofc, *rho_cp[lev+1]); + + // buffer patch of fine level + MultiFab rhofb(charge_buf[lev+1]->boxArray(), + charge_buf[lev+1]->DistributionMap(), ncomp, ng); + applyFilter(rhofb, *charge_buf[lev+1]); + + MultiFab::Copy(rhofb, rhofc, 0, 0, ncomp, ng); + mf.ParallelAdd(rhofb, 0, 0, ncomp, ng, IntVect::TheZeroVector(), period); + + rhofc.SumBoundary(period); + MultiFab::Copy(*rho_cp[lev+1], rhofc, 0, 0, ncomp, 0); + } + else if (use_filter) // but no buffer + { IntVect ng = rho_cp[lev+1]->nGrowVect(); ng += 1; MultiFab rf(rho_cp[lev+1]->boxArray(), rho_cp[lev+1]->DistributionMap(), ncomp, ng); @@ -754,7 +803,21 @@ WarpX::AddRhoFromFineLevelandSumBoundary(int lev, int icomp, int ncomp) mf.ParallelAdd(rf, 0, 0, ncomp, ng, IntVect::TheZeroVector(), period); rf.SumBoundary(0, ncomp, period); MultiFab::Copy(*rho_cp[lev+1], rf, 0, icomp, ncomp, 0); - } else { + } + else if (charge_buf[lev+1]) // but no filter + { + MultiFab::Copy(*charge_buf[lev+1], + *rho_cp[lev+1], 0, 0, + rho_cp[lev+1]->nComp(), + rho_cp[lev+1]->nGrow()); + mf.ParallelAdd(*charge_buf[lev+1], 0, 0, + rho_cp[lev+1]->nComp(), + charge_buf[lev+1]->nGrowVect(), IntVect::TheZeroVector(), + period); + rho_cp[lev+1]->SumBoundary(period); + } + else // no filter, no buffer + { mf.ParallelAdd(*rho_cp[lev+1], icomp, 0, ncomp, rho_cp[lev+1]->nGrowVect(), IntVect::TheZeroVector(), period); diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index 4a305b0f370..a40f93718a9 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -806,7 +806,7 @@ WarpX::PushParticlesandDepose (int lev, Real cur_time) *Bfield_aux[lev][0],*Bfield_aux[lev][1],*Bfield_aux[lev][2], *current_fp[lev][0],*current_fp[lev][1],*current_fp[lev][2], current_buf[lev][0].get(), current_buf[lev][1].get(), current_buf[lev][2].get(), - rho_fp[lev].get(), + 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]); diff --git a/Source/WarpXParticleContainer.H b/Source/WarpXParticleContainer.H index dd94b52978f..9f311c84c08 100644 --- a/Source/WarpXParticleContainer.H +++ b/Source/WarpXParticleContainer.H @@ -95,7 +95,7 @@ public: const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz, amrex::MultiFab& jx, amrex::MultiFab& jy, amrex::MultiFab& jz, amrex::MultiFab* cjx, amrex::MultiFab* cjy, amrex::MultiFab* cjz, - amrex::MultiFab* rho, + 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) = 0; diff --git a/Source/WarpXRegrid.cpp b/Source/WarpXRegrid.cpp index 360033635c5..744dd18e8e1 100644 --- a/Source/WarpXRegrid.cpp +++ b/Source/WarpXRegrid.cpp @@ -183,6 +183,14 @@ WarpX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa current_buf[lev][idim] = std::move(pmf); } } + if (charge_buf[lev]) + { + const IntVect& ng = charge_buf[lev]->nGrowVect(); + auto pmf = std::unique_ptr(new MultiFab(charge_buf[lev]->boxArray(), + dm, 1, ng)); + // pmf->ParallelCopy(*charge_buf[lev][idim], 0, 0, 1, ng, ng); + charge_buf[lev] = std::move(pmf); + } if (current_buffer_masks[lev]) { const IntVect& ng = current_buffer_masks[lev]->nGrowVect(); From c627cc8ff5c9b706e6980958cc693688e9f8aa48 Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Mon, 24 Sep 2018 13:12:32 -0700 Subject: [PATCH 22/43] only turn on the current buffers if dive cleaning is on. --- Source/WarpX.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 2c7bdecc55e..5f2c61ef578 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -663,7 +663,9 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d current_buf[lev][0].reset( new MultiFab(amrex::convert(cba,jx_nodal_flag),dm,1,ngJ)); current_buf[lev][1].reset( new MultiFab(amrex::convert(cba,jy_nodal_flag),dm,1,ngJ)); current_buf[lev][2].reset( new MultiFab(amrex::convert(cba,jz_nodal_flag),dm,1,ngJ)); - charge_buf[lev].reset( new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,2,ngRho)); + if (do_dive_cleaning) { + charge_buf[lev].reset( new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,2,ngRho)); + } current_buffer_masks[lev].reset( new iMultiFab(ba, dm, 1, 1) ); } } From fdaf457daeef55ee5a55241b666fe3930ec7d5d1 Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Mon, 24 Sep 2018 13:30:00 -0700 Subject: [PATCH 23/43] don't name this argument since it isn't used --- Source/LaserParticleContainer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/LaserParticleContainer.cpp b/Source/LaserParticleContainer.cpp index abb0f0fe5b8..92b64b28cdb 100644 --- a/Source/LaserParticleContainer.cpp +++ b/Source/LaserParticleContainer.cpp @@ -280,7 +280,7 @@ LaserParticleContainer::Evolve (int lev, const MultiFab&, const MultiFab&, const MultiFab&, MultiFab& jx, MultiFab& jy, MultiFab& jz, MultiFab*, MultiFab*, MultiFab*, - MultiFab* rho, MultiFab* crho, + MultiFab* rho, MultiFab*, const MultiFab*, const MultiFab*, const MultiFab*, const MultiFab*, const MultiFab*, const MultiFab*, Real t, Real dt) From 482f763dabfb957ceb9b181fc754cd239ee596fa Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Mon, 24 Sep 2018 13:30:13 -0700 Subject: [PATCH 24/43] also initialize the crho to zero --- Source/ParticleContainer.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/Source/ParticleContainer.cpp b/Source/ParticleContainer.cpp index ac7b3800635..66a5736ab0a 100644 --- a/Source/ParticleContainer.cpp +++ b/Source/ParticleContainer.cpp @@ -220,6 +220,7 @@ MultiParticleContainer::Evolve (int lev, 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); From cc64b39c7c46ab3ad016d720b821af0a7061637b Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Mon, 24 Sep 2018 13:50:30 -0700 Subject: [PATCH 25/43] fix typo --- Source/PhysicalParticleContainer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp index 346c18037e6..fe52916421d 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/PhysicalParticleContainer.cpp @@ -968,7 +968,7 @@ PhysicalParticleContainer::Evolve (int lev, const long ny = 0; const long nz = rholen[1]-1-2*ngRho; #endif - warpx_charge_deposition(data_ptr, &np, + warpx_charge_deposition(data_ptr, &np_current, xp.data(), yp.data(), zp.data(), wp.data(), &this->charge, &xyzmin[0], &xyzmin[1], &xyzmin[2], From cac001d5e8a662b68c574ae7a3c29cc149d99f8a Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Mon, 24 Sep 2018 15:23:38 -0700 Subject: [PATCH 26/43] fix nic corrector for non-subcycling --- Source/ParticleContainer.H | 4 ++-- Source/PhysicalParticleContainer.cpp | 12 +++++------ Source/WarpXInitData.cpp | 31 ++++++++++++++++------------ 3 files changed, 26 insertions(+), 21 deletions(-) diff --git a/Source/ParticleContainer.H b/Source/ParticleContainer.H index f643823527c..9618edf696c 100644 --- a/Source/ParticleContainer.H +++ b/Source/ParticleContainer.H @@ -177,8 +177,8 @@ public: // Number of coefficients for the stencil of the NCI corrector. // The stencil is applied in the z direction only. static constexpr int nstencilz_fdtd_nci_corr=5; - std::array fdtd_nci_stencilz_ex; - std::array fdtd_nci_stencilz_by; + amrex::Vector > fdtd_nci_stencilz_ex; + amrex::Vector > fdtd_nci_stencilz_by; protected: diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp index 6f8b05a70b3..91ec9e76136 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/PhysicalParticleContainer.cpp @@ -776,7 +776,7 @@ PhysicalParticleContainer::Evolve (int lev, WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ex), BL_TO_FORTRAN_ANYD(filtered_Ex), BL_TO_FORTRAN_ANYD(Ex[pti]), - mypc.fdtd_nci_stencilz_ex.data(), + mypc.fdtd_nci_stencilz_ex[lev].data(), &nstencilz_fdtd_nci_corr); exfab = &filtered_Ex; @@ -784,7 +784,7 @@ PhysicalParticleContainer::Evolve (int lev, WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ez), BL_TO_FORTRAN_ANYD(filtered_Ez), BL_TO_FORTRAN_ANYD(Ez[pti]), - mypc.fdtd_nci_stencilz_by.data(), + mypc.fdtd_nci_stencilz_by[lev].data(), &nstencilz_fdtd_nci_corr); ezfab = &filtered_Ez; @@ -792,7 +792,7 @@ PhysicalParticleContainer::Evolve (int lev, WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_By), BL_TO_FORTRAN_ANYD(filtered_By), BL_TO_FORTRAN_ANYD(By[pti]), - mypc.fdtd_nci_stencilz_by.data(), + mypc.fdtd_nci_stencilz_by[lev].data(), &nstencilz_fdtd_nci_corr); byfab = &filtered_By; @@ -801,7 +801,7 @@ PhysicalParticleContainer::Evolve (int lev, WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ey), BL_TO_FORTRAN_ANYD(filtered_Ey), BL_TO_FORTRAN_ANYD(Ey[pti]), - mypc.fdtd_nci_stencilz_ex.data(), + mypc.fdtd_nci_stencilz_ex[lev].data(), &nstencilz_fdtd_nci_corr); eyfab = &filtered_Ey; @@ -809,7 +809,7 @@ PhysicalParticleContainer::Evolve (int lev, WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bx), BL_TO_FORTRAN_ANYD(filtered_Bx), BL_TO_FORTRAN_ANYD(Bx[pti]), - mypc.fdtd_nci_stencilz_by.data(), + mypc.fdtd_nci_stencilz_by[lev].data(), &nstencilz_fdtd_nci_corr); bxfab = &filtered_Bx; @@ -817,7 +817,7 @@ PhysicalParticleContainer::Evolve (int lev, WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bz), BL_TO_FORTRAN_ANYD(filtered_Bz), BL_TO_FORTRAN_ANYD(Bz[pti]), - mypc.fdtd_nci_stencilz_ex.data(), + mypc.fdtd_nci_stencilz_ex[lev].data(), &nstencilz_fdtd_nci_corr); bzfab = &filtered_Bz; #endif diff --git a/Source/WarpXInitData.cpp b/Source/WarpXInitData.cpp index cb38038293c..8ca4a940918 100644 --- a/Source/WarpXInitData.cpp +++ b/Source/WarpXInitData.cpp @@ -134,20 +134,25 @@ WarpX::InitNCICorrector () { if (warpx_use_fdtd_nci_corr()) { - const Geometry& gm = Geom(finest_level); - const Real* dx = gm.CellSize(); - const int l_lower_order_in_v = warpx_l_lower_order_in_v(); - amrex::Real dz, cdtodz; - if (AMREX_SPACEDIM == 3){ - dz = dx[2]; - }else{ - dz = dx[1]; + mypc->fdtd_nci_stencilz_ex.resize(max_level+1); + mypc->fdtd_nci_stencilz_by.resize(max_level+1); + for (int lev = 0; lev <= max_level; ++lev) + { + const Geometry& gm = Geom(lev); + const Real* dx = gm.CellSize(); + const int l_lower_order_in_v = warpx_l_lower_order_in_v(); + amrex::Real dz, cdtodz; + if (AMREX_SPACEDIM == 3){ + dz = dx[2]; + }else{ + dz = dx[1]; + } + cdtodz = PhysConst::c * dt[lev] / dz; + WRPX_PXR_NCI_CORR_INIT( (mypc->fdtd_nci_stencilz_ex)[lev].data(), + (mypc->fdtd_nci_stencilz_by)[lev].data(), + mypc->nstencilz_fdtd_nci_corr, cdtodz, + l_lower_order_in_v); } - cdtodz = PhysConst::c * dt[finest_level] / dz; - WRPX_PXR_NCI_CORR_INIT( mypc->fdtd_nci_stencilz_ex.data(), - mypc->fdtd_nci_stencilz_by.data(), - mypc->nstencilz_fdtd_nci_corr, cdtodz, - l_lower_order_in_v); } } From 5aa0d85086c548fbe3965ecbe04e6686b67833f0 Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Mon, 24 Sep 2018 17:18:20 -0700 Subject: [PATCH 27/43] MultiFab::Copy -> MultiFab::Add --- Source/WarpXComm.cpp | 20 +++++++++----------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/Source/WarpXComm.cpp b/Source/WarpXComm.cpp index 031b4d14263..a6868a2c9ff 100644 --- a/Source/WarpXComm.cpp +++ b/Source/WarpXComm.cpp @@ -694,7 +694,7 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev) current_buf[lev+1][idim]->DistributionMap(), 1, ng); applyFilter(jfb, *current_buf[lev+1][idim]); - MultiFab::Copy(jfb, jfc, 0, 0, 1, ng); + MultiFab::Add(jfb, jfc, 0, 0, 1, ng); mf.ParallelAdd(jfb, 0, 0, 1, ng, IntVect::TheZeroVector(), period); jfc.SumBoundary(period); @@ -773,22 +773,21 @@ WarpX::AddRhoFromFineLevelandSumBoundary(int lev, int icomp, int ncomp) rho_fp[lev]->DistributionMap(), ncomp, 0); mf.setVal(0.0); - if (use_filter && charge_buf[lev]) + if (use_filter && charge_buf[lev+1]) { // coarse patch of fine level IntVect ng = rho_cp[lev+1]->nGrowVect(); - const int ncomp = rho_cp[lev+1]->nComp(); ng += 1; MultiFab rhofc(rho_cp[lev+1]->boxArray(), rho_cp[lev+1]->DistributionMap(), ncomp, ng); - applyFilter(rhofc, *rho_cp[lev+1]); + applyFilter(rhofc, *rho_cp[lev+1], icomp, 0, ncomp); // buffer patch of fine level MultiFab rhofb(charge_buf[lev+1]->boxArray(), charge_buf[lev+1]->DistributionMap(), ncomp, ng); - applyFilter(rhofb, *charge_buf[lev+1]); + applyFilter(rhofb, *charge_buf[lev+1], icomp, 0, ncomp); - MultiFab::Copy(rhofb, rhofc, 0, 0, ncomp, ng); + MultiFab::Add(rhofb, rhofc, 0, 0, ncomp, ng); mf.ParallelAdd(rhofb, 0, 0, ncomp, ng, IntVect::TheZeroVector(), period); rhofc.SumBoundary(period); @@ -807,14 +806,13 @@ WarpX::AddRhoFromFineLevelandSumBoundary(int lev, int icomp, int ncomp) else if (charge_buf[lev+1]) // but no filter { MultiFab::Copy(*charge_buf[lev+1], - *rho_cp[lev+1], 0, 0, - rho_cp[lev+1]->nComp(), + *rho_cp[lev+1], icomp, icomp, ncomp, rho_cp[lev+1]->nGrow()); - mf.ParallelAdd(*charge_buf[lev+1], 0, 0, - rho_cp[lev+1]->nComp(), + mf.ParallelAdd(*charge_buf[lev+1], icomp, 0, + ncomp, charge_buf[lev+1]->nGrowVect(), IntVect::TheZeroVector(), period); - rho_cp[lev+1]->SumBoundary(period); + rho_cp[lev+1]->SumBoundary(icomp, ncomp, period); } else // no filter, no buffer { From 8fc2f552c13de666f8c94b82ed8a2f0e6ca7bd00 Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Mon, 24 Sep 2018 17:37:49 -0700 Subject: [PATCH 28/43] MultiFab::Copy -> MultiFab::Add --- Source/WarpXComm.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/WarpXComm.cpp b/Source/WarpXComm.cpp index 33e26e51a8e..0604fbc5c34 100644 --- a/Source/WarpXComm.cpp +++ b/Source/WarpXComm.cpp @@ -667,7 +667,7 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev) current_buf[lev+1][idim]->DistributionMap(), 1, ng); applyFilter(jfb, *current_buf[lev+1][idim]); - MultiFab::Copy(jfb, jfc, 0, 0, 1, ng); + MultiFab::Add(jfb, jfc, 0, 0, 1, ng); mf.ParallelAdd(jfb, 0, 0, 1, ng, IntVect::TheZeroVector(), period); jfc.SumBoundary(period); From 7eda00a9e1704e1f8d1ece8148b2394a06955f36 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Mon, 24 Sep 2018 18:01:57 -0700 Subject: [PATCH 29/43] add missing FillBoundary --- Source/WarpXEvolve.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index 4a305b0f370..6dc33ac8f2c 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -92,6 +92,8 @@ WarpX::EvolveEM (int numsteps) } else { // Beyond one step, we have E^{n} and B^{n}. // Particles have p^{n-1/2} and x^{n}. + FillBoundaryE(); + FillBoundaryB(); UpdateAuxilaryData(); } From 1eac7fba545a14d6ec2d4fe70fa5e226c53ebecc Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Tue, 25 Sep 2018 10:31:13 -0700 Subject: [PATCH 30/43] change deposit charge lambda function for consistency --- Source/PhysicalParticleContainer.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp index 18188ed32c0..c1a98935408 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/PhysicalParticleContainer.cpp @@ -942,7 +942,7 @@ PhysicalParticleContainer::Evolve (int lev, long lvect = 8; - auto depositCharge = [&] (MultiFab* rhomf, int icomp) + auto depositCharge = [&] (MultiFab* rhomf, MultiFab* crhomf, int icomp) { long ngRho = rhomf->nGrow(); Real* data_ptr; @@ -1021,7 +1021,7 @@ PhysicalParticleContainer::Evolve (int lev, &WarpX::nox,&WarpX::noy,&WarpX::noz, &lvect, &WarpX::charge_deposition_algo); - FArrayBox& crhofab = (*crho)[pti]; + FArrayBox& crhofab = (*crhomf)[pti]; const int ncomp = 1; amrex_atomic_accumulate_fab(BL_TO_FORTRAN_3D(local_rho), @@ -1029,7 +1029,7 @@ PhysicalParticleContainer::Evolve (int lev, } }; - if (rho) depositCharge(rho,0); + if (rho) depositCharge(rho, crho, 0); if (! do_not_push) { @@ -1229,7 +1229,7 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_STOP(blp_copy); } - if (rho) depositCharge(rho,1); + if (rho) depositCharge(rho, crho, 1); if (cost) { const Box& tbx = pti.tilebox(); From 1a0ee5fa4155895020a9ebc7bdae1679b0739f66 Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Tue, 25 Sep 2018 15:17:32 -0700 Subject: [PATCH 31/43] use correct species names in the boosted frame diagnostics --- Source/ParticleContainer.H | 2 ++ Source/WarpXBoostedFrameDiagnostic.cpp | 17 ++++++++++------- 2 files changed, 12 insertions(+), 7 deletions(-) diff --git a/Source/ParticleContainer.H b/Source/ParticleContainer.H index efa8f0c9f42..1b955d55f0c 100644 --- a/Source/ParticleContainer.H +++ b/Source/ParticleContainer.H @@ -180,6 +180,8 @@ public: amrex::Vector > fdtd_nci_stencilz_ex; amrex::Vector > fdtd_nci_stencilz_by; + std::vector GetSpeciesNames() const { return species_names; } + protected: // Particle container types diff --git a/Source/WarpXBoostedFrameDiagnostic.cpp b/Source/WarpXBoostedFrameDiagnostic.cpp index f3bb1eb8bcf..96903f99204 100644 --- a/Source/WarpXBoostedFrameDiagnostic.cpp +++ b/Source/WarpXBoostedFrameDiagnostic.cpp @@ -54,6 +54,9 @@ void BoostedFrameDiagnostic::Flush(const Geometry& geom) VisMF::Header::Version current_version = VisMF::GetHeaderVersion(); VisMF::SetHeaderVersion(amrex::VisMF::Header::NoFabHeader_v1); + auto & mypc = WarpX::GetInstance().GetPartContainer(); + const std::vector species_names = mypc.GetSpeciesNames(); + for (int i = 0; i < N_snapshots_; ++i) { int i_lab = (snapshots_[i].current_z_lab - snapshots_[i].zmin_lab) / dz_lab_; @@ -84,11 +87,9 @@ void BoostedFrameDiagnostic::Flush(const Geometry& geom) } if (WarpX::do_boosted_frame_particles) { - auto & mypc = WarpX::GetInstance().GetPartContainer(); for (int j = 0; j < mypc.nSpecies(); ++j) { std::stringstream part_ss; - std::string species_name = "/particle" + std::to_string(j) + "/"; - part_ss << snapshots_[i].file_name + species_name; + part_ss << snapshots_[i].file_name + "/" + species_names[j] + "/"; writeParticleData(particles_buffer_[i][j], part_ss.str(), i_lab); } particles_buffer_[i].clear(); @@ -115,6 +116,8 @@ writeLabFrameData(const MultiFab* cell_centered_data, const Real zlo_boost = domain_z_boost.lo(boost_direction_); const Real zhi_boost = domain_z_boost.hi(boost_direction_); + const std::vector species_names = mypc.GetSpeciesNames(); + for (int i = 0; i < N_snapshots_; ++i) { const Real old_z_boost = snapshots_[i].current_z_boost; snapshots_[i].updateCurrentZPositions(t_boost, @@ -201,8 +204,7 @@ writeLabFrameData(const MultiFab* cell_centered_data, if (WarpX::do_boosted_frame_particles) { for (int j = 0; j < mypc.nSpecies(); ++j) { std::stringstream part_ss; - std::string species_name = "/particle" + std::to_string(j) + "/"; - part_ss << snapshots_[i].file_name + species_name; + part_ss << snapshots_[i].file_name + "/" + species_names[j] + "/"; writeParticleData(particles_buffer_[i][j], part_ss.str(), i_lab); } particles_buffer_[i].clear(); @@ -311,7 +313,7 @@ LabSnapShot(Real t_lab_in, Real zmin_lab_in, file_name = Concatenate("lab_frame_data/snapshot", file_num, 5); if (ParallelDescriptor::IOProcessor()) { - + if (!UtilCreateDirectory(file_name, 0755)) CreateDirectoryFailed(file_name); @@ -323,11 +325,12 @@ LabSnapShot(Real t_lab_in, Real zmin_lab_in, } auto & mypc = WarpX::GetInstance().GetPartContainer(); + const std::vector species_names = mypc.GetSpeciesNames(); int nspecies = mypc.nSpecies(); const std::string particles_prefix = "particle"; for(int i = 0; i < nspecies; ++i) { - const std::string fullpath = file_name + "/" + particles_prefix + std::to_string(i); + const std::string fullpath = file_name + "/" + species_names[i]; if (!UtilCreateDirectory(fullpath, 0755)) CreateDirectoryFailed(fullpath); } From 3e1b5b3e1b86fb2561e1ca0b0d8169644db99ab9 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Wed, 26 Sep 2018 09:34:54 -0700 Subject: [PATCH 32/43] patch for PML F --- Source/WarpX.H | 1 + Source/WarpXComm.cpp | 34 ++++++++++++++++++++++++---------- Source/WarpXPML.H | 1 + Source/WarpXPML.cpp | 12 +++++++++--- 4 files changed, 35 insertions(+), 13 deletions(-) diff --git a/Source/WarpX.H b/Source/WarpX.H index aced4a8e6e8..25a4a7bb654 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -254,6 +254,7 @@ private: void FillBoundaryB (int lev, PatchType patch_type); void FillBoundaryE (int lev, PatchType patch_type); + void FillBoundaryF (int lev, PatchType patch_type); void EvolveB (int lev, PatchType patch_type, amrex::Real dt); void EvolveE (int lev, PatchType patch_type, amrex::Real dt); diff --git a/Source/WarpXComm.cpp b/Source/WarpXComm.cpp index 93432daef29..2066b497529 100644 --- a/Source/WarpXComm.cpp +++ b/Source/WarpXComm.cpp @@ -324,22 +324,36 @@ WarpX::FillBoundaryB (int lev, PatchType patch_type) } void -WarpX::FillBoundaryF(int lev) +WarpX::FillBoundaryF (int lev) { - const auto& period = Geom(lev).periodicity(); + FillBoundaryF(lev, PatchType::fine); + if (lev > 0) FillBoundaryF(lev, PatchType::coarse); +} - if (do_pml && pml[lev]->ok()) +void +WarpX::FillBoundaryF (int lev, PatchType patch_type) +{ + if (patch_type == PatchType::fine && F_fp[lev]) { - ExchangeWithPmlF(lev); - pml[lev]->FillBoundaryF(); - } - - if (F_fp[lev]) F_fp[lev]->FillBoundary(period); + if (do_pml && pml[lev]->ok()) + { + pml[lev]->ExchangeF(patch_type, F_fp[lev].get()); + pml[lev]->FillBoundaryF(patch_type); + } - if (lev > 0) + const auto& period = Geom(lev).periodicity(); + F_fp[lev]->FillBoundary(period); + } + else if (patch_type == PatchType::coarse && F_cp[lev]) { + if (do_pml && pml[lev]->ok()) + { + pml[lev]->ExchangeF(patch_type, F_cp[lev].get()); + pml[lev]->FillBoundaryF(patch_type); + } + const auto& cperiod = Geom(lev-1).periodicity(); - if (F_cp[lev]) F_cp[lev]->FillBoundary(cperiod); + F_cp[lev]->FillBoundary(cperiod); } } diff --git a/Source/WarpXPML.H b/Source/WarpXPML.H index b45ea4091f0..841075366fd 100644 --- a/Source/WarpXPML.H +++ b/Source/WarpXPML.H @@ -126,6 +126,7 @@ public: void FillBoundaryF (); void FillBoundaryE (PatchType patch_type); void FillBoundaryB (PatchType patch_type); + void FillBoundaryF (PatchType patch_type); bool ok () const { return m_ok; } diff --git a/Source/WarpXPML.cpp b/Source/WarpXPML.cpp index 04b3c5bebf2..ea8c44080d8 100644 --- a/Source/WarpXPML.cpp +++ b/Source/WarpXPML.cpp @@ -691,13 +691,19 @@ PML::FillBoundaryB (PatchType patch_type) void PML::FillBoundaryF () { - if (pml_F_fp && pml_F_fp->nGrowVect().max() > 0) + FillBoundaryF(PatchType::fine); + FillBoundaryF(PatchType::coarse); +} + +void +PML::FillBoundaryF (PatchType patch_type) +{ + if (patch_type == PatchType::fine && pml_F_fp && pml_F_fp->nGrowVect().max() > 0) { const auto& period = m_geom->periodicity(); pml_F_fp->FillBoundary(period); } - - if (pml_F_cp && pml_F_cp->nGrowVect().max() > 0) + else if (patch_type == PatchType::coarse && pml_F_cp && pml_F_cp->nGrowVect().max() > 0) { const auto& period = m_cgeom->periodicity(); pml_F_cp->FillBoundary(period); From 0daf6cc76b230fe90351b49a4dbe713d284639a0 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Wed, 26 Sep 2018 09:40:20 -0700 Subject: [PATCH 33/43] add FillBoundaryF to subcycling --- Source/WarpXEvolve.cpp | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index 46504dd2c54..ca443a70ba0 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -283,6 +283,7 @@ WarpX::OneStep_sub1 (Real curtime) EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]); EvolveF(fine_lev, PatchType::fine, 0.5*dt[fine_lev], DtType::FirstHalf); FillBoundaryB(fine_lev, PatchType::fine); + FillBoundaryF(fine_lev, PatchType::fine); EvolveE(fine_lev, PatchType::fine, dt[fine_lev]); FillBoundaryE(fine_lev, PatchType::fine); @@ -306,6 +307,7 @@ WarpX::OneStep_sub1 (Real curtime) EvolveB(fine_lev, PatchType::coarse, dt[fine_lev]); EvolveF(fine_lev, PatchType::coarse, dt[fine_lev], DtType::FirstHalf); FillBoundaryB(fine_lev, PatchType::coarse); + FillBoundaryF(fine_lev, PatchType::coarse); EvolveE(fine_lev, PatchType::coarse, dt[fine_lev]); FillBoundaryE(fine_lev, PatchType::coarse); @@ -313,7 +315,8 @@ WarpX::OneStep_sub1 (Real curtime) EvolveB(coarse_lev, PatchType::fine, 0.5*dt[coarse_lev]); EvolveF(coarse_lev, PatchType::fine, 0.5*dt[coarse_lev], DtType::FirstHalf); FillBoundaryB(coarse_lev, PatchType::fine); - + FillBoundaryF(coarse_lev, PatchType::fine); + EvolveE(coarse_lev, PatchType::fine, 0.5*dt[coarse_lev]); FillBoundaryE(coarse_lev, PatchType::fine); @@ -332,6 +335,7 @@ WarpX::OneStep_sub1 (Real curtime) EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]); EvolveF(fine_lev, PatchType::fine, 0.5*dt[fine_lev], DtType::FirstHalf); FillBoundaryB(fine_lev, PatchType::fine); + FillBoundaryF(fine_lev, PatchType::fine); EvolveE(fine_lev, PatchType::fine, dt[fine_lev]); FillBoundaryE(fine_lev, PatchType::fine); @@ -345,6 +349,7 @@ WarpX::OneStep_sub1 (Real curtime) } FillBoundaryB(fine_lev, PatchType::fine); + FillBoundaryF(fine_lev, PatchType::fine); // v) RestoreCurrent(coarse_lev); @@ -364,6 +369,7 @@ WarpX::OneStep_sub1 (Real curtime) } FillBoundaryB(fine_lev, PatchType::coarse); + FillBoundaryF(fine_lev, PatchType::coarse); EvolveE(coarse_lev, PatchType::fine, 0.5*dt[coarse_lev]); FillBoundaryE(coarse_lev, PatchType::fine); From ab2dc42b3b039680d81873c2863ded064e54c482 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Wed, 26 Sep 2018 11:38:49 -0700 Subject: [PATCH 34/43] need to make sure there are at least two ghost cells when there is moving window --- Source/WarpX.cpp | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index a2bb7b54bb5..a8f7c3ecc7f 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -526,6 +526,15 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d int ngJy = ngy_tmp; int ngJz = ngz_tmp; + if (do_moving_window) { + ngx = std::max(ngx,2); + ngy = std::max(ngy,2); + ngz = std::max(ngz,2); + ngJx = std::max(ngJx,2); + ngJy = std::max(ngJy,2); + ngJz = std::max(ngJz,2); + } + #if (AMREX_SPACEDIM == 3) IntVect ngE(ngx,ngy,ngz); IntVect ngJ(ngJx,ngJy,ngJz); From 80090d493e667e22855fd6804d1c15fc0b00eed0 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Wed, 26 Sep 2018 17:14:14 -0700 Subject: [PATCH 35/43] use fines_level instead of max_level in moving window --- Source/WarpXEvolve.cpp | 1 + Source/WarpXMove.cpp | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index ca443a70ba0..71407ee66eb 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -59,6 +59,7 @@ WarpX::EvolveEM (int numsteps) #ifdef WARPX_USE_PSATD amrex::Abort("LoadBalance for PSATD: TODO"); #endif + if (step > 0 && (step+1) % load_balance_int == 0) { LoadBalance(); diff --git a/Source/WarpXMove.cpp b/Source/WarpXMove.cpp index d8d26c5ff74..17fa6db2b20 100644 --- a/Source/WarpXMove.cpp +++ b/Source/WarpXMove.cpp @@ -54,13 +54,13 @@ WarpX::MoveWindow (bool move_j) new_lo[dir] = current_lo[dir] + num_shift_base * dx[dir]; new_hi[dir] = current_hi[dir] + num_shift_base * dx[dir]; RealBox new_box(new_lo, new_hi); - geom[0].ProbDomain(new_box); + Geometry::ProbDomain(new_box); int num_shift = num_shift_base; int num_shift_crse = num_shift; // Shift the mesh fields - for (int lev = 0; lev <= max_level; ++lev) { + for (int lev = 0; lev <= finest_level; ++lev) { if (lev > 0) { num_shift_crse = num_shift; From 026e892bde22645bba96688a088bf30e3a936c40 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Wed, 26 Sep 2018 18:15:59 -0700 Subject: [PATCH 36/43] fix bug in load balance --- Source/WarpXEvolve.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index 71407ee66eb..ce3f5617a08 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -432,7 +432,7 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real dt) } MultiFab* cost = costs[lev].get(); - const IntVect& rr = (lev < finestLevel()) ? refRatio(lev) : IntVect::TheUnitVector(); + const IntVect& rr = (lev > 0) ? refRatio(lev-1) : IntVect::TheUnitVector(); // Loop through the grids, and over the tiles within each grid #ifdef _OPENMP @@ -557,7 +557,7 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real dt) } MultiFab* cost = costs[lev].get(); - const IntVect& rr = (lev < finestLevel()) ? refRatio(lev) : IntVect::TheUnitVector(); + const IntVect& rr = (lev > 0) ? refRatio(lev-1) : IntVect::TheUnitVector(); // Loop through the grids, and over the tiles within each grid #ifdef _OPENMP From da6e0bfa554c3155b5a1320718713126b3061d3a Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Wed, 26 Sep 2018 20:53:53 -0700 Subject: [PATCH 37/43] recreate current_store in regrid --- Source/WarpXRegrid.cpp | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/Source/WarpXRegrid.cpp b/Source/WarpXRegrid.cpp index 744dd18e8e1..d21e2ad9c39 100644 --- a/Source/WarpXRegrid.cpp +++ b/Source/WarpXRegrid.cpp @@ -63,6 +63,14 @@ WarpX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa // pmf->Redistribute(*current_fp[lev][idim], 0, 0, 1, ng); current_fp[lev][idim] = std::move(pmf); } + if (current_store[lev][idim]) + { + const IntVect& ng = current_store[lev][idim]->nGrowVect(); + auto pmf = std::unique_ptr(new MultiFab(current_store[lev][idim]->boxArray(), + dm, 1, ng)); + // no need to redistribute + current_store[lev][idim] = std::move(pmf); + } } if (F_fp[lev] != nullptr) { From 194e94ff2edfa206644868b7a5e1b4c866dec1e8 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Fri, 28 Sep 2018 09:34:20 -0700 Subject: [PATCH 38/43] more walltime info --- Source/WarpXEvolve.cpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index ce3f5617a08..a0fe641a35b 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -47,6 +47,7 @@ WarpX::EvolveEM (int numsteps) Real walltime, walltime_start = amrex::second(); for (int step = istep[0]; step < numsteps_max && cur_time < stop_time; ++step) { + Real walltime_beg_step = amrex::second(); // Start loop on time steps amrex::Print() << "\nSTEP " << step+1 << " starts ...\n"; @@ -146,9 +147,11 @@ WarpX::EvolveEM (int numsteps) amrex::Print()<< "STEP " << step+1 << " ends." << " TIME = " << cur_time << " DT = " << dt[0] << "\n"; - walltime = amrex::second() - walltime_start; + Real walltime_end_step = amrex::second(); + walltime = walltime_end_step - walltime_start; amrex::Print()<< "Walltime = " << walltime - << " s; Avg. per step = " << walltime/(step+1) << " s\n"; + << " s; This step = " << walltime_end_step-walltime_beg_step + << " s; Avg. per step = " << walltime/(step+1) << " s\n"; // sync up time for (int i = 0; i <= max_level; ++i) { From afc99d266555c40496d4ff9a033c7702c123d897 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 10 Oct 2018 17:02:14 -0700 Subject: [PATCH 39/43] Fix merge error --- Source/WarpX.H | 11 +---------- 1 file changed, 1 insertion(+), 10 deletions(-) diff --git a/Source/WarpX.H b/Source/WarpX.H index 03bb697357f..32df63e1916 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -236,7 +236,7 @@ private: /// Advance the simulation by numsteps steps, electromagnetic case. /// void EvolveEM(int numsteps); - + void FillBoundaryB (int lev, PatchType patch_type); void FillBoundaryE (int lev, PatchType patch_type); void FillBoundaryF (int lev, PatchType patch_type); @@ -261,15 +261,6 @@ private: void AddRhoFromFineLevelandSumBoundary (int lev, int icomp, int ncomp); void NodalSyncRho (int lev, PatchType patch_type, int icomp, int ncomp); - void FillBoundaryB (int lev, PatchType patch_type); - void FillBoundaryE (int lev, PatchType patch_type); - void FillBoundaryF (int lev, PatchType patch_type); - - void EvolveB (int lev, PatchType patch_type, amrex::Real dt); - void EvolveE (int lev, PatchType patch_type, amrex::Real dt); - void EvolveF (int lev, PatchType patch_type, amrex::Real dt, DtType dt_type); - void DampPML (int lev, PatchType patch_type); - #ifdef WARPX_DO_ELECTROSTATIC /// /// Advance the simulation by numsteps steps, electrostatic case. From 05d31ae6ccff4b9be787d6ab77526fce02f505b3 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Mon, 15 Oct 2018 10:35:23 -0700 Subject: [PATCH 40/43] Add script for subcycling --- Examples/Tests/subcycling/inputs.2d | 154 ++++++++++++++++++++++++++++ Regression/WarpX-tests.ini | 13 +++ 2 files changed, 167 insertions(+) create mode 100644 Examples/Tests/subcycling/inputs.2d diff --git a/Examples/Tests/subcycling/inputs.2d b/Examples/Tests/subcycling/inputs.2d new file mode 100644 index 00000000000..f1e2951c98c --- /dev/null +++ b/Examples/Tests/subcycling/inputs.2d @@ -0,0 +1,154 @@ +# stop_time = 1.12966840205e-10 +amrex.v=1 +max_step = 2500 +amr.n_cell = 64 256 + +amr.max_grid_size = 4096 +amr.blocking_factor = 16 + +# Maximum level in hierarchy (for now must be 0, i.e., one level in total) +amr.max_level = 1 +amr.plot_file = "plotfiles/plt" +amr.plot_int = 200 + +warpx.fine_tag_lo = -2.e-6 -15.e-6 +warpx.fine_tag_hi = 2.e-6 -7.e-6 + + +# Geometry +geometry.coord_sys = 0 # 0: Cartesian +geometry.is_periodic = 0 0 # Is periodic? +geometry.prob_lo = -30.e-6 -20.e-6 # physical domain +geometry.prob_hi = 30.e-6 0.e-6 + +# Verbosity +warpx.verbose = 1 +warpx.plot_raw_fields = 1 +warpx.do_dive_cleaning = 0 +warpx.use_filter = 1 +warpx.do_pml = 1 +warpx.do_subcycling = 0 +warpx.refine_plasma = 0 +warpx.plot_raw_fields = 1 +warpx.plot_raw_fields_guards = 1 +warpx.plot_finepatch = 1 +warpx.plot_crsepatch = 1 +warpx.n_current_deposition_buffer = 0 +warpx.n_field_gather_buffer = 0 + +# Algorithms +algo.current_deposition = 0 +algo.charge_deposition = 0 +algo.field_gathering = 0 +algo.particle_pusher = 0 +algo.maxwell_fdtd_solver = "ckc" + +# CFL +warpx.cfl = .9999 +particles.nspecies = 4 +particles.species_names = driver beam plasma_e plasma_p + +# interpolation +interpolation.nox = 3 +interpolation.noy = 3 +interpolation.noz = 3 + +# +# The driver species information +# + +driver.charge = -q_e +driver.mass = m_e +driver.injection_style = "gaussian_beam" +driver.x_rms = 3.e-6 +driver.y_rms = 3.e-6 +driver.z_rms = .2e-6 +driver.x_m = 0. +driver.y_m = 0. +driver.z_m = -3.e-6 +driver.npart = 10000 +driver.q_tot = -1.e-10 +driver.profile = "constant" +driver.density = 5.e24 # number of particles per m^3 +driver.momentum_distribution_type = "gaussian" +driver.ux_m = 0.0 +driver.uy_m = 0.0 +driver.uz_m = 1e6 +driver.uz_th = 0. +#driver.ux_th = 2. +#driver.uy_th = 2. +#driver.uz_th = 20000. +driver.zinject_plane = 0. +driver.rigid_advance = true +driver.projected = true +driver.focused = false + +# +# The beam species information +# + +beam.charge = -q_e +beam.mass = m_e +beam.injection_style = "gaussian_beam" +beam.x_rms = .1e-6 +beam.y_rms = .1e-6 +beam.z_rms = .2e-6 +beam.x_m = 0. +beam.y_m = 0. +beam.z_m = -11.e-6 +beam.npart = 10000 +beam.q_tot = -1.e-12 +beam.profile = "constant" +beam.density = 1.e24 +beam.momentum_distribution_type = "gaussian" +beam.ux_m = 0.0 +beam.uy_m = 0.0 +beam.uz_m = 20. +beam.u_th = 0. + +# +# The plasma species information +# + +plasma_e.charge = -q_e +plasma_e.mass = m_e +plasma_e.injection_style = "NUniformPerCell" +plasma_e.xmin = -15.e-6 +plasma_e.xmax = 15.e-6 +plasma_e.ymin = -15.e-6 +plasma_e.ymax = 15.e-6 +plasma_e.zmin = 0.e-6 +plasma_e.profile = "constant" +plasma_e.density = 1.e25 +plasma_e.num_particles_per_cell_each_dim = 2 2 +plasma_e.momentum_distribution_type = "constant" +plasma_e.ux = 0.0 +plasma_e.uy = 0.0 +plasma_e.uz = 0.0 + + +plasma_p.charge = q_e +plasma_p.mass = m_p +plasma_p.injection_style = "NUniformPerCell" +plasma_p.xmin = -15.e-6 +plasma_p.xmax = 15.e-6 +plasma_p.ymin = -15.e-6 +plasma_p.ymax = 15.e-6 +plasma_p.zmin = 0.e-6 +plasma_p.profile = "constant" +plasma_p.density = 1.e25 +plasma_p.num_particles_per_cell_each_dim = 2 2 +plasma_p.momentum_distribution_type = "constant" +plasma_p.ux = 0.0 +plasma_p.uy = 0.0 +plasma_p.uz = 0.0 + +# Moving window +warpx.do_moving_window = 1 +warpx.moving_window_dir = z +warpx.moving_window_v = 1. # in units of the speed of light + +# Particle Injection +warpx.do_plasma_injection = 1 +warpx.num_injected_species = 2 +warpx.injected_plasma_species = 2 3 diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index c6b8e378f17..6b42e3211af 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -263,6 +263,19 @@ doVis = 0 compareParticles = 1 particleTypes = electrons beam +[subcyclingMR] +buildDir = . +inputFile = Examples/Tests/subcycling/inputs.2d +dim = 2 +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 2 +compileTest = 0 +doVis = 0 +compareParticles = 0 + [LaserAccelerationMR] buildDir = . inputFile = Examples/Physics_applications/laser_acceleration/inputs.2d From abbee70a6c55a671e89d90c26b41a3d005fd63d7 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Tue, 16 Oct 2018 17:15:20 -0700 Subject: [PATCH 41/43] Add comments --- Source/WarpX.cpp | 17 ++++++++++++----- Source/WarpXComm.cpp | 35 +++++++++++++++++++++++++++++++++++ 2 files changed, 47 insertions(+), 5 deletions(-) diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 9d19732eb9b..1891bc7cb86 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -263,13 +263,13 @@ WarpX::ReadParameters () const std::string msg = "Unknown moving_window_dir: "+s; amrex::Abort(msg.c_str()); } - + moving_window_x = geom[0].ProbLo(moving_window_dir); - + pp.get("moving_window_v", moving_window_v); moving_window_v *= PhysConst::c; } - + pp.query("do_plasma_injection", do_plasma_injection); if (do_plasma_injection) { pp.get("num_injected_species", num_injected_species); @@ -284,10 +284,10 @@ WarpX::ReadParameters () current_injection_position = geom[0].ProbLo(moving_window_dir); } } - + pp.query("do_boosted_frame_diagnostic", do_boosted_frame_diagnostic); if (do_boosted_frame_diagnostic) { - + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(gamma_boost > 1.0, "gamma_boost must be > 1 to use the boosted frame diagnostic."); @@ -510,6 +510,9 @@ WarpX::ClearLevel (int lev) void WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& dm) { + // When using subcycling, the particles on the fine level perform two pushes + // before being redistributed ; therefore, we need one extra guard cell + // (the particles may move by 2*c*dt) const int ngx_tmp = (maxLevel() > 0 && do_subcycling == 1) ? WarpX::nox+1 : WarpX::nox; const int ngy_tmp = (maxLevel() > 0 && do_subcycling == 1) ? WarpX::noy+1 : WarpX::noy; const int ngz_tmp = (maxLevel() > 0 && do_subcycling == 1) ? WarpX::noz+1 : WarpX::noz; @@ -518,6 +521,8 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d // jx, jy, jz and rho have the same number of ghost cells. // E and B have the same number of ghost cells as j and rho if NCI filter is not used, // but different number of ghost cells in z-direction if NCI filter is used. + // The number of cells should be even, in order to easily perform the + // interpolation from fine grid to coarse grid. int ngx = (ngx_tmp % 2) ? ngx_tmp+1 : ngx_tmp; // Always even number int ngy = (ngy_tmp % 2) ? ngy_tmp+1 : ngy_tmp; // Always even number int ngz_nonci = (ngz_tmp % 2) ? ngz_tmp+1 : ngz_tmp; // Always even number @@ -533,6 +538,8 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d int ngJy = ngy_tmp; int ngJz = ngz_tmp; + // When calling the moving window (with one level of refinement), we shift + // the fine grid by 2 cells ; therefore, we need at least 2 guard cells if (do_moving_window) { ngx = std::max(ngx,2); ngy = std::max(ngy,2); diff --git a/Source/WarpXComm.cpp b/Source/WarpXComm.cpp index 3d7321913ca..f225f10153a 100644 --- a/Source/WarpXComm.cpp +++ b/Source/WarpXComm.cpp @@ -504,6 +504,9 @@ WarpX::SyncCurrent () } } +/** \brief Fills the values of the current on the coarse patch by + * averaging the values of the current of the fine patch (on the same level). + */ void WarpX::SyncCurrent (const std::array& fine, const std::array< amrex::MultiFab*,3>& crse, @@ -651,6 +654,9 @@ WarpX::SyncRho (amrex::Vector >& rhof, } } +/** \brief Fills the values of the charge density on the coarse patch by + * averaging the values of the charge density of the fine patch (on the same level). + */ void WarpX::SyncRho (const MultiFab& fine, MultiFab& crse, int ref_ratio) { @@ -679,6 +685,9 @@ WarpX::SyncRho (const MultiFab& fine, MultiFab& crse, int ref_ratio) } } +/** \brief Fills the values of the current on the coarse patch by + * averaging the values of the current of the fine patch (on the same level). + */ void WarpX::RestrictCurrentFromFineToCoarsePatch (int lev) { @@ -717,6 +726,19 @@ WarpX::ApplyFilterandSumBoundaryJ (int lev, PatchType patch_type) } } +/* /brief Update the currents of `lev` by adding the currents from particles +* that are in the mesh refinement patches at `lev+1` +* +* More precisely, apply filter and sum boundaries for the current of: +* - the fine patch of `lev` +* - the coarse patch of `lev+1` (same resolution) +* - the buffer regions of the coarse patch of `lev+1` (i.e. for particules +* that are within the mesh refinement patch, but do not deposit on the +* mesh refinement patch because they are too close to the boundary) +* +* Then update the fine patch of `lev` by adding the currents for the coarse +* patch (and buffer region) of `lev+1` +*/ void WarpX::AddCurrentFromFineLevelandSumBoundary (int lev) { @@ -814,6 +836,19 @@ WarpX::ApplyFilterandSumBoundaryRho (int lev, PatchType patch_type, int icomp, i } } +/* /brief Update the charge density of `lev` by adding the charge density from particles +* that are in the mesh refinement patches at `lev+1` +* +* More precisely, apply filter and sum boundaries for the charge density of: +* - the fine patch of `lev` +* - the coarse patch of `lev+1` (same resolution) +* - the buffer regions of the coarse patch of `lev+1` (i.e. for particules +* that are within the mesh refinement patch, but do not deposit on the +* mesh refinement patch because they are too close to the boundary) +* +* Then update the fine patch of `lev` by adding the charge density for the coarse +* patch (and buffer region) of `lev+1` +*/ void WarpX::AddRhoFromFineLevelandSumBoundary(int lev, int icomp, int ncomp) { From 8cd1ff6e08ccf0a93457f7bbeea592fc21388e57 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 17 Oct 2018 14:45:37 -0700 Subject: [PATCH 42/43] Add more comments --- Source/WarpX.cpp | 8 +++++++- Source/WarpXEvolve.cpp | 32 +++++++++++++++++++++++++++----- 2 files changed, 34 insertions(+), 6 deletions(-) diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 1891bc7cb86..4da208aac85 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -522,7 +522,7 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d // E and B have the same number of ghost cells as j and rho if NCI filter is not used, // but different number of ghost cells in z-direction if NCI filter is used. // The number of cells should be even, in order to easily perform the - // interpolation from fine grid to coarse grid. + // interpolation from coarse grid to fine grid. int ngx = (ngx_tmp % 2) ? ngx_tmp+1 : ngx_tmp; // Always even number int ngy = (ngy_tmp % 2) ? ngy_tmp+1 : ngy_tmp; // Always even number int ngz_nonci = (ngz_tmp % 2) ? ngz_tmp+1 : ngz_tmp; // Always even number @@ -534,6 +534,8 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d ngz = ngz_nonci; } + // J is only interpolated from fine to coarse (not coarse to fine) + // and therefore does not need to be even. int ngJx = ngx_tmp; int ngJy = ngy_tmp; int ngJz = ngz_tmp; @@ -687,6 +689,8 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d Efield_cax[lev][2].reset( new MultiFab(amrex::convert(cba,Ez_nodal_flag),dm,1,ngE)); gather_buffer_masks[lev].reset( new iMultiFab(ba, dm, 1, 1) ); + // Gather buffer masks have 1 ghost cell, because of the fact + // that particles may move by more than one cell when using subcycling. } if (n_current_deposition_buffer > 0) { @@ -697,6 +701,8 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d charge_buf[lev].reset( new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,2,ngRho)); } current_buffer_masks[lev].reset( new iMultiFab(ba, dm, 1, 1) ); + // Current buffer masks have 1 ghost cell, because of the fact + // that particles may move by more than one cell when using subcycling. } } diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index e3af72cf4f3..69302929731 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -222,6 +222,10 @@ WarpX::EvolveEM (int numsteps) } } +/* /brief Perform one PIC iteration, without subcycling +* i.e. all levels/patches use the same timestep (that of the finest level) +* for the field advance and particle pusher. +*/ void WarpX::OneStep_nosub (Real cur_time) { @@ -266,6 +270,21 @@ WarpX::OneStep_nosub (Real cur_time) #endif } +/* /brief Perform one PIC iteration, with subcycling +* i.e. The fine patch uses a smaller timestep (and steps more often) +* than the coarse patch, for the field advance and particle pusher. +* +* This version of subcycling only works for 2 levels and with a refinement +* ratio of 2. +* The particles and fields of the fine patch are pushed twice +* (with dt[coarse]/2) in this routine. +* The particles of the coarse patch and mother grid are pushed only once +* (with dt[coarse]). The fields on the coarse patch and mother grid +* are pushed in a way which is equivalent to pushing once only, with +* a current which is the average of the coarse + fine current at the 2 +* steps of the fine grid. +* +*/ void WarpX::OneStep_sub1 (Real curtime) { @@ -275,7 +294,7 @@ WarpX::OneStep_sub1 (Real curtime) const int fine_lev = 1; const int coarse_lev = 0; - // i) + // i) Push particles and fields on the fine patch (first fine step) PushParticlesandDepose(fine_lev, curtime); RestrictCurrentFromFineToCoarsePatch(fine_lev); RestrictRhoFromFineToCoarsePatch(fine_lev); @@ -302,7 +321,9 @@ WarpX::OneStep_sub1 (Real curtime) FillBoundaryB(fine_lev, PatchType::fine); - // ii) + // ii) Push particles on the coarse patch and mother grid. + // Push the fields on the coarse patch and mother grid + // by only half a coarse step (first half) PushParticlesandDepose(coarse_lev, curtime); StoreCurrent(coarse_lev); AddCurrentFromFineLevelandSumBoundary(coarse_lev); @@ -324,10 +345,10 @@ WarpX::OneStep_sub1 (Real curtime) EvolveE(coarse_lev, PatchType::fine, 0.5*dt[coarse_lev]); FillBoundaryE(coarse_lev, PatchType::fine); - // iii) + // iii) Get auxiliary fields on the fine grid, at dt[fine_lev] UpdateAuxilaryData(); - // iv) + // iv) Push particles and fields on the fine patch (second fine step) PushParticlesandDepose(fine_lev, curtime+dt[fine_lev]); RestrictCurrentFromFineToCoarsePatch(fine_lev); RestrictRhoFromFineToCoarsePatch(fine_lev); @@ -355,7 +376,8 @@ WarpX::OneStep_sub1 (Real curtime) FillBoundaryB(fine_lev, PatchType::fine); FillBoundaryF(fine_lev, PatchType::fine); - // v) + // v) Push the fields on the coarse patch and mother grid + // by only half a coarse step (second half) RestoreCurrent(coarse_lev); AddCurrentFromFineLevelandSumBoundary(coarse_lev); AddRhoFromFineLevelandSumBoundary(coarse_lev, 1, 1); From c3d396d52db7c0438db3c06c716e0a5f9387b3c2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maxence=20Th=C3=A9venet?= Date: Wed, 17 Oct 2018 18:27:58 -0700 Subject: [PATCH 43/43] minor comment --- Source/WarpX.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 4da208aac85..ae39f073c52 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -542,6 +542,7 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d // When calling the moving window (with one level of refinement), we shift // the fine grid by 2 cells ; therefore, we need at least 2 guard cells + // on level 1. This may not be necessary for level 0. if (do_moving_window) { ngx = std::max(ngx,2); ngy = std::max(ngy,2);