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

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