order parameter now correctly passed down to buildFDMatrix; upwind operators available for order 1 to 5; buildFDMatrix still needs modification to support more than order 4

Author: Daniel Ruprecht

examples/boussinesq_2d_imex/build2DFDMatrix.py

@@ -4,11 +4,17 @@
 import scipy.sparse as sp
 from buildFDMatrix import getMatrix, getUpwindMatrix, getBCLeft, getBCRight
 
-def get2DUpwindMatrix(N, dx):
+#
+#
+#
+def get2DUpwindMatrix(N, dx, order):
 
-  Dx = getUpwindMatrix( N[0], dx)
+  Dx = getUpwindMatrix( N[0], dx, order)
   return sp.kron( Dx, sp.eye(N[1]), format="csr" )
 
+#
+#
+#
 def get2DMesh(N, x_b, z_b, bc_hor, bc_ver):
   assert np.size(N)==2, 'N needs to be an array with two entries: N[0]=Nx and N[1]=Nz'
   assert np.size(x_b)==2, 'x_b needs to be an array with two entries: x_b[0] = left boundary, x_b[1] = right boundary'
@@ -39,12 +45,16 @@ def get2DMesh(N, x_b, z_b, bc_hor, bc_ver):
   xx, zz = np.meshgrid(x,z,indexing="ij")
   return xx, zz, h
 
-def get2DMatrix(N, h, bc_hor, bc_ver):
+#
+#
+#
+def get2DMatrix(N, h, bc_hor, bc_ver, order):
+
   assert np.size(N)==2, 'N needs to be an array with two entries: N[0]=Nx and N[1]=Nz'
   assert np.size(h)==2, 'h needs to be an array with two entries: h[0]=dx and h[1]=dz'
 
-  Ax = getMatrix( N[0], h[0], bc_hor[0], bc_hor[1])
-  Az = getMatrix( N[1], h[1], bc_ver[0], bc_ver[1])
+  Ax = getMatrix( N[0], h[0], bc_hor[0], bc_hor[1], order)
+  Az = getMatrix( N[1], h[1], bc_ver[0], bc_ver[1], order)
 
   Dx = sp.kron( Ax, sp.eye(N[1]), format="csr")
   Dz = sp.kron( sp.eye(N[0]), Az, format="csr")
@@ -83,4 +93,4 @@ def getBCVertical( value, N, dz, bc_ver):
   bu = getBCRight( value[1], N[1], dz, bc_ver[1])
   bu = np.kron( np.ones(N[0]), bu)
 
-  return bd, bu
+  return bd, bu

examples/boussinesq_2d_imex/buildBoussinesq2DMatrix.py

@@ -4,7 +4,7 @@
 
 def getBoussinesq2DUpwindMatrix(N, dx, u_adv, order):
 
-  Dx   = get2DUpwindMatrix(N, dx)
+  Dx   = get2DUpwindMatrix(N, dx, order)
 
   # Note: In the equations it is u_t + u_adv* D_x u = ... so in order to comply with the form u_t = M u,
   # add a minus sign in front of u_adv
@@ -19,10 +19,12 @@ def getBoussinesq2DUpwindMatrix(N, dx, u_adv, order):
   return sp.csc_matrix(M)
 
 def getBoussinesq2DMatrix(N, h, bc_hor, bc_ver, c_s, Nfreq, order):
-  Dx_u, Dz_u = get2DMatrix(N, h, bc_hor[0], bc_ver[0])
-  Dx_w, Dz_w = get2DMatrix(N, h, bc_hor[1], bc_ver[1])
-  Dx_b, Dz_b = get2DMatrix(N, h, bc_hor[2], bc_ver[2])
-  Dx_p, Dz_p = get2DMatrix(N, h, bc_hor[3], bc_ver[3])
+
+  Dx_u, Dz_u = get2DMatrix(N, h, bc_hor[0], bc_ver[0], order)
+  Dx_w, Dz_w = get2DMatrix(N, h, bc_hor[1], bc_ver[1], order)
+  Dx_b, Dz_b = get2DMatrix(N, h, bc_hor[2], bc_ver[2], order)
+  Dx_p, Dz_p = get2DMatrix(N, h, bc_hor[3], bc_ver[3], order)
+
   Id_N = sp.eye(N[0]*N[1])
 
   Zero = np.zeros((N[0]*N[1],N[0]*N[1]))
