aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Source/LaserParticleContainer.H2
-rw-r--r--Source/LaserParticleContainer.cpp2
-rw-r--r--Source/ParticleContainer.H2
-rw-r--r--Source/ParticleContainer.cpp5
-rw-r--r--Source/PhysicalParticleContainer.H1
-rw-r--r--Source/PhysicalParticleContainer.cpp116
-rw-r--r--Source/RigidInjectedParticleContainer.H1
-rw-r--r--Source/RigidInjectedParticleContainer.cpp4
-rw-r--r--Source/WarpX.H4
-rw-r--r--Source/WarpX.cpp6
-rw-r--r--Source/WarpXComm.cpp31
-rw-r--r--Source/WarpXEvolve.cpp2
-rw-r--r--Source/WarpXParticleContainer.H2
-rw-r--r--Source/WarpXRegrid.cpp8
14 files changed, 140 insertions, 46 deletions
diff --git a/Source/LaserParticleContainer.H b/Source/LaserParticleContainer.H
index b5dcf72a6..8fe012fc0 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 10a9117df..92b64b28c 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*,
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 55d064cd2..8a4033795 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 c40e0c1bb..02de0a4db 100644
--- a/Source/ParticleContainer.cpp
+++ b/Source/ParticleContainer.cpp
@@ -197,7 +197,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)
@@ -209,9 +209,10 @@ 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, 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 1f3e28c9c..52e59713b 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 403606bfb..ddb066301 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)
@@ -923,6 +923,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
//
@@ -936,46 +938,94 @@ 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;
- const int *rholen;
- FArrayBox& rhofab = (*rhomf)[pti];
Box tile_box = convert(pti.tilebox(), IntVect::TheUnitVector());
- Box grown_box;
- const std::array<Real, 3>& 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<Real, 3>& 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;
+ 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);
+ warpx_charge_deposition(data_ptr, &np_current,
+ 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);
+ 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<Real,3>& 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;
+#else
+ const long nx = rholen[0]-1-2*ngRho;
+ const long ny = 0;
+ const long nz = rholen[1]-1-2*ngRho;
+#endif
+
+ 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 = (*crhomf)[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);
+ if (rho) depositCharge(rho, crho, 0);
if (! do_not_push)
{
@@ -1123,8 +1173,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<Real, 3>& xyzmin = xyzmin_tile;
tbx.grow(ngJ);
@@ -1236,7 +1284,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();
diff --git a/Source/RigidInjectedParticleContainer.H b/Source/RigidInjectedParticleContainer.H
index 18e9743ce..f2440ad7a 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 a3e2968e4..96bdc59a4 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 2fd758159..19b3d74ad 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -379,8 +379,10 @@ private:
amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3 > > Bfield_cax;
amrex::Vector<std::unique_ptr<amrex::iMultiFab> > current_buffer_masks;
amrex::Vector<std::unique_ptr<amrex::iMultiFab> > gather_buffer_masks;
- // If current deposition buffers are used
+
+ // If charge/current deposition buffers are used
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_buf;
+ amrex::Vector<std::unique_ptr<amrex::MultiFab> > charge_buf;
// div E cleaning
int do_dive_cleaning = 0;
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index ea54db2b7..9ca17c759 100644
--- a/Source/WarpX.cpp
+++ b/Source/WarpX.cpp
@@ -161,6 +161,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();
@@ -644,6 +647,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));
+ 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, 0) );
}
}
diff --git a/Source/WarpXComm.cpp b/Source/WarpXComm.cpp
index 435a430b0..23ce0ce0b 100644
--- a/Source/WarpXComm.cpp
+++ b/Source/WarpXComm.cpp
@@ -499,7 +499,8 @@ WarpX::SyncRho (amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rhof,
Vector<std::unique_ptr<MultiFab> > rho_f_g(finest_level+1);
Vector<std::unique_ptr<MultiFab> > rho_c_g(finest_level+1);
-
+ Vector<std::unique_ptr<MultiFab> > rho_buf_g(finest_level+1);
+
if (WarpX::use_filter) {
for (int lev = 0; lev <= finest_level; ++lev) {
const int ncomp = rhof[lev]->nComp();
@@ -521,6 +522,18 @@ WarpX::SyncRho (amrex::Vector<std::unique_ptr<amrex::MultiFab> >& 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
@@ -537,7 +550,14 @@ WarpX::SyncRho (amrex::Vector<std::unique_ptr<amrex::MultiFab> >& 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
@@ -556,6 +576,13 @@ WarpX::SyncRho (amrex::Vector<std::unique_ptr<amrex::MultiFab> >& 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
diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp
index 68eb229ae..be5152fc2 100644
--- a/Source/WarpXEvolve.cpp
+++ b/Source/WarpXEvolve.cpp
@@ -688,7 +688,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 8eaff7241..2138a07a5 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 360033635..744dd18e8 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<MultiFab>(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();
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-21 19:15:27 +0000