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