GitHub

Deep-Water

In this notebook, we compare several models for the propagation of surface gravity waves over an infinite fluid layer. Specifically, we consider

Initialization

Import package

using WaterWaves1Dimport LinearAlgebra.norm

Define initial data

η(x) = exp.(-x .^ 2); # Gaussian initial data for the surface deformation.v(x) = zero(x)      # we set the initial velocity as zero to avoid inconsistencies among different models.init = Init(η, v);

The quadratic spectral model (WW2)

We investigate here the spectral model with quadratic nonlinearity (WW2) proposed by Dommermuth and Yue, West et al., and Craig and Sulem. We show that the model exhibits high-frequency instabilities, unless a suitable regularization is included.

First experiment with a small number of modes

param_stable = (    # Physical parameters. Variables are non-dimensionalized as in Lannes, The water waves problem, isbn:978-0-8218-9470-5    ϵ = 1 / 4,    # nonlinearity dimensionless parameter    μ = Inf,    # inifinite-depth case    # Numerical parameters    N = 2^9,    # number of collocation points    L = 2 * π,    # half-length of the numerical tank (-L,L)    T = 5,      # final time of computation    dt = 0.001,  # timestep);WW2_stable = Problem(WWn(param_stable, dealias = 1), init, param_stable)solve!(WW2_stable);
┌ Info: Now solving the initial-value problem WW2
│ with timestep dt=0.001, final time T=5.0,
└ and N=512 collocation points.

Plot solution at final time

using Plotsplot(WW2_stable; var = [:surface, :fourier])
Example block output

Second experiment with more modes

param_unstable = (    # Physical parameters. Variables are non-dimensionalized as in Lannes, The water waves problem, isbn:978-0-8218-9470-5    ϵ = 1 / 4,    # nonlinearity dimensionless parameter    μ = Inf,    # inifinite-depth case    # Numerical parameters    N = 2^11,   # number of collocation points    L = 2 * π,    # half-length of the numerical tank (-L,L)    T = 1.5,      # final time of computation    dt = 0.001,  # timestep);WW2_unstable = Problem(WWn(param_unstable, dealias = 1), init, param_unstable)solve!(WW2_unstable);
┌ Info: Now solving the initial-value problem WW2
│ with timestep dt=0.001, final time T=1.5,
└ and N=2048 collocation points.

Plot solution just before breakup

plot(WW2_unstable; T = 1.1, var = [:surface, :fourier])
Example block output

The (very) high-frequency component of the flow, which is generated by machine-precision rounding errors, is rapidly amplified and produces spurious numerical oscillations. This phenomenon is not due to aliasing, since a dealiasing filter (adapted to quadratic nonlinearities) has been used. Although we use an explicit time integrator (RK4) is not due to CFL-type conditions, since it is stable when diminishing the timestep.

The rectified quadratic spectral model

We use here the rectified WW2 model introduced by Duchêne and Melinand

WW2_rectified = Problem(WWn(param_unstable, dealias = 1, δ = 0.01), init, param_unstable, label = "rWW2") # δ is the strength of the regularizationsolve!(WW2_rectified);
┌ Info: Now solving the initial-value problem rWW2
│ with timestep dt=0.001, final time T=1.5,
└ and N=2048 collocation points.

Plot solution at final time

plot(WW2_rectified, var = [:surface, :fourier])
Example block output

The solution shows no more sign of instabilities. And this is so for much higher number of modes (provided the timestep is sufficiently small).

Plot the difference between the solutions to the rectified and non-rectified WW2 models.

plot([WW2_rectified, WW2_unstable], T = 1, var = [:difference, :difference_fourier])
Example block output

By plotting the difference, we see that the rectification suppresses high-frequency instabilities without affecting the main profile of the wave determined by small wavenumbers.

Disintegration of a heap of water

We now compare

with the "exact" water waves equation, with an initial data consisting of an excess of mass, with zero initial velocity.

param_heap = (    # Physical parameters. Variables are non-dimensionalized as in Lannes, The water waves problem, isbn:978-0-8218-9470-5    ϵ = 1 / 4,    # nonlinearity dimensionless parameter    μ = Inf,    # infinite-depth case    # Numerical parameters    N = 2^9,   # number of collocation points    L = 2 * π,    # half-length of the numerical tank (-L,L)    T = 5,      # final time of computation    dt = 0.001,  # timestep)init = Init(x -> exp.(-x .^ 2), x -> zero(x));WW_heap = Problem(WaterWaves(param_heap, dealias = 1, verbose = false), init, param_heap);rWW2_heap = Problem(WWn(param_heap, dealias = 1, δ = 0.02), init, param_heap, label = "rWW2") # δ is the strength of the regularizationMatsuno_heap = Problem(Matsuno(param_heap, dealias = 1), init, param_heap)AkersNicholls_heap = Problem(AkersNicholls(param_heap, dealias = 1), init, param_heap)solve!(WW_heap);solve!(rWW2_heap);solve!(Matsuno_heap);solve!(AkersNicholls_heap);plot([WW_heap, rWW2_heap, Matsuno_heap, AkersNicholls_heap])plt = plot([(WW_heap, rWW2_heap), (WW_heap, Matsuno_heap), (WW_heap, AkersNicholls_heap)], var = :difference)display(plt)WWη, WWv, WWx, = solution(WW_heap)rWW2η, = solution(rWW2_heap, x = WWx);errWW2 = norm(WWη - rWW2η) / sqrt(length(WWx))Matη, = solution(Matsuno_heap, x = WWx);errMat = norm(WWη - Matη) / sqrt(length(WWx))ANη, = solution(AkersNicholls_heap, x = WWx);errAN = norm(WWη - ANη) / sqrt(length(WWx))print("Errors (in l²). rectified WW2: $errWW2, Matsuno: $errMat, Akers&Nicholls: $errAN")
┌ Info: Now solving the initial-value problem water waves
│ with timestep dt=0.001, final time T=5.0,
└ and N=512 collocation points.

┌ Info: Now solving the initial-value problem rWW2
│ with timestep dt=0.001, final time T=5.0,
└ and N=512 collocation points.

┌ Info: Now solving the initial-value problem Matsuno
│ with timestep dt=0.001, final time T=5.0,
└ and N=512 collocation points.

┌ Info: Now solving the initial-value problem Akers-Nicholls
│ with timestep dt=0.001, final time T=5.0,
└ and N=512 collocation points.

Errors (in l²). rectified WW2: 0.00200876373700181, Matsuno: 0.000918672651668852, Akers&Nicholls: 0.0012203806493777008

Random initial data

We now compare

with the "exact" water waves equation, with an initial data consisting of a random initial surface perturbation, with zero initial velocity.

param_rand = (    # Physical parameters. Variables are non-dimensionalized as in Lannes, The water waves problem, isbn:978-0-8218-9470-5    ϵ = 1 / 4,    # nonlinearity dimensionless parameter    μ = Inf,    # infinite-depth case    # Numerical parameters    N = 2^10,   # number of collocation points    L = 2 * π,    # half-length of the numerical tank (-L,L)    T = 1,      # final time of computation    dt = 0.001,  # timestep)mesh = Mesh(param_rand);x = mesh.x;init = Random(x; L = 2, a = (1, 0)); # we set the initial velocity as zero to avoid possible inconsistencies among different models.WW_rand = Problem(WaterWaves(param_rand, dealias = 1, verbose = false), init, param_rand);rWW2_rand = Problem(WWn(param_rand, dealias = 1, δ = 0.05), init, param_rand, label = "rWW2");Matsuno_rand = Problem(Matsuno(param_rand, dealias = 1), init, param_rand);AkersNicholls_rand = Problem(AkersNicholls(param_rand, dealias = 1), init, param_rand);solve!(WW_rand);solve!(rWW2_rand);solve!(Matsuno_rand);solve!(AkersNicholls_rand);plot([WW_rand, rWW2_rand, Matsuno_rand, AkersNicholls_rand])plt = plot([(WW_rand, rWW2_rand), (WW_rand, Matsuno_rand), (WW_rand, AkersNicholls_rand)], var = :difference)display(plt)WWη, WWv, WWx, = solution(WW_rand)rWW2η, = solution(rWW2_rand, x = WWx);errWW2 = norm(WWη - rWW2η) / sqrt(length(WWx))Matη, = solution(Matsuno_rand, x = WWx);errMat = norm(WWη - Matη) / sqrt(length(WWx))ANη, = solution(AkersNicholls_rand, x = WWx);errAN = norm(WWη - ANη) / sqrt(length(WWx))print("Errors (in l²). rectified WW2: $errWW2, Matsuno: $errMat, Akers&Nicholls: $errAN")
┌ Info: Now solving the initial-value problem water waves
│ with timestep dt=0.001, final time T=1.0,
└ and N=1024 collocation points.

┌ Info: Now solving the initial-value problem rWW2
│ with timestep dt=0.001, final time T=1.0,
└ and N=1024 collocation points.

┌ Info: Now solving the initial-value problem Matsuno
│ with timestep dt=0.001, final time T=1.0,
└ and N=1024 collocation points.

┌ Info: Now solving the initial-value problem Akers-Nicholls
│ with timestep dt=0.001, final time T=1.0,
└ and N=1024 collocation points.

Errors (in l²). rectified WW2: 0.0001451131388095949, Matsuno: 0.00025198532431141576, Akers&Nicholls: 0.0001646263107059437

This page was generated using Literate.jl.