DNDSR 0.1.0.dev1+gcd065ad
Distributed Numeric Data Structure for CFV
Loading...
Searching...
No Matches
EulerEvaluator_EvaluateRHS.hxx
Go to the documentation of this file.
1/** @file EulerEvaluator_EvaluateRHS.hxx
2 * @brief Template implementation of EulerEvaluator::EvaluateRHS, the main spatial
3 * right-hand side evaluation for the compressible Navier-Stokes / Euler equations.
4 *
5 * Covers inviscid flux accumulation over internal and boundary faces, viscous flux
6 * computation, RANS source terms, mass-force and rotating-frame source terms,
7 * boundary integration recording, and optional direct 2nd-order reconstruction modes.
8 */
9#pragma once
10#include "EulerEvaluator.hpp"
11#include "DNDS/CppUtils/ScopedValueAlternator.hpp"
12#include <fmt/core.h>
13
14namespace DNDS::Euler
15{
16 static const auto model = NS_SA;
17 /**
18 * @details
19 * about RHS:
20 * with topology fixed, RHS is dependent on:
21 * flux:
22 * dofs: u_L u_R, urec_L, urec_R,
23 * geoms: dbv_l, dbv_r, detJacobian_f, uNorm_f,
24 */
25#define IF_NOT_NOREC (1)
26 DNDS_SWITCH_INTELLISENSE(
27 // the real definition
28 template <EulerModel model>
29 ,
30 // the intellisense friendly definition
31 template <>)
32 /** @brief Evaluate the spatial right-hand side (RHS) of the semi-discrete equations.
33 *
34 * This is the core spatial operator. It performs the following steps:
35 * 1. Loop over internal faces: reconstruct L/R states, compute inviscid numerical flux
36 * via the selected Riemann solver, and accumulate face contributions to cell RHS.
37 * 2. Loop over boundary faces: generate ghost boundary values, compute boundary flux.
38 * 3. Loop over cells: compute viscous flux, source terms (RANS, mass force, rotating frame),
39 * and add volume contributions.
40 * 4. Track faces where reduced-order reconstruction is used for robustness.
41 * 5. Optionally use direct 2nd-order reconstruction methods for multigrid.
42 *
43 * @param rhs Cell RHS residual (output, zeroed then accumulated).
44 * @param JSource Source-term Jacobian diagonal block (output for implicit).
45 * @param u Conservative variable DOF array.
46 * @param uRecUnlim Unlimited reconstruction coefficients.
47 * @param uRec Limited reconstruction coefficients.
48 * @param uRecBeta Per-cell reconstruction limiter coefficient.
49 * @param cellRHSAlpha Per-cell RHS scaling factor for positivity preservation.
50 * @param onlyOnHalfAlpha If true, evaluate RHS only on cells with alpha < 1.
51 * @param t Current simulation time.
52 * @param flags Bitfield flags controlling viscosity, integration recording,
53 * direct 2nd-order reconstruction, and other options.
54 */
55 void EulerEvaluator<model>::EvaluateRHS(
56 ArrayDOFV<nVarsFixed> &rhs,
57 JacobianDiagBlock<nVarsFixed> &JSource,
58 ArrayDOFV<nVarsFixed> &u,
59 ArrayRECV<nVarsFixed> &uRecUnlim,
60 ArrayRECV<nVarsFixed> &uRec,
61 ArrayDOFV<1> &uRecBeta,
62 ArrayDOFV<1> &cellRHSAlpha,
63 bool onlyOnHalfAlpha,
64 real t,
65 uint64_t flags)
66 {
67 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
68 using namespace Geom;
69 DNDS_MPI_InsertCheck(u.father->getMPI(), "EvaluateRHS 1");
70 int cnvars = nVars;
71
72 rhs.setConstant(0.0);
73 if (settings.useSourceGradFixGG)
74 for (auto &v : gradUFix)
75 v.setZero();
76 nFaceReducedOrder = 0;
77
78 const bool ignoreVis = flags & RHS_Ignore_Viscosity;
79 const bool dontUpdateIntegration = flags & RHS_Dont_Update_Integration;
80 const bool dontUpdateBndFlux = flags & RHS_Dont_Record_Bud_Flux;
81 const bool direct2ndRec = flags & RHS_Direct_2nd_Rec;
82 const bool direct2ndRec1stConv = flags & RHS_Direct_2nd_Rec_1st_Conv;
83 const bool direct2ndUseLimiter = flags & RHS_Direct_2nd_Rec_use_limiter;
84 const bool direct2ndRec_already_have_uGradBufNoLim = flags & RHS_Direct_2nd_Rec_already_have_uGradBufNoLim;
85 const bool recoverIncFScale = flags & RHS_Recover_IncFScale;
86
87 DNDS_assert(direct2ndRec1stConv ? direct2ndRec : true);
88 auto rsType = settings.rsType;
89
90 // recover the upwind factor if using direct2ndRec
91 std::unique_ptr<ScopedValueAlternator<real>> p_rsIncFScaleAlternator;
92 if (recoverIncFScale)
93 p_rsIncFScaleAlternator = std::make_unique<ScopedValueAlternator<real>>(settings.rsIncFScale, 1.0);
94
95#ifdef USE_MG_O1_LLF_FLUX
96 if (direct2ndRec)
97 rsType = Gas::Roe_M2; // to LLF flux
98#endif
99
100 TU fluxWallSumLocal;
101 fluxWallSumLocal.setZero(cnvars);
102 if (!dontUpdateIntegration)
103 {
104 fluxWallSum.setZero(cnvars);
105 for (Geom::t_index i = Geom::BC_ID_DEFAULT_MAX; i < pBCHandler->size(); i++) // init code, consider adding to ctor
106 {
107 if (pBCHandler->GetFlagFromIDSoft(i, "integrationOpt") == 0)
108 continue;
109 if (!bndIntegrations.count(i))
110 {
111 auto intOpt = pBCHandler->GetFlagFromIDSoft(i, "integrationOpt");
112 bndIntegrations.emplace(i, IntegrationRecorder(mesh->getMPI(), intOpt == 1 ? nVars : nVars + 2));
113 }
114 }
115 for (auto &v : bndIntegrations)
116 v.second.Reset();
117 }
118
119 // direct2ndRec = true;
120 if (direct2ndRec)
121 {
122 typename TVFV::template TFBoundary<nVarsFixed>
123 FBoundary = [&](const TU &UL, const TU &UMean, index iCell, index iFace, int ig,
124 const Geom::tPoint &normOut, const Geom::tPoint &pPhy, const Geom::t_index bType) -> TU
125 {
126 TVec normOutV = normOut(Seq012);
127 Eigen::Matrix<real, dim, dim> normBase = Geom::NormBuildLocalBaseV<dim>(normOutV);
128 bool compressed = false;
129 TU ULfixed = this->CompressRecPart(
130 UMean,
131 UL - UMean,
132 compressed);
133 return this->generateBoundaryValue(ULfixed, UMean, iCell, iFace, ig, normOutV, normBase, pPhy, t, bType, true, 1);
134 };
135 if (!direct2ndRec_already_have_uGradBufNoLim)
136 {
137 vfv->DoReconstruction2ndGrad(uGradBufNoLim, u, FBoundary, settings.direct2ndRecMethod);
138 uGradBufNoLim.trans.startPersistentPull(); // this can be safely put before LimiterUGradCall
139 }
140 // uGradBuf = uGradBufNoLim;
141 if (!direct2ndRec1stConv)
142 this->LimiterUGrad(u, uGradBufNoLim, uGradBuf, direct2ndUseLimiter ? LIMITER_UGRAD_No_Flags : LIMITER_UGRAD_Disable_Shock_Limiter);
143 if (!direct2ndRec1stConv)
144 uGradBuf.trans.startPersistentPull();
145 if (!direct2ndRec_already_have_uGradBufNoLim)
146 uGradBufNoLim.trans.waitPersistentPull(); // todo: utilize the uGradBufNoLim to match the implicit reconstruction's "useViscousLimited": false
147 if (!direct2ndRec1stConv)
148 uGradBuf.trans.waitPersistentPull();
149 else
150 uGradBuf = uGradBufNoLim;
151 }
152
153#if defined(DNDS_DIST_MT_USE_OMP)
154 static std::vector<TU> faceFluxBuf;
155 if (faceFluxBuf.size() < mesh->NumFaceProc())
156 {
157 faceFluxBuf.resize(mesh->NumFaceProc(), TU::Zero(nVars));
158 }
159#endif
160
161 auto faceOp = [&](index iFace) {
162
163 };
164
165 double t0 = MPI_Wtime();
166
167#if defined(DNDS_DIST_MT_USE_OMP)
168# pragma omp declare reduction(TUAdd:TU : omp_out += omp_in) initializer(omp_priv = omp_orig)
169# pragma omp parallel for schedule(runtime) reduction(TUAdd : fluxWallSumLocal)
170#endif
171 for (index iFace = 0; iFace < mesh->NumFaceProc(); iFace++)
172 {
173 faceOp(iFace);
174 auto f2c = mesh->face2cell[iFace];
175 Elem::Quadrature gFace = direct2ndRec ? vfv->GetFaceQuadO1(iFace) : vfv->GetFaceQuad(iFace);
176 Eigen::Matrix<real, nVarsFixed, 1, Eigen::ColMajor> fluxEs(cnvars, 1);
177
178 fluxEs.setZero();
179
180 // auto f2n = mesh->face2node[iFace];
181 // Geom::tSmallCoords coords;
182 // mesh->LoadCoords(f2n, coords);
183
184 Geom::Elem::SummationNoOp noOp;
185 auto faceBndID = mesh->GetFaceZone(iFace);
186 auto faceBCType = pBCHandler->GetTypeFromID(faceBndID);
187 TU_Batch ULxyV, URxyV;
188 ULxyV.resize(u.father->MatRowSize(), gFace.GetNumPoints());
189 URxyV.resizeLike(ULxyV);
190 TDiffU_Batch DiffUxyV, DiffUxyPrimV;
191 DiffUxyV.resize(dim * gFace.GetNumPoints(), u.father->MatRowSize());
192 DiffUxyPrimV.resizeLike(DiffUxyV);
193 TVec_Batch unitNormV, vgXYV;
194 unitNormV.resize(dim, gFace.GetNumPoints()), vgXYV.resizeLike(unitNormV);
195
196 TVec unitNormCent = vfv->GetFaceNorm(iFace, -1)(Seq012);
197 TMat normBaseCent = Geom::NormBuildLocalBaseV<dim>(unitNormCent);
198 TU ULMeanXy = u[f2c[0]];
199 this->UFromCell2Face(ULMeanXy, iFace, f2c[0], 0);
200 TU URMeanXy;
201 if (f2c[1] != UnInitIndex)
202 {
203 URMeanXy = u[f2c[1]];
204 this->UFromCell2Face(URMeanXy, iFace, f2c[1], 1);
205 }
206 else
207 {
208 URMeanXy = generateBoundaryValue(
209 ULMeanXy, ULMeanXy, f2c[0], iFace, -1,
210 unitNormCent,
211 normBaseCent,
212 vfv->GetFaceQuadraturePPhys(iFace, -1),
213 t,
214 mesh->GetFaceZone(iFace), false, 0);
215 }
216 TU_Batch FLFix, FRFix;
217#ifdef USE_FLUX_BALANCE_TERM // todo: decide if the flfix and frfix arguments for fluxface should be completely deleted
218 FLFix.setZero(cnvars, gFace.GetNumPoints()), FRFix.setZero(cnvars, gFace.GetNumPoints()); // todo: finish in faceFlux
219#endif
220 TDiffU faceGradFix;
221 if (settings.useSourceGradFixGG)
222 faceGradFix.setZero(Eigen::NoChange, u.father->MatRowSize());
223
224 gFace.IntegrationSimple(
225 noOp,
226 [&](decltype(noOp) &finc, int iG, real w)
227 {
228 int iGQ = direct2ndRec ? -1 : iG;
229 // finc.resizeLike(fluxEs);
230 int nDiff = vfv->GetFaceAtr(iFace).NDIFF;
231 TVec unitNorm = vfv->GetFaceNorm(iFace, iGQ)(Seq012);
232 TMat normBase = Geom::NormBuildLocalBaseV<dim>(unitNorm);
233 PerformanceTimer::Instance().StartTimer(PerformanceTimer::LimiterB);
234
235 TU ULxy = u[f2c[0]];
236 if (direct2ndRec && !direct2ndRec1stConv)
237 ULxy += uGradBuf[f2c[0]].transpose() * (vfv->GetFaceQuadraturePPhysFromCell(iFace, f2c[0], 0, -1) - vfv->GetCellQuadraturePPhys(f2c[0], -1))(SeqG012);
238 else if (!direct2ndRec1stConv)
239 ULxy += (vfv->GetIntPointDiffBaseValue(f2c[0], iFace, 0, iGQ, std::array<int, 1>{0}, 1) *
240 uRec[f2c[0]])
241 .transpose() *
242 IF_NOT_NOREC;
243 this->UFromCell2Face(ULxy, iFace, f2c[0], 0);
244 TU ULxyUnlim;
245 if (&uRecUnlim != &uRec && !direct2ndRec)
246 {
247 ULxyUnlim = u[f2c[0]];
248 ULxyUnlim += (vfv->GetIntPointDiffBaseValue(f2c[0], iFace, 0, iGQ, std::array<int, 1>{0}, 1) *
249 uRecUnlim[f2c[0]])
250 .transpose() *
251 IF_NOT_NOREC;
252 this->UFromCell2Face(ULxyUnlim, iFace, f2c[0], 0);
253 }
254 else
255 ULxyUnlim = ULxy;
256
257 TU URxy, URxyUnlim;
258#ifndef DNDS_FV_EULEREVALUATOR_IGNORE_VISCOUS_TERM
259 TDiffU GradULxy, GradURxy;
260 GradULxy.resize(Eigen::NoChange, cnvars);
261 GradURxy.resize(Eigen::NoChange, cnvars);
262 GradULxy.setZero(), GradURxy.setZero();
263
264 if (direct2ndRec && !direct2ndRec1stConv)
265 GradULxy(SeqG012, EigenAll) = uGradBuf[f2c[0]];
266 else if (!direct2ndRec1stConv)
267 {
268 if constexpr (gDim == 2)
269 GradULxy({0, 1}, EigenAll) =
270 vfv->GetIntPointDiffBaseValue(f2c[0], iFace, 0, iGQ, std::array<int, 2>{1, 2}, 3) *
271 uRecUnlim[f2c[0]] * IF_NOT_NOREC; // 2d here
272 else
273 GradULxy({0, 1, 2}, EigenAll) =
274 vfv->GetIntPointDiffBaseValue(f2c[0], iFace, 0, iGQ, std::array<int, 3>{1, 2, 3}, 4) *
275 uRecUnlim[f2c[0]] * IF_NOT_NOREC; // 3d here
276 }
277 this->DiffUFromCell2Face(GradULxy, iFace, f2c[0], 0);
278 if (ignoreVis)
279 GradULxy *= 0.;
280
281#endif
282 real minVol = vfv->GetCellVol(f2c[0]);
283 // DNDS_MPI_InsertCheck(u.father->getMPI(), "RHS inner 2");
284 real distBary = veryLargeReal;
285 real distBaryPerp = veryLargeReal;
286
287 if (f2c[1] != UnInitIndex)
288 {
289 URxy = u[f2c[1]];
290 if (direct2ndRec && !direct2ndRec1stConv)
291 URxy += uGradBuf[f2c[1]].transpose() * (vfv->GetFaceQuadraturePPhysFromCell(iFace, f2c[1], 1, -1) - vfv->GetCellQuadraturePPhys(f2c[1], -1))(SeqG012);
292 else if (!direct2ndRec1stConv)
293 URxy += (vfv->GetIntPointDiffBaseValue(f2c[1], iFace, 1, iGQ, std::array<int, 1>{0}, 1) *
294 uRec[f2c[1]])
295 .transpose() *
296 IF_NOT_NOREC;
297 this->UFromCell2Face(URxy, iFace, f2c[1], 1);
298 if (&uRecUnlim != &uRec && !direct2ndRec)
299 {
300 URxyUnlim = u[f2c[1]];
301 URxyUnlim += (vfv->GetIntPointDiffBaseValue(f2c[1], iFace, 1, iGQ, std::array<int, 1>{0}, 1) *
302 uRecUnlim[f2c[1]])
303 .transpose() *
304 IF_NOT_NOREC;
305 this->UFromCell2Face(URxyUnlim, iFace, f2c[1], 1);
306 }
307 else
308 URxyUnlim = URxy;
309
310#ifndef DNDS_FV_EULEREVALUATOR_IGNORE_VISCOUS_TERM
311
312 if (direct2ndRec && !direct2ndRec1stConv)
313 GradURxy(SeqG012, EigenAll) = uGradBuf[f2c[1]];
314 else if (!direct2ndRec1stConv)
315 {
316 if constexpr (gDim == 2)
317 GradURxy({0, 1}, EigenAll) =
318 vfv->GetIntPointDiffBaseValue(f2c[1], iFace, 1, iGQ, std::array<int, 2>{1, 2}, 3) *
319 uRecUnlim[f2c[1]] * IF_NOT_NOREC; // 2d here
320 else
321 GradURxy({0, 1, 2}, EigenAll) =
322 vfv->GetIntPointDiffBaseValue(f2c[1], iFace, 1, iGQ, std::array<int, 3>{1, 2, 3}, 4) *
323 uRecUnlim[f2c[1]] * IF_NOT_NOREC; // 3d here
324 }
325 this->DiffUFromCell2Face(GradURxy, iFace, f2c[1], 1);
326 if (ignoreVis)
327 GradURxy *= 0.;
328#endif
329 minVol = std::min(minVol, vfv->GetCellVol(f2c[1]));
330 distBary = (vfv->GetOtherCellBaryFromCell(f2c[0], f2c[1], iFace) - vfv->GetCellBary(f2c[0])).norm();
331 distBaryPerp =
332 std::abs(
333 (vfv->GetOtherCellBaryFromCell(f2c[0], f2c[1], iFace) -
334 vfv->GetCellBary(f2c[0]))
335 .dot(unitNorm));
336 }
337 else if (true) // is bc
338 {
339 URxy = generateBoundaryValue(
340 ULxy, ULMeanXy, f2c[0], iFace, iGQ,
341 unitNorm,
342 normBase,
343 vfv->GetFaceQuadraturePPhys(iFace, iGQ),
344 t,
345 mesh->GetFaceZone(iFace), true, 0);
346 if (&uRecUnlim != &uRec && false) //! disabled now, as ULxyUnlim may have negative pressure failing in generateBV
347 URxyUnlim = generateBoundaryValue(
348 ULxyUnlim, ULMeanXy, f2c[0], iFace, iGQ,
349 unitNorm,
350 normBase,
351 vfv->GetFaceQuadraturePPhys(iFace, iGQ),
352 t,
353 mesh->GetFaceZone(iFace), true, 0);
354 else
355 URxyUnlim = URxy;
356#ifndef DNDS_FV_EULEREVALUATOR_IGNORE_VISCOUS_TERM
357 GradURxy = GradULxy;
358#endif
359 distBary = (vfv->GetFaceQuadraturePPhysFromCell(iFace, f2c[0], 0, -1) - vfv->GetCellBary(f2c[0])).norm() * 2.;
360 distBaryPerp =
361 std::abs(
362 (vfv->GetFaceQuadraturePPhysFromCell(iFace, f2c[0], 0, -1) -
363 vfv->GetCellBary(f2c[0]))
364 .dot(unitNorm)) *
365 2.;
366 }
367 PerformanceTimer::Instance().StopTimer(PerformanceTimer::LimiterB);
368
369 real distGRP = minVol / vfv->GetFaceArea(iFace) * 2;
370 distGRP = std::max(std::min(distBaryPerp, distGRP * 2), distGRP * 0.25); //! USING REAL GEOMETRICAL
371 if (direct2ndRec1stConv)
372 distGRP = distBary;
373 if (settings.noGRPOnWall && !direct2ndRec1stConv)
374 distGRP += (faceBCType == EulerBCType::BCWall ||
375 faceBCType == EulerBCType::BCWallIsothermal)
376 ? veryLargeReal
377 : 0.0;
378
379 distGRP += faceBCType == EulerBCType::BCWallInvis ? veryLargeReal : 0.0;
380 distGRP += faceBCType == EulerBCType::BCSym ? veryLargeReal : 0.0;
381 TU UMeanXy = 0.5 * (ULxy + URxy);
382
383#ifndef DNDS_FV_EULEREVALUATOR_IGNORE_VISCOUS_TERM
384 TDiffU GradUMeanXy = (GradURxy + GradULxy) * 0.5 +
385 (1.0 / distGRP) *
386 (unitNorm * (URxyUnlim - ULxyUnlim).transpose());
387 if (!GradUMeanXy.allFinite())
388 {
389 std::cout << "GradUMeanXy\n"
390 << GradUMeanXy << "\n";
391 std::cout << "GradULxy\n"
392 << GradULxy << "\n";
393 std::cout << "GradURxy\n"
394 << GradURxy << "\n";
395 std::cout << std::endl;
396
397 DNDS_assert(false);
398 }
399
400 if (ignoreVis)
401 GradUMeanXy *= 0.;
402
403 TDiffU GradUMeanXyPrim;
404 real gamma = settings.idealGasProperty.gamma;
405 if (settings.usePrimGradInVisFlux)
406 {
407 TDiffU GradULxyPrim, GradURxyPrim;
408 GradULxyPrim.resizeLike(GradURxy), GradURxyPrim.resizeLike(GradURxy);
409 Gas::GradientCons2Prim_IdealGas<dim>(ULxy, GradULxy, GradULxyPrim, gamma);
410 Gas::GradientCons2Prim_IdealGas<dim>(URxy, GradURxy, GradURxyPrim, gamma);
411 TU URxyPrim(cnvars), ULxyPrim(cnvars);
412 Gas::IdealGasThermalConservative2Primitive<dim>(ULxy, ULxyPrim, gamma);
413 Gas::IdealGasThermalConservative2Primitive<dim>(URxy, URxyPrim, gamma);
414
415 GradUMeanXyPrim = (GradURxyPrim + GradULxyPrim) * 0.5 +
416 (1.0 / distGRP) *
417 (unitNorm * (URxyPrim - ULxyPrim).transpose());
418 }
419 else
420 Gas::GradientCons2Prim_IdealGas(UMeanXy, GradUMeanXy, GradUMeanXyPrim, gamma);
421
422#else
423 TDiffU GradUMeanXy;
424#endif
425 if (settings.useSourceGradFixGG)
426 {
427 faceGradFix += 0.5 * unitNorm * (URxy - ULxy).transpose() * ((direct2ndRec ? vfv->GetFaceArea(iFace) / vfv->GetFaceParamArea(iFace) : vfv->GetFaceJacobiDet(iFace, iG)) * w);
428 }
429 if (faceBCType == EulerBCType::BCWallInvis ||
430 // faceBCType == EulerBCType::BCIn ||
431 faceBCType == EulerBCType::BCOut ||
432 faceBCType == EulerBCType::BCFar ||
433 faceBCType == EulerBCType::BCSpecial ||
434 faceBCType == EulerBCType::BCSym)
435 GradUMeanXy *= 0, GradUMeanXyPrim *= 0; // force no viscid flux
436
437 if (!GradUMeanXy.allFinite())
438 {
439 std::cout << GradURxy << std::endl;
440 std::cout << GradULxy << std::endl;
441 std::cout << distGRP << std::endl;
442 std::cout << f2c[0] << " " << f2c[1] << " " << mesh->NumCell() << " " << mesh->NumCellProc() << std::endl;
443 std::cout << uRec[f2c[0]] << std::endl;
444 std::cout << "-----------------------------------\n";
445 if (f2c[1] != UnInitIndex)
446 std::cout << uRec[f2c[1]] << std::endl;
447 std::cout << "-----------------------------------\n";
448 std::cout << u[f2c[0]].transpose() << std::endl;
449 if (f2c[1] != UnInitIndex)
450 std::cout << u[f2c[1]].transpose() << std::endl;
451 DNDS_assert(false);
452 }
453
454 if ((faceBCType == EulerBCType::BCWall ||
455 faceBCType == EulerBCType::BCWallIsothermal) &&
456 settings.noRsOnWall)
457 {
458 TU ULc = ULxy;
459 TU ULcPrim;
460 Gas::IdealGasThermalConservative2Primitive<dim>(ULc, ULcPrim, settings.idealGasProperty.gamma);
461 ULcPrim(Seq123).setZero();
462 Gas::IdealGasThermalPrimitive2Conservative<dim>(ULcPrim, ULc, settings.idealGasProperty.gamma);
463 if (faceBCType == EulerBCType::BCWallIsothermal)
464 {
465 real temp = pBCHandler->GetValueFromID(mesh->GetFaceZone(iFace))(0);
466 TU ULcPrim;
467 ULcPrim.resizeLike(ULc);
468 Gas::IdealGasThermalConservative2Primitive<dim>(ULc, ULcPrim, settings.idealGasProperty.gamma);
469 DNDS_assert(ULcPrim(0) > 0 && temp > 0);
470 DNDS_assert_info(ULcPrim(0) > 0 && ULcPrim(I4) > 0 && temp > 0, fmt::format("{}, {}, {}", ULcPrim(0), ULcPrim(I4), temp));
471 // real newPressure = ULcPrim(0) * settings.idealGasProperty.Rgas * temp;
472 // ULcPrim(I4) = newPressure;
473 real newDensity = ULcPrim(I4) / temp / settings.idealGasProperty.Rgas;
474 ULcPrim(0) = newDensity;
475 Gas::IdealGasThermalPrimitive2Conservative<dim>(ULcPrim, URxy, settings.idealGasProperty.gamma);
476 }
477 ULxy = ULc;
478 URxy = ULc;
479 }
480
481 auto seqC = Eigen::seq(iG * dim, iG * dim + dim - 1);
482 ULxyV(EigenAll, iG) = ULxy;
483 URxyV(EigenAll, iG) = URxy;
484 if (!ignoreVis)
485 {
486 DiffUxyV(seqC, EigenAll) = GradUMeanXy;
487 DiffUxyPrimV(seqC, EigenAll) = GradUMeanXyPrim;
488 }
489 unitNormV(EigenAll, iG) = unitNorm;
490 vgXYV(EigenAll, iG) = GetFaceVGrid(iFace, iGQ);
491 });
492 TReal_Batch lam0V, lam123V, lam4V;
493
494 TU_Batch fincC;
495 fluxFace(
496 ULxyV, URxyV,
497 ULMeanXy, URMeanXy,
498 DiffUxyV, DiffUxyPrimV,
499 unitNormV,
500 vgXYV,
501 unitNormCent,
502 GetFaceVGrid(iFace, -1),
503 FLFix,
504 FRFix,
505 fincC,
506 lam0V, lam123V, lam4V,
507 mesh->GetFaceZone(iFace),
508 rsType,
509 iFace, ignoreVis);
510 if (mesh->getMPI().rank == 0)
511 {
512 // std::cout << fincC << std::endl;
513 // std::cout << ULxyV << std::endl;
514 // std::cout << "======" << std::endl;
515 // DNDS_assert(false);
516 }
517
518 gFace.IntegrationSimple(
519 fluxEs,
520 [&](decltype(fluxEs) &finc, int iG)
521 {
522 finc.resizeLike(fluxEs);
523 finc(EigenAll, 0) = fincC(EigenAll, iG);
524 finc *= (direct2ndRec ? vfv->GetFaceArea(iFace) / vfv->GetFaceParamArea(iFace) : vfv->GetFaceJacobiDet(iFace, iG)); // !don't forget this
525 });
526
527 if (settings.useRoeJacobian)
528 {
529 Eigen::Vector<real, 3> lamEs;
530 lamEs.setZero();
531 gFace.IntegrationSimple(
532 lamEs,
533 [&](decltype(lamEs) &finc, int iG)
534 {
535 finc(0) = lam0V(iG);
536 finc(1) = lam123V(iG);
537 finc(2) = lam4V(iG);
538 finc *= (direct2ndRec ? vfv->GetFaceArea(iFace) / vfv->GetFaceParamArea(iFace) : vfv->GetFaceJacobiDet(iFace, iG)); // !don't forget this
539 });
540 lamEs /= vfv->GetFaceArea(iFace);
541 lambdaFace0[iFace] = lamEs(0);
542 lambdaFace123[iFace] = lamEs(1);
543 lambdaFace4[iFace] = lamEs(2);
544 }
545
546#if defined(DNDS_DIST_MT_USE_OMP)
547 faceFluxBuf.at(iFace) = fluxEs(EigenAll, 0);
548#else
549 // ! original code why alphaFace not used?
550 TU fluxIncL = fluxEs(EigenAll, 0);
551 TU fluxIncR = -fluxEs(EigenAll, 0);
552
553 this->UFromFace2Cell(fluxIncL, iFace, f2c[0], 0);
554 if (f2c[1] != UnInitIndex)
555 this->UFromFace2Cell(fluxIncR, iFace, f2c[1], 1); // periodic back to cell
556 // real alphaFace = cellRHSAlpha[f2c[0]](0);
557 // if (f2c[1] != UnInitIndex)
558 // alphaFace = std::min(alphaFace, cellRHSAlpha[f2c[1]](0));
559
560 rhs[f2c[0]] += fluxIncL / vfv->GetCellVol(f2c[0]);
561 if (f2c[1] != UnInitIndex)
562 rhs[f2c[1]] += fluxIncR / vfv->GetCellVol(f2c[1]);
563#endif
564
565 if (settings.useSourceGradFixGG)
566#if defined(DNDS_DIST_MT_USE_OMP)
567// todo: save face value to buffer
568# pragma omp critical
569#endif
570 {
571 TDiffU faceGradFixL{faceGradFix}, faceGradFixR{faceGradFix};
572 if (f2c[0] < mesh->NumCell())
573 this->DiffUFromCell2Face(faceGradFixL, iFace, f2c[0], 0, true), gradUFix[f2c[0]] += faceGradFixL;
574 if (f2c[1] < mesh->NumCell() && f2c[1] != UnInitIndex)
575 this->DiffUFromCell2Face(faceGradFixR, iFace, f2c[1], 1, true), gradUFix[f2c[1]] += faceGradFixR;
576 }
577
578 // record bc flux and tangential bc flux
579 if (!dontUpdateBndFlux)
580 if (f2c[1] == UnInitIndex)
581 {
582 DNDS_assert(mesh->face2bndM.find(iFace) != mesh->face2bndM.end());
583 fluxBnd.at(mesh->face2bndM[iFace]) = fluxEs(EigenAll, 0) / vfv->GetFaceArea(iFace);
584 TVec fluxBndForceTInt;
585 fluxBndForceTInt.setZero();
586 gFace.IntegrationSimple(
587 fluxBndForceTInt,
588 [&](decltype(fluxBndForceTInt) &finc, int iG)
589 {
590 TU fcur = fincC(EigenAll, iG);
591 TVec ncur = unitNormV(EigenAll, iG);
592 finc = fcur(Seq123);
593 finc -= ncur * (ncur.dot(finc));
594 finc *= (direct2ndRec ? vfv->GetFaceArea(iFace) / vfv->GetFaceParamArea(iFace) : vfv->GetFaceJacobiDet(iFace, iG)); // !don't forget this
595 });
596 fluxBndForceT.at(mesh->face2bndM[iFace]) = fluxBndForceTInt / vfv->GetFaceArea(iFace);
597 }
598
599 // integrate BCWall flux
600 if (!dontUpdateIntegration)
601 if (faceBCType == EulerBCType::BCWall || // TODO: update to general
602 faceBCType == EulerBCType::BCWallIsothermal ||
603 (faceBCType == EulerBCType::BCWallInvis && settings.idealGasProperty.muGas < 1e-99))
604 {
605 fluxWallSumLocal -= fluxEs(EigenAll, 0);
606 if (iFace >= mesh->NumFace())
607 DNDS_assert(false);
608 }
609
610 // integrations
611 if (!dontUpdateIntegration)
612#if defined(DNDS_DIST_MT_USE_OMP) // todo: use reduction
613# pragma omp critical
614#endif
615 {
616 if (pBCHandler->GetFlagFromIDSoft(mesh->GetFaceZone(iFace), "integrationOpt") == 1)
617 {
618 bndIntegrations.at(mesh->GetFaceZone(iFace)).Add(-fluxEs(EigenAll, 0), vfv->GetFaceArea(iFace));
619 }
620 if (pBCHandler->GetFlagFromIDSoft(mesh->GetFaceZone(iFace), "integrationOpt") == 2)
621 {
622 TU uL = u[f2c[0]];
623 if (settings.frameConstRotation.enabled)
624 this->TransformURotatingFrame(uL, vfv->GetFaceQuadraturePPhys(iFace, -1), 1);
625 TU uLPrim = uL;
626 auto gamma = settings.idealGasProperty.gamma;
627 Gas::IdealGasThermalConservative2Primitive<dim>(uL, uLPrim, gamma);
628 Eigen::Vector<real, Eigen::Dynamic> vInt;
629 vInt.resize(nVars + 2);
630 vInt(Eigen::seq(0, nVars - 1)) = uL;
631
632 auto [p0, T0] = Gas::IdealGasThermalPrimitiveGetP0T0<dim>(uLPrim, gamma, settings.idealGasProperty.Rgas);
633 vInt(nVars) = p0, vInt(nVars + 1) = T0;
634 vInt(0) = 1;
635 bndIntegrations.at(mesh->GetFaceZone(iFace)).Add(vInt * fluxEs(0, 0), fluxEs(0, 0));
636 }
637 }
638 }
639
640#if defined(DNDS_DIST_MT_USE_OMP)
641# pragma omp parallel for schedule(static)
642 for (int iPart = 0; iPart < mesh->NLocalParts(); iPart++)
643 for (index iCell = mesh->LocalPartStart(iPart); iCell < mesh->LocalPartEnd(iPart); iCell++)
644 {
645 auto c2f = mesh->cell2face[iCell];
646 for (int ic2f = 0; ic2f < c2f.size(); ic2f++)
647 {
648 index iFace = c2f[ic2f];
649 int if2c = mesh->CellIsFaceBack(iCell, iFace) ? 0 : 1;
650 TU fluxFaceC = faceFluxBuf[iFace] * (if2c ? -1 : 1);
651 this->UFromFace2Cell(fluxFaceC, iFace, iCell, if2c);
652
653 rhs[iCell] += fluxFaceC / vfv->GetCellVol(iCell);
654 if (mesh->face2cell(iFace, 1 - if2c) == iCell) // check for self-facing face
655 {
656 TU fluxFaceC = faceFluxBuf[iFace] * (if2c ? -1 : 1) * (-1);
657 this->UFromFace2Cell(fluxFaceC, iFace, iCell, 1 - if2c); // use 1-if2c to force use the other periodic state
658 rhs[iCell] += fluxFaceC / vfv->GetCellVol(iCell);
659 }
660 }
661 }
662#endif
663 double t1 = MPI_Wtime();
664
665 DNDS_MPI_InsertCheck(u.father->getMPI(), "EvaluateRHS After Flux");
666
667 auto cellOp = [&](index iCell) {
668
669 };
670
671 if (!settings.ignoreSourceTerm)
672 {
673 JSource.clearValues();
674#if defined(DNDS_DIST_MT_USE_OMP)
675# pragma omp parallel for schedule(guided)
676#endif
677 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
678 {
679 cellOp(iCell);
680 auto gCell = direct2ndRec ? vfv->GetCellQuadO1(iCell) : vfv->GetCellQuad(iCell);
681
682 TDiffU cellGrad2nd;
683 if (settings.ransSource2nd || settings.source2nd)
684 {
685 cellGrad2nd.setZero(Eigen::NoChange, u[iCell].size());
686 TU uC = u[iCell];
687 for (index iFace : mesh->cell2face[iCell])
688 {
689 index iCellOther = mesh->CellFaceOther(iCell, iFace);
690 TVec uNorm = vfv->GetFaceNormFromCell(iFace, iCell, -1, -1)(Seq012) *
691 (mesh->CellIsFaceBack(iCell, iFace) ? 1 : -1);
692 TU uR;
693 if (iCellOther != UnInitIndex)
694 uR = u[iCellOther], this->UFromOtherCell(uR, iFace, iCell, iCellOther, -1);
695 else
696 uR = generateBoundaryValue(
697 uC, uC, iCell, iFace, -1,
698 uNorm, Geom::NormBuildLocalBaseV<dim>(uNorm),
699 vfv->GetFaceQuadraturePPhys(iFace, -1),
700 t, mesh->GetFaceZone(iFace), true, 0);
701 cellGrad2nd += vfv->GetFaceArea(iFace) * uNorm * (uR - uC).transpose() * 0.5;
702 }
703 cellGrad2nd /= vfv->GetCellVol(iCell);
704 }
705
706 Eigen::Matrix<real, nVarsFixed, Eigen::Dynamic> sourceV; // now includes sourcejacobian diag
707 sourceV.setZero(cnvars, JSource.isBlock() ? cnvars + 1 : 2);
708
709 Geom::Elem::SummationNoOp noOp;
710
711 TDiffU cellGradFix;
712 if (settings.useSourceGradFixGG)
713 cellGradFix = gradUFix[iCell] / vfv->GetCellVol(iCell);
714
715 gCell.IntegrationSimple(
716 sourceV,
717 [&](decltype(sourceV) &finc, int iG)
718 {
719 int iGQ = direct2ndRec ? -1 : iG;
720 TDiffU GradU;
721 GradU.resize(Eigen::NoChange, cnvars);
722 GradU.setZero();
723 PerformanceTimer::Instance().StartTimer(PerformanceTimer::LimiterB);
724 if (direct2ndRec) // should use limited version here or not?
725 GradU(SeqG012, EigenAll) = uGradBufNoLim[iCell];
726 else if (settings.source2nd)
727 GradU = cellGrad2nd;
728 else
729 {
730 if constexpr (gDim == 2)
731 GradU({0, 1}, EigenAll) =
732 vfv->GetIntPointDiffBaseValue(iCell, -1, -1, iGQ, std::array<int, 2>{1, 2}, 3) *
733 uRecUnlim[iCell] * IF_NOT_NOREC; // 2d specific
734 else
735 GradU({0, 1, 2}, EigenAll) =
736 vfv->GetIntPointDiffBaseValue(iCell, -1, -1, iGQ, std::array<int, 3>{1, 2, 3}, 4) *
737 uRecUnlim[iCell] * IF_NOT_NOREC; // 3d specific
738 if (settings.useSourceGradFixGG)
739 GradU += cellGradFix;
740 if (settings.ransSource2nd)
741 {
742 if constexpr (Traits::hasSA)
743 GradU(EigenAll, I4 + 1) = cellGrad2nd(EigenAll, I4 + 1);
744 if constexpr (Traits::has2EQ)
745 GradU(EigenAll, {I4 + 1, I4 + 2}) = cellGrad2nd(EigenAll, {I4 + 1, I4 + 2});
746 }
747 }
748
749 TU ULxy = u[iCell];
750 if (!direct2ndRec)
751 ULxy +=
752 (vfv->GetIntPointDiffBaseValue(iCell, -1, -1, iGQ, std::array<int, 1>{0}, 1) *
753 uRec[iCell])
754 .transpose() *
755 IF_NOT_NOREC;
756
757 PerformanceTimer::Instance().StopTimer(PerformanceTimer::LimiterB);
758
759 // bool compressed = false;
760 // ULxy = CompressRecPart(u[iCell], ULxy, compressed); //! do not forget the mean value
761
762 finc.resizeLike(sourceV);
763 TJacobianU jac;
764 finc(EigenAll, 0) =
765 source(
766 ULxy,
767 GradU,
768 vfv->GetCellQuadraturePPhys(iCell, iGQ), jac,
769 iCell, iGQ, 0);
770 TU sourceJDiag =
771 source(
772 ULxy,
773 GradU,
774 vfv->GetCellQuadraturePPhys(iCell, iGQ), jac,
775 iCell, iGQ, JSource.isBlock() ? 2 : 1);
776 if (JSource.isBlock())
777 finc(EigenAll, Eigen::seq(Eigen::fix<1>, EigenLast)) = jac;
778 else
779 finc(EigenAll, 1) = sourceJDiag;
780
781 finc *= direct2ndRec ? vfv->GetCellVol(iCell) / vfv->GetCellParamVol(iCell) : vfv->GetCellJacobiDet(iCell, iG); //! don't forget this
782 if (finc.hasNaN() || (!finc.allFinite()))
783 {
784 std::cout << finc.transpose() << std::endl;
785 std::cout << ULxy.transpose() << std::endl;
786 std::cout << GradU << std::endl;
787 DNDS_assert(false);
788 }
789 });
790 sourceV *= cellRHSAlpha[iCell](0) / vfv->GetCellVol(iCell); // becomes mean value
791 rhs[iCell] += sourceV(EigenAll, 0);
792 if (JSource.isBlock())
793 JSource.getBlock(iCell) = sourceV(EigenAll, Eigen::seq(Eigen::fix<1>, EigenLast));
794 else
795 JSource.getDiag(iCell) = sourceV(EigenAll, 1);
796
797 // if (iCell == 18195)
798 // {
799 // std::cout << rhs[iCell].transpose() << std::endl;
800 // }
801 }
802 }
803
804 // quick aux: reduce the wall flux sum
805 if (!dontUpdateIntegration)
806 MPI::Allreduce(fluxWallSumLocal.data(), fluxWallSum.data(), fluxWallSum.size(), DNDS_MPI_REAL, MPI_SUM, u.father->getMPI().comm);
807 if (!dontUpdateIntegration)
808 for (auto &i : bndIntegrations)
809 i.second.Reduce();
810 DNDS_MPI_InsertCheck(u.father->getMPI(), "EvaluateRHS -1");
811
812 double t2 = MPI_Wtime();
813
814 // if (u.father->getMPI().rank == 0)
815 // {
816 // std::cout << fmt::format("ti01 [{}] ti12 [{}]", t1 - t0, t2 - t1) << std::endl;
817 // }
818 }
819}
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.
IF_NOT_NOREC
#define IF_NOT_NOREC
Definition EulerEvaluator_EvaluateRHS.hxx:25
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
ScopedValueAlternator.hpp
RAII helper that temporarily overrides a variable and restores its original value when the scope exit...
DNDS::Euler::EulerEvaluator::EvaluateRHS
void EvaluateRHS(ArrayDOFV< nVarsFixed > &rhs, JacobianDiagBlock< nVarsFixed > &JSource, ArrayDOFV< nVarsFixed > &u, ArrayRECV< nVarsFixed > &uRecUnlim, ArrayRECV< nVarsFixed > &uRec, ArrayDOFV< 1 > &uRecBeta, ArrayDOFV< 1 > &cellRHSAlpha, bool onlyOnHalfAlpha, real t, uint64_t flags=RHS_No_Flags)
Evaluate the spatial right-hand side of the semi-discrete system.
DNDS::PerformanceTimer::LimiterB
@ LimiterB
Limiter sub-phase B.
Definition Profiling.hpp:35
DNDS::PerformanceTimer::Instance
static PerformanceTimer & Instance()
Access the process-wide singleton.
Definition Profiling.cpp:9
DNDS::PerformanceTimer::StartTimer
void StartTimer(TimerType t)
Record the current wall time in the "start" slot for timer t.
Definition Profiling.cpp:15
DNDS::PerformanceTimer::StopTimer
void StopTimer(TimerType t)
Add (now - start) to the accumulated time for timer t.
Definition Profiling.cpp:25
DNDS::EulerP::TU
Eigen::Vector< real, nVarsFlow > TU
Fixed-size 5-component conservative state vector (rho, rhoU, rhoV, rhoW, E).
Definition EulerP.hpp:33
DNDS::EulerP::TDiffU
Eigen::Matrix< real, 3, nVarsFlow > TDiffU
Fixed-size 3x5 spatial gradient of the conservative state (one row per spatial dimension).
Definition EulerP.hpp:34
DNDS::Euler::Gas::GradientCons2Prim_IdealGas
void GradientCons2Prim_IdealGas(const TU &U, const TGradU &GradU, TGradUPrim &GradUPrim, real gamma)
Converts the gradient of conservative variables to the gradient of primitive variables for an ideal g...
Definition Gas.hpp:1678
DNDS::Euler::NS_SA
@ NS_SA
NS + Spalart-Allmaras, 2D geometry (6 vars).
Definition Euler.hpp:877
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::BCSym
@ BCSym
Symmetry plane boundary.
Definition EulerBC.hpp:46
DNDS::Euler::BCFar
@ BCFar
Far-field (characteristic-based) boundary.
Definition EulerBC.hpp:38
DNDS::Geom::t_index
int32_t t_index
Definition Geometric.hpp:6
DNDS::Geom::tPoint
Eigen::Vector3d tPoint
Definition Geometric.hpp:9
DNDS::UnInitIndex
DNDS_CONSTANT const index UnInitIndex
Sentinel "not initialised" index value (= INT64_MIN).
Definition Defines.hpp:176
DNDS::index
int64_t index
Global row / DOF index type (signed 64-bit; handles multi-billion-cell meshes).
Definition Defines.hpp:107
DNDS::real
double real
Canonical floating-point scalar used throughout DNDSR (double precision).
Definition Defines.hpp:105
DNDS::veryLargeReal
DNDS_CONSTANT const real veryLargeReal
Catch-all upper bound ("practically infinity") for physical scalars.
Definition Defines.hpp:190
v
Eigen::Matrix< real, 5, 1 > v
Definition test_ArrayDOF.cpp:468