CoreFlows/examples/Python/StationaryDiffusionEquation/StationaryDiffusionEquation_3DFV_StructuredTetrahedra.py

   1 #!/usr/bin/env python
   2 # -*-coding:utf-8 -*
   3 #===============================================================================================================================
   4 # Name        : Résolution VF de l'équation de Poisson 3D -\triangle T = f sur un cube 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 tétraédrique
   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)*sin(pi*z)
  10 #================================================================================================================================
  11
  12 import CoreFlows as cf
  13 import cdmath as cm
  14 from math import sin, pi
  15
  16 def StationaryDiffusionEquation_3DFV_StructuredTetrahedra():
  17         spaceDim = 3;
  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         zinf = 0.0;
  25         zsup = 1.0;
  26         nx = 5;
  27         ny = 5;
  28         nz = 5;
  29         M=cm.Mesh(xinf,xsup,nx,yinf,ysup,ny,zinf,zsup,nz,1)#Regular tetrahedral mesh
  30         # set the limit field for each boundary
  31         eps=1e-6;
  32         M.setGroupAtPlan(xsup,0,eps,"Bord1")
  33         M.setGroupAtPlan(xinf,0,eps,"Bord2")
  34         M.setGroupAtPlan(ysup,1,eps,"Bord3")
  35         M.setGroupAtPlan(yinf,1,eps,"Bord4")
  36         M.setGroupAtPlan(zsup,2,eps,"Bord5")
  37         M.setGroupAtPlan(zinf,2,eps,"Bord6")
  38
  39         print "Built a regular 3D tetrahedral mesh from ", nx,"x" ,ny,"x" ,nz, " cubic cells"
  40
  41         FEComputation=False
  42         myProblem = cf.StationaryDiffusionEquation(spaceDim,FEComputation);
  43         myProblem.setMesh(M);
  44
  45         # set the limit value for each boundary
  46         T1=0;
  47         T2=0;
  48         T3=0;
  49         T4=0;
  50         T5=0;
  51         T6=0;
  52
  53         myProblem.setDirichletBoundaryCondition("Bord1",T1)
  54         myProblem.setDirichletBoundaryCondition("Bord2",T2)
  55         myProblem.setDirichletBoundaryCondition("Bord3",T3)
  56         myProblem.setDirichletBoundaryCondition("Bord4",T4)
  57         myProblem.setDirichletBoundaryCondition("Bord5",T5)
  58         myProblem.setDirichletBoundaryCondition("Bord6",T6)
  59
  60         #Set the right hand side function
  61         my_RHSfield = cm.Field("RHS_field", cm.CELLS, M, 1)
  62         for i in range(M.getNumberOfCells()):
  63                 Ci= M.getCell(i)
  64                 x = Ci.x()
  65                 y = Ci.y()
  66                 z = Ci.z()
  67
  68                 my_RHSfield[i]=2*pi*pi*sin(pi*x)*sin(pi*y)*sin(pi*z)#mettre la fonction definie au second membre de l'edp
  69
  70         myProblem.setHeatPowerField(my_RHSfield)
  71         myProblem.setLinearSolver(cf.GMRES,cf.ILU);
  72
  73         # name of result file
  74         fileName = "StationaryDiffusion_3DFV_StructuredTetrahedra";
  75
  76         # computation parameters
  77         myProblem.setFileName(fileName);
  78
  79         # Run the computation
  80         myProblem.initialize();
  81         print("Running python "+ fileName );
  82
  83         ok = myProblem.solveStationaryProblem();
  84         if (not ok):
  85                 print( "Python simulation of " + fileName + "  failed ! " );
  86                 pass
  87         else:
  88                 print( "Python simulation of " + fileName + " is successful !" );
  89                 ####################### Postprocessing #########################
  90                 my_ResultField = myProblem.getOutputTemperatureField()
  91                 #The following formulas use the fact that the exact solution is equal the right hand side divided by 3*pi*pi
  92                 max_abs_sol_exacte=max(my_RHSfield.max(),-my_RHSfield.min())/(3*pi*pi)
  93                 max_sol_num=my_ResultField.max()
  94                 min_sol_num=my_ResultField.min()
  95                 erreur_abs=0
  96                 for i in range(M.getNumberOfCells()) :
  97                         if  erreur_abs < abs(my_RHSfield[i]/(3*pi*pi) - my_ResultField[i]) :
  98                                 erreur_abs = abs(my_RHSfield[i]/(3*pi*pi) - my_ResultField[i])
  99
 100                 print("Absolute error = max(| exact solution - numerical solution |) = ",erreur_abs )
 101                 print("Relative error = max(| exact solution - numerical solution |)/max(| exact solution |) = ",erreur_abs/max_abs_sol_exacte)
 102                 print ("Maximum numerical solution = ", max_sol_num, " Minimum numerical solution = ", min_sol_num)
 103
 104                 assert erreur_abs/max_abs_sol_exacte <1.
 105         pass
 106
 107         print( "------------ !!! End of calculation !!! -----------" );
 108
 109         myProblem.terminate();
 110         return ok
 111
 112 if __name__ == """__main__""":
 113         StationaryDiffusionEquation_3DFV_StructuredTetrahedra()