CDMATH/tests/examples/Poisson3DVF/FiniteVolumes3DPoisson_CUBE.py

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