1 import FiniteVolumes2DDiffusion_SQUARE
2 import matplotlib.pyplot as plt
4 from math import log10, sqrt
7 convergence_synthesis=dict(FiniteVolumes2DDiffusion_SQUARE.test_desc)
9 def test_validation2DVF_cross_triangles():
11 ### 2D FV cross triangles mesh
12 #meshList=[5,9,15,21,31]
13 meshList=['squareWithCrossTriangles_0','squareWithCrossTriangles_1','squareWithCrossTriangles_2','squareWithCrossTriangles_3']#,'squareWithCrossTriangles_4'
14 mesh_path='../../../ressources/2DCrossTriangles/'
15 meshType="Regular_cross_triangles"
17 nbMeshes=len(meshList)
18 error_tab=[0]*nbMeshes
19 mesh_size_tab=[0]*nbMeshes
20 mesh_name='squareWithCrossTriangles'
21 diag_data=[0]*nbMeshes
24 curv_abs=np.linspace(0,sqrt(2),resolution+1)
27 # Storing of numerical errors, mesh sizes and diagonal values
28 for filename in meshList:
30 #my_mesh=cdmath.Mesh(0,1,nx,0,1,nx*nx)#Use script provided by Adrien to split each quadrangle in 4 triangles
31 #error_tab[i], mesh_size_tab[i], diag_data[i], min_sol_num, max_sol_num, time_tab[i] =FiniteVolumes2DDiffusion_SQUARE.solve(my_mesh,str(nx)+'x'+str(nx),resolution,meshType,testColor)
32 error_tab[i], mesh_size_tab[i], diag_data[i], min_sol_num, max_sol_num, time_tab[i] =FiniteVolumes2DDiffusion_SQUARE.solve_file(mesh_path+filename,resolution,meshType,testColor)
33 assert min_sol_num>-0.01
34 assert max_sol_num<1.53
35 plt.plot(curv_abs, diag_data[i], label= str(mesh_size_tab[i]) + ' cells')
36 error_tab[i]=log10(error_tab[i])
37 time_tab[i]=log10(time_tab[i])
38 mesh_size_tab[i] = 0.5*log10(mesh_size_tab[i])
43 # Plot over diagonal line
45 plt.xlabel('Position on diagonal line')
46 plt.ylabel('Value on diagonal line')
47 plt.title('Plot over diagonal line for finite volumes \n for the diffusion equation on 2D cross triangles meshes')
48 plt.savefig(mesh_name+"_2DDiffusionFV_PlotOverDiagonalLine.png", dpi='figure')
50 # Least square linear regression
51 # Find the best a,b such that f(x)=ax+b best approximates the convergence curve
52 # The vector X=(a,b) solves a symmetric linear system AX=B with A=(a1,a2\\a2,a3), B=(b1,b2)
53 a1=np.dot(mesh_size_tab,mesh_size_tab)
54 a2=np.sum(mesh_size_tab)
56 b1=np.dot(error_tab,mesh_size_tab)
60 assert det!=0, 'test_validation2DVF_cross_triangles() : Make sure you use distinct meshes and at least two meshes'
64 print "FV for diffusion on 2D cross triangles meshes : scheme order is ", -a
65 assert abs(a+0.035)<0.001
67 # Plot of convergence curve
69 plt.plot(mesh_size_tab, error_tab, label='log(|numerical-exact|)')
70 plt.plot(mesh_size_tab, a*np.array(mesh_size_tab)+b,label='least square slope : '+'%.3f' % a)
72 plt.plot(mesh_size_tab, error_tab)
73 plt.xlabel('log(sqrt(number of cells))')
74 plt.ylabel('log(error)')
75 plt.title('Convergence of finite volumes for \n the diffusion equation on 2D cross triangles meshes')
76 plt.savefig(mesh_name+"_2DDiffusionFV_ConvergenceCurve.png", dpi='figure')
78 # Plot of computational time
80 plt.plot(mesh_size_tab, time_tab, label='log(cpu time)')
82 plt.xlabel('log(sqrt(number of cells))')
83 plt.ylabel('log(cpu time)')
84 plt.title('Computational time of finite volumes \n for the diffusion equation on 2D cross triangles meshes')
85 plt.savefig(mesh_name+"_2DDiffusionFV_ComputationalTime.png", dpi='figure')
89 convergence_synthesis["Mesh_names"]=meshList
90 convergence_synthesis["Mesh_type"]=meshType
91 convergence_synthesis["Mesh_path"]=mesh_path
92 convergence_synthesis["Mesh_description"]=mesh_name
93 convergence_synthesis["Mesh_sizes"]=[10**x for x in mesh_size_tab]
94 convergence_synthesis["Space_dimension"]=2
95 convergence_synthesis["Mesh_dimension"]=2
96 convergence_synthesis["Mesh_cell_type"]="cross triangles"
97 convergence_synthesis["Errors"]=[10**x for x in error_tab]
98 convergence_synthesis["Scheme_order"]=-a
99 convergence_synthesis["Test_color"]=testColor
100 convergence_synthesis["Computational_time"]=end-start
102 with open('Convergence_Diffusion_2DVF_'+mesh_name+'.json', 'w') as outfile:
103 json.dump(convergence_synthesis, outfile)
105 if __name__ == """__main__""":
106 test_validation2DVF_cross_triangles()