CDMATH/tests/examples/Poisson2DVF_DISK_StiffBC/FiniteVolumes2DPoisson_DISK_StiffBC.py

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