convergence_no_an_Pgl.py 3.35 KB
 Davide Torlo committed Apr 21, 2021 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()``````