CoCalc -- EnthalpySolver.F90

GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/elmerice/Solvers/EnthalpySolver.F90
³²⁰⁴ views
1
!/*****************************************************************************/
2
! *  Elmer/Ice, a glaciological add-on to Elmer
3
! *  http://elmerice.elmerfem.org
4
! * 
5
! *  This program is free software; you can redistribute it and/or
6
! *  modify it under the terms of the GNU General Public License
7
! *  as published by the Free Software Foundation; either version 2
8
! *  of the License, or (at your option) any later version.
9
! * 
10
! *  This program is distributed in the hope that it will be useful,
11
! *  but WITHOUT ANY WARRANTY; without even the implied warranty of
12
! *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13
! *  GNU General Public License for more details.
14
! *
15
! *  You should have received a copy of the GNU General Public License
16
! *  along with this program (in file fem/GPL-2); if not, write to the 
17
! *  Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 
18
! *  Boston, MA 02110-1301, USA.
19
! *
20
! *****************************************************************************/
21
!
22
!/******************************************************************************
23
! *
24
! *  Module containing a solver for enthalpy equation in glaciology
25
! *
26
! ******************************************************************************
27
! *
28
! *  Authors: Adrien Gilbert, 
29
! *
30
! *  Original Date: Jun 2014
31
! *
32
! *****************************************************************************/
33

34
!------------------------------------------------------------------------------
35
   SUBROUTINE EnthalpySolver( Model,Solver,Timestep,TransientSimulation )
36
!------------------------------------------------------------------------------
37
     USE DiffuseConvective
38
     USE DiffuseConvectiveGeneral
39
     USE Differentials
40
     USE Radiation
41
     USE MaterialModels
42
     USE Adaptive
43
     USE DefUtils
44

45
!------------------------------------------------------------------------------
46
     IMPLICIT NONE
47
!------------------------------------------------------------------------------
48
     INTEGER, PARAMETER :: PHASE_SPATIAL_1 = 1
49
     INTEGER, PARAMETER :: PHASE_SPATIAL_2 = 2
50
     INTEGER, PARAMETER :: PHASE_TEMPORAL  = 3
51
    
52
     TYPE(Model_t)  :: Model
53
     TYPE(Solver_t), TARGET :: Solver
54

55
     LOGICAL :: TransientSimulation
56
     REAL(KIND=dp) :: Timestep
57
!------------------------------------------------------------------------------
58
!    Local variables
59
!------------------------------------------------------------------------------
60
     TYPE(Matrix_t), POINTER :: StiffMatrix
61

62
     INTEGER :: i,j,k,l,m,n,nd,t,tt,iter,k1,k2,body_id,eq_id,istat,LocalNodes,bf_id
63

64
     TYPE(Nodes_t)   :: ElementNodes
65
     TYPE(Element_t),POINTER :: Element,Parent,RadiationElement
66

67
     REAL(KIND=dp) :: RelativeChange, &
68
           Norm,PrevNorm,Text,S,C,C1,Emissivity,StefanBoltzmann, &
69
           ReferencePressure=0.0d0, SpecificHeatRatio
70

71
     CHARACTER(LEN=MAX_NAME_LEN) :: RadiationFlag,ConvectionFlag
72

73
     INTEGER :: PhaseChangeModel
74
     CHARACTER(LEN=MAX_NAME_LEN) :: PhaseModel, StabilizeFlag, VarName
75

76
     INTEGER, POINTER :: NodeIndexes(:)
77
     LOGICAL :: Stabilize = .TRUE., Bubbles = .TRUE., UseBubbles,NewtonLinearization = .FALSE., &
78
         Found, GotIt, HeatFluxBC, HeatGapBC, GotMeltPoint, IsRadiation, InfBC, UnFoundFatal=.TRUE.
79
! Which compressibility model is used
80
     CHARACTER(LEN=MAX_NAME_LEN) :: CompressibilityFlag, ConvectionField
81
     INTEGER :: CompressibilityModel
82

83
     LOGICAL :: AllocationsDone = .FALSE.,PhaseSpatial=.FALSE., &
84
        PhaseChange=.FALSE., CheckLatentHeatRelease=.FALSE., FirstTime, &
85
        SmartHeaterControl, IntegralHeaterControl, HeaterControlLocal, SmartTolReached=.FALSE., &
86
        TransientHeaterControl, SmartHeaterAverage, ConstantBulk, SaveBulk, &
87
	TransientAssembly
88
     LOGICAL, POINTER :: SmartHeaters(:), IntegralHeaters(:)
89

90
     TYPE(Variable_t), POINTER :: TempSol,FlowSol,HeatSol,CurrentSol, MeshSol, DensitySol
91
     TYPE(ValueList_t), POINTER :: Equation,Material,SolverParams,BodyForce,BC,Constants
92

93
     INTEGER, POINTER :: EnthalpyPerm(:),FlowPerm(:),CurrentPerm(:),MeshPerm(:)
94

95
     INTEGER :: NSDOFs,NewtonIter,NonlinearIter,MDOFs, &
96
         SmartHeaterBC, SmartHeaterNode, DoneTime=0
97
     REAL(KIND=dp) :: NonlinearTol,NewtonTol,SmartTol,Relax, &
98
            SaveRelax,dt,dt0,CumulativeTime, VisibleFraction, PowerScaling=1.0, PrevPowerScaling=1.0, &
99
            PowerRelax, PowerTimeScale, PowerSensitivity, xave, yave, Normal(3), &
100
	    dist, mindist, ControlPoint(3)
101

102
     REAL(KIND=dp), POINTER :: Enthalpy_h(:),PrevEnthalpy_h(:),FlowSolution(:), &
103
       ElectricCurrent(:), PhaseChangeIntervals(:,:),ForceVector(:), &
104
       PrevSolution(:), HC(:), Hwrk(:,:,:),MeshVelocity(:), XX(:), YY(:),ForceHeater(:),&
105
       RealWork(:,:)
106

107
     REAL(KIND=dp), ALLOCATABLE :: vals(:)
108
     REAL(KIND=dp) :: Jx,Jy,Jz,JAbs, Power, MeltPoint, IntHeatSource
109

110
     REAL(KIND=dp) :: Tref,Tm,P,hm,hi,L_heat,A_cap,B_cap,Ptriple,Psurf,beta
111

112
     CHARACTER(LEN=MAX_NAME_LEN) :: PressureName,WaterName,PhaseEnthName,TempName
113

114
     TYPE(Variable_t), POINTER :: PressureVariable,PhaseChangeEnthVar,WaterVar,TemphomoVar
115
     REAL(KIND=dp), POINTER :: PressureValues(:),PhaseChangeEnthValues(:)
116
     INTEGER, POINTER :: PressurePerm(:),PhaseChangeEnthPerm(:)
117

118
      INTEGER, ALLOCATABLE, SAVE :: Indexes(:), SaveIndexes(:)
119

120
     REAL(KIND=dp), ALLOCATABLE :: MASS(:,:), &
121
       STIFF(:,:), LOAD(:), HeatConductivity(:,:,:), &
122
       FORCE(:), U(:), V(:), W(:), MU(:,:),TimeForce(:), &
123
       Density(:), LatentHeat(:), HeatTransferCoeff(:), &
124
       HeatCapacity(:),WaterDiffusivity(:), Enthalpy(:), Viscosity(:), LocalEnthalpy_h(:), &
125
       NodalEmissivity(:), ElectricConductivity(:), Permeability(:), Work(:), C0(:), &
126
       Pressure(:), dPressuredt(:), GasConstant(:),AText(:), HeaterArea(:), &
127
       HeaterTarget(:), HeaterScaling(:), HeaterDensity(:), HeaterSource(:), &
128
       HeatExpansionCoeff(:), ReferenceEnthalpy_h(:), PressureCoeff(:), &
129
       PhaseVelocity(:,:), HeatConductivityIso(:), &
130
       PerfusionRate(:), PerfusionDensity(:), PerfusionHeatCapacity(:), PerfusionRefEnthalpy_h(:)
131

132
     SAVE U, V, W, MU, MASS, STIFF, LOAD, PressureCoeff, &
133
       FORCE, ElementNodes, HeatConductivity, HeatCapacity,WaterDiffusivity, HeatTransferCoeff, &
134
       Enthalpy, Density, LatentHeat, PhaseVelocity, AllocationsDone, Viscosity, TimeForce, &
135
       LocalNodes, LocalEnthalpy_h, Work, ElectricConductivity, &
136
       NodalEmissivity, Permeability, C0, dPressuredt, Pressure, &
137
       GasConstant,AText,Hwrk, XX, YY, ForceHeater, Power, HeaterArea, HeaterTarget, &
138
       HeaterScaling, HeaterDensity, HeaterSource, SmartHeaters, IntegralHeaters, SmartTolReached,    &
139
       ReferenceEnthalpy_h, HeatExpansionCoeff, PrevPowerScaling, PowerScaling, &
140
       MeltPoint, DoneTime, SmartHeaterNode, SmartHeaterBC, SmartHeaterAverage, &
141
       HeatConductivityIso, &
142
       PerfusionRate, PerfusionDensity, PerfusionHeatCapacity, PerfusionRefEnthalpy_h
143

144

145
     INTERFACE
146
        SUBROUTINE EnthalpySolver_Boundary_Residual( Model,Edge,Mesh,Quant,Perm,Gnorm,Indicator)
147
          USE Types
148
          TYPE(Element_t), POINTER :: Edge
149
          TYPE(Model_t) :: Model
150
          TYPE(Mesh_t), POINTER :: Mesh
151
          REAL(KIND=dp) :: Quant(:), Indicator(2), Gnorm
152
          INTEGER :: Perm(:)
153
        END SUBROUTINE EnthalpySolver_Boundary_Residual
154

155
        SUBROUTINE EnthalpySolver_Edge_Residual( Model,Edge,Mesh,Quant,Perm,Indicator)
156
          USE Types
157
          TYPE(Element_t), POINTER :: Edge
158
          TYPE(Model_t) :: Model
159
          TYPE(Mesh_t), POINTER :: Mesh
160
          REAL(KIND=dp) :: Quant(:), Indicator(2)
161
          INTEGER :: Perm(:)
162
        END SUBROUTINE EnthalpySolver_Edge_Residual
163

164
        SUBROUTINE EnthalpySolver_Inside_Residual( Model,Element,Mesh,Quant,Perm, Fnorm, Indicator)
165
          USE Types
166
          TYPE(Element_t), POINTER :: Element
167
          TYPE(Model_t) :: Model
168
          TYPE(Mesh_t), POINTER :: Mesh
169
          REAL(KIND=dp) :: Quant(:), Indicator(2), Fnorm
170
          INTEGER :: Perm(:)
171
        END SUBROUTINE EnthalpySolver_Inside_Residual
172
     END INTERFACE
173

174
  REAL(KIND=dp) :: at,at0,totat,st,totst,t1
175

176
!------------------------------------------------------------------------------
177
!    The View and Gebhart factors may change. If this is necessary, this is 
178
!    done within this subroutine. The routine is called in the
179
!    start as it may affect the matrix topology.
180
!    Newton lineariarization option is needed only when there is radiation.
181
!------------------------------------------------------------------------------
182
     IsRadiation = ListCheckPresentAnyBC( Model,'Radiation')
183
     IF( IsRadiation ) THEN
184
       CALL RadiationFactors( Solver, .FALSE.)
185
     END IF
186

187
!------------------------------------------------------------------------------
188
!    Get variables needed for solution
189
!------------------------------------------------------------------------------
190

191
     IF ( .NOT. ASSOCIATED( Solver % Matrix ) ) RETURN
192

193
     StiffMatrix => GetMatrix()
194
     ForceVector => Solver % Matrix % RHS
195

196
     TempSol => Solver % Variable
197
     EnthalpyPerm    => TempSol % Perm
198
     Enthalpy_h => TempSol % Values
199
     VarName = GetVarName( TempSol ) 
200
  
201
     LocalNodes = COUNT( EnthalpyPerm > 0 )
202
     IF ( LocalNodes <= 0 ) RETURN
203

204
     SolverParams => GetSolverParams()
205
     ConvectionField = GetString( SolverParams, 'Enthalpy_h Convection Field', Found )
206

207
     IF ( Found ) THEN
208
       FlowSol => VariableGet( Solver % Mesh % Variables, ConvectionField )
209
     ELSE
210
       FlowSol => VariableGet( Solver % Mesh % Variables, 'Flow Solution' )
211
     END IF
212

213
     IF ( ASSOCIATED( FlowSol ) ) THEN
214
       FlowPerm     => FlowSol % Perm
215
       NSDOFs       =  FlowSol % DOFs
216
       FlowSolution => FlowSol % Values
217
     END IF
218

219
     DensitySol => VariableGet( Solver % Mesh % Variables, 'Enthalpy Density' )
220

221
     ! Check whether we have some heater controls. This will affect initialization stuff. 
222
     SmartHeaterControl = ListCheckPresentAnyBodyForce( Model,'Smart Heater Control')
223
     IntegralHeaterControl = ListCheckPresentAnyBodyForce( Model,'Integral Heat Source')
224
   
225
!------------------------------------------------------------------------------
226
!    Allocate some permanent storage, this is done first time only
227
!------------------------------------------------------------------------------
228
     IF ( .NOT. AllocationsDone .OR. Solver % Mesh % Changed ) THEN
229
        N = Solver % Mesh % MaxElementDOFs
230

231
        IF ( AllocationsDone ) THEN
232
          DEALLOCATE(  &
233
                 U, V, W, MU,           &
234
                 Pressure,              &
235
                 dPressureDt,           &
236
                 PressureCoeff,        &
237
                 ElementNodes % x,      &
238
                 ElementNodes % y,      &
239
                 ElementNodes % z,      &
240
                 Density,Work,          &
241
                 LatentHeat,            &
242
                 PhaseVelocity,         &
243
                 ElectricConductivity,  &
244
                 Permeability,          &
245
                 Viscosity,C0,          &
246
                 HeatTransferCoeff,     &
247
                 HeatExpansionCoeff,    &
248
                 ReferenceEnthalpy_h,  &
249
                 MASS,       &
250
                 LocalEnthalpy_h,      &
251
                 HeatCapacity,Enthalpy,WaterDiffusivity, &
252
                 NodalEmissivity,       &
253
                 GasConstant, AText,    &
254
                 HeatConductivity,      &
255
                 STIFF,LOAD,            &
256
                 Indexes, SaveIndexes,  &
257
                 FORCE, TimeForce,      &
258
                 HeatConductivityIso,   &
259
                 PerfusionRate,         &
260
                 PerfusionDensity,      &
261
                 PerfusionHeatCapacity, &
262
                 PerfusionRefEnthalpy_h )
263
       END IF
264

265
       ALLOCATE( &
266
                 Indexes(N), SaveIndexes(N),           &
267
                 U( N ),   V( N ),  W( N ),            &
268
                 MU( 3,N ),                            &
269
                 Pressure( N ),                        &
270
                 dPressureDt( N ),                     &
271
                 PressureCoeff( N ),                   &
272
                 ElementNodes % x( N ),                &
273
                 ElementNodes % y( N ),                &
274
                 ElementNodes % z( N ),                &
275
                 Density( N ),Work( N ),               &
276
                 LatentHeat( N ),                      &
277
                 PhaseVelocity(3, N),                  &
278
                 ElectricConductivity(N),              &
279
                 Permeability(N),                      &
280
                 Viscosity(N),C0(N),                   &
281
                 HeatTransferCoeff( N ),               &
282
                 HeatExpansionCoeff( N ),              &
283
                 ReferenceEnthalpy_h( N ),            &
284
                 MASS(  2*N,2*N ),                     &
285
                 LocalEnthalpy_h( N ),                &
286
                 HeatCapacity( N ),Enthalpy( N ), WaterDiffusivity( N ),     &
287
                 NodalEmissivity( N ),                 &
288
                 GasConstant( N ),AText( N ),          &
289
                 HeatConductivity( 3,3,N ),            &
290
                 HeatConductivityIso( N ),             &
291
                 STIFF( 2*N,2*N ),LOAD( N ), &
292
                 FORCE( 2*N ), TimeForce(2*N), &
293
                 PerfusionRate( N ),     &
294
                 PerfusionDensity( N ),      &
295
                 PerfusionHeatCapacity( N ), &
296
                 PerfusionRefEnthalpy_h( N ), &
297
                 STAT=istat )
298

299
       IF ( istat /= 0 ) THEN
300
         CALL Fatal( 'HeatSolve', 'Memory allocation error' )
301
       END IF
302

303

304
       IF ( SmartHeaterControl .OR. IntegralHeaterControl) THEN          
305
          n = Model % NumberOfBodyForces
306
          IF ( AllocationsDone ) DEALLOCATE( HeaterArea, HeaterDensity, HeaterSource, &
307
               HeaterScaling, HeaterTarget, SmartHeaters, IntegralHeaters )
308
          ALLOCATE( HeaterArea(n), HeaterDensity(n), HeaterSource(n), &
309
               HeaterScaling(n), HeaterTarget(n), SmartHeaters(n), &
310
               IntegralHeaters(n) )
311
          IF ( istat /= 0 ) THEN
312
             CALL Fatal( 'HeatSolve', 'Memory allocation error' )
313
          END IF
314
          SmartHeaters = .FALSE.
315
          IntegralHeaters = .FALSE.
316
       END IF
317
       
318
       IF( SmartHeaterControl ) THEN
319
          IF ( AllocationsDone ) DEALLOCATE( XX, YY, ForceHeater  )
320
          n = SIZE( Enthalpy_h )
321
          ALLOCATE( XX( n ), YY(n), ForceHeater( n ), STAT=istat )
322
          IF ( istat /= 0 ) THEN
323
             CALL Fatal( 'HeatSolve', 'Memory allocation error' )
324
          END IF
325
          XX = 0.0d0 
326
          YY = 0.0d0
327
          ForceHeater = 0.0d0
328
       END IF
329
       
330
       NULLIFY( Hwrk )
331
       AllocationsDone = .TRUE.
332
    END IF
333

334
!------------------------------------------------------------------------------
335
!    Do some additional initialization, and go for it
336
!------------------------------------------------------------------------------
337
     dt = Timestep
338
     Constants => GetConstants()
