diff options
Diffstat (limited to 'Source/WarpXIO.cpp')
-rw-r--r-- | Source/WarpXIO.cpp | 398 |
1 files changed, 144 insertions, 254 deletions
diff --git a/Source/WarpXIO.cpp b/Source/WarpXIO.cpp index bdef4f7df..e22ef95c7 100644 --- a/Source/WarpXIO.cpp +++ b/Source/WarpXIO.cpp @@ -658,77 +658,12 @@ WarpX::WritePlotFile () const } else { - MultiFab Ex_fine(Efield_aux[lev][0]->boxArray(), dmap[lev], 1, ngrow); - MultiFab Ey_fine(Efield_aux[lev][1]->boxArray(), dmap[lev], 1, ngrow); - MultiFab Ez_fine(Efield_aux[lev][2]->boxArray(), dmap[lev], 1, ngrow); - - Ex_fine.setVal(0.0); Ey_fine.setVal(0.0); Ez_fine.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(Ex_fine); 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 - - Ex_fine[mfi].plus(efab[0], Ex_fine[mfi].box(), Ex_fine[mfi].box(), 0, 0, 1); - Ey_fine[mfi].plus(efab[1], Ey_fine[mfi].box(), Ey_fine[mfi].box(), 0, 0, 1); - Ez_fine[mfi].plus(efab[2], Ez_fine[mfi].box(), Ez_fine[mfi].box(), 0, 0, 1); - } - } - -#if (AMREX_SPACEDIM == 3) - srcmf[0] = &Ex_fine; - srcmf[1] = &Ey_fine; - srcmf[2] = &Ez_fine; -#elif (AMREX_SPACEDIM == 2) - srcmf[0] = &Ex_fine;; - srcmf[1] = &Ez_fine;; -#endif - + std::array<std::unique_ptr<MultiFab>, 3> E = getInterpolatedE(lev); + PackPlotDataPtrs(srcmf, E); amrex::average_edge_to_cellcenter(*mf[lev], dcomp, srcmf); - - VisMF::Write(*mf[lev], "test"); - #if (AMREX_SPACEDIM == 2) MultiFab::Copy(*mf[lev], *mf[lev], dcomp+1, dcomp+2, 1, ngrow); - amrex::average_node_to_cellcenter(*mf[lev], dcomp+1, Ey_fine, 0, 1); + amrex::average_node_to_cellcenter(*mf[lev], dcomp+1, *E[1], 0, 1); #endif } if (lev == 0) @@ -746,76 +681,12 @@ WarpX::WritePlotFile () const } else { - - MultiFab Bx_fine(Bfield_aux[lev][0]->boxArray(), dmap[lev], 1, ngrow); - MultiFab By_fine(Bfield_aux[lev][1]->boxArray(), dmap[lev], 1, ngrow); - MultiFab Bz_fine(Bfield_aux[lev][2]->boxArray(), dmap[lev], 1, ngrow); - - Bx_fine.setVal(0.0); By_fine.setVal(0.0); Bz_fine.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(Bx_fine); 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 - - Bx_fine[mfi].plus(bfab[0], Bx_fine[mfi].box(), Bx_fine[mfi].box(), 0, 0, 1); - By_fine[mfi].plus(bfab[1], By_fine[mfi].box(), By_fine[mfi].box(), 0, 0, 1); - Bz_fine[mfi].plus(bfab[2], Bz_fine[mfi].box(), Bz_fine[mfi].box(), 0, 0, 1); - } - } - -#if (AMREX_SPACEDIM == 3) - srcmf[0] = &Bx_fine; - srcmf[1] = &By_fine; - srcmf[2] = &Bz_fine; -#elif (AMREX_SPACEDIM == 2) - srcmf[0] = &Bx_fine; - srcmf[1] = &Bz_fine; -#endif - + std::array<std::unique_ptr<MultiFab>, 3> B = getInterpolatedB(lev); + PackPlotDataPtrs(srcmf, B); amrex::average_face_to_cellcenter(*mf[lev], dcomp, srcmf); #if (AMREX_SPACEDIM == 2) MultiFab::Copy(*mf[lev], *mf[lev], dcomp+1, dcomp+2, 1, ngrow); - MultiFab::Copy(*mf[lev], By_fine, 0, dcomp+1, 1, ngrow); + MultiFab::Copy(*mf[lev], *B[1], 0, dcomp+1, 1, ngrow); #endif } if (lev == 0) @@ -939,9 +810,8 @@ WarpX::WritePlotFile () const } // Plot coarse patch - if (plot_crsepatch) { - const int ngrow = 0; - + if (plot_crsepatch) + { if (lev == 0) { const DistributionMapping& dm = DistributionMap(lev); @@ -971,122 +841,15 @@ WarpX::WritePlotFile () const 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 { - MultiFab Ex_fine(Efield_aux[lev][0]->boxArray(), dmap[lev], 1, ngrow); - MultiFab Ey_fine(Efield_aux[lev][1]->boxArray(), dmap[lev], 1, ngrow); - MultiFab Ez_fine(Efield_aux[lev][2]->boxArray(), dmap[lev], 1, ngrow); - - Ex_fine.setVal(0.0); Ey_fine.setVal(0.0); Ez_fine.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(Ex_fine); 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 - - Ex_fine[mfi].plus(efab[0], Ex_fine[mfi].box(), Ex_fine[mfi].box(), 0, 0, 1); - Ey_fine[mfi].plus(efab[1], Ey_fine[mfi].box(), Ey_fine[mfi].box(), 0, 0, 1); - Ez_fine[mfi].plus(efab[2], Ez_fine[mfi].box(), Ez_fine[mfi].box(), 0, 0, 1); - } - } - - MultiFab Bx_fine(Bfield_aux[lev][0]->boxArray(), dmap[lev], 1, ngrow); - MultiFab By_fine(Bfield_aux[lev][1]->boxArray(), dmap[lev], 1, ngrow); - MultiFab Bz_fine(Bfield_aux[lev][2]->boxArray(), dmap[lev], 1, ngrow); - - Bx_fine.setVal(0.0); By_fine.setVal(0.0); Bz_fine.setVal(0.0); - - const Real* dx = Geom(lev-1).CellSize(); -#ifdef _OPEMP -#pragma omp parallel -#endif - { - std::array<FArrayBox,3> bfab; - for (MFIter mfi(Bx_fine); 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 - - Bx_fine[mfi].plus(bfab[0], Bx_fine[mfi].box(), Bx_fine[mfi].box(), 0, 0, 1); - By_fine[mfi].plus(bfab[1], By_fine[mfi].box(), By_fine[mfi].box(), 0, 0, 1); - Bz_fine[mfi].plus(bfab[2], Bz_fine[mfi].box(), Bz_fine[mfi].box(), 0, 0, 1); - } - } - - - VisMF::Write(Ex_fine, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ex_cp")); - VisMF::Write(Ey_fine, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ey_cp")); - VisMF::Write(Ez_fine, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ez_cp")); - VisMF::Write(Bx_fine, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Bx_cp")); - VisMF::Write(By_fine, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "By_cp")); - VisMF::Write(Bz_fine, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Bz_cp")); + 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")); } } } @@ -1251,6 +1014,133 @@ WarpX::WriteJobInfo (const std::string& dir) const } } +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; +} + void WarpX::PackPlotDataPtrs (Vector<const MultiFab*>& pmf, const std::array<std::unique_ptr<MultiFab>,3>& data) |