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 carré 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((x*x-y*y)/(2*x*y))
   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"]="Square"
  30 test_desc["Part_of_mesh_convergence_analysis"]=True
  31
  32 def solve(my_mesh,filename,resolution, meshName, testColor):
  33     start = time.time()
  34     test_desc["Mesh_name"]=meshName
  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((x*x-y*y)/(2*x*y))
  63         if x<-eps or x>1+eps or y<-eps or y>1+eps :
  64             print("!!! Warning Mesh ", meshName," !!! Cell is not in the unit square.",", eps=",eps, ", x= ",x, ", y= ",y)
  65             #raise ValueError("!!! Domain should be the unit square.")
  66         if abs(x*y) > eps :
  67             my_ExactSol[i] = atan((x*x-y*y)/(2*x*y))
  68         elif x**2-y**2>0 :
  69             my_ExactSol[i] = pi/2
  70         elif x**2-y**2<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                 #Robust calculation of atan((x*x-y*y)/(2*x*y))
 107                 if x<-eps or x>1+eps or y<-eps or y>1+eps :
 108                     print("!!! Warning Mesh ", meshName," !!! Node is not in the unit square.",", eps=",eps, ", x= ",x, ", y= ",y)
 109                     #raise ValueError("!!! Domain should be the unit square.")
 110                 if abs(x*y) > eps :
 111                     RHS[i]+= atan((x*x-y*y)/(2*x*y))*coeff
 112                 elif x**2-y**2>0 :
 113                     RHS[i]+= pi/2*coeff
 114                 elif x**2-y**2<0 :
 115                     RHS[i]+= -pi/2*coeff
 116                 else : #x=0
 117                     RHS[i]+= 0
 118             Rigidite.addValue(i,i,coeff) # terme diagonal
 119
 120     print("Linear system matrix building done")
 121
 122     # Résolution du système linéaire
 123     #=================================
 124     LS=cdmath.LinearSolver(Rigidite,RHS,100,1.E-6,"CG","LU")
 125     LS.setComputeConditionNumber()
 126     SolSyst=LS.solve()
 127
 128     print( "Preconditioner used : ", LS.getNameOfPc() )
 129     print( "Number of iterations used : ", LS.getNumberOfIter() )
 130     print("Linear system solved")
 131
 132     test_desc["Linear_solver_algorithm"]=LS.getNameOfMethod()
 133     test_desc["Linear_solver_preconditioner"]=LS.getNameOfPc()
 134     test_desc["Linear_solver_precision"]=LS.getTolerance()
 135     test_desc["Linear_solver_maximum_iterations"]=LS.getNumberMaxOfIter()
 136     test_desc["Linear_system_max_actual_iterations_number"]=LS.getNumberOfIter()
 137     test_desc["Linear_system_max_actual_error"]=LS.getResidu()
 138     test_desc["Linear_system_max_actual_condition number"]=LS.getConditionNumber()
 139
 140     # Création du champ résultat
 141     #===========================
 142     my_ResultField = cdmath.Field("ResultField", cdmath.CELLS, my_mesh, 1)
 143     for i in range(nbCells):
 144         my_ResultField[i]=SolSyst[i];
 145     #sauvegarde sur le disque dur du résultat dans un fichier paraview
 146     my_ResultField.writeVTK("FiniteVolumes2DPoissonStiffBC_SQUARE_"+meshName+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_SQUARE_"+meshName+str(nbCells)+'_0.vtu',"ResultField",'CELLS',"FiniteVolumes2DPoissonStiffBC_SQUARE_"+meshName+str(nbCells))
 167
 168     test_desc["Computational_time_taken_by_run"]=end-start
 169     test_desc["Absolute_error"]=l2_error
 170     test_desc["Relative_error"]=l2_error/l2_norm_sol_exacte
 171
 172     with open('test_PoissonStiffBC'+str(my_mesh.getMeshDimension())+'D_VF_'+"SQUARE_"+str(nbCells)+ "Cells.json", 'w') as outfile:
 173         json.dump(test_desc, outfile)
 174
 175     return l2_error/l2_norm_sol_exacte, nbCells, diag_data, my_ResultField.min(), my_ResultField.max(), end - start
 176
 177
 178 def solve_file( filename,resolution, meshName, testColor):
 179     my_mesh = cdmath.Mesh(filename+".med")
 180     return solve(my_mesh, filename,resolution, meshName, testColor)
 181
 182 if __name__ == """__main__""":
 183     solve("diskWithTriangles",100,"Unstructured_triangles","Green")