DNDSR 0.1.0.dev1+gcd065ad
Distributed Numeric Data Structure for CFV
Loading...
Searching...
No Matches
VariationalReconstruction_Reconstruction.hxx
Go to the documentation of this file.
1#pragma once
2
3#include "VariationalReconstruction.hpp"
4
5namespace DNDS
6{
7 namespace CFV
8 {
9
10 static const int dim = 3;
11 static const int nVarsFixed = 6;
12
13 DNDS_SWITCH_INTELLISENSE(
14 template <int dim>
15 template <int nVarsFixed>
16 ,
17 template <>
18 template <>)
19 auto VariationalReconstruction<dim>::GetBoundaryRHS(tURec<nVarsFixed> &uRec,
20 tUDof<nVarsFixed> &u,
21 index iCell, index iFace,
22 const TFBoundary<nVarsFixed> &FBoundary)
23 {
24 using namespace Geom;
25 using namespace Geom::Elem;
26 static const auto Seq012 = Eigen::seq(Eigen::fix<0>, Eigen::fix<dim - 1>);
27
28 Eigen::Matrix<real, Eigen::Dynamic, nVarsFixed> BCC;
29 BCC.setZero(uRec[iCell].rows(), uRec[iCell].cols());
30 auto faceID = mesh->GetFaceZone(iFace);
31 DNDS_assert(FaceIDIsExternalBC(faceID));
32
33 auto qFace = this->GetFaceQuad(iFace);
34 if (settings.intOrderVRBCIsSame() || settings.functionalSettings.greenGauss1Weight != 0)
35 qFace.IntegrationSimple(
36 BCC,
37 [&](auto &vInc, int iG)
38 {
39 RowVectorXR dbv =
40 this->GetIntPointDiffBaseValue(iCell, iFace, -1, iG, std::array<int, 1>{0}, 1);
41 Eigen::Vector<real, nVarsFixed> uBL =
42 (dbv *
43 uRec[iCell])
44 .transpose();
45 uBL += u[iCell]; //! need fixing?
46 Eigen::Vector<real, nVarsFixed> uBV =
47 FBoundary(
48 uBL,
49 u[iCell], iCell, iFace, iG,
50 this->GetFaceNorm(iFace, iG),
51 this->GetFaceQuadraturePPhysFromCell(iFace, iCell, -1, iG), faceID);
52 Eigen::RowVector<real, nVarsFixed> uIncBV = (uBV - u[iCell]).transpose();
53 if (settings.intOrderVRBCIsSame())
54 vInc = this->FFaceFunctional(dbv, uIncBV, iFace, iCell, iCell);
55 else
56 vInc.resizeLike(BCC), vInc.setZero();
57
58 if (settings.functionalSettings.greenGauss1Weight != 0)
59 {
60 // DNDS_assert(false); // not yet implemented
61 vInc += (settings.functionalSettings.greenGauss1Bias * this->GetGreenGauss1WeightOnCell(iCell) *
62 this->coefficients_.matrixAHalf_GG[iCell].transpose() * this->GetFaceNorm(iFace, iG)(Seq012)) *
63 uIncBV;
64 }
65 vInc *= this->GetFaceJacobiDet(iFace, iG);
66 // std::cout << baseWeight_.faceWeight[iFace].transpose() << std::endl;
67 });
68 if (!settings.intOrderVRBCIsSame())
69 {
70 auto qFace = Quadrature(mesh->GetFaceElement(iFace), settings.intOrderVRBCValue());
71 tSmallCoords coords;
72 mesh->GetCoordsOnFace(iFace, coords);
73 qFace.Integration(
74 BCC,
75 [&](auto &vInc, int __xxx_iG, const tPoint &pParam, const Elem::tD01Nj &DiNj) { // todo: cache these for bnd: pPhy JDet norm and dbv
76 BndVRPointCache &bndCacheEntry = baseWeight_.bndVRCaches.at(iFace).at(__xxx_iG);
77 auto &dbv = bndCacheEntry.D0Bj;
78 auto &np = bndCacheEntry.norm;
79 auto &JDet = bndCacheEntry.JDet;
80 auto &pPhy = bndCacheEntry.PPhy;
81
82 Eigen::Vector<real, nVarsFixed> uBL = (dbv * uRec[iCell]).transpose();
83 uBL += u[iCell];
84 Eigen::Vector<real, nVarsFixed> uBV =
85 FBoundary(
86 uBL,
87 u[iCell], iCell, iFace, -2,
88 np,
89 this->GetFacePointFromCell(iFace, iCell, -1, pPhy), faceID);
90 Eigen::RowVector<real, nVarsFixed> uIncBV = (uBV - u[iCell]).transpose();
91 vInc = this->FFaceFunctional(dbv, uIncBV, iFace, iCell, iCell);
92 vInc *= JDet;
93 });
94 }
95 return BCC;
96 }
97
98 template <int dim>
99 template <int nVarsFixed>
100 auto VariationalReconstruction<dim>::GetBoundaryRHSDiff(tURec<nVarsFixed> &uRec,
101 tURec<nVarsFixed> &uRecDiff,
102 tUDof<nVarsFixed> &u,
103 index iCell, index iFace,
104 const TFBoundaryDiff<nVarsFixed> &FBoundaryDiff)
105 {
106 using namespace Geom;
107 using namespace Geom::Elem;
108 static const auto Seq012 = Eigen::seq(Eigen::fix<0>, Eigen::fix<dim - 1>);
109
110 Eigen::Matrix<real, Eigen::Dynamic, nVarsFixed> BCC;
111 BCC.setZero(uRec[iCell].rows(), uRec[iCell].cols());
112 auto faceID = mesh->GetFaceZone(iFace);
113 DNDS_assert(FaceIDIsExternalBC(faceID));
114
115 auto qFace = this->GetFaceQuad(iFace);
116 if (settings.intOrderVRBCIsSame() || settings.functionalSettings.greenGauss1Weight != 0)
117 qFace.IntegrationSimple(
118 BCC,
119 [&](auto &vInc, int iG)
120 {
121 RowVectorXR dbv =
122 this->GetIntPointDiffBaseValue(iCell, iFace, -1, iG, std::array<int, 1>{0}, 1);
123 Eigen::Vector<real, nVarsFixed> uBL = (dbv * uRec[iCell]).transpose();
124 uBL += u[iCell]; //! need fixing?
125 Eigen::Vector<real, nVarsFixed> uBLDiff = (dbv * uRecDiff[iCell]).transpose();
126 Eigen::Vector<real, nVarsFixed>
127 uBV = FBoundaryDiff(
128 uBL, uBLDiff,
129 u[iCell], iCell, iFace, iG,
130 this->GetFaceNorm(iFace, iG),
131 this->GetFaceQuadraturePPhysFromCell(iFace, iCell, -1, iG), faceID);
132 Eigen::RowVector<real, nVarsFixed> uIncBV = uBV.transpose();
133 if (settings.intOrderVRBCIsSame())
134 vInc = this->FFaceFunctional(dbv, uIncBV, iFace, iCell, iCell);
135 else
136 vInc.resizeLike(BCC), vInc.setZero();
137 // std::cout << baseWeight_.faceWeight[iFace].transpose() << std::endl;
138 if (settings.functionalSettings.greenGauss1Weight != 0)
139 {
140 // DNDS_assert(false); // not yet implemented
141 vInc += (settings.functionalSettings.greenGauss1Bias * this->GetGreenGauss1WeightOnCell(iCell) *
142 this->coefficients_.matrixAHalf_GG[iCell].transpose() * this->GetFaceNorm(iFace, iG)(Seq012)) *
143 uIncBV;
144 }
145 vInc *= this->GetFaceJacobiDet(iFace, iG);
146 });
147 if (!settings.intOrderVRBCIsSame())
148 {
149 auto qFace = Quadrature(mesh->GetFaceElement(iFace), settings.intOrderVRBCValue());
150 tSmallCoords coords;
151 mesh->GetCoordsOnFace(iFace, coords);
152 qFace.Integration(
153 BCC,
154 [&](auto &vInc, int __xxx_iG, const tPoint &pParam, const Elem::tD01Nj &DiNj)
155 {
156 BndVRPointCache &bndCacheEntry = baseWeight_.bndVRCaches.at(iFace).at(__xxx_iG);
157 auto &dbv = bndCacheEntry.D0Bj;
158 auto &np = bndCacheEntry.norm;
159 auto &JDet = bndCacheEntry.JDet;
160 auto &pPhy = bndCacheEntry.PPhy;
161 Eigen::Vector<real, nVarsFixed> uBL = (dbv * uRec[iCell]).transpose();
162 uBL += u[iCell];
163 Eigen::Vector<real, nVarsFixed> uBLDiff = (dbv * uRecDiff[iCell]).transpose();
164 Eigen::Vector<real, nVarsFixed> uBV =
165 FBoundaryDiff(
166 uBL, uBLDiff,
167 u[iCell], iCell, iFace, -2,
168 np,
169 this->GetFacePointFromCell(iFace, iCell, -1, pPhy), faceID);
170 Eigen::RowVector<real, nVarsFixed> uIncBV = uBV.transpose();
171 vInc = this->FFaceFunctional(dbv, uIncBV, iFace, iCell, iCell);
172 vInc *= JDet;
173 });
174 }
175 return BCC;
176 }
177
178 DNDS_SWITCH_INTELLISENSE(
179 template <int dim>
180 template <int nVarsFixed>
181 ,
182 template <>
183 template <>)
184 void VariationalReconstruction<dim>::DoReconstruction2ndGrad(
185 tUGrad<nVarsFixed, dim> &uRec,
186 tUDof<nVarsFixed> &u,
187 const TFBoundary<nVarsFixed> &FBoundary,
188 int method)
189 {
190 using namespace Geom;
191 using namespace Geom::Elem;
192
193 static const auto Seq012 = Eigen::seq(Eigen::fix<0>, Eigen::fix<dim - 1>); // note this is gDim!
194 static const auto Seq123 = Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>);
195
196 if (method == 1 || method == 11)
197 {
198 // simple gauss rule
199#if defined(DNDS_DIST_MT_USE_OMP)
200# pragma omp parallel for schedule(runtime)
201#endif
202 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
203 {
204 int nVars = u[iCell].size();
205 auto c2f = mesh->cell2face[iCell];
206 Eigen::Matrix<real, nVarsFixed, dim> grad;
207 grad.setZero(nVars, dim);
208 MatrixXR mGG;
209 Eigen::Matrix<real, Eigen::Dynamic, nVarsFixed> bGG;
210 if (method == 11)
211 {
212 mGG.setZero(dim + c2f.size(), dim + c2f.size());
213 bGG.setZero(dim + c2f.size(), u[0].rows());
214 mGG(Seq012, Seq012).setIdentity(); // for GGMP
215 // mGG(Seq012, Seq012) *= this->GetCellVol(iCell);
216 }
217
218 for (int ic2f = 0; ic2f < c2f.size(); ic2f++)
219 {
220 index iFace = c2f[ic2f];
221 index iCellOther = CellFaceOther(iCell, iFace);
222 auto faceID = mesh->GetFaceZone(iFace);
223 int if2c = CellIsFaceBack(iCell, iFace) ? 0 : 1;
224 auto faceCent = this->GetFaceQuadraturePPhysFromCell(iFace, iCell, if2c, -1);
225 real lThis = (faceCent - this->GetCellBary(iCell)).norm();
226 if (settings.subs2ndOrderGGScheme == 0)
227 lThis = 1;
228 real lThat = lThis;
229 if (iCellOther != UnInitIndex)
230 {
231 lThat = (this->GetOtherCellBaryFromCell(iCell, iCellOther, iFace) - faceCent).norm();
232 if (settings.subs2ndOrderGGScheme == 0)
233 lThat = 1;
234 Eigen::RowVector<real, nVarsFixed> uOther = u[iCellOther].transpose();
235 this->ApplyPeriodicTransform(if2c, faceID, uOther);
236 Eigen::Vector<real, dim> uNorm = this->GetFaceNormFromCell(iFace, iCell, -1, -1)(Seq012) * (CellIsFaceBack(iCell, iFace) ? 1. : -1.);
237 Eigen::Matrix<real, nVarsFixed, dim> gradInc =
238 (uOther.transpose() - u[iCell]) *
239 lThis / (lThis + lThat + verySmallReal) * this->GetFaceArea(iFace) * uNorm.transpose();
240 grad += gradInc;
241 if (method == 11)
242 {
243 mGG(Seq012, dim + ic2f) = -this->GetFaceArea(iFace) / (this->GetCellVol(iCell) + verySmallReal) * uNorm;
244 mGG(dim + ic2f, Seq012) =
245 (this->GetCellBary(iCell)(Seq012) + this->GetOtherCellBaryFromCell(iCell, iCellOther, iFace)(Seq012) -
246 2 * faceCent(Seq012))
247 .transpose();
248 mGG(dim + ic2f, dim + ic2f) = 2;
249 bGG(dim + ic2f, EigenAll) = uOther + u[iCell].transpose();
250 }
251 }
252 else
253 {
254 auto faceID = mesh->GetFaceZone(iFace);
255 DNDS_assert(FaceIDIsExternalBC(faceID));
256
257 Eigen::Vector<real, nVarsFixed> uBL = u[iCell]; //! need fixing?
258 Eigen::Vector<real, nVarsFixed> uBV =
259 FBoundary(
260 uBL,
261 u[iCell], iCell, iFace, -1,
262 this->GetFaceNorm(iFace, -1),
263 faceCent, faceID);
264 Eigen::Vector<real, dim> uNorm = this->GetFaceNorm(iFace, -1)(Seq012);
265 grad += (uBV - u[iCell]) * 0.5 * this->GetFaceArea(iFace) * uNorm.transpose();
266 if (method == 11)
267 {
268 mGG(Seq012, dim + ic2f) = -this->GetFaceArea(iFace) / (this->GetCellVol(iCell) + verySmallReal) * uNorm;
269 // Eigen::Vector<real, dim> BaryOther = this->GetCellBary(iCell)(Seq012) +
270 // 2 * uNorm * uNorm.dot(faceCent(Seq012) - this->GetCellBary(iCell)(Seq012));
271 Eigen::Vector<real, dim> BaryOther = 2 * faceCent(Seq012) - this->GetCellBary(iCell)(Seq012);
272 mGG(dim + ic2f, Seq012) = (this->GetCellBary(iCell)(Seq012) + BaryOther -
273 2 * faceCent(Seq012))
274 .transpose();
275 mGG(dim + ic2f, dim + ic2f) = 2;
276 bGG(dim + ic2f, EigenAll) = (uBV + u[iCell]).transpose();
277 }
278 }
279 }
280
281 grad /= GetCellVol(iCell);
282
283 if (method == 11)
284 {
285
286 // auto mGGLU = mGG.colPivHouseholderQr();
287 auto mGGLU = mGG.fullPivLu();
288 DNDS_assert(mGGLU.isInvertible());
289 Eigen::Matrix<real, Eigen::Dynamic, nVarsFixed> xGG = mGGLU.solve(bGG);
290 grad = xGG(Seq012, EigenAll).transpose();
291
292 if (mesh->GetCellElement(iCell).type == Elem::Pyramid5)
293 {
294 std::cout << "pyramid5\n";
295 std::cout << mGG << std::endl;
296 std::cout << bGG << std::endl;
297 auto svd = mGG.bdcSvd();
298 std::cout << svd.singularValues().transpose() << std::endl;
299 std::cout << svd.singularValues().maxCoeff() / svd.singularValues().minCoeff() << std::endl;
300 }
301 if (mesh->GetCellElement(iCell).type == Elem::Tet4)
302 {
303 std::cout << "Tet4\n";
304 std::cout << mGG << std::endl;
305 std::cout << bGG << std::endl;
306 auto svd = mGG.bdcSvd();
307 std::cout << svd.singularValues().transpose() << std::endl;
308 std::cout << svd.singularValues().maxCoeff() / svd.singularValues().minCoeff() << std::endl;
309 }
310 }
311
312 // tPoint cellBary = GetCellBary(iCell);
313 // Eigen::MatrixXd vvv = (grad * (Geom::RotZ(90) * cellBary)(Eigen::seq(0, dim - 1))).transpose();
314 // std::cout << cellBary.transpose() << " ---- " << vvv << std::endl;
315 // DNDS_assert(vvv(0) < 1e-10);
316
317 uRec[iCell] = grad.transpose();
318 }
319 }
320 else if (method == 2) //! warning, periodic not implemented here
321 {
322#if defined(DNDS_DIST_MT_USE_OMP)
323# pragma omp parallel for schedule(runtime)
324#endif
325 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
326 {
327 int nVars = u[iCell].size();
328 auto c2f = mesh->cell2face[iCell];
329 Eigen::Matrix<real, nVarsFixed, dim> grad;
330 grad.setZero(nVars, dim);
331 Eigen::Matrix<real, Eigen::Dynamic, dim> dcs;
332 Eigen::Matrix<real, Eigen::Dynamic, nVarsFixed> dus;
333 dcs.resize(c2f.size(), dim);
334 dus.resize(c2f.size(), nVars);
335
336 for (int ic2f = 0; ic2f < c2f.size(); ic2f++)
337 {
338 index iFace = c2f[ic2f];
339 index iCellOther = CellFaceOther(iCell, iFace);
340
341 if (iCellOther != UnInitIndex)
342 {
343 dcs(ic2f, EigenAll) = (GetOtherCellBaryFromCell(iCell, iCellOther, iFace) - GetCellBary(iCell))(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim - 1>))
344 .transpose();
345 dus(ic2f, EigenAll) = (u[iCellOther] - u[iCell]).transpose();
346 }
347 else
348 {
349 auto faceID = mesh->GetFaceZone(iFace);
350 DNDS_assert(FaceIDIsExternalBC(faceID));
351
352 Eigen::Vector<real, nVarsFixed> uBL;
353 uBL.setZero();
354 uBL += u[iCell]; //! need fixing?
355 Eigen::Vector<real, nVarsFixed> uBV =
356 FBoundary(
357 uBL,
358 u[iCell], iCell, iFace, -1,
359 this->GetFaceNorm(iFace, -1),
360 this->GetFaceQuadraturePPhysFromCell(iFace, iCell, -1, -1), faceID);
361 dus(ic2f, EigenAll) = (uBV - u[iCell]).transpose();
362
363 dcs(ic2f, EigenAll) =
364 ((GetCellBary(iCell) - GetFaceQuadraturePPhys(iFace, -1)).dot(GetFaceNorm(iFace, -1)) * (-2.) * GetFaceNorm(iFace, -1))(
365 Eigen::seq(Eigen::fix<0>, Eigen::fix<dim - 1>))
366 .transpose();
367 }
368 }
369 // Eigen::MatrixXd m;
370 // m.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve()
371 auto svd = dcs.template bdcSvd<Eigen::ComputeFullU | Eigen::ComputeFullV>();
372 grad = svd.solve(dus).transpose();
373 uRec[iCell] = grad.transpose();
374 }
375 }
376 else
377 {
378 DNDS_assert_info(false, "NO SUCH 2nd rec METHOD");
379 }
380 }
381
382 DNDS_SWITCH_INTELLISENSE(
383 template <int dim>
384 template <int nVarsFixed>
385 ,
386 template <>
387 template <>)
388 void VariationalReconstruction<dim>::DoReconstruction2nd(
389 tURec<nVarsFixed> &uRec,
390 tUDof<nVarsFixed> &u,
391 const TFBoundary<nVarsFixed> &FBoundary,
392 int method,
393 const std::vector<int> &mask)
394 {
395 using namespace Geom;
396 using namespace Geom::Elem;
397 // using TU = Eigen::Vector<real, nVarsFixed>;
398 // using TU_Batch = Eigen::Matrix<real, nVarsFixed, Eigen::Dynamic>;
399 int nVars = u.father->MatRowSize();
400 auto vfv = this;
401
402 tUGrad<nVarsFixed, dim> uGrad;
403 this->BuildUGrad(uGrad, u.father->MatRowSize(), false, false); // no son or trans makes this fully local operation
404
405 static const auto Seq012 = Eigen::seq(Eigen::fix<0>, Eigen::fix<dim - 1>); // note this is gDim!
406 static const auto Seq123 = Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>);
407
408 this->DoReconstruction2ndGrad(uGrad, u, FBoundary, method);
409
410 // // in place barth limiting
411 // for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
412 // {
413 // auto c2f = mesh->cell2face[iCell];
414
415 // TU_Batch uFaceInc;
416 // uFaceInc.setZero(nVars, c2f.size() * 2); // j < c2f.size(): faceInc; j > c2f.size(): baryInc
417 // TU uOtherMin = u[iCell];
418 // TU uOtherMax = u[iCell];
419 // for (rowsize ic2f = 0; ic2f < c2f.size(); ic2f++)
420 // {
421 // index iFace = c2f[ic2f];
422 // uFaceInc(EigenAll, ic2f) =
423 // uGrad[iCell].transpose() *
424 // (vfv->GetFaceQuadraturePPhysFromCell(iFace, iCell, -1, -1) - vfv->GetCellQuadraturePPhys(iCell, -1))(Seq012);
425 // index iCellOther = mesh->CellFaceOther(iCell, iFace);
426 // if (iCellOther != UnInitIndex)
427 // {
428 // uOtherMin = uOtherMin.array().min(u[iCellOther].array());
429 // uOtherMax = uOtherMin.array().max(u[iCellOther].array());
430 // }
431 // }
432
433 // TU uFaceIncMax = uFaceInc.array().rowwise().maxCoeff();
434 // TU uFaceIncMin = uFaceInc.array().rowwise().minCoeff();
435 // TU alpha0;
436 // alpha0.setConstant(nVars, 1.0);
437 // alpha0 = alpha0.array().min(((uOtherMax - u[iCell]).array().abs() / (uFaceIncMax.array().abs() + verySmallReal)));
438 // alpha0 = alpha0.array().min(((uOtherMin - u[iCell]).array().abs() / (uFaceIncMin.array().abs() + verySmallReal)));
439 // uGrad[iCell].array().rowwise() *= alpha0.array().transpose();
440
441 // // kill the axis-cells
442 // for (rowsize ic2f = 0; ic2f < c2f.size(); ic2f++)
443 // {
444 // index iFace = c2f[ic2f];
445 // if (this->axisFaces.count(iFace))
446 // uGrad[iCell] *= 0;
447 // }
448 // }
449
450#if defined(DNDS_DIST_MT_USE_OMP)
451# pragma omp parallel for schedule(static)
452#endif
453 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
454
455 // #if defined(DNDS_DIST_MT_USE_OMP)
456 // #pragma omp parallel for schedule(static)
457 // #endif
458 // for (int iPart = 0; iPart < mesh->NLocalParts(); iPart++)
459 // for (index iCell = mesh->LocalPartStart(iPart); iCell < mesh->LocalPartEnd(iPart); iCell++)
460 {
461 int nVars = u[iCell].size();
462 auto c2f = mesh->cell2face[iCell];
463 Eigen::Matrix<real, nVarsFixed, dim> grad = uGrad[iCell].transpose();
464
465 Eigen::Matrix<real, dim, dim> d1bv;
466 d1bv = this->GetIntPointDiffBaseValue(
467 iCell, -1, -1, -1, Seq123, dim + 1)(EigenAll, Seq012);
468 auto lud1bv = d1bv.partialPivLu();
469 if (lud1bv.rcond() > 1e9)
470 std::cout << "Large Cond " << lud1bv.rcond() << std::endl;
471 if (mask.size() == 0)
472 {
473 uRec[iCell].setZero();
474 uRec[iCell](Seq012, EigenAll) = lud1bv.solve(grad.transpose());
475 }
476 else
477 {
478 uRec[iCell](EigenAll, mask).setZero();
479 uRec[iCell](Seq012, mask) = lud1bv.solve(grad.transpose())(EigenAll, mask);
480 }
481
482 }
483 }
484
485 template <int dim>
486 template <int nVarsFixed>
487 void VariationalReconstruction<dim>::DoReconstructionIter(
488 tURec<nVarsFixed> &uRec,
489 tURec<nVarsFixed> &uRecNew,
490 tUDof<nVarsFixed> &u,
491 const TFBoundary<nVarsFixed> &FBoundary,
492 bool putIntoNew,
493 bool recordInc,
494 bool uRecIsZero)
495 {
496 using namespace Geom;
497 using namespace Geom::Elem;
498 using namespace Geom::Base;
499 static const auto Seq012 = Eigen::seq(Eigen::fix<0>, Eigen::fix<dim - 1>);
500 int maxNDOF = GetNDof<dim>(settings.maxOrder);
501 if (recordInc)
502 DNDS_assert_info(putIntoNew, "the -RHS must be put into uRecNew");
503 if (settings.maxOrder == 1 && settings.subs2ndOrder != 0)
504 {
505 if (recordInc)
506 this->DoReconstruction2nd<nVarsFixed>(uRecNew, u, (FBoundary), settings.subs2ndOrder, std::vector<int>());
507 else if (putIntoNew)
508 this->DoReconstruction2nd<nVarsFixed>(uRecNew, u, (FBoundary), settings.subs2ndOrder, std::vector<int>());
509 else
510 this->DoReconstruction2nd<nVarsFixed>(uRec, u, (FBoundary), settings.subs2ndOrder, std::vector<int>());
511 return;
512 }
513
514 // #if defined(DNDS_DIST_MT_USE_OMP)
515 // #pragma omp parallel for schedule(static)
516 // #endif
517 // for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
518
519#if defined(DNDS_DIST_MT_USE_OMP)
520# pragma omp parallel for schedule(static)
521#endif
522 for (int iPart = 0; iPart < mesh->NLocalParts(); iPart++)
523 for (index iCell = mesh->LocalPartStart(iPart); iCell < mesh->LocalPartEnd(iPart); iCell++)
524 {
525 real relax = GetCellAtr(iCell).relax;
526
527 if (recordInc)
528 {
529 if (uRecIsZero)
530 uRecNew[iCell].setZero();
531 else
532 uRecNew[iCell] = uRec[iCell];
533 }
534 else if (settings.SORInstead)
535 uRec[iCell] = uRec[iCell] * ((recordInc ? 0 : 1) - relax);
536 else
537 uRecNew[iCell] = uRec[iCell] * ((recordInc ? 0 : 1) - relax);
538
539 auto c2f = mesh->cell2face[iCell];
540 auto matrixAAInvBRow = coefficients_.matrixAAInvB[iCell];
541 auto vectorAInvBRow = coefficients_.vectorAInvB[iCell];
542 for (int ic2f = 0; ic2f < c2f.size(); ic2f++)
543 {
544 index iFace = c2f[ic2f];
545 index iCellOther = CellFaceOther(iCell, iFace);
546 auto faceID = mesh->GetFaceZone(iFace);
547 int if2c = CellIsFaceBack(iCell, iFace) ? 0 : 1;
548 if (iCellOther != UnInitIndex)
549 {
550 Eigen::RowVector<real, nVarsFixed> uOther = u[iCellOther];
551 Eigen::Matrix<real, Eigen::Dynamic, nVarsFixed> uRecOther = uRec[iCellOther];
552 this->ApplyPeriodicTransform(if2c, faceID, uOther, uRecOther);
553 if (recordInc)
554 {
555 if (uRecIsZero)
556 uRecNew[iCell] -=
557 (vectorAInvBRow[ic2f] * (uOther - u[iCell].transpose()));
558 else
559 uRecNew[iCell] -=
560 (matrixAAInvBRow[ic2f + 1] * uRecOther +
561 vectorAInvBRow[ic2f] * (uOther - u[iCell].transpose()));
562 }
563 else if (settings.SORInstead)
564 uRec[iCell] +=
565 relax *
566 (matrixAAInvBRow[ic2f + 1] * uRecOther +
567 vectorAInvBRow[ic2f] * (uOther - u[iCell].transpose()));
568 else
569 uRecNew[iCell] +=
570 relax *
571 (matrixAAInvBRow[ic2f + 1] * uRecOther +
572 vectorAInvBRow[ic2f] * (uOther - u[iCell].transpose()));
573 }
574 else
575 {
576 Eigen::Matrix<real, Eigen::Dynamic, nVarsFixed> BCC = GetBoundaryRHS(uRec, u, iCell, iFace, FBoundary);
577 // BCC *= 0;
578 if (recordInc)
579 {
580 if (uRecIsZero)
581 ;
582 else
583 uRecNew[iCell] -=
584 matrixAAInvBRow[0] * BCC;
585 }
586 else if (settings.SORInstead && !recordInc)
587 uRec[iCell] +=
588 relax * matrixAAInvBRow[0] * BCC;
589 else
590 uRecNew[iCell] +=
591 relax * matrixAAInvBRow[0] * BCC;
592 }
593 }
594 if ((!uRecNew[iCell].allFinite()) || (!uRec[iCell].allFinite()))
595 {
596 DNDS_assert(false);
597 }
598 }
599
600 if (!recordInc)
601 {
602 if (putIntoNew && settings.SORInstead)
603#if defined(DNDS_DIST_MT_USE_OMP)
604# pragma omp parallel for schedule(static)
605#endif
606 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
607 uRecNew[iCell].swap(uRec[iCell]);
608 if ((!putIntoNew) && (!settings.SORInstead))
609#if defined(DNDS_DIST_MT_USE_OMP)
610# pragma omp parallel for schedule(static)
611#endif
612 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
613 uRec[iCell].swap(uRecNew[iCell]);
614 }
615 }
616
617 template <int dim>
618 template <int nVarsFixed>
619 void VariationalReconstruction<dim>::DoReconstructionIterDiff(
620 tURec<nVarsFixed> &uRec,
621 tURec<nVarsFixed> &uRecDiff,
622 tURec<nVarsFixed> &uRecNew,
623 tUDof<nVarsFixed> &u,
624 const TFBoundaryDiff<nVarsFixed> &FBoundaryDiff)
625 {
626 using namespace Geom;
627 using namespace Geom::Elem;
628 using namespace Geom::Base;
629 int maxNDOF = GetNDof<dim>(settings.maxOrder);
630 static const auto Seq012 = Eigen::seq(Eigen::fix<0>, Eigen::fix<dim - 1>);
631
632 if (settings.maxOrder == 1 && settings.subs2ndOrder != 0)
633 {
634#if defined(DNDS_DIST_MT_USE_OMP)
635# pragma omp parallel for schedule(static)
636#endif
637 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
638 uRecNew[iCell].setZero();
639 return;
640 }
641#if defined(DNDS_DIST_MT_USE_OMP)
642# pragma omp parallel for schedule(runtime)
643#endif
644 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
645 {
646 uRecNew[iCell] = uRecDiff[iCell];
647
648 auto c2f = mesh->cell2face[iCell];
649 auto matrixAAInvBRow = coefficients_.matrixAAInvB[iCell];
650 for (int ic2f = 0; ic2f < c2f.size(); ic2f++)
651 {
652 index iFace = c2f[ic2f];
653 index iCellOther = CellFaceOther(iCell, iFace);
654 auto faceID = mesh->GetFaceZone(iFace);
655 int if2c = CellIsFaceBack(iCell, iFace) ? 0 : 1;
656 if (iCellOther != UnInitIndex)
657 {
658 Eigen::Matrix<real, Eigen::Dynamic, nVarsFixed> uRecOtherDiff = uRecDiff[iCellOther];
659 this->ApplyPeriodicTransform(if2c, faceID, uRecOtherDiff);
660 uRecNew[iCell] -= matrixAAInvBRow[ic2f + 1] * uRecOtherDiff; // mind the sign
661 }
662 else
663 {
664 Eigen::Matrix<real, Eigen::Dynamic, nVarsFixed> BCC = GetBoundaryRHSDiff(uRec, uRecDiff, u, iCell, iFace, FBoundaryDiff);
665 // BCC *= 0;
666 uRecNew[iCell] -= matrixAAInvBRow[0] * BCC; // mind the sign
667 }
668 }
669 }
670 }
671
672 template <int dim>
673 template <int nVarsFixed>
674 void VariationalReconstruction<dim>::
675 DoReconstructionIterSOR(
676 tURec<nVarsFixed> &uRec,
677 tURec<nVarsFixed> &uRecInc,
678 tURec<nVarsFixed> &uRecNew,
679 tUDof<nVarsFixed> &u,
680 const TFBoundaryDiff<nVarsFixed> &FBoundaryDiff,
681 bool reverse)
682 {
683 using namespace Geom;
684 using namespace Geom::Elem;
685 using namespace Geom::Base;
686 int maxNDOF = GetNDof<dim>(settings.maxOrder);
687 static const auto Seq012 = Eigen::seq(Eigen::fix<0>, Eigen::fix<dim - 1>);
688
689 for (index iScan = 0; iScan < mesh->NumCell(); iScan++)
690 {
691 index iCell = iScan;
692 if (reverse)
693 iCell = mesh->NumCell() - 1 - iCell;
694 real relax = GetCellAtr(iCell).relax;
695
696 uRecNew[iCell] = (1 - relax) * uRecNew[iCell] + uRecInc[iCell];
697
698 auto c2f = mesh->cell2face[iCell];
699 auto matrixAAInvBRow = coefficients_.matrixAAInvB[iCell];
700 for (int ic2f = 0; ic2f < c2f.size(); ic2f++)
701 {
702 index iFace = c2f[ic2f];
703 index iCellOther = CellFaceOther(iCell, iFace);
704 auto faceID = mesh->GetFaceZone(iFace);
705 int if2c = CellIsFaceBack(iCell, iFace) ? 0 : 1;
706 if (iCellOther != UnInitIndex)
707 {
708 Eigen::Matrix<real, Eigen::Dynamic, nVarsFixed> uRecOtherNew = uRecNew[iCellOther];
709 this->ApplyPeriodicTransform(if2c, faceID, uRecOtherNew);
710 uRecNew[iCell] += relax * matrixAAInvBRow[ic2f + 1] * uRecOtherNew; // mind the sign
711 }
712 else
713 {
714 Eigen::Matrix<real, Eigen::Dynamic, nVarsFixed> BCC = GetBoundaryRHSDiff(uRec, uRecNew, u, iCell, iFace, FBoundaryDiff);
715 // BCC *= 0;
716 uRecNew[iCell] += relax * matrixAAInvBRow[0] * BCC; // mind the sign
717 }
718 }
719 }
720 }
721 }
722}
DNDS_assert_info
#define DNDS_assert_info(expr, info)
Debug-only assertion with an extra std::string info message.
Definition Errors.hpp:113
DNDS_assert
#define DNDS_assert(expr)
Debug-only assertion (compiled out when DNDS_NDEBUG is defined). Prints the expression + file/line + ...
Definition Errors.hpp:108
DNDS_SWITCH_INTELLISENSE
#define DNDS_SWITCH_INTELLISENSE(real, intellisense)
Definition Macros.hpp:86
DNDS::ArrayDof
Primary solver state container: an ArrayEigenMatrix pair with MPI-collective vector-space operations.
Definition ArrayDOF.hpp:172
DNDS::CFV::VariationalReconstruction
The VR class that provides any information needed in high-order CFV.
Definition VariationalReconstruction.hpp:110
DNDS::CFV::VariationalReconstruction::DoReconstruction2ndGrad
void DoReconstruction2ndGrad(tUGrad< nVarsFixed, dim > &uRec, tUDof< nVarsFixed > &u, const TFBoundary< nVarsFixed > &FBoundary, int method)
fallback reconstruction method, explicit 2nd order FV reconstruction
DNDS::CFV::VariationalReconstruction::DoReconstructionIter
void DoReconstructionIter(tURec< nVarsFixed > &uRec, tURec< nVarsFixed > &uRecNew, tUDof< nVarsFixed > &u, const TFBoundary< nVarsFixed > &FBoundary, bool putIntoNew=false, bool recordInc=false, bool uRecIsZero=false)
iterative variational reconstruction method
Definition VariationalReconstruction_Reconstruction.hxx:487
DNDS::CFV::VariationalReconstruction::DoReconstruction2nd
void DoReconstruction2nd(tURec< nVarsFixed > &uRec, tUDof< nVarsFixed > &u, const TFBoundary< nVarsFixed > &FBoundary, int method, const std::vector< int > &mask)
fallback reconstruction method, explicit 2nd order FV reconstruction
DNDS::CFV::VariationalReconstruction::DoReconstructionIterSOR
void DoReconstructionIterSOR(tURec< nVarsFixed > &uRec, tURec< nVarsFixed > &uRecInc, tURec< nVarsFixed > &uRecNew, tUDof< nVarsFixed > &u, const TFBoundaryDiff< nVarsFixed > &FBoundaryDiff, bool reverse=false)
do a SOR iteration from uRecNew, with uRecInc as the RHSterm of Block-Jacobi preconditioned system
Definition VariationalReconstruction_Reconstruction.hxx:675
DNDS::CFV::VariationalReconstruction::DoReconstructionIterDiff
void DoReconstructionIterDiff(tURec< nVarsFixed > &uRec, tURec< nVarsFixed > &uRecDiff, tURec< nVarsFixed > &uRecNew, tUDof< nVarsFixed > &u, const TFBoundaryDiff< nVarsFixed > &FBoundaryDiff)
puts into uRecNew with Mat * uRecDiff; uses the Block-jacobi preconditioned reconstruction system as ...
Definition VariationalReconstruction_Reconstruction.hxx:619
DNDS::Geom::tSmallCoords
Eigen::Matrix< real, 3, Eigen::Dynamic > tSmallCoords
Definition Geometric.hpp:50
DNDS::Geom::FaceIDIsExternalBC
DNDS_DEVICE_CALLABLE bool FaceIDIsExternalBC(t_index id)
Definition BoundaryCondition.hpp:80
DNDS
the host side operators are provided as implemented
Definition BenchmarkFiniteVolume.cpp:7
DNDS::UnInitIndex
DNDS_CONSTANT const index UnInitIndex
Sentinel "not initialised" index value (= INT64_MIN).
Definition Defines.hpp:176
DNDS::verySmallReal
DNDS_CONSTANT const real verySmallReal
Catch-all lower bound ("effectively zero").
Definition Defines.hpp:194
DNDS::MatrixXR
Eigen::Matrix< real, Eigen::Dynamic, Eigen::Dynamic > MatrixXR
Column-major dynamic Eigen matrix of reals (default layout).
Definition Defines.hpp:205
DNDS::index
int64_t index
Global row / DOF index type (signed 64-bit; handles multi-billion-cell meshes).
Definition Defines.hpp:107
DNDS::RowVectorXR
Eigen::RowVector< real, Eigen::Dynamic > RowVectorXR
Dynamic row-vector of reals.
Definition Defines.hpp:209
DNDS::real
double real
Canonical floating-point scalar used throughout DNDSR (double precision).
Definition Defines.hpp:105