DNDSR 0.1.0.dev1+gcd065ad
Distributed Numeric Data Structure for CFV
Loading...
Searching...
No Matches
VariationalReconstruction_LimiterProcedure.hxx
Go to the documentation of this file.
1#pragma once
3#include "Limiters.hpp"
4
5namespace DNDS::CFV
6{
7 /**
8 * @brief Dispatches the biway limiter function selected by limiterBiwayAlter.
9 *
10 * Consolidates the switch(limiterBiwayAlter) block that appeared identically
11 * in both DoLimiterWBAP_C and DoLimiterWBAP_3.
12 */
13 template <int dim, int nVarsFixed, typename Tin1, typename Tin2, typename Tout>
14 inline void DispatchBiwayLimiter(int limiterBiwayAlter,
15 const Tin1 &uThis, const Tin2 &uOther,
16 Tout &uOut, real n)
17 {
18 switch (limiterBiwayAlter)
19 {
20 case 0:
21 FWBAP_L2_Biway(uThis, uOther, uOut, n);
22 break;
23 case 1:
24 FMINMOD_Biway(uThis, uOther, uOut, n);
25 break;
26 case 2:
27 FWBAP_L2_Biway_PolynomialNorm<dim, nVarsFixed>(uThis, uOther, uOut, n);
28 break;
29 case 3:
30 FMEMM_Biway_PolynomialNorm<dim, nVarsFixed>(uThis, uOther, uOut, n);
31 break;
32 case 4:
33 FWBAP_L2_Cut_Biway(uThis, uOther, uOut, n);
34 break;
35 default:
36 DNDS_assert_info(false, "no such limiterBiwayAlter code!");
37 }
38 }
39
40 template <int dim>
41 template <int nVarsFixed, int nVarsSee>
44 const std::array<int, nVarsSee> &varsSee)
45 {
46 using namespace Geom;
47 static const int maxNDiff = dim == 2 ? 10 : 20;
48
49#if defined(DNDS_DIST_MT_USE_OMP)
50# pragma omp parallel for schedule(runtime)
51#endif
52 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
53 {
54 // int NRecDOF = cellAtr[iCell].NDOF - 1; // ! not good ! TODO
55
56 auto c2f = mesh->cell2face[iCell];
57 Eigen::Matrix<real, nVarsSee, 2> IJIISIsum;
58 IJIISIsum.setZero(nVarsSee, 2);
59 for (int ic2f = 0; ic2f < c2f.size(); ic2f++)
60 {
62 index iCellOther = this->CellFaceOther(iCell, iFace);
63 auto gFace = this->GetFaceQuadO1(iFace);
64 decltype(IJIISIsum) IJIISI;
65 IJIISI.setZero(nVarsSee, 2);
66 gFace.IntegrationSimple(
67 IJIISI,
68 [&](auto &finc, int ig)
69 {
70 int nDiff = GetFaceAtr(iFace).NDIFF;
71 // int nDiff = 1;
72 tPoint unitNorm = faceMeanNorm[iFace];
73
74 Eigen::Matrix<real, Eigen::Dynamic, nVarsSee, Eigen::DontAlign, maxNDiff, nVarsSee>
76 uRecVal.setZero(), uRecValJump.setZero();
77 uRecValL = this->GetIntPointDiffBaseValue(iCell, iFace, -1, -1, Eigen::seq(0, nDiff - 1)) *
78 uRec[iCell](EigenAll, varsSee);
79 uRecValL(0, EigenAll) += u[iCell](varsSee).transpose();
80
82 {
83 uRecValR = this->GetIntPointDiffBaseValue(iCellOther, iFace, -1, -1, Eigen::seq(0, nDiff - 1)) *
84 uRec[iCellOther](EigenAll, varsSee);
85 uRecValR(0, EigenAll) += u[iCellOther](varsSee).transpose();
86 uRecVal = (uRecValL + uRecValR) * 0.5;
87 uRecValJump = (uRecValL - uRecValR) * 0.5;
88 }
89
90 Eigen::Matrix<real, nVarsSee, nVarsSee> IJI, ISI;
91 IJI = FFaceFunctional(uRecValJump, uRecValJump, iFace, iCell, iCellOther);
92 ISI = FFaceFunctional(uRecVal, uRecVal, iFace, iCell, iCellOther);
93
94 finc(EigenAll, 0) = IJI.diagonal();
95 finc(EigenAll, 1) = ISI.diagonal();
96
97 finc *= GetFaceArea(iFace); // don't forget this
98 });
100 }
101 Eigen::Vector<real, nVarsSee> smoothIndicator =
102 (IJIISIsum(EigenAll, 0).array() /
103 (IJIISIsum(EigenAll, 1).array() + verySmallReal))
104 .matrix();
105 real sImax = smoothIndicator.array().abs().maxCoeff();
106 si(iCell, 0) = std::sqrt(sImax) * sqr(settings.maxOrder);
107 }
108 }
109
110 template <int dim>
111 template <int nVarsFixed>
114 const Eigen::Vector<real, nVarsFixed> &varsSee,
116 {
117 using namespace Geom;
118 static const int maxNDiff = dim == 2 ? 10 : 20;
119 int nVars = u.RowSize();
120
121#if defined(DNDS_DIST_MT_USE_OMP)
122# pragma omp parallel for schedule(runtime)
123#endif
124 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
125 {
126 // int NRecDOF = cellAtr[iCell].NDOF - 1; // ! not good ! TODO
127
128 auto c2f = mesh->cell2face[iCell];
129 Eigen::Matrix<real, nVarsFixed, 2> IJIISIsum;
130 IJIISIsum.setZero(nVars, 2);
131 for (int ic2f = 0; ic2f < c2f.size(); ic2f++)
132 {
133 index iFace = c2f[ic2f];
134 index iCellOther = this->CellFaceOther(iCell, iFace);
135 auto gFace = this->GetFaceQuadO1(iFace);
136 decltype(IJIISIsum) IJIISI;
137 // if (iCellOther != UnInitIndex)
138 // {
139 // uRec[iCell].setConstant(1);
140 // uRec[iCellOther].setConstant(0);
141 // u[iCell].setConstant(1);
142 // u[iCellOther].setConstant(0);
143 // }
144 IJIISI.setZero(nVars, 2);
145 gFace.IntegrationSimple(
146 IJIISI,
147 [&](auto &finc, int ig)
148 {
149 tPoint unitNorm = faceMeanNorm[iFace];
150
151 Eigen::Matrix<real, 1, nVarsFixed>
152 uRecVal(1, nVarsFixed), uRecValL(1, nVarsFixed), uRecValR(1, nVarsFixed), uRecValJump(1, nVarsFixed);
153 uRecVal.setZero(1, nVars), uRecValJump.setZero(1, nVars);
154 uRecValL = this->GetIntPointDiffBaseValue(iCell, iFace, -1, -1, std::array<int, 1>{0}, 1) *
155 uRec[iCell];
156 uRecValL(0, EigenAll) += u[iCell].transpose();
158
159 if (iCellOther != UnInitIndex)
160 {
161 uRecValR = this->GetIntPointDiffBaseValue(iCellOther, iFace, -1, -1, std::array<int, 1>{0}, 1) *
162 uRec[iCellOther];
163 uRecValR(0, EigenAll) += u[iCellOther].transpose();
165 uRecVal = (uRecValL + uRecValR) * 0.5;
166 uRecValJump = (uRecValL - uRecValR) * 0.5;
167 }
168
169 for (int i = 0; i < nVars; i++)
170 {
171 finc(i, 0) = FFaceFunctional(uRecValJump(EigenAll, {i}), uRecValJump(EigenAll, {i}), iFace, iCell, iCellOther)(0, 0);
172 finc(i, 1) = FFaceFunctional(uRecVal(EigenAll, {i}), uRecVal(EigenAll, {i}), iFace, iCell, iCellOther)(0, 0);
173 }
174 finc *= GetFaceArea(iFace); // don't forget this
175 });
176 IJIISIsum += IJIISI;
177 }
178 Eigen::Vector<real, nVarsFixed> smoothIndicator =
179 (IJIISIsum(EigenAll, 0).array() /
180 (IJIISIsum(EigenAll, 1).array() + verySmallReal))
181 .matrix();
182 smoothIndicator.array() *= varsSee.array();
183 real sImax = smoothIndicator.array().abs().maxCoeff();
184 si(iCell, 0) = std::sqrt(sImax) * sqr(settings.maxOrder);
185 }
186 }
187
188 template <int dim>
189 template <int nVarsFixed>
192 tURec<nVarsFixed> &uRec,
193 tURec<nVarsFixed> &uRecNew,
196 bool ifAll,
198 bool putIntoNew)
199 {
200 using namespace Geom;
201
202#if defined(DNDS_DIST_MT_USE_OMP)
203# pragma omp parallel for schedule(runtime)
204#endif
205 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
206 {
207 if ((!ifAll) &&
208 si(iCell, 0) < settings.smoothThreshold)
209 {
210 uRecNew[iCell] = uRec[iCell]; //! no lim need to copy !!!!
211 continue;
212 }
213 index NRecDOF = GetCellAtr(iCell).NDOF - 1;
214 auto c2f = mesh->cell2face[iCell];
215 std::vector<Eigen::Matrix<real, Eigen::Dynamic, nVarsFixed, 0, maxRecDOF>> uFaces(c2f.size());
216 for (int ic2f = 0; ic2f < c2f.size(); ic2f++)
217 {
218 // * safety initialization
219 index iFace = c2f[ic2f];
220 index iCellOther = this->CellFaceOther(iCell, iFace);
221 if (iCellOther != UnInitIndex)
222 {
223 uFaces[ic2f].resizeLike(uRec[iCellOther]);
224 }
225 }
226
227 int cPOrder = settings.maxOrder;
228 for (; cPOrder >= 1; cPOrder--)
229 {
231
232 std::vector<Eigen::Array<real, Eigen::Dynamic, nVarsFixed, 0, maxRecDOFBatch>>
233 uOthers;
234 Eigen::Array<real, Eigen::Dynamic, nVarsFixed, 0, maxRecDOFBatch>
235 uC = uRec[iCell](
236 Eigen::seq(
237 LimStart,
238 LimEnd),
239 EigenAll);
240 uOthers.reserve(maxNeighbour);
241 uOthers.push_back(uC); // using uC centered
242 // DNDS_MPI_InsertCheck(mpi, "HereAAC");
243 for (int ic2f = 0; ic2f < c2f.size(); ic2f++)
244 {
245 index iFace = c2f[ic2f];
246 auto f2c = mesh->face2cell[iFace];
247 index iCellOther = this->CellFaceOther(iCell, iFace);
248 index iCellAtFace = f2c[0] == iCell ? 0 : 1;
249
250 if (iCellOther != UnInitIndex)
251 {
252 index NRecDOFOther = GetCellAtr(iCellOther).NDOF - 1;
254 if (NRecDOFLim < (LimEnd + 1))
255 continue; // reserved for p-adaption
256 // if (!(ifUseLimiter[iCell] & 0x0000000FU))
257 // continue;
258
259 tPoint unitNorm = faceMeanNorm[iFace];
260
261 const auto &matrixSecondaryThis =
262 this->GetMatrixSecondary(iCell, iFace, -1);
263
264 const auto &matrixSecondaryOther =
265 this->GetMatrixSecondary(iCellOther, iFace, -1);
266
267 Eigen::Matrix<real, Eigen::Dynamic, nVarsFixed, 0, maxRecDOF>
268 uOtherOther = uRec[iCellOther](Eigen::seq(0, NRecDOFLim - 1), EigenAll);
269
270 if (LimEnd < uOtherOther.rows() - 1) // successive SR
271 uOtherOther(Eigen::seq(LimEnd + 1, NRecDOFLim - 1), EigenAll) =
272 matrixSecondaryOther(Eigen::seq(LimEnd + 1, NRecDOFLim - 1), Eigen::seq(LimEnd + 1, NRecDOFLim - 1)) *
273 uFaces[ic2f](Eigen::seq(LimEnd + 1, NRecDOFLim - 1), EigenAll);
274
275 Eigen::Matrix<real, Eigen::Dynamic, nVarsFixed, 0, maxRecDOFBatch>
276 uOtherIn =
277 matrixSecondaryThis(Eigen::seq(LimStart, LimEnd), EigenAll) * uOtherOther;
278
279 Eigen::Matrix<real, Eigen::Dynamic, nVarsFixed, 0, maxRecDOFBatch>
280 uThisIn =
281 uC.matrix();
282
283 // 2 eig space :
284 auto uR = iCellAtFace ? u[iCell] : u[iCellOther];
285 auto uL = iCellAtFace ? u[iCellOther] : u[iCell];
286
288 uThisIn = (FM(uL, uR, unitNorm, uThisIn));
289
290 Eigen::Array<real, Eigen::Dynamic, nVarsFixed, 0, maxRecDOFBatch>
292
293 real n = settings.WBAP_nStd;
295 settings.limiterBiwayAlter,
296 uThisIn.array(), uOtherIn.array(), uLimOutArray, 1);
297
298 // to phys space
299 uLimOutArray = (FMI(uL, uR, unitNorm, uLimOutArray.matrix())).array();
300
301 uFaces[ic2f](Eigen::seq(LimStart, LimEnd), EigenAll) = uLimOutArray.matrix();
302 uOthers.push_back(uLimOutArray);
303 }
304 else
305 {
306 }
307 }
308 Eigen::Array<real, Eigen::Dynamic, nVarsFixed, 0, maxRecDOFBatch>
310
311 real n = settings.WBAP_nStd;
312 if (settings.normWBAP)
313 FWBAP_L2_Multiway_Polynomial2D(uOthers, uOthers.size(), uLimOutArray, n); // TODO: add 3D version here!
314 else
316
317 uRecNew[iCell](
318 Eigen::seq(
319 LimStart,
320 LimEnd),
321 EigenAll) = uLimOutArray.matrix();
322 }
323 }
324 uRecNew.trans.startPersistentPull();
325 uRecNew.trans.waitPersistentPull();
326 if (!putIntoNew)
327 {
328#if defined(DNDS_DIST_MT_USE_OMP)
329# pragma omp parallel for schedule(runtime)
330#endif
331 for (index iCell = 0; iCell < mesh->NumCellProc(); iCell++) // mind the edge
332 uRec[iCell] = uRecNew[iCell];
333 }
334 }
335
336 template <int dim>
337 template <int nVarsFixed>
340 tURec<nVarsFixed> &uRec,
341 tURec<nVarsFixed> &uRecNew,
344 bool ifAll,
346 bool putIntoNew)
347 {
348 using namespace Geom;
349
350 int cPOrder = settings.maxOrder;
351#if defined(DNDS_DIST_MT_USE_OMP)
352# pragma omp parallel for schedule(runtime)
353#endif
354 for (index iCell = 0; iCell < mesh->NumCellProc(); iCell++) // mind the edge
355 uRecNew[iCell] = uRec[iCell];
356 for (; cPOrder >= 1; cPOrder--)
357 {
359#if defined(DNDS_DIST_MT_USE_OMP)
360# pragma omp parallel for schedule(runtime)
361#endif
362 for (index iCell = 0; iCell < mesh->NumCellProc(); iCell++) // mind the edge
363 uRecBuf[iCell] = uRecNew[iCell];
364#if defined(DNDS_DIST_MT_USE_OMP)
365# pragma omp parallel for schedule(runtime)
366#endif
367 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
368 {
369 if ((!ifAll) &&
370 si(iCell, 0) < settings.smoothThreshold)
371 {
372 // uRecNew[iCell] = uRecBuf[iCell]; //! no copy for 3wbap!
373 continue;
374 }
375 index NRecDOF = GetCellAtr(iCell).NDOF - 1;
376 auto c2f = mesh->cell2face[iCell];
377 // std::vector<Eigen::Matrix<real, Eigen::Dynamic, nVarsFixed, 0, maxRecDOF>> uFaces(c2f.size());
378 for (int ic2f = 0; ic2f < c2f.size(); ic2f++)
379 {
380 // * safety initialization
381 index iFace = c2f[ic2f];
382 index iCellOther = this->CellFaceOther(iCell, iFace);
383 if (iCellOther != UnInitIndex)
384 {
385 // uFaces[ic2f].resizeLike(uRec[iCellOther]);
386 }
387 }
388
389 std::vector<Eigen::Array<real, Eigen::Dynamic, nVarsFixed, 0, maxRecDOFBatch>>
390 uOthers;
391 Eigen::Array<real, Eigen::Dynamic, nVarsFixed, 0, maxRecDOFBatch>
392 uC = uRecBuf[iCell](
393 Eigen::seq(
394 LimStart,
395 LimEnd),
396 EigenAll);
397 uOthers.reserve(maxNeighbour);
398 uOthers.push_back(uC); // using uC centered
399 // DNDS_MPI_InsertCheck(mpi, "HereAAC");
400 for (int ic2f = 0; ic2f < c2f.size(); ic2f++)
401 {
402 index iFace = c2f[ic2f];
403 auto f2c = mesh->face2cell[iFace];
404 index iCellOther = this->CellFaceOther(iCell, iFace);
405 index iCellAtFace = f2c[0] == iCell ? 0 : 1;
406
407 if (iCellOther != UnInitIndex)
408 {
409 index NRecDOFOther = GetCellAtr(iCellOther).NDOF - 1;
411 if (NRecDOFLim < (LimEnd + 1))
412 continue; // reserved for p-adaption
413 // if (!(ifUseLimiter[iCell] & 0x0000000FU))
414 // continue;
415
416 tPoint unitNorm = faceMeanNorm[iFace];
417
418 const auto &matrixSecondaryThis =
419 this->GetMatrixSecondary(iCell, iFace, -1);
420
421 const auto &matrixSecondaryOther =
422 this->GetMatrixSecondary(iCellOther, iFace, -1);
423
424 Eigen::Matrix<real, Eigen::Dynamic, nVarsFixed, 0, maxRecDOF>
425 uOtherOther = uRecBuf[iCellOther](Eigen::seq(0, NRecDOFLim - 1), EigenAll);
426
427 // if (LimEnd < uOtherOther.rows() - 1) // successive SR
428 // uOtherOther(Eigen::seq(LimEnd + 1, NRecDOFLim - 1), EigenAll) =
429 // matrixSecondaryOther(Eigen::seq(LimEnd + 1, NRecDOFLim - 1), Eigen::seq(LimEnd + 1, NRecDOFLim - 1)) *
430 // uFaces[ic2f](Eigen::seq(LimEnd + 1, NRecDOFLim - 1), EigenAll);
431
432 Eigen::Matrix<real, Eigen::Dynamic, nVarsFixed, 0, maxRecDOFBatch>
433 uOtherIn =
434 matrixSecondaryThis(Eigen::seq(LimStart, LimEnd), EigenAll) * uOtherOther;
435
436 Eigen::Matrix<real, Eigen::Dynamic, nVarsFixed, 0, maxRecDOFBatch>
437 uThisIn =
438 uC.matrix();
439
440 // 2 eig space :
441 auto uR = iCellAtFace ? u[iCell] : u[iCellOther];
442 auto uL = iCellAtFace ? u[iCellOther] : u[iCell];
443
446
447 Eigen::Array<real, Eigen::Dynamic, nVarsFixed, 0, maxRecDOFBatch>
449
450 real n = settings.WBAP_nStd;
451
453 settings.limiterBiwayAlter,
454 uThisIn.array(), uOtherIn.array(), uLimOutArray, 1);
455
456 // to phys space
457 uLimOutArray = FMI(uL, uR, unitNorm, uLimOutArray.matrix()).array();
458
459 // uFaces[ic2f](Eigen::seq(LimStart, LimEnd), EigenAll) = uLimOutArray.matrix();
460 uOthers.push_back(uLimOutArray);
461 }
462 else
463 {
464 }
465 }
466 Eigen::Array<real, Eigen::Dynamic, nVarsFixed, 0, maxRecDOFBatch>
468
469 real n = settings.WBAP_nStd;
470 if (settings.normWBAP)
471 FWBAP_L2_Multiway_Polynomial2D(uOthers, uOthers.size(), uLimOutArray, n);
472
473 else
475
476 uRecNew[iCell](
477 Eigen::seq(
478 LimStart,
479 LimEnd),
480 EigenAll) = uLimOutArray.matrix();
481 }
482 uRecNew.trans.startPersistentPull();
483 uRecNew.trans.waitPersistentPull();
484 }
485 if (!putIntoNew)
486 {
487#if defined(DNDS_DIST_MT_USE_OMP)
488# pragma omp parallel for schedule(runtime)
489#endif
490 for (index iCell = 0; iCell < mesh->NumCellProc(); iCell++) // mind the edge
491 uRec[iCell] = uRecNew[iCell];
492 }
493 }
494}
#define DNDS_assert_info(expr, info)
Debug-only assertion with an extra std::string info message.
Definition Errors.hpp:113
void DoLimiterWBAP_3(tUDof< nVarsFixed > &u, tURec< nVarsFixed > &uRec, tURec< nVarsFixed > &uRecNew, tURec< nVarsFixed > &uRecBuf, tScalarPair &si, bool ifAll, const tFMEig< nVarsFixed > &FM, const tFMEig< nVarsFixed > &FMI, bool putIntoNew=false)
FM(uLeft,uRight,norm) gives vsize * vsize mat of Left Eigen Vectors.
void DoLimiterWBAP_C(tUDof< nVarsFixed > &u, tURec< nVarsFixed > &uRec, tURec< nVarsFixed > &uRecNew, tURec< nVarsFixed > &uRecBuf, tScalarPair &si, bool ifAll, const tFMEig< nVarsFixed > &FM, const tFMEig< nVarsFixed > &FMI, bool putIntoNew=false)
FM(uLeft,uRight,norm) gives vsize * vsize mat of Left Eigen Vectors.
void DoCalculateSmoothIndicatorV1(tScalarPair &si, tURec< nVarsFixed > &uRec, tUDof< nVarsFixed > &u, const Eigen::Vector< real, nVarsFixed > &varsSee, const TFPost< nVarsFixed > &FPost)
void DoCalculateSmoothIndicator(tScalarPair &si, tURec< nVarsFixed > &uRec, tUDof< nVarsFixed > &u, const std::array< int, nVarsSee > &varsSee)
void DispatchBiwayLimiter(int limiterBiwayAlter, const Tin1 &uThis, const Tin2 &uOther, Tout &uOut, real n)
Dispatches the biway limiter function selected by limiterBiwayAlter.
void FWBAP_L2_Cut_Biway(const Tin1 &u1, const Tin2 &u2, Tout &uOut, real n)
Definition Limiters.hpp:483
void FWBAP_L2_Multiway(const TinOthers &uOthers, int Nother, Tout &uOut, real n1=1.0)
input vector<Eigen::Array-like>
Definition Limiters.hpp:361
void FMINMOD_Biway(const Tin1 &u1, const Tin2 &u2, Tout &uOut, real n)
input eigen arrays
Definition Limiters.hpp:503
void FWBAP_L2_Biway(const Tin1 &u1, const Tin2 &u2, Tout &uOut, real n)
input eigen arrays
Definition Limiters.hpp:432
DNDS_CONSTANT const index UnInitIndex
Sentinel "not initialised" index value (= INT64_MIN).
Definition Defines.hpp:176
DNDS_CONSTANT const real verySmallReal
Catch-all lower bound ("effectively zero").
Definition Defines.hpp:194
int64_t index
Global row / DOF index type (signed 64-bit; handles multi-billion-cell meshes).
Definition Defines.hpp:107
constexpr T sqr(const T &a)
a * a, constexpr. Works for all arithmetic types.
Definition Defines.hpp:517
double real
Canonical floating-point scalar used throughout DNDSR (double precision).
Definition Defines.hpp:105
Eigen::Vector3d n(1.0, 0.0, 0.0)