When given the task to change my explicit ocean model to an implicit version, I had assumed that I would need to invert the ‘A’ matrix to solve Ax=B, and then do some matrix math. I found a few web sites that warn against that because it is possibly a lot more computationally intensive than solving directly. John Cook’s blog explains more about why.
Python’s NumPy has linalg.solve(A, B), which returns the ‘x’ array
The next assignment was to rewrite the sliding plates using an implicit scheme. Here is a video of the animation produced using Visual Python arrows.
Colored according to speed. The initial velocity of the fluid is zero and the upper plate moves at a constant speed.
The visualization setup using Visual Python looks like this:
#SET UP VISUALIZATION SCENE
scene = display(title='Velocity profile', width=800, height=800, center=(u/2,y_divs/2,0),
autoscale=True, fullscreen=False)
base = box(pos=(u/2,-0.5,-2), length=u, height=0.4, width=4)
plate = box(pos=(u/2,y_divs+1.5,-2), length=u, height=0.4, width=4)
quiver = [] #TO HOLD ALL THE VELOCITY VECTOR ARROWS
quiver.append( arrow(pos=(0,0,0), axis=(BOUNDS[0,0],0,0),
shaftwidth=0.5, color=color.blue) )
for each in xrange(y_divs):
quiver.append( arrow(pos=(0,each+1,0), axis=(UM[each,0],0,0),
shaftwidth=0.5, color=color.blue) )
quiver.append( arrow(pos=(0,y_divs+1,0), axis=(BOUNDS[-1,0],0,0),
shaftwidth=0.5, color=color.red) )
scene.autoscale = False
The main loop with the visualization update looks like this:
for t in arange(t_divs-1): #START FROM 1, SKIP ZERO
B = zeros(y_divs)
B[:] = -1 * UM[:,t]
B[ 0] -= BOUNDS[ 0,t+1] * alpha
B[-1] -= BOUNDS[-1,t+1] * alpha
UMnext = linalg.solve( A, B )
UM[:,t+1] = UMnext[:]
for y in arange(y_divs):
#UPDATE VECTOR LENGTH BASED ON SPEED
quiver[y+1].axis = (UM[y,t+1],0,0)
#CHANGE ARROW COLOR BASED ON SPEED: Blue is slow, Red is fast
quiver[y+1].color = (UM[y,t+1]/u,0.3,1.-UM[y,t+1]/u)
Completed the sliding plates assignment.
Assignment: Simulate the vertical velocity profile of a liquid between two infinitely long plates separated by 40mm. The bottom plate remains fixed, but the top plates begins moving at a constant 40m/s in the +x direction. Using the explicit finite difference method.
Visual python makes it so easy to make a nice graphical simulation. I use arrows along the profile that change from blue to red as they increase in length to show the velocity.
.index() is the builtin function for finding the index of something in a list or tuple. This does not work with NumPy arrays! NP arrays rely on where() to return an array of indices that match the condition.
—————————
for-else statement
‘else’ can be used after a for-loop. But here, the ‘else’ statement only executes if the for-loop completes all iterations. Think of it as using the for-loop to search for something, or expecting a break from the loop. If nothing is resolved or no break occurs then the else statement can be used for notification. The else can be used after a try-except in a similar way; executing only when nothing breaks down, throwing an exception.
—————————
how to prevent a new-line character in a print() statement
Use a comma at the end. You can have consecutive print() statements print to the same line simply by adding a comma at the end of each line.
—————————
Use a dictionary for formatting a long output string
%20s will format a string in the list after a text. But when using this formatting method, the order must match in the supplied list. %(keyword)20s or %(keyword)10f with a supplied dictionary can avoid the ordering issue and having to have the same variable or value in a list multiple times.
—————————
use an asterisk * in front of a functions input argument to accept an infinite number of args.
def MyFunction( *args ): This will accept any number of arguments and create a list to use within the function.
—————————
ternary operation or inline if (iif)
In Java you can use a iif like this ( x == 2 ? y = 4 : y = 0 ).
( condition ? if true do this : if false do this)
I found that Python also has a version of this but ordered differently:
(do this) if (condition) else (do this if false)
Not the same order as java but still pretty straight forward.
This is what I have so far for a 1D Shallow Water model. Based on the lecture from class and also an ocean modeling book. This version has no friction parameter.
For the first entry of my programming activities record , I will start by posting my current script template. I plan on periodically updating this as I learn more about Python, but I guessing it won’t change much from this format. Maybe I should put it on a page instead of post.
_____________template.py_______________
#!/usr/bin/env python
# -*- coding: utf-8 -*-
“””
Created on %(date)s
@author: Ripley6811
@abstract:
@input: none
@output: none
@acknowledgments:
“””
#from numpy import * #imports ndarray(), arange(), zeros(), ones(), and other array methods
#from visual import * #imports NumPy.*, SciPy.*, and Visual objects (sphere, box, etc.)
#import matplotlib.pyplot as plot #use plot(x,y) and show() for plotting
#from pylab import * #imports NumPy.*, SciPy.*, and matplotlib.*
def main():
“””Description of main()”””
if __name__ == ‘__main__’:
main()
”’######################################
linspace(0,2,9) > array of 9 numbers from 0 to 2
zeros( (10,10,10) ) > a 3-d array of 0.0 (pass a list)
ones( (3,4,5) ) > a 3-d array of 1.0’s
”’######################################
The Python Path 
