diff options
Diffstat (limited to 'Source')
| -rw-r--r-- | Source/Diagnostics/FieldIO.H | 55 | ||||
| -rw-r--r-- | Source/Diagnostics/FieldIO.cpp | 288 | ||||
| -rw-r--r-- | Source/Diagnostics/WarpXIO.cpp | 303 |
3 files changed, 398 insertions, 248 deletions
diff --git a/Source/Diagnostics/FieldIO.H b/Source/Diagnostics/FieldIO.H index 1d94c07ce..121fb9374 100644 --- a/Source/Diagnostics/FieldIO.H +++ b/Source/Diagnostics/FieldIO.H @@ -20,6 +20,61 @@ AverageAndPackScalarField( MultiFab& mf_avg, const int dcomp, const int ngrow ); void +WriteRawField( const MultiFab& F, const DistributionMapping& dm, + const std::string& filename, + const std::string& level_prefix, + const std::string& field_name, + const int lev, const bool plot_guards ); + +void +WriteZeroRawField( const MultiFab& F, const DistributionMapping& dm, + const std::string& filename, + const std::string& level_prefix, + const std::string& field_name, + const int lev, const int ng ); + +void +WriteCoarseScalar( const std::string field_name, + const std::unique_ptr<MultiFab>& F_cp, + const std::unique_ptr<MultiFab>& F_fp, + const DistributionMapping& dm, + const std::string& filename, + const std::string& level_prefix, + const int lev, const bool plot_guards, + const int r_ratio, const Real* dx, const int icomp=0 ); + +void +WriteCoarseVector( const std::string field_name, + const std::unique_ptr<MultiFab>& Fx_cp, + const std::unique_ptr<MultiFab>& Fy_cp, + const std::unique_ptr<MultiFab>& Fz_cp, + const std::unique_ptr<MultiFab>& Fx_fp, + const std::unique_ptr<MultiFab>& Fy_fp, + const std::unique_ptr<MultiFab>& Fz_fp, + const DistributionMapping& dm, + const std::string& filename, + const std::string& level_prefix, + const int lev, const bool plot_guards, + const int r_ratio, const Real* dx ); + +std::unique_ptr<MultiFab> +getInterpolatedScalar( + const MultiFab& F_cp, const MultiFab& F_fp, + const DistributionMapping& dm, const int r_ratio, + const Real* dx, const int ngrow ); + +std::array<std::unique_ptr<MultiFab>, 3> +getInterpolatedVector( + const std::unique_ptr<MultiFab>& Fx_cp, + const std::unique_ptr<MultiFab>& Fy_cp, + const std::unique_ptr<MultiFab>& Fz_cp, + const std::unique_ptr<MultiFab>& Fx_fp, + const std::unique_ptr<MultiFab>& Fy_fp, + const std::unique_ptr<MultiFab>& Fz_fp, + const DistributionMapping& dm, + const int r_ratio, const Real* dx, + const int ngrow ); + coarsenCellCenteredFields( Vector<MultiFab>& coarse_mf, Vector<Geometry>& coarse_geom, const Vector<MultiFab>& source_mf, const Vector<Geometry>& source_geom, diff --git a/Source/Diagnostics/FieldIO.cpp b/Source/Diagnostics/FieldIO.cpp index f8fb86fb2..965324d02 100644 --- a/Source/Diagnostics/FieldIO.cpp +++ b/Source/Diagnostics/FieldIO.cpp @@ -2,6 +2,8 @@ #include <WarpX.H> #include <FieldIO.H> +#include <AMReX_FillPatchUtil_F.H> + using namespace amrex; void @@ -327,3 +329,289 @@ coarsenCellCenteredFields( average_down(source_mf[lev], coarse_mf[lev], 0, ncomp, IntVect(coarse_ratio)); } }; + + +/** \brief TODO + */ +void +WriteRawField( const MultiFab& F, const DistributionMapping& dm, + const std::string& filename, + const std::string& level_prefix, + const std::string& field_name, + const int lev, const bool plot_guards ) +{ + std::string prefix = amrex::MultiFabFileFullPrefix(lev, + filename, level_prefix, field_name); + + if (plot_guards) { + // Dump original MultiFab F + 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); + VisMF::Write(tmpF, prefix); + } + +} + +/** \brief TODO + */ +void +WriteZeroRawField( const MultiFab& F, const DistributionMapping& dm, + const std::string& filename, + const std::string& level_prefix, + const std::string& field_name, + const int lev, const int ng ) +{ + std::string prefix = amrex::MultiFabFileFullPrefix(lev, + filename, level_prefix, field_name); + + MultiFab tmpF(F.boxArray(), dm, 1, ng); + tmpF.setVal(0.); + VisMF::Write(tmpF, prefix); +} + +/** \brief TODO Writes same size as fp for compatibility with yt + */ +void +WriteCoarseScalar( const std::string field_name, + const std::unique_ptr<MultiFab>& F_cp, + const std::unique_ptr<MultiFab>& F_fp, + const DistributionMapping& dm, + const std::string& filename, + const std::string& level_prefix, + const int lev, const bool plot_guards, + const int r_ratio, const Real* dx, const int icomp ) +{ + int ng = 0; + if (plot_guards) ng = F_fp->nGrow(); + + if (lev == 0) { + // No coarse field for level 0: instead write a MultiFab + // filled with 0, with the same number of cells as the _fp field + WriteZeroRawField( *F_fp, dm, filename, level_prefix, field_name+"_cp", lev, ng ); + } else { + // Create an alias to the component `icomp` of F_cp + MultiFab F_comp(*F_cp, amrex::make_alias, 1, icomp); + auto F = getInterpolatedScalar( F_comp, *F_fp, dm, r_ratio, dx, ng ); + WriteRawField( *F, dm, filename, level_prefix, field_name+"_cp", lev, plot_guards ); + } +} + +/** \brief TODO Writes same size as fp for compatibility with yt + */ +void +WriteCoarseVector( const std::string field_name, + const std::unique_ptr<MultiFab>& Fx_cp, + const std::unique_ptr<MultiFab>& Fy_cp, + const std::unique_ptr<MultiFab>& Fz_cp, + const std::unique_ptr<MultiFab>& Fx_fp, + const std::unique_ptr<MultiFab>& Fy_fp, + const std::unique_ptr<MultiFab>& Fz_fp, + const DistributionMapping& dm, + const std::string& filename, + const std::string& level_prefix, + const int lev, const bool plot_guards, + const int r_ratio, const Real* dx ) +{ + int ng = 0; + if (plot_guards) ng = Fx_fp->nGrow(); + + if (lev == 0) { + // No coarse field for level 0: instead write a MultiFab + // filled with 0, with the same number of cells as the _fp field + WriteZeroRawField( *Fx_fp, dm, filename, level_prefix, field_name+"x_cp", lev, ng ); + WriteZeroRawField( *Fy_fp, dm, filename, level_prefix, field_name+"y_cp", lev, ng ); + WriteZeroRawField( *Fz_fp, dm, filename, level_prefix, field_name+"z_cp", lev, ng ); + } else { + auto F = getInterpolatedVector( Fx_cp, Fy_cp, Fz_cp, Fx_fp, Fy_fp, Fz_fp, + dm, r_ratio, dx, ng ); + WriteRawField( *F[0], dm, filename, level_prefix, field_name+"x_cp", lev, plot_guards ); + WriteRawField( *F[1], dm, filename, level_prefix, field_name+"y_cp", lev, plot_guards ); + WriteRawField( *F[2], dm, filename, level_prefix, field_name+"z_cp", lev, plot_guards ); + } +} + +/** \brief TODO + * + */ +std::unique_ptr<MultiFab> +getInterpolatedScalar( + const MultiFab& F_cp, const MultiFab& F_fp, + const DistributionMapping& dm, const int r_ratio, + const Real* dx, const int ngrow ) +{ + // Prepare the structure that will contain the returned fields + std::unique_ptr<MultiFab> interpolated_F; + interpolated_F.reset( new MultiFab(F_fp.boxArray(), dm, 1, ngrow) ); + interpolated_F->setVal(0.); + + // Loop through the boxes and interpolate the values from the _cp data + const int use_limiter = 0; +#ifdef _OPEMP +#pragma omp parallel +#endif + { + FArrayBox ffab; // Temporary array ; contains interpolated fields + for (MFIter mfi(*interpolated_F); mfi.isValid(); ++mfi) + { + Box ccbx = mfi.fabbox(); + ccbx.enclosedCells(); + ccbx.coarsen(r_ratio).refine(r_ratio); // so that ccbx is coarsenable + + const FArrayBox& cfab = (F_cp)[mfi]; + ffab.resize(amrex::convert(ccbx,(F_fp)[mfi].box().type())); + + // - Fully centered + if ( F_fp.is_cell_centered() ){ + amrex_interp_cc_bfield(ccbx.loVect(), ccbx.hiVect(), + BL_TO_FORTRAN_ANYD(ffab), + BL_TO_FORTRAN_ANYD(cfab), + &r_ratio, &use_limiter); + // - Edge centered, in the same way as E on a Yee grid + } else if ( F_fp.is_nodal() ){ + amrex_interp_nd_efield(ccbx.loVect(), ccbx.hiVect(), + BL_TO_FORTRAN_ANYD(ffab), + BL_TO_FORTRAN_ANYD(cfab), + &r_ratio); + } else { + amrex::Abort("Unknown field staggering."); + } + + // Add temporary array to the returned structure + const Box& bx = (*interpolated_F)[mfi].box(); + (*interpolated_F)[mfi].plus(ffab, bx, bx, 0, 0, 1); + } + } + return interpolated_F; +} + +/** \brief TODO + * + */ +std::array<std::unique_ptr<MultiFab>, 3> +getInterpolatedVector( + const std::unique_ptr<MultiFab>& Fx_cp, + const std::unique_ptr<MultiFab>& Fy_cp, + const std::unique_ptr<MultiFab>& Fz_cp, + const std::unique_ptr<MultiFab>& Fx_fp, + const std::unique_ptr<MultiFab>& Fy_fp, + const std::unique_ptr<MultiFab>& Fz_fp, + const DistributionMapping& dm, const int r_ratio, + const Real* dx, const int ngrow ) +{ + + // Prepare the structure that will contain the returned fields + std::array<std::unique_ptr<MultiFab>, 3> interpolated_F; + interpolated_F[0].reset( new MultiFab(Fx_fp->boxArray(), dm, 1, ngrow) ); + interpolated_F[1].reset( new MultiFab(Fy_fp->boxArray(), dm, 1, ngrow) ); + interpolated_F[2].reset( new MultiFab(Fz_fp->boxArray(), dm, 1, ngrow) ); + for (int i=0; i<3; i++) interpolated_F[i]->setVal(0.); + + // Loop through the boxes and interpolate the values from the _cp data + const int use_limiter = 0; +#ifdef _OPEMP +#pragma omp parallel +#endif + { + std::array<FArrayBox,3> ffab; // Temporary array ; contains interpolated fields + for (MFIter mfi(*interpolated_F[0]); mfi.isValid(); ++mfi) + { + Box ccbx = mfi.fabbox(); + ccbx.enclosedCells(); + ccbx.coarsen(r_ratio).refine(r_ratio); // so that ccbx is coarsenable + + const FArrayBox& cxfab = (*Fx_cp)[mfi]; + const FArrayBox& cyfab = (*Fy_cp)[mfi]; + const FArrayBox& czfab = (*Fz_cp)[mfi]; + ffab[0].resize(amrex::convert(ccbx,(*Fx_fp)[mfi].box().type())); + ffab[1].resize(amrex::convert(ccbx,(*Fy_fp)[mfi].box().type())); + ffab[2].resize(amrex::convert(ccbx,(*Fz_fp)[mfi].box().type())); + + // - Face centered, in the same way as B on a Yee grid + if ( (*Fx_fp)[mfi].box().type() == IntVect{AMREX_D_DECL(1,0,0)} ){ +#if (AMREX_SPACEDIM == 3) + amrex_interp_div_free_bfield(ccbx.loVect(), ccbx.hiVect(), + BL_TO_FORTRAN_ANYD(ffab[0]), + BL_TO_FORTRAN_ANYD(ffab[1]), + BL_TO_FORTRAN_ANYD(ffab[2]), + BL_TO_FORTRAN_ANYD(cxfab), + BL_TO_FORTRAN_ANYD(cyfab), + BL_TO_FORTRAN_ANYD(czfab), + dx, &r_ratio, &use_limiter); +#else + amrex_interp_div_free_bfield(ccbx.loVect(), ccbx.hiVect(), + BL_TO_FORTRAN_ANYD(ffab[0]), + BL_TO_FORTRAN_ANYD(ffab[2]), + BL_TO_FORTRAN_ANYD(cxfab), + BL_TO_FORTRAN_ANYD(czfab), + dx, &r_ratio, &use_limiter); + amrex_interp_cc_bfield(ccbx.loVect(), ccbx.hiVect(), + BL_TO_FORTRAN_ANYD(ffab[1]), + BL_TO_FORTRAN_ANYD(cyfab), + &r_ratio, &use_limiter); +#endif + // - Edge centered, in the same way as E on a Yee grid + } else if ( (*Fx_fp)[mfi].box().type() == IntVect{AMREX_D_DECL(0,1,1)} ){ +#if (AMREX_SPACEDIM == 3) + amrex_interp_efield(ccbx.loVect(), ccbx.hiVect(), + BL_TO_FORTRAN_ANYD(ffab[0]), + BL_TO_FORTRAN_ANYD(ffab[1]), + BL_TO_FORTRAN_ANYD(ffab[2]), + BL_TO_FORTRAN_ANYD(cxfab), + BL_TO_FORTRAN_ANYD(cyfab), + BL_TO_FORTRAN_ANYD(czfab), + &r_ratio, &use_limiter); +#else + amrex_interp_efield(ccbx.loVect(), ccbx.hiVect(), + BL_TO_FORTRAN_ANYD(ffab[0]), + BL_TO_FORTRAN_ANYD(ffab[2]), + BL_TO_FORTRAN_ANYD(cxfab), + BL_TO_FORTRAN_ANYD(czfab), + &r_ratio,&use_limiter); + amrex_interp_nd_efield(ccbx.loVect(), ccbx.hiVect(), + BL_TO_FORTRAN_ANYD(ffab[1]), + BL_TO_FORTRAN_ANYD(cyfab), + &r_ratio); +#endif + } else { + amrex::Abort("Unknown field staggering."); + } + + // Add temporary array to the returned structure + for (int i = 0; i < 3; ++i) { + const Box& bx = (*interpolated_F[i])[mfi].box(); + (*interpolated_F[i])[mfi].plus(ffab[i], bx, bx, 0, 0, 1); + } + } + } + return interpolated_F; +} + +std::array<std::unique_ptr<MultiFab>, 3> WarpX::getInterpolatedE(int lev) const +{ + + const int ngrow = 0; + const DistributionMapping& dm = DistributionMap(lev); + const Real* dx = Geom(lev-1).CellSize(); + const int r_ratio = refRatio(lev-1)[0]; + + return getInterpolatedVector( + Efield_cp[lev][0], Efield_cp[lev][1], Efield_cp[lev][2], + Efield_fp[lev][0], Efield_fp[lev][1], Efield_fp[lev][2], + dm, r_ratio, dx, ngrow ); +} + +std::array<std::unique_ptr<MultiFab>, 3> WarpX::getInterpolatedB(int lev) const +{ + const int ngrow = 0; + const DistributionMapping& dm = DistributionMap(lev); + const Real* dx = Geom(lev-1).CellSize(); + const int r_ratio = refRatio(lev-1)[0]; + + return getInterpolatedVector( + Bfield_cp[lev][0], Bfield_cp[lev][1], Bfield_cp[lev][2], + Bfield_fp[lev][0], Bfield_fp[lev][1], Bfield_fp[lev][2], + dm, r_ratio, dx, ngrow ); +} diff --git a/Source/Diagnostics/WarpXIO.cpp b/Source/Diagnostics/WarpXIO.cpp index 0d65b2796..ee1991c21 100644 --- a/Source/Diagnostics/WarpXIO.cpp +++ b/Source/Diagnostics/WarpXIO.cpp @@ -1,4 +1,3 @@ - #include <AMReX_MultiFabUtil.H> #include <AMReX_PlotFileUtil.H> #include <AMReX_FillPatchUtil_F.H> @@ -181,7 +180,6 @@ WarpX::WriteCheckPointFile() const VisMF::SetHeaderVersion(current_version); } - void WarpX::InitFromCheckpoint () { @@ -514,130 +512,66 @@ WarpX::WritePlotFile () const const int nlevels = finestLevel()+1; for (int lev = 0; lev < nlevels; ++lev) { - const std::string raw_plotfilename = plotfilename + "/raw_fields"; - // Plot auxilary patch - if (plot_raw_fields_guards) { - VisMF::Write(*Efield_aux[lev][0], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ex_aux")); - VisMF::Write(*Efield_aux[lev][1], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ey_aux")); - VisMF::Write(*Efield_aux[lev][2], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ez_aux")); - VisMF::Write(*Bfield_aux[lev][0], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Bx_aux")); - VisMF::Write(*Bfield_aux[lev][1], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "By_aux")); - VisMF::Write(*Bfield_aux[lev][2], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Bz_aux")); - } else { - const DistributionMapping& dm = DistributionMap(lev); - MultiFab Ex(Efield_aux[lev][0]->boxArray(), dm, 1, 0); - MultiFab Ey(Efield_aux[lev][1]->boxArray(), dm, 1, 0); - MultiFab Ez(Efield_aux[lev][2]->boxArray(), dm, 1, 0); - MultiFab Bx(Bfield_aux[lev][0]->boxArray(), dm, 1, 0); - MultiFab By(Bfield_aux[lev][1]->boxArray(), dm, 1, 0); - MultiFab Bz(Bfield_aux[lev][2]->boxArray(), dm, 1, 0); - MultiFab::Copy(Ex, *Efield_aux[lev][0], 0, 0, 1, 0); - MultiFab::Copy(Ey, *Efield_aux[lev][1], 0, 0, 1, 0); - MultiFab::Copy(Ez, *Efield_aux[lev][2], 0, 0, 1, 0); - MultiFab::Copy(Bx, *Bfield_aux[lev][0], 0, 0, 1, 0); - MultiFab::Copy(By, *Bfield_aux[lev][1], 0, 0, 1, 0); - MultiFab::Copy(Bz, *Bfield_aux[lev][2], 0, 0, 1, 0); - VisMF::Write(Ex, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ex_aux")); - VisMF::Write(Ey, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ey_aux")); - VisMF::Write(Ez, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ez_aux")); - VisMF::Write(Bx, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Bx_aux")); - VisMF::Write(By, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "By_aux")); - VisMF::Write(Bz, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Bz_aux")); - } - - // Plot fine patch + const std::unique_ptr<MultiFab> empty_ptr; + const std::string raw_pltname = plotfilename + "/raw_fields"; + const DistributionMapping& dm = DistributionMap(lev); + + // Auxiliary patch + WriteRawField( *Efield_aux[lev][0], dm, raw_pltname, level_prefix, "Ex_aux", lev, plot_raw_fields_guards); + WriteRawField( *Efield_aux[lev][1], dm, raw_pltname, level_prefix, "Ey_aux", lev, plot_raw_fields_guards); + WriteRawField( *Efield_aux[lev][2], dm, raw_pltname, level_prefix, "Ez_aux", lev, plot_raw_fields_guards); + WriteRawField( *Bfield_aux[lev][0], dm, raw_pltname, level_prefix, "Bx_aux", lev, plot_raw_fields_guards); + WriteRawField( *Bfield_aux[lev][1], dm, raw_pltname, level_prefix, "By_aux", lev, plot_raw_fields_guards); + WriteRawField( *Bfield_aux[lev][2], dm, raw_pltname, level_prefix, "Bz_aux", lev, plot_raw_fields_guards); + + // Fine patch if (plot_finepatch) { - if (plot_raw_fields_guards) { - VisMF::Write(*Efield_fp[lev][0], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ex_fp")); - VisMF::Write(*Efield_fp[lev][1], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ey_fp")); - VisMF::Write(*Efield_fp[lev][2], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ez_fp")); - VisMF::Write(*Bfield_fp[lev][0], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Bx_fp")); - VisMF::Write(*Bfield_fp[lev][1], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "By_fp")); - VisMF::Write(*Bfield_fp[lev][2], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Bz_fp")); - VisMF::Write(*current_fp[lev][0], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "jx_fp")); - VisMF::Write(*current_fp[lev][1], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "jy_fp")); - VisMF::Write(*current_fp[lev][2], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "jz_fp")); - } else { - const DistributionMapping& dm = DistributionMap(lev); - MultiFab Ex(Efield_fp[lev][0]->boxArray(), dm, 1, 0); - MultiFab Ey(Efield_fp[lev][1]->boxArray(), dm, 1, 0); - MultiFab Ez(Efield_fp[lev][2]->boxArray(), dm, 1, 0); - MultiFab Bx(Bfield_fp[lev][0]->boxArray(), dm, 1, 0); - MultiFab By(Bfield_fp[lev][1]->boxArray(), dm, 1, 0); - MultiFab Bz(Bfield_fp[lev][2]->boxArray(), dm, 1, 0); - MultiFab jx(current_fp[lev][0]->boxArray(), dm, 1, 0); - MultiFab jy(current_fp[lev][1]->boxArray(), dm, 1, 0); - MultiFab jz(current_fp[lev][2]->boxArray(), dm, 1, 0); - MultiFab::Copy(Ex, *Efield_fp[lev][0], 0, 0, 1, 0); - MultiFab::Copy(Ey, *Efield_fp[lev][1], 0, 0, 1, 0); - MultiFab::Copy(Ez, *Efield_fp[lev][2], 0, 0, 1, 0); - MultiFab::Copy(Bx, *Bfield_fp[lev][0], 0, 0, 1, 0); - MultiFab::Copy(By, *Bfield_fp[lev][1], 0, 0, 1, 0); - MultiFab::Copy(Bz, *Bfield_fp[lev][2], 0, 0, 1, 0); - MultiFab::Copy(jx, *current_fp[lev][0], 0, 0, 1, 0); - MultiFab::Copy(jy, *current_fp[lev][1], 0, 0, 1, 0); - MultiFab::Copy(jz, *current_fp[lev][2], 0, 0, 1, 0); - VisMF::Write(Ex, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ex_fp")); - VisMF::Write(Ey, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ey_fp")); - VisMF::Write(Ez, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ez_fp")); - VisMF::Write(Bx, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Bx_fp")); - VisMF::Write(By, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "By_fp")); - VisMF::Write(Bz, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Bz_fp")); - VisMF::Write(jx, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "jx_fp")); - VisMF::Write(jy, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "jy_fp")); - VisMF::Write(jz, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "jz_fp")); + WriteRawField( *Efield_fp[lev][0], dm, raw_pltname, level_prefix, "Ex_fp", lev, plot_raw_fields_guards); + WriteRawField( *Efield_fp[lev][1], dm, raw_pltname, level_prefix, "Ey_fp", lev, plot_raw_fields_guards); + WriteRawField( *Efield_fp[lev][2], dm, raw_pltname, level_prefix, "Ez_fp", lev, plot_raw_fields_guards); + WriteRawField( *current_fp[lev][0], dm, raw_pltname, level_prefix, "jx_fp", lev, plot_raw_fields_guards); + WriteRawField( *current_fp[lev][1], dm, raw_pltname, level_prefix, "jy_fp", lev, plot_raw_fields_guards); + WriteRawField( *current_fp[lev][2], dm, raw_pltname, level_prefix, "jz_fp", lev, plot_raw_fields_guards); + WriteRawField( *Bfield_fp[lev][0], dm, raw_pltname, level_prefix, "Bx_fp", lev, plot_raw_fields_guards); + WriteRawField( *Bfield_fp[lev][1], dm, raw_pltname, level_prefix, "By_fp", lev, plot_raw_fields_guards); + WriteRawField( *Bfield_fp[lev][2], dm, raw_pltname, level_prefix, "Bz_fp", lev, plot_raw_fields_guards); + 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); + WriteRawField( rho_new, dm, raw_pltname, level_prefix, "rho_fp", lev, plot_raw_fields_guards); } } - // Plot coarse patch - if (plot_crsepatch) - { - if (lev == 0) - { - const DistributionMapping& dm = DistributionMap(lev); - MultiFab Ex(Efield_aux[lev][0]->boxArray(), dm, 1, 0); - MultiFab Ey(Efield_aux[lev][1]->boxArray(), dm, 1, 0); - MultiFab Ez(Efield_aux[lev][2]->boxArray(), dm, 1, 0); - MultiFab Bx(Bfield_aux[lev][0]->boxArray(), dm, 1, 0); - MultiFab By(Bfield_aux[lev][1]->boxArray(), dm, 1, 0); - MultiFab Bz(Bfield_aux[lev][2]->boxArray(), dm, 1, 0); - - Ex.setVal(0.0); Ey.setVal(0.0); Ez.setVal(0.0); - Bx.setVal(0.0); By.setVal(0.0); Bz.setVal(0.0); - - VisMF::Write(Ex, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ex_cp")); - VisMF::Write(Ey, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ey_cp")); - VisMF::Write(Ez, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ez_cp")); - VisMF::Write(Bx, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Bx_cp")); - VisMF::Write(By, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "By_cp")); - VisMF::Write(Bz, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Bz_cp")); - } else { - - if (plot_raw_fields_guards) { - VisMF::Write(*Efield_cp[lev][0], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ex_cp")); - VisMF::Write(*Efield_cp[lev][1], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ey_cp")); - VisMF::Write(*Efield_cp[lev][2], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ez_cp")); - VisMF::Write(*Bfield_cp[lev][0], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Bx_cp")); - VisMF::Write(*Bfield_cp[lev][1], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "By_cp")); - VisMF::Write(*Bfield_cp[lev][2], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Bz_cp")); - } else { - std::array<std::unique_ptr<MultiFab>, 3> E = getInterpolatedE(lev); - std::array<std::unique_ptr<MultiFab>, 3> B = getInterpolatedB(lev); - - VisMF::Write(*E[0], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ex_cp")); - VisMF::Write(*E[1], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ey_cp")); - VisMF::Write(*E[2], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ez_cp")); - VisMF::Write(*B[0], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Bx_cp")); - VisMF::Write(*B[1], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "By_cp")); - VisMF::Write(*B[2], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Bz_cp")); - } - } + // Coarse path + if (plot_crsepatch) { + const Real* dx = Geom(lev-1).CellSize(); + const int r_ratio = refRatio(lev-1)[0]; + WriteCoarseVector( "E", + Efield_cp[lev][0], Efield_cp[lev][1], Efield_cp[lev][2], + Efield_fp[lev][0], Efield_fp[lev][1], Efield_fp[lev][2], + dm, raw_pltname, level_prefix, lev, plot_raw_fields_guards, + r_ratio, dx ); + WriteCoarseVector( "B", + Bfield_cp[lev][0], Bfield_cp[lev][1], Bfield_cp[lev][2], + Bfield_fp[lev][0], Bfield_fp[lev][1], Bfield_fp[lev][2], + dm, raw_pltname, level_prefix, lev, plot_raw_fields_guards, + r_ratio, dx ); + WriteCoarseVector( "j", + current_cp[lev][0], current_cp[lev][1], current_cp[lev][2], + current_fp[lev][0], current_fp[lev][1], current_fp[lev][2], + dm, raw_pltname, level_prefix, lev, plot_raw_fields_guards, + r_ratio, dx ); + if (F_cp[lev]) WriteCoarseScalar( + "F", F_cp[lev], F_fp[lev], + dm, raw_pltname, level_prefix, lev, + plot_raw_fields_guards, r_ratio, dx ); + if (plot_rho) WriteCoarseScalar( + "rho", rho_cp[lev], rho_fp[lev], + dm, raw_pltname, level_prefix, lev, + plot_raw_fields_guards, r_ratio, dx, 1 ); + // Use the component 1 of `rho_cp`, i.e. rho_new for time synchronization } - - if (F_fp[lev]) { - VisMF::Write(*F_fp[lev], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "F_fp")); - } - } } @@ -793,130 +727,3 @@ WarpX::WriteJobInfo (const std::string& dir) const jobInfoFile.close(); } } - -std::array<std::unique_ptr<MultiFab>, 3> WarpX::getInterpolatedE(int lev) const -{ - - const int ngrow = 0; - - std::array<std::unique_ptr<MultiFab>, 3> interpolated_E; - for (int i = 0; i < 3; ++i) { - interpolated_E[i].reset( new MultiFab(Efield_aux[lev][i]->boxArray(), dmap[lev], 1, ngrow) ); - interpolated_E[i]->setVal(0.0); - } - - const int r_ratio = refRatio(lev-1)[0]; - const int use_limiter = 0; -#ifdef _OPEMP -#pragma omp parallel -#endif - { - std::array<FArrayBox,3> efab; - for (MFIter mfi(*interpolated_E[0]); mfi.isValid(); ++mfi) - { - Box ccbx = mfi.fabbox(); - ccbx.enclosedCells(); - ccbx.coarsen(r_ratio).refine(r_ratio); // so that ccbx is coarsenable - - const FArrayBox& cxfab = (*Efield_cp[lev][0])[mfi]; - const FArrayBox& cyfab = (*Efield_cp[lev][1])[mfi]; - const FArrayBox& czfab = (*Efield_cp[lev][2])[mfi]; - - efab[0].resize(amrex::convert(ccbx,Ex_nodal_flag)); - efab[1].resize(amrex::convert(ccbx,Ey_nodal_flag)); - efab[2].resize(amrex::convert(ccbx,Ez_nodal_flag)); - -#if (AMREX_SPACEDIM == 3) - amrex_interp_efield(ccbx.loVect(), ccbx.hiVect(), - BL_TO_FORTRAN_ANYD(efab[0]), - BL_TO_FORTRAN_ANYD(efab[1]), - BL_TO_FORTRAN_ANYD(efab[2]), - BL_TO_FORTRAN_ANYD(cxfab), - BL_TO_FORTRAN_ANYD(cyfab), - BL_TO_FORTRAN_ANYD(czfab), - &r_ratio,&use_limiter); -#else - amrex_interp_efield(ccbx.loVect(), ccbx.hiVect(), - BL_TO_FORTRAN_ANYD(efab[0]), - BL_TO_FORTRAN_ANYD(efab[2]), - BL_TO_FORTRAN_ANYD(cxfab), - BL_TO_FORTRAN_ANYD(czfab), - &r_ratio,&use_limiter); - amrex_interp_nd_efield(ccbx.loVect(), ccbx.hiVect(), - BL_TO_FORTRAN_ANYD(efab[1]), - BL_TO_FORTRAN_ANYD(cyfab), - &r_ratio); -#endif - - for (int i = 0; i < 3; ++i) { - const Box& bx = (*interpolated_E[i])[mfi].box(); - (*interpolated_E[i])[mfi].plus(efab[i], bx, bx, 0, 0, 1); - } - } - } - - return interpolated_E; -} - -std::array<std::unique_ptr<MultiFab>, 3> WarpX::getInterpolatedB(int lev) const -{ - const int ngrow = 0; - - std::array<std::unique_ptr<MultiFab>, 3> interpolated_B; - for (int i = 0; i < 3; ++i) { - interpolated_B[i].reset( new MultiFab(Bfield_aux[lev][i]->boxArray(), dmap[lev], 1, ngrow) ); - interpolated_B[i]->setVal(0.0); - } - - const Real* dx = Geom(lev-1).CellSize(); - const int r_ratio = refRatio(lev-1)[0]; - const int use_limiter = 0; -#ifdef _OPEMP -#pragma omp parallel -#endif - { - std::array<FArrayBox,3> bfab; - for (MFIter mfi(*interpolated_B[0]); mfi.isValid(); ++mfi) - { - Box ccbx = mfi.fabbox(); - ccbx.enclosedCells(); - ccbx.coarsen(r_ratio).refine(r_ratio); // so that ccbx is coarsenable - - const FArrayBox& cxfab = (*Bfield_cp[lev][0])[mfi]; - const FArrayBox& cyfab = (*Bfield_cp[lev][1])[mfi]; - const FArrayBox& czfab = (*Bfield_cp[lev][2])[mfi]; - - bfab[0].resize(amrex::convert(ccbx,Bx_nodal_flag)); - bfab[1].resize(amrex::convert(ccbx,By_nodal_flag)); - bfab[2].resize(amrex::convert(ccbx,Bz_nodal_flag)); - -#if (AMREX_SPACEDIM == 3) - amrex_interp_div_free_bfield(ccbx.loVect(), ccbx.hiVect(), - BL_TO_FORTRAN_ANYD(bfab[0]), - BL_TO_FORTRAN_ANYD(bfab[1]), - BL_TO_FORTRAN_ANYD(bfab[2]), - BL_TO_FORTRAN_ANYD(cxfab), - BL_TO_FORTRAN_ANYD(cyfab), - BL_TO_FORTRAN_ANYD(czfab), - dx, &r_ratio, &use_limiter); -#else - amrex_interp_div_free_bfield(ccbx.loVect(), ccbx.hiVect(), - BL_TO_FORTRAN_ANYD(bfab[0]), - BL_TO_FORTRAN_ANYD(bfab[2]), - BL_TO_FORTRAN_ANYD(cxfab), - BL_TO_FORTRAN_ANYD(czfab), - dx, &r_ratio, &use_limiter); - amrex_interp_cc_bfield(ccbx.loVect(), ccbx.hiVect(), - BL_TO_FORTRAN_ANYD(bfab[1]), - BL_TO_FORTRAN_ANYD(cyfab), - &r_ratio, &use_limiter); -#endif - - for (int i = 0; i < 3; ++i) { - const Box& bx = (*interpolated_B[i])[mfi].box(); - (*interpolated_B[i])[mfi].plus(bfab[i], bx, bx, 0, 0, 1); - } - } - } - return interpolated_B; -} |