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

   1 import WaveSystemUpwind
   2 import matplotlib.pyplot as plt
   3 import numpy as np
   4 from math import log10, sqrt
   5 import sys
   6 import time, json
   7
   8 def test_validation2DWaveSystemUpwindDelaunayTriangles(bctype,scaling):
   9     start = time.time()
  10     #### 2D Delaunay triangles mesh
  11     meshList=['squareWithTriangles_1','squareWithTriangles_2','squareWithTriangles_3','squareWithTriangles_4','squareWithTriangles_5']
  12     meshType="Unstructured_Delaunay_triangles"
  13     testColor="Green"
  14     nbMeshes=len(meshList)
  15     error_p_tab=[0]*nbMeshes
  16     error_u_tab=[0]*nbMeshes
  17     mesh_size_tab=[0]*nbMeshes
  18     mesh_path='../../../ressources/2DTriangles/'
  19     mesh_name='squareWithDelaunayTriangles'
  20     diag_data_press=[0]*nbMeshes
  21     diag_data_vel=[0]*nbMeshes
  22     time_tab=[0]*nbMeshes
  23     t_final=[0]*nbMeshes
  24     ndt_final=[0]*nbMeshes
  25     max_vel=[0]*nbMeshes
  26     resolution=100
  27     curv_abs=np.linspace(0,sqrt(2),resolution+1)
  28
  29     plt.close('all')
  30     i=0
  31     cfl=0.5
  32     # Storing of numerical errors, mesh sizes and diagonal values
  33     for filename in meshList:
  34         error_p_tab[i], error_u_tab[i], mesh_size_tab[i], t_final[i], ndt_final[i], max_vel[i], diag_data_press[i], diag_data_vel[i], time_tab[i] =WaveSystemUpwind.solve_file(mesh_path+filename, mesh_name, resolution,scaling,meshType,testColor,cfl,bctype)
  35         assert max_vel[i]>0.94 and max_vel[i]<1
  36         error_p_tab[i]=log10(error_p_tab[i])
  37         error_u_tab[i]=log10(error_u_tab[i])
  38         time_tab[i]=log10(time_tab[i])
  39         i=i+1
  40
  41     end = time.time()
  42
  43     # Plot over diagonal line
  44     for i in range(nbMeshes):
  45         plt.plot(curv_abs, diag_data_press[i], label= str(mesh_size_tab[i]) + ' cells')
  46     plt.legend()
  47     plt.xlabel('Position on diagonal line')
  48     plt.ylabel('Pressure on diagonal line')
  49     plt.title('Plot over diagonal line for stationary wave system \n with upwind scheme on 2D Delaunay triangles meshes')
  50     plt.savefig(mesh_name+'_Pressure_2DWaveSystemUpwind_'+"PlotOverDiagonalLine.png")
  51     plt.close()
  52
  53     plt.clf()
  54     for i in range(nbMeshes):
  55         plt.plot(curv_abs, diag_data_vel[i],   label= str(mesh_size_tab[i]) + ' cells')
  56     plt.legend()
  57     plt.xlabel('Position on diagonal line')
  58     plt.ylabel('Velocity on diagonal line')
  59     plt.title('Plot over diagonal line for the stationary wave system \n with upwind scheme on 2D Delaunay triangles meshes')
  60     plt.savefig(mesh_name+"_Velocity_2DWaveSystemUpwind_"+"PlotOverDiagonalLine.png")
  61     plt.close()
  62
  63     # Plot of number of time steps
  64     plt.close()
  65     plt.plot(mesh_size_tab, ndt_final, label='Number of time step to reach stationary regime')
  66     plt.legend()
  67     plt.xlabel('Number of cells')
  68     plt.ylabel('Max time steps for stationary regime')
  69     plt.title('Number of times steps required for the stationary Wave System \n with upwind scheme on 2D Delaunay triangles meshes')
  70     plt.savefig(mesh_name+"_2DWaveSystemUpwind_"+"TimeSteps.png")
  71
  72     # Plot of number of stationary time
  73     plt.close()
  74     plt.plot(mesh_size_tab, t_final, label='Time where stationary regime is reached')
  75     plt.legend()
  76     plt.xlabel('Number of cells')
  77     plt.ylabel('Max time for stationary regime')
  78     plt.title('Simulated time for the stationary Wave System \n with upwind scheme on 2D Delaunay triangles meshes')
  79     plt.savefig(mesh_name+"_2DWaveSystemUpwind_"+"TimeFinal.png")
  80
  81     # Plot of number of maximal velocity norm
  82     plt.close()
  83     plt.plot(mesh_size_tab, max_vel, label='Maximum velocity norm')
  84     plt.legend()
  85     plt.xlabel('Number of cells')
  86     plt.ylabel('Max velocity norm')
  87     plt.title('Maximum velocity norm for the stationary Wave System \n with upwind scheme on 2D Delaunay triangles meshes')
  88     plt.savefig(mesh_name+"_2DWaveSystemUpwind_"+"MaxVelNorm.png")
  89
  90     for i in range(nbMeshes):
  91         mesh_size_tab[i]=0.5*log10(mesh_size_tab[i])
  92
  93     # Least square linear regression
  94     # Find the best a,b such that f(x)=ax+b best approximates the convergence curve
  95     # The vector X=(a,b) solves a symmetric linear system AX=B with A=(a1,a2\\a2,a3), B=(b1,b2)
  96     a1=np.dot(mesh_size_tab,mesh_size_tab)
  97     a2=np.sum(mesh_size_tab)
  98     a3=nbMeshes
  99
 100     det=a1*a3-a2*a2
 101     assert det!=0, 'test_validation2DWaveSystemUpwindTrianglesFV() : Make sure you use distinct meshes and at least two meshes'
 102
 103     b1p=np.dot(error_p_tab,mesh_size_tab)
 104     b2p=np.sum(error_p_tab)
 105     ap=( a3*b1p-a2*b2p)/det
 106     bp=(-a2*b1p+a1*b2p)/det
 107
 108     print("FV upwind on 2D triangle meshes : scheme order for pressure is ", -ap)
 109
 110     b1u=np.dot(error_u_tab,mesh_size_tab)
 111     b2u=np.sum(error_u_tab)
 112     au=( a3*b1u-a2*b2u)/det
 113     bu=(-a2*b1u+a1*b2u)/det
 114
 115     print("FV upwind on 2D triangle meshes : scheme order for velocity is ", -au)
 116
 117     # Plot of convergence curves
 118     plt.close()
 119     plt.plot(mesh_size_tab, error_p_tab, label='|error on stationary pressure|')
 120     plt.legend()
 121     plt.xlabel('1/2 log(number of cells)')
 122     plt.ylabel('log(error p)')
 123     plt.title('Convergence of finite volumes for the stationary Wave System \n with upwind scheme on 2D Delaunay triangles meshes')
 124     plt.savefig(mesh_name+"_Pressure_2DWaveSystemUpwind_"+"ConvergenceCurve.png")
 125
 126     plt.close()
 127     plt.plot(mesh_size_tab, error_u_tab, label='|error on stationary velocity|')
 128     plt.legend()
 129     plt.xlabel('1/2 log(number of cells)')
 130     plt.ylabel('log(error u)')
 131     plt.title('Convergence of finite volumes for the stationary Wave System \n with upwind scheme on 2D Delaunay triangles meshes')
 132     plt.savefig(mesh_name+"_Velocity_2DWaveSystemUpwind_"+"ConvergenceCurve.png")
 133
 134     # Plot of computational time
 135     plt.close()
 136     plt.plot(mesh_size_tab, time_tab, label='log(cpu time)')
 137     plt.legend()
 138     plt.xlabel('1/2 log(number of cells)')
 139     plt.ylabel('log(cpu time)')
 140     plt.title('Computational time of finite volumes for the stationary Wave System \n with upwind scheme on 2D Delaunay triangles meshes')
 141     plt.savefig(mesh_name+"_2DWaveSystemUpwind_ComputationalTime.png")
 142
 143     plt.close('all')
 144
 145     convergence_synthesis={}
 146
 147     convergence_synthesis["PDE_model"]="Wave system"
 148     convergence_synthesis["PDE_is_stationary"]=False
 149     convergence_synthesis["PDE_search_for_stationary_solution"]=True
 150     convergence_synthesis["Numerical_method_name"]="Upwind"
 151     convergence_synthesis["Numerical_method_space_discretization"]="Finite volumes"
 152     convergence_synthesis["Numerical_method_time_discretization"]="Implicit"
 153     convergence_synthesis["Initial_data"]="Constant pressure, divergence free velocity"
 154     convergence_synthesis["Boundary_conditions"]="Periodic"
 155     convergence_synthesis["Numerical_parameter_cfl"]=cfl
 156     convergence_synthesis["Space_dimension"]=2
 157     convergence_synthesis["Mesh_dimension"]=2
 158     convergence_synthesis["Mesh_names"]=meshList
 159     convergence_synthesis["Mesh_type"]=meshType
 160     convergence_synthesis["Mesh_path"]=mesh_path
 161     convergence_synthesis["Mesh_description"]=mesh_name
 162     convergence_synthesis["Mesh_sizes"]=mesh_size_tab
 163     convergence_synthesis["Mesh_cell_type"]="Triangles"
 164     convergence_synthesis["Numerical_error_velocity"]=error_u_tab
 165     convergence_synthesis["Numerical_error_pressure"]=error_p_tab
 166     convergence_synthesis["Max_vel_norm"]=max_vel
 167     convergence_synthesis["Final_time"]=t_final
 168     convergence_synthesis["Final_time_step"]=ndt_final
 169     convergence_synthesis["Scheme_order"]=min(-au,-ap)
 170     convergence_synthesis["Scheme_order_vel"]=-au
 171     convergence_synthesis["Scheme_order_press"]=-ap
 172     convergence_synthesis["Scaling_preconditioner"]="None"
 173     convergence_synthesis["Test_color"]=testColor
 174     convergence_synthesis["Computational_time"]=end-start
 175
 176     with open('Convergence_WaveSystem_2DFV_Upwind_'+mesh_name+'.json', 'w') as outfile:
 177         json.dump(convergence_synthesis, outfile)
 178
 179 if __name__ == """__main__""":
 180     if len(sys.argv) >2 :
 181         bctype = sys.argv[1]
 182         scaling = int(sys.argv[2])
 183         test_validation2DWaveSystemUpwindDelaunayTriangles(bctype,scaling)
 184     else :
 185         raise ValueError("test_validation2DWaveSystemUpwindTriangles.py expects a mesh file name and a scaling parameter")