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