19 BL_PROFILE(
"Integrator::Integrator()");
27 pp_query(
"restart", restart_file_cell);
28 pp_query(
"restart_cell", restart_file_cell);
29 pp_query(
"restart_node", restart_file_node);
39 std::vector<std::string> ignore;
42 for (
unsigned int i = 0; i < ignore.size(); i++)
45 pp.remove(ignore[i].c_str());
70 nsubsteps.resize(maxLevel() + 1, 1);
71 int cnt = pp.countval(
"nsubsteps");
73 if (cnt == maxLevel()) {
75 nsubsteps.insert(nsubsteps.begin(), 1);
82 for (
int lev = 1; lev <= maxLevel(); ++lev) nsubsteps[lev] = nsubsteps_all;
85 Util::Abort(
INFO,
"number of nsubsteps input must equal either 1 or amr.max_level");
87 for (
int lev = 1; lev <= maxLevel(); ++lev)
88 nsubsteps[lev] = MaxRefRatio(lev - 1);
108 for (
int ilev = 0; ilev < maxLevel(); ++ilev)
110 std::string strlo =
"lo" + std::to_string(ilev + 1);
111 std::string strhi =
"hi" + std::to_string(ilev + 1);
116 amrex::Vector<int> lodata, hidata;
119 amrex::IntVect lo(
AMREX_D_DECL(lodata[0], lodata[1], lodata[2]));
120 amrex::IntVect hi(
AMREX_D_DECL(hidata[0], hidata[1], hidata[2]));
122 explicitmesh.box.push_back(amrex::Box(lo, hi));
127 int nlevs_max = maxLevel() + 1;
129 istep.resize(nlevs_max, 0);
131 t_new.resize(nlevs_max, 0.0);
132 t_old.resize(nlevs_max, -1.e100);
133 SetTimestep(timestep);
143 Integrator::~Integrator()
145 BL_PROFILE(
"Integrator::~Integrator");
146 if (amrex::ParallelDescriptor::IOProcessor())
152 BL_PROFILE(
"Integrator::SetTimestep");
153 int nlevs_max = maxLevel() + 1;
154 timestep = _timestep;
155 dt.resize(nlevs_max, 1.e100);
157 for (
int i = 1; i < nlevs_max; i++)
158 dt[i] = dt[i - 1] / (amrex::Real)nsubsteps[i];
160 void Integrator::SetPlotInt(
int a_plot_int)
162 BL_PROFILE(
"Integrator::SetPlotInt");
163 plot_int = a_plot_int;
171 Integrator::MakeNewLevelFromCoarse(
int lev, amrex::Real time,
const amrex::BoxArray& cgrids,
const amrex::DistributionMapping& dm)
173 BL_PROFILE(
"Integrator::MakeNewLevelFromCoarse");
175 for (
int n = 0; n < cell.number_of_fabs; n++)
177 const int ncomp = (*cell.fab_array[n])[lev - 1]->nComp();
178 const int nghost = (*cell.fab_array[n])[lev - 1]->nGrow();
180 (*cell.fab_array[n])[lev].reset(
new amrex::MultiFab(cgrids, dm, ncomp, nghost));
182 (*cell.fab_array[n])[lev]->setVal(0.0);
184 FillCoarsePatch(lev, time, *cell.fab_array[n], *cell.physbc_array[n], 0, ncomp);
187 amrex::BoxArray ngrids = cgrids;
188 ngrids.convert(amrex::IntVect::TheNodeVector());
190 for (
int n = 0; n < node.number_of_fabs; n++)
192 const int ncomp = (*node.fab_array[n])[lev - 1]->nComp();
193 const int nghost = (*node.fab_array[n])[lev - 1]->nGrow();
195 (*node.fab_array[n])[lev].reset(
new amrex::MultiFab(ngrids, dm, ncomp, nghost));
196 (*node.fab_array[n])[lev]->setVal(0.0);
198 FillCoarsePatch(lev, time, *node.fab_array[n], *node.physbc_array[n], 0, ncomp);
201 for (
unsigned int n = 0; n < m_basefields.size(); n++)
203 m_basefields[n]->MakeNewLevelFromCoarse(lev, time, cgrids, dm);
205 for (
unsigned int n = 0; n < m_basefields_cell.size(); n++)
207 m_basefields_cell[n]->MakeNewLevelFromCoarse(lev, time, cgrids, dm);
220 Integrator::RemakeLevel(
int lev,
222 const amrex::BoxArray& cgrids,
223 const amrex::DistributionMapping& dm)
225 BL_PROFILE(
"Integrator::RemakeLevel");
226 for (
int n = 0; n < cell.number_of_fabs; n++)
228 const int ncomp = (*cell.fab_array[n])[lev]->nComp();
229 const int nghost = (*cell.fab_array[n])[lev]->nGrow();
231 amrex::MultiFab new_state(cgrids, dm, ncomp, nghost);
233 new_state.setVal(0.0);
234 FillPatch(lev, time, *cell.fab_array[n], new_state, *cell.physbc_array[n], 0);
235 std::swap(new_state, *(*cell.fab_array[n])[lev]);
238 amrex::BoxArray ngrids = cgrids;
239 ngrids.convert(amrex::IntVect::TheNodeVector());
241 for (
int n = 0; n < node.number_of_fabs; n++)
243 const int ncomp = (*node.fab_array[n])[lev]->nComp();
244 const int nghost = (*node.fab_array[n])[lev]->nGrow();
246 amrex::MultiFab new_state(ngrids, dm, ncomp, nghost);
248 new_state.setVal(0.0);
249 FillPatch(lev, time, *node.fab_array[n], new_state, *node.physbc_array[n], 0);
250 std::swap(new_state, *(*node.fab_array[n])[lev]);
253 for (
unsigned int n = 0; n < m_basefields_cell.size(); n++)
255 m_basefields_cell[n]->RemakeLevel(lev, time, cgrids, dm);
257 for (
unsigned int n = 0; n < m_basefields.size(); n++)
259 m_basefields[n]->RemakeLevel(lev, time, cgrids, dm);
270 Integrator::ClearLevel(
int lev)
272 BL_PROFILE(
"Integrator::ClearLevel");
273 for (
int n = 0; n < cell.number_of_fabs; n++)
275 (*cell.fab_array[n])[lev].reset(
nullptr);
277 for (
int n = 0; n < node.number_of_fabs; n++)
279 (*node.fab_array[n])[lev].reset(
nullptr);
295 AddField<Set::Scalar, Set::Hypercube::Cell>(new_fab, new_bc, ncomp, nghost, name, writeout,
true);
301 AddField<Set::Scalar, Set::Hypercube::Cell>(new_fab,
nullptr, ncomp, 0, name, writeout,
true);
307 AddField<Set::Scalar, Set::Hypercube::Node>(new_fab, new_bc, ncomp, nghost, name, writeout,
true);
313 AddField<Set::Scalar, Set::Hypercube::Node>(new_fab,
nullptr, ncomp, nghost, name, writeout,
true);
320 Integrator::RegisterIntegratedVariable(
Set::Scalar *integrated_variable, std::string name,
bool extensive)
322 BL_PROFILE(
"Integrator::RegisterIntegratedVariable");
323 thermo.vars.push_back(integrated_variable);
324 thermo.names.push_back(name);
325 thermo.extensives.push_back(extensive);
330 Integrator::CountCells(
int lev)
332 BL_PROFILE(
"Integrator::CountCells");
333 const int N = grids[lev].size();
337 for (
int i = 0; i < N; ++i)
339 cnt += grids[lev][i].numPts();
346 Integrator::FillPatch(
int lev, amrex::Real time,
348 amrex::MultiFab& destination_mf,
351 BL_PROFILE(
"Integrator::FillPatch");
355 amrex::Vector<amrex::MultiFab*> smf;
356 smf.push_back(source_mf[lev].get());
357 amrex::Vector<amrex::Real> stime;
358 stime.push_back(time);
361 amrex::FillPatchSingleLevel(destination_mf,
367 destination_mf.nComp(),
374 amrex::Vector<amrex::MultiFab*> cmf, fmf;
375 cmf.push_back(source_mf[lev - 1].get());
376 fmf.push_back(source_mf[lev].get());
377 amrex::Vector<amrex::Real> ctime, ftime;
378 ctime.push_back(time);
379 ftime.push_back(time);
383 amrex::Interpolater* mapper;
385 if (destination_mf.boxArray().ixType() == amrex::IndexType::TheNodeType())
386 mapper = &amrex::node_bilinear_interp;
388 mapper = &amrex::cell_cons_interp;
390 amrex::Vector<amrex::BCRec> bcs(destination_mf.nComp(), physbc.
GetBCRec());
391 amrex::FillPatchTwoLevels(destination_mf, time, cmf, ctime, fmf, ftime,
392 0, icomp, destination_mf.nComp(), geom[lev - 1], geom[lev],
405 Integrator::FillCoarsePatch(
int lev,
412 BL_PROFILE(
"Integrator::FillCoarsePatch");
413 AMREX_ASSERT(lev > 0);
414 amrex::Vector<amrex::MultiFab*> cmf;
415 cmf.push_back(mf[lev - 1].get());
416 amrex::Vector<amrex::Real> ctime;
417 ctime.push_back(time);
421 amrex::Interpolater* mapper;
422 if (mf[lev]->boxArray().ixType() == amrex::IndexType::TheNodeType())
423 mapper = &amrex::node_bilinear_interp;
425 mapper = &amrex::cell_cons_interp;
427 amrex::Vector<amrex::BCRec> bcs(ncomp, physbc.
GetBCRec());
428 amrex::InterpFromCoarseLevel(*mf[lev], time, *cmf[0], 0, icomp, ncomp, geom[lev - 1], geom[lev],
436 Integrator::ErrorEst(
int lev, amrex::TagBoxArray& tags, amrex::Real time,
int ngrow)
438 BL_PROFILE(
"Integrator::ErrorEst");
439 TagCellsForRefinement(lev, tags, time, ngrow);
444 Integrator::InitData()
446 BL_PROFILE(
"Integrator::InitData");
448 if (restart_file_cell ==
"" && restart_file_node ==
"")
450 const amrex::Real time = 0.0;
451 InitFromScratch(time);
453 for (
int lev = finest_level - 1; lev >= 0; --lev)
455 if (lev < max_level) regrid(lev, 0.0);
456 for (
int n = 0; n < cell.number_of_fabs; n++)
457 amrex::average_down(*(*cell.fab_array[n])[lev + 1], *(*cell.fab_array[n])[lev],
458 geom[lev + 1], geom[lev],
459 0, (*cell.fab_array[n])[lev]->nComp(), refRatio(lev));
461 SetFinestLevel(finest_level);
463 if (restart_file_cell !=
"")
465 Restart(restart_file_cell,
false);
467 if (restart_file_node !=
"")
469 Restart(restart_file_node,
true);
472 if (plot_int > 0 || plot_dt > 0.0) {
478 Integrator::Restart(
const std::string dirname,
bool a_nodal)
480 BL_PROFILE(
"Integrator::Restart");
482 if (a_nodal && node.fab_array.size() == 0)
487 if (!a_nodal && cell.fab_array.size() == 0)
493 std::string
filename = dirname +
"/Header";
494 std::string chkptfilename = dirname +
"/Checkpoint";
495 amrex::VisMF::IO_Buffer io_buffer(amrex::VisMF::GetIOBufferSize());
496 amrex::Vector<char> fileCharPtr, chkptfileCharPtr;
497 amrex::ParallelDescriptor::ReadAndBcastFile(
filename, fileCharPtr);
498 amrex::ParallelDescriptor::ReadAndBcastFile(chkptfilename, chkptfileCharPtr);
499 std::string fileCharPtrString(fileCharPtr.dataPtr());
500 std::string chkptfileCharPtrString(chkptfileCharPtr.dataPtr());
501 std::istringstream is(fileCharPtrString, std::istringstream::in);
502 std::istringstream chkpt_is(chkptfileCharPtrString, std::istringstream::in);
504 std::string line, word;
507 std::getline(is, line);
512 std::getline(is, line); tmp_numfabs = std::stoi(line);
514 std::vector<std::string> tmp_name_array;
516 for (
int i = 0; i < tmp_numfabs; i++)
518 std::getline(is, line);
519 tmp_name_array.push_back(line);
523 std::getline(is, line);
528 std::getline(is, line); tmp_time = std::stof(line);
Util::Message(
INFO,
"Current time: ", tmp_time);
529 for (
int i = 0; i < max_level + 1; i++)
531 t_new[i] = tmp_time; t_old[i] = tmp_time;
536 std::getline(is, line); tmp_max_level = std::stoi(line);
Util::Message(
INFO,
"Max AMR level: ", line);
537 if (tmp_max_level != max_level)
538 Util::Abort(
INFO,
"The max level specified (", max_level,
") does not match the max level in the restart file (", tmp_max_level,
")");
539 finest_level = tmp_max_level;
551 std::getline(is, line);
553 if (tmp_iters.size() != (
unsigned int)(max_level + 1))
Util::Abort(
INFO,
"Error reading in interation counts: line = ", line);
554 for (
int lev = 0; lev <= max_level; lev++) { istep[lev] = std::stoi(tmp_iters[lev]);
Util::Message(
INFO,
"Iter on level ", lev,
" = ", istep[lev]); }
556 amrex::Vector<amrex::MultiFab> tmpdata(tmp_max_level + 1);
559 if (a_nodal)
for (
unsigned int i = 0; i < node.fab_array.size(); i++) total_ncomp += node.ncomp_array[i];
560 else for (
unsigned int i = 0; i < cell.fab_array.size(); i++) total_ncomp += cell.ncomp_array[i];
562 int total_nghost = a_nodal ? 0 : cell.nghost_array[0];
564 for (
int lev = 0; lev <= max_level; lev++)
566 amrex::BoxArray tmp_ba;
567 tmp_ba.readFrom(chkpt_is);
568 SetBoxArray(lev, tmp_ba);
569 amrex::DistributionMapping tmp_dm(tmp_ba, amrex::ParallelDescriptor::NProcs());
570 SetDistributionMap(lev, tmp_dm);
574 amrex::BoxArray ngrids = grids[lev];
575 ngrids.convert(amrex::IntVect::TheNodeVector());
576 tmpdata[lev].define(ngrids, dmap[lev], total_ncomp, total_nghost);
580 tmpdata[lev].define(grids[lev], dmap[lev], total_ncomp, total_nghost);
582 amrex::VisMF::Read(tmpdata[lev],
583 amrex::MultiFabFileFullPrefix(lev, dirname,
"Level_",
"Cell"));
587 for (
int i = 0; i < node.number_of_fabs; i++)
589 amrex::BoxArray ngrids = grids[lev];
590 ngrids.convert(amrex::IntVect::TheNodeVector());
591 (*node.fab_array[i])[lev].reset(
new amrex::MultiFab(ngrids, dmap[lev], node.ncomp_array[i], node.nghost_array[i]));
592 (*node.fab_array[i])[lev]->setVal(0.);
595 for (
int i = 0; i < cell.number_of_fabs; i++)
596 (*cell.fab_array[i])[lev].reset(
new amrex::MultiFab(grids[lev], dmap[lev], cell.ncomp_array[i], cell.nghost_array[i]));
597 for (
int i = 0; i < tmp_numfabs; i++)
602 for (
int j = 0; j < node.number_of_fabs; j++)
604 if (tmp_name_array[i] == node.name_array[j])
607 Util::Message(
INFO,
"Initializing ", node.name_array[j],
"; nghost=", node.nghost_array[j],
" with ", tmp_name_array[i]);
608 amrex::MultiFab::Copy(*((*node.fab_array[j])[lev]).get(), tmpdata[lev], i, 0, 1, total_nghost);
610 for (
int k = 0; k < node.ncomp_array[j]; k++)
612 if (tmp_name_array[i] == amrex::Concatenate(node.name_array[j], k + 1, 3))
615 Util::Message(
INFO,
"Initializing ", node.name_array[j],
"[", k,
"]; ncomp=", node.ncomp_array[j],
"; nghost=", node.nghost_array[j],
" with ", tmp_name_array[i]);
616 amrex::MultiFab::Copy(*((*node.fab_array[j])[lev]).get(), tmpdata[lev], i, k, 1, total_nghost);
624 for (
int j = 0; j < cell.number_of_fabs; j++)
625 for (
int k = 0; k < cell.ncomp_array[j]; k++)
627 if (tmp_name_array[i] == amrex::Concatenate(cell.name_array[j], k + 1, 3))
630 Util::Message(
INFO,
"Initializing ", cell.name_array[j],
"[", k,
"]; ncomp=", cell.ncomp_array[j],
"; nghost=", cell.nghost_array[j],
" with ", tmp_name_array[i]);
631 amrex::MultiFab::Copy(*((*cell.fab_array[j])[lev]).get(), tmpdata[lev], i, k, 1, cell.nghost_array[j]);
635 if (!match)
Util::Warning(
INFO,
"Fab ", tmp_name_array[i],
" is in the restart file, but there is no fab with that name here.");
638 for (
unsigned int n = 0; n < m_basefields_cell.size(); n++)
640 m_basefields_cell[n]->MakeNewLevelFromScratch(lev, t_new[lev], grids[lev], dmap[lev]);
642 for (
unsigned int n = 0; n < m_basefields.size(); n++)
644 m_basefields[n]->MakeNewLevelFromScratch(lev, t_new[lev], grids[lev], dmap[lev]);
648 SetFinestLevel(max_level);
652 Integrator::MakeNewLevelFromScratch(
int lev, amrex::Real
t,
const amrex::BoxArray& cgrids,
653 const amrex::DistributionMapping& dm)
655 BL_PROFILE(
"Integrator::MakeNewLevelFromScratch");
656 for (
int n = 0; n < cell.number_of_fabs; n++)
658 (*cell.fab_array[n])[lev].reset(
new amrex::MultiFab(cgrids, dm, cell.ncomp_array[n], cell.nghost_array[n]));
659 (*cell.fab_array[n])[lev]->setVal(0.0);
662 amrex::BoxArray ngrids = cgrids;
663 ngrids.convert(amrex::IntVect::TheNodeVector());
664 for (
int n = 0; n < node.number_of_fabs; n++)
666 (*node.fab_array[n])[lev].reset(
new amrex::MultiFab(ngrids, dm, node.ncomp_array[n], node.nghost_array[n]));
667 (*node.fab_array[n])[lev]->setVal(0.0);
669 for (
unsigned int n = 0; n < m_basefields_cell.size(); n++)
671 m_basefields_cell[n]->MakeNewLevelFromScratch(lev,
t, cgrids, dm);
673 for (
unsigned int n = 0; n < m_basefields.size(); n++)
675 m_basefields[n]->MakeNewLevelFromScratch(lev,
t, cgrids, dm);
679 t_old[lev] =
t - dt[lev];
683 for (
int n = 0; n < cell.number_of_fabs; n++)
685 cell.physbc_array[n]->define(geom[lev]);
686 cell.physbc_array[n]->FillBoundary(*(*cell.fab_array[n])[lev], 0, 0,
t, 0);
700 for (
unsigned int n = 0; n < m_basefields_cell.size(); n++)
702 m_basefields_cell[n]->FillBoundary(lev,
t);
704 for (
unsigned int n = 0; n < m_basefields.size(); n++)
706 m_basefields[n]->FillBoundary(lev,
t);
710 std::vector<std::string>
711 Integrator::PlotFileName(
int lev, std::string prefix)
const
713 BL_PROFILE(
"Integrator::PlotFileName");
714 std::vector<std::string> name;
715 name.push_back(plot_file +
"/" + prefix +
"/");
716 name.push_back(amrex::Concatenate(
"", lev, 5));
721 Integrator::WritePlotFile(
bool initial)
const
723 WritePlotFile(t_new[0], istep, initial,
"");
726 Integrator::WritePlotFile(std::string prefix,
Set::Scalar time,
int step)
const
728 amrex::Vector<int> istep(max_level + 1, step);
729 WritePlotFile(time, istep,
false, prefix);
733 Integrator::WritePlotFile(
Set::Scalar time, amrex::Vector<int> iter,
bool initial, std::string prefix)
const
735 BL_PROFILE(
"Integrator::WritePlotFile");
736 int nlevels = finest_level + 1;
737 if (max_plot_level >= 0) nlevels = std::min(nlevels, max_plot_level);
739 int ccomponents = 0, ncomponents = 0, bfcomponents_cell = 0, bfcomponents = 0;
740 amrex::Vector<std::string> cnames, nnames, bfnames_cell, bfnames;
741 for (
int i = 0; i < cell.number_of_fabs; i++)
743 if (!cell.writeout_array[i])
continue;
744 ccomponents += cell.ncomp_array[i];
745 if (cell.ncomp_array[i] > 1)
746 for (
int j = 1; j <= cell.ncomp_array[i]; j++)
747 cnames.push_back(amrex::Concatenate(cell.name_array[i], j, 3));
749 cnames.push_back(cell.name_array[i]);
751 for (
int i = 0; i < node.number_of_fabs; i++)
753 if (!node.writeout_array[i])
continue;
754 ncomponents += node.ncomp_array[i];
755 if (node.ncomp_array[i] > 1)
756 for (
int j = 1; j <= node.ncomp_array[i]; j++)
757 nnames.push_back(amrex::Concatenate(node.name_array[i], j, 3));
759 nnames.push_back(node.name_array[i]);
761 for (
unsigned int i = 0; i < m_basefields_cell.size(); i++)
763 if (m_basefields_cell[i]->writeout)
765 bfcomponents_cell += m_basefields_cell[i]->NComp();
766 for (
int j = 0; j < m_basefields_cell[i]->NComp(); j++)
767 bfnames_cell.push_back(m_basefields_cell[i]->Name(j));
770 for (
unsigned int i = 0; i < m_basefields.size(); i++)
772 if (m_basefields[i]->writeout)
774 bfcomponents += m_basefields[i]->NComp();
775 for (
int j = 0; j < m_basefields[i]->NComp(); j++)
776 bfnames.push_back(m_basefields[i]->Name(j));
780 amrex::Vector<amrex::MultiFab> cplotmf(nlevels), nplotmf(nlevels);
782 bool do_cell_plotfile = (ccomponents + bfcomponents_cell > 0 || (ncomponents + bfcomponents > 0 && cell.all)) && cell.any;
783 bool do_node_plotfile = (ncomponents + bfcomponents > 0 || (ccomponents + bfcomponents_cell > 0 && node.all)) && node.any;
785 for (
int ilev = 0; ilev < nlevels; ++ilev)
787 if (do_cell_plotfile)
789 int ncomp = ccomponents + bfcomponents_cell;
790 if (cell.all) ncomp += ncomponents + bfcomponents;
791 cplotmf[ilev].define(grids[ilev], dmap[ilev], ncomp, 0);
794 for (
int i = 0; i < cell.number_of_fabs; i++)
796 if (!cell.writeout_array[i])
continue;
797 if ((*cell.fab_array[i])[ilev]->contains_nan())
799 if (abort_on_nan)
Util::Abort(
INFO, cnames[i],
" contains nan (i=", i,
")");
802 if ((*cell.fab_array[i])[ilev]->contains_inf())
804 if (abort_on_nan)
Util::Abort(
INFO, cnames[i],
" contains inf (i=", i,
")");
807 amrex::MultiFab::Copy(cplotmf[ilev], *(*cell.fab_array[i])[ilev], 0, n, cell.ncomp_array[i], 0);
808 n += cell.ncomp_array[i];
810 for (
unsigned int i = 0; i < m_basefields_cell.size(); i++)
812 if (m_basefields_cell[i]->writeout)
814 m_basefields_cell[i]->Copy(ilev, cplotmf[ilev], n, 0);
815 n += m_basefields_cell[i]->NComp();
821 for (
int i = 0; i < node.number_of_fabs; i++)
823 if (!node.writeout_array[i])
continue;
824 if ((*node.fab_array[i])[ilev]->contains_nan())
Util::Abort(
INFO, nnames[i],
" contains nan (i=", i,
")");
825 if ((*node.fab_array[i])[ilev]->contains_inf())
Util::Abort(
INFO, nnames[i],
" contains inf (i=", i,
")");
826 amrex::average_node_to_cellcenter(cplotmf[ilev], n, *(*node.fab_array[i])[ilev], 0, node.ncomp_array[i], 0);
827 n += node.ncomp_array[i];
830 if (bfcomponents > 0)
832 amrex::BoxArray ngrids = grids[ilev];
833 ngrids.convert(amrex::IntVect::TheNodeVector());
834 amrex::MultiFab bfplotmf(ngrids, dmap[ilev], bfcomponents, 0);
836 for (
unsigned int i = 0; i < m_basefields.size(); i++)
838 if (m_basefields[i]->writeout)
840 m_basefields[i]->Copy(ilev, bfplotmf, ctr, 0);
841 ctr += m_basefields[i]->NComp();
844 amrex::average_node_to_cellcenter(cplotmf[ilev], n, bfplotmf, 0, bfcomponents);
850 if (do_node_plotfile)
852 amrex::BoxArray ngrids = grids[ilev];
853 ngrids.convert(amrex::IntVect::TheNodeVector());
854 int ncomp = ncomponents + bfcomponents;
855 if (node.all) ncomp += ccomponents + bfcomponents_cell;
856 nplotmf[ilev].define(ngrids, dmap[ilev], ncomp, 0);
859 for (
int i = 0; i < node.number_of_fabs; i++)
861 if (!node.writeout_array[i])
continue;
862 if ((*node.fab_array[i])[ilev]->contains_nan())
Util::Warning(
INFO, nnames[i],
" contains nan (i=", i,
"). Resetting to zero.");
863 if ((*node.fab_array[i])[ilev]->contains_inf())
Util::Abort(
INFO, nnames[i],
" contains inf (i=", i,
")");
864 amrex::MultiFab::Copy(nplotmf[ilev], *(*node.fab_array[i])[ilev], 0, n, node.ncomp_array[i], 0);
865 n += node.ncomp_array[i];
867 for (
unsigned int i = 0; i < m_basefields.size(); i++)
869 if (m_basefields[i]->writeout)
871 m_basefields[i]->Copy(ilev, nplotmf[ilev], n, 0);
872 n += m_basefields[i]->NComp();
878 for (
int i = 0; i < cell.number_of_fabs; i++)
880 if (!cell.writeout_array[i])
continue;
881 if ((*cell.fab_array[i])[ilev]->contains_nan())
Util::Warning(
INFO, cnames[i],
" contains nan (i=", i,
"). Resetting to zero.");
882 if ((*cell.fab_array[i])[ilev]->contains_inf())
Util::Abort(
INFO, cnames[i],
" contains inf (i=", i,
")");
883 if ((*cell.fab_array[i])[ilev]->nGrow() == 0)
885 if (initial)
Util::Warning(
INFO, cnames[i],
" has no ghost cells and will not be included in nodal output");
889 n += cell.ncomp_array[i];
892 if (bfcomponents_cell > 0)
894 amrex::BoxArray cgrids = grids[ilev];
895 amrex::MultiFab bfplotmf(cgrids, dmap[ilev], bfcomponents_cell, 0);
897 for (
unsigned int i = 0; i < m_basefields_cell.size(); i++)
899 if (m_basefields_cell[i]->writeout)
901 m_basefields_cell[i]->Copy(ilev, bfplotmf, ctr, 0);
902 ctr += m_basefields_cell[i]->NComp();
906 n += bfcomponents_cell;
912 std::vector<std::string> plotfilename = PlotFileName(istep[0], prefix);
913 if (initial) plotfilename[1] = plotfilename[1] +
"init";
915 if (do_cell_plotfile)
917 amrex::Vector<std::string> allnames = cnames;
918 allnames.insert(allnames.end(), bfnames_cell.begin(), bfnames_cell.end());
920 allnames.insert(allnames.end(), nnames.begin(), nnames.end());
921 allnames.insert(allnames.end(), bfnames.begin(), bfnames.end());
923 WriteMultiLevelPlotfile(plotfilename[0] + plotfilename[1] +
"cell", nlevels, amrex::GetVecOfConstPtrs(cplotmf), allnames,
924 Geom(), time, iter, refRatio());
926 std::ofstream chkptfile;
927 chkptfile.open(plotfilename[0] + plotfilename[1] +
"cell/Checkpoint");
928 for (
int i = 0; i <= max_level; i++) boxArray(i).writeOn(chkptfile);
932 if (do_node_plotfile)
934 amrex::Vector<std::string> allnames = nnames;
935 allnames.insert(allnames.end(), bfnames.begin(), bfnames.end());
936 if (node.all) allnames.insert(allnames.end(), cnames.begin(), cnames.end());
937 WriteMultiLevelPlotfile(plotfilename[0] + plotfilename[1] +
"node", nlevels, amrex::GetVecOfConstPtrs(nplotmf), allnames,
938 Geom(), time, iter, refRatio());
940 std::ofstream chkptfile;
941 chkptfile.open(plotfilename[0] + plotfilename[1] +
"node/Checkpoint");
942 for (
int i = 0; i <= max_level; i++) boxArray(i).writeOn(chkptfile);
946 if (amrex::ParallelDescriptor::IOProcessor())
948 std::ofstream coutfile, noutfile;
951 if (do_cell_plotfile) coutfile.open(plot_file +
"/celloutput.visit", std::ios_base::out);
952 if (do_node_plotfile) noutfile.open(plot_file +
"/nodeoutput.visit", std::ios_base::out);
956 if (do_cell_plotfile) coutfile.open(plot_file +
"/celloutput.visit", std::ios_base::app);
957 if (do_node_plotfile) noutfile.open(plot_file +
"/nodeoutput.visit", std::ios_base::app);
959 if (do_cell_plotfile) coutfile << plotfilename[1] +
"cell" +
"/Header" << std::endl;
960 if (do_node_plotfile) noutfile << plotfilename[1] +
"node" +
"/Header" << std::endl;
967 BL_PROFILE(
"Integrator::Evolve");
968 amrex::Real cur_time = t_new[0];
969 int last_plot_file_step = 0;
971 for (
int step = istep[0]; step < max_step && cur_time < stop_time; ++step)
973 if (amrex::ParallelDescriptor::IOProcessor()) {
974 std::cout <<
"\nSTEP " << step + 1 <<
" starts ..." << std::endl;
978 TimeStepBegin(cur_time, step);
979 if (integrate_variables_before_advance) IntegrateVariables(cur_time, step);
980 TimeStep(lev, cur_time, iteration);
981 if (integrate_variables_after_advance) IntegrateVariables(cur_time, step);
982 TimeStepComplete(cur_time, step);
985 if (amrex::ParallelDescriptor::IOProcessor()) {
986 std::cout <<
"STEP " << step + 1 <<
" ends."
987 <<
" TIME = " << cur_time <<
" DT = " << dt[0]
992 for (
int lev = 0; lev <= finest_level; ++lev) {
993 t_new[lev] = cur_time;
996 if (plot_int > 0 && (step + 1) % plot_int == 0) {
997 last_plot_file_step = step + 1;
1001 else if (std::fabs(std::remainder(cur_time, plot_dt)) < 0.5 * dt[0])
1003 last_plot_file_step = step + 1;
1008 if (cur_time >= stop_time - 1.e-6 * dt[0])
break;
1010 if (plot_int > 0 && istep[0] > last_plot_file_step) {
1016 Integrator::IntegrateVariables(amrex::Real time,
int step)
1018 BL_PROFILE(
"Integrator::IntegrateVariables");
1019 if (!thermo.number)
return;
1021 if ((thermo.interval > 0 && (step) % thermo.interval == 0) ||
1022 ((thermo.dt > 0.0) && (std::fabs(std::remainder(time, plot_dt)) < 0.5 * dt[0])))
1025 for (
int i = 0; i < thermo.number; i++)
1027 if (thermo.extensives[i]) *thermo.vars[i] = 0;
1031 for (
int ilev = 0; ilev < max_level; ilev++)
1033 const amrex::BoxArray& cfba = amrex::coarsen(grids[ilev + 1], refRatio(ilev));
1036 #pragma omp parallel
1038 for (amrex::MFIter mfi(grids[ilev], dmap[ilev],
true); mfi.isValid(); ++mfi)
1040 const amrex::Box& box = mfi.tilebox();
1041 const amrex::BoxArray& comp = amrex::complementIn(box, cfba);
1043 for (
int i = 0; i < comp.size(); i++)
1045 Integrate(ilev, time, step,
1053 #pragma omp parallel
1055 for (amrex::MFIter mfi(grids[max_level], dmap[max_level],
true); mfi.isValid(); ++mfi)
1057 const amrex::Box& box = mfi.tilebox();
1058 Integrate(max_level, time, step, mfi, box);
1063 for (
int i = 0; i < thermo.number; i++)
1065 if (thermo.extensives[i])
1066 amrex::ParallelDescriptor::ReduceRealSum(*thermo.vars[i]);
1069 if (amrex::ParallelDescriptor::IOProcessor() &&
1071 (thermo.plot_int > 0 && step % thermo.plot_int == 0) ||
1072 (thermo.plot_dt > 0.0 && std::fabs(std::remainder(time, thermo.plot_dt)) < 0.5 * dt[0])
1075 std::ofstream outfile;
1078 outfile.open(plot_file +
"/thermo.dat", std::ios_base::out);
1080 for (
int i = 0; i < thermo.number; i++)
1081 outfile <<
"\t" << thermo.names[i];
1082 outfile << std::endl;
1084 else outfile.open(plot_file +
"/thermo.dat", std::ios_base::app);
1086 for (
int i = 0; i < thermo.number; i++)
1087 outfile <<
"\t" << *thermo.vars[i];
1088 outfile << std::endl;
1096 Integrator::TimeStep(
int lev, amrex::Real time,
int )
1098 BL_PROFILE(
"Integrator::TimeStep");
1099 if (base_regrid_int <= 0 || istep[0] % base_regrid_int == 0)
1101 if (regrid_int > 0 || base_regrid_int > 0)
1103 static amrex::Vector<int> last_regrid_step(max_level + 1, 0);
1106 if (lev < max_level && istep[lev] > last_regrid_step[lev])
1108 if (istep[lev] % regrid_int == 0)
1110 regrid(lev, time,
false);
1115 SetFinestLevel(finest_level);
1117 if (Verbose() && amrex::ParallelDescriptor::IOProcessor()) {
1118 std::cout <<
"[Level " << lev
1119 <<
" step " << istep[lev] + 1 <<
"] ";
1120 std::cout <<
"ADVANCE with dt = "
1125 for (
int n = 0; n < cell.number_of_fabs; n++)
1126 FillPatch(lev, time, *cell.fab_array[n], *(*cell.fab_array[n])[lev], *cell.physbc_array[n], 0);
1127 for (
int n = 0; n < node.number_of_fabs; n++)
1128 FillPatch(lev, time, *node.fab_array[n], *(*node.fab_array[n])[lev], *node.physbc_array[n], 0);
1129 for (
unsigned int n = 0; n < m_basefields_cell.size(); n++)
1130 if (m_basefields_cell[n]->evolving) m_basefields_cell[n]->FillPatch(lev, time);
1131 for (
unsigned int n = 0; n < m_basefields.size(); n++)
1132 if (m_basefields[n]->evolving) m_basefields[n]->FillPatch(lev, time);
1134 Advance(lev, time, dt[lev]);
1137 if (Verbose() && amrex::ParallelDescriptor::IOProcessor())
1139 std::cout <<
"[Level " << lev
1140 <<
" step " << istep[lev] <<
"] ";
1141 std::cout <<
"Advanced "
1147 if (lev < finest_level)
1149 for (
int i = 1; i <= nsubsteps[lev + 1]; ++i)
1150 TimeStep(lev + 1, time + (i - 1) * dt[lev + 1], i);
1152 for (
int n = 0; n < cell.number_of_fabs; n++)
1154 amrex::average_down(*(*cell.fab_array[n])[lev + 1], *(*cell.fab_array[n])[lev],
1155 geom[lev + 1], geom[lev],
1156 0, (*cell.fab_array[n])[lev]->nComp(), refRatio(lev));
1158 for (
int n = 0; n < node.number_of_fabs; n++)
1160 amrex::average_down(*(*node.fab_array[n])[lev + 1], *(*node.fab_array[n])[lev],
1161 0, (*node.fab_array[n])[lev]->nComp(), refRatio(lev));
1163 for (
unsigned int n = 0; n < m_basefields_cell.size(); n++)
1165 if (m_basefields_cell[n]->evolving)
1166 m_basefields_cell[n]->AverageDown(lev, refRatio(lev));
1168 for (
unsigned int n = 0; n < m_basefields.size(); n++)
1170 if (m_basefields[n]->evolving)
1171 m_basefields[n]->AverageDown(lev, refRatio(lev));