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

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