converge.py

   1 import math
   2 import numpy as np
   3 import matplotlib as mpl
   4 import matplotlib.pyplot as plt
   5 import matplotlib.mlab as mlab
   6 from matplotlib.pyplot import *
   7 from pylab import *
   8
   9 DATA_RK = "../rk_conv.dat"
  10 DATA_EU = "../eu_conv.dat"
  11 FOUT = "../converge.png"
  12
  13 xlabel = r'$h$'
  14 ylabel = r'$\delta y$'
  15 res = (1440/80,900/80) # default dpi is 80
  16
  17 ######################################################################
  18
  19 def mean_var(param,var):
  20         unique_param = unique(param)
  21         list_var = [ [] for DUMMYVAR in range(len(unique_param)) ]
  22         ret_avg = []
  23         for i in range(len(unique_param)):
  24                 print "now:", unique_param[i]
  25                 for j in range(len(var)):
  26                         if (param[j] == unique_param[i]):
  27                                 list_var[i].append(var[j])
  28                                 #assert len(list_var[i]) == len(param/unique_param)
  29                 ret_avg.append(float(sum(list_var[i])) /
  30                                 len(list_var[i]))
  31                 #print list_var[i], ret_avg[i]
  32                 #assert list_var[i] == ret_avg[i]
  33         return ret_avg
  34
  35 print "reading data..."
  36
  37 rk = mlab.csv2rec(DATA_RK, delimiter='\t', comments='#')
  38 eu = mlab.csv2rec(DATA_EU, delimiter='\t', comments='#')
  39
  40 h_rk = rk.r
  41 h_eu = eu.r
  42 x_rk = rk.x
  43 x_eu = eu.x
  44
  45 h_min = h_rk.min()
  46 h_max = h_rk.max()
  47
  48 print "calculating mean..."
  49
  50 rk_mean = mean_var(h_rk, x_rk)
  51 eu_mean = mean_var(h_eu, x_eu)
  52
  53 #print rk_mean
  54
  55
  56 fig = plt.figure(figsize=res)
  57
  58 #ax1 = fig.add_subplot(211)
  59 ax2 = fig.add_subplot(111)
  60
  61 #ax1.plot(unique(h_rk), rk_mean, marker=',', ls='')
  62 #ax1.plot(unique(h_eu), eu_mean, marker=',', ls='')
  63
  64 #ax1.axhline(x_rk[0], color='r')
  65 #ax1.axhline(x_eu[0], color='r')
  66
  67 print "calculating error..."
  68
  69 err_rk = []
  70 err_eu = []
  71
  72 for i in range(len(rk_mean)):
  73         tmp = np.abs(rk_mean[0] - rk_mean[i])
  74         err_rk.append(tmp)
  75
  76 for i in range(len(eu_mean)):
  77         tmp = np.abs(rk_mean[0] - eu_mean[i])
  78         err_eu.append(tmp)
  79
  80 print "plotting..."
  81 3
  82 ax2.loglog(unique(h_rk), err_rk, marker=',', ls='', ms=0.1)
  83 ax2.loglog(unique(h_eu), err_eu, marker=',', ls='', ms=0.1)
  84
  85 pad = h_max
  86
  87 ax2.set_xlabel(xlabel)
  88 ax2.set_ylabel(ylabel)
  89 ax2.set_xlim(h_min,h_max)
  90 ax2.set_ylim(min(err_eu),max(err_eu))
  91
  92 #ax1.set_xlim()
  93 #ax1.set_ylim(-3.3,-3.2)
  94 plt.savefig(FOUT)