41#ifndef NEKTAR_LIB_UTILITIES_TIME_INTEGRATION_DIRK_TIME_INTEGRATION_SCHEME
42#define NEKTAR_LIB_UTILITIES_TIME_INTEGRATION_DIRK_TIME_INTEGRATION_SCHEME
44#define LUE LIB_UTILITIES_EXPORT
58 std::vector<NekDouble> freeParams)
63 "Diagonally Implicit Runge Kutta integration scheme bad order "
65 std::to_string(order));
71 ASSERTL1((variant ==
"" || variant ==
"ES5" || variant ==
"ES6"),
72 "DIRK Time integration scheme bad variant: " + variant +
74 "Must blank, 'ES5', or 'ES6'");
93 std::string variant,
size_t order, std::vector<NekDouble> freeParams)
97 variant, order, freeParams);
107 phase->m_order = order;
109 std::string(
"DIRKOrder" + std::to_string(phase->m_order));
111 phase->m_numsteps = 1;
112 phase->m_numstages = phase->m_order;
126 switch (phase->m_order)
133 phase->m_A[0][0][0] = 1.0;
135 phase->m_B[0][0][0] = 1.0;
147 phase->m_A[0][0][0] = lambda;
148 phase->m_A[0][1][0] = 1.0 - lambda;
149 phase->m_A[0][1][1] = lambda;
151 phase->m_B[0][0][0] = 1.0 - lambda;
152 phase->m_B[0][0][1] = lambda;
160 constexpr NekDouble lambda = 0.4358665215;
162 phase->m_A[0][0][0] = lambda;
163 phase->m_A[0][1][0] = 0.5 * (1.0 - lambda);
164 phase->m_A[0][2][0] =
165 0.25 * (-6.0 * lambda * lambda + 16.0 * lambda - 1.0);
167 phase->m_A[0][1][1] = lambda;
169 phase->m_A[0][2][1] =
170 0.25 * (6.0 * lambda * lambda - 20.0 * lambda + 5.0);
171 phase->m_A[0][2][2] = lambda;
173 phase->m_B[0][0][0] =
174 0.25 * (-6.0 * lambda * lambda + 16.0 * lambda - 1.0);
175 phase->m_B[0][0][1] =
176 0.25 * (6.0 * lambda * lambda - 20.0 * lambda + 5.0);
177 phase->m_B[0][0][2] = lambda;
182 phase->m_numMultiStepValues = 1;
183 phase->m_numMultiStepImplicitDerivs = 0;
184 phase->m_numMultiStepExplicitDerivs = 0;
186 phase->m_timeLevelOffset[0] = 0;
188 phase->CheckAndVerify();
193 size_t order, std::vector<NekDouble> freeParams)
196 phase->m_order = order;
198 std::string(
"DIRKOrder" + std::to_string(phase->m_order) + variant);
200 phase->m_numsteps = 1;
201 char const &stage = variant.back();
202 phase->m_numstages = std::atoi(&stage);
217 switch (phase->m_order)
225 std::string(
"DIRKOrder3_ES" +
226 std::to_string(phase->m_numstages) +
230 (freeParams.size()) ? freeParams[0] : 9.0 / 40.0;
232 phase->m_A[0][0][0] = 0.0;
233 phase->m_A[0][1][0] = lambda;
234 phase->m_A[0][2][0] = 9.0 * (1.0 + ConstSqrt2) / 80.0;
235 phase->m_A[0][3][0] =
236 (22.0 + 15.0 * ConstSqrt2) / (80.0 * (1.0 + ConstSqrt2));
237 phase->m_A[0][4][0] = (2398.0 + 1205.0 * ConstSqrt2) /
238 (2835.0 * (4.0 + 3.0 * ConstSqrt2));
240 phase->m_A[0][1][1] = phase->m_A[0][1][0];
241 phase->m_A[0][2][1] = phase->m_A[0][2][0];
242 phase->m_A[0][3][1] = phase->m_A[0][3][0];
243 phase->m_A[0][4][1] = phase->m_A[0][4][0];
245 phase->m_A[0][2][2] = lambda;
246 phase->m_A[0][3][2] = -7.0 / (40.0 * (1.0 + ConstSqrt2));
247 phase->m_A[0][4][2] = -2374.0 * (1.0 + 2.0 * ConstSqrt2) /
248 (2835.0 * (5.0 + 3.0 * ConstSqrt2));
250 phase->m_A[0][3][3] = lambda;
251 phase->m_A[0][4][3] = 5827.0 / 7560.0;
253 phase->m_A[0][4][4] = lambda;
255 phase->m_B[0][0][0] = phase->m_A[0][4][0];
256 phase->m_B[0][0][1] = phase->m_A[0][4][1];
257 phase->m_B[0][0][2] = phase->m_A[0][4][2];
258 phase->m_B[0][0][3] = phase->m_A[0][4][3];
259 phase->m_B[0][0][4] = phase->m_A[0][4][4];
269 std::string(
"DIRKOrder4_ES" +
270 std::to_string(phase->m_numstages) +
274 (freeParams.size()) ? freeParams[0] : 1.0 / 4.0;
276 phase->m_A[0][0][0] = 0.0;
277 phase->m_A[0][1][0] = lambda;
278 phase->m_A[0][2][0] = (1.0 - ConstSqrt2) / 8.0;
279 phase->m_A[0][3][0] = (5.0 - 7.0 * ConstSqrt2) / 64.0;
280 phase->m_A[0][4][0] =
281 (-13796.0 - 54539.0 * ConstSqrt2) / 125000.0;
282 phase->m_A[0][5][0] = (1181.0 - 987.0 * ConstSqrt2) / 13782.0;
284 phase->m_A[0][1][1] = phase->m_A[0][1][0];
285 phase->m_A[0][2][1] = phase->m_A[0][2][0];
286 phase->m_A[0][3][1] = phase->m_A[0][3][0];
287 phase->m_A[0][4][1] = phase->m_A[0][4][0];
288 phase->m_A[0][5][1] = phase->m_A[0][5][0];
290 phase->m_A[0][2][2] = lambda;
291 phase->m_A[0][3][2] = 7.0 * (1.0 + ConstSqrt2) / 32.0;
292 phase->m_A[0][4][2] =
293 (506605.0 + 132109.0 * ConstSqrt2) / 437500.0;
294 phase->m_A[0][5][2] =
295 47.0 * (-267.0 + 1783.0 * ConstSqrt2) / 273343.0;
297 phase->m_A[0][3][3] = lambda;
298 phase->m_A[0][4][3] =
299 166.0 * (-97.0 + 376.0 * ConstSqrt2) / 109375.0;
300 phase->m_A[0][5][3] =
301 -16.0 * (-22922.0 + 3525.0 * ConstSqrt2) / 571953.0;
303 phase->m_A[0][4][4] = lambda;
304 phase->m_A[0][5][4] =
305 -15625.0 * (97.0 + 376.0 * ConstSqrt2) / 90749876.0;
307 phase->m_A[0][5][5] = lambda;
309 phase->m_B[0][0][0] = phase->m_A[0][5][0];
310 phase->m_B[0][0][1] = phase->m_A[0][5][1];
311 phase->m_B[0][0][2] = phase->m_A[0][5][2];
312 phase->m_B[0][0][3] = phase->m_A[0][5][3];
313 phase->m_B[0][0][4] = phase->m_A[0][5][4];
314 phase->m_B[0][0][5] = phase->m_A[0][5][5];
319 ASSERTL0(
false, std::string(
"ESDIRK of order" +
320 std::to_string(phase->m_order) +
326 phase->m_numMultiStepValues = 1;
327 phase->m_numMultiStepImplicitDerivs = 0;
328 phase->m_numMultiStepExplicitDerivs = 0;
330 phase->m_timeLevelOffset[0] = 0;
332 phase->CheckAndVerify();
338 return std::string(
"DIRK");
354 std::vector<NekDouble> freeParams)
357 boost::ignore_unused(variant, order);
361 [[maybe_unused]] std::string variant, [[maybe_unused]]
size_t order,
362 std::vector<NekDouble> freeParams)
382 std::vector<NekDouble> freeParams)
385 boost::ignore_unused(variant, order);
389 [[maybe_unused]] std::string variant, [[maybe_unused]]
size_t order,
390 std::vector<NekDouble> freeParams)
410 std::vector<NekDouble> freeParams)
413 boost::ignore_unused(variant, order);
417 [[maybe_unused]] std::string variant, [[maybe_unused]]
size_t order,
418 std::vector<NekDouble> freeParams)
438 std::vector<NekDouble> freeParams)
441 boost::ignore_unused(variant, order);
445 [[maybe_unused]] std::string variant, [[maybe_unused]]
size_t order,
446 std::vector<NekDouble> freeParams)
450 "ES5", 3, freeParams);
466 std::vector<NekDouble> freeParams)
469 boost::ignore_unused(variant, order);
473 [[maybe_unused]] std::string variant, [[maybe_unused]]
size_t order,
474 std::vector<NekDouble> freeParams)
478 "ES6", 4, freeParams);
#define ASSERTL0(condition, msg)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
static std::string className
static TimeIntegrationSchemeSharedPtr create(std::string variant, size_t order, std::vector< NekDouble > freeParams)
DIRKOrder1TimeIntegrationScheme(std::string variant, size_t order, std::vector< NekDouble > freeParams)
static std::string TimeIntegrationMethodLookupId
static TimeIntegrationSchemeSharedPtr create(std::string variant, size_t order, std::vector< NekDouble > freeParams)
static std::string className
DIRKOrder2TimeIntegrationScheme(std::string variant, size_t order, std::vector< NekDouble > freeParams)
static std::string TimeIntegrationMethodLookupId
DIRKOrder3_ES5TimeIntegrationScheme(std::string variant, size_t order, std::vector< NekDouble > freeParams)
static TimeIntegrationSchemeSharedPtr create(std::string variant, size_t order, std::vector< NekDouble > freeParams)
static std::string TimeIntegrationMethodLookupId
static std::string className
static std::string TimeIntegrationMethodLookupId
DIRKOrder3TimeIntegrationScheme(std::string variant, size_t order, std::vector< NekDouble > freeParams)
static TimeIntegrationSchemeSharedPtr create(std::string variant, size_t order, std::vector< NekDouble > freeParams)
static std::string className
static TimeIntegrationSchemeSharedPtr create(std::string variant, size_t order, std::vector< NekDouble > freeParams)
static std::string className
static std::string TimeIntegrationMethodLookupId
DIRKOrder4_ES6TimeIntegrationScheme(std::string variant, size_t order, std::vector< NekDouble > freeParams)
static LUE void SetupSchemeDataESDIRK(TimeIntegrationAlgorithmGLMSharedPtr &phase, std::string variant, size_t order, std::vector< NekDouble > freeParams)
static LUE void SetupSchemeData(TimeIntegrationAlgorithmGLMSharedPtr &phase, size_t order)
LUE NekDouble v_GetTimeStability() const override
LUE std::string v_GetName() const override
DIRKTimeIntegrationScheme(std::string variant, size_t order, std::vector< NekDouble > freeParams)
~DIRKTimeIntegrationScheme() override
static TimeIntegrationSchemeSharedPtr create(std::string variant, size_t order, std::vector< NekDouble > freeParams)
static std::string className
Base class for GLM time integration schemes.
TimeIntegrationAlgorithmGLMVector m_integration_phases
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
std::shared_ptr< TimeIntegrationAlgorithmGLM > TimeIntegrationAlgorithmGLMSharedPtr
@ eDiagonallyImplicit
Diagonally implicit scheme (e.g. the DIRK schemes)
std::vector< TimeIntegrationAlgorithmGLMSharedPtr > TimeIntegrationAlgorithmGLMVector
std::shared_ptr< TimeIntegrationScheme > TimeIntegrationSchemeSharedPtr
scalarT< T > sqrt(scalarT< T > in)