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

   1 import WaveSystemPStag
   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
   9 def test_validation2DWaveSystemPStagTriangles(scaling):
  10     start = time.time()
  11     #### 2D Delaunay triangles meshes
  12     meshList=['squareWithTriangles_1','squareWithTriangles_2','squareWithTriangles_3','squareWithTriangles_4']
  13     meshType="Unstructured_Delaunay_triangles"
  14     testColor="Green"
  15     nbMeshes=len(meshList)
  16     mesh_size_tab=[0]*nbMeshes
  17     mesh_path='../../../ressources/2DTriangles/'
  18     mesh_name='squareWithTriangles'
  19     resolution=100
  20     curv_abs=np.linspace(0,sqrt(2),resolution+1)
  21
  22     error_p_tab=[0]*nbMeshes
  23     error_u_tab=[0]*nbMeshes
  24     diag_data_press=[0]*nbMeshes
  25     diag_data_vel=[0]*nbMeshes
  26     time_tab=[0]*nbMeshes
  27     t_final=[0]*nbMeshes
  28     ndt_final=[0]*nbMeshes
  29     max_vel=[0]*nbMeshes
  30     cond_number=[0]*nbMeshes
  31
  32     plt.close('all')
  33     i=0
  34     cfl=0.5
  35     # Storing of numerical errors, mesh sizes and diagonal values
  36     for filename in meshList:
  37         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], cond_number[i] = WaveSystemPStag.solve_file(mesh_path+filename, mesh_name, resolution,scaling,meshType,testColor,cfl,"Periodic")
  38         assert max_vel[i]>0.3 and max_vel[i]<1.1
  39         error_p_tab[i]=log10(error_p_tab[i])
  40         error_u_tab[i]=log10(error_u_tab[i])
  41         time_tab[i]=log10(time_tab[i])
  42         i=i+1
  43
  44     end = time.time()
  45
  46     # Plot over diagonal line
  47     for i in range(nbMeshes):
  48         if(scaling==0):
  49             plt.plot(curv_abs, diag_data_press[i], label= str(mesh_size_tab[i]) + ' cells - no scaling')
  50         else:
  51             plt.plot(curv_abs, diag_data_press[i], label= str(mesh_size_tab[i]) + ' cells - with scaling')
  52     plt.legend()
  53     plt.xlabel('Position on diagonal line')
  54     plt.ylabel('Pressure on diagonal line')
  55     plt.title('Plot over diagonal line for stationary wave system \n with PStagggered scheme on 2D Delaunay triangles meshes')
  56     plt.savefig(mesh_name+'_Pressure_2DWaveSystemPStag_'+"scaling"+str(scaling)+"_PlotOverDiagonalLine.png")
  57     plt.close()
  58
  59     plt.clf()
  60     for i in range(nbMeshes):
  61         if(scaling==0):
  62             plt.plot(curv_abs, diag_data_vel[i],   label= str(mesh_size_tab[i]) + ' cells - no scaling')
  63         else:
  64             plt.plot(curv_abs, diag_data_vel[i],   label= str(mesh_size_tab[i]) + ' cells - with scaling')
  65     plt.legend()
  66     plt.xlabel('Position on diagonal line')
  67     plt.ylabel('Velocity on diagonal line')
  68     plt.title('Plot over diagonal line for the stationary wave system \n with PStagggered scheme on 2D Delaunay triangles meshes')
  69     plt.savefig(mesh_name+"_Velocity_2DWaveSystemPStag_"+"scaling"+str(scaling)+"_PlotOverDiagonalLine.png")
  70     plt.close()
  71
  72     # Plot of number of time steps
  73     plt.close()
  74     if(scaling==0):
  75         plt.plot(mesh_size_tab, ndt_final, label='Number of time step to reach stationary regime - no scaling')
  76     else:
  77         plt.plot(mesh_size_tab, ndt_final, label='Number of time step to reach stationary regime - with scaling')
  78     plt.legend()
  79     plt.xlabel('Number of cells')
  80     plt.ylabel('Max time steps for stationary regime')
  81     plt.title('Number of times steps required for the stationary Wave System \n with PStagggered scheme on 2D Delaunay triangles meshes')
  82     plt.savefig(mesh_name+"_2DWaveSystemPStag_"+"scaling"+str(scaling)+"_TimeSteps.png")
  83
  84     # Plot of number of stationary time
  85     plt.close()
  86     if(scaling==0):
  87         plt.plot(mesh_size_tab, t_final, label='Time where stationary regime is reached - no scaling')
  88     else:
  89         plt.plot(mesh_size_tab, t_final, label='Time where stationary regime is reached - with scaling')
  90     plt.legend()
  91     plt.xlabel('number of cells')
  92     plt.ylabel('Max time for stationary regime')
  93     plt.title('Simulated time for the stationary Wave System \n with PStagggered scheme on 2D Delaunay triangles meshes')
  94     plt.savefig(mesh_name+"_2DWaveSystemPStag_"+"scaling"+str(scaling)+"_FinalTime.png")
  95
  96     # Plot of number of maximal velocity norm
  97     plt.close()
  98     if(scaling==0):
  99         plt.plot(mesh_size_tab, max_vel, label='Maximum velocity norm - no scaling')
 100     else:
 101         plt.plot(mesh_size_tab, max_vel, label='Maximum velocity norm - with scaling')
 102     plt.legend()
 103     plt.xlabel('Number of cells')
 104     plt.ylabel('Max velocity norm')
 105     plt.title('Maximum velocity norm for the stationary Wave System \n with PStagggered scheme on 2D Delaunay triangles meshes')
 106     plt.savefig(mesh_name+"_2DWaveSystemPStag_"+"scaling"+str(scaling)+"_MaxVelNorm.png")
 107
 108     # Plot of condition number
 109     plt.close()
 110     if(scaling==0):
 111         plt.plot(mesh_size_tab, cond_number, label='Condition number - no scaling')
 112     else:
 113         plt.plot(mesh_size_tab, cond_number, label='Condition number - with scaling')
 114     plt.legend()
 115     plt.xlabel('Number of cells')
 116     plt.ylabel('Condition number')
 117     plt.title('Condition number for the stationary Wave System \n with PStagggered scheme on 2D square meshes')
 118     plt.savefig(mesh_name+"_2DWaveSystemPStag_"+"scaling"+str(scaling)+"_condition_number.png")
 119
 120     for i in range(nbMeshes):
 121         mesh_size_tab[i]=0.5*log10(mesh_size_tab[i])
 122
 123     # Least square linear regression
 124     # Find the best a,b such that f(x)=ax+b best approximates the convergence curve
 125     # The vector X=(a,b) solves a symmetric linear system AX=B with A=(a1,a2\\a2,a3), B=(b1,b2)
 126     a1=np.dot(mesh_size_tab,mesh_size_tab)
 127     a2=np.sum(mesh_size_tab)
 128     a3=nbMeshes
 129
 130     det=a1*a3-a2*a2
 131     assert det!=0, 'test_validation2DWaveSystemFVPStag() : Make sure you use distinct meshes and at least two meshes'
 132
 133     b1p=np.dot(error_p_tab,mesh_size_tab)
 134     b2p=np.sum(error_p_tab)
 135     ap=( a3*b1p-a2*b2p)/det
 136     bp=(-a2*b1p+a1*b2p)/det
 137
 138     print("FV PStag on 2D Delaunay triangles meshes : scheme order for pressure is ", -ap)
 139
 140     b1u=np.dot(error_u_tab,mesh_size_tab)
 141     b2u=np.sum(error_u_tab)
 142     au=( a3*b1u-a2*b2u)/det
 143     bu=(-a2*b1u+a1*b2u)/det
 144
 145     if(scaling==0):
 146         print("FV PStag on 2D Delaunay triangles meshes : scheme order for velocity without scaling is ", -au)
 147     else:
 148         print("FV PStag on 2D Delaunay triangles meshes : scheme order for velocity with scaling is ", -au)
 149
 150     # Plot of convergence curves
 151     plt.close()
 152     if(scaling==0):
 153         plt.plot(mesh_size_tab, error_p_tab, label='|error on stationary pressure| - no scaling')
 154     else:
 155         plt.plot(mesh_size_tab, error_p_tab, label='|error on stationary pressure| - with scaling')
 156     plt.legend()
 157     plt.xlabel('1/2 log(number of cells)')
 158     plt.ylabel('|error p|')
 159     plt.title('Convergence of finite volumes for the stationary Wave System \n with PStagggered scheme on 2D Delaunay triangles meshes')
 160     plt.savefig(mesh_name+"_Pressure_2DWaveSystemPStag_"+"scaling"+str(scaling)+"_ConvergenceCurve.png")
 161
 162     plt.close()
 163     if(scaling==0):
 164         plt.plot(mesh_size_tab, error_u_tab, label='log(|error on stationary velocity|) - no scaling')
 165     else:
 166         plt.plot(mesh_size_tab, error_u_tab, label='log(|error on stationary velocity|) - with scaling')
 167     plt.legend()
 168     plt.xlabel('1/2 log(number of cells)')
 169     plt.ylabel('|error u|')
 170     plt.title('Convergence of finite volumes for the stationary Wave System \n with PStagggered scheme on 2D Delaunay triangles meshes')
 171     plt.savefig(mesh_name+"_Velocity_2DWaveSystemPStag_"+"scaling"+str(scaling)+"_ConvergenceCurve.png")
 172
 173     # Plot of computational time
 174     plt.close()
 175     if(scaling==0):
 176         plt.plot(mesh_size_tab, time_tab, label='log(cpu time) - no scaling')
 177     else:
 178         plt.plot(mesh_size_tab, time_tab, label='log(cpu time) - with scaling')
 179     plt.legend()
 180     plt.xlabel('1/2 log(number of cells)')
 181     plt.ylabel('log(cpu time)')
 182     plt.title('Computational time of finite volumes for the stationary Wave System \n with PStagggered scheme on 2D Delaunay triangles meshes')
 183     plt.savefig(mesh_name+"2DWaveSystemPStag_"+"scaling"+str(scaling)+"_ComputationalTime.png")
 184
 185     plt.close('all')
 186
 187     convergence_synthesis={}
 188
 189     convergence_synthesis["PDE_model"]="Wave system"
 190     convergence_synthesis["PDE_is_stationary"]=False
 191     convergence_synthesis["PDE_search_for_stationary_solution"]=True
 192     if(scaling==0):
 193         convergence_synthesis["Numerical_method_name"]="PStag no scaling"
 194     else:
 195         convergence_synthesis["Numerical_method_name"]="PStag scaling"
 196     convergence_synthesis["Numerical_method_space_discretization"]="Finite volumes"
 197     convergence_synthesis["Numerical_method_time_discretization"]="Implicit"
 198     convergence_synthesis["Initial_data"]="Constant pressure, divergence free velocity"
 199     convergence_synthesis["Boundary_conditions"]="Periodic"
 200     convergence_synthesis["Numerical_parameter_cfl"]=cfl
 201     convergence_synthesis["Space_dimension"]=2
 202     convergence_synthesis["Mesh_dimension"]=2
 203     convergence_synthesis["Mesh_names"]=meshList
 204     convergence_synthesis["Mesh_type"]=meshType
 205     convergence_synthesis["Mesh_path"]=mesh_path
 206     convergence_synthesis["Mesh_description"]=mesh_name
 207     convergence_synthesis["Mesh_sizes"]=mesh_size_tab
 208     convergence_synthesis["Mesh_cell_type"]="Triangles"
 209     convergence_synthesis["Numerical_error_velocity"]=error_u_tab
 210     convergence_synthesis["Numerical_error_pressure"]=error_p_tab
 211     convergence_synthesis["Max_vel_norm"]=max_vel
 212     convergence_synthesis["Final_time"]=t_final
 213     convergence_synthesis["Final_time_step"]=ndt_final
 214     convergence_synthesis["Scheme_order"]=-au
 215     convergence_synthesis["Scheme_order_vel"]=-au
 216     convergence_synthesis["Scheme_order_press"]=-ap
 217     convergence_synthesis["Scaling_preconditioner"]=scaling
 218     convergence_synthesis["Condition_numbers"]=cond_number
 219     convergence_synthesis["Test_color"]=testColor
 220     convergence_synthesis["Computational_time"]=end-start
 221
 222     with open('Convergence_WaveSystem_2DFV_PStag_'+mesh_name+'.json', 'w') as outfile:
 223         json.dump(convergence_synthesis, outfile)
 224
 225 if __name__ == """__main__""":
 226     if len(sys.argv) >1 :
 227         scaling = int(sys.argv[1])
 228         test_validation2DWaveSystemPStagTriangles(scaling)
 229     else :
 230         test_validation2DWaveSystemPStagTriangles(2)
 231