/home/runner/work/kynema-sgf/kynema-sgf/src/ocean_waves/relaxation_zones/stokes_waves_K.H Source File

Kynema-SGF API: /home/runner/work/kynema-sgf/kynema-sgf/src/ocean_waves/relaxation_zones/stokes_waves_K.H Source File
Kynema-SGF API v0.1.0
CFD solver for wind plant simulations
Loading...
Searching...
No Matches
stokes_waves_K.H
Go to the documentation of this file.
1#ifndef STOKES_WAVES_K_H
2#define STOKES_WAVES_K_H
3
4#include <algorithm>
5#include <cmath>
6#include <numbers>
7#include "AMReX_FArrayBox.H"
8#include "src/utilities/math_ops.H"
9
10using namespace amrex::literals;
11
12// Stokes waves theory adapted from
13// Fenton, J., Fifth Order Stokes Theory for Steady Waves
14// Journal of Waterway, Port, Coastal and Ocean Engineering, 1985, 111, 216-234
15
16// Updated Table 1 from Fenton 1985 found in
17// Fenton, J., Nonlinear Wave Theories
18// The Sea Vol. 9 Ocean Engineering Science, 1990
19// https://johndfenton.com/Papers/Fenton90b-Nonlinear-wave-theories.pdf
20
21// Relevant results are summarized in
22// Kinnas S. A., Notes on fifth order gravity wave theory
23// https://www.caee.utexas.edu/prof/kinnas/ce358/oenotes/kinnas_stokes11.pdf
24
25namespace kynema_sgf::ocean_waves::relaxation_zones {
26
27// Compute wavelength as a function of wave period, water depth, and g
28AMREX_FORCE_INLINE amrex::Real stokes_wave_length(
29 const amrex::Real T,
30 const amrex::Real d,
31 const amrex::Real H,
32 const int order,
33 const amrex::Real g,
34 const amrex::Real tol,
35 const int iter_max)
36{
37 // Calculate constants and derivatives that do not change with iteration
38 const amrex::Real omega = 2.0_rt * std::numbers::pi_v<amrex::Real> / T;
39 const amrex::Real depsdk = H / 2.0_rt;
40
41 // Begin Newton-Raphson iterations
42 int iter = 0;
43 // First guess is first-order
44 amrex::Real k = (omega * omega) / g;
45 // Cannot skip loop
46 amrex::Real f = tol + 1.0_rt;
47 while (std::abs(f) > tol && iter < iter_max) {
48 // Calculate current constants
49
50 // Exponential definition of S = sech(2kd)
51 const amrex::Real S = 2.0_rt * std::exp(2.0_rt * k * d) /
52 (std::exp(4.0_rt * k * d) + 1.0_rt);
53 const amrex::Real C = 1.0_rt - S;
54 const amrex::Real eps = k * H / 2.0_rt;
55 const amrex::Real C0 = std::sqrt(std::tanh(k * d));
56 const amrex::Real C2 =
57 C0 * (2.0_rt + 7.0_rt * S * S) / (4.0_rt * C * C);
58 const amrex::Real numterm_C4 =
59 (4.0_rt + (32.0_rt * S) - (116.0_rt * utils::powi(S, 2)) -
60 (400.0_rt * utils::powi(S, 3)) - (71.0_rt * utils::powi(S, 4)) +
61 (146.0_rt * utils::powi(S, 5)));
62 const amrex::Real C4 = C0 * numterm_C4 / (32.0_rt * utils::powi(C, 5));
63 // Calculate pure derivates
64 const amrex::Real dSdk = -2.0_rt * d * std::sinh(2.0_rt * k * d) /
65 utils::powi(std::cosh(2.0_rt * k * d), 2);
66 const amrex::Real dCdk = -dSdk;
67 const amrex::Real dC0dk =
68 d / (2.0_rt * C0 * utils::powi(std::cosh(k * d), 2));
69 // Calculate derivatives with products
70 const amrex::Real dC2dk =
71 (4.0_rt * utils::powi(C, 2) *
72 (dC0dk * (2.0_rt + 7.0_rt * utils::powi(S, 2)) +
73 C0 * 14.0_rt * S * dSdk) -
74 C0 * (2.0_rt + 7.0_rt * utils::powi(S, 2)) * 8.0_rt * C * dCdk) /
75 (16.0_rt * utils::powi(C, 4));
76 const amrex::Real dC4dk =
77 (32.0_rt * utils::powi(C, 5) *
78 (dC0dk * numterm_C4 +
79 C0 * (32.0_rt * dSdk - 232.0_rt * S * dSdk -
80 1200.0_rt * utils::powi(S, 2) * dSdk -
81 284.0_rt * utils::powi(S, 3) * dSdk +
82 730.0_rt * utils::powi(S, 4))) -
83 C0 * numterm_C4 * 160.0_rt * utils::powi(C, 4) * dCdk) /
84 (1024.0_rt * utils::powi(C, 10));
85
86 // Calculate derivative for loop convergence
87 amrex::Real dfdk =
88 (g * utils::powi(C0, 2)) + (g * k * 2.0_rt * C0 * dC0dk);
89 // Add additional terms depending on order
90 if (order >= 2) {
91 dfdk +=
92 (g * (2.0_rt * C0 * utils::powi(eps, 2) * C2 +
93 utils::powi(eps, 4) * utils::powi(C2, 2))) +
94 (g * k *
95 (2.0_rt * dC0dk * utils::powi(eps, 2) * C2 +
96 2.0_rt * C0 *
97 (2.0_rt * eps * depsdk * C2 +
98 utils::powi(eps, 2) * dC2dk) +
99 4.0_rt * utils::powi(eps, 3) * utils::powi(C2, 2) * depsdk +
100 utils::powi(eps, 4) * 2.0_rt * C2 * dC2dk));
101 }
102 if (order >= 4) {
103 dfdk +=
104 (g * (2.0_rt * utils::powi(eps, 4) * C0 * C4 +
105 2.0_rt * utils::powi(eps, 6) * C2 * C4 +
106 utils::powi(eps, 8) * utils::powi(C4, 2))) +
107 (g * k *
108 (8.0_rt * utils::powi(eps, 3) * depsdk * C0 * C4 +
109 2.0_rt * utils::powi(eps, 4) * (C0 * dC4dk + C4 * dC0dk) +
110 12.0_rt * utils::powi(eps, 5) * depsdk * C2 * C4 +
111 2.0_rt * utils::powi(eps, 6) * (C2 * dC4dk + C4 * dC2dk) +
112 8.0_rt * utils::powi(eps, 7) * depsdk * utils::powi(C4, 2) +
113 utils::powi(eps, 8) * 2.0_rt * C4 * dC4dk));
114 }
115 k = k - (f / dfdk);
116
117 iter += 1;
118 f = g * k * utils::powi(C0, 2);
119 // Add additional terms depending on order
120 if (order >= 2) {
121 f += g * k *
122 (2.0_rt * C0 * utils::powi(eps, 2) * C2 +
123 utils::powi(eps, 4) * utils::powi(C2, 2));
124 }
125 if (order >= 4) {
126 f += g * k *
127 (2.0_rt * utils::powi(eps, 4) * C0 * C4 +
128 2.0_rt * utils::powi(eps, 6) * C2 * C4 +
129 utils::powi(eps, 8) * utils::powi(C4, 2));
130 }
131 // Subtract omega^2 to measure convergence
132 f -= omega * omega;
133 }
134
135 if (k < tol) {
136 // Return negative wavelength if faulty result
137 return -1;
138 }
139 if (std::isnan(k)) {
140 return -2;
141 }
142 if (iter == iter_max) {
143 return -3;
144 }
145 // Return wavelength calculated from wavenumber
146 return 2.0_rt * std::numbers::pi_v<amrex::Real> / k;
147}
148
149AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void stokes_coefficients(
150 int stokes_order,
151 amrex::Real wavenumber,
152 amrex::Real water_depth,
153 amrex::Real& c0,
154 amrex::Real& a11,
155 amrex::Real& a22,
156 amrex::Real& b22,
157 amrex::Real& c2,
158 amrex::Real& a31,
159 amrex::Real& a33,
160 amrex::Real& b31,
161 amrex::Real& a42,
162 amrex::Real& a44,
163 amrex::Real& b42,
164 amrex::Real& b44,
165 amrex::Real& c4,
166 amrex::Real& a51,
167 amrex::Real& a53,
168 amrex::Real& a55,
169 amrex::Real& b53,
170 amrex::Real& b55)
171{
172
173 amrex::Real kd = wavenumber * water_depth;
174 kd = amrex::min<amrex::Real>(kd, 50.0_rt * std::numbers::pi_v<amrex::Real>);
175
176 // Exponential definition of S = sech(2kd)
177 amrex::Real S =
178 2.0_rt * std::exp(2.0_rt * kd) / (std::exp(4.0_rt * kd) + 1.0_rt);
179 amrex::Real C = 1.0_rt - S;
180 amrex::Real Sh = std::sinh(kd);
181 amrex::Real Th = std::tanh(kd);
182 // Exponential definition of coth(kd)
183 amrex::Real CTh =
184 (1.0_rt + std::exp(-2.0_rt * kd)) / (1.0_rt - std::exp(-2 * kd));
185
186 c0 = std::sqrt(Th);
187 a11 = 1.0_rt / std::sinh(kd); // Hyperbolic cosecant
188 // Second order coefficients
189 a22 = 3.0_rt * utils::powi(S, 2) / (2.0_rt * utils::powi(C, 2));
190 b22 = CTh * (1.0_rt + 2.0_rt * S) / (2.0_rt * C);
191 c2 = std::sqrt(Th) * (2.0_rt + 7.0_rt * utils::powi(S, 2)) /
192 (4.0_rt * utils::powi(C, 2));
193 if (stokes_order == 2) {
194 return;
195 }
196
197 // Third order coefficients
198 a31 = (-4.0_rt - 20.0_rt * S + 10.0_rt * utils::powi(S, 2) -
199 13.0_rt * utils::powi(S, 3)) /
200 (8.0_rt * Sh * utils::powi(C, 3));
201 a33 = (-2.0_rt * utils::powi(S, 2) + 11.0_rt * utils::powi(S, 3)) /
202 (8.0_rt * Sh * utils::powi(C, 3));
203 b31 = -3.0_rt *
204 (1.0_rt + 3.0_rt * S + 3.0_rt * utils::powi(S, 2) +
205 2.0_rt * utils::powi(S, 3)) /
206 (8.0_rt * utils::powi(C, 3));
207 if (stokes_order == 3) {
208 return;
209 }
210
211 // Fourth order coefficients
212 a42 = (12.0_rt * S - 14.0_rt * utils::powi(S, 2) -
213 264.0_rt * utils::powi(S, 3) - 45.0_rt * utils::powi(S, 4) -
214 13.0_rt * utils::powi(S, 5)) /
215 (24.0_rt * utils::powi(C, 5));
216 a44 = (10.0_rt * utils::powi(S, 3) - 174.0_rt * utils::powi(S, 4) +
217 291.0_rt * utils::powi(S, 5) + 278.0_rt * utils::powi(S, 6)) /
218 (48.0_rt * (3.0_rt + 2.0_rt * S) * utils::powi(C, 5));
219 b42 = CTh *
220 (6.0_rt - 26.0_rt * S - 182.0_rt * utils::powi(S, 2) -
221 204.0_rt * utils::powi(S, 3) - 25.0_rt * utils::powi(S, 4) +
222 26.0_rt * utils::powi(S, 5)) /
223 (6.0_rt * (3.0_rt + 2.0_rt * S) * utils::powi(C, 4));
224 b44 = CTh *
225 (24.0_rt + 92.0_rt * S + 122.0_rt * utils::powi(S, 2) +
226 66.0_rt * utils::powi(S, 3) + 67.0_rt * utils::powi(S, 4) +
227 34.0_rt * utils::powi(S, 5)) /
228 (24.0_rt * (3.0_rt + 2.0_rt * S) * utils::powi(C, 4));
229 c4 = std::sqrt(Th) *
230 (4.0_rt + 32.0_rt * S - 116.0_rt * utils::powi(S, 2) -
231 400.0_rt * utils::powi(S, 3) - 71.0_rt * utils::powi(S, 4) +
232 146.0_rt * utils::powi(S, 5)) /
233 (32.0_rt * utils::powi(C, 5));
234 if (stokes_order == 4) {
235 return;
236 }
237
238 // Fifth order coefficients
239 a51 = (-1184.0_rt + 32.0_rt * S + 13232.0_rt * utils::powi(S, 2) +
240 21712.0_rt * utils::powi(S, 3) + 20940.0_rt * utils::powi(S, 4) +
241 12554.0_rt * utils::powi(S, 5) - 500.0_rt * utils::powi(S, 6) -
242 3341.0_rt * utils::powi(S, 7) - 670.0_rt * utils::powi(S, 8)) /
243 (64.0_rt * Sh * (3.0_rt + 2.0_rt * S) * (4.0_rt + S) *
244 utils::powi(C, 6));
245 a53 = (4.0_rt * S + 105.0_rt * utils::powi(S, 2) +
246 198.0_rt * utils::powi(S, 3) - 1376.0_rt * utils::powi(S, 4) -
247 1302.0_rt * utils::powi(S, 5) - 117.0_rt * utils::powi(S, 6) +
248 58.0_rt * utils::powi(S, 7)) /
249 (32.0_rt * Sh * (3.0_rt + 2.0_rt * S) * utils::powi(C, 6));
250 a55 = (-6.0_rt * utils::powi(S, 3) + 272.0_rt * utils::powi(S, 4) -
251 1552.0_rt * utils::powi(S, 5) + 852.0_rt * utils::powi(S, 6) +
252 2029.0_rt * utils::powi(S, 7) + 430.0_rt * utils::powi(S, 8)) /
253 (64.0_rt * Sh * (3.0_rt + 2.0_rt * S) * (4.0_rt + S) *
254 utils::powi(C, 6));
255 b53 = 9.0_rt *
256 (132.0_rt + 17.0_rt * S - 2216.0_rt * utils::powi(S, 2) -
257 5897.0_rt * utils::powi(S, 3) - 6292.0_rt * utils::powi(S, 4) -
258 2687.0_rt * utils::powi(S, 5) + 194.0_rt * utils::powi(S, 6) +
259 467.0_rt * utils::powi(S, 7) + 82.0_rt * utils::powi(S, 8)) /
260 (128.0_rt * (3.0_rt + 2.0_rt * S) * (4.0_rt + S) * utils::powi(C, 6));
261 b55 = 5.0_rt *
262 (300.0_rt + 1579.0_rt * S + 3176.0_rt * utils::powi(S, 2) +
263 2949.0_rt * utils::powi(S, 3) + 1188.0_rt * utils::powi(S, 4) +
264 675.0_rt * utils::powi(S, 5) + 1326.0_rt * utils::powi(S, 6) +
265 827.0_rt * utils::powi(S, 7) + 130.0_rt * utils::powi(S, 8)) /
266 (384.0_rt * (3.0_rt + 2.0_rt * S) * (4.0_rt + S) * utils::powi(C, 6));
267 if (stokes_order == 5) {
268 return;
269 }
270
271 if (stokes_order > 5 || stokes_order < 2) {
272 amrex::Abort(
273 "invalid stokes order specified. It should be between 2,3,4 or 5 ");
274 }
275}
276
277// Based on Fenton 1985
278AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void stokes_waves(
279 int stokes_order,
280 amrex::Real wavelength,
281 amrex::Real water_depth,
282 amrex::Real wave_height,
283 amrex::Real zsl,
284 amrex::Real g,
285 amrex::Real x,
286 amrex::Real z,
287 amrex::Real time,
288 amrex::Real phase_offset,
289 amrex::Real hyp_lo_limit,
290 amrex::Real hyp_hi_limit,
291 amrex::Real& eta,
292 amrex::Real& u_w,
293 amrex::Real& v_w,
294 amrex::Real& w_w)
295{
296 const amrex::Real wavenumber =
297 2.0_rt * std::numbers::pi_v<amrex::Real> / wavelength;
298
299 // some parameters
300 amrex::Real c0{0.0_rt};
301 amrex::Real a11{0.0_rt}, a22{0.0_rt}, b22{0.0_rt}, c2{0.0_rt};
302 amrex::Real a31{0.0_rt}, a33{0.0_rt}, b31{0.0_rt}, a42{0.0_rt}, a44{0.0_rt};
303 amrex::Real b42{0.0_rt}, b44{0.0_rt}, c4{0.0_rt};
304 amrex::Real a51{0.0_rt}, a53{0.0_rt}, a55{0.0_rt}, b53{0.0_rt}, b55{0.0_rt};
305
306 stokes_coefficients(
307 stokes_order, wavenumber, water_depth, c0, a11, a22, b22, c2, a31, a33,
308 b31, a42, a44, b42, b44, c4, a51, a53, a55, b53, b55);
309
310 const amrex::Real eps = wavenumber * wave_height / 2.0_rt; // Steepness (ka)
311 const amrex::Real c =
312 (c0 + utils::powi(eps, 2) * c2 + utils::powi(eps, 4) * c4) *
313 std::sqrt(g / wavenumber);
314
315 const amrex::Real omega = c * wavenumber;
316 const amrex::Real phase = (wavenumber * x) - (omega * time) - phase_offset;
317
318 eta = ((eps * std::cos(phase) // first order term
319 + utils::powi(eps, 2) * b22 *
320 std::cos(2.0_rt * phase) // second order term
321 + utils::powi(eps, 3) * b31 *
322 (std::cos(phase) - std::cos(3.0_rt * phase)) +
323 utils::powi(eps, 4) * (b42 * std::cos(2.0_rt * phase) +
324 b44 * std::cos(4.0_rt * phase)) +
325 utils::powi(eps, 5) * (-(b53 + b55) * std::cos(phase) +
326 b53 * std::cos(3.0_rt * phase) +
327 b55 * std::cos(5.0_rt * phase))) /
328 wavenumber) +
329 zsl;
330
331 // Compute velocities components using Eq.(21) Eq.(23) from Kinnas
332 // https://www.caee.utexas.edu/prof/kinnas/ce358/oenotes/kinnas_stokes11.pdf
333 // Define coefficients using Eq.(19)
334
335 const int MAX_ORDER = 5;
336 amrex::GpuArray<amrex::Real, MAX_ORDER> a;
337 if (stokes_order == 2) {
338 a[0] = a11;
339 a[1] = a22;
340 }
341 if (stokes_order == 3) {
342 a[0] = a11 + ((eps * eps) * a31);
343 a[1] = a22;
344 a[2] = a33;
345 }
346 if (stokes_order == 4) {
347 a[0] = a11 + ((eps * eps) * a31);
348 a[1] = a22 + ((eps * eps) * a42);
349 a[2] = a33;
350 a[3] = a44;
351 }
352 if (stokes_order == 5) {
353 a[0] = a11 + ((eps * eps) * a31) + (utils::powi(eps, 4) * a51);
354 a[1] = a22 + ((eps * eps) * a42);
355 a[2] = a33 + ((eps * eps) * a53);
356 a[3] = a44;
357 a[4] = a55;
358 }
359
360 u_w = 0.0_rt;
361 v_w = 0.0_rt;
362 w_w = 0.0_rt;
363 for (int n = 1; n <= stokes_order; ++n) {
364 // Ensure finite values of velocity for large kd
365 const amrex::Real nkdz = n * wavenumber * (water_depth + (z - zsl));
366 amrex::Real cosh_nkdz = std::cosh(nkdz);
367 amrex::Real sinh_nkdz = std::sinh(nkdz);
368 cosh_nkdz = amrex::min<amrex::Real>(
369 hyp_hi_limit, amrex::max<amrex::Real>(hyp_lo_limit, cosh_nkdz));
370 sinh_nkdz = amrex::min<amrex::Real>(
371 hyp_hi_limit, amrex::max<amrex::Real>(hyp_lo_limit, sinh_nkdz));
372 // Contributions to velocity
373 u_w += utils::powi(eps, n) * static_cast<amrex::Real>(n) * a[n - 1] *
374 cosh_nkdz * std::cos(static_cast<amrex::Real>(n) * phase);
375 w_w += utils::powi(eps, n) * static_cast<amrex::Real>(n) * a[n - 1] *
376 sinh_nkdz * std::sin(static_cast<amrex::Real>(n) * phase);
377 }
378 u_w *= (c0 * std::sqrt(g / wavenumber));
379 w_w *= (c0 * std::sqrt(g / wavenumber));
380}
381
382} // namespace kynema_sgf::ocean_waves::relaxation_zones
383#endif
math_ops.H
kynema_sgf::ocean_waves::relaxation_zones
Definition relaxation_zones_ops.cpp:14
kynema_sgf::ocean_waves::relaxation_zones::stokes_wave_length
AMREX_FORCE_INLINE amrex::Real stokes_wave_length(const amrex::Real T, const amrex::Real d, const amrex::Real H, const int order, const amrex::Real g, const amrex::Real tol, const int iter_max)
Definition stokes_waves_K.H:28
kynema_sgf::ocean_waves::relaxation_zones::stokes_coefficients
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:149
kynema_sgf::ocean_waves::relaxation_zones::stokes_waves
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 hyp_lo_limit, amrex::Real hyp_hi_limit, amrex::Real &eta, amrex::Real &u_w, amrex::Real &v_w, amrex::Real &w_w)
Definition stokes_waves_K.H:278
kynema_sgf::utils::powi
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real powi(const amrex::Real &x, const int &y)
Helper function to do pow() with integer exponents and output amrex::Real.
Definition math_ops.H:18
Generated by doxygen 1.14.0