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

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