convergence_no_an_Pgl.py 3.35 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
#!/usr/bin/python
import numpy as np 
from numpy import linalg as LA
import scipy
import matplotlib.pyplot as pl
import sys 
import itertools
#import exact_sol as ex

# run the convergence without an exact solution
# files saved in columns x, h, u
# ns and the file are below and must be changed accordingly

def remove_extra(x, u):
   ind=[]
   for i in range(len(x)-1):
      if x[i]==x[i+1]:
          ind.append(i)
   x=np.delete(x,ind,0)
   u=np.delete(u,ind,0)
   return x,u

pl.rc("text", usetex=True)
pl.rc("font", family="serif")

# Grids for convergence study
#nNodes = [5,10,20,40,80,160,320,640,1280,2560]
#nNodes = [5,10,20,40,80,160,320,640,1280]
#nNodes = [5,10,20,40,80,160,320,640]
#nNodes = [5,10,20,40,80,160,320]
#nNodes = [5,10,20,40,80,160]
#nNodes = [5,10,20,40,80]
#nNodes = [5,10,20,40]
nNodes = [8,16,32,64,128]#,256,512, 1024]#,2048,4096]
#nNodes = [25,50,100,200]#,400,800]
# Test name
test = "test"

# scheme
scheme = "scheme4"

# Scheme order
orders = [1, 2,3]#,4]

marker = itertools.cycle(('s','o','*'))

varnames = ['h','u']

for nvar, var in enumerate(varnames):
		 # Initialize plots
		 fig_err_L1 = pl.figure()
		 ax_err_L1  = fig_err_L1.add_subplot(1,1,1)
		
		 # Initialize array to store errors
		 err_L1 = np.zeros((np.size(orders),np.size(nNodes)-1))
		 dx     = np.zeros((np.size(orders),np.size(nNodes)-1))

		 # Initialize array to store orders
		 order_L1 = np.zeros((np.size(orders),np.size(nNodes)-1))
		
		 for k, order in enumerate(orders):
		    xs =[]
		    num_sols = []
		    for i, N in enumerate(nNodes):

		       # Load numerical solution
		       fname = "Pgl"+str(order)+"/"+str(test)+"_"+scheme+"_"+str(N)+".dat"  #here the folder and filenames
		       x, num_sol = np.loadtxt(fname,delimiter=None,usecols=(0,nvar+1), unpack=True)
		       x, num_sol = remove_extra(x,num_sol)
		       xs.append(x)
		       num_sols.append(num_sol)

		    for i in range(len(nNodes)-1):
		       err_L1[k,i]=LA.norm(num_sols[i][::order]-num_sols[i+1][::2*order],ord=2)/float(np.sqrt(nNodes[i]))

		       # Compute exact solution
	#         ex_sol = ex.exact_sw_smooth(x)
	 
		       # Compute errors
		       #err_L1[k,i] = LA.norm(num_sol-ex_sol[nvar],ord=np.inf)
	#         err_L1[k,i] = LA.norm(num_sol-ex_sol[nvar],ord=1)/float(N)
		       #err_L1[k,i] = LA.norm(num_sol-ex_sol[nvar],ord=2)/float(np.sqrt(N))
		       #err_L1[k,i] = LA.norm((num_sol-ex_sol[nvar])*dxvec,ord=1)
		       
		       if (i > 0):
		          order_L1[k,i] = -( np.log(err_L1[k][i])-np.log(err_L1[k][i-1]) ) / ( np.log(nNodes[i])-np.log(nNodes[i-1]) )
		    
		    print("Variable : ", var)
		    print("L1-errors for order ", order, " :")
		    print("Computed errors : ")
		    print(err_L1[k,:])
		    print("Computed orders : ")
		    print(order_L1[k,:])

		    # Plot figures

		    # error vs mesh size
		    ax_err_L1.loglog(nNodes[1:],err_L1[k,:],"-"+marker.next(),label="Pgl"+str(order))
		    ref_err = [err_L1[k,1]*nNodes[1]**(float(order+1))*nNodes[i]**(-float(order+1)) for i in range(np.size(nNodes[1:]))]
		    ax_err_L1.loglog(nNodes[1:],ref_err,"--",label="order "+str(order+1))
		    ax_err_L1.set_xlabel(r"$N$")
		    ax_err_L1.set_ylabel(r"$L_2$-error")
		    pl.xlim(nNodes[1]*0.7,nNodes[-1]*1.3)
		    ax_err_L1.legend(loc=0)
		    pl.title("Convergence of "+var)      

     # Save plots
		 fig_err_L1.savefig(scheme+"_"+var+"_convergence_no_an_Pgl.pdf",format="pdf")

pl.show()