Numerical methods

LogRange(kmin,kmax,M)

Create M logarithmic grid points between kmin and kmax.

area_ratio(yi,yu,yl,yr)

Area of a box (x,y) ∈ [xl,xr] × [yi,yu] that is below the line passing by (xl,yl) and (xr,yr), divided by the total area of the box.

Note: does not depend on xl and xr.

area_ratio_grid(x, y, f)

For each cell of the grid x × y, return the area under f divided by the total area of the box.

Note: convolution allows to integrate over y < f.

area_ratio_logbins(kkh, kkz, logλz, f)

Same as area_ratio_grid but for logarithmic grids.

• kkh: horizontal wave vectors grid points
• kkz: vertical wave vectors grid points
• logλz: log(λz) where λz is the logarithmic increment of the vertical wave vectors grid
• f: values of f vs kkh

Note: convolution allows to integrate over kkz < f.

change_mesh!(Nknew,Nkold;interp_scheeme=BS_interp)

Fill Nknew wave action spectrum from interpolation interp_scheeme of Nkold.

meshgrid(x, y)

Create 2D grids from the 1D vectors x and y.

BS_interp(kk::AbstractVector)

Allocate nodes for BSplines interpolation. See val_nk.

bilin_interp_khkz(kkh::AbstractVector, kkz::AbstractVector)

Allocate coefficients for bilinear interpolation and define interpolation method. See val_nk.

bilinlog_interp_khkz(kkh::AbstractVector, kkz::AbstractVector)

Allocate coefficients for bilinear interpolation in log scale and define interpolation method. See val_nk.

cpow_interp_khkz(kkh::AbstractVector, kkz::AbstractVector)

Allocate coefficients for cpow interpolation and define interpolation method. See val_nk.

exp_interp_khkz(kkh::AbstractVector, kkz::AbstractVector)

Allocate coefficients for exponential interpolation and define interpolation method. See val_nk.

lin_interp(kk::AbstractVector)

Allocate coefficients for linear interpolation and define interpolation method. See val_nk.

linlog_interp(kk::AbstractVector)

Allocate coefficients for linear interpolation in log scale (power law) and define interpolation method. See val_nk.

powGauss_interp(kk::AbstractVector)

Allocate coefficients for power-law Gaussian interpolation. See val_nk.

powexp_interp(kk::AbstractVector)

Allocate coefficients for power-law exponential interpolation. See val_nk.

bilin_interp_vs(kh,kz,khinf,khsup,kzinf,kzsup,n1,n2,n3,n4)

Return the bilinear interpolation of $n_{\bf k}$ at point (kh,kz) computed using the points (khinf,kzinf,n1), (khsup,kzinf,n2), (khsup,kzsup,n3) and (khinf,kzsup,n4).

lin_interp_vs(k,k1,k2,n1,n2)

Return the linear interpolation of $n_k$ at point k computed using the points (k1,n1) and (k2,n2)

$n_k = n_1 \frac{k_2-k}{k_2-k_1} + n_2 \frac{k-k_1}{k_2-k_1}.$

update_coeff_interp!(interp::BS_interp,Nk::wave_spectrum)

Compute the coefficients for B-spline interpolation. The method is given by the BSplineKit package.

• Nk: Structure containing the wave action wave_spectrum
• interp: structure where the coefficients are stored

update_coeff_interp!(::bilin_interp_khkz, Nk::wave_spectrum_khkz)

Compute the coefficients $C_0$, $\alpha_h$, $\alpha_z$ and $\beta$ for bilinear interpolation such that $n_{\bf k} = C_0 - \alpha_h k_h - \alpha_z k_z - \beta k_h k_z$.

• Nk: Structure containing the wave action wave_spectrum
• interp: structure where the coefficients are stored

update_coeff_interp!(::bilinlog_interp_khkz, Nk::wave_spectrum_khkz)

Compute the coefficients $C_0$, $\alpha_h$, $\alpha_z$ and $\beta$ for bilinear interpolation in logarithmic scale (power law) such that $n_{\bf k} = e^{C_0 - \alpha_h \log k_h - \alpha_z \log k_z - \beta \log k_h \log k_z}$.

• Nk: Structure containing the wave action wave_spectrum
• interp: structure where the coefficients are stored

update_coeff_interp!(::cpow_interp_khkz, Nk::wave_spectrum_khkz)

Compute the coefficients $C_0$, $\alpha_h$, $\alpha_z$ and $\beta$ for bilinear interpolation such that $n_{\bf k} = \beta + C_0 k_h^{-\alpha_h} k_z^{-\alpha_z k_z}$.

