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