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

   1 import WaveSystemCentered
   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_validation2DWaveSystemSourceCenteredCheckerboard(scaling):
   9     start = time.time()
  10     #### 2D checkerboard mesh
  11     meshList=['checkerboard_5x5','checkerboard_9x9','checkerboard_17x17','checkerboard_33x33','checkerboard_65x65']
  12     meshType="Regular checkerboard"
  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/2DCheckerboard/'
  19     mesh_name='squareWithCheckerboard'
  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     cond_number=[0]*nbMeshes
  29
  30     plt.close('all')
  31     i=0
  32     cfl=0.5
  33     # Storing of numerical errors, mesh sizes and diagonal values
  34     for filename in meshList:
  35         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] =WaveSystemCentered.solve_file(mesh_path+filename, mesh_name, resolution,scaling,meshType,testColor,cfl,"Periodic",True)
  36         assert max_vel[i]>0.001 and max_vel[i]<0.01
  37         error_p_tab[i]=log10(error_p_tab[i])
  38         error_u_tab[i]=log10(error_u_tab[i])
  39         time_tab[i]=log10(time_tab[i])
  40         i=i+1
  41
  42     end = time.time()
  43
  44     # Plot over diagonal line
  45     for i in range(nbMeshes):
  46         if(scaling==0):
  47             plt.plot(curv_abs, diag_data_press[i], label= str(mesh_size_tab[i]) + ' cells - no scaling')
  48         else:
  49             plt.plot(curv_abs, diag_data_press[i], label= str(mesh_size_tab[i]) + ' cells - with scaling')
  50     plt.legend()
  51     plt.xlabel('Position on diagonal line')
  52     plt.ylabel('Pressure on diagonal line')
  53     plt.title('Plot over diagonal line for stationary Wave System with source term \n with centered scheme on 2D checkerboard meshes')
  54     plt.savefig(mesh_name+'_Pressure_2DWaveSystemSourceCentered_'+"scaling"+str(scaling)+"_PlotOverDiagonalLine.png")
  55     plt.close()
  56
  57     # Plot of number of time steps
  58     plt.close()
  59     if(scaling==0):
  60         plt.plot(mesh_size_tab, ndt_final, label='Number of time step to reach stationary regime - no scaling')
  61     else:
  62         plt.plot(mesh_size_tab, ndt_final, label='Number of time step to reach stationary regime - with scaling')
  63     plt.legend()
  64     plt.xlabel('Number of cells')
  65     plt.ylabel('Max time steps for stationary regime')
  66     plt.title('Number of times steps required for the stationary Wave System \n with centered scheme on 2D checkerboard meshes')
  67     plt.savefig(mesh_name+"_2DWaveSystemSourceCentered_"+"scaling"+str(scaling)+"_TimeSteps.png")
  68
  69     # Plot of number of stationary time
  70     plt.close()
  71     if(scaling==0):
  72         plt.plot(mesh_size_tab, t_final, label='Time where stationary regime is reached - no scaling')
  73     else:
  74         plt.plot(mesh_size_tab, t_final, label='Time where stationary regime is reached - with scaling')
  75     plt.legend()
  76     plt.xlabel('Number of cells')
  77     plt.ylabel('Max time for stationary regime')
  78     plt.title('Simulated time for the stationary Wave System \n with centered scheme on 2D checkerboard meshes')
  79     plt.savefig(mesh_name+"_2DWaveSystemSourceCentered_"+"scaling"+str(scaling)+"_TimeFinal.png")
  80
  81     # Plot of number of maximal velocity norm
  82     plt.close()
  83     if(scaling==0):
  84         plt.plot(mesh_size_tab, max_vel, label='Maximum velocity norm - no scaling')
  85     else:
  86         plt.plot(mesh_size_tab, max_vel, label='Maximum velocity norm - with scaling')
  87     plt.legend()
  88     plt.xlabel('Number of cells')
  89     plt.ylabel('Max velocity norm')
  90     plt.title('Maximum velocity norm for the stationary Wave System \n with centered scheme on 2D checkerboard meshes')
  91     plt.savefig(mesh_name+"_2DWaveSystemSourceCentered_"+"scaling"+str(scaling)+"_MaxVelNorm.png")
  92
  93     # Plot of condition number
  94     plt.close()
  95     if(scaling==0):
  96         plt.plot(mesh_size_tab, cond_number, label='Condition number - no scaling')
  97     else:
  98         plt.plot(mesh_size_tab, cond_number, label='Condition number - with scaling')
  99     plt.legend()
 100     plt.xlabel('Number of cells')
 101     plt.ylabel('Condition number')
 102     plt.title('Condition number for the stationary Wave System \n with centered scheme on 2D square meshes')
 103     plt.savefig(mesh_name+"_2DWaveSystemSourceCentered_"+"scaling"+str(scaling)+"_condition_number.png")
 104
 105     plt.clf()
 106     for i in range(nbMeshes):
 107         if(scaling==0):
 108             plt.plot(curv_abs, diag_data_vel[i],   label= str(mesh_size_tab[i]) + ' cells - no scaling')
 109         else:
 110             plt.plot(curv_abs, diag_data_vel[i],   label= str(mesh_size_tab[i]) + ' cells - with scaling')
 111     plt.legend()
 112     plt.xlabel('Position on diagonal line')
 113     plt.ylabel('Velocity on diagonal line')
 114     plt.title('Plot over diagonal line for the stationary Wave System with source term \n with centered scheme on 2D checkerboard meshes')
 115     plt.savefig(mesh_name+"_Velocity_2DWaveSystemSourceCentered_"+"scaling"+str(scaling)+"_PlotOverDiagonalLine.png")
 116     plt.close()
 117
 118     for i in range(nbMeshes):
 119         mesh_size_tab[i]=0.5*log10(mesh_size_tab[i])
 120
 121     # Least square linear regression
 122     # Find the best a,b such that f(x)=ax+b best approximates the convergence curve
 123     # The vector X=(a,b) solves a symmetric linear system AX=B with A=(a1,a2\\a2,a3), B=(b1,b2)
 124     a1=np.dot(mesh_size_tab,mesh_size_tab)
 125     a2=np.sum(mesh_size_tab)
 126     a3=nbMeshes
 127
 128     det=a1*a3-a2*a2
 129     assert det!=0, 'test_validation2DWaveSystemSourceCenteredCheckerboard() : Make sure you use distinct meshes and at least two meshes'
 130
 131     b1p=np.dot(error_p_tab,mesh_size_tab)
 132     b2p=np.sum(error_p_tab)
 133     ap=( a3*b1p-a2*b2p)/det
 134     bp=(-a2*b1p+a1*b2p)/det
 135
 136     if(scaling==0):
 137         print( "FV Centered on 2D checkerboard meshes : scheme order for pressure without scaling is ", -ap)
 138     else:
 139         print( "FV Centered on 2D checkerboard meshes : scheme order for pressure with    scaling is ", -ap)
 140
 141     b1u=np.dot(error_u_tab,mesh_size_tab)
 142     b2u=np.sum(error_u_tab)
 143     au=( a3*b1u-a2*b2u)/det
 144     bu=(-a2*b1u+a1*b2u)/det
 145
 146     if(scaling==0):
 147         print( "FV Centered on 2D checkerboard meshes : scheme order for velocity without scaling is ", -au)
 148     else:
 149         print( "FV Centered on 2D checkerboard meshes : scheme order for velocity with    scaling is ", -au)
 150
 151     # Plot of convergence curves
 152     plt.close()
 153     if(scaling==0):
 154         plt.plot(mesh_size_tab, error_p_tab, label='|error on stationary pressure| - no scaling')
 155     else:
 156         plt.plot(mesh_size_tab, error_p_tab, label='|error on stationary pressure| - with scaling')
 157     plt.legend()
 158     plt.xlabel('1/2 log(Number of cells)')
 159     plt.ylabel('log(Error p)')
 160     plt.title('Convergence of finite volumes for the stationary Wave System \n with centered scheme on 2D checkerboard meshes')
 161     plt.savefig(mesh_name+"_Pressure_2DWaveSystemSourceCentered_"+"scaling"+str(scaling)+"_ConvergenceCurve.png")
 162
 163     plt.close()
 164     if(scaling==0):
 165         plt.plot(mesh_size_tab, error_u_tab, label='log(|error on stationary velocity|) - no scaling')
 166     else:
 167         plt.plot(mesh_size_tab, error_u_tab, label='log(|error on stationary velocity|) - with scaling')
 168     plt.legend()
 169     plt.xlabel('1/2 log(Number of cells)')
 170     plt.ylabel('log(Error u)')
 171     plt.title('Convergence of finite volumes for the stationary Wave System \n with centered scheme on 2D checkerboard meshes')
 172     plt.savefig(mesh_name+"_Velocity_2DWaveSystemSourceCentered_"+"scaling"+str(scaling)+"_ConvergenceCurve.png")
 173
 174     # Plot of computational time
 175     plt.close()
 176     if(scaling==0):
 177         plt.plot(mesh_size_tab, time_tab, label='log(cpu time) - no scaling')
 178     else:
 179         plt.plot(mesh_size_tab, time_tab, label='log(cpu time) - with scaling')
 180     plt.legend()
 181     plt.xlabel('1/2 log(Number of cells)')
 182     plt.ylabel('log(cpu time)')
 183     plt.title('Computational time of finite volumes for the stationary Wave System \n with centered scheme on 2D checkerboard meshes')
 184     plt.savefig(mesh_name+"_2DWaveSystemSourceCentered_"+"scaling"+str(scaling)+"_ComputationalTime.png")
 185
 186     plt.close('all')
 187
 188     convergence_synthesis={}
 189
 190     convergence_synthesis["PDE_model"]="Wave system"
 191     convergence_synthesis["PDE_is_stationary"]=False
 192     convergence_synthesis["PDE_search_for_stationary_solution"]=True
 193     convergence_synthesis["Numerical_method_name"]="Centered"
 194     convergence_synthesis["Numerical_method_space_discretization"]="Finite volumes"
 195     convergence_synthesis["Numerical_method_time_discretization"]="Implicit"
 196     convergence_synthesis["Initial_data"]="Constant pressure, divergence free velocity"
 197     convergence_synthesis["Boundary_conditions"]="Periodic"
 198     convergence_synthesis["Numerical_parameter_cfl"]=cfl
 199     convergence_synthesis["Space_dimension"]=2
 200     convergence_synthesis["Mesh_dimension"]=2
 201     convergence_synthesis["Mesh_names"]=meshList
 202     convergence_synthesis["Mesh_type"]=meshType
 203     convergence_synthesis["Mesh_path"]=mesh_path
 204     convergence_synthesis["Mesh_description"]=mesh_name
 205     convergence_synthesis["Mesh_sizes"]=mesh_size_tab
 206     convergence_synthesis["Mesh_cell_type"]="Squares"
 207     convergence_synthesis["Numerical_error_velocity"]=error_u_tab
 208     convergence_synthesis["Numerical_error_pressure"]=error_p_tab
 209     convergence_synthesis["Max_vel_norm"]=max_vel
 210     convergence_synthesis["Final_time"]=t_final
 211     convergence_synthesis["Final_time_step"]=ndt_final
 212     convergence_synthesis["Scheme_order"]=min(-au,-ap)
 213     convergence_synthesis["Scheme_order_vel"]=-au
 214     convergence_synthesis["Scheme_order_press"]=-ap
 215     convergence_synthesis["Scaling_preconditioner"]="None"
 216     convergence_synthesis["Test_color"]=testColor
 217     convergence_synthesis["Computational_time"]=end-start
 218
 219     with open('Convergence_WaveSystemSource_2DFV_Centered_'+mesh_name+'.json', 'w') as outfile:
 220         json.dump(convergence_synthesis, outfile)
 221
 222 if __name__ == """__main__""":
 223     if len(sys.argv) >1 :
 224         scaling = int(sys.argv[1])
 225         test_validation2DWaveSystemSourceCenteredCheckerboard(scaling)
 226     else :
 227         test_validation2DWaveSystemSourceCenteredCheckerboard(2)