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