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_cc_vel(finest_level+1);
195 for (
int lev = 0; lev <= finest_level; ++lev) {
196 mf_cc_vel[lev].define(grids[lev], dmap[lev], AMREX_SPACEDIM, IntVect(1,1,0));
197 mf_cc_vel[lev].setVal(0.0_rt);
198 average_face_to_cellcenter(mf_cc_vel[lev],0,
200 mf_cc_vel[lev].FillBoundary(geom[lev].periodicity());
204 amrex::Interpolater* mapper = &cell_cons_interp;
206 for (
int lev = 1; lev <= finest_level; ++lev) {
207 Vector<MultiFab*> fmf = {&(mf_cc_vel[lev]), &(mf_cc_vel[lev])};
208 Vector<Real> ftime = {
t_new[lev],
t_new[lev]};
209 Vector<MultiFab*> cmf = {&mf_cc_vel[lev-1], &mf_cc_vel[lev-1]};
210 Vector<Real> ctime = {
t_new[lev],
t_new[lev]};
213 amrex::FillPatchTwoLevels(mf_cc_vel[lev],
t_new[lev], cmf, ctime, fmf, ftime,
214 0, 0, AMREX_SPACEDIM, geom[lev-1], geom[lev],
221 for (
int lev = 0; lev <= finest_level; ++lev)
227 for (
int i = 0; i <
NCONS; ++i) {
230 amrex::Abort(
"Found while writing output: Cons (salt, temp, or scalar, etc) contains nan or inf");
232 MultiFab::Copy(mf[lev],*
cons_new[lev],i,mf_comp,1,ngrow_vars);
239 if (mf_cc_vel[lev].contains_nan(0,1) || mf_cc_vel[lev].contains_inf(0,1)) {
240 amrex::Abort(
"Found while writing output: u velocity contains nan or inf");
242 MultiFab::Copy(mf[lev], mf_cc_vel[lev], 0, mf_comp, 1, 0);
246 if (mf_cc_vel[lev].contains_nan(1,1) || mf_cc_vel[lev].contains_inf(1,1)) {
247 amrex::Abort(
"Found while writing output: v velocity contains nan or inf");
249 MultiFab::Copy(mf[lev], mf_cc_vel[lev], 1, mf_comp, 1, 0);
253 if (mf_cc_vel[lev].contains_nan(2,1) || mf_cc_vel[lev].contains_inf(2,1)) {
254 amrex::Abort(
"Found while writing output: z velocity contains nan or inf");
256 MultiFab::Copy(mf[lev], mf_cc_vel[lev], 2, mf_comp, 1, 0);
261 auto calculate_derived = [&](
const std::string& der_name,
265 MultiFab dmf(mf[lev], make_alias, mf_comp, 1);
267#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
269 for (MFIter mfi(dmf, TilingIfNotGPU()); mfi.isValid(); ++mfi)
271 const Box& bx = mfi.tilebox();
272 auto& dfab = dmf[mfi];
274 if (der_name ==
"vorticity") {
275 auto const& sfab = mf_cc_vel[lev][mfi];
276 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);
278 auto const& sfab = (*
cons_new[lev])[mfi];
279 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);
291 Real dx = Geom()[lev].CellSizeArray()[0];
292 Real dy = Geom()[lev].CellSizeArray()[1];
299 MultiFab dmf(mf[lev], make_alias, mf_comp, AMREX_SPACEDIM);
301#pragma omp parallel if (Gpu::notInLaunchRegion())
303 for (MFIter mfi(dmf, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
304 const Box& bx = mfi.tilebox();
305 const Array4<Real> loc_arr = dmf.array(mfi);
306 const Array4<Real const> zp_arr =
vec_z_phys_nd[lev]->const_array(mfi);
308 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) {
309 loc_arr(i,j,k,0) = (i+0.5_rt) * dx;
310 loc_arr(i,j,k,1) = (j+0.5_rt) * dy;
311 loc_arr(i,j,k,2) = 0.125_rt * (zp_arr(i,j ,k ) + zp_arr(i+1,j ,k ) +
312 zp_arr(i,j+1,k ) + zp_arr(i+1,j+1,k ) +
313 zp_arr(i,j ,k+1) + zp_arr(i+1,j ,k+1) +
314 zp_arr(i,j+1,k+1) + zp_arr(i+1,j+1,k+1) );
317 mf_comp += AMREX_SPACEDIM;
320#ifdef REMORA_USE_PARTICLES
321 const auto& particles_namelist( particleData.getNames() );
322 for (ParticlesNamesVector::size_type i = 0; i < particles_namelist.size(); i++) {
324 MultiFab temp_dat(mf[lev].boxArray(), mf[lev].DistributionMap(), 1, 0);
326 particleData[particles_namelist[i]]->Increment(temp_dat, lev);
327 MultiFab::Copy(mf[lev], temp_dat, 0, mf_comp, 1, 0);
332 Vector<std::string> particle_mesh_plot_names(0);
333 particleData.GetMeshPlotVarNames( particle_mesh_plot_names );
334 for (
int i = 0; i < particle_mesh_plot_names.size(); i++) {
335 std::string plot_var_name(particle_mesh_plot_names[i]);
337 MultiFab temp_dat(mf[lev].boxArray(), mf[lev].DistributionMap(), 1, 1);
339 particleData.GetMeshPlotVar(plot_var_name, temp_dat, lev);
340 MultiFab::Copy(mf[lev], temp_dat, 0, mf_comp, 1, 0);
347 Real dz = Geom()[lev].CellSizeArray()[2];
348 int N = Geom()[lev].Domain().size()[2];
351#pragma omp parallel if (Gpu::notInLaunchRegion())
353 for (MFIter mfi(mf_nd[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
355 const Box& bx = mfi.tilebox();
356 Array4<Real> mf_arr = mf_nd[lev].array(mfi);
357 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) {
358 mf_arr(i,j,k,2) = mf_arr(i,j,k,2) + (N-k) * dz;
366 if (finest_level == 0)
369 amrex::Print() <<
"Writing plotfile " << plotfilename <<
"\n";
371 GetVecOfConstPtrs(mf),
372 GetVecOfConstPtrs(mf_nd),
377#ifdef REMORA_USE_PARTICLES
378 particleData.Checkpoint(plotfilename);
381#ifdef REMORA_USE_HDF5
383 amrex::Print() <<
"Writing plotfile " << plotfilename+
"d01.h5" <<
"\n";
384 WriteMultiLevelPlotfileHDF5(plotfilename, finest_level+1,
385 GetVecOfConstPtrs(mf),
390 amrex::Abort(
"User specified unknown plot_filetype");
395 Vector<IntVect> r2(finest_level);
396 Vector<Geometry> g2(finest_level+1);
397 Vector<MultiFab> mf2(finest_level+1);
399 mf2[0].define(grids[0], dmap[0], ncomp_mf, 0);
402 MultiFab::Copy(mf2[0],mf[0],0,0,mf[0].nComp(),0);
405 Array<int,AMREX_SPACEDIM> periodicity =
406 {Geom()[0].isPeriodic(0),Geom()[0].isPeriodic(1),Geom()[0].isPeriodic(2)};
407 g2[0].define(Geom()[0].Domain(),&(Geom()[0].ProbDomain()),0,periodicity.data());
410 r2[0] = IntVect(1,1,ref_ratio[0][0]);
411 for (
int lev = 1; lev <= finest_level; ++lev) {
415 r2[lev-1][2] = r2[lev-2][2] * ref_ratio[lev-1][0];
418 mf2[lev].define(refine(grids[lev],r2[lev-1]), dmap[lev], ncomp_mf, 0);
421 Box d2(Geom()[lev].Domain());
422 d2.refine(r2[lev-1]);
424 g2[lev].define(d2,&(Geom()[lev].ProbDomain()),0,periodicity.data());
429 amrex::Vector<amrex::BCRec> null_dom_bcs;
430 null_dom_bcs.resize(mf2[0].nComp());
431 for (
int n = 0; n < mf2[0].nComp(); n++) {
432 for (
int dir = 0; dir < AMREX_SPACEDIM; dir++) {
439 for (
int lev = 1; lev <= finest_level; ++lev) {
440 Interpolater* mapper_c = &pc_interp;
441 InterpFromCoarseLevel(mf2[lev],
t_new[lev], mf[lev],
442 0, 0, mf2[lev].nComp(),
445 r2[lev-1], mapper_c, null_dom_bcs, 0);
449 Vector<IntVect> rr(finest_level);
450 for (
int lev = 0; lev < finest_level; ++lev) {
451 rr[lev] = IntVect(ref_ratio[lev][0],ref_ratio[lev][1],ref_ratio[lev][0]);
454 WriteMultiLevelPlotfile(plotfilename, finest_level+1, GetVecOfConstPtrs(mf2), varnames,
458#ifdef REMORA_USE_PARTICLES
459 particleData.Checkpoint(plotfilename);
463 for (
int lev = 0; lev <= finest_level; ++lev) {
483 const Vector<const MultiFab*>& mf,
484 const Vector<const MultiFab*>& mf_nd,
485 const Vector<std::string>& varnames,
487 const Vector<int>& level_steps,
488 const std::string &versionName,
489 const std::string &levelPrefix,
490 const std::string &mfPrefix,
491 const Vector<std::string>& extra_dirs)
const
493 BL_PROFILE(
"WriteMultiLevelPlotfileWithBathymetry()");
495 BL_ASSERT(nlevels <= mf.size());
496 BL_ASSERT(nlevels <= ref_ratio.size()+1);
497 BL_ASSERT(nlevels <= level_steps.size());
498 BL_ASSERT(mf[0]->nComp() == varnames.size());
500 bool callBarrier(
false);
501 PreBuildDirectorHierarchy(plotfilename, levelPrefix, nlevels, callBarrier);
502 if (!extra_dirs.empty()) {
503 for (
const auto& d : extra_dirs) {
504 const std::string ed = plotfilename+
"/"+d;
505 PreBuildDirectorHierarchy(ed, levelPrefix, nlevels, callBarrier);
508 ParallelDescriptor::Barrier();
510 if (ParallelDescriptor::MyProc() == ParallelDescriptor::NProcs()-1) {
511 Vector<BoxArray> boxArrays(nlevels);
512 for(
int level(0); level < boxArrays.size(); ++level) {
513 boxArrays[level] = mf[level]->boxArray();
517 VisMF::IO_Buffer io_buffer(VisMF::IO_Buffer_Size);
518 std::string HeaderFileName(plotfilename +
"/Header");
519 std::ofstream HeaderFile;
520 HeaderFile.rdbuf()->pubsetbuf(io_buffer.dataPtr(), io_buffer.size());
521 HeaderFile.open(HeaderFileName.c_str(), std::ofstream::out |
522 std::ofstream::trunc |
523 std::ofstream::binary);
524 if( ! HeaderFile.good()) FileOpenFailed(HeaderFileName);
526 time, level_steps, versionName,
527 levelPrefix, mfPrefix);
530 if (AsyncOut::UseAsyncOut()) {
531 AsyncOut::Submit(std::move(f));
537 std::string mf_nodal_prefix =
"Nu_nd";
538 for (
int level = 0; level <= finest_level; ++level)
540 if (AsyncOut::UseAsyncOut()) {
541 VisMF::AsyncWrite(*mf[level],
542 MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mfPrefix),
544 VisMF::AsyncWrite(*mf_nd[level],
545 MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mf_nodal_prefix),
548 const MultiFab* data;
549 std::unique_ptr<MultiFab> mf_tmp;
550 if (mf[level]->nGrowVect() != 0) {
551 mf_tmp = std::make_unique<MultiFab>(mf[level]->boxArray(),
552 mf[level]->DistributionMap(),
553 mf[level]->nComp(), 0, MFInfo(),
554 mf[level]->Factory());
555 MultiFab::Copy(*mf_tmp, *mf[level], 0, 0, mf[level]->nComp(), 0);
560 VisMF::Write(*data , MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mfPrefix));
561 VisMF::Write(*mf_nd[level], MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mf_nodal_prefix));
580 const Vector<BoxArray> &bArray,
581 const Vector<std::string> &varnames,
583 const Vector<int> &level_steps,
584 const std::string &versionName,
585 const std::string &levelPrefix,
586 const std::string &mfPrefix)
const
588 BL_ASSERT(nlevels <= bArray.size());
589 BL_ASSERT(nlevels <= ref_ratio.size()+1);
590 BL_ASSERT(nlevels <= level_steps.size());
592 HeaderFile.precision(17);
595 HeaderFile << versionName <<
'\n';
597 HeaderFile << varnames.size() <<
'\n';
599 for (
int ivar = 0; ivar < varnames.size(); ++ivar) {
600 HeaderFile << varnames[ivar] <<
"\n";
602 HeaderFile << AMREX_SPACEDIM <<
'\n';
603 HeaderFile << time <<
'\n';
604 HeaderFile << finest_level <<
'\n';
605 for (
int i = 0; i < AMREX_SPACEDIM; ++i) {
606 HeaderFile << geom[0].ProbLo(i) <<
' ';
609 for (
int i = 0; i < AMREX_SPACEDIM; ++i) {
610 HeaderFile << geom[0].ProbHi(i) <<
' ';
613 for (
int i = 0; i < finest_level; ++i) {
614 HeaderFile << ref_ratio[i][0] <<
' ';
617 for (
int i = 0; i <= finest_level; ++i) {
618 HeaderFile << geom[i].Domain() <<
' ';
621 for (
int i = 0; i <= finest_level; ++i) {
622 HeaderFile << level_steps[i] <<
' ';
625 for (
int i = 0; i <= finest_level; ++i) {
626 for (
int k = 0; k < AMREX_SPACEDIM; ++k) {
627 HeaderFile << geom[i].CellSize()[k] <<
' ';
631 HeaderFile << (int) geom[0].
Coord() <<
'\n';
634 for (
int level = 0; level <= finest_level; ++level) {
635 HeaderFile << level <<
' ' << bArray[level].size() <<
' ' << time <<
'\n';
636 HeaderFile << level_steps[level] <<
'\n';
638 const IntVect& domain_lo = geom[level].Domain().smallEnd();
639 for (
int i = 0; i < bArray[level].size(); ++i)
644 const Box& b = shift(bArray[level][i], -domain_lo);
645 RealBox loc = RealBox(b, geom[level].CellSize(), geom[level].ProbLo());
646 for (
int n = 0; n < AMREX_SPACEDIM; ++n) {
647 HeaderFile << loc.lo(n) <<
' ' << loc.hi(n) <<
'\n';
651 HeaderFile << MultiFabHeaderPath(level, levelPrefix, mfPrefix) <<
'\n';
653 HeaderFile <<
"1" <<
"\n";
654 HeaderFile <<
"3" <<
"\n";
655 HeaderFile <<
"amrexvec_nu_x" <<
"\n";
656 HeaderFile <<
"amrexvec_nu_y" <<
"\n";
657 HeaderFile <<
"amrexvec_nu_z" <<
"\n";
658 std::string mf_nodal_prefix =
"Nu_nd";
659 for (
int level = 0; level <= finest_level; ++level) {
660 HeaderFile << MultiFabHeaderPath(level, levelPrefix, mf_nodal_prefix) <<
'\n';
671 for (MFIter mfi(*
cons_new[lev],
false); mfi.isValid(); ++mfi) {
672 Box gbx1 = mfi.growntilebox(IntVect(
NGROW+1,
NGROW+1,0));
673 Box ubx = mfi.grownnodaltilebox(0,IntVect(
NGROW,
NGROW,0));
674 Box vbx = mfi.grownnodaltilebox(1,IntVect(
NGROW,
NGROW,0));
676 Array4<Real>
const& Zt_avg1 =
vec_Zt_avg1[lev]->array(mfi);
677 Array4<Real>
const& ubar =
vec_ubar[lev]->array(mfi);
678 Array4<Real>
const& vbar =
vec_vbar[lev]->array(mfi);
679 Array4<Real>
const& xvel =
xvel_new[lev]->array(mfi);
680 Array4<Real>
const& yvel =
yvel_new[lev]->array(mfi);
684 Array4<Real const>
const& mskr =
vec_mskr[lev]->array(mfi);
685 Array4<Real const>
const& msku =
vec_msku[lev]->array(mfi);
686 Array4<Real const>
const& mskv =
vec_mskv[lev]->array(mfi);
688 ParallelFor(makeSlab(gbx1,2,0), [=] AMREX_GPU_DEVICE (
int i,
int j,
int )
691 Zt_avg1(i,j,0) = fill_value;
694 ParallelFor(gbx1, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
697 temp(i,j,k) = fill_value;
698 salt(i,j,k) = fill_value;
701 ParallelFor(makeSlab(ubx,2,0), 3, [=] AMREX_GPU_DEVICE (
int i,
int j,
int ,
int n)
703 if (!msku(i,j,0) && ubar(i,j,0)==fill_where) {
704 ubar(i,j,0,n) = fill_value;
707 ParallelFor(makeSlab(vbx,2,0), 3, [=] AMREX_GPU_DEVICE (
int i,
int j,
int ,
int n)
709 if (!mskv(i,j,0) && vbar(i,j,0)==fill_where) {
710 vbar(i,j,0,n) = fill_value;
713 ParallelFor(ubx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
715 if (!msku(i,j,0) && xvel(i,j,k)==fill_where) {
716 xvel(i,j,k) = fill_value;
719 ParallelFor(vbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
721 if (!mskv(i,j,0) && yvel(i,j,k)==fill_where) {
722 yvel(i,j,k) = fill_value;
726 Gpu::streamSynchronize();