/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_H_
2#define STOKES_WAVES_H_
3
4#include <AMReX_FArrayBox.H>
5#include <cmath>
6#include "AMReX_REAL.H"
7
8using namespace amrex::literals;
9
10// Stokes waves theory adapted from
11// Fenton, J., Fifth Order Stokes Theory for Steady Waves
12// Journal of Waterway, Port, Coastal and Ocean Engineering, 1985, 111, 216-234
13
14// Updated Table 1 from Fenton 1985 found in
15// Fenton, J., Nonlinear Wave Theories
16// The Sea Vol. 9 Ocean Engineering Science, 1990
17// https://johndfenton.com/Papers/Fenton90b-Nonlinear-wave-theories.pdf
18
19// Relevant results are summarized in
20// Kinnas S. A., Notes on fifth order gravity wave theory
21// https://www.caee.utexas.edu/prof/kinnas/ce358/oenotes/kinnas_stokes11.pdf
22
23namespace amr_wind::ocean_waves::relaxation_zones {
24
25// Compute wavelength as a function of wave period, water depth, and g
26AMREX_FORCE_INLINE amrex::Real stokes_wave_length(
27 const amrex::Real T,
28 const amrex::Real d,
29 const amrex::Real H,
30 const int order,
31 const amrex::Real g,
32 const amrex::Real tol,
33 const int iter_max)
34{
35 // Calculate constants and derivatives that do not change with iteration
36 const amrex::Real omega = 2.0_rt * static_cast<amrex::Real>(M_PI) / T;
37 const amrex::Real depsdk = H / 2.0_rt;
38
39 // Begin Newton-Raphson iterations
40 int iter = 0;
41 // First guess is first-order
42 amrex::Real k = (omega * omega) / g;
43 // Cannot skip loop
44 amrex::Real f = tol + 1.0_rt;
45 while (std::abs(f) > tol && iter < iter_max) {
46 // Calculate current constants
47
48 // Exponential definition of S = sech(2kd)
49 const amrex::Real S = 2.0_rt * std::exp(2.0_rt * k * d) /
50 (std::exp(4.0_rt * k * d) + 1.0_rt);
51 const amrex::Real C = 1.0_rt - S;
52 const amrex::Real eps = k * H / 2.0_rt;
53 const amrex::Real C0 = std::sqrt(std::tanh(k * d));
54 const amrex::Real C2 =
55 C0 * (2.0_rt + 7.0_rt * S * S) / (4.0_rt * C * C);
56 const amrex::Real numterm_C4 =
57 (4.0_rt + 32.0_rt * S - 116.0_rt * std::pow(S, 2.0_rt) -
58 400.0_rt * std::pow(S, 3.0_rt) - 71.0_rt * std::pow(S, 4.0_rt) +
59 146.0_rt * std::pow(S, 5.0_rt));
60 const amrex::Real C4 =
61 C0 * numterm_C4 / (32.0_rt * std::pow(C, 5.0_rt));
62 // Calculate pure derivates
63 const amrex::Real dSdk = -2.0_rt * d * std::sinh(2.0_rt * k * d) /
64 std::pow(std::cosh(2.0_rt * k * d), 2.0_rt);
65 const amrex::Real dCdk = -dSdk;
66 const amrex::Real dC0dk =
67 d / (2.0_rt * C0 * std::pow(std::cosh(k * d), 2.0_rt));
68 // Calculate derivatives with products
69 const amrex::Real dC2dk =
70 (4.0_rt * std::pow(C, 2.0_rt) *
71 (dC0dk * (2.0_rt + 7.0_rt * std::pow(S, 2.0_rt)) +
72 C0 * 14.0_rt * S * dSdk) -
73 C0 * (2.0_rt + 7.0_rt * std::pow(S, 2.0_rt)) * 8.0_rt * C * dCdk) /
74 (16.0_rt * std::pow(C, 4.0_rt));
75 const amrex::Real dC4dk =
76 (32.0_rt * std::pow(C, 5.0_rt) *
77 (dC0dk * numterm_C4 +
78 C0 * (32.0_rt * dSdk - 232.0_rt * S * dSdk -
79 1200.0_rt * std::pow(S, 2.0_rt) * dSdk -
80 284.0_rt * std::pow(S, 3.0_rt) * dSdk +
81 730.0_rt * std::pow(S, 4.0_rt))) -
82 C0 * numterm_C4 * 160.0_rt * std::pow(C, 4.0_rt) * dCdk) /
83 (1024.0_rt * std::pow(C, 10.0_rt));
84
85 // Calculate derivative for loop convergence
86 amrex::Real dfdk =
87 g * std::pow(C0, 2.0_rt) + g * k * 2.0_rt * C0 * dC0dk;
88 // Add additional terms depending on order
89 if (order >= 2) {
90 dfdk += g * (2.0_rt * C0 * std::pow(eps, 2.0_rt) * C2 +
91 std::pow(eps, 4.0_rt) * std::pow(C2, 2.0_rt)) +
92 g * k *
93 (2.0_rt * dC0dk * std::pow(eps, 2.0_rt) * C2 +
94 2.0_rt * C0 *
95 (2.0_rt * eps * depsdk * C2 +
96 std::pow(eps, 2.0_rt) * dC2dk) +
97 4.0_rt * std::pow(eps, 3.0_rt) * std::pow(C2, 2.0_rt) *
98 depsdk +
99 std::pow(eps, 4.0_rt) * 2.0_rt * C2 * dC2dk);
100 }
101 if (order >= 4) {
102 dfdk += g * (2.0_rt * std::pow(eps, 4.0_rt) * C0 * C4 +
103 2.0_rt * std::pow(eps, 6.0_rt) * C2 * C4 +
104 std::pow(eps, 8.0_rt) * std::pow(C4, 2.0_rt)) +
105 g * k *
106 (8.0_rt * std::pow(eps, 3.0_rt) * depsdk * C0 * C4 +
107 2.0_rt * std::pow(eps, 4.0_rt) *
108 (C0 * dC4dk + C4 * dC0dk) +
109 12.0_rt * std::pow(eps, 5.0_rt) * depsdk * C2 * C4 +
110 2.0_rt * std::pow(eps, 6.0_rt) *
111 (C2 * dC4dk + C4 * dC2dk) +
112 8.0_rt * std::pow(eps, 7.0_rt) * depsdk *
113 std::pow(C4, 2.0_rt) +
114 std::pow(eps, 8.0_rt) * 2.0_rt * C4 * dC4dk);
115 }
116 k = k - f / dfdk;
117
118 iter += 1;
119 f = g * k * std::pow(C0, 2.0_rt);
120 // Add additional terms depending on order
121 if (order >= 2) {
122 f += g * k *
123 (2.0_rt * C0 * std::pow(eps, 2.0_rt) * C2 +
124 std::pow(eps, 4.0_rt) * std::pow(C2, 2.0_rt));
125 }
126 if (order >= 4) {
127 f += g * k *
128 (2.0_rt * std::pow(eps, 4.0_rt) * C0 * C4 +
129 2.0_rt * std::pow(eps, 6.0_rt) * C2 * C4 +
130 std::pow(eps, 8.0_rt) * std::pow(C4, 2.0_rt));
131 }
132 // Subtract omega^2 to measure convergence
133 f -= omega * omega;
134 }
135
136 if (k < tol) {
137 // Return negative wavelength if faulty result
138 return -1;
139 }
140 if (std::isnan(k)) {
141 return -2;
142 }
143 if (iter == iter_max) {
144 return -3;
145 }
146 // Return wavelength calculated from wavenumber
147 return 2.0_rt * static_cast<amrex::Real>(M_PI) / k;
148}
149
150AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void stokes_coefficients(
151 int stokes_order,
152 amrex::Real wavenumber,
153 amrex::Real water_depth,
154 amrex::Real& c0,
155 amrex::Real& a11,
156 amrex::Real& a22,
157 amrex::Real& b22,
158 amrex::Real& c2,
159 amrex::Real& a31,
160 amrex::Real& a33,
161 amrex::Real& b31,
162 amrex::Real& a42,
163 amrex::Real& a44,
164 amrex::Real& b42,
165 amrex::Real& b44,
166 amrex::Real& c4,
167 amrex::Real& a51,
168 amrex::Real& a53,
169 amrex::Real& a55,
170 amrex::Real& b53,
171 amrex::Real& b55)
172{
173
174 amrex::Real kd = wavenumber * water_depth;
175 if (kd > 50.0_rt * static_cast<amrex::Real>(M_PI)) {
176 kd = 50.0_rt * static_cast<amrex::Real>(M_PI);
177 }
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 * static_cast<amrex::Real>(M_PI) / 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 if (nkdz > 50.0_rt * static_cast<amrex::Real>(M_PI)) {
371 nkdz = 50.0_rt * static_cast<amrex::Real>(M_PI);
372 }
373 u_w += std::pow(eps, static_cast<amrex::Real>(n)) *
374 static_cast<amrex::Real>(n) * a[n - 1] * std::cosh(nkdz) *
375 std::cos(static_cast<amrex::Real>(n) * phase);
376 w_w += std::pow(eps, static_cast<amrex::Real>(n)) *
377 static_cast<amrex::Real>(n) * a[n - 1] * std::sinh(nkdz) *
378 std::sin(static_cast<amrex::Real>(n) * phase);
379 }
380 u_w *= (c0 * std::sqrt(g / wavenumber));
381 w_w *= (c0 * std::sqrt(g / wavenumber));
382}
383
384} // namespace amr_wind::ocean_waves::relaxation_zones
385#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:150
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:26
Generated by doxygen 1.14.0