sandbox/M1EMN/Exemples/maree_bretagne.c

    Marées proche Atlantique, Manche et Mer du Nord

    A temptative to compute tides around France based on Tsunami example, with a given M2 wave (only Moon contribution) at the left of the domain in Atlantic ocean the first of August 2017 (time in TU+1). It needs 3 days to obtain the forced state.

    Solver setup

    The following headers specify that we use the Saint-Venant solver together with (dynamic) terrain reconstruction.

    #include "spherical.h"
    #include "saint-venant.h"
    #include "terrain.h"

    We then define a few useful macros and constants.

    #define ETAE     1e-2 // error on free surface elevation (1 cm)
    #define HMAXE    5e-2 // error on maximum free surface elevation (5 cm)
    double z0=0,tmax;

    When using a quadtree (i.e. adaptive) discretisation, we want to start with the coarsest grid, otherwise we directly refine to the maximum level. Note that 1 << n is C for 2^n.

    int main()
    {
    #define MAXLEVEL 10   // as The domain is 25 degrees hence 25*30=750 points ($2^9=512$). 9 / 6
    #define MINLEVEL 7
    
      #if QUADTREE
      //  grid points to start with
      N = 1 << MINLEVEL;
      #else // Cartesian
      //   grid points
      N = 1 << MAXLEVEL;
      #endif

    Here we setup the domain geometry. For the moment Basilisk only supports square domains. For this example we have to use degrees as horizontal units because that is what the topographic database uses (eventually coordinate mappings will give more flexibility). We set the size of the box L0 and the coordinates of the lower-left corner (X0,Y0).

      Radius = 6371220.;
      // the domain is 25 degrees squared
      L0 = 25;
      // centered on 94,8 longitude,latitude
      X0 = -2 - L0/2.;
      Y0 = 48. - L0/2.;

    G is the acceleration of gravity required by the Saint-Venant solver. This is the only dimensional parameter. We rescale it so that time is in minutes, horizontal distances in degrees and vertical distances in metres. This is a trick to circumvent the current lack of coordinate mapping.

    units : meters in z but degrees en x,y, so 60 Miles 40075e3/360.= 111km which is 60 times one Nauticale Miles 1852m.

      // acceleration of gravity in degrees^2/min^2/m
      //G = 9.81*sq(mtd)*sq(60.);
      //mtd = 360./40000e3; 
      G=9.81*sq(60.);
      // time in minutes, one day 24*60s computation during half a week
      tmax=22*60*4;//(24*60)*7;

    We then call the run() method of the Saint-Venant solver to perform the integration.

      run();
    }

    We declare and allocate another scalar field which will be used to store the maximum wave elevation reached over time.

    scalar hmax[];

    Boundary conditions

    We set the normal velocity component on the left, right and bottom boundaries to a “radiation condition” with a reference sealevel of zero. see http://basilisk.fr/src/elevation.h#radiation-boundary-conditions

    The top boundary is always “dry” in this example so can be left alone. Note that the sign is important and needs to reflect the orientation of the boundary. At first we just put the period of M2 tide of the Moon (a complete computation needs the Sun AS2) The period of M2 is 12.4206h = 12h24min14s (745.236 minutes, 44714.2 s).

    the amplitude 2.00 and phase shift 320 of the velocity at the entrance of the domain are adjusted to fit the M2 tide as soon as possible after the first of August 2017 (this is a trick).

    u.n[right]  = + radiation(0);
    u.n[bottom] = - radiation(0);
    u.n[top]    = + radiation(0);
    u.n[left]   = - radiation(2.00*sin(2*pi*(t+320)/(745.236)));

    Adaptation

    Here we define an auxilliary function which we will use several times in what follows. Again we have two #if…#else branches selecting whether the simulation is being run on an (adaptive) quadtree or a (static) Cartesian grid.

    We want to adapt according to two criteria: an estimate of the error on the free surface position – to track the wave in time – and an estimate of the error on the maximum wave height hmax – to make sure that the final maximum wave height field is properly resolved.

    We first define a temporary field (in the automatic variable η) which we set to h+z_b but only for “wet” cells. If we used h+z_b everywhere (i.e. the default \eta provided by the Saint-Venant solver) we would also refine the dry topography, which is not useful.

    int adapt() {
      #if QUADTREE
      scalar eta[];
      foreach()
        eta[] = h[] > dry ? h[] + zb[] : z0;
      boundary ({eta});

    We can now use wavelet adaptation on the list of scalars {η,hmax} with thresholds {ETAE,HMAXE}. The compiler is not clever enough yet and needs to be told explicitly that this is a list of doubles, hence the (double[]) type casting.

    The function then returns the number of cells refined.

      astats s = adapt_wavelet ({eta, hmax}, (double[]){ETAE,HMAXE}, MAXLEVEL, MINLEVEL);
      // fprintf (stderr, "# refined %d cells, coarsened %d cells\n", s.nf, s.nc);
      return s.nf;
      #else // Cartesian
      return 0;
      #endif
    }

    Initial conditions

    We first specify the terrain database to use to reconstruct the topography z_b. This KDT database needs to be built beforehand. See the xyz2kdt manual for explanations on how to do this.

    The horizontal grid spacing is 2-minutes of latitude and longitude (1 minute of latitude = 1.853 km at the Equator). The vertical precision is 1 meter. The domain is 25 degrees hence 25*30=750 points (2^9=512).

    The next line tells the Saint-Venant solver to conserve water surface elevation rather than volume when adapting the mesh.

    event init (i = 0)
    {
      terrain (zb, "/home/basilisk/terrain/etopo2", NULL);
      conserve_elevation();

    The initial still water surface is at z=0 so that the water depth h is…

      foreach(){ 
        h[] = max(0.0, 0.0 - zb[]);
      }
      boundary ({h});
    }

    At each timestep

    We use a simple implicit scheme to implement quadratic bottom friction i.e. \displaystyle \frac{\partial \mathbf{u}}{\partial t} = - C_f|\mathbf{u}|\frac{\mathbf{u}}{h} with C_f=10^{-4}.

    event friction (i++) {
      double Cf=1e-4;
      foreach() {
        double a = h[] < dry ? HUGE : 1. + Cf*dt*norm(u)/(h[]);
        foreach_dimension()
          u.x[] /= a;

    That is also where we update hmax.

        if (h[] > dry && h[] + zb[] > hmax[])
          hmax[] = h[] + zb[];
      }
      boundary ({hmax, u});
    }

    Add now we compute the Coriolis force \displaystyle \frac{\partial \mathbf{u}}{\partial t} = -2 \omega \sin \lambda \times \mathbf{u} with \lambda around 48° (this is y ), and sideral period 23h56min 04s = 1436.07 min.

    The splited derivative is solved in a implicit way: \displaystyle \frac{\partial u}{\partial t} = 2 \omega \sin \lambda v, \;\; \;\; \frac{\partial v}{\partial t} = -2 \omega \sin \lambda u define \Omega = 2 \omega \sin \lambda, write Z=u+iv so \frac{\partial Z}{\partial t} = -i \Omega Z so \displaystyle Z^{n+1}= \frac{Z^{n}}{1+ i \Omega \Delta t} which gives \displaystyle u^{n+1} = \frac{u^n+(\Omega \Delta t) v^n}{1 + (\Omega \Delta t)^2)},\;\; v^{n+1} = \frac{v^n -(\Omega \Delta t) u^n }{1 + (\Omega \Delta t)^2)}

    event coriolis (i++) {
      double Omeg,Ro,u1,v1;
        foreach() {
          Omeg=2*sin(y*pi/180.)*(2*pi/1436.07);
          Ro=1 + sq(Omeg*dt);
          u1=(u.x[] + u.y[]*(Omeg*dt))/Ro; 
          v1=(u.y[] - u.x[]*(Omeg*dt))/Ro; 
          u.x[]=u1;
          u.y[]=v1;
        }
       boundary ({u.x,u.y}); 
    }

    Movies

    This is done every minute (t++). The static variable fp is NULL when the simulation starts and is kept between calls (that is what static means). The first time the event is called we set fp to a ppm2mpeg pipe. This will convert the stream of PPM images into an mpeg video using ffmpeg externally.

    We use the mask option of output_ppm() to mask out the dry topography. Any part of the image for which m[] is negative (i.e. for which etam[] < zb[]) will be masked out.

    event movies (t+=10) {
      static FILE * fp = popen ("ppm2mpeg > eta.mpg", "w");
      scalar m[], etam[];
      foreach() {
        etam[] = eta[]*(h[] > dry);
        m[] = etam[] - zb[];
      }
      boundary ({m, etam});

    save only the last period in a film (play the film in loop!), and save one image to display

     // if(t>(tmax - 745.236))
        output_ppm (etam, fp, mask = m,min = -3, max = 3, n = 1024,  linear = true);
     
        if(t==2500) {
      output_ppm (etam, file = "eta.png", mask = m,min = -5, max = 5, n = 512,  linear = true);
      output_ppm (hmax, file = "hmax.png", mask = m,min = -5, max = 5, n = 512,  linear = true);    
    }

    We also use the box option to only output a subset of the domain (defined by the lower-left, upper-right coordinates).

      static FILE * fp2 = popen ("ppm2mpeg > eta-zoom.mpg", "w");
      output_ppm (etam, fp2                 , mask = m, min = -3 , max = 3, n = 1024, linear = true,box = {{-4.8,48},{2.,51.3}});   
      output_ppm (etam, file ="eta-zoom.mp4", mask = m, min = -3 , max = 3, n = 1024, linear = true,box = {{-4.8,48},{2.,51.3}});  
    
      //if((t>=tmax-2*745.236)&&(t<tmax))
      //{output_ppm (etam, file ="eta-SM.mp4",   min = -3 , max = 3, n = 1024, linear = true,box = {{-4.14,48.46},{1.23,49.8}}); }  
        
      if(t>0)
         output_ppm (etam, file ="eta-SM.mp4",   min = -5 , max = 5, n = 512, linear = true, box = {{-4.14,48.5},{-1.5,50}}); 
     
        if(t==1000) {
        output_ppm (etam, file = "eta-zoom.png", mask = m,min = -3, max = 3, n = 256,  linear = true,box = {{-4.8,48},{2.,51.3}});
    }

    And repeat the operation for the level of refinement…

      static FILE * fp1 = popen ("ppm2mpeg > level.mpg", "w");
      scalar l = etam;
      foreach()
        l[] = level;
      output_ppm (l, fp1, min = MINLEVEL, max = MAXLEVEL, n = 512);
          if(t==1000) {
      output_ppm (l, file = "level.png", mask = m, min = MINLEVEL, max = MAXLEVEL, n = 256, linear = true);
    }
    
     foreach()
        l[] = level;
        boundary ({l}); 
        output_ppm (l, n = 1024,  min = 0.,   max= 12.,
          box = {{0,-1},{10,2.5}},
          file = "level.mp4")  ;

    …and for the process id for parallel runs.

     // fprintf (stderr, " heur= %g, i.e. jour=%g \n",t/60,t/60/24.);
     
    }

    Tide gauges

    We define a list of file names, locations and descriptions and use the output_gauges() function to output timeseries (for each timestep) of \eta for each location.

    Gauge gauges[] = {
      // file   lon      lat         description
      {"manche.txt", 0,  50, "#dans la manche"},
      {"saintmalo.txt", -2.140962, 49, "#SM"},
      {"rotterdam.txt", 4.095950 , 52.092836, "#Rott"},
      {"brest.txt", -4.806819, 48.214563, "#Brst"},
      {"houat.txt", -2.708353,47.406863,"#Houat"},
      {"large.txt", -14,  48, "#large"},
      {NULL}
    };
    
    event gauges1 (t += 15; t <= tmax) output_gauges (gauges, {eta});

    Adaptivity

    And finally we apply our adapt() function at every timestep.

    event do_adapt (i++) adapt();

    Run

    compilation

    qcc  -g -O2 marees_bretagne.c -o marees_bretagne ~/basilisk/basilisk/src/kdt/kdt.o -lm 
    ./marees_bretagne
     make maree_bretagne.tst;make maree_bretagne/plots;make maree_bretagne.c.html

    Results images

    After completion this will give the following animation

    Animation of the wave elevation. Dark blue is -2 metres and less. Dark red is +2 metres and more.

    Animation of the wave elevation. Dark blue is -2 metres and less. Dark red is +2 metres and more.

    maximum elevation of tide, note that Saint Malo and Bristol bay are indeed high tide places. Note the green spot corresponing to the amphidromic points north of France. Max h

    a zoom in “La Manche”

    Animation of the wave elevation. Dark blue is -2 metres and less. Dark red is +2 metres and more.

    Animation of the wave elevation. Dark blue is -2 metres and less. Dark red is +2 metres and more.

    Animation of the wave elevation. Dark blue is -2 metres and less. Dark red is +2 metres and more.

    Animation of the wave elevation Saint Malo

    adaptation level

    Animation of the level of refinement. Dark blue is 7 (128 pt 25*60./2**7=11 M/cel) and dark red is 9 (512 pt, precision 3M/cel)

    Animation of the level of refinement. Dark blue is 7 (128 pt 25*60./2**7=11 M/cel) and dark red is 9 (512 pt, precision 3M/cel)

    Results tides

    Value of the tide in Brest from SHOM the 1st of August 2017 (time in TU+1) \displaystyle Z(t)= Z_0 + AS_2 \cos(GS_2- WS_2 t) + AM_2 \cos(GS_2 - WM_2 t) with Moon parameters amplitude AM_2=2.16, phase GS_2=5642.320805440571, period WM_2=0.5058680493365497/60 which is 12h25min14s (low tide to hight tide 6h12m37s) 44714.1643s

    and sun parameters WS_2=0.5235987755982988/60 which is 12 hours, (low tide to hight tide 6h) AS_2 =0.755m the amplitude GS_2 = 3.1764992386296798 the phase shift

    and mean height Z_0=4.45m

    From Table des Marées des Grands Ports du Monde, SHOM 1984 \displaystyle \frac{\text{Z0}}{100} + \frac{\text{AM2} \cos \left(\frac{1}{180} \pi \left(-\text{GM2}+2 \left(\text{h0}+\left(\frac{\text{heure}}{24}+13270\right) \text{hp0}\right)-2 \left(\left(\frac{\text{heure}}{24}+13270\right) \text{sp0}+\text{s0}\right)+30 \text{heure}\right)\right)+\text{AS2} \cos \left(\frac{1}{180} \pi (30 \text{heure}-\text{GS2})\right)}{1000} with the parameters for Brest \displaystyle \{\text{Z0}\to 445,\text{AM2}\to 2160,\text{GM2}\to 142,\text{AS2}\to 755,\text{GS2}\to 182\} and astronomical parameters \displaystyle \{\text{s0}\to 78.16,\text{sp0}\to 13.1764,\text{h0}\to 279.82,\text{hp0}\to 0.985647\}

    here we just plot the tide from SHOM

    set xlabel "t in min"
    set ylabel "h in m"
    Z0=4.45
    AM2=2.16
    AS2 =0.755
    h1(x)=AS2*cos(3.1764992386296798 - 0.5235987755982988*x/60) 
    h2(x)=AM2*cos(5642.320805440571 - 0.5058680493365497*x/60)
    p[0:24*60]Z0+h1(x)+h2(x) t'maree Brest SHOM 2 modes' w l,  Z0+h2(x) t'onde M2 uniq.'
    maree a Brest deux premiers modes le 1er Aout 2017 pendant un jour (script)

    maree a Brest deux premiers modes le 1er Aout 2017 pendant un jour (script)

    Plot of the tide in Brest and height at the boundary condition, the phase shift and amplitude of the left BC have been adjusted to fit the tide in Brest

    set xlabel "t in min"
    set ylabel "h in m" 
    p[0:]'large.txt' w lp,'brest.txt' w lp,   h2(x) t'onde M2 uniq.'
    au large (condition d’entree) (script)

    au large (condition d’entree) (script)

    Once, the phase shift and amplitude of the left BC have been adjusted in radiation BC, we check that the computed M2 tide in Saint Malo is not so far from the SHOM the 1st of August 2017 \displaystyle Z(t)= Z_0 + AM_2 \cos(GS_2 - WM_2 t) Z_0=6.71 m with Moon parameters amplitude AM_2=3.68, phase is GS_2=5643.47, period WM_2=0.5058680493365497/60
    (we remove the Z_0 in the plot as we start by a mean constant level which takes into account the initial depth in the tiopography)

    set xlabel "t in min"
    set ylabel "h in m" 
    sm2(x)= 3.68*cos(5643.47 - 0.505868*x/60)
    p[0:][:7]'brest.txt' w lp, h2(x) t' Brest onde M2 uniq.','saintmalo.txt' w lp,sm2(x) t' SM onde M2 uniq.'
    Brest calculé et Brest SHOM et Saint Malo Calculé et Saint Malo SHOM (script)

    Brest calculé et Brest SHOM et Saint Malo Calculé et Saint Malo SHOM (script)

    set xlabel "t in min"
    set ylabel "h in m"
    sm2(x)= 3.68*cos(5643.47 - 0.505868*x/60)
    p[64*60:88*60][:7]'brest.txt' w lp, h2(x) t' Brest onde M2 uniq.','saintmalo.txt' w lp,sm2(x) t' SM onde M2 uniq.'
    Brest calculé et Brest SHOM et Saint Malo Calculé et Saint Malo SHOM (script)

    Brest calculé et Brest SHOM et Saint Malo Calculé et Saint Malo SHOM (script)

    Plot of tide in several harbours; note amplification in Saint Malo

    set key left
    set xlabel "t in min"
    set ylabel "h in m"
    p[0:][:4]'saintmalo.txt' w lp,'rotterdam.txt' w l,'manche.txt' w l,h2(x) t'M2','brest.txt' w lp,'houat.txt' w l 
    (script)

    (script)

    Links

    bibliography

    notes

    bug de formule en local sous macosx

     
     make maree_bretagne/plots;make maree_bretagne.c.html
     
     sed -i -e 's/\\)//g' maree_bretagne.c.html;sed -i -e 's/\\(//g' maree_bretagne.c.html
     sed -i -e 's/\\\[//g' maree_bretagne.c.html
     sed -i -e 's/\\\]//g' maree_bretagne.c.html;

    reste le problème des films à résoudre

    [Animation] of the wave elevation. LIEN TOUT BETE
    LIEN TOUT BETE