4 views (last 30 days)

Show older comments

This code plots an idealized cross section of a delta (geologic feature). It works when the shoreline (SH) transgression (lines 71-74) is commented out, however this doesn't allow the shoreline to transgress (which doesn't happen in real life). The lines that allow for shoreline transgression work only when line 55 is commented (Z).

Any thoughts on how I can remedy this error message? My original idea is to have an initial delta present at the start of the model run, but I am not sure how to write that into the code. Perhaps I need to adjust where the initial shoreline is to do that?

Thank you!

%% Input parameter values

Rab=0.3; %Sediment input

A=0.2; %sea level amplitude

B=1; %sea level frequency

%% setting up the domain

nx=1000; %length of delta

p = 2; %starting position

pos=p;

dx=.01;

%% Time vector

dt=.00005;

Tmax=10;

t=0:dt:Tmax;

nt=length(t);

%% Initialize time vectors

sc=zeros(1,nt); % cycle Shoreline position

rc=zeros(1,nt); % cycle Alluvial-basement position

sm=zeros(1,nt); % Fixed SL Shoreline position

rm=zeros(1,nt); % Fixed SL Alluvial-basement position

Zp=zeros(1,nt); % SL position

Vc=zeros(1,nt); % Volume

%% Initialize space vectors

x=zeros(1,nx);

%cycle values

Fc=zeros(1,nx); %sediment flux

Hc=zeros(1,nx); %Enthalpy

Hnewc=zeros(1,nx);

hc=zeros(1,nx); %height above current sea-level

% To define the dimension of F(flux)

for i=1:p

x(i)=(i-(p-0.5))*dx;

hc(i) = -x(i);

hm(i) = -x(i);

end

for i=p:nx

x(i)=(i-(p-0.5))*dx;

end

Z = 0;

Z0=1;

Fc(1) = Rab;

%Cycle time loop

for j=1:nt

%% Sea Level

Z = Z0 + A*sin(j*dt*B);

Zp(j) = Z;

sc(j)=(pos-p)*dx - Hc(pos)*dx/(-Z); %shoreline position

for i=p:nx-1

Fc(i) = min((hc(i) - hc(i+1))/dx, Hc(i)*dx/dt+Fc(i-1));

%% Exner equation

Hnewc(i) = Hc(i) + dt/dx*(Fc(i-1)-Fc(i));

%tracks SH regression

if Hnewc(pos)> Z

pos = pos + 1;

end

%tracks SH transregression

while Hc(pos-1) < Z

pos = pos-1;

end

Hc(i)=Hnewc(i);

hc(i) = max(Hc(i)-Z,0);

end

%% Counter

if mod(j,.5/dt) == 0

disp(['t=', num2str(j*dt)]);

end

%Checks if trajectories hit boundary of domain

if j>1000 %skips roughness at first few steps

if pos == nx-1 || rc(j)/dx >= p || abs(rc(j) - rc(j-1)) >= 100*dx

disp('end of domain reached')

break

end

end

end

plot(t,sc, 'b')

hold on

plot(t,Zp)

hold on

Walter Roberson
on 22 Feb 2020

Vimal Rathod
on 24 Feb 2020

Hi,

Using the following lines decreases the pos value to less than 1, Thus the MATLAB throws an error. From your given code, all the Hc values have 0 value and Z value is 1 (Thus the Error). Commenting the following loop makes sure the error doesn't occur and if you comment the line 55 (Where Z is initialized) Z will always have a value of zero (Thus the loop won't run) which is same as commenting the loop.

%tracks SH transregression

while Hc(pos-1) < Z

pos = pos-1;

end

You could put a certain boundary check that to make sure the pos value doesn't fall behind the boundary.

You could also debug your code, check the values by using breakpoints and verifying if your values are appropriate.

Refer to the following link to know more about debugging using breakpoints.

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

Start Hunting!