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