aboutsummaryrefslogtreecommitdiff
path: root/Source
diff options
context:
space:
mode:
Diffstat (limited to 'Source')
-rw-r--r--Source/Diagnostics/FieldIO.H1
-rw-r--r--Source/Diagnostics/FieldIO.cpp184
-rw-r--r--Source/Diagnostics/WarpXIO.cpp23
-rw-r--r--Source/Evolve/WarpXEvolveEM.cpp27
-rw-r--r--Source/FortranInterface/WarpX_f.H3
-rw-r--r--Source/FortranInterface/WarpX_picsar.F9059
-rw-r--r--Source/Initialization/PlasmaInjector.cpp6
-rw-r--r--Source/Parallelization/WarpXComm.cpp90
-rw-r--r--Source/Parallelization/WarpXRegrid.cpp66
-rw-r--r--Source/Particles/MultiParticleContainer.cpp2
-rw-r--r--Source/Particles/PhysicalParticleContainer.H1
-rw-r--r--Source/Particles/PhysicalParticleContainer.cpp11
-rw-r--r--Source/Particles/RigidInjectedParticleContainer.cpp3
-rw-r--r--Source/Particles/WarpXParticleContainer.H3
-rw-r--r--Source/Particles/WarpXParticleContainer.cpp27
-rw-r--r--Source/WarpX.H4
-rw-r--r--Source/WarpX.cpp112
17 files changed, 435 insertions, 187 deletions
diff --git a/Source/Diagnostics/FieldIO.H b/Source/Diagnostics/FieldIO.H
index 1a3b45580..24fd6abb6 100644
--- a/Source/Diagnostics/FieldIO.H
+++ b/Source/Diagnostics/FieldIO.H
@@ -15,6 +15,7 @@ PackPlotDataPtrs (Vector<const MultiFab*>& pmf,
void
AverageAndPackVectorField( MultiFab& mf_avg,
const std::array< std::unique_ptr<MultiFab>, 3 >& vector_field,
+ const DistributionMapping& dm,
const int dcomp, const int ngrow );
void
diff --git a/Source/Diagnostics/FieldIO.cpp b/Source/Diagnostics/FieldIO.cpp
index 5c523804e..a1fd61b25 100644
--- a/Source/Diagnostics/FieldIO.cpp
+++ b/Source/Diagnostics/FieldIO.cpp
@@ -204,6 +204,20 @@ WriteOpenPMDFields( const std::string& filename,
}
#endif // WARPX_USE_OPENPMD
+void
+ConstructTotalRZField(std::array< std::unique_ptr<MultiFab>, 3 >& mf_total,
+ const std::array< std::unique_ptr<MultiFab>, 3 >& vector_field)
+{
+ // Sum over the real components, giving quantity at theta=0
+ MultiFab::Copy(*mf_total[0], *vector_field[0], 0, 0, 1, vector_field[0]->nGrowVect());
+ MultiFab::Copy(*mf_total[1], *vector_field[1], 0, 0, 1, vector_field[1]->nGrowVect());
+ MultiFab::Copy(*mf_total[2], *vector_field[2], 0, 0, 1, vector_field[2]->nGrowVect());
+ for (int ic=2 ; ic < vector_field[0]->nComp() ; ic += 2) {
+ MultiFab::Add(*mf_total[0], *vector_field[0], ic, 0, 1, vector_field[0]->nGrowVect());
+ MultiFab::Add(*mf_total[1], *vector_field[1], ic, 0, 1, vector_field[1]->nGrowVect());
+ MultiFab::Add(*mf_total[2], *vector_field[2], ic, 0, 1, vector_field[2]->nGrowVect());
+ }
+}
void
PackPlotDataPtrs (Vector<const MultiFab*>& pmf,
@@ -228,6 +242,7 @@ PackPlotDataPtrs (Vector<const MultiFab*>& pmf,
void
AverageAndPackVectorField( MultiFab& mf_avg,
const std::array< std::unique_ptr<MultiFab>, 3 >& vector_field,
+ const DistributionMapping& dm,
const int dcomp, const int ngrow )
{
// The object below is temporary, and is needed because
@@ -256,23 +271,101 @@ AverageAndPackVectorField( MultiFab& mf_avg,
// - Face centered, in the same way as B on a Yee grid
} else if ( vector_field[0]->is_nodal(0) ){
- PackPlotDataPtrs(srcmf, vector_field);
- amrex::average_face_to_cellcenter( mf_avg, dcomp, srcmf, ngrow);
+ // Note that average_face_to_cellcenter operates only on the number of
+ // arrays equal to the number of dimensions. So, for 2D, PackPlotDataPtrs
+ // packs in the x and z (or r and z) arrays, which are then cell averaged.
+ // The Copy code then copies the z from the 2nd to the 3rd field,
+ // and copies over directly the y (or theta) component (which is
+ // already cell centered).
+ if (vector_field[0]->nComp() > 1) {
+#ifdef WARPX_RZ
+ // When there are more than one components, the total
+ // fields needs to be constructed in temporary MultiFabs
+ // Note that mf_total is declared in the same way as
+ // vector_field so that it can be passed into PackPlotDataPtrs.
+ std::array<std::unique_ptr<MultiFab>,3> mf_total;
+ mf_total[0].reset(new MultiFab(vector_field[0]->boxArray(), dm, 1, vector_field[0]->nGrowVect()));
+ mf_total[1].reset(new MultiFab(vector_field[1]->boxArray(), dm, 1, vector_field[1]->nGrowVect()));
+ mf_total[2].reset(new MultiFab(vector_field[2]->boxArray(), dm, 1, vector_field[2]->nGrowVect()));
+ ConstructTotalRZField(mf_total, vector_field);
+ PackPlotDataPtrs(srcmf, mf_total);
+ amrex::average_face_to_cellcenter( mf_avg, dcomp, srcmf, ngrow);
+ MultiFab::Copy( mf_avg, mf_avg, dcomp+1, dcomp+2, 1, ngrow);
+ MultiFab::Copy( mf_avg, *mf_total[1], 0, dcomp+1, 1, ngrow);
+ // Also add the real and imaginary parts of each mode.
+ /* Not supported yet
+ for (int i=0 ; i < vector_field[0]->nComp() ; i++) {
+ MultiFab v_comp0(*vector_field[0], amrex::make_alias, i, 1);
+ MultiFab v_comp1(*vector_field[1], amrex::make_alias, i, 1);
+ MultiFab v_comp2(*vector_field[2], amrex::make_alias, i, 1);
+ srcmf[0] = &v_comp0;
+ srcmf[1] = &v_comp2;
+ int id = dcomp + 3*(i + 1);
+ amrex::average_face_to_cellcenter( mf_avg, id, srcmf, ngrow);
+ MultiFab::Copy( mf_avg, mf_avg, id+1, id+2, 1, ngrow);
+ MultiFab::Copy( mf_avg, v_comp1, 0, id+1, 1, ngrow);
+ }
+ */
+#else
+ amrex::Abort("AverageAndPackVectorField not implemented for ncomp > 1");
+#endif
+ } else {
+ PackPlotDataPtrs(srcmf, vector_field);
+ amrex::average_face_to_cellcenter( mf_avg, dcomp, srcmf, ngrow);
#if (AMREX_SPACEDIM == 2)
- MultiFab::Copy( mf_avg, mf_avg, dcomp+1, dcomp+2, 1, ngrow);
- MultiFab::Copy( mf_avg, *vector_field[1], 0, dcomp+1, 1, ngrow);
+ MultiFab::Copy( mf_avg, mf_avg, dcomp+1, dcomp+2, 1, ngrow);
+ MultiFab::Copy( mf_avg, *vector_field[1], 0, dcomp+1, 1, ngrow);
#endif
+ }
// - Edge centered, in the same way as E on a Yee grid
} else if ( !vector_field[0]->is_nodal(0) ){
- PackPlotDataPtrs(srcmf, vector_field);
- amrex::average_edge_to_cellcenter( mf_avg, dcomp, srcmf, ngrow);
+ // See comment above, though here, the y (or theta) component
+ // has node centering.
+ if (vector_field[0]->nComp() > 1) {
+#ifdef WARPX_RZ
+ // When there are more than one components, the total
+ // fields needs to be constructed in temporary MultiFabs
+ // Note that mf_total is declared in the same way as
+ // vector_field so that it can be passed into PackPlotDataPtrs.
+ std::array<std::unique_ptr<MultiFab>,3> mf_total;
+ mf_total[0].reset(new MultiFab(vector_field[0]->boxArray(), dm, 1, vector_field[0]->nGrowVect()));
+ mf_total[1].reset(new MultiFab(vector_field[1]->boxArray(), dm, 1, vector_field[1]->nGrowVect()));
+ mf_total[2].reset(new MultiFab(vector_field[2]->boxArray(), dm, 1, vector_field[2]->nGrowVect()));
+ ConstructTotalRZField(mf_total, vector_field);
+ PackPlotDataPtrs(srcmf, mf_total);
+ amrex::average_edge_to_cellcenter( mf_avg, dcomp, srcmf, ngrow);
+ MultiFab::Copy( mf_avg, mf_avg, dcomp+1, dcomp+2, 1, ngrow);
+ amrex::average_node_to_cellcenter( mf_avg, dcomp+1,
+ *mf_total[1], 0, 1, ngrow);
+ // Also add the real and imaginary parts of each mode.
+ /* Not supported yet
+ for (int i=0 ; i < vector_field[0]->nComp() ; i++) {
+ MultiFab v_comp0(*vector_field[0], amrex::make_alias, i, 1);
+ MultiFab v_comp1(*vector_field[1], amrex::make_alias, i, 1);
+ MultiFab v_comp2(*vector_field[2], amrex::make_alias, i, 1);
+ srcmf[0] = &v_comp0;
+ srcmf[1] = &v_comp2;
+ int id = dcomp + 3*(i + 1);
+ amrex::average_edge_to_cellcenter( mf_avg, id, srcmf, ngrow);
+ MultiFab::Copy( mf_avg, mf_avg, id+1, id+2, 1, ngrow);
+ amrex::average_node_to_cellcenter( mf_avg, id+1,
+ v_comp1, 0, 1, ngrow);
+ }
+ */
+#else
+ amrex::Abort("AverageAndPackVectorField not implemented for ncomp > 1");
+#endif
+ } else {
+ PackPlotDataPtrs(srcmf, vector_field);
+ amrex::average_edge_to_cellcenter( mf_avg, dcomp, srcmf, ngrow);
#if (AMREX_SPACEDIM == 2)
- MultiFab::Copy( mf_avg, mf_avg, dcomp+1, dcomp+2, 1, ngrow);
- amrex::average_node_to_cellcenter( mf_avg, dcomp+1,
- *vector_field[1], 0, 1, ngrow);
+ MultiFab::Copy( mf_avg, mf_avg, dcomp+1, dcomp+2, 1, ngrow);
+ amrex::average_node_to_cellcenter( mf_avg, dcomp+1,
+ *vector_field[1], 0, 1, ngrow);
#endif
+ }
} else {
amrex::Abort("Unknown staggering.");
@@ -305,6 +398,35 @@ AverageAndPackScalarField( MultiFab& mf_avg,
}
}
+/** \brief Add variable names to the list.
+ * If there are more that one mode, add the
+ * name of the total field and then the
+ * names of the real and imaginary parts of each
+ * mode.
+ */
+void
+AddToVarNames( Vector<std::string>& varnames,
+ std::string name, std::string suffix,
+ int n_rz_azimuthal_modes ) {
+ auto coords = {"x", "y", "z"};
+ for(auto coord:coords) varnames.push_back(name+coord+suffix);
+ /* Not supported yet
+ auto coordsRZ = {"r", "theta", "z"};
+ if (n_rz_azimuthal_modes > 1) {
+ // Note that the names are added in the same order as the fields
+ // are packed in AverageAndPackVectorField
+ for (int i = 0 ; i < n_rz_azimuthal_modes ; i++) {
+ for(auto coord:coordsRZ) {
+ varnames.push_back(name + coord + suffix + std::to_string(i) + "_real");
+ }
+ for(auto coord:coordsRZ) {
+ varnames.push_back(name + coord + suffix + std::to_string(i) + "_imag");
+ }
+ }
+ }
+ */
+}
+
/** \brief Write the different fields that are meant for output,
* into the vector of MultiFab `mf_avg` (one MultiFab per level)
* after averaging them to the cell centers.
@@ -313,6 +435,14 @@ void
WarpX::AverageAndPackFields ( Vector<std::string>& varnames,
amrex::Vector<MultiFab>& mf_avg, const int ngrow) const
{
+ // Factor to account for quantities that have multiple components.
+ // If n_rz_azimuthal_modes > 1, allow space for total field and the real and
+ // imaginary part of each node. For now, also include the
+ // imaginary part of mode 0 for code symmetry, even though
+ // it is always zero.
+ int modes_factor = 1;
+ if (n_rz_azimuthal_modes > 1) modes_factor = 2*n_rz_azimuthal_modes + 1;
+
// Count how many different fields should be written (ncomp)
const int ncomp = fields_to_plot.size()
+ static_cast<int>(plot_finepatch)*6
@@ -333,17 +463,17 @@ WarpX::AverageAndPackFields ( Vector<std::string>& varnames,
// Allocate temp MultiFab with 3 components
mf_tmp_E = MultiFab(grids[lev], dmap[lev], 3, ngrow);
// Fill MultiFab mf_tmp_E with averaged E
- AverageAndPackVectorField(mf_tmp_E, Efield_aux[lev], 0, ngrow);
+ AverageAndPackVectorField(mf_tmp_E, Efield_aux[lev], dmap[lev], 0, ngrow);
}
// Same for B
if (is_in_vector(fields_to_plot, {"Bx", "By", "Bz"} )){
mf_tmp_B = MultiFab(grids[lev], dmap[lev], 3, ngrow);
- AverageAndPackVectorField(mf_tmp_B, Bfield_aux[lev], 0, ngrow);
+ AverageAndPackVectorField(mf_tmp_B, Bfield_aux[lev], dmap[lev], 0, ngrow);
}
// Same for J
if (is_in_vector(fields_to_plot, {"jx", "jy", "jz"} )){
mf_tmp_J = MultiFab(grids[lev], dmap[lev], 3, ngrow);
- AverageAndPackVectorField(mf_tmp_J, current_fp[lev], 0, ngrow);
+ AverageAndPackVectorField(mf_tmp_J, current_fp[lev], dmap[lev], 0, ngrow);
}
int dcomp;
@@ -437,12 +567,12 @@ WarpX::AverageAndPackFields ( Vector<std::string>& varnames,
}
if (plot_finepatch)
{
- AverageAndPackVectorField( mf_avg[lev], Efield_fp[lev], dcomp, ngrow );
- if(lev==0) for(auto name:{"Ex_fp","Ey_fp","Ez_fp"}) varnames.push_back(name);
- dcomp += 3;
- AverageAndPackVectorField( mf_avg[lev], Bfield_fp[lev], dcomp, ngrow );
- if(lev==0) for(auto name:{"Bx_fp","By_fp","Bz_fp"}) varnames.push_back(name);
- dcomp += 3;
+ AverageAndPackVectorField( mf_avg[lev], Efield_fp[lev], dmap[lev], dcomp, ngrow );
+ if (lev == 0) AddToVarNames(varnames, "E", "_fp", n_rz_azimuthal_modes);
+ dcomp += 3*modes_factor;
+ AverageAndPackVectorField( mf_avg[lev], Bfield_fp[lev], dmap[lev], dcomp, ngrow );
+ if (lev == 0) AddToVarNames(varnames, "B", "_fp", n_rz_azimuthal_modes);
+ dcomp += 3*modes_factor;
}
if (plot_crsepatch)
@@ -455,11 +585,11 @@ WarpX::AverageAndPackFields ( Vector<std::string>& varnames,
{
if (do_nodal) amrex::Abort("TODO: do_nodal && plot_crsepatch");
std::array<std::unique_ptr<MultiFab>, 3> E = getInterpolatedE(lev);
- AverageAndPackVectorField( mf_avg[lev], E, dcomp, ngrow );
+ AverageAndPackVectorField( mf_avg[lev], E, dmap[lev], dcomp, ngrow );
}
- if(lev==0) for(auto name:{"Ex_cp","Ey_cp","Ez_cp"}) varnames.push_back(name);
- dcomp += 3;
+ if (lev == 0) AddToVarNames(varnames, "E", "_cp", n_rz_azimuthal_modes);
+ dcomp += 3*modes_factor;
// now the magnetic field
if (lev == 0)
@@ -470,10 +600,10 @@ WarpX::AverageAndPackFields ( Vector<std::string>& varnames,
{
if (do_nodal) amrex::Abort("TODO: do_nodal && plot_crsepatch");
std::array<std::unique_ptr<MultiFab>, 3> B = getInterpolatedB(lev);
- AverageAndPackVectorField( mf_avg[lev], B, dcomp, ngrow );
+ AverageAndPackVectorField( mf_avg[lev], B, dmap[lev], dcomp, ngrow );
}
- if(lev==0) for(auto name:{"Bx_cp","By_cp","Bz_cp"}) varnames.push_back(name);
- dcomp += 3;
+ if (lev == 0) AddToVarNames(varnames, "B", "_cp", n_rz_azimuthal_modes);
+ dcomp += 3*modes_factor;
}
if (costs[0] != nullptr)
@@ -535,8 +665,8 @@ WriteRawField( const MultiFab& F, const DistributionMapping& dm,
VisMF::Write(F, prefix);
} else {
// Copy original MultiFab into one that does not have guard cells
- MultiFab tmpF( F.boxArray(), dm, 1, 0);
- MultiFab::Copy(tmpF, F, 0, 0, 1, 0);
+ MultiFab tmpF( F.boxArray(), dm, F.nComp(), 0);
+ MultiFab::Copy(tmpF, F, 0, 0, F.nComp(), 0);
VisMF::Write(tmpF, prefix);
}
@@ -558,7 +688,7 @@ WriteZeroRawField( const MultiFab& F, const DistributionMapping& dm,
std::string prefix = amrex::MultiFabFileFullPrefix(lev,
filename, level_prefix, field_name);
- MultiFab tmpF(F.boxArray(), dm, 1, ng);
+ MultiFab tmpF(F.boxArray(), dm, F.nComp(), ng);
tmpF.setVal(0.);
VisMF::Write(tmpF, prefix);
}
diff --git a/Source/Diagnostics/WarpXIO.cpp b/Source/Diagnostics/WarpXIO.cpp
index 38399bf9e..2ea8a7fd7 100644
--- a/Source/Diagnostics/WarpXIO.cpp
+++ b/Source/Diagnostics/WarpXIO.cpp
@@ -415,20 +415,28 @@ WarpX::GetCellCenteredData() {
Vector<std::unique_ptr<MultiFab> > cc(finest_level+1);
+ // Factor to account for quantities that have multiple components.
+ // If n_rz_azimuthal_modes > 1, allow space for total field and the real and
+ // imaginary part of each node. For now, also include the
+ // imaginary part of mode 0 for code symmetry, even though
+ // it is always zero.
+ int modes_factor = 1;
+ if (n_rz_azimuthal_modes > 1) modes_factor = 2*n_rz_azimuthal_modes + 1;
+
for (int lev = 0; lev <= finest_level; ++lev)
{
cc[lev].reset( new MultiFab(grids[lev], dmap[lev], nc, ng) );
int dcomp = 0;
// first the electric field
- AverageAndPackVectorField( *cc[lev], Efield_aux[lev], dcomp, ng );
- dcomp += 3;
+ AverageAndPackVectorField( *cc[lev], Efield_aux[lev], dmap[lev], dcomp, ng );
+ dcomp += 3*modes_factor;
// then the magnetic field
- AverageAndPackVectorField( *cc[lev], Bfield_aux[lev], dcomp, ng );
- dcomp += 3;
+ AverageAndPackVectorField( *cc[lev], Bfield_aux[lev], dmap[lev], dcomp, ng );
+ dcomp += 3*modes_factor;
// then the current density
- AverageAndPackVectorField( *cc[lev], current_fp[lev], dcomp, ng );
- dcomp += 3;
+ AverageAndPackVectorField( *cc[lev], current_fp[lev], dmap[lev], dcomp, ng );
+ dcomp += 3*modes_factor;
// then the charge density
const std::unique_ptr<MultiFab>& charge_density = mypc->GetChargeDensity(lev);
AverageAndPackScalarField( *cc[lev], *charge_density, dcomp, ng );
@@ -582,7 +590,8 @@ WarpX::WritePlotFile () const
if (F_fp[lev]) WriteRawField( *F_fp[lev], dm, raw_pltname, level_prefix, "F_fp", lev, plot_raw_fields_guards);
if (plot_rho) {
// Use the component 1 of `rho_fp`, i.e. rho_new for time synchronization
- MultiFab rho_new(*rho_fp[lev], amrex::make_alias, 1, 1);
+ // If nComp > 1, this is the upper half of the list of components.
+ MultiFab rho_new(*rho_fp[lev], amrex::make_alias, rho_fp[lev]->nComp()/2, rho_fp[lev]->nComp()/2);
WriteRawField( rho_new, dm, raw_pltname, level_prefix, "rho_fp", lev, plot_raw_fields_guards);
}
}
diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp
index 170a8f604..a532a22fd 100644
--- a/Source/Evolve/WarpXEvolveEM.cpp
+++ b/Source/Evolve/WarpXEvolveEM.cpp
@@ -349,7 +349,7 @@ WarpX::OneStep_sub1 (Real curtime)
RestrictRhoFromFineToCoarsePatch(fine_lev);
ApplyFilterandSumBoundaryJ(fine_lev, PatchType::fine);
NodalSyncJ(fine_lev, PatchType::fine);
- ApplyFilterandSumBoundaryRho(fine_lev, PatchType::fine, 0, 2);
+ ApplyFilterandSumBoundaryRho(fine_lev, PatchType::fine, 0, 2*n_rz_azimuthal_modes);
NodalSyncRho(fine_lev, PatchType::fine, 0, 2);
EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]);
@@ -376,7 +376,7 @@ WarpX::OneStep_sub1 (Real curtime)
PushParticlesandDepose(coarse_lev, curtime);
StoreCurrent(coarse_lev);
AddCurrentFromFineLevelandSumBoundary(coarse_lev);
- AddRhoFromFineLevelandSumBoundary(coarse_lev, 0, 1);
+ AddRhoFromFineLevelandSumBoundary(coarse_lev, 0, n_rz_azimuthal_modes);
EvolveB(fine_lev, PatchType::coarse, dt[fine_lev]);
EvolveF(fine_lev, PatchType::coarse, dt[fine_lev], DtType::FirstHalf);
@@ -403,7 +403,7 @@ WarpX::OneStep_sub1 (Real curtime)
RestrictRhoFromFineToCoarsePatch(fine_lev);
ApplyFilterandSumBoundaryJ(fine_lev, PatchType::fine);
NodalSyncJ(fine_lev, PatchType::fine);
- ApplyFilterandSumBoundaryRho(fine_lev, PatchType::fine, 0, 2);
+ ApplyFilterandSumBoundaryRho(fine_lev, PatchType::fine, 0, 2*n_rz_azimuthal_modes);
NodalSyncRho(fine_lev, PatchType::fine, 0, 2);
EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]);
@@ -429,7 +429,7 @@ WarpX::OneStep_sub1 (Real curtime)
// by only half a coarse step (second half)
RestoreCurrent(coarse_lev);
AddCurrentFromFineLevelandSumBoundary(coarse_lev);
- AddRhoFromFineLevelandSumBoundary(coarse_lev, 1, 1);
+ AddRhoFromFineLevelandSumBoundary(coarse_lev, n_rz_azimuthal_modes, n_rz_azimuthal_modes);
EvolveE(fine_lev, PatchType::coarse, dt[fine_lev]);
FillBoundaryE(fine_lev, PatchType::coarse);
@@ -506,8 +506,23 @@ WarpX::ComputeDt ()
if (maxwell_fdtd_solver_id == 0) {
// CFL time step Yee solver
#ifdef WARPX_DIM_RZ
- // Derived semi-analytically by R. Lehe
- deltat = cfl * 1./( std::sqrt((1+0.2105)/(dx[0]*dx[0]) + 1./(dx[1]*dx[1])) * PhysConst::c );
+ // In the rz case, the Courant limit has been evaluated
+ // semi-analytically by R. Lehe, and resulted in the following
+ // coefficients. For an explanation, see (not officially published)
+ // www.normalesup.org/~lehe/Disp_relation_Circ.pdf
+ // NB : Here the coefficient for m=1 as compared to this document,
+ // as it was observed in practice that this coefficient was not
+ // high enough (The simulation became unstable).
+ Real multimode_coeffs[6] = { 0.2105, 1.0, 3.5234, 8.5104, 15.5059, 24.5037 };
+ Real multimode_alpha;
+ if (n_rz_azimuthal_modes < 7) {
+ // Use the table of the coefficients
+ multimode_alpha = multimode_coeffs[n_rz_azimuthal_modes-1];
+ } else {
+ // Use a realistic extrapolation
+ multimode_alpha = (n_rz_azimuthal_modes - 1)*(n_rz_azimuthal_modes - 1) - 0.4;
+ }
+ deltat = cfl * 1./( std::sqrt((1+multimode_alpha)/(dx[0]*dx[0]) + 1./(dx[1]*dx[1])) * PhysConst::c );
#else
deltat = cfl * 1./( std::sqrt(AMREX_D_TERM( 1./(dx[0]*dx[0]),
+ 1./(dx[1]*dx[1]),
diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H
index aac23f781..6dff3469d 100644
--- a/Source/FortranInterface/WarpX_f.H
+++ b/Source/FortranInterface/WarpX_f.H
@@ -80,6 +80,7 @@ extern "C"
amrex::Real* jx, const long* jx_ng, const int* jx_ntot,
amrex::Real* jy, const long* jy_ng, const int* jy_ntot,
amrex::Real* jz, const long* jz_ng, const int* jz_ntot,
+ const long* n_rz_azimuthal_modes,
const long* np,
const amrex::Real* xp, const amrex::Real* yp, const amrex::Real* zp,
const amrex::Real* uxp, const amrex::Real* uyp,const amrex::Real* uzp,
@@ -144,6 +145,7 @@ extern "C"
const int* xlo, const int* xhi,
const int* ylo, const int* yhi,
const int* zlo, const int* zhi,
+ const long* n_rz_azimuthal_modes,
BL_FORT_FAB_ARG_3D(ex),
BL_FORT_FAB_ARG_3D(ey),
BL_FORT_FAB_ARG_3D(ez),
@@ -164,6 +166,7 @@ extern "C"
const int* xlo, const int* xhi,
const int* ylo, const int* yhi,
const int* zlo, const int* zhi,
+ const long* n_rz_azimuthal_modes,
const BL_FORT_FAB_ARG_3D(ex),
const BL_FORT_FAB_ARG_3D(ey),
const BL_FORT_FAB_ARG_3D(ez),
diff --git a/Source/FortranInterface/WarpX_picsar.F90 b/Source/FortranInterface/WarpX_picsar.F90
index 34084d753..3a9f8f41e 100644
--- a/Source/FortranInterface/WarpX_picsar.F90
+++ b/Source/FortranInterface/WarpX_picsar.F90
@@ -74,17 +74,25 @@ contains
!> @param[in] charge_depo_algo algorithm choice for the charge deposition
!>
subroutine warpx_current_deposition( &
- jx,jx_ng,jx_ntot,jy,jy_ng,jy_ntot,jz,jz_ng,jz_ntot, &
+ jx,jx_ng,jx_ntot,jy,jy_ng,jy_ntot,jz,jz_ng,jz_ntot,n_rz_azimuthal_modes, &
np,xp,yp,zp,uxp,uyp,uzp,gaminv,w,q,xmin,ymin,zmin,dt,dx,dy,dz,nox,noy,noz,&
l_nodal,lvect,current_depo_algo) &
bind(C, name="warpx_current_deposition")
integer, intent(in) :: jx_ntot(AMREX_SPACEDIM), jy_ntot(AMREX_SPACEDIM), jz_ntot(AMREX_SPACEDIM)
integer(c_long), intent(in) :: jx_ng, jy_ng, jz_ng
+ integer(c_long), intent(IN) :: n_rz_azimuthal_modes
integer(c_long), intent(IN) :: np
integer(c_long), intent(IN) :: nox,noy,noz
integer(c_int), intent(in) :: l_nodal
+
+#ifdef WARPX_RZ
+ real(amrex_real), intent(IN OUT):: jx(jx_ntot(1),jx_ntot(2),2,n_rz_azimuthal_modes)
+ real(amrex_real), intent(IN OUT):: jy(jy_ntot(1),jy_ntot(2),2,n_rz_azimuthal_modes)
+ real(amrex_real), intent(IN OUT):: jz(jz_ntot(1),jz_ntot(2),2,n_rz_azimuthal_modes)
+#else
real(amrex_real), intent(IN OUT):: jx(*), jy(*), jz(*)
+#endif
real(amrex_real), intent(IN) :: q
real(amrex_real), intent(IN) :: dx,dy,dz
real(amrex_real), intent(IN) :: dt
@@ -96,6 +104,13 @@ contains
integer(c_long), intent(IN) :: current_depo_algo
logical(pxr_logical) :: pxr_l_nodal
+#ifdef WARPX_RZ
+ logical(pxr_logical) :: l_particles_weight = .true.
+ integer(c_long) :: type_rz_depose = 1
+ complex(amrex_real), allocatable, dimension(:,:,:) :: jr_c, jt_c, jz_c
+ integer :: alloc_status
+#endif
+
! Compute the number of valid cells and guard cells
integer(c_long) :: jx_nvalid(AMREX_SPACEDIM), jy_nvalid(AMREX_SPACEDIM), jz_nvalid(AMREX_SPACEDIM), &
jx_nguards(AMREX_SPACEDIM), jy_nguards(AMREX_SPACEDIM), jz_nguards(AMREX_SPACEDIM)
@@ -118,13 +133,53 @@ contains
nox,noy,noz,pxr_l_nodal,current_depo_algo)
! Dimension 2
#elif (AMREX_SPACEDIM==2)
- CALL WRPX_PXR_CURRENT_DEPOSITION( &
+#ifdef WARPX_RZ
+ if (n_rz_azimuthal_modes > 1) then
+
+ allocate(jr_c(jx_ntot(1),jx_ntot(2),n_rz_azimuthal_modes), &
+ jt_c(jy_ntot(1),jy_ntot(2),n_rz_azimuthal_modes), &
+ jz_c(jz_ntot(1),jz_ntot(2),n_rz_azimuthal_modes), stat=alloc_status)
+ if (alloc_status /= 0) then
+ print*,"Error: warpx_current_deposition: complex arrays could not be allocated"
+ stop
+ endif
+
+ jr_c = 0._amrex_real
+ jt_c = 0._amrex_real
+ jz_c = 0._amrex_real
+
+ CALL pxr_depose_jrjtjz_esirkepov_n_2d_circ( &
+ jr_c,jx_nguards,jx_nvalid, &
+ jt_c,jy_nguards,jy_nvalid, &
+ jz_c,jz_nguards,jz_nvalid, &
+ n_rz_azimuthal_modes, &
+ np,xp,yp,zp,uxp,uyp,uzp,gaminv,w,q, &
+ xmin,zmin,dt,dx,dz, &
+ nox,noz,l_particles_weight,type_rz_depose)
+
+ jx(:,:,1,:) = jx(:,:,1,:) + real(jr_c)
+ jx(:,:,2,:) = jx(:,:,2,:) + aimag(jr_c)
+ jy(:,:,1,:) = jy(:,:,1,:) + real(jt_c)
+ jy(:,:,2,:) = jy(:,:,2,:) + aimag(jt_c)
+ jz(:,:,1,:) = jz(:,:,1,:) + real(jz_c)
+ jz(:,:,2,:) = jz(:,:,2,:) + aimag(jz_c)
+
+ deallocate(jr_c)
+ deallocate(jt_c)
+ deallocate(jz_c)
+
+ else
+#endif
+ CALL WRPX_PXR_CURRENT_DEPOSITION( &
jx,jx_nguards,jx_nvalid, &
jy,jy_nguards,jy_nvalid, &
jz,jz_nguards,jz_nvalid, &
np,xp,yp,zp,uxp,uyp,uzp,gaminv,w,q, &
xmin,zmin,dt,dx,dz,nox,noz,pxr_l_nodal, &
lvect,current_depo_algo)
+#ifdef WARPX_RZ
+ endif
+#endif
#endif
end subroutine warpx_current_deposition
diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp
index 541999789..5a5c190c6 100644
--- a/Source/Initialization/PlasmaInjector.cpp
+++ b/Source/Initialization/PlasmaInjector.cpp
@@ -143,9 +143,11 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name)
parseDensity(pp);
parseMomentum(pp);
} else if (part_pos_s == "nuniformpercell") {
- num_particles_per_cell_each_dim.resize(3);
+ // Note that for RZ, three numbers are expected, r, theta, and z.
+ // For 2D, only two are expected. The third is overwritten with 1.
+ num_particles_per_cell_each_dim.assign(3, 1);
pp.getarr("num_particles_per_cell_each_dim", num_particles_per_cell_each_dim);
-#if ( AMREX_SPACEDIM == 2 )
+#if ( AMREX_SPACEDIM == 2 ) && !defined(WARPX_RZ)
num_particles_per_cell_each_dim[2] = 1;
#endif
// Construct InjectorPosition from InjectorPositionRegular.
diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp
index 9d85783b0..0ca1e8a5d 100644
--- a/Source/Parallelization/WarpXComm.cpp
+++ b/Source/Parallelization/WarpXComm.cpp
@@ -59,24 +59,24 @@ WarpX::UpdateAuxilaryData ()
// B field
{
- MultiFab dBx(Bfield_cp[lev][0]->boxArray(), dm, 1, ng);
- MultiFab dBy(Bfield_cp[lev][1]->boxArray(), dm, 1, ng);
- MultiFab dBz(Bfield_cp[lev][2]->boxArray(), dm, 1, ng);
+ MultiFab dBx(Bfield_cp[lev][0]->boxArray(), dm, Bfield_cp[lev][0]->nComp(), ng);
+ MultiFab dBy(Bfield_cp[lev][1]->boxArray(), dm, Bfield_cp[lev][1]->nComp(), ng);
+ MultiFab dBz(Bfield_cp[lev][2]->boxArray(), dm, Bfield_cp[lev][2]->nComp(), ng);
dBx.setVal(0.0);
dBy.setVal(0.0);
dBz.setVal(0.0);
- dBx.ParallelCopy(*Bfield_aux[lev-1][0], 0, 0, 1, ng, ng, crse_period);
- dBy.ParallelCopy(*Bfield_aux[lev-1][1], 0, 0, 1, ng, ng, crse_period);
- dBz.ParallelCopy(*Bfield_aux[lev-1][2], 0, 0, 1, ng, ng, crse_period);
+ dBx.ParallelCopy(*Bfield_aux[lev-1][0], 0, 0, Bfield_aux[lev-1][0]->nComp(), ng, ng, crse_period);
+ dBy.ParallelCopy(*Bfield_aux[lev-1][1], 0, 0, Bfield_aux[lev-1][1]->nComp(), ng, ng, crse_period);
+ dBz.ParallelCopy(*Bfield_aux[lev-1][2], 0, 0, Bfield_aux[lev-1][2]->nComp(), ng, ng, crse_period);
if (Bfield_cax[lev][0])
{
- MultiFab::Copy(*Bfield_cax[lev][0], dBx, 0, 0, 1, ng);
- MultiFab::Copy(*Bfield_cax[lev][1], dBy, 0, 0, 1, ng);
- MultiFab::Copy(*Bfield_cax[lev][2], dBz, 0, 0, 1, ng);
+ MultiFab::Copy(*Bfield_cax[lev][0], dBx, 0, 0, Bfield_cax[lev][0]->nComp(), ng);
+ MultiFab::Copy(*Bfield_cax[lev][1], dBy, 0, 0, Bfield_cax[lev][1]->nComp(), ng);
+ MultiFab::Copy(*Bfield_cax[lev][2], dBz, 0, 0, Bfield_cax[lev][2]->nComp(), ng);
}
- MultiFab::Subtract(dBx, *Bfield_cp[lev][0], 0, 0, 1, ng);
- MultiFab::Subtract(dBy, *Bfield_cp[lev][1], 0, 0, 1, ng);
- MultiFab::Subtract(dBz, *Bfield_cp[lev][2], 0, 0, 1, ng);
+ MultiFab::Subtract(dBx, *Bfield_cp[lev][0], 0, 0, Bfield_cp[lev][0]->nComp(), ng);
+ MultiFab::Subtract(dBy, *Bfield_cp[lev][1], 0, 0, Bfield_cp[lev][1]->nComp(), ng);
+ MultiFab::Subtract(dBz, *Bfield_cp[lev][2], 0, 0, Bfield_cp[lev][2]->nComp(), ng);
const Real* dx = Geom(lev-1).CellSize();
const int refinement_ratio = refRatio(lev-1)[0];
@@ -134,24 +134,24 @@ WarpX::UpdateAuxilaryData ()
// E field
{
- MultiFab dEx(Efield_cp[lev][0]->boxArray(), dm, 1, ng);
- MultiFab dEy(Efield_cp[lev][1]->boxArray(), dm, 1, ng);
- MultiFab dEz(Efield_cp[lev][2]->boxArray(), dm, 1, ng);
+ MultiFab dEx(Efield_cp[lev][0]->boxArray(), dm, Efield_cp[lev][0]->nComp(), ng);
+ MultiFab dEy(Efield_cp[lev][1]->boxArray(), dm, Efield_cp[lev][1]->nComp(), ng);
+ MultiFab dEz(Efield_cp[lev][2]->boxArray(), dm, Efield_cp[lev][2]->nComp(), ng);
dEx.setVal(0.0);
dEy.setVal(0.0);
dEz.setVal(0.0);
- dEx.ParallelCopy(*Efield_aux[lev-1][0], 0, 0, 1, ng, ng, crse_period);
- dEy.ParallelCopy(*Efield_aux[lev-1][1], 0, 0, 1, ng, ng, crse_period);
- dEz.ParallelCopy(*Efield_aux[lev-1][2], 0, 0, 1, ng, ng, crse_period);
+ dEx.ParallelCopy(*Efield_aux[lev-1][0], 0, 0, Efield_aux[lev-1][0]->nComp(), ng, ng, crse_period);
+ dEy.ParallelCopy(*Efield_aux[lev-1][1], 0, 0, Efield_aux[lev-1][1]->nComp(), ng, ng, crse_period);
+ dEz.ParallelCopy(*Efield_aux[lev-1][2], 0, 0, Efield_aux[lev-1][2]->nComp(), ng, ng, crse_period);
if (Efield_cax[lev][0])
{
- MultiFab::Copy(*Efield_cax[lev][0], dEx, 0, 0, 1, ng);
- MultiFab::Copy(*Efield_cax[lev][1], dEy, 0, 0, 1, ng);
- MultiFab::Copy(*Efield_cax[lev][2], dEz, 0, 0, 1, ng);
+ MultiFab::Copy(*Efield_cax[lev][0], dEx, 0, 0, Efield_cax[lev][0]->nComp(), ng);
+ MultiFab::Copy(*Efield_cax[lev][1], dEy, 0, 0, Efield_cax[lev][1]->nComp(), ng);
+ MultiFab::Copy(*Efield_cax[lev][2], dEz, 0, 0, Efield_cax[lev][2]->nComp(), ng);
}
- MultiFab::Subtract(dEx, *Efield_cp[lev][0], 0, 0, 1, ng);
- MultiFab::Subtract(dEy, *Efield_cp[lev][1], 0, 0, 1, ng);
- MultiFab::Subtract(dEz, *Efield_cp[lev][2], 0, 0, 1, ng);
+ MultiFab::Subtract(dEx, *Efield_cp[lev][0], 0, 0, Efield_cp[lev][0]->nComp(), ng);
+ MultiFab::Subtract(dEy, *Efield_cp[lev][1], 0, 0, Efield_cp[lev][1]->nComp(), ng);
+ MultiFab::Subtract(dEz, *Efield_cp[lev][2], 0, 0, Efield_cp[lev][2]->nComp(), ng);
const int refinement_ratio = refRatio(lev-1)[0];
#ifdef _OPEMP
@@ -199,8 +199,8 @@ WarpX::UpdateAuxilaryData ()
FArrayBox& aux = (*Efield_aux[lev][idim])[mfi];
FArrayBox& fp = (*Efield_fp[lev][idim])[mfi];
const Box& bx = aux.box();
- aux.copy(fp, bx, 0, bx, 0, 1);
- aux.plus(efab[idim], bx, bx, 0, 0, 1);
+ aux.copy(fp, bx, 0, bx, 0, Efield_fp[lev][idim]->nComp());
+ aux.plus(efab[idim], bx, bx, 0, 0, Efield_fp[lev][idim]->nComp());
}
}
}
@@ -409,7 +409,7 @@ WarpX::SyncCurrent (const std::array<const amrex::MultiFab*,3>& fine,
ffab.resize(fbx);
fbx &= (*fine[idim])[mfi].box();
ffab.setVal(0.0);
- ffab.copy((*fine[idim])[mfi], fbx, 0, fbx, 0, 1);
+ ffab.copy((*fine[idim])[mfi], fbx, 0, fbx, 0, fine[idim]->nComp());
WRPX_SYNC_CURRENT(bx.loVect(), bx.hiVect(),
BL_TO_FORTRAN_ANYD((*crse[idim])[mfi]),
BL_TO_FORTRAN_ANYD(ffab),
@@ -505,11 +505,11 @@ WarpX::ApplyFilterandSumBoundaryJ (int lev, PatchType patch_type)
if (use_filter) {
IntVect ng = j[idim]->nGrowVect();
ng += bilinear_filter.stencil_length_each_dir-1;
- MultiFab jf(j[idim]->boxArray(), j[idim]->DistributionMap(), 1, ng);
+ MultiFab jf(j[idim]->boxArray(), j[idim]->DistributionMap(), j[idim]->nComp(), ng);
bilinear_filter.ApplyStencil(jf, *j[idim]);
- WarpXSumGuardCells(*(j[idim]), jf, period);
+ WarpXSumGuardCells(*(j[idim]), jf, period, 0, (j[idim])->nComp());
} else {
- WarpXSumGuardCells(*(j[idim]), period);
+ WarpXSumGuardCells(*(j[idim]), period, 0, (j[idim])->nComp());
}
}
}
@@ -539,7 +539,7 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev)
const auto& period = Geom(lev).periodicity();
for (int idim = 0; idim < 3; ++idim) {
MultiFab mf(current_fp[lev][idim]->boxArray(),
- current_fp[lev][idim]->DistributionMap(), 1, 0);
+ current_fp[lev][idim]->DistributionMap(), current_fp[lev][idim]->nComp(), 0);
mf.setVal(0.0);
if (use_filter && current_buf[lev+1][idim])
{
@@ -547,18 +547,18 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev)
IntVect ng = current_cp[lev+1][idim]->nGrowVect();
ng += bilinear_filter.stencil_length_each_dir-1;
MultiFab jfc(current_cp[lev+1][idim]->boxArray(),
- current_cp[lev+1][idim]->DistributionMap(), 1, ng);
+ current_cp[lev+1][idim]->DistributionMap(), current_cp[lev+1][idim]->nComp(), ng);
bilinear_filter.ApplyStencil(jfc, *current_cp[lev+1][idim]);
// buffer patch of fine level
MultiFab jfb(current_buf[lev+1][idim]->boxArray(),
- current_buf[lev+1][idim]->DistributionMap(), 1, ng);
+ current_buf[lev+1][idim]->DistributionMap(), current_buf[lev+1][idim]->nComp(), ng);
bilinear_filter.ApplyStencil(jfb, *current_buf[lev+1][idim]);
- MultiFab::Add(jfb, jfc, 0, 0, 1, ng);
- mf.ParallelAdd(jfb, 0, 0, 1, ng, IntVect::TheZeroVector(), period);
+ MultiFab::Add(jfb, jfc, 0, 0, current_buf[lev+1][idim]->nComp(), ng);
+ mf.ParallelAdd(jfb, 0, 0, current_buf[lev+1][idim]->nComp(), ng, IntVect::TheZeroVector(), period);
- WarpXSumGuardCells(*current_cp[lev+1][idim], jfc, period);
+ WarpXSumGuardCells(*current_cp[lev+1][idim], jfc, period, 0, current_cp[lev+1][idim]->nComp());
}
else if (use_filter) // but no buffer
{
@@ -566,29 +566,29 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev)
IntVect ng = current_cp[lev+1][idim]->nGrowVect();
ng += bilinear_filter.stencil_length_each_dir-1;
MultiFab jf(current_cp[lev+1][idim]->boxArray(),
- current_cp[lev+1][idim]->DistributionMap(), 1, ng);
+ current_cp[lev+1][idim]->DistributionMap(), current_cp[lev+1][idim]->nComp(), ng);
bilinear_filter.ApplyStencil(jf, *current_cp[lev+1][idim]);
- mf.ParallelAdd(jf, 0, 0, 1, ng, IntVect::TheZeroVector(), period);
- WarpXSumGuardCells(*current_cp[lev+1][idim], jf, period);
+ mf.ParallelAdd(jf, 0, 0, current_cp[lev+1][idim]->nComp(), ng, IntVect::TheZeroVector(), period);
+ WarpXSumGuardCells(*current_cp[lev+1][idim], jf, period, 0, current_cp[lev+1][idim]->nComp());
}
else if (current_buf[lev+1][idim]) // but no filter
{
MultiFab::Copy(*current_buf[lev+1][idim],
- *current_cp [lev+1][idim], 0, 0, 1,
+ *current_cp [lev+1][idim], 0, 0, current_buf[lev+1][idim]->nComp(),
current_cp[lev+1][idim]->nGrow());
- mf.ParallelAdd(*current_buf[lev+1][idim], 0, 0, 1,
+ mf.ParallelAdd(*current_buf[lev+1][idim], 0, 0, current_buf[lev+1][idim]->nComp(),
current_buf[lev+1][idim]->nGrowVect(), IntVect::TheZeroVector(),
period);
- WarpXSumGuardCells(*(current_cp[lev+1][idim]), period);
+ WarpXSumGuardCells(*(current_cp[lev+1][idim]), period, 0, current_cp[lev+1][idim]->nComp());
}
else // no filter, no buffer
{
- mf.ParallelAdd(*current_cp[lev+1][idim], 0, 0, 1,
+ mf.ParallelAdd(*current_cp[lev+1][idim], 0, 0, current_cp[lev+1][idim]->nComp(),
current_cp[lev+1][idim]->nGrowVect(), IntVect::TheZeroVector(),
period);
- WarpXSumGuardCells(*(current_cp[lev+1][idim]), period);
+ WarpXSumGuardCells(*(current_cp[lev+1][idim]), period, 0, current_cp[lev+1][idim]->nComp());
}
- MultiFab::Add(*current_fp[lev][idim], mf, 0, 0, 1, 0);
+ MultiFab::Add(*current_fp[lev][idim], mf, 0, 0, current_fp[lev+1][idim]->nComp(), 0);
}
NodalSyncJ(lev+1, PatchType::coarse);
}
diff --git a/Source/Parallelization/WarpXRegrid.cpp b/Source/Parallelization/WarpXRegrid.cpp
index 8d7873041..eb119d4a2 100644
--- a/Source/Parallelization/WarpXRegrid.cpp
+++ b/Source/Parallelization/WarpXRegrid.cpp
@@ -46,21 +46,21 @@ WarpX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa
{
const IntVect& ng = Bfield_fp[lev][idim]->nGrowVect();
auto pmf = std::unique_ptr<MultiFab>(new MultiFab(Bfield_fp[lev][idim]->boxArray(),
- dm, 1, ng));
- pmf->Redistribute(*Bfield_fp[lev][idim], 0, 0, 1, ng);
+ dm, Bfield_fp[lev][idim]->nComp(), ng));
+ pmf->Redistribute(*Bfield_fp[lev][idim], 0, 0, Bfield_fp[lev][idim]->nComp(), ng);
Bfield_fp[lev][idim] = std::move(pmf);
}
{
const IntVect& ng = Efield_fp[lev][idim]->nGrowVect();
auto pmf = std::unique_ptr<MultiFab>(new MultiFab(Efield_fp[lev][idim]->boxArray(),
- dm, 1, ng));
- pmf->Redistribute(*Efield_fp[lev][idim], 0, 0, 1, ng);
+ dm, Efield_fp[lev][idim]->nComp(), ng));
+ pmf->Redistribute(*Efield_fp[lev][idim], 0, 0, Efield_fp[lev][idim]->nComp(), ng);
Efield_fp[lev][idim] = std::move(pmf);
}
{
const IntVect& ng = current_fp[lev][idim]->nGrowVect();
auto pmf = std::unique_ptr<MultiFab>(new MultiFab(current_fp[lev][idim]->boxArray(),
- dm, 1, ng));
+ dm, current_fp[lev][idim]->nComp(), ng));
current_fp[lev][idim] = std::move(pmf);
current_fp_owner_masks[lev][idim] = std::move(current_fp[lev][idim]->OwnerMask(period));
}
@@ -68,7 +68,7 @@ WarpX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa
{
const IntVect& ng = current_store[lev][idim]->nGrowVect();
auto pmf = std::unique_ptr<MultiFab>(new MultiFab(current_store[lev][idim]->boxArray(),
- dm, 1, ng));
+ dm, current_store[lev][idim]->nComp(), ng));
// no need to redistribute
current_store[lev][idim] = std::move(pmf);
}
@@ -77,8 +77,8 @@ WarpX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa
if (F_fp[lev] != nullptr) {
const IntVect& ng = F_fp[lev]->nGrowVect();
auto pmf = std::unique_ptr<MultiFab>(new MultiFab(F_fp[lev]->boxArray(),
- dm, 1, ng));
- pmf->Redistribute(*F_fp[lev], 0, 0, 1, ng);
+ dm, F_fp[lev]->nComp(), ng));
+ pmf->Redistribute(*F_fp[lev], 0, 0, F_fp[lev]->nComp(), ng);
F_fp[lev] = std::move(pmf);
}
@@ -96,8 +96,8 @@ WarpX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa
if (lev == 0)
{
for (int idim = 0; idim < 3; ++idim) {
- Bfield_aux[lev][idim].reset(new MultiFab(*Bfield_fp[lev][idim], amrex::make_alias, 0, 1));
- Efield_aux[lev][idim].reset(new MultiFab(*Efield_fp[lev][idim], amrex::make_alias, 0, 1));
+ Bfield_aux[lev][idim].reset(new MultiFab(*Bfield_fp[lev][idim], amrex::make_alias, 0, Bfield_aux[lev][idim]->nComp()));
+ Efield_aux[lev][idim].reset(new MultiFab(*Efield_fp[lev][idim], amrex::make_alias, 0, Efield_aux[lev][idim]->nComp()));
}
} else {
for (int idim=0; idim < 3; ++idim)
@@ -105,15 +105,15 @@ WarpX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa
{
const IntVect& ng = Bfield_aux[lev][idim]->nGrowVect();
auto pmf = std::unique_ptr<MultiFab>(new MultiFab(Bfield_aux[lev][idim]->boxArray(),
- dm, 1, ng));
- // pmf->Redistribute(*Bfield_aux[lev][idim], 0, 0, 1, ng);
+ dm, Bfield_aux[lev][idim]->nComp(), ng));
+ // pmf->Redistribute(*Bfield_aux[lev][idim], 0, 0, Bfield_aux[lev][idim]->nComp(), ng);
Bfield_aux[lev][idim] = std::move(pmf);
}
{
const IntVect& ng = Efield_aux[lev][idim]->nGrowVect();
auto pmf = std::unique_ptr<MultiFab>(new MultiFab(Efield_aux[lev][idim]->boxArray(),
- dm, 1, ng));
- // pmf->Redistribute(*Efield_aux[lev][idim], 0, 0, 1, ng);
+ dm, Efield_aux[lev][idim]->nComp(), ng));
+ // pmf->Redistribute(*Efield_aux[lev][idim], 0, 0, Efield_aux[lev][idim]->nComp(), ng);
Efield_aux[lev][idim] = std::move(pmf);
}
}
@@ -127,21 +127,21 @@ WarpX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa
{
const IntVect& ng = Bfield_cp[lev][idim]->nGrowVect();
auto pmf = std::unique_ptr<MultiFab>(new MultiFab(Bfield_cp[lev][idim]->boxArray(),
- dm, 1, ng));
- pmf->Redistribute(*Bfield_cp[lev][idim], 0, 0, 1, ng);
+ dm, Bfield_cp[lev][idim]->nComp(), ng));
+ pmf->Redistribute(*Bfield_cp[lev][idim], 0, 0, Bfield_cp[lev][idim]->nComp(), ng);
Bfield_cp[lev][idim] = std::move(pmf);
}
{
const IntVect& ng = Efield_cp[lev][idim]->nGrowVect();
auto pmf = std::unique_ptr<MultiFab>(new MultiFab(Efield_cp[lev][idim]->boxArray(),
- dm, 1, ng));
- pmf->Redistribute(*Efield_cp[lev][idim], 0, 0, 1, ng);
+ dm, Efield_cp[lev][idim]->nComp(), ng));
+ pmf->Redistribute(*Efield_cp[lev][idim], 0, 0, Efield_cp[lev][idim]->nComp(), ng);
Efield_cp[lev][idim] = std::move(pmf);
}
{
const IntVect& ng = current_cp[lev][idim]->nGrowVect();
auto pmf = std::unique_ptr<MultiFab>( new MultiFab(current_cp[lev][idim]->boxArray(),
- dm, 1, ng));
+ dm, current_cp[lev][idim]->nComp(), ng));
current_cp[lev][idim] = std::move(pmf);
current_cp_owner_masks[lev][idim] = std::move(
current_cp[lev][idim]->OwnerMask(cperiod));
@@ -151,8 +151,8 @@ WarpX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa
if (F_cp[lev] != nullptr) {
const IntVect& ng = F_cp[lev]->nGrowVect();
auto pmf = std::unique_ptr<MultiFab>(new MultiFab(F_cp[lev]->boxArray(),
- dm, 1, ng));
- pmf->Redistribute(*F_cp[lev], 0, 0, 1, ng);
+ dm, F_cp[lev]->nComp(), ng));
+ pmf->Redistribute(*F_cp[lev], 0, 0, F_cp[lev]->nComp(), ng);
F_cp[lev] = std::move(pmf);
}
@@ -173,24 +173,24 @@ WarpX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa
{
const IntVect& ng = Bfield_cax[lev][idim]->nGrowVect();
auto pmf = std::unique_ptr<MultiFab>(new MultiFab(Bfield_cax[lev][idim]->boxArray(),
- dm, 1, ng));
- // pmf->ParallelCopy(*Bfield_cax[lev][idim], 0, 0, 1, ng, ng);
+ dm, Bfield_cax[lev][idim]->nComp(), ng));
+ // pmf->ParallelCopy(*Bfield_cax[lev][idim], 0, 0, Bfield_cax[lev][idim]->nComp(), ng, ng);
Bfield_cax[lev][idim] = std::move(pmf);
}
if (Efield_cax[lev][idim])
{
const IntVect& ng = Efield_cax[lev][idim]->nGrowVect();
auto pmf = std::unique_ptr<MultiFab>(new MultiFab(Efield_cax[lev][idim]->boxArray(),
- dm, 1, ng));
- // pmf->ParallelCopy(*Efield_cax[lev][idim], 0, 0, 1, ng, ng);
+ dm, Efield_cax[lev][idim]->nComp(), ng));
+ // pmf->ParallelCopy(*Efield_cax[lev][idim], 0, 0, Efield_cax[lev][idim]->nComp(), ng, ng);
Efield_cax[lev][idim] = std::move(pmf);
}
if (current_buf[lev][idim])
{
const IntVect& ng = current_buf[lev][idim]->nGrowVect();
auto pmf = std::unique_ptr<MultiFab>(new MultiFab(current_buf[lev][idim]->boxArray(),
- dm, 1, ng));
- // pmf->ParallelCopy(*current_buf[lev][idim], 0, 0, 1, ng, ng);
+ dm, current_buf[lev][idim]->nComp(), ng));
+ // pmf->ParallelCopy(*current_buf[lev][idim], 0, 0, current_buf[lev][idim]->nComp(), ng, ng);
current_buf[lev][idim] = std::move(pmf);
}
}
@@ -198,24 +198,24 @@ WarpX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa
{
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);
+ dm, charge_buf[lev]->nComp(), ng));
+ // pmf->ParallelCopy(*charge_buf[lev][idim], 0, 0, charge_buf[lev]->nComp(), ng, ng);
charge_buf[lev] = std::move(pmf);
}
if (current_buffer_masks[lev])
{
const IntVect& ng = current_buffer_masks[lev]->nGrowVect();
auto pmf = std::unique_ptr<iMultiFab>(new iMultiFab(current_buffer_masks[lev]->boxArray(),
- dm, 1, ng));
- // pmf->ParallelCopy(*current_buffer_masks[lev], 0, 0, 1, ng, ng);
+ dm, current_buffer_masks[lev]->nComp(), ng));
+ // pmf->ParallelCopy(*current_buffer_masks[lev], 0, 0, current_buffer_masks[lev]->nComp(), ng, ng);
current_buffer_masks[lev] = std::move(pmf);
}
if (gather_buffer_masks[lev])
{
const IntVect& ng = gather_buffer_masks[lev]->nGrowVect();
auto pmf = std::unique_ptr<iMultiFab>(new iMultiFab(gather_buffer_masks[lev]->boxArray(),
- dm, 1, ng));
- // pmf->ParallelCopy(*gather_buffer_masks[lev], 0, 0, 1, ng, ng);
+ dm, gather_buffer_masks[lev]->nComp(), ng));
+ // pmf->ParallelCopy(*gather_buffer_masks[lev], 0, 0, gather_buffer_masks[lev]->nComp(), ng, ng);
gather_buffer_masks[lev] = std::move(pmf);
}
}
diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp
index 982e04e39..6a1294ddf 100644
--- a/Source/Particles/MultiParticleContainer.cpp
+++ b/Source/Particles/MultiParticleContainer.cpp
@@ -274,7 +274,7 @@ MultiParticleContainer::GetChargeDensity (int lev, bool local)
std::unique_ptr<MultiFab> rho = allcontainers[0]->GetChargeDensity(lev, true);
for (unsigned i = 1, n = allcontainers.size(); i < n; ++i) {
std::unique_ptr<MultiFab> rhoi = allcontainers[i]->GetChargeDensity(lev, true);
- MultiFab::Add(*rho, *rhoi, 0, 0, 1, rho->nGrow());
+ MultiFab::Add(*rho, *rhoi, 0, 0, rho->nComp(), rho->nGrow());
}
if (!local) {
const Geometry& gm = allcontainers[0]->Geom(lev);
diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H
index b80619733..ebd4eac2b 100644
--- a/Source/Particles/PhysicalParticleContainer.H
+++ b/Source/Particles/PhysicalParticleContainer.H
@@ -53,6 +53,7 @@ public:
amrex::FArrayBox const * byfab,
amrex::FArrayBox const * bzfab,
const int ngE, const int e_is_nodal,
+ const long n_rz_azimuthal_modes,
const long offset,
const long np_to_gather,
int thread_num,
diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
index d10390204..c44a9a8b7 100644
--- a/Source/Particles/PhysicalParticleContainer.cpp
+++ b/Source/Particles/PhysicalParticleContainer.cpp
@@ -899,7 +899,8 @@ PhysicalParticleContainer::FieldGather (int lev,
int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal();
FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
&exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab,
- Ex.nGrow(), e_is_nodal, 0, np, thread_num, lev, lev);
+ Ex.nGrow(), e_is_nodal, WarpX::n_rz_azimuthal_modes,
+ 0, np, thread_num, lev, lev);
if (cost) {
const Box& tbx = pti.tilebox();
@@ -1175,7 +1176,8 @@ PhysicalParticleContainer::Evolve (int lev,
BL_PROFILE_VAR_START(blp_pxr_fg);
FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
exfab, eyfab, ezfab, bxfab, byfab, bzfab,
- Ex.nGrow(), e_is_nodal, 0, np_gather, thread_num, lev, lev);
+ Ex.nGrow(), e_is_nodal, WarpX::n_rz_azimuthal_modes,
+ 0, np_gather, thread_num, lev, lev);
if (np_gather < np)
{
@@ -1248,6 +1250,7 @@ PhysicalParticleContainer::Evolve (int lev,
cexfab, ceyfab, cezfab,
cbxfab, cbyfab, cbzfab,
cEx->nGrow(), e_is_nodal,
+ WarpX::n_rz_azimuthal_modes,
nfine_gather, np-nfine_gather,
thread_num, lev, lev-1);
}
@@ -1591,7 +1594,8 @@ PhysicalParticleContainer::PushP (int lev, Real dt,
int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal();
FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
&exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab,
- Ex.nGrow(), e_is_nodal, 0, np, thread_num, lev, lev);
+ Ex.nGrow(), e_is_nodal, WarpX::n_rz_azimuthal_modes,
+ 0, np, thread_num, lev, lev);
// This wraps the momentum advance so that inheritors can modify the call.
// Extract pointers to the different particle quantities
@@ -1831,6 +1835,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti,
const int ngE, const int e_is_nodal,
const long offset,
const long np_to_gather,
+ const long n_rz_azimuthal_modes,
int thread_num,
int lev,
int gather_lev)
diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp
index 36cb9d224..695955e15 100644
--- a/Source/Particles/RigidInjectedParticleContainer.cpp
+++ b/Source/Particles/RigidInjectedParticleContainer.cpp
@@ -426,7 +426,8 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt,
int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal();
FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
&exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab,
- Ex.nGrow(), e_is_nodal, 0, np, thread_num, lev, lev);
+ Ex.nGrow(), e_is_nodal, WarpX::n_rz_azimuthal_modes,
+ 0, np, thread_num, lev, lev);
// Save the position and momenta, making copies
auto uxp_save = uxp;
diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H
index ac5b47ada..ee263c27b 100644
--- a/Source/Particles/WarpXParticleContainer.H
+++ b/Source/Particles/WarpXParticleContainer.H
@@ -42,6 +42,9 @@ namespace ParticleStringNames
{"Bx", PIdx::Bx },
{"By", PIdx::By },
{"Bz", PIdx::Bz }
+#ifdef WARPX_RZ
+ ,{"theta", PIdx::theta}
+#endif
};
}
diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp
index befa5cfed..0ff80941d 100644
--- a/Source/Particles/WarpXParticleContainer.cpp
+++ b/Source/Particles/WarpXParticleContainer.cpp
@@ -356,9 +356,9 @@ WarpXParticleContainer::DepositCurrentFortran(WarpXParIter& pti,
tby.grow(ngJ);
tbz.grow(ngJ);
- local_jx[thread_num].resize(tbx);
- local_jy[thread_num].resize(tby);
- local_jz[thread_num].resize(tbz);
+ local_jx[thread_num].resize(tbx, jx->nComp());
+ local_jy[thread_num].resize(tby, jy->nComp());
+ local_jz[thread_num].resize(tbz, jz->nComp());
Real* jx_ptr = local_jx[thread_num].dataPtr();
Real* jy_ptr = local_jy[thread_num].dataPtr();
@@ -381,6 +381,7 @@ WarpXParticleContainer::DepositCurrentFortran(WarpXParIter& pti,
jx_ptr, &ngJ, jxntot.getVect(),
jy_ptr, &ngJ, jyntot.getVect(),
jz_ptr, &ngJ, jzntot.getVect(),
+ &WarpX::n_rz_azimuthal_modes,
&np_to_depose,
m_xp[thread_num].dataPtr() + offset,
m_yp[thread_num].dataPtr() + offset,
@@ -401,9 +402,9 @@ WarpXParticleContainer::DepositCurrentFortran(WarpXParIter& pti,
BL_PROFILE_VAR_START(blp_accumulate);
// CPU, tiling: atomicAdd local_jx into jx
// (same for jx and jz)
- (*jx)[pti].atomicAdd(local_jx[thread_num], tbx, tbx, 0, 0, 1);
- (*jy)[pti].atomicAdd(local_jy[thread_num], tby, tby, 0, 0, 1);
- (*jz)[pti].atomicAdd(local_jz[thread_num], tbz, tbz, 0, 0, 1);
+ (*jx)[pti].atomicAdd(local_jx[thread_num], tbx, tbx, 0, 0, local_jx[thread_num].nComp());
+ (*jy)[pti].atomicAdd(local_jy[thread_num], tby, tby, 0, 0, local_jy[thread_num].nComp());
+ (*jz)[pti].atomicAdd(local_jz[thread_num], tbz, tbz, 0, 0, local_jz[thread_num].nComp());
BL_PROFILE_VAR_STOP(blp_accumulate);
#endif
}
@@ -585,23 +586,23 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp,
tilebox.grow(ngRho);
+ const int nc = (rhomf->nComp() == 1 ? 1 : rhomf->nComp()/2);
+
#ifdef AMREX_USE_GPU
// No tiling on GPU: rho_arr points to the full rho array.
- MultiFab rhoi(*rho, amrex::make_alias, icomp, 1);
+ MultiFab rhoi(*rho, amrex::make_alias, icomp*nc, nc);
Array4<Real> const& rho_arr = rhoi.array(pti);
#else
// Tiling is on: rho_arr points to local_rho[thread_num]
const Box tb = amrex::convert(tilebox, IntVect::TheUnitVector());
- local_rho[thread_num].resize(tb);
+ local_rho[thread_num].resize(tb, nc);
// local_rho[thread_num] is set to zero
local_rho[thread_num].setVal(0.0);
Array4<Real> const& rho_arr = local_rho[thread_num].array();
#endif
- // GPU, no tiling: deposit directly in rho
- // CPU, tiling: deposit into local_rho
Real* AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset;
Real* AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset;
@@ -629,7 +630,7 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp,
#ifndef AMREX_USE_GPU
BL_PROFILE_VAR_START(blp_accumulate);
- (*rho)[pti].atomicAdd(local_rho[thread_num], tb, tb, 0, icomp, 1);
+ (*rho)[pti].atomicAdd(local_rho[thread_num], tb, tb, 0, icomp*nc, nc);
BL_PROFILE_VAR_STOP(blp_accumulate);
#endif
@@ -683,7 +684,7 @@ WarpXParticleContainer::DepositCharge (Vector<std::unique_ptr<MultiFab> >& rho,
BoxArray coarsened_fine_BA = fine_BA;
coarsened_fine_BA.coarsen(m_gdb->refRatio(lev));
- MultiFab coarsened_fine_data(coarsened_fine_BA, fine_dm, 1, 0);
+ MultiFab coarsened_fine_data(coarsened_fine_BA, fine_dm, rho[lev+1]->nComp(), 0);
coarsened_fine_data.setVal(0.0);
IntVect ratio(AMREX_D_DECL(2, 2, 2)); // FIXME
@@ -712,7 +713,7 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local)
const int ng = WarpX::nox;
- auto rho = std::unique_ptr<MultiFab>(new MultiFab(nba,dm,1,ng));
+ auto rho = std::unique_ptr<MultiFab>(new MultiFab(nba,dm,WarpX::ncomps,ng));
rho->setVal(0.0);
#ifdef _OPENMP
diff --git a/Source/WarpX.H b/Source/WarpX.H
index 927cc1f32..cc44541ee 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -95,6 +95,10 @@ public:
static long noy;
static long noz;
+ // Number of modes for the RZ multimode version
+ static long n_rz_azimuthal_modes;
+ static long ncomps;
+
static bool use_fdtd_nci_corr;
static int l_lower_order_in_v;
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index a6a7a3399..2d6df625f 100644
--- a/Source/WarpX.cpp
+++ b/Source/WarpX.cpp
@@ -45,6 +45,9 @@ long WarpX::field_gathering_algo;
long WarpX::particle_pusher_algo;
int WarpX::maxwell_fdtd_solver_id;
+long WarpX::n_rz_azimuthal_modes = 1;
+long WarpX::ncomps = 1;
+
long WarpX::nox = 1;
long WarpX::noy = 1;
long WarpX::noz = 1;
@@ -499,6 +502,10 @@ WarpX::ReadParameters ()
// Use same shape factors in all directions, for gathering
l_lower_order_in_v = false;
}
+
+ // Only needs to be set with WARPX_RZ, otherwise defaults to 1.
+ pp.query("n_rz_azimuthal_modes", n_rz_azimuthal_modes);
+
}
{
@@ -768,20 +775,31 @@ void
WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm,
const IntVect& ngE, const IntVect& ngJ, const IntVect& ngRho, int ngF)
{
+
+#if defined WARPX_RZ
+ if (n_rz_azimuthal_modes > 1) {
+ // There is a real and imaginary component for each mode
+ ncomps = n_rz_azimuthal_modes*2;
+ } else {
+ // With only mode 0, only reals are used
+ ncomps = 1;
+ }
+#endif
+
//
// The fine patch
//
- Bfield_fp[lev][0].reset( new MultiFab(amrex::convert(ba,Bx_nodal_flag),dm,1,ngE));
- Bfield_fp[lev][1].reset( new MultiFab(amrex::convert(ba,By_nodal_flag),dm,1,ngE));
- Bfield_fp[lev][2].reset( new MultiFab(amrex::convert(ba,Bz_nodal_flag),dm,1,ngE));
+ Bfield_fp[lev][0].reset( new MultiFab(amrex::convert(ba,Bx_nodal_flag),dm,ncomps,ngE));
+ Bfield_fp[lev][1].reset( new MultiFab(amrex::convert(ba,By_nodal_flag),dm,ncomps,ngE));
+ Bfield_fp[lev][2].reset( new MultiFab(amrex::convert(ba,Bz_nodal_flag),dm,ncomps,ngE));
- Efield_fp[lev][0].reset( new MultiFab(amrex::convert(ba,Ex_nodal_flag),dm,1,ngE));
- Efield_fp[lev][1].reset( new MultiFab(amrex::convert(ba,Ey_nodal_flag),dm,1,ngE));
- Efield_fp[lev][2].reset( new MultiFab(amrex::convert(ba,Ez_nodal_flag),dm,1,ngE));
+ Efield_fp[lev][0].reset( new MultiFab(amrex::convert(ba,Ex_nodal_flag),dm,ncomps,ngE));
+ Efield_fp[lev][1].reset( new MultiFab(amrex::convert(ba,Ey_nodal_flag),dm,ncomps,ngE));
+ Efield_fp[lev][2].reset( new MultiFab(amrex::convert(ba,Ez_nodal_flag),dm,ncomps,ngE));
- current_fp[lev][0].reset( new MultiFab(amrex::convert(ba,jx_nodal_flag),dm,1,ngJ));
- current_fp[lev][1].reset( new MultiFab(amrex::convert(ba,jy_nodal_flag),dm,1,ngJ));
- current_fp[lev][2].reset( new MultiFab(amrex::convert(ba,jz_nodal_flag),dm,1,ngJ));
+ current_fp[lev][0].reset( new MultiFab(amrex::convert(ba,jx_nodal_flag),dm,ncomps,ngJ));
+ current_fp[lev][1].reset( new MultiFab(amrex::convert(ba,jy_nodal_flag),dm,ncomps,ngJ));
+ current_fp[lev][2].reset( new MultiFab(amrex::convert(ba,jz_nodal_flag),dm,ncomps,ngJ));
const auto& period = Geom(lev).periodicity();
current_fp_owner_masks[lev][0] = std::move(current_fp[lev][0]->OwnerMask(period));
@@ -790,25 +808,25 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
if (do_dive_cleaning || plot_rho)
{
- rho_fp[lev].reset(new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,2,ngRho));
+ rho_fp[lev].reset(new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,2*ncomps,ngRho));
rho_fp_owner_masks[lev] = std::move(rho_fp[lev]->OwnerMask(period));
}
if (do_subcycling == 1 && lev == 0)
{
- current_store[lev][0].reset( new MultiFab(amrex::convert(ba,jx_nodal_flag),dm,1,ngJ));
- current_store[lev][1].reset( new MultiFab(amrex::convert(ba,jy_nodal_flag),dm,1,ngJ));
- current_store[lev][2].reset( new MultiFab(amrex::convert(ba,jz_nodal_flag),dm,1,ngJ));
+ current_store[lev][0].reset( new MultiFab(amrex::convert(ba,jx_nodal_flag),dm,ncomps,ngJ));
+ current_store[lev][1].reset( new MultiFab(amrex::convert(ba,jy_nodal_flag),dm,ncomps,ngJ));
+ current_store[lev][2].reset( new MultiFab(amrex::convert(ba,jz_nodal_flag),dm,ncomps,ngJ));
}
if (do_dive_cleaning)
{
- F_fp[lev].reset (new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,1, ngF));
+ F_fp[lev].reset (new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,ncomps, ngF));
}
#ifdef WARPX_USE_PSATD
else
{
- rho_fp[lev].reset(new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,2,ngRho));
+ rho_fp[lev].reset(new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,2*ncomps,ngRho));
rho_fp_owner_masks[lev] = std::move(rho_fp[lev]->OwnerMask(period));
}
if (fft_hybrid_mpi_decomposition == false){
@@ -834,19 +852,19 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
if (lev == 0)
{
for (int idir = 0; idir < 3; ++idir) {
- Efield_aux[lev][idir].reset(new MultiFab(*Efield_fp[lev][idir], amrex::make_alias, 0, 1));
- Bfield_aux[lev][idir].reset(new MultiFab(*Bfield_fp[lev][idir], amrex::make_alias, 0, 1));
+ Efield_aux[lev][idir].reset(new MultiFab(*Efield_fp[lev][idir], amrex::make_alias, 0, ncomps));
+ Bfield_aux[lev][idir].reset(new MultiFab(*Bfield_fp[lev][idir], amrex::make_alias, 0, ncomps));
}
}
else
{
- Bfield_aux[lev][0].reset( new MultiFab(amrex::convert(ba,Bx_nodal_flag),dm,1,ngE));
- Bfield_aux[lev][1].reset( new MultiFab(amrex::convert(ba,By_nodal_flag),dm,1,ngE));
- Bfield_aux[lev][2].reset( new MultiFab(amrex::convert(ba,Bz_nodal_flag),dm,1,ngE));
+ Bfield_aux[lev][0].reset( new MultiFab(amrex::convert(ba,Bx_nodal_flag),dm,ncomps,ngE));
+ Bfield_aux[lev][1].reset( new MultiFab(amrex::convert(ba,By_nodal_flag),dm,ncomps,ngE));
+ Bfield_aux[lev][2].reset( new MultiFab(amrex::convert(ba,Bz_nodal_flag),dm,ncomps,ngE));
- Efield_aux[lev][0].reset( new MultiFab(amrex::convert(ba,Ex_nodal_flag),dm,1,ngE));
- Efield_aux[lev][1].reset( new MultiFab(amrex::convert(ba,Ey_nodal_flag),dm,1,ngE));
- Efield_aux[lev][2].reset( new MultiFab(amrex::convert(ba,Ez_nodal_flag),dm,1,ngE));
+ Efield_aux[lev][0].reset( new MultiFab(amrex::convert(ba,Ex_nodal_flag),dm,ncomps,ngE));
+ Efield_aux[lev][1].reset( new MultiFab(amrex::convert(ba,Ey_nodal_flag),dm,ncomps,ngE));
+ Efield_aux[lev][2].reset( new MultiFab(amrex::convert(ba,Ez_nodal_flag),dm,ncomps,ngE));
}
//
@@ -858,19 +876,19 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
cba.coarsen(refRatio(lev-1));
// Create the MultiFabs for B
- Bfield_cp[lev][0].reset( new MultiFab(amrex::convert(cba,Bx_nodal_flag),dm,1,ngE));
- Bfield_cp[lev][1].reset( new MultiFab(amrex::convert(cba,By_nodal_flag),dm,1,ngE));
- Bfield_cp[lev][2].reset( new MultiFab(amrex::convert(cba,Bz_nodal_flag),dm,1,ngE));
+ Bfield_cp[lev][0].reset( new MultiFab(amrex::convert(cba,Bx_nodal_flag),dm,ncomps,ngE));
+ Bfield_cp[lev][1].reset( new MultiFab(amrex::convert(cba,By_nodal_flag),dm,ncomps,ngE));
+ Bfield_cp[lev][2].reset( new MultiFab(amrex::convert(cba,Bz_nodal_flag),dm,ncomps,ngE));
// Create the MultiFabs for E
- Efield_cp[lev][0].reset( new MultiFab(amrex::convert(cba,Ex_nodal_flag),dm,1,ngE));
- Efield_cp[lev][1].reset( new MultiFab(amrex::convert(cba,Ey_nodal_flag),dm,1,ngE));
- Efield_cp[lev][2].reset( new MultiFab(amrex::convert(cba,Ez_nodal_flag),dm,1,ngE));
+ Efield_cp[lev][0].reset( new MultiFab(amrex::convert(cba,Ex_nodal_flag),dm,ncomps,ngE));
+ Efield_cp[lev][1].reset( new MultiFab(amrex::convert(cba,Ey_nodal_flag),dm,ncomps,ngE));
+ Efield_cp[lev][2].reset( new MultiFab(amrex::convert(cba,Ez_nodal_flag),dm,ncomps,ngE));
// Create the MultiFabs for the current
- current_cp[lev][0].reset( new MultiFab(amrex::convert(cba,jx_nodal_flag),dm,1,ngJ));
- current_cp[lev][1].reset( new MultiFab(amrex::convert(cba,jy_nodal_flag),dm,1,ngJ));
- current_cp[lev][2].reset( new MultiFab(amrex::convert(cba,jz_nodal_flag),dm,1,ngJ));
+ current_cp[lev][0].reset( new MultiFab(amrex::convert(cba,jx_nodal_flag),dm,ncomps,ngJ));
+ current_cp[lev][1].reset( new MultiFab(amrex::convert(cba,jy_nodal_flag),dm,ncomps,ngJ));
+ current_cp[lev][2].reset( new MultiFab(amrex::convert(cba,jz_nodal_flag),dm,ncomps,ngJ));
const auto& cperiod = Geom(lev-1).periodicity();
current_cp_owner_masks[lev][0] = std::move(current_cp[lev][0]->OwnerMask(cperiod));
@@ -878,17 +896,17 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
current_cp_owner_masks[lev][2] = std::move(current_cp[lev][2]->OwnerMask(cperiod));
if (do_dive_cleaning || plot_rho){
- rho_cp[lev].reset(new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,2,ngRho));
+ rho_cp[lev].reset(new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,2*ncomps,ngRho));
rho_cp_owner_masks[lev] = std::move(rho_cp[lev]->OwnerMask(cperiod));
}
if (do_dive_cleaning)
{
- F_cp[lev].reset (new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,1, ngF));
+ F_cp[lev].reset (new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,ncomps, ngF));
}
#ifdef WARPX_USE_PSATD
else
{
- rho_cp[lev].reset(new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,2,ngRho));
+ rho_cp[lev].reset(new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,2*ncomps,ngRho));
rho_cp_owner_masks[lev] = std::move(rho_cp[lev]->OwnerMask(cperiod));
}
#endif
@@ -904,28 +922,28 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
if (n_field_gather_buffer > 0) {
// Create the MultiFabs for B
- Bfield_cax[lev][0].reset( new MultiFab(amrex::convert(cba,Bx_nodal_flag),dm,1,ngE));
- Bfield_cax[lev][1].reset( new MultiFab(amrex::convert(cba,By_nodal_flag),dm,1,ngE));
- Bfield_cax[lev][2].reset( new MultiFab(amrex::convert(cba,Bz_nodal_flag),dm,1,ngE));
+ Bfield_cax[lev][0].reset( new MultiFab(amrex::convert(cba,Bx_nodal_flag),dm,ncomps,ngE));
+ Bfield_cax[lev][1].reset( new MultiFab(amrex::convert(cba,By_nodal_flag),dm,ncomps,ngE));
+ Bfield_cax[lev][2].reset( new MultiFab(amrex::convert(cba,Bz_nodal_flag),dm,ncomps,ngE));
// Create the MultiFabs for E
- Efield_cax[lev][0].reset( new MultiFab(amrex::convert(cba,Ex_nodal_flag),dm,1,ngE));
- Efield_cax[lev][1].reset( new MultiFab(amrex::convert(cba,Ey_nodal_flag),dm,1,ngE));
- Efield_cax[lev][2].reset( new MultiFab(amrex::convert(cba,Ez_nodal_flag),dm,1,ngE));
+ Efield_cax[lev][0].reset( new MultiFab(amrex::convert(cba,Ex_nodal_flag),dm,ncomps,ngE));
+ Efield_cax[lev][1].reset( new MultiFab(amrex::convert(cba,Ey_nodal_flag),dm,ncomps,ngE));
+ Efield_cax[lev][2].reset( new MultiFab(amrex::convert(cba,Ez_nodal_flag),dm,ncomps,ngE));
- gather_buffer_masks[lev].reset( new iMultiFab(ba, dm, 1, 1) );
+ gather_buffer_masks[lev].reset( new iMultiFab(ba, dm, ncomps, 1) );
// Gather buffer masks have 1 ghost cell, because of the fact
// that particles may move by more than one cell when using subcycling.
}
if (n_current_deposition_buffer > 0) {
- 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));
+ current_buf[lev][0].reset( new MultiFab(amrex::convert(cba,jx_nodal_flag),dm,ncomps,ngJ));
+ current_buf[lev][1].reset( new MultiFab(amrex::convert(cba,jy_nodal_flag),dm,ncomps,ngJ));
+ current_buf[lev][2].reset( new MultiFab(amrex::convert(cba,jz_nodal_flag),dm,ncomps,ngJ));
if (do_dive_cleaning || plot_rho) {
- charge_buf[lev].reset( new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,2,ngRho));
+ charge_buf[lev].reset( new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,2*ncomps,ngRho));
}
- current_buffer_masks[lev].reset( new iMultiFab(ba, dm, 1, 1) );
+ current_buffer_masks[lev].reset( new iMultiFab(ba, dm, ncomps, 1) );
// Current buffer masks have 1 ghost cell, because of the fact
// that particles may move by more than one cell when using subcycling.
}
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-11 17:44:17 +0000