diff options
Diffstat (limited to 'Source')
-rw-r--r-- | Source/PhysicalParticleContainer.H | 10 | ||||
-rw-r--r-- | Source/PhysicalParticleContainer.cpp | 108 |
2 files changed, 37 insertions, 81 deletions
diff --git a/Source/PhysicalParticleContainer.H b/Source/PhysicalParticleContainer.H index 98d15e22c..40a8a6068 100644 --- a/Source/PhysicalParticleContainer.H +++ b/Source/PhysicalParticleContainer.H @@ -42,14 +42,14 @@ public: virtual void PostRestart () override {} - // Inject 'ppc' particles per cell in Box 'part_box' - void AddParticles (int lev, const amrex::Box& part_box); + // Inject particles in Box 'part_box' + void AddParticles (int lev, amrex::Box part_box = amrex::Box()); protected: - void InitNDiagPerCell (); - void InitNRandomUniformPerCell (); - void InitNRandomNormal (); + void AddNDiagPerCell (int lev, amrex::Box part_box = amrex::Box()); + void AddNRandomUniformPerCell (int lev, amrex::Box part_box = amrex::Box()); + void AddNRandomNormal (int lev, amrex::Box part_box = amrex::Box()); std::string species_name; std::unique_ptr<PlasmaInjector> plasma_injector; diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp index 21fd95ba7..bc42b36a7 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/PhysicalParticleContainer.cpp @@ -16,12 +16,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp } void PhysicalParticleContainer::InitData() { - if (plasma_injector->injection_style == "ndiagpercell") - InitNDiagPerCell(); - else if (plasma_injector->injection_style == "nrandomuniformpercell") - InitNRandomUniformPerCell(); - else if (plasma_injector->injection_style == "nrandomnormal") - InitNRandomNormal(); + AddParticles(0); // Note - only one level right now } void @@ -34,15 +29,16 @@ PhysicalParticleContainer::AllocData () } void -PhysicalParticleContainer::InitNRandomUniformPerCell () { - BL_PROFILE("PhysicalParticleContainer::InitNRandomPerCell()"); +PhysicalParticleContainer::AddNRandomUniformPerCell (int lev, Box part_box) { + BL_PROFILE("PhysicalParticleContainer::AddNRandomPerCell()"); charge = plasma_injector->getCharge(); mass = plasma_injector->getMass(); - const int lev = 0; const Geometry& geom = Geom(lev); const Real* dx = geom.CellSize(); + + if (!part_box.ok()) part_box = geom.Domain(); Real scale_fac; int n_part_per_cell = plasma_injector->numParticlesPerCell(); @@ -62,13 +58,16 @@ PhysicalParticleContainer::InitNRandomUniformPerCell () { for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) { const Box& tile_box = mfi.tilebox(); - RealBox tile_real_box { tile_box, dx, geom.ProbLo() }; + const Box& intersectBox = tile_box & part_box; + if (!intersectBox.ok()) continue; + + RealBox tile_real_box { intersectBox, dx, geom.ProbLo() }; const int grid_id = mfi.index(); const int tile_id = mfi.LocalTileIndex(); - const auto& boxlo = tile_box.smallEnd(); - for (IntVect iv = tile_box.smallEnd(); iv <= tile_box.bigEnd(); tile_box.next(iv)) + const auto& boxlo = intersectBox.smallEnd(); + for (IntVect iv = intersectBox.smallEnd(); iv <= intersectBox.bigEnd(); intersectBox.next(iv)) { for (int i_part=0; i_part<n_part_per_cell;i_part++) { @@ -104,22 +103,23 @@ PhysicalParticleContainer::InitNRandomUniformPerCell () { } void -PhysicalParticleContainer::InitNRandomNormal () { - BL_PROFILE("PhysicalParticleContainer::InitNRandomNormal()"); - amrex::Abort("PhysicalParticleContainer::InitNRandomNormal() not implemented yet."); +PhysicalParticleContainer::AddNRandomNormal (int lev, Box part_box) { + BL_PROFILE("PhysicalParticleContainer::AddNRandomNormal()"); + amrex::Abort("PhysicalParticleContainer::AddNRandomNormal() not implemented yet."); } void -PhysicalParticleContainer::InitNDiagPerCell () { - BL_PROFILE("PhysicalParticleContainer::InitNDiagPerCell()"); +PhysicalParticleContainer::AddNDiagPerCell (int lev, Box part_box) { + BL_PROFILE("PhysicalParticleContainer::AddNDiagPerCell()"); charge = plasma_injector->getCharge(); mass = plasma_injector->getMass(); - const int lev = 0; const Geometry& geom = Geom(lev); const Real* dx = geom.CellSize(); + if (!part_box.ok()) part_box = geom.Domain(); + Real scale_fac; int n_part_per_cell = plasma_injector->numParticlesPerCell(); @@ -133,13 +133,16 @@ PhysicalParticleContainer::InitNDiagPerCell () { attribs.fill(0.0); for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) { const Box& tile_box = mfi.tilebox(); - RealBox tile_real_box { tile_box, dx, geom.ProbLo() }; + const Box& intersectBox = tile_box & part_box; + if (!intersectBox.ok()) continue; + + RealBox tile_real_box { intersectBox, dx, geom.ProbLo() }; const int grid_id = mfi.index(); const int tile_id = mfi.LocalTileIndex(); - const auto& boxlo = tile_box.smallEnd(); - for (IntVect iv = tile_box.smallEnd(); iv <= tile_box.bigEnd(); tile_box.next(iv)) + const auto& boxlo = intersectBox.smallEnd(); + for (IntVect iv = intersectBox.smallEnd(); iv <= intersectBox.bigEnd(); intersectBox.next(iv)) { for (int i_part=0; i_part<n_part_per_cell;i_part++) { @@ -153,7 +156,7 @@ PhysicalParticleContainer::InitNDiagPerCell () { Real y = 0.0; Real z = tile_real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; #endif - + if (plasma_injector->insideBounds(x, y, z)) { Real weight; std::array<Real, 3> u; @@ -279,61 +282,14 @@ PhysicalParticleContainer::FieldGather (int lev, } void -PhysicalParticleContainer::AddParticles (int lev, const Box& part_box) +PhysicalParticleContainer::AddParticles (int lev, Box part_box) { - const Geometry& geom = Geom(lev); - const Real* dx = geom.CellSize(); - - Real scale_fac; - int ppc = plasma_injector->numParticlesPerCell(); - -#if BL_SPACEDIM==3 - scale_fac = dx[0]*dx[1]*dx[2]/ppc; -#elif BL_SPACEDIM==2 - scale_fac = dx[0]*dx[1]/ppc; -#endif - - std::array<Real,PIdx::nattribs> 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()) - { - RealBox real_box { intersectBox, dx, geom.ProbLo() }; - - const int grid_id = mfi.index(); - const int tile_id = mfi.LocalTileIndex(); - - const IntVect& boxlo = tile_box.smallEnd(); - for (IntVect iv = boxlo; iv <= tile_box.bigEnd(); tile_box.next(iv)) - { - for (int i_part=0; i_part < ppc; i_part++) - { - Real particle_shift = (0.5+i_part)/ppc; -#if (BL_SPACEDIM == 3) - Real x = real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0]; - Real y = real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; - Real z = real_box.lo(2) + (iv[2]-boxlo[2] + particle_shift)*dx[2]; -#elif (BL_SPACEDIM == 2) - Real x = real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0]; - Real y = 0.0; - Real z = real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; -#endif - Real weight; - std::array<Real, 3> u; - weight = plasma_injector->getDensity(x, y, z) * scale_fac; - plasma_injector->getMomentum(u); - attribs[PIdx::w ] = weight; - attribs[PIdx::ux] = u[0]; - attribs[PIdx::uy] = u[1]; - attribs[PIdx::uz] = u[2]; - AddOneParticle(lev, grid_id, tile_id, x, y, z, attribs); - } - } - } - } + if (plasma_injector->injection_style == "ndiagpercell") + AddNDiagPerCell(lev, part_box); + else if (plasma_injector->injection_style == "nrandomuniformpercell") + AddNRandomUniformPerCell(lev, part_box); + else if (plasma_injector->injection_style == "nrandomnormal") + AddNRandomNormal(lev, part_box); } void |