The function "solver" was closed with an 'end', but at least one other function definition was not. All functions in a script must be closed with an 'end'.

20 views (last 30 days)
%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

Answers (1)

Navya Singam
Navya Singam on 1 Dec 2021
Hi,
"end" statement can be used to terminate the declared function. "end" statement at the end of function is optional, but if a file contains multiple functions and one of the function is terminated with "end", in that case all the other functions should also be terminated with "end" statement.
The error states that only the "solver" function is terminated with end,
The function "solver" was closed with an 'end', but at least one other function definition was not. All functions in a script must be closed with an 'end'.
It can be resolved by terminating all the declared functions with "end" statement.
Refer to end documentation for more information on usage of "end" keyword in functions.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!