update

onedimtest.py

@@ -45,7 +45,7 @@ def F(Q, DQ, d):
     return F_
 
 
-nx = 100
+nx = 5000
 L = [1.]
 tf = .5
 
@@ -76,18 +76,18 @@ from pypde import pde_solver
 #initial states
 plt.plot(Q0)
 plt.gca().legend(('rho','rho1','vx'))
-plt.ylim((-.8,2))
+plt.ylim((-.8,2.5))
 plt.savefig('{}/rho{:0>4d}'.format('results',0) )
 plt.close()
 
 '''SOLVER'''
-out = pde_solver(Q0, tf, L, F=F, flux='rusanov', order=3,boundaryTypes='periodic',ndt=100)
+out = pde_solver(Q0, tf, L, F=F, flux='rusanov', order=3,boundaryTypes='periodic',ndt=100,stiff=False)
 
 
 for n in range(out.shape[0]):
   if not out[n, :, 0].all() < 1e-10:
     plt.plot(out[n, :, :])
     plt.gca().legend(('rho','rho1','vx'))
-    plt.ylim((-.8,2))
+    plt.ylim((-.8,2.5))
     plt.savefig('{}/rho{:0>4d}'.format('results',n+1) )
     plt.close()

twodimtest.py

@@ -53,8 +53,8 @@ def F(Q, DQ, d):
     F_[0] = u[d]
     F_[1] = r1 * u[d] / r
 
-    F_[1] += - tau * dr1[d] + tau * r1 / r * dr[d]  # this diffusion makes it difficult
-    F_[1] +=  np.array([-1,1])[d]*( -tau * dr1[d2] + tau * r1 / r * dr[d2] )*a
+    #F_[1] += - tau * dr1[d] + tau * r1 / r * dr[d]  # this diffusion makes it difficult
+    #F_[1] +=  np.array([-1,1])[d]*( -tau * dr1[d2] + tau * r1 / r * dr[d2] )*a
 
 
     F_[2:4] = u * u[d] / r   # u u[d] = u_i u_j
@@ -63,11 +63,14 @@ def F(Q, DQ, d):
     return F_
 
 
-nx,ny = 100, 100
+nx,ny = 200, 200
 
 
 L = [2 * pi, 2 * pi]
-tf = 1.5
+
+L = [1., 1.]
+
+tf = 1.
 
 
 
@@ -96,8 +99,12 @@ cx,cy = np.array(L)*0.5
 
 for i in range(nx):
   for j in range(ny):
-    x = (i + 0.5) * L[0] / nx
-    y = (j + 0.5) * L[1] / ny
+    #x = (i + 0.5) * L[0] / nx
+    #y = (j + 0.5) * L[1] / ny
+
+    x = (i ) * L[0] / ( nx -1 )
+    y = (j ) * L[1] / ( ny -1 )
+
 
     if 0.2*L[0] < np.sqrt( (x-cx)**2+(y-cy)**2 )  < 0.4*L[1]:
         v = np.array([-(x-cx),-(y-cy)])/np.sqrt((x-cx)**2+(y-cy)**2)*0.1
@@ -105,7 +112,7 @@ for i in range(nx):
         v = np.array([0.,0.])
 
     Q0[i, j, 0] = 1.0 #- 0.1 * cos(2*y) * cos(2*x)
-    Q0[i, j, 1] = 0.6 #0.1 * np.sin(2*pi*y) * np.sin(2*pi*x) + 0.5
+    Q0[i, j, 1] = 0.4 #0.1 * np.sin(2*pi*y) * np.sin(2*pi*x) + 0.5
     Q0[i, j, 2] = v[0]
     Q0[i, j, 3] = v[1]
 #'''
@@ -116,7 +123,8 @@ outx[0] = Q0
 savevtk2(0,outx,L,'results',0)
 
 '''SOLVER'''
-out = pde_solver(Q0, tf, L, F=F, flux='rusanov', order=3,boundaryTypes='periodic', stiff=False,ndt=100)
+out = pde_solver(Q0, tf, L, F=F, flux='rusanov',
+        order=4,boundaryTypes='periodic', stiff=False,ndt=100)
 
 #out = np.vstack((out,[Q0]))#add initial condition
 for n in range(out.shape[0]):