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

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