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

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