## Answers - Visualization of Diffusion Model Results

In the April 1997 tutorial there were many tasks for you to perform. This document shows some script files and plots for solving these tasks. It is important to note that there are often many different ways to solve the problems - it is likely that your script files will use a different approach

To run these script files, you must save them as Matlab m-files. One way is to highlight the section of code, copy it to the clipboard, paste it into a new document and save it as a file (such as diff1b.m). The other approach would be to save the entire file diffvian.htm as a file and then edit out each script file. The only problem with this approach is in dealing with the less-than (<) and greater-than (>) symbols. Because of how HTML interprets these symbols, I have replaced all < and > symbols with the text strings &lt; and &gt;. You will have to edit the source code and replace the string &lt; with < (and &gt; with >)

Copy your program program diff1.m to diff1b.m and modify it as shown below.
```      % DIFF1B.M    1-Dimensional Diffusion Solution
%
% Tom Huber, March 1997
%
ncells = 25;       % Number of Cells in the Array
NTimes = 20;       % Number of Time steps to perform

FracLeft  = 1/10;   % Fraction diffusing to the left
FracSame  = 3/10;   % Fraction which remains in same cell
FracRight = 6/10;   % Fraction diffusing to the right

Values = zeros(ncells,1);  % Initialize the Array to have zero concentration
VMax = 100;                % Maximum Concentration
Values(ncells/4) = VMax;   % Initially set the concentration at the midpoint
AllValues = Values;

for i=1:NTimes  % Perform a total of NTimes time steps
NewValues = Values*FracSame;  % New values initially a fraction of original
for j=1:ncells  % For each cell except leftmost or rightmost
if j>1
NewValues(j-1) = NewValues(j-1)+Values(j)*FracLeft;  % Diffuse Left
end
if j<ncells
NewValues(j+1) = NewValues(j+1)+Values(j)*FracRight; % Diffuse Right
end
end % for j=1...

Values = NewValues;     % The updated values become the current values
AllValues = [AllValues Values];
end  % for i=1...
```

1. Make a script file valsurf.m that will plot and label the surface
```          % VALSURF.M   Plot the surface from the diffusion model
clf;                  % Clear the current figure
surf(AllValues)       % Plot the surface
colormap(jet)         % Change the colormap
view(0,90)            % Set the viewpoint to be above
colorbar              % Show a colorbar which has the values
xlabel('Time Step')   % Label the axes
ylabel('Position')
```
2. Extracting regions of the data set Values=AllValues(10,:) and Values=AllValues(:,5) the results are shown below
 Value = AllValues(10,:) This can be interpreted as the concentration at the 10th location. (Click on image for an enlarged view) Value = AllValues(:,5) This can be interpreted as the concentration after the 5th time step
3. Make an M-File plotloc.m to plot the results at a selected location
```      % plotloc.m   Plot Results at some location
%   Uses the array AllValues and plots the results
%
loc = input('Enter Cell Location to Display ');
Values = AllValues(loc,:);
plot(Values)
xlabel('Time Step')
ylabel('Concentration')
title(['Cell Number ' num2str(loc)])
[ValMax LocMax] = max(Values);
disp(['Max Value = ' num2str(ValMax) ' At Time Step ' num2str(LocMax)])```
4. When does the crest (maximum) occur for a position 5 cells downstream from the initial concentration? What about 15 cells downstream? For 5 cells downstream from maximum (at Location 10), the crest value is 22.6 and occurs at time step 7. For 15 cells downstream, we need to extend the time range in diff1b.m by setting NTimes=35. This results in a crest value of 11.4 at time step 27. Does doubling the initial concentration double the crest at these locations? Yes
5. Running the model with different left, right, and same fractions, we can view the results using plotloc2.m below.
 Comparing the original run (Model 1) with left, same, right fractions of (1/10, 3/10, 6/10) to a run (Model 2) with fractions (0.5/10, 3/10, 6.5/10) we see that the peak in the second model comes earlier and is higher.
