Skip to content

Commit

Permalink
add geometry for recording plane and ensure particles copied to buffe…
Browse files Browse the repository at this point in the history 
…r are made valid
  • Loading branch information
@RevathiJambunathan
RevathiJambunathan committed Feb 7, 2024
1 parent 62c8c40 commit 7e8d149
Show file tree
Hide file tree
Showing 3 changed files with 80 additions and 9 deletions.
4 changes: 4 additions & 0 deletions Source/Diagnostics/ComputeDiagFunctors/RecordingPlaneParticleFunctor.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -126,6 +126,10 @@ struct PlaneCrossingTime
#else
amrex::ignore_unused(xp, yp, zp);
#endif
// The particle may have crossed the computational domain and could be set to have a neg id
// But since the particle has crossed the recording plane location, this particle must be stored
// in the buffer
if (dst.m_aos[i_dst].id() < 0) dst.m_aos[i_dst].id() = dst.m_aos[i_dst].id()*-1;
dst.m_rdata[PIdx::w][i_dst] = m_wpnew[i_src];
dst.m_rdata[PIdx::ux][i_dst] = uxp;
dst.m_rdata[PIdx::uy][i_dst] = uyp;
Expand Down
4 changes: 4 additions & 0 deletions Source/Diagnostics/RecordingPlaneDiagnostics.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,12 +69,16 @@ private:

void ClearParticlesBuffer (int i_buffer);

void DefineSnapshotGeometry(const int i_buffer, const int lev);

int m_buffer_size = 256;
int m_slice_counter = 0;
amrex::Real m_tmin = std::numeric_limits<amrex::Real>::lowest();
amrex::Real m_tmax = std::numeric_limits<amrex::Real>::lowest();
amrex::Box m_buffer_box;
int m_flush_counter = 0;
amrex::Vector< amrex::Vector< amrex::Geometry> > m_geom_snapshot;
amrex::Vector<int> m_snapshot_geometry_defined;

[[nodiscard]] bool GetZSliceInDomain (int lev) const;
bool m_last_timeslice_filled = false;
Expand Down
81 changes: 72 additions & 9 deletions Source/Diagnostics/RecordingPlaneDiagnostics.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -133,6 +133,13 @@ RecordingPlaneDiagnostics::ReadParameters ()
m_particles_buffer.resize(m_num_buffers);
m_totalParticles_in_buffer.resize(m_num_buffers);
m_totalParticles_flushed_already.resize(m_num_buffers);
m_geom_snapshot.resize(m_num_buffers);
m_snapshot_geometry_defined.resize(m_num_buffers);
for(int i=0; i < m_num_buffers; ++i ){
int nlev_output = 1;
m_geom_snapshot[i].resize(nlev_output);
m_snapshot_geometry_defined[i] = 0;
}
}

bool
Expand Down Expand Up @@ -185,7 +192,7 @@ RecordingPlaneDiagnostics::InitializeBufferData (int /*i_buffer*/, int lev, bool
m_buffer_box = domain;
amrex::BoxArray diag_ba;
diag_ba.define(m_buffer_box);
amrex::BoxArray ba = diag_ba.maxSize(256);
amrex::BoxArray ba = diag_ba.maxSize(m_buffer_size);
amrex::IntVect typ(1);
typ[WARPX_ZINDEX] = 0; // nodal except in z-direction
ba.convert(typ);
Expand All @@ -206,6 +213,9 @@ RecordingPlaneDiagnostics::PrepareFieldDataForOutput ()
f->PrepareFunctorData(0, ZSliceInDomain, m_station_loc, m_buffer_box,
m_slice_counter, m_buffer_size);
}
for( int i_buffer=0; i_buffer < m_num_buffers; ++i_buffer ) {
DefineSnapshotGeometry(i_buffer,lev);
}
}
}

Expand Down Expand Up @@ -277,19 +287,21 @@ RecordingPlaneDiagnostics::PrepareParticleDataForOutput ()
{
const int lev = 0;
auto& warpx = WarpX::GetInstance();

for (int i = 0; i < m_particles_buffer.size(); ++i) {
DefineSnapshotGeometry(i,lev);
for (int isp = 0; isp < m_all_particle_functors.size(); ++isp )
{
amrex::Box domain = (warpx.boxArray(lev)).minimalBox();
domain.setSmall(WARPX_ZINDEX, 0);
domain.setBig(WARPX_ZINDEX, (m_buffer_size - 1));
domain.setBig(WARPX_ZINDEX, (1));
amrex::BoxArray diag_ba;
diag_ba.define(m_buffer_box);
amrex::BoxArray ba = diag_ba.maxSize(256);
diag_ba.define(domain);
amrex::BoxArray ba = diag_ba.maxSize(m_buffer_size);
amrex::DistributionMapping dmap(ba);
m_particles_buffer[i][isp]->SetParticleBoxArray(lev, ba);
m_particles_buffer[i][isp]->SetParticleDistributionMap(lev, dmap);
m_particles_buffer[i][isp]->SetParticleGeometry(lev, warpx.Geom(0));
m_particles_buffer[i][isp]->SetParticleGeometry(lev, m_geom_snapshot[0][0]);
}
}
for (int isp = 0; isp < m_all_particle_functors.size(); ++isp)
Expand All @@ -298,6 +310,29 @@ RecordingPlaneDiagnostics::PrepareParticleDataForOutput ()
}
}