339
     IF( IsRadiation ) THEN
340
       StefanBoltzmann = ListGetConstReal( Model % Constants, &
341
                     'Stefan Boltzmann' )
342
     END IF
343

344
!------------------------------------------------------------------------------
345
     Stabilize = GetLogical( SolverParams,'Stabilize',Found )
346

347
     UseBubbles = GetLogical( SolverParams,'Bubbles',Found )
348
     IF ( .NOT.Found ) UseBubbles = .TRUE.
349

350
     StabilizeFlag = GetString( SolverParams, &
351
          'Stabilization Method',Found )
352

353
     SELECT CASE(StabilizeFlag)
354
     CASE('vms')
355
       Stabilize = .FALSE.
356
       UseBubbles= .FALSE.
357
     CASE('stabilized')
358
       Stabilize = .TRUE.
359
       UseBubbles = .FALSE.
360
     CASE('bubbles')
361
       Stabilize = .FALSE.
362
       UseBubbles = .TRUE.
363
     END SELECT
364

365
     NonlinearIter = GetInteger(   SolverParams, &
366
                     'Nonlinear System Max Iterations', Found )
367
     IF ( .NOT.Found ) NonlinearIter = 1
368

369
     NonlinearTol  = GetConstReal( SolverParams, &
370
                     'Nonlinear System Convergence Tolerance',    Found )
371

372
     IF( IsRadiation ) THEN
373
       NewtonTol     = GetConstReal( SolverParams, &
374
                      'Nonlinear System Newton After Tolerance',  Found )
375
       NewtonIter    = GetInteger(   SolverParams, &
376
                      'Nonlinear System Newton After Iterations', Found )
377
     ELSE
378
       NewtonTol = 1.0_dp
379
       NewtonIter =  0
380
     END IF
381
     IF ( NewtonIter == 0) NewtonLinearization = .TRUE.
382

383
     Relax = GetCReal( SolverParams, &
384
               'Nonlinear System Relaxation Factor',Found )
385
     IF ( .NOT.Found ) Relax = 1
386

387
     TransientAssembly = TransientSimulation
388
     dt0 = ListGetConstReal(SolverParams,'Steady State Transition Timestep',Found)
389
     IF(.NOT. Found) dt0 = ListGetConstReal(SolverParams,'Smart Heater Time Scale',Found)
390

391
     IF(Found .AND. dt > dt0) TransientAssembly = .FALSE.
392

393

394
!------------------------------------------------------------------------------
395

396
     TransientHeaterControl = .FALSE.
397
     IF(SmartHeaterControl) THEN
398

399
       ! Mark the smart heaters 
400
       SmartHeaters = .FALSE.
401
       bf_id = 0
402
       DO i = 1,Model % NumberOfBodyForces
403
         IF( ListGetLogical( Model % BodyForces(i) % Values, &
404
             'Smart Heater Control', Found ) ) THEN
405
           SmartHeaters(i) = .TRUE.	     
406
           bf_id = i
407
         END IF
408
       END DO
409

410
       ! Find the BC that controls the heater 
411
       ! If not found assume that smart heater is related to phase change 
412
       MeltPoint = GetCReal( Model % BodyForces(bf_id) % Values,&
413
           'Smart Heater Enthalpy_h',GotMeltPoint)           
414
              
415
       SmartHeaterAverage = .FALSE.
416
       SmartHeaterNode = ListGetInteger( Model % BodyForces(bf_id) % Values,&
417
           'Smart Heater Control Node',GotIt) 
418
       IF(.NOT. GotIt) THEN
419
         RealWork => ListGetConstRealArray( Model % BodyForces(bf_id) % Values,&
420
             'Smart Heater Control Point',GotIt) 
421
         IF( GotIt ) THEN
422
           ControlPoint(1:3) = RealWork(1:3,1)
423
           
424
           mindist = HUGE( mindist )
425
           DO l=1,Model % NumberOfNodes
426
             IF( EnthalpyPerm(l) == 0 ) CYCLE
427
             
428
             jx = Model % Mesh % Nodes % x(l)
429
             jy = Model % Mesh % Nodes % y(l)
430
             jz = Model % Mesh % Nodes % z(l)
431
             
432
             dist = (ControlPoint(1)-jx)**2 + (ControlPoint(2)-jy)**2 + (ControlPoint(3)-jz)**2
433
             IF( dist < mindist ) THEN
434
               mindist = dist
435
               SmartHeaterNode = l
436
             END IF
437
           END DO
438
         END IF
439

440
         WRITE(Message,*) 'Found Control Point at distance:',SQRT(mindist)
441
         CALL Info('HeatSolve',Message)
442
         WRITE(Message,*) 'Control Point Index:',SmartHeaterNode
443
         CALL Info('HeatSolve',Message)        
444
       END IF
445
       
446
       IF( .NOT. GotMeltPoint .OR. SmartHeaterNode == 0) THEN
447
         GotIt = .FALSE.
448
         Found = .FALSE.
449
         SmartHeaterBC = 0
450
         
451
         DO i=1,Model % NumberOfBCs
452
           GotIt = ListGetLogical( Model % BCs(i) % Values,'Smart Heater Boundary', Found ) 
453
           IF(GotIt) THEN
454
             SmartHeaterBC = i
455
             EXIT
456
           END IF
457
         END DO
458
         IF(.NOT. GotIt) THEN
459
           DO i=1,Model % NumberOfBCs
460
             GotIt = ListGetLogical( Model % BCs(i) % Values,'Phase Change', Found ) 
461
             IF(GotIt) THEN
462
               SmartHeaterBC = i
463
               EXIT
464
             END IF
465
           END DO
466
         END IF
467
         IF(SmartHeaterBC == 0) THEN
468
           CALL Fatal('HeatSolve','Smart Heater Boundary / Phase Change is undefined')
469
         END IF
470
         
471
         MeltPoint = GetCReal( Model % BCs(SmartHeaterBC) % Values,&
472
             'Smart Heater Enthalpy_h',Found)
473
         IF(.NOT. Found) THEN
474
           DO k=1, Model % NumberOfMaterials
475
             MeltPoint = GetCReal( Model % Materials(k) % Values, &
476
                 'Melting Point', Found )
477
             IF(Found) EXIT
478
           END DO
479
           IF(.NOT. Found) THEN
480
             CALL Fatal('HeatSolver','Smart Heater Enthalpy_h / Melting Point is undefined')
481
           END IF
482
         END IF
483
         
484
         ! Find the node related to Enthalpy_h control 
485
         SmartHeaterAverage = ListGetLogical(Solver % Values,'Smart Heater Average', Found)
486
         IF(.NOT. SmartHeaterAverage) THEN
487
           jx = -HUGE(jx)
488
           DO k = Model % Mesh % NumberOfBulkElements + 1, &
489
               Model % Mesh % NumberOfBulkElements + Model % Mesh % NumberOfBoundaryElements
490
             
491
             Element => Model % Mesh % Elements(k)
492
             
493
             IF ( Element % BoundaryInfo % Constraint == SmartHeaterBC ) THEN
494
               DO l=1,Element % TYPE % NumberOfNodes
495
                 IF ( Model % Mesh % Nodes % x(Element % NodeIndexes(l)) >= jx ) THEN
496
                   j = Element % NodeIndexes(l) 
497
                   jx = Model % Mesh % Nodes % x(Element % NodeIndexes(l))
498
                 END IF
499
               END DO
500
             END IF
501
           END DO
502
           SmartHeaterNode = j
503
         END IF
504
       END IF
505

506
        SmartTol  = GetConstReal( SolverParams, &
507
             'Smart Heater Control After Tolerance',  Found )
508
        IF(.NOT. Found) THEN
509
          SmartTolReached = .TRUE.
510
          SmartTol = 1.0
511
        END IF   
512
     
513
        PowerTimeScale = ListGetConstReal(Solver % Values, &
514
             'Smart Heater Time Scale',Found)
515

516
        IF(TransientSimulation .AND. dt < PowerTimeScale) THEN
517
           TransientHeaterControl = .TRUE.
518
           CALL Info( 'HeatSolve', 'Using Transient Heater Control')
519
        ELSE
520
           TransientHeaterControl = .FALSE.
521
           CALL Info( 'HeatSolve', 'Using Steady-state Heater Control')
522
        END IF
523
        
524
        IF(Solver % DoneTime /= DoneTime) THEN
525
           PrevPowerScaling = PowerScaling
526
           DoneTime = Solver % DoneTime
527
        END IF
528
     END IF
529
    
530
     IF( IntegralHeaterControl) THEN
531
        CALL Info( 'HeatSolve', 'Using Integral Heater Control')       
532
        IntegralHeaters = .FALSE.
533
        DO i = 1,Model % NumberOfBodyForces
534
           IntegralHeaters(i) = ListCheckPresent( Model % BodyForces(i) % Values, &
535
                'Integral Heat Source')
536
        END DO
537
     END IF
538

539
!------------------------------------------------------------------------------
540

541
     ConstantBulk = GetLogical( SolverParams, 'Constant Bulk System', Found )
542
     SaveBulk = ConstantBulk .OR. GetLogical( SolverParams, 'Save Bulk System', Found )
543
     SaveBulk = ConstantBulk .OR. GetLogical( SolverParams, 'Calculate Loads', Found )
544

545
!------------------------------------------------------------------------------
546

547
     SaveRelax = Relax
548
     CumulativeTime = 0.0d0
549

550
!------------------------------------------------------------------------------
551
     FirstTime = .TRUE.
552
     ALLOCATE(PrevSolution(LocalNodes))
553
     
554
     DO WHILE( CumulativeTime < Timestep-1.0d-12 .OR. .NOT. TransientSimulation )
555
!------------------------------------------------------------------------------
556
!    The first time around this has been done by the caller...
557
!------------------------------------------------------------------------------
558
     IF ( TransientSimulation .AND. .NOT.FirstTime ) &
559
       CALL InitializeTimestep(Solver)
560
     FirstTime = .FALSE.
561
!------------------------------------------------------------------------------
562
!    Save current solution
563
!------------------------------------------------------------------------------
564
     PrevSolution = Enthalpy_h(1:LocalNodes)
565
     IF ( TransientSimulation ) THEN
566
       PrevEnthalpy_h => Solver % Variable % PrevValues(:,1)
567
     END IF
568
!------------------------------------------------------------------------------
569

570
     totat = 0.0d0
571
     totst = 0.0d0
572

573
     Norm = Solver % Variable % Norm
574

575

576
!Calculate Phase change enthalpy ==================================================================
577
!==================================================================================================
578

579
	  A_cap = GetConstReal(Model % Constants, "Enthalpy Heat Capacity A",GotIt)
580
      IF (.NOT.GotIt) THEN
