<< . .

. 43
( : 51)



. . >>

for (j=1;j<=lattice_size_y;j++)
{jdown[j]=j-1;
jup[j]=j+1;}
jdown[1]=lattice_size_y;
jup[lattice_size_y]=1;
}

Another common trick we will use involves time stepping. In the case of the linear diffusion
equation with constant diffusivity, the stability criterion is t < ( x)2 /(2D ). In nonlinear and/or
coupled equations in which the diffusivity changes through time as the solution evolves, the stability
criterion is sometimes dif¬cult to determine. Stability can be enforced by varying the time step to
limit the maximum change in the variable or variables that are being modeled. This approach is
crude but it works. The downside to this approach is that it requires some trial and error in order
to determine how big a step is allowed before the solution becomes unstable. In the example code
below, maxchange is set equal to the largest change allowed in any value of topo[i][j] during a single
time step (determined by trial and error). If the largest change exceeds that maximum value, the
values of topo revert to the values at the last time step topoold and the time step is reduced by a
factor of 2. If the largest change is smaller than one tenth of the maximum value then the time step
is increased by a factor of 1.2. In this way, the time step is increased and decreased dynamically in
order to maintain similar maximum changes in the function or functions we are modeling. We will
use this scheme in a number of example codes in this book.
max=0;
for (i=1;i<lattice_size_x;i++)
{deltah=timestep/deltax*(flux[i]-flux[iup[i]]);
topo[i]+=deltah;
if (fabs(deltah)>max) max=fabs(deltah);}
elapsedtime+=timestep;
A1.2 FTCS METHOD FOR 2D DIFFUSION EQUATION 225


if (max>maxchange)
{elapsedtime-=timestep;
timestep/=2.0;
for (i=1;i<=lattice_size_x;i++)
topo[i]=topoold[i];}
else
{if (max<0.1*maxchange) timestep*=1.2;}



A1.2 FTCS method for 2D diffusion equation

The sample code below solved the 2D diffusion equation using the FTCS method. The code accepts an
input grid of type float, size 300 — 300, and resolution 10 m/pixel from the local directory, calculates
the solution for t = duration, and then outputs the result to an ASCII ¬le in the local directory.


int lattice_size_x,lattice_size_y,*iup,*idown,*jup,*jdown;

main()
{ FILE *fp1,*fp2;
float delta,**topo,**topoold,D,duration,timestep;
int i,j,t,nsteps;

fp1=fopen("inputdem","r");
fp2=fopen("outputdem","w");
lattice_size_x=300;
lattice_size_y=300;
delta=10.0; /* m */
D=1.0; /* m^2/kyr */
duration=100.0; /* kyr */
setupgridneighbors();
topo=matrix(1,lattice_size_x,1,lattice_size_y);
topoold=matrix(1,lattice_size_x,1,lattice_size_y);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
{fscanf(fp1,"%f",&topo[i][j]);
topoold[i][j]=topo[i][j];}
timestep=0.5*delta*delta/(2*D);
nsteps=(int)(duration/timestep);
for (t=1;t<=nsteps;t++)
{for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
topo[i][j]+=timestep*D/(delta*delta)*(topoold[iup[i]][j]+topoold[idown[i]][j]
+topoold[i][jup[j]]+topoold[i][jdown[j]]-4*topoold[i][j]);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
topoold[i][j]=topo[i][j];}
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
226 CODES FOR SOLVING THE DIFFUSION EQUATION


fprintf(fp2,"%f\n",topo[i][j]);
fclose(fp1);
fclose(fp2);
}

As a concrete example, we consider the case of terrace diffusion described in Section 2.3.4. The
following code can be used to reproduce Figure 2.23b, given the input ¬le rillmask, available from
the CUP website. This code also makes use of an ˜˜ef¬ciency mask™™ that identi¬es the grid points and
their neighbors that are actively undergoing erosion. This speeds up the code early in the simulation
when erosion is localized close to the rills.


#define PI 3.141592653589793

