Reworked structure of argument-passing to operators.

operators.py

@@ -3,8 +3,10 @@
 import scipy.sparse as sparse
 
 
-def dyOp(nx, ny):
+def dyOp():
     # TODO
+    nx = config.nx
+    ny = config.ny
     dy = config.dy
 
     plus1 = np.ones((ny,))
@@ -24,35 +26,45 @@ def dyOp(nx, ny):
     return sparse.csr_matrix(sparse.diags([plus1[:-1], mid, min1[:-1]], [1, 0, -1]) / (2.0 * dy))
 
 
-def dyOpTemp(nx, ny, bcDirArray):
-    # Seems to work
+def dyOpTemp(bcDirArray):
+    # the dT/dy operator for the temperature field. Returns matrix and RHS. To calculate the dT/dy field one would
+    # perform A @ T - RHS
+
+    nx = config.nx
+    ny = config.ny
     dy = config.dy
 
+    # First construct the right hand side from one-sided difference formulas for y = 0 & y = ymax and the Dirichlet BC.
     rhs = np.zeros((nx * ny,))
     for i in range(nx):
         rhs[0 + ny * i] = bcDirArray[1] / dy  # y = 0
         rhs[ny * (i + 1) - 1] = - bcDirArray[0] / dy  # y = ymax
     rhs = np.expand_dims(rhs, 1)  # Make it an explicit column vector
 
+    # We use the central difference formula for all the inner derivatives ...
     plus1 = np.ones((ny,))
-    plus1[-1] = 0
-    plus1[0] = 2
     mid = np.zeros((ny,))
-    # mid[0] = -2
-    # mid[-1] = 2
     min1 = -np.ones((ny,))
+
+    # ... and again the one-sided difference formulas for y = 0 & y = ymax
+    plus1[-1] = 0
+    plus1[0] = 2
     min1[-1] = 0
     min1[-2] = -2
 
+    # Repeat along x = 0, 1, 2, ...
     plus1 = np.tile(plus1, (nx,))
     mid = np.tile(mid, (nx,))
     min1 = np.tile(min1, (nx,))
 
+    # Construct a sparse matrix from diagonals
     return sparse.csr_matrix(sparse.diags([plus1[:-1], mid, min1[:-1]], [1, 0, -1]) / (2.0 * dy)), rhs
 
 
-def dxOpTemp(nx, ny, bcNeuArray):
+def dxOpTemp(bcNeuArray):
     # Seems to work
+    nx = config.nx
+    ny = config.ny
     dx = config.dx
 
     rhs = np.zeros((nx * ny,))
@@ -75,7 +87,9 @@ def dxOpTemp(nx, ny, bcNeuArray):
     return sparse.csr_matrix(sparse.diags([plus1[:], mid, min1[:]], [ny, 0, -ny]) / (2.0 * dx)), rhs
 
 
-def dxOpStream(nx, ny):
+def dxOpStream():
+    nx = config.nx
+    ny = config.ny
     dx = config.dx
 
     plus1 = np.ones((ny,))
@@ -95,8 +109,10 @@ def dxOpStream(nx, ny):
     return sparse.csr_matrix(sparse.diags([plus1[:], mid, min1[:]], [ny, 0, -ny]) / (2.0 * dx))
 
 
-def dlOpStream(nx, ny):
+def dlOpStream():
     # Streamline Laplace operator
+    nx = config.nx
+    ny = config.ny
     dx = config.dx
     dy = config.dy
 
@@ -142,9 +158,11 @@ def dlOpStream(nx, ny):
     return sparse.csr_matrix(sparse.diags([mid, plus1y, min1y, plus1x, min1x], [0, 1, -1, ny, -ny]))
 
 
-def dlOpTemp(nx, ny, bcDirArray, bcNeuArray):
+def dlOpTemp(bcDirArray, bcNeuArray):
     # TODO
     # Streamline Laplace operator
+    nx = config.nx
+    ny = config.ny
     dx = config.dx
     dy = config.dy