CDMATH/tests/examples/Poisson2DVF_DISK/FiniteVolumes2DPoisson_DISK.py

   1 # -*-coding:utf-8 -*
   2 #===============================================================================================================================
   3 # Name        : Résolution VF de l'équation de Poisson -\triangle u = f sur un disque avec conditions aux limites de Dirichlet u=0
   4 # Author      : Michaël Ndjinga
   5 # Copyright   : CEA Saclay 2018
   6 # Description : Utilisation de la méthode des volumes finis avec champs u et f discrétisés aux cellules d'un maillage quelconque
   7 #                               Création et sauvegarde du champ résultant ainsi que du champ second membre en utilisant CDMATH
   8 #               Comparaison de la solution numérique avec la solution exacte u=-(r-1)*r**2*sin(theta)
   9 #================================================================================================================================
  10
  11 import cdmath
  12 from math import sin, sqrt, atan2
  13 import numpy as np
  14 import matplotlib.pyplot as plt
  15 import PV_routines
  16 import VTK_routines
  17 import sys
  18
  19 if len(sys.argv) >1 :#non rectangular mesh
  20     my_mesh = cdmath.Mesh(sys.argv[1])
  21 else :
  22     raise ValueError("Give an input mesh of the disk")
  23
  24 nbCells = my_mesh.getNumberOfCells()
  25
  26 if( my_mesh.getSpaceDimension()!=2 or my_mesh.getMeshDimension()!=2) :
  27     raise ValueError("Wrong space or mesh dimension : space and mesh dimensions should be 2")
  28
  29 print "Mesh loading done"
  30 print "Number of cells  = ", nbCells
  31
  32 #Discrétisation du second membre et extraction du nb max de voisins d'une cellule
  33 #================================================================================
  34 my_RHSfield = cdmath.Field("RHS_field", cdmath.CELLS, my_mesh, 1)
  35 my_ExactSol = cdmath.Field("Exact_field", cdmath.CELLS, my_mesh, 1)
  36
  37 maxNbNeighbours=0#This is to determine the number of non zero coefficients in the sparse finite element rigidity matrix
  38 #parcours des cellules pour discrétisation du second membre et extraction du nb max de voisins d'une cellule
  39 for i in range(nbCells):
  40         Ci = my_mesh.getCell(i)
  41         x = Ci.x()
  42         y = Ci.y()
  43
  44         r=sqrt(x*x+y*y)
  45         theta=atan2(y,x)
  46
  47         my_RHSfield[i]=(8*r-3)*sin(theta)#mettre la fonction definie au second membre de l'edp
  48         my_ExactSol[i]=-(r-1)*r*r*sin(theta)#mettre la solution exacte de l'edp
  49         # compute maximum number of neighbours
  50         maxNbNeighbours= max(1+Ci.getNumberOfFaces(),maxNbNeighbours)
  51
  52 print("Right hand side discretisation done")
  53 print "Maximum number of neighbours = ", maxNbNeighbours
  54
  55 # Construction de la matrice et du vecteur second membre du système linéaire
  56 #===========================================================================
  57 Rigidite=cdmath.SparseMatrixPetsc(nbCells,nbCells,maxNbNeighbours)# warning : third argument is max number of non zero coefficients per line of the matrix
  58 RHS=cdmath.Vector(nbCells)
  59 #Parcours des cellules du domaine
  60 for i in range(nbCells):
  61         RHS[i]=my_RHSfield[i] #la valeur moyenne du second membre f dans la cellule i
  62         Ci=my_mesh.getCell(i)
  63         for j in range(Ci.getNumberOfFaces()):# parcours des faces voisinnes
  64                 Fj=my_mesh.getFace(Ci.getFaceId(j))
  65                 if not Fj.isBorder():
  66                         k=Fj.getCellId(0)
  67                         if k==i :
  68                                 k=Fj.getCellId(1)
  69                         Ck=my_mesh.getCell(k)
  70                         distance=Ci.getBarryCenter().distance(Ck.getBarryCenter())
  71                         coeff=Fj.getMeasure()/Ci.getMeasure()/distance
  72                         Rigidite.setValue(i,k,-coeff) # terme extradiagonal
  73                 else:
  74                         coeff=Fj.getMeasure()/Ci.getMeasure()/Ci.getBarryCenter().distance(Fj.getBarryCenter())
  75                         #For the particular case where the mesh boundary does not coincide with the domain boundary
  76                         x=Fj.getBarryCenter().x()
  77                         y=Fj.getBarryCenter().y()
  78                         r=sqrt(x*x+y*y)
  79                         theta=atan2(y,x)
  80                         RHS[i]+=coeff*(-(r-1)*r*r*sin(theta))#mettre ici la solution exacte de l'edp
  81                 Rigidite.addValue(i,i,coeff) # terme diagonal
  82
  83 print("Linear system matrix building done")
  84
  85 # Résolution du système linéaire
  86 #=================================
  87 LS=cdmath.LinearSolver(Rigidite,RHS,500,1.E-6,"GMRES","ILU")
  88 SolSyst=LS.solve()
  89
  90 print "Preconditioner used : ", LS.getNameOfPc()
  91 print "Number of iterations used : ", LS.getNumberOfIter()
  92 print "Final residual : ", LS.getResidu()
  93 print "Linear system solved"
  94
  95 # Création du champ résultat
  96 #===========================
  97 my_ResultField = cdmath.Field("ResultField", cdmath.CELLS, my_mesh, 1)
  98 for i in range(nbCells):
  99     my_ResultField[i]=SolSyst[i];
 100 #sauvegarde sur le disque dur du résultat dans un fichier paraview
 101 my_ResultField.writeVTK("FiniteVolumes2DPoisson_DISK_square_ResultField")
 102
 103 #Postprocessing :
 104 #===============
 105 # save 2D picture
 106 PV_routines.Save_PV_data_to_picture_file("FiniteVolumes2DPoisson_DISK_square_ResultField"+'_0.vtu',"ResultField",'CELLS',"FiniteVolumes2DPoisson_DISK_square_ResultField")
 107
 108 # extract and plot diagonal values
 109 resolution=100
 110 curv_abs=np.linspace(0,sqrt(2),resolution+1)
 111 diag_data=VTK_routines.Extract_field_data_over_line_to_numpyArray(my_ResultField,[0,1,0],[1,0,0], resolution)
 112 plt.legend()
 113 plt.xlabel('Position on diagonal line')
 114 plt.ylabel('Value on diagonal line')
 115 if len(sys.argv) >1 :
 116     plt.title('Plot over diagonal line for finite Volumes \n for Laplace operator on a 2D disk mesh '+my_mesh.getName())
 117     plt.plot(curv_abs, diag_data, label= str(nbCells)+ ' cells mesh')
 118     plt.savefig("FiniteVolumes2DPoisson_DISK_square_ResultField_"+str(nbCells)+ '_cells'+"_PlotOverDiagonalLine.png")
 119 else :
 120     plt.title('Plot over diagonal line for finite Volumes \n for Laplace operator on a 2D disk with a rectangular grid')
 121     plt.plot(curv_abs, diag_data, label= str(nx) +'x'+str(ny)+ ' cells mesh')
 122     plt.savefig("FiniteVolumes2DPoisson_DISK_square_ResultField_"+str(nx) +'x'+str(ny)+ '_cells'+"_PlotOverDiagonalLine.png")
 123
 124 print("Numerical solution of 2D Poisson equation on a disk using finite volumes done")
 125
 126 #Calcul de l'erreur commise par rapport à la solution exacte
 127 #===========================================================
 128 max_abs_sol_exacte=max(my_ExactSol.max(),-my_ExactSol.min())
 129 max_sol_num=my_ResultField.max()
 130 min_sol_num=my_ResultField.min()
 131 erreur_abs=0
 132 for i in range(nbCells) :
 133     if erreur_abs < abs(my_ExactSol[i] - my_ResultField[i]) :
 134         erreur_abs = abs(my_ExactSol[i] - my_ResultField[i])
 135
 136 print("Absolute error = max(| exact solution - numerical solution |) = ",erreur_abs )
 137 print("Relative error = max(| exact solution - numerical solution |)/max(| exact solution |) = ",erreur_abs/max_abs_sol_exacte)
 138 print("Maximum numerical solution = ", max_sol_num, " Minimum numerical solution = ", min_sol_num)
 139 print("Maximum exact solution = ", my_ExactSol.max(), " Minimum exact solution = ", my_ExactSol.min() )
 140
 141 assert erreur_abs/max_abs_sol_exacte <1.