• Nk: Structure containing the wave action wave_spectrum
• interp: structure where the coefficients are stored

update_coeff_interp!(::exp_interp_khkz, Nk::wave_spectrum_khkz)

Compute the coefficients $C_0$, $\alpha_h$, $\alpha_z$ and $\beta$ for exponential interpolation such that $n_{\bf k} = e^{C_0 - \alpha_h k_h - \alpha_z k_z - \beta k_h k_z}$.

• Nk: Structure containing the wave action wave_spectrum
• interp: structure where the coefficients are stored

update_coeff_interp!(::lin_interp,Nk::wave_spectrum)

Compute the coefficients $\alpha$ and $\beta$ for linear interpolation such that $n_{\bf k} = \alpha + \beta k$.

• Nk: Structure containing the wave action wave_spectrum
• interp: structure where the coefficients are stored

update_coeff_interp!(::linlog_interp,Nk::wave_spectrum)

Compute the coefficients $\alpha$ and $C_0$ for linear interpolation in logarithmic scale (power law) such that $n_{\bf k} = C_0 k^{-\alpha} = C_0 e^{- \alpha \log k}$.

• Nk: Structure containing the wave action wave_spectrum
• interp: structure where the coefficients are stored

update_coeff_interp!(::powexp_interp,Nk::wave_spectrum)

Compute the coefficients $C_0$, $\alpha$ and $\beta$ for power law Gaussian interpolation such that $n_{\bf k} = C_0 k^{-\alpha} e^{-\beta k^2}$.

• Nk: Structure containing the wave action wave_spectrum
• interp: structure where the coefficients are stored

update_coeff_interp!(::powexp_interp,Nk::wave_spectrum)

Compute the coefficients $C_0$, $\alpha$ and $\beta$ for power law exponential interpolation such that $n_{\bf k} = C_0 k^{-\alpha} e^{-\beta k}$.

• Nk: Structure containing the wave action wave_spectrum
• interp: structure where the coefficients are stored

val_nk(::bilin_interp_khkz, Nk::wave_spectrum_khkz, kh, kz)

Return the value of Nk.nk interpolated at (kh, kz) (with bilin_interp_khkz). Coefficients should be updated.

See also update_coeff_interp!

val_nk(::bilinlog_interp_khkz, Nk::wave_spectrum_khkz, kh, kz)

Return the value of Nk.nk interpolated at (kh, kz) (with bilinlog_interp_khkz). Coefficients should be updated.

See also update_coeff_interp!

val_nk(::cpow_interp_khkz, Nk::wave_spectrum_khkz, kh, kz)

Return the value of Nk.nk interpolated at (kh, kz) (with cpow_interp_khkz). Coefficients should be updated.

See also update_coeff_interp!

val_nk(::exp_interp_khkz, Nk::wave_spectrum_khkz, kh, kz)

Return the value of Nk.nk interpolated at (kh, kz) (with exp_interp_khkz). Coefficients should be updated.

See also update_coeff_interp!

val_nk(::lin_interp,Nk::waveaction,k)

Return the value of Nk.nk interpolated at k (with lin_interp). Coefficients should be updated.

See also update_coeff_interp!

val_nk(::linlog_interp,Nk::waveaction,k)

Return the value of Nk.nk interpolated at k (with linlog_interp). Coefficients should be updated.

See also update_coeff_interp!

val_nk(interp::powGauss_interp, Nk::wave_spectrum, k)

Return the value of Nk.nk interpolated at k (with powGauss_interp). Coefficients should be updated.

See also update_coeff_interp!

val_nk(::powexp_interp,Nk::waveaction,k)

Return the value of Nk.nk interpolated at k (with powexp_interp). Coefficients should be updated.

See also update_coeff_interp!

integrate_with_cpow_khkz()

Select the use of the cpow integration with log bins khkz method. See integrate.

integrate_with_log_bins_khkz()

Select the use of the trapezoidal integration with log bins khkz method. See integrate.

clean_waveaction!(Nk::waveaction)

Set the minimum of the waveaction spectrum Nk.nk to min_nk. By default min_nk=0 so it removes negative values.

cumintegrate!(CumInt::AbstractVector, M::Int, λ::Float64, logλ::Float64, kk::Vector{Float64}, F::Vector{Float64})

Compute cumulative integral of F using log bins and trapezoidal rule. The result is stored in CumInt.

• M: number of grid points
• λ: logarithmic increment of the grid
• logλ: log(λ)
• kk: grid

See also integrate_with_log_bins.

cumintegrate!(CumInt::AbstractVector, Sk::wave_spectrum)

