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