DNDSR 0.1.0.dev1+gcd065ad
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 = 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 };
1640
1641 if (Mode == 0)
1642 {
1643 sourceCaller(0);
1644 }
1645 else if (Mode == 1)
1646 {
1647 sourceCaller(1);
1648 }
1649 else if (Mode == 2)
1650 {
1651 sourceCaller(1);
1652 jacobian += retInc.asDiagonal(); //! TODO: make really block jacobian
1653 }
1654 ret += retInc;
1655 }
1656 else if constexpr (Traits::has2EQ)
1657 {
1658 real pMean, asqrMean, Hmean;
1659 real gamma = settings.idealGasProperty.gamma;
1660 Gas::IdealGasThermal(UMeanXy(I4), UMeanXy(0), (UMeanXy(Seq123) / UMeanXy(0)).squaredNorm(),
1661 gamma, pMean, asqrMean, Hmean);
1662 // ! refvalue:
1663 real muRef = settings.idealGasProperty.muGas;
1664 real T = pMean / ((gamma - 1) / gamma * settings.idealGasProperty.CpGas * UMeanXy(0));
1665
1666 real mufPhy, muf;
1667 mufPhy = muf = muEff(UMeanXy, T);
1668
1669 TU retInc;
1670 retInc.setZero(UMeanXy.size());
1671
1672 TU UMeanXyFixed = UMeanXy;
1673
1674 // if constexpr (Traits::has2EQ)
1675 // for (auto f : mesh->cell2face[iCell])
1676 // if (pBCHandler->GetTypeFromID(mesh->GetFaceZone(f)) == BCWall)
1677 // {
1678 // real d1 = dWallC;
1679 // real rhoOmegaaaWall = mufPhy / sqr(d1) * 800;
1680 // UMeanXyFixed(I4 + 2) = rhoOmegaaaWall;
1681 // }
1682
1683 auto sourceCaller = [&](int mode)
1684 {
1685 if (settings.ransModel == RANSModel::RANS_KOSST)
1686 RANS::GetSource_SST<dim>(UMeanXyFixed, DiffUxy, mufPhy, dWallC, vfv->GetCellMaxLenScale(iCell) * settings.SADESScale, retInc, mode);
1687 else if (settings.ransModel == RANSModel::RANS_KOWilcox)
1688 RANS::GetSource_KOWilcox<dim>(UMeanXyFixed, DiffUxy, mufPhy, dWallC, retInc, mode);
1689 else if (settings.ransModel == RANSModel::RANS_RKE)
1690 RANS::GetSource_RealizableKe<dim>(UMeanXyFixed, DiffUxy, mufPhy, dWallC, retInc, mode);
1691 };
1692
1693 if (Mode == 0)
1694 {
1695 sourceCaller(0);
1696 }
1697 else if (Mode == 1)
1698 {
1699 sourceCaller(1);
1700 }
1701 else if (Mode == 2)
1702 {
1703 sourceCaller(1);
1704 jacobian += retInc.asDiagonal(); //! TODO: make really block jacobian
1705 }
1706 ret += retInc;
1707 }
1708 else
1709 {
1710 DNDS_assert(false);
1711 }
1712 // if (Mode == 1)
1713 // std::cout << ret.transpose() << std::endl;
1714 return ret;
1715 }
1716
1717 DNDS_SWITCH_INTELLISENSE(
1718 template <EulerModel model>
1719 ,
1720 template <>)
1721 /** @brief Dispatcher for boundary ghost-value generation, routing to per-BC-type handlers.
1722 *
1723 * Selects the appropriate boundary value generator based on the BC type (farfield,
1724 * wall, inviscid wall, outflow, inflow, total-condition inflow, symmetry, etc.)
1725 * and returns the ghost (right) state for the boundary face.
1726 *
1727 * @param ULxy Left (interior) state; may be modified by some handlers (e.g., fixUL).
1728 * @param ULMeanXy Left cell mean state.
1729 * @param iCell Owner cell index.
1730 * @param iFace Face index.
1731 * @param iG Quadrature point index (-1 for face center).
1732 * @param uNorm Outward unit normal vector.
1733 * @param normBase Local coordinate frame built from the normal.
1734 * @param pPhysics Physical coordinates of the quadrature point.
1735 * @param t Current simulation time.
1736 * @param btype Boundary zone identifier.
1737 * @param fixUL If true, handler may modify ULxy (e.g., enforce zero SA at wall).
1738 * @param geomMode Geometry mode for normal evaluation.
1739 * @param linMode If 1, return linearized boundary perturbation (for Jacobian).
1740 * @return Ghost (right) state for the boundary face.
1741 */
1742 typename EulerEvaluator<model>::TU EulerEvaluator<model>::generateBoundaryValue(
1743 TU &ULxy, //! warning, possible that UL is also modified
1744 const TU &ULMeanXy,
1745 index iCell, index iFace, int iG,
1746 const TVec &uNorm,
1747 const TMat &normBase,
1748 const Geom::tPoint &pPhysics,
1749 real t,
1750 Geom::t_index btype,
1751 bool fixUL,
1752 int geomMode, int linMode)
1753 {
1754 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
1755 DNDS_assert(iG >= -2);
1756
1757 TU URxy;
1758 URxy.resizeLike(ULxy);
1759 auto bTypeEuler = pBCHandler->GetTypeFromID(btype);
1760
1761 // TODO: for all linMode == 1: implement more precise linearized BC
1762
1763 if (bTypeEuler == EulerBCType::BCSpecial ||
1764 bTypeEuler == EulerBCType::BCFar ||
1765 bTypeEuler == EulerBCType::BCOutP)
1766 {
1767 if (linMode == 1)
1768 {
1769 URxy.setZero();
1770 return URxy;
1771 }
1772 DNDS_assert(ULxy(0) > 0);
1773 if (btype == Geom::BC_ID_DEFAULT_FAR ||
1774 bTypeEuler == EulerBCType::BCFar ||
1775 bTypeEuler == EulerBCType::BCOutP)
1776 {
1777 URxy = generateBV_FarField(ULxy, ULMeanXy, iCell, iFace, iG,
1778 uNorm, normBase, pPhysics, t, btype, fixUL, geomMode);
1779 }
1780 else
1781 {
1782 URxy = generateBV_SpecialFar(ULxy, ULMeanXy, iCell, iFace, iG,
1783 uNorm, normBase, pPhysics, t, btype);
1784 }
1785 }
1786 else if (bTypeEuler == EulerBCType::BCWallInvis ||
1787 bTypeEuler == EulerBCType::BCSym)
1788 {
1789 URxy = generateBV_InviscidWall(ULxy, ULMeanXy, iCell, iFace, iG,
1790 uNorm, normBase, pPhysics, t, btype);
1791 if (linMode == 1)
1792 return URxy;
1793 }
1794 else if (bTypeEuler == EulerBCType::BCWall ||
1795 bTypeEuler == EulerBCType::BCWallIsothermal)
1796 {
1797 URxy = generateBV_ViscousWall(ULxy, ULMeanXy, iCell, iFace, iG,
1798 uNorm, normBase, pPhysics, t, btype, fixUL, linMode);
1799 if (linMode == 1)
1800 return URxy;
1801 }
1802 else if (bTypeEuler == EulerBCType::BCOut)
1803 {
1804 URxy = generateBV_Outflow(ULxy, ULMeanXy, iCell, iFace, iG,
1805 uNorm, normBase, pPhysics, t, btype);
1806 if (linMode == 1)
1807 {
1808 URxy.setZero();
1809 return URxy;
1810 }
1811 }
1812 else if (bTypeEuler == EulerBCType::BCIn)
1813 {
1814 if (linMode == 1)
1815 {
1816 URxy.setZero();
1817 return URxy;
1818 }
1819 URxy = generateBV_Inflow(ULxy, ULMeanXy, iCell, iFace, iG,
1820 uNorm, normBase, pPhysics, t, btype);
1821 }
1822 else if (bTypeEuler == EulerBCType::BCInPsTs)
1823 {
1824 if (linMode == 1)
1825 {
1826 URxy.setZero();
1827 return URxy;
1828 }
1829 URxy = generateBV_TotalConditionInflow(ULxy, ULMeanXy, iCell, iFace, iG,
1830 uNorm, normBase, pPhysics, t, btype);
1831 }
1832 else
1833 {
1834 DNDS_assert(false);
1835 }
1836 return URxy;
1837 }
1838
1839 /**************************************************************************
1840 * Per-BC-type handlers
1841 **************************************************************************/
1842
1843 DNDS_SWITCH_INTELLISENSE(
1844 template <EulerModel model>
1845 ,
1846 template <>)
1847 /** @brief Generate ghost state for farfield and back-pressure outlet boundary conditions.
1848 *
1849 * Uses Riemann-invariant-based switching: supersonic outflow extrapolates the
1850 * interior state; subsonic outflow imposes farfield pressure; subsonic inflow
1851 * imposes farfield values but retains interior pressure; supersonic inflow uses
1852 * the full farfield state. Supports anchor-point and radial-profile pressure corrections,
1853 * CL-driver AoA rotation, and rotating-frame transformations.
1854 *
1855 * @param ULxy Left (interior) state.
1856 * @param ULMeanXy Left cell mean state.
1857 * @param iCell Owner cell index.
1858 * @param iFace Face index.
1859 * @param iG Quadrature point index.
1860 * @param uNorm Outward unit normal.
1861 * @param normBase Local coordinate frame.
1862 * @param pPhysics Physical coordinates.
1863 * @param t Current simulation time.
1864 * @param btype Boundary zone identifier.
1865 * @param fixUL If true, may modify ULxy.
1866 * @param geomMode Geometry mode.
1867 * @return Ghost (right) state for the farfield face.
1868 */
1869 typename EulerEvaluator<model>::TU EulerEvaluator<model>::generateBV_FarField(
1870 TU &ULxy, const TU &ULMeanXy,
1871 index iCell, index iFace, int iG,
1872 const TVec &uNorm, const TMat &normBase,
1873 const Geom::tPoint &pPhysics, real t,
1874 Geom::t_index btype, bool fixUL, int geomMode)
1875 {
1876 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
1877 auto bTypeEuler = pBCHandler->GetTypeFromID(btype);
1878 TU URxy;
1879 URxy.resizeLike(ULxy);
1880
1881 TU far = btype >= Geom::BC_ID_DEFAULT_MAX
1882 ? pBCHandler->GetValueFromID(btype)
1883 : TU(settings.farFieldStaticValue);
1884 if (pCLDriver)
1885 far(Seq123) = pCLDriver->GetAOARotation()(Seq012, Seq012) * far(Seq123);
1886
1887 if (bTypeEuler == EulerBCType::BCFar)
1888 {
1889 if (settings.frameConstRotation.enabled && pBCHandler->GetFlagFromID(btype, "frameOpt") != 0)
1890 far(Seq123) = (Geom::RotateAxis(-settings.frameConstRotation.vOmega() * t) * Geom::ToThreeDim<dim>(far(Seq123)))(Seq012);
1891 }
1892
1893 TU ULxyStatic = ULxy;
1894 if (settings.frameConstRotation.enabled) // to static frame velocity
1895 TransformURotatingFrame(ULxyStatic, pPhysics, 1);
1896
1897 real un = ULxy(Seq123).dot(uNorm) / ULxy(0); // using relative velo for in/out judgement
1898 real gamma = settings.idealGasProperty.gamma;
1899 real asqr, H, p;
1900 Gas::IdealGasThermal(ULxyStatic(I4), ULxyStatic(0), (ULxyStatic(Seq123) / ULxyStatic(0)).squaredNorm(), gamma, p, asqr, H);
1901
1902 DNDS_assert(asqr >= 0);
1903 real a = std::sqrt(asqr);
1904
1905 auto vg = this->GetFaceVGrid(iFace, iG, pPhysics);
1906 real vgN = vg.dot(uNorm);
1907
1908 if (un - vgN - a > 0) // full outflow
1909 {
1910 URxy = ULxyStatic;
1911 }
1912 else if (un - vgN > 0) // subsonic out
1913 {
1914 TU farPrimitive, ULxyPrimitive;
1915 farPrimitive.resizeLike(ULxyStatic);
1916 ULxyPrimitive.resizeLike(URxy);
1917 Gas::IdealGasThermalConservative2Primitive<dim>(far, farPrimitive, gamma);
1918 Gas::IdealGasThermalConservative2Primitive<dim>(ULxyStatic, ULxyPrimitive, gamma);
1919 if (bTypeEuler == EulerBCType::BCOutP && pBCHandler->GetFlagFromID(btype, "anchorOpt") == 1)
1920 {
1921 {
1922 TU anchorPointRel = ULxy;
1923 if (anchorRecorders.count(btype)) // if doesn't have anchor value yet, use UL as anchor
1924 anchorPointRel = anchorRecorders.at(btype).val;
1925 TU anchorPointRelPrimitive = anchorPointRel;
1926 Gas::IdealGasThermalConservative2Primitive<dim>(anchorPointRel, anchorPointRelPrimitive, gamma);
1927 farPrimitive(I4) += std::max(ULxyPrimitive(I4) - anchorPointRelPrimitive(I4), -0.95 * farPrimitive(I4));
1928 }
1929 }
1930 if (bTypeEuler == EulerBCType::BCOutP && pBCHandler->GetFlagFromID(btype, "anchorOpt") == 2)
1931 {
1932 real pInc = 0;
1933 if (profileRecorders.count(btype))
1934 pInc = profileRecorders.at(btype).GetPlain(settings.frameConstRotation.rVec(pPhysics).norm())(I4);
1935 farPrimitive(I4) += std::max(pInc, -0.95 * farPrimitive(I4));
1936 }
1937 ULxyPrimitive(I4) = farPrimitive(I4); // using far pressure
1938 Gas::IdealGasThermalPrimitive2Conservative<dim>(ULxyPrimitive, URxy, gamma);
1939 }
1940 else if (un - vgN + a > 0) // subsonic in
1941 {
1942 TU farPrimitive, ULxyPrimitive;
1943 farPrimitive.resizeLike(ULxyStatic);
1944 ULxyPrimitive.resizeLike(URxy);
1945 Gas::IdealGasThermalConservative2Primitive<dim>(far, farPrimitive, gamma);
1946 Gas::IdealGasThermalConservative2Primitive<dim>(ULxyStatic, ULxyPrimitive, gamma);
1947 farPrimitive(I4) = ULxyPrimitive(I4); // using inner pressure
1948 Gas::IdealGasThermalPrimitive2Conservative<dim>(farPrimitive, URxy, gamma);
1949 }
1950 else // full inflow
1951 {
1952 URxy = far;
1953 }
1954 if (settings.frameConstRotation.enabled) // to rotating frame velocity
1955 TransformURotatingFrame(URxy, pPhysics, -1);
1956 return URxy;
1957 }
1958
1959 DNDS_SWITCH_INTELLISENSE(
1960 template <EulerModel model>
1961 ,
1962 template <>)
1963 /** @brief Generate ghost state for special farfield BCs (DMR, Rayleigh-Taylor, etc.).
1964 *
1965 * Handles built-in special boundary conditions identified by reserved BC IDs:
1966 * double Mach reflection (DMR), Rayleigh-Taylor instability, and user-specified
1967 * special options. Uses Riemann-invariant switching like the standard farfield.
1968 *
1969 * @param ULxy Left (interior) state.
1970 * @param ULMeanXy Left cell mean state.
1971 * @param iCell Owner cell index.
1972 * @param iFace Face index.
1973 * @param iG Quadrature point index.
1974 * @param uNorm Outward unit normal.
1975 * @param normBase Local coordinate frame.
1976 * @param pPhysics Physical coordinates.
1977 * @param t Current simulation time.
1978 * @param btype Boundary zone identifier (special reserved ID).
1979 * @return Ghost (right) state for the special farfield face.
1980 */
1981 typename EulerEvaluator<model>::TU EulerEvaluator<model>::generateBV_SpecialFar(
1982 TU &ULxy, const TU &ULMeanXy,
1983 index iCell, index iFace, int iG,
1984 const TVec &uNorm, const TMat &normBase,
1985 const Geom::tPoint &pPhysics, real t,
1986 Geom::t_index btype)
1987 {
1988 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
1989 TU URxy;
1990 URxy.resizeLike(ULxy);
1991
1992 if (btype == Geom::BC_ID_DEFAULT_SPECIAL_DMR_FAR) // DMR
1993 {
1994 DNDS_assert(dim > 1);
1995 URxy = settings.farFieldStaticValue;
1996 real uShock = 10;
1997 if constexpr (dim == 3)
1998 {
1999 if (((pPhysics(0) - uShock / std::sin(pi / 3) * t - 1. / 6.) -
2000 pPhysics(1) / std::tan(pi / 3)) > 0)
2001 URxy({0, 1, 2, 3, 4}) = Eigen::Vector<real, 5>{1.4, 0, 0, 0, 2.5};
2002 else
2003 URxy({0, 1, 2, 3, 4}) = Eigen::Vector<real, 5>{8, 57.157676649772960, -33, 0, 5.635e2};
2004 }
2005 else
2006 {
2007 if (((pPhysics(0) - uShock / std::sin(pi / 3) * t - 1. / 6.) -
2008 pPhysics(1) / std::tan(pi / 3)) > 0)
2009 URxy({0, 1, 2, 3}) = Eigen::Vector<real, 4>{1.4, 0, 0, 2.5};
2010 else
2011 URxy({0, 1, 2, 3}) = Eigen::Vector<real, 4>{8, 57.157676649772960, -33, 5.635e2};
2012 }
2013 }
2014 else if (btype == Geom::BC_ID_DEFAULT_SPECIAL_RT_FAR) // Rayleigh-Taylor
2015 {
2016 DNDS_assert(dim > 1);
2017 Eigen::VectorXd far = settings.farFieldStaticValue;
2018 real gamma = settings.idealGasProperty.gamma;
2019 real un = ULxy(Seq123).dot(uNorm) / ULxy(0);
2020 real vsqr = (ULxy(Seq123) / ULxy(0)).squaredNorm();
2021 real asqr, H, p;
2022 Gas::IdealGasThermal(ULxy(I4), ULxy(0), vsqr, gamma, p, asqr, H);
2023
2024 DNDS_assert(asqr >= 0);
2025 real a = std::sqrt(asqr);
2026 real v = -0.025 * a * cos(pPhysics(0) * 8 * pi);
2027
2028 if (pPhysics(1) < 0.5)
2029 {
2030 real rho = 2;
2031 real p = 1;
2032 far(0) = rho;
2033 far(1) = 0;
2034 far(2) = rho * v;
2035 far(I4) = 0.5 * rho * sqr(v) + p / (gamma - 1);
2036 }
2037 else
2038 {
2039 real rho = 1;
2040 real p = 2.5;
2041 far(0) = rho;
2042 far(1) = 0;
2043 far(2) = rho * v;
2044 far(I4) = 0.5 * rho * sqr(v) + p / (gamma - 1);
2045 }
2046
2047 if (un - a > 0)
2048 URxy = ULxy;
2049 else if (un > 0)
2050 {
2051 TU farPrimitive, ULxyPrimitive;
2052 farPrimitive.resizeLike(ULxy);
2053 ULxyPrimitive.resizeLike(URxy);
2054 Gas::IdealGasThermalConservative2Primitive<dim>(far, farPrimitive, gamma);
2055 Gas::IdealGasThermalConservative2Primitive<dim>(ULxy, ULxyPrimitive, gamma);
2056 ULxyPrimitive(I4) = farPrimitive(I4);
2057 Gas::IdealGasThermalPrimitive2Conservative<dim>(ULxyPrimitive, URxy, gamma);
2058 }
2059 else if (un + a > 0)
2060 {
2061 TU farPrimitive, ULxyPrimitive;
2062 farPrimitive.resizeLike(ULxy);
2063 ULxyPrimitive.resizeLike(URxy);
2064 Gas::IdealGasThermalConservative2Primitive<dim>(far, farPrimitive, gamma);
2065 Gas::IdealGasThermalConservative2Primitive<dim>(ULxy, ULxyPrimitive, gamma);
2066 farPrimitive(I4) = ULxyPrimitive(I4);
2067 Gas::IdealGasThermalPrimitive2Conservative<dim>(farPrimitive, URxy, gamma);
2068 }
2069 else
2070 URxy = far;
2071 }
2072 else if (btype == Geom::BC_ID_DEFAULT_SPECIAL_IV_FAR) // Isentropic Vortex
2073 {
2074 real chi = 5;
2075 real gamma = settings.idealGasProperty.gamma;
2076 real xc = 5 + t;
2077 real yc = 5 + t;
2078 real r = std::sqrt(sqr(pPhysics(0) - xc) + sqr(pPhysics(1) - yc));
2079 real dT = -(gamma - 1) / (8 * gamma * sqr(pi)) * sqr(chi) * std::exp(1 - sqr(r));
2080 real dux = chi / 2 / pi * std::exp((1 - sqr(r)) / 2) * -(pPhysics(1) - xc);
2081 real duy = chi / 2 / pi * std::exp((1 - sqr(r)) / 2) * +(pPhysics(0) - yc);
2082 real T = dT + 1;
2083 real ux = dux + 1;
2084 real uy = duy + 1;
2085 real S = 1;
2086 real rho = std::pow(T / S, 1 / (gamma - 1));
2087 real pVal = T * rho;
2088 real E = 0.5 * (sqr(ux) + sqr(uy)) * rho + pVal / (gamma - 1);
2089
2090 URxy.setZero();
2091 URxy(0) = rho;
2092 URxy(1) = rho * ux;
2093 URxy(2) = rho * uy;
2094 URxy(dim + 1) = E;
2095 }
2096 else if (btype == Geom::BC_ID_DEFAULT_SPECIAL_2DRiemann_FAR) // 2D Riemann
2097 {
2098 real gamma = settings.idealGasProperty.gamma;
2099 real bdL = 0.0, bdR = 1.0, bdD = 0.0, bdU = 1.0;
2100
2101 real phi1 = -0.663324958071080;
2102 real phi2 = -0.422115882408869;
2103 real location = 0.8;
2104 real p1 = location + phi1 * t;
2105 real p2 = location + phi2 * t;
2106 real rho, u, v, pre;
2107 TU ULxyPrimitive;
2108 ULxyPrimitive.resizeLike(ULxy);
2109
2110 Gas::IdealGasThermalConservative2Primitive<dim>(ULxy, ULxyPrimitive, gamma);
2111 real rhoL = ULxyPrimitive(0);
2112 real uL = ULxyPrimitive(1);
2113 real vL = ULxyPrimitive(2);
2114 real preL = ULxyPrimitive(I4);
2115 TU farPrimitive = ULxyPrimitive;
2116
2117 static const real bTol = 1e-9;
2118 if (std::abs(pPhysics(0) - bdL) < bTol)
2119 { // left, phi2
2120 if (pPhysics(1) <= p2)
2121 rho = 0.137992831541219, u = 1.206045378311055, v = 1.206045378311055, pre = 0.029032258064516;
2122 else
2123 rho = 0.532258064516129, u = 1.206045378311055, v = 0.0, pre = 0.3;
2124 }
2125 else if (std::abs(pPhysics(0) - bdR) < bTol)
2126 { // right, phi1
2127 if (pPhysics(1) <= p1)
2128 rho = rhoL, u = -uL, v = vL, pre = preL;
2129 else
2130 rho = rhoL, u = -uL, v = vL, pre = preL;
2131 }
2132 else if (std::abs(pPhysics(1) - bdU) < bTol)
2133 { // up, phi1
2134 if (pPhysics(0) <= p1)
2135 rho = rhoL, u = uL, v = -vL, pre = preL;
2136 else
2137 rho = rhoL, u = uL, v = -vL, pre = preL;
2138 }
2139 else if (std::abs(pPhysics(1) - bdD) < bTol)
2140 { // down, phi2
2141 if (pPhysics(0) <= p2)
2142 rho = 0.137992831541219, u = 1.206045378311055, v = 1.206045378311055, pre = 0.029032258064516;
2143 else
2144 rho = 0.532258064516129, u = 0.0, v = 1.206045378311055, pre = 0.3;
2145 }
2146 else
2147 {
2148 rho = u = v = pre = std::nan("1");
2149 DNDS_assert(false);
2150 }
2151 farPrimitive(0) = rho;
2152 farPrimitive(1) = u, farPrimitive(2) = v;
2153 farPrimitive(I4) = pre;
2154 Gas::IdealGasThermalPrimitive2Conservative<dim>(farPrimitive, URxy, gamma);
2155 }
2156 else if (pBCHandler->GetFlagFromID(btype, "specialOpt") == 3001) // Noh
2157 {
2158 TU farPrimitive;
2159 Gas::IdealGasThermalConservative2Primitive<dim>(settings.farFieldStaticValue, farPrimitive, settings.idealGasProperty.gamma);
2160 real pInf = farPrimitive(I4);
2161 real r = pPhysics.norm();
2162 TVec velo = -pPhysics(Seq012) / (r + smallReal);
2163 real rho = sqr(1. + t / (r + smallReal));
2164 farPrimitive(0) = rho;
2165 farPrimitive(Seq123) = velo;
2166 farPrimitive(I4) = pInf;
2167 Gas::IdealGasThermalPrimitive2Conservative<dim>(farPrimitive, URxy, settings.idealGasProperty.gamma);
2168 }
2169 else
2170 DNDS_assert_info(false, fmt::format(
2171 "btype [{}] or bTypeEuler [{}] or specialOpt [{}] is not supported",
2172 btype, to_string(pBCHandler->GetTypeFromID(btype)),
2173 pBCHandler->GetFlagFromIDSoft(btype, "specialOpt")));
2174 return URxy;
2175 }
2176
2177 DNDS_SWITCH_INTELLISENSE(
2178 template <EulerModel model>
2179 ,
2180 template <>)
2181 /** @brief Generate ghost state for inviscid wall (slip wall) and symmetry boundary conditions.
2182 *
2183 * Mirrors the velocity component normal to the wall while preserving the
2184 * tangential component, density, and energy. Handles rotating-frame transformations.
2185 *
2186 * @param ULxy Left (interior) state.
2187 * @param ULMeanXy Left cell mean state.
2188 * @param iCell Owner cell index.
2189 * @param iFace Face index.
2190 * @param iG Quadrature point index.
2191 * @param uNorm Outward unit normal.
2192 * @param normBase Local coordinate frame.
2193 * @param pPhysics Physical coordinates.
2194 * @param t Current simulation time.
2195 * @param btype Boundary zone identifier.
2196 * @return Ghost (right) state with reflected normal velocity.
2197 */
2198 typename EulerEvaluator<model>::TU EulerEvaluator<model>::generateBV_InviscidWall(
2199 TU &ULxy, const TU &ULMeanXy,
2200 index iCell, index iFace, int iG,
2201 const TVec &uNorm, const TMat &normBase,
2202 const Geom::tPoint &pPhysics, real t,
2203 Geom::t_index btype)
2204 {
2205 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
2206 TU URxy = ULxy;
2207 if (settings.frameConstRotation.enabled)
2208 this->TransformURotatingFrame_ABS_VELO(URxy, pPhysics, -1);
2209 URxy(Seq123) -= 2 * URxy(Seq123).dot(uNorm) * uNorm; // mirrored!
2210 if (settings.frameConstRotation.enabled)
2211 this->TransformURotatingFrame_ABS_VELO(URxy, pPhysics, 1);
2212 return URxy;
2213 }
2214
2215 DNDS_SWITCH_INTELLISENSE(
2216 template <EulerModel model>
2217 ,
2218 template <>)
2219 /** @brief Generate ghost state for viscous (no-slip) wall boundary conditions.
2220 *
2221 * Reverses the velocity to enforce zero velocity at the wall (ghost mirroring).
2222 * For isothermal walls, adjusts density to match the prescribed wall temperature.
2223 * Handles SA turbulence variable (negated or zeroed at wall) and k-omega/k-epsilon
2224 * wall boundary conditions including the realizable k-epsilon wall epsilon formula.
2225 *
2226 * @param ULxy Left (interior) state (may be modified when fixUL is true for SA).
2227 * @param ULMeanXy Left cell mean state.
2228 * @param iCell Owner cell index.
2229 * @param iFace Face index.
2230 * @param iG Quadrature point index.
2231 * @param uNorm Outward unit normal.
2232 * @param normBase Local coordinate frame.
2233 * @param pPhysics Physical coordinates.
2234 * @param t Current simulation time.
2235 * @param btype Boundary zone identifier.
2236 * @param fixUL If true, may zero SA variables in ULxy at the wall.
2237 * @param linMode If 1, return linearized ghost perturbation for Jacobian.
2238 * @return Ghost (right) state with no-slip wall condition.
2239 */
2240 typename EulerEvaluator<model>::TU EulerEvaluator<model>::generateBV_ViscousWall(
2241 TU &ULxy, const TU &ULMeanXy,
2242 index iCell, index iFace, int iG,
2243 const TVec &uNorm, const TMat &normBase,
2244 const Geom::tPoint &pPhysics, real t,
2245 Geom::t_index btype, bool fixUL, int linMode)
2246 {
2247 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
2248 auto bTypeEuler = pBCHandler->GetTypeFromID(btype);
2249 TU URxy = ULxy;
2250 if (true) // physical mode
2251 {
2252 if (settings.frameConstRotation.enabled && pBCHandler->GetFlagFromID(btype, "frameOpt") == 0)
2253 this->TransformURotatingFrame_ABS_VELO(URxy, pPhysics, -1);
2254 if (settings.frameConstRotation.enabled && pBCHandler->GetFlagFromID(btype, "frameOpt") != 0)
2255 this->TransformURotatingFrame(URxy, pPhysics, 1);
2256 URxy(Seq123) *= -1;
2257 if (settings.frameConstRotation.enabled && pBCHandler->GetFlagFromID(btype, "frameOpt") == 0)
2258 this->TransformURotatingFrame_ABS_VELO(URxy, pPhysics, 1);
2259 if (settings.frameConstRotation.enabled && pBCHandler->GetFlagFromID(btype, "frameOpt") != 0)
2260 this->TransformURotatingFrame(URxy, pPhysics, -1);
2261 }
2262
2263 if (linMode == 1)
2264 return URxy;
2265
2266 if (bTypeEuler == EulerBCType::BCWallIsothermal)
2267 {
2268 real temp = pBCHandler->GetValueFromID(btype)(0);
2269 TU URxyPrim;
2270 URxyPrim.resizeLike(ULxy);
2271 Gas::IdealGasThermalConservative2Primitive<dim>(URxy, URxyPrim, settings.idealGasProperty.gamma);
2272 DNDS_assert_info(URxyPrim(0) > 0 && URxyPrim(I4) > 0 && temp > 0, fmt::format("{}, {}, {}", URxyPrim(0), URxyPrim(I4), temp));
2273 real newDensity = URxyPrim(I4) / temp / settings.idealGasProperty.Rgas;
2274 URxyPrim(0) = newDensity;
2275 Gas::IdealGasThermalPrimitive2Conservative<dim>(URxyPrim, URxy, settings.idealGasProperty.gamma);
2276 }
2277 if constexpr (Traits::hasSA)
2278 {
2279 URxy(I4 + 1) *= -1;
2280#ifdef USE_FIX_ZERO_SA_NUT_AT_WALL
2281 if (fixUL)
2282 ULxy(I4 + 1) = URxy(I4 + 1) = 0; //! modifing UL
2283#endif
2284 }
2285 if constexpr (Traits::has2EQ)
2286 {
2287 URxy({I4 + 1, I4 + 2}) *= -1;
2288#ifdef USE_FIX_ZERO_SA_NUT_AT_WALL
2289 // if (fixUL)
2290 // ULxy({I4 + 1, I4 + 2}).setZero(), URxy({I4 + 1, I4 + 2}).setZero();
2291#endif
2292 if (settings.ransModel == RANSModel::RANS_RKE)
2293 { // BC for RealizableKe
2294 real d1 = dWall[iCell].mean();
2295 real k1 = ULMeanXy(I4 + 1) / ULMeanXy(0);
2296
2297 real pMean, asqrMean, Hmean;
2298 real gamma = settings.idealGasProperty.gamma;
2299 Gas::IdealGasThermal(ULMeanXy(I4), ULMeanXy(0), (ULMeanXy(Seq123) / ULMeanXy(0)).squaredNorm(),
2300 gamma, pMean, asqrMean, Hmean);
2301 real T = pMean / ((gamma - 1) / gamma * settings.idealGasProperty.CpGas * ULMeanXy(0));
2302 real mufPhy1 = muEff(ULMeanXy, T);
2303 real epsWall = 2 * mufPhy1 / ULMeanXy(0) * k1 / sqr(d1);
2304 URxy(I4 + 2) = 2 * epsWall * ULxy(0) - ULxy(I4 + 2);
2305 }
2306 if (settings.ransModel == RANSModel::RANS_KOSST ||
2307 settings.ransModel == RANSModel::RANS_KOWilcox)
2308 { // BC for SST or KOWilcox
2309 real d1 = dWall[iCell].mean();
2310 real pMean, asqrMean, Hmean;
2311 real gamma = settings.idealGasProperty.gamma;
2312 Gas::IdealGasThermal(ULMeanXy(I4), ULMeanXy(0), (ULMeanXy(Seq123) / ULMeanXy(0)).squaredNorm(),
2313 gamma, pMean, asqrMean, Hmean);
2314 real T = pMean / ((gamma - 1) / gamma * settings.idealGasProperty.CpGas * ULMeanXy(0));
2315 real mufPhy1 = muEff(ULMeanXy, T);
2316
2317 real rhoOmegaaaWall = mufPhy1 / sqr(d1) * RANS::kWallOmegaCoeff;
2318 URxy(I4 + 2) = 2 * rhoOmegaaaWall - ULxy(I4 + 2);
2319 }
2320 }
2321 return URxy;
2322 }
2323
2324 DNDS_SWITCH_INTELLISENSE(
2325 template <EulerModel model>
2326 ,
2327 template <>)
2328 /** @brief Generate ghost state for supersonic/extrapolation outflow BC.
2329 *
2330 * Simply returns the interior state as the ghost value, allowing all waves to
2331 * exit the domain without reflection.
2332 *
2333 * @param ULxy Left (interior) state.
2334 * @param ULMeanXy Left cell mean state.
2335 * @param iCell Owner cell index.
2336 * @param iFace Face index.
2337 * @param iG Quadrature point index.
2338 * @param uNorm Outward unit normal.
2339 * @param normBase Local coordinate frame.
2340 * @param pPhysics Physical coordinates.
2341 * @param t Current simulation time.
2342 * @param btype Boundary zone identifier.
2343 * @return Ghost state equal to the interior state (full extrapolation).
2344 */
2345 typename EulerEvaluator<model>::TU EulerEvaluator<model>::generateBV_Outflow(
2346 TU &ULxy, const TU &ULMeanXy,
2347 index iCell, index iFace, int iG,
2348 const TVec &uNorm, const TMat &normBase,
2349 const Geom::tPoint &pPhysics, real t,
2350 Geom::t_index btype)
2351 {
2352 return ULxy;
2353 }
2354
2355 DNDS_SWITCH_INTELLISENSE(
2356 template <EulerModel model>
2357 ,
2358 template <>)
2359 /** @brief Generate ghost state for prescribed inflow boundary condition.
2360 *
2361 * Sets the ghost state to the BC-prescribed conservative values. Applies
2362 * CL-driver AoA rotation and rotating-frame transformation if applicable.
2363 *
2364 * @param ULxy Left (interior) state (unused for pure inflow).
2365 * @param ULMeanXy Left cell mean state.
2366 * @param iCell Owner cell index.
2367 * @param iFace Face index.
2368 * @param iG Quadrature point index.
2369 * @param uNorm Outward unit normal.
2370 * @param normBase Local coordinate frame.
2371 * @param pPhysics Physical coordinates.
2372 * @param t Current simulation time.
2373 * @param btype Boundary zone identifier.
2374 * @return Ghost state set to prescribed inflow values.
2375 */
2376 typename EulerEvaluator<model>::TU EulerEvaluator<model>::generateBV_Inflow(
2377 TU &ULxy, const TU &ULMeanXy,
2378 index iCell, index iFace, int iG,
2379 const TVec &uNorm, const TMat &normBase,
2380 const Geom::tPoint &pPhysics, real t,
2381 Geom::t_index btype)
2382 {
2383 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
2384 TU URxy = pBCHandler->GetValueFromID(btype);
2385 if (pCLDriver)
2386 URxy(Seq123) = pCLDriver->GetAOARotation()(Seq012, Seq012) * URxy(Seq123);
2387 // Note: removed dead code that checked bTypeEuler == BCFar (unreachable in BCIn branch)
2388 if (settings.frameConstRotation.enabled)
2389 TransformURotatingFrame(URxy, pPhysics, -1);
2390 return URxy;
2391 }
2392
2393 DNDS_SWITCH_INTELLISENSE(
2394 template <EulerModel model>
2395 ,
2396 template <>)
2397 /** @brief Generate ghost state for total-condition (stagnation) inflow BC.
2398 *
2399 * Given prescribed stagnation pressure and temperature, computes the static
2400 * conditions from the interior velocity magnitude using isentropic relations.
2401 * The flow direction is taken from the BC-prescribed direction vector.
2402 * Applies CL-driver AoA rotation and rotating-frame transformation.
2403 *
2404 * @param ULxy Left (interior) state (used for velocity magnitude).
2405 * @param ULMeanXy Left cell mean state.
2406 * @param iCell Owner cell index.
2407 * @param iFace Face index.
2408 * @param iG Quadrature point index.
2409 * @param uNorm Outward unit normal.
2410 * @param normBase Local coordinate frame.
2411 * @param pPhysics Physical coordinates.
2412 * @param t Current simulation time.
2413 * @param btype Boundary zone identifier.
2414 * @return Ghost state computed from total conditions and interior velocity.
2415 */
2416 typename EulerEvaluator<model>::TU EulerEvaluator<model>::generateBV_TotalConditionInflow(
2417 TU &ULxy, const TU &ULMeanXy,
2418 index iCell, index iFace, int iG,
2419 const TVec &uNorm, const TMat &normBase,
2420 const Geom::tPoint &pPhysics, real t,
2421 Geom::t_index btype)
2422 {
2423 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
2424 TU URxy;
2425 URxy.resizeLike(ULxy);
2426
2427 TU ULxyStatic = ULxy;
2428 if (settings.frameConstRotation.enabled)
2429 TransformURotatingFrame(ULxyStatic, pPhysics, 1);
2430 TU ULxyPrimitive;
2431 ULxyPrimitive.resizeLike(ULxy);
2432 real gamma = settings.idealGasProperty.gamma;
2433 Gas::IdealGasThermalConservative2Primitive<dim>(ULxyStatic, ULxyPrimitive, gamma);
2434 TVec v = ULxyStatic(Seq123).array() / ULxyStatic(0);
2435 real vSqr = v.squaredNorm();
2436 {
2437 TU farPrimitive = pBCHandler->GetValueFromID(btype); // primitive passive scalar components like Nu
2438
2439 real pStag = pBCHandler->GetValueFromID(btype)(0);
2440 real tStag = pBCHandler->GetValueFromID(btype)(1);
2441 vSqr = std::min(vSqr, tStag * 2 * settings.idealGasProperty.CpGas * 0.95);
2442 real tStatic = tStag - 0.5 * vSqr / (settings.idealGasProperty.CpGas);
2443 real pStatic = pStag * std::pow(tStatic / tStag, gamma / (gamma - 1));
2444 real rStatic = pStatic / (settings.idealGasProperty.Rgas * tStatic);
2445 farPrimitive(0) = rStatic;
2446 farPrimitive(Seq123) = pBCHandler->GetValueFromID(btype)(Seq234).normalized() * std::sqrt(vSqr);
2447 farPrimitive(I4) = pStatic;
2448 Gas::IdealGasThermalPrimitive2Conservative<dim>(farPrimitive, URxy, gamma);
2449 }
2450 if (pCLDriver)
2451 URxy(Seq123) = pCLDriver->GetAOARotation()(Seq012, Seq012) * URxy(Seq123);
2452 if (settings.frameConstRotation.enabled)
2453 TransformURotatingFrame(URxy, pPhysics, -1);
2454 return URxy;
2455 }
2456
2457 template <EulerModel model>
2458 /** @brief Register available output scalar fields in the OutputPicker.
2459 *
2460 * Populates the output picker with named lambda functions that extract per-cell
2461 * scalar quantities (conservative components, reconstruction coefficients,
2462 * limiter values, wall distance, cell volume, turbulent viscosity, etc.)
2463 * for use in post-processing output.
2464 *
2465 * @param op OutputPicker to populate with named field extractors (output).
2466 * @param dataRefs References to the solution arrays (u, uRec, betaPP, alphaPP).
2467 */
2468 void EulerEvaluator<model>::InitializeOutputPicker(OutputPicker &op, OutputOverlapDataRefs dataRefs)
2469 {
2470 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
2471
2472 auto &eval = *this;
2473 auto &u = dataRefs.u;
2474 auto &uRec = dataRefs.uRec;
2475 auto &betaPP = dataRefs.betaPP;
2476 auto &alphaPP = dataRefs.alphaPP;
2477
2478 OutputPicker::tMap outMap;
2479 // outMap["R"] = [&](index iCell)
2480 // { return u[iCell](0); };
2481 outMap["RU"] = [&](index iCell)
2482 { return u[iCell](1); };
2483 outMap["RV"] = [&](index iCell)
2484 { return u[iCell](2); };
2485 outMap["RV"] = [&](index iCell)
2486 { return u[iCell](I4 - 1); };
2487 outMap["RE"] = [&](index iCell)
2488 { return u[iCell](I4); };
2489 outMap["R_REC_1"] = [&](index iCell)
2490 { return uRec[iCell](0, 0); };
2491 outMap["RU_REC_1"] = [&](index iCell)
2492 { return uRec[iCell](1, 0); }; // TODO: to be continued to...
2493
2494 // pps:
2495 outMap["betaPP"] = [&](index iCell)
2496 { return betaPP[iCell](0); };
2497 outMap["alphaPP"] = [&](index iCell)
2498 { return alphaPP[iCell](0); };
2499 outMap["ACond"] = [&](index iCell)
2500 {
2501 auto AI = vfv->GetCellRecMatAInv(iCell);
2502 Eigen::MatrixXd AIInv = AI;
2503 real aCond = HardEigen::EigenLeastSquareInverse(AI, AIInv);
2504 return aCond;
2505 };
2506 outMap["dWall"] = [&](index iCell)
2507 {
2508 return eval.dWall.at(iCell).mean();
2509 };
2510 outMap["minJacobiDetRel"] = [&](index iCell)
2511 {
2512 auto eCell = mesh->GetCellElement(iCell);
2513 auto qCell = vfv->GetCellQuad(iCell);
2514 real minDetJac = veryLargeReal;
2515 for (int iG = 0; iG < qCell.GetNumPoints(); iG++)
2516 minDetJac = std::min(vfv->GetCellJacobiDet(iCell, iG), minDetJac);
2517 return minDetJac * Geom::Elem::ParamSpaceVol(eCell.GetParamSpace()) / vfv->GetCellVol(iCell);
2518 };
2519 outMap["cellVolume"] = [&](index iCell)
2520 {
2521 return vfv->GetCellVol(iCell);
2522 };
2523 outMap["mut"] = [&](index iCell)
2524 {
2525 real mut = 0;
2526 if (model == NS_2EQ || model == NS_2EQ_3D)
2527 {
2528 TU Uxy = u[iCell];
2529 TDiffU GradU;
2530 GradU.resize(Eigen::NoChange, nVars);
2531 GradU.setZero();
2532 if constexpr (gDim == 2)
2533 GradU({0, 1}, EigenAll) =
2534 vfv->GetIntPointDiffBaseValue(iCell, -1, -1, -1, std::array<int, 2>{1, 2}, 3) *
2535 uRec[iCell]; // 2d specific
2536 else
2537 GradU({0, 1, 2}, EigenAll) =
2538 vfv->GetIntPointDiffBaseValue(iCell, -1, -1, -1, std::array<int, 3>{1, 2, 3}, 4) *
2539 uRec[iCell]; // 3d specific
2540 real pMean, asqrMean, Hmean;
2541 real gamma = settings.idealGasProperty.gamma;
2542 auto ULMeanXy = Uxy;
2543 Gas::IdealGasThermal(ULMeanXy(I4), ULMeanXy(0), (ULMeanXy(Seq123) / ULMeanXy(0)).squaredNorm(),
2544 gamma, pMean, asqrMean, Hmean);
2545 // ! refvalue:
2546 real muRef = settings.idealGasProperty.muGas;
2547 real T = pMean / ((gamma - 1) / gamma * settings.idealGasProperty.CpGas * ULMeanXy(0));
2548 real mufPhy = muEff(ULMeanXy, T);
2549 if (settings.ransModel == RANSModel::RANS_KOSST)
2550 mut = RANS::GetMut_SST<dim>(Uxy, GradU, mufPhy, dWall[iCell].mean());
2551 else if (settings.ransModel == RANSModel::RANS_KOWilcox)
2552 mut = RANS::GetMut_KOWilcox<dim>(Uxy, GradU, mufPhy, dWall[iCell].mean());
2553 else if (settings.ransModel == RANSModel::RANS_RKE)
2554 mut = RANS::GetMut_RealizableKe<dim>(Uxy, GradU, mufPhy, dWall[iCell].mean());
2555 }
2556
2557 return mut;
2558 };
2559
2560 op.setMap(outMap);
2561 }
2562}
DNDS_MAKE_SSP
#define DNDS_MAKE_SSP(ssp,...)
Definition Defines.hpp:212
DNDS_RESTRICT
#define DNDS_RESTRICT
Definition Defines.hpp:977
DNDS_assert_info
#define DNDS_assert_info(expr, info)
Debug-only assertion with an extra std::string info message.
Definition Errors.hpp:113
DNDS_assert
#define DNDS_assert(expr)
Debug-only assertion (compiled out when DNDS_NDEBUG is defined). Prints the expression + file/line + ...
Definition Errors.hpp:108
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:320
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:2468
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:77
DNDS::Geom::Elem::GetQuadPatches
auto GetQuadPatches(Quadrature &q)
Definition Quadrature.hpp:185
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:203
DNDS::MPI::Barrier
MPI_int Barrier(MPI_Comm comm)
Wrapper over MPI_Barrier.
Definition MPI.cpp:248
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:679
DNDS::UnInitIndex
DNDS_CONSTANT const index UnInitIndex
Sentinel "not initialised" index value (= INT64_MIN).
Definition Defines.hpp:176
DNDS::MPISerialDo
void MPISerialDo(const MPIInfo &mpi, F f)
Execute f on each rank serially, in rank order.
Definition MPI.hpp:698
DNDS::verySmallReal
DNDS_CONSTANT const real verySmallReal
Catch-all lower bound ("effectively zero").
Definition Defines.hpp:194
DNDS::pi
DNDS_CONSTANT const real pi
π in double precision (matches DNDS_E_PI macro).
Definition Defines.hpp:199
DNDS::sign
constexpr real sign(real a)
Signum function: +1, 0, or -1.
Definition Defines.hpp:532
DNDS::MatrixXR
Eigen::Matrix< real, Eigen::Dynamic, Eigen::Dynamic > MatrixXR
Column-major dynamic Eigen matrix of reals (default layout).
Definition Defines.hpp:205
DNDS::index
int64_t index
Global row / DOF index type (signed 64-bit; handles multi-billion-cell meshes).
Definition Defines.hpp:107
DNDS::sqr
constexpr T sqr(const T &a)
a * a, constexpr. Works for all arithmetic types.
Definition Defines.hpp:517
DNDS::real
double real
Canonical floating-point scalar used throughout DNDSR (double precision).
Definition Defines.hpp:105
DNDS::smallReal
DNDS_CONSTANT const real smallReal
Loose lower bound (for iterative-solver tolerances etc.).
Definition Defines.hpp:196
DNDS::log
std::ostream & log()
Return the current DNDSR log stream (either std::cout or the installed file).
Definition Defines.cpp:49
DNDS::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:190
DNDS::DNDS_MPI_REAL
const MPI_Datatype DNDS_MPI_REAL
MPI datatype matching real (= MPI_REAL8).
Definition MPI.hpp:92
DNDS::ArrayTransformerType::Type
ArrayTransformer< typename TArray::value_type, TArray::rs, TArray::rm, TArray::al > Type
Definition ArrayTransformer.hpp:1480
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:62
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_ODE.cpp:196
n
Eigen::Vector3d n(1.0, 0.0, 0.0)