CoreFlows/examples/Python/Convergence/StationaryDiffusion/2DFV_Dirichlet_Hexagons/convergence_StationaryDiffusion_2DFV_Dirichlet_Hexagons.py

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