Skip to content

Commit

Permalink
Add fixed_ppc option to flux injection
Browse files Browse the repository at this point in the history
  • Loading branch information
@dpgrote
dpgrote committed Nov 9, 2023
1 parent a074014 commit 9d48854
Show file tree
Hide file tree
Showing 3 changed files with 66 additions and 13 deletions.
2 changes: 2 additions & 0 deletions Source/Initialization/PlasmaInjector.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -125,6 +125,8 @@ public:
// When compiled in cylindrical geometry, 0 = radial, 1 = azimuthal, 2 = z
int flux_normal_axis;
int flux_direction; // -1 for left, +1 for right
bool fixed_ppc_is_specified = false;
int fixed_ppc;

bool radially_weighted = true;

Expand Down
6 changes: 6 additions & 0 deletions Source/Initialization/PlasmaInjector.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -535,6 +535,12 @@ void PlasmaInjector::parseFlux (ParmParseWithOptionalGroup& pp_species_dot_sourc
utils::parser::makeParser(str_flux_function,{"x","y","z","t"}));
h_inj_flux.reset(new InjectorFlux((InjectorFluxParser*)nullptr,
flux_parser->compile<4>()));
} else if (flux_prof_s == "fixed_ppc") {
fixed_ppc_is_specified = true;
pp_species_dot_source.getWithParser("fixed_ppc", fixed_ppc);
pp_species_dot_source.getWithParser("flux", flux);
// Construct InjectorFlux with InjectorFluxConstant.
h_inj_flux.reset(new InjectorFlux((InjectorFluxConstant*)nullptr, flux));
} else {
SpeciesUtils::StringParseAbortMessage("Flux profile type", flux_prof_s);
}
Expand Down
71 changes: 58 additions & 13 deletions Source/Particles/PhysicalParticleContainer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1535,9 +1535,9 @@ PhysicalParticleContainer::AddPlasmaFlux (const PlasmaInjector * plasma_injector
Real scale_fac = 0._rt;
// Scale particle weight by the area of the emitting surface, within one cell
#if defined(WARPX_DIM_3D)
scale_fac = dx[0]*dx[1]*dx[2]/dx[plasma_injector->flux_normal_axis]/num_ppc_real;
scale_fac = dx[0]*dx[1]*dx[2]/dx[plasma_injector->flux_normal_axis]/num_ppc_real*dt;
#elif defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ)
scale_fac = dx[0]*dx[1]/num_ppc_real;
scale_fac = dx[0]*dx[1]/num_ppc_real*dt;
// When emission is in the r direction, the emitting surface is a cylinder.
// The factor 2*pi*r is added later below.
if (plasma_injector->flux_normal_axis == 0) scale_fac /= dx[0];
Expand All @@ -1547,7 +1547,7 @@ PhysicalParticleContainer::AddPlasmaFlux (const PlasmaInjector * plasma_injector
// When emission is in the theta direction (flux_normal_axis == 1),
// the emitting surface is a rectangle, within the plane of the simulation
#elif defined(WARPX_DIM_1D_Z)
scale_fac = dx[0]/num_ppc_real;
scale_fac = dx[0]/num_ppc_real*dt;
if (plasma_injector->flux_normal_axis == 2) scale_fac /= dx[0];
#endif

Expand Down Expand Up @@ -1595,15 +1595,15 @@ PhysicalParticleContainer::AddPlasmaFlux (const PlasmaInjector * plasma_injector
info.SetDynamic(true);
#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
#endif
for (MFIter mfi = MakeMFIter(0, info); mfi.isValid(); ++mfi)
for (WarpXParIter pti(*this, level_zero); pti.isValid(); ++pti)
{
if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
{
amrex::Gpu::synchronize();
}
auto wt = static_cast<amrex::Real>(amrex::second());

const Box& tile_box = mfi.tilebox();
const Box& tile_box = pti.tilebox();
const RealBox tile_realbox = WarpX::getRealBox(tile_box, 0);

// Find the cells of part_realbox that overlap with tile_realbox
Expand Down Expand Up @@ -1670,8 +1670,8 @@ PhysicalParticleContainer::AddPlasmaFlux (const PlasmaInjector * plasma_injector
continue; // Go to the next tile
}

const int grid_id = mfi.index();
const int tile_id = mfi.LocalTileIndex();
const int grid_id = pti.index();
const int tile_id = pti.LocalTileIndex();

const GpuArray<Real,AMREX_SPACEDIM> overlap_corner
{AMREX_D_DECL(overlap_realbox.lo(0),
Expand All @@ -1681,20 +1681,57 @@ PhysicalParticleContainer::AddPlasmaFlux (const PlasmaInjector * plasma_injector
// count the number of particles that each cell in overlap_box could add
Gpu::DeviceVector<int> counts(overlap_box.numPts(), 0);
Gpu::DeviceVector<int> offset(overlap_box.numPts());
auto pcounts = counts.data();
int * pcounts = counts.data();
const amrex::IntVect lrrfac = rrfac;
Box fine_overlap_box; // default Box is NOT ok().
if (refine_injection) {
fine_overlap_box = overlap_box & amrex::shift(fine_injection_box, -shifted);
}

bool fixed_ppc_is_specified = plasma_injector->fixed_ppc_is_specified;
int fixed_ppc = plasma_injector->fixed_ppc;
if (fixed_ppc_is_specified) {
// count the current number of particles in the cells
const auto np = pti.numParticles();
const auto GetPosition = GetParticlePosition(pti);
#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) || defined(WARPX_DIM_3D)
const amrex::Real xmin = overlap_corner[0];
const amrex::Real dxi = 1.0_rt/dx[0];
#endif
#if defined(WARPX_DIM_3D)
const amrex::Real ymin = overlap_corner[1];
const amrex::Real dyi = 1.0_rt/dx[1];
#endif
const amrex::Real zmin = overlap_corner[WARPX_ZINDEX];
const amrex::Real dzi = 1.0_rt/dx[WARPX_ZINDEX];
// Note that lbound(overlap_box) is always 0
amrex::ParallelFor(np,
[=] AMREX_GPU_HOST_DEVICE (long ip)
{
ParticleReal xp, yp, zp;
GetPosition(ip, xp, yp, zp);
#if defined(WARPX_DIM_3D)
const int i = static_cast<int>((xp - xmin)*dxi);
const int j = static_cast<int>((yp - ymin)*dyi);
const int k = static_cast<int>((zp - zmin)*dzi);
#elif defined(WARPX_DIM_1D_Z)
const int i = static_cast<int>((zp - zmin)*dzi);
#endif
const IntVect iv(AMREX_D_DECL(i, j, k));

if (overlap_box.contains(iv)) {
auto index = overlap_box.index(iv);
amrex::Gpu::Atomic::AddNoRet( pcounts+index, 1);
}
});
}

amrex::ParallelForRNG(overlap_box, [=] AMREX_GPU_DEVICE (int i, int j, int k, amrex::RandomEngine const& engine) noexcept
{
const IntVect iv(AMREX_D_DECL(i, j, k));
auto lo = getCellCoords(overlap_corner, dx, {0._rt, 0._rt, 0._rt}, iv);
auto hi = getCellCoords(overlap_corner, dx, {1._rt, 1._rt, 1._rt}, iv);

const int num_ppc_int = static_cast<int>(num_ppc_real + amrex::Random(engine));

if (flux_pos->overlapsWith(lo, hi))
{
auto index = overlap_box.index(iv);
Expand All @@ -1704,6 +1741,14 @@ PhysicalParticleContainer::AddPlasmaFlux (const PlasmaInjector * plasma_injector
} else {
r = 1;
}
int num_ppc_int;
if (fixed_ppc_is_specified) {
num_ppc_int = std::max(0, fixed_ppc - pcounts[index]);
amrex::Print() << "num_ppc_int " << num_ppc_int << " " << fixed_ppc << " " << pcounts[index] << std::endl;
} else {
num_ppc_int = static_cast<int>(num_ppc_real + amrex::Random(engine));
}

pcounts[index] = num_ppc_int*r;
}
#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
Expand Down Expand Up @@ -1953,7 +1998,7 @@ PhysicalParticleContainer::AddPlasmaFlux (const PlasmaInjector * plasma_injector
// For cylindrical emission (flux_normal_axis==0
// or flux_normal_axis==2), the emission surface depends on
// the radius ; thus, the calculation is finalized here
Real t_weight = flux * scale_fac * dt;
Real t_weight = flux * scale_fac;
if (loc_flux_normal_axis != 1) {
if (radially_weighted) {
t_weight *= 2._rt*MathConst::pi*radial_position;
Expand All @@ -1966,7 +2011,7 @@ PhysicalParticleContainer::AddPlasmaFlux (const PlasmaInjector * plasma_injector
}
const Real weight = t_weight;
#else
const Real weight = flux * scale_fac * dt;
const Real weight = flux * scale_fac;
#endif
pa[PIdx::w ][ip] = weight;
pa[PIdx::ux][ip] = pu.x;
Expand Down Expand Up @@ -2000,7 +2045,7 @@ PhysicalParticleContainer::AddPlasmaFlux (const PlasmaInjector * plasma_injector
if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
{
wt = static_cast<amrex::Real>(amrex::second()) - wt;
amrex::HostDevice::Atomic::Add( &(*cost)[mfi.index()], wt);
amrex::HostDevice::Atomic::Add( &(*cost)[pti.index()], wt);
}
}

Expand Down

0 comments on commit 9d48854

Please sign in to comment.