CoreFlows/examples/Python/DiffusionEquation/DiffusionEquation_2DSpherical.py

   1 #!/usr/bin/env python3
   2 # -*-coding:utf-8 -*
   3
   4 import sys
   5 import solverlab
   6
   7 #===============================================================================================================================
   8 # Name        : Simulation of a 2D heat equation
   9 # Description : Test solving the diffusion of the temperature T in a solid
  10 #               \rho cp dT/dt-\lambda\Delta T=\Phi + \lambda_{sf} (T_{fluid}-T_{solid})
  11 #               Neumann or Dirichlet boundary conditions
  12 #               Finite elements or finite volumes
  13 # Author      : Michaël Ndjinga
  14 # Copyright   : CEA Saclay 2021
  15 #================================================================================================================================
  16
  17
  18 def DiffusionEquation_2DSpherical(FECalculation, fileName):
  19
  20         """ Description : Test solving the diffusion of the temperature T in a solid (default is Uranium).
  21                 Equation : Thermal diffusion equation  \rho cp dT/dt-\lambda\Delta T=\Phi + \lambda_{sf} (T_{fluid}-T_{solid})
  22                         Heat capacity cp, density \rho, and conductivity \lambda of the solid MUST be defined
  23                         The solid may be refrigerated by a fluid with temperature T_{solid} transfer coefficient \lambda_{sf} using functions setFluidTemperature and setHeatTransfertCoeff
  24                         The solid may receive some extra heat power due to nuclear fissions using function setHeatSource
  25         """
  26         #Space dimension of the problem
  27         spaceDim=2
  28
  29     # Mandatory physical values
  30         solid_specific_heat=300# specific heat capacity, default value 300
  31         solid_density=10000# density, default value 10000
  32         solid_conductivity=5# conductivity, default value 5
  33
  34         myProblem = solverlab.DiffusionEquation(spaceDim,FECalculation,solid_density,solid_specific_heat,solid_conductivity);
  35
  36         # Definition of field support parameter
  37         if( FECalculation):
  38                 supportOfField=solverlab.NODES
  39         else:
  40                 supportOfField=solverlab.CELLS
  41
  42     # Set the mesh and initial data
  43         initial_data_inputfile="../resources/meshSquare";
  44         initial_data_fieldName="Solid temperature";
  45         print("Loading unstructured mesh and initial data", " in file ", initial_data_inputfile )
  46         initial_data_time_iteration=0# default value is 0
  47         initial_data_time_sub_iteration=0# default value is 0
  48         initial_data_time_meshLevel=0# default value is 0
  49         myProblem.setInitialField(initial_data_inputfile, initial_data_fieldName, initial_data_time_iteration, initial_data_time_sub_iteration, initial_data_time_meshLevel, supportOfField)
  50
  51     #### Optional physical values (default value is zero) : fluid temperature field, heat transfert coefficient, heat power field
  52         # Loading and setting fluid temperature field
  53         fluid_temperature_inputfile="../resources/meshSquare";
  54         fluid_temperature_fieldName="Fluid temperature";
  55         fluid_temperature_time_iteration=0# default value is 0
  56         fluid_temperature_time_sub_iteration=0# default value is 0
  57         fluid_temperature_meshLevel=0# default value is 0
  58         print("Loading field :", fluid_temperature_fieldName, " in file ", fluid_temperature_inputfile)
  59         fluidTemperatureField=solverlab.Field(fluid_temperature_inputfile, supportOfField, fluid_temperature_fieldName, fluid_temperature_time_iteration, fluid_temperature_time_sub_iteration, fluid_temperature_meshLevel)
  60         myProblem.setFluidTemperatureField(fluidTemperatureField)
  61         # Setting heat transfert coefficient
  62         heatTransfertCoeff=1000.;#fluid/solid exchange coefficient, default value is 0
  63         myProblem.setHeatTransfertCoeff(heatTransfertCoeff);
  64         # Loading heat power field
  65         heat_power_inputfile="../resources/meshSquare";
  66         heat_power_fieldName="Heat power";
  67         heat_power_time_iteration=0# default value is 0
  68         heat_power_time_sub_iteration=0# default value is 0
  69         heat_power_meshLevel=0# default value is 0
  70         print("Loading field :", heat_power_fieldName, " in file ", heat_power_inputfile)
  71         heatPowerField=solverlab.Field(heat_power_inputfile, supportOfField, heat_power_fieldName, heat_power_time_iteration, heat_power_time_sub_iteration, heat_power_meshLevel)
  72         myProblem.setHeatPowerField(heatPowerField)
  73
  74     # the boundary conditions :
  75         if( FECalculation):
  76                 boundaryGroupNames=myProblem.getMesh().getNameOfNodeGroups()
  77                 print(len(boundaryGroupNames), " Boundary Node Group detected : ", boundaryGroupNames)
  78         else:
  79                 boundaryGroupNames=myProblem.getMesh().getNameOfFaceGroups()
  80                 print(len(boundaryGroupNames), " Boundary Face Group detected : ", boundaryGroupNames)
  81
  82         # for each boundary we load the boundary field (replace by a loop over the boundaries)
  83         boundary1_type=solverlab.NeumannDiffusion
  84         boundary1_inputfile="../resources/meshSquare";
  85         boundary1_fieldName="Left temperature";
  86         boundary1_time_iteration=0# default value is 0
  87         boundary1_time_sub_iteration=0# default value is 0
  88         boundary1_meshLevel=0# default value is 0
  89         print("Boundary ", boundaryGroupNames[3], ", loading field :", boundary1_fieldName, " in file ", boundary1_inputfile)
  90         boundary1Field=solverlab.Field(boundary1_inputfile, supportOfField, boundary1_fieldName, boundary1_time_iteration, boundary1_time_sub_iteration, boundary1_meshLevel)
  91         boundary2_type=solverlab.DirichletDiffusion
  92         boundary2_inputfile="../resources/meshSquare";
  93         boundary2_fieldName="Right temperature";
  94         boundary2_time_iteration=0# default value is 0
  95         boundary2_time_sub_iteration=0# default value is 0
  96         boundary2_meshLevel=0# default value is 0
  97         print("Boundary ", boundaryGroupNames[2], ", loading field :", boundary2_fieldName, " in file ", boundary2_inputfile)
  98         boundary2Field=solverlab.Field(boundary2_inputfile, supportOfField, boundary2_fieldName, boundary2_time_iteration, boundary2_time_sub_iteration, boundary2_meshLevel)
  99         boundary3_type=solverlab.NeumannDiffusion
 100         boundary3_inputfile="../resources/meshSquare";
 101         boundary3_fieldName="Top temperature";
 102         boundary3_time_iteration=0# default value is 0
 103         boundary3_time_sub_iteration=0# default value is 0
 104         boundary3_meshLevel=0# default value is 0
 105         print("Boundary ", boundaryGroupNames[4], ", loading field :", boundary3_fieldName, " in file ", boundary3_inputfile)
 106         boundary3Field=solverlab.Field(boundary3_inputfile, supportOfField, boundary3_fieldName, boundary3_time_iteration, boundary3_time_sub_iteration, boundary3_meshLevel)
 107         boundary4_type=solverlab.DirichletDiffusion
 108         boundary4_inputfile="../resources/meshSquare";
 109         boundary4_fieldName="Bottom temperature";
 110         boundary4_time_iteration=0# default value is 0
 111         boundary4_time_sub_iteration=0# default value is 0
 112         boundary4_meshLevel=0# default value is 0
 113         print("Boundary ", boundaryGroupNames[1], ", loading field :", boundary4_fieldName, " in file ", boundary4_inputfile)
 114         boundary4Field=solverlab.Field(boundary4_inputfile, supportOfField, boundary4_fieldName, boundary4_time_iteration, boundary4_time_sub_iteration, boundary4_meshLevel)
 115
 116         # for each boundary we need to know if we want a Neumann or a Dirichlet boundary condition
 117         if boundary1_type==solverlab.NeumannDiffusion :
 118                 myProblem.setNeumannBoundaryCondition("Left", boundary1Field)
 119         elif boundary1_type==solverlab.DirichletDiffusion :
 120                 myProblem.setDirichletBoundaryCondition("Left", boundary1Field)
 121         if boundary2_type==solverlab.NeumannDiffusion :
 122                 myProblem.setNeumannBoundaryCondition("Right", boundary2Field)
 123         elif boundary2_type==solverlab.DirichletDiffusion :
 124                 myProblem.setDirichletBoundaryCondition("Right", boundary2Field)
 125         if boundary3_type==solverlab.NeumannDiffusion :
 126                 myProblem.setNeumannBoundaryCondition("Top", boundary3Field)
 127         elif boundary3_type==solverlab.DirichletDiffusion :
 128                 myProblem.setDirichletBoundaryCondition("Top", boundary3Field);
 129         if boundary4_type==solverlab.NeumannDiffusion :
 130                 myProblem.setNeumannBoundaryCondition("Bottom", boundary4Field)
 131         elif boundary4_type==solverlab.DirichletDiffusion :
 132                 myProblem.setDirichletBoundaryCondition("Bottom", boundary4Field);
 133
 134     # set the numerical method
 135         myProblem.setTimeScheme( solverlab.Explicit)# Otherwise solverlab.Implicit
 136         max_nb_its_lin_solver = 50
 137         myProblem.setLinearSolver(solverlab.GMRES, solverlab.ILU, max_nb_its_lin_solver );
 138
 139     # computation parameters
 140         MaxNbOfTimeStep = 3 ;# default value is 10
 141         freqSave = 1;# default value is 1
 142         cfl = 0.95;# default value is 1
 143         maxTime = 100000000;# default value is 10
 144         precision = 1e-6;# default value is 1e-6
 145         result_directory="."# default value = current directory
 146
 147         myProblem.setCFL(cfl);
 148         myProblem.setPrecision(precision);
 149         myProblem.setMaxNbOfTimeStep(MaxNbOfTimeStep);
 150         myProblem.setTimeMax(maxTime);
 151         myProblem.setFreqSave(freqSave);
 152         myProblem.setFileName(fileName);
 153         myProblem.setResultDirectory(result_directory)
 154         myProblem.setSaveFileFormat(solverlab.MED)#default value is solverlab.VTK
 155
 156     # evolution
 157         myProblem.initialize();
 158         print("Running python test "+ fileName );
 159
 160         ok = myProblem.run();
 161         if (ok):
 162                 print( "Python simulation " + fileName + " is successful !" );
 163                 pass
 164         else:
 165                 print( "Python simulation " + fileName + "  failed ! " );
 166                 pass
 167
 168         print( "------------ End of simulation !!! -----------" );
 169
 170         myProblem.terminate();
 171         return ok
 172
 173 if __name__ == """__main__""":
 174     if len(sys.argv) >1 :
 175         FECalculation = bool(int(sys.argv[1]))
 176         # name of result file
 177         if( FECalculation):
 178                 fileName = "2DSpherical_FE";# default value is ""
 179         else:
 180                 fileName = "2DSpherical_FV";# default value is ""
 181         DiffusionEquation_2DSpherical(FECalculation, fileName)
 182     else :
 183         raise ValueError("DiffusionEquation_2DSpherical : missing one argument")