PDE-Impl.hpp
Go to the documentation of this file.
1 // This code is based on Jet framework.
2 // Copyright (c) 2018 Doyub Kim
3 // CubbyFlow is voxel-based fluid simulation engine for computer games.
4 // Copyright (c) 2020 CubbyFlow Team
5 // Core Part: Chris Ohk, Junwoo Hwang, Jihong Sin, Seungwoo Yoo
6 // AI Part: Dongheon Cho, Minseo Kim
7 // We are making my contributions/submissions to this project solely in our
8 // personal capacity and are not conveying any rights to any intellectual
9 // property of any third parties.
10 
11 #ifndef CUBBYFLOW_PDE_IMPL_HPP
12 #define CUBBYFLOW_PDE_IMPL_HPP
13 
14 #include <Core/Math/MathUtils.hpp>
15 
16 namespace CubbyFlow
17 {
18 template <typename T>
19 std::array<T, 2> Upwind1(T* d0, T dx)
20 {
21  const T invdx = 1 / dx;
22  std::array<T, 2> dfx{};
23 
24  dfx[0] = invdx * (d0[1] - d0[0]);
25  dfx[1] = invdx * (d0[2] - d0[1]);
26 
27  return dfx;
28 }
29 
30 template <typename T>
31 T Upwind1(T* d0, T dx, bool isDirectionPositive)
32 {
33  const T invdx = 1 / dx;
34  return isDirectionPositive ? (invdx * (d0[1] - d0[0]))
35  : invdx * (d0[2] - d0[1]);
36 }
37 
38 template <typename T>
39 T CD2(T* d0, T dx)
40 {
41  const T hinvdx = 0.5f / dx;
42  return hinvdx * (d0[2] - d0[0]);
43 }
44 
45 template <typename T>
46 std::array<T, 2> ENO3(T* d0, T dx)
47 {
48  const T invdx = 1 / dx;
49  const T hinvdx = invdx / 2;
50  const T tinvdx = invdx / 3;
51  T d1[6], d2[5], d3[2];
52  T c, cstar;
53  int kstar;
54  std::array<T, 2> dfx{};
55 
56  d1[0] = invdx * (d0[1] - d0[0]);
57  d1[1] = invdx * (d0[2] - d0[1]);
58  d1[2] = invdx * (d0[3] - d0[2]);
59  d1[3] = invdx * (d0[4] - d0[3]);
60  d1[4] = invdx * (d0[5] - d0[4]);
61  d1[5] = invdx * (d0[6] - d0[5]);
62 
63  d2[0] = hinvdx * (d1[1] - d1[0]);
64  d2[1] = hinvdx * (d1[2] - d1[1]);
65  d2[2] = hinvdx * (d1[3] - d1[2]);
66  d2[3] = hinvdx * (d1[4] - d1[3]);
67  d2[4] = hinvdx * (d1[5] - d1[4]);
68 
69  for (int k = 0; k < 2; ++k)
70  {
71  if (std::fabs(d2[k + 1]) < std::fabs(d2[k + 2]))
72  {
73  c = d2[k + 1];
74  kstar = k - 1;
75  d3[0] = tinvdx * (d2[k + 1] - d2[k]);
76  d3[1] = tinvdx * (d2[k + 2] - d2[k + 1]);
77  }
78  else
79  {
80  c = d2[k + 2];
81  kstar = k;
82  d3[0] = tinvdx * (d2[k + 2] - d2[k + 1]);
83  d3[1] = tinvdx * (d2[k + 3] - d2[k + 2]);
84  }
85 
86  if (std::fabs(d3[0]) < std::fabs(d3[1]))
87  {
88  cstar = d3[0];
89  }
90  else
91  {
92  cstar = d3[1];
93  }
94 
95  T dq1 = d1[k + 2];
96  T dq2 = c * (2 * (1 - k) - 1) * dx;
97  T dq3 = cstar * (3 * Square(1 - kstar) - 6 * (1 - kstar) + 2) * dx * dx;
98 
99  dfx[k] = dq1 + dq2 + dq3;
100  }
101 
102  return dfx;
103 }
104 
105 template <typename T>
106 T ENO3(T* d0, T dx, bool isDirectionPositive)
107 {
108  const T invdx = 1 / dx;
109  const T hinvdx = invdx / 2;
110  const T tinvdx = invdx / 3;
111  T d1[6], d2[5], d3[2];
112  T c, cstar;
113  int kstar;
114 
115  d1[0] = invdx * (d0[1] - d0[0]);
116  d1[1] = invdx * (d0[2] - d0[1]);
117  d1[2] = invdx * (d0[3] - d0[2]);
118  d1[3] = invdx * (d0[4] - d0[3]);
119  d1[4] = invdx * (d0[5] - d0[4]);
120  d1[5] = invdx * (d0[6] - d0[5]);
121 
122  d2[0] = hinvdx * (d1[1] - d1[0]);
123  d2[1] = hinvdx * (d1[2] - d1[1]);
124  d2[2] = hinvdx * (d1[3] - d1[2]);
125  d2[3] = hinvdx * (d1[4] - d1[3]);
126  d2[4] = hinvdx * (d1[5] - d1[4]);
127 
128  int k = isDirectionPositive ? 0 : 1;
129 
130  if (std::fabs(d2[k + 1]) < std::fabs(d2[k + 2]))
131  {
132  c = d2[k + 1];
133  kstar = k - 1;
134  d3[0] = tinvdx * (d2[k + 1] - d2[k]);
135  d3[1] = tinvdx * (d2[k + 2] - d2[k + 1]);
136  }
137  else
138  {
139  c = d2[k + 2];
140  kstar = k;
141  d3[0] = tinvdx * (d2[k + 2] - d2[k + 1]);
142  d3[1] = tinvdx * (d2[k + 3] - d2[k + 2]);
143  }
144 
145  if (std::fabs(d3[0]) < std::fabs(d3[1]))
146  {
147  cstar = d3[0];
148  }
149  else
150  {
151  cstar = d3[1];
152  }
153 
154  T dq1 = d1[k + 2];
155  T dq2 = c * (2 * (1 - k) - 1) * dx;
156  T dq3 = cstar * (3 * Square(1 - kstar) - 6 * (1 - kstar) + 2) * dx * dx;
157 
158  return dq1 + dq2 + dq3;
159 }
160 
161 template <typename T>
162 std::array<T, 2> WENO5(T* v, T h, T eps)
163 {
164  static const T c13 = T(1.0 / 3.0), c14 = T(0.25), c16 = T(1.0 / 6.0);
165  static const T c56 = T(5.0 / 6.0), c76 = T(7.0 / 6.0), c116 = T(11.0 / 6.0);
166  static const T c1312 = T(13.0 / 12.0);
167 
168  const T hinv = T(1) / h;
169  std::array<T, 2> dfx{};
170  T vdev[5];
171 
172  for (int k = 0; k < 2; ++k)
173  {
174  if (k == 0)
175  {
176  for (int m = 0; m < 5; ++m)
177  {
178  vdev[m] = (v[m + 1] - v[m]) * hinv;
179  }
180  }
181  else
182  {
183  for (int m = 0; m < 5; ++m)
184  {
185  vdev[m] = (v[6 - m] - v[5 - m]) * hinv;
186  }
187  }
188 
189  const T phix1 = vdev[0] * c13 - vdev[1] * c76 + vdev[2] * c116;
190  const T phix2 = -vdev[1] * c16 + vdev[2] * c56 + vdev[3] * c13;
191  const T phix3 = vdev[2] * c13 + vdev[3] * c56 - vdev[4] * c16;
192 
193  T s1 = c1312 * Square(vdev[0] - 2 * vdev[1] + vdev[2]) +
194  c14 * Square(vdev[0] - 4 * vdev[1] + 3 * vdev[2]);
195  T s2 = c1312 * Square(vdev[1] - 2 * vdev[2] + vdev[3]) +
196  c14 * Square(vdev[1] - vdev[3]);
197  T s3 = c1312 * Square(vdev[2] - 2 * vdev[3] + vdev[4]) +
198  c14 * Square(3 * vdev[2] - 4 * vdev[3] + vdev[4]);
199 
200  T alpha1 = T(0.1 / Square(s1 + eps));
201  T alpha2 = T(0.6 / Square(s2 + eps));
202  T alpha3 = T(0.3 / Square(s3 + eps));
203 
204  T sum = alpha1 + alpha2 + alpha3;
205 
206  dfx[k] = (alpha1 * phix1 + alpha2 * phix2 + alpha3 * phix3) / sum;
207  }
208 
209  return dfx;
210 }
211 
212 template <typename T>
213 T WENO5(T* v, T h, bool isDirectionPositive, T eps)
214 {
215  static const T c13 = T(1.0 / 3.0), c14 = T(0.25), c16 = T(1.0 / 6.0);
216  static const T c56 = T(5.0 / 6.0), c76 = T(7.0 / 6.0), c116 = T(11.0 / 6.0);
217  static const T c1312 = T(13.0 / 12.0);
218 
219  const T hinv = T(1) / h;
220  T vdev[5];
221 
222  if (isDirectionPositive)
223  {
224  for (int m = 0; m < 5; ++m)
225  {
226  vdev[m] = (v[m + 1] - v[m]) * hinv;
227  }
228  }
229  else
230  {
231  for (int m = 0; m < 5; ++m)
232  {
233  vdev[m] = (v[6 - m] - v[5 - m]) * hinv;
234  }
235  }
236 
237  const T phix1 = vdev[0] * c13 - vdev[1] * c76 + vdev[2] * c116;
238  const T phix2 = -vdev[1] * c16 + vdev[2] * c56 + vdev[3] * c13;
239  const T phix3 = vdev[2] * c13 + vdev[3] * c56 - vdev[4] * c16;
240 
241  T s1 = c1312 * Square(vdev[0] - 2 * vdev[1] + vdev[2]) +
242  c14 * Square(vdev[0] - 4 * vdev[1] + 3 * vdev[2]);
243  T s2 = c1312 * Square(vdev[1] - 2 * vdev[2] + vdev[3]) +
244  c14 * Square(vdev[1] - vdev[3]);
245  T s3 = c1312 * Square(vdev[2] - 2 * vdev[3] + vdev[4]) +
246  c14 * Square(3 * vdev[2] - 4 * vdev[3] + vdev[4]);
247 
248  T alpha1 = T(0.1 / Square(s1 + eps));
249  T alpha2 = T(0.6 / Square(s2 + eps));
250  T alpha3 = T(0.3 / Square(s3 + eps));
251 
252  T sum = alpha1 + alpha2 + alpha3;
253 
254  return (alpha1 * phix1 + alpha2 * phix2 + alpha3 * phix3) / sum;
255 }
256 } // namespace CubbyFlow
257 
258 #endif
CubbyFlow::ENO3
std::array< T, 2 > ENO3(T *d0, T dx)
3rd order ENO. d0[3] is the origin.
Definition: PDE-Impl.hpp:46
CubbyFlow::Square
std::enable_if_t< std::is_arithmetic< T >::value, T > Square(T x)
Returns the square of x.
Definition: MathUtils-Impl.hpp:154
MathUtils.hpp
CubbyFlow::WENO5
std::array< T, 2 > WENO5(T *v, T h, T eps)
5th order WENO. d0[3] is the origin.
Definition: PDE-Impl.hpp:162
CubbyFlow::CD2
T CD2(T *d0, T dx)
2nd order central differencing. d0[1] is the origin.
Definition: PDE-Impl.hpp:39
CubbyFlow
Definition: pybind11Utils.hpp:20
CubbyFlow::Upwind1
std::array< T, 2 > Upwind1(T *d0, T dx)
1st order upwind differencing. d0[1] is the origin.
Definition: PDE-Impl.hpp:19
Generated on Wed Jul 27 2022 14:50:10 for CubbyFlow by   doxygen 1.8.13