src/test/nonlinear.c

    Non-linear geostrophic adjustment

    For a circular vortex defined by a tangential velocity u_\theta(r), the radial height/pressure profile is a solution of \displaystyle h'_0(r) = \frac{u_\theta}{g}\left(2\Omega + \frac{u_\theta}{r}\right) with \Omega the angular velocity. For this test case we take \displaystyle u_\theta(r) = (r < 0.4)\epsilon(1 + \cos(\pi(r - 0.2)/0.2))/2 The control parameters are the Froude number U/\sqrt{gH} and the Rossby number Ro=U/\Omega L. We set the Froude number to 0.1 and consider Rossby numbers 0.1 and \infty (no rotation). In the case without rotation the errors reflect only the accuracy of the momentum advection terms. With rotation, the errors also depend on the accuracy of the discretisation of the geostrophic balance.

    The figures below illustrate the evolution of the errors on free surface height/pressure for the different solvers, with and without rotation.

    set xlabel 'Time'
    set ylabel 'Maximum relative error'
    set logscale y
    plot 'log' index 'F0 = 0' u 1:2 w l t 'C grid', \
         '../nonlinear-ml/log' index 'F0 = 0' u 1:2 w l t 'multilayer'
    Evolution of the maximum relative error on free-surface height. Ro = \infty. (script)

    Evolution of the maximum relative error on free-surface height. Ro = \infty. (script)

    plot 'log' index 'F0 = 0.1' u 1:2 w l t 'C grid', \
         '../nonlinear-ml/log' index 'F0 = 0.1' u 1:2 w l t 'multilayer'
    Evolution of the maximum relative error on free-surface height. Ro = 0.1. (script)

    Evolution of the maximum relative error on free-surface height. Ro = 0.1. (script)

    The total energy is exactly conserved for the C-grid scheme and very well conserved for the multilayer scheme.

    set ylabel 'Normalised total energy'
    unset logscale
    set key bottom left
    plot 'log' index 'F0 = 0.1' u 1:3 w l t 'C grid (Ro = 0.1)', \
         '../nonlinear-ml/log' index 'F0 = 0.1' u 1:3 w l	     \
           t 'multilayer (Ro = 0.1)',			     \
         'log' index 'F0 = 0' u 1:3 w l t 'C grid (Ro = infty)', \
         '../nonlinear-ml/log' index 'F0 = 0' u 1:3 w l t 'multilayer (Ro = infty)'
    Evolution of the total energy (script)

    Evolution of the total energy (script)

    See also

    #include <gsl/gsl_integration.h>
    #pragma autolink -lgsl -lgslcblas
    #include "grid/multigrid.h"
    #if ML
    double F0 = 0.;
    # define F0() F0
    # include "layered/hydro.h"
    # include "layered/implicit.h"
    #else
    # include "atmosphere.h"
    #endif
    
    int main()
    {
      // coordinates of lower-left corner
      origin (-0.5, -0.5);
      // number of grid points
      init_grid (64);
      // size of the box
      size (1.);
      // acceleration of gravity
      G = 1.;
      // Coriolis parameter
      F0 = 0.1;
      // CFL number: the C-grid model is unstable for larger CFL
      CFL = 0.25;
      for (F0 = 0.; F0 <= 0.1; F0 += 0.1) {
        fprintf (stderr, "# F0 = %g\n", F0);
        run();
        fprintf (stderr, "\n\n");
      }
    }
    
    /* ---------------- Initial conditions ------------------- */
    
    #define H0 1.
    #define FROUDE 0.1
    
    double vtheta (double r) {
      return FROUDE*(r < 0.4)*(1. + cos((r - 0.2)/0.2*M_PI))/2.;
    }
    
    double h0p (double r, void * p) {
      double vt = vtheta(r);
      return vt*(F0 + vt/r)/G;
    }
    
    double h0 (double r) {
      gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);
      double result, error;
      gsl_function F;
      F.function = &h0p;
      gsl_integration_qags (&F, 0, r, 0, 1e-7, 1000, w, &result, &error);
      gsl_integration_workspace_free (w);
      return result;
    }
    
    scalar h1[];
    
    event init (i = 0)
    {
    #if ML
      CFL_H = 1.;
    #endif
      foreach()
        h1[] = h[] = (H0 + h0(sqrt (x*x + y*y)));
    #if ML
      foreach() {
        double r = sqrt (x*x + y*y), vt = vtheta(r);
        u.x[] = - vt*y/r;
        u.y[] =   vt*x/r;
      }
    #else
      foreach_face(x)
        u.x[] = - vtheta(sqrt (x*x + y*y))*y/sqrt (x*x + y*y);
      foreach_face(y)
        u.y[] =   vtheta(sqrt (x*x + y*y))*x/sqrt (x*x + y*y);
    #endif
    }
    
    /* ------------------ Output helper functions --------------- */
    
    scalar e[];
    
    double error()
    {
      double max = 0.;
      foreach(reduction(max:max)) {
        e[] = fabs (h1[]  - h[]);
        if (e[] > max) max = e[];
      }
      return max/(normf(h1).max - H0);
    }
    
    double energy()
    {
      double se = 0.;
      foreach(reduction(+:se)) {
        double ke = (sq(u.x[] + u.x[1,0]) + sq(u.y[] + u.y[0,1]))/8.;
        se += (h[]*ke + G*h[]*(h[]/2. + zb[]))*sq(Delta);
      }
      return se;
    }
    
    event logfile (i += 10; t <= 5.)
    {
      static double e0 = 0.;
      if (i == 0)
        e0 = energy();
      fprintf (stderr, "%g %g %.12g\n", t, error(), energy()/e0);
    }
    
    event plots (t = end)
    {
      output_ppm (e, file = "e.png", spread = -1, n = 256);
    }