Block-Structured AMR Software Framework
Loading...
Searching...
No Matches
AMReX_SundialsIntegrator.H
Go to the documentation of this file.
1#ifndef AMREX_SUNDIALS_INTEGRATOR_H
2#define AMREX_SUNDIALS_INTEGRATOR_H
3
4#include <functional>
5
6#include <AMReX_Config.H>
7#include <AMReX_REAL.H>
8#include <AMReX_Vector.H>
9#include <AMReX_ParmParse.H>
10#include <AMReX_IntegratorBase.H>
11#include <AMReX_NVector_MultiFab.H>
12#include <AMReX_Sundials.H>
13
14#include <nvector/nvector_manyvector.h>
15#include <sunnonlinsol/sunnonlinsol_fixedpoint.h>
16#include <sunlinsol/sunlinsol_spgmr.h>
17#include <arkode/arkode_arkstep.h>
18#include <arkode/arkode_mristep.h>
19
20namespace amrex {
21
22struct SundialsUserData {
23 // ERK or DIRK right-hand side function
24 // EX-MRI or IM-MRI slow right-hand side function
25 std::function<int(amrex::Real, N_Vector, N_Vector, void*)> f;
26
27 // ImEx-RK right-hand side functions
28 // ImEx-MRI slow right-hand side functions
29 std::function<int(amrex::Real, N_Vector, N_Vector, void*)> fi;
30 std::function<int(amrex::Real, N_Vector, N_Vector, void*)> fe;
31
32 // MRI fast time scale right-hand side function
33 std::function<int(amrex::Real, N_Vector, N_Vector, void*)> ff;
34
35 // Post stage and step actions
36 std::function<int(amrex::Real, N_Vector, void*)> post_stage;
37 std::function<int(amrex::Real, N_Vector, void*)> post_step;
38
39 // Post fast stage and step actions
40 std::function<int(amrex::Real, N_Vector, void*)> post_fast_stage;
41 std::function<int(amrex::Real, N_Vector, void*)> post_fast_step;
42};
43
44namespace SundialsUserFun {
45 static int f (amrex::Real t, N_Vector y_data, N_Vector y_rhs, void *user_data) {
46 SundialsUserData* udata = static_cast<SundialsUserData*>(user_data);
47 return udata->f(t, y_data, y_rhs, user_data);
48 }
49
50 static int fi (amrex::Real t, N_Vector y_data, N_Vector y_rhs, void *user_data) {
51 SundialsUserData* udata = static_cast<SundialsUserData*>(user_data);
52 return udata->fi(t, y_data, y_rhs, user_data);
53 }
54
55 static int fe (amrex::Real t, N_Vector y_data, N_Vector y_rhs, void *user_data) {
56 SundialsUserData* udata = static_cast<SundialsUserData*>(user_data);
57 return udata->fe(t, y_data, y_rhs, user_data);
58 }
59
60 static int ff (amrex::Real t, N_Vector y_data, N_Vector y_rhs, void *user_data) {
61 SundialsUserData* udata = static_cast<SundialsUserData*>(user_data);
62 return udata->ff(t, y_data, y_rhs, user_data);
63 }
64
65 static int post_stage (amrex::Real t, N_Vector y_data, void *user_data) {
66 SundialsUserData* udata = static_cast<SundialsUserData*>(user_data);
67 return udata->post_stage(t, y_data, user_data);
68 }
69
70 static int post_step (amrex::Real t, N_Vector y_data, void *user_data) {
71 SundialsUserData* udata = static_cast<SundialsUserData*>(user_data);
72 return udata->post_step(t, y_data, user_data);
73 }
74
75 static int post_fast_stage (amrex::Real t, N_Vector y_data, void *user_data) {
76 SundialsUserData* udata = static_cast<SundialsUserData*>(user_data);
77 return udata->post_fast_stage(t, y_data, user_data);
78 }
79
80 static int post_fast_step (amrex::Real t, N_Vector y_data, void *user_data) {
81 SundialsUserData* udata = static_cast<SundialsUserData*>(user_data);
82 return udata->post_fast_step(t, y_data, user_data);
83 }
84}
85
86template<class T>
87class SundialsIntegrator : public IntegratorBase<T>
88{
89private:
90 using BaseT = IntegratorBase<T>;
91
92 // Method type: ERK, DIRK, IMEX-RK, EX-MRI, IM-MRI, IMEX-MRI
93 std::string type = "ERK";
94
95 // Use SUNDIALS default methods
96 std::string method = "DEFAULT"; // ERK, DIRK, or slow method with MRI
97 std::string method_e = "DEFAULT"; // ERK in IMEX-RK
98 std::string method_i = "DEFAULT"; // DIRK in IMEX-RK
99
100 // Fast method type (ERK or DIRK) and method
101 std::string fast_type = "ERK";
102 std::string fast_method = "DEFAULT";
103
104 // Nonlinear solver
105 std::string nonlinear_solver = "Newton";
106 int max_nonlinear_iters = 0;
107
108 std::string fast_nonlinear_solver = "Newton";
109 int fast_max_nonlinear_iters = 0;
110
111 // Linear solver
112 std::string linear_solver = "GMRES";
113 int max_linear_iters = 0;
114
115 std::string fast_linear_solver = "GMRES";
116 int fast_max_linear_iters = 0;
117
118 // SUNDIALS package flags, set based on type
119 bool use_ark = false;
120 bool use_mri = false;
121
122 // structure for interfacing with user-supplied functions
123 SundialsUserData udata;
124
125 // SUNDIALS context
126 //
127 // We should probably use context created by amrex:sundials::Initialize but
128 // that context is not MPI-aware
129 ::sundials::Context sunctx;
130
131 // Single rate or slow time scale
132 void *arkode_mem = nullptr;
133 SUNLinearSolver LS = nullptr;
134 SUNNonlinearSolver NLS = nullptr;
135
136 // Fast time scale
137 void *arkode_fast_mem = nullptr;
138 MRIStepInnerStepper fast_stepper = nullptr;
139 SUNLinearSolver fast_LS = nullptr;
140 SUNNonlinearSolver fast_NLS = nullptr;
141
142 // Integrator stop time
143 bool set_stop_time = false;
144 amrex::Real stop_time = 0.0;
145
146 // Max steps between returns
147 amrex::Long max_num_steps = 0;
148
149 void initialize_parameters ()
150 {
151 amrex::ParmParse pp("integration.sundials");
152
153 pp.query("type", type);
154 pp.query("method", method);
155 pp.query("method_e", method_e);
156 pp.query("method_i", method_i);
157
158 pp.query("fast_type", fast_type);
159 pp.query("fast_method", fast_method);
160
161 if (type == "ERK" || type == "DIRK" || type == "IMEX-RK") {
162 use_ark = true;
163 }
164 else if (type == "EX-MRI" || type == "IM-MRI" || type == "IMEX-MRI") {
165 use_mri = true;
166 }
167 else {
168 std::string msg("Unknown method type: ");
169 msg += type;
170 amrex::Error(msg.c_str());
171 }
172
173 pp.query("nonlinear_solver", nonlinear_solver);
174 pp.query("max_nonlinear_iters", max_nonlinear_iters);
175
176 pp.query("fast_nonlinear_solver", fast_nonlinear_solver);
177 pp.query("fast_max_nonlinear_iters", fast_max_nonlinear_iters);
178
179 pp.query("linear_solver", linear_solver);
180 pp.query("max_linear_iters", max_linear_iters);
181
182 pp.query("fast_linear_solver", fast_linear_solver);
183 pp.query("fast_max_linear_iters", fast_max_linear_iters);
184
185 set_stop_time = pp.query("stop_time", stop_time);
186
187 pp.query("max_num_steps", max_num_steps);
188 }
189
190 void SetupRK (amrex::Real time, N_Vector y_data)
191 {
192 if (amrex::Verbose()) { amrex::Print() << "Using SUNDIALS time integrator\n"; }
193 int flag = 0;
194
195 // Create integrator and select method
196 if (type == "ERK") {
197 if (amrex::Verbose()) { amrex::Print() << "ERK method: " << method << "\n"; }
198 arkode_mem = ARKStepCreate(SundialsUserFun::f, nullptr, time, y_data, sunctx);
199 AMREX_ALWAYS_ASSERT(arkode_mem != nullptr);
200 if (method != "DEFAULT") {
201 flag = ARKStepSetTableName(arkode_mem, "ARKODE_DIRK_NONE", method.c_str());
202 AMREX_ALWAYS_ASSERT(flag == 0);
203 }
204 }
205 else if (type == "DIRK") {
206 if (amrex::Verbose()) { amrex::Print() << "DIRK method: " << method << "\n"; }
207 arkode_mem = ARKStepCreate(nullptr, SundialsUserFun::f, time, y_data, sunctx);
208 AMREX_ALWAYS_ASSERT(arkode_mem != nullptr);
209 if (method != "DEFAULT") {
210 flag = ARKStepSetTableName(arkode_mem, method.c_str(), "ARKODE_ERK_NONE");
211 AMREX_ALWAYS_ASSERT(flag == 0);
212 }
213 }
214 else if (type == "IMEX-RK") {
215 if (amrex::Verbose()) { amrex::Print() << "IMEX-RK method: " << method_i << " and "
216 << method_e << "\n"; }
217 arkode_mem = ARKStepCreate(SundialsUserFun::fe, SundialsUserFun::fi, time, y_data, sunctx);
218 AMREX_ALWAYS_ASSERT(arkode_mem != 0);
219 if (method_e != "DEFAULT" && method_i != "DEFAULT")
220 {
221 flag = ARKStepSetTableName(arkode_mem, method_i.c_str(), method_e.c_str());
222 AMREX_ALWAYS_ASSERT(flag == 0);
223 }
224 }
225
226 // Attach structure with user-supplied function wrappers
227 flag = ARKStepSetUserData(arkode_mem, &udata);
228 AMREX_ALWAYS_ASSERT(flag == 0);
229
230 // Set integrator tolerances
231 if (BaseT::use_adaptive_time_step || type == "DIRK" || type == "IMEX-RK") {
232 if (amrex::Verbose()) {
233 amrex::Print() << "Relative tolerance: " << BaseT::rel_tol << "\n";
234 amrex::Print() << "Absolute tolerance: " << BaseT::abs_tol << "\n";
235 }
236 flag = ARKStepSStolerances(arkode_mem, BaseT::rel_tol, BaseT::abs_tol);
237 AMREX_ALWAYS_ASSERT(flag == 0);
238 }
239
240 // Create and attach linear solver for implicit methods
241 if (type == "DIRK" || type == "IMEX-RK") {
242 if (amrex::Verbose()) {
243 amrex::Print() << "Nonlinear solver: " << nonlinear_solver << "\n";
244 amrex::Print() << "Max nonlinear iters: " << max_nonlinear_iters << "\n";
245 }
246 if (nonlinear_solver == "fixed-point") {
247 NLS = SUNNonlinSol_FixedPoint(y_data, 0, sunctx);
248 AMREX_ALWAYS_ASSERT(NLS != nullptr);
249 flag = ARKStepSetNonlinearSolver(arkode_mem, NLS);
250 AMREX_ALWAYS_ASSERT(flag == 0);
251 }
252 flag = ARKStepSetMaxNonlinIters(arkode_mem, max_nonlinear_iters);
253 AMREX_ALWAYS_ASSERT(flag == 0);
254
255 if (nonlinear_solver == "Newton") {
256 if (amrex::Verbose()) {
257 amrex::Print() << "Linear solver: " << linear_solver << "\n";
258 amrex::Print() << "Max linear iters: " << max_linear_iters << "\n";
259 }
260 LS = SUNLinSol_SPGMR(y_data, SUN_PREC_NONE, max_linear_iters, sunctx);
261 AMREX_ALWAYS_ASSERT(LS != nullptr);
262 flag = ARKStepSetLinearSolver(arkode_mem, LS, nullptr);
263 AMREX_ALWAYS_ASSERT(flag == 0);
264 }
265 }
266
267 // Set post stage and step function
268 flag = ARKStepSetPostprocessStageFn(arkode_mem, SundialsUserFun::post_stage);
269 AMREX_ALWAYS_ASSERT(flag == 0);
270 flag = ARKStepSetPostprocessStepFn(arkode_mem, SundialsUserFun::post_step);
271 AMREX_ALWAYS_ASSERT(flag == 0);
272
273 // Set a stop time
274 if (set_stop_time) {
275 if (amrex::Verbose()) { amrex::Print() << "Stop time: " << stop_time << "\n"; }
276 flag = ARKStepSetStopTime(arkode_mem, stop_time);
277 AMREX_ALWAYS_ASSERT(flag == 0);
278 }
279
280 // Set max number of steps between returns
281 ARKStepSetMaxNumSteps(arkode_mem, max_num_steps);
282 AMREX_ALWAYS_ASSERT(flag == 0);
283 }
284
285 void SetupMRI (amrex::Real time, N_Vector y_data)
286 {
287 if (amrex::Verbose()) { amrex::Print() << "Using SUNDIALS multirate time integrator\n"; }
288 int flag = 0;
289
290 // Create the fast integrator and select method
291 if (fast_type == "ERK") {
292 if (amrex::Verbose()) { amrex::Print() << "Fast ERK method: " << fast_method << "\n"; }
293 arkode_fast_mem = ARKStepCreate(SundialsUserFun::ff, nullptr, time, y_data, sunctx);
294 AMREX_ALWAYS_ASSERT(arkode_fast_mem != nullptr);
295 if (fast_method != "DEFAULT") {
296 flag = ARKStepSetTableName(arkode_fast_mem, "ARKODE_DIRK_NONE", fast_method.c_str());
297 AMREX_ALWAYS_ASSERT(flag == 0);
298 }
299 }
300 else if (fast_type == "DIRK") {
301 if (amrex::Verbose()) { amrex::Print() << "Fast DIRK method: " << fast_method << "\n"; }
302 arkode_fast_mem = ARKStepCreate(nullptr, SundialsUserFun::ff, time, y_data, sunctx);
303 AMREX_ALWAYS_ASSERT(arkode_fast_mem != nullptr);
304 if (fast_method != "DEFAULT") {
305 flag = ARKStepSetTableName(arkode_fast_mem, fast_method.c_str(), "ARKODE_ERK_NONE");
306 AMREX_ALWAYS_ASSERT(flag == 0);
307 }
308
309 if (amrex::Verbose()) {
310 amrex::Print() << "Fast nonlinear solver: " << fast_nonlinear_solver << "\n";
311 amrex::Print() << "Fast max nonlinear iters: " << fast_max_nonlinear_iters << "\n";
312 }
313 if (fast_nonlinear_solver == "fixed-point") {
314 fast_NLS = SUNNonlinSol_FixedPoint(y_data, 0, sunctx);
315 AMREX_ALWAYS_ASSERT(fast_NLS != nullptr);
316 flag = ARKStepSetNonlinearSolver(arkode_mem, fast_NLS);
317 AMREX_ALWAYS_ASSERT(flag == 0);
318 }
319 flag = ARKStepSetMaxNonlinIters(arkode_mem, fast_max_nonlinear_iters);
320 AMREX_ALWAYS_ASSERT(flag == 0);
321
322 if (fast_nonlinear_solver == "Newton") {
323 if (amrex::Verbose()) {
324 amrex::Print() << "Linear solver: " << fast_linear_solver << "\n";
325 amrex::Print() << "Max linear iters: " << fast_max_linear_iters << "\n";
326 }
327 fast_LS = SUNLinSol_SPGMR(y_data, SUN_PREC_NONE, fast_max_linear_iters, sunctx);
328 AMREX_ALWAYS_ASSERT(fast_LS != nullptr);
329 flag = ARKStepSetLinearSolver(arkode_mem, fast_LS, nullptr);
330 AMREX_ALWAYS_ASSERT(flag == 0);
331 }
332 }
333
334 // Attach structure with user-supplied function wrappers
335 flag = ARKStepSetUserData(arkode_fast_mem, &udata);
336 AMREX_ALWAYS_ASSERT(flag == 0);
337
338 // Set integrator tolerances
339 if (BaseT::use_adaptive_fast_time_step || fast_type == "DIRK" || fast_type == "IMEX-RK") {
340 if (amrex::Verbose()) {
341 amrex::Print() << "Fast relative tolerance: " << BaseT::fast_rel_tol << "\n";
342 amrex::Print() << "Fast absolute tolerance: " << BaseT::fast_abs_tol << "\n";
343 }
344 flag = ARKStepSStolerances(arkode_fast_mem, BaseT::fast_rel_tol, BaseT::fast_abs_tol);
345 AMREX_ALWAYS_ASSERT(flag == 0);
346 }
347
348 // Set post stage and step function
349 flag = ARKStepSetPostprocessStageFn(arkode_fast_mem, SundialsUserFun::post_fast_stage);
350 AMREX_ALWAYS_ASSERT(flag == 0);
351 flag = ARKStepSetPostprocessStepFn(arkode_fast_mem, SundialsUserFun::post_fast_step);
352 AMREX_ALWAYS_ASSERT(flag == 0);
353
354 // Set max number of steps between returns
355 ARKStepSetMaxNumSteps(arkode_fast_mem, max_num_steps);
356 AMREX_ALWAYS_ASSERT(flag == 0);
357
358 // Wrap fast integrator as an inner stepper
359 flag = ARKStepCreateMRIStepInnerStepper(arkode_fast_mem, &fast_stepper);
360 AMREX_ALWAYS_ASSERT(flag == 0);
361
362 // Create slow integrator
363 if (type == "EX-MRI") {
364 if (amrex::Verbose()) { amrex::Print() << "EX-MRI method: " << method << "\n"; }
365 arkode_mem = MRIStepCreate(SundialsUserFun::f, nullptr, time, y_data,
366 fast_stepper, sunctx);
367 AMREX_ALWAYS_ASSERT(arkode_mem != nullptr);
368 }
369 else if (type == "IM-MRI") {
370 if (amrex::Verbose()) { amrex::Print() << "IM-MRI method: " << method << "\n"; }
371 arkode_mem = MRIStepCreate(nullptr, SundialsUserFun::f, time, y_data,
372 fast_stepper, sunctx);
373 AMREX_ALWAYS_ASSERT(arkode_mem != nullptr);
374 }
375 else if (type == "IMEX-MRI") {
376 if (amrex::Verbose()) { amrex::Print() << "IMEX-MRI method: " << method << "\n"; }
377 arkode_mem = MRIStepCreate(SundialsUserFun::fe, SundialsUserFun::fi,
378 time, y_data, fast_stepper, sunctx);
379 AMREX_ALWAYS_ASSERT(arkode_mem != nullptr);
380 }
381
382 // Set method
383 if (method != "DEFAULT") {
384 MRIStepCoupling MRIC = MRIStepCoupling_LoadTableByName(method.c_str());
385 AMREX_ALWAYS_ASSERT(MRIC != nullptr);
386 flag = MRIStepSetCoupling(arkode_mem, MRIC);
387 AMREX_ALWAYS_ASSERT(flag == 0);
388 MRIStepCoupling_Free(MRIC);
389 }
390
391 // Attach structure with user-supplied function wrappers
392 flag = MRIStepSetUserData(arkode_mem, &udata);
393 AMREX_ALWAYS_ASSERT(flag == 0);
394
395 // Set integrator tolerances
396 if (BaseT::use_adaptive_time_step || fast_type == "IM-MRI" || fast_type == "IMEX-MRI") {
397 if (amrex::Verbose()) {
398 amrex::Print() << "Relative tolerance: " << BaseT::rel_tol << "\n";
399 amrex::Print() << "Absolute tolerance: " << BaseT::abs_tol << "\n";
400 }
401 flag = MRIStepSStolerances(arkode_mem, BaseT::rel_tol, BaseT::abs_tol);
402 AMREX_ALWAYS_ASSERT(flag == 0);
403 }
404
405 // Create and attach linear solver
406 if (type == "IM-MRI" || type == "IMEX-MRI") {
407 if (amrex::Verbose()) {
408 amrex::Print() << "Nonlinear solver: " << nonlinear_solver << "\n";
409 amrex::Print() << "Max nonlinear iters: " << max_nonlinear_iters << "\n";
410 }
411 if (nonlinear_solver == "fixed-point") {
412 NLS = SUNNonlinSol_FixedPoint(y_data, 0, sunctx);
413 AMREX_ALWAYS_ASSERT(NLS != nullptr);
414 flag = ARKStepSetNonlinearSolver(arkode_mem, NLS);
415 AMREX_ALWAYS_ASSERT(flag == 0);
416 }
417 flag = ARKStepSetMaxNonlinIters(arkode_mem, max_nonlinear_iters);
418 AMREX_ALWAYS_ASSERT(flag == 0);
419
420 if (nonlinear_solver == "Newton") {
421 if (amrex::Verbose()) {
422 amrex::Print() << "Linear solver: " << linear_solver << "\n";
423 amrex::Print() << "Max linear iters: " << max_linear_iters << "\n";
424 }
425 LS = SUNLinSol_SPGMR(y_data, SUN_PREC_NONE, max_linear_iters, sunctx);
426 AMREX_ALWAYS_ASSERT(LS != nullptr);
427 flag = ARKStepSetLinearSolver(arkode_mem, LS, nullptr);
428 AMREX_ALWAYS_ASSERT(flag == 0);
429 }
430 }
431
432 // Set post stage and step function
433 flag = MRIStepSetPostprocessStageFn(arkode_mem, SundialsUserFun::post_stage);
434 AMREX_ALWAYS_ASSERT(flag == 0);
435 flag = MRIStepSetPostprocessStepFn(arkode_mem, SundialsUserFun::post_step);
436 AMREX_ALWAYS_ASSERT(flag == 0);
437
438 // Set a stop time
439 if (set_stop_time) {
440 if (amrex::Verbose()) { amrex::Print() << "Stop time: " << stop_time << "\n"; }
441 flag = ARKStepSetStopTime(arkode_mem, stop_time);
442 AMREX_ALWAYS_ASSERT(flag == 0);
443 }
444
445 // Set max number of steps between returns
446 ARKStepSetMaxNumSteps(arkode_mem, max_num_steps);
447 AMREX_ALWAYS_ASSERT(flag == 0);
448 }
449
450 // -------------------------------------
451 // Vector<MultiFab> / N_Vector Utilities
452 // -------------------------------------
453
454 // Utility to unpack a SUNDIALS ManyVector into a vector of MultiFabs
455 void unpack_vector (N_Vector y_data, amrex::Vector<amrex::MultiFab>& S_data)
456 {
457 const int num_vecs = N_VGetNumSubvectors_ManyVector(y_data);
458 S_data.resize(num_vecs);
459
460 for(int i = 0; i < num_vecs; i++)
461 {
462 S_data.at(i) = amrex::MultiFab(*amrex::sundials::getMFptr(N_VGetSubvector_ManyVector(y_data, i)),
463 amrex::make_alias,
464 0,
465 amrex::sundials::getMFptr(N_VGetSubvector_ManyVector(y_data, i))->nComp());
466 }
467 };
468
469 // Utility to wrap vector of MultiFabs as a SUNDIALS ManyVector
470 N_Vector wrap_data (amrex::Vector<amrex::MultiFab>& S_data)
471 {
472 auto get_length = [&](int index) -> sunindextype {
473 auto* p_mf = &S_data[index];
474 return p_mf->nComp() * (p_mf->boxArray()).numPts();
475 };
476
477 sunindextype NV_len = S_data.size();
478 N_Vector* NV_array = new N_Vector[NV_len];
479
480 for (int i = 0; i < NV_len; ++i) {
481 NV_array[i] = amrex::sundials::N_VMake_MultiFab(get_length(i),
482 &S_data[i], &sunctx);
483 }
484
485 N_Vector y_data = N_VNew_ManyVector(NV_len, NV_array, sunctx);
486
487 delete[] NV_array;
488
489 return y_data;
490 };
491
492 // Utility to wrap vector of MultiFabs as a SUNDIALS ManyVector
493 N_Vector copy_data (const amrex::Vector<amrex::MultiFab>& S_data)
494 {
495 auto get_length = [&](int index) -> sunindextype {
496 auto* p_mf = &S_data[index];
497 return p_mf->nComp() * (p_mf->boxArray()).numPts();
498 };
499
500 sunindextype NV_len = S_data.size();
501 N_Vector* NV_array = new N_Vector[NV_len];
502
503 for (int i = 0; i < NV_len; ++i) {
504 NV_array[i] = amrex::sundials::N_VNew_MultiFab(get_length(i),
505 S_data[i].boxArray(),
506 S_data[i].DistributionMap(),
507 S_data[i].nComp(),
508 S_data[i].nGrow(),
509 &sunctx);
510
511 MultiFab::Copy(*amrex::sundials::getMFptr(NV_array[i]),
512 S_data[i],
513 0,
514 0,
515 S_data[i].nComp(),
516 S_data[i].nGrow());
517 }
518
519 N_Vector y_data = N_VNew_ManyVector(NV_len, NV_array, sunctx);
520
521 delete[] NV_array;
522
523 return y_data;
524 };
525
526 // -----------------------------
527 // MultiFab / N_Vector Utilities
528 // -----------------------------
529
530 // Utility to unpack a SUNDIALS Vector into a MultiFab
531 void unpack_vector (N_Vector y_data, amrex::MultiFab& S_data)
532 {
533 S_data = amrex::MultiFab(*amrex::sundials::getMFptr(y_data),
534 amrex::make_alias,
535 0,
536 amrex::sundials::getMFptr(y_data)->nComp());
537 };
538
539 // Utility to wrap a MultiFab as a SUNDIALS Vector
540 N_Vector wrap_data (amrex::MultiFab& S_data)
541 {
542 return amrex::sundials::N_VMake_MultiFab(S_data.nComp() * S_data.boxArray().numPts(),
543 &S_data, &sunctx);
544 };
545
546 // Utility to wrap a MultiFab as a SUNDIALS Vector
547 N_Vector copy_data (const amrex::MultiFab& S_data)
548 {
549 N_Vector y_data = amrex::sundials::N_VNew_MultiFab(S_data.nComp() * S_data.boxArray().numPts(),
550 S_data.boxArray(),
551 S_data.DistributionMap(),
552 S_data.nComp(),
553 S_data.nGrow(),
554 &sunctx);
555
556 MultiFab::Copy(*amrex::sundials::getMFptr(y_data),
557 S_data,
558 0,
559 0,
560 S_data.nComp(),
561 S_data.nGrow());
562
563 return y_data;
564 };
565
566public:
567 SundialsIntegrator () {}
568
569 SundialsIntegrator (const T& S_data, const amrex::Real time = 0.0)
570 {
571 initialize(S_data, time);
572 }
573
574 void initialize (const T& S_data, const amrex::Real time = 0.0)
575 {
576 initialize_parameters();
577 MPI_Comm mpi_comm = ParallelContext::CommunicatorSub();
578#if defined(SUNDIALS_VERSION_MAJOR) && (SUNDIALS_VERSION_MAJOR < 7)
579# ifdef AMREX_USE_MPI
580 sunctx = ::sundials::Context(&mpi_comm);
581# else
582 sunctx = ::sundials::Context(nullptr);
583# endif
584#else
585# ifdef AMREX_USE_MPI
586 sunctx = ::sundials::Context(mpi_comm);
587# else
588 sunctx = ::sundials::Context(SUN_COMM_NULL);
589# endif
590#endif
591
592 // Right-hand side function wrappers
593 udata.f = [&](amrex::Real rhs_time, N_Vector y_data, N_Vector y_rhs,
594 void * /* user_data */) -> int {
595
596 T S_data;
597 unpack_vector(y_data, S_data);
598
599 T S_rhs;
600 unpack_vector(y_rhs, S_rhs);
601
602 BaseT::Rhs(S_rhs, S_data, rhs_time);
603
604 return 0;
605 };
606
607 udata.fi = [&](amrex::Real rhs_time, N_Vector y_data, N_Vector y_rhs,
608 void * /* user_data */) -> int {
609
610 T S_data;
611 unpack_vector(y_data, S_data);
612
613 T S_rhs;
614 unpack_vector(y_rhs, S_rhs);
615
616 BaseT::RhsIm(S_rhs, S_data, rhs_time);
617
618 return 0;
619 };
620
621 udata.fe = [&](amrex::Real rhs_time, N_Vector y_data, N_Vector y_rhs,
622 void * /* user_data */) -> int {
623
624 T S_data;
625 unpack_vector(y_data, S_data);
626
627 T S_rhs;
628 unpack_vector(y_rhs, S_rhs);
629
630 BaseT::RhsEx(S_rhs, S_data, rhs_time);
631
632 return 0;
633 };
634
635 udata.ff = [&](amrex::Real rhs_time, N_Vector y_data, N_Vector y_rhs,
636 void * /* user_data */) -> int {
637
638 T S_data;
639 unpack_vector(y_data, S_data);
640
641 T S_rhs;
642 unpack_vector(y_rhs, S_rhs);
643
644 BaseT::RhsFast(S_rhs, S_data, rhs_time);
645
646 return 0;
647 };
648
649 udata.post_stage = [&](amrex::Real time, N_Vector y_data,
650 void * /* user_data */) -> int {
651
652 T S_data;
653 unpack_vector(y_data, S_data);
654
655 BaseT::post_stage_action(S_data, time);
656
657 return 0;
658 };
659
660 udata.post_step = [&](amrex::Real time, N_Vector y_data,
661 void * /* user_data */) -> int {
662
663 T S_data;
664 unpack_vector(y_data, S_data);
665
666 BaseT::post_step_action(S_data, time);
667
668 return 0;
669 };
670
671 udata.post_fast_stage = [&](amrex::Real time, N_Vector y_data,
672 void * /* user_data */) -> int {
673
674 T S_data;
675 unpack_vector(y_data, S_data);
676
677 BaseT::post_fast_stage_action(S_data, time);
678
679 return 0;
680 };
681
682 udata.post_fast_step = [&](amrex::Real time, N_Vector y_data,
683 void * /* user_data */) -> int {
684
685 T S_data;
686 unpack_vector(y_data, S_data);
687
688 BaseT::post_fast_step_action(S_data, time);
689
690 return 0;
691 };
692
693 N_Vector y_data = copy_data(S_data); // ideally just wrap and ignore const
694
695 if (use_ark) {
696 SetupRK(time, y_data);
697 }
698 else if (use_mri)
699 {
700 SetupMRI(time, y_data);
701 }
702
703 N_VDestroy(y_data);
704 }
705
706 virtual ~SundialsIntegrator () {
707 // Print integrator statistics
708 if (amrex::Verbose()) {
709 if (type == "EX-MRI" || type == "IM-MRI" || type == "IMEX-MRI") {
710 amrex::Print() << "Slow Time Integrator Stats\n";
711 if (ParallelDescriptor::IOProcessor()) {
712 MRIStepPrintAllStats(arkode_mem, stdout, SUN_OUTPUTFORMAT_TABLE);
713 }
714 amrex::Print() << "Fast Time Integrator Stats\n";
715 if (ParallelDescriptor::IOProcessor()) {
716 ARKStepPrintAllStats(arkode_fast_mem, stdout, SUN_OUTPUTFORMAT_TABLE);
717 }
718 } else {
719 amrex::Print() << "Time Integrator Stats\n";
720 if (ParallelDescriptor::IOProcessor()) {
721 ARKStepPrintAllStats(arkode_mem, stdout, SUN_OUTPUTFORMAT_TABLE);
722 }
723 }
724 }
725
726 // Clean up allocated memory
727 SUNLinSolFree(LS);
728 SUNLinSolFree(fast_LS);
729 SUNNonlinSolFree(NLS);
730 SUNNonlinSolFree(fast_NLS);
731 MRIStepInnerStepper_Free(&fast_stepper);
732 MRIStepFree(&arkode_fast_mem);
733 ARKStepFree(&arkode_mem);
734 }
735
736 amrex::Real advance (T& S_old, T& S_new, amrex::Real time, const amrex::Real dt) override
737 {
738 amrex::Real tout = time + dt;
739 amrex::Real tret;
740
741 N_Vector y_old = wrap_data(S_old);
742 N_Vector y_new = wrap_data(S_new);
743
744 if (use_ark) {
745 ARKStepReset(arkode_mem, time, y_old); // should probably resize
746 ARKStepSetFixedStep(arkode_mem, dt);
747 int flag = ARKStepEvolve(arkode_mem, tout, y_new, &tret, ARK_ONE_STEP);
748 AMREX_ALWAYS_ASSERT(flag >= 0);
749 }
750 else if (use_mri) {
751 MRIStepReset(arkode_mem, time, y_old); // should probably resize -- need to resize inner stepper
752 MRIStepSetFixedStep(arkode_mem, dt);
753 int flag = MRIStepEvolve(arkode_mem, tout, y_new, &tret, ARK_ONE_STEP);
754 AMREX_ALWAYS_ASSERT(flag >= 0);
755 } else {
756 Error("SUNDIALS integrator type not specified.");
757 }
758
759 N_VDestroy(y_old);
760 N_VDestroy(y_new);
761
762 return dt;
763 }
764
765 void evolve (T& S_out, const amrex::Real time_out) override
766 {
767 int flag = 0; // SUNDIALS return status
768 amrex::Real time_ret; // SUNDIALS return time
769
770 N_Vector y_out = wrap_data(S_out);
771
772 if (use_ark) {
773 if (!BaseT::use_adaptive_time_step) {
774 ARKStepSetFixedStep(arkode_mem, BaseT::time_step);
775 }
776 flag = ARKStepEvolve(arkode_mem, time_out, y_out, &time_ret, ARK_NORMAL);
777 AMREX_ALWAYS_ASSERT(flag >= 0);
778 }
779 else if (use_mri) {
780 if (!BaseT::use_adaptive_time_step) {
781 MRIStepSetFixedStep(arkode_mem, BaseT::time_step);
782 }
783 if (!BaseT::use_adaptive_fast_time_step) {
784 ARKStepSetFixedStep(arkode_fast_mem, BaseT::fast_time_step);
785 }
786 flag = MRIStepEvolve(arkode_mem, time_out, y_out, &time_ret, ARK_NORMAL);
787 AMREX_ALWAYS_ASSERT(flag >= 0);
788 } else {
789 Error("SUNDIALS integrator type not specified.");
790 }
791
792 N_VDestroy(y_out);
793 }
794
795 void time_interpolate (const T& /* S_new */, const T& /* S_old */, amrex::Real /* timestep_fraction */, T& /* data */) override {}
796
797 void map_data (std::function<void(T&)> /* Map */) override {}
798};
799
800}
801
802#endif
AMREX_ALWAYS_ASSERT
#define AMREX_ALWAYS_ASSERT(EX)
Definition AMReX_BLassert.H:50
pp
amrex::ParmParse pp
Input file parser instance for the given namespace.
Definition AMReX_HypreIJIface.cpp:15
AMReX_IntegratorBase.H
AMReX_NVector_MultiFab.H
AMReX_ParmParse.H
AMReX_REAL.H
AMReX_Sundials.H
AMReX_Vector.H
amrex::BoxArray::numPts
Long numPts() const noexcept
Returns the total number of cells contained in all boxes in the BoxArray.
Definition AMReX_BoxArray.cpp:394
amrex::FabArrayBase::nGrow
int nGrow(int direction=0) const noexcept
Return the grow factor that defines the region of definition.
Definition AMReX_FabArrayBase.H:78
amrex::FabArrayBase::DistributionMap
const DistributionMapping & DistributionMap() const noexcept
Return constant reference to associated DistributionMapping.
Definition AMReX_FabArrayBase.H:131
amrex::FabArrayBase::nComp
int nComp() const noexcept
Return number of variables (aka components) associated with each point.
Definition AMReX_FabArrayBase.H:83
amrex::FabArrayBase::boxArray
const BoxArray & boxArray() const noexcept
Return a constant reference to the BoxArray that defines the valid region associated with this FabArr...
Definition AMReX_FabArrayBase.H:95
amrex::IntegratorBase
Definition AMReX_IntegratorBase.H:164
amrex::IntegratorBase::use_adaptive_fast_time_step
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:246
amrex::IntegratorBase::fast_rel_tol
amrex::Real fast_rel_tol
Relative tolerance for adaptive time stepping at the fast time scale (Real)
Definition AMReX_IntegratorBase.H:278
amrex::IntegratorBase::rel_tol
amrex::Real rel_tol
Relative tolerance for adaptive time stepping (Real)
Definition AMReX_IntegratorBase.H:267
amrex::IntegratorBase::post_fast_stage_action
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:218
amrex::IntegratorBase::fast_abs_tol
amrex::Real fast_abs_tol
Absolute tolerance for adaptive time stepping at the fast time scale (Real)
Definition AMReX_IntegratorBase.H:284
amrex::IntegratorBase::fast_time_step
amrex::Real fast_time_step
Current integrator fast time scale time step size with multirate methods (Real)
Definition AMReX_IntegratorBase.H:252
amrex::IntegratorBase::RhsEx
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:194
amrex::IntegratorBase::post_step_action
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:212
amrex::IntegratorBase::post_fast_step_action
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:224
amrex::IntegratorBase::RhsIm
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:188
amrex::IntegratorBase::Rhs
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:182
amrex::IntegratorBase::use_adaptive_time_step
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:230
amrex::IntegratorBase::post_stage_action
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:206
amrex::IntegratorBase::time_step
amrex::Real time_step
Current integrator time step size (Real)
Definition AMReX_IntegratorBase.H:235
amrex::IntegratorBase::RhsFast
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:200
amrex::IntegratorBase::abs_tol
amrex::Real abs_tol
Absolute tolerance for adaptive time stepping (Real)
Definition AMReX_IntegratorBase.H:272
amrex::MultiFab
A collection (stored as an array) of FArrayBox objects.
Definition AMReX_MultiFab.H:40
amrex::MultiFab::Copy
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
amrex::ParmParse
Parse Parameters From Command Line and Input Files.
Definition AMReX_ParmParse.H:346
amrex::ParmParse::query
int query(std::string_view name, bool &ref, int ival=FIRST) const
Same as querykth() but searches for the last occurrence of name.
Definition AMReX_ParmParse.cpp:1450
amrex::Print
This class provides the user with a few print options.
Definition AMReX_Print.H:35
amrex::SundialsIntegrator
Definition AMReX_SundialsIntegrator.H:88
amrex::SundialsIntegrator::time_interpolate
void time_interpolate(const T &, const T &, amrex::Real, T &) override
Definition AMReX_SundialsIntegrator.H:795
amrex::SundialsIntegrator::SundialsIntegrator
SundialsIntegrator()
Definition AMReX_SundialsIntegrator.H:567
amrex::SundialsIntegrator::initialize
void initialize(const T &S_data, const amrex::Real time=0.0)
Definition AMReX_SundialsIntegrator.H:574
amrex::SundialsIntegrator::advance
amrex::Real advance(T &S_old, T &S_new, amrex::Real time, const amrex::Real dt) override
Take a single time step from (time, S_old) to (time + dt, S_new) with the given step size.
Definition AMReX_SundialsIntegrator.H:736
amrex::SundialsIntegrator::evolve
void evolve(T &S_out, const amrex::Real time_out) override
Evolve the current (internal) integrator state to time_out.
Definition AMReX_SundialsIntegrator.H:765
amrex::SundialsIntegrator::~SundialsIntegrator
virtual ~SundialsIntegrator()
Definition AMReX_SundialsIntegrator.H:706
amrex::SundialsIntegrator::SundialsIntegrator
SundialsIntegrator(const T &S_data, const amrex::Real time=0.0)
Definition AMReX_SundialsIntegrator.H:569
amrex::SundialsIntegrator::map_data
void map_data(std::function< void(T &)>) override
Definition AMReX_SundialsIntegrator.H:797
amrex::Vector
This class is a thin wrapper around std::vector. Unlike vector, Vector::operator[] provides bound che...
Definition AMReX_Vector.H:28
amrex::Vector::size
Long size() const noexcept
Definition AMReX_Vector.H:53
amrex::Real
amrex_real Real
Floating Point Type for Fields.
Definition AMReX_REAL.H:79
amrex::Long
amrex_long Long
Definition AMReX_INT.H:30
amrex::ParallelDescriptor::IOProcessor
bool IOProcessor() noexcept
Is this CPU the I/O Processor? To get the rank number, call IOProcessorNumber()
Definition AMReX_ParallelDescriptor.H:289
amrex::ParallelContext::CommunicatorSub
MPI_Comm CommunicatorSub() noexcept
sub-communicator for current frame
Definition AMReX_ParallelContext.H:70
amrex::SundialsUserFun::fi
static int fi(amrex::Real t, N_Vector y_data, N_Vector y_rhs, void *user_data)
Definition AMReX_SundialsIntegrator.H:50
amrex::SundialsUserFun::fe
static int fe(amrex::Real t, N_Vector y_data, N_Vector y_rhs, void *user_data)
Definition AMReX_SundialsIntegrator.H:55
amrex::SundialsUserFun::post_fast_step
static int post_fast_step(amrex::Real t, N_Vector y_data, void *user_data)
Definition AMReX_SundialsIntegrator.H:80
amrex::SundialsUserFun::f
static int f(amrex::Real t, N_Vector y_data, N_Vector y_rhs, void *user_data)
Definition AMReX_SundialsIntegrator.H:45
amrex::SundialsUserFun::post_step
static int post_step(amrex::Real t, N_Vector y_data, void *user_data)
Definition AMReX_SundialsIntegrator.H:70
amrex::SundialsUserFun::ff
static int ff(amrex::Real t, N_Vector y_data, N_Vector y_rhs, void *user_data)
Definition AMReX_SundialsIntegrator.H:60
amrex::SundialsUserFun::post_stage
static int post_stage(amrex::Real t, N_Vector y_data, void *user_data)
Definition AMReX_SundialsIntegrator.H:65
amrex::SundialsUserFun::post_fast_stage
static int post_fast_stage(amrex::Real t, N_Vector y_data, void *user_data)
Definition AMReX_SundialsIntegrator.H:75
amrex::mpidatatypes::MPI_Comm
int MPI_Comm
Definition AMReX_ccse-mpi.H:51
amrex::sundials::N_VMake_MultiFab
N_Vector N_VMake_MultiFab(sunindextype length, amrex::MultiFab *v_mf, ::sundials::Context *sunctx)
Definition AMReX_NVector_MultiFab.cpp:103
amrex::sundials::getMFptr
amrex::MultiFab *& getMFptr(N_Vector v)
Definition AMReX_NVector_MultiFab.cpp:228
amrex::sundials::N_VNew_MultiFab
N_Vector N_VNew_MultiFab(sunindextype length, const amrex::BoxArray &ba, const amrex::DistributionMapping &dm, sunindextype nComp, sunindextype nGhost, ::sundials::Context *sunctx)
Definition AMReX_NVector_MultiFab.cpp:78
amrex
Definition AMReX_Amr.cpp:49
amrex::make_alias
@ make_alias
Definition AMReX_MakeType.H:7
amrex::nComp
int nComp(FabArrayBase const &fa)
Definition AMReX_FabArrayBase.cpp:2854
amrex::DistributionMap
DistributionMapping const & DistributionMap(FabArrayBase const &fa)
Definition AMReX_FabArrayBase.cpp:2869
amrex::Error
void Error(const std::string &msg)
Print out message to cerr and exit via amrex::Abort().
Definition AMReX.cpp:224
amrex::Verbose
int Verbose() noexcept
Definition AMReX.cpp:169
amrex::int
const int[]
Definition AMReX_BLProfiler.cpp:1664
amrex::boxArray
BoxArray const & boxArray(FabArrayBase const &fa)
Definition AMReX_FabArrayBase.cpp:2864
amrex::SundialsUserData
Definition AMReX_SundialsIntegrator.H:22
amrex::SundialsUserData::fi
std::function< int(amrex::Real, N_Vector, N_Vector, void *)> fi
Definition AMReX_SundialsIntegrator.H:29
amrex::SundialsUserData::post_fast_stage
std::function< int(amrex::Real, N_Vector, void *)> post_fast_stage
Definition AMReX_SundialsIntegrator.H:40
amrex::SundialsUserData::post_step
std::function< int(amrex::Real, N_Vector, void *)> post_step
Definition AMReX_SundialsIntegrator.H:37
amrex::SundialsUserData::f
std::function< int(amrex::Real, N_Vector, N_Vector, void *)> f
Definition AMReX_SundialsIntegrator.H:25
amrex::SundialsUserData::post_stage
std::function< int(amrex::Real, N_Vector, void *)> post_stage
Definition AMReX_SundialsIntegrator.H:36
amrex::SundialsUserData::ff
std::function< int(amrex::Real, N_Vector, N_Vector, void *)> ff
Definition AMReX_SundialsIntegrator.H:33
amrex::SundialsUserData::post_fast_step
std::function< int(amrex::Real, N_Vector, void *)> post_fast_step
Definition AMReX_SundialsIntegrator.H:41
amrex::SundialsUserData::fe
std::function< int(amrex::Real, N_Vector, N_Vector, void *)> fe
Definition AMReX_SundialsIntegrator.H:30