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