1d acoustic now running correctly and showing expected convergence

Author: Daniel Ruprecht

examples/acoustic_1d_imex/ProblemClass.py

@@ -17,11 +17,7 @@
 from pySDC.Problem import ptype
 from pySDC.datatype_classes.mesh import mesh, rhs_imex_mesh
 
-# Sharpclaw imports
-from clawpack import pyclaw
-from clawpack import riemann
-
-from getFDMatrix import getFDMatrix
+from buildWave1DMatrix import getWave1DMatrix, getWave1DAdvectionMatrix
 
 def u_initial(x):
     return np.sin(x)
@@ -60,29 +56,13 @@ def __init__(self, cparams, dtype_u, dtype_f):
         # invoke super init, passing number of dofs, dtype_u and dtype_f
         super(acoustic_1d_imex,self).__init__(self.nvars,dtype_u,dtype_f)
 
-        riemann_solver              = riemann.advection_1D # NOTE: This uses the FORTRAN kernels of clawpack
-        self.solver                 = pyclaw.SharpClawSolver1D(riemann_solver)
-        self.solver.weno_order      = 5
-        self.solver.time_integrator = 'Euler' # Remove later
-        self.solver.kernel_language = 'Fortran'
-        self.solver.bc_lower[0]     = pyclaw.BC.periodic
-        self.solver.bc_upper[0]     = pyclaw.BC.periodic
-        self.solver.cfl_max         = 1.0
-        assert self.solver.is_valid()
-
-        x = pyclaw.Dimension(0.0, 2*np.pi, self.nvars[1], name='x')
-        self.domain = pyclaw.Domain(x)
-        self.state  = pyclaw.State(self.domain, self.solver.num_eqn)
-        self.mesh   = self.state.grid.x.centers
+        self.mesh   = np.linspace(0.0, 2.0*np.pi, self.nvars[1], endpoint=False)
         self.dx     = self.mesh[1] - self.mesh[0]
-        self.A      = -self.cs*getFDMatrix(self.nvars[1], self.order_adv, self.dx)
-        
-        self.state.problem_data['u'] = self.cadv
 
-        solution = pyclaw.Solution(self.state, self.domain)
-        self.solver.setup(solution)
-
-
+        self.Dx     = -self.cadv*getWave1DAdvectionMatrix(self.nvars[1], self.dx, self.order_adv)
+        self.Id, A  = getWave1DMatrix(self.nvars[1], self.dx, ['periodic','periodic'], ['periodic','periodic'])
+        self.A = -self.cs*A
+                
     def solve_system(self,rhs,factor,u0,t):
         """
         Simple linear solver for (I-dtA)u = rhs
@@ -97,9 +77,7 @@ def solve_system(self,rhs,factor,u0,t):
             solution as mesh
         """
 
-        M1 = sp.hstack( (sp.eye(self.nvars[1]), -factor*self.A) )
-        M2 = sp.hstack( (-factor*self.A, sp.eye(self.nvars[1])) )
-        M  = sp.vstack( (M1, M2) )
+        M  = self.Id - factor*self.A
 
         b = np.concatenate( (rhs.values[0,:], rhs.values[1,:]) )
 
@@ -124,26 +102,12 @@ def __eval_fexpl(self,u,t):
         """
 
 
+        b = np.concatenate( (u.values[0,:], u.values[1,:]) )
+        sol = self.Dx.dot(b)
+        
         fexpl        = mesh(self.nvars)
+        fexpl.values[0,:], fexpl.values[1,:] = np.split(sol, 2)
 
-        # Copy values of u into pyClaw state object
-        self.state.q[0,:] = u.values[0,:]
-
-        # Evaluate right hand side
-        tmp = self.solver.dqdt(self.state)
-        fexpl.values[0,:] = tmp.reshape(self.nvars[1:])
-
-        # Copy values of u into pyClaw state object
-        self.state.q[0,:] = u.values[1,:]
-
-        # Evaluate right hand side
-        tmp = self.solver.dqdt(self.state)
-        fexpl.values[1,:] = tmp.reshape(self.nvars[1:])
-        
-        
-        # DEBUGGING
-        # fexpl.values[0,:] = 0.0*self.mesh
-        # fexpl.values[1,:] = 0.0*self.mesh
         return fexpl
 
 
@@ -159,9 +123,11 @@ def __eval_fimpl(self,u,t):
             implicit part of RHS
         """
 
