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

   1 #!/usr/bin/env python
   2 # -*-coding:utf-8 -*
   3 #===============================================================================================================================
   4 # Name        : Résolution VF de l'équation de Poisson 2D -\triangle T = f sur un carré avec conditions aux limites de Neumann T=0
   5 # Author      : Michaël Ndjinga
   6 # Copyright   : CEA Saclay 2019
   7 # Description : Utilisation de la méthode des volumes finis avec champs T et f discrétisés aux cellules d'un maillage triangulaire
   8 #                               Création et sauvegarde du champ résultant ainsi que du champ second membre en utilisant CDMATH
   9 #               Comparaison de la solution numérique avec la solution exacte T=cos(pi*x)*cos(pi*y)
  10 #================================================================================================================================
  11
  12 import CoreFlows as cf
  13 import cdmath as cm
  14 from math import cos, pi
  15
  16 def StationaryDiffusionEquation_2DFV_StructuredSquares_Neumann():
  17         spaceDim = 2;
  18         # Prepare for the mesh
  19         print("Building mesh " );
  20         xinf = 0.0;
  21         xsup = 1.0;
  22         yinf = 0.0;
  23         ysup = 1.0;
  24         nx=30;
  25         ny=30;
  26         M=cm.Mesh(xinf,xsup,nx,yinf,ysup,ny)#Regular square mesh
  27         # set the limit field for each boundary
  28         eps=1e-6;
  29         M.setGroupAtPlan(xsup,0,eps,"Bord1")
  30         M.setGroupAtPlan(xinf,0,eps,"Bord2")
  31         M.setGroupAtPlan(ysup,1,eps,"Bord3")
  32         M.setGroupAtPlan(yinf,1,eps,"Bord4")
  33
  34         print "Built a regular 2D square mesh with ", nx,"x" ,ny, " cells"
  35
  36         FEComputation=False
  37         myProblem = cf.StationaryDiffusionEquation(spaceDim,FEComputation);
  38         myProblem.setMesh(M);
  39
  40         # set the limit value for each boundary
  41         T1=0;
  42         T2=0;
  43         T3=0;
  44         T4=0;
  45
  46         myProblem.setNeumannBoundaryCondition("Bord1")
  47         myProblem.setNeumannBoundaryCondition("Bord2")
  48         myProblem.setNeumannBoundaryCondition("Bord3")
  49         myProblem.setNeumannBoundaryCondition("Bord4")
  50
  51         #Set the right hand side function
  52         my_RHSfield = cm.Field("RHS_field", cm.CELLS, M, 1)
  53         for i in range(M.getNumberOfCells()):
  54                 Ci= M.getCell(i)
  55                 x = Ci.x()
  56                 y = Ci.y()
  57
  58                 my_RHSfield[i]=2*pi*pi*cos(pi*x)*cos(pi*y)#mettre la fonction definie au second membre de l'edp
  59
  60         myProblem.setHeatPowerField(my_RHSfield)
  61         myProblem.setLinearSolver(cf.GMRES,cf.ILU);
  62
  63         # name of result file
  64         fileName = "StationnaryDiffusion_2DFV_StructuredSquares_Neumann";
  65
  66         # computation parameters
  67         myProblem.setFileName(fileName);
  68
  69         # Run the computation
  70         myProblem.initialize();
  71         print("Running python "+ fileName );
  72
  73         ok = myProblem.solveStationaryProblem();
  74         if (not ok):
  75                 print( "Python simulation of " + fileName + "  failed ! " );
  76                 pass
  77         else:
  78                 print( "Python simulation of " + fileName + " is successful !" );
  79                 ####################### Postprocessing #########################
  80                 my_ResultField = myProblem.getOutputTemperatureField()
  81                 #The following formulas use the fact that the exact solution is equal the right hand side divided by 2*pi*pi
  82                 max_abs_sol_exacte=max(my_RHSfield.max(),-my_RHSfield.min())/(2*pi*pi)
  83                 max_sol_num=my_ResultField.max()
  84                 min_sol_num=my_ResultField.min()
  85                 erreur_abs=0
  86                 for i in range(M.getNumberOfCells()) :
  87                         if erreur_abs < abs(my_RHSfield[i]/(2*pi*pi) - my_ResultField[i]) :
  88                                 erreur_abs = abs(my_RHSfield[i]/(2*pi*pi) - my_ResultField[i])
  89
  90                 print("Absolute error = max(| exact solution - numerical solution |) = ",erreur_abs )
  91                 print("Relative error = max(| exact solution - numerical solution |)/max(| exact solution |) = ",erreur_abs/max_abs_sol_exacte)
  92                 print ("Maximum numerical solution = ", max_sol_num, " Minimum numerical solution = ", min_sol_num)
  93
  94                 assert erreur_abs/max_abs_sol_exacte <1.
  95         pass
  96
  97         print( "------------ !!! End of calculation !!! -----------" );
  98
  99         myProblem.terminate();
 100         return ok
 101
 102 if __name__ == """__main__""":
 103         StationaryDiffusionEquation_2DFV_StructuredSquares_Neumann()