# Green-Naghdi soliton

The Green-Naghdi system of equations admits solitary wave solutions of the form $\eta \left(\psi \right)=a{h}_{0}{\mathrm{\text{sech}}}^{2}\left(\psi \frac{\sqrt{3a{h}_{0}}}{2{h}_{0}\sqrt{{h}_{0}\left(1+a\right)}}\right)$ $u\left(\psi \right)=\frac{c\eta }{{h}_{0}+\eta }$ with $\psi =x-ct$ and the soliton velocity ${c}^{2}=g\left(1+a\right){h}_{0}$.

``````#include "grid/multigrid1D.h"
#include "green-naghdi.h"``````

The domain is 700 metres long, the acceleration of gravity is 10 m/s2. We need to set the dispersion parameter ${\alpha }_{d}$ to one. We compute the solution in one dimension for a number of grid points varying between 128 and 1024.

``````int main()
{
X0 = -200.;
L0 = 700.;
G = 10.;
alpha_d = 1.;
for (N = 128; N <= 1024; N *= 2)
run();
}``````

We follow Le Métayer et al, 2010 (section 6.1.2) for the values of ${h}_{0}$ and $a$.

``````double h0 = 10, a = 0.21;

double sech2 (double x) {
double a = 2./(exp(x) + exp(-x));
return a*a;
}

double soliton (double x, double t)
{
double c = sqrt(G*(1. + a)*h0), ψ = x - c*t;
double k = sqrt(3.*a*h0)/(2.*h0*sqrt(h0*(1. + a)));
return a*h0*sech2 (k*ψ);
}

event init (i = 0)
{
double c = sqrt(G*(1. + a)*h0);
foreach() {
double η = soliton (x, t);
h[] = h0 + η;
u.x[] = c*η/(h0 + η);
}
}``````

We output the profiles and reference solution at regular intervals.

``````event output (t = {0,7.3,14.6,21.9,29.2}) {
if (N == 256) {
foreach()
fprintf (stdout, "%g %g %g %g\n", x, h[] - h0, u.x[], soliton (x, t));
fprintf (stdout, "\n");
}
}``````

We compute the error between the theoretical and numerical solutions at $t=29.2$.

``````event error (t = end) {
scalar e[];
foreach()
e[] = h[] - h0 - soliton (x, t);
fprintf (stderr, "%d %g\n", N, normf(e).max/(a*h0));
}``````

The method has a second-order rate of convergence as expected.