Alsvinn  0.5.3
The fast FVM simulator with UQ support
WENOF2.hpp
Go to the documentation of this file.
1 /* Copyright (c) 2018 ETH Zurich, Kjetil Olsen Lye
2  * This program is free software: you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation, either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program. If not, see <http://www.gnu.org/licenses/>.
14  */
15 
16 #pragma once
17 #include "alsfvm/types.hpp"
20 #include <algorithm>
21 
22 namespace alsfvm {
23 namespace reconstruction {
31 template<class Equation>
32 class WENOF2 {
33 public:
34  __device__ __host__ static void reconstruct(Equation eq,
35  typename Equation::ConstViews& in,
36  size_t x, size_t y, size_t z,
37  typename Equation::Views& leftView,
38  typename Equation::Views& rightView,
39  bool xDir, bool yDir, bool zDir) {
40  const real BOUND = 1.9;
41  const size_t indexOut = leftView.index(x, y, z);
42  const size_t indexRight = leftView.index(x + xDir, y + yDir, z + zDir);
43  const size_t indexLeft = leftView.index(x - xDir, y - yDir, z - zDir);
44  const real i0 = eq.getWeight(in, indexOut);
45  const real b0 = square(eq.getWeight(in, indexRight) - i0);
46  const real b1 = square(i0 - eq.getWeight(in, indexLeft));
47  const real a0Left = 1 / (3 * (ALSFVM_WENO_EPSILON + b0) *
48  (ALSFVM_WENO_EPSILON + b0));
49  const real a1Left = 2 / (3 * (ALSFVM_WENO_EPSILON + b1) *
50  (ALSFVM_WENO_EPSILON + b1));
51  const real w0Left = a0Left / (a0Left + a1Left);
52  const real w1Left = a1Left / (a0Left + a1Left);
53 
54  const real a0Right = 2 / (3 * (ALSFVM_WENO_EPSILON + b0) *
55  (ALSFVM_WENO_EPSILON + b0));
56  const real a1Right = 1 / (3 * (ALSFVM_WENO_EPSILON + b1) *
57  (ALSFVM_WENO_EPSILON + b1));
58  const real w0Right = a0Right / (a0Right + a1Right);
59  const real w1Right = a1Right / (a0Right + a1Right);
60 
61  typename Equation::PrimitiveVariables inLeft =
62  eq.computePrimitiveVariables(eq.fetchConservedVariables(in, indexLeft));
63 
64  typename Equation::PrimitiveVariables inMiddle =
65  eq.computePrimitiveVariables(eq.fetchConservedVariables(in, indexOut));
66 
67  typename Equation::PrimitiveVariables inRight =
68  eq.computePrimitiveVariables(eq.fetchConservedVariables(in, indexRight));
69 
70  typename Equation::PrimitiveVariables dPLeft = w1Left * (inMiddle - inLeft) +
71  w0Left * (inRight - inMiddle);
72 
73 
74  dPLeft.rho = fmax(-BOUND * inMiddle.rho, fmin(BOUND * inMiddle.rho,
75  dPLeft.rho));
76  dPLeft.p = fmax(-BOUND * inMiddle.p, fmin(BOUND * inMiddle.p, dPLeft.p));
77 
78  real LLLeft = 0.125 * inMiddle.rho * (dPLeft.u.dot(dPLeft.u)) -
79  0.5 * fmin(real(0.0), dPLeft.rho * (inMiddle.u.dot(dPLeft.u))) +
80  0.5 * dPLeft.rho * dPLeft.rho * dPLeft.u.dot(dPLeft.u) / inMiddle.rho;
81 
82 
83  real R = inMiddle.p / (eq.getGamma() - 1);
84  real aijkLeft = 0.5 * std::sqrt(R / fmax(R, LLLeft));
85 
86 
87 
88 
89 
90  typename Equation::ConservedVariables left = eq.computeConserved(
91  inMiddle - aijkLeft * dPLeft);
92 
93 
94  typename Equation::PrimitiveVariables dPRight = w1Right * (inMiddle - inLeft) +
95  w0Right * (inRight - inMiddle);
96 
97 
98  dPRight.rho = fmax(-BOUND * inMiddle.rho, fmin(BOUND * inMiddle.rho,
99  dPRight.rho));
100  dPRight.p = fmax(-BOUND * inMiddle.p, fmin(BOUND * inMiddle.p, dPRight.p));
101 
102  real LLRight = 0.125 * inMiddle.rho * (dPRight.u.dot(dPRight.u)) -
103  0.5 * fmin(real(0.0), dPRight.rho * (inMiddle.u.dot(dPRight.u))) +
104  0.5 * dPRight.rho * dPRight.rho * dPRight.u.dot(dPRight.u) / inMiddle.rho;
105 
106  real aijkRight = 0.5 * std::sqrt(R / fmax(R, LLRight));
107 
108  typename Equation::ConservedVariables right = eq.computeConserved(
109  inMiddle + aijkRight * dPRight);
110  eq.setViewAt(leftView, indexOut, left);
111  eq.setViewAt(rightView, indexOut, right);
112 
113  }
114 
116  return 2;
117  }
118 };
119 } // namespace alsfvm
120 } // namespace reconstruction
alsfvm::shared_ptr< reconstruction::Reconstruction > & reconstruction
Definition: NumericalFluxFactory.cpp:101
#define ALSFVM_WENO_EPSILON
Definition: WENOCoefficients.hpp:21
#define __host__
Definition: types.hpp:46
__device__ static __host__ int getNumberOfGhostCells()
Definition: WENOF2.hpp:115
double real
Definition: types.hpp:65
__device__ static __host__ void reconstruct(Equation eq, typename Equation::ConstViews &in, size_t x, size_t y, size_t z, typename Equation::Views &leftView, typename Equation::Views &rightView, bool xDir, bool yDir, bool zDir)
Definition: WENOF2.hpp:34
Definition: WENOF2.hpp:32
#define __device__
Definition: types.hpp:45
Various utility functions to implement the tecno flux.
Definition: types.hpp:30
__host__ __device__ real square(const real &x)
Definition: types.hpp:159