Merge pull request #34 from danielru/setup/fwsw_acoustic1d

Setup/fwsw acoustic1d

Author: pancetta

examples/acoustic_1d_imex/HookClass.py

@@ -20,7 +20,7 @@ def __init__(self):
 
         # add figure object for further use
         #        self.fig = plt.figure(figsize=(18,6))
-        self.fig = plt.figure(figsize=(9,9))
+        #self.fig = plt.figure(figsize=(9,9))
 
 
     def dump_step(self,status):
@@ -31,17 +31,17 @@ def dump_step(self,status):
             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=-1.0, ymax=1.0)
-          plt.axes().set_aspect('equal')
+          plt.axes().set_ylim(ymin=-0.1, ymax=1.1)
+          #plt.axes().set_aspect('equal')
           plt.xlabel('x')
-          plt.ylabel('u')
+          plt.ylabel('p')
           #plt.tight_layout()
           plt.show(block=False)
           plt.pause(0.00001)

examples/acoustic_1d_imex/ProblemClass.py

@@ -56,12 +56,12 @@ 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)
 
-        self.mesh   = np.linspace(0.0, 2.0*np.pi, self.nvars[1], endpoint=False)
+        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
+        self.A      = -self.cs*A
 
     def solve_system(self,rhs,factor,u0,t):
         """
@@ -160,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/playground.py

@@ -30,19 +30,19 @@
     lparams['restol'] = 1E-10
 
     sparams = {}
-    sparams['maxiter'] = 20
+    sparams['maxiter'] = 2
 
     # setup parameters "in time"
-    t0 = 0
-    dt = 0.01
-    Tend = 100*dt
-
+    t0 = 0.0
+    Tend = 3.0
+    dt = Tend/float(154)
+    
     # This comes as read-in for the problem class
     pparams = {}
-    pparams['nvars'] = [(2,500)]
-    pparams['cadv']  = 0.1
+    pparams['nvars'] = [(2,512)]
+    pparams['cadv']  = 0.0
     pparams['cs']    = 1.0
-    pparams['order_adv'] = 4
+    pparams['order_adv'] = 5
 
     # This comes as read-in for the transfer operations
     tparams = {}
@@ -55,7 +55,7 @@
     description['dtype_u']           = mesh
     description['dtype_f']           = rhs_imex_mesh
     description['collocation_class'] = collclass.CollGaussLobatto
-    description['num_nodes']         = 4
+    description['num_nodes']         = 2
     description['sweeper_class']     = imex_1st_order
     description['level_params']      = lparams
     description['hook_class']        = plot_solution
@@ -79,11 +79,17 @@
 
     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.legend()
+    sigma_0 = 0.1
+    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, uex.values[1,:],  '+', color='r', label='p (exact)')
+    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.legend(loc=2)
     plt.xlim([0, 1])
-    plt.ylim([-1, 1])
-    plt.show()
+    plt.ylim([-0.1, 1.1])
+    #plt.show()
+    plt.gcf().savefig('fwsw-sdc-K'+str(sparams['maxiter'])+'-M'+str(description['num_nodes'])+'.pdf', bbox_inches='tight')