# Title: Tsunami runup onto a complex three-dimensional beach
#
# Description:
#
# This example is a classical validation test case for tsunami
# models. It was proposed at the \htmladdnormallinkfoot{"Third
# international workshop on long-wave runup
# models"}{http://isec.nacse.org/workshop/2004\_cornell/bmark2.html}. It
# is based on experimental data obtained in a wave tank in order to
# understand the extreme runups observed near the village of Monai
# during the 1993 Okushiri tsunami.
#
# The animation in Figure \ref{monai} gives a general idea of the
# geometry and time evolution of the modelled tsunami. The bathymetry
# data and channel geometry matches that used in the experimental wave
# tank. The water surface is forced on the open boundary with the
# experimentally-imposed waveform (outside the field of view on the
# right-hand-side in the animation). The initial dryout as well as
# extreme runup in the narrow central valley are clearly visible in
# the animation, as well as wave reflections from the boundaries and
# dimples in the water surface caused by underwater vortices.
#
# \begin{figure}[htbp]
# \caption{\label{monai}Animation of the water surface (light blue)
# and bathymetry (white) for the Monai tsunami.}
# \begin{center}
# \htmladdnormallinkfoot{\includegraphics[width=0.8\hsize]{monai.eps}}{monai.mpg}
# \end{center}
# \end{figure}
#
# Another view of the process together with the adaptive mesh used to
# resolve the flow is given in Figure \ref{topview}. This figure can
# be compared with that of LeVeque and George (Figure 4 of
# \cite{leveque2006}).
#
# \begin{figure}[htbp]
# \caption{\label{topview}Evolution of the free surface elevation
# (left column) and corresponding adaptive mesh (right
# column). Contour lines of the topography are represented. The areas
# in white in the left column are dry. The wet areas are coloured
# according to free surface elevation relative to the unperturbed
# water surface. A ``jet'' colour scale is used with a maximum value
# (dark red) of $+5$ cm and a minimum value (dark blue) of $-5$ cm.}
# \begin{center}
# \begin{tabular}{cc}
# \includegraphics[width=0.5\hsize]{fig4-10.eps} &
# \includegraphics[width=0.5\hsize]{mesh-10.eps} \\
# \multicolumn{2}{c}{$t = 10$ s} \\
# \includegraphics[width=0.5\hsize]{fig4-12.eps} &
# \includegraphics[width=0.5\hsize]{mesh-12.eps} \\
# \multicolumn{2}{c}{$t = 12$ s} \\
# \includegraphics[width=0.5\hsize]{fig4-14.eps} &
# \includegraphics[width=0.5\hsize]{mesh-14.eps} \\
# \multicolumn{2}{c}{$t = 14$ s} \\
# \includegraphics[width=0.5\hsize]{fig4-16.eps} &
# \includegraphics[width=0.5\hsize]{mesh-16.eps} \\
# \multicolumn{2}{c}{$t = 16$ s} \\
# \includegraphics[width=0.5\hsize]{fig4-18.eps} &
# \includegraphics[width=0.5\hsize]{mesh-18.eps} \\
# \multicolumn{2}{c}{$t = 18$ s} \\
# \includegraphics[width=0.5\hsize]{fig4-20.eps} &
# \includegraphics[width=0.5\hsize]{mesh-20.eps} \\
# \multicolumn{2}{c}{$t = 20$ s}
# \end{tabular}
# \end{center}
# \end{figure}
#
# Experimental data includes free-surface elevation timeseries at
# several locations (see Figure \ref{topview}, $t = 18$ s). Figures
# \ref{p5}, \ref{p7} and \ref{p9} give comparisons of the experimental
# and numerical timeseries.
#
# \begin{figure}[htbp]
# \caption{\label{p5}Time-series of free-surface elevation measured
# and calculated at the location of probe 5.}
# \begin{center}
# \includegraphics[width=0.8\hsize]{p5.eps}
# \end{center}
# \end{figure}
#
# \begin{figure}[htbp]
# \caption{\label{p7}Time-series of free-surface elevation measured
# and calculated at the location of probe 7.}
# \begin{center}
# \includegraphics[width=0.8\hsize]{p7.eps}
# \end{center}
# \end{figure}
#
# \begin{figure}[htbp]
# \caption{\label{p9}Time-series of free-surface elevation measured
# and calculated at the location of probe 9.}
# \begin{center}
# \includegraphics[width=0.8\hsize]{p9.eps}
# \end{center}
# \end{figure}
#
# Finally, the animation in Figure \ref{comparison} gives a comparison
# between the overhead video of the experiment and the corresponding
# view of the simulation.
#
# \begin{figure}[htbp]
# \caption{\label{comparison}Comparison between the video of the
# experiment (left) and the simulation results (right).}
# \begin{center}
# \htmladdnormallinkfoot{\includegraphics[width=0.8\hsize]{comparison.eps}}{comparison.mp4}
# \end{center}
# \end{figure}
#
# More details on this simulation and the method used is given in
# \cite{popinet2011}.
#
# Author: St\'ephane Popinet
# Command: sh monai.sh
# Required files: monai.sh 3D.gfv leveque.gfv mesh.gfv overhead.gfv probe.gfv
# Version: 110107
# Running time: 50 minutes
# Generated files: monai.mpg monai.eps fig4-10.eps mesh-10.eps fig4-12.eps mesh-12.eps fig4-14.eps mesh-14.eps fig4-16.eps mesh-16.eps fig4-18.eps mesh-18.eps fig4-20.eps mesh-20.eps p5.eps p7.eps p9.eps comparison.eps comparison.mp4

