%FVM1D - FOR WET BED PROBLEMS ONLY
%by Brett F. Sanders (and pieces of code from Scott F. Bradford)
%
%This is a very simple 1D solver of the shallow-water equations
%that uses the Hancock predictor-corrector time-stepping scheme,
%the MUSCL method of slope limiting and variable reconstruction
%and Roe's approximate Riemann solver to compute fluxes.
%
%The code is kept as simple as possible to emphasize the basic
%flow of logic. To account for problems involving a dry bed,
%one needs to add a number of "if" statements to avoid
%division by zero. This makes the code pretty messy and
%therefore these lines have been omitted.
%
%Note that the code can either be run in a first order or
%second order accurate mode. The user can select from
%several limiters to see how these impact the solution.
%
%To run this program, copy all the .m files into a directory
%and run "fvm1d.m" by either typing "fvm1d" at the matlab
%command prompt or pushing the execute button in the matlab
%text editor.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%------------------Script edited by Adam M Gould USQ------------------%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%--------------------------Initial Setup------------------------------%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all
format long
global grav
global tol_h
%%%%%%%%%%%%%%%%%%%%%%%%%% Set up grid %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nc=200; %number of cells
nf=nc+1; %number of edges
L=1000; %length of channel
dx=L/nc; %length of cell
x=0:dx:L; %array of edge coordinates
xc=dx/2:dx:L-dx/2; %array of cell center coordinates
nm=0.2; %Manning n
%can include bed slope So
%%%%%%%%%%%% Set up time marching and output interval %%%%%%%%%%%%%%
dt=0.05; %time step (s)
nt=2000000; %number of time steps
ntplot=5000; %plot interval (number of time steps)
So=2/1000;
%%%%%%%%%%%%% Define bed elevation at faces z=f(x) %%%%%%%%%%%%%%%
%if using So t need this section
%z=[2:-.04:0 0:.04:2]
%z=zeros(size(x)); %flat bed - can enter own function here, z=f(x)
z=So*(L-x);
%zc=So*(L-xc);
%%%%%%%%%%%%%% Compute bed slope %%%%%%%%%%%%%%%%%%
%j = jth cell
for j=1:nc,
zc(j)=0.5d0*(z(j+1)+z(j)); %elevation of cell center
zmin(j)=min([z(j) z(j+1)]);%lowest point in the cell
zmax(j)=max([z(j) z(j+1)]);%highest point in the cell
hcrit(j)=0.5*(zmax(j)-zmin(j));
dz(j)=z(j+1)-z(j);%change in elavation at the cell
%So=dz/dx;
end
% for i=1:nc, %for cell 1 through to number of cells
% dz(i)=z(i+1)-z(i); %dimensions of length
% zc(i)=0.5d0*(z(i+1)+z(i)); %elevation of cell center
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%----------------------Set parameter values-----------------------------%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
grav=9.806;
tol_h=1.d-4; %height to define whether cell is wet or dry
%%%%%--------------------Set attributes of solver------------------%%%%%%
iorder=2; %1=first order scheme, 2=second order scheme
beta=6; %controls limiter used by model
%Notes on limiters
%beta=1 => Minmod
%beta=2 => Superbee
%beta=3 => Fromm scheme, predicts oscillations at
%sharp fronts
%beta=4 => Van Leer
%beta=5 => Van Albada
%beta=6 => Double Minmod
xo=L/20; %position along channel where wave starts
etalo=.2; %height on left hand side of the wave
etaro=-3; %height on right hand side of the wave
ulo=.1; %initial velocity left side of wave (U/S)
uro=.1; %initial velocity right side of wave (D/S)
for j=1:nc,
if (xc(j) < xo), %%If cell is before wave then
%eta(j)=etalo; %%height from datum = initial height
left
%%of wave
eta(j)=etalo+zc(j);
if (eta(j) > zmax(j)), %%if height from datum is higher
%%than the highest bed elevation
h(j)=eta(j)-zc(j); %%height of water = initial height
from
%%datum - height of bed
elseif (eta(j) > zmin(j));%if height from datum is smaller than
%%highest cell point but bigger than
%%the lowest cell point
h(j)=0.5*(eta(j)-zmin(j));%h = 1/2 height of water from datum
%%- bed elevation
else
h(j)=0; %else water = 0
end
if (h(j) > tol_h) %if water height is bigger than threshold
u(j)=ulo;%%velocity = velocity of left side of wave
else
u(j)=0;%if not velocity is zero
h(j)=0;
end
else %%if cell is after wave then
eta(j)=etaro; %%height from datum = initial height on
if (eta(j)-zc(j))<0;
eta(j)=zc(j);
h(j)=0;
end
%%right side of wave
if (eta(j) > zmax(j)), %%if height from datum is higher than
%%the highest bed elevation
h(j)=eta(j)-zc(j); %%height of water = initial height from
%%datum - height of bed
elseif ((eta(j)-zc(j)) > 0);%zmin(j));%if height
%from datum is smaller than
%%highest cell point but bigger than
%%the lowest cell point
h(j)=0.5*(eta(j)-zmin(j));%h = 1/2 height of water from datum
%- bed elevation
else
h(j)=0; %%if not water height is zero
end
if (h(j) > tol_h), %if water height is bigger than threshold
u(j)=ulo;%%velocity = velocity of left side of wave%I believe
%%this should be u(i)=uro;
else
u(j)=0;
h(j)=0;%if not water height is zero
end
end
if eta(j) < zmin(j), %if water height from datum is smaller than the
eta(j) = zmin(j); %%bed height there is no water at this point
%%and eta=bed height
end
end
%%%%%%%%%%%%%%%%%%%%%%%----Initialize arrays----%%%%%%%%%%%%%%%%%%%%%%%%%
uh=h.*u; %%height of water * velocity ????flow rate
matrix
deta=zeros(size(eta)); %%Initialize change in height from datum
du=zeros(size(u)); %%Initialize change in velocity
dh=zeros(size(h));
duh=zeros(size(uh));
t=0; %start time
for i=1:nt, %Begin time-marching loop
if (iorder == 2), %for second order accuracy only
deta = limiter(nc,beta,eta,h,tol_h);
du = limiter(nc,beta,u,h,tol_h);
%%we now have found the
%%difference in velocity and
%%difference in
%%height at the current moment in time
[etap, hp, up, Uhp, Sfnew]=predictor_d...
(nc,eta,h,uh,deta,du,dz,zc,dt,dx,nm,So,u);
%%we now have height and velocity
%%at t+(1/2dt)
[F, amax, Hplot] = fluxes_d...
(grav,nf,etap,up,Uhp,z,deta,du,duh,dz,zc); %Compute fluxes
%%these fluxes are the momentum and mass
%%fluxes at each cell interface
else
[F, amax, Hplot] = fluxes_d(grav,nf,eta,u,uh,z,deta,du,duh,dz,zc);
end
[uh, h, u] = corrector_d(nc,h,uh,F,S,dx,dt,So,nm,h,u);
eta=h+zc; %compute new free surface height
e=eta+0.5*u.^2/grav; %compute energy in units of length (head)
t=t+dt;
%%%%%%%%%%------------Check the courant number-----------%%%%%%%%%%%%
cr=amax*dt/dx; %compute Courant number
fprintf(1,'%g %d\n',i,cr)
if (cr > 1) %Stops program if Courant number exceeds one.
break
end
%%%%%%%%%%%%%---------------plotting loop---------------%%%%%%%%%%%%%
if (mod(i,ntplot) == 0),
hold on
figure('Position',get(0,'ScreenSize'))
subplot(3,1,1)
plot(xc,e,'r-',xc,h+zc,'b-',xc,zc,'k-')
ylabel('$\eta$ (m)','Interpreter','Latex')
axis([0 L -2 10])
legend('Energy','Free Surface','Bed')
subplot(3,1,2)
plot(xc,u,'b-')
xlabel('$x$ (m)','Interpreter','Latex')
ylabel('$u$ (m/s)','Interpreter','Latex')
axis([0 L -5 10])
legend('Velocity')
subplot(3,1,3)
plot(xc,uh,'b-')
xlabel('$x$ (m)','Interpreter','Latex')
ylabel('q (m^3/s)')
axis([0 L -2 10])
legend('Discharge')
%pause(0.005)
saveas(gcf, num2str(t), 'bmp')
close
end
end
%E.2 limiter
function df = limiter(nc,beta,f,h,tol_h) %%deta or du = (nc,beta,eta) or (nc,beta,u)
df(nc)=f(nc)-f(nc-1); %%Ghost cell boundary Extrapolation
df(1)=f(2)-f(1); %%Ghost cell boundary extrapolations
for i=2:nc-1, %%for rest of cells
df1=f(i+1)-f(i); %%forward difference
df2=f(i)-f(i-1); %%backwards difference
df(i)=limit(df1,df2,beta); %%Call limit function
end
%E.3 limit
%%Only really concerned with beta = 2
function f=limit(d1,d2,beta)
%
% 0 = first order, 1-2 = beta, 3= fromm, 4=vanleer, 5=vanalbada, 6=double
% minmod
if beta==0, %first order
f=0;
elseif beta >= 1 & beta <= 2%beta: minmod (beta=1) and superbee (beta=2)
if(d1*d2 < 0),
f=0;
else
s=sign(d1);
%%sign means>0=1 0=0 <0=-1
a=abs(d1);
%%absolute values
b=abs(d2);
%%absolute values
f=s*min(max([a b]),beta*min([a b]));
%%where f = df which is deta or du
%%df = sign(forward
%%diff)*the min of(max(forward diff, back diff),
%beta*min(foward diff,back diff))
end
elseif beta == 3 %Fromm
f=0.5*(d1+d2);
elseif beta == 4 %vanleer
if(d1*d2 <= 0),
f=0;
else
f=2*d1*d2/(d1+d2);
end
elseif beta == 5 %vanalbada
eps=1.e-20;
f=(d1*(d2*d2+eps)+d2*(d1*d1+eps))/(d1*d1+d2*d2+2*eps);
elseif beta == 6 %double minmod
if(d1*d2 < 0),
f=0;
else
s=sign(d1);
a=abs(d1);
b=abs(d2);
c=0.5*(a+b);
f=s*min([2*a 2*b c]);
end
end
%E.4 predictor
%%The prediction step is used to predict future velocity and future height
%the future is dt/2
%%This is know as the Hancock's method described on pg323 of
%%high resolution and non-oscillatory
%%note these equations are eq12(A) and eq13(v)without friction and bedslope
%%terms
function [etap, hp, up, Uhp,Sfnew]=predictor_d...
(nc,eta,h,Uh,deta,dU,dz,zc,dt,dx,nm,So,u);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%----------------------Variable Declaration---------------------------%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
b=2;%width of rectangular channel
global grav
global tol_h
dh=deta-dz;
a=1;
Sfnew=zeros(nc,1);
for a = 1:10;
for j=1:1:nc,
etap(j)=eta(j)-0.5*dt/dx*dU(j); %find the eta at t+1/2
hp(j)=etap(j)-zc(j);%find the height at t+1/2
%Sfold
if (h(j) > tol_h),%if water height is bigger than threshold
A(j)=h(j)*b;%Area for rectangle
P(j)=h(j)*2+b;%Perimeter for rectangle
R(j)=A(j)/P(j);% Hydraulic radius
Sfold(j)=nm*nm*u(j)*abs(u(j))/((R(j))^(4/3));%sf at current time
%than g*(manning n)^2 /h(at cell j)
else%if dry
Sfold(j)=0; %if not water height is zero
h(j)=0;
A(j)=0;
P(j)=0;
R(j)=0;
end
%find Velocity @ t+1/2
U_temp(j)=u(j)-((0.5*dt/dx)*((grav*dh(j))-(u(j)*dU(j))))...
+ (0.5*grav*dt*(So-(0.5*Sfold(j))-(0.5*Sfnew(j))));
if (hp(j) > tol_h),%if water height is bigger than threshold
Ap(j)=hp(j)*b;
Pp(j)=hp(j)*2+b;
Rp(j)=Ap(j)/Pp(j);
Sfnew(j)=(nm*nm*U_temp(j)*abs(U_temp(j)))/((Rp(j))^(4/3));
%than g*(manning n)^2 /h(at cell j)
else
Snew(j)=0; %if not water height is zero
hp(j)=0;
Ap(j)=0;
Pp(j)=0;
Rp(j)=0;
end
%refind velocity at t+1/2
up(j)=u(j)-(0.5*dt/dx)*(grav*dh(j)-u(j)*dU(j))...
+ 0.5*grav*dt*(So-(0.5*Sfold(j))-(0.5*Sfnew(j)));
Uhp(j)=up(j)*hp(j);%flow rate
if (hp(j) < tol_h),%if water height is smaller than threshold
up(j)=0;
Uhp(j)=0;
end
end
a=a+1;
end
%E.5 fluxes
function [F, amax0,Hplot] = fluxes_d(grav,nf,eta,u,uh,z,deta,du,duh,dz,zc)
global tol_h %dry bed threshhold
%I'm using a 2d solver in 1d, so I'm setting sn=0 and cn=1 for all cases
sn=0;
cn=1;
vl=0; %%%these are zero because 1D solver
vr=0; %%%These are zero because 1D solver
q=0.05; %flow rate
hr=eta(1)-0.5*deta(1)-(z(1)); %extrapolation for the ghost cell
if hr < tol_h, %if height is under tol_h set to 0.
hr=0;
end
F(1,1)=q; %mass flux
F(1,2)=q^2/hr+0.5*grav*hr^2; %momentum Flux
%%eq18 pg324
%Right Boundary : model as wall
hl=eta(nf-1)+0.5*deta(nf-1)-z(nf); %extrapolation for the ghost cell
if hl < tol_h,%if height is under tol_h set to 0.
hl=0;
end
ul=u(nf-1)+0.5*du(nf-1); %extrapolation for the ghost cell
[fdum, amax]=solver(hl,hl,ul,-ul,vl,vr,sn,cn); %solid wall boundary
F(nf,1)=fdum(1);
F(nf,2)=fdum(2);
% ul=((u(nf-1))+0.5*du(nf-1)); Free flowing boundary conditions
% F(nf,1)=hl*ul;
% F(nf,2)=ul^2*hl+0.5*grav*hl^2;
amax0=0;
%same as above just for every internal cell
for i=2:nf-1,
hl=eta(i-1)+0.5*deta(i-1)-(z(i));
if hl < tol_h, %if calculated h is below tol it is set to 0
hl=0;
end
if hl > tol_h,
ul=(uh(i-1)+0.5*duh(i-1))/hl; %find subcritical velocity for cell
%with atleast hl wet
if (ul*ul/(grav*hl)) > 1, %Froude number
ul=u(i-1)+0.5*du(i-1); %Super critical
end
else
hl = max([hl 0]);
ul = 0;
end
hr=eta(i)-0.5*deta(i)-(z(i));
if hr < tol_h,%if calculated h is below tol it is set to 0
hr=0;
end
if hr > tol_h,
ur=(uh(i)-0.5*duh(i))/hr;%find subcritical velocity for cell
%with atleast hr wet
if (ur*ur/(grav*hr)) > 1, %Froude number
ur=u(i)-0.5*du(i); %Super critical
end
else
hr = max([hr 0]);
ur = 0;
end
if (hl > tol_h || hr > tol_h), %If the water is above both call
solver
[fdum, amax]=solver(hl,hr,ul,ur,vl,vr,sn,cn);%find fluxes
F(i,1)=fdum(1);%Mass Flux
F(i,2)=fdum(2);%Momentum flux
amax0=max([amax0 amax]);
else
F(i,1)=0;%Dry Bed
F(i,2)=0;%Dry Bed
end
end
%E.6 solver
%%fluxes are calculated as described within pg290 and 291 of Finite
volume
%%model for shallow water flooding of arbitrary topography
function [F, amax]=solver(hl,hr,ul,ur,vl,vr,sn,cn)
global grav
%Compute Roe averages
%the majority of this code is used to set up variables to simplify the
end
%equations that documented below
duml=hl^0.5;
dumr=hr^0.5;
cl=(grav*hl)^0.5;
cr=(grav*hr)^0.5;
hhat=duml*dumr;
uhat=(duml*ul + dumr*ur)/(duml+dumr);
vhat=(duml*vl + dumr*vr)/(duml+dumr);
chat=(0.5*grav*(hl+hr))^0.5;
uperp=uhat*cn+vhat*sn;
dh=hr-hl;
du=ur-ul;
dv=vr-vl;
dupar=-du*sn+dv*cn;
duperp=du*cn+dv*sn;
dW=[0.5*(dh-hhat*duperp/chat); hhat*dupar; 0.5*(dh+hhat*duperp/chat)];
%%eq12
uperpl=ul*cn+vl*sn;
uperpr=ur*cn+vr*sn;
al1=uperpl-cl;
al3=uperpl+cl;
ar1=uperpr-cr;
ar3=uperpr+cr;
R=[1 0 1;
uhat-chat*cn -sn uhat+chat*cn;
vhat-chat*sn cn vhat+chat*sn];
%%eq11
da1=max([0 2*(ar1-al1)]);
da3=max([0 2*(ar3-al3)]);
a1=abs(uperp-chat);
a2=abs(uperp);
a3=abs(uperp+chat);
%Critical flow fix
if a1 < da1,
a1=0.5*(a1*a1/da1+da1);
end
if a3 < da3,
a3=0.5*(a3*a3/da3+da3);
end
%Compute interface flux
A=diag([a1 a2 a3]);
%%eq9
FL=[uperpl*hl; ul*uperpl*hl + 0.5*grav*hl*hl*cn;...
vl*uperpl*hl + 0.5*grav*hl*hl*sn];
%%eq7
FR=[uperpr*hr; ur*uperpr*hr + 0.5*grav*hr*hr*cn;...
vr*uperpr*hr + 0.5*grav*hr*hr*sn];
%%eq7
F=0.5*(FL + FR - R*A*dW);
%%eq8
amax=chat+abs(uperp);
%%just below eq9
%E.7 corrector
function [uhnew, hnew, unew]=corrector_d...
(nc,hold,uhold,F,dx,dt,So,nm,h,uold)
%%$----------------------Variable Declaration--------------------------%%
b=2;%channel width
global tol_h
global grav
%%$--------------------------------------------------------------------%%
a=0;
Sfnew=zeros(50);
for a = 1:10
for i=1:1:nc,
hnew(i)=hold(i)-dt/dx*(F(i+1,1)-F(i,1));%determine new height at t+1
Anew(i)=hnew(i)*b;%Area
if hnew(i)<tol_h%dry bed check
hnew(i)=0;
Anew(i)=0;
end
