-
Notifications
You must be signed in to change notification settings - Fork 12
/
simulation.hpp
2680 lines (2277 loc) · 99.4 KB
/
simulation.hpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#ifndef SIMULATION_HPP_ // NOLINT
#define SIMULATION_HPP_
//==============================================================================
// TwoMomentRad - a radiation transport library for patch-based AMR codes
// Copyright 2020 Benjamin Wibking.
// Released under the MIT license. See LICENSE file included in the GitHub repo.
//==============================================================================
/// \file simulation.cpp
/// \brief Implements classes and functions to organise the overall setup,
/// timestepping, solving, and I/O of a simulation.
// c++ headers
#include <cassert>
#include <cmath>
#include <csignal>
#include <cstdint>
#include <cstdio>
#if __has_include(<filesystem>)
#include <filesystem>
#elif __has_include(<experimental/filesystem>)
#include <experimental/filesystem>
namespace std
{
namespace filesystem = experimental::filesystem;
}
#endif
#include <fstream>
#include <iomanip>
#include <iostream>
#include <limits>
#include <memory>
#include <optional>
#include <ostream>
#include <stdexcept>
#include <variant>
// library headers
#include "AMReX.H"
#include "AMReX_AmrCore.H"
#include "AMReX_Array.H"
#include "AMReX_Array4.H"
#include "AMReX_AsyncOut.H"
#include "AMReX_BCRec.H"
#include "AMReX_BLassert.H"
#include "AMReX_DistributionMapping.H"
#include "AMReX_Extension.H"
#include "AMReX_FArrayBox.H"
#include "AMReX_FillPatchUtil.H"
#include "AMReX_FillPatcher.H"
#include "AMReX_GpuQualifiers.H"
#include "AMReX_INT.H"
#include "AMReX_IndexType.H"
#include "AMReX_IntVect.H"
#include "AMReX_Interpolater.H"
#include "AMReX_MultiFabUtil.H"
#include "AMReX_ParallelDescriptor.H"
#include "AMReX_REAL.H"
#include "AMReX_SPACE.H"
#include "AMReX_Vector.H"
#include "AMReX_VisMF.H"
#include "AMReX_YAFluxRegister.H"
#include <AMReX_Geometry.H>
#include <AMReX_MultiFab.H>
#include <AMReX_ParmParse.H>
#include <AMReX_PlotFileUtil.H>
#include <AMReX_Print.H>
#include <AMReX_Utility.H>
#include <fmt/core.h>
#include <yaml-cpp/yaml.h>
#ifdef AMREX_PARTICLES
#include "CICParticles.hpp"
#include <AMReX_AmrParticles.H>
#include <AMReX_Particles.H>
#endif
#if AMREX_SPACEDIM == 3
#include "AMReX_OpenBC.H"
#endif
#ifdef AMREX_USE_ASCENT
#include <AMReX_Conduit_Blueprint.H>
#include <ascent.hpp>
#endif
// internal headers
#include "fundamental_constants.H"
#include "grid.hpp"
#include "physics_info.hpp"
#include "physics_numVars.hpp"
#ifdef QUOKKA_USE_OPENPMD
#include "openPMD.hpp"
#endif
#define USE_YAFLUXREGISTER
#ifdef AMREX_USE_ASCENT
using namespace conduit;
using namespace ascent;
#endif
enum class ParticleStep { BeforePoissonSolve, AfterPoissonSolve };
using variant_t = std::variant<amrex::Real, std::string>;
namespace YAML
{
template <typename T> struct as_if<T, std::optional<T>> {
explicit as_if(const Node &node_) : node(node_) {}
const Node &node; // NOLINT(cppcoreguidelines-avoid-const-or-ref-data-members)
auto operator()() const -> std::optional<T>
{
std::optional<T> val;
T t;
if ((node.m_pNode != nullptr) && convert<T>::decode(node, t)) {
val = std::move(t);
}
return val;
}
};
// There is already a std::string partial specialisation,
// so we need a full specialisation here
template <> struct as_if<std::string, std::optional<std::string>> {
explicit as_if(const Node &node_) : node(node_) {}
const Node &node; // NOLINT(cppcoreguidelines-avoid-const-or-ref-data-members)
auto operator()() const -> std::optional<std::string>
{
std::optional<std::string> val;
std::string t;
if ((node.m_pNode != nullptr) && convert<std::string>::decode(node, t)) {
val = std::move(t);
}
return val;
}
};
} // namespace YAML
enum class FillPatchType { fillpatch_class, fillpatch_function };
// Main simulation class; solvers should inherit from this
template <typename problem_t> class AMRSimulation : public amrex::AmrCore
{
public:
amrex::Real maxDt_ = std::numeric_limits<double>::max(); // no limit by default
amrex::Real initDt_ = std::numeric_limits<double>::max(); // no limit by default
amrex::Real constantDt_ = 0.0;
amrex::Vector<int> istep; // which step?
amrex::Vector<int> nsubsteps; // how many substeps on each level?
amrex::Vector<amrex::Real> tNew_; // for state_new_cc_
amrex::Vector<amrex::Real> tOld_; // for state_old_cc_
amrex::Vector<amrex::Real> dt_; // timestep for each level
amrex::Vector<int> reductionFactor_; // timestep reduction factor for each level
amrex::Real stopTime_ = 1.0; // default
amrex::Real cflNumber_ = 0.3; // default
amrex::Real dtToleranceFactor_ = 1.1; // default
amrex::Long cycleCount_ = 0;
amrex::Long maxTimesteps_ = 1e4; // default
amrex::Long maxWalltime_ = 0; // default: no limit
int ascentInterval_ = -1; // -1 == no in-situ renders with Ascent
int plotfileInterval_ = -1; // -1 == no output
int projectionInterval_ = -1; // -1 == no output
int statisticsInterval_ = -1; // -1 == no output
amrex::Real plotTimeInterval_ = -1.0; // time interval for plt file
amrex::Real checkpointTimeInterval_ = -1.0; // time interval for checkpoints
int checkpointInterval_ = -1; // -1 == no output
int amrInterpMethod_ = 1; // 0 == piecewise constant, 1 == lincc_interp
amrex::Real reltolPoisson_ = 1.0e-5; // default
amrex::Real abstolPoisson_ = 1.0e-5; // default (scaled by minimum RHS value)
int doPoissonSolve_ = 0; // 1 == self-gravity enabled, 0 == disabled
amrex::Vector<amrex::MultiFab> phi;
amrex::Real densityFloor_ = 0.0; // default
amrex::Real tempCeiling_ = std::numeric_limits<double>::max(); // default
amrex::Real tempFloor_ = 0.0; // default
amrex::Real speedCeiling_ = std::numeric_limits<double>::max(); // default
std::unordered_map<std::string, variant_t> simulationMetadata_;
// constructor
explicit AMRSimulation(amrex::Vector<amrex::BCRec> &BCs_cc, amrex::Vector<amrex::BCRec> &BCs_fc) : BCs_cc_(BCs_cc), BCs_fc_(BCs_fc) { initialize(); }
explicit AMRSimulation(amrex::Vector<amrex::BCRec> &BCs_cc) : BCs_cc_(BCs_cc), BCs_fc_(builtin_BCs_fc(BCs_cc)) { initialize(); }
inline auto builtin_BCs_fc(amrex::Vector<amrex::BCRec> &BCs_cc) -> amrex::Vector<amrex::BCRec>
{
amrex::Vector<amrex::BCRec> BCs_fc(Physics_Indices<problem_t>::nvarPerDim_fc);
if (Physics_Traits<problem_t>::is_hydro_enabled) {
AMREX_ALWAYS_ASSERT(Physics_Indices<problem_t>::nvarPerDim_fc == 1);
// set boundary conditions for face velocities (used ONLY for tracer particles)
for (int i = 0; i < AMREX_SPACEDIM; ++i) {
// lower boundary
if (BCs_cc[Physics_Indices<problem_t>::hydroFirstIndex].lo(i) == amrex::BCType::int_dir) {
BCs_fc[Physics_Indices<problem_t>::velFirstIndex].setLo(i, amrex::BCType::int_dir);
} else {
BCs_fc[Physics_Indices<problem_t>::velFirstIndex].setLo(i, amrex::BCType::foextrap);
}
// upper boundary
if (BCs_cc[Physics_Indices<problem_t>::hydroFirstIndex].hi(i) == amrex::BCType::int_dir) {
BCs_fc[Physics_Indices<problem_t>::velFirstIndex].setHi(i, amrex::BCType::int_dir);
} else {
BCs_fc[Physics_Indices<problem_t>::velFirstIndex].setHi(i, amrex::BCType::foextrap);
}
}
}
return BCs_fc;
}
void initialize();
void PerformanceHints();
void readParameters();
void setInitialConditions();
void setInitialConditionsAtLevel_cc(int level, amrex::Real time);
void setInitialConditionsAtLevel_fc(int level, amrex::Real time);
void evolve();
void computeTimestep();
auto computeTimestepAtLevel(int lev) -> amrex::Real;
void AverageFCToCC(amrex::MultiFab &mf_cc, const amrex::MultiFab &mf_fc, int idim, int dstcomp_start, int srccomp_start, int srccomp_total,
int nGrow) const;
virtual void computeMaxSignalLocal(int level) = 0;
virtual auto computeExtraPhysicsTimestep(int lev) -> amrex::Real = 0;
virtual void advanceSingleTimestepAtLevel(int lev, amrex::Real time, amrex::Real dt_lev, int ncycle) = 0;
virtual void preCalculateInitialConditions() = 0;
virtual void setInitialConditionsOnGrid(quokka::grid grid_elem) = 0;
virtual void setInitialConditionsOnGridFaceVars(quokka::grid grid_elem) = 0;
virtual void createInitialParticles() = 0;
virtual void computeAfterTimestep() = 0;
virtual void computeAfterEvolve(amrex::Vector<amrex::Real> &initSumCons) = 0;
virtual void fillPoissonRhsAtLevel(amrex::MultiFab &rhs, int lev) = 0;
virtual void applyPoissonGravityAtLevel(amrex::MultiFab const &phi, int lev, amrex::Real dt) = 0;
// compute derived variables
virtual void ComputeDerivedVar(int lev, std::string const &dname, amrex::MultiFab &mf, int ncomp) const = 0;
// compute projected vars
[[nodiscard]] virtual auto ComputeProjections(int dir) const -> std::unordered_map<std::string, amrex::BaseFab<amrex::Real>> = 0;
// compute statistics
virtual auto ComputeStatistics() -> std::map<std::string, amrex::Real> = 0;
// fix-up any unphysical states created by AMR operations
// (e.g., caused by the flux register or from interpolation)
virtual void FixupState(int level) = 0;
// tag cells for refinement
void ErrorEst(int lev, amrex::TagBoxArray &tags, amrex::Real time, int ngrow) override = 0;
// Make a new level using provided BoxArray and DistributionMapping
void MakeNewLevelFromCoarse(int lev, amrex::Real time, const amrex::BoxArray &ba, const amrex::DistributionMapping &dm) override;
// Remake an existing level using provided BoxArray and DistributionMapping
void RemakeLevel(int lev, amrex::Real time, const amrex::BoxArray &ba, const amrex::DistributionMapping &dm) override;
// Delete level data
void ClearLevel(int lev) override;
// Make a new level from scratch using provided BoxArray and
// DistributionMapping
void MakeNewLevelFromScratch(int lev, amrex::Real time, const amrex::BoxArray &ba, const amrex::DistributionMapping &dm) override;
// AMR utility functions
template <typename PreInterpHook, typename PostInterpHook>
void fillBoundaryConditions(amrex::MultiFab &S_filled, amrex::MultiFab &state, int lev, amrex::Real time, quokka::centering cen, quokka::direction dir,
PreInterpHook const &pre_interp, PostInterpHook const &post_interp, FillPatchType fptype = FillPatchType::fillpatch_class);
template <typename PreInterpHook, typename PostInterpHook>
void FillPatchWithData(int lev, amrex::Real time, amrex::MultiFab &mf, amrex::Vector<amrex::MultiFab *> &coarseData,
amrex::Vector<amrex::Real> &coarseTime, amrex::Vector<amrex::MultiFab *> &fineData, amrex::Vector<amrex::Real> &fineTime,
int icomp, int ncomp, amrex::Vector<amrex::BCRec> &BCs, quokka::centering &cen, FillPatchType fptype,
PreInterpHook const &pre_interp, PostInterpHook const &post_interp);
static void InterpHookNone(amrex::MultiFab &mf, int scomp, int ncomp);
virtual void FillPatch(int lev, amrex::Real time, amrex::MultiFab &mf, int icomp, int ncomp, quokka::centering cen, quokka::direction dir,
FillPatchType fptype);
auto getAmrInterpolaterCellCentered() -> amrex::MFInterpolater *;
auto getAmrInterpolaterFaceCentered() -> amrex::Interpolater *;
void FillCoarsePatch(int lev, amrex::Real time, amrex::MultiFab &mf, int icomp, int ncomp, amrex::Vector<amrex::BCRec> &BCs, quokka::centering cen,
quokka::direction dir);
void GetData(int lev, amrex::Real time, amrex::Vector<amrex::MultiFab *> &data, amrex::Vector<amrex::Real> &datatime, quokka::centering cen,
quokka::direction dir);
void AverageDown();
void AverageDownTo(int crse_lev);
void timeStepWithSubcycling(int lev, amrex::Real time, int iteration);
void calculateGpotAllLevels();
void gravAccelAllLevels(amrex::Real dt);
void ellipticSolveAllLevels(amrex::Real dt);
void incrementFluxRegisters(amrex::MFIter &mfi, amrex::YAFluxRegister *fr_as_crse, amrex::YAFluxRegister *fr_as_fine,
std::array<amrex::FArrayBox, AMREX_SPACEDIM> &fluxArrays, int lev, amrex::Real dt_lev);
void incrementFluxRegisters(amrex::YAFluxRegister *fr_as_crse, amrex::YAFluxRegister *fr_as_fine,
std::array<amrex::MultiFab, AMREX_SPACEDIM> &fluxArrays, int lev, amrex::Real dt_lev);
// boundary condition
AMREX_GPU_DEVICE static void setCustomBoundaryConditions(const amrex::IntVect &iv, amrex::Array4<amrex::Real> const &dest, int dcomp, int numcomp,
amrex::GeometryData const &geom, amrex::Real time, const amrex::BCRec *bcr, int bcomp,
int orig_comp); // template specialized by problem generator
// boundary condition
AMREX_GPU_DEVICE static void setCustomBoundaryConditionsFaceVar(const amrex::IntVect &iv, amrex::Array4<amrex::Real> const &dest, int dcomp,
int numcomp, amrex::GeometryData const &geom, amrex::Real time, const amrex::BCRec *bcr,
int bcomp,
int orig_comp); // template specialized by problem generator
template <typename ReduceOp, typename F> auto computePlaneProjection(F const &user_f, int dir) const -> amrex::BaseFab<amrex::Real>;
// I/O functions
[[nodiscard]] auto PlotFileName(int lev) const -> std::string;
[[nodiscard]] auto CustomPlotFileName(const char *base, int lev) const -> std::string;
[[nodiscard]] auto GetPlotfileVarNames() const -> amrex::Vector<std::string>;
[[nodiscard]] auto PlotFileMF(int included_ghosts) -> amrex::Vector<amrex::MultiFab>;
[[nodiscard]] auto PlotFileMFAtLevel(int lev, int included_ghosts) -> amrex::MultiFab;
void WriteMetadataFile(std::string const &MetadataFileName) const;
void ReadMetadataFile(std::string const &chkfilename);
void WriteStatisticsFile();
void WritePlotFile();
void WriteProjectionPlotfile() const;
void WriteCheckpointFile() const;
void SetLastCheckpointSymlink(std::string const &checkpointname) const;
void ReadCheckpointFile();
auto getWalltime() -> amrex::Real;
void setChkFile(std::string const &chkfile_number);
[[nodiscard]] auto getOldMF_fc() const -> amrex::Vector<amrex::Array<amrex::MultiFab, AMREX_SPACEDIM>> const &;
[[nodiscard]] auto getNewMF_fc() const -> amrex::Vector<amrex::Array<amrex::MultiFab, AMREX_SPACEDIM>> const &;
// particle functions
void kickParticlesAllLevels(amrex::Real dt);
void driftParticlesAllLevels(amrex::Real dt);
#ifdef AMREX_USE_ASCENT
void AscentCustomActions(conduit::Node const &blueprintMesh);
void RenderAscent();
#endif
protected:
amrex::Vector<amrex::BCRec> BCs_cc_; // on level 0
amrex::Vector<amrex::BCRec> BCs_fc_; // on level 0
amrex::Vector<amrex::MultiFab> state_old_cc_;
amrex::Vector<amrex::MultiFab> state_new_cc_;
amrex::Vector<amrex::Array<amrex::MultiFab, AMREX_SPACEDIM>> state_old_fc_;
amrex::Vector<amrex::Array<amrex::MultiFab, AMREX_SPACEDIM>> state_new_fc_;
amrex::Vector<amrex::MultiFab> max_signal_speed_; // needed to compute CFL timestep
// flux registers: store fluxes at coarse-fine interface for synchronization
// this will be sized "nlevs_max+1"
// NOTE: the flux register associated with flux_reg[lev] is associated with
// the lev/lev-1 interface (and has grid spacing associated with lev-1)
// therefore flux_reg[0] and flux_reg[nlevs_max] are never actually used in
// the reflux operation
amrex::Vector<std::unique_ptr<amrex::YAFluxRegister>> flux_reg_;
// This is for fillpatch during timestepping, but not for regridding.
amrex::Vector<std::unique_ptr<amrex::FillPatcher<amrex::MultiFab>>> fillpatcher_;
// Nghost = number of ghost cells for each array
int nghost_cc_ = 4; // PPM needs nghost >= 3, PPM+flattening needs nghost >= 4
int nghost_fc_ = Physics_Traits<problem_t>::is_mhd_enabled ? 4 : 2; // 4 needed for MHD, otherwise only 2 for tracer particles
amrex::Vector<std::string> componentNames_cc_;
amrex::Vector<std::string> componentNames_fc_;
amrex::Vector<std::string> derivedNames_;
bool areInitialConditionsDefined_ = false;
/// output parameters
std::string plot_file{"plt"}; // plotfile prefix
std::string chk_file{"chk"}; // checkpoint prefix
std::string stats_file{"history.txt"}; // statistics filename
/// input parameters (if >= 0 we restart from a checkpoint)
std::string restart_chkfile;
/// AMR-specific parameters
int regrid_int = 2; // regrid interval (number of coarse steps)
int do_reflux = 1; // 1 == reflux, 0 == no reflux
int do_subcycle = 1; // 1 == subcycle, 0 == no subcyle
int suppress_output = 0; // 1 == show timestepping, 0 == do not output each timestep
// performance metrics
amrex::Long cellUpdates_ = 0;
amrex::Vector<amrex::Long> cellUpdatesEachLevel_;
// gravity
amrex::Real Gconst_ = C::Gconst; // gravitational constant G
// tracer particles
#ifdef AMREX_PARTICLES
void InitParticles(); // create tracer particles
void InitCICParticles(); // create CIC particles
int do_tracers = 0;
int do_cic_particles = 0;
std::unique_ptr<amrex::AmrTracerParticleContainer> TracerPC;
std::unique_ptr<quokka::CICParticleContainer> CICParticles;
#endif
// external objects
#ifdef AMREX_USE_ASCENT
Ascent ascent_;
#endif
};
template <typename problem_t> void AMRSimulation<problem_t>::setChkFile(std::string const &chkfile_number) { restart_chkfile = chkfile_number; }
template <typename problem_t> auto AMRSimulation<problem_t>::getOldMF_fc() const -> const amrex::Vector<amrex::Array<amrex::MultiFab, AMREX_SPACEDIM>> &
{
return state_old_fc_;
}
template <typename problem_t> auto AMRSimulation<problem_t>::getNewMF_fc() const -> const amrex::Vector<amrex::Array<amrex::MultiFab, AMREX_SPACEDIM>> &
{
return state_new_fc_;
}
template <typename problem_t> void AMRSimulation<problem_t>::initialize()
{
BL_PROFILE("AMRSimulation::initialize()");
readParameters();
// print derived vars
if (!derivedNames_.empty()) {
amrex::Print() << "Using derived variables:\n";
for (auto const &name : derivedNames_) {
amrex::Print() << "\t" << name << "\n";
}
amrex::Print() << "\n";
}
int nlevs_max = max_level + 1;
istep.resize(nlevs_max, 0);
nsubsteps.resize(nlevs_max, 1);
if (do_subcycle == 1) {
for (int lev = 1; lev <= max_level; ++lev) {
nsubsteps[lev] = MaxRefRatio(lev - 1);
}
}
tNew_.resize(nlevs_max, 0.0);
tOld_.resize(nlevs_max, -1.e100);
dt_.resize(nlevs_max, 1.e100);
reductionFactor_.resize(nlevs_max, 1);
state_new_cc_.resize(nlevs_max);
state_old_cc_.resize(nlevs_max);
if constexpr (Physics_Indices<problem_t>::nvarTotal_fc > 0) {
state_new_fc_.resize(nlevs_max);
state_old_fc_.resize(nlevs_max);
}
max_signal_speed_.resize(nlevs_max);
flux_reg_.resize(nlevs_max + 1);
fillpatcher_.resize(nlevs_max + 1);
cellUpdatesEachLevel_.resize(nlevs_max, 0);
// check that grids will be properly nested on each level
// (this is necessary since FillPatch only fills from non-ghost cells on
// lev-1)
auto checkIsProperlyNested = [this](int const lev, amrex::IntVect const &blockingFactor) {
return amrex::ProperlyNested(refRatio(lev - 1), blockingFactor, nghost_cc_, amrex::IndexType::TheCellType(), &amrex::cell_cons_interp);
};
for (int lev = 1; lev <= max_level; ++lev) {
if (!checkIsProperlyNested(lev, blocking_factor[lev])) {
// level lev is not properly nested
amrex::Print() << "Blocking factor is too small for proper grid nesting! "
"Increase blocking factor to >= ceil(nghost,ref_ratio)*ref_ratio."
<< std::endl; // NOLINT(performance-avoid-endl)
amrex::Abort("Grids not properly nested!");
}
}
#ifdef AMREX_USE_ASCENT
// initialize Ascent
conduit::Node ascent_options;
ascent_options["mpi_comm"] = MPI_Comm_c2f(amrex::ParallelContext::CommunicatorSub());
ascent_.open(ascent_options);
#endif
}
template <typename problem_t> void AMRSimulation<problem_t>::PerformanceHints()
{
// Check requested MPI ranks and available boxes
for (int ilev = 0; ilev <= finestLevel(); ++ilev) {
const amrex::Long nboxes = boxArray(ilev).size();
if (amrex::ParallelDescriptor::NProcs() > nboxes) {
amrex::Print() << "\n[Warning] [Performance] Too many resources / too little work!\n"
<< " It looks like you requested more compute resources than "
<< " the number of boxes of cells available on level " << ilev << " (" << nboxes << "). "
<< "You started with (" << amrex::ParallelDescriptor::NProcs() << ") MPI ranks, so ("
<< amrex::ParallelDescriptor::NProcs() - nboxes << ") rank(s) will have no work on this level.\n"
#ifdef AMREX_USE_GPU
<< " On GPUs, consider using 1-8 boxes per GPU per level that "
"together fill each GPU's memory sufficiently.\n"
#endif
<< "\n";
}
}
// check that blocking_factor and max_grid_size are set to reasonable values
#ifdef AMREX_USE_GPU
const int recommended_blocking_factor = 32;
const int recommended_max_grid_size = 128;
#else
const int recommended_blocking_factor = 16;
const int recommended_max_grid_size = 64;
#endif
int min_blocking_factor = INT_MAX;
int min_max_grid_size = INT_MAX;
for (int ilev = 0; ilev <= finestLevel(); ++ilev) {
min_blocking_factor = std::min(min_blocking_factor, blocking_factor[ilev].min());
min_max_grid_size = std::min(min_max_grid_size, max_grid_size[ilev].min());
}
if (min_blocking_factor < recommended_blocking_factor) {
amrex::Print() << "\n[Warning] [Performance] The grid blocking factor (" << min_blocking_factor
<< ") is too small for reasonable performance. It should be 32 (or "
"greater) when running on GPUs, and 16 (or greater) when running on "
"CPUs.\n";
}
if (min_max_grid_size < recommended_max_grid_size) {
amrex::Print() << "\n[Warning] [Performance] The maximum grid size (" << min_max_grid_size
<< ") is too small for reasonable performance. It should be "
"128 (or greater) when running on GPUs, and 64 (or "
"greater) when running on CPUs.\n";
}
#ifdef QUOKKA_USE_OPENPMD
// warning about face-centered variables and OpenPMD outputs
if constexpr (Physics_Indices<problem_t>::nvarTotal_fc > 0) {
amrex::Print() << "\n[Warning] [I/O] Plotfiles will ONLY contain cell-centered averages of face-centered variables!"
<< " Support for outputting face-centered variables for openPMD is not yet implemented.\n";
}
#endif
}
template <typename problem_t> void AMRSimulation<problem_t>::readParameters()
{
BL_PROFILE("AMRSimulation::readParameters()");
// ParmParse reads inputs from the *.inputs file
amrex::ParmParse pp;
// Default nsteps = 1e4
pp.query("max_timesteps", maxTimesteps_);
// Default CFL number == 0.3, set to whatever is in the file
pp.query("cfl", cflNumber_);
// Default AMR interpolation method == lincc_interp
pp.query("amr_interpolation_method", amrInterpMethod_);
// Default stopping time
pp.query("stop_time", stopTime_);
// Default ascent render interval
pp.query("ascent_interval", ascentInterval_);
// Default output interval
pp.query("plotfile_interval", plotfileInterval_);
// Default projection interval
pp.query("projection_interval", projectionInterval_);
// Default statistics interval
pp.query("statistics_interval", statisticsInterval_);
// Default Time interval
pp.query("plottime_interval", plotTimeInterval_);
// Default Time interval
pp.query("checkpointtime_interval", checkpointTimeInterval_);
// Default checkpoint interval
pp.query("checkpoint_interval", checkpointInterval_);
// Default do_reflux = 1
pp.query("do_reflux", do_reflux);
// Default do_subcycle = 1
pp.query("do_subcycle", do_subcycle);
// Default do_tracers = 0 (turns on/off tracer particles)
pp.query("do_tracers", do_tracers);
// Default do_cic_particles = 0 (turns on/off CIC particles)
pp.query("do_cic_particles", do_cic_particles);
// Default suppress_output = 0
pp.query("suppress_output", suppress_output);
// specify this on the command-line in order to restart from a checkpoint
// file
pp.query("restartfile", restart_chkfile);
// Specify derived variables to save to plotfiles
pp.queryarr("derived_vars", derivedNames_);
// re-grid interval
pp.query("regrid_interval", regrid_int);
// read density floor in g cm^-3
pp.query("density_floor", densityFloor_);
// read temperature floor in K
pp.query("temperature_floor", tempFloor_);
// read temperature ceiling in K
pp.query("temperature_ceiling", tempCeiling_);
// read speed ceiling in cm s^-1
pp.query("speed_ceiling", speedCeiling_);
// specify maximum walltime in HH:MM:SS format
std::string maxWalltimeInput;
pp.query("max_walltime", maxWalltimeInput);
// convert to seconds
int hours = 0;
int minutes = 0;
int seconds = 0;
int nargs = std::sscanf(maxWalltimeInput.c_str(), "%d:%d:%d", &hours, &minutes, &seconds); // NOLINT
if (nargs == 3) {
maxWalltime_ = 3600 * hours + 60 * minutes + seconds;
amrex::Print() << fmt::format("Setting walltime limit to {} hours, {} minutes, {} seconds.\n", hours, minutes, seconds);
}
// set gravity runtime parameters
{
const amrex::ParmParse hpp("gravity");
hpp.query("Gconst", Gconst_);
}
}
template <typename problem_t> void AMRSimulation<problem_t>::setInitialConditions()
{
BL_PROFILE("AMRSimulation::setInitialConditions()");
if (restart_chkfile.empty()) {
// start simulation from the beginning
const amrex::Real time = 0.0;
InitFromScratch(time);
AverageDown();
#ifdef AMREX_PARTICLES
if (do_tracers != 0) {
InitParticles();
}
if (do_cic_particles != 0) {
InitCICParticles();
}
#endif
if (checkpointInterval_ > 0) {
WriteCheckpointFile();
}
} else {
// restart from a checkpoint
ReadCheckpointFile();
}
calculateGpotAllLevels();
// abort if amrex.async_out=1, it is currently broken
if (amrex::AsyncOut::UseAsyncOut()) {
amrex::Print() << "[ERROR] [FATAL] AsyncOut is currently broken! If you want to "
"run with AsyncOut anyway (THIS MAY CAUSE DATA CORRUPTION), comment "
"out this line in src/simulation.hpp. Aborting."
<< std::endl; // NOLINT(performance-avoid-endl)
amrex::Abort();
}
#ifdef AMREX_USE_ASCENT
if (ascentInterval_ > 0) {
RenderAscent();
}
#endif
if (plotfileInterval_ > 0) {
WritePlotFile();
}
if (projectionInterval_ > 0) {
WriteProjectionPlotfile();
}
if (statisticsInterval_ > 0) {
WriteStatisticsFile();
}
// ensure that there are enough boxes per MPI rank
PerformanceHints();
}
template <typename problem_t> auto AMRSimulation<problem_t>::computeTimestepAtLevel(int lev) -> amrex::Real
{
// compute CFL timestep on level 'lev'
BL_PROFILE("AMRSimulation::computeTimestepAtLevel()");
// compute hydro timestep on level 'lev'
computeMaxSignalLocal(lev);
const amrex::Real domain_signal_max = max_signal_speed_[lev].norminf();
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const &dx = geom[lev].CellSizeArray();
const amrex::Real dx_min = std::min({AMREX_D_DECL(dx[0], dx[1], dx[2])});
const amrex::Real hydro_dt = cflNumber_ * (dx_min / domain_signal_max);
// compute timestep due to extra physics on level 'lev'
const amrex::Real extra_physics_dt = computeExtraPhysicsTimestep(lev);
// return minimum timestep
return std::min(hydro_dt, extra_physics_dt);
}
template <typename problem_t> void AMRSimulation<problem_t>::computeTimestep()
{
BL_PROFILE("AMRSimulation::computeTimestep()");
// compute candidate timestep dt_tmp on each level
amrex::Vector<amrex::Real> dt_tmp(finest_level + 1);
for (int level = 0; level <= finest_level; ++level) {
dt_tmp[level] = computeTimestepAtLevel(level);
}
// limit change in timestep on each level
constexpr amrex::Real change_max = 1.1;
for (int level = 0; level <= finest_level; ++level) {
dt_tmp[level] = std::min(dt_tmp[level], change_max * dt_[level]);
}
// set default subcycling pattern
if (do_subcycle == 1) {
for (int lev = 1; lev <= max_level; ++lev) {
nsubsteps[lev] = MaxRefRatio(lev - 1);
reductionFactor_[lev] = 1; // reset additional subcycling
}
}
// compute root level timestep given nsubsteps
amrex::Real dt_0 = dt_tmp[0];
amrex::Long n_factor = 1;
for (int level = 0; level <= finest_level; ++level) {
n_factor *= nsubsteps[level];
dt_0 = std::min(dt_0, static_cast<amrex::Real>(n_factor) * dt_tmp[level]);
dt_0 = std::min(dt_0, maxDt_); // limit to maxDt_
if (tNew_[level] == 0.0) { // first timestep
dt_0 = std::min(dt_0, initDt_);
}
if (constantDt_ > 0.0) { // use constant timestep if set
dt_0 = constantDt_;
}
}
// compute global timestep assuming no subcycling
amrex::Real dt_global = dt_tmp[0];
for (int level = 0; level <= finest_level; ++level) {
dt_global = std::min(dt_global, dt_tmp[level]);
dt_global = std::min(dt_global, maxDt_); // limit to maxDt_
if (tNew_[level] == 0.0) { // special case: first timestep
dt_global = std::min(dt_global, initDt_);
}
if (constantDt_ > 0.0) { // special case: constant timestep
dt_global = constantDt_;
}
}
// compute work estimate for subcycling
amrex::Long n_factor_work = 1;
amrex::Long work_subcycling = 0;
for (int level = 0; level <= finest_level; ++level) {
n_factor_work *= nsubsteps[level];
work_subcycling += n_factor_work * CountCells(level);
}
// compute work estimate for non-subcycling
amrex::Long total_cells = 0;
for (int level = 0; level <= finest_level; ++level) {
total_cells += CountCells(level);
}
const amrex::Real work_nonsubcycling = static_cast<amrex::Real>(total_cells) * (dt_0 / dt_global);
if (work_nonsubcycling <= static_cast<amrex::Real>(work_subcycling)) {
// use global timestep on this coarse step
if (verbose) {
const amrex::Real ratio = work_nonsubcycling / static_cast<amrex::Real>(work_subcycling);
amrex::Print() << "\t>> Using global timestep on this coarse step (estimated work ratio: " << ratio << ").\n";
}
for (int lev = 1; lev <= max_level; ++lev) {
nsubsteps[lev] = 1;
}
}
// Limit dt to avoid overshooting stop_time
const amrex::Real eps = 1.e-3 * dt_0;
if (tNew_[0] + dt_0 > stopTime_ - eps) {
dt_0 = stopTime_ - tNew_[0];
}
// assign timesteps on each level
dt_[0] = dt_0;
for (int level = 1; level <= finest_level; ++level) {
dt_[level] = dt_[level - 1] / nsubsteps[level];
}
}
template <typename problem_t> auto AMRSimulation<problem_t>::getWalltime() -> amrex::Real
{
const static amrex::Real start_time = amrex::ParallelDescriptor::second(); // initialized on first call
const amrex::Real time = amrex::ParallelDescriptor::second();
return time - start_time;
}
template <typename problem_t> void AMRSimulation<problem_t>::evolve()
{
BL_PROFILE("AMRSimulation::evolve()");
AMREX_ALWAYS_ASSERT(areInitialConditionsDefined_);
amrex::Real cur_time = tNew_[0];
#ifdef AMREX_USE_ASCENT
int last_ascent_step = 0;
#endif
int last_projection_step = 0;
int last_statistics_step = 0;
int last_plot_file_step = 0;
double next_plot_file_time = plotTimeInterval_;
double next_chk_file_time = checkpointTimeInterval_;
int last_chk_file_step = 0;
const int ncomp_cc = Physics_Indices<problem_t>::nvarTotal_cc;
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const &dx0 = geom[0].CellSizeArray();
amrex::Real const vol = AMREX_D_TERM(dx0[0], *dx0[1], *dx0[2]);
amrex::Vector<amrex::Real> init_sum_cons(ncomp_cc);
for (int n = 0; n < ncomp_cc; ++n) {
const int lev = 0;
init_sum_cons[n] = state_new_cc_[lev].sum(n) * vol;
}
getWalltime(); // initialize start_time
// Main time loop
for (int step = istep[0]; step < maxTimesteps_ && cur_time < stopTime_; ++step) {
if (suppress_output == 0) {
amrex::Print() << "\nCoarse STEP " << step + 1 << " at t = " << cur_time << " (" << (cur_time / stopTime_) * 100. << "%) starts ..."
<< '\n';
}
amrex::ParallelDescriptor::Barrier(); // synchronize all MPI ranks
computeTimestep();
// do particle leapfrog (first kick at time t)
kickParticlesAllLevels(dt_[0]);
// hyperbolic advance over all levels
// (N.B. when AMR is enabled, regridding may happen during this function!)
int lev = 0; // coarsest level
const int iteration = 1; // this is the first call to advance level 'lev'
timeStepWithSubcycling(lev, cur_time, iteration);
// drift particles from t to (t + dt)
// N.B.: MUST be done *before* Poisson solve at new time!
driftParticlesAllLevels(dt_[0]);
// elliptic solve over entire AMR grid (post-timestep)
ellipticSolveAllLevels(dt_[0]);
// do particle leapfrog (second kick at t + dt)
kickParticlesAllLevels(dt_[0]);
cur_time += dt_[0];
++cycleCount_;
computeAfterTimestep();
// sync up time (to avoid roundoff error)
for (lev = 0; lev <= finest_level; ++lev) {
AMREX_ALWAYS_ASSERT(std::abs((tNew_[lev] - cur_time) / cur_time) < 1e-10);
tNew_[lev] = cur_time;
}
#ifdef AMREX_USE_ASCENT
if (ascentInterval_ > 0 && (step + 1) % ascentInterval_ == 0) {
last_ascent_step = step + 1;
RenderAscent();
}
#endif
if (statisticsInterval_ > 0 && (step + 1) % statisticsInterval_ == 0) {
last_statistics_step = step + 1;
WriteStatisticsFile();
}
if (plotfileInterval_ > 0 && (step + 1) % plotfileInterval_ == 0) {
last_plot_file_step = step + 1;
WritePlotFile();
}
if (projectionInterval_ > 0 && (step + 1) % projectionInterval_ == 0) {
last_projection_step = step + 1;
WriteProjectionPlotfile();
}
// Writing Plot files at time intervals
if (plotTimeInterval_ > 0 && next_plot_file_time <= cur_time) {
next_plot_file_time += plotTimeInterval_;
WritePlotFile();
}
if (checkpointTimeInterval_ > 0 && next_chk_file_time <= cur_time) {
next_chk_file_time += checkpointTimeInterval_;
WriteCheckpointFile();
}
if (checkpointInterval_ > 0 && (step + 1) % checkpointInterval_ == 0) {
last_chk_file_step = step + 1;
WriteCheckpointFile();
}
if (cur_time >= stopTime_ - 1.e-6 * dt_[0]) {
// we have reached stopTime_
break;
}
if (maxWalltime_ > 0 && getWalltime() > 0.9 * maxWalltime_) {
// we have exceeded 90% of maxWalltime_
break;
}
}
amrex::Real elapsed_sec = getWalltime();
// compute reference solution (if it's a test problem)
computeAfterEvolve(init_sum_cons);
// compute conservation error
for (int n = 0; n < ncomp_cc; ++n) {
amrex::Real const final_sum = state_new_cc_[0].sum(n) * vol;
amrex::Real const abs_err = (final_sum - init_sum_cons[n]);
amrex::Print() << "Initial " << componentNames_cc_[n] << " = " << init_sum_cons[n] << '\n';
amrex::Print() << "\tabsolute conservation error = " << abs_err << '\n';
if (init_sum_cons[n] != 0.0) {
amrex::Real const rel_err = abs_err / init_sum_cons[n];
amrex::Print() << "\trelative conservation error = " << rel_err << '\n';
}
amrex::Print() << '\n';
}
// compute zone-cycles/sec
const int IOProc = amrex::ParallelDescriptor::IOProcessorNumber();
amrex::ParallelDescriptor::ReduceRealMax(elapsed_sec, IOProc);
const double microseconds_per_update = 1.0e6 * elapsed_sec / cellUpdates_;
const double megaupdates_per_second = 1.0 / microseconds_per_update;
amrex::Print() << "Performance figure-of-merit: " << microseconds_per_update << " μs/zone-update [" << megaupdates_per_second << " Mupdates/s]\n";
for (int lev = 0; lev <= max_level; ++lev) {
amrex::Print() << "Zone-updates on level " << lev << ": " << cellUpdatesEachLevel_[lev] << "\n";
}
amrex::Print() << '\n';
// write final checkpoint
if (checkpointInterval_ > 0 && istep[0] > last_chk_file_step) {
WriteCheckpointFile();
}
// write final plotfile
if (plotfileInterval_ > 0 && istep[0] > last_plot_file_step) {
WritePlotFile();
}
// write final projection
if (projectionInterval_ > 0 && istep[0] > last_projection_step) {
WriteProjectionPlotfile();
}
// write final statistics
if (statisticsInterval_ > 0 && istep[0] > last_statistics_step) {
WriteStatisticsFile();
}
#ifdef AMREX_USE_ASCENT
// close Ascent
ascent_.close();
#endif
}
template <typename problem_t> void AMRSimulation<problem_t>::calculateGpotAllLevels()
{
#if AMREX_SPACEDIM == 3
if (doPoissonSolve_) {
if (do_subcycle == 1) { // not supported
amrex::Abort("Poisson solve is not support when AMR subcycling is enabled! You must set do_subcycle = 0.");
}
BL_PROFILE_REGION("GravitySolver");
// set up elliptic solve object
amrex::OpenBCSolver poissonSolver(Geom(0, finest_level), boxArray(0, finest_level), DistributionMap(0, finest_level));
if (verbose) {
poissonSolver.setVerbose(true);