Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Improvements to the structure of stokes_waves_K.H #1403

Merged
merged 24 commits into from
Dec 19, 2024
Merged

Improvements to the structure of stokes_waves_K.H #1403

Show file tree
Hide file tree
Changes from 15 commits
Commits
Show all changes
24 commits
Select commit Hold shift + click to select a range
6b4d1a9
adding expression for C constant in stokes_coefficients()
mpolimeno Dec 11, 2024
4a06f9a
started replacing '1-S' with C in the relevant places
mpolimeno Dec 11, 2024
f1deec1
Merge branch 'main' into stokes_waves_K_header_improvements
mpolimeno Dec 11, 2024
de70a86
adding comments with relevant references for Stokes wave theory
mpolimeno Dec 12, 2024
467e5f1
Merge branch 'main' into stokes_waves_K_header_improvements
jrood-nrel Dec 12, 2024
5b4688f
made definition of S consistent across functions; changed 1-S to C in…
mpolimeno Dec 12, 2024
4b4d9c1
removing unused coefficients
mpolimeno Dec 12, 2024
397d508
adding upper bound for inputs of cosh and sinh to ensure finite veloc…
mpolimeno Dec 17, 2024
17adbca
draft unit test for free surface profile
mpolimeno Dec 17, 2024
d3e6035
fix gpu compilation issue
mpolimeno Dec 17, 2024
f4128fe
more gpu-related fixes
mpolimeno Dec 17, 2024
7bafa75
Merge branch 'main' into stokes_waves_K_header_improvements
mpolimeno Dec 17, 2024
f5135a4
using snake_case for variable names when appropriate, following amrwi…
mpolimeno Dec 17, 2024
360f271
draft unit test for velocity components
mpolimeno Dec 17, 2024
5f27f99
expanding unit test for velocity components to be able to test all re…
mpolimeno Dec 18, 2024
aa287d5
temporary fix
mpolimeno Dec 18, 2024
7005841
updated formulation for the computation of the velocity components
mpolimeno Dec 19, 2024
246b93e
fix for gpu compilation
mpolimeno Dec 19, 2024
cd422ec
Merge branch 'main' into stokes_waves_K_header_improvements
mpolimeno Dec 19, 2024
aa0d01f
Update unit_tests/ocean_waves/test_wave_theories.cpp
mpolimeno Dec 19, 2024
82f56af
Update unit_tests/ocean_waves/test_wave_theories.cpp
mpolimeno Dec 19, 2024
430d985
Update unit_tests/ocean_waves/test_wave_theories.cpp
mpolimeno Dec 19, 2024
7d6c416
constexpr for variables whose value will not change
mpolimeno Dec 19, 2024
aa9afde
Merge branch 'main' into stokes_waves_K_header_improvements
mbkuhn Dec 19, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
97 changes: 52 additions & 45 deletions amr-wind/ocean_waves/relaxation_zones/stokes_waves_K.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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));
Expand Down Expand Up @@ -131,8 +143,6 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void stokes_coefficients(
amrex::Real& a22,
amrex::Real& b22,
amrex::Real& c2,
amrex::Real& d2,
amrex::Real& e2,
amrex::Real& a31,
amrex::Real& a33,
amrex::Real& b31,
Expand All @@ -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,
Expand All @@ -155,59 +163,54 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void stokes_coefficients(
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));
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 (StokesOrder == 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));
(8 * std::pow(C, 3));
if (StokesOrder == 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));
(32 * std::pow(C, 5));
if (StokesOrder == 4) {
return;
}
Expand All @@ -217,26 +220,26 @@ 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));
(384 * (3 + 2 * S) * (4 + S) * std::pow(C, 6));
if (StokesOrder == 5) {
return;
}
Expand All @@ -247,7 +250,7 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void stokes_coefficients(
}
}

AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real stokes_sinh_sin(
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real stokes_sinh_sin(
int m,
int n,
amrex::Real phase,
Expand All @@ -266,6 +269,11 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real stokes_sinh_sin(
amrex::Real zsl,
amrex::Real z)
{
amrex::Real nkdz = n * wavenumber * (waterdepth + (z - zsl));
if (nkdz > 50. * M_PI) {
nkdz = 50. * M_PI;
}

amrex::Real D = 0.0;
if (m == 1 && n == 1) {
D = a11;
Expand Down Expand Up @@ -295,12 +303,11 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real stokes_sinh_sin(
D = a55;
}

return std::pow(eps, m) * D * n * wavenumber *
std::sinh(n * wavenumber * (waterdepth + (z - zsl))) *
return std::pow(eps, m) * D * n * wavenumber * std::sinh(nkdz) *
std::sin(n * phase);
}

AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real stokes_cosh_cos(
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real stokes_cosh_cos(
int m,
int n,
amrex::Real phase,
Expand All @@ -319,6 +326,11 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real stokes_cosh_cos(
amrex::Real zsl,
amrex::Real z)
{
amrex::Real nkdz = n * wavenumber * (waterdepth + (z - zsl));
if (nkdz > 50. * M_PI) {
nkdz = 50. * M_PI;
}

amrex::Real D = 0.0;
if (m == 1 && n == 1) {
D = a11;
Expand Down Expand Up @@ -348,13 +360,12 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real stokes_cosh_cos(
D = a55;
}

return std::pow(eps, m) * D * n * wavenumber *
std::cosh(n * wavenumber * (waterdepth + (z - zsl))) *
return std::pow(eps, m) * D * n * wavenumber * std::cosh(nkdz) *
std::cos(n * phase);
}

// Based on Fenton 1985
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void stokes_waves(
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void stokes_waves(
int StokesOrder,
amrex::Real wavelength,
amrex::Real waterdepth,
Expand All @@ -374,22 +385,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);
StokesOrder, wavenumber, waterdepth, 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 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;
Expand Down
Loading