Merge pull request #37 from danielru/boussinesq/selectable_order

Boussinesq/selectable order

Author: pancetta

examples/boussinesq_2d_imex/ProblemClass.py

@@ -62,6 +62,7 @@ def __init__(self, cparams, dtype_u, dtype_f):
         assert 'Nfreq' in cparams
         assert 'x_bounds' in cparams
         assert 'z_bounds' in cparams
+        assert 'order_upw' in cparams
         assert 'order' in cparams
         assert 'gmres_maxiter' in cparams
         assert 'gmres_restart' in cparams
@@ -82,7 +83,7 @@ def __init__(self, cparams, dtype_u, dtype_f):
         self.xx, self.zz, self.h = get2DMesh(self.N, self.x_bounds, self.z_bounds, self.bc_hor[0], self.bc_ver[0])
 
         self.Id, self.M = getBoussinesq2DMatrix(self.N, self.h, self.bc_hor, self.bc_ver, self.c_s, self.Nfreq, self.order)
-        self.D_upwind   = getBoussinesq2DUpwindMatrix( self.N, self.h[0], self.u_adv , self.order)
+        self.D_upwind   = getBoussinesq2DUpwindMatrix( self.N, self.h[0], self.u_adv , self.order_upw)
 
         self.logger = logging()
 

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,78 +43,145 @@ def getUpwindMatrix(N, dx):
 
   return sp.csc_matrix( coeff*(1.0/dx)*la.circulant(first_col) )
 
-def getMatrix(N, dx, bc_left, bc_right):
-  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)
+def getMatrix(N, dx, bc_left, bc_right, order):
 
   assert bc_left in ['periodic','neumann','dirichlet'], "Unknown type of BC"
+
+  if order==2:
+    stencil = [-1.0, 0.0, 1.0]
+    range   = [-1, 0, 1]
+    coeff   = 1.0/2.0
+  elif order==4:
+    stencil = [1.0, -8.0, 0.0, 8.0, -1.0]
+    range   = [ -2,   -1,   0,   1,    2]
+    coeff   = 1.0/12.0
+  else:
+    sys.exit('Order '+str(order)+' not implemented')
+
+  A       = sp.diags(stencil, range, shape=(N,N))
+  A       = sp.lil_matrix(A)
+  
+  #
+  # Periodic boundary conditions
+  #
   if bc_left in ['periodic']:
     assert bc_right in ['periodic'], "Periodic BC can only be selected for both sides simultaneously"
 
   if bc_left in ['periodic']:
-    A[0,N-2] = stencil[0]
-    A[0,N-1] = stencil[1]
-    A[1,N-1] = stencil[0]
+    if order==2:
+      A[0,N-1] = stencil[0]
 
-  if bc_right in ['periodic']:
-    A[N-2,0] = stencil[4]
-    A[N-1,0] = stencil[3]
-    A[N-1,1] = stencil[4]
+    elif order==4:
+      A[0,N-2] = stencil[0]
+      A[0,N-1] = stencil[1]
+      A[1,N-1] = stencil[0]
 
+  if bc_right in ['periodic']:
+    if order==2:
+      A[N-1,0] = stencil[2]
+    elif order==4:
+      A[N-2,0] = stencil[4]
+      A[N-1,0] = stencil[3]
+      A[N-1,1] = stencil[4]
+
+  #
+  # Neumann boundary conditions
+  #
   if bc_left in ['neumann']:
     A[0,:] = np.zeros(N)
-    A[0,0] = -8.0
-    A[0,1] =  8.0
-    A[1,0] = -8.0 + 4.0/3.0
-    A[1,1] = -1.0/3.0
+    if order==2:
+      A[0,0] = -4.0/3.0
+      A[0,1] = 4.0/3.0
+    elif order==4:
+      A[0,0] = -8.0
+      A[0,1] =  8.0
+      A[1,0] = -8.0 + 4.0/3.0
+      A[1,1] = -1.0/3.0
 
   if bc_right in ['neumann']:
     A[N-1,:]   = np.zeros(N)
-    A[N-2,N-1] = 8.0 - 4.0/3.0
-    A[N-2,N-2] = 1.0/3.0
-    A[N-1,N-1] =  8.0
-    A[N-1,N-2] = -8.0
-
+    if order==2:
+      A[N-1,N-2] = -4.0/3.0
+      A[N-1,N-1] = 4.0/3.0
+    elif order==4:
+      A[N-2,N-1] = 8.0 - 4.0/3.0
+      A[N-2,N-2] = 1.0/3.0
+      A[N-1,N-1] =  8.0
+      A[N-1,N-2] = -8.0
+
+  #
+  # Dirichlet boundary conditions
+  #
   if bc_left in ['dirichlet']:
-    A[0,:] = np.zeros(N)
-    A[0,1] = 6.0
+    # For order==2, nothing to do here
+    if order==4:
+      A[0,:] = np.zeros(N)
+      A[0,1] = 6.0
 
   if bc_right in ['dirichlet']:
-    A[N-1,:]   = np.zeros(N)
-    A[N-1,N-2] = -6.0
-
-  A = 1.0/(12.0*dx)*A
+    # For order==2, nothing to do here
+    if order==4:
+      A[N-1,:]   = np.zeros(N)
+      A[N-1,N-2] = -6.0
 
+  
+  A = coeff*(1.0/dx)*A
   return sp.csc_matrix(A)
 
-def getBCLeft(value, N, dx, type):
+#
+#
+#
+def getBCLeft(value, N, dx, type, order):
 
   assert type in ['periodic','neumann','dirichlet'], "Unknown type of BC"
 
+  if order==2:
+    coeff = 1.0/2.0
+  elif order==4:
+    coeff = 1.0/12.0
+
   b = np.zeros(N)
   if type in ['dirichlet']:
-    b[0] = -6.0*value
-    b[1] =  1.0*value
+    if order==2:
+      b[0] = -value;
+    elif order==4:
+      b[0] = -6.0*value
+      b[1] =  1.0*value
 
   if type in ['neumann']:
-    b[0] = 4.0*dx*value
-    b[1] = -(2.0/3.0)*dx*value
+    if order==2:
+      b[0] = (2.0/3.0)*dx*value
+    elif order==4:
+      b[0] = 4.0*dx*value
+      b[1] = -(2.0/3.0)*dx*value
 
-  return (1.0/(12.0*dx))*b
+  return coeff*(1.0/dx)*b
 
-def getBCRight(value, N, dx, type):
+#
+#
+#
+def getBCRight(value, N, dx, type, order):
 
   assert type in ['periodic','neumann','dirichlet'], "Unknown type of BC"
 
+  if order==2:
+    coeff = 1.0/2.0
+  elif order==4:
+    coeff = 1.0/12.0
+
   b = np.zeros(N)
   if type in ['dirichlet']:
-    b[N-2] = -1.0*value
-    b[N-1] =  6.0*value
+    if order==2:
+      b[N-1] = value
+    elif order==4:
+      b[N-2] = -1.0*value
+      b[N-1] =  6.0*value
 
   if type in ['neumann']:
-    b[N-2] = -(2.0/3.0)*dx*value
-    b[N-1] =  4.0*dx*value
+    if order==2:
+      b[N-1] = (2.0/3.0)*dx*value
+    elif order==4:
+      b[N-2] = -(2.0/3.0)*dx*value
+      b[N-1] =  4.0*dx*value
 
-  return (1.0/(12.0*dx))*b
+  return coeff*(1.0/dx)*b