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_validation2DWaveSystemSourceCentered_brickwall(scaling):
  10     start = time.time()
  11     #### 2D triangular meshes
  12     meshList=['squareWithBrickWall_1','squareWithBrickWall_2','squareWithBrickWall_3','squareWithBrickWall_4']
  13     meshType="Unstructured brick wall"
  14     testColor="Green"
  15     nbMeshes=len(meshList)
  16     mesh_size_tab=[0]*nbMeshes
  17     mesh_path='../../../ressources/2DBrickWall/'
  18     mesh_name='squareWithBrickWall'
  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",True)
  38         assert max_vel[i]>0.001 and max_vel[i]<0.01
  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 with source term \n with centered scheme on 2D triangular meshes')
  56     plt.savefig(mesh_name+'_Pressure_2DWaveSystemSourceCentered_'+"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 with source term \n with centered scheme on 2D triangular meshes')
  69     plt.savefig(mesh_name+"_Velocity_2DWaveSystemSourceCentered_"+"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 with source term \n with centered scheme on 2D triangular meshes')
  82     plt.savefig(mesh_name+"_2DWaveSystemSourceCentered_"+"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 with source term \n with centered scheme on 2D triangular meshes')
  94     plt.savefig(mesh_name+"_2DWaveSystemSourceCentered_"+"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 with source term \n with centered scheme on 2D triangular meshes')
 106     plt.savefig(mesh_name+"_2DWaveSystemSourceCentered_"+"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 with source term \n with centered scheme on 2D square meshes')
 118     plt.savefig(mesh_name+"_2DWaveSystemSourceCentered_"+"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_validation2DWaveSystemBrickWallFVSourceCentered() : Make sure you use distinct meshes and at least two meshes'
 132
 133     b1u=np.dot(error_u_tab,mesh_size_tab)
 134     b2u=np.sum(error_u_tab)
 135     au=( a3*b1u-a2*b2u)/det
 136     bu=(-a2*b1u+a1*b2u)/det
 137
 138     if(scaling==0):
 139         print( "FVCentered on 2D triangular meshes : scheme order for velocity without scaling is ", -au)
 140     else:
 141         print( "FVCentered on 2D triangular meshes : scheme order for velocity with    scaling is ", -au)
 142
 143     # Plot of convergence curves
 144     plt.close()
 145     if(scaling==0):
 146         plt.plot(mesh_size_tab, error_p_tab, label='|error on stationary pressure| - no scaling')
 147     else:
 148         plt.plot(mesh_size_tab, error_p_tab, label='|error on stationary pressure| - with scaling')
 149     plt.legend()
 150     plt.xlabel('1/2 log(Number of cells)')
 151     plt.ylabel('log(|error p|)')
 152     plt.title('Convergence of finite volumes for the stationary Wave System with source term \n with centered scheme on 2D triangular meshes')
 153     plt.savefig(mesh_name+"_Pressure_2DWaveSystemSourceCentered_"+"scaling"+str(scaling)+"_ConvergenceCurve.png")
 154
 155     plt.close()
 156     if(scaling==0):
 157         plt.plot(mesh_size_tab, error_u_tab, label='log(|error on stationary velocity|) - no scaling')
 158     else:
 159         plt.plot(mesh_size_tab, error_u_tab, label='log(|error on stationary velocity|) - with scaling')
 160     plt.legend()
 161     plt.xlabel('1/2 log(Number of cells)')
 162     plt.ylabel('log(|error u|)')
 163     plt.title('Convergence of finite volumes for the stationary Wave System with source term \n with centered scheme on 2D triangular meshes')
 164     plt.savefig(mesh_name+"_Velocity_2DWaveSystemSourceCentered_"+"scaling"+str(scaling)+"_ConvergenceCurve.png")
 165
 166     # Plot of computational time
 167     plt.close()
 168     if(scaling==0):
 169         plt.plot(mesh_size_tab, time_tab, label='log(cpu time) - no scaling')
 170     else:
 171         plt.plot(mesh_size_tab, time_tab, label='log(cpu time) - with scaling')
 172     plt.legend()
 173     plt.xlabel('1/2 log(Number of cells)')
 174     plt.ylabel('log(cpu time)')
 175     plt.title('Computational time of finite volumes for the stationary Wave System with source term \n with centered scheme on 2D triangular meshes')
 176     plt.savefig(mesh_name+"2DWaveSystemSourceCentered_"+"scaling"+str(scaling)+"_ComputationalTime.png")
 177
 178     plt.close('all')
 179
 180     convergence_synthesis={}
 181
 182     convergence_synthesis["PDE_model"]="Wave system"
 183     convergence_synthesis["PDE_is_stationary"]=False
 184     convergence_synthesis["PDE_search_for_stationary_solution"]=True
 185     if(scaling==0):
 186         convergence_synthesis["Numerical_method_name"]="Centered scheme no scaling"
 187     else:
 188         convergence_synthesis["Numerical_method_name"]="Centered scheme scaling"
 189     convergence_synthesis["Numerical_method_space_discretization"]="Finite volumes"
 190     convergence_synthesis["Numerical_method_time_discretization"]="Implicit"
 191     convergence_synthesis["Initial_data"]="Constant pressure, divergence free velocity"
 192     convergence_synthesis["Boundary_conditions"]="Periodic"
 193     convergence_synthesis["Numerical_parameter_cfl"]=cfl
 194     convergence_synthesis["Space_dimension"]=2
 195     convergence_synthesis["Mesh_dimension"]=2
 196     convergence_synthesis["Mesh_names"]=meshList
 197     convergence_synthesis["Mesh_type"]=meshType
 198     convergence_synthesis["Mesh_path"]=mesh_path
 199     convergence_synthesis["Mesh_description"]=mesh_name
 200     convergence_synthesis["Mesh_sizes"]=mesh_size_tab
 201     convergence_synthesis["Mesh_cell_type"]="BrickWall"
 202     convergence_synthesis["Numerical_error_velocity"]=error_u_tab
 203     convergence_synthesis["Numerical_error_pressure"]=error_p_tab
 204     convergence_synthesis["Max_vel_norm"]=max_vel
 205     convergence_synthesis["Final_time"]=t_final
 206     convergence_synthesis["Final_time_step"]=ndt_final
 207     convergence_synthesis["Scheme_order"]=-au
 208     convergence_synthesis["Scheme_order_vel"]=-au
 209     convergence_synthesis["Scaling_preconditioner"]=scaling
 210     convergence_synthesis["Condition_numbers"]=cond_number
 211     convergence_synthesis["Test_color"]=testColor
 212     convergence_synthesis["Computational_time"]=end-start
 213
 214     with open('Convergence_WaveSystemSource_2DFV_Centered_'+mesh_name+'.json', 'w') as outfile:
 215         json.dump(convergence_synthesis, outfile)
 216
 217 if __name__ == """__main__""":
 218     if len(sys.argv) >1 :
 219         scaling = int(sys.argv[1])
 220         test_validation2DWaveSystemSourceCentered_brickwall(scaling)
 221     else :
 222         test_validation2DWaveSystemSourceCentered_brickwall(2)
 223