29 const amrex::Real tol,
33 const amrex::Real omega = 2.0 * M_PI / T;
34 const amrex::Real depsdk = H / 2.0;
39 amrex::Real k = (omega * omega) / g;
41 amrex::Real f = tol + 1.0;
42 while (std::abs(f) > tol && iter < iter_max) {
47 2. * std::exp(2. * k * d) / (std::exp(4. * k * d) + 1.);
48 const amrex::Real C = 1.0 - S;
49 const amrex::Real eps = k * H / 2.0;
50 const amrex::Real C0 = std::sqrt(std::tanh(k * d));
51 const amrex::Real C2 = C0 * (2.0 + 7.0 * S * S) / (4.0 * C * C);
52 const amrex::Real numterm_C4 =
53 (4.0 + 32.0 * S - 116.0 * std::pow(S, 2) - 400.0 * std::pow(S, 3) -
54 71.0 * std::pow(S, 4) + 146.0 * std::pow(S, 5));
55 const amrex::Real C4 = C0 * numterm_C4 / (32.0 * std::pow(C, 5));
57 const amrex::Real dSdk = -2.0 * d * std::sinh(2.0 * k * d) /
58 std::pow(std::cosh(2.0 * k * d), 2);
59 const amrex::Real dCdk = -dSdk;
60 const amrex::Real dC0dk =
61 d / (2.0 * C0 * std::pow(std::cosh(k * d), 2));
63 const amrex::Real dC2dk =
65 (dC0dk * (2.0 + 7.0 * std::pow(S, 2)) + C0 * 14.0 * S * dSdk) -
66 C0 * (2.0 + 7.0 * std::pow(S, 2)) * 8.0 * C * dCdk) /
67 (16.0 * std::pow(C, 4));
68 const amrex::Real dC4dk =
69 (32.0 * std::pow(C, 5) *
71 C0 * (32.0 * dSdk - 232.0 * S * dSdk -
72 1200 * std::pow(S, 2) * dSdk -
73 284.0 * std::pow(S, 3) * dSdk + 730 * std::pow(S, 4))) -
74 C0 * numterm_C4 * 160.0 * std::pow(C, 4) * dCdk) /
75 (1024.0 * std::pow(C, 10));
78 amrex::Real dfdk = g * std::pow(C0, 2) + g * k * 2.0 * C0 * dC0dk;
82 g * (2.0 * C0 * std::pow(eps, 2) * C2 +
83 std::pow(eps, 4) * std::pow(C2, 2)) +
85 (2.0 * dC0dk * std::pow(eps, 2) * C2 +
87 (2.0 * eps * depsdk * C2 + std::pow(eps, 2) * dC2dk) +
88 4.0 * std::pow(eps, 3) * std::pow(C2, 2) * depsdk +
89 std::pow(eps, 4) * 2.0 * C2 * dC2dk);
92 dfdk += g * (2.0 * std::pow(eps, 4) * C0 * C4 +
93 2.0 * std::pow(eps, 6) * C2 * C4 +
94 std::pow(eps, 8) * std::pow(C4, 2)) +
96 (8.0 * std::pow(eps, 3) * depsdk * C0 * C4 +
97 2.0 * std::pow(eps, 4) * (C0 * dC4dk + C4 * dC0dk) +
98 12 * std::pow(eps, 5) * depsdk * C2 * C4 +
99 2.0 * std::pow(eps, 6) * (C2 * dC4dk + C4 * dC2dk) +
100 8.0 * std::pow(eps, 7) * depsdk * std::pow(C4, 2) +
101 std::pow(eps, 8) * 2.0 * C4 * dC4dk);
106 f = g * k * std::pow(C0, 2);
110 (2.0 * C0 * std::pow(eps, 2) * C2 +
111 std::pow(eps, 4) * std::pow(C2, 2));
115 (2.0 * std::pow(eps, 4) * C0 * C4 +
116 2.0 * std::pow(eps, 6) * C2 * C4 +
117 std::pow(eps, 8) * std::pow(C4, 2));
130 if (iter == iter_max) {
134 return 2.0 * M_PI / k;
139 amrex::Real wavenumber,
140 amrex::Real water_depth,
161 amrex::Real kd = wavenumber * water_depth;
162 if (kd > 50. * M_PI) {
167 amrex::Real S = 2. * std::exp(2. * kd) / (std::exp(4. * kd) + 1.);
168 amrex::Real C = 1.0 - S;
169 amrex::Real Sh = std::sinh(kd);
170 amrex::Real Th = std::tanh(kd);
172 amrex::Real CTh = (1. + std::exp(-2 * kd)) / (1. - std::exp(-2 * kd));
175 a11 = 1. / std::sinh(kd);
177 a22 = 3. * std::pow(S, 2) / (2 * std::pow(C, 2));
178 b22 = CTh * (1 + 2. * S) / (2 * C);
179 c2 = std::sqrt(Th) * (2 + 7 * std::pow(S, 2)) / (4 * std::pow(C, 2));
180 if (stokes_order == 2) {
185 a31 = (-4 - 20 * S + 10 * std::pow(S, 2) - 13 * std::pow(S, 3)) /
186 (8 * Sh * std::pow(C, 3));
188 (-2 * std::pow(S, 2) + 11 * std::pow(S, 3)) / (8 * Sh * std::pow(C, 3));
189 b31 = -3 * (1 + 3 * S + 3 * std::pow(S, 2) + 2 * std::pow(S, 3)) /
190 (8 * std::pow(C, 3));
191 if (stokes_order == 3) {
196 a42 = (12 * S - 14 * std::pow(S, 2) - 264 * std::pow(S, 3) -
197 45 * std::pow(S, 4) - 13 * std::pow(S, 5)) /
198 (24 * std::pow(C, 5));
199 a44 = (10 * std::pow(S, 3) - 174 * std::pow(S, 4) + 291 * std::pow(S, 5) +
200 278 * std::pow(S, 6)) /
201 (48 * (3 + 2 * S) * std::pow(C, 5));
203 (6 - 26 * S - 182 * std::pow(S, 2) - 204 * std::pow(S, 3) -
204 25 * std::pow(S, 4) + 26 * std::pow(S, 5)) /
205 (6 * (3 + 2 * S) * std::pow(C, 4));
207 (24 + 92 * S + 122 * std::pow(S, 2) + 66 * std::pow(S, 3) +
208 67 * std::pow(S, 4) + 34 * std::pow(S, 5)) /
209 (24 * (3 + 2 * S) * std::pow(C, 4));
211 (4 + 32 * S - 116 * std::pow(S, 2) - 400 * std::pow(S, 3) -
212 71 * std::pow(S, 4) + 146 * std::pow(S, 5)) /
213 (32 * std::pow(C, 5));
214 if (stokes_order == 4) {
220 (-1184 + 32 * S + 13232 * std::pow(S, 2) + 21712 * std::pow(S, 3) +
221 20940 * std::pow(S, 4) + 12554 * std::pow(S, 5) -
222 500 * std::pow(S, 6) - 3341 * std::pow(S, 7) - 670 * std::pow(S, 8)) /
223 (64 * Sh * (3 + 2 * S) * (4 + S) * std::pow(C, 6));
224 a53 = (4 * S + 105 * pow(S, 2) + 198 * std::pow(S, 3) -
225 1376 * std::pow(S, 4) - 1302 * std::pow(S, 5) -
226 117 * std::pow(S, 6) + 58 * std::pow(S, 7)) /
227 (32 * Sh * (3 + 2 * S) * std::pow(C, 6));
229 (-6 * std::pow(S, 3) + 272 * std::pow(S, 4) - 1552 * std::pow(S, 5) +
230 852 * std::pow(S, 6) + 2029 * std::pow(S, 7) + 430 * std::pow(S, 8)) /
231 (64 * Sh * (3 + 2 * S) * (4 + S) * std::pow(C, 6));
233 (132 + 17 * S - 2216 * std::pow(S, 2) - 5897 * std::pow(S, 3) -
234 6292 * std::pow(S, 4) - 2687 * std::pow(S, 5) +
235 194 * std::pow(S, 6) + 467 * std::pow(S, 7) + 82 * std::pow(S, 8)) /
236 (128 * (3 + 2 * S) * (4 + S) * std::pow(C, 6));
239 (300 + 1579 * S + 3176 * std::pow(S, 2) + 2949 * std::pow(S, 3) +
240 1188 * std::pow(S, 4) + 675 * std::pow(S, 5) + 1326 * std::pow(S, 6) +
241 827 * std::pow(S, 7) + 130 * std::pow(S, 8)) /
242 (384 * (3 + 2 * S) * (4 + S) * std::pow(C, 6));
243 if (stokes_order == 5) {
247 if (stokes_order > 5 || stokes_order < 2) {
249 "invalid stokes order specified. It should be between 2,3,4 or 5 ");
256 amrex::Real wavelength,
257 amrex::Real water_depth,
258 amrex::Real wave_height,
264 amrex::Real phase_offset,
270 const amrex::Real wavenumber = 2.0 * M_PI / wavelength;
274 amrex::Real a11{0.0}, a22{0.0}, b22{0.0}, c2{0.0};
275 amrex::Real a31{0.0}, a33{0.0}, b31{0.0}, a42{0.0}, a44{0.0};
276 amrex::Real b42{0.0}, b44{0.0}, c4{0.0};
277 amrex::Real a51{0.0}, a53{0.0}, a55{0.0}, b53{0.0}, b55{0.0};
280 stokes_order, wavenumber, water_depth, c0, a11, a22, b22, c2, a31, a33,
281 b31, a42, a44, b42, b44, c4, a51, a53, a55, b53, b55);
283 const amrex::Real eps = wavenumber * wave_height / 2.;
284 const amrex::Real c = (c0 + std::pow(eps, 2) * c2 + std::pow(eps, 4) * c4) *
285 std::sqrt(g / wavenumber);
287 const amrex::Real omega = c * wavenumber;
288 const amrex::Real phase = wavenumber * x - omega * time - phase_offset;
290 eta = (eps * std::cos(phase)
291 + std::pow(eps, 2) * b22 * std::cos(2. * phase)
292 + std::pow(eps, 3) * b31 * (std::cos(phase) - std::cos(3. * phase)) +
294 (b42 * std::cos(2. * phase) + b44 * std::cos(4. * phase)) +
296 (-(b53 + b55) * std::cos(phase) + b53 * std::cos(3 * phase) +
297 b55 * std::cos(5 * phase))) /
305 const int MAX_ORDER = 5;
306 amrex::GpuArray<amrex::Real, MAX_ORDER> a;
307 if (stokes_order == 2) {
311 if (stokes_order == 3) {
312 a[0] = a11 + (eps * eps) * a31;
316 if (stokes_order == 4) {
317 a[0] = a11 + (eps * eps) * a31;
318 a[1] = a22 + (eps * eps) * a42;
322 if (stokes_order == 5) {
323 a[0] = a11 + (eps * eps) * a31 + std::pow(eps, 4) * a51;
324 a[1] = a22 + (eps * eps) * a42;
325 a[2] = a33 + (eps * eps) * a53;
333 for (
int n = 1; n <= stokes_order; ++n) {
336 amrex::Real nkdz = n * wavenumber * (water_depth + (z - zsl));
337 if (nkdz > 50. * M_PI) {
340 u_w += std::pow(eps, n) * n * a[n - 1] * std::cosh(nkdz) *
342 w_w += std::pow(eps, n) * n * a[n - 1] * std::sinh(nkdz) *
345 u_w *= (c0 * std::sqrt(g / wavenumber));
346 w_w *= (c0 * std::sqrt(g / wavenumber));
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void stokes_coefficients(int stokes_order, amrex::Real wavenumber, amrex::Real water_depth, amrex::Real &c0, amrex::Real &a11, amrex::Real &a22, amrex::Real &b22, amrex::Real &c2, amrex::Real &a31, amrex::Real &a33, amrex::Real &b31, amrex::Real &a42, amrex::Real &a44, amrex::Real &b42, amrex::Real &b44, amrex::Real &c4, amrex::Real &a51, amrex::Real &a53, amrex::Real &a55, amrex::Real &b53, amrex::Real &b55)
Definition stokes_waves_K.H:137
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void stokes_waves(int stokes_order, amrex::Real wavelength, amrex::Real water_depth, amrex::Real wave_height, amrex::Real zsl, amrex::Real g, amrex::Real x, amrex::Real z, amrex::Real time, amrex::Real phase_offset, amrex::Real &eta, amrex::Real &u_w, amrex::Real &v_w, amrex::Real &w_w)
Definition stokes_waves_K.H:254