main()
{ FILE *fp1,*fp2;
int length_x,length_y,delta,lattice_size_x,lattice_size_y,**rills,**efficiencymask,i,j;
float deltah,simulationlength,elapsedtime,timestep,dissectionrelief,kappa,depositslope;
float depositaspect,northwestdatum,**topo,**topoold,**initialtopo;

fp1=fopen("rillmask","r");
fp2=fopen("terracediffusion2d","w");
length_x=300; /* m */
length_y=300;
delta=1.0; /* m */
dissectionrelief=10.0;
kappa=1.0; /* m^2/kyr */
depositslope=0.06; /* m/m */
depositaspect=300*PI/180; /* surface drains SE; degrees converted to radians */
northwestdatum=0.0; /* elevations are relative to NE corner = 0 m */
lattice_size_x=length_x/delta;
lattice_size_y=length_y/delta;
topo=matrix(1,lattice_size_x,1,lattice_size_y);
topoold=matrix(1,lattice_size_x,1,lattice_size_y);
initialtopo=matrix(1,lattice_size_x,1,lattice_size_y);
rills=imatrix(1,lattice_size_x,1,lattice_size_y);
efficiencymask=imatrix(1,lattice_size_x,1,lattice_size_y);
timestep=0.5*delta*delta/(2*kappa);
simulationlength=50; /* kyr */
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
{fscanf(fp1,"%d",&rills[i][j]);
if (rills[i][j]>0) rills[i][j]=0; else rills[i][j]=1;
efficiencymask[i][j]=1;
if (i<lattice_size_x) if (rills[i+1][j]==1) efficiencymask[i][j]=0;
if (i>1) if (rills[i-1][j]==1) efficiencymask[i][j]=0;
if (j<lattice_size_y) if (rills[i][j+1]==1) efficiencymask[i][j]=0;
if (j>1) if (rills[i][j-1]==1) efficiencymask[i][j]=0;
topo[i][j]=northwestdatum-(i-1)*delta*depositslope*cos(depositaspect)+
(j-1)*delta*depositslope*sin(depositaspect);
A1.3 FTCS METHOD FOR 1D NONLINEAR DIFFUSION EQUATION 227


if (rills[i][j]==1) topo[i][j]-=dissectionrelief;
topoold[i][j]=topo[i][j];
initialtopo[i][j]=topo[i][j];}
elapsedtime=0.0;
while (elapsedtime<=simulationlength)
{printf("%f %f\n",elapsedtime,timestep);
for (j=2;j<=lattice_size_y-1;j++)
for (i=2;i<=lattice_size_x-1;i++)
if (efficiencymask[i][j]==0)
{if (rills[i][j]==0)
deltah=timestep*kappa/(delta*delta)*(topoold[i+1][j]+topoold[i-1][j]+
topoold[i][j+1]+topoold[i][j-1]-4*topoold[i][j]);
else
deltah=0;
topo[i][j]+=deltah;
if (fabs(topo[i][j]-initialtopo[i][j])>0.01*dissectionrelief)
{efficiencymask[i+1][j]=0;
efficiencymask[i-1][j]=0;
efficiencymask[i][j+1]=0;
efficiencymask[i][j-1]=0;}}
elapsedtime+=timestep;
i=1;
for (j=1;j<=lattice_size_y;j++)
topo[i][j]=topo[i+1][j]+depositslope*cos(depositaspect)*delta;
i=lattice_size_x;
for (j=1;j<=lattice_size_y;j++)
topo[i][j]=topo[i-1][j]-depositslope*cos(depositaspect)*delta;
j=1;
for (i=1;i<=lattice_size_x;i++)
topo[i][j]=topo[i][j+1]-depositslope*sin(depositaspect)*delta;
j=lattice_size_y;
for (i=1;i<=lattice_size_x;i++)
topo[i][j]=topo[i][j-1]+depositslope*sin(depositaspect)*delta;
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
topoold[i][j]=topo[i][j];}
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
fprintf(fp2,"%f\n",topo[i][j]);
fclose(fp1);
fclose(fp2);
}



A1.3 FTCS method for 1D nonlinear diffusion equation

The code presented below is used to solve for the evolution of a 1D hillslope according to the
nonlinear diffusion equation (described in Section 2.3.3) subject to instantaneous base level drop.
228 CODES FOR SOLVING THE DIFFUSION EQUATION


