Block-Structured AMR Software Framework
AMReX_IntegratorBase.H
Go to the documentation of this file.
1 #ifndef AMREX_INTEGRATOR_BASE_H
2 #define AMREX_INTEGRATOR_BASE_H
3 #include <AMReX_Config.H>
4 #include <AMReX_REAL.H>
5 #include <AMReX_Vector.H>
6 #include <AMReX_MultiFab.H>
7 
8 #if defined(AMREX_PARTICLES)
9 #include <AMReX_Particles.H>
10 #endif
11 
12 #include <functional>
13 #include <type_traits>
14 
15 namespace amrex {
16 
17 template<class T, typename Tv = void> struct IntegratorOps;
18 
19 #if defined(AMREX_PARTICLES)
20 template<class T>
21 struct IntegratorOps<T, std::enable_if_t<std::is_base_of_v<amrex::ParticleContainerBase, T> > >
22 {
23 
24  static void CreateLike (amrex::Vector<std::unique_ptr<T> >& V, const T& Other)
25  {
26  // Emplace a new T in V with the same size as Other and get a reference
27  V.emplace_back(std::make_unique<T>(Other.Geom(0), Other.ParticleDistributionMap(0), Other.ParticleBoxArray(0)));
28  T& pc = *V[V.size()-1];
29 
30  // We want the particles to have all the same position, cpu, etc.
31  // as in Other, so do a copy from Other to our new particle container.
32  Copy(pc, Other);
33  }
34 
35  static void Copy (T& Y, const T& Other)
36  {
37  // Copy the contents of Other into Y
38  const bool local = true;
39  Y.copyParticles(Other, local);
40  }
41 
42  static void Saxpy (T& Y, const amrex::Real a, T& X)
43  {
44  // Calculate Y += a * X using a particle-level saxpy function supplied by the particle container T
47 
48  int lev = 0;
49  TParIter pty(Y, lev);
50  TParIter ptx(X, lev);
51 
52  auto checkValid = [&]() -> bool {
53  bool pty_v = pty.isValid();
54  bool ptx_v = ptx.isValid();
55  AMREX_ASSERT(pty_v == ptx_v);
56  return pty_v && ptx_v;
57  };
58 
59  auto ptIncrement = [&](){ ++pty; ++ptx; };
60 
61 #ifdef AMREX_USE_OMP
62 #pragma omp parallel if (Gpu::notInLaunchRegion())
63 #endif
64  for (; checkValid(); ptIncrement())
65  {
66  const int npy = pty.numParticles();
67  const int npx = ptx.numParticles();
68  AMREX_ALWAYS_ASSERT(npy == npx);
69 
70  ParticleType* psy = &(pty.GetArrayOfStructs()[0]);
71  ParticleType* psx = &(ptx.GetArrayOfStructs()[0]);
72 
73  auto particle_apply_rhs = T::particle_apply_rhs;
74 
75  amrex::ParallelFor ( npy, [=] AMREX_GPU_DEVICE (int i) {
76  ParticleType& py = psy[i];
77  const ParticleType& px = psx[i];
78  particle_apply_rhs(py, a, px);
79  });
80  }
81  }
82 
83 };
84 #endif
85 
86 template<class T>
87 struct IntegratorOps<T, std::enable_if_t<std::is_same_v<amrex::Vector<amrex::MultiFab>, T> > >
88 {
89 
90  static void CreateLike (amrex::Vector<std::unique_ptr<T> >& V, const T& Other, bool Grow = false)
91  {
92  // Emplace a new T in V with the same size as Other
93  V.emplace_back(std::make_unique<T>());
94  for (auto const& other_mf : Other) {
95  IntVect nGrow = Grow ? other_mf.nGrowVect() : IntVect(0);
96  V.back()->push_back(amrex::MultiFab(other_mf.boxArray(), other_mf.DistributionMap(), other_mf.nComp(), nGrow));
97  }
98  }
99 
100  static void Copy (T& Y, const T& Other, const Vector<int>& scomp={}, const Vector<int>& ncomp={}, bool Grow = true)
101  {
102  // Copy the contents of Other into Y
103  const int size = Y.size();
104  bool specify_components = !scomp.empty() && ncomp.size() == scomp.size();
105  for (int i = 0; i < size; ++i) {
106  IntVect nGrow = Grow ? Other[i].nGrowVect() : IntVect(0);
107  const int iscomp = specify_components ? scomp[i] : 0;
108  const int incomp = specify_components ? ncomp[i] : Other[i].nComp();
109  if (incomp > 0) {
110  amrex::MultiFab::Copy(Y[i], Other[i], iscomp, iscomp, incomp, nGrow);
111  }
112  }
113  }
114 
115  static void Saxpy (T& Y, const amrex::Real a, const T& X, const Vector<int>& scomp={}, const Vector<int>& ncomp={}, bool Grow = false)
116  {
117  // Calculate Y += a * X
118  const int size = Y.size();
119  bool specify_components = !scomp.empty() && ncomp.size() == scomp.size();
120  for (int i = 0; i < size; ++i) {
121  IntVect nGrow = Grow ? X[i].nGrowVect() : IntVect(0);
122  const int iscomp = specify_components ? scomp[i] : 0;
123  const int incomp = specify_components ? ncomp[i] : X[i].nComp();
124  if (incomp > 0) {
125  amrex::MultiFab::Saxpy(Y[i], a, X[i], iscomp, iscomp, incomp, nGrow);
126  }
127  }
128  }
129 
130 };
131 
132 template<class T>
133 struct IntegratorOps<T, std::enable_if_t<std::is_same_v<amrex::MultiFab, T> > >
134 {
135 
136  static void CreateLike (amrex::Vector<std::unique_ptr<T> >& V, const T& Other, bool Grow = false)
137  {
138  // Emplace a new T in V with the same size as Other
139  IntVect nGrow = Grow ? Other.nGrowVect() : IntVect(0);
140  V.emplace_back(std::make_unique<T>(Other.boxArray(), Other.DistributionMap(), Other.nComp(), nGrow));
141  }
142 
143  static void Copy (T& Y, const T& Other, const int scomp=0, const int ncomp=-1, bool Grow = true)
144  {
145  // Copy the contents of Other into Y
146  IntVect nGrow = Grow ? Other.nGrowVect() : IntVect(0);
147  const int mf_ncomp = ncomp > 0 ? ncomp : Other.nComp();
148  amrex::MultiFab::Copy(Y, Other, scomp, scomp, mf_ncomp, nGrow);
149  }
150 
151  static void Saxpy (T& Y, const amrex::Real a, const T& X, const int scomp=0, const int ncomp=-1, bool Grow = false)
152  {
153  // Calculate Y += a * X
154  IntVect nGrow = Grow ? X.nGrowVect() : IntVect(0);
155  const int mf_ncomp = ncomp > 0 ? ncomp : X.nComp();
156  amrex::MultiFab::Saxpy(Y, a, X, scomp, scomp, mf_ncomp, nGrow);
157  }
158 
159 };
160 
161 template<class T>
163 {
164 protected:
168  std::function<void(T& rhs, T& state, const amrex::Real time)> Rhs;
169 
174  std::function<void(T& rhs, T& state, const amrex::Real time)> RhsIm;
175 
180  std::function<void(T& rhs, T& state, const amrex::Real time)> RhsEx;
181 
186  std::function<void(T& rhs, T& state, const amrex::Real time)> RhsFast;
187 
192  std::function<void (T&, amrex::Real)> post_stage_action;
193 
198  std::function<void (T&, amrex::Real)> post_step_action;
199 
204  std::function<void (T&, amrex::Real)> post_fast_stage_action;
205 
210  std::function<void (T&, amrex::Real)> post_fast_step_action;
211 
217 
221  amrex::Real time_step;
222 
226  amrex::Real previous_time_step;
227 
233 
238  amrex::Real fast_time_step;
239 
243  amrex::Long num_steps = 0;
244 
248  int max_steps = 500;
249 
253  amrex::Real rel_tol = 1.0e-4;
254 
258  amrex::Real abs_tol = 1.0e-9;
259 
264  amrex::Real fast_rel_tol = 1.0e-4;
265 
270  amrex::Real fast_abs_tol = 1.0e-9;
271 
272 
273 public:
274  IntegratorBase () = default;
275 
276  IntegratorBase (const T& /* S_data */) {}
277 
278  virtual ~IntegratorBase () = default;
279 
280  void set_rhs (std::function<void(T&, T&, const amrex::Real)> F)
281  {
282  Rhs = F;
283  }
284 
285  void set_imex_rhs (std::function<void(T&, T&, const amrex::Real)> Fi,
286  std::function<void(T&, T&, const amrex::Real)> Fe)
287  {
288  RhsIm = Fi;
289  RhsEx = Fe;
290  }
291 
292  void set_fast_rhs (std::function<void(T&, T&, const amrex::Real)> F)
293  {
294  RhsFast = F;
295  }
296 
297  void set_post_stage_action (std::function<void (T&, amrex::Real)> A)
298  {
299  post_stage_action = A;
300  }
301 
302  void set_post_step_action (std::function<void (T&, amrex::Real)> A)
303  {
304  post_step_action = A;
305  }
306 
307  void set_post_fast_stage_action (std::function<void (T&, amrex::Real)> A)
308  {
310  }
311 
312  void set_post_fast_step_action (std::function<void (T&, amrex::Real)> A)
313  {
315  }
316 
317  amrex::Real get_time_step ()
318  {
319  return time_step;
320  }
321 
322  void set_time_step (amrex::Real dt)
323  {
324  time_step = dt;
325  use_adaptive_time_step = false;
326  }
327 
329  {
330  use_adaptive_time_step = true;
331  }
332 
333  void set_fast_time_step (amrex::Real dt)
334  {
335  fast_time_step = dt;
337  }
338 
340  {
342  }
343 
344  void set_max_steps (int steps)
345  {
346  max_steps = steps;
347  }
348 
349  void set_tolerances (amrex::Real rtol, amrex::Real atol)
350  {
351  rel_tol = rtol;
352  abs_tol = atol;
353  }
354 
355  void set_fast_tolerances (amrex::Real rtol, amrex::Real atol)
356  {
357  fast_rel_tol = rtol;
358  fast_abs_tol = atol;
359  }
360 
365  virtual amrex::Real advance (T& S_old, T& S_new, amrex::Real time,
366  amrex::Real dt) = 0;
367 
371  virtual void evolve (T& S_out, const amrex::Real time_out) = 0;
372 
373  virtual void time_interpolate (const T& S_new, const T& S_old,
374  amrex::Real timestep_fraction, T& data) = 0;
375 
376  virtual void map_data (std::function<void(T&)> Map) = 0;
377 };
378 
379 }
380 
381 #endif
#define AMREX_ASSERT(EX)
Definition: AMReX_BLassert.H:38
#define AMREX_ALWAYS_ASSERT(EX)
Definition: AMReX_BLassert.H:50
#define AMREX_GPU_DEVICE
Definition: AMReX_GpuQualifiers.H:18
Definition: AMReX_IntegratorBase.H:163
void set_post_step_action(std::function< void(T &, amrex::Real)> A)
Definition: AMReX_IntegratorBase.H:302
void set_rhs(std::function< void(T &, T &, const amrex::Real)> F)
Definition: AMReX_IntegratorBase.H:280
void set_adaptive_step()
Definition: AMReX_IntegratorBase.H:328
bool use_adaptive_fast_time_step
Flag to enable/disable adaptive time stepping at the fast time scale in multirate methods (bool)
Definition: AMReX_IntegratorBase.H:232
amrex::Real fast_rel_tol
Relative tolerance for adaptive time stepping at the fast time scale (Real)
Definition: AMReX_IntegratorBase.H:264
amrex::Real rel_tol
Relative tolerance for adaptive time stepping (Real)
Definition: AMReX_IntegratorBase.H:253
void set_fast_tolerances(amrex::Real rtol, amrex::Real atol)
Definition: AMReX_IntegratorBase.H:355
void set_post_fast_stage_action(std::function< void(T &, amrex::Real)> A)
Definition: AMReX_IntegratorBase.H:307
std::function< void(T &, amrex::Real)> post_fast_stage_action
The post_stage_action function is called by the integrator on the computed stage just after it is com...
Definition: AMReX_IntegratorBase.H:204
virtual void time_interpolate(const T &S_new, const T &S_old, amrex::Real timestep_fraction, T &data)=0
void set_adaptive_fast_step()
Definition: AMReX_IntegratorBase.H:339
virtual amrex::Real advance(T &S_old, T &S_new, amrex::Real time, amrex::Real dt)=0
Take a single time step from (time, S_old) to (time + dt, S_new) with the given step size.
amrex::Real get_time_step()
Definition: AMReX_IntegratorBase.H:317
amrex::Real fast_abs_tol
Absolute tolerance for adaptive time stepping at the fast time scale (Real)
Definition: AMReX_IntegratorBase.H:270
amrex::Real fast_time_step
Current integrator fast time scale time step size with multirate methods (Real)
Definition: AMReX_IntegratorBase.H:238
std::function< void(T &rhs, T &state, const amrex::Real time)> RhsEx
RhsEx is the explicit right-hand-side function an ImEx integrator will use.
Definition: AMReX_IntegratorBase.H:180
std::function< void(T &, amrex::Real)> post_step_action
The post_step_action function is called by the integrator on the computed state just after it is comp...
Definition: AMReX_IntegratorBase.H:198
void set_post_stage_action(std::function< void(T &, amrex::Real)> A)
Definition: AMReX_IntegratorBase.H:297
std::function< void(T &, amrex::Real)> post_fast_step_action
The post_step_action function is called by the integrator on the computed state just after it is comp...
Definition: AMReX_IntegratorBase.H:210
int max_steps
Max number of internal steps before an error is returned (Long)
Definition: AMReX_IntegratorBase.H:248
virtual void evolve(T &S_out, const amrex::Real time_out)=0
Evolve the current (internal) integrator state to time_out.
std::function< void(T &rhs, T &state, const amrex::Real time)> RhsIm
RhsIm is the implicit right-hand-side function an ImEx integrator will use.
Definition: AMReX_IntegratorBase.H:174
virtual void map_data(std::function< void(T &)> Map)=0
void set_post_fast_step_action(std::function< void(T &, amrex::Real)> A)
Definition: AMReX_IntegratorBase.H:312
void set_imex_rhs(std::function< void(T &, T &, const amrex::Real)> Fi, std::function< void(T &, T &, const amrex::Real)> Fe)
Definition: AMReX_IntegratorBase.H:285
void set_max_steps(int steps)
Definition: AMReX_IntegratorBase.H:344
void set_time_step(amrex::Real dt)
Definition: AMReX_IntegratorBase.H:322
void set_tolerances(amrex::Real rtol, amrex::Real atol)
Definition: AMReX_IntegratorBase.H:349
std::function< void(T &rhs, T &state, const amrex::Real time)> Rhs
Rhs is the right-hand-side function the integrator will use.
Definition: AMReX_IntegratorBase.H:168
virtual ~IntegratorBase()=default
amrex::Long num_steps
Number of integrator time steps (Long)
Definition: AMReX_IntegratorBase.H:243
bool use_adaptive_time_step
Flag to enable/disable adaptive time stepping in single rate methods or at the slow time scale in mul...
Definition: AMReX_IntegratorBase.H:216
void set_fast_rhs(std::function< void(T &, T &, const amrex::Real)> F)
Definition: AMReX_IntegratorBase.H:292
std::function< void(T &, amrex::Real)> post_stage_action
The post_stage_action function is called by the integrator on the computed stage just after it is com...
Definition: AMReX_IntegratorBase.H:192
void set_fast_time_step(amrex::Real dt)
Definition: AMReX_IntegratorBase.H:333
IntegratorBase(const T &)
Definition: AMReX_IntegratorBase.H:276
amrex::Real time_step
Current integrator time step size (Real)
Definition: AMReX_IntegratorBase.H:221
std::function< void(T &rhs, T &state, const amrex::Real time)> RhsFast
RhsFast is the fast timescale right-hand-side function a multirate integrator will use.
Definition: AMReX_IntegratorBase.H:186
amrex::Real previous_time_step
Step size of the last completed step (Real)
Definition: AMReX_IntegratorBase.H:226
amrex::Real abs_tol
Absolute tolerance for adaptive time stepping (Real)
Definition: AMReX_IntegratorBase.H:258
A collection (stored as an array) of FArrayBox objects.
Definition: AMReX_MultiFab.H:38
static void Saxpy(MultiFab &dst, Real a, const MultiFab &src, int srccomp, int dstcomp, int numcomp, int nghost)
dst += a*src
Definition: AMReX_MultiFab.cpp:345
static void Copy(MultiFab &dst, const MultiFab &src, int srccomp, int dstcomp, int numcomp, int nghost)
Copy from src to dst including nghost ghost cells. The two MultiFabs MUST have the same underlying Bo...
Definition: AMReX_MultiFab.cpp:193
Definition: AMReX_ParIter.H:113
This class is a thin wrapper around std::vector. Unlike vector, Vector::operator[] provides bound che...
Definition: AMReX_Vector.H:27
Long size() const noexcept
Definition: AMReX_Vector.H:50
AMREX_GPU_HOST_DEVICE Long size(T const &b) noexcept
integer version
Definition: AMReX_GpuRange.H:26
Definition: AMReX_Amr.cpp:49
std::enable_if_t< std::is_integral_v< T > > ParallelFor(TypeList< CTOs... > ctos, std::array< int, sizeof...(CTOs)> const &runtime_options, T N, F &&f)
Definition: AMReX_CTOParallelForImpl.H:200
void Saxpy(MF &dst, typename MF::value_type a, MF const &src, int scomp, int dcomp, int ncomp, IntVect const &nghost)
dst += a * src
Definition: AMReX_FabArrayUtility.H:1646
void Copy(FabArray< DFAB > &dst, FabArray< SFAB > const &src, int srccomp, int dstcomp, int numcomp, int nghost)
Definition: AMReX_FabArray.H:179
IntVectND< AMREX_SPACEDIM > IntVect
Definition: AMReX_BaseFwd.H:30
static void CreateLike(amrex::Vector< std::unique_ptr< T > > &V, const T &Other, bool Grow=false)
Definition: AMReX_IntegratorBase.H:136
static void Saxpy(T &Y, const amrex::Real a, const T &X, const int scomp=0, const int ncomp=-1, bool Grow=false)
Definition: AMReX_IntegratorBase.H:151
static void Copy(T &Y, const T &Other, const int scomp=0, const int ncomp=-1, bool Grow=true)
Definition: AMReX_IntegratorBase.H:143
static void Saxpy(T &Y, const amrex::Real a, const T &X, const Vector< int > &scomp={}, const Vector< int > &ncomp={}, bool Grow=false)
Definition: AMReX_IntegratorBase.H:115
static void CreateLike(amrex::Vector< std::unique_ptr< T > > &V, const T &Other, bool Grow=false)
Definition: AMReX_IntegratorBase.H:90
static void Copy(T &Y, const T &Other, const Vector< int > &scomp={}, const Vector< int > &ncomp={}, bool Grow=true)
Definition: AMReX_IntegratorBase.H:100
Definition: AMReX_IntegratorBase.H:17
The struct used to store particles.
Definition: AMReX_Particle.H:295