new function to run bdf2, implicit euler and trapezoidal rule for comparison; fixed wave length in fwsw sdc initial data

Author: Daniel Ruprecht

examples/acoustic_1d_imex/ProblemClass.py

@@ -161,7 +161,7 @@ def u_exact(self,t):
         """
 
         sigma_0 = 0.1
-        k       = 7.2*np.pi
+        k       = 7.0*2.0*np.pi
         x_0     = 0.75
         x_1     = 0.25
 

examples/acoustic_1d_imex/playground.py

@@ -13,8 +13,8 @@
 from pySDC.Stats import grep_stats, sort_stats
 
 # Sharpclaw imports
-from clawpack import pyclaw
-from clawpack import riemann
+#from clawpack import pyclaw
+#from clawpack import riemann
 from matplotlib import pyplot as plt
 
 
@@ -40,7 +40,7 @@
     # This comes as read-in for the problem class
     pparams = {}
     pparams['nvars'] = [(2,512)]
-    pparams['cadv']  = 0.0
+    pparams['cadv']  = 0.1
     pparams['cs']    = 1.0
     pparams['order_adv'] = 5
 
@@ -83,13 +83,14 @@
     x_0     = 0.75
 
     #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, uend.values[1,:], '-', color='b', label='SDC')
     #plt.plot(P.mesh, uex.values[1,:],  '+', color='r', label='p (exact)')
-    plt.plot(P.mesh, uend.values[1,:], '-', color='b', linewidth=2.0, label='p (SDC)')
+    #plt.plot(P.mesh, uend.values[1,:], '-', color='b', linewidth=2.0, label='p (SDC)')
     p_slow = np.exp(-np.square(P.mesh-x_0)/(sigma_0*sigma_0))
-    plt.plot(P.mesh, p_slow, '+', color='r', markersize=4, label='slow mode')
+    #plt.plot(P.mesh, p_slow, '-', color='r', markersize=4, label='slow mode')
     plt.legend(loc=2)
     plt.xlim([0, 1])
     plt.ylim([-0.1, 1.1])
+    fig.gca().grid()
     #plt.show()
     plt.gcf().savefig('fwsw-sdc-K'+str(sparams['maxiter'])+'-M'+str(description['num_nodes'])+'.pdf', bbox_inches='tight')

examples/acoustic_1d_imex/run_standard.py

@@ -0,0 +1,87 @@
+import numpy as np
+import scipy.sparse as sp
+import scipy.sparse.linalg as LA
+from matplotlib import pyplot as plt
+from buildWave1DMatrix import getWave1DMatrix, getWave1DAdvectionMatrix
+
+sigma_0 = 0.1
+k       = 7.0*2*np.pi
+x_0     = 0.75
+x_1     = 0.25
+multi_scale = 1.0
+
+nvars = 512
+cs    = 1.0
+cadv  = 0.0
+order = 4
+
+def u(x,t):
+  u0 = np.exp(-np.square( np.mod( mesh- cs*t, 1.0 ) -x_0 )/(sigma_0*sigma_0)) + multi_scale*np.exp(-np.square( np.mod( mesh -cs*t, 1.0 ) -x_1 )/(sigma_0*sigma_0))*np.cos(k*( np.mod( mesh-cs*t, 1.0 ))/sigma_0)
+  p0 = u0
+  return u0, p0
+
+Tend   = 3.0
+Nsteps = 154
+dt = Tend/float(Nsteps)
+
+mesh   = np.linspace(0.0, 1.0, nvars, endpoint=False)
+dx     = mesh[1] - mesh[0]
+Dx     = -cadv*getWave1DAdvectionMatrix(nvars, dx, order)
+Id, A  = getWave1DMatrix(nvars, dx, ['periodic','periodic'], ['periodic','periodic'])
+A      = -cs*A
+
+M_ieuler = Id - dt*(A + Dx)
+
+M_bdf    = Id - (2.0/3.0)*dt*(A + Dx)
+
+alpha    = 0.5
+M_trap   = Id - alpha*dt*(A+Dx)
+B_trap   = Id + (1-alpha)*dt*(A+Dx)
+ 
+u0, p0 = u(mesh, 0.0)
+y0_ie  = np.concatenate( (u0, p0) )
+y0_tp  = y0_ie
+
+y0_bdf = y0_ie
+
+fig = plt.figure(figsize=(8,8))
+
+for i in range(0,Nsteps):
+
+  # implicit Euler step
+  ynew_ie = LA.spsolve(M_ieuler, y0_ie)
+
+  # trapezoidal rule step
+  b_tp    = B_trap.dot(y0_tp)
+  ynew_tp = LA.spsolve(M_trap, b_tp)
+
+  # BDF-2 scheme
+  if i==0:
+    ynew_bdf = LA.spsolve(M_ieuler, y0_bdf)
+  else:
+    b_bdf    = (4.0/3.0)*y0_bdf - (1.0/3.0)*ym1_bdf
+    ynew_bdf = LA.spsolve(M_bdf, b_bdf)
+
+  unew_ie, pnew_ie = np.split(ynew_ie, 2)
+  unew_tp, pnew_tp = np.split(ynew_tp, 2)
+  unew_bdf, pnew_bdf = np.split(ynew_bdf, 2)
+  uex, pex = u(mesh, float(i+1)*dt)
+
+  fig.gca().clear()
+  #plt.plot(mesh, pnew_bdf, 'b', label='BDF-2')
+  plt.plot(mesh, pnew_tp, 'r', label='Trapezoidal')
+  plt.plot(mesh, pex, 'k', label='Exact')
+  fig.gca().set_xlim([0, 1.0])
+  fig.gca().set_ylim([-0.5, 1.1])
+  fig.gca().legend(loc=3)
+  fig.gca().grid()
+  plt.draw()
+  plt.pause(0.0001)  
+  #if i==0:
+  #  plt.gcf().savefig('initial.pdf', bbox_inches='tight')
+  y0_ie   = ynew_ie
+  y0_tp   = ynew_tp
+  ym1_bdf = y0_bdf
+  y0_bdf  = ynew_bdf
+#plt.show()
+#plt.gcf().savefig('final.pdf', bbox_inches='tight')