main()
{ FILE *fp1;
int *iup,*idown,printed,outputinterval,lattice_size_x,i;
float kappa,Sc,slope,den,*flux,duration,deltax,deltah,elapsedtime,max,timestep;
float *topo,*topoold,maxchange;

fp1=fopen("nonlinbasedrop","w");
lattice_size_x=102;
deltax=1.0; /* m */
kappa=1.0; /* m^2/kyr */
Sc=1.0; /* m/m */
iup=ivector(1,lattice_size_x);
idown=ivector(1,lattice_size_x);
topo=vector(1,lattice_size_x);
topoold=vector(1,lattice_size_x);
flux=vector(1,lattice_size_x);
duration=1000.0; /* kyr */
outputinterval=1000.0;
maxchange=0.001;
for (i=1;i<=lattice_size_x;i++)
{topo[i]=1.0;
if ((i==lattice_size_x)||(i==lattice_size_x-1)||(i==lattice_size_x-2)) topo[i]=0.0;
topoold[i]=topo[i];
flux[i]=0;
iup[i]=i+1;idown[i]=i-1;}
iup[lattice_size_x]=lattice_size_x;idown[1]=1;
timestep=0.0001*deltax*deltax;
elapsedtime=0.0;printed=0;
while (elapsedtime<duration)
{printed=elapsedtime/outputinterval;
for (i=1;i<=lattice_size_x;i++)
{slope=(topoold[idown[i]]-topoold[i])/deltax;
den=1-fabs(slope)*fabs(slope)/(Sc*Sc);
if (den<0.01) den=0.01;
flux[i]=kappa*slope/den;
if (i==1) flux[i]=0.0;
if (i==lattice_size_x) flux[i]=flux[idown[i]];
if (flux[i]<flux[idown[i]]) flux[i]=flux[idown[i]];}
max=0;
for (i=1;i<lattice_size_x;i++)
{deltah=timestep/deltax*(flux[i]-flux[iup[i]]);
topo[i]+=deltah;
if (fabs(deltah)>max) max=fabs(deltah);}
elapsedtime+=timestep;
if (max>maxchange)
{elapsedtime-=timestep;
timestep/=2.0;
for (i=1;i<lattice_size_x;i++)
topo[i]=topoold[i];}
A1.4 ADI METHOD FOR SOLVING THE 2D DIFFUSION EQUATION 229


else
{if (max<0.1*maxchange) timestep*=1.2;}
for (i=1;i<lattice_size_x;i++)
topoold[i]=topo[i];
if ((int)(elapsedtime/outputinterval)>printed)
{printf("%f %f\n",elapsedtime,timestep);
for (i=1;i<=lattice_size_x;i++)
fprintf(fp1,"%d %f\n",i,topo[i]);
printed=(int)(elapsedtime/outputinterval);}}
fclose(fp1);
}


A1.4 ADI method for solving the 2D diffusion equation

The following code accepts a 30 m/pixel US Geological Survey DEM covering Hanaupah Canyon from
¬le hanaupahtopo and uses the ADI technique to solve the diffusion equation on hillslopes. In the
hillslopediffusion subroutine, the diffusion equation is solved for each row and each column in
an alternating fashion using the tridag routine of Numerical Recipes, which inverts a tridiagonal
matrix. Hillslopes are de¬ned to be pixels with a contributing area less than 0.1 km2 . Channel pixels
act as ¬xed boundary conditions for the channel pixels. In this example, 10 Myr of diffusive evo-
lution is simulated with a diffusivity value of D = 10 m2 /kyr using time steps of 1 Myr. The ¬le
hanaupahtopo can be downloaded from the CUP website. The code below makes use of two functions:
fillinpitsandflats and mfdflowroute introduced in Appendix 2.

float D,delta,thresholdarea,*ax,*bx,*cx,*ux,*rx,*ay,*by,*cy,*uy,*ry;
float **topoold,*topovec,**topo,**flow,**flow1,**flow2,**flow3,**flow4,**flow5;
float **flow6,**flow7,**flow8;
int *topovecind,lattice_size_x,lattice_size_y,*iup,*idown,*jup,*jdown;

