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

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