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

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