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! *
35! ****************************************************************************/
36
37!> \ingroup Solvers
38!> \{
39
40!------------------------------------------------------------------------------
41!> Initialization of the main solver: AdvectionDiffusionSolver
42!------------------------------------------------------------------------------
43   SUBROUTINE FlowSolver_init( Model,Solver,Timestep,TransientSimulation )
44!------------------------------------------------------------------------------
45     USE DefUtils
46     IMPLICIT NONE
47!------------------------------------------------------------------------------
48     TYPE(Solver_t) :: Solver          !< Linear & nonlinear equation solver options
49     TYPE(Model_t), TARGET :: Model    !< All model information (mesh, materials, BCs, etc...)
50     REAL(KIND=dp) :: Timestep         !< Timestep size for time dependent simulations
51     LOGICAL :: TransientSimulation    !< Steady state or transient simulation
52!------------------------------------------------------------------------------
53!    Local variables
54!------------------------------------------------------------------------------
55     TYPE(ValueList_t), POINTER :: Params
56
57     Params => GetSolverParams()
58     CALL ListAddInteger( Params,'Time Derivative Order',1 )
59
60   END SUBROUTINE FlowSolver_init
61
62
63!------------------------------------------------------------------------------
64!> Solver for the Navier-Stokes equation in various different coordinate systems.
65!------------------------------------------------------------------------------
66   SUBROUTINE FlowSolver( Model,Solver,dt,TransientSimulation)
67!------------------------------------------------------------------------------
68    USE NavierStokes
69    USE NavierStokesGeneral
70    USE NavierStokesCylindrical
71    USE Adaptive
72    USE DefUtils
73    USE FreeSurface
74    USE ElementDescription, ONLY: GetEdgeMap
75!------------------------------------------------------------------------------
76    IMPLICIT NONE
77
78     TYPE(Model_t) :: Model
79     TYPE(Solver_t), TARGET :: Solver
80
81     REAL(KIND=dp) :: dt
82     LOGICAL :: TransientSimulation
83!------------------------------------------------------------------------------
84!    Local variables
85!------------------------------------------------------------------------------
86     TYPE(Matrix_t),POINTER :: StiffMatrix
87
88     INTEGER :: i,j,k,n,nb,nd,t,iter,LocalNodes,istat,q,m
89
90     TYPE(ValueList_t),POINTER :: Material, BC, BodyForce, Equation
91     TYPE(Nodes_t) :: ElementNodes
92     TYPE(Element_t),POINTER :: Element
93
94     REAL(KIND=dp) :: RelativeChange,UNorm,Gravity(3),AngularVelocity(3), &
95       Tdiff,s,Relaxation,NewtonTol,NewtonUBound,NonlinearTol, &
96       ReferencePressure=0.0, SpecificHeatRatio, &
97       PseudoCompressibilityScale=1.0, FreeSTol, res
98
99     INTEGER :: NSDOFs,NewtonIter,NewtonMaxIter,NonlinearIter,FreeSIter
100
101     TYPE(Variable_t), POINTER :: DensitySol, TimeVar
102     TYPE(Variable_t), POINTER :: FlowSol, TempSol, MeshSol
103
104     INTEGER, POINTER :: FlowPerm(:),TempPerm(:), MeshPerm(:)
105     REAL(KIND=dp), POINTER :: FlowSolution(:), Temperature(:), &
106        gWork(:,:), ForceVector(:), LayerThickness(:), &
107           SurfaceRoughness(:),MeshVelocity(:)
108
109     REAL(KIND=dp), POINTER :: TempPrev(:)
110     REAL(KIND=DP), POINTER :: Pwrk(:,:,:)
111
112     LOGICAL :: Stabilize,NewtonLinearization = .FALSE., GotForceBC, GotIt, &
113                  MBFlag, Convect  = .TRUE., NormalTangential, RelaxBefore, &
114                  divDiscretization, GradPDiscretization, ComputeFree=.FALSE., &
115                  Transient, Rotating, AnyRotating, OutOfPlaneFlow=.FALSE.,&
116                  RecheckNewton=.FALSE.
117
118! Which compressibility model is used
119     CHARACTER(LEN=MAX_NAME_LEN) :: CompressibilityFlag, StabilizeFlag, VarName
120     CHARACTER(LEN=MAX_NAME_LEN) :: LocalCoords, FlowModel
121     INTEGER :: CompressibilityModel, ModelCoords, ModelDim, NoActive
122     INTEGER :: body_id,bf_id,eq_id,DIM
123     INTEGER :: MidEdgeNodes(12), BrickFaceMap(6,4)
124     INTEGER, POINTER :: NodeIndexes(:), Indexes(:)
125     INTEGER, POINTER :: EdgeMap(:,:)
126
127
128     INTEGER, SAVE :: Timestep, SaveTimestep=-1
129     REAL(KIND=dp), ALLOCATABLE, SAVE :: pDensity0(:), pDensity1(:)
130!
131     LOGICAL :: AllocationsDone = .FALSE., FreeSurfaceFlag, &
132         PseudoPressureExists, PseudoCompressible, Bubbles, P2P1, &
133         Porous =.FALSE., PotentialForce=.FALSE., Hydrostatic=.FALSE., &
134         MagneticForce =.FALSE., UseLocalCoords, PseudoPressureUpdate, &
135         AllIncompressible
136
137
138     REAL(KIND=dp),ALLOCATABLE :: MASS(:,:),STIFF(:,:), LoadVector(:,:), &
139       Viscosity(:),FORCE(:), TimeForce(:), PrevDensity(:),Density(:),   &
140       U(:),V(:),W(:),MU(:),MV(:),MW(:), Pressure(:),Alpha(:),Beta(:),   &
141       ExtPressure(:),PrevPressure(:), HeatExpansionCoeff(:),            &
142       ReferenceTemperature(:), Permeability(:),Mx(:),My(:),Mz(:),       &
143       LocalTemperature(:), GasConstant(:), HeatCapacity(:),             &
144       LocalTempPrev(:),SlipCoeff(:,:), PseudoCompressibility(:),        &
145       PseudoPressure(:), Drag(:,:), PotentialField(:),    &
146       PotentialCoefficient(:)
147
148     SAVE U,V,W,MASS,STIFF,LoadVector,Viscosity, TimeForce,FORCE,ElementNodes,  &
149       Alpha,Beta,ExtPressure,Pressure,PrevPressure, PrevDensity,Density,       &
150       AllocationsDone,LocalNodes, HeatExpansionCoeff,ReferenceTemperature,     &
151       Permeability,Mx,My,Mz,LayerThickness, SlipCoeff, SurfaceRoughness,       &
152       LocalTemperature, GasConstant, HeatCapacity, LocalTempPrev,MU,MV,MW,     &
153       PseudoCompressibilityScale, PseudoCompressibility, PseudoPressure,       &
154       PseudoPressureExists, Drag, PotentialField, PotentialCoefficient, &
155       ComputeFree, Indexes
156
157      REAL(KIND=dp) :: at,at0,at1,totat,st,totst
158!------------------------------------------------------------------------------
159
160     INTERFACE
161        FUNCTION FlowBoundaryResidual( Model,Edge,Mesh,Quant,Perm,Gnorm ) RESULT(Indicator)
162          USE Types
163          TYPE(Element_t), POINTER :: Edge
164          TYPE(Model_t) :: Model
165          TYPE(Mesh_t), POINTER :: Mesh
166          REAL(KIND=dp) :: Quant(:), Indicator(2), Gnorm
167          INTEGER :: Perm(:)
168        END FUNCTION FlowBoundaryResidual
169
170        FUNCTION FlowEdgeResidual( Model,Edge,Mesh,Quant,Perm ) RESULT(Indicator)
171          USE Types
172          TYPE(Element_t), POINTER :: Edge
173          TYPE(Model_t) :: Model
174          TYPE(Mesh_t), POINTER :: Mesh
175          REAL(KIND=dp) :: Quant(:), Indicator(2)
176          INTEGER :: Perm(:)
177        END FUNCTION FlowEdgeResidual
178
179        FUNCTION FlowInsideResidual( Model,Element,Mesh,Quant,Perm,Fnorm ) RESULT(Indicator)
180          USE Types
181          TYPE(Element_t), POINTER :: Element
182          TYPE(Model_t) :: Model
183          TYPE(Mesh_t), POINTER :: Mesh
184          REAL(KIND=dp) :: Quant(:), Indicator(2), Fnorm
185          INTEGER :: Perm(:)
186        END FUNCTION FlowInsideResidual
187     END INTERFACE
188!------------------------------------------------------------------------------
189
190!------------------------------------------------------------------------------
191!    Get variables needed for solving the system
192!------------------------------------------------------------------------------
193     CALL Info('FlowSolver','Solving the Navier-Stokes equations',Level=6)
194
195
196     IF ( .NOT. ASSOCIATED( Solver % Matrix ) ) RETURN
197
198     CALL DefaultStart()
199
200!    Check for local coordinate system
201
202     LocalCoords = GetString( Solver % Values, 'Solver Coordinate System', &
203          UseLocalCoords )
204
205     IF ( UseLocalCoords ) THEN
206        ModelCoords = Coordinates
207        ModelDim = Model % DIMENSION
208        SELECT CASE ( LocalCoords )
209           CASE( 'cartesian 2d' )
210              Coordinates = 1
211              Model % DIMENSION = 2
212              CALL Info( 'FlowSolve', 'Solver Coordinate System is cartesian 2d', LEVEL=31 )
213           CASE( 'cartesian 3d' )
214              Coordinates = 1
215              Model % DIMENSION = 3
216              CALL Info( 'FlowSolve', 'Solver Coordinate System is cartesian 3d', LEVEL=31 )
217           CASE( 'axi symmetric' )
218              Coordinates = 4
219              Model % DIMENSION = 2
220              CALL Info( 'FlowSolve', 'Solver Coordinate System is axi symmetric', LEVEL=31 )
221           CASE( 'cylindric symmetric' )
222              Coordinates = 3
223              Model % DIMENSION = 3
224              CALL Info( 'FlowSolve', 'Solver Coordinate System is cylindric symmetric', LEVEL=31 )
225           CASE DEFAULT
226              CALL Warn( 'FlowSolve', 'Solver coordinate system not recognised, using original' )
227        END SELECT
228     END IF
229
230     ! check for Flow model, one of 'full', 'no convection', 'stokes':
231     ! ---------------------------------------------------------------
232     Transient = TransientSimulation
233     Convect = .TRUE.
234     FlowModel = GetString( GetSolverParams(), 'Flow Model', Gotit )
235
236     SELECT CASE(FlowModel)
237     CASE('no convection')
238       Convect = .FALSE.
239     CASE('stokes')
240       Convect = .FALSE.
241       Transient = .FALSE.
242     CASE DEFAULT
243       FlowModel = 'full'
244     END SELECT
245
246     DIM = CoordinateSystemDimension()
247
248     FlowSol => Solver % Variable
249     NSDOFs         =  FlowSol % DOFs
250     FlowPerm       => FlowSol % Perm
251     FlowSolution   => FlowSol % Values
252
253     VarName = GetVarName(FlowSol)
254
255     LocalNodes = COUNT( FlowPerm > 0 )
256     IF ( LocalNodes <= 0 ) RETURN
257
258     TempSol => VariableGet( Solver % Mesh % Variables, 'Temperature' )
259     IF ( ASSOCIATED( TempSol ) ) THEN
260       TempPerm     => TempSol % Perm
261       Temperature  => TempSol % Values
262       IF( Transient ) THEN
263         IF ( ASSOCIATED(TempSol % PrevValues) ) TempPrev => TempSol % PrevValues(:,1)
264       END IF
265     END IF
266
267     MeshSol => VariableGet( Solver % Mesh % Variables, 'Mesh Velocity')
268     NULLIFY( MeshVelocity )
269     IF ( ASSOCIATED( MeshSol ) ) THEN
270       MeshPerm     => MeshSol % Perm
271       MeshVelocity => MeshSol % Values
272     END IF
273
274     DensitySol => VariableGet( Solver % Mesh % Variables, 'Density' )
275
276!------------------------------------------------------------------------------
277
278     StiffMatrix => Solver % Matrix
279     ForceVector => StiffMatrix % RHS
280     UNorm = Solver % Variable % Norm
281
282!------------------------------------------------------------------------------
283!     Allocate some permanent storage, this is done first time only
284!------------------------------------------------------------------------------
285
286     IF ( .NOT.AllocationsDone .OR. Solver % MeshChanged ) THEN
287
288       N = Solver % Mesh % MaxElementDOFs
289
290       IF( AllocationsDone ) THEN
291          DEALLOCATE(                               &
292               U,  V,  W,                           &
293               MU, MV, MW,                          &
294               Indexes,                             &
295               Pressure,                            &
296               PrevPressure,                        &
297               PseudoCompressibility,               &
298               PrevDensity,Density,                 &
299               LayerThickness,                      &
300               SurfaceRoughness,                    &
301               Permeability,                        &
302               Mx,My,Mz,                            &
303               SlipCoeff, Drag,                     &
304               TimeForce,FORCE, Viscosity,          &
305               MASS,  STIFF,                        &
306               HeatExpansionCoeff,                  &
307               GasConstant, HeatCapacity,           &
308               ReferenceTemperature,                &
309               LocalTempPrev, LocalTemperature,     &
310               PotentialField, PotentialCoefficient, &
311               LoadVector, Alpha, Beta, &
312               ExtPressure, STAT=istat )
313       END IF
314
315       ALLOCATE( U(N),  V(N),  W(N),                     &
316                 MU(N), MV(N), MW(N),                    &
317                 Indexes( N ),                           &
318                 Pressure( N ),                          &
319                 PrevPressure( N ),                      &
320                 PseudoCompressibility( N ),             &
321                 PrevDensity(N),Density( N ),            &
322                 LayerThickness(N),                      &
323                 SurfaceRoughness(N),                    &
324                 Permeability(N),                        &
325                 Mx(N),My(N),Mz(N),                      &
326                 SlipCoeff(3,N), Drag(3,N),              &
327                 TimeForce( 2*NSDOFs*N ),                &
328                 FORCE( 2*NSDOFs*N ), Viscosity( N ), &
329                 MASS(  2*NSDOFs*N,2*NSDOFs*N ),&
330                 STIFF( 2*NSDOFs*N,2*NSDOFs*N ),&
331                 HeatExpansionCoeff(N),                  &
332                 GasConstant( N ), HeatCapacity( N ),    &
333                 ReferenceTemperature(N),                &
334                 LocalTempPrev(N), LocalTemperature(N),  &
335                 PotentialField( N ), PotentialCoefficient( N ), &
336                 LoadVector( 4,N ), Alpha( N ), Beta( N ), &
337                 ExtPressure( N ), STAT=istat )
338
339       Drag = 0.0d0
340       NULLIFY(Pwrk)
341
342       PseudoPressureExists = .FALSE.
343       AllIncompressible = .TRUE.
344
345       DO k=1,Model % NumberOfMaterials
346         Material => Model % Materials(k) % Values
347         CompressibilityFlag = ListGetString( Material, &
348             'Compressibility Model', GotIt)
349         IF (.NOT. gotIt ) CYCLE
350         IF( CompressibilityFlag == 'artificial compressible') THEN
351           PseudoPressureExists = .TRUE.
352           AllIncompressible = .FALSE.
353         ELSE IF ( CompressibilityFlag == 'incompressible' ) THEN
354           CONTINUE
355         ELSE
356           AllIncompressible = .FALSE.
357         END IF
358       END DO
359
360       IF ( PseudoPressureExists ) THEN
361          IF ( AllocationsDone ) THEN
362             DEALLOCATE( PseudoPressure )
363          END IF
364          n = SIZE( FlowSolution ) / NSDOFs
365          ALLOCATE( PseudoPressure(n),STAT=istat )
366       END IF
367
368!------------------------------------------------------------------------------
369!     This hack is needed  cause of the fluctuating pressure levels
370!------------------------------------------------------------------------------
371       IF( AllIncompressible ) THEN
372         CALL Info('FlowSolve','Enforcing relative pressure relaxation',Level=8)
373         CALL ListAddNewLogical( Solver % Values,'Relative Pressure Relaxation',.TRUE.)
374       END IF
375
376       IF ( istat /= 0 ) THEN
377         CALL Fatal( 'FlowSolve','Memory allocation error, Aborting.' )
378       END IF
379
380!------------------------------------------------------------------------------
381
382       AllocationsDone = .TRUE.
383     END IF
384!------------------------------------------------------------------------------
385
386
387     TimeVar => VariableGet( Solver % Mesh % Variables, 'Timestep')
388     Timestep = NINT(Timevar % Values(1))
389     IF ( SaveTimestep /= Timestep ) THEN
390       IF ( ALLOCATED(pDensity0) ) pDensity0 = pDensity1
391       SaveTimestep=Timestep
392     END IF
393
394!------------------------------------------------------------------------------
395!    Do some additional initialization, and go for it
396!------------------------------------------------------------------------------
397
398     gWork => ListGetConstRealArray( Model % Constants,'Gravity',GotIt)
399     IF ( GotIt ) THEN
400       Gravity = gWork(1:3,1)*gWork(4,1)
401     ELSE
402       Gravity    =  0.00D0
403       Gravity(2) = -9.81D0
404     END IF
405
406     AnyRotating = ListCheckPresentAnyBodyForce(Model,'Angular Velocity') .OR. &
407         ListCheckPresentAnyBodyForce(Model,'Angular Velocity 1') .OR. &
408         ListCheckPresentAnyBodyForce(Model,'Angular Velocity 2') .OR. &
409         ListCheckPresentAnyBodyForce(Model,'Angular Velocity 3')
410
411!------------------------------------------------------------------------------
412     P2P1 = .FALSE.
413     Bubbles   = ListGetLogical( Solver % Values,'Bubbles',GotIt )
414     Stabilize = ListGetLogical( Solver % Values,'Stabilize',GotIt )
415
416     StabilizeFlag = ListGetString( Solver % Values, &
417           'Stabilization Method', GotIt )
418     IF ( .NOT. GotIt ) THEN
419       IF ( Stabilize ) THEN
420          StabilizeFlag = 'stabilized'
421       ELSE IF ( Bubbles  ) THEN
422          StabilizeFlag = 'bubbles'
423       ELSE
424          StabilizeFlag = 'stabilized'
425       END IF
426     ELSE
427       IF (StabilizeFlag == 'p2/p1' .OR. StabilizeFlag == 'p2p1') THEN
428         P2P1 = .TRUE.
429         Bubbles = .FALSE.
430         Stabilize = .FALSE.
431       END IF
432     END IF
433
434     IF ( StabilizeFlag == 'bubbles' ) Bubbles = .TRUE.
435
436     DivDiscretization = ListGetLogical( Solver % Values, &
437              'Div Discretization', GotIt )
438
439     GradPDiscretization = ListGetLogical( Solver % Values, &
440              'Gradp Discretization', GotIt )
441
442     NonlinearTol = ListGetConstReal( Solver % Values, &
443        'Nonlinear System Convergence Tolerance',minv=0.0d0 )
444
445     NewtonTol = ListGetConstReal( Solver % Values, &
446        'Nonlinear System Newton After Tolerance', minv=0.0d0 )
447
448     !Option to switch back to picard if convergence exceeds certain tolerance
449     NewtonUBound = ListGetConstReal( Solver % Values, &
450        'Nonlinear System Newton Max Tolerance', GotIt )
451     IF(GotIt) RecheckNewton = .TRUE.
452
453     NewtonIter = ListGetInteger( Solver % Values, &
454        'Nonlinear System Newton After Iterations', minv=0 )
455     IF ( NewtonIter == 0 ) NewtonLinearization = .TRUE.
456
457     !Option to switch back to picard after NewtonMaxIter iterations
458     NewtonMaxIter = ListGetInteger( Solver % Values, &
459        'Nonlinear System Newton Max Iterations', GotIt )
460     RecheckNewton = RecheckNewton .OR. GotIt
461
462     IF (GetLogical( GetSolverParams(), &
463         'Nonlinear System Reset Newton',  GotIt)) NewtonLinearization=.FALSE.
464
465     NonlinearIter = ListGetInteger( Solver % Values, &
466        'Nonlinear System Max Iterations', minv=0 )
467
468     IF ( .NOT. ListCheckPresent( Solver % Values, &
469        'Nonlinear System Norm Dofs' ) ) THEN
470       CALL ListAddInteger( Solver % Values, 'Nonlinear System Norm DOFs', NSDOFs-1 )
471     END IF
472
473     FreeSTol = ListGetConstReal( Solver % Values, &
474        'Free Surface After Tolerance', GotIt, minv=0.0d0 )
475     IF ( .NOT. GotIt ) FreeSTol = HUGE(1.0d0)
476
477     FreeSIter = ListGetInteger( Solver % Values, &
478        'Free Surface After Iterations', GotIt, minv=0 )
479     IF ( .NOT. GotIt ) FreeSIter = 0
480
481!------------------------------------------------------------------------------
482!    Check if free surfaces present
483!------------------------------------------------------------------------------
484     FreeSurfaceFlag = .FALSE.
485     DO i=1,Model % NumberOfBCs
486       FreeSurfaceFlag = FreeSurfaceFlag .OR. ListGetLogical( &
487          Model % BCs(i) % Values,'Free Surface', GotIt )
488       IF ( FreeSurfaceFlag ) EXIT
489     END DO
490
491     CALL CheckCircleBoundary()
492!------------------------------------------------------------------------------
493
494     totat = 0.0d0
495     totst = 0.0d0
496
497
498     ! Initialize the pressure to be used in artificial compressibility
499     IF(PseudoPressureExists) THEN
500       PseudoPressure = FlowSolution(NSDOFs:SIZE(FlowSolution):NSDOFs)
501
502       WRITE(Message,'(A,T25,E15.4)') 'PseudoPressure mean: ',&
503           SUM(PseudoPressure)/SIZE(PseudoPressure)
504       CALL Info('FlowSolve',Message,Level=5)
505
506       PseudoCompressibilityScale = ListGetConstReal( Model % Simulation, &
507           'Artificial Compressibility Scaling',GotIt)
508
509       IF(.NOT.GotIt) PseudoCompressibilityScale = 1.0
510       IF(Transient) THEN
511         PseudoCompressibilityScale = PseudoCompressibilityScale / dt
512       END IF
513       PseudoPressureUpdate = ListGetLogical( Model % Simulation, &
514           'Pseudo Pressure Update',GotIt)
515       IF (.NOT.GotIt) PseudoPressureUpdate = .FALSE.
516     END IF
517
518     DO iter=1,NonlinearIter
519
520       IF (PseudoPressureExists .AND. PseudoPressureUpdate) &
521          PseudoPressure = FlowSolution(NSDOFs:SIZE(FlowSolution):NSDOFs)
522
523       at  = CPUTime()
524       at0 = RealTime()
525       at1 = RealTime()
526
527       CALL Info( 'FlowSolve', ' ', Level=4 )
528       CALL Info( 'FlowSolve', ' ', Level=4 )
529       CALL Info( 'FlowSolve', '-------------------------------------', Level=4 )
530       WRITE( Message, * ) 'NAVIER-STOKES ITERATION', iter
531       CALL Info( 'FlowSolve',Message, Level=4 )
532       CALL Info( 'FlowSolve','-------------------------------------', Level=4 )
533       CALL Info( 'FlowSolve', ' ', Level=4 )
534       CALL Info( 'FlowSolve','Starting Assembly...', Level=4 )
535
536!------------------------------------------------------------------------------
537       !CALL InitializeToZero( StiffMatrix, ForceVector )
538       CALL DefaultInitialize()
539!------------------------------------------------------------------------------
540
541       bf_id   = -1
542       body_id = -1
543
544       CALL StartAdvanceOutput( 'FlowSolve', 'Assembly: ' )
545       NoActive = GetNOFActive()
546
547       DO t = 1,NoActive
548
549         CALL AdvanceOutput( t, NoActive )
550!
551         Element => GetActiveElement(t)
552         NodeIndexes => Element % NodeIndexes
553
554!------------------------------------------------------------------------------
555
556         IF ( Element % BodyId /= body_id ) THEN
557           body_id = Element % BodyId
558
559           eq_id = ListGetInteger( Model % Bodies(body_id) % Values,'Equation', &
560                   minv=1, maxv=Model % NumberOfEquations )
561           Equation => Model % Equations(eq_id) % Values
562
563           bf_id = ListGetInteger( Model % Bodies(body_id) % Values, &
564              'Body Force', gotIt, 1, Model % NumberOfBodyForces )
565           IF( bf_id > 0 ) THEN
566             BodyForce => Model % BodyForces(bf_id) % Values
567           END IF
568
569           IF ( FlowModel == 'full' ) THEN
570             Convect = ListGetLogical( Equation,'NS Convect', GotIt )
571             IF ( .NOT. GotIt ) Convect = .TRUE.
572           ENDIF
573
574           k = ListGetInteger( Model % Bodies(body_id) % Values, 'Material', &
575                  minv=1, maxv=Model % NumberOfMaterials )
576           Material => Model % Materials(k) % Values
577
578!------------------------------------------------------------------------------
579           CompressibilityFlag = ListGetString( Material, &
580               'Compressibility Model', GotIt)
581           IF ( .NOT.GotIt ) CompressibilityModel = Incompressible
582           PseudoCompressible = .FALSE.
583!------------------------------------------------------------------------------
584           SELECT CASE( CompressibilityFlag )
585!------------------------------------------------------------------------------
586             CASE( 'incompressible' )
587               CompressibilityModel = Incompressible
588
589             CASE( 'perfect gas', 'perfect gas equation 1' )
590               CompressibilityModel = PerfectGas1
591
592             CASE( 'thermal' )
593               CompressibilityModel = Thermal
594
595             CASE( 'user defined' )
596               CompressibilityModel = UserDefined1
597
598             CASE( 'pressure dependent' )
599               CompressibilityModel = UserDefined2
600
601             CASE( 'artificial compressible' )
602               CompressibilityModel = Incompressible
603               PseudoCompressible = .TRUE.
604
605             CASE DEFAULT
606               CompressibilityModel = Incompressible
607!------------------------------------------------------------------------------
608           END SELECT
609!------------------------------------------------------------------------------
610
611           Gotit = .FALSE.
612           IF ( bf_id > 0 ) THEN
613             MagneticForce = ListGetLogical( BodyForce,'Lorentz Force', gotIt )
614             Hydrostatic = ListGetLogical( BodyForce,'Hydrostatic Pressure',gotIt )
615           END IF
616           IF ( .NOT. GotIt ) THEN
617             Hydrostatic = ListGetLogical( Equation,'Hydrostatic Pressure',gotIt )
618           END IF
619!------------------------------------------------------------------------------
620
621           Rotating = .FALSE.
622           IF( bf_id > 0 ) THEN
623             gWork => ListGetConstRealArray( BodyForce,'Angular Velocity',GotIt)
624             IF ( GotIt ) THEN
625               IF( Coordinates == Cartesian ) THEN
626                 AngularVelocity = gWork(1:3,1)
627                 Rotating = .TRUE.
628               ELSE
629                 CALL Fatal('FlowSolve','Rotating coordinate implemented only for cartesian coodinates')
630               END IF
631             ELSE
632               AngularVelocity = 0.0_dp
633             END IF
634           END IF
635         END IF
636!------------------------------------------------------------------------------
637
638         n = GetElementNOFNodes()
639         nb = GetElementNOFBDOFs()
640         nd = GetElementDOFs( Indexes )
641
642         CALL GetElementNodes( ElementNodes )
643
644         SELECT CASE( NSDOFs )
645           CASE(3)
646             U(1:nd) = FlowSolution(NSDOFs*FlowPerm(Indexes(1:nd))-2)
647             V(1:nd) = FlowSolution(NSDOFs*FlowPerm(Indexes(1:nd))-1)
648             W(1:nd) = 0.0_dp
649             IF (bf_id > 0 ) THEN
650               W(1:n)  = ListGetReal(BodyForce,'Out Of Plane Velocity',&
651                    n, NodeIndexes(1:n),OutOfPlaneFlow)
652               IF (.NOT.OutOfPlaneFlow) &
653                    W(1:n) = 0.0_dp
654             END IF
655           CASE(4)
656             U(1:nd) = FlowSolution(NSDOFs*FlowPerm(Indexes(1:nd))-3)
657             V(1:nd) = FlowSolution(NSDOFs*FlowPerm(Indexes(1:nd))-2)
658             W(1:nd) = FlowSolution(NSDOFs*FlowPerm(Indexes(1:nd))-1)
659         END SELECT
660
661         MU(1:nd) = 0.0d0
662         MV(1:nd) = 0.0d0
663         MW(1:nd) = 0.0d0
664         IF ( ASSOCIATED( MeshVelocity ) ) THEN
665            SELECT CASE( MeshSol % DOFs )
666            CASE(2)
667               IF ( ALL( MeshPerm( Indexes(1:nd) ) > 0 ) ) THEN
668                  MU(1:nd) = MeshVelocity(2*MeshPerm(Indexes(1:nd))-1)
669                  MV(1:nd) = MeshVelocity(2*MeshPerm(Indexes(1:nd))-0)
670               END IF
671
672            CASE(3)
673               IF ( ALL( MeshPerm( NodeIndexes ) > 0 ) ) THEN
674                  MU(1:nd) = MeshVelocity(3*MeshPerm(Indexes(1:nd))-2)
675                  MV(1:nd) = MeshVelocity(3*MeshPerm(Indexes(1:nd))-1)
676                  MW(1:nd) = MeshVelocity(3*MeshPerm(Indexes(1:nd))-0)
677               END IF
678            END SELECT
679         END IF
680
681         LocalTemperature = 0.0d0
682         LocalTempPrev    = 0.0d0
683         IF ( ASSOCIATED( TempSol ) ) THEN
684            IF ( ALL( TempPerm(NodeIndexes) > 0 ) ) THEN
685               LocalTemperature(1:nd) = Temperature( TempPerm(Indexes(1:nd)) )
686               IF ( Transient .AND. CompressibilityModel /= Incompressible) THEN
687                 LocalTempPrev(1:nd) = TempPrev( TempPerm(Indexes(1:nd)) )
688               END IF
689            END IF
690         END IF
691         ReferencePressure = 0.0d0
692
693         PrevDensity = 0.0d0
694         Density = 0.0d0
695!------------------------------------------------------------------------------
696         SELECT CASE( CompressibilityModel )
697!------------------------------------------------------------------------------
698           CASE( Incompressible )
699!------------------------------------------------------------------------------
700             Pressure(1:nd) = FlowSolution( NSDOFs*FlowPerm(Indexes(1:nd)) )
701             Density(1:n) = GetReal( Material, 'Density' )
702
703             IF(PseudoCompressible) THEN
704               Pressure(1:n) = GetReal( Material,'Artificial Pressure', GotIt )
705               IF(.NOT. GotIt) THEN
706                 Pressure(1:nd) = PseudoPressure(FlowPerm(Indexes(1:nd)))
707               ELSE
708                 Pressure(n+1:nd) = 0.0d0
709               END IF
710               PseudoCompressibility(1:n) = PseudoCompressibilityScale * &
711                   GetReal(Material,'Artificial Compressibility', GotIt )
712               IF(.NOT. gotIt) PseudoCompressibility(1:n) = 0.0d0
713             END IF
714
715!------------------------------------------------------------------------------
716           CASE( PerfectGas1 )
717
718              ! Use  ReferenceTemperature in .sif file for fixed temperature
719              ! field. At the moment can not have both fixed T ideal gas and
720              ! Boussinesq force:
721              !-------------------------------------------------------------
722              IF ( .NOT. ASSOCIATED( TempSol ) ) THEN
723                 LocalTemperature(1:n) = GetReal( Material, &
724                   'Reference Temperature' )
725                 LocalTempPrev = LocalTemperature
726              END IF
727
728              HeatCapacity(1:n) = GetReal( Material, 'Heat Capacity', GotIt )
729
730
731              ! Read Specific Heat Ratio:
732              !--------------------------
733              SpecificHeatRatio = ListGetConstReal( Material, &
734                     'Specific Heat Ratio', GotIt )
735              IF ( .NOT.GotIt ) SpecificHeatRatio = 5.d0/3.d0
736
737
738              ! For an ideal gas, \gamma, c_p and R are really a constant
739              ! GasConstant is an array only since HeatCapacity formally is
740              !------------------------------------------------------------
741              GasConstant(1:n) = ( SpecificHeatRatio - 1.d0 ) *  &
742                   HeatCapacity(1:n) / SpecificHeatRatio
743
744
745              ! For ideal gases take pressure deviation p_d as the
746              ! dependent variable: p = p_0 + p_d
747              ! Read p_0
748              !---------------------------------------------------
749              ReferencePressure = ListGetConstReal( Material, &
750                      'Reference Pressure', GotIt )
751              IF ( .NOT.GotIt ) ReferencePressure = 0.0d0
752
753              Pressure(1:nd) = FlowSolution(NSDOFs*FlowPerm(Indexes(1:nd)))
754              IF ( Transient ) THEN
755                PrevPressure(1:nd) = Solver % Variable % PrevValues( &
756                          NSDOFs*FlowPerm(Indexes(1:nd)),1 )
757              END IF
758              Density(1:n) = ( Pressure(1:n) + ReferencePressure ) / &
759                 ( GasConstant(1:n) * LocalTemperature(1:n) )
760
761           CASE( UserDefined1 )
762             Pressure(1:nd) = FlowSolution(NSDOFs*FlowPerm(Indexes(1:nd)) )
763             IF ( ASSOCIATED( DensitySol ) ) THEN
764               Density(1:nd) = DensitySol % Values( DensitySol % Perm(Indexes(1:nd)) )
765               IF ( Transient ) THEN
766                  PrevDensity(1:nd) = DensitySol % PrevValues( &
767                       DensitySol % Perm(Indexes(1:nd)),1)
768                END IF
769             ELSE
770               Density(1:n) = GetReal( Material,'Density' )
771               IF ( Transient ) THEN
772                 IF (.NOT.ALLOCATED(pDensity0))  THEN
773                   ALLOCATE(pDensity0(LocalNodes), &
774                            pDensity1(LocalNodes))
775                 END IF
776
777                 IF ( Timestep==1 ) &
778                     pDensity0(Indexes(1:n)) = Density(1:n)
779                 pDensity1(Indexes(1:n)) = Density(1:n)
780                 PrevDensity(1:n) = pDensity0(Indexes(1:n))
781               END IF
782             END IF
783
784           CASE( UserDefined2 )
785             Density(1:n) = GetReal( Material,'Density' )
786             Pressure(1:nd) = FlowSolution(NSDOFs*FlowPerm(Indexes(1:nd)))
787
788           CASE( Thermal )
789             Pressure(1:n) = FlowSolution(NSDOFs*FlowPerm(Indexes(1:nd)))
790
791             HeatExpansionCoeff(1:n) = GetReal( Material, &
792               'Heat Expansion Coefficient' )
793
794             ReferenceTemperature(1:n) = GetReal( Material, &
795               'Reference Temperature' )
796
797             Density(1:n) = GetReal( Material,'Density' )
798
799             IF( Transient ) THEN
800               PrevDensity(1:n) = Density(1:n) * ( 1 - HeatExpansionCoeff(1:n)  * &
801                  ( LocalTempPrev(1:n) - ReferenceTemperature(1:n) ) )
802             END IF
803             Density(1:n) = Density(1:n) * ( 1 - HeatExpansionCoeff(1:n)  * &
804                  ( LocalTemperature(1:n) - ReferenceTemperature(1:n) ) )
805!------------------------------------------------------------------------------
806         END SELECT
807!------------------------------------------------------------------------------
808
809!------------------------------------------------------------------------------
810!        Read in porous media defs
811!------------------------------------------------------------------------------
812         Porous = ListGetLogical( Material,'Porous Media', GotIt)
813         IF(Porous) THEN
814           CALL GetRealArray( Material,  Pwrk,'Porous Resistivity',GotIt)
815
816           IF( .NOT. GotIt ) THEN
817             Drag( 1,1:n) = GetReal( Material,'Porous Resistivity 1',GotIt )
818	     Drag( 2,1:n) = GetReal( Material,'Porous Resistivity 2',GotIt )
819             IF( NSDOFs -1 > 2 ) THEN
820   	       Drag( 3,1:n) = GetReal( Material,'Porous Resistivity 3',GotIt )
821             END IF
822           ELSE IF ( SIZE(Pwrk,1) == 1 ) THEN
823             DO i=1,NSDOFs-1
824               Drag( i,1:n ) = Pwrk( 1,1,1:n )
825             END DO
826           ELSE
827             DO i=1,MIN(NSDOFs,SIZE(Pwrk,1))
828               Drag(i,1:n) = Pwrk(i,1,1:n)
829             END DO
830           END IF
831         END IF
832
833!------------------------------------------------------------------------------
834!        Viscosity = Laminar viscosity
835!------------------------------------------------------------------------------
836         Viscosity(1:n) = GetReal( Material,'Viscosity' )
837
838!------------------------------------------------------------------------------
839!        Set body forces, if any
840!------------------------------------------------------------------------------
841         LoadVector = 0.0D0
842
843         IF ( bf_id > 0 ) THEN
844           HeatExpansionCoeff   = 0.0D0
845           ReferenceTemperature = 0.0D0
846
847!------------------------------------------------------------------------------
848!          Boussinesq body force & gravity
849!------------------------------------------------------------------------------
850           IF ( ListGetLogical( BodyForce,'Boussinesq',gotIt) ) THEN
851
852             HeatExpansionCoeff(1:n) = GetReal( Material, &
853                 'Heat Expansion Coefficient' )
854
855             ReferenceTemperature(1:n) = GetReal( Material, &
856                 'Reference Temperature' )
857
858             DO i=1,n
859               k = TempPerm(NodeIndexes(i))
860               IF ( k > 0 ) THEN
861                 IF ( Hydrostatic ) THEN
862                   Tdiff = 1 - HeatExpansionCoeff(i) * &
863                      (Temperature(k) - ReferenceTemperature(i))
864
865                   IF ( Tdiff <= 0.0D0 ) THEN
866                      CALL Warn( 'FlowSolve','Zero or negative density.' )
867                   END IF
868                 ELSE
869                   Tdiff = -HeatExpansionCoeff(i) * &
870                               (Temperature(k) - ReferenceTemperature(i))
871                 END IF
872
873                 LoadVector(1,i)   = Gravity(1) * Tdiff
874                 LoadVector(2,i)   = Gravity(2) * Tdiff
875                 IF ( NSDOFs > 3 ) THEN
876                   LoadVector(3,i) = Gravity(3) * Tdiff
877                 END IF
878               END IF
879             END DO
880           ELSE IF ( Hydrostatic ) THEN
881             LoadVector(1,1:n)   = Gravity(1)
882             LoadVector(2,1:n)   = Gravity(2)
883             IF ( NSDOFs > 3 ) LoadVector(3,1:n) = Gravity(3)
884           END IF
885!------------------------------------------------------------------------------
886           LoadVector(1,1:n) = LoadVector(1,1:n) + ListGetReal( BodyForce, &
887               'Flow Bodyforce 1',n,NodeIndexes,gotIt )
888
889           LoadVector(2,1:n) = LoadVector(2,1:n) + ListGetReal( BodyForce, &
890               'Flow Bodyforce 2',n,NodeIndexes,gotIt )
891
892           IF ( NSDOFs > 3 ) THEN
893             LoadVector(3,1:n) = LoadVector(3,1:n) + ListGetReal( BodyForce, &
894                 'Flow Bodyforce 3',n,NodeIndexes,gotIt )
895           END IF
896
897!------------------------------------------------------------------------------
898
899           PotentialForce = ListGetLogical( BodyForce,'Potential Force',gotIt)
900           IF(PotentialForce) THEN
901             PotentialField(1:n) = ListGetReal( BodyForce, &
902                 'Potential Field',n,NodeIndexes)
903             PotentialCoefficient(1:n) = ListGetReal( BodyForce, &
904                 'Potential Coefficient',n,NodeIndexes)
905           END IF
906
907
908
909!------------------------------------------------------------------------------
910         END IF ! of body forces
911
912!------------------------------------------------------------------------------
913!
914! NOTE: LoadVector is multiplied by density inside *Navier* routines
915!
916         IF ( Transient ) THEN
917           SELECT CASE( CompressibilityModel )
918           CASE( PerfectGas1 )
919             IF ( ASSOCIATED( TempSol ) ) THEN
920               DO i=1,n
921                 k = TempPerm(NodeIndexes(i))
922                 IF ( k > 0 ) THEN
923                    LoadVector(NSDOFs,i) = LoadVector(NSDOFs,i) + &
924                      ( Temperature(k) - TempPrev(k) ) / dt
925                 END IF
926               END DO
927             END IF
928           CASE( UserDefined1, Thermal )
929              DO i=1,n
930                LoadVector(NSDOFs,i) = LoadVector(NSDOFs,i) - &
931                  ( Density(i) - PrevDensity(i) ) / (Density(i)*dt)
932              END DO
933           END SELECT
934         END IF
935
936!------------------------------------------------------------------------------
937!        Get element local stiffness & mass matrices
938!------------------------------------------------------------------------------
939         SELECT CASE(Coordinates)
940         CASE( Cartesian )
941!------------------------------------------------------------------------------
942           SELECT CASE( CompressibilityModel )
943!------------------------------------------------------------------------------
944             CASE( Incompressible,PerfectGas1,UserDefined1,UserDefined2,Thermal)
945!------------------------------------------------------------------------------
946! Density needed for steady-state, also pressure for transient
947!------------------------------------------------------------------------------
948               CALL NavierStokesCompose( MASS,STIFF,FORCE, LoadVector, &
949                   Viscosity,Density,U,V,W,MU,MV,MW,ReferencePressure+Pressure(1:n), &
950                   LocalTemperature, Convect, StabilizeFlag, CompressibilityModel, &
951                   PseudoCompressible, PseudoCompressibility, GasConstant, Porous, &
952                   Drag, PotentialForce, PotentialField, PotentialCoefficient,  &
953                   MagneticForce, Rotating, AngularVelocity, DivDiscretization, &
954                   GradPDiscretization, NewtonLinearization, Element,n,ElementNodes)
955!------------------------------------------------------------------------------
956           END SELECT
957!------------------------------------------------------------------------------
958
959         CASE( Cylindric,CylindricSymmetric,AxisSymmetric )
960! Same comments as Cartesian
961!------------------------------------------------------------------------------
962           SELECT CASE( CompressibilityModel )
963!------------------------------------------------------------------------------
964             CASE( Incompressible,PerfectGas1)
965!------------------------------------------------------------------------------
966               CALL NavierStokesCylindricalCompose( &
967                   MASS,STIFF,FORCE, &
968                   LoadVector, Viscosity,Density,U,V,W,MU,MV,MW, &
969                   ReferencePressure+Pressure(1:n),LocalTemperature,&
970                   Convect, StabilizeFlag, CompressibilityModel /= Incompressible, &
971                   PseudoCompressible, PseudoCompressibility, GasConstant, Porous, Drag, &
972                   PotentialForce, PotentialField, PotentialCoefficient, &
973                   MagneticForce, divDiscretization, gradpDiscretization, NewtonLinearization, &
974                   Element,n,ElementNodes )
975!------------------------------------------------------------------------------
976             CASE DEFAULT
977               CALL Fatal('FlowSolver','Missing compressibility model in cylindrical coordinates')
978
979
980           END SELECT
981!------------------------------------------------------------------------------
982
983!------------------------------------------------------------------------------
984         CASE DEFAULT
985!------------------------------------------------------------------------------
986
987           SELECT CASE( CompressibilityModel )
988!------------------------------------------------------------------------------
989           CASE( Incompressible,PerfectGas1)
990
991             CALL NavierStokesGeneralCompose( &
992                 MASS,STIFF,FORCE, &
993                 LoadVector, Viscosity,Density,U,V,W,MU,MV,MW,Stabilize, &
994                 NewtonLinearization,Element,n,ElementNodes )
995
996           CASE DEFAULT
997             CALL Fatal('FlowSolver','Missing compressibility model in general coordinates')
998
999           END SELECT
1000
1001!------------------------------------------------------------------------------
1002         END SELECT
1003!------------------------------------------------------------------------------
1004!        If time dependent simulation, add mass matrix to global
1005!        matrix and global RHS vector
1006!------------------------------------------------------------------------------
1007         IF ( CompressibilityModel /= Incompressible .AND. &
1008                 StabilizeFlag == 'stabilized' ) THEN
1009            Bubbles = .TRUE.
1010            StabilizeFlag = 'bubbles'
1011         END IF
1012         IF ( Element % TYPE % BasisFunctionDegree <= 1 .AND. P2P1 ) THEN
1013            Bubbles = .TRUE.
1014            StabilizeFlag = 'bubbles'
1015         END IF
1016
1017         IF ( nb==0 .AND. Bubbles ) nb = n
1018
1019         TimeForce = 0.0_dp
1020         IF ( Transient ) THEN
1021!------------------------------------------------------------------------------
1022!          NOTE: the following will replace STIFF and FORCE
1023!          with the combined information
1024!------------------------------------------------------------------------------
1025           CALL Default1stOrderTime( MASS, STIFF, FORCE )
1026         END IF
1027
1028         IF ( nb > 0 ) THEN
1029            CALL NSCondensate( nd, nb, NSDOFs-1, STIFF, FORCE, TimeForce )
1030         END IF
1031
1032!------------------------------------------------------------------------------
1033!        Add local stiffness matrix and force vector to global matrix & vector
1034!------------------------------------------------------------------------------
1035         CALL DefaultUpdateEquations( STIFF, FORCE )
1036!------------------------------------------------------------------------------
1037      END DO
1038!------------------------------------------------------------------------------
1039
1040      CALL DefaultFinishBulkAssembly()
1041
1042      CALL Info( 'FlowSolve', 'Assembly done', Level=4 )
1043
1044!------------------------------------------------------------------------------
1045!     Neumann & Newton boundary conditions
1046!------------------------------------------------------------------------------
1047      NoActive = GetNOFBoundaryElements()
1048
1049      DO t = 1,NoActive
1050
1051        Element => GetBoundaryElement(t)
1052        IF ( .NOT. ActiveBoundaryElement() ) CYCLE
1053
1054        n = GetElementNOFNodes()
1055
1056        CALL GetElementNodes( ElementNodes )
1057        NodeIndexes => Element % NodeIndexes
1058
1059        BC => GetBC()
1060        IF ( .NOT. ASSOCIATED(BC) ) CYCLE
1061
1062!------------------------------------------------------------------------------
1063        GotForceBC = GetLogical( BC, 'Flow Force BC',gotIt )
1064        IF ( .NOT. gotIt ) GotForceBC = .TRUE.
1065
1066        IF ( GotForceBC ) THEN
1067          LoadVector  = 0.0d0
1068          Alpha       = 0.0d0
1069          ExtPressure = 0.0d0
1070          Beta        = 0.0d0
1071          SlipCoeff   = 0.0d0
1072          STIFF = 0.0d0
1073          FORCE = 0.0d0
1074
1075!------------------------------------------------------------------------------
1076!         (at the moment the following is done...)
1077!         BC: \tau \cdot n = \alpha n +  @\beta/@t + R_k u_k + F
1078!------------------------------------------------------------------------------
1079!------------------------------------------------------------------------------
1080!         normal force BC: \tau\cdot n = \alpha n
1081!------------------------------------------------------------------------------
1082          IF ( GetLogical( BC, 'Free Surface',gotIt) ) THEN
1083            Alpha(1:n) = GetReal( BC,'Surface Tension Coefficient', gotIt )
1084          END IF
1085
1086          ExtPressure(1:n) = GetReal( BC, 'External Pressure', GotForceBC )
1087          IF(.NOT. GotForceBC) ExtPressure(1:n) = GetReal( BC, 'Normal Pressure', GotForceBC )
1088!------------------------------------------------------------------------------
1089!         tangential force BC:
1090!         \tau\cdot n = @\beta/@t (tangential derivative of something)
1091!------------------------------------------------------------------------------
1092
1093          IF ( ASSOCIATED( TempSol ) ) THEN
1094            Beta(1:n) = GetReal( BC, &
1095                'Surface Tension Expansion Coefficient',gotIt )
1096
1097            IF ( gotIt ) THEN
1098              DO j=1,n
1099                k = TempPerm( NodeIndexes(j) )
1100                IF ( k>0 ) Beta(j) = 1.0_dp - Beta(j) * Temperature(k)
1101              END DO
1102              Beta(1:n) = Beta(1:n) * GetReal(BC, 'Surface Tension Coefficient' )
1103            ELSE
1104              Beta(1:n) = GetReal( BC,'Surface Tension Coefficient', gotIt )
1105            END IF
1106          END IF
1107
1108!------------------------------------------------------------------------------
1109!         force in given direction BC: \tau\cdot n = F
1110!------------------------------------------------------------------------------
1111
1112          LoadVector(1,1:n) =  GetReal( BC, 'Pressure 1', GotIt )
1113          LoadVector(2,1:n) =  GetReal( BC, 'Pressure 2', GotIt )
1114          LoadVector(3,1:n) =  GetReal( BC, 'Pressure 3', GotIt )
1115          LoadVector(4,1:n) =  GetReal( BC, 'Mass Flux', GotIt )
1116
1117!------------------------------------------------------------------------------
1118!         slip boundary condition BC: \tau\cdot n = R_k u_k
1119!------------------------------------------------------------------------------
1120
1121          SlipCoeff = 0.0d0
1122          SlipCoeff(1,1:n) =  GetReal( BC, 'Slip Coefficient 1',GotIt )
1123          SlipCoeff(2,1:n) =  GetReal( BC, 'Slip Coefficient 2',GotIt )
1124          SlipCoeff(3,1:n) =  GetReal( BC, 'Slip Coefficient 3',GotIt )
1125
1126          NormalTangential = GetLogical( BC, &
1127                 'Normal-Tangential Velocity', GotIt )
1128
1129          IF (.NOT.GotIt) THEN
1130            NormalTangential = GetLogical( BC, &
1131                   'Normal-Tangential '//GetVarName(Solver % Variable), GotIt )
1132          END IF
1133!------------------------------------------------------------------------------
1134          SELECT CASE( CurrentCoordinateSystem() )
1135          CASE( Cartesian )
1136
1137            CALL NavierStokesBoundary(  STIFF, FORCE, &
1138             LoadVector, Alpha, Beta, ExtPressure, SlipCoeff, NormalTangential,   &
1139                Element, n, ElementNodes )
1140
1141         CASE( Cylindric, CylindricSymmetric,  AxisSymmetric )
1142
1143            CALL NavierStokesCylindricalBoundary( STIFF, &
1144             FORCE, LoadVector, Alpha, Beta, ExtPressure, SlipCoeff, &
1145                 NormalTangential, Element, n, ElementNodes)
1146
1147         CASE DEFAULT
1148
1149            CALL NavierStokesGeneralBoundary( STIFF, &
1150             FORCE, LoadVector, Alpha, Beta, ExtPressure, SlipCoeff, &
1151                Element, n, ElementNodes)
1152
1153         END SELECT
1154
1155!------------------------------------------------------------------------------
1156
1157          IF ( GetLogical( BC, 'Wall Law',GotIt ) ) THEN
1158            !/*
1159            ! * TODO: note that the following is not really valid, the
1160            ! * pointer to the Material structure is from the remains
1161            ! * of the last of the bulk elements.
1162            ! */
1163            Density(1:n)   = GetParentMatProp( 'Density' )
1164            Viscosity(1:n) = GetParentMatProp( 'Viscosity' )
1165
1166            LayerThickness(1:n) = GetReal( BC, 'Boundary Layer Thickness' )
1167            SurfaceRoughness(1:n) = GetReal( BC, 'Surface Roughness',GotIt )
1168
1169            SELECT CASE( NSDOFs )
1170              CASE(3)
1171                U(1:n) = FlowSolution( NSDOFs*FlowPerm(NodeIndexes)-2 )
1172                V(1:n) = FlowSolution( NSDOFs*FlowPerm(NodeIndexes)-1 )
1173                W(1:n) = 0.0d0
1174
1175              CASE(4)
1176                U(1:n) = FlowSolution( NSDOFs*FlowPerm(NodeIndexes)-3 )
1177                V(1:n) = FlowSolution( NSDOFs*FlowPerm(NodeIndexes)-2 )
1178                W(1:n) = FlowSolution( NSDOFs*FlowPerm(NodeIndexes)-1 )
1179            END SELECT
1180
1181            CALL NavierStokesWallLaw( STIFF,FORCE,     &
1182              LayerThickness,SurfaceRoughness,Viscosity,Density,U,V,W, &
1183                     Element,n, ElementNodes )
1184          ELSE IF ( GetLogical( BC, 'VMS Wall', GotIt ) ) THEN
1185            Density(1:n)   = GetParentMatProp( 'Density' )
1186            Viscosity(1:n) = GetParentMatProp( 'Viscosity' )
1187
1188            LayerThickness(1:n) = GetReal( BC, 'Boundary Layer Thickness', GotIt )
1189            SurfaceRoughness(1:n) = GetReal( BC, 'Surface Roughness',GotIt )
1190
1191            SELECT CASE( NSDOFs )
1192              CASE(3)
1193                U(1:n) = FlowSolution( NSDOFs*FlowPerm(NodeIndexes)-2 )
1194                V(1:n) = FlowSolution( NSDOFs*FlowPerm(NodeIndexes)-1 )
1195                W(1:n) = 0.0d0
1196
1197              CASE(4)
1198                U(1:n) = FlowSolution( NSDOFs*FlowPerm(NodeIndexes)-3 )
1199                V(1:n) = FlowSolution( NSDOFs*FlowPerm(NodeIndexes)-2 )
1200                W(1:n) = FlowSolution( NSDOFs*FlowPerm(NodeIndexes)-1 )
1201            END SELECT
1202            CALL VMSWalls( STIFF,FORCE,     &
1203              LayerThickness,SurfaceRoughness,Viscosity,Density,U,V,W, &
1204                     Element,n, ElementNodes )
1205
1206          END IF
1207!------------------------------------------------------------------------------
1208
1209          IF ( Transient ) THEN
1210            MASS = 0.0d0
1211            CALL Default1stOrderTime( MASS, STIFF, FORCE )
1212          END IF
1213
1214!------------------------------------------------------------------------------
1215!         Add local stiffness matrix and force vector to
1216!         global matrix & vector
1217!--------------------------------------------------------------------------
1218          CALL DefaultUpdateEquations( STIFF, FORCE )
1219!------------------------------------------------------------------------------
1220        END IF
1221      END DO
1222!------------------------------------------------------------------------------
1223      !
1224      ! IMPLEMENT NOSLIP WALL BC CODE:
1225      ! ------------------------------
1226      DO i=1,Model % NumberOFBCs
1227        BC => Model % BCs(i) % Values
1228        IF ( GetLogical(  BC, 'Noslip wall BC', gotit ) ) THEN
1229          IF ( VarName  == 'flow solution' ) THEN
1230            CALL ListAddConstReal( BC, 'Velocity 1', 0.0_dp )
1231            CALL ListAddConstReal( BC, 'Velocity 2', 0.0_dp )
1232            IF ( NSDOFs>3 ) CALL ListAddConstReal( BC, 'Velocity 3', 0.0_dp )
1233          ELSE
1234            DO j=1,NSDOFs-1
1235              CALL ListAddConstReal( BC, ComponentName( &
1236                 Solver % Variable % Name, j), 0.0_dp )
1237            END DO
1238          END IF
1239        END IF
1240      END DO
1241
1242      CALL DefaultFinishBoundaryAssembly()
1243      CALL DefaultFinishAssembly()
1244
1245!------------------------------------------------------------------------------
1246!     Dirichlet boundary conditions
1247!------------------------------------------------------------------------------
1248      CALL DefaultDirichletBCs()
1249      CALL Info( 'FlowSolve', 'Dirichlet conditions done', Level=4 )
1250
1251!------------------------------------------------------------------------------
1252!     Solve the system and check for convergence
1253!------------------------------------------------------------------------------
1254      at = CPUTime() - at
1255      st = CPUTime()
1256
1257      Unorm = DefaultSolve()
1258
1259      st = CPUTIme()-st
1260      totat = totat + at
1261      totst = totst + st
1262      WRITE(Message,'(a,i4,a,F8.2,F8.2)') 'iter: ',iter,' Assembly: (s)', at, totat
1263      CALL Info( 'FlowSolve', Message, Level=4 )
1264      WRITE(Message,'(a,i4,a,F8.2,F8.2)') 'iter: ',iter,' Solve:    (s)', st, totst
1265      CALL Info( 'FlowSolve', Message, Level=4 )
1266
1267      n = NSDOFs * LocalNodes
1268
1269
1270!------------------------------------------------------------------------------
1271
1272      WRITE( Message, * ) 'Result Norm     : ',Solver % Variable % Norm
1273      CALL Info( 'FlowSolve', Message, Level=4 )
1274      WRITE( Message, * ) 'Relative Change : ',Solver % Variable % NonlinChange
1275      CALL Info( 'FlowSolve', Message, Level=4 )
1276
1277      RelativeChange = Solver % Variable % NonlinChange
1278      IF ( RelativeChange < NewtonTol .OR. &
1279             iter > NewtonIter ) NewtonLinearization = .TRUE.
1280
1281      IF ( RecheckNewton .AND. NewtonLinearization .AND. (RelativeChange > NewtonUBound)) THEN
1282        NewtonLinearization = .FALSE.
1283	CALL Info('FlowSolve', 'Newton tolerance exceeded, switching back to picard', Level=6)
1284      END IF
1285
1286      IF ( RecheckNewton .AND. NewtonLinearization .AND. (iter >= NewtonMaxIter)) THEN
1287        NewtonLinearization = .FALSE.
1288	CALL Info('FlowSolve', 'Newton iteration limit exceeded, switching back to picard', Level=6)
1289      END IF
1290
1291      IF ( RecheckNewton .AND. (RelativeChange > NewtonUBound) .AND. NewtonLinearization ) THEN
1292        NewtonLinearization = .FALSE.
1293	CALL Info('FlowSolve', 'Newton tolerance exceeded, switching back to picard', Level=6)
1294      END IF
1295
1296      IF ( RelativeChange < NonLinearTol .AND. Iter<NonlinearIter ) EXIT
1297
1298!------------------------------------------------------------------------------
1299!     If free surfaces in model, this will move the nodal points
1300!------------------------------------------------------------------------------
1301      IF ( FreeSurfaceFlag ) THEN
1302
1303        IF ( RelativeChange < FreeSTol .OR. &
1304             iter > FreeSIter ) ComputeFree = .TRUE.
1305
1306        IF ( ComputeFree ) THEN
1307          Relaxation = GetCReal( Solver % Values, &
1308             'Free Surface Relaxation Factor', GotIt )
1309
1310          IF ( .NOT.GotIt ) Relaxation = 1.0d0
1311
1312          MBFlag = GetLogical(GetSolverParams(), 'Internal Move Boundary', GotIt)
1313          IF ( MBFlag .OR. .NOT. GotIt ) CALL MoveBoundary( Model, Relaxation )
1314        END IF
1315      END IF
1316!------------------------------------------------------------------------------
1317    END DO  ! of nonlinear iteration
1318
1319    IF ( P2P1 ) THEN
1320      !----------------------------------------------------------------------------------------
1321      ! Replace the zero pressure solution at the nodes which are not needed in the linear
1322      ! pressure approximation by the interpolated values for right visualization:
1323      !----------------------------------------------------------------------------------------
1324      NoActive = GetNOFActive()
1325
1326      DO t=1,NoActive
1327        ! First the midedge nodes:
1328        Element => GetActiveElement(t)
1329        IF ( Element % TYPE % BasisFunctionDegree <= 1 ) CYCLE
1330
1331        nd = GetElementDOFs( Indexes )
1332        k = GetElementFamily()
1333        EdgeMap => GetEdgeMap(k)
1334        SELECT CASE(k)
1335        CASE (3)
1336          MidEdgeNodes(1:3) = (/ 4, 5, 6 /)
1337        CASE (4)
1338          MidEdgeNodes(1:4) = (/ 5, 6, 7, 8 /)
1339        CASE (5)
1340          MidEdgeNodes(1:6) = (/ 5, 6, 7, 8, 9, 10 /)
1341        CASE (6)
1342          MidEdgeNodes(1:8) = (/ 6, 7, 8, 9, 10, 11, 12, 13 /)
1343        CASE (7)
1344          MidEdgeNodes(1:9) = (/ 7, 8, 9, 10, 11, 12, 13, 14, 15 /)
1345        CASE (8)
1346          MidEdgeNodes(1:12) = (/ 9, 10, 11, 12, 17, 18, 19, 20, 13, 14, 15, 16 /)
1347        END SELECT
1348
1349        DO q=1,SIZE(EdgeMap,1)
1350          m = (dim+1) * Solver % Variable % Perm(Indexes(MidEdgeNodes(q)))
1351          i = (dim+1) * Solver % Variable % Perm(Indexes(EdgeMap(q,1)))
1352          j = (dim+1) * Solver % Variable % Perm(Indexes(EdgeMap(q,2)))
1353          Solver % Variable % Values(m) = 0.5d0 * ( Solver % Variable % Values(i) + &
1354              Solver % Variable % Values(j) )
1355        END DO
1356
1357        ! The pressure at the midface nodes for 409 elements:
1358        IF (k==4 .AND. nd==9) THEN
1359          res = 0.0d0
1360          DO q=1,4
1361            i = (dim+1) * Solver % Variable % Perm(Indexes(q))
1362            res = res + Solver % Variable % Values(i)
1363          END DO
1364          m = (dim+1) * Solver % Variable % Perm(Indexes(9))
1365          Solver % Variable % Values(m) = 0.25d0 * res
1366        END IF
1367
1368        ! The pressure at the midpoint and at the midface nodes for 827 elements:
1369        IF (k==8 .AND. nd==27) THEN
1370          BrickFaceMap(1,:) = (/ 1,2,6,5 /)
1371          BrickFaceMap(2,:) = (/ 2,3,7,6 /)
1372          BrickFaceMap(3,:) = (/ 4,3,7,8 /)
1373          BrickFaceMap(4,:) = (/ 1,4,8,5 /)
1374          BrickFaceMap(5,:) = (/ 1,2,3,4 /)
1375          BrickFaceMap(6,:) = (/ 5,6,7,8 /)
1376          DO j=1,6
1377            res = 0.0d0
1378            DO q=1,4
1379              i = (dim+1) * Solver % Variable % Perm(Indexes(BrickFaceMap(j,q)))
1380              res = res + Solver % Variable % Values(i)
1381            END DO
1382            m = (dim+1) * Solver % Variable % Perm(Indexes(20+j))
1383            Solver % Variable % Values(m) = 0.25d0 * res
1384          END DO
1385
1386          res = 0.0d0
1387          DO q=1,8
1388            i = (dim+1) * Solver % Variable % Perm(Indexes(q))
1389            res = res + Solver % Variable % Values(i)
1390          END DO
1391          m = (dim+1) * Solver % Variable % Perm(Indexes(27))
1392          Solver % Variable % Values(m) = 0.125d0 * res
1393        END IF
1394
1395      END DO
1396    END IF
1397
1398    IF (ListGetLogical(Solver % Values,'Adaptive Mesh Refinement',GotIt)) &
1399      CALL RefineMesh( Model,Solver,FlowSolution,FlowPerm, &
1400         FlowInsideResidual, FlowEdgeResidual, FlowBoundaryResidual )
1401
1402!------------------------------------------------------------------------------
1403    CALL CheckCircleBoundary()
1404!------------------------------------------------------------------------------
1405
1406    IF ( UseLocalCoords ) THEN
1407       Coordinates = ModelCoords
1408       Model % DIMENSION = ModelDim
1409    END IF
1410
1411CONTAINS
1412!------------------------------------------------------------------------------
1413
1414
1415!------------------------------------------------------------------------------
1416   SUBROUTINE CheckCircleBoundary()
1417!------------------------------------------------------------------------------
1418      REAL(KIND=dp) :: x,y,phi,x0,y0,r
1419      LOGICAL :: GotIt
1420      INTEGER :: i,j,k,l
1421!------------------------------------------------------------------------------
1422
1423      l = 0
1424      DO i=1,Model % NumberOfBCs
1425         IF ( .NOT.ListgetLogical( Model % BCs(i) % Values, &
1426                  'Circle Boundary', GotIt ) ) CYCLE
1427
1428         x0 = ListGetConstReal( Model % BCs(i) % Values, 'Circle X', GotIt )
1429         IF ( .NOT. GotIt ) x0 = 0.0d0
1430
1431         y0 = ListGetConstReal( Model % BCs(i) % Values, 'Circle Y', GotIt )
1432         IF ( .NOT. GotIt ) y0 = 0.0d0
1433
1434         R  = ListGetConstReal( Model % BCs(i) % Values, 'Circle R', GotIt )
1435         IF ( .NOT. GotIt ) R = 1.0d0
1436
1437         DO j=Solver % Mesh % NumberOfBulkElements+1, &
1438            Solver % Mesh % NumberOfBulkElements+ &
1439               Solver % Mesh % NumberOfBoundaryElements
1440            Element => Solver % Mesh % Elements(j)
1441            IF ( Element % BoundaryInfo % Constraint &
1442                 /= Model % BCs(i) % Tag ) CYCLE
1443
1444            n = Element % TYPE % NumberOfNodes
1445            NodeIndexes => Element % NodeIndexes
1446            DO k=1,n
1447               x = Solver % Mesh % Nodes % x(NodeIndexes(k)) - x0
1448               y = Solver % Mesh % Nodes % y(NodeIndexes(k)) - y0
1449
1450               phi = ATAN2( y,x )
1451               x = R * COS( phi )
1452               y = R * SIN( phi )
1453
1454               Solver % Mesh % Nodes % x(NodeIndexes(k)) = x + x0
1455               Solver % Mesh % Nodes % y(NodeIndexes(k)) = y + y0
1456            END DO
1457            l = l + 1
1458        END DO
1459     END DO
1460
1461     IF ( l > 0 ) THEN
1462        WRITE( Message, * ) 'Elements on Circle', l
1463        CALL Info( 'FlowSolve', Message, Level=6 )
1464     END IF
1465!------------------------------------------------------------------------------
1466   END SUBROUTINE CheckCircleBoundary
1467!------------------------------------------------------------------------------
1468
1469!------------------------------------------------------------------------------
1470  END SUBROUTINE FlowSolver
1471!------------------------------------------------------------------------------
1472
1473!------------------------------------------------------------------------------
1474!> Compute the residual of the Navier-Stokes equation for the boundary elements.
1475!------------------------------------------------------------------------------
1476  FUNCTION FlowBoundaryResidual( Model, Edge, Mesh, &
1477        Quant, Perm, Gnorm ) RESULT( Indicator )
1478!------------------------------------------------------------------------------
1479     USE DefUtils
1480     IMPLICIT NONE
1481!------------------------------------------------------------------------------
1482     TYPE(Model_t) :: Model
1483     INTEGER :: Perm(:)
1484     REAL(KIND=dp) :: Quant(:), Indicator(2), Gnorm
1485     TYPE( Mesh_t ), POINTER    :: Mesh
1486     TYPE( Element_t ), POINTER :: Edge
1487!------------------------------------------------------------------------------
1488
1489     TYPE(Nodes_t) :: Nodes, EdgeNodes
1490     TYPE(Element_t), POINTER :: Element
1491
1492     INTEGER :: i,k,n,l,t,bc,DIM,DOFs,Pn,En
1493     LOGICAL :: stat, GotIt, Compressible
1494
1495     REAL(KIND=dp) :: Grad(3,3), Grad1(3,3), Stress(3,3), Normal(3), ForceSolved(3), &
1496                      EdgeLength,u, v, w, s, detJ
1497     REAL(KIND=dp) :: SqrtMetric, Metric(3,3), Symb(3,3,3), dSymb(3,3,3,3)
1498
1499     REAL(KIND=dp), ALLOCATABLE :: EdgeBasis(:), dEdgeBasisdx(:,:), Basis(:),dBasisdx(:,:)
1500     REAL(KIND=dp), ALLOCATABLE ::  x(:), y(:), z(:), ExtPressure(:)
1501     REAL(KIND=dp), ALLOCATABLE :: Temperature(:), Tension(:), SlipCoeff(:,:)
1502     REAL(KIND=dp), ALLOCATABLE :: Velocity(:,:), Pressure(:), Force(:,:), NodalViscosity(:)
1503
1504     REAL(KIND=dp) :: Residual(3), ResidualNorm, Viscosity, Slip, Dir(3)
1505
1506     TYPE(Variable_t), POINTER  :: TempSol
1507     TYPE(ValueList_t), POINTER :: Material
1508     TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff
1509!------------------------------------------------------------------------------
1510
1511!    Initialize:
1512!    -----------
1513     Indicator = 0.0d0
1514     Gnorm     = 0.0d0
1515
1516     Metric = 0.0d0
1517     DO i=1,3
1518        Metric(i,i) = 1.0d0
1519     END DO
1520
1521     SELECT CASE( CurrentCoordinateSystem() )
1522        CASE( AxisSymmetric, CylindricSymmetric )
1523           DIM = 3
1524        CASE DEFAULT
1525           DIM = CoordinateSystemDimension()
1526     END SELECT
1527
1528     DOFs = DIM + 1
1529     IF ( CurrentCoordinateSystem() == AxisSymmetric ) DOFs = DOFs-1
1530!
1531!    --------------------------------------------------
1532     Element => Edge % BoundaryInfo % Left
1533
1534     IF ( .NOT. ASSOCIATED( Element ) ) THEN
1535
1536        Element => Edge % BoundaryInfo % Right
1537
1538     ELSE IF ( ANY( Perm( Element % NodeIndexes ) <= 0 ) ) THEN
1539
1540        Element => Edge % BoundaryInfo % Right
1541
1542     END IF
1543
1544     IF ( .NOT. ASSOCIATED( Element ) ) RETURN
1545     IF ( ANY( Perm( Element % NodeIndexes ) <= 0 ) ) RETURN
1546
1547     En = Edge % TYPE % NumberOfNodes
1548     Pn = Element % TYPE % NumberOfNodes
1549
1550     ALLOCATE( EdgeNodes % x(En), EdgeNodes % y(En), EdgeNodes % z(En) )
1551
1552     EdgeNodes % x = Mesh % Nodes % x(Edge % NodeIndexes)
1553     EdgeNodes % y = Mesh % Nodes % y(Edge % NodeIndexes)
1554     EdgeNodes % z = Mesh % Nodes % z(Edge % NodeIndexes)
1555
1556     ALLOCATE( Nodes % x(Pn), Nodes % y(Pn), Nodes % z(Pn) )
1557
1558     Nodes % x = Mesh % Nodes % x(Element % NodeIndexes)
1559     Nodes % y = Mesh % Nodes % y(Element % NodeIndexes)
1560     Nodes % z = Mesh % Nodes % z(Element % NodeIndexes)
1561
1562     ALLOCATE( EdgeBasis(En), dEdgeBasisdx(En,3), Basis(Pn), dBasisdx(Pn,3), &
1563      x(En), y(En), z(En), ExtPressure(En), Temperature(Pn), Tension(En),    &
1564      SlipCoeff(3,En), Velocity(3,Pn), Pressure(Pn), Force(3,En), NodalViscosity(En) )
1565
1566     DO l = 1,En
1567       DO k = 1,Pn
1568          IF ( Edge % NodeIndexes(l) == Element % NodeIndexes(k) ) THEN
1569             x(l) = Element % TYPE % NodeU(k)
1570             y(l) = Element % TYPE % NodeV(k)
1571             z(l) = Element % TYPE % NodeW(k)
1572             EXIT
1573          END IF
1574       END DO
1575     END DO
1576!
1577!    Integrate square of residual over boundary element:
1578!    ---------------------------------------------------
1579
1580     Indicator    = 0.0d0
1581     EdgeLength   = 0.0d0
1582     ResidualNorm = 0.0d0
1583
1584     DO bc=1,Model % NumberOfBCs
1585        IF ( Edge % BoundaryInfo % Constraint /= Model % BCs(bc) % Tag ) CYCLE
1586
1587!       IF ( .NOT. ListGetLogical( Model % BCs(bc) % Values, &
1588!                 'Flow Force BC', gotIt ) ) CYCLE
1589!
1590!       Get material parameters:
1591!       ------------------------
1592
1593        k = ListGetInteger(Model % Bodies(Element % BodyId) % Values,'Material', &
1594                     minv=1, maxv=Model % NumberOfMaterials )
1595        Material => Model % Materials(k) % Values
1596
1597        NodalViscosity(1:En) = ListGetReal( Material, &
1598                 'Viscosity', En, Edge % NodeIndexes, GotIt )
1599
1600        Compressible = .FALSE.
1601        IF ( ListGetString( Material, 'Compressibility Model', GotIt ) == &
1602               'perfect gas equation 1' ) Compressible = .TRUE.
1603
1604
1605!       Given traction:
1606!       ---------------
1607        Force = 0.0d0
1608
1609        Force(1,1:En) = ListGetReal( Model % BCs(bc) % Values, &
1610            'Pressure 1', En, Edge % NodeIndexes, GotIt )
1611
1612        Force(2,1:En) = ListGetReal( Model % BCs(bc) % Values, &
1613            'Pressure 2', En, Edge % NodeIndexes, GotIt )
1614
1615        Force(3,1:En) = ListGetReal( Model % BCs(bc) % Values, &
1616            'Pressure 3', En, Edge % NodeIndexes, GotIt )
1617
1618!
1619!       Force in normal direction:
1620!       ---------------------------
1621        ExtPressure(1:En) = ListGetReal( Model % BCs(bc) % Values, &
1622          'External Pressure', En, Edge % NodeIndexes, GotIt )
1623
1624!
1625!       Slip BC condition:
1626!       ------------------
1627        SlipCoeff = 0.0d0
1628        SlipCoeff(1,1:En) =  ListGetReal( Model % BCs(bc) % Values, &
1629             'Slip Coefficient 1',En,Edge % NodeIndexes,GotIt )
1630
1631        SlipCoeff(2,1:En) =  ListGetReal( Model % BCs(bc) % Values, &
1632             'Slip Coefficient 2',En,Edge % NodeIndexes,GotIt )
1633
1634        SlipCoeff(3,1:En) =  ListGetReal( Model % BCs(bc) % Values, &
1635             'Slip Coefficient 3',En,Edge % NodeIndexes,GotIt )
1636
1637!
1638!       Surface tension induced by temperature gradient (or otherwise):
1639!       ---------------------------------------------------------------
1640        TempSol => VariableGet( Mesh % Variables, 'Temperature', .TRUE. )
1641
1642        IF ( ASSOCIATED( TempSol ) ) THEN
1643          Tension(1:En) = ListGetReal( Model % BCs(bc) % Values, &
1644           'Surface Tension Expansion Coefficient',En,Edge % NodeIndexes,gotIt )
1645
1646           IF ( gotIt ) THEN
1647              DO n=1,En
1648                 k = TempSol % Perm( Edge % NodeIndexes(n) )
1649                 IF (k>0) Tension(n) = 1.0d0 - Tension(n) * TempSol % Values(k)
1650              END DO
1651
1652              Tension(1:En) = Tension(1:En) * ListGetReal( &
1653                 Model % BCs(bc) % Values,'Surface Tension Coefficient', &
1654                               En, Edge % NodeIndexes )
1655           ELSE
1656              Tension(1:En) = ListGetReal( &
1657                  Model % BCs(bc) % Values,'Surface Tension Coefficient', &
1658                         En, Edge % NodeIndexes,gotIt )
1659           END IF
1660        ELSE
1661           Tension(1:En) = ListGetReal( &
1662               Model % BCs(bc) % Values,'Surface Tension Coefficient', &
1663                      En, Edge % NodeIndexes,gotIt )
1664        END IF
1665
1666!
1667!       If dirichlet BC for velocity in any direction given,
1668!       nullify force in that directon:
1669!       ------------------------------------------------------------------
1670        Dir = 1
1671        s = ListGetConstReal( Model % BCs(bc) % Values, 'Velocity 1', GotIt )
1672        IF ( GotIt ) Dir(1) = 0
1673
1674        s = ListGetConstReal( Model % BCs(bc) % Values, 'Velocity 2', GotIt )
1675        IF ( GotIt ) Dir(2) = 0
1676
1677        s = ListGetConstReal( Model % BCs(bc) % Values, 'Velocity 3', GotIt )
1678        IF ( GotIt ) Dir(3) = 0
1679
1680!
1681!       Elementwise nodal solution:
1682!       ---------------------------
1683        Velocity = 0.0d0
1684        DO k=1,DOFs-1
1685           Velocity(k,1:Pn) = Quant(DOFs*Perm(Element % NodeIndexes)-DOFs + k)
1686        END DO
1687        Pressure(1:Pn) = Quant( DOFs*Perm(Element % NodeIndexes) )
1688
1689!       do the integration:
1690!       -------------------
1691        EdgeLength   = 0.0d0
1692        ResidualNorm = 0.0d0
1693
1694        IntegStuff = GaussPoints( Edge )
1695
1696        DO t=1,IntegStuff % n
1697           u = IntegStuff % u(t)
1698           v = IntegStuff % v(t)
1699           w = IntegStuff % w(t)
1700
1701           stat = ElementInfo( Edge, EdgeNodes, u, v, w, detJ, &
1702               EdgeBasis, dEdgeBasisdx )
1703
1704           IF ( CurrentCoordinateSystem() == Cartesian ) THEN
1705              s = IntegStuff % s(t) * detJ
1706           ELSE
1707              u = SUM( EdgeBasis(1:En) * EdgeNodes % x(1:En) )
1708              v = SUM( EdgeBasis(1:En) * EdgeNodes % y(1:En) )
1709              w = SUM( EdgeBasis(1:En) * EdgeNodes % z(1:En) )
1710
1711              CALL CoordinateSystemInfo( Metric, SqrtMetric, &
1712                          Symb, dSymb, u, v, w )
1713
1714              s = IntegStuff % s(t) * detJ * SqrtMetric
1715           END IF
1716
1717           Normal = NormalVector( Edge, EdgeNodes, u, v, .TRUE. )
1718
1719           u = SUM( EdgeBasis(1:En) * x(1:En) )
1720           v = SUM( EdgeBasis(1:En) * y(1:En) )
1721           w = SUM( EdgeBasis(1:En) * z(1:En) )
1722
1723           stat = ElementInfo( Element, Nodes, u, v, w, detJ, &
1724              Basis, dBasisdx )
1725
1726           Viscosity = SUM( NodalViscosity(1:En) * EdgeBasis(1:En) )
1727
1728           Residual = 0.0d0
1729!
1730!          Given force at the integration point:
1731!          -------------------------------------
1732           Residual = Residual + MATMUL( Force(:,1:En), EdgeBasis(1:En) ) - &
1733                 SUM( ExtPressure(1:En) * EdgeBasis(1:En) ) * Normal
1734
1735!
1736!          Slip velocity BC:
1737!          -----------------
1738           DO i=1,DIM
1739              Slip = SUM( SlipCoeff(i,1:En) * EdgeBasis(1:En) )
1740              Residual(i) = Residual(i) - &
1741                   Slip * SUM( Velocity(i,1:Pn) * Basis(1:Pn) )
1742           END DO
1743
1744!
1745!          Tangential tension force:
1746!          -------------------------
1747           DO i=1,DIM
1748              Residual(i) = Residual(i) + &
1749                   SUM( dEdgeBasisdx(1:En,i) * Tension(1:En) )
1750           END DO
1751
1752!
1753!          Force given by the computed solution:
1754!          -------------------------------------
1755!
1756!          Stress tensor on the boundary:
1757!          ------------------------------
1758           Grad = MATMUL( Velocity(:,1:Pn), dBasisdx(1:Pn,:) )
1759
1760           IF ( CurrentCoordinateSystem() /= Cartesian ) THEN
1761              Grad1 = Grad
1762              DO i=1,DIM
1763                 DO k=1,DIM
1764                    DO l=1,DIM
1765                       Grad1(i,k) = Grad1(i,k) - &
1766                          Symb(k,l,i) * SUM ( Velocity(l,1:Pn) * Basis(1:Pn) )
1767                    END DO
1768                 END DO
1769              END DO
1770
1771              Grad = 0.0d0
1772              DO i=1,DIM
1773                 DO k=1,DIM
1774                    DO l=1,DIM
1775                       Grad(i,k) = Grad(i,k) + Metric(k,l) * Grad1(i,l)
1776                    END DO
1777                 END DO
1778              END DO
1779           END IF
1780
1781           Stress = Viscosity * ( Grad + TRANSPOSE(Grad) )
1782           Stress = Stress - Metric * SUM( Pressure(1:Pn) * Basis(1:Pn) )
1783
1784           IF ( Compressible ) THEN
1785              IF ( CurrentCoordinateSystem() == Cartesian ) THEN
1786                 DO i=1,DIM
1787                    DO k=1,DIM
1788                       Stress(i,i) = Stress(i,i) - &
1789                           (2.0d0/3.0d0) * Viscosity * Grad(k,k)
1790                    END DO
1791                 END DO
1792              ELSE
1793                 DO i=1,DIM
1794                    DO k=1,DIM
1795                       DO l=1,DIM
1796                          Stress(i,k) = Stress(i,k) - &
1797                             Metric(i,k) * (2.0d0/3.0d0) * Viscosity * Grad(l,l)
1798                       END DO
1799                    END DO
1800                 END DO
1801              END IF
1802           END IF
1803
1804           ForceSolved = MATMUL(Stress,Normal)
1805           Residual = Residual - ForceSolved * Dir
1806
1807           EdgeLength = EdgeLength + s
1808
1809           IF ( CurrentCoordinateSystem() == Cartesian ) THEN
1810              Gnorm = Gnorm + s * SUM( ForceSolved**2 )
1811              ResidualNorm = ResidualNorm + s * SUM( Residual(1:DIM) ** 2 )
1812           ELSE
1813              CALL InvertMatrix( Metric,3 )
1814              DO i=1,DIM
1815                 DO k=1,DIM
1816                    ResidualNorm = ResidualNorm + &
1817                            s * Metric(i,k) * Residual(i) * Residual(k)
1818                    Gnorm = GNorm + s * Metric(i,k) * &
1819                                        ForceSolved(i) * ForceSolved(k)
1820                 END DO
1821              END DO
1822           END IF
1823        END DO
1824        EXIT
1825     END DO
1826
1827     IF ( CoordinateSystemDimension() == 3 ) EdgeLength = SQRT(EdgeLength)
1828     Indicator = EdgeLength * ResidualNorm
1829
1830     DEALLOCATE( Nodes % x, Nodes % y, Nodes % z)
1831     DEALLOCATE( EdgeNodes % x, EdgeNodes % y, EdgeNodes % z)
1832
1833     DEALLOCATE( EdgeBasis, dEdgeBasisdx, Basis, dBasisdx, x, y, z,   &
1834      ExtPressure, Temperature, Tension,SlipCoeff, Velocity, Pressure, &
1835      Force, NodalViscosity )
1836!------------------------------------------------------------------------------
1837  END FUNCTION FlowBoundaryResidual
1838!------------------------------------------------------------------------------
1839
1840
1841!------------------------------------------------------------------------------
1842!> Compute the residual of the Navier-Stokes equation for the edge elements.
1843!------------------------------------------------------------------------------
1844  FUNCTION FlowEdgeResidual( Model,Edge,Mesh,Quant,Perm ) RESULT( Indicator )
1845!------------------------------------------------------------------------------
1846     USE DefUtils
1847     IMPLICIT NONE
1848
1849     TYPE(Model_t) :: Model
1850     INTEGER :: Perm(:)
1851     REAL(KIND=dp) :: Quant(:), Indicator(2)
1852     TYPE( Mesh_t ), POINTER    :: Mesh
1853     TYPE( Element_t ), POINTER :: Edge
1854!------------------------------------------------------------------------------
1855
1856     TYPE(Nodes_t) :: Nodes, EdgeNodes
1857     TYPE(Element_t), POINTER :: Element
1858
1859     INTEGER :: i,j,k,l,n,t,DIM,DOFs,En,Pn
1860     LOGICAL :: stat, GotIt
1861
1862     REAL(KIND=dp) :: SqrtMetric, Metric(3,3), Symb(3,3,3), dSymb(3,3,3,3)
1863     REAL(KIND=dp) :: Stress(3,3,2), Jump(3), Viscosity
1864
1865     REAL(KIND=dp) :: Grad(3,3), Grad1(3,3), Normal(3)
1866
1867     REAL(KIND=dp), ALLOCATABLE :: NodalViscosity(:), x(:), y(:), z(:), &
1868               EdgeBasis(:), Basis(:), dBasisdx(:,:)
1869     REAL(KIND=dp), ALLOCATABLE :: Velocity(:,:), Pressure(:)
1870
1871     REAL(KIND=dp) :: ResidualNorm, EdgeLength, u, v, w, s, detJ
1872
1873     TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff
1874!------------------------------------------------------------------------------
1875
1876!    Initialize:
1877!    -----------
1878     SELECT CASE( CurrentCoordinateSystem() )
1879        CASE( AxisSymmetric, CylindricSymmetric )
1880           DIM = 3
1881        CASE DEFAULT
1882           DIM = CoordinateSystemDimension()
1883     END SELECT
1884
1885     DOFs = DIM + 1
1886     IF ( CurrentCoordinateSystem() == AxisSymmetric ) DOFs = DOFs - 1
1887
1888     Metric = 0.0d0
1889     DO i = 1,3
1890        Metric(i,i) = 1.0d0
1891     END DO
1892
1893     Grad = 0.0d0
1894!
1895!    ---------------------------------------------
1896
1897     Element => Edge % BoundaryInfo % Left
1898     n = Element % TYPE % NumberOfNodes
1899
1900     Element => Edge % BoundaryInfo % Right
1901     n = MAX( n, Element % TYPE % NumberOfNodes )
1902
1903     ALLOCATE( Nodes % x(n), Nodes % y(n), Nodes % z(n) )
1904
1905     En = Edge % TYPE % NumberOfNodes
1906     ALLOCATE( EdgeNodes % x(En), EdgeNodes % y(En), EdgeNodes % z(En) )
1907
1908     EdgeNodes % x = Mesh % Nodes % x(Edge % NodeIndexes)
1909     EdgeNodes % y = Mesh % Nodes % y(Edge % NodeIndexes)
1910     EdgeNodes % z = Mesh % Nodes % z(Edge % NodeIndexes)
1911
1912     ALLOCATE( NodalViscosity(En), x(En), y(En), z(En), EdgeBasis(En), &
1913           Basis(n), dBasisdx(n,3), Velocity(3,n), Pressure(n) )
1914
1915!    Integrate square of jump over edge:
1916!    ------------------------------------
1917     ResidualNorm = 0.0d0
1918     EdgeLength   = 0.0d0
1919     Indicator    = 0.0d0
1920
1921     IntegStuff = GaussPoints( Edge )
1922
1923     DO t=1,IntegStuff % n
1924
1925        u = IntegStuff % u(t)
1926        v = IntegStuff % v(t)
1927        w = IntegStuff % w(t)
1928
1929        stat = ElementInfo( Edge, EdgeNodes, u, v, w, detJ, &
1930             EdgeBasis, dBasisdx )
1931
1932        Normal = NormalVector( Edge, EdgeNodes, u, v, .FALSE. )
1933
1934        IF ( CurrentCoordinateSystem() == Cartesian ) THEN
1935           s = IntegStuff % s(t) * detJ
1936        ELSE
1937           u = SUM( EdgeBasis(1:En) * EdgeNodes % x(1:En) )
1938           v = SUM( EdgeBasis(1:En) * EdgeNodes % y(1:En) )
1939           w = SUM( EdgeBasis(1:En) * EdgeNodes % z(1:En) )
1940
1941           CALL CoordinateSystemInfo( Metric, SqrtMetric, &
1942                       Symb, dSymb, u, v, w )
1943           s = IntegStuff % s(t) * detJ * SqrtMetric
1944        END IF
1945
1946        Stress = 0.0d0
1947        DO i = 1,2
1948           IF ( i==1 ) THEN
1949              Element => Edge % BoundaryInfo % Left
1950           ELSE
1951              Element => Edge % BoundaryInfo % Right
1952           END IF
1953
1954           IF ( ANY( Perm( Element % NodeIndexes ) <= 0 ) ) CYCLE
1955
1956           Pn = Element % TYPE % NumberOfNodes
1957           Nodes % x(1:Pn) = Mesh % Nodes % x(Element % NodeIndexes)
1958           Nodes % y(1:Pn) = Mesh % Nodes % y(Element % NodeIndexes)
1959           Nodes % z(1:Pn) = Mesh % Nodes % z(Element % NodeIndexes)
1960
1961           DO j = 1,En
1962              DO k = 1,Pn
1963                 IF ( Edge % NodeIndexes(j) == Element % NodeIndexes(k) ) THEN
1964                    x(j) = Element % TYPE % NodeU(k)
1965                    y(j) = Element % TYPE % NodeV(k)
1966                    z(j) = Element % TYPE % NodeW(k)
1967                    EXIT
1968                 END IF
1969              END DO
1970           END DO
1971
1972           u = SUM( EdgeBasis(1:En) * x(1:En) )
1973           v = SUM( EdgeBasis(1:En) * y(1:En) )
1974           w = SUM( EdgeBasis(1:En) * z(1:En) )
1975
1976           stat = ElementInfo( Element, Nodes, u, v, w, detJ, &
1977               Basis, dBasisdx )
1978
1979           k = ListGetInteger( Model % Bodies( Element % BodyId) % Values, 'Material', &
1980                            minv=1, maxv=Model % NumberOfMaterials )
1981
1982           NodalViscosity(1:En) = ListGetReal( &
1983               Model % Materials(k) % Values, 'Viscosity', &
1984                    En, Edge % NodeIndexes, GotIt )
1985
1986           Viscosity = SUM( NodalViscosity(1:En) * EdgeBasis(1:En) )
1987!
1988!          Elementwise nodal solution:
1989!          ---------------------------
1990           Velocity = 0.0d0
1991           DO k=1,DOFs-1
1992              Velocity(k,1:Pn) = Quant(DOFs*Perm(Element % NodeIndexes)-DOFs+k)
1993           END DO
1994           Pressure(1:Pn) = Quant( DOFs*Perm(Element % NodeIndexes) )
1995!
1996!          Stress tensor on the edge:
1997!          --------------------------
1998           Grad = MATMUL( Velocity(:,1:Pn), dBasisdx(1:Pn,:) )
1999
2000           IF ( CurrentCoordinateSystem() /= Cartesian ) THEN
2001              Grad1 = Grad
2002              DO j=1,DIM
2003                 DO k=1,DIM
2004                    DO l=1,DIM
2005                       Grad1(j,k) = Grad1(j,k) - &
2006                          Symb(k,l,j) * SUM ( Velocity(l,1:Pn) * Basis(1:Pn) )
2007                    END DO
2008                 END DO
2009              END DO
2010
2011              Grad = 0.0d0
2012              DO j=1,DIM
2013                 DO k=1,DIM
2014                    DO l=1,DIM
2015                       Grad(j,k) = Grad(j,k) + Metric(k,l) * Grad1(j,l)
2016                    END DO
2017                 END DO
2018              END DO
2019           END IF
2020
2021           Stress(:,:,i) = Viscosity * ( Grad + TRANSPOSE(Grad) )
2022
2023           IF ( CurrentCoordinateSystem() == Cartesian ) THEN
2024              DO j=1,DIM
2025                 Stress(j,j,i) = Stress(j,j,i) - SUM( Pressure(1:Pn) * Basis(1:Pn))
2026                 DO k=1,DIM
2027                    Stress(j,j,i) = Stress(j,j,i) - (2.0d0/3.0d0)*Viscosity*Grad(k,k)
2028                 END DO
2029              END DO
2030           ELSE
2031              DO j=1,DIM
2032                 DO k=1,DIM
2033                    Stress(j,k,i) = Stress(j,k,i) - &
2034                           Metric(j,k) * SUM( Pressure(1:Pn) * Basis(1:Pn) )
2035
2036                    DO l=1,DIM
2037                       Stress(j,k,i) = Stress(j,k,i) - &
2038                           Metric(j,k) * (2.0d0/3.0d0) * Viscosity * Grad(l,l)
2039                    END DO
2040                 END DO
2041              END DO
2042           END IF
2043
2044        END DO
2045
2046        EdgeLength = EdgeLength + s
2047
2048        Jump = MATMUL( ( Stress(:,:,1) - Stress(:,:,2)), Normal )
2049
2050        IF ( CurrentCoordinateSystem() == Cartesian ) THEN
2051           ResidualNorm = ResidualNorm + s * SUM( Jump(1:DIM) ** 2 )
2052        ELSE
2053           CALL InvertMatrix( Metric,3 )
2054           DO i=1,DIM
2055              DO j=1,DIM
2056                 ResidualNorm = ResidualNorm + s*Metric(i,j)*Jump(i)*Jump(j)
2057              END DO
2058           END DO
2059        END IF
2060     END DO
2061
2062     Indicator = EdgeLength * ResidualNorm
2063
2064     DEALLOCATE( Nodes % x, Nodes % y, Nodes % z)
2065     DEALLOCATE( EdgeNodes % x, EdgeNodes % y, EdgeNodes % z)
2066
2067     DEALLOCATE( NodalViscosity, x, y, z, EdgeBasis, &
2068           Basis, dBasisdx, Velocity, Pressure )
2069!------------------------------------------------------------------------------
2070  END FUNCTION FlowEdgeResidual
2071!------------------------------------------------------------------------------
2072
2073
2074!------------------------------------------------------------------------------
2075!> Compute the residual of the Navier-Stokes equation for the bulk elements.
2076!------------------------------------------------------------------------------
2077   FUNCTION FlowInsideResidual( Model, Element,  &
2078          Mesh, Quant, Perm, Fnorm ) RESULT( Indicator )
2079!------------------------------------------------------------------------------
2080     USE DefUtils
2081!------------------------------------------------------------------------------
2082     IMPLICIT NONE
2083!------------------------------------------------------------------------------
2084     TYPE(Model_t) :: Model
2085     INTEGER :: Perm(:)
2086     REAL(KIND=dp) :: Quant(:), Indicator(2), FNorm
2087     TYPE( Mesh_t ), POINTER    :: Mesh
2088     TYPE( Element_t ), POINTER :: Element
2089!------------------------------------------------------------------------------
2090
2091     TYPE(Nodes_t) :: Nodes
2092
2093     INTEGER :: i,j,k,l,m,n,t,DIM,DOFs
2094
2095     LOGICAL :: stat, GotIt, Compressible, Convect
2096
2097     TYPE( Variable_t ), POINTER :: Var
2098
2099     REAL(KIND=dp) :: SqrtMetric, Metric(3,3), Symb(3,3,3), dSymb(3,3,3,3)
2100     REAL(KIND=dp) :: Density, Viscosity,u, v, w, s, detJ
2101     REAL(KIND=dp) :: Residual(4), ResidualNorm, Area, ReferencePressure, dt
2102
2103     REAL(KIND=dp), ALLOCATABLE :: NodalViscosity(:), NodalDensity(:), &
2104            Basis(:),  dBasisdx(:,:), ddBasisddx(:,:,:)
2105     REAL(KIND=dp),ALLOCATABLE :: Velocity(:,:), Pressure(:)
2106     REAL(KIND=dp),ALLOCATABLE :: PrevVelo(:,:), PrevPres(:)
2107     REAL(KIND=dp),ALLOCATABLE :: Temperature(:), NodalForce(:,:)
2108     REAL(KIND=dp),ALLOCATABLE :: HeatCapacity(:), ReferenceTemperature(:), &
2109                      HeatExpansionCoeff(:)
2110
2111     REAL(KIND=dp) :: SpecificHeatRatio
2112
2113     REAL(KIND=dp), POINTER :: Gravity(:,:)
2114
2115     TYPE(ValueList_t), POINTER :: Material
2116
2117     TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff
2118!------------------------------------------------------------------------------
2119
2120!    Initialize:
2121!    -----------
2122     Indicator = 0.0d0
2123     FNorm = 0.0d0
2124
2125     IF ( ANY( Perm( Element % NodeIndexes ) <= 0 ) ) RETURN
2126
2127     Metric = 0.0d0
2128     DO i=1,3
2129        Metric(i,i) = 1.0d0
2130     END DO
2131
2132     SELECT CASE( CurrentCoordinateSystem() )
2133        CASE( AxisSymmetric, CylindricSymmetric )
2134           DIM = 3
2135        CASE DEFAULT
2136           DIM = CoordinateSystemDimension()
2137     END SELECT
2138
2139     DOFs = DIM + 1
2140     IF ( CurrentCoordinateSystem() == AxisSymmetric ) DOFs = DOFs-1
2141!
2142!    Element nodal points:
2143!    ---------------------
2144     n = Element % TYPE % NumberOfNodes
2145
2146     ALLOCATE( Nodes % x(n), Nodes % y(n), Nodes % z(n) )
2147     Nodes % x = Mesh % Nodes % x(Element % NodeIndexes)
2148     Nodes % y = Mesh % Nodes % y(Element % NodeIndexes)
2149     Nodes % z = Mesh % Nodes % z(Element % NodeIndexes)
2150
2151     ALLOCATE( NodalViscosity(n), NodalDensity(n), Basis(n), dBasisdx(n,3), &
2152        ddBasisddx(n,3,3), Velocity(3,n), Pressure(n), PrevVelo(3,n),  &
2153        PrevPres(n), Temperature(n), NodalForce(4,n), HeatCapacity(n), &
2154        ReferenceTemperature(n), HeatExpansionCoeff(n) )
2155
2156!
2157!    Material parameters: density, viscosity, etc.
2158!    ----------------------------------------------
2159     k = ListGetInteger( Model % Bodies(Element % BodyId) % Values, 'Material', &
2160                  minv=1, maxv=Model % NumberOfMaterials )
2161
2162     Material => Model % Materials(k) % Values
2163
2164     NodalDensity(1:n) = ListGetReal( &
2165         Material, 'Density', n, Element % NodeIndexes, GotIt )
2166
2167     NodalViscosity(1:n) = ListGetReal( &
2168         Material, 'Viscosity', n, Element % NodeIndexes, GotIt )
2169
2170     k = ListGetInteger( Model % Bodies(Element % BodyId) % Values,'Equation', &
2171                      minv=1, maxv=Model % NumberOfEquations   )
2172
2173     Convect = ListGetLogical( Model % Equations(k) % Values, &
2174                   'NS Convect', GotIt )
2175     IF ( .NOT. GotIt ) Convect = .TRUE.
2176!
2177!    Elementwise nodal solution:
2178!    ---------------------------
2179     Velocity = 0.0d0
2180     DO k=1,DOFs-1
2181        Velocity(k,1:n) = Quant( DOFs*Perm(Element % NodeIndexes)-DOFs+k )
2182     END DO
2183     Pressure(1:n) = Quant( DOFs*Perm(Element % NodeIndexes) )
2184
2185!
2186!    Check for time dep.
2187!    -------------------
2188     PrevPres(1:n)     = Pressure(1:n)
2189     PrevVelo(1:3,1:n) = Velocity(1:3,1:n)
2190
2191     dt = Model % Solver % dt
2192
2193     IF ( ListGetString( Model % Simulation, 'Simulation Type') == 'transient' ) THEN
2194        Var => VariableGet( Model % Variables, 'Flow Solution', .TRUE. )
2195
2196        PrevVelo = 0.0d0
2197        DO k=1,DOFs-1
2198           PrevVelo(k,1:n) = &
2199              Var % PrevValues(DOFs*Var % Perm(Element % NodeIndexes)-DOFs+k,1)
2200        END DO
2201        PrevPres(1:n)=Var % PrevValues(DOFs*Var % Perm(Element % NodeIndexes),1)
2202     END IF
2203
2204
2205!
2206!    Check for compressible flow equations:
2207!    --------------------------------------
2208     Compressible = .FALSE.
2209
2210     IF (  ListGetString( Material, 'Compressibility Model', GotIt ) == &
2211                      'perfect gas equation 1' ) THEN
2212
2213        Compressible = .TRUE.
2214
2215        Var => VariableGet( Mesh % Variables, 'Temperature', .TRUE. )
2216        IF ( ASSOCIATED( Var ) ) THEN
2217           Temperature(1:n) = &
2218               Var % Values( Var % Perm(Element % NodeIndexes) )
2219        ELSE
2220           Temperature(1:n) = ListGetReal( Material, &
2221               'Reference Temperature',n,Element % NodeIndexes )
2222        END IF
2223
2224        SpecificHeatRatio = ListGetConstReal( Material, &
2225                  'Specific Heat Ratio' )
2226
2227        ReferencePressure = ListGetConstReal( Material, &
2228                   'Reference Pressure' )
2229
2230        HeatCapacity(1:n) = ListGetReal( Material, &
2231                      'Heat Capacity',n,Element % NodeIndexes )
2232
2233        NodalDensity(1:n) =  (Pressure(1:n) + ReferencePressure) * SpecificHeatRatio / &
2234              ( (SpecificHeatRatio - 1) * HeatCapacity(1:n) * Temperature(1:n) )
2235     END IF
2236!
2237!    Body Forces:
2238!    ------------
2239!
2240     k = ListGetInteger( Model % Bodies(Element % BodyId) % Values, &
2241       'Body Force', GotIt, 1, Model % NumberOfBodyForces )
2242
2243     NodalForce = 0.0d0
2244
2245     IF ( GotIt .AND. k > 0  ) THEN
2246!
2247!       Boussinesq approximation of heat expansion for
2248!       incompressible flow equations:
2249!
2250!       Density for the force term equals to
2251!
2252!       \rho = rho_0 (1-\beta(T-T_0)),
2253!
2254!       where \beta is the  heat expansion  coefficient,
2255!       T temperature and \rho_0 and T_0 correspond to
2256!       stress free state. Otherwise density is assumed
2257!       constant.
2258!       ----------------------------------------------
2259        IF (ListGetLogical(Model % BodyForces(k) % Values,'Boussinesq',GotIt)) THEN
2260
2261           Var => VariableGet( Mesh % Variables, 'Temperature', .TRUE. )
2262           IF ( ASSOCIATED( Var ) ) THEN
2263              Temperature(1:n) = &
2264                  Var % Values( Var % Perm(Element % NodeIndexes) )
2265
2266              HeatExpansionCoeff(1:n) = ListGetReal( Material, &
2267                 'Heat Expansion Coefficient',n,Element % NodeIndexes )
2268
2269              ReferenceTemperature(1:n) = ListGetReal( Material, &
2270                 'Reference Temperature',n,Element % NodeIndexes )
2271
2272              Gravity => ListGetConstRealArray( Model % Constants, &
2273                             'Gravity' )
2274
2275              k = ListGetInteger( Model % Bodies(Element % BodyId) % Values,'Equation', &
2276                        minv=1, maxv=Model % NumberOfEquations )
2277
2278              IF ( ListGetLogical( Model % Equations(k) % Values, &
2279                            'Hydrostatic Pressure', GotIt) ) THEN
2280                 DO i=1,DIM
2281                    NodalForce(i,1:n) = ( 1 - HeatExpansionCoeff(1:n) * &
2282                       ( Temperature(1:n) - ReferenceTemperature(1:n) ) ) * &
2283                            Gravity(i,1) * Gravity(4,1)
2284                 END DO
2285              ELSE
2286                 DO i=1,DIM
2287                    NodalForce(i,1:n) = ( -HeatExpansionCoeff(1:n) * &
2288                       ( Temperature(1:n) - ReferenceTemperature(1:n) ) ) * &
2289                            Gravity(i,1) * Gravity(4,1)
2290                 END DO
2291              END IF
2292           END IF
2293        END IF
2294
2295!
2296!       Given external force:
2297!       ---------------------
2298        NodalForce(1,1:n) = NodalForce(1,1:n) + ListGetReal( &
2299             Model % BodyForces(k) % Values, 'Flow BodyForce 1', &
2300                  n, Element % NodeIndexes, GotIt )
2301
2302        NodalForce(2,1:n) = NodalForce(2,1:n) + ListGetReal( &
2303             Model % BodyForces(k) % Values, 'Flow BodyForce 2', &
2304                  n, Element % NodeIndexes, GotIt )
2305
2306        NodalForce(3,1:n) = NodalForce(3,1:n) + ListGetReal( &
2307             Model % BodyForces(k) % Values, 'Flow BodyForce 3', &
2308                  n, Element % NodeIndexes, GotIt )
2309     END IF
2310!
2311!    Integrate square of residual over element:
2312!    ------------------------------------------
2313     ResidualNorm = 0.0d0
2314     Area = 0.0d0
2315
2316     IntegStuff = GaussPoints( Element )
2317
2318     DO t=1,IntegStuff % n
2319        u = IntegStuff % u(t)
2320        v = IntegStuff % v(t)
2321        w = IntegStuff % w(t)
2322
2323        stat = ElementInfo( Element, Nodes, u, v, w, detJ, &
2324            Basis, dBasisdx, ddBasisddx, .TRUE. )
2325
2326        IF ( CurrentCoordinateSystem() == Cartesian ) THEN
2327           s = IntegStuff % s(t) * detJ
2328        ELSE
2329           u = SUM( Basis(1:n) * Nodes % x(1:n) )
2330           v = SUM( Basis(1:n) * Nodes % y(1:n) )
2331           w = SUM( Basis(1:n) * Nodes % z(1:n) )
2332
2333           CALL CoordinateSystemInfo( Metric, SqrtMetric, &
2334                      Symb, dSymb, u, v, w )
2335           s = IntegStuff % s(t) * detJ * SqrtMetric
2336        END IF
2337
2338        Density   = SUM( NodalDensity(1:n)   * Basis(1:n) )
2339        Viscosity = SUM( NodalViscosity(1:n) * Basis(1:n) )
2340!
2341!       Residual of the navier-stokes equations:
2342!
2343!       or more generally:
2344!
2345!       ----------------------------------------------------------
2346!
2347        Residual = 0.0d0
2348        DO i=1,DIM
2349!
2350!          given force:
2351!          -------------
2352           Residual(i) = -Density * SUM( NodalForce(i,1:n) * Basis(1:n) )
2353
2354           IF ( CurrentCoordinateSystem() == Cartesian ) THEN
2355!             + grad(p):
2356!             ----------
2357              Residual(i) = Residual(i) + SUM( Pressure(1:n) * dBasisdx(1:n,i) )
2358
2359              DO j=1,DIM
2360!
2361!                - 2 ( \mu \epsilon^{ij} )_{,j}:
2362!                -------------------------------
2363                 Residual(i) = Residual(i) - Viscosity * &
2364                     SUM( Velocity(i,1:n) * ddBasisddx(1:n,j,j) )
2365
2366                 Residual(i) = Residual(i) - &
2367                      SUM( NodalViscosity(1:n) * dBasisdx(1:n,j) ) * &
2368                          SUM( Velocity(i,1:n) * dBasisdx(1:n,j) )
2369
2370                  Residual(i) = Residual(i) - Viscosity * &
2371                      SUM( Velocity(j,1:n) * ddBasisddx(1:n,i,j) )
2372
2373                  Residual(i) = Residual(i) - &
2374                      SUM( NodalViscosity(1:n) * dBasisdx(1:n,j) ) * &
2375                          SUM( Velocity(j,1:n) * dBasisdx(1:n,i) )
2376
2377                  IF ( Compressible ) THEN
2378!
2379!                    + (2/3) grad(\mu div(u)):
2380!                    -------------------------
2381                     Residual(i) = Residual(i) + &
2382                        Viscosity * ( 2.0d0 / 3.0d0 ) * &
2383                           SUM( Velocity(j,1:n) * ddBasisddx(1:n,j,i) )
2384
2385                     Residual(i) = Residual(i) + &
2386                         SUM( NodalViscosity(1:n) * dBasisdx(1:n,i) ) * &
2387                             SUM( Velocity(j,1:n) * dBasisdx(1:n,j) )
2388
2389                  END IF
2390              END DO
2391
2392              IF ( Convect ) THEN
2393!
2394!                + \rho * (@u/@t + u.grad(u)):
2395!                -----------------------------
2396                 Residual(i) = Residual(i) + Density *  &
2397                     SUM((Velocity(i,1:n)-PrevVelo(i,1:n))*Basis(1:n)) / dt
2398
2399                 DO j=1,DIM
2400                    Residual(i) = Residual(i) + &
2401                        Density * SUM( Velocity(j,1:n) * Basis(1:n) ) * &
2402                            SUM( Velocity(i,1:n) * dBasisdx(1:n,j) )
2403                 END DO
2404              END IF
2405           ELSE
2406!             + g^{ij}p_{,j}:
2407!             ---------------
2408              DO j=1,DIM
2409                 Residual(i) = Residual(i) + Metric(i,j) * &
2410                      SUM( Pressure(1:n) * dBasisdx(1:n,i) )
2411              END DO
2412
2413!             - g^{jk} (\mu u^i_{,k})_{,j}):
2414!             ------------------------------
2415              DO j=1,DIM
2416                 DO k=1,DIM
2417                    Residual(i) = Residual(i) -   &
2418                         Metric(j,k) * Viscosity * &
2419                         SUM( Velocity(i,1:n) * ddBasisddx(1:n,j,k) )
2420
2421                    DO l=1,DIM
2422                       Residual(i) = Residual(i) +  &
2423                            Metric(j,k) * Viscosity * Symb(j,k,l) * &
2424                            SUM( Velocity(i,1:n) * dBasisdx(1:n,l) )
2425
2426                       Residual(i) = Residual(i) -  &
2427                            Metric(j,k) * Viscosity * Symb(l,j,i) * &
2428                            SUM( Velocity(l,1:n) * dBasisdx(1:n,k) )
2429
2430                       Residual(i) = Residual(i) -  &
2431                            Metric(j,k) * Viscosity * Symb(l,k,i) * &
2432                            SUM( Velocity(l,1:n) * dBasisdx(1:n,j) )
2433
2434                       Residual(i) = Residual(i) -  &
2435                            Metric(j,k) * Viscosity * dSymb(l,j,i,k) * &
2436                            SUM( Velocity(l,1:n) * Basis(1:n) )
2437
2438                       DO m=1,DIM
2439                          Residual(i) = Residual(i) - Metric(j,k) * Viscosity *&
2440                                  Symb(m,k,i) * Symb(l,j,m) * &
2441                                        SUM( Velocity(l,1:n) * Basis(1:n) )
2442
2443                          Residual(i) = Residual(i) + Metric(j,k) * Viscosity *&
2444                                  Symb(j,k,m) * Symb(l,m,i) * &
2445                                        SUM( Velocity(l,1:n) * Basis(1:n) )
2446                       END DO
2447                    END DO
2448                 END DO
2449              END DO
2450
2451!             - g^{ik} (\mu u^j_{,k})_{,j}):
2452!             ------------------------------
2453              DO j=1,DIM
2454                 DO k=1,DIM
2455                    Residual(i) = Residual(i) -   &
2456                         Metric(i,k) * Viscosity * &
2457                         SUM( Velocity(j,1:n) * ddBasisddx(1:n,j,k) )
2458
2459                    DO l=1,DIM
2460                       Residual(i) = Residual(i) +  &
2461                            Metric(i,k) * Viscosity * Symb(j,k,l) * &
2462                            SUM( Velocity(j,1:n) * dBasisdx(1:n,l) )
2463
2464                       Residual(i) = Residual(i) -  &
2465                            Metric(i,k) * Viscosity * Symb(l,j,j) * &
2466                            SUM( Velocity(l,1:n) * dBasisdx(1:n,k) )
2467
2468                       Residual(i) = Residual(i) -  &
2469                            Metric(i,k) * Viscosity * Symb(l,k,j) * &
2470                            SUM( Velocity(l,1:n) * dBasisdx(1:n,j) )
2471
2472                       Residual(i) = Residual(i) -  &
2473                            Metric(i,k) * Viscosity * dSymb(l,j,j,k) * &
2474                            SUM( Velocity(l,1:n) * Basis(1:n) )
2475
2476                       DO m=1,DIM
2477                          Residual(i) = Residual(i) - Metric(i,k) * Viscosity *&
2478                                  Symb(m,k,j) * Symb(l,j,m) * &
2479                                        SUM( Velocity(l,1:n) * Basis(1:n) )
2480
2481                          Residual(i) = Residual(i) + Metric(i,k) * Viscosity *&
2482                                  Symb(j,k,m) * Symb(l,m,j) * &
2483                                        SUM( Velocity(l,1:n) * Basis(1:n) )
2484                       END DO
2485                    END DO
2486                 END DO
2487              END DO
2488
2489              IF ( Convect ) THEN
2490!
2491!                + \rho * (@u/@t + u^j u^i_{,j}):
2492!                --------------------------------
2493                 Residual(i) = Residual(i) + Density *  &
2494                     SUM((Velocity(i,1:n)-PrevVelo(i,1:n))*Basis(1:n)) / dt
2495
2496                 DO j=1,DIM
2497                    Residual(i) = Residual(i) + &
2498                         Density * SUM( Velocity(j,1:n) * Basis(1:n) ) * &
2499                         SUM( Velocity(i,1:n) * dBasisdx(1:n,j) )
2500
2501                    DO k=1,DIM
2502                       Residual(i) = Residual(i) + &
2503                            Density * SUM( Velocity(j,1:n) * Basis(1:n) ) * &
2504                            Symb(j,k,i) * SUM( Velocity(k,1:n) * Basis(1:n) )
2505                    END DO
2506                 END DO
2507              END IF
2508           END IF
2509        END DO
2510
2511!
2512!       Continuity equation:
2513!       --------------------
2514        IF ( CurrentCoordinateSystem() == Cartesian ) THEN
2515!
2516!          + \rho * div(u):
2517!          ----------------
2518           DO j=1,DIM
2519              Residual(DIM+1) = Residual(DIM+1) + &
2520                   Density * SUM( Velocity(j,1:n) * dBasisdx(1:n,j) )
2521           END DO
2522
2523           IF ( Compressible ) THEN
2524!
2525!             + u.grad(\rho):
2526!             ----------------
2527              DO j=1,DIM
2528                 Residual(DIM+1) = Residual(DIM+1) + &
2529                      SUM( Velocity(j,1:n) * Basis(1:n) ) *  &
2530                           SUM( NodalDensity(1:n) * dBasisdx(1:n,j) )
2531              END DO
2532           END IF
2533        ELSE
2534!
2535!          + \rho * u^j_{,j}:
2536!          ------------------
2537           DO j=1,DIM
2538              Residual(DIM+1) = Residual(DIM+1) + &
2539                   Density * SUM( Velocity(j,1:n) * dBasisdx(1:n,j) )
2540
2541              DO k=1,DIM
2542                 Residual(DIM+1) = Residual(DIM+1) + Density * &
2543                      Symb(k,j,j) * SUM( Velocity(k,1:n) * Basis(1:n) )
2544              END DO
2545           END DO
2546
2547           IF ( Compressible ) THEN
2548!
2549!             + u^j \rho_{,j}:
2550!             ----------------
2551              DO j=1,DIM
2552                 Residual(DIM+1) = Residual(DIM+1) + &
2553                      SUM( Velocity(j,1:n) * Basis(1:n) ) *  &
2554                      SUM( NodalDensity(1:n) * dBasisdx(1:n,j) )
2555              END DO
2556           END IF
2557        END IF
2558
2559        DO i=1,DIM
2560           FNorm = FNorm + s * (Density * SUM(NodalForce(i,1:n)*Basis(1:n))**2)
2561        END DO
2562        Area = Area + s
2563
2564        IF ( CurrentCoordinateSystem() == Cartesian ) THEN
2565           ResidualNorm = ResidualNorm + &
2566             s * (Element % hK**2 * SUM(Residual(1:dim)**2) + Residual(dim+1)**2 )
2567        ELSE
2568           CALL InvertMatrix( Metric,3 )
2569           DO i=1,dim
2570              DO j=1,dim
2571                 ResidualNorm = ResidualNorm + &
2572                    s * Element % hK **2 * Metric(i,j) * Residual(i) * Residual(j)
2573              END DO
2574           END DO
2575           ResidualNorm = ResidualNorm + s * Residual(dim+1)**2
2576        END IF
2577     END DO
2578
2579!    FNorm = Area * FNorm
2580     Indicator = ResidualNorm
2581
2582     DEALLOCATE( NodalViscosity, NodalDensity, Basis, dBasisdx,           &
2583        ddBasisddx, Velocity, Pressure, PrevVelo, PrevPres, Temperature,   &
2584        NodalForce, HeatCapacity, ReferenceTemperature, HeatExpansionCoeff,&
2585        Nodes % x, Nodes % y, Nodes % z )
2586!------------------------------------------------------------------------------
2587  END FUNCTION FlowInsideResidual
2588!------------------------------------------------------------------------------
2589
2590!> \}
2591