CoreFlows/examples/Python/TransportEquation/TransportEquation_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 transport equation
   9 # Description : dh/dt + \div h \vec v  = \Phi + \lambda_{sf} (T_{solid} - T_{fluid})
  10 #               Neumann or Dirichlet boundary conditions
  11 #               Finite volumes or finite elements
  12 # Author      : Michaël Ndjinga
  13 # Copyright   : CEA Saclay 2021
  14 #================================================================================================================================
  15
  16
  17 def TransportEquation_2DSpherical(  fileName):
  18
  19         """ Description : Test solving the transport of the enthalpy h in a fluid (water or steam).
  20                 Equation : Linear transport equation  dh/dt + \div h \vec v  = \lambda_{sf} (T_{solid} - T_{fluid})
  21                                         Phase (water or steam), pressure of the fluid MUST be defined
  22                         The fluid may be heated by a solid with temperature T_{solid} and transfer coefficient \lambda_{sf} using functions setRodTemperature and setHeatTransfertCoeff
  23                         The solid may receive some extra heat power due to nuclear fissions using function setHeatSource
  24         """
  25         #Space dimension of the problem
  26         spaceDim=2
  27
  28     # Mandatory physical values
  29         transport_velocity=[5,5]# Constant velocity vector of dimension spaceDim
  30         fluid_phase=solverlab.Water# ou bien solverlab.Air
  31         pressure_regime=solverlab.around155bars600K# ou bien solverlab.around1bar300KTransport
  32         myProblem = solverlab.TransportEquation(fluid_phase,pressure_regime,transport_velocity);
  33
  34     # Set the mesh and initial data
  35         initial_data_inputfile="../resources/meshSquare";
  36         initial_data_fieldName="Fluid temperature";
  37         print("Loading unstructured mesh and initial data", " in file ", initial_data_inputfile )
  38         initial_data_time_iteration=0# default value is 0
  39         initial_data_time_sub_iteration=0# default value is 0
  40         initial_data_time_meshLevel=0# default value is 0
  41         myProblem.setInitialField(initial_data_inputfile, initial_data_fieldName, initial_data_time_iteration, initial_data_time_sub_iteration, initial_data_time_meshLevel, solverlab.CELLS)
  42
  43     #### Optional physical values (default value is zero) : fluid temperature field, heat transfert coefficient, heat power field
  44         # Loading and setting fluid temperature field
  45         solid_temperature_inputfile="../resources/meshSquare";
  46         solid_temperature_fieldName="Fluid temperature";
  47         solid_temperature_time_iteration=0# default value is 0
  48         solid_temperature_time_sub_iteration=0# default value is 0
  49         solid_temperature_meshLevel=0# default value is 0
  50         print("Loading field :", solid_temperature_fieldName, " in file ", solid_temperature_inputfile)
  51         solidTemperatureField=solverlab.Field(solid_temperature_inputfile, solverlab.CELLS, solid_temperature_fieldName, solid_temperature_time_iteration, solid_temperature_time_sub_iteration, solid_temperature_meshLevel)
  52         myProblem.setRodTemperatureField(solidTemperatureField)
  53         # Setting heat transfert coefficient
  54         heatTransfertCoeff=1000.;#fluid/solid exchange coefficient, default value is 0
  55         myProblem.setHeatTransfertCoeff(heatTransfertCoeff);
  56         # Loading heat power field
  57         heat_power_inputfile="../resources/meshSquare";
  58         heat_power_fieldName="Heat power";
  59         heat_power_time_iteration=0# default value is 0
  60         heat_power_time_sub_iteration=0# default value is 0
  61         heat_power_meshLevel=0# default value is 0
  62         print("Loading field :", heat_power_fieldName, " in file ", heat_power_inputfile)
  63         heatPowerField=solverlab.Field(heat_power_inputfile, solverlab.CELLS, heat_power_fieldName, heat_power_time_iteration, heat_power_time_sub_iteration, heat_power_meshLevel)
  64         myProblem.setHeatPowerField(heatPowerField)
  65
  66     # the boundary conditions :
  67         boundaryGroupNames=myProblem.getMesh().getNameOfFaceGroups()
  68         print(len(boundaryGroupNames), " Boundary Face Group detected : ", boundaryGroupNames)
  69
  70         # for each boundary we load the boundary field (replace by a loop over the boundaries)
  71         boundary1_type=solverlab.NeumannTransport
  72         boundary1_inputfile="../resources/meshSquare";
  73         boundary1_fieldName="Left temperature";
  74         boundary1_time_iteration=0# default value is 0
  75         boundary1_time_sub_iteration=0# default value is 0
  76         boundary1_meshLevel=0# default value is 0
  77         print("Boundary ", boundaryGroupNames[3], ", loading field :", boundary1_fieldName, " in file ", boundary1_inputfile)
  78         boundary1Field=solverlab.Field(boundary1_inputfile, solverlab.CELLS, boundary1_fieldName, boundary1_time_iteration, boundary1_time_sub_iteration, boundary1_meshLevel)
  79         boundary2_type=solverlab.DirichletTransport
  80         boundary2_inputfile="../resources/meshSquare";
  81         boundary2_fieldName="Right temperature";
  82         boundary2_time_iteration=0# default value is 0
  83         boundary2_time_sub_iteration=0# default value is 0
  84         boundary2_meshLevel=0# default value is 0
  85         print("Boundary ", boundaryGroupNames[2], ", loading field :", boundary2_fieldName, " in file ", boundary2_inputfile)
  86         boundary2Field=solverlab.Field(boundary2_inputfile, solverlab.CELLS, boundary2_fieldName, boundary2_time_iteration, boundary2_time_sub_iteration, boundary2_meshLevel)
  87         boundary3_type=solverlab.NeumannTransport
  88         boundary3_inputfile="../resources/meshSquare";
  89         boundary3_fieldName="Top temperature";
  90         boundary3_time_iteration=0# default value is 0
  91         boundary3_time_sub_iteration=0# default value is 0
  92         boundary3_meshLevel=0# default value is 0
  93         print("Boundary ", boundaryGroupNames[4], ", loading field :", boundary3_fieldName, " in file ", boundary3_inputfile)
  94         boundary3Field=solverlab.Field(boundary3_inputfile, solverlab.CELLS, boundary3_fieldName, boundary3_time_iteration, boundary3_time_sub_iteration, boundary3_meshLevel)
  95         boundary4_type=solverlab.DirichletTransport
  96         boundary4_inputfile="../resources/meshSquare";
  97         boundary4_fieldName="Bottom temperature";
  98         boundary4_time_iteration=0# default value is 0
  99         boundary4_time_sub_iteration=0# default value is 0
 100         boundary4_meshLevel=0# default value is 0
 101         print("Boundary ", boundaryGroupNames[1], ", loading field :", boundary4_fieldName, " in file ", boundary4_inputfile)
 102         boundary4Field=solverlab.Field(boundary4_inputfile, solverlab.CELLS, boundary4_fieldName, boundary4_time_iteration, boundary4_time_sub_iteration, boundary4_meshLevel)
 103
 104         # for each boundary we need to know if we want a Neumann or a Dirichlet boundary condition
 105         if boundary1_type==solverlab.NeumannTransport :
 106                 myProblem.setNeumannBoundaryCondition("Left", boundary1Field)
 107         elif boundary1_type==solverlab.DirichletTransport :
 108                 myProblem.setDirichletBoundaryCondition("Left", boundary1Field)
 109         if boundary2_type==solverlab.NeumannTransport :
 110                 myProblem.setNeumannBoundaryCondition("Right", boundary2Field)
 111         elif boundary2_type==solverlab.DirichletTransport :
 112                 myProblem.setDirichletBoundaryCondition("Right", boundary2Field)
 113         if boundary3_type==solverlab.NeumannTransport :
 114                 myProblem.setNeumannBoundaryCondition("Top", boundary3Field)
 115         elif boundary3_type==solverlab.DirichletTransport :
 116                 myProblem.setDirichletBoundaryCondition("Top", boundary3Field);
 117         if boundary4_type==solverlab.NeumannTransport :
 118                 myProblem.setNeumannBoundaryCondition("Bottom", boundary4Field)
 119         elif boundary4_type==solverlab.DirichletTransport :
 120                 myProblem.setDirichletBoundaryCondition("Bottom", boundary4Field);
 121
 122     # set the numerical method
 123         myProblem.setTimeScheme( solverlab.Explicit)# Otherwise solverlab.Implicit
 124         max_nb_its_lin_solver = 50
 125         myProblem.setLinearSolver(solverlab.GMRES, solverlab.ILU, max_nb_its_lin_solver );
 126
 127     # computation parameters
 128         MaxNbOfTimeStep = 3 ;# default value is 10
 129         freqSave = 1;# default value is 1
 130         cfl = 0.95;# default value is 1
 131         maxTime = 100000000;# default value is 10
 132         precision = 1e-6;# default value is 1e-6
 133         result_directory="."# default value = current directory
 134
 135         myProblem.setCFL(cfl);
 136         myProblem.setPrecision(precision);
 137         myProblem.setMaxNbOfTimeStep(MaxNbOfTimeStep);
 138         myProblem.setTimeMax(maxTime);
 139         myProblem.setFreqSave(freqSave);
 140         myProblem.setFileName(fileName);
 141         myProblem.setResultDirectory(result_directory)
 142         myProblem.setSaveFileFormat(solverlab.MED)#default value is solverlab.VTK
 143
 144     # evolution
 145         myProblem.initialize();
 146         print("Running python test "+ fileName );
 147
 148         ok = myProblem.run();
 149         if (ok):
 150                 print( "Python simulation " + fileName + " is successful !" );
 151                 pass
 152         else:
 153                 print( "Python simulation " + fileName + "  failed ! " );
 154                 pass
 155
 156         print( "------------ End of simulation !!! -----------" );
 157
 158         myProblem.terminate();
 159         return ok
 160
 161 if __name__ == """__main__""":
 162         # name of result file
 163         fileName = "2DSpherical";# default value is ""
 164         TransportEquation_2DSpherical( fileName)