# sandbox/fcor_test.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 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187  /* Western Boundary Current intensification This is using Stommel's western boundary current as a test case Stommel, Henry (April 1948). "The Westward Intensification of Wind-Driven Ocean Currents". Transactions, American Geophysical Union 29 (2): 202–206. doi:10.1029/tr029i002p00202 See http://ido.at.fcen.uba.ar/index_archivos/Stommel_1948.pdf A sinusoidal wind with wind stress in the x direction of $1.7 \times 10^{-4} \cos(\pi y / L0)$ m/s blows in the x direction over a square pool of L0 = 1000 km side length with constant depth 2 km. There is linear friction with coefficient 2.5e-3. A beta-plane approximation is made with f=-9.5e-5+1.7e-11*y Which corresponds to the origin being at 41 degrees South.  */ #include "grid/cartesian.h" #include "storm_surge.h" #define MAXLEVEL 8 #define MINLEVEL 4 #define ETAE 1e-8 double t1=8*604800.; /** Here we can set standard parameters in Basilisk */ int main() { /** Here we setup the domain geometry. For the moment Basilisk only supports square domains. This case uses metres east and north. We set the size of the box L0 and the origin to be the lower-left corner (X0,Y0). In this case we are assuming a square 'pool' of length 1000 km. */ // the domain is L0 = 1000000.; /** G is the acceleration of gravity required by the Saint-Venant solver. This is the only dimensional parameter.. */ G = 9.81; #if QUADTREE // 32^2 grid points to start with init_grid( 1 << MINLEVEL ); #else // Cartesian // 1024^2 grid points init_grid( 1 << MAXLEVEL ); #endif run(); } /** ## 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](http://en.wikipedia.org/wiki/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[] : 0; boundary ({eta}); astats s = adapt_wavelet ({eta}, (double[]){ETAE}, MAXLEVEL, MINLEVEL); fprintf (stderr, "# refined %d cells, coarsened %d cells\n", s.nf, s.nc); return s.nf; #else // Cartesian return 0; #endif } event initiate(i=0) { foreach() { fcor[]=-9.5e-5+1.7e-11*y; zb[] = -2000.; h[] = max(0., - zb[]); ts.x[] = -1.7e-4*cos(pi*y/L0); ts.y[] = 0.; } boundary ({h,zb}); } /** We want the simulation to stop when we are close to steady state. To do this we store the h field of the previous timestep in an auxilliary variable hn. */ scalar hn[]; event init_hn (i = 0) { foreach() { hn[] = h[]; } } /** We output running statistics to the standard error. */ event logfile (i+=10; t eta.mpg", "w"); scalar m[], etam[]; foreach() { etam[] = eta[]*(h[] > dry); m[] = etam[] - zb[]; } boundary ({m, etam}); output_ppm (etam, fp, min = -0.01, max = 0.01 , n = 512, linear = true); #if QUADTREE static FILE * fp1 = NULL; if (!fp1) fp1 = popen ("ppm2mpeg > level.mpg", "w"); scalar l = etam; foreach() l[] = level; output_ppm (l, fp1, min = MINLEVEL, max = MAXLEVEL, n = 512); #endif } /** ## Adaptivity We apply our adapt() function at every timestep. */ event do_adapt (i++) adapt();