Compute cumulative integral of Sk.Nk using log bins and trapezoidal rule. The result is stored in CumInt.

cumintegraterev!(CumInt::AbstractVector, M::Int, λ::Float64, logλ::Float64, kk::Vector{Float64}, F::Vector{Float64})

Compute cumulative integral using log bins in reverse order.

cumintegraterev!(CumInt::AbstractVector, Sk::wave_spectrum)

Compute cumulative integral using log bins in reverse order.

integrate(Nk::wave_spectrum, imin::Int=1, imax::Int=-1)

Return the integral of Nk.nk using (simple) trapezoidal rule from index imin to imax. By default, it integrates over the whole domain.

integrate_with_grid(Sk::field_grid_1D, imin::Int=1, imax::Int=-1)

Return the integral of Sk.F using (simple) trapezoidal rule on grid Sk.kk between index imin and imax.

integrate(integ::integrate_with_cpow_khkz, Nk::wave_spectrum_khkz, ihmin::Int=1, ihmax::Int=-1, izmin::Int=1, izmax::Int=-1)

Return the integral of Nk.nk assuming $n_{\bf k} = \beta + C_0 k_h^{-\alpha_h} k_z^{-\alpha_z k_z}$ in each mesh from index ihmin to ihmax and izmin to izmax. By default, it integrates over the whole domain.

See also cpow_interp_khkz.

integrate(integ::integrate_with_log_bins_khkz, Nk::wave_spectrum_khkz, ihmin::Int=1, ihmax::Int=-1, izmin::Int=1, izmax::Int=-1)

Return the integral of Nk.nk using log bins and trapezoidal rule from index ihmin to ihmax and izmin to izmax. By default, it integrates over the whole domain.

Namely $I = \int\limits_{k_{h{\rm min}}}^{k_{h{\rm max}}} \int\limits_{k_{z{\rm min}}}^{k_{z{\rm max}}} ~ F(k_h,k_z) ~ \mathrm{d}k_h \mathrm{d}k_z = \int\limits_{i_{h{\rm min}}}^{i_{h{\rm max}}} \int\limits_{i_{z{\rm min}}}^{i_{z{\rm max}}} ~ F(k_h=\lambda_h^{i_h}, k_z=\lambda_z^{i_z}) ~ \lambda_h^{i_h} \lambda_z^{i_z} \ln \lambda_h \ln \lambda_z ~ \mathrm{d} i_h \mathrm{d} i_z$.

integrate_quad_with_log_bins(Nk::wave_spectrum, imin::Int=1, imax::Int=-1)

Return the integral of Nk.nk using log bins and Simpson rule from index imin to imax. By default, it integrates over the whole domain.

integrate_with_grid(x, F, imin::Int=1, imax::Int=-1)

Return the integral of F using (simple) trapezoidal rule on grid x between index imin and imax.

integrate_with_log_bins(Nk::wave_spectrum, imin::Int=1, imax::Int=-1)

Return the integral of Nk.nk using log bins and trapezoidal rule from index imin to imax. By default, it integrates over the whole domain.

integrate_with_log_bins(M::Int, λ::Float64, logλ::Float64, kk::Vector{Float64}, F::Vector{Float64}, imin::Int=1, imax::Int=-1)

Return the integral of F using log bins and trapezoidal rule from index imin to imax.

• M: number of grid points
• λ: logarithmic increment of the grid
• logλ: log(λ)
• kk: grid

By default, it integrates over the whole domain. Namely $I = \int\limits_{k_{\rm min}}^{k_{\rm max}} ~ F(k) ~ \mathrm{d}k = \int\limits_{i_{\rm min}}^{i_{\rm max}} ~ F(k=\lambda^{i} ) \lambda^{i} \ln \lambda ~ \mathrm{d} i$. It uses trapezoidal rule.

integrate_with_log_bins_segment(Fa, Fb, ka, kb)

Return the integral of F using trapezoidal log bins from ka to kb.

• Fa: value of F at ka
• Fb: value of F at kb

AB2_RK2_step()

Set the use of split Second-order-Exponential-Adams-Bashforth and Runge-Kutta 2. See advance!.

AB2_RK2arrays

Structure for split Second-order-Exponential-Adams-Bashforth and Runge-Kutta 2

AB_Euler_step()

Set the use of split Exponential-Adams-Bashforth and Euler. See advance!.

ETD2_step()

Set the use of ETD2. See advance!.

ETD2arrays

Structure for ETD2

ETD4_step()

Set the use of ETD4. See advance!.

ETD4arrays

