Alamo
Integrator.cpp
Go to the documentation of this file.
1///
2/// \file Integrator.cpp
3/// \brief Compute the volume integral of two multiplied Fourier series
4///
5
6#include "Integrator.H"
7#include "IO/FileNameParse.H"
8#include "IO/ParmParse.H"
9#include "Util/Util.H"
10#include "Unit/Unit.H"
11#include <numeric>
12
13
14
15namespace Integrator
16{
17
18Integrator::Integrator() : amrex::AmrCore()
19{
20 IO::ParmParse pp;
21 Integrator::Parse (*this, pp);
22}
23void
24Integrator::Parse(Integrator &value, IO::ParmParse &pp)
25{
26 BL_PROFILE("Integrator::Integrator()");
27 {
28 // These are basic parameters that are, in
29 // general, common to all Alamo simulations.
30 // Number of iterations before ending (default is maximum possible int)
31 pp.query_default("max_step", value.max_step, 2147483647);
32 // Simulation time before ending
33 pp.query_required("stop_time", value.stop_time, Unit::Time());
34 // Nominal timestep on amrlev = 0
35 pp.query_required("timestep", value.timestep, Unit::Time());
36 pp.query_default("restart", value.restart_file_cell,""); // Name of restart file to read from
37 pp.query_default("restart_cell", value.restart_file_cell,""); // Name of cell-fab restart file to read from
38 pp.query_default("restart_node", value.restart_file_node,""); // Name of node-fab restart file to read from
39 }
40#ifdef AMREX_USE_HDF5
41 {
42 // These parameters are only used when Alamo
43 // is compiled with HDF5 support.
44 // The level of compression to use with HDF5's zlib compression algorithm (0-9)
45 pp.query_default("compression_level", value.compression_level, "6");
46 }
47#endif
48 {
49 // These are parameters that are specific to
50 // the AMR/regridding part of the code.
51 pp.query_default("amr.regrid_int", value.regrid_int, 2); // Regridding interval in step numbers
52 pp.query_default("amr.base_regrid_int", value.base_regrid_int, 0); // Regridding interval based on coarse level only
53 pp.query_default("amr.plot_int", value.plot_int, -1); // Interval (in timesteps) between plotfiles (Default negative value will cause the plot interval to be ignored.)
54 pp.query_default("amr.plot_dt", value.plot_dt, "-1.0", Unit::Time()); // Interval (in simulation time) between plotfiles (Default negative value will cause the plot dt to be ignored.)
55
56
57 // Output file: see IO::FileNameParse for wildcards and variable substitution
58 pp.query_default("amr.plot_file", value.plot_file, "output");
59
60 pp.query_default("amr.cell.all", value.cell.all, false); // Turn on to write all output in cell fabs (default: off)
61 pp.query_default("amr.cell.any", value.cell.any, true); // Turn off to prevent any cell based output (default: on)
62 pp.query_default("amr.node.all", value.node.all, false); // Turn on to write all output in node fabs (default: off)
63 pp.query_default("amr.node.any", value.node.any, true); // Turn off to prevent any node based output (default: on)
64
65 pp.query_default("amr.abort_on_nan",value.abort_on_nan, true); // Abort if a plotfile contains nan or inf.
66
67 Util::Assert(INFO, TEST(!(!value.cell.any && value.cell.all)));
68 Util::Assert(INFO, TEST(!(!value.node.any && value.node.all)));
69
70 pp.query_default("amr.max_plot_level", value.max_plot_level, -1); // Specify a maximum level of refinement for output files (NO REFINEMENT)
71
72 pp.query_default("amr.print_ghost_nodes", value.print_ghost_nodes, 0); // include ghost nodes in output
73 pp.query_default("amr.print_ghost_cells", value.print_ghost_cells, 0); // include ghost cells in output
74
75 IO::FileNameParse(value.plot_file);
76
77 value.nsubsteps.resize(value.maxLevel() + 1, 1);
78 int cnt = pp.countval("amr.nsubsteps");
79 if (cnt != 0)
80 if (cnt == value.maxLevel()) {
81 pp.queryarr("amr.nsubsteps", value.nsubsteps); // Number of substeps to take on each level (default: 2)
82 value.nsubsteps.insert(value.nsubsteps.begin(), 1);
83 value.nsubsteps.pop_back();
84 }
85 else if (cnt == 1)
86 {
87 int nsubsteps_all;
88 pp.query_required("amr.nsubsteps", nsubsteps_all);// Number of substeps to take on each level (set all levels to this value)
89 for (int lev = 1; lev <= value.maxLevel(); ++lev) value.nsubsteps[lev] = nsubsteps_all;
90 }
91 else
92 Util::Abort(INFO, "number of nsubsteps input must equal either 1 or amr.max_level");
93 else
94 for (int lev = 1; lev <= value.maxLevel(); ++lev)
95 value.nsubsteps[lev] = value.MaxRefRatio(lev - 1);
96 }
97 {
98 // activate dynamic CFL-based timestep
99 pp.query_default("dynamictimestep.on",value.dynamictimestep.on,false);
100 if (value.dynamictimestep.on)
101 {
102 // how much information to print
103 pp.query_validate("dynamictimestep.verbose",value.dynamictimestep.verbose,{0,1});
104 // number of previous timesteps for rolling average
105 pp.query_default("dynamictimestep.nprevious",value.dynamictimestep.nprevious,5);
106 // dynamic teimstep CFL condition
107 pp.query_default("dynamictimestep.cfl",value.dynamictimestep.cfl,1.0);
108 // minimum timestep size allowed shen stepping dynamically
109 pp.query_default("dynamictimestep.min",value.dynamictimestep.min,value.timestep);
110 // maximum timestep size allowed shen stepping dynamically
111 pp.query_default("dynamictimestep.max",value.dynamictimestep.max,value.timestep);
112
113 Util::AssertException(INFO,TEST(value.dynamictimestep.max >= value.dynamictimestep.min));
114 }
115 }
116 {
117 // Information on how to generate thermodynamic
118 // data (to show up in thermo.dat)
119 value.thermo.interval = 1; // Default: integrate every time.
120 pp.query_default("amr.thermo.int", value.thermo.interval, 1); // Integration interval (1)
121 pp.query_default("amr.thermo.plot_int", value.thermo.plot_int, -1); // Interval (in timesteps) between writing (Default negative value will cause the plot interval to be ignored.)
122 pp.query_default("amr.thermo.plot_dt", value.thermo.plot_dt, "-1.0", Unit::Time()); // Interval (in simulation time) between writing (Default negative value will cause the plot dt to be ignored.)
123 }
124
125 {
126 // Instead of using AMR, prescribe an explicit, user-defined
127 // set of grids to work on. This is pretty much always used
128 // for testing purposes only.
129 pp.query_default("explicitmesh.on", value.explicitmesh.on, 0); // Use explicit mesh instead of AMR
130 if (value.explicitmesh.on)
131 {
132 for (int ilev = 0; ilev < value.maxLevel(); ++ilev)
133 {
134 std::string strlo = "explicitmesh.lo" + std::to_string(ilev + 1);
135 std::string strhi = "explicitmesh.hi" + std::to_string(ilev + 1);
136
137 Util::Assert(INFO, TEST(pp.contains(strlo.c_str())));
138 Util::Assert(INFO, TEST(pp.contains(strhi.c_str())));
139
140 amrex::Vector<int> lodata, hidata;
141 pp.queryarr(strlo.c_str(), lodata);
142 pp.queryarr(strhi.c_str(), hidata);
143 amrex::IntVect lo(AMREX_D_DECL(lodata[0], lodata[1], lodata[2]));
144 amrex::IntVect hi(AMREX_D_DECL(hidata[0], hidata[1], hidata[2]));
145
146 value.explicitmesh.box.push_back(amrex::Box(lo, hi));
147 }
148 }
149 }
150
151
152 {
153 //
154 // These parameters are NOT USED!
155 // They are shadow paramters for AMReX::TimeIntegration.
156 // The purpose here is:
157 // (1) to set default values
158 // (2) to prevent "unused input" warnings when integration is parsed, since it is parsed in the
159 // AMReX convention, not the Alamo convention.
160 //
161
162 std::string str;
163 // Type of time integration to use (see amrex::TimeIntegrator for more details)
164 pp.query_validate("integration.type", str, {"ForwardEuler","RungeKutta"});
165 if (str == "RungeKutta")
166 {
167 int type;
168 // If RungeKutta specified, which order to use (3=SSPRK3, 4=RK4)
169 pp.query_validate("integration.rk.type", type, {1,2,3,4});
170 }
171 }
172
173 int nlevs_max = value.maxLevel() + 1;
174
175 value.istep.resize(nlevs_max, 0);
176
177 value.t_new.resize(nlevs_max, 0.0);
178 value.t_old.resize(nlevs_max, -1.e100);
179 value.SetTimestep(value.timestep);
180
181 value.plot_file = Util::GetFileName();
182 IO::WriteMetaData(value.plot_file, IO::Status::Running, 0);
183}
184
185// Destructor
186Integrator::~Integrator()
187{
188 if (Util::finalized)
189 {
190 std::cout << "!! ERROR !! Integrator destructor called after alamo has been finalized." << std::endl;
191 std::cout << " Behavior occurring after this is undefined." << std::endl;
192 std::abort();
193 }
194
195 // Close out the metadata file and mark completed.
196 IO::WriteMetaData(plot_file, IO::Status::Complete);
197
198 // De-initialize all of the base fields and clear the arrays.
199 for (unsigned int i = 0; i < m_basefields.size(); i++) delete m_basefields[i];
200 for (unsigned int i = 0; i < m_basefields_cell.size(); i++) delete m_basefields_cell[i];
201 m_basefields.clear();
202 m_basefields_cell.clear();
203}
204
205void Integrator::SetTimestep(Set::Scalar _timestep)
206{
207 BL_PROFILE("Integrator::SetTimestep");
208 int nlevs_max = maxLevel() + 1;
209 timestep = _timestep;
210 dt.resize(nlevs_max, 1.e100);
211 dt[0] = timestep;
212 for (int i = 1; i < nlevs_max; i++)
213 dt[i] = dt[i - 1] / (amrex::Real)nsubsteps[i];
214}
215void Integrator::SetPlotInt(int a_plot_int)
216{
217 BL_PROFILE("Integrator::SetPlotInt");
218 plot_int = a_plot_int;
219}
220
221/// \fn Integrator::MakeNewLevelFromCoarse
222/// \brief Wrapper to call FillCoarsePatch
223/// \note **THIS OVERRIDES A PURE VIRTUAL METHOD - DO NOT CHANGE**
224///
225void
226Integrator::MakeNewLevelFromCoarse(int lev, amrex::Real time, const amrex::BoxArray& cgrids, const amrex::DistributionMapping& dm)
227{
228 BL_PROFILE("Integrator::MakeNewLevelFromCoarse");
229
230 for (int n = 0; n < cell.number_of_fabs; n++)
231 {
232 const int ncomp = (*cell.fab_array[n])[lev - 1]->nComp();
233 const int nghost = (*cell.fab_array[n])[lev - 1]->nGrow();
234
235 (*cell.fab_array[n])[lev].reset(new amrex::MultiFab(cgrids, dm, ncomp, nghost));
236
237 (*cell.fab_array[n])[lev]->setVal(0.0);
238
239 FillCoarsePatch(lev, time, *cell.fab_array[n], *cell.physbc_array[n], 0, ncomp);
240 }
241
242 amrex::BoxArray ngrids = cgrids;
243 ngrids.convert(amrex::IntVect::TheNodeVector());
244
245 for (int n = 0; n < node.number_of_fabs; n++)
246 {
247 const int ncomp = (*node.fab_array[n])[lev - 1]->nComp();
248 const int nghost = (*node.fab_array[n])[lev - 1]->nGrow();
249
250 (*node.fab_array[n])[lev].reset(new amrex::MultiFab(ngrids, dm, ncomp, nghost));
251 (*node.fab_array[n])[lev]->setVal(0.0);
252
253 FillCoarsePatch(lev, time, *node.fab_array[n], *node.physbc_array[n], 0, ncomp);
254 }
255
256 for (unsigned int n = 0; n < m_basefields.size(); n++)
257 {
258 m_basefields[n]->MakeNewLevelFromCoarse(lev, time, cgrids, dm);
259 }
260 for (unsigned int n = 0; n < m_basefields_cell.size(); n++)
261 {
262 m_basefields_cell[n]->MakeNewLevelFromCoarse(lev, time, cgrids, dm);
263 }
264
265 Regrid(lev, time);
266}
267
268
269///
270/// RESETS ALL MULTIFABS AT A GIVEN LEVEL
271///
272/// (OVERRIDES PURE VIRTUAL METHOD - DO NOT CHANGE)
273///
274void
275Integrator::RemakeLevel(int lev, ///<[in] AMR Level
276 amrex::Real time, ///<[in] Simulation time
277 const amrex::BoxArray &cgrids, ///<[in] Coarse grids
278 const amrex::DistributionMapping &dm ///[in] Distribution mapping
279 )
280{
281 BL_PROFILE("Integrator::RemakeLevel");
282 for (int n = 0; n < cell.number_of_fabs; n++)
283 {
284 const int ncomp = (*cell.fab_array[n])[lev]->nComp();
285 const int nghost = (*cell.fab_array[n])[lev]->nGrow();
286
287 amrex::MultiFab new_state(cgrids, dm, ncomp, nghost);
288
289 new_state.setVal(0.0);
290 FillPatch(lev, time, *cell.fab_array[n], new_state, *cell.physbc_array[n], 0);
291 std::swap(new_state, *(*cell.fab_array[n])[lev]);
292 }
293
294 amrex::BoxArray ngrids = cgrids;
295 ngrids.convert(amrex::IntVect::TheNodeVector());
296
297 for (int n = 0; n < node.number_of_fabs; n++)
298 {
299 const int ncomp = (*node.fab_array[n])[lev]->nComp();
300 const int nghost = (*node.fab_array[n])[lev]->nGrow();
301
302 amrex::MultiFab new_state(ngrids, dm, ncomp, nghost);
303
304 new_state.setVal(0.0);
305 FillPatch(lev, time, *node.fab_array[n], new_state, *node.physbc_array[n], 0);
306 std::swap(new_state, *(*node.fab_array[n])[lev]);
307 }
308
309 for (unsigned int n = 0; n < m_basefields_cell.size(); n++)
310 {
311 m_basefields_cell[n]->RemakeLevel(lev, time, cgrids, dm);
312 }
313 for (unsigned int n = 0; n < m_basefields.size(); n++)
314 {
315 m_basefields[n]->RemakeLevel(lev, time, cgrids, dm);
316 }
317 Regrid(lev, time);
318}
319
320//
321// DELETE EVERYTHING
322//
323// (OVERRIDES PURE VIRTUAL METHOD - DO NOT CHANGE)
324//
325void
326Integrator::ClearLevel(int lev)
327{
328 BL_PROFILE("Integrator::ClearLevel");
329 for (int n = 0; n < cell.number_of_fabs; n++)
330 {
331 (*cell.fab_array[n])[lev].reset(nullptr);
332 }
333 for (int n = 0; n < node.number_of_fabs; n++)
334 {
335 (*node.fab_array[n])[lev].reset(nullptr);
336 }
337}
338
339//
340//
341//
342
343
344
345
346
347void
348Integrator::RegisterNewFab(Set::Field<Set::Scalar>& new_fab, BC::BC<Set::Scalar>* new_bc, int ncomp, int nghost, std::string name, bool writeout, bool evolving, std::vector<std::string> suffix)
349{
350 //Util::Warning(INFO, "RegisterNewFab is depricated. Please replace with AddField");
351 AddField<Set::Scalar, Set::Hypercube::Cell>(new_fab, new_bc, ncomp, nghost, name, writeout, evolving, suffix);
352}
353void
354Integrator::RegisterNewFab(Set::Field<Set::Scalar>& new_fab, int ncomp, std::string name, bool writeout, bool evolving, std::vector<std::string> suffix)
355{
356 //Util::Warning(INFO, "RegisterNewFab is depricated. Please replace with AddField");
357 AddField<Set::Scalar, Set::Hypercube::Cell>(new_fab, nullptr, ncomp, 0, name, writeout, evolving, suffix);
358}
359void
360Integrator::RegisterNodalFab(Set::Field<Set::Scalar>& new_fab, BC::BC<Set::Scalar>* new_bc, int ncomp, int nghost, std::string name, bool writeout, bool evolving, std::vector<std::string> suffix)
361{
362 //Util::Warning(INFO, "RegisterNodalFab is depricated. Please replace with AddField");
363 AddField<Set::Scalar, Set::Hypercube::Node>(new_fab, new_bc, ncomp, nghost, name, writeout, evolving,suffix);
364}
365void
366Integrator::RegisterNodalFab(Set::Field<Set::Scalar>& new_fab, int ncomp, int nghost, std::string name, bool writeout, bool evolving, std::vector<std::string> suffix)
367{
368 //Util::Warning(INFO, "RegisterNodalFab is depricated. Please replace with AddField");
369 AddField<Set::Scalar, Set::Hypercube::Node>(new_fab, nullptr, ncomp, nghost, name, writeout, evolving,suffix);
370}
371
372
373
374
375void // CUSTOM METHOD - CHANGEABLE
376Integrator::RegisterIntegratedVariable(Set::Scalar *integrated_variable, std::string name, bool extensive)
377{
378 BL_PROFILE("Integrator::RegisterIntegratedVariable");
379 thermo.vars.push_back(integrated_variable);
380 thermo.names.push_back(name);
381 thermo.extensives.push_back(extensive);
382 thermo.number++;
383}
384
385long // CUSTOM METHOD - CHANGEABLE
386Integrator::CountCells(int lev)
387{
388 BL_PROFILE("Integrator::CountCells");
389 const int N = grids[lev].size();
390
391 long cnt = 0;
392
393 for (int i = 0; i < N; ++i)
394 {
395 cnt += grids[lev][i].numPts();
396 }
397
398 return cnt;
399}
400
401void // CUSTOM METHOD - CHANGEABLE
402Integrator::FillPatch(int lev, amrex::Real time,
403 amrex::Vector<std::unique_ptr<amrex::MultiFab>>& source_mf,
404 amrex::MultiFab& destination_mf,
405 BC::BC<Set::Scalar>& physbc, int icomp)
406{
407 BL_PROFILE("Integrator::FillPatch");
408 if (lev == 0)
409 {
410
411 amrex::Vector<amrex::MultiFab*> smf;
412 smf.push_back(source_mf[lev].get());
413 amrex::Vector<amrex::Real> stime;
414 stime.push_back(time);
415
416 physbc.define(geom[lev]);
417 amrex::FillPatchSingleLevel(destination_mf, // Multifab
418 time, // time
419 smf, // Vector<MultiFab*> &smf (CONST)
420 stime, // Vector<Real> &stime (CONST)
421 0, // scomp - Source component
422 icomp, // dcomp - Destination component
423 destination_mf.nComp(), // ncomp - Number of components
424 geom[lev], // Geometry (CONST)
425 physbc,
426 0); // BC
427 }
428 else
429 {
430 amrex::Vector<amrex::MultiFab*> cmf, fmf;
431 cmf.push_back(source_mf[lev - 1].get());
432 fmf.push_back(source_mf[lev].get());
433 amrex::Vector<amrex::Real> ctime, ftime;
434 ctime.push_back(time);
435 ftime.push_back(time);
436
437 physbc.define(geom[lev]);
438
439 amrex::Interpolater* mapper;
440
441 if (destination_mf.boxArray().ixType() == amrex::IndexType::TheNodeType())
442 mapper = &amrex::node_bilinear_interp;
443 else
444 mapper = &amrex::cell_cons_interp;
445
446 amrex::Vector<amrex::BCRec> bcs(destination_mf.nComp(), physbc.GetBCRec()); // todo
447 amrex::FillPatchTwoLevels(destination_mf, time, cmf, ctime, fmf, ftime,
448 0, icomp, destination_mf.nComp(), geom[lev - 1], geom[lev],
449 physbc, 0,
450 physbc, 0,
451 refRatio(lev - 1),
452 mapper, bcs, 0);
453 }
454}
455
456/// \fn Integrator::FillCoarsePatch
457/// \brief Fill a fab at current level with the data from one level up
458///
459/// \note This is a custom method and is changeable
460void
461Integrator::FillCoarsePatch(int lev, ///<[in] AMR level
462 amrex::Real time, ///<[in] Simulatinon time
463 Set::Field<Set::Scalar>& mf, ///<[in] Fab to fill
464 BC::BC<Set::Scalar>& physbc, ///<[in] BC object applying to Fab
465 int icomp, ///<[in] start component
466 int ncomp) ///<[in] end component (i.e. applies to components `icomp`...`ncomp`)
467{
468 BL_PROFILE("Integrator::FillCoarsePatch");
469 AMREX_ASSERT(lev > 0);
470 amrex::Vector<amrex::MultiFab*> cmf;
471 cmf.push_back(mf[lev - 1].get());
472 amrex::Vector<amrex::Real> ctime;
473 ctime.push_back(time);
474
475 physbc.define(geom[lev]);
476
477 amrex::Interpolater* mapper;
478 if (mf[lev]->boxArray().ixType() == amrex::IndexType::TheNodeType())
479 mapper = &amrex::node_bilinear_interp;
480 else
481 mapper = &amrex::cell_cons_interp;
482
483 amrex::Vector<amrex::BCRec> bcs(ncomp, physbc.GetBCRec());
484 amrex::InterpFromCoarseLevel(*mf[lev], time, *cmf[0], 0, icomp, ncomp, geom[lev - 1], geom[lev],
485 physbc, 0,
486 physbc, 0,
487 refRatio(lev - 1),
488 mapper, bcs, 0);
489}
490
491void
492Integrator::ErrorEst(int lev, amrex::TagBoxArray& tags, amrex::Real time, int ngrow)
493{
494 BL_PROFILE("Integrator::ErrorEst");
495 TagCellsForRefinement(lev, tags, time, ngrow);
496}
497
498
499void
500Integrator::InitData()
501{
502 BL_PROFILE("Integrator::InitData");
503
504 if (restart_file_cell == "" && restart_file_node == "")
505 {
506 const amrex::Real time = 0.0;
507 InitFromScratch(time);
508
509 for (int lev = finest_level - 1; lev >= 0; --lev)
510 {
511 if (lev < max_level) regrid(lev, 0.0);
512 for (int n = 0; n < cell.number_of_fabs; n++)
513 amrex::average_down(*(*cell.fab_array[n])[lev + 1], *(*cell.fab_array[n])[lev],
514 geom[lev + 1], geom[lev],
515 0, (*cell.fab_array[n])[lev]->nComp(), refRatio(lev));
516 }
517 SetFinestLevel(finest_level);
518 }
519 if (restart_file_cell != "")
520 {
521 Restart(restart_file_cell, false);
522 }
523 if (restart_file_node != "")
524 {
525 Restart(restart_file_node, true);
526 }
527
528 if (plot_int > 0 || plot_dt > 0.0) {
529 WritePlotFile();
530 }
531}
532
533void
534Integrator::Restart(const std::string dirname, bool a_nodal)
535{
536 BL_PROFILE("Integrator::Restart");
537
538 if (a_nodal && node.fab_array.size() == 0)
539 {
540 Util::Message(INFO, "Nothing here for nodal fabs");
541 return;
542 }
543 if (!a_nodal && cell.fab_array.size() == 0)
544 {
545 Util::Message(INFO, "Nothing here for cell-based fabs");
546 return;
547 }
548
549 std::string filename = dirname + "/Header";
550 std::string chkptfilename = dirname + "/Checkpoint";
551 amrex::VisMF::IO_Buffer io_buffer(amrex::VisMF::GetIOBufferSize());
552 amrex::Vector<char> fileCharPtr, chkptfileCharPtr;
553 amrex::ParallelDescriptor::ReadAndBcastFile(filename, fileCharPtr);
554 amrex::ParallelDescriptor::ReadAndBcastFile(chkptfilename, chkptfileCharPtr);
555 std::string fileCharPtrString(fileCharPtr.dataPtr());
556 std::string chkptfileCharPtrString(chkptfileCharPtr.dataPtr());
557 std::istringstream is(fileCharPtrString, std::istringstream::in);
558 std::istringstream chkpt_is(chkptfileCharPtrString, std::istringstream::in);
559
560 std::string line, word;
561
562 // Get version
563 std::getline(is, line);
564 Util::Message(INFO, "Version: ", line);
565
566 // Get number of fabs
567 int tmp_numfabs;
568 std::getline(is, line); tmp_numfabs = std::stoi(line);
569 Util::Message(INFO, "number of fabs:", tmp_numfabs);
570 std::vector<std::string> tmp_name_array;
571
572 for (int i = 0; i < tmp_numfabs; i++)
573 {
574 std::getline(is, line);
575 tmp_name_array.push_back(line);
576 }
577
578 // Dimension?
579 std::getline(is, line);
580 Util::Warning(INFO, "Dimension: " + line);
581
582 // Current time
583 Set::Scalar tmp_time = 0.0;
584 std::getline(is, line); tmp_time = std::stof(line); Util::Message(INFO, "Current time: ", tmp_time);
585 for (int i = 0; i < max_level + 1; i++)
586 {
587 t_new[i] = tmp_time; t_old[i] = tmp_time;
588 }
589
590 // AMR level
591 int tmp_max_level;
592 std::getline(is, line); tmp_max_level = std::stoi(line); Util::Message(INFO, "Max AMR level: ", line);
593 if (tmp_max_level > max_level)
594 Util::Abort(INFO, "The max level specified (", max_level, ") is smaller than the finest level in the restart file (", tmp_max_level, ")");
595 finest_level = tmp_max_level;
596 // Geometry ?
597 std::getline(is, line); Util::Message(INFO, "Input geometry: ", line);
598 std::getline(is, line); Util::Message(INFO, " ", line);
599
600 // Mesh refinement ratio?
601 std::getline(is, line); Util::Message(INFO, "Mesh refinement ratio: ", line);
602
603 // Domain
604 std::getline(is, line); Util::Warning(INFO, "Domain: ", line);
605
606 // Domain
607 std::getline(is, line);
608 std::vector<std::string> tmp_iters = Util::String::Split(line);
609 if (tmp_iters.size() != (unsigned int)(finest_level + 1)) Util::Abort(INFO, "Error reading in interation counts: line = ", line);
610 for (int lev = 0; lev <= finest_level; lev++) { istep[lev] = std::stoi(tmp_iters[lev]); Util::Message(INFO, "Iter on level ", lev, " = ", istep[lev]); }
611
612 amrex::Vector<amrex::MultiFab> tmpdata(tmp_max_level + 1);
613 int total_ncomp = 0;
614
615 if (a_nodal)
616 {
617 for (unsigned int i = 0; i < node.fab_array.size(); i++)
618 if (node.writeout_array[i])
619 total_ncomp += node.ncomp_array[i];
620 }
621 else
622 {
623 for (unsigned int i = 0; i < cell.fab_array.size(); i++)
624 if (cell.writeout_array[i])
625 total_ncomp += cell.ncomp_array[i];
626 }
627
628 int total_nghost = a_nodal ? 0 : cell.nghost_array[0];
629
630 for (int lev = 0; lev <= finest_level; lev++)
631 {
632 amrex::BoxArray tmp_ba;
633 tmp_ba.readFrom(chkpt_is);
634 SetBoxArray(lev, tmp_ba);
635 amrex::DistributionMapping tmp_dm(tmp_ba, amrex::ParallelDescriptor::NProcs());
636 SetDistributionMap(lev, tmp_dm);
637
638 if (a_nodal)
639 {
640 amrex::BoxArray ngrids = grids[lev];
641 ngrids.convert(amrex::IntVect::TheNodeVector());
642 //tmpdata[lev].define(ngrids, dmap[lev], total_ncomp, total_nghost);
643 }
644 else
645 {
646 //tmpdata[lev].define(grids[lev], dmap[lev], total_ncomp, total_nghost);
647 }
648 Util::Message(INFO,max_level);
649 Util::Message(INFO,finest_level);
650 Util::Message(INFO,lev,dirname);
651 Util::Message(INFO,lev,grids[lev]);
652 Util::Message(INFO,amrex::MultiFabFileFullPrefix(lev, dirname, "Level_", "Cell"));
653 Util::Message(INFO,total_ncomp);
654 Util::Message(INFO,total_nghost);
655 amrex::VisMF::Read(tmpdata[lev],
656 amrex::MultiFabFileFullPrefix(lev, dirname, "Level_", "Cell"));
657
658 if (a_nodal)
659 for (int i = 0; i < node.number_of_fabs; i++)
660 {
661 amrex::BoxArray ngrids = grids[lev];
662 ngrids.convert(amrex::IntVect::TheNodeVector());
663 (*node.fab_array[i])[lev].reset(new amrex::MultiFab(ngrids, dmap[lev], node.ncomp_array[i], node.nghost_array[i]));
664 (*node.fab_array[i])[lev]->setVal(0.);
665 }
666 else
667 for (int i = 0; i < cell.number_of_fabs; i++)
668 (*cell.fab_array[i])[lev].reset(new amrex::MultiFab(grids[lev], dmap[lev], cell.ncomp_array[i], cell.nghost_array[i]));
669 for (int i = 0; i < tmp_numfabs; i++)
670 {
671 bool match = false;
672 if (a_nodal)
673 {
674 for (int j = 0; j < node.number_of_fabs; j++)
675 {
676 if (tmp_name_array[i] == node.name_array[i][j])
677 {
678 match = true;
679 Util::Message(INFO, "Initializing ", node.name_array[i][j], "; nghost=", node.nghost_array[j], " with ", tmp_name_array[i]);
680 amrex::MultiFab::Copy(*((*node.fab_array[j])[lev]).get(), tmpdata[lev], i, 0, 1, total_nghost);
681 }
682 for (int k = 0; k < node.ncomp_array[j]; k++)
683 {
684 if (tmp_name_array[i] == node.name_array[j][k])
685 {
686 match = true;
687 Util::Message(INFO, "Initializing ", node.name_array[j][k], "; ncomp=", node.ncomp_array[j], "; nghost=", node.nghost_array[j], " with ", tmp_name_array[i]);
688 amrex::MultiFab::Copy(*((*node.fab_array[j])[lev]).get(), tmpdata[lev], i, k, 1, total_nghost);
689 }
690 }
691 Util::RealFillBoundary(*((*node.fab_array[j])[lev]).get(), geom[lev]);
692 }
693 }
694 else
695 {
696 for (int j = 0; j < cell.number_of_fabs; j++)
697 {
698 for (int k = 0; k < cell.ncomp_array[j]; k++)
699 {
700 if (tmp_name_array[i] == cell.name_array[j][k])
701 {
702 match = true;
703 Util::Message(INFO, "Initializing ", cell.name_array[j][k], "; ncomp=", cell.ncomp_array[j], "; nghost=", cell.nghost_array[j], " with ", tmp_name_array[i]);
704 amrex::MultiFab::Copy(*((*cell.fab_array[j])[lev]).get(), tmpdata[lev], i, k, 1, 0 /*cell.nghost_array[j]*/);
705 }
706 }
707 Util::RealFillBoundary(*(*cell.fab_array[j])[lev].get(), geom[lev]);
708 }
709 }
710 if (!match) Util::Warning(INFO, "Fab ", tmp_name_array[i], " is in the restart file, but there is no fab with that name here.");
711 }
712
713 for (unsigned int n = 0; n < m_basefields_cell.size(); n++)
714 {
715 m_basefields_cell[n]->MakeNewLevelFromScratch(lev, t_new[lev], grids[lev], dmap[lev]);
716 }
717 for (unsigned int n = 0; n < m_basefields.size(); n++)
718 {
719 m_basefields[n]->MakeNewLevelFromScratch(lev, t_new[lev], grids[lev], dmap[lev]);
720 }
721
722
723 for (int n = 0; n < cell.number_of_fabs; n++)
724 {
725 if (cell.writeout_array[n])
726 FillPatch(lev, t_new[lev], *cell.fab_array[n], *(*cell.fab_array[n])[lev], *cell.physbc_array[n], 0);
727 }
728
729 }
730
731 SetFinestLevel(max_level);
732}
733
734void
735Integrator::MakeNewLevelFromScratch(int lev, amrex::Real t, const amrex::BoxArray& cgrids,
736 const amrex::DistributionMapping& dm)
737{
738 BL_PROFILE("Integrator::MakeNewLevelFromScratch");
739 for (int n = 0; n < cell.number_of_fabs; n++)
740 {
741 (*cell.fab_array[n])[lev].reset(new amrex::MultiFab(cgrids, dm, cell.ncomp_array[n], cell.nghost_array[n]));
742 (*cell.fab_array[n])[lev]->setVal(0.0);
743 }
744
745 amrex::BoxArray ngrids = cgrids;
746 ngrids.convert(amrex::IntVect::TheNodeVector());
747 for (int n = 0; n < node.number_of_fabs; n++)
748 {
749 (*node.fab_array[n])[lev].reset(new amrex::MultiFab(ngrids, dm, node.ncomp_array[n], node.nghost_array[n]));
750 (*node.fab_array[n])[lev]->setVal(0.0);
751 }
752 for (unsigned int n = 0; n < m_basefields_cell.size(); n++)
753 {
754 m_basefields_cell[n]->MakeNewLevelFromScratch(lev, t, cgrids, dm);
755 }
756 for (unsigned int n = 0; n < m_basefields.size(); n++)
757 {
758 m_basefields[n]->MakeNewLevelFromScratch(lev, t, cgrids, dm);
759 }
760
761 t_new[lev] = t;
762 t_old[lev] = t - dt[lev];
763
764 Initialize(lev);
765
766 for (int n = 0; n < cell.number_of_fabs; n++)
767 {
768 cell.physbc_array[n]->define(geom[lev]);
769 cell.physbc_array[n]->FillBoundary(*(*cell.fab_array[n])[lev], 0, 0, t, 0);
770 }
771
772 //for (int n = 0 ; n < node.number_of_fabs; n++)
773 //{
774 // bcnothing->define(geom[lev]);
775 // for (amrex::MFIter mfi(*(*node.fab_array[n])[lev],true); mfi.isValid(); ++mfi)
776 // {
777 // amrex::BaseFab<Set::Scalar> &patch = (*(*node.fab_array[n])[lev])[mfi];
778 // const amrex::Box& box = mfi.tilebox();
779 // bcnothing->FillBoundary(patch,box,0,0,0,t);
780 // }
781 //}
782
783 for (unsigned int n = 0; n < m_basefields_cell.size(); n++)
784 {
785 m_basefields_cell[n]->FillBoundary(lev, t);
786 }
787 for (unsigned int n = 0; n < m_basefields.size(); n++)
788 {
789 m_basefields[n]->FillBoundary(lev, t);
790 }
791}
792
793std::vector<std::string>
794Integrator::PlotFileName(int lev, std::string prefix) const
795{
796 BL_PROFILE("Integrator::PlotFileName");
797 std::vector<std::string> name;
798 name.push_back(plot_file + "/" + prefix + "/");
799 name.push_back(amrex::Concatenate("", lev, 5));
800 return name;
801}
802
803void
804Integrator::WritePlotFile(bool initial) const
805{
806 WritePlotFile(t_new[0], istep, initial, "");
807}
808void
809Integrator::WritePlotFile(std::string prefix, Set::Scalar time, int step) const
810{
811 amrex::Vector<int> istep(max_level + 1, step);
812 WritePlotFile(time, istep, false, prefix);
813}
814
815void
816Integrator::WritePlotFile(Set::Scalar time, amrex::Vector<int> iter, bool initial, std::string prefix) const
817{
818 BL_PROFILE("Integrator::WritePlotFile");
819 int nlevels = finest_level + 1;
820 if (max_plot_level >= 0) nlevels = std::min(nlevels, max_plot_level);
821
822 int ccomponents = 0, ncomponents = 0, bfcomponents_cell = 0, bfcomponents = 0;
823 amrex::Vector<std::string> cnames, nnames, bfnames_cell, bfnames;
824 for (int i = 0; i < cell.number_of_fabs; i++)
825 {
826 if (!cell.writeout_array[i]) continue;
827 ccomponents += cell.ncomp_array[i];
828 if (cell.ncomp_array[i] > 1)
829 for (int j = 0; j < cell.ncomp_array[i]; j++)
830 cnames.push_back(cell.name_array[i][j]);
831 else
832 cnames.push_back(cell.name_array[i][0]);
833 }
834 for (int i = 0; i < node.number_of_fabs; i++)
835 {
836 if (!node.writeout_array[i]) continue;
837 ncomponents += node.ncomp_array[i];
838 if (node.ncomp_array[i] > 1)
839 for (int j = 0; j < node.ncomp_array[i]; j++)
840 nnames.push_back(node.name_array[i][j]);
841 else
842 nnames.push_back(node.name_array[i][0]);
843 }
844 for (unsigned int i = 0; i < m_basefields_cell.size(); i++)
845 {
846 if (m_basefields_cell[i]->writeout)
847 {
848 bfcomponents_cell += m_basefields_cell[i]->NComp();
849 for (int j = 0; j < m_basefields_cell[i]->NComp(); j++)
850 bfnames_cell.push_back(m_basefields_cell[i]->Name(j));
851 }
852 }
853 for (unsigned int i = 0; i < m_basefields.size(); i++)
854 {
855 if (m_basefields[i]->writeout)
856 {
857 bfcomponents += m_basefields[i]->NComp();
858 for (int j = 0; j < m_basefields[i]->NComp(); j++)
859 bfnames.push_back(m_basefields[i]->Name(j));
860 }
861 }
862
863 amrex::Vector<amrex::MultiFab> cplotmf(nlevels), nplotmf(nlevels);
864
865 bool do_cell_plotfile = (ccomponents + bfcomponents_cell > 0 || (ncomponents + bfcomponents > 0 && cell.all)) && cell.any;
866 bool do_node_plotfile = (ncomponents + bfcomponents > 0 || (ccomponents + bfcomponents_cell > 0 && node.all)) && node.any;
867
868 for (int ilev = 0; ilev < nlevels; ++ilev)
869 {
870 if (do_cell_plotfile)
871 {
872 int ncomp = ccomponents + bfcomponents_cell;
873 if (cell.all) ncomp += ncomponents + bfcomponents;
874 amrex::BoxArray cgrids_ghost = grids[ilev];
875 cgrids_ghost.grow(print_ghost_cells);
876 cplotmf[ilev].define(cgrids_ghost, dmap[ilev], ncomp, 0);
877
878 int n = 0;
879 int cnames_cnt = 0;
880 for (int i = 0; i < cell.number_of_fabs; i++)
881 {
882 if (!cell.writeout_array[i]) continue;
883 if ((*cell.fab_array[i])[ilev]->contains_nan())
884 {
885 if (abort_on_nan) Util::Abort(INFO, cnames[cnames_cnt], " contains nan (i=", i, ")");
886 else Util::Warning(INFO, cnames[cnames_cnt], " contains nan (i=", i, ")");
887 }
888 if ((*cell.fab_array[i])[ilev]->contains_inf())
889 {
890 if (abort_on_nan) Util::Abort(INFO, cnames[cnames_cnt], " contains inf (i=", i, ")");
891 else Util::Warning(INFO, cnames[cnames_cnt], " contains inf (i=", i, ")");
892 }
893 cnames_cnt++;
894 amrex::MultiFab::Copy(cplotmf[ilev], *(*cell.fab_array[i])[ilev], 0, n, cell.ncomp_array[i], 0);
895 n += cell.ncomp_array[i];
896 }
897 for (unsigned int i = 0; i < m_basefields_cell.size(); i++)
898 {
899 if (m_basefields_cell[i]->writeout)
900 {
901 m_basefields_cell[i]->Copy(ilev, cplotmf[ilev], n, 0);
902 n += m_basefields_cell[i]->NComp();
903 }
904 }
905
906 if (cell.all)
907 {
908 int nnames_cnt = 0;
909 for (int i = 0; i < node.number_of_fabs; i++)
910 {
911 if (!node.writeout_array[i]) continue;
912 if ((*node.fab_array[i])[ilev]->contains_nan())
913 {
914 if (abort_on_nan) Util::Abort(INFO, nnames[nnames_cnt], " contains nan (i=", i, ")");
915 else Util::Warning(INFO, nnames[nnames_cnt], " contains nan (i=", i, ")");
916 }
917 if ((*node.fab_array[i])[ilev]->contains_inf())
918 {
919 if (abort_on_nan) Util::Abort(INFO, nnames[nnames_cnt], " contains inf (i=", i, ")");
920 else Util::Warning(INFO, nnames[nnames_cnt], " contains inf (i=", i, ")");
921 }
922 nnames_cnt++;
923 amrex::average_node_to_cellcenter(cplotmf[ilev], n, *(*node.fab_array[i])[ilev], 0, node.ncomp_array[i], 0);
924 n += node.ncomp_array[i];
925 }
926 if (bfcomponents > 0)
927 {
928 amrex::BoxArray ngrids = grids[ilev];
929 ngrids.convert(amrex::IntVect::TheNodeVector());
930 amrex::MultiFab bfplotmf(ngrids, dmap[ilev], bfcomponents, 0);
931 int ctr = 0;
932 for (unsigned int i = 0; i < m_basefields.size(); i++)
933 {
934 if (m_basefields[i]->writeout)
935 {
936 m_basefields[i]->Copy(ilev, bfplotmf, ctr, 0);
937 ctr += m_basefields[i]->NComp();
938 }
939 }
940 amrex::average_node_to_cellcenter(cplotmf[ilev], n, bfplotmf, 0, bfcomponents);
941 n += bfcomponents;
942 }
943 }
944 }
945
946 if (do_node_plotfile)
947 {
948 amrex::BoxArray ngrids = grids[ilev];
949 ngrids.convert(amrex::IntVect::TheNodeVector());
950 int ncomp = ncomponents + bfcomponents;
951 if (node.all) ncomp += ccomponents + bfcomponents_cell;
952
953 amrex::BoxArray ngrids_ghost = ngrids;
954 ngrids_ghost.grow(print_ghost_nodes);
955
956 nplotmf[ilev].define(ngrids_ghost, dmap[ilev], ncomp, 0);
957
958 int n = 0;
959 for (int i = 0; i < node.number_of_fabs; i++)
960 {
961 if (!node.writeout_array[i]) continue;
962 if ((*node.fab_array[i])[ilev]->contains_nan()) Util::Warning(INFO, nnames[i], " contains nan (i=", i, "). Resetting to zero.");
963 if ((*node.fab_array[i])[ilev]->contains_inf()) Util::Abort(INFO, nnames[i], " contains inf (i=", i, ")");
964 amrex::MultiFab::Copy(nplotmf[ilev], *(*node.fab_array[i])[ilev], 0, n, node.ncomp_array[i], 0);
965 n += node.ncomp_array[i];
966 }
967 for (unsigned int i = 0; i < m_basefields.size(); i++)
968 {
969 if (m_basefields[i]->writeout)
970 {
971 m_basefields[i]->Copy(ilev, nplotmf[ilev], n, 0);
972 n += m_basefields[i]->NComp();
973 }
974 }
975
976 if (node.all)
977 {
978 for (int i = 0; i < cell.number_of_fabs; i++)
979 {
980 if (!cell.writeout_array[i]) continue;
981 if ((*cell.fab_array[i])[ilev]->contains_nan()) Util::Warning(INFO, cnames[i], " contains nan (i=", i, "). Resetting to zero.");
982 if ((*cell.fab_array[i])[ilev]->contains_inf()) Util::Abort(INFO, cnames[i], " contains inf (i=", i, ")");
983 if ((*cell.fab_array[i])[ilev]->nGrow() < 1)
984 {
985 if (initial) Util::Warning(INFO, cnames[i], " has no ghost cells and will not be included in nodal output");
986 continue;
987 }
988 Util::AverageCellcenterToNode(nplotmf[ilev], n, *(*cell.fab_array[i])[ilev], 0, cell.ncomp_array[i]);
989 n += cell.ncomp_array[i];
990 }
991
992 if (bfcomponents_cell > 0)
993 {
994 amrex::BoxArray cgrids = grids[ilev];
995 amrex::MultiFab bfplotmf(cgrids, dmap[ilev], bfcomponents_cell, 0);
996 int ctr = 0;
997 for (unsigned int i = 0; i < m_basefields_cell.size(); i++)
998 {
999 if (m_basefields_cell[i]->writeout)
1000 {
1001 m_basefields_cell[i]->Copy(ilev, bfplotmf, ctr, 0);
1002 ctr += m_basefields_cell[i]->NComp();
1003 }
1004 }
1005 Util::AverageCellcenterToNode(nplotmf[ilev], n, bfplotmf, 0, bfcomponents_cell);
1006 n += bfcomponents_cell;
1007 }
1008 }
1009 }
1010 }
1011
1012 std::vector<std::string> plotfilename = PlotFileName(istep[0], prefix);
1013 if (initial) plotfilename[1] = plotfilename[1] + "init";
1014
1015 if (do_cell_plotfile)
1016 {
1017 amrex::Vector<std::string> allnames = cnames;
1018 allnames.insert(allnames.end(), bfnames_cell.begin(), bfnames_cell.end());
1019 if (cell.all) {
1020 allnames.insert(allnames.end(), nnames.begin(), nnames.end());
1021 allnames.insert(allnames.end(), bfnames.begin(), bfnames.end());
1022 }
1023 const std::string base = plotfilename[0] + plotfilename[1] + "cell";
1024#ifdef AMREX_USE_HDF5
1025 WriteMultiLevelPlotfileHDF5(base, nlevels, amrex::GetVecOfConstPtrs(cplotmf), allnames,
1026 Geom(), time, iter, refRatio(),"ZLIB@" + compression_level);
1027 const std::string chkptfilename = base + ".Checkpoint";
1028#else
1029 WriteMultiLevelPlotfile(base, nlevels, amrex::GetVecOfConstPtrs(cplotmf), allnames,
1030 Geom(), time, iter, refRatio());
1031 const std::string chkptfilename = base + "/Checkpoint";
1032#endif
1033 std::ofstream chkptfile(chkptfilename);
1034 if (!chkptfile.good()) amrex::FileOpenFailed(chkptfilename);
1035 for (int i = 0; i <= max_level; i++) boxArray(i).writeOn(chkptfile);
1036 chkptfile.close();
1037 if (chkptfile.fail()) amrex::FileOpenFailed(chkptfilename);
1038 }
1039
1040 if (do_node_plotfile)
1041 {
1042 amrex::Vector<std::string> allnames = nnames;
1043 allnames.insert(allnames.end(), bfnames.begin(), bfnames.end());
1044 if (node.all) allnames.insert(allnames.end(), cnames.begin(), cnames.end());
1045 const std::string base = plotfilename[0] + plotfilename[1] + "node";
1046#ifdef AMREX_USE_HDF5
1047 WriteMultiLevelPlotfileHDF5(base, nlevels, amrex::GetVecOfConstPtrs(nplotmf), allnames,
1048 Geom(), time, iter, refRatio(),"ZLIB@" + compression_level);
1049 const std::string chkptfilename = base + ".Checkpoint";
1050#else
1051 WriteMultiLevelPlotfile(base, nlevels, amrex::GetVecOfConstPtrs(nplotmf), allnames,
1052 Geom(), time, iter, refRatio());
1053 const std::string chkptfilename = base + "/Checkpoint";
1054#endif
1055 std::ofstream chkptfile(chkptfilename);
1056 if (!chkptfile.good()) amrex::FileOpenFailed(chkptfilename);
1057 for (int i = 0; i <= max_level; i++) boxArray(i).writeOn(chkptfile);
1058 chkptfile.close();
1059 if (chkptfile.fail()) amrex::FileOpenFailed(chkptfilename);
1060 }
1061
1062 if (amrex::ParallelDescriptor::IOProcessor())
1063 {
1064 std::ofstream coutfile, noutfile;
1065 if (istep[0] == 0)
1066 {
1067 if (do_cell_plotfile) coutfile.open(plot_file + "/celloutput.visit", std::ios_base::out);
1068 if (do_node_plotfile) noutfile.open(plot_file + "/nodeoutput.visit", std::ios_base::out);
1069 }
1070 else
1071 {
1072 if (do_cell_plotfile) coutfile.open(plot_file + "/celloutput.visit", std::ios_base::app);
1073 if (do_node_plotfile) noutfile.open(plot_file + "/nodeoutput.visit", std::ios_base::app);
1074 }
1075#ifdef AMREX_USE_HDF5
1076 const std::string header_suffix = ".h5";
1077#else
1078 const std::string header_suffix = "/Header";
1079#endif
1080 if (do_cell_plotfile) coutfile << plotfilename[1] + "cell" + header_suffix << std::endl;
1081 if (do_node_plotfile) noutfile << plotfilename[1] + "node" + header_suffix << std::endl;
1082 }
1083}
1084
1085void
1086Integrator::Evolve()
1087{
1088 BL_PROFILE("Integrator::Evolve");
1089 amrex::Real cur_time = t_new[0];
1090 int last_plot_file_step = 0;
1091
1092 for (int step = istep[0]; step < max_step && cur_time < stop_time; ++step)
1093 {
1094 if (amrex::ParallelDescriptor::IOProcessor()) {
1095 std::cout << "\nSTEP " << step + 1 << " starts ..." << std::endl;
1096 }
1097 int lev = 0;
1098 int iteration = 1;
1099 TimeStepBegin(cur_time, step);
1100 if (integrate_variables_before_advance) IntegrateVariables(cur_time, step);
1101 TimeStep(lev, cur_time, iteration);
1102 if (integrate_variables_after_advance) IntegrateVariables(cur_time, step);
1103 TimeStepComplete(cur_time, step);
1104 cur_time += dt[0];
1105
1106 if (amrex::ParallelDescriptor::IOProcessor()) {
1107 std::cout << "STEP " << step + 1 << " ends."
1108 << " TIME = " << cur_time << " DT = " << dt[0]
1109 << std::endl;
1110 }
1111
1112 // sync up time
1113 for (int lev = 0; lev <= finest_level; ++lev) {
1114 t_new[lev] = cur_time;
1115 }
1116
1117 if (plot_int > 0 && (step + 1) % plot_int == 0) {
1118 last_plot_file_step = step + 1;
1119 WritePlotFile();
1120 IO::WriteMetaData(plot_file, IO::Status::Running, (int)(100.0 * cur_time / stop_time));
1121 }
1122 else if (std::fabs(std::remainder(cur_time, plot_dt)) < 0.5 * dt[0])
1123 {
1124 last_plot_file_step = step + 1;
1125 WritePlotFile();
1126 IO::WriteMetaData(plot_file, IO::Status::Running, (int)(100.0 * cur_time / stop_time));
1127 }
1128
1129 if (cur_time >= stop_time - 1.e-6 * dt[0]) break;
1130 }
1131 if (plot_int > 0 && istep[0] > last_plot_file_step) {
1132 WritePlotFile();
1133 }
1134}
1135
1136void
1137Integrator::IntegrateVariables(amrex::Real time, int step)
1138{
1139 BL_PROFILE("Integrator::IntegrateVariables");
1140 if (!thermo.number) return;
1141
1142 if ((thermo.interval > 0 && (step) % thermo.interval == 0) ||
1143 ((thermo.dt > 0.0) && (std::fabs(std::remainder(time, plot_dt)) < 0.5 * dt[0])))
1144 {
1145 // Zero out all variables
1146 for (int i = 0; i < thermo.number; i++)
1147 {
1148 if (thermo.extensives[i]) *thermo.vars[i] = 0;
1149 }
1150
1151 // All levels except the finest
1152 for (int ilev = 0; ilev < max_level; ilev++)
1153 {
1154 const amrex::BoxArray& cfba = amrex::coarsen(grids[ilev + 1], refRatio(ilev));
1155
1156#ifdef OMP
1157#pragma omp parallel
1158#endif
1159 for (amrex::MFIter mfi(grids[ilev], dmap[ilev], true); mfi.isValid(); ++mfi)
1160 {
1161 const amrex::Box& box = mfi.tilebox();
1162 const amrex::BoxArray& comp = amrex::complementIn(box, cfba);
1163
1164 for (int i = 0; i < comp.size(); i++)
1165 {
1166 Integrate(ilev, time, step,
1167 mfi, comp[i]);
1168 }
1169 }
1170 }
1171 // Now do the finest level
1172 {
1173#ifdef OMP
1174#pragma omp parallel
1175#endif
1176 for (amrex::MFIter mfi(grids[max_level], dmap[max_level], true); mfi.isValid(); ++mfi)
1177 {
1178 const amrex::Box& box = mfi.tilebox();
1179 Integrate(max_level, time, step, mfi, box);
1180 }
1181 }
1182
1183 // Sum up across all processors
1184 for (int i = 0; i < thermo.number; i++)
1185 {
1186 if (thermo.extensives[i])
1187 amrex::ParallelDescriptor::ReduceRealSum(*thermo.vars[i]);
1188 }
1189 }
1190 if (amrex::ParallelDescriptor::IOProcessor() &&
1191 (
1192 (thermo.plot_int > 0 && step % thermo.plot_int == 0) ||
1193 (thermo.plot_dt > 0.0 && std::fabs(std::remainder(time, thermo.plot_dt)) < 0.5 * dt[0])
1194 ))
1195 {
1196 std::ofstream outfile;
1197 if (step == 0)
1198 {
1199 outfile.open(plot_file + "/thermo.dat", std::ios_base::out);
1200 outfile << "time";
1201 for (int i = 0; i < thermo.number; i++)
1202 outfile << "\t" << thermo.names[i];
1203 outfile << std::endl;
1204 }
1205 else outfile.open(plot_file + "/thermo.dat", std::ios_base::app);
1206 outfile << time;
1207 for (int i = 0; i < thermo.number; i++)
1208 outfile << "\t" << *thermo.vars[i];
1209 outfile << std::endl;
1210 outfile.close();
1211 }
1212
1213}
1214
1215
1216void
1217Integrator::TimeStep(int lev, amrex::Real time, int /*iteration*/)
1218{
1219 BL_PROFILE("Integrator::TimeStep");
1220 if (base_regrid_int <= 0 || istep[0] % base_regrid_int == 0)
1221 {
1222 if (regrid_int > 0 || base_regrid_int > 0) // We may need to regrid
1223 {
1224 static amrex::Vector<int> last_regrid_step(max_level + 1, 0);
1225
1226 // regrid doesn't change the base level, so we don't regrid on max_level
1227 if (lev < max_level && istep[lev] > last_regrid_step[lev])
1228 {
1229 if (istep[lev] % regrid_int == 0)
1230 {
1231 regrid(lev, time, false);
1232 }
1233 }
1234 }
1235 }
1236 SetFinestLevel(finest_level);
1237
1238 if (Verbose() && amrex::ParallelDescriptor::IOProcessor()) {
1239 std::cout << "[Level " << lev
1240 << " step " << istep[lev] + 1 << "] ";
1241 std::cout << "ADVANCE with dt = "
1242 << dt[lev]
1243 << std::endl;
1244 }
1245
1246 for (int n = 0; n < cell.number_of_fabs; n++)
1247 if (cell.evolving_array[n]) FillPatch(lev, time, *cell.fab_array[n], *(*cell.fab_array[n])[lev], *cell.physbc_array[n], 0);
1248 for (int n = 0; n < node.number_of_fabs; n++)
1249 if (node.evolving_array[n]) FillPatch(lev, time, *node.fab_array[n], *(*node.fab_array[n])[lev], *node.physbc_array[n], 0);
1250 for (unsigned int n = 0; n < m_basefields_cell.size(); n++)
1251 if (m_basefields_cell[n]->evolving) m_basefields_cell[n]->FillPatch(lev, time);
1252 for (unsigned int n = 0; n < m_basefields.size(); n++)
1253 if (m_basefields[n]->evolving) m_basefields[n]->FillPatch(lev, time);
1254
1255 Advance(lev, time, dt[lev]);
1256 ++istep[lev];
1257
1258 if (Verbose() && amrex::ParallelDescriptor::IOProcessor())
1259 {
1260 std::cout << "[Level " << lev
1261 << " step " << istep[lev] << "] ";
1262 std::cout << "Advanced "
1263 << CountCells(lev)
1264 << " cells"
1265 << std::endl;
1266 }
1267
1268 if (lev < finest_level)
1269 {
1270 for (int i = 1; i <= nsubsteps[lev + 1]; ++i)
1271 TimeStep(lev + 1, time + (i - 1) * dt[lev + 1], i);
1272
1273 for (int n = 0; n < cell.number_of_fabs; n++)
1274 {
1275 amrex::average_down(*(*cell.fab_array[n])[lev + 1], *(*cell.fab_array[n])[lev],
1276 geom[lev + 1], geom[lev],
1277 0, (*cell.fab_array[n])[lev]->nComp(), refRatio(lev));
1278 }
1279 for (int n = 0; n < node.number_of_fabs; n++)
1280 {
1281 amrex::average_down(*(*node.fab_array[n])[lev + 1], *(*node.fab_array[n])[lev],
1282 0, (*node.fab_array[n])[lev]->nComp(), refRatio(lev));
1283 }
1284 for (unsigned int n = 0; n < m_basefields_cell.size(); n++)
1285 {
1286 if (m_basefields_cell[n]->evolving)
1287 m_basefields_cell[n]->AverageDown(lev, refRatio(lev));
1288 }
1289 for (unsigned int n = 0; n < m_basefields.size(); n++)
1290 {
1291 if (m_basefields[n]->evolving)
1292 m_basefields[n]->AverageDown(lev, refRatio(lev));
1293 }
1294
1295 }
1296}
1297}
t
std::time_t t
Definition FileNameParse.cpp:9
TEST
#define TEST(x)
Definition Util.H:25
INFO
#define INFO
Definition Util.H:24
BC::BC
Definition BC.H:43
BC::BC::define
void define(const amrex::Geometry &a_geom)
Definition BC.H:48
BC::BC::GetBCRec
virtual amrex::BCRec GetBCRec()=0
IO::ParmParse::contains
bool contains(std::string name)
Definition ParmParse.H:173
IO::ParmParse::queryarr
int queryarr(std::string name, std::vector< T > &value)
Definition ParmParse.H:535
IO::ParmParse::query_required
int query_required(std::string name, T &value)
Definition ParmParse.H:249
IO::ParmParse::query_default
int query_default(std::string name, T &value, T defaultvalue)
Definition ParmParse.H:293
IO::ParmParse::query_validate
int query_validate(std::string name, int &value, std::vector< int > possibleintvals)
Definition ParmParse.H:336
Integrator::Integrator::WritePlotFile
void WritePlotFile(bool initial=false) const
Definition Integrator.cpp:804
Integrator::Integrator::Integrator
Integrator()
This is the constructor for the intetgrator class, which reads timestep information,...
Definition Integrator.cpp:18
Integrator::Integrator::PlotFileName
std::vector< std::string > PlotFileName(int lev, std::string prefix="") const
Definition Integrator.cpp:794
Integrator::Integrator::FillCoarsePatch
void FillCoarsePatch(int lev, amrex::Real time, Set::Field< Set::Scalar > &mf, BC::BC< Set::Scalar > &physbc, int icomp, int ncomp)
Fill a fab at current level with the data from one level up.
Definition Integrator.cpp:461
Integrator::Integrator::RegisterNodalFab
void RegisterNodalFab(Set::Field< Set::Scalar > &new_fab, int ncomp, int nghost, std::string name, bool writeout, bool evolving=true, std::vector< std::string > suffix={})
Add a new node-based scalar field.
Definition Integrator.cpp:366
Integrator::Integrator::RegisterNewFab
void RegisterNewFab(Set::Field< Set::Scalar > &new_fab, BC::BC< Set::Scalar > *new_bc, int ncomp, int nghost, std::string name, bool writeout, bool evolving=true, std::vector< std::string > suffix={})
Add a new cell-based scalar field.
Definition Integrator.cpp:348
Integrator::Integrator::Restart
void Restart(std::string restartfile, bool a_node=false)
Read in output from previous simulation and start simulation at that point - Not currently tested.
Definition Integrator.cpp:534
Integrator::Integrator::TimeStep
void TimeStep(int lev, amrex::Real time, int iteration)
Timestep marching.
Definition Integrator.cpp:1217
Integrator::Integrator::MakeNewLevelFromScratch
virtual void MakeNewLevelFromScratch(int lev, amrex::Real time, const amrex::BoxArray &ba, const amrex::DistributionMapping &dm) override
Definition Integrator.cpp:735
Integrator::Integrator::ClearLevel
virtual void ClearLevel(int lev) override
Definition Integrator.cpp:326
Integrator::Integrator::CountCells
long CountCells(int lev)
Simple utility to count cells.
Definition Integrator.cpp:386
Integrator::Integrator::FillPatch
void FillPatch(int lev, amrex::Real time, amrex::Vector< std::unique_ptr< amrex::MultiFab > > &source_mf, amrex::MultiFab &destination_multifab, BC::BC< Set::Scalar > &physbc, int icomp)
This is the function that is responsible for updating patch data.
Definition Integrator.cpp:402
Integrator::Integrator::~Integrator
virtual ~Integrator()
Virtual destructure; make sure delete any pointers that you create here.
Definition Integrator.cpp:186
Integrator::Integrator::InitData
void InitData()
Front-end method to initialize simulation on all levels.
Definition Integrator.cpp:500
Integrator::Integrator::Evolve
void Evolve()
Front-end method to start simulation.
Definition Integrator.cpp:1086
Integrator::Integrator::SetPlotInt
void SetPlotInt(int plot_int)
Utility to set the frequency (in timesteps) of plotfile dumping.
Definition Integrator.cpp:215
Integrator::Integrator::ErrorEst
virtual void ErrorEst(int lev, amrex::TagBoxArray &tags, amrex::Real time, int ngrow) override
Definition Integrator.cpp:492
Integrator::Integrator::IntegrateVariables
void IntegrateVariables(Set::Scalar cur_time, int step)
Definition Integrator.cpp:1137
Integrator::Integrator::RegisterIntegratedVariable
void RegisterIntegratedVariable(Set::Scalar *integrated_variable, std::string name, bool extensive=true)
Register a variable to be integrated over the spatial domain using the Integrate function.
Definition Integrator.cpp:376
Integrator::Integrator::SetTimestep
void SetTimestep(Set::Scalar _timestep)
Utility to set the coarse-grid timestep.
Definition Integrator.cpp:205
Integrator::Integrator::MakeNewLevelFromCoarse
virtual void MakeNewLevelFromCoarse(int lev, amrex::Real time, const amrex::BoxArray &ba, const amrex::DistributionMapping &dm) override
Wrapper to call FillCoarsePatch.
Definition Integrator.cpp:226
Integrator::Integrator::RemakeLevel
virtual void RemakeLevel(int lev, amrex::Real time, const amrex::BoxArray &ba, const amrex::DistributionMapping &dm) override
RESETS ALL MULTIFABS AT A GIVEN LEVEL.
Definition Integrator.cpp:275
Integrator::Integrator::Parse
static void Parse(Integrator &, IO::ParmParse &)
Definition Integrator.cpp:24
IO::WriteMetaData
void WriteMetaData(std::string plot_file, Status status, int per)
Definition WriteMetaData.cpp:17
IO::FileNameParse
void FileNameParse(std::string &filename)
Internal function to do processing of the file name.
Definition FileNameParse.cpp:12
IO::Complete
@ Complete
Definition WriteMetaData.H:31
IO::Running
@ Running
Definition WriteMetaData.H:30
Integrator
Collection of numerical integrator objects.
Definition AllenCahn.H:43
Integrator::AddField< Set::Scalar, Set::Hypercube::Node >
ALAMO_SINGLE_DEFINITION void Integrator::AddField< Set::Scalar, Set::Hypercube::Node >(Set::Field< Set::Scalar > &new_field, BC::BC< Set::Scalar > *new_bc, int ncomp, int nghost, std::string name, bool writeout, bool evolving, std::vector< std::string > suffix)
Definition Integrator.H:524
Integrator::AddField< Set::Scalar, Set::Hypercube::Cell >
ALAMO_SINGLE_DEFINITION void Integrator::AddField< Set::Scalar, Set::Hypercube::Cell >(Set::Field< Set::Scalar > &new_field, BC::BC< Set::Scalar > *new_bc, int ncomp, int nghost, std::string name, bool writeout, bool evolving, std::vector< std::string > suffix)
Definition Integrator.H:484
Set::Scalar
amrex::Real Scalar
Definition Base.H:19
Set::Vector
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, 1 > Vector
Definition Base.H:21
Util::String::Split
AMREX_FORCE_INLINE std::vector< std::string > Split(std::string &str, const char delim=' ')
Definition String.H:138
Util::GetFileName
std::string GetFileName()
Definition Util.cpp:32
Util::AssertException
AMREX_FORCE_INLINE void AssertException(std::string file, std::string func, int line, std::string smt, bool pass, Args const &... args)
Definition Util.H:254
Util::finalized
bool finalized
Definition Util.cpp:30
Util::Abort
void Abort(const char *msg)
Definition Util.cpp:268
Util::AverageCellcenterToNode
void AverageCellcenterToNode(amrex::MultiFab &node_mf, const int &dcomp, const amrex::MultiFab &cell_mf, const int &scomp, const int &ncomp)
Definition Util.cpp:407
Util::Assert
AMREX_FORCE_INLINE void Assert(std::string file, std::string func, int line, std::string smt, bool pass, Args const &... args)
Definition Util.H:58
Util::Warning
void Warning(std::string file, std::string func, int line, Args const &... args)
Definition Util.H:202
Util::Message
void Message(std::string file, std::string func, int line, Args const &... args)
Definition Util.H:129
Util::RealFillBoundary
AMREX_FORCE_INLINE void RealFillBoundary(amrex::FabArray< amrex::BaseFab< T > > &a_mf, const amrex::Geometry &, const int nghost=2)
Definition Util.H:339
Unit::Time
static Unit Time()
Definition Unit.H:199
Generated by doxygen 1.9.8