+        b = np.concatenate( (u.values[0,:], u.values[1,:]) )
+        sol = self.A.dot(b)
+        
         fimpl             = mesh(self.nvars,val=0)
-        fimpl.values[0,:] = self.A.dot(u.values[1,:])
-        fimpl.values[1,:] = self.A.dot(u.values[0,:])
+        fimpl.values[0,:], fimpl.values[1,:] = np.split(sol, 2)
 
         return fimpl
 

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 "+order+" not implemented."
+    print "Order "+str(order)+" not implemented."
 
   first_col = np.zeros(N)
 

examples/acoustic_1d_imex/plotconvdata.py

@@ -15,10 +15,6 @@
   nsteps = np.append(nsteps, int(items[1]))
   error  = np.append(error,  float(items[2]))
 
-print order
-print nsteps
-print error
-
 assert np.size(order)==np.size(nsteps), 'Found different number of entries in order and nsteps'
 assert np.size(nsteps)==np.size(error), 'Found different number of entries in nsteps and error'
 
@@ -31,11 +27,11 @@
 order_plot  = np.zeros(3)
 
 for ii in range(0,3):
-  order_plot[ii] = order[3*ii]
+  order_plot[ii] = order[N*ii]
   for jj in range(0,N):
-    error_plot[ii,jj]  = error[3*ii+jj]
-    nsteps_plot[ii,jj] = nsteps[3*ii+jj]
-    convline[ii,jj]    = error_plot[ii,0]*(float(nsteps_plot[ii,0])/float(nsteps_plot[ii,jj]))**order[3*ii+jj]
+    error_plot[ii,jj]  = error[N*ii+jj]
+    nsteps_plot[ii,jj] = nsteps[N*ii+jj]
+    convline[ii,jj]    = error_plot[ii,0]*(float(nsteps_plot[ii,0])/float(nsteps_plot[ii,jj]))**order_plot[ii]
 
 color = [ 'r', 'b', 'g' ]
 fig = plt.figure(figsize=(8,8))
@@ -44,6 +40,9 @@
   plt.loglog(nsteps_plot[ii,:], convline[ii,:], '-', color=color[ii])
 
 plt.legend()
-plt.xlabel(r'$N_t$')
+plt.xlabel(r'Number of time step $N_t$')
 plt.ylabel('Relative error')
+plt.xlim([0.9*np.min(nsteps_plot), 1.1*np.max(nsteps_plot)])
 plt.show()
+fig.savefig('sdc_fwsw_convergence.pdf',bbox_inches='tight')
+

examples/acoustic_1d_imex/runconvergence.py

@@ -34,11 +34,11 @@
 
     # This comes as read-in for the problem class
     pparams = {}
-    pparams['nvars']     = [(2,500)]
-    pparams['cadv']      = 0.1
+    pparams['nvars']     = [(2,250)]
+    pparams['cadv']      = 0.05
     pparams['cs']        = 1.0
-    pparams['order_adv'] = 4
-    
+    pparams['order_adv'] = 5
+
     # This comes as read-in for the transfer operations
     tparams = {}
     tparams['finter'] = True
@@ -56,20 +56,26 @@
     #description['transfer_class'] = mesh_to_mesh_1d
     #description['transfer_params'] = tparams
 
-    Nsteps = [20, 40, 50, 80]
-    order  = 2
+    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 = 5.5
+    Tend = 1.0
 
     if order==2:
       file = open('conv-data.txt', 'w')
     else:
       file = open('conv-data.txt', 'a')
 
-    description['num_nodes'] = (order+2)/2
+    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