some clean-up and docstringing

Author: Michael Zingale

advection/advectiveFluxes.py

@@ -1,74 +1,90 @@
-"""
-solve
+import numpy
 
-   a  + u a  + v a  = 0
-    t      x      y
+import mesh.reconstruction_f as reconstruction_f
 
-"""
+def unsplitFluxes(my_data, dt, scalar_name):
+    """
+    Construct the fluxes through the interfaces for the linear advection 
+    equation:
 
-import numpy
+      a  + u a  + v a  = 0
+       t      x      y
 
-import mesh.reconstruction_f as reconstruction_f
+    We use a second-order (piecewise linear) unsplit Godunov method 
+    (following Colella 1990).
 
-def unsplitFluxes(my_data, dt, a):
+    In the pure advection case, there is no Riemann problem we need to
+    solve -- we just simply do upwinding.  So there is only one 'state'
+    at each interface, and the zone the information comes from depends
+    on the sign of the velocity.
+    
+    Our convection is that the fluxes are going to be defined on the 
+    left edge of the computational zones
+
+            
+     |             |             |             |
+     |             |             |             |
+    -+------+------+------+------+------+------+--
+     |     i-1     |      i      |     i+1     |
+    
+              a_l,i  a_r,i   a_l,i+1
+            
+    
+    a_r,i and a_l,i+1 are computed using the information in
+    zone i,j.
+
+    Parameters
+    ----------
+    my_data : CellCenterData2d object
+        The data object containing the grid and advective scalar that
+        we are advecting.
+    dt : float
+        The timestep we are advancing through.
+    scalar_name : str
+        The name of the variable contained in my_data that we are
+        advecting
+
+    Returns
+    -------
+    out : ndarray, ndarray
+        The fluxes on the x- and y-interfaces
+
+    """
 
     my_grid = my_data.grid
     rp = my_data.rp
 
+    a = my_data.get_var(scalar_name)
+
     # get the advection velocities
     u = rp.get_param("advection.u")
     v = rp.get_param("advection.v")
 
     cx = u*dt/my_grid.dx
     cy = v*dt/my_grid.dy
-    
+
+    qx = my_grid.qx
+    qy = my_grid.qy    
 
     #--------------------------------------------------------------------------
     # monotonized central differences
     #--------------------------------------------------------------------------
 
     limiter = rp.get_param("advection.limiter")
-    if (limiter == 0):
+    if limiter == 0:
         limitFunc = reconstruction_f.nolimit
-    elif (limiter == 1):
+    elif limiter == 1:
         limitFunc = reconstruction_f.limit2
     else:
         limitFunc = reconstruction_f.limit4
 
-    ldelta_a = limitFunc(1, a, my_grid.qx, my_grid.qy, my_grid.ng)
+    ldelta_a = limitFunc(1, a, qx, qy, my_grid.ng)
 
-
-
-    # in the pure advection case, there is no Riemann problem we need to
-    # solve -- we just simply do upwinding.  So there is only one 'state'
-    # at each interface, and the zone the information comes from depends
-    # on the sign of the velocity -- this should be vectorized.
-    a_x = numpy.zeros((my_grid.qx,my_grid.qy), dtype=numpy.float64)
-    
-    """
-    
-    the fluxes are going to be defined on the left edge of the
-    computational zones
-
-            
-     |             |             |             |
-     |             |             |             |
-    -+------+------+------+------+------+------+--
-     |     i-1     |      i      |     i+1     |
-    
-              a_l,i  a_r,i   a_l,i+1
-            
-    
-    a_r,i and a_l,i+1 are computed using the information in
-    zone i,j.
+    a_x = my_grid.scratch_array()
 
-    """
-
-    qx = my_grid.qx
-    qy = my_grid.qy
 
     # upwind
-    if (u < 0):
+    if u < 0:
         # a_x[i,j] = a[i,j] - 0.5*(1.0 + cx)*ldelta_a[i,j]
         a_x[:,:] = a[:,:] - 0.5*(1.0 + cx)*ldelta_a[:,:]
     else:
