SALOME platform Git repositories - tools/solverlab.git/blob

   1 import FiniteVolumes2DPoisson_SQUARE
   2 import matplotlib.pyplot as plt
   3 import numpy as np
   4 from math import log10, sqrt
   5 import time, json
   6
   7 convergence_synthesis=dict(FiniteVolumes2DPoisson_SQUARE.test_desc)
   8
   9 def test_validation2DVF_flat_cross_triangles():
  10     start = time.time()
  11     ### 2D FV flat cross triangles mesh
  12     #meshList=[5,9,15,21,31]
  13     meshList=['squareWithFlatCrossTriangles_0','squareWithFlatCrossTriangles_1','squareWithFlatCrossTriangles_2','squareWithFlatCrossTriangles_3']#,'squareWithFlatCrossTriangles_4'
  14     mesh_path='../../../ressources/2DFlatCrossTriangles/'
  15     meshType="Regular_flat_cross_triangles"
  16     testColor="Green"
  17     nbMeshes=len(meshList)
  18     error_tab=[0]*nbMeshes
  19     mesh_size_tab=[0]*nbMeshes
  20     mesh_name='squareWithFlatCrossTriangles'
  21     diag_data=[0]*nbMeshes
  22     time_tab=[0]*nbMeshes
  23     resolution=100
  24     curv_abs=np.linspace(0,sqrt(2),resolution+1)
  25     plt.close('all')
  26     i=0
  27     # Storing of numerical errors, mesh sizes and diagonal values
  28     for filename in meshList:
  29     #for nx 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] =FiniteVolumes2DPoisson_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] =FiniteVolumes2DPoisson_SQUARE.solve_file(mesh_path+filename,resolution,meshType,testColor)
  33         assert min_sol_num>-0.01
  34         assert max_sol_num<1.4
  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])
  39         i=i+1
  40
  41     end = time.time()
  42
  43     # Plot over diagonal line
  44     plt.legend()
  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 Laplace operator on 2D flat cross triangles meshes')
  48     plt.savefig(mesh_name+"_2DPoissonFV_PlotOverDiagonalLine.png")
  49
  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)
  55     a3=nbMeshes
  56     b1=np.dot(error_tab,mesh_size_tab)
  57     b2=np.sum(error_tab)
  58
  59     det=a1*a3-a2*a2
  60     assert det!=0, 'test_validation2DVF_flat_triangles() : Make sure you use distinct meshes and at least two meshes'
  61     a=( a3*b1-a2*b2)/det
  62     b=(-a2*b1+a1*b2)/det
  63
  64     print( "FV on 2D flat cross triangles mesh : scheme order is ", -a )
  65     assert abs(a-0.12)<0.1
  66
  67     # Plot of convergence curve
  68     plt.close()
  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)
  71     plt.legend()
  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 Laplace operator on 2D flat cross triangles meshes')
  76     plt.savefig(mesh_name+"_2DPoissonFV_ConvergenceCurve.png")
  77
  78     # Plot of computational time
  79     plt.close()
  80     plt.plot(mesh_size_tab, time_tab, label='log(cpu time)')
  81     plt.legend()
  82     plt.xlabel('log(sqrt(number of cells))')
  83     plt.ylabel('log(cpu time)')
  84     plt.title('Computational time of finite volumes \n for Laplace operator on 2D flat cross triangles meshes')
  85     plt.savefig(mesh_name+"_2DPoissonFV_ComputationalTime.png")
  86
  87     plt.close('all')
  88
  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"]="flat 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
 101
 102     with open('Convergence_Poisson_2DVF_'+mesh_name+'.json', 'w') as outfile:
 103         json.dump(convergence_synthesis, outfile)
 104
 105     import os
 106     os.system("jupyter-nbconvert --to notebook --execute Convergence_Poisson_FV5_SQUARE_flat_triangles.ipynb")
 107     os.system("jupyter-nbconvert --to html Convergence_Poisson_FV5_SQUARE_flat_triangles.ipynb")
 108     os.system("jupyter-nbconvert --to pdf Convergence_Poisson_FV5_SQUARE_flat_triangles.ipynb")
 109
 110 if __name__ == """__main__""":
 111     test_validation2DVF_flat_cross_triangles()