```      % plotloc2.m   Plot Results at some location
%   Uses the array AllValues and plots the results
%
loc = input('Enter Cell Location to Display ');
clf
Value1 = AllValue1(loc,:);
plot(Value1,'g-')
Values = AllValues(loc,:);
hold on
plot(Values,'r-')
xlabel('Time Step')
ylabel('Concentration')
title(['Cell Number ' num2str(loc)])
legend('Model 1','Model 2')
[ValMax LocMax] = max(Value1);
disp(['Model 1  Max Value = ' num2str(ValMax)  ...
' At Time Step ' num2str(LocMax)])
[ValMax LocMax] = max(Values);
disp(['Model 2  Max Value = ' num2str(ValMax)  ...
'  At Time Step ' num2str(LocMax)])

```
6. To view the results as a movie, you can use the file Valmovie.m below
```      % Valmovie.m   Plot Results for successive times
[NCells NSteps] = size(AllValues);
VMax = max(max(AllValues));
clf
for i=1:NSteps
Values = AllValues(:,i);
plot(Values,'r-')
xlabel('Location')
ylabel('Concentration')
title(['Time Step ' num2str(i)])
axis([1 NCells 0 VMax])
drawnow
end
```
Another variation on Valmovie.m (which I call fstmovie.m uses Matlab's ability to do animation using HandleGraphics. In this program, instead of displaying the axes & titles each time, only the Y-Data values are changed on the graph. The resulting movie executes much faster than the previous version
```      % Fstmovie.m   Plot Results for successive times
%   Uses HandleGraphics for faster animation of the movie
[NCells NSteps] = size(AllValues);
VMax = max(max(AllValues));
TDelay = .1;   % This is the minimum amount of time between movies
clf;
plot ([],[]);
axis([1 npts 0 VMax]);
axis(AXIS);
x = 1:NCells;  % Initialize X Data for 1 to NCells
%     Draw the line and store its handle in LineData
LineData = line('color','y','erase','xor','xdata',x,'ydata',x);
xlabel('Position')
ylabel('Concentration')
%     Display Iteration number and store its handle in TextData
TextData = text(npts/10,0.9*VMax,'Iteration: 0','Color','r');
set(TextData,'erase','Background')
drawnow

for i=1:NSteps
Time0 = clock;  % Get the current time
Values = AllValues(:,i);
set(LineData,'ydata',Values); % Update the Y Values on the Graph
set(TextData,'string',['Iteration: ' num2str(i)]) % Update text
drawnow

while etime(clock,Time0)<TDelay  % Pauses at least .1 second
% disp('This pauses on a fast CPU')
end %   while ...
end   % for i...
```

1. The first task has a reservoir that "melts" at some rate (this is similar to the model melting.m which we did at the beginning of the Tutorial on Matlab Scripts.
 Comparing the runs shown here, one notes that a slower melting rate leads to a lower peak value, but also leads to a longer flooding time (the peak is wider)
```      % DIFF1MLT.M    1-Dimensional Diffusion Solution with a Melting Reservoir
%
% Tom Huber, March 1997
%
ncells = 25;       % Number of Cells in the Array
NTimes = 35;       % Number of Time steps to perform

FracLeft  = 1/10;   % Fraction diffusing to the left
FracSame  = 3/10;   % Fraction which remains in same cell
FracRight = 6/10;   % Fraction diffusing to the right

Values = zeros(ncells,1);  % Initialize the Array to have zero concentration
VMax = 100;                % Maximum Concentration
MeltFrac = 0.1;            % Fraction which "Melts" each day

VRemain = VMax;
Values(ncells/4) = MeltFrac*VRemain; % Initially let some fraction melt
VRemain = (1-MeltFrac)*VRemain;   % This is fraction remaining after step
AllValues = Values;

for i=1:NTimes  % Perform a total of NTimes time steps
% First we let the current water flow downstrem
NewValues = Values*FracSame;  % New values initially a fraction of original
for j=1:ncells  % For each cell except leftmost or rightmost
if j>1
NewValues(j-1) = NewValues(j-1)+Values(j)*FracLeft;  % Diffuse Left
end
if j<ncells
NewValues(j+1) = NewValues(j+1)+Values(j)*FracRight; % Diffuse Right
end
end % for j=1...
Values = NewValues;     % The updated values become the current values

%  Now we will add some new melted water to the river
Melted = MeltFrac*VRemain;   % This is how much melted
VRemain = VRemain - Melted;  % This is how much remains
Values(ncells/4) = Values(ncells/4) + Melted;  % Add the melted

AllValues = [AllValues Values];
end  % for i=1...```
2. The next task is to have a reservoir at every location. This is accomplished by making VRemain and Meleted vectors. In the example below, the reservoir of 100 is replaced by 5 reservoirs of 20 (the melting fraction is 0.4 in both cases).
 Comparing the runs shown here, one notes that an extended reservoir (5 cells of 20 each - solid curve) leads to a lower peak value, but also leads to a longer flooding time than a single reservoir
```      % DIFF1RES.M    1-Dimensional Diffusion Solution with a Reservoir
%               at each location
%
% Tom Huber, March 1997
%
ncells = 25;       % Number of Cells in the Array
NTimes = 35;       % Number of Time steps to perform

FracLeft  = 1/10;   % Fraction diffusing to the left
FracSame  = 3/10;   % Fraction which remains in same cell
FracRight = 6/10;   % Fraction diffusing to the right

MeltFrac = 0.4;            % Fraction which "Melts" each day

VRemain = zeros(ncells,1);      % Initialize all reservoirs to zero

%  Now load 5 reservoirs (locations 4-8)
%  with 20 instead of 1 reservoir with 100
VRemain(4:8) = 20*ones(5,1);

Values = MeltFrac*VRemain;      % Initially let some fraction melt
VRemain = (1-MeltFrac)*VRemain;   % This is fraction remaining after step
AllValues = Values;

for i=1:NTimes  % Perform a total of NTimes time steps
% First we let the current water flow downstrem
NewValues = Values*FracSame;  % New values initially a fraction of original
for j=1:ncells  % For each cell except leftmost or rightmost
if j>1
NewValues(j-1) = NewValues(j-1)+Values(j)*FracLeft;  % Diffuse Left
end
if j<ncells
NewValues(j+1) = NewValues(j+1)+Values(j)*FracRight; % Diffuse Right
end
end % for j=1...
Values = NewValues;     % The updated values become the current values

%  Now we will add some new melted water to the river
Melted = MeltFrac*VRemain;   % This is how much melted
VRemain = VRemain - Melted;  % This is how much remains
Values = Values + Melted;  % Add the melted

AllValues = [AllValues Values];
end  % for i=1...
```
3. Finally, the left, right and same fractions can be treated as arrays and a constriction can be placed some location (in this case cell 10). After running diff1vec.m, view the results using valsurf and valmovie
 Comparing the runs shown here, one notes that the crest value and the length of flooding is significantly higher at the location of the constriction (in this case the constriction and the plot are both at location 10). Comparing the runs shown here, one notes that the crest value is significantly lower and occurs later downstream of the constriction (in this case the constriction is at location 10 and the plots is at location 10). The surface plot (from valsurf) also allows one to see the effect of the constriction be noticing that diffusion occurs very slowly beyond cell 10
```      % DIFF1VEC.M    1-Dimensional Diffusion Solution with a Reservoir
%               at each location and the diffusion fractions
%               stored in vectors
%
% Tom Huber, March 1997
%
ncells = 25;       % Number of Cells in the Array
NTimes = 35;       % Number of Time steps to perform

FracLeft  = 1/10*ones(ncells,1);  % Fraction diffusing to the left
FracSame  = 3/10*ones(ncells,1);  % Fraction which remains in same cell
FracRight = 6/10*ones(ncells,1);  % Fraction diffusing to the right

% Make a constriction at one cell (cell number 10)
FracLeft(10) = 2/10;
FracSame(10) = 6/10;
FracRight(10) = 2/10;

MeltFrac = 0.4;            % Fraction which "Melts" each day

VRemain = zeros(ncells,1);      % Initialize all reservoirs to zero

%  Now load 5 reservoirs (locations 4-8)
%  with 20 instead of 1 reservoir with 100
VRemain(4:8) = 20*ones(5,1);

Values = MeltFrac*VRemain;      % Initially let some fraction melt
VRemain = (1-MeltFrac)*VRemain;   % This is fraction remaining after step
AllValues = Values;

for i=1:NTimes  % Perform a total of NTimes time steps
% First we let the current water flow downstrem
% NOTE: Sinc FracSame is now a vector, we use the .* operator
NewValues = Values.*FracSame;  % New values initially a fraction of original

% NOTE: We must use subsets of the vectors of length ncells-1 and
%  add them either to the left or right

NewValues(1:ncells-1)=NewValues(1:ncells-1)+FracLeft(2:ncells).*Values(2:ncells);
NewValues(2:ncells)=NewValues(2:ncells)+FracRight(1:ncells-1).*Values(1:ncells-1);
Values = NewValues;     % The updated values become the current values

%  Now we will add some new melted water to the river
Melted = MeltFrac*VRemain;   % This is how much melted
VRemain = VRemain - Melted;  % This is how much remains
Values = Values + Melted;  % Add the melted

AllValues = [AllValues Values];
end  % for i=1...
```

Electronic Copy: http://physics.gac.edu/~huber/envision/tutor2/diffvian.htm
Created: 4-APR-1997 by Tom Huber, Physics Department, Gustavus Adolphus College.