Line data Source code
1 : //
2 : // Pure abstract class for managing data structures, time integration (with substepping),
3 : // mesh refinement, and I/O.
4 : //
5 : // Native input file parameters:
6 : //
7 : // .. code-block:: make
8 : //
9 : // max_step = [maximum number of timesteps]
10 : // stop_time = [maximum simulation time]
11 : // timestep = [time step for coarsest level]
12 : //
13 : // amr.regrid_int = [number of timesteps between regridding]
14 : // amr.plot_int = [number of timesteps between dumping output]
15 : // amr.plot_file = [base name of output directory]
16 : //
17 : // amr.nsubsteps = [number of temporal substeps at each level. This can be
18 : // either a single int (which is then applied to every refinement
19 : // level) or an array of ints (equal to amr.max_level)
20 : // corresponding to the refinement for each level.]
21 : //
22 : // Inherited input file parameters (from amrex AmrMesh class):
23 : //
24 : // .. code-block:: make
25 : //
26 : // amr.v = [verbosity level]
27 : // amr.max_level = [maximum level of refinement]
28 : // amr.n_proper =
29 : // amr.grid_eff =
30 : // amr.n_error_buff =
31 : // amr.ref_ratio_vect = [refinement ratios in each direction]
32 : // amr.ref_ratio = [refinement ratio in all directions (cannot be used with ref_ratio_vect)]
33 : // amr.max_grid_x =
34 : // amr.blocking_factor =
35 : // amr.n_cell = [number of cells on coarsest level]
36 : // amr.refine_grid_layout =
37 : // amr.check_input =
38 : //
39 :
40 : #ifndef INTEGRATOR_INTEGRATOR_H
41 : #define INTEGRATOR_INTEGRATOR_H
42 :
43 : #include <chrono>
44 : #include <ctime>
45 : #include <string>
46 : #include <limits>
47 : #include <memory>
48 :
49 : #ifdef _OPENMP
50 : #include <omp.h>
51 : #endif
52 :
53 : #include <AMReX_ParallelDescriptor.H>
54 : #include <AMReX_ParmParse.H>
55 : #include <AMReX_MultiFabUtil.H>
56 : #include <AMReX_FillPatchUtil.H>
57 : #include <AMReX_BC_TYPES.H>
58 : #include <AMReX_AmrCore.H>
59 : #include <AMReX_FluxRegister.H>
60 : #include <AMReX_Utility.H>
61 : #include <AMReX_PlotFileUtil.H>
62 :
63 : #include "Set/Set.H"
64 : #include "BC/BC.H"
65 : #include "BC/Nothing.H"
66 : #include "IO/WriteMetaData.H"
67 : #include "BaseField.H"
68 :
69 : /// \brief Collection of numerical integrator objects
70 : namespace Integrator
71 : {
72 :
73 : class Integrator
74 : : public amrex::AmrCore
75 : {
76 : public:
77 :
78 : /// \function Integrator
79 : /// \brief Constructor
80 : ///
81 : /// Does the following things:
82 : /// - Read in simulation TIME(STEP) information
83 : /// - Read in simulation output and AMR information
84 : /// - Initalize timestep substep information
85 : /// - Create a clean directory
86 : /// For derived classes this **must** be called for the derived constructor. For instance: `code`
87 : /// ```cpp
88 : /// class MyDerivedClass : Integrator
89 : /// {
90 : /// MyDerivedClass() : Integrator() { ... }
91 : /// ...
92 : /// }
93 : /// ```
94 : Integrator();
95 :
96 : virtual ~Integrator();
97 :
98 : /// \fn FrontData
99 : /// \brief Front-end method to initialize simulation
100 : void InitData();
101 :
102 : void Restart(std::string restartfile, bool a_node = false);
103 :
104 : /// \fn Evolve
105 : /// \brief Front-end method to start simulation
106 : void Evolve();
107 :
108 : void SetFilename(std::string _plot_file) { plot_file = _plot_file; };
109 : std::string GetFilename() { return plot_file; };
110 :
111 2196 : void regrid(int lbase, Set::Scalar time, bool initial = false) override
112 : {
113 2196 : if (!explicitmesh.on)
114 2188 : AmrCore::regrid(lbase, time, initial);
115 2196 : }
116 :
117 31 : void InitFromScratch(Set::Scalar time)
118 : {
119 31 : if (!explicitmesh.on) AmrCore::InitFromScratch(time);
120 : else
121 : {
122 : // Generate the coarse level mesh.
123 : {
124 7 : finest_level = 0;
125 14 : const amrex::BoxArray& ba = MakeBaseGrids();
126 14 : amrex::DistributionMapping dm(ba);
127 7 : const auto old_num_setdm = num_setdm;
128 7 : const auto old_num_setba = num_setba;
129 7 : MakeNewLevelFromScratch(0, time, ba, dm);
130 7 : if (old_num_setba == num_setba) SetBoxArray(0, ba);
131 7 : if (old_num_setdm == num_setdm) SetDistributionMap(0, dm);
132 : }
133 : // Generate subsequent level meshes based on user input
134 13 : for (int ilev = 0; ilev < maxLevel(); ++ilev)
135 : {
136 6 : finest_level = ilev + 1;
137 12 : amrex::BoxArray grids(explicitmesh.box[ilev]);
138 6 : ChopGrids(ilev + 1, grids, amrex::ParallelDescriptor::NProcs());
139 12 : amrex::DistributionMapping dmap(grids);
140 6 : SetBoxArray(ilev + 1, grids);
141 6 : SetDistributionMap(ilev + 1, dmap);
142 6 : MakeNewLevelFromScratch(ilev + 1, time, grids, dmap);
143 : }
144 :
145 : }
146 31 : }
147 :
148 : protected:
149 :
150 : /// \fn Initialize
151 : /// \brief Apply initial conditions
152 : ///
153 : /// You **must** override this function to inherit this class.
154 : /// This function is called before the simulation begins, and is where
155 : /// initial conditions should be applied.
156 : virtual void Initialize(int lev ///<[in] AMR Level
157 : ) = 0;
158 :
159 : /// \fn Advance
160 : /// \brief Perform computation
161 : ///
162 : /// You **must** override this function to inherit this class.
163 : /// Advance is called every time(sub)step, and implements the evolution of
164 : /// the system in time.
165 : ///
166 : virtual void Advance(int lev, ///<[in] AMR Level
167 : amrex::Real time, ///< [in] System time
168 : amrex::Real dt ///< [in] Timestep for this level
169 : ) = 0;
170 :
171 : /// \fn TagCellsForRefinement
172 : /// \brief Tag cells where mesh refinement is needed
173 : ///
174 : /// You **must** override this function to inherit this class.
175 : /// Advance is called every time(sub)step, and implements the evolution of
176 : /// the system in time.
177 : ///
178 : virtual void TagCellsForRefinement(int lev, amrex::TagBoxArray& tags, amrex::Real time,
179 : int ngrow) = 0;
180 :
181 : /// \fn TimeStepBegin
182 : /// \brief Run another system calculation (e.g. implicit solve) before integration step
183 : ///
184 : /// This function is called at the beginning of every timestep. This function can be used
185 : /// to complete additional global solves, e.g. a MLMG implicit solve.
186 : ///
187 : /// Overriding is optional; the default is to do nothing.
188 : ///
189 1606 : virtual void TimeStepBegin(Set::Scalar /*time*/, int /*iter*/) {};
190 :
191 : /// \fn TimeStepComplete
192 : /// \brief Run another system calculation (e.g. implicit solve) after integration step
193 : ///
194 : /// This function is called at the end of every timestep. This function can be used
195 : /// to complete additional global solves, e.g. a MLMG implicit solve.
196 : ///
197 : /// Overriding is optional; the default is to do nothing.
198 : ///
199 2032 : virtual void TimeStepComplete(Set::Scalar /*time*/, int /*iter*/) {};
200 :
201 : /// \fn Integrate
202 : /// \brief Perform an integration to compute integrated quantities
203 : ///
204 : /// This is a function that is called by `Integrator` to update the variables registered in
205 : /// RegisterIntegratedVariable.
206 : /// The following variables are used:
207 : /// - amrlev: current amr level
208 : /// - time: current simulation time
209 : /// - iter: current simulation iteration
210 : /// - mfi: current MFIter object (used to get FArrayBox from MultiFab)
211 : /// - box: Use this box (not mfi.tilebox). This box covers only cells on this level that are
212 : /// not also on a finer level.
213 0 : virtual void Integrate(int /*amrlev*/, Set::Scalar /*time*/, int /*iter*/,
214 : const amrex::MFIter&/*mfi*/, const amrex::Box&/*box*/)
215 : {
216 0 : if (thermo.number > 0)
217 0 : Util::Warning(INFO, "integrated variables registered, but no integration implemented!");
218 0 : }
219 :
220 409 : virtual void Regrid(int /* amrlev */, Set::Scalar /* time */)
221 409 : {}
222 :
223 :
224 :
225 :
226 : /// \fn RegisterNewFab
227 : /// \brief Register a field variable for AMR with this class
228 : void RegisterNewFab(Set::Field<Set::Scalar>& new_fab, BC::BC<Set::Scalar>* new_bc, int ncomp, int nghost, std::string name, bool writeout);
229 : void RegisterNewFab(Set::Field<Set::Scalar>& new_fab, int ncomp, std::string name, bool writeout);
230 : void RegisterNodalFab(Set::Field<Set::Scalar>& new_fab, int ncomp, int nghost, std::string name, bool writeout);
231 : void RegisterNodalFab(Set::Field<Set::Scalar>& new_fab, BC::BC<Set::Scalar>* new_bc, int ncomp, int nghost, std::string name, bool writeout);
232 : template<class T>
233 : void RegisterGeneralFab(Set::Field<T>& new_fab, int ncomp, int nghost, bool evolving = true);
234 : template<class T>
235 : void RegisterGeneralFab(Set::Field<T>& new_fab, int ncomp, int nghost, std::string a_name, bool evolving = true);
236 : template<class T>
237 : void RegisterGeneralFab(Set::Field<T>& new_fab, int ncomp, int nghost, bool writeout, std::string a_name, bool evolving = true);
238 :
239 : template<class T, int d>
240 : void AddField(Set::Field<T>& new_field, BC::BC<T>* new_bc, int ncomp, int nghost, std::string, bool writeout, bool evolving);
241 :
242 70407 : void SetFinestLevel(const int a_finestlevel)
243 : {
244 232949 : for (unsigned int i = 0; i < cell.fab_array.size(); i++)
245 162542 : cell.fab_array[i]->finest_level = a_finestlevel;
246 125200 : for (unsigned int i = 0; i < node.fab_array.size(); i++)
247 54793 : node.fab_array[i]->finest_level = a_finestlevel;
248 70407 : for (unsigned int i = 0; i < m_basefields_cell.size(); i++)
249 0 : m_basefields_cell[i]->SetFinestLevel(finest_level);
250 74227 : for (unsigned int i = 0; i < m_basefields.size(); i++)
251 3820 : m_basefields[i]->SetFinestLevel(finest_level);
252 70407 : }
253 :
254 : void RegisterIntegratedVariable(Set::Scalar* integrated_variable, std::string name, bool extensive=true);
255 :
256 : void SetTimestep(Set::Scalar _timestep);
257 : void SetPlotInt(int plot_int);
258 1 : void SetThermoInt(int a_thermo_int) { thermo.interval = a_thermo_int; }
259 : void SetThermoPlotInt(int a_thermo_plot_int) { thermo.plot_int = a_thermo_plot_int; }
260 : void SetStopTime(Set::Scalar a_stop_time) { stop_time = a_stop_time; }
261 :
262 : amrex::Vector<amrex::Real> t_new; ///< Keep track of current old simulation time on each level
263 : amrex::Vector<int> istep; ///< Keep track of where each level is
264 : // PLOT FILES
265 : std::string plot_file{ "plt" }; ///< Plotfile name
266 :
267 : private:
268 : virtual void MakeNewLevelFromScratch(int lev, amrex::Real time, const amrex::BoxArray& ba,
269 : const amrex::DistributionMapping& dm) override;
270 : virtual void MakeNewLevelFromCoarse(int lev, amrex::Real time, const amrex::BoxArray& ba,
271 : const amrex::DistributionMapping& dm) override;
272 : virtual void RemakeLevel(int lev, amrex::Real time, const amrex::BoxArray& ba,
273 : const amrex::DistributionMapping& dm) override;
274 : virtual void ClearLevel(int lev) override;
275 : virtual void ErrorEst(int lev, amrex::TagBoxArray& tags, amrex::Real time, int ngrow) override;
276 :
277 :
278 : void FillPatch(int lev, amrex::Real time,
279 : amrex::Vector<std::unique_ptr<amrex::MultiFab>>& source_mf,
280 : amrex::MultiFab& destination_multifab,
281 : BC::BC<Set::Scalar>& physbc,
282 : int icomp);
283 : long CountCells(int lev);
284 : void TimeStep(int lev, amrex::Real time, int iteration);
285 : void FillCoarsePatch(int lev, amrex::Real time, Set::Field<Set::Scalar>& mf, BC::BC<Set::Scalar>& physbc, int icomp, int ncomp);
286 : void GetData(const int lev, const amrex::Real time, amrex::Vector<amrex::MultiFab*>& data, amrex::Vector<amrex::Real>& datatime);
287 :
288 : std::vector<std::string> PlotFileName(int lev, std::string prefix = "") const;
289 : protected:
290 : void IntegrateVariables(Set::Scalar cur_time, int step);
291 : void WritePlotFile(bool initial = false) const;
292 : void WritePlotFile(std::string prefix, Set::Scalar time, int step) const;
293 : void WritePlotFile(Set::Scalar time, amrex::Vector<int> iter, bool initial = false, std::string prefix = "") const;
294 :
295 : //
296 : // MEMBER VARIABLES
297 : //
298 :
299 : // TIME (STEP) KEEPINGamrex::Vector<std::unique_ptr<amrex::MultiFab> >
300 : protected:
301 : amrex::Real timestep = NAN; ///< Timestep for the base level of refinement
302 : private:
303 : amrex::Vector<amrex::Real> dt; ///< Timesteps for each level of refinement
304 : amrex::Vector<int> nsubsteps; ///< how many substeps on each level?
305 : int max_plot_level = -1;
306 :
307 : amrex::Vector<amrex::Real> t_old;///< Keep track of current old simulation time on each level
308 : int max_step = std::numeric_limits<int>::max(); ///< Maximum allowable timestep
309 : amrex::Real tstart = 0; ///< Default start time (default: 0)
310 : amrex::Real stop_time = NAN; ///< Default stop time
311 :
312 : protected:
313 : bool integrate_variables_before_advance = true;
314 : bool integrate_variables_after_advance = false;
315 :
316 : protected:
317 : struct {
318 : int number_of_fabs = 0;
319 : std::vector<Set::Field<Set::Scalar>*> fab_array;
320 : std::vector<int> ncomp_array;
321 : std::vector<int> nghost_array;
322 : std::vector<std::string> name_array;
323 : std::vector<BC::BC<Set::Scalar>*> physbc_array;
324 : std::vector<bool> writeout_array;
325 : bool any = true;
326 : bool all = false;
327 : } node;
328 :
329 : struct {
330 : int number_of_fabs = 0;
331 : std::vector<Set::Field<Set::Scalar>*> fab_array;
332 : std::vector<int> ncomp_array;
333 : std::vector<int> nghost_array;
334 : std::vector<std::string> name_array;
335 : std::vector<BC::BC<Set::Scalar>*> physbc_array;
336 : std::vector<bool> writeout_array;
337 : bool any = true;
338 : bool all = false;
339 : } cell;
340 :
341 : std::vector<BaseField*> m_basefields;
342 : std::vector<BaseField*> m_basefields_cell;
343 :
344 : BC::Nothing bcnothing;
345 :
346 : // KEEP TRACK OF ALL INTEGRATED VARIABLES
347 : struct {
348 : int interval = -1;
349 : Set::Scalar dt = NAN;
350 : int plot_int = -1;
351 : Set::Scalar plot_dt = NAN;
352 : int number = 0;
353 : std::vector<Set::Scalar*> vars;
354 : std::vector<std::string> names;
355 : std::vector<bool> extensives;
356 : } thermo;
357 :
358 : // REGRIDDING
359 : int regrid_int = -1; ///< Determine how often to regrid (default: 2)
360 : int base_regrid_int = -1; ///< Determine how often to regrid based on coarse level only (default: 0)
361 :
362 : std::string restart_file_cell = "";
363 : std::string restart_file_node = "";
364 :
365 : struct {
366 : int on = 0;
367 : std::vector<amrex::Box> box;
368 : } explicitmesh;
369 :
370 : protected:
371 : int plot_int = -1; ///< How frequently to dump plot file (default: never)
372 : Set::Scalar plot_dt = -1.0;
373 :
374 : int abort_on_nan = true;
375 : };
376 :
377 :
378 : template<>
379 : ALAMO_SINGLE_DEFINITION
380 40 : void Integrator::AddField<Set::Scalar, Set::Hypercube::Cell>(Set::Field<Set::Scalar>& new_field,
381 : BC::BC<Set::Scalar>* new_bc,
382 : int ncomp,
383 : int nghost,
384 : std::string name,
385 : bool writeout,
386 : bool /*evolving*/)
387 : {
388 40 : int nlevs_max = maxLevel() + 1;
389 40 : new_field.resize(nlevs_max);
390 40 : cell.fab_array.push_back(&new_field);
391 40 : if (new_bc != nullptr) cell.physbc_array.push_back(new_bc);
392 0 : else cell.physbc_array.push_back(&bcnothing);
393 40 : cell.ncomp_array.push_back(ncomp);
394 40 : cell.nghost_array.push_back(nghost);
395 40 : cell.name_array.push_back(name);
396 40 : cell.writeout_array.push_back(writeout);
397 40 : cell.number_of_fabs++;
398 40 : }
399 :
400 : template<>
401 : ALAMO_SINGLE_DEFINITION
402 21 : void Integrator::AddField<Set::Scalar, Set::Hypercube::Node>(Set::Field<Set::Scalar>& new_field,
403 : BC::BC<Set::Scalar>* new_bc,
404 : int ncomp,
405 : int nghost,
406 : std::string name,
407 : bool writeout,
408 : bool /*evolving*/)
409 : {
410 : BL_PROFILE("Integrator::RegisterNodalFab");
411 42 : Util::Assert(INFO, TEST(new_bc == nullptr));
412 21 : int nlevs_max = maxLevel() + 1;
413 21 : new_field.resize(nlevs_max);
414 21 : node.fab_array.push_back(&new_field);
415 21 : node.physbc_array.push_back(&bcnothing);
416 21 : node.ncomp_array.push_back(ncomp);
417 21 : node.nghost_array.push_back(nghost);
418 21 : node.name_array.push_back(name);
419 21 : node.writeout_array.push_back(writeout);
420 21 : node.number_of_fabs++;
421 21 : }
422 :
423 : template<class T, int d>
424 : ALAMO_SINGLE_DEFINITION
425 98 : void Integrator::AddField(Set::Field<T>& new_field, BC::BC<T>* new_bc, int ncomp, int nghost, std::string name, bool writeout, bool evolving)
426 : {
427 : if (d == Set::Hypercube::Node)
428 : {
429 196 : Util::Assert(INFO, TEST(new_bc == nullptr));
430 98 : int nlevs_max = maxLevel() + 1;
431 98 : new_field.resize(nlevs_max);
432 98 : m_basefields.push_back(new Field<T>(new_field, geom, refRatio(), ncomp, nghost));
433 98 : m_basefields.back()->evolving = evolving;
434 98 : m_basefields.back()->writeout = writeout;
435 98 : m_basefields.back()->setName(name);
436 98 : m_basefields.back()->evolving = evolving;
437 98 : m_basefields.back()->m_gridtype = Set::Hypercube::Node;
438 : }
439 : else if (d == Set::Hypercube::Cell)
440 : {
441 : int nlevs_max = maxLevel() + 1;
442 : new_field.resize(nlevs_max);
443 : m_basefields_cell.push_back(new Field<T>(new_field, geom, refRatio(), ncomp, nghost));
444 : m_basefields_cell.back()->evolving = evolving;
445 : m_basefields_cell.back()->writeout = writeout;
446 : m_basefields_cell.back()->setName(name);
447 : m_basefields_cell.back()->evolving = evolving;
448 : if (new_bc) m_basefields_cell.back()->setBC(new_bc);
449 : else m_basefields_cell.back()->setBC(&bcnothing);
450 : m_basefields_cell.back()->m_gridtype = Set::Hypercube::Cell;
451 : }
452 : else
453 : {
454 : Util::Abort(INFO, "Only node and cell based fields can be added at this time");
455 : }
456 98 : }
457 :
458 :
459 :
460 : template<class T>
461 : ALAMO_SINGLE_DEFINITION
462 19 : void Integrator::RegisterGeneralFab(Set::Field<T>& new_fab, int ncomp, int nghost, bool evolving)
463 : {
464 : //Util::Warning(INFO, "RegisterGeneralFab is depricated. Please replace with AddField");
465 19 : AddField<T, Set::Hypercube::Node>(new_fab, nullptr, ncomp, nghost, "", true, evolving);
466 19 : }
467 : template<class T>
468 : ALAMO_SINGLE_DEFINITION
469 3 : void Integrator::RegisterGeneralFab(Set::Field<T>& new_fab, int ncomp, int nghost, std::string a_name, bool evolving)
470 : {
471 : //Util::Warning(INFO, "RegisterGeneralFab is depricated. Please replace with AddField");
472 3 : AddField<T, Set::Hypercube::Node>(new_fab, nullptr, ncomp, nghost, a_name, true, evolving);
473 3 : }
474 : template<class T>
475 : AMREX_ATTRIBUTE_WEAK
476 76 : void Integrator::RegisterGeneralFab(Set::Field<T>& new_fab, int ncomp, int nghost, bool writeout, std::string a_name, bool evolving)
477 : {
478 : //Util::Warning(INFO, "RegisterGeneralFab is depricated. Please replace with AddField");
479 76 : AddField<T, Set::Hypercube::Node>(new_fab, nullptr, ncomp, nghost, a_name, writeout, evolving);
480 76 : }
481 :
482 : }
483 : #endif
|
