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

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