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