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