void hillslopediffusion()
{
int i,j,count;
float term1;

count=0;
while (count<5)
{count++;
for (i=1;i<=lattice_size_x;i++)
for (j=1;j<=lattice_size_y;j++)
topoold[i][j]=topo[i][j];
for (i=1;i<=lattice_size_x;i++)
{for (j=1;j<=lattice_size_y;j++)
{term1=D*timestep/(delta*delta);
if (flow[i][j]<thresholdarea)
{ay[j]=-term1;
cy[j]=-term1;
by[j]=4*term1+1;
ry[j]=term1*(topo[iup[i]][j]+topo[idown[i]][j])+topoold[i][j];}
230 CODES FOR SOLVING THE DIFFUSION EQUATION


else
{by[j]=1;
ay[j]=0;
cy[j]=0;
ry[j]=topoold[i][j];}
if (j==1)
{by[j]=1;
cy[j]=0;
ry[j]=topoold[i][j];}
if (j==lattice_size_y)
{by[j]=1;
ay[j]=0;
ry[j]=topoold[i][j];}}
tridag(ay,by,cy,ry,uy,lattice_size_y);
for (j=1;j<=lattice_size_y;j++)
topo[i][j]=uy[j];}
for (i=1;i<=lattice_size_x;i++)
for (j=1;j<=lattice_size_y;j++)
topoold[i][j]=topo[i][j];
for (j=1;j<=lattice_size_y;j++)
{for (i=1;i<=lattice_size_x;i++)
{term1=D*timestep/(delta*delta);
if (flow[i][j]<thresholdarea)
{ax[i]=-term1;
cx[i]=-term1;
bx[i]=4*term1+1;
rx[i]=term1*(topo[i][jup[j]]+topo[i][jdown[j]])+topoold[i][j];}
else
{bx[i]=1;
ax[i]=0;
cx[i]=0;
rx[i]=topoold[i][j];}
if (i==1)
{bx[i]=1;
cx[i]=0;
rx[i]=topoold[i][j];}
if (i==lattice_size_x)
{bx[i]=1;
ax[i]=0;
rx[i]=topoold[i][j];}}
tridag(ax,bx,cx,rx,ux,lattice_size_x);
for (i=1;i<=lattice_size_x;i++)
topo[i][j]=ux[i];}}
}

main()
{ FILE *fp1,*fp2;
float time,timestep,duration;
int i,j,t,dum;
A1.4 ADI METHOD FOR SOLVING THE 2D DIFFUSION EQUATION 231


fp1=fopen("hanaupahtopo","r");
fp2=fopen("hanaupahtopodiffuse","w");
lattice_size_x=640;
lattice_size_y=335;
delta=30.0; /* m */
D=10.0; /* m^2/kyr */
thresholdarea=0.1; /* km */
duration=10000.0; /* kyr */
timestep=1000.0; /* kyr */
setupgridneighbors();
/* grids needed for matrix inversion */
ax=vector(1,lattice_size_x);
ay=vector(1,lattice_size_y);
bx=vector(1,lattice_size_x);
by=vector(1,lattice_size_y);
cx=vector(1,lattice_size_x);
cy=vector(1,lattice_size_y);
ux=vector(1,lattice_size_x);
uy=vector(1,lattice_size_y);
rx=vector(1,lattice_size_x);
ry=vector(1,lattice_size_y);
/* other grids */
topo=matrix(1,lattice_size_x,1,lattice_size_y);
topovec=vector(1,lattice_size_x*lattice_size_y);
topovecind=ivector(1,lattice_size_x*lattice_size_y);
topoold=matrix(1,lattice_size_x,1,lattice_size_y);
flow=matrix(1,lattice_size_x,1,lattice_size_y);
flow1=matrix(1,lattice_size_x,1,lattice_size_y);
flow2=matrix(1,lattice_size_x,1,lattice_size_y);
flow3=matrix(1,lattice_size_x,1,lattice_size_y);
flow4=matrix(1,lattice_size_x,1,lattice_size_y);
flow5=matrix(1,lattice_size_x,1,lattice_size_y);
flow6=matrix(1,lattice_size_x,1,lattice_size_y);
flow7=matrix(1,lattice_size_x,1,lattice_size_y);
flow8=matrix(1,lattice_size_x,1,lattice_size_y);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
{flow[i][j]=(delta/1000)*(delta/1000); /* km */
fscanf(fp1,"%d",&dum);
topo[i][j]=dum;}
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
fillinpitsandflats(i,j);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
topovec[(j-1)*lattice_size_x+i]=topo[i][j];
indexx(lattice_size_x*lattice_size_y,topovec,topovecind);
t=lattice_size_x*lattice_size_y+1;
while (t>1)
232 CODES FOR SOLVING THE DIFFUSION EQUATION


{t--;
i=(topovecind[t])%lattice_size_x;
if (i==0) i=lattice_size_x;
j=(topovecind[t])/lattice_size_x+1;
if (i==lattice_size_x) j--;
mfdflowroute(i,j);}
time=0;
while (time<duration)
{time+=timestep;
hillslopediffusion();}
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
fprintf(fp2,"%f\n",topo[i][j]);
fclose(fp1);

<< . .

. 43
( : 51)



. . >>