DNDSR 0.1.0.dev1+gcd065ad
Distributed Numeric Data Structure for CFV
Loading...
Searching...
No Matches
EulerSolver.hxx
Go to the documentation of this file.
1/** @file EulerSolver.hxx
2 * @brief Template implementation of EulerSolver::RunImplicitEuler, the main implicit
3 * time-marching loop, along with linear solver dispatch (solveLinear),
4 * preconditioner application (doPrecondition), and running environment
5 * initialization (InitializeRunningEnvironment).
6 *
7 * Covers CFL ramping, reconstruction update, RHS evaluation, LUSGS/GMRES
8 * implicit time stepping, residual monitoring, CL-driver integration,
9 * checkpoint/restart writing, time averaging, and output scheduling.
10 */
11#pragma once
12
13// #ifndef __DNDS_REALLY_COMPILING__
14// #define __DNDS_REALLY_COMPILING__
15// #define __DNDS_REALLY_COMPILING__HEADER_ON__
16// #endif
17#include "EulerSolver.hpp"
18#include "DNDS/EigenUtil.hpp"
19#include "Solver/ODE.hpp"
20#include "Solver/Linear.hpp"
21#include "SpecialFields.hpp"
22// #ifdef __DNDS_REALLY_COMPILING__HEADER_ON__
23// #undef __DNDS_REALLY_COMPILING__
24// #endif
25// #include "fmt/ranges.h"
26
27namespace DNDS::Euler
28{
29 // /***************/ // IDE mode:
30 // static const auto model = NS_SA;
31 // template <>
32 // void EulerSolver<model>::RunImplicitEuler()
33 // /***************/ // IDE mode;
34 static const auto model = NS_SA;
35 DNDS_SWITCH_INTELLISENSE(
36 // the real definition
37 template <EulerModel model>
38 ,
39 // the intellisense friendly definition
40 template <>)
41 /** @brief Main implicit time-marching loop for the compressible Navier-Stokes solver.
42 *
43 * Performs the following at each time step:
44 * 1. CFL ramping with exponential growth and linear fallback on divergence.
45 * 2. Reconstruction update: VR coefficient computation, gradient limiting,
46 * positivity-preserving beta evaluation.
47 * 3. Local time step (dTau) estimation from face eigenvalues.
48 * 4. RHS evaluation via EvaluateRHS.
49 * 5. Implicit time stepping using either LU-SGS sweeps or GMRES with LU preconditioner.
50 * 6. Positivity-preserving alpha limiter for the update.
51 * 7. Component-wise L1 residual monitoring with convergence checks.
52 * 8. CL-driver integration for angle-of-attack adjustment.
53 * 9. Checkpoint/restart writing on schedule or SIGUSR1 signal.
54 * 10. Time averaging accumulation.
55 * 11. VTK output scheduling based on time or iteration.
56 */
57 void EulerSolver<model>::RunImplicitEuler()
58 {
59 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
60 using namespace std::literals;
61
62 auto runningEnvironment = RunningEnvironment();
63 InitializeRunningEnvironment(runningEnvironment);
64 DNDS_EULERSOLVER_RUNNINGENV_GET_REF_LIST
65
66 DNDS_MPI_InsertCheck(mpi, "Implicit 1 nvars " + std::to_string(nVars));
67 /*******************************************************/
68 /* CHECK MESH */
69 /*******************************************************/
70 auto hashCoord = mesh->coords.hash();
71 if (mpi.rank == 0)
72 {
73 log() << "Mesh coord hash is: [" << std::hex << hashCoord << std::dec << "]" << std::endl
74 << std::scientific;
75 }
76 if (config.timeMarchControl.partitionMeshOnly)
77 {
78 if (mpi.rank == 0)
79 {
80 log() << "Mesh Is not altered; partitioning done" << std::endl;
81 }
82 return;
83 }
84
85 /*******************************************************/
86 /* DIRECT FACTORIZATION */
87 /*******************************************************/
88
89 mesh->ObtainLocalFactFillOrdering(*eval.symLU, config.linearSolverControl.directPrecControl);
90 mesh->ObtainSymmetricSymbolicFactorization(*eval.symLU, config.linearSolverControl.directPrecControl);
91 if (config.linearSolverControl.jacobiCode == 2) // do lu on mean-value jacobian
92 {
93 DNDS_MAKE_SSP(JLocalLU, eval.symLU, nVars);
94 }
95
96 // fmt::print("pEval is {}", (void*)(pEval.get()));
97
98 /*******************************************************/
99 /* INITIALIZE U */
100 /*******************************************************/
101
102 eval.InitializeUDOF(u);
103 if (config.timeAverageControl.enabled)
104 wAveraged.setConstant(0.0);
105 if (config.timeMarchControl.useRestart)
106 {
107
108 if (!config.restartState.lastRestartFile.empty())
109 {
110 DNDS_assert(config.restartState.iStep >= 1);
111 ReadRestart(config.restartState.lastRestartFile);
112 }
113 if (!config.restartState.otherRestartFile.empty())
114 ReadRestartOtherSolver(config.restartState.otherRestartFile, config.restartState.otherRestartStoreDim);
115 }
116 OutputPicker outputPicker;
117 eval.InitializeOutputPicker(outputPicker, {u, uRec, betaPP, alphaPP});
118 addOutList = outputPicker.getSubsetList(config.dataIOControl.outCellScalarNames);
119
120 /*******************************************************/
121 /* TEMPORARY Us */
122 /*******************************************************/
123
124 auto initUDOF = [&](ArrayDOFV<nVarsFixed> &uu)
125 { vfv->BuildUDof(uu, nVars); };
126 auto initUREC = [&](ArrayRECV<nVarsFixed> &uu)
127 { vfv->BuildURec(uu, nVars); };
128
129#define DNDS_EULER_SOLVER_GET_TEMP_UDOF(name) \
130 auto __p##name = uPool.getAllocInit(initUDOF); \
131 auto &name = *__p##name;
132
133 /*******************************************************/
134 /* SOLVER MAJOR START */
135 /*******************************************************/
136 tstart = MPI_Wtime();
137 tstartInternal = tstart;
138
139 Timer().clearAllTimer();
140
141 DNDS_MPI_InsertCheck(mpi, "Implicit 2 nvars " + std::to_string(nVars));
142
143 /*******************************************************/
144 /* DEFINE LAMBDAS */
145 /*******************************************************/
146
147 auto frhsOuter =
148 [&](
149 ArrayDOFV<nVarsFixed> &crhs,
150 ArrayDOFV<nVarsFixed> &cx,
151 ArrayDOFV<1> &dTau,
152 int iter, real ct, int uPos, int reconstructionFlag)
153 {
154 cx.trans.startPersistentPull();
155 cx.trans.waitPersistentPull(); // for hermite3
156 auto &uRecC = config.timeMarchControl.timeMarchIsTwoStage() && uPos == 1 ? uRec1 : uRec;
157 auto &JSourceC = config.timeMarchControl.timeMarchIsTwoStage() && uPos == 1 ? JSource1 : JSource;
158 auto &uRecIncC = config.timeMarchControl.timeMarchIsTwoStage() && uPos == 1 ? uRecInc1 : uRecInc;
159 auto &alphaPPC = config.timeMarchControl.timeMarchIsTwoStage() && uPos == 1 ? alphaPP1 : alphaPP;
160 auto &betaPPC = config.timeMarchControl.timeMarchIsTwoStage() && uPos == 1 ? betaPP1 : betaPP;
161 // if (mpi.rank == 0)
162 // std::cout << uRecC.father.get() << std::endl;
163 typename TVFV::template TFBoundary<nVarsFixed>
164 FBoundary = [&](const TU &UL, const TU &UMean, index iCell, index iFace, int ig,
165 const Geom::tPoint &normOut, const Geom::tPoint &pPhy, const Geom::t_index bType) -> TU
166 {
167 TVec normOutV = normOut(Seq012);
168 Eigen::Matrix<real, dim, dim> normBase = Geom::NormBuildLocalBaseV<dim>(normOutV);
169 bool compressed = false;
170 TU ULfixed = eval.CompressRecPart(
171 UMean,
172 UL - UMean,
173 compressed);
174 return eval.generateBoundaryValue(ULfixed, UMean, iCell, iFace, ig, normOutV, normBase, pPhy, tSimu + ct * curDtImplicit, bType, true, 1);
175 };
176 typename TVFV::template TFBoundaryDiff<nVarsFixed>
177 FBoundaryDiff = [&](const TU &UL, const TU &dU, const TU &UMean, index iCell, index iFace, int ig,
178 const Geom::tPoint &normOut, const Geom::tPoint &pPhy, const Geom::t_index bType) -> TU
179 {
180 TVec normOutV = normOut(Seq012);
181 Eigen::Matrix<real, dim, dim> normBase = Geom::NormBuildLocalBaseV<dim>(normOutV);
182 bool compressed = false;
183 TU ULfixed = eval.CompressRecPart(
184 UMean,
185 UL - UMean,
186 compressed);
187 TU ULfixedPlus = eval.CompressRecPart(
188 UMean,
189 UL - UMean + dU,
190 compressed);
191 return eval.generateBoundaryValue(ULfixedPlus, UMean, iCell, iFace, ig, normOutV, normBase, pPhy, tSimu + ct * curDtImplicit, bType, true, 1) -
192 eval.generateBoundaryValue(ULfixed, UMean, iCell, iFace, ig, normOutV, normBase, pPhy, tSimu + ct * curDtImplicit, bType, true, 1);
193 };
194
195 eval.FixUMaxFilter(cx);
196 if ((iter > config.convergenceControl.nAnchorUpdateStart &&
197 (iter - config.convergenceControl.nAnchorUpdateStart - 1) % config.convergenceControl.nAnchorUpdate == 0))
198 {
199 eval.updateBCAnchors(cx, uRecC);
200 eval.updateBCProfiles(cx, uRecC);
201 eval.updateBCProfilesPressureRadialEq();
202 }
203 // cx.trans.startPersistentPull();
204 // cx.trans.waitPersistentPull();
205
206 // for (index iCell = 0; iCell < uOld.size(); iCell++)
207 // uOld[iCell].m() = uRec[iCell].m();
208 if (!reconstructionFlag)
209 {
210 betaPPC.setConstant(1.0);
211 alphaPP_tmp.setConstant(1.0);
212 uRecNew.setConstant(0.0);
213 eval.EvaluateRHS(crhs, JSourceC, cx, uRecNew, uRecNew, betaPPC, alphaPP_tmp, false, tSimu + ct * curDtImplicit,
214 TEval::RHS_Ignore_Viscosity); // TODO: test with viscosity // TODO: RHS_Direct_2nd_Rec_1st_Conv?
215 // vfv->DoReconstruction2nd(uRecOld, cx, FBoundary, 1, std::vector<int>());
216 // eval.EvaluateRHS(crhs, JSourceC, cx, uRecOld, uRecNew, betaPPC, alphaPP_tmp, false, tSimu + ct * curDtImplicit,
217 // 0); // TEval::RHS_Ignore_Viscosity
218 return;
219 }
220
221 DNDS_MPI_InsertCheck(mpi, " Lambda RHS: StartRec");
222 int nRec = (gradIsZero ? config.implicitReconstructionControl.nRecMultiplyForZeroedGrad : 1) *
223 config.implicitReconstructionControl.nInternalRecStep;
224 if (step <= config.implicitReconstructionControl.zeroRecForSteps ||
225 iter <= config.implicitReconstructionControl.zeroRecForStepsInternal)
226 {
227 nRec = 0;
228 uRec.setConstant(0.0);
229 }
230 real recIncBase = 0;
231 Timer().StartTimer(PerformanceTimer::Reconstruction);
232 if (config.implicitReconstructionControl.storeRecInc)
233 uRecOld = uRecC;
234
235 if (config.implicitReconstructionControl.useExplicit)
236 { // pass
237 }
238 else if (config.implicitReconstructionControl.recLinearScheme == 0)
239 {
240 for (int iRec = 1; iRec <= nRec; iRec++)
241 {
242 if (nRec > 1)
243 uRecNew1 = uRecC;
244
245 vfv->DoReconstructionIter(
246 uRecC, uRecNew, cx,
247 FBoundary,
248 false);
249
250 uRecC.trans.startPersistentPull();
251 uRecC.trans.waitPersistentPull();
252
253 if (nRec > 1)
254 {
255 uRecNew1 -= uRecC;
256 real recInc = uRecNew1.norm2();
257 if (iRec == 1)
258 recIncBase = recInc;
259
260 if (iRec % config.implicitReconstructionControl.nRecConsolCheck == 0)
261 {
262 if (mpi.rank == 0)
263 log() << iRec << " Rec inc: " << recIncBase << " -> " << recInc << std::endl;
264 }
265 if (recInc < recIncBase * config.implicitReconstructionControl.recThreshold)
266 break;
267 }
268 }
269 }
270 else if (config.implicitReconstructionControl.recLinearScheme == 1)
271 {
272 Eigen::VectorXd meanScale;
273 if (config.implicitReconstructionControl.gmresRecScale == 1)
274 {
275 meanScale = eval.settings.refU;
276 meanScale(Seq123).setConstant(std::sqrt(meanScale(0) * meanScale(I4))); //! using consistent rho U scale
277 }
278 else
279 meanScale.setConstant(nVars, 1.0);
280 // meanScale(0) = 10;
281 TU meanScaleInv = (meanScale.array() + verySmallReal).inverse();
282
283 int nGMRESrestartAll{0};
284 real gmresResidualB = 0;
285 for (int iRec = 1; iRec <= nRec; iRec++)
286 {
287 vfv->DoReconstructionIter(
288 uRecC, uRecNew, cx,
289 FBoundary,
290 true, true);
291 uRecNew.trans.startPersistentPull();
292 uRecNew.trans.waitPersistentPull();
293 // uRec.setConstant(0.0);
294 uRecNew1 = uRecNew;
295 if (iRec == 1)
296 gmresResidualB = uRecNew1.norm2();
297
298 bool gmresConverge =
299 gmresRec->solve(
300 [&](ArrayRECV<nVarsFixed> &x, ArrayRECV<nVarsFixed> &Ax)
301 {
302 vfv->DoReconstructionIterDiff(uRec, x, Ax, cx, FBoundaryDiff);
303 Ax.trans.startPersistentPull();
304 Ax.trans.waitPersistentPull();
305 },
306 [&](ArrayRECV<nVarsFixed> &x, ArrayRECV<nVarsFixed> &MLx)
307 {
308 MLx = x; // initial value; for the input is mostly a good estimation
309 // MLx no need to comm
310 // vfv->DoReconstructionIterSOR(uRecC, x, MLx, cx, FBoundaryDiff, false); //! causes loss of accuracy; why?????
311 },
312 [&](ArrayRECV<nVarsFixed> &a, ArrayRECV<nVarsFixed> &b) -> real
313 {
314 return (a.dotV(b).array() * meanScaleInv.transpose().array().square()).sum(); //! dim balancing here
315 },
316 uRecNew, uRecNew1, config.implicitReconstructionControl.nGmresIter,
317 [&](uint32_t i, real res, real resB) -> bool
318 {
319 if (i > 0)
320 {
321 nGMRESrestartAll++;
322 if (mpi.rank == 0 &&
323 (nGMRESrestartAll % config.implicitReconstructionControl.nRecConsolCheck == 0))
324 {
325 auto fmt = log().flags();
326 log() << std::scientific;
327 log() << "GMRES for Rec: " << iRec << " Restart " << i << " " << resB << " -> " << res << std::endl;
328 log().setf(fmt);
329 }
330 }
331 return res < gmresResidualB * config.implicitReconstructionControl.recThreshold;
332 });
333 uRecC.addTo(uRecNew1, -1);
334 if (gmresConverge)
335 break;
336 }
337 uRecC.trans.startPersistentPull();
338 uRecC.trans.waitPersistentPull();
339 }
340 else if (config.implicitReconstructionControl.recLinearScheme == 2)
341 {
342 Eigen::Array<real, 1, Eigen::Dynamic> resB;
343 int nPCGIterAll{0};
344 auto &pcgRecC = config.timeMarchControl.timeMarchIsTwoStage() && uPos == 1 ? pcgRec1 : pcgRec;
345 auto &uRecBC = config.timeMarchControl.timeMarchIsTwoStage() && uPos == 1 ? uRecB1 : uRecB;
346
347 if (iter <= 2 //! consecutive pcg is bad in 0012, using separate pcg
348 || pcgRecC->getPHistorySize() >= config.implicitReconstructionControl.fpcgMaxPHistory)
349 pcgRecC->reset();
350 uRecNew = uRecBC; // last version of B
351 vfv->DoReconstructionIter(
352 uRecC, uRecBC, cx,
353 FBoundary, true, true, true); // puts rhs of reconstruction A^-1b into uRecBC
354 if (iter > 2)
355 {
356 if (config.implicitReconstructionControl.fpcgResetScheme == 0)
357 {
358 Eigen::RowVectorXd bPrevSqr = uRecNew.dotV(uRecNew);
359 uRecNew.addTo(uRecBC, -1.0);
360 Eigen::RowVectorXd bIncSqr = uRecNew.dotV(uRecNew);
361 real maxPortion = (bIncSqr.array() / (bPrevSqr.array() + smallReal)).sqrt().maxCoeff();
362 if (maxPortion >= config.implicitReconstructionControl.fpcgResetThres)
363 {
364 if (mpi.rank == 0 && config.implicitReconstructionControl.fpcgResetReport > 0)
365 log() << "FPCG force reset at portion " << fmt::format("{:.4g}", maxPortion) << std::endl;
366 pcgRecC->reset();
367 }
368 }
369 else if (config.implicitReconstructionControl.fpcgResetScheme == 1)
370 pcgRecC->reset();
371 else if (config.implicitReconstructionControl.fpcgResetScheme == 2)
372 ;
373 else
374 DNDS_assert_info(false, "invalid fpcgResetScheme");
375 }
376 // pcgRecC->reset(); // using separate pcg
377 for (int iRec = 1; iRec <= nRec; iRec++)
378 {
379 bool pcgConverged = pcgRecC->solve(
380 [&](ArrayRECV<nVarsFixed> &x, ArrayRECV<nVarsFixed> &Ax)
381 {
382 vfv->DoReconstructionIterDiff(uRec, x, Ax, cx, FBoundaryDiff);
383 Ax.trans.startPersistentPull();
384 Ax.trans.waitPersistentPull();
385 vfv->MatrixAMult(Ax, Ax);
386 },
387 [&](ArrayRECV<nVarsFixed> &x, ArrayRECV<nVarsFixed> &Mx)
388 {
389 vfv->MatrixAMult(x, Mx);
390 },
391 [&](ArrayRECV<nVarsFixed> &x, ArrayRECV<nVarsFixed> &res)
392 {
393 vfv->DoReconstructionIter(
394 x, res, cx,
395 FBoundary, true, true);
396 res.trans.startPersistentPull();
397 res.trans.waitPersistentPull();
398 res *= -1;
399 },
400 [&](ArrayRECV<nVarsFixed> &a, ArrayRECV<nVarsFixed> &b)
401 {
402 return (a.dotV(b)).array();
403 },
404 uRecC, config.implicitReconstructionControl.nGmresIter,
405 [&](uint32_t i, const Eigen::Array<real, 1, Eigen::Dynamic> &res, const Eigen::Array<real, 1, Eigen::Dynamic> &res1) -> bool
406 {
407 if (i == 1 && iRec == 1)
408 resB = res1;
409 if (i > 0)
410 {
411 nPCGIterAll++;
412 if (mpi.rank == 0 &&
413 (nPCGIterAll % config.implicitReconstructionControl.nRecConsolCheck == 0))
414 {
415 auto fmt = log().flags();
416 log() << std::scientific << std::setprecision(3);
417 log() << "PCG for Rec: " << iRec << " Restart " << i << " [" << resB << "] -> [" << res1 << "]" << std::endl;
418 log().setf(fmt);
419 }
420 }
421 return (res1 / resB).maxCoeff() < config.implicitReconstructionControl.recThreshold;
422 });
423 if (pcgConverged)
424 break;
425 }
426 uRecC.trans.startPersistentPull();
427 uRecC.trans.waitPersistentPull();
428 }
429 else
430 DNDS_assert_info(false, "no such recLinearScheme");
431 if ((model == NS_SA || model == NS_SA_3D) && eval.settings.ransForce2nd)
432 {
433 std::vector<int> mask;
434 mask.resize(1);
435 mask[0] = 5;
436 vfv->DoReconstruction2nd(uRec, u, FBoundary, 1, mask);
437 }
438 if (model == NS_2EQ || model == NS_2EQ_3D)
439 {
440 std::vector<int> mask;
441 mask.resize(2);
442 mask[0] = 5;
443 mask[1] = 6;
444 vfv->DoReconstruction2nd(uRec, u, FBoundary, 1, mask);
445 }
446 Timer().StopTimer(PerformanceTimer::Reconstruction);
447 if (gradIsZero)
448 {
449 uRec = uRecC;
450 if (config.timeMarchControl.timeMarchIsTwoStage())
451 uRec1 = uRecC;
452 gradIsZero = false;
453 }
454
455 Timer().StartTimer(PerformanceTimer::Limiter);
456 // for (index iCell = 0; iCell < uOld.size(); iCell++)
457 // uRecC[iCell].m() -= uOld[iCell].m();
458
459 DNDS_MPI_InsertCheck(mpi, " Lambda RHS: StartLim");
460 if (!config.implicitReconstructionControl.useExplicit && config.limiterControl.useLimiter)
461 {
462 // vfv->ReconstructionWBAPLimitFacial(
463 // cx, uRecC, uRecNew, uF0, uF1, ifUseLimiter,
464
465 using tLimitBatch = typename TVFV::template tLimitBatch<nVarsFixed>;
466
467 auto fML = [&](const TU &UL, const TU &UR, const Geom::tPoint &n,
468 const Eigen::Ref<tLimitBatch> &data) -> tLimitBatch
469 {
470 Timer().StartTimer(PerformanceTimer::LimiterA);
471 Eigen::Vector<real, I4 + 1> UC = (UL + UR)(Seq01234) * 0.5;
472 Eigen::Matrix<real, dim, dim> normBase = Geom::NormBuildLocalBaseV<dim>(n(Seq012));
473 UC(Seq123) = normBase.transpose() * UC(Seq123);
474
475 auto M = Gas::IdealGas_EulerGasLeftEigenVector<dim>(UC, eval.settings.idealGasProperty.gamma);
476 M(EigenAll, Seq123) *= normBase.transpose();
477
478 Eigen::Matrix<real, nVarsFixed, nVarsFixed> ret(nVars, nVars);
479 ret.setIdentity();
480 ret(Seq01234, Seq01234) = M;
481 Timer().StopTimer(PerformanceTimer::LimiterA);
482 return (ret * data.transpose()).transpose();
483 // return real(1);
484 };
485 auto fMR = [&](const TU &UL, const TU &UR, const Geom::tPoint &n,
486 const Eigen::Ref<tLimitBatch> &data) -> tLimitBatch
487 {
488 Timer().StartTimer(PerformanceTimer::LimiterA);
489 Eigen::Vector<real, I4 + 1> UC = (UL + UR)(Seq01234) * 0.5;
490 Eigen::Matrix<real, dim, dim> normBase = Geom::NormBuildLocalBaseV<dim>(n(Seq012));
491 UC(Seq123) = normBase.transpose() * UC(Seq123);
492
493 auto M = Gas::IdealGas_EulerGasRightEigenVector<dim>(UC, eval.settings.idealGasProperty.gamma);
494 M(Seq123, EigenAll) = normBase * M(Seq123, EigenAll);
495
496 Eigen::Matrix<real, nVarsFixed, nVarsFixed> ret(nVars, nVars);
497 ret.setIdentity();
498 ret(Seq01234, Seq01234) = M;
499
500 Timer().StopTimer(PerformanceTimer::LimiterA);
501 return (ret * data.transpose()).transpose();
502 // return real(1);
503 };
504 if (config.limiterControl.smoothIndicatorProcedure == 0)
505 vfv->template DoCalculateSmoothIndicator<nVarsFixed, 2>(
506 ifUseLimiter, (uRecC), (u),
507 std::array<int, 2>{0, I4});
508 else if (config.limiterControl.smoothIndicatorProcedure == 1)
509 vfv->template DoCalculateSmoothIndicatorV1<nVarsFixed>(
510 ifUseLimiter, (uRecC), (u),
511 TU::Ones(nVars),
512 [&](Eigen::Matrix<real, 1, nVarsFixed> &v)
513 {
514 TU prim;
515 TU cons;
516 cons = v.transpose();
517 if (cons(0) < verySmallReal)
518 {
519 v.setConstant(-veryLargeReal * v(I4));
520 return;
521 }
522 Gas::IdealGasThermalConservative2Primitive<dim>(cons, prim, eval.settings.idealGasProperty.gamma);
523 v.setConstant(prim(I4));
524 return;
525 });
526
527 else
528 {
529 DNDS_assert(false);
530 }
531 if (config.limiterControl.limiterProcedure == 1)
532 vfv->template DoLimiterWBAP_C<nVarsFixed>(
533 (cx),
534 (uRecC),
535 (uRecLimited),
536 (uRecNew),
537 ifUseLimiter,
538 iter < config.limiterControl.nPartialLimiterStartLocal && step < config.limiterControl.nPartialLimiterStart,
539 fML, fMR, true);
540 else if (config.limiterControl.limiterProcedure == 0)
541 vfv->template DoLimiterWBAP_3<nVarsFixed>(
542 (cx),
543 (uRecC),
544 (uRecLimited),
545 (uRecNew),
546 ifUseLimiter,
547 iter < config.limiterControl.nPartialLimiterStartLocal && step < config.limiterControl.nPartialLimiterStart,
548 fML, fMR, true);
549 else
550 {
551 DNDS_assert(false);
552 }
553 // uRecLimited.trans.startPersistentPull();
554 // uRecLimited.trans.waitPersistentPull();
555 }
556 Timer().StopTimer(PerformanceTimer::Limiter);
557 if (config.implicitReconstructionControl.storeRecInc)
558 {
559 uRecIncC = uRecC;
560 uRecIncC -= uRecOld; //! uRecIncC now stores uRecIncrement
561 }
562 if (config.implicitReconstructionControl.storeRecInc && config.implicitReconstructionControl.dampRecIncDTau)
563 {
564 ArrayDOFV<1> damper;
565 ArrayDOFV<1> damper1;
566 vfv->BuildUDof(damper, 1);
567 vfv->BuildUDof(damper1, 1);
568 damper = dTau;
569 damper1 = dTau;
570 damper1 += curDtImplicit;
571 damper /= damper1;
572 // for (auto &v : damper)
573 // v = v / (curDtImplicit + v); //! warning: teleported value here
574
575 uRecIncC *= damper;
576 uRecC = uRecOld;
577 uRecC += uRecIncC;
578 }
579 if (config.limiterControl.preserveLimited && config.limiterControl.useLimiter)
580 uRecC = uRecLimited;
581
582 // uRecC.trans.startPersistentPull(); //! this also need to update!
583 // uRecC.trans.waitPersistentPull();
584
585 // }
586
587 DNDS_MPI_InsertCheck(mpi, " Lambda RHS: StartEval");
588 eval.setPassiveDiscardSource(iter <= 0);
589
590 if (iter == 1)
591 alphaPPC.setConstant(1.0); // make RHS un-disturbed
592 alphaPP_tmp.setConstant(1.0); // make RHS un-disturbed
593 if (!config.implicitReconstructionControl.useExplicit && config.limiterControl.usePPRecLimiter)
594 {
595 Timer().StartTimer(PerformanceTimer::Positivity);
596 nLimBeta = 0;
597 minBeta = 1;
598 if (!config.limiterControl.useLimiter)
599 uRecLimited = uRecC;
600 eval.EvaluateURecBeta(cx, uRecLimited, betaPPC, nLimBeta, minBeta,
601 config.limiterControl.ppRecLimiterCompressToMean
602 ? TEval::EvaluateURecBeta_COMPRESS_TO_MEAN
603 : TEval::EvaluateURecBeta_DEFAULT); //*cx instead of u!
604 if (nLimBeta)
605 if (mpi.rank == 0 &&
606 (config.outputControl.consoleOutputEveryFix == 1 || config.outputControl.consoleOutputEveryFix == 3))
607 {
608 log() << std::scientific << std::setprecision(config.outputControl.nPrecisionConsole)
609 << "PPRecLimiter: nLimBeta [" << nLimBeta << "]"
610 << " minBeta[" << minBeta << "]" << std::endl;
611 }
612 uRecLimited.trans.startPersistentPull();
613 betaPPC.trans.startPersistentPull();
614 uRecLimited.trans.waitPersistentPull();
615 betaPPC.trans.waitPersistentPull();
616 Timer().StopTimer(PerformanceTimer::Positivity);
617 }
618
619 Timer().StartTimer(PerformanceTimer::RHS);
620 if (config.implicitReconstructionControl.useExplicit)
621 eval.EvaluateRHS(crhs, JSourceC, cx, uRecC /* dummy*/, uRecC /* dummy*/,
622 betaPPC /* dummy*/, alphaPP_tmp /* dummy*/, false, tSimu + ct * curDtImplicit,
623 TEval::RHS_Direct_2nd_Rec | (config.limiterControl.useLimiter ? TEval::RHS_Direct_2nd_Rec_use_limiter : TEval::RHS_No_Flags));
624 else if (config.limiterControl.useLimiter || config.limiterControl.usePPRecLimiter) // todo: opt to using limited for uRecUnlim
625 eval.EvaluateRHS(crhs, JSourceC, cx, config.limiterControl.useViscousLimited ? uRecLimited : uRecC, uRecLimited,
626 betaPPC, alphaPP_tmp, false, tSimu + ct * curDtImplicit);
627 else
628 eval.EvaluateRHS(crhs, JSourceC, cx, uRecC, uRecC,
629 betaPPC, alphaPP_tmp, false, tSimu + ct * curDtImplicit);
630
631 crhs.trans.startPersistentPull();
632 crhs.trans.waitPersistentPull();
633 if (getNVars(model) > (I4 + 1) && iter <= config.others.nFreezePassiveInner)
634 {
635 for (int i = 0; i < crhs.Size(); i++)
636 crhs[i](Eigen::seq(I4 + 1, EigenLast)).setZero();
637 // if (mpi.rank == 0)
638 // std::cout << "Freezing all passive" << std::endl;
639 }
640 Timer().StopTimer(PerformanceTimer::RHS);
641
642 DNDS_MPI_InsertCheck(mpi, " Lambda RHS: End");
643 };
644
645 auto frhs =
646 [&](
647 ArrayDOFV<nVarsFixed> &crhs,
648 ArrayDOFV<nVarsFixed> &cx,
649 ArrayDOFV<1> &dTau,
650 int iter, real ct, int uPos)
651 {
652 if (config.timeMarchControl.timeMarchIsTwoStage() && uPos == 2) // special for two stage (HM3):
653 {
654 // ! mind that JSource, betaPP uses as uPos == 0
655 return eval.EvaluateRHS(crhs, JSource, cx, uRecNew, uRecNew, betaPP, alphaPP_tmp, false, tSimu + ct * curDtImplicit,
656 TEval::RHS_Direct_2nd_Rec | TEval::RHS_Dont_Record_Bud_Flux | TEval::RHS_Dont_Update_Integration |
657 TEval::RHS_Direct_2nd_Rec_already_have_uGradBufNoLim | //! uGradBufNoLim already existent in fdtau
658 (config.limiterControl.useLimiter ? TEval::RHS_Direct_2nd_Rec_use_limiter : TEval::RHS_No_Flags));
659 //! note: in HM3 IV test for TPMG, this O1 version makes convergence slower compared to the O2 one, why?
660 // static const int use_1st_conv = 1;
661 // static const int use_1st_conv_ignore_vis = 0;
662 // return eval.EvaluateRHS(crhs, JSource, cx, uRecNew, uRecNew, betaPP, alphaPP_tmp, false, tSimu + ct * curDtImplicit,
663 // (TEval::RHS_Direct_2nd_Rec_1st_Conv * use_1st_conv) |
664 // TEval::RHS_Direct_2nd_Rec |
665 // TEval::RHS_Dont_Record_Bud_Flux |
666 // TEval::RHS_Dont_Update_Integration |
667 // (TEval::RHS_Ignore_Viscosity * use_1st_conv_ignore_vis) |
668 // TEval::RHS_Direct_2nd_Rec_already_have_uGradBufNoLim | //! uGradBufNoLim already existent in fdtau
669 // (config.limiterControl.useLimiter ? TEval::RHS_Direct_2nd_Rec_use_limiter : TEval::RHS_No_Flags));
670 }
671 return frhsOuter(crhs, cx, dTau, iter, ct, uPos, 1); // reconstructionFlag == 1
672 };
673
674 auto fdtau = [&](ArrayDOFV<nVarsFixed> &cx, ArrayDOFV<1> &dTau, real alphaDiag, int uPos)
675 {
676 eval.FixUMaxFilter(cx);
677 cx.trans.startPersistentPull(); //! this also need to update!
678 cx.trans.waitPersistentPull();
679 // uRec.trans.startPersistentPull();
680 // uRec.trans.waitPersistentPull();
681 auto &uRecC = config.timeMarchControl.timeMarchIsTwoStage() && uPos == 1 ? uRec1 : uRec;
682 eval.EvaluateDt(dTau, cx, uRecC, CFLNow, curDtMin, 1e100, config.implicitCFLControl.useLocalDt, tSimu);
683 for (int iS = 1; iS <= config.implicitCFLControl.nSmoothDTau; iS++)
684 {
685 // ArrayDOFV<1> dTauNew = dTau; //TODO: copying is still unusable; consider doing copiers on the level of ArrayDOFV and ArrayRecV
686 dTauTmp = dTau;
687 dTau.trans.startPersistentPull();
688 dTau.trans.waitPersistentPull();
689 eval.MinSmoothDTau(dTau, dTauTmp);
690 dTau = dTauTmp;
691 }
692
693 dTau *= 1. / alphaDiag;
694 };
695
696 auto fincrement = [&](
697 ArrayDOFV<nVarsFixed> &cx,
698 ArrayDOFV<nVarsFixed> &cxInc,
699 real alpha, int uPos)
700 {
701 Timer().StartTimer(PerformanceTimer::Positivity);
702 auto &alphaPPC = config.timeMarchControl.timeMarchIsTwoStage() && uPos == 1 ? alphaPP1 : alphaPP;
703 auto &betaPPC = config.timeMarchControl.timeMarchIsTwoStage() && uPos == 1 ? betaPP1 : betaPP;
704 auto &uRecC = config.timeMarchControl.timeMarchIsTwoStage() && uPos == 1 ? uRec1 : uRec;
705 auto &JSourceC = config.timeMarchControl.timeMarchIsTwoStage() && uPos == 1 ? JSource1 : JSource;
706 nLimInc = 0;
707 alphaMinInc = 1;
708 eval.EvaluateCellRHSAlpha(cx, uRecC, betaPPC, cxInc, alphaPP_tmp, nLimInc, alphaMinInc, config.timeMarchControl.incrementPPRelax,
709 2, TEval::EvaluateCellRHSAlpha_DEFAULT);
710 if (nLimInc)
711 if (mpi.rank == 0 &&
712 (config.outputControl.consoleOutputEveryFix == 1 || config.outputControl.consoleOutputEveryFix == 2))
713 {
714 log() << std::scientific << std::setw(3) << TermColor::Red;
715 log() << "PPIncrementLimiter: nIncrementRes[" << nLimInc << "] minAlpha [" << alphaMinInc << "]"
716 << TermColor::Reset
717 << std::endl;
718 }
719 auto pUTemp = uPool.getAllocInit(initUDOF);
720 auto &uTemp = *pUTemp;
721
722 uTemp = cxInc;
723 uTemp *= alphaPP_tmp;
724 //*directadd
725 // cx += uTemp;
726 //*fixing add
727 if (model == NS_2EQ || model == NS_2EQ_3D)
728 {
729 if (config.implicitCFLControl.RANSRelax != 1)
730 for (index i = 0; i < uTemp.Size(); i++)
731 uTemp[i]({I4, I4 + 1}) *= config.implicitCFLControl.RANSRelax;
732 }
733 Timer().StopTimer(PerformanceTimer::Positivity);
734 eval.AddFixedIncrement(cx, uTemp, alpha);
735 eval.AssertMeanValuePP(cx, true);
736
737 // eval.AddFixedIncrement(cx, cxInc, alpha);
738 };
739
740 auto fsolve = [&](ArrayDOFV<nVarsFixed> &cx, ArrayDOFV<nVarsFixed> &cres, ArrayDOFV<nVarsFixed> &resOther, ArrayDOFV<1> &dTau,
741 real dt, real alphaDiag, ArrayDOFV<nVarsFixed> &cxInc, int iter, real ct, int uPos)
742 {
743 {
744 DNDS_EULER_SOLVER_GET_TEMP_UDOF(rhsTemp)
745 DNDS_EULER_SOLVER_GET_TEMP_UDOF(uTemp)
746 rhsTemp = cres;
747 eval.CentralSmoothResidual(rhsTemp, cres, uTemp);
748 }
749
750 auto &JDC = config.timeMarchControl.timeMarchIsTwoStage() && uPos == 1 ? JD1 : JD;
751 auto &JSourceC = config.timeMarchControl.timeMarchIsTwoStage() && uPos == 1 ? JSource1 : JSource;
752 auto &uRecC = config.timeMarchControl.timeMarchIsTwoStage() && uPos == 1 ? uRec1 : uRec;
753 auto &uRecIncC = config.timeMarchControl.timeMarchIsTwoStage() && uPos == 1 ? uRecInc1 : uRecInc;
754 auto &alphaPPC = config.timeMarchControl.timeMarchIsTwoStage() && uPos == 1 ? alphaPP1 : alphaPP;
755 auto &betaPPC = config.timeMarchControl.timeMarchIsTwoStage() && uPos == 1 ? betaPP1 : betaPP;
756 bool isTPMGLevel = config.timeMarchControl.timeMarchIsTwoStage() && uPos == 2;
757
758 Timer().StartTimer(PerformanceTimer::Positivity);
759 if (config.timeMarchControl.rhsFPPMode == 1 || config.timeMarchControl.rhsFPPMode == 11)
760 {
761 DNDS_EULER_SOLVER_GET_TEMP_UDOF(rhsTemp)
762 // ! experimental: bad now ?
763 rhsTemp = cres;
764 rhsTemp *= dTau;
765 rhsTemp *= config.timeMarchControl.rhsFPPScale;
766 index nLimFRes = 0;
767 real alphaMinFRes = 1;
768 eval.EvaluateCellRHSAlpha(cx, uRecC, betaPPC, rhsTemp, alphaPP_tmp, nLimFRes, alphaMinFRes, config.timeMarchControl.rhsFPPRelax,
769 2, config.timeMarchControl.rhsFPPMode == 1 ? TEval::EvaluateCellRHSAlpha_DEFAULT : TEval::EvaluateCellRHSAlpha_MIN_IF_NOT_ONE);
770 if (nLimFRes)
771 if (mpi.rank == 0)
772 {
773 log() << std::scientific << std::setprecision(config.outputControl.nPrecisionConsole) << std::setw(3) << TermColor::Red;
774 log() << "PPFResLimiter: nLimFRes[" << nLimFRes << "] minAlpha [" << alphaMinFRes << "]" << TermColor::Reset << std::endl;
775 }
776
777 cres *= alphaPP_tmp;
778 }
779 else if (config.timeMarchControl.rhsFPPMode == 2)
780 {
781 DNDS_EULER_SOLVER_GET_TEMP_UDOF(rhsTemp)
782 rhsTemp = cres;
783 rhsTemp *= dTau;
784 rhsTemp *= config.timeMarchControl.rhsFPPScale;
785 index nLimFRes = 0;
786 real alphaMinFRes = 1;
787 eval.EvaluateCellRHSAlpha(cx, uRecC, betaPPC, rhsTemp, alphaPP_tmp, nLimFRes, alphaMinFRes, config.timeMarchControl.rhsFPPRelax,
788 2, TEval::EvaluateCellRHSAlpha_DEFAULT);
789 if (nLimFRes)
790 if (mpi.rank == 0)
791 {
792 log() << std::scientific << std::setw(3) << TermColor::Red;
793 log() << "PPFResLimiter: nLimFRes[" << nLimFRes << "] minAlpha [" << alphaMinFRes << "]" << TermColor::Reset << std::endl;
794 }
795 alphaPP_tmp.setMaxWith(smallReal); // dTau cannot be zero
796 dTauTmp = dTau;
797 dTauTmp *= alphaPP_tmp;
798 }
799 auto &dTauC = config.timeMarchControl.rhsFPPMode == 2 ? dTauTmp : dTau;
800 Timer().StopTimer(PerformanceTimer::Positivity);
801
802 if (config.limiterControl.useLimiter) // uses urec value
803 eval.LUSGSMatrixInit(JDC, JSourceC,
804 dTauC, dt, alphaDiag,
805 cx, uRecLimited,
806 0,
807 tSimu);
808 else
809 eval.LUSGSMatrixInit(JDC, JSourceC,
810 dTauC, dt, alphaDiag,
811 cx, uRecC,
812 0,
813 tSimu);
814
815 // for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
816 // {
817 // cres[iCell] = eval.CompressInc(cx[iCell], cres[iCell] * dTau[iCell]) / dTau[iCell];
818 // }
819 cxInc.setConstant(0.0);
820 this->solveLinear(alphaDiag, tSimu, cres, cx, cxInc, uRecC, uRecIncC,
821 JDC, *gmres, !isTPMGLevel ? 0 : 1); //! here we borrow PMG's level1 setting into TPMG
822 // cxInc: in: full increment from previous level; out: full increment form current level
823 const auto solve_multigrid = [&](TDof &x_upper, TDof &xIncBuf, TDof &rhsBuf, const TDof &resOther, int mgLevelInit, int mgLevelMax)
824 {
825 uRecNew.setConstant(0.0);
826 // templated lambda recursion:
827 // https://stackoverflow.com/questions/2067988/how-to-make-a-recursive-lambda
828 // std::function<void(TDof &, TDof &, int, int)> solve_multigrid_impl;
829 // solve_multigrid_impl = [&](TDof &x_base, TDof &cxInc, int mgLevel, int mgLevelMax)
830
831 // [&frhs, &fincrement, &fdtau, &initUDOF,
832 // &eval, &config = config, &dTauC, &cres, &resOther,
833 // &mpi = mpi, &JSourceTmp = JSourceTmp, &JDTmp = JDTmp,
834 // &uPool = uPool, &uRecNew = uRecNew, &uRec = uRec, &betaPP = betaPP, &alphaPP = alphaPP,
835 // &gmres, tSimu, solveLinear = solveLinear,
836 // alphaDiag, dt,
837 // iter, ct, uPos]
838 auto solve_multigrid_impl = [&](TDof &x_upper, const TDof &rhs_upper, const TDof &resOther_upper, int mgLevel, int mgLevelMax, auto &solve_multigrid_impl_ref) -> void
839 {
840 static const int use_1st_conv = 1;
841
842 static const int use_1st_conv_ignore_vis =
843#ifdef USE_MG_O1_NO_VISCOUS
844 1;
845#else
846 0;
847#endif
848 DNDS_assert(mgLevel > 0 && mgLevel <= mgLevelMax);
849 std::string level_key = std::to_string(mgLevel);
850 DNDS_assert(config.linearSolverControl.coarseGridLinearSolverControlList.at(level_key).multiGridNIterPost >= 0);
851 int curMGIter = config.linearSolverControl.coarseGridLinearSolverControlList.at(level_key).multiGridNIter +
852 config.linearSolverControl.coarseGridLinearSolverControlList.at(level_key).multiGridNIterPost;
853 if (curMGIter < 0)
854 curMGIter = config.linearSolverControl.multiGridLPInnerNIter;
855 DNDS_EULER_SOLVER_GET_TEMP_UDOF(uMG1)
856 // DNDS_EULER_SOLVER_GET_TEMP_UDOF(uMG1Init)
857 // DNDS_EULER_SOLVER_GET_TEMP_UDOF(rhsInitCurMG1)
858 DNDS_EULER_SOLVER_GET_TEMP_UDOF(resOtherCurMG)
859 DNDS_EULER_SOLVER_GET_TEMP_UDOF(rhsTemp)
860 // if (mgLevel == 2)
861 // return;
862 uMG1 = x_upper; // projection
863 // std::cout << "here0" << std::endl;
864 resOtherCurMG = rhs_upper; // projection
865 resOtherCurMG *= alphaDiag;
866 resOtherCurMG += resOther_upper; // projection
867 // res0CurMG.addTo(uMG1, -1. / dt);
868 // uMG1Init = uMG1;
869
870 if (config.linearSolverControl.coarseGridLinearSolverControlList.at(level_key).centralSmoothInputResidual > 0)
871 { // make resOtherCurMG <- proj(res of upper grid) and do smoothing
872 resOtherCurMG.addTo(uMG1, -1. / dt);
873 DNDS_EULER_SOLVER_GET_TEMP_UDOF(rhsTemp1)
874 rhsTemp = resOtherCurMG;
875 eval.CentralSmoothResidual(resOtherCurMG, rhsTemp, rhsTemp1, config.linearSolverControl.coarseGridLinearSolverControlList.at(level_key).centralSmoothInputResidual);
876 resOtherCurMG = rhsTemp;
877 resOtherCurMG.addTo(uMG1, 1. / dt);
878 }
879
880 // from uMG1 to rhsTemp - JSourceTmp
881 auto call_evaluate_rhs = [&]()
882 {
883 if (mgLevel == 1)
884 eval.EvaluateRHS(rhsTemp, JSourceTmp, uMG1,
885 config.limiterControl.useViscousLimited ? uRecNew : uRec /*dummy*/, uRec /*dummy*/,
886 betaPP /*dummy*/, alphaPP /*dummy*/, false, tSimu + dt * ct,
887 TEval::RHS_Direct_2nd_Rec | TEval::RHS_Dont_Record_Bud_Flux | TEval::RHS_Dont_Update_Integration |
888 TEval::RHS_Direct_2nd_Rec_already_have_uGradBufNoLim | //! uGradBufNoLim already existent in fdtau
889 (config.limiterControl.useLimiter ? TEval::RHS_Direct_2nd_Rec_use_limiter : TEval::RHS_No_Flags) |
890 TEval::RHS_Recover_IncFScale);
891 else if (mgLevel == 2)
892 eval.EvaluateRHS(rhsTemp, JSourceTmp, uMG1,
893 config.limiterControl.useViscousLimited ? uRecNew : uRec /*dummy*/, uRec /*dummy*/,
894 betaPP /*dummy*/, alphaPP /*dummy*/, false, tSimu + dt * ct,
895 (TEval::RHS_Direct_2nd_Rec_1st_Conv * use_1st_conv) |
896 TEval::RHS_Direct_2nd_Rec |
897 TEval::RHS_Dont_Record_Bud_Flux |
898 TEval::RHS_Dont_Update_Integration |
899 (TEval::RHS_Ignore_Viscosity * use_1st_conv_ignore_vis) |
900 TEval::RHS_Direct_2nd_Rec_already_have_uGradBufNoLim | //! uGradBufNoLim already existent in fdtau
901 (config.limiterControl.useLimiter ? TEval::RHS_Direct_2nd_Rec_use_limiter : TEval::RHS_No_Flags) |
902 TEval::RHS_Recover_IncFScale);
903 else
904 DNDS_assert(false);
905 };
906
907 for (int iIterMG = 1; iIterMG <= curMGIter; iIterMG++)
908 {
909
910 // if (curMGIter > 1 && mgLevel == mgLevelMax) // this is used for checking lusgs-1lusgs == 2xlusgs
911 {
912 if (iIterMG > 1)
913 fdtau(uMG1, dTauC, alphaDiag, uPos); //! warning! dTauC is overwritten
914 call_evaluate_rhs();
915 rhsTemp.trans.startPersistentPull();
916 rhsTemp.trans.waitPersistentPull();
917 if (iIterMG == 1)
918 {
919 // rhsInitCurMG1 = rhsTemp;
920 resOtherCurMG.addTo(rhsTemp, -alphaDiag);
921 }
922 // if (mgLevel < mgLevelMax && iIterMG == 1) // pre smoother coarser grid call
923 // solve_multigrid_impl_ref(uMG1, rhsTemp, resOtherCurMG, mgLevel + 1, mgLevelMax, solve_multigrid_impl_ref);
924 rhsTemp *= alphaDiag;
925 rhsTemp += resOtherCurMG;
926 rhsTemp.addTo(uMG1, -1. / dt);
927 // todo: add rhsfpphere
928 // rhsTemp === alphaDiag * rhs(cur) - alphaDiag * rhs(at_step_1) - uMG1 / dt + uMG1Init / dt
929 }
930 // else
931 // {
932 // rhsTemp = resOtherCurMG;
933 // rhsTemp.addTo(uMG1, -1. / dt);
934 // }
935
936 eval.LUSGSMatrixInit(JDTmp, JSourceTmp, dTauC, dt, alphaDiag, uMG1, uRecNew, 0, tSimu);
937
938 if (iIterMG % config.linearSolverControl.multiGridLPInnerNSee == 0)
939 {
940 Eigen::VectorFMTSafe<real, -1> resNorm;
941 eval.EvaluateNorm(resNorm, rhsTemp);
942 if (mpi.rank == 0)
943 log() << fmt::format("MG Level LP [{}] iter [{}] res [{:.3e}]", mgLevel, iIterMG, resNorm.transpose()) << std::endl;
944 }
945 xIncBuf.setConstant(0.0);
946 solveLinear(alphaDiag, tSimu, rhsTemp, uMG1, xIncBuf, uRecNew, uRecNew,
947 JDTmp, *gmres, mgLevel);
948 fincrement(uMG1, xIncBuf, 1.0, uPos);
949 // solve_multigrid_impl(x_base, cxInc, mgLevel + 1, mgLevelMax);
950 // cxInc *= -1;
951 // std::cout << "here" << std::endl;
952
953 if (mgLevel < mgLevelMax &&
954 iIterMG == config.linearSolverControl.coarseGridLinearSolverControlList.at(level_key).multiGridNIter) // post smoother coarser grid call
955 {
956 call_evaluate_rhs();
957 solve_multigrid_impl_ref(uMG1, rhsTemp, resOtherCurMG, mgLevel + 1, mgLevelMax, solve_multigrid_impl_ref);
958 }
959 }
960
961 x_upper = uMG1; // interpolate
962
963 // alphaDiag *rhsTemp(x + xinc) - xinc / dt == alphaDiag *rhsTemp(x) - res_of_first
964 };
965
966 fdtau(x_upper, dTauC, alphaDiag, uPos); // warning: fdtau resets lambda01234, crucial if useRoeJacobian
967 frhs(rhsBuf, x_upper, dTauC, iter, ct, uPos);
968 rhsBuf.trans.startPersistentPull();
969 rhsBuf.trans.waitPersistentPull();
970 solve_multigrid_impl(x_upper, rhsBuf, resOther, mgLevelInit, mgLevelMax, solve_multigrid_impl);
971 };
972
973 if (config.linearSolverControl.multiGridLP >= 1 && iter > config.linearSolverControl.multiGridLPStartIter)
974 {
975 DNDS_assert(config.linearSolverControl.multiGridLP <= 2);
976 DNDS_EULER_SOLVER_GET_TEMP_UDOF(cxTemp)
977 DNDS_EULER_SOLVER_GET_TEMP_UDOF(resTemp)
978 cxTemp = cx;
979 fincrement(cxTemp, cxInc, 1.0, uPos);
980 // eval.settings.useRoeJacobian = true;
981 //! overwrites cxInc, cxTemp and resTemp do not overwrite cres! (as we might use cres for evaluation of convergence)
982 solve_multigrid(cxTemp, cxInc, resTemp, resOther, 1, config.linearSolverControl.multiGridLP);
983 // eval.settings.useRoeJacobian = false;
984 cxInc = cxTemp;
985 cxInc -= cx; // TODO: renew fsolve to produce cxNew instead of cxInc!!!
986 }
987 // eval.FixIncrement(cx, cxInc);
988 // !freeze something
989 if (getNVars(model) > I4 + 1 && iter <= config.others.nFreezePassiveInner)
990 {
991 for (int i = 0; i < cres.Size(); i++)
992 cxInc[i](Eigen::seq(I4 + 1, EigenLast)).setZero();
993 // if (mpi.rank == 0)
994 // std::cout << "Freezing all passive" << std::endl;
995 }
996 };
997
998 auto fsolveNest = [&](
999 ArrayDOFV<nVarsFixed> &cx,
1000 ArrayDOFV<nVarsFixed> &cx1,
1001 ArrayDOFV<nVarsFixed> &crhs,
1002 ArrayDOFV<1> &dTau,
1003 const std::vector<real> &Coefs, // coefs are dU * c[0] + dt * c[1] * (I/(dt * c[2]) - JMid) * (I/(dt * c[3]) - J)
1004 real dt, real alphaDiag, ArrayDOFV<nVarsFixed> &cxInc, int iter, int uPos)
1005 {
1006 {
1007 DNDS_EULER_SOLVER_GET_TEMP_UDOF(uTemp)
1008 DNDS_EULER_SOLVER_GET_TEMP_UDOF(rhsTemp)
1009
1010 crhs.trans.startPersistentPull();
1011 crhs.trans.waitPersistentPull();
1012 rhsTemp = crhs;
1013 eval.CentralSmoothResidual(rhsTemp, crhs, uTemp);
1014 }
1015
1016 cxInc.setConstant(0.0);
1017 auto &JDC = config.timeMarchControl.timeMarchIsTwoStage() && uPos == 1 ? JD1 : JD;
1018 auto &JSourceC = config.timeMarchControl.timeMarchIsTwoStage() && uPos == 1 ? JSource1 : JSource;
1019 auto &uRecC = config.timeMarchControl.timeMarchIsTwoStage() && uPos == 1 ? uRec1 : uRec;
1020 // TODO: use "update spectral radius" procedure? or force update in fsolve
1021 eval.EvaluateDt(dTau, cx1, uRecC, CFLNow, curDtMin, 1e100, config.implicitCFLControl.useLocalDt, tSimu);
1022 dTau *= Coefs[2];
1023 eval.LUSGSMatrixInit(JD1, JSource1,
1024 dTau, dt * Coefs[2], alphaDiag,
1025 cx1, uRec,
1026 0,
1027 tSimu);
1028 eval.EvaluateDt(dTau, cx, uRecC, CFLNow, curDtMin, 1e100, config.implicitCFLControl.useLocalDt, tSimu);
1029 dTau *= Coefs[3] * veryLargeReal;
1030 eval.LUSGSMatrixInit(JD, JSource,
1031 dTau, dt * Coefs[3], alphaDiag,
1032 cx, uRec,
1033 0,
1034 tSimu);
1035
1036 // for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
1037 // {
1038 // crhs[iCell] = eval.CompressInc(cx[iCell], crhs[iCell] * dTau[iCell]) / dTau[iCell];
1039 // }
1040
1041 if (config.linearSolverControl.gmresCode == 0 || config.linearSolverControl.gmresCode == 2)
1042 {
1043 // //! LUSGS
1044 DNDS_EULER_SOLVER_GET_TEMP_UDOF(uTemp)
1045
1046 eval.UpdateLUSGSForward(alphaDiag, tSimu, crhs, cx1, cxInc, JD1, cxInc);
1047 cxInc.trans.startPersistentPull();
1048 cxInc.trans.waitPersistentPull();
1049 eval.UpdateLUSGSBackward(alphaDiag, tSimu, crhs, cx1, cxInc, JD1, cxInc);
1050 cxInc.trans.startPersistentPull();
1051 cxInc.trans.waitPersistentPull();
1052 uTemp = cxInc;
1053 eval.UpdateLUSGSForward(alphaDiag, tSimu, uTemp, cx, cxInc, JD, cxInc);
1054 cxInc.trans.startPersistentPull();
1055 cxInc.trans.waitPersistentPull();
1056 eval.UpdateLUSGSBackward(alphaDiag, tSimu, uTemp, cx, cxInc, JD, cxInc);
1057 cxInc.trans.startPersistentPull();
1058 cxInc.trans.waitPersistentPull();
1059 cxInc *= 1. / (dt * Coefs[1]);
1060 }
1061
1062 if (config.linearSolverControl.gmresCode != 0)
1063 {
1064 DNDS_EULER_SOLVER_GET_TEMP_UDOF(uTemp)
1065 // ! GMRES
1066 // ! for gmres solver: A * uinc = rhsinc, rhsinc is average value insdead of cumulated on vol
1067 gmres->solve(
1068 [&](decltype(u) &x, decltype(u) &Ax)
1069 {
1070 eval.LUSGSMatrixVec(alphaDiag, tSimu, cx, x, JD, Ax);
1071 Ax.trans.startPersistentPull();
1072 Ax.trans.waitPersistentPull();
1073 uTemp = Ax;
1074 eval.LUSGSMatrixVec(alphaDiag, tSimu, cx1, uTemp, JD1, Ax);
1075 Ax.trans.startPersistentPull();
1076 Ax.trans.waitPersistentPull();
1077 Ax *= dt * Coefs[1];
1078 Ax.addTo(x, Coefs[0]);
1079 },
1080 [&](decltype(u) &x, decltype(u) &MLx)
1081 {
1082 // x as rhs, and MLx as uinc
1083 eval.UpdateLUSGSForward(alphaDiag, tSimu, x, cx1, MLx, JD1, MLx);
1084 MLx.trans.startPersistentPull();
1085 MLx.trans.waitPersistentPull();
1086 eval.UpdateLUSGSBackward(alphaDiag, tSimu, x, cx1, MLx, JD1, MLx);
1087 MLx.trans.startPersistentPull();
1088 MLx.trans.waitPersistentPull();
1089 uTemp = MLx;
1090 eval.UpdateLUSGSForward(alphaDiag, tSimu, uTemp, cx, MLx, JD, MLx);
1091 MLx.trans.startPersistentPull();
1092 MLx.trans.waitPersistentPull();
1093 eval.UpdateLUSGSBackward(alphaDiag, tSimu, uTemp, cx, MLx, JD, MLx);
1094 MLx.trans.startPersistentPull();
1095 MLx.trans.waitPersistentPull();
1096 MLx *= 1. / (dt * Coefs[1]);
1097 },
1098 [&](decltype(u) &a, decltype(u) &b) -> real
1099 {
1100 return a.dot(b); //! need dim balancing here
1101 },
1102 crhs, cxInc, config.linearSolverControl.nGmresIter,
1103 [&](uint32_t i, real res, real resB) -> bool
1104 {
1105 if (i > 0)
1106 {
1107 if (mpi.rank == 0)
1108 {
1109 // log() << std::scientific;
1110 // log() << "GMRES: " << i << " " << resB << " -> " << res << std::endl;
1111 }
1112 }
1113 return false;
1114 });
1115 }
1116 // eval.FixIncrement(cx, cxInc);
1117 // !freeze something
1118 if (getNVars(model) > I4 + 1 && iter <= config.others.nFreezePassiveInner)
1119 {
1120 for (int i = 0; i < crhs.Size(); i++)
1121 cxInc[i](Eigen::seq(I4 + 1, EigenLast)).setZero();
1122 // if (mpi.rank == 0)
1123 // std::cout << "Freezing all passive" << std::endl;
1124 }
1125 };
1126
1127 auto falphaLimSource = [&](
1128 ArrayDOFV<nVarsFixed> &v,
1129 int uPos)
1130 {
1131 auto &alphaPPC = config.timeMarchControl.timeMarchIsTwoStage() && uPos == 1 ? alphaPP1 : alphaPP;
1132 for (index i = 0; i < v.Size(); i++)
1133 v[i] *= alphaPPC[i](0);
1134 };
1135
1136 auto fresidualIncPP = [&](
1137 ArrayDOFV<nVarsFixed> &cx,
1138 ArrayDOFV<nVarsFixed> &xPrev,
1139 ArrayDOFV<nVarsFixed> &crhs,
1140 ArrayDOFV<nVarsFixed> &rhsIncPart,
1141 const std::function<void()> &renewRhsIncPart,
1142 real ct,
1143 int uPos)
1144 {
1145 Timer().StartTimer(PerformanceTimer::Positivity);
1146 auto &alphaPPC = config.timeMarchControl.timeMarchIsTwoStage() && uPos == 1 ? alphaPP1 : alphaPP;
1147 auto &betaPPC = config.timeMarchControl.timeMarchIsTwoStage() && uPos == 1 ? betaPP1 : betaPP;
1148 auto &uRecC = config.timeMarchControl.timeMarchIsTwoStage() && uPos == 1 ? uRec1 : uRec;
1149 auto &JSourceC = config.timeMarchControl.timeMarchIsTwoStage() && uPos == 1 ? JSource1 : JSource;
1150 renewRhsIncPart(); // un-fixed now
1151 // rhsIncPart.trans.startPersistentPull();
1152 // rhsIncPart.trans.waitPersistentPull(); //seems not needed
1153 eval.EvaluateCellRHSAlpha(xPrev, uRecC, betaPPC, rhsIncPart, alphaPP_tmp, nLimAlpha, minAlpha, 1.,
1154 2, TEval::EvaluateCellRHSAlpha_MIN_IF_NOT_ONE);
1155 alphaPP_tmp.trans.startPersistentPull();
1156 alphaPP_tmp.trans.waitPersistentPull();
1157 if (nLimAlpha)
1158 if (mpi.rank == 0 &&
1159 (config.outputControl.consoleOutputEveryFix == 1 || config.outputControl.consoleOutputEveryFix == 4))
1160 {
1161 log() << std::scientific << std::setprecision(config.outputControl.nPrecisionConsole)
1162 << "PPResidualLimiter: nLimAlpha [" << nLimAlpha << "]"
1163 << " minAlpha[" << minAlpha << "]" << std::endl;
1164 }
1165 alphaPPC = alphaPP_tmp;
1166 if (config.limiterControl.useLimiter || config.limiterControl.usePPRecLimiter)
1167 eval.EvaluateRHS(crhs, JSourceC, cx, config.limiterControl.useViscousLimited ? uRecLimited : uRecC, uRecLimited,
1168 betaPPC, alphaPPC, false, tSimu + ct * curDtImplicit);
1169 else
1170 eval.EvaluateRHS(crhs, JSourceC, cx, uRecC, uRecC,
1171 betaPPC, alphaPPC, false, tSimu + ct * curDtImplicit);
1172 // rhs now last-fixed
1173 crhs.trans.startPersistentPull();
1174 crhs.trans.waitPersistentPull();
1175
1176 /********************************************/
1177 // first fix
1178 {
1179 renewRhsIncPart(); // now last-fixed
1180 // rhsIncPart.trans.startPersistentPull();
1181 // rhsIncPart.trans.waitPersistentPull(); //seems not needed
1182 eval.EvaluateCellRHSAlphaExpansion(xPrev, uRecC, betaPPC, rhsIncPart, alphaPP_tmp, nLimAlpha, minAlpha);
1183 alphaPP_tmp.trans.startPersistentPull();
1184 alphaPP_tmp.trans.waitPersistentPull();
1185 if (nLimAlpha)
1186 if (mpi.rank == 0 &&
1187 (config.outputControl.consoleOutputEveryFix == 1 || config.outputControl.consoleOutputEveryFix == 4))
1188 {
1189 log() << std::scientific << std::setprecision(config.outputControl.nPrecisionConsole)
1190 << "PPResidualLimiter - first expand: nLimAlpha [" << nLimAlpha << "]"
1191 << " minAlpha[" << minAlpha << "]" << std::endl;
1192 }
1193 alphaPPC = alphaPP_tmp;
1194 if (config.limiterControl.useLimiter || config.limiterControl.usePPRecLimiter)
1195 eval.EvaluateRHS(crhs, JSourceC, cx, config.limiterControl.useViscousLimited ? uRecLimited : uRecC, uRecLimited,
1196 betaPPC, alphaPPC, false, tSimu + ct * curDtImplicit);
1197 else
1198 eval.EvaluateRHS(crhs, JSourceC, cx, uRecC, uRecC,
1199 betaPPC, alphaPPC, false, tSimu + ct * curDtImplicit);
1200 crhs.trans.startPersistentPull();
1201 crhs.trans.waitPersistentPull();
1202 }
1203 Timer().StopTimer(PerformanceTimer::Positivity);
1204 };
1205
1206 auto fstop = [&](int iter, ArrayDOFV<nVarsFixed> &cres, int iStep) -> bool
1207 {
1208 return functor_fstop(iter, cres, iStep, runningEnvironment);
1209 };
1210
1211 // fmainloop gets the time-variant residual norm,
1212 // handles the output / log nested loops,
1213 // integrates physical time tsimu
1214 // and finally decides if break time loop
1215 auto fmainloop = [&]() -> bool
1216 {
1217 return functor_fmainloop(runningEnvironment);
1218 };
1219
1220 /***************************************************************************************************************************************
1221 .___ ___. ___ __ .__ __. __ ______ ______ .______
1222 | \/ | / \ | | | \ | | | | / __ \ / __ \ | _ \
1223 | \ / | / ^ \ | | | \| | | | | | | | | | | | | |_) |
1224 | |\/| | / /_\ \ | | | . ` | | | | | | | | | | | | ___/
1225 | | | | / _____ \ | | | |\ | | `----.| `--' | | `--' | | |
1226 |__| |__| /__/ \__\ |__| |__| \__| |_______| \______/ \______/ | _|
1227 ***************************************************************************************************************************************/
1228 // step 0 extra:
1229 if (config.outputControl.dataOutAtInit)
1230 {
1231 PrintData(
1232 config.dataIOControl.outPltName + "_" + output_stamp + "_" + "00000",
1233 "",
1234 [&](index iCell)
1235 { return ode->getLatestRHS()[iCell](0); },
1236 addOutList,
1237 eval, tSimu);
1238 eval.PrintBCProfiles(config.dataIOControl.outPltName + "_" + output_stamp + "_" + "00000",
1239 u, uRec);
1240 }
1241 if (config.outputControl.restartOutAtInit)
1242 {
1243 PrintRestart(config.dataIOControl.getOutRestartName() + "_" + output_stamp + "_" + "00000");
1244 }
1245
1246 for (step = 1; step <= config.timeMarchControl.nTimeStep; step++)
1247 {
1248 DNDS_MPI_InsertCheck(mpi, "Implicit Step");
1249 ifOutT = false;
1250 real curDtImplicitOld = curDtImplicit;
1251 curDtImplicit = config.timeMarchControl.dtImplicit;
1252 CFLNow = config.implicitCFLControl.CFL;
1253 fdtau(u, dTauTmp, 1., 0); // generates a curDtMin / CFLNow value as a CFL=1 dt value
1254 curDtImplicit = std::min(curDtMin / CFLNow * config.timeMarchControl.dtCFLLimitScale, curDtImplicit); // limits dt by CFL
1255
1256 if (config.timeMarchControl.useDtPPLimit)
1257 {
1258 DNDS_EULER_SOLVER_GET_TEMP_UDOF(uTemp)
1259 DNDS_EULER_SOLVER_GET_TEMP_UDOF(rhsTemp)
1260 Timer().StartTimer(PerformanceTimer::PositivityOuter);
1261 dTauTmp.setConstant(curDtImplicit * config.timeMarchControl.dtPPLimitScale); //? used as damper here, appropriate?
1262 frhsOuter(rhsTemp, u, dTauTmp, 1, 0.0, 0, 0); //* trick: use 0th order reconstruction RHS for dt PP limiting
1263 uTemp = u;
1264 rhsTemp *= curDtImplicit * config.timeMarchControl.dtPPLimitScale;
1265 index nLim = 0;
1266 real minLim = 1;
1267 eval.EvaluateCellRHSAlpha(u, uRec, alphaPP, rhsTemp, alphaPP_tmp, nLim, minLim, config.timeMarchControl.dtPPLimitRelax,
1268 1, TEval::EvaluateCellRHSAlpha_DEFAULT); // using compress = 1, only minLim is used as output
1269 if (nLim)
1270 curDtImplicit = std::min(curDtImplicit, minLim * curDtImplicit);
1271 if (curDtImplicitHistory.size() && curDtImplicit > curDtImplicitHistory.back() * config.timeMarchControl.dtIncreaseLimit)
1272 {
1273 curDtImplicit = curDtImplicitHistory.back() * config.timeMarchControl.dtIncreaseLimit;
1274 }
1275 if (mpi.rank == 0 && nLim)
1276 {
1277 log() << "##################################################################" << std::endl;
1278 log() << fmt::format("At Step [{:d}] t [{:.8g}] Changing dt to [{}], using PP", step, tSimu, curDtImplicit) << std::endl;
1279 log() << "##################################################################" << std::endl;
1280 }
1281 Timer().StopTimer(PerformanceTimer::PositivityOuter);
1282 }
1283 curDtImplicit = std::max(curDtImplicit, config.timeMarchControl.dtImplicitMin);
1284
1285 DNDS_assert(curDtImplicit > 0);
1286 if (tSimu + curDtImplicit >= nextTout - curDtImplicit * smallReal) // limits dt by output nodes
1287 {
1288 ifOutT = true;
1289 curDtImplicit = std::max(0.0, nextTout - tSimu);
1290 }
1291
1292 if (config.timeMarchControl.useImplicitPP)
1293 {
1294 switch (config.timeMarchControl.odeCode)
1295 {
1296 case 1: // BDF2
1297 std::dynamic_pointer_cast<ODE::ImplicitBDFDualTimeStep<ArrayDOFV<nVarsFixed>, ArrayDOFV<1>>>(ode)
1298 ->StepPP(
1299 u, uIncBufODE,
1300 frhs,
1301 fdtau,
1302 fsolve,
1303 config.convergenceControl.nTimeStepInternal,
1304 fstop, fincrement,
1305 curDtImplicit + verySmallReal,
1306 falphaLimSource,
1307 fresidualIncPP);
1308 break;
1309 case 102: // VBDFPP
1310 {
1311 index nLimAlpha;
1312 real minAlpha;
1313 auto odeVBDF = std::dynamic_pointer_cast<ODE::ImplicitVBDFDualTimeStep<ArrayDOFV<nVarsFixed>, ArrayDOFV<1>>>(ode);
1314 odeVBDF->LimitDt_StepPPV2(
1315 u, [&](ArrayDOFV<nVarsFixed> &u, ArrayDOFV<nVarsFixed> &uInc) -> real
1316 {
1317 eval.EvaluateCellRHSAlpha(u, uRec, alphaPP, uInc, alphaPP_tmp, nLimAlpha, minAlpha, 1.,
1318 1, TEval::EvaluateCellRHSAlpha_MIN_IF_NOT_ONE); // using compress == 1, only minAlpha is used
1319 return minAlpha; },
1320 curDtImplicit, 2); // curDtImplicit modified
1321 if (curDtImplicit > curDtImplicitOld)
1322 {
1323 dtIncreaseCounter++;
1324 if (dtIncreaseCounter > config.timeMarchControl.dtIncreaseAfterCount)
1325 curDtImplicit = std::min(curDtImplicitOld * config.timeMarchControl.dtIncreaseLimit, curDtImplicit);
1326 else
1327 curDtImplicit = curDtImplicitOld;
1328 }
1329 else
1330 {
1331 dtIncreaseCounter = 0;
1332 }
1333
1334 if (mpi.rank == 0)
1335 {
1336 log() << "##################################################################" << std::endl;
1337 log() << fmt::format("At Step [{:d}] t [{:.8g}] Changing dt to {}", step, tSimu, curDtImplicit) << std::endl;
1338 log() << "##################################################################" << std::endl;
1339 }
1340 odeVBDF->StepPP(
1341 u, uIncBufODE,
1342 frhs,
1343 fdtau,
1344 fsolve,
1345 config.convergenceControl.nTimeStepInternal,
1346 fstop, fincrement,
1347 curDtImplicit + verySmallReal,
1348 falphaLimSource,
1349 fresidualIncPP);
1350 }
1351 break;
1352
1353 default:
1354 DNDS_assert_info(false, "unsupported odeCode for PP!");
1355 break;
1356 }
1357 }
1358 else if (config.timeMarchControl.timeMarchIsTwoStage() && false)
1359 std::dynamic_pointer_cast<ODE::ImplicitHermite3SimpleJacobianDualStep<ArrayDOFV<nVarsFixed>, ArrayDOFV<1>>>(ode)
1360 ->StepNested(
1361 u, uIncBufODE,
1362 frhs,
1363 fdtau,
1364 fsolve, fsolveNest,
1365 config.convergenceControl.nTimeStepInternal,
1366 fstop, fincrement,
1367 curDtImplicit + verySmallReal);
1368 else
1369 ode
1370 ->Step(
1371 u, uIncBufODE,
1372 frhs,
1373 fdtau,
1374 fsolve,
1375 config.convergenceControl.nTimeStepInternal,
1376 fstop, fincrement,
1377 curDtImplicit + verySmallReal);
1378 curDtImplicitHistory.push_back(curDtImplicit);
1379 if (fmainloop())
1380 break;
1381 }
1382 }
1383
1384 DNDS_SWITCH_INTELLISENSE(
1385 // the real definition
1386 template <EulerModel model>
1387 ,
1388 // the intellisense friendly definition
1389 template <>)
1390 /** @brief Dispatch the linear solve for implicit time stepping.
1391 *
1392 * Solves the linearized system using either SGS sweeps (gmresCode=0) or
1393 * FGMRES with preconditioner (gmresCode=1/2). Supports multi-grid levels
1394 * with per-level configuration of iteration counts and solver parameters.
1395 *
1396 * @param alphaDiag Diagonal scaling factor for the implicit operator.
1397 * @param t Current simulation time.
1398 * @param cres RHS residual (right-hand side of the linear system).
1399 * @param cx Current solution DOF.
1400 * @param cxInc Solution increment (input/output).
1401 * @param uRecC Reconstruction coefficients.
1402 * @param uRecIncC Reconstruction increment (for SGS-with-rec mode).
1403 * @param JDC Diagonal Jacobian block.
1404 * @param gmres GMRES solver object.
1405 * @param gridLevel Multigrid level (0 = finest).
1406 */
1407 void EulerSolver<model>::solveLinear(
1408 real alphaDiag, real t,
1409 TDof &cres, TDof &cx, TDof &cxInc, TRec &uRecC, TRec uRecIncC,
1410 JacobianDiagBlock<nVarsFixed> &JDC, tGMRES_u &gmres, int gridLevel)
1411 {
1412 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
1413 auto &eval = *pEval;
1414 auto initUDOF = [&](ArrayDOFV<nVarsFixed> &uu)
1415 { vfv->BuildUDof(uu, nVars); };
1416
1417 TU sgsRes(nVars), sgsRes0(nVars);
1418 bool inputIsZero{true}, hasLUDone{false};
1419
1420 typename TVFV::template TFBoundary<nVarsFixed>
1421 FBoundary = [&](const TU &UL, const TU &UMean, index iCell, index iFace, int iG,
1422 const Geom::tPoint &normOut, const Geom::tPoint &pPhy, const Geom::t_index bType) -> TU
1423 {
1424 TU UR = UL;
1425 UR.setZero();
1426 return UR;
1427 };
1428 DNDS_assert(gridLevel <= config.linearSolverControl.coarseGridLinearSolverControlList.size());
1429 std::string level_key = std::to_string(gridLevel);
1430 auto gmresCode =
1431 gridLevel > 0
1432 ? config.linearSolverControl.coarseGridLinearSolverControlList.at(level_key).gmresCode
1433 : config.linearSolverControl.gmresCode;
1434 bool initWithLastURecInc =
1435 gridLevel > 0
1436 ? false
1437 : config.linearSolverControl.initWithLastURecInc;
1438 int sgsWithRec =
1439 gridLevel > 0
1440 ? 0
1441 : config.linearSolverControl.sgsWithRec;
1442 int sgsIter =
1443 gridLevel > 0
1444 ? config.linearSolverControl.coarseGridLinearSolverControlList.at(level_key).sgsIter
1445 : config.linearSolverControl.sgsIter;
1446 int nSgsConsoleCheck =
1447 gridLevel > 0
1448 ? config.linearSolverControl.coarseGridLinearSolverControlList.at(level_key).nSgsConsoleCheck
1449 : config.linearSolverControl.nSgsConsoleCheck;
1450 int gmresScale =
1451 gridLevel > 0
1452 ? config.linearSolverControl.coarseGridLinearSolverControlList.at(level_key).gmresScale
1453 : config.linearSolverControl.gmresScale;
1454 int nGmresIter =
1455 gridLevel > 0
1456 ? config.linearSolverControl.coarseGridLinearSolverControlList.at(level_key).nGmresIter
1457 : config.linearSolverControl.nGmresIter;
1458 int nGmresConsoleCheck =
1459 gridLevel > 0
1460 ? config.linearSolverControl.coarseGridLinearSolverControlList.at(level_key).nGmresConsoleCheck
1461 : config.linearSolverControl.nGmresConsoleCheck;
1462 if (gmresCode == 0 || gmresCode == 2)
1463 {
1464 // //! LUSGS
1465
1466 if (initWithLastURecInc)
1467 {
1468 DNDS_assert(config.implicitReconstructionControl.storeRecInc);
1469 eval.UpdateSGSWithRec(alphaDiag, t, cres, cx, uRecC, cxInc, uRecIncC, JDC, true, sgsRes);
1470 // for (index iCell = 0; iCell < uRecIncC.Size(); iCell++)
1471 // std::cout << "-------\n"
1472 // << uRecIncC[iCell] << std::endl;
1473 cxInc.trans.startPersistentPull();
1474 cxInc.trans.waitPersistentPull();
1475 eval.UpdateSGSWithRec(alphaDiag, t, cres, cx, uRecC, cxInc, uRecIncC, JDC, false, sgsRes);
1476 cxInc.trans.startPersistentPull();
1477 cxInc.trans.waitPersistentPull();
1478 }
1479 else
1480 {
1481 DNDS_EULER_SOLVER_GET_TEMP_UDOF(uTemp)
1482 doPrecondition(alphaDiag, t, cres, cx, cxInc, uTemp, JDC, sgsRes, inputIsZero, hasLUDone, gridLevel);
1483 }
1484
1485 if (sgsWithRec != 0)
1486 uRecNew.setConstant(0.0);
1487 for (int iterSGS = 1; iterSGS <= sgsIter; iterSGS++)
1488 {
1489 if (sgsWithRec != 0)
1490 {
1491 vfv->DoReconstructionIter(
1492 uRecNew, uRecNew1, cxInc,
1493 FBoundary,
1494 false);
1495 uRecNew.trans.startPersistentPull();
1496 uRecNew.trans.waitPersistentPull();
1497 eval.UpdateSGSWithRec(alphaDiag, t, cres, cx, uRecC, cxInc, uRecNew, JDC, true, sgsRes);
1498 cxInc.trans.startPersistentPull();
1499 cxInc.trans.waitPersistentPull();
1500 eval.UpdateSGSWithRec(alphaDiag, t, cres, cx, uRecC, cxInc, uRecNew, JDC, false, sgsRes);
1501 cxInc.trans.startPersistentPull();
1502 cxInc.trans.waitPersistentPull();
1503 }
1504 else
1505 {
1506 DNDS_EULER_SOLVER_GET_TEMP_UDOF(uTemp)
1507 doPrecondition(alphaDiag, t, cres, cx, cxInc, uTemp, JDC, sgsRes, inputIsZero, hasLUDone, gridLevel);
1508 }
1509 if (iterSGS == 1)
1510 sgsRes0 = sgsRes;
1511
1512 if (mpi.rank == 0 && iterSGS % nSgsConsoleCheck == 0)
1513 log() << std::scientific << "SGS1 " << std::to_string(iterSGS)
1514 << " [" << sgsRes0.transpose() << "] ->"
1515 << " [" << sgsRes.transpose() << "] " << std::endl;
1516 }
1517 }
1518 Eigen::VectorXd meanScale;
1519 if (gmresScale == 1)
1520 {
1521 meanScale = eval.settings.refU;
1522 meanScale(Seq123).setConstant(std::sqrt(meanScale(0) * meanScale(I4))); //! using consistent rho U scale
1523 // meanScale(I4) = sqr(meanScale(1)) / (meanScale(0) + verySmallReal);
1524 // meanScale(0) = 0.01;
1525 // meanScale(Seq123).setConstant(0.1);
1526 // meanScale(I4) = 1;
1527 }
1528 else if (gmresScale == 2)
1529 {
1530 eval.EvaluateNorm(meanScale, cx, 1, true, true);
1531 meanScale(Seq123).setConstant(meanScale(Seq123).norm());
1532 meanScale(I4) = sqr(meanScale(1)) / (meanScale(0) + verySmallReal);
1533 }
1534 else
1535 meanScale.setOnes(nVars);
1536 // meanScale(0) = 10;
1537 TU meanScaleInv = (meanScale.array() + verySmallReal).inverse();
1538
1539 if (gmresCode != 0)
1540 {
1541 DNDS_EULER_SOLVER_GET_TEMP_UDOF(uTemp)
1542 // ! GMRES
1543 // ! for gmres solver: A * uinc = rhsinc, rhsinc is average value insdead of cumulated on vol
1544 gmres.solve(
1545 [&](TDof &x, TDof &Ax)
1546 {
1547 eval.LUSGSMatrixVec(alphaDiag, t, cx, x, JDC, Ax);
1548 Ax.trans.startPersistentPull();
1549 Ax.trans.waitPersistentPull();
1550 },
1551 [&](TDof &x, TDof &MLx)
1552 {
1553 // x as rhs, and MLx as uinc
1554 MLx.setConstant(0.0), inputIsZero = true; //! start as zero
1555 doPrecondition(alphaDiag, t, x, cx, MLx, uTemp, JDC, sgsRes, inputIsZero, hasLUDone, gridLevel);
1556 for (int i = 0; i < sgsIter; i++)
1557 {
1558 doPrecondition(alphaDiag, t, x, cx, MLx, uTemp, JDC, sgsRes, inputIsZero, hasLUDone, gridLevel);
1559 }
1560 },
1561 [&](TDof &a, TDof &b) -> real
1562 {
1563 return a.dot(b, meanScaleInv.array(), meanScaleInv.array());
1564 },
1565 cres, cxInc, nGmresIter,
1566 [&](uint32_t i, real res, real resB) -> bool
1567 {
1568 if (i > 0 && i % nGmresConsoleCheck == 0)
1569 {
1570 if (mpi.rank == 0)
1571 {
1572 log() << std::scientific;
1573 log() << "GMRES: " << i << " " << resB << " -> " << res << std::endl;
1574 }
1575 }
1576 return false;
1577 });
1578 }
1579 }
1580
1581 DNDS_SWITCH_INTELLISENSE(
1582 // the real definition
1583 template <EulerModel model>
1584 ,
1585 // the intellisense friendly definition
1586 template <>)
1587 /** @brief Apply the preconditioner for the GMRES linear solver.
1588 *
1589 * Depending on jacobiCode: applies symmetric SGS sweeps (code 0/1) or
1590 * local LU factorization solve (code 2). Manages the inputIsZero and
1591 * hasLUDone state flags across GMRES restarts.
1592 *
1593 * @param alphaDiag Diagonal scaling factor.
1594 * @param t Current simulation time.
1595 * @param cres Right-hand side for preconditioning.
1596 * @param cx Current solution DOF.
1597 * @param cxInc Preconditioned result (output).
1598 * @param uTemp Temporary DOF buffer.
1599 * @param JDC Diagonal Jacobian block.
1600 * @param sgsRes SGS residual for convergence tracking.
1601 * @param inputIsZero Tracks whether cxInc is zero (updated in-place).
1602 * @param hasLUDone Tracks whether LU factorization has been computed (updated).
1603 * @param gridLevel Multigrid level.
1604 */
1605 void EulerSolver<model>::doPrecondition(real alphaDiag, real t,
1606 TDof &cres, TDof &cx, TDof &cxInc, TDof &uTemp,
1607 JacobianDiagBlock<nVarsFixed> &JDC, TU &sgsRes, bool &inputIsZero, bool &hasLUDone, int gridLevel)
1608 {
1609 DNDS_assert(pEval);
1610 auto &eval = *pEval;
1611 auto initUDOF = [&](ArrayDOFV<nVarsFixed> &uu)
1612 { vfv->BuildUDof(uu, nVars); };
1613 // static int nCall{0};
1614 // nCall++;
1615 // if (mpi.rank == 0)
1616 // std::cout << "doPrecondition nCall " << nCall << fmt::format(" {} ", hasLUDone) << std::endl;
1617 std::string level_key = std::to_string(gridLevel);
1618 int jacobiCode =
1619 gridLevel > 0
1620 ? config.linearSolverControl.coarseGridLinearSolverControlList.at(level_key).jacobiCode
1621 : config.linearSolverControl.jacobiCode;
1622
1623 if (jacobiCode <= 1)
1624 {
1625 bool useGS = jacobiCode == 1;
1626 eval.UpdateSGS(alphaDiag, t, cres, cx, cxInc, uTemp, JDC, true, useGS, sgsRes);
1627 cxInc = uTemp;
1628 cxInc.trans.startPersistentPull();
1629 cxInc.trans.waitPersistentPull();
1630 eval.UpdateSGS(alphaDiag, t, cres, cx, cxInc, uTemp, JDC, false, useGS, sgsRes);
1631 cxInc = uTemp;
1632 cxInc.trans.startPersistentPull();
1633 cxInc.trans.waitPersistentPull();
1634 // eval.UpdateLUSGSForward(alphaDiag, cres, cx, cxInc, JDC, cxInc);
1635 // cxInc.trans.startPersistentPull();
1636 // cxInc.trans.waitPersistentPull();
1637 // eval.UpdateLUSGSBackward(alphaDiag, cres, cx, cxInc, JDC, cxInc);
1638 // cxInc.trans.startPersistentPull();
1639 // cxInc.trans.waitPersistentPull();
1640 inputIsZero = false;
1641 }
1642 else if (jacobiCode == 2)
1643 {
1644 DNDS_EULER_SOLVER_GET_TEMP_UDOF(rhsTemp)
1645 DNDS_assert_info(config.linearSolverControl.directPrecControl.useDirectPrec, "need to use config.linearSolverControl.directPrecControl.useDirectPrec first !");
1646 if (!hasLUDone)
1647 eval.LUSGSMatrixToJacobianLU(alphaDiag, t, cx, JDC, *JLocalLU), hasLUDone = true;
1648 for (int iii = 0; iii < 2; iii++)
1649 {
1650 eval.LUSGSMatrixSolveJacobianLU(alphaDiag, t, cres, cx, cxInc, uTemp, rhsTemp, JDC, *JLocalLU, inputIsZero, sgsRes);
1651 uTemp.SwapDataFatherSon(cxInc);
1652 // cxInc = uTemp;
1653 cxInc.trans.startPersistentPull();
1654 cxInc.trans.waitPersistentPull();
1655 inputIsZero = false;
1656 }
1657 }
1658 }
1659
1660 DNDS_SWITCH_INTELLISENSE(
1661 // the real definition
1662 template <EulerModel model>
1663 ,
1664 // the intellisense friendly definition
1665 template <>)
1666 /** @brief Initialize the running environment for the time-marching loop.
1667 *
1668 * Sets up the ODE solver, error logging stream, temporary arrays, signal handlers,
1669 * output data structures, and applies steady-mode overrides when configured.
1670 * Must be called before entering the main time-stepping loop.
1671 *
1672 * @param runningEnvironment Running environment structure to populate (output).
1673 */
1674 void EulerSolver<model>::InitializeRunningEnvironment(EulerSolver<model>::RunningEnvironment &runningEnvironment)
1675 {
1676 if (config.timeMarchControl.partitionMeshOnly)
1677 {
1678 return;
1679 }
1680 // mind we need to get ptr-to-actual-eval into env.pEval,
1681 // before assigning the ref (ptr), or the ptr is null
1682 runningEnvironment.pEval = pEval;
1683 DNDS_EULERSOLVER_RUNNINGENV_GET_REF_LIST
1684
1685 /*******************************************************/
1686 /* LOGERR STREAM */
1687 /*******************************************************/
1688 logErr = LogErrInitialize();
1689
1690 /*******************************************************/
1691 /* PICK ODE SOLVER */
1692 /*******************************************************/
1693
1694 auto buildDOF = [&](ArrayDOFV<nVarsFixed> &data)
1695 {
1696 vfv->BuildUDof(data, nVars);
1697 };
1698 auto buildScalar = [&](ArrayDOFV<1> &data)
1699 {
1700 vfv->BuildUDof(data, 1);
1701 };
1702 // std::cout << fmt::format("nVars {}, here50", nVars);
1703 if (config.timeMarchControl.steadyQuit)
1704 {
1705 if (mpi.rank == 0)
1706 log() << "Using steady!" << std::endl;
1707 config.timeMarchControl.odeCode = 1; // To bdf;
1708 config.timeMarchControl.nTimeStep = 1;
1709 config.timeMarchControl.dtImplicit = 1e100;
1710 config.timeMarchControl.useDtPPLimit = false;
1711 config.timeMarchControl.dtCFLLimitScale = 1e110;
1712 config.outputControl.tDataOut = 1e300; // no t-out steps
1713 }
1714 switch (config.timeMarchControl.odeCode)
1715 {
1716 case 0: // esdirk4
1717 if (mpi.rank == 0)
1718 log() << "=== ODE: ESDIRK4 " << std::endl;
1719 ode = std::make_shared<ODE::ImplicitSDIRK4DualTimeStep<ArrayDOFV<nVarsFixed>, ArrayDOFV<1>>>(
1720 mesh->NumCell(),
1721 buildDOF, buildScalar,
1722 1); // 1 for esdirk
1723 break;
1724 case 101: // sdirk4
1725 if (mpi.rank == 0)
1726 log() << "=== ODE: SSP-SDIRK4 " << std::endl;
1727 ode = std::make_shared<ODE::ImplicitSDIRK4DualTimeStep<ArrayDOFV<nVarsFixed>, ArrayDOFV<1>>>(
1728 mesh->NumCell(),
1729 buildDOF, buildScalar,
1730 0);
1731 break;
1732 case 202: // esdirk3
1733 case 203: // trapz
1734 case 204: // esdirk2
1735 if (mpi.rank == 0)
1736 log() << "=== ODE: " + std::vector<std::string>{"ESDIRK3", "Trapz", "ESDIRK2"}.at(config.timeMarchControl.odeCode - 202) << std::endl;
1737 ode = std::make_shared<ODE::ImplicitSDIRK4DualTimeStep<ArrayDOFV<nVarsFixed>, ArrayDOFV<1>>>(
1738 mesh->NumCell(),
1739 buildDOF, buildScalar,
1740 config.timeMarchControl.odeCode - 200); // 202 -> 2
1741 break;
1742 case 1: // BDF2
1743 case 103: // Backward Euler
1744 if (mpi.rank == 0 && config.timeMarchControl.odeCode == 1)
1745 log() << "=== ODE: BDF2 " << std::endl;
1746 if (mpi.rank == 0 && config.timeMarchControl.odeCode == 103)
1747 log() << "=== ODE: Backward Euler " << std::endl;
1748 ode = std::make_shared<ODE::ImplicitVBDFDualTimeStep<ArrayDOFV<nVarsFixed>, ArrayDOFV<1>>>(
1749 mesh->NumCell(),
1750 buildDOF, buildScalar,
1751 config.timeMarchControl.odeCode == 1 ? 2 : 1);
1752 break;
1753 case 102: // VBDF2
1754 if (mpi.rank == 0)
1755 log() << "=== ODE: VBDF2 " << std::endl;
1756 ode = std::make_shared<ODE::ImplicitVBDFDualTimeStep<ArrayDOFV<nVarsFixed>, ArrayDOFV<1>>>(
1757 mesh->NumCell(),
1758 buildDOF, buildScalar,
1759 2);
1760 break;
1761 case 2: // SSPRK
1762 if (mpi.rank == 0)
1763 log() << "=== ODE: SSPRK4 " << std::endl;
1764 ode = std::make_shared<ODE::ExplicitSSPRK3TimeStepAsImplicitDualTimeStep<ArrayDOFV<nVarsFixed>, ArrayDOFV<1>>>(
1765 mesh->NumCell(),
1766 buildDOF, buildScalar,
1767 false); // TODO: add local stepping options
1768 break;
1769 case 401: // H3S
1770 DNDS_assert(config.timeMarchControl.timeMarchIsTwoStage());
1771 if (mpi.rank == 0)
1772 log() << "=== ODE: Hermite3 (simple jacobian) " << std::endl;
1773 ode = std::make_shared<ODE::ImplicitHermite3SimpleJacobianDualStep<ArrayDOFV<nVarsFixed>, ArrayDOFV<1>>>(
1774 mesh->NumCell(),
1775 buildDOF, buildScalar,
1776 config.timeMarchControl.odeSetting1 == 0 ? 0.55 : config.timeMarchControl.odeSetting1,
1777 std::round(config.timeMarchControl.odeSetting2), // Backward Euler Starter
1778 0, // method
1779 config.timeMarchControl.odeSetting3 == 0 ? 0.9146 : config.timeMarchControl.odeSetting3, // thetaM1
1780 0.0, // thetaM2 = 0
1781 std::round(config.timeMarchControl.odeSetting4) // mask
1782 );
1783 break;
1784 case 411: // H3S U2R2
1785 case 412: // H3S U2R1
1786 case 413: // H3S U3R1
1787 DNDS_assert(config.timeMarchControl.timeMarchIsTwoStage());
1788 if (mpi.rank == 0)
1789 log() << "=== ODE: Hermite3 (simple jacobian) " + fmt::format("Mask {}", config.timeMarchControl.odeCode - 411)
1790 << std::endl;
1791 ode = std::make_shared<ODE::ImplicitHermite3SimpleJacobianDualStep<ArrayDOFV<nVarsFixed>, ArrayDOFV<1>>>(
1792 mesh->NumCell(),
1793 buildDOF, buildScalar,
1794 config.timeMarchControl.odeSetting1 == 0 ? 0.55 : config.timeMarchControl.odeSetting1,
1795 std::round(config.timeMarchControl.odeSetting2), // Backward Euler Starter
1796 0, // method
1797 config.timeMarchControl.odeSetting3 == 0 ? 0.9146 : config.timeMarchControl.odeSetting3, // thetaM1
1798 config.timeMarchControl.odeSetting4 == 0 ? 0.0 : config.timeMarchControl.odeSetting4, // thetaM2
1799 config.timeMarchControl.odeCode - 411 // mask
1800 );
1801 break;
1802 default:
1803 DNDS_assert_info(false, "no such ode code");
1804 }
1805 if (config.timeMarchControl.odeSettingsExtra.is_object())
1806 ode->SetExtraParams(config.timeMarchControl.odeSettingsExtra);
1807
1808 if (config.timeMarchControl.useImplicitPP)
1809 {
1810 DNDS_assert(config.timeMarchControl.odeCode == 1 || config.timeMarchControl.odeCode == 102);
1811 }
1812
1813 // std::cout << fmt::format("nVars {}, here100", nVars);
1814
1815 /*******************************************************/
1816 /* INIT GMRES AND PCG */
1817 /*******************************************************/
1818
1819 if (config.linearSolverControl.gmresCode == 1 ||
1820 config.linearSolverControl.gmresCode == 2)
1821 gmres = std::make_unique<tGMRES_u>(
1822 config.linearSolverControl.nGmresSpace,
1823 [&](decltype(u) &data)
1824 {
1825 vfv->BuildUDof(data, nVars);
1826 });
1827
1828 if (config.implicitReconstructionControl.recLinearScheme == 1)
1829 gmresRec = std::make_unique<tGMRES_uRec>(
1830 config.implicitReconstructionControl.nGmresSpace,
1831 [&](decltype(uRec) &data)
1832 {
1833 vfv->BuildURec(data, nVars);
1834 });
1835
1836 if (config.implicitReconstructionControl.recLinearScheme == 2)
1837 pcgRec = std::make_unique<tPCG_uRec>(
1838 [&](decltype(uRec) &data)
1839 {
1840 vfv->BuildURec(data, nVars);
1841 });
1842
1843 if (config.implicitReconstructionControl.recLinearScheme == 2 || config.timeMarchControl.timeMarchIsTwoStage())
1844 pcgRec1 = std::make_unique<tPCG_uRec>(
1845 [&](decltype(uRec) &data)
1846 {
1847 vfv->BuildURec(data, nVars);
1848 });
1849
1850 /*******************************************************/
1851 /* SOLVER MAJOR START */
1852 /*******************************************************/
1853 tstart = MPI_Wtime();
1854 tstartInternal = tstart;
1855 stepCount = 0;
1856
1857 resBaseC.resize(nVars);
1858 resBaseCInternal.resize(nVars);
1859 resBaseC.setConstant(config.convergenceControl.res_base);
1860
1861 tSimu = 0.0;
1862 tAverage = 0.0;
1863 nextTout = std::min(config.outputControl.tDataOut, config.timeMarchControl.tEnd); // ensures the destination time output
1864 nextStepOut = config.outputControl.nDataOut;
1865 nextStepOutC = config.outputControl.nDataOutC;
1866 nextStepRestart = config.outputControl.nRestartOut;
1867 nextStepRestartC = config.outputControl.nRestartOutC;
1868 nextStepOutAverage = config.outputControl.nTimeAverageOut;
1869 nextStepOutAverageC = config.outputControl.nTimeAverageOutC;
1870
1871 CFLNow = config.implicitCFLControl.CFL;
1872 ifOutT = false;
1873 curDtMin = veryLargeReal;
1874 curDtImplicit = config.timeMarchControl.dtImplicit;
1875 step = 0;
1876 gradIsZero = true;
1877
1878 dtIncreaseCounter = 0;
1879 }
1880}
DNDS_MAKE_SSP
#define DNDS_MAKE_SSP(ssp,...)
Definition Defines.hpp:212
EigenUtil.hpp
Eigen extensions: to_string, an fmt-safe wrapper, and fmt formatter specialisations for dense Eigen m...
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
EulerSolver.hpp
Top-level solver orchestrator for compressible Navier-Stokes / Euler simulations.
DNDS_EULERSOLVER_RUNNINGENV_GET_REF_LIST
#define DNDS_EULERSOLVER_RUNNINGENV_GET_REF_LIST
Definition EulerSolver.hpp:1300
DNDS_EULER_SOLVER_GET_TEMP_UDOF
#define DNDS_EULER_SOLVER_GET_TEMP_UDOF(name)
DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
#define DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
Declares compile-time Eigen index sequences for sub-vector access into conservative state vectors.
Definition Euler.hpp:56
DNDS_MPI_InsertCheck
#define DNDS_MPI_InsertCheck(mpi, info)
Debug helper: barrier + print "CHECK <info> RANK r @ fn @ file:line".
Definition MPI.hpp:320
DNDS_SWITCH_INTELLISENSE
#define DNDS_SWITCH_INTELLISENSE(real, intellisense)
Definition Macros.hpp:86
SpecialFields.hpp
Analytic isentropic-vortex solutions for inviscid accuracy verification.
DNDS::Euler::EulerSolver::RunImplicitEuler
void RunImplicitEuler()
Run the main implicit time-marching loop.
DNDS::Euler::EulerSolver::solveLinear
void solveLinear(real alphaDiag, real t, TDof &cres, TDof &cx, TDof &cxInc, TRec &uRecC, TRec uRecIncC, JacobianDiagBlock< nVarsFixed > &JDC, tGMRES_u &gmres, int gridLevel)
Solve the implicit linear system at a given grid level.
DNDS::Euler::EulerSolver::InitializeRunningEnvironment
void InitializeRunningEnvironment(RunningEnvironment &env)
Populate a RunningEnvironment with allocated solvers, loggers, and initial state.
DNDS::Euler::EulerSolver::doPrecondition
void doPrecondition(real alphaDiag, real t, TDof &crhs, TDof &cx, TDof &cxInc, TDof &uTemp, JacobianDiagBlock< nVarsFixed > &JDC, TU &sgsRes, bool &inputIsZero, bool &hasLUDone, int gridLevel)
Apply the preconditioner (SGS sweeps or ILU) to a right-hand side vector.
DNDS::PerformanceTimer::clearAllTimer
void clearAllTimer()
Zero every timer slot.
Definition Profiling.cpp:69
DNDS::PerformanceTimer::Positivity
@ Positivity
Positivity preservation.
Definition Profiling.hpp:45
DNDS::PerformanceTimer::Limiter
@ Limiter
Slope / variable limiter.
Definition Profiling.hpp:33
DNDS::PerformanceTimer::PositivityOuter
@ PositivityOuter
Outer-iteration positivity.
Definition Profiling.hpp:46
DNDS::PerformanceTimer::LimiterA
@ LimiterA
Limiter sub-phase A.
Definition Profiling.hpp:34
DNDS::PerformanceTimer::Reconstruction
@ Reconstruction
Variational reconstruction.
Definition Profiling.hpp:31
DNDS::PerformanceTimer::RHS
@ RHS
Total RHS evaluation.
Definition Profiling.hpp:29
DNDS::PerformanceTimer::StartTimer
void StartTimer(TimerType t)
Record the current wall time in the "start" slot for timer t.
Definition Profiling.cpp:15
DNDS::PerformanceTimer::StopTimer
void StopTimer(TimerType t)
Add (now - start) to the accumulated time for timer t.
Definition Profiling.cpp:25
DNDS::EulerP::TU
Eigen::Vector< real, nVarsFlow > TU
Fixed-size 5-component conservative state vector (rho, rhoU, rhoV, rhoW, E).
Definition EulerP.hpp:33
DNDS::Euler::NS_2EQ_3D
@ NS_2EQ_3D
NS + 2-equation RANS, 3D geometry (7 vars).
Definition Euler.hpp:882
DNDS::Euler::NS_SA_3D
@ NS_SA_3D
NS + Spalart-Allmaras, 3D geometry (6 vars).
Definition Euler.hpp:880
DNDS::Euler::NS_SA
@ NS_SA
NS + Spalart-Allmaras, 2D geometry (6 vars).
Definition Euler.hpp:877
DNDS::Euler::NS_2EQ
@ NS_2EQ
NS + 2-equation RANS, 2D geometry (7 vars).
Definition Euler.hpp:881
DNDS::Geom::t_index
int32_t t_index
Definition Geometric.hpp:6
DNDS::Geom::tPoint
Eigen::Vector3d tPoint
Definition Geometric.hpp:9
DNDS::TermColor::Reset
constexpr std::string_view Reset
ANSI escape: reset all attributes.
Definition Defines.hpp:883
DNDS::TermColor::Red
constexpr std::string_view Red
ANSI escape: bright red foreground.
Definition Defines.hpp:869
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::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
DNDS::smallReal
DNDS_CONSTANT const real smallReal
Loose lower bound (for iterative-solver tolerances etc.).
Definition Defines.hpp:196
DNDS::log
std::ostream & log()
Return the current DNDSR log stream (either std::cout or the installed file).
Definition Defines.cpp:49
DNDS::Timer
PerformanceTimer & Timer()
Short-hand accessor to the PerformanceTimer singleton.
Definition Profiling.hpp:136
DNDS::veryLargeReal
DNDS_CONSTANT const real veryLargeReal
Catch-all upper bound ("practically infinity") for physical scalars.
Definition Defines.hpp:190
test_basic_fv.mpi
mpi
Definition test_basic_fv.py:487
Eigen::MatrixFMTSafe
Eigen::Matrix wrapper that hides begin/end from fmt.
Definition EigenUtil.hpp:62
v
Eigen::Matrix< real, 5, 1 > v
Definition test_ArrayDOF.cpp:468
Ax
tVec Ax(NCells)
x
tVec x(NCells)
b
tVec b(NCells)
alpha
real alpha
Definition test_GasThermo.cpp:454
dU
Eigen::Vector< real, 5 > dU
Definition test_GasThermo.cpp:452
ode
std::cout<< "[ESDIRK2-sc4] err1="<< std::scientific<< res.err_coarse<< " err2="<< res.err_fine<< " order="<< std::fixed<< res.order<< std::endl;CHECK(res.order > 1.8);CHECK(res.order< 2.5);}TEST_CASE("VBDF k=1: 1st-order on oscillator"){ ODE::ImplicitVBDFDualTimeStep< Vec2, Vec2 > ode(1, v2init, v2init, 1)
res
auto res
Definition test_ODE.cpp:196
n
Eigen::Vector3d n(1.0, 0.0, 0.0)