From d34d9fa39f5612b596d3dd5abcfa7d014ad81d88 Mon Sep 17 00:00:00 2001 From: michael Date: Sun, 21 Mar 2021 15:14:20 +0100 Subject: [PATCH] Plot pressure convergence curve only if no zero singularity --- ...est_validation2DWaveSystemUpwindSquares.py | 33 ++++++++++--------- 1 file changed, 18 insertions(+), 15 deletions(-) diff --git a/CDMATH/tests/validation/WaveSystem_stationary/2DWaveSystemUpwindSquares/test_validation2DWaveSystemUpwindSquares.py b/CDMATH/tests/validation/WaveSystem_stationary/2DWaveSystemUpwindSquares/test_validation2DWaveSystemUpwindSquares.py index 4028b07..ec7ca5b 100755 --- a/CDMATH/tests/validation/WaveSystem_stationary/2DWaveSystemUpwindSquares/test_validation2DWaveSystemUpwindSquares.py +++ b/CDMATH/tests/validation/WaveSystem_stationary/2DWaveSystemUpwindSquares/test_validation2DWaveSystemUpwindSquares.py @@ -41,8 +41,11 @@ def test_validation2DWaveSystemUpwindSquares(bctype,scaling): 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,"square"+str(nx)+'x'+str(nx), resolution,scaling,meshType,testColor,cfl,bctype) #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) assert max_vel[i]>0.0001 and max_vel[i]<1. - print( "error_p_tab[i]=", error_p_tab[i], " zero leads to singularity in plotting of convergence curve") - #error_p_tab[i]=log10(error_p_tab[i]) + if(error_p_tab[i]>0): + error_p_tab[i]=log10(error_p_tab[i]) + else: + print( "error_p_tab[i]=", error_p_tab[i], " zero leads to singularity in plotting of convergence curve") + error_p_tab[i]=-16 error_u_tab[i]=log10(error_u_tab[i]) time_tab[i]=log10(time_tab[i]) i=i+1 @@ -110,10 +113,10 @@ def test_validation2DWaveSystemUpwindSquares(bctype,scaling): assert det!=0, 'test_validation2DWaveSystemUpwind_squares() : Make sure you use distinct meshes and at least two meshes' #zero pressure leads to singularity - # b1p=np.dot(error_p_tab,mesh_size_tab) - # b2p=np.sum(error_p_tab) - # ap=( a3*b1p-a2*b2p)/det - # bp=(-a2*b1p+a1*b2p)/det + b1p=np.dot(error_p_tab,mesh_size_tab) + b2p=np.sum(error_p_tab) + ap=( a3*b1p-a2*b2p)/det + bp=(-a2*b1p+a1*b2p)/det # print("FV upwind on 2D square meshes : scheme order for pressure is ", -ap) @@ -125,13 +128,13 @@ def test_validation2DWaveSystemUpwindSquares(bctype,scaling): print("FV upwind on 2D square meshes : scheme order for velocity is ", -au) # Plot of convergence curves - # plt.close() - # plt.plot(mesh_size_tab, error_p_tab, label='|error on stationary pressure|') - # plt.legend() - # plt.xlabel('1/2 log(number of cells)') - # plt.ylabel('|error p|') - # plt.title('Convergence of finite volumes for the stationary Wave System \n with upwind scheme on 2D square meshes') - # plt.savefig(mesh_name+"_Pressure_2DWaveSystemUpwind_"+"ConvergenceCurve.png") + plt.close() + plt.plot(mesh_size_tab, error_p_tab, label='|error on stationary pressure|') + plt.legend() + plt.xlabel('1/2 log(number of cells)') + plt.ylabel('|error p|') + plt.title('Convergence of finite volumes for the stationary Wave System \n with upwind scheme on 2D square meshes') + plt.savefig(mesh_name+"_Pressure_2DWaveSystemUpwind_"+"ConvergenceCurve.png") plt.close() plt.plot(mesh_size_tab, error_u_tab, label='log(|error on stationary velocity|)') @@ -176,9 +179,9 @@ def test_validation2DWaveSystemUpwindSquares(bctype,scaling): convergence_synthesis["Max_vel_norm"]=max_vel convergence_synthesis["Final_time"]=t_final convergence_synthesis["Final_time_step"]=ndt_final - convergence_synthesis["Scheme_order"]=-au#min(-au,-ap) + convergence_synthesis["Scheme_order"]=min(-au,-ap) convergence_synthesis["Scheme_order_vel"]=-au - #convergence_synthesis["Scheme_order_press"]=-ap + convergence_synthesis["Scheme_order_press"]=-ap convergence_synthesis["Scaling_preconditioner"]="None" convergence_synthesis["Test_color"]=testColor convergence_synthesis["Computational_time"]=end-start -- 2.39.2