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_3DFE_Dirichlet_RegularTetrahedra():
  14     start = time.time()
  15     ### 3D FE Regular Tetrahedra meshes
  16     method = 'FE'
  17     BC = 'Dirichlet'
  18     meshList=[5,20,50,70]
  19     mesh_name='cubeWithRegularTetrahedra'
  20     meshType="Regular_Tetrahedra"
  21     nbMeshes=len(meshList)
  22     error_tab=[0]*nbMeshes
  23     mesh_size_tab=[0]*nbMeshes
  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="Orange (not order 2)"
  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,0,1,nx,1)
  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
  36                 assert min_sol_num>=0.
  37                 assert max_sol_num<1.
  38                 plt.plot(curv_abs, diag_data[i], label= str(mesh_size_tab[i]) + ' cells')
  39                 error_tab[i]=log10(error_tab[i])
  40                 time_tab[i]=log10(time_tab[i])
  41                 mesh_size_tab[i] = 0.5*log10(mesh_size_tab[i])
  42                 i=i+1
  43     end = time.time()
  44
  45
  46     # Plot over diagonal line
  47     plt.legend()
  48     plt.xlabel('Position on diagonal line')
  49     plt.ylabel('Value on diagonal line')
  50     plt.title('Plot over diagonal line for finite elements for Poisson problem \n on 3D regular tetrahedra with Dirichlet BC')
  51     plt.savefig(mesh_name+"_3DFE_StatDiffusion_Dirichlet_PlotOverDiagonalLine.png")
  52
  53     # Least square linear regression
  54     # Find the best a,b such that f(x)=ax+b best approximates the convergence curve
  55     # The vector X=(a,b) solves a symmetric linear system AX=B with A=(a1,a2\\a2,a3), B=(b1,b2)
  56     a1=np.dot(mesh_size_tab,mesh_size_tab)
  57     a2=np.sum(mesh_size_tab)
  58     a3=nbMeshes
  59     b1=np.dot(error_tab,mesh_size_tab)
  60     b2=np.sum(error_tab)
  61
  62     det=a1*a3-a2*a2
  63     assert det!=0, 'convergence_StationaryDiffusion_3DFE_Dirichlet_regularTetrahedra() : Make sure you use distinct meshes and at least two meshes'
  64     a=( a3*b1-a2*b2)/det
  65     b=(-a2*b1+a1*b2)/det
  66
  67     print "FE for diffusion on 3D regular tetrahedron meshes: scheme order is ", -a
  68     assert abs(a+1.34)<0.001
  69
  70     # Plot of convergence curve
  71     plt.close()
  72     plt.plot(mesh_size_tab, error_tab, label='log(|numerical-exact|)')
  73     plt.plot(mesh_size_tab, a*np.array(mesh_size_tab)+b,label='least square slope : '+'%.3f' % a)
  74     plt.legend()
  75     plt.plot(mesh_size_tab, error_tab)
  76     plt.xlabel('log(sqrt(number of cells))')
  77     plt.ylabel('log(error)')
  78     plt.title('Convergence of finite elements for Poisson problem \n on 3D regular tetrahedra meshes with Dirichlet BC')
  79     plt.savefig(mesh_name+"_3DFE_StatDiffusion_Dirichlet_ConvergenceCurve.png")
  80
  81     # Plot of computational time
  82     plt.close()
  83     plt.plot(mesh_size_tab, time_tab, label='log(cpu time)')
  84     plt.legend()
  85     plt.xlabel('log(sqrt(number of cells))')
  86     plt.ylabel('log(cpu time)')
  87     plt.title('Computational time of finite elements for Poisson problem \n on 3D regular tetrahedra meshes with Dirichlet BC')
  88     plt.savefig(mesh_name+"_3DFE_StatDiffusion_Dirichlet_ComputationalTime.png")
  89
  90     plt.close('all')
  91
  92     convergence_synthesis["Mesh_names"]=meshList
  93     convergence_synthesis["Mesh_type"]=meshType
  94     #convergence_synthesis["Mesh_path"]=mesh_path
  95     convergence_synthesis["Mesh_description"]=mesh_name
  96     convergence_synthesis["Mesh_sizes"]=[10**x for x in mesh_size_tab]
  97     convergence_synthesis["Space_dim"]=3
  98     convergence_synthesis["Mesh_dim"]=3
  99     convergence_synthesis["Mesh_cell_type"]="Tetrahedron"
 100     convergence_synthesis["Errors"]=[10**x for x in error_tab]
 101     convergence_synthesis["Scheme_order"]=round(-a,4)
 102     convergence_synthesis["Test_color"]=testColor
 103     convergence_synthesis["PDE_model"]='Poisson'
 104     convergence_synthesis["Num_method"]=method
 105     convergence_synthesis["Bound_cond"]=BC
 106     convergence_synthesis["Comput_time"]=round(end-start,3)
 107
 108     with open('Convergence_Poisson_3DFE_'+mesh_name+'.json', 'w') as outfile:
 109         json.dump(convergence_synthesis, outfile)
 110
 111
 112 if __name__ == """__main__""":
 113     convergence_StationaryDiffusion_3DFE_Dirichlet_RegularTetrahedra()