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