DNDSR 0.1.0.dev1+gcd065ad
Distributed Numeric Data Structure for CFV
Loading...
Searching...
No Matches
EulerEvaluator.hxx
Go to the documentation of this file.
1/** @file EulerEvaluator.hxx
2 * @brief Template implementations of EulerEvaluator methods for implicit time-stepping,
3 * LU-SGS preconditioning, DOF initialization, boundary value generation,
4 * reconstruction limiting, norm evaluation, and BC profile updates.
5 */
6#pragma once
7
8#include "DNDS/Defines.hpp" // for correct DNDS_SWITCH_INTELLISENSE
9#include "EulerEvaluator.hpp"
10#include "DNDS/HardEigen.hpp"
11#include "SpecialFields.hpp"
12#include "DNDS/ExprtkWrapper.hpp"
13
14namespace DNDS::Euler
15{
16
17 static const auto model = NS_SA; // to be hidden by template params
18
19 DNDS_SWITCH_INTELLISENSE(
20 template <EulerModel model>, )
21 /** @brief Assemble diagonal Jacobian block for each cell for LU-SGS preconditioning.
22 *
23 * Computes the spectral radius contribution from face eigenvalues and accumulates
24 * into JDiag. When settings.useRoeJacobian is enabled, also computes Roe flux
25 * Jacobian block contributions including boundary face linearization.
26 *
27 * @param JDiag Diagonal Jacobian block storage (output, cleared then filled).
28 * @param JSource Source-term Jacobian diagonal block (must match JDiag block mode).
29 * @param dTau Local pseudo-time step per cell.
30 * @param dt Global physical time step.
31 * @param alphaDiag Diagonal scaling factor for the spectral radius.
32 * @param u Conservative variable DOF array.
33 * @param uRec Reconstruction coefficients.
34 * @param jacobianCode Jacobian assembly mode (must be 0).
35 * @param t Current simulation time.
36 */
37 void EulerEvaluator<model>::LUSGSMatrixInit(
38 JacobianDiagBlock<nVarsFixed> &JDiag,
39 JacobianDiagBlock<nVarsFixed> &JSource,
40 ArrayDOFV<1> &dTau, real dt, real alphaDiag,
41 ArrayDOFV<nVarsFixed> &u,
42 ArrayRECV<nVarsFixed> &uRec,
43 int jacobianCode,
44 real t)
45 {
46 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
47 // TODO: for code0: flux jacobian with lambdaFace, and source jacobian with integration, only diagpart dealt with
48 DNDS_assert(JDiag.isBlock() == JSource.isBlock());
49 DNDS_assert(jacobianCode == 0);
50 if (settings.useRoeJacobian)
51 DNDS_assert(JDiag.isBlock());
52 JDiag.clearValues();
53 int cnvars = nVars;
54#if defined(DNDS_DIST_MT_USE_OMP)
55# pragma omp parallel for schedule(runtime)
56#endif
57 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
58 {
59 auto c2f = mesh->cell2face[iCell];
60
61 // LUSGS diag part
62 real fpDivisor = 1.0 / dTau[iCell](0) + 1.0 / dt;
63 for (int ic2f = 0; ic2f < c2f.size(); ic2f++)
64 {
65 index iFace = c2f[ic2f];
66 fpDivisor += (0.5 * alphaDiag) * vfv->GetFaceArea(iFace) * lambdaFace[iFace] / vfv->GetCellVol(iCell);
67 if (!settings.useRoeJacobian)
68 continue;
69 // roe term jacobi
70 auto f2c = mesh->face2cell[iFace];
71 index iCellOther = f2c[0] == iCell ? f2c[1] : f2c[0];
72 rowsize iCellAtFace = f2c[0] == iCell ? 0 : 1;
73 TVec unitNorm = vfv->GetFaceNormFromCell(iFace, iCell, iCellAtFace, -1)(Seq012) *
74 (iCellAtFace ? -1 : 1); // faces out
75 TU uj;
76 if (iCellOther == UnInitIndex) // handle BC
77 {
78 TU UL = u[iCell];
79 uj = this->generateBoundaryValue(UL, u[iCell], iCell, iFace, -1,
80 unitNorm,
81 Geom::NormBuildLocalBaseV<dim>(unitNorm),
82 vfv->GetFaceQuadraturePPhys(iFace, -1),
83 t,
84 mesh->GetFaceZone(iFace),
85 false, 0);
86 }
87 else
88 uj = u[iCellOther];
89 if (iCellOther != UnInitIndex)
90 this->UFromOtherCell(uj, iFace, iCell, iCellOther, iCellAtFace);
91 TJacobianU jacII = fluxJacobian0_Right_Times_du_AsMatrix( // unitnorm and uj are both respect with this cell
92 u[iCell], uj,
93 unitNorm, GetFaceVGridFromCell(iFace, iCell, iCellAtFace, -1),
94 mesh->GetFaceZone(iFace),
95 lambdaFace[iFace], lambdaFaceC[iFace], lambdaFaceVis[iFace],
96 iCellAtFace ? lambdaFace4[iFace] : lambdaFace0[iFace], lambdaFace123[iFace], iCellAtFace ? lambdaFace0[iFace] : lambdaFace4[iFace],
97 // swap lambda0 and lambda4 if iCellAtFace==1
98 true, +1, 1); // for this is diff(uthis) not diff(uthat)
99 if (iCellOther == UnInitIndex) // handle BC
100 {
101 TJacobianU JBC;
102 JBC.resize(nVars, nVars);
103 JBC.setIdentity();
104 for (int i = 0; i < nVars; i++)
105 {
106 TU VE = JBC(EigenAll, i);
107 JBC(EigenAll, i) = this->generateBoundaryValue(VE, u[iCell], iCell, iFace, -1,
108 unitNorm,
109 Geom::NormBuildLocalBaseV<dim>(unitNorm),
110 vfv->GetFaceQuadraturePPhys(iFace, -1),
111 t,
112 mesh->GetFaceZone(iFace),
113 false, 0, /*linMode=*/1);
114 }
115 TJacobianU jacIJ = fluxJacobian0_Right_Times_du_AsMatrix( // unitnorm and uj are both respect with this cell
116 uj, u[iCell],
117 unitNorm, GetFaceVGridFromCell(iFace, iCell, iCellAtFace, -1),
118 mesh->GetFaceZone(iFace),
119 lambdaFace[iFace], lambdaFaceC[iFace], lambdaFaceVis[iFace],
120 iCellAtFace ? lambdaFace4[iFace] : lambdaFace0[iFace], lambdaFace123[iFace], iCellAtFace ? lambdaFace0[iFace] : lambdaFace4[iFace],
121 // swap lambda0 and lambda4 if iCellAtFace==1
122 true, -1, 0);
123 JDiag.getBlock(iCell) += (0.5 * alphaDiag) * vfv->GetFaceArea(iFace) / vfv->GetCellVol(iCell) * (jacIJ * JBC);
124 }
125 JDiag.getBlock(iCell) += (0.5 * alphaDiag) * vfv->GetFaceArea(iFace) / vfv->GetCellVol(iCell) * jacII;
126 // std::cout << "JacII\n";
127 // std::cout << jacII << "\n";
128 // std::cout << lambdaFace[iFace] << std::endl;
129 }
130 if (settings.useRoeJacobian)
131 JDiag.getBlock(iCell).diagonal().array() += 1.0 / dTau[iCell](0) + 1.0 / dt; // time term!!
132 else
133 {
134 if (JDiag.isBlock())
135 JDiag.getBlock(iCell).diagonal().setConstant(fpDivisor);
136 else
137 JDiag.getDiag(iCell).setConstant(fpDivisor);
138 }
139
140 // std::cout << fpDivisor << std::endl;
141
142 // jacobian diag
143
144 if (!settings.ignoreSourceTerm)
145 {
146 if (JDiag.isBlock())
147 JDiag.getBlock(iCell) += alphaDiag * JSource.getBlock(iCell);
148 else
149 JDiag.getDiag(iCell) += alphaDiag * JSource.getDiag(iCell);
150 }
151
152 // jacobianCellInv[iCell] = jacobianCell[iCell].partialPivLu().inverse();
153
154 // std::cout << "jacobian Diag\n"
155 // << jacobianCell[iCell] << std::endl;
156 // std::cout << dTau[iCell] << "\n";
157 }
158 // exit(-1);
159 }
160
161 DNDS_SWITCH_INTELLISENSE(
162 template <EulerModel model>, )
163 /** @brief Compute the implicit matrix-vector product A*uInc for GMRES.
164 *
165 * Evaluates the action of the implicit operator (diagonal block + off-diagonal
166 * flux Jacobian contributions) on the increment vector, storing the result in AuInc.
167 *
168 * @param alphaDiag Diagonal scaling factor.
169 * @param t Current simulation time.
170 * @param u Conservative variable DOF array.
171 * @param uInc Increment vector to multiply.
172 * @param JDiag Diagonal Jacobian block.
173 * @param AuInc Result of the matrix-vector product (output).
174 */
175 void EulerEvaluator<model>::LUSGSMatrixVec(
176 real alphaDiag,
177 real t,
178 ArrayDOFV<nVarsFixed> &u,
179 ArrayDOFV<nVarsFixed> &uInc,
180 JacobianDiagBlock<nVarsFixed> &JDiag,
181 ArrayDOFV<nVarsFixed> &AuInc)
182 {
183 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
184 DNDS_MPI_InsertCheck(u.father->getMPI(), "LUSGSMatrixVec 1");
185 int cnvars = nVars;
186
187 auto cellOp = [&](index iCell) __attribute__((always_inline)) {
188
189 };
190
191#if defined(DNDS_DIST_MT_USE_OMP)
192# pragma omp parallel for schedule(static)
193#endif
194 for (int iPart = 0; iPart < mesh->NLocalParts(); iPart++)
195 for (index iScan = mesh->LocalPartStart(iPart); iScan < mesh->LocalPartEnd(iPart); iScan++)
196 {
197 index iCell = iScan;
198 cellOp(iCell);
199 // iCell = (*vfv->SOR_iScan2iCell)[iCell];//TODO: add rb-sor
200 auto c2f = mesh->cell2face[iCell];
201 TU uIncNewBuf(cnvars);
202 uIncNewBuf.setZero(); // norhs
203 auto uINCi = uInc[iCell];
204
205 if (uINCi.hasNaN())
206 {
207 std::cout << uINCi << std::endl;
208 DNDS_assert(false);
209 }
210
211 for (int ic2f = 0; ic2f < c2f.size(); ic2f++)
212 {
213 index iFace = c2f[ic2f];
214 auto f2c = mesh->face2cell[iFace];
215 index iCellOther = f2c[0] == iCell ? f2c[1] : f2c[0];
216 index iCellAtFace = f2c[0] == iCell ? 0 : 1;
217 TVec unitNorm = vfv->GetFaceNormFromCell(iFace, iCell, iCellAtFace, -1)(Seq012) *
218 (iCellAtFace ? -1 : 1); // faces out
219 if (iCellOther != UnInitIndex)
220 {
221
222 if (true)
223 {
224 TU uINCj = uInc[iCellOther];
225 TU uj = u[iCellOther];
226 this->UFromOtherCell(uINCj, iFace, iCell, iCellOther, iCellAtFace);
227 this->UFromOtherCell(uj, iFace, iCell, iCellOther, iCellAtFace);
228 TU fInc;
229 {
230
231 // fInc = fluxJacobian0_Right(
232 // u[iCellOther],
233 // unitNorm,
234 // BoundaryType::Inner) *
235 // uInc[iCellOther]; //! always inner here
236 fInc = fluxJacobian0_Right_Times_du(
237 uj, u[iCell],
238 unitNorm, GetFaceVGridFromCell(iFace, iCell, iCellAtFace, -1),
239 Geom::BC_ID_INTERNAL, uINCj,
240 lambdaFace[iFace], lambdaFaceC[iFace], lambdaFaceVis[iFace],
241 iCellAtFace ? lambdaFace4[iFace] : lambdaFace0[iFace], lambdaFace123[iFace], iCellAtFace ? lambdaFace0[iFace] : lambdaFace4[iFace],
242 // swap lambda0 and lambda4 if iCellAtFace==1
243 settings.useRoeJacobian); //! always inner here
244 }
245
246 uIncNewBuf -= (0.5 * alphaDiag) * vfv->GetFaceArea(iFace) / vfv->GetCellVol(iCell) *
247 (fInc);
248 if (uIncNewBuf.hasNaN() || (!uIncNewBuf.allFinite()))
249 {
250 std::cout
251 << fInc.transpose() << std::endl
252 << uInc[iCellOther].transpose() << std::endl;
253 DNDS_assert(!(uIncNewBuf.hasNaN() || (!uIncNewBuf.allFinite())));
254 }
255 }
256 }
257 }
258 // uIncNewBuf /= fpDivisor;
259 // uIncNew[iCell] = uIncNewBuf;
260 AuInc[iCell] = JDiag.MatVecLeft(iCell, uInc[iCell]) - uIncNewBuf;
261
262 auto AuIncI = AuInc[iCell];
263 if (AuIncI.hasNaN())
264 {
265 std::cout << AuIncI.transpose() << std::endl
266 << uINCi.transpose() << std::endl
267 << u[iCell].transpose() << std::endl
268 << JDiag.getValue(iCell) << std::endl
269 << iCell << std::endl;
270 DNDS_assert(!AuInc[iCell].hasNaN());
271 }
272 }
273 DNDS_MPI_InsertCheck(u.father->getMPI(), "LUSGSMatrixVec -1");
274 }
275
276 DNDS_SWITCH_INTELLISENSE(
277 template <EulerModel model>, )
278 /** @brief Assemble and factorize the local LU Jacobian from diagonal and off-diagonal blocks.
279 *
280 * Populates the sparse local LU structure with diagonal blocks from JDiag and
281 * off-diagonal Roe flux Jacobian entries for local-partition neighbor cells,
282 * then performs in-place LU decomposition.
283 *
284 * @param alphaDiag Diagonal scaling factor for face flux Jacobian.
285 * @param t Current simulation time.
286 * @param u Conservative variable DOF array.
287 * @param JDiag Diagonal Jacobian block.
288 * @param jacLU Local LU factorization structure (output).
289 */
290 void EulerEvaluator<model>::LUSGSMatrixToJacobianLU(
291 real alphaDiag, real t,
292 ArrayDOFV<nVarsFixed> &u,
293 JacobianDiagBlock<nVarsFixed> &JDiag,
294 JacobianLocalLU<nVarsFixed> &jacLU)
295 {
296 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
297 DNDS_MPI_InsertCheck(u.father->getMPI(), "LUSGSMatrixToJacobianLU 1");
298 int cnvars = nVars;
299 jacLU.setZero();
300#if defined(DNDS_DIST_MT_USE_OMP)
301# pragma omp parallel for schedule(static)
302#endif
303 for (int iPart = 0; iPart < mesh->NLocalParts(); iPart++)
304 for (index iScan = mesh->LocalPartStart(iPart); iScan < mesh->LocalPartEnd(iPart); iScan++)
305 {
306 index iCell = iScan;
307 auto c2f = mesh->cell2face[iCell];
308 for (int ic2f = 0; ic2f < c2f.size(); ic2f++)
309 {
310 index iFace = c2f[ic2f];
311 auto f2c = mesh->face2cell[iFace];
312 index iCellOther = f2c[0] == iCell ? f2c[1] : f2c[0];
313 rowsize iCellAtFace = f2c[0] == iCell ? 0 : 1;
314 TVec unitNorm = vfv->GetFaceNormFromCell(iFace, iCell, iCellAtFace, -1)(Seq012) *
315 (iCellAtFace ? -1 : 1); // faces out
316 if (iCellOther != UnInitIndex && iCellOther != iCell &&
317 iCellOther < mesh->LocalPartEnd(iPart) && iCellOther >= mesh->LocalPartStart(iPart))
318 {
319 TU uj = u[iCellOther];
320 this->UFromOtherCell(uj, iFace, iCell, iCellOther, iCellAtFace);
321 int iC2CInLocal = -1;
322 for (int ic2c = 0; ic2c < mesh->cell2cellFaceVLocalParts[iCell].size(); ic2c++)
323 if (iCellOther == mesh->cell2cellFaceVLocalParts[iCell][ic2c])
324 iC2CInLocal = ic2c; // TODO: pre-search this
325 DNDS_assert(iC2CInLocal != -1);
326 TJacobianU jacIJ;
327 {
328 jacIJ = fluxJacobian0_Right_Times_du_AsMatrix( // unitnorm and uj are both respect with this cell
329 uj, u[iCell],
330 unitNorm, GetFaceVGridFromCell(iFace, iCell, iCellAtFace, -1),
331 Geom::BC_ID_INTERNAL,
332 lambdaFace[iFace], lambdaFaceC[iFace], lambdaFaceVis[iFace],
333 iCellAtFace ? lambdaFace4[iFace] : lambdaFace0[iFace], lambdaFace123[iFace], iCellAtFace ? lambdaFace0[iFace] : lambdaFace4[iFace],
334 // swap lambda0 and lambda4 if iCellAtFace==1
335 settings.useRoeJacobian); //! always inner here
336 }
337 auto faceID = mesh->GetFaceZone(iFace);
338 mesh->CellOtherCellPeriodicHandle(
339 iFace, iCellAtFace,
340 [&]()
341 { jacIJ(EigenAll, Seq123) =
342 mesh->periodicInfo.TransVectorBack<dim, nVarsFixed>(
343 jacIJ(EigenAll, Seq123).transpose(), faceID)
344 .transpose(); },
345 [&]()
346 { jacIJ(EigenAll, Seq123) =
347 mesh->periodicInfo.TransVector<dim, nVarsFixed>(
348 jacIJ(EigenAll, Seq123).transpose(), faceID)
349 .transpose(); });
350 jacLU.LDU(iCell, symLU->cell2cellFaceVLocal2FullRowPos[iCell][iC2CInLocal]) =
351 (0.5 * alphaDiag) * vfv->GetFaceArea(iFace) / vfv->GetCellVol(iCell) * jacIJ;
352 }
353 }
354 jacLU.GetDiag(iCell) = JDiag.getValue(iCell);
355 }
356 // TODO: make below OMP-ed
357 jacLU.InPlaceDecompose();
358 DNDS_MPI_InsertCheck(u.father->getMPI(), "LUSGSMatrixToJacobianLU -1");
359 }
360
361 DNDS_SWITCH_INTELLISENSE(
362 template <EulerModel model>, template <>)
363 /** @brief [[deprecated]] Forward LU-SGS sweep. Use UpdateSGS with uIncIsZero=true instead.
364 *
365 * Performs a forward Gauss-Seidel sweep using diagonal inversion and lower-triangular
366 * off-diagonal flux Jacobian contributions. Retained for callers that alias uInc==uIncNew.
367 *
368 * @param alphaDiag Diagonal scaling factor.
369 * @param t Current simulation time.
370 * @param rhs Right-hand side residual.
371 * @param u Conservative variable DOF array.
372 * @param uInc Current increment (input, may alias uIncNew).
373 * @param JDiag Diagonal Jacobian block.
374 * @param uIncNew Updated increment (output, may alias uInc).
375 */
376 void EulerEvaluator<model>::UpdateLUSGSForward(
377 real alphaDiag, real t,
378 ArrayDOFV<nVarsFixed> &rhs,
379 ArrayDOFV<nVarsFixed> &u,
380 ArrayDOFV<nVarsFixed> &uInc,
381 JacobianDiagBlock<nVarsFixed> &JDiag,
382 ArrayDOFV<nVarsFixed> &uIncNew)
383 {
384 // Deprecated: delegate to UpdateSGS with uIncIsZero.
385 // Note: callers pass uInc == uIncNew (aliased), but UpdateSGS requires
386 // uInc != uIncNew. We keep the original implementation for now to
387 // preserve this aliasing behavior. Migrate callers to use UpdateSGS
388 // with separate buffers instead.
389 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
390 DNDS_MPI_InsertCheck(u.father->getMPI(), "UpdateLUSGSForward 1");
391 int cnvars = nVars;
392 JDiag.GetInvert();
393 index nCellDist = mesh->NumCell();
394 for (index iScan = 0; iScan < nCellDist; iScan++)
395 {
396 index iCell = iScan;
397
398 auto c2f = mesh->cell2face[iCell];
399 TU uIncNewBuf(nVars);
400 auto RHSI = rhs[iCell];
401 uIncNewBuf = RHSI;
402
403 for (int ic2f = 0; ic2f < c2f.size(); ic2f++)
404 {
405 index iFace = c2f[ic2f];
406 auto f2c = mesh->face2cell[iFace];
407 index iCellOther = f2c[0] == iCell ? f2c[1] : f2c[0];
408 index iCellAtFace = f2c[0] == iCell ? 0 : 1;
409 TVec unitNorm = vfv->GetFaceNormFromCell(iFace, iCell, iCellAtFace, -1)(Seq012) *
410 (iCellAtFace ? -1 : 1); // faces out
411 if (iCellOther != UnInitIndex)
412 {
413 index iScanOther = iCellOther;
414 if (iScanOther < iScan)
415 {
416 TU uINCj = uInc[iCellOther];
417 TU uj = u[iCellOther];
418 this->UFromOtherCell(uINCj, iFace, iCell, iCellOther, iCellAtFace);
419 this->UFromOtherCell(uj, iFace, iCell, iCellOther, iCellAtFace);
420
421 TU fInc = fluxJacobian0_Right_Times_du(
422 uj, u[iCell],
423 unitNorm, GetFaceVGridFromCell(iFace, iCell, iCellAtFace, -1),
424 Geom::BC_ID_INTERNAL, uINCj,
425 lambdaFace[iFace], lambdaFaceC[iFace], lambdaFaceVis[iFace],
426 iCellAtFace ? lambdaFace4[iFace] : lambdaFace0[iFace], lambdaFace123[iFace], iCellAtFace ? lambdaFace0[iFace] : lambdaFace4[iFace],
427 settings.useRoeJacobian);
428
429 uIncNewBuf -= (0.5 * alphaDiag) * vfv->GetFaceArea(iFace) / vfv->GetCellVol(iCell) *
430 (fInc);
431
432 if ((!uIncNewBuf.allFinite()))
433 {
434 std::cout << RHSI.transpose() << std::endl
435 << fInc.transpose() << std::endl
436 << uINCj.transpose() << std::endl;
437 DNDS_assert(false);
438 }
439 }
440 }
441 }
442 auto uIncNewI = uIncNew[iCell];
443 uIncNewI = JDiag.MatVecLeftInvert(iCell, uIncNewBuf);
444
445 if (uIncNewI.hasNaN())
446 {
447 std::cout << uIncNewI.transpose() << std::endl
448 << uIncNewBuf.transpose() << std::endl
449 << JDiag.getValue(iCell) << std::endl
450 << iCell << std::endl;
451 DNDS_assert(!uIncNew[iCell].hasNaN());
452 }
453 }
454 DNDS_MPI_InsertCheck(u.father->getMPI(), "UpdateLUSGSForward -1");
455 }
456
457 DNDS_SWITCH_INTELLISENSE(
458 template <EulerModel model>, template <>)
459 /** @brief [[deprecated]] Backward LU-SGS sweep. Use UpdateSGS instead.
460 *
461 * Performs a backward Gauss-Seidel sweep using diagonal inversion and upper-triangular
462 * off-diagonal flux Jacobian contributions. Adds the correction to uIncNew in-place.
463 *
464 * @param alphaDiag Diagonal scaling factor.
465 * @param t Current simulation time.
466 * @param rhs Right-hand side residual.
467 * @param u Conservative variable DOF array.
468 * @param uInc Current increment (input).
469 * @param JDiag Diagonal Jacobian block.
470 * @param uIncNew Updated increment (input/output, correction added).
471 */
472 void EulerEvaluator<model>::UpdateLUSGSBackward(
473 real alphaDiag, real t,
474 ArrayDOFV<nVarsFixed> &rhs,
475 ArrayDOFV<nVarsFixed> &u,
476 ArrayDOFV<nVarsFixed> &uInc,
477 JacobianDiagBlock<nVarsFixed> &JDiag,
478 ArrayDOFV<nVarsFixed> &uIncNew)
479 {
480 // Deprecated: see UpdateLUSGSForward comment.
481 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
482 DNDS_MPI_InsertCheck(u.father->getMPI(), "UpdateLUSGSBackward 1");
483 int cnvars = nVars;
484 JDiag.GetInvert();
485 index nCellDist = mesh->NumCell();
486 for (index iScan = nCellDist - 1; iScan >= 0; iScan--)
487 {
488 index iCell = iScan;
489
490 auto c2f = mesh->cell2face[iCell];
491 TU uIncNewBuf(cnvars);
492 uIncNewBuf.setZero();
493
494 for (int ic2f = 0; ic2f < c2f.size(); ic2f++)
495 {
496 index iFace = c2f[ic2f];
497 auto f2c = mesh->face2cell[iFace];
498 index iCellOther = f2c[0] == iCell ? f2c[1] : f2c[0];
499 index iCellAtFace = f2c[0] == iCell ? 0 : 1;
500 TVec unitNorm = vfv->GetFaceNormFromCell(iFace, iCell, iCellAtFace, -1)(Seq012) *
501 (iCellAtFace ? -1 : 1); // faces out
502 if (iCellOther != UnInitIndex)
503 {
504 index iScanOther = iCellOther;
505 if (iScanOther > iScan)
506 {
507 TU uINCj = uInc[iCellOther];
508 TU uj = u[iCellOther];
509 this->UFromOtherCell(uINCj, iFace, iCell, iCellOther, iCellAtFace);
510 this->UFromOtherCell(uj, iFace, iCell, iCellOther, iCellAtFace);
511
512 TU fInc = fluxJacobian0_Right_Times_du(
513 uj, u[iCell],
514 unitNorm, GetFaceVGridFromCell(iFace, iCell, iCellAtFace, -1),
515 Geom::BC_ID_INTERNAL, uINCj,
516 lambdaFace[iFace], lambdaFaceC[iFace], lambdaFaceVis[iFace],
517 iCellAtFace ? lambdaFace4[iFace] : lambdaFace0[iFace], lambdaFace123[iFace], iCellAtFace ? lambdaFace0[iFace] : lambdaFace4[iFace],
518 settings.useRoeJacobian);
519
520 uIncNewBuf -= (0.5 * alphaDiag) * vfv->GetFaceArea(iFace) / vfv->GetCellVol(iCell) *
521 (fInc);
522 }
523 }
524 }
525 auto uIncNewI = uIncNew[iCell];
526 uIncNewI += JDiag.MatVecLeftInvert(iCell, uIncNewBuf);
527 }
528 DNDS_MPI_InsertCheck(u.father->getMPI(), "UpdateLUSGSBackward -1");
529 }
530
531 template <EulerModel model>
532 /** @brief Symmetric Gauss-Seidel (SGS) sweep for implicit preconditioning.
533 *
534 * Performs a single forward or backward sweep. When uIncIsZero is true, skips
535 * off-diagonal terms involving the zeroed increment in the first forward pass
536 * to save computation. Accumulates the global L1 increment change into sumInc.
537 *
538 * @param alphaDiag Diagonal scaling factor.
539 * @param t Current simulation time.
540 * @param rhs Right-hand side residual.
541 * @param u Conservative variable DOF array.
542 * @param uInc Current increment (input, must not alias uIncNew).
543 * @param uIncNew Updated increment (output).
544 * @param JDiag Diagonal Jacobian block.
545 * @param forward If true, sweep from cell 0 to N-1; otherwise N-1 to 0.
546 * @param gsUpdate If true, use Gauss-Seidel update (use latest values); otherwise Jacobi.
547 * @param sumInc Global L1 sum of increment changes (output, MPI-reduced).
548 * @param uIncIsZero If true, skip off-diagonal terms for zero-initialized increment.
549 */
550 void EulerEvaluator<model>::UpdateSGS(
551 real alphaDiag, real t,
552 ArrayDOFV<nVarsFixed> &rhs,
553 ArrayDOFV<nVarsFixed> &u,
554 ArrayDOFV<nVarsFixed> &uInc,
555 ArrayDOFV<nVarsFixed> &uIncNew,
556 JacobianDiagBlock<nVarsFixed> &JDiag,
557 bool forward, bool gsUpdate, TU &sumInc,
558 bool uIncIsZero)
559 {
560 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
561 DNDS_MPI_InsertCheck(u.father->getMPI(), "UpdateSGS 1");
562 DNDS_assert(&uInc != &uIncNew);
563 int cnvars = nVars;
564 JDiag.GetInvert();
565 const index nCellDist = mesh->NumCell();
566 sumInc.setZero(cnvars);
567
568#if defined(DNDS_DIST_MT_USE_OMP)
569 auto cellOp = [&](index iScanStart, index iScanEnd)
570 {
571#else
572 const index iScanStart{0}, iScanEnd(nCellDist);
573#endif
574 for (index iScan = iScanStart; iScan < iScanEnd; iScan++)
575 {
576 index iCell = forward ? iScan : iScanEnd - 1 - (iScan - iScanStart);
577
578 auto c2f = mesh->cell2face[iCell];
579 TU uIncNewBuf(nVars);
580 auto RHSI = rhs[iCell];
581 // std::cout << rhs[iCell](0) << std::endl;
582 uIncNewBuf = RHSI;
583
584 for (int ic2f = 0; ic2f < c2f.size(); ic2f++)
585 {
586 index iFace = c2f[ic2f];
587 auto btype = mesh->GetFaceZone(iFace);
588 auto f2c = mesh->face2cell[iFace];
589 index iCellOther = f2c[0] == iCell ? f2c[1] : f2c[0];
590 index iCellAtFace = f2c[0] == iCell ? 0 : 1;
591 TVec unitNorm = vfv->GetFaceNormFromCell(iFace, iCell, iCellAtFace, -1)(Seq012) *
592 (iCellAtFace ? -1 : 1); // faces out
593 if (iCellOther != UnInitIndex)
594 {
595 if (iCell != iCellOther)
596 {
597 TU fInc;
598 bool iCellOtherIsThisPart = iScanStart <= iCellOther && iCellOther < iScanEnd;
599 bool gsUseNew = gsUpdate && iCellOtherIsThisPart && (forward ? (iCellOther < iCell) : (iCellOther > iCell));
600
601 // When uInc is known to be zero, skip flux for not-yet-processed
602 // neighbours whose increment is still zero (LUSGS optimisation).
603 if (uIncIsZero && !gsUseNew)
604 continue;
605
606 TU uINCj = gsUseNew ? uIncNew[iCellOther] : uInc[iCellOther];
607 TU uj = u[iCellOther];
608 this->UFromOtherCell(uINCj, iFace, iCell, iCellOther, iCellAtFace);
609 this->UFromOtherCell(uj, iFace, iCell, iCellOther, iCellAtFace);
610
611 {
612 fInc = fluxJacobian0_Right_Times_du(
613 uj, u[iCell],
614 unitNorm, GetFaceVGridFromCell(iFace, iCell, iCellAtFace, -1),
615 Geom::BC_ID_INTERNAL, uINCj,
616 lambdaFace[iFace], lambdaFaceC[iFace], lambdaFaceVis[iFace],
617 iCellAtFace ? lambdaFace4[iFace] : lambdaFace0[iFace], lambdaFace123[iFace], iCellAtFace ? lambdaFace0[iFace] : lambdaFace4[iFace],
618 // swap lambda0 and lambda4 if iCellAtFace==1
619 settings.useRoeJacobian); //! always inner here
620 }
621
622 uIncNewBuf -= (0.5 * alphaDiag) * vfv->GetFaceArea(iFace) / vfv->GetCellVol(iCell) *
623 (fInc);
624
625 if ((!uIncNewBuf.allFinite()))
626 {
627 std::cout << RHSI.transpose() << std::endl
628 << fInc.transpose() << std::endl
629 << uINCj.transpose() << std::endl;
630 DNDS_assert(false);
631 }
632 }
633 }
634 // else if (pBCHandler->GetTypeFromID(btype) == BCWall || pBCHandler->GetTypeFromID(btype) == BCWallIsothermal)
635 // {
636 // TMat normBase = Geom::NormBuildLocalBaseV<dim>(unitNorm);
637 // Geom::tPoint pPhysics = vfv->GetFaceQuadraturePPhysFromCell(iFace, iCell, iCellAtFace, -1);
638 // TU uThis = u[iCell];
639 // TU uINCj = uInc[iCell];
640 // //! using t = 0 in generateBoudnaryValue!
641 // TU uj = generateBoundaryValue(uThis, uThis, iCell, iFace, -1, unitNorm, normBase, pPhysics, 0, btype, false, 0);
642 // uINCj(Seq123) *= -1;
643
644 // if (model == NS_SA || model == NS_SA_3D)
645 // uINCj(I4 + 1) = 0;
646 // if (model == NS_2EQ || model == NS_2EQ_3D)
647 // uINCj({I4 + 1, I4 + 2}).setZero();
648 // TU fInc = fluxJacobian0_Right_Times_du(
649 // uj, u[iCell],
650 // unitNorm, GetFaceVGridFromCell(iFace, iCell, iCellAtFace, -1),
651 // Geom::BC_ID_INTERNAL, uINCj,
652 // lambdaFace[iFace], lambdaFaceC[iFace], lambdaFaceVis[iFace],
653 // iCellAtFace ? lambdaFace4[iFace] : lambdaFace0[iFace], lambdaFace123[iFace], iCellAtFace ? lambdaFace0[iFace] : lambdaFace4[iFace],
654 // // swap lambda0 and lambda4 if iCellAtFace==1
655 // settings.useRoeJacobian); //! treat as inner here
656
657 // uIncNewBuf -= (0.5 * alphaDiag) * vfv->GetFaceArea(iFace) / vfv->GetCellVol(iCell) *
658 // (fInc);
659 // }
660 }
661 auto uIncNewI = uIncNew[iCell];
662 TU uIncOld = uIncNewI;
663
664 uIncNewI = JDiag.MatVecLeftInvert(iCell, uIncNewBuf);
665 sumInc.array() += (uIncNewI - uIncOld).array().abs();
666
667 if (uIncNewI.hasNaN())
668 {
669 std::cout << uIncNewI.transpose() << std::endl
670 << uIncNewBuf.transpose() << std::endl
671 << JDiag.getValue(iCell) << std::endl
672 << iCell << std::endl;
673 DNDS_assert(!uInc[iCell].hasNaN());
674 }
675 }
676#if defined(DNDS_DIST_MT_USE_OMP)
677 };
678#endif
679
680#if defined(DNDS_DIST_MT_USE_OMP)
681# pragma omp declare reduction(TUAdd:TU : omp_out += omp_in) initializer(omp_priv = omp_orig)
682# pragma omp parallel for schedule(static) reduction(TUAdd : sumInc)
683 for (int iPart = 0; iPart < mesh->NLocalParts(); iPart++)
684 cellOp(mesh->LocalPartStart(iPart), mesh->LocalPartEnd(iPart));
685#endif
686
687 TU sumIncAll(cnvars);
688 // std::abort();
689 MPI::Allreduce(sumInc.data(), sumIncAll.data(), sumInc.size(), DNDS_MPI_REAL, MPI_SUM, rhs.father->getMPI().comm);
690 sumInc = sumIncAll;
691 DNDS_MPI_InsertCheck(u.father->getMPI(), "UpdateSGS -1");
692 // exit(-1);
693 }
694
695 template <EulerModel model>
696 /** @brief SGS sweep with reconstruction-level off-diagonal terms.
697 *
698 * Extends UpdateSGS by including reconstruction increment contributions
699 * (uRecInc) in the off-diagonal flux Jacobian product, enabling higher-order
700 * implicit coupling.
701 *
702 * @param alphaDiag Diagonal scaling factor.
703 * @param t Current simulation time.
704 * @param rhs Right-hand side residual.
705 * @param u Conservative variable DOF array.
706 * @param uRec Reconstruction coefficients.
707 * @param uInc Current DOF increment.
708 * @param uRecInc Current reconstruction increment.
709 * @param JDiag Diagonal Jacobian block.
710 * @param forward Sweep direction: true for forward, false for backward.
711 * @param sumInc Global L1 sum of increment changes (output, MPI-reduced).
712 */
713 void EulerEvaluator<model>::UpdateSGSWithRec(
714 real alphaDiag, real t,
715 ArrayDOFV<nVarsFixed> &rhs,
716 ArrayDOFV<nVarsFixed> &u,
717 ArrayRECV<nVarsFixed> &uRec,
718 ArrayDOFV<nVarsFixed> &uInc,
719 ArrayRECV<nVarsFixed> &uRecInc,
720 JacobianDiagBlock<nVarsFixed> &JDiag,
721 bool forward, TU &sumInc)
722 {
723 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
724 DNDS_MPI_InsertCheck(u.father->getMPI(), "UpdateSGS 1");
725 int cnvars = nVars;
726 JDiag.GetInvert();
727 index nCellDist = mesh->NumCell();
728 sumInc.setZero(cnvars);
729#if defined(DNDS_DIST_MT_USE_OMP)
730# pragma omp declare reduction(TUAdd:TU : omp_out += omp_in) initializer(omp_priv = omp_orig)
731# pragma omp parallel for schedule(static) reduction(TUAdd : sumInc)
732#endif
733 for (index iScan = 0; iScan < nCellDist; iScan++)
734 {
735 index iCell = iScan;
736 iCell = forward ? iScan : nCellDist - 1 - iScan; // TODO: add rb-sor
737
738 auto c2f = mesh->cell2face[iCell];
739 TU uIncNewBuf(nVars);
740 auto RHSI = rhs[iCell];
741 // std::cout << rhs[iCell](0) << std::endl;
742 uIncNewBuf = RHSI;
743
744 for (int ic2f = 0; ic2f < c2f.size(); ic2f++)
745 {
746 index iFace = c2f[ic2f];
747 auto f2c = mesh->face2cell[iFace];
748 index iCellOther = f2c[0] == iCell ? f2c[1] : f2c[0];
749 index iCellAtFace = f2c[0] == iCell ? 0 : 1;
750 TVec unitNorm = vfv->GetFaceNormFromCell(iFace, iCell, iCellAtFace, -1)(Seq012) *
751 (iCellAtFace ? -1 : 1); // faces out
752 if (iCellOther != UnInitIndex)
753 {
754 index iScanOther = forward ? iCellOther : nCellDist - 1 - iCellOther; // TODO: add rb-sor
755 if (iCell != iCellOther)
756 {
757 TU fInc, fIncS;
758 auto uINCj = uInc[iCellOther];
759 {
760 fInc = fluxJacobian0_Right_Times_du(
761 u[iCellOther], u[iCell], //! TODO periodic here
762 unitNorm, GetFaceVGridFromCell(iFace, iCell, iCellAtFace, -1),
763 Geom::BC_ID_INTERNAL, uINCj,
764 lambdaFace[iFace], lambdaFaceC[iFace], lambdaFaceVis[iFace],
765 iCellAtFace ? lambdaFace4[iFace] : lambdaFace0[iFace], lambdaFace123[iFace], iCellAtFace ? lambdaFace0[iFace] : lambdaFace4[iFace],
766 // swap lambda0 and lambda4 if iCellAtFace==1
767 settings.useRoeJacobian); //! always inner here
768 }
769 {
770 TU uRecSLInc =
771 (vfv->GetIntPointDiffBaseValue(iCell, iFace, iCellAtFace, -1, std::array<int, 1>{0}, 1) *
772 uRecInc[iCell])
773 .transpose();
774 TU uRecSRInc =
775 (vfv->GetIntPointDiffBaseValue(iCellOther, iFace, 1 - iCellAtFace, -1, std::array<int, 1>{0}, 1) *
776 uRecInc[iCellOther]) //! TODO periodic here
777 .transpose();
778 TU fIncSL = fluxJacobianC_Right_Times_du(u[iCell], unitNorm, GetFaceVGridFromCell(iFace, iCell, iCellAtFace, -1), Geom::BC_ID_INTERNAL, uRecSLInc);
779 TU fIncSR = fluxJacobianC_Right_Times_du(u[iCellOther], unitNorm, GetFaceVGridFromCell(iFace, iCell, iCellAtFace, -1), Geom::BC_ID_INTERNAL, uRecSRInc);
780 fIncS = fIncSL + fIncSR + lambdaFaceC[iFace] * (uRecSLInc - uRecSRInc);
781 }
782
783 uIncNewBuf -= (0.5 * alphaDiag) * vfv->GetFaceArea(iFace) / vfv->GetCellVol(iCell) *
784 (fInc);
785 uIncNewBuf -= (0.5 * alphaDiag) * vfv->GetFaceArea(iFace) / vfv->GetCellVol(iCell) *
786 (fIncS);
787
788 if ((!uIncNewBuf.allFinite()))
789 {
790 std::cout << RHSI.transpose() << std::endl
791 << fInc.transpose() << std::endl
792 << uINCj.transpose() << std::endl;
793 DNDS_assert(false);
794 }
795 }
796 }
797 }
798 auto uIncNewI = uInc[iCell];
799 TU uIncOld = uIncNewI;
800
801 uIncNewI = JDiag.MatVecLeftInvert(iCell, uIncNewBuf);
802 sumInc.array() += (uIncNewI - uIncOld).array().abs();
803
804 if (uIncNewI.hasNaN())
805 {
806 std::cout << uIncNewI.transpose() << std::endl
807 << uIncNewBuf.transpose() << std::endl
808 << JDiag.getValue(iCell) << std::endl
809 << iCell << std::endl;
810 DNDS_assert(!uInc[iCell].hasNaN());
811 }
812 }
813 TU sumIncAll(cnvars);
814 MPI::Allreduce(sumInc.data(), sumIncAll.data(), sumInc.size(), DNDS_MPI_REAL, MPI_SUM, rhs.father->getMPI().comm);
815 DNDS_MPI_InsertCheck(u.father->getMPI(), "UpdateSGS -1");
816 }
817
818 DNDS_SWITCH_INTELLISENSE(
819 template <EulerModel model>, template <>)
820 /** @brief Solve the local LU-factorized Jacobian system for the increment.
821 *
822 * Applies the pre-factorized local LU solve (from LUSGSMatrixToJacobianLU) to
823 * compute the new increment from the RHS residual and off-diagonal contributions.
824 * Supports accumulation onto a nonzero initial increment when uIncIsZero is false.
825 *
826 * @param alphaDiag Diagonal scaling factor.
827 * @param t Current simulation time.
828 * @param rhs Right-hand side residual.
829 * @param u Conservative variable DOF array.
830 * @param uInc Current increment (input).
831 * @param uIncNew Updated increment (output, must not alias uInc).
832 * @param bBuf Temporary buffer for RHS assembly.
833 * @param JDiag Diagonal Jacobian block.
834 * @param jacLU Pre-factorized local LU structure.
835 * @param uIncIsZero If true, uInc is zero and can be skipped in accumulation.
836 * @param sumInc Component-wise L1 norm of the new increment (output).
837 */
838 void EulerEvaluator<model>::LUSGSMatrixSolveJacobianLU(
839 real alphaDiag, real t,
840 ArrayDOFV<nVarsFixed> &rhs,
841 ArrayDOFV<nVarsFixed> &u,
842 ArrayDOFV<nVarsFixed> &uInc,
843 ArrayDOFV<nVarsFixed> &uIncNew,
844 ArrayDOFV<nVarsFixed> &bBuf,
845 JacobianDiagBlock<nVarsFixed> &JDiag,
846 JacobianLocalLU<nVarsFixed> &jacLU,
847 bool uIncIsZero,
848 TU &sumInc)
849 {
850 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
851 DNDS_MPI_InsertCheck(u.father->getMPI(), "LUSGSMatrixSolveJacobianLU 1");
852 int cnvars = nVars;
853 index nCellDist = mesh->NumCell();
854 sumInc.setZero(cnvars);
855 //- revert to old LU solve here
856 if (!uIncIsZero)
857 this->LUSGSMatrixVec(alphaDiag, t, u, uInc, JDiag, bBuf); // bBuf = Ax
858 else
859 bBuf.setConstant(0.0);
860#if defined(DNDS_DIST_MT_USE_OMP)
861# pragma omp parallel for schedule(static)
862#endif
863 for (int iPart = 0; iPart < mesh->NLocalParts(); iPart++)
864 for (index iScan = mesh->LocalPartStart(iPart); iScan < mesh->LocalPartEnd(iPart); iScan++)
865 // update the ghost part (non proc-block) rhs
866 {
867 index iCell = iScan;
868 auto c2f = mesh->cell2face[iCell];
869 // std::cout << uIncIsZero << " " << uInc[iCell].norm() << " " << bBuf[iCell] << std::endl;
870 bBuf[iCell] = rhs[iCell] - bBuf[iCell]; // bBuf = b - Ax
871
872 // bBuf[iCell] = rhs[iCell]; //+ revert to old LU solve here
873 // if (uIncIsZero) //+ revert to old LU solve here
874 continue;
875
876 for (int ic2f = 0; ic2f < c2f.size(); ic2f++)
877 {
878 index iFace = c2f[ic2f];
879 auto f2c = mesh->face2cell[iFace];
880 index iCellOther = f2c[0] == iCell ? f2c[1] : f2c[0];
881 index iCellAtFace = f2c[0] == iCell ? 0 : 1;
882 TVec unitNorm = vfv->GetFaceNormFromCell(iFace, iCell, iCellAtFace, -1)(Seq012) *
883 (iCellAtFace ? -1 : 1); // faces out
884 if (iCellOther != UnInitIndex && iCell != iCellOther
885 // if is a ghost neighbour
886 && iCellOther >= mesh->NumCell())
887 {
888 TU fInc;
889 TU uINCj = uInc[iCellOther];
890 TU uj = u[iCellOther];
891 this->UFromOtherCell(uINCj, iFace, iCell, iCellOther, iCellAtFace);
892 this->UFromOtherCell(uj, iFace, iCell, iCellOther, iCellAtFace);
893 {
894
895 fInc = fluxJacobian0_Right_Times_du(
896 uj, u[iCell],
897 unitNorm, GetFaceVGridFromCell(iFace, iCell, iCellAtFace, -1),
898 Geom::BC_ID_INTERNAL, uINCj,
899 lambdaFace[iFace], lambdaFaceC[iFace], lambdaFaceVis[iFace],
900 iCellAtFace ? lambdaFace4[iFace] : lambdaFace0[iFace], lambdaFace123[iFace], iCellAtFace ? lambdaFace0[iFace] : lambdaFace4[iFace],
901 // swap lambda0 and lambda4 if iCellAtFace==1
902 settings.useRoeJacobian); //! always inner here
903 }
904
905 bBuf[iCell] -= (0.5 * alphaDiag) * vfv->GetFaceArea(iFace) / vfv->GetCellVol(iCell) *
906 (fInc);
907 }
908 }
909 // TU uIncOld = uIncNew[iCell];
910 // uIncNew[iCell] = JDiag[iCell].array().inverse() * bBuf[iCell].array();
911 // sumInc.array() += (uIncNew[iCell] - uIncOld).array().abs();
912 }
913 jacLU.Solve(bBuf, uIncNew); // top-diagonal solve
914 sumInc = uIncNew.componentWiseNorm1();
915
916 //- revert to old LU solve here
917 if (!uIncIsZero)
918 uIncNew += uInc;
919
920 DNDS_assert(uIncNew.father.get() != uInc.father.get()); // no aliasing
921 // TU sumIncAll(cnvars);
922 // MPI::Allreduce(sumInc.data(), sumIncAll.data(), sumInc.size(), DNDS_MPI_REAL, MPI_SUM, rhs.father->getMPI().comm);
923 // sumInc = sumIncAll;
924
925 DNDS_MPI_InsertCheck(u.father->getMPI(), "LUSGSMatrixSolveJacobianLU -1");
926 // exit(-1);
927 }
928
929 DNDS_SWITCH_INTELLISENSE(
930 template <EulerModel model>, )
931 /** @brief Initialize conservative variable DOF from farfield and special-field initializers.
932 *
933 * Sets all cells to the farfield static value, then applies model-specific
934 * initialization (SA wall-distance scaling), box/plane/exprtk region overrides,
935 * and built-in special fields (isentropic vortex, Rayleigh-Taylor, Taylor-Green,
936 * double Mach reflection, etc.) based on settings.
937 *
938 * @param u Conservative variable DOF array (output, fully initialized).
939 */
940 void EulerEvaluator<model>::InitializeUDOF(ArrayDOFV<nVarsFixed> &u)
941 {
942 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
943 Eigen::VectorXd initConstVal = this->settings.farFieldStaticValue;
944 u.setConstant(initConstVal);
945 if (model == EulerModel::NS_SA || model == NS_SA_3D)
946 {
947#if defined(DNDS_DIST_MT_USE_OMP)
948# pragma omp parallel for schedule(runtime)
949#endif
950 for (int iCell = 0; iCell < mesh->NumCell(); iCell++)
951 {
952 auto c2f = mesh->cell2face[iCell];
953 for (int ic2f = 0; ic2f < c2f.size(); ic2f++)
954 {
955 index iFace = c2f[ic2f];
956 if (pBCHandler->GetTypeFromID(mesh->GetFaceZone(iFace)) == EulerBCType::BCWall ||
957 pBCHandler->GetTypeFromID(mesh->GetFaceZone(iFace)) == EulerBCType::BCWallIsothermal)
958 u[iCell](I4 + 1) *= 1.0; // ! not fixing first layer!
959 }
960 }
961 }
962 if (model == EulerModel::NS_2EQ || model == NS_2EQ_3D)
963 {
964 if (settings.ransModel == RANSModel::RANS_KOSST ||
965 settings.ransModel == RANSModel::RANS_KOWilcox)
966#if defined(DNDS_DIST_MT_USE_OMP)
967# pragma omp parallel for schedule(runtime)
968#endif
969 for (int iCell = 0; iCell < mesh->NumCell(); iCell++)
970 {
971 auto c2f = mesh->cell2face[iCell];
972 real d = dWall.at(iCell).mean();
973 // for SST or KOWilcox
974 for (int ic2f = 0; ic2f < c2f.size(); ic2f++)
975 // if (pBCHandler->GetTypeFromID(mesh->GetFaceZone(c2f[ic2f])) == EulerBCType::BCWall ||
976 // pBCHandler->GetTypeFromID(mesh->GetFaceZone(c2f[ic2f])) == EulerBCType::BCWallIsothermal)
977 {
978 real pMean, asqrMean, Hmean;
979 real gamma = settings.idealGasProperty.gamma;
980 Gas::IdealGasThermal(u[iCell](I4), u[iCell](0), (u[iCell](Seq123) / u[iCell](0)).squaredNorm(),
981 gamma, pMean, asqrMean, Hmean);
982 real muRef = settings.idealGasProperty.muGas;
983 real T = pMean / ((gamma - 1) / gamma * settings.idealGasProperty.CpGas * u[iCell](0));
984 real mufPhy1 = muEff(u[iCell], T);
985 real rhoOmegaaaWall = mufPhy1 / sqr(d) * RANS::kWallOmegaCoeff * 0.1;
986
987 real rhoOmegaaaNew = std::max(rhoOmegaaaWall, u[iCell](I4 + 2));
988 real rhoOmegaaaOld = u[iCell](I4 + 2);
989 // u[iCell](I4 + 2) = rhoOmegaaaNew;
990 // u[iCell](I4 + 1) = rhoOmegaaaNew / rhoOmegaaaOld * u[iCell](I4 + 1);
991 }
992 }
993 }
994
995 switch (settings.specialBuiltinInitializer)
996 {
997 case 1: // for RT problem
998 DNDS_assert(model == NS || model == NS_2D || model == NS_3D);
999 if constexpr (model == NS || model == NS_2D)
1000#if defined(DNDS_DIST_MT_USE_OMP)
1001# pragma omp parallel for schedule(runtime)
1002#endif
1003 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
1004 {
1005 Geom::tPoint pos = vfv->GetCellBary(iCell);
1006 real gamma = settings.idealGasProperty.gamma;
1007 real rho = 2;
1008 real p = 1 + 2 * pos(1);
1009 if (pos(1) >= 0.5)
1010 {
1011 rho = 1;
1012 p = 1.5 + pos(1);
1013 }
1014 real v = -0.025 * sqrt(gamma * p / rho) * std::cos(8 * pi * pos(0));
1015 if constexpr (dim == 3)
1016 u[iCell] = Eigen::Vector<real, 5>{rho, 0, rho * v, 0, 0.5 * rho * sqr(v) + p / (gamma - 1)};
1017 else
1018 u[iCell] = Eigen::Vector<real, 4>{rho, 0, rho * v, 0.5 * rho * sqr(v) + p / (gamma - 1)};
1019 }
1020 else if constexpr (model == NS_3D)
1021#if defined(DNDS_DIST_MT_USE_OMP)
1022# pragma omp parallel for schedule(runtime)
1023#endif
1024 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
1025 {
1026 Geom::tPoint pos = vfv->GetCellBary(iCell);
1027 real gamma = settings.idealGasProperty.gamma;
1028 real rho = 2;
1029 real p = 1 + 2 * pos(1);
1030 if (pos(1) >= 0.5)
1031 {
1032 rho = 1;
1033 p = 1.5 + pos(1);
1034 }
1035 real v = -0.025 * sqrt(gamma * p / rho) * std::cos(8 * pi * pos(0)) * std::cos(8 * pi * pos(2));
1036 u[iCell] = Eigen::Vector<real, 5>{rho, 0, rho * v, 0, 0.5 * rho * sqr(v) + p / (gamma - 1)};
1037 }
1038 break;
1039 case 2: // for IV10 problem
1040 case 203: // for IV10 problem with PP
1041 case 202:
1042 DNDS_assert(model == NS || model == NS_2D || model == NS_EX);
1043 if constexpr (model == NS || model == NS_2D || model == NS_EX)
1044#if defined(DNDS_DIST_MT_USE_OMP)
1045# pragma omp parallel for schedule(runtime)
1046#endif
1047 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
1048 {
1049 Geom::tPoint pos = vfv->GetCellBary(iCell);
1050 auto c2n = mesh->cell2node[iCell];
1051 auto gCell = vfv->GetCellQuad(iCell);
1052 TU um;
1053 um.resizeLike(u[iCell]);
1054 um.setZero();
1055 gCell.IntegrationSimple(
1056 um,
1057 [&](TU &inc, int ig)
1058 {
1059 Geom::tPoint pPhysics = vfv->GetCellQuadraturePPhys(iCell, ig);
1060 if (settings.specialBuiltinInitializer == 2)
1061 inc = SpecialFields::IsentropicVortex10(*this, pPhysics, 0, nVars, 5);
1062 else if (settings.specialBuiltinInitializer == 203)
1063 inc = SpecialFields::IsentropicVortex10(*this, pPhysics, 0, nVars, 10.0828);
1064 else if (settings.specialBuiltinInitializer == 202)
1065 inc = SpecialFields::IsentropicVortex30(*this, pPhysics, 0, nVars);
1066 else
1067 DNDS_assert(false);
1068 inc *= vfv->GetCellJacobiDet(iCell, ig); // don't forget this
1069 });
1070 u[iCell] = um / vfv->GetCellVol(iCell); // mean value
1071 }
1072 break;
1073 case 201: // for IVCent problem
1074 DNDS_assert(model == NS || model == NS_2D);
1075 if constexpr (model == NS || model == NS_2D)
1076#if defined(DNDS_DIST_MT_USE_OMP)
1077# pragma omp parallel for schedule(runtime)
1078#endif
1079 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
1080 {
1081 Geom::tPoint pos = vfv->GetCellBary(iCell);
1082 auto c2n = mesh->cell2node[iCell];
1083 auto gCell = vfv->GetCellQuad(iCell);
1084 TU um;
1085 um.resizeLike(u[iCell]);
1086 um.setZero();
1087 gCell.IntegrationSimple(
1088 um,
1089 [&](TU &inc, int ig)
1090 {
1091 Geom::tPoint pPhysics = vfv->GetCellQuadraturePPhys(iCell, ig);
1092 inc = SpecialFields::IsentropicVortexCent(*this, pPhysics, 0, nVars);
1093 inc *= vfv->GetCellJacobiDet(iCell, ig); // don't forget this
1094 });
1095 u[iCell] = um / vfv->GetCellVol(iCell); // mean value
1096 }
1097 break;
1098 case 3: // for taylor-green vortex problem
1099 DNDS_assert(model == NS_3D);
1100 if constexpr (model == NS_3D)
1101 {
1102#if defined(DNDS_DIST_MT_USE_OMP)
1103# pragma omp parallel for schedule(runtime)
1104#endif
1105 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
1106 {
1107 Geom::tPoint pos = vfv->GetCellBary(iCell);
1108 real M0 = 0.1;
1109 real gamma = settings.idealGasProperty.gamma;
1110 auto c2n = mesh->cell2node[iCell];
1111 auto gCell = vfv->GetCellQuad(iCell);
1112 TU um;
1113 um.resizeLike(u[iCell]);
1114 um.setZero();
1115 // Eigen::MatrixXd coords;
1116 // mesh->GetCoords(c2n, coords);
1117 gCell.IntegrationSimple(
1118 um,
1119 [&](TU &inc, int ig)
1120 {
1121 // std::cout << coords<< std::endl << std::endl;
1122 // std::cout << DiNj << std::endl;
1123 Geom::tPoint pPhysics = vfv->GetCellQuadraturePPhys(iCell, ig);
1124 real x{pPhysics(0)}, y{pPhysics(1)}, z{pPhysics(2)};
1125 real ux = std::sin(x) * std::cos(y) * std::cos(z);
1126 real uy = -std::cos(x) * std::sin(y) * std::cos(z);
1127 real p = 1. / (gamma * sqr(M0)) + 1. / 16 * ((std::cos(2 * x) + std::cos(2 * y)) * (2 + std::cos(2 * z)));
1128 real rho = gamma * sqr(M0) * p;
1129 real E = 0.5 * (sqr(ux) + sqr(uy)) * rho + p / (gamma - 1);
1130
1131 // std::cout << T << " " << rho << std::endl;
1132 inc.setZero();
1133 inc(0) = rho;
1134 inc(1) = rho * ux;
1135 inc(2) = rho * uy;
1136 inc(dim + 1) = E;
1137
1138 inc *= vfv->GetCellJacobiDet(iCell, ig); // don't forget this
1139 });
1140 u[iCell] = um / vfv->GetCellVol(iCell); // mean value
1141 }
1142 }
1143 break;
1144 case 3001: // for nol problem
1145 DNDS_assert(model == NS_3D);
1146#if defined(DNDS_DIST_MT_USE_OMP)
1147# pragma omp parallel for schedule(runtime)
1148#endif
1149 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
1150 {
1151 Geom::tPoint pos = vfv->GetCellBary(iCell);
1152 real M0 = 0.1;
1153 real gamma = settings.idealGasProperty.gamma;
1154 auto c2n = mesh->cell2node[iCell];
1155 auto gCell = vfv->GetCellQuad(iCell);
1156 TU um;
1157 um.resizeLike(u[iCell]);
1158 um.setZero();
1159 // Eigen::MatrixXd coords;
1160 // mesh->GetCoords(c2n, coords);
1161 gCell.IntegrationSimple(
1162 um,
1163 [&](TU &inc, int ig)
1164 {
1165 // std::cout << coords<< std::endl << std::endl;
1166 // std::cout << DiNj << std::endl;
1167 Geom::tPoint pPhysics = vfv->GetCellQuadraturePPhys(iCell, ig);
1168 TU farPrimitive;
1169 Gas::IdealGasThermalConservative2Primitive<dim>(settings.farFieldStaticValue, farPrimitive, settings.idealGasProperty.gamma);
1170 real pInf = farPrimitive(I4);
1171 real r = pPhysics.norm();
1172 TVec velo = -pPhysics(Seq012) / (r + smallReal);
1173 farPrimitive(I4) = pInf;
1174 farPrimitive(0) = 1;
1175 farPrimitive(Seq123) = velo;
1176
1177 Gas::IdealGasThermalPrimitive2Conservative<dim>(farPrimitive, inc, settings.idealGasProperty.gamma);
1178 inc *= vfv->GetCellJacobiDet(iCell, ig); // don't forget this
1179 });
1180 u[iCell] = um / vfv->GetCellVol(iCell); // mean value
1181 }
1182 break;
1183 case 0:
1184 break;
1185 default:
1186 log() << "Wrong specialBuiltinInitializer" << std::endl;
1187 DNDS_assert(false);
1188 break;
1189 }
1190
1191 if (settings.frameConstRotation.enabled)
1192 {
1193#if defined(DNDS_DIST_MT_USE_OMP)
1194# pragma omp parallel for schedule(runtime)
1195#endif
1196 for (int iCell = 0; iCell < mesh->NumCell(); iCell++)
1197 {
1198 TU ui = u[iCell];
1199 TransformURotatingFrame(ui, vfv->GetCellQuadraturePPhys(iCell, -1), -1);
1200 u[iCell] = ui;
1201 }
1202 }
1203
1204 // Box
1205 for (auto &i : settings.boxInitializers)
1206 {
1207#if defined(DNDS_DIST_MT_USE_OMP)
1208# pragma omp parallel for schedule(runtime)
1209#endif
1210 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
1211 {
1212 Geom::tPoint pos = vfv->GetCellBary(iCell);
1213 if (pos(0) > i.x0 && pos(0) < i.x1 &&
1214 pos(1) > i.y0 && pos(1) < i.y1 &&
1215 pos(2) > i.z0 && pos(2) < i.z1)
1216 {
1217 u[iCell] = i.v;
1218 }
1219 }
1220 }
1221
1222 // Plane
1223 for (auto &i : settings.planeInitializers)
1224 {
1225#if defined(DNDS_DIST_MT_USE_OMP)
1226# pragma omp parallel for schedule(runtime)
1227#endif
1228 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
1229 {
1230 Geom::tPoint pos = vfv->GetCellBary(iCell);
1231 if (pos(0) * i.a + pos(1) * i.b + pos(2) * i.c + i.h > 0)
1232 {
1233 // std::cout << pos << std::endl << i.a << i.b << std::endl << i.h <<std::endl;
1234 // DNDS_assert(false);
1235 u[iCell] = i.v;
1236 }
1237 }
1238 }
1239
1240 for (auto &i : settings.exprtkInitializers)
1241 {
1242 auto exprStr = i.GetExpr();
1243 ExprtkWrapperEvaluator exprtkEval;
1244 exprtkEval.AddScalar("inRegion");
1245 exprtkEval.AddScalar("iCell");
1246 exprtkEval.AddVector("x", dim);
1247 exprtkEval.AddVector("UPrim", nVars);
1248 exprtkEval.Compile(exprStr);
1249#if defined(DNDS_DIST_MT_USE_OMP)
1250# pragma omp parallel for schedule(runtime)
1251#endif
1252 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
1253 {
1254 Geom::tPoint pos = vfv->GetCellBary(iCell);
1255 auto c2n = mesh->cell2node[iCell];
1256 auto gCell = vfv->GetCellQuad(iCell);
1257 TU um;
1258 um.resizeLike(u[iCell]);
1259 um.setZero();
1260 bool allIn = true;
1261 bool someIn = false;
1262 gCell.IntegrationSimple(
1263 um,
1264 [&](TU &inc, int ig)
1265 {
1266 Geom::tPoint pPhysics = vfv->GetCellQuadraturePPhys(iCell, ig);
1267
1268 exprtkEval.Var("inRegion") = 0;
1269 exprtkEval.Var("iCell") = real(iCell);
1270 exprtkEval.VarVec("x", 0) = pPhysics(0);
1271 exprtkEval.VarVec("x", 1) = pPhysics(1);
1272 if constexpr (dim == 3)
1273 exprtkEval.VarVec("x", 2) = pPhysics(2);
1274
1275 TU uPrimitive;
1276 Gas::IdealGasThermalConservative2Primitive<dim>(u[iCell], uPrimitive, settings.idealGasProperty.gamma);
1277 for (int i = 0; i < nVars; i++)
1278 exprtkEval.VarVec("UPrim", i) = uPrimitive(i);
1279
1280 real ret = exprtkEval.Evaluate();
1281
1282 DNDS_check_throw_info(ret == 0.0, "return of \n" + exprStr + "\nis non-zero");
1283 allIn = allIn && exprtkEval.Var("inRegion");
1284 someIn = someIn || exprtkEval.Var("inRegion");
1285
1286 if (exprtkEval.Var("inRegion"))
1287 for (int i = 0; i < nVars; i++)
1288 uPrimitive(i) = exprtkEval.VarVec("UPrim", i);
1289 Gas::IdealGasThermalPrimitive2Conservative<dim>(uPrimitive, inc, settings.idealGasProperty.gamma);
1290 if (!inc.allFinite())
1291 {
1292 std::ostringstream oss0, oss1, oss2;
1293 oss0 << uPrimitive.transpose();
1294 oss1 << inc.transpose();
1295 oss2 << "x=" << pPhysics.transpose();
1296 DNDS_check_throw_info(
1297 false,
1298 fmt::format("Got invalid exprtk product: [{}] -> [{}]\n{}\n {}\n",
1299 oss0.str(), oss1.str(), oss2.str(), exprStr));
1300 }
1301 inc *= vfv->GetCellJacobiDet(iCell, ig); // don't forget this
1302 });
1303 // if (allIn)
1304 if (someIn)
1305 u[iCell] = um / vfv->GetCellVol(iCell); // mean value
1306 }
1307 }
1308 }
1309
1310 template <EulerModel model>
1311 /** @brief Apply maximum-value spatial filter to conservative variables.
1312 *
1313 * Placeholder for a spatial filter Jacobian approach. Currently returns immediately
1314 * without modification.
1315 *
1316 * @param u Conservative variable DOF array (unmodified).
1317 */
1318 void EulerEvaluator<model>::FixUMaxFilter(ArrayDOFV<nVarsFixed> &u)
1319 {
1320 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
1321 // TODO: make spacial filter jacobian
1322 return; // ! nofix shortcut
1323 }
1324
1325 template <EulerModel model>
1326 /** @brief Accumulate time-averaged solution field using running weighted average.
1327 *
1328 * Updates wAveraged as a running weighted average: wAveraged = (tCur*wAveraged + dt*w) / (tCur+dt).
1329 *
1330 * @param w Current primitive variable field.
1331 * @param wAveraged Running time-averaged field (input/output).
1332 * @param dt Time step weight for current sample.
1333 * @param tCur Accumulated averaging time (input/output, incremented by dt).
1334 */
1335 void EulerEvaluator<model>::TimeAverageAddition(ArrayDOFV<nVarsFixed> &w, ArrayDOFV<nVarsFixed> &wAveraged, real dt, real &tCur)
1336 {
1337 wAveraged *= (tCur / (tCur + dt + verySmallReal));
1338 wAveraged.addTo(w, (dt + verySmallReal) / (tCur + dt + verySmallReal));
1339 tCur += dt + verySmallReal;
1340 }
1341
1342 template <EulerModel model>
1343 /** @brief Convert cell mean values from conservative to primitive variables.
1344 *
1345 * @param u Conservative variable DOF array (input).
1346 * @param w Primitive variable DOF array (output).
1347 */
1348 void EulerEvaluator<model>::MeanValueCons2Prim(ArrayDOFV<nVarsFixed> &u, ArrayDOFV<nVarsFixed> &w)
1349 {
1350#if defined(DNDS_DIST_MT_USE_OMP)
1351# pragma omp parallel for schedule(static)
1352#endif
1353 for (index iCell = 0; iCell < u.Size(); iCell++)
1354 {
1355 real gamma = settings.idealGasProperty.gamma;
1356 TU out;
1357 Gas::IdealGasThermalConservative2Primitive<dim>(u[iCell], out, gamma);
1358 w[iCell] = out;
1359 }
1360 }
1361
1362 template <EulerModel model>
1363 /** @brief Convert cell mean values from primitive to conservative variables.
1364 *
1365 * @param w Primitive variable DOF array (input).
1366 * @param u Conservative variable DOF array (output).
1367 */
1368 void EulerEvaluator<model>::MeanValuePrim2Cons(ArrayDOFV<nVarsFixed> &w, ArrayDOFV<nVarsFixed> &u)
1369 {
1370#if defined(DNDS_DIST_MT_USE_OMP)
1371# pragma omp parallel for schedule(static)
1372#endif
1373 for (index iCell = 0; iCell < w.Size(); iCell++)
1374 {
1375 real gamma = settings.idealGasProperty.gamma;
1376 TU out;
1377 Gas::IdealGasThermalPrimitive2Conservative<dim>(w[iCell], out, gamma);
1378 u[iCell] = out;
1379 }
1380 }
1381
1382 template <EulerModel model>
1383 /** @brief Evaluate the LP norm of the RHS residual (MPI-global).
1384 *
1385 * For P < 3, computes component-wise LP norm; for P >= 3, computes L-infinity norm.
1386 * Optionally volume-weighted and/or averaged.
1387 *
1388 * @param res Result vector of nVars component norms (output).
1389 * @param rhs Right-hand side residual.
1390 * @param P Norm exponent (1 for L1, 2 for L2, >=3 for L-infinity).
1391 * @param volWise If true, weight by cell volume.
1392 * @param average If true, divide by total weight.
1393 */
1394 void EulerEvaluator<model>::EvaluateNorm(Eigen::Vector<real, -1> &res, ArrayDOFV<nVarsFixed> &rhs, index P, bool volWise, bool average)
1395 {
1396 res.resize(nVars);
1397 if (P < 3)
1398 {
1399 TU resc;
1400 resc.setZero(nVars);
1401 real rescBase{0};
1402#if defined(DNDS_DIST_MT_USE_OMP)
1403# pragma omp declare reduction(TUAdd:TU : omp_out += omp_in) initializer(omp_priv = omp_orig)
1404# pragma omp parallel for schedule(static) reduction(TUAdd : resc)
1405#endif
1406 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
1407 {
1408 if (rhs[iCell].hasNaN() || (!rhs[iCell].allFinite()))
1409 {
1410 std::cout << rhs[iCell] << std::endl;
1411 DNDS_assert(false);
1412 }
1413 if (volWise)
1414 resc += rhs[iCell].array().abs().pow(P).matrix() * vfv->GetCellVol(iCell), rescBase += vfv->GetCellVol(iCell);
1415 else
1416 resc += rhs[iCell].array().abs().pow(P).matrix(), rescBase += 1;
1417 }
1418 MPI::Allreduce(resc.data(), res.data(), res.size(), DNDS_MPI_REAL, MPI_SUM, rhs.father->getMPI().comm);
1419 if (average)
1420 MPI::AllreduceOneReal(rescBase, MPI_SUM, rhs.father->getMPI());
1421 res = res.array().pow(1.0 / P).matrix();
1422 if (average)
1423 res *= 1. / (rescBase + verySmallReal);
1424 // std::cout << res << std::endl;
1425 }
1426 else
1427 {
1428 TU resc;
1429 resc.resizeLike(rhs[0]);
1430 resc.setZero();
1431#if defined(DNDS_DIST_MT_USE_OMP)
1432# pragma omp declare reduction(TUAdd:TU : omp_out = omp_out.array().max(omp_in.array())) initializer(omp_priv = omp_orig)
1433# pragma omp parallel for schedule(static) reduction(TUAdd : resc)
1434#endif
1435 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
1436 resc = resc.array().max(rhs[iCell].array().abs()).matrix();
1437 MPI::Allreduce(resc.data(), res.data(), res.size(), DNDS_MPI_REAL, MPI_MAX, rhs.father->getMPI().comm);
1438 }
1439 }
1440
1441 template <class TU>
1442 struct TU_P_Reduction
1443 {
1444 TU _v;
1445 index P{1};
1446
1447 void setZero(int n)
1448 {
1449 _v.setZero(n);
1450 }
1451
1452 void reduce(const TU_P_Reduction<TU> &R)
1453 {
1454 _v = ((P < 3) ? TU(_v + R._v) : TU((_v.array().max(R._v.array())).matrix()));
1455 }
1456 void setP(index nP) { P = nP; }
1457 };
1458
1459 template <EulerModel model>
1460 /** @brief Evaluate the LP norm of the reconstructed solution, optionally against a reference field.
1461 *
1462 * Integrates the reconstructed solution (or its error against FCompareField) over
1463 * quadrature points in each cell and reduces globally.
1464 *
1465 * @param res Result vector of nVars component norms (output).
1466 * @param u Conservative variable DOF array.
1467 * @param uRec Reconstruction coefficients.
1468 * @param P Norm exponent (1 for L1, 2 for L2, >=3 for L-infinity).
1469 * @param compare If true, compute error against FCompareField.
1470 * @param FCompareField Reference field function(point, time) -> TU.
1471 * @param FCompareFieldWeight Weight function(point, time) -> TU for error weighting.
1472 * @param t Current simulation time.
1473 */
1474 void EulerEvaluator<model>::EvaluateRecNorm(
1475 Eigen::Vector<real, -1> &res,
1476 ArrayDOFV<nVarsFixed> &u,
1477 ArrayRECV<nVarsFixed> &uRec,
1478 index P,
1479 bool compare,
1480 const tFCompareField &FCompareField,
1481 const tFCompareFieldWeight &FCompareFieldWeight,
1482 real t)
1483 {
1484 res.resize(nVars);
1485 TU_P_Reduction<TU> resc;
1486 resc.setZero(nVars);
1487 resc.setP(P);
1488#if defined(DNDS_DIST_MT_USE_OMP)
1489# pragma omp declare reduction(TUAdd : TU_P_Reduction<TU> : omp_out.reduce(omp_in)) initializer(omp_priv = omp_orig)
1490# pragma omp parallel for schedule(static) reduction(TUAdd : resc)
1491#endif
1492 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
1493 {
1494 auto qCell = vfv->GetCellQuad(iCell);
1495 TU rescCell;
1496 rescCell.setZero(nVars);
1497 qCell.IntegrationSimple(
1498 rescCell,
1499 [&](auto &inc, int iG)
1500 {
1501 TU uR = u[iCell] + (vfv->GetIntPointDiffBaseValue(iCell, -1, -1, iG, 0, 1) * uRec[iCell]).transpose();
1502 if (compare)
1503 {
1504 Geom::tPoint pPhysics = vfv->GetCellQuadraturePPhys(iCell, iG);
1505 uR -= FCompareField(pPhysics, t);
1506 uR *= FCompareFieldWeight(pPhysics, t);
1507 }
1508 if (P >= 3)
1509 resc._v = resc._v.array().max(uR.array().abs());
1510 inc = uR.array().abs().pow(P);
1511 inc *= vfv->GetCellJacobiDet(iCell, iG);
1512 });
1513 if (P < 3)
1514 resc._v += rescCell;
1515 }
1516 if (P > 3)
1517 MPI::Allreduce(resc._v.data(), res.data(), res.size(), DNDS_MPI_REAL, MPI_MAX, u.father->getMPI().comm);
1518 else
1519 {
1520 MPI::Allreduce(resc._v.data(), res.data(), res.size(), DNDS_MPI_REAL, MPI_SUM, u.father->getMPI().comm);
1521 res = res.array().pow(1.0 / P).matrix();
1522 }
1523 }
1524
1525 DNDS_SWITCH_INTELLISENSE(
1526 template <EulerModel model>, template <>)
1527 /** @brief Apply positivity-preserving gradient limiter to reconstructed solution gradients.
1528 *
1529 * Limits the gradient so that density and internal energy remain above threshold
1530 * values at all face quadrature points (Zhang-Shu style for gradients).
1531 * Optionally disables shock-based limiting via flags.
1532 *
1533 * @param u Conservative variable DOF array.
1534 * @param uGrad Input gradient array.
1535 * @param uGradNew Limited gradient array (output).
1536 * @param flags Bitfield flags (e.g., LIMITER_UGRAD_Disable_Shock_Limiter).
1537 */
1538 void EulerEvaluator<model>::LimiterUGrad(
1539 ArrayDOFV<nVarsFixed> &u, ArrayGRADV<nVarsFixed, gDim> &uGrad, ArrayGRADV<nVarsFixed, gDim> &uGradNew,
1540 uint64_t flags)
1541 {
1542 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
1543 static const real safetyRatio = 1 - 1e-2;
1544 static const real E_lb_eps = smallReal;
1545
1546 bool disable_shock_limiter = flags & LIMITER_UGRAD_Disable_Shock_Limiter;
1547
1548#if defined(DNDS_DIST_MT_USE_OMP)
1549# pragma omp parallel for schedule(runtime)
1550#endif
1551 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
1552 {
1553 uGradNew[iCell] = uGrad[iCell];
1554 auto c2f = mesh->cell2face[iCell];
1555
1556 TU_Batch uFaceInc;
1557 uFaceInc.setZero(nVars, c2f.size() * 2); // j < c2f.size(): faceInc; j > c2f.size(): baryInc
1558 TU uOtherMin = u[iCell];
1559 TU uOtherMax = u[iCell];
1560 auto fEInternal = [](const TU &u) -> real
1561 { return u(I4) - 0.5 * u(Seq123).squaredNorm() / (u(0) + verySmallReal); };
1562 real eOtherMin, eOtherMax;
1563 eOtherMin = eOtherMax = fEInternal(u[iCell]);
1564 for (rowsize ic2f = 0; ic2f < c2f.size(); ic2f++)
1565 {
1566 index iFace = c2f[ic2f];
1567 uFaceInc(EigenAll, ic2f) =
1568 uGrad[iCell].transpose() *
1569 (vfv->GetFaceQuadraturePPhysFromCell(iFace, iCell, -1, -1) - vfv->GetCellQuadraturePPhys(iCell, -1))(SeqG012);
1570 index iCellOther = mesh->CellFaceOther(iCell, iFace);
1571 if (iCellOther != UnInitIndex)
1572 {
1573 uOtherMin = uOtherMin.array().min(u[iCellOther].array());
1574 uOtherMax = uOtherMin.array().max(u[iCellOther].array());
1575 eOtherMin = std::min(eOtherMin, fEInternal(u[iCellOther]));
1576 eOtherMax = std::max(eOtherMax, fEInternal(u[iCellOther]));
1577 // ! adding bary value if the other cell exists
1578 // uFaceInc(EigenAll, ic2f + c2f.size()) =
1579 // uGrad[iCell].transpose() *
1580 // (vfv->GetOtherCellBaryFromCell(iCell, iCellOther, -1) - vfv->GetCellQuadraturePPhys(iCell, -1))(SeqG012);
1581 }
1582 }
1583
1584 TU uFaceIncMax = uFaceInc.array().rowwise().maxCoeff();
1585 TU uFaceIncMin = uFaceInc.array().rowwise().minCoeff();
1586 if (!disable_shock_limiter)
1587 {
1588 TU alpha0;
1589 alpha0.setConstant(nVars, 1.0);
1590 alpha0 = alpha0.array().min(((uOtherMax - u[iCell]).array().abs() / (uFaceIncMax.array().abs() + verySmallReal)));
1591 alpha0 = alpha0.array().min(((uOtherMin - u[iCell]).array().abs() / (uFaceIncMin.array().abs() + verySmallReal)));
1592 uGradNew[iCell].array().rowwise() *= alpha0.array().transpose();
1593 uFaceInc.array().colwise() *= alpha0.array();
1594 }
1595
1596 // start PP
1597 real alphaPP_Rho = 1.0;
1598 if (disable_shock_limiter) // do rho PP first
1599 {
1600 alphaPP_Rho = std::min(alphaPP_Rho, u[iCell][0] / (std::abs(uFaceIncMin(0)) + smallReal * u[iCell][0]));
1601 if (alphaPP_Rho < 1.0)
1602 alphaPP_Rho *= safetyRatio;
1603 uGradNew[iCell].array() *= alphaPP_Rho;
1604 uFaceInc.array() *= alphaPP_Rho;
1605 }
1606
1607 TU_Batch uFaceAlpha0 = uFaceInc.colwise() + u[iCell];
1608 for (int j = 0; j < uFaceAlpha0.cols(); j++)
1609 DNDS_assert(uFaceAlpha0(0, j) > 0);
1610 real minEFace =
1611 (uFaceAlpha0(I4, EigenAll).array() -
1612 0.5 * uFaceAlpha0(Seq123, EigenAll).colwise().squaredNorm().array() / (uFaceAlpha0(0, EigenAll).array() + verySmallReal))
1613 .minCoeff();
1614 real eC = fEInternal(u[iCell]);
1615 real deltaEFaceMin = minEFace - eC;
1616 real alphaPP_E = 1.0;
1617 if (deltaEFaceMin < 0)
1618 alphaPP_E = std::min(alphaPP_E, std::abs(eC * (1 - E_lb_eps)) / (verySmallReal - deltaEFaceMin));
1619 if (alphaPP_E < 1.0)
1620 alphaPP_E *= safetyRatio;
1621 uGradNew[iCell](EigenAll, Seq01234) *= alphaPP_E;
1622 }
1623 }
1624
1625 DNDS_SWITCH_INTELLISENSE(
1626 template <EulerModel model>, template <>)
1627 /** @brief Evaluate the positivity-preserving reconstruction limiter coefficient (beta) per cell.
1628 *
1629 * Computes a scalar beta in [0,1] for each cell such that the reconstructed
1630 * solution at all face and (optionally) volume quadrature points maintains
1631 * positive density and pressure. Uses bisection to find the maximum allowable beta.
1632 * Reports the total number of limited cells and the global minimum beta.
1633 *
1634 * @param u Conservative variable DOF array.
1635 * @param uRec Reconstruction coefficients (may be modified for order reduction).
1636 * @param uRecBeta Per-cell limiter coefficient (output).
1637 * @param nLim Total number of limited cells across all MPI ranks (output).
1638 * @param betaMin Global minimum beta value (output).
1639 * @param flag Evaluation mode flags.
1640 */
1641 void EulerEvaluator<model>::EvaluateURecBeta(
1642 ArrayDOFV<nVarsFixed> &u,
1643 ArrayRECV<nVarsFixed> &uRec,
1644 ArrayDOFV<1> &uRecBeta, index &nLim, real &betaMin,
1645 int flag)
1646 {
1647 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
1648 static const real safetyRatio = 1 - 1e-5;
1649 static const real minRatio = 0.5;
1650 real rhoEps = smallReal * settings.refUPrim(0) * 1e-1;
1651 real pEps = smallReal * settings.refUPrim(I4) * 1e-1;
1652 real betaCutOff = 1e-3;
1653 bool restrictOnVolPoints = (!settings.ignoreSourceTerm) || settings.forceVolURecBeta;
1654
1655 if (settings.ppEpsIsRelaxed)
1656 {
1657 real rhoMin = veryLargeReal;
1658 real pMin = veryLargeReal;
1659#if defined(DNDS_DIST_MT_USE_OMP)
1660# pragma omp parallel for schedule(runtime) reduction(min : rhoMin, pMin)
1661#endif
1662 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
1663 {
1664 TU UPrim;
1665 Gas::IdealGasThermalConservative2Primitive<dim>(u[iCell], UPrim, settings.idealGasProperty.gamma);
1666 rhoMin = std::min(rhoMin, UPrim(0));
1667 pMin = std::min(pMin, UPrim(I4));
1668 }
1669 MPI::AllreduceOneReal(rhoMin, MPI_MIN, mesh->getMPI());
1670 MPI::AllreduceOneReal(pMin, MPI_MIN, mesh->getMPI());
1671 rhoEps = std::min(rhoEps, minRatio * rhoMin);
1672 pEps = std::min(pEps, minRatio * pMin);
1673 }
1674
1675 index nLimLocal = 0;
1676 real minBetaLocal = 1;
1677#if defined(DNDS_DIST_MT_USE_OMP)
1678# pragma omp parallel for schedule(runtime) reduction(+ : nLimLocal) reduction(min : minBetaLocal)
1679#endif
1680 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
1681 {
1682 auto gCell = vfv->GetCellQuad(iCell);
1683 int nPoint = restrictOnVolPoints ? gCell.GetNumPoints() : 0;
1684 auto c2f = mesh->cell2face[iCell];
1685 for (int ic2f = 0; ic2f < c2f.size(); ic2f++)
1686 nPoint += vfv->GetFaceQuad(c2f[ic2f]).GetNumPoints();
1687 /***********/
1688 MatrixXR quadBase;
1689 quadBase.resize(nPoint, vfv->GetCellAtr(iCell).NDOF - 1);
1690 nPoint = 0;
1691 if (restrictOnVolPoints)
1692 {
1693 for (int iG = 0; iG < gCell.GetNumPoints(); iG++)
1694 quadBase(iG, EigenAll) = vfv->GetIntPointDiffBaseValue(iCell, -1, -1, iG, 0, 1);
1695 nPoint += gCell.GetNumPoints();
1696 }
1697 for (int ic2f = 0; ic2f < c2f.size(); ic2f++)
1698 {
1699 auto gFace = vfv->GetFaceQuad(c2f[ic2f]);
1700 for (int iG = 0; iG < gFace.GetNumPoints(); iG++)
1701 quadBase(nPoint + iG, EigenAll) = vfv->GetIntPointDiffBaseValue(iCell, c2f[ic2f], -1, iG, 0, 1);
1702 nPoint += gFace.GetNumPoints();
1703 }
1704 /***********/
1705 DNDS_assert_info(u[iCell](0) >= rhoEps, fmt::format("rhoMean {}, {}", u[iCell](0), rhoEps));
1706 real gamma = settings.idealGasProperty.gamma;
1707 real pCent = (gamma - 1) * (u[iCell](I4) - 0.5 * u[iCell](Seq123).squaredNorm() / u[iCell](0));
1708 DNDS_assert_info(pCent >= pEps, fmt::format("pMean {}, {}", pCent, pEps));
1709
1710 // alter uRec if necessary
1711 int curOrder = vfv->GetCellOrder(iCell);
1712 Eigen::Matrix<real, Eigen::Dynamic, nVarsFixed> uRecBase = uRec[iCell];
1713 Eigen::Matrix<real, Eigen::Dynamic, nVarsFixed> recBase; // * has to call checkRecBaseGood() to hold valid value
1714 auto checkRecBaseGood = [&]()
1715 {
1716 recBase = (quadBase * uRecBase).rowwise() + u[iCell].transpose();
1717 if (recBase(EigenAll, 0).minCoeff() < rhoEps) // TODO: add relaxation to eps values
1718 return false;
1719 if constexpr (Traits::hasSA)
1720 if (recBase(EigenAll, I4 + 1).minCoeff() < rhoEps)
1721 return false;
1722 if constexpr (Traits::has2EQ)
1723 if (recBase(EigenAll, I4 + 1).minCoeff() < rhoEps || recBase(EigenAll, I4 + 2).minCoeff() < rhoEps)
1724 return false;
1725 Eigen::Vector<real, Eigen::Dynamic> ek =
1726 0.5 * (recBase(EigenAll, Seq123).array().square().rowwise().sum()) / recBase(EigenAll, 0).array();
1727 Eigen::Vector<real, Eigen::Dynamic> eInternalS = (recBase(EigenAll, I4) - ek);
1728 if (eInternalS.minCoeff() < pEps)
1729 return false;
1730 return true;
1731 };
1732 if (checkRecBaseGood())
1733 {
1734 uRecBeta[iCell](0) = 1;
1735 continue; //! early exit, reconstruction is good it self
1736 }
1737 if (flag & EvaluateURecBeta_COMPRESS_TO_MEAN)
1738 curOrder = 1;
1739 while (curOrder > 0)
1740 {
1741 uRecBase = vfv->template DownCastURecOrder<nVarsFixed>(curOrder, iCell, uRec, 0);
1742 if (checkRecBaseGood())
1743 break;
1744 uRec[iCell] = uRecBase; // uRec[iCell] could be altered
1745 curOrder--;
1746 }
1747
1748 Eigen::Matrix<real, Eigen::Dynamic, nVarsFixed>
1749 recInc = quadBase * (uRec[iCell] - uRecBase);
1750 Eigen::Vector<real, Eigen::Dynamic> rhoS = recInc(EigenAll, 0) + recBase(EigenAll, 0);
1751 Eigen::Index rhoMinIdx;
1752 real rhoMin = rhoS.minCoeff(&rhoMinIdx);
1753 real theta1 = 1;
1754 if (rhoMin < rhoEps)
1755 for (int iG = 0; iG < rhoS.size(); iG++)
1756 if (recInc(iG, 0) < 0) // negative increment
1757 theta1 = std::min(theta1, (recBase(iG, 0) - rhoEps) / (-recInc(iG, 0) + verySmallReal));
1758#ifdef USE_NS_SA_NUT_REDUCED_ORDER
1759 if constexpr (Traits::hasSA)
1760 {
1761 static real v1Eps = smallReal * settings.refUPrim(I4 + 1);
1762 Eigen::Vector<real, Eigen::Dynamic> v1S = recInc(EigenAll, I4 + 1) + recBase(EigenAll, I4 + 1);
1763 real v1Min = v1S.minCoeff();
1764 if (v1Min < v1Eps)
1765 for (int iG = 0; iG < rhoS.size(); iG++)
1766 if (recInc(iG, I4 + 1) < 0) // negative increment
1767 theta1 = std::min(theta1, (recBase(iG, I4 + 1) - v1Eps) / (-recInc(iG, I4 + 1) + verySmallReal))
1768 // * 0 // to gain fully reduced order
1769 ;
1770 }
1771#endif
1772 if constexpr (Traits::has2EQ)
1773 {
1774 real thetaC = 1;
1775 static real v1Eps = smallReal * settings.refUPrim(I4 + 1);
1776 static real v2Eps = smallReal * settings.refUPrim(I4 + 2);
1777 Eigen::Vector<real, Eigen::Dynamic> v1S = recInc(EigenAll, I4 + 1) + recBase(EigenAll, I4 + 1);
1778 Eigen::Vector<real, Eigen::Dynamic> v2S = recInc(EigenAll, I4 + 2) + recBase(EigenAll, I4 + 2);
1779 real v1Min = v1S.minCoeff();
1780 real v2Min = v2S.minCoeff();
1781 if (v1Min < v1Eps)
1782 for (int iG = 0; iG < rhoS.size(); iG++)
1783 if (recInc(iG, I4 + 1) < 0) // negative increment
1784 thetaC = std::min(thetaC, (recBase(iG, I4 + 1) - v1Eps) / (-recInc(iG, I4 + 1) + verySmallReal));
1785 if (v2Min < v2Eps)
1786 for (int iG = 0; iG < rhoS.size(); iG++)
1787 if (recInc(iG, I4 + 2) < 0) // negative increment
1788 thetaC = std::min(thetaC, (recBase(iG, I4 + 2) - v2Eps) / (-recInc(iG, I4 + 2) + verySmallReal));
1789 // theta1 = std::min(theta1, thetaC);
1790 if (thetaC < 1) // 2eq's pp not disturbing main flow
1791 {
1792 uRec[iCell](EigenAll, {I4 + 1, I4 + 2}) *= thetaC * 0.9;
1793 }
1794 }
1795
1796 if constexpr (Traits::hasRANS)
1797 recInc(EigenAll, Seq01234) *= theta1; // to leave SA unchanged
1798 else
1799 recInc *= theta1;
1800 Eigen::Matrix<real, Eigen::Dynamic, nVarsFixed>
1801 recVRhoG = recInc + recBase;
1802
1803 Eigen::Vector<real, Eigen::Dynamic> ek = 0.5 * (recVRhoG(EigenAll, Seq123).array().square().rowwise().sum()) / recVRhoG(EigenAll, 0).array();
1804 Eigen::Vector<real, Eigen::Dynamic> eInternalS = recVRhoG(EigenAll, I4) - ek;
1805 real thetaP = 1.0;
1806
1807 if (pCent <= 2 * pEps)
1808 thetaP = 0;
1809 else
1810 for (int iG = 0; iG < rhoS.size(); iG++)
1811 {
1812 if (eInternalS(iG) < 2 * pEps / (gamma - 1))
1813 {
1814 real thetaThis = Gas::IdealGasGetCompressionRatioPressure<dim, 0, nVarsFixed>(
1815 recBase(iG, EigenAll).transpose(), recInc(iG, EigenAll).transpose(),
1816 pEps);
1817 thetaP = std::min(thetaP, thetaThis);
1818 }
1819 }
1820
1821 uRecBeta[iCell](0) = theta1 * thetaP;
1822 // if (uRecBeta[iCell](0) < 1)
1823 // uRecBeta[iCell](0) *= 1 - 1e-8;
1824 // uRecBeta[iCell](0) = std::min(theta1, thetaP);
1825 if (uRecBeta[iCell](0) < 1)
1826 {
1827 nLimLocal++;
1828 // uRecBeta[iCell](0) *= uRecBeta[iCell](0) < 0.99 ? 0. : 0.99; //! for safety
1829 uRecBeta[iCell](0) *= std::pow(uRecBeta[iCell](0), static_cast<int>(std::round(settings.uRecBetaCompressPower))) * safetyRatio;
1830 minBetaLocal = std::min(uRecBeta[iCell](0), minBetaLocal);
1831 if (uRecBeta[iCell](0) < betaCutOff)
1832 uRecBeta[iCell](0) = 0;
1833 }
1834 if (uRecBeta[iCell](0) < 0)
1835 {
1836 std::cout << fmt::format("theta1 {}, thetaP {}", theta1, thetaP) << std::endl;
1837 DNDS_assert(false);
1838 }
1839 if (uRecBeta[iCell](0) < 1)
1840 uRec[iCell] = (uRec[iCell] - uRecBase) * uRecBeta[iCell](0) + uRecBase;
1841
1842 // validation:
1843 recInc = quadBase * uRec[iCell];
1844 recVRhoG = recInc.rowwise() + u[iCell].transpose();
1845 ek = 0.5 * (recVRhoG(EigenAll, Seq123).array().square().rowwise().sum()) / recVRhoG(EigenAll, 0).array();
1846 eInternalS = (recVRhoG(EigenAll, I4) - ek);
1847 for (int iG = 0; iG < eInternalS.size(); iG++)
1848 {
1849 if (eInternalS(iG) < pEps)
1850 {
1851 std::cout << std::scientific;
1852 std::cout << eInternalS.transpose() << std::endl;
1853 std::cout << curOrder << std::endl;
1854 std::cout << fmt::format("{} {} {}", theta1, thetaP, uRecBeta[iCell](0)) << std::endl;
1855 std::cout << u[iCell] << std::endl;
1856 std::cout << recInc.transpose() << std::endl;
1857 DNDS_assert(false);
1858 }
1859 }
1860 }
1861 MPI::Allreduce(&nLimLocal, &nLim, 1, DNDS_MPI_INDEX, MPI_SUM, u.father->getMPI().comm);
1862 MPI::Allreduce(&minBetaLocal, &betaMin, 1, DNDS_MPI_REAL, MPI_MIN, u.father->getMPI().comm);
1863 }
1864
1865 template <EulerModel model>
1866 bool EulerEvaluator<model>::AssertMeanValuePP(
1867 ArrayDOFV<nVarsFixed> &u, bool panic)
1868 {
1869 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
1870 real rhoEps = smallReal * settings.refUPrim(0) * 1e-1;
1871 real pEps = smallReal * settings.refUPrim(I4) * 1e-1;
1872 if (settings.ppEpsIsRelaxed)
1873 rhoEps *= 0, pEps *= 0;
1874 bool ret{true};
1875#if defined(DNDS_DIST_MT_USE_OMP)
1876# pragma omp parallel for schedule(runtime)
1877#endif
1878 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
1879 {
1880 real gamma = settings.idealGasProperty.gamma;
1881 real alphaRho = 1;
1882 if (u[iCell](0) < rhoEps)
1883 {
1884 if (panic)
1885 DNDS_assert_info(
1886 false,
1887 fmt::format(
1888 "AssertMeanValuePP Failed on cell {} rho\n",
1889 iCell) +
1890 fmt::format(
1891 " eps={}, value={}",
1892 rhoEps, u[iCell](0)));
1893#if defined(DNDS_DIST_MT_USE_OMP)
1894# pragma omp critical
1895#endif
1896 ret = false;
1897 }
1898 real rhoEi = u[iCell](I4) - 0.5 * u[iCell](Seq123).squaredNorm() / u[iCell](0);
1899 if (rhoEi < pEps / (gamma - 1))
1900 {
1901 if (panic)
1902 DNDS_assert_info(
1903 false,
1904 fmt::format(
1905 "AssertMeanValuePP Failed on cell {} rhoEi\n",
1906 iCell) +
1907 fmt::format(
1908 " eps={}, value={}",
1909 pEps / (gamma - 1), rhoEi));
1910#if defined(DNDS_DIST_MT_USE_OMP)
1911# pragma omp critical
1912#endif
1913 ret = false;
1914 }
1915 }
1916
1917 return ret;
1918 }
1919
1920 template <EulerModel model>
1921 /** @brief Evaluate per-cell RHS scaling factor (alpha) for positivity-preserving time stepping.
1922 *
1923 * Determines the largest alpha in [0,1] such that u + alpha*res remains physically
1924 * realizable (positive density and pressure) at all quadrature points. Used for
1925 * under-relaxation to preserve positivity during explicit or implicit updates.
1926 *
1927 * @param u Conservative variable DOF array.
1928 * @param uRec Reconstruction coefficients.
1929 * @param uRecBeta Per-cell reconstruction limiter coefficient.
1930 * @param res RHS residual to be scaled.
1931 * @param cellRHSAlpha Per-cell alpha scaling factor (output).
1932 * @param nLim Total number of limited cells (output, MPI-reduced).
1933 * @param alphaMin Global minimum alpha (output, MPI-reduced).
1934 * @param relax Relaxation factor applied to alpha.
1935 * @param compress Compression mode for reconstruction part.
1936 * @param flag Evaluation mode flags (e.g., EvaluateCellRHSAlpha_MIN_ALL).
1937 */
1938 void EulerEvaluator<model>::EvaluateCellRHSAlpha(
1939 ArrayDOFV<nVarsFixed> &u,
1940 ArrayRECV<nVarsFixed> &uRec,
1941 ArrayDOFV<1> &uRecBeta,
1942 ArrayDOFV<nVarsFixed> &res,
1943 ArrayDOFV<1> &cellRHSAlpha, index &nLim, real &alphaMin,
1944 real relax, int compress,
1945 int flag)
1946 {
1947 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
1948 real rhoEps = smallReal * settings.refUPrim(0) * 1e-1;
1949 real pEps = smallReal * settings.refUPrim(I4) * 1e-1;
1950 static const real safetyRatio = 1 - 1e-5;
1951 static const real minRatio = 0.5;
1952
1953 if (settings.ppEpsIsRelaxed)
1954 {
1955 pEps *= 0, rhoEps *= 0;
1956 DNDS_assert_info(relax < 1, "Relaxed eps only for using relaxation in alpha");
1957 }
1958
1959 index nLimLocal = 0;
1960 real alphaMinLocal = 1;
1961#if defined(DNDS_DIST_MT_USE_OMP)
1962# pragma omp parallel for schedule(runtime) reduction(+ : nLimLocal) reduction(min : alphaMinLocal)
1963#endif
1964 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
1965 {
1966 real gamma = settings.idealGasProperty.gamma;
1967 real alphaRho = 1;
1968 TU inc = res[iCell];
1969 DNDS_assert(u[iCell](0) >= rhoEps);
1970 real relaxedRho = rhoEps * relax + (u[iCell](0)) * (1 - relax);
1971 if (inc(0) < 0) // not < rhoEps!!!
1972 alphaRho = std::min(1.0, (u[iCell](0) - relaxedRho) / (-inc(0) - smallReal * inc(0)));
1973 DNDS_assert(alphaRho >= 0 && alphaRho <= 1);
1974 if constexpr (Traits::hasSA)
1975 {
1976 // ** ! currently do not mass - fix for SA
1977 // static real v1Eps = smallReal * settings.refUPrim(I4 + 1);
1978 // if (inc(I4 + 1) < 0)
1979 // alphaRho = std::min(alphaRho,
1980 // (u[iCell](I4 + 1) - v1Eps) / (-inc(I4 + 1) - smallReal * inc(I4 + 1)));
1981 // // use exp down:
1982 // if (inc(I4 + 1) + u[iCell](I4 + 1) < v1Eps)
1983 // {
1984 // DNDS_assert(inc(I4 + 1) < 0);
1985 // real declineV = inc(I4 + 1) / (u[iCell](I4 + 1) + 1e-6);
1986 // real newu5 = u[iCell](I4 + 1) * std::exp(declineV);
1987 // alphaRho = std::min(alphaRho,
1988 // (u[iCell](I4 + 1) - newu5) / (-inc(I4 + 1) - smallReal * inc(I4 + 1)));
1989 // }
1990 }
1991 if constexpr (Traits::has2EQ)
1992 {
1993 // static real v1Eps = smallReal * settings.refUPrim(I4 + 1);
1994 // static real v2Eps = smallReal * settings.refUPrim(I4 + 2);
1995 // if (inc(I4 + 1) < 0)
1996 // alphaRho = std::min(alphaRho,
1997 // (u[iCell](I4 + 1) - v1Eps) / (-inc(I4 + 1) - smallReal * inc(I4 + 1)));
1998 // if (inc(I4 + 2) < 0)
1999 // alphaRho = std::min(alphaRho,
2000 // (u[iCell](I4 + 2) - v2Eps) / (-inc(I4 + 2) - smallReal * inc(I4 + 2)));
2001 }
2002
2003 inc *= alphaRho;
2004
2005 TU uNew = u[iCell] + inc;
2006 real pNew = (uNew(I4) - 0.5 * uNew(Seq123).squaredNorm() / uNew(0)) * (gamma - 1);
2007 real pOld = (u[iCell](I4) - 0.5 * u[iCell](Seq123).squaredNorm() / u[iCell](0)) * (gamma - 1);
2008 real relaxedP = pEps;
2009 if (pNew < pOld)
2010 relaxedP = pEps + (pOld - pEps) * (1 - relax);
2011
2012 real alphaP = 1;
2013 if (pNew < relaxedP)
2014 {
2015 // todo: use high order accurate (add control switch)
2016 real alphaC = Gas::IdealGasGetCompressionRatioPressure<dim, 0, nVarsFixed>(
2017 u[iCell], inc, relaxedP / (gamma - 1));
2018 alphaP = std::min(alphaP, alphaC);
2019 }
2020 cellRHSAlpha[iCell](0) = alphaRho * alphaP;
2021 // cellRHSAlpha[iCell](0) = std::min(alphaRho, alphaP);
2022 if (cellRHSAlpha[iCell](0) < 1)
2023 {
2024 cellRHSAlpha[iCell](0) = std::pow(cellRHSAlpha[iCell](0), compress * static_cast<int>(std::round(settings.uRecAlphaCompressPower)));
2025 nLimLocal++,
2026 cellRHSAlpha[iCell] *= safetyRatio,
2027 alphaMinLocal = std::min(alphaMinLocal, cellRHSAlpha[iCell](0));
2028 } //! for safety
2029 }
2030 MPI::Allreduce(&nLimLocal, &nLim, 1, DNDS_MPI_INDEX, MPI_SUM, u.father->getMPI().comm);
2031 MPI::Allreduce(&alphaMinLocal, &alphaMin, 1, DNDS_MPI_REAL, MPI_MIN, u.father->getMPI().comm);
2032 if (flag & EvaluateCellRHSAlpha_MIN_IF_NOT_ONE)
2033#if defined(DNDS_DIST_MT_USE_OMP)
2034# pragma omp parallel for schedule(static)
2035#endif
2036 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
2037 {
2038 if (cellRHSAlpha[iCell](0) < 1)
2039 cellRHSAlpha[iCell](0) = alphaMin;
2040 }
2041 if (flag & EvaluateCellRHSAlpha_MIN_ALL)
2042#if defined(DNDS_DIST_MT_USE_OMP)
2043# pragma omp parallel for schedule(static)
2044#endif
2045 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
2046 cellRHSAlpha[iCell](0) = alphaMin;
2047 }
2048
2049 template <EulerModel model>
2050 /** @brief Expand the per-cell alpha limiter to neighboring cells for smoothness.
2051 *
2052 * Cells adjacent to limited cells (alpha < 1) are also given reduced alpha
2053 * values to create a smooth transition region, preventing isolated sharp
2054 * limiting boundaries.
2055 *
2056 * @param u Conservative variable DOF array.
2057 * @param uRec Reconstruction coefficients.
2058 * @param uRecBeta Per-cell reconstruction limiter coefficient.
2059 * @param res RHS residual.
2060 * @param cellRHSAlpha Per-cell alpha (input/output, expanded to neighbors).
2061 * @param nLim Additional limited cells from expansion (output, accumulated).
2062 * @param alphaMin Global minimum alpha (input).
2063 */
2064 void EulerEvaluator<model>::EvaluateCellRHSAlphaExpansion(
2065 ArrayDOFV<nVarsFixed> &u,
2066 ArrayRECV<nVarsFixed> &uRec,
2067 ArrayDOFV<1> &uRecBeta,
2068 ArrayDOFV<nVarsFixed> &res,
2069 ArrayDOFV<1> &cellRHSAlpha, index &nLim, real alphaMin)
2070 {
2071 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
2072 real rhoEps = smallReal * settings.refUPrim(0) * 1e-1;
2073 real pEps = smallReal * settings.refUPrim(I4) * 1e-1;
2074
2075 auto cellIsHalfAlpha = [&](index iCell) -> bool // iCell should be internal
2076 {
2077 bool ret = false;
2078 if (cellRHSAlpha[iCell](0) == 1.0)
2079 {
2080 auto c2f = mesh->cell2face[iCell];
2081 for (int ic2f = 0; ic2f < c2f.size(); ic2f++)
2082 {
2083 index iCellOther = vfv->CellFaceOther(iCell, c2f[ic2f]);
2084 if (iCellOther != UnInitIndex)
2085 if (cellRHSAlpha[iCellOther](0) != 1.0)
2086 ret = true;
2087 }
2088 }
2089 return ret;
2090 };
2091
2092 auto cellAdjAlphaMin = [&](index iCell) -> real // iCell should be internal
2093 {
2094 real ret = 1;
2095 auto c2f = mesh->cell2face[iCell];
2096 for (int ic2f = 0; ic2f < c2f.size(); ic2f++)
2097 {
2098 index iCellOther = vfv->CellFaceOther(iCell, c2f[ic2f]);
2099 if (iCellOther != UnInitIndex)
2100 ret = std::min(ret, cellRHSAlpha[iCellOther](0));
2101 }
2102 return ret;
2103 };
2104
2105 // std::vector<index> InterCells;
2106
2107 // for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
2108 // if (cellIsHalfAlpha(iCell))
2109 // InterCells.emplace_back(iCell);
2110
2111 index nLimLocal = 0;
2112 index nLimAdd = 0;
2113 // for (index iCell : InterCells)
2114 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
2115 {
2116 real gamma = settings.idealGasProperty.gamma;
2117 TU inc = res[iCell];
2118
2119 TU uNew = u[iCell] + inc;
2120 real pNew = (uNew(I4) - 0.5 * uNew(Seq123).squaredNorm() / uNew(0)) * (gamma - 1);
2121
2122 if (pNew < pEps || uNew(0) < rhoEps)
2123 {
2124 // cellRHSAlpha[iCell](0) = cellAdjAlphaMin(iCell);
2125 // DNDS_assert(cellRHSAlpha[iCell](0) == alphaMin);
2126 cellRHSAlpha[iCell](0) = alphaMin;
2127 }
2128 if constexpr (Traits::hasSA)
2129 {
2130 // ** ! currently do not mass - fix for SA
2131 // static real v1Eps = smallReal * settings.refUPrim(I4 + 1);
2132 // if (uNew(I4 + 1) < v1Eps)
2133 // cellRHSAlpha[iCell](0) = alphaMin;
2134 }
2135 if constexpr (Traits::has2EQ)
2136 {
2137 // static real v1Eps = smallReal * settings.refUPrim(I4 + 1);
2138 // static real v2Eps = smallReal * settings.refUPrim(I4 + 2);
2139 // if (uNew(I4 + 1) < v1Eps)
2140 // cellRHSAlpha[iCell](0) = alphaMin;
2141 // if (uNew(I4 + 2) < v2Eps)
2142 // cellRHSAlpha[iCell](0) = alphaMin;
2143 }
2144 }
2145 MPI::Allreduce(&nLimLocal, &nLimAdd, 1, DNDS_MPI_INDEX, MPI_SUM, u.father->getMPI().comm);
2146 nLim += nLimAdd;
2147 }
2148
2149 template <EulerModel model>
2150 /** @brief Smooth the local pseudo-time step by taking the minimum of the cell value and its neighbor average.
2151 *
2152 * Prevents large time-step jumps between adjacent cells by capping each cell's
2153 * dTau at the weighted average of its neighbors' values.
2154 *
2155 * @param dTau Input local pseudo-time step per cell.
2156 * @param dTauNew Smoothed pseudo-time step (output).
2157 */
2158 void EulerEvaluator<model>::MinSmoothDTau(
2159 ArrayDOFV<1> &dTau, ArrayDOFV<1> &dTauNew)
2160 {
2161 real smootherCentWeight = 1;
2162#if defined(DNDS_DIST_MT_USE_OMP)
2163# pragma omp parallel for schedule(runtime)
2164#endif
2165 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
2166 {
2167 auto c2f = mesh->cell2face[iCell];
2168 real nAdj = 0.;
2169 real dTMean = 0.;
2170 for (index iFace : c2f)
2171 {
2172 index iCellOther = vfv->CellFaceOther(iCell, iFace);
2173 if (iCellOther != UnInitIndex)
2174 {
2175 nAdj += 1.;
2176 dTMean += dTau[iCellOther](0);
2177 }
2178 }
2179 dTMean += nAdj * smootherCentWeight * dTau[iCell](0);
2180 dTMean /= nAdj * (1 + smootherCentWeight);
2181 dTauNew[iCell](0) = std::min(dTau[iCell](0), dTMean);
2182 }
2183 }
2184
2185 template <EulerModel model>
2186 /** @brief Update boundary condition profiles from current solution for profile-anchored BCs.
2187 *
2188 * For BC zones with anchorOpt==2, integrates the solution at boundary faces into
2189 * a 1D radial profile recorder. Used for radial equilibrium pressure BCs in
2190 * turbomachinery applications.
2191 *
2192 * @param u Conservative variable DOF array.
2193 * @param uRec Reconstruction coefficients.
2194 */
2195 void EulerEvaluator<model>::updateBCProfiles(ArrayDOFV<nVarsFixed> &u, ArrayRECV<nVarsFixed> &uRec)
2196 {
2197 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
2198 for (Geom::t_index i = Geom::BC_ID_DEFAULT_MAX; i < pBCHandler->size(); i++) // init code, consider adding to ctor
2199 {
2200 if (pBCHandler->GetFlagFromIDSoft(i, "anchorOpt") != 2)
2201 continue;
2202 if (!profileRecorders.count(i))
2203 {
2204 real RMin = veryLargeReal;
2205 real RMax = -veryLargeReal;
2206 profileRecorders.emplace(std::make_pair(i, OneDimProfile<nVarsFixed>(mesh->getMPI())));
2207 for (index iBnd = 0; iBnd < mesh->NumBnd(); iBnd++)
2208 {
2209 index iFace = mesh->bnd2faceV.at(iBnd);
2210 if (iFace < 0) // remember that some iBnd do not have iFace (for periodic case)
2211 continue;
2212 auto f2c = mesh->face2cell[iFace];
2213 auto gFace = vfv->GetFaceQuad(iFace);
2214
2215 Geom::Elem::SummationNoOp noOp;
2216 auto faceBndID = mesh->GetFaceZone(iFace);
2217 auto faceBCType = pBCHandler->GetTypeFromID(faceBndID);
2218 if (faceBndID == i)
2219 {
2220 Geom::tSmallCoords coo;
2221 mesh->GetCoordsOnFace(iFace, coo);
2222 for (int ic = 0; ic < coo.cols(); ic++)
2223 {
2224 real r = settings.frameConstRotation.rVec(coo(EigenAll, ic)).norm();
2225 RMin = std::min(r, RMin);
2226 RMax = std::max(r, RMax);
2227 }
2228 }
2229 }
2230 MPI::AllreduceOneReal(RMin, MPI_MIN, mesh->getMPI());
2231 MPI::AllreduceOneReal(RMax, MPI_MAX, mesh->getMPI());
2232 auto vExtra = pBCHandler->GetValueExtraFromID(i);
2233 // * valueExtra[0] == nDiv
2234 // * valueExtra[1] == divMethod
2235 // * valueExtra[2] == d0
2236 // * valueExtra[3] == printInfo
2237 index nDiv = vExtra.size() >= 1 ? vExtra(0) : 10;
2238 index divMethod = vExtra.size() >= 2 ? vExtra(1) : 0; // TODO: implement other distributions
2239 real divd0 = vExtra.size() >= 3 ? vExtra(2) : veryLargeReal;
2240
2241 if (divMethod == 0)
2242 profileRecorders.at(i).GenerateUniform(std::max(nDiv, index(10)), nVars, RMin, RMax);
2243 else
2244 profileRecorders.at(i).GenerateTanh(std::max(nDiv, index(10)), nVars, RMin, RMax, divd0);
2245 }
2246 }
2247 for (auto &v : profileRecorders)
2248 v.second.SetZero();
2249 for (index iBnd = 0; iBnd < mesh->NumBnd(); iBnd++)
2250 {
2251 index iFace = mesh->bnd2faceV.at(iBnd);
2252 if (iFace < 0) // remember that some iBnd do not have iFace (for periodic case)
2253 continue;
2254 auto f2c = mesh->face2cell[iFace];
2255 auto gFace = vfv->GetFaceQuad(iFace);
2256
2257 Geom::Elem::SummationNoOp noOp;
2258 auto faceBndID = mesh->GetFaceZone(iFace);
2259 auto faceBCType = pBCHandler->GetTypeFromID(faceBndID);
2260
2261 if (pBCHandler->GetFlagFromIDSoft(faceBndID, "anchorOpt") != 2)
2262 continue;
2263
2264 real RMin = veryLargeReal;
2265 real RMax = -veryLargeReal;
2266
2267 Geom::tSmallCoords coo;
2268 mesh->GetCoordsOnFace(iFace, coo);
2269 for (int ic = 0; ic < coo.cols(); ic++)
2270 {
2271 real r = settings.frameConstRotation.rVec(coo(EigenAll, ic)).norm();
2272 RMin = std::min(r, RMin);
2273 RMax = std::max(r, RMax);
2274 }
2275 TU valIn = u[f2c[0]];
2276 valIn(Seq123) = settings.frameConstRotation.rtzFrame(vfv->GetFaceQuadraturePPhys(iFace, -1)).transpose()(Seq012, Seq012) * valIn(Seq123); // to rtz frame
2277#ifndef USE_ABS_VELO_IN_ROTATION
2278 valIn(2) += valIn(0) * settings.frameConstRotation.Omega() * settings.frameConstRotation.rVec(vfv->GetFaceQuadraturePPhys(iFace, -1)).norm(); // to static value
2279#endif
2280 // std::cout << valIn.transpose() << std::endl;
2281 // std::cout << RMin << " " << RMax << " " << vfv->GetFaceArea(iFace) << std::endl;
2282 profileRecorders.at(faceBndID).AddSimpleInterval(valIn, vfv->GetFaceArea(iFace), RMin, RMax);
2283 }
2284 for (auto &v : profileRecorders)
2285 v.second.Reduce();
2286 }
2287
2288 template <EulerModel model>
2289 /** @brief Compute radial equilibrium pressure distribution from BC profiles.
2290 *
2291 * Integrates the centrifugal pressure gradient (rho * vt^2 / r) along the
2292 * radial direction in boundary profiles to obtain the radial equilibrium
2293 * pressure increment for turbomachinery outlet BCs.
2294 */
2295 void EulerEvaluator<model>::updateBCProfilesPressureRadialEq()
2296 {
2297 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
2298 for (auto &v : profileRecorders)
2299 {
2300 if (pBCHandler->GetFlagFromIDSoft(v.first, "anchorOpt") == 2)
2301 {
2302 v.second.v.array().rowwise() /= (v.second.div.array() + verySmallReal);
2303 v.second.div.setConstant(1.);
2304 v.second.v(I4, 0) = 0;
2305 for (index i = 1; i < v.second.Size(); i++)
2306 {
2307 real vt0 = v.second.v(2, i - 1) / v.second.v(0, i - 1);
2308 real vt1 = v.second.v(2, i) / v.second.v(0, i);
2309 real l0 = v.second.Len(i - 1);
2310 real l1 = v.second.Len(i);
2311 real ldist = 0.5 * (l0 + l1);
2312 real vtm = (vt0 * l0 + vt1 * l1) / (l0 + l1);
2313 real rhom = (v.second.v(0, i - 1) * l0 + v.second.v(0, i) * l1) / (l0 + l1);
2314 real rc = v.second.nodes[i];
2315 v.second.v(I4, i) = v.second.v(I4, i - 1) + rhom * sqr(vtm) / rc * ldist;
2316 }
2317 if (mesh->getMPI().rank == 0)
2318 {
2319 // std::cout << "nodes";
2320 // for (auto vv : v.second.nodes)
2321 // std::cout << vv << " ";
2322 // std::cout << "\n";
2323 auto vExtra = pBCHandler->GetValueExtraFromID(v.first);
2324 int showMethod = vExtra.size() >= 4 ? vExtra(3) : 0;
2325 if (showMethod)
2326 log() << fmt::format("EulerEvaluator<model>::updateBCProfilesPressureRadialEq: p rise: [{:.3e}]", v.second.v(I4, EigenLast))
2327 << std::endl;
2328 }
2329 }
2330 }
2331 }
2332}
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
DNDS_check_throw_info
#define DNDS_check_throw_info(expr, info)
Same as DNDS_check_throw but attaches a user-supplied info message to the thrown std::runtime_error.
Definition Errors.hpp:96
EulerEvaluator.hpp
Core finite-volume evaluator for compressible Navier-Stokes / Euler equations.
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
ExprtkWrapper.hpp
Runtime mathematical expression compiler and evaluator wrapping exprtk.
HardEigen.hpp
Robust linear-algebra primitives for small / moderately-sized Eigen matrices that are more numericall...
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::ArrayDof::min
real min()
Global minimum across all entries (collective).
Definition ArrayDOF.hpp:313
DNDS::Euler::ArrayDOFV
MPI-parallel distributed array of per-cell Degrees-of-Freedom (conservative variable vectors).
Definition Euler.hpp:93
DNDS::Euler::ArrayDOFV::addTo
void addTo(const t_self &R, real r)
Scaled addition: this[i] += r * R[i] for every cell (axpy-style).
Definition Euler.hpp:153
DNDS::Euler::ArrayDOFV::setConstant
void setConstant(real R)
Set all DOF entries to a uniform scalar value.
Definition Euler.hpp:103
DNDS::Euler::ArrayDOFV::componentWiseNorm1
Eigen::Vector< real, nVarsFixed > componentWiseNorm1()
Compute the per-component global L1 norm.
Definition Euler.hpp:269
DNDS::Euler::ArrayGRADV
MPI-parallel distributed array of per-cell gradient matrices.
Definition Euler.hpp:643
DNDS::Euler::ArrayRECV
MPI-parallel distributed array of per-cell reconstruction coefficient matrices.
Definition Euler.hpp:401
DNDS::Euler::EulerEvaluator
Core finite-volume evaluator for compressible Navier-Stokes / Euler equations.
Definition EulerEvaluator.hpp:78
DNDS::Euler::EulerEvaluator::MinSmoothDTau
void MinSmoothDTau(ArrayDOFV< 1 > &dTau, ArrayDOFV< 1 > &dTauNew)
Smooth the local time step across neighboring cells.
Definition EulerEvaluator.hxx:2158
DNDS::Euler::EulerEvaluator::EvaluateRecNorm
void EvaluateRecNorm(Eigen::Vector< real, -1 > &res, ArrayDOFV< nVarsFixed > &u, ArrayRECV< nVarsFixed > &uRec, index P=1, bool compare=false, const tFCompareField &FCompareField=[](const Geom::tPoint &p, real t) { return TU::Zero();}, const tFCompareFieldWeight &FCompareFieldWeight=[](const Geom::tPoint &p, real t) { return 1.0;}, real t=0)
Compute the reconstruction error norm (optionally against an analytical field).
Definition EulerEvaluator.hxx:1474
DNDS::Euler::EulerEvaluator::AssertMeanValuePP
bool AssertMeanValuePP(ArrayDOFV< nVarsFixed > &u, bool panic)
Assert that all cell-mean values are physically realizable.
Definition EulerEvaluator.hxx:1866
DNDS::Euler::EulerEvaluator::UpdateLUSGSBackward
void UpdateLUSGSBackward(real alphaDiag, real t, ArrayDOFV< nVarsFixed > &rhs, ArrayDOFV< nVarsFixed > &u, ArrayDOFV< nVarsFixed > &uInc, JacobianDiagBlock< nVarsFixed > &JDiag, ArrayDOFV< nVarsFixed > &uIncNew)
Deprecated: use UpdateSGS with uIncIsZero=true instead. Kept for backward compatibility....
DNDS::Euler::EulerEvaluator::FixUMaxFilter
void FixUMaxFilter(ArrayDOFV< nVarsFixed > &u)
Clip extreme conserved-variable values to prevent overflow.
Definition EulerEvaluator.hxx:1318
DNDS::Euler::EulerEvaluator::UpdateSGS
void UpdateSGS(real alphaDiag, real t, ArrayDOFV< nVarsFixed > &rhs, ArrayDOFV< nVarsFixed > &u, ArrayDOFV< nVarsFixed > &uInc, ArrayDOFV< nVarsFixed > &uIncNew, JacobianDiagBlock< nVarsFixed > &JDiag, bool forward, bool gsUpdate, TU &sumInc, bool uIncIsZero=false)
Symmetric Gauss-Seidel update for the implicit linear system.
Definition EulerEvaluator.hxx:550
DNDS::Euler::EulerEvaluator::EvaluateURecBeta
void EvaluateURecBeta(ArrayDOFV< nVarsFixed > &u, ArrayRECV< nVarsFixed > &uRec, ArrayDOFV< 1 > &uRecBeta, index &nLim, real &betaMin, int flag)
Evaluate the positivity-preserving reconstruction limiter (beta).
DNDS::Euler::EulerEvaluator::MeanValuePrim2Cons
void MeanValuePrim2Cons(ArrayDOFV< nVarsFixed > &w, ArrayDOFV< nVarsFixed > &u)
Convert cell-mean primitive variables to conservative variables.
Definition EulerEvaluator.hxx:1368
DNDS::Euler::EulerEvaluator::tFCompareField
std::function< TU(const Geom::tPoint &, real)> tFCompareField
Callback type for analytical comparison field.
Definition EulerEvaluator.hpp:601
DNDS::Euler::EulerEvaluator::UpdateLUSGSForward
void UpdateLUSGSForward(real alphaDiag, real t, ArrayDOFV< nVarsFixed > &rhs, ArrayDOFV< nVarsFixed > &u, ArrayDOFV< nVarsFixed > &uInc, JacobianDiagBlock< nVarsFixed > &JDiag, ArrayDOFV< nVarsFixed > &uIncNew)
Deprecated: use UpdateSGS with uIncIsZero=true instead. Kept for backward compatibility....
DNDS::Euler::EulerEvaluator::EvaluateNorm
void EvaluateNorm(Eigen::Vector< real, -1 > &res, ArrayDOFV< nVarsFixed > &rhs, index P=1, bool volWise=false, bool average=false)
Compute the norm of the RHS residual vector.
Definition EulerEvaluator.hxx:1394
DNDS::Euler::EulerEvaluator::updateBCProfilesPressureRadialEq
void updateBCProfilesPressureRadialEq()
Update boundary profiles using radial-equilibrium pressure extrapolation.
Definition EulerEvaluator.hxx:2295
DNDS::Euler::EulerEvaluator::MeanValueCons2Prim
void MeanValueCons2Prim(ArrayDOFV< nVarsFixed > &u, ArrayDOFV< nVarsFixed > &w)
Convert cell-mean conservative variables to primitive variables.
Definition EulerEvaluator.hxx:1348
DNDS::Euler::EulerEvaluator::tFCompareFieldWeight
std::function< real(const Geom::tPoint &, real)> tFCompareFieldWeight
Callback type for comparison weighting function.
Definition EulerEvaluator.hpp:602
DNDS::Euler::EulerEvaluator::LUSGSMatrixSolveJacobianLU
void LUSGSMatrixSolveJacobianLU(real alphaDiag, real t, ArrayDOFV< nVarsFixed > &rhs, ArrayDOFV< nVarsFixed > &u, ArrayDOFV< nVarsFixed > &uInc, ArrayDOFV< nVarsFixed > &uIncNew, ArrayDOFV< nVarsFixed > &bBuf, JacobianDiagBlock< nVarsFixed > &JDiag, JacobianLocalLU< nVarsFixed > &jacLU, bool uIncIsZero, TU &sumInc)
Solve the implicit linear system using the pre-factored local LU.
DNDS::Euler::EulerEvaluator::LimiterUGrad
void LimiterUGrad(ArrayDOFV< nVarsFixed > &u, ArrayGRADV< nVarsFixed, gDim > &uGrad, ArrayGRADV< nVarsFixed, gDim > &uGradNew, uint64_t flags=LIMITER_UGRAD_No_Flags)
Apply slope limiter to the gradient field.
DNDS::Euler::EulerEvaluator::EvaluateCellRHSAlpha
void EvaluateCellRHSAlpha(ArrayDOFV< nVarsFixed > &u, ArrayRECV< nVarsFixed > &uRec, ArrayDOFV< 1 > &uRecBeta, ArrayDOFV< nVarsFixed > &res, ArrayDOFV< 1 > &cellRHSAlpha, index &nLim, real &alphaMin, real relax, int compress=1, int flag=0)
Compute the positivity-preserving RHS scaling factor (alpha) per cell.
Definition EulerEvaluator.hxx:1938
DNDS::Euler::EulerEvaluator::updateBCProfiles
void updateBCProfiles(ArrayDOFV< nVarsFixed > &u, ArrayRECV< nVarsFixed > &uRec)
Update boundary 1-D profiles from the current solution.
Definition EulerEvaluator.hxx:2195
DNDS::Euler::EulerEvaluator::TimeAverageAddition
void TimeAverageAddition(ArrayDOFV< nVarsFixed > &w, ArrayDOFV< nVarsFixed > &wAveraged, real dt, real &tCur)
Accumulate time-averaged primitive variables for unsteady statistics.
Definition EulerEvaluator.hxx:1335
DNDS::Euler::EulerEvaluator::EvaluateCellRHSAlphaExpansion
void EvaluateCellRHSAlphaExpansion(ArrayDOFV< nVarsFixed > &u, ArrayRECV< nVarsFixed > &uRec, ArrayDOFV< 1 > &uRecBeta, ArrayDOFV< nVarsFixed > &res, ArrayDOFV< 1 > &cellRHSAlpha, index &nLim, real alphaMin)
Expand a previously computed cellRHSAlpha toward 1 where safe.
Definition EulerEvaluator.hxx:2064
DNDS::Euler::EulerEvaluator::UpdateSGSWithRec
void UpdateSGSWithRec(real alphaDiag, real t, ArrayDOFV< nVarsFixed > &rhs, ArrayDOFV< nVarsFixed > &u, ArrayRECV< nVarsFixed > &uRec, ArrayDOFV< nVarsFixed > &uInc, ArrayRECV< nVarsFixed > &uRecInc, JacobianDiagBlock< nVarsFixed > &JDiag, bool forward, TU &sumInc)
SGS sweep coupled with reconstruction update.
Definition EulerEvaluator.hxx:713
DNDS::Euler::JacobianDiagBlock
Per-cell diagonal or block-diagonal Jacobian storage for implicit time stepping.
Definition EulerJacobian.hpp:43
DNDS::Euler::JacobianDiagBlock::MatVecLeftInvert
TU MatVecLeftInvert(index iCell, TV v)
Left-multiply a vector by the inverse Jacobian block for cell iCell.
Definition EulerJacobian.hpp:215
DNDS::Euler::JacobianDiagBlock::getValue
tComponent getValue(index iCell) const
Return the Jacobian as a full matrix for cell iCell (either mode).
Definition EulerJacobian.hpp:126
DNDS::Euler::JacobianDiagBlock::GetInvert
void GetInvert()
Compute and cache the inverse of every cell's Jacobian block.
Definition EulerJacobian.hpp:153
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::Gas::IdealGasThermal
void IdealGasThermal(real E, real rho, real vSqr, real gamma, real &p, real &asqr, real &H)
Thin wrapper delegating to IdealGas::IdealGasThermal.
Definition Gas.hpp:190
DNDS::Euler::SpecialFields::IsentropicVortexCent
auto IsentropicVortexCent(EulerEvaluator< model > &eval, const Geom::tPoint &x, real t, int cnVars)
Stationary isentropic vortex centered at the origin with zero background flow.
Definition SpecialFields.hpp:148
DNDS::Euler::SpecialFields::IsentropicVortex30
auto IsentropicVortex30(EulerEvaluator< model > &eval, const Geom::tPoint &x, real t, int cnVars)
Analytic isentropic vortex on a [-10,20] x [-10,20] periodic domain (period 30).
Definition SpecialFields.hpp:93
DNDS::Euler::SpecialFields::IsentropicVortex10
auto IsentropicVortex10(EulerEvaluator< model > &eval, const Geom::tPoint &x, real t, int cnVars, real chi)
Analytic isentropic vortex on a [0,10] x [0,10] periodic domain.
Definition SpecialFields.hpp:42
DNDS::Euler::RANS_KOSST
@ RANS_KOSST
Menter k-omega SST two-equation model.
Definition Euler.hpp:903
DNDS::Euler::RANS_KOWilcox
@ RANS_KOWilcox
Wilcox k-omega two-equation model.
Definition Euler.hpp:902
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
@ NS
Navier-Stokes, 2D geometry, 3D physics (5 vars).
Definition Euler.hpp:876
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::Euler::NS_3D
@ NS_3D
Navier-Stokes, 3D geometry, 3D physics (5 vars).
Definition Euler.hpp:879
DNDS::Euler::NS_EX
@ NS_EX
Extended NS, 2D geometry, dynamic nVars (Eigen::Dynamic).
Definition Euler.hpp:883
DNDS::Euler::NS_2D
@ NS_2D
NS with 2D physics (2 velocity components, 4 vars). The only model with dim=2.
Definition Euler.hpp:878
DNDS::Euler::BCWallIsothermal
@ BCWallIsothermal
Isothermal wall; requires wall temperature in the BC value vector.
Definition EulerBC.hpp:41
DNDS::Euler::BCWall
@ BCWall
No-slip viscous wall boundary.
Definition EulerBC.hpp:39
DNDS::Geom::t_index
int32_t t_index
Definition Geometric.hpp:6
DNDS::Geom::tPoint
Eigen::Vector3d tPoint
Definition Geometric.hpp:9
DNDS::Geom::tSmallCoords
Eigen::Matrix< real, 3, Eigen::Dynamic > tSmallCoords
Definition Geometric.hpp:50
DNDS::MPI::AllreduceOneReal
void AllreduceOneReal(real &v, MPI_Op op, const MPIInfo &mpi)
Single-scalar Allreduce helper for reals (in-place, count = 1).
Definition MPI.hpp:673
DNDS::MPI::Allreduce
MPI_int Allreduce(const void *sendbuf, void *recvbuf, MPI_int count, MPI_Datatype datatype, MPI_Op op, MPI_Comm comm)
Wrapper over MPI_Allreduce.
Definition MPI.cpp:203
DNDS::DNDS_MPI_INDEX
const MPI_Datatype DNDS_MPI_INDEX
MPI datatype matching index (= MPI_INT64_T).
Definition MPI.hpp:90
DNDS::UnInitIndex
DNDS_CONSTANT const index UnInitIndex
Sentinel "not initialised" index value (= INT64_MIN).
Definition Defines.hpp:176
DNDS::verySmallReal
DNDS_CONSTANT const real verySmallReal
Catch-all lower bound ("effectively zero").
Definition Defines.hpp:194
DNDS::rowsize
int32_t rowsize
Row-width / per-row element-count type (signed 32-bit).
Definition Defines.hpp:109
DNDS::pi
DNDS_CONSTANT const real pi
π in double precision (matches DNDS_E_PI macro).
Definition Defines.hpp:199
DNDS::MatrixXR
Eigen::Matrix< real, Eigen::Dynamic, Eigen::Dynamic > MatrixXR
Column-major dynamic Eigen matrix of reals (default layout).
Definition Defines.hpp:205
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::veryLargeReal
DNDS_CONSTANT const real veryLargeReal
Catch-all upper bound ("practically infinity") for physical scalars.
Definition Defines.hpp:190
DNDS::DNDS_MPI_REAL
const MPI_Datatype DNDS_MPI_REAL
MPI datatype matching real (= MPI_REAL8).
Definition MPI.hpp:92
DNDS::ArrayPair::father
ssp< TArray > father
Owned-side array (must be resized before ghost setup).
Definition ArrayPair.hpp:170
DNDS::ArrayPair::Size
index Size() const
Combined row count (father->Size() + son->Size()).
Definition ArrayPair.hpp:292
DNDS::Direct::LocalLUBase::Solve
void Solve(tVec &b, tVec &result)
Definition Direct.hpp:485
DNDS::Euler::JacobianLocalLU
Complete LU factorization of the cell-local Jacobian for GMRES preconditioning.
Definition EulerJacobian.hpp:273
DNDS::Euler::OneDimProfile
One-dimensional profile for inlet/outlet boundary data.
Definition EulerBC.hpp:768
DNDS::Euler::TU_P_Reduction
Definition EulerEvaluator.hxx:1443
DNDS::Euler::TU_P_Reduction::reduce
void reduce(const TU_P_Reduction< TU > &R)
Definition EulerEvaluator.hxx:1452
DNDS::Euler::TU_P_Reduction::P
index P
Definition EulerEvaluator.hxx:1445
DNDS::Euler::TU_P_Reduction::setP
void setP(index nP)
Definition EulerEvaluator.hxx:1456
DNDS::Euler::TU_P_Reduction::_v
TU _v
Definition EulerEvaluator.hxx:1444
DNDS::Euler::TU_P_Reduction::setZero
void setZero(int n)
Definition EulerEvaluator.hxx:1447
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
r
tVec r(NCells)
z
tVec z(NCells)
x
tVec x(NCells)
velo
Eigen::Vector3d velo
Definition test_GasThermo.cpp:282
res
auto res
Definition test_ODE.cpp:196
n
Eigen::Vector3d n(1.0, 0.0, 0.0)