diff options
Diffstat (limited to 'Source')
26 files changed, 1583 insertions, 336 deletions
diff --git a/Source/Diagnostics/FieldIO.cpp b/Source/Diagnostics/FieldIO.cpp index e3d44d1fc..b1181f22f 100644 --- a/Source/Diagnostics/FieldIO.cpp +++ b/Source/Diagnostics/FieldIO.cpp @@ -677,7 +677,7 @@ getInterpolatedScalar( if ( F_fp.is_nodal() ){ IntVect refinement_vector{AMREX_D_DECL(r_ratio, r_ratio, r_ratio)}; node_bilinear_interp.interp(cfab, 0, ffab, 0, 1, - finebx, refinement_vector, {}, {}, {}, 0, 0); + finebx, refinement_vector, {}, {}, {}, 0, 0, RunOn::Cpu); } else { amrex::Abort("Unknown field staggering."); } diff --git a/Source/Diagnostics/Make.package b/Source/Diagnostics/Make.package index 3b3b3d202..dfd947d53 100644 --- a/Source/Diagnostics/Make.package +++ b/Source/Diagnostics/Make.package @@ -5,6 +5,8 @@ CEXE_sources += FieldIO.cpp CEXE_headers += FieldIO.H CEXE_headers += BoostedFrameDiagnostic.H CEXE_headers += ElectrostaticIO.cpp +CEXE_headers += SliceDiagnostic.H +CEXE_sources += SliceDiagnostic.cpp INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Diagnostics VPATH_LOCATIONS += $(WARPX_HOME)/Source/Diagnostics diff --git a/Source/Diagnostics/SliceDiagnostic.H b/Source/Diagnostics/SliceDiagnostic.H new file mode 100644 index 000000000..31eea83be --- /dev/null +++ b/Source/Diagnostics/SliceDiagnostic.H @@ -0,0 +1,42 @@ +#ifndef WARPX_SliceDiagnostic_H_ +#define WARPX_SliceDiagnostic_H_ + +#include <AMReX_VisMF.H> +#include <AMReX_PlotFileUtil.H> +#include <AMReX_ParallelDescriptor.H> +#include <AMReX_Geometry.H> + +#include <WarpX.H> +#include <AMReX_FArrayBox.H> +#include <AMReX_IArrayBox.H> +#include <AMReX_Vector.H> +#include <AMReX_BLassert.H> +#include <AMReX_MultiFabUtil.H> +#include <AMReX_MultiFabUtil_C.H> + + + +std::unique_ptr<amrex::MultiFab> CreateSlice( const amrex::MultiFab& mf, + const amrex::Vector<amrex::Geometry> &dom_geom, + amrex::RealBox &slice_realbox, + amrex::IntVect &slice_cr_ratio ); + +void CheckSliceInput( const amrex::RealBox real_box, + amrex::RealBox &slice_cc_nd_box, amrex::RealBox &slice_realbox, + amrex::IntVect &slice_cr_ratio, amrex::Vector<amrex::Geometry> dom_geom, + amrex::IntVect const SliceType, amrex::IntVect &slice_lo, + amrex::IntVect &slice_hi, amrex::IntVect &interp_lo); + +void InterpolateSliceValues( amrex::MultiFab& smf, + amrex::IntVect interp_lo, amrex::RealBox slice_realbox, + amrex::Vector<amrex::Geometry> geom, int ncomp, int nghost, + amrex::IntVect slice_lo, amrex::IntVect slice_hi, + amrex::IntVect SliceType, const amrex::RealBox real_box); + +void InterpolateLo( const amrex::Box& bx, amrex::FArrayBox &fabox, + amrex::IntVect slice_lo, amrex::Vector<amrex::Geometry> geom, + int idir, amrex::IntVect IndType, amrex::RealBox slice_realbox, + int srccomp, int ncomp, int nghost, const amrex::RealBox real_box); + +#endif + diff --git a/Source/Diagnostics/SliceDiagnostic.cpp b/Source/Diagnostics/SliceDiagnostic.cpp new file mode 100644 index 000000000..994f990c6 --- /dev/null +++ b/Source/Diagnostics/SliceDiagnostic.cpp @@ -0,0 +1,475 @@ +#include "SliceDiagnostic.H" +#include <AMReX_MultiFabUtil.H> +#include <AMReX_PlotFileUtil.H> +#include <AMReX_FillPatchUtil_F.H> + +#include <WarpX.H> + +using namespace amrex; + + +/* \brief + * The functions creates the slice for diagnostics based on the user-input. + * The slice can be 1D, 2D, or 3D and it inherts the index type of the underlying data. + * The implementation assumes that the slice is aligned with the coordinate axes. + * The input parameters are modified if the user-input does not comply with requirements of coarsenability or if the slice extent is not contained within the simulation domain. + * First a slice multifab (smf) with cell size equal to that of the simulation grid is created such that it extends from slice.dim_lo to slice.dim_hi and shares the same index space as the source multifab (mf) + * The values are copied from src mf to dst smf using amrex::ParallelCopy + * If interpolation is required, then on the smf, using data points stored in the ghost cells, the data in interpolated. + * If coarsening is required, then a coarse slice multifab is generated (cs_mf) and the + * values of the refined slice (smf) is averaged down to obtain the coarse slice. + * \param mf is the source multifab containing the field data + * \param dom_geom is the geometry of the domain and used in the function to obtain the + * CellSize of the underlying grid. + * \param slice_realbox defines the extent of the slice + * \param slice_cr_ratio provides the coarsening ratio for diagnostics + */ + +std::unique_ptr<MultiFab> +CreateSlice( const MultiFab& mf, const Vector<Geometry> &dom_geom, + RealBox &slice_realbox, IntVect &slice_cr_ratio ) +{ + std::unique_ptr<MultiFab> smf; + std::unique_ptr<MultiFab> cs_mf; + + Vector<int> slice_ncells(AMREX_SPACEDIM); + int nghost = 1; + int nlevels = dom_geom.size(); + int ncomp = (mf).nComp(); + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( nlevels==1, + "Slice diagnostics does not work with mesh refinement yet (TO DO)."); + + const auto conversionType = (mf).ixType(); + IntVect SliceType(AMREX_D_DECL(0,0,0)); + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim ) + { + SliceType[idim] = conversionType.nodeCentered(idim); + } + + const RealBox& real_box = dom_geom[0].ProbDomain(); + RealBox slice_cc_nd_box; + int slice_grid_size = 32; + + bool interpolate = false; + bool coarsen = false; + + // same index space as domain // + IntVect slice_lo(AMREX_D_DECL(0,0,0)); + IntVect slice_hi(AMREX_D_DECL(1,1,1)); + IntVect interp_lo(AMREX_D_DECL(0,0,0)); + + CheckSliceInput(real_box, slice_cc_nd_box, slice_realbox, slice_cr_ratio, + dom_geom, SliceType, slice_lo, + slice_hi, interp_lo); + // Determine if interpolation is required and number of cells in slice // + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + + // Flag for interpolation if required // + if ( interp_lo[idim] == 1) { + interpolate = 1; + } + + // For the case when a dimension is reduced // + if ( ( slice_hi[idim] - slice_lo[idim]) == 1) { + slice_ncells[idim] = 1; + } + else { + slice_ncells[idim] = ( slice_hi[idim] - slice_lo[idim] + 1 ) + / slice_cr_ratio[idim]; + + int refined_ncells = slice_hi[idim] - slice_lo[idim] + 1 ; + if ( slice_cr_ratio[idim] > 1) { + coarsen = true; + + // modify slice_grid_size if >= refines_cells // + if ( slice_grid_size >= refined_ncells ) { + slice_grid_size = refined_ncells - 1; + } + + } + } + } + + // Slice generation with index type inheritance // + Box slice(slice_lo, slice_hi); + + Vector<BoxArray> sba(1); + sba[0].define(slice); + sba[0].maxSize(slice_grid_size); + + // Distribution mapping for slice can be different from that of domain // + Vector<DistributionMapping> sdmap(1); + sdmap[0] = DistributionMapping{sba[0]}; + + smf.reset(new MultiFab(amrex::convert(sba[0],SliceType), sdmap[0], + ncomp, nghost)); + + // Copy data from domain to slice that has same cell size as that of // + // the domain mf. src and dst have the same number of ghost cells // + smf->ParallelCopy(mf, 0, 0, ncomp,nghost,nghost); + + // inteprolate if required on refined slice // + if (interpolate == 1 ) { + InterpolateSliceValues( *smf, interp_lo, slice_cc_nd_box, dom_geom, + ncomp, nghost, slice_lo, slice_hi, SliceType, real_box); + } + + + if (coarsen == false) { + return smf; + } + else if ( coarsen == true ) { + Vector<BoxArray> crse_ba(1); + crse_ba[0] = sba[0]; + crse_ba[0].coarsen(slice_cr_ratio); + + AMREX_ALWAYS_ASSERT(crse_ba[0].size() == sba[0].size()); + + cs_mf.reset( new MultiFab(amrex::convert(crse_ba[0],SliceType), + sdmap[0], ncomp,nghost)); + + MultiFab& mfSrc = *smf; + MultiFab& mfDst = *cs_mf; + + MFIter mfi_dst(mfDst); + for (MFIter mfi(mfSrc); mfi.isValid(); ++mfi) { + + FArrayBox& Src_fabox = mfSrc[mfi]; + + const Box& Dst_bx = mfi_dst.validbox(); + FArrayBox& Dst_fabox = mfDst[mfi_dst]; + + int scomp = 0; + int dcomp = 0; + + IntVect cctype(AMREX_D_DECL(0,0,0)); + if( SliceType==cctype ) { + amrex::amrex_avgdown(Dst_bx, Dst_fabox, Src_fabox, dcomp, scomp, + ncomp, slice_cr_ratio); + } + IntVect ndtype(AMREX_D_DECL(1,1,1)); + if( SliceType == ndtype ) { + amrex::amrex_avgdown_nodes(Dst_bx, Dst_fabox, Src_fabox, dcomp, + scomp, ncomp, slice_cr_ratio); + } + if( SliceType == WarpX::Ex_nodal_flag ) { + amrex::amrex_avgdown_edges(Dst_bx, Dst_fabox, Src_fabox, dcomp, + scomp, ncomp, slice_cr_ratio, 0); + } + if( SliceType == WarpX::Ey_nodal_flag) { + amrex::amrex_avgdown_edges(Dst_bx, Dst_fabox, Src_fabox, dcomp, + scomp, ncomp, slice_cr_ratio, 1); + } + if( SliceType == WarpX::Ez_nodal_flag ) { + amrex::amrex_avgdown_edges(Dst_bx, Dst_fabox, Src_fabox, dcomp, + scomp, ncomp, slice_cr_ratio, 2); + } + if( SliceType == WarpX::Bx_nodal_flag) { + amrex::amrex_avgdown_faces(Dst_bx, Dst_fabox, Src_fabox, dcomp, + scomp, ncomp, slice_cr_ratio, 0); + } + if( SliceType == WarpX::By_nodal_flag ) { + amrex::amrex_avgdown_faces(Dst_bx, Dst_fabox, Src_fabox, dcomp, + scomp, ncomp, slice_cr_ratio, 1); + } + if( SliceType == WarpX::Bz_nodal_flag ) { + amrex::amrex_avgdown_faces(Dst_bx, Dst_fabox, Src_fabox, dcomp, + scomp, ncomp, slice_cr_ratio, 2); + } + + if ( mfi_dst.isValid() ) { + ++mfi_dst; + } + + } + return cs_mf; + + } + amrex::Abort("Should not hit this return statement."); + return smf; +} + + +/* \brief + * This function modifies the slice input parameters under certain conditions. + * The coarsening ratio, slice_cr_ratio is modified if the input is not an exponent of 2. + * for example, if the coarsening ratio is 3, 5 or 6, which is not an exponent of 2, + * then the value of coarsening ratio is modified to the nearest exponent of 2. + * The default value for coarsening ratio is 1. + * slice_realbox.lo and slice_realbox.hi are set equal to the simulation domain lo and hi + * if for the user-input for the slice lo and hi coordinates are outside the domain. + * If the slice_realbox.lo and slice_realbox.hi coordinates do not align with the data + * points and the number of cells in that dimension is greater than 1, and if the extent of + * the slice in that dimension is not coarsenable, then the value lo and hi coordinates are + * shifted to the nearest coarsenable point to include some extra data points in the slice. + * If slice_realbox.lo==slice_realbox.hi, then that dimension has only one cell and no + * modifications are made to the value. If the lo and hi do not align with a data point, + * then it is flagged for interpolation. + * \param real_box a Real box defined for the underlying domain. + * \param slice_realbox a Real box for defining the slice dimension. + * \param slice_cc_nd_box a Real box for defining the modified lo and hi of the slice + * such that the coordinates align with the underlying data points. + * If the dimension is reduced to have only one cell, the slice_realbox is not modified and * instead the values are interpolated to the coordinate from the nearest data points. + * \param slice_cr_ratio contains values of the coarsening ratio which may be modified + * if the input values do not satisfy coarsenability conditions. + * \param slice_lo and slice_hi are the index values of the slice + * \param interp_lo are set to 0 or 1 if they are flagged for interpolation. + * The slice shares the same index space as that of the simulation domain. + */ + + +void +CheckSliceInput( const RealBox real_box, RealBox &slice_cc_nd_box, + RealBox &slice_realbox, IntVect &slice_cr_ratio, + Vector<Geometry> dom_geom, IntVect const SliceType, + IntVect &slice_lo, IntVect &slice_hi, IntVect &interp_lo) +{ + + IntVect slice_lo2(AMREX_D_DECL(0,0,0)); + for ( int idim = 0; idim < AMREX_SPACEDIM; ++idim) + { + // Modify coarsening ratio if the input value is not an exponent of 2 for AMR // + if ( slice_cr_ratio[idim] > 0 ) { + int log_cr_ratio = floor ( log2( double(slice_cr_ratio[idim]))); + slice_cr_ratio[idim] = exp2( double(log_cr_ratio) ); + } + + //// Default coarsening ratio is 1 // + // Modify lo if input is out of bounds // + if ( slice_realbox.lo(idim) < real_box.lo(idim) ) { + slice_realbox.setLo( idim, real_box.lo(idim)); + amrex::Print() << " slice lo is out of bounds. " << + " Modified it in dimension " << idim << + " to be aligned with the domain box\n"; + } + + // Modify hi if input in out od bounds // + if ( slice_realbox.hi(idim) > real_box.hi(idim) ) { + slice_realbox.setHi( idim, real_box.hi(idim)); + amrex::Print() << " slice hi is out of bounds." << + " Modified it in dimension " << idim << + " to be aligned with the domain box\n"; + } + + // Factor to ensure index values computation depending on index type // + double fac = ( 1.0 - SliceType[idim] )*dom_geom[0].CellSize(idim)*0.5; + // if dimension is reduced to one cell length // + if ( slice_realbox.hi(idim) - slice_realbox.lo(idim) <= 0) + { + slice_cc_nd_box.setLo( idim, slice_realbox.lo(idim) ); + slice_cc_nd_box.setHi( idim, slice_realbox.hi(idim) ); + + if ( slice_cr_ratio[idim] > 1) slice_cr_ratio[idim] = 1; + + // check for interpolation -- compute index lo with floor and ceil + if ( slice_cc_nd_box.lo(idim) - real_box.lo(idim) >= fac ) { + slice_lo[idim] = floor( ( (slice_cc_nd_box.lo(idim) + - (real_box.lo(idim) + fac ) ) + / dom_geom[0].CellSize(idim)) + fac * 1E-10); + slice_lo2[idim] = ceil( ( (slice_cc_nd_box.lo(idim) + - (real_box.lo(idim) + fac) ) + / dom_geom[0].CellSize(idim)) - fac * 1E-10 ); + } + else { + slice_lo[idim] = round( (slice_cc_nd_box.lo(idim) + - (real_box.lo(idim) ) ) + / dom_geom[0].CellSize(idim)); + slice_lo2[idim] = ceil((slice_cc_nd_box.lo(idim) + - (real_box.lo(idim) ) ) + / dom_geom[0].CellSize(idim) ); + } + + // flag for interpolation -- if reduced dimension location // + // does not align with data point // + if ( slice_lo[idim] == slice_lo2[idim]) { + if ( slice_cc_nd_box.lo(idim) - real_box.lo(idim) < fac ) { + interp_lo[idim] = 1; + } + } + else { + interp_lo[idim] = 1; + } + + // ncells = 1 if dimension is reduced // + slice_hi[idim] = slice_lo[idim] + 1; + } + else + { + // moving realbox.lo and reabox.hi to nearest coarsenable grid point // + int index_lo = floor(((slice_realbox.lo(idim) + 1E-10 + - (real_box.lo(idim))) / dom_geom[0].CellSize(idim))); + int index_hi = ceil(((slice_realbox.hi(idim) - 1E-10 + - (real_box.lo(idim))) / dom_geom[0].CellSize(idim))); + + bool modify_cr = true; + + while ( modify_cr == true) { + int lo_new = index_lo; + int hi_new = index_hi; + int mod_lo = index_lo % slice_cr_ratio[idim]; + int mod_hi = index_hi % slice_cr_ratio[idim]; + modify_cr = false; + + // To ensure that the index.lo is coarsenable // + if ( mod_lo > 0) { + lo_new = index_lo - mod_lo; + } + // To ensure that the index.hi is coarsenable // + if ( mod_hi > 0) { + hi_new = index_hi + (slice_cr_ratio[idim] - mod_hi); + } + + // If modified index.hi is > baselinebox.hi, move the point // + // to the previous coarsenable point // + if ( (hi_new * dom_geom[0].CellSize(idim)) + > real_box.hi(idim) - real_box.lo(idim) + dom_geom[0].CellSize(idim)*0.01 ) + { + hi_new = index_hi - mod_hi; + } + + if ( (hi_new - lo_new) == 0 ){ + amrex::Print() << " Diagnostic Warning :: "; + amrex::Print() << " Coarsening ratio "; + amrex::Print() << slice_cr_ratio[idim] << " in dim "<< idim; + amrex::Print() << "is leading to zero cells for slice."; + amrex::Print() << " Thus reducing cr_ratio by half.\n"; + + slice_cr_ratio[idim] = slice_cr_ratio[idim]/2; + modify_cr = true; + } + + if ( modify_cr == false ) { + index_lo = lo_new; + index_hi = hi_new; + } + slice_lo[idim] = index_lo; + slice_hi[idim] = index_hi - 1; // since default is cell-centered + } + slice_realbox.setLo( idim, index_lo * dom_geom[0].CellSize(idim) + + real_box.lo(idim) ); + slice_realbox.setHi( idim, index_hi * dom_geom[0].CellSize(idim) + + real_box.lo(idim) ); + slice_cc_nd_box.setLo( idim, slice_realbox.lo(idim) + fac ); + slice_cc_nd_box.setHi( idim, slice_realbox.hi(idim) - fac ); + } + } +} + + +/* \brief + * This function is called if the coordinates of the slice do not align with data points + * \param interp_lo is an IntVect which is flagged as 1, if interpolation + is required in the dimension. + */ +void +InterpolateSliceValues(MultiFab& smf, IntVect interp_lo, RealBox slice_realbox, + Vector<Geometry> geom, int ncomp, int nghost, + IntVect slice_lo, IntVect slice_hi, IntVect SliceType, + const RealBox real_box) +{ + for (MFIter mfi(smf); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const auto IndType = smf.ixType(); + const auto lo = amrex::lbound(bx); + const auto hi = amrex::ubound(bx); + FArrayBox& fabox = smf[mfi]; + + for ( int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + if ( interp_lo[idim] == 1 ) { + InterpolateLo( bx, fabox, slice_lo, geom, idim, SliceType, + slice_realbox, 0, ncomp, nghost, real_box); + } + } + } + +} + +void +InterpolateLo(const Box& bx, FArrayBox &fabox, IntVect slice_lo, + Vector<Geometry> geom, int idir, IntVect IndType, + RealBox slice_realbox, int srccomp, int ncomp, + int nghost, const RealBox real_box ) +{ + auto fabarr = fabox.array(); + const auto lo = amrex::lbound(bx); + const auto hi = amrex::ubound(bx); + double fac = ( 1.0-IndType[idir] )*geom[0].CellSize(idir) * 0.5; + int imin = slice_lo[idir]; + double minpos = imin*geom[0].CellSize(idir) + fac + real_box.lo(idir); + double maxpos = (imin+1)*geom[0].CellSize(idir) + fac + real_box.lo(idir); + double slice_minpos = slice_realbox.lo(idir) ; + + switch (idir) { + case 0: + { + if ( imin >= lo.x && imin <= lo.x) { + for (int n = srccomp; n < srccomp + ncomp; ++n) { + for (int k = lo.z; k <= hi.z; ++k) { + for (int j = lo.y; j <= hi.y; ++j) { + for (int i = lo.x; i <= hi.x; ++i) { + double minval = fabarr(i,j,k,n); + double maxval = fabarr(i+1,j,k,n); + double ratio = (maxval - minval) / (maxpos - minpos); + double xdiff = slice_minpos - minpos; + double newval = minval + xdiff * ratio; + fabarr(i,j,k,n) = newval; + } + } + } + } + } + break; + } + case 1: + { + if ( imin >= lo.y && imin <= lo.y) { + for (int n = srccomp; n < srccomp+ncomp; ++n) { + for (int k = lo.z; k <= hi.z; ++k) { + for (int j = lo.y; j <= hi.y; ++j) { + for (int i = lo.x; i <= hi.x; ++i) { + double minval = fabarr(i,j,k,n); + double maxval = fabarr(i,j+1,k,n); + double ratio = (maxval - minval) / (maxpos - minpos); + double xdiff = slice_minpos - minpos; + double newval = minval + xdiff * ratio; + fabarr(i,j,k,n) = newval; + } + } + } + } + } + break; + } + case 2: + { + if ( imin >= lo.z && imin <= lo.z) { + for (int n = srccomp; n < srccomp+ncomp; ++n) { + for (int k = lo.z; k <= hi.z; ++k) { + for (int j = lo.y; j <= hi.y; ++j) { + for (int i = lo.x; i <= hi.x; ++i) { + double minval = fabarr(i,j,k,n); + double maxval = fabarr(i,j,k+1,n); + double ratio = (maxval - minval) / (maxpos - minpos); + double xdiff = slice_minpos - minpos; + double newval = minval + xdiff * ratio; + fabarr(i,j,k,n) = newval; + } + } + } + } + } + break; + } + + } + +} + + + + + + + diff --git a/Source/Diagnostics/WarpXIO.cpp b/Source/Diagnostics/WarpXIO.cpp index 6d646dd42..7e7ed278d 100644 --- a/Source/Diagnostics/WarpXIO.cpp +++ b/Source/Diagnostics/WarpXIO.cpp @@ -11,6 +11,8 @@ #include <AMReX_AmrMeshInSituBridge.H> #endif +#include "SliceDiagnostic.H" + #ifdef AMREX_USE_ASCENT #include <ascent.hpp> #include <AMReX_Conduit_Blueprint.H> @@ -771,3 +773,94 @@ WarpX::WriteJobInfo (const std::string& dir) const jobInfoFile.close(); } } + + +/* \brief + * The slice is ouput using visMF and can be visualized used amrvis. + */ +void +WarpX::WriteSlicePlotFile () const +{ + if (F_fp[0] ) { + VisMF::Write( (*F_slice[0]), "vismf_F_slice"); + } + + if (rho_fp[0]) { + VisMF::Write( (*rho_slice[0]), "vismf_rho_slice"); + } + + VisMF::Write( (*Efield_slice[0][0]), amrex::Concatenate("vismf_Ex_slice_",istep[0])); + VisMF::Write( (*Efield_slice[0][1]), amrex::Concatenate("vismf_Ey_slice_",istep[0])); + VisMF::Write( (*Efield_slice[0][2]), amrex::Concatenate("vismf_Ez_slice_",istep[0])); + VisMF::Write( (*Bfield_slice[0][0]), amrex::Concatenate("vismf_Bx_slice_",istep[0])); + VisMF::Write( (*Bfield_slice[0][1]), amrex::Concatenate("vismf_By_slice_",istep[0])); + VisMF::Write( (*Bfield_slice[0][2]), amrex::Concatenate("vismf_Bz_slice_",istep[0])); + VisMF::Write( (*current_slice[0][0]), amrex::Concatenate("vismf_jx_slice_",istep[0])); + VisMF::Write( (*current_slice[0][1]), amrex::Concatenate("vismf_jy_slice_",istep[0])); + VisMF::Write( (*current_slice[0][2]), amrex::Concatenate("vismf_jz_slice_",istep[0])); + +} + + +void +WarpX::InitializeSliceMultiFabs () +{ + + int nlevels = Geom().size(); + + F_slice.resize(nlevels); + rho_slice.resize(nlevels); + current_slice.resize(nlevels); + Efield_slice.resize(nlevels); + Bfield_slice.resize(nlevels); + +} + + +// To generate slice that inherits index type of underlying data // +void +WarpX::SliceGenerationForDiagnostics () +{ + + Vector<Geometry> dom_geom; + dom_geom = Geom(); + + if (F_fp[0] ) { + F_slice[0] = CreateSlice( *F_fp[0].get(), dom_geom, slice_realbox, + slice_cr_ratio ); + } + if (rho_fp[0]) { + rho_slice[0] = CreateSlice( *rho_fp[0].get(), dom_geom, slice_realbox, + slice_cr_ratio ); + } + + for (int idim = 0; idim < 3; ++idim) { + Efield_slice[0][idim] = CreateSlice( *Efield_fp[0][idim].get(), + dom_geom, slice_realbox, slice_cr_ratio ); + Bfield_slice[0][idim] = CreateSlice( *Bfield_fp[0][idim].get(), + dom_geom, slice_realbox, slice_cr_ratio ); + current_slice[0][idim] = CreateSlice( *current_fp[0][idim].get(), + dom_geom, slice_realbox, slice_cr_ratio ); + } + + +} + + +void +WarpX::ClearSliceMultiFabs () +{ + + F_slice.clear(); + rho_slice.clear(); + current_slice.clear(); + Efield_slice.clear(); + Bfield_slice.clear(); + F_slice.shrink_to_fit(); + rho_slice.shrink_to_fit(); + current_slice.shrink_to_fit(); + Efield_slice.shrink_to_fit(); + Bfield_slice.shrink_to_fit(); + +} + diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index ad7c7d840..a5d68e4f9 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -132,6 +132,8 @@ WarpX::EvolveEM (int numsteps) bool to_make_plot = (plot_int > 0) && ((step+1) % plot_int == 0); + // slice generation // + bool to_make_slice_plot = (slice_plot_int > 0) && ( (step+1)% slice_plot_int == 0); bool do_insitu = ((step+1) >= insitu_start) && (insitu_int > 0) && ((step+1) % insitu_int == 0); @@ -176,7 +178,8 @@ WarpX::EvolveEM (int numsteps) myBFD->writeLabFrameData(cell_centered_data.get(), *mypc, geom[0], cur_time, dt[0]); } - if (to_make_plot || do_insitu) + // slice gen // + if (to_make_plot || do_insitu || to_make_slice_plot) { FillBoundaryE(); FillBoundaryB(); @@ -192,11 +195,19 @@ WarpX::EvolveEM (int numsteps) last_insitu_step = step+1; if (to_make_plot) - WritePlotFile(); + WritePlotFile(); + + if (to_make_slice_plot) + { + InitializeSliceMultiFabs (); + SliceGenerationForDiagnostics(); + WriteSlicePlotFile(); + ClearSliceMultiFabs (); + } if (do_insitu) UpdateInSitu(); - } + } if (check_int > 0 && (step+1) % check_int == 0) { last_check_file_step = step+1; diff --git a/Source/Filter/BilinearFilter.H b/Source/Filter/BilinearFilter.H index 9bade3c77..4fcfec14b 100644 --- a/Source/Filter/BilinearFilter.H +++ b/Source/Filter/BilinearFilter.H @@ -1,30 +1,17 @@ -#include <AMReX_MultiFab.H> -#include <AMReX_CudaContainers.H> +#include <Filter.H> -#ifndef WARPX_FILTER_H_ -#define WARPX_FILTER_H_ +#ifndef WARPX_BILIN_FILTER_H_ +#define WARPX_BILIN_FILTER_H_ -class BilinearFilter +class BilinearFilter : public Filter { public: - BilinearFilter () = default; + + BilinearFilter() = default; void ComputeStencils(); - void ApplyStencil(amrex::MultiFab& dstmf, const amrex::MultiFab& srcmf, int scomp=0, int dcomp=0, int ncomp=10000); amrex::IntVect npass_each_dir; - amrex::IntVect stencil_length_each_dir; - - // public for cuda - void Filter(const amrex::Box& tbx, - amrex::Array4<amrex::Real const> const& tmp, - amrex::Array4<amrex::Real > const& dst, - int scomp, int dcomp, int ncomp); - -private: - - amrex::Gpu::ManagedVector<amrex::Real> stencil_x, stencil_y, stencil_z; - amrex::Dim3 slen; }; -#endif +#endif // #ifndef WARPX_BILIN_FILTER_H_ diff --git a/Source/Filter/BilinearFilter.cpp b/Source/Filter/BilinearFilter.cpp index f6acaa5df..68ebde687 100644 --- a/Source/Filter/BilinearFilter.cpp +++ b/Source/Filter/BilinearFilter.cpp @@ -68,173 +68,3 @@ void BilinearFilter::ComputeStencils(){ slen.z = 1; #endif } - - -#ifdef AMREX_USE_CUDA - -void -BilinearFilter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dcomp, int ncomp) -{ - BL_PROFILE("BilinearFilter::ApplyStencil()"); - ncomp = std::min(ncomp, srcmf.nComp()); - - for (MFIter mfi(dstmf); mfi.isValid(); ++mfi) - { - const auto& src = srcmf.array(mfi); - const auto& dst = dstmf.array(mfi); - const Box& tbx = mfi.growntilebox(); - const Box& gbx = amrex::grow(tbx,stencil_length_each_dir-1); - - // tmpfab has enough ghost cells for the stencil - FArrayBox tmp_fab(gbx,ncomp); - Elixir tmp_eli = tmp_fab.elixir(); // Prevent the tmp data from being deleted too early - auto const& tmp = tmp_fab.array(); - - // Copy values in srcfab into tmpfab - const Box& ibx = gbx & srcmf[mfi].box(); - AMREX_PARALLEL_FOR_4D ( gbx, ncomp, i, j, k, n, - { - if (ibx.contains(IntVect(AMREX_D_DECL(i,j,k)))) { - tmp(i,j,k,n) = src(i,j,k,n+scomp); - } else { - tmp(i,j,k,n) = 0.0; - } - }); - - // Apply filter - Filter(tbx, tmp, dst, 0, dcomp, ncomp); - } -} - -void BilinearFilter::Filter (const Box& tbx, - Array4<Real const> const& tmp, - Array4<Real > const& dst, - int scomp, int dcomp, int ncomp) -{ - amrex::Real const* AMREX_RESTRICT sx = stencil_x.data(); - amrex::Real const* AMREX_RESTRICT sy = stencil_y.data(); - amrex::Real const* AMREX_RESTRICT sz = stencil_z.data(); - Dim3 slen_local = slen; - AMREX_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n, - { - Real d = 0.0; - - for (int iz=0; iz < slen_local.z; ++iz){ - for (int iy=0; iy < slen_local.y; ++iy){ - for (int ix=0; ix < slen_local.x; ++ix){ -#if (AMREX_SPACEDIM == 3) - Real sss = sx[ix]*sy[iy]*sz[iz]; -#else - Real sss = sx[ix]*sz[iy]; -#endif -#if (AMREX_SPACEDIM == 3) - d += sss*( tmp(i-ix,j-iy,k-iz,scomp+n) - +tmp(i+ix,j-iy,k-iz,scomp+n) - +tmp(i-ix,j+iy,k-iz,scomp+n) - +tmp(i+ix,j+iy,k-iz,scomp+n) - +tmp(i-ix,j-iy,k+iz,scomp+n) - +tmp(i+ix,j-iy,k+iz,scomp+n) - +tmp(i-ix,j+iy,k+iz,scomp+n) - +tmp(i+ix,j+iy,k+iz,scomp+n)); -#else - d += sss*( tmp(i-ix,j-iy,k,scomp+n) - +tmp(i+ix,j-iy,k,scomp+n) - +tmp(i-ix,j+iy,k,scomp+n) - +tmp(i+ix,j+iy,k,scomp+n)); -#endif - } - } - } - - dst(i,j,k,dcomp+n) = d; - }); -} - -#else - -void -BilinearFilter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dcomp, int ncomp) -{ - BL_PROFILE("BilinearFilter::ApplyStencil()"); - ncomp = std::min(ncomp, srcmf.nComp()); -#ifdef _OPENMP -#pragma omp parallel -#endif - { - FArrayBox tmpfab; - for (MFIter mfi(dstmf,true); mfi.isValid(); ++mfi){ - const auto& srcfab = srcmf[mfi]; - auto& dstfab = dstmf[mfi]; - const Box& tbx = mfi.growntilebox(); - const Box& gbx = amrex::grow(tbx,stencil_length_each_dir-1); - // tmpfab has enough ghost cells for the stencil - tmpfab.resize(gbx,ncomp); - tmpfab.setVal(0.0, gbx, 0, ncomp); - // Copy values in srcfab into tmpfab - const Box& ibx = gbx & srcfab.box(); - tmpfab.copy(srcfab, ibx, scomp, ibx, 0, ncomp); - // Apply filter - Filter(tbx, tmpfab.array(), dstfab.array(), 0, dcomp, ncomp); - } - } -} - -void BilinearFilter::Filter (const Box& tbx, - Array4<Real const> const& tmp, - Array4<Real > const& dst, - int scomp, int dcomp, int ncomp) -{ - const auto lo = amrex::lbound(tbx); - const auto hi = amrex::ubound(tbx); - // tmp and dst are of type Array4 (Fortran ordering) - amrex::Real const* AMREX_RESTRICT sx = stencil_x.data(); - amrex::Real const* AMREX_RESTRICT sy = stencil_y.data(); - amrex::Real const* AMREX_RESTRICT sz = stencil_z.data(); - for (int n = 0; n < ncomp; ++n) { - // Set dst value to 0. - for (int k = lo.z; k <= hi.z; ++k) { - for (int j = lo.y; j <= hi.y; ++j) { - for (int i = lo.x; i <= hi.x; ++i) { - dst(i,j,k,dcomp+n) = 0.0; - } - } - } - // 3 nested loop on 3D stencil - for (int iz=0; iz < slen.z; ++iz){ - for (int iy=0; iy < slen.y; ++iy){ - for (int ix=0; ix < slen.x; ++ix){ -#if (AMREX_SPACEDIM == 3) - Real sss = sx[ix]*sy[iy]*sz[iz]; -#else - Real sss = sx[ix]*sz[iy]; -#endif - // 3 nested loop on 3D array - for (int k = lo.z; k <= hi.z; ++k) { - for (int j = lo.y; j <= hi.y; ++j) { - AMREX_PRAGMA_SIMD - for (int i = lo.x; i <= hi.x; ++i) { -#if (AMREX_SPACEDIM == 3) - dst(i,j,k,dcomp+n) += sss*(tmp(i-ix,j-iy,k-iz,scomp+n) - +tmp(i+ix,j-iy,k-iz,scomp+n) - +tmp(i-ix,j+iy,k-iz,scomp+n) - +tmp(i+ix,j+iy,k-iz,scomp+n) - +tmp(i-ix,j-iy,k+iz,scomp+n) - +tmp(i+ix,j-iy,k+iz,scomp+n) - +tmp(i-ix,j+iy,k+iz,scomp+n) - +tmp(i+ix,j+iy,k+iz,scomp+n)); -#else - dst(i,j,k,dcomp+n) += sss*(tmp(i-ix,j-iy,k,scomp+n) - +tmp(i+ix,j-iy,k,scomp+n) - +tmp(i-ix,j+iy,k,scomp+n) - +tmp(i+ix,j+iy,k,scomp+n)); -#endif - } - } - } - } - } - } - } -} - -#endif diff --git a/Source/Filter/Filter.H b/Source/Filter/Filter.H new file mode 100644 index 000000000..eaa0498c9 --- /dev/null +++ b/Source/Filter/Filter.H @@ -0,0 +1,43 @@ +#include <AMReX_MultiFab.H> +#include <AMReX_CudaContainers.H> + +#ifndef WARPX_FILTER_H_ +#define WARPX_FILTER_H_ + +class Filter +{ +public: + Filter () = default; + + // Apply stencil on MultiFab. + // Guard cells are handled inside this function + void ApplyStencil(amrex::MultiFab& dstmf, + const amrex::MultiFab& srcmf, int scomp=0, + int dcomp=0, int ncomp=10000); + + // Apply stencil on a FabArray. + void ApplyStencil (amrex::FArrayBox& dstfab, + const amrex::FArrayBox& srcfab, const amrex::Box& tbx, + int scomp=0, int dcomp=0, int ncomp=10000); + + // public for cuda + void DoFilter(const amrex::Box& tbx, + amrex::Array4<amrex::Real const> const& tmp, + amrex::Array4<amrex::Real > const& dst, + int scomp, int dcomp, int ncomp); + + // In 2D, stencil_length_each_dir = {length(stencil_x), length(stencil_z)} + amrex::IntVect stencil_length_each_dir; + +protected: + // Stencil along each direction. + // in 2D, stencil_y is not initialized. + amrex::Gpu::ManagedVector<amrex::Real> stencil_x, stencil_y, stencil_z; + // Length of each stencil. + // In 2D, slen = {length(stencil_x), length(stencil_z), 1} + amrex::Dim3 slen; + +private: + +}; +#endif // #ifndef WARPX_FILTER_H_ diff --git a/Source/Filter/Filter.cpp b/Source/Filter/Filter.cpp new file mode 100644 index 000000000..f8a2dd6c2 --- /dev/null +++ b/Source/Filter/Filter.cpp @@ -0,0 +1,257 @@ +#include <WarpX.H> +#include <Filter.H> + +#ifdef _OPENMP +#include <omp.h> +#endif + +using namespace amrex; + +#ifdef AMREX_USE_CUDA + +/* \brief Apply stencil on MultiFab (GPU version, 2D/3D). + * \param dstmf Destination MultiFab + * \param srcmf source MultiFab + * \param scomp first component of srcmf on which the filter is applied + * \param dcomp first component of dstmf on which the filter is applied + * \param ncomp Number of components on which the filter is applied. + */ +void +Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dcomp, int ncomp) +{ + BL_PROFILE("BilinearFilter::ApplyStencil(MultiFab)"); + ncomp = std::min(ncomp, srcmf.nComp()); + + for (MFIter mfi(dstmf); mfi.isValid(); ++mfi) + { + const auto& src = srcmf.array(mfi); + const auto& dst = dstmf.array(mfi); + const Box& tbx = mfi.growntilebox(); + const Box& gbx = amrex::grow(tbx,stencil_length_each_dir-1); + + // tmpfab has enough ghost cells for the stencil + FArrayBox tmp_fab(gbx,ncomp); + Elixir tmp_eli = tmp_fab.elixir(); // Prevent the tmp data from being deleted too early + auto const& tmp = tmp_fab.array(); + + // Copy values in srcfab into tmpfab + const Box& ibx = gbx & srcmf[mfi].box(); + AMREX_PARALLEL_FOR_4D ( gbx, ncomp, i, j, k, n, + { + if (ibx.contains(IntVect(AMREX_D_DECL(i,j,k)))) { + tmp(i,j,k,n) = src(i,j,k,n+scomp); + } else { + tmp(i,j,k,n) = 0.0; + } + }); + + // Apply filter + DoFilter(tbx, tmp, dst, 0, dcomp, ncomp); + } +} + +/* \brief Apply stencil on FArrayBox (GPU version, 2D/3D). + * \param dstfab Destination FArrayBox + * \param srcmf source FArrayBox + * \param tbx Grown box on which srcfab is defined. + * \param scomp first component of srcfab on which the filter is applied + * \param dcomp first component of dstfab on which the filter is applied + * \param ncomp Number of components on which the filter is applied. + */ +void +Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, + const Box& tbx, int scomp, int dcomp, int ncomp) +{ + BL_PROFILE("BilinearFilter::ApplyStencil(FArrayBox)"); + ncomp = std::min(ncomp, srcfab.nComp()); + const auto& src = srcfab.array(); + const auto& dst = dstfab.array(); + const Box& gbx = amrex::grow(tbx,stencil_length_each_dir-1); + + // tmpfab has enough ghost cells for the stencil + FArrayBox tmp_fab(gbx,ncomp); + Elixir tmp_eli = tmp_fab.elixir(); // Prevent the tmp data from being deleted too early + auto const& tmp = tmp_fab.array(); + + // Copy values in srcfab into tmpfab + const Box& ibx = gbx & srcfab.box(); + AMREX_PARALLEL_FOR_4D ( gbx, ncomp, i, j, k, n, + { + if (ibx.contains(IntVect(AMREX_D_DECL(i,j,k)))) { + tmp(i,j,k,n) = src(i,j,k,n+scomp); + } else { + tmp(i,j,k,n) = 0.0; + } + }); + + // Apply filter + DoFilter(tbx, tmp, dst, 0, dcomp, ncomp); +} + +/* \brief Apply stencil (2D/3D, CPU/GPU) + */ +void Filter::DoFilter (const Box& tbx, + Array4<Real const> const& tmp, + Array4<Real > const& dst, + int scomp, int dcomp, int ncomp) +{ + amrex::Real const* AMREX_RESTRICT sx = stencil_x.data(); + amrex::Real const* AMREX_RESTRICT sy = stencil_y.data(); + amrex::Real const* AMREX_RESTRICT sz = stencil_z.data(); + Dim3 slen_local = slen; + AMREX_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n, + { + Real d = 0.0; + + for (int iz=0; iz < slen_local.z; ++iz){ + for (int iy=0; iy < slen_local.y; ++iy){ + for (int ix=0; ix < slen_local.x; ++ix){ +#if (AMREX_SPACEDIM == 3) + Real sss = sx[ix]*sy[iy]*sz[iz]; +#else + Real sss = sx[ix]*sz[iy]; +#endif +#if (AMREX_SPACEDIM == 3) + d += sss*( tmp(i-ix,j-iy,k-iz,scomp+n) + +tmp(i+ix,j-iy,k-iz,scomp+n) + +tmp(i-ix,j+iy,k-iz,scomp+n) + +tmp(i+ix,j+iy,k-iz,scomp+n) + +tmp(i-ix,j-iy,k+iz,scomp+n) + +tmp(i+ix,j-iy,k+iz,scomp+n) + +tmp(i-ix,j+iy,k+iz,scomp+n) + +tmp(i+ix,j+iy,k+iz,scomp+n)); +#else + d += sss*( tmp(i-ix,j-iy,k,scomp+n) + +tmp(i+ix,j-iy,k,scomp+n) + +tmp(i-ix,j+iy,k,scomp+n) + +tmp(i+ix,j+iy,k,scomp+n)); +#endif + } + } + } + + dst(i,j,k,dcomp+n) = d; + }); +} + +#else + +/* \brief Apply stencil on MultiFab (CPU version, 2D/3D). + * \param dstmf Destination MultiFab + * \param srcmf source MultiFab + * \param scomp first component of srcmf on which the filter is applied + * \param dcomp first component of dstmf on which the filter is applied + * \param ncomp Number of components on which the filter is applied. + */ +void +Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dcomp, int ncomp) +{ + BL_PROFILE("BilinearFilter::ApplyStencil()"); + ncomp = std::min(ncomp, srcmf.nComp()); +#ifdef _OPENMP +#pragma omp parallel +#endif + { + FArrayBox tmpfab; + for (MFIter mfi(dstmf,true); mfi.isValid(); ++mfi){ + const auto& srcfab = srcmf[mfi]; + auto& dstfab = dstmf[mfi]; + const Box& tbx = mfi.growntilebox(); + const Box& gbx = amrex::grow(tbx,stencil_length_each_dir-1); + // tmpfab has enough ghost cells for the stencil + tmpfab.resize(gbx,ncomp); + tmpfab.setVal(0.0, gbx, 0, ncomp); + // Copy values in srcfab into tmpfab + const Box& ibx = gbx & srcfab.box(); + tmpfab.copy(srcfab, ibx, scomp, ibx, 0, ncomp); + // Apply filter + DoFilter(tbx, tmpfab.array(), dstfab.array(), 0, dcomp, ncomp); + } + } +} + +/* \brief Apply stencil on FArrayBox (CPU version, 2D/3D). + * \param dstfab Destination FArrayBox + * \param srcmf source FArrayBox + * \param tbx Grown box on which srcfab is defined. + * \param scomp first component of srcfab on which the filter is applied + * \param dcomp first component of dstfab on which the filter is applied + * \param ncomp Number of components on which the filter is applied. + */ +void +Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, + const Box& tbx, int scomp, int dcomp, int ncomp) +{ + BL_PROFILE("BilinearFilter::ApplyStencil(FArrayBox)"); + ncomp = std::min(ncomp, srcfab.nComp()); + FArrayBox tmpfab; + const Box& gbx = amrex::grow(tbx,stencil_length_each_dir-1); + // tmpfab has enough ghost cells for the stencil + tmpfab.resize(gbx,ncomp); + tmpfab.setVal(0.0, gbx, 0, ncomp); + // Copy values in srcfab into tmpfab + const Box& ibx = gbx & srcfab.box(); + tmpfab.copy(srcfab, ibx, scomp, ibx, 0, ncomp); + // Apply filter + DoFilter(tbx, tmpfab.array(), dstfab.array(), 0, dcomp, ncomp); +} + +void Filter::DoFilter (const Box& tbx, + Array4<Real const> const& tmp, + Array4<Real > const& dst, + int scomp, int dcomp, int ncomp) +{ + const auto lo = amrex::lbound(tbx); + const auto hi = amrex::ubound(tbx); + // tmp and dst are of type Array4 (Fortran ordering) + amrex::Real const* AMREX_RESTRICT sx = stencil_x.data(); + amrex::Real const* AMREX_RESTRICT sy = stencil_y.data(); + amrex::Real const* AMREX_RESTRICT sz = stencil_z.data(); + for (int n = 0; n < ncomp; ++n) { + // Set dst value to 0. + for (int k = lo.z; k <= hi.z; ++k) { + for (int j = lo.y; j <= hi.y; ++j) { + for (int i = lo.x; i <= hi.x; ++i) { + dst(i,j,k,dcomp+n) = 0.0; + } + } + } + // 3 nested loop on 3D stencil + for (int iz=0; iz < slen.z; ++iz){ + for (int iy=0; iy < slen.y; ++iy){ + for (int ix=0; ix < slen.x; ++ix){ +#if (AMREX_SPACEDIM == 3) + Real sss = sx[ix]*sy[iy]*sz[iz]; +#else + Real sss = sx[ix]*sz[iy]; +#endif + // 3 nested loop on 3D array + for (int k = lo.z; k <= hi.z; ++k) { + for (int j = lo.y; j <= hi.y; ++j) { + AMREX_PRAGMA_SIMD + for (int i = lo.x; i <= hi.x; ++i) { +#if (AMREX_SPACEDIM == 3) + dst(i,j,k,dcomp+n) += sss*(tmp(i-ix,j-iy,k-iz,scomp+n) + +tmp(i+ix,j-iy,k-iz,scomp+n) + +tmp(i-ix,j+iy,k-iz,scomp+n) + +tmp(i+ix,j+iy,k-iz,scomp+n) + +tmp(i-ix,j-iy,k+iz,scomp+n) + +tmp(i+ix,j-iy,k+iz,scomp+n) + +tmp(i-ix,j+iy,k+iz,scomp+n) + +tmp(i+ix,j+iy,k+iz,scomp+n)); +#else + dst(i,j,k,dcomp+n) += sss*(tmp(i-ix,j-iy,k,scomp+n) + +tmp(i+ix,j-iy,k,scomp+n) + +tmp(i-ix,j+iy,k,scomp+n) + +tmp(i+ix,j+iy,k,scomp+n)); +#endif + } + } + } + } + } + } + } +} + +#endif // #ifdef AMREX_USE_CUDA diff --git a/Source/Filter/Make.package b/Source/Filter/Make.package index 36e74bfb0..8b8e9b50b 100644 --- a/Source/Filter/Make.package +++ b/Source/Filter/Make.package @@ -1,5 +1,9 @@ +CEXE_headers += Filter.H +CEXE_sources += Filter.cpp CEXE_sources += BilinearFilter.cpp CEXE_headers += BilinearFilter.H +CEXE_sources += NCIGodfreyFilter.cpp +CEXE_headers += NCIGodfreyFilter.H INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Filter VPATH_LOCATIONS += $(WARPX_HOME)/Source/Filter diff --git a/Source/Filter/NCIGodfreyFilter.H b/Source/Filter/NCIGodfreyFilter.H new file mode 100644 index 000000000..a53039dfa --- /dev/null +++ b/Source/Filter/NCIGodfreyFilter.H @@ -0,0 +1,30 @@ +#include <Filter.H> + +#ifndef WARPX_GODFREY_FILTER_H_ +#define WARPX_GODFREY_FILTER_H_ + +enum class godfrey_coeff_set { Ex_Ey_Bz=0, Bx_By_Ez=1 }; + +class NCIGodfreyFilter : public Filter +{ +public: + + NCIGodfreyFilter () = default; + + NCIGodfreyFilter(godfrey_coeff_set coeff_set_, amrex::Real cdtodz_, amrex::Real l_lower_order_in_v_); + + void ComputeStencils(); + + void getGodfreyCoeffs(godfrey_coeff_set coeff_set_in); + + static constexpr int stencil_width = 4; + +private: + + godfrey_coeff_set coeff_set; + amrex::Real cdtodz; + int l_lower_order_in_v; + +}; + +#endif // #ifndef WARPX_GODFREY_FILTER_H_ diff --git a/Source/Filter/NCIGodfreyFilter.cpp b/Source/Filter/NCIGodfreyFilter.cpp new file mode 100644 index 000000000..3f589260a --- /dev/null +++ b/Source/Filter/NCIGodfreyFilter.cpp @@ -0,0 +1,78 @@ +#include <WarpX.H> +#include <NCIGodfreyFilter.H> +#include <NCIGodfreyTables.H> + +#ifdef _OPENMP +#include <omp.h> +#endif + +using namespace amrex; + +NCIGodfreyFilter::NCIGodfreyFilter(godfrey_coeff_set coeff_set_, amrex::Real cdtodz_, amrex::Real l_lower_order_in_v_){ + // Store parameters into class data members + coeff_set = coeff_set_; + cdtodz = cdtodz_; + l_lower_order_in_v = l_lower_order_in_v_; + // NCI Godfrey filter has fixed size, and is applied along z only. +#if (AMREX_SPACEDIM == 3) + stencil_length_each_dir = {1,1,5}; + slen = {1,1,5}; +#else + stencil_length_each_dir = {1,5}; + slen = {1,5,1}; +#endif +} + +void NCIGodfreyFilter::ComputeStencils(){ + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( +#if ( AMREX_SPACEDIM == 3 ) + slen.z==5, +#else + slen.y==5, +#endif + "ERROR: NCI filter requires 5 points in z"); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(l_lower_order_in_v==1, + "ERROR: NCI corrector requires l_lower_order_in_v=1, i.e., Galerkin scheme"); + + // Interpolate coefficients from the table, and store into prestencil. + int index = tab_length*cdtodz; + index = min(index, tab_length-2); + index = max(index, 0); + Real weight_right = cdtodz - index/tab_length; + Real prestencil[4]; + for(int i=0; i<tab_width; i++){ + if (coeff_set == godfrey_coeff_set::Ex_Ey_Bz) { + prestencil[i] = (1-weight_right)*table_nci_godfrey_Ex_Ey_Bz[index ][i] + + weight_right*table_nci_godfrey_Ex_Ey_Bz[index+1][i]; + } else if (coeff_set == godfrey_coeff_set::Bx_By_Ez) { + prestencil[i] = (1-weight_right)*table_nci_godfrey_Bx_By_Ez[index ][i] + + weight_right*table_nci_godfrey_Bx_By_Ez[index+1][i]; + } + } + + // Compute stencil_z + stencil_z.resize( 5 ); + stencil_z[0] = (256 + 128*prestencil[0] + 96*prestencil[1] + 80*prestencil[2] + 70*prestencil[3]) / 256; + stencil_z[1] = -( 64*prestencil[0] + 64*prestencil[1] + 60*prestencil[2] + 56*prestencil[3]) / 256; + stencil_z[2] = ( 16*prestencil[1] + 24*prestencil[2] + 28*prestencil[3]) / 256; + stencil_z[3] = -( 4*prestencil[2] + 8*prestencil[3]) / 256; + stencil_z[4] = ( 1*prestencil[3]) / 256; + + // Compute stencil_x and stencil_y (no filter in these directions, + // so only 1 coeff, equal to 1) + stencil_x.resize(1); + stencil_x[0] = 1.; +#if (AMREX_SPACEDIM == 3) + stencil_y.resize(1); + stencil_y[0] = 1.; +#endif + + // Due to the way Filter::DoFilter() is written, + // coefficient 0 has to be /2 + stencil_x[0] /= 2.; +#if (AMREX_SPACEDIM == 3) + stencil_y[0] /= 2.; +#endif + stencil_z[0] /= 2.; + +} diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H index 52996a60a..1053ace89 100644 --- a/Source/FortranInterface/WarpX_f.H +++ b/Source/FortranInterface/WarpX_f.H @@ -39,10 +39,6 @@ #define WRPX_PUSH_LEAPFROG warpx_push_leapfrog_3d #define WRPX_PUSH_LEAPFROG_POSITIONS warpx_push_leapfrog_positions_3d -#define WRPX_PXR_NCI_CORR_INIT init_godfrey_filter_coeffs -#define WRPX_PXR_GODFREY_FILTER apply_godfrey_filter_z_3d - - #elif (AMREX_SPACEDIM == 2) #define WRPX_COMPUTE_DIVB warpx_compute_divb_2d @@ -66,9 +62,6 @@ #define WRPX_PUSH_LEAPFROG warpx_push_leapfrog_2d #define WRPX_PUSH_LEAPFROG_POSITIONS warpx_push_leapfrog_positions_2d -#define WRPX_PXR_NCI_CORR_INIT init_godfrey_filter_coeffs -#define WRPX_PXR_GODFREY_FILTER apply_godfrey_filter_z_2d - #ifdef WARPX_RZ #define WRPX_COMPUTE_DIVE warpx_compute_dive_rz #else @@ -455,14 +448,6 @@ extern "C" const BL_FORT_FAB_ARG_ANYD(fine), const int* ncomp); - void WRPX_PXR_NCI_CORR_INIT(amrex::Real*, amrex::Real*, const int, - const amrex::Real, const int); - - void WRPX_PXR_GODFREY_FILTER (const int* lo, const int* hi, - amrex_real* fou, const int* olo, const int* ohi, - const amrex_real* fin, const int* ilo, const int* ihi, - const amrex_real* stencil, const int* nsten); - #ifdef WARPX_USE_PSATD void warpx_fft_mpi_init (int fcomm); void warpx_fft_domain_decomp (int* warpx_local_nz, int* warpx_local_z0, diff --git a/Source/FortranInterface/WarpX_picsar.F90 b/Source/FortranInterface/WarpX_picsar.F90 index c17e8861b..12d541b08 100644 --- a/Source/FortranInterface/WarpX_picsar.F90 +++ b/Source/FortranInterface/WarpX_picsar.F90 @@ -117,7 +117,7 @@ contains pxr_ll4symtry = ll4symtry .eq. 1 pxr_l_lower_order_in_v = l_lower_order_in_v .eq. 1 pxr_l_nodal = l_nodal .eq. 1 - + exg_nguards = ixyzmin - exg_lo eyg_nguards = ixyzmin - eyg_lo ezg_nguards = ixyzmin - ezg_lo @@ -218,6 +218,11 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n IF ((nox.eq.1).and.(noy.eq.1).and.(noz.eq.1)) THEN CALL depose_rho_vecHVv2_1_1_1(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,nx,ny,nz,& nxguard,nyguard,nzguard,lvect) + + ELSE IF ((nox.eq.2).and.(noy.eq.2).and.(noz.eq.2)) THEN + CALL depose_rho_vecHVv2_2_2_2(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,nx,ny,nz,& + nxguard,nyguard,nzguard,lvect) + ELSE CALL pxr_depose_rho_n(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,nx,ny,nz,& nxguard,nyguard,nzguard,nox,noy,noz, & @@ -268,7 +273,7 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n #ifdef WARPX_RZ integer(c_long) :: type_rz_depose = 1 #endif - + ! Compute the number of valid cells and guard cells integer(c_long) :: rho_nvalid(AMREX_SPACEDIM), rho_nguards(AMREX_SPACEDIM) rho_nvalid = rho_ntot - 2*rho_ng @@ -387,7 +392,7 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n real(amrex_real), intent(IN OUT):: jx(*), jy(*), jz(*) real(amrex_real), intent(IN) :: rmin, dr -#ifdef WARPX_RZ +#ifdef WARPX_RZ integer(c_long) :: type_rz_depose = 1 #endif ! Compute the number of valid cells and guard cells diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 23637ec97..0f33d1a0f 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -7,6 +7,7 @@ #include <WarpX.H> #include <WarpX_f.H> #include <BilinearFilter.H> +#include <NCIGodfreyFilter.H> #ifdef BL_USE_SENSEI_INSITU #include <AMReX_AmrMeshInSituBridge.H> @@ -161,8 +162,6 @@ WarpX::InitNCICorrector () { if (WarpX::use_fdtd_nci_corr) { - mypc->fdtd_nci_stencilz_ex.resize(max_level+1); - mypc->fdtd_nci_stencilz_by.resize(max_level+1); for (int lev = 0; lev <= max_level; ++lev) { const Geometry& gm = Geom(lev); @@ -174,10 +173,15 @@ WarpX::InitNCICorrector () dz = dx[1]; } cdtodz = PhysConst::c * dt[lev] / dz; - WRPX_PXR_NCI_CORR_INIT( (mypc->fdtd_nci_stencilz_ex)[lev].data(), - (mypc->fdtd_nci_stencilz_by)[lev].data(), - mypc->nstencilz_fdtd_nci_corr, cdtodz, - WarpX::l_lower_order_in_v); + + // Initialize Godfrey filters + // Same filter for fields Ex, Ey and Bz + nci_godfrey_filter_exeybz[lev].reset( new NCIGodfreyFilter(godfrey_coeff_set::Ex_Ey_Bz, cdtodz, WarpX::l_lower_order_in_v) ); + // Same filter for fields Bx, By and Ez + nci_godfrey_filter_bxbyez[lev].reset( new NCIGodfreyFilter(godfrey_coeff_set::Bx_By_Ez, cdtodz, WarpX::l_lower_order_in_v) ); + // Compute Godfrey filters stencils + nci_godfrey_filter_exeybz[lev]->ComputeStencils(); + nci_godfrey_filter_bxbyez[lev]->ComputeStencils(); } } } diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 513f30b99..aaf7ead0e 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -181,17 +181,6 @@ public: void UpdateContinuousInjectionPosition(amrex::Real dt) const; int doContinuousInjection() const; - // - // Parameters for the Cherenkov corrector in the FDTD solver. - // Both stencils are calculated ar runtime. - // - // Number of coefficients for the stencil of the NCI corrector. - // The stencil is applied in the z direction only. - static constexpr int nstencilz_fdtd_nci_corr=5; - - amrex::Vector<amrex::Array<amrex::Real, nstencilz_fdtd_nci_corr> > fdtd_nci_stencilz_ex; - amrex::Vector<amrex::Array<amrex::Real, nstencilz_fdtd_nci_corr> > fdtd_nci_stencilz_by; - std::vector<std::string> GetSpeciesNames() const { return species_names; } PhysicalParticleContainer& GetPCtmp () { return *pc_tmp; } diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 6d618c096..9d39ec2f9 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -8,8 +8,6 @@ using namespace amrex; -constexpr int MultiParticleContainer::nstencilz_fdtd_nci_corr; - MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) { ReadParameters(); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 212084e64..ab37bf6ff 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1132,8 +1132,9 @@ PhysicalParticleContainer::Evolve (int lev, const std::array<Real,3>& dx = WarpX::CellSize(lev); const std::array<Real,3>& cdx = WarpX::CellSize(std::max(lev-1,0)); - const auto& mypc = WarpX::GetInstance().GetPartContainer(); - const int nstencilz_fdtd_nci_corr = mypc.nstencilz_fdtd_nci_corr; + // Get instances of NCI Godfrey filters + const auto& nci_godfrey_filter_exeybz = WarpX::GetInstance().nci_godfrey_filter_exeybz; + const auto& nci_godfrey_filter_bxbyez = WarpX::GetInstance().nci_godfrey_filter_bxbyez; BL_ASSERT(OnSameGrids(lev,jx)); @@ -1165,7 +1166,7 @@ PhysicalParticleContainer::Evolve (int lev, { Real wt = amrex::second(); - const Box& box = pti.validbox(); + const Box& box = pti.validbox(); auto& attribs = pti.GetAttribs(); @@ -1182,7 +1183,7 @@ PhysicalParticleContainer::Evolve (int lev, const long np = pti.numParticles(); - // Data on the grid + // Data on the grid FArrayBox const* exfab = &(Ex[pti]); FArrayBox const* eyfab = &(Ey[pti]); FArrayBox const* ezfab = &(Ez[pti]); @@ -1190,6 +1191,7 @@ PhysicalParticleContainer::Evolve (int lev, FArrayBox const* byfab = &(By[pti]); FArrayBox const* bzfab = &(Bz[pti]); + Elixir exeli, eyeli, ezeli, bxeli, byeli, bzeli; if (WarpX::use_fdtd_nci_corr) { #if (AMREX_SPACEDIM == 2) @@ -1201,54 +1203,43 @@ PhysicalParticleContainer::Evolve (int lev, static_cast<int>(WarpX::noz)}); #endif - // both 2d and 3d + // Filter Ex (Both 2D and 3D) filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ex), - BL_TO_FORTRAN_ANYD(filtered_Ex), - BL_TO_FORTRAN_ANYD(Ex[pti]), - mypc.fdtd_nci_stencilz_ex[lev].data(), - &nstencilz_fdtd_nci_corr); + // Safeguard for GPU + exeli = filtered_Ex.elixir(); + // Apply filter on Ex, result stored in filtered_Ex + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ex, Ex[pti], filtered_Ex.box()); + // Update exfab reference exfab = &filtered_Ex; + // Filter Ez filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ez), - BL_TO_FORTRAN_ANYD(filtered_Ez), - BL_TO_FORTRAN_ANYD(Ez[pti]), - mypc.fdtd_nci_stencilz_by[lev].data(), - &nstencilz_fdtd_nci_corr); + ezeli = filtered_Ez.elixir(); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Ez, Ez[pti], filtered_Ez.box()); ezfab = &filtered_Ez; + // Filter By filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_By), - BL_TO_FORTRAN_ANYD(filtered_By), - BL_TO_FORTRAN_ANYD(By[pti]), - mypc.fdtd_nci_stencilz_by[lev].data(), - &nstencilz_fdtd_nci_corr); + byeli = filtered_By.elixir(); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_By, By[pti], filtered_By.box()); byfab = &filtered_By; - #if (AMREX_SPACEDIM == 3) + // Filter Ey filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ey), - BL_TO_FORTRAN_ANYD(filtered_Ey), - BL_TO_FORTRAN_ANYD(Ey[pti]), - mypc.fdtd_nci_stencilz_ex[lev].data(), - &nstencilz_fdtd_nci_corr); + eyeli = filtered_Ey.elixir(); + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ey, Ey[pti], filtered_Ey.box()); eyfab = &filtered_Ey; + // Filter Bx filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bx), - BL_TO_FORTRAN_ANYD(filtered_Bx), - BL_TO_FORTRAN_ANYD(Bx[pti]), - mypc.fdtd_nci_stencilz_by[lev].data(), - &nstencilz_fdtd_nci_corr); + bxeli = filtered_Bx.elixir(); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Bx, Bx[pti], filtered_Bx.box()); bxfab = &filtered_Bx; + // Filter Bz filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bz), - BL_TO_FORTRAN_ANYD(filtered_Bz), - BL_TO_FORTRAN_ANYD(Bz[pti]), - mypc.fdtd_nci_stencilz_ex[lev].data(), - &nstencilz_fdtd_nci_corr); + bzeli = filtered_Bz.elixir(); + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Bz, Bz[pti], filtered_Bz.box()); bzfab = &filtered_Bz; #endif } @@ -1424,53 +1415,43 @@ PhysicalParticleContainer::Evolve (int lev, static_cast<int>(WarpX::noz)}); #endif - // both 2d and 3d + // Filter Ex (both 2D and 3D) filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ex), - BL_TO_FORTRAN_ANYD(filtered_Ex), - BL_TO_FORTRAN_ANYD((*cEx)[pti]), - mypc.fdtd_nci_stencilz_ex[lev-1].data(), - &nstencilz_fdtd_nci_corr); + // Safeguard for GPU + exeli = filtered_Ex.elixir(); + // Apply filter on Ex, result stored in filtered_Ex + nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Ex, (*cEx)[pti], filtered_Ex.box()); + // Update exfab reference cexfab = &filtered_Ex; + // Filter Ez filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ez), - BL_TO_FORTRAN_ANYD(filtered_Ez), - BL_TO_FORTRAN_ANYD((*cEz)[pti]), - mypc.fdtd_nci_stencilz_by[lev-1].data(), - &nstencilz_fdtd_nci_corr); + ezeli = filtered_Ez.elixir(); + nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_Ez, (*cEz)[pti], filtered_Ez.box()); cezfab = &filtered_Ez; + + // Filter By filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_By), - BL_TO_FORTRAN_ANYD(filtered_By), - BL_TO_FORTRAN_ANYD((*cBy)[pti]), - mypc.fdtd_nci_stencilz_by[lev-1].data(), - &nstencilz_fdtd_nci_corr); + byeli = filtered_By.elixir(); + nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_By, (*cBy)[pti], filtered_By.box()); cbyfab = &filtered_By; - #if (AMREX_SPACEDIM == 3) + // Filter Ey filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ey), - BL_TO_FORTRAN_ANYD(filtered_Ey), - BL_TO_FORTRAN_ANYD((*cEy)[pti]), - mypc.fdtd_nci_stencilz_ex[lev-1].data(), - &nstencilz_fdtd_nci_corr); + eyeli = filtered_Ey.elixir(); + nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Ey, (*cEy)[pti], filtered_Ey.box()); ceyfab = &filtered_Ey; + // Filter Bx filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bx), - BL_TO_FORTRAN_ANYD(filtered_Bx), - BL_TO_FORTRAN_ANYD((*cBx)[pti]), - mypc.fdtd_nci_stencilz_by[lev-1].data(), - &nstencilz_fdtd_nci_corr); + bxeli = filtered_Bx.elixir(); + nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_Bx, (*cBx)[pti], filtered_Bx.box()); cbxfab = &filtered_Bx; + // Filter Bz filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bz), - BL_TO_FORTRAN_ANYD(filtered_Bz), - BL_TO_FORTRAN_ANYD((*cBz)[pti]), - mypc.fdtd_nci_stencilz_ex[lev-1].data(), - &nstencilz_fdtd_nci_corr); + bzeli = filtered_Bz.elixir(); + nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Bz, (*cBz)[pti], filtered_Bz.box()); cbzfab = &filtered_Bz; #endif } diff --git a/Source/Utils/Make.package b/Source/Utils/Make.package index d8e2d2dab..cd335dbcf 100644 --- a/Source/Utils/Make.package +++ b/Source/Utils/Make.package @@ -4,6 +4,9 @@ CEXE_sources += WarpXTagging.cpp CEXE_sources += WarpXUtil.cpp CEXE_headers += WarpXConst.H CEXE_headers += WarpXUtil.H +CEXE_headers += WarpXAlgorithmSelection.H +CEXE_sources += WarpXAlgorithmSelection.cpp +CEXE_headers += NCIGodfreyTables.H INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Utils VPATH_LOCATIONS += $(WARPX_HOME)/Source/Utils diff --git a/Source/Utils/NCIGodfreyTables.H b/Source/Utils/NCIGodfreyTables.H new file mode 100644 index 000000000..8cb105aa0 --- /dev/null +++ b/Source/Utils/NCIGodfreyTables.H @@ -0,0 +1,224 @@ +#include <AMReX_AmrCore.H> + +#ifndef WARPX_GODFREY_COEFF_TABLE_H_ +#define WARPX_GODFREY_COEFF_TABLE_H_ + +// Table width. This is related to the stencil length +const int tab_width = 4; +// table length. Each line correspond to 1 value of cdtodz +// (here 101 values). +const int tab_length = 101; + +// Table of coefficient for Ex, Ey abd Bz +// We typically interpolate between two lines +const amrex::Real table_nci_godfrey_Ex_Ey_Bz[tab_length][tab_width]{ + -2.47536,2.04288,-0.598163,0.0314711, + -2.47536,2.04288,-0.598163,0.0314711, + -2.47545,2.04309,-0.598307,0.0315029, + -2.4756,2.04342,-0.598549,0.0315558, + -2.47581,2.0439,-0.598886,0.0316298, + -2.47608,2.0445,-0.59932,0.031725, + -2.47641,2.04525,-0.59985,0.0318412, + -2.4768,2.04612,-0.600477,0.0319785, + -2.47725,2.04714,-0.6012,0.0321367, + -2.47776,2.04829,-0.602019,0.0323158, + -2.47833,2.04957,-0.602934,0.0325158, + -2.47896,2.05099,-0.603944,0.0327364, + -2.47965,2.05254,-0.605051,0.0329777, + -2.4804,2.05423,-0.606253,0.0332396, + -2.48121,2.05606,-0.60755,0.0335218, + -2.48208,2.05802,-0.608942,0.0338243, + -2.48301,2.06012,-0.610429,0.0341469, + -2.48401,2.06235,-0.61201,0.0344895, + -2.48506,2.06471,-0.613685,0.0348519, + -2.48618,2.06721,-0.615453,0.0352339, + -2.48735,2.06984,-0.617314,0.0356353, + -2.48859,2.07261,-0.619268,0.0360559, + -2.48988,2.0755,-0.621312,0.0364954, + -2.49123,2.07853,-0.623447,0.0369536, + -2.49265,2.08169,-0.625672,0.0374302, + -2.49412,2.08498,-0.627986,0.0379248, + -2.49565,2.0884,-0.630386,0.0384372, + -2.49724,2.09194,-0.632873,0.0389669, + -2.49888,2.09561,-0.635443,0.0395135, + -2.50058,2.09939,-0.638096,0.0400766, + -2.50234,2.1033,-0.640829,0.0406557, + -2.50415,2.10732,-0.64364,0.0412502, + -2.50601,2.11145,-0.646526,0.0418594, + -2.50791,2.1157,-0.649485,0.0424828, + -2.50987,2.12004,-0.652512,0.0431196, + -2.51187,2.12448,-0.655604,0.0437688, + -2.51392,2.12901,-0.658756,0.0444297, + -2.516,2.13363,-0.661964,0.0451011, + -2.51812,2.13832,-0.665221,0.0457818, + -2.52027,2.14308,-0.668521,0.0464705, + -2.52244,2.14789,-0.671856,0.0471658, + -2.52464,2.15274,-0.675218,0.0478658, + -2.52684,2.15762,-0.678596,0.0485687, + -2.52906,2.16251,-0.68198,0.0492723, + -2.53126,2.16738,-0.685355,0.049974, + -2.53345,2.17222,-0.688706,0.0506708, + -2.53561,2.177,-0.692015,0.0513594, + -2.53773,2.18168,-0.69526,0.0520359, + -2.53978,2.18623,-0.698416,0.0526955, + -2.54175,2.19059,-0.701452,0.053333, + -2.5436,2.19471,-0.704331,0.0539417, + -2.54531,2.19852,-0.70701,0.0545141, + -2.54683,2.20193,-0.709433,0.0550409, + -2.5481,2.20483,-0.711533,0.0555106, + -2.54906,2.20709,-0.713224,0.0559094, + -2.54963,2.20852,-0.714397,0.0562198, + -2.54968,2.20888,-0.714907,0.0564196, + -2.54905,2.20785,-0.714562,0.0564797, + -2.54751,2.20496,-0.713094,0.0563618, + -2.54472,2.19955,-0.710118,0.0560124, + -2.54014,2.19058,-0.705048,0.0553544, + -2.53286,2.1763,-0.69693,0.0542684, + -2.52115,2.15344,-0.684027,0.05255, + -2.50098,2.11466,-0.66255,0.0497817, + -2.45797,2.03459,-0.620099,0.0446889, + -2.28371,1.72254,-0.465905,0.0283268, + -2.4885,2.04899,-0.599292,0.0390466, + -2.1433,1.36735,-0.220924,-0.00215633, + -2.4943,2.07019,-0.610552,0.035166, + -2.84529,2.77303,-1.00018,0.0724884, + -2.72242,2.51888,-0.847226,0.0509964, + -2.65633,2.3744,-0.750392,0.0326366, + -2.59601,2.23412,-0.646421,0.00868027, + -2.51477,2.0369,-0.491066,-0.0306397, + -2.35935,1.65155,-0.178971,-0.112713, + -1.84315,0.361693,0.876104,-0.393844, + -2.65422,2.39262,-0.789663,0.0516265, + -3.46529,4.42354,-2.45543,0.497097, + -3.15747,3.65311,-1.824,0.328432, + -3.04694,3.37613,-1.59668,0.267631, + -2.99205,3.23814,-1.48302,0.237103, + -2.96075,3.15894,-1.41733,0.219317, + -2.94172,3.11028,-1.37649,0.20811, + -2.92994,3.07962,-1.35025,0.200755, + -2.92283,3.06054,-1.33338,0.195859, + -2.91894,3.04938,-1.3229,0.192637, + -2.91736,3.04394,-1.31702,0.190612, + -2.91753,3.04278,-1.31456,0.189477, + -2.91905,3.04494,-1.31475,0.189026, + -2.92165,3.04973,-1.31705,0.189117, + -2.92512,3.05667,-1.32105,0.189646, + -2.92933,3.06539,-1.32646,0.190538, + -2.93416,3.07562,-1.33308,0.191735, + -2.93952,3.08715,-1.34072,0.193194, + -2.94535,3.09982,-1.34925,0.194881, + -2.95159,3.11349,-1.35858,0.196769, + -2.9582,3.12805,-1.36861,0.198838, + -2.96514,3.14342,-1.37929,0.201068, + -2.97239,3.15953,-1.39055,0.203448, + -2.97991,3.17632,-1.40234,0.205964, + -2.98769,3.19374,-1.41463,0.208607 + }; + +// Table of coefficient for Bx, By and Ez +// We typically interpolate between two lines +const amrex::Real table_nci_godfrey_Bx_By_Ez[tab_length][tab_width]{ + -2.80862,2.80104,-1.14615,0.154077, + -2.80862,2.80104,-1.14615,0.154077, + -2.80851,2.80078,-1.14595,0.154027, + -2.80832,2.80034,-1.14561,0.153945, + -2.80807,2.79973,-1.14514,0.153829, + -2.80774,2.79894,-1.14454,0.15368, + -2.80733,2.79798,-1.1438,0.153498, + -2.80685,2.79685,-1.14292,0.153284, + -2.8063,2.79554,-1.14192,0.153036, + -2.80568,2.79405,-1.14077,0.152756, + -2.80498,2.79239,-1.1395,0.152443, + -2.80421,2.79056,-1.13809,0.152098, + -2.80337,2.78856,-1.13656,0.151721, + -2.80246,2.78638,-1.13488,0.151312, + -2.80147,2.78404,-1.13308,0.150871, + -2.80041,2.78152,-1.13115,0.150397, + -2.79927,2.77882,-1.12908,0.149893, + -2.79807,2.77596,-1.12689,0.149356, + -2.79679,2.77292,-1.12456,0.148789, + -2.79543,2.76972,-1.12211,0.14819, + -2.79401,2.76634,-1.11953,0.14756, + -2.79251,2.76279,-1.11681,0.1469, + -2.79094,2.75907,-1.11397,0.146208, + -2.78929,2.75517,-1.111,0.145486, + -2.78757,2.7511,-1.10789,0.144733, + -2.78578,2.74686,-1.10466,0.14395, + -2.78391,2.74245,-1.1013,0.143137, + -2.78196,2.73786,-1.09781,0.142293, + -2.77994,2.73309,-1.09419,0.141419, + -2.77784,2.72814,-1.09043,0.140514, + -2.77566,2.72301,-1.08654,0.139578, + -2.7734,2.7177,-1.08252,0.138612, + -2.77106,2.7122,-1.07836,0.137614, + -2.76864,2.70651,-1.07406,0.136586, + -2.76613,2.70062,-1.06962,0.135525, + -2.76353,2.69453,-1.06503,0.134432, + -2.76084,2.68824,-1.0603,0.133307, + -2.75806,2.68173,-1.05541,0.132148, + -2.75518,2.675,-1.05037,0.130954, + -2.75219,2.66804,-1.04516,0.129725, + -2.7491,2.66084,-1.03978,0.12846, + -2.7459,2.65339,-1.03423,0.127156, + -2.74257,2.64566,-1.02848,0.125813, + -2.73912,2.63765,-1.02254,0.124428, + -2.73552,2.62934,-1.01638,0.122999, + -2.73178,2.62069,-1.01,0.121523, + -2.72787,2.61169,-1.00337,0.119996, + -2.72379,2.6023,-0.996479,0.118417, + -2.71951,2.59248,-0.989294,0.116778, + -2.71501,2.58218,-0.981786,0.115076, + -2.71026,2.57135,-0.97392,0.113303, + -2.70524,2.55991,-0.965651,0.111453, + -2.69989,2.54778,-0.956922,0.109514, + -2.69416,2.53484,-0.947666,0.107476, + -2.68799,2.52096,-0.937795,0.105324, + -2.68129,2.50596,-0.927197,0.103039, + -2.67394,2.48959,-0.915724,0.100597, + -2.66578,2.47153,-0.903179,0.097968, + -2.65657,2.4513,-0.889283,0.0951084, + -2.64598,2.42824,-0.873638,0.0919592, + -2.63347,2.40127,-0.855632,0.0884325, + -2.61813,2.36864,-0.834261,0.0843898, + -2.59821,2.32701,-0.807691,0.0795876, + -2.56971,2.26887,-0.77188,0.0735132, + -2.51823,2.16823,-0.713448,0.0645399, + -2.33537,1.8294,-0.533852,0.0409941, + -2.53143,2.14818,-0.670502,0.053982, + -2.17737,1.43641,-0.259095,0.00101255, + -2.51929,2.12931,-0.654743,0.0452381, + -2.86122,2.82221,-1.05039,0.0894636, + -2.72908,2.54506,-0.87834,0.0626188, + -2.6536,2.37954,-0.7665,0.0409117, + -2.58374,2.21923,-0.649738,0.0146791, + -2.49284,2.00346,-0.48457,-0.0255348, + -2.32762,1.60337,-0.1698,-0.105287, + -1.80149,0.316787,0.855414,-0.369652, + -2.60242,2.28418,-0.721378,0.040091, + -3.40335,4.25157,-2.29817,0.449834, + -3.0852,3.47341,-1.67791,0.28982, + -2.9642,3.17856,-1.44399,0.229852, + -2.89872,3.01966,-1.31861,0.197945, + -2.85668,2.91811,-1.23894,0.17783, + -2.82679,2.84621,-1.18287,0.163785, + -2.80401,2.79167,-1.14058,0.153278, + -2.78577,2.74819,-1.10706,0.145015, + -2.77061,2.7122,-1.07946,0.138267, + -2.75764,2.68152,-1.05606,0.132589, + -2.74627,2.65475,-1.03575,0.127695, + -2.73612,2.63093,-1.01777,0.123395, + -2.72692,2.6094,-1.00159,0.119553, + -2.71846,2.58968,-0.986841,0.116074, + -2.71061,2.57142,-0.973239,0.112887, + -2.70323,2.55434,-0.960573,0.109937, + -2.69626,2.53824,-0.948678,0.107185, + -2.68962,2.52294,-0.937429,0.104598, + -2.68327,2.50833,-0.926722,0.102151, + -2.67714,2.4943,-0.916477,0.0998223, + -2.67122,2.48076,-0.906627,0.0975966, + -2.66546,2.46764,-0.897118,0.0954599, + -2.65985,2.45489,-0.887903,0.0934011, + -2.65437,2.44244,-0.878945,0.0914107 + }; + +#endif // #ifndef WARPX_GODFREY_COEFF_TABLE_H_ diff --git a/Source/Utils/WarpXAlgorithmSelection.H b/Source/Utils/WarpXAlgorithmSelection.H new file mode 100644 index 000000000..3fb23698a --- /dev/null +++ b/Source/Utils/WarpXAlgorithmSelection.H @@ -0,0 +1,57 @@ +#ifndef UTILS_WARPXALGORITHMSELECTION_H_ +#define UTILS_WARPXALGORITHMSELECTION_H_ + +#include <AMReX_ParmParse.H> +#include <string> + +struct MaxwellSolverAlgo { + // These numbers corresponds to the algorithm code in WarpX's + // `warpx_push_bvec` and `warpx_push_evec_f` + enum { + Yee = 0, + CKC = 1 + }; +}; + +struct ParticlePusherAlgo { + // These numbers correspond to the algorithm code in WarpX's + // `warpx_particle_pusher` + enum { + Boris = 0, + Vay = 1 + }; +}; + +struct CurrentDepositionAlgo { + enum { + // These numbers corresponds to the algorithm code in PICSAR's + // `depose_jxjyjz_generic` and `depose_jxjyjz_generic_2d` + Direct = 3, + DirectVectorized = 2, + EsirkepovNonOptimized = 1, + Esirkepov = 0 + }; +}; + +struct ChargeDepositionAlgo { + // These numbers corresponds to the algorithm code in WarpX's + // `warpx_charge_deposition` function + enum { + Vectorized = 0, + Standard = 1 + }; +}; + +struct GatheringAlgo { + // These numbers corresponds to the algorithm code in PICSAR's + // `geteb3d_energy_conserving_generic` function + enum { + Vectorized = 0, + Standard = 1 + }; +}; + +int +GetAlgorithmInteger( amrex::ParmParse& pp, const char* pp_search_key ); + +#endif // UTILS_WARPXALGORITHMSELECTION_H_ diff --git a/Source/Utils/WarpXAlgorithmSelection.cpp b/Source/Utils/WarpXAlgorithmSelection.cpp new file mode 100644 index 000000000..89802e064 --- /dev/null +++ b/Source/Utils/WarpXAlgorithmSelection.cpp @@ -0,0 +1,94 @@ +#include <WarpXAlgorithmSelection.H> +#include <map> +#include <algorithm> + +// Define dictionary with correspondance between user-input strings, +// and corresponding integer for use inside the code (e.g. in PICSAR). + +const std::map<std::string, int> maxwell_solver_algo_to_int = { + {"yee", MaxwellSolverAlgo::Yee }, +#ifndef WARPX_RZ // Not available in RZ + {"ckc", MaxwellSolverAlgo::CKC }, +#endif + {"default", MaxwellSolverAlgo::Yee } +}; + +const std::map<std::string, int> particle_pusher_algo_to_int = { + {"boris", ParticlePusherAlgo::Boris }, + {"vay", ParticlePusherAlgo::Vay }, + {"default", ParticlePusherAlgo::Boris } +}; + +const std::map<std::string, int> current_deposition_algo_to_int = { + {"esirkepov", CurrentDepositionAlgo::Esirkepov }, + {"direct", CurrentDepositionAlgo::Direct }, +#if (!defined AMREX_USE_GPU)&&(AMREX_SPACEDIM == 3) // Only available on CPU and 3D + {"direct-vectorized", CurrentDepositionAlgo::DirectVectorized }, +#endif + {"default", CurrentDepositionAlgo::Esirkepov } +}; + +const std::map<std::string, int> charge_deposition_algo_to_int = { + {"standard", ChargeDepositionAlgo::Standard }, +#if (!defined AMREX_USE_GPU)&&(AMREX_SPACEDIM == 3) // Only available on CPU and 3D + {"vectorized", ChargeDepositionAlgo::Vectorized }, + {"default", ChargeDepositionAlgo::Vectorized } +#else + {"default", ChargeDepositionAlgo::Standard } +#endif +}; + +const std::map<std::string, int> gathering_algo_to_int = { + {"standard", GatheringAlgo::Standard }, +#ifndef AMREX_USE_GPU // Only available on CPU + {"vectorized", GatheringAlgo::Vectorized }, + {"default", GatheringAlgo::Vectorized } +#else + {"default", GatheringAlgo::Standard } +#endif +}; + + +int +GetAlgorithmInteger( amrex::ParmParse& pp, const char* pp_search_key ){ + + // Read user input ; use "default" if it is not found + std::string algo = "default"; + pp.query( pp_search_key, algo ); + // Convert to lower case + std::transform(algo.begin(), algo.end(), algo.begin(), ::tolower); + + // Pick the right dictionary + std::map<std::string, int> algo_to_int; + if (pp_search_key == "maxwell_fdtd_solver") { + algo_to_int = maxwell_solver_algo_to_int; + } else if (pp_search_key == "particle_pusher") { + algo_to_int = particle_pusher_algo_to_int; + } else if (pp_search_key == "current_deposition") { + algo_to_int = current_deposition_algo_to_int; + } else if (pp_search_key == "charge_deposition") { + algo_to_int = charge_deposition_algo_to_int; + } else if (pp_search_key == "field_gathering") { + algo_to_int = gathering_algo_to_int; + } else { + std::string pp_search_string = pp_search_key; + amrex::Abort("Unknown algorithm type: " + pp_search_string); + } + + // Check if the user-input is a valid key for the dictionary + if (algo_to_int.count(algo) == 0){ + // Not a valid key ; print error message + std::string pp_search_string = pp_search_key; + std::string error_message = "Invalid string for algo." + pp_search_string + + ": " + algo + ".\nThe valid values are:\n"; + for ( const auto &valid_pair : algo_to_int ) { + if (valid_pair.first != "default"){ + error_message += " - " + valid_pair.first + "\n"; + } + } + amrex::Abort(error_message); + } + + // If the input is a valid key, return the value + return algo_to_int[algo]; +} diff --git a/Source/Utils/WarpXMovingWindow.cpp b/Source/Utils/WarpXMovingWindow.cpp index a0ab1f26f..ae781f9aa 100644 --- a/Source/Utils/WarpXMovingWindow.cpp +++ b/Source/Utils/WarpXMovingWindow.cpp @@ -65,6 +65,26 @@ WarpX::MoveWindow (bool move_j) ResetProbDomain(RealBox(new_lo, new_hi)); + // Moving slice coordinates - lo and hi - with moving window // + // slice box is modified only if slice diagnostics is initialized in input // + if ( slice_plot_int > 0 ) + { + Real new_slice_lo[AMREX_SPACEDIM]; + Real new_slice_hi[AMREX_SPACEDIM]; + const Real* current_slice_lo = slice_realbox.lo(); + const Real* current_slice_hi = slice_realbox.hi(); + for ( int i = 0; i < AMREX_SPACEDIM; i++) { + new_slice_lo[i] = current_slice_lo[i]; + new_slice_hi[i] = current_slice_hi[i]; + } + int num_shift_base_slice = static_cast<int> ((moving_window_x - + current_slice_lo[dir]) / cdx[dir]); + new_slice_lo[dir] = current_slice_lo[dir] + num_shift_base_slice*cdx[dir]; + new_slice_hi[dir] = current_slice_hi[dir] + num_shift_base_slice*cdx[dir]; + slice_realbox.setLo(new_slice_lo); + slice_realbox.setHi(new_slice_hi); + } + int num_shift = num_shift_base; int num_shift_crse = num_shift; diff --git a/Source/WarpX.H b/Source/WarpX.H index 35b072142..2580c5320 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -23,6 +23,7 @@ #include <PML.H> #include <BoostedFrameDiagnostic.H> #include <BilinearFilter.H> +#include <NCIGodfreyFilter.H> #ifdef WARPX_USE_PSATD #include <fftw3.h> @@ -146,7 +147,9 @@ public: static amrex::IntVect filter_npass_each_dir; BilinearFilter bilinear_filter; - + amrex::Vector< std::unique_ptr<NCIGodfreyFilter> > nci_godfrey_filter_exeybz; + amrex::Vector< std::unique_ptr<NCIGodfreyFilter> > nci_godfrey_filter_bxbyez; + static int num_mirrors; amrex::Vector<amrex::Real> mirror_z; amrex::Vector<amrex::Real> mirror_z_width; @@ -242,6 +245,12 @@ public: static int do_moving_window; static int moving_window_dir; + // slice generation // + void InitializeSliceMultiFabs (); + void SliceGenerationForDiagnostics (); + void WriteSlicePlotFile () const; + void ClearSliceMultiFabs (); + protected: //! Tagging cells for refinement @@ -553,6 +562,17 @@ private: bool is_synchronized = true; + //Slice Parameters + int slice_max_grid_size; + int slice_plot_int = -1; + amrex::RealBox slice_realbox; + amrex::IntVect slice_cr_ratio; + amrex::Vector< std::unique_ptr<amrex::MultiFab> > F_slice; + amrex::Vector< std::unique_ptr<amrex::MultiFab> > rho_slice; + amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_slice; + amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Efield_slice; + amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Bfield_slice; + #ifdef WARPX_USE_PSATD // Store fields in real space on the dual grid (i.e. the grid for the FFT push of the fields) // This includes data for the FFT guard cells (between FFT groups) diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 3d7f7dcc5..c2cf97f30 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -17,6 +17,7 @@ #include <WarpXConst.H> #include <WarpXWrappers.h> #include <WarpXUtil.H> +#include <WarpXAlgorithmSelection.H> #ifdef BL_USE_SENSEI_INSITU #include <AMReX_AmrMeshInSituBridge.H> @@ -35,11 +36,11 @@ Vector<int> WarpX::boost_direction = {0,0,0}; int WarpX::do_compute_max_step_from_zmax = 0; Real WarpX::zmax_plasma_to_compute_max_step = 0.; -long WarpX::current_deposition_algo = 3; -long WarpX::charge_deposition_algo = 0; -long WarpX::field_gathering_algo = 1; -long WarpX::particle_pusher_algo = 0; -int WarpX::maxwell_fdtd_solver_id = 0; +long WarpX::current_deposition_algo; +long WarpX::charge_deposition_algo; +long WarpX::field_gathering_algo; +long WarpX::particle_pusher_algo; +int WarpX::maxwell_fdtd_solver_id; long WarpX::nox = 1; long WarpX::noy = 1; @@ -118,7 +119,7 @@ WarpX::ResetInstance () { delete m_instance; m_instance = nullptr; -} +} WarpX::WarpX () { @@ -228,6 +229,11 @@ WarpX::WarpX () #ifdef BL_USE_SENSEI_INSITU insitu_bridge = nullptr; #endif + + // NCI Godfrey filters can have different stencils + // at different levels (the stencil depends on c*dt/dz) + nci_godfrey_filter_exeybz.resize(nlevs_max); + nci_godfrey_filter_bxbyez.resize(nlevs_max); } WarpX::~WarpX () @@ -276,10 +282,10 @@ WarpX::ReadParameters () ReadBoostedFrameParameters(gamma_boost, beta_boost, boost_direction); - // pp.query returns 1 if argument zmax_plasma_to_compute_max_step is + // pp.query returns 1 if argument zmax_plasma_to_compute_max_step is // specified by the user, 0 otherwise. - do_compute_max_step_from_zmax = - pp.query("zmax_plasma_to_compute_max_step", + do_compute_max_step_from_zmax = + pp.query("zmax_plasma_to_compute_max_step", zmax_plasma_to_compute_max_step); pp.queryarr("B_external", B_external); @@ -318,7 +324,7 @@ WarpX::ReadParameters () "gamma_boost must be > 1 to use the boosted frame diagnostic."); pp.query("lab_data_directory", lab_data_directory); - + std::string s; pp.get("boost_direction", s); AMREX_ALWAYS_ASSERT_WITH_MESSAGE( (s == "z" || s == "Z"), @@ -476,29 +482,12 @@ WarpX::ReadParameters () } { - ParmParse pp("algo"); - pp.query("current_deposition", current_deposition_algo); - pp.query("charge_deposition", charge_deposition_algo); - pp.query("field_gathering", field_gathering_algo); - pp.query("particle_pusher", particle_pusher_algo); - std::string s_solver = ""; - pp.query("maxwell_fdtd_solver", s_solver); - std::transform(s_solver.begin(), - s_solver.end(), - s_solver.begin(), - ::tolower); - // if maxwell_fdtd_solver is specified, set the value - // of maxwell_fdtd_solver_id accordingly. - // Otherwise keep the default value maxwell_fdtd_solver_id=0 - if (s_solver != "") { - if (s_solver == "yee") { - maxwell_fdtd_solver_id = 0; - } else if (s_solver == "ckc") { - maxwell_fdtd_solver_id = 1; - } else { - amrex::Abort("Unknown FDTD Solver type " + s_solver); - } - } + ParmParse pp("algo"); + current_deposition_algo = GetAlgorithmInteger(pp, "current_deposition"); + charge_deposition_algo = GetAlgorithmInteger(pp, "charge_deposition"); + field_gathering_algo = GetAlgorithmInteger(pp, "field_gathering"); + particle_pusher_algo = GetAlgorithmInteger(pp, "particle_pusher"); + maxwell_fdtd_solver_id = GetAlgorithmInteger(pp, "maxwell_fdtd_solver"); } #ifdef WARPX_USE_PSATD @@ -527,6 +516,34 @@ WarpX::ReadParameters () pp.query("config", insitu_config); pp.query("pin_mesh", insitu_pin_mesh); } + + // for slice generation // + { + ParmParse pp("slice"); + amrex::Vector<Real> slice_lo(AMREX_SPACEDIM); + amrex::Vector<Real> slice_hi(AMREX_SPACEDIM); + Vector<int> slice_crse_ratio(AMREX_SPACEDIM); + // set default slice_crse_ratio // + for (int idim=0; idim < AMREX_SPACEDIM; ++idim ) + { + slice_crse_ratio[idim] = 1; + } + pp.queryarr("dom_lo",slice_lo,0,AMREX_SPACEDIM); + pp.queryarr("dom_hi",slice_hi,0,AMREX_SPACEDIM); + pp.queryarr("coarsening_ratio",slice_crse_ratio,0,AMREX_SPACEDIM); + pp.query("plot_int",slice_plot_int); + slice_realbox.setLo(slice_lo); + slice_realbox.setHi(slice_hi); + slice_cr_ratio = IntVect(AMREX_D_DECL(1,1,1)); + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) + { + if (slice_crse_ratio[idim] > 1 ) { + slice_cr_ratio[idim] = slice_crse_ratio[idim]; + } + } + + } + } // This is a virtual function. @@ -628,7 +645,7 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d int ngz_nonci = (ngz_tmp % 2) ? ngz_tmp+1 : ngz_tmp; // Always even number int ngz; if (WarpX::use_fdtd_nci_corr) { - int ng = ngz_tmp + (mypc->nstencilz_fdtd_nci_corr-1); + int ng = ngz_tmp + NCIGodfreyFilter::stencil_width; ngz = (ng % 2) ? ng+1 : ng; } else { ngz = ngz_nonci; @@ -833,8 +850,6 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm if (load_balance_int > 0) { costs[lev].reset(new MultiFab(ba, dm, 1, 0)); } - - } std::array<Real,3> |