Structure for ETD4

Euler_step()

Set the use of a simple Euler. See advance!.

RK2_step()

Set the use of a Runge-Kutta 2. See advance!.

RK2arrays

Structure for standard Runge-Kutta 2

RK4_step()

Set the use of 4th-order Runge-Kutta 4, no forcing nor dissipation. See advance!.

RK4arrays

Structure for standard Runge-Kutta 4, no forcing nor dissipation

adaptative_time_step(Run, dtmin, dtmax, dt, cmin=0.05, cmax=0.5)

Return the adaptated time step based on the wave action spectrum.

• Run: run structure.
• dtmin: minimal time step.
• dtmax: maximal time step.
• dt: current time step.
• cmin: lower threshold for the ratio of the nonlinear time step and current time step.
• cmax: upper threshold for the ratio of the nonlinear time step and current time step.

See also get_T_nonlinear.

advance!(::Euler_step, Run, dt)

Make one time step with Euler_step.

• Run: run structure
• dt: time step

The time advancement is simply $n_{\bf k}(t+\mathrm{d}t) = \left[ St_{\rm k} + f_{\bf k} - d_{\bf k} n_{\bf k} \right] ~ \mathrm{d}t$.

advance!(::RK2_step, Run, dt)

Make one time step with RK2_step.

• Run: run structure
• dt: time step

The time advancement is done with Runge-Kutta 2 method.

advance!(::RK4_step, Run, dt)

Make one time step with RK4_step.

• Run: run structure
• dt: time step

The time advancement is done with Runge-Kutta 4 method.

advance!(::AB2_RK2_step, Run, dt)

Make one time step with AB2_RK2_step.

• Run: run structure
• dt: time step

The collision integral is decomposed as $St_{\bf k} = - \gamma_{\bf k} n_{\bf k} + \eta_{\bf k}$. Time advancement is done with a splitting method, treating the term $-(\gamma_{\bf k} + d_{\bf k}) n_{\bf k}$ implicitly and the term $\eta_{\bf k} + f_{\bf k}$ explicitly with Runge-Kutta 2 method.

Note: The Run structure must have the coefficients $\gamma_{\bf k}$ and $\eta_{\bf k}$ that must be computed in St_k!. This scheeme is not implemented for all physical systems.

advance!(::AB_Euler_step, Run, dt)

Make one time step with AB_Euler_step.

• Run: run structure
• dt: time step

The collision integral is decomposed as $St_{\bf k} = - \gamma_{\bf k} n_{\bf k} + \eta_{\bf k}$. Time advancement is done with a splitting method, treating the term $-(\gamma_{\bf k} + d_{\bf k}) n_{\bf k}$ implicitly and the term $\eta_{\bf k} + f_{\bf k}$ explicitly with Euler method.

Note: The Run structure must have the coefficients $\gamma_{\bf k}$ and $\eta_{\bf k}$ that must be computed in St_k!. This scheeme is not implemented for all physical systems.

advance!(::ETD2_step, Run, dt)

Make one time step with ETD2_step.

• Run: run structure
• dt: time step

advance!(::ETD4_step, Run, dt)

Make one time step with ETD4_step.

• Run: run structure
• dt: time step

get_T_nonlinear(Run)

Non linear time based on the wave action spectrum

$\tau_{\rm nl} = \frac{1}{\max\limits_{\bf k} \left| \frac{St_{\bf k}}{n_{\bf k}} \right| }.$

We use only the energetic modes with $n_{\bf k} > 10^{-50}$.

Note: We use $\tau_{\rm nl}$ to fix the adaptative time step (see adaptative_time_step).

init_temporal_scheeme!(::Euler_step, Run, dt)

Initialization of the temporal scheeme.

• First argument is the time-stepping method: Euler_step, AB_Euler_step, AB2_RK2_step, RK2_step, ETD2_step, RK4_step, ETD4_step.
• Run: run structure
• dt: time step

NewtonRaphson(f, x0; ϵ=1e-9, tol=1e-12, maxIter = 1000)

Find the zero of the function f using the Newton-Raphson algorithm.

• x0: starting point
• ϵ: spacing used to estimate derivative
• tol: tolerence
• maxIter: maximal number of iterations

Adapted from here.

NewtonRaphson_with_derivative(f, fp, x0, tol=1e-12, maxIter = 1000)

Same as NewtonRaphson, but using the explicit derivative of f, which is the argument fp, instead of a numerical estimate.

cardano(p::Float64, q::Float64)

Assuming that p > 0 and q < 0, it returns the positive real root of $x^3 +$ p $x +$ q.