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