Visualization of Diffusion Model Results

In the diffusion model we have developed, we had only one way of visualizing the results, namely a series of plots of values versus time. Below we will use some of Matlab's powerful features to enhance our ability to understand and visualize the results.

In particular, we are going to modify the previous program so it stores all the profiles in a 2 dimensional array where each column represents the values (for example, the volume of water in each cell) at a subsequent time step. Then we will be able to use this result in a number of ways:
• Viewing it as a surface, we can see the trends in the flow.
• We can slice it so we can see the values versus position at any time, or the values versus time at any position.
• We can take successive slices and view them as a "movie" to watch the flow
First, copy your program program diff1.m to diff1b.m. Modify the program diff1b.m by doing the following:
1. Add a statement AllValues=Values; right before the first for loop. This will store the initial distribution in the array AllValues
2. Remove all of the plotting statements from near the end of the program. In their place, we will append the updated column of values onto the array of previous values by including the statement AllValues=[AllValues Values]; (this is similar to what we did in Steve McKelvey's tutorial on population growth)
Now when we run the program, it will create the array AllValues where successive columns of the array have the values versus position.

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')
```
Experiment with other views of the same surface by typing in a command such as view(150,30) (the first arguement is the rotation angle of the XY plane around the Z axis, and the second arguement is the view angle above the XY plane.
2. We can extract a row or column from the AllValues matrix using, for example Values=AllValues(10,:) to extract the 10th row or Values=AllValues(:,5) to extract the 5th column. Plot the resulting Values and interpret the results.
3. Make an M-File plotloc.m which will ask the user for a cell number (hint: use the input statement) and then will extract the values versus time for that location from AllValues. Make a plot of this (including labels), and also print out (using the disp command) the maximum value and time step when it occurs (to do this, use the max statement. In particular, the command
`          [ValMax LocMax] = max(Values);`
will store the maximum value of the array Values in the variable ValMax and the location in the array (the element number) in LocMax
4. Use the program above to answer questions such as: When does the crest (maximum) occur for a position 5 cells downstream from the initial concentration? What about 15 cells downstream? Does doubling the initial concentration double the crest at these locations?
5. Perform the command AllValues1=AllValues (so that the current array is stored in this variable). Modify your program diff1b.m in some way (for example, changing the fractions to left, right and same) and re-run the program. Now modify your M-file plotloc.m from the previous task to extract and plot (on the same graph) both the values from AllValues and AllValues1. This will allow you to directly compare the results of two different models.
6. Make a script file valmovie.m which extracts successive columns from AllValues and plots these one after another. This "movie" will be similar to the output from our original program diff1.m (from the previous handout). You could even plot both AllValues and AllValues1 on the same graph.
At this point, we have a method for running the model to produce an output file and then to visualize the results using the programs valsurf, plotloc and valmovie. We can now add additional details to our model to make it closer to the actual situation that would exist in a river. You may want to work as a group on some of the following tasks which extend the model.

1. We modeled our source as coming at one time from one location (this would be a good model for a sudden pollution spill at some position). A more realistic model for snow melting would be a reservoir which melts and puts this melting snow into the river. Make a new version of your program which does the following:
• Defines an initial reservoir of material
• During the iteration for each day (the for i=1:NTimes loop) calculate the outflow from the reservoir and add this to the Values array at the appropriate spot (such as ncells/4).
• Proceed with the rest of the calculations as before.
• The resulting AllValues array can be viewed with the same programs valsurf, plotloc and valmovie.
• How does the flooding compare with the results from an "instant" melting?
2. The above example could be extended to have a reservoir located at every location (a more realistic model where there would be some amount of snow at every point along the river) and then the program would do melting and outflow at each location. To do this, make the reservoir an array and calculate the flow into the river for every cell.
3. Another modification would be to store the fractions left, right and same in an array (one element for each cell). These could originally be set to the original values, but later you could set one cell with fractions left, right and same to be 3/10, 2/10 and 5/10 respectively. This would model a region where the flow slowed substantially (such as if there was a dam that limited the flow of water in the river) and see how this changes the resulting flooding both upstream and downstream of the dam.

Answers to these problems can be found in http://physics.gac.edu/~huber/envision/tutor2/diffvis.htm
Electronic Copy: http://physics.gac.edu/~huber/envision/tutor2/diffvis.htm
Revised: 6-June-1997 by Tom Huber, Physics Department, Gustavus Adolphus College.