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