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

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