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