2 import matplotlib.pyplot as plt
4 from math import log10, sqrt
9 def test_validation2DWaveSystemSourcePStagFlatCrossTriangles(scaling):
11 #### 2D flat cross triangle meshes
12 meshList=['squareWithFlatCrossTriangles_00','squareWithFlatCrossTriangles_0']#,'squareWithFlatCrossTriangles_1','squareWithFlatCrossTriangles_2','squareWithFlatCrossTriangles_3'
13 meshType="Unstructured_flat_cross_triangles"
14 testColor="Orange (no stationary found on large meshes)"
15 nbMeshes=len(meshList)
16 mesh_size_tab=[0]*nbMeshes
17 mesh_path='../../../ressources/2DFlatCrossTriangles/'
18 mesh_name='squareWithFlatCrossTriangles'
20 curv_abs=np.linspace(0,sqrt(2),resolution+1)
22 error_p_tab=[0]*nbMeshes
23 error_u_tab=[0]*nbMeshes
24 diag_data_press=[0]*nbMeshes
25 diag_data_vel=[0]*nbMeshes
28 ndt_final=[0]*nbMeshes
30 cond_number=[0]*nbMeshes
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] = WaveSystemPStag.solve_file(mesh_path+filename, mesh_name, resolution,scaling,meshType,testColor,cfl,"Periodic",True)
38 assert max_vel[i]>0.001 and max_vel[i]<0.01
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])
46 # Plot over diagonal line
47 for i in range(nbMeshes):
49 plt.plot(curv_abs, diag_data_press[i], label= str(mesh_size_tab[i]) + ' cells - no scaling')
51 plt.plot(curv_abs, diag_data_press[i], label= str(mesh_size_tab[i]) + ' cells - with scaling')
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 with source term \n with PStagggered scheme on 2D flat cross triangle meshes')
56 plt.savefig(mesh_name+'_Pressure_2DWaveSystemSourcePStag_'+"scaling"+str(scaling)+"_PlotOverDiagonalLine.png")
60 for i in range(nbMeshes):
62 plt.plot(curv_abs, diag_data_vel[i], label= str(mesh_size_tab[i]) + ' cells - no scaling')
64 plt.plot(curv_abs, diag_data_vel[i], label= str(mesh_size_tab[i]) + ' cells - with scaling')
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 with source term \n with PStagggered scheme on 2D flat cross triangle meshes')
69 plt.savefig(mesh_name+"_Velocity_2DWaveSystemSourcePStag_"+"scaling"+str(scaling)+"_PlotOverDiagonalLine.png")
72 # Plot of number of time steps
75 plt.plot(mesh_size_tab, ndt_final, label='Number of time step to reach stationary regime - no scaling')
77 plt.plot(mesh_size_tab, ndt_final, label='Number of time step to reach stationary regime - with scaling')
79 plt.xlabel('Number of cells')
80 plt.ylabel('Max time steps for stationary regime')
81 plt.title('Number of times steps required for the stationary Wave System \n with PStagggered scheme on 2D flat cross triangle meshes')
82 plt.savefig(mesh_name+"_2DWaveSystemSquareSourcePStag_"+"scaling"+str(scaling)+"_TimeSteps.png")
84 # Plot of number of stationary time
87 plt.plot(mesh_size_tab, t_final, label='Time where stationary regime is reached - no scaling')
89 plt.plot(mesh_size_tab, t_final, label='Time where stationary regime is reached - with scaling')
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 PStagggered scheme on 2D flat cross triangle meshes')
94 plt.savefig(mesh_name+"_2DWaveSystemSourcePStag_"+"scaling"+str(scaling)+"_FinalTime.png")
96 # Plot of number of maximal velocity norm
99 plt.plot(mesh_size_tab, max_vel, label='Maximum velocity norm - no scaling')
101 plt.plot(mesh_size_tab, max_vel, label='Maximum velocity norm - with scaling')
103 plt.xlabel('Number of cells')
104 plt.ylabel('Max velocity norm')
105 plt.title('Maximum velocity norm for the stationary Wave System \n with PStagggered scheme on 2D flat cross triangle meshes')
106 plt.savefig(mesh_name+"_2DWaveSystemSourcePStag_"+"scaling"+str(scaling)+"_MaxVelNorm.png")
108 # Plot of condition number
111 plt.plot(mesh_size_tab, cond_number, label='Condition number - no scaling')
113 plt.plot(mesh_size_tab, cond_number, label='Condition number - with scaling')
115 plt.xlabel('Number of cells')
116 plt.ylabel('Condition number')
117 plt.title('Condition number for the stationary Wave System \n with PStagggered scheme on 2D square meshes')
118 plt.savefig(mesh_name+"_2DWaveSystemSourcePStag_"+"scaling"+str(scaling)+"_condition_number.png")
120 for i in range(nbMeshes):
121 mesh_size_tab[i]=0.5*log10(mesh_size_tab[i])
123 # Least square linear regression
124 # Find the best a,b such that f(x)=ax+b best approximates the convergence curve
125 # The vector X=(a,b) solves a symmetric linear system AX=B with A=(a1,a2\\a2,a3), B=(b1,b2)
126 a1=np.dot(mesh_size_tab,mesh_size_tab)
127 a2=np.sum(mesh_size_tab)
131 assert det!=0, 'test_validation2DWaveSystemSourceFVPStag() : Make sure you use distinct meshes and at least two meshes'
133 b1u=np.dot(error_u_tab,mesh_size_tab)
134 b2u=np.sum(error_u_tab)
135 au=( a3*b1u-a2*b2u)/det
136 bu=(-a2*b1u+a1*b2u)/det
139 print( "FVPStag on 2D flat cross triangle meshes : scheme order for velocity without scaling is ", -au)
141 print( "FVPStag on 2D flat cross triangle meshes : scheme order for velocity with scaling is ", -au)
143 # Plot of convergence curves
146 plt.plot(mesh_size_tab, error_p_tab, label='|error on stationary pressure| - no scaling')
148 plt.plot(mesh_size_tab, error_p_tab, label='|error on stationary pressure| - with scaling')
150 plt.xlabel('1/2 log(number of cells)')
151 plt.ylabel('|error p|')
152 plt.title('Convergence of finite volumes for the stationary Wave System \n with PStagggered scheme on 2D flat cross triangle meshes')
153 plt.savefig(mesh_name+"_Pressure_2DWaveSystemSourcePStag_"+"scaling"+str(scaling)+"_ConvergenceCurve.png")
157 plt.plot(mesh_size_tab, error_u_tab, label='log(|error on stationary velocity|) - no scaling')
159 plt.plot(mesh_size_tab, error_u_tab, label='log(|error on stationary velocity|) - with scaling')
161 plt.xlabel('1/2 log(number of cells)')
162 plt.ylabel('|error u|')
163 plt.title('Convergence of finite volumes for the stationary Wave System \n with PStagggered scheme on 2D flat cross triangle meshes')
164 plt.savefig(mesh_name+"_Velocity_2DWaveSystemSourcePStag_"+"scaling"+str(scaling)+"_ConvergenceCurve.png")
166 # Plot of computational time
169 plt.plot(mesh_size_tab, time_tab, label='log(cpu time) - no scaling')
171 plt.plot(mesh_size_tab, time_tab, label='log(cpu time) - with scaling')
173 plt.xlabel('1/2 log(number of cells)')
174 plt.ylabel('log(cpu time)')
175 plt.title('Computational time of finite volumes for the stationary Wave System \n with PStagggered scheme on 2D flat cross triangle meshes')
176 plt.savefig(mesh_name+"2DWaveSystemSourcePStag_"+"scaling"+str(scaling)+"_ComputationalTime.png")
180 convergence_synthesis={}
182 convergence_synthesis["PDE_model"]="Wave system"
183 convergence_synthesis["PDE_is_stationary"]=False
184 convergence_synthesis["PDE_search_for_stationary_solution"]=True
186 convergence_synthesis["Numerical_method_name"]="PStag no scaling"
188 convergence_synthesis["Numerical_method_name"]="PStag scaling"
189 convergence_synthesis["Numerical_method_space_discretization"]="Finite volumes"
190 convergence_synthesis["Numerical_method_time_discretization"]="Implicit"
191 convergence_synthesis["Initial_data"]="Constant pressure, divergence free velocity"
192 convergence_synthesis["Boundary_conditions"]="Periodic"
193 convergence_synthesis["Numerical_parameter_cfl"]=cfl
194 convergence_synthesis["Space_dimension"]=2
195 convergence_synthesis["Mesh_dimension"]=2
196 convergence_synthesis["Mesh_names"]=meshList
197 convergence_synthesis["Mesh_type"]=meshType
198 convergence_synthesis["Mesh_path"]=mesh_path
199 convergence_synthesis["Mesh_description"]=mesh_name
200 convergence_synthesis["Mesh_sizes"]=mesh_size_tab
201 convergence_synthesis["Mesh_cell_type"]="Triangles"
202 convergence_synthesis["Numerical_error_velocity"]=error_u_tab
203 convergence_synthesis["Numerical_error_pressure"]=error_p_tab
204 convergence_synthesis["Max_vel_norm"]=max_vel
205 convergence_synthesis["Final_time"]=t_final
206 convergence_synthesis["Final_time_step"]=ndt_final
207 convergence_synthesis["Scheme_order"]=-au
208 convergence_synthesis["Scheme_order_vel"]=-au
209 convergence_synthesis["Scaling_preconditioner"]=scaling
210 convergence_synthesis["Condition_numbers"]=cond_number
211 convergence_synthesis["Test_color"]=testColor
212 convergence_synthesis["Computational_time"]=end-start
214 with open('Convergence_WaveSystemSource_2DFV_PStag_'+mesh_name+'.json', 'w') as outfile:
215 json.dump(convergence_synthesis, outfile)
217 if __name__ == """__main__""":
218 if len(sys.argv) >1 :
219 scaling = int(sys.argv[1])
220 test_validation2DWaveSystemSourcePStagFlatCrossTriangles(scaling)
222 test_validation2DWaveSystemSourcePStagFlatCrossTriangles(2)