Merge pull request #38 from danielru/boussinesq/selectable_order

Boussinesq/selectable order

Author: pancetta

examples/SWFW/playground.py

@@ -19,18 +19,21 @@
     logger = Log.setup_custom_logger('root')
 
     num_procs = 1
+    
+    N_s = 100
+    N_f = 125
 
     # This comes as read-in for the level class
     lparams = {}
     lparams['restol'] = 1E-12
 
     sparams = {}
-    sparams['maxiter'] = 2
+    sparams['maxiter'] = 4
 
     # This comes as read-in for the problem class
     pparams = {}
-    pparams['lambda_s'] = 1j*np.linspace(0,3,100)
-    pparams['lambda_f'] = 1j*np.linspace(0,8,100)
+    pparams['lambda_s'] = 1j*np.linspace(0.0, 4.0, N_s)
+    pparams['lambda_f'] = 1j*np.linspace(0.0, 8.0, N_f)
     pparams['u0'] = 1
 
     # Fill description dictionary for easy hierarchy creation
@@ -40,7 +43,8 @@
     description['dtype_u'] = mesh
     description['dtype_f'] = rhs_imex_mesh
     description['collocation_class'] = collclass.CollGaussLobatto
-    description['num_nodes'] = [2]
+    description['num_nodes'] = [3]
+    description['do_LU'] = False
     description['sweeper_class'] = imex_1st_order
     description['level_params'] = lparams
 
@@ -62,15 +66,19 @@
 
     uex = P.u_exact(Tend)
 
-    fig = plt.figure(figsize=(8,8))
-    plt.pcolor(pparams['lambda_s'].imag, pparams['lambda_f'].imag, np.absolute(uend.values).T,vmin=1,vmax=1.01)
+    fig = plt.figure(figsize=(4,4))
+    #pcol = plt.pcolor(pparams['lambda_s'].imag, pparams['lambda_f'].imag, np.absolute(uend.values).T, vmin=0.99, vmax=1.01)
+    pcol = plt.pcolor(pparams['lambda_s'].imag, pparams['lambda_f'].imag, np.absolute(uend.values).T, vmin=1.0, vmax=2.0)
     # plt.pcolor(np.imag(uend.values))
-    plt.colorbar()
-    plt.xlabel('$\Delta t \lambda_{slow}$', fontsize=18, labelpad=20)
+    pcol.set_edgecolor('face')
+    plt.plot([0, 4], [0, 4], color='w', linewidth=2.5)
+    #plt.colorbar()
+    plt.gca().set_xticks([0.0, 1.0, 2.0, 3.0])
+    plt.xlabel('$\Delta t \lambda_{slow}$', fontsize=18, labelpad=0.0)
     plt.ylabel('$\Delta t \lambda_{fast}$', fontsize=18)
 
-    plt.show()
-
+    #plt.show()
+    fig.savefig('sdc-fwsw-stability-K'+str(sparams['maxiter'])+'-M'+str(description['num_nodes'][0])+'.pdf', bbox_inches='tight')
     print('error at time %s: %s' %(Tend,np.linalg.norm(uex.values-uend.values,np.inf)/np.linalg.norm(uex.values,
                                                                                                      np.inf)))
 

examples/acoustic_1d_imex/HookClass.py

@@ -0,0 +1,49 @@
+from __future__ import division
+from pySDC.Hooks import hooks
+from pySDC.Stats import stats
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from matplotlib import pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+from matplotlib import cm
+from matplotlib.ticker import LinearLocator, FormatStrFormatter
+
+class plot_solution(hooks):
+
+    def __init__(self):
+        """
+        Initialization of output
+        """
+        super(plot_solution,self).__init__()
+
+        # add figure object for further use
+        #        self.fig = plt.figure(figsize=(18,6))
+        #self.fig = plt.figure(figsize=(9,9))
+
+
+    def dump_step(self,status):
+        """
+        Overwrite standard dump per step
+
+        Args:
+            status: status object per step
+        """
+        super(plot_solution,self).dump_step(status)
+
+        if False:
+          yplot = self.level.uend.values
+          xx    = self.level.prob.mesh
+          self.fig.clear()
+          plt.plot(xx, yplot[0,:])
+          plt.axes().set_xlim(xmin = xx[0], xmax = np.max(xx))
+          plt.axes().set_ylim(ymin=-0.1, ymax=1.1)
+          #plt.axes().set_aspect('equal')
+          plt.xlabel('x')
+          plt.ylabel('p')
+          #plt.tight_layout()
+          plt.show(block=False)
+          plt.pause(0.00001)
+
+        return None

