src/examples/tides.c

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    95
    96
    97
    98
    99
    100
    101
    102
    103
    104
    105
    106
    107
    108
    109
    110
    111
    112
    113
    
    #include "saint-venant.h"
    #include "terrain.h"
    
    #define MAXLEVEL 9
    #define MINLEVEL 5
    
    // metres to degrees
    double mtd = 360./40075e3;
    
    int main()
    {
      // 512^2 grid points
      init_grid (1 << MAXLEVEL);
      // the domain is 5 degrees squared
      size (5.);
      // centered on 174,-40.8 longitude,latitude
      origin (174 - L0/2., -40.8 - L0/2.);
      /* rescale G so that time is in minutes, horizontal length scales in
         degrees and vertical length scales in metres */
      G = 9.81*sq(mtd)*sq(60.);
      run();
    }
    
    // M2 tidal frequency. The period is 12h25 (745 minutes).
    #define M2F (2.*M_PI/745.)
    
    // "radiation" boundary conditions on left,right,top,bottom
    u.n[left]   = - radiation (  sin(M2F*t));
    u.n[top]    =   radiation (  sin(M2F*t));
    u.n[right]  = + radiation (- sin(M2F*t));
    u.n[bottom] = - radiation (- sin(M2F*t));
    
    event init (i = 0)
    {
      // use several databases for topography zb
      terrain (zb, 
    	   "/home/popinet/terrain/niwa", 
    	   "/home/popinet/terrain/geo",
    	   "/home/popinet/terrain/etopo2", 
    	   NULL);
      // ensure that water level is conserved during refinement/coarsening
      // the default is to conserve volume
      conserve_elevation();
      // sealevel at z = 0
      foreach()
        h[] = max(0., - zb[]);
    }
    
    // every timestep
    event logfile (i++) {
      stats s = statsf (h);
      norm n = normf (u.x);
      if (i == 0)
        fprintf (stderr, "t i h.min h.max h.sum u.x.rms u.x.max dt\n");
      fprintf (stderr, "%g %d %g %g %g %g %g %g\n", t, i, s.min, s.max, s.sum, 
    	   n.rms, n.max, dt);
    
      foreach() {
        // quadratic bottom friction, coefficient 1e-4 (dimensionless)
        double a = h[] < dry ? HUGE : 1. + 1e-4*dt*norm(u)/(h[]*mtd);
        foreach_dimension()
          u.x[] /= a;
      }
      boundary ((scalar *){u});
    }
    
    // snapshots every hour
    event snapshots (t += 60; t <= 6000) {
      printf ("file: t-%g\n", t);
      output_field ({h, zb, eta, u}, stdout, n = N, linear = true);
    }
    
    // movies every 10 minutes
    event movies (t += 10) {
      scalar m[], etam[];
      foreach() {
        etam[] = eta[]*(h[] > dry);
        m[] = etam[] - zb[];
      }
      boundary ({m, etam});
      output_ppm (etam, mask = m, min = -1, max = 1, n = 512, linear = true,
    	      file = "eta.mp4");
    
      scalar vort = etam;
      vorticity (u, vort);
      output_ppm (vort, mask = m, // min = -1e-2, max = 1e-2, 
    	      n = 512, linear = true, file = "vort.mp4");
    
      scalar l = etam;
      foreach()
        l[] = level;
      output_ppm (l, min = MINLEVEL, max = MAXLEVEL, n = 512, file = "level.mp4");
    }
    
    // tide gauges
    
    Gauge gauges[] = {
      {"well", 96.88,  -12.13, "Wellington, New Zealand"},
      {NULL}
    };
    
    event gauges1 (i++) output_gauges (gauges, {eta});
    
    event adapt (i++) {
      scalar eta[];
      foreach()
        eta[] = h[] > dry ? h[] + zb[] : 0;
      boundary ({eta});
    
      double cmax = 1e-2;
      astats s = adapt_wavelet ({eta}, (double[]){cmax}, MAXLEVEL, MINLEVEL);
      fprintf (stderr, "# refined %d cells, coarsened %d cells\n", s.nf, s.nc);
    }