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, ") does not match the max 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)(max_level + 1)) Util::Abort(INFO, "Error reading in interation counts: line = ", line);
594 0 : for (int lev = 0; lev <= max_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) for (unsigned int i = 0; i < node.fab_array.size(); i++) total_ncomp += node.ncomp_array[i];
600 0 : else for (unsigned int i = 0; i < cell.fab_array.size(); i++) total_ncomp += cell.ncomp_array[i];
601 :
602 0 : int total_nghost = a_nodal ? 0 : cell.nghost_array[0];
603 :
604 0 : for (int lev = 0; lev <= max_level; lev++)
605 : {
606 0 : amrex::BoxArray tmp_ba;
607 0 : tmp_ba.readFrom(chkpt_is);
608 0 : SetBoxArray(lev, tmp_ba);
609 0 : amrex::DistributionMapping tmp_dm(tmp_ba, amrex::ParallelDescriptor::NProcs());
610 0 : SetDistributionMap(lev, tmp_dm);
611 :
612 0 : if (a_nodal)
613 : {
614 0 : amrex::BoxArray ngrids = grids[lev];
615 0 : ngrids.convert(amrex::IntVect::TheNodeVector());
616 0 : tmpdata[lev].define(ngrids, dmap[lev], total_ncomp, total_nghost);
617 0 : }
618 : else
619 : {
620 0 : tmpdata[lev].define(grids[lev], dmap[lev], total_ncomp, total_nghost);
621 : }
622 0 : amrex::VisMF::Read(tmpdata[lev],
623 0 : amrex::MultiFabFileFullPrefix(lev, dirname, "Level_", "Cell"));
624 :
625 :
626 0 : if (a_nodal)
627 0 : for (int i = 0; i < node.number_of_fabs; i++)
628 : {
629 0 : amrex::BoxArray ngrids = grids[lev];
630 0 : ngrids.convert(amrex::IntVect::TheNodeVector());
631 0 : (*node.fab_array[i])[lev].reset(new amrex::MultiFab(ngrids, dmap[lev], node.ncomp_array[i], node.nghost_array[i]));
632 0 : (*node.fab_array[i])[lev]->setVal(0.);
633 0 : }
634 : else
635 0 : for (int i = 0; i < cell.number_of_fabs; i++)
636 0 : (*cell.fab_array[i])[lev].reset(new amrex::MultiFab(grids[lev], dmap[lev], cell.ncomp_array[i], cell.nghost_array[i]));
637 0 : for (int i = 0; i < tmp_numfabs; i++)
638 : {
639 0 : bool match = false;
640 0 : if (a_nodal)
641 : {
642 0 : for (int j = 0; j < node.number_of_fabs; j++)
643 : {
644 0 : if (tmp_name_array[i] == node.name_array[i][j])
645 : {
646 0 : match = true;
647 0 : Util::Message(INFO, "Initializing ", node.name_array[i][j], "; nghost=", node.nghost_array[j], " with ", tmp_name_array[i]);
648 0 : amrex::MultiFab::Copy(*((*node.fab_array[j])[lev]).get(), tmpdata[lev], i, 0, 1, total_nghost);
649 : }
650 0 : for (int k = 0; k < node.ncomp_array[j]; k++)
651 : {
652 0 : if (tmp_name_array[i] == node.name_array[j][k])
653 : {
654 0 : match = true;
655 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]);
656 0 : amrex::MultiFab::Copy(*((*node.fab_array[j])[lev]).get(), tmpdata[lev], i, k, 1, total_nghost);
657 : }
658 : }
659 0 : Util::RealFillBoundary(*((*node.fab_array[j])[lev]).get(), geom[lev]);
660 : }
661 : }
662 : else
663 : {
664 0 : for (int j = 0; j < cell.number_of_fabs; j++)
665 0 : for (int k = 0; k < cell.ncomp_array[j]; k++)
666 : {
667 0 : if (tmp_name_array[i] == cell.name_array[j][k])
668 : {
669 0 : match = true;
670 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]);
671 0 : amrex::MultiFab::Copy(*((*cell.fab_array[j])[lev]).get(), tmpdata[lev], i, k, 1, cell.nghost_array[j]);
672 : }
673 : }
674 : }
675 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.");
676 : }
677 :
678 0 : for (unsigned int n = 0; n < m_basefields_cell.size(); n++)
679 : {
680 0 : m_basefields_cell[n]->MakeNewLevelFromScratch(lev, t_new[lev], grids[lev], dmap[lev]);
681 : }
682 0 : for (unsigned int n = 0; n < m_basefields.size(); n++)
683 : {
684 0 : m_basefields[n]->MakeNewLevelFromScratch(lev, t_new[lev], grids[lev], dmap[lev]);
685 : }
686 0 : }
687 :
688 0 : SetFinestLevel(max_level);
689 0 : }
690 :
691 : void
692 139 : Integrator::MakeNewLevelFromScratch(int lev, amrex::Real t, const amrex::BoxArray& cgrids,
693 : const amrex::DistributionMapping& dm)
694 : {
695 : BL_PROFILE("Integrator::MakeNewLevelFromScratch");
696 464 : for (int n = 0; n < cell.number_of_fabs; n++)
697 : {
698 325 : (*cell.fab_array[n])[lev].reset(new amrex::MultiFab(cgrids, dm, cell.ncomp_array[n], cell.nghost_array[n]));
699 325 : (*cell.fab_array[n])[lev]->setVal(0.0);
700 : }
701 :
702 139 : amrex::BoxArray ngrids = cgrids;
703 139 : ngrids.convert(amrex::IntVect::TheNodeVector());
704 226 : for (int n = 0; n < node.number_of_fabs; n++)
705 : {
706 87 : (*node.fab_array[n])[lev].reset(new amrex::MultiFab(ngrids, dm, node.ncomp_array[n], node.nghost_array[n]));
707 87 : (*node.fab_array[n])[lev]->setVal(0.0);
708 : }
709 139 : for (unsigned int n = 0; n < m_basefields_cell.size(); n++)
710 : {
711 0 : m_basefields_cell[n]->MakeNewLevelFromScratch(lev, t, cgrids, dm);
712 : }
713 534 : for (unsigned int n = 0; n < m_basefields.size(); n++)
714 : {
715 395 : m_basefields[n]->MakeNewLevelFromScratch(lev, t, cgrids, dm);
716 : }
717 :
718 139 : t_new[lev] = t;
719 139 : t_old[lev] = t - dt[lev];
720 :
721 139 : Initialize(lev);
722 :
723 464 : for (int n = 0; n < cell.number_of_fabs; n++)
724 : {
725 325 : cell.physbc_array[n]->define(geom[lev]);
726 325 : cell.physbc_array[n]->FillBoundary(*(*cell.fab_array[n])[lev], 0, 0, t, 0);
727 : }
728 :
729 : //for (int n = 0 ; n < node.number_of_fabs; n++)
730 : //{
731 : // bcnothing->define(geom[lev]);
732 : // for (amrex::MFIter mfi(*(*node.fab_array[n])[lev],true); mfi.isValid(); ++mfi)
733 : // {
734 : // amrex::BaseFab<Set::Scalar> &patch = (*(*node.fab_array[n])[lev])[mfi];
735 : // const amrex::Box& box = mfi.tilebox();
736 : // bcnothing->FillBoundary(patch,box,0,0,0,t);
737 : // }
738 : //}
739 :
740 139 : for (unsigned int n = 0; n < m_basefields_cell.size(); n++)
741 : {
742 0 : m_basefields_cell[n]->FillBoundary(lev, t);
743 : }
744 534 : for (unsigned int n = 0; n < m_basefields.size(); n++)
745 : {
746 395 : m_basefields[n]->FillBoundary(lev, t);
747 : }
748 139 : }
749 :
750 : std::vector<std::string>
751 250 : Integrator::PlotFileName(int lev, std::string prefix) const
752 : {
753 : BL_PROFILE("Integrator::PlotFileName");
754 250 : std::vector<std::string> name;
755 250 : name.push_back(plot_file + "/" + prefix + "/");
756 250 : name.push_back(amrex::Concatenate("", lev, 5));
757 250 : return name;
758 0 : }
759 :
760 : void
761 250 : Integrator::WritePlotFile(bool initial) const
762 : {
763 250 : WritePlotFile(t_new[0], istep, initial, "");
764 250 : }
765 : void
766 0 : Integrator::WritePlotFile(std::string prefix, Set::Scalar time, int step) const
767 : {
768 0 : amrex::Vector<int> istep(max_level + 1, step);
769 0 : WritePlotFile(time, istep, false, prefix);
770 0 : }
771 :
772 : void
773 250 : Integrator::WritePlotFile(Set::Scalar time, amrex::Vector<int> iter, bool initial, std::string prefix) const
774 : {
775 : BL_PROFILE("Integrator::WritePlotFile");
776 250 : int nlevels = finest_level + 1;
777 250 : if (max_plot_level >= 0) nlevels = std::min(nlevels, max_plot_level);
778 :
779 250 : int ccomponents = 0, ncomponents = 0, bfcomponents_cell = 0, bfcomponents = 0;
780 250 : amrex::Vector<std::string> cnames, nnames, bfnames_cell, bfnames;
781 1226 : for (int i = 0; i < cell.number_of_fabs; i++)
782 : {
783 976 : if (!cell.writeout_array[i]) continue;
784 731 : ccomponents += cell.ncomp_array[i];
785 731 : if (cell.ncomp_array[i] > 1)
786 737 : for (int j = 0; j < cell.ncomp_array[i]; j++)
787 518 : cnames.push_back(cell.name_array[i][j]);
788 : else
789 512 : cnames.push_back(cell.name_array[i][0]);
790 : }
791 335 : for (int i = 0; i < node.number_of_fabs; i++)
792 : {
793 85 : if (!node.writeout_array[i]) continue;
794 85 : ncomponents += node.ncomp_array[i];
795 85 : if (node.ncomp_array[i] > 1)
796 18 : for (int j = 0; j < node.ncomp_array[i]; j++)
797 12 : nnames.push_back(node.name_array[i][j]);
798 : else
799 79 : nnames.push_back(node.name_array[i][0]);
800 : }
801 250 : for (unsigned int i = 0; i < m_basefields_cell.size(); i++)
802 : {
803 0 : if (m_basefields_cell[i]->writeout)
804 : {
805 0 : bfcomponents_cell += m_basefields_cell[i]->NComp();
806 0 : for (int j = 0; j < m_basefields_cell[i]->NComp(); j++)
807 0 : bfnames_cell.push_back(m_basefields_cell[i]->Name(j));
808 : }
809 : }
810 677 : for (unsigned int i = 0; i < m_basefields.size(); i++)
811 : {
812 427 : if (m_basefields[i]->writeout)
813 : {
814 427 : bfcomponents += m_basefields[i]->NComp();
815 2057 : for (int j = 0; j < m_basefields[i]->NComp(); j++)
816 1630 : bfnames.push_back(m_basefields[i]->Name(j));
817 : }
818 : }
819 :
820 750 : amrex::Vector<amrex::MultiFab> cplotmf(nlevels), nplotmf(nlevels);
821 :
822 250 : bool do_cell_plotfile = (ccomponents + bfcomponents_cell > 0 || (ncomponents + bfcomponents > 0 && cell.all)) && cell.any;
823 250 : bool do_node_plotfile = (ncomponents + bfcomponents > 0 || (ccomponents + bfcomponents_cell > 0 && node.all)) && node.any;
824 :
825 953 : for (int ilev = 0; ilev < nlevels; ++ilev)
826 : {
827 703 : if (do_cell_plotfile)
828 : {
829 653 : int ncomp = ccomponents + bfcomponents_cell;
830 653 : if (cell.all) ncomp += ncomponents + bfcomponents;
831 653 : cplotmf[ilev].define(grids[ilev], dmap[ilev], ncomp, 0);
832 :
833 653 : int n = 0;
834 653 : int cnames_cnt = 0;
835 2480 : for (int i = 0; i < cell.number_of_fabs; i++)
836 : {
837 1827 : if (!cell.writeout_array[i]) continue;
838 1188 : if ((*cell.fab_array[i])[ilev]->contains_nan())
839 : {
840 0 : if (abort_on_nan) Util::Abort(INFO, cnames[cnames_cnt], " contains nan (i=", i, ")");
841 0 : else Util::Warning(INFO, cnames[cnames_cnt], " contains nan (i=", i, ")");
842 : }
843 1188 : if ((*cell.fab_array[i])[ilev]->contains_inf())
844 : {
845 0 : if (abort_on_nan) Util::Abort(INFO, cnames[cnames_cnt], " contains inf (i=", i, ")");
846 0 : else Util::Warning(INFO, cnames[cnames_cnt], " contains inf (i=", i, ")");
847 : }
848 1188 : cnames_cnt++;
849 1188 : amrex::MultiFab::Copy(cplotmf[ilev], *(*cell.fab_array[i])[ilev], 0, n, cell.ncomp_array[i], 0);
850 1188 : n += cell.ncomp_array[i];
851 : }
852 653 : for (unsigned int i = 0; i < m_basefields_cell.size(); i++)
853 : {
854 0 : if (m_basefields_cell[i]->writeout)
855 : {
856 0 : m_basefields_cell[i]->Copy(ilev, cplotmf[ilev], n, 0);
857 0 : n += m_basefields_cell[i]->NComp();
858 : }
859 : }
860 :
861 653 : if (cell.all)
862 : {
863 108 : int nnames_cnt = 0;
864 216 : for (int i = 0; i < node.number_of_fabs; i++)
865 : {
866 108 : if (!node.writeout_array[i]) continue;
867 108 : if ((*node.fab_array[i])[ilev]->contains_nan())
868 : {
869 0 : if (abort_on_nan) Util::Abort(INFO, nnames[nnames_cnt], " contains nan (i=", i, ")");
870 0 : else Util::Warning(INFO, nnames[nnames_cnt], " contains nan (i=", i, ")");
871 : }
872 108 : if ((*node.fab_array[i])[ilev]->contains_inf())
873 : {
874 0 : if (abort_on_nan) Util::Abort(INFO, nnames[nnames_cnt], " contains inf (i=", i, ")");
875 0 : else Util::Warning(INFO, nnames[nnames_cnt], " contains inf (i=", i, ")");
876 : }
877 108 : nnames_cnt++;
878 108 : amrex::average_node_to_cellcenter(cplotmf[ilev], n, *(*node.fab_array[i])[ilev], 0, node.ncomp_array[i], 0);
879 108 : n += node.ncomp_array[i];
880 : }
881 108 : if (bfcomponents > 0)
882 : {
883 108 : amrex::BoxArray ngrids = grids[ilev];
884 108 : ngrids.convert(amrex::IntVect::TheNodeVector());
885 108 : amrex::MultiFab bfplotmf(ngrids, dmap[ilev], bfcomponents, 0);
886 108 : int ctr = 0;
887 672 : for (unsigned int i = 0; i < m_basefields.size(); i++)
888 : {
889 564 : if (m_basefields[i]->writeout)
890 : {
891 564 : m_basefields[i]->Copy(ilev, bfplotmf, ctr, 0);
892 564 : ctr += m_basefields[i]->NComp();
893 : }
894 : }
895 108 : amrex::average_node_to_cellcenter(cplotmf[ilev], n, bfplotmf, 0, bfcomponents);
896 108 : n += bfcomponents;
897 108 : }
898 : }
899 : }
900 :
901 703 : if (do_node_plotfile)
902 : {
903 175 : amrex::BoxArray ngrids = grids[ilev];
904 175 : ngrids.convert(amrex::IntVect::TheNodeVector());
905 175 : int ncomp = ncomponents + bfcomponents;
906 175 : if (node.all) ncomp += ccomponents + bfcomponents_cell;
907 175 : nplotmf[ilev].define(ngrids, dmap[ilev], ncomp, 0);
908 :
909 175 : int n = 0;
910 345 : for (int i = 0; i < node.number_of_fabs; i++)
911 : {
912 170 : if (!node.writeout_array[i]) continue;
913 170 : if ((*node.fab_array[i])[ilev]->contains_nan()) Util::Warning(INFO, nnames[i], " contains nan (i=", i, "). Resetting to zero.");
914 170 : if ((*node.fab_array[i])[ilev]->contains_inf()) Util::Abort(INFO, nnames[i], " contains inf (i=", i, ")");
915 170 : amrex::MultiFab::Copy(nplotmf[ilev], *(*node.fab_array[i])[ilev], 0, n, node.ncomp_array[i], 0);
916 170 : n += node.ncomp_array[i];
917 : }
918 994 : for (unsigned int i = 0; i < m_basefields.size(); i++)
919 : {
920 819 : if (m_basefields[i]->writeout)
921 : {
922 819 : m_basefields[i]->Copy(ilev, nplotmf[ilev], n, 0);
923 819 : n += m_basefields[i]->NComp();
924 : }
925 : }
926 :
927 175 : if (node.all)
928 : {
929 129 : for (int i = 0; i < cell.number_of_fabs; i++)
930 : {
931 67 : if (!cell.writeout_array[i]) continue;
932 62 : if ((*cell.fab_array[i])[ilev]->contains_nan()) Util::Warning(INFO, cnames[i], " contains nan (i=", i, "). Resetting to zero.");
933 62 : if ((*cell.fab_array[i])[ilev]->contains_inf()) Util::Abort(INFO, cnames[i], " contains inf (i=", i, ")");
934 62 : if ((*cell.fab_array[i])[ilev]->nGrow() == 0)
935 : {
936 0 : if (initial) Util::Warning(INFO, cnames[i], " has no ghost cells and will not be included in nodal output");
937 0 : continue;
938 : }
939 62 : Util::AverageCellcenterToNode(nplotmf[ilev], n, *(*cell.fab_array[i])[ilev], 0, cell.ncomp_array[i]);
940 62 : n += cell.ncomp_array[i];
941 : }
942 :
943 62 : if (bfcomponents_cell > 0)
944 : {
945 0 : amrex::BoxArray cgrids = grids[ilev];
946 0 : amrex::MultiFab bfplotmf(cgrids, dmap[ilev], bfcomponents_cell, 0);
947 0 : int ctr = 0;
948 0 : for (unsigned int i = 0; i < m_basefields_cell.size(); i++)
949 : {
950 0 : if (m_basefields_cell[i]->writeout)
951 : {
952 0 : m_basefields_cell[i]->Copy(ilev, bfplotmf, ctr, 0);
953 0 : ctr += m_basefields_cell[i]->NComp();
954 : }
955 : }
956 0 : Util::AverageCellcenterToNode(nplotmf[ilev], n, bfplotmf, 0, bfcomponents_cell);
957 0 : n += bfcomponents_cell;
958 0 : }
959 : }
960 175 : }
961 : }
962 :
963 250 : std::vector<std::string> plotfilename = PlotFileName(istep[0], prefix);
964 250 : if (initial) plotfilename[1] = plotfilename[1] + "init";
965 :
966 250 : if (do_cell_plotfile)
967 : {
968 204 : amrex::Vector<std::string> allnames = cnames;
969 204 : allnames.insert(allnames.end(), bfnames_cell.begin(), bfnames_cell.end());
970 204 : if (cell.all) {
971 35 : allnames.insert(allnames.end(), nnames.begin(), nnames.end());
972 35 : allnames.insert(allnames.end(), bfnames.begin(), bfnames.end());
973 : }
974 1428 : WriteMultiLevelPlotfile(plotfilename[0] + plotfilename[1] + "cell", nlevels, amrex::GetVecOfConstPtrs(cplotmf), allnames,
975 : Geom(), time, iter, refRatio());
976 :
977 204 : std::ofstream chkptfile;
978 204 : chkptfile.open(plotfilename[0] + plotfilename[1] + "cell/Checkpoint");
979 864 : for (int i = 0; i <= max_level; i++) boxArray(i).writeOn(chkptfile);
980 204 : chkptfile.close();
981 204 : }
982 :
983 250 : if (do_node_plotfile)
984 : {
985 86 : amrex::Vector<std::string> allnames = nnames;
986 86 : allnames.insert(allnames.end(), bfnames.begin(), bfnames.end());
987 86 : if (node.all) allnames.insert(allnames.end(), cnames.begin(), cnames.end());
988 602 : WriteMultiLevelPlotfile(plotfilename[0] + plotfilename[1] + "node", nlevels, amrex::GetVecOfConstPtrs(nplotmf), allnames,
989 : Geom(), time, iter, refRatio());
990 :
991 86 : std::ofstream chkptfile;
992 86 : chkptfile.open(plotfilename[0] + plotfilename[1] + "node/Checkpoint");
993 265 : for (int i = 0; i <= max_level; i++) boxArray(i).writeOn(chkptfile);
994 86 : chkptfile.close();
995 86 : }
996 :
997 250 : if (amrex::ParallelDescriptor::IOProcessor())
998 : {
999 250 : std::ofstream coutfile, noutfile;
1000 250 : if (istep[0] == 0)
1001 : {
1002 37 : if (do_cell_plotfile) coutfile.open(plot_file + "/celloutput.visit", std::ios_base::out);
1003 37 : if (do_node_plotfile) noutfile.open(plot_file + "/nodeoutput.visit", std::ios_base::out);
1004 : }
1005 : else
1006 : {
1007 213 : if (do_cell_plotfile) coutfile.open(plot_file + "/celloutput.visit", std::ios_base::app);
1008 213 : if (do_node_plotfile) noutfile.open(plot_file + "/nodeoutput.visit", std::ios_base::app);
1009 : }
1010 250 : if (do_cell_plotfile) coutfile << plotfilename[1] + "cell" + "/Header" << std::endl;
1011 250 : if (do_node_plotfile) noutfile << plotfilename[1] + "node" + "/Header" << std::endl;
1012 250 : }
1013 250 : }
1014 :
1015 : void
1016 37 : Integrator::Evolve()
1017 : {
1018 : BL_PROFILE("Integrator::Evolve");
1019 37 : amrex::Real cur_time = t_new[0];
1020 37 : int last_plot_file_step = 0;
1021 :
1022 12040 : for (int step = istep[0]; step < max_step && cur_time < stop_time; ++step)
1023 : {
1024 12040 : if (amrex::ParallelDescriptor::IOProcessor()) {
1025 12040 : std::cout << "\nSTEP " << step + 1 << " starts ..." << std::endl;
1026 : }
1027 12040 : int lev = 0;
1028 12040 : int iteration = 1;
1029 12040 : TimeStepBegin(cur_time, step);
1030 12040 : if (integrate_variables_before_advance) IntegrateVariables(cur_time, step);
1031 12040 : TimeStep(lev, cur_time, iteration);
1032 12040 : if (integrate_variables_after_advance) IntegrateVariables(cur_time, step);
1033 12040 : TimeStepComplete(cur_time, step);
1034 12040 : cur_time += dt[0];
1035 :
1036 12040 : if (amrex::ParallelDescriptor::IOProcessor()) {
1037 12040 : std::cout << "STEP " << step + 1 << " ends."
1038 12040 : << " TIME = " << cur_time << " DT = " << dt[0]
1039 12040 : << std::endl;
1040 : }
1041 :
1042 : // sync up time
1043 46111 : for (int lev = 0; lev <= finest_level; ++lev) {
1044 34071 : t_new[lev] = cur_time;
1045 : }
1046 :
1047 12040 : if (plot_int > 0 && (step + 1) % plot_int == 0) {
1048 133 : last_plot_file_step = step + 1;
1049 133 : WritePlotFile();
1050 133 : IO::WriteMetaData(plot_file, IO::Status::Running, (int)(100.0 * cur_time / stop_time));
1051 : }
1052 11907 : else if (std::fabs(std::remainder(cur_time, plot_dt)) < 0.5 * dt[0])
1053 : {
1054 76 : last_plot_file_step = step + 1;
1055 76 : WritePlotFile();
1056 76 : IO::WriteMetaData(plot_file, IO::Status::Running, (int)(100.0 * cur_time / stop_time));
1057 : }
1058 :
1059 12040 : if (cur_time >= stop_time - 1.e-6 * dt[0]) break;
1060 : }
1061 37 : if (plot_int > 0 && istep[0] > last_plot_file_step) {
1062 4 : WritePlotFile();
1063 : }
1064 37 : }
1065 :
1066 : void
1067 12040 : Integrator::IntegrateVariables(amrex::Real time, int step)
1068 : {
1069 : BL_PROFILE("Integrator::IntegrateVariables");
1070 12040 : if (!thermo.number) return;
1071 :
1072 480 : if ((thermo.interval > 0 && (step) % thermo.interval == 0) ||
1073 0 : ((thermo.dt > 0.0) && (std::fabs(std::remainder(time, plot_dt)) < 0.5 * dt[0])))
1074 : {
1075 : // Zero out all variables
1076 4146 : for (int i = 0; i < thermo.number; i++)
1077 : {
1078 3666 : if (thermo.extensives[i]) *thermo.vars[i] = 0;
1079 : }
1080 :
1081 : // All levels except the finest
1082 696 : for (int ilev = 0; ilev < max_level; ilev++)
1083 : {
1084 216 : const amrex::BoxArray& cfba = amrex::coarsen(grids[ilev + 1], refRatio(ilev));
1085 :
1086 : #ifdef OMP
1087 : #pragma omp parallel
1088 : #endif
1089 6205 : for (amrex::MFIter mfi(grids[ilev], dmap[ilev], true); mfi.isValid(); ++mfi)
1090 : {
1091 5989 : const amrex::Box& box = mfi.tilebox();
1092 5989 : const amrex::BoxArray& comp = amrex::complementIn(box, cfba);
1093 :
1094 13579 : for (int i = 0; i < comp.size(); i++)
1095 : {
1096 7590 : Integrate(ilev, time, step,
1097 15180 : mfi, comp[i]);
1098 : }
1099 6205 : }
1100 216 : }
1101 : // Now do the finest level
1102 : {
1103 : #ifdef OMP
1104 : #pragma omp parallel
1105 : #endif
1106 5864 : for (amrex::MFIter mfi(grids[max_level], dmap[max_level], true); mfi.isValid(); ++mfi)
1107 : {
1108 5384 : const amrex::Box& box = mfi.tilebox();
1109 5384 : Integrate(max_level, time, step, mfi, box);
1110 480 : }
1111 : }
1112 :
1113 : // Sum up across all processors
1114 4146 : for (int i = 0; i < thermo.number; i++)
1115 : {
1116 3666 : if (thermo.extensives[i])
1117 3666 : amrex::ParallelDescriptor::ReduceRealSum(*thermo.vars[i]);
1118 : }
1119 : }
1120 960 : if (amrex::ParallelDescriptor::IOProcessor() &&
1121 : (
1122 480 : (thermo.plot_int > 0 && step % thermo.plot_int == 0) ||
1123 69 : (thermo.plot_dt > 0.0 && std::fabs(std::remainder(time, thermo.plot_dt)) < 0.5 * dt[0])
1124 : ))
1125 : {
1126 414 : std::ofstream outfile;
1127 414 : if (step == 0)
1128 : {
1129 8 : outfile.open(plot_file + "/thermo.dat", std::ios_base::out);
1130 8 : outfile << "time";
1131 73 : for (int i = 0; i < thermo.number; i++)
1132 65 : outfile << "\t" << thermo.names[i];
1133 8 : outfile << std::endl;
1134 : }
1135 406 : else outfile.open(plot_file + "/thermo.dat", std::ios_base::app);
1136 414 : outfile << time;
1137 3731 : for (int i = 0; i < thermo.number; i++)
1138 3317 : outfile << "\t" << *thermo.vars[i];
1139 414 : outfile << std::endl;
1140 414 : outfile.close();
1141 414 : }
1142 :
1143 : }
1144 :
1145 :
1146 : void
1147 110612 : Integrator::TimeStep(int lev, amrex::Real time, int /*iteration*/)
1148 : {
1149 : BL_PROFILE("Integrator::TimeStep");
1150 110612 : if (base_regrid_int <= 0 || istep[0] % base_regrid_int == 0)
1151 : {
1152 102942 : if (regrid_int > 0 || base_regrid_int > 0) // We may need to regrid
1153 : {
1154 102851 : static amrex::Vector<int> last_regrid_step(max_level + 1, 0);
1155 :
1156 : // regrid doesn't change the base level, so we don't regrid on max_level
1157 102816 : if (lev < max_level && istep[lev] > last_regrid_step[lev])
1158 : {
1159 45735 : if (istep[lev] % regrid_int == 0)
1160 : {
1161 19967 : regrid(lev, time, false);
1162 : }
1163 : }
1164 : }
1165 : }
1166 110612 : SetFinestLevel(finest_level);
1167 :
1168 110612 : if (Verbose() && amrex::ParallelDescriptor::IOProcessor()) {
1169 0 : std::cout << "[Level " << lev
1170 0 : << " step " << istep[lev] + 1 << "] ";
1171 0 : std::cout << "ADVANCE with dt = "
1172 0 : << dt[lev]
1173 0 : << std::endl;
1174 : }
1175 :
1176 410536 : for (int n = 0; n < cell.number_of_fabs; n++)
1177 299924 : if (cell.evolving_array[n]) FillPatch(lev, time, *cell.fab_array[n], *(*cell.fab_array[n])[lev], *cell.physbc_array[n], 0);
1178 119530 : for (int n = 0; n < node.number_of_fabs; n++)
1179 8918 : if (node.evolving_array[n]) FillPatch(lev, time, *node.fab_array[n], *(*node.fab_array[n])[lev], *node.physbc_array[n], 0);
1180 110612 : for (unsigned int n = 0; n < m_basefields_cell.size(); n++)
1181 0 : if (m_basefields_cell[n]->evolving) m_basefields_cell[n]->FillPatch(lev, time);
1182 114334 : for (unsigned int n = 0; n < m_basefields.size(); n++)
1183 3722 : if (m_basefields[n]->evolving) m_basefields[n]->FillPatch(lev, time);
1184 :
1185 110612 : Advance(lev, time, dt[lev]);
1186 110612 : ++istep[lev];
1187 :
1188 110612 : if (Verbose() && amrex::ParallelDescriptor::IOProcessor())
1189 : {
1190 0 : std::cout << "[Level " << lev
1191 0 : << " step " << istep[lev] << "] ";
1192 0 : std::cout << "Advanced "
1193 0 : << CountCells(lev)
1194 0 : << " cells"
1195 0 : << std::endl;
1196 : }
1197 :
1198 110612 : if (lev < finest_level)
1199 : {
1200 147858 : for (int i = 1; i <= nsubsteps[lev + 1]; ++i)
1201 98572 : TimeStep(lev + 1, time + (i - 1) * dt[lev + 1], i);
1202 :
1203 156564 : for (int n = 0; n < cell.number_of_fabs; n++)
1204 : {
1205 214556 : amrex::average_down(*(*cell.fab_array[n])[lev + 1], *(*cell.fab_array[n])[lev],
1206 107278 : geom[lev + 1], geom[lev],
1207 214556 : 0, (*cell.fab_array[n])[lev]->nComp(), refRatio(lev));
1208 : }
1209 53124 : for (int n = 0; n < node.number_of_fabs; n++)
1210 : {
1211 3838 : amrex::average_down(*(*node.fab_array[n])[lev + 1], *(*node.fab_array[n])[lev],
1212 7676 : 0, (*node.fab_array[n])[lev]->nComp(), refRatio(lev));
1213 : }
1214 49286 : for (unsigned int n = 0; n < m_basefields_cell.size(); n++)
1215 : {
1216 0 : if (m_basefields_cell[n]->evolving)
1217 0 : m_basefields_cell[n]->AverageDown(lev, refRatio(lev));
1218 : }
1219 50080 : for (unsigned int n = 0; n < m_basefields.size(); n++)
1220 : {
1221 794 : if (m_basefields[n]->evolving)
1222 0 : m_basefields[n]->AverageDown(lev, refRatio(lev));
1223 : }
1224 :
1225 : }
1226 110612 : }
1227 : }
|
