update finite difference matrix builders

Daniel Ruprecht · Daniel Ruprecht · commit 26890cf6cb4c · 2015-12-03T09:00:35.000Z
diff --git a/examples/acoustic_1d_imex/buildFDMatrix.py b/examples/acoustic_1d_imex/buildFDMatrix.py
@@ -30,7 +30,7 @@ def getHorizontalDx(N, dx, order):
     coeff    = 1.0/60.0
     zero_pos = 5
   else:
-    print "Order "+str(order)+" not implemented."
+    print "Order "+order+" not implemented."
 
   first_col = np.zeros(N)
   
@@ -42,79 +42,143 @@ def getHorizontalDx(N, dx, order):
 
   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
+
+  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
diff --git a/examples/acoustic_1d_imex/buildWave1DMatrix.py b/examples/acoustic_1d_imex/buildWave1DMatrix.py
@@ -4,12 +4,14 @@
 import scipy.sparse as sp
 from buildFDMatrix import getMatrix, getHorizontalDx, getBCLeft, getBCRight
 
+wave_order = 4
+
 def getWave1DMatrix(N, dx, bc_left, bc_right):
   
   Id = sp.eye(2*N)
 
-  D_u  = getMatrix(N, dx, bc_left[0], bc_right[0])
-  D_p  = getMatrix(N, dx, bc_left[1], bc_right[1])
+  D_u  = getMatrix(N, dx, bc_left[0], bc_right[0], wave_order)
+  D_p  = getMatrix(N, dx, bc_left[1], bc_right[1], wave_order)
   Zero = np.zeros((N,N))
   M1 = sp.hstack((Zero, D_p), format="csc")
   M2 = sp.hstack((D_u, Zero), format="csc")
@@ -25,11 +27,11 @@ def getWave1DAdvectionMatrix(N, dx, order):
   return sp.csc_matrix(M)
 
 def getWaveBCLeft(value, N, dx, bc_left):
-  bu = getBCLeft(value[0],  N, dx, bc_left[0])
-  bp = getBCLeft(value[1],  N, dx, bc_left[1])
+  bu = getBCLeft(value[0],  N, dx, bc_left[0], wave_order)
+  bp = getBCLeft(value[1],  N, dx, bc_left[1], wave_order)
   return np.concatenate((bp, bu))
 
 def getWaveBCRight(value, N, dx, bc_right):
-  bu = getBCRight(value[0], N, dx, bc_right[0])
-  bp = getBCRight(value[1], N, dx, bc_right[1])
+  bu = getBCRight(value[0], N, dx, bc_right[0], wave_order)
+  bp = getBCRight(value[1], N, dx, bc_right[1], wave_order)
   return np.concatenate((bp, bu))
diff --git a/examples/acoustic_1d_imex/plotconvdata.py b/examples/acoustic_1d_imex/plotconvdata.py
@@ -1,7 +1,7 @@
 import numpy as np
 from matplotlib import pyplot as plt
 
-fs = 18
+fs     = 18
 order  = np.array([])
 nsteps = np.array([])
 error  = np.array([])
diff --git a/examples/acoustic_1d_imex/runconvergence.py b/examples/acoustic_1d_imex/runconvergence.py
@@ -3,7 +3,7 @@
 
 import numpy as np
 
-from ProblemClass import acoustic_1d_imex
+from ProblemClass_conv import acoustic_1d_imex
 from examples.acoustic_1d_imex.HookClass import plot_solution
 
 from pySDC.datatype_classes.mesh import mesh, rhs_imex_mesh
@@ -12,9 +12,6 @@
 from pySDC import Log
 from pySDC.Stats import grep_stats, sort_stats
 
-# Sharpclaw imports
-from clawpack import pyclaw
-from clawpack import riemann
 from matplotlib import pyplot as plt
 
 
@@ -36,7 +33,7 @@
     pparams = {}
     pparams['nvars']     = [(2,250)]
     pparams['cadv']      = 0.05
-    pparams['cs']        = 1.0
+    pparams['cs']        = 1.00
     pparams['order_adv'] = 5
     pparams['waveno']            = 5
 
@@ -50,67 +47,68 @@
     description['problem_params']    = pparams
     description['dtype_u']           = mesh
     description['dtype_f']           = rhs_imex_mesh
