CoreFlows/examples/Python/StationaryDiffusionEquation_2DEF.py

   1 #!/usr/bin/env python
   2 # -*-coding:utf-8 -*
   3 #===============================================================================================================================
   4 # Name        : Résolution EF de l'équation de Poisson 2D -\triangle T = f 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 éléménts finis P1 avec champs T et f discrétisés aux noeuds d'un maillage triangulaire
   8 #                       Création et sauvegarde du champ résultant ainsi que du champ second membre en utilisant la librairie 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
  16 def StationaryDiffusionEquation_2DEF_StructuredTriangles():
  17         spaceDim = 2;
  18         # Prepare for the mesh
  19         print("Building mesh " );
  20         xinf = 0 ;
  21         xsup=1.0;
  22         yinf=0.0;
  23         ysup=1.0;
  24         nx=20;
  25         ny=20;
  26         M=cm.Mesh(xinf,xsup,nx,yinf,ysup,ny,0)#Regular triangular 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 triangular 2D mesh from a square mesh with ", nx,"x" ,ny, " cells"
  35
  36         FEComputation=True
  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.setDirichletBoundaryCondition("Bord1",T1)
  47         myProblem.setDirichletBoundaryCondition("Bord2",T2)
  48         myProblem.setDirichletBoundaryCondition("Bord3",T3)
  49         myProblem.setDirichletBoundaryCondition("Bord4",T4)
  50
  51         #Set the right hand side function
  52         my_RHSfield = cm.Field("RHS_field", cm.NODES, M, 1)
  53         for i in range(M.getNumberOfNodes()):
  54                 Ni= M.getNode(i)
  55                 x = Ni.x()
  56                 y = Ni.y()
  57
  58                 my_RHSfield[i]=2*pi*pi*sin(pi*x)*sin(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_2DEF_StructuredTriangles";
  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                 ####################### Postprocessing #########################
  79                 my_ResultField = myProblem.getOutputTemperatureField()
  80                 #The following formulas use the fact that the exact solution is equal the right hand side divided by 2*pi*pi
  81                 max_abs_sol_exacte=max(my_RHSfield.max(),-my_RHSfield.min())/(2*pi*pi)
  82                 max_sol_num=my_ResultField.max()
  83                 min_sol_num=my_ResultField.min()
  84                 erreur_abs=0
  85                 for i in range(M.getNumberOfNodes()) :
  86                         if erreur_abs < abs(my_RHSfield[i]/(2*pi*pi) - my_ResultField[i]) :
  87                                 erreur_abs = abs(my_RHSfield[i]/(2*pi*pi) - my_ResultField[i])
  88
  89                 print("Absolute error = max(| exact solution - numerical solution |) = ",erreur_abs )
  90                 print("Relative error = max(| exact solution - numerical solution |)/max(| exact solution |) = ",erreur_abs/max_abs_sol_exacte)
  91                 print ("Maximum numerical solution = ", max_sol_num, " Minimum numerical solution = ", min_sol_num)
  92
  93                 assert erreur_abs/max_abs_sol_exacte <1.
  94         pass
  95
  96         print( "------------ !!! End of calculation !!! -----------" );
  97
  98         myProblem.terminate();
  99         return ok
 100
 101 if __name__ == """__main__""":
 102     StationaryDiffusionEquation_2DEF_StructuredTriangles()