# Title: Geostrophic adjustment # # Description: # # We consider the geostrophic adjustment problem studied by # Dupont {\cite{dupont}} and Le Roux et al {\cite{leroux98}}. A Gaussian bump # $\eta ( x, y ) = \eta_0 e^{^{- \frac{x^2 + y^2}{R^2}}}$ # is initialised in a 1000$\times$1000 km, 1000 m deep square basin. A reduced # gravity $g = 0.01$ m/s$^2$ is used to approximate a 10 m-thick stratified surface # layer. On an $f$-plane the corresponding geostrophic velocities are given by # \begin{eqnarray*} # u ( x, y ) & = & \frac{2 g \eta_0 y}{f_0 R^2} e^{- \frac{x^2 + y^2}{R^2}},\\ # v ( x, y ) & = & - \frac{2 g \eta_0 x}{f_0 R^2} e^{- \frac{x^2 + y^2}{R^2}}, # \end{eqnarray*} # where $f_0$ is the Coriolis parameter. Following Dupont we set $f_0 = 1.0285 # \times 10^{- 4}$ s$^{- 1}$, $R = 100$ km, $\eta_0 = 599.5$ m which gives a # maximum geostrophic velocity of 0.5 m/s. # # In the context of the linearised shallow-water equations, the geostrophic # balance is an exact solution which should be preserved by the numerical # method. In practice, this would require an exact numerical balance between # terms computed very differently: the pressure gradient and the Coriolis terms # in the momentum equation. If this numerical balance is not exact, the # numerical solution will adjust toward numerical equilibrium through the # emission of gravity-wave noise which should not affect the stability of the # solution. This problem is thus a good test of both the overall accuracy of the # numerical scheme and its stability properties when dealing with # inertia--gravity waves. We note in particular that a standard A-grid # discretisation would develop a strong computational-mode instability in this # case. Also, as studied by Leroux et al, an inappropriate choice of # finite-element basis functions will result in growing gravity-wave noise. # # \begin{figure}[htbp] # \caption{\label{geo-error}Evolution of the maximum error on the surface height for the # geostrophic adjustment problem.} # \begin{center} # \includegraphics[width=\hsize]{geo_error.eps} # \end{center} # \end{figure} # # \begin{figure}[htbp] # \caption{\label{geo-error1}Evolution of the surface-height error field. (a) $t =$1.157 # days, (b) $t = 2.315$ days, (c) $t =$3.472 days, (d) $t =$4.630 days, (e) $t # =$17.361 days.} # \begin{center} # \begin{tabular}{ccccc} # \includegraphics[width=0.18\hsize]{error-100.eps} & # \includegraphics[width=0.18\hsize]{error-200.eps} & # \includegraphics[width=0.18\hsize]{error-300.eps} & # \includegraphics[width=0.18\hsize]{error-400.eps} & # \includegraphics[width=0.18\hsize]{error-1500.eps} \\ # (a) & (b) & (c) & (d) & (e) # \end{tabular} # \end{center} # \end{figure} # # Figures \ref{geo-error} and \ref{geo-error1} summarise the results obtained # when running the geostrophic adjustment problem on a $64 \times 64$ uniform # grid with a timestep $\Delta t = 1000$ s. The maximum error on the height # field (Figure \ref{geo-error}) is small even after 18 days. After a strong # initial transient corresponding to the emission of gravity waves, the error # reaches a minimum at day 3 and then slowly grows with time with modulations # due to the reflexions of the initial gravity waves on the domain boundaries. # As illustrated on figure \ref{geo-error1}, this growth is not due to any # instability of the solution but to the slow decrease of the maximum amplitude # of the Gaussian bump due to numerical energy dissipation. # # Author: St\'ephane Popinet # Command: sh geo.sh geo.gfs # Version: 100323 # Required files: geo.sh geo.gfv e.ref # Running time: 3 minutes # Generated files: geo_error.eps error-100.eps error-200.eps error-300.eps error-400.eps error-1500.eps # 1 0 GfsOcean GfsBox GfsGEdge {} { # dt = 1000 s Time { iend = 1580 dtmax = 0.10285 } Refine 6 # Lx = Ly = 1000 km # H0 = 1000 m # g = 0.01 m/s^2 PhysicalParams { g = 9.4534734306584e-4 } AdvectionParams { scheme = none } ApproxProjectionParams { tolerance = 1e-6 } Solid (z + 1.) Init {} { # e-folding radius = 100 km # umax = 0.5 m/s = sqrt(200)*exp(-1/2) U = (5.667583815e-4*200.*y*exp (-100.*(x*x + y*y))) V = (- 5.667583815e-4*200.*x*exp (-100.*(x*x + y*y))) P = (5.667583815e-4*exp (-100.*(x*x + y*y))) } # f0 = 1.0285e-4 s-1 SourceCoriolis 1 OutputErrorNorm { istep = 1 } { awk '{print $3/1.0285e-4/3600./24. " "$9*1000e6*1.0285e-4*1.0285e-4/0.01;}' > e } { v = P } { s = (5.667583815e-4*exp (-100.*(x*x + y*y))) unbiased = 1 v = E } OutputSimulation { istart = 100 iend = 500 istep = 100 } stdout EventScript { istart = 100 iend = 500 istep = 100 } { echo "Save error-$GfsIter.eps { format = EPS }" } OutputSimulation { istart = 1500 } stdout EventScript { istart = 1500 } { echo "Save error-$GfsIter.eps { format = EPS }"} EventScript { start = end } { cat <