void
RecordingPlaneDiagnostics::DefineSnapshotGeometry (const int i_buffer, const int lev)
{
if (m_snapshot_geometry_defined[i_buffer]) { return; }

if (lev == 0) {
auto& warpx = WarpX::GetInstance();
amrex::Box domain = (warpx.boxArray(lev)).minimalBox();
domain.setSmall(WARPX_ZINDEX,0);
domain.setBig(WARPX_ZINDEX, (domain.smallEnd(WARPX_ZINDEX)+1));
amrex::Real zlo = m_station_loc - warpx.Geom(lev).CellSize(WARPX_ZINDEX)*1.1;
amrex::Real zhi = m_station_loc + warpx.Geom(lev).CellSize(WARPX_ZINDEX)*1.1;
amrex::RealBox realdom = warpx.Geom(lev).ProbDomain();
realdom.setLo(WARPX_ZINDEX,zlo);
realdom.setHi(WARPX_ZINDEX,zhi);
amrex::Vector<int> periodicity(AMREX_SPACEDIM, 0);
m_geom_snapshot[i_buffer][lev].define(domain, &realdom, amrex::CoordSys::cartesian, periodicity.data());
} else if (lev > 0) {
m_geom_snapshot[i_buffer][lev] = amrex::refine(m_geom_snapshot[i_buffer][lev-1], WarpX::RefRatio(lev-1));
}
}


void
RecordingPlaneDiagnostics::Flush (int i_buffer, bool /* force_flush */)
{
Expand All @@ -320,6 +355,10 @@ RecordingPlaneDiagnostics::Flush (int i_buffer, bool /* force_flush */)
}
}
amrex::ParallelDescriptor::Barrier();
amrex::Vector<amrex::BoxArray> vba;
amrex::Vector<amrex::DistributionMapping> vdmap;
amrex::Vector<amrex::Geometry> vgeom;
amrex::Vector<amrex::IntVect> vrefratio;
for (int lev = 0; lev < 1; ++lev) {
std::string buffer_string = amrex::Concatenate("buffer-",m_flush_counter,1);
amrex::Print() << " Writing station buffer " << buffer_string << "\n";
Expand Down Expand Up @@ -377,20 +416,44 @@ RecordingPlaneDiagnostics::Flush (int i_buffer, bool /* force_flush */)
// particles buffer is saved in pinned particle container
bool const use_pinned_pc = true;
auto const& warpx = WarpX::GetInstance();
const amrex::Geometry& geom = warpx.Geom(0); // not used in WriteToFile
bool const plot_raw_fields = false;
bool const plot_raw_fields_guards = false;

if (nparticles > 0) {

const int nlevels = m_particles_buffer[i_buffer][0]->numLevels();
for (int lev = 0 ; lev < nlevels; ++lev) {
// Store BoxArray, dmap, geometry, and refratio for every level
vba.push_back(m_particles_buffer[i_buffer][0]->ParticleBoxArray(lev));
vdmap.push_back(m_particles_buffer[i_buffer][0]->ParticleDistributionMap(lev));
vgeom.push_back(m_particles_buffer[i_buffer][0]->ParticleGeom(lev));
if (lev < nlevels - 1) {
vrefratio.push_back(m_particles_buffer[i_buffer][0]->GetParGDB()->refRatio(lev));
}
}

amrex::Box domain = (m_particles_buffer[i_buffer][0]->ParticleBoxArray(0)).minimalBox();
domain.setSmall(WARPX_ZINDEX, (domain.smallEnd(WARPX_ZINDEX)));
domain.setBig(WARPX_ZINDEX, (domain.bigEnd(WARPX_ZINDEX)));
amrex::BoxArray diag_ba;
diag_ba.define(domain);
m_particles_buffer[i_buffer][0]->SetParticleBoxArray(0,diag_ba);
for (int isp = 0; isp < m_particles_buffer.at(i_buffer).size(); ++isp) {
m_particles_buffer[i_buffer][isp]->SetParGDB(vgeom[0], vdmap[0], diag_ba);
m_particles_buffer[i_buffer][isp]->Redistribute();
}
m_particles_buffer[i_buffer][0]->SetParticleBoxArray(0,vba.back());
m_particles_buffer[i_buffer][0]->SetParGDB(vgeom[0], vdmap[0], vba.back());

m_flush_format->WriteToFile(m_varnames, {}, m_geom_output[i_buffer],
warpx.getistep(), warpx.gett_new(0), m_output_species[i_buffer], nlev_output, filename,
m_file_min_digits, plot_raw_fields, plot_raw_fields_guards, use_pinned_pc,
isBTD, i_buffer, m_flush_counter, maxBTDBuffers, geom,
isBTD, i_buffer, m_flush_counter, maxBTDBuffers, m_geom_snapshot[i_buffer][0],
isLastBTD, m_totalParticles_flushed_already[i_buffer]);
}

for (int isp = 0; isp < m_particles_buffer[0].size(); ++isp) {
m_totalParticles_flushed_already[i_buffer][isp] += int(m_particles_buffer[i_buffer][isp]->TotalNumberOfParticles());
m_totalParticles_flushed_already[i_buffer][isp] += m_totalParticles_in_buffer[i_buffer][isp];
m_totalParticles_in_buffer[i_buffer][isp] = 0;
}

Expand All @@ -409,7 +472,7 @@ RecordingPlaneDiagnostics::Flush (int i_buffer, bool /* force_flush */)
// slice
for (int lev = 0; lev < 1; ++lev) {
auto& out = m_mf_output[i_buffer][lev];
amrex::Box destbox = WarpX::GetInstance().Geom(lev).Domain();
amrex::Box destbox = m_geom_snapshot[i_buffer][lev].Domain();
destbox.convert(out.ixType()).setRange(WARPX_ZINDEX,0);
amrex::NonLocalBC::MultiBlockIndexMapping dtos;
dtos.offset[WARPX_ZINDEX] = -(m_buffer_size-1);
Expand Down

0 comments on commit 7e8d149

Please sign in to comment.