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

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