DNDSR 0.2.1
Distributed Numeric Data Structure for CFV
Loading...
Searching...
No Matches
EulerEvaluator_EvaluateDt.hxx
Go to the documentation of this file.
1/** @file EulerEvaluator_EvaluateDt.hxx
2 * @brief Template implementations of EulerEvaluator methods for local time-step estimation,
3 * wall distance computation (CGAL AABB and Poisson-based), positivity-preserving
4 * limiters, face eigenvalue estimation, numerical flux evaluation, source terms,
5 * boundary value generation, and output variable registration.
6 */
7#pragma once
8#include "EulerEvaluator.hpp"
9#include "RANS_ke.hpp"
10#include "Solver/Linear.hpp"
11
12#define CGAL_DISABLE_ROUNDING_MATH_CHECK // for valgrind
13#include <CGAL/Simple_cartesian.h>
14#include <CGAL/AABB_tree.h>
15#include <CGAL/AABB_traits_3.h>
16#include <CGAL/AABB_triangle_primitive_3.h>
17#undef CGAL_DISABLE_ROUNDING_MATH_CHECK
18
19namespace DNDS::Euler
20{
21
22 static const auto model = NS_SA;
23 DNDS_SWITCH_INTELLISENSE(
24 template <EulerModel model>
25 ,
26 template <>)
27 /** @brief Collect wall boundary face triangulations for AABB-based wall distance computation.
28 *
29 * Iterates over wall and isothermal-wall boundary faces, decomposes them into triangles
30 * (handling Line, Tri, and Quad element types), and appends to the output vector.
31 * When useQuadPatches is true, uses quadrature-point-based sub-patches instead of
32 * vertex-based triangulation.
33 *
34 * @param useQuadPatches If true, create triangle patches from quadrature points.
35 * @param trianglesOut Collected triangles as 3x3 matrices (columns = vertices) (output).
36 */
37 void EulerEvaluator<model>::GetWallDist_CollectTriangles(
38 bool useQuadPatches,
39 std::vector<Eigen::Matrix<real, 3, 3>> &trianglesOut)
40 {
41 for (index iBnd = 0; iBnd < mesh->NumBnd(); iBnd++)
42 {
43 if (pBCHandler->GetTypeFromID(mesh->GetBndZone(iBnd)) == EulerBCType::BCWall ||
44 pBCHandler->GetTypeFromID(mesh->GetBndZone(iBnd)) == EulerBCType::BCWallIsothermal)
45 {
46 index iFace = mesh->bnd2faceV[iBnd];
47 auto elem = mesh->GetFaceElement(iFace);
48 auto quad = vfv->GetFaceQuad(iFace);
49 if (!useQuadPatches)
50 {
51 if (elem.type == Geom::Elem::ElemType::Line2 || elem.type == Geom::Elem::ElemType::Line3)
52 {
53 Geom::tSmallCoords coords;
54 mesh->GetCoordsOnFace(iFace, coords);
55 Eigen::Matrix<real, 3, 3> tri;
56 tri(EigenAll, 0) = coords(EigenAll, 0);
57 tri(EigenAll, 1) = coords(EigenAll, 1);
58 tri(EigenAll, 2) = coords(EigenAll, 1) + Geom::tPoint{0., 0., vfv->GetFaceArea(iFace)};
59 trianglesOut.push_back(tri);
60 }
61 else if (elem.type == Geom::Elem::ElemType::Tri3 || elem.type == Geom::Elem::ElemType::Tri6)
62 {
63 Geom::tSmallCoords coords;
64 mesh->GetCoordsOnFace(iFace, coords);
65 Eigen::Matrix<real, 3, 3> tri;
66 tri(EigenAll, 0) = coords(EigenAll, 0);
67 tri(EigenAll, 1) = coords(EigenAll, 1);
68 tri(EigenAll, 2) = coords(EigenAll, 2);
69 trianglesOut.push_back(tri);
70 }
71 else if (elem.type == Geom::Elem::ElemType::Quad4 || elem.type == Geom::Elem::ElemType::Quad9)
72 {
73 Geom::tSmallCoords coords;
74 mesh->GetCoordsOnFace(iFace, coords);
75 Eigen::Matrix<real, 3, 3> tri;
76 tri(EigenAll, 0) = coords(EigenAll, 0);
77 tri(EigenAll, 1) = coords(EigenAll, 1);
78 tri(EigenAll, 2) = coords(EigenAll, 2);
79 trianglesOut.push_back(tri);
80 tri(EigenAll, 0) = coords(EigenAll, 0);
81 tri(EigenAll, 1) = coords(EigenAll, 2);
82 tri(EigenAll, 2) = coords(EigenAll, 3);
83 trianglesOut.push_back(tri);
84 }
85 else
86 {
87 DNDS_assert_info(false, "This elem not implemented yet for GetWallDist()");
88 }
89 }
90 else
91 {
92 auto qPatches = Geom::Elem::GetQuadPatches(quad);
93 for (auto &qPatch : qPatches)
94 {
95 Eigen::Matrix<real, 3, 3> tri;
96 Geom::tSmallCoords coords;
97 mesh->GetCoordsOnFace(iFace, coords);
98 for (int iV = 0; iV < 3; iV++)
99 if (qPatch[iV] > 0)
100 tri(EigenAll, iV) = coords(EigenAll, qPatch[iV] - 1);
101 else if (qPatch[iV] < 0)
102 tri(EigenAll, iV) = vfv->GetFaceQuadraturePPhys(iFace, -qPatch[iV] - 1);
103 else
104 tri(EigenAll, iV) = coords(EigenAll, 1) + Geom::tPoint{0., 0., vfv->GetFaceArea(iFace)};
105 trianglesOut.push_back(tri);
106 }
107 }
108 }
109 }
110 }
111
112 DNDS_SWITCH_INTELLISENSE(
113 template <EulerModel model>
114 ,
115 template <>)
116 /** @brief Compute face-averaged wall distances from cell wall distances.
117 *
118 * For each face, averages the wall distance of its two adjacent cells (or uses
119 * the single owner cell for boundary faces) and stores in dWallFace.
120 */
121 void EulerEvaluator<model>::GetWallDist_ComputeFaceDistances()
122 {
123 dWallFace.resize(mesh->NumFaceProc());
124 for (index iFace = 0; iFace < mesh->NumFaceProc(); iFace++)
125 {
126 auto f2c = mesh->face2cell[iFace];
127 real facialDist = dWall.at(f2c[0]).mean();
128 if (f2c[1] != UnInitIndex)
129 facialDist = 0.5 * (facialDist + dWall.at(f2c[1]).mean());
130 dWallFace[iFace] = facialDist;
131 }
132 }
133
134 DNDS_SWITCH_INTELLISENSE(
135 template <EulerModel model>
136 ,
137 template <>)
138 /** @brief Compute geometric wall distance using CGAL AABB tree (all-to-all broadcast).
139 *
140 * Gathers all wall triangles globally, builds a CGAL AABB tree, and queries the
141 * nearest distance for each cell's quadrature points. Uses execution serialization
142 * (wallDistExection) to limit concurrent memory usage. Clamps result to minWallDist.
143 */
144 void EulerEvaluator<model>::GetWallDist_AABB()
145 {
146 using TriArray = ArrayEigenMatrix<3, 3>;
147 ssp<TriArray> TrianglesLocal, TrianglesFull;
148 DNDS_MAKE_SSP(TrianglesLocal, mesh->getMPI());
149 DNDS_MAKE_SSP(TrianglesFull, mesh->getMPI());
150
151 bool useQuadPatches = (settings.wallDistScheme == 1);
152 std::vector<Eigen::Matrix<real, 3, 3>> Triangles;
153 GetWallDist_CollectTriangles(useQuadPatches, Triangles);
154
155 TrianglesLocal->Resize(Triangles.size(), 3, 3);
156 for (index i = 0; i < TrianglesLocal->Size(); i++)
157 (*TrianglesLocal)[i] = Triangles[i];
158 Triangles.clear();
159 ArrayTransformerType<TriArray>::Type TrianglesTransformer;
160 TrianglesTransformer.setFatherSon(TrianglesLocal, TrianglesFull);
161 TrianglesTransformer.createFatherGlobalMapping();
162
163 std::vector<index> pullingSet;
164 pullingSet.resize(TrianglesTransformer.pLGlobalMapping->globalSize());
165 for (index i = 0; i < pullingSet.size(); i++)
166 pullingSet[i] = i;
167 TrianglesTransformer.createGhostMapping(pullingSet);
168 TrianglesTransformer.createMPITypes();
169 TrianglesTransformer.pullOnce();
170 if (mesh->coords.father->getMPI().rank == 0)
171 log() << fmt::format("=== EulerEvaluator<model>::GetWallDist() with minWallDist = {:.4e}, ", settings.minWallDist)
172 << " To search in " << TrianglesFull->Size() << std::endl;
173
174 auto executeSearch = [&]()
175 {
176 log() << fmt::format("Start search rank [{}]", mesh->getMPI().rank) << std::endl;
177 typedef CGAL::Simple_cartesian<double> K;
178 typedef K::FT FT;
179 typedef K::Point_3 Point;
180 typedef K::Triangle_3 Triangle;
181 typedef std::vector<Triangle>::iterator Iterator;
182 typedef CGAL::AABB_triangle_primitive_3<K, Iterator> Primitive;
183 typedef CGAL::AABB_traits_3<K, Primitive> AABB_triangle_traits;
184 typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
185
186 std::vector<Triangle> triangles;
187 triangles.reserve(TrianglesFull->Size());
188
189 for (index i = 0; i < TrianglesFull->Size(); i++)
190 {
191 Point p0((*TrianglesFull)[i](0, 0), (*TrianglesFull)[i](1, 0), (*TrianglesFull)[i](2, 0));
192 Point p1((*TrianglesFull)[i](0, 1), (*TrianglesFull)[i](1, 1), (*TrianglesFull)[i](2, 1));
193 Point p2((*TrianglesFull)[i](0, 2), (*TrianglesFull)[i](1, 2), (*TrianglesFull)[i](2, 2));
194 triangles.push_back(Triangle(p0, p1, p2));
195 }
196 TrianglesLocal->Resize(0, 3, 3);
197 TrianglesFull->Resize(0, 3, 3);
198 double minDist = veryLargeReal;
199 this->dWall.resize(mesh->NumCellProc());
200
201 if (!triangles.empty())
202 {
203 Tree tree(triangles.begin(), triangles.end());
204
205 for (index iCell = 0; iCell < mesh->NumCellProc(); iCell++)
206 {
207 auto quadCell = vfv->GetCellQuad(iCell);
208 dWall[iCell].resize(quadCell.GetNumPoints());
209 for (int ig = 0; ig < quadCell.GetNumPoints(); ig++)
210 {
211 auto p = vfv->GetCellQuadraturePPhys(iCell, ig);
212 Point pQ(p[0], p[1], p[2]);
213 FT sqd = tree.squared_distance(pQ);
214 dWall[iCell][ig] = std::sqrt(sqd);
215 if (dWall[iCell][ig] < minDist)
216 minDist = dWall[iCell][ig];
217 dWall[iCell][ig] = std::max(settings.minWallDist, dWall[iCell][ig]);
218 }
219 }
220 }
221 else
222 {
223 for (index iCell = 0; iCell < mesh->NumCellProc(); iCell++)
224 {
225 auto quadCell = vfv->GetCellQuad(iCell);
226 dWall[iCell].setConstant(quadCell.GetNumPoints(), std::pow(veryLargeReal, 1. / 4.));
227 }
228 }
229 log() << fmt::format("[{}] MinDist: ", mesh->getMPI().rank) << minDist << std::endl;
230 };
231 if (settings.wallDistExection == 1)
232 MPISerialDo(mesh->getMPI(), [&]()
233 { executeSearch(); });
234 else if (settings.wallDistExection > 1)
235 for (int i = 0; i < settings.wallDistExection; i++)
236 {
237 if (mesh->getMPI().rank % settings.wallDistExection == i)
238 executeSearch();
239 MPI::Barrier(mesh->getMPI().comm);
240 }
241 else
242 executeSearch();
243 }
244
245 DNDS_SWITCH_INTELLISENSE(
246 template <EulerModel model>
247 ,
248 template <>)
249 /** @brief Compute geometric wall distance using batched CGAL AABB tree queries.
250 *
251 * Builds a local AABB tree from wall triangles (using quadrature patches), then
252 * performs batched distributed queries where cells are processed in configurable
253 * load-sized chunks. Supports refinement mode (wallDistScheme==20) for selective
254 * re-querying of cells below wallDistRefineMax.
255 */
256 void EulerEvaluator<model>::GetWallDist_BatchedAABB()
257 {
258 typedef CGAL::Simple_cartesian<double> K;
259 typedef K::FT FT;
260 typedef K::Point_3 Point;
261 typedef K::Triangle_3 Triangle;
262 typedef std::vector<Triangle>::iterator Iterator;
263 typedef CGAL::AABB_triangle_primitive_3<K, Iterator> Primitive;
264 typedef CGAL::AABB_traits_3<K, Primitive> AABB_triangle_traits;
265 typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
266
267 // Build local AABB tree from wall boundary triangles using quad patches.
268 std::vector<Triangle> triangles;
269 triangles.reserve(mesh->NumBnd() * 8 + 8);
270 {
271 std::vector<Eigen::Matrix<real, 3, 3>> triEigen;
272 GetWallDist_CollectTriangles(/*useQuadPatches=*/true, triEigen);
273 for (auto &tri : triEigen)
274 {
275 Point p0(tri(0, 0), tri(1, 0), tri(2, 0));
276 Point p1(tri(0, 1), tri(1, 1), tri(2, 1));
277 Point p2(tri(0, 2), tri(1, 2), tri(2, 2));
278 triangles.push_back(Triangle(p0, p1, p2));
279 }
280 }
281
282 if (triangles.empty())
283 {
284 triangles.push_back(Triangle(
285 Point(veryLargeReal, veryLargeReal, veryLargeReal),
286 Point(veryLargeReal, veryLargeReal, veryLargeReal),
287 Point(veryLargeReal, veryLargeReal, veryLargeReal)));
288 }
289 MPISerialDo(mesh->getMPI(),
290 [&]()
291 {
292 log() << fmt::format("[{},{}] ", mesh->getMPI().rank, triangles.size());
293 if (mesh->getMPI().rank % 10 == 0)
294 log() << "\n";
295 log().flush();
296 });
297 if (mesh->getMPI().rank == 0)
298 log() << std::endl;
299 Tree tree(triangles.begin(), triangles.end());
300
301 if (settings.wallDistScheme == 2)
302 dWall.resize(mesh->NumCellProc());
303 index iCellBase = 0;
304 index nProcessed = 0;
305 int cellLoadNum = std::max(1, static_cast<int>(std::ceil(settings.wallDistCellLoadSize / real(mesh->getMPI().size))));
306 if (mesh->coords.father->getMPI().rank == 0)
307 log() << fmt::format("=== EulerEvaluator<model>::GetWallDist() using cellLoadNum = {}, ", cellLoadNum)
308 << std::endl;
309
310 if (settings.wallDistScheme == 20)
311 {
312 index nRefine = 0;
313 for (auto &ds : dWall)
314 for (auto d : ds)
315 if (d <= settings.wallDistRefineMax)
316 nRefine++;
317 MPI::AllreduceOneIndex(nRefine, MPI_SUM, mesh->getMPI());
318 if (mesh->coords.father->getMPI().rank == 0)
319 log() << fmt::format("=== EulerEvaluator<model>::GetWallDist() to refine {}, ", nRefine)
320 << std::endl;
321 }
322
323 real t0 = MPI_Wtime();
324 for (index iIter = 0;; iIter++)
325 {
326 std::vector<Geom::tPoint> pnts;
327 pnts.reserve(cellLoadNum * 64);
328 for (int iCLoad = 0; iCLoad < cellLoadNum; iCLoad++)
329 {
330 index iCell = iCellBase + iCLoad;
331 if (iCell < mesh->NumCellProc())
332 {
333 auto quadCell = vfv->GetCellQuad(iCell);
334 for (int ig = 0; ig < quadCell.GetNumPoints(); ig++)
335 {
336 if (settings.wallDistScheme == 20)
337 if (dWall[iCell][ig] > settings.wallDistRefineMax)
338 continue;
339 auto p = vfv->GetCellQuadraturePPhys(iCell, ig);
340 pnts.push_back(p);
341 }
342 }
343 }
344
345 using PntArray = ArrayEigenMatrix<3, 1>;
346 ssp<PntArray> PntArrayLocal, PntArrayFull;
347 DNDS_MAKE_SSP(PntArrayLocal, mesh->getMPI());
348 DNDS_MAKE_SSP(PntArrayFull, mesh->getMPI());
349 PntArrayLocal->Resize(pnts.size(), 3, 1);
350 for (size_t i = 0; i < pnts.size(); i++)
351 (*PntArrayLocal)[i] = pnts[i];
352 ArrayTransformerType<PntArray>::Type PntTransformer;
353 PntTransformer.setFatherSon(PntArrayLocal, PntArrayFull);
354 PntTransformer.createFatherGlobalMapping();
355 std::vector<index> pullingSet;
356 pullingSet.resize(PntTransformer.pLGlobalMapping->globalSize());
357 if (!pullingSet.size())
358 break;
359 if (mesh->getMPI().rank == 0)
360 log() << fmt::format("=== EulerEvaluator<model>::GetWallDist() iter [{}], nProcessed [{}], t [{:.6g}] ",
361 iIter, nProcessed, MPI_Wtime() - t0)
362 << std::endl;
363 for (index i = 0; i < pullingSet.size(); i++)
364 pullingSet[i] = i;
365 PntTransformer.createGhostMapping(pullingSet);
366 PntTransformer.createMPITypes();
367 PntTransformer.pullOnce();
368 if (mesh->getMPI().rank == 0)
369 log() << fmt::format("=== EulerEvaluator<model>::GetWallDist() iter [{}], pullOnce done, t [{:.6g}] ",
370 iIter, MPI_Wtime() - t0)
371 << std::endl;
372
373 std::vector<real> distQueryFull(PntArrayFull->Size(), veryLargeReal);
374 for (int iQ = 0; iQ < PntArrayFull->Size(); iQ++)
375 {
376 Point pQ((*PntArrayFull)[iQ][0], (*PntArrayFull)[iQ][1], (*PntArrayFull)[iQ][2]);
377 FT sqd = tree.squared_distance(pQ);
378 distQueryFull[iQ] = std::sqrt(sqd);
379 }
380 if (mesh->getMPI().rank == 0)
381 log() << fmt::format("=== EulerEvaluator<model>::GetWallDist() iter [{}], query done, t [{:.6g}] ",
382 iIter, MPI_Wtime() - t0)
383 << std::endl;
384
385 std::vector<real> distQueryFullReduced(PntArrayFull->Size(), veryLargeReal);
386 MPI::Allreduce(distQueryFull.data(), distQueryFullReduced.data(), distQueryFull.size(), DNDS_MPI_REAL, MPI_MIN, mesh->getMPI().comm);
387
388 if (mesh->getMPI().rank == 0)
389 log() << fmt::format("=== EulerEvaluator<model>::GetWallDist() iter [{}], reduce done, t [{:.6g}] ",
390 iIter, MPI_Wtime() - t0)
391 << std::endl;
392 index iQLoad = 0;
393 for (int iCLoad = 0; iCLoad < cellLoadNum; iCLoad++)
394 {
395 index iCell = iCellBase + iCLoad;
396 if (iCell < mesh->NumCellProc())
397 {
398 auto quadCell = vfv->GetCellQuad(iCell);
399 if (settings.wallDistScheme == 2)
400 dWall[iCell].resize(quadCell.GetNumPoints());
401 for (int ig = 0; ig < quadCell.GetNumPoints(); ig++)
402 {
403 if (settings.wallDistScheme == 20)
404 if (dWall[iCell][ig] > settings.wallDistRefineMax)
405 continue;
406 dWall[iCell][ig] = distQueryFullReduced.at(
407 PntTransformer.pLGhostMapping->ghostStart.at(mesh->getMPI().rank) +
408 iQLoad);
409 iQLoad++;
410 }
411 }
412 }
413
414 iCellBase += cellLoadNum;
415 nProcessed += pullingSet.size();
416 }
417 }
418
419 DNDS_SWITCH_INTELLISENSE(
420 template <EulerModel model>
421 ,
422 template <>)
423 /** @brief Compute wall distance by solving a Poisson equation with pseudo-time stepping.
424 *
425 * Solves a p-Poisson equation (grad^2 phi = -1) with Dirichlet BC (phi=0) on
426 * wall boundaries and Neumann BC elsewhere, using Jacobi or GMRES iterations.
427 * The wall distance is recovered from the solution as d = |grad(phi)| + sqrt(|grad(phi)|^2 + 2*phi).
428 * Configurable via settings: wallDistNJacobiSweep, wallDistIter, wallDistPoissonP, etc.
429 */
430 void EulerEvaluator<model>::GetWallDist_Poisson()
431 {
432 int nSweep = settings.wallDistNJacobiSweep;
433 int nIter = settings.wallDistIter;
434 int nIterSee = 10;
435 real smoothBias = 0.0;
436 int nGmresSubspace = 10;
437 int nGmresRestart = 1;
438 int useGmres = settings.wallDistLinSolver;
439 real resTh = settings.wallDistResTol;
440
441 int nPPoissonStartIter = settings.wallDistIterStart;
442 int nPPoisson = settings.wallDistPoissonP;
443
444 real dTauScale = settings.wallDistDTauScale;
445
446 ArrayDOFV<1> phi, rPhi, dPhi, dPhiNew, phiTmp, dTauInv;
447 ArrayDOFV<3> diffPhi;
448 vfv->BuildUDof(phi, 1);
449 vfv->BuildUDof(rPhi, 1);
450 vfv->BuildUDof(dPhi, 1);
451 vfv->BuildUDof(dPhiNew, 1);
452 vfv->BuildUDof(phiTmp, 1);
453 vfv->BuildUDof(dTauInv, 1);
454 vfv->BuildUDof(diffPhi, 3);
455 phi.setConstant(0.0);
456 dTauInv.setConstant(0.0);
457
458 std::vector<std::vector<real>> coefs;
459 coefs.resize(mesh->NumCell());
460 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
461 coefs.at(iCell).resize(mesh->cell2face[iCell].size() + 1);
462
463 std::vector<MatrixXR> mGGs;
464 mGGs.resize(mesh->NumCell());
465
466 // get mGG coefs
467 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
468 {
469 Geom::tPoint bary = vfv->GetCellQuadraturePPhys(iCell, -1);
470 auto &mGG = mGGs.at(iCell);
471 mGG.setZero(3 + mesh->cell2face[iCell].size(), 3 + mesh->cell2face[iCell].size());
472 mGG({0, 1, 2}, {0, 1, 2}) = Eigen::Matrix3d::Identity();
473 real sumFaceArea = 0.;
474 for (int ic2f = 0; ic2f < mesh->cell2face[iCell].size(); ic2f++)
475 {
476 index iFace = mesh->cell2face[iCell][ic2f];
477 index iCellOther = mesh->CellFaceOther(iCell, iFace);
478 int if2c = mesh->CellIsFaceBack(iCell, iFace) ? 0 : 1;
479 Geom::tPoint uNormOut = vfv->GetFaceNormFromCell(iFace, iCell, if2c, -1) * (if2c ? -1 : 1);
480 auto faceBndID = mesh->GetFaceZone(iFace);
481 auto faceBCType = pBCHandler->GetTypeFromID(faceBndID);
482 Geom::tPoint baryOther = bary;
483 Geom::tPoint bFace = vfv->GetFaceQuadraturePPhysFromCell(iFace, iCell, -1, -1);
484 if (iCellOther != UnInitIndex)
485 {
486 baryOther = vfv->GetOtherCellPointFromCell(
487 iCell, iCellOther, iFace,
488 vfv->GetCellQuadraturePPhys(iCellOther, -1));
489 }
490 else
491 {
492 DNDS_assert(faceBCType != BCUnknown);
493 DNDS_assert(!Geom::FaceIDIsInternal(faceBndID));
494 baryOther = bFace * 2 - bary;
495 }
496 mGG({0, 1, 2}, 3 + ic2f) = -uNormOut * vfv->GetFaceArea(iFace) / vfv->GetCellVol(iCell);
497 mGG(3 + ic2f, {0, 1, 2}) = (baryOther + bary - 2 * bFace).transpose();
498 mGG(3 + ic2f, 3 + ic2f) = 2;
499 sumFaceArea += vfv->GetFaceArea(iFace);
500 }
501 diffPhi[iCell] /= vfv->GetCellVol(iCell);
502 auto mGGLU = mGG.fullPivLu();
503 DNDS_assert_info(mGGLU.isInvertible(), fmt::format("[[{}]]\n", Eigen::MatrixFMTSafe<real, -1, -1>(mGG)));
504 MatrixXR mGGInv = mGGLU.inverse();
505 mGG = mGGInv;
506
507 real LCell = vfv->GetCellVol(iCell) / sumFaceArea;
508 dTauInv[iCell](0) = 1. / (1. / std::pow(LCell, 2) * dTauScale);
509 }
510 dTauInv.trans.startPersistentPull();
511 dTauInv.trans.waitPersistentPull();
512
513 auto getDiffPhi = [&](ArrayDOFV<1> &phi, ArrayDOFV<3> &diffPhi)
514 {
515 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
516 {
517 Geom::tPoint bary = vfv->GetCellQuadraturePPhys(iCell, -1);
518 diffPhi[iCell].setZero();
519 Eigen::Vector<real, Eigen::Dynamic> bGG;
520 bGG.setZero(3 + mesh->cell2face[iCell].size());
521 for (int ic2f = 0; ic2f < mesh->cell2face[iCell].size(); ic2f++)
522 {
523 index iFace = mesh->cell2face[iCell][ic2f];
524 index iCellOther = mesh->CellFaceOther(iCell, iFace);
525 int if2c = mesh->CellIsFaceBack(iCell, iFace) ? 0 : 1;
526 Geom::tPoint uNormOut = vfv->GetFaceNormFromCell(iFace, iCell, if2c, -1) * (if2c ? -1 : 1);
527 auto faceBndID = mesh->GetFaceZone(iFace);
528 auto faceBCType = pBCHandler->GetTypeFromID(faceBndID);
529 real phiOther = phi[iCell](0);
530 Geom::tPoint baryOther = bary;
531 Geom::tPoint bFace = vfv->GetFaceQuadraturePPhysFromCell(iFace, iCell, -1, -1);
532 if (iCellOther != UnInitIndex)
533 {
534 phiOther = phi[iCellOther](0);
535 baryOther = vfv->GetOtherCellPointFromCell(
536 iCell, iCellOther, iFace,
537 vfv->GetCellQuadraturePPhys(iCellOther, -1));
538 }
539 else
540 {
541 DNDS_assert(faceBCType != BCUnknown);
542 DNDS_assert(!Geom::FaceIDIsInternal(faceBndID));
543 if (faceBCType == BCWall || faceBCType == BCWallIsothermal)
544 phiOther = -phi[iCell](0);
545 baryOther = bFace * 2 - bary;
546 }
547 real distThis = (bFace - bary).norm();
548 real distOther = (bFace - baryOther).norm();
549 real phiFace = (distOther * phi[iCell](0) + distThis * phiOther) / (distOther + distThis) - phi[iCell](0);
550 diffPhi[iCell] += uNormOut * vfv->GetFaceArea(iFace) * phiFace;
551 bGG(3 + ic2f) = phiOther + phi[iCell](0);
552 }
553 diffPhi[iCell] /= vfv->GetCellVol(iCell);
554 diffPhi[iCell] = (mGGs.at(iCell) * bGG)({0, 1, 2}); // comment to use traditional GG
555 }
556 };
557
558 int pPoissonCur = 2;
559
560 auto rhsPhi = [&](ArrayDOFV<1> &phi, ArrayDOFV<3> &diffPhi, ArrayDOFV<1> &rhs, std::vector<std::vector<real>> &coefs, bool updateCoefs)
561 {
562 const real supressRec = 1.0;
563 getDiffPhi(phi, diffPhi);
564 diffPhi.trans.startPersistentPull();
565 diffPhi.trans.waitPersistentPull();
566 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
567 {
568 rhs[iCell](0) = 1.;
569 Geom::tPoint bary = vfv->GetCellQuadraturePPhys(iCell, -1);
570 if (updateCoefs)
571 coefs.at(iCell).at(0) = 0.;
572 for (int ic2f = 0; ic2f < mesh->cell2face[iCell].size(); ic2f++)
573 {
574 index iFace = mesh->cell2face[iCell][ic2f];
575 index iCellOther = mesh->CellFaceOther(iCell, iFace);
576 int if2c = mesh->CellIsFaceBack(iCell, iFace) ? 0 : 1;
577 Geom::tPoint uNormOut = vfv->GetFaceNormFromCell(iFace, iCell, if2c, -1) * (if2c ? -1 : 1);
578 auto faceBndID = mesh->GetFaceZone(iFace);
579 auto faceBCType = pBCHandler->GetTypeFromID(faceBndID);
580 real phiOther = phi[iCell](0);
581 Geom::tPoint baryOther = bary;
582 Geom::tPoint bFace = vfv->GetFaceQuadraturePPhysFromCell(iFace, iCell, -1, -1);
583 real phiThisFace = phi[iCell](0) + (bFace - bary).dot(diffPhi[iCell]) * supressRec;
584 real phiOtherFace = phiThisFace;
585 real diffPhiNormThis = diffPhi[iCell].dot(uNormOut) * supressRec;
586 real diffPhiNormOther = diffPhiNormThis;
587 Geom::tPoint diffPhiFace = diffPhi[iCell] * supressRec;
588 if (iCellOther != UnInitIndex)
589 {
590 phiOther = phi[iCellOther](0);
591 baryOther = vfv->GetOtherCellPointFromCell(
592 iCell, iCellOther, iFace,
593 vfv->GetCellQuadraturePPhys(iCellOther, -1));
594 phiOtherFace = phiOther + (bFace - baryOther).dot(diffPhi[iCellOther]) * supressRec;
595 diffPhiNormOther = diffPhi[iCellOther].dot(uNormOut) * supressRec; //! todo: periodic!!
596 diffPhiFace = 0.5 * (diffPhiFace + diffPhi[iCellOther] * supressRec);
597 }
598 else
599 {
600 DNDS_assert(faceBCType != BCUnknown);
601 DNDS_assert(!Geom::FaceIDIsInternal(faceBndID));
602 if (faceBCType == BCWall || faceBCType == BCWallIsothermal)
603 phiOther = -phi[iCell](0), phiOtherFace = -phiThisFace;
604 else
605 {
606 diffPhiNormOther = -diffPhiNormThis;
607 diffPhiFace.setZero();
608 }
609 baryOther = bFace * 2 - bary;
610 }
611 // real dist = (baryOther - bary).norm();
612 real dist = std::abs((baryOther - bary).dot(uNormOut));
613
614 real diffPhiNorm = ((phiOtherFace - phiThisFace) * 1.0 / dist + 0.5 * (diffPhiNormOther + diffPhiNormThis));
615 diffPhiFace += ((phiOtherFace - phiThisFace) * 1.0 / dist) * uNormOut;
616 real diffPhiFaceMag = diffPhiFace.norm();
617
618 // rhs[iCell](0) += (phiOther - phi[iCell](0)) * 1.0 / dist * vfv->GetFaceArea(iFace) / vfv->GetCellVol(iCell);
619 rhs[iCell](0) += std::pow(diffPhiFaceMag, pPoissonCur - 2) * diffPhiNorm * vfv->GetFaceArea(iFace) / vfv->GetCellVol(iCell);
620 if (updateCoefs)
621 {
622 coefs.at(iCell).at(0) -= 1.0 / dist * vfv->GetFaceArea(iFace) / vfv->GetCellVol(iCell) * std::pow(diffPhiFaceMag, pPoissonCur - 2) * (pPoissonCur - 1);
623 coefs.at(iCell).at(1 + ic2f) = 1.0 / dist * vfv->GetFaceArea(iFace) / vfv->GetCellVol(iCell) * std::pow(diffPhiFaceMag, pPoissonCur - 2) * (pPoissonCur - 1);
624 }
625 }
626 }
627 };
628
629 auto solveDphi = [&](ArrayDOFV<1> &rhsPhi, ArrayDOFV<1> &dphi, ArrayDOFV<1> &dphiNew, std::vector<std::vector<real>> &coefs, ArrayDOFV<1> &dTauInv)
630 {
631 dphi.setConstant(0.0);
632 for (int iSweep = 1; iSweep <= nSweep; iSweep++)
633 {
634 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
635 {
636 dphiNew[iCell] = rhsPhi[iCell];
637 for (int ic2f = 0; ic2f < mesh->cell2face[iCell].size(); ic2f++)
638 {
639 index iFace = mesh->cell2face[iCell][ic2f];
640 index iCellOther = mesh->CellFaceOther(iCell, iFace);
641 if (iCellOther != UnInitIndex)
642 dphiNew[iCell] += coefs[iCell][ic2f + 1] * dphi[iCellOther];
643 }
644 dphiNew[iCell] /= -coefs[iCell][0] + dTauInv[iCell](0);
645 }
646 dphiNew.trans.startPersistentPull();
647 dphiNew.trans.waitPersistentPull();
648 dphi = dphiNew;
649 }
650 };
651 Linear::GMRES_LeftPreconditioned<ArrayDOFV<1>> gmres(
652 nGmresSubspace,
653 [&](ArrayDOFV<1> &v)
654 {
655 vfv->BuildUDof(v, 1);
656 });
657 for (; pPoissonCur <= nPPoisson; pPoissonCur += 2)
658 {
659
660 real incNormBase{0};
661 real resNormBase{0};
662 for (int iIter = 1; iIter <= nIter; iIter++)
663 {
664 rhsPhi(phi, diffPhi, rPhi, coefs, true);
665 real normR = rPhi.norm2();
666 real normPhi = phi.norm2();
667 dPhi.setConstant(0.0);
668
669 if (useGmres == 0)
670 {
671 solveDphi(rPhi, dPhi, dPhiNew, coefs, dTauInv);
672 }
673 else
674 {
675 gmres.solve(
676 [&](ArrayDOFV<1> &x, ArrayDOFV<1> &Ax)
677 {
678 real normX = x.norm2();
679 real ratio = normPhi * 1e-7 / (normX + normPhi * 1e-7 + verySmallReal) + 1e-12;
680 dPhiNew = phi;
681 rhsPhi(dPhiNew, diffPhi, phiTmp, coefs, false);
682 dPhiNew.addTo(x, ratio);
683 rhsPhi(dPhiNew, diffPhi, Ax, coefs, false);
684 Ax.addTo(phiTmp, -1.0);
685 Ax *= -1. / ratio;
686 phiTmp = x;
687 phiTmp *= dTauInv;
688 Ax += phiTmp;
689 Ax.trans.startPersistentPull();
690 Ax.trans.waitPersistentPull();
691
692 // for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
693 // {
694 // Ax[iCell] = -coefs[iCell][0] * x[iCell];
695 // for (int ic2f = 0; ic2f < mesh->cell2face[iCell].size(); ic2f++)
696 // {
697 // index iFace = mesh->cell2face[iCell][ic2f];
698 // index iCellOther = mesh->CellFaceOther(iCell, iFace);
699 // if (iCellOther != UnInitIndex)
700 // Ax[iCell] += -coefs[iCell][ic2f + 1] * x[iCellOther];
701 // }
702 // }
703 // Ax.trans.startPersistentPull();
704 // Ax.trans.waitPersistentPull();
705 },
706 [&](ArrayDOFV<1> &x, ArrayDOFV<1> &Mx)
707 {
708 Mx.setConstant(0.0);
709 solveDphi(x, Mx, dPhiNew, coefs, dTauInv);
710 // Mx.trans.startPersistentPull();
711 // Mx.trans.waitPersistentPull();
712 },
713 [&](ArrayDOFV<1> &a, ArrayDOFV<1> &b) -> real
714 {
715 return a.dot(b);
716 },
717 rPhi,
718 dPhi,
719 nGmresRestart,
720 [&](int iRestart, real res, real resB)
721 {
722 // if (res < resB * 1e-6)
723 // return true;
724 return false;
725 });
726 }
727
728 real incNorm = dPhi.norm2();
729 DNDS_assert(std::isfinite(incNorm));
730 phi += dPhi;
731 phi.trans.startPersistentPull();
732 phi.trans.waitPersistentPull();
733 if (iIter == 1)
734 incNormBase = incNorm;
735 resNormBase = std::max(resNormBase, normR);
736 bool nowExit = false;
737 if ((normR < resNormBase * resTh * std::pow(0.1, pPoissonCur - 2)) && iIter > nPPoissonStartIter)
738 nowExit = true;
739 if (iIter % nIterSee == 0 || iIter == nIter || nowExit)
740 if (phi.father->getMPI().rank == 0)
741 log() << fmt::format("EulerEvaluator<model>::GetWallDist(): poisson [p={}] [{}] inc: [{:.4e}] -> [{:.4e}], res: [{:.4e}] -> [{:.4e}]",
742 pPoissonCur, iIter, incNormBase, incNorm, resNormBase, normR)
743 << std::endl;
744 if (nowExit)
745 break;
746 }
747 getDiffPhi(phi, diffPhi);
748 diffPhi.trans.startPersistentPull();
749 diffPhi.trans.waitPersistentPull();
750 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
751 {
752 Geom::tPoint gradPhi;
753 gradPhi.setZero();
754 gradPhi = diffPhi[iCell];
755 real nD = 1;
756 for (int ic2f = 0; ic2f < mesh->cell2face[iCell].size(); ic2f++)
757 {
758 index iFace = mesh->cell2face[iCell][ic2f];
759 index iCellOther = mesh->CellFaceOther(iCell, iFace);
760 if (iCellOther != UnInitIndex)
761 {
762 gradPhi += diffPhi[iCell];
763 gradPhi += diffPhi[iCellOther] * smoothBias;
764 nD += 1. + smoothBias;
765 }
766 }
767 gradPhi /= nD;
768 real gradPhiNorm = gradPhi.norm();
769 real dEst = std::pow(pPoissonCur / real(pPoissonCur - 1) * phi[iCell](0) + std::pow(gradPhiNorm, pPoissonCur), real(pPoissonCur - 1) / pPoissonCur) - std::pow(gradPhiNorm, pPoissonCur - 1);
770 // dPhi[iCell](0) = phi[iCell](0);
771 dPhi[iCell](0) = dEst;
772 }
773 dPhi.trans.startPersistentPull();
774 dPhi.trans.waitPersistentPull();
775 phi = dPhi;
776 }
777
778 auto minval = dPhi.min();
779
780 dWall.resize(mesh->NumCellProc());
781 for (index iCell = 0; iCell < mesh->NumCellProc(); iCell++)
782 {
783 auto quadCell = vfv->GetCellQuad(iCell);
784 dWall[iCell].resize(quadCell.GetNumPoints());
785 }
786
787 for (index iCell = 0; iCell < mesh->NumCellProc(); iCell++)
788 dWall.at(iCell).setConstant(std::max(dPhi[iCell](0), settings.minWallDist));
789 if (phi.father->getMPI().rank == 0)
790 log() << fmt::format("EulerEvaluator<model>::GetWallDist(): poisson min dist [{}]", minval(0)) << std::endl;
791 }
792
793 DNDS_SWITCH_INTELLISENSE(
794 template <EulerModel model>
795 ,
796 template <>)
797 /** @brief Dispatcher for wall distance computation, selecting method based on settings.
798 *
799 * Routes to GetWallDist_AABB (schemes 0,1,20), GetWallDist_BatchedAABB (schemes 2,20),
800 * or GetWallDist_Poisson (scheme 3), then computes face-averaged distances.
801 */
802 void EulerEvaluator<model>::GetWallDist()
803 {
804 if (settings.wallDistScheme == 0 || settings.wallDistScheme == 1 || settings.wallDistScheme == 20)
805 GetWallDist_AABB();
806 if (settings.wallDistScheme == 2 || settings.wallDistScheme == 20)
807 GetWallDist_BatchedAABB();
808 if (settings.wallDistScheme == 3)
809 GetWallDist_Poisson();
810 GetWallDist_ComputeFaceDistances();
811 }
812
813 // Eigen::Vector<real, -1> EulerEvaluator::CompressRecPart(
814 // const Eigen::Vector<real, -1> &umean,
815 // const Eigen::Vector<real, -1> &uRecInc)1
816
817 //! evaluates dt and facial spectral radius
818 DNDS_SWITCH_INTELLISENSE(
819 template <EulerModel model>
820 ,
821 template <>)
822 /** @brief Estimate the local time step and face spectral radii for CFL-based time stepping.
823 *
824 * For each face, evaluates L/R states using 2nd-order reconstruction, computes the
825 * convective and viscous eigenvalue estimates, and stores face spectral radii
826 * (lambdaFace, lambdaFace0..4, lambdaFaceC, lambdaFaceVis). Then for each cell,
827 * sums the face spectral radii to obtain a CFL-limited local time step, optionally
828 * clamped by MaxDt or made uniform (global minimum).
829 *
830 * @param dt Local time step per cell (output).
831 * @param u Conservative variable DOF array.
832 * @param uRec Reconstruction coefficients.
833 * @param CFL CFL number for time step scaling.
834 * @param dtMinall Global minimum time step (output, MPI-reduced).
835 * @param MaxDt Maximum allowed time step (upper clamp).
836 * @param UseLocaldt If true, use local time stepping; otherwise set all cells to dtMinall.
837 * @param t Current simulation time.
838 * @param flags Bitfield flags (e.g., DT_Dont_update_lambda01234).
839 */
840 void EulerEvaluator<model>::EvaluateDt(
841 ArrayDOFV<1> &dt,
842 ArrayDOFV<nVarsFixed> &u,
843 ArrayRECV<nVarsFixed> &uRec,
844 real CFL, real &dtMinall, real MaxDt,
845 bool UseLocaldt,
846 real t,
847 uint64_t flags)
848 {
849 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
850 DNDS_MPI_InsertCheck(u.father->getMPI(), "EvaluateDt 1");
851
852 bool dont_update_lambda01234 = flags & DT_Dont_update_lambda01234;
853
854 { // 2nd order reconstruction
855 typename TVFV::template TFBoundary<nVarsFixed>
856 FBoundary = [&](const TU &UL, const TU &UMean, index iCell, index iFace, int ig,
857 const Geom::tPoint &normOut, const Geom::tPoint &pPhy, const Geom::t_index bType) -> TU
858 {
859 TVec normOutV = normOut(Seq012);
860 Eigen::Matrix<real, dim, dim> normBase = Geom::NormBuildLocalBaseV<dim>(normOutV);
861 bool compressed = false;
862 TU ULfixed = this->CompressRecPart(
863 UMean,
864 UL - UMean,
865 compressed);
866 return this->generateBoundaryValue(ULfixed, UMean, iCell, iFace, ig, normOutV, normBase, pPhy, t, bType, true, 1);
867 };
868 vfv->DoReconstruction2ndGrad(uGradBufNoLim, u, FBoundary, settings.direct2ndRecMethod);
869 uGradBufNoLim.trans.startPersistentPull();
870 uGradBufNoLim.trans.waitPersistentPull();
871 }
872#if defined(DNDS_DIST_MT_USE_OMP)
873# pragma omp parallel for schedule(runtime)
874#endif
875 for (index iFace = 0; iFace < mesh->NumFaceProc(); iFace++)
876 {
877 auto f2c = mesh->face2cell[iFace];
878 TVec unitNorm = vfv->GetFaceNorm(iFace, -1)(Seq012);
879
880 index iCellL = f2c[0];
881 TU UL = u[iCellL];
882 this->UFromCell2Face(UL, iFace, f2c[0], 0);
883 TU uMean = UL;
884 real pL, asqrL, HL, pR, asqrR, HR;
885 TVec vL = UL(Seq123) / UL(0);
886 TVec vR = vL;
887 Gas::IdealGasThermal(UL(I4), UL(0), vL.squaredNorm(),
888 settings.idealGasProperty.gamma,
889 pL, asqrL, HL);
890 pR = pL, HR = HL, asqrR = asqrL;
891 if (f2c[1] != UnInitIndex)
892 {
893 TU UR = u[f2c[1]];
894 this->UFromCell2Face(UR, iFace, f2c[1], 1);
895 uMean = (uMean + UR) * 0.5;
896 vR = UR(Seq123) / UR(0);
897 Gas::IdealGasThermal(UR(I4), UR(0), vR.squaredNorm(),
898 settings.idealGasProperty.gamma,
899 pR, asqrR, HR);
900 }
901 TDiffU GradULxy, GradURxy;
902 GradULxy.resize(Eigen::NoChange, nVars);
903 GradURxy.resize(Eigen::NoChange, nVars);
904 GradULxy.setZero(), GradURxy.setZero();
905 GradULxy(SeqG012, EigenAll) = uGradBuf[f2c[0]];
906 // if constexpr (gDim == 2)
907 // GradULxy({0, 1}, EigenAll) =
908 // vfv->GetIntPointDiffBaseValue(f2c[0], iFace, 0, -1, std::array<int, 2>{1, 2}, 3) *
909 // uRec[f2c[0]]; // 2d here
910 // else
911 // GradULxy({0, 1, 2}, EigenAll) =
912 // vfv->GetIntPointDiffBaseValue(f2c[0], iFace, 0, -1, std::array<int, 3>{1, 2, 3}, 4) *
913 // uRec[f2c[0]]; // 3d here
914 this->DiffUFromCell2Face(GradULxy, iFace, f2c[0], 0);
915 GradURxy = GradULxy;
916 if (f2c[1] != UnInitIndex)
917 {
918 GradURxy(SeqG012, EigenAll) = uGradBuf[f2c[1]];
919 // if constexpr (gDim == 2)
920 // GradURxy({0, 1}, EigenAll) =
921 // vfv->GetIntPointDiffBaseValue(f2c[1], iFace, 1, -1, std::array<int, 2>{1, 2}, 3) *
922 // uRec[f2c[1]]; // 2d here
923 // else
924 // GradURxy({0, 1, 2}, EigenAll) =
925 // vfv->GetIntPointDiffBaseValue(f2c[1], iFace, 1, -1, std::array<int, 3>{1, 2, 3}, 4) *
926 // uRec[f2c[1]]; // 3d here
927 this->DiffUFromCell2Face(GradURxy, iFace, f2c[1], 1);
928 }
929 TDiffU GradUMeanXY = (GradURxy + GradULxy) / 2;
930
931 DNDS_assert(uMean(0) > 0);
932 TVec veloMean = (uMean(Seq123).array() / uMean(0)).matrix();
933 // real veloNMean = veloMean.dot(unitNorm); // original
934 real veloNMean = 0.5 * (vL + vR).dot(unitNorm); // paper
935 real veloNL = vL.dot(unitNorm);
936 real veloNR = vR.dot(unitNorm);
937 real vgN = this->GetFaceVGrid(iFace, -1).dot(unitNorm);
938
939 // real ekFixRatio = 0.001;
940 // Eigen::Vector3d velo = uMean({1, 2, 3}) / uMean(0);
941 // real vsqr = velo.squaredNorm();
942 // real Ek = vsqr * 0.5 * uMean(0);
943 // real Efix = Ek * ekFixRatio;
944 // real e = uMean(4) - Ek;
945 // if (e < 0)
946 // e = 0.5 * Efix;
947 // else if (e < Efix)
948 // e = (e * e + Efix * Efix) / (2 * Efix);
949 // uMean(4) = Ek + e;
950
951 real pMean, asqrMean, HMean;
952 Gas::IdealGasThermal(uMean(I4), uMean(0), veloMean.squaredNorm(),
953 settings.idealGasProperty.gamma,
954 pMean, asqrMean, HMean);
955
956 pMean = (pL + pR) * 0.5;
957 real aMean = sqrt(settings.idealGasProperty.gamma * pMean / uMean(0)); // paper
958
959 // DNDS_assert(asqrMean >= 0);
960 // real aMean = std::sqrt(asqrMean); // original
961 real lambdaConvection = std::abs(veloNMean - vgN) + aMean;
962 DNDS_assert_info(
963 asqrL >= 0 && asqrR >= 0,
964 fmt::format(" mean value violates PP! asqr: [{} {}]", asqrL, asqrR));
965 real aL = std::sqrt(asqrL);
966 real aR = std::sqrt(asqrR);
967 lambdaConvection = std::max(aL + std::abs(veloNL - vgN), aR + std::abs(veloNR - vgN));
968
969 // ! refvalue:
970 real muRef = settings.idealGasProperty.muGas;
971
972 real gamma = settings.idealGasProperty.gamma;
973 real T = pMean / ((gamma - 1) / gamma * settings.idealGasProperty.CpGas * uMean(0));
974 real muf = muEff(uMean, T);
975 real muPhy = muf;
976 real muTurb = this->getMuTur(uMean, GradUMeanXY, muRef, muf, iFace);
977 muf += muTurb;
978
979 real lamVis = muf / uMean(0) *
980 std::max(4. / 3., gamma / settings.idealGasProperty.prGas);
981
982 real lamFace = lambdaConvection * vfv->GetFaceArea(iFace);
983
984 real area = vfv->GetFaceArea(iFace);
985 real areaSqr = area * area;
986 real volR = vfv->GetCellVol(iCellL);
987 // lambdaCell[iCellL] += lamFace + 2 * lamVis * areaSqr / fv->GetCellVol(iCellL);
988 if (f2c[1] != UnInitIndex) // can't be non local
989 // lambdaCell[f2c[1]] += lamFace + 2 * lamVis * areaSqr / fv->volumeLocal[f2c[1]],
990 volR = vfv->GetCellVol(f2c[1]);
991
992 lambdaFace[iFace] = lambdaConvection + lamVis * area * (1. / vfv->GetCellVol(iCellL) + 1. / volR);
993 lambdaFaceC[iFace] = std::abs(veloNMean - vgN) + lamVis * area * (1. / vfv->GetCellVol(iCellL) + 1. / volR); // passive part
994 lambdaFaceVis[iFace] = lamVis * area * (1. / vfv->GetCellVol(iCellL) + 1. / volR);
995 if (!dont_update_lambda01234)
996 lambdaFace0[iFace] = lambdaFace123[iFace] = lambdaFace4[iFace] = lambdaConvection;
997
998 // if (f2c[0] == 10756)
999 // {
1000 // std::cout << "----Lambdas" << std::setprecision(16) << iFace << std::endl;
1001 // std::cout << lambdaConvection << std::endl;
1002 // std::cout << lambdaFaceVis[iFace] << std::endl;
1003 // std::cout << veloNMean << " " << aMean << std::endl;
1004 // std::cout << gamma << " " << pMean << " " << uMean(0) << std::endl;
1005 // }
1006
1007 // put to cell part to be OMP compliant
1008 // lambdaCell[iCellL] += lambdaFace[iFace] * vfv->GetFaceArea(iFace);
1009 // if (f2c[1] != UnInitIndex) // can't be non local
1010 // lambdaCell[f2c[1]] += lambdaFace[iFace] * vfv->GetFaceArea(iFace);
1011
1012 deltaLambdaFace[iFace] = std::max<real>({std::abs((vR - vL).dot(unitNorm) + aR - aL), std::abs((vR - vL).dot(unitNorm) - aR + aL), std::abs(aL - aR)});
1013 }
1014 real dtMin = veryLargeReal;
1015#if defined(DNDS_DIST_MT_USE_OMP)
1016# pragma omp parallel for schedule(runtime) reduction(min : dtMin)
1017#endif
1018 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
1019 {
1020 real lambdaCellC = 0;
1021 for (auto iFace : mesh->cell2face[iCell])
1022 lambdaCellC += lambdaFace[iFace] * vfv->GetFaceArea(iFace);
1023 // std::cout << fv->GetCellVol(iCell) << " " << (lambdaCellC) << " " << CFL << std::endl;
1024 // exit(0);
1025 dt[iCell](0) = std::min(CFL * vfv->GetCellVol(iCell) * vfv->GetCellSmoothScaleRatio(iCell) / (lambdaCellC + 1e-100), MaxDt);
1026 dtMin = std::min(dtMin, dt[iCell](0));
1027 deltaLambdaCell[iCell](0) = 0;
1028 for (auto iFace : mesh->cell2face[iCell])
1029 deltaLambdaCell[iCell](0) = std::max(deltaLambdaCell[iCell](0), deltaLambdaFace[iFace]);
1030 // if (iCell == 10756)
1031 // {
1032 // std::cout << std::endl;
1033 // }
1034 }
1035
1036 deltaLambdaCell.trans.startPersistentPull();
1037 MPI::Allreduce(&dtMin, &dtMinall, 1, DNDS_MPI_REAL, MPI_MIN, u.father->getMPI().comm);
1038 deltaLambdaCell.trans.waitPersistentPull();
1039
1040 // if (uRec.father->getMPI().rank == 0)
1041 // std::cout << "dt min is " << dtMinall << std::endl;
1042 if (!UseLocaldt)
1043 {
1044 dt.setConstant(dtMinall);
1045 }
1046 // if (uRec.father->getMPI().rank == 0)
1047 // log() << "dt: " << dtMin << std::endl;
1048 }
1049
1050 DNDS_SWITCH_INTELLISENSE(
1051 // the real definition
1052 template <EulerModel model>
1053 ,
1054 // the intellisense friendly definition
1055 template <>)
1056 /** @brief Evaluate the numerical flux at a single face for a batch of quadrature points.
1057 *
1058 * Computes the inviscid and viscous numerical flux at a face given batched L/R states,
1059 * gradients, and face geometry. Dispatches to the configured Riemann solver type.
1060 * Outputs left/right viscous corrections (FLfix, FRfix), the total flux increment (finc),
1061 * and per-quadrature-point eigenvalue estimates (lam0V, lam123V, lam4V).
1062 *
1063 * @param ULxy Batched left states at quadrature points.
1064 * @param URxy Batched right states at quadrature points.
1065 * @param ULMeanXy Left cell mean state.
1066 * @param URMeanXy Right cell mean state.
1067 * @param DiffUxy Batched conservative gradient at face.
1068 * @param DiffUxyPrim Batched primitive gradient at face.
1069 * @param unitNorm Batched outward unit normals at quadrature points.
1070 * @param vgXY Batched grid velocity at quadrature points.
1071 * @param unitNormC Face-center unit normal.
1072 * @param vgC Face-center grid velocity.
1073 * @param FLfix Left viscous correction (output).
1074 * @param FRfix Right viscous correction (output).
1075 * @param finc Total flux increment (output).
1076 * @param lam0V Eigenvalue estimate for wave 0 (output).
1077 * @param lam123V Eigenvalue estimate for waves 1-3 (output).
1078 * @param lam4V Eigenvalue estimate for wave 4 (output).
1079 * @param btype Boundary type index of the face.
1080 * @param rsType Riemann solver type to use.
1081 * @param iFace Face index.
1082 * @param ignoreVis If true, skip viscous flux computation.
1083 */
1084 void EulerEvaluator<model>::fluxFace(
1085 const TU_Batch &ULxy,
1086 const TU_Batch &URxy,
1087 const TU &ULMeanXy,
1088 const TU &URMeanXy,
1089 const TDiffU_Batch &DiffUxy,
1090 const TDiffU_Batch &DiffUxyPrim,
1091 const TVec_Batch &unitNorm,
1092 const TVec_Batch &vgXY,
1093 const TVec &unitNormC,
1094 const TVec &vgC,
1095 TU_Batch &FLfix,
1096 TU_Batch &FRfix,
1097 TU_Batch &finc,
1098 TReal_Batch &lam0V, TReal_Batch &lam123V, TReal_Batch &lam4V,
1099 Geom::t_index btype,
1100 typename Gas::RiemannSolverType rsType,
1101 index iFace, bool ignoreVis)
1102 {
1103 finc.resizeLike(ULxy);
1104 DNDS_assert(&DiffUxy != &DiffUxyPrim);
1105 DNDS_assert(&ULMeanXy != &URMeanXy);
1106 DNDS_assert(&DiffUxy != &DiffUxyPrim);
1107 DNDS_assert(&unitNorm != &vgXY);
1108 DNDS_assert(&unitNormC != &vgC);
1109 DNDS_assert(&FLfix != &FRfix);
1110 DNDS_assert(&lam0V != &lam123V);
1111 DNDS_assert(&lam123V != &lam4V); // aliasing check not complete here
1112
1113 auto fluxFace_impl =
1114 [&](
1115 const TU_Batch *DNDS_RESTRICT p_ULxy,
1116 const TU_Batch *DNDS_RESTRICT p_URxy,
1117 const TU *DNDS_RESTRICT p_ULMeanXy,
1118 const TU *DNDS_RESTRICT p_URMeanXy,
1119 const TDiffU_Batch *DNDS_RESTRICT p_DiffUxy,
1120 const TDiffU_Batch *DNDS_RESTRICT p_DiffUxyPrim,
1121 const TVec_Batch *DNDS_RESTRICT p_unitNorm,
1122 const TVec_Batch *DNDS_RESTRICT p_vgXY,
1123 const TVec *DNDS_RESTRICT p_unitNormC,
1124 const TVec *DNDS_RESTRICT p_vgC,
1125 TU_Batch *DNDS_RESTRICT p_FLfix,
1126 TU_Batch *DNDS_RESTRICT p_FRfix,
1127 TReal_Batch *DNDS_RESTRICT p_lam0V, TReal_Batch *DNDS_RESTRICT p_lam123V, TReal_Batch *DNDS_RESTRICT p_lam4V,
1128 Geom::t_index btype,
1129 typename Gas::RiemannSolverType rsType)
1130 // clang-format off
1131 {
1132 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
1133 const TU_Batch &ULxy = *p_ULxy;
1134 const TU_Batch &URxy = *p_URxy;
1135 const TU &ULMeanXy = *p_ULMeanXy;
1136 const TU &URMeanXy = *p_URMeanXy;
1137 const TDiffU_Batch &DiffUxy = *p_DiffUxy;
1138 const TDiffU_Batch &DiffUxyPrim = *p_DiffUxyPrim;
1139 const TVec_Batch &unitNorm = *p_unitNorm;
1140 const TVec_Batch &vgXY = *p_vgXY;
1141 const TVec &unitNormC = *p_unitNormC;
1142 const TVec &vgC = *p_vgC;
1143 TU_Batch &FLfix = *p_FLfix;
1144 TU_Batch &FRfix = *p_FRfix;
1145 TReal_Batch &lam0V = *p_lam0V;
1146 TReal_Batch &lam123V = *p_lam123V;
1147 TReal_Batch &lam4V = *p_lam4V;
1148
1149 int nB = ULxy.cols();
1150
1151 TU_Batch UMeanXy = 0.5 * (ULxy + URxy);
1152
1153 // PerformanceTimer::Instance().StartTimer(PerformanceTimer::LimiterB);
1154 real muRef = settings.idealGasProperty.muGas;
1155 real TRef = settings.idealGasProperty.TRef;
1156
1157 auto f2c = mesh->face2cell[iFace];
1158 real dLambda = deltaLambdaCell[f2c[0]](0);
1159 if (f2c[1] != UnInitIndex)
1160 dLambda = std::max(dLambda, deltaLambdaCell[f2c[1]](0));
1161 real fixScale = settings.rsFixScale;
1162 real incFScale = settings.rsIncFScale;
1163 /** viscous flux **/
1164 TU_Batch visFluxV;
1165 visFluxV.resizeLike(ULxy);
1166 for (int iB = 0; iB < nB; iB++)
1167 {
1168 real pMean, asqrMean, Hmean;
1169 real gamma = settings.idealGasProperty.gamma;
1170 TU UMeanXYC = UMeanXy(EigenAll, iB);
1171 auto seqC = Eigen::seq(iB * dim, iB * dim + dim - 1);
1172 TDiffU DiffUxyC = DiffUxy(seqC, EigenAll);
1173 TDiffU DiffUxyPrimC = DiffUxyPrim(seqC, EigenAll);
1174 TVec uNormC = unitNorm(EigenAll, iB);
1175 Gas::IdealGasThermal(UMeanXYC(I4), UMeanXYC(0), (UMeanXYC(Seq123) / UMeanXYC(0)).squaredNorm(),
1176 gamma, pMean, asqrMean, Hmean);
1177 DNDS_assert_info(pMean > 0, fmt::format("{}, {}, {}", UMeanXYC(I4), UMeanXYC(0), (UMeanXYC(Seq123) / UMeanXYC(0)).squaredNorm()));
1178 real T = pMean / ((gamma - 1) / gamma * settings.idealGasProperty.CpGas * UMeanXYC(0));
1179 real mufPhy, muf;
1180 muf = muEff(UMeanXYC, T);
1181 mufPhy = muf;
1182 // PerformanceTimer::Instance().StopTimer(PerformanceTimer::LimiterB);
1183#ifndef DNDS_FV_EULEREVALUATOR_IGNORE_VISCOUS_TERM
1184 if (!ignoreVis)
1185 {
1186 real muTurb = this->getMuTur(UMeanXYC, DiffUxyC, muRef, muf, iFace); // TODO: make this accept primitive gradients instead
1187 muf += muTurb;
1188
1189 real k = settings.idealGasProperty.CpGas * muTurb / 0.9 +
1190 settings.idealGasProperty.CpGas * mufPhy / settings.idealGasProperty.prGas;
1191 TU VisFlux;
1192 VisFlux.resizeLike(ULMeanXy);
1193 VisFlux.setZero();
1194 Gas::ViscousFlux_IdealGas<dim>(
1195 UMeanXYC, DiffUxyPrimC, uNormC, pBCHandler->GetTypeFromID(btype) == EulerBCType::BCWall,
1196 settings.idealGasProperty.gamma,
1197 muf, muTurb / (muf + verySmallReal), settings.ransUseQCR,
1198 k,
1199 settings.idealGasProperty.CpGas,
1200 VisFlux);
1201
1202 this->visFluxTurVariable(UMeanXYC, DiffUxyPrimC, muRef, mufPhy, muTurb, uNormC, iFace, VisFlux);
1203 if (pBCHandler->GetTypeFromID(btype) == EulerBCType::BCWallInvis ||
1204 pBCHandler->GetTypeFromID(btype) == EulerBCType::BCSym)
1205 {
1206 VisFlux *= 0.0;
1207 }
1208 visFluxV(EigenAll, iB) = VisFlux;
1209 }
1210#endif
1211 if (!isfinite(pMean) || !isfinite(pMean) || !isfinite(pMean))
1212 {
1213 std::cout << T << std::endl;
1214 std::cout << muf << std::endl;
1215 std::cout << pMean << std::endl;
1216 DNDS_assert(false);
1217 }
1218 }
1219
1220 auto exitFun = [&]()
1221 {
1222 std::cout << "face at" << vfv->GetFaceQuadraturePPhys(iFace, -1) << '\n';
1223 std::cout << "UL" << ULxy.transpose() << '\n';
1224 std::cout << "UR" << URxy.transpose() << std::endl;
1225 std::cout << "ULM" << ULMeanXy.transpose() << '\n';
1226 std::cout << "URM" << URMeanXy.transpose() << std::endl;
1227 };
1228
1229 if (pBCHandler->GetTypeFromID(btype) == EulerBCType::BCWallInvis ||
1230 pBCHandler->GetTypeFromID(btype) == EulerBCType::BCSym)
1231 {
1232 // ? normal invert here?
1233 }
1234
1235 lam0V.resize(nB);
1236 lam123V.resize(nB);
1237 lam4V.resize(nB);
1238
1239 auto RSWrapper_XY =
1240 [&](Gas::RiemannSolverType rsType,
1241 auto &&UL, auto &&UR, auto &&ULm, auto &&URm, auto &&vg, auto &&n,
1242 real gamma, auto &&finc, real dLambda, real fixScale, real incFScale,
1243 real &lam0, real &lam123, real &lam4)
1244 {
1245 Gas::InviscidFlux_IdealGas_Dispatcher<dim>(rsType, UL, UR, ULm, URm, vg, n, gamma, finc, dLambda, fixScale, incFScale,exitFun, lam0, lam123, lam4);
1246 };
1247
1248 // TU_Batch finc1;
1249 // finc1.resizeLike(ULxy);
1250 if (settings.rsTypeWall != Gas::UnknownRS &&
1251 (pBCHandler->GetTypeFromID(btype) == EulerBCType::BCWall ||
1252 pBCHandler->GetTypeFromID(btype) == BCWallIsothermal))
1253 {
1254 rsType = settings.rsTypeWall;
1255 }
1256
1257 auto runRsOnNorm = [&]()
1258 {
1259 if (settings.rsMeanValueEig != 0 &&
1260 (rsType >= Gas::Roe_M1 && rsType <= Gas::Roe_M5))
1261 {
1262 real lam0{0}, lam123{0}, lam4{0};
1263 Gas::InviscidFlux_IdealGas_Batch_Dispatcher<dim>(
1264 rsType,
1265 ULxy, URxy, ULMeanXy, URMeanXy, vgXY, vgC, unitNorm, unitNormC,
1266 settings.idealGasProperty.gamma, finc, dLambda, fixScale, incFScale,
1267 exitFun, lam0, lam123, lam4);
1268 lam0V.setConstant(lam0);
1269 lam123V.setConstant(lam123);
1270 lam4V.setConstant(lam4);
1271 }
1272 else
1273 for (int iB = 0; iB < nB; iB++)
1274 {
1275 RSWrapper_XY(rsType, ULxy(EigenAll, iB), URxy(EigenAll, iB), ULMeanXy, URMeanXy,
1276 vgXY(EigenAll, iB), unitNorm(EigenAll, iB),
1277 settings.idealGasProperty.gamma, finc(EigenAll, iB), dLambda, fixScale, incFScale,
1278 lam0V(iB), lam123V(iB), lam4V(iB));
1279 }
1280 };
1281
1282 if (settings.rsRotateScheme == 0)
1283 {
1284 runRsOnNorm();
1285 }
1286 else
1287 {
1288 TVec veloL = ULMeanXy(Seq123) / ULMeanXy(0);
1289 TVec veloR = URMeanXy(Seq123) / URMeanXy(0);
1290 TVec diffVelo = veloR - veloL;
1291 real diffVeloN = diffVelo.norm();
1292 real veloLN = veloL.norm();
1293 real veloRN = veloR.norm();
1294 if (diffVeloN < (smallReal * 10) * (veloLN + veloRN) || diffVeloN < std::sqrt(verySmallReal))
1295 runRsOnNorm();
1296 else
1297 {
1298 TVec N1;
1299 if (settings.rsRotateScheme == 1)
1300 N1 = diffVelo / diffVeloN;
1301 else if (settings.rsRotateScheme == 2)
1302 N1 = unitNormC;
1303 else
1304 DNDS_assert(false);
1305 DNDS_assert_info(std::abs(N1.norm() - 1) < 1e-5, fmt::format("{}", diffVeloN));
1306 N1 *= sign(N1.dot(unitNormC));
1307 TReal_Batch N1Proj = N1.transpose() * unitNorm;
1308
1309 TVec_Batch N2 = unitNorm - N1 * N1Proj;
1310 TReal_Batch N2Proj = N2.colwise().norm().array().max(verySmallReal * 10);
1311 N2.array().rowwise() /= N2Proj.array();
1312
1313 real N1ProjC = N1.dot(unitNormC);
1314 TVec N2C = unitNormC - N1 * N1ProjC;
1315 real N2CProj = std::max(N2C.norm(), verySmallReal * 10);
1316 N2C /= N2CProj;
1317
1318 TVec_Batch N1B;
1319 N1B.resizeLike(N2);
1320 N1B.colwise() = N1;
1321
1322 if (false)
1323 {
1324 std::cout << unitNorm << "\n";
1325 std::cout << "N1" << "\n";
1326 std::cout << N1.transpose() << "\n";
1327 std::cout << N1B.transpose() << "\n";
1328 std::cout << N1Proj << "\n";
1329 std::cout << "N2" << "\n";
1330 std::cout << N2.transpose() << "\n";
1331 std::cout << N2Proj << "\n";
1332 std::cout << N2C.transpose() << "\n";
1333 std::cout << N2CProj << "\n";
1334 std::cout << std::endl;
1335 }
1336
1337 TReal_Batch lam4V1, lam0V1, lam123V1;
1338 lam0V1.resizeLike(lam0V);
1339 lam4V1.resizeLike(lam0V);
1340 lam123V1.resizeLike(lam0V);
1341
1342 TU_Batch F1;
1343 F1.resizeLike(finc);
1344
1345 auto rsTypeAux = settings.rsTypeAux == Gas::UnknownRS ? rsType : settings.rsTypeAux;
1346
1347 if (settings.rsMeanValueEig != 0 &&
1348 (rsTypeAux >= Gas::Roe_M1 && rsTypeAux <= Gas::Roe_M5))
1349 {
1350 real lam0{0}, lam123{0}, lam4{0};
1351 Gas::InviscidFlux_IdealGas_Batch_Dispatcher<dim>(
1352 rsTypeAux,
1353 ULxy, URxy, ULMeanXy, URMeanXy, vgXY, vgC, N1B, N1,
1354 settings.idealGasProperty.gamma, F1, dLambda, fixScale, incFScale,
1355 exitFun, lam0, lam123, lam4);
1356 lam0V1.setConstant(lam0);
1357 lam123V1.setConstant(lam123);
1358 lam4V1.setConstant(lam4);
1359 }
1360 else
1361 for (int iB = 0; iB < nB; iB++)
1362 {
1363 RSWrapper_XY(rsTypeAux, ULxy(EigenAll, iB), URxy(EigenAll, iB), ULMeanXy, URMeanXy,
1364 vgXY(EigenAll, iB), N1,
1365 settings.idealGasProperty.gamma, F1(EigenAll, iB), dLambda, fixScale, incFScale,
1366 lam0V1(iB), lam123V1(iB), lam4V1(iB));
1367 }
1368
1369 if (settings.rsMeanValueEig != 0 &&
1370 (rsType >= Gas::Roe_M1 && rsType <= Gas::Roe_M5))
1371 {
1372 real lam0{0}, lam123{0}, lam4{0};
1373 Gas::InviscidFlux_IdealGas_Batch_Dispatcher<dim>(
1374 rsType,
1375 ULxy, URxy, ULMeanXy, URMeanXy, vgXY, vgC, N2, N2C,
1376 settings.idealGasProperty.gamma, finc, dLambda, fixScale, incFScale,
1377 exitFun, lam0, lam123, lam4);
1378 lam0V.setConstant(lam0);
1379 lam123V.setConstant(lam123);
1380 lam4V.setConstant(lam4);
1381 }
1382 else
1383 for (int iB = 0; iB < nB; iB++)
1384 {
1385 RSWrapper_XY(rsType, ULxy(EigenAll, iB), URxy(EigenAll, iB), ULMeanXy, URMeanXy,
1386 vgXY(EigenAll, iB), N2(EigenAll, iB),
1387 settings.idealGasProperty.gamma, finc(EigenAll, iB), dLambda, fixScale, incFScale,
1388 lam0V(iB), lam123V(iB), lam4V(iB));
1389 }
1390
1391 finc.array().rowwise() *= N2Proj.array();
1392 finc.array() += F1.array().rowwise() * N1Proj.array();
1393
1394 TReal_Batch N12ProjSum = N1Proj.array() + N2Proj.array();
1395 lam0V.array() *= N2Proj.array() / N12ProjSum.array();
1396 lam4V.array() *= N2Proj.array() / N12ProjSum.array();
1397 lam123V.array() *= N2Proj.array() / N12ProjSum.array();
1398 lam0V.array() += N1Proj.array() * lam0V1.array() / N12ProjSum.array();
1399 lam4V.array() += N1Proj.array() * lam4V1.array() / N12ProjSum.array();
1400 lam123V.array() += N1Proj.array() * lam123V1.array() / N12ProjSum.array(); // todo: fix these
1401
1402 lam0V = lam0V1;
1403 lam123V = lam123V1;
1404 lam4V = lam4V1;
1405 }
1406 }
1407
1408#ifndef USE_ENTROPY_FIXED_LAMBDA_IN_SA
1409 lam123 = (std::abs(UL(1) / UL(0) - vg(0)) + std::abs(UR(1) / UR(0) - vg(0))) * 0.5; //! high fix
1410 // lam123 = std::abs(UL(1) / UL(0) + UR(1) / UR(0)) * 0.5 - vg(0); //! low fix
1411#endif
1412
1413 if constexpr (Traits::hasSA)
1414 {
1415 // real lambdaFaceCC = sqrt(std::abs(asqrMean)) + std::abs((UL(1) / UL(0) - vg(0)) + (UR(1) / UR(0) - vg(0))) * 0.5;
1416 Eigen::RowVector<real, -1> lambdaFaceCC = lam123V; //! using velo instead of velo + a
1417 if (settings.ransEigScheme == 1)
1418 lambdaFaceCC = lambdaFaceCC.array().max(lam0V.array()).max(lam4V.array());
1419 auto vnR = ((URxy(Seq123, EigenAll).array().rowwise() / URxy(0, EigenAll).array() - vgXY.array()) * unitNorm.array()).colwise().sum();
1420 auto vnL = ((ULxy(Seq123, EigenAll).array().rowwise() / ULxy(0, EigenAll).array() - vgXY.array()) * unitNorm.array()).colwise().sum();
1421 finc(I4 + 1, EigenAll) =
1422 ((vnL * ULxy(I4 + 1, EigenAll).array() + vnR * URxy(I4 + 1, EigenAll).array()) -
1423 (URxy(I4 + 1, EigenAll).array() - ULxy(I4 + 1, EigenAll).array()) * lambdaFaceCC.array()) *
1424 0.5;
1425 }
1426 if constexpr (Traits::has2EQ)
1427 {
1428 Eigen::RowVector<real, -1> lambdaFaceCC = lam123V; //! using velo instead of velo + a
1429 if (settings.ransEigScheme == 1)
1430 lambdaFaceCC = lambdaFaceCC.array().max(lam0V.array()).max(lam4V.array());
1431 auto vnR = ((URxy(Seq123, EigenAll).array().rowwise() / URxy(0, EigenAll).array() - vgXY.array()) * unitNorm.array()).colwise().sum();
1432 auto vnL = ((ULxy(Seq123, EigenAll).array().rowwise() / ULxy(0, EigenAll).array() - vgXY.array()) * unitNorm.array()).colwise().sum();
1433 finc(I4 + 1, EigenAll) =
1434 ((vnL * ULxy(I4 + 1, EigenAll).array() + vnR * URxy(I4 + 1, EigenAll).array()) -
1435 (URxy(I4 + 1, EigenAll).array() - ULxy(I4 + 1, EigenAll).array()) * lambdaFaceCC.array()) *
1436 0.5;
1437 finc(I4 + 2, EigenAll) =
1438 ((vnL * ULxy(I4 + 2, EigenAll).array() + vnR * URxy(I4 + 2, EigenAll).array()) -
1439 (URxy(I4 + 2, EigenAll).array() - ULxy(I4 + 2, EigenAll).array()) * lambdaFaceCC.array()) *
1440 0.5;
1441 finc(1, EigenAll).array() += UMeanXy(I4 + 1, EigenAll).array() * (2. / 3.); //! k's normal stress
1442 finc(I4, EigenAll).array() += UMeanXy(I4 + 1, EigenAll).array() * (2. / 3.) * UMeanXy(1, EigenAll).array() / UMeanXy(0, EigenAll).array();
1443 }
1444 if constexpr (Traits::isExtended)
1445 {
1446 // real lambdaFaceCC = sqrt(std::abs(asqrMean)) + std::abs((UL(1) / UL(0) - vg(0)) + (UR(1) / UR(0) - vg(0))) * 0.5;
1447 Eigen::RowVector<real, -1> lambdaFaceCC = lam123V; //! using velo instead of velo + a
1448 if (settings.ransEigScheme == 1)
1449 lambdaFaceCC = lambdaFaceCC.array().max(lam0V.array()).max(lam4V.array());
1450 auto vnR = ((URxy(Seq123, EigenAll).array().rowwise() / URxy(0, EigenAll).array() - vgXY.array()) * unitNorm.array()).colwise().sum();
1451 auto vnL = ((ULxy(Seq123, EigenAll).array().rowwise() / ULxy(0, EigenAll).array() - vgXY.array()) * unitNorm.array()).colwise().sum();
1452 finc(SeqI52Last, EigenAll) =
1453 ((vnL * ULxy(SeqI52Last, EigenAll).array() + vnR * URxy(SeqI52Last, EigenAll).array()) -
1454 (URxy(SeqI52Last, EigenAll).array() - ULxy(SeqI52Last, EigenAll).array()) * lambdaFaceCC.array()) *
1455 0.5;
1456 }
1457
1458#ifndef DNDS_FV_EULEREVALUATOR_IGNORE_VISCOUS_TERM
1459 if (!ignoreVis)
1460 finc -= visFluxV;
1461#endif
1462
1463 if (finc.hasNaN() || (!finc.allFinite()))
1464 {
1465 std::cout << "finc\n"
1466 << finc.transpose() << "\n";
1467 std::cout << "visFluxV\n"
1468 << visFluxV << "\n";
1469 std::cout << "lam0V\n"
1470 << lam0V << "\n";
1471 std::cout << "lam123V\n"
1472 << lam0V << "\n";
1473 std::cout << "lam4V\n"
1474 << lam0V << "\n";
1475 std::cout << "ULxy\n"
1476 << ULxy.transpose() << "\n";
1477 std::cout << "URxy\n"
1478 << URxy.transpose() << "\n";
1479 std::cout << "UMeanXy\n"
1480 << UMeanXy.transpose() << "\n";
1481 std::cout << "DiffUxy\n"
1482 << DiffUxy << "\n";
1483 std::cout << "DiffUxyPrim\n"
1484 << DiffUxyPrim << "\n";
1485 std::cout << "unitNorm\n"
1486 << unitNorm << "\n";
1487 std::cout << "btype\n"
1488 << btype << std::endl;
1489 DNDS_assert(false);
1490 }
1491
1492 finc *= -1.0;
1493 };
1494 // clang-format on
1495 return fluxFace_impl(
1496 &ULxy,
1497 &URxy,
1498 &ULMeanXy,
1499 &URMeanXy,
1500 &DiffUxy,
1501 &DiffUxyPrim,
1502 &unitNorm,
1503 &vgXY,
1504 &unitNormC,
1505 &vgC,
1506 &FLfix,
1507 &FRfix,
1508 &lam0V, &lam123V, &lam4V,
1509 btype,
1510 rsType);
1511 }
1512
1513 DNDS_SWITCH_INTELLISENSE(
1514 template <EulerModel model>
1515 ,
1516 template <>)
1517 /** @brief Compute source terms for a cell at a quadrature point.
1518 *
1519 * Evaluates volume source terms including: constant mass force, rotating-frame
1520 * fictitious forces (Coriolis and centrifugal), axisymmetric source terms,
1521 * SA turbulence model source, and k-omega/k-epsilon RANS model sources.
1522 *
1523 * @param UMeanXy Cell mean conservative state.
1524 * @param DiffUxy Gradient of conservative variables at the quadrature point.
1525 * @param pPhy Physical coordinates of the quadrature point.
1526 * @param jacobian Source-term Jacobian matrix (output when Mode==2).
1527 * @param iCell Cell index.
1528 * @param ig Quadrature point index within the cell (-1 for cell center).
1529 * @param Mode 0: return source vector; 1: return diagonal Jacobian; 2: return full Jacobian.
1530 * @return Source term vector.
1531 */
1532 typename EulerEvaluator<model>::TU EulerEvaluator<model>::source(
1533 const TU &UMeanXy,
1534 const TDiffU &DiffUxy,
1535 const Geom::tPoint &pPhy,
1536 TJacobianU &jacobian,
1537 index iCell,
1538 index ig,
1539 int Mode) // mode =0: source; mode = 1, diagJacobi; mode = 2,
1540 {
1541 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
1542 TU ret;
1543 ret.resizeLike(UMeanXy);
1544 ret.setZero();
1545 real dWallC;
1546 if (ig < 0)
1547 dWallC = dWall[iCell].mean();
1548 else
1549 dWallC = dWall[iCell][ig];
1550 if (Mode == 2)
1551 jacobian.setZero(UMeanXy.size(), UMeanXy.size());
1552#ifdef DNDS_FV_EULEREVALUATOR_SOURCE_TERM_ZERO
1553 return ret;
1554#endif
1555 if (Mode == 0)
1556 {
1557 ret(Seq123) += settings.constMassForce(Seq012) * UMeanXy(0);
1558 ret(I4) += settings.constMassForce(Seq012).dot(UMeanXy(Seq123));
1559 }
1560 if (Mode == 2)
1561 {
1562 jacobian(I4, Seq123) -= settings.constMassForce(Seq012);
1563 }
1564#ifdef USE_ABS_VELO_IN_ROTATION
1565 if (settings.frameConstRotation.enabled)
1566 {
1567 if (Mode == 0 || Mode == 2)
1568 ret(Seq123) += -settings.frameConstRotation.vOmega().cross(Geom::ToThreeDim<dim>(UMeanXy(Seq123)))(Seq012);
1569 if (Mode == 2)
1570 jacobian(Seq123, Seq123) -= Geom::CrossVecToMat(-settings.frameConstRotation.vOmega())(Seq012, Seq012);
1571 }
1572#else
1573 if (settings.frameConstRotation.enabled)
1574 {
1575 Geom::tPoint radi = pPhy - settings.frameConstRotation.center;
1576 Geom::tPoint radiR = radi - settings.frameConstRotation.axis * (settings.frameConstRotation.axis.dot(radi));
1577 TVec mvolForce = (radiR * sqr(settings.frameConstRotation.Omega()) * UMeanXy(0))(Seq012);
1578 mvolForce += -2.0 * settings.frameConstRotation.vOmega().cross(Geom::ToThreeDim<dim>(UMeanXy(Seq123)))(Seq012);
1579 if (Mode == 0)
1580 {
1581 ret(Seq123) += mvolForce;
1582 ret(I4) += mvolForce.dot(UMeanXy(Seq123)) / UMeanXy(0);
1583 }
1584 if (Mode == 2)
1585 {
1586 TMat dmvolForceDrhov = Geom::CrossVecToMat(-2 * settings.frameConstRotation.vOmega())(Seq012, Seq012);
1587 jacobian(Seq123, Seq123) -= dmvolForceDrhov;
1588 jacobian(I4, Seq123) -= mvolForce + dmvolForceDrhov.transpose() * UMeanXy(Seq123) / UMeanXy(0);
1589 jacobian(I4, 0) -= -mvolForce.dot(UMeanXy(Seq123)) / sqr(UMeanXy(0));
1590 }
1591 }
1592#endif
1593 if (axisSymmetric)
1594 {
1595 TU uPrim;
1596 uPrim.resizeLike(UMeanXy);
1597 Gas::IdealGasThermalConservative2Primitive(UMeanXy, uPrim, settings.idealGasProperty.gamma);
1598 if (Mode == 0)
1599 ret(2) += uPrim(I4) / std::max(verySmallReal, pPhy(1)); // y direction force
1600 if (Mode == 1)
1601 ; // not implementing axisSymmetric jacobian addition
1602 if (Mode == 2)
1603 ; // not implementing axisSymmetric jacobian addition
1604 }
1605 if constexpr (Traits::isPlainNS)
1606 {
1607 }
1608 else if constexpr (Traits::hasSA)
1609 {
1610
1611 real pMean, asqrMean, Hmean;
1612 real gamma = settings.idealGasProperty.gamma;
1613 Gas::IdealGasThermal(UMeanXy(I4), UMeanXy(0), (UMeanXy(Seq123) / UMeanXy(0)).squaredNorm(),
1614 gamma, pMean, asqrMean, Hmean);
1615 // ! refvalue:
1616 real muRef = settings.idealGasProperty.muGas;
1617 real T = pMean / ((gamma - 1) / gamma * settings.idealGasProperty.CpGas * UMeanXy(0));
1618
1619 real mufPhy, muf;
1620 mufPhy = muf = muEff(UMeanXy, T);
1621
1622 real d = std::min(dWallC, std::pow(veryLargeReal, 1. / 6.));
1623 TU retInc;
1624 retInc.setZero(UMeanXy.size());
1625
1626 //! DES mesh lengths should be cached!
1627 real hMax = vfv->GetCellMaxLenScale(iCell);
1628 real cWall = settings.SADESScale > 100.0 ? 1.0 : 0.15;
1629 real lLES = hMax * settings.SADESScale;
1630 //! missing hWallNormal!
1631 lLES = std::min(lLES, std::max({d * cWall, hMax * cWall}));
1632 auto sourceCaller = [&](int mode)
1633 {
1634 RANS::GetSource_SA<dim>(UMeanXy, DiffUxy, settings.idealGasProperty.muGas, muf,
1635 gamma,
1636 d, lLES, hMax, settings.SADESMode,
1637 retInc,
1638 settings.ransSARotCorrection, mode,
1639 settings.SAVersion);
1640 };
1641
1642 if (Mode == 0)
1643 {
1644 sourceCaller(0);
1645 }
1646 else if (Mode == 1)
1647 {
1648 sourceCaller(1);
1649 }
1650 else if (Mode == 2)
1651 {
1652 sourceCaller(1);
1653 jacobian += retInc.asDiagonal(); //! TODO: make really block jacobian
1654 }
1655 ret += retInc;
1656 }
1657 else if constexpr (Traits::has2EQ)
1658 {
1659 real pMean, asqrMean, Hmean;
1660 real gamma = settings.idealGasProperty.gamma;
1661 Gas::IdealGasThermal(UMeanXy(I4), UMeanXy(0), (UMeanXy(Seq123) / UMeanXy(0)).squaredNorm(),
1662 gamma, pMean, asqrMean, Hmean);
1663 // ! refvalue:
1664 real muRef = settings.idealGasProperty.muGas;
1665 real T = pMean / ((gamma - 1) / gamma * settings.idealGasProperty.CpGas * UMeanXy(0));
1666
1667 real mufPhy, muf;
1668 mufPhy = muf = muEff(UMeanXy, T);
1669
1670 TU retInc;
1671 retInc.setZero(UMeanXy.size());
1672
1673 TU UMeanXyFixed = UMeanXy;
1674
1675 // if constexpr (Traits::has2EQ)
1676 // for (auto f : mesh->cell2face[iCell])
1677 // if (pBCHandler->GetTypeFromID(mesh->GetFaceZone(f)) == BCWall)
1678 // {
1679 // real d1 = dWallC;
1680 // real rhoOmegaaaWall = mufPhy / sqr(d1) * 800;
1681 // UMeanXyFixed(I4 + 2) = rhoOmegaaaWall;
1682 // }
1683
1684 auto sourceCaller = [&](int mode)
1685 {
1686 if (settings.ransModel == RANSModel::RANS_KOSST)
1687 RANS::GetSource_SST<dim>(UMeanXyFixed, DiffUxy, mufPhy, dWallC, vfv->GetCellMaxLenScale(iCell) * settings.SADESScale, retInc, mode);
1688 else if (settings.ransModel == RANSModel::RANS_KOWilcox)
1689 RANS::GetSource_KOWilcox<dim>(UMeanXyFixed, DiffUxy, mufPhy, dWallC, retInc, mode);
1690 else if (settings.ransModel == RANSModel::RANS_RKE)
1691 RANS::GetSource_RealizableKe<dim>(UMeanXyFixed, DiffUxy, mufPhy, dWallC, retInc, mode);
1692 };
1693
1694 if (Mode == 0)
1695 {
1696 sourceCaller(0);
1697 }
1698 else if (Mode == 1)
1699 {
1700 sourceCaller(1);
1701 }
1702 else if (Mode == 2)
1703 {
1704 sourceCaller(1);
1705 jacobian += retInc.asDiagonal(); //! TODO: make really block jacobian
1706 }
1707 ret += retInc;
1708 }
1709 else
1710 {
1711 DNDS_assert(false);
1712 }
1713 // if (Mode == 1)
1714 // std::cout << ret.transpose() << std::endl;
1715 return ret;
1716 }
1717
1718 DNDS_SWITCH_INTELLISENSE(
1719 template <EulerModel model>
1720 ,
1721 template <>)
1722 /** @brief Dispatcher for boundary ghost-value generation, routing to per-BC-type handlers.
1723 *
1724 * Selects the appropriate boundary value generator based on the BC type (farfield,
1725 * wall, inviscid wall, outflow, inflow, total-condition inflow, symmetry, etc.)
1726 * and returns the ghost (right) state for the boundary face.
1727 *
1728 * @param ULxy Left (interior) state; may be modified by some handlers (e.g., fixUL).
1729 * @param ULMeanXy Left cell mean state.
1730 * @param iCell Owner cell index.
1731 * @param iFace Face index.
1732 * @param iG Quadrature point index (-1 for face center).
1733 * @param uNorm Outward unit normal vector.
1734 * @param normBase Local coordinate frame built from the normal.
1735 * @param pPhysics Physical coordinates of the quadrature point.
1736 * @param t Current simulation time.
1737 * @param btype Boundary zone identifier.
1738 * @param fixUL If true, handler may modify ULxy (e.g., enforce zero SA at wall).
1739 * @param geomMode Geometry mode for normal evaluation.
1740 * @param linMode If 1, return linearized boundary perturbation (for Jacobian).
1741 * @return Ghost (right) state for the boundary face.
1742 */
1743 typename EulerEvaluator<model>::TU EulerEvaluator<model>::generateBoundaryValue(
1744 TU &ULxy, //! warning, possible that UL is also modified
1745 const TU &ULMeanXy,
1746 index iCell, index iFace, int iG,
1747 const TVec &uNorm,
1748 const TMat &normBase,
1749 const Geom::tPoint &pPhysics,
1750 real t,
1751 Geom::t_index btype,
1752 bool fixUL,
1753 int geomMode, int linMode)
1754 {
1755 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
1756 DNDS_assert(iG >= -2);
1757
1758 TU URxy;
1759 URxy.resizeLike(ULxy);
1760 auto bTypeEuler = pBCHandler->GetTypeFromID(btype);
1761
1762 // TODO: for all linMode == 1: implement more precise linearized BC
1763
1764 if (bTypeEuler == EulerBCType::BCSpecial ||
1765 bTypeEuler == EulerBCType::BCFar ||
1766 bTypeEuler == EulerBCType::BCOutP)
1767 {
1768 if (linMode == 1)
1769 {
1770 URxy.setZero();
1771 return URxy;
1772 }
1773 DNDS_assert(ULxy(0) > 0);
1774 if (btype == Geom::BC_ID_DEFAULT_FAR ||
1775 bTypeEuler == EulerBCType::BCFar ||
1776 bTypeEuler == EulerBCType::BCOutP)
1777 {
1778 URxy = generateBV_FarField(ULxy, ULMeanXy, iCell, iFace, iG,
1779 uNorm, normBase, pPhysics, t, btype, fixUL, geomMode);
1780 }
1781 else
1782 {
1783 URxy = generateBV_SpecialFar(ULxy, ULMeanXy, iCell, iFace, iG,
1784 uNorm, normBase, pPhysics, t, btype);
1785 }
1786 }
1787 else if (bTypeEuler == EulerBCType::BCWallInvis ||
1788 bTypeEuler == EulerBCType::BCSym)
1789 {
1790 URxy = generateBV_InviscidWall(ULxy, ULMeanXy, iCell, iFace, iG,
1791 uNorm, normBase, pPhysics, t, btype);
1792 if (linMode == 1)
1793 return URxy;
1794 }
1795 else if (bTypeEuler == EulerBCType::BCWall ||
1796 bTypeEuler == EulerBCType::BCWallIsothermal)
1797 {
1798 URxy = generateBV_ViscousWall(ULxy, ULMeanXy, iCell, iFace, iG,
1799 uNorm, normBase, pPhysics, t, btype, fixUL, linMode);
1800 if (linMode == 1)
1801 return URxy;
1802 }
1803 else if (bTypeEuler == EulerBCType::BCOut)
1804 {
1805 URxy = generateBV_Outflow(ULxy, ULMeanXy, iCell, iFace, iG,
1806 uNorm, normBase, pPhysics, t, btype);
1807 if (linMode == 1)
1808 {
1809 URxy.setZero();
1810 return URxy;
1811 }
1812 }
1813 else if (bTypeEuler == EulerBCType::BCIn)
1814 {
1815 if (linMode == 1)
1816 {
1817 URxy.setZero();
1818 return URxy;
1819 }
1820 URxy = generateBV_Inflow(ULxy, ULMeanXy, iCell, iFace, iG,
1821 uNorm, normBase, pPhysics, t, btype);
1822 }
1823 else if (bTypeEuler == EulerBCType::BCInPsTs)
1824 {
1825 if (linMode == 1)
1826 {
1827 URxy.setZero();
1828 return URxy;
1829 }
1830 URxy = generateBV_TotalConditionInflow(ULxy, ULMeanXy, iCell, iFace, iG,
1831 uNorm, normBase, pPhysics, t, btype);
1832 }
1833 else
1834 {
1835 DNDS_assert(false);
1836 }
1837 return URxy;
1838 }
1839
1840 /**************************************************************************
1841 * Per-BC-type handlers
1842 **************************************************************************/
1843
1844 DNDS_SWITCH_INTELLISENSE(
1845 template <EulerModel model>
1846 ,
1847 template <>)
1848 /** @brief Generate ghost state for farfield and back-pressure outlet boundary conditions.
1849 *
1850 * Uses Riemann-invariant-based switching: supersonic outflow extrapolates the
1851 * interior state; subsonic outflow imposes farfield pressure; subsonic inflow
1852 * imposes farfield values but retains interior pressure; supersonic inflow uses
1853 * the full farfield state. Supports anchor-point and radial-profile pressure corrections,
1854 * CL-driver AoA rotation, and rotating-frame transformations.
1855 *
1856 * @param ULxy Left (interior) state.
1857 * @param ULMeanXy Left cell mean state.
1858 * @param iCell Owner cell index.
1859 * @param iFace Face index.
1860 * @param iG Quadrature point index.
1861 * @param uNorm Outward unit normal.
1862 * @param normBase Local coordinate frame.
1863 * @param pPhysics Physical coordinates.
1864 * @param t Current simulation time.
1865 * @param btype Boundary zone identifier.
1866 * @param fixUL If true, may modify ULxy.
1867 * @param geomMode Geometry mode.
1868 * @return Ghost (right) state for the farfield face.
1869 */
1870 typename EulerEvaluator<model>::TU EulerEvaluator<model>::generateBV_FarField(
1871 TU &ULxy, const TU &ULMeanXy,
1872 index iCell, index iFace, int iG,
1873 const TVec &uNorm, const TMat &normBase,
1874 const Geom::tPoint &pPhysics, real t,
1875 Geom::t_index btype, bool fixUL, int geomMode)
1876 {
1877 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
1878 auto bTypeEuler = pBCHandler->GetTypeFromID(btype);
1879 TU URxy;
1880 URxy.resizeLike(ULxy);
1881
1882 TU far = btype >= Geom::BC_ID_DEFAULT_MAX
1883 ? pBCHandler->GetValueFromID(btype)
1884 : TU(settings.farFieldStaticValue);
1885 if (pCLDriver)
1886 far(Seq123) = pCLDriver->GetAOARotation()(Seq012, Seq012) * far(Seq123);
1887
1888 if (bTypeEuler == EulerBCType::BCFar)
1889 {
1890 if (settings.frameConstRotation.enabled && pBCHandler->GetFlagFromID(btype, "frameOpt") != 0)
1891 far(Seq123) = (Geom::RotateAxis(-settings.frameConstRotation.vOmega() * t) * Geom::ToThreeDim<dim>(far(Seq123)))(Seq012);
1892 }
1893
1894 TU ULxyStatic = ULxy;
1895 if (settings.frameConstRotation.enabled) // to static frame velocity
1896 TransformURotatingFrame(ULxyStatic, pPhysics, 1);
1897
1898 real un = ULxy(Seq123).dot(uNorm) / ULxy(0); // using relative velo for in/out judgement
1899 real gamma = settings.idealGasProperty.gamma;
1900 real asqr, H, p;
1901 Gas::IdealGasThermal(ULxyStatic(I4), ULxyStatic(0), (ULxyStatic(Seq123) / ULxyStatic(0)).squaredNorm(), gamma, p, asqr, H);
1902
1903 DNDS_assert(asqr >= 0);
1904 real a = std::sqrt(asqr);
1905
1906 auto vg = this->GetFaceVGrid(iFace, iG, pPhysics);
1907 real vgN = vg.dot(uNorm);
1908
1909 if (un - vgN - a > 0) // full outflow
1910 {
1911 URxy = ULxyStatic;
1912 }
1913 else if (un - vgN > 0) // subsonic out
1914 {
1915 TU farPrimitive, ULxyPrimitive;
1916 farPrimitive.resizeLike(ULxyStatic);
1917 ULxyPrimitive.resizeLike(URxy);
1918 Gas::IdealGasThermalConservative2Primitive<dim>(far, farPrimitive, gamma);
1919 Gas::IdealGasThermalConservative2Primitive<dim>(ULxyStatic, ULxyPrimitive, gamma);
1920 if (bTypeEuler == EulerBCType::BCOutP && pBCHandler->GetFlagFromID(btype, "anchorOpt") == 1)
1921 {
1922 {
1923 TU anchorPointRel = ULxy;
1924 if (anchorRecorders.count(btype)) // if doesn't have anchor value yet, use UL as anchor
1925 anchorPointRel = anchorRecorders.at(btype).val;
1926 TU anchorPointRelPrimitive = anchorPointRel;
1927 Gas::IdealGasThermalConservative2Primitive<dim>(anchorPointRel, anchorPointRelPrimitive, gamma);
1928 farPrimitive(I4) += std::max(ULxyPrimitive(I4) - anchorPointRelPrimitive(I4), -0.95 * farPrimitive(I4));
1929 }
1930 }
1931 if (bTypeEuler == EulerBCType::BCOutP && pBCHandler->GetFlagFromID(btype, "anchorOpt") == 2)
1932 {
1933 real pInc = 0;
1934 if (profileRecorders.count(btype))
1935 pInc = profileRecorders.at(btype).GetPlain(settings.frameConstRotation.rVec(pPhysics).norm())(I4);
1936 farPrimitive(I4) += std::max(pInc, -0.95 * farPrimitive(I4));
1937 }
1938 ULxyPrimitive(I4) = farPrimitive(I4); // using far pressure
1939 Gas::IdealGasThermalPrimitive2Conservative<dim>(ULxyPrimitive, URxy, gamma);
1940 }
1941 else if (un - vgN + a > 0) // subsonic in
1942 {
1943 TU farPrimitive, ULxyPrimitive;
1944 farPrimitive.resizeLike(ULxyStatic);
1945 ULxyPrimitive.resizeLike(URxy);
1946 Gas::IdealGasThermalConservative2Primitive<dim>(far, farPrimitive, gamma);
1947 Gas::IdealGasThermalConservative2Primitive<dim>(ULxyStatic, ULxyPrimitive, gamma);
1948 farPrimitive(I4) = ULxyPrimitive(I4); // using inner pressure
1949 Gas::IdealGasThermalPrimitive2Conservative<dim>(farPrimitive, URxy, gamma);
1950 }
1951 else // full inflow
1952 {
1953 URxy = far;
1954 }
1955 if (settings.frameConstRotation.enabled) // to rotating frame velocity
1956 TransformURotatingFrame(URxy, pPhysics, -1);
1957 return URxy;
1958 }
1959
1960 DNDS_SWITCH_INTELLISENSE(
1961 template <EulerModel model>
1962 ,
1963 template <>)
1964 /** @brief Generate ghost state for special farfield BCs (DMR, Rayleigh-Taylor, etc.).
1965 *
1966 * Handles built-in special boundary conditions identified by reserved BC IDs:
1967 * double Mach reflection (DMR), Rayleigh-Taylor instability, and user-specified
1968 * special options. Uses Riemann-invariant switching like the standard farfield.
1969 *
1970 * @param ULxy Left (interior) state.
1971 * @param ULMeanXy Left cell mean state.
1972 * @param iCell Owner cell index.
1973 * @param iFace Face index.
1974 * @param iG Quadrature point index.
1975 * @param uNorm Outward unit normal.
1976 * @param normBase Local coordinate frame.
1977 * @param pPhysics Physical coordinates.
1978 * @param t Current simulation time.
1979 * @param btype Boundary zone identifier (special reserved ID).
1980 * @return Ghost (right) state for the special farfield face.
1981 */
1982 typename EulerEvaluator<model>::TU EulerEvaluator<model>::generateBV_SpecialFar(
1983 TU &ULxy, const TU &ULMeanXy,
1984 index iCell, index iFace, int iG,
1985 const TVec &uNorm, const TMat &normBase,
1986 const Geom::tPoint &pPhysics, real t,
1987 Geom::t_index btype)
1988 {
1989 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
1990 TU URxy;
1991 URxy.resizeLike(ULxy);
1992
1993 if (btype == Geom::BC_ID_DEFAULT_SPECIAL_DMR_FAR) // DMR
1994 {
1995 DNDS_assert(dim > 1);
1996 URxy = settings.farFieldStaticValue;
1997 real uShock = 10;
1998 if constexpr (dim == 3)
1999 {
2000 if (((pPhysics(0) - uShock / std::sin(pi / 3) * t - 1. / 6.) -
2001 pPhysics(1) / std::tan(pi / 3)) > 0)
2002 URxy({0, 1, 2, 3, 4}) = Eigen::Vector<real, 5>{1.4, 0, 0, 0, 2.5};
2003 else
2004 URxy({0, 1, 2, 3, 4}) = Eigen::Vector<real, 5>{8, 57.157676649772960, -33, 0, 5.635e2};
2005 }
2006 else
2007 {
2008 if (((pPhysics(0) - uShock / std::sin(pi / 3) * t - 1. / 6.) -
2009 pPhysics(1) / std::tan(pi / 3)) > 0)
2010 URxy({0, 1, 2, 3}) = Eigen::Vector<real, 4>{1.4, 0, 0, 2.5};
2011 else
2012 URxy({0, 1, 2, 3}) = Eigen::Vector<real, 4>{8, 57.157676649772960, -33, 5.635e2};
2013 }
2014 }
2015 else if (btype == Geom::BC_ID_DEFAULT_SPECIAL_RT_FAR) // Rayleigh-Taylor
2016 {
2017 DNDS_assert(dim > 1);
2018 Eigen::VectorXd far = settings.farFieldStaticValue;
2019 real gamma = settings.idealGasProperty.gamma;
2020 real un = ULxy(Seq123).dot(uNorm) / ULxy(0);
2021 real vsqr = (ULxy(Seq123) / ULxy(0)).squaredNorm();
2022 real asqr, H, p;
2023 Gas::IdealGasThermal(ULxy(I4), ULxy(0), vsqr, gamma, p, asqr, H);
2024
2025 DNDS_assert(asqr >= 0);
2026 real a = std::sqrt(asqr);
2027 real v = -0.025 * a * cos(pPhysics(0) * 8 * pi);
2028
2029 if (pPhysics(1) < 0.5)
2030 {
2031 real rho = 2;
2032 real p = 1;
2033 far(0) = rho;
2034 far(1) = 0;
2035 far(2) = rho * v;
2036 far(I4) = 0.5 * rho * sqr(v) + p / (gamma - 1);
2037 }
2038 else
2039 {
2040 real rho = 1;
2041 real p = 2.5;
2042 far(0) = rho;
2043 far(1) = 0;
2044 far(2) = rho * v;
2045 far(I4) = 0.5 * rho * sqr(v) + p / (gamma - 1);
2046 }
2047
2048 if (un - a > 0)
2049 URxy = ULxy;
2050 else if (un > 0)
2051 {
2052 TU farPrimitive, ULxyPrimitive;
2053 farPrimitive.resizeLike(ULxy);
2054 ULxyPrimitive.resizeLike(URxy);
2055 Gas::IdealGasThermalConservative2Primitive<dim>(far, farPrimitive, gamma);
2056 Gas::IdealGasThermalConservative2Primitive<dim>(ULxy, ULxyPrimitive, gamma);
2057 ULxyPrimitive(I4) = farPrimitive(I4);
2058 Gas::IdealGasThermalPrimitive2Conservative<dim>(ULxyPrimitive, URxy, gamma);
2059 }
2060 else if (un + a > 0)
2061 {
2062 TU farPrimitive, ULxyPrimitive;
2063 farPrimitive.resizeLike(ULxy);
2064 ULxyPrimitive.resizeLike(URxy);
2065 Gas::IdealGasThermalConservative2Primitive<dim>(far, farPrimitive, gamma);
2066 Gas::IdealGasThermalConservative2Primitive<dim>(ULxy, ULxyPrimitive, gamma);
2067 farPrimitive(I4) = ULxyPrimitive(I4);
2068 Gas::IdealGasThermalPrimitive2Conservative<dim>(farPrimitive, URxy, gamma);
2069 }
2070 else
2071 URxy = far;
2072 }
2073 else if (btype == Geom::BC_ID_DEFAULT_SPECIAL_IV_FAR) // Isentropic Vortex
2074 {
2075 real chi = 5;
2076 real gamma = settings.idealGasProperty.gamma;
2077 real xc = 5 + t;
2078 real yc = 5 + t;
2079 real r = std::sqrt(sqr(pPhysics(0) - xc) + sqr(pPhysics(1) - yc));
2080 real dT = -(gamma - 1) / (8 * gamma * sqr(pi)) * sqr(chi) * std::exp(1 - sqr(r));
2081 real dux = chi / 2 / pi * std::exp((1 - sqr(r)) / 2) * -(pPhysics(1) - xc);
2082 real duy = chi / 2 / pi * std::exp((1 - sqr(r)) / 2) * +(pPhysics(0) - yc);
2083 real T = dT + 1;
2084 real ux = dux + 1;
2085 real uy = duy + 1;
2086 real S = 1;
2087 real rho = std::pow(T / S, 1 / (gamma - 1));
2088 real pVal = T * rho;
2089 real E = 0.5 * (sqr(ux) + sqr(uy)) * rho + pVal / (gamma - 1);
2090
2091 URxy.setZero();
2092 URxy(0) = rho;
2093 URxy(1) = rho * ux;
2094 URxy(2) = rho * uy;
2095 URxy(dim + 1) = E;
2096 }
2097 else if (btype == Geom::BC_ID_DEFAULT_SPECIAL_2DRiemann_FAR) // 2D Riemann
2098 {
2099 real gamma = settings.idealGasProperty.gamma;
2100 real bdL = 0.0, bdR = 1.0, bdD = 0.0, bdU = 1.0;
2101
2102 real phi1 = -0.663324958071080;
2103 real phi2 = -0.422115882408869;
2104 real location = 0.8;
2105 real p1 = location + phi1 * t;
2106 real p2 = location + phi2 * t;
2107 real rho, u, v, pre;
2108 TU ULxyPrimitive;
2109 ULxyPrimitive.resizeLike(ULxy);
2110
2111 Gas::IdealGasThermalConservative2Primitive<dim>(ULxy, ULxyPrimitive, gamma);
2112 real rhoL = ULxyPrimitive(0);
2113 real uL = ULxyPrimitive(1);
2114 real vL = ULxyPrimitive(2);
2115 real preL = ULxyPrimitive(I4);
2116 TU farPrimitive = ULxyPrimitive;
2117
2118 static const real bTol = 1e-9;
2119 if (std::abs(pPhysics(0) - bdL) < bTol)
2120 { // left, phi2
2121 if (pPhysics(1) <= p2)
2122 rho = 0.137992831541219, u = 1.206045378311055, v = 1.206045378311055, pre = 0.029032258064516;
2123 else
2124 rho = 0.532258064516129, u = 1.206045378311055, v = 0.0, pre = 0.3;
2125 }
2126 else if (std::abs(pPhysics(0) - bdR) < bTol)
2127 { // right, phi1
2128 if (pPhysics(1) <= p1)
2129 rho = rhoL, u = -uL, v = vL, pre = preL;
2130 else
2131 rho = rhoL, u = -uL, v = vL, pre = preL;
2132 }
2133 else if (std::abs(pPhysics(1) - bdU) < bTol)
2134 { // up, phi1
2135 if (pPhysics(0) <= p1)
2136 rho = rhoL, u = uL, v = -vL, pre = preL;
2137 else
2138 rho = rhoL, u = uL, v = -vL, pre = preL;
2139 }
2140 else if (std::abs(pPhysics(1) - bdD) < bTol)
2141 { // down, phi2
2142 if (pPhysics(0) <= p2)
2143 rho = 0.137992831541219, u = 1.206045378311055, v = 1.206045378311055, pre = 0.029032258064516;
2144 else
2145 rho = 0.532258064516129, u = 0.0, v = 1.206045378311055, pre = 0.3;
2146 }
2147 else
2148 {
2149 rho = u = v = pre = std::nan("1");
2150 DNDS_assert(false);
2151 }
2152 farPrimitive(0) = rho;
2153 farPrimitive(1) = u, farPrimitive(2) = v;
2154 farPrimitive(I4) = pre;
2155 Gas::IdealGasThermalPrimitive2Conservative<dim>(farPrimitive, URxy, gamma);
2156 }
2157 else if (pBCHandler->GetFlagFromID(btype, "specialOpt") == 3001) // Noh
2158 {
2159 TU farPrimitive;
2160 Gas::IdealGasThermalConservative2Primitive<dim>(settings.farFieldStaticValue, farPrimitive, settings.idealGasProperty.gamma);
2161 real pInf = farPrimitive(I4);
2162 real r = pPhysics.norm();
2163 TVec velo = -pPhysics(Seq012) / (r + smallReal);
2164 real rho = sqr(1. + t / (r + smallReal));
2165 farPrimitive(0) = rho;
2166 farPrimitive(Seq123) = velo;
2167 farPrimitive(I4) = pInf;
2168 Gas::IdealGasThermalPrimitive2Conservative<dim>(farPrimitive, URxy, settings.idealGasProperty.gamma);
2169 }
2170 else
2171 DNDS_assert_info(false, fmt::format(
2172 "btype [{}] or bTypeEuler [{}] or specialOpt [{}] is not supported",
2173 btype, to_string(pBCHandler->GetTypeFromID(btype)),
2174 pBCHandler->GetFlagFromIDSoft(btype, "specialOpt")));
2175 return URxy;
2176 }
2177
2178 DNDS_SWITCH_INTELLISENSE(
2179 template <EulerModel model>
2180 ,
2181 template <>)
2182 /** @brief Generate ghost state for inviscid wall (slip wall) and symmetry boundary conditions.
2183 *
2184 * Mirrors the velocity component normal to the wall while preserving the
2185 * tangential component, density, and energy. Handles rotating-frame transformations.
2186 *
2187 * @param ULxy Left (interior) state.
2188 * @param ULMeanXy Left cell mean state.
2189 * @param iCell Owner cell index.
2190 * @param iFace Face index.
2191 * @param iG Quadrature point index.
2192 * @param uNorm Outward unit normal.
2193 * @param normBase Local coordinate frame.
2194 * @param pPhysics Physical coordinates.
2195 * @param t Current simulation time.
2196 * @param btype Boundary zone identifier.
2197 * @return Ghost (right) state with reflected normal velocity.
2198 */
2199 typename EulerEvaluator<model>::TU EulerEvaluator<model>::generateBV_InviscidWall(
2200 TU &ULxy, const TU &ULMeanXy,
2201 index iCell, index iFace, int iG,
2202 const TVec &uNorm, const TMat &normBase,
2203 const Geom::tPoint &pPhysics, real t,
2204 Geom::t_index btype)
2205 {
2206 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
2207 TU URxy = ULxy;
2208 if (settings.frameConstRotation.enabled)
2209 this->TransformURotatingFrame_ABS_VELO(URxy, pPhysics, -1);
2210 URxy(Seq123) -= 2 * URxy(Seq123).dot(uNorm) * uNorm; // mirrored!
2211 if (settings.frameConstRotation.enabled)
2212 this->TransformURotatingFrame_ABS_VELO(URxy, pPhysics, 1);
2213 return URxy;
2214 }
2215
2216 DNDS_SWITCH_INTELLISENSE(
2217 template <EulerModel model>
2218 ,
2219 template <>)
2220 /** @brief Generate ghost state for viscous (no-slip) wall boundary conditions.
2221 *
2222 * Reverses the velocity to enforce zero velocity at the wall (ghost mirroring).
2223 * For isothermal walls, adjusts density to match the prescribed wall temperature.
2224 * Handles SA turbulence variable (negated or zeroed at wall) and k-omega/k-epsilon
2225 * wall boundary conditions including the realizable k-epsilon wall epsilon formula.
2226 *
2227 * @param ULxy Left (interior) state (may be modified when fixUL is true for SA).
2228 * @param ULMeanXy Left cell mean state.
2229 * @param iCell Owner cell index.
2230 * @param iFace Face index.
2231 * @param iG Quadrature point index.
2232 * @param uNorm Outward unit normal.
2233 * @param normBase Local coordinate frame.
2234 * @param pPhysics Physical coordinates.
2235 * @param t Current simulation time.
2236 * @param btype Boundary zone identifier.
2237 * @param fixUL If true, may zero SA variables in ULxy at the wall.
2238 * @param linMode If 1, return linearized ghost perturbation for Jacobian.
2239 * @return Ghost (right) state with no-slip wall condition.
2240 */
2241 typename EulerEvaluator<model>::TU EulerEvaluator<model>::generateBV_ViscousWall(
2242 TU &ULxy, const TU &ULMeanXy,
2243 index iCell, index iFace, int iG,
2244 const TVec &uNorm, const TMat &normBase,
2245 const Geom::tPoint &pPhysics, real t,
2246 Geom::t_index btype, bool fixUL, int linMode)
2247 {
2248 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
2249 auto bTypeEuler = pBCHandler->GetTypeFromID(btype);
2250 TU URxy = ULxy;
2251 if (true) // physical mode
2252 {
2253 if (settings.frameConstRotation.enabled && pBCHandler->GetFlagFromID(btype, "frameOpt") == 0)
2254 this->TransformURotatingFrame_ABS_VELO(URxy, pPhysics, -1);
2255 if (settings.frameConstRotation.enabled && pBCHandler->GetFlagFromID(btype, "frameOpt") != 0)
2256 this->TransformURotatingFrame(URxy, pPhysics, 1);
2257 URxy(Seq123) *= -1;
2258 if (settings.frameConstRotation.enabled && pBCHandler->GetFlagFromID(btype, "frameOpt") == 0)
2259 this->TransformURotatingFrame_ABS_VELO(URxy, pPhysics, 1);
2260 if (settings.frameConstRotation.enabled && pBCHandler->GetFlagFromID(btype, "frameOpt") != 0)
2261 this->TransformURotatingFrame(URxy, pPhysics, -1);
2262 }
2263
2264 if (linMode == 1)
2265 return URxy;
2266
2267 if (bTypeEuler == EulerBCType::BCWallIsothermal)
2268 {
2269 real temp = pBCHandler->GetValueFromID(btype)(0);
2270 TU URxyPrim;
2271 URxyPrim.resizeLike(ULxy);
2272 Gas::IdealGasThermalConservative2Primitive<dim>(URxy, URxyPrim, settings.idealGasProperty.gamma);
2273 DNDS_assert_info(URxyPrim(0) > 0 && URxyPrim(I4) > 0 && temp > 0, fmt::format("{}, {}, {}", URxyPrim(0), URxyPrim(I4), temp));
2274 real newDensity = URxyPrim(I4) / temp / settings.idealGasProperty.Rgas;
2275 URxyPrim(0) = newDensity;
2276 Gas::IdealGasThermalPrimitive2Conservative<dim>(URxyPrim, URxy, settings.idealGasProperty.gamma);
2277 }
2278 if constexpr (Traits::hasSA)
2279 {
2280 URxy(I4 + 1) *= -1;
2281#ifdef USE_FIX_ZERO_SA_NUT_AT_WALL
2282 if (fixUL)
2283 ULxy(I4 + 1) = URxy(I4 + 1) = 0; //! modifing UL
2284#endif
2285 }
2286 if constexpr (Traits::has2EQ)
2287 {
2288 URxy({I4 + 1, I4 + 2}) *= -1;
2289#ifdef USE_FIX_ZERO_SA_NUT_AT_WALL
2290 // if (fixUL)
2291 // ULxy({I4 + 1, I4 + 2}).setZero(), URxy({I4 + 1, I4 + 2}).setZero();
2292#endif
2293 if (settings.ransModel == RANSModel::RANS_RKE)
2294 { // BC for RealizableKe
2295 real d1 = dWall[iCell].mean();
2296 real k1 = ULMeanXy(I4 + 1) / ULMeanXy(0);
2297
2298 real pMean, asqrMean, Hmean;
2299 real gamma = settings.idealGasProperty.gamma;
2300 Gas::IdealGasThermal(ULMeanXy(I4), ULMeanXy(0), (ULMeanXy(Seq123) / ULMeanXy(0)).squaredNorm(),
2301 gamma, pMean, asqrMean, Hmean);
2302 real T = pMean / ((gamma - 1) / gamma * settings.idealGasProperty.CpGas * ULMeanXy(0));
2303 real mufPhy1 = muEff(ULMeanXy, T);
2304 real epsWall = 2 * mufPhy1 / ULMeanXy(0) * k1 / sqr(d1);
2305 URxy(I4 + 2) = 2 * epsWall * ULxy(0) - ULxy(I4 + 2);
2306 }
2307 if (settings.ransModel == RANSModel::RANS_KOSST ||
2308 settings.ransModel == RANSModel::RANS_KOWilcox)
2309 { // BC for SST or KOWilcox
2310 real d1 = dWall[iCell].mean();
2311 real pMean, asqrMean, Hmean;
2312 real gamma = settings.idealGasProperty.gamma;
2313 Gas::IdealGasThermal(ULMeanXy(I4), ULMeanXy(0), (ULMeanXy(Seq123) / ULMeanXy(0)).squaredNorm(),
2314 gamma, pMean, asqrMean, Hmean);
2315 real T = pMean / ((gamma - 1) / gamma * settings.idealGasProperty.CpGas * ULMeanXy(0));
2316 real mufPhy1 = muEff(ULMeanXy, T);
2317
2318 real rhoOmegaaaWall = mufPhy1 / sqr(d1) * RANS::kWallOmegaCoeff;
2319 URxy(I4 + 2) = 2 * rhoOmegaaaWall - ULxy(I4 + 2);
2320 }
2321 }
2322 return URxy;
2323 }
2324
2325 DNDS_SWITCH_INTELLISENSE(
2326 template <EulerModel model>
2327 ,
2328 template <>)
2329 /** @brief Generate ghost state for supersonic/extrapolation outflow BC.
2330 *
2331 * Simply returns the interior state as the ghost value, allowing all waves to
2332 * exit the domain without reflection.
2333 *
2334 * @param ULxy Left (interior) state.
2335 * @param ULMeanXy Left cell mean state.
2336 * @param iCell Owner cell index.
2337 * @param iFace Face index.
2338 * @param iG Quadrature point index.
2339 * @param uNorm Outward unit normal.
2340 * @param normBase Local coordinate frame.
2341 * @param pPhysics Physical coordinates.
2342 * @param t Current simulation time.
2343 * @param btype Boundary zone identifier.
2344 * @return Ghost state equal to the interior state (full extrapolation).
2345 */
2346 typename EulerEvaluator<model>::TU EulerEvaluator<model>::generateBV_Outflow(
2347 TU &ULxy, const TU &ULMeanXy,
2348 index iCell, index iFace, int iG,
2349 const TVec &uNorm, const TMat &normBase,
2350 const Geom::tPoint &pPhysics, real t,
2351 Geom::t_index btype)
2352 {
2353 return ULxy;
2354 }
2355
2356 DNDS_SWITCH_INTELLISENSE(
2357 template <EulerModel model>
2358 ,
2359 template <>)
2360 /** @brief Generate ghost state for prescribed inflow boundary condition.
2361 *
2362 * Sets the ghost state to the BC-prescribed conservative values. Applies
2363 * CL-driver AoA rotation and rotating-frame transformation if applicable.
2364 *
2365 * @param ULxy Left (interior) state (unused for pure inflow).
2366 * @param ULMeanXy Left cell mean state.
2367 * @param iCell Owner cell index.
2368 * @param iFace Face index.
2369 * @param iG Quadrature point index.
2370 * @param uNorm Outward unit normal.
2371 * @param normBase Local coordinate frame.
2372 * @param pPhysics Physical coordinates.
2373 * @param t Current simulation time.
2374 * @param btype Boundary zone identifier.
2375 * @return Ghost state set to prescribed inflow values.
2376 */
2377 typename EulerEvaluator<model>::TU EulerEvaluator<model>::generateBV_Inflow(
2378 TU &ULxy, const TU &ULMeanXy,
2379 index iCell, index iFace, int iG,
2380 const TVec &uNorm, const TMat &normBase,
2381 const Geom::tPoint &pPhysics, real t,
2382 Geom::t_index btype)
2383 {
2384 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
2385 TU URxy = pBCHandler->GetValueFromID(btype);
2386 if (pCLDriver)
2387 URxy(Seq123) = pCLDriver->GetAOARotation()(Seq012, Seq012) * URxy(Seq123);
2388 // Note: removed dead code that checked bTypeEuler == BCFar (unreachable in BCIn branch)
2389 if (settings.frameConstRotation.enabled)
2390 TransformURotatingFrame(URxy, pPhysics, -1);
2391 return URxy;
2392 }
2393
2394 DNDS_SWITCH_INTELLISENSE(
2395 template <EulerModel model>
2396 ,
2397 template <>)
2398 /** @brief Generate ghost state for total-condition (stagnation) inflow BC.
2399 *
2400 * Given prescribed stagnation pressure and temperature, computes the static
2401 * conditions from the interior velocity magnitude using isentropic relations.
2402 * The flow direction is taken from the BC-prescribed direction vector.
2403 * Applies CL-driver AoA rotation and rotating-frame transformation.
2404 *
2405 * @param ULxy Left (interior) state (used for velocity magnitude).
2406 * @param ULMeanXy Left cell mean state.
2407 * @param iCell Owner cell index.
2408 * @param iFace Face index.
2409 * @param iG Quadrature point index.
2410 * @param uNorm Outward unit normal.
2411 * @param normBase Local coordinate frame.
2412 * @param pPhysics Physical coordinates.
2413 * @param t Current simulation time.
2414 * @param btype Boundary zone identifier.
2415 * @return Ghost state computed from total conditions and interior velocity.
2416 */
2417 typename EulerEvaluator<model>::TU EulerEvaluator<model>::generateBV_TotalConditionInflow(
2418 TU &ULxy, const TU &ULMeanXy,
2419 index iCell, index iFace, int iG,
2420 const TVec &uNorm, const TMat &normBase,
2421 const Geom::tPoint &pPhysics, real t,
2422 Geom::t_index btype)
2423 {
2424 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
2425 TU URxy;
2426 URxy.resizeLike(ULxy);
2427
2428 TU ULxyStatic = ULxy;
2429 if (settings.frameConstRotation.enabled)
2430 TransformURotatingFrame(ULxyStatic, pPhysics, 1);
2431 TU ULxyPrimitive;
2432 ULxyPrimitive.resizeLike(ULxy);
2433 real gamma = settings.idealGasProperty.gamma;
2434 Gas::IdealGasThermalConservative2Primitive<dim>(ULxyStatic, ULxyPrimitive, gamma);
2435 TVec v = ULxyStatic(Seq123).array() / ULxyStatic(0);
2436 real vSqr = v.squaredNorm();
2437 {
2438 TU farPrimitive = pBCHandler->GetValueFromID(btype); // primitive passive scalar components like Nu
2439
2440 real pStag = pBCHandler->GetValueFromID(btype)(0);
2441 real tStag = pBCHandler->GetValueFromID(btype)(1);
2442 vSqr = std::min(vSqr, tStag * 2 * settings.idealGasProperty.CpGas * 0.95);
2443 real tStatic = tStag - 0.5 * vSqr / (settings.idealGasProperty.CpGas);
2444 real pStatic = pStag * std::pow(tStatic / tStag, gamma / (gamma - 1));
2445 real rStatic = pStatic / (settings.idealGasProperty.Rgas * tStatic);
2446 farPrimitive(0) = rStatic;
2447 farPrimitive(Seq123) = pBCHandler->GetValueFromID(btype)(Seq234).normalized() * std::sqrt(vSqr);
2448 farPrimitive(I4) = pStatic;
2449 Gas::IdealGasThermalPrimitive2Conservative<dim>(farPrimitive, URxy, gamma);
2450 }
2451 if (pCLDriver)
2452 URxy(Seq123) = pCLDriver->GetAOARotation()(Seq012, Seq012) * URxy(Seq123);
2453 if (settings.frameConstRotation.enabled)
2454 TransformURotatingFrame(URxy, pPhysics, -1);
2455 return URxy;
2456 }
2457
2458 template <EulerModel model>
2459 /** @brief Register available output scalar fields in the OutputPicker.
2460 *
2461 * Populates the output picker with named lambda functions that extract per-cell
2462 * scalar quantities (conservative components, reconstruction coefficients,
2463 * limiter values, wall distance, cell volume, turbulent viscosity, etc.)
2464 * for use in post-processing output.
2465 *
2466 * @param op OutputPicker to populate with named field extractors (output).
2467 * @param dataRefs References to the solution arrays (u, uRec, betaPP, alphaPP).
2468 */
2469 void EulerEvaluator<model>::InitializeOutputPicker(OutputPicker &op, OutputOverlapDataRefs dataRefs)
2470 {
2471 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
2472
2473 auto &eval = *this;
2474 auto &u = dataRefs.u;
2475 auto &uRec = dataRefs.uRec;
2476 auto &betaPP = dataRefs.betaPP;
2477 auto &alphaPP = dataRefs.alphaPP;
2478
2479 OutputPicker::tMap outMap;
2480 // outMap["R"] = [&](index iCell)
2481 // { return u[iCell](0); };
2482 outMap["RU"] = [&](index iCell)
2483 { return u[iCell](1); };
2484 outMap["RV"] = [&](index iCell)
2485 { return u[iCell](2); };
2486 outMap["RV"] = [&](index iCell)
2487 { return u[iCell](I4 - 1); };
2488 outMap["RE"] = [&](index iCell)
2489 { return u[iCell](I4); };
2490 outMap["R_REC_1"] = [&](index iCell)
2491 { return uRec[iCell](0, 0); };
2492 outMap["RU_REC_1"] = [&](index iCell)
2493 { return uRec[iCell](1, 0); }; // TODO: to be continued to...
2494
2495 // pps:
2496 outMap["betaPP"] = [&](index iCell)
2497 { return betaPP[iCell](0); };
2498 outMap["alphaPP"] = [&](index iCell)
2499 { return alphaPP[iCell](0); };
2500 outMap["ACond"] = [&](index iCell)
2501 {
2502 auto AI = vfv->GetCellRecMatAInv(iCell);
2503 Eigen::MatrixXd AIInv = AI;
2504 real aCond = HardEigen::EigenLeastSquareInverse(AI, AIInv);
2505 return aCond;
2506 };
2507 outMap["dWall"] = [&](index iCell)
2508 {
2509 return eval.dWall.at(iCell).mean();
2510 };
2511 outMap["minJacobiDetRel"] = [&](index iCell)
2512 {
2513 auto eCell = mesh->GetCellElement(iCell);
2514 auto qCell = vfv->GetCellQuad(iCell);
2515 real minDetJac = veryLargeReal;
2516 for (int iG = 0; iG < qCell.GetNumPoints(); iG++)
2517 minDetJac = std::min(vfv->GetCellJacobiDet(iCell, iG), minDetJac);
2518 return minDetJac * Geom::Elem::ParamSpaceVol(eCell.GetParamSpace()) / vfv->GetCellVol(iCell);
2519 };
2520 outMap["cellVolume"] = [&](index iCell)
2521 {
2522 return vfv->GetCellVol(iCell);
2523 };
2524 outMap["mut"] = [&](index iCell)
2525 {
2526 real mut = 0;
2527 if (model == NS_2EQ || model == NS_2EQ_3D)
2528 {
2529 TU Uxy = u[iCell];
2530 TDiffU GradU;
2531 GradU.resize(Eigen::NoChange, nVars);
2532 GradU.setZero();
2533 if constexpr (gDim == 2)
2534 GradU({0, 1}, EigenAll) =
2535 vfv->GetIntPointDiffBaseValue(iCell, -1, -1, -1, std::array<int, 2>{1, 2}, 3) *
2536 uRec[iCell]; // 2d specific
2537 else
2538 GradU({0, 1, 2}, EigenAll) =
2539 vfv->GetIntPointDiffBaseValue(iCell, -1, -1, -1, std::array<int, 3>{1, 2, 3}, 4) *
2540 uRec[iCell]; // 3d specific
2541 real pMean, asqrMean, Hmean;
2542 real gamma = settings.idealGasProperty.gamma;
2543 auto ULMeanXy = Uxy;
2544 Gas::IdealGasThermal(ULMeanXy(I4), ULMeanXy(0), (ULMeanXy(Seq123) / ULMeanXy(0)).squaredNorm(),
2545 gamma, pMean, asqrMean, Hmean);
2546 // ! refvalue:
2547 real muRef = settings.idealGasProperty.muGas;
2548 real T = pMean / ((gamma - 1) / gamma * settings.idealGasProperty.CpGas * ULMeanXy(0));
2549 real mufPhy = muEff(ULMeanXy, T);
2550 if (settings.ransModel == RANSModel::RANS_KOSST)
2551 mut = RANS::GetMut_SST<dim>(Uxy, GradU, mufPhy, dWall[iCell].mean());
2552 else if (settings.ransModel == RANSModel::RANS_KOWilcox)
2553 mut = RANS::GetMut_KOWilcox<dim>(Uxy, GradU, mufPhy, dWall[iCell].mean());
2554 else if (settings.ransModel == RANSModel::RANS_RKE)
2555 mut = RANS::GetMut_RealizableKe<dim>(Uxy, GradU, mufPhy, dWall[iCell].mean());
2556 }
2557
2558 return mut;
2559 };
2560
2561 op.setMap(outMap);
2562 }
2563}
DNDS_MAKE_SSP
#define DNDS_MAKE_SSP(ssp,...)
Definition Defines.hpp:217
DNDS_RESTRICT
#define DNDS_RESTRICT
Definition Defines.hpp:982
DNDS_assert_info
#define DNDS_assert_info(expr, info)
Debug-only assertion with an extra std::string info message.
Definition Errors.hpp:117
DNDS_assert
#define DNDS_assert(expr)
Debug-only assertion (compiled out when DNDS_NDEBUG is defined). Prints the expression + file/line + ...
Definition Errors.hpp:112
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
DNDS_MPI_InsertCheck
#define DNDS_MPI_InsertCheck(mpi, info)
Debug helper: barrier + print "CHECK <info> RANK r @ fn @ file:line".
Definition MPI.hpp:343
DNDS_SWITCH_INTELLISENSE
#define DNDS_SWITCH_INTELLISENSE(real, intellisense)
Definition Macros.hpp:86
RANS_ke.hpp
RANS two-equation turbulence model implementations for the DNDSR CFD solver.
DNDS::Euler::EulerEvaluator::GetWallDist
void GetWallDist()
Compute wall distance for all cells (dispatcher).
DNDS::Euler::EulerEvaluator::InitializeOutputPicker
void InitializeOutputPicker(OutputPicker &op, OutputOverlapDataRefs dataRefs)
Initialize an OutputPicker with field callbacks for VTK/HDF5 output.
Definition EulerEvaluator_EvaluateDt.hxx:2469
DNDS::Euler::EulerEvaluator::TU
Eigen::VectorFMTSafe< real, nVarsFixed > TU
Conservative variable vector (nVarsFixed components).
Definition EulerEvaluator.hpp:94
DNDS::Euler::EulerEvaluator::EvaluateDt
void EvaluateDt(ArrayDOFV< 1 > &dt, ArrayDOFV< nVarsFixed > &u, ArrayRECV< nVarsFixed > &uRec, real CFL, real &dtMinall, real MaxDt, bool UseLocaldt, real t, uint64_t flags=DT_No_Flags)
Estimate the local or global time step for each cell.
DNDS::Euler::EulerEvaluator::generateBoundaryValue
TU generateBoundaryValue(TU &ULxy, const TU &ULMeanXy, index iCell, index iFace, int iG, const TVec &uNorm, const TMat &normBase, const Geom::tPoint &pPhysics, real t, Geom::t_index btype, bool fixUL=false, int geomMode=0, int linMode=0)
Generate the ghost (boundary) state for a boundary face.
DNDS::Euler::EulerEvaluator::source
TU source(const TU &UMeanXy, const TDiffU &DiffUxy, const Geom::tPoint &pPhy, TJacobianU &jacobian, index iCell, index ig, int Mode)
Compute the source term at a cell quadrature point.
DNDS::Euler::EulerEvaluator::fluxFace
void fluxFace(const TU_Batch &ULxy, const TU_Batch &URxy, const TU &ULMeanXy, const TU &URMeanXy, const TDiffU_Batch &DiffUxy, const TDiffU_Batch &DiffUxyPrim, const TVec_Batch &unitNorm, const TVec_Batch &vg, const TVec &unitNormC, const TVec &vgC, TU_Batch &FLfix, TU_Batch &FRfix, TU_Batch &finc, TReal_Batch &lam0V, TReal_Batch &lam123V, TReal_Batch &lam4V, Geom::t_index btype, typename Gas::RiemannSolverType rsType, index iFace, bool ignoreVis)
Compute the numerical flux at a face (batched over quadrature points).
DNDS::Euler::OutputPicker
Registry that maps named cell-scalar output fields to getter lambdas.
Definition EulerBC.hpp:716
DNDS::Euler::OutputPicker::tMap
std::map< std::string, tCellScalarFGet > tMap
Name-to-getter mapping type.
Definition EulerBC.hpp:718
DNDS::Euler::Gas::RiemannSolverType
RiemannSolverType
Selects the approximate Riemann solver and its entropy-fix variant.
Definition Gas.hpp:62
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::Gas::IdealGasThermalConservative2Primitive
void IdealGasThermalConservative2Primitive(const TCons &U, TPrim &prim, real gamma)
Converts conservative variables to primitive variables for an ideal gas.
Definition Gas.hpp:279
DNDS::Euler::RANS_KOSST
@ RANS_KOSST
Menter k-omega SST two-equation model.
Definition Euler.hpp:903
DNDS::Euler::RANS_RKE
@ RANS_RKE
Realizable k-epsilon two-equation model.
Definition Euler.hpp:904
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
@ 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::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::Euler::BCWallInvis
@ BCWallInvis
Inviscid slip wall boundary.
Definition EulerBC.hpp:40
DNDS::Euler::BCSpecial
@ BCSpecial
Test-case-specific boundary (Riemann, DMR, isentropic vortex, RT).
Definition EulerBC.hpp:47
DNDS::Euler::BCOut
@ BCOut
Supersonic outflow (extrapolation) boundary.
Definition EulerBC.hpp:42
DNDS::Euler::BCUnknown
@ BCUnknown
Uninitialized / invalid sentinel.
Definition EulerBC.hpp:37
DNDS::Euler::BCInPsTs
@ BCInPsTs
Subsonic inflow with prescribed total pressure and temperature.
Definition EulerBC.hpp:45
DNDS::Euler::BCSym
@ BCSym
Symmetry plane boundary.
Definition EulerBC.hpp:46
DNDS::Euler::BCFar
@ BCFar
Far-field (characteristic-based) boundary.
Definition EulerBC.hpp:38
DNDS::Euler::BCOutP
@ BCOutP
Back-pressure (subsonic) outflow boundary.
Definition EulerBC.hpp:43
DNDS::Euler::BCIn
@ BCIn
Supersonic inflow (fully prescribed state) boundary.
Definition EulerBC.hpp:44
DNDS::Geom::Elem::ParamSpaceVol
DNDS_DEVICE_CALLABLE constexpr t_real ParamSpaceVol(ParamSpace ps)
Definition Elements.hpp:101
DNDS::Geom::Elem::GetQuadPatches
auto GetQuadPatches(Quadrature &q)
Definition Quadrature.hpp:186
DNDS::Geom::t_index
int32_t t_index
Definition Geometric.hpp:6
DNDS::Geom::FaceIDIsInternal
DNDS_DEVICE_CALLABLE bool FaceIDIsInternal(t_index id)
Definition BoundaryCondition.hpp:85
DNDS::Geom::RotateAxis
tGPoint RotateAxis(const tPoint &axis)
Definition Geometric.hpp:183
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::Geom::CrossVecToMat
tGPoint CrossVecToMat(const tPoint &axn)
Definition Geometric.hpp:195
DNDS::HardEigen::EigenLeastSquareInverse
real EigenLeastSquareInverse(const Eigen::MatrixXd &A, Eigen::MatrixXd &AI, real svdTol)
Moore-Penrose pseudoinverse via SVD, dropping singular values below svdTol (relative to the largest)....
Definition HardEigen.cpp:13
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:202
DNDS::MPI::Barrier
MPI_int Barrier(MPI_Comm comm)
Wrapper over MPI_Barrier.
Definition MPI.cpp:247
DNDS::MPI::AllreduceOneIndex
void AllreduceOneIndex(index &v, MPI_Op op, const MPIInfo &mpi)
Single-scalar Allreduce helper for indices (in-place, count = 1).
Definition MPI.hpp:727
DNDS::UnInitIndex
DNDS_CONSTANT const index UnInitIndex
Sentinel "not initialised" index value (= INT64_MIN).
Definition Defines.hpp:181
DNDS::MPISerialDo
void MPISerialDo(const MPIInfo &mpi, F f)
Execute f on each rank serially, in rank order.
Definition MPI.hpp:749
DNDS::verySmallReal
DNDS_CONSTANT const real verySmallReal
Catch-all lower bound ("effectively zero").
Definition Defines.hpp:199
DNDS::pi
DNDS_CONSTANT const real pi
π in double precision (matches DNDS_E_PI macro).
Definition Defines.hpp:204
DNDS::sign
constexpr real sign(real a)
Signum function: +1, 0, or -1.
Definition Defines.hpp:537
DNDS::MatrixXR
Eigen::Matrix< real, Eigen::Dynamic, Eigen::Dynamic > MatrixXR
Column-major dynamic Eigen matrix of reals (default layout).
Definition Defines.hpp:210
DNDS::index
int64_t index
Global row / DOF index type (signed 64-bit; handles multi-billion-cell meshes).
Definition Defines.hpp:112
DNDS::sqr
constexpr T sqr(const T &a)
a * a, constexpr. Works for all arithmetic types.
Definition Defines.hpp:522
DNDS::real
double real
Canonical floating-point scalar used throughout DNDSR (double precision).
Definition Defines.hpp:110
DNDS::smallReal
DNDS_CONSTANT const real smallReal
Loose lower bound (for iterative-solver tolerances etc.).
Definition Defines.hpp:201
DNDS::log
std::ostream & log()
Return the current DNDSR log stream (either std::cout or the installed file).
Definition Defines.cpp:50
DNDS::to_string
std::string to_string(const Eigen::DenseBase< dir > &v, int precision=5, bool scientific=false)
Render an Eigen::DenseBase to a string via operator<<.
Definition EigenUtil.hpp:29
DNDS::veryLargeReal
DNDS_CONSTANT const real veryLargeReal
Catch-all upper bound ("practically infinity") for physical scalars.
Definition Defines.hpp:195
DNDS::DNDS_MPI_REAL
const MPI_Datatype DNDS_MPI_REAL
MPI datatype matching real (= MPI_REAL8).
Definition MPI.hpp:108
DNDS::ArrayTransformerType::Type
ArrayTransformer< typename TArray::value_type, TArray::rs, TArray::rm, TArray::al > Type
Definition ArrayTransformer.hpp:1490
DNDS::Euler::EulerEvaluator::OutputOverlapDataRefs
References to arrays needed by the output data picker.
Definition EulerEvaluator.hpp:1971
DNDS::Euler::EulerEvaluator::OutputOverlapDataRefs::alphaPP
ArrayDOFV< 1 > & alphaPP
PP RHS limiter alpha.
Definition EulerEvaluator.hpp:1975
DNDS::Euler::EulerEvaluator::OutputOverlapDataRefs::uRec
ArrayRECV< nVarsFixed > & uRec
Reconstruction coefficients.
Definition EulerEvaluator.hpp:1973
DNDS::Euler::EulerEvaluator::OutputOverlapDataRefs::betaPP
ArrayDOFV< 1 > & betaPP
PP reconstruction limiter beta.
Definition EulerEvaluator.hpp:1974
DNDS::Euler::EulerEvaluator::OutputOverlapDataRefs::u
ArrayDOFV< nVarsFixed > & u
Cell-centered conservative DOFs.
Definition EulerEvaluator.hpp:1972
Eigen::MatrixFMTSafe
Eigen::Matrix wrapper that hides begin/end from fmt.
Definition EigenUtil.hpp:66
v
Eigen::Matrix< real, 5, 1 > v
Definition test_ArrayDOF.cpp:468
Ax
tVec Ax(NCells)
r
tVec r(NCells)
x
tVec x(NCells)
b
tVec b(NCells)
vg
Eigen::Vector3d vg(0, 0, 0)
velo
Eigen::Vector3d velo
Definition test_GasThermo.cpp:282
res
auto res
Definition test_MeshConnectivity_Interpolate.cpp:54
i
res i
Definition test_MeshConnectivity_Interpolate.cpp:205
mesh
DistributedHex3D mesh
Definition test_MeshConnectivity_Interpolate.cpp:929
p
const tPoint const tPoint const tPoint & p
Definition test_Reconstruction.cpp:484
n
Eigen::Vector3d n(1.0, 0.0, 0.0)