CDMATH/tests/examples/Poisson2DVF/FiniteVolumes2DPoisson_SQUARE.py

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