1 function waterwave
  2 
  3  
  4 n = 64;                  % grid size
  5 g = 9.8;                 % gravitational constant
  6 dt = 0.01;               % hardwired timestep
  7 dx = 1.0;
  8 dy = 1.0;
  9 nplotstep = 8;           % plot interval
 10 ndrops = 5;              % maximum number of drops
 11 dropstep = 500;          % drop interval
 12 D = droplet(1.5,21);     % simulate a water drop
 13  
 14 % Initialize graphics
 15  
 16 [surfplot,top,start,stop] = initgraphics(n);
 17  
 18 % Outer loop, restarts.
 19  
 20 while get(stop,'value') == 0
 21    set(start,'value',0)
 22    
 23    H = ones(n+2,n+2);   U = zeros(n+2,n+2);  V = zeros(n+2,n+2);
 24    Hx = zeros(n+1,n+1); Ux = zeros(n+1,n+1); Vx = zeros(n+1,n+1);
 25    Hy = zeros(n+1,n+1); Uy = zeros(n+1,n+1); Vy = zeros(n+1,n+1);
 26    ndrop = ceil(rand*ndrops);
 27    nstep = 0;
 28  
 29    % Inner loop, time steps.
 30  
 31    while get(start,'value')==0 && get(stop,'value')==0
 32        nstep = nstep + 1;
 33  
 34        % Random water drops
 35        if mod(nstep,dropstep) == 0 && nstep <= ndrop*dropstep
 36            w = size(D,1);
 37            i = ceil(rand*(n-w))+(1:w);
 38            j = ceil(rand*(n-w))+(1:w);
 39            H(i,j) = H(i,j) + rand*D;
 40        end
 41      
 42        % Reflective boundary conditions
 43        H(:,1) = H(:,2);      U(:,1) = U(:,2);       V(:,1) = -V(:,2);
 44        H(:,n+2) = H(:,n+1);  U(:,n+2) = U(:,n+1);   V(:,n+2) = -V(:,n+1);
 45        H(1,:) = H(2,:);      U(1,:) = -U(2,:);      V(1,:) = V(2,:);
 46        H(n+2,:) = H(n+1,:);  U(n+2,:) = -U(n+1,:);  V(n+2,:) = V(n+1,:);
 47  
 48        % First half step
 49    
 50        % x direction
 51        i = 1:n+1;
 52        j = 1:n;
 53    
 54        % height
 55        Hx(i,j) = (H(i+1,j+1)+H(i,j+1))/2 - dt/(2*dx)*(U(i+1,j+1)-U(i,j+1));
 56    
 57        % x momentum
 58        Ux(i,j) = (U(i+1,j+1)+U(i,j+1))/2 -  ...
 59                  dt/(2*dx)*((U(i+1,j+1).^2./H(i+1,j+1) + g/2*H(i+1,j+1).^2) - ...
 60                             (U(i,j+1).^2./H(i,j+1) + g/2*H(i,j+1).^2));
 61    
 62        % y momentum
 63        Vx(i,j) = (V(i+1,j+1)+V(i,j+1))/2 - ...
 64                  dt/(2*dx)*((U(i+1,j+1).*V(i+1,j+1)./H(i+1,j+1)) - ...
 65                             (U(i,j+1).*V(i,j+1)./H(i,j+1)));
 66        
 67        % y direction
 68        i = 1:n;
 69        j = 1:n+1;
 70    
 71        % height
 72        Hy(i,j) = (H(i+1,j+1)+H(i+1,j))/2 - dt/(2*dy)*(V(i+1,j+1)-V(i+1,j));
 73    
 74        % x momentum
 75        Uy(i,j) = (U(i+1,j+1)+U(i+1,j))/2 - ...
 76                  dt/(2*dy)*((V(i+1,j+1).*U(i+1,j+1)./H(i+1,j+1)) - ...
 77                             (V(i+1,j).*U(i+1,j)./H(i+1,j)));
 78        % y momentum
 79        Vy(i,j) = (V(i+1,j+1)+V(i+1,j))/2 - ...
 80                  dt/(2*dy)*((V(i+1,j+1).^2./H(i+1,j+1) + g/2*H(i+1,j+1).^2) - ...
 81                             (V(i+1,j).^2./H(i+1,j) + g/2*H(i+1,j).^2));
 82    
 83        % Second half step
 84        i = 2:n+1;
 85        j = 2:n+1;
 86    
 87        % height
 88        H(i,j) = H(i,j) - (dt/dx)*(Ux(i,j-1)-Ux(i-1,j-1)) - ...
 89                          (dt/dy)*(Vy(i-1,j)-Vy(i-1,j-1));
 90        % x momentum
 91        U(i,j) = U(i,j) - (dt/dx)*((Ux(i,j-1).^2./Hx(i,j-1) + g/2*Hx(i,j-1).^2) - ...
 92                          (Ux(i-1,j-1).^2./Hx(i-1,j-1) + g/2*Hx(i-1,j-1).^2)) ...
 93                        - (dt/dy)*((Vy(i-1,j).*Uy(i-1,j)./Hy(i-1,j)) - ...
 94                          (Vy(i-1,j-1).*Uy(i-1,j-1)./Hy(i-1,j-1)));
 95        % y momentum
 96        V(i,j) = V(i,j) - (dt/dx)*((Ux(i,j-1).*Vx(i,j-1)./Hx(i,j-1)) - ...
 97                          (Ux(i-1,j-1).*Vx(i-1,j-1)./Hx(i-1,j-1))) ...
 98                        - (dt/dy)*((Vy(i-1,j).^2./Hy(i-1,j) + g/2*Hy(i-1,j).^2) - ...
 99                          (Vy(i-1,j-1).^2./Hy(i-1,j-1) + g/2*Hy(i-1,j-1).^2));
100    
101        % Update plot
102        if mod(nstep,nplotstep) == 0
103           C = abs(U(i,j)) + abs(V(i,j));  % Color shows momemtum
104           t = nstep*dt;
105           tv = norm(C,'fro');
106           set(surfplot,'zdata',H(i,j),'cdata',C);
107           set(top,'string',sprintf('t = %6.2f,  tv = %6.2f',t,tv))
108           drawnow
109        end
110       
111        if all(all(isnan(H))), break, end  % Unstable, restart
112    end
113 end
114 close(gcf)
115  
116 % ------------------------------------
117  
118 function D = droplet(height,width)
119 % DROPLET  2D Gaussian
120 % D = droplet(height,width)
121    [x,y] = ndgrid(-1:(2/(width-1)):1);
122    D = height*exp(-5*(x.^2+y.^2));
123  
124 % ------------------------------------
125  
126 function [surfplot,top,start,stop] = initgraphics(n);
127 % INITGRAPHICS  Initialize graphics for waterwave.
128 % [surfplot,top,start,stop] = initgraphics(n)
129 % returns handles to a surface plot, its title, and two uicontrol toggles.
130  
131    clf
132    shg
133    set(gcf,'numbertitle','off','name','Shallow_water')
134    x = (0:n-1)/(n-1);
135    surfplot = surf(x,x,ones(n,n),zeros(n,n));
136    grid off
137    axis([0 1 0 1 -1 3])
138    caxis([-1 1])
139    shading faceted
140    c = (1:64)'/64;
141    cyan = [0*c c c];
142    colormap(cyan)
143    top = title('Click start');
144    start = uicontrol('position',[20 20 80 20],'style','toggle','string','start');
145    stop = uicontrol('position',[120 20 80 20],'style','toggle','string','stop');