PeriDyno 1.2.1
Loading...
Searching...
No Matches
NumericalScheme.h
Go to the documentation of this file.
1
16#pragma once
17
18#include "Vector/Vector4D.h"
19
20namespace dyno
21{
22#define POSITIVTY Real(0.1)
23
24 //A stable solution to calculate u and v
25 template<typename Real, typename Coord4D>
26 inline DYN_FUNC void ComputeVelocity(Real& u, Real& v, Coord4D gp)
27 {
28 Real EPSILON = 0.000001f;
29
30 Real h = maximum(gp.x, 0.0f);
31 Real h4 = h * h * h * h;
32 u = sqrtf(2.0f) * h * gp.y / (sqrtf(h4 + maximum(h4, EPSILON)));
33 v = sqrtf(2.0f) * h * gp.z / (sqrtf(h4 + maximum(h4, EPSILON)));
34 }
35
36 //A first-order upwind scheme to calculate -ghB'
37 template<typename Real, typename Coord4D>
38 DYN_FUNC Real FirstOrderUpwindPotential(Coord4D gpl, Coord4D gpr, Real GRAVITY, Real dt)
39 {
40 Real EPSILON = 0.000001f;
41
42 Real hl = maximum(gpl.x, 0.0f);
43 Real hl4 = hl * hl * hl * hl;
44 Real ul = sqrtf(2.0f) * hl * gpl.y / (sqrtf(hl4 + maximum(hl4, EPSILON)));
45 Real vl = sqrtf(2.0f) * hl * gpl.z / (sqrtf(hl4 + maximum(hl4, EPSILON)));
46
47 Real hr = maximum(gpr.x, 0.0f);
48 Real hr4 = hr * hr * hr * hr;
49 Real ur = sqrtf(2.0f) * hr * gpr.y / (sqrtf(hr4 + maximum(hr4, EPSILON)));
50 Real vr = sqrtf(2.0f) * hr * gpr.z / (sqrtf(hr4 + maximum(hr4, EPSILON)));
51
52 //Potential
53 Real wl = hl + gpl.w;
54 Real wr = hr + gpr.w;
55 Real potential = 0;
56 //Left
57 {
58 Real wl_diff = maximum(wl - gpr.w, Real(0));
59 potential -= GRAVITY * hl * wl_diff;
60 }
61 // Right
62 {
63 Real wr_diff = maximum(wr - gpl.w, Real(0));
64 potential += GRAVITY * hr * wr_diff;
65 }
66
67 return - dt * potential;
68 }
69
70 //A first-order upwind scheme to calculate (hu, huu, huv, 0)
71 template<typename Real, typename Coord4D>
72 DYN_FUNC Coord4D FirstOrderUpwindX(Coord4D gpl, Coord4D gpr, Real GRAVITY, Real dt)
73 {
74 Real EPSILON = 0.000001f;
75
76 Real hl = maximum(gpl.x, 0.0f);
77 Real hl4 = hl * hl * hl * hl;
78 Real ul = sqrtf(2.0f) * hl * gpl.y / (sqrtf(hl4 + maximum(hl4, EPSILON)));
79 Real vl = sqrtf(2.0f) * hl * gpl.z / (sqrtf(hl4 + maximum(hl4, EPSILON)));
80
81 Real hr = maximum(gpr.x, 0.0f);
82 Real hr4 = hr * hr * hr * hr;
83 Real ur = sqrtf(2.0f) * hr * gpr.y / (sqrtf(hr4 + maximum(hr4, EPSILON)));
84 Real vr = sqrtf(2.0f) * hr * gpr.z / (sqrtf(hr4 + maximum(hr4, EPSILON)));
85
86 Real hu = 0.5 * (hl * ul + hr * ur);
87 Real hm = 0.5f * (hl + hr);
88 Real hm4 = hm * hm * hm * hm;
89 Real um = sqrtf(2.0f) * hm * hu / (sqrtf(hm4 + maximum(hm4, EPSILON)));
90 Real vm = 0.5f * (vl + vr);
91
92 //Advection
93 Real h_hat = um >= 0 ? maximum(hl - maximum(gpr.w - gpl.w, 0.0f), 0.0f) : maximum(hr - maximum(gpl.w - gpr.w, 0.0f), 0.0f);
94 Real uh_hat = h_hat * um;// -dt * potential;
95 Real vh_hat = h_hat * vm;
96
97 Real alpha = (um >= 0 ? hl : hr) / maximum(dt * uh_hat, EPSILON);
98 Coord4D flux = minimum(alpha, POSITIVTY) * Coord4D(uh_hat, uh_hat * um, vh_hat * um, 0.0f);
99
100 return flux;
101 }
102
103 //A first-order upwind scheme to calculate (hv, huv, hvv, 0)
104 template<typename Real, typename Coord4D>
105 DYN_FUNC Coord4D FirstOrderUpwindY(Coord4D gpl, Coord4D gpr, Real GRAVITY, Real dt)
106 {
107 Real EPSILON = 0.000001f;
108
109 Real hl = maximum(gpl.x, 0.0f);
110 Real hl4 = hl * hl * hl * hl;
111 Real ul = sqrtf(2.0f) * hl * gpl.y / (sqrtf(hl4 + maximum(hl4, EPSILON)));
112 Real vl = sqrtf(2.0f) * hl * gpl.z / (sqrtf(hl4 + maximum(hl4, EPSILON)));
113
114 Real hr = maximum(gpr.x, 0.0f);
115 Real hr4 = hr * hr * hr * hr;
116 Real ur = sqrtf(2.0f) * hr * gpr.y / (sqrtf(hr4 + maximum(hr4, EPSILON)));
117 Real vr = sqrtf(2.0f) * hr * gpr.z / (sqrtf(hr4 + maximum(hr4, EPSILON)));
118
119 Real hv = 0.5 * (hl * vl + hr * vr);
120 Real hm = 0.5f * (hl + hr);
121 Real hm4 = hm * hm * hm * hm;
122 Real um = 0.5f * (ul + ur);
123 Real vm = sqrtf(2.0f) * hm * hv / (sqrtf(hm4 + maximum(hm4, EPSILON)));
124
125 //Advection
126 Real h_hat = vm >= 0 ? maximum(hl - maximum(gpr.w - gpl.w, 0.0f), 0.0f) : maximum(hr - maximum(gpl.w - gpr.w, 0.0f), 0.0f);
127 Real uh_hat = h_hat * um;
128 Real vh_hat = h_hat * vm;// -dt * potential;
129
130 Real alpha = (vm >= 0 ? hl : hr) / maximum(dt * vh_hat, EPSILON);
131 Coord4D flux = minimum(alpha, POSITIVTY) * Coord4D(vh_hat, uh_hat * vm, vh_hat * vm, 0.0f);
132
133 return flux;
134 }
135
136 template<typename Real, typename Coord4D>
137 DYN_FUNC Coord4D CentralUpwindX(Coord4D gpl, Coord4D gpr, Real GRAVITY)
138 {
139 Real EPSILON = 0.00001f;
140
141 Real h = maximum(0.5f * (gpl.x + gpr.x), 0.0f);
142 Real b = 0.5f * (gpl.w + gpr.w);
143
144 Real hl = maximum(gpl.x, 0.0f);
145 Real hl4 = hl * hl * hl * hl;
146 Real ul = sqrtf(2.0f) * hl * gpl.y / (sqrtf(hl4 + maximum(hl4, Real(EPSILON))));
147
148 Real hr = maximum(gpr.x, 0.0f);
149 Real hr4 = hr * hr * hr * hr;
150 Real ur = sqrtf(2.0f) * hr * gpr.y / (sqrtf(hr4 + maximum(hr4, Real(EPSILON))));
151
152 if (hl < EPSILON && hr < EPSILON)
153 {
154 return Coord4D(0.0f);
155 }
156
157 Real a_plus;
158 Real a_minus;
159 a_plus = maximum(maximum(Real(ul + sqrtf(GRAVITY * (gpl.x/*+gpl.w*/))), Real(ur + sqrtf(GRAVITY * (gpr.x/*+gpr.w*/)))), Real(0));
160 a_minus = minimum(minimum(Real(ul - sqrtf(GRAVITY * (gpl.x/*+gpl.w*/))), Real(ur - sqrtf(GRAVITY * (gpr.x/*+gpr.w*/)))), Real(0));
161
162 Coord4D delta_U = gpr - gpl;
163 if (gpl.x > EPSILON && gpr.x > EPSILON) {
164 delta_U.x += delta_U.w;
165 }
166
167 delta_U.w = 0.0f;
168
169 Coord4D Fl = Coord4D(gpl.y, gpl.y * ul, gpl.z * ul, 0.0f);
170 Coord4D Fr = Coord4D(gpr.y, gpr.y * ur, gpr.z * ur, 0.0f);
171
172 Coord4D re = (a_plus * Fl - a_minus * Fr) / (a_plus - a_minus + EPSILON) + a_plus * a_minus / (a_plus - a_minus + EPSILON) * delta_U;
173
174 if (ul == 0 && ur == 0) {//abs(ul) <EPSILON && abs(ur) <EPSILON
175 re.x = 0;
176 re.y = 0;
177 re.z = 0;
178 }
179
180 return re;
181 }
182
183 template<typename Real, typename Coord4D>
184 DYN_FUNC Coord4D CentralUpwindY(Coord4D gpl, Coord4D gpr, Real GRAVITY)
185 {
186 Real EPSILON = 0.00001f;
187
188 Real hl = maximum(gpl.x, 0.0f);
189 Real hl4 = hl * hl * hl * hl;
190 Real vl = sqrtf(2.0f) * hl * gpl.z / (sqrtf(hl4 + max(hl4, EPSILON)));
191
192 Real hr = maximum(gpr.x, 0.0f);
193 Real hr4 = hr * hr * hr * hr;
194 Real vr = sqrtf(2.0f) * hr * gpr.z / (sqrtf(hr4 + max(hr4, EPSILON)));
195
196 if (hl < EPSILON && hr < EPSILON)
197 {
198 return Coord4D(0.0f);
199 }
200
201 Real a_plus = maximum(maximum(Real(vl + sqrtf(GRAVITY * (gpl.x/* + gpl.w*/))), Real(vr + sqrtf(GRAVITY * (gpr.x/* + gpr.w*/)))), Real(0));
202 Real a_minus = minimum(minimum(Real(vl - sqrtf(GRAVITY * (gpl.x/* + gpl.w*/))), Real(vr - sqrtf(GRAVITY * (gpr.x/* + gpr.w*/)))), Real(0));
203
204 Real b = 0.5f * (gpl.w + gpr.w);
205
206 Coord4D delta_U = gpr - gpl;
207 if (gpl.x > EPSILON && gpr.x > EPSILON)
208 {
209 delta_U.x += delta_U.w;
210 }
211 delta_U.w = 0.0f;
212
213 Coord4D Fl = Coord4D(gpl.z, gpl.y * vl, gpl.z * vl, 0.0f);
214 Coord4D Fr = Coord4D(gpr.z, gpr.y * vr, gpr.z * vr, 0.0f);
215
216 Coord4D re = (a_plus * Fl - a_minus * Fr) / (a_plus - a_minus) + a_plus * a_minus / (a_plus - a_minus) * delta_U;
217
218 if (vl == 0 && vr == 0)
219 {
220 re.x = 0;
221 re.y = 0;
222 re.z = 0;
223 }
224 return re;
225 }
226}
POSITIVTY
#define POSITIVTY
Definition NumericalScheme.h:22
Real
double Real
Definition Typedef.inl:23
dyno
This is an implementation of AdditiveCCD based on peridyno.
Definition Array.h:25
dyno::FirstOrderUpwindX
DYN_FUNC Coord4D FirstOrderUpwindX(Coord4D gpl, Coord4D gpr, Real GRAVITY, Real dt)
Definition NumericalScheme.h:72
dyno::ComputeVelocity
DYN_FUNC void ComputeVelocity(Real &u, Real &v, Coord4D gp)
Definition NumericalScheme.h:26
dyno::EPSILON
constexpr Real EPSILON
Definition Typedef.inl:42
dyno::FirstOrderUpwindPotential
DYN_FUNC Real FirstOrderUpwindPotential(Coord4D gpl, Coord4D gpr, Real GRAVITY, Real dt)
Definition NumericalScheme.h:38
dyno::CentralUpwindY
DYN_FUNC Coord4D CentralUpwindY(Coord4D gpl, Coord4D gpr, Real GRAVITY)
Definition NumericalScheme.h:184
dyno::minimum
DYN_FUNC T minimum(const T &v0, const T &v1)
Definition SimpleMath.h:120
dyno::FirstOrderUpwindY
DYN_FUNC Coord4D FirstOrderUpwindY(Coord4D gpl, Coord4D gpr, Real GRAVITY, Real dt)
Definition NumericalScheme.h:105
dyno::maximum
DYN_FUNC T maximum(const T &v0, const T &v1)
Definition SimpleMath.h:160
dyno::CentralUpwindX
DYN_FUNC Coord4D CentralUpwindX(Coord4D gpl, Coord4D gpr, Real GRAVITY)
Definition NumericalScheme.h:137
max
#define max(x, y)
Definition svd3_cuda.h:41
Generated by doxygen 1.13.2