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

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