149 Vector<std::string> varnames;
152 const int ncomp_mf = varnames.size();
153 const auto ngrow_vars = IntVect(
NGROW-1,
NGROW-1,0);
161 for (
int lev = 0; lev <= finest_level; ++lev) {
168 Real fill_value = 0.0_rt;
169 for (
int lev = 0; lev <= finest_level; ++lev) {
174 Vector<MultiFab> mf(finest_level+1);
175 for (
int lev = 0; lev <= finest_level; ++lev) {
176 mf[lev].define(grids[lev], dmap[lev], ncomp_mf, ngrow_vars);
180 Vector<MultiFab> mf_nd(finest_level+1);
181 for (
int lev = 0; lev <= finest_level; ++lev) {
182 BoxArray nodal_grids(grids[lev]); nodal_grids.surroundingNodes();
183 mf_nd[lev].define(nodal_grids, dmap[lev], AMREX_SPACEDIM, 0);
184 mf_nd[lev].setVal(0.);
188 Vector<MultiFab> mf_u(finest_level+1);
189 Vector<MultiFab> mf_v(finest_level+1);
190 Vector<MultiFab> mf_w(finest_level+1);
192 for (
int lev = 0; lev <= finest_level; ++lev) {
193 BoxArray grid_stag_u(grids[lev]); grid_stag_u.surroundingNodes(0);
194 BoxArray grid_stag_v(grids[lev]); grid_stag_v.surroundingNodes(1);
195 BoxArray grid_stag_w(grids[lev]); grid_stag_w.surroundingNodes(2);
196 mf_u[lev].define(grid_stag_u, dmap[lev], 1, 0);
197 mf_v[lev].define(grid_stag_v, dmap[lev], 1, 0);
198 mf_w[lev].define(grid_stag_w, dmap[lev], 1, 0);
199 MultiFab::Copy(mf_u[lev],*
xvel_new[lev],0,0,1,0);
200 MultiFab::Copy(mf_v[lev],*
yvel_new[lev],0,0,1,0);
201 MultiFab::Copy(mf_w[lev],*
zvel_new[lev],0,0,1,0);
206 Vector<MultiFab> mf_cc_vel(finest_level+1);
213 for (
int lev = 0; lev <= finest_level; ++lev) {
214 mf_cc_vel[lev].define(grids[lev], dmap[lev], AMREX_SPACEDIM, IntVect(1,1,0));
215 mf_cc_vel[lev].setVal(0.0_rt);
216 average_face_to_cellcenter(mf_cc_vel[lev],0,
218 mf_cc_vel[lev].FillBoundary(geom[lev].periodicity());
222 amrex::Interpolater* mapper = &cell_cons_interp;
224 for (
int lev = 1; lev <= finest_level; ++lev) {
225 Vector<MultiFab*> fmf = {&(mf_cc_vel[lev]), &(mf_cc_vel[lev])};
226 Vector<Real> ftime = {
t_new[lev],
t_new[lev]};
227 Vector<MultiFab*> cmf = {&mf_cc_vel[lev-1], &mf_cc_vel[lev-1]};
228 Vector<Real> ctime = {
t_new[lev],
t_new[lev]};
231 amrex::FillPatchTwoLevels(mf_cc_vel[lev],
t_new[lev], cmf, ctime, fmf, ftime,
232 0, 0, AMREX_SPACEDIM, geom[lev-1], geom[lev],
239 for (
int lev = 0; lev <= finest_level; ++lev)
245 for (
int i = 0; i <
NCONS; ++i) {
248 amrex::Abort(
"Found while writing output: Cons (salt, temp, or scalar, etc) contains nan or inf");
250 MultiFab::Copy(mf[lev],*
cons_new[lev],i,mf_comp,1,ngrow_vars);
257 if (mf_cc_vel[lev].contains_nan(0,1) || mf_cc_vel[lev].contains_inf(0,1)) {
258 amrex::Abort(
"Found while writing output: u velocity contains nan or inf");
260 MultiFab::Copy(mf[lev], mf_cc_vel[lev], 0, mf_comp, 1, 0);
264 if (mf_cc_vel[lev].contains_nan(1,1) || mf_cc_vel[lev].contains_inf(1,1)) {
265 amrex::Abort(
"Found while writing output: v velocity contains nan or inf");
267 MultiFab::Copy(mf[lev], mf_cc_vel[lev], 1, mf_comp, 1, 0);
271 if (mf_cc_vel[lev].contains_nan(2,1) || mf_cc_vel[lev].contains_inf(2,1)) {
272 amrex::Abort(
"Found while writing output: z velocity contains nan or inf");
274 MultiFab::Copy(mf[lev], mf_cc_vel[lev], 2, mf_comp, 1, 0);
279 auto calculate_derived = [&](
const std::string& der_name,
283 MultiFab dmf(mf[lev], make_alias, mf_comp, 1);
285#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
287 for (MFIter mfi(dmf, TilingIfNotGPU()); mfi.isValid(); ++mfi)
289 const Box& bx = mfi.tilebox();
290 auto& dfab = dmf[mfi];
292 if (der_name ==
"vorticity") {
293 auto const& sfab = mf_cc_vel[lev][mfi];
294 der_function(bx, dfab, 0, 1, sfab,
vec_pm[lev]->const_array(mfi),
vec_pn[lev]->const_array(mfi), Geom(lev),
t_new[0],
nullptr, lev);
296 auto const& sfab = (*
cons_new[lev])[mfi];
297 der_function(bx, dfab, 0, 1, sfab,
vec_pm[lev]->const_array(mfi),
vec_pn[lev]->const_array(mfi), Geom(lev),
t_new[0],
nullptr, lev);
309 Real dx = Geom()[lev].CellSizeArray()[0];
310 Real dy = Geom()[lev].CellSizeArray()[1];
317 MultiFab dmf(mf[lev], make_alias, mf_comp, AMREX_SPACEDIM);
319#pragma omp parallel if (Gpu::notInLaunchRegion())
321 for (MFIter mfi(dmf, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
322 const Box& bx = mfi.tilebox();
323 const Array4<Real> loc_arr = dmf.array(mfi);
324 const Array4<Real const> zp_arr =
vec_z_phys_nd[lev]->const_array(mfi);
326 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) {
327 loc_arr(i,j,k,0) = (i+0.5_rt) * dx;
328 loc_arr(i,j,k,1) = (j+0.5_rt) * dy;
329 loc_arr(i,j,k,2) = 0.125_rt * (zp_arr(i,j ,k ) + zp_arr(i+1,j ,k ) +
330 zp_arr(i,j+1,k ) + zp_arr(i+1,j+1,k ) +
331 zp_arr(i,j ,k+1) + zp_arr(i+1,j ,k+1) +
332 zp_arr(i,j+1,k+1) + zp_arr(i+1,j+1,k+1) );
335 mf_comp += AMREX_SPACEDIM;
338#ifdef REMORA_USE_PARTICLES
339 const auto& particles_namelist( particleData.getNames() );
340 for (ParticlesNamesVector::size_type i = 0; i < particles_namelist.size(); i++) {
342 MultiFab temp_dat(mf[lev].boxArray(), mf[lev].DistributionMap(), 1, 0);
344 particleData[particles_namelist[i]]->Increment(temp_dat, lev);
345 MultiFab::Copy(mf[lev], temp_dat, 0, mf_comp, 1, 0);
350 Vector<std::string> particle_mesh_plot_names(0);
351 particleData.GetMeshPlotVarNames( particle_mesh_plot_names );
352 for (
int i = 0; i < particle_mesh_plot_names.size(); i++) {
353 std::string plot_var_name(particle_mesh_plot_names[i]);
355 MultiFab temp_dat(mf[lev].boxArray(), mf[lev].DistributionMap(), 1, 1);
357 particleData.GetMeshPlotVar(plot_var_name, temp_dat, lev);
358 MultiFab::Copy(mf[lev], temp_dat, 0, mf_comp, 1, 0);
365 Real dz = Geom()[lev].CellSizeArray()[2];
366 int N = Geom()[lev].Domain().size()[2];
369#pragma omp parallel if (Gpu::notInLaunchRegion())
371 for (MFIter mfi(mf_nd[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
373 const Box& bx = mfi.tilebox();
374 Array4<Real> mf_arr = mf_nd[lev].array(mfi);
375 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) {
376 mf_arr(i,j,k,2) = mf_arr(i,j,k,2) + (N-k) * dz;
384 if (finest_level == 0)
387 amrex::Print() <<
"Writing plotfile " << plotfilename <<
"\n";
389 GetVecOfConstPtrs(mf),
390 GetVecOfConstPtrs(mf_nd),
391 GetVecOfConstPtrs(mf_u),
392 GetVecOfConstPtrs(mf_v),
393 GetVecOfConstPtrs(mf_w),
398#ifdef REMORA_USE_PARTICLES
399 particleData.Checkpoint(plotfilename);
402#ifdef REMORA_USE_HDF5
404 amrex::Print() <<
"Writing plotfile " << plotfilename+
"d01.h5" <<
"\n";
405 WriteMultiLevelPlotfileHDF5(plotfilename, finest_level+1,
406 GetVecOfConstPtrs(mf),
411 amrex::Abort(
"User specified unknown plot_filetype");
416 Vector<IntVect> r2(finest_level);
417 Vector<Geometry> g2(finest_level+1);
418 Vector<MultiFab> mf2(finest_level+1);
420 mf2[0].define(grids[0], dmap[0], ncomp_mf, 0);
423 MultiFab::Copy(mf2[0],mf[0],0,0,mf[0].nComp(),0);
426 Array<int,AMREX_SPACEDIM> periodicity =
427 {Geom()[0].isPeriodic(0),Geom()[0].isPeriodic(1),Geom()[0].isPeriodic(2)};
428 g2[0].define(Geom()[0].Domain(),&(Geom()[0].ProbDomain()),0,periodicity.data());
431 r2[0] = IntVect(1,1,ref_ratio[0][0]);
432 for (
int lev = 1; lev <= finest_level; ++lev) {
436 r2[lev-1][2] = r2[lev-2][2] * ref_ratio[lev-1][0];
439 mf2[lev].define(refine(grids[lev],r2[lev-1]), dmap[lev], ncomp_mf, 0);
442 Box d2(Geom()[lev].Domain());
443 d2.refine(r2[lev-1]);
445 g2[lev].define(d2,&(Geom()[lev].ProbDomain()),0,periodicity.data());
450 amrex::Vector<amrex::BCRec> null_dom_bcs;
451 null_dom_bcs.resize(mf2[0].nComp());
452 for (
int n = 0; n < mf2[0].nComp(); n++) {
453 for (
int dir = 0; dir < AMREX_SPACEDIM; dir++) {
460 for (
int lev = 1; lev <= finest_level; ++lev) {
461 Interpolater* mapper_c = &pc_interp;
462 InterpFromCoarseLevel(mf2[lev],
t_new[lev], mf[lev],
463 0, 0, mf2[lev].nComp(),
466 r2[lev-1], mapper_c, null_dom_bcs, 0);
470 Vector<IntVect> rr(finest_level);
471 for (
int lev = 0; lev < finest_level; ++lev) {
472 rr[lev] = IntVect(ref_ratio[lev][0],ref_ratio[lev][1],ref_ratio[lev][0]);
475 WriteMultiLevelPlotfile(plotfilename, finest_level+1, GetVecOfConstPtrs(mf2), varnames,
479#ifdef REMORA_USE_PARTICLES
480 particleData.Checkpoint(plotfilename);
484 for (
int lev = 0; lev <= finest_level; ++lev) {
504 const Vector<const MultiFab*>& mf,
505 const Vector<const MultiFab*>& mf_nd,
506 const Vector<const MultiFab*>& mf_u,
507 const Vector<const MultiFab*>& mf_v,
508 const Vector<const MultiFab*>& mf_w,
509 const Vector<std::string>& varnames,
511 const Vector<int>& level_steps,
512 const std::string &versionName,
513 const std::string &levelPrefix,
514 const std::string &mfPrefix,
515 const Vector<std::string>& extra_dirs)
const
517 BL_PROFILE(
"WriteMultiLevelPlotfileWithBathymetry()");
519 BL_ASSERT(nlevels <= mf.size());
520 BL_ASSERT(nlevels <= ref_ratio.size()+1);
521 BL_ASSERT(nlevels <= level_steps.size());
522 BL_ASSERT(mf[0]->nComp() == varnames.size());
524 bool callBarrier(
false);
525 PreBuildDirectorHierarchy(plotfilename, levelPrefix, nlevels, callBarrier);
526 if (!extra_dirs.empty()) {
527 for (
const auto& d : extra_dirs) {
528 const std::string ed = plotfilename+
"/"+d;
529 PreBuildDirectorHierarchy(ed, levelPrefix, nlevels, callBarrier);
532 ParallelDescriptor::Barrier();
534 if (ParallelDescriptor::MyProc() == ParallelDescriptor::NProcs()-1) {
535 Vector<BoxArray> boxArrays(nlevels);
536 for(
int level(0); level < boxArrays.size(); ++level) {
537 boxArrays[level] = mf[level]->boxArray();
541 VisMF::IO_Buffer io_buffer(VisMF::IO_Buffer_Size);
542 std::string HeaderFileName(plotfilename +
"/Header");
543 std::ofstream HeaderFile;
544 HeaderFile.rdbuf()->pubsetbuf(io_buffer.dataPtr(), io_buffer.size());
545 HeaderFile.open(HeaderFileName.c_str(), std::ofstream::out |
546 std::ofstream::trunc |
547 std::ofstream::binary);
548 if( ! HeaderFile.good()) FileOpenFailed(HeaderFileName);
550 time, level_steps, versionName,
551 levelPrefix, mfPrefix);
554 if (AsyncOut::UseAsyncOut()) {
555 AsyncOut::Submit(std::move(f));
561 std::string mf_nodal_prefix =
"Nu_nd";
562 std::string mf_uface_prefix =
"UonXFace";
563 std::string mf_vface_prefix =
"VonYFace";
564 std::string mf_wface_prefix =
"WonZFace";
566 for (
int level = 0; level <= finest_level; ++level)
568 if (AsyncOut::UseAsyncOut()) {
569 VisMF::AsyncWrite(*mf[level],
570 MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mfPrefix),
572 VisMF::AsyncWrite(*mf_nd[level],
573 MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mf_nodal_prefix),
576 VisMF::AsyncWrite(*mf_u[level],
577 MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mf_uface_prefix),
579 VisMF::AsyncWrite(*mf_v[level],
580 MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mf_vface_prefix),
582 VisMF::AsyncWrite(*mf_w[level],
583 MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mf_wface_prefix),
586 }
else {
const MultiFab* data;
587 std::unique_ptr<MultiFab> mf_tmp;
588 if (mf[level]->nGrowVect() != 0) {
589 mf_tmp = std::make_unique<MultiFab>(mf[level]->boxArray(),
590 mf[level]->DistributionMap(),
591 mf[level]->nComp(), 0, MFInfo(),
592 mf[level]->Factory());
593 MultiFab::Copy(*mf_tmp, *mf[level], 0, 0, mf[level]->nComp(), 0);
598 VisMF::Write(*data , MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mfPrefix));
599 VisMF::Write(*mf_nd[level], MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mf_nodal_prefix));
601 VisMF::Write(*mf_u[level], MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mf_uface_prefix));
602 VisMF::Write(*mf_v[level], MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mf_vface_prefix));
603 VisMF::Write(*mf_w[level], MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mf_wface_prefix));
623 const Vector<BoxArray> &bArray,
624 const Vector<std::string> &varnames,
626 const Vector<int> &level_steps,
627 const std::string &versionName,
628 const std::string &levelPrefix,
629 const std::string &mfPrefix)
const
631 BL_ASSERT(nlevels <= bArray.size());
632 BL_ASSERT(nlevels <= ref_ratio.size()+1);
633 BL_ASSERT(nlevels <= level_steps.size());
635 HeaderFile.precision(17);
638 HeaderFile << versionName <<
'\n';
640 HeaderFile << varnames.size() <<
'\n';
642 for (
int ivar = 0; ivar < varnames.size(); ++ivar) {
643 HeaderFile << varnames[ivar] <<
"\n";
645 HeaderFile << AMREX_SPACEDIM <<
'\n';
646 HeaderFile << time <<
'\n';
647 HeaderFile << finest_level <<
'\n';
648 for (
int i = 0; i < AMREX_SPACEDIM; ++i) {
649 HeaderFile << geom[0].ProbLo(i) <<
' ';
652 for (
int i = 0; i < AMREX_SPACEDIM; ++i) {
653 HeaderFile << geom[0].ProbHi(i) <<
' ';
656 for (
int i = 0; i < finest_level; ++i) {
657 HeaderFile << ref_ratio[i][0] <<
' ';
660 for (
int i = 0; i <= finest_level; ++i) {
661 HeaderFile << geom[i].Domain() <<
' ';
664 for (
int i = 0; i <= finest_level; ++i) {
665 HeaderFile << level_steps[i] <<
' ';
668 for (
int i = 0; i <= finest_level; ++i) {
669 for (
int k = 0; k < AMREX_SPACEDIM; ++k) {
670 HeaderFile << geom[i].CellSize()[k] <<
' ';
674 HeaderFile << (int) geom[0].
Coord() <<
'\n';
677 for (
int level = 0; level <= finest_level; ++level) {
678 HeaderFile << level <<
' ' << bArray[level].size() <<
' ' << time <<
'\n';
679 HeaderFile << level_steps[level] <<
'\n';
681 const IntVect& domain_lo = geom[level].Domain().smallEnd();
682 for (
int i = 0; i < bArray[level].size(); ++i)
687 const Box& b = shift(bArray[level][i], -domain_lo);
688 RealBox loc = RealBox(b, geom[level].CellSize(), geom[level].ProbLo());
689 for (
int n = 0; n < AMREX_SPACEDIM; ++n) {
690 HeaderFile << loc.lo(n) <<
' ' << loc.hi(n) <<
'\n';
694 HeaderFile << MultiFabHeaderPath(level, levelPrefix, mfPrefix) <<
'\n';
696 HeaderFile <<
"1" <<
"\n";
697 HeaderFile <<
"3" <<
"\n";
698 HeaderFile <<
"amrexvec_nu_x" <<
"\n";
699 HeaderFile <<
"amrexvec_nu_y" <<
"\n";
700 HeaderFile <<
"amrexvec_nu_z" <<
"\n";
701 std::string mf_nodal_prefix =
"Nu_nd";
702 for (
int level = 0; level <= finest_level; ++level) {
703 HeaderFile << MultiFabHeaderPath(level, levelPrefix, mf_nodal_prefix) <<
'\n';
714 for (MFIter mfi(*
cons_new[lev],
false); mfi.isValid(); ++mfi) {
715 Box gbx1 = mfi.growntilebox(IntVect(
NGROW+1,
NGROW+1,0));
716 Box ubx = mfi.grownnodaltilebox(0,IntVect(
NGROW,
NGROW,0));
717 Box vbx = mfi.grownnodaltilebox(1,IntVect(
NGROW,
NGROW,0));
719 Array4<Real>
const& Zt_avg1 =
vec_Zt_avg1[lev]->array(mfi);
720 Array4<Real>
const& ubar =
vec_ubar[lev]->array(mfi);
721 Array4<Real>
const& vbar =
vec_vbar[lev]->array(mfi);
722 Array4<Real>
const& xvel =
xvel_new[lev]->array(mfi);
723 Array4<Real>
const& yvel =
yvel_new[lev]->array(mfi);
727 Array4<Real const>
const& mskr =
vec_mskr[lev]->array(mfi);
728 Array4<Real const>
const& msku =
vec_msku[lev]->array(mfi);
729 Array4<Real const>
const& mskv =
vec_mskv[lev]->array(mfi);
731 ParallelFor(makeSlab(gbx1,2,0), [=] AMREX_GPU_DEVICE (
int i,
int j,
int )
734 Zt_avg1(i,j,0) = fill_value;
737 ParallelFor(gbx1, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
740 temp(i,j,k) = fill_value;
741 salt(i,j,k) = fill_value;
744 ParallelFor(makeSlab(ubx,2,0), 3, [=] AMREX_GPU_DEVICE (
int i,
int j,
int ,
int n)
746 if (!msku(i,j,0) && ubar(i,j,0)==fill_where) {
747 ubar(i,j,0,n) = fill_value;
750 ParallelFor(makeSlab(vbx,2,0), 3, [=] AMREX_GPU_DEVICE (
int i,
int j,
int ,
int n)
752 if (!mskv(i,j,0) && vbar(i,j,0)==fill_where) {
753 vbar(i,j,0,n) = fill_value;
756 ParallelFor(ubx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
758 if (!msku(i,j,0) && xvel(i,j,k)==fill_where) {
759 xvel(i,j,k) = fill_value;
762 ParallelFor(vbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
764 if (!mskv(i,j,0) && yvel(i,j,k)==fill_where) {
765 yvel(i,j,k) = fill_value;
769 Gpu::streamSynchronize();