DNDSR 0.1.0.dev1+gcd065ad
Distributed Numeric Data Structure for CFV
Loading...
Searching...
No Matches
ODE.hpp
Go to the documentation of this file.
1#pragma once
2#include "DNDS/Defines.hpp"
3#include "DNDS/MPI.hpp"
4#include "Scalar.hpp"
5#include "DNDS/JsonUtil.hpp"
6
7namespace DNDS::ODE
8{
9 template <class TDATA, class TDTAU>
10 class ImplicitDualTimeStep
11 {
12 public:
13 using Frhs = std::function<void(TDATA &, TDATA &, TDTAU &, int, real, int)>;
14 using Fdt = std::function<void(TDATA &, TDTAU &, real, int)>;
15 // x, res, resOther(=res-alpha*rhs), dTau, dt, alpha, xinc, iter, ct, stage
16 using Fsolve = std::function<void(TDATA &, TDATA &, TDATA &, TDTAU &, real, real, TDATA &, int, real, int)>;
17 using Fstop = std::function<bool(int, TDATA &, int)>;
18 using Fincrement = std::function<void(TDATA &, TDATA &, real, int)>;
19
20 virtual void Step(TDATA &x, TDATA &xinc, const Frhs &frhs, const Fdt &fdt, const Fsolve &fsolve,
21 int maxIter, const Fstop &fstop, const Fincrement &fincrement, real dt) = 0;
22 virtual ~ImplicitDualTimeStep() = default;
23
24 virtual TDATA &getLatestRHS() = 0;
25
26 virtual TDATA &getRHS(int i) = 0;
27
28 virtual TDATA &getRES(int i) = 0;
29
30 virtual void SetExtraParams(const nlohmann::ordered_json &j)
31 {
32 for (auto &[k, v] : j.items())
33 DNDS_assert_info(false, "no extra params to handle! key is " + k);
34 };
35 };
36
37 template <class TDATA, class TDTAU>
38 class ImplicitEulerDualTimeStep : public ImplicitDualTimeStep<TDATA, TDTAU>
39 {
40 public:
41 using tBase = ImplicitDualTimeStep<TDATA, TDTAU>;
42 using Frhs = typename tBase::Frhs;
43 using Fdt = typename tBase::Fdt;
44 using Fsolve = typename tBase::Fsolve;
45 using Fstop = typename tBase::Fstop;
46 using Fincrement = typename tBase::Fincrement;
47
48 TDTAU dTau;
49 std::vector<TDATA> rhsbuf;
50 TDATA rhs;
51 TDATA resOther;
52 TDATA xLast;
53 TDATA xInc;
54 index DOF;
55
56 /**
57 * @brief mind that NDOF is the dof of dt
58 * finit(TDATA& data)
59 */
60 template <class Finit, class FinitDtau>
61 ImplicitEulerDualTimeStep(
62 index NDOF, Finit &&finit = [](TDATA &) {}, FinitDtau &&finitDtau = [](TDTAU &) {}) : DOF(NDOF)
63 {
64 rhsbuf.resize(1);
65 for (auto &i : rhsbuf)
66 finit(i);
67 finit(rhs);
68 finit(resOther);
69 finit(xLast);
70 finit(xInc);
71 finitDtau(dTau);
72 }
73
74 /**
75 * @brief
76 * frhs(TDATA &rhs, TDATA &x)
77 * fdt(TDTAU& dTau)
78 * fsolve(TDATA &x, TDATA &rhs, TDTAU& dTau, real dt, real alphaDiag, TDATA &xinc)
79 * bool fstop(int iter, TDATA &xinc, int iInternal)
80 */
81 virtual void Step(TDATA &x, TDATA &xinc, const Frhs &frhs, const Fdt &fdt, const Fsolve &fsolve,
82 int maxIter, const Fstop &fstop, const Fincrement &fincrement, real dt) override
83 {
84 xLast = x;
85 for (int iter = 1; iter <= maxIter; iter++)
86 {
87 fdt(x, dTau, 1, 0);
88
89 frhs(rhsbuf[0], x, dTau, iter, 1.0, 0);
90 rhs = xLast;
91 rhs *= 1.0 / dt;
92 resOther = rhs;
93 rhs.addTo(x, -1. / dt);
94 rhs += rhsbuf[0]; // crhs = rhs + (x_i - x_j) / dt
95
96 fsolve(x, rhs, resOther, dTau, dt, 1.0, xinc, iter, 1.0, 0);
97 fincrement(x, xinc, 1.0, 0);
98
99 if (fstop(iter, rhs, 1))
100 break;
101 }
102 }
103
104 virtual ~ImplicitEulerDualTimeStep() = default;
105
106 virtual TDATA &getLatestRHS() override
107 {
108 return rhsbuf[0];
109 }
110
111 TDATA &getRHS(int i) override
112 {
113 return rhsbuf[0];
114 }
115
116 TDATA &getRES(int i) override
117 {
118 return rhs;
119 }
120 };
121
122 template <class TDATA, class TDTAU>
123 class ImplicitSDIRK4DualTimeStep : public ImplicitDualTimeStep<TDATA, TDTAU>
124 {
125
126 Eigen::Matrix<real, -1, -1> butcherA;
127 Eigen::Vector<real, -1> butcherC;
128 Eigen::RowVector<real, -1> butcherB;
129 int nInnerStage = 3;
130 int schemeC = 0;
131 bool explicitFirst = false;
132 int hasLastEndPointR = 0;
133 int latestStage = 0;
134
135 public:
136 using tBase = ImplicitDualTimeStep<TDATA, TDTAU>;
137 using Frhs = typename tBase::Frhs;
138 using Fdt = typename tBase::Fdt;
139 using Fsolve = typename tBase::Fsolve;
140 using Fstop = typename tBase::Fstop;
141 using Fincrement = typename tBase::Fincrement;
142
143 TDTAU dTau;
144 std::vector<TDATA> rhsbuf;
145 TDATA rhs;
146 TDATA resOther;
147 TDATA xLast;
148 TDATA xIncPrev;
149 index DOF;
150
151 /**
152 * @brief mind that NDOF is the dof of dt
153 * finit(TDATA& data)
154 */
155 template <class Finit, class FinitDtau>
156 ImplicitSDIRK4DualTimeStep(
157 index NDOF, Finit &&finit = [](TDATA &) {}, FinitDtau &&finitDtau = [](TDTAU &) {}, int schemeCode = 0) : DOF(NDOF)
158 {
159
160 schemeC = schemeCode;
161
162 if (schemeCode == 0)
163 {
164#define _zeta 0.128886400515
165 nInnerStage = 3;
166 butcherA.resize(nInnerStage, nInnerStage);
167 butcherC.resize(nInnerStage);
168 butcherB.resize(nInnerStage);
169
170 butcherA << _zeta, 0, 0,
171 0.5 - _zeta, _zeta, 0,
172 2 * _zeta, 1 - 4 * _zeta, _zeta;
173 butcherC << _zeta, 0.5, 1 - _zeta;
174 butcherB << 1. / (6 * sqr(2 * _zeta - 1)),
175 (4 * sqr(_zeta) - 4 * _zeta + 2. / 3.) / sqr(2 * _zeta - 1),
176 1. / (6 * sqr(2 * _zeta - 1));
177#undef _zeta
178 }
179 else if (schemeCode == 1)
180 {
181 nInnerStage = 6;
182 explicitFirst = true;
183 butcherA.resize(nInnerStage, nInnerStage);
184 butcherC.resize(nInnerStage);
185 butcherB.resize(nInnerStage);
186
187 butcherA << verySmallReal, 0, 0, 0, 0, 0,
188 0.25, 0.25, 0, 0, 0, 0,
189 0.137776, -0.055776, 0.25, 0, 0, 0,
190 0.1446368660269822, -0.2239319076133447, 0.4492950415863626, 0.25, 0, 0,
191 0.09825878328356477, -0.5915442428196704, 0.8101210205756877, 0.283164405707806, 0.25, 0,
192 0.1579162951616714, 0, 0.1867589405240008, 0.6805652953093346, -0.2752405309950067, 0.25;
193 butcherB = butcherA(EigenLast, EigenAll);
194 butcherC << 0, 0.5, 0.332, 0.62, 0.849999966747388, 1;
195 }
196 else if (schemeCode == 2) // esdirk3
197 {
198 nInnerStage = 4;
199 explicitFirst = true;
200
201 butcherA.resize(nInnerStage, nInnerStage);
202 butcherC.resize(nInnerStage);
203 butcherB.resize(nInnerStage);
204 double alphaC = double(1767732205903) / double(4055673282236);
205
206 butcherA << verySmallReal, 0, 0, 0,
207 alphaC, alphaC, 0, 0,
208 real(2746238789719) / real(10658868560708), real(-640167445237) / real(6845629431997), alphaC, 0,
209 real(1471266399579) / real(7840856788654), real(-4482444167858) / real(7529755066697), real(11266239266428) / real(11593286722821), alphaC;
210
211 butcherB = butcherA(EigenLast, EigenAll);
212 butcherC = butcherA.rowwise().sum();
213 butcherC(0) = 0;
214 butcherC(3) = 1;
215 }
216 else if (schemeCode == 3) // trapezoid rule
217 {
218 nInnerStage = 2;
219 explicitFirst = true;
220 butcherA.resize(nInnerStage, nInnerStage);
221 butcherC.resize(nInnerStage);
222 butcherB.resize(nInnerStage);
223 butcherA << verySmallReal, 0,
224 0.5, 0.5;
225 butcherB = butcherA(EigenLast, EigenAll);
226 butcherC << 0, 1;
227 }
228 else if (schemeCode == 4) // esdirk2
229 {
230 nInnerStage = 3;
231 explicitFirst = true;
232 butcherA.resize(nInnerStage, nInnerStage);
233 butcherC.resize(nInnerStage);
234 butcherB.resize(nInnerStage);
235
236 real gamma = 1. - std::sqrt(2.0) * 0.5;
237 real b2 = (1. - 2 * gamma) / (4 * gamma);
238 double alphaC = double(1767732205903) / double(4055673282236);
239
240 butcherA << verySmallReal, 0, 0,
241 gamma, gamma, 0,
242 1 - b2 - gamma, b2, gamma;
243
244 butcherB = butcherA(EigenLast, EigenAll);
245 butcherC = butcherA.rowwise().sum();
246 butcherC(0) = 0;
247 butcherC(2) = 1;
248 }
249 else
250 {
251 DNDS_assert(false);
252 }
253
254 rhsbuf.resize(nInnerStage);
255 for (auto &i : rhsbuf)
256 finit(i);
257 finit(rhs);
258 finit(resOther);
259 finit(xLast);
260 finit(xIncPrev);
261 finitDtau(dTau);
262 }
263
264 /**
265 * @brief
266 * frhs(TDATA &rhs, TDATA &x)
267 * fdt(TDTAU& dTau)
268 * fsolve(TDATA &x, TDATA &rhs, TDTAU& dTau, real dt, real alphaDiag, TDATA &xinc)
269 * bool fstop(int iter, TDATA &xinc, int iInternal)
270 */
271 virtual void Step(TDATA &x, TDATA &xinc, const Frhs &frhs, const Fdt &fdt, const Fsolve &fsolve,
272 int maxIter, const Fstop &fstop, const Fincrement &fincrement, real dt) override
273 {
274 xLast = x;
275 for (int iB = 0; iB < nInnerStage; iB++)
276 {
277 x = xLast;
278 xIncPrev.setConstant(0.0);
279 int iter = 1;
280 for (; iter <= maxIter; iter++)
281 {
282 if (explicitFirst && iB == 0) // for esdirk first frhs evaluation
283 {
284 if (hasLastEndPointR)
285 rhsbuf[0] = rhsbuf[nInnerStage - 1];
286 else
287 frhs(rhsbuf[0], x, dTau, INT_MAX, butcherC(0), 0), latestStage = 0;
288 break;
289 }
290 fdt(x, dTau, butcherA(iB, iB), 0);
291
292 frhs(rhsbuf[iB], x, dTau, iter, butcherC(iB), 0), latestStage = iB;
293
294 // //!test explicit
295 // rhs = rhsbuf[iB];
296 // rhs *= dTau;
297 // xinc = rhs;
298 // x += xinc;
299 // if (fstop(iter, xinc, iB + 1))
300 // break;
301 // continue;
302 // //! test explicit
303
304 // rhsbuf[0] = rhs;
305 rhs = xLast;
306 rhs *= 1.0 / dt;
307 for (int jB = 0; jB < iB; jB++)
308 rhs.addTo(rhsbuf[jB], butcherA(iB, jB)); // crhs = rhs + (x_i - x_j) / dt
309 resOther = rhs;
310 rhs.addTo(x, -1. / dt);
311 rhs.addTo(rhsbuf[iB], butcherA(iB, iB));
312
313 fsolve(x, rhs, resOther, dTau, dt, butcherA(iB, iB), xinc, iter, butcherC(iB), 0);
314 // x += xinc;
315 fincrement(x, xinc, 1.0, 0);
316 // x.addTo(xIncPrev, -0.5);
317
318 xIncPrev = xinc;
319
320 if (fstop(iter, rhs, iB + 1))
321 break;
322 if (explicitFirst && iB == 0) // for esdirk
323 break;
324
325 // TODO: add time dependent rhs
326 }
327 if (iter > maxIter)
328 fstop(iter, rhs, iB + 1);
329 }
330 if (explicitFirst) // for esdirk
331 {
332 hasLastEndPointR = 1;
333 return;
334 }
335 x = xLast;
336 // for (int jB = 0; jB < nInnerStage; jB++)
337 // fincrement(x, rhsbuf[jB], butcherB(jB) * dt, 0); //!bad here
338 for (int jB = 0; jB < nInnerStage; jB++)
339 x.addTo(rhsbuf[jB], butcherB(jB) * dt);
340 }
341
342 virtual TDATA &getLatestRHS() override
343 {
344 return rhsbuf[latestStage];
345 }
346
347 TDATA &getRHS(int i) override
348 {
349 DNDS_assert(i >= 0 && i < rhsbuf.size());
350 return rhsbuf[i];
351 }
352
353 TDATA &getRES(int i) override
354 {
355 // todo: enable on-demanc calculation
356 return rhs;
357 }
358
359 virtual ~ImplicitSDIRK4DualTimeStep() = default;
360 };
361
362 template <class TDATA, class TDTAU>
363 class ImplicitBDFDualTimeStep : public ImplicitDualTimeStep<TDATA, TDTAU>
364 {
365
366 static const Eigen::Matrix<real, 4, 5> BDFCoefs;
367
368 public:
369 using tBase = ImplicitDualTimeStep<TDATA, TDTAU>;
370 using Frhs = typename tBase::Frhs;
371 using Fdt = typename tBase::Fdt;
372 using Fsolve = typename tBase::Fsolve;
373 using Fstop = typename tBase::Fstop;
374 using Fincrement = typename tBase::Fincrement;
375
376 TDTAU dTau;
377 std::vector<TDATA> xPrevs;
378 Eigen::VectorXd dtPrevs;
379 std::vector<TDATA> rhsbuf;
380 TDATA rhs;
381 TDATA resOther;
382 TDATA xLast;
383 TDATA xIncPrev;
384 TDATA resInc;
385
386 index DOF;
387 index cnPrev;
388 index prevStart;
389 index kBDF;
390
391 /**
392 * @brief mind that NDOF is the dof of dt
393 * finit(TDATA& data)
394 */
395 template <class Finit, class FinitDtau>
396 ImplicitBDFDualTimeStep(
397 index NDOF, Finit &&finit = [](TDATA &) {}, FinitDtau &&finitDtau = [](TDTAU &) {},
398 index k = 2) : DOF(NDOF), cnPrev(0), prevStart(k - 2), kBDF(k)
399 {
400 assert(k > 0 && k <= 4);
401
402 xPrevs.resize(k - 1);
403 dtPrevs.resize(k - 1);
404 for (auto &i : xPrevs)
405 finit(i);
406 rhsbuf.resize(1);
407 finit(rhsbuf[0]);
408 finit(rhs);
409 finit(resOther);
410 finit(resInc);
411 finit(xLast);
412 finit(xIncPrev);
413 finitDtau(dTau);
414 }
415
416 /**
417 * @brief
418 * frhs(TDATA &rhs, TDATA &x)
419 * fdt(TDTAU& dTau)
420 * fsolve(TDATA &x, TDATA &rhs, TDTAU& dTau, real dt, real alphaDiag, TDATA &xinc)
421 * bool fstop(int iter, TDATA &xinc, int iInternal)
422 */
423 virtual void Step(TDATA &x, TDATA &xinc, const Frhs &frhs, const Fdt &fdt, const Fsolve &fsolve,
424 int maxIter, const Fstop &fstop, const Fincrement &fincrement, real dt) override
425 {
426 index kCurrent = cnPrev + 1;
427 index prevSiz = kBDF - 1;
428 for (index iPrev = 0; iPrev < cnPrev; iPrev++)
429 assert(prevSiz && std::abs(dtPrevs[mod(iPrev + prevStart, prevSiz)] - dt) < dt * 1e-8);
430
431 xLast = x;
432 // x = xLast;
433 xIncPrev.setConstant(0.0);
434 int iter = 1;
435 for (; iter <= maxIter; iter++)
436 {
437 fdt(x, dTau, BDFCoefs(kCurrent - 1, 0), 0);
438 frhs(rhsbuf[0], x, dTau, iter, 1.0, 0);
439
440 rhs.setConstant(0.0);
441 rhs.addTo(xLast, BDFCoefs(kCurrent - 1, 1) / dt);
442 // std::cout << "add " << BDFCoefs(kCurrent - 1, 1) << " " << "last" << std::endl;
443 if (prevSiz)
444 for (index iPrev = 0; iPrev < cnPrev; iPrev++)
445 {
446 // std::cout << "add " << BDFCoefs(kCurrent - 1, 2 + iPrev) <<" " << mod(iPrev + prevStart, prevSiz) << std::endl;
447 rhs.addTo(xPrevs[mod(iPrev + prevStart, prevSiz)], BDFCoefs(kCurrent - 1, 2 + iPrev) / dt);
448 }
449 resOther = rhs;
450 rhs.addTo(x, -1. / dt);
451 rhs.addTo(rhsbuf[0], BDFCoefs(kCurrent - 1, 0));
452 fsolve(x, rhs, resOther, dTau, dt, BDFCoefs(kCurrent - 1, 0), xinc, iter, 1.0, 0);
453 //* xinc = (I/dtau-A*alphaDiag)\rhs
454
455 // std::cout << "BDF::\n";
456 // std::cout << kCurrent << " " << cnPrev<<" " << BDFCoefs(kCurrent - 1, 0) << std::endl;
457
458 // x += xinc;
459 fincrement(x, xinc, 1.0, 0);
460 // x.addTo(xIncPrev, -0.5);
461
462 xIncPrev = xinc;
463
464 if (fstop(iter, rhs, 1))
465 break;
466 }
467 if (iter > maxIter)
468 fstop(iter, rhs, 1);
469 if (prevSiz)
470 {
471 prevStart = mod(prevStart - 1, prevSiz);
472 // std::cout << dtPrevs.size() << " " << prevStart << std::endl;
473 xPrevs[prevStart] = xLast;
474 dtPrevs[prevStart] = dt;
475 cnPrev = std::min(cnPrev + 1, prevSiz);
476 }
477 }
478
479 using FresidualIncPP = std::function<void(
480 TDATA &, // cx
481 TDATA &, // xPrev
482 TDATA &, // crhs
483 TDATA &, // rhsIncPart,
484 const std::function<void()> &, // renewRhsIncPart
485 real, // ct
486 int // uPos
487 )>;
488 using FalphaLimSource = std::function<void(
489 TDATA &, // source
490 int // uPos
491 )>;
492
493 void StepPP(TDATA &x, TDATA &xinc, const Frhs &frhs, const Fdt &fdt, const Fsolve &fsolve,
494 int maxIter, const Fstop &fstop, const Fincrement &fincrement, real dt,
495 const FalphaLimSource &falphaLimSource,
496 const FresidualIncPP &fresidualIncPP)
497 {
498 index kCurrent = cnPrev + 1;
499 index prevSiz = kBDF - 1;
500 for (index iPrev = 0; iPrev < cnPrev; iPrev++)
501 assert(prevSiz && std::abs(dtPrevs[mod(iPrev + prevStart, prevSiz)] - dt) < dt * 1e-8);
502
503 xLast = x;
504 // x = xLast;
505 xIncPrev.setConstant(0.0);
506 int iter = 1;
507 for (; iter <= maxIter; iter++)
508 {
509 fdt(x, dTau, BDFCoefs(kCurrent - 1, 0), 0);
510
511 frhs(rhsbuf[0], x, dTau, iter, 1.0, 0);
512 fresidualIncPP(
513 x, xLast, rhsbuf[0], resInc,
514 [&]()
515 {
516 resInc.setConstant(0.0);
517 // resInc.addTo(x, -1. / dt);
518 resInc.addTo(xLast, (BDFCoefs(kCurrent - 1, 1) - 1) / dt);
519 // std::cout << "add " << BDFCoefs(kCurrent - 1, 1) << " " << "last" << std::endl;
520 if (prevSiz)
521 for (index iPrev = 0; iPrev < cnPrev; iPrev++)
522 {
523 // std::cout << "add " << BDFCoefs(kCurrent - 1, 2 + iPrev) <<" " << mod(iPrev + prevStart, prevSiz) << std::endl;
524 resInc.addTo(xPrevs[mod(iPrev + prevStart, prevSiz)], BDFCoefs(kCurrent - 1, 2 + iPrev) / dt);
525 }
526 falphaLimSource(resInc, 0); // non-rhs part of residual, fixed with alpha too
527 resInc.addTo(rhsbuf[0], BDFCoefs(kCurrent - 1, 0));
528 resInc *= dt; // so that equation is resInc == x - xLast
529 },
530 1.0,
531 0);
532
533 rhs.setConstant(0.0);
534 // rhsbuf .addTo(x, -1. / dt);
535 rhs.addTo(xLast, (BDFCoefs(kCurrent - 1, 1) - 1) / dt);
536 // std::cout << "add " << BDFCoefs(kCurrent - 1, 1) << " " << "last" << std::endl;
537 if (prevSiz)
538 for (index iPrev = 0; iPrev < cnPrev; iPrev++)
539 {
540 // std::cout << "add " << BDFCoefs(kCurrent - 1, 2 + iPrev) <<" " << mod(iPrev + prevStart, prevSiz) << std::endl;
541 rhs.addTo(xPrevs[mod(iPrev + prevStart, prevSiz)], BDFCoefs(kCurrent - 1, 2 + iPrev) / dt);
542 }
543 falphaLimSource(rhs, 0); // non-rhs part of residual, fixed with alpha too
544 rhs.addTo(xLast, 1 / dt);
545 resOther = rhs;
546 rhs.addTo(x, -1. / dt);
547 rhs.addTo(rhs, BDFCoefs(kCurrent - 1, 0));
548 fsolve(x, rhs, resOther, dTau, dt, BDFCoefs(kCurrent - 1, 0), xinc, iter, 1.0, 0);
549 //* xinc = (I/dtau-A*alphaDiag)\rhs
550
551 // std::cout << "BDF::\n";
552 // std::cout << kCurrent << " " << cnPrev<<" " << BDFCoefs(kCurrent - 1, 0) << std::endl;
553
554 // x += xinc;
555 fincrement(x, xinc, 1.0, 0);
556 // x.addTo(xIncPrev, -0.5);
557
558 xIncPrev = xinc;
559
560 if (fstop(iter, rhs, 1))
561 break;
562 }
563 if (iter > maxIter)
564 fstop(iter, rhs, 1);
565 if (prevSiz)
566 {
567 prevStart = mod(prevStart - 1, prevSiz);
568 // std::cout << dtPrevs.size() << " " << prevStart << std::endl;
569 xPrevs[prevStart] = xLast;
570 dtPrevs[prevStart] = dt;
571 cnPrev = std::min(cnPrev + 1, prevSiz);
572 }
573 }
574
575 virtual TDATA &getLatestRHS() override
576 {
577 return rhsbuf[0];
578 }
579
580 TDATA &getRHS(int i) override
581 {
582 return rhsbuf[0];
583 }
584
585 TDATA &getRES(int i) override
586 {
587 return rhs;
588 }
589
590 virtual ~ImplicitBDFDualTimeStep() = default;
591 };
592
593 template <class TDATA, class TDTAU>
594 const Eigen::Matrix<real, 4, 5> ImplicitBDFDualTimeStep<TDATA, TDTAU>::BDFCoefs{
595 {1. / 1., 1. / 1., std::nan("1"), std::nan("1"), std::nan("1")},
596 {2. / 3., 4. / 3., -1. / 3., std::nan("1"), std::nan("1")},
597 {6. / 11., 18. / 11., -9. / 11., 2. / 11., std::nan("1")},
598 {12. / 25., 48. / 25., -36. / 25., 16. / 25., -3. / 25.}};
599
600 template <class TDATA, class TDTAU>
601 class ImplicitVBDFDualTimeStep : public ImplicitDualTimeStep<TDATA, TDTAU>
602 {
603
604 public:
605 using tBase = ImplicitDualTimeStep<TDATA, TDTAU>;
606 using Frhs = typename tBase::Frhs;
607 using Fdt = typename tBase::Fdt;
608 using Fsolve = typename tBase::Fsolve;
609 using Fstop = typename tBase::Fstop;
610 using Fincrement = typename tBase::Fincrement;
611
612 TDTAU dTau;
613 std::vector<TDATA> xPrevs;
614 Eigen::VectorXd dtPrevs;
615 std::vector<TDATA> rhsbuf;
616 TDATA rhs;
617 TDATA resOther;
618 TDATA xLast;
619 TDATA xBase;
620 TDATA xIncPrev;
621 TDATA resInc;
622
623 index DOF;
624 index cnPrev;
625 index prevStart;
626 index kBDF;
627
628 private:
629 index kCurrent = 1;
630 index prevSiz = 0;
631 Eigen::Vector<real, Eigen::Dynamic> BDFCoefs;
632
633 public:
634 /**
635 * @brief mind that NDOF is the dof of dt
636 * finit(TDATA& data)
637 */
638 template <class Finit, class FinitDtau>
639 ImplicitVBDFDualTimeStep(
640 index NDOF, Finit &&finit = [](TDATA &) {}, FinitDtau &&finitDtau = [](TDTAU &) {},
641 index k = 2) : DOF(NDOF), cnPrev(0), prevStart(k - 2), kBDF(k)
642 {
643 assert(k > 0 && k <= 4);
644
645 xPrevs.resize(k - 1);
646 dtPrevs.resize(k - 1);
647 for (auto &i : xPrevs)
648 finit(i);
649 rhsbuf.resize(1);
650 finit(rhsbuf[0]);
651 finit(rhs);
652 finit(resOther);
653 finit(resInc);
654 finit(xLast);
655 finit(xBase);
656 finit(xIncPrev);
657 finitDtau(dTau);
658 DNDS_assert(k <= 2);
659 }
660
661 void VBDFFrontMatters(real dt)
662 {
663 kCurrent = cnPrev + 1;
664 prevSiz = kBDF - 1;
665 // for (index iPrev = 0; iPrev < cnPrev; iPrev++)
666 // assert(prevSiz && std::abs(dtPrevs[mod(iPrev + prevStart, prevSiz)] - dt) < dt * 1e-8);
667 BDFCoefs.resize(kCurrent + 1);
668
669 switch (kCurrent)
670 {
671 case 1:
672 BDFCoefs(0) = BDFCoefs(1) = 1;
673 break;
674 case 2:
675 {
676 real phi = dt / (dtPrevs[mod(0 + prevStart, prevSiz)] + dt);
677 real Rt = dt / dtPrevs[mod(0 + prevStart, prevSiz)];
678 BDFCoefs(0) = 1 / (1 + phi); // C
679 BDFCoefs(1) = 1 + Rt * phi / (1 + phi); // -A
680 BDFCoefs(2) = -Rt * phi / (1 + phi); // -B
681 }
682 break;
683
684 default:
685 DNDS_assert(false);
686 break;
687 }
688 }
689
690 /**
691 * @brief
692 * frhs(TDATA &rhs, TDATA &x)
693 * fdt(std::vector<real>& dTau)
694 * fsolve(TDATA &x, TDATA &rhs, std::vector<real>& dTau, real dt, real alphaDiag, TDATA &xinc)
695 * bool fstop(int iter, TDATA &xinc, int iInternal)
696 */
697 virtual void Step(TDATA &x, TDATA &xinc, const Frhs &frhs, const Fdt &fdt, const Fsolve &fsolve,
698 int maxIter, const Fstop &fstop, const Fincrement &fincrement, real dt) override
699 {
700 VBDFFrontMatters(dt);
701 xLast = x;
702 // x = xLast;
703 xIncPrev.setConstant(0.0);
704 int iter = 1;
705 for (; iter <= maxIter; iter++)
706 {
707 fdt(x, dTau, BDFCoefs(0), 0);
708
709 frhs(rhsbuf[0], x, dTau, iter, 1.0, 0);
710
711 rhs.setConstant(0.0);
712 rhs.addTo(xLast, BDFCoefs(1) / dt);
713 // std::cout << "add " << BDFCoefs(1) << " " << "last" << std::endl;
714 if (prevSiz)
715 for (index iPrev = 0; iPrev < cnPrev; iPrev++)
716 {
717 // std::cout << "add " << BDFCoefs(2 + iPrev) <<" " << mod(iPrev + prevStart, prevSiz) << std::endl;
718 rhs.addTo(xPrevs[mod(iPrev + prevStart, prevSiz)], BDFCoefs(2 + iPrev) / dt);
719 }
720 resOther = rhs;
721 rhs.addTo(x, -1. / dt);
722 rhs.addTo(rhsbuf[0], BDFCoefs(0));
723 fsolve(x, rhs, resOther, dTau, dt, BDFCoefs(0), xinc, iter, 1.0, 0);
724 //* xinc = (I/dtau-A*alphaDiag)\rhs
725
726 // std::cout << "BDF::\n";
727 // std::cout << kCurrent << " " << cnPrev<<" " << BDFCoefs(0) << std::endl;
728
729 // x += xinc;
730 fincrement(x, xinc, 1.0, 0);
731 // x.addTo(xIncPrev, -0.5);
732
733 xIncPrev = xinc;
734
735 if (fstop(iter, rhs, 1))
736 break;
737 }
738 if (iter > maxIter)
739 fstop(iter, rhs, 1);
740 if (prevSiz)
741 {
742 prevStart = mod(prevStart - 1, prevSiz);
743 // std::cout << dtPrevs.size() << " " << prevStart << std::endl;
744 xPrevs[prevStart] = xLast;
745 dtPrevs[prevStart] = dt;
746 cnPrev = std::min(cnPrev + 1, prevSiz);
747 }
748 }
749
750 using FresidualIncPP = std::function<void(
751 TDATA &, // cx
752 TDATA &, // xPrev
753 TDATA &, // crhs
754 TDATA &, // rhsIncPart,
755 const std::function<void()> &, // renewRhsIncPart
756 real, // ct
757 int // uPos
758 )>;
759 using FalphaLimSource = std::function<void(
760 TDATA &, // source
761 int // uPos
762 )>;
763 using FlimitDtBDF = std::function<real(
764 TDATA &, // base u,
765 TDATA & // to be limited incu
766 )>;
767
768 void LimitDt_StepPPV2(TDATA &xIn, const FlimitDtBDF &flimitDtBDF, real &dtOut, real maxIncrease = 2)
769 {
770 VBDFFrontMatters(dtOut); // using a wished dt
771 switch (kCurrent)
772 {
773 case 1:
774 // nothing
775 break;
776 case 2:
777 {
778 DNDS_assert(prevSiz == 1);
779 xLast = xIn;
780 xLast.addTo(xPrevs[mod(0 + prevStart, prevSiz)], -1.0);
781 xLast *= -BDFCoefs(2); // max value
782 real limitingV = flimitDtBDF(xIn, xLast);
783 // std::cout << fmt::format("limitingV {}, B {}", limitingV, -BDFCoefs(2)) << std::endl;
784 if (limitingV >= 1)
785 break;
786 DNDS_assert(limitingV >= 0);
787 real dtm1 = dtPrevs[mod(0 + prevStart, prevSiz)];
788 auto fB = [=](real dt)
789 {
790 real Rt = dt / dtm1;
791 real phi = dt / (dt + dtm1);
792 return Rt * phi / (1 + phi);
793 };
794 real targetB = limitingV * -BDFCoefs(2);
795 real dtNew = dtm1;
796 if (targetB >= fB(dtm1 * maxIncrease))
797 dtNew = dtm1 * maxIncrease; // max increase value
798 else
799 dtNew = Scalar::BisectSolveLower(fB, 0., dtm1 * maxIncrease, targetB * 0.99, 20);
800
801 dtOut = std::min(dtNew, dtOut);
802 }
803 break;
804
805 default:
806 DNDS_assert(false);
807 break;
808 }
809 }
810
811 void StepPP(TDATA &x, TDATA &xinc, const Frhs &frhs, const Fdt &fdt, const Fsolve &fsolve,
812 int maxIter, const Fstop &fstop, const Fincrement &fincrement, real dt,
813 const FalphaLimSource &falphaLimSource,
814 const FresidualIncPP &fresidualIncPP)
815 {
816 this->VBDFFrontMatters(dt);
817
818 xLast = x;
819 // x = xLast;
820 xIncPrev.setConstant(0.0);
821 int iter = 1;
822 for (; iter <= maxIter; iter++)
823 {
824 fdt(x, dTau, BDFCoefs(0), 0);
825
826 frhs(rhs, x, dTau, iter, 1.0, 0);
827
828 xBase.setConstant(0.0);
829 xBase.addTo(xLast, (BDFCoefs(1) - 0) / dt); // 0 because xBase includes xLast/dt part
830 if (prevSiz)
831 for (index iPrev = 0; iPrev < cnPrev; iPrev++)
832 {
833 // std::cout << "add " << BDFCoefs(2 + iPrev) <<" " << mod(iPrev + prevStart, prevSiz) << std::endl;
834 xBase.addTo(xPrevs[mod(iPrev + prevStart, prevSiz)], BDFCoefs(2 + iPrev) / dt);
835 }
836 fresidualIncPP(
837 x, xBase, rhsbuf[0], resInc,
838 [&]()
839 {
840 resInc.setConstant(0.0);
841 // ***** excluded with VBDF's ts limiting
842 // // resInc.addTo(x, -1. / dt);
843 // resInc.addTo(xLast, (BDFCoefs(1) - 1) / dt);
844 // // std::cout << "add " << BDFCoefs( 1) << " " << "last" << std::endl;
845 // if (prevSiz)
846 // for (index iPrev = 0; iPrev < cnPrev; iPrev++)
847 // {
848 // // std::cout << "add " << BDFCoefs(2 + iPrev) <<" " << mod(iPrev + prevStart, prevSiz) << std::endl;
849 // resInc.addTo(xPrevs[mod(iPrev + prevStart, prevSiz)], BDFCoefs(2 + iPrev) / dt);
850 // }
851 // falphaLimSource(resInc, 0); // non-rhs part of residual, fixed with alpha too
852 // *****
853 resInc.addTo(rhsbuf[0], BDFCoefs(0));
854 resInc *= dt; // so that equation is resInc == x - xLast
855 },
856 1.0,
857 0);
858
859 rhs = xBase;
860 // rhs.addTo(xLast, 1 / dt); // 0 because xBase includes xLast/dt part
861 resOther = rhs;
862 rhs.addTo(x, -1. / dt);
863 rhs.addTo(rhsbuf[0], BDFCoefs(0));
864
865 fsolve(x, rhs, resOther, dTau, dt, BDFCoefs(0), xinc, iter, 1.0, 0);
866 //* xinc = (I/dtau-A*alphaDiag)\rhs
867
868 // std::cout << "BDF::\n";
869 // std::cout << kCurrent << " " << cnPrev<<" " << BDFCoefs(0) << std::endl;
870
871 // x += xinc;
872 fincrement(x, xinc, 1.0, 0);
873 // x.addTo(xIncPrev, -0.5);
874
875 xIncPrev = xinc;
876
877 if (fstop(iter, rhs, 1))
878 break;
879 }
880 if (iter > maxIter)
881 fstop(iter, rhs, 1);
882 if (prevSiz)
883 {
884 prevStart = mod(prevStart - 1, prevSiz);
885 // std::cout << dtPrevs.size() << " " << prevStart << std::endl;
886 xPrevs[prevStart] = xLast;
887 dtPrevs[prevStart] = dt;
888 cnPrev = std::min(cnPrev + 1, prevSiz);
889 }
890 }
891
892 virtual TDATA &getLatestRHS() override
893 {
894 return rhsbuf[0];
895 }
896
897 TDATA &getRHS(int i) override
898 {
899 return rhsbuf[0];
900 }
901
902 TDATA &getRES(int i) override
903 {
904 return rhs;
905 }
906
907 virtual ~ImplicitVBDFDualTimeStep() = default;
908 };
909
910 /*******************************************************************************************/
911 /* */
912 /* */
913 /* */
914 /*******************************************************************************************/
915
916 template <class TDATA, class TDTAU>
917 class ImplicitHermite3SimpleJacobianDualStep : public ImplicitDualTimeStep<TDATA, TDTAU>
918 {
919 int hasLastEndPointR = 0;
920
921 public:
922 using tBase = ImplicitDualTimeStep<TDATA, TDTAU>;
923 using Frhs = typename tBase::Frhs;
924 using Fdt = typename tBase::Fdt;
925 using Fsolve = typename tBase::Fsolve;
926 using Fstop = typename tBase::Fstop;
927 using Fincrement = typename tBase::Fincrement;
928 using FsolveNest = std::function<void(
929 TDATA &, TDATA &, TDATA &,
930 TDTAU &, const std::vector<real> &,
931 real, real, TDATA &, int, int)>;
932
933 TDTAU dTau, dTauMid;
934 TDATA xMid, rhsMid, rhsFull, resOther;
935 std::vector<TDATA> rhsbuf;
936 TDATA xLast, xMG0, xMG;
937 TDATA xIncPrev;
938 index DOF;
939 index cnPrev;
940 TDATA xIncDamper;
941 TDATA xIncDamper2;
942
943 Eigen::Vector<real, 4> cInter;
944 Eigen::Vector<real, 3> wInteg;
945 int curSolveMethod = 0;
946 int nStartIter = 0;
947 real thetaM1 = 0.9146;
948 real thetaM2 = 0.0;
949 real thetaMMG = 1.0;
950 real coefIncMidMG = 1.0;
951 real alphaHM3 = 0.5;
952 int maskHM3 = 0;
953 int maskHM3Exe = 0;
954
955 int nMG = 0;
956
957 TDATA xPrev;
958 real dtPrev = 0;
959 int prevSize = 0;
960
961 /**
962 * @brief mind that NDOF is the dof of dt
963 * finit(TDATA& data)
964 */
965 template <class Finit, class FinitDtau>
966 ImplicitHermite3SimpleJacobianDualStep(
967 index NDOF, Finit &&finit = [](TDATA &) {}, FinitDtau &&finitDtau = [](TDTAU &) {},
968 real alpha = 0.55, int nCurSolveMethod = 0, int nnStartIter = 0,
969 real thetaM1n = 0.9146, real thetaM2n = 0.0, int mask = 0,
970 int nMGn = 4)
971 : DOF(NDOF),
972 cnPrev(0),
973 curSolveMethod(nCurSolveMethod),
974 nStartIter(nnStartIter),
975 thetaM1(thetaM1n),
976 thetaM2(thetaM2n),
977 alphaHM3(alpha),
978 maskHM3(mask),
979 nMG(nMGn)
980 {
981 rhsbuf.resize(3);
982 finit(rhsbuf[0]);
983 finit(rhsbuf[1]);
984 finit(rhsbuf[2]);
985 finit(xMid);
986 finit(rhsMid);
987 finit(rhsFull);
988 finit(resOther);
989 finit(xLast);
990 finit(xIncPrev);
991 finit(xMG0);
992 finit(xMG);
993 finit(xIncDamper);
994 finit(xIncDamper2);
995 finit(xPrev);
996 finitDtau(dTau);
997 finitDtau(dTauMid);
998
999 DNDS_assert_info(mask == 0 || mask == 1 || mask == 2, "mask not supported");
1000 SetCoefs(1);
1001 }
1002
1003 virtual void SetExtraParams(const nlohmann::ordered_json &j) override
1004 {
1005 for (auto &[k, v] : j.items())
1006 {
1007 if (k == "nMG")
1008 {
1009 DNDS_assert_info(v.is_number(), "need a number for nMG");
1010 nMG = v;
1011 }
1012 else if (k == "thetaMMG")
1013 {
1014 DNDS_assert_info(v.is_number(), "need a number for thetaMMG");
1015 thetaMMG = v;
1016 }
1017 else if (k == "coefIncMidMG")
1018 {
1019 DNDS_assert_info(v.is_number(), "need a number for coefIncMidMG");
1020 coefIncMidMG = v;
1021 }
1022 else
1023 {
1024 DNDS_assert_info(false, "no such key for HM3: " + k);
1025 }
1026 }
1027 }
1028
1029 void SetCoefs(real hR1 = 1)
1030 {
1031 real alpha = alphaHM3;
1032 int mask = maskHM3;
1033 maskHM3Exe = mask;
1034 assert(alpha > 0 && alpha < 1);
1035 cInter.setZero();
1036 switch (mask)
1037 {
1038 case 0: // U2R2
1039 {
1040 cInter[0] = (alpha * alpha) * -3.0 + (alpha * alpha * alpha) * 2.0 + 1.0;
1041 cInter[1] = (alpha * alpha) * 3.0 - (alpha * alpha * alpha) * 2.0;
1042 cInter[2] = alpha - (alpha * alpha) * 2.0 + alpha * alpha * alpha;
1043 cInter[3] = -alpha * alpha + alpha * alpha * alpha;
1044 }
1045 break;
1046 case 1: // U2R1
1047 {
1048 cInter[0] = alpha * -2.0 + alpha * alpha + 1.0;
1049 cInter[1] = alpha * 2.0 - alpha * alpha;
1050 cInter[2] = 0;
1051 cInter[3] = -alpha + alpha * alpha;
1052 }
1053 break;
1054 case 2: // U3R1
1055 {
1056 if (prevSize == 1 && (hR1 < 1e2 && hR1 > 1e-2))
1057 {
1058 // note that the meaning of cInter[2] is changed!
1059 cInter[2] = -(alpha * pow(alpha - 1.0, 2.0) * 1.0 / pow(hR1 + 1.0, 2.0)) / hR1;
1060 cInter[0] = ((alpha + hR1) * pow(alpha - 1.0, 2.0)) / hR1;
1061 cInter[1] = alpha * 1.0 / pow(hR1 + 1.0, 2.0) * (alpha * 3.0 + hR1 * 3.0 - alpha * (hR1 * hR1) - (alpha * alpha) * hR1 - (alpha * alpha) * 2.0 + (hR1 * hR1) * 2.0);
1062 cInter[3] = (alpha * (alpha + hR1) * (alpha - 1.0)) / (hR1 + 1.0);
1063 }
1064 else
1065 {
1066 maskHM3Exe = 1;
1067 alpha = 0.25;
1068 cInter[0] = alpha * -2.0 + alpha * alpha + 1.0;
1069 cInter[1] = alpha * 2.0 - alpha * alpha;
1070 cInter[2] = 0;
1071 cInter[3] = -alpha + alpha * alpha;
1072 }
1073 }
1074 break;
1075 default:
1076 DNDS_assert(false);
1077 break;
1078 }
1079
1080 wInteg[0] = (-1.0 / 6.0) / alpha + 1.0 / 2.0;
1081 wInteg[1] = (-1.0 / 6.0) / (alpha * (alpha - 1.0));
1082 wInteg[2] = 1.0 / (alpha * 6.0 - 6.0) + 1.0 / 2.0;
1083 }
1084
1085 /**
1086 * @brief
1087 * frhs(TDATA &rhs, TDATA &x)
1088 * fdt(TDTAU& dTau)
1089 * fsolve(TDATA &x, TDATA &rhs, TDTAU& dTau, real dt, real alphaDiag, TDATA &xinc)
1090 * bool fstop(int iter, TDATA &xinc, int iInternal)
1091 */
1092 virtual void Step(TDATA &x, TDATA &xinc, const Frhs &frhs, const Fdt &fdt, const Fsolve &fsolve,
1093 int maxIter, const Fstop &fstop, const Fincrement &fincrement, real dt) override
1094 {
1095 SetCoefs(dtPrev / (dt + verySmallReal));
1096 xLast = x;
1097 xMid = x;
1098
1099 if (hasLastEndPointR)
1100 {
1101 rhsbuf[0] = rhsbuf[1];
1102 // dTau = dTau here
1103 }
1104 else
1105 {
1106 fdt(xLast, dTau, 1.0, 0);
1107 frhs(rhsbuf[0], xLast, dTau, INT_MAX, 0.0, 0);
1108 }
1109 rhsbuf[1] = rhsbuf[0];
1110 rhsbuf[2] = rhsbuf[0];
1111 dTauMid = dTau;
1112
1113 xIncPrev.setConstant(0.0);
1114 int iter = 1;
1115
1116 int method = curSolveMethod;
1117 bool stepIsRealU3R1 = prevSize >= 1 && maskHM3 == 2;
1118
1119 for (; iter <= maxIter; iter++)
1120 {
1121
1122 if (iter < nStartIter)
1123 {
1124 fdt(x, dTau, 1.0, 0);
1125 frhs(rhsbuf[1], x, dTau, iter, 1.0, 0);
1126 rhsFull = xLast;
1127 rhsFull *= 1.0 / dt;
1128 resOther = rhsFull;
1129 rhsFull.addTo(x, -1. / dt);
1130 rhsFull += rhsbuf[1];
1131 fsolve(x, rhsFull, resOther, dTau, dt, 1.0, xinc, iter, 1.0, 0);
1132 }
1133 else
1134 {
1135 if (method == 0)
1136 {
1137 rhsMid.setConstant(0.0);
1138 real thetaCur = thetaM1;
1139 {
1140 if (maskHM3Exe == 1 && maskHM3 == 2)
1141 thetaCur = 1; // for U2R1 filling first step of U3R1
1142 rhsMid.addTo(xLast, (cInter(0) + thetaCur) / dt);
1143 rhsMid.addTo(x, (cInter(1) - thetaCur) / dt);
1144 if (stepIsRealU3R1)
1145 rhsMid.addTo(xPrev, cInter(2) / dt);
1146 rhsMid.addTo(rhsbuf[0], (stepIsRealU3R1 ? 0.0 : cInter(2)) + thetaCur * wInteg(0));
1147 rhsMid.addTo(rhsbuf[1], cInter(3) + thetaCur * wInteg(2));
1148 resOther = rhsMid;
1149 rhsMid.addTo(xMid, -1. / dt);
1150 rhsMid.addTo(rhsbuf[2], thetaCur * wInteg(1));
1151
1152 fsolve(xMid, rhsMid, resOther, dTauMid, dt, std::abs(thetaCur * wInteg(1)), xinc, iter, alphaHM3, 1);
1153 }
1154
1155 fincrement(xMid, xinc, 1.0, 1);
1156 fdt(xMid, dTauMid, 1.0, 1);
1157 frhs(rhsbuf[2], xMid, dTauMid, iter, alphaHM3, 1);
1158
1159 rhsFull.setConstant(0.0);
1160 real theta2Scale = 1 - thetaM2 * cInter(1);
1161 DNDS_assert_info(theta2Scale > 0, fmt::format("thetaM2 {}, cInter(1) {}", thetaM2, cInter(1)));
1162 theta2Scale = 1. / theta2Scale;
1163 bool theta2LimitSitu = thetaM2 < -1e5;
1164 real xLastCoef = !theta2LimitSitu
1165 ? (theta2Scale * (1. + thetaM2 * cInter(0)) / dt)
1166 : (-cInter(0) / cInter(1) / dt);
1167 real xMidCoef = !theta2LimitSitu
1168 ? (-theta2Scale * thetaM2 * 1. / dt)
1169 : (1 / cInter(1) * 1. / dt);
1170 real xPrevCoef = !theta2LimitSitu
1171 ? (theta2Scale * thetaM2 * cInter(2) / dt)
1172 : (-1. / cInter(1) * cInter(2) / dt);
1173 real rLastCoef = theta2Scale * (wInteg(0) + thetaM2 * (stepIsRealU3R1 ? 0.0 : cInter(2)));
1174 if (theta2LimitSitu)
1175 rLastCoef = -1. / cInter(1) * (stepIsRealU3R1 ? 0.0 : cInter(2));
1176 real rMidCoef = theta2LimitSitu ? 0.0 : theta2Scale * wInteg(1);
1177 real rCoef = !theta2LimitSitu ? (theta2Scale * (wInteg(2) + thetaM2 * cInter(3)))
1178 : (-1. / cInter(1) * cInter(3));
1179 // std::cout << xLastCoef << " "
1180 // << xMidCoef << " "
1181 // << xPrevCoef << " "
1182 // << rLastCoef << " "
1183 // << rMidCoef << " "
1184 // << rCoef << " " << std::endl;
1185 {
1186 rhsFull.addTo(xLast, xLastCoef);
1187 if (thetaM2)
1188 {
1189 // rhsFull.addTo(x, thetaM2 * cInter(1) / dt); // need to add to diag part!!
1190 rhsFull.addTo(xMid, xMidCoef);
1191 if (stepIsRealU3R1)
1192 rhsFull.addTo(xPrev, xPrevCoef);
1193 }
1194
1195 rhsFull.addTo(rhsbuf[0], rLastCoef);
1196 rhsFull.addTo(rhsbuf[2], rMidCoef);
1197
1198 resOther = rhsFull;
1199 rhsFull.addTo(x, -1. / dt);
1200 rhsFull.addTo(rhsbuf[1], rCoef);
1201
1202 fsolve(x, rhsFull, resOther, dTau, dt, rCoef, xinc, iter, 1.0, 0);
1203 }
1204
1205 fincrement(x, xinc, 1.0, 0);
1206 fdt(x, dTau, 1.0, 0);
1207 frhs(rhsbuf[1], x, dTau, iter, 1.0, 0);
1208
1209 // residual #1 + #0 * thetaMMG
1210 {
1211 rhsFull.setConstant(0.0);
1212 rhsFull.addTo(xLast, (1. + thetaMMG * cInter(0)) / dt);
1213 rhsFull.addTo(x, (thetaMMG * cInter(1)) / dt);
1214 rhsFull.addTo(xMid, -thetaMMG * 1. / dt);
1215 if (stepIsRealU3R1)
1216 rhsFull.addTo(xPrev, thetaMMG * cInter(2) / dt);
1217 rhsFull.addTo(rhsbuf[0], wInteg(0) + thetaMMG * (stepIsRealU3R1 ? 0.0 : cInter(2)));
1218 rhsFull.addTo(rhsbuf[2], wInteg(1));
1219 rhsFull.addTo(x, -1. / dt);
1220 rhsFull.addTo(rhsbuf[1], wInteg(2) + thetaMMG * cInter(3));
1221 }
1222 // std::cout << thetaMMG << std::endl;
1223
1224 // * START pMG part
1225 if (nMG)
1226 {
1227 {
1228 // resOther.setConstant(0.0);
1229 // resOther.addTo(xLast, (1. + thetaMMG * cInter(0)) / dt);
1230 // resOther.addTo(x, (thetaMMG * cInter(1)) / dt);
1231 // resOther.addTo(xMid, -thetaMMG * 1. / dt);
1232 // if (stepIsRealU3R1)
1233 // resOther.addTo(xPrev, thetaMMG * cInter(2) / dt);
1234 // resOther.addTo(rhsbuf[0], wInteg(0) + thetaMMG * (stepIsRealU3R1 ? 0.0 : cInter(2)));
1235 // resOther.addTo(rhsbuf[2], wInteg(1));
1236 // resOther.addTo(x, -1. / dt);
1237 // resOther.addTo(rhsbuf[1], wInteg(2) + thetaMMG * cInter(3));
1238 }
1239 resOther = rhsFull;
1240 xMG0 = x;
1241 xMG0 *= 0.5;
1242 xMG0.addTo(xMid, 0.5);
1243 xMG = xMG0;
1244 // F_trapz == (xLast - x) / dt + (rLast + r) * 0.5
1245 // F_trapz - F_trapz_0 == (x0 - x) / dt + (r - r_0) * 0.5
1246 real mgTrapzAlpha = 1.0;
1247 for (int iMG = 1; iMG <= nMG; iMG++)
1248 {
1249 // use upos == 2 for lower order spacial frhs
1250 fdt(xMG, dTau, 1.0, /*upos=*/2);
1251 frhs(rhsbuf[1], xMG, dTau, iter, 1.0, /*upos=*/2); // pos = 0 for original RHS
1252 if (iMG == 1 and /*1 smoother shortcut*/ nMG > 1)
1253 {
1254 resOther.addTo(rhsbuf[1], -mgTrapzAlpha);
1255 resOther.addTo(xMG, 1. / dt);
1256 }
1257 rhsMid = resOther;
1258 if (/*1 smoother shortcut*/ nMG > 1)
1259 {
1260 rhsMid.addTo(rhsbuf[1], mgTrapzAlpha);
1261 rhsMid.addTo(xMG, -1. / dt);
1262 }
1263
1264 fsolve(xMG, rhsMid, resOther, dTau, dt, mgTrapzAlpha, xinc, iter, 1.0, /*upos=*/2);
1265 fincrement(xMG, xinc, 1.0, /*upos=*/2);
1266 }
1267 xinc = xMG;
1268 xinc.addTo(xMG0, -1.0);
1269 fincrement(x, xinc, 1.0, 0);
1270 if (coefIncMidMG)
1271 fincrement(xMid, xinc, coefIncMidMG, 1);
1272 fdt(xMid, dTauMid, 1.0, 1);
1273 frhs(rhsbuf[2], xMid, dTauMid, iter, alphaHM3, 1);
1274
1275 fdt(x, dTau, 1.0, 0);
1276 frhs(rhsbuf[1], x, dTau, iter, 1.0, 0);
1277 }
1278 }
1279 else
1280 {
1281 DNDS_assert(false);
1282 }
1283 }
1284
1285 xIncPrev = xinc;
1286
1287 if (fstop(iter, rhsFull, 1))
1288 if (iter >= nStartIter)
1289 break;
1290 }
1291 if (iter > maxIter)
1292 fstop(iter, rhsFull, 1);
1293
1294 hasLastEndPointR = 1;
1295
1296 prevSize = 1;
1297 xPrev = xLast;
1298 dtPrev = dt;
1299 }
1300
1301 void StepNested(TDATA &x, TDATA &xinc, const Frhs &frhs, const Fdt &fdt, const Fsolve &fsolve, const FsolveNest &fsolveN,
1302 int maxIter, const Fstop &fstop, const Fincrement &fincrement, real dt)
1303 {
1304 xLast = x;
1305 frhs(rhsbuf[0], x, dTau, 0, 1.0, 0);
1306
1307 xIncPrev.setConstant(0.0);
1308 int iter = 1;
1309 for (; iter <= maxIter; iter++)
1310 {
1311
1312 if (iter < nStartIter)
1313 {
1314 fdt(x, dTau, 1.0, 0);
1315 frhs(rhsbuf[1], x, dTau, iter, 1.0, 0);
1316 rhsFull = xLast;
1317 rhsFull *= 1.0 / dt;
1318 resOther = rhsFull;
1319 rhsFull.addTo(x, -1. / dt);
1320 rhsFull += rhsbuf[1];
1321 fsolve(x, rhsFull, resOther, dTau, dt, 1.0, xinc, iter, 1.0, 0);
1322 }
1323 else
1324 {
1325 fdt(x, dTau, 1.0, 0);
1326
1327 frhs(rhsbuf[1], x, dTau, iter, 1.0, 0);
1328 xMid.setConstant(0.0);
1329 xMid.addTo(xLast, cInter[0]);
1330 xMid.addTo(x, cInter[1]);
1331 fincrement(xMid, rhsbuf[0], cInter[2] * dt, 1);
1332 fincrement(xMid, rhsbuf[1], cInter[3] * dt, 1);
1333 frhs(rhsMid, xMid, dTau, iter, 1.0, 1);
1334 rhsFull.setConstant(0.0);
1335 rhsFull.addTo(rhsbuf[0], wInteg[0]);
1336 rhsFull.addTo(rhsMid, wInteg[1]);
1337 rhsFull.addTo(rhsbuf[1], wInteg[2]);
1338 rhsFull.addTo(x, -1. / dt);
1339 rhsFull.addTo(xLast, 1. / dt);
1340
1341 {
1342 // damping
1343 // xIncDamper = xLast;
1344 // xIncDamper.addTo(x, -1.);
1345 // xIncDamper.setAbs();
1346 // xIncDamper += 1e-100;
1347 // xIncDamper2 = xIncPrev;
1348 // xIncDamper2.setAbs();
1349 // xIncDamper2 += xIncDamper;
1350 // xIncDamper /= xIncDamper2;
1351 // rhsFull *= xIncDamper;
1352 }
1353 {
1354 // fdt(x, dTau, 1.0); // TODO: use "update spectral radius" procedure? or force update in fsolve
1355 // for(auto &v : dTau)
1356 // v *= -cInter[3] / cInter[1];
1357 // fsolve(x, rhsFull, dTau, -dt * cInter[3] / cInter[1],
1358 // 1.0, xinc, iter);
1359 // rhsFull = xinc;
1360 // fdt(xMid, dTau, 1.0);
1361 // for (auto &v : dTau)
1362 // v = veryLargeReal;
1363 // fsolve(xMid, rhsFull, dTau, -dt * cInter[3] * wInteg[1] / wInteg[2],
1364 // 1.0, xinc, iter);
1365 // xinc *= -1. / (2 * dt * cInter[3] * wInteg[1]);
1366 fsolveN(x, xMid, rhsFull, dTau,
1367 std::vector<real>{
1368 // cInter[1] * wInteg[2] / (cInter[3] * dt),
1369 // -(cInter[3] * wInteg[1]),
1370 // -cInter[3] * wInteg[1] / wInteg[2],
1371 // -cInter[3] / cInter[1],
1372 0. / dt,
1373 1,
1374 1,
1375 1,
1376 },
1377 dt, 1.0, xinc, iter, 0);
1378 }
1379 }
1380
1381 //** xinc = (I/dtau-A*alphaDiag)\rhs
1382
1383 // x += xinc;
1384 {
1385 // xIncDamper = xLast;
1386 // xIncDamper.addTo(x, -1.);
1387 // xIncDamper.setAbs();
1388 // xIncDamper += 1e-100;
1389 // xIncDamper2 = xIncPrev;
1390 // xIncDamper2.setAbs();
1391 // xIncDamper2 += xIncDamper;
1392 // xIncDamper /= xIncDamper2;
1393 // xinc *= xIncDamper;
1394 }
1395
1396 fincrement(x, xinc, 1.0, 0);
1397 // x.addTo(xIncPrev, -0.5);
1398
1399 xIncPrev = xinc;
1400
1401 if (fstop(iter, xinc, 1))
1402 if (iter >= nStartIter)
1403 break;
1404 }
1405 if (iter > maxIter)
1406 fstop(iter, xinc, 1);
1407 }
1408
1409 virtual TDATA &getLatestRHS() override
1410 {
1411 return rhsbuf[1];
1412 }
1413
1414 TDATA &getRHS(int i) override
1415 {
1416 // note that 0 for tn, 1 for tn+1, 2 or tn+c2
1417 DNDS_assert(i >= 0 && i < rhsbuf.size());
1418 return rhsbuf[i];
1419 }
1420
1421 TDATA &getRES(int i) override
1422 {
1423 DNDS_assert(i >= 0 && i < 2);
1424 // todo: enable on-demanc calculation
1425 return i == 0 ? rhsFull : rhsMid;
1426 }
1427
1428 virtual ~ImplicitHermite3SimpleJacobianDualStep() = default;
1429 };
1430
1431 /*******************************************************************************************/
1432 /* */
1433 /* */
1434 /* */
1435 /*******************************************************************************************/
1436
1437 template <class TDATA, class TDTAU>
1438 class ExplicitSSPRK3TimeStepAsImplicitDualTimeStep : public ImplicitDualTimeStep<TDATA, TDTAU>
1439 {
1440 public:
1441 using tBase = ImplicitDualTimeStep<TDATA, TDTAU>;
1442 using Frhs = typename tBase::Frhs;
1443 using Fdt = typename tBase::Fdt;
1444 using Fsolve = typename tBase::Fsolve;
1445 using Fstop = typename tBase::Fstop;
1446 using Fincrement = typename tBase::Fincrement;
1447
1448 TDTAU dTau;
1449 std::vector<TDATA> rhsbuf;
1450 TDATA rhs;
1451 TDATA xLast;
1452 TDATA xInc;
1453 index DOF;
1454 bool localDtStepping{false};
1455
1456 template <class Finit, class FinitDtau>
1457 ExplicitSSPRK3TimeStepAsImplicitDualTimeStep(
1458 index NDOF, Finit &&finit = [](TDATA &) {}, FinitDtau &&finitDtau = [](TDTAU &) {},
1459 bool nLocalDtStepping = false)
1460 : DOF(NDOF), localDtStepping(nLocalDtStepping)
1461 {
1462
1463 rhsbuf.resize(3);
1464 for (auto &i : rhsbuf)
1465 finit(i);
1466 finit(rhs);
1467 finit(xLast);
1468 finit(xInc);
1469 finitDtau(dTau);
1470 }
1471
1472 /*!
1473
1474
1475 @brief fsolve, maxIter, fstop are omitted here
1476 */
1477 virtual void Step(TDATA &x, TDATA &xinc, const Frhs &frhs, const Fdt &fdt, const Fsolve &fsolve,
1478 int maxIter, const Fstop &fstop, const Fincrement &fincrement, real dt) override
1479 {
1480
1481 fdt(x, dTau, 1.0, 0); // always gets dTau for CFL evaluation
1482 xLast = x;
1483 // MPI::Barrier(MPI_COMM_WORLD);
1484 // std::cout << "fucked" << std::endl;
1485
1486 frhs(rhs, x, dTau, 1, 0.5, 0);
1487 rhsbuf[0] = rhs;
1488 if (localDtStepping)
1489 rhs *= dTau;
1490 else
1491 rhs *= dt;
1492
1493 // x += rhs;
1494 fincrement(x, rhs, 1.0, 0);
1495
1496 frhs(rhs, x, dTau, 1, 1, 0);
1497 rhsbuf[1] = rhs;
1498 if (localDtStepping)
1499 rhs *= dTau;
1500 else
1501 rhs *= dt;
1502 x *= 0.25;
1503 x.addTo(xLast, 0.75);
1504 // x.addTo(rhs, 0.25);
1505 fincrement(x, rhs, 0.25, 0);
1506
1507 frhs(rhs, x, dTau, 1, 0.25, 0);
1508 rhsbuf[2] = rhs;
1509 if (localDtStepping)
1510 rhs *= dTau;
1511 else
1512 rhs *= dt;
1513 x *= 2. / 3.;
1514 x.addTo(xLast, 1. / 3.);
1515 // x.addTo(rhs, 2. / 3.);
1516 fincrement(x, rhs, 2. / 3., 0);
1517 }
1518
1519 virtual TDATA &getLatestRHS() override
1520 {
1521 return rhsbuf[0];
1522 }
1523
1524 TDATA &getRHS(int i) override
1525 {
1526 DNDS_assert(i >= 0 && i < rhsbuf.size());
1527 return rhsbuf[i];
1528 }
1529
1530 TDATA &getRES(int i) override
1531 {
1532 // todo: return a zero (for this is explicit)
1533 return rhs;
1534 }
1535 };
1536}
Defines.hpp
Core type aliases, constants, and metaprogramming utilities for the DNDS framework.
DNDS_assert_info
#define DNDS_assert_info(expr, info)
Debug-only assertion with an extra std::string info message.
Definition Errors.hpp:113
DNDS_assert
#define DNDS_assert(expr)
Debug-only assertion (compiled out when DNDS_NDEBUG is defined). Prints the expression + file/line + ...
Definition Errors.hpp:108
JsonUtil.hpp
JSON-to-Eigen conversion utilities and nlohmann_json helper macros.
MPI.hpp
MPI wrappers: MPIInfo, collective operations, type mapping, CommStrategy.
_zeta
#define _zeta
DNDS::ODE::ExplicitSSPRK3TimeStepAsImplicitDualTimeStep::ExplicitSSPRK3TimeStepAsImplicitDualTimeStep
ExplicitSSPRK3TimeStepAsImplicitDualTimeStep(index NDOF, Finit &&finit=[](TDATA &) {}, FinitDtau &&finitDtau=[](TDTAU &) {}, bool nLocalDtStepping=false)
Definition ODE.hpp:1457
DNDS::ODE::ExplicitSSPRK3TimeStepAsImplicitDualTimeStep::Step
virtual void Step(TDATA &x, TDATA &xinc, const Frhs &frhs, const Fdt &fdt, const Fsolve &fsolve, int maxIter, const Fstop &fstop, const Fincrement &fincrement, real dt) override
fsolve, maxIter, fstop are omitted here
Definition ODE.hpp:1477
DNDS::ODE::ImplicitBDFDualTimeStep::getLatestRHS
virtual TDATA & getLatestRHS() override
Definition ODE.hpp:575
DNDS::ODE::ImplicitBDFDualTimeStep::ImplicitBDFDualTimeStep
ImplicitBDFDualTimeStep(index NDOF, Finit &&finit=[](TDATA &) {}, FinitDtau &&finitDtau=[](TDTAU &) {}, index k=2)
mind that NDOF is the dof of dt finit(TDATA& data)
Definition ODE.hpp:396
DNDS::ODE::ImplicitBDFDualTimeStep::Fsolve
typename tBase::Fsolve Fsolve
Definition ODE.hpp:372
DNDS::ODE::ImplicitBDFDualTimeStep::StepPP
void StepPP(TDATA &x, TDATA &xinc, const Frhs &frhs, const Fdt &fdt, const Fsolve &fsolve, int maxIter, const Fstop &fstop, const Fincrement &fincrement, real dt, const FalphaLimSource &falphaLimSource, const FresidualIncPP &fresidualIncPP)
Definition ODE.hpp:493
DNDS::ODE::ImplicitBDFDualTimeStep::Fdt
typename tBase::Fdt Fdt
Definition ODE.hpp:371
DNDS::ODE::ImplicitBDFDualTimeStep::Step
virtual void Step(TDATA &x, TDATA &xinc, const Frhs &frhs, const Fdt &fdt, const Fsolve &fsolve, int maxIter, const Fstop &fstop, const Fincrement &fincrement, real dt) override
frhs(TDATA &rhs, TDATA &x) fdt(TDTAU& dTau) fsolve(TDATA &x, TDATA &rhs, TDTAU& dTau,...
Definition ODE.hpp:423
DNDS::ODE::ImplicitBDFDualTimeStep::FresidualIncPP
std::function< void(TDATA &, TDATA &, TDATA &, TDATA &, const std::function< void()> &, real, int)> FresidualIncPP
Definition ODE.hpp:487
DNDS::ODE::ImplicitBDFDualTimeStep::Frhs
typename tBase::Frhs Frhs
Definition ODE.hpp:370
DNDS::ODE::ImplicitBDFDualTimeStep::Fincrement
typename tBase::Fincrement Fincrement
Definition ODE.hpp:374
DNDS::ODE::ImplicitBDFDualTimeStep::getRES
TDATA & getRES(int i) override
Definition ODE.hpp:585
DNDS::ODE::ImplicitBDFDualTimeStep::xPrevs
std::vector< TDATA > xPrevs
Definition ODE.hpp:377
DNDS::ODE::ImplicitBDFDualTimeStep::FalphaLimSource
std::function< void(TDATA &, int)> FalphaLimSource
Definition ODE.hpp:491
DNDS::ODE::ImplicitBDFDualTimeStep::Fstop
typename tBase::Fstop Fstop
Definition ODE.hpp:373
DNDS::ODE::ImplicitBDFDualTimeStep::~ImplicitBDFDualTimeStep
virtual ~ImplicitBDFDualTimeStep()=default
DNDS::ODE::ImplicitBDFDualTimeStep::rhsbuf
std::vector< TDATA > rhsbuf
Definition ODE.hpp:379
DNDS::ODE::ImplicitBDFDualTimeStep::getRHS
TDATA & getRHS(int i) override
Definition ODE.hpp:580
DNDS::ODE::ImplicitDualTimeStep::Frhs
std::function< void(TDATA &, TDATA &, TDTAU &, int, real, int)> Frhs
Definition ODE.hpp:13
DNDS::ODE::ImplicitDualTimeStep::getLatestRHS
virtual TDATA & getLatestRHS()=0
DNDS::ODE::ImplicitDualTimeStep::~ImplicitDualTimeStep
virtual ~ImplicitDualTimeStep()=default
DNDS::ODE::ImplicitDualTimeStep::Fsolve
std::function< void(TDATA &, TDATA &, TDATA &, TDTAU &, real, real, TDATA &, int, real, int)> Fsolve
Definition ODE.hpp:16
DNDS::ODE::ImplicitDualTimeStep::getRHS
virtual TDATA & getRHS(int i)=0
DNDS::ODE::ImplicitDualTimeStep::SetExtraParams
virtual void SetExtraParams(const nlohmann::ordered_json &j)
Definition ODE.hpp:30
DNDS::ODE::ImplicitDualTimeStep::Fdt
std::function< void(TDATA &, TDTAU &, real, int)> Fdt
Definition ODE.hpp:14
DNDS::ODE::ImplicitDualTimeStep::Step
virtual void Step(TDATA &x, TDATA &xinc, const Frhs &frhs, const Fdt &fdt, const Fsolve &fsolve, int maxIter, const Fstop &fstop, const Fincrement &fincrement, real dt)=0
DNDS::ODE::ImplicitDualTimeStep::getRES
virtual TDATA & getRES(int i)=0
DNDS::ODE::ImplicitDualTimeStep::Fincrement
std::function< void(TDATA &, TDATA &, real, int)> Fincrement
Definition ODE.hpp:18
DNDS::ODE::ImplicitDualTimeStep::Fstop
std::function< bool(int, TDATA &, int)> Fstop
Definition ODE.hpp:17
DNDS::ODE::ImplicitEulerDualTimeStep::rhsbuf
std::vector< TDATA > rhsbuf
Definition ODE.hpp:49
DNDS::ODE::ImplicitEulerDualTimeStep::Fsolve
typename tBase::Fsolve Fsolve
Definition ODE.hpp:44
DNDS::ODE::ImplicitEulerDualTimeStep::Fstop
typename tBase::Fstop Fstop
Definition ODE.hpp:45
DNDS::ODE::ImplicitEulerDualTimeStep::getRHS
TDATA & getRHS(int i) override
Definition ODE.hpp:111
DNDS::ODE::ImplicitEulerDualTimeStep::getRES
TDATA & getRES(int i) override
Definition ODE.hpp:116
DNDS::ODE::ImplicitEulerDualTimeStep::Step
virtual void Step(TDATA &x, TDATA &xinc, const Frhs &frhs, const Fdt &fdt, const Fsolve &fsolve, int maxIter, const Fstop &fstop, const Fincrement &fincrement, real dt) override
frhs(TDATA &rhs, TDATA &x) fdt(TDTAU& dTau) fsolve(TDATA &x, TDATA &rhs, TDTAU& dTau,...
Definition ODE.hpp:81
DNDS::ODE::ImplicitEulerDualTimeStep::~ImplicitEulerDualTimeStep
virtual ~ImplicitEulerDualTimeStep()=default
DNDS::ODE::ImplicitEulerDualTimeStep::getLatestRHS
virtual TDATA & getLatestRHS() override
Definition ODE.hpp:106
DNDS::ODE::ImplicitEulerDualTimeStep::ImplicitEulerDualTimeStep
ImplicitEulerDualTimeStep(index NDOF, Finit &&finit=[](TDATA &) {}, FinitDtau &&finitDtau=[](TDTAU &) {})
mind that NDOF is the dof of dt finit(TDATA& data)
Definition ODE.hpp:61
DNDS::ODE::ImplicitEulerDualTimeStep::Frhs
typename tBase::Frhs Frhs
Definition ODE.hpp:42
DNDS::ODE::ImplicitEulerDualTimeStep::Fdt
typename tBase::Fdt Fdt
Definition ODE.hpp:43
DNDS::ODE::ImplicitEulerDualTimeStep::Fincrement
typename tBase::Fincrement Fincrement
Definition ODE.hpp:46
DNDS::ODE::ImplicitHermite3SimpleJacobianDualStep::StepNested
void StepNested(TDATA &x, TDATA &xinc, const Frhs &frhs, const Fdt &fdt, const Fsolve &fsolve, const FsolveNest &fsolveN, int maxIter, const Fstop &fstop, const Fincrement &fincrement, real dt)
Definition ODE.hpp:1301
DNDS::ODE::ImplicitHermite3SimpleJacobianDualStep::ImplicitHermite3SimpleJacobianDualStep
ImplicitHermite3SimpleJacobianDualStep(index NDOF, Finit &&finit=[](TDATA &) {}, FinitDtau &&finitDtau=[](TDTAU &) {}, real alpha=0.55, int nCurSolveMethod=0, int nnStartIter=0, real thetaM1n=0.9146, real thetaM2n=0.0, int mask=0, int nMGn=4)
mind that NDOF is the dof of dt finit(TDATA& data)
Definition ODE.hpp:966
DNDS::ODE::ImplicitHermite3SimpleJacobianDualStep::SetExtraParams
virtual void SetExtraParams(const nlohmann::ordered_json &j) override
Definition ODE.hpp:1003
DNDS::ODE::ImplicitHermite3SimpleJacobianDualStep::Fincrement
typename tBase::Fincrement Fincrement
Definition ODE.hpp:927
DNDS::ODE::ImplicitHermite3SimpleJacobianDualStep::getLatestRHS
virtual TDATA & getLatestRHS() override
Definition ODE.hpp:1409
DNDS::ODE::ImplicitHermite3SimpleJacobianDualStep::FsolveNest
std::function< void(TDATA &, TDATA &, TDATA &, TDTAU &, const std::vector< real > &, real, real, TDATA &, int, int)> FsolveNest
Definition ODE.hpp:931
DNDS::ODE::ImplicitHermite3SimpleJacobianDualStep::Step
virtual void Step(TDATA &x, TDATA &xinc, const Frhs &frhs, const Fdt &fdt, const Fsolve &fsolve, int maxIter, const Fstop &fstop, const Fincrement &fincrement, real dt) override
frhs(TDATA &rhs, TDATA &x) fdt(TDTAU& dTau) fsolve(TDATA &x, TDATA &rhs, TDTAU& dTau,...
Definition ODE.hpp:1092
DNDS::ODE::ImplicitSDIRK4DualTimeStep::Fincrement
typename tBase::Fincrement Fincrement
Definition ODE.hpp:141
DNDS::ODE::ImplicitSDIRK4DualTimeStep::getLatestRHS
virtual TDATA & getLatestRHS() override
Definition ODE.hpp:342
DNDS::ODE::ImplicitSDIRK4DualTimeStep::getRES
TDATA & getRES(int i) override
Definition ODE.hpp:353
DNDS::ODE::ImplicitSDIRK4DualTimeStep::getRHS
TDATA & getRHS(int i) override
Definition ODE.hpp:347
DNDS::ODE::ImplicitSDIRK4DualTimeStep::ImplicitSDIRK4DualTimeStep
ImplicitSDIRK4DualTimeStep(index NDOF, Finit &&finit=[](TDATA &) {}, FinitDtau &&finitDtau=[](TDTAU &) {}, int schemeCode=0)
mind that NDOF is the dof of dt finit(TDATA& data)
Definition ODE.hpp:156
DNDS::ODE::ImplicitSDIRK4DualTimeStep::~ImplicitSDIRK4DualTimeStep
virtual ~ImplicitSDIRK4DualTimeStep()=default
DNDS::ODE::ImplicitSDIRK4DualTimeStep::Fsolve
typename tBase::Fsolve Fsolve
Definition ODE.hpp:139
DNDS::ODE::ImplicitSDIRK4DualTimeStep::Step
virtual void Step(TDATA &x, TDATA &xinc, const Frhs &frhs, const Fdt &fdt, const Fsolve &fsolve, int maxIter, const Fstop &fstop, const Fincrement &fincrement, real dt) override
frhs(TDATA &rhs, TDATA &x) fdt(TDTAU& dTau) fsolve(TDATA &x, TDATA &rhs, TDTAU& dTau,...
Definition ODE.hpp:271
DNDS::ODE::ImplicitSDIRK4DualTimeStep::rhsbuf
std::vector< TDATA > rhsbuf
Definition ODE.hpp:144
DNDS::ODE::ImplicitSDIRK4DualTimeStep::Fstop
typename tBase::Fstop Fstop
Definition ODE.hpp:140
DNDS::ODE::ImplicitSDIRK4DualTimeStep::Frhs
typename tBase::Frhs Frhs
Definition ODE.hpp:137
DNDS::ODE::ImplicitVBDFDualTimeStep::Fstop
typename tBase::Fstop Fstop
Definition ODE.hpp:609
DNDS::ODE::ImplicitVBDFDualTimeStep::StepPP
void StepPP(TDATA &x, TDATA &xinc, const Frhs &frhs, const Fdt &fdt, const Fsolve &fsolve, int maxIter, const Fstop &fstop, const Fincrement &fincrement, real dt, const FalphaLimSource &falphaLimSource, const FresidualIncPP &fresidualIncPP)
Definition ODE.hpp:811
DNDS::ODE::ImplicitVBDFDualTimeStep::Fsolve
typename tBase::Fsolve Fsolve
Definition ODE.hpp:608
DNDS::ODE::ImplicitVBDFDualTimeStep::ImplicitVBDFDualTimeStep
ImplicitVBDFDualTimeStep(index NDOF, Finit &&finit=[](TDATA &) {}, FinitDtau &&finitDtau=[](TDTAU &) {}, index k=2)
mind that NDOF is the dof of dt finit(TDATA& data)
Definition ODE.hpp:639
DNDS::ODE::ImplicitVBDFDualTimeStep::Step
virtual void Step(TDATA &x, TDATA &xinc, const Frhs &frhs, const Fdt &fdt, const Fsolve &fsolve, int maxIter, const Fstop &fstop, const Fincrement &fincrement, real dt) override
frhs(TDATA &rhs, TDATA &x) fdt(std::vector<real>& dTau) fsolve(TDATA &x, TDATA &rhs,...
Definition ODE.hpp:697
DNDS::ODE::ImplicitVBDFDualTimeStep::FlimitDtBDF
std::function< real(TDATA &, TDATA &)> FlimitDtBDF
Definition ODE.hpp:766
DNDS::ODE::ImplicitVBDFDualTimeStep::getRES
TDATA & getRES(int i) override
Definition ODE.hpp:902
DNDS::ODE::ImplicitVBDFDualTimeStep::getRHS
TDATA & getRHS(int i) override
Definition ODE.hpp:897
DNDS::ODE::ImplicitVBDFDualTimeStep::FresidualIncPP
std::function< void(TDATA &, TDATA &, TDATA &, TDATA &, const std::function< void()> &, real, int)> FresidualIncPP
Definition ODE.hpp:758
DNDS::ODE::ImplicitVBDFDualTimeStep::Fdt
typename tBase::Fdt Fdt
Definition ODE.hpp:607
DNDS::ODE::ImplicitVBDFDualTimeStep::FalphaLimSource
std::function< void(TDATA &, int)> FalphaLimSource
Definition ODE.hpp:762
DNDS::ODE::ImplicitVBDFDualTimeStep::getLatestRHS
virtual TDATA & getLatestRHS() override
Definition ODE.hpp:892
DNDS::ODE::ImplicitVBDFDualTimeStep::rhsbuf
std::vector< TDATA > rhsbuf
Definition ODE.hpp:615
DNDS::ODE::ImplicitVBDFDualTimeStep::LimitDt_StepPPV2
void LimitDt_StepPPV2(TDATA &xIn, const FlimitDtBDF &flimitDtBDF, real &dtOut, real maxIncrease=2)
Definition ODE.hpp:768
DNDS::ODE::ImplicitVBDFDualTimeStep::Frhs
typename tBase::Frhs Frhs
Definition ODE.hpp:606
DNDS::ODE::ImplicitVBDFDualTimeStep::xPrevs
std::vector< TDATA > xPrevs
Definition ODE.hpp:613
DNDS::ODE::ImplicitVBDFDualTimeStep::~ImplicitVBDFDualTimeStep
virtual ~ImplicitVBDFDualTimeStep()=default
DNDS::ODE::ImplicitVBDFDualTimeStep::Fincrement
typename tBase::Fincrement Fincrement
Definition ODE.hpp:610
DNDS::Scalar::BisectSolveLower
real BisectSolveLower(TF &&F, real v0, real v1, real fTarget, int maxIter)
Definition Scalar.hpp:10
DNDS::verySmallReal
DNDS_CONSTANT const real verySmallReal
Catch-all lower bound ("effectively zero").
Definition Defines.hpp:194
DNDS::index
int64_t index
Global row / DOF index type (signed 64-bit; handles multi-billion-cell meshes).
Definition Defines.hpp:107
DNDS::mod
constexpr T mod(T a, T b)
Mathematical modulo that always returns a non-negative result. Unlike % in C++ where (-1) % 3 == -1,...
Definition Defines.hpp:558
DNDS::sqr
constexpr T sqr(const T &a)
a * a, constexpr. Works for all arithmetic types.
Definition Defines.hpp:517
DNDS::real
double real
Canonical floating-point scalar used throughout DNDSR (double precision).
Definition Defines.hpp:105
v
Eigen::Matrix< real, 5, 1 > v
Definition test_ArrayDOF.cpp:468
x
tVec x(NCells)
alpha
real alpha
Definition test_GasThermo.cpp:454