SALOME platform Git repositories - tools/solverlab.git/blob

   1 # -*-coding:utf-8 -*
   2 #===============================================================================================================================
   3 # Name        : Résolution VF de l'équation de Poisson 2D -\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 triangulaire
   7 #                       Création et sauvegarde du champ résultant ainsi que du champ second membre en utilisant la librairie CDMATH
   8 #               Comparaison de la solution numérique avec la solution exacte u=atan(2*x/(x*x+y*y-1))
   9 #================================================================================================================================
  10
  11 import cdmath
  12 from math import atan, pi
  13 import time, json
  14 import PV_routines
  15 import VTK_routines
  16
  17 test_desc={}
  18 test_desc["Initial_data"]="None"
  19 test_desc["Boundary_conditions"]="Dirichlet"
  20 test_desc["Global_name"]="FV simulation of the 2D Poisson equation"
  21 test_desc["Global_comment"]="Condition limite discontinue, Maillage triangulaire"
  22 test_desc["PDE_model"]="Poisson"
  23 test_desc["PDE_is_stationary"]=True
  24 test_desc["PDE_search_for_stationary_solution"]=False
  25 test_desc["Numerical_method_name"]="VF5"
  26 test_desc["Numerical_method_space_discretization"]="Finite elements"
  27 test_desc["Numerical_method_time_discretization"]="None"
  28 test_desc["Mesh_is_unstructured"]=True
  29 test_desc["Geometry"]="Disk"
  30 test_desc["Part_of_mesh_convergence_analysis"]=True
  31
  32 def solve(my_mesh,filename,resolution, meshType, testColor):
  33     start = time.time()
  34     test_desc["Mesh_type"]=meshType
  35     test_desc["Test_color"]=testColor
  36
  37     nbCells = my_mesh.getNumberOfCells()
  38
  39     if( my_mesh.getSpaceDimension()!=2 or my_mesh.getMeshDimension()!=2) :
  40         raise ValueError("Wrong space or mesh dimension : space and mesh dimensions should be 2")
  41
  42     test_desc["Space_dimension"]=my_mesh.getSpaceDimension()
  43     test_desc["Mesh_dimension"]=my_mesh.getMeshDimension()
  44     test_desc["Mesh_number_of_elements"]=my_mesh.getNumberOfCells()
  45     test_desc["Mesh_cell_type"]=my_mesh.getElementTypesNames()
  46
  47     print("Mesh groups done")
  48     print("Number of cells  = ", nbCells)
  49
  50     #Discrétisation du second membre et extraction du nb max de voisins d'une cellule
  51     #================================================================================
  52     my_ExactSol = cdmath.Field("Exact_field", cdmath.CELLS, my_mesh, 1)
  53     maxNbNeighbours=0#This is to determine the number of non zero coefficients in the sparse finite element rigidity matrix
  54     eps=1e-6#For coarse meshes
  55
  56     #parcours des cellules pour discrétisation du second membre et extraction du nb max de voisins d'une cellule
  57     for i in range(nbCells):
  58         Ci = my_mesh.getCell(i)
  59         x = Ci.x()
  60         y = Ci.y()
  61
  62         #Robust calculation of atan(2x/(x**2+y**2-1)
  63         if x**2+y**2-1 > eps :
  64             print("!!! Warning Mesh ",meshType," !!! Cell is not in the unit disk.",", eps=",eps, ", x**2+y**2 - 1=",x**2+y**2 - 1)
  65             #raise ValueError("Exact solution computation !!! Domain should be the unit disk.")
  66         if x**2+y**2-1 < -eps :
  67             my_ExactSol[i] = atan(2*x/(x**2+y**2-1))
  68         elif x>0 : #x**2+y**2-1>=0
  69             my_ExactSol[i] = -pi/2
  70         elif x<0 : #x**2+y**2-1>=0
  71             my_ExactSol[i] =  pi/2
  72         else : #x=0
  73             my_ExactSol[i] = 0
  74
  75         # compute maximum number of neighbours
  76         maxNbNeighbours= max(1+Ci.getNumberOfFaces(),maxNbNeighbours)
  77
  78     test_desc["Mesh_max_number_of_neighbours"]=maxNbNeighbours
  79
  80     print("Right hand side discretisation done")
  81     print("Max nb of neighbours=", maxNbNeighbours)
  82
  83     # Construction de la matrice et du vecteur second membre du système linéaire
  84     #===========================================================================
  85     Rigidite=cdmath.SparseMatrixPetsc(nbCells,nbCells,maxNbNeighbours) # warning : third argument is maximum number of non zero coefficients per line of the matrix
  86     RHS=cdmath.Vector(nbCells)
  87
  88     #Parcours des cellules du domaine
  89     for i in range(nbCells):
  90         Ci=my_mesh.getCell(i)
  91         for j in range(Ci.getNumberOfFaces()):# parcours des faces voisinnes
  92             Fj=my_mesh.getFace(Ci.getFaceId(j))
  93             if not Fj.isBorder():
  94                 k=Fj.getCellId(0)
  95                 if k==i :
  96                     k=Fj.getCellId(1)
  97                 Ck=my_mesh.getCell(k)
  98                 distance=Ci.getBarryCenter().distance(Ck.getBarryCenter())
  99                 coeff=Fj.getMeasure()/Ci.getMeasure()/distance
 100                 Rigidite.setValue(i,k,-coeff) # terme extradiagonal
 101             else:
 102                 coeff=Fj.getMeasure()/Ci.getMeasure()/Ci.getBarryCenter().distance(Fj.getBarryCenter())
 103                 #For the particular case where the mesh boundary does not coincide with the domain boundary
 104                 x=Fj.getBarryCenter().x()
 105                 y=Fj.getBarryCenter().y()
 106                 if x**2+y**2-1 > eps :
 107                     print("!!! Warning Mesh ", meshType," !!! Face is not in the unit disk.",", eps=",eps, ", x**2+y**2 - 1=",x**2+y**2 - 1)
 108                     #raise ValueError("!!! Domain should be the unit disk.")
 109                 if x**2+y**2-1 < -eps :
 110                     RHS[i]+= coeff*atan(2*x/(x**2+y**2-1))
 111                 elif x>0 : #x**2+y**2-1>=0
 112                     RHS[i]+= coeff*(-pi/2)
 113                 elif x<0 : #x**2+y**2-1>=0
 114                     RHS[i]+= coeff*pi/2
 115                 else : #x=0
 116                     RHS[i]+=  0
 117             Rigidite.addValue(i,i,coeff) # terme diagonal
 118
 119     print("Linear system matrix building done")
 120
 121     # Résolution du système linéaire
 122     #=================================
 123     LS=cdmath.LinearSolver(Rigidite,RHS,500,1.E-6,"CG","LU")
 124     LS.setComputeConditionNumber()
 125     SolSyst=LS.solve()
 126
 127     print( "Preconditioner used : ", LS.getNameOfPc() )
 128     print( "Number of iterations used : ", LS.getNumberOfIter() )
 129     print("Linear system solved")
 130
 131     test_desc["Linear_solver_algorithm"]=LS.getNameOfMethod()
 132     test_desc["Linear_solver_preconditioner"]=LS.getNameOfPc()
 133     test_desc["Linear_solver_precision"]=LS.getTolerance()
 134     test_desc["Linear_solver_maximum_iterations"]=LS.getNumberMaxOfIter()
 135     test_desc["Linear_system_max_actual_iterations_number"]=LS.getNumberOfIter()
 136     test_desc["Linear_system_max_actual_error"]=LS.getResidu()
 137     test_desc["Linear_system_max_actual_condition number"]=LS.getConditionNumber()
 138
 139     # Création du champ résultat
 140     #===========================
 141     my_ResultField = cdmath.Field("ResultField", cdmath.CELLS, my_mesh, 1)
 142     for i in range(nbCells):
 143         my_ResultField[i]=SolSyst[i];
 144     #sauvegarde sur le disque dur du résultat dans un fichier paraview
 145     my_ResultField.writeVTK("FiniteVolumes2DPoissonStiffBC_DISK_"+meshType+str(nbCells))
 146     my_ExactSol.writeVTK("ExactSol2DPoissonStiffBC_DISK_"+meshType+str(nbCells))
 147
 148     print("Numerical solution of 2D Poisson equation on a disk using finite volumes done")
 149
 150     end = time.time()
 151
 152     #Calcul de l'erreur commise par rapport à la solution exacte
 153     #===========================================================
 154     l2_norm_sol_exacte=my_ExactSol.normL2()[0]
 155     l2_error = (my_ExactSol - my_ResultField).normL2()[0]
 156
 157     print("L2 relative error = norm( exact solution - numerical solution )/norm( exact solution ) = ", l2_error/l2_norm_sol_exacte)
 158     print("Maximum numerical solution = ", my_ResultField.max(), " Minimum numerical solution = ", my_ResultField.min())
 159     print("Maximum exact solution = ", my_ExactSol.max(), " Minimum exact solution = ", my_ExactSol.min())
 160
 161     #Postprocessing :
 162     #================
 163         # Extraction of the diagonal data
 164     diag_data=VTK_routines.Extract_field_data_over_line_to_numpyArray(my_ResultField,[0,-1,0],[0,1,0], resolution)
 165     # save 2D picture
 166     PV_routines.Save_PV_data_to_picture_file("FiniteVolumes2DPoissonStiffBC_DISK_"+meshType+str(nbCells)+'_0.vtu',"ResultField",'CELLS',"FiniteVolumes2DPoissonStiffBC_DISK_"+meshType+str(nbCells))
 167     PV_routines.Save_PV_data_to_picture_file("ExactSol2DPoissonStiffBC_DISK_"+meshType+str(nbCells)+'_0.vtu',"Exact_field",'CELLS',"ExactSol2DPoissonStiffBC_DISK_"+meshType+str(nbCells))
 168
 169     test_desc["Computational_time_taken_by_run"]=end-start
 170     test_desc["Absolute_error"]=l2_error
 171     test_desc["Relative_error"]=l2_error/l2_norm_sol_exacte
 172
 173     with open('test_PoissonStiffBC'+str(my_mesh.getMeshDimension())+'D_VF_'+"DISK_"+str(nbCells)+ "Cells.json", 'w') as outfile:
 174         json.dump(test_desc, outfile)
 175
 176     return l2_error/l2_norm_sol_exacte, nbCells, diag_data, my_ResultField.min(), my_ResultField.max(), end - start
 177
 178
 179 def solve_file( filename,resolution, meshType, testColor):
 180     my_mesh = cdmath.Mesh(filename+".med")
 181     return solve(my_mesh, filename,resolution, meshType, testColor)
 182
 183 if __name__ == """__main__""":
 184     solve("diskWithTriangles",100,"Unstructured_triangles","Green")