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