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

   1 import WaveSystemCentered
   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_validation2DWaveSystemCenteredFlatCrossTriangles(scaling):
  10     start = time.time()
  11     #### 2D flat cross triangle meshes
  12     meshList=['squareWithFlatCrossTriangles_00','squareWithFlatCrossTriangles_0','squareWithFlatCrossTriangles_1']#,'squareWithFlatCrossTriangles_2','squareWithFlatCrossTriangles_3'
  13     meshType="Unstructured_flat_cross_triangles"
  14     testColor="Green"
  15     nbMeshes=len(meshList)
  16     mesh_size_tab=[0]*nbMeshes
  17     mesh_path='../../../ressources/2DFlatCrossTriangles/'
  18     mesh_name='squareWithFlatCrossTriangles'
  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] = WaveSystemCentered.solve_file(mesh_path+filename, mesh_name, resolution,scaling,meshType,testColor,cfl,"Periodic")
  38         assert max_vel[i]>0.94 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 centered scheme on 2D flat cross triangular meshes')
  56     plt.savefig(mesh_name+'_Pressure_2DWaveSystemCentered_'+"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 centered scheme on 2D flat cross triangle meshes')
  69     plt.savefig(mesh_name+"_Velocity_2DWaveSystemCentered_"+"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 centered scheme on 2D flat cross triangle meshes')
  82     plt.savefig(mesh_name+"_2DWaveSystemCentered_"+"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 centered scheme on 2D flat cross triangle meshes')
  94     plt.savefig(mesh_name+"_2DWaveSystemCentered_"+"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 centered scheme on 2D flat cross triangle meshes')
 106     plt.savefig(mesh_name+"_2DWaveSystemCentered_"+"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 centered scheme on 2D square meshes')
 118     plt.savefig(mesh_name+"_2DWaveSystemCentered_"+"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_validation2DWaveSystemFVCentered() : 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     if(scaling==0):
 139         print("FV Centered on 2D flat cross triangle meshes : scheme order for pressure without scaling is ", -ap)
 140     else:
 141         print("FV Centered on 2D flat cross triangle meshes : scheme order for pressure with scaling is ", -ap)
 142
 143     b1u=np.dot(error_u_tab,mesh_size_tab)
 144     b2u=np.sum(error_u_tab)
 145     au=( a3*b1u-a2*b2u)/det
 146     bu=(-a2*b1u+a1*b2u)/det
 147
 148     if(scaling==0):
 149         print("FV Centered on 2D flat cross triangle meshes : scheme order for velocity without scaling is ", -au)
 150     else:
 151         print("FV Centered on 2D flat cross triangle meshes : scheme order for velocity with scaling is ", -au)
 152
 153     # Plot of convergence curves
 154     plt.close()
 155     if(scaling==0):
 156         plt.plot(mesh_size_tab, error_p_tab, label='|error on stationary pressure| - no scaling')
 157     else:
 158         plt.plot(mesh_size_tab, error_p_tab, label='|error on stationary pressure| - with scaling')
 159     plt.legend()
 160     plt.xlabel('1/2 log(Number of cells)')
 161     plt.ylabel('log(|error p|)')
 162     plt.title('Convergence of finite volumes for the stationary Wave System \n with centered scheme on 2D flat cross triangle meshes')
 163     plt.savefig(mesh_name+"_Pressure_2DWaveSystemCentered_"+"scaling"+str(scaling)+"_ConvergenceCurve.png")
 164
 165     plt.close()
 166     if(scaling==0):
 167         plt.plot(mesh_size_tab, error_u_tab, label='log(|error on stationary velocity|) - no scaling')
 168     else:
 169         plt.plot(mesh_size_tab, error_u_tab, label='log(|error on stationary velocity|) - with scaling')
 170     plt.legend()
 171     plt.xlabel('1/2 log(Number of cells)')
 172     plt.ylabel('log(|error u|)')
 173     plt.title('Convergence of finite volumes for the stationary Wave System \n with centered scheme on 2D flat cross triangle meshes')
 174     plt.savefig(mesh_name+"_Velocity_2DWaveSystemCentered_"+"scaling"+str(scaling)+"_ConvergenceCurve.png")
 175
 176     # Plot of computational time
 177     plt.close()
 178     if(scaling==0):
 179         plt.plot(mesh_size_tab, time_tab, label='log(cpu time) - no scaling')
 180     else:
 181         plt.plot(mesh_size_tab, time_tab, label='log(cpu time) - with scaling')
 182     plt.legend()
 183     plt.xlabel('1/2 log(Number of cells)')
 184     plt.ylabel('log(cpu time)')
 185     plt.title('Computational time of finite volumes for the stationary Wave System \n with centered scheme on 2D flat cross triangle meshes')
 186     plt.savefig(mesh_name+"2DWaveSystemCentered_"+"scaling"+str(scaling)+"_ComputationalTime.png")
 187
 188     plt.close('all')
 189
 190     convergence_synthesis={}
 191
 192     convergence_synthesis["PDE_model"]="Wave system"
 193     convergence_synthesis["PDE_is_stationary"]=False
 194     convergence_synthesis["PDE_search_for_stationary_solution"]=True
 195     if(scaling==0):
 196         convergence_synthesis["Numerical_method_name"]="Centered scheme no scaling"
 197     else:
 198         convergence_synthesis["Numerical_method_name"]="Centered scheme scaling"
 199     convergence_synthesis["Numerical_method_space_discretization"]="Finite volumes"
 200     convergence_synthesis["Numerical_method_time_discretization"]="Implicit"
 201     convergence_synthesis["Initial_data"]="Constant pressure, divergence free velocity"
 202     convergence_synthesis["Boundary_conditions"]="Periodic"
 203     convergence_synthesis["Numerical_parameter_cfl"]=cfl
 204     convergence_synthesis["Space_dimension"]=2
 205     convergence_synthesis["Mesh_dimension"]=2
 206     convergence_synthesis["Mesh_names"]=meshList
 207     convergence_synthesis["Mesh_type"]=meshType
 208     convergence_synthesis["Mesh_path"]=mesh_path
 209     convergence_synthesis["Mesh_description"]=mesh_name
 210     convergence_synthesis["Mesh_sizes"]=mesh_size_tab
 211     convergence_synthesis["Mesh_cell_type"]="Triangles"
 212     convergence_synthesis["Numerical_error_velocity"]=error_u_tab
 213     convergence_synthesis["Numerical_error_pressure"]=error_p_tab
 214     convergence_synthesis["Max_vel_norm"]=max_vel
 215     convergence_synthesis["Final_time"]=t_final
 216     convergence_synthesis["Final_time_step"]=ndt_final
 217     convergence_synthesis["Scheme_order"]=-au
 218     convergence_synthesis["Scheme_order_vel"]=-au
 219     convergence_synthesis["Scaling_preconditioner"]=scaling
 220     convergence_synthesis["Condition_numbers"]=cond_number
 221     convergence_synthesis["Test_color"]=testColor
 222     convergence_synthesis["Computational_time"]=end-start
 223
 224     with open('Convergence_WaveSystem_2DFV_Centered_'+mesh_name+'.json', 'w') as outfile:
 225         json.dump(convergence_synthesis, outfile)
 226
 227 if __name__ == """__main__""":
 228     if len(sys.argv) >1 :
 229         scaling = int(sys.argv[1])
 230         test_validation2DWaveSystemCenteredFlatCrossTriangles(scaling)
 231     else :
 232         test_validation2DWaveSystemCenteredFlatCrossTriangles(2)
 233