---
tags: HW4.2
---
# ***Computational Hydraulics***
## **Group 01:** N88071053 洪豪男 N88091011 黃彥鈞
## **HW4.2:** Please check the speed of the hydraulic jump by the m-code “code_20191106a_SWE_jump.m”. Does it satisfy the shock speed condition mentioned in Unit 02?

%
% Shallow Water Equations (conservative form)
%                   (NOC scheme + TVD minmod)
%    dh/dt + d(hu)/dx = 0   (NOC scheme, Burger equation)
% d(hu)/dt + d/dx (hu^2 + 0.5 g h^2 cos(theta)) = g h sin(theta) - drag
%
%  drag = g * h * cos(theta) * tan(delta)
%
%  u(:,1) :  h
%  u(:,2) :  hu
%  u(:,3) :  u = hu/h
%
clear all;
%
xx = -1:0.1:14;   % computational domain
nx = length(xx);  % number of grid (cell)
%
dx = xx(2) - xx(1); % size of cell
CFL = 0.4;          % CFL number
%
TVD_on = 2;  % 0: without  1: Minmod  2: superbee TVD recon
%
incAng = 0.0/180*pi;    % Inclination angel
dlta = 0.0/180*pi;      % friction angle
grav   = 9.8;           % gravitational accelration
Hcut = 1.0e-8;          % for computing the velocity (1.0e-8)
%
%dt = 0.01;
%dtdx = dt/dx;
%
%   Initial condition
% 
u0 = zeros(nx,3);
h0 = 0.5;
Cpos = 10.0;
R = 1.5;
Vel_0 = 1.0;
%
%   Initial condition (parabolic cap)
%
for ii = 1:length(xx)
    if  xx(ii) >= Cpos
        u0(ii,1) =  h0;  % 0.5
        u0(ii,3) =  0.1; % u
    else
        u0(ii,1) =  0.1*h0;
        u0(ii,3) =  1.0; % u
    end
    u0(ii,1) = u0(ii,1) + 0.0;
    u0(ii,2) =  u0(ii,1)*u0(ii,3); % hu
end
%
% Plot the initial condition
% 
figure(1);
axes('position',[0.1 0.1 0.8 0.4]);
plot(xx,u0(:,1)); ylabel('h');
axis([-1 14 -0.1 0.8]);
axes('position',[0.1 0.56 0.8 0.4]);
plot(xx,u0(:,3)); ylabel('u');
axis([-1 14 -0.1 1.8]);
%
% pause;
%
u = u0;
uone = u0;  % at j+1/2 (for first step)
unew = u0;  % at j
du = zeros(size(u));
dupls = zeros(size(u));
dumns = zeros(size(u));
%
F_p = zeros(nx,2);   % flux at RHS
F_m = zeros(nx,2);   % flux at LHS
%
TotalV = sum(u(:,1))*dx;
TimeTotal = 0.0;
%
figure(2);
for nstep = 1:50
    % compute the time step
    hmax = max(u(:,1));
    umax = max(abs(u(:,3)));
    a_max = max(umax+sqrt(grav*hmax*cos(incAng)));  % maximum wave speed
    dt = CFL * dx/a_max;
%  First step
    for ii = 1:nx-1
        F_p(ii,:) = fkt_flux(u(ii+1,:),grav,incAng);
        F_m(ii,:) = fkt_flux(u(ii  ,:),grav,incAng);
    end
    F_p(nx,:) = F_p(nx-1,:);
    F_m(nx,:) = F_m(nx-1,:);
%   Compute the value of du (first step)
    for mm = 1:3
    for ii = 2:nx-1
        dupls = u(ii+1,mm) - u(ii  ,mm);
        dumns = u(ii  ,mm) - u(ii-1,mm);
        if TVD_on == 1
            du(ii,mm) = fkt_Minmod(dupls,dumns); % Minmod limiter
        elseif TVD_on == 2
            du(ii,mm) = fkt_Superbee2(dupls,dumns); % Superbee limiter
        else
            du(ii,mm) = 0.0;
        end
    end
    end
