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

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