DNDSR 0.2.1
Distributed Numeric Data Structure for CFV
Loading...
Searching...
No Matches
EulerSolver_Init.hxx
Go to the documentation of this file.
1/** @file EulerSolver_Init.hxx
2 * @brief Template implementation of EulerSolver::ReadMeshAndInitialize, the full
3 * solver initialization pipeline from mesh reading to evaluator setup.
4 *
5 * Covers CGNS/OpenFOAM mesh reading, periodic boundary deduplication, mesh
6 * partitioning (ParMetis), O1-to-O2 elevation, h-refinement bisection, boundary
7 * mesh extraction for wall distance, VFV (VariationalReconstruction) construction,
8 * BC handler configuration, wall distance computation, evaluator initialization,
9 * restart loading, and initial VTK output.
10 */
11#pragma once
12
13#include "EulerSolver.hpp"
14#include "SpecialFields.hpp"
15#include "DNDS/OMP.hpp"
16#include "Geom/Mesh/Mesh_Helpers.hpp"
17
18namespace DNDS::Euler
19{
20 static const auto model = NS_SA;
21 DNDS_SWITCH_INTELLISENSE(template <EulerModel model>, )
22 /** @brief Read the mesh, partition it, build solver data structures, and initialize the evaluator.
23 *
24 * Complete initialization pipeline:
25 * 1. Read mesh from CGNS or OpenFOAM format (serial, parallel, or distributed mode).
26 * 2. Handle periodic boundary deduplication and mesh topology (cell2cell, node2cell).
27 * 3. Optionally elevate mesh order (O1 to O2) and apply h-refinement bisection.
28 * 4. Partition with ParMetis and redistribute across MPI ranks.
29 * 5. Build ghost layers, adjacency, and coordinate structures.
30 * 6. Extract boundary mesh for wall distance computation.
31 * 7. Construct VFV (VariationalReconstruction) with configured settings.
32 * 8. Configure BC handlers from JSON settings.
33 * 9. Compute wall distance (CGAL or Poisson).
34 * 10. Initialize the EulerEvaluator (arrays, face data, etc.).
35 * 11. Load restart if configured.
36 * 12. Write initial VTK output.
37 */
38 void EulerSolver<model>::ReadMeshAndInitialize()
39 {
40 DNDS_MPI_InsertCheck(mpi, "ReadMeshAndInitialize 1 nvars " + std::to_string(nVars));
41 output_stamp = getTimeStamp(mpi);
42 if (!config.dataIOControl.uniqueStamps)
43 output_stamp = "";
44 if (mpi.rank == 0)
45 log() << "=== Got Time Stamp: [" << output_stamp << "] ===" << std::endl;
46 // Debug::MPIDebugHold(mpi);
47
48 int gDimLocal = gDim; //! or else the linker breaks down here (with clang++ or g++, -g -O0,2; c++ non-optimizer bug?)
49
50 auto &BCHandler = *pBCHandler;
51
52 DNDS_MAKE_SSP(mesh, mpi, gDimLocal);
53 DNDS_MAKE_SSP(meshBnd, mpi, gDimLocal - 1);
54
55 DNDS_MAKE_SSP(reader, mesh, 0);
56 DNDS_MAKE_SSP(readerBnd, meshBnd, 0);
57 DNDS_assert(config.dataIOControl.readMeshMode == 0 || config.dataIOControl.readMeshMode == 1 || config.dataIOControl.readMeshMode == 2);
58 DNDS_assert(config.dataIOControl.outPltMode == 0 || config.dataIOControl.outPltMode == 1);
59 mesh->periodicInfo.translation[1].map() = config.boundaryDefinition.PeriodicTranslation1;
60 mesh->periodicInfo.translation[2].map() = config.boundaryDefinition.PeriodicTranslation2;
61 mesh->periodicInfo.translation[3].map() = config.boundaryDefinition.PeriodicTranslation3;
62 mesh->periodicInfo.rotationCenter[1].map() = config.boundaryDefinition.PeriodicRotationCent1;
63 mesh->periodicInfo.rotationCenter[2].map() = config.boundaryDefinition.PeriodicRotationCent2;
64 mesh->periodicInfo.rotationCenter[3].map() = config.boundaryDefinition.PeriodicRotationCent3;
65 mesh->periodicInfo.rotation[1].map() =
66 Geom::RotZ(config.boundaryDefinition.PeriodicRotationEulerAngles1[2]) *
67 Geom::RotY(config.boundaryDefinition.PeriodicRotationEulerAngles1[1]) *
68 Geom::RotX(config.boundaryDefinition.PeriodicRotationEulerAngles1[0]);
69 mesh->periodicInfo.rotation[2].map() =
70 Geom::RotZ(config.boundaryDefinition.PeriodicRotationEulerAngles2[2]) *
71 Geom::RotY(config.boundaryDefinition.PeriodicRotationEulerAngles2[1]) *
72 Geom::RotX(config.boundaryDefinition.PeriodicRotationEulerAngles2[0]);
73 mesh->periodicInfo.rotation[3].map() =
74 Geom::RotZ(config.boundaryDefinition.PeriodicRotationEulerAngles3[2]) *
75 Geom::RotY(config.boundaryDefinition.PeriodicRotationEulerAngles3[1]) *
76 Geom::RotX(config.boundaryDefinition.PeriodicRotationEulerAngles3[0]);
77
78 if (config.dataIOControl.readMeshMode == 0)
79 {
80 if (config.dataIOControl.meshFormat == 1)
81 {
82 reader->ReadFromOpenFOAMAndConvertSerial(
83 config.dataIOControl.meshFile,
84 config.bcNameMapping,
85 [&](const std::string &name)
86 -> Geom::t_index
87 { return BCHandler.GetIDFromName(name); });
88 reader->Deduplicate1to1Periodic(config.boundaryDefinition.periodicTolerance);
89 reader->BuildCell2Cell();
90 reader->MeshPartitionCell2Cell(config.dataIOControl.meshPartitionOptions);
91 reader->PartitionReorderToMeshCell2Cell();
92
93 mesh->RecoverNode2CellAndNode2Bnd();
94 mesh->RecoverCell2CellAndBnd2Cell();
95 mesh->BuildGhostPrimary();
96 mesh->AdjGlobal2LocalPrimary();
97 mesh->AdjGlobal2LocalN2CB();
98 if (config.dataIOControl.meshElevation == 1)
99 {
100 DNDS::ssp<DNDS::Geom::UnstructuredMesh> meshO2;
101 DNDS_MAKE_SSP(meshO2, mpi, gDimLocal);
102 meshO2->BuildO2FromO1Elevation(*mesh);
103 std::swap(meshO2, mesh);
104
105 reader->mesh = mesh;
106 mesh->RecoverNode2CellAndNode2Bnd();
107 mesh->RecoverCell2CellAndBnd2Cell();
108 mesh->BuildGhostPrimary();
109 mesh->AdjGlobal2LocalPrimary();
110 mesh->AdjGlobal2LocalN2CB();
111 }
112 DNDS_assert(config.dataIOControl.meshDirectBisect <= 4);
113 for (int iter = 1; iter <= config.dataIOControl.meshDirectBisect; iter++)
114 {
115 DNDS::ssp<DNDS::Geom::UnstructuredMesh> meshO2;
116 DNDS_MAKE_SSP(meshO2, mpi, gDimLocal);
117 meshO2->BuildO2FromO1Elevation(*mesh);
118
119 meshO2->RecoverNode2CellAndNode2Bnd();
120 meshO2->RecoverCell2CellAndBnd2Cell();
121 meshO2->BuildGhostPrimary();
122 DNDS::ssp<DNDS::Geom::UnstructuredMesh> meshO1B;
123 DNDS_MAKE_SSP(meshO1B, mpi, gDimLocal);
124 meshO1B->BuildBisectO1FormO2(*meshO2);
125
126 std::swap(meshO1B, mesh);
127 reader->mesh = mesh;
128 mesh->RecoverNode2CellAndNode2Bnd();
129 mesh->RecoverCell2CellAndBnd2Cell();
130 mesh->BuildGhostPrimary();
131 mesh->AdjGlobal2LocalPrimary();
132 mesh->AdjGlobal2LocalN2CB();
133 index nCell = mesh->NumCellGlobal();
134 index nNode = mesh->NumNodeGlobal();
135 if (mesh->getMPI().rank == 0)
136 {
137 log() << fmt::format("Mesh Direct Bisect {} done, nCell [{}], nNode [{}]", iter, nCell, nNode) << std::endl;
138 }
139 }
140 }
141 else
142 {
143 Geom::ReadMeshFromCGNS(
144 mesh, reader,
145 config.dataIOControl.meshFile,
146 config.dataIOControl.meshPartitionOptions,
147 config.boundaryDefinition.periodicTolerance,
148 config.dataIOControl.meshElevation,
149 config.dataIOControl.meshDirectBisect,
150 [&](const std::string &name) -> Geom::t_index
151 { return BCHandler.GetIDFromName(name); });
152 }
153 }
154 else if (config.dataIOControl.readMeshMode == 1)
155 {
156 std::string meshOutName = Geom::MeshH5Path(
157 config.dataIOControl.meshFile, mpi.size,
158 config.dataIOControl.meshElevation,
159 config.dataIOControl.meshDirectBisect);
160 Geom::ReadMeshFromH5Parallel(
161 *mesh,
162 Serializer::SerializerFactory(config.dataIOControl.meshPartitionedReaderType),
163 meshOutName);
164 }
165 else if (config.dataIOControl.readMeshMode == 2)
166 {
167 std::string meshOutName = Geom::MeshH5Path(
168 config.dataIOControl.meshFile, mpi.size,
169 config.dataIOControl.meshElevation,
170 config.dataIOControl.meshDirectBisect);
171 Geom::ReadMeshFromH5(
172 *mesh,
173 Serializer::SerializerFactory(config.dataIOControl.meshPartitionedReaderType),
174 meshOutName,
175 config.dataIOControl.meshPartitionOptions);
176 }
177 OMP::set_full_parallel_OMP();
178 {
179 Geom::PrepareMeshOptions prepOpts;
180 prepOpts.reorderCells = config.dataIOControl.meshReorderCells;
181#ifdef DNDS_USE_OMP
182 prepOpts.reorderParts = std::max(get_env_DNDS_DIST_OMP_NUM_THREADS(), 1);
183#else
184 prepOpts.reorderParts = 1;
185#endif
186 prepOpts.buildSerialOut = (config.dataIOControl.outPltMode == 0);
187 Geom::PrepareMesh(*mesh, *reader, prepOpts);
188 }
189
190 // note: node2cell may contain negative entries for ghost nodes
191 // whose cell neighbors are not in the local partition.
192 // The old code checked iCell >= 0 here but only for NumNode()
193 // (father nodes), not NumNodeProc() (father + son).
194
195 if (config.dataIOControl.meshElevation == 1 && config.dataIOControl.readMeshMode == 0)
196 {
197 mesh->elevationInfo.nIter = config.dataIOControl.meshElevationIter;
198 mesh->elevationInfo.nSearch = config.dataIOControl.meshElevationNSearch;
199 mesh->elevationInfo.RBFRadius = config.dataIOControl.meshElevationRBFRadius;
200 mesh->elevationInfo.RBFPower = config.dataIOControl.meshElevationRBFPower;
201 mesh->elevationInfo.kernel = config.dataIOControl.meshElevationRBFKernel;
202 mesh->elevationInfo.MaxIncludedAngle = config.dataIOControl.meshElevationMaxIncludedAngle;
203 mesh->elevationInfo.refDWall = config.dataIOControl.meshElevationRefDWall;
204 mesh->ElevatedNodesGetBoundarySmooth(
205 [&](Geom::t_index bndId)
206 {
207 auto bType = pBCHandler->GetTypeFromID(bndId);
208 if (bType == BCWall || bType == BCWallIsothermal)
209 return true;
210 if (config.dataIOControl.meshElevationBoundaryMode == 1 &&
211 (bType == BCWallInvis || bType == BCSym))
212 return true;
213 return false;
214 });
215 if (config.dataIOControl.meshElevationInternalSmoother == 0)
216 mesh->ElevatedNodesSolveInternalSmooth();
217 else if (config.dataIOControl.meshElevationInternalSmoother == 1)
218 mesh->ElevatedNodesSolveInternalSmoothV1();
219 else if (config.dataIOControl.meshElevationInternalSmoother == 2)
220 mesh->ElevatedNodesSolveInternalSmoothV2();
221 else if (config.dataIOControl.meshElevationInternalSmoother == -1)
222 {
223 if (mpi.rank == 0)
224 log() << " WARNING !!! Not Smoothing internal, abandoning boundary smooth displacements" << std::endl;
225 }
226 else
227 DNDS_assert(false);
228 }
229
230 if (config.dataIOControl.meshBuildWallDist)
231 {
232 mesh->BuildNodeWallDist(
233 [pBCHandler = pBCHandler](Geom::t_index id)
234 {
235 auto euler_bc_type = pBCHandler->GetTypeFromID(id);
236 return static_cast<bool>(euler_bc_type == Euler::BCWall || euler_bc_type == Euler::BCWallIsothermal);
237 },
238 config.dataIOControl.meshWallDistOptions);
239 }
240
241 if (config.timeMarchControl.partitionMeshOnly)
242 {
243 std::string meshOutName = Geom::MeshH5Path(
244 config.dataIOControl.meshFile, mpi.size,
245 config.dataIOControl.meshElevation,
246 config.dataIOControl.meshDirectBisect);
247 Geom::SerializeMesh(*mesh, meshOutName, config.dataIOControl.meshPartitionedWriter);
248 return; //** mesh preprocess only (without transformation)
249 }
250
251 Geom::BuildBndMesh(*mesh, *meshBnd, *readerBnd,
252 config.dataIOControl.outPltMode == 0);
253
254 if (config.dataIOControl.meshRotZ != 0.0)
255 {
256 real rz = config.dataIOControl.meshRotZ / 180.0 * pi;
257 Eigen::Matrix3d Rz{
258 {std::cos(rz), -std::sin(rz), 0},
259 {std::sin(rz), std::cos(rz), 0},
260 {0, 0, 1},
261 };
262 mesh->TransformCoords(
263 [&](const Geom::tPoint &p)
264 { return Rz * p; });
265 meshBnd->TransformCoords(
266 [&](const Geom::tPoint &p)
267 { return Rz * p; });
268
269 for (auto &r : mesh->periodicInfo.rotation)
270 r.map() = Rz * r.map() * Rz.transpose();
271 for (auto &p : mesh->periodicInfo.rotationCenter)
272 p.map() = Rz * p.map();
273 for (auto &p : mesh->periodicInfo.translation)
274 p.map() = Rz * p.map();
275 // @todo //! todo: alter the rotation and translation in periodicInfo mesh->periodicInfo
276 }
277 if (config.dataIOControl.meshScale != 1.0)
278 {
279 auto scale = config.dataIOControl.meshScale;
280 mesh->TransformCoords(
281 [&](const Geom::tPoint &p)
282 { return p * scale; });
283 meshBnd->TransformCoords(
284 [&](const Geom::tPoint &p)
285 { return p * scale; });
286
287 for (auto &i : mesh->periodicInfo.translation)
288 i.map() *= scale;
289 for (auto &i : mesh->periodicInfo.rotationCenter)
290 i.map() *= scale;
291 }
292 if (config.dataIOControl.rectifyNearPlane)
293 {
294 auto fTrans = [&](const Geom::tPoint &p)
295 {
296 Geom::tPoint ret = p;
297 if (config.dataIOControl.rectifyNearPlane & 1)
298 if (std::abs(ret(0)) < config.dataIOControl.rectifyNearPlaneThres)
299 ret(0) = 0;
300 if (config.dataIOControl.rectifyNearPlane & 2)
301 if (std::abs(ret(1)) < config.dataIOControl.rectifyNearPlaneThres)
302 ret(1) = 0;
303 if (config.dataIOControl.rectifyNearPlane & 4)
304 if (std::abs(ret(2)) < config.dataIOControl.rectifyNearPlaneThres)
305 ret(2) = 0;
306 return ret;
307 };
308 mesh->TransformCoords(fTrans);
309 meshBnd->TransformCoords(fTrans);
310 }
311 { //* symBnd's rectifying: ! altering mesh
312 for (index iB = 0; iB < mesh->NumBnd(); iB++)
313 {
314 index iFace = mesh->bnd2faceV.at(iB);
315 auto bndID = mesh->bndElemInfo(iB, 0).zone;
316 EulerBCType bndType = pBCHandler->GetTypeFromID(bndID);
317 if (bndType == BCSym)
318 {
319 auto rectifyOpt = pBCHandler->GetFlagFromID(bndID, "rectifyOpt");
320 if (rectifyOpt >= 1 && rectifyOpt <= 3)
321 for (auto iNode : mesh->bnd2node[iB])
322 mesh->coords[iNode](rectifyOpt - 1) = 0.0;
323 }
324 }
325 mesh->coords.trans.pullOnce();
326 for (index iB = 0; iB < meshBnd->NumCell(); iB++)
327 {
328 auto bndID = meshBnd->cellElemInfo(iB, 0).zone;
329 EulerBCType bndType = pBCHandler->GetTypeFromID(bndID);
330 if (bndType == BCSym)
331 {
332 auto rectifyOpt = pBCHandler->GetFlagFromID(bndID, "rectifyOpt");
333 if (rectifyOpt >= 1 && rectifyOpt <= 3)
334 for (auto iNode : meshBnd->cell2node[iB])
335 meshBnd->coords[iNode](rectifyOpt - 1) = 0.0;
336 }
337 }
338 }
339 /// @todo //todo: upgrade to optional
340 if (config.dataIOControl.outPltMode == 0)
341 reader->coordSerialOutTrans.pullOnce(),
342 readerBnd->coordSerialOutTrans.pullOnce();
343
344 mesh->RecreatePeriodicNodes();
345 mesh->BuildVTKConnectivity();
346 meshBnd->RecreatePeriodicNodes();
347 meshBnd->BuildVTKConnectivity();
348
349 // *demo code
350 // mesh->PrintMeshCGNS("test.cgns", [&](Geom::t_index i)
351 // { return BCHandler.GetNameFormID(i); }, BCHandler.GetAllNames());
352
353 DNDS_MAKE_SSP(vfv, mpi, mesh);
354 vfv->parseSettings(config.vfvSettings);
355 if (vfv->SetAxisSymmetric(config.others.axisSymmetric) && mpi.rank == 0)
356 log() << "EulerSolver === Using Axis Symmetric" << std::endl;
357 TEval::InitializeFV(mesh, vfv, pBCHandler);
358 if (config.others.printRecMatrix)
359 {
360 auto serializerP = config.others.recMatrixWriter.BuildSerializer(mpi);
361 std::string fName = config.dataIOControl.outPltName + "_RecMatrix";
362 auto [fNameMod, partPath] = config.others.recMatrixWriter.ModifyFilePath(fName, mpi, "part_%d", true);
363 serializerP->OpenFile(partPath, false);
364 vfv->WriteSerializeRecMatrix(serializerP);
365 serializerP->CloseFile();
366 }
367
368 vfv->BuildUDof(u, nVars);
369 vfv->BuildUDof(uIncBufODE, nVars);
370 if (config.timeAverageControl.enabled)
371 {
372 vfv->BuildUDof(wAveraged, nVars);
373 vfv->BuildUDof(uAveraged, nVars);
374 }
375 uPool.resizeInit(3,
376 [&](ArrayDOFV<nVarsFixed> &uu)
377 { vfv->BuildUDof(uu, nVars); });
378
379 vfv->BuildURec(uRec, nVars);
380 if (config.timeMarchControl.timeMarchIsTwoStage())
381 vfv->BuildURec(uRec1, nVars);
382 vfv->BuildURec(uRecLimited, nVars);
383 vfv->BuildURec(uRecNew, nVars);
384 vfv->BuildURec(uRecNew1, nVars);
385 vfv->BuildURec(uRecB, nVars);
386 vfv->BuildURec(uRecB1, nVars);
387 vfv->BuildURec(uRecOld, nVars);
388 vfv->BuildScalar(ifUseLimiter);
389 vfv->BuildUDof(betaPP, 1);
390 vfv->BuildUDof(alphaPP, 1);
391 vfv->BuildUDof(alphaPP_tmp, 1);
392 vfv->BuildUDof(dTauTmp, 1);
393 betaPP.setConstant(1.0);
394 alphaPP.setConstant(1.0);
395 if (config.timeMarchControl.timeMarchIsTwoStage())
396 {
397 vfv->BuildUDof(betaPP1, 1);
398 vfv->BuildUDof(alphaPP1, 1);
399 betaPP1.setConstant(1.0);
400 alphaPP1.setConstant(1.0);
401 }
402 if (config.implicitReconstructionControl.storeRecInc)
403 {
404 vfv->BuildURec(uRecInc, nVars);
405 if (config.timeMarchControl.timeMarchIsTwoStage())
406 vfv->BuildURec(uRecInc1, nVars);
407 }
408
409 DNDS_MPI_InsertCheck(mpi, "ReadMeshAndInitialize 2 nvars " + std::to_string(nVars));
410 /*******************************/
411 // initialize pEval
412 DNDS_MAKE_SSP(pEval, mesh, vfv, pBCHandler, config.eulerSettings, nVars);
413 EulerEvaluator<model> &eval = *pEval;
414
415 JD.SetModeAndInit(eval.settings.useScalarJacobian ? 0 : 1, nVars, u);
416 JSource.SetModeAndInit(eval.settings.useScalarJacobian ? 0 : 1, nVars, u);
417 if (config.timeMarchControl.timeMarchIsTwoStage())
418 {
419 JD1.SetModeAndInit(eval.settings.useScalarJacobian ? 0 : 1, nVars, u);
420 JSource1.SetModeAndInit(eval.settings.useScalarJacobian ? 0 : 1, nVars, u);
421 }
422 JDTmp.SetModeAndInit(eval.settings.useScalarJacobian ? 0 : 1, nVars, u);
423 JSourceTmp.SetModeAndInit(eval.settings.useScalarJacobian ? 0 : 1, nVars, u);
424 /*******************************/
425 // ** initialize output Array
426
427 DNDS_MPI_InsertCheck(mpi, "ReadMeshAndInitialize 3 nvars " + std::to_string(nVars));
428
429 // update output number
430 DNDS_assert(config.dataIOControl.outCellScalarNames.size() < 128);
431 nOUTS += config.dataIOControl.outCellScalarNames.size();
432
433 DNDS_assert(config.dataIOControl.outAtCellData || config.dataIOControl.outAtPointData);
434 DNDS_assert(config.dataIOControl.outPltVTKFormat || config.dataIOControl.outPltTecplotFormat || config.dataIOControl.outPltVTKHDFFormat);
435 DNDS_MAKE_SSP(outDistBnd, mpi);
436 outDistBnd->Resize(meshBnd->NumCell(), nOUTSBnd);
437
438 if (config.dataIOControl.outAtCellData)
439 {
440 DNDS_MAKE_SSP(outDist, mpi);
441 outDist->Resize(mesh->NumCell(), nOUTS);
442 }
443 if (config.dataIOControl.outAtPointData)
444 {
445 DNDS_MAKE_SSP(outDistPoint, mpi);
446 outDistPoint->Resize(mesh->NumNode(), nOUTSPoint);
447 DNDS_assert(nOUTSPoint >= nVars);
448
449 outDistPointPair.father = outDistPoint;
450 DNDS_MAKE_SSP(outDistPointPair.son, mpi);
451 outDistPointPair.TransAttach();
452 outDistPointPair.trans.BorrowGGIndexing(mesh->coords.trans);
453 outDistPointPair.trans.createMPITypes();
454 outDistPointPair.trans.initPersistentPull();
455 }
456
457 if (config.dataIOControl.outPltMode == 0)
458 {
459 //! serial mesh specific output method
460 DNDS_MAKE_SSP(outSerialBnd, mpi);
461 outDist2SerialTransBnd.setFatherSon(outDistBnd, outSerialBnd);
462 DNDS_assert(readerBnd->mode == Geom::MeshReaderMode::SerialOutput);
463 outDist2SerialTransBnd.BorrowGGIndexing(readerBnd->cell2nodeSerialOutTrans);
464 outDist2SerialTransBnd.createMPITypes();
465 outDist2SerialTransBnd.initPersistentPull();
466
467 if (config.dataIOControl.outAtCellData)
468 {
469 DNDS_MAKE_SSP(outSerial, mpi);
470 outDist2SerialTrans.setFatherSon(outDist, outSerial);
471 DNDS_assert(reader->mode == Geom::MeshReaderMode::SerialOutput);
472 outDist2SerialTrans.BorrowGGIndexing(reader->cell2nodeSerialOutTrans);
473 outDist2SerialTrans.createMPITypes();
474 outDist2SerialTrans.initPersistentPull();
475 }
476 if (config.dataIOControl.outAtPointData)
477 {
478 DNDS_MAKE_SSP(outSerialPoint, mpi);
479 outDist2SerialTransPoint.setFatherSon(outDistPoint, outSerialPoint);
480 DNDS_assert(reader->mode == Geom::MeshReaderMode::SerialOutput);
481 outDist2SerialTransPoint.BorrowGGIndexing(reader->coordSerialOutTrans);
482 outDist2SerialTransPoint.createMPITypes();
483 outDist2SerialTransPoint.initPersistentPull();
484 }
485 }
486 DNDS_MPI_InsertCheck(mpi, "ReadMeshAndInitialize -1 nvars " + std::to_string(nVars));
487 }
488
489 DNDS_SWITCH_INTELLISENSE(template <EulerModel model>, )
490 bool EulerSolver<model>::functor_fstop(int iter, ArrayDOFV<nVarsFixed> &cres, int iStep, RunningEnvironment &runningEnvironment)
491 {
492 using namespace std::literals;
493 DNDS_EULERSOLVER_RUNNINGENV_GET_REF_LIST
494 iterAll++;
495 // auto &uRecC = config.timeMarchControl.timeMarchIsTwoStage() && uPos == 1 ? uRec1 : uRec;
496
497 // auto &uRecC = config.timeMarchControl.timeMarchIsTwoStage() && uPos == 1 ? uRec1 : uRec;
498
499 bool useRHSasResBase = !config.timeMarchControl.steadyQuit && config.convergenceControl.resBaseType == 1;
500
501 Eigen::VectorFMTSafe<real, -1> res(nVars);
502 eval.EvaluateNorm(res, cres, config.convergenceControl.normOrd, config.convergenceControl.useVolWiseResidual);
503 if (config.convergenceControl.mergeMultiResidual == 1 && config.timeMarchControl.timeMarchIsTwoStage())
504 {
505 Eigen::VectorFMTSafe<real, -1> res1(nVars);
506 eval.EvaluateNorm(res1, ode->getRES(1), config.convergenceControl.normOrd, config.convergenceControl.useVolWiseResidual);
507 res += res1;
508 res *= 0.5;
509 }
510
511 Eigen::VectorFMTSafe<real, -1> resBaseNorm = res;
512
513 if (iter == 1)
514 resBaseCInternal.setZero();
515
516 if (useRHSasResBase)
517 {
518 if (iter == 1)
519 {
520 Eigen::VectorFMTSafe<real, -1> rhsBaseNorm(nVars);
521 eval.EvaluateNorm(rhsBaseNorm, ode->getLatestRHS(), 1, config.convergenceControl.useVolWiseResidual);
522 resBaseCInternal = rhsBaseNorm;
523 }
524 // {
525 // resBaseNorm = resBaseCInternal; // using Only RHS no res
526 // }
527 }
528
529 // if (iter == 1 && iStep == 1)
530 // * using 1st rk step for reference
531
532 resBaseCInternal = resBaseCInternal.array().max(resBaseNorm.array()); //! using max !
533
534 Eigen::VectorFMTSafe<real, -1> resRel = (res.array() / (resBaseCInternal.array() + verySmallReal)).matrix();
535 bool ifStop = resRel(0) < config.convergenceControl.rhsThresholdInternal; // ! using only rho's residual
536 if (config.convergenceControl.useCLDriver)
537 ifStop = ifStop && eval.pCLDriver->ConvergedAtTarget();
538 if (iter < config.convergenceControl.nTimeStepInternalMin)
539 ifStop = false;
540 auto [CLCur, CDCur, AOACur] = eval.CLDriverGetIntegrationUpdate(iter);
541 if (iter % config.outputControl.nConsoleCheckInternal == 0 || iter > config.convergenceControl.nTimeStepInternal || ifStop)
542 {
543 double tWall = MPI_Wtime();
544 real telapsed = MPI_Wtime() - tstartInternal;
545 bool useCollectiveTimer = config.outputControl.useCollectiveTimer;
546 real tcomm = Timer().getTimerColOrLoc(PerformanceTimer::Comm, mpi, useCollectiveTimer);
547 real tLimiterA = Timer().getTimerColOrLoc(PerformanceTimer::LimiterA, mpi, useCollectiveTimer);
548 real tLimiterB = Timer().getTimerColOrLoc(PerformanceTimer::LimiterB, mpi, useCollectiveTimer);
549 real trhs = Timer().getTimerColOrLoc(PerformanceTimer::RHS, mpi, useCollectiveTimer);
550 real trec = Timer().getTimerColOrLoc(PerformanceTimer::Reconstruction, mpi, useCollectiveTimer);
551 real tLim = Timer().getTimerColOrLoc(PerformanceTimer::Limiter, mpi, useCollectiveTimer);
552 real tPP = Timer().getTimerColOrLoc(PerformanceTimer::Positivity, mpi, useCollectiveTimer);
553 auto [telapsedM, telapsedS] = tInternalStats["t"].update(telapsed).get();
554 auto [tcommM, tcommS] = tInternalStats["c"].update(tcomm).get();
555 auto [trhsM, trhsS] = tInternalStats["r"].update(trhs).get();
556 auto [trecM, trecS] = tInternalStats["v"].update(trec).get();
557 auto [tLimM, tLimS] = tInternalStats["l"].update(tLim).get();
558 auto [tPPrM, tPPrS] = tInternalStats["p"].update(tPP).get();
559
560 if (mpi.rank == 0)
561 {
562 auto fmt = log().flags();
563 std::string formatStringMain = "";
564 for (auto &s : config.outputControl.consoleMainOutputFormatInternal)
565 formatStringMain += s;
566 log() << fmt::format(formatStringMain +
567 " "s +
568 (config.outputControl.consoleOutputMode == 1
569 ? "WallFlux {termYellow}{wallFlux:.6e}{termReset} CL,CD,AoA [{CLCur:.6e},{CDCur:.6e},{AOACur:.2e}]"s
570 : ""s),
571 DNDS_FMT_ARG(step),
572 DNDS_FMT_ARG(iStep),
573 DNDS_FMT_ARG(iter),
574 DNDS_FMT_ARG(iterAll),
575 fmt::arg("resRel", resRel.transpose()),
576 fmt::arg("wallFlux", eval.fluxWallSum.transpose()),
577 DNDS_FMT_ARG(tSimu),
578 DNDS_FMT_ARG(curDtImplicit),
579 DNDS_FMT_ARG(curDtMin),
580 DNDS_FMT_ARG(CFLNow),
581 DNDS_FMT_ARG(nLimInc),
582 DNDS_FMT_ARG(alphaMinInc),
583 DNDS_FMT_ARG(nLimBeta),
584 DNDS_FMT_ARG(minBeta),
585 DNDS_FMT_ARG(nLimAlpha),
586 DNDS_FMT_ARG(minAlpha),
587 DNDS_FMT_ARG(telapsed), DNDS_FMT_ARG(telapsedM),
588 DNDS_FMT_ARG(trec), DNDS_FMT_ARG(trecM),
589 DNDS_FMT_ARG(trhs), DNDS_FMT_ARG(trhsM),
590 DNDS_FMT_ARG(tcomm), DNDS_FMT_ARG(tcommM),
591 DNDS_FMT_ARG(tLim), DNDS_FMT_ARG(tLimM),
592 DNDS_FMT_ARG(tPP), DNDS_FMT_ARG(tPPrM),
593 DNDS_FMT_ARG(tLimiterA),
594 DNDS_FMT_ARG(tLimiterB),
595 DNDS_FMT_ARG(tWall),
596 DNDS_FMT_ARG(CLCur),
597 DNDS_FMT_ARG(CDCur),
598 DNDS_FMT_ARG(AOACur),
599 fmt::arg("termRed", TermColor::Red),
600 fmt::arg("termBlue", TermColor::Blue),
601 fmt::arg("termGreen", TermColor::Green),
602 fmt::arg("termCyan", TermColor::Cyan),
603 fmt::arg("termYellow", TermColor::Yellow),
604 fmt::arg("termBold", TermColor::Bold),
605 fmt::arg("termReset", TermColor::Reset));
606 log() << std::endl;
607 log().setf(fmt);
608
609 // std::string delimC = " ";
610 // logErr
611 // << std::left
612 // << step << delimC
613 // << std::left
614 // << iter << delimC
615 // << std::left
616 // << std::setprecision(config.outputControl.nPrecisionLog) << std::scientific
617 // << res.transpose() << delimC
618 // << tSimu << delimC
619 // << curDtMin << delimC
620 // << real(eval.nFaceReducedOrder) << delimC
621 // << eval.fluxWallSum.transpose() << delimC
622 // << (nLimInc) << delimC << (alphaMinInc) << delimC
623 // << (nLimBeta) << delimC << (minBeta) << delimC
624 // << (nLimAlpha) << delimC << (minAlpha) << delimC
625 // << std::endl;
626
627 // std::vector<std::string> logfileOutputTitles{
628 // "step", "iStep", "iterAll", "iter", "tSimu",
629 // "res", "curDtImplicit", "curDtMin", "CFLNow",
630 // "nLimInc", "alphaMinInc",
631 // "nLimBeta", "minBeta",
632 // "nLimAlpha", "minAlpha",
633 // "tWall", "telapsed", "trec", "trhs", "tcomm", "tLim", "tLimiterA", "tLimiterB",
634 // "fluxWall", "CL", "CD", "AoA"};
635 auto &fluxWall = eval.fluxWallSum;
636 auto &logErrVal = std::get<1>(logErr);
637#define DNDS_FILL_IN_LOG_ERR_VAL(v) FillLogValue(logErrVal, #v, v)
638 DNDS_FILL_IN_LOG_ERR_VAL(step);
639 DNDS_FILL_IN_LOG_ERR_VAL(iStep);
640 DNDS_FILL_IN_LOG_ERR_VAL(iterAll);
641 DNDS_FILL_IN_LOG_ERR_VAL(iter);
642 DNDS_FILL_IN_LOG_ERR_VAL(tSimu);
643 DNDS_FILL_IN_LOG_ERR_VAL(res);
644 DNDS_FILL_IN_LOG_ERR_VAL(curDtImplicit);
645 DNDS_FILL_IN_LOG_ERR_VAL(curDtMin);
646 DNDS_FILL_IN_LOG_ERR_VAL(CFLNow);
647
648 DNDS_FILL_IN_LOG_ERR_VAL(nLimInc);
649 DNDS_FILL_IN_LOG_ERR_VAL(alphaMinInc);
650 DNDS_FILL_IN_LOG_ERR_VAL(nLimBeta);
651 DNDS_FILL_IN_LOG_ERR_VAL(minBeta);
652 DNDS_FILL_IN_LOG_ERR_VAL(nLimAlpha);
653 DNDS_FILL_IN_LOG_ERR_VAL(minAlpha);
654
655 DNDS_FILL_IN_LOG_ERR_VAL(tWall);
656 DNDS_FILL_IN_LOG_ERR_VAL(telapsed);
657 DNDS_FILL_IN_LOG_ERR_VAL(trec);
658 DNDS_FILL_IN_LOG_ERR_VAL(trhs);
659 DNDS_FILL_IN_LOG_ERR_VAL(tcomm);
660 DNDS_FILL_IN_LOG_ERR_VAL(tLim);
661 DNDS_FILL_IN_LOG_ERR_VAL(tLimiterA);
662 DNDS_FILL_IN_LOG_ERR_VAL(tLimiterB);
663
664 DNDS_FILL_IN_LOG_ERR_VAL(fluxWall);
665 real CL{CLCur}, CD{CDCur}, AoA(AOACur);
666 DNDS_FILL_IN_LOG_ERR_VAL(CL);
667 DNDS_FILL_IN_LOG_ERR_VAL(CD);
668 DNDS_FILL_IN_LOG_ERR_VAL(AoA);
669#undef DNDS_FILL_IN_LOG_ERR_VAL
670
671 std::get<0>(logErr)->WriteLine(std::get<1>(logErr), config.outputControl.nPrecisionLog);
672
673 FillLogValue(logErrVal, "step", step);
674
675 eval.ConsoleOutputBndIntegrations();
676 eval.BndIntegrationLogWriteLine(
677 config.dataIOControl.getOutLogName() + "_" + output_stamp,
678 step, iStep, iter);
679 }
680 tstartInternal = MPI_Wtime();
681 Timer().clearAllTimer();
682 }
683
684 if (iter % config.outputControl.nDataOutInternal == 0)
685 {
686 eval.FixUMaxFilter(u);
687 PrintData(
688 config.dataIOControl.outPltName + "_" + output_stamp + "_" + std::to_string(step) + "_" + std::to_string(iter),
689 config.dataIOControl.outPltName + "_" + output_stamp + "_" + std::to_string(step), // internal series
690 [&](index iCell)
691 { return ode->getLatestRHS()[iCell](0); },
692 addOutList,
693 eval, tSimu);
694 eval.PrintBCProfiles(config.dataIOControl.outPltName + "_" + output_stamp + "_" + std::to_string(step),
695 u, uRec);
696 }
697 if ((iter % config.outputControl.nDataOutCInternal == 0) &&
698 !(config.outputControl.lazyCoverDataOutput && (iter % config.outputControl.nDataOutInternal == 0)))
699 {
700 eval.FixUMaxFilter(u);
701 PrintData(
702 config.dataIOControl.outPltName + "_" + output_stamp + "_" + "C",
703 "",
704 [&](index iCell)
705 { return ode->getLatestRHS()[iCell](0); },
706 addOutList,
707 eval, tSimu);
708 eval.PrintBCProfiles(config.dataIOControl.outPltName + "_" + output_stamp + "_" + "C",
709 u, uRec);
710 }
711 if (iter % config.outputControl.nRestartOutInternal == 0)
712 {
713 config.restartState.iStep = step;
714 config.restartState.iStepInternal = iter;
715 PrintRestart(config.dataIOControl.getOutRestartName() + "_" + output_stamp + "_" + std::to_string(step) + "_" + std::to_string(iter));
716 }
717 if ((iter % config.outputControl.nRestartOutCInternal == 0) &&
718 !(config.outputControl.lazyCoverDataOutput && (iter % config.outputControl.nRestartOutInternal == 0)))
719 {
720 config.restartState.iStep = step;
721 config.restartState.iStepInternal = iter;
722 PrintRestart(config.dataIOControl.getOutRestartName() + "_" + output_stamp + "_" + "C");
723 }
724 if (iter >= config.implicitCFLControl.nCFLRampStart && iter <= config.implicitCFLControl.nCFLRampLength + config.implicitCFLControl.nCFLRampStart)
725 {
726 real inter = real(iter - config.implicitCFLControl.nCFLRampStart) / config.implicitCFLControl.nCFLRampLength;
727 real logCFL = std::log(config.implicitCFLControl.CFL) + (std::log(config.implicitCFLControl.CFLRampEnd / config.implicitCFLControl.CFL) * inter);
728 CFLNow = std::exp(logCFL);
729 }
730 if (ifStop || iter > config.convergenceControl.nTimeStepInternal) //! TODO: reconstruct the framework of ODE-top-level-control
731 {
732 CFLNow = config.implicitCFLControl.CFL;
733 }
734 // return resRel.maxCoeff() < config.convergenceControl.rhsThresholdInternal;
735 return ifStop;
736 }
737
738 DNDS_SWITCH_INTELLISENSE(template <EulerModel model>, )
739 bool EulerSolver<model>::functor_fmainloop(RunningEnvironment &runningEnvironment)
740 {
741 using namespace std::literals;
742 DNDS_EULERSOLVER_RUNNINGENV_GET_REF_LIST
743 auto initUDOF = [&](ArrayDOFV<nVarsFixed> &uu)
744 { vfv->BuildUDof(uu, nVars); };
745 auto initUREC = [&](ArrayRECV<nVarsFixed> &uu)
746 { vfv->BuildURec(uu, nVars); };
747#define DNDS_EULER_SOLVER_GET_TEMP_UDOF(name) \
748 auto __p##name = uPool.getAllocInit(initUDOF); \
749 auto &name = *__p##name;
750
751 tSimu += curDtImplicit;
752 if (ifOutT)
753 tSimu = nextTout;
754 Eigen::VectorFMTSafe<real, -1> res(nVars);
755 eval.EvaluateNorm(res, ode->getLatestRHS(), 1, config.convergenceControl.useVolWiseResidual);
756 if (stepCount == 0 && resBaseC.norm() == 0)
757 resBaseC = res;
758
759 if (config.timeAverageControl.enabled)
760 {
761 DNDS_EULER_SOLVER_GET_TEMP_UDOF(uTemp)
762 eval.MeanValueCons2Prim(u, uTemp); // could use time-step-mean-u instead of latest-u
763 eval.TimeAverageAddition(uTemp, wAveraged, curDtImplicit, tAverage);
764 }
765
766 real CLCur{0.0}, CDCur{0.0}, AOACur{0.0};
767
768 if (step % config.outputControl.nConsoleCheck == 0)
769 {
770 double tWall = MPI_Wtime();
771 real telapsed = MPI_Wtime() - tstart;
772 bool useCollectiveTimer = config.outputControl.useCollectiveTimer;
773 real tcomm = Timer().getTimerColOrLoc(PerformanceTimer::Comm, mpi, useCollectiveTimer);
774 real tLimiterA = Timer().getTimerColOrLoc(PerformanceTimer::LimiterA, mpi, useCollectiveTimer);
775 real tLimiterB = Timer().getTimerColOrLoc(PerformanceTimer::LimiterB, mpi, useCollectiveTimer);
776 real trhs = Timer().getTimerColOrLoc(PerformanceTimer::RHS, mpi, useCollectiveTimer);
777 real trec = Timer().getTimerColOrLoc(PerformanceTimer::Reconstruction, mpi, useCollectiveTimer);
778 real tLim = Timer().getTimerColOrLoc(PerformanceTimer::Limiter, mpi, useCollectiveTimer);
779
780 tcomm = tInternalStats["c"].update(tcomm).getSum();
781 trhs = tInternalStats["r"].update(trhs).getSum();
782 trec = tInternalStats["v"].update(trec).getSum();
783 tLim = tInternalStats["l"].update(tLim).getSum();
784 auto tPPr = tInternalStats["p"].getSum() + Timer().getTimerColOrLoc(PerformanceTimer::PositivityOuter, mpi, useCollectiveTimer);
785 if (mpi.rank == 0)
786 {
787 auto format = log().flags();
788 std::string formatStringMain = "";
789 for (auto &s : config.outputControl.consoleMainOutputFormat)
790 formatStringMain += s;
791 log() << fmt::format(formatStringMain +
792 " "s +
793 (config.outputControl.consoleOutputMode == 1
794 ? "WallFlux {termYellow}{wallFlux:.6e}{termReset}"s
795 : ""s),
796 DNDS_FMT_ARG(step),
797 // DNDS_FMT_ARG(iStep),
798 // DNDS_FMT_ARG(iter),
799 fmt::arg("resRel", (res.array() / (resBaseC.array() + verySmallReal)).transpose()),
800 fmt::arg("wallFlux", eval.fluxWallSum.transpose()),
801 DNDS_FMT_ARG(tSimu),
802 DNDS_FMT_ARG(curDtImplicit),
803 DNDS_FMT_ARG(curDtMin),
804 DNDS_FMT_ARG(CFLNow),
805 DNDS_FMT_ARG(nLimInc),
806 DNDS_FMT_ARG(alphaMinInc),
807 DNDS_FMT_ARG(nLimBeta),
808 DNDS_FMT_ARG(minBeta),
809 DNDS_FMT_ARG(nLimAlpha),
810 DNDS_FMT_ARG(minAlpha),
811 DNDS_FMT_ARG(telapsed),
812 DNDS_FMT_ARG(trec),
813 DNDS_FMT_ARG(trhs),
814 DNDS_FMT_ARG(tcomm),
815 DNDS_FMT_ARG(tLim),
816 DNDS_FMT_ARG(tPPr),
817 DNDS_FMT_ARG(tLimiterA),
818 DNDS_FMT_ARG(tLimiterB),
819 DNDS_FMT_ARG(tWall),
820 fmt::arg("termRed", TermColor::Red),
821 fmt::arg("termBlue", TermColor::Blue),
822 fmt::arg("termGreen", TermColor::Green),
823 fmt::arg("termCyan", TermColor::Cyan),
824 fmt::arg("termYellow", TermColor::Yellow),
825 fmt::arg("termBold", TermColor::Bold),
826 fmt::arg("termReset", TermColor::Reset));
827 log() << std::endl;
828 log().setf(format);
829 auto &fluxWall = eval.fluxWallSum;
830 auto &logErrVal = std::get<1>(logErr);
831 int iStep{-1}, iter{-1};
832#define DNDS_FILL_IN_LOG_ERR_VAL(v) FillLogValue(logErrVal, #v, v)
833 DNDS_FILL_IN_LOG_ERR_VAL(step);
834 DNDS_FILL_IN_LOG_ERR_VAL(iStep);
835 DNDS_FILL_IN_LOG_ERR_VAL(iter);
836 DNDS_FILL_IN_LOG_ERR_VAL(tSimu);
837 DNDS_FILL_IN_LOG_ERR_VAL(res);
838 DNDS_FILL_IN_LOG_ERR_VAL(curDtImplicit);
839 DNDS_FILL_IN_LOG_ERR_VAL(curDtMin);
840 DNDS_FILL_IN_LOG_ERR_VAL(CFLNow);
841
842 DNDS_FILL_IN_LOG_ERR_VAL(nLimInc);
843 DNDS_FILL_IN_LOG_ERR_VAL(alphaMinInc);
844 DNDS_FILL_IN_LOG_ERR_VAL(nLimBeta);
845 DNDS_FILL_IN_LOG_ERR_VAL(minBeta);
846 DNDS_FILL_IN_LOG_ERR_VAL(nLimAlpha);
847 DNDS_FILL_IN_LOG_ERR_VAL(minAlpha);
848
849 DNDS_FILL_IN_LOG_ERR_VAL(tWall);
850 DNDS_FILL_IN_LOG_ERR_VAL(telapsed);
851 DNDS_FILL_IN_LOG_ERR_VAL(trec);
852 DNDS_FILL_IN_LOG_ERR_VAL(trhs);
853 DNDS_FILL_IN_LOG_ERR_VAL(tcomm);
854 DNDS_FILL_IN_LOG_ERR_VAL(tLim);
855 DNDS_FILL_IN_LOG_ERR_VAL(tLimiterA);
856 DNDS_FILL_IN_LOG_ERR_VAL(tLimiterB);
857
858 DNDS_FILL_IN_LOG_ERR_VAL(fluxWall);
859 real CL{CLCur}, CD{CDCur}, AoA(AOACur);
860 DNDS_FILL_IN_LOG_ERR_VAL(CL);
861 DNDS_FILL_IN_LOG_ERR_VAL(CD);
862 DNDS_FILL_IN_LOG_ERR_VAL(AoA);
863#undef DNDS_FILL_IN_LOG_ERR_VAL
864 std::get<0>(logErr)->WriteLine(std::get<1>(logErr), config.outputControl.nPrecisionLog);
865
866 eval.ConsoleOutputBndIntegrations();
867 eval.BndIntegrationLogWriteLine(
868 config.dataIOControl.getOutLogName() + "_" + output_stamp,
869 step, -1, -1);
870 }
871 tstart = MPI_Wtime();
872 Timer().clearAllTimer();
873 for (auto &s : tInternalStats)
874 s.second.clear();
875 }
876 if (step == nextStepOutC)
877 {
878 if (!(config.outputControl.lazyCoverDataOutput && (step == nextStepOut)))
879 {
880 eval.FixUMaxFilter(u);
881 PrintData(
882 config.dataIOControl.outPltName + "_" + output_stamp + "_" + "C",
883 "",
884 [&](index iCell)
885 { return ode->getLatestRHS()[iCell](0); },
886 addOutList,
887 eval, tSimu);
888 eval.PrintBCProfiles(config.dataIOControl.outPltName + "_" + output_stamp + "_" + "C",
889 u, uRec);
890 }
891 nextStepOutC += config.outputControl.nDataOutC;
892 }
893 if (step == nextStepOut)
894 {
895 eval.FixUMaxFilter(u);
896 PrintData(
897 config.dataIOControl.outPltName + "_" + output_stamp + "_" + std::to_string(step),
898 config.dataIOControl.outPltName + "_" + output_stamp, // physical ts series
899 [&](index iCell)
900 { return ode->getLatestRHS()[iCell](0); },
901 addOutList,
902 eval, tSimu);
903 eval.PrintBCProfiles(config.dataIOControl.outPltName + "_" + output_stamp + "_" + std::to_string(step),
904 u, uRec);
905 nextStepOut += config.outputControl.nDataOut;
906 }
907 if (step == nextStepOutAverageC)
908 {
909 if (!(config.outputControl.lazyCoverDataOutput && (step == nextStepOutAverage)))
910 {
911 DNDS_assert(config.timeAverageControl.enabled);
912 eval.MeanValuePrim2Cons(wAveraged, uAveraged);
913 eval.FixUMaxFilter(uAveraged);
914 PrintData(
915 config.dataIOControl.outPltName + "_TimeAveraged_" + output_stamp + "_" + "C",
916 "",
917 [&](index iCell)
918 { return ode->getLatestRHS()[iCell](0); },
919 addOutList,
920 eval, tSimu,
921 PrintDataTimeAverage);
922 }
923 nextStepOutAverageC += config.outputControl.nTimeAverageOutC;
924 }
925 if (step == nextStepOutAverage)
926 {
927 DNDS_assert(config.timeAverageControl.enabled);
928 eval.MeanValuePrim2Cons(wAveraged, uAveraged);
929 eval.FixUMaxFilter(uAveraged);
930 PrintData(
931 config.dataIOControl.outPltName + "_TimeAveraged_" + output_stamp + "_" + std::to_string(step),
932 config.dataIOControl.outPltName + "_TimeAveraged_" + output_stamp, // time average series
933 [&](index iCell)
934 { return ode->getLatestRHS()[iCell](0); },
935 addOutList,
936 eval, tSimu,
937 PrintDataTimeAverage);
938 nextStepOutAverage += config.outputControl.nTimeAverageOut;
939 }
940 if (step == nextStepRestartC)
941 {
942 if (!(config.outputControl.lazyCoverDataOutput && (step == nextStepRestart)))
943 {
944 config.restartState.iStep = step;
945 config.restartState.iStepInternal = -1;
946 PrintRestart(config.dataIOControl.getOutRestartName() + "_" + output_stamp + "_" + "C");
947 }
948 nextStepRestartC += config.outputControl.nRestartOutC;
949 }
950 if (step == nextStepRestart)
951 {
952 config.restartState.iStep = step;
953 config.restartState.iStepInternal = -1;
954 PrintRestart(config.dataIOControl.getOutRestartName() + "_" + output_stamp + "_" + std::to_string(step));
955 nextStepRestart += config.outputControl.nRestartOut;
956 }
957 if (ifOutT)
958 {
959 eval.FixUMaxFilter(u);
960 PrintData(
961 config.dataIOControl.outPltName + "_" + output_stamp + "_" + "t_" + std::to_string(nextTout),
962 config.dataIOControl.outPltName + "_" + output_stamp, // physical ts series
963 [&](index iCell)
964 { return ode->getLatestRHS()[iCell](0); },
965 addOutList,
966 eval, tSimu);
967 eval.PrintBCProfiles(config.dataIOControl.outPltName + "_" + output_stamp + "_" + "t_" + std::to_string(nextTout),
968 u, uRec);
969 nextTout += config.outputControl.tDataOut;
970 if (nextTout >= config.timeMarchControl.tEnd)
971 nextTout = config.timeMarchControl.tEnd;
972 }
973 if ((eval.settings.specialBuiltinInitializer == 2 ||
974 eval.settings.specialBuiltinInitializer == 203) &&
975 (step % config.outputControl.nConsoleCheck == 0)) // IV problem special: reduction on solution
976 {
977 auto FVal = [&](const Geom::tPoint &p, real t)
978 {
979 switch (eval.settings.specialBuiltinInitializer)
980 {
981 case 203:
982 return SpecialFields::IsentropicVortex10(eval, p, t, nVars, 10.0828);
983 default:
984 case 2:
985 return SpecialFields::IsentropicVortex10(eval, p, t, nVars, 5);
986 }
987 };
988 auto FWeight = [&](const Geom::tPoint &p, real t) -> real
989 {
990 // real xyOrig = t;
991 // real xCC = float_mod(p(0) - xyOrig, 10);
992 // real yCC = float_mod(p(1) - xyOrig, 10);
993 // return std::abs(xCC - 5.0) <= 2 && std::abs(yCC - 5.0) <= 2 ? 1.0 : 0.0;
994 return 1.0;
995 };
996 Eigen::Vector<real, -1> err1, errInf;
997 eval.EvaluateRecNorm(
998 err1, u, uRec, 1, true,
999 FVal, FWeight,
1000 tSimu);
1001 eval.EvaluateRecNorm(
1002 errInf, u, uRec, 1000, true,
1003 FVal, FWeight,
1004 tSimu);
1005
1006 if (mpi.rank == 0)
1007 {
1008 log() << "=== Mean Error IV: [" << std::scientific
1009 << std::setprecision(config.outputControl.nPrecisionConsole + 4) << err1(0) << ", "
1010 << err1(0) / vfv->GetGlobalVol() << ", "
1011 << errInf(0)
1012 << "]" << std::endl;
1013 }
1014 }
1015 if (config.implicitReconstructionControl.zeroGrads)
1016 uRec.setConstant(0.0), gradIsZero = true;
1017
1018 stepCount++;
1019
1020 return tSimu >= config.timeMarchControl.tEnd;
1021 }
1022}
DNDS_MAKE_SSP
#define DNDS_MAKE_SSP(ssp,...)
Definition Defines.hpp:217
DNDS_FMT_ARG
#define DNDS_FMT_ARG(V)
Definition Defines.hpp:43
DNDS_assert
#define DNDS_assert(expr)
Debug-only assertion (compiled out when DNDS_NDEBUG is defined). Prints the expression + file/line + ...
Definition Errors.hpp:112
EulerSolver.hpp
Top-level solver orchestrator for compressible Navier-Stokes / Euler simulations.
DNDS_EULERSOLVER_RUNNINGENV_GET_REF_LIST
#define DNDS_EULERSOLVER_RUNNINGENV_GET_REF_LIST
Definition EulerSolver.hpp:1328
DNDS_EULER_SOLVER_GET_TEMP_UDOF
#define DNDS_EULER_SOLVER_GET_TEMP_UDOF(name)
DNDS_FILL_IN_LOG_ERR_VAL
#define DNDS_FILL_IN_LOG_ERR_VAL(v)
DNDS_MPI_InsertCheck
#define DNDS_MPI_InsertCheck(mpi, info)
Debug helper: barrier + print "CHECK <info> RANK r @ fn @ file:line".
Definition MPI.hpp:343
DNDS_SWITCH_INTELLISENSE
#define DNDS_SWITCH_INTELLISENSE(real, intellisense)
Definition Macros.hpp:86
Mesh_Helpers.hpp
Unified mesh helpers: read, prepare, build boundary, serialize.
OMP.hpp
Helpers for OpenMP thread-count configuration at solver entry points.
SpecialFields.hpp
Analytic isentropic-vortex solutions for inviscid accuracy verification.
DNDS::Euler::ArrayDOFV
MPI-parallel distributed array of per-cell Degrees-of-Freedom (conservative variable vectors).
Definition Euler.hpp:93
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::settings
EulerEvaluatorSettings< model > settings
Physics and numerics settings for this evaluator.
Definition EulerEvaluator.hpp:220
DNDS::Euler::EulerSolver
Top-level solver orchestrator for compressible Navier-Stokes / Euler equations.
Definition EulerSolver.hpp:75
DNDS::PerformanceTimer::clearAllTimer
void clearAllTimer()
Zero every timer slot.
Definition Profiling.cpp:69
DNDS::PerformanceTimer::Positivity
@ Positivity
Positivity preservation.
Definition Profiling.hpp:45
DNDS::PerformanceTimer::Limiter
@ Limiter
Slope / variable limiter.
Definition Profiling.hpp:33
DNDS::PerformanceTimer::PositivityOuter
@ PositivityOuter
Outer-iteration positivity.
Definition Profiling.hpp:46
DNDS::PerformanceTimer::LimiterA
@ LimiterA
Limiter sub-phase A.
Definition Profiling.hpp:34
DNDS::PerformanceTimer::Comm
@ Comm
Catch-all MPI comm.
Definition Profiling.hpp:37
DNDS::PerformanceTimer::Reconstruction
@ Reconstruction
Variational reconstruction.
Definition Profiling.hpp:31
DNDS::PerformanceTimer::LimiterB
@ LimiterB
Limiter sub-phase B.
Definition Profiling.hpp:35
DNDS::PerformanceTimer::RHS
@ RHS
Total RHS evaluation.
Definition Profiling.hpp:29
DNDS::PerformanceTimer::getTimerColOrLoc
real getTimerColOrLoc(T t, const MPIInfo &mpi, bool col)
Either getTimerCollective (when col == true) or getTimer.
Definition Profiling.hpp:86
DNDS::Euler::SpecialFields::IsentropicVortex10
auto IsentropicVortex10(EulerEvaluator< model > &eval, const Geom::tPoint &x, real t, int cnVars, real chi)
Analytic isentropic vortex on a [0,10] x [0,10] periodic domain.
Definition SpecialFields.hpp:42
DNDS::Euler::NS_SA
@ NS_SA
NS + Spalart-Allmaras, 2D geometry (6 vars).
Definition Euler.hpp:877
DNDS::Euler::EulerBCType
EulerBCType
Boundary condition type identifiers for compressible flow solvers.
Definition EulerBC.hpp:36
DNDS::Euler::BCWallIsothermal
@ BCWallIsothermal
Isothermal wall; requires wall temperature in the BC value vector.
Definition EulerBC.hpp:41
DNDS::Euler::BCWall
@ BCWall
No-slip viscous wall boundary.
Definition EulerBC.hpp:39
DNDS::Euler::BCWallInvis
@ BCWallInvis
Inviscid slip wall boundary.
Definition EulerBC.hpp:40
DNDS::Euler::BCSym
@ BCSym
Symmetry plane boundary.
Definition EulerBC.hpp:46
DNDS::Geom::ReadMeshFromCGNS
void ReadMeshFromCGNS(ssp< UnstructuredMesh > &mesh, ssp< UnstructuredMeshSerialRW > &reader, const std::string &meshFile, const PartitionOptions &partOpts={}, real periodicTol=1e-9, int elevation=0, int bisect=0, std::function< t_index(const std::string &)> nameMapper=nullptr)
Read a CGNS mesh, partition, and optionally elevate/bisect.
Definition Mesh_Helpers.hpp:61
DNDS::Geom::BuildBndMesh
void BuildBndMesh(UnstructuredMesh &mesh, UnstructuredMesh &meshBnd, UnstructuredMeshSerialRW &readerBnd, bool buildSerialOut=true)
Extract the boundary surface mesh from a solver-ready volume mesh.
Definition Mesh_Helpers.hpp:248
DNDS::Geom::t_index
int32_t t_index
Definition Geometric.hpp:6
DNDS::Geom::SerialOutput
@ SerialOutput
Definition Mesh.hpp:1081
DNDS::Geom::SerializeMesh
void SerializeMesh(UnstructuredMesh &mesh, const std::string &outputPath, Serializer::SerializerFactory &factory)
Write a partitioned mesh to H5 for later distributed read.
Definition Mesh_Helpers.hpp:271
DNDS::Geom::tPoint
Eigen::Vector3d tPoint
Definition Geometric.hpp:9
DNDS::Geom::ReadMeshFromH5
void ReadMeshFromH5(UnstructuredMesh &mesh, Serializer::SerializerFactory factory, const std::string &h5Path, const PartitionOptions &partOpts={})
Read a mesh from DNDSR H5 with even-split + ParMetis repartition.
Definition Mesh_Helpers.hpp:135
DNDS::Geom::RotX
tGPoint RotX(real theta)
Definition Geometric.hpp:166
DNDS::Geom::RotZ
tGPoint RotZ(real theta)
Definition Geometric.hpp:158
DNDS::Geom::MeshH5Path
std::string MeshH5Path(const std::string &base, int mpiSize, int elevation=0, int bisect=0)
Build the conventional H5 filename for a partitioned mesh.
Definition Mesh_Helpers.hpp:296
DNDS::Geom::PrepareMesh
void PrepareMesh(UnstructuredMesh &mesh, UnstructuredMeshSerialRW &reader, const PrepareMeshOptions &opts={})
Prepare a distributed mesh for solver use.
Definition Mesh_Helpers.hpp:204
DNDS::Geom::ReadMeshFromH5Parallel
void ReadMeshFromH5Parallel(UnstructuredMesh &mesh, Serializer::SerializerFactory factory, const std::string &h5Path)
Read a pre-partitioned mesh from H5 (exact np match required).
Definition Mesh_Helpers.hpp:161
DNDS::Geom::RotY
tGPoint RotY(real theta)
Definition Geometric.hpp:174
DNDS::OMP::set_full_parallel_OMP
void set_full_parallel_OMP()
Configure OpenMP and Eigen threading for the "distributed OMP" pattern.
Definition OMP.hpp:19
DNDS::TermColor::Blue
constexpr std::string_view Blue
ANSI escape: bright blue foreground.
Definition Defines.hpp:880
DNDS::TermColor::Cyan
constexpr std::string_view Cyan
ANSI escape: bright cyan foreground.
Definition Defines.hpp:884
DNDS::TermColor::Reset
constexpr std::string_view Reset
ANSI escape: reset all attributes.
Definition Defines.hpp:888
DNDS::TermColor::Green
constexpr std::string_view Green
ANSI escape: bright green foreground.
Definition Defines.hpp:876
DNDS::TermColor::Yellow
constexpr std::string_view Yellow
ANSI escape: bright yellow foreground.
Definition Defines.hpp:878
DNDS::TermColor::Red
constexpr std::string_view Red
ANSI escape: bright red foreground.
Definition Defines.hpp:874
DNDS::TermColor::Bold
constexpr std::string_view Bold
ANSI escape: bold.
Definition Defines.hpp:890
DNDS::verySmallReal
DNDS_CONSTANT const real verySmallReal
Catch-all lower bound ("effectively zero").
Definition Defines.hpp:199
DNDS::pi
DNDS_CONSTANT const real pi
π in double precision (matches DNDS_E_PI macro).
Definition Defines.hpp:204
DNDS::index
int64_t index
Global row / DOF index type (signed 64-bit; handles multi-billion-cell meshes).
Definition Defines.hpp:112
DNDS::get_env_DNDS_DIST_OMP_NUM_THREADS
int get_env_DNDS_DIST_OMP_NUM_THREADS()
Read the DNDSR-specific DNDS_DIST_OMP_NUM_THREADS override, falling back to get_env_OMP_NUM_THREADS.
Definition Defines.cpp:142
DNDS::ssp
std::shared_ptr< T > ssp
Shortened alias for std::shared_ptr used pervasively in DNDSR.
Definition Defines.hpp:143
DNDS::real
double real
Canonical floating-point scalar used throughout DNDSR (double precision).
Definition Defines.hpp:110
DNDS::getTimeStamp
std::string getTimeStamp(const MPIInfo &mpi)
Format a human-readable timestamp using the calling rank as context.
Definition MPI.cpp:115
DNDS::log
std::ostream & log()
Return the current DNDSR log stream (either std::cout or the installed file).
Definition Defines.cpp:50
DNDS::Timer
PerformanceTimer & Timer()
Short-hand accessor to the PerformanceTimer singleton.
Definition Profiling.hpp:142
DNDS::Euler::EulerSolver::RunningEnvironment
Mutable state bundle for the time-marching loop.
Definition EulerSolver.hpp:1282
DNDS::Geom::PrepareMeshOptions
Options for PrepareMesh.
Definition Mesh_Helpers.hpp:187
DNDS::Geom::PrepareMeshOptions::buildSerialOut
bool buildSerialOut
build serial output arrays
Definition Mesh_Helpers.hpp:190
DNDS::Geom::PrepareMeshOptions::reorderCells
int reorderCells
0 = natural, 1 = reorder
Definition Mesh_Helpers.hpp:188
DNDS::Geom::PrepareMeshOptions::reorderParts
int reorderParts
nParts for ReorderLocalCells
Definition Mesh_Helpers.hpp:189
DNDS::Serializer::SerializerFactory
Config-backed factory selecting between JSON and HDF5 serializers.
Definition SerializerFactory.hpp:24
DistributedHex3D::cellElemInfo
tElemInfoArrayPair cellElemInfo
Definition SyntheticMeshBuilders.hpp:614
Eigen::MatrixFMTSafe
Eigen::Matrix wrapper that hides begin/end from fmt.
Definition EigenUtil.hpp:66
r
tVec r(NCells)
res
auto res
Definition test_MeshConnectivity_Interpolate.cpp:54
i
res i
Definition test_MeshConnectivity_Interpolate.cpp:205
mesh
DistributedHex3D mesh
Definition test_MeshConnectivity_Interpolate.cpp:929
nCell
const DNDS::index nCell
Definition test_MeshReorder.cpp:926
nNode
const DNDS::index nNode
Definition test_MeshReorder.cpp:927
ode
std::cout<< "[ESDIRK2-sc4] err1="<< std::scientific<< res.err_coarse<< " err2="<< res.err_fine<< " order="<< std::fixed<< res.order<< std::endl;CHECK(res.order > 1.8);CHECK(res.order< 2.5);}TEST_CASE("VBDF k=1: 1st-order on oscillator"){ ODE::ImplicitVBDFDualTimeStep< Vec2, Vec2 > ode(1, v2init, v2init, 1)
p
const tPoint const tPoint const tPoint & p
Definition test_Reconstruction.cpp:484