@@ -77,13 +93,13 @@ def unsplitFluxes(my_data, dt, a):
 
 
     # y-direction
-    ldelta_a = limitFunc(2, a, my_grid.qx, my_grid.qy, my_grid.ng)
+    ldelta_a = limitFunc(2, a, qx, qy, my_grid.ng)
 
-    a_y = numpy.zeros((my_grid.qx,my_grid.qy), dtype=numpy.float64)
+    a_y = my_grid.scratch_array()
 
 
     # upwind
-    if (v < 0):
+    if v < 0:
         # a_y[i,j] = a[i,j] - 0.5*(1.0 + cy)*ldelta_a[i,j]
         a_y[:,:] = a[:,:] - 0.5*(1.0 + cy)*ldelta_a[:,:]
     else:
@@ -96,8 +112,8 @@ def unsplitFluxes(my_data, dt, a):
     F_xt = u*a_x
     F_yt = v*a_y
 
-    F_x = numpy.zeros((my_grid.qx,my_grid.qy), dtype=numpy.float64)
-    F_y = numpy.zeros((my_grid.qx,my_grid.qy), dtype=numpy.float64)
+    F_x = my_grid.scratch_array()
+    F_y = my_grid.scratch_array()
 
     # the zone where we grab the transverse flux derivative from
     # depends on the sign of the advective velocity
@@ -107,7 +123,6 @@ def unsplitFluxes(my_data, dt, a):
     else:
         mx = -1
 
-
     if v <= 0:
         my = 0
     else:
@@ -118,29 +133,26 @@ def unsplitFluxes(my_data, dt, a):
     dtdy2 = 0.5*dt/my_grid.dy
 
     i = my_grid.ilo
-    while (i <= my_grid.ihi+1):
+    while i <= my_grid.ihi+1:
 
         j = my_grid.jlo
-        while (j <= my_grid.jhi+1):
+        while j <= my_grid.jhi+1:
 
             F_x[i,j] = u*(a_x[i,j] - dtdy2*(F_yt[i+mx,j+1] - F_yt[i+mx,j]))
 
             j += 1
         i += 1
 
     i = my_grid.ilo
-    while (i <= my_grid.ihi+1):
+    while i <= my_grid.ihi+1:
 
         j = my_grid.jlo
-        while (j <= my_grid.jhi+1):
+        while j <= my_grid.jhi+1:
 
             F_y[i,j] = v*(a_y[i,j] - dtdx2*(F_xt[i+1,j+my] - F_xt[i,j+my]))
 
             j += 1
         i += 1
 
 
-    return [F_x, F_y]
-
-
-    
+    return F_x, F_y

advection/evolve.py

@@ -3,13 +3,10 @@
 def evolve(my_data, dt):
     """ evolve the advection equations through one timestep """
 
-    
-    dens = my_data.get_var("density")
-
     dtdx = dt/my_data.grid.dx
     dtdy = dt/my_data.grid.dy
 
-    [flux_x, flux_y] =  unsplitFluxes(my_data, dt, dens)
+    flux_x, flux_y =  unsplitFluxes(my_data, dt, "density")
 
     """
     do the differencing for the fluxes now.  Here, we use slices so we
@@ -22,9 +19,9 @@ def evolve(my_data, dt):
 
     qx = my_data.grid.qx
     qy = my_data.grid.qy
+
+    dens = my_data.get_var("density")
+
     dens[0:qx-1,0:qy-1] = dens[0:qx-1,0:qy-1] + \
         dtdx*(flux_x[0:qx-1,0:qy-1] - flux_x[1:qx,0:qy-1]) + \
-        dtdy*(flux_y[0:qx-1,0:qy-1] - flux_y[0:qx-1,1:qy])
-            
-
-    
+        dtdy*(flux_y[0:qx-1,0:qy-1] - flux_y[0:qx-1,1:qy])