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

   1 # -*-coding:utf-8 -*
   2 #===============================================================================================================================
   3 # Name        : Résolution VF de l'équation de Diffusion anisotrope -(\partial_{xx} u + K \partial_{yy} 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 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 import time, json
  13 from math import sin, pi
  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 Diffusion equation"
  21 test_desc["Global_comment"]="2 points FV diffusion scheme"
  22 test_desc["PDE_model"]="Diffusion"
  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 volumes"
  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, 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     dim=my_mesh.getSpaceDimension()
  42
  43     test_desc["Space_dimension"]=my_mesh.getSpaceDimension()
  44     test_desc["Mesh_dimension"]=my_mesh.getMeshDimension()
  45     test_desc["Mesh_number_of_elements"]=my_mesh.getNumberOfCells()
  46     test_desc["Mesh_cell_type"]=my_mesh.getElementTypes()
  47
  48     print("Mesh groups done")
  49     print("Number of cells  = ", nbCells)
  50
  51     #Discrétisation du second membre et extraction du nb max de voisins d'une cellule
  52     #================================================================================
  53     K=10^4
  54     my_RHSfield = cdmath.Field("RHS_field", cdmath.CELLS, my_mesh, 1)
  55     maxNbNeighbours=0#This is to determine the number of non zero coefficients in the sparse finite element rigidity matrix
  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         my_RHSfield[i]=(1+K)*pi*pi*sin(pi*x)*sin(pi*y)#mettre la fonction definie au second membre de l edp
  63         # compute maximum number of neighbours
  64         maxNbNeighbours= max(1+Ci.getNumberOfFaces(),maxNbNeighbours)
  65
  66     test_desc["Mesh_max_number_of_neighbours"]=maxNbNeighbours
  67
  68     print("Right hand side discretisation done")
  69     print("Max nb 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 maximum number of non zero coefficients per line of the matrix
  74     RHS=cdmath.Vector(nbCells)
  75     normal=cdmath.Vector(dim)
  76     #Parcours des cellules du domaine
  77     for i in range(nbCells):
  78         RHS[i]=my_RHSfield[i] #la valeur moyenne du second membre f dans la cellule i
  79         Ci=my_mesh.getCell(i)
  80         for j in range(Ci.getNumberOfFaces()):# parcours des faces voisinnes
  81             Fj=my_mesh.getFace(Ci.getFaceId(j))
  82             for idim in range(dim) :
  83                 normal[idim] = Ci.getNormalVector(j, idim);#normale sortante
  84             if not Fj.isBorder():
  85                 k=Fj.getCellId(0)
  86                 if k==i :
  87                     k=Fj.getCellId(1)
  88                 Ck=my_mesh.getCell(k)
  89                 distance=Ci.getBarryCenter().distance(Ck.getBarryCenter())
  90                 coeff=Fj.getMeasure()/Ci.getMeasure()/distance*(normal[0]*normal[0] + K*normal[1]*normal[1])
  91                 Rigidite.setValue(i,k,-coeff) # terme extradiagonal
  92             else:
  93                 coeff=Fj.getMeasure()/Ci.getMeasure()/Ci.getBarryCenter().distance(Fj.getBarryCenter())*(normal[0]*normal[0] + K*normal[1]*normal[1])
  94                 #For the particular case where the mesh boundary does not coincide with the domain boundary
  95                 x=Fj.getBarryCenter().x()
  96                 y=Fj.getBarryCenter().y()
  97                 RHS[i]+=coeff*sin(pi*x)*sin(pi*y)#mettre ici la condition limite du problème de Dirichlet
  98             Rigidite.addValue(i,i,coeff) # terme diagonal
  99
 100     print("Linear system matrix building done")
 101
 102     # Résolution du système linéaire
 103     #=================================
 104     LS=cdmath.LinearSolver(Rigidite,RHS,500,1.E-6,"CG","LU")
 105     LS.setComputeConditionNumber()
 106     SolSyst=LS.solve()
 107
 108     print "Preconditioner used : ", LS.getNameOfPc()
 109     print "Number of iterations used : ", LS.getNumberOfIter()
 110     print("Linear system solved")
 111
 112     test_desc["Linear_solver_algorithm"]=LS.getNameOfMethod()
 113     test_desc["Linear_solver_preconditioner"]=LS.getNameOfPc()
 114     test_desc["Linear_solver_precision"]=LS.getTolerance()
 115     test_desc["Linear_solver_maximum_iterations"]=LS.getNumberMaxOfIter()
 116     test_desc["Linear_system_max_actual_iterations_number"]=LS.getNumberOfIter()
 117     test_desc["Linear_system_max_actual_error"]=LS.getResidu()
 118     test_desc["Linear_system_max_actual_condition number"]=LS.getConditionNumber()
 119
 120     # Création du champ résultat
 121     #===========================
 122     my_ResultField = cdmath.Field("ResultField", cdmath.CELLS, my_mesh, 1)
 123     for i in range(nbCells):
 124         my_ResultField[i]=SolSyst[i];
 125     #sauvegarde sur le disque dur du résultat dans un fichier paraview
 126     my_ResultField.writeVTK("FiniteVolumes2DDiffusion_SQUARE_"+meshType+str(nbCells))
 127
 128     print("Numerical solution of 2D Diffusion equation on a square using finite elements done")
 129
 130     end = time.time()
 131
 132     #Calcul de l'erreur commise par rapport à la solution exacte
 133     #===========================================================
 134     #The following formulas use the fact that the exact solution is equal the right hand side divided by (1+K)*pi*pi
 135     max_abs_sol_exacte=max(my_RHSfield.max(),-my_RHSfield.min())/((1+K)*pi*pi)
 136     max_sol_num=my_ResultField.max()
 137     min_sol_num=my_ResultField.min()
 138     erreur_abs=0
 139     for i in range(nbCells) :
 140         if  erreur_abs < abs(my_RHSfield[i]/((1+K)*pi*pi) - my_ResultField[i]) :
 141             erreur_abs = abs(my_RHSfield[i]/((1+K)*pi*pi) - my_ResultField[i])
 142
 143     print("Relative error = max(| exact solution - numerical solution |)/max(| exact solution |) = ",erreur_abs/max_abs_sol_exacte)
 144     print ("Maximum numerical solution = ", max_sol_num, " Minimum numerical solution = ", min_sol_num)
 145
 146     #Postprocessing :
 147     #================
 148         # Extraction of the diagonal data
 149     diag_data=VTK_routines.Extract_field_data_over_line_to_numpyArray(my_ResultField,[0,1,0],[1,0,0], resolution)
 150     # save 2D picture
 151     PV_routines.Save_PV_data_to_picture_file("FiniteVolumes2DDiffusion_SQUARE_"+meshType+str(nbCells)+'_0.vtu',"ResultField",'CELLS',"FiniteVolumes2DDiffusion_SQUARE_"+meshType+str(nbCells))
 152
 153     test_desc["Computational_time_taken_by_run"]=end-start
 154     test_desc["Absolute_error"]=erreur_abs
 155     test_desc["Relative_error"]=erreur_abs/max_abs_sol_exacte
 156
 157     with open('test_Diffusion'+str(my_mesh.getMeshDimension())+'D_VF_'+str(nbCells)+ "Cells.json", 'w') as outfile:
 158         json.dump(test_desc, outfile)
 159
 160     return erreur_abs/max_abs_sol_exacte, nbCells, diag_data, min_sol_num, max_sol_num, end - start
 161
 162
 163 def solve_file( filename,resolution, meshType, testColor):
 164     my_mesh = cdmath.Mesh(filename+".med")
 165     return solve(my_mesh, filename,resolution, meshType, testColor)
 166
 167 if __name__ == """__main__""":
 168         mesh51 = cdmath.Mesh(0,1,51,0,1,51)
 169         solve(mesh51,'51',100,"Regular_squares","Green")