DNDSR 0.1.0.dev1+gcd065ad
Distributed Numeric Data Structure for CFV
Loading...
Searching...
No Matches
Limiters.hpp
Go to the documentation of this file.
1#pragma once
2#include "DNDS/Defines.hpp"
3
4namespace DNDS::CFV
5{
6 /**
7 * @brief Computes the weighted polynomial squared norm of each column of theta.
8 *
9 * The weights follow from the L2 inner product structure of the polynomial
10 * basis: pure-order terms get weight 1, mixed terms get reduced weights
11 * proportional to the combinatorial mixing of coordinate directions.
12 *
13 * For dim=2:
14 * nRows=2: theta(0)^2 + theta(1)^2
15 * nRows=3: theta(0)^2 + theta(1)^2 * 0.5 + theta(2)^2
16 * nRows=4: theta(0)^2 + theta(1)^2 * (1/3) + theta(2)^2 * (1/3) + theta(3)^2
17 *
18 * For dim=3:
19 * nRows=3: uniform weight 1
20 * nRows=6: first 3 weight 1, next 3 weight 0.5
21 * nRows=10: first 3 weight 1, next 6 weight 1/3, last 1 weight 1/6
22 *
23 * @tparam dim Spatial dimension (2 or 3)
24 * @param theta Array of shape (nRows, nCols), normalized polynomial coefficients
25 * @return Eigen::ArrayXd of shape (nCols) containing the weighted squared norms
26 */
27 template <int dim, typename ThetaArray>
28 inline Eigen::ArrayXd PolynomialSquaredNorm(const ThetaArray &theta)
29 {
30 Eigen::ArrayXd result;
31 if constexpr (dim == 2)
32 {
33 switch (theta.rows())
34 {
35 case 2:
36 result =
37 theta(0, EigenAll).square() +
38 theta(1, EigenAll).square();
39 break;
40 case 3:
41 result =
42 theta(0, EigenAll).square() +
43 theta(1, EigenAll).square() * 0.5 +
44 theta(2, EigenAll).square();
45 break;
46 case 4:
47 result =
48 theta(0, EigenAll).square() +
49 theta(1, EigenAll).square() * (1. / 3.) +
50 theta(2, EigenAll).square() * (1. / 3.) +
51 theta(3, EigenAll).square();
52 break;
53 default:
54 DNDS_assert(false);
55 break;
56 }
57 }
58 else
59 {
60 switch (theta.rows())
61 {
62 case 3:
63 result =
64 theta(0, EigenAll).square() +
65 theta(1, EigenAll).square() +
66 theta(2, EigenAll).square();
67 break;
68 case 6:
69 result =
70 theta(0, EigenAll).square() +
71 theta(1, EigenAll).square() +
72 theta(2, EigenAll).square() +
73 theta(3, EigenAll).square() * 0.5 +
74 theta(4, EigenAll).square() * 0.5 +
75 theta(5, EigenAll).square() * 0.5;
76 break;
77 case 10:
78 result =
79 theta(0, EigenAll).square() +
80 theta(1, EigenAll).square() +
81 theta(2, EigenAll).square() +
82 theta(3, EigenAll).square() * (1. / 3.) +
83 theta(4, EigenAll).square() * (1. / 3.) +
84 theta(5, EigenAll).square() * (1. / 3.) +
85 theta(6, EigenAll).square() * (1. / 3.) +
86 theta(7, EigenAll).square() * (1. / 3.) +
87 theta(8, EigenAll).square() * (1. / 3.) +
88 theta(9, EigenAll).square() * (1. / 6.);
89 break;
90 default:
91 DNDS_assert(false);
92 break;
93 }
94 }
95 return result;
96 }
97
98 /**
99 * @brief Computes a weighted polynomial dot product between theta1 and theta2.
100 *
101 * Uses the same weight structure as PolynomialSquaredNorm<dim>.
102 * Only used for dim=2 in FMEMM_Multiway_Polynomial2D.
103 */
104 template <int dim, typename ThetaArray1, typename ThetaArray2>
105 inline Eigen::ArrayXd PolynomialDotProduct(const ThetaArray1 &theta1, const ThetaArray2 &theta2)
106 {
107 Eigen::ArrayXd result;
108 static_assert(dim == 2, "PolynomialDotProduct only implemented for dim=2");
109 switch (theta1.rows())
110 {
111 case 2:
112 result =
113 theta1(0, EigenAll) * theta2(0, EigenAll) +
114 theta1(1, EigenAll) * theta2(1, EigenAll);
115 break;
116 case 3:
117 result =
118 theta1(0, EigenAll) * theta2(0, EigenAll) +
119 theta1(1, EigenAll) * theta2(1, EigenAll) * 0.5 +
120 theta1(2, EigenAll) * theta2(2, EigenAll);
121 break;
122 case 4:
123 result =
124 theta1(0, EigenAll) * theta2(0, EigenAll) +
125 theta1(1, EigenAll) * theta2(1, EigenAll) * (1. / 3.) +
126 theta1(2, EigenAll) * theta2(2, EigenAll) * (1. / 3.) +
127 theta1(3, EigenAll) * theta2(3, EigenAll);
128 break;
129 default:
130 DNDS_assert(false);
131 break;
132 }
133 return result;
134 }
135 /**
136 * @brief input vector<Eigen::Array-like>
137 */
138 template <typename TinOthers, typename Tout>
139 static inline void FWBAP_L2_Multiway_Polynomial2D(const TinOthers &uOthers, int Nother, Tout &uOut, real n1 = 1)
140 {
141 using namespace DNDS;
142 static const int p = 4;
143 static const real verySmallReal_pDiP = std::pow(verySmallReal, 1.0 / p);
144
145 Eigen::ArrayXXd uUp; //* copy!
146 uUp.resizeLike(uOthers[0]);
147 uUp.setZero();
148 Eigen::ArrayXd uDown;
149 uDown.resize(uOthers[0].cols());
150 uDown.setZero();
151 Eigen::ArrayXXd uMax = uUp + verySmallReal;
152 for (int iOther = 0; iOther < Nother; iOther++)
153 uMax = uMax.max(uOthers[iOther].abs());
154 uMax.rowwise() = uMax.colwise().maxCoeff();
155 uOut = uMax;
156
157 for (int iOther = 0; iOther < Nother; iOther++)
158 {
159 Eigen::ArrayXd thetaNorm;
160 Eigen::ArrayXXd theta = uOthers[iOther] / uMax;
161 thetaNorm = PolynomialSquaredNorm<2>(theta);
162 thetaNorm += verySmallReal_pDiP;
163 thetaNorm = thetaNorm.pow(-p / 2);
164
165 real exn = iOther ? 1.0 : n1;
166 uDown += thetaNorm * exn;
167 uUp += theta.rowwise() * thetaNorm.transpose() * exn;
168 }
169
170 uOut *= uUp.rowwise() / (uDown.transpose() + verySmallReal);
171
172 // * Do cut off
173 // for (int iOther = 0; iOther < Nother; iOther++)
174 // {
175 // Eigen::ArrayXd uDotuOut;
176 // switch (uOut.rows())
177 // {
178 // case 2:
179 // uDotuOut =
180 // uOthers[iOther](0, EigenAll) * uOut(0, EigenAll) +
181 // uOthers[iOther](1, EigenAll) * uOut(1, EigenAll);
182 // break;
183 // case 3:
184 // uDotuOut =
185 // uOthers[iOther](0, EigenAll) * uOut(0, EigenAll) +
186 // uOthers[iOther](1, EigenAll) * uOut(1, EigenAll) * 0.5 +
187 // uOthers[iOther](2, EigenAll) * uOut(2, EigenAll);
188 // break;
189 // case 4:
190 // uDotuOut =
191 // uOthers[iOther](0, EigenAll) * uOut(0, EigenAll) +
192 // uOthers[iOther](1, EigenAll) * uOut(1, EigenAll) * (1. / 3.) +
193 // uOthers[iOther](2, EigenAll) * uOut(2, EigenAll) * (1. / 3.) +
194 // uOthers[iOther](3, EigenAll) * uOut(3, EigenAll);
195 // break;
196
197 // default:
198 // DNDS_assert(false);
199 // break;
200 // }
201
202 // uOut.rowwise() *= 0.5 * (uDotuOut.sign().transpose() + 1);
203 // }
204 // // * Do cut off
205
206 if (uOut.hasNaN())
207 {
208 std::cout << "Limiter FWBAP_L2_Multiway Failed" << std::endl;
209 std::cout << uMax.transpose() << std::endl;
210 std::cout << uUp.transpose() << std::endl;
211 std::cout << uDown.transpose() << std::endl;
212 abort();
213 }
214 }
215
216 /**
217 * @brief input vector<Eigen::Array-like>
218 */
219 template <typename Tcenter, typename TinOthers, typename Tout>
220 static inline void FMEMM_Multiway_Polynomial2D(const Tcenter &u, const TinOthers &uOthers, int Nother, Tout &uOut, real n1 = 1)
221 {
222 using namespace DNDS;
223 static const int p = 4;
224
225 Eigen::ArrayXXd umax = u.abs();
226 umax.rowwise() = umax.colwise().maxCoeff() + verySmallReal;
227
228 Eigen::ArrayXXd theta0 = u / umax;
229 Eigen::ArrayXXd thetaMinNorm = theta0;
230
231 for (int iOther = 0; iOther < Nother; iOther++)
232 {
233 Eigen::ArrayXd thetaMinNormNorm;
234 Eigen::ArrayXd thetaNorm;
235 Eigen::ArrayXXd theta = uOthers[iOther] / umax;
236 Eigen::ArrayXd theta0DotTheta;
237 thetaNorm = PolynomialSquaredNorm<2>(theta);
238 thetaMinNormNorm = PolynomialSquaredNorm<2>(thetaMinNorm);
239 theta0DotTheta = PolynomialDotProduct<2>(theta, theta0);
240 Eigen::ArrayXd selection = (thetaNorm - thetaMinNormNorm).sign() * 0.5 + 0.5; //! need eliminate one side?
241 thetaMinNorm = theta.rowwise() * (1 - selection).transpose() +
242 thetaMinNorm.rowwise() * selection.transpose();
243 // //! cutting
244 // theta0 = theta0.rowwise() * (theta0DotTheta.sign() + 1).transpose() * 0.5;
245 }
246 Eigen::ArrayXd thetaNorm;
247 Eigen::ArrayXd thetaMinNormNorm;
248 thetaNorm = PolynomialSquaredNorm<2>(theta0);
249 thetaMinNormNorm = PolynomialSquaredNorm<2>(thetaMinNorm);
250 Eigen::ArrayXd replaceLoc = ((thetaNorm / (thetaMinNormNorm + verySmallReal)).sqrt() - 1).max(verySmallReal);
251 // Eigen::ArrayXd replaceFactor = 2 - (-replaceLoc).exp();
252 // Eigen::ArrayXd replaceFactor = 2 - (replaceLoc * p + 1).pow(-1. / p);
253 Eigen::ArrayXd replaceFactor = replaceLoc * 0 + 1;
254 // Eigen::ArrayXd replaceFactor = 1+(1 - (replaceLoc * p + 1).pow(-1. / p)) / (replaceLoc/10+1 );
255
256 replaceFactor = (replaceFactor - 1) / replaceLoc;
257
258 // !safety?
259 Eigen::ArrayXd ifReplace = (thetaNorm - thetaMinNormNorm).sign() * 0.5 + 0.5;
260 replaceFactor = ifReplace * replaceFactor + (1 - ifReplace);
261
262 uOut = u.rowwise() * replaceFactor.transpose() + (thetaMinNorm * umax).rowwise() * (1 - replaceFactor).transpose();
263
264 if (uOut.hasNaN())
265 {
266 std::cout << "Limiter FMEMM_L2_Multiway Failed" << std::endl;
267 std::cout << umax.transpose() << std::endl;
268 std::cout << uOut.transpose() << std::endl;
269 std::cout << replaceFactor << std::endl;
270 std::cout << replaceLoc << std::endl;
271 abort();
272 }
273 }
274
275 /**
276 * @brief input vector<Eigen::Array-like>
277 */
278 template <typename TinOthers, typename Tout>
279 static inline void FWBAP_L2_Multiway_PolynomialOrth(const TinOthers &uOthers, int Nother, Tout &uOut, real n1 = 1)
280 {
281 using namespace DNDS;
282 static const int p = 4;
283 static const real verySmallReal_pDiP = std::pow(verySmallReal, 1.0 / p);
284
285 Eigen::ArrayXXd uUp; //* copy!
286 uUp.resizeLike(uOthers[0]);
287 uUp.setZero();
288 Eigen::ArrayXd uDown;
289 uDown.resize(uOthers[0].cols());
290 uDown.setZero();
291 Eigen::ArrayXXd uMax = uUp + verySmallReal;
292 for (int iOther = 0; iOther < Nother; iOther++)
293 uMax = uMax.max(uOthers[iOther].abs());
294 uMax.rowwise() = uMax.colwise().maxCoeff();
295 uOut = uMax;
296
297 for (int iOther = 0; iOther < Nother; iOther++)
298 {
299 Eigen::ArrayXd thetaNorm;
300 Eigen::ArrayXXd theta = uOthers[iOther] / uMax;
301 thetaNorm = (theta * theta).colwise().sum();
302 thetaNorm += verySmallReal_pDiP;
303 thetaNorm = thetaNorm.pow(-p / 2);
304
305 real exn = iOther ? 1.0 : n1;
306 uDown += thetaNorm * exn;
307 uUp += theta.rowwise() * thetaNorm.transpose() * exn;
308 }
309 // std::cout << uUp << std::endl;
310 // std::cout << uDown << std::endl;
311 uOut *= uUp.rowwise() / (uDown.transpose() + verySmallReal);
312 if (uOut.hasNaN())
313 {
314 std::cout << "Limiter FWBAP_L2_Multiway Failed" << std::endl;
315 std::cout << uMax.transpose() << std::endl;
316 std::cout << uUp.transpose() << std::endl;
317 std::cout << uDown.transpose() << std::endl;
318 abort();
319 }
320 }
321
322 namespace _Limiters_Internal
323 {
324 inline real W12n1(real u1, real u2)
325 {
326 real n = 1.0;
327 real p = 4.0;
328
329 real theta1 = std::pow(u1 / (u2 + signP(u2) * 1e-12), p - 1.0);
330 real theta2 = std::pow(u1 / (u2 + signP(u2) * 1e-12), p);
331
332 return u1 * (n + theta1) / (n + theta2);
333 }
334
335 inline real W12center(real *u, const int J, real n)
336 {
337 std::vector<real> thetaBuf(J);
338 real *theta = thetaBuf.data();
339 for (int ii = 0; ii < J; ++ii)
340 {
341 theta[ii] = (u[0] + signP(u[0]) * 1e-12) / (u[ii] + signP(u[ii]) * 1e-12);
342 }
343
344 real p = 4.0;
345 real sumLocal1 = n;
346 real sumLocal2 = n;
347 for (int ii = 0; ii < J; ++ii)
348 {
349 sumLocal1 += std::pow(theta[ii], (p - 1.0));
350 sumLocal2 += std::pow(theta[ii], p);
351 }
352
353 return u[0] * sumLocal1 / (sumLocal2 + 1e-12);
354 }
355 }
356
357 /**
358 * @brief input vector<Eigen::Array-like>
359 */
360 template <typename TinOthers, typename Tout>
361 inline void FWBAP_L2_Multiway(const TinOthers &uOthers, int Nother, Tout &uOut, real n1 = 1.0)
362 {
363 // PerformanceTimer::Instance().StartTimer(PerformanceTimer::LimiterA);
364 static const int p = 4;
365 static const real verySmallReal_pDiP = std::pow(verySmallReal, 1.0 / p);
366
367 Tout uUp; //* copy!
368 uUp.resizeLike(uOthers[0]);
369 uUp.setZero();
370 Tout uDown = uUp; //* copy!
371 Tout uCent = uOthers[0];
372 // Tout uMax = uDown + verySmallReal;
373
374 // for (int iOther = 0; iOther < Nother; iOther++)
375 // uMax = uMax.max(uOthers[iOther].abs());
376 // uOut = uMax;
377
378 for (int iOther = 1; iOther < Nother; iOther++)
379 {
380 Tout thetaInverse = uOthers[0] / (uOthers[iOther].abs() + verySmallReal_pDiP) * uOthers[iOther].sign();
381 uDown += thetaInverse.square() * thetaInverse.square();
382 uUp += thetaInverse.cube();
383 uCent *= uOthers[iOther].sign().abs();
384 }
385 uOut = uCent * (uUp + n1) / (uDown + n1); // currently fast version
386
387 // if (uOut.hasNaN())
388 // {
389 // std::cout << "Limiter FWBAP_L2_Multiway Failed" << std::endl;
390 // std::cout << uMax.transpose() << std::endl;
391 // std::cout << uUp.transpose() << std::endl;
392 // std::cout << uDown.transpose() << std::endl;
393 // abort();
394 // }
395
396 /************************/ //! safe version
397 // uOut.resizeLike(uOthers[0]);
398 // static const int maxNeighbour = 7;
399 // DNDS_assert(uOthers.size() <= maxNeighbour);
400 // real theta[maxNeighbour];
401
402 // for (int idof = 0; idof < uOthers[0].cols(); idof++)
403 // for (int irec = 0; irec < uOthers[0].rows(); irec++)
404 // {
405 // real u0 = uOthers[0](irec, idof);
406 // for (int ii = 0; ii < uOthers.size(); ++ii)
407 // {
408 // real uother = uOthers[ii](irec, idof);
409 // theta[ii] = (u0 + signM(u0) * 1e-12) /
410 // (uother + signM(uother) * 1e-12);
411 // }
412
413 // static const real p = 4.0;
414 // real sumLocal1 = n1;
415 // real sumLocal2 = n1;
416 // for (int ii = 1; ii < uOthers.size(); ++ii)
417 // {
418 // sumLocal1 += std::pow(theta[ii], (p - 1.0));
419 // sumLocal2 += std::pow(theta[ii], p);
420 // }
421
422 // uOut(irec, idof) = u0 * sumLocal1 / (sumLocal2 + 1e-12);
423 // }
424
425 // PerformanceTimer::Instance().StopTimer(PerformanceTimer::LimiterA);
426 }
427
428 /**
429 * @brief input eigen arrays
430 */
431 template <typename Tin1, typename Tin2, typename Tout>
432 inline void FWBAP_L2_Biway(const Tin1 &u1, const Tin2 &u2, Tout &uOut, real n)
433 {
434 // PerformanceTimer::Instance().StartTimer(PerformanceTimer::LimiterA);
435 static const int p = 4;
436 static const real verySmallReal_pDiP = std::pow(verySmallReal, 1.0 / p);
437 // // static const real n = 10.0;
438 // static const real verySmallReal_pDiP = std::pow(verySmallReal, 1.0 / p);
439 // auto uMax = u1.abs().max(u2.abs()) + verySmallReal_pDiP;
440 // auto u1p = (u1 / uMax).pow(p);
441 // auto u2p = (u2 / uMax).pow(p);
442 // // std::cout << u1 << std::endl;
443
444 // uOut = (u1p * u2 + n * u2p * u1) / ((u1p + n * u2p) + verySmallReal);
445 // // uOut *= (u1.sign() + u2.sign()).abs() * 0.5; //! cutting below zero!!!
446 // // std::cout << u2 << std::endl;
447
448 ///////////
449 Tout frac = (u1 / (u2.abs() + verySmallReal_pDiP) * u2.sign());
450
451 // for (index i = 0; i < frac.size(); i++)
452 // {
453 // auto &v = frac(i);
454 // if (v == 0.)
455 // v = 1. / verySmallReal_pDiP;
456 // }
457
458 // auto theta1 = frac.pow(p - 1);
459 // auto theta2 = frac.pow(p);
460
461 auto theta1 = frac.cube();
462 auto theta2 = frac.square() * frac.square();
463
464 uOut = u1 * (n + theta1) / (n + theta2); // currently fast version
465 uOut *= u2.sign().abs();
466 ///////////
467 /************************/ //! safe version
468 // uOut.resizeLike(u1);
469 // for (int idof = 0; idof < u1.cols(); idof++)
470 // for (int irec = 0; irec < u1.rows(); irec++)
471 // {
472 // real u1c = u1(irec, idof);
473 // real u2c = u2(irec, idof);
474 // real frac = (u1c) / (u2c + signM(u2c) * 1e-12);
475 // real theta1 = std::pow(frac, p - 1);
476 // real theta2 = std::pow(frac, p);
477 // uOut(irec, idof) = u1c * (n + theta1) / (n + theta2);
478 // }
479 // PerformanceTimer::Instance().StopTimer(PerformanceTimer::LimiterA);
480 }
481
482 template <typename Tin1, typename Tin2, typename Tout>
483 inline void FWBAP_L2_Cut_Biway(const Tin1 &u1, const Tin2 &u2, Tout &uOut, real n)
484 {
485 // PerformanceTimer::Instance().StartTimer(PerformanceTimer::LimiterA);
486 static const int p = 4;
487 static const real verySmallReal_pDiP = std::pow(verySmallReal, 1.0 / p);
488
489 Tout frac = (u1 / (u2.abs() + verySmallReal_pDiP) * u2.sign());
490
491 auto theta1 = frac.cube();
492 auto theta2 = frac.square() * frac.square();
493
494 uOut = u1 * (n + theta1) / (n + theta2); // currently fast version
495 uOut *= u2.sign().abs() * (0.5 * (u1.sign() + u2.sign()).abs());
496 ///////////
497 }
498
499 /**
500 * @brief input eigen arrays
501 */
502 template <typename Tin1, typename Tin2, typename Tout>
503 inline void FMINMOD_Biway(const Tin1 &u1, const Tin2 &u2, Tout &uOut, real n)
504 {
505 uOut = u1.abs().min(u2.abs()) * (u1.sign() + u2.sign()) * 0.5;
506 }
507
508 /**
509 * @brief input eigen arrays
510 */
511 template <typename Tin1, typename Tin2, typename Tout>
512 inline void FVanLeer_Biway(const Tin1 &u1, const Tin2 &u2, Tout &uOut, real n)
513 {
514 uOut = (u1.sign() + u2.sign()) * (u1 * u2).abs() / ((u1 + u2).abs() + verySmallReal);
515 }
516
517 template <int dim, int nVarsFixed, typename Tin1, typename Tin2, typename Tout>
518 inline void FWBAP_L2_Biway_PolynomialNorm(const Tin1 &u1, const Tin2 &u2, Tout &uOut, real n)
519 {
520 static const int p = 4;
521 // static const real n = 10.0;
522 static const real verySmallReal_pDiP = std::pow(verySmallReal, 1.0 / p);
523 static const int nDofMax = dim == 2 ? 6 : 10; //! hard code max 4thorder
524 Eigen::Array<real, Eigen::Dynamic, nVarsFixed, Eigen::AutoAlign, nDofMax> uMax = u1.abs().max(u2.abs()) + verySmallReal_pDiP;
525 uMax.rowwise() = uMax.colwise().maxCoeff();
526 Eigen::ArrayXd u1p, u2p;
527 Eigen::Array<real, Eigen::Dynamic, nVarsFixed, Eigen::AutoAlign, nDofMax> theta1 = u1 / uMax;
528 Eigen::Array<real, Eigen::Dynamic, nVarsFixed, Eigen::AutoAlign, nDofMax> theta2 = u2 / uMax;
529 u1p = PolynomialSquaredNorm<dim>(theta1);
530 u2p = PolynomialSquaredNorm<dim>(theta2);
531 // u1p = u1p.pow(p / 2);
532 // u2p = u2p.pow(p / 2);
533
534 u1p = u1p.square();
535 u2p = u2p.square();
536
537 uOut = (u2.rowwise() * u1p.transpose() + n * (u1.rowwise() * u2p.transpose())).rowwise() / ((u1p + n * u2p) + verySmallReal).transpose();
538 }
539
540 template <int dim, int nVarsFixed, typename Tin1, typename Tin2, typename Tout>
541 inline void FMEMM_Biway_PolynomialNorm(const Tin1 &u1, const Tin2 &u2, Tout &uOut, real n)
542 {
543 static const int p = 4;
544 // static const real n = 10.0;
545 static const real verySmallReal_pDiP = std::pow(verySmallReal, 1.0 / p);
546 static const int nDofMax = dim == 2 ? 6 : 10; //! hard code max 4thorder
547 Eigen::Array<real, Eigen::Dynamic, nVarsFixed, Eigen::AutoAlign, nDofMax> uMax = u1.abs().max(u2.abs()) + verySmallReal_pDiP;
548 uMax.rowwise() = uMax.colwise().maxCoeff();
549 Eigen::ArrayXd u1p, u2p;
550 Eigen::Array<real, Eigen::Dynamic, nVarsFixed, Eigen::AutoAlign, nDofMax> theta1 = u1 / uMax;
551 Eigen::Array<real, Eigen::Dynamic, nVarsFixed, Eigen::AutoAlign, nDofMax> theta2 = u2 / uMax;
552 u1p = PolynomialSquaredNorm<dim>(theta1);
553 u2p = PolynomialSquaredNorm<dim>(theta2);
554 u1p = u1p.sqrt();
555 u2p = u2p.sqrt();
556
557 Eigen::ArrayXd replaceLoc = (u1p / (u2p + verySmallReal) - 1).max(verySmallReal);
558 // Eigen::ArrayXd replaceFactor = 2 - (-replaceLoc).exp();
559 // Eigen::ArrayXd replaceFactor = 2 - (replaceLoc * p + 1).pow(-1. / p);
560 Eigen::ArrayXd replaceFactor = replaceLoc * 0 + 1;
561 // Eigen::ArrayXd replaceFactor = 1+(1 - (replaceLoc * p + 1).pow(-1. / p)) / (replaceLoc/10+1 );
562
563 replaceFactor = (replaceFactor - 1) / replaceLoc;
564
565 // !safety?
566 Eigen::ArrayXd ifReplace = (u1p - u2p).sign() * 0.5 + 0.5;
567 // replaceFactor = ifReplace * replaceFactor + (1 - ifReplace);
568 // ! fast version:
569 replaceFactor = (1 - ifReplace);
570
571 uOut = u1.rowwise() * replaceFactor.transpose() + u2.rowwise() * (1 - replaceFactor).transpose();
572 // //! cutting
573 // uOut = uOut.rowwise() * (u1u2.sign() + 1).transpose() * 0.5;
574 }
575
576 template <typename Tin1, typename Tin2, typename Tout>
577 inline void FWBAP_L2_Biway_PolynomialOrth(const Tin1 &u1, const Tin2 &u2, Tout &uOut, real n)
578 {
579 static const int p = 4;
580 // static const real n = 10.0;
581 static const real verySmallReal_pDiP = std::pow(verySmallReal, 1.0 / p);
582 Eigen::ArrayXXd uMax = u1.abs().max(u2.abs()) + verySmallReal_pDiP;
583 uMax.rowwise() = uMax.colwise().maxCoeff();
584 Eigen::ArrayXd u1p, u2p;
585 Eigen::ArrayXXd theta1 = u1 / uMax;
586 Eigen::ArrayXXd theta2 = u2 / uMax;
587 u1p = (theta1 * theta1).colwise().sum();
588 u2p = (theta2 * theta2).colwise().sum();
589 u1p = u1p.pow(p / 2);
590 u2p = u2p.pow(p / 2);
591
592 uOut = (u2.rowwise() * u1p.transpose() + n * (u1.rowwise() * u2p.transpose())).rowwise() / ((u1p + n * u2p) + verySmallReal).transpose();
593 }
594
595}
Defines.hpp
Core type aliases, constants, and metaprogramming utilities for the DNDS framework.
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::CFV::_Limiters_Internal::W12n1
real W12n1(real u1, real u2)
Definition Limiters.hpp:324
DNDS::CFV::_Limiters_Internal::W12center
real W12center(real *u, const int J, real n)
Definition Limiters.hpp:335
DNDS::CFV::FWBAP_L2_Cut_Biway
void FWBAP_L2_Cut_Biway(const Tin1 &u1, const Tin2 &u2, Tout &uOut, real n)
Definition Limiters.hpp:483
DNDS::CFV::FWBAP_L2_Multiway
void FWBAP_L2_Multiway(const TinOthers &uOthers, int Nother, Tout &uOut, real n1=1.0)
input vector<Eigen::Array-like>
Definition Limiters.hpp:361
DNDS::CFV::FMEMM_Biway_PolynomialNorm
void FMEMM_Biway_PolynomialNorm(const Tin1 &u1, const Tin2 &u2, Tout &uOut, real n)
Definition Limiters.hpp:541
DNDS::CFV::FMINMOD_Biway
void FMINMOD_Biway(const Tin1 &u1, const Tin2 &u2, Tout &uOut, real n)
input eigen arrays
Definition Limiters.hpp:503
DNDS::CFV::FWBAP_L2_Biway_PolynomialOrth
void FWBAP_L2_Biway_PolynomialOrth(const Tin1 &u1, const Tin2 &u2, Tout &uOut, real n)
Definition Limiters.hpp:577
DNDS::CFV::FVanLeer_Biway
void FVanLeer_Biway(const Tin1 &u1, const Tin2 &u2, Tout &uOut, real n)
input eigen arrays
Definition Limiters.hpp:512
DNDS::CFV::FWBAP_L2_Biway_PolynomialNorm
void FWBAP_L2_Biway_PolynomialNorm(const Tin1 &u1, const Tin2 &u2, Tout &uOut, real n)
Definition Limiters.hpp:518
DNDS::CFV::PolynomialSquaredNorm
Eigen::ArrayXd PolynomialSquaredNorm(const ThetaArray &theta)
Computes the weighted polynomial squared norm of each column of theta.
Definition Limiters.hpp:28
DNDS::CFV::FWBAP_L2_Biway
void FWBAP_L2_Biway(const Tin1 &u1, const Tin2 &u2, Tout &uOut, real n)
input eigen arrays
Definition Limiters.hpp:432
DNDS::CFV::PolynomialDotProduct
Eigen::ArrayXd PolynomialDotProduct(const ThetaArray1 &theta1, const ThetaArray2 &theta2)
Computes a weighted polynomial dot product between theta1 and theta2.
Definition Limiters.hpp:105
DNDS
the host side operators are provided as implemented
Definition BenchmarkFiniteVolume.cpp:7
DNDS::verySmallReal
DNDS_CONSTANT const real verySmallReal
Catch-all lower bound ("effectively zero").
Definition Defines.hpp:194
DNDS::sign
constexpr real sign(real a)
Signum function: +1, 0, or -1.
Definition Defines.hpp:532
DNDS::real
double real
Canonical floating-point scalar used throughout DNDSR (double precision).
Definition Defines.hpp:105
DNDS::signP
constexpr real signP(real a)
"Signum, biased toward +1": treats 0 as positive.
Definition Defines.hpp:544
result
auto result
Definition test_Limiters.cpp:110
n
Eigen::Vector3d n(1.0, 0.0, 0.0)