CDMATH/tests/examples/Poisson1DVF/FiniteVolumes1DPoisson.py

   1 # -*-coding:utf-8 -*
   2 #===============================================================================================================================
   3 # Name        : Résolution VF de l'équation de Poisson 1D -\triangle u = f avec conditions aux limites de Dirichlet u=0
   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 1D 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)
   9 #================================================================================================================================
  10
  11 import cdmath
  12 from math import sin, pi
  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 uniform mesh
  20     my_mesh = cdmath.Mesh(sys.argv[1])
  21 else :   #rectangular mesh
  22 #Création d'un maillage uniforme du segment [0,1], définition des bords
  23 #======================================================================
  24         nx=100
  25         my_mesh = cdmath.Mesh(0,1,nx)
  26
  27 if( my_mesh.getSpaceDimension()!=1 or my_mesh.getMeshDimension()!=1) :
  28     raise ValueError("Wrong space or mesh dimension : space and mesh dimensions should be 1")
  29
  30 eps=1e-6
  31 my_mesh.setGroupAtPlan(0.,0,eps,"DirichletBorder")#Bord GAUCHE
  32 my_mesh.setGroupAtPlan(1.,0,eps,"DirichletBorder")#Bord DROIT
  33
  34 nbCells = my_mesh.getNumberOfCells()
  35
  36 print( "Mesh loading/building done")
  37 print( "Number of cells  = ", nbCells)
  38
  39 #Discrétisation du second membre et extraction du nb max de voisins d'une cellule
  40 #================================================================================
  41 my_RHSfield = cdmath.Field("RHS_field", cdmath.CELLS, my_mesh, 1)
  42 maxNbNeighbours=0#This is to determine the number of non zero coefficients in the sparse finite element rigidity matrix
  43 #parcours des cellules pour discrétisation du second membre et extraction du nb max de voisins d'une cellule
  44 for i in range(nbCells):
  45         Ci = my_mesh.getCell(i)
  46         x = Ci.x()
  47
  48         my_RHSfield[i]=pi*pi*sin(pi*x)#mettre la fonction definie au second membre 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( "Max nb 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                         RHS[i]+=coeff*sin(pi*x)#mettre ici  la solution exacte de l'edp
  78                 Rigidite.addValue(i,i,coeff) # terme diagonal
  79
  80 print("Linear system matrix building done")
  81
  82 # Résolution du système linéaire
  83 #=================================
  84 LS=cdmath.LinearSolver(Rigidite,RHS,100,1.E-6,"GMRES","ILU")
  85 SolSyst=LS.solve()
  86
  87 print( "Preconditioner used : ", LS.getNameOfPc() )
  88 print( "Number of iterations used : ", LS.getNumberOfIter() )
  89 print( "Final residual : ", LS.getResidu() )
  90 print( "Linear system solved")
  91
  92 # Création du champ résultat
  93 #===========================
  94 my_ResultField = cdmath.Field("ResultField", cdmath.CELLS, my_mesh, 1)
  95 for i in range(nbCells):
  96     my_ResultField[i]=SolSyst[i];
  97 #sauvegarde sur le disque dur du résultat dans un fichier paraview
  98 my_ResultField.writeVTK("FiniteVolumes1DPoisson_ResultField")
  99
 100 #Postprocessing :
 101 #===============
 102 # save 1D picture
 103 PV_routines.Save_PV_data_to_picture_file("FiniteVolumes1DPoisson_ResultField"+'_0.vtu',"ResultField",'CELLS',"FiniteVolumes1DPoisson_ResultField")
 104
 105 # extract and plot diagonal values
 106 resolution=100
 107 curv_abs=linspace(0,1,resolution+1)
 108 diag_data=VTK_routines.Extract_field_data_over_line_to_numpyArray(my_ResultField,[0,0,0],[1,0,0], resolution)
 109 plt.legend()
 110 plt.xlabel('Position')
 111 plt.ylabel('Value')
 112 if len(sys.argv) >1 :
 113     plt.title('Finite Volumes \n for Laplace operator in 1D'+my_mesh.getName())
 114     plt.plot(curv_abs, diag_data, label= str(nbCells)+ ' cells mesh')
 115     plt.savefig("FiniteVolumes1DPoisson_ResultField_"+str(nbCells)+ '_cells'+".png")
 116 else :
 117     plt.title('Finite Volumes \n for Laplace operator on a 1D regular grid')
 118     plt.plot(curv_abs, diag_data, label= str(nx) + ' cells mesh')
 119     plt.savefig("FiniteVolumes1DPoisson_ResultField_"+str(nx) + '_cells'+".png")
 120
 121 print("Numerical solution of 1D Poisson equation using finite volumes done")
 122
 123 #Calcul de l'erreur commise par rapport à la solution exacte
 124 #===========================================================
 125 #The following formulas use the fact that the exact solution is equal the right hand side divided by pi*pi
 126 max_abs_sol_exacte=max(my_RHSfield.max(),-my_RHSfield.min())/(pi*pi)
 127 max_sol_num=my_ResultField.max()
 128 min_sol_num=my_ResultField.min()
 129 erreur_abs=0
 130 for i in range(nbCells) :
 131     if erreur_abs <  abs(my_RHSfield[i]/(pi*pi) - my_ResultField[i]) :
 132         erreur_abs = abs(my_RHSfield[i]/(pi*pi) - my_ResultField[i])
 133
 134 print("Absolute error = max(| exact solution - numerical solution |) = ",erreur_abs )
 135 print("Relative error = max(| exact solution - numerical solution |)/max(| exact solution |) = ",erreur_abs/max_abs_sol_exacte)
 136 print("Maximum numerical solution = ", max_sol_num, " Minimum numerical solution = ", min_sol_num)
 137 print("Maximum exact solution = ", my_RHSfield.max()/(pi*pi), " Minimum exact solution = ", my_RHSfield.min()/(pi*pi))
 138
 139 assert erreur_abs/max_abs_sol_exacte <1.