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

   1 #!/usr/bin/env python
   2 # -*-coding:utf-8 -*
   3
   4 import validationStationaryDiffusionEquation
   5 import cdmath as cm
   6 import matplotlib.pyplot as plt
   7 import numpy as np
   8 from math import log10, sqrt
   9 import time, json
  10
  11 convergence_synthesis=dict(validationStationaryDiffusionEquation.test_desc)
  12
  13 def convergence_StationaryDiffusion_2DFV_Neumann_Squares():
  14     start = time.time()
  15     ### 2D FV Square mesh
  16     method='FV'
  17     BC = 'Neumann'
  18     meshList=[5,20,50,100, 200,400]
  19     nbMeshes=len(meshList)
  20     error_tab=[0]*nbMeshes
  21     mesh_size_tab=[0]*nbMeshes
  22     mesh_name='squareWithSquares'
  23     meshType="RegularSquares"
  24     diag_data=[0]*nbMeshes
  25     time_tab=[0]*nbMeshes
  26     resolution=100
  27     curv_abs=np.linspace(0,sqrt(2),resolution+1)
  28     plt.close('all')
  29     i=0
  30     testColor="Green, despite singular matrix"
  31     # Storing of numerical errors, mesh sizes and diagonal values
  32     for nx in meshList:
  33                 my_mesh=cm.Mesh(0,1,nx,0,1,nx)
  34                 error_tab[i], mesh_size_tab[i], diag_data[i], min_sol_num, max_sol_num, time_tab[i] =validationStationaryDiffusionEquation.SolveStationaryDiffusionEquation(my_mesh,resolution,meshType,method,BC)
  35                 assert min_sol_num>     -1.
  36                 assert max_sol_num<1.
  37                 plt.plot(curv_abs, diag_data[i], label= str(mesh_size_tab[i]) + ' cells')
  38                 error_tab[i]=log10(error_tab[i])
  39                 time_tab[i]=log10(time_tab[i])
  40                 mesh_size_tab[i] = 0.5*log10(mesh_size_tab[i])
  41                 i=i+1
  42
  43     end = time.time()
  44
  45     # Plot over diagonal line
  46     plt.legend()
  47     plt.xlabel('Position on diagonal line')
  48     plt.ylabel('Value on diagonal line')
  49     plt.title('Plot over diagonal line for finite volumes with Neumann BC \n for Poisson problem on 2D square meshes')
  50     plt.savefig(mesh_name+"_2DFV_StatDiffusion_Neumann_PlotOverDiagonalLine.png")
  51
  52     # Least square linear regression
  53     # Find the best a,b such that f(x)=ax+b best approximates the convergence curve
  54     # The vector X=(a,b) solves a symmetric linear system AX=B with A=(a1,a2\\a2,a3), B=(b1,b2)
  55     a1=np.dot(mesh_size_tab,mesh_size_tab)
  56     a2=np.sum(mesh_size_tab)
  57     a3=nbMeshes
  58     b1=np.dot(error_tab,mesh_size_tab)
  59     b2=np.sum(error_tab)
  60
  61     det=a1*a3-a2*a2
  62     assert det!=0, 'convergence_StationaryDiffusion_2DFV_Neumann_Squares() : Make sure you use distinct meshes and at least two meshes'
  63     a=( a3*b1-a2*b2)/det
  64     b=(-a2*b1+a1*b2)/det
  65
  66     print "FV for diffusion on 2D square meshes: scheme order is ", -a
  67     assert abs(a+2)<0.1
  68
  69     # Plot of convergence curve
  70     plt.close()
  71     plt.plot(mesh_size_tab, error_tab, label='log(|numerical-exact|)')
  72     plt.plot(mesh_size_tab, a*np.array(mesh_size_tab)+b,label='least square slope : '+'%.3f' % a)
  73     plt.legend()
  74     plt.plot(mesh_size_tab, error_tab)
  75     plt.xlabel('log(sqrt(number of cells))')
  76     plt.ylabel('log(error)')
  77     plt.title('Convergence of finite volumes with Neumann BC \n for Poisson problem  on 2D square meshes')
  78     plt.savefig(mesh_name+"_2DFV_StatDiffusion_Neumann_ConvergenceCurve.png")
  79
  80     # Plot of computational time
  81     plt.close()
  82     plt.plot(mesh_size_tab, time_tab, label='log(cpu time)')
  83     plt.legend()
  84     plt.xlabel('log(sqrt(number of cells))')
  85     plt.ylabel('log(cpu time)')
  86     plt.title('Computational time of finite volumes with Neumann BC \n for Poisson problem on 2D square StrucMeshes')
  87     plt.savefig(mesh_name+"_2DFV_StatDiffusion_Neumann_ComputationalTime.png")
  88
  89     plt.close('all')
  90
  91     convergence_synthesis["Mesh_names"]=meshList
  92     convergence_synthesis["Mesh_type"]=meshType
  93     #convergence_synthesis["Mesh_path"]=mesh_path
  94     convergence_synthesis["Mesh_description"]=mesh_name
  95     convergence_synthesis["Mesh_sizes"]=[10**x for x in mesh_size_tab]
  96     convergence_synthesis["Space_dim"]=2
  97     convergence_synthesis["Mesh_dim"]=2
  98     convergence_synthesis["Mesh_cell_type"]="Squares"
  99     convergence_synthesis["Errors"]=[10**x for x in error_tab]
 100     convergence_synthesis["Scheme_order"]=round(-a,4)
 101     convergence_synthesis["Test_color"]=testColor
 102     convergence_synthesis["PDE_model"]='Poisson'
 103     convergence_synthesis["Num_method"]=method
 104     convergence_synthesis["Bound_cond"]=BC
 105     convergence_synthesis["Comput_time"]=round(end-start,3)
 106
 107     with open('Convergence_Poisson_2DFV_'+mesh_name+'.json', 'w') as outfile:
 108         json.dump(convergence_synthesis, outfile)
 109 if __name__ == """__main__""":
 110         convergence_StationaryDiffusion_2DFV_Neumann_Squares()