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

   1 import TransportEquation1DUpwindExplicit
   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
   9 def test_validation1DTransportEquationUpwindExplicit(cfl,isSmooth):
  10     start = time.time()
  11     #### 1D regular grid
  12     meshList=[10,20,50,100,200, 400]
  13     meshType="1D regular grid"
  14     testColor="Green"
  15     nbMeshes=len(meshList)
  16     mesh_size_tab=meshList
  17     mesh_name='RegularGrid'
  18
  19     a=0.  ;  b=1.
  20     error_u_tab=[0]*nbMeshes
  21     sol_u=[0]*nbMeshes
  22     total_var_u=[0]*nbMeshes
  23     min_u=[0]*nbMeshes
  24     max_u=[0]*nbMeshes
  25     time_tab=[0]*nbMeshes
  26
  27     plt.close('all')
  28     i=0
  29
  30     # Storing of numerical errors, mesh sizes and solution
  31     for nx in meshList:
  32         min_u[i], max_u[i], sol_u[i], total_var_u[i], error_u_tab[i], time_tab[i] = TransportEquation1DUpwindExplicit.solve(nx,cfl,a,b,isSmooth)
  33         assert max_u[i]>-1e-5 and max_u[i]<1+1e-5
  34         error_u_tab[i]=log10(error_u_tab[i])
  35         time_tab[i]=log10(time_tab[i])
  36         i=i+1
  37
  38     end = time.time()
  39
  40     # Plot of solution
  41     plt.close()
  42     for i in range(nbMeshes):
  43         plt.plot(np.linspace(a,b,meshList[i]), sol_u[i],   label= str(mesh_size_tab[i]) + ' cells')
  44     plt.xlabel('x')
  45     plt.ylabel('u')
  46     plt.title('Plot of the numerical solution of the transport equation \n with explicit upwind scheme on a 1D regular grid')
  47     plt.ylim( - 0.1, 1.1 )
  48     plt.legend()
  49     if(isSmooth):
  50         plt.savefig(mesh_name+"_1DTransportEquation_UpwindExplicit_CFL"+str(cfl)+"_Smooth_PlotOfSolution.png")
  51     else:
  52         plt.savefig(mesh_name+"_1DTransportEquation_UpwindExplicit_CFL"+str(cfl)+"_Stiff_PlotOfSolution.png")
  53     plt.close()
  54
  55     # Plot of maximal value
  56     plt.close()
  57     plt.plot(mesh_size_tab, max_u, label='Maximum value')
  58     plt.legend()
  59     plt.xlabel('Number of cells')
  60     plt.ylabel('Max |u|')
  61     plt.title('Maximum solution norm of the transport equation \n with explicit upwind scheme on a 1D regular grid')
  62     if(isSmooth):
  63         plt.savefig(mesh_name+"_1DTransportEquationUpwindExplicit_CFL"+str(cfl)+"_Smooth_MaxSolution.png")
  64     else:
  65         plt.savefig(mesh_name+"_1DTransportEquationUpwindExplicit_CFL"+str(cfl)+"_Stiff_MaxSolution.png")
  66
  67     # Plot of total variation
  68     plt.close()
  69     plt.plot(mesh_size_tab, total_var_u, label='Total variation')
  70     plt.legend()
  71     plt.xlabel('Number of cells')
  72     plt.ylabel('Var(u)')
  73     plt.title('Total variation for the transport equation \n with explicit upwind scheme on a 1D regular grid')
  74     if(isSmooth):
  75         plt.savefig(mesh_name+"_1DTransportEquationUpwindExplicit_CFL"+str(cfl)+"_Smooth_TotalVariation.png")
  76     else:
  77         plt.savefig(mesh_name+"_1DTransportEquationUpwindExplicit_CFL"+str(cfl)+"_Stiff_TotalVariation.png")
  78
  79     for i in range(nbMeshes):
  80         mesh_size_tab[i]=log10(mesh_size_tab[i])
  81
  82     # Least square linear regression
  83     # Find the best a,b such that f(x)=ax+b best approximates the convergence curve
  84     # The vector X=(a,b) solves a symmetric linear system AX=B with A=(a1,a2\\a2,a3), B=(b1,b2)
  85     a1=np.dot(mesh_size_tab,mesh_size_tab)
  86     a2=np.sum(mesh_size_tab)
  87     a3=nbMeshes
  88
  89     det=a1*a3-a2*a2
  90     assert det!=0, 'test_validation1DTransportEquationUpwindExplicit() : Make sure you use distinct meshes and at least two meshes'
  91
  92     b1u=np.dot(error_u_tab,mesh_size_tab)
  93     b2u=np.sum(error_u_tab)
  94     a=( a3*b1u-a2*b2u)/det
  95     b=(-a2*b1u+a1*b2u)/det
  96
  97     print("Explicit Upwind scheme for Transport Equation on 1D regular grid : scheme order is ", -a)
  98
  99     assert -a>0.48 and -a<1.02
 100
 101     # Plot of convergence curve
 102     plt.close()
 103     plt.plot(mesh_size_tab, error_u_tab, label='log(|error|)')
 104     plt.plot(mesh_size_tab, a*np.array(mesh_size_tab)+b,label='straight line with slope : '+'%.3f' % a)
 105     plt.legend()
 106     plt.xlabel('log(Number of cells)')
 107     plt.ylabel('log(|error u|)')
 108     plt.title('Convergence of finite volumes for the transport equation \n with explicit upwind scheme on a 1D regular grid')
 109     if(isSmooth):
 110         plt.savefig(mesh_name+"1DTransportEquationUpwindExplicit_CFL"+str(cfl)+"_Smooth_ConvergenceCurve.png")
 111     else:
 112         plt.savefig(mesh_name+"1DTransportEquationUpwindExplicit_CFL"+str(cfl)+"_Stiff_ConvergenceCurve.png")
 113
 114     # Plot of computational time
 115     plt.close()
 116     plt.plot(mesh_size_tab, time_tab, label='log(cpu time)')
 117     plt.legend()
 118     plt.xlabel('log(Number of cells)')
 119     plt.ylabel('log(cpu time)')
 120     plt.title('Computational time of finite volumes for the transport equation \n with explicit upwind scheme on a 1D regular grid')
 121     if(isSmooth):
 122         plt.savefig(mesh_name+"1DTransportEquationUpwindExplicit_CFL"+str(cfl)+"_Smooth_ComputationalTime.png")
 123     else:
 124         plt.savefig(mesh_name+"1DTransportEquationUpwindExplicit_CFL"+str(cfl)+"_Stiff_ComputationalTime.png")
 125
 126     plt.close('all')
 127
 128     convergence_synthesis={}
 129
 130     convergence_synthesis["PDE_model"]="Transport_Equation"
 131     convergence_synthesis["PDE_is_stationary"]=False
 132     convergence_synthesis["PDE_search_for_stationary_solution"]=True
 133     convergence_synthesis["Numerical_method_name"]="Upwind scheme"
 134     convergence_synthesis["Numerical_method_space_discretization"]="Finite volumes"
 135     convergence_synthesis["Numerical_method_time_discretization"]="Explicit"
 136     convergence_synthesis["Initial_data"]="sine"
 137     convergence_synthesis["Boundary_conditions"]="Periodic"
 138     convergence_synthesis["Numerical_parameter_cfl"]=cfl
 139     convergence_synthesis["Space_dimension"]=2
 140     convergence_synthesis["Mesh_dimension"]=2
 141     convergence_synthesis["Mesh_names"]=meshList
 142     convergence_synthesis["Mesh_type"]=meshType
 143     convergence_synthesis["Mesh_description"]=mesh_name
 144     convergence_synthesis["Mesh_sizes"]=mesh_size_tab
 145     convergence_synthesis["Mesh_cell_type"]="1D regular grid"
 146     convergence_synthesis["Numerical_ersolution"]=max_u
 147     convergence_synthesis["Scheme_order"]=-a
 148     convergence_synthesis["Test_color"]=testColor
 149     convergence_synthesis["Computational_time"]=end-start
 150
 151     with open('Convergence_1DTransportEquationUpwindExplicit_'+mesh_name+'.json', 'w') as outfile:
 152         json.dump(convergence_synthesis, outfile)
 153
 154 if __name__ == """__main__""":
 155     if len(sys.argv) >2 :
 156         cfl = float(sys.argv[1])
 157         isSmooth = bool(int(sys.argv[2]))
 158         test_validation1DTransportEquationUpwindExplicit(cfl,isSmooth)
 159     else :
 160         test_validation1DTransportEquationUpwindExplicit(0.99,True)
 161