diff options
Diffstat (limited to 'Source')
27 files changed, 716 insertions, 275 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 540a968a2..d4d46f1bd 100644 --- a/Source/Diagnostics/FieldIO.cpp +++ b/Source/Diagnostics/FieldIO.cpp @@ -211,6 +211,22 @@ WriteOpenPMDFields( const std::string& filename, } #endif // WARPX_USE_OPENPMD +#ifdef WARPX_DIM_RZ +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=1 ; 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()); + } +} +#endif void PackPlotDataPtrs (Vector<const MultiFab*>& pmf, @@ -235,6 +251,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 @@ -263,23 +280,72 @@ 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_DIM_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); +#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_DIM_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); +#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."); @@ -312,6 +378,15 @@ AverageAndPackScalarField( MultiFab& mf_avg, } } +/** \brief Add variable names to the list. + */ +void +AddToVarNames (Vector<std::string>& varnames, + std::string name, std::string suffix) { + auto coords = {"x", "y", "z"}; + for(auto coord:coords) varnames.push_back(name+coord+suffix); +} + /** \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. @@ -340,17 +415,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; @@ -444,11 +519,11 @@ 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); + AverageAndPackVectorField( mf_avg[lev], Efield_fp[lev], dmap[lev], dcomp, ngrow ); + if (lev == 0) AddToVarNames(varnames, "E", "_fp"); 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); + AverageAndPackVectorField( mf_avg[lev], Bfield_fp[lev], dmap[lev], dcomp, ngrow ); + if (lev == 0) AddToVarNames(varnames, "B", "_fp"); dcomp += 3; } @@ -462,10 +537,10 @@ 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); + if (lev == 0) AddToVarNames(varnames, "E", "_cp"); dcomp += 3; // now the magnetic field @@ -477,9 +552,9 @@ 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); + if (lev == 0) AddToVarNames(varnames, "B", "_cp"); dcomp += 3; } @@ -542,8 +617,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); } @@ -565,7 +640,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 9da0348d6..f79079b0e 100644 --- a/Source/Diagnostics/WarpXIO.cpp +++ b/Source/Diagnostics/WarpXIO.cpp @@ -421,13 +421,13 @@ WarpX::GetCellCenteredData() { int dcomp = 0; // first the electric field - AverageAndPackVectorField( *cc[lev], Efield_aux[lev], dcomp, ng ); + AverageAndPackVectorField( *cc[lev], Efield_aux[lev], dmap[lev], dcomp, ng ); dcomp += 3; // then the magnetic field - AverageAndPackVectorField( *cc[lev], Bfield_aux[lev], dcomp, ng ); + AverageAndPackVectorField( *cc[lev], Bfield_aux[lev], dmap[lev], dcomp, ng ); dcomp += 3; // then the current density - AverageAndPackVectorField( *cc[lev], current_fp[lev], dcomp, ng ); + AverageAndPackVectorField( *cc[lev], current_fp[lev], dmap[lev], dcomp, ng ); dcomp += 3; // then the charge density const std::unique_ptr<MultiFab>& charge_density = mypc->GetChargeDensity(lev); @@ -584,7 +584,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 594ea55ee..8f800665d 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -363,7 +363,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*ncomps); NodalSyncRho(fine_lev, PatchType::fine, 0, 2); EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]); @@ -390,7 +390,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, ncomps); EvolveB(fine_lev, PatchType::coarse, dt[fine_lev]); EvolveF(fine_lev, PatchType::coarse, dt[fine_lev], DtType::FirstHalf); @@ -417,7 +417,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, ncomps); NodalSyncRho(fine_lev, PatchType::fine, 0, 2); EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]); @@ -443,7 +443,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, ncomps, ncomps); EvolveE(fine_lev, PatchType::coarse, dt[fine_lev]); FillBoundaryE(fine_lev, PatchType::coarse); @@ -520,8 +520,22 @@ 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. + // 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/FieldSolver/SpectralSolver/Make.package b/Source/FieldSolver/SpectralSolver/Make.package index b0ee54095..205b363f2 100644 --- a/Source/FieldSolver/SpectralSolver/Make.package +++ b/Source/FieldSolver/SpectralSolver/Make.package @@ -1,4 +1,4 @@ -CEXE_headers += WarpX_Complex.H +CEXE_headers += WarpX_ComplexForFFT.H CEXE_headers += SpectralSolver.H CEXE_sources += SpectralSolver.cpp CEXE_headers += SpectralFieldData.H diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H index 6a2446981..01ca11083 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -1,7 +1,7 @@ #ifndef WARPX_SPECTRAL_FIELD_DATA_H_ #define WARPX_SPECTRAL_FIELD_DATA_H_ -#include <WarpX_Complex.H> +#include <WarpX_ComplexForFFT.H> #include <SpectralKSpace.H> #include <AMReX_MultiFab.H> diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H index ae7124b2e..a73356dca 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H @@ -1,7 +1,7 @@ #ifndef WARPX_SPECTRAL_K_SPACE_H_ #define WARPX_SPECTRAL_K_SPACE_H_ -#include <WarpX_Complex.H> +#include <WarpX_ComplexForFFT.H> #include <AMReX_BoxArray.H> #include <AMReX_LayoutData.H> diff --git a/Source/FieldSolver/SpectralSolver/WarpX_ComplexForFFT.H b/Source/FieldSolver/SpectralSolver/WarpX_ComplexForFFT.H new file mode 100644 index 000000000..43c7a1950 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/WarpX_ComplexForFFT.H @@ -0,0 +1,21 @@ +#ifndef WARPX_COMPLEXFORFFT_H_ +#define WARPX_COMPLEXFORFFT_H_ + +#include <WarpX_Complex.H> + +// Check the complex type on GPU/CPU +#ifdef AMREX_USE_GPU + +#include <cufft.h> +static_assert( sizeof(Complex) == sizeof(cuDoubleComplex), + "The complex types in WarpX and cuFFT do not match."); + +#else + +#include <fftw3.h> +static_assert( sizeof(Complex) == sizeof(fftw_complex), + "The complex types in WarpX and FFTW do not match."); + +#endif + +#endif //WARPX_COMPLEXFORFFT_H_ diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 11d598b98..b3d8e9d76 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -175,18 +175,19 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) warpx_push_bz_nodal(j,k,l,Bzfab,Exfab,Eyfab,dtsdx,dtsdy); }); } else if (WarpX::maxwell_fdtd_solver_id == 0) { + const long nmodes = n_rz_azimuthal_modes; amrex::ParallelFor(tbx, tby, tbz, [=] AMREX_GPU_DEVICE (int j, int k, int l) { - warpx_push_bx_yee(j,k,l,Bxfab,Eyfab,Ezfab,dtsdy,dtsdz); + warpx_push_bx_yee(j,k,l,Bxfab,Eyfab,Ezfab,dtsdx,dtsdy,dtsdz,dxinv,xmin,nmodes); }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { - warpx_push_by_yee(j,k,l,Byfab,Exfab,Ezfab,dtsdx,dtsdz); + warpx_push_by_yee(j,k,l,Byfab,Exfab,Ezfab,dtsdx,dtsdz,nmodes); }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { - warpx_push_bz_yee(j,k,l,Bzfab,Exfab,Eyfab,dtsdx,dtsdy,dxinv,xmin); + warpx_push_bz_yee(j,k,l,Bzfab,Exfab,Eyfab,dtsdx,dtsdy,dxinv,xmin,nmodes); }); } else if (WarpX::maxwell_fdtd_solver_id == 1) { Real betaxy, betaxz, betayx, betayz, betazx, betazy; @@ -366,18 +367,19 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) warpx_push_ez_nodal(j,k,l,Ezfab,Bxfab,Byfab,jzfab,mu_c2_dt,dtsdx_c2,dtsdy_c2); }); } else { + const long nmodes = n_rz_azimuthal_modes; amrex::ParallelFor(tex, tey, tez, [=] AMREX_GPU_DEVICE (int j, int k, int l) { - warpx_push_ex_yee(j,k,l,Exfab,Byfab,Bzfab,jxfab,mu_c2_dt,dtsdy_c2,dtsdz_c2); + warpx_push_ex_yee(j,k,l,Exfab,Byfab,Bzfab,jxfab,mu_c2_dt,dtsdx_c2,dtsdy_c2,dtsdz_c2,dxinv,xmin,nmodes); }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { - warpx_push_ey_yee(j,k,l,Eyfab,Bxfab,Bzfab,jyfab,mu_c2_dt,dtsdx_c2,dtsdz_c2,xmin); + warpx_push_ey_yee(j,k,l,Eyfab,Bxfab,Bzfab,jyfab,Exfab,mu_c2_dt,dtsdx_c2,dtsdz_c2,xmin,nmodes); }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { - warpx_push_ez_yee(j,k,l,Ezfab,Bxfab,Byfab,jzfab,mu_c2_dt,dtsdx_c2,dtsdy_c2,dxinv,xmin); + warpx_push_ez_yee(j,k,l,Ezfab,Bxfab,Byfab,jzfab,mu_c2_dt,dtsdx_c2,dtsdy_c2,dxinv,xmin,nmodes); }); } @@ -647,6 +649,8 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu const Real rmin = xyzmin[0]; const int irmin = lo.x; + const long nmodes = n_rz_azimuthal_modes; + // Rescale current in r-z mode since the inverse volume factor was not // included in the current deposition. amrex::ParallelFor(tbr, tbt, tbz, @@ -656,13 +660,29 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu // to the cells above the axis. // Note that Jr(i==0) is at 1/2 dr. if (rmin == 0. && 0 <= i && i < ngJ) { - Jr_arr(i,j,0) -= Jr_arr(-1-i,j,0); + Jr_arr(i,j,0,0) -= Jr_arr(-1-i,j,0,0); } // Apply the inverse volume scaling // Since Jr is not cell centered in r, no need for distinction // between on axis and off-axis factors const amrex::Real r = std::abs(rmin + (i - irmin + 0.5)*dr); - Jr_arr(i,j,0) /= (2.*MathConst::pi*r); + Jr_arr(i,j,0,0) /= (2.*MathConst::pi*r); + + for (int imode=1 ; imode < nmodes ; imode++) { + const Real ifact = ( (imode%2) == 0 ? +1. : -1.); + // Wrap the current density deposited in the guard cells around + // to the cells above the axis. + // Note that Jr(i==0) is at 1/2 dr. + if (rmin == 0. && 0 <= i && i < ngJ) { + Jr_arr(i,j,0,2*imode-1) -= ifact*Jr_arr(-1-i,j,0,2*imode-1); + Jr_arr(i,j,0,2*imode) -= ifact*Jr_arr(-1-i,j,0,2*imode); + } + // Apply the inverse volume scaling + // Since Jr is not cell centered in r, no need for distinction + // between on axis and off-axis factors + Jr_arr(i,j,0,2*imode-1) /= (2.*MathConst::pi*r); + Jr_arr(i,j,0,2*imode) /= (2.*MathConst::pi*r); + } }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { @@ -670,16 +690,37 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu // to the cells above the axis. // Jt is located on the boundary if (rmin == 0. && 0 < i && i <= ngJ) { - Jt_arr(i,j,0) += Jt_arr(-i,j,0); + Jt_arr(i,j,0,0) += Jt_arr(-i,j,0,0); } // Apply the inverse volume scaling // Jt is forced to zero on axis. const amrex::Real r = std::abs(rmin + (i - irmin)*dr); if (r == 0.) { - Jt_arr(i,j,0) = 0.; + Jt_arr(i,j,0,0) = 0.; } else { - Jt_arr(i,j,0) /= (2.*MathConst::pi*r); + Jt_arr(i,j,0,0) /= (2.*MathConst::pi*r); + } + + for (int imode=1 ; imode < nmodes ; imode++) { + const Real ifact = ( (imode%2) == 0 ? +1. : -1.); + // Wrap the current density deposited in the guard cells around + // to the cells above the axis. + // Jt is located on the boundary + if (rmin == 0. && 0 < i && i <= ngJ) { + Jt_arr(i,j,0,2*imode-1) += ifact*Jt_arr(-i,j,0,2*imode-1); + Jt_arr(i,j,0,2*imode) += ifact*Jt_arr(-i,j,0,2*imode); + } + + // Apply the inverse volume scaling + // Jt is forced to zero on axis. + if (r == 0.) { + Jt_arr(i,j,0,2*imode-1) = 0.; + Jt_arr(i,j,0,2*imode) = 0.; + } else { + Jt_arr(i,j,0,2*imode-1) /= (2.*MathConst::pi*r); + Jt_arr(i,j,0,2*imode) /= (2.*MathConst::pi*r); + } } }, [=] AMREX_GPU_DEVICE (int i, int j, int k) @@ -688,17 +729,39 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu // to the cells above the axis. // Jz is located on the boundary if (rmin == 0. && 0 < i && i <= ngJ) { - Jz_arr(i,j,0) += Jz_arr(-i,j,0); + Jz_arr(i,j,0,0) += Jz_arr(-i,j,0,0); } // Apply the inverse volume scaling const amrex::Real r = std::abs(rmin + (i - irmin)*dr); if (r == 0.) { // Verboncoeur JCP 164, 421-427 (2001) : corrected volume on axis - Jz_arr(i,j,0) /= (MathConst::pi*dr/3.); + Jz_arr(i,j,0,0) /= (MathConst::pi*dr/3.); } else { - Jz_arr(i,j,0) /= (2.*MathConst::pi*r); + Jz_arr(i,j,0,0) /= (2.*MathConst::pi*r); } + + for (int imode=1 ; imode < nmodes ; imode++) { + const Real ifact = ( (imode%2) == 0 ? +1. : -1.); + // Wrap the current density deposited in the guard cells around + // to the cells above the axis. + // Jz is located on the boundary + if (rmin == 0. && 0 < i && i <= ngJ) { + Jz_arr(i,j,0,2*imode-1) += ifact*Jz_arr(-i,j,0,2*imode-1); + Jz_arr(i,j,0,2*imode) += ifact*Jz_arr(-i,j,0,2*imode); + } + + // Apply the inverse volume scaling + if (r == 0.) { + // Verboncoeur JCP 164, 421-427 (2001) : corrected volume on axis + Jz_arr(i,j,0,2*imode-1) /= (MathConst::pi*dr/3.); + Jz_arr(i,j,0,2*imode) /= (MathConst::pi*dr/3.); + } else { + Jz_arr(i,j,0,2*imode-1) /= (2.*MathConst::pi*r); + Jz_arr(i,j,0,2*imode) /= (2.*MathConst::pi*r); + } + } + }); } } diff --git a/Source/FieldSolver/WarpX_FDTD.H b/Source/FieldSolver/WarpX_FDTD.H index 8f91036cc..232a84e8e 100644 --- a/Source/FieldSolver/WarpX_FDTD.H +++ b/Source/FieldSolver/WarpX_FDTD.H @@ -8,36 +8,75 @@ using namespace amrex; AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void warpx_push_bx_yee(int i, int j, int k, Array4<Real> const& Bx, Array4<Real const> const& Ey, Array4<Real const> const& Ez, - Real dtsdy, Real dtsdz) + Real dtsdx, Real dtsdy, Real dtsdz, Real dxinv, Real rmin, const long nmodes) { #if defined WARPX_DIM_3D Bx(i,j,k) += - dtsdy * (Ez(i,j+1,k ) - Ez(i,j,k)) + dtsdz * (Ey(i,j ,k+1) - Ey(i,j,k)); -#elif (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) - // Note that the 2D Cartesian and RZ versions are the same +#elif (defined WARPX_DIM_XZ) Bx(i,j,0) += + dtsdz * (Ey(i,j+1,0) - Ey(i,j,0)); +#if (defined WARPX_DIM_RZ) + if (i != 0 || rmin != 0.) { + Bx(i,j,0,0) += + dtsdz * (Ey(i,j+1,0,0) - Ey(i,j,0,0)); + } else { + Bx(i,j,0,0) = 0.; + } + for (int imode=1 ; imode < nmodes ; imode++) { + if (i == 0 && rmin == 0) { + if (imode == 1) { + // For the mode m = 1, the bulk equation diverges on axis + // (due to the 1/r terms). The following expressions regularize + // these divergences by assuming, on axis : + // Ez/r = 0/r + dEz/dr + Bx(i,j,0,2*imode-1) = Bx(i,j,0,2*imode-1) - imode*dtsdx*Ez(i+1,j,0,2*imode) & + + dtsdz*(Ey(i,j+1,0,2*imode-1) - Ey(i,j,0,2*imode-1)); + Bx(i,j,0,2*imode) = Bx(i,j,0,2*imode) + imode*dtsdx*Ez(i+1,j,0,2*imode-1) & + + dtsdz*(Ey(i,j+1,0,2*imode) - Ey(i,j,0,2*imode)); + } else { + Bx(i,j,0,2*imode-1) = 0.; + Bx(i,j,0,2*imode) = 0.; + } + } else { + // Br(i,j,m) = Br(i,j,m) + I*m*dt*Ez(i,j,m)/r + dtsdz*(Et(i,j+1,m) - Et(i,j,m)) + const Real r = rmin*dxinv + i; + Bx(i,j,0,2*imode-1) = Bx(i,j,0,2*imode-1) - imode*dtsdx*Ez(i,j,0,2*imode)/r & + + dtsdz*(Ey(i,j+1,0,2*imode-1) - Ey(i,j,0,2*imode-1)); + Bx(i,j,0,2*imode) = Bx(i,j,0,2*imode) + imode*dtsdx*Ez(i,j,0,2*imode-1)/r & + + dtsdz*(Ey(i,j+1,0,2*imode) - Ey(i,j,0,2*imode)); + } + } +#endif #endif } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void warpx_push_by_yee(int i, int j, int k, Array4<Real> const& By, Array4<Real const> const& Ex, Array4<Real const> const& Ez, - Real dtsdx, Real dtsdz) + Real dtsdx, Real dtsdz, const long nmodes) { #if defined WARPX_DIM_3D By(i,j,k) += + dtsdx * (Ez(i+1,j,k ) - Ez(i,j,k)) - dtsdz * (Ex(i ,j,k+1) - Ex(i,j,k)); #elif (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) - // Note that the 2D Cartesian and RZ versions are the same - By(i,j,0) += + dtsdx * (Ez(i+1,j ,0) - Ez(i,j,0)) - - dtsdz * (Ex(i ,j+1,0) - Ex(i,j,0)); + // Note that the 2D Cartesian and RZ mode 0 are the same + By(i,j,0,0) += + dtsdx * (Ez(i+1,j ,0,0) - Ez(i,j,0,0)) + - dtsdz * (Ex(i ,j+1,0,0) - Ex(i,j,0,0)); +#if (defined WARPX_DIM_RZ) + for (int imode=1 ; imode < nmodes ; imode++) { + // Bt(i,j,m) = Bt(i,j,m) + dtsdr*(Ez(i+1,j,m) - Ez(i,j,m)) - dtsdz*(Er(i,j+1,m) - Er(i,j,m)) + By(i,j,0,2*imode-1) += + dtsdx*(Ez(i+1,j ,0,2*imode-1) - Ez(i,j,0,2*imode-1)) + - dtsdz*(Ex(i ,j+1,0,2*imode-1) - Ex(i,j,0,2*imode-1)); + By(i,j,0,2*imode) += + dtsdx*(Ez(i+1,j ,0,2*imode) - Ez(i,j,0,2*imode)) + - dtsdz*(Ex(i ,j+1,0,2*imode) - Ex(i,j,0,2*imode)); + } +#endif #endif } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void warpx_push_bz_yee(int i, int j, int k, Array4<Real> const& Bz, Array4<Real const> const& Ex, Array4<Real const> const& Ey, - Real dtsdx, Real dtsdy, Real dxinv, Real rmin) + Real dtsdx, Real dtsdy, Real dxinv, Real rmin, const long nmodes) { #if defined WARPX_DIM_3D Bz(i,j,k) += - dtsdx * (Ey(i+1,j ,k) - Ey(i,j,k)) @@ -45,46 +84,94 @@ void warpx_push_bz_yee(int i, int j, int k, Array4<Real> const& Bz, #elif defined WARPX_DIM_XZ Bz(i,j,0) += - dtsdx * (Ey(i+1,j,0) - Ey(i,j,0)); #elif defined WARPX_DIM_RZ + const Real r = rmin*dxinv + i + 0.5; const Real ru = 1. + 0.5/(rmin*dxinv + i + 0.5); const Real rd = 1. - 0.5/(rmin*dxinv + i + 0.5); - Bz(i,j,0) += - dtsdx * (ru*Ey(i+1,j,0) - rd*Ey(i,j,0)); + Bz(i,j,0,0) += - dtsdx*(ru*Ey(i+1,j,0,0) - rd*Ey(i,j,0,0)); + for (int imode=1 ; imode < nmodes ; imode++) { + // Bz(i,j,m) = Bz(i,j,m) - dtsdr*(ru*Et(i+1,j,m) - rd*Et(i,j,m)) - I*m*dt*Er(i,j,m)/r + Bz(i,j,0,2*imode-1) += - dtsdx*(ru*Ey(i+1,j,0,2*imode-1) - rd*Ey(i,j,0,2*imode-1)) + imode*dtsdx*Ex(i,j,0,2*imode)/r; + Bz(i,j,0,2*imode) += - dtsdx*(ru*Ey(i+1,j,0,2*imode) - rd*Ey(i,j,0,2*imode)) - imode*dtsdx*Ex(i,j,0,2*imode-1)/r; + } #endif } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void warpx_push_ex_yee(int i, int j, int k, Array4<Real> const& Ex, Array4<Real const> const& By, Array4<Real const> const& Bz, Array4<Real const> const& Jx, - Real mu_c2_dt, Real dtsdy_c2, Real dtsdz_c2) + Real mu_c2_dt, Real dtsdx_c2, Real dtsdy_c2, Real dtsdz_c2, Real dxinv, Real rmin, const long nmodes) { #if defined WARPX_DIM_3D Ex(i,j,k) += + dtsdy_c2 * (Bz(i,j,k) - Bz(i,j-1,k )) - dtsdz_c2 * (By(i,j,k) - By(i,j ,k-1)) - mu_c2_dt * Jx(i,j,k); #elif (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) - // Note that the 2D Cartesian and RZ versions are the same - Ex(i,j,0) += - dtsdz_c2 * (By(i,j,0) - By(i,j-1,0)) - - mu_c2_dt * Jx(i,j,0); + // Note that the 2D Cartesian and RZ mode 0 are the same + Ex(i,j,0,0) += - dtsdz_c2 * (By(i,j,0,0) - By(i,j-1,0,0)) + - mu_c2_dt * Jx(i,j,0,0); +#if (defined WARPX_DIM_RZ) + const Real r = rmin*dxinv+ i + 0.5; + for (int imode=1 ; imode < nmodes ; imode++) { + // Er(i,j,m) = Er(i,j,m) - I*m*dt*Bz(i,j,m)/r - dtsdz*(Bt(i,j,m) - Bt(i,j-1,m)) - mudt*Jr(i,j,m) + Ex(i,j,0,2*imode-1) += - dtsdz_c2*(By(i,j,0,2*imode-1) - By(i,j-1,0,2*imode-1)) + imode*dtsdx_c2*Bz(i,j,0,2*imode)/r + - mu_c2_dt*Jx(i,j,0,2*imode-1); + Ex(i,j,0,2*imode) += - dtsdz_c2*(By(i,j,0,2*imode) - By(i,j-1,0,2*imode)) - imode*dtsdx_c2*Bz(i,j,0,2*imode-1)/r + - mu_c2_dt*Jx(i,j,0,2*imode); + } +#endif #endif } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void warpx_push_ey_yee(int i, int j, int k, Array4<Real> const& Ey, Array4<Real const> const& Bx, Array4<Real const> const& Bz, Array4<Real const> const& Jy, - Real mu_c2_dt, Real dtsdx_c2, Real dtsdz_c2, Real rmin) + Array4<Real> const& Ex, + Real mu_c2_dt, Real dtsdx_c2, Real dtsdz_c2, Real rmin, const long nmodes) { #if defined WARPX_DIM_3D Ey(i,j,k) += - dtsdx_c2 * (Bz(i,j,k) - Bz(i-1,j,k)) + dtsdz_c2 * (Bx(i,j,k) - Bx(i,j,k-1)) - mu_c2_dt * Jy(i,j,k); -#elif (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) - // 2D Cartesian and RZ are the same, except that the axis is skipped with RZ -#ifdef WARPX_DIM_RZ - if (i != 0 || rmin != 0.) -#endif - { +#elif (defined WARPX_DIM_XZ) Ey(i,j,0) += - dtsdx_c2 * (Bz(i,j,0) - Bz(i-1,j,0)) + dtsdz_c2 * (Bx(i,j,0) - Bx(i,j-1,0)) - mu_c2_dt * Jy(i,j,0); +#elif (defined WARPX_DIM_RZ) + if (i != 0 || rmin != 0.) { + Ey(i,j,0,0) += - dtsdx_c2 * (Bz(i,j,0,0) - Bz(i-1,j,0,0)) + + dtsdz_c2 * (Bx(i,j,0,0) - Bx(i,j-1,0,0)) + - mu_c2_dt * Jy(i,j,0,0); + } else { + Ey(i,j,0,0) = 0.; + } + for (int imode=1 ; imode < nmodes ; imode++) { + if (i == 0 && rmin == 0) { + if (imode == 1) { + // The bulk equation could in principle be used here since it does not diverge + // on axis. However, it typically gives poor results e.g. for the propagation + // of a laser pulse (the field is spuriously reduced on axis). For this reason + // a modified on-axis condition is used here: we use the fact that + // Etheta(r=0,m=1) should equal -iEr(r=0,m=1), for the fields Er and Et to be + // independent of theta at r=0. Now with linear interpolation: + // Er(r=0,m=1) = 0.5*[Er(r=dr/2,m=1) + Er(r=-dr/2,m=1)] + // And using the rule applying for the guards cells + // Er(r=-dr/2,m=1) = Er(r=dr/2,m=1). Thus: Et(i,j,m) = -i*Er(i,j,m) + Ey(i,j,0,2*imode-1) = Ex(i,j,0,2*imode); + Ey(i,j,0,2*imode) = -Ex(i,j,0,2*imode-1); + } else { + // Etheta should remain 0 on axis, for modes different than m=1 + Ey(i,j,0,2*imode-1) = 0.; + Ey(i,j,0,2*imode) = 0.; + } + } else { + // Et(i,j,m) = Et(i,j,m) - dtsdr*(Bz(i,j,m) - Bz(i-1,j,m)) + dtsdz*(Br(i,j,m) - Br(i,j-1,m)) - mudt*Jt(i,j,m) + Ey(i,j,0,2*imode-1) += - dtsdx_c2 * (Bz(i,j,0,2*imode-1) - Bz(i-1,j,0,2*imode-1)) + + dtsdz_c2 * (Bx(i,j,0,2*imode-1) - Bx(i,j-1,0,2*imode-1)) + - mu_c2_dt * Jy(i,j,0,2*imode-1); + Ey(i,j,0,2*imode) += - dtsdx_c2 * (Bz(i,j,0,2*imode) - Bz(i-1,j,0,2*imode)) + + dtsdz_c2 * (Bx(i,j,0,2*imode) - Bx(i,j-1,0,2*imode)) + - mu_c2_dt * Jy(i,j,0,2*imode); + } } #endif } @@ -92,7 +179,7 @@ void warpx_push_ey_yee(int i, int j, int k, Array4<Real> const& Ey, AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void warpx_push_ez_yee(int i, int j, int k, Array4<Real> const& Ez, Array4<Real const> const& Bx, Array4<Real const> const& By, Array4<Real const> const& Jz, - Real mu_c2_dt, Real dtsdx_c2, Real dtsdy_c2, Real dxinv, Real rmin) + Real mu_c2_dt, Real dtsdx_c2, Real dtsdy_c2, Real dxinv, Real rmin, const long nmodes) { #if defined WARPX_DIM_3D Ez(i,j,k) += + dtsdx_c2 * (By(i,j,k) - By(i-1,j ,k)) @@ -105,11 +192,28 @@ void warpx_push_ez_yee(int i, int j, int k, Array4<Real> const& Ez, if (i != 0 || rmin != 0.) { const Real ru = 1. + 0.5/(rmin*dxinv + i); const Real rd = 1. - 0.5/(rmin*dxinv + i); - Ez(i,j,0) += + dtsdx_c2 * (ru*By(i,j,0) - rd*By(i-1,j,0)) - - mu_c2_dt * Jz(i,j,0); + Ez(i,j,0,0) += + dtsdx_c2 * (ru*By(i,j,0,0) - rd*By(i-1,j,0,0)) + - mu_c2_dt * Jz(i,j,0,0); } else { - Ez(i,j,0) += + 4.*dtsdx_c2 * By(i,j,0) - - mu_c2_dt * Jz(i,j,0); + Ez(i,j,0,0) += + 4.*dtsdx_c2 * By(i,j,0,0) + - mu_c2_dt * Jz(i,j,0,0); + } + for (int imode=1 ; imode < nmodes ; imode++) { + if (i == 0 && rmin == 0) { + Ez(i,j,0,2*imode-1) = 0.; + Ez(i,j,0,2*imode) = 0.; + } else { + const Real r = rmin*dxinv + i + 0.5; + const Real ru = 1. + 0.5/(rmin*dxinv + i); + const Real rd = 1. - 0.5/(rmin*dxinv + i); + // Ez(i,j,m) = Ez(i,j,m) + dtsdr*(ru*Bt(i,j,m) - rd*Bt(i-1,j,m)) + I*m*dt*Br(i,j,m)/r - mudt*Jz(i,j,m) + Ez(i,j,0,2*imode-1) += + dtsdx_c2 * (ru*By(i,j,0,2*imode-1) - rd*By(i-1,j,0,2*imode-1)) + - imode*dtsdx_c2*Bx(i,j,0,2*imode)/r + - mu_c2_dt * Jz(i,j,0,2*imode-1); + Ez(i,j,0,2*imode) += + dtsdx_c2 * (ru*By(i,j,0,2*imode) - rd*By(i-1,j,0,2*imode)) + + imode*dtsdx_c2*Bx(i,j,0,2*imode-1)/r + - mu_c2_dt * Jz(i,j,0,2*imode); + } } #endif } diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H index 60a9817a2..26da42606 100644 --- a/Source/FortranInterface/WarpX_f.H +++ b/Source/FortranInterface/WarpX_f.H @@ -63,6 +63,7 @@ extern "C" { #endif + // Laser pusher void warpx_gaussian_laser( const long* np, amrex::Real* Xp, amrex::Real* Yp, amrex::Real* t, amrex::Real* wavelength, amrex::Real* e_max, amrex::Real* waist, diff --git a/Source/Initialization/InjectorPosition.H b/Source/Initialization/InjectorPosition.H index 19bb092dd..0e43f889e 100644 --- a/Source/Initialization/InjectorPosition.H +++ b/Source/Initialization/InjectorPosition.H @@ -33,7 +33,7 @@ struct InjectorPositionRegular { int nx = ref_fac*ppc.x; int ny = ref_fac*ppc.y; -#if (AMREX_SPACEDIM == 3) +#if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ) int nz = ref_fac*ppc.z; #else int nz = 1; diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp index 541999789..a65899902 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 WARPX_DIM_XZ 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 feb9d3564..88adbc147 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -62,24 +62,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]; @@ -137,24 +137,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 @@ -202,8 +202,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()); } } } @@ -418,7 +418,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), @@ -514,11 +514,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()); } } } @@ -548,7 +548,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]) { @@ -556,18 +556,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 { @@ -575,29 +575,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/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index b879970ef..c7dfde75a 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -2,6 +2,7 @@ #define CURRENTDEPOSITION_H_ #include "ShapeFactors.H" +#include <WarpX_Complex.H> /* \brief Current Deposition for thread thread_num * /param xp, yp, zp : Pointer to arrays of particle positions. @@ -169,7 +170,7 @@ void doDepositionShapeN(const amrex::Real * const xp, } /* \brief Esirkepov Current Deposition for thread thread_num - * /param xp, yp, zp : Pointer to arrays of particle positions. + * \param xp, yp, zp : Pointer to arrays of particle positions. * \param wp : Pointer to array of particle weights. * \param uxp uyp uzp : Pointer to arrays of particle momentum. * \param ion_lev : Pointer to array of particle ionization level. This is @@ -184,7 +185,8 @@ void doDepositionShapeN(const amrex::Real * const xp, * \param dx : 3D cell size * \param xyzmin : Physical lower bounds of domain. * \param lo : Index lower bounds of domain. - * /param q : species charge. + * \param q : species charge. + * \param n_rz_azimuthal_modes: Number of azimuthal modes when using RZ geometry */ template <int depos_order> void doEsirkepovDepositionShapeN (const amrex::Real * const xp, @@ -203,7 +205,8 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, const std::array<amrex::Real,3>& dx, const std::array<amrex::Real, 3> xyzmin, const amrex::Dim3 lo, - const amrex::Real q) + const amrex::Real q, + const long n_rz_azimuthal_modes) { // Whether ion_lev is a null pointer (do_ionization=0) or a real pointer // (do_ionization=1) @@ -230,6 +233,10 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, const amrex::Real invvol = 1.0/(dx[0]*dx[2]); #endif +#if (defined WARPX_DIM_RZ) + const Complex I = Complex{0., 1.}; +#endif + const amrex::Real clightsq = 1.0/PhysConst::c/PhysConst::c; // Loop over particles and deposit into Jx_arr, Jy_arr and Jz_arr @@ -255,11 +262,42 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, // computes current and old position in grid units #if (defined WARPX_DIM_RZ) - const amrex::Real r_new = std::sqrt(xp[ip]*xp[ip] + yp[ip]*yp[ip]); - const amrex::Real r_old = std::sqrt((xp[ip] - dt*uxp[ip]*gaminv)*(xp[ip] - dt*uxp[ip]*gaminv) + - (yp[ip] - dt*uyp[ip]*gaminv)*(yp[ip] - dt*uyp[ip]*gaminv)); - const amrex::Real x_new = (r_new - xmin)*dxi; - const amrex::Real x_old = (r_old - xmin)*dxi; + const amrex::Real xp_mid = xp[ip] - 0.5*dt*uxp[ip]*gaminv; + const amrex::Real yp_mid = yp[ip] - 0.5*dt*uyp[ip]*gaminv; + const amrex::Real xp_old = xp[ip] - dt*uxp[ip]*gaminv; + const amrex::Real yp_old = yp[ip] - dt*uyp[ip]*gaminv; + const amrex::Real rp_new = std::sqrt(xp[ip]*xp[ip] + yp[ip]*yp[ip]); + const amrex::Real rp_mid = std::sqrt(xp_mid*xp_mid + yp_mid*yp_mid); + const amrex::Real rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old); + amrex::Real costheta_new, sintheta_new; + if (rp_new > 0.) { + costheta_new = xp[ip]/rp_new; + sintheta_new = yp[ip]/rp_new; + } else { + costheta_new = 1.; + sintheta_new = 0.; + } + amrex::Real costheta_mid, sintheta_mid; + if (rp_mid > 0.) { + costheta_mid = xp_mid/rp_mid; + sintheta_mid = yp_mid/rp_mid; + } else { + costheta_mid = 1.; + sintheta_mid = 0.; + } + amrex::Real costheta_old, sintheta_old; + if (rp_old > 0.) { + costheta_old = xp_old/rp_old; + sintheta_old = yp_old/rp_old; + } else { + costheta_old = 1.; + sintheta_old = 0.; + } + const Complex xy_new0 = Complex{costheta_new, sintheta_new}; + const Complex xy_mid0 = Complex{costheta_mid, sintheta_mid}; + const Complex xy_old0 = Complex{costheta_old, sintheta_old}; + const amrex::Real x_new = (rp_new - xmin)*dxi; + const amrex::Real x_old = (rp_old - xmin)*dxi; #else const amrex::Real x_new = (xp[ip] - xmin)*dxi; const amrex::Real x_old = x_new - dtsdx0*uxp[ip]*gaminv; @@ -272,16 +310,7 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, const amrex::Real z_old = z_new - dtsdz0*uzp[ip]*gaminv; #if (defined WARPX_DIM_RZ) - amrex::Real costheta; - amrex::Real sintheta; - if (r_new > 0.) { - costheta = xp[ip]/r_new; - sintheta = yp[ip]/r_new; - } else { - costheta = 1.; - sintheta = 0.; - } - const amrex::Real vy = (-uxp[ip]*sintheta + uyp[ip]*costheta)*gaminv; + const amrex::Real vy = (-uxp[ip]*sintheta_mid + uyp[ip]*costheta_mid)*gaminv; #elif (defined WARPX_DIM_XZ) const amrex::Real vy = uyp[ip]*gaminv; #endif @@ -363,25 +392,61 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, amrex::Real sdxi = 0.; for (int i=dil; i<=depos_order+1-diu; i++) { sdxi += wqx*(sx_old[i] - sx_new[i])*(sz_new[k] + 0.5*(sz_old[k] - sz_new[k])); - amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0), sdxi); + amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdxi); +#if (defined WARPX_DIM_RZ) + Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta} + for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) { + // The factor 2 comes from the normalization of the modes + const Complex djr_cmplx = 2.*sdxi*xy_mid; + amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djr_cmplx.real()); + amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djr_cmplx.imag()); + xy_mid = xy_mid*xy_mid0; + } +#endif } } for (int k=dkl; k<=depos_order+2-dku; k++) { for (int i=dil; i<=depos_order+2-diu; i++) { const amrex::Real sdyj = wq*vy*invvol*((sz_new[k] + 0.5*(sz_old[k] - sz_new[k]))*sx_new[i] + (0.5*sz_new[k] + 1./3.*(sz_old[k] - sz_new[k]))*(sx_old[i] - sx_new[i])); - amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0), sdyj); + amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdyj); +#if (defined WARPX_DIM_RZ) + Complex xy_new = xy_new0; + Complex xy_mid = xy_mid0; + Complex xy_old = xy_old0; + // Throughout the following loop, xy_ takes the value e^{i m theta_} + for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) { + // The factor 2 comes from the normalization of the modes + // The minus sign comes from the different convention with respect to Davidson et al. + const Complex djt_cmplx = -2.*I*(i_new-1 + i + xmin*dxi)*wq*invdtdx/(double)imode* + (sx_new[i]*sz_new[k]*(xy_new - xy_mid) + sx_old[i]*sz_old[k]*(xy_mid - xy_old)); + amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djt_cmplx.real()); + amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djt_cmplx.imag()); + xy_new = xy_new*xy_new0; + xy_mid = xy_mid*xy_mid0; + xy_old = xy_old*xy_old0; + } +#endif } } for (int i=dil; i<=depos_order+2-diu; i++) { amrex::Real sdzk = 0.; for (int k=dkl; k<=depos_order+1-dku; k++) { sdzk += wqz*(sz_old[k] - sz_new[k])*(sx_new[i] + 0.5*(sx_old[i] - sx_new[i])); - amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0), sdzk); + amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdzk); +#if (defined WARPX_DIM_RZ) + Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta} + for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) { + // The factor 2 comes from the normalization of the modes + const Complex djz_cmplx = 2.*sdzk*xy_mid; + amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djz_cmplx.real()); + amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djz_cmplx.imag()); + xy_mid = xy_mid*xy_mid0; + } +#endif } } - #endif } ); diff --git a/Source/Particles/Gather/FieldGather.H b/Source/Particles/Gather/FieldGather.H index 1978e8141..d8d1d78ef 100644 --- a/Source/Particles/Gather/FieldGather.H +++ b/Source/Particles/Gather/FieldGather.H @@ -2,9 +2,10 @@ #define FIELDGATHER_H_ #include "ShapeFactors.H" +#include <WarpX_Complex.H> /* \brief Field gather for particles handled by thread thread_num - * /param xp, yp, zp : Pointer to arrays of particle positions. + * \param xp, yp, zp : Pointer to arrays of particle positions. * \param Exp, Eyp, Ezp: Pointer to array of electric field on particles. * \param Bxp, Byp, Bzp: Pointer to array of magnetic field on particles. * \param ex_arr ey_arr: Array4 of current density, either full array or tile. @@ -15,6 +16,7 @@ * \param xyzmin : Physical lower bounds of domain. * \param lo : Index lower bounds of domain. * \param stagger_shift: 0 if nodal, 0.5 if staggered. + * \param n_rz_azimuthal_modes: Number of azimuthal modes when using RZ geometry */ template <int depos_order, int lower_in_v> void doGatherShapeN(const amrex::Real * const xp, @@ -33,7 +35,8 @@ void doGatherShapeN(const amrex::Real * const xp, const std::array<amrex::Real, 3>& dx, const std::array<amrex::Real, 3> xyzmin, const amrex::Dim3 lo, - const amrex::Real stagger_shift) + const amrex::Real stagger_shift, + const long n_rz_azimuthal_modes) { const amrex::Real dxi = 1.0/dx[0]; const amrex::Real dzi = 1.0/dx[2]; @@ -56,8 +59,8 @@ void doGatherShapeN(const amrex::Real * const xp, // x direction // Get particle position #ifdef WARPX_DIM_RZ - const amrex::Real r = std::sqrt(xp[ip]*xp[ip] + yp[ip]*yp[ip]); - const amrex::Real x = (r - xmin)*dxi; + const amrex::Real rp = std::sqrt(xp[ip]*xp[ip] + yp[ip]*yp[ip]); + const amrex::Real x = (rp - xmin)*dxi; #else const amrex::Real x = (xp[ip]-xmin)*dxi; #endif @@ -103,7 +106,7 @@ void doGatherShapeN(const amrex::Real * const xp, for (int iz=0; iz<=depos_order; iz++){ for (int ix=0; ix<=depos_order; ix++){ Eyp[ip] += sx[ix]*sz[iz]* - ey_arr(lo.x+j+ix, lo.y+l+iz, 0); + ey_arr(lo.x+j+ix, lo.y+l+iz, 0, 0); } } // Gather field on particle Exp[i] from field on grid ex_arr @@ -111,9 +114,9 @@ void doGatherShapeN(const amrex::Real * const xp, for (int iz=0; iz<=depos_order; iz++){ for (int ix=0; ix<=depos_order-lower_in_v; ix++){ Exp[ip] += sx0[ix]*sz[iz]* - ex_arr(lo.x+j0+ix, lo.y+l +iz, 0); + ex_arr(lo.x+j0+ix, lo.y+l +iz, 0, 0); Bzp[ip] += sx0[ix]*sz[iz]* - bz_arr(lo.x+j0+ix, lo.y+l +iz, 0); + bz_arr(lo.x+j0+ix, lo.y+l +iz, 0, 0); } } // Gather field on particle Ezp[i] from field on grid ez_arr @@ -121,30 +124,79 @@ void doGatherShapeN(const amrex::Real * const xp, for (int iz=0; iz<=depos_order-lower_in_v; iz++){ for (int ix=0; ix<=depos_order; ix++){ Ezp[ip] += sx[ix]*sz0[iz]* - ez_arr(lo.x+j+ix, lo.y+l0 +iz, 0); + ez_arr(lo.x+j+ix, lo.y+l0 +iz, 0, 0); Bxp[ip] += sx[ix]*sz0[iz]* - bx_arr(lo.x+j+ix, lo.y+l0 +iz, 0); + bx_arr(lo.x+j+ix, lo.y+l0 +iz, 0, 0); } } // Gather field on particle Byp[i] from field on grid by_arr for (int iz=0; iz<=depos_order-lower_in_v; iz++){ for (int ix=0; ix<=depos_order-lower_in_v; ix++){ Byp[ip] += sx0[ix]*sz0[iz]* - by_arr(lo.x+j0+ix, lo.y+l0+iz, 0); + by_arr(lo.x+j0+ix, lo.y+l0+iz, 0, 0); } } #ifdef WARPX_DIM_RZ - // Convert Exp and Eyp (which are actually Er and Etheta) to Ex and Ey + amrex::Real costheta; amrex::Real sintheta; - if (r > 0.) { - costheta = xp[ip]/r; - sintheta = yp[ip]/r; + if (rp > 0.) { + costheta = xp[ip]/rp; + sintheta = yp[ip]/rp; } else { costheta = 1.; sintheta = 0.; } + const Complex xy0 = Complex{costheta, -sintheta}; + Complex xy = xy0; + + for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) { + + // Gather field on particle Eyp[i] from field on grid ey_arr + for (int iz=0; iz<=depos_order; iz++){ + for (int ix=0; ix<=depos_order; ix++){ + const amrex::Real dEy = (+ ey_arr(lo.x+j+ix, lo.y+l+iz, 0, 2*imode-1)*xy.real() + - ey_arr(lo.x+j+ix, lo.y+l+iz, 0, 2*imode)*xy.imag()); + Eyp[ip] += sx[ix]*sz[iz]*dEy; + } + } + // Gather field on particle Exp[i] from field on grid ex_arr + // Gather field on particle Bzp[i] from field on grid bz_arr + for (int iz=0; iz<=depos_order; iz++){ + for (int ix=0; ix<=depos_order-lower_in_v; ix++){ + const amrex::Real dEx = (+ ex_arr(lo.x+j0+ix, lo.y+l +iz, 0, 2*imode-1)*xy.real() + - ex_arr(lo.x+j0+ix, lo.y+l +iz, 0, 2*imode)*xy.imag()); + Exp[ip] += sx0[ix]*sz[iz]*dEx; + const amrex::Real dBz = (+ bz_arr(lo.x+j0+ix, lo.y+l +iz, 0, 2*imode-1)*xy.real() + - bz_arr(lo.x+j0+ix, lo.y+l +iz, 0, 2*imode)*xy.imag()); + Bzp[ip] += sx0[ix]*sz[iz]*dBz; + } + } + // Gather field on particle Ezp[i] from field on grid ez_arr + // Gather field on particle Bxp[i] from field on grid bx_arr + for (int iz=0; iz<=depos_order-lower_in_v; iz++){ + for (int ix=0; ix<=depos_order; ix++){ + const amrex::Real dEz = (+ ez_arr(lo.x+j+ix, lo.y+l0 +iz, 0, 2*imode-1)*xy.real() + - ez_arr(lo.x+j+ix, lo.y+l0 +iz, 0, 2*imode)*xy.imag()); + Ezp[ip] += sx[ix]*sz0[iz]*dEz; + const amrex::Real dBx = (+ bx_arr(lo.x+j+ix, lo.y+l0 +iz, 0, 2*imode-1)*xy.real() + - bx_arr(lo.x+j+ix, lo.y+l0 +iz, 0, 2*imode)*xy.imag()); + Bxp[ip] += sx[ix]*sz0[iz]*dBx; + } + } + // Gather field on particle Byp[i] from field on grid by_arr + for (int iz=0; iz<=depos_order-lower_in_v; iz++){ + for (int ix=0; ix<=depos_order-lower_in_v; ix++){ + const amrex::Real dBy = (+ by_arr(lo.x+j0+ix, lo.y+l0+iz, 0, 2*imode-1)*xy.real() + - by_arr(lo.x+j0+ix, lo.y+l0+iz, 0, 2*imode)*xy.imag()); + Byp[ip] += sx0[ix]*sz0[iz]*dBy; + } + } + xy = xy*xy0; + } + + // Convert Exp and Eyp (which are actually Er and Etheta) to Ex and Ey const amrex::Real Exp_save = Exp[ip]; Exp[ip] = costheta*Exp[ip] - sintheta*Eyp[ip]; Eyp[ip] = costheta*Eyp[ip] + sintheta*Exp_save; diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 143f4b7f9..7803bdae1 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -257,7 +257,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.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 015482e3f..2dc25e6fa 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -160,12 +160,12 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, npart /= 4; } for (long i = 0; i < npart; ++i) { -#if ( AMREX_SPACEDIM == 3 | WARPX_DIM_RZ) +#if (defined WARPX_DIM_3D) || (WARPX_DIM_RZ) Real weight = q_tot/npart/charge; Real x = distx(mt); Real y = disty(mt); Real z = distz(mt); -#elif ( AMREX_SPACEDIM == 2 ) +#elif (defined WARPX_DIM_XZ) Real weight = q_tot/npart/charge/y_rms; Real x = distx(mt); Real y = 0.; @@ -332,6 +332,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) Real density_max = plasma_injector->density_max; #ifdef WARPX_DIM_RZ + const long nmodes = WarpX::n_rz_azimuthal_modes; bool radially_weighted = plasma_injector->radially_weighted; #endif @@ -489,7 +490,12 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) #else Real x = overlap_corner[0] + (iv[0]+r.x)*dx[0]; Real y = 0.0; +#if defined WARPX_DIM_XZ Real z = overlap_corner[1] + (iv[1]+r.y)*dx[1]; +#elif defined WARPX_DIM_RZ + // Note that for RZ, r.y will be theta + Real z = overlap_corner[1] + (iv[1]+r.z)*dx[1]; +#endif #endif #if (AMREX_SPACEDIM == 3) @@ -510,9 +516,16 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) Real yb = y; #ifdef WARPX_DIM_RZ - // Replace the x and y, choosing the angle randomly. + // Replace the x and y, setting an angle theta. // These x and y are used to get the momentum and density - Real theta = 2.*MathConst::pi*amrex::Random(); + Real theta; + if (nmodes == 1) { + // With only 1 mode, the angle doesn't matter so + // choose it randomly. + theta = 2.*MathConst::pi*amrex::Random(); + } else { + theta = 2.*MathConst::pi*r.y; + } x = xb*std::cos(theta); y = xb*std::sin(theta); #endif @@ -903,7 +916,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, + 0, np, thread_num, lev, lev); if (cost) { const Box& tbx = pti.tilebox(); @@ -1201,7 +1215,8 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_START(blp_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, + 0, np_gather, thread_num, lev, lev); if (np_gather < np) { @@ -1635,7 +1650,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, + 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 @@ -1928,7 +1944,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, np_to_gather, dx, - xyzmin, lo, stagger_shift); + xyzmin, lo, stagger_shift, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 2){ doGatherShapeN<2,1>(xp, yp, zp, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, @@ -1936,7 +1952,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, np_to_gather, dx, - xyzmin, lo, stagger_shift); + xyzmin, lo, stagger_shift, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 3){ doGatherShapeN<3,1>(xp, yp, zp, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, @@ -1944,7 +1960,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, np_to_gather, dx, - xyzmin, lo, stagger_shift); + xyzmin, lo, stagger_shift, WarpX::n_rz_azimuthal_modes); } } else { if (WarpX::nox == 1){ @@ -1954,7 +1970,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, np_to_gather, dx, - xyzmin, lo, stagger_shift); + xyzmin, lo, stagger_shift, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 2){ doGatherShapeN<2,0>(xp, yp, zp, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, @@ -1962,7 +1978,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, np_to_gather, dx, - xyzmin, lo, stagger_shift); + xyzmin, lo, stagger_shift, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 3){ doGatherShapeN<3,0>(xp, yp, zp, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, @@ -1970,7 +1986,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, np_to_gather, dx, - xyzmin, lo, stagger_shift); + xyzmin, lo, stagger_shift, WarpX::n_rz_azimuthal_modes); } } } diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index 0c94d1e33..4893b3294 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -409,7 +409,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, + 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 7393f7301..c39eec9dc 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -51,6 +51,9 @@ namespace ParticleStringNames {"Bx", PIdx::Bx }, {"By", PIdx::By }, {"Bz", PIdx::Bz } +#ifdef WARPX_DIM_RZ + ,{"theta", PIdx::theta} +#endif }; } diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 80d1caa0f..4e80374d8 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -336,9 +336,9 @@ WarpXParticleContainer::DepositCurrent(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()); // local_jx[thread_num] is set to zero local_jx[thread_num].setVal(0.0); @@ -373,17 +373,20 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, doEsirkepovDepositionShapeN<1>( xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q); + jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q, + WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 2){ doEsirkepovDepositionShapeN<2>( xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q); + jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q, + WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 3){ doEsirkepovDepositionShapeN<3>( xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q); + jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q, + WarpX::n_rz_azimuthal_modes); } } else { if (WarpX::nox == 1){ @@ -412,9 +415,9 @@ WarpXParticleContainer::DepositCurrent(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, jx->nComp()); + (*jy)[pti].atomicAdd(local_jy[thread_num], tby, tby, 0, 0, jy->nComp()); + (*jz)[pti].atomicAdd(local_jz[thread_num], tbz, tbz, 0, 0, jz->nComp()); BL_PROFILE_VAR_STOP(blp_accumulate); #endif } @@ -471,15 +474,17 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, tilebox.grow(ngRho); + const int nc = (rho->nComp() == 1 ? 1 : rho->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); @@ -515,7 +520,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 @@ -569,7 +574,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 @@ -598,7 +603,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/Python/WarpXWrappers.cpp b/Source/Python/WarpXWrappers.cpp index 3ed4830f5..a60efd498 100644 --- a/Source/Python/WarpXWrappers.cpp +++ b/Source/Python/WarpXWrappers.cpp @@ -29,7 +29,7 @@ namespace for (int j = 0; j < AMREX_SPACEDIM; ++j) { (*shapes)[shapesize*i+j] = mf[mfi].box().length(j); } - if (mf.nComp() > 1) (*shapes)[shapesize*i+2] = mf.nComp(); + if (mf.nComp() > 1) (*shapes)[shapesize*i+AMREX_SPACEDIM] = mf.nComp(); } return data; } diff --git a/Source/Utils/Make.package b/Source/Utils/Make.package index 389b3670a..8c071a08b 100644 --- a/Source/Utils/Make.package +++ b/Source/Utils/Make.package @@ -7,6 +7,7 @@ CEXE_headers += WarpXUtil.H CEXE_headers += WarpXAlgorithmSelection.H CEXE_sources += WarpXAlgorithmSelection.cpp CEXE_headers += NCIGodfreyTables.H +CEXE_headers += WarpX_Complex.H CEXE_headers += IonizationEnergiesTable.H INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Utils diff --git a/Source/FieldSolver/SpectralSolver/WarpX_Complex.H b/Source/Utils/WarpX_Complex.H index c898c5baa..f04e134c7 100644 --- a/Source/FieldSolver/SpectralSolver/WarpX_Complex.H +++ b/Source/Utils/WarpX_Complex.H @@ -7,20 +7,15 @@ #ifdef AMREX_USE_GPU #include <thrust/complex.h> -#include <cufft.h> using Complex = thrust::complex<amrex::Real>; -static_assert( sizeof(Complex) == sizeof(cuDoubleComplex), - "The complex types in WarpX and cuFFT do not match."); #else #include <complex> -#include <fftw3.h> using Complex = std::complex<amrex::Real>; -static_assert( sizeof(Complex) == sizeof(fftw_complex), - "The complex types in WarpX and FFTW do not match."); #endif + static_assert(sizeof(Complex) == sizeof(amrex::Real[2]), "Unexpected complex type."); diff --git a/Source/WarpX.H b/Source/WarpX.H index 36d2a8f35..9236b8e20 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -94,6 +94,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 517fb2332..95826c075 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -44,6 +44,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; @@ -517,6 +520,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_DIM_RZ, otherwise defaults to 1. + pp.query("n_rz_azimuthal_modes", n_rz_azimuthal_modes); + } { @@ -782,20 +789,30 @@ 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_DIM_RZ + // With RZ multimode, there is a real and imaginary component + // for each mode, except mode 0 which is purely real + // Component 0 is mode 0. + // Odd components are the real parts. + // Even components are the imaginary parts. + ncomps = n_rz_azimuthal_modes*2 - 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)); @@ -804,25 +821,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){ @@ -848,19 +865,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)); } // @@ -872,19 +889,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)); @@ -892,17 +909,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)); } if (fft_hybrid_mpi_decomposition == false){ @@ -933,28 +950,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 (rho_cp[lev]) { - 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. } |