Simulation of a tsunami

This example also shows how to use a image to define the initial thickness.

You MUST use the double upwind scheme for tsunami simulation
(doubleUpwind = 1 – a simple scheme would smooth the results).

download the image here and add it in the same folder than the input file :

 Input file :

%A simple example to simulate a collapse in 2D
%basic instruction to write a VolcFlow file
nrow = 100;   %equations are in 2D. At least 3 lignes (or columns) should be defined, even for 1D simulation
ncol = 100; %number of column of the simulation
dx_horiz = 5;        %space step in x (in m)
dy_horiz = dx_horiz; %space step in y
dt=0.05;           %time step (in s)
tstep_adjust = 0;  %0 = the time step is fixed by the user, 1 = its varies in time using a CFL condition
doubleUpwind = 1; %To use the double-upwind scheme. Must be used for tsunami simulation
dtplot = 0.1;        %time step for plot
tmax =  6;        %duration of the simulation
g=9.81;                  %gravity
cohesion = 0;         %cohesion or yield strength (depending of the rheology used)
rho = 1000;              %density (kg/m3)
representation = ‘cla; hrep=h; hrep(h<=0)=NaN; surf(x,y,z+h,hrep-hini);  axis equal;pause(0.1);colormap(mycol);’; %function used for the visualisation of results (at each dtplot)
f_avi =  »;  %avi file done with the plots.  » = no movie is done
f_data =  »; %save various variables at each dtplot.  » = no data is saved
bound_cond =  »;
%vectors x and y – not needed by VolcFlow (for the visualisation of results only)
x = [dx_horiz/2:dx_horiz:ncol*dx_horiz-dx_horiz/2];
y = [dy_horiz/2:dy_horiz:nrow*dy_horiz-dy_horiz/2];
%definition of the topography
z = 50*sin(y’/100)*sin(x/50);
%definition of the initial thickness
h=30-z;    %the thickness is defined from a giver elevation
h(h<0)=0;  %to avoid the negative values of the water depth
hini = h;  %initial thickness
dh = (255-double(imread(‘volcflow.bmp’)))/10; %use an image for the initail wave
h=h+dh;
%initial velocity – 0 here
u_xx = z(1:nrow, 1:ncol-1)*0;
u_xy = z(1:nrow, 1:ncol-1)*0;
u_yy = z(1:nrow-1, 1:ncol)*0;
u_yx = z(1:nrow-1, 1:ncol)*0;
mycol = cool;