-    description['collocation_class'] = collclass.CollGaussLobatto
+    description['collocation_class'] = collclass.CollGaussRadau_Right
     description['sweeper_class']     = imex_1st_order
     description['level_params']      = lparams
     description['hook_class']        = plot_solution
     #description['transfer_class'] = mesh_to_mesh_1d
     #description['transfer_params'] = tparams
     
-    Nsteps = [3, 4, 6, 8, 10, 12, 15, 18, 20]
-    order  = 4
-    error  = np.zeros(np.size(Nsteps))
-  
-    # setup parameters "in time"
-    t0   = 0
-    Tend = 1.0
+    Nsteps = [65]
     
-    if order==2:
-      file = open('conv-data.txt', 'w')
-    else:
-      file = open('conv-data.txt', 'a')
-
-    if order==2:
-      description['num_nodes'] = 2
-    elif order==3:
-      description['num_nodes'] = 3
-    elif order==4:
-      description['num_nodes'] = 3
+    for order in [4]:
 
-    sparams['maxiter'] = order
+      error  = np.zeros(np.size(Nsteps))
     
-     # quickly generate block of steps
-    MS = mp.generate_steps(num_procs,sparams,description)
-    
-    for ii in range(0,np.size(Nsteps)):
-    
-      dt = Tend/float(Nsteps[ii])
+      # setup parameters "in time"
+      t0   = 0
+      Tend = 1.0
+      
+      if order==2:
+        file = open('conv-data.txt', 'w')
+      else:
+        file = open('conv-data.txt', 'a')
+
+      if order==2:
+        description['num_nodes'] = 2
+      elif order==3:
+        description['num_nodes'] = 3
+      elif order==4:
+        description['num_nodes'] = 3
+
+      sparams['maxiter'] = order
+      
+       # quickly generate block of steps
+      MS = mp.generate_steps(num_procs,sparams,description)
+      for ii in range(0,np.size(Nsteps)):
+      
+        dt = Tend/float(Nsteps[ii])
 
-      # get initial values on finest level
-      P = MS[0].levels[0].prob
-      uinit = P.u_exact(t0)
+        # get initial values on finest level
+        P = MS[0].levels[0].prob
+        uinit = P.u_exact(t0)
 
-      # call main function to get things done...
-      uend,stats = mp.run_pfasst(MS, u0=uinit, t0=t0, dt=dt, Tend=Tend)
+        # call main function to get things done...
+        uend,stats = mp.run_pfasst(MS, u0=uinit, t0=t0, dt=dt, Tend=Tend)
 
-      # compute exact solution and compare
-      uex = P.u_exact(Tend)
+        # compute exact solution and compare
+        uex = P.u_exact(Tend)
 
-      error[ii] = np.linalg.norm(uex.values-uend.values,np.inf)/np.linalg.norm(uex.values,np.inf)
-      file.write(str(order)+"    "+str(Nsteps[ii])+"    "+str(error[ii])+"\n")
+        error[ii] = np.linalg.norm(uex.values-uend.values,np.inf)/np.linalg.norm(uex.values,np.inf)
+        file.write(str(order)+"    "+str(Nsteps[ii])+"    "+str(error[ii])+"\n")
 
-    file.close()
+      file.close()
 
     if np.size(Nsteps)==1:
       fig = plt.figure(figsize=(8,8))
 
-      plt.plot(P.state.grid.x.centers, uex.values[0,:],  '+', color='b', label='u (exact)')
-      plt.plot(P.state.grid.x.centers, uend.values[0,:], '-', color='b', label='u (SDC)')
-      plt.plot(P.state.grid.x.centers, uex.values[1,:],  '+', color='r', label='p (exact)')
-      plt.plot(P.state.grid.x.centers, uend.values[1,:], '-', color='r', label='p (SDC)')
+      plt.plot(P.mesh, uex.values[0,:],  '+', color='b', label='u (exact)')
+      plt.plot(P.mesh, uend.values[0,:], '-', color='b', label='u (SDC)')
+      plt.plot(P.mesh, uex.values[1,:],  '+', color='r', label='p (exact)')
+      plt.plot(P.mesh, uend.values[1,:], '-', color='r', label='p (SDC)')
       plt.legend()
-      plt.xlim([0, 1])
-      plt.ylim([-1, 1])
+      plt.xlim([0.0, 1.0])
+      plt.ylim([-1.0, 1.0])
       plt.show()
     else:
       for ii in range(0,np.size(Nsteps)):
