/home/runner/work/amr-wind/amr-wind/amr-wind/ocean_waves/relaxation_zones/stokes_waves_K.H Source File

AMR-Wind API: /home/runner/work/amr-wind/amr-wind/amr-wind/ocean_waves/relaxation_zones/stokes_waves_K.H Source File
AMR-Wind 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 "AMReX_REAL.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 amr_wind::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 * std::pow(S, 2.0_rt)) -
60 (400.0_rt * std::pow(S, 3.0_rt)) -
61 (71.0_rt * std::pow(S, 4.0_rt)) +
62 (146.0_rt * std::pow(S, 5.0_rt)));
63 const amrex::Real C4 =
64 C0 * numterm_C4 / (32.0_rt * std::pow(C, 5.0_rt));
65 // Calculate pure derivates
66 const amrex::Real dSdk = -2.0_rt * d * std::sinh(2.0_rt * k * d) /
67 std::pow(std::cosh(2.0_rt * k * d), 2.0_rt);
68 const amrex::Real dCdk = -dSdk;
69 const amrex::Real dC0dk =
70 d / (2.0_rt * C0 * std::pow(std::cosh(k * d), 2.0_rt));
71 // Calculate derivatives with products
72 const amrex::Real dC2dk =
73 (4.0_rt * std::pow(C, 2.0_rt) *
74 (dC0dk * (2.0_rt + 7.0_rt * std::pow(S, 2.0_rt)) +
75 C0 * 14.0_rt * S * dSdk) -
76 C0 * (2.0_rt + 7.0_rt * std::pow(S, 2.0_rt)) * 8.0_rt * C * dCdk) /
77 (16.0_rt * std::pow(C, 4.0_rt));
78 const amrex::Real dC4dk =
79 (32.0_rt * std::pow(C, 5.0_rt) *
80 (dC0dk * numterm_C4 +
81 C0 * (32.0_rt * dSdk - 232.0_rt * S * dSdk -
82 1200.0_rt * std::pow(S, 2.0_rt) * dSdk -
83 284.0_rt * std::pow(S, 3.0_rt) * dSdk +
84 730.0_rt * std::pow(S, 4.0_rt))) -
85 C0 * numterm_C4 * 160.0_rt * std::pow(C, 4.0_rt) * dCdk) /
86 (1024.0_rt * std::pow(C, 10.0_rt));
87
88 // Calculate derivative for loop convergence
89 amrex::Real dfdk =
90 (g * std::pow(C0, 2.0_rt)) + (g * k * 2.0_rt * C0 * dC0dk);
91 // Add additional terms depending on order
92 if (order >= 2) {
93 dfdk += (g * (2.0_rt * C0 * std::pow(eps, 2.0_rt) * C2 +
94 std::pow(eps, 4.0_rt) * std::pow(C2, 2.0_rt))) +
95 (g * k *
96 (2.0_rt * dC0dk * std::pow(eps, 2.0_rt) * C2 +
97 2.0_rt * C0 *
98 (2.0_rt * eps * depsdk * C2 +
99 std::pow(eps, 2.0_rt) * dC2dk) +
100 4.0_rt * std::pow(eps, 3.0_rt) * std::pow(C2, 2.0_rt) *
101 depsdk +
102 std::pow(eps, 4.0_rt) * 2.0_rt * C2 * dC2dk));
103 }
104 if (order >= 4) {
105 dfdk +=
106 (g * (2.0_rt * std::pow(eps, 4.0_rt) * C0 * C4 +
107 2.0_rt * std::pow(eps, 6.0_rt) * C2 * C4 +
108 std::pow(eps, 8.0_rt) * std::pow(C4, 2.0_rt))) +
109 (g * k *
110 (8.0_rt * std::pow(eps, 3.0_rt) * depsdk * C0 * C4 +
111 2.0_rt * std::pow(eps, 4.0_rt) * (C0 * dC4dk + C4 * dC0dk) +
112 12.0_rt * std::pow(eps, 5.0_rt) * depsdk * C2 * C4 +
113 2.0_rt * std::pow(eps, 6.0_rt) * (C2 * dC4dk + C4 * dC2dk) +
114 8.0_rt * std::pow(eps, 7.0_rt) * depsdk *
115 std::pow(C4, 2.0_rt) +
116 std::pow(eps, 8.0_rt) * 2.0_rt * C4 * dC4dk));
117 }
118 k = k - (f / dfdk);
119
120 iter += 1;
121 f = g * k * std::pow(C0, 2.0_rt);
122 // Add additional terms depending on order
123 if (order >= 2) {
124 f += g * k *
125 (2.0_rt * C0 * std::pow(eps, 2.0_rt) * C2 +
126 std::pow(eps, 4.0_rt) * std::pow(C2, 2.0_rt));
127 }
128 if (order >= 4) {
129 f += g * k *
130 (2.0_rt * std::pow(eps, 4.0_rt) * C0 * C4 +
131 2.0_rt * std::pow(eps, 6.0_rt) * C2 * C4 +
132 std::pow(eps, 8.0_rt) * std::pow(C4, 2.0_rt));
133 }
134 // Subtract omega^2 to measure convergence
135 f -= omega * omega;
136 }
137
138 if (k < tol) {
139 // Return negative wavelength if faulty result
140 return -1;
141 }
142 if (std::isnan(k)) {
143 return -2;
144 }
145 if (iter == iter_max) {
146 return -3;
147 }
148 // Return wavelength calculated from wavenumber
149 return 2.0_rt * std::numbers::pi_v<amrex::Real> / k;
150}
151
152AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void stokes_coefficients(
153 int stokes_order,
154 amrex::Real wavenumber,
155 amrex::Real water_depth,
156 amrex::Real& c0,
157 amrex::Real& a11,
158 amrex::Real& a22,
159 amrex::Real& b22,
160 amrex::Real& c2,
161 amrex::Real& a31,
162 amrex::Real& a33,
163 amrex::Real& b31,
164 amrex::Real& a42,
165 amrex::Real& a44,
166 amrex::Real& b42,
167 amrex::Real& b44,
168 amrex::Real& c4,
169 amrex::Real& a51,
170 amrex::Real& a53,
171 amrex::Real& a55,
172 amrex::Real& b53,
173 amrex::Real& b55)
174{
175
176 amrex::Real kd = wavenumber * water_depth;
177 kd = amrex::min<amrex::Real>(kd, 50.0_rt * std::numbers::pi_v<amrex::Real>);
178
179 // Exponential definition of S = sech(2kd)
180 amrex::Real S =
181 2.0_rt * std::exp(2.0_rt * kd) / (std::exp(4.0_rt * kd) + 1.0_rt);
182 amrex::Real C = 1.0_rt - S;
183 amrex::Real Sh = std::sinh(kd);
184 amrex::Real Th = std::tanh(kd);
185 // Exponential definition of coth(kd)
186 amrex::Real CTh =
187 (1.0_rt + std::exp(-2.0_rt * kd)) / (1.0_rt - std::exp(-2 * kd));
188
189 c0 = std::sqrt(Th);
190 a11 = 1.0_rt / std::sinh(kd); // Hyperbolic cosecant
191 // Second order coefficients
192 a22 = 3.0_rt * std::pow(S, 2.0_rt) / (2.0_rt * std::pow(C, 2.0_rt));
193 b22 = CTh * (1.0_rt + 2.0_rt * S) / (2.0_rt * C);
194 c2 = std::sqrt(Th) * (2.0_rt + 7.0_rt * std::pow(S, 2.0_rt)) /
195 (4.0_rt * std::pow(C, 2.0_rt));
196 if (stokes_order == 2) {
197 return;
198 }
199
200 // Third order coefficients
201 a31 = (-4.0_rt - 20.0_rt * S + 10.0_rt * std::pow(S, 2.0_rt) -
202 13.0_rt * std::pow(S, 3.0_rt)) /
203 (8.0_rt * Sh * std::pow(C, 3.0_rt));
204 a33 = (-2.0_rt * std::pow(S, 2.0_rt) + 11.0_rt * std::pow(S, 3.0_rt)) /
205 (8.0_rt * Sh * std::pow(C, 3.0_rt));
206 b31 = -3.0_rt *
207 (1.0_rt + 3.0_rt * S + 3.0_rt * std::pow(S, 2.0_rt) +
208 2.0_rt * std::pow(S, 3.0_rt)) /
209 (8.0_rt * std::pow(C, 3.0_rt));
210 if (stokes_order == 3) {
211 return;
212 }
213
214 // Fourth order coefficients
215 a42 = (12.0_rt * S - 14.0_rt * std::pow(S, 2.0_rt) -
216 264.0_rt * std::pow(S, 3.0_rt) - 45.0_rt * std::pow(S, 4.0_rt) -
217 13.0_rt * std::pow(S, 5.0_rt)) /
218 (24.0_rt * std::pow(C, 5.0_rt));
219 a44 = (10.0_rt * std::pow(S, 3.0_rt) - 174.0_rt * std::pow(S, 4.0_rt) +
220 291.0_rt * std::pow(S, 5.0_rt) + 278.0_rt * std::pow(S, 6.0_rt)) /
221 (48.0_rt * (3.0_rt + 2.0_rt * S) * std::pow(C, 5.0_rt));
222 b42 = CTh *
223 (6.0_rt - 26.0_rt * S - 182.0_rt * std::pow(S, 2.0_rt) -
224 204.0_rt * std::pow(S, 3.0_rt) - 25.0_rt * std::pow(S, 4.0_rt) +
225 26.0_rt * std::pow(S, 5.0_rt)) /
226 (6.0_rt * (3.0_rt + 2.0_rt * S) * std::pow(C, 4.0_rt));
227 b44 = CTh *
228 (24.0_rt + 92.0_rt * S + 122.0_rt * std::pow(S, 2.0_rt) +
229 66.0_rt * std::pow(S, 3.0_rt) + 67.0_rt * std::pow(S, 4.0_rt) +
230 34.0_rt * std::pow(S, 5.0_rt)) /
231 (24.0_rt * (3.0_rt + 2.0_rt * S) * std::pow(C, 4.0_rt));
232 c4 = std::sqrt(Th) *
233 (4.0_rt + 32.0_rt * S - 116.0_rt * std::pow(S, 2.0_rt) -
234 400.0_rt * std::pow(S, 3.0_rt) - 71.0_rt * std::pow(S, 4.0_rt) +
235 146.0_rt * std::pow(S, 5.0_rt)) /
236 (32.0_rt * std::pow(C, 5.0_rt));
237 if (stokes_order == 4) {
238 return;
239 }
240
241 // Fifth order coefficients
242 a51 = (-1184.0_rt + 32.0_rt * S + 13232.0_rt * std::pow(S, 2.0_rt) +
243 21712.0_rt * std::pow(S, 3.0_rt) + 20940.0_rt * std::pow(S, 4.0_rt) +
244 12554.0_rt * std::pow(S, 5.0_rt) - 500.0_rt * std::pow(S, 6.0_rt) -
245 3341.0_rt * std::pow(S, 7.0_rt) - 670.0_rt * std::pow(S, 8.0_rt)) /
246 (64.0_rt * Sh * (3.0_rt + 2.0_rt * S) * (4.0_rt + S) *
247 std::pow(C, 6.0_rt));
248 a53 = (4.0_rt * S + 105.0_rt * std::pow(S, 2.0_rt) +
249 198.0_rt * std::pow(S, 3.0_rt) - 1376.0_rt * std::pow(S, 4.0_rt) -
250 1302.0_rt * std::pow(S, 5.0_rt) - 117.0_rt * std::pow(S, 6.0_rt) +
251 58.0_rt * std::pow(S, 7.0_rt)) /
252 (32.0_rt * Sh * (3.0_rt + 2.0_rt * S) * std::pow(C, 6.0_rt));
253 a55 = (-6.0_rt * std::pow(S, 3.0_rt) + 272.0_rt * std::pow(S, 4.0_rt) -
254 1552.0_rt * std::pow(S, 5.0_rt) + 852.0_rt * std::pow(S, 6.0_rt) +
255 2029.0_rt * std::pow(S, 7.0_rt) + 430.0_rt * std::pow(S, 8.0_rt)) /
256 (64.0_rt * Sh * (3.0_rt + 2.0_rt * S) * (4.0_rt + S) *
257 std::pow(C, 6.0_rt));
258 b53 =
259 9.0_rt *
260 (132.0_rt + 17.0_rt * S - 2216.0_rt * std::pow(S, 2.0_rt) -
261 5897.0_rt * std::pow(S, 3.0_rt) - 6292.0_rt * std::pow(S, 4.0_rt) -
262 2687.0_rt * std::pow(S, 5.0_rt) + 194.0_rt * std::pow(S, 6.0_rt) +
263 467.0_rt * std::pow(S, 7.0_rt) + 82.0_rt * std::pow(S, 8.0_rt)) /
264 (128.0_rt * (3.0_rt + 2.0_rt * S) * (4.0_rt + S) * std::pow(C, 6.0_rt));
265 b55 =
266 5.0_rt *
267 (300.0_rt + 1579.0_rt * S + 3176.0_rt * std::pow(S, 2.0_rt) +
268 2949.0_rt * std::pow(S, 3.0_rt) + 1188.0_rt * std::pow(S, 4.0_rt) +
269 675.0_rt * std::pow(S, 5.0_rt) + 1326.0_rt * std::pow(S, 6.0_rt) +
270 827.0_rt * std::pow(S, 7.0_rt) + 130.0_rt * std::pow(S, 8.0_rt)) /
271 (384.0_rt * (3.0_rt + 2.0_rt * S) * (4.0_rt + S) * std::pow(C, 6.0_rt));
272 if (stokes_order == 5) {
273 return;
274 }
275
276 if (stokes_order > 5 || stokes_order < 2) {
277 amrex::Abort(
278 "invalid stokes order specified. It should be between 2,3,4 or 5 ");
279 }
280}
281
282// Based on Fenton 1985
283AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void stokes_waves(
284 int stokes_order,
285 amrex::Real wavelength,
286 amrex::Real water_depth,
287 amrex::Real wave_height,
288 amrex::Real zsl,
289 amrex::Real g,
290 amrex::Real x,
291 amrex::Real z,
292 amrex::Real time,
293 amrex::Real phase_offset,
294 amrex::Real& eta,
295 amrex::Real& u_w,
296 amrex::Real& v_w,
297 amrex::Real& w_w)
298{
299 const amrex::Real wavenumber =
300 2.0_rt * std::numbers::pi_v<amrex::Real> / wavelength;
301
302 // some parameters
303 amrex::Real c0{0.0_rt};
304 amrex::Real a11{0.0_rt}, a22{0.0_rt}, b22{0.0_rt}, c2{0.0_rt};
305 amrex::Real a31{0.0_rt}, a33{0.0_rt}, b31{0.0_rt}, a42{0.0_rt}, a44{0.0_rt};
306 amrex::Real b42{0.0_rt}, b44{0.0_rt}, c4{0.0_rt};
307 amrex::Real a51{0.0_rt}, a53{0.0_rt}, a55{0.0_rt}, b53{0.0_rt}, b55{0.0_rt};
308
309 stokes_coefficients(
310 stokes_order, wavenumber, water_depth, c0, a11, a22, b22, c2, a31, a33,
311 b31, a42, a44, b42, b44, c4, a51, a53, a55, b53, b55);
312
313 const amrex::Real eps = wavenumber * wave_height / 2.0_rt; // Steepness (ka)
314 const amrex::Real c =
315 (c0 + std::pow(eps, 2.0_rt) * c2 + std::pow(eps, 4.0_rt) * c4) *
316 std::sqrt(g / wavenumber);
317
318 const amrex::Real omega = c * wavenumber;
319 const amrex::Real phase = (wavenumber * x) - (omega * time) - phase_offset;
320
321 eta = ((eps * std::cos(phase) // first order term
322 + std::pow(eps, 2.0_rt) * b22 *
323 std::cos(2.0_rt * phase) // second order term
324 + std::pow(eps, 3.0_rt) * b31 *
325 (std::cos(phase) - std::cos(3.0_rt * phase)) +
326 std::pow(eps, 4.0_rt) * (b42 * std::cos(2.0_rt * phase) +
327 b44 * std::cos(4.0_rt * phase)) +
328 std::pow(eps, 5.0_rt) * (-(b53 + b55) * std::cos(phase) +
329 b53 * std::cos(3.0_rt * phase) +
330 b55 * std::cos(5.0_rt * phase))) /
331 wavenumber) +
332 zsl;
333
334 // Compute velocities components using Eq.(21) Eq.(23) from Kinnas
335 // https://www.caee.utexas.edu/prof/kinnas/ce358/oenotes/kinnas_stokes11.pdf
336 // Define coefficients using Eq.(19)
337
338 const int MAX_ORDER = 5;
339 amrex::GpuArray<amrex::Real, MAX_ORDER> a;
340 if (stokes_order == 2) {
341 a[0] = a11;
342 a[1] = a22;
343 }
344 if (stokes_order == 3) {
345 a[0] = a11 + ((eps * eps) * a31);
346 a[1] = a22;
347 a[2] = a33;
348 }
349 if (stokes_order == 4) {
350 a[0] = a11 + ((eps * eps) * a31);
351 a[1] = a22 + ((eps * eps) * a42);
352 a[2] = a33;
353 a[3] = a44;
354 }
355 if (stokes_order == 5) {
356 a[0] = a11 + ((eps * eps) * a31) + (std::pow(eps, 4.0_rt) * a51);
357 a[1] = a22 + ((eps * eps) * a42);
358 a[2] = a33 + ((eps * eps) * a53);
359 a[3] = a44;
360 a[4] = a55;
361 }
362
363 u_w = 0.0_rt;
364 v_w = 0.0_rt;
365 w_w = 0.0_rt;
366 for (int n = 1; n <= stokes_order; ++n) {
367 // Upper bound for deep water case
368 // This ensure finite values of velocity for large kd's
369 amrex::Real nkdz = n * wavenumber * (water_depth + (z - zsl));
370 nkdz = amrex::min<amrex::Real>(
371 nkdz, 50.0_rt * std::numbers::pi_v<amrex::Real>);
372 u_w += std::pow(eps, static_cast<amrex::Real>(n)) *
373 static_cast<amrex::Real>(n) * a[n - 1] * std::cosh(nkdz) *
374 std::cos(static_cast<amrex::Real>(n) * phase);
375 w_w += std::pow(eps, static_cast<amrex::Real>(n)) *
376 static_cast<amrex::Real>(n) * a[n - 1] * std::sinh(nkdz) *
377 std::sin(static_cast<amrex::Real>(n) * phase);
378 }
379 u_w *= (c0 * std::sqrt(g / wavenumber));
380 w_w *= (c0 * std::sqrt(g / wavenumber));
381}
382
383} // namespace amr_wind::ocean_waves::relaxation_zones
384#endif
amr_wind::ocean_waves::relaxation_zones
Definition relaxation_zones_ops.cpp:14
amr_wind::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:152
amr_wind::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 &eta, amrex::Real &u_w, amrex::Real &v_w, amrex::Real &w_w)
Definition stokes_waves_K.H:283
amr_wind::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
Generated by doxygen 1.14.0