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 def test_validation2DWaveSystemPStagHexagons(bctype,scaling):
  11     start = time.time()
  12     #### 2D hexagonal mesh
  13     meshList=['squareWithHexagons_1','squareWithHexagons_2','squareWithHexagons_3','squareWithHexagons_4']#,'squareWithHexagons_5'
  14     meshType="Regular hexagons"
  15     testColor="Orange : Pb with wall BC (no mass conservation)"
  16     nbMeshes=len(meshList)
  17     error_p_tab=[0]*nbMeshes
  18     error_u_tab=[0]*nbMeshes
  19     mesh_size_tab=[0]*nbMeshes
  20     mesh_path='../../../ressources/2DHexagons/'
  21     mesh_name='squareWithHexagons'
  22     diag_data_press=[0]*nbMeshes
  23     diag_data_vel=[0]*nbMeshes
  24     time_tab=[0]*nbMeshes
  25     t_final=[0]*nbMeshes
  26     ndt_final=[0]*nbMeshes
  27     max_vel=[0]*nbMeshes
  28     cond_number=[0]*nbMeshes
  29     resolution=100
  30     curv_abs=np.linspace(0,sqrt(2),resolution+1)
  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] =WaveSystemPStag.solve_file(mesh_path+filename, mesh_name, resolution,scaling, meshType,testColor,cfl,bctype)
  38         assert max_vel[i]>0.0001 and max_vel[i]<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         plt.plot(curv_abs, diag_data_press[i], label= str(mesh_size_tab[i]) + ' cells')
  49     plt.legend()
  50     plt.xlabel('Position on diagonal line')
  51     plt.ylabel('Pressure on diagonal line')
  52     plt.title('Plot over diagonal line for stationary wave system \n with PStag scheme on 2D hexagonal meshes')
  53     plt.savefig(mesh_name+'_Pressure_2DWaveSystemPStag_'+"PlotOverDiagonalLine.png")
  54     plt.close()
  55
  56     plt.clf()
  57     for i in range(nbMeshes):
  58         plt.plot(curv_abs, diag_data_vel[i],   label= str(mesh_size_tab[i]) + ' cells')
  59     plt.legend()
  60     plt.xlabel('Position on diagonal line')
  61     plt.ylabel('Velocity on diagonal line')
  62     plt.title('Plot over diagonal line for the stationary wave system \n with PStag scheme on 2D hexagonal meshes')
  63     plt.savefig(mesh_name+"_Velocity_2DWaveSystemPStag_"+"PlotOverDiagonalLine.png")
  64     plt.close()
  65
  66     # Plot of number of time steps
  67     plt.close()
  68     plt.plot(mesh_size_tab, ndt_final, label='Number of time step to reach stationary regime')
  69     plt.legend()
  70     plt.xlabel('Number of cells')
  71     plt.ylabel('Max time steps for stationary regime')
  72     plt.title('Number of times steps required for the stationary Wave System \n with PStag scheme on 2D hexagonal meshes')
  73     plt.savefig(mesh_name+"_2DWaveSystemPStag_"+"TimeSteps.png")
  74
  75     # Plot of number of stationary time
  76     plt.close()
  77     plt.plot(mesh_size_tab, t_final, label='Time where stationary regime is reached')
  78     plt.legend()
  79     plt.xlabel('Number of cells')
  80     plt.ylabel('Max time for stationary regime')
  81     plt.title('Simulated time for the stationary Wave System \n with PStag scheme on 2D hexagonal meshes')
  82     plt.savefig(mesh_name+"_2DWaveSystemPStag_"+"TimeFinal.png")
  83
  84     # Plot of number of maximal velocity norm
  85     plt.close()
  86     plt.plot(mesh_size_tab, max_vel, label='Maximum velocity norm')
  87     plt.legend()
  88     plt.xlabel('Number of cells')
  89     plt.ylabel('Max velocity norm')
  90     plt.title('Maximum velocity norm for the stationary Wave System \n with PStag scheme on 2D hexagonal meshes')
  91     plt.savefig(mesh_name+"_2DWaveSystemPStag_"+"MaxVelNorm.png")
  92
  93     for i in range(nbMeshes):
  94         mesh_size_tab[i]=0.5*log10(mesh_size_tab[i])
  95
  96     # Least square linear regression
  97     # Find the best a,b such that f(x)=ax+b best approximates the convergence curve
  98     # The vector X=(a,b) solves a symmetric linear system AX=B with A=(a1,a2\\a2,a3), B=(b1,b2)
  99     a1=np.dot(mesh_size_tab,mesh_size_tab)
 100     a2=np.sum(mesh_size_tab)
 101     a3=nbMeshes
 102
 103     det=a1*a3-a2*a2
 104     assert det!=0, 'test_validation2DWaveSystemPStagHexagonsFV() : Make sure you use distinct meshes and at least two meshes'
 105
 106     b1p=np.dot(error_p_tab,mesh_size_tab)
 107     b2p=np.sum(error_p_tab)
 108     ap=( a3*b1p-a2*b2p)/det
 109     bp=(-a2*b1p+a1*b2p)/det
 110
 111     print("FV PStag on 2D hexagonal meshes : scheme order for pressure is ", -ap)
 112
 113     b1u=np.dot(error_u_tab,mesh_size_tab)
 114     b2u=np.sum(error_u_tab)
 115     au=( a3*b1u-a2*b2u)/det
 116     bu=(-a2*b1u+a1*b2u)/det
 117
 118     print("FV PStag on 2D hexagonal meshes : scheme order for velocity is ", -au)
 119
 120     # Plot of convergence curves
 121     plt.close()
 122     plt.plot(mesh_size_tab, error_p_tab, label='|error on stationary pressure|')
 123     plt.legend()
 124     plt.xlabel('1/2 log(number of cells)')
 125     plt.ylabel('log(error p)')
 126     plt.title('Convergence of finite volumes for the stationary Wave System \n with PStag scheme on 2D hexagonal meshes')
 127     plt.savefig(mesh_name+"_Pressure_2DWaveSystemPStag_"+"ConvergenceCurve.png")
 128
 129     plt.close()
 130     plt.plot(mesh_size_tab, error_u_tab, label='|error on stationary velocity|')
 131     plt.legend()
 132     plt.xlabel('1/2 log(number of cells)')
 133     plt.ylabel('log(error u)')
 134     plt.title('Convergence of finite volumes for the stationary Wave System \n with PStag scheme on 2D hexagonal meshes')
 135     plt.savefig(mesh_name+"_Velocity_2DWaveSystemPStag_"+"ConvergenceCurve.png")
 136
 137     # Plot of computational time
 138     plt.close()
 139     plt.plot(mesh_size_tab, time_tab, label='log(cpu time)')
 140     plt.legend()
 141     plt.xlabel('1/2 log(number of cells)')
 142     plt.ylabel('log(cpu time)')
 143     plt.title('Computational time of finite volumes for the stationary Wave System \n with PStag scheme on 2D hexagonal meshes')
 144     plt.savefig(mesh_name+"_2DWaveSystemPStag_ComputationalTime.png")
 145
 146     plt.close('all')
 147
 148     convergence_synthesis={}
 149
 150     convergence_synthesis["PDE_model"]="Wave system"
 151     convergence_synthesis["PDE_is_stationary"]=False
 152     convergence_synthesis["PDE_search_for_stationary_solution"]=True
 153     convergence_synthesis["Numerical_method_name"]="PStag"
 154     convergence_synthesis["Numerical_method_space_discretization"]="Finite volumes"
 155     convergence_synthesis["Numerical_method_time_discretization"]="Implicit"
 156     convergence_synthesis["Initial_data"]="Constant pressure, divergence free velocity"
 157     convergence_synthesis["Boundary_conditions"]="Periodic"
 158     convergence_synthesis["Numerical_parameter_cfl"]=cfl
 159     convergence_synthesis["Space_dimension"]=2
 160     convergence_synthesis["Mesh_dimension"]=2
 161     convergence_synthesis["Mesh_names"]=meshList
 162     convergence_synthesis["Mesh_type"]=meshType
 163     convergence_synthesis["Mesh_path"]=mesh_path
 164     convergence_synthesis["Mesh_description"]=mesh_name
 165     convergence_synthesis["Mesh_sizes"]=mesh_size_tab
 166     convergence_synthesis["Mesh_cell_type"]="Hexagons"
 167     convergence_synthesis["Numerical_error_velocity"]=error_u_tab
 168     convergence_synthesis["Numerical_error_pressure"]=error_p_tab
 169     convergence_synthesis["Max_vel_norm"]=max_vel
 170     convergence_synthesis["Final_time"]=t_final
 171     convergence_synthesis["Final_time_step"]=ndt_final
 172     convergence_synthesis["Scheme_order"]=min(-au,-ap)
 173     convergence_synthesis["Scheme_order_vel"]=-au
 174     convergence_synthesis["Scheme_order_press"]=-ap
 175     convergence_synthesis["Scaling_preconditioner"]="None"
 176     convergence_synthesis["Test_color"]=testColor
 177     convergence_synthesis["Computational_time"]=end-start
 178
 179     with open('Convergence_WaveSystem_2DFV_PStag_'+mesh_name+'.json', 'w') as outfile:
 180         json.dump(convergence_synthesis, outfile)
 181
 182 if __name__ == """__main__""":
 183     if len(sys.argv) >2 :
 184         bctype = sys.argv[1]
 185         scaling = int(sys.argv[2])
 186         test_validation2DWaveSystemPStagHexagons(bctype,scaling)
 187     else :
 188         raise ValueError("test_validation2DWaveSystemPStagHexagons.py expects a mesh file name and a scaling parameter")