# below this depth the flow is considered "dry"
Define DRY 1e-4

# use the GfsRiver Saint-Venant solver
# shift the origin of the reference box to (0,0)
2 1 GfsRiver GfsBox GfsGEdge { x = 0.5 y = 0.5 } {
    # the domain is 3.402 m X 6.804 m
    # units for time are seconds
    PhysicalParams { L = 3.402 g = 9.81 }
    Refine 6

    # maintain the Zb1 variable using the GTS surface
    VariableFunction Zb1 bathy.gts

    # the initial water level is at z = 0, so the depth P is...
    Init {} {
        P = MAX (0., -Zb1)
    }

    # use a Sweby limiter rather than the default minmod which is too
    # dissipative
    AdvectionParams { gradient = gfs_center_sweby_gradient }

    # adapt down to 9 levels based on the slope of the (wet)
    # free-surface and with a tolerance of 1 mm
    AdaptGradient { istep = 1 } {
        cmax = 1e-3
        cfactor = 2
        maxlevel = 9
        minlevel = 6
    } (P < DRY ? 0. : P + Zb)

    # at each timestep
    Init { istep = 1 } {
	# Add a "shelf" to simulate the wall on the right-hand-side boundary
	Zb = (x > 5.448 ? 0.13535 : Zb1)
        # read in the experimental timeseries in variable 'input'
	input = input.cgd
        # implicit quadratic bottom friction with coefficient 1e-3
	U = (P > DRY ? U/(1. + dt*1e-3*Velocity/P) : 0.)
	V = (P > DRY ? V/(1. + dt*1e-3*Velocity/P) : 0.)
        P = (P > DRY ? P : 0.)
    }

    Time { end = 22.5 }
    OutputTime { istep = 10 } stderr
    OutputTiming { start = end } stderr
    OutputSimulation { step = 1 } stdout
    OutputSimulation { step = 1 } sim-%g.gfs
    EventScript { step = 1 } { gzip -f sim-*.gfs }

    # output data at probe locations
    OutputLocation { istep = 1 } input 1e-3 1.7 0
    OutputLocation { istep = 1 } p5 4.521 1.196 0
    OutputLocation { istep = 1 } p7 4.521 1.696 0
    OutputLocation { istep = 1 } p9 4.521 2.196 0

    # kinetic energy
    OutputScalarSum { istep = 1 } ke { v = (P > 0. ? Velocity2*P : 0.) }
    # free-surface elevation
    OutputScalarStats { istep = 1 } p { v = (Zb > 0. ? P : P + Zb) }
    OutputScalarNorm { istep = 1 } u { v = Velocity }
    OutputTime { istep = 1 } balance
    OutputBalance { istep = 1 } balance

    # generate movies
    GModule gfsview
    OutputView { start = 9 step = 0.0416 } { ppm2mpeg -s 640x480 > monai.mpg } { 
	width = 1280 height = 960
    } 3D.gfv
    OutputView { start = 14.63 end = 19.5 step = 0.033333333 } { 
	ppm2mpeg -s 400x600 > overhead.mpg 
    } { 
	width = 800 height = 1200
    } overhead.gfv
} {
    dry = DRY
}
GfsBox {
    # use 'subcritical' boundary condition to impose the experimental
    # water level on the left boundary
    left = Boundary { BcSubcritical U (input - Zb) }
    top = Boundary
    bottom = Boundary
}
GfsBox {
    top = Boundary
    bottom = Boundary
}
1 2 right