From e9d166c0ded8a6186ea4bc43069b52e819597fa1 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Tue, 24 Oct 2017 17:18:06 -0700 Subject: Started implementation of boosted injection --- Source/PhysicalParticleContainer.cpp | 149 +++++++++++++++++++---------------- 1 file changed, 83 insertions(+), 66 deletions(-) (limited to 'Source/PhysicalParticleContainer.cpp') diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp index 53a465820..9a616bad6 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/PhysicalParticleContainer.cpp @@ -26,9 +26,9 @@ void PhysicalParticleContainer::InitData() } } -void +void PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, - Real x_rms, Real y_rms, Real z_rms, + Real x_rms, Real y_rms, Real z_rms, Real q_tot, long npart) { const Geometry& geom = m_gdb->Geom(0); @@ -76,7 +76,7 @@ PhysicalParticleContainer::AddParticles (int lev, Box part_box) BL_PROFILE("PhysicalParticleContainer::AddParticles()"); if (plasma_injector->add_single_particle) { - AddNParticles(lev, 1, + AddNParticles(lev, 1, &(plasma_injector->single_particle_pos[0]), &(plasma_injector->single_particle_pos[1]), &(plasma_injector->single_particle_pos[2]), @@ -91,12 +91,12 @@ PhysicalParticleContainer::AddParticles (int lev, Box part_box) AddGaussianBeam(plasma_injector->x_m, plasma_injector->y_m, plasma_injector->z_m, - plasma_injector->x_rms, + plasma_injector->x_rms, plasma_injector->y_rms, - plasma_injector->z_rms, - plasma_injector->q_tot, + plasma_injector->z_rms, + plasma_injector->q_tot, plasma_injector->npart); - + return; } @@ -116,35 +116,35 @@ PhysicalParticleContainer::AddParticles (int lev, Box part_box) scale_fac = dx[0]*dx[2]/num_ppc; #endif -#ifdef _OPENMP +#ifdef _OPENMP // First touch all tiles in the map in serial for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) { const int grid_id = mfi.index(); - const int tile_id = mfi.LocalTileIndex(); + const int tile_id = mfi.LocalTileIndex(); GetParticles(lev)[std::make_pair(grid_id, tile_id)]; } #endif - + #ifdef _OPENMP #pragma omp parallel if (not WarpX::serialize_ics) #endif - { + { std::array attribs; attribs.fill(0.0); - + for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) { const Box& tile_box = mfi.tilebox(); const Box& intersectBox = tile_box & part_box; if (!intersectBox.ok()) continue; - - const std::array& tile_corner = + + const std::array& tile_corner = WarpX::LowerCorner(intersectBox, lev); - + const int grid_id = mfi.index(); const int tile_id = mfi.LocalTileIndex(); - + const auto& boxlo = intersectBox.smallEnd(); - for (IntVect iv = intersectBox.smallEnd(); + for (IntVect iv = intersectBox.smallEnd(); iv <= intersectBox.bigEnd(); intersectBox.next(iv)) { for (int i_part=0; i_part r; @@ -159,11 +159,28 @@ PhysicalParticleContainer::AddParticles (int lev, Box part_box) Real z = tile_corner[2] + (iv[1]-boxlo[1] + r[2])*dx[2]; #endif if (plasma_injector->insideBounds(x, y, z)) { - Real weight; + Real dens; std::array u; - plasma_injector->getMomentum(u, x, y, z); - weight = plasma_injector->getDensity(x, y, z) * scale_fac; - attribs[PIdx::w ] = weight; + if (WarpX::gamma_boost == 1.){ + // Lab-frame simulation + plasma_injector->getMomentum(u, x, y, z); + dens = plasma_injector->getDensity(x, y, z); + } else { + // Boosted-frame simulation: call `getMomentum` + // and `getDensity` with lab-frame Parameters + // (Assumes that the plasma has a low velocity, + // and that the boost is along z) + Real t = WarpX::t_new[lev]; + Real v_boost = WarpX::beta_boost*PhysConst::c; + Real z_lab = WarpX::gamma_boost*( z - v_boost*t ); + plasma_injector->getMomentum(u, x, y, z_lab); + dens = plasma_injector->getDensity(x, y, z); + // Perform Lorentz transform + // (Assumes that the plasma has a low velocity) + u[2] = WarpX::gamma_boost * ( u[2] - v_boost ); + dens = WarpX::gamma_boost * dens; + } + attribs[PIdx::w ] = dens * scale_fac; attribs[PIdx::ux] = u[0]; attribs[PIdx::uy] = u[1]; attribs[PIdx::uz] = u[2]; @@ -223,11 +240,11 @@ FieldGatherES (const amrex::Vector, #endif WRPX_INTERPOLATE_CIC(particles.data(), nstride, np, - Exp.data(), Eyp.data(), -#if BL_SPACEDIM == 3 + Exp.data(), Eyp.data(), +#if BL_SPACEDIM == 3 Ezp.data(), #endif - exfab.dataPtr(), eyfab.dataPtr(), + exfab.dataPtr(), eyfab.dataPtr(), #if BL_SPACEDIM == 3 ezfab.dataPtr(), #endif @@ -247,7 +264,7 @@ FieldGatherES (const amrex::Vector, #if BL_SPACEDIM == 3 MultiFab coarse_Ez(coarsened_fine_BA, fine_dm, 1, 1); #endif - + coarse_Ex.copy(*E[0][0], 0, 0, 1, 1, 1); coarse_Ey.copy(*E[0][1], 0, 0, 1, 1, 1); #if BL_SPACEDIM == 3 @@ -293,35 +310,35 @@ FieldGatherES (const amrex::Vector, if (lev == 0) { WRPX_INTERPOLATE_CIC(particles.data(), nstride, np, - Exp.data(), Eyp.data(), -#if BL_SPACEDIM == 3 + Exp.data(), Eyp.data(), +#if BL_SPACEDIM == 3 Ezp.data(), #endif - exfab.dataPtr(), eyfab.dataPtr(), + exfab.dataPtr(), eyfab.dataPtr(), #if BL_SPACEDIM == 3 ezfab.dataPtr(), #endif - box.loVect(), box.hiVect(), plo, dx, &ng); + box.loVect(), box.hiVect(), plo, dx, &ng); } else { - + const FArrayBox& exfab_coarse = coarse_Ex[pti]; const FArrayBox& eyfab_coarse = coarse_Ey[pti]; #if BL_SPACEDIM == 3 const FArrayBox& ezfab_coarse = coarse_Ez[pti]; -#endif +#endif const Box& coarse_box = coarsened_fine_BA[pti]; const Real* coarse_dx = Geom(0).CellSize(); - + WRPX_INTERPOLATE_CIC_TWO_LEVELS(particles.data(), nstride, np, - Exp.data(), Eyp.data(), -#if BL_SPACEDIM == 3 + Exp.data(), Eyp.data(), +#if BL_SPACEDIM == 3 Ezp.data(), #endif - exfab.dataPtr(), eyfab.dataPtr(), + exfab.dataPtr(), eyfab.dataPtr(), #if BL_SPACEDIM == 3 ezfab.dataPtr(), #endif - box.loVect(), box.hiVect(), dx, + box.loVect(), box.hiVect(), dx, exfab_coarse.dataPtr(), eyfab_coarse.dataPtr(), #if BL_SPACEDIM == 3 ezfab_coarse.dataPtr(), @@ -439,9 +456,9 @@ PhysicalParticleContainer::EvolveES (const Vectorcharge, &this->mass, &dt, - prob_domain.lo(), prob_domain.hi()); + prob_domain.lo(), prob_domain.hi()); } } } @@ -488,23 +505,23 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_NS("PICSAR::ParticlePush", blp_pxr_pp); BL_PROFILE_VAR_NS("PICSAR::CurrentDeposition", blp_pxr_cd); BL_PROFILE_VAR_NS("PPC::Evolve::Accumulate", blp_accumulate); - + const std::array& dx = WarpX::CellSize(lev); - + // WarpX assumes the same number of guard cells for Ex, Ey, Ez, Bx, By, Bz long ngE = Ex.nGrow(); // WarpX assumes the same number of guard cells for Jx, Jy, Jz long ngJ = jx.nGrow(); - + long ngRho = (rho) ? rho->nGrow() : 0; - + long ngRhoDeposit = (WarpX::use_filter) ? ngRho +1 : ngRho; long ngJDeposit = (WarpX::use_filter) ? ngJ +1 : ngJ; BL_ASSERT(OnSameGrids(lev,Ex)); MultiFab* cost = WarpX::getCosts(lev); - + #ifdef _OPENMP #pragma omp parallel #endif @@ -512,7 +529,7 @@ PhysicalParticleContainer::Evolve (int lev, Vector xp, yp, zp, giv; FArrayBox local_rho, local_jx, local_jy, local_jz; FArrayBox filtered_rho, filtered_jx, filtered_jy, filtered_jz; - + for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { Real wt = ParallelDescriptor::second(); @@ -600,7 +617,7 @@ PhysicalParticleContainer::Evolve (int lev, &this->charge, &xyzmin[0], &xyzmin[1], &xyzmin[2], &dx[0], &dx[1], &dx[2], &nx, &ny, &nz, - &ngRhoDeposit, &ngRhoDeposit, &ngRhoDeposit, + &ngRhoDeposit, &ngRhoDeposit, &ngRhoDeposit, &WarpX::nox,&WarpX::noy,&WarpX::noz, &lvect, &WarpX::charge_deposition_algo); @@ -627,7 +644,7 @@ PhysicalParticleContainer::Evolve (int lev, BL_TO_FORTRAN_3D(rhofab), ncomp); } } - + if (! do_not_push) { // @@ -653,7 +670,7 @@ PhysicalParticleContainer::Evolve (int lev, &ll4symtry, &l_lower_order_in_v, &lvect_fieldgathe, &WarpX::field_gathering_algo); BL_PROFILE_VAR_STOP(blp_pxr_fg); - + // // Particle Push // @@ -665,7 +682,7 @@ PhysicalParticleContainer::Evolve (int lev, &this->charge, &this->mass, &dt, &WarpX::particle_pusher_algo); BL_PROFILE_VAR_STOP(blp_pxr_pp); - + // // Current Deposition onto fine patch // @@ -677,18 +694,18 @@ PhysicalParticleContainer::Evolve (int lev, Box tby = convert(pti.tilebox(), WarpX::jy_nodal_flag); Box tbz = convert(pti.tilebox(), WarpX::jz_nodal_flag); Box gtbx, gtby, gtbz; - + const std::array& xyzmin = xyzmin_tile; - + tbx.grow(ngJ); tby.grow(ngJ); tbz.grow(ngJ); - + if (WarpX::use_filter) { gtbx = tbx; gtbx.grow(1); - + gtby = tby; gtby.grow(1); @@ -707,11 +724,11 @@ PhysicalParticleContainer::Evolve (int lev, 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(); @@ -730,7 +747,7 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_STOP(blp_pxr_cd); - BL_PROFILE_VAR_START(blp_accumulate); + BL_PROFILE_VAR_START(blp_accumulate); const int ncomp = 1; if (WarpX::use_filter) { @@ -767,29 +784,29 @@ PhysicalParticleContainer::Evolve (int lev, filtered_jz.loVect(), filtered_jz.hiVect(), ncomp); - + amrex_atomic_accumulate_fab(BL_TO_FORTRAN_3D(filtered_jx), BL_TO_FORTRAN_3D(jxfab), ncomp); - + amrex_atomic_accumulate_fab(BL_TO_FORTRAN_3D(filtered_jy), BL_TO_FORTRAN_3D(jyfab), ncomp); - + amrex_atomic_accumulate_fab(BL_TO_FORTRAN_3D(filtered_jz), BL_TO_FORTRAN_3D(jzfab), ncomp); - + } else { - + 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); - + // // copy particle data back // -- cgit v1.2.3 From 39736ee382cd8b124b34dfe6d7eb283352aecb3a Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Mon, 6 Nov 2017 15:38:06 -0800 Subject: Inject particles in integer number of cells --- Source/PhysicalParticleContainer.H | 8 ++--- Source/PhysicalParticleContainer.cpp | 60 +++++++++++++++++++++++++++--------- Source/WarpX.H | 5 ++- Source/WarpX.cpp | 14 ++++++--- Source/WarpXEvolve.cpp | 52 +++++++++++++------------------ Source/WarpXMove.cpp | 60 ++++++++++++++++++++++++++++++------ 6 files changed, 135 insertions(+), 64 deletions(-) (limited to 'Source/PhysicalParticleContainer.cpp') diff --git a/Source/PhysicalParticleContainer.H b/Source/PhysicalParticleContainer.H index 17a6f0774..e14ca5f21 100644 --- a/Source/PhysicalParticleContainer.H +++ b/Source/PhysicalParticleContainer.H @@ -29,7 +29,7 @@ public: const amrex::MultiFab& Bz) final; virtual void EvolveES (const amrex::Vector, 3> >& E, - amrex::Vector >& rho, + amrex::Vector >& rho, amrex::Real t, amrex::Real dt) override; virtual void Evolve (int lev, @@ -49,12 +49,12 @@ public: virtual void PostRestart () final {} // Inject particles in Box 'part_box' - void AddParticles (int lev, amrex::Box part_box = amrex::Box()); + void AddParticles (int lev, amrex::RealBox part_box = amrex::RealBox()); void AddGaussianBeam(amrex::Real x_m, amrex::Real y_m, amrex::Real z_m, - amrex::Real x_rms, amrex::Real y_rms, amrex::Real z_rms, + amrex::Real x_rms, amrex::Real y_rms, amrex::Real z_rms, amrex::Real q_tot, long npart); - + protected: std::string species_name; diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp index 9a616bad6..f95f26831 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/PhysicalParticleContainer.cpp @@ -71,7 +71,7 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, } void -PhysicalParticleContainer::AddParticles (int lev, Box part_box) +PhysicalParticleContainer::AddParticles (int lev, RealBox part_realbox) { BL_PROFILE("PhysicalParticleContainer::AddParticles()"); @@ -102,8 +102,9 @@ PhysicalParticleContainer::AddParticles (int lev, Box part_box) if ( not plasma_injector->doInjection() ) return; + // If no part_realbox is provided, initialize particles in the whole domain const Geometry& geom = Geom(lev); - if (!part_box.ok()) part_box = geom.Domain(); + if (!part_realbox.ok()) part_realbox = geom.ProbDomain(); int num_ppc = plasma_injector->num_particles_per_cell; @@ -132,32 +133,61 @@ PhysicalParticleContainer::AddParticles (int lev, Box part_box) std::array attribs; attribs.fill(0.0); + // Loop through the tiles for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) { - const Box& tile_box = mfi.tilebox(); - const Box& intersectBox = tile_box & part_box; - if (!intersectBox.ok()) continue; - const std::array& tile_corner = - WarpX::LowerCorner(intersectBox, lev); + const Box& tile_box = mfi.tilebox(); + const RealBox tile_realbox = WarpX::getRealBox(tile_box, lev); + + // Find the cells of part_box that overlap with tile_realbox + // If there is no overlap, just go to the next tile in the loop + RealBox overlap_realbox; + Box overlap_box; + bool no_overlap = 0; + for (int dir=0; dir part_realbox.lo(dir) ) { + overlap_realbox.setHi( dir, part_realbox.lo(dir) + + std::ceil( (tile_realbox.hi(dir) - part_realbox.lo(dir))/dx[dir]) * dx[dir]); + } else { + no_overlap = 1; break; + } + // Count the number of cells in this direction in overlap_realbox + overlap_box.setSmall( dir, 0 ); + overlap_box.setBig( dir, + int( (overlap_realbox.hi(dir)-overlap_realbox.lo(dir))/dx[dir] )); + } + if (no_overlap == 1) continue; // Go to the next tile const int grid_id = mfi.index(); const int tile_id = mfi.LocalTileIndex(); - const auto& boxlo = intersectBox.smallEnd(); - for (IntVect iv = intersectBox.smallEnd(); - iv <= intersectBox.bigEnd(); intersectBox.next(iv)) { + // Loop through the cells of overlap_realbox and inject + // the corresponding particles + const auto& overlap_corner = overlap_realbox.lo(); + for (IntVect iv = overlap_box.smallEnd(); + iv <= overlap_box.bigEnd(); overlap_box.next(iv)) { for (int i_part=0; i_part r; plasma_injector->getPositionUnitBox(r, i_part); #if ( BL_SPACEDIM == 3 ) - Real x = tile_corner[0] + (iv[0]-boxlo[0] + r[0])*dx[0]; - Real y = tile_corner[1] + (iv[1]-boxlo[1] + r[1])*dx[1]; - Real z = tile_corner[2] + (iv[2]-boxlo[2] + r[2])*dx[2]; + Real x = overlap_corner[0] + (iv[0] + r[0])*dx[0]; + Real y = overlap_corner[1] + (iv[1] + r[1])*dx[1]; + Real z = overlap_corner[2]+ (iv[2] + r[2])*dx[2]; #elif ( BL_SPACEDIM == 2 ) - Real x = tile_corner[0] + (iv[0]-boxlo[0] + r[0])*dx[0]; + Real x = overlap_corner[0] + (iv[0] + r[0])*dx[0]; Real y = 0.; - Real z = tile_corner[2] + (iv[1]-boxlo[1] + r[2])*dx[2]; + Real z = overlap_corner[2] + (iv[1] + r[2])*dx[2]; #endif + // If the new particle is not inside the tile box, + // go to the next generated particle. + if(!tile_realbox.contains({x, y, z})) continue; + if (plasma_injector->insideBounds(x, y, z)) { Real dens; std::array u; diff --git a/Source/WarpX.H b/Source/WarpX.H index 6c2b4f402..372fb9955 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -140,6 +140,7 @@ public: const amrex::Vector, 3> >& E); static std::array CellSize (int lev); + static amrex::RealBox getRealBox(const amrex::Box& bx, int lev); static std::array LowerCorner (const amrex::Box& bx, int lev); static std::array UpperCorner (const amrex::Box& bx, int lev); @@ -229,7 +230,7 @@ private: void InitPML (); void ComputePMLFactors (); - void InjectPlasma(int num_shift, int dir); + void InjectPlasma(int lev, amrex::RealBox particleBox); void WriteWarpXHeader(const std::string& name) const; void WriteJobInfo (const std::string& dir) const; @@ -325,6 +326,8 @@ private: int moving_window_dir = -1; amrex::Real moving_window_x = std::numeric_limits::max(); amrex::Real moving_window_v = std::numeric_limits::max(); + Real current_injection_position; + Real last_injection_time; // Plasma injection parameters int do_plasma_injection = 0; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index c856121d0..c4b89a3e4 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -423,11 +423,18 @@ WarpX::CellSize (int lev) #endif } -std::array -WarpX::LowerCorner(const Box& bx, int lev) +amrex::RealBox +WarpX::getRealBox(const Box& bx, int lev) { const auto& gm = GetInstance().Geom(lev); const RealBox grid_box{bx, gm.CellSize(), gm.ProbLo()}; + return( grid_box ); +} + +std::array +WarpX::LowerCorner(const Box& bx, int lev) +{ + const RealBox grid_box = getRealBox( bx, lev ); const Real* xyzmin = grid_box.lo(); #if (BL_SPACEDIM == 3) return { xyzmin[0], xyzmin[1], xyzmin[2] }; @@ -439,8 +446,7 @@ WarpX::LowerCorner(const Box& bx, int lev) std::array WarpX::UpperCorner(const Box& bx, int lev) { - const auto& gm = GetInstance().Geom(lev); - const RealBox grid_box{bx, gm.CellSize(), gm.ProbLo()}; + const RealBox grid_box = getRealBox( bx, lev ); const Real* xyzmax = grid_box.hi(); #if (BL_SPACEDIM == 3) return { xyzmax[0], xyzmax[1], xyzmax[2] }; diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index 86b7bc2ce..36657dc03 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -42,12 +42,12 @@ WarpX::EvolveES(int numsteps) { // nodal storage for the electrostatic case const int num_levels = max_level + 1; Vector > rhoNodal(num_levels); - Vector > phiNodal(num_levels); - Vector, 3> > eFieldNodal(num_levels); + Vector > phiNodal(num_levels); + Vector, 3> > eFieldNodal(num_levels); const int ng = 1; for (int lev = 0; lev <= max_level; lev++) { BoxArray nba = boxArray(lev); - nba.surroundingNodes(); + nba.surroundingNodes(); rhoNodal[lev].reset(new MultiFab(nba, dmap[lev], 1, ng)); phiNodal[lev].reset(new MultiFab(nba, dmap[lev], 1, 2)); @@ -59,12 +59,12 @@ WarpX::EvolveES(int numsteps) { const int lev = 0; for (int step = istep[0]; step < numsteps_max && cur_time < stop_time; ++step) { - + // Start loop on time steps amrex::Print() << "\nSTEP " << step+1 << " starts ...\n"; - - // At initialization, particles have p^{n-1/2} and x^{n-1/2}. - // Beyond one step, particles have p^{n-1/2} and x^{n}. + + // At initialization, particles have p^{n-1/2} and x^{n-1/2}. + // Beyond one step, particles have p^{n-1/2} and x^{n}. if (is_synchronized) { // on first step, push X by 0.5*dt mypc->PushXES(0.5*dt[lev]); @@ -88,17 +88,17 @@ WarpX::EvolveES(int numsteps) { computePhi(rhoNodal, phiNodal); computeE(eFieldNodal, phiNodal); - + if (cur_time + dt[0] >= stop_time - 1.e-3*dt[0] || step == numsteps_max-1) { // on last step, push by only 0.5*dt to synchronize all at n+1/2 mypc->PushXES(-0.5*dt[lev]); is_synchronized = true; - } + } mypc->Redistribute(); - + ++istep[0]; - + cur_time += dt[0]; bool to_make_plot = (plot_int > 0) && ((step+1) % plot_int == 0); @@ -152,7 +152,7 @@ WarpX::EvolveEM (int numsteps) Real cur_time = t_new[0]; static int last_plot_file_step = 0; static int last_check_file_step = 0; - + int numsteps_max; if (numsteps < 0) { // Note that the default argument is numsteps = -1 numsteps_max = max_step; @@ -231,7 +231,7 @@ WarpX::EvolveEM (int numsteps) } else { EvolveE(dt[0], DtType::Full); // We now have E^{n+1} } - + for (int lev = 0; lev <= max_level; ++lev) { ++istep[lev]; } @@ -289,13 +289,13 @@ WarpX::EvolveEM (int numsteps) FillBoundaryE(); FillBoundaryB(); UpdateAuxilaryData(); - + for (int lev = 0; lev <= finest_level; ++lev) { mypc->FieldGather(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]); } - + WritePlotFile(); } @@ -327,7 +327,7 @@ WarpX::EvolveB (int lev, Real dt, DtType typ) 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) { @@ -362,7 +362,7 @@ WarpX::EvolveB (int lev, Real dt, DtType typ) 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 WRPX_PXR_PUSH_BVEC( tbx.loVect(), tbx.hiVect(), @@ -405,7 +405,7 @@ WarpX::EvolveB (int lev, Real dt, DtType typ) 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(), @@ -479,7 +479,7 @@ WarpX::EvolveE (int lev, Real dt, DtType typ) 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 @@ -552,7 +552,7 @@ WarpX::EvolveE (int lev, Real dt, DtType typ) 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( tex.loVect(), tex.hiVect(), tey.loVect(), tey.hiVect(), @@ -697,19 +697,10 @@ WarpX::ComputeDt () } void -WarpX::InjectPlasma (int num_shift, int dir) +WarpX::InjectPlasma (int lev, RealBox particleBox) { if(do_plasma_injection) { - const int lev = 0; - - // particleBox encloses the cells where we generate particles - Box particleBox = geom[lev].Domain(); - int domainLength = particleBox.length(dir); - int sign = (num_shift < 0) ? -1 : 1; - particleBox.shift(dir, sign*(domainLength - std::abs(num_shift))); - particleBox &= geom[lev].Domain(); - for (int i = 0; i < num_injected_species; ++i) { int ispecies = injected_plasma_species[i]; WarpXParticleContainer& pc = mypc->GetParticleContainer(ispecies); @@ -718,4 +709,3 @@ WarpX::InjectPlasma (int num_shift, int dir) } } } - diff --git a/Source/WarpXMove.cpp b/Source/WarpXMove.cpp index 566205e5c..8c9fcdf71 100644 --- a/Source/WarpXMove.cpp +++ b/Source/WarpXMove.cpp @@ -1,5 +1,6 @@ #include +#include using namespace amrex; @@ -33,14 +34,15 @@ WarpX::MoveWindow (bool move_j) int num_shift = num_shift_base; int num_shift_crse = num_shift; + + // Shift the mesh fields for (int lev = 0; lev <= max_level; ++lev) { - + if (lev > 0) { num_shift_crse = num_shift; num_shift *= refRatio(lev-1)[dir]; } - // shift the mesh fields for (int dim = 0; dim < 3; ++dim) { shiftMF(*Bfield_fp[lev][dim], geom[lev], num_shift, dir); @@ -57,7 +59,7 @@ WarpX::MoveWindow (bool move_j) if (do_pml && pml[lev]->ok()) { const std::array& pml_B = pml[lev]->GetB_fp(); - const std::array& pml_E = pml[lev]->GetE_fp(); + const std::array& pml_E = pml[lev]->GetE_fp(); shiftMF(*pml_B[dim], geom[lev], num_shift, dir); shiftMF(*pml_E[dim], geom[lev], num_shift, dir); if (do_dive_cleaning) { @@ -84,7 +86,7 @@ WarpX::MoveWindow (bool move_j) if (do_pml && pml[lev]->ok()) { const std::array& pml_B = pml[lev]->GetB_cp(); - const std::array& pml_E = pml[lev]->GetE_cp(); + const std::array& pml_E = pml[lev]->GetE_cp(); shiftMF(*pml_B[dim], geom[lev-1], num_shift_crse, dir); shiftMF(*pml_E[dim], geom[lev-1], num_shift_crse, dir); if (do_dive_cleaning) { @@ -95,8 +97,48 @@ WarpX::MoveWindow (bool move_j) } } } - - InjectPlasma(num_shift_base, dir); + + // Continuously inject plasma in new cells (by default only on level 0) + if (WarpX::do_plasma_injection) { + const int lev = 0; + + if (WarpX::gamma_boost > 1){ + // In boosted-frame simulations, the plasma has moved since the last + // call to this function, and injection position needs to be updated + current_injection_position -= WarpX::beta_boost * + WarpX::boost_direction[dir] * PhysConst::c * dt[0]; + } + + // particleBox encloses the cells where we generate particles + // (only injects particles in an integer number of cells, + // for correct particle spacing) + RealBox particleBox = geom[lev].ProbDomain(); + Real new_injection_position; + if (moving_window_v > 0){ + // Forward-moving window + Real dx = geom[lev].CellSize(dir); + new_injection_position = current_injection_position + + std::floor( (geom[lev].ProbHi(dir) - current_injection_position)/dx ) * dx; + } else { + // Backward-moving window + Real dx = geom[lev].CellSize(dir); + new_injection_position = current_injection_position - + std::floor( (current_injection_position - geom[lev].ProbLo(dir))/dx) * dx; + } + // Modify the corresponding bounds of the particleBox + if (num_shift > 0) { + particleBox.setLo( dir, current_injection_position ); + particleBox.setHi( dir, new_injection_position ); + } else { + particleBox.setLo( dir, new_injection_position ); + particleBox.setHi( dir, current_injection_position ); + } + // Perform the injection of new particles in particleBox + if (particleBox.ok()) { + InjectPlasma( lev, particleBox); + current_injection_position = new_injection_position; + } + } } void @@ -106,13 +148,13 @@ WarpX::shiftMF(MultiFab& mf, const Geometry& geom, int num_shift, int dir) const DistributionMapping& dm = mf.DistributionMap(); const int nc = mf.nComp(); const int ng = mf.nGrow(); - + BL_ASSERT(ng >= num_shift); - + MultiFab tmpmf(ba, dm, nc, ng); MultiFab::Copy(tmpmf, mf, 0, 0, nc, ng); tmpmf.FillBoundary(geom.periodicity()); - + // Make a box that covers the region that the window moved into const IndexType& typ = ba.ixType(); const Box& domainBox = geom.Domain(); -- cgit v1.2.3 From ce96584c41470573f15bf5d29d20bdc6f52dddc1 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Mon, 6 Nov 2017 16:31:03 -0800 Subject: Modify access to time --- Source/PhysicalParticleContainer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'Source/PhysicalParticleContainer.cpp') diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp index f95f26831..e3fa62dc3 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/PhysicalParticleContainer.cpp @@ -200,7 +200,7 @@ PhysicalParticleContainer::AddParticles (int lev, RealBox part_realbox) // and `getDensity` with lab-frame Parameters // (Assumes that the plasma has a low velocity, // and that the boost is along z) - Real t = WarpX::t_new[lev]; + Real t = WarpX::GetInstance().gett_new(lev); Real v_boost = WarpX::beta_boost*PhysConst::c; Real z_lab = WarpX::gamma_boost*( z - v_boost*t ); plasma_injector->getMomentum(u, x, y, z_lab); -- cgit v1.2.3 From 7d88ced5ad8448a459e0b1753a378c9305c15448 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sat, 11 Nov 2017 15:16:08 -0800 Subject: Corrected injection so that the first cell is not duplicated --- Source/PhysicalParticleContainer.cpp | 9 +++++---- Source/WarpXMove.cpp | 6 +++--- 2 files changed, 8 insertions(+), 7 deletions(-) (limited to 'Source/PhysicalParticleContainer.cpp') diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp index e3fa62dc3..65d1b1381 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/PhysicalParticleContainer.cpp @@ -143,17 +143,18 @@ PhysicalParticleContainer::AddParticles (int lev, RealBox part_realbox) // If there is no overlap, just go to the next tile in the loop RealBox overlap_realbox; Box overlap_box; + Real ncells_adjust; bool no_overlap = 0; for (int dir=0; dir part_realbox.lo(dir) ) { - overlap_realbox.setHi( dir, part_realbox.lo(dir) + - std::ceil( (tile_realbox.hi(dir) - part_realbox.lo(dir))/dx[dir]) * dx[dir]); + ncells_adjust = std::floor( (part_realbox.hi(dir) - tile_realbox.hi(dir))/dx[dir] ); + overlap_realbox.setHi( dir, part_realbox.hi(dir) - std::max(ncells_adjust, 0.) * dx[dir]); } else { no_overlap = 1; break; } diff --git a/Source/WarpXMove.cpp b/Source/WarpXMove.cpp index 8c9fcdf71..3080ff442 100644 --- a/Source/WarpXMove.cpp +++ b/Source/WarpXMove.cpp @@ -114,7 +114,7 @@ WarpX::MoveWindow (bool move_j) // for correct particle spacing) RealBox particleBox = geom[lev].ProbDomain(); Real new_injection_position; - if (moving_window_v > 0){ + if (moving_window_v >= 0){ // Forward-moving window Real dx = geom[lev].CellSize(dir); new_injection_position = current_injection_position + @@ -126,7 +126,7 @@ WarpX::MoveWindow (bool move_j) std::floor( (current_injection_position - geom[lev].ProbLo(dir))/dx) * dx; } // Modify the corresponding bounds of the particleBox - if (num_shift > 0) { + if (moving_window_v >= 0) { particleBox.setLo( dir, current_injection_position ); particleBox.setHi( dir, new_injection_position ); } else { @@ -134,7 +134,7 @@ WarpX::MoveWindow (bool move_j) particleBox.setHi( dir, current_injection_position ); } // Perform the injection of new particles in particleBox - if (particleBox.ok()) { + if (particleBox.ok() and (current_injection_position != new_injection_position)){ InjectPlasma( lev, particleBox); current_injection_position = new_injection_position; } -- cgit v1.2.3 From ebf7a53e792841fb446e6a7bf96d81c26b1beaf6 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sat, 11 Nov 2017 15:55:13 -0800 Subject: Add documentation and suppress unused function --- Source/PhysicalParticleContainer.H | 3 ++- Source/PhysicalParticleContainer.cpp | 21 +++++++++++++++++++-- Source/WarpXEvolve.cpp | 14 -------------- Source/WarpXMove.cpp | 8 +++++++- 4 files changed, 28 insertions(+), 18 deletions(-) (limited to 'Source/PhysicalParticleContainer.cpp') diff --git a/Source/PhysicalParticleContainer.H b/Source/PhysicalParticleContainer.H index e14ca5f21..cb221c426 100644 --- a/Source/PhysicalParticleContainer.H +++ b/Source/PhysicalParticleContainer.H @@ -49,7 +49,8 @@ public: virtual void PostRestart () final {} // Inject particles in Box 'part_box' - void AddParticles (int lev, amrex::RealBox part_box = amrex::RealBox()); + void AddParticles (int lev); + void AddPlasma(int lev, amrex::RealBox part_realbox = amrex::RealBox()); void AddGaussianBeam(amrex::Real x_m, amrex::Real y_m, amrex::Real z_m, amrex::Real x_rms, amrex::Real y_rms, amrex::Real z_rms, diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp index 65d1b1381..871dc6a6b 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/PhysicalParticleContainer.cpp @@ -71,7 +71,7 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, } void -PhysicalParticleContainer::AddParticles (int lev, RealBox part_realbox) +PhysicalParticleContainer::AddParticles (int lev) { BL_PROFILE("PhysicalParticleContainer::AddParticles()"); @@ -100,7 +100,24 @@ PhysicalParticleContainer::AddParticles (int lev, RealBox part_realbox) return; } - if ( not plasma_injector->doInjection() ) return; + if ( plasma_injector->doInjection() ) { + AddPlasma( lev ); + } +} + +/** + * Create new macroparticles for this species, with a fixed + * number of particles per cell (in the cells of `part_realbox`). + * The new particles are only created inside the intersection of `part_realbox` + * with the local grid for the current proc. + * @param lev the index of the refinement level + * @param part_realbox the box in which new particles should be created + * (this box should correspond to an integer number of cells in each direction, + * but its boundaries need not be aligned with the actual cells of the simulation) + */ +void +PhysicalParticleContainer::AddPlasma(int lev, RealBox part_realbox ) +{ // If no part_realbox is provided, initialize particles in the whole domain const Geometry& geom = Geom(lev); diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index 36657dc03..5bd90bf7d 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -695,17 +695,3 @@ WarpX::ComputeDt () dt[0] = const_dt; } } - -void -WarpX::InjectPlasma (int lev, RealBox particleBox) -{ - if(do_plasma_injection) - { - for (int i = 0; i < num_injected_species; ++i) { - int ispecies = injected_plasma_species[i]; - WarpXParticleContainer& pc = mypc->GetParticleContainer(ispecies); - auto& ppc = dynamic_cast(pc); - ppc.AddParticles(lev, particleBox); - } - } -} diff --git a/Source/WarpXMove.cpp b/Source/WarpXMove.cpp index 3080ff442..2d7a4fa2a 100644 --- a/Source/WarpXMove.cpp +++ b/Source/WarpXMove.cpp @@ -135,7 +135,13 @@ WarpX::MoveWindow (bool move_j) } // Perform the injection of new particles in particleBox if (particleBox.ok() and (current_injection_position != new_injection_position)){ - InjectPlasma( lev, particleBox); + for (int i = 0; i < num_injected_species; ++i) { + int ispecies = injected_plasma_species[i]; + WarpXParticleContainer& pc = mypc->GetParticleContainer(ispecies); + auto& ppc = dynamic_cast(pc); + ppc.AddPlasma(lev, particleBox); + } + // Update the injection position current_injection_position = new_injection_position; } } -- cgit v1.2.3 From 0af58798392aa5250475414c9cfd7c882ce2879b Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 15 Nov 2017 08:24:59 -0800 Subject: Allow compilation for both 3D and 2D --- Source/PhysicalParticleContainer.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) (limited to 'Source/PhysicalParticleContainer.cpp') diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp index 871dc6a6b..848878f16 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/PhysicalParticleContainer.cpp @@ -204,7 +204,11 @@ PhysicalParticleContainer::AddPlasma(int lev, RealBox part_realbox ) #endif // If the new particle is not inside the tile box, // go to the next generated particle. - if(!tile_realbox.contains({x, y, z})) continue; +#if ( BL_SPACEDIM == 3 ) + if(!tile_realbox.contains( {x, y, z} )) continue; +#elif ( BL_SPACEDIM == 2 ) + if(!tile_realbox.contains( {x, z} )) continue; +#endif if (plasma_injector->insideBounds(x, y, z)) { Real dens; -- cgit v1.2.3 From 813b2fce7451c5169fe60c2cb20dba8ade6e53eb Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sun, 19 Nov 2017 17:20:56 -0800 Subject: Fix 2D injection --- Source/PhysicalParticleContainer.cpp | 17 +++++++++++------ Source/WarpXMove.cpp | 1 + 2 files changed, 12 insertions(+), 6 deletions(-) (limited to 'Source/PhysicalParticleContainer.cpp') diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp index 848878f16..95e04f345 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/PhysicalParticleContainer.cpp @@ -125,7 +125,7 @@ PhysicalParticleContainer::AddPlasma(int lev, RealBox part_realbox ) int num_ppc = plasma_injector->num_particles_per_cell; - const std::array& dx = WarpX::CellSize(lev); + const Real* dx = geom.CellSize(); Real scale_fac; #if BL_SPACEDIM==3 @@ -150,6 +150,10 @@ PhysicalParticleContainer::AddPlasma(int lev, RealBox part_realbox ) std::array attribs; attribs.fill(0.0); + std::cout << "0 " << part_realbox.lo(0) << " " << part_realbox.hi(0) << " " << dx[0] << std::endl; + std::cout << "1 " << part_realbox.lo(1) << " " << part_realbox.hi(1) << " " << dx[1] << std::endl; + std::cout << "2 " << part_realbox.lo(2) << " " << part_realbox.hi(2) << " " << dx[2] << std::endl; + // Loop through the tiles for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) { @@ -162,6 +166,7 @@ PhysicalParticleContainer::AddPlasma(int lev, RealBox part_realbox ) Box overlap_box; Real ncells_adjust; bool no_overlap = 0; + for (int dir=0; dirinsideBounds(x, y, z)) { diff --git a/Source/WarpXMove.cpp b/Source/WarpXMove.cpp index 4a638fe19..1fea06cc3 100644 --- a/Source/WarpXMove.cpp +++ b/Source/WarpXMove.cpp @@ -145,6 +145,7 @@ WarpX::MoveWindow (bool move_j) ppc.AddPlasma(lev, particleBox); } // Update the injection position + std::cout << "Moving window " << dir << " " << current_injection_position << " " << new_injection_position << std::endl; current_injection_position = new_injection_position; } } -- cgit v1.2.3 From 0bd9ab04c20a1d91f1d956b9109ad942a804631d Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Mon, 20 Nov 2017 09:56:27 -0800 Subject: Fix another bug in the injection --- Source/PhysicalParticleContainer.cpp | 6 +----- Source/WarpXMove.cpp | 1 - 2 files changed, 1 insertion(+), 6 deletions(-) (limited to 'Source/PhysicalParticleContainer.cpp') diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp index 95e04f345..9687b30b3 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/PhysicalParticleContainer.cpp @@ -131,7 +131,7 @@ PhysicalParticleContainer::AddPlasma(int lev, RealBox part_realbox ) #if BL_SPACEDIM==3 scale_fac = dx[0]*dx[1]*dx[2]/num_ppc; #elif BL_SPACEDIM==2 - scale_fac = dx[0]*dx[2]/num_ppc; + scale_fac = dx[0]*dx[1]/num_ppc; #endif #ifdef _OPENMP @@ -150,10 +150,6 @@ PhysicalParticleContainer::AddPlasma(int lev, RealBox part_realbox ) std::array attribs; attribs.fill(0.0); - std::cout << "0 " << part_realbox.lo(0) << " " << part_realbox.hi(0) << " " << dx[0] << std::endl; - std::cout << "1 " << part_realbox.lo(1) << " " << part_realbox.hi(1) << " " << dx[1] << std::endl; - std::cout << "2 " << part_realbox.lo(2) << " " << part_realbox.hi(2) << " " << dx[2] << std::endl; - // Loop through the tiles for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) { diff --git a/Source/WarpXMove.cpp b/Source/WarpXMove.cpp index 1fea06cc3..4a638fe19 100644 --- a/Source/WarpXMove.cpp +++ b/Source/WarpXMove.cpp @@ -145,7 +145,6 @@ WarpX::MoveWindow (bool move_j) ppc.AddPlasma(lev, particleBox); } // Update the injection position - std::cout << "Moving window " << dir << " " << current_injection_position << " " << new_injection_position << std::endl; current_injection_position = new_injection_position; } } -- cgit v1.2.3 From d41b7adebf9222272f6a390e82a49912bce070eb Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sun, 26 Nov 2017 21:52:30 -0800 Subject: Correct error in number of cells --- Source/PhysicalParticleContainer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'Source/PhysicalParticleContainer.cpp') diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp index 9687b30b3..49e9b7184 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/PhysicalParticleContainer.cpp @@ -179,7 +179,7 @@ PhysicalParticleContainer::AddPlasma(int lev, RealBox part_realbox ) // Count the number of cells in this direction in overlap_realbox overlap_box.setSmall( dir, 0 ); overlap_box.setBig( dir, - int( (overlap_realbox.hi(dir)-overlap_realbox.lo(dir))/dx[dir] )); + int( round((overlap_realbox.hi(dir)-overlap_realbox.lo(dir))/dx[dir] )) - 1); } if (no_overlap == 1) continue; // Go to the next tile -- cgit v1.2.3 From 289e66d7cbc4082bd51b745e0b9a4dc950930e68 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Mon, 27 Nov 2017 15:54:57 -0800 Subject: Take into account comments from @atmyers --- Source/PhysicalParticleContainer.cpp | 2 +- Source/WarpX.H | 2 -- 2 files changed, 1 insertion(+), 3 deletions(-) (limited to 'Source/PhysicalParticleContainer.cpp') diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp index 49e9b7184..e5f6b0c82 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/PhysicalParticleContainer.cpp @@ -186,7 +186,7 @@ PhysicalParticleContainer::AddPlasma(int lev, RealBox part_realbox ) const int grid_id = mfi.index(); const int tile_id = mfi.LocalTileIndex(); - // Loop through the cells of overlap_realbox and inject + // Loop through the cells of overlap_box and inject // the corresponding particles const auto& overlap_corner = overlap_realbox.lo(); for (IntVect iv = overlap_box.smallEnd(); diff --git a/Source/WarpX.H b/Source/WarpX.H index ab1cfef2f..c194e1038 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -236,8 +236,6 @@ private: void InitPML (); void ComputePMLFactors (); - void InjectPlasma(int lev, amrex::RealBox particleBox); - void InitDiagnostics (); void WriteWarpXHeader(const std::string& name) const; -- cgit v1.2.3