581
         CALL WARN('EnthalpySolver', 'No Keyword >Enthalpy Heat Capacity A< &
582
         & defined in model constants. Using >7.253< as default (Paterson, 1994).')
583
         A_cap = 7.253
584
      END IF
585
	  
586
	  B_cap = GetConstReal(Model % Constants, "Enthalpy Heat Capacity B",GotIt)
587
      IF (.NOT.GotIt) THEN
588
         CALL WARN('EnthalpySolver', 'No Keyword >Enthalpy Heat Capacity B< & 
589
         & defined in model constants. Using >146.3< as default (Paterson, 1994).')
590
         B_cap = 146.3
591
      END IF
592

593
  PressureName = GetString(Constants,'Pressure Variable', GotIt)
594
      IF (.NOT.GotIt) THEN
595
         CALL WARN('EnthalpySolver', 'No Keyword >Pressure Variable< defined. Using >Pressure< as default.')
596
         WRITE(PressureName,'(A)') 'Pressure'
597
      ELSE
598
         WRITE(Message,'(a,a)') 'Variable Name for pressure: ', PressureName
599
         CALL INFO('EnthalpySolve',Message,Level=12)
600
      END IF
601

602
      Tref = GetConstReal(Model % Constants, "T_ref_enthalpy",GotIt)
603
      IF (.NOT.GotIt) THEN
604
         CALL WARN('EnthalpySolver', 'No Keyword >Reference Enthalpy< defined. Using >200K< as default.')
605
         Tref = 200.0
606
      END IF
607
	  
608
	  beta = GetConstReal(Model % Constants, "beta_clapeyron",GotIt)
609
      IF (.NOT.GotIt) THEN
610
         CALL WARN('EnthalpySolver', 'No Keyword >beta_clapeyron< defined. Using >9.74 10-2 K Mpa-1< as default.')
611
         beta = 0.0974
612
      END IF
613
	  
614
	  Psurf = GetConstReal(Model % Constants, "P_surf",GotIt)
615
      IF (.NOT.GotIt) THEN
616
         CALL WARN('EnthalpySolver', 'No Keyword >Psurf< defined. Using >1.013 10-1 MPa < as default.')
617
         Psurf = 0.1013
618
      END IF
619
	  
620
	  Ptriple = GetConstReal(Model % Constants,"P_triple",GotIt)
621
      IF (.NOT.GotIt) THEN
622
         CALL WARN('EnthalpySolver', 'No Keyword >P_triple< defined. Using >0.061173 MPa < as default.')
623
         Ptriple = 0.061173
624
      END IF
625

626
      PressureVariable => VariableGet(Solver % Mesh %Variables,PressureName)
627
      IF ( ASSOCIATED( PressureVariable ) ) THEN
628
       PressurePerm    => PressureVariable % Perm
629
       PressureValues => PressureVariable % Values
630
      END IF
631

632
 PhaseEnthName = GetString(SolverParams , 'Exported Variable 1', Found )
633

634
     IF (.NOT.Found) &
635
        CALL FATAL('EnthalpySolver','No value > Exported Variable 1 (enthalpy phase change) < found in Solver')
636

637
      WaterVar => VariableGet(Solver % Mesh % Variables, WaterName)
638
      TemphomoVar => VariableGet(Solver % Mesh % Variables,TempName)
639

640
      PhaseChangeEnthVar => &
641
               VariableGet(Solver % Mesh % Variables,PhaseEnthName)
642
      IF ( ASSOCIATED( PhaseChangeEnthVar ) ) THEN
643
         PhaseChangeEnthPerm => PhaseChangeEnthVar % Perm    
644
         PhaseChangeEnthValues => PhaseChangeEnthVar % Values  
645
      END IF
646

647
do i=1,model % NumberOfNodes
648
  P=PressureValues (PressurePerm(i))+Psurf
649
  Tm=273.16-beta*(P-Ptriple)! Clapeyron
650
  PhaseChangeEnthValues (PhaseChangeEnthPerm(i)) = A_cap*0.5*Tm*Tm+B_cap*Tm-Tref*(A_cap*0.5*Tref+B_cap) ! use CP(T)=B_cap+A_cap*T
651
enddo
652

653
!=====================================================================================================================
654
!=====================================================================================================================
655

656

657
     DO iter=1,NonlinearIter
658
       at  = CPUTime()
659
       at0 = RealTime()
660

661
       CALL Info( 'HeatSolve', ' ', Level=4 )
662
       CALL Info( 'HeatSolve', ' ', Level=4 )
663
       CALL Info( 'HeatSolve', '-------------------------------------',Level=4 )
664
       WRITE( Message,* ) 'Enthalpy_h ITERATION', iter
665
       CALL Info( 'HeatSolve', Message, Level=4 )
666
       CALL Info( 'HeatSolve', '-------------------------------------',Level=4 )
667
       CALL Info( 'HeatSolve', ' ', Level=4 )
668
       CALL Info( 'HeatSolve', 'Starting Assembly...', Level=4 )
669

670
       IF ( ConstantBulk .AND. ASSOCIATED(Solver % Matrix % BulkValues) ) THEN
671
         Solver % Matrix % Values = Solver % Matrix % BulkValues
672
         Solver % Matrix % RHS = Solver % Matrix % BulkRHS
673
         GOTO 1000
674
       END IF
675

676
!------------------------------------------------------------------------------
677
       CALL DefaultInitialize()
678
!------------------------------------------------------------------------------
679

680

681

682

683
       IF ( SmartHeaterControl .OR. IntegralHeaterControl ) THEN
684
          IF( SmartHeaterControl) ForceHeater = 0.0d0
685
          HeaterArea = 0.0d0
686
          HeaterSource = 0.0d0
687
          HeaterScaling = 1.0d0
688
          HeaterDensity = 0.0d0
689
          HeaterTarget = 0.0d0
690
          HeaterControlLocal = .FALSE.
691

692
          DO t=1,Solver % NumberOfActiveElements             
693
             Element => GetActiveElement(t)             
694
             BodyForce => GetBodyForce()
695
             
696
             IF ( .NOT. ASSOCIATED( BodyForce ) ) CYCLE
697
             bf_id = GetBodyForceId()
698
             
699
             IF( .NOT. (SmartHeaters(bf_id) .OR. IntegralHeaters(bf_id) ) ) CYCLE
700

701
             n = GetElementNOFNodes()
702

703
             Material => GetMaterial()
704

705
             Density(1:n) = GetReal( Material, 'Enthalpy Density' )
706
             Load(1:n) = GetReal( BodyForce, 'Heat Source' )
707

708
             s = ElementArea( Solver % Mesh, Element, n )
709

710
             IF( CurrentCoordinateSystem() == AxisSymmetric .OR. &
711
                  CurrentCoordinateSystem() == CylindricSymmetric ) s = 2 * PI * s
712

713
             HeaterSource(bf_id) = HeaterSource(bf_id) + s * SUM(Density(1:n) * Load(1:n)) / n
714
             HeaterArea(bf_id) = HeaterArea(bf_id) + s
715
             HeaterDensity(bf_id) = HeaterDensity(bf_id) + s * SUM( Density(1:n) ) / n
716
          END DO
717

718
          DO i = 1,Model % NumberOfBodyForces
719
             IF( IntegralHeaters(i) .OR. SmartHeaters(i) ) THEN
720
                HeaterDensity(i) = HeaterDensity(i) / HeaterArea(i)
721
             END IF
722
             IF(IntegralHeaters(i)) THEN
723
                HeaterTarget(i) = GetCReal(  Model % BodyForces(i) % Values, &
724
                     'Integral Heat Source', Found )
725
                HeaterScaling(i) = HeaterTarget(i) / HeaterSource(i)
726
             END IF
727
          END DO
728
       END IF
729

730
!------------------------------------------------------------------------------
731
       body_id = -1
732
       NULLIFY(Material)
733
!------------------------------------------------------------------------------
734
!      Bulk elements
735
!------------------------------------------------------------------------------
736
       CALL StartAdvanceOutput( 'HeatSolve', 'Assembly:' )
737
       DO t=1,Solver % NumberOfActiveElements
738

739
         CALL AdvanceOutput(t,GetNOFActive())
740
!------------------------------------------------------------------------------
741
!        Check if this element belongs to a body where Enthalpy_h 
742
!        should be calculated
743
!------------------------------------------------------------------------------
744
         Element => GetActiveElement(t)
745

746
!------------------------------------------------------------------------------
747
         IF ( Element % BodyId /= body_id ) THEN
748
!------------------------------------------------------------------------------
749
           Equation => GetEquation()
750
           ConvectionFlag = GetString( Equation, 'Convection', Found )
751

752
           Material => GetMaterial()
753
!------------------------------------------------------------------------------
754
           CompressibilityFlag = GetString( Material, &
755
                 'Compressibility Model', Found)
756
           IF ( .NOT.Found ) CompressibilityModel = Incompressible
757

758
           SELECT CASE( CompressibilityFlag )
759

760
             CASE( 'incompressible' )
761
               CompressibilityModel = Incompressible
762

763
             CASE( 'user defined' )
764
               CompressibilityModel = UserDefined1
765

766
             CASE( 'perfect gas', 'perfect gas equation 1' )
767
               CompressibilityModel = PerfectGas1
768

769
             CASE( 'thermal' )
770
               CompressibilityModel = Thermal
771

772
             CASE DEFAULT
773
               CompressibilityModel = Incompressible
774
           END SELECT
775
!------------------------------------------------------------------------------
776

777
           PhaseModel = GetString( Equation, 'Phase Change Model',Found )
778
           IF(.NOT. Found) PhaseModel = GetString( Material, 'Phase Change Model',Found )
779

780
           PhaseChange = Found .AND. (PhaseModel(1:4) /= 'none')
781
           IF ( PhaseChange ) THEN
782
              CheckLatentHeatRelease = GetLogical( Equation, &
783
                   'Check Latent Heat Release',Found )
784
           END IF
785
         END IF
786
!------------------------------------------------------------------------------
787

788
         n = GetElementNOFNodes()
789
         CALL GetElementNodes( ElementNodes )
790

791
         CALL GetScalarLocalSolution( LocalEnthalpy_h )
792
!------------------------------------------------------------------------------
793
!        Get element material parameters
794
!------------------------------------------------------------------------------
795
         HeatCapacity(1:n) = 1.0 
796
 	 WaterDiffusivity(1:n) = GetReal( Material, 'Enthalpy Water Diffusivity', Found )
797

798
IF (.NOT.Found) THEN
799
         CALL FATAL('EnthalpySolver', 'No value for >Enthalpy Water Diffusivity< defined in material.')
800
ENDIF
801
	
802

803
         CALL ListGetRealArray( Material,'Enthalpy Heat Diffusivity',Hwrk,n, &
804
                      Element % NodeIndexes )
805
					 
806
         HeatConductivity = 0.0d0
807
         IF ( SIZE(Hwrk,1) == 1 ) THEN
808
           DO i=1,3
809
		DO j=1,n
810
			IF (PhaseChangeEnthValues(PhaseChangeEnthPerm(Element % NodeIndexes(j)))<&
811
				&Enthalpy_h(EnthalpyPerm(Element % NodeIndexes(j)))) then
812
             			HeatConductivity( i,i,j ) = WaterDiffusivity(j)
813
			ELSE
814
				HeatConductivity( i,i,j ) = Hwrk( 1,1,j )
815
			ENDIF
816
		ENDDO
817
           END DO
818
         ELSE IF ( SIZE(Hwrk,2) == 1 ) THEN
819
           DO i=1,MIN(3,SIZE(Hwrk,1))
820
		DO j=1,n
821
			IF (PhaseChangeEnthValues(PhaseChangeEnthPerm(Element % NodeIndexes(j)))<&
822
				&Enthalpy_h(EnthalpyPerm(Element % NodeIndexes(j)))) then
823
				HeatConductivity(i,i,j) = WaterDiffusivity(j)
824
			ELSE
825
			        HeatConductivity(i,i,j) = Hwrk(i,1,j)
826
			ENDIF
827
		ENDDO
828
           END DO
829
         ELSE
830
           DO i=1,MIN(3,SIZE(Hwrk,1))
831
             DO j=1,MIN(3,SIZE(Hwrk,2))
832
		DO k=1,n
833
			IF (PhaseChangeEnthValues(PhaseChangeEnthPerm(Element % NodeIndexes(k)))<&
834
				&Enthalpy_h(EnthalpyPerm(Element % NodeIndexes(k)))) then
835
				HeatConductivity( i,j,k ) = WaterDiffusivity(k)
836
			ELSE
837
				HeatConductivity( i,j,k ) = Hwrk(i,j,k)
838
			ENDIF
839

840
		ENDDO
841
             END DO
842
           END DO
843
         END IF
844
!------------------------------------------------------------------------------
845

846
         IF ( CompressibilityModel == PerfectGas1 ) THEN
847

848
           ! Read Specific Heat Ratio:
849
           !--------------------------
850
           SpecificHeatRatio = GetConstReal( Material, &
851
               'Specific Heat Ratio', Found )
852
           IF ( .NOT.Found ) SpecificHeatRatio = 5.d0/3.d0
853

854
           ! For an ideal gas, \gamma, c_p and R are really a constant
855
           ! GasConstant is an array only since HeatCapacity formally is:
856
           !-------------------------------------------------------------
857
           GasConstant(1:n) = ( SpecificHeatRatio - 1.d0 ) * &
858
               HeatCapacity(1:n) / SpecificHeatRatio
859

860
           PressureCoeff(1:n) = GetReal( Material, 'Pressure Coefficient', Found )
861
           IF ( .NOT. Found ) PressureCoeff(1:n) = 1.0_dp
862
         ELSE IF ( CompressibilityModel == Thermal ) THEN
863
           ReferenceEnthalpy_h(1:n) = GetReal( Material, 'Reference Enthalpy_h' )
864
           HeatExpansionCoeff(1:n) = GetReal( Material, 'Heat Expansion Coefficient' )
865

866
           Density(1:n) = GetReal( Material,'Enthalpy Density' )
867
           Density(1:n) = Density(1:n) * ( 1 - HeatExpansionCoeff(1:n)  * &
868
                ( LocalEnthalpy_h(1:n) - ReferenceEnthalpy_h(1:n) ) )
869

870
           PressureCoeff(1:n) = GetReal( Material, 'Pressure Coefficient', Found )
871
           IF ( .NOT. Found ) &
872
             PressureCoeff(1:n) = LocalEnthalpy_h(1:n) * HeatExpansionCoeff(1:n) / &
873
                   ( 1-HeatExpansionCoeff(1:n)*( &
874
                               LocalEnthalpy_h(1:n)-ReferenceEnthalpy_h(1:n)) )
875
         ELSE IF ( CompressibilityModel == UserDefined1 ) THEN
876
           IF ( ASSOCIATED( DensitySol ) ) THEN
877
             CALL GetScalarLocalSolution( Density, 'Enthalpy Density' ) 
878
           ELSE
879
             Density(1:n) = GetReal( Material,'Enthalpy Density' )
880
           END IF
881
           PressureCoeff(1:n) = GetReal( Material, 'Pressure Coefficient', Found )
882
           IF ( .NOT. Found ) PressureCoeff(1:n) = 0.0_dp
883
         ELSE
884
           PressureCoeff(1:n) = GetReal( Material, 'Pressure Coefficient', Found )
885
           IF ( .NOT. Found ) PressureCoeff(1:n) = 0.0_dp
886
           Density(1:n) = GetReal( Material, 'Enthalpy Density' )
887
         END IF
888

889
!------------------------------------------------------------------------------
890
! Take pressure deviation p_d as the dependent variable, p = p_0 + p_d
891
! for PerfectGas, read p_0
892
!------------------------------------------------------------------------------
893
         IF ( CompressibilityModel /= Incompressible ) THEN
894
           ReferencePressure = ListGetConstReal( Material, &
895
               'Reference Pressure', Found,UnFoundFatal=UnFoundFatal)
896
            !Previous default value: ReferencePressure = 0.0d0
897
         END IF
898
!------------------------------------------------------------------------------
899

900
         HeaterControlLocal = .FALSE.
901
         Load = 0.0D0
902
         Pressure = 0.0d0
903
         dPressuredt = 0.0d0
904
!------------------------------------------------------------------------------
905
!        Check for convection model
906
!------------------------------------------------------------------------------
907
         C1 = 1.0D0
908
         U = 0._dp
909
         V = 0._dp
910
         W = 0._dp
911

912
         MU = 0.0d0
913
         CALL GetVectorLocalSolution( MU, 'Mesh Velocity' )
914

915
         IF ( ConvectionFlag == 'constant' ) THEN
916

917
           U(1:n) = GetReal( Material, 'Convection Velocity 1', Found )
918
           IF ( .NOT. Found ) &
919
              U(1:n) = GetReal( Equation, 'Convection Velocity 1', Found )
920
           V(1:n) = GetReal( Material, 'Convection Velocity 2', Found )
921
           IF ( .NOT. Found ) &
922
             V(1:n) = GetReal( Equation, 'Convection Velocity 2', Found )
923
           W(1:n) = GetReal( Material, 'Convection Velocity 3', Found )
924
           IF ( .NOT. Found ) &
925
             W(1:n) = GetReal( Equation, 'Convection Velocity 3', Found )
926

927
         ELSE IF ( ConvectionFlag == 'computed' .AND. &
928
              ASSOCIATED(FlowSolution) ) THEN
929
           DO i=1,n
930
             k = FlowPerm(Element % NodeIndexes(i))
931
             IF ( k > 0 ) THEN
932
!------------------------------------------------------------------------------
933
               Pressure(i) = FlowSolution(NSDOFs*k) + ReferencePressure
934
               SELECT CASE( CompressibilityModel )
935
                 CASE( PerfectGas1 )
936
                   Density(i)  = Pressure(i) / &
937
                       ( GasConstant(i) * LocalEnthalpy_h(i) )
938
               END SELECT
939
               IF ( TransientSimulation ) THEN
940
                 dPressureDt(i) = ( FlowSolution(NSDOFs*k) - &
941
                     FlowSol % PrevValues(NSDOFs*k,1) ) / dt
942
               END IF
943
!------------------------------------------------------------------------------
944

945
               SELECT CASE( NSDOFs )
946
               CASE(3)
947
                 U(i) = FlowSolution( NSDOFs*k-2 )
948
                 V(i) = FlowSolution( NSDOFs*k-1 )
949
                 W(i) = 0.0D0
950

951
               CASE(4)
952
                 U(i) = FlowSolution( NSDOFs*k-3 )
953
                 V(i) = FlowSolution( NSDOFs*k-2 )
954
                 W(i) = FlowSolution( NSDOFs*k-1 )
955
               END SELECT
956
             ELSE
957
               U(i) = 0.0d0
958
               V(i) = 0.0d0
959
               W(i) = 0.0d0
960
             END IF
961
           END DO
962
         ELSE IF (ConvectionFlag=='computed' ) THEN
963
           CALL Warn( 'HeatSolver', 'Convection model specified ' //  &
964
                 'but no associated flow field present?' )
965
         ELSE
966
           IF ( ALL(MU==0) ) C1 = 0.0D0 
967
         END IF
968
!------------------------------------------------------------------------------
969
!        Check if modelling Phase Change with Eulerian approach 
970
!------------------------------------------------------------------------------
971
         PhaseSpatial = .FALSE.
972
         IF (  PhaseChange ) THEN
973
           CALL EffectiveHeatCapacity()
974
         ELSE
975
           HeatCapacity(1:n) = Density(1:n) * HeatCapacity(1:n)
976
         END IF
977

978
         Viscosity = 0.0d0
979
!------------------------------------------------------------------------------
980
!        Add body forces, if any
981
!------------------------------------------------------------------------------
982
         BodyForce => GetBodyForce()
983
         IF ( ASSOCIATED( BodyForce ) ) THEN
984
           bf_id = GetBodyForceId()
985
!------------------------------------------------------------------------------
986
!          Frictional viscous heating
987
!------------------------------------------------------------------------------
988
           IF ( GetLogical( BodyForce, 'Friction Heat',Found) ) THEN
989
              Viscosity(1:n) = GetReal( Material,'Viscosity' )
990
           END IF
991
!------------------------------------------------------------------------------
992
!          Given heat source
993
!------------------------------------------------------------------------------
994
           Load(1:n) = Density(1:n) *  GetReal( BodyForce, 'Heat Source', Found )
995
           
996
           IF ( SmartHeaterControl .AND. NewtonLinearization .AND. SmartTolReached) THEN
997
              IF(  SmartHeaters(bf_id) ) THEN
998
               HeaterControlLocal = .TRUE.
999
               IF( TransientHeaterControl ) THEN
1000
                 Load(1:n) = PrevPowerScaling * Load(1:n)
1001
                 HeaterScaling(bf_id) = PrevPowerScaling
1002
               END IF
1003
             END IF
1004
           END IF
1005

1006
           IF ( IntegralHeaterControl ) THEN
1007
              IF( IntegralHeaters(bf_id) ) THEN
1008
                 Load(1:n) = Load(1:n) * HeaterScaling(bf_id) 
1009
              END IF
1010
           END IF
1011
         END IF
1012
	
1013
         C0 = 0.0_dp
1014
!------------------------------------------------------------------------------
1015
! Note at this point HeatCapacity = \rho * c_p OR \rho * (c_p - R)
1016
! and C1 = 0 (diffusion) or 1 (convection)
1017
!------------------------------------------------------------------------------
1018
	 
1019
!------------------------------------------------------------------------------
1020
!          Perfusion (added as suggested by Matthias Zenker)
1021
!------------------------------------------------------------------------------
1022
         PerfusionRate(1:n) = GetReal( BodyForce, 'Perfusion Rate', Found )
1023
         
1024
         IF ( Found ) THEN
1025
           PerfusionRefEnthalpy_h(1:n) = GetReal( BodyForce, 'Perfusion Reference Enthalpy_h' )
1026
           PerfusionDensity(1:n) = GetReal( BodyForce, 'Perfusion Density' )
1027
           PerfusionHeatCapacity(1:n) = GetReal( BodyForce, 'Perfusion Heat Capacity' )
1028
           C0(1:n) = PerfusionHeatCapacity(1:n) * PerfusionRate(1:n) * PerfusionDensity(1:n) 
1029
           Load(1:n) = Load(1:n) + C0(1:n) * PerfusionRefEnthalpy_h(1:n)           
1030
         END IF
1031

1032
!------------------------------------------------------------------------------
1033
!        Get element local matrices, and RHS vectors
1034
!------------------------------------------------------------------------------
1035
         IF ( CurrentCoordinateSystem() == Cartesian ) THEN
1036
!------------------------------------------------------------------------------
1037
           CALL DiffuseConvectiveCompose( &
1038
               MASS, STIFF, FORCE, LOAD, &
1039
               HeatCapacity, C0, C1*HeatCapacity(1:n), HeatConductivity, &
1040
               PhaseSpatial, LocalEnthalpy_h, Enthalpy, U, V, W, &
1041
               MU(1,1:n),MU(2,1:n),MU(3,1:n), Viscosity, Density, Pressure, &
1042
               dPressureDt, PressureCoeff, CompressibilityModel /= Incompressible, &
1043
               Stabilize, UseBubbles, Element, n, ElementNodes )
1044

1045
!------------------------------------------------------------------------------
1046
         ELSE
1047
!------------------------------------------------------------------------------
1048
           CALL DiffuseConvectiveGenCompose( &
1049
               MASS, STIFF, FORCE, LOAD, &
1050
               HeatCapacity, C0, C1*HeatCapacity(1:n), HeatConductivity, &
1051
               PhaseSpatial, LocalEnthalpy_h, Enthalpy, U, V, W, &
1052
               MU(1,1:n),MU(2,1:n),MU(3,1:n), Viscosity, Density, Pressure, &
1053
               dPressureDt, PressureCoeff, CompressibilityModel /= Incompressible, &
1054
               Stabilize, Element, n, ElementNodes )
1055
!------------------------------------------------------------------------------
1056
         END IF
1057
!------------------------------------------------------------------------------
1058

1059
         IF ( HeaterControlLocal .AND. .NOT. TransientHeaterControl) THEN
1060

1061
           IF ( TransientAssembly .AND. .NOT. ConstantBulk ) THEN
1062
             CALL Default1stOrderTime( MASS, STIFF, FORCE )
1063
           END IF
1064

1065
           CALL UpdateGlobalEquations( Solver % Matrix, STIFF, &
1066
               ForceHeater, FORCE, n, 1, EnthalpyPerm(Element % NodeIndexes) )
1067
         ELSE
1068
            Bubbles = UseBubbles .AND. .NOT.Stabilize .AND. &
1069
            ( ConvectionFlag == 'computed' .OR. ConvectionFlag == 'constant' )
1070
            
1071
!------------------------------------------------------------------------------
1072
!           If time dependent simulation add mass matrix to stiff matrix
1073
!------------------------------------------------------------------------------
1074
            TimeForce  = 0.0_dp
1075
            IF ( TransientAssembly ) THEN
1076
               IF ( ConstantBulk ) THEN
1077
                 CALL DefaultUpdateMass( MASS )
1078
               ELSE
1079
                 CALL Default1stOrderTime( MASS,STIFF,FORCE )
1080
               END IF
1081
            ELSE IF ( Solver % NOFEigenValues>0 ) THEN
1082
              CALL DefaultUpdateDamp(MASS)
1083
            END IF
1084
!------------------------------------------------------------------------------
1085
!           Update global matrices from local matrices
1086
!------------------------------------------------------------------------------
1087
            IF (  Bubbles ) THEN
1088
               CALL Condensate( N, STIFF, FORCE, TimeForce )
1089
            END IF
1090

1091
            CALL DefaultUpdateEquations( STIFF, FORCE )
1092
         END IF
1093
!------------------------------------------------------------------------------
1094
      END DO     !  Bulk elements
1095
!------------------------------------------------------------------------------
1096

1097
      CALL DefaultFinishBulkAssembly()
1098

1099

1100
1000  CONTINUE
1101

1102
!------------------------------------------------------------------------------
1103
!     Neumann & Newton boundary conditions
1104
!------------------------------------------------------------------------------
1105
      DO t=1, Solver % Mesh % NumberOfBoundaryElements
1106
        Element => GetBoundaryElement(t)
1107
        IF ( .NOT. ActiveBoundaryElement() ) CYCLE
1108

1109
        n = GetElementNOFNodes()
1110

1111
        ! Check that the dimension of element is suitable for fluxes
1112
        IF( .NOT. PossibleFluxElement(Element) ) CYCLE
1113

1114
        BC => GetBC()
1115
        IF ( .NOT. ASSOCIATED(BC) ) CYCLE
1116

1117
        ! This checks whether there are any Dirichlet conditions on the 
1118
        ! smart heater boundary. If there are the r.h.s. must be zero as 
1119
        ! there can possibly not be any effect on Enthalpy_h.
1120
        !-----------------------------------------------------------------
1121
	IF ( HeaterControlLocal .AND. .NOT. TransientHeaterControl) THEN
1122
          IF( ListCheckPresent(BC, Varname) ) THEN
1123
             nd = GetElementDOFs(Indexes)
1124
             ForceHeater(EnthalpyPerm(Indexes(1:nd))) = 0.0_dp
1125
          END IF
1126
        END IF
1127

1128
        HeatFluxBC = GetLogical( BC, 'Enthalpy Heat Flux BC', Found )
1129
        IF ( Found .AND. .NOT. HeatFluxBC ) CYCLE
1130

1131
        HeatGapBC = ListGetLogical( BC, 'Heat Gap', Found )
1132
        CALL AddHeatFluxBC()
1133

1134
        IF ( HeatGapBC ) THEN
1135
          CALL FindGapIndexes( Element, Indexes, n )
1136
          SaveIndexes(1:n) = Element % NodeIndexes
1137
          Element % NodeIndexes = Indexes(1:n)
1138
          CALL AddHeatFluxBC()
1139
          Element % NodeIndexes = SaveIndexes(1:n)
1140
        END IF
1141

1142
      END DO   ! Neumann & Newton BCs
1143
!------------------------------------------------------------------------------
1144

1145
      IF ( TransientSimulation .AND. ConstantBulk ) CALL AddGlobalTime()
1146

1147
      CALL DefaultFinishAssembly()
1148
      CALL Info( 'HeatSolve', 'Assembly done', Level=4 )
1149

1150
      CALL DefaultDirichletBCs()
1151

1152
!------------------------------------------------------------------------------
1153

1154
!------------------------------------------------------------------------------
1155
!     Solve the system and check for convergence
1156
!------------------------------------------------------------------------------
1157
      at = CPUTime() - at
1158
      st = CPUTime()
1159

1160
      PrevNorm = Norm
1161

1162
      IF(SmartHeaterControl .AND. NewtonLinearization .AND. SmartTolReached) THEN
1163
      
1164
        IF(.NOT. TransientHeaterControl) THEN
1165
 
1166
          CALL ListAddLogical(SolverParams, &
1167
              'Skip Compute Nonlinear Change',.TRUE.)
1168

1169
          Relax = GetCReal( SolverParams, &
1170
              'Nonlinear System Relaxation Factor', Found )
1171
          
1172
          IF ( Found .AND. Relax /= 1.0d0 ) THEN
1173
            CALL ListAddConstReal( Solver % Values, &
1174
                'Nonlinear System Relaxation Factor', 1.0d0 )
1175
          ELSE
1176
            Relax = 1.0d0
1177
          END IF          
1178

1179
          CALL SolveSystem( Solver % Matrix, ParMatrix, &
1180
              ForceHeater, XX, Norm, 1, Solver )
1181
         
1182
          CALL SolveSystem( Solver % Matrix, ParMatrix, &
1183
              Solver % Matrix % RHS, YY, Norm, 1, Solver )
1184

1185
          CALL ListAddLogical(SolverParams,'Skip Compute Nonlinear Change',.FALSE.)
1186
        ELSE          
1187
          CALL SolveSystem( Solver % Matrix, ParMatrix, &
1188
              Solver % Matrix % RHS, Enthalpy_h, Norm, 1, Solver )
1189
          YY = Enthalpy_h
1190
        END IF
1191

1192
        IF(.NOT. SmartHeaterAverage) THEN
1193
          xave = XX(EnthalpyPerm(SmartHeaterNode))
1194
          yave = YY(EnthalpyPerm(SmartHeaterNode))
1195
        ELSE          
1196
          xave = 0.0d0
1197
          yave = 0.0d0
1198
          j = 0
1199
          
1200
          DO k = Model % Mesh % NumberOfBulkElements + 1, &
1201
              Model % Mesh % NumberOfBulkElements + Model % Mesh % NumberOfBoundaryElements            
1202

1203
            Element => Model % Mesh % Elements(k)            
1204
            IF ( Element % BoundaryInfo % Constraint == SmartHeaterBC ) THEN
1205
              l = Element % TYPE % NumberOfNodes
1206
              j = j + l
1207
              xave = xave + SUM( XX(EnthalpyPerm(Element % NodeIndexes)) )
1208
              yave = yave + SUM( YY(EnthalpyPerm(Element % NodeIndexes)) )
1209
            END IF
1210
          END DO
1211
          xave = xave / j
1212
          yave = yave / j 
1213
          CALL ListAddConstReal(Model % Simulation,'res: Smart Heater Enthalpy_h',yave)
1214
        END IF
1215

1216
        IF(.NOT. TransientHeaterControl) THEN
1217
          IF ( ASSOCIATED(Solver % Variable % NonlinValues) ) THEN
1218
            Solver % Variable % NonlinValues = Enthalpy_h
1219
          END IF
1220

1221
          PowerScaling = (MeltPoint - yave) / xave 
1222
          Enthalpy_h = YY + PowerScaling * XX
1223

1224
          ! The change is computed separately for the controlled Enthalpy_h field
1225
          !-----------------------------------------------------------------------
1226
          CALL ComputeChange(Solver,.FALSE.,LocalNodes,Enthalpy_h)
1227
          Norm = Solver % Variable % Norm
1228

1229
        END IF
1230

1231
        IF(dt > PowerTimeScale) THEN
1232
          IF ( Relax /= 1.0d0 ) THEN
1233
            CALL ListAddConstReal( Solver % Values,  &
1234
                'Nonlinear System Relaxation Factor', Relax )
1235
          END IF
1236
        END IF
1237
      ELSE
1238
        Norm = DefaultSolve()
1239
      END IF
1240

1241

1242
      IF( SmartHeaterControl .OR. IntegralHeaterControl) THEN
1243
         
1244
         CALL ListAddConstReal(Model % Simulation,'res: Heater Power Scaling',PowerScaling)
1245
         
1246
         CALL Info( 'HeatSolve', 'Heater Control Information', Level=4 )
1247
         DO i=1,Model % NumberOfBodyForces
1248
            IF( .NOT. (SmartHeaters(i) .OR. IntegralHeaters(i))) CYCLE
1249
            IF( SmartHeaters(i) )  HeaterScaling(i) = PowerScaling
1250

1251
            WRITE( Message, '(A,T35,I15)' ) 'Heater for body: ', i
1252
            CALL Info( 'HeatSolve', Message, Level=4 )
1253
            IF(SmartHeaters(i)) WRITE( Message, '(A,T35,A)' ) 'Heater type:','Smart heater'
1254
            IF(IntegralHeaters(i)) WRITE( Message, '(A,T35,A)' ) 'Heater type:','Integral heater'
1255
            CALL Info( 'HeatSolve', Message, Level=4 )
1256

1257
            WRITE( Message,'(A,T35,ES15.4)' ) 'Heater Volume (m^3): ', HeaterArea(i)
1258
            CALL Info( 'HeatSolve', Message, Level=4 )
1259
            s = HeaterSource(i) * HeaterScaling(i)
1260
            WRITE( Message,'(A,T35,ES15.4)' ) 'Heater Power (W): ', s
1261
            CALL Info( 'HeatSolve', Message, Level=4 )
1262

1263
            WRITE( Message,'(A,T35,ES15.4)' ) 'Heater scaling: ', HeaterScaling(i)
1264
            CALL Info( 'HeatSolve', Message, Level=4 )
1265
            WRITE( Message, '(A,T35,ES15.4)' ) 'Heater Power Density (W/kg): ', s/(HeaterDensity(i) * HeaterArea(i))
1266
            CALL Info( 'HeatSolve', Message, Level=4 )
1267
            
1268
            IF( SmartHeaters(i)) CALL ListAddConstReal(Model % Simulation,'res: Heater Power Density',&
1269
                 s/(HeaterDensity(i) * HeaterArea(i)))
1270
         END DO
1271
      END IF
1272

1273

1274
      st = CPUTIme()-st
1275
      totat = totat + at
1276
      totst = totst + st
1277
      WRITE(Message,'(a,i4,a,F8.2,F8.2)') 'iter: ',iter,' Assembly: (s)', at, totat
1278
      CALL Info( 'HeatSolve', Message, Level=4 )
1279
      WRITE(Message,'(a,i4,a,F8.2,F8.2)') 'iter: ',iter,' Solve:    (s)', st, totst
1280
      CALL Info( 'HeatSolve', Message, Level=4 )
1281
!------------------------------------------------------------------------------
1282
!     If modelling phase change (and if requested by the user), check if any
1283
!     node has jumped over the phase change interval, and if so, reduce
1284
!     timestep and or relaxation and recompute.
1285
!------------------------------------------------------------------------------
1286
      IF (PhaseChange .AND. CheckLatentHeatRelease ) THEN
1287
!------------------------------------------------------------------------------
1288
        IF ( CheckLatentHeat() ) THEN
1289
          Enthalpy_h(1:LocalNodes) = PrevSolution
1290
          Norm = PrevNorm
1291

1292
          IF ( TransientSimulation ) THEN
1293
            dt = dt / 2
1294
            Solver % dt = dt
1295
            WRITE( Message, * ) &
1296
                  'Latent heat release check: reducing timestep to: ',dt
1297
            CALL Info( 'HeatSolve', Message, Level=4 )
1298
          ELSE
1299
            Relax = Relax / 2
1300
            CALL  ListAddConstReal( Solver % Values,  &
1301
                 'Nonlinear System Relaxation Factor', Relax )
1302
            WRITE( Message, * ) &
1303
                 'Latent heat release check: reducing relaxation to: ',Relax
1304
            CALL Info( 'HeatSolve', Message, Level=4 )
1305
          END IF
1306

1307
          CYCLE
1308
        END IF
1309
        IF ( .NOT.TransientSimulation ) PrevSolution=Enthalpy_h(1:LocalNodes)
1310
      END IF
1311
!------------------------------------------------------------------------------
1312
     
1313
      RelativeChange = Solver % Variable % NonlinChange
1314

1315
      WRITE( Message, * ) 'Result Norm   : ',Norm
1316
      CALL Info( 'HeatSolve', Message, Level=4 )
1317
      WRITE( Message, * ) 'Relative Change : ',RelativeChange
1318
      CALL Info( 'HeatSolve', Message, Level=4 )
1319

1320
      IF ( RelativeChange < NewtonTol .OR. iter >= NewtonIter ) &
1321
               NewtonLinearization = .TRUE.
1322
      IF ( RelativeChange < NonlinearTol .AND. &
1323
          (.NOT. SmartHeaterControl .OR. SmartTolReached)) EXIT
1324

1325
      IF(SmartHeaterControl) THEN
1326
        IF ( RelativeChange < SmartTol ) THEN
1327
          SmartTolReached = .TRUE.
1328
          YY = Enthalpy_h
1329
        END IF
1330
      END IF
1331
      
1332
!------------------------------------------------------------------------------
1333
    END DO ! of the nonlinear iteration
1334
!------------------------------------------------------------------------------
1335
    IF(TransientHeaterControl) THEN
1336
      PowerRelax = GetCReal(Solver % Values,'Smart Heater Relaxation Factor', GotIt)
1337
      IF(.NOT. GotIt) PowerRelax = 1.0_dp
1338
      PowerSensitivity = ListGetConstReal(Solver % Values,'Smart Heater Power Sensivity',GotIt,UnFoundFatal=UnFoundFatal)
1339
      !Previous default value: PowerSensitivity = 4.0_dp
1340
      PowerScaling = PowerScaling * (1 + PowerSensitivity * PowerRelax * (MeltPoint/yave - 1.0d0) ) 
1341

1342
      IF( ListGetLogical( Solver % Values,'Smart Heater Transient Speedup',GotIt ) ) THEN
1343
        Enthalpy_h = Enthalpy_h * (1 + PowerRelax * (MeltPoint/yave - 1.0d0)   )     
1344
      END IF
1345
      YY = Enthalpy_h
1346
    END IF
1347

1348

1349
!------------------------------------------------------------------------------
1350
!   Compute temperature homologous and water content
1351
!------------------------------------------------------------------------------
1352

1353
 WaterName = GetString(SolverParams , 'Exported Variable 2', Found )
1354

1355
     IF (.NOT.Found) &
1356
        CALL FATAL('EnthalpySolver','No value > Exported Variable 2 (water content) < found in Solver')
1357

1358
 TempName = GetString(SolverParams , 'Exported Variable 3', Found )
1359

1360
     IF (.NOT.Found) &
1361
        CALL FATAL('EnthalpySolver','No value > Exported Variable 3 (temperature) < found in Solver')
1362

1363
      WaterVar => VariableGet(Solver % Mesh % Variables, WaterName)
1364
      TemphomoVar => VariableGet(Solver % Mesh % Variables,TempName)
1365

1366
      L_heat = GetConstReal(Model % Constants, "L_heat",GotIt)
1367
      IF (.NOT.GotIt) THEN
1368
         CALL WARN('EnthalpySolver', 'No Keyword >L_heat< defined in model constants. Using >334000Jkg-1< as default.')
1369
         L_heat = 334000.0
1370
      END IF
1371
	  
1372
do i=1,Model % NumberOfNodes
1373

1374
hi = Enthalpy_h (EnthalpyPerm (i) )
1375
hm = PhaseChangeEnthValues (PhaseChangeEnthPerm(i))
1376

1377
  if (hi<hm) then
1378
    WaterVar % values ( WaterVar % perm (i) ) = 0.0
1379
    TemphomoVar % values ( TemphomoVar % perm (i) ) = &
1380
    & (-B_cap+(B_cap**2+A_cap*2*(A_cap*0.5*Tref**2+B_cap*Tref+hi))**0.5 ) / A_cap - 273.16 ! Use CP(T)=B_cap+A_cap*T
1381
  else
1382
    TemphomoVar % values ( TemphomoVar % perm (i) ) =  &
1383
    & (-B_cap+(B_cap**2+A_cap*2*(A_cap*0.5*Tref**2+B_cap*Tref+hm))**0.5 ) / A_cap - 273.16 ! Use CP(T)=B_cap+A_cap*T
1384
    WaterVar % values ( WaterVar % perm (i) ) = (hi-hm)/L_heat
1385
  endif
1386

1387

1388

1389
enddo
1390

1391
!------------------------------------------------------------------------------
1392
!   Compute cumulative time done by now and time remaining
1393
!------------------------------------------------------------------------------
1394
    IF ( .NOT. TransientSimulation ) EXIT
1395
    CumulativeTime = CumulativeTime + dt
1396
    dt = Timestep - CumulativeTime
1397

1398
   END DO ! time interval
1399
   Solver % dt = Timestep
1400

1401
!------------------------------------------------------------------------------
1402
   CALL  ListAddConstReal( Solver % Values,  &
1403
        'Nonlinear System Relaxation Factor', SaveRelax )
1404
!------------------------------------------------------------------------------
1405

1406
   DEALLOCATE( PrevSolution )
1407

1408
   IF ( ListGetLogical( Solver % Values, 'Adaptive Mesh Refinement', Found ) ) THEN
1409
     IF ( .NOT. ListGetLogical( Solver % Values, 'Library Adaptivity', Found ) ) THEN
1410
       CALL RefineMesh( Model,Solver,Enthalpy_h,EnthalpyPerm, &
1411
            EnthalpySolver_Inside_Residual, EnthalpySolver_Edge_Residual, EnthalpySolver_Boundary_Residual )
1412
     END IF
1413
   END IF
1414

1415
CONTAINS
1416

1417

1418

1419
!------------------------------------------------------------------------------
1420
   SUBROUTINE AddHeatFluxBC()
1421
!------------------------------------------------------------------------------
1422
      CALL GetElementNodes( ElementNodes )
1423

1424
      HeatTransferCoeff = 0.0D0
1425
      LOAD  = 0.0D0
1426
!------------------------------------------------------------------------------
1427
!     BC: -k@T/@n = \epsilon\sigma(T^4 - Text^4)
1428
!------------------------------------------------------------------------------
1429
      RadiationFlag = GetString( BC, 'Radiation', Found )
1430

1431
      IF ( Found .AND. RadiationFlag(1:4) /= 'none' ) THEN
1432
        NodalEmissivity(1:n) = GetReal(BC, 'Emissivity', Found)
1433
        IF(.NOT. Found) THEN
1434
           NodalEmissivity(1:n) = GetParentMatProp( 'Emissivity' )
1435
        END IF
1436
        Emissivity = SUM( NodalEmissivity(1:n) ) / n
1437

1438
!------------------------------------------------------------------------------
1439
        IF (  RadiationFlag(1:9) == 'idealized' ) THEN
1440
          AText(1:n) = GetReal( BC, 'Radiation External Enthalpy_h',Found )
1441
          IF(.NOT. Found) AText(1:n) = GetReal( BC, 'External Enthalpy_h' )
1442
        ELSE
1443
          CALL DiffuseGrayRadiation( Model, Solver, Element, & 
1444
              Enthalpy_h, EnthalpyPerm, ForceVector, VisibleFraction, Text)
1445

1446
          IF( GetLogical( BC, 'Radiation Boundary Open', Found) ) THEN
1447
            AText(1:n) = GetReal( BC, 'Radiation External Enthalpy_h',Found )
1448
            IF(.NOT. Found) AText(1:n) = GetReal( BC, 'External Enthalpy_h' )
1449
            IF( VisibleFraction >= 1.0_dp ) THEN
1450
              Atext(1:n) = Text
1451
            ELSE
1452
              Atext(1:n) = ( (1 - VisibleFraction) * Atext(1:n)**4 + &
1453
                  VisibleFraction * Text**4 ) ** 0.25_dp
1454
            END IF
1455
          ELSE
1456
            AText(1:n) = Text
1457
          END IF
1458
        END IF
1459
!------------------------------------------------------------------------------
1460
!       Add our own contribution to surface Enthalpy_h (and external
1461
!       if using linear type iteration or idealized radiation)
1462
!------------------------------------------------------------------------------
1463
        DO j=1,n
1464
          k = EnthalpyPerm(Element % NodeIndexes(j))
1465
          Text = AText(j)
1466

1467
          IF ( .NOT. HeatGapBC .AND. NewtonLinearization ) THEN
1468
             HeatTransferCoeff(j) = Emissivity * 4*Enthalpy_h(k)**3 * &
1469
                               StefanBoltzmann
1470
             LOAD(j) = Emissivity*(3*Enthalpy_h(k)**4+Text**4) * &
1471
                               StefanBoltzmann
1472
          ELSE
1473
             HeatTransferCoeff(j) = Emissivity * (Enthalpy_h(k)**3 + &
1474
             Enthalpy_h(k)**2*Text+Enthalpy_h(k)*Text**2 + Text**3) * &
1475
                               StefanBoltzmann 
1476
             LOAD(j) = HeatTransferCoeff(j) * Text
1477
          END IF
1478
        END DO
1479
      END IF  ! of radition
1480
!------------------------------------------------------------------------------
1481

1482
      Work(1:n)  = GetReal( BC, 'Heat Transfer Coefficient',Found )
1483
      IF ( Found ) THEN
1484
       AText(1:n) = GetReal( BC, 'External Enthalpy_h',Found )
1485
       DO j=1,n
1486
!------------------------------------------------------------------------------
1487
!         BC: -k@T/@n = \alpha(T - Text)
1488
!------------------------------------------------------------------------------
1489
          k = EnthalpyPerm(Element % NodeIndexes(j))
1490
          LOAD(j) = LOAD(j) + Work(j) * AText(j)
1491
          HeatTransferCoeff(j) = HeatTransferCoeff(j) + Work(j)
1492
        END DO
1493
      END IF
1494

1495
!------------------------------------------------------------------------------
1496
!     BC: -k@T/@n = (rho*L)*v.n 
1497
!     Heating related to pulling is possible only in ss cases where pull velocity
1498
!     is desrcibed.
1499
!------------------------------------------------------------------------------
1500

1501
      IF( GetLogical( BC, 'Phase Change',Found ) ) THEN
1502
         PhaseVelocity(1,1:n) = GetReal( BC,'Phase Velocity 1', Found  )
1503
         PhaseVelocity(2,1:n) = GetReal( BC,'Phase Velocity 2', Found  )
1504
         PhaseVelocity(3,1:n) = GetReal( BC,'Phase Velocity 3', Found  )
1505
  
1506
         ! Ensure that the latent heat and density come from the same side
1507
         LatentHeat(1:n) = GetParentMatProp( 'Latent Heat', &
1508
              UElement = Element, UParent = Parent )
1509
         IF(.NOT. ASSOCIATED(Parent) ) THEN
1510
           CALL Warn('HeatSolve','Parent not associated')
1511
         ELSE
1512
           k = GetInteger(Model % Bodies(Parent % BodyId) % Values,'Material')
1513
           Density(1:n) = GetReal( Model % Materials(k) % Values, 'Enthalpy Density' )
1514
         END IF
1515

1516
         ! This could be rather put as a new type of BC into the assembly routine and 
1517
         ! then the Normal could be taken at the proper Gaussian integration points. 
1518
         Normal = NormalVector( Element, ElementNodes, 0.0_dp, 0.0_dp, .TRUE. )
1519

1520
         DO i=1,n
1521
            LOAD(i) = LOAD(i) + &
1522
                 LatentHeat(i) * Density(i) * SUM( Normal(1:3) * PhaseVelocity(1:3,i))
1523
         END DO
1524
      END IF
1525

1526
!------------------------------------------------------------------------------
1527
!     BC: -k@T/@n = g
1528
!------------------------------------------------------------------------------
1529
      LOAD(1:n) = LOAD(1:n) +  GetReal( BC, 'Enthalpy Heat Flux', Found )
1530
      WaterDiffusivity(1:n) = GetReal( Material, 'Enthalpy Water Diffusivity', Found )
1531

1532
      InfBC = ListGetLogical( BC,'Infinity BC '//TRIM(VarName),GotIt)
1533
      IF( InfBC ) THEN
1534
        AText(1:n) = GetReal( BC,'Infinity BC '//TRIM(VarName)//' Offset',GotIt)
1535
        ! currently only isotropic heat conductivity supported
1536
        HeatConductivityIso(1:n) = GetParentMatProp('Enthalpy Heat Diffusivity',Element,GotIt)
1537

1538
		DO k=1,n
1539
			IF (PhaseChangeEnthValues(PhaseChangeEnthPerm(Element % NodeIndexes(k)))<&
1540
				&Enthalpy_h(EnthalpyPerm(Element % NodeIndexes(k)))) then
1541
				HeatConductivityIso(k) = WaterDiffusivity(k)
1542
			ENDIF
1543
		ENDDO
1544

1545

1546
        IF(.NOT. GotIt) THEN
1547
          CALL Fatal( 'HeatSolver','Could not find > Heat Conductivity < for parent!' )           
1548
        END IF
1549
      END IF
1550

1551
!------------------------------------------------------------------------------
1552
!     Get element matrix and rhs due to boundary conditions ...
1553
!------------------------------------------------------------------------------
1554
      IF ( CurrentCoordinateSystem() == Cartesian ) THEN
1555
        CALL DiffuseConvectiveBoundary( STIFF,FORCE, &
1556
            LOAD,HeatTransferCoeff,InfBC,HeatConductivityIso,AText(1:n),&
1557
            Element,n,ElementNodes )
1558
      ELSE
1559
        IF( InfBC ) THEN
1560
          CALL Fatal('HeatSolver','Infinity BC not implemented only for cartersian case!')
1561
        END IF
1562
        CALL DiffuseConvectiveGenBoundary(STIFF,FORCE,&
1563
            LOAD,HeatTransferCoeff,Element,n,ElementNodes ) 
1564
      END IF
1565

1566
!------------------------------------------------------------------------------
1567
!     Update global matrices from local matrices
1568
!------------------------------------------------------------------------------
1569
      IF ( TransientAssembly .AND. .NOT. ConstantBulk ) THEN
1570
        MASS = 0.d0
1571
        CALL Default1stOrderTime( MASS, STIFF, FORCE )
1572
      END IF
1573

1574
      IF ( HeatGapBC ) &
1575
        CALL AddHeatGap( Solver, Element, STIFF, EnthalpyPerm)
1576

1577
      CALL DefaultUpdateEquations( STIFF, FORCE )
1578
!------------------------------------------------------------------------------
1579
    END SUBROUTINE AddHeatFluxBC
1580
!------------------------------------------------------------------------------
1581

1582

1583
!------------------------------------------------------------------------------
1584
    SUBROUTINE AddGlobalTime()
1585
!------------------------------------------------------------------------------
1586
      INTEGER :: i,j,k,l,n
1587
      REAL(KIND=dp) :: FORCE(1)
1588
      REAL(KIND=dp), POINTER :: SaveValues(:) => NULL()
1589
      SAVE STIFF, MASS, X
1590
      REAL(KIND=dp), ALLOCATABLE :: STIFF(:,:),MASS(:,:), X(:,:)
1591

1592
      IF ( .NOT.ASSOCIATED(Solver % Variable % Values, SaveValues) ) THEN
1593
         IF ( ALLOCATED(STIFF) ) DEALLOCATE( STIFF,MASS,X )
1594
         n = 0
1595
         DO i=1,Solver % Matrix % NumberOfRows
1596
           n = MAX( n,Solver % Matrix % Rows(i+1)-Solver % Matrix % Rows(i) )
1597
         END DO
1598
         k = SIZE(Solver % Variable % PrevValues,2)
1599
         ALLOCATE( STIFF(1,n),MASS(1,n),X(n,k) )
1600
 
1601
         SaveValues => Solver % Variable % Values
1602
      END IF
1603

1604
      DO i=1,Solver % Matrix % NumberOFRows
1605
        n = 0
1606
        DO j=Solver % Matrix % Rows(i),Solver % Matrix % Rows(i+1)-1
1607
          n=n+1
1608
          STIFF(1,n) = Solver % Matrix % Values(j)
1609
          MASS(1,n)  = Solver % Matrix % MassValues(j)
1610
          X(n,:) = Solver % Variable % PrevValues(Solver % Matrix % Cols(j),:)
1611
        END DO
1612
        FORCE(1) = Solver % Matrix % RHS(i)
1613
        Solver % Matrix % Force(i,1) = FORCE(1)
1614
        k = MIN( Solver % DoneTime, Solver % Order )
1615
        CALL BDFLocal( n, dt, MASS, STIFF, FORCE, X, k )
1616

1617
        n = 0
1618
        DO j=Solver % Matrix % Rows(i),Solver % Matrix % Rows(i+1)-1
1619
           n=n+1
1620
          Solver % Matrix % Values(j) = STIFF(1,n)
1621
        END DO
1622
        Solver % Matrix % RHS(i) = FORCE(1)
1623
      END DO
1624
!------------------------------------------------------------------------------
1625
    END SUBROUTINE AddGlobalTime
1626
!------------------------------------------------------------------------------
1627

1628
!------------------------------------------------------------------------------
1629
    SUBROUTINE DiffuseGrayRadiation( Model, Solver, Element,  &
1630
      Enthalpy_h, EnthalpyPerm, ForceVector,AngleFraction, Text)
1631
!------------------------------------------------------------------------------
1632
      TYPE(Model_t)  :: Model
1633
      TYPE(Solver_t) :: Solver
1634
      TYPE(Element_t), POINTER :: Element
1635
      INTEGER :: EnthalpyPerm(:)
1636
      REAL(KIND=dp) :: Enthalpy_h(:), ForceVector(:)
1637
      REAL(KIND=dp) :: AngleFraction, Text
1638
!------------------------------------------------------------------------------
1639
      REAL(KIND=dp) :: Area, Asum
1640
      INTEGER :: i,j,k,l,m,ImplicitFactors
1641
      INTEGER, POINTER :: ElementList(:)
1642
!------------------------------------------------------------------------------
1643
!     If linear iteration compute radiation load
1644
!------------------------------------------------------------------------------
1645

1646

1647
      Asum = 0.0d0
1648
      IF ( .NOT. NewtonLinearization ) THEN
1649
        Text = ComputeRadiationLoad( Model, Solver % Mesh, Element, &
1650
                 Enthalpy_h, EnthalpyPerm, Emissivity, AngleFraction)
1651
      ELSE   !  Full Newton-Raphson solver
1652
!------------------------------------------------------------------------------
1653
!       Go through surfaces (j) this surface (i) is getting
1654
!       radiated from.
1655
!------------------------------------------------------------------------------
1656

1657
        Area  = ElementArea( Solver % Mesh, Element, n )
1658
        ElementList => Element % BoundaryInfo % RadiationFactors % Elements
1659

1660
        DO j=1,Element % BoundaryInfo % RadiationFactors % NumberOfFactors
1661

1662
          RadiationElement => Solver % Mesh % Elements( ElementList(j) )
1663

1664
          Text = ComputeRadiationCoeff(Model,Solver % Mesh,Element,j) / ( Area )
1665
          Asum = Asum + Text
1666
!------------------------------------------------------------------------------
1667
!         Gebhart factors are given elementwise at the center
1668
!         of the element, so take average of nodal Enthalpy_hs
1669
!         (or integrate over surface j)
1670
!------------------------------------------------------------------------------
1671

1672
          k = RadiationElement % TYPE % NumberOfNodes
1673
          ImplicitFactors = Element % BoundaryInfo % RadiationFactors % NumberOfImplicitFactors
1674
          IF(ImplicitFactors == 0) &
1675
              ImplicitFactors = Element % BoundaryInfo % RadiationFactors % NumberOfFactors
1676

1677
          IF(j <= ImplicitFactors) THEN
1678
            
1679
            S = (SUM( Enthalpy_h( EnthalpyPerm( RadiationElement % &
1680
                NodeIndexes))**4 )/k )**(1.0d0/4.0d0)
1681
!------------------------------------------------------------------------------
1682
!         Linearization of the G_jiT^4_j term
1683
!------------------------------------------------------------------------------
1684
            HeatTransferCoeff(1:n) = -4 * Text * S**3 * StefanBoltzmann
1685
            LOAD(1:n) = -3 * Text * S**4 * StefanBoltzmann
1686
!------------------------------------------------------------------------------
1687
!         Integrate the contribution of surface j over surface i
1688
!         and add to global matrix
1689
!------------------------------------------------------------------------------
1690
            CALL IntegOverA( STIFF, FORCE, LOAD, &
1691
                HeatTransferCoeff, Element, n, k, ElementNodes ) 
1692
            
1693
            IF ( TransientAssembly ) THEN
1694
              MASS = 0.d0
1695
              CALL Add1stOrderTime( MASS, STIFF, &
1696
                  FORCE,dt,n,1,EnthalpyPerm(Element % NodeIndexes),Solver )
1697
            END IF
1698
            
1699
            DO m=1,n
1700
              k1 = EnthalpyPerm( Element % NodeIndexes(m) )
1701
              DO l=1,k
1702
                k2 = EnthalpyPerm( RadiationElement % NodeIndexes(l) )
1703
                CALL AddToMatrixElement( StiffMatrix,k1, &
1704
                    k2,STIFF(m,l) )
1705
              END DO
1706
              ForceVector(k1) = ForceVector(k1) + FORCE(m)
1707
            END DO
1708

1709
          ELSE
1710

1711
            S = (SUM( Enthalpy_h( EnthalpyPerm( RadiationElement % &
1712
                NodeIndexes))**4 )/k )
1713
            
1714
            HeatTransferCoeff(1:n) = 0.0d0
1715
            LOAD(1:n) = Text * S * StefanBoltzmann
1716
            
1717
            CALL IntegOverA( STIFF, FORCE, LOAD, &
1718
                HeatTransferCoeff, Element, n, k, ElementNodes ) 
1719
            
1720
            DO m=1,n
1721
              k1 = EnthalpyPerm( Element % NodeIndexes(m) )
1722
              ForceVector(k1) = ForceVector(k1) + FORCE(m)
1723
            END DO
1724
            
1725
          END IF 
1726

1727
        END DO
1728

1729
!------------------------------------------------------------------------------
1730
!       We have already added all external Enthalpy_h contributions
1731
!       to the matrix for the Newton type iteration
1732
!------------------------------------------------------------------------------
1733
        AngleFraction = Asum / Emissivity
1734
        Text = 0.0
1735

1736
      END IF  !  of newton-raphson
1737

1738
    END SUBROUTINE DiffuseGrayRadiation
1739
!------------------------------------------------------------------------------
1740

1741

1742
!------------------------------------------------------------------------------
1743
    SUBROUTINE EffectiveHeatCapacity()
1744
      LOGICAL :: Found, Specific
1745
      REAL(KIND=dp), ALLOCATABLE :: dT(:)
1746

1747
!------------------------------------------------------------------------------
1748
!     See if Enthalpy_h gradient indside the element is large enough 
1749
!     to use  the c_p = SQRT( (dH/dx)^2 / (dT/dx)^2 ), otherwise
1750
!     use c_p = dH/dT, or if in time dependent simulation, use
1751
!     c_p = (dH/dt) / (dT/dt), if requested. 
1752
!------------------------------------------------------------------------------
1753

1754
      SELECT CASE(PhaseModel)
1755
!------------------------------------------------------------------------------
1756
        CASE( 'spatial 1' )
1757
          PhaseChangeModel = PHASE_SPATIAL_1
1758
!------------------------------------------------------------------------------
1759

1760
        CASE( 'spatial 2' )
1761
!------------------------------------------------------------------------------
1762
! Check if local variation of Enthalpy_h is large enough to actually use the
1763
! Spatial 2 model. Should perhaps be scaled to element size (or actually
1764
! compute the gradient, but this will do for now...).
1765
!------------------------------------------------------------------------------
1766
          s = MAXVAL(LocalEnthalpy_h(1:n))-MINVAL(LocalEnthalpy_h(1:n))
1767
          IF ( s < AEPS ) THEN
1768
            PhaseChangeModel = PHASE_SPATIAL_1
1769
          ELSE
1770
            PhaseChangeModel = PHASE_SPATIAL_2
1771
          END IF
1772

1773
!------------------------------------------------------------------------------
1774
! Note that here HeatCapacity is miused for saving dT.
1775
!------------------------------------------------------------------------------
1776
        CASE('temporal')
1777
          IF ( TransientSimulation )  THEN
1778
            ALLOCATE( dT(n) )
1779
            dT(1:n) = Enthalpy_h(EnthalpyPerm(Element % NodeIndexes)) - &
1780
                     PrevEnthalpy_h(EnthalpyPerm(Element % NodeIndexes))
1781

1782
            IF ( ANY(ABS(dT(1:n)) < AEPS) ) THEN
1783
              PhaseChangeModel = PHASE_SPATIAL_1
1784
            ELSE
1785
              PhaseChangeModel = PHASE_TEMPORAL
1786
            END IF
1787
          ELSE
1788
             PhaseChangeModel = PHASE_SPATIAL_1
1789
          END IF
1790

1791
!------------------------------------------------------------------------------
1792
        CASE DEFAULT
1793
          PhaseChangeModel = PHASE_SPATIAL_1
1794

1795
      END SELECT
1796
!------------------------------------------------------------------------------
1797

1798
      PhaseSpatial = ( PhaseChangeModel == PHASE_SPATIAL_2 )
1799
      Specific = ListCheckPresent( Material,'Specific Enthalpy')
1800

1801
!-----------------------------------------------------------------------------
1802
      SELECT CASE( PhaseChangeModel )
1803

1804
!------------------------------------------------------------------------------
1805
! This phase change model is available only for some type of real entries 
1806
! that have an implemented analytical derivation rule.
1807
!-----------------------------------------------------------------------------
1808
      CASE( PHASE_SPATIAL_1 )
1809
        HeatCapacity(1:n) = ListGetReal( Material, &
1810
             'Effective Heat Capacity', n,Element % NodeIndexes, Found )
1811
        IF ( .NOT. Found ) THEN
1812
          IF( Specific ) THEN
1813
            HeatCapacity(1:n) = ListGetDerivValue( Material, &
1814
                'Specific Enthalpy', n,Element % NodeIndexes )
1815
            HeatCapacity(1:n) = Density(1:n) * HeatCapacity(1:n)
1816
          ELSE
1817
            HeatCapacity(1:n) = ListGetDerivValue( Material, &
1818
                'Enthalpy', n,Element % NodeIndexes )            
1819
          END IF
1820
        END IF
1821
          
1822
!---------------------------------------------------------------------------------------
1823
! Note that for the 'spatial 2' model the evaluation of c_p is done in each integration
1824
! point and thus Enthalphy and PhaseSpatial flag are used instead of HeatCapacity directly.
1825
!-----------------------------------------------------------------------------------------
1826
      CASE( PHASE_SPATIAL_2 )
1827
        IF( Specific ) THEN
1828
          Enthalpy(1:n) = ListGetReal(Material,'Specific Enthalpy',n,Element % NodeIndexes)
1829
          Enthalpy(1:n) = Density(1:n) * Enthalpy(1:n)
1830
        ELSE
1831
          Enthalpy(1:n) = ListGetReal(Material,'Enthalpy',n,Element % NodeIndexes)          
1832
        END IF
1833
          
1834
!------------------------------------------------------------------------------
1835
      CASE( PHASE_TEMPORAL )
1836
        ! When retrieving the value of enthalphy on the previous timestep 
1837
        ! the relevant entries of the Enthalpy_h solution in the global vector
1838
        ! are tampered in order to make the ListGetReal command work as wanted. 
1839
        ! 1) Values at current Enthalpy_h     
1840
        !------------------------------------------------------------------------
1841
        IF( Specific ) THEN
1842
          Work(1:n) = ListGetReal( Material,'Specific Enthalpy',n,Element % NodeIndexes )
1843
        ELSE
1844
          Work(1:n) = ListGetReal( Material,'Enthalpy',n,Element % NodeIndexes )
1845
        END IF
1846

1847
        ! 2) Values at previous Enthalpy_h
1848
        Enthalpy_h(EnthalpyPerm(Element % NodeIndexes)) = & 
1849
            PrevEnthalpy_h(EnthalpyPerm(Element % NodeIndexes)) 
1850

1851
        IF( Specific ) THEN
1852
          Work(1:n) = Work(1:n) - ListGetReal( Material,'Specific Enthalpy', &
1853
              n,Element % NodeIndexes )          
1854
          HeatCapacity(1:n) = Density(1:n) * Work(1:n) / dT(1:n)
1855
       ELSE
1856
          Work(1:n) = Work(1:n) - ListGetReal( Material,'Enthalpy', &
1857
              n,Element % NodeIndexes )
1858
          HeatCapacity(1:n) = Work(1:n) / dT(1:n)
1859
        END IF
1860

1861
        ! Revert to current Enthalpy_h
1862
        Enthalpy_h(EnthalpyPerm(Element % NodeIndexes)) = & 
1863
            PrevEnthalpy_h(EnthalpyPerm(Element % NodeIndexes)) + dT(1:n)
1864

1865
!------------------------------------------------------------------------------
1866
      END SELECT
1867
!------------------------------------------------------------------------------
1868
    END SUBROUTINE EffectiveHeatCapacity
1869
!------------------------------------------------------------------------------
1870

1871

1872

1873
!------------------------------------------------------------------------------
1874
    FUNCTION CheckLatentHeat() RESULT(Failure)
1875
!------------------------------------------------------------------------------
1876
      LOGICAL :: Failure, PhaseChange, CheckLatentHeatRelease
1877
      INTEGER :: t, eq_id, body_id
1878
      CHARACTER(LEN=MAX_NAME_LEN) :: PhaseModel
1879
!------------------------------------------------------------------------------
1880

1881
      Failure = .FALSE.
1882
!------------------------------------------------------------------------------
1883
      DO t=1,Solver % Mesh % NumberOfBulkElements
1884
!------------------------------------------------------------------------------
1885
!       Check if this element belongs to a body where Enthalpy_h 
1886
!       has been calculated
1887
!------------------------------------------------------------------------------
1888
        Element => Solver % Mesh % Elements(t)
1889

1890
        NodeIndexes => Element % NodeIndexes
1891
        IF ( ANY( EnthalpyPerm( NodeIndexes ) <= 0 ) ) CYCLE
1892

1893
        body_id = Element % Bodyid
1894
        eq_id = ListGetInteger( Model % Bodies(body_id) % Values, &
1895
            'Equation', minv=1, maxv=Model % NumberOfEquations )
1896

1897
        PhaseModel = ListGetString( Model % Equations(eq_id) % Values, &
1898
                          'Phase Change Model',Found )
1899

1900
        PhaseChange = Found .AND. (PhaseModel(1:4) /= 'none')
1901

1902
        IF ( PhaseChange ) THEN
1903
          CheckLatentHeatRelease = ListGetLogical(Model % Equations(eq_id) % &
1904
                    Values, 'Check Latent Heat Release',Found )
1905
        END IF
1906
        IF ( .NOT. ( PhaseChange .AND. CheckLatentHeatRelease ) ) CYCLE
1907

1908
        n = Element % TYPE % NumberOfNodes
1909
!------------------------------------------------------------------------------
1910
!       Set the current element pointer in the model structure to
1911
!       reflect the element being processed
1912
!------------------------------------------------------------------------------
1913
        Model % CurrentElement => Element
1914
!------------------------------------------------------------------------------
1915
!------------------------------------------------------------------------------
1916
!       Get element material parameters
1917
!------------------------------------------------------------------------------
1918
        k = ListGetInteger( Model % Bodies(body_id) % Values,'Material', &
1919
                minv=1, maxv=Model % NumberOfMaterials )
1920
        Material => Model % Materials(k) % Values
1921

1922
        PhaseChangeIntervals => ListGetConstRealArray( Material, &
1923
                        'Phase Change Intervals' )
1924

1925
        DO k=1,n
1926
          i = EnthalpyPerm( NodeIndexes(k) )
1927
          DO j=1,SIZE(PhaseChangeIntervals,2)
1928
            IF ( ( Enthalpy_h(i)  < PhaseChangeIntervals(1,j) .AND. &
1929
                   PrevSolution(i) > PhaseChangeIntervals(2,j) ).OR. &
1930
                 ( Enthalpy_h(i)  > PhaseChangeIntervals(2,j) .AND. &
1931
                   PrevSolution(i) < PhaseChangeIntervals(1,j) )  ) THEN
1932
              Failure = .TRUE.
1933
              EXIT
1934
            END IF
1935
          END DO
1936
          IF ( Failure ) EXIT
1937
        END DO
1938
        IF ( Failure ) EXIT
1939
      END DO
1940
!------------------------------------------------------------------------------
1941
    END FUNCTION CheckLatentHeat
1942
!------------------------------------------------------------------------------
1943

1944

1945

1946
!------------------------------------------------------------------------------
1947
   SUBROUTINE IntegOverA( BoundaryMatrix, BoundaryVector, &
1948
     LOAD, NodalAlpha, Element, n, m, Nodes )
1949
!------------------------------------------------------------------------------
1950
     REAL(KIND=dp) :: BoundaryMatrix(:,:),BoundaryVector(:), &
1951
                    LOAD(:),NodalAlpha(:)
1952

1953
     TYPE(Nodes_t)   :: Nodes
1954
     TYPE(Element_t) :: Element
1955

1956
     INTEGER :: n,  m
1957

1958
     REAL(KIND=dp) :: Basis(n)
1959
     REAL(KIND=dp) :: dBasisdx(n,3),SqrtElementMetric
1960

1961
     REAL(KIND=dp) :: u,v,w,s,x,y,z
1962
     REAL(KIND=dp) :: Force,Alpha
1963
     REAL(KIND=dp), POINTER :: U_Integ(:),V_Integ(:),W_Integ(:),S_Integ(:)
1964

1965
     INTEGER :: i,t,q,p,N_Integ
1966

1967
     TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff
1968

1969
     LOGICAL :: stat
1970
!------------------------------------------------------------------------------
1971

1972
     BoundaryVector = 0.0D0
1973
     BoundaryMatrix = 0.0D0
1974
!------------------------------------------------------------------------------
1975
!    Integration stuff
1976
!------------------------------------------------------------------------------
1977
     IntegStuff = GaussPoints( Element )
1978
     U_Integ => IntegStuff % u
1979
     V_Integ => IntegStuff % v
1980
     W_Integ => IntegStuff % w
1981
     S_Integ => IntegStuff % s
1982
     N_Integ =  IntegStuff % n
1983

1984
!------------------------------------------------------------------------------
1985
!   Now we start integrating
1986
!------------------------------------------------------------------------------
1987
     DO t=1,N_Integ
1988
       u = U_Integ(t)
1989
       v = V_Integ(t)
1990
       w = W_Integ(t)
1991
!------------------------------------------------------------------------------
1992
!     Basis function values & derivatives at the integration point
1993
!------------------------------------------------------------------------------
1994
       stat = ElementInfo( Element,Nodes,u,v,w,SqrtElementMetric, &
1995
                  Basis,dBasisdx )
1996

1997
       s = SqrtElementMetric * S_Integ(t)
1998
!------------------------------------------------------------------------------
1999
!      Coordinatesystem dependent info
2000
!------------------------------------------------------------------------------
2001
       IF ( CurrentCoordinateSystem() /= Cartesian ) THEN
2002
         x = SUM( Nodes % x(1:n)*Basis )
2003
         y = SUM( Nodes % y(1:n)*Basis )
2004
         z = SUM( Nodes % z(1:n)*Basis )
2005
         s = s * CoordinateSqrtMetric( x,y,z )
2006
       END IF
2007
!------------------------------------------------------------------------------
2008
       Force = SUM( LOAD(1:n) * Basis )
2009
       Alpha = SUM( NodalAlpha(1:n) * Basis )
2010

2011
       DO p=1,N
2012
         DO q=1,M
2013
           BoundaryMatrix(p,q) = BoundaryMatrix(p,q) + &
2014
                  s * Alpha * Basis(p) / m
2015
         END DO
2016
       END DO
2017

2018
       DO p=1,N
2019
         BoundaryVector(p) = BoundaryVector(p) + s * Force * Basis(p)
2020
       END DO
2021
     END DO
2022
   END SUBROUTINE IntegOverA
2023
!------------------------------------------------------------------------------
2024

2025

2026

2027
!------------------------------------------------------------------------------
2028
    SUBROUTINE FindGapIndexes( Element, Indexes, n )
2029
!------------------------------------------------------------------------------
2030
      TYPE(Element_t) :: Element
2031
      INTEGER :: n,Indexes(:)
2032
!------------------------------------------------------------------------------
2033
      TYPE(Element_t), POINTER :: Parent,Left,Right
2034
      INTEGER :: i,j,k,l
2035
      REAL(KIND=dp) :: x0,y0,z0,x,y,z
2036
!------------------------------------------------------------------------------
2037
      Left  => Element % BoundaryInfo % Left
2038
      Right => Element % BoundaryInfo % Right
2039

2040
      IF ( .NOT.ASSOCIATED(Left) .OR. .NOT.ASSOCIATED(Right) ) RETURN
2041

2042
      l = 0
2043
      DO i=1,n
2044
        Parent => Left
2045
        k = Element % NodeIndexes(i)
2046

2047
        IF ( ANY( Parent % NodeIndexes == k ) ) &
2048
          Parent => Right
2049

2050
        x0 = ElementNodes % x(i)
2051
        y0 = ElementNodes % y(i)
2052
        z0 = ElementNodes % z(i)
2053
        DO j=1,Parent % TYPE % NumberOfNodes
2054
          k = Parent % NodeIndexes(j)
2055
          x = Solver % Mesh % Nodes % x(k) - x0
2056
          y = Solver % Mesh % Nodes % y(k) - y0
2057
          z = Solver % Mesh % Nodes % z(k) - z0
2058
          IF ( x**2 + y**2 + z**2 < AEPS ) EXIT
2059
        END DO
2060
        Indexes(i) = k
2061
      END DO
2062
!------------------------------------------------------------------------------
2063
    END SUBROUTINE FindGapIndexes
2064
!------------------------------------------------------------------------------
2065

2066

2067
!------------------------------------------------------------------------------
2068
    SUBROUTINE AddHeatGap( Solver, Element, STIFF, EnthalpyPerm )
2069
!------------------------------------------------------------------------------
2070
      TYPE(Solver_t) :: Solver
2071
      REAL(KIND=dp) :: STIFF(:,:)
2072
      INTEGER :: EnthalpyPerm(:)
2073
      TYPE(Element_t) :: Element
2074
!------------------------------------------------------------------------------
2075
      TYPE(Element_t), POINTER :: Parent,Left,Right
2076
      INTEGER :: i,j,k,l, Ind(n)
2077
      REAL(KIND=dp) :: x0,y0,z0,x,y,z
2078
!------------------------------------------------------------------------------
2079
      CALL FindGapIndexes( Element, Ind, n )
2080
      DO i=1,n
2081
        DO j=1,n
2082
          k = EnthalpyPerm( Element % NodeIndexes(i) )
2083
          l = EnthalpyPerm( Ind(j) )
2084
          IF ( k > 0 .AND. l > 0 ) THEN
2085
            CALL AddToMatrixElement( Solver % Matrix,k,l,-STIFF(i,j) )
2086
          END IF
2087
        END DO
2088
      END DO
2089
!------------------------------------------------------------------------------
2090
    END SUBROUTINE AddHeatGap
2091
!------------------------------------------------------------------------------
2092

2093

2094
!------------------------------------------------------------------------------
2095
  END SUBROUTINE EnthalpySolver
2096
!------------------------------------------------------------------------------
2097

2098

2099
!------------------------------------------------------------------------------
2100
  SUBROUTINE EnthalpySolver_Boundary_Residual( Model, Edge, Mesh, Quant, Perm,Gnorm, Indicator )
2101
!------------------------------------------------------------------------------
2102
     USE DefUtils
2103
     USE Radiation
2104

2105
     IMPLICIT NONE
2106
!------------------------------------------------------------------------------
2107
     TYPE(Model_t) :: Model
2108
     INTEGER :: Perm(:)
2109
     REAL(KIND=dp) :: Quant(:), Indicator(2), Gnorm
2110
     TYPE( Mesh_t ), POINTER    :: Mesh
2111
     TYPE( Element_t ), POINTER :: Edge
2112
!------------------------------------------------------------------------------
2113

2114
     TYPE(Nodes_t) :: Nodes, EdgeNodes
2115
     TYPE(Element_t), POINTER :: Element, Bndry
2116

2117
     INTEGER :: i,j,k,n,l,t,DIM,Pn,En
2118
     LOGICAL :: stat, Found
2119

2120
     REAL(KIND=dp), POINTER :: Hwrk(:,:,:)
2121

2122
     REAL(KIND=dp) :: SqrtMetric, Metric(3,3), Symb(3,3,3), dSymb(3,3,3,3)
2123

2124
     REAL(KIND=dp), ALLOCATABLE :: NodalConductivity(:), ExtEnthalpy_h(:), &
2125
       TransferCoeff(:), EdgeBasis(:), Basis(:), x(:), y(:), z(:), &
2126
       dBasisdx(:,:), Enthalpy_h(:), Flux(:), NodalEmissivity(:)
2127

2128
     REAL(KIND=dp) :: Conductivity, Emissivity, StefanBoltzmann
2129

2130
     REAL(KIND=dp) :: Grad(3,3), Normal(3), EdgeLength, gx, gy, gz
2131

2132
     REAL(KIND=dp) :: u, v, w, s, detJ
2133

2134
     REAL(KIND=dp) :: Source, Residual, ResidualNorm, Area
2135

2136
     TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff
2137

2138
     LOGICAL :: First = .TRUE., Dirichlet
2139
     SAVE Hwrk, First
2140
!------------------------------------------------------------------------------
2141

2142
!    Initialize:
2143
!    -----------
2144
     IF ( First ) THEN
2145
        First = .FALSE.
2146
        NULLIFY( Hwrk )
2147
     END IF
2148

2149
     Indicator = 0.0d0
2150
     Gnorm     = 0.0d0
2151

2152
     Metric = 0.0d0
2153
     DO i=1,3
2154
        Metric(i,i) = 1.0d0
2155
     END DO
2156

2157
     SELECT CASE( CurrentCoordinateSystem() )
2158
        CASE( AxisSymmetric, CylindricSymmetric )
2159
           dim = 3
2160
        CASE DEFAULT
2161
           dim = CoordinateSystemDimension()
2162
     END SELECT
2163
!    
2164
!    ---------------------------------------------
2165

2166
     Element => Edge % BoundaryInfo % Left
2167

2168
     IF ( .NOT. ASSOCIATED( Element ) ) THEN
2169
        Element => Edge % BoundaryInfo % Right
2170
     ELSE IF ( ANY( Perm( Element % NodeIndexes ) <= 0 ) ) THEN
2171
        Element => Edge % BoundaryInfo % Right
2172
     END IF
2173

2174
     IF ( .NOT. ASSOCIATED( Element ) ) RETURN
2175
     IF ( ANY( Perm( Element % NodeIndexes ) <= 0 ) ) RETURN
2176

2177
     En = Edge % TYPE % NumberOfNodes
2178
     Pn = Element % TYPE % NumberOfNodes
2179

2180
     ALLOCATE( EdgeNodes % x(En), EdgeNodes % y(En), EdgeNodes % z(En) )
2181

2182
     EdgeNodes % x = Mesh % Nodes % x(Edge % NodeIndexes)
2183
     EdgeNodes % y = Mesh % Nodes % y(Edge % NodeIndexes)
2184
     EdgeNodes % z = Mesh % Nodes % z(Edge % NodeIndexes)
2185

2186
     ALLOCATE( Nodes % x(Pn), Nodes % y(Pn), Nodes % z(Pn) )
2187

2188
     Nodes % x = Mesh % Nodes % x(Element % NodeIndexes)
2189
     Nodes % y = Mesh % Nodes % y(Element % NodeIndexes)
2190
     Nodes % z = Mesh % Nodes % z(Element % NodeIndexes)
2191

2192
     ALLOCATE( Enthalpy_h(Pn), Basis(Pn), ExtEnthalpy_h(En), &
2193
        TransferCoeff(En), x(En), y(En), z(En), EdgeBasis(En), &
2194
        dBasisdx(Pn,3), NodalConductivity(En), Flux(En), &
2195
        NodalEmissivity(En) ) 
2196

2197
     DO l = 1,En
2198
       DO k = 1,Pn
2199
          IF ( Edge % NodeIndexes(l) == Element % NodeIndexes(k) ) THEN
2200
             x(l) = Element % TYPE % NodeU(k)
2201
             y(l) = Element % TYPE % NodeV(k)
2202
             z(l) = Element % TYPE % NodeW(k)
2203
             EXIT
2204
          END IF
2205
       END DO
2206
     END DO
2207
!
2208
!    Integrate square of residual over boundary element:
2209
!    ---------------------------------------------------
2210

2211
     Indicator    = 0.0d0
2212
     EdgeLength   = 0.0d0
2213
     ResidualNorm = 0.0d0
2214

2215
     DO j=1,Model % NumberOfBCs
2216
        IF ( Edge % BoundaryInfo % Constraint /= Model % BCs(j) % Tag ) CYCLE
2217

2218
!       IF ( .NOT. ListGetLogical( Model % BCs(j) % Values, &
2219
!                 'Enthalpy Heat Flux BC', Found ) ) CYCLE
2220

2221
!
2222
!       Check if dirichlet BC given:
2223
!       ----------------------------
2224
        s = ListGetConstReal( Model % BCs(j) % Values,'Enthalpy_h',Dirichlet )
2225

2226
!       Get various flux bc options:
2227
!       ----------------------------
2228

2229
!       ...given flux:
2230
!       --------------
2231
        Flux(1:En) = ListGetReal( Model % BCs(j) % Values, &
2232
          'Enthalpy Heat Flux', En, Edge % NodeIndexes, Found )
2233

2234
!       ...convective heat transfer:
2235
!       ----------------------------
2236
        TransferCoeff(1:En) =  ListGetReal( Model % BCs(j) % Values, &
2237
          'Heat Transfer Coefficient', En, Edge % NodeIndexes, Found )
2238

2239
        ExtEnthalpy_h(1:En) = ListGetReal( Model % BCs(j) % Values, &
2240
          'External Enthalpy_h', En, Edge % NodeIndexes, Found )
2241

2242
!       ...black body radiation:
2243
!       ------------------------
2244
        Emissivity      = 0.0d0
2245
        StefanBoltzmann = 0.0d0
2246

2247
        SELECT CASE(ListGetString(Model % BCs(j) % Values,'Radiation',Found))
2248
           !------------------
2249
           CASE( 'idealized' )
2250
           !------------------
2251

2252
              NodalEmissivity(1:En) = ListGetReal( Model % BCs(j) % Values, &
2253
                   'Emissivity', En, Edge % NodeIndexes, Found)
2254
              IF(.NOT. Found) THEN
2255
                 NodalEmissivity(1:En) = GetParentMatProp( 'Emissivity', Edge)
2256
              END IF
2257
              Emissivity = SUM( NodalEmissivity(1:En)) / En
2258

2259
              StefanBoltzMann = &
2260
                    ListGetConstReal( Model % Constants,'Stefan Boltzmann' )
2261

2262
           !---------------------
2263
           CASE( 'diffuse gray' )
2264
           !---------------------
2265

2266
              NodalEmissivity(1:En) = ListGetReal( Model % BCs(j) % Values, &
2267
                   'Emissivity', En, Edge % NodeIndexes, Found)
2268
              IF(.NOT. Found) THEN
2269
                 NodalEmissivity(1:En) = GetParentMatProp( 'Emissivity', Edge)
2270
              END IF
2271
              Emissivity = SUM( NodalEmissivity(1:En)) / En
2272

2273
              StefanBoltzMann = &
2274
                    ListGetConstReal( Model % Constants,'Stefan Boltzmann' )
2275

2276
              ExtEnthalpy_h(1:En) =  ComputeRadiationLoad( Model, &
2277
                      Mesh, Edge, Quant, Perm, Emissivity )
2278
        END SELECT
2279

2280
!       get material parameters:
2281
!       ------------------------
2282
        k = ListGetInteger(Model % Bodies(Element % BodyId) % Values,'Material', &
2283
                    minv=1, maxv=Model % NumberOFMaterials)
2284

2285
        CALL ListGetRealArray( Model % Materials(k) % Values, &
2286
               'Enthalpy Heat Diffusivity', Hwrk, En, Edge % NodeIndexes )
2287

2288
        NodalConductivity( 1:En ) = Hwrk( 1,1,1:En )
2289

2290
!       elementwise nodal solution:
2291
!       ---------------------------
2292
        Enthalpy_h(1:Pn) = Quant( Perm(Element % NodeIndexes) )
2293

2294
!       do the integration:
2295
!       -------------------
2296
        EdgeLength   = 0.0d0
2297
        ResidualNorm = 0.0d0
2298

2299
        IntegStuff = GaussPoints( Edge )
2300

2301
        DO t=1,IntegStuff % n
2302
           u = IntegStuff % u(t)
2303
           v = IntegStuff % v(t)
2304
           w = IntegStuff % w(t)
2305

2306
           stat = ElementInfo( Edge, EdgeNodes, u, v, w, detJ, &
2307
               EdgeBasis, dBasisdx )
2308

2309
           Normal = NormalVector( Edge, EdgeNodes, u, v, .TRUE. )
2310

2311
           IF ( CurrentCoordinateSystem() == Cartesian ) THEN
2312
              s = IntegStuff % s(t) * detJ
2313
           ELSE
2314
              gx = SUM( EdgeBasis(1:En) * EdgeNodes % x(1:En) )
2315
              gy = SUM( EdgeBasis(1:En) * EdgeNodes % y(1:En) )
2316
              gz = SUM( EdgeBasis(1:En) * EdgeNodes % z(1:En) )
2317
      
2318
              CALL CoordinateSystemInfo( Metric, SqrtMetric, &
2319
                         Symb, dSymb, gx, gy, gz )
2320

2321
              s = IntegStuff % s(t) * detJ * SqrtMetric
2322
           END IF
2323

2324
!
2325
!          Integration point in parent element local
2326
!          coordinates:
2327
!          -----------------------------------------
2328
           u = SUM( EdgeBasis(1:En) * x(1:En) )
2329
           v = SUM( EdgeBasis(1:En) * y(1:En) )
2330
           w = SUM( EdgeBasis(1:En) * z(1:En) )
2331

2332
           stat = ElementInfo( Element, Nodes, u, v, w, detJ, &
2333
                 Basis, dBasisdx )
2334
!
2335
!          Heat conductivity at the integration point:
2336
!          --------------------------------------------
2337
           Conductivity = SUM( NodalConductivity(1:En) * EdgeBasis(1:En) )
2338
!
2339
!          given flux at integration point:
2340
!          --------------------------------
2341
           Residual = -SUM( Flux(1:En) * EdgeBasis(1:En) )
2342

2343
!          convective ...:
2344
!          ----------------
2345
           Residual = Residual + SUM(TransferCoeff(1:En) * EdgeBasis(1:En)) * &
2346
                     ( SUM( Enthalpy_h(1:Pn) * Basis(1:Pn) ) - &
2347
                       SUM( ExtEnthalpy_h(1:En) * EdgeBasis(1:En) ) )
2348

2349
!          black body radiation...:
2350
!          -------------------------
2351
           Residual = Residual + &
2352
                Emissivity * StefanBoltzmann * &
2353
                     ( SUM( Enthalpy_h(1:Pn) * Basis(1:Pn) ) ** 4 - &
2354
                       SUM( ExtEnthalpy_h(1:En) * EdgeBasis(1:En) ) ** 4 )
2355

2356
!          flux given by the computed solution, and 
2357
!          force norm for scaling the residual:
2358
!          -----------------------------------------
2359
           IF ( CurrentCoordinateSystem() == Cartesian ) THEN
2360
              DO k=1,DIM
2361
                 Residual = Residual + Conductivity  * &
2362
                    SUM( dBasisdx(1:Pn,k) * Enthalpy_h(1:Pn) ) * Normal(k)
2363

2364
                 Gnorm = Gnorm + s * (Conductivity * &
2365
                       SUM(dBasisdx(1:Pn,k) * Enthalpy_h(1:Pn)) * Normal(k))**2
2366
              END DO
2367
           ELSE
2368
              DO k=1,DIM
2369
                 DO l=1,DIM
2370
                    Residual = Residual + Metric(k,l) * Conductivity  * &
2371
                       SUM( dBasisdx(1:Pn,k) * Enthalpy_h(1:Pn) ) * Normal(l)
2372

2373
                    Gnorm = Gnorm + s * (Metric(k,l) * Conductivity * &
2374
                      SUM(dBasisdx(1:Pn,k) * Enthalpy_h(1:Pn) ) * Normal(l))**2
2375
                 END DO
2376
              END DO
2377
           END IF
2378

2379
           EdgeLength   = EdgeLength + s
2380
           IF ( .NOT. Dirichlet ) THEN
2381
              ResidualNorm = ResidualNorm + s * Residual ** 2
2382
           END IF
2383
        END DO
2384
        EXIT
2385
     END DO
2386

2387
     IF ( CoordinateSystemDimension() == 3 ) THEN
2388
        EdgeLength = SQRT(EdgeLength)
2389
     END IF
2390

2391
!    Gnorm = EdgeLength * Gnorm
2392
     Indicator = EdgeLength * ResidualNorm
2393

2394
     DEALLOCATE( Nodes % x, Nodes % y, Nodes % z)
2395
     DEALLOCATE( EdgeNodes % x, EdgeNodes % y, EdgeNodes % z)
2396

2397
     DEALLOCATE( Enthalpy_h, Basis, ExtEnthalpy_h, TransferCoeff,  &
2398
        x, y, z, EdgeBasis, dBasisdx, NodalConductivity, Flux, &
2399
        NodalEmissivity ) 
2400
!------------------------------------------------------------------------------
2401
  END SUBROUTINE EnthalpySolver_Boundary_Residual
2402
!------------------------------------------------------------------------------
2403

2404

2405

2406
!------------------------------------------------------------------------------
2407
  SUBROUTINE EnthalpySolver_Edge_Residual( Model, Edge, Mesh, Quant, Perm, Indicator )
2408
!------------------------------------------------------------------------------
2409
     USE DefUtils
2410
     IMPLICIT NONE
2411

2412
     TYPE(Model_t) :: Model
2413
     INTEGER :: Perm(:)
2414
     REAL(KIND=dp) :: Quant(:), Indicator(2)
2415
     TYPE( Mesh_t ), POINTER    :: Mesh
2416
     TYPE( Element_t ), POINTER :: Edge
2417
!------------------------------------------------------------------------------
2418

2419
     TYPE(Nodes_t) :: Nodes, EdgeNodes
2420
     TYPE(Element_t), POINTER :: Element, Bndry
2421

2422
     INTEGER :: i,j,k,l,n,t,DIM,En,Pn
2423
     LOGICAL :: stat, Found
2424
     REAL(KIND=dp), POINTER :: Hwrk(:,:,:)
2425

2426
     REAL(KIND=dp) :: SqrtMetric, Metric(3,3), Symb(3,3,3), dSymb(3,3,3,3)
2427

2428
     REAL(KIND=dp), ALLOCATABLE :: NodalConductivity(:), x(:), y(:), z(:), &
2429
            EdgeBasis(:), Basis(:), dBasisdx(:,:), Enthalpy_h(:)
2430

2431
     REAL(KIND=dp) :: Grad(3,3), Normal(3), EdgeLength, Jump, Conductivity
2432

2433
     REAL(KIND=dp) :: u, v, w, s, detJ
2434

2435
     REAL(KIND=dp) :: Residual, ResidualNorm, Area
2436

2437
     TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff
2438

2439
     LOGICAL :: First = .TRUE.
2440
     SAVE Hwrk, First
2441
!------------------------------------------------------------------------------
2442

2443
     !    Initialize:
2444
     !    -----------
2445
     IF ( First ) THEN
2446
        First = .FALSE.
2447
        NULLIFY( Hwrk )
2448
     END IF
2449

2450
     SELECT CASE( CurrentCoordinateSystem() )
2451
        CASE( AxisSymmetric, CylindricSymmetric )
2452
           DIM = 3
2453
        CASE DEFAULT
2454
           DIM = CoordinateSystemDimension()
2455
     END SELECT
2456

2457
     Metric = 0.0d0
2458
     DO i = 1,3
2459
        Metric(i,i) = 1.0d0
2460
     END DO
2461

2462
     Grad = 0.0d0
2463
!
2464
!    ---------------------------------------------
2465

2466
     Element => Edge % BoundaryInfo % Left
2467
     n = Element % TYPE % NumberOfNodes
2468

2469
     Element => Edge % BoundaryInfo % Right
2470
     n = MAX( n, Element % TYPE % NumberOfNodes )
2471

2472
     ALLOCATE( Nodes % x(n), Nodes % y(n), Nodes % z(n) )
2473

2474
     En = Edge % TYPE % NumberOfNodes
2475
     ALLOCATE( EdgeNodes % x(En), EdgeNodes % y(En), EdgeNodes % z(En) )
2476

2477
     EdgeNodes % x = Mesh % Nodes % x(Edge % NodeIndexes)
2478
     EdgeNodes % y = Mesh % Nodes % y(Edge % NodeIndexes)
2479
     EdgeNodes % z = Mesh % Nodes % z(Edge % NodeIndexes)
2480

2481
     ALLOCATE( NodalConductivity(En), EdgeBasis(En), Basis(n), &
2482
        dBasisdx(n,3), x(En), y(En), z(En), Enthalpy_h(n) )
2483

2484
!    Integrate square of jump over edge:
2485
!    -----------------------------------
2486
     ResidualNorm = 0.0d0
2487
     EdgeLength   = 0.0d0
2488
     Indicator    = 0.0d0
2489

2490
     IntegStuff = GaussPoints( Edge )
2491

2492
     DO t=1,IntegStuff % n
2493

2494
        u = IntegStuff % u(t)
2495
        v = IntegStuff % v(t)
2496
        w = IntegStuff % w(t)
2497

2498
        stat = ElementInfo( Edge, EdgeNodes, u, v, w, detJ, &
2499
             EdgeBasis, dBasisdx )
2500

2501
        Normal = NormalVector( Edge, EdgeNodes, u, v, .FALSE. )
2502

2503
        IF ( CurrentCoordinateSystem() == Cartesian ) THEN
2504
           s = IntegStuff % s(t) * detJ
2505
        ELSE
2506
           u = SUM( EdgeBasis(1:En) * EdgeNodes % x(1:En) )
2507
           v = SUM( EdgeBasis(1:En) * EdgeNodes % y(1:En) )
2508
           w = SUM( EdgeBasis(1:En) * EdgeNodes % z(1:En) )
2509

2510
           CALL CoordinateSystemInfo( Metric, SqrtMetric, &
2511
                      Symb, dSymb, u, v, w )
2512
           s = IntegStuff % s(t) * detJ * SqrtMetric
2513
        END IF
2514

2515
        ! 
2516
        ! Compute flux over the edge as seen by elements
2517
        ! on both sides of the edge:
2518
        ! ----------------------------------------------
2519
        DO i = 1,2
2520
           SELECT CASE(i)
2521
              CASE(1)
2522
                 Element => Edge % BoundaryInfo % Left
2523
              CASE(2)
2524
                 Element => Edge % BoundaryInfo % Right
2525
           END SELECT
2526
!
2527
!          Can this really happen (maybe it can...)  ?      
2528
!          -------------------------------------------
2529
           IF ( ANY( Perm( Element % NodeIndexes ) <= 0 ) ) CYCLE
2530
!
2531
!          Next, get the integration point in parent
2532
!          local coordinates:
2533
!          -----------------------------------------
2534
           Pn = Element % TYPE % NumberOfNodes
2535

2536
           DO j = 1,En
2537
              DO k = 1,Pn
2538
                 IF ( Edge % NodeIndexes(j) == Element % NodeIndexes(k) ) THEN
2539
                    x(j) = Element % TYPE % NodeU(k)
2540
                    y(j) = Element % TYPE % NodeV(k)
2541
                    z(j) = Element % TYPE % NodeW(k)
2542
                    EXIT
2543
                 END IF
2544
              END DO
2545
           END DO
2546

2547
           u = SUM( EdgeBasis(1:En) * x(1:En) )
2548
           v = SUM( EdgeBasis(1:En) * y(1:En) )
2549
           w = SUM( EdgeBasis(1:En) * z(1:En) )
2550
!
2551
!          Get parent element basis & derivatives at the integration point:
2552
!          -----------------------------------------------------------------
2553
           Nodes % x(1:Pn) = Mesh % Nodes % x(Element % NodeIndexes)
2554
           Nodes % y(1:Pn) = Mesh % Nodes % y(Element % NodeIndexes)
2555
           Nodes % z(1:Pn) = Mesh % Nodes % z(Element % NodeIndexes)
2556

2557
           stat = ElementInfo( Element, Nodes, u, v, w, detJ, &
2558
             Basis, dBasisdx )
2559
!
2560
!          Material parameters:
2561
!          --------------------
2562
           k = ListGetInteger( Model % Bodies( &
2563
                    Element % BodyId) % Values, 'Material', &
2564
                     minv=1, maxv=Model % NumberOFMaterials )
2565

2566
           CALL ListGetRealArray( Model % Materials(k) % Values, &
2567
                   'Enthalpy Heat Diffusivity', Hwrk,En, Edge % NodeIndexes )
2568

2569
           NodalConductivity( 1:En ) = Hwrk( 1,1,1:En )
2570
           Conductivity = SUM( NodalConductivity(1:En) * EdgeBasis(1:En) )
2571
!
2572
!          Enthalpy_h at element nodal points:
2573
!          ------------------------------------
2574
           Enthalpy_h(1:Pn) = Quant( Perm(Element % NodeIndexes) )
2575
!
2576
!          Finally, the flux:
2577
!          ------------------
2578
           DO j=1,DIM
2579
              Grad(j,i) = Conductivity * SUM( dBasisdx(1:Pn,j) * Enthalpy_h(1:Pn) )
2580
           END DO
2581
        END DO
2582

2583
!       Compute squre of the flux jump:
2584
!       -------------------------------   
2585
        EdgeLength  = EdgeLength + s
2586
        Jump = 0.0d0
2587
        DO k=1,DIM
2588
           IF ( CurrentCoordinateSystem() == Cartesian ) THEN
2589
              Jump = Jump + (Grad(k,1) - Grad(k,2)) * Normal(k)
2590
           ELSE
2591
              DO l=1,DIM
2592
                 Jump = Jump + &
2593
                       Metric(k,l) * (Grad(k,1) - Grad(k,2)) * Normal(l)
2594
              END DO
2595
           END IF
2596
        END DO
2597
        ResidualNorm = ResidualNorm + s * Jump ** 2
2598
     END DO
2599

2600
     IF ( CoordinateSystemDimension() == 3 ) THEN
2601
        EdgeLength = SQRT(EdgeLength)
2602
     END IF
2603
     Indicator = EdgeLength * ResidualNorm
2604

2605
     DEALLOCATE( Nodes % x, Nodes % y, Nodes % z, x, y, z)
2606
     DEALLOCATE( EdgeNodes % x, EdgeNodes % y, EdgeNodes % z)
2607
     DEALLOCATE( NodalConductivity, EdgeBasis, Basis, dBasisdx, Enthalpy_h)
2608

2609
!------------------------------------------------------------------------------
2610
  END SUBROUTINE EnthalpySolver_Edge_Residual
2611
!------------------------------------------------------------------------------
2612

2613

2614
!------------------------------------------------------------------------------
2615
   SUBROUTINE EnthalpySolver_Inside_Residual( Model, Element, Mesh, &
2616
        Quant, Perm, Fnorm, Indicator )
2617
!------------------------------------------------------------------------------
2618
     USE DefUtils
2619
!------------------------------------------------------------------------------
2620
     IMPLICIT NONE
2621
!------------------------------------------------------------------------------
2622
     TYPE(Model_t) :: Model
2623
     INTEGER :: Perm(:)
2624
     REAL(KIND=dp) :: Quant(:), Indicator(2), Fnorm
2625
     TYPE( Mesh_t ), POINTER    :: Mesh
2626
     TYPE( Element_t ), POINTER :: Element
2627
!------------------------------------------------------------------------------
2628

2629
     TYPE(Nodes_t) :: Nodes
2630

2631
     INTEGER :: i,j,k,l,n,t,DIM
2632

2633
     LOGICAL :: stat, Found, Compressible
2634
     TYPE( Variable_t ), POINTER :: Var
2635

2636
     REAL(KIND=dp), POINTER :: Hwrk(:,:,:)
2637

2638
     REAL(KIND=dp) :: SqrtMetric, Metric(3,3), Symb(3,3,3), dSymb(3,3,3,3)
2639

2640
     REAL(KIND=dp), ALLOCATABLE :: NodalDensity(:)
2641
     REAL(KIND=dp), ALLOCATABLE :: NodalCapacity(:)
2642
     REAL(KIND=dp), ALLOCATABLE :: x(:), y(:), z(:)
2643
     REAL(KIND=dp), ALLOCATABLE :: NodalConductivity(:)
2644
     REAL(KIND=dp), ALLOCATABLE :: Velo(:,:), Pressure(:)
2645
     REAL(KIND=dp), ALLOCATABLE :: NodalSource(:), Enthalpy_h(:), PrevTemp(:)
2646
     REAL(KIND=dp), ALLOCATABLE :: Basis(:), dBasisdx(:,:), ddBasisddx(:,:,:)
2647

2648
     REAL(KIND=dp) :: u, v, w, s, detJ, Density, Capacity
2649

2650
     REAL(KIND=dp) :: SpecificHeatRatio, ReferencePressure, dt
2651
     REAL(KIND=dp) :: Source, Residual, ResidualNorm, Area, Conductivity
2652

2653
     TYPE( ValueList_t ), POINTER :: Material
2654

2655
     TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff
2656

2657
     LOGICAL :: First = .TRUE.
2658
     SAVE Hwrk, First
2659
!------------------------------------------------------------------------------
2660

2661
!    Initialize:
2662
!    -----------
2663
     Indicator = 0.0d0
2664
     Fnorm     = 0.0d0
2665
!
2666
!    Check if this eq. computed in this element:
2667
!    -------------------------------------------
2668
     IF ( ANY( Perm( Element % NodeIndexes ) <= 0 ) ) RETURN
2669

2670
     IF ( First ) THEN
2671
        First = .FALSE.
2672
        NULLIFY( Hwrk )
2673
     END IF
2674

2675
     Metric = 0.0d0
2676
     DO i=1,3
2677
        Metric(i,i) = 1.0d0
2678
     END DO
2679

2680
     SELECT CASE( CurrentCoordinateSystem() )
2681
        CASE( AxisSymmetric, CylindricSymmetric )
2682
           DIM = 3
2683
        CASE DEFAULT
2684
           DIM = CoordinateSystemDimension()
2685
     END SELECT
2686
!
2687
!    Element nodal points:
2688
!    ---------------------
2689
     n = Element % TYPE % NumberOfNodes
2690

2691
     ALLOCATE( Nodes % x(n), Nodes % y(n), Nodes % z(n) )
2692

2693
     Nodes % x = Mesh % Nodes % x(Element % NodeIndexes)
2694
     Nodes % y = Mesh % Nodes % y(Element % NodeIndexes)
2695
     Nodes % z = Mesh % Nodes % z(Element % NodeIndexes)
2696

2697
     ALLOCATE( NodalDensity(n), NodalCapacity(n), NodalConductivity(n),       &
2698
         Velo(3,n), Pressure(n), NodalSource(n), Enthalpy_h(n), PrevTemp(n), &
2699
         Basis(n), dBasisdx(n,3), ddBasisddx(n,3,3) )
2700
!
2701
!    Elementwise nodal solution:
2702
!    ---------------------------
2703
     Enthalpy_h(1:n) = Quant( Perm(Element % NodeIndexes) )
2704
!
2705
!    Check for time dep.
2706
!    -------------------
2707
     PrevTemp(1:n) = Enthalpy_h(1:n)
2708
     dt = Model % Solver % dt
2709
     IF ( ListGetString( Model % Simulation, 'Simulation Type') == 'transient' ) THEN
2710
        Var => VariableGet( Model % Variables, 'Enthalpy_h', .TRUE. )
2711
        PrevTemp(1:n) = Var % PrevValues(Var % Perm(Element % NodeIndexes),1)
2712
     END IF
2713
!
2714
!    Material parameters: conductivity, heat capacity and density
2715
!    -------------------------------------------------------------
2716
     k = ListGetInteger( Model % Bodies(Element % BodyId) % Values, 'Material', &
2717
                     minv=1, maxv=Model % NumberOfMaterials )
2718

2719
     Material => Model % Materials(k) % Values
2720

2721
     CALL ListGetRealArray( Material, &
2722
                  'Enthalpy Heat Diffusivity', Hwrk,n, Element % NodeIndexes )
2723

2724
     NodalConductivity( 1:n ) = Hwrk( 1,1,1:n )
2725

2726
     NodalDensity(1:n) = ListGetReal( Material, &
2727
            'Enthalpy Density', n, Element % NodeIndexes, Found )
2728

2729
     NodalCapacity(1:n) = 1.0 
2730
!
2731
!    Check for compressible flow equations:
2732
!    --------------------------------------
2733
     Compressible = .FALSE.
2734

2735
     IF (  ListGetString( Material, 'Compressibility Model', Found ) == &
2736
                 'perfect gas equation 1' ) THEN
2737

2738
        Compressible = .TRUE.
2739

2740
        Pressure = 0.0d0
2741
        Var => VariableGet( Mesh % Variables, 'Pressure', .TRUE. )
2742
        IF ( ASSOCIATED( Var ) ) THEN
2743
           Pressure(1:n) = &
2744
               Var % Values( Var % Perm(Element % NodeIndexes) )
2745
        END IF
2746

2747
        ReferencePressure = ListGetConstReal( Material, &
2748
                   'Reference Pressure' )
2749

2750
        SpecificHeatRatio = ListGetConstReal( Material, &
2751
                   'Specific Heat Ratio' )
2752

2753
        NodalDensity(1:n) =  (Pressure(1:n) + ReferencePressure) * SpecificHeatRatio / &
2754
              ( (SpecificHeatRatio - 1) * NodalCapacity(1:n) * Enthalpy_h(1:n) )
2755
     END IF
2756
!
2757
!    Get (possible) convection velocity at the nodes of the element:
2758
!    ----------------------------------------------------------------
2759
     k = ListGetInteger( Model % Bodies(Element % BodyId) % Values, 'Equation', &
2760
                minv=1, maxv=Model % NumberOFEquations )
2761

2762
     Velo = 0.0d0
2763
     SELECT CASE( ListGetString( Model % Equations(k) % Values, &
2764
                         'Convection', Found ) )
2765

2766
        !-----------------
2767
        CASE( 'constant' )
2768
        !-----------------
2769

2770
           Velo(1,1:n) = ListGetReal( Material, &
2771
              'Convection Velocity 1', n, Element % NodeIndexes, Found )
2772

2773
           Velo(2,1:n) = ListGetReal( Material, &
2774
              'Convection Velocity 2', n, Element % NodeIndexes, Found )
2775

2776
           Velo(3,1:n) = ListGetReal( Material, &
2777
              'Convection Velocity 3', n, Element % NodeIndexes, Found )
2778

2779
        !-----------------
2780
        CASE( 'computed' )
2781
        !-----------------
2782

2783
           Var => VariableGet( Mesh % Variables, 'Velocity 1', .TRUE. )
2784
           IF ( ASSOCIATED( Var ) ) THEN
2785
              IF ( ALL( Var % Perm( Element % NodeIndexes ) > 0 ) ) THEN
2786
                 Velo(1,1:n) = Var % Values(Var % Perm(Element % NodeIndexes))
2787
   
2788
                 Var => VariableGet( Mesh % Variables, 'Velocity 2', .TRUE. )
2789
                 IF ( ASSOCIATED( Var ) ) &
2790
                    Velo(2,1:n) = Var % Values( &
2791
                              Var % Perm(Element % NodeIndexes ) )
2792
   
2793
                 Var => VariableGet( Mesh % Variables, 'Velocity 3', .TRUE. )
2794
                 IF ( ASSOCIATED( Var ) ) &
2795
                    Velo(3,1:n) = Var % Values( &
2796
                             Var % Perm( Element % NodeIndexes ) )
2797
              END IF
2798
           END IF
2799

2800
     END SELECT
2801

2802
!
2803
!    Heat source:
2804
!    ------------
2805
!
2806
     k = ListGetInteger( &
2807
         Model % Bodies(Element % BodyId) % Values,'Body Force',Found, &
2808
                 1, Model % NumberOFBodyForces)
2809

2810
     NodalSource = 0.0d0
2811
     IF ( Found .AND. k > 0  ) THEN
2812
        NodalSource(1:n) = ListGetReal( Model % BodyForces(k) % Values, &
2813
               'Heat Source', n, Element % NodeIndexes, Found )
2814
     END IF
2815

2816
!
2817
!    Integrate square of residual over element:
2818
!    ------------------------------------------
2819

2820
     ResidualNorm = 0.0d0
2821
     Area = 0.0d0
2822

2823
     IntegStuff = GaussPoints( Element )
2824

2825
     DO t=1,IntegStuff % n
2826
        u = IntegStuff % u(t)
2827
        v = IntegStuff % v(t)
2828
        w = IntegStuff % w(t)
2829

2830
        stat = ElementInfo( Element, Nodes, u, v, w, detJ, &
2831
            Basis, dBasisdx, ddBasisddx, .TRUE., .FALSE. )
2832

2833
        IF ( CurrentCoordinateSystem() == Cartesian ) THEN
2834
           s = IntegStuff % s(t) * detJ
2835
        ELSE
2836
           u = SUM( Basis(1:n) * Nodes % x(1:n) )
2837
           v = SUM( Basis(1:n) * Nodes % y(1:n) )
2838
           w = SUM( Basis(1:n) * Nodes % z(1:n) )
2839

2840
           CALL CoordinateSystemInfo( Metric, SqrtMetric, &
2841
                       Symb, dSymb, u, v, w )
2842
           s = IntegStuff % s(t) * detJ * SqrtMetric
2843
        END IF
2844

2845
        Capacity     = SUM( NodalCapacity(1:n) * Basis(1:n) )
2846
        Density      = SUM( NodalDensity(1:n) * Basis(1:n) )
2847
        Conductivity = SUM( NodalConductivity(1:n) * Basis(1:n) )
2848
!
2849
!       Residual of the convection-diffusion (heat) equation:
2850
!        R = \rho * c_p * (@T/@t + u.grad(T)) - &
2851
!            div(C grad(T)) + p div(u) - h,
2852
!       ---------------------------------------------------
2853
!
2854
!       or more generally:
2855
!
2856
!        R = \rho * c_p * (@T/@t + u^j T_{,j}) - &
2857
!          g^{jk} (C T_{,j}}_{,k} + p div(u) - h
2858
!       ---------------------------------------------------
2859
!
2860
        Residual = -Density * SUM( NodalSource(1:n) * Basis(1:n) )
2861

2862
        IF ( CurrentCoordinateSystem() == Cartesian ) THEN
2863
           DO j=1,DIM
2864
!
2865
!             - grad(C).grad(T):
2866
!             --------------------
2867
!
2868
              Residual = Residual - &
2869
                 SUM( Enthalpy_h(1:n) * dBasisdx(1:n,j) ) * &
2870
                 SUM( NodalConductivity(1:n) * dBasisdx(1:n,j) )
2871

2872
!
2873
!             - C div(grad(T)):
2874
!             -------------------
2875
!
2876
              Residual = Residual - Conductivity * &
2877
                 SUM( Enthalpy_h(1:n) * ddBasisddx(1:n,j,j) )
2878
           END DO
2879
        ELSE
2880
           DO j=1,DIM
2881
              DO k=1,DIM
2882
!
2883
!                - g^{jk} C_{,k}T_{j}:
2884
!                ---------------------
2885
!
2886
                 Residual = Residual - Metric(j,k) * &
2887
                    SUM( Enthalpy_h(1:n) * dBasisdx(1:n,j) ) * &
2888
                    SUM( NodalConductivity(1:n) * dBasisdx(1:n,k) )
2889

2890
!
2891
!                - g^{jk} C T_{,jk}:
2892
!                -------------------
2893
!
2894
                 Residual = Residual - Metric(j,k) * Conductivity * &
2895
                    SUM( Enthalpy_h(1:n) * ddBasisddx(1:n,j,k) )
2896
!
2897
!                + g^{jk} C {_jk^l} T_{,l}:
2898
!                ---------------------------
2899
                 DO l=1,DIM
2900
                    Residual = Residual + Metric(j,k) * Conductivity * &
2901
                      Symb(j,k,l) * SUM( Enthalpy_h(1:n) * dBasisdx(1:n,l) )
2902
                 END DO
2903
              END DO
2904
           END DO
2905
        END IF
2906

2907
!       + \rho * c_p * (@T/@t + u.grad(T)):
2908
!       -----------------------------------
2909
        Residual = Residual + Density * Capacity *  &
2910
           SUM((Enthalpy_h(1:n)-PrevTemp(1:n))*Basis(1:n)) / dt
2911

2912
        DO j=1,DIM
2913
           Residual = Residual + &
2914
              Density * Capacity * SUM( Velo(j,1:n) * Basis(1:n) ) * &
2915
                    SUM( Enthalpy_h(1:n) * dBasisdx(1:n,j) )
2916
        END DO
2917

2918

2919
        IF ( Compressible ) THEN
2920
!
2921
!          + p div(u) or p u^j_{,j}:
2922
!          -------------------------
2923
!
2924
           DO j=1,DIM
2925
              Residual = Residual + &
2926
                 SUM( Pressure(1:n) * Basis(1:n) ) * &
2927
                      SUM( Velo(j,1:n) * dBasisdx(1:n,j) )
2928

2929
              IF ( CurrentCoordinateSystem() /= Cartesian ) THEN
2930
                 DO k=1,DIM
2931
                    Residual = Residual + &
2932
                       SUM( Pressure(1:n) * Basis(1:n) ) * &
2933
                           Symb(j,k,j) * SUM( Velo(k,1:n) * Basis(1:n) )
2934
                 END DO
2935
              END IF
2936
           END DO
2937
        END IF
2938

2939
!
2940
!       Compute also force norm for scaling the residual:
2941
!       -------------------------------------------------
2942
        DO i=1,DIM
2943
           Fnorm = Fnorm + s * ( Density * &
2944
             SUM( NodalSource(1:n) * Basis(1:n) ) ) ** 2
2945
        END DO
2946

2947
        Area = Area + s
2948
        ResidualNorm = ResidualNorm + s *  Residual ** 2
2949
     END DO
2950

2951
!    Fnorm = Element % hk**2 * Fnorm
2952
     Indicator = Element % hK**2 * ResidualNorm
2953

2954
     DEALLOCATE( NodalDensity, NodalCapacity, NodalConductivity,    &
2955
         Velo, Pressure, NodalSource, Enthalpy_h, PrevTemp, Basis, &
2956
         dBasisdx, ddBasisddx )
2957
!------------------------------------------------------------------------------
2958
  END SUBROUTINE EnthalpySolver_Inside_Residual
2959
!------------------------------------------------------------------------------
2960

2961

2962
Product

Resources

Company