examples/acoustic_1d_imex/ProblemClass.py

@@ -17,14 +17,10 @@
 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(2.0*np.pi*x)
+    return np.sin(x)
 #    return np.exp(-0.5*(x-0.5)**2/0.1**2)
 
 class acoustic_1d_imex(ptype):
@@ -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, 1.0, 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, 1.0, 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
 
@@ -194,9 +160,16 @@ def u_exact(self,t):
             exact solution
         """
 
+        sigma_0 = 0.1
+        k       = 7.2*np.pi
+        x_0     = 0.75
+        x_1     = 0.25
+        
         me             = mesh(self.nvars)
-        me.values[0,:] = 0.5*u_initial(self.mesh - (self.cadv + self.cs)*t) + 0.5*u_initial(self.mesh - (self.cadv - self.cs)*t)
-        me.values[1,:] = 0.5*u_initial(self.mesh - (self.cadv + self.cs)*t) - 0.5*u_initial(self.mesh - (self.cadv - self.cs)*t)
+        #me.values[0,:] = 0.5*u_initial(self.mesh - (self.cadv + self.cs)*t) + 0.5*u_initial(self.mesh - (self.cadv - self.cs)*t)
+        #me.values[1,:] = 0.5*u_initial(self.mesh - (self.cadv + self.cs)*t) - 0.5*u_initial(self.mesh - (self.cadv - self.cs)*t)
+        me.values[0,:] = np.exp(-np.square(self.mesh-x_0-self.cs*t)/(sigma_0*sigma_0)) + np.exp(-np.square(self.mesh-x_1-self.cs*t)/(sigma_0*sigma_0))*np.cos(k*(self.mesh-self.cs*t)/sigma_0)
+        me.values[1,:] = me.values[0,:]
         return me
 
 

examples/acoustic_1d_imex/buildFDMatrix.py

@@ -0,0 +1,120 @@
+import numpy as np
+import scipy.sparse as sp
+import scipy.linalg as la
+
+# Only for periodic BC because we have advection only in x direction
+def getHorizontalDx(N, dx, order):
+  
+  if order==1:
+    stencil    = [-1.0, 1.0]
+    zero_pos   = 2
+    coeff      = 1.0
+  
+  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:
+    print "Order "+str(order)+" not implemented."
+
+  first_col = np.zeros(N)
+  
+  # Because we need to specific first column (not row) in circulant, flip stencil array
+  first_col[0:np.size(stencil)] = np.flipud(stencil)
+
+  # Circulant shift of coefficient column so that entry number zero_pos becomes first entry
+  first_col = np.roll(first_col, -np.size(stencil)+zero_pos, axis=0)
+
+  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)
+
+  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"
+  
+  if bc_left in ['periodic']:
+    A[0,N-2] = stencil[0]
+    A[0,N-1] = stencil[1]
+    A[1,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]
+
+  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 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 bc_left in ['dirichlet']:
+    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
+
+  return sp.csc_matrix(A)
+
+def getBCLeft(value, N, dx, type):
+
+  assert type in ['periodic','neumann','dirichlet'], "Unknown type of BC"
+
+  b = np.zeros(N)
+  if type in ['dirichlet']:
+    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
+
+  return (1.0/(12.0*dx))*b
+
+def getBCRight(value, N, dx, type):
+
+  assert type in ['periodic','neumann','dirichlet'], "Unknown type of BC"
+
+  b = np.zeros(N)
+  if type in ['dirichlet']:
+    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
+
+  return (1.0/(12.0*dx))*b

examples/acoustic_1d_imex/buildWave1DMatrix.py

@@ -0,0 +1,35 @@
+import sys
+sys.path.append('../')
+import numpy as np
+import scipy.sparse as sp
+from buildFDMatrix import getMatrix, getHorizontalDx, getBCLeft, getBCRight
+
+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])
+  Zero = np.zeros((N,N))
+  M1 = sp.hstack((Zero, D_p), format="csc")
+  M2 = sp.hstack((D_u, Zero), format="csc")
+  M   = sp.vstack((M1, M2),  format="csc")
+  return sp.csc_matrix(Id), sp.csc_matrix(M)
+
+def getWave1DAdvectionMatrix(N, dx, order):
+  Dx   = getHorizontalDx(N, dx, order)
+  Zero = np.zeros((N,N))
+  M1   = sp.hstack((Dx, Zero), format="csc")
+  M2   = sp.hstack((Zero, Dx), format="csc")
+  M    = sp.vstack((M1, M2), format="csc")
+  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])
+  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])
+  return np.concatenate((bp, bu))