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! *  Module for the solution of incompressible Navier-Stokes equation.
27! *  Utilizes multithreading and vectorization features initially introduced by Mikko Byckling.
28! *  Replaces partly the legacy solver FlowSolve which is not optimized.
29! *
30! *  Authors: Mika Malinen, Juhani Kataja, Juha Ruokolainen, Peter Råback
31! *  Email:   elmeradm@csc.fi
32! *  Web:     http://www.csc.fi/elmer
33! *  Address: CSC - IT Center for Science Ltd.
34! *           Keilaranta 14
35! *           02101 Espoo, Finland
36! *
37! *  Created: 28.01.2019
38! *
39!/*****************************************************************************/
40
41MODULE IncompressibleLocalForms
42
43  USE DefUtils
44
45  REAL(KIND=dp), ALLOCATABLE, SAVE :: bx(:), bxprev(:)
46
47CONTAINS
48
49!------------------------------------------------------------------------------
50! Assemble local finite element matrix for a single bulk element and glue
51! it to the global matrix. Routine is vectorized and multithreaded.
52!------------------------------------------------------------------------------
53  SUBROUTINE LocalBulkMatrix(Element, n, nd, ntot, dim, &
54       DivCurlForm, GradPVersion, SpecificLoad, StokesFlow, &
55       dt, LinearAssembly, nb, Newton, Transient, InitHandles, SchurSolver )
56!------------------------------------------------------------------------------
57    USE LinearForms
58    IMPLICIT NONE
59
60    TYPE(Element_t), POINTER, INTENT(IN) :: Element
61    INTEGER, INTENT(IN) :: n, nd, ntot, dim, nb
62    LOGICAL, INTENT(IN) :: DivCurlForm, GradPVersion, SpecificLoad, StokesFlow
63    REAL(KIND=dp), INTENT(IN) :: dt
64    LOGICAL, INTENT(IN) :: LinearAssembly, Newton, Transient, InitHandles
65    TYPE(Solver_t), POINTER :: SchurSolver
66!------------------------------------------------------------------------------
67    TYPE(GaussIntegrationPoints_t) :: IP
68    TYPE(Nodes_t) :: Nodes
69
70    LOGICAL :: Stat, Found
71
72    REAL(KIND=dp), TARGET :: MASS(ntot*(dim+1),ntot*(dim+1)), &
73        STIFF(ntot*(dim+1),ntot*(dim+1)), FORCE(ntot*(dim+1))
74    REAL(KIND=dp) :: NodalSol(dim+1,ntot)
75    REAL(KIND=dp) :: PrevNodalSol(dim+1,ntot)
76    REAL(KIND=dp) :: s, rho
77
78    REAL(KIND=dp), ALLOCATABLE, SAVE :: BasisVec(:,:), dBasisdxVec(:,:,:), DetJVec(:), &
79        rhoVec(:), VeloPresVec(:,:), loadAtIpVec(:,:), VelocityMass(:,:), &
80        PressureMass(:,:), ForcePart(:), &
81        weight_a(:), weight_b(:), weight_c(:), tauVec(:), PrevTempVec(:), PrevPressureVec(:), &
82        VeloVec(:,:), PresVec(:), GradVec(:,:,:)
83    REAL(KIND=dp), POINTER :: muVec(:), LoadVec(:)
84    REAL(KIND=dp), ALLOCATABLE :: muDerVec0(:),g(:,:,:),StrainRateVec(:,:,:)
85
86    REAL(kind=dp) :: stifford(ntot,ntot,dim+1,dim+1), muder, jacord(ntot,ntot,dim+1,dim+1), &
87                       JAC(ntot*(dim+1),ntot*(dim+1) ), jsum
88
89    INTEGER :: t, i, j, k, ii, p, q, ngp, allocstat
90    INTEGER, SAVE :: elemdim
91    INTEGER :: DOFs
92
93    TYPE(ValueHandle_t), SAVE :: Dens_h, Load_h(3)
94
95!DIR$ ATTRIBUTES ALIGN:64 :: BasisVec, dBasisdxVec, DetJVec, rhoVec, VeloPresVec, loadAtIpVec
96!DIR$ ATTRIBUTES ALIGN:64 :: MASS, STIFF, FORCE, weight_a, weight_b, weight_c
97!$OMP THREADPRIVATE(BasisVec, dBasisdxVec, DetJVec, rhoVec, VeloPresVec, loadAtIpVec, ElemDim )
98!$OMP THREADPRIVATE(VelocityMass, PressureMass, ForcePart, Weight_a, weight_b, weight_c)
99!$OMP THREADPRIVATE(tauVec, PrevTempVec, PrevPressureVec, VeloVec, PresVec, GradVec, Nodes)
100
101    SAVE Nodes
102!------------------------------------------------------------------------------
103
104    CALL GetElementNodesVec( Nodes )
105    STIFF = 0._dp
106    MASS  = 0._dp
107    FORCE = 0._dp
108    JAC   = 0._dp
109    JacOrd = 0._dp
110    stifford = 0._dp
111
112    DOFs = dim + 1
113
114    ! Numerical integration:
115    !-----------------------
116    IP = GaussPointsAdapt( Element, PReferenceElement = .TRUE. )
117    ngp = IP % n
118
119    ElemDim = Element % Type % Dimension
120
121    ! Storage size depending ngp
122    !-------------------------------------------------------------------------------
123
124    ! Deallocate storage if needed
125    IF (ALLOCATED(BasisVec)) THEN
126      IF (SIZE(BasisVec,1) < ngp .OR. SIZE(BasisVec,2) < ntot) &
127          DEALLOCATE(BasisVec,dBasisdxVec, DetJVec, rhoVec, VeloVec, PresVec, &
128          LoadAtIpVec, weight_a, weight_b, weight_c, tauVec, PrevTempVec, &
129          PrevPressureVec, VeloPresVec, GradVec)
130    END IF
131
132    ! Allocate storage if needed
133    IF (.NOT. ALLOCATED(BasisVec)) THEN
134      ALLOCATE(BasisVec(ngp,ntot), dBasisdxVec(ngp,ntot,3), DetJVec(ngp), &
135          rhoVec(ngp), VeloVec(ngp, dim), PresVec(ngp), velopresvec(ngp,dofs), LoadAtIpVec(ngp,dim+1), &
136          weight_a(ngp), weight_b(ngp), weight_c(ngp), tauVec(ngp), PrevTempVec(ngp), &
137          PrevPressureVec(ngp), GradVec(ngp,dim,dim), &
138          STAT=allocstat)
139      IF (allocstat /= 0) THEN
140        CALL Fatal('IncompressibleNSSolver::LocalBulkMatrix','Local storage allocation failed')
141      END IF
142    END IF
143
144
145    ! Deallocate storage (ntot) if needed
146    IF (ALLOCATED(VelocityMass)) THEN
147      IF(SIZE(VelocityMass,1) < ntot ) THEN
148        DEALLOCATE(VelocityMass, PressureMass, ForcePart)
149      END IF
150    END IF
151
152    ! Allocate storage (ntot) if needed
153    IF(.NOT. ALLOCATED(VelocityMass)) THEN
154      ALLOCATE(VelocityMass(ntot,ntot), PressureMass(ntot, ntot), &
155          ForcePart(ntot))
156    END IF
157
158    IF (Newton) THEN
159      ALLOCATE(muDerVec0(ngp), g(ngp,ntot,dim), StrainRateVec(ngp,dim,dim))
160      muDerVec0 = 0._dp
161    END IF
162
163    IF( InitHandles ) THEN
164      CALL ListInitElementKeyword( Dens_h,'Material','Density')
165      CALL ListInitElementKeyword( Load_h(1),'Body Force','Flow Bodyforce 1')
166      CALL ListInitElementKeyword( Load_h(2),'Body Force','Flow Bodyforce 2')
167      CALL ListInitElementKeyword( Load_h(3),'Body Force','Flow Bodyforce 3')
168    END IF
169
170    ! We assume constant density so far:
171    !-----------------------------------
172    rho = ListGetElementReal( Dens_h, Element = Element )
173
174    ! Get the previous elementwise velocity-pressure iterate:
175    !--------------------------------------------------------
176    IF ( LinearAssembly ) THEN
177      VeloPresVec = 0._dp
178    ELSE
179      CALL GetLocalSolution( NodalSol )
180      IF (nb > 0 .AND. Transient .AND. .NOT. StokesFlow) &
181         CALL GetLocalSolution(PrevNodalSol, tStep=-1)
182    END IF
183
184    VelocityMass = 0.0d0
185
186    ! Vectorized basis functions
187    stat = ElementInfoVec(Element, nodes, ngp, IP % U, IP % V, &
188        IP % W, detJvec, SIZE(basisVec, 2), BasisVec, dBasisdxVec)
189
190    ! Weights at integration points
191    DO t = 1, ngp
192      DetJVec(t) = DetJVec(t) * IP % s(t)
193    END DO
194
195    ! Velocity and pressure from previous iteration at integration points
196    IF(.NOT. LinearAssembly) THEN
197      CALL DGEMM('N', 'T', ngp, DOFs, n, 1._dp, BasisVec, ngp, NodalSol, &
198          dofs, 0._dp, VeloPresVec, ngp)
199    END IF
200
201    ! Return the effective viscosity. Currently only non-newtonian models supported.
202    muvec => EffectiveViscosityVec( ngp, BasisVec, dBasisdxVec, Element, NodalSol, &
203              muDerVec0, Newton,  InitHandles )
204
205    ! Rho
206    rhovec(1:ngp) = rho
207
208    ! Flow bodyforce if present
209    LoadAtIpVec = 0._dp
210    DO i=1,dim
211      LoadVec => ListGetElementRealVec( Load_h(i), ngp, BasisVec, Element, Found )
212      IF( Found ) THEN
213        IF (SpecificLoad) THEN
214          LoadAtIpVec(1:ngp,i) = LoadVec(1:ngp)
215        ELSE
216          LoadAtIpVec(1:ngp,i) = rho * LoadVec(1:ngp)
217        END IF
218      END IF
219    END DO
220
221    IF ( Newton ) THEN
222
223      DO i = 1,dim
224        DO j = 1,dim
225          GradVec(1:ngp, i, j) = MATMUL(dBasisdxVec(1:ngp,1:ntot,j),nodalsol(i,1:ntot))
226        END DO
227      END DO
228
229      IF( .NOT. StokesFlow ) THEN
230        DO i = 1,dim
231          LoadAtIpVec(1:ngp, i) = LoadAtIpVec(1:ngp, i) + rhovec(1:ngp) * &
232               SUM(gradvec(1:ngp,i,1:dim)*velopresvec(1:ngp,1:dim),2)
233        END DO
234      END IF
235
236      IF (ANY(muDerVec0(1:ngp)/=0)) THEN
237        DO i = 1,dim
238          DO j = 1,dim
239            StrainRateVec(1:ngp,i,j) = ( GradVec(1:ngp,i,j) + GradVec(1:ngp,j,i) ) / 2
240          END DO
241        END DO
242
243        muDerVec0(1:ngp) = muderVec0(1:ngp)*detJVec(1:ngp)*8
244        DO i=1,dim
245          DO q = 1,ntot
246            g(1:ngp,q,i) = SUM(StrainRateVec(1:ngp,i,:)*dBasisdxvec(1:ngp,q,1:dim),2)
247          END DO
248        END DO
249
250        DO i=1,dim
251          DO j=1,dim
252            CALL LinearForms_udotv(ngp,ntot,dim,g(:,:,j),g(:,:,i),mudervec0,jacord(:,:,j,i))
253          END DO
254        END DO
255      END IF
256    END IF
257
258    IF (DivCurlForm) THEN
259      weight_a(1:ngp) = muVec(1:ngp) * detJVec(1:ngp)
260      weight_b(1:ngp) = -muVec(1:ngp) * detJVec(1:ngp)
261
262      ! The following assumes that the bulk viscosity of the fluid vanishes:
263      weight_c(1:ngp) =  4.0_dp / 3.0_dp * muVec(1:ngp) * detJVec(1:ngp)
264
265      ! curl-curl part and div-div parts
266      SELECT CASE(dim)
267      CASE(3) ! {{{
268        i = 1; j = 1
269        ! curl-curl
270        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
271            dbasisdxvec(:, :, 2), dbasisdxvec(:,:,2), weight_a, stifford(:,:,i,j))
272        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
273            dbasisdxvec(:, :, 3), dbasisdxvec(:,:,3), weight_a, stifford(:,:,i,j))
274        ! div-div
275        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
276            dbasisdxvec(:, :, 1), dbasisdxvec(:,:, 1), weight_c, stifford(:,:,i,j))
277
278        i=1;j=2
279        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
280            dbasisdxvec(:, :, 2), dbasisdxvec(:,:,1), weight_b, stifford(:,:,i,j))
281        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
282            dbasisdxvec(:, :, 1), dbasisdxvec(:,:, 2), weight_c, stifford(:,:,i,j))
283
284        i=1;j=3
285        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
286            dbasisdxvec(:, :, 3), dbasisdxvec(:,:,1), weight_b, stifford(:,:,i,j))
287        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
288            dbasisdxvec(:, :, 1), dbasisdxvec(:,:, 3), weight_c, stifford(:,:,i,j))
289
290        i=2;j=1
291        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
292            dbasisdxvec(:, :, 1), dbasisdxvec(:,:,2), weight_b, stifford(:,:,i,j))
293        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
294            dbasisdxvec(:, :, 2), dbasisdxvec(:,:, 1), weight_c, stifford(:,:,i,j))
295
296        i=2;j=2
297        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
298            dbasisdxvec(:, :, 1), dbasisdxvec(:,:,1), weight_a, stifford(:,:,i,j))
299        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
300            dbasisdxvec(:, :, 3), dbasisdxvec(:,:,3), weight_a, stifford(:,:,i,j))
301        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
302            dbasisdxvec(:, :, 2), dbasisdxvec(:,:, 2), weight_c, stifford(:,:,i,j))
303
304        i=2;j=3
305        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
306            dbasisdxvec(:, :, 3), dbasisdxvec(:,:,2), weight_b, stifford(:,:,i,j))
307        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
308            dbasisdxvec(:, :, 2), dbasisdxvec(:,:, 3), weight_c, stifford(:,:,i,j))
309
310        i=3;j=1
311        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
312            dbasisdxvec(:, :, 1), dbasisdxvec(:,:,3), weight_b, stifford(:,:,i,j))
313        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
314            dbasisdxvec(:, :, 3), dbasisdxvec(:,:, 1), weight_c, stifford(:,:,i,j))
315
316        i=3;j=2
317        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
318            dbasisdxvec(:, :, 2), dbasisdxvec(:,:,3), weight_b, stifford(:,:,i,j))
319        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
320            dbasisdxvec(:, :, 3), dbasisdxvec(:,:, 2), weight_c, stifford(:,:,i,j))
321
322        i=3;j=3
323        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
324            dbasisdxvec(:, :, 1), dbasisdxvec(:,:,1), weight_a, stifford(:,:,i,j))
325        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
326            dbasisdxvec(:, :, 2), dbasisdxvec(:,:,2), weight_a, stifford(:,:,i,j))
327        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
328            dbasisdxvec(:, :, 3), dbasisdxvec(:,:, 3), weight_c, stifford(:,:,i,j))
329        ! }}}
330      CASE(2)  ! {{{
331        i = 1; j = 1
332        ! curl-curl
333        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
334            dbasisdxvec(:, :, 2), dbasisdxvec(:,:,2), weight_a, stifford(:,:,i,j))
335        ! div-div
336        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
337            dbasisdxvec(:, :, 1), dbasisdxvec(:,:, 1), weight_c, stifford(:,:,i,j))
338
339        i = 1; j = 2
340        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
341            dbasisdxvec(:, :, 2), dbasisdxvec(:,:,1), weight_b, stifford(:,:,i,j))
342        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
343            dbasisdxvec(:, :, 1), dbasisdxvec(:,:, 2), weight_c, stifford(:,:,i,j))
344
345        i = 2; j = 1
346        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
347            dbasisdxvec(:, :, 1), dbasisdxvec(:,:,2), weight_b, stifford(:,:,i,j))
348        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
349            dbasisdxvec(:, :, 2), dbasisdxvec(:,:, 1), weight_c, stifford(:,:,i,j))
350
351        i = 2; j = 2
352        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
353            dbasisdxvec(:, :, 1), dbasisdxvec(:,:,1), weight_a, stifford(:,:,i,j))
354        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
355            dbasisdxvec(:, :, 2), dbasisdxvec(:,:, 2), weight_c, stifford(:,:,i,j))
356        ! }}}
357      END SELECT
358
359    ELSE !DivCurlForm
360
361      weight_a(1:ngp) = muVec(1:ngp) * detJvec(1:ngp)
362
363      ! The following assumes that the bulk viscosity of the fluid vanishes:
364!     weight_c(1:ngp) = -2.0_dp / 3.0_dp * muVec(1:ngp) * detJVec(1:ngp)
365
366      DO i=1,dim
367        DO j=1,dim
368          CALL LinearForms_UdotV(ngp, ntot, elemdim, &
369              dBasisdxVec(1:ngp,1:ntot,j), dBasisdxVec(1:ngp,1:ntot,j), weight_a, stifford(1:ntot,1:ntot,i,i))
370
371          CALL LinearForms_UdotV(ngp, ntot, elemdim, &
372              dBasisdxVec(1:ngp,1:ntot,j), dBasisdxVec(1:ngp,1:ntot,i), weight_a, stifford(1:ntot,1:ntot,i,j))
373
374!         CALL LinearForms_UdotV(ngp, ntot, elemdim, &
375!             dBasisdxVec(1:ngp,1:ntot,i), dBasisdxVec(1:ngp,1:ntot,j), weight_c, stifford(1:ntot,1:ntot,i,j))
376        END DO
377      END DO
378    END IF
379
380    IF (GradPVersion) THEN
381       ! b(u,q) = (u, grad q) part
382      DO i = 1, dim
383        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
384            BasisVec, dbasisdxvec(:,:,i), detJVec, stifford(:,:,i,dofs))
385        StiffOrd(:,:,dofs,i) = transpose(stifford(:,:,i,dofs))
386      END DO
387    ELSE
388       DO i = 1, dim
389         CALL LinearForms_UdotV(ngp, ntot, elemdim, &
390            dBasisdxVec(:, :, i), BasisVec, -detJVec, StiffOrd(:,:,i,dofs))
391        StiffOrd(:,:,dofs,i) = transpose(stifford(:,:,i,dofs))
392      END DO
393    END IF
394
395    ! Masses (use symmetry)
396    ! Compute bilinear form G=G+(alpha u, u) = u .dot. (grad u)
397    IF ( .NOT. StokesFlow ) THEN
398      CALL LinearForms_UdotU(ngp, ntot, elemdim, BasisVec, DetJVec, VelocityMass, rhovec)
399
400      ! Scatter to the usual local mass matrix
401      DO i = 1, dim
402        mass(i::dofs, i::dofs) = mass(i::dofs, i::dofs) + VelocityMass(1:ntot, 1:ntot)
403      END DO
404
405      !mass(dofs::dofs, dofs::dofs) = mass(dofs::dofs, dofs::dofs) + PressureMass(1:ntot,1:ntot)
406
407      ! These loop unrolls look bad, maybe do nicer weight precomputation?
408      weight_a(1:ngp) = rhovec(1:ngp) * veloPresVec(1:ngp,1)
409      weight_b(1:ngp) = rhovec(1:ngp) * veloPresVec(1:ngp,2)
410      weight_c(1:ngp) = rhovec(1:ngp) * veloPresVec(1:ngp,3)
411      DO i = 1, dim
412        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
413            basisvec, dbasisdxvec(:,:,1), detJvec, stifford(:,:,i,i), weight_a)
414        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
415            basisvec, dbasisdxvec(:,:,2), detJvec, stifford(:,:,i,i), weight_b)
416        CALL LinearForms_UdotV(ngp, ntot, elemdim, &
417            basisvec, dbasisdxvec(:,:,3), detJvec, stifford(:,:,i,i), weight_c)
418      END DO
419
420      IF ( Newton ) THEN
421        DO i = 1, dim
422          DO j = 1, dim
423            CALL LinearForms_UdotV(ngp, ntot, elemdim, &
424                basisvec, basisvec, detJvec, stifford(:,:,i,j), rhovec*gradvec(:,i,j))
425          END DO
426        END DO
427      END IF
428    END IF
429
430    ! add loads
431    DO i = 1,dim+1
432      ForcePart = 0._dp
433      CALL LinearForms_UdotF(ngp, ntot, basisVec, detJVec, LoadAtIpVec(:,i), ForcePart)
434      FORCE(i::dofs) = ForcePart(1:ntot)
435    END DO
436
437    DO i = 1, DOFS
438      DO j = 1, DOFS
439        Stiff(i::DOFS, j::DOFS) = StiffOrd(1:ntot, 1:ntot, i,j)
440      END DO
441    END DO
442
443    IF ( Newton ) THEN
444BLOCK
445     REAL(KIND=dp) :: SOL(ntot*(dim+1))
446
447     SOL=0._dp
448     DO i = 1, DOFS
449       DO j = 1, DOFS
450         JAC(i::DOFS, j::DOFS) = JacOrd(1:ntot, 1:ntot, i,j)
451       END DO
452       SOL(i::DOFs) = NodalSol(i,1:ntot)
453     END DO
454
455     STIFF = STIFF + JAC
456     FORCE = FORCE + MATMUL(JAC,SOL)
457END BLOCK
458   END IF
459
460    IF(StokesFlow) THEN
461      IF ( nb>0 ) THEN
462        CALL LCondensate(nd, nb, dim, MASS, STIFF, FORCE)
463      ELSE
464        DO p = n+1,ntot
465          i = DOFs * p
466          FORCE(i)   = 0._dp
467          STIFF(i,:) = 0._dp
468          STIFF(:,i) = 0._dp
469          STIFF(i,i) = 1._dp
470        END DO
471      END IF
472
473    ELSE IF (nb > 0 .AND. nd==n .AND. Transient) THEN
474      !-------------------------------------------------------------------------
475      ! This branch is primarily intended to handle the (enhanced) MINI element
476      ! approximation together with the static condensation for the velocity
477      ! bubbles. The subroutine LCondensate constructs the time derivative of
478      ! the bubble-augmented part and performs the static condensation.
479      !-------------------------------------------------------------------------
480      CALL LCondensate(nd, nb, dim, MASS, STIFF, FORCE, PrevNodalSol, &
481                   NodalSol, Element % ElementIndex)
482    ELSE
483      !-------------------------------------------------------------------------
484      ! The cases handled here include the MINI element approximation with the
485      ! velocity bubbles left in the global system and P2/Q2-P1/Q1 approximation.
486      ! First, enforce P1/Q1 pressure approximation by setting Dirichlet
487      ! constraints for unused dofs:
488      !-------------------------------------------------------------------------
489      DO p = n+1,ntot
490        i = DOFs * p
491        FORCE(i)   = 0._dp
492        MASS(:,i)  = 0._dp
493        MASS(i,:)  = 0._dp
494        STIFF(i,:) = 0._dp
495        STIFF(:,i) = 0._dp
496        STIFF(i,i) = 1._dp
497      END DO
498
499      IF ( Transient ) THEN
500        CALL Default1stOrderTime( MASS, STIFF, FORCE )
501      END IF
502      IF (nb > 0) THEN
503        IF (Transient) THEN
504          CALL LCondensate(nd, nb, dim, MASS, STIFF, FORCE, &
505              PrevNodalSol, NodalSol, Element % ElementIndex)
506        ELSE
507          CALL LCondensate(nd, nb, dim, MASS, STIFF, FORCE)
508        END IF
509      END IF
510    END IF
511
512    CALL DefaultUpdateEquations( STIFF, FORCE, UElement=Element, VecAssembly=.TRUE.)
513
514    IF( ASSOCIATED( SchurSolver ) ) THEN
515      ! Preconditioner for pressure block when using block preconditioning
516      weight_a(1:ngp) = -1.0_dp / muvec(1:ngp) * detJVec(1:ngp)
517      PressureMass = 0.0_dp
518      FORCE = 0.0_dp
519      CALL LinearForms_UdotU(ngp, nd, elemdim, BasisVec, weight_a, PressureMass)
520      CALL DefaultUpdateEquations( PressureMass, FORCE, UElement=Element, &
521                 Usolver = SchurSolver, VecAssembly = .TRUE.)
522    END IF
523
524
525
526!------------------------------------------------------------------------------
527
528  CONTAINS
529
530
531    FUNCTION EffectiveViscosityVec( ngp, BasisVec, dBasisdxVec, Element, NodalSol, &
532        ViscDerVec, ViscNewton, InitHandles ) RESULT ( EffViscVec )
533
534      INTEGER :: ngp
535      REAL(KIND=dp) :: BasisVec(:,:), dBasisdxVec(:,:,:)
536      TYPE(Element_t), POINTER :: Element
537      REAL(KIND=dp) :: NodalSol(:,:), ViscDerVec(:)
538      LOGICAL :: InitHandles , ViscNewton
539      REAL(KIND=dp), POINTER  :: EffViscVec(:)
540
541      LOGICAL :: Found
542      CHARACTER(LEN=MAX_NAME_LEN) :: ViscModel
543      TYPE(ValueHandle_t), SAVE :: Visc_h, ViscModel_h, ViscExp_h, ViscCritical_h, &
544          ViscNominal_h, ViscDiff_h, ViscTrans_h, ViscYasuda_h, ViscGlenExp_h, ViscGlenFactor_h, &
545          ViscArrSet_h, ViscArr_h, ViscTLimit_h, ViscRate1_h, ViscRate2_h, ViscEne1_h, ViscEne2_h, &
546          ViscTemp_h
547      REAL(KIND=dp), SAVE :: R
548      REAL(KIND=dp) :: c1, c2, c3, c4, Ehf, Temp, Tlimit, ArrheniusFactor, A1, A2, Q1, Q2, ViscCond
549      LOGICAL, SAVE :: ConstantVisc = .FALSE., Visited = .FALSE.
550      REAL(KIND=dp), ALLOCATABLE, SAVE :: ss(:), s(:), ArrheniusFactorVec(:)
551      REAL(KIND=dp), POINTER, SAVE :: ViscVec0(:), ViscVec(:), TempVec(:), EhfVec(:)
552
553!$OMP THREADPRIVATE(ss,s,ViscVec0,ViscVec,ArrheniusFactorVec)
554
555      IF(InitHandles ) THEN
556        CALL Info('EffectiveViscosityVec','Initializing handles for viscosity models',Level=8)
557
558        CALL ListInitElementKeyword( Visc_h,'Material','Viscosity')
559        CALL ListInitElementKeyword( ViscModel_h,'Material','Viscosity Model')
560
561        IF( ListGetElementSomewhere( ViscModel_h) ) THEN
562          ViscCond = ListGetCReal( CurrentModel % Solver % Values,&
563              'Newtonian Viscosity Condition',Found )
564          ConstantVisc = ( Found .AND. ViscCond > 0.0_dp )
565
566          IF( ListGetLogical( CurrentModel % Solver % Values,&
567              'Constant-Viscosity Start', Found) ) ConstantVisc = (.NOT. Visited )
568
569          CALL ListInitElementKeyword( ViscExp_h,'Material','Viscosity Exponent')
570          CALL ListInitElementKeyword( ViscCritical_h,'Material','Critical Shear Rate')
571          CALL ListInitElementKeyword( ViscNominal_h,'Material','Nominal Shear Rate')
572          CALL ListInitElementKeyword( ViscDiff_h,'Material','Viscosity Difference')
573          CALL ListInitElementKeyword( ViscTrans_h,'Material','Viscosity Transition')
574          CALL ListInitElementKeyword( ViscYasuda_h,'Material','Yasuda Exponent')
575
576          ! Do these initializations for glen's model only
577          IF ( ListCompareElementAnyString( ViscModel_h,'glen') ) THEN
578            CALL ListInitElementKeyword( ViscGlenExp_h,'Material','Glen Exponent',DefRValue=3.0_dp)
579            CALL ListInitElementKeyword( ViscGlenFactor_h,'Material','Glen Enhancement Factor',DefRValue=1.0_dp)
580            CALL ListInitElementKeyword( ViscArrSet_h,'Material','Set Arrhenius Factor',DefLValue=.FALSE.)
581            CALL ListInitElementKeyword( ViscArr_h,'Material','Arrhenius Factor')
582            CALL ListInitElementKeyword( ViscTLimit_h,'Material','Limit Temperature',DefRValue=-10.0_dp)
583            CALL ListInitElementKeyword( ViscRate1_h,'Material','Rate Factor 1',DefRValue=3.985d-13)
584            CALL ListInitElementKeyword( ViscRate2_h,'Material','Rate Factor 2',DefRValue=1.916d3)
585            CALL ListInitElementKeyword( ViscEne1_h,'Material','Activation Energy 1',DefRValue=60.0d03)
586            CALL ListInitElementKeyword( ViscEne2_h,'Material','Activation Energy 2',DefRValue=139.0d03)
587            CALL ListInitElementKeyword( ViscTemp_h,'Material','Relative Temperature')
588
589            IF (.NOT.ListCheckPresentAnyMaterial( CurrentModel,'Glen Allow Old Keywords')) THEN
590              IF( ListCheckPresentAnyMaterial( CurrentModel,'Constant Temperature') ) THEN
591                CALL Fatal('EffectiveViscosityVec','Replace >Constant Temperature< with >Relative Temperature<')
592              END IF
593              IF( ListCheckPresentAnyMaterial( CurrentModel,'Temperature Field Variable') ) THEN
594                CALL Fatal('EffectiveViscosityVec','Replace >Temperature Field Variable< with >Relative Temperature<')
595              END IF
596            END IF
597            IF( ListCheckPresentAnyMaterial( CurrentModel,'Glen Enhancement Factor Function')  ) THEN
598              CALL Fatal('EffectiveViscosityVec','No Glen function API yet!')
599            END IF
600            R = GetConstReal( CurrentModel % Constants,'Gas Constant',Found)
601            IF (.NOT.Found) R = 8.314_dp
602          END IF
603        END IF
604
605        Visited = .TRUE.
606      END IF
607
608      ViscVec0 => ListGetElementRealVec( Visc_h, ngp, BasisVec, Element )
609
610      ViscModel = ListGetElementString( ViscModel_h, Element, Found )
611      IF( .NOT. Found ) THEN
612        ! Return the plain viscosity
613        EffViscVec => ViscVec0
614        RETURN
615      END IF
616
617      ! Initialize derivative of viscosity for when newtonian linearization is used
618      IF( ViscNewton ) THEN
619        ViscDerVec(1:ngp) = 0.0_dp
620      END IF
621
622      ! This reverts the viscosity model to linear
623      IF( ConstantVisc ) THEN
624        EffViscVec => ViscVec0
625        RETURN
626      END IF
627
628      ! Deallocate too small storage if needed
629      IF (ALLOCATED(ss)) THEN
630        IF (SIZE(ss) < ngp ) DEALLOCATE(ss, s, ViscVec, ArrheniusFactorVec )
631      END IF
632
633      ! Allocate storage if needed
634      IF (.NOT. ALLOCATED(ss)) THEN
635        ALLOCATE(ss(ngp),s(ngp),ViscVec(ngp),ArrheniusFactorVec(ngp),STAT=allocstat)
636        IF (allocstat /= 0) THEN
637          CALL Fatal('IncompressibleNSSolver::LocalBulkMatrix','Local storage allocation failed')
638        END IF
639      END IF
640
641      ! For non-newtonian models compute the viscosity here
642      EffViscVec => ViscVec
643
644      ! Calculate the strain rate velocity at all integration points
645      ss(1:ngp) = 0.0_dp
646      DO i=1,dim
647        DO j=1,dim
648          s(1:ngp) = MATMUL( dBasisdxVec(1:ngp,1:ntot,i), nodalsol(j,1:ntot) ) + &
649              MATMUL( dBasisdxVec(1:ngp,1:ntot,j), nodalsol(i,1:ntot) )
650          ss(1:ngp) = ss(1:ngp) + s(1:ngp)**2
651        END DO
652      END DO
653      ss(1:ngp) = 0.5_dp * ss(1:ngp)
654
655
656
657
658      SELECT CASE( ViscModel )
659
660      CASE('glen')
661        c2 = ListGetElementReal( ViscGlenExp_h,Element=Element,Found=Found)
662
663        ! the second invariant is not taken from the strain rate tensor,
664        ! but rather 2*strain rate tensor (that's why we divide by 4 = 2**2)
665        s(1:ngp) = ss(1:ngp)/4.0_dp
666
667        c3 = ListGetElementReal( ViscCritical_h,Element=Element,Found=Found)
668        IF( Found ) THEN
669          c3 = c3**2
670          WHERE( s(1:ngp) < c3 ) s(1:ngp) = c3
671        END IF
672
673        IF( ListGetElementLogical( ViscArrSet_h,Element,Found=Found) ) THEN
674          ArrheniusFactor = ListGetElementReal( ViscArr_h,Element=Element)
675          ViscVec(1:ngp) = 0.5_dp * (ArrheniusFactor)**(-1.0_dp/c2) * s(1:ngp)**(((1.0_dp/c2)-1.0_dp)/2.0_dp);
676
677          IF( ViscNewton ) THEN
678            WHERE( s(1:ngp) > c3 ) ViscDerVec(1:ngp) = 0.5_dp * ArrheniusFactor**(-1.0_dp/c2) &
679                * ((1.0_dp/c2)-1.0_dp)/2.0_dp * s(1:ngp)**(((1.0_dp/c2)-1.0_dp)/2.0_dp - 1.0_dp)/4.0_dp
680          END IF
681        ELSE
682          ! lets for the time being have this hardcoded
683          Tlimit = ListGetElementReal( ViscTlimit_h,Element=Element)
684          A1 = ListGetElementReal( ViscRate1_h,Element=Element)
685          A2 = ListGetElementReal( ViscRate2_h,Element=Element)
686          Q1 = ListGetElementReal( ViscEne1_h,Element=Element)
687          Q2 = ListGetElementReal( ViscEne2_h,Element=Element)
688#if 1
689          ! WHERE is faster than DO + IF
690          TempVec => ListGetElementRealVec( ViscTemp_h, ngp, BasisVec, Element )
691
692          WHERE( TempVec(1:ngp ) < Tlimit )
693            ArrheniusFactorVec(1:ngp) = A1 * EXP( -Q1/(R * (273.15_dp + TempVec(1:ngp))))
694          ELSE WHERE( TempVec(1:ngp) > 0.0_dp )
695            ArrheniusFactorVec(1:ngp) = A2 * EXP( -Q2/(R * (273.15_dp)))
696          ELSE WHERE
697            ArrheniusFactorVec(1:ngp) = A2 * EXP( -Q2/(R * (273.15_dp + TempVec(1:ngp))))
698          END WHERE
699
700          EhfVec => ListGetElementRealVec( ViscGlenFactor_h, ngp, BasisVec,Element=Element )
701          ViscVec(1:ngp) = 0.5_dp * (EhFVec(1:ngp) * ArrheniusFactorVec(1:ngp))**(-1.0_dp/c2) * &
702              s(1:ngp)**(((1.0_dp/c2)-1.0_dp)/2.0_dp);
703
704          IF( ViscNewton ) THEN
705            WHERE( s(1:ngp) > c3 )
706              ViscDerVec(1:ngp) = 0.5_dp * (  EhFVec(1:ngp) * ArrheniusFactorVec(1:ngp))**(-1.0_dp/c2) &
707                    * ((1.0_dp/c2)-1.0_dp)/2.0_dp * s(1:ngp)**(((1.0_dp/c2)-1.0_dp)/2.0_dp - 1.0_dp)/4.0_dp
708            END WHERE
709          END IF
710#else
711          DO i=1,ngp
712            Temp = ListGetElementReal( ViscTemp_h,BasisVec(i,:),Element=Element,&
713                Found=Found,GaussPoint=i)
714
715            IF( Temp <= Tlimit ) THEN
716              ArrheniusFactor = A1 * EXP( -Q1/(R * (273.15_dp + Temp)))
717            ELSE IF((Tlimit < Temp ) .AND. (Temp <= 0.0_dp)) THEN
718              ArrheniusFactor = A2 * EXP( -Q2/(R * (273.15_dp + Temp)))
719            ELSE
720              ArrheniusFactor = A2 * EXP( -Q2/(R * (273.15_dp)))
721              CALL Info('EffectiveViscosityVec',&
722                  'Positive Temperature detected in Glen - limiting to zero!', Level = 12)
723            END IF
724
725            Ehf = ListGetElementReal( ViscGlenFactor_h,BasisVec(i,:),Element=Element,&
726                Found=Found,GaussPoint=i)
727
728            ! compute the effective viscosity
729            ViscVec(i) = 0.5_dp * (EhF * ArrheniusFactor)**(-1.0_dp/c2) * s(i)**(((1.0_dp/c2)-1.0_dp)/2.0_dp);
730
731            IF( ViscNewton ) THEN
732              IF( s(i) > c3 ) THEN
733                ViscDerVec(i) = 0.5_dp * (  EhF * ArrheniusFactor)**(-1.0_dp/c2) &
734                    * ((1.0_dp/c2)-1.0_dp)/2.0_dp * s(i)**(((1.0_dp/c2)-1.0_dp)/2.0_dp - 1.0_dp)/4.0_dp
735              END IF
736            END IF
737
738          END DO
739#endif
740
741        END IF
742
743
744      CASE('power law')
745        c2 = ListGetElementReal( ViscExp_h,Element=Element)
746
747        c3 = ListGetElementReal( ViscCritical_h,Element=Element,Found=Found)
748        IF( Found ) THEN
749          c3 = c3**2
750          WHERE( ss(1:ngp) < c3 ) ss(1:ngp) = c3
751        END IF
752
753        ViscVec(1:ngp) = ViscVec0(1:ngp) * ss(1:ngp)**((c2-1)/2)
754
755        IF (ViscNewton ) THEN
756          WHERE(ss(1:ngp) /= 0) ViscDerVec(1:ngp) = &
757              ViscVec0(1:ngp) * (c2-1)/2 * ss(1:ngp)**((c2-1)/2-1)
758        END IF
759
760        c4 = ListGetElementReal( ViscNominal_h,Element=Element,Found=Found)
761        IF( Found ) THEN
762          ViscVec(1:ngp) = ViscVec(1:ngp) / c4**(c2-1)
763          IF (ViscNewton ) THEN
764            ViscDerVec(1:ngp) = ViscDerVec(1:ngp) / c4**(c2-1)
765          END IF
766        END IF
767
768      CASE('power law too')
769        c2 = ListGetElementReal( ViscExp_h,Element=Element)
770        ViscVec(1:ngp) = ViscVec0(1:ngp)**(-1/c2)* ss(1:ngp)**(-(c2-1)/(2*c2)) / 2
771
772        IF (ViscNewton ) THEN
773          ViscDerVec(1:ngp) = ViscVec0(1:ngp)**(-1/c2)*(-(c2-1)/(2*c2))*ss(1:ngp)*(-(c2-1)/(2*c2)-1) / 2
774        END IF
775
776      CASE ('carreau')
777        c1 = ListGetElementReal( ViscDiff_h,Element=Element)
778        c2 = ListGetElementReal( ViscExp_h,Element=Element)
779        c3 = ListGetElementReal( ViscTrans_h,Element=Element)
780        c4 = ListGetElementReal( ViscYasuda_h,Element=Element,Found=Found)
781        IF( Found ) THEN
782          ViscVec(1:ngp) = ViscVec0(1:ngp) + c1 * (1 + c3**c4*ss(1:ngp)**(c4/2))**((c2-1)/c4)
783
784          IF( ViscNewton ) THEN
785            ViscDerVec(1:ngp) = c1*(1+c3**c4*ss(1:ngp)**(c4/2))**((c2-1)/c4-1)*(c2-1)/2*c3**c4*&
786                ss(1:ngp)**(c4/2-1)
787          END IF
788        ELSE
789          ViscVec(1:ngp) = ViscVec0(1:ngp) + c1 * (1 + c3*c3*ss(1:ngp))**((c2-1)/2)
790
791          IF( ViscNewton ) THEN
792            ViscDerVec(1:ngp) = c1*(c2-1)/2*c3**2*(1+c3**2*ss(1:ngp))**((c2-1)/2-1)
793          END IF
794        END IF
795
796      CASE ('cross')
797        c1 = ListGetElementReal( ViscDiff_h,Element=Element)
798        c2 = ListGetElementReal( ViscExp_h,Element=Element)
799        c3 = ListGetElementReal( ViscTrans_h,Element=Element)
800
801        ViscVec(1:ngp) = ViscVec0(1:ngp) + c1 / (1 + c3*ss(1:ngp)**(c2/2))
802
803        IF( ViscNewton ) THEN
804          ViscDerVec(1:ngp) = -c1*c3*ss(1:ngp)**(c2/2)*c2 / (2*(1+c3*ss(1:ngp)**(c2/2))**2*ss(1:ngp))
805        END IF
806
807      CASE ('powell eyring')
808        c1 = ListGetElementReal( ViscDiff_h,Element=Element)
809        c2 = ListGetElementReal( ViscTrans_h,Element=Element)
810
811        s(1:ngp) = SQRT(ss(1:ngp))
812
813        IF( ViscNewton ) THEN
814          WHERE( c2*s(1:ngp) < 1.0d-5 )
815            ViscVec(1:ngp) = ViscVec0(1:ngp) + c1
816            ViscDerVec(1:ngp) = 0.0_dp
817          ELSE WHERE
818            ViscVec(1:ngp) = ViscVec0(1:ngp) + c1 * LOG(c2*s(1:ngp)+SQRT(c2*c2*ss(1:ngp)+1))/(c2*ss(1:ngp))
819            ViscDerVec(1:ngp) = c1*(c2/(2*s(1:ngp))+c2**2/(2*SQRT(c2**2*ss(1:ngp)+1)))/ &
820                ((c2*s(1:ngp)+SQRT(c2*ss(1:ngp)+1))*c2*s(1:ngp)) - &
821                c1*LOG(c2*s(1:ngp)+SQRT(c2**2*ss(1:ngp)+1))/(c2*s(1:ngp)**3)/2
822          END WHERE
823        ELSE
824          WHERE( c2*s(1:ngp) < 1.0d-5 )
825            ViscVec(1:ngp) = ViscVec0(1:ngp) + c1
826          ELSE WHERE
827            ViscVec(1:ngp) = ViscVec0(1:ngp) + c1 * LOG(c2*s(1:ngp)+SQRT(c2*c2*ss(1:ngp)+1))/(c2*ss(1:ngp))
828          END WHERE
829        END IF
830
831      CASE DEFAULT
832        CALL Fatal('EffectiveViscosityVec','Unknown material model')
833
834      END SELECT
835
836
837    END FUNCTION EffectiveViscosityVec
838
839
840
841    !------------------------------------------------------------------------------
842    ! A special subroutine for performing static condensation when the velocity
843    ! (the first dim components of the solution) is augmented by a bubble part.
844    ! The speciality is that the bubble part at the previous time level
845    ! is retrieved to approximate the first time derivative in a consistent manner.
846    ! This version works only with the BDF(1) method, as higher-order versions have
847    ! not yet been implemented.
848    !------------------------------------------------------------------------------
849    SUBROUTINE LCondensate( N, nb, dim, M, K, F, xprev, x, Element_id )
850      !------------------------------------------------------------------------------
851      USE LinearAlgebra
852      IMPLICIT NONE
853
854      INTEGER, INTENT(IN) :: N   ! The number of retained DOFs per scalar field
855      INTEGER, INTENT(IN) :: nb  ! The number of eliminated DOFs per scalar field
856      INTEGER, INTENT(IN) :: dim ! The number of bubble-augmented fields
857      REAL(KIND=dp), INTENT(IN) :: M(:,:)
858      REAL(KIND=dp), INTENT(INOUT) :: K(:,:), F(:)
859      REAL(KIND=dp), INTENT(IN), OPTIONAL :: x(:,:)     ! The solution without
860      ! bubbles
861      REAL(KIND=dp), INTENT(IN), OPTIONAL :: xprev(:,:) ! The previous solution
862      ! without bubbles
863      INTEGER, INTENT(IN), OPTIONAL :: Element_id       ! The element identifier
864
865      ! This subroutine accesses also the real-valued arrays bx(:) and bxprev(:)
866      ! that contain coefficients of the bubble basis functions
867
868  !------------------------------------------------------------------------------
869      LOGICAL :: ComputeBubblePart
870      REAL(KIND=dp) :: Kbb(nb*dim,nb*dim), Fb(nb*dim)
871      REAL(KIND=dp) :: Kbl(nb*dim,n*(dim+1)), Klb(n*(dim+1), nb*dim)
872
873      REAL(KIND=dp) :: xl(n*(dim+1)), xlprev((n+nb)*(dim+1))
874
875      INTEGER :: DOFs
876      INTEGER :: p, q, Cdofs((dim+1)*n), Bdofs(dim*nb)
877  !------------------------------------------------------------------------------
878      ComputeBubblePart = PRESENT(x) .AND. PRESENT(xprev) .AND. PRESENT(Element_id)
879
880      DOFs = dim + 1
881
882      ! Vectorize the input array x and
883      ! create xlprev that contains the full previous solution including
884      ! the bubble DOFs. First insert the DOFs that are retained:
885      q = 0
886      DO p = 1,n
887        DO i = 1,DOFs
888          q = q + 1
889          Cdofs(q) = DOFs*(p-1) + i
890        END DO
891      END DO
892
893      ! Then the DOFs of the bubble part:
894      q = 0
895      DO p = 1,nb
896        DO i = 1,dim
897          q = q + 1
898          Bdofs(q) = DOFs*(p-1) + i + n*DOFs
899        END DO
900      END DO
901
902      ! The following only works for the BDF(1) method:
903      IF (ComputeBubblePart) THEN
904        xlprev = 0
905        q = 0
906        DO p = 1,n
907          DO i = 1,DOFs
908            q = q + 1
909            xl(q) = x(i,p)
910            xlprev(cdofs(q)) = xprev(i,p) ! cdofs identity mapping?
911          END DO
912        END DO
913
914        q = 0
915        DO p = 1,nb
916          DO i = 1,dim
917            q = q + 1
918            xlprev(bdofs(q)) = bxprev((Element_id-1)*dim*nb+q)
919          END DO
920        END DO
921
922        K = K + M / dt
923        F = F + MATMUL(M,xlprev) / dt
924      END IF
925
926      Kbb = K(Bdofs,Bdofs)
927      Kbl = K(Bdofs,Cdofs)
928      Klb = K(Cdofs,Bdofs)
929      Fb  = F(Bdofs)
930
931      CALL InvertMatrix( Kbb,Nb*dim )
932
933      F(cdofs) = F(cdofs) - MATMUL( Klb, MATMUL( Kbb, Fb ) )
934      K(cdofs,cdofs) = &
935          K(cdofs,cdofs) - MATMUL( Klb, MATMUL( Kbb,Kbl ) )
936
937      ! The bubble part evaluated for the current solution candidate:
938      IF (ComputeBubblePart) bx((Element_id-1)*dim*nb+1:Element_id*dim*nb) = &
939          MATMUL(Kbb,Fb-MATMUL(Kbl,xl))
940      !------------------------------------------------------------------------------
941    END SUBROUTINE LCondensate
942    !------------------------------------------------------------------------------
943
944  END SUBROUTINE LocalBulkMatrix
945
946
947!------------------------------------------------------------------------------
948! Assemble local finite element matrix for a single boundary element and glue
949! it to the global matrix.
950!------------------------------------------------------------------------------
951  SUBROUTINE LocalBoundaryMatrix( Element, n, nd, dim, InitHandles)
952!------------------------------------------------------------------------------
953    IMPLICIT NONE
954
955    TYPE(Element_t), POINTER, INTENT(IN) :: Element
956    INTEGER, INTENT(IN) :: n, nd, dim
957    LOGICAL, INTENT(INOUT) :: InitHandles
958!------------------------------------------------------------------------------
959    TYPE(GaussIntegrationPoints_t) :: IP
960    REAL(KIND=dp), TARGET :: STIFF(nd*(dim+1),nd*(dim+1)), FORCE(nd*(dim+1))
961    REAL(KIND=dp), ALLOCATABLE :: Basis(:)
962    INTEGER :: c,i,j,k,l,p,q,t,ngp
963    LOGICAL :: NormalTangential, HaveSlip, HaveForce, HavePres, Found, Stat
964    REAL(KIND=dp) :: ExtPressure, s, detJ
965    REAL(KIND=dp) :: SlipCoeff(3), SurfaceTraction(3), Normal(3), Tangent(3), Tangent2(3), Vect(3)
966    TYPE(Nodes_t), SAVE :: Nodes
967    TYPE(ValueHandle_t), SAVE :: ExtPressure_h, SurfaceTraction_h, SlipCoeff_h, &
968                NormalTangential_h, NormalTangentialVelo_h
969
970    SAVE Basis
971
972!------------------------------------------------------------------------------
973
974    IF( InitHandles ) THEN
975      CALL ListInitElementKeyword( ExtPressure_h,'Boundary Condition','Normal Surface Traction')
976      IF( .NOT. ListGetElementSomewhere( ExtPressure_h) ) THEN
977        CALL ListInitElementKeyword( ExtPressure_h,'Boundary Condition','External Pressure')
978      END IF
979      CALL ListInitElementKeyword( SurfaceTraction_h,'Boundary Condition','Surface Traction',InitVec3D=.TRUE.)
980      CALL ListInitElementKeyword( SlipCoeff_h,'Boundary Condition','Slip Coefficient',InitVec3D=.TRUE.)
981
982      CALL ListInitElementKeyword( NormalTangentialVelo_h,'Boundary Condition',&
983             'Normal-Tangential Velocity' )
984
985      CALL ListInitElementKeyword( NormalTangential_h,'Boundary Condition',&
986          'Normal-Tangential '//GetVarName(CurrentModel % Solver % Variable) )
987
988      InitHandles = .FALSE.
989    END IF
990
991    IF( ALLOCATED( Basis ) ) THEN
992      IF( SIZE( Basis ) < nd ) THEN
993        DEALLOCATE( Basis )
994      END IF
995    END IF
996
997    IF( .NOT. ALLOCATED( Basis ) ) THEN
998      ALLOCATE( Basis(nd) )
999    END IF
1000
1001
1002    CALL GetElementNodes( Nodes )
1003    STIFF = 0.0d0
1004    FORCE = 0.0d0
1005    c = dim + 1
1006
1007    ! Numerical integration:
1008    !-----------------------
1009    IP = GaussPoints( Element )
1010    ngp = IP % n
1011
1012    NormalTangential = ListGetElementLogical( NormalTangentialVelo_h, Element, Found )
1013    IF (.NOT.Found) THEN
1014        NormalTangential = ListGetElementLogical( NormalTangential_h, Element, Found )
1015    END IF
1016
1017    DO t=1,ngp
1018!------------------------------------------------------------------------------
1019!    Basis function values & derivatives at the integration point
1020!------------------------------------------------------------------------------
1021      stat = ElementInfo( Element,Nodes,IP % u(t),IP % v(t),IP % w(t), detJ, Basis )
1022
1023      s = detJ * IP % s(t)
1024
1025      ! Slip coefficient
1026      !----------------------------------
1027      SlipCoeff = ListGetElementReal3D( SlipCoeff_h, Basis, Element, HaveSlip, GaussPoint = t )
1028
1029      ! Given force on a boundary componentwise
1030      !----------------------------------------
1031      SurfaceTraction = ListGetElementReal3D( SurfaceTraction_h, Basis, Element, HaveForce, GaussPoint = t )
1032
1033      ! Given force to the normal direction
1034      !------------------------------------
1035      ExtPressure = ListGetElementReal( ExtPressure_h, Basis, Element, HavePres, GaussPoint = t )
1036
1037      ! Nothing to do, exit the routine
1038      IF(.NOT. (HaveSlip .OR. HaveForce .OR. HavePres)) RETURN
1039
1040      ! Calculate normal vector only if needed
1041      IF( HavePres .OR. NormalTangential ) THEN
1042        Normal = NormalVector( Element, Nodes, IP % u(t),IP %v(t),.TRUE. )
1043      END IF
1044
1045      ! Project external pressure to the normal direction
1046      IF( HavePres ) THEN
1047        IF( NormalTangential ) THEN
1048          SurfaceTraction(1) = SurfaceTraction(1) + ExtPressure
1049        ELSE
1050          SurfaceTraction = SurfaceTraction + ExtPressure * Normal
1051        END IF
1052        HaveForce = .TRUE.
1053      END IF
1054
1055      ! Calculate directions for N-T system
1056      IF( NormalTangential ) THEN
1057        SELECT CASE( dim )
1058        CASE(2)
1059          Tangent(1) =  Normal(2)
1060          Tangent(2) = -Normal(1)
1061          Tangent(3) =  0.0_dp
1062          Tangent2   =  0.0_dp
1063        CASE(3)
1064          CALL TangentDirections( Normal, Tangent, Tangent2 )
1065        END SELECT
1066      END IF
1067
1068      ! Assemble the slip coefficients to the stiffness matrix
1069      IF( HaveSlip ) THEN
1070        IF ( NormalTangential ) THEN
1071          DO i=1,dim
1072            SELECT CASE(i)
1073            CASE(1)
1074              Vect = Normal
1075            CASE(2)
1076              Vect = Tangent
1077            CASE(3)
1078              Vect = Tangent2
1079            END SELECT
1080
1081            DO p=1,nd
1082              DO q=1,nd
1083                DO j=1,dim
1084                  DO k=1,dim
1085                    STIFF( (p-1)*c+j,(q-1)*c+k ) = &
1086                        STIFF( (p-1)*c+j,(q-1)*c+k ) + &
1087                        s * SlipCoeff(i) * Basis(q) * Basis(p) * Vect(j) * Vect(k)
1088                  END DO
1089                END DO
1090              END DO
1091            END DO
1092          END DO
1093        ELSE
1094          DO p=1,nd
1095            DO q=1,nd
1096              DO i=1,dim
1097                STIFF( (p-1)*c+i,(q-1)*c+i ) = &
1098                    STIFF( (p-1)*c+i,(q-1)*c+i ) + &
1099                    s * SlipCoeff(i) * Basis(q) * Basis(p)
1100              END DO
1101            END DO
1102          END DO
1103        END IF
1104      END IF
1105
1106      ! Assemble given forces to r.h.s.
1107      IF( HaveForce .OR. HavePres ) THEN
1108        IF ( NormalTangential ) THEN
1109          DO i=1,dim
1110            SELECT CASE(i)
1111            CASE(1)
1112              Vect = Normal
1113            CASE(2)
1114              Vect = Tangent
1115            CASE(3)
1116              Vect = Tangent2
1117            END SELECT
1118
1119            DO q=1,nd
1120              DO j=1,dim
1121                l = (q-1)*c + j
1122                FORCE(l) = FORCE(l) + s * Basis(q) * SurfaceTraction(i) * Vect(j)
1123              END DO
1124            END DO
1125          END DO
1126        ELSE
1127          DO i=1,dim
1128            DO q=1,nd
1129              k = (q-1)*c + i
1130              FORCE(k) = FORCE(k) + s * Basis(q) * SurfaceTraction(i)
1131            END DO
1132          END DO
1133        END IF
1134      END IF
1135    END DO
1136
1137    CALL DefaultUpdateEquations( STIFF, FORCE )
1138
1139  END SUBROUTINE LocalBoundaryMatrix
1140
1141
1142
1143
1144END MODULE IncompressibleLocalForms
1145
1146
1147!------------------------------------------------------------------------------
1148SUBROUTINE IncompressibleNSSolver_Init0(Model, Solver, dt, Transient)
1149!------------------------------------------------------------------------------
1150  USE DefUtils
1151  IMPLICIT NONE
1152!------------------------------------------------------------------------------
1153  TYPE(Model_t) :: Model
1154  TYPE(Solver_t) :: Solver
1155  REAL(KIND=dp) :: dt
1156  LOGICAL :: Transient
1157!------------------------------------------------------------------------------
1158  CALL ListAddNewString(GetSolverParams(),'Element', &
1159    'p:1 -tri b:1 -tetra b:1 -quad b:3 -brick b:4 -prism b:4 -pyramid b:4')
1160!------------------------------------------------------------------------------
1161END SUBROUTINE IncompressibleNSSolver_Init0
1162!------------------------------------------------------------------------------
1163
1164
1165!------------------------------------------------------------------------------
1166SUBROUTINE IncompressibleNSSolver_init(Model, Solver, dt, Transient)
1167!------------------------------------------------------------------------------
1168  USE DefUtils
1169  IMPLICIT NONE
1170!------------------------------------------------------------------------------
1171  TYPE(Model_t) :: Model
1172  TYPE(Solver_t) :: Solver
1173  REAL(KIND=dp) :: dt
1174  LOGICAL :: Transient
1175!------------------------------------------------------------------------------
1176  TYPE(ValueList_t), POINTER :: Params
1177  LOGICAL :: Found
1178  INTEGER :: dim
1179  CHARACTER(*), PARAMETER :: Caller = 'IncompressibleNSSolver_init'
1180!------------------------------------------------------------------------------
1181  Params => GetSolverParams()
1182
1183  IF( ListCheckPresentAnyBC( Model, 'Pressure 1' ) ) THEN
1184    CALL Fatal( Caller,'Use >Surface Traction 1< instead of >Pressure 1<')
1185  END IF
1186  IF( ListCheckPresentAnyBC( Model, 'Pressure 2' ) ) THEN
1187    CALL Fatal( Caller,'Use >Surface Traction 3< instead of >Pressure 2<')
1188  END IF
1189  IF( ListCheckPresentAnyBC( Model, 'Pressure 3' ) ) THEN
1190    CALL Fatal( Caller,'Use >Surface Traction 3< instead of >Pressure 3<')
1191  END IF
1192
1193  dim = CoordinateSystemDimension()
1194
1195  IF ( dim == 2 ) THEN
1196    CALL ListAddNewString(Params, 'Variable', &
1197        'Flow Solution[Velocity:2 Pressure:1]')
1198  ELSE
1199    CALL ListAddNewString(Params, 'Variable', &
1200        'Flow Solution[Velocity:3 Pressure:1]')
1201  END IF
1202
1203  ! Study only velocity components in linear system
1204  CALL ListAddNewInteger(Params, 'Nonlinear System Norm DOFs', dim )
1205
1206  ! This should be true to incompressible flows where pressure level is not uniquely determined
1207  CALL ListAddNewLogical(Params, 'Relative Pressure Relaxation', .TRUE. )
1208
1209  ! Automate the choice for the variational formulation:
1210  CALL ListAddNewLogical(Params, 'GradP Discretization', .FALSE.)
1211  CALL ListAddNewLogical(Params, 'Div-Curl Discretization', .FALSE.)
1212
1213
1214  ! It makes sense to eliminate the bubbles to save memory and time
1215  CALL ListAddNewLogical(Params, 'Bubbles in Global System', .FALSE.)
1216
1217  ! The recovery of transient bubble DOFs is done such that at least two
1218  ! iterations within the same time step are needed. The multiple solutions are
1219  ! ensured by making at least two nonlinear iterations. However, if "steady
1220  ! state" iterations are performed, computing just one nonlinear iterate
1221  ! is sufficient.
1222  IF ( .NOT. ListGetLogical(Params, 'Bubbles In Global System', Found) ) THEN
1223    CALL ListAddNewInteger(Params, 'Nonlinear System Min Iterations', 2)
1224  END IF
1225
1226  ! Create solver related to variable "iterpres" when using block preconditioning
1227  ! There keyword ensure that the matrix is truly used in the library version of the
1228  ! block solver.
1229  IF( ListGetLogical( Params,'Block Preconditioner',Found ) ) THEN
1230    CALL ListAddNewString( Params,'Block Matrix Schur Variable','schur')
1231  END IF
1232
1233!------------------------------------------------------------------------------
1234END SUBROUTINE IncompressibleNSSolver_Init
1235!------------------------------------------------------------------------------
1236
1237!------------------------------------------------------------------------------
1238SUBROUTINE IncompressibleNSSolver(Model, Solver, dt, Transient)
1239!------------------------------------------------------------------------------
1240  USE DefUtils
1241  USE IncompressibleLocalForms
1242  USE MainUtils
1243
1244
1245  IMPLICIT NONE
1246!------------------------------------------------------------------------------
1247  TYPE(Solver_t) :: Solver
1248  TYPE(Model_t) :: Model
1249  REAL(KIND=dp) :: dt
1250  LOGICAL :: Transient
1251!------------------------------------------------------------------------------
1252! Local variables
1253!------------------------------------------------------------------------------
1254  TYPE(Element_t), POINTER :: Element
1255  TYPE(ValueList_t), POINTER :: Params
1256  TYPE(Mesh_t), POINTER :: Mesh
1257  TYPE(GaussIntegrationPoints_t) :: IP
1258
1259  INTEGER :: Element_id
1260  INTEGER :: i, n, nb, nd, nbdofs, dim, Active, maxiter, iter
1261  INTEGER :: stimestep = -1
1262
1263  REAL(KIND=dp) :: Norm
1264
1265  LOGICAL :: AllocationsDone = .FALSE., Found, StokesFlow, BlockPrec
1266  LOGICAL :: GradPVersion, DivCurlForm, SpecificLoad, InitHandles,InitBCHandles
1267
1268  TYPE(Solver_t), POINTER, SAVE :: SchurSolver => Null()
1269
1270  CHARACTER(*), PARAMETER :: Caller = 'IncompressibleNSSolver'
1271
1272  SAVE AllocationsDone, stimestep
1273
1274!------------------------------------------------------------------------------
1275! Local variables to be accessed by the contained subroutines:
1276!------------------------------------------------------------------------------
1277  LOGICAL :: LinearAssembly, Newton
1278!------------------------------------------------------------------------------
1279
1280  CALL DefaultStart()
1281
1282  dim = CoordinateSystemDimension()
1283  Mesh => GetMesh()
1284
1285  !-----------------------------------------------------------------------------
1286  ! Allocate some permanent storage, this is done first time only:
1287  !-----------------------------------------------------------------------------
1288  IF (.NOT. AllocationsDone .AND. Transient .AND. &
1289      Mesh % MaxBDOFs > 0) THEN
1290    !
1291    ! Allocate arrays having a sufficient size for listing all bubble entries of
1292    ! the velocity solution (current and previous). These are needed in order to
1293    ! evaluate the time derivative of the bubble part.
1294    !
1295    nbdofs = Mesh % MaxBDOFs*dim*GetNOFActive()
1296    ALLOCATE(bx(nbdofs), bxprev(nbdofs));
1297    bx=0.0_dp; bxprev=0.0_dp
1298
1299    AllocationsDone = .TRUE.
1300  END IF
1301
1302  ! Check if the previous bubble part (bxprev) needs to be updated:
1303  IF (Transient .AND. GetTimestep() /= stimestep .AND. &
1304      Mesh % MaxBDOFs > 0) THEN
1305    bxprev = bx
1306    stimestep = GetTimestep()
1307  END IF
1308
1309  Params => GetSolverParams()
1310
1311  !-----------------------------------------------------------------------------
1312  ! Output the number of integration points as information.
1313  ! This in not fully informative if several element types are present.
1314  !-----------------------------------------------------------------------------
1315  Element => Mesh % Elements( Solver % ActiveElements(1) )
1316  IP = GaussPointsAdapt( Element, PReferenceElement = .TRUE. )
1317  CALL Info('IncompressibleNSSolver', &
1318      'Number of 1st integration points: '//TRIM(I2S(IP % n)), Level=5)
1319
1320  !-----------------------------------------------------------------------------
1321  ! Set the flags/parameters which define how the system is assembled:
1322  !-----------------------------------------------------------------------------
1323  LinearAssembly = GetLogical(Params, 'Linear Equation', Found )
1324  StokesFlow = GetLogical(Params, 'Stokes Flow', Found )
1325  GradPVersion = GetLogical(Params, 'GradP Discretization', Found)
1326  DivCurlForm = GetLogical(Params, 'Div-Curl Discretization', Found)
1327  BlockPrec = GetLogical(Params,'Block Preconditioner',Found )
1328  SpecificLoad = GetLogical(Params,'Specific Load',Found)
1329
1330  Maxiter = GetInteger(Params, 'Nonlinear system max iterations', Found)
1331  IF (.NOT.Found) Maxiter = 1
1332  !-----------------------------------------------------------------------------
1333
1334  IF (DivCurlForm) CALL Info(Caller, 'The div-curl form is used for the viscous terms')
1335  IF (GradPVersion) CALL Info(Caller, 'The pressure gradient is not integrated by parts')
1336  IF (BlockPrec) CALL Info(Caller,'Creating pressure block for block preconditioner')
1337
1338  IF(BlockPrec ) THEN
1339    ! Create solver that only acts as a container for the shcur complement
1340    ! matrix used in the block preconditioning solver of the library.
1341    IF( .NOT. ASSOCIATED( SchurSolver ) ) THEN
1342      SchurSolver => CreateChildSolver( Solver,'schur', 1 )
1343    END IF
1344  END IF
1345
1346
1347
1348  DO iter=1,maxiter
1349
1350    CALL Info(Caller,'--------------------------------------------------------', Level=4)
1351    WRITE( Message,'(A,I4)') 'Nonlinear iteration:', Iter
1352    CALL Info(Caller, Message, Level=4)
1353    CALL Info(Caller,'--------------------------------------------------------', Level=4)
1354
1355    Active = GetNOFActive()
1356    CALL DefaultInitialize()
1357    IF (ASSOCIATED(SchurSolver)) THEN
1358      CALL DefaultInitialize(USolver=SchurSolver)
1359    END IF
1360
1361    Newton = GetNewtonActive( Solver )
1362
1363    DO Element_id=1,1
1364      Element => GetActiveElement(Element_id)
1365      n  = GetElementNOFNodes(Element)
1366      !
1367      ! When the number of bubbles is obtained with the Update=.TRUE. flag,
1368      ! we need to call GetElementNOFBDOFs before calling GetElementNOFDOFs.
1369      !
1370      nb = GetElementNOFBDOFs(Element, Update=.TRUE.)
1371      nd = GetElementNOFDOFs(Element)
1372
1373      ! Get element local matrix and rhs vector:
1374      !-----------------------------------------
1375      CALL LocalBulkMatrix(Element, n, nd, nd+nb, dim,  DivCurlForm, GradPVersion, &
1376          SpecificLoad, StokesFlow, dt, LinearAssembly, nb, Newton, Transient, .TRUE., &
1377          SchurSolver )
1378    END DO
1379
1380    !$OMP PARALLEL SHARED(Active, dim, SpecificLoad, StokesFlow, &
1381    !$OMP                 DivCurlForm, GradPVersion, &
1382    !$OMP                 dt, LinearAssembly, Newton, Transient, SchurSolver ) &
1383    !$OMP PRIVATE(Element, Element_id, n, nd, nb)  DEFAULT(None)
1384    !$OMP DO
1385    DO Element_id=2,Active
1386      Element => GetActiveElement(Element_id)
1387      n  = GetElementNOFNodes(Element)
1388      nb = GetElementNOFBDOFs(Element, Update=.TRUE.)
1389      nd = GetElementNOFDOFs(Element)
1390
1391      ! Get element local matrix and rhs vector:
1392      !-----------------------------------------
1393      CALL LocalBulkMatrix(Element, n, nd, nd+nb, dim,  DivCurlForm, GradPVersion, &
1394          SpecificLoad, StokesFlow, dt, LinearAssembly, nb, Newton, Transient, .FALSE.,&
1395          SchurSolver )
1396    END DO
1397    !$OMP END DO
1398    !$OMP END PARALLEL
1399
1400    CALL DefaultFinishBulkAssembly()
1401
1402    Active = GetNOFBoundaryElements()
1403    InitBCHandles = .TRUE.
1404    DO Element_id=1,Active
1405      Element => GetBoundaryElement(Element_id)
1406      IF (ActiveBoundaryElement()) THEN
1407        n  = GetElementNOFNodes()
1408        nd = GetElementNOFDOFs()
1409
1410        IF ( GetElementFamily() == 1 ) CYCLE
1411
1412        ! Get element local matrix and rhs vector:
1413        !-----------------------------------------
1414        CALL LocalBoundaryMatrix(Element, n, nd, dim, InitBCHandles )
1415        InitBCHandles = .FALSE.
1416      END IF
1417    END DO
1418
1419    CALL DefaultFinishBoundaryAssembly()
1420
1421    CALL DefaultFinishAssembly()
1422    CALL DefaultDirichletBCs()
1423    IF(ASSOCIATED(SchurSolver)) CALL DefaultDirichletBCs(USolver=SchurSolver)
1424
1425    Norm = DefaultSolve()
1426
1427    IF ( Solver % Variable % NonlinConverged == 1 ) EXIT
1428  END DO
1429
1430  CALL DefaultFinish()
1431
1432  CALL Info( Caller,'All done',Level=10)
1433!------------------------------------------------------------------------------
1434END SUBROUTINE IncompressibleNSSolver
1435!------------------------------------------------------------------------------
1436