CoreFlows/examples/Python/StationaryDiffusionEquation/StationaryDiffusionEquation_2DFV_EquilateralTriangles.py

   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 Dirichlet 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 de triangules équilatéraux
   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=sin(pi*x)*sin(pi*y)
  10 #================================================================================================================================
  11
  12 import CoreFlows as cf
  13 import cdmath as cm
  14 from math import sin, pi
  15 import os
  16
  17 def StationaryDiffusionEquation_2DFV_EquilateralTriangles():
  18         spaceDim = 2;
  19         # Prepare for the mesh
  20         print("Loading mesh " );
  21         M=cm.Mesh('../resources/squareWithEquilateralTriangles20.med')#Equilateral triangular mesh
  22
  23         print( "Loaded 2D equilateral triangle mesh with ", M.getNumberOfCells(), " cells")
  24
  25         # set the limit field
  26         boundaryFaces = M.getBoundaryFaceIds()
  27         boundaryValues = {}
  28         print("Setting Dirichlet boundary values")
  29         for i in boundaryFaces :
  30                 Fi=M.getFace(i)
  31                 x=Fi.x()
  32                 y=Fi.y()
  33                 boundaryValues[i] = sin(pi*x)*sin(pi*y)
  34
  35         FEComputation=False
  36         myProblem = cf.StationaryDiffusionEquation(spaceDim,FEComputation);
  37         myProblem.setMesh(M);
  38         myProblem.setDirichletValues(cf.MapIntDouble(boundaryValues))
  39
  40         #Set the right hand side function
  41         my_RHSfield = cm.Field("RHS_field", cm.CELLS, M, 1)
  42         for i in range(M.getNumberOfCells()):
  43                 Ci= M.getCell(i)
  44                 x = Ci.x()
  45                 y = Ci.y()
  46
  47                 my_RHSfield[i]=2*pi*pi*sin(pi*x)*sin(pi*y)#mettre la fonction definie au second membre de l'edp
  48
  49         myProblem.setHeatPowerField(my_RHSfield)
  50         myProblem.setLinearSolver(cf.GMRES,cf.ILU);
  51
  52         # name of result file
  53         fileName = "StationnaryDiffusion_2DFV_StructuredTriangles";
  54
  55         # computation parameters
  56         myProblem.setFileName(fileName);
  57
  58         # Run the computation
  59         myProblem.initialize();
  60         print("Running python "+ fileName );
  61
  62         ok = myProblem.solveStationaryProblem();
  63         if (not ok):
  64                 print( "Python simulation of " + fileName + "  failed ! " );
  65                 pass
  66         else:
  67                 print( "Python simulation of " + fileName + " is successful !" );
  68                 ####################### Postprocessing #########################
  69                 my_ResultField = myProblem.getOutputTemperatureField()
  70                 #The following formulas use the fact that the exact solution is equal the right hand side divided by 2*pi*pi
  71                 max_abs_sol_exacte=max(my_RHSfield.max(),-my_RHSfield.min())/(2*pi*pi)
  72                 max_sol_num=my_ResultField.max()
  73                 min_sol_num=my_ResultField.min()
  74                 erreur_abs=0
  75                 for i in range(M.getNumberOfCells()) :
  76                         if erreur_abs < abs(my_RHSfield[i]/(2*pi*pi) - my_ResultField[i]) :
  77                                 erreur_abs = abs(my_RHSfield[i]/(2*pi*pi) - my_ResultField[i])
  78
  79                 print("Absolute error = max(| exact solution - numerical solution |) = ",erreur_abs )
  80                 print("Relative error = max(| exact solution - numerical solution |)/max(| exact solution |) = ",erreur_abs/max_abs_sol_exacte)
  81                 print("Maximum numerical solution = ", max_sol_num, " Minimum numerical solution = ", min_sol_num)
  82
  83                 assert erreur_abs/max_abs_sol_exacte <1.
  84                 pass
  85
  86         print( "------------ !!! End of calculation !!! -----------" );
  87
  88         myProblem.terminate();
  89         return ok
  90
  91 if __name__ == """__main__""":
  92         StationaryDiffusionEquation_2DFV_EquilateralTriangles()