can now also choose between order 2 and 4 for non-advection terms; simple test for PFASST on eight processors with operator coarsening indicates reasonable convergence

Author: Daniel Ruprecht

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/buildFDMatrix.py

@@ -31,7 +31,7 @@ def getUpwindMatrix(N, dx, order):
     coeff    = 1.0/60.0
     zero_pos = 5
   else:
-    sys.exit("Order "+str(order)+" not implemented.")
+    sys.exit('Order '+str(order)+' not implemented')
 
   first_col = np.zeros(N)
 
@@ -44,78 +44,144 @@ def getUpwindMatrix(N, dx, order):
   return sp.csc_matrix( coeff*(1.0/dx)*la.circulant(first_col) )
 
 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 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

examples/boussinesq_2d_imex/playground.py

@@ -24,7 +24,7 @@
     # set global logger (remove this if you do not want the output at all)
     logger = Log.setup_custom_logger('root')
 
-    num_procs = 1
+    num_procs = 8
 
     # This comes as read-in for the level class
     lparams = {}
@@ -48,16 +48,17 @@
 
     # This comes as read-in for the problem class
     pparams = {}
-    pparams['nvars']    = [(4,450,30)]
+    pparams['nvars']    = [(4,450,30), (4,450,30)]
     pparams['u_adv']    = 0.02
     pparams['c_s']      = 0.3
     pparams['Nfreq']    = 0.01
     pparams['x_bounds'] = [(-150.0, 150.0)]
     pparams['z_bounds'] = [(   0.0,  10.0)]
-    pparams['order']    = [1] # [fine_level, coarse_level]
-    pparams['gmres_maxiter'] = [50]
-    pparams['gmres_restart'] = 20
-    pparams['gmres_tol']     = 1e-6
+    pparams['order']    = [4, 2] # [fine_level, coarse_level]
+    pparams['order_upw'] = [5, 1]
+    pparams['gmres_maxiter'] = [50, 50]
+    pparams['gmres_restart'] = [20, 20]
+    pparams['gmres_tol']     = [1e-6, 1e-6]
 
     # This comes as read-in for the transfer operations
     tparams = {}
@@ -73,8 +74,8 @@
     description['sweeper_class']     = imex_1st_order
     description['level_params']      = lparams
     description['hook_class']        = plot_solution
-    #description['transfer_class']    = mesh_to_mesh_2d
-    #description['transfer_params']   = tparams
+    description['transfer_class']    = mesh_to_mesh_2d
+    description['transfer_params']   = tparams
 
     # quickly generate block of steps
     MS = mp.generate_steps(num_procs,sparams,description)