diff --git a/amr-wind/ocean_waves/relaxation_zones/stokes_waves_K.H b/amr-wind/ocean_waves/relaxation_zones/stokes_waves_K.H index 20cf57a4c6..2d729de2d5 100644 --- a/amr-wind/ocean_waves/relaxation_zones/stokes_waves_K.H +++ b/amr-wind/ocean_waves/relaxation_zones/stokes_waves_K.H @@ -8,6 +8,15 @@ // Fenton, J., Fifth Order Stokes Theory for Steady Waves // Journal of Waterway, Port, Coastal and Ocean Engineering, 1985, 111, 216-234 +// Updated Table 1 from Fenton 1985 found in +// Fenton, J., Nonlinear Wave Theories +// The Sea Vol. 9 Ocean Engineering Science, 1990 +// https://johndfenton.com/Papers/Fenton90b-Nonlinear-wave-theories.pdf + +// Relevant results are summarized in +// Kinnas S. A., Notes on fifth order gravity wave theory +// https://www.caee.utexas.edu/prof/kinnas/ce358/oenotes/kinnas_stokes11.pdf + namespace amr_wind::ocean_waves::relaxation_zones { // Compute wavelength as a function of wave period, water depth, and g @@ -32,7 +41,10 @@ AMREX_FORCE_INLINE amrex::Real stokes_wave_length( amrex::Real f = tol + 1.0; while (std::abs(f) > tol && iter < iter_max) { // Calculate current constants - const amrex::Real S = 1.0 / std::cosh(2.0 * k * d); + + // Exponential definition of S = sech(2kd) + const amrex::Real S = + 2. * std::exp(2. * k * d) / (std::exp(4. * k * d) + 1.); const amrex::Real C = 1.0 - S; const amrex::Real eps = k * H / 2.0; const amrex::Real C0 = std::sqrt(std::tanh(k * d)); @@ -123,16 +135,14 @@ AMREX_FORCE_INLINE amrex::Real stokes_wave_length( } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void stokes_coefficients( - int StokesOrder, + int stokes_order, amrex::Real wavenumber, - amrex::Real waterdepth, + amrex::Real water_depth, amrex::Real& c0, amrex::Real& a11, amrex::Real& a22, amrex::Real& b22, amrex::Real& c2, - amrex::Real& d2, - amrex::Real& e2, amrex::Real& a31, amrex::Real& a33, amrex::Real& b31, @@ -141,8 +151,6 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void stokes_coefficients( amrex::Real& b42, amrex::Real& b44, amrex::Real& c4, - amrex::Real& d4, - amrex::Real& e4, amrex::Real& a51, amrex::Real& a53, amrex::Real& a55, @@ -150,65 +158,60 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void stokes_coefficients( amrex::Real& b55) { - amrex::Real kd = wavenumber * waterdepth; + amrex::Real kd = wavenumber * water_depth; if (kd > 50. * M_PI) { kd = 50 * M_PI; } + // Exponential definition of S = sech(2kd) amrex::Real S = 2. * std::exp(2. * kd) / (std::exp(4. * kd) + 1.); + amrex::Real C = 1.0 - S; amrex::Real Sh = std::sinh(kd); amrex::Real Th = std::tanh(kd); + // Exponential definition of coth(kd) amrex::Real CTh = (1. + std::exp(-2 * kd)) / (1. - std::exp(-2 * kd)); c0 = std::sqrt(Th); a11 = 1. / std::sinh(kd); // Hyperbolic cosecant // Second order coefficients - a22 = 3. * std::pow(S, 2) / (2 * std::pow(1 - S, 2)); - b22 = CTh * (1 + 2. * S) / (2 * (1 - S)); - c2 = std::sqrt(Th) * (2 + 7 * std::pow(S, 2)) / (4 * std::pow(1 - S, 2)); - d2 = -std::sqrt(CTh) / 2.; - e2 = Th * (2 + 2 * S + 5 * std::pow(S, 2)) / (4 * std::pow(1 - S, 2)); - if (StokesOrder == 2) { + a22 = 3. * std::pow(S, 2) / (2 * std::pow(C, 2)); + b22 = CTh * (1 + 2. * S) / (2 * C); + c2 = std::sqrt(Th) * (2 + 7 * std::pow(S, 2)) / (4 * std::pow(C, 2)); + if (stokes_order == 2) { return; } // Third order coefficients a31 = (-4 - 20 * S + 10 * std::pow(S, 2) - 13 * std::pow(S, 3)) / - (8 * Sh * std::pow(1 - S, 3)); - a33 = (-2 * std::pow(S, 2) + 11 * std::pow(S, 3)) / - (8 * Sh * std::pow(1 - S, 3)); + (8 * Sh * std::pow(C, 3)); + a33 = + (-2 * std::pow(S, 2) + 11 * std::pow(S, 3)) / (8 * Sh * std::pow(C, 3)); b31 = -3 * (1 + 3 * S + 3 * std::pow(S, 2) + 2 * std::pow(S, 3)) / - (8 * std::pow(1 - S, 3)); - if (StokesOrder == 3) { + (8 * std::pow(C, 3)); + if (stokes_order == 3) { return; } // Fourth order coefficients a42 = (12 * S - 14 * std::pow(S, 2) - 264 * std::pow(S, 3) - 45 * std::pow(S, 4) - 13 * std::pow(S, 5)) / - (24 * std::pow(1 - S, 5)); + (24 * std::pow(C, 5)); a44 = (10 * std::pow(S, 3) - 174 * std::pow(S, 4) + 291 * std::pow(S, 5) + 278 * std::pow(S, 6)) / - (48 * (3 + 2 * S) * std::pow(1 - S, 5)); + (48 * (3 + 2 * S) * std::pow(C, 5)); b42 = CTh * (6 - 26 * S - 182 * std::pow(S, 2) - 204 * std::pow(S, 3) - 25 * std::pow(S, 4) + 26 * std::pow(S, 5)) / - (6 * (3 + 2 * S) * std::pow(1 - S, 4)); + (6 * (3 + 2 * S) * std::pow(C, 4)); b44 = CTh * (24 + 92 * S + 122 * std::pow(S, 2) + 66 * std::pow(S, 3) + 67 * std::pow(S, 4) + 34 * std::pow(S, 5)) / - (24 * (3 + 2 * S) * std::pow(1 - S, 4)); + (24 * (3 + 2 * S) * std::pow(C, 4)); c4 = std::sqrt(Th) * (4 + 32 * S - 116 * std::pow(S, 2) - 400 * std::pow(S, 3) - 71 * std::pow(S, 4) + 146 * std::pow(S, 5)) / - (32 * std::pow(1 - S, 5)); - d4 = std::sqrt(CTh) * (2 + 4 * S + std::pow(S, 2) + 2 * std::pow(S, 3)) / - (8 * std::pow(1 - S, 3)); - e4 = Th * - (8 + 12 * S - 152 * std::pow(S, 2) - 308 * std::pow(S, 3) - - 42 * std::pow(S, 4) + 77 * std::pow(S, 5)) / - (32 * std::pow(1 - S, 5)); - if (StokesOrder == 4) { + (32 * std::pow(C, 5)); + if (stokes_order == 4) { return; } @@ -217,148 +220,42 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void stokes_coefficients( (-1184 + 32 * S + 13232 * std::pow(S, 2) + 21712 * std::pow(S, 3) + 20940 * std::pow(S, 4) + 12554 * std::pow(S, 5) - 500 * std::pow(S, 6) - 3341 * std::pow(S, 7) - 670 * std::pow(S, 8)) / - (64 * Sh * (3 + 2 * S) * (4 + S) * std::pow(1 - S, 6)); + (64 * Sh * (3 + 2 * S) * (4 + S) * std::pow(C, 6)); a53 = (4 * S + 105 * pow(S, 2) + 198 * std::pow(S, 3) - 1376 * std::pow(S, 4) - 1302 * std::pow(S, 5) - 117 * std::pow(S, 6) + 58 * std::pow(S, 7)) / - (32 * Sh * (3 + 2 * S) * std::pow(1 - S, 6)); + (32 * Sh * (3 + 2 * S) * std::pow(C, 6)); a55 = (-6 * std::pow(S, 3) + 272 * std::pow(S, 4) - 1552 * std::pow(S, 5) + 852 * std::pow(S, 6) + 2029 * std::pow(S, 7) + 430 * std::pow(S, 8)) / - (64 * Sh * (3 + 2 * S) * (4 + S) * std::pow(1 - S, 6)); + (64 * Sh * (3 + 2 * S) * (4 + S) * std::pow(C, 6)); b53 = 9 * (132 + 17 * S - 2216 * std::pow(S, 2) - 5897 * std::pow(S, 3) - 6292 * std::pow(S, 4) - 2687 * std::pow(S, 5) + 194 * std::pow(S, 6) + 467 * std::pow(S, 7) + 82 * std::pow(S, 8)) / - (128 * (3 + 2 * S) * (4 + S) * std::pow(1 - S, 6)); + (128 * (3 + 2 * S) * (4 + S) * std::pow(C, 6)); b55 = 5 * (300 + 1579 * S + 3176 * std::pow(S, 2) + 2949 * std::pow(S, 3) + 1188 * std::pow(S, 4) + 675 * std::pow(S, 5) + 1326 * std::pow(S, 6) + 827 * std::pow(S, 7) + 130 * std::pow(S, 8)) / - (384 * (3 + 2 * S) * (4 + S) * std::pow(1 - S, 6)); - if (StokesOrder == 5) { + (384 * (3 + 2 * S) * (4 + S) * std::pow(C, 6)); + if (stokes_order == 5) { return; } - if (StokesOrder > 5 || StokesOrder < 2) { + if (stokes_order > 5 || stokes_order < 2) { amrex::Abort( "invalid stokes order specified. It should be between 2,3,4 or 5 "); } } -AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real stokes_sinh_sin( - int m, - int n, - amrex::Real phase, - amrex::Real a11, - amrex::Real a22, - amrex::Real a31, - amrex::Real a33, - amrex::Real a42, - amrex::Real a44, - amrex::Real a51, - amrex::Real a53, - amrex::Real a55, - amrex::Real eps, - amrex::Real wavenumber, - amrex::Real waterdepth, - amrex::Real zsl, - amrex::Real z) -{ - amrex::Real D = 0.0; - if (m == 1 && n == 1) { - D = a11; - } - if (m == 2 && n == 2) { - D = a22; - } - if (m == 3 && n == 1) { - D = a31; - } - if (m == 3 && n == 3) { - D = a33; - } - if (m == 4 && n == 2) { - D = a42; - } - if (m == 4 && n == 4) { - D = a44; - } - if (m == 5 && n == 1) { - D = a51; - } - if (m == 5 && n == 3) { - D = a53; - } - if (m == 5 && n == 5) { - D = a55; - } - - return std::pow(eps, m) * D * n * wavenumber * - std::sinh(n * wavenumber * (waterdepth + (z - zsl))) * - std::sin(n * phase); -} - -AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real stokes_cosh_cos( - int m, - int n, - amrex::Real phase, - amrex::Real a11, - amrex::Real a22, - amrex::Real a31, - amrex::Real a33, - amrex::Real a42, - amrex::Real a44, - amrex::Real a51, - amrex::Real a53, - amrex::Real a55, - amrex::Real eps, - amrex::Real wavenumber, - amrex::Real waterdepth, - amrex::Real zsl, - amrex::Real z) -{ - amrex::Real D = 0.0; - if (m == 1 && n == 1) { - D = a11; - } - if (m == 2 && n == 2) { - D = a22; - } - if (m == 3 && n == 1) { - D = a31; - } - if (m == 3 && n == 3) { - D = a33; - } - if (m == 4 && n == 2) { - D = a42; - } - if (m == 4 && n == 4) { - D = a44; - } - if (m == 5 && n == 1) { - D = a51; - } - if (m == 5 && n == 3) { - D = a53; - } - if (m == 5 && n == 5) { - D = a55; - } - - return std::pow(eps, m) * D * n * wavenumber * - std::cosh(n * wavenumber * (waterdepth + (z - zsl))) * - std::cos(n * phase); -} - // Based on Fenton 1985 -AMREX_GPU_DEVICE AMREX_FORCE_INLINE void stokes_waves( - int StokesOrder, +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void stokes_waves( + int stokes_order, amrex::Real wavelength, - amrex::Real waterdepth, - amrex::Real waveheight, + amrex::Real water_depth, + amrex::Real wave_height, amrex::Real zsl, amrex::Real g, amrex::Real x, @@ -374,22 +271,18 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE void stokes_waves( // some parameters amrex::Real c0{0.0}; - amrex::Real a11{0.0}, a22{0.0}, b22{0.0}, c2{0.0}, d2{0.0}, e2{0.0}; + amrex::Real a11{0.0}, a22{0.0}, b22{0.0}, c2{0.0}; amrex::Real a31{0.0}, a33{0.0}, b31{0.0}, a42{0.0}, a44{0.0}; - amrex::Real b42{0.0}, b44{0.0}, c4{0.0}, d4{0.0}, e4{0.0}; + amrex::Real b42{0.0}, b44{0.0}, c4{0.0}; amrex::Real a51{0.0}, a53{0.0}, a55{0.0}, b53{0.0}, b55{0.0}; stokes_coefficients( - StokesOrder, wavenumber, waterdepth, c0, a11, a22, b22, c2, d2, e2, a31, - a33, b31, a42, a44, b42, b44, c4, d4, e4, a51, a53, a55, b53, b55); + stokes_order, wavenumber, water_depth, c0, a11, a22, b22, c2, a31, a33, + b31, a42, a44, b42, b44, c4, a51, a53, a55, b53, b55); - const amrex::Real eps = wavenumber * waveheight / 2.; // Steepness (ka) + const amrex::Real eps = wavenumber * wave_height / 2.; // Steepness (ka) const amrex::Real c = (c0 + std::pow(eps, 2) * c2 + std::pow(eps, 4) * c4) * std::sqrt(g / wavenumber); - // const amrex::Real Q = - // c * waterdepth * std::sqrt(std::pow(k, 3) / g) + - // d2 * std::pow(eps, 2) + - // d4 * std::pow(eps, 4) * std::sqrt(g / std::pow(k, 3)); const amrex::Real omega = c * wavenumber; const amrex::Real phase = wavenumber * x - omega * time - phase_offset; @@ -405,69 +298,52 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE void stokes_waves( wavenumber + zsl; - // Compute cosines - const amrex::Real cc11 = stokes_cosh_cos( - 1, 1, phase, a11, a22, a31, a33, a42, a44, a51, a53, a55, eps, - wavenumber, waterdepth, zsl, z); - const amrex::Real cc22 = stokes_cosh_cos( - 2, 2, phase, a11, a22, a31, a33, a42, a44, a51, a53, a55, eps, - wavenumber, waterdepth, zsl, z); - const amrex::Real cc31 = stokes_cosh_cos( - 3, 1, phase, a11, a22, a31, a33, a42, a44, a51, a53, a55, eps, - wavenumber, waterdepth, zsl, z); - const amrex::Real cc33 = stokes_cosh_cos( - 3, 3, phase, a11, a22, a31, a33, a42, a44, a51, a53, a55, eps, - wavenumber, waterdepth, zsl, z); - const amrex::Real cc42 = stokes_cosh_cos( - 4, 2, phase, a11, a22, a31, a33, a42, a44, a51, a53, a55, eps, - wavenumber, waterdepth, zsl, z); - const amrex::Real cc44 = stokes_cosh_cos( - 4, 4, phase, a11, a22, a31, a33, a42, a44, a51, a53, a55, eps, - wavenumber, waterdepth, zsl, z); - const amrex::Real cc51 = stokes_cosh_cos( - 5, 1, phase, a11, a22, a31, a33, a42, a44, a51, a53, a55, eps, - wavenumber, waterdepth, zsl, z); - const amrex::Real cc53 = stokes_cosh_cos( - 5, 3, phase, a11, a22, a31, a33, a42, a44, a51, a53, a55, eps, - wavenumber, waterdepth, zsl, z); - const amrex::Real cc55 = stokes_cosh_cos( - 5, 5, phase, a11, a22, a31, a33, a42, a44, a51, a53, a55, eps, - wavenumber, waterdepth, zsl, z); + // Compute velocities components using Eq.(21) Eq.(23) from Kinnas + // https://www.caee.utexas.edu/prof/kinnas/ce358/oenotes/kinnas_stokes11.pdf + // Define coefficients using Eq.(19) - // Compute sines - const amrex::Real ss11 = stokes_sinh_sin( - 1, 1, phase, a11, a22, a31, a33, a42, a44, a51, a53, a55, eps, - wavenumber, waterdepth, zsl, z); - const amrex::Real ss22 = stokes_sinh_sin( - 2, 2, phase, a11, a22, a31, a33, a42, a44, a51, a53, a55, eps, - wavenumber, waterdepth, zsl, z); - const amrex::Real ss31 = stokes_sinh_sin( - 3, 1, phase, a11, a22, a31, a33, a42, a44, a51, a53, a55, eps, - wavenumber, waterdepth, zsl, z); - const amrex::Real ss33 = stokes_sinh_sin( - 3, 3, phase, a11, a22, a31, a33, a42, a44, a51, a53, a55, eps, - wavenumber, waterdepth, zsl, z); - const amrex::Real ss42 = stokes_sinh_sin( - 4, 2, phase, a11, a22, a31, a33, a42, a44, a51, a53, a55, eps, - wavenumber, waterdepth, zsl, z); - const amrex::Real ss44 = stokes_sinh_sin( - 4, 4, phase, a11, a22, a31, a33, a42, a44, a51, a53, a55, eps, - wavenumber, waterdepth, zsl, z); - const amrex::Real ss51 = stokes_sinh_sin( - 5, 1, phase, a11, a22, a31, a33, a42, a44, a51, a53, a55, eps, - wavenumber, waterdepth, zsl, z); - const amrex::Real ss53 = stokes_sinh_sin( - 5, 3, phase, a11, a22, a31, a33, a42, a44, a51, a53, a55, eps, - wavenumber, waterdepth, zsl, z); - const amrex::Real ss55 = stokes_sinh_sin( - 5, 5, phase, a11, a22, a31, a33, a42, a44, a51, a53, a55, eps, - wavenumber, waterdepth, zsl, z); + const int MAX_ORDER = 5; + amrex::GpuArray a; + if (stokes_order == 2) { + a[0] = a11; + a[1] = a22; + } + if (stokes_order == 3) { + a[0] = a11 + (eps * eps) * a31; + a[1] = a22; + a[2] = a33; + } + if (stokes_order == 4) { + a[0] = a11 + (eps * eps) * a31; + a[1] = a22 + (eps * eps) * a42; + a[2] = a33; + a[3] = a44; + } + if (stokes_order == 5) { + a[0] = a11 + (eps * eps) * a31 + std::pow(eps, 4) * a51; + a[1] = a22 + (eps * eps) * a42; + a[2] = a33 + (eps * eps) * a53; + a[3] = a44; + a[4] = a55; + } - u_w = c0 * std::sqrt(g / std::pow(wavenumber, 3)) * - (cc11 + cc22 + cc31 + cc33 + cc42 + cc44 + cc51 + cc53 + cc55); - v_w = 0.0; - w_w = c0 * std::sqrt(g / std::pow(wavenumber, 3)) * - (ss11 + ss22 + ss31 + ss33 + ss42 + ss44 + ss51 + ss53 + ss55); + u_w = 0.; + v_w = 0.; + w_w = 0.; + for (int n = 1; n <= stokes_order; ++n) { + // Upper bound for deep water case + // This ensure finite values of velocity for large kd's + amrex::Real nkdz = n * wavenumber * (water_depth + (z - zsl)); + if (nkdz > 50. * M_PI) { + nkdz = 50. * M_PI; + } + u_w += std::pow(eps, n) * n * a[n - 1] * std::cosh(nkdz) * + std::cos(n * phase); + w_w += std::pow(eps, n) * n * a[n - 1] * std::sinh(nkdz) * + std::sin(n * phase); + } + u_w *= (c0 * std::sqrt(g / wavenumber)); + w_w *= (c0 * std::sqrt(g / wavenumber)); } } // namespace amr_wind::ocean_waves::relaxation_zones diff --git a/unit_tests/ocean_waves/test_wave_theories.cpp b/unit_tests/ocean_waves/test_wave_theories.cpp index 50d1b0274f..65034f077e 100644 --- a/unit_tests/ocean_waves/test_wave_theories.cpp +++ b/unit_tests/ocean_waves/test_wave_theories.cpp @@ -10,17 +10,17 @@ class WaveTheoriesTest : public MeshTest TEST_F(WaveTheoriesTest, StokesWaves) { - amrex::Real CoeffTol = 1e-4; + constexpr amrex::Real coeff_tol = 1e-4; - amrex::Real wavenumber = 2.0; - amrex::Real waterdepth = 0.376991; - int StokesOrder = 5; - amrex::Real c0, a11, a22, b22, c2, d2, e2, a31, a33, b31, a42; - amrex::Real a44, b42, b44, c4, d4, e4, a51, a53, a55, b53, b55; + constexpr amrex::Real wavenumber = 2.0; + constexpr amrex::Real water_depth = 0.376991; + constexpr int stokes_order = 5; + amrex::Real c0, a11, a22, b22, c2, a31, a33, b31, a42; + amrex::Real a44, b42, b44, c4, a51, a53, a55, b53, b55; amr_wind::ocean_waves::relaxation_zones::stokes_coefficients( - StokesOrder, wavenumber, waterdepth, c0, a11, a22, b22, c2, d2, e2, a31, - a33, b31, a42, a44, b42, b44, c4, d4, e4, a51, a53, a55, b53, b55); + stokes_order, wavenumber, water_depth, c0, a11, a22, b22, c2, a31, a33, + b31, a42, a44, b42, b44, c4, a51, a53, a55, b53, b55); // Gold coefficient values taken from table 2 of // Fenton, J. Fifth Order Stokes Theory for Steady Waves @@ -45,33 +45,253 @@ TEST_F(WaveTheoriesTest, StokesWaves) const amrex::Real gold_C0 = 0.798448; const amrex::Real gold_C2 = 1.940215; const amrex::Real gold_C4 = -12.970403; - const amrex::Real gold_D2 = -0.626215; - const amrex::Real gold_D4 = 3.257104; - const amrex::Real gold_E2 = 1.781926; - const amrex::Real gold_E4 = -11.573657; - - EXPECT_NEAR(gold_A11, a11, CoeffTol); - EXPECT_NEAR(gold_A22, a22, CoeffTol); - EXPECT_NEAR(gold_A31, a31, CoeffTol); - EXPECT_NEAR(gold_A33, a33, CoeffTol); - EXPECT_NEAR(gold_A42, a42, CoeffTol); - EXPECT_NEAR(gold_A44, a44, CoeffTol); - EXPECT_NEAR(gold_A51, a51, CoeffTol); - EXPECT_NEAR(gold_A53, a53, CoeffTol); - EXPECT_NEAR(gold_A55, a55, CoeffTol); - EXPECT_NEAR(gold_B22, b22, CoeffTol); - EXPECT_NEAR(gold_B31, b31, CoeffTol); - EXPECT_NEAR(gold_B42, b42, CoeffTol); - EXPECT_NEAR(gold_B44, b44, CoeffTol); - EXPECT_NEAR(gold_B53, b53, CoeffTol); - EXPECT_NEAR(gold_B55, b55, CoeffTol); - EXPECT_NEAR(gold_C0, c0, CoeffTol); - EXPECT_NEAR(gold_C2, c2, CoeffTol); - EXPECT_NEAR(gold_C4, c4, CoeffTol); - EXPECT_NEAR(gold_D2, d2, CoeffTol); - EXPECT_NEAR(gold_D4, d4, CoeffTol); - EXPECT_NEAR(gold_E2, e2, CoeffTol); - EXPECT_NEAR(gold_E4, e4, CoeffTol); + + EXPECT_NEAR(gold_A11, a11, coeff_tol); + EXPECT_NEAR(gold_A22, a22, coeff_tol); + EXPECT_NEAR(gold_A31, a31, coeff_tol); + EXPECT_NEAR(gold_A33, a33, coeff_tol); + EXPECT_NEAR(gold_A42, a42, coeff_tol); + EXPECT_NEAR(gold_A44, a44, coeff_tol); + EXPECT_NEAR(gold_A51, a51, coeff_tol); + EXPECT_NEAR(gold_A53, a53, coeff_tol); + EXPECT_NEAR(gold_A55, a55, coeff_tol); + EXPECT_NEAR(gold_B22, b22, coeff_tol); + EXPECT_NEAR(gold_B31, b31, coeff_tol); + EXPECT_NEAR(gold_B42, b42, coeff_tol); + EXPECT_NEAR(gold_B44, b44, coeff_tol); + EXPECT_NEAR(gold_B53, b53, coeff_tol); + EXPECT_NEAR(gold_B55, b55, coeff_tol); + EXPECT_NEAR(gold_C0, c0, coeff_tol); + EXPECT_NEAR(gold_C2, c2, coeff_tol); + EXPECT_NEAR(gold_C4, c4, coeff_tol); +} + +TEST_F(WaveTheoriesTest, StokesWavesFreeSurfaceProfile) +{ + constexpr amrex::Real tol = 1e-4; + + constexpr amrex::Real g = 9.81; + constexpr int stokes_order = 5; + // wavenumber k and water_depth d chosen so that kd = 0.758932 + // to match value of column 3 from table 2 in + // Fenton, J. Fifth Order Stokes Theory for Steady Waves + // Journal of Waterway, Port, Coastal and Ocean Engineering, 1985, 111, + // 216-234 + amrex::Real wavenumber = 2.0; + amrex::Real water_depth = 0.376991; + amrex::Real wavelength = 2. * M_PI / wavenumber; + amrex::Real wave_height = 0.1; + amrex::Real zsl = 0.; + amrex::Real x = 0.; + amrex::Real z = -0.25; + amrex::Real phase_offset = 0.; + amrex::Real time = 0.; + amrex::Real eta = 0.; + amrex::Real u_w = 0.; + amrex::Real v_w = 0.; + amrex::Real w_w = 0.; + + amr_wind::ocean_waves::relaxation_zones::stokes_waves( + stokes_order, wavelength, water_depth, wave_height, zsl, g, x, z, time, + phase_offset, eta, u_w, v_w, w_w); + + // Coefficients values taken from column 3 of table 2 of + // Fenton, J. Fifth Order Stokes Theory for Steady Waves + // Journal of Waterway, Port, Coastal and Ocean Engineering, 1985, 111, + // 216-234 + amrex::Real B22 = 2.502414; + amrex::Real B31 = -5.731666; + amrex::Real B42 = -32.407508; + amrex::Real B44 = 14.033758; + amrex::Real B53 = -103.44536875; + amrex::Real B55 = 37.200027; + + amrex::Real eps = wavenumber * wave_height / 2.; + amrex::Real S = 2. * std::exp(2. * wavenumber * water_depth) / + (std::exp(4. * wavenumber * water_depth) + 1.); + amrex::Real C = 1.0 - S; + amrex::Real C0 = std::sqrt(std::tanh(wavenumber * water_depth)); + amrex::Real C2 = C0 * (2 + 7 * std::pow(S, 2)) / (4 * std::pow(C, 2)); + amrex::Real C4 = C0 * + (4 + 32 * S - 116 * std::pow(S, 2) - 400 * std::pow(S, 3) - + 71 * std::pow(S, 4) + 146 * std::pow(S, 5)) / + (32 * std::pow(C, 5)); + amrex::Real wave_speed = + (C0 + std::pow(eps, 2) * C2 + std::pow(eps, 4) * C4) * + std::sqrt(g / wavenumber); + + amrex::Real omega = wave_speed * wavenumber; + amrex::Real phase = wavenumber * x - omega * time - phase_offset; + + // Check against Eq. (14) from Fenton 1985 + amrex::Real eta_theory = + (eps * std::cos(phase) + std::pow(eps, 2) * B22 * std::cos(2. * phase) + + std::pow(eps, 3) * B31 * (std::cos(phase) - std::cos(3. * phase)) + + std::pow(eps, 4) * + (B42 * std::cos(2. * phase) + B44 * std::cos(4. * phase)) + + std::pow(eps, 5) * + (-(B53 + B55) * std::cos(phase) + B53 * std::cos(3 * phase) + + B55 * std::cos(5 * phase))) / + wavenumber + + zsl; + + EXPECT_NEAR(eta, eta_theory, tol); + + // Re-evaluate with new set of coefficients + // Deep-water limit (k*d->\infty) + water_depth = 100; + wave_height = 0.16; + wavenumber = 0.156; + wavelength = 2. * M_PI / wavenumber; + zsl = 0.; + x = 4.; + z = 0.; + phase_offset = M_PI; + time = 2.7; + eta = 0.; + u_w = 0.; + v_w = 0.; + w_w = 0.; + + amr_wind::ocean_waves::relaxation_zones::stokes_waves( + stokes_order, wavelength, water_depth, wave_height, zsl, g, x, z, time, + phase_offset, eta, u_w, v_w, w_w); + + // Coefficients values taken from column 1 of table 2 of + // Fenton, J. Fifth Order Stokes Theory for Steady Waves + // Journal of Waterway, Port, Coastal and Ocean Engineering, 1985, 111, + // 216-234 + B22 = 0.5; + B31 = -0.375; + B42 = 0.3333333; + B44 = 0.3333333; + B53 = 0.7734375; + B55 = 0.3255208; + + eps = wavenumber * wave_height / 2.; + // Coefficients computed analytically by taking limit kd -> \infty and + // simplified accordingly + // Note that in this limit S = 0 and C = 1 and thus they are omitted here + C0 = 1.0; + C2 = 0.5; + C4 = 0.125; + wave_speed = (C0 + std::pow(eps, 2) * C2 + std::pow(eps, 4) * C4) * + std::sqrt(g / wavenumber); + + omega = wave_speed * wavenumber; + phase = wavenumber * x - omega * time - phase_offset; + + // Matches Eq. (18) from Fenton 1985 + eta_theory = + (eps * std::cos(phase) + std::pow(eps, 2) * B22 * std::cos(2. * phase) + + std::pow(eps, 3) * B31 * (std::cos(phase) - std::cos(3. * phase)) + + std::pow(eps, 4) * + (B42 * std::cos(2. * phase) + B44 * std::cos(4. * phase)) + + std::pow(eps, 5) * + (-(B53 + B55) * std::cos(phase) + B53 * std::cos(3 * phase) + + B55 * std::cos(5 * phase))) / + wavenumber + + zsl; + + EXPECT_NEAR(eta, eta_theory, tol); +} + +TEST_F(WaveTheoriesTest, StokesWavesVelocityComponents) +{ + constexpr amrex::Real tol = 1e-4; + + constexpr amrex::Real g = 9.81; + constexpr int stokes_order = 5; + // wavenumber k and water_depth d chosen so that kd = 0.758932 + // to match value of column 3 from table 2 in + // Fenton, J. Fifth Order Stokes Theory for Steady Waves + // Journal of Waterway, Port, Coastal and Ocean Engineering, 1985, 111, + // 216-234 + constexpr amrex::Real wavenumber = 2.0; + constexpr amrex::Real water_depth = 0.376991; + constexpr amrex::Real wavelength = 2. * M_PI / wavenumber; + constexpr amrex::Real wave_height = 0.1; + constexpr amrex::Real zsl = 0.; + constexpr amrex::Real x = 0.; + constexpr amrex::Real z = -0.25; + constexpr amrex::Real phase_offset = 0.; + constexpr amrex::Real time = 0.; + amrex::Real eta = 0.; + amrex::Real u_w = 0.; + amrex::Real v_w = 0.; + amrex::Real w_w = 0.; + + amr_wind::ocean_waves::relaxation_zones::stokes_waves( + stokes_order, wavelength, water_depth, wave_height, zsl, g, x, z, time, + phase_offset, eta, u_w, v_w, w_w); + + // Coefficients values taken from column 3 of table 2 of + // Fenton, J. Fifth Order Stokes Theory for Steady Waves + // Journal of Waterway, Port, Coastal and Ocean Engineering, 1985, 111, + // 216-234 + const amrex::Real A11 = 1.208490; + const amrex::Real A22 = 0.799840; + const amrex::Real A31 = -9.105340; + const amrex::Real A33 = 0.368275; + const amrex::Real A42 = -12.196150; + const amrex::Real A44 = 0.058723; + const amrex::Real A51 = 108.46831725; + const amrex::Real A53 = -6.941756; + const amrex::Real A55 = -0.074979; + + const amrex::Real eps = wavenumber * wave_height / 2.; + const amrex::Real S = 2. * std::exp(2. * wavenumber * water_depth) / + (std::exp(4. * wavenumber * water_depth) + 1.); + const amrex::Real C = 1.0 - S; + const amrex::Real C0 = std::sqrt(std::tanh(wavenumber * water_depth)); + const amrex::Real C2 = C0 * (2 + 7 * std::pow(S, 2)) / (4 * std::pow(C, 2)); + const amrex::Real C4 = + C0 * + (4 + 32 * S - 116 * std::pow(S, 2) - 400 * std::pow(S, 3) - + 71 * std::pow(S, 4) + 146 * std::pow(S, 5)) / + (32 * std::pow(C, 5)); + const amrex::Real wave_speed = + (C0 + std::pow(eps, 2) * C2 + std::pow(eps, 4) * C4) * + std::sqrt(g / wavenumber); + + const amrex::Real omega = wave_speed * wavenumber; + const amrex::Real phase = wavenumber * x - omega * time - phase_offset; + + // Compare with theoretical Results from Kinnas + // https://www.sciencedirect.com/science/article/pii/S0029801817306066 + // Define coefficients using Eq.(19) + amrex::Vector a(stokes_order); + a[0] = A11 + (eps * eps) * A31 + std::pow(eps, 4) * A51; + a[1] = A22 + (eps * eps) * A42; + a[2] = A33 + (eps * eps) * A53; + a[3] = A44; + a[4] = A55; + + // Horizontal velocity from Eq.(21) and vertical velocity from Eq.(23) in + // Kinnas + // Initialize to first order terms + amrex::Real horizontal_velocity = + eps * a[0] * std::cosh(wavenumber * (water_depth + (z - zsl))) * + std::cos(phase); + amrex::Real vertical_velocity = + eps * a[0] * std::sinh(wavenumber * (water_depth + (z - zsl))) * + std::sin(phase); + + for (int n = 1; n < stokes_order; ++n) { + horizontal_velocity += + std::pow(eps, n + 1) * (n + 1) * a[n] * + std::cosh((n + 1) * wavenumber * (water_depth + (z - zsl))) * + std::cos((n + 1) * phase); + vertical_velocity += + std::pow(eps, n + 1) * (n + 1) * a[n] * + std::sinh((n + 1) * wavenumber * (water_depth + (z - zsl))) * + std::sin((n + 1) * phase); + } + horizontal_velocity *= (C0 * std::sqrt(g / wavenumber)); + vertical_velocity *= (C0 * std::sqrt(g / wavenumber)); + + EXPECT_NEAR(u_w, horizontal_velocity, tol); + EXPECT_NEAR(w_w, vertical_velocity, tol); } TEST_F(WaveTheoriesTest, StokesWaveLength) @@ -93,14 +313,14 @@ TEST_F(WaveTheoriesTest, StokesWaveLength) wave_period, water_depth, wave_height, wave_order, g, tol_lambda, iter_max); - const amrex::Real k_Newton = 2.0 * M_PI / lambda; + const amrex::Real k_newton = 2.0 * M_PI / lambda; // Compare with expected wavenumber from theory k = omega^2/g, where omega = // 2Pi/wave_period const amrex::Real k_theory = (2.0 * M_PI / wave_period) * (2.0 * M_PI / wave_period) / g; - EXPECT_NEAR(k_Newton, k_theory, 1e-8); + EXPECT_NEAR(k_newton, k_theory, 1e-8); // Check wave theory wave_height = 0.2;