1!/*****************************************************************************/
2! *
3! *  Elmer, A Finite Element Software for Multiphysical Problems
4! *
5! *  Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland
6! *
7! *  This program is free software; you can redistribute it and/or
8! *  modify it under the terms of the GNU General Public License
9! *  as published by the Free Software Foundation; either version 2
10! *  of the License, or (at your option) any later version.
11! *
12! *  This program is distributed in the hope that it will be useful,
13! *  but WITHOUT ANY WARRANTY; without even the implied warranty of
14! *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15! *  GNU General Public License for more details.
16! *
17! *  You should have received a copy of the GNU General Public License
18! *  along with this program (in file fem/GPL-2); if not, write to the
19! *  Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
20! *  Boston, MA 02110-1301, USA.
21! *
22! *****************************************************************************/
23!
24! ******************************************************************************
25! *
26! *  Authors: Juha Ruokolainen
27! *  Email:   Juha.Ruokolainen@csc.fi
28! *  Web:     http://www.csc.fi/elmer
29! *  Address: CSC - IT Center for Science Ltd.
30! *           Keilaranta 14
31! *           02101 Espoo, Finland
32! *
33! *  Original Date: 08 Jun 1997
34! *  Modified by: Jussi Heikonen, Ville Savolainen, Peter Raback
35! *  Modification date: 2.9.2006
36! *
37! *****************************************************************************/
38
39
40!------------------------------------------------------------------------------
41!> Initialization for the primary solver: PhaseChangeSolve
42!> \ingroup Solvers
43!------------------------------------------------------------------------------
44SUBROUTINE PhaseChangeSolve_Init( Model,Solver,dt,TransientSimulation)
45!------------------------------------------------------------------------------
46    USE DefUtils
47    IMPLICIT NONE
48!------------------------------------------------------------------------------
49    TYPE(Model_t)  :: Model
50    TYPE(Solver_t), TARGET :: Solver
51    LOGICAL ::  TransientSimulation
52    REAL(KIND=dp) :: dt
53!------------------------------------------------------------------------------
54    LOGICAL :: Found
55    TYPE(ValueList_t), POINTER :: Params
56    CHARACTER(LEN=MAX_NAME_LEN) :: VariableName
57
58    Params => GetSolverParams()
59    VariableName = GetString(Params,'Variable')
60    CALL ListAddString( Params,NextFreeKeyword('Exported Variable ',Params), &
61          '-nooutput '//TRIM(ComponentName(VariableName))//'Move' )
62
63    IF (ListGetLogical(Params,'Use Average Velocity',Found)) &
64      CALL ListAddString( Params,NextFreeKeyword('Exported Variable ',Params), &
65          '-nooutput '//TRIM(ComponentName(VariableName))//'MoveAve' )
66!------------------------------------------------------------------------------
67END SUBROUTINE PhaseChangeSolve_Init
68!------------------------------------------------------------------------------
69
70
71!------------------------------------------------------------------------------
72!>  Solve the free surface in the phase change problem using Lagrangian techniques.
73!> \deprecated This had been replaced by separate versions for transient and steady state phase change.
74!> \ingroup Solvers
75!------------------------------------------------------------------------------
76SUBROUTINE PhaseChangeSolve( Model,Solver,dt,TransientSimulation )
77!------------------------------------------------------------------------------
78  USE DefUtils
79
80  IMPLICIT NONE
81!------------------------------------------------------------------------------
82  TYPE(Model_t)  :: Model
83  TYPE(Solver_t), TARGET :: Solver
84  LOGICAL ::  TransientSimulation
85  REAL(KIND=dp) :: dt
86!------------------------------------------------------------------------------
87!    Local variables
88!------------------------------------------------------------------------------
89  TYPE(Element_t), POINTER :: CurrentElement, Parent, Element
90  TYPE(Variable_t), POINTER :: SurfSol, TempSol, HelpSol, HelpSol2
91  TYPE(Nodes_t) :: Nodes, PNodes
92  TYPE(GaussIntegrationPoints_t) :: IntegStuff
93  TYPE(Matrix_t),POINTER  :: StiffMatrix
94  TYPE(ValueList_t), POINTER :: Material
95  TYPE(Solver_t), POINTER :: PSolver
96
97  REAL(KIND=dp) :: Normal(3), u, v, w, UPull(3), PrevUpull(3), &
98      Density, Update, MaxUpdate, MaxTempDiff, Relax, LocalRelax, AverageRelax, &
99      AverageRelax2, surf, xx, yy, r, detJ, Temp, MeltPoint, NonlinearTol, &
100      Temp1,Temp2,Temppi,dxmin,dymin, xmin, xmax, d, HeatCond, s, tave, prevtave, cvol, clim, ccum=1, &
101      tabs, prevtabs, volabs, prevvolabs, CoordMin(3), CoordMax(3), RelativeChange, &
102      dTdz, ttemp, NewtonAfterTol, area, volume, prevvolume, Eps, &
103      Norm, maxds, maxds0, ds, FluxCorrect = 1.0, dpos, &
104      pos0=0.0, prevpos0, ThermalInertiaFactor, Coeff, OrigNorm, &
105      SteadyDt, AveragingDt, Trip_Temp, MeanSurface, PullRelax, &
106      SpeedUp, PrevNorm
107  REAL(KIND=dp), POINTER :: Surface(:), Temperature(:), ForceVector(:),  &
108      x(:), y(:), z(:), Basis(:), dBasisdx(:,:), NodalTemp(:), &
109      Conductivity(:), LatentHeat(:), TempDiff(:), Taverage(:), Taverage2(:), &
110      Normals(:), Weights(:), SurfaceMove(:), SurfaceMoveAve(:)
111  REAL (KIND=dp), ALLOCATABLE :: PrevTemp(:), IsoSurf(:,:)
112
113  REAL(KIND=dp), ALLOCATABLE :: &
114      LocalStiffMatrix(:,:), LocalForceVector(:), LocalMassMatrix(:,:)
115
116  INTEGER :: i,j,k,t,n,nn,pn,DIM,kl,kr,l, bc, Trip_node, NoBNodes, NonlinearIter, istat, &
117       NElems,ElementCode,Next,Vertex,ii,imin,NewtonAfterIter,Node, iter, LiquidInd, Visited = -1, &
118       SubroutineVisited = 0, NormalDirection, TangentDirection, CoordMini(3), CoordMaxi(3), &
119       Axis_node
120  INTEGER, POINTER :: NodeIndexes(:),TempPerm(:),SurfPerm(:),NormalsPerm(:)
121
122  LOGICAL :: Stat, FirstTime = .TRUE., Newton = .FALSE., UseTaverage, Debug, &
123      PullControl, PullVelocitySet = .FALSE., IsoSurfAllocated, AllocationsDone = .FALSE., &
124      TransientAlgo = .FALSE., SteadyAlgo = .FALSE., UseLoads, &
125      AverageNormal, SurfaceVelocitySet = .FALSE., TriplePointFixed
126  LOGICAL, POINTER :: NodeDone(:)
127  CHARACTER(LEN=MAX_NAME_LEN) :: VariableName, str
128
129  SAVE FirstTime, Trip_node, NoBNodes, SubroutineVisited, prevpos0, &
130      PrevTemp, Newton, ccum, NormalDirection, TangentDirection, MeltPoint, &
131      ForceVector, FluxCorrect, NodeDone, Eps, PullControl, &
132      Visited, Nodes, PNodes, NodalTemp, Conductivity, LatentHeat, &
133      TempDiff, AllocationsDone, LocalStiffMatrix, LocalForceVector, LocalMassMatrix, &
134      x, y, z, Basis, dBasisdx, norm, Axis_node, PullVelocitySet, &
135      Normals, Weights, NormalsPerm, AverageNormal, SurfaceMove, SurfaceMoveAve, &
136      SurfaceVelocitySet, CoordMax, CoordMin, CoordMaxi, CoordMini, UPull
137
138  !------------------------------------------------------------------------------
139  ! Decide which kind of algorithm to use for the current timestep size
140  !------------------------------------------------------------------------------
141
142  TransientAlgo = TransientSimulation
143  IF(TransientAlgo) THEN
144    SteadyDt = ListGetConstReal(Solver % Values,'Steady Transition Timestep',Stat)
145    IF(Stat) TransientAlgo = (dt < SteadyDt)
146  END IF
147  SteadyAlgo = .NOT. TransientAlgo
148
149  CALL Info('PhaseChangeSolve',                  '--------------------------------------------')
150  IF(SteadyAlgo) CALL Info('PhaseChangeSolve',   'Using steady algorithm to find the isotherm')
151  IF(TransientAlgo) CALL Info('PhaseChangeSolve','Using transient algorithm for surface update')
152  CALL Info('PhaseChangeSolve',                  '--------------------------------------------')
153
154  SubroutineVisited = SubroutineVisited + 1
155
156  DIM = CoordinateSystemDimension()
157  IF(DIM /= 2) THEN
158    CALL Fatal('PhaseChangeSolve','Implemented only in 2D')
159  END IF
160
161!------------------------------------------------------------------------------
162! Get variables needed for solution
163!------------------------------------------------------------------------------
164
165  PSolver => Solver
166  SurfSol  => Solver % Variable
167  Surface  => SurfSol % Values
168  SurfPerm => SurfSol % Perm
169  IF(.NOT. ASSOCIATED (Surface) .OR. ALL(SurfPerm <= 0) ) THEN
170    CALL Fatal('PhaseChangeSolve','Surface field needed for Phase Change')
171  END IF
172
173  VariableName = ListGetString( Solver % Values, 'Phase Change Variable', Stat )
174  IF(Stat) THEN
175    TempSol => VariableGet( Solver % Mesh % Variables, TRIM(VariableName) )
176  ELSE
177    TempSol => VariableGet( Solver % Mesh % Variables, 'Temperature' )
178  END IF
179  TempPerm    => TempSol % Perm
180  Temperature => TempSol % Values
181  IF(.NOT. ASSOCIATED (Temperature) .OR. ALL(TempPerm <= 0) ) THEN
182    CALL Fatal('PhaseChangeSolve','Temperature field needed for Phase Change')
183  END IF
184
185  Relax = GetCReal( Solver % Values,  &
186      'Nonlinear System Relaxation Factor', stat )
187  IF ( .NOT. stat ) Relax = 1.0d0
188  NonlinearIter = ListGetInteger( Solver % Values, &
189      'Nonlinear System Max Iterations', stat )
190  IF ( .NOT. stat ) NonlinearIter = 1
191  NonlinearTol  = ListGetConstReal( Solver % Values, &
192      'Nonlinear System Convergence Tolerance', stat )
193
194  PullControl = ListGetLogical( Solver % Values,'Pull Rate Control',stat)
195  TriplePointFixed = ListGetLogical( Solver % Values,'Triple Point Fixed',stat)
196
197!---------------------------------------------------------------------------------
198! The first time the main axis of the free surface is determined
199! and some permanent vectors related to the surface are allocated.
200!---------------------------------------------------------------------------------
201
202  IF(FirstTime) THEN
203    UPull = 0.0
204    NoBNodes = 0
205    CoordMax = -HUGE(CoordMax)
206    CoordMin = HUGE(CoordMin)
207
208    DO k=1, Model % Mesh % NumberOfNodes
209      IF( SurfPerm(k) <= 0) CYCLE
210      NoBnodes = NoBnodes + 1
211
212      DO j=1,DIM
213        IF(j==1) xx = Model % Mesh % Nodes % x(k)
214        IF(j==2) xx = Model % Mesh % Nodes % y(k)
215        IF(j==3) xx = Model % Mesh % Nodes % z(k)
216        IF(xx > CoordMax(j)) THEN
217          CoordMax(j) = xx
218          CoordMaxi(j) = k
219        END IF
220        IF(xx < CoordMin(j)) THEN
221          CoordMin(j) = xx
222          CoordMini(j) = k
223        END IF
224      END DO
225    END DO
226
227    ! Direction of minimum change
228    j = 1
229    DO i=1,DIM
230      IF(CoordMax(i)-CoordMin(i) < CoordMax(j)-CoordMin(j)) THEN
231        j = i
232      END IF
233    END DO
234    NormalDirection = j
235
236    ! Direction of maximum change
237    j = 1
238    DO i=1,DIM
239      IF(CoordMax(i)-CoordMin(i) > CoordMax(j)-CoordMin(j)) THEN
240        j = i
241      END IF
242    END DO
243    TangentDirection = j
244
245    ! Find triple point:
246    !-------------------
247    Trip_node = CoordMaxi(TangentDirection)
248    Axis_node = CoordMini(TangentDirection)
249
250    n = Solver % Mesh % MaxElementNodes
251
252    ALLOCATE( Nodes % x(n), Nodes % y(n), Nodes % z(n), &
253        PNodes % x(n), PNodes % y(n), PNodes % z(n), &
254        x(n), y(n), z(n), Basis(n), dBasisdx(n,3), NodalTemp(n), &
255        Conductivity(n), LatentHeat(n), TempDiff(n), &
256        LocalStiffMatrix(n,n), LocalForceVector(n), LocalMassMatrix(n,n), &
257        PrevTemp(NobNodes), &
258        STAT=istat)
259    IF ( istat /= 0 ) CALL Fatal( 'PhaseChangeSolver', 'Memory allocation error 1.' )
260
261    Nodes % x = 0.0d0
262    Nodes % y = 0.0d0
263    Nodes % z = 0.0d0
264    PrevTemp = 0.0d0
265
266    IF( SteadyAlgo ) THEN
267      ALLOCATE( NodeDone( NobNodes ), STAT=istat)
268      IF ( istat /= 0 ) CALL Fatal( 'PhaseChangeSolver', 'Memory allocation error 3.' )
269    END IF
270
271    VariableName = ListGetString( Solver % Values, 'Normal Variable', Stat )
272    IF(Stat) THEN
273      HelpSol => VariableGet( Solver % Mesh % Variables, TRIM(VariableName), ThisOnly=.TRUE. )
274    ELSE
275      HelpSol => VariableGet( Solver % Mesh % Variables, 'Normals',ThisOnly=.TRUE. )
276    END IF
277    IF(ASSOCIATED(HelpSol)) THEN
278      Normals => HelpSol % Values
279      NormalsPerm => HelpSol % Perm
280      AverageNormal = .TRUE.
281    ELSE
282      AverageNormal = .FALSE.
283    END IF
284
285    HelpSol => VariableGet( Solver % Mesh % Variables, &
286         TRIM(ComponentName(Solver % Variable))//'Move', stat)
287    IF(ASSOCIATED(HelpSol)) THEN
288      SurfaceMove => HelpSol % Values
289    ELSE
290      ALLOCATE( SurfaceMove(SIZE(Surface)), STAT=istat)
291      IF ( istat /= 0 ) CALL Fatal( 'PhaseChangeSolver', 'Memory allocation error.' )
292      SurfaceMove = 0.0d0
293      CALL VariableAdd( Solver % Mesh % Variables, Solver % Mesh, &
294          PSolver,TRIM(ComponentName(Solver % Variable))//'Move',1, &
295              SurfaceMove,SurfPerm,Output=.FALSE.)
296    END IF
297
298    IF(ListGetLogical(Solver % Values,'Use Average Velocity',Stat)) THEN
299      HelpSol => VariableGet( Solver % Mesh % Variables, &
300             TRIM(ComponentName(Solver % Variable))//'MoveAve', stat)
301      IF(ASSOCIATED(HelpSol)) THEN
302        SurfaceMoveAve => HelpSol % Values
303      ELSE
304        ALLOCATE( SurfaceMoveAve(SIZE(Surface)), STAT=istat)
305        IF ( istat /= 0 ) CALL Fatal( 'PhaseChangeSolver', 'Memory allocation error.' )
306        SurfaceMoveAve = 0.0d0
307        CALL VariableAdd( Solver % Mesh % Variables, Solver % Mesh, &
308            PSolver,TRIM(ComponentName(Solver % Variable))//'MoveAve', &
309               1,SurfaceMoveAve,SurfPerm,Output=.FALSE.)
310      END IF
311    END IF
312
313    AllocationsDone = .TRUE.
314  END IF
315
316  ! Find triple point temperature
317  !-------------------------------
318
319  Trip_Temp =  Temperature( TempPerm(Trip_node) )
320
321
322  i =  ListGetInteger( Solver % Values,'Passive Steps',Stat)
323  IF( i >= SubroutineVisited) GOTO 200
324
325  ! Find the constant material parameters
326  !--------------------------------------
327
328  CurrentElement => Solver % Mesh % Elements(Solver % ActiveElements (1))
329  n = CurrentElement % TYPE % NumberOfNodes
330  NodeIndexes => CurrentElement % NodeIndexes
331  k = ListGetInteger(Model % Bodies(CurrentElement % BodyId) % Values,'Material')
332  Material => Model % Materials(k) % Values
333  MeltPoint = ListGetConstReal( Material,'Melting Point' )
334  Density = ListGetConstReal( Material, 'Density' )
335
336  ! The first velocity should always be set
337  IF ( .NOT. PullVelocitySet ) THEN
338    Upull(1) = ListGetConstReal( Material, 'Convection Velocity 1', Stat )
339    Upull(2) = ListGetConstReal( Material, 'Convection Velocity 2', Stat )
340    Upull(3) = 0.0
341    PullVelocitySet = .TRUE.
342  END IF
343
344!--------------------------------------------------------------------
345! The transient algorithm
346! In the transient algorhitm a heat flux over the interface is computed and
347! it is assumed to be used solely in the melting of the solid into liquid.
348! This melting speed gives an estimate for the melting speed that may be
349! improved by iteration.
350!--------------------------------------------------------------------
351
352  IF ( TransientAlgo ) THEN
353    StiffMatrix => Solver % Matrix
354    ForceVector => Solver % Matrix % RHS
355
356    UseLoads = ListGetLogical( Solver % Values,'Use Heat Load',Stat)
357    IF(UseLoads) THEN
358      HelpSol => VariableGet( Solver % Mesh % Variables, 'Nodal Heat Load' )
359    END IF
360
361    PrevUpull = Upull
362    PrevTemp = SurfaceMove
363
364    ! nonlinear iteration could only be associated if normal is computed from the solution
365    DO iter = 1, NonlinearIter
366
367      ! First solve the velocity field
368      CALL InitializeToZero( StiffMatrix, ForceVector )
369
370      IF(UseLoads) THEN
371        DO t=1,Solver % Mesh % NumberOfNodes
372          i = SurfPerm(t)
373          IF( i <= 0) CYCLE
374          ForceVector(i) = HelpSol % Values(HelpSol % Perm(t))
375        END DO
376      END IF
377
378      DO t = 1, Solver % NumberOfActiveElements
379        CurrentElement => Solver % Mesh % Elements(Solver % ActiveElements (t))
380        n = CurrentElement % TYPE % NumberOfNodes
381        NodeIndexes => CurrentElement % NodeIndexes
382
383        k = ListGetInteger(Model % Bodies(CurrentElement % BodyId) % Values,'Material')
384        Material => Model % Materials(k) % Values
385
386        LatentHeat(1:n) = ListGetReal( Material, 'Latent Heat', n, NodeIndexes )
387        Conductivity(1:n) = ListGetReal( Material,'Heat Conductivity', n, NodeIndexes )
388        TempDiff(1:n) = Temperature( TempPerm(NodeIndexes) ) - MeltPoint
389
390        CALL VelocityLocalMatrix( LocalStiffMatrix, LocalMassMatrix, LocalForceVector,&
391            CurrentElement, n, Nodes )
392
393        CALL UpdateGlobalEquations( StiffMatrix, LocalStiffMatrix, &
394            ForceVector, LocalForceVector, n, 1, SurfPerm(NodeIndexes) )
395      END DO
396
397      CALL FinishAssembly( Solver, ForceVector )
398
399      ! No Dirihtlet conditions here since
400      ! One should not really try to force the phase change at some point,
401      ! rather use feedback to tune the pull velocity
402
403      CALL SolveSystem( StiffMatrix, ParMatrix, ForceVector, SurfaceMove, Norm, 1, Solver )
404      RelativeChange = Solver % Variable % NonlinChange
405
406      WRITE( Message, * ) 'Result Norm     : ',Norm
407      CALL Info( 'PhaseChangeSolve', Message, Level=4 )
408
409      WRITE( Message, * ) 'Relative Change : ',RelativeChange
410      CALL Info( 'PhaseChangeSolve', Message, Level=4 )
411
412      IF ( RelativeChange < NonLinearTol ) EXIT
413    END DO
414
415    IF(ListGetLogical(Solver % Values,'Use Average Velocity',Stat)) THEN
416      IF(SurfaceVelocitySet) THEN
417        LocalRelax = GetCReal(Solver % Values,'Velocity Averaging Factor')
418        SurfaceMoveAve = LocalRelax * SurfaceMove + (1-LocalRelax) * SurfaceMoveAve
419      ELSE
420        SurfaceMoveAve = SurfaceMove
421        SurfaceVelocitySet = .TRUE.
422      END IF
423      LocalRelax = GetCReal(Solver % Values,'Velocity Relaxation Factor',Stat)
424      IF(.NOT. Stat) LocalRelax = 1.0d0
425      SurfaceMove = LocalRelax * SurfaceMove + (1-LocalRelax) * SurfaceMoveAve
426    ELSE
427      LocalRelax = GetCReal(Solver % Values,'Velocity Relaxation Factor',Stat)
428      IF(.NOT. Stat) LocalRelax = 1.0d0
429      SurfaceMove = LocalRelax * SurfaceMove + (1-LocalRelax) * PrevTemp
430    END IF
431
432    IF(PullControl .OR. TriplePointFixed) THEN
433      Upull(NormalDirection) = -SurfaceMove(SurfPerm(Trip_node))
434      IF(PullControl) THEN
435        WRITE(Message,*) 'Pull velocity: ', Upull(NormalDirection)
436        CALL Info('PhaseChangeSolve',Message)
437      END IF
438    END IF
439
440
441    ! Then solve the corresponding update in displacement field
442    SpeedUp = ListGetConstReal( solver % Values,'Transient SpeedUp',Stat)
443    IF(.NOT. Stat) SpeedUp = 1.0d0
444    !-----------------------------------------------------------------------------------------
445    ! If no smoothing is applied directly to the free surface then its update may be computed
446    ! directly from the velocity field.
447    !-----------------------------------------------------------------------------------------
448    IF( ListGetConstReal(Solver % Values,'Surface Smoothing Factor') < 1.0d-20) THEN
449      Surface = Surface + SpeedUp * dt * ( SurfaceMove + Upull(NormalDirection))
450    ELSE
451      CALL InitializeToZero( StiffMatrix, ForceVector )
452
453      DO t = 1, Solver % NumberOfActiveElements
454        CurrentElement => Solver % Mesh % Elements(Solver % ActiveElements (t))
455        n = CurrentElement % TYPE % NumberOfNodes
456        NodeIndexes => CurrentElement % NodeIndexes
457
458        CALL SurfaceLocalMatrix( LocalStiffMatrix, LocalMassMatrix, LocalForceVector,&
459            CurrentElement, n, Nodes )
460
461        CALL Add1stOrderTime( LocalMassMatrix, LocalStiffMatrix, &
462            LocalForceVector, dt, n, 1, SurfPerm(NodeIndexes), Solver )
463
464        CALL UpdateGlobalEquations( StiffMatrix, LocalStiffMatrix, &
465            ForceVector, LocalForceVector, n, 1, SurfPerm(NodeIndexes) )
466      END DO
467
468      CALL FinishAssembly( Solver, ForceVector )
469      CALL SolveSystem( StiffMatrix, ParMatrix, ForceVector, Surface, Norm, 1, Solver )
470    END IF
471
472    IF(PullControl .OR. TriplePointFixed) THEN
473      IF(Visited /= Solver % DoneTime) THEN
474        IF(Solver % DoneTime == 1) THEN
475          prevpos0 = 0.0
476        ELSE
477          prevpos0 = pos0
478        END IF
479        Visited = Solver % DoneTime
480      END IF
481
482      pos0 = prevpos0 + UPull(NormalDirection) * dt
483
484      IF(PullControl) THEN
485        ! This sets the maximum position of the crystal side that is still
486        ! aligned with the pull direction.
487        CALL FindPullBoundary()
488      END IF
489
490    END IF
491  END IF
492
493!--------------------------------------------------------------------
494! Transient algorithm end
495!--------------------------------------------------------------------
496
497
498!----------------------------------------------------------------------------
499! The Steady State the simulation is based on a geometric determination of the
500! isotherm. The solution may be accelerated using local or global Newton
501! type of iteration. It is activated only after the solution is quite accurate.
502!-----------------------------------------------------------------------------
503
504  IF( SteadyAlgo ) THEN
505
506    UseTaverage = ListGetLogical( Solver % Values,'Average Temperature', stat )
507    AveragingDt = ListGetConstReal(Solver % Values,'Steady Averaging Timestep',Stat)
508
509    IF(UseTaverage) THEN
510      HelpSol => VariableGet( Solver % Mesh % Variables, 'Taverage' )
511      IF(.NOT. ASSOCIATED (HelpSol)) THEN
512        ALLOCATE( Taverage( Model % NumberOfNodes ), STAT=istat )
513        IF ( istat /= 0 ) CALL Fatal( 'PhaseChangeSolver', 'Memory allocation error.' )
514        CALL VariableAdd( Solver % Mesh % Variables, Solver % Mesh, &
515            PSolver, 'Taverage', 1, Taverage, TempPerm)
516        HelpSol => VariableGet( Solver % Mesh % Variables, 'Taverage' )
517      END IF
518
519      NULLIFY(HelpSol2)
520      AverageRelax2 = GetCReal(Solver % Values,'Temperature Averaging Factor',stat)
521      IF(stat) THEN
522        HelpSol2 => VariableGet( Solver % Mesh % Variables, 'Taverage Slow' )
523        IF(.NOT. ASSOCIATED (HelpSol2)) THEN
524          ALLOCATE( Taverage2( Model % NumberOfNodes ), STAT=istat )
525          IF ( istat /= 0 ) CALL Fatal( 'PhaseChangeSolver', 'Memory allocation error.' )
526          CALL VariableAdd( Solver % Mesh % Variables, Solver % Mesh, &
527              PSolver, 'Taverage Slow', 1, Taverage2, TempPerm)
528          HelpSol2 => VariableGet( Solver % Mesh % Variables, 'Taverage Slow' )
529        END IF
530      END IF
531
532      IF(dt >= AveragingDt) THEN
533        HelpSol % Values = TempSol % Values
534        IF(ASSOCIATED(HelpSol2)) HelpSol2 % Values = TempSol % Values
535      ELSE
536        AverageRelax = GetCReal( Solver % Values,'Temperature Relaxation Factor')
537        IF(.NOT. ASSOCIATED(HelpSol2)) THEN
538          HelpSol % Values = (1.0d0-AverageRelax) * HelpSol % Values + AverageRelax * TempSol % Values
539        ELSE
540          HelpSol2 % Values = (1.0d0-AverageRelax2) * HelpSol2 % Values + AverageRelax2 * TempSol % Values
541          HelpSol % Values = (1.0d0-AverageRelax) * HelpSol2 % Values + AverageRelax * TempSol % Values
542        END IF
543        PRINT *,'temperature set to average temperature'
544        Temperature => HelpSol % Values
545      END IF
546    END IF
547
548    NewtonAfterIter = ListGetInteger( Solver % Values, &
549        'Nonlinear System Newton After Iterations', stat )
550    IF ( stat .AND. SubroutineVisited > NewtonAfterIter ) Newton = .TRUE.
551
552    NewtonAfterTol = ListGetConstReal( Solver % Values, &
553        'Nonlinear System Newton After Tolerance', stat )
554
555    IF(Newton) THEN
556      CALL Info( 'PhaseChangeSolve','Steady state newton formulation', Level=4 )
557    ELSE
558      CALL Info( 'PhaseChangeSolve','Steady state isoterm formulation', Level=4 )
559    END IF
560
561! Find the melting point:
562!-----------------------
563
564    IF( TriplePointFixed ) THEN
565      MeltPoint = Temperature( TempPerm(Trip_node) )
566      WRITE(Message,*) 'Melting point set to triple point temperature',MeltPoint,Trip_node
567      CALL Info('PhaseChangeSolve',Message)
568    END IF
569
570
571! Create the isotherm for T=T_0:
572!------------------------------
573
574    IF(.NOT. Newton) THEN
575
576      IsoSurfAllocated = .FALSE.
577      xmin = HUGE(xmin)
578      xmax = -HUGE(xmax)
579
580100   NElems = 0
581
582      DO t=1,Solver % Mesh % NumberOfBulkElements
583
584        CurrentElement => Solver % Mesh % Elements(t)
585
586        n = CurrentElement % TYPE % NumberOfNodes
587        NodeIndexes => CurrentElement % NodeIndexes
588        ElementCode = CurrentElement % TYPE % ElementCode
589
590        k = ListGetInteger(Model % Bodies(CurrentElement % BodyId) % Values,'Material')
591        IF (.NOT. (ListGetLogical(Model % Materials(k) % Values, 'Solid', stat) .OR. &
592            ListGetLogical(Model % Materials(k) % Values, 'Liquid', stat) )) CYCLE
593
594        IF(ElementCode < 300) CYCLE
595
596        IF ( ANY( TempPerm( NodeIndexes(1:n) ) <= 0 ) )  CYCLE
597        TempDiff(1:n) = Temperature(TempPerm(NodeIndexes(1:n))) - MeltPoint
598
599        IF( ALL ( TempDiff(1:n) < 0.0 ) ) CYCLE
600        IF( ALL ( TempDiff(1:n) > 0.0 ) ) CYCLE
601
602        Nodes % x(1:n) = Solver % Mesh % Nodes % x(NodeIndexes)
603        Nodes % y(1:n) = Solver % Mesh % Nodes % y(NodeIndexes)
604
605        Vertex = ElementCode/100
606        n=0
607        DO nn=1,Vertex
608          next  = MODULO(nn,Vertex) + 1
609
610          temp1 = TempDiff(nn)
611          temp2 = TempDiff(next)
612
613          IF ( ( (temp1 < 0.0) .AND. (0.0 <= temp2) ) .OR. &
614              ( (temp2 <= 0.0) .AND. (0.0 < temp1) ) ) THEN
615
616            n = n + 1
617
618            IF ( n <= 2 ) THEN
619              NElems = NElems + 1
620              IF(IsoSurfAllocated) THEN
621                IsoSurf(NElems,1) = Nodes % x(nn) + &
622                    temp1 * ((Nodes % x(next) - Nodes % x(nn)) / (temp1-temp2))
623                IsoSurf(NElems,2) = Nodes % y(nn) + &
624                    temp1 * ((Nodes % y(next) - Nodes % y(nn)) / (temp1-temp2))
625
626                xmin = MIN( IsoSurf(Nelems,1), xmin )
627                xmax = MAX( IsoSurf(Nelems,1), xmax )
628              END IF
629            ELSE
630              CALL Warn('PhaseChangeSolve','Wiggly Isotherm')
631              WRITE(Message,*) 'nodeindexes',nodeindexes(1:Vertex)
632              CALL Warn('PhaseChangeSolve',Message)
633              WRITE(Message,*) 'temperature',Temperature(TempPerm(NodeIndexes))
634              CALL Warn('PhaseChangeSolve',Message)
635              PRINT *,'Nodes % x',Nodes % x(n)
636              PRINT *,'Nodes % y',Nodes % y(n)
637            END IF
638          END IF
639        END DO
640
641        IF ( n == 1 ) THEN
642          NElems = NElems - 1
643          CYCLE
644        END IF
645
646        IF (IsoSurfAllocated .AND. n == 2) THEN
647          IF ( IsoSurf(Nelems-1,TangentDirection) > IsoSurf(Nelems,TangentDirection) ) THEN
648            Temppi = IsoSurf(Nelems-1,1)
649            IsoSurf(Nelems-1,1) = IsoSurf(Nelems,1)
650            IsoSurf(Nelems,1) = Temppi
651            Temppi = IsoSurf(Nelems-1,2)
652            IsoSurf(Nelems-1,2) = IsoSurf(Nelems,2)
653            IsoSurf(Nelems,2) = Temppi
654          END IF
655        END IF
656
657      END DO
658
659      IF(Nelems == 0) CALL Fatal('PhaseChangeSolve','Isotherm is empty thus cannot map phase change surface')
660
661      IF(.NOT. IsoSurfAllocated) THEN
662        ALLOCATE( IsoSurf(Nelems+1,2))
663        IsoSurfAllocated = .TRUE.
664        WRITE(Message,*) 'Isotherm created with number of segments',Nelems
665        CALL Info('PhaseChangeSolve',Message)
666        GOTO 100
667      END IF
668    END IF
669
670    area = 0.0
671    volume = 0.0
672    tave = 0.0
673    tabs = 0.0
674    volabs = 0.0
675    NodeDone = .FALSE.
676    MaxTempDiff = 0.0
677
678
679    DO t = 1, Solver % NumberOfActiveElements
680
681      CurrentElement => Solver % Mesh % Elements(Solver % ActiveElements(t))
682      n = CurrentElement % TYPE % NumberOfNodes
683      NodeIndexes => CurrentElement % NodeIndexes
684
685      ElementCode = CurrentElement % TYPE % ElementCode
686      IF(ElementCode < 200) CYCLE
687      IF(ElementCode > 203) THEN
688        CALL Fatal('PhaseChangeSolve','Implemented only for elements 202 and 203!')
689        CYCLE
690      END IF
691
692      Nodes % x(1:n) = Solver % Mesh % Nodes % x(NodeIndexes)
693      Nodes % y(1:n) = Solver % Mesh % Nodes % y(NodeIndexes)
694      Nodes % z(1:n) = Solver % Mesh % Nodes % z(NodeIndexes)
695
696      TempDiff(1:n) = Temperature( TempPerm(NodeIndexes(1:n)) ) - MeltPoint
697      MaxTempDiff = MAX(MaxTempDiff, MAXVAL(ABS(TempDiff(1:n))))
698
699      DO nn=1,n
700
701        k = SurfPerm(NodeIndexes(nn))
702        IF ( NodeDone(k) ) CYCLE
703        NodeDone(k) = .TRUE.
704
705        IF( TriplePointFixed .AND. NodeIndexes(nn) == trip_node) THEN
706          SurfaceMove(k) = 0.0d0
707          PRINT *,'triple point fixed by construction'
708          CYCLE
709        END IF
710
711        IF ( nn == 3 ) THEN
712          Surface(k)=0.5*SUM(Surface(SurfPerm(NodeIndexes(1:2))))
713          Update = 0.0
714        END IF
715
716        IF(TangentDirection == 1) THEN
717          xx = Nodes % x(nn)
718          yy = Nodes % y(nn)
719        ELSE
720          xx = Nodes % y(nn)
721          yy = Nodes % x(nn)
722        END IF
723
724
725        IF ( .NOT. Newton ) THEN
726
727          ! Find the the contour element that has the x-coordinate in closest to the that of the
728          ! free surface
729
730          Eps = 1.0d-6 * ( xmax - xmin )
731
732          dxmin = HUGE(dxmin)
733          dymin = HUGE(dymin)
734          stat = .FALSE.
735
736          DO i=1,Nelems-1,2
737
738            IF ( (xx > IsoSurf(i,TangentDirection) - Eps) .AND. (xx < IsoSurf(i+1,TangentDirection) + Eps)) THEN
739              dxmin = 0.0
740              d = MIN( ABS(yy - IsoSurf(i,NormalDirection)), ABS(yy - IsoSurf(i+1,NormalDirection)) )
741
742              ! Punish for overlapping the boundaries
743              d = d + MAX(0.0d0, IsoSurf(i,TangentDirection) - xx)
744              d = d + MAX(0.0d0, xx - IsoSurf(i+1,TangentDirection) )
745
746              IF(d <= dymin) THEN
747                stat = .TRUE.
748                dymin = d
749                imin = i
750              END IF
751            END IF
752
753            IF(.NOT. stat) THEN
754              d = MIN( ABS(xx - IsoSurf(i,TangentDirection)), ABS(xx - IsoSurf(i+1,TangentDirection)) )
755              IF (d <= dxmin) THEN
756                dxmin = d
757                imin = i
758              END IF
759            END IF
760          END DO
761
762          i = imin
763
764          ! There may be a problem if the boundary cannot be mapped on an isotherm
765          IF (.NOT. stat) THEN
766            IF(dxmin > 1.0d-2* ABS(IsoSurf(i,TangentDirection)- IsoSurf(i+1,TangentDirection))) THEN
767              CALL Warn('PhaseChangeSolve','Isotherm error?')
768              WRITE(Message,*) 'Nodeindexes',NodeIndexes(nn)
769              CALL Warn('PhaseChangeSolve',Message)
770              WRITE(Message,*) 'x:',xx,' y:',yy
771              CALL Warn('PhaseChangeSolve',Message)
772              WRITE(Message,*) 'dxmin:',dxmin,' dymin:',dymin
773              CALL Warn('PhaseChangeSolve',Message)
774            END IF
775          END IF
776
777          IF ( ABS( IsoSurf(i+1,TangentDirection) - IsoSurf(i,TangentDirection) ) > AEPS ) THEN
778            Update = IsoSurf(i,NormalDirection) + ( xx - IsoSurf(i,TangentDirection) ) * &
779                ( IsoSurf(i+1,NormalDirection) - IsoSurf(i,NormalDirection) ) / &
780                ( IsoSurf(i+1,TangentDirection) - IsoSurf(i,TangentDirection) ) - yy
781          ELSE
782            Update = 0.5d0 * ( IsoSurf(i,NormalDirection) + IsoSurf(i+1,NormalDirection) ) - yy
783          END IF
784
785        END IF
786
787        TTemp = Temperature( TempPerm(NodeIndexes(nn)) )
788
789        IF ( Newton ) THEN
790          dTdz = TTemp - PrevTemp(k)
791          IF ( ABS(dTdz) < AEPS ) THEN
792            CALL Warn( 'PhaseChangeSolve', 'Very small temperature update.' )
793            dTdz = 1
794          END IF
795          Update = SurfaceMove(k) * ( MeltPoint - TTemp ) / dTdz
796        END IF
797
798        ! This enforcing is rather than by setting meltpoint to triple point temperature
799        ! IF ( NodeIndexes(nn) == Trip_node ) Update = 0
800
801        PrevTemp(k) = TTemp
802        SurfaceMove(k) = Update
803      END DO
804
805      IntegStuff = GaussPoints( CurrentElement )
806      DO i=1,IntegStuff % n
807
808        u = IntegStuff % u(i)
809        v = IntegStuff % v(i)
810        w = IntegStuff % w(i)
811
812        stat = ElementInfo( CurrentElement, Nodes, u, v, w, detJ, Basis, dBasisdx )
813
814        s = IntegStuff % s(i) * detJ
815
816        IF ( CurrentCoordinateSystem() /= Cartesian ) THEN
817          s = s * SUM(Basis(1:n) * Nodes % x(1:n)) * 2.0 * PI
818        END IF
819
820        area = area + S
821        volume = volume + S * SUM(Basis(1:n) * SurfaceMove(SurfPerm(NodeIndexes(1:n))))
822        tave = tave + S * SUM(Basis(1:n) * TempDiff(1:n) )
823        volabs = volabs + S * SUM(Basis(1:n) * ABS(SurfaceMove(SurfPerm(NodeIndexes(1:n)))))
824        tabs = tabs + S * SUM(Basis(1:n) * ABS(TempDiff(1:n)) )
825      END DO
826
827    END DO
828
829
830!    IF(.NOT. UseTAverage .AND. dt < SteadyDt) THEN
831!      LocaRelax = Relax * GetCReal( Solver % Values,'Temperature Relaxation Factor')
832!    ELSE
833      LocalRelax = Relax
834!    END IF
835
836    ! There are several different acceleration methods which are mainly inactive
837    tave = tave / area
838    tabs = tabs / area
839    volume = volume / area
840    volabs = volabs / area
841
842    i = ListGetInteger(Solver % Values,'Lumped Newton After Iterations', Stat)
843    IF(Stat .AND. SubroutineVisited > i) THEN
844
845      j = ListGetInteger(Solver % Values,'Lumped Newton Mode', Stat)
846      SELECT CASE( j )
847      CASE( 1 )
848        cvol = 0.5*(prevtave+tave)/(prevtave-tave)
849
850      CASE( 2 )
851        cvol = 0.5*(prevvolabs+volabs)/(prevvolabs-volabs)
852
853      CASE( 3 )
854        cvol = 0.5*(prevtabs+tabs)/(prevtabs-tabs)
855
856      CASE DEFAULT
857        cvol = 0.5*(prevvolume+volume)/(prevvolume-volume)
858
859      END SELECT
860
861      IF(cvol < 0.0) THEN
862        cvol = 1.0
863        ccum = 1.0
864      END IF
865
866      clim = ListGetConstReal(Solver % Values,'Lumped Newton Limiter', Stat)
867      IF(.NOT. Stat) clim = 100.0
868      cvol = MIN(clim,cvol)
869      cvol = MAX(1.0/clim,cvol)
870
871      ccum = ccum * cvol
872
873      WRITE(Message,*) 'Lumped Newton relaxation: ', ccum
874      CALL Info('PhaseChangeSolve',Message)
875
876      LocalRelax = LocalRelax * ccum
877    END IF
878
879    SurfaceMove = LocalRelax * SurfaceMove
880    Surface = Surface + SurfaceMove
881
882    dpos = SurfaceMove(SurfPerm(Trip_node))
883
884    MaxUpdate = MAXVAL(ABS(SurfaceMove)) / MAXVAL(ABS(Surface))
885
886    IF ( ABS(MaxUpdate) < NewtonAfterTol ) Newton = .TRUE.
887    WRITE(Message,*) 'Maximum surface update: ', MaxUpdate
888    CALL Info('PhaseChangeSolve',Message)
889
890    WRITE(Message,*) 'Maximum temperature difference: ', MaxTempDiff
891    CALL Info('PhaseChangeSolve',Message)
892
893    Norm = SQRT( SUM( Surface**2 ) / SIZE( Surface ) )
894    WRITE( Message, * ) 'Result Norm     : ',Norm
895    CALL Info( 'PhaseChangeSolve', Message, Level=4 )
896    Solver % Variable % Norm = Norm
897
898    prevvolume = volume
899    prevtave = tave
900    prevtabs = tabs
901    prevvolabs = volabs
902
903    IF(IsoSurfAllocated) THEN
904      IsoSurfAllocated = .FALSE.
905      DEALLOCATE(IsoSurf)
906    END IF
907  END IF
908
909!--------------------------------------------------------------------
910! Steady algorithm end
911!--------------------------------------------------------------------
912
913
914200 CALL ListAddConstReal(Model % Simulation,'res: Triple point temperature',Trip_temp)
915  CALL ListAddConstReal( Model % Simulation,'res: triple point movement',dpos)
916  CALL ListAddConstReal(Model % Simulation,'res: Pull Position',pos0)
917
918  IF(PullControl) THEN
919    IF(NormalDirection == 1) THEN
920      CALL ListAddConstReal( Model % Simulation,'res: Pull Velocity 1',UPull(NormalDirection))
921    ELSE
922      CALL ListAddConstReal( Model % Simulation,'res: Pull Velocity 2',UPull(NormalDirection))
923    END IF
924  END IF
925
926  FirstTime = .FALSE.
927
928!------------------------------------------------------------------------------
929CONTAINS
930
931!------------------------------------------------------------------------------
932  SUBROUTINE VelocityLocalMatrix( StiffMatrix, MassMatrix, ForceVector,&
933      Element, nCoord, Nodes )
934
935    ! external variables:
936    REAL(KIND=dp) :: StiffMatrix(:,:), MassMatrix(:,:), ForceVector(:)
937    INTEGER :: nCoord
938    TYPE(Nodes_t) :: Nodes
939    TYPE(Element_t), POINTER :: Element
940
941    ! internal variables:
942    TYPE(Nodes_t) :: PNodes
943    TYPE(Element_t), POINTER :: Parent
944    REAL(KIND=dp) :: Basis(3*nCoord),dBasisdx(3*nCoord,3), &
945        X,Y,Z,U,V,W,S,detJ, TGrad(3,3),Flux,pu,pv,pw,pull, LocalHeat, &
946        NodalTemp(3*nCoord), xx(10),yy(10),zz(10),  NodalSurf(nCoord), NodalNormal(2,nCoord)
947    REAL(KIND=dp) :: Velo, Ny, Normal(3), StabFactor, StabCoeff, DerSurf, xcoord
948    LOGICAL :: Stat
949    INTEGER :: i,j,k,l,t,p,q, n, pn
950    TYPE(GaussIntegrationPoints_t) :: IntegStuff
951
952!------------------------------------------------------------------------------
953
954    ForceVector = 0.0d0
955    StiffMatrix = 0.0d0
956    MassMatrix  = 0.0d0
957    n = nCoord
958
959    StabFactor = ListGetConstReal(Solver % Values,'Velocity Smoothing Factor')
960
961    Nodes % x(1:n) = Solver % Mesh % Nodes % x(NodeIndexes)
962    Nodes % y(1:n) = Solver % Mesh % Nodes % y(NodeIndexes)
963    NodalSurf(1:n) = Surface( SurfPerm(NodeIndexes) )
964
965    IF(AverageNormal) THEN
966      NodalNormal(1,1:n) = Normals(2*NormalsPerm(NodeIndexes(1:n))-1)
967      NodalNormal(2,1:n) = Normals(2*NormalsPerm(NodeIndexes(1:n)))
968    END IF
969
970    IF(TransientAlgo) THEN
971      ALLOCATE( PNodes % x(10), PNodes % y(10), PNodes % z(10) )
972    END IF
973
974!      Numerical integration:
975!      ----------------------
976
977    IntegStuff = GaussPoints( Element )
978
979    DO t = 1,IntegStuff % n
980
981      u = IntegStuff % u(t)
982      v = IntegStuff % v(t)
983      w = IntegStuff % w(t)
984      s = IntegStuff % s(t)
985
986!        Basis function values & derivatives at the integration point:
987!        -------------------------------------------------------------
988      stat = ElementInfo( Element,Nodes,U,V,W,detJ,Basis,dBasisdx)
989      s = s * detJ
990      xcoord = SUM( Nodes % x(1:nCoord) * Basis(1:nCoord) )
991
992!      IF ( CurrentCoordinateSystem() /= Cartesian ) THEN
993!        s = s * xcoord
994!      END IF
995
996      IF(AverageNormal) THEN
997        Normal(1) = SUM( Basis(1:n) * NodalNormal(1,1:n))
998        Normal(2) = SUM( Basis(1:n) * NodalNormal(2,1:n))
999        Normal(3) = 0.0d0
1000      ELSE
1001        Normal = NormalVector( Element, Nodes, u, v, .TRUE. )
1002      END IF
1003
1004      LocalHeat = SUM(Basis(1:n) * LatentHeat(1:n))
1005      StabCoeff = StabFactor * LocalHeat * Density
1006
1007
1008      IF(UseLoads) THEN
1009        ! do nothing, loads already computed
1010
1011      ELSE
1012        ! Compute the flux from normal derivaties
1013        TGrad = 0.0d0
1014        l = 0
1015        DO i=1,2
1016
1017          IF( i == 1) THEN
1018            Parent => Element % BoundaryInfo % Left
1019          ELSE
1020            Parent => Element % BoundaryInfo % Right
1021          END IF
1022
1023          k = ListGetInteger(Model % Bodies(Parent % BodyId) % Values,'Material')
1024          IF (ListGetLogical(Model % Materials(k) % Values,'Solid',stat)) THEN
1025            IF(l == 2) CALL Fatal('PhaseChangeSolve','Both materials cannot be solid!')
1026            l = 2
1027          ELSE
1028            IF(l == 1) CALL Fatal('PhaseChangeSolve','Both materials cannot be liquid!')
1029            l = 1
1030          END IF
1031
1032          pn = Parent % TYPE % NumberOfNodes
1033          k = ListGetInteger(Model % Bodies(Parent % BodyId) % Values,'Material')
1034
1035          Conductivity(1:pn) = ListGetReal( Model % Materials(k) % Values, &
1036              'Heat Conductivity', pn, Parent % NodeIndexes )
1037          NodalTemp(1:pn) = Temperature( TempPerm(Parent % NodeIndexes) )
1038
1039          stat = ElementInfo( Element, Nodes, u, v, w, detJ, Basis, dBasisdx)
1040
1041          !           Calculate the basis functions for the parent element:
1042          !           -----------------------------------------------------
1043          DO j = 1,n
1044            DO k = 1,pn
1045              IF( NodeIndexes(j) == Parent % NodeIndexes(k) ) THEN
1046                xx(j) = Parent % TYPE % NodeU(k)
1047                yy(j) = Parent % TYPE % NodeV(k)
1048                zz(j) = Parent % TYPE % NodeW(k)
1049                EXIT
1050              END IF
1051            END DO
1052          END DO
1053
1054          pu = SUM( Basis(1:n) * xx(1:n) )
1055          pv = SUM( Basis(1:n) * yy(1:n) )
1056          pw = SUM( Basis(1:n) * zz(1:n) )
1057
1058          PNodes % x(1:pn) = Solver % Mesh % Nodes % x(Parent % NodeIndexes)
1059          PNodes % y(1:pn) = Solver % Mesh % Nodes % y(Parent % NodeIndexes)
1060          PNodes % z(1:pn) = Solver % Mesh % Nodes % z(Parent % NodeIndexes)
1061
1062          stat = ElementInfo( Parent, PNodes, pu, pv, pw, detJ, Basis, dBasisdx )
1063
1064          DO j=1,DIM
1065            TGrad(l,j) = SUM( Conductivity(1:pn) * Basis(1:pn) ) * &
1066                SUM( dBasisdx(1:pn,j) * NodalTemp(1:pn) )
1067          END DO
1068        END DO
1069
1070        Flux = SUM( (TGrad(1,:) - TGrad(2,:)) *  Normal)
1071        stat = ElementInfo( Element,Nodes,U, V, W, detJ, Basis, dBasisdx )
1072      END IF
1073
1074      ! Assembly the matrix
1075      DO p=1,n
1076        DO q=1,n
1077          StiffMatrix(p,q) = StiffMatrix(p,q) + s * Basis(p) * Basis(q) * &
1078              Normal(NormalDirection) * Density * LocalHeat
1079          IF(NodeIndexes(p) /= Trip_node) THEN
1080            StiffMatrix(p,q) = StiffMatrix(p,q) + &
1081                s * StabCoeff * dBasisdx(q,TangentDirection) * dBasisdx(p,TangentDirection)
1082          END IF
1083
1084          ! BC for tipple node
1085!          IF(NodeIndexes(p) == Trip_node) THEN
1086!            StiffMatrix(p,q) = StiffMatrix(p,q) - &
1087!                StabCoeff * Basis(p) * dBasisdx(q,TangentDirection)
1088!          END IF
1089        END DO
1090
1091        ! transient part of heat flux
1092        IF(.NOT. (UseLoads)) THEN
1093          ForceVector(p) = ForceVector(p) - s * Basis(p) * Flux
1094        END IF
1095
1096      END DO
1097
1098    END DO
1099
1100    IF(TransientAlgo) DEALLOCATE( PNodes % x, PNodes % y, PNodes % z )
1101
1102!------------------------------------------------------------------------------
1103  END SUBROUTINE VelocityLocalMatrix
1104!------------------------------------------------------------------------------
1105
1106
1107!------------------------------------------------------------------------------
1108  SUBROUTINE SurfaceLocalMatrix( StiffMatrix, MassMatrix, ForceVector,&
1109      Element, nCoord, Nodes )
1110
1111    ! external variables:
1112    REAL(KIND=dp) :: StiffMatrix(:,:), MassMatrix(:,:), ForceVector(:)
1113    INTEGER :: nCoord
1114    TYPE(Nodes_t) :: Nodes
1115    TYPE(Element_t), POINTER :: Element
1116
1117    ! internal variables:
1118    REAL(KIND=dp) :: Basis(3*nCoord),dBasisdx(3*nCoord,3), &
1119        X,Y,Z,U,V,W,S,detJ,pull, NodalVelo(nCoord), xcoord, NodalNormal(2,nCoord)
1120    REAL(KIND=dp) :: Velo, StabFactor, StabCoeff
1121    LOGICAL :: Stat
1122    INTEGER :: i,j,k,l,t,p,q, n, pn
1123    TYPE(GaussIntegrationPoints_t) :: IntegStuff
1124
1125!------------------------------------------------------------------------------
1126
1127    ForceVector = 0.0d0
1128    StiffMatrix = 0.0d0
1129    MassMatrix  = 0.0d0
1130    n = nCoord
1131
1132    StabFactor = ListGetConstReal(Solver % Values,'Surface Smoothing Factor')
1133    StabCoeff = StabFactor / dt
1134
1135    Nodes % x(1:n) = Solver % Mesh % Nodes % x(NodeIndexes)
1136    Nodes % y(1:n) = Solver % Mesh % Nodes % y(NodeIndexes)
1137    NodalVelo(1:n) = SurfaceMove( SurfPerm(NodeIndexes) )
1138
1139    IF(AverageNormal) THEN
1140      NodalNormal(1,1:n) = Normals(2*NormalsPerm(NodeIndexes(1:n))-1)
1141      NodalNormal(2,1:n) = Normals(2*NormalsPerm(NodeIndexes(1:n)))
1142    END IF
1143
1144
1145!      Numerical integration:
1146!      ----------------------
1147
1148    IntegStuff = GaussPoints( Element )
1149
1150    DO t = 1,IntegStuff % n
1151
1152      u = IntegStuff % u(t)
1153      v = IntegStuff % v(t)
1154      w = IntegStuff % w(t)
1155      s = IntegStuff % s(t)
1156
1157!        Basis function values & derivatives at the integration point:
1158!        -------------------------------------------------------------
1159      stat = ElementInfo( Element,Nodes,U,V,W,detJ,Basis,dBasisdx)
1160      s = s * detJ
1161      xcoord = SUM( Nodes % x(1:nCoord) * Basis(1:nCoord) )
1162
1163      IF(AverageNormal) THEN
1164        Normal(1) = SUM( Basis(1:n) * NodalNormal(1,1:n) )
1165        Normal(2) = SUM( Basis(1:n) * NodalNormal(2,1:n) )
1166        Normal(3) = 0.0d0
1167      ELSE
1168        Normal = NormalVector( Element, Nodes, u, v, .TRUE. )
1169      END IF
1170
1171      Velo = SUM( Basis(1:n) * NodalVelo(1:n)) + Upull(NormalDirection)
1172      Velo = SpeedUp * Velo
1173
1174      DO p=1,n
1175        DO q=1,n
1176          MassMatrix(p,q) = MassMatrix(p,q) + s * Basis(p) * Basis(q)
1177          StiffMatrix(p,q) = StiffMatrix(p,q) + &
1178              s * StabCoeff * dBasisdx(q,TangentDirection) * dBasisdx(p,TangentDirection)
1179
1180          ! BC for tipple node
1181          IF(NodeIndexes(p) == Trip_node) THEN
1182            StiffMatrix(p,q) = StiffMatrix(p,q) - &
1183                StabCoeff * Basis(p) * dBasisdx(q,TangentDirection)
1184          END IF
1185        END DO
1186        ForceVector(p) = ForceVector(p) + s * Basis(p) * Velo
1187      END DO
1188
1189    END DO
1190
1191!------------------------------------------------------------------------------
1192  END SUBROUTINE SurfaceLocalMatrix
1193!------------------------------------------------------------------------------
1194
1195
1196!------------------------------------------------------------------------------
1197    SUBROUTINE FindPullBoundary()
1198
1199      REAL (KIND=dp) :: Ybot, Ymax, Ytop, Xtrip, x, y
1200      INTEGER :: t,k,i,n
1201      LOGICAL :: PullBoundary
1202
1203      IF(TangentDirection /= 1) CALL Fatal('PhaseChangeSolve',&
1204          'FindPullBoundary implemented only for lateral boundaries!')
1205
1206      PullBoundary = .FALSE.
1207      Ybot = Solver % Mesh % Nodes % y(Trip_node)
1208      Ymax = MAXVAL(Solver % Mesh % Nodes % y)
1209      Ytop = Ymax
1210      Xtrip = Solver % Mesh % Nodes % x(Trip_node)
1211
1212      DO t = Solver % Mesh % NumberOfBulkElements + 1, &
1213          Solver % Mesh % NumberOfBulkElements + Solver % Mesh % NumberOfBoundaryElements
1214
1215        CurrentElement => Solver % Mesh % Elements(t)
1216        Model % CurrentElement => CurrentElement
1217        n = CurrentElement % TYPE % NumberOfNodes
1218        NodeIndexes => CurrentElement % NodeIndexes
1219
1220        DO k=1, Model % NumberOfBCs
1221          IF ( Model % BCs(k) % Tag /= CurrentElement % BoundaryInfo % Constraint ) CYCLE
1222          IF( ListGetLogical(Model % BCs(k) % Values,'Pull Boundary',stat ) ) THEN
1223            PullBoundary = .TRUE.
1224            DO i = 1,n
1225              x = Solver % Mesh % Nodes % x(NodeIndexes(i))
1226              y = Solver % Mesh % Nodes % y(NodeIndexes(i))
1227              IF(y > Ybot .AND. ABS(x-Xtrip) > 1.0e-6 * (Ymax-Ybot)) Ytop = y
1228            END DO
1229          END IF
1230        END DO
1231      END DO
1232
1233      IF (PullBoundary) THEN
1234        CALL ListAddConstReal(Model % Simulation,'res: Triple point position',Ybot)
1235        CALL ListAddConstReal(Model % Simulation,'res: Full pull position',Ytop)
1236      END IF
1237
1238    END SUBROUTINE FindPullBoundary
1239
1240!------------------------------------------------------------------------------
1241  END SUBROUTINE PhaseChangeSolve
1242!------------------------------------------------------------------------------
1243
1244
1245!-------------------------------------------------------------------------------
1246  FUNCTION MeltingHeat(Model, Node, t) RESULT(Flux)
1247!-------------------------------------------------------------------------------
1248! This subroutine computes the heat flux resulting from solidification
1249! This
1250!-------------------------------------------------------------------------------
1251  USE Types
1252  USE Lists
1253  USE ElementDescription
1254  IMPLICIT NONE
1255!-------------------------------------------------------------------------------
1256  TYPE(Model_t) :: Model
1257  INTEGER:: Node
1258  REAL (KIND=dp):: Flux,t
1259!-------------------------------------------------------------------------------
1260  TYPE(Variable_t), POINTER :: NormalSol
1261  INTEGER:: k,n,i
1262  INTEGER, POINTER :: NodeIndexes(:), NormalsPerm(:)
1263  REAL (KIND=dp):: NodeLatentHeat, Density, NormalPull, u, v
1264  REAL (KIND=dp):: UPull(3) = (/ 0,0,0 /), Normal(3), ElemLatentHeat(4)
1265  REAL (KIND=dp), POINTER :: Normals(:)
1266  LOGICAL:: stat, NormalExist = .FALSE., Visited = .FALSE.
1267  TYPE(Nodes_t) :: Nodes
1268  TYPE(Element_t), POINTER :: CurrentElement, Parent
1269
1270!------------------------------------------------------------------------------
1271
1272  SAVE NormalExist, Normals, NormalsPerm, Nodes
1273
1274
1275  IF(.NOT. Visited) THEN
1276    NormalSol  => VariableGet( Model % Variables, 'Normals',ThisOnly=.TRUE. )
1277    IF(ASSOCIATED(NormalSol)) THEN
1278      NormalsPerm => NormalSol % Perm
1279      Normals => NormalSol % Values
1280      NormalExist = .TRUE.
1281    ELSE
1282      n = Model % Mesh % MaxElementNodes
1283      ALLOCATE( Nodes % x(n), Nodes % y(n), Nodes % z(n) )
1284    END IF
1285    Visited = .TRUE.
1286  END IF
1287
1288  CurrentElement => Model % CurrentElement
1289  NodeIndexes => CurrentElement % NodeIndexes
1290  n = CurrentElement % TYPE % NumberOfNodes
1291
1292  k = ListGetInteger(Model % Bodies(CurrentElement % BodyId) % Values,'Material')
1293  ElemLatentHeat(1:n) = ListGetReal( Model % Materials(k) % Values, 'Latent Heat', n, NodeIndexes )
1294
1295  DO i=1,n
1296    IF(NodeIndexes(i) == Node) EXIT
1297  END DO
1298  IF(NodeIndexes(i) /= Node) CALL Fatal('PhaseProcs','Node not found')
1299  NodeLatentHeat = ElemLatentHeat(i)
1300
1301  UPull = 0.0
1302  UPull(1) = ListGetConstReal(Model % Simulation,'res: Pull Velocity 1',stat)
1303  IF(.NOT. stat) UPull(1) = ListGetConstReal(Model % Materials(k) % Values,'Convection Velocity 1',stat)
1304
1305  UPull(2) = ListGetConstReal(Model % Simulation,'res: Pull Velocity 2',stat)
1306  IF(.NOT. stat) UPull(2) = ListGetConstReal(Model % Materials(k) % Values,'Convection Velocity 2',stat)
1307
1308  Parent => CurrentElement % BoundaryInfo % Left
1309  k = ListGetInteger(Model % Bodies(Parent % BodyId) % Values,'Material')
1310  IF (ListGetLogical(Model % Materials(k) % Values, 'Solid', stat)) THEN
1311    Density = ListGetConstReal( Model % Materials(k) % Values, 'Density' )
1312  ELSE
1313    Parent => CurrentElement % BoundaryInfo % Right
1314    k = ListGetInteger(Model % Bodies(Parent % BodyId) % Values,'Material')
1315    Density = ListGetConstReal( Model % Materials(k) % Values, 'Density' )
1316  END IF
1317
1318  IF(.NOT. ASSOCIATED(Parent)) CALL Fatal('PhaseProcs','Parent not found')
1319!------------------------------------------------------------------------------
1320
1321  IF( NormalExist ) THEN
1322    Normal(1) = Normals( 2*NormalsPerm( Node ) - 1 )
1323    Normal(2) = Normals( 2*NormalsPerm( Node ))
1324    Normal(3) = 0.0d0
1325  ELSE
1326    Nodes % x(1:n) = Model % Nodes % x(NodeIndexes)
1327    Nodes % y(1:n) = Model % Nodes % y(NodeIndexes)
1328    Nodes % z(1:n) = Model % Nodes % z(NodeIndexes)
1329
1330    ! For line segments the normal in the center is usually sufficient
1331    u = 0.0d0
1332    v = 0.0d0
1333
1334    ! If inner boundary, Normal Target Body should be defined for the boundary
1335    ! (if not, material density will be used to determine then normal direction
1336    ! and should be defined for bodies on both sides):
1337
1338    Normal = NormalVector( CurrentElement, Nodes, u, v, .TRUE. )
1339  END IF
1340
1341  NormalPull = SUM( Normal * UPull )
1342
1343  Flux = NodeLatentHeat * Density * NormalPull
1344
1345!------------------------------------------------------------------------------
1346END FUNCTION MeltingHeat
1347!------------------------------------------------------------------------------
1348