Merge pull request #48 from danielru/fwsw_sdc_acoustic1d

Iteration and convergence plots for fwsw SDC

Author: pancetta

examples/acoustic_1d_imex/HookClass.py

@@ -22,6 +22,8 @@ def __init__(self):
         #        self.fig = plt.figure(figsize=(18,6))
         #self.fig = plt.figure(figsize=(9,9))
 
+    def dump_sweep(self,status):
+        return None
 
     def dump_step(self,status):
         """

examples/acoustic_1d_imex/ProblemClass_conv.py

@@ -0,0 +1,168 @@
+r"""
+  One-dimensional IMEX acoustic-advection
+  =========================
+  
+  Integrate the linear 1D acoustic-advection problem:
+  
+  .. math::
+  u_t + U u_x + c p_x & = 0 \\
+  p_t + U p_x + c u_x & = 0.
+  
+"""
+
+import numpy as np
+import scipy.sparse as sp
+import scipy.sparse.linalg as LA
+
+from pySDC.Problem import ptype
+from pySDC.datatype_classes.mesh import mesh, rhs_imex_mesh
+
+from buildWave1DMatrix import getWave1DMatrix, getWave1DAdvectionMatrix
+
+def u_initial(x, k):
+    return np.sin(k*2.0*np.pi*x)
+
+class acoustic_1d_imex(ptype):
+    """
+    Example implementing the forced 1D heat equation with Dirichlet-0 BC in [0,1]
+
+    Attributes:
+      solver: Sharpclaw solver
+      state:  Sharclaw state
+      domain: Sharpclaw domain
+    """
+
+    def __init__(self, cparams, dtype_u, dtype_f):
+        """
+        Initialization routine
+
+        Args:
+            cparams: custom parameters for the example
+            dtype_u: particle data type (will be passed parent class)
+            dtype_f: acceleration data type (will be passed parent class)
+        """
+
+        # these parameters will be used later, so assert their existence
+        assert 'nvars' in cparams
+        assert 'cs' in cparams
+        assert 'cadv' in cparams
+        assert 'order_adv' in cparams
+        assert 'waveno' in cparams
+        
+        # add parameters as attributes for further reference
+        for k,v in cparams.items():
+            setattr(self,k,v)
+
+        # invoke super init, passing number of dofs, dtype_u and dtype_f
+        super(acoustic_1d_imex,self).__init__(self.nvars,dtype_u,dtype_f)
+        
+        self.mesh   = np.linspace(0.0, 1.0, self.nvars[1], endpoint=False)
+        self.dx     = self.mesh[1] - self.mesh[0]
+        
+        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
+
+        Args:
+            rhs: right-hand side for the nonlinear system
+            factor: abbrev. for the node-to-node stepsize (or any other factor required)
+            u0: initial guess for the iterative solver (not used here so far)
+            t: current time (e.g. for time-dependent BCs)
+
+        Returns:
+            solution as mesh
+        """
+        
+        M  = self.Id - factor*self.A
+        
+        b = np.concatenate( (rhs.values[0,:], rhs.values[1,:]) )
+        
+        sol = LA.spsolve(M, b)
+
+        me = mesh(self.nvars)
+        me.values[0,:], me.values[1,:] = np.split(sol, 2)
+        
+        return me
+
+
+    def __eval_fexpl(self,u,t):
+        """
+        Helper routine to evaluate the explicit part of the RHS
+
+        Args:
+            u: current values (not used here)
+            t: current time
+
+        Returns:
+            explicit part of RHS
+        """
+
+
+        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)
+
+        return fexpl
+
+
+    def __eval_fimpl(self,u,t):
+        """
+        Helper routine to evaluate the implicit part of the RHS
+
+        Args:
+            u: current values
+            t: current time (not used here)
+
+        Returns:
+            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,:], fimpl.values[1,:] = np.split(sol, 2)
+        
+        return fimpl
+
+
+    def eval_f(self,u,t):
+        """
+        Routine to evaluate both parts of the RHS
+
+        Args:
+            u: current values
+            t: current time
+
+        Returns:
+            the RHS divided into two parts
+        """
+
+        f = rhs_imex_mesh(self.nvars)
+        f.impl = self.__eval_fimpl(u,t)
+        f.expl = self.__eval_fexpl(u,t)
+        return f
+
+    def u_exact(self,t):
+        """
+        Routine to compute the exact solution at time t
+
+        Args:
+            t: current time
+
+        Returns:
+            exact solution
+        """
+        
+        me             = mesh(self.nvars)
+        me.values[0,:] = 0.5*u_initial( self.mesh - (self.cadv + self.cs)*t , self.waveno) - 0.5*u_initial( self.mesh - (self.cadv - self.cs)*t , self.waveno)
+        me.values[1,:] = 0.5*u_initial( self.mesh - (self.cadv + self.cs)*t , self.waveno) + 0.5*u_initial( self.mesh - (self.cadv - self.cs)*t , self.waveno)
+        return me
+
+

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