%
    for ii = 2:nx-1
        htmp  = 0.5*(u(ii,1)+0.25*du(ii,1)+u(ii+1,1)-0.25*du(ii+1,1));
        hutmp = 0.5*(u(ii,2)+0.25*du(ii,2)+u(ii+1,2)-0.25*du(ii+1,2));
        if htmp < 0.0
            htmp = 0.0; hutmp = 0.0;
        end
        %utmp  = hutmp/(htmp+Hcut);
        utmp = 0.5*(u(ii,3)+0.25*du(ii,3)+u(ii+1,3)-0.25*du(ii+1,3));
        
        uone(ii,1) = htmp - dt/dx* (F_p(ii,1)-F_m(ii,1));
        uone(ii,2) = hutmp ...
            - dt/dx* (F_p(ii,2)-F_m(ii,2)) + grav*htmp*sin(incAng) ...
            - dt*fkt_drag2(htmp,utmp,grav,incAng,dlta);
    end
    uone(:,3) = uone(:,2)./(uone(:,1)+Hcut);  % compute the velocity
%   B.C.
    uone(1:2,:)  = u0(1:2,:);
    uone(nx-1:nx,:) = u0(nx-1:nx,:);
%
    TimeTotal = TimeTotal + dt
%
%   Second step
%
    for ii = 2:nx
        F_p(ii,:) = fkt_flux(uone(ii  ,:),grav,incAng);
        F_m(ii,:) = fkt_flux(uone(ii-1,:),grav,incAng);
    end
    F_p(1,:) = F_p(2,:);
    F_m(1,:) = F_m(2,:);
%   Compute the value of du  (2nd step)
    for mm = 1:3
    for ii = 2:nx-1
        dupls = uone(ii+1,mm) - uone(ii  ,mm);
        dumns = uone(ii  ,mm) - uone(ii-1,mm);
        if TVD_on == 1
            du(ii,mm) = fkt_Minmod(dupls,dumns); % Minmod limiter
        elseif TVD_on == 2
            du(ii,mm) = fkt_Superbee2(dupls,dumns); % Superbee limiter
        else
            du(ii,mm) = 0.0;
        end
    end
    end
%
    for ii = 2:nx
        htmp  = 0.5*(uone(ii,1)-0.25*du(ii,1)+uone(ii-1,1)+0.25*du(ii-1,1));
        hutmp = 0.5*(uone(ii,2)-0.25*du(ii,2)+uone(ii-1,2)+0.25*du(ii-1,2));
        if htmp < 0.0
            htmp = 0.0; hutmp = 0.0;
        end
        %utmp  = hutmp/(htmp+Hcut);
        utmp = 0.5*(uone(ii-1,3)+0.25*du(ii-1,3)+uone(ii,3)-0.25*du(ii,3));
        
        unew(ii,1) = htmp - dt/dx* (F_p(ii,1)-F_m(ii,1));
        unew(ii,2) = hutmp ...
            - dt/dx* (F_p(ii,2)-F_m(ii,2)) + grav*htmp*sin(incAng)...
            - dt*fkt_drag2(htmp,utmp,grav,incAng,dlta);
    end
    unew(:,3) = unew(:,2)./(unew(:,1)+Hcut);  % compute the velocity
%   B.C.
    unew(1:2,:)  = u0(1:2,:);
    unew(nx-1:nx,:) = u0(nx-1:nx,:);
%
    TimeTotal = TimeTotal + dt
    TotalV = sum(unew(:,1))*dx
%
    u = unew;
%{
    axes('position',[0.1 0.1 0.8 0.4]);
    plot(xx,u0(:,1),'--r',xx,u(:,1),'-b'); hold on;
    ylabel('h');
    stepTXT = ['step No. = ',num2str(nstep)];
    legend('initial',stepTXT);
    timeTXT = ['time = ',num2str(TimeTotal),' s'];
    text(6,0.5,timeTXT);
    axis([-1 14 -0.1 0.6]);
    axes('position',[0.1 0.56 0.8 0.4]);
    plot(xx,u(:,2)); ylabel('hu');
    axis([-1 14 -5. 15.]); hold off;
    pause(0.5);
    %}
    plot(xx,u0(:,3),'--r',xx,u(:,3),'-b'); 
    xlabel('x');
    ylabel('u');
    axis([0 14 -2 3]);
    timeTXT = ['time = ',num2str(TimeTotal),' s'];
    text(5,2.4,timeTXT);
    stepTXT = ['step No. = ',num2str(nstep)];
    text(5,2.2,stepTXT);
    pause(0.5);
end
    
%


* We see that a curve of discontinuity is a shock curve satisfies the Rankine-Hugoniot jump condition and the entropy condition for that solution u.