CoreFlows/examples/Python/StationaryDiffusionEquation_2DEF_Neumann.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 Neumann 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=-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_2DEF_StructuredTriangles_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=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 boundary condition for each boundary
  41         myProblem.setNeumannBoundaryCondition("Bord1")
  42         myProblem.setNeumannBoundaryCondition("Bord2")
  43         myProblem.setNeumannBoundaryCondition("Bord3")
  44         myProblem.setNeumannBoundaryCondition("Bord4")
  45
  46         #Set the right hand side function
  47         my_RHSfield = cm.Field("RHS_field", cm.NODES, M, 1)
  48         for i in range(M.getNumberOfNodes()):
  49                 Ni= M.getNode(i)
  50                 x = Ni.x()
  51                 y = Ni.y()
  52
  53                 my_RHSfield[i]=2*pi*pi*cos(pi*x)*cos(pi*y)#mettre la fonction definie au second membre de l'edp
  54
  55         myProblem.setHeatPowerField(my_RHSfield)
  56         myProblem.setLinearSolver(cf.GMRES,cf.ILU);
  57
  58     # name of result file
  59         fileName = "StationnaryDiffusion_2DEF_StructuredTriangles_Neumann";
  60
  61     # computation parameters
  62         myProblem.setFileName(fileName);
  63
  64     # Run the computation
  65         myProblem.initialize();
  66         print("Running python "+ fileName );
  67
  68         ok = myProblem.solveStationaryProblem();
  69         if (not ok):
  70                 print( "Python simulation of " + fileName + "  failed ! " );
  71                 pass
  72         else:
  73                 ####################### Postprocessing #########################
  74                 my_ResultField = myProblem.getOutputTemperatureField()
  75                 #The following formulas use the fact that the exact solution is equal the right hand side divided by 2*pi*pi
  76                 max_abs_sol_exacte=max(my_RHSfield.max(),-my_RHSfield.min())/(2*pi*pi)
  77                 max_sol_num=my_ResultField.max()
  78                 min_sol_num=my_ResultField.min()
  79                 erreur_abs=0
  80                 for i in range(M.getNumberOfNodes()) :
  81                         if erreur_abs < abs(my_RHSfield[i]/(2*pi*pi) - my_ResultField[i]) :
  82                                 erreur_abs = abs(my_RHSfield[i]/(2*pi*pi) - my_ResultField[i])
  83
  84                 print("Absolute error = max(| exact solution - numerical solution |) = ",erreur_abs )
  85                 print("Relative error = max(| exact solution - numerical solution |)/max(| exact solution |) = ",erreur_abs/max_abs_sol_exacte)
  86                 print ("Maximum numerical solution = ", max_sol_num, " Minimum numerical solution = ", min_sol_num)
  87
  88                 assert erreur_abs/max_abs_sol_exacte <1.
  89         pass
  90
  91         print( "------------ !!! End of calculation !!! -----------" );
  92
  93         myProblem.terminate();
  94         return ok
  95
  96 if __name__ == """__main__""":
  97     StationaryDiffusionEquation_2DEF_StructuredTriangles_Neumann()