@@ -53,4 +55,4 @@ def getBoussinesqBCHorizontal(value, N, dx, bc_hor):
   return b_left, b_right
 
 def getBoussinesqBCVertical():
-  return 0.0
+  return 0.0

examples/boussinesq_2d_imex/buildFDMatrix.py

@@ -1,29 +1,38 @@
 import numpy as np
 import scipy.sparse as sp
 import scipy.linalg as la
+import sys
+
 # Only for periodic BC because we have advection only in x direction
-def getUpwindMatrix(N, dx):
-     
-  #stencil    = [-1.0, 1.0]
-  #zero_pos = 2
-  #coeff      = 1.0
-  
-  #stencil    = [1.0, -4.0, 3.0]
-  #coeff      = 1.0/2.0
-  #zero_pos   = 3
-  
-  #stencil    = [1.0, -6.0, 3.0, 2.0]
-  #coeff      = 1.0/6.0
-  #zero_pos   = 3
+def getUpwindMatrix(N, dx, order):
 
-  #stencil  = [-5.0, 30.0, -90.0, 50.0, 15.0]
-  #coeff    = 1.0/60.0
-  #zero_pos = 4
+  if order==1:
+    stencil    = [-1.0, 1.0]
+    zero_pos   = 2
+    coeff      = 1.0
 
-  stencil = [3.0, -20.0, 60.0, -120.0, 65.0, 12.0]
-  coeff   = 1.0/60.0
-  zero_pos = 5
+  elif order==2:
+    stencil    = [1.0, -4.0, 3.0]
+    coeff      = 1.0/2.0
+    zero_pos   = 3
+
+  elif order==3:
+    stencil    = [1.0, -6.0, 3.0, 2.0]
+    coeff      = 1.0/6.0
+    zero_pos   = 3
+
+  elif order==4:
+    stencil  = [-5.0, 30.0, -90.0, 50.0, 15.0]
+    coeff    = 1.0/60.0
+    zero_pos = 4
 
+  elif order==5:
+    stencil  = [3.0, -20.0, 60.0, -120.0, 65.0, 12.0]
+    coeff    = 1.0/60.0
+    zero_pos = 5
+  else:
+    sys.exit("Order "+str(order)+" not implemented.")
+
   first_col = np.zeros(N)
 
   # Because we need to specific first column (not row) in circulant, flip stencil array
@@ -34,13 +43,14 @@ def getUpwindMatrix(N, dx):
 
   return sp.csc_matrix( coeff*(1.0/dx)*la.circulant(first_col) )
 
-def getMatrix(N, dx, bc_left, bc_right):
+def getMatrix(N, dx, bc_left, bc_right, order):
   stencil = [1.0, -8.0, 0.0, 8.0, -1.0]
   range   = [ -2,   -1,   0,   1,    2]
   A       = sp.diags(stencil, range, shape=(N,N))
   A       = sp.lil_matrix(A)
 
   assert bc_left in ['periodic','neumann','dirichlet'], "Unknown type of BC"
+  
   if bc_left in ['periodic']:
     assert bc_right in ['periodic'], "Periodic BC can only be selected for both sides simultaneously"
 

examples/boussinesq_2d_imex/playground.py

@@ -40,8 +40,8 @@
 
     # setup parameters "in time"
     t0     = 0
-    Tend   = 3000
-    Nsteps =  500
+    #Tend   = 3000
+    #Nsteps =  500
     Tend   = 30
     Nsteps =  5
     dt = Tend/float(Nsteps)
@@ -54,7 +54,7 @@
     pparams['Nfreq']    = 0.01
     pparams['x_bounds'] = [(-150.0, 150.0)]
     pparams['z_bounds'] = [(   0.0,  10.0)]
-    pparams['order']    = [0] # [fine_level, coarse_level]
+    pparams['order']    = [1] # [fine_level, coarse_level]
     pparams['gmres_maxiter'] = [50]
     pparams['gmres_restart'] = 20
     pparams['gmres_tol']     = 1e-6
@@ -105,4 +105,4 @@
 
     #extract_stats = grep_stats(stats,iter=-1,type='residual')
     #sortedlist_stats = sort_stats(extract_stats,sortby='step')
-    #print(extract_stats,sortedlist_stats)
+    #print(extract_stats,sortedlist_stats)