DNDSR 0.1.0.dev1+gcd065ad
Distributed Numeric Data Structure for CFV
Loading...
Searching...
No Matches
EulerEvaluator.hpp
Go to the documentation of this file.
1/**
2 * @file EulerEvaluator.hpp
3 * @brief Core finite-volume evaluator for compressible Navier-Stokes / Euler equations.
4 *
5 * Provides the EulerEvaluator class template, parameterized by EulerModel, which
6 * encapsulates the spatial discretization of the compressible Navier-Stokes equations
7 * using Compact Finite Volume (CFV) methods with variational reconstruction.
8 *
9 * Responsibilities include:
10 * - Right-hand side (RHS) evaluation of the semi-discrete system
11 * - Local time step estimation
12 * - Implicit Jacobian assembly (LU-SGS / SGS / GMRES preconditioning)
13 * - Boundary condition ghost-value generation for all supported BC types
14 * - Wall distance computation (AABB tree, batched AABB, p-Poisson)
15 * - Positivity-preserving limiters and increment compression
16 * - RANS turbulence model source terms (SA, k-omega SST, k-omega Wilcox, Realizable k-e)
17 * - Periodic boundary handling and rotating-frame transformations
18 *
19 * Supported model specializations (via EulerModel enum):
20 * NS, NS_2D, NS_SA, NS_SA_3D, NS_2EQ, NS_2EQ_3D, NS_3D
21 *
22 * @see EulerSolver.hpp Top-level solver orchestrating time marching
23 * @see EulerBC.hpp Boundary condition handler
24 * @see Gas.hpp Gas thermodynamic and Riemann solver routines
25 */
26#pragma once
27
28// #ifndef __DNDS_REALLY_COMPILING__
29// #define __DNDS_REALLY_COMPILING__
30// #define __DNDS_REALLY_COMPILING__HEADER_ON__
31// #endif
32
33#include "Gas.hpp"
34#include "Geom/Mesh.hpp"
35#include "CFV/VariationalReconstruction.hpp"
36#include "DNDS/JsonUtil.hpp"
37#include "Euler.hpp"
38#include "EulerBC.hpp"
39#include "EulerJacobian.hpp"
40#include "EulerEvaluatorSettings.hpp"
41#include "DNDS/SerializerBase.hpp"
42#include "RANS_ke.hpp"
43
44// #ifdef __DNDS_REALLY_COMPILING__HEADER_ON__
45// #undef __DNDS_REALLY_COMPILING__
46// #endif
47
48#include "DNDS/JsonUtil.hpp"
49#include "fmt/core.h"
50#include <iomanip>
51#include <functional>
52
53// #define DNDS_FV_EULEREVALUATOR_SOURCE_TERM_ZERO
54// // #define DNDS_FV_EULEREVALUATOR_IGNORE_SOURCE_TERM
55// // #define DNDS_FV_EULEREVALUATOR_IGNORE_VISCOUS_TERM
56
57// // #ifdef DNDS_FV_EULEREVALUATOR_IGNORE_SOURCE_TERM // term dependency
58// // // #define DNDS_FV_EULEREVALUATOR_USE_SCALAR_JACOBIAN
59// // #endif
60
61namespace DNDS::Euler
62{
63
64 /**
65 * @brief Core finite-volume evaluator for compressible Navier-Stokes / Euler equations.
66 *
67 * Implements the spatial discretization of the compressible N-S equations on
68 * unstructured meshes using Compact Finite Volume (CFV) with variational
69 * reconstruction. Handles inviscid and viscous flux evaluation, source terms
70 * (including RANS turbulence models), implicit Jacobian assembly, boundary
71 * condition ghost-value generation, and positivity-preserving fixes.
72 *
73 * @tparam model EulerModel enum selecting the equation set and spatial dimension
74 * (NS, NS_2D, NS_SA, NS_SA_3D, NS_2EQ, NS_2EQ_3D, NS_3D).
75 */
76 template <EulerModel model = NS>
77 class EulerEvaluator
78 {
79 public:
80 using Traits = EulerModelTraits<model>; ///< Compile-time traits: hasSA, has2EQ, etc.
81 static const int nVarsFixed = getnVarsFixed(model); ///< Number of conserved variables (compile-time fixed).
82 static const int dim = getDim_Fixed(model); ///< Spatial dimension (2 or 3).
83 static const int gDim = getGeomDim_Fixed(model); ///< Geometric dimension of the mesh.
84 static const auto I4 = dim + 1; ///< Index of the energy equation (= dim+1).
85
86 static const int MaxBatch = 16; ///< Maximum batch size for vectorized quadrature-point evaluation.
87 /// @brief Compute the maximum batch-multiplied column count for batched Eigen matrices.
88 static constexpr int MaxBatchMult(int n) { return MaxBatch > 0 ? (n * MaxBatch) : Eigen::Dynamic; }
89
90 typedef Eigen::VectorFMTSafe<real, dim> TVec; ///< Spatial vector (dim components).
91 typedef Eigen::MatrixFMTSafe<real, dim, Eigen::Dynamic, Eigen::ColMajor, dim, MaxBatch> TVec_Batch; ///< Batch of spatial vectors.
92 typedef Eigen::MatrixFMTSafe<real, dim, dim> TMat; ///< Spatial matrix (dim x dim).
93 typedef Eigen::MatrixFMTSafe<real, dim, Eigen::Dynamic, Eigen::ColMajor, dim, MaxBatchMult(3)> TMat_Batch; ///< Batch of spatial matrices.
94 typedef Eigen::VectorFMTSafe<real, nVarsFixed> TU; ///< Conservative variable vector (nVarsFixed components).
95 typedef Eigen::MatrixFMTSafe<real, nVarsFixed, Eigen::Dynamic, Eigen::ColMajor, nVarsFixed, MaxBatch> TU_Batch; ///< Batch of conservative variable vectors.
96 typedef Eigen::MatrixFMTSafe<real, 1, Eigen::Dynamic, Eigen::RowMajor, 1, MaxBatch> TReal_Batch; ///< Batch of scalar values.
97 typedef Eigen::MatrixFMTSafe<real, nVarsFixed, nVarsFixed> TJacobianU; ///< Jacobian matrix (nVars x nVars) of the flux w.r.t. conserved variables.
98 typedef Eigen::MatrixFMTSafe<real, dim, nVarsFixed> TDiffU; ///< Gradient of conserved variables (dim x nVars).
99 typedef Eigen::MatrixFMTSafe<real, Eigen::Dynamic, nVarsFixed, Eigen::ColMajor, MaxBatchMult(3)> TDiffU_Batch; ///< Batch of gradient matrices.
100 typedef Eigen::MatrixFMTSafe<real, nVarsFixed, dim> TDiffUTransposed; ///< Transposed gradient (nVars x dim).
101 typedef ArrayDOFV<nVarsFixed> TDof; ///< Cell-centered DOF array (mean values).
102 typedef ArrayRECV<nVarsFixed> TRec; ///< Reconstruction coefficient array (per-cell polynomial coefficients).
103 typedef ArrayRECV<1> TScalar; ///< Scalar reconstruction coefficient array.
104
105 typedef CFV::VariationalReconstruction<gDim> TVFV; ///< Variational reconstruction type for this geometric dimension.
106 typedef ssp<CFV::VariationalReconstruction<gDim>> TpVFV; ///< Shared pointer to the variational reconstruction object.
107 typedef ssp<BoundaryHandler<model>> TpBCHandler; ///< Shared pointer to the boundary condition handler.
108
109 public:
110 /**
111 * @brief Initialize the finite-volume infrastructure on the mesh.
112 *
113 * Sets periodic transformations on the VFV object, constructs cell/face
114 * geometric metrics, builds reconstruction base functions and weights
115 * (with BC-type-dependent Dirichlet/Neumann weighting), and computes
116 * reconstruction coefficients.
117 *
118 * @param mesh Shared pointer to the unstructured mesh.
119 * @param vfv Shared pointer to the variational reconstruction object.
120 * @param pBCHandler Shared pointer to the boundary condition handler (used
121 * for per-face reconstruction weight selection).
122 */
123 static void InitializeFV(ssp<Geom::UnstructuredMesh> &mesh, TpVFV vfv, TpBCHandler pBCHandler)
124 {
125 vfv->SetPeriodicTransformations(
126 [mesh](auto u, Geom::t_index id) //! important caveat! using & captures mesh shared_ptr as reference, lost if inbound pointer is destoried!!
127 {
128 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
129 u(EigenAll, Seq123) = mesh->periodicInfo.TransVector<dim, Eigen::Dynamic>(u(EigenAll, Seq123).transpose(), id).transpose();
130 },
131 [mesh](auto u, Geom::t_index id)
132 {
133 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
134 u(EigenAll, Seq123) = mesh->periodicInfo.TransVectorBack<dim, Eigen::Dynamic>(u(EigenAll, Seq123).transpose(), id).transpose();
135 });
136
137 vfv->ConstructMetrics();
138 vfv->ConstructBaseAndWeight(
139 [&](Geom::t_index id, int iOrder) -> real
140 {
141 auto type = pBCHandler->GetTypeFromID(id);
142 if (type == BCSpecial || type == BCOut)
143 return 0;
144 if (type == BCFar) // use Dirichlet type
145 return iOrder ? 0. : 1.;
146 if (type == BCWallInvis || type == BCSym)
147 return iOrder ? 0. : 1.;
148 if (Geom::FaceIDIsPeriodic(id))
149 return iOrder ? 1. : 1.; //! treat as real internal
150 // others: use Dirichlet type
151 return iOrder ? 0. : 1.;
152 });
153 vfv->ConstructRecCoeff();
154 }
155
156 public:
157 // static const int gdim = 2; //* geometry dim
158
159 private:
160 int nVars = 5; ///< Runtime number of conserved variables (may differ from nVarsFixed for dynamic models).
161
162 bool passiveDiscardSource = false; ///< When true, discard source terms for passive scalar equations.
163
164 public:
165 /// @brief Enable or disable discarding of passive-scalar source terms.
166 void setPassiveDiscardSource(bool n) { passiveDiscardSource = n; }
167
168 private:
169 int axisSymmetric = 0; ///< Non-zero if the solver operates in axisymmetric mode (2D axisymmetric).
170
171 public:
172 /// @brief Return the axisymmetric mode flag.
173 int GetAxisSymmetric() { return axisSymmetric; }
174
175 private:
176 public:
177 ssp<Geom::UnstructuredMesh> mesh; ///< Shared pointer to the unstructured mesh.
178 ssp<CFV::VariationalReconstruction<gDim>> vfv; //! gDim -> 3 for intellisense //!tmptmp ///< Variational reconstruction object.
179 ssp<BoundaryHandler<model>> pBCHandler; ///< Boundary condition handler.
180 int kAv = 0; ///< Artificial viscosity polynomial order (maxOrder + 1).
181
182 // buffer for fdtau
183 // std::vector<real> lambdaCell;
184 std::vector<real> lambdaFace; ///< Per-face spectral radius (inviscid + viscous combined).
185 std::vector<real> lambdaFaceC; ///< Per-face convective spectral radius.
186 std::vector<real> lambdaFaceVis; ///< Per-face viscous spectral radius.
187 std::vector<real> lambdaFace0; ///< Per-face eigenvalue |u·n| (contact wave).
188 std::vector<real> lambdaFace123; ///< Per-face eigenvalue |u·n + a| (acoustic wave).
189 std::vector<real> lambdaFace4; ///< Per-face eigenvalue |u·n - a| (acoustic wave).
190 std::vector<real> deltaLambdaFace; ///< Per-face spectral radius difference for implicit diagonal.
191 ArrayDOFV<1> deltaLambdaCell; ///< Per-cell accumulated spectral radius difference.
192
193 // grad fix
194 std::vector<TDiffU> gradUFix; ///< Green-Gauss gradient correction buffer for source term stabilization.
195
196 // wall distance
197 std::vector<Eigen::Vector<real, Eigen::Dynamic>> dWall; ///< Per-cell wall distance (one value per cell node/quadrature point).
198 std::vector<real> dWallFace; ///< Per-face wall distance (interpolated from cell values).
199
200 // maps from bc id to various objects
201 std::map<Geom::t_index, AnchorPointRecorder<nVarsFixed>> anchorRecorders; ///< Per-BC anchor point recorders for profile extraction.
202 std::map<Geom::t_index, OneDimProfile<nVarsFixed>> profileRecorders; ///< Per-BC 1-D profile recorders.
203 std::map<Geom::t_index, IntegrationRecorder> bndIntegrations; ///< Per-BC boundary flux/force integration accumulators.
204 std::map<Geom::t_index, std::ofstream> bndIntegrationLogs; ///< Per-BC log file streams for integration output.
205
206 std::set<Geom::t_index> cLDriverBndIDs; ///< Boundary IDs driven by the CL (lift-coefficient) driver.
207 std::unique_ptr<CLDriver> pCLDriver; ///< Lift-coefficient driver for AoA adaptation (null if unused).
208
209 // ArrayVDOF<25> dRdUrec;
210 // ArrayVDOF<25> dRdb;
211 ArrayGRADV<nVarsFixed, gDim> uGradBuf, uGradBufNoLim; ///< Gradient buffers: limited and unlimited.
212
213 Eigen::Vector<real, -1> fluxWallSum; ///< Accumulated wall flux integral (force on wall boundaries).
214 std::vector<TU> fluxBnd; ///< Per-boundary-face flux values.
215 std::vector<TVec> fluxBndForceT; ///< Per-boundary-face tangential force.
216 index nFaceReducedOrder = 0; ///< Count of faces where reconstruction order was reduced.
217
218 ssp<Direct::SerialSymLUStructure> symLU; ///< Symmetric LU structure for direct preconditioner.
219
220 EulerEvaluatorSettings<model> settings; ///< Physics and numerics settings for this evaluator.
221
222 /**
223 * @brief Construct an EulerEvaluator and initialize all internal buffers.
224 *
225 * Allocates per-face spectral-radius buffers, per-boundary flux arrays,
226 * gradient correction arrays (if enabled), wall distance fields (for RANS),
227 * and the symmetric LU structure for direct preconditioning. Also validates
228 * CL-driver boundary configuration.
229 *
230 * @param Nmesh Shared pointer to the unstructured mesh.
231 * @param Nvfv Shared pointer to the variational reconstruction object.
232 * @param npBCHandler Shared pointer to the boundary condition handler.
233 * @param nSettings Physics and numerics settings.
234 * @param n_nVars Runtime number of conserved variables (default from model).
235 */
236 EulerEvaluator(const decltype(mesh) &Nmesh, const decltype(vfv) &Nvfv, const decltype(pBCHandler) &npBCHandler,
237 const EulerEvaluatorSettings<model> &nSettings,
238 int n_nVars = getNVars(model))
239 : nVars(n_nVars), axisSymmetric(Nvfv->GetAxisSymmetric()), mesh(Nmesh), vfv(Nvfv), pBCHandler(npBCHandler), kAv(Nvfv->getSettings().maxOrder + 1), settings(nSettings)
240 {
241 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
242 if (getNVars(model) == DynamicSize)
243 DNDS_assert_info(nVars >= getDim_Fixed(model) + 2, "nVars too small");
244 else
245 DNDS_assert_info(nVars == getNVars(model), "do not change the nVars for this model");
246
247 vfv->BuildUGrad(uGradBuf, nVars);
248 vfv->BuildUGrad(uGradBufNoLim, nVars);
249
250 if (axisSymmetric)
251 DNDS_assert_info(!settings.ignoreSourceTerm, "you have set source term, do not use ignoreSourceTerm! ");
252
253 // lambdaCell.resize(mesh->NumCellProc()); // but only dist part are used, ghost part to not judge for it in facial iter
254 lambdaFace.resize(mesh->NumFaceProc());
255 lambdaFaceC.resize(mesh->NumFaceProc());
256 lambdaFaceVis.resize(lambdaFace.size());
257 lambdaFace0.resize(mesh->NumFaceProc(), 0.);
258 lambdaFace123.resize(mesh->NumFaceProc(), 0.);
259 lambdaFace4.resize(mesh->NumFaceProc(), 0.);
260
261 deltaLambdaFace.resize(lambdaFace.size());
262 vfv->BuildUDof(deltaLambdaCell, 1);
263
264 if (settings.useSourceGradFixGG)
265 {
266 gradUFix.resize(mesh->NumCell());
267 for (auto &v : gradUFix)
268 v.resize(Eigen::NoChange, nVars);
269 }
270
271 fluxBnd.resize(mesh->NumBnd());
272 for (auto &v : fluxBnd)
273 v.resize(nVars);
274 fluxBndForceT.resize(mesh->NumBnd());
275
276 this->GetWallDist();
277
278 if (model == NS_2EQ || model == NS_2EQ_3D)
279 {
280 TU farPrim = settings.farFieldStaticValue;
281 real gamma = settings.idealGasProperty.gamma;
282 Gas::IdealGasThermalConservative2Primitive<dim>(settings.farFieldStaticValue, farPrim, gamma);
283 real T = farPrim(I4) / ((gamma - 1) / gamma * settings.idealGasProperty.CpGas * farPrim(0));
284 // auto [rhs0, rhs] = RANS::SolveZeroGradEquilibrium<dim>(settings.farFieldStaticValue, this->muEff(settings.farFieldStaticValue, T));
285 // if(mesh->getMPI().rank == 0)
286 // log()
287 // << "EulerEvaluator===EulerEvaluator: got 2EQ init for farFieldStaticValue: " << settings.farFieldStaticValue.transpose() << "\n"
288 // << fmt::format(" [{:.3e} -> {:.3e}] ", rhs0, rhs) << std::endl;
289 }
290
291 symLU = std::make_shared<Direct::SerialSymLUStructure>(mesh->getMPI(), mesh->NumCell());
292
293 for (auto &name : settings.cLDriverBCNames)
294 {
295 auto bcID = pBCHandler->GetIDFromName(name);
296 cLDriverBndIDs.insert(bcID);
297 if (bcID >= Geom::BC_ID_DEFAULT_MAX)
298 DNDS_assert_info(pBCHandler->GetFlagFromIDSoft(bcID, "integrationOpt") == 1,
299 "you have to set integrationOption == 1 to make this bc available by CLDriver");
300 else
301 DNDS_assert_info(bcID == Geom::BC_ID_DEFAULT_WALL || bcID == Geom::BC_ID_DEFAULT_WALL_INVIS,
302 "default bc must be WALL or WALL_INVIS for CLDriver");
303 }
304 if (cLDriverBndIDs.size())
305 pCLDriver = std::make_unique<CLDriver>(settings.cLDriverSettings);
306 }
307
308 /**
309 * @brief Compute wall distance for all cells (dispatcher).
310 *
311 * Selects the wall distance algorithm based on the wallDistScheme setting
312 * (AABB tree, batched AABB, p-Poisson, or combination) and populates dWall
313 * and dWallFace arrays used by RANS turbulence models.
314 */
315 void GetWallDist();
316
317 private:
318 /// Collect wall-boundary triangles for CGAL AABB queries.
319 /// @param useQuadPatches If true, triangulate using quadrature patches (scheme 1);
320 /// otherwise use element vertices directly (scheme 0/20).
321 void GetWallDist_CollectTriangles(bool useQuadPatches,
322 std::vector<Eigen::Matrix<real, 3, 3>> &trianglesOut);
323
324 /// Wall distance via global CGAL AABB tree (wallDistScheme 0, 1, 20 first pass).
325 void GetWallDist_AABB();
326
327 /// Wall distance via batched per-rank CGAL AABB queries (wallDistScheme 2, 20 refine pass).
328 void GetWallDist_BatchedAABB();
329
330 /// Wall distance via p-Poisson equation (wallDistScheme 3).
331 void GetWallDist_Poisson();
332
333 /// Compute dWallFace from dWall (shared postprocess).
334 void GetWallDist_ComputeFaceDistances();
335
336 public:
337
338 /******************************************************/
339 static const uint64_t DT_No_Flags = 0x0ull; ///< No flags for EvaluateDt.
340 static const uint64_t DT_Dont_update_lambda01234 = 0x1ull << 0; ///< Skip recomputation of per-face eigenvalues lambda0/123/4.
341 /**
342 * @brief Estimate the local or global time step for each cell.
343 *
344 * Computes the local pseudo-time step dTau based on CFL number, spectral radii
345 * of inviscid and viscous fluxes, and optionally enforces a global minimum.
346 * Updates per-face eigenvalue arrays (lambda0, lambda123, lambda4) unless
347 * DT_Dont_update_lambda01234 is set.
348 *
349 * @param[out] dt Per-cell time step array (overwritten).
350 * @param[in] u Cell-centered conservative variable DOFs.
351 * @param[in] uRec Reconstruction coefficients.
352 * @param[in] CFL CFL number.
353 * @param[out] dtMinall Global minimum time step (MPI-reduced).
354 * @param[in] MaxDt Upper bound on the time step.
355 * @param[in] UseLocaldt If true, use local (per-cell) time stepping.
356 * @param[in] t Current simulation time.
357 * @param[in] flags Bitwise combination of DT_* flags.
358 */
359 void EvaluateDt(
360 ArrayDOFV<1> &dt,
361 ArrayDOFV<nVarsFixed> &u,
362 ArrayRECV<nVarsFixed> &uRec,
363 real CFL, real &dtMinall, real MaxDt,
364 bool UseLocaldt,
365 real t,
366 uint64_t flags = DT_No_Flags);
367
368 /// @name RHS evaluation flags (bitwise OR combinable)
369 /// @{
370 static const uint64_t RHS_No_Flags = 0x0ull; ///< Default: full RHS evaluation.
371 static const uint64_t RHS_Ignore_Viscosity = 0x1ull << 0; ///< Skip viscous flux contribution.
372 static const uint64_t RHS_Dont_Update_Integration = 0x1ull << 1; ///< Skip boundary integration accumulation.
373 static const uint64_t RHS_Dont_Record_Bud_Flux = 0x1ull << 2; ///< Skip recording per-boundary flux.
374 static const uint64_t RHS_Direct_2nd_Rec = 0x1ull << 8; ///< Use direct 2nd-order reconstruction.
375 static const uint64_t RHS_Direct_2nd_Rec_1st_Conv = 0x1ull << 9; ///< 2nd-order rec with 1st-order convection.
376 static const uint64_t RHS_Direct_2nd_Rec_use_limiter = 0x1ull << 10; ///< Apply limiter when using direct 2nd rec.
377 static const uint64_t RHS_Direct_2nd_Rec_already_have_uGradBufNoLim = 0x1ull << 11; ///< uGradBufNoLim is already computed.
378 static const uint64_t RHS_Recover_IncFScale = 0x1ull << 12; ///< Recover incremental face scaling.
379 /// @}
380
381 /**
382 * @brief Evaluate the spatial right-hand side of the semi-discrete system.
383 *
384 * Computes inviscid flux (Riemann solver), viscous flux, and source terms
385 * over all internal and boundary faces, accumulating cell residuals. Also
386 * records boundary force/flux integrations and updates reduced-order face counts.
387 *
388 * @param[out] rhs Cell residual array (overwritten).
389 * @param[out] JSource Diagonal Jacobian block from source terms.
390 * @param[in] u Cell-centered conservative DOFs.
391 * @param[in] uRecUnlim Unlimited reconstruction coefficients.
392 * @param[in] uRec Limited reconstruction coefficients.
393 * @param[in] uRecBeta Per-cell reconstruction compression factor (PP limiter).
394 * @param[in] cellRHSAlpha Per-cell RHS scaling factor (PP limiter).
395 * @param[in] onlyOnHalfAlpha If true, evaluate only cells with alpha < 1.
396 * @param[in] t Current simulation time.
397 * @param[in] flags Bitwise combination of RHS_* flags.
398 */
399 void EvaluateRHS(
400 ArrayDOFV<nVarsFixed> &rhs,
401 JacobianDiagBlock<nVarsFixed> &JSource,
402 ArrayDOFV<nVarsFixed> &u,
403 ArrayRECV<nVarsFixed> &uRecUnlim,
404 ArrayRECV<nVarsFixed> &uRec,
405 ArrayDOFV<1> &uRecBeta,
406 ArrayDOFV<1> &cellRHSAlpha,
407 bool onlyOnHalfAlpha,
408 real t,
409 uint64_t flags = RHS_No_Flags);
410
411 /**
412 * @brief Assemble the diagonal blocks of the implicit Jacobian for LU-SGS / SGS.
413 *
414 * Computes J_diag = (V/dTau + alphaDiag * sum_faces(spectral_radius)) * I + J_source
415 * for each cell, where V is cell volume and dTau is the local pseudo-time step.
416 *
417 * @param[out] JDiag Per-cell diagonal Jacobian block (overwritten).
418 * @param[in] JSource Source-term Jacobian contribution to diagonal.
419 * @param[in] dTau Per-cell local pseudo-time step.
420 * @param[in] dt Physical time step (for dual time stepping).
421 * @param[in] alphaDiag Diagonal scaling factor for implicit relaxation.
422 * @param[in] u Cell-centered conservative DOFs.
423 * @param[in] uRec Reconstruction coefficients.
424 * @param[in] jacobianCode Controls Jacobian approximation: 0=scalar, 1=analytical flux Jacobian.
425 * @param[in] t Current simulation time.
426 */
427 void LUSGSMatrixInit(
428 JacobianDiagBlock<nVarsFixed> &JDiag,
429 JacobianDiagBlock<nVarsFixed> &JSource,
430 ArrayDOFV<1> &dTau, real dt, real alphaDiag,
431 ArrayDOFV<nVarsFixed> &u,
432 ArrayRECV<nVarsFixed> &uRec,
433 int jacobianCode,
434 real t);
435
436 /**
437 * @brief Compute the matrix-vector product A * uInc for the implicit system.
438 *
439 * Evaluates the action of the approximate Jacobian on an increment vector,
440 * used as the matvec operation inside GMRES.
441 *
442 * @param[in] alphaDiag Diagonal scaling factor.
443 * @param[in] t Current simulation time.
444 * @param[in] u Cell-centered conservative DOFs.
445 * @param[in] uInc Increment vector to multiply.
446 * @param[in] JDiag Pre-assembled diagonal Jacobian blocks.
447 * @param[out] AuInc Result of A * uInc (overwritten).
448 */
449 void LUSGSMatrixVec(
450 real alphaDiag,
451 real t,
452 ArrayDOFV<nVarsFixed> &u,
453 ArrayDOFV<nVarsFixed> &uInc,
454 JacobianDiagBlock<nVarsFixed> &JDiag,
455 ArrayDOFV<nVarsFixed> &AuInc);
456
457 /**
458 * @brief Build the local LU factorization of the Jacobian for direct solve.
459 *
460 * Assembles the full local Jacobian (including off-diagonal face coupling)
461 * and factorizes it, storing the result in jacLU for subsequent direct solves.
462 *
463 * @param[in] alphaDiag Diagonal scaling factor.
464 * @param[in] t Current simulation time.
465 * @param[in] u Cell-centered conservative DOFs.
466 * @param[in] JDiag Pre-assembled diagonal Jacobian blocks.
467 * @param[out] jacLU Local LU factorization result (overwritten).
468 */
469 void LUSGSMatrixToJacobianLU(
470 real alphaDiag,
471 real t,
472 ArrayDOFV<nVarsFixed> &u,
473 JacobianDiagBlock<nVarsFixed> &JDiag,
474 JacobianLocalLU<nVarsFixed> &jacLU);
475
476 /**
477 * @brief Deprecated: use UpdateSGS with uIncIsZero=true instead.
478 * Kept for backward compatibility. Delegates to UpdateSGS.
479 */
480 [[deprecated("Use UpdateSGS with uIncIsZero=true instead")]]
481 void UpdateLUSGSForward(
482 real alphaDiag,
483 real t,
484 ArrayDOFV<nVarsFixed> &rhs,
485 ArrayDOFV<nVarsFixed> &u,
486 ArrayDOFV<nVarsFixed> &uInc,
487 JacobianDiagBlock<nVarsFixed> &JDiag,
488 ArrayDOFV<nVarsFixed> &uIncNew);
489
490 /**
491 * @brief Deprecated: use UpdateSGS with uIncIsZero=true instead.
492 * Kept for backward compatibility. Delegates to UpdateSGS.
493 */
494 [[deprecated("Use UpdateSGS with uIncIsZero=true instead")]]
495 void UpdateLUSGSBackward(
496 real alphaDiag,
497 real t,
498 ArrayDOFV<nVarsFixed> &rhs,
499 ArrayDOFV<nVarsFixed> &u,
500 ArrayDOFV<nVarsFixed> &uInc,
501 JacobianDiagBlock<nVarsFixed> &JDiag,
502 ArrayDOFV<nVarsFixed> &uIncNew);
503
504 /**
505 * @brief Symmetric Gauss-Seidel update for the implicit linear system.
506 *
507 * @param forward Scan ascending (true) or descending (false).
508 * @param gsUpdate Use Gauss-Seidel update (read from uIncNew for already-processed cells).
509 * @param uIncIsZero If true, uInc is assumed zero for not-yet-processed cells,
510 * so their flux contribution is skipped (LUSGS optimisation).
511 * @param sumInc Output: accumulated absolute increment for convergence tracking.
512 */
513 void UpdateSGS(
514 real alphaDiag,
515 real t,
516 ArrayDOFV<nVarsFixed> &rhs,
517 ArrayDOFV<nVarsFixed> &u,
518 ArrayDOFV<nVarsFixed> &uInc,
519 ArrayDOFV<nVarsFixed> &uIncNew,
520 JacobianDiagBlock<nVarsFixed> &JDiag,
521 bool forward, bool gsUpdate, TU &sumInc,
522 bool uIncIsZero = false);
523
524 /**
525 * @brief Solve the implicit linear system using the pre-factored local LU.
526 *
527 * @param[in] alphaDiag Diagonal scaling factor.
528 * @param[in] t Current simulation time.
529 * @param[in] rhs Right-hand side residual.
530 * @param[in] u Cell-centered conservative DOFs.
531 * @param[in] uInc Current increment (input guess).
532 * @param[out] uIncNew Updated increment (output).
533 * @param[out] bBuf Buffer for intermediate right-hand side assembly.
534 * @param[in] JDiag Diagonal Jacobian blocks.
535 * @param[in] jacLU Pre-factored local LU decomposition.
536 * @param[in] uIncIsZero If true, skip contributions from zero-increment cells.
537 * @param[out] sumInc Accumulated absolute increment for convergence monitoring.
538 */
539 void LUSGSMatrixSolveJacobianLU(
540 real alphaDiag,
541 real t,
542 ArrayDOFV<nVarsFixed> &rhs,
543 ArrayDOFV<nVarsFixed> &u,
544 ArrayDOFV<nVarsFixed> &uInc,
545 ArrayDOFV<nVarsFixed> &uIncNew,
546 ArrayDOFV<nVarsFixed> &bBuf,
547 JacobianDiagBlock<nVarsFixed> &JDiag,
548 JacobianLocalLU<nVarsFixed> &jacLU,
549 bool uIncIsZero,
550 TU &sumInc);
551
552 /**
553 * @brief SGS sweep coupled with reconstruction update.
554 *
555 * Performs a single SGS sweep that simultaneously updates the conservative
556 * increment (uInc) and the reconstruction increment (uRecInc).
557 *
558 * @param[in] alphaDiag Diagonal scaling factor.
559 * @param[in] t Current simulation time.
560 * @param[in] rhs Right-hand side residual.
561 * @param[in] u Cell-centered conservative DOFs.
562 * @param[in] uRec Reconstruction coefficients.
563 * @param[in] uInc Current increment for conservative variables.
564 * @param[in] uRecInc Current reconstruction increment.
565 * @param[in] JDiag Diagonal Jacobian blocks.
566 * @param[in] forward Sweep direction: ascending (true) or descending (false).
567 * @param[out] sumInc Accumulated absolute increment for convergence monitoring.
568 */
569 void UpdateSGSWithRec(
570 real alphaDiag,
571 real t,
572 ArrayDOFV<nVarsFixed> &rhs,
573 ArrayDOFV<nVarsFixed> &u,
574 ArrayRECV<nVarsFixed> &uRec,
575 ArrayDOFV<nVarsFixed> &uInc,
576 ArrayRECV<nVarsFixed> &uRecInc,
577 JacobianDiagBlock<nVarsFixed> &JDiag,
578 bool forward, TU &sumInc);
579
580 // void UpdateLUSGSForwardWithRec(
581 // real alphaDiag,
582 // ArrayDOFV<nVarsFixed> &rhs,
583 // ArrayDOFV<nVarsFixed> &u,
584 // ArrayRECV<nVarsFixed> &uRec,
585 // ArrayDOFV<nVarsFixed> &uInc,
586 // ArrayRECV<nVarsFixed> &uRecInc,
587 // ArrayDOFV<nVarsFixed> &JDiag,
588 // ArrayDOFV<nVarsFixed> &uIncNew);
589
590 /// @brief Clip extreme conserved-variable values to prevent overflow.
591 void FixUMaxFilter(ArrayDOFV<nVarsFixed> &u);
592
593 /// @brief Accumulate time-averaged primitive variables for unsteady statistics.
594 void TimeAverageAddition(ArrayDOFV<nVarsFixed> &w, ArrayDOFV<nVarsFixed> &wAveraged, real dt, real &tCur);
595
596 /// @brief Convert cell-mean conservative variables to primitive variables.
597 void MeanValueCons2Prim(ArrayDOFV<nVarsFixed> &u, ArrayDOFV<nVarsFixed> &w);
598 /// @brief Convert cell-mean primitive variables to conservative variables.
599 void MeanValuePrim2Cons(ArrayDOFV<nVarsFixed> &w, ArrayDOFV<nVarsFixed> &u);
600
601 using tFCompareField = std::function<TU(const Geom::tPoint &, real)>; ///< Callback type for analytical comparison field.
602 using tFCompareFieldWeight = std::function<real(const Geom::tPoint &, real)>; ///< Callback type for comparison weighting function.
603
604 /**
605 * @brief Compute the norm of the RHS residual vector.
606 *
607 * @param[out] res Norm result per variable (resized to nVars).
608 * @param[in] rhs Cell residual array.
609 * @param[in] P Norm order (1 = L1, 2 = L2, etc.).
610 * @param[in] volWise If true, weight by cell volume.
611 * @param[in] average If true, divide by total volume/count.
612 */
613 void EvaluateNorm(Eigen::Vector<real, -1> &res, ArrayDOFV<nVarsFixed> &rhs, index P = 1, bool volWise = false, bool average = false);
614
615 /**
616 * @brief Compute the reconstruction error norm (optionally against an analytical field).
617 *
618 * @param[out] res Norm result per variable.
619 * @param[in] u Cell-centered conservative DOFs.
620 * @param[in] uRec Reconstruction coefficients.
621 * @param[in] P Norm order.
622 * @param[in] compare If true, compute error against FCompareField.
623 * @param[in] FCompareField Analytical field callback.
624 * @param[in] FCompareFieldWeight Weighting callback.
625 * @param[in] t Current simulation time.
626 */
627 void EvaluateRecNorm(
628 Eigen::Vector<real, -1> &res,
629 ArrayDOFV<nVarsFixed> &u,
630 ArrayRECV<nVarsFixed> &uRec,
631 index P = 1,
632 bool compare = false,
633 const tFCompareField &FCompareField = [](const Geom::tPoint &p, real t)
634 { return TU::Zero(); },
635 const tFCompareFieldWeight &FCompareFieldWeight = [](const Geom::tPoint &p, real t)
636 { return 1.0; },
637 real t = 0);
638
639 /// @name Limiter flags
640 /// @{
641 static const uint64_t LIMITER_UGRAD_No_Flags = 0x0ull; ///< Default limiter flags.
642 static const uint64_t LIMITER_UGRAD_Disable_Shock_Limiter = 0x1ull << 0; ///< Disable shock-detecting component of the limiter.
643 /// @}
644
645 /**
646 * @brief Apply slope limiter to the gradient field.
647 *
648 * Limits the reconstructed gradient (uGrad) to produce a monotonicity-preserving
649 * gradient (uGradNew). Supports WBAP and CWBAP limiter variants.
650 *
651 * @param[in] u Cell-centered conservative DOFs.
652 * @param[in] uGrad Input (unlimited) gradient array.
653 * @param[out] uGradNew Output limited gradient array.
654 * @param[in] flags Bitwise combination of LIMITER_UGRAD_* flags.
655 */
656 void LimiterUGrad(ArrayDOFV<nVarsFixed> &u, ArrayGRADV<nVarsFixed, gDim> &uGrad, ArrayGRADV<nVarsFixed, gDim> &uGradNew,
657 uint64_t flags = LIMITER_UGRAD_No_Flags);
658
659 static const int EvaluateURecBeta_DEFAULT = 0x00; ///< Default: evaluate beta without compression.
660 static const int EvaluateURecBeta_COMPRESS_TO_MEAN = 0x01; ///< Compress reconstruction toward cell mean to enforce positivity.
661 /**
662 * @brief Evaluate the positivity-preserving reconstruction limiter (beta).
663 *
664 * For each cell, computes the maximum compression factor beta such that
665 * u_mean + beta * uRec remains physically realizable (positive density and pressure).
666 *
667 * @param[in] u Cell-centered conservative DOFs.
668 * @param[in] uRec Reconstruction coefficients.
669 * @param[out] uRecBeta Per-cell compression factor in [0,1].
670 * @param[out] nLim Number of cells where beta < 1.
671 * @param[out] betaMin Global minimum beta value.
672 * @param[in] flag EvaluateURecBeta_DEFAULT or EvaluateURecBeta_COMPRESS_TO_MEAN.
673 */
674 void EvaluateURecBeta(
675 ArrayDOFV<nVarsFixed> &u,
676 ArrayRECV<nVarsFixed> &uRec,
677 ArrayDOFV<1> &uRecBeta, index &nLim, real &betaMin, int flag);
678
679 /**
680 * @brief Assert that all cell-mean values are physically realizable.
681 *
682 * Checks that density > 0 and internal energy > 0 for all cells.
683 *
684 * @param[in] u Cell-centered conservative DOFs.
685 * @param[in] panic If true, abort on first violation; otherwise just report.
686 * @return true if all cells pass the check.
687 */
688 bool AssertMeanValuePP(
689 ArrayDOFV<nVarsFixed> &u, bool panic);
690
691 static const int EvaluateCellRHSAlpha_DEFAULT = 0x00; ///< Default alpha evaluation mode.
692 static const int EvaluateCellRHSAlpha_MIN_IF_NOT_ONE = 0x01; ///< Take min(alpha, prev) only if prev != 1.
693 static const int EvaluateCellRHSAlpha_MIN_ALL = 0x02; ///< Always take min(alpha, prev).
694 /**
695 * @brief Compute the positivity-preserving RHS scaling factor (alpha) per cell.
696 *
697 * Determines the maximum safe scaling alpha in [0,1] such that
698 * u + alpha * res remains physically realizable.
699 *
700 * @param[in] u Cell-centered conservative DOFs.
701 * @param[in] uRec Reconstruction coefficients.
702 * @param[in] uRecBeta Per-cell reconstruction beta from PP limiter.
703 * @param res Incremental residual (the RHS increment to scale).
704 * @param[out] cellRHSAlpha Per-cell alpha factor.
705 * @param[out] nLim Number of cells where alpha < 1.
706 * @param[out] alphaMin Global minimum alpha.
707 * @param[in] relax Relaxation factor applied to alpha.
708 * @param[in] compress Compression mode (1=compress, 0=clamp).
709 * @param[in] flag EvaluateCellRHSAlpha_* mode flag.
710 */
711 void EvaluateCellRHSAlpha(
712 ArrayDOFV<nVarsFixed> &u,
713 ArrayRECV<nVarsFixed> &uRec,
714 ArrayDOFV<1> &uRecBeta,
715 ArrayDOFV<nVarsFixed> &res,
716 ArrayDOFV<1> &cellRHSAlpha, index &nLim, real &alphaMin,
717 real relax, int compress = 1,
718 int flag = 0);
719
720 /**
721 * @brief Expand a previously computed cellRHSAlpha toward 1 where safe.
722 *
723 * @param res Incremental residual fixed previously.
724 * @param cellRHSAlpha Limiting factor evaluated previously (expanded in-place).
725 * @param[out] nLim Number of cells still limited after expansion.
726 * @param[out] alphaMin Global minimum alpha after expansion.
727 */
728 void EvaluateCellRHSAlphaExpansion(
729 ArrayDOFV<nVarsFixed> &u,
730 ArrayRECV<nVarsFixed> &uRec,
731 ArrayDOFV<1> &uRecBeta,
732 ArrayDOFV<nVarsFixed> &res,
733 ArrayDOFV<1> &cellRHSAlpha, index &nLim, real alphaMin);
734
735 /// @brief Smooth the local time step across neighboring cells.
736 void MinSmoothDTau(
737 ArrayDOFV<1> &dTau, ArrayDOFV<1> &dTauNew);
738
739 /******************************************************/
740
741 /**
742 * @brief Compute effective molecular viscosity using the configured viscosity model.
743 *
744 * Supports constant viscosity (muModel=0), Sutherland's law (muModel=1),
745 * and density-proportional viscosity (muModel=2).
746 *
747 * @param U Conservative state vector.
748 * @param T Temperature.
749 * @return Effective molecular dynamic viscosity.
750 */
751 real muEff(const TU &U, real T) // TODO: more than sutherland law
752 {
753
754 switch (settings.idealGasProperty.muModel)
755 {
756 case 0:
757 return settings.idealGasProperty.muGas;
758 case 1:
759 {
760 real TRel = T / settings.idealGasProperty.TRef;
761 return settings.idealGasProperty.muGas *
762 TRel * std::sqrt(TRel) *
763 (settings.idealGasProperty.TRef + settings.idealGasProperty.CSutherland) /
764 (T + settings.idealGasProperty.CSutherland);
765 }
766 break;
767 case 2:
768 {
769 return settings.idealGasProperty.muGas * U(0);
770 }
771 break;
772 default:
773 DNDS_assert_info(false, "No such muModel");
774 }
775 return std::nan("0");
776 }
777
778 /**
779 * @brief Compute turbulent eddy viscosity at a face.
780 *
781 * Dispatches to the appropriate RANS model (SA, k-omega SST, k-omega Wilcox,
782 * or Realizable k-epsilon) based on settings.ransModel.
783 *
784 * @param uMean Cell-mean conservative state.
785 * @param GradUMeanXY Gradient of conservative variables in physical coordinates.
786 * @param muRef Reference dynamic viscosity scaling.
787 * @param muf Molecular (physical) viscosity at the face.
788 * @param iFace Face index (for wall distance lookup).
789 * @return Turbulent eddy viscosity mu_t.
790 */
791 real getMuTur(const TU &uMean, const TDiffU &GradUMeanXY, real muRef, real muf, index iFace)
792 {
793 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
794 real muTur = 0;
795 if constexpr (Traits::hasSA)
796 {
797 real Chi = uMean(I4 + 1) * muRef / muf;
798#ifdef USE_NS_SA_NEGATIVE_MODEL
799 if (Chi < 10)
800 Chi = 0.05 * std::log(1 + std::exp(20 * Chi));
801#endif
802 real Chi3 = cube(Chi);
803 real fnu1 = Chi3 / (Chi3 + std::pow(RANS::SA::cnu1, 3));
804 muTur = muf * std::max((Chi * fnu1), 0.0);
805 }
806 if constexpr (Traits::has2EQ)
807 {
808 real mut = 0;
809 if (settings.ransModel == RANSModel::RANS_KOSST)
810 mut = RANS::GetMut_SST<dim>(uMean, GradUMeanXY, muf, dWallFace[iFace]);
811 else if (settings.ransModel == RANSModel::RANS_KOWilcox)
812 mut = RANS::GetMut_KOWilcox<dim>(uMean, GradUMeanXY, muf, dWallFace[iFace]);
813 else if (settings.ransModel == RANSModel::RANS_RKE)
814 mut = RANS::GetMut_RealizableKe<dim>(uMean, GradUMeanXY, muf, dWallFace[iFace]);
815 muTur = mut;
816 }
817 return muTur;
818 }
819
820 /**
821 * @brief Compute viscous flux contribution from turbulence model variables.
822 *
823 * Adds the turbulent diffusion flux for SA (nuTilde) or two-equation (k, omega/epsilon)
824 * variables, and the turbulent normal-stress correction to the momentum and energy fluxes.
825 *
826 * @param UMeanXYC Cell-mean state in physical coordinates.
827 * @param DiffUxyPrimC Gradient of primitive variables in physical coordinates.
828 * @param muRef Reference viscosity scaling.
829 * @param mufPhy Physical molecular viscosity.
830 * @param muTur Turbulent eddy viscosity.
831 * @param uNormC Face outward unit normal.
832 * @param iFace Face index (for wall distance lookup).
833 * @param VisFlux Viscous flux vector (accumulated in-place).
834 */
835 void visFluxTurVariable(const TU &UMeanXYC, const TDiffU &DiffUxyPrimC,
836 real muRef, real mufPhy, real muTur, const TVec &uNormC, index iFace, TU &VisFlux)
837 {
838 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
839 if constexpr (Traits::hasSA)
840 {
841 real sigma = RANS::SA::sigma;
842 real fn = 1;
843#ifdef USE_NS_SA_NEGATIVE_MODEL
844 if (UMeanXYC(I4 + 1) < 0)
845 {
846 real Chi = UMeanXYC(I4 + 1) * muRef / mufPhy;
847 fn = (RANS::SA::cn1 + std::pow(Chi, 3)) / (RANS::SA::cn1 - std::pow(Chi, 3));
848 }
849#endif
850 VisFlux(I4 + 1) = DiffUxyPrimC(Seq012, {I4 + 1}).dot(uNormC) * (mufPhy + UMeanXYC(I4 + 1) * muRef * fn) / sigma;
851
852 real tauPressure = DiffUxyPrimC(Seq012, Seq123).trace() * (2. / 3.) * (muTur); //! SA's normal stress
853 tauPressure *= 0; // !not standard SA, abandoning
854 VisFlux(Seq123) -= tauPressure * uNormC;
855 VisFlux(I4) -= tauPressure * UMeanXYC(Seq123).dot(uNormC) / UMeanXYC(0);
856 }
857 if constexpr (Traits::has2EQ)
858 {
859 if (settings.ransModel == RANSModel::RANS_KOSST)
860 RANS::GetVisFlux_SST<dim>(UMeanXYC, DiffUxyPrimC, uNormC, muTur, dWallFace[iFace], mufPhy, VisFlux);
861 else if (settings.ransModel == RANSModel::RANS_KOWilcox)
862 RANS::GetVisFlux_KOWilcox<dim>(UMeanXYC, DiffUxyPrimC, uNormC, muTur, dWallFace[iFace], mufPhy, VisFlux);
863 else if (settings.ransModel == RANSModel::RANS_RKE)
864 RANS::GetVisFlux_RealizableKe<dim>(UMeanXYC, DiffUxyPrimC, uNormC, muTur, dWallFace[iFace], mufPhy, VisFlux);
865 }
866 }
867
868 /**
869 * @brief Transform a conservative state vector from cell frame to face frame for periodic BCs.
870 *
871 * Applies the periodic rotation/translation to the momentum components when
872 * the face is a periodic boundary and the cell is on the donor side.
873 *
874 * @param[in,out] u Conservative state vector (modified in-place).
875 * @param[in] iFace Face index.
876 * @param[in] iCell Cell index.
877 * @param[in] if2c Face-to-cell local index (0=back, 1=front, <0=auto-detect).
878 */
879 void UFromCell2Face(TU &u, index iFace, index iCell, rowsize if2c)
880 {
881 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
882 if (!mesh->isPeriodic)
883 return;
884 auto faceID = mesh->GetFaceZone(iFace);
885 if (!Geom::FaceIDIsPeriodic(faceID))
886 return;
887 if (if2c < 0)
888 if2c = vfv->CellIsFaceBack(iCell, iFace) ? 0 : 1;
889 if (if2c == 1 && Geom::FaceIDIsPeriodicMain(faceID))
890 u(Seq123) = mesh->periodicInfo.TransVectorBack(Eigen::Vector<real, dim>{u(Seq123)}, faceID);
891 if (if2c == 1 && Geom::FaceIDIsPeriodicDonor(faceID))
892 u(Seq123) = mesh->periodicInfo.TransVector(Eigen::Vector<real, dim>{u(Seq123)}, faceID);
893 }
894
895 /// @brief Inverse of UFromCell2Face: transform from face frame back to cell frame.
896 void UFromFace2Cell(TU &u, index iFace, index iCell, rowsize if2c)
897 {
898 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
899 if (!mesh->isPeriodic)
900 return;
901 auto faceID = mesh->GetFaceZone(iFace);
902 if (!Geom::FaceIDIsPeriodic(faceID))
903 return;
904 if (if2c < 0)
905 if2c = vfv->CellIsFaceBack(iCell, iFace) ? 0 : 1;
906 if (if2c == 1 && Geom::FaceIDIsPeriodicMain(faceID))
907 u(Seq123) = mesh->periodicInfo.TransVector(Eigen::Vector<real, dim>{u(Seq123)}, faceID);
908 if (if2c == 1 && Geom::FaceIDIsPeriodicDonor(faceID))
909 u(Seq123) = mesh->periodicInfo.TransVectorBack(Eigen::Vector<real, dim>{u(Seq123)}, faceID);
910 }
911
912 /// @brief Transform a state from a neighbor cell across a periodic face.
913 void UFromOtherCell(TU &u, index iFace, index iCell, index iCellOther, rowsize if2c)
914 {
915 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
916 auto faceID = mesh->GetFaceZone(iFace);
917 if (if2c < 0)
918 if2c = vfv->CellIsFaceBack(iCell, iFace) ? 0 : 1;
919 mesh->CellOtherCellPeriodicHandle(
920 iFace, if2c,
921 [&]()
922 { u(Seq123) = mesh->periodicInfo.TransVector(Eigen::Vector<real, dim>{u(Seq123)}, faceID); },
923 [&]()
924 { u(Seq123) = mesh->periodicInfo.TransVectorBack(Eigen::Vector<real, dim>{u(Seq123)}, faceID); });
925 }
926
927 /// @brief Transform a gradient tensor from cell frame to face frame for periodic BCs.
928 void DiffUFromCell2Face(TDiffU &u, index iFace, index iCell, rowsize if2c, bool reverse = false)
929 {
930 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
931 if (!mesh->isPeriodic)
932 return;
933 auto faceID = mesh->GetFaceZone(iFace);
934 if (!Geom::FaceIDIsPeriodic(faceID))
935 return;
936 if (if2c < 0)
937 if2c = vfv->CellIsFaceBack(iCell, iFace) ? 0 : 1;
938 if ((if2c == 1 && Geom::FaceIDIsPeriodicMain(faceID) && !reverse) ||
939 (if2c == 1 && Geom::FaceIDIsPeriodicDonor(faceID) && reverse))
940 {
941 u(Seq012, EigenAll) = mesh->periodicInfo.TransVectorBack<dim, nVarsFixed>(u(Seq012, EigenAll), faceID);
942 u(EigenAll, Seq123) = mesh->periodicInfo.TransVectorBack<dim, dim>(u(EigenAll, Seq123).transpose(), faceID).transpose();
943 }
944 if ((if2c == 1 && Geom::FaceIDIsPeriodicDonor(faceID) && !reverse) ||
945 (if2c == 1 && Geom::FaceIDIsPeriodicMain(faceID) && reverse))
946 {
947 u(Seq012, EigenAll) = mesh->periodicInfo.TransVector<dim, nVarsFixed>(u(Seq012, EigenAll), faceID);
948 u(EigenAll, Seq123) = mesh->periodicInfo.TransVector<dim, dim>(u(EigenAll, Seq123).transpose(), faceID).transpose();
949 }
950 }
951
952 /**
953 * @brief Compute the numerical flux at a face (batched over quadrature points).
954 *
955 * Evaluates the Riemann-solver-based inviscid flux and (optionally) the viscous
956 * flux at all quadrature points of a face simultaneously. Supports Roe, HLLC,
957 * Lax-Friedrichs, and other Riemann solvers selected by rsType.
958 *
959 * @param[in] ULxy Left state at quadrature points (batched).
960 * @param[in] URxy Right state at quadrature points (batched).
961 * @param[in] ULMeanXy Left cell mean state.
962 * @param[in] URMeanXy Right cell mean state.
963 * @param[in] DiffUxy Gradient of conservative variables at quad points.
964 * @param[in] DiffUxyPrim Gradient of primitive variables at quad points.
965 * @param[in] unitNorm Face outward unit normals at quad points.
966 * @param[in] vg Grid velocity at quad points (for ALE / rotating frame).
967 * @param[in] unitNormC Face-center outward unit normal.
968 * @param[in] vgC Grid velocity at face center.
969 * @param[out] FLfix Left-biased flux correction (batched).
970 * @param[out] FRfix Right-biased flux correction (batched).
971 * @param[out] finc Numerical flux increment (batched).
972 * @param[out] lam0V Eigenvalue |u·n| at quad points.
973 * @param[out] lam123V Eigenvalue |u·n + a| at quad points.
974 * @param[out] lam4V Eigenvalue |u·n - a| at quad points.
975 * @param[in] btype Boundary type (UnInitIndex for internal faces).
976 * @param[in] rsType Riemann solver type (Roe, HLLC, LF, etc.).
977 * @param[in] iFace Face index.
978 * @param[in] ignoreVis If true, skip viscous flux computation.
979 */
980 void fluxFace(
981 const TU_Batch &ULxy,
982 const TU_Batch &URxy,
983 const TU &ULMeanXy,
984 const TU &URMeanXy,
985 const TDiffU_Batch &DiffUxy,
986 const TDiffU_Batch &DiffUxyPrim,
987 const TVec_Batch &unitNorm,
988 const TVec_Batch &vg,
989 const TVec &unitNormC,
990 const TVec &vgC,
991 TU_Batch &FLfix,
992 TU_Batch &FRfix,
993 TU_Batch &finc,
994 TReal_Batch &lam0V, TReal_Batch &lam123V, TReal_Batch &lam4V,
995 Geom::t_index btype,
996 typename Gas::RiemannSolverType rsType,
997 index iFace, bool ignoreVis);
998
999 /**
1000 * @brief Compute the source term at a cell quadrature point.
1001 *
1002 * Evaluates body force, axisymmetric, rotating-frame, and RANS turbulence
1003 * model source terms. Optionally computes the source Jacobian for implicit methods.
1004 *
1005 * @param[in] UMeanXy Conservative state at the quadrature point.
1006 * @param[in] DiffUxy Gradient of conservative variables.
1007 * @param[in] pPhy Physical coordinates of the quadrature point.
1008 * @param[out] jacobian Source Jacobian matrix (populated if Mode=1 or 2).
1009 * @param[in] iCell Cell index.
1010 * @param[in] ig Quadrature point index within the cell.
1011 * @param[in] Mode 0=source vector only, 1=diagonal Jacobian, 2=full Jacobian.
1012 * @return Source term vector.
1013 */
1014 TU source(
1015 const TU &UMeanXy,
1016 const TDiffU &DiffUxy,
1017 const Geom::tPoint &pPhy,
1018 TJacobianU &jacobian,
1019 index iCell,
1020 index ig,
1021 int Mode) // mode =0: source; mode = 1, diagJacobi; mode = 2,
1022 ;
1023
1024 /**
1025 * @brief inviscid flux approx jacobian (flux term not reconstructed / no riemann)
1026 *
1027 */
1028 auto fluxJacobian0_Right(
1029 TU &UR,
1030 const TVec &uNorm,
1031 Geom::t_index btype)
1032 {
1033 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
1034 DNDS_assert(dim == 3); // only for 3D!!!!!!!!
1035 const TU &U = UR;
1036 const TVec &n = uNorm;
1037
1038 real rhoun = n.dot(U({1, 2, 3}));
1039 real rhousqr = U({1, 2, 3}).squaredNorm();
1040 real gamma = settings.idealGasProperty.gamma;
1041 TJacobianU subFdU;
1042 subFdU.resize(nVars, nVars);
1043
1044 subFdU.setZero();
1045 subFdU(0, 1) = n(1 - 1);
1046 subFdU(0, 2) = n(2 - 1);
1047 subFdU(0, 3) = n(3 - 1);
1048 subFdU(1, 0) = -1.0 / (U(1 - 1) * U(1 - 1)) * U(2 - 1) * rhoun + (1.0 / (U(1 - 1) * U(1 - 1)) * n(1 - 1) * (gamma - 1.0) * (rhousqr - U(1 - 1) * U(5 - 1) * 2.0)) / 2.0 + (U(5 - 1) * n(1 - 1) * (gamma - 1.0)) / U(1 - 1);
1049 subFdU(1, 1) = (rhoun + U(2 - 1) * n(1 - 1) * 2.0 - U(2 - 1) * gamma * n(1 - 1)) / U(1 - 1);
1050 subFdU(1, 2) = (U(2 - 1) * n(2 - 1)) / U(1 - 1) - (U(3 - 1) * n(1 - 1) * (gamma - 1.0)) / U(1 - 1);
1051 subFdU(1, 3) = (U(2 - 1) * n(3 - 1)) / U(1 - 1) - (U(4 - 1) * n(1 - 1) * (gamma - 1.0)) / U(1 - 1);
1052 subFdU(1, 4) = n(1 - 1) * (gamma - 1.0);
1053 subFdU(2, 0) = -1.0 / (U(1 - 1) * U(1 - 1)) * U(3 - 1) * rhoun + (1.0 / (U(1 - 1) * U(1 - 1)) * n(2 - 1) * (gamma - 1.0) * (rhousqr - U(1 - 1) * U(5 - 1) * 2.0)) / 2.0 + (U(5 - 1) * n(2 - 1) * (gamma - 1.0)) / U(1 - 1);
1054 subFdU(2, 1) = (U(3 - 1) * n(1 - 1)) / U(1 - 1) - (U(2 - 1) * n(2 - 1) * (gamma - 1.0)) / U(1 - 1);
1055 subFdU(2, 2) = (rhoun + U(3 - 1) * n(2 - 1) * 2.0 - U(3 - 1) * gamma * n(2 - 1)) / U(1 - 1);
1056 subFdU(2, 3) = (U(3 - 1) * n(3 - 1)) / U(1 - 1) - (U(4 - 1) * n(2 - 1) * (gamma - 1.0)) / U(1 - 1);
1057 subFdU(2, 4) = n(2 - 1) * (gamma - 1.0);
1058 subFdU(3, 0) = -1.0 / (U(1 - 1) * U(1 - 1)) * U(4 - 1) * rhoun + (1.0 / (U(1 - 1) * U(1 - 1)) * n(3 - 1) * (gamma - 1.0) * (rhousqr - U(1 - 1) * U(5 - 1) * 2.0)) / 2.0 + (U(5 - 1) * n(3 - 1) * (gamma - 1.0)) / U(1 - 1);
1059 subFdU(3, 1) = (U(4 - 1) * n(1 - 1)) / U(1 - 1) - (U(2 - 1) * n(3 - 1) * (gamma - 1.0)) / U(1 - 1);
1060 subFdU(3, 2) = (U(4 - 1) * n(2 - 1)) / U(1 - 1) - (U(3 - 1) * n(3 - 1) * (gamma - 1.0)) / U(1 - 1);
1061 subFdU(3, 3) = (rhoun + U(4 - 1) * n(3 - 1) * 2.0 - U(4 - 1) * gamma * n(3 - 1)) / U(1 - 1);
1062 subFdU(3, 4) = n(3 - 1) * (gamma - 1.0);
1063 subFdU(4, 0) = 1.0 / (U(1 - 1) * U(1 - 1) * U(1 - 1)) * rhoun * (-rhousqr + (U(2 - 1) * U(2 - 1)) * gamma + (U(3 - 1) * U(3 - 1)) * gamma + (U(4 - 1) * U(4 - 1)) * gamma - U(1 - 1) * U(5 - 1) * gamma);
1064 subFdU(4, 1) = 1.0 / (U(1 - 1) * U(1 - 1)) * n(1 - 1) * (-rhousqr + (U(2 - 1) * U(2 - 1)) * gamma + (U(3 - 1) * U(3 - 1)) * gamma + (U(4 - 1) * U(4 - 1)) * gamma - U(1 - 1) * U(5 - 1) * gamma * 2.0) * (-1.0 / 2.0) - 1.0 / (U(1 - 1) * U(1 - 1)) * U(2 - 1) * rhoun * (gamma - 1.0);
1065 subFdU(4, 2) = 1.0 / (U(1 - 1) * U(1 - 1)) * n(2 - 1) * (-rhousqr + (U(2 - 1) * U(2 - 1)) * gamma + (U(3 - 1) * U(3 - 1)) * gamma + (U(4 - 1) * U(4 - 1)) * gamma - U(1 - 1) * U(5 - 1) * gamma * 2.0) * (-1.0 / 2.0) - 1.0 / (U(1 - 1) * U(1 - 1)) * U(3 - 1) * rhoun * (gamma - 1.0);
1066 subFdU(4, 3) = 1.0 / (U(1 - 1) * U(1 - 1)) * n(3 - 1) * (-rhousqr + (U(2 - 1) * U(2 - 1)) * gamma + (U(3 - 1) * U(3 - 1)) * gamma + (U(4 - 1) * U(4 - 1)) * gamma - U(1 - 1) * U(5 - 1) * gamma * 2.0) * (-1.0 / 2.0) - 1.0 / (U(1 - 1) * U(1 - 1)) * U(4 - 1) * rhoun * (gamma - 1.0);
1067 subFdU(4, 4) = (gamma * rhoun) / U(1 - 1);
1068
1069 real un = rhoun / U(0);
1070
1071 if constexpr (Traits::hasSA)
1072 {
1073 subFdU(5, 5) = un;
1074 subFdU(5, 0) = -un * U(5) / U(0);
1075 subFdU(5, 1) = n(0) * U(5) / U(0);
1076 subFdU(5, 2) = n(1) * U(5) / U(0);
1077 subFdU(5, 3) = n(2) * U(5) / U(0);
1078 }
1079 if constexpr (Traits::has2EQ)
1080 {
1081 subFdU(5, 5) = un;
1082 subFdU(5, 0) = -un * U(5) / U(0);
1083 subFdU(5, 1) = n(0) * U(5) / U(0);
1084 subFdU(5, 2) = n(1) * U(5) / U(0);
1085 subFdU(5, 3) = n(2) * U(5) / U(0);
1086 subFdU(6, 6) = un;
1087 subFdU(6, 0) = -un * U(6) / U(0);
1088 subFdU(6, 1) = n(0) * U(6) / U(0);
1089 subFdU(6, 2) = n(1) * U(6) / U(0);
1090 subFdU(6, 3) = n(2) * U(6) / U(0);
1091 }
1092 return subFdU;
1093 }
1094
1095 /**
1096 * @brief inviscid flux approx jacobian (flux term not reconstructed / no riemann)
1097 * if lambdaMain == veryLargeReal, then use lambda0~4 for roe-flux type jacobian
1098 *
1099 */
1100 TU fluxJacobian0_Right_Times_du(
1101 const TU &U,
1102 const TU &UOther,
1103 const TVec &n,
1104 const TVec &vg,
1105 Geom::t_index btype,
1106 const TU &dU,
1107 real lambdaMain, real lambdaC, real lambdaVis,
1108 real lambda0, real lambda123, real lambda4,
1109 int useRoeTerm, int incFsign = -1, int omitF = 0)
1110 {
1111 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
1112 real gamma = settings.idealGasProperty.gamma;
1113 TVec velo = U(Seq123) / U(0);
1114 real p, H, asqr;
1115 Gas::IdealGasThermal(U(I4), U(0), velo.squaredNorm(), gamma, p, asqr, H);
1116 TVec dVelo;
1117 real dp;
1118 Gas::IdealGasUIncrement<dim>(U, dU, velo, gamma, dVelo, dp);
1119 TU dF(U.size());
1120 if (omitF == 0)
1121 Gas::GasInviscidFluxFacialIncrement<dim>(
1122 U, dU,
1123 n,
1124 velo, dVelo, vg,
1125 dp, p,
1126 dF);
1127 else
1128 dF.setZero();
1129 if (useRoeTerm == 0)
1130 dF += incFsign * lambdaMain * dU;
1131 else
1132 {
1133 TVec veloRoe;
1134 real vsqrRoe{0}, aRoe{0}, asqrRoe{0}, HRoe{0};
1135 TU uMean(U.size());
1136 Gas::GetRoeAverage<dim>(U, UOther, gamma, veloRoe, vsqrRoe, aRoe, asqrRoe, HRoe, uMean);
1137 {
1138 // TVec dVeloRoe;
1139 // real dpRoe;
1140 // Gas::IdealGasUIncrement<dim>(uMean, dU, veloRoe, gamma, dVeloRoe, dpRoe);
1141 // Gas::GasInviscidFluxFacialIncrement<dim>(
1142 // uMean, dU,
1143 // n,
1144 // veloRoe, dVeloRoe, vg,
1145 // dp, p,
1146 // dF);
1147 }
1148 // lambda0 = lambda123 = lambda4 = aRoe + std::abs(veloRoe.dot(n));
1149 Gas::RoeFluxIncFDiff<dim>(dU, n, veloRoe, vsqrRoe, aRoe, asqrRoe, HRoe,
1150 incFsign * lambda0, incFsign * lambda123, incFsign * lambda4, gamma,
1151 dF);
1152 dF += incFsign * lambdaVis * dU;
1153 }
1154 //! now dF(U, dU) (GasInviscidFluxFacialIncrement) part actually treats the SeqI52Last part (as passive scalars)
1155 //! the RANS models and eulerEX use the same transport jacobian form
1156
1157 return dF;
1158 }
1159
1160 TJacobianU fluxJacobian0_Right_Times_du_AsMatrix(
1161 const TU &U,
1162 const TU &UOther,
1163 const TVec &n,
1164 const TVec &vg,
1165 Geom::t_index btype,
1166 real lambdaMain, real lambdaC, real lambdaVis,
1167 real lambda0, real lambda123, real lambda4,
1168 int useRoeTerm, int incFsign = -1, int omitF = 0)
1169 { // TODO: optimize this
1170 TJacobianU J;
1171 J.resize(nVars, nVars);
1172 J.setIdentity();
1173 for (int i = 0; i < nVars; i++)
1174 J(EigenAll, i) = fluxJacobian0_Right_Times_du(
1175 U, UOther, n, vg, btype, J(EigenAll, i),
1176 lambdaMain, lambdaC, lambdaVis, lambda0, lambda123, lambda4,
1177 useRoeTerm, incFsign, omitF);
1178 // TODO: for eulerEX, use scalar for SeqI52Last part
1179 return J;
1180 }
1181
1182 TU fluxJacobianC_Right_Times_du(
1183 const TU &U,
1184 const TVec &n,
1185 const TVec &vg,
1186 Geom::t_index btype,
1187 const TU &dU)
1188 {
1189 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
1190 real gamma = settings.idealGasProperty.gamma;
1191 TVec velo = U(Seq123) / U(0);
1192 real p, H, asqr;
1193 Gas::IdealGasThermal(U(I4), U(0), velo.squaredNorm(), gamma, p, asqr, H);
1194 TVec dVelo;
1195 real dp;
1196 Gas::IdealGasUIncrement<dim>(U, dU, velo, gamma, dVelo, dp);
1197 TU dF(U.size());
1198 Gas::GasInviscidFluxFacialIncrement<dim>(
1199 U, dU,
1200 n,
1201 velo, dVelo, vg,
1202 dp, p,
1203 dF);
1204 return dF;
1205 }
1206
1207 /// @brief Get the grid velocity at a face quadrature point (rotating frame).
1208 TVec GetFaceVGrid(index iFace, index iG)
1209 {
1210 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
1211 TVec ret;
1212 ret.setZero();
1213#ifdef USE_ABS_VELO_IN_ROTATION
1214 if (settings.frameConstRotation.enabled)
1215 ret += settings.frameConstRotation.vOmega().cross(vfv->GetFaceQuadraturePPhys(iFace, iG) - settings.frameConstRotation.center)(Seq012);
1216#endif
1217 return ret;
1218 }
1219
1220 /// @brief Get the grid velocity at a face quadrature point (with explicit physical point).
1221 TVec GetFaceVGrid(index iFace, index iG, const Geom::tPoint &pPhy)
1222 {
1223 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
1224 TVec ret;
1225 ret.setZero();
1226#ifdef USE_ABS_VELO_IN_ROTATION
1227 if (settings.frameConstRotation.enabled)
1228 ret += settings.frameConstRotation.vOmega().cross(
1229 (iG >= -1 ? vfv->GetFaceQuadraturePPhys(iFace, iG) : pPhy) - settings.frameConstRotation.center)(Seq012);
1230#endif
1231 return ret;
1232 }
1233
1234 /// @brief Get the grid velocity at a face quadrature point from a specific cell's perspective.
1235 TVec GetFaceVGridFromCell(index iFace, index iCell, int if2c, index iG)
1236 {
1237 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
1238 TVec ret;
1239 ret.setZero();
1240#ifdef USE_ABS_VELO_IN_ROTATION
1241 if (settings.frameConstRotation.enabled)
1242 ret += settings.frameConstRotation.vOmega().cross(vfv->GetFaceQuadraturePPhysFromCell(iFace, iCell, if2c, iG) - settings.frameConstRotation.center)(Seq012);
1243#endif
1244 return ret;
1245 }
1246
1247 /// @brief Transform momentum in-place between inertial and rotating frame (velocity only).
1248 void TransformVelocityRotatingFrame(TU &U, const Geom::tPoint &pPhysics, int direction)
1249 {
1250 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
1251 U(Seq123) += direction * settings.frameConstRotation.vOmega().cross(pPhysics - settings.frameConstRotation.center)(Seq012) * U(0);
1252 }
1253
1254 /// @brief Transform full conservative state (momentum + total energy) between frames (relative velocity formulation).
1255 void TransformURotatingFrame(TU &U, const Geom::tPoint &pPhysics, int direction)
1256 {
1257 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
1258#ifndef USE_ABS_VELO_IN_ROTATION
1259 U(I4) -= U(Seq123).squaredNorm() / (2 * U(0));
1260 U(Seq123) += direction * settings.frameConstRotation.vOmega().cross(pPhysics - settings.frameConstRotation.center)(Seq012) * U(0);
1261 U(I4) += U(Seq123).squaredNorm() / (2 * U(0));
1262#endif
1263 }
1264
1265 /// @brief Transform full conservative state for the absolute-velocity rotating frame formulation.
1266 void TransformURotatingFrame_ABS_VELO(TU &U, const Geom::tPoint &pPhysics, int direction)
1267 {
1268 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
1269#ifdef USE_ABS_VELO_IN_ROTATION
1270 U(I4) -= U(Seq123).squaredNorm() / (2 * U(0));
1271 U(Seq123) += direction * settings.frameConstRotation.vOmega().cross(pPhysics - settings.frameConstRotation.center)(Seq012) * U(0);
1272 U(I4) += U(Seq123).squaredNorm() / (2 * U(0));
1273#endif
1274 }
1275
1276 /// @brief Update boundary anchor point recorders from current solution.
1277 void updateBCAnchors(ArrayDOFV<nVarsFixed> &u, ArrayRECV<nVarsFixed> &uRec)
1278 {
1279 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
1280 for (Geom::t_index i = Geom::BC_ID_DEFAULT_MAX; i < pBCHandler->size(); i++) // init code, consider adding to ctor
1281 {
1282 if (pBCHandler->GetFlagFromIDSoft(i, "anchorOpt") == 0)
1283 continue;
1284 if (!anchorRecorders.count(i))
1285 anchorRecorders.emplace(std::make_pair(i, AnchorPointRecorder<nVarsFixed>(mesh->getMPI())));
1286 }
1287 for (auto &v : anchorRecorders)
1288 v.second.Reset();
1289 for (index iBnd = 0; iBnd < mesh->NumBnd(); iBnd++)
1290 {
1291 index iFace = mesh->bnd2faceV.at(iBnd);
1292 if (iFace < 0) // remember that some iBnd do not have iFace (for periodic case)
1293 continue;
1294 auto f2c = mesh->face2cell[iFace];
1295 auto gFace = vfv->GetFaceQuad(iFace);
1296
1297 Geom::Elem::SummationNoOp noOp;
1298 auto faceBndID = mesh->GetFaceZone(iFace);
1299 auto faceBCType = pBCHandler->GetTypeFromID(faceBndID);
1300
1301 if (pBCHandler->GetFlagFromIDSoft(faceBndID, "anchorOpt") == 0)
1302 continue;
1303 gFace.IntegrationSimple(
1304 noOp,
1305 [&](auto finc, int iG)
1306 {
1307 TU ULxy = u[f2c[0]];
1308 ULxy += (vfv->GetIntPointDiffBaseValue(f2c[0], iFace, 0, iG, std::array<int, 1>{0}, 1) *
1309 uRec[f2c[0]])
1310 .transpose();
1311 this->UFromCell2Face(ULxy, iFace, f2c[0], 0);
1312 real dist = vfv->GetFaceQuadraturePPhys(iFace, iG).norm();
1313 if (pBCHandler->GetValueExtraFromID(faceBndID).size() >= 3)
1314 {
1315 Geom::tPoint vOrig = pBCHandler->GetValueExtraFromID(faceBndID)({0, 1, 2});
1316 dist = (vfv->GetFaceQuadraturePPhys(iFace, iG) - vOrig).norm();
1317 }
1318 anchorRecorders.at(faceBndID).AddAnchor(ULxy, dist);
1319 });
1320 }
1321 for (auto &v : anchorRecorders)
1322 v.second.ObtainAnchorMPI();
1323 }
1324
1325 /// @brief Update boundary 1-D profiles from the current solution.
1326 void updateBCProfiles(ArrayDOFV<nVarsFixed> &u, ArrayRECV<nVarsFixed> &uRec);
1327
1328 /// @brief Update boundary profiles using radial-equilibrium pressure extrapolation.
1329 void updateBCProfilesPressureRadialEq();
1330
1331 /**
1332 * @brief Generate the ghost (boundary) state for a boundary face.
1333 *
1334 * Dispatches to type-specific BC handlers (far-field, wall, outflow, inflow,
1335 * symmetry, etc.) based on btype. Used by the RHS evaluator to obtain the
1336 * right state at boundary faces for the Riemann solver.
1337 *
1338 * @param[in,out] ULxy Left (interior) state; may be modified for certain BCs.
1339 * @param[in] ULMeanXy Left cell-mean state (base state for linearized mode).
1340 * @param[in] iCell Cell index adjacent to the boundary face.
1341 * @param[in] iFace Face index.
1342 * @param[in] iG Quadrature point index (< -1 for arbitrary location).
1343 * @param[in] uNorm Face outward unit normal.
1344 * @param[in] normBase Orthonormal basis with uNorm as first column.
1345 * @param[in] pPhysics Physical coordinates of the evaluation point.
1346 * @param[in] t Current simulation time.
1347 * @param[in] btype Boundary type ID.
1348 * @param[in] fixUL If true, do not modify the left state.
1349 * @param[in] geomMode Geometry evaluation mode (0=standard).
1350 * @param[in] linMode Linearization mode (0=nonlinear, nonzero=linearized about ULMeanXy).
1351 * @return Ghost (right) state for the Riemann solver.
1352 */
1353 TU generateBoundaryValue(
1354 TU &ULxy, //! warning, possible that UL is also modified
1355 const TU &ULMeanXy,
1356 index iCell, index iFace, int iG,
1357 const TVec &uNorm,
1358 const TMat &normBase,
1359 const Geom::tPoint &pPhysics,
1360 real t,
1361 Geom::t_index btype,
1362 bool fixUL = false,
1363 int geomMode = 0,
1364 int linMode = 0);
1365
1366 private:
1367 /// @name Per-BC-type handlers (called by generateBoundaryValue)
1368 /// @{
1369
1370 /// Characteristic-based far-field / outflow-pressure BC (BCFar, BCOutP, BCSpecial far-field).
1371 TU generateBV_FarField(
1372 TU &ULxy, const TU &ULMeanXy,
1373 index iCell, index iFace, int iG,
1374 const TVec &uNorm, const TMat &normBase,
1375 const Geom::tPoint &pPhysics, real t,
1376 Geom::t_index btype, bool fixUL, int geomMode);
1377
1378 /// Analytical / special far-field BCs (DMR, RT, IV, 2D Riemann, Noh).
1379 TU generateBV_SpecialFar(
1380 TU &ULxy, const TU &ULMeanXy,
1381 index iCell, index iFace, int iG,
1382 const TVec &uNorm, const TMat &normBase,
1383 const Geom::tPoint &pPhysics, real t,
1384 Geom::t_index btype);
1385
1386 /// Inviscid wall / symmetry BC (BCWallInvis, BCSym).
1387 TU generateBV_InviscidWall(
1388 TU &ULxy, const TU &ULMeanXy,
1389 index iCell, index iFace, int iG,
1390 const TVec &uNorm, const TMat &normBase,
1391 const Geom::tPoint &pPhysics, real t,
1392 Geom::t_index btype);
1393
1394 /// Viscous no-slip wall BC (BCWall, BCWallIsothermal), including RANS wall treatment.
1395 TU generateBV_ViscousWall(
1396 TU &ULxy, const TU &ULMeanXy,
1397 index iCell, index iFace, int iG,
1398 const TVec &uNorm, const TMat &normBase,
1399 const Geom::tPoint &pPhysics, real t,
1400 Geom::t_index btype, bool fixUL, int linMode);
1401
1402 /// Supersonic / extrapolation outflow BC (BCOut).
1403 TU generateBV_Outflow(
1404 TU &ULxy, const TU &ULMeanXy,
1405 index iCell, index iFace, int iG,
1406 const TVec &uNorm, const TMat &normBase,
1407 const Geom::tPoint &pPhysics, real t,
1408 Geom::t_index btype);
1409
1410 /// Prescribed inflow BC (BCIn).
1411 TU generateBV_Inflow(
1412 TU &ULxy, const TU &ULMeanXy,
1413 index iCell, index iFace, int iG,
1414 const TVec &uNorm, const TMat &normBase,
1415 const Geom::tPoint &pPhysics, real t,
1416 Geom::t_index btype);
1417
1418 /// Total pressure / total temperature inflow BC (BCInPsTs).
1419 TU generateBV_TotalConditionInflow(
1420 TU &ULxy, const TU &ULMeanXy,
1421 index iCell, index iFace, int iG,
1422 const TVec &uNorm, const TMat &normBase,
1423 const Geom::tPoint &pPhysics, real t,
1424 Geom::t_index btype);
1425 /// @}
1426
1427 public:
1428
1429 /// @brief Write boundary profile data to CSV files (rank 0 only).
1430 void PrintBCProfiles(const std::string &name, ArrayDOFV<nVarsFixed> &u, ArrayRECV<nVarsFixed> &uRec)
1431 {
1432 this->updateBCProfiles(u, uRec);
1433 if (mesh->getMPI().rank != 0)
1434 return; //! only 0 needs to write
1435 for (auto &[id, bcProfile] : profileRecorders)
1436 {
1437 std::string fname = name + "_bc[" + pBCHandler->GetNameFormID(id) + "]_profile.csv";
1438 std::filesystem::path outFile{fname};
1439 std::filesystem::create_directories(outFile.parent_path() / ".");
1440 std::ofstream fout(fname);
1441 DNDS_assert_info(fout, fmt::format("failed to open [{}]", fname));
1442 bcProfile.OutProfileCSV(fout);
1443 }
1444 }
1445
1446 /// @brief Print boundary integration results (force/flux) to console on rank 0.
1447 void ConsoleOutputBndIntegrations()
1448 {
1449 for (auto &i : bndIntegrations)
1450 {
1451 auto intOpt = pBCHandler->GetFlagFromIDSoft(i.first, "integrationOpt");
1452 if (mesh->getMPI().rank == 0)
1453 {
1454 Eigen::VectorFMTSafe<real, Eigen::Dynamic> vPrint = i.second.v;
1455 if (intOpt == 2)
1456 vPrint(Eigen::seq(nVars, nVars + 1)) /= i.second.div;
1457 log() << fmt::format("Bnd [{}] integarted values option [{}] : {:.5e}",
1458 pBCHandler->GetNameFormID(i.first),
1459 intOpt, vPrint.transpose())
1460 << std::endl;
1461 }
1462 }
1463 }
1464
1465 /// @brief Append a line to the per-BC boundary integration CSV log files.
1466 void BndIntegrationLogWriteLine(const std::string &name, index step, index stage, index iter)
1467 {
1468 if (mesh->getMPI().rank != 0)
1469 return; //! only 0 needs to write
1470 for (auto &[id, bndInt] : bndIntegrations)
1471 {
1472 auto intOpt = pBCHandler->GetFlagFromIDSoft(id, "integrationOpt");
1473 if (!bndIntegrationLogs.count(id))
1474 {
1475 std::string fname = name + "_bc[" + pBCHandler->GetNameFormID(id) + "]_integrationLog.csv";
1476 bndIntegrationLogs.emplace(std::make_pair(id, std::ofstream(fname)));
1477 DNDS_assert_info(bndIntegrationLogs.at(id), fmt::format("failed to open [{}]", fname));
1478 bndIntegrationLogs.at(id) << "step, stage, iter";
1479 for (int i = 0; i < bndInt.v.size(); i++)
1480 bndIntegrationLogs.at(id) << ", F" << std::to_string(i);
1481 bndIntegrationLogs.at(id) << "\n";
1482 }
1483 Eigen::Vector<real, Eigen::Dynamic> vPrint = bndInt.v;
1484 if (intOpt == 2)
1485 vPrint(Eigen::seq(nVars, nVars + 1)) /= bndInt.div;
1486 bndIntegrationLogs.at(id) << step << ", " << stage << ", " << iter << std::setprecision(16) << std::scientific;
1487 for (auto &val : vPrint)
1488 bndIntegrationLogs.at(id) << ", " << val;
1489 bndIntegrationLogs.at(id) << std::endl;
1490 }
1491 }
1492
1493 /**
1494 * @brief Query the CL driver for current lift/drag coefficients and update AoA.
1495 *
1496 * Collects boundary force integrals from CL-driven boundaries, projects them
1497 * onto lift/drag directions, and feeds the result to the CLDriver for AoA adaptation.
1498 *
1499 * @param iter Current iteration count.
1500 * @return Tuple of (CL, CD, AoA) after the driver update.
1501 */
1502 std::tuple<real, real, real> CLDriverGetIntegrationUpdate(index iter)
1503 {
1504 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
1505 if (!pCLDriver)
1506 return {0.0, 0.0, 0.0};
1507 TU fluxBndCur;
1508 fluxBndCur.setZero(nVars);
1509 for (auto bcID : cLDriverBndIDs)
1510 {
1511 if (bcID >= Geom::BC_ID_DEFAULT_MAX)
1512 {
1513 fluxBndCur += bndIntegrations.at(bcID).v;
1514 }
1515 else
1516 fluxBndCur += this->fluxWallSum;
1517 }
1518 Geom::tPoint fluxFaceForce;
1519 fluxFaceForce.setZero();
1520 fluxFaceForce(Seq012) = fluxBndCur(Seq123);
1521 fluxFaceForce = pCLDriver->GetAOARotation().transpose() * fluxFaceForce;
1522 real CLCur = fluxFaceForce.dot(pCLDriver->GetCL0Direction()) * pCLDriver->GetForce2CoeffRatio();
1523 real CDCur = fluxFaceForce.dot(pCLDriver->GetCD0Direction()) * pCLDriver->GetForce2CoeffRatio();
1524 pCLDriver->Update(iter, CLCur, mesh->getMPI());
1525 real AOACur = pCLDriver->GetAOA();
1526 return {CLCur, CDCur, AOACur};
1527 }
1528
1529 /**
1530 * @brief Compress a reconstruction increment to maintain positivity.
1531 *
1532 * Given a cell mean and reconstruction increment, returns umean + uRecInc
1533 * with the increment clamped so that density, internal energy, and (for RANS)
1534 * turbulent variables remain non-negative.
1535 *
1536 * @param umean Cell-mean conservative state.
1537 * @param uRecInc Reconstruction polynomial increment at an evaluation point.
1538 * @param compressed Set to true if compression was applied.
1539 * @return Compressed reconstructed state.
1540 */
1541 inline TU CompressRecPart(
1542 const TU &umean,
1543 const TU &uRecInc,
1544 bool &compressed)
1545 {
1546 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
1547 // if (umean(0) + uRecInc(0) < 0)
1548 // {
1549 // std::cout << umean.transpose() << std::endl
1550 // << uRecInc.transpose() << std::endl;
1551 // DNDS_assert(false);
1552 // }
1553 // return umean + uRecInc; // ! no compress shortcut
1554 // return umean; // ! 0th order shortcut
1555
1556 // // * Compress Method
1557 // real compressT = 0.00001;
1558 // real eFixRatio = 0.00001;
1559 // Eigen::Vector<real, 5> ret;
1560
1561 // real compress = 1.0;
1562 // if ((umean(0) + uRecInc(0)) < umean(0) * compressT)
1563 // compress *= umean(0) * (1 - compressT) / uRecInc(0);
1564
1565 // ret = umean + uRecInc * compress;
1566
1567 // real Ek = ret({1, 2, 3}).squaredNorm() * 0.5 / (verySmallReal + ret(0));
1568 // real eT = eFixRatio * Ek;
1569 // real e = ret(4) - Ek;
1570 // if (e < 0)
1571 // e = eT * 0.5;
1572 // else if (e < eT)
1573 // e = (e * e + eT * eT) / (2 * eT);
1574 // ret(4) = e + Ek;
1575 // // * Compress Method
1576
1577 // TU ret = umean + uRecInc;
1578 // real eK = ret(Seq123).squaredNorm() * 0.5 / (verySmallReal + std::abs(ret(0)));
1579 // real e = ret(I4) - eK;
1580 // if (e <= 0 || ret(0) <= 0)
1581 // ret = umean, compressed = true;
1582 // if constexpr (Traits::hasSA)
1583 // if (ret(I4 + 1) < 0)
1584 // ret = umean, compressed = true;
1585
1586 bool rhoFixed = false;
1587 bool eFixed = false;
1588 TU ret = umean + uRecInc;
1589 // do rho fix
1590 if (ret(0) < 0)
1591 {
1592 // rhoFixed = true;
1593 // TVec veloOld = umean(Seq123) / umean(0);
1594 // real eOld = umean(I4) - 0.5 * veloOld.squaredNorm() * umean(0);
1595 // ret(0) = umean(0) * std::exp(uRecInc(0) / (umean(0) + verySmallReal));
1596 // ret(Seq123) = veloOld * ret(0);
1597 // ret(I4) = eOld + 0.5 * veloOld.squaredNorm() * ret(0);
1598
1599 ret = umean;
1600 compressed = true;
1601 }
1602
1603 real eK = ret(Seq123).squaredNorm() * 0.5 / (verySmallReal + ret(0));
1604 real e = ret(I4) - eK;
1605 if (e < 0)
1606 {
1607 // eFixed = true;
1608 // real eOld = umean(I4) - eK;
1609 // real eNew = eOld * std::exp(eOld / (umean(I4) - eK));
1610 // ret(I4) = eNew + eK;
1611
1612 ret = umean;
1613
1614 compressed = true;
1615 }
1616
1617#ifdef USE_NS_SA_NUT_REDUCED_ORDER
1618 if constexpr (Traits::hasSA)
1619 if (ret(I4 + 1) < 0)
1620 ret(I4 + 1) = umean(I4 + 1), compressed = true;
1621#endif
1622 if constexpr (Traits::has2EQ)
1623 {
1624 if (ret(I4 + 1) < 0)
1625 ret(I4 + 1) = umean(I4 + 1), compressed = true;
1626 if (ret(I4 + 2) < 0)
1627 ret(I4 + 2) = umean(I4 + 2), compressed = true;
1628 }
1629
1630 return ret;
1631 }
1632
1633 /**
1634 * @brief Compress a solution increment to maintain positivity.
1635 *
1636 * Given the current state u and an update increment uInc, returns a modified
1637 * increment that ensures u + result has positive density, positive internal
1638 * energy, and (for RANS) non-negative turbulent variables. Uses exponential
1639 * decay clamping for density and a quadratic solve for internal energy.
1640 *
1641 * @param u Current conservative state.
1642 * @param uInc Proposed increment.
1643 * @return Modified (compressed) increment safe to add to u.
1644 */
1645 inline TU CompressInc(
1646 const TU &u,
1647 const TU &uInc)
1648 {
1649 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
1650 real rhoEps = smallReal * settings.refUPrim(0) * 1e-1;
1651 real pEps = smallReal * settings.refUPrim(I4) * 1e-1;
1652 if (settings.ppEpsIsRelaxed)
1653 {
1654 pEps *= verySmallReal, rhoEps *= verySmallReal;
1655 }
1656
1657 TU ret = uInc;
1658
1659 /** A intuitive fix **/ //! need positive perserving technique!
1660 DNDS_assert(u(0) > 0);
1661 if (u(0) + ret(0) <= rhoEps)
1662 {
1663 real declineV = ret(0) / (u(0) + verySmallReal);
1664 real newrho = u(0) * std::exp(declineV);
1665 // newrho = std::max(newrho, rhoEps);
1666 newrho = rhoEps;
1667 ret *= (newrho - u(0)) / (ret(0) - verySmallReal);
1668 // std::cout << (newrho - u(0)) / (ret(0) + verySmallReal) << std::endl;
1669 // DNDS_assert(false);
1670 }
1671 real ekOld = 0.5 * u(Seq123).squaredNorm() / (u(0) + verySmallReal);
1672 real rhoEinternal = u(I4) - ekOld;
1673 DNDS_assert(rhoEinternal > 0);
1674 real ek = 0.5 * (u(Seq123) + ret(Seq123)).squaredNorm() / (u(0) + ret(0) + verySmallReal);
1675 real rhoEinternalNew = u(I4) + ret(I4) - ek;
1676 if (rhoEinternalNew <= pEps)
1677 {
1678 real declineV = (rhoEinternalNew - rhoEinternal) / (rhoEinternal + verySmallReal);
1679 real newrhoEinteralNew = (std::exp(declineV) + verySmallReal) * rhoEinternal;
1680 real gamma = settings.idealGasProperty.gamma;
1681 // newrhoEinteralNew = std::max(pEps / (gamma - 1), newrhoEinteralNew);
1682 newrhoEinteralNew = pEps / (gamma - 1);
1683 real c0 = 2 * u(I4) * u(0) - u(Seq123).squaredNorm() - 2 * u(0) * newrhoEinteralNew;
1684 real c1 = 2 * u(I4) * ret(0) + 2 * u(0) * ret(I4) - 2 * u(Seq123).dot(ret(Seq123)) - 2 * ret(0) * newrhoEinteralNew;
1685 real c2 = 2 * ret(I4) * ret(0) - ret(Seq123).squaredNorm();
1686 real deltaC = sqr(c1) - 4 * c0 * c2;
1687 DNDS_assert(deltaC > 0);
1688 real alphaL = (-std::sqrt(deltaC) - c1) / (2 * c2);
1689 real alphaR = (std::sqrt(deltaC) - c1) / (2 * c2);
1690 // if (c2 > 0)
1691 // DNDS_assert(alphaL > 0);
1692 // DNDS_assert(alphaR > 0);
1693 // DNDS_assert(alphaL < 1);
1694 // if (c2 < 0)
1695 // DNDS_assert(alphaR < 1);
1696 real alpha = std::min((c2 > 0 ? alphaL : alphaL), 1.);
1697 alpha = std::max(0., alpha);
1698 ret *= alpha * (0.99);
1699
1700 real decay = 1 - 1e-1;
1701 for (int iter = 0; iter < 1000; iter++)
1702 {
1703 ek = 0.5 * (u(Seq123) + ret(Seq123)).squaredNorm() / (u(0) + ret(0) + verySmallReal);
1704 if (ret(I4) + u(I4) - ek < newrhoEinteralNew)
1705 ret *= decay, alpha *= decay;
1706 else
1707 break;
1708 }
1709
1710 ek = 0.5 * (u(Seq123) + ret(Seq123)).squaredNorm() / (u(0) + ret(0) + verySmallReal);
1711
1712 if (ret(I4) + u(I4) - ek < newrhoEinteralNew * 0.5)
1713 {
1714 std::cout << std::scientific << std::setprecision(5);
1715 std::cout << u(0) << " " << ret(0) << std::endl;
1716 std::cout << rhoEinternalNew << " " << rhoEinternal << std::endl;
1717 std::cout << declineV << std::endl;
1718 std::cout << newrhoEinteralNew << std::endl;
1719 std::cout << ret(I4) + u(I4) - ek << std::endl;
1720 std::cout << alpha << std::endl;
1721 DNDS_assert(false);
1722 }
1723 }
1724
1725 /** A intuitive fix **/
1726// #define USE_NS_SA_ALLOW_NEGATIVE_MEAN
1727#ifndef USE_NS_SA_ALLOW_NEGATIVE_MEAN
1728 if constexpr (Traits::hasSA)
1729 {
1730 if (u(I4 + 1) + ret(I4 + 1) < 0)
1731 {
1732 // std::cout << "Fixing SA inc " << std::endl;
1733
1734 DNDS_assert(u(I4 + 1) >= 0); //! might be bad using gmres, add this to gmres inc!
1735 real declineV = ret(I4 + 1) / (u(I4 + 1) + 1e-6);
1736 real newu5 = u(I4 + 1) * std::exp(declineV);
1737 // ! refvalue:
1738 real muRef = settings.idealGasProperty.muGas;
1739 newu5 = std::max(1e-6, newu5);
1740 ret(I4 + 1) = newu5 - u(I4 + 1);
1741 }
1742 }
1743#endif
1744 if constexpr (Traits::has2EQ)
1745 {
1746 if (u(I4 + 1) + ret(I4 + 1) < 0)
1747 {
1748 // std::cout << "Fixing KE inc " << std::endl;
1749
1750 DNDS_assert(u(I4 + 1) >= 0); //! might be bad using gmres, add this to gmres inc!
1751 real declineV = ret(I4 + 1) / (u(I4 + 1) + 1e-6);
1752 real newu5 = u(I4 + 1) * std::exp(declineV);
1753 // ! refvalue:
1754 real muRef = settings.idealGasProperty.muGas;
1755 // newu5 = std::max(1e-10, newu5);
1756 ret(I4 + 1) = newu5 - u(I4 + 1);
1757 }
1758
1759 if (u(I4 + 2) + ret(I4 + 2) < 0)
1760 {
1761 // std::cout << "Fixing KE inc " << std::endl;
1762
1763 DNDS_assert(u(I4 + 2) >= 0); //! might be bad using gmres, add this to gmres inc!
1764 real declineV = ret(I4 + 2) / (u(I4 + 2) + 1e-6);
1765 real newu5 = u(I4 + 2) * std::exp(declineV);
1766 // ! refvalue:
1767 real muRef = settings.idealGasProperty.muGas;
1768 // newu5 = std::max(1e-10, newu5);
1769 ret(I4 + 2) = newu5 - u(I4 + 2);
1770 }
1771 }
1772
1773 return ret;
1774 }
1775
1776 /// @brief Apply CompressInc to every cell, modifying cxInc in-place.
1777 void FixIncrement(
1778 ArrayDOFV<nVarsFixed> &cx,
1779 ArrayDOFV<nVarsFixed> &cxInc, real alpha = 1.0)
1780 {
1781 for (index iCell = 0; iCell < cxInc.Size(); iCell++)
1782 cxInc[iCell] = this->CompressInc(cx[iCell], cxInc[iCell] * alpha);
1783 }
1784
1785 /**
1786 * @brief Add a positivity-compressed increment to the solution.
1787 *
1788 * For each cell, compresses the increment via CompressInc, then adds it to cx.
1789 * Reports the global minimum compression factor via MPI reduction.
1790 *
1791 * @param[in,out] cx Solution array (updated in-place).
1792 * @param[in] cxInc Increment array.
1793 * @param[in] alpha Scaling factor applied to the increment before compression.
1794 */
1795 void AddFixedIncrement(
1796 ArrayDOFV<nVarsFixed> &cx,
1797 ArrayDOFV<nVarsFixed> &cxInc, real alpha = 1.0)
1798 {
1799 DNDS_FV_EULEREVALUATOR_GET_FIXED_EIGEN_SEQS
1800 real alpha_fix_min = 1.0;
1801 for (index iCell = 0; iCell < cxInc.Size(); iCell++)
1802 {
1803 TU compressedInc = this->CompressInc(cx[iCell], cxInc[iCell] * alpha);
1804 real newAlpha = std::abs(compressedInc(0)) /
1805 (std::abs((cxInc[iCell] * alpha)(0)));
1806 if (std::abs((cxInc[iCell] * alpha)(0)) < verySmallReal)
1807 newAlpha = 1.; //! old inc could be zero, so compresion alpha is always 1
1808 alpha_fix_min = std::min(
1809 alpha_fix_min,
1810 newAlpha);
1811 // if (newAlpha < 1.0 - 1e-14)
1812 // std::cout << "KL\n"
1813 // << std::scientific << std::setprecision(5)
1814 // << this->CompressInc(cx[iCell], cxInc[iCell] * alpha).transpose() << "\n"
1815 // << cxInc[iCell].transpose() * alpha << std::endl;
1816 cx[iCell] += compressedInc;
1817 // wall fix in not needed
1818 // if (model == NS_2EQ || model == NS_2EQ_3D)
1819 // if (iCell < mesh->NumCell())
1820 // for (auto f : mesh->cell2face[iCell])
1821 // if (pBCHandler->GetTypeFromID(mesh->GetFaceZone(f)) == BCWall || pBCHandler->GetTypeFromID(mesh->GetFaceZone(f)) == BCWallIsothermal)
1822 // { // for SST or KOWilcox
1823 // TVec uNorm = vfv->GetFaceNorm(f, -1)(Seq012);
1824 // real vt = (cx[iCell](Seq123) - cx[iCell](Seq123).dot(uNorm) * uNorm).norm() / cx[iCell](0);
1825 // // cx[iCell](I4 + 1) = sqr(vt) * cx[iCell](0) * 1; // k = v_tang ^2 in sublayer, Wilcox book
1826 // // cx[iCell](I4 + 1) *= 0;
1827
1828 // real d1 = dWall[iCell].mean();
1829 // // cx[iCell](I4 + 1) = 0.; // superfix, actually works
1830 // // real d1 = dWall[iCell].minCoeff();
1831 // real pMean, asqrMean, Hmean;
1832 // real gamma = settings.idealGasProperty.gamma;
1833 // auto ULMeanXy = cx[iCell];
1834 // Gas::IdealGasThermal(ULMeanXy(I4), ULMeanXy(0), (ULMeanXy(Seq123) / ULMeanXy(0)).squaredNorm(),
1835 // gamma, pMean, asqrMean, Hmean);
1836 // // ! refvalue:
1837 // real muRef = settings.idealGasProperty.muGas;
1838 // real T = pMean / ((gamma - 1) / gamma * settings.idealGasProperty.CpGas * ULMeanXy(0));
1839 // real mufPhy1 = muEff(ULMeanXy, T);
1840
1841 // real rhoOmegaaaWall = mufPhy1 / sqr(d1) * 800;
1842 // // cx[iCell](I4 + 2) = rhoOmegaaaWall * 0.5; // this is bad
1843 // }
1844
1845 if (model == NS_2EQ || model == NS_2EQ_3D)
1846 { // for SST or KOWilcox
1847 if (settings.ransModel == RANSModel::RANS_KOSST ||
1848 settings.ransModel == RANSModel::RANS_KOWilcox)
1849 cx[iCell](I4 + 2) = std::max(cx[iCell](I4 + 2), settings.RANSBottomLimit * settings.farFieldStaticValue(I4 + 2));
1850 }
1851 }
1852 real alpha_fix_min_c = alpha_fix_min;
1853 MPI::Allreduce(&alpha_fix_min_c, &alpha_fix_min, 1, DNDS_MPI_REAL, MPI_MIN, cx.father->getMPI().comm);
1854 if (alpha_fix_min < 1.0)
1855 if (cx.father->getMPI().rank == 0)
1856 log() << TermColor::Magenta << "Increment fixed " << std::scientific << std::setprecision(5) << alpha_fix_min << TermColor::Reset << std::endl;
1857 }
1858
1859 // void AddFixedIncrement(
1860 // ArrayDOFV<nVarsFixed> &cx,
1861 // ArrayDOFV<nVarsFixed> &cxInc, real alpha = 1.0)
1862 // {
1863 // real alpha_fix_min = 1.0;
1864 // for (index iCell = 0; iCell < cxInc.Size(); iCell++)
1865 // {
1866 // TU compressedInc = this->CompressInc(cx[iCell], cxInc[iCell] * alpha);
1867 // real newAlpha = std::abs(compressedInc(0)) /
1868 // (std::abs((cxInc[iCell] * alpha)(0)));
1869 // if (std::abs((cxInc[iCell] * alpha)(0)) < verySmallReal)
1870 // newAlpha = 1.; //! old inc could be zero, so compresion alpha is always 1
1871 // alpha_fix_min = std::min(
1872 // alpha_fix_min,
1873 // newAlpha);
1874 // // if (newAlpha < 1.0 - 1e-14)
1875 // // std::cout << "KL\n"
1876 // // << std::scientific << std::setprecision(5)
1877 // // << this->CompressInc(cx[iCell], cxInc[iCell] * alpha).transpose() << "\n"
1878 // // << cxInc[iCell].transpose() * alpha << std::endl;
1879 // // cx[iCell] += compressedInc;
1880 // }
1881 // real alpha_fix_min_c = alpha_fix_min;
1882 // MPI::Allreduce(&alpha_fix_min_c, &alpha_fix_min, 1, DNDS_MPI_REAL, MPI_MIN, cx.father->getMPI().comm);
1883 // if (alpha_fix_min < 1.0)
1884 // if (cx.father->getMPI().rank == 0)
1885 // std::cout << "Increment fixed " << std::scientific << std::setprecision(5) << alpha_fix_min << std::endl;
1886
1887 // // for (index iCell = 0; iCell < cxInc.Size(); iCell++)
1888 // // cx[iCell] += alpha_fix_min * alpha * cxInc[iCell];
1889 // index nFixed = 0;
1890 // for (index iCell = 0; iCell < cxInc.Size(); iCell++)
1891 // {
1892 // TU compressedInc = this->CompressInc(cx[iCell], cxInc[iCell] * alpha);
1893 // real newAlpha = std::abs(compressedInc(0)) /
1894 // (std::abs((cxInc[iCell] * alpha)(0)));
1895 // if (std::abs((cxInc[iCell] * alpha)(0)) < verySmallReal)
1896 // newAlpha = 1.; //! old inc could be zero, so compresion alpha is always 1
1897
1898 // // if (newAlpha < 1.0 - 1e-14)
1899 // // std::cout << "KL\n"
1900 // // << std::scientific << std::setprecision(5)
1901 // // << this->CompressInc(cx[iCell], cxInc[iCell] * alpha).transpose() << "\n"
1902 // // << cxInc[iCell].transpose() * alpha << std::endl;
1903 // cx[iCell] += (newAlpha < 1 ? nFixed ++,alpha_fix_min : 1) * alpha * cxInc[iCell];
1904 // }
1905 // index nFixed_c = nFixed;
1906 // MPI::Allreduce(&nFixed_c, &nFixed, 1, DNDS_MPI_INDEX, MPI_SUM, cx.father->getMPI().comm);
1907 // if (alpha_fix_min < 1.0)
1908 // if (cx.father->getMPI().rank == 0)
1909 // std::cout << "Increment fixed number " << nFixed_c << std::endl;
1910 // }
1911
1912 /**
1913 * @brief Apply Laplacian smoothing to a residual field.
1914 *
1915 * Iteratively smooths r into rs using weighted neighbor averaging.
1916 * Used to stabilize central-difference residuals on coarse grids.
1917 *
1918 * @param[in] r Input residual.
1919 * @param[in,out] rs Smoothed residual (output; also used as work buffer).
1920 * @param[out] rtemp Temporary buffer (same size as r).
1921 * @param[in] nStep Number of smoothing passes (0 = use settings.nCentralSmoothStep).
1922 */
1923 void CentralSmoothResidual(ArrayDOFV<nVarsFixed> &r, ArrayDOFV<nVarsFixed> &rs, ArrayDOFV<nVarsFixed> &rtemp, int nStep = 0)
1924 {
1925 for (int iterS = 1; iterS <= (nStep > 0 ? nStep : settings.nCentralSmoothStep); iterS++)
1926 {
1927 real epsC = settings.centralSmoothEps;
1928#if defined(DNDS_DIST_MT_USE_OMP)
1929# pragma omp parallel for schedule(runtime)
1930#endif
1931 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
1932 {
1933 real div = 1.;
1934 TU vC = r[iCell];
1935 auto c2f = mesh->cell2face[iCell];
1936 for (int ic2f = 0; ic2f < c2f.size(); ic2f++)
1937 {
1938 index iFace = c2f[ic2f];
1939 index iCellOther = vfv->CellFaceOther(iCell, iFace);
1940 if (iCellOther != UnInitIndex)
1941 {
1942 div += epsC;
1943 vC += epsC * rs[iCellOther];
1944 }
1945 }
1946 rtemp[iCell] = vC / div;
1947 }
1948 rs = rtemp;
1949 rs.trans.startPersistentPull();
1950 rs.trans.waitPersistentPull();
1951 }
1952 }
1953
1954 /**
1955 * @brief Set initial conservative DOF values for all cells.
1956 *
1957 * Populates u with the far-field state (or a problem-specific initial condition)
1958 * based on the evaluator settings. Handles rotating-frame velocity transformations.
1959 *
1960 * @param[out] u Cell-centered DOF array to initialize.
1961 */
1962 void InitializeUDOF(ArrayDOFV<nVarsFixed> &u);
1963
1964 /**
1965 * @brief References to arrays needed by the output data picker.
1966 *
1967 * Groups the solution, reconstruction, and PP-limiter arrays for
1968 * use by InitializeOutputPicker to set up probe/output callbacks.
1969 */
1970 struct OutputOverlapDataRefs
1971 {
1972 ArrayDOFV<nVarsFixed> &u; ///< Cell-centered conservative DOFs.
1973 ArrayRECV<nVarsFixed> &uRec; ///< Reconstruction coefficients.
1974 ArrayDOFV<1> &betaPP; ///< PP reconstruction limiter beta.
1975 ArrayDOFV<1> &alphaPP; ///< PP RHS limiter alpha.
1976 };
1977
1978 /// @brief Initialize an OutputPicker with field callbacks for VTK/HDF5 output.
1979 void InitializeOutputPicker(OutputPicker &op, OutputOverlapDataRefs dataRefs);
1980 };
1981}
1982
1983#define DNDS_EulerEvaluator_INS_EXTERN(model, ext) \
1984 namespace DNDS::Euler \
1985 { \
1986 ext template void EulerEvaluator<model>::LUSGSMatrixInit( \
1987 JacobianDiagBlock<nVarsFixed> &JDiag, \
1988 JacobianDiagBlock<nVarsFixed> &JSource, \
1989 ArrayDOFV<1> &dTau, real dt, real alphaDiag, \
1990 ArrayDOFV<nVarsFixed> &u, \
1991 ArrayRECV<nVarsFixed> &uRec, \
1992 int jacobianCode, \
1993 real t); \
1994 \
1995 ext template void EulerEvaluator<model>::LUSGSMatrixVec( \
1996 real alphaDiag, \
1997 real t, \
1998 ArrayDOFV<nVarsFixed> &u, \
1999 ArrayDOFV<nVarsFixed> &uInc, \
2000 JacobianDiagBlock<nVarsFixed> &JDiag, \
2001 ArrayDOFV<nVarsFixed> &AuInc); \
2002 \
2003 ext template void EulerEvaluator<model>::LUSGSMatrixToJacobianLU( \
2004 real alphaDiag, \
2005 real t, \
2006 ArrayDOFV<nVarsFixed> &u, \
2007 JacobianDiagBlock<nVarsFixed> &JDiag, \
2008 JacobianLocalLU<nVarsFixed> &jacLU); \
2009 \
2010 ext template void EulerEvaluator<model>::UpdateLUSGSForward( \
2011 real alphaDiag, \
2012 real t, \
2013 ArrayDOFV<nVarsFixed> &rhs, \
2014 ArrayDOFV<nVarsFixed> &u, \
2015 ArrayDOFV<nVarsFixed> &uInc, \
2016 JacobianDiagBlock<nVarsFixed> &JDiag, \
2017 ArrayDOFV<nVarsFixed> &uIncNew); \
2018 \
2019 ext template void EulerEvaluator<model>::UpdateLUSGSBackward( \
2020 real alphaDiag, \
2021 real t, \
2022 ArrayDOFV<nVarsFixed> &rhs, \
2023 ArrayDOFV<nVarsFixed> &u, \
2024 ArrayDOFV<nVarsFixed> &uInc, \
2025 JacobianDiagBlock<nVarsFixed> &JDiag, \
2026 ArrayDOFV<nVarsFixed> &uIncNew); \
2027 \
2028 ext template void EulerEvaluator<model>::UpdateSGS( \
2029 real alphaDiag, \
2030 real t, \
2031 ArrayDOFV<nVarsFixed> &rhs, \
2032 ArrayDOFV<nVarsFixed> &u, \
2033 ArrayDOFV<nVarsFixed> &uInc, \
2034 ArrayDOFV<nVarsFixed> &uIncNew, \
2035 JacobianDiagBlock<nVarsFixed> &JDiag, \
2036 bool forward, bool gsUpdate, TU &sumInc, \
2037 bool uIncIsZero); \
2038 ext template void EulerEvaluator<model>::UpdateSGSWithRec( \
2039 real alphaDiag, \
2040 real t, \
2041 ArrayDOFV<nVarsFixed> &rhs, \
2042 ArrayDOFV<nVarsFixed> &u, \
2043 ArrayRECV<nVarsFixed> &uRec, \
2044 ArrayDOFV<nVarsFixed> &uInc, \
2045 ArrayRECV<nVarsFixed> &uRecInc, \
2046 JacobianDiagBlock<nVarsFixed> &JDiag, \
2047 bool forward, TU &sumInc); \
2048 \
2049 ext template void EulerEvaluator<model>::LUSGSMatrixSolveJacobianLU( \
2050 real alphaDiag, \
2051 real t, \
2052 ArrayDOFV<nVarsFixed> &rhs, \
2053 ArrayDOFV<nVarsFixed> &u, \
2054 ArrayDOFV<nVarsFixed> &uInc, \
2055 ArrayDOFV<nVarsFixed> &uIncNew, \
2056 ArrayDOFV<nVarsFixed> &bBuf, \
2057 JacobianDiagBlock<nVarsFixed> &JDiag, \
2058 JacobianLocalLU<nVarsFixed> &jacLU, \
2059 bool uIncIsZero, \
2060 TU &sumInc); \
2061 \
2062 ext template void EulerEvaluator<model>::InitializeUDOF(ArrayDOFV<nVarsFixed> &u); \
2063 \
2064 ext template void EulerEvaluator<model>::FixUMaxFilter( \
2065 ArrayDOFV<nVarsFixed> &u); \
2066 \
2067 ext template void EulerEvaluator<model>::TimeAverageAddition( \
2068 ArrayDOFV<nVarsFixed> &w, ArrayDOFV<nVarsFixed> &wAveraged, real dt, real &tCur); \
2069 ext template void EulerEvaluator<model>::MeanValueCons2Prim( \
2070 ArrayDOFV<nVarsFixed> &u, ArrayDOFV<nVarsFixed> &w); \
2071 ext template void EulerEvaluator<model>::MeanValuePrim2Cons( \
2072 ArrayDOFV<nVarsFixed> &w, ArrayDOFV<nVarsFixed> &u); \
2073 \
2074 ext template void EulerEvaluator<model>::EvaluateNorm( \
2075 Eigen::Vector<real, -1> &res, ArrayDOFV<nVarsFixed> &rhs, index P, bool volWise, bool average); \
2076 \
2077 ext template void EulerEvaluator<model>::EvaluateRecNorm( \
2078 Eigen::Vector<real, -1> &res, \
2079 ArrayDOFV<nVarsFixed> &u, \
2080 ArrayRECV<nVarsFixed> &uRec, \
2081 index P, \
2082 bool compare, \
2083 const tFCompareField &FCompareField, \
2084 const tFCompareFieldWeight &FCompareFieldWeight, \
2085 real t); \
2086 \
2087 ext template void EulerEvaluator<model>::LimiterUGrad( \
2088 ArrayDOFV<nVarsFixed> &u, ArrayGRADV<nVarsFixed, gDim> &uGrad, ArrayGRADV<nVarsFixed, gDim> &uGradNew, \
2089 uint64_t flags); \
2090 \
2091 ext template void EulerEvaluator<model>::EvaluateURecBeta( \
2092 ArrayDOFV<nVarsFixed> &u, \
2093 ArrayRECV<nVarsFixed> &uRec, \
2094 ArrayDOFV<1> &uRecBeta, index &nLim, real &betaMin, int flag); \
2095 \
2096 ext template bool EulerEvaluator<model>::AssertMeanValuePP( \
2097 ArrayDOFV<nVarsFixed> &u, bool panic); \
2098 \
2099 ext template void EulerEvaluator<model>::EvaluateCellRHSAlpha( \
2100 ArrayDOFV<nVarsFixed> &u, \
2101 ArrayRECV<nVarsFixed> &uRec, \
2102 ArrayDOFV<1> &uRecBeta, \
2103 ArrayDOFV<nVarsFixed> &rhs, \
2104 ArrayDOFV<1> &cellRHSAlpha, index &nLim, real &alphaMin, real relax, \
2105 int compress, \
2106 int flag); \
2107 \
2108 ext template void EulerEvaluator<model>::EvaluateCellRHSAlphaExpansion( \
2109 ArrayDOFV<nVarsFixed> &u, \
2110 ArrayRECV<nVarsFixed> &uRec, \
2111 ArrayDOFV<1> &uRecBeta, \
2112 ArrayDOFV<nVarsFixed> &res, \
2113 ArrayDOFV<1> &cellRHSAlpha, index &nLim, real alphaMin); \
2114 ext template void EulerEvaluator<model>::MinSmoothDTau( \
2115 ArrayDOFV<1> &dTau, ArrayDOFV<1> &dTauNew); \
2116 ext template void EulerEvaluator<model>::updateBCProfiles(ArrayDOFV<nVarsFixed> &u, ArrayRECV<nVarsFixed> &uRec); \
2117 ext template void EulerEvaluator<model>::updateBCProfilesPressureRadialEq(); \
2118 }
2119
2120DNDS_EulerEvaluator_INS_EXTERN(NS, extern);
2121DNDS_EulerEvaluator_INS_EXTERN(NS_2D, extern);
2122DNDS_EulerEvaluator_INS_EXTERN(NS_SA, extern);
2123DNDS_EulerEvaluator_INS_EXTERN(NS_2EQ, extern);
2124DNDS_EulerEvaluator_INS_EXTERN(NS_3D, extern);
2125DNDS_EulerEvaluator_INS_EXTERN(NS_SA_3D, extern);
2126DNDS_EulerEvaluator_INS_EXTERN(NS_2EQ_3D, extern);
2127
2128#define DNDS_EulerEvaluator_EvaluateDt_INS_EXTERN(model, ext) \
2129 namespace DNDS::Euler \
2130 { \
2131 ext template void EulerEvaluator<model>::GetWallDist(); \
2132 ext template void EulerEvaluator<model>::GetWallDist_CollectTriangles( \
2133 bool, std::vector<Eigen::Matrix<real, 3, 3>> &); \
2134 ext template void EulerEvaluator<model>::GetWallDist_AABB(); \
2135 ext template void EulerEvaluator<model>::GetWallDist_BatchedAABB(); \
2136 ext template void EulerEvaluator<model>::GetWallDist_Poisson(); \
2137 ext template void EulerEvaluator<model>::GetWallDist_ComputeFaceDistances(); \
2138 ext template void EulerEvaluator<model>::EvaluateDt( \
2139 ArrayDOFV<1> &dt, \
2140 ArrayDOFV<nVarsFixed> &u, \
2141 ArrayRECV<nVarsFixed> &uRec, \
2142 real CFL, real &dtMinall, real MaxDt, \
2143 bool UseLocaldt, \
2144 real t, \
2145 uint64_t flags); \
2146 ext template void \
2147 EulerEvaluator<model>::fluxFace( \
2148 const TU_Batch &ULxy, \
2149 const TU_Batch &URxy, \
2150 const TU &ULMeanXy, \
2151 const TU &URMeanXy, \
2152 const TDiffU_Batch &DiffUxy, \
2153 const TDiffU_Batch &DiffUxyPrim, \
2154 const TVec_Batch &unitNorm, \
2155 const TVec_Batch &vg, \
2156 const TVec &unitNormC, \
2157 const TVec &vgC, \
2158 TU_Batch &FLfix, \
2159 TU_Batch &FRfix, \
2160 TU_Batch &finc, \
2161 TReal_Batch &lam0V, TReal_Batch &lam123V, TReal_Batch &lam4V, \
2162 Geom::t_index btype, \
2163 typename Gas::RiemannSolverType rsType, \
2164 index iFace, bool ignoreVis); \
2165 ext template \
2166 typename EulerEvaluator<model>::TU \
2167 EulerEvaluator<model>::source( \
2168 const TU &UMeanXy, \
2169 const TDiffU &DiffUxy, \
2170 const Geom::tPoint &pPhy, \
2171 TJacobianU &jacobian, \
2172 index iCell, \
2173 index ig, \
2174 int Mode); \
2175 ext template \
2176 typename EulerEvaluator<model>::TU \
2177 EulerEvaluator<model>::generateBoundaryValue( \
2178 TU &ULxy, \
2179 const TU &ULMeanXy, \
2180 index iCell, index iFace, int iG, \
2181 const TVec &uNorm, \
2182 const TMat &normBase, \
2183 const Geom::tPoint &pPhysics, \
2184 real t, \
2185 Geom::t_index btype, \
2186 bool fixUL, \
2187 int geomMode, int linMode); \
2188 ext template typename EulerEvaluator<model>::TU EulerEvaluator<model>::generateBV_FarField( \
2189 TU &, const TU &, index, index, int, const TVec &, const TMat &, \
2190 const Geom::tPoint &, real, Geom::t_index, bool, int); \
2191 ext template typename EulerEvaluator<model>::TU EulerEvaluator<model>::generateBV_SpecialFar( \
2192 TU &, const TU &, index, index, int, const TVec &, const TMat &, \
2193 const Geom::tPoint &, real, Geom::t_index); \
2194 ext template typename EulerEvaluator<model>::TU EulerEvaluator<model>::generateBV_InviscidWall( \
2195 TU &, const TU &, index, index, int, const TVec &, const TMat &, \
2196 const Geom::tPoint &, real, Geom::t_index); \
2197 ext template typename EulerEvaluator<model>::TU EulerEvaluator<model>::generateBV_ViscousWall( \
2198 TU &, const TU &, index, index, int, const TVec &, const TMat &, \
2199 const Geom::tPoint &, real, Geom::t_index, bool, int); \
2200 ext template typename EulerEvaluator<model>::TU EulerEvaluator<model>::generateBV_Outflow( \
2201 TU &, const TU &, index, index, int, const TVec &, const TMat &, \
2202 const Geom::tPoint &, real, Geom::t_index); \
2203 ext template typename EulerEvaluator<model>::TU EulerEvaluator<model>::generateBV_Inflow( \
2204 TU &, const TU &, index, index, int, const TVec &, const TMat &, \
2205 const Geom::tPoint &, real, Geom::t_index); \
2206 ext template typename EulerEvaluator<model>::TU EulerEvaluator<model>::generateBV_TotalConditionInflow( \
2207 TU &, const TU &, index, index, int, const TVec &, const TMat &, \
2208 const Geom::tPoint &, real, Geom::t_index); \
2209 ext template void EulerEvaluator<model>::InitializeOutputPicker(OutputPicker &op, OutputOverlapDataRefs dataRefs); \
2210 }
2211
2212DNDS_EulerEvaluator_EvaluateDt_INS_EXTERN(NS, extern);
2213DNDS_EulerEvaluator_EvaluateDt_INS_EXTERN(NS_2D, extern);
2214DNDS_EulerEvaluator_EvaluateDt_INS_EXTERN(NS_SA, extern);
2215DNDS_EulerEvaluator_EvaluateDt_INS_EXTERN(NS_2EQ, extern);
2216DNDS_EulerEvaluator_EvaluateDt_INS_EXTERN(NS_3D, extern);
2217DNDS_EulerEvaluator_EvaluateDt_INS_EXTERN(NS_SA_3D, extern);
2218DNDS_EulerEvaluator_EvaluateDt_INS_EXTERN(NS_2EQ_3D, extern);
2219
2220#define DNDS_EulerEvaluator_EvaluateRHS_INS_EXTERN(model, ext) \
2221 namespace DNDS::Euler \
2222 { \
2223 ext template void EulerEvaluator<model>::EvaluateRHS( \
2224 ArrayDOFV<nVarsFixed> &rhs, \
2225 JacobianDiagBlock<nVarsFixed> &JSource, \
2226 ArrayDOFV<nVarsFixed> &u, \
2227 ArrayRECV<nVarsFixed> &uRecUnlim, \
2228 ArrayRECV<nVarsFixed> &uRec, \
2229 ArrayDOFV<1> &uRecBeta, \
2230 ArrayDOFV<1> &cellRHSAlpha, \
2231 bool onlyOnHalfAlpha, \
2232 real t, \
2233 uint64_t flags); \
2234 }
2235
2236DNDS_EulerEvaluator_EvaluateRHS_INS_EXTERN(NS, extern);
2237DNDS_EulerEvaluator_EvaluateRHS_INS_EXTERN(NS_2D, extern);
2238DNDS_EulerEvaluator_EvaluateRHS_INS_EXTERN(NS_SA, extern);
2239DNDS_EulerEvaluator_EvaluateRHS_INS_EXTERN(NS_2EQ, extern);
2240DNDS_EulerEvaluator_EvaluateRHS_INS_EXTERN(NS_3D, extern);
2241DNDS_EulerEvaluator_EvaluateRHS_INS_EXTERN(NS_SA_3D, extern);
2242DNDS_EulerEvaluator_EvaluateRHS_INS_EXTERN(NS_2EQ_3D, extern);
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
EulerBC.hpp
Boundary condition types, handlers, integration recorders, and 1-D profile utilities for the compress...
EulerEvaluatorSettings.hpp
Complete solver configuration for the Euler/Navier-Stokes evaluator.
DNDS_EulerEvaluator_INS_EXTERN
#define DNDS_EulerEvaluator_INS_EXTERN(model, ext)
Definition EulerEvaluator.hpp:1983
DNDS_EulerEvaluator_EvaluateDt_INS_EXTERN
#define DNDS_EulerEvaluator_EvaluateDt_INS_EXTERN(model, ext)
Definition EulerEvaluator.hpp:2128
DNDS_EulerEvaluator_EvaluateRHS_INS_EXTERN
#define DNDS_EulerEvaluator_EvaluateRHS_INS_EXTERN(model, ext)
Definition EulerEvaluator.hpp:2220
EulerJacobian.hpp
Jacobian storage and factorization structures for implicit time stepping.
Euler.hpp
Core type definitions, enumerations, and distributed array wrappers for compressible Navier-Stokes / ...
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
Gas.hpp
Ideal-gas Riemann solvers, flux functions, and thermodynamic utilities for the compressible Euler / N...
JsonUtil.hpp
JSON-to-Eigen conversion utilities and nlohmann_json helper macros.
RANS_ke.hpp
RANS two-equation turbulence model implementations for the DNDSR CFD solver.
SerializerBase.hpp
Base types and abstract interface for array serialization.
DNDS::CFV::VariationalReconstruction
The VR class that provides any information needed in high-order CFV.
Definition VariationalReconstruction.hpp:110
DNDS::Euler::ArrayDOFV
MPI-parallel distributed array of per-cell Degrees-of-Freedom (conservative variable vectors).
Definition Euler.hpp:93
DNDS::Euler::ArrayGRADV
MPI-parallel distributed array of per-cell gradient matrices.
Definition Euler.hpp:643
DNDS::Euler::ArrayRECV
MPI-parallel distributed array of per-cell reconstruction coefficient matrices.
Definition Euler.hpp:401
DNDS::Euler::EulerEvaluator
Core finite-volume evaluator for compressible Navier-Stokes / Euler equations.
Definition EulerEvaluator.hpp:78
DNDS::Euler::EulerEvaluator::TJacobianU
Eigen::MatrixFMTSafe< real, nVarsFixed, nVarsFixed > TJacobianU
Jacobian matrix (nVars x nVars) of the flux w.r.t. conserved variables.
Definition EulerEvaluator.hpp:97
DNDS::Euler::EulerEvaluator::TU_Batch
Eigen::MatrixFMTSafe< real, nVarsFixed, Eigen::Dynamic, Eigen::ColMajor, nVarsFixed, MaxBatch > TU_Batch
Batch of conservative variable vectors.
Definition EulerEvaluator.hpp:95
DNDS::Euler::EulerEvaluator::TMat
Eigen::MatrixFMTSafe< real, dim, dim > TMat
Spatial matrix (dim x dim).
Definition EulerEvaluator.hpp:92
DNDS::Euler::EulerEvaluator::TMat_Batch
Eigen::MatrixFMTSafe< real, dim, Eigen::Dynamic, Eigen::ColMajor, dim, MaxBatchMult(3)> TMat_Batch
Batch of spatial matrices.
Definition EulerEvaluator.hpp:93
DNDS::Euler::EulerEvaluator::nFaceReducedOrder
index nFaceReducedOrder
Count of faces where reconstruction order was reduced.
Definition EulerEvaluator.hpp:216
DNDS::Euler::EulerEvaluator::RHS_Ignore_Viscosity
static const uint64_t RHS_Ignore_Viscosity
Skip viscous flux contribution.
Definition EulerEvaluator.hpp:371
DNDS::Euler::EulerEvaluator::TDiffUTransposed
Eigen::MatrixFMTSafe< real, nVarsFixed, dim > TDiffUTransposed
Transposed gradient (nVars x dim).
Definition EulerEvaluator.hpp:100
DNDS::Euler::EulerEvaluator::RHS_Recover_IncFScale
static const uint64_t RHS_Recover_IncFScale
Recover incremental face scaling.
Definition EulerEvaluator.hpp:378
DNDS::Euler::EulerEvaluator::setPassiveDiscardSource
void setPassiveDiscardSource(bool n)
Enable or disable discarding of passive-scalar source terms.
Definition EulerEvaluator.hpp:166
DNDS::Euler::EulerEvaluator::TpBCHandler
ssp< BoundaryHandler< model > > TpBCHandler
Shared pointer to the boundary condition handler.
Definition EulerEvaluator.hpp:107
DNDS::Euler::EulerEvaluator::InitializeFV
static void InitializeFV(ssp< Geom::UnstructuredMesh > &mesh, TpVFV vfv, TpBCHandler pBCHandler)
Initialize the finite-volume infrastructure on the mesh.
Definition EulerEvaluator.hpp:123
DNDS::Euler::EulerEvaluator::MinSmoothDTau
void MinSmoothDTau(ArrayDOFV< 1 > &dTau, ArrayDOFV< 1 > &dTauNew)
Smooth the local time step across neighboring cells.
Definition EulerEvaluator.hxx:2158
DNDS::Euler::EulerEvaluator::LUSGSMatrixInit
void LUSGSMatrixInit(JacobianDiagBlock< nVarsFixed > &JDiag, JacobianDiagBlock< nVarsFixed > &JSource, ArrayDOFV< 1 > &dTau, real dt, real alphaDiag, ArrayDOFV< nVarsFixed > &u, ArrayRECV< nVarsFixed > &uRec, int jacobianCode, real t)
Assemble the diagonal blocks of the implicit Jacobian for LU-SGS / SGS.
DNDS::Euler::EulerEvaluator::DT_Dont_update_lambda01234
static const uint64_t DT_Dont_update_lambda01234
Skip recomputation of per-face eigenvalues lambda0/123/4.
Definition EulerEvaluator.hpp:340
DNDS::Euler::EulerEvaluator::EvaluateRecNorm
void EvaluateRecNorm(Eigen::Vector< real, -1 > &res, ArrayDOFV< nVarsFixed > &u, ArrayRECV< nVarsFixed > &uRec, index P=1, bool compare=false, const tFCompareField &FCompareField=[](const Geom::tPoint &p, real t) { return TU::Zero();}, const tFCompareFieldWeight &FCompareFieldWeight=[](const Geom::tPoint &p, real t) { return 1.0;}, real t=0)
Compute the reconstruction error norm (optionally against an analytical field).
Definition EulerEvaluator.hxx:1474
DNDS::Euler::EulerEvaluator::GetWallDist
void GetWallDist()
Compute wall distance for all cells (dispatcher).
DNDS::Euler::EulerEvaluator::deltaLambdaCell
ArrayDOFV< 1 > deltaLambdaCell
Per-cell accumulated spectral radius difference.
Definition EulerEvaluator.hpp:191
DNDS::Euler::EulerEvaluator::TDiffU
Eigen::MatrixFMTSafe< real, dim, nVarsFixed > TDiffU
Gradient of conserved variables (dim x nVars).
Definition EulerEvaluator.hpp:98
DNDS::Euler::EulerEvaluator::RHS_Dont_Record_Bud_Flux
static const uint64_t RHS_Dont_Record_Bud_Flux
Skip recording per-boundary flux.
Definition EulerEvaluator.hpp:373
DNDS::Euler::EulerEvaluator::GetFaceVGridFromCell
TVec GetFaceVGridFromCell(index iFace, index iCell, int if2c, index iG)
Get the grid velocity at a face quadrature point from a specific cell's perspective.
Definition EulerEvaluator.hpp:1235
DNDS::Euler::EulerEvaluator::dWallFace
std::vector< real > dWallFace
Per-face wall distance (interpolated from cell values).
Definition EulerEvaluator.hpp:198
DNDS::Euler::EulerEvaluator::AssertMeanValuePP
bool AssertMeanValuePP(ArrayDOFV< nVarsFixed > &u, bool panic)
Assert that all cell-mean values are physically realizable.
Definition EulerEvaluator.hxx:1866
DNDS::Euler::EulerEvaluator::BndIntegrationLogWriteLine
void BndIntegrationLogWriteLine(const std::string &name, index step, index stage, index iter)
Append a line to the per-BC boundary integration CSV log files.
Definition EulerEvaluator.hpp:1466
DNDS::Euler::EulerEvaluator::CLDriverGetIntegrationUpdate
std::tuple< real, real, real > CLDriverGetIntegrationUpdate(index iter)
Query the CL driver for current lift/drag coefficients and update AoA.
Definition EulerEvaluator.hpp:1502
DNDS::Euler::EulerEvaluator::RHS_Direct_2nd_Rec
static const uint64_t RHS_Direct_2nd_Rec
Use direct 2nd-order reconstruction.
Definition EulerEvaluator.hpp:374
DNDS::Euler::EulerEvaluator::muEff
real muEff(const TU &U, real T)
Compute effective molecular viscosity using the configured viscosity model.
Definition EulerEvaluator.hpp:751
DNDS::Euler::EulerEvaluator::TReal_Batch
Eigen::MatrixFMTSafe< real, 1, Eigen::Dynamic, Eigen::RowMajor, 1, MaxBatch > TReal_Batch
Batch of scalar values.
Definition EulerEvaluator.hpp:96
DNDS::Euler::EulerEvaluator::UpdateLUSGSBackward
void UpdateLUSGSBackward(real alphaDiag, real t, ArrayDOFV< nVarsFixed > &rhs, ArrayDOFV< nVarsFixed > &u, ArrayDOFV< nVarsFixed > &uInc, JacobianDiagBlock< nVarsFixed > &JDiag, ArrayDOFV< nVarsFixed > &uIncNew)
Deprecated: use UpdateSGS with uIncIsZero=true instead. Kept for backward compatibility....
DNDS::Euler::EulerEvaluator::fluxBndForceT
std::vector< TVec > fluxBndForceT
Per-boundary-face tangential force.
Definition EulerEvaluator.hpp:215
DNDS::Euler::EulerEvaluator::RHS_Direct_2nd_Rec_already_have_uGradBufNoLim
static const uint64_t RHS_Direct_2nd_Rec_already_have_uGradBufNoLim
uGradBufNoLim is already computed.
Definition EulerEvaluator.hpp:377
DNDS::Euler::EulerEvaluator::EvaluateURecBeta_DEFAULT
static const int EvaluateURecBeta_DEFAULT
Default: evaluate beta without compression.
Definition EulerEvaluator.hpp:659
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::nVarsFixed
static const int nVarsFixed
Number of conserved variables (compile-time fixed).
Definition EulerEvaluator.hpp:81
DNDS::Euler::EulerEvaluator::lambdaFace123
std::vector< real > lambdaFace123
Per-face eigenvalue |u·n + a| (acoustic wave).
Definition EulerEvaluator.hpp:188
DNDS::Euler::EulerEvaluator::lambdaFace0
std::vector< real > lambdaFace0
Per-face eigenvalue |u·n| (contact wave).
Definition EulerEvaluator.hpp:187
DNDS::Euler::EulerEvaluator::CompressInc
TU CompressInc(const TU &u, const TU &uInc)
Compress a solution increment to maintain positivity.
Definition EulerEvaluator.hpp:1645
DNDS::Euler::EulerEvaluator::vfv
ssp< CFV::VariationalReconstruction< gDim > > vfv
Definition EulerEvaluator.hpp:178
DNDS::Euler::EulerEvaluator::deltaLambdaFace
std::vector< real > deltaLambdaFace
Per-face spectral radius difference for implicit diagonal.
Definition EulerEvaluator.hpp:190
DNDS::Euler::EulerEvaluator::InitializeUDOF
void InitializeUDOF(ArrayDOFV< nVarsFixed > &u)
Set initial conservative DOF values for all cells.
DNDS::Euler::EulerEvaluator::TU
Eigen::VectorFMTSafe< real, nVarsFixed > TU
Conservative variable vector (nVarsFixed components).
Definition EulerEvaluator.hpp:94
DNDS::Euler::EulerEvaluator::uGradBuf
ArrayGRADV< nVarsFixed, gDim > uGradBuf
Definition EulerEvaluator.hpp:211
DNDS::Euler::EulerEvaluator::EvaluateURecBeta_COMPRESS_TO_MEAN
static const int EvaluateURecBeta_COMPRESS_TO_MEAN
Compress reconstruction toward cell mean to enforce positivity.
Definition EulerEvaluator.hpp:660
DNDS::Euler::EulerEvaluator::LIMITER_UGRAD_No_Flags
static const uint64_t LIMITER_UGRAD_No_Flags
Default limiter flags.
Definition EulerEvaluator.hpp:641
DNDS::Euler::EulerEvaluator::GetFaceVGrid
TVec GetFaceVGrid(index iFace, index iG)
Get the grid velocity at a face quadrature point (rotating frame).
Definition EulerEvaluator.hpp:1208
DNDS::Euler::EulerEvaluator::cLDriverBndIDs
std::set< Geom::t_index > cLDriverBndIDs
Boundary IDs driven by the CL (lift-coefficient) driver.
Definition EulerEvaluator.hpp:206
DNDS::Euler::EulerEvaluator::anchorRecorders
std::map< Geom::t_index, AnchorPointRecorder< nVarsFixed > > anchorRecorders
Per-BC anchor point recorders for profile extraction.
Definition EulerEvaluator.hpp:201
DNDS::Euler::EulerEvaluator::profileRecorders
std::map< Geom::t_index, OneDimProfile< nVarsFixed > > profileRecorders
Per-BC 1-D profile recorders.
Definition EulerEvaluator.hpp:202
DNDS::Euler::EulerEvaluator::GetAxisSymmetric
int GetAxisSymmetric()
Return the axisymmetric mode flag.
Definition EulerEvaluator.hpp:173
DNDS::Euler::EulerEvaluator::CompressRecPart
TU CompressRecPart(const TU &umean, const TU &uRecInc, bool &compressed)
Compress a reconstruction increment to maintain positivity.
Definition EulerEvaluator.hpp:1541
DNDS::Euler::EulerEvaluator::FixUMaxFilter
void FixUMaxFilter(ArrayDOFV< nVarsFixed > &u)
Clip extreme conserved-variable values to prevent overflow.
Definition EulerEvaluator.hxx:1318
DNDS::Euler::EulerEvaluator::TVec_Batch
Eigen::MatrixFMTSafe< real, dim, Eigen::Dynamic, Eigen::ColMajor, dim, MaxBatch > TVec_Batch
Batch of spatial vectors.
Definition EulerEvaluator.hpp:91
DNDS::Euler::EulerEvaluator::FixIncrement
void FixIncrement(ArrayDOFV< nVarsFixed > &cx, ArrayDOFV< nVarsFixed > &cxInc, real alpha=1.0)
Apply CompressInc to every cell, modifying cxInc in-place.
Definition EulerEvaluator.hpp:1777
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::Euler::EulerEvaluator::TVFV
CFV::VariationalReconstruction< gDim > TVFV
Variational reconstruction type for this geometric dimension.
Definition EulerEvaluator.hpp:105
DNDS::Euler::EulerEvaluator::I4
static const auto I4
Index of the energy equation (= dim+1).
Definition EulerEvaluator.hpp:84
DNDS::Euler::EulerEvaluator::settings
EulerEvaluatorSettings< model > settings
Physics and numerics settings for this evaluator.
Definition EulerEvaluator.hpp:220
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::lambdaFaceC
std::vector< real > lambdaFaceC
Per-face convective spectral radius.
Definition EulerEvaluator.hpp:185
DNDS::Euler::EulerEvaluator::UpdateSGS
void UpdateSGS(real alphaDiag, real t, ArrayDOFV< nVarsFixed > &rhs, ArrayDOFV< nVarsFixed > &u, ArrayDOFV< nVarsFixed > &uInc, ArrayDOFV< nVarsFixed > &uIncNew, JacobianDiagBlock< nVarsFixed > &JDiag, bool forward, bool gsUpdate, TU &sumInc, bool uIncIsZero=false)
Symmetric Gauss-Seidel update for the implicit linear system.
Definition EulerEvaluator.hxx:550
DNDS::Euler::EulerEvaluator::EvaluateURecBeta
void EvaluateURecBeta(ArrayDOFV< nVarsFixed > &u, ArrayRECV< nVarsFixed > &uRec, ArrayDOFV< 1 > &uRecBeta, index &nLim, real &betaMin, int flag)
Evaluate the positivity-preserving reconstruction limiter (beta).
DNDS::Euler::EulerEvaluator::RHS_Direct_2nd_Rec_1st_Conv
static const uint64_t RHS_Direct_2nd_Rec_1st_Conv
2nd-order rec with 1st-order convection.
Definition EulerEvaluator.hpp:375
DNDS::Euler::EulerEvaluator::MaxBatchMult
static constexpr int MaxBatchMult(int n)
Compute the maximum batch-multiplied column count for batched Eigen matrices.
Definition EulerEvaluator.hpp:88
DNDS::Euler::EulerEvaluator::TpVFV
ssp< CFV::VariationalReconstruction< gDim > > TpVFV
Shared pointer to the variational reconstruction object.
Definition EulerEvaluator.hpp:106
DNDS::Euler::EulerEvaluator::kAv
int kAv
Artificial viscosity polynomial order (maxOrder + 1).
Definition EulerEvaluator.hpp:180
DNDS::Euler::EulerEvaluator::fluxWallSum
Eigen::Vector< real, -1 > fluxWallSum
Accumulated wall flux integral (force on wall boundaries).
Definition EulerEvaluator.hpp:213
DNDS::Euler::EulerEvaluator::RHS_Direct_2nd_Rec_use_limiter
static const uint64_t RHS_Direct_2nd_Rec_use_limiter
Apply limiter when using direct 2nd rec.
Definition EulerEvaluator.hpp:376
DNDS::Euler::EulerEvaluator::LUSGSMatrixToJacobianLU
void LUSGSMatrixToJacobianLU(real alphaDiag, real t, ArrayDOFV< nVarsFixed > &u, JacobianDiagBlock< nVarsFixed > &JDiag, JacobianLocalLU< nVarsFixed > &jacLU)
Build the local LU factorization of the Jacobian for direct solve.
DNDS::Euler::EulerEvaluator::pCLDriver
std::unique_ptr< CLDriver > pCLDriver
Lift-coefficient driver for AoA adaptation (null if unused).
Definition EulerEvaluator.hpp:207
DNDS::Euler::EulerEvaluator::TDof
ArrayDOFV< nVarsFixed > TDof
Cell-centered DOF array (mean values).
Definition EulerEvaluator.hpp:101
DNDS::Euler::EulerEvaluator::EvaluateCellRHSAlpha_DEFAULT
static const int EvaluateCellRHSAlpha_DEFAULT
Default alpha evaluation mode.
Definition EulerEvaluator.hpp:691
DNDS::Euler::EulerEvaluator::MeanValuePrim2Cons
void MeanValuePrim2Cons(ArrayDOFV< nVarsFixed > &w, ArrayDOFV< nVarsFixed > &u)
Convert cell-mean primitive variables to conservative variables.
Definition EulerEvaluator.hxx:1368
DNDS::Euler::EulerEvaluator::TransformURotatingFrame
void TransformURotatingFrame(TU &U, const Geom::tPoint &pPhysics, int direction)
Transform full conservative state (momentum + total energy) between frames (relative velocity formula...
Definition EulerEvaluator.hpp:1255
DNDS::Euler::EulerEvaluator::AddFixedIncrement
void AddFixedIncrement(ArrayDOFV< nVarsFixed > &cx, ArrayDOFV< nVarsFixed > &cxInc, real alpha=1.0)
Add a positivity-compressed increment to the solution.
Definition EulerEvaluator.hpp:1795
DNDS::Euler::EulerEvaluator::gradUFix
std::vector< TDiffU > gradUFix
Green-Gauss gradient correction buffer for source term stabilization.
Definition EulerEvaluator.hpp:194
DNDS::Euler::EulerEvaluator::LIMITER_UGRAD_Disable_Shock_Limiter
static const uint64_t LIMITER_UGRAD_Disable_Shock_Limiter
Disable shock-detecting component of the limiter.
Definition EulerEvaluator.hpp:642
DNDS::Euler::EulerEvaluator::PrintBCProfiles
void PrintBCProfiles(const std::string &name, ArrayDOFV< nVarsFixed > &u, ArrayRECV< nVarsFixed > &uRec)
Write boundary profile data to CSV files (rank 0 only).
Definition EulerEvaluator.hpp:1430
DNDS::Euler::EulerEvaluator::uGradBufNoLim
ArrayGRADV< nVarsFixed, gDim > uGradBufNoLim
Gradient buffers: limited and unlimited.
Definition EulerEvaluator.hpp:211
DNDS::Euler::EulerEvaluator::pBCHandler
ssp< BoundaryHandler< model > > pBCHandler
gDim -> 3 for intellisense //!tmptmp ///< Variational reconstruction object.
Definition EulerEvaluator.hpp:179
DNDS::Euler::EulerEvaluator::UFromFace2Cell
void UFromFace2Cell(TU &u, index iFace, index iCell, rowsize if2c)
Inverse of UFromCell2Face: transform from face frame back to cell frame.
Definition EulerEvaluator.hpp:896
DNDS::Euler::EulerEvaluator::tFCompareField
std::function< TU(const Geom::tPoint &, real)> tFCompareField
Callback type for analytical comparison field.
Definition EulerEvaluator.hpp:601
DNDS::Euler::EulerEvaluator::DT_No_Flags
static const uint64_t DT_No_Flags
No flags for EvaluateDt.
Definition EulerEvaluator.hpp:339
DNDS::Euler::EulerEvaluator::UpdateLUSGSForward
void UpdateLUSGSForward(real alphaDiag, real t, ArrayDOFV< nVarsFixed > &rhs, ArrayDOFV< nVarsFixed > &u, ArrayDOFV< nVarsFixed > &uInc, JacobianDiagBlock< nVarsFixed > &JDiag, ArrayDOFV< nVarsFixed > &uIncNew)
Deprecated: use UpdateSGS with uIncIsZero=true instead. Kept for backward compatibility....
DNDS::Euler::EulerEvaluator::EvaluateNorm
void EvaluateNorm(Eigen::Vector< real, -1 > &res, ArrayDOFV< nVarsFixed > &rhs, index P=1, bool volWise=false, bool average=false)
Compute the norm of the RHS residual vector.
Definition EulerEvaluator.hxx:1394
DNDS::Euler::EulerEvaluator::bndIntegrationLogs
std::map< Geom::t_index, std::ofstream > bndIntegrationLogs
Per-BC log file streams for integration output.
Definition EulerEvaluator.hpp:204
DNDS::Euler::EulerEvaluator::symLU
ssp< Direct::SerialSymLUStructure > symLU
Symmetric LU structure for direct preconditioner.
Definition EulerEvaluator.hpp:218
DNDS::Euler::EulerEvaluator::getMuTur
real getMuTur(const TU &uMean, const TDiffU &GradUMeanXY, real muRef, real muf, index iFace)
Compute turbulent eddy viscosity at a face.
Definition EulerEvaluator.hpp:791
DNDS::Euler::EulerEvaluator::mesh
ssp< Geom::UnstructuredMesh > mesh
Shared pointer to the unstructured mesh.
Definition EulerEvaluator.hpp:177
DNDS::Euler::EulerEvaluator::gDim
static const int gDim
Geometric dimension of the mesh.
Definition EulerEvaluator.hpp:83
DNDS::Euler::EulerEvaluator::UFromCell2Face
void UFromCell2Face(TU &u, index iFace, index iCell, rowsize if2c)
Transform a conservative state vector from cell frame to face frame for periodic BCs.
Definition EulerEvaluator.hpp:879
DNDS::Euler::EulerEvaluator::CentralSmoothResidual
void CentralSmoothResidual(ArrayDOFV< nVarsFixed > &r, ArrayDOFV< nVarsFixed > &rs, ArrayDOFV< nVarsFixed > &rtemp, int nStep=0)
Apply Laplacian smoothing to a residual field.
Definition EulerEvaluator.hpp:1923
DNDS::Euler::EulerEvaluator::updateBCProfilesPressureRadialEq
void updateBCProfilesPressureRadialEq()
Update boundary profiles using radial-equilibrium pressure extrapolation.
Definition EulerEvaluator.hxx:2295
DNDS::Euler::EulerEvaluator::MeanValueCons2Prim
void MeanValueCons2Prim(ArrayDOFV< nVarsFixed > &u, ArrayDOFV< nVarsFixed > &w)
Convert cell-mean conservative variables to primitive variables.
Definition EulerEvaluator.hxx:1348
DNDS::Euler::EulerEvaluator::dim
static const int dim
Spatial dimension (2 or 3).
Definition EulerEvaluator.hpp:82
DNDS::Euler::EulerEvaluator::fluxJacobian0_Right_Times_du
TU fluxJacobian0_Right_Times_du(const TU &U, const TU &UOther, const TVec &n, const TVec &vg, Geom::t_index btype, const TU &dU, real lambdaMain, real lambdaC, real lambdaVis, real lambda0, real lambda123, real lambda4, int useRoeTerm, int incFsign=-1, int omitF=0)
inviscid flux approx jacobian (flux term not reconstructed / no riemann) if lambdaMain == veryLargeRe...
Definition EulerEvaluator.hpp:1100
DNDS::Euler::EulerEvaluator::GetFaceVGrid
TVec GetFaceVGrid(index iFace, index iG, const Geom::tPoint &pPhy)
Get the grid velocity at a face quadrature point (with explicit physical point).
Definition EulerEvaluator.hpp:1221
DNDS::Euler::EulerEvaluator::LUSGSMatrixVec
void LUSGSMatrixVec(real alphaDiag, real t, ArrayDOFV< nVarsFixed > &u, ArrayDOFV< nVarsFixed > &uInc, JacobianDiagBlock< nVarsFixed > &JDiag, ArrayDOFV< nVarsFixed > &AuInc)
Compute the matrix-vector product A * uInc for the implicit system.
DNDS::Euler::EulerEvaluator::EvaluateCellRHSAlpha_MIN_ALL
static const int EvaluateCellRHSAlpha_MIN_ALL
Always take min(alpha, prev).
Definition EulerEvaluator.hpp:693
DNDS::Euler::EulerEvaluator::visFluxTurVariable
void visFluxTurVariable(const TU &UMeanXYC, const TDiffU &DiffUxyPrimC, real muRef, real mufPhy, real muTur, const TVec &uNormC, index iFace, TU &VisFlux)
Compute viscous flux contribution from turbulence model variables.
Definition EulerEvaluator.hpp:835
DNDS::Euler::EulerEvaluator::tFCompareFieldWeight
std::function< real(const Geom::tPoint &, real)> tFCompareFieldWeight
Callback type for comparison weighting function.
Definition EulerEvaluator.hpp:602
DNDS::Euler::EulerEvaluator::lambdaFaceVis
std::vector< real > lambdaFaceVis
Per-face viscous spectral radius.
Definition EulerEvaluator.hpp:186
DNDS::Euler::EulerEvaluator::lambdaFace
std::vector< real > lambdaFace
Per-face spectral radius (inviscid + viscous combined).
Definition EulerEvaluator.hpp:184
DNDS::Euler::EulerEvaluator::ConsoleOutputBndIntegrations
void ConsoleOutputBndIntegrations()
Print boundary integration results (force/flux) to console on rank 0.
Definition EulerEvaluator.hpp:1447
DNDS::Euler::EulerEvaluator::LUSGSMatrixSolveJacobianLU
void LUSGSMatrixSolveJacobianLU(real alphaDiag, real t, ArrayDOFV< nVarsFixed > &rhs, ArrayDOFV< nVarsFixed > &u, ArrayDOFV< nVarsFixed > &uInc, ArrayDOFV< nVarsFixed > &uIncNew, ArrayDOFV< nVarsFixed > &bBuf, JacobianDiagBlock< nVarsFixed > &JDiag, JacobianLocalLU< nVarsFixed > &jacLU, bool uIncIsZero, TU &sumInc)
Solve the implicit linear system using the pre-factored local LU.
DNDS::Euler::EulerEvaluator::fluxJacobian0_Right_Times_du_AsMatrix
TJacobianU fluxJacobian0_Right_Times_du_AsMatrix(const TU &U, const TU &UOther, const TVec &n, const TVec &vg, Geom::t_index btype, real lambdaMain, real lambdaC, real lambdaVis, real lambda0, real lambda123, real lambda4, int useRoeTerm, int incFsign=-1, int omitF=0)
Definition EulerEvaluator.hpp:1160
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::TransformVelocityRotatingFrame
void TransformVelocityRotatingFrame(TU &U, const Geom::tPoint &pPhysics, int direction)
Transform momentum in-place between inertial and rotating frame (velocity only).
Definition EulerEvaluator.hpp:1248
DNDS::Euler::EulerEvaluator::TVec
Eigen::VectorFMTSafe< real, dim > TVec
Spatial vector (dim components).
Definition EulerEvaluator.hpp:90
DNDS::Euler::EulerEvaluator::LimiterUGrad
void LimiterUGrad(ArrayDOFV< nVarsFixed > &u, ArrayGRADV< nVarsFixed, gDim > &uGrad, ArrayGRADV< nVarsFixed, gDim > &uGradNew, uint64_t flags=LIMITER_UGRAD_No_Flags)
Apply slope limiter to the gradient field.
DNDS::Euler::EulerEvaluator::fluxJacobianC_Right_Times_du
TU fluxJacobianC_Right_Times_du(const TU &U, const TVec &n, const TVec &vg, Geom::t_index btype, const TU &dU)
Definition EulerEvaluator.hpp:1182
DNDS::Euler::EulerEvaluator::EvaluateCellRHSAlpha
void EvaluateCellRHSAlpha(ArrayDOFV< nVarsFixed > &u, ArrayRECV< nVarsFixed > &uRec, ArrayDOFV< 1 > &uRecBeta, ArrayDOFV< nVarsFixed > &res, ArrayDOFV< 1 > &cellRHSAlpha, index &nLim, real &alphaMin, real relax, int compress=1, int flag=0)
Compute the positivity-preserving RHS scaling factor (alpha) per cell.
Definition EulerEvaluator.hxx:1938
DNDS::Euler::EulerEvaluator::updateBCAnchors
void updateBCAnchors(ArrayDOFV< nVarsFixed > &u, ArrayRECV< nVarsFixed > &uRec)
Update boundary anchor point recorders from current solution.
Definition EulerEvaluator.hpp:1277
DNDS::Euler::EulerEvaluator::TScalar
ArrayRECV< 1 > TScalar
Scalar reconstruction coefficient array.
Definition EulerEvaluator.hpp:103
DNDS::Euler::EulerEvaluator::bndIntegrations
std::map< Geom::t_index, IntegrationRecorder > bndIntegrations
Per-BC boundary flux/force integration accumulators.
Definition EulerEvaluator.hpp:203
DNDS::Euler::EulerEvaluator::updateBCProfiles
void updateBCProfiles(ArrayDOFV< nVarsFixed > &u, ArrayRECV< nVarsFixed > &uRec)
Update boundary 1-D profiles from the current solution.
Definition EulerEvaluator.hxx:2195
DNDS::Euler::EulerEvaluator::TransformURotatingFrame_ABS_VELO
void TransformURotatingFrame_ABS_VELO(TU &U, const Geom::tPoint &pPhysics, int direction)
Transform full conservative state for the absolute-velocity rotating frame formulation.
Definition EulerEvaluator.hpp:1266
DNDS::Euler::EulerEvaluator::TDiffU_Batch
Eigen::MatrixFMTSafe< real, Eigen::Dynamic, nVarsFixed, Eigen::ColMajor, MaxBatchMult(3)> TDiffU_Batch
Batch of gradient matrices.
Definition EulerEvaluator.hpp:99
DNDS::Euler::EulerEvaluator::EulerEvaluator
EulerEvaluator(const decltype(mesh) &Nmesh, const decltype(vfv) &Nvfv, const decltype(pBCHandler) &npBCHandler, const EulerEvaluatorSettings< model > &nSettings, int n_nVars=getNVars(model))
Construct an EulerEvaluator and initialize all internal buffers.
Definition EulerEvaluator.hpp:236
DNDS::Euler::EulerEvaluator::RHS_Dont_Update_Integration
static const uint64_t RHS_Dont_Update_Integration
Skip boundary integration accumulation.
Definition EulerEvaluator.hpp:372
DNDS::Euler::EulerEvaluator::DiffUFromCell2Face
void DiffUFromCell2Face(TDiffU &u, index iFace, index iCell, rowsize if2c, bool reverse=false)
Transform a gradient tensor from cell frame to face frame for periodic BCs.
Definition EulerEvaluator.hpp:928
DNDS::Euler::EulerEvaluator::TimeAverageAddition
void TimeAverageAddition(ArrayDOFV< nVarsFixed > &w, ArrayDOFV< nVarsFixed > &wAveraged, real dt, real &tCur)
Accumulate time-averaged primitive variables for unsteady statistics.
Definition EulerEvaluator.hxx:1335
DNDS::Euler::EulerEvaluator::fluxBnd
std::vector< TU > fluxBnd
Per-boundary-face flux values.
Definition EulerEvaluator.hpp:214
DNDS::Euler::EulerEvaluator::EvaluateCellRHSAlpha_MIN_IF_NOT_ONE
static const int EvaluateCellRHSAlpha_MIN_IF_NOT_ONE
Take min(alpha, prev) only if prev != 1.
Definition EulerEvaluator.hpp:692
DNDS::Euler::EulerEvaluator::lambdaFace4
std::vector< real > lambdaFace4
Per-face eigenvalue |u·n - a| (acoustic wave).
Definition EulerEvaluator.hpp:189
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::UFromOtherCell
void UFromOtherCell(TU &u, index iFace, index iCell, index iCellOther, rowsize if2c)
Transform a state from a neighbor cell across a periodic face.
Definition EulerEvaluator.hpp:913
DNDS::Euler::EulerEvaluator::EvaluateCellRHSAlphaExpansion
void EvaluateCellRHSAlphaExpansion(ArrayDOFV< nVarsFixed > &u, ArrayRECV< nVarsFixed > &uRec, ArrayDOFV< 1 > &uRecBeta, ArrayDOFV< nVarsFixed > &res, ArrayDOFV< 1 > &cellRHSAlpha, index &nLim, real alphaMin)
Expand a previously computed cellRHSAlpha toward 1 where safe.
Definition EulerEvaluator.hxx:2064
DNDS::Euler::EulerEvaluator::TRec
ArrayRECV< nVarsFixed > TRec
Reconstruction coefficient array (per-cell polynomial coefficients).
Definition EulerEvaluator.hpp:102
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::EulerEvaluator::dWall
std::vector< Eigen::Vector< real, Eigen::Dynamic > > dWall
Per-cell wall distance (one value per cell node/quadrature point).
Definition EulerEvaluator.hpp:197
DNDS::Euler::EulerEvaluator::RHS_No_Flags
static const uint64_t RHS_No_Flags
Default: full RHS evaluation.
Definition EulerEvaluator.hpp:370
DNDS::Euler::EulerEvaluator::fluxJacobian0_Right
auto fluxJacobian0_Right(TU &UR, const TVec &uNorm, Geom::t_index btype)
inviscid flux approx jacobian (flux term not reconstructed / no riemann)
Definition EulerEvaluator.hpp:1028
DNDS::Euler::EulerEvaluator::UpdateSGSWithRec
void UpdateSGSWithRec(real alphaDiag, real t, ArrayDOFV< nVarsFixed > &rhs, ArrayDOFV< nVarsFixed > &u, ArrayRECV< nVarsFixed > &uRec, ArrayDOFV< nVarsFixed > &uInc, ArrayRECV< nVarsFixed > &uRecInc, JacobianDiagBlock< nVarsFixed > &JDiag, bool forward, TU &sumInc)
SGS sweep coupled with reconstruction update.
Definition EulerEvaluator.hxx:713
DNDS::Euler::JacobianDiagBlock
Per-cell diagonal or block-diagonal Jacobian storage for implicit time stepping.
Definition EulerJacobian.hpp:43
DNDS::Euler::OutputPicker
Registry that maps named cell-scalar output fields to getter lambdas.
Definition EulerBC.hpp:716
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::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_2EQ
@ NS_2EQ
NS + 2-equation RANS, 2D geometry (7 vars).
Definition Euler.hpp:881
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::Geom::FaceIDIsPeriodic
DNDS_DEVICE_CALLABLE bool FaceIDIsPeriodic(t_index id)
Definition BoundaryCondition.hpp:95
DNDS::Geom::FaceIDIsPeriodicDonor
DNDS_DEVICE_CALLABLE bool FaceIDIsPeriodicDonor(t_index id)
Definition BoundaryCondition.hpp:130
DNDS::Geom::FaceIDIsPeriodicMain
DNDS_DEVICE_CALLABLE bool FaceIDIsPeriodicMain(t_index id)
Definition BoundaryCondition.hpp:123
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::TermColor::Magenta
constexpr std::string_view Magenta
ANSI escape: bright magenta foreground.
Definition Defines.hpp:877
DNDS::TermColor::Reset
constexpr std::string_view Reset
ANSI escape: reset all attributes.
Definition Defines.hpp:883
DNDS::cube
constexpr T cube(const T &a)
a * a * a, constexpr.
Definition Defines.hpp:525
DNDS::UnInitIndex
DNDS_CONSTANT const index UnInitIndex
Sentinel "not initialised" index value (= INT64_MIN).
Definition Defines.hpp:176
DNDS::verySmallReal
DNDS_CONSTANT const real verySmallReal
Catch-all lower bound ("effectively zero").
Definition Defines.hpp:194
DNDS::rowsize
int32_t rowsize
Row-width / per-row element-count type (signed 32-bit).
Definition Defines.hpp:109
DNDS::DynamicSize
DNDS_CONSTANT const rowsize DynamicSize
Template parameter flag: "row width is set at runtime but uniform".
Definition Defines.hpp:277
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::ssp
std::shared_ptr< T > ssp
Shortened alias for std::shared_ptr used pervasively in DNDSR.
Definition Defines.hpp:138
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::DNDS_MPI_REAL
const MPI_Datatype DNDS_MPI_REAL
MPI datatype matching real (= MPI_REAL8).
Definition MPI.hpp:92
DNDS::ArrayPair::father
ssp< TArray > father
Owned-side array (must be resized before ghost setup).
Definition ArrayPair.hpp:170
DNDS::ArrayPair::Size
index Size() const
Combined row count (father->Size() + son->Size()).
Definition ArrayPair.hpp:292
DNDS::ArrayPair::trans
TTrans trans
Ghost-communication engine bound to father and son.
Definition ArrayPair.hpp:175
DNDS::Euler::AnchorPointRecorder
Finds the cell closest to a specified anchor point across all MPI ranks.
Definition EulerBC.hpp:647
DNDS::Euler::EulerEvaluatorSettings
Master configuration struct for the compressible Euler/Navier-Stokes evaluator.
Definition EulerEvaluatorSettings.hpp:44
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
DNDS::Euler::EulerModelTraits
Compile-time traits for EulerModel variants.
Definition Euler.hpp:1086
DNDS::Euler::EulerModelTraits::hasSA
static constexpr bool hasSA
True for Spalart-Allmaras models (NS_SA, NS_SA_3D).
Definition Euler.hpp:1095
DNDS::Euler::EulerModelTraits::has2EQ
static constexpr bool has2EQ
True for 2-equation RANS models (NS_2EQ, NS_2EQ_3D).
Definition Euler.hpp:1097
DNDS::Euler::JacobianLocalLU
Complete LU factorization of the cell-local Jacobian for GMRES preconditioning.
Definition EulerJacobian.hpp:273
Eigen::MatrixFMTSafe
Eigen::Matrix wrapper that hides begin/end from fmt.
Definition EigenUtil.hpp:62
v
Eigen::Matrix< real, 5, 1 > v
Definition test_ArrayDOF.cpp:468
r
tVec r(NCells)
alpha
real alpha
Definition test_GasThermo.cpp:454
vg
Eigen::Vector3d vg(0, 0, 0)
velo
Eigen::Vector3d velo
Definition test_GasThermo.cpp:282
dU
Eigen::Vector< real, 5 > dU
Definition test_GasThermo.cpp:452
res
auto res
Definition test_ODE.cpp:196
n
Eigen::Vector3d n(1.0, 0.0, 0.0)