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