1!/*****************************************************************************/ 2! * 3! * Elmer, A Finite Element Software for Multiphysical Problems 4! * 5! * Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland 6! * 7! * This program is free software; you can redistribute it and/or 8! * modify it under the terms of the GNU General Public License 9! * as published by the Free Software Foundation; either version 2 10! * of the License, or (at your option) any later version. 11! * 12! * This program is distributed in the hope that it will be useful, 13! * but WITHOUT ANY WARRANTY; without even the implied warranty of 14! * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15! * GNU General Public License for more details. 16! * 17! * You should have received a copy of the GNU General Public License 18! * along with this program (in file fem/GPL-2); if not, write to the 19! * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 20! * Boston, MA 02110-1301, USA. 21! * 22! *****************************************************************************/ 23! 24!/****************************************************************************** 25! * 26! * Authors: Juha Ruokolainen, Mika Malinen 27! * Email: mika.malinen@csc.fi & Juha.Ruokolainen@csc.fi 28! * Web: http://www.csc.fi/elmer 29! * Address: CSC - IT Center for Science Ltd. 30! * Keilaranta 14 31! * 02101 Espoo, Finland 32! * 33! * Original Date: 2012-01-30 34! * 35! *****************************************************************************/ 36 37!------------------------------------------------------------------------------ 38SUBROUTINE StokesSolver_Init0(Model, Solver, dt, Transient) 39!------------------------------------------------------------------------------ 40 USE DefUtils 41 USE SolverUtils 42 USE ElementUtils 43 44 IMPLICIT NONE 45!------------------------------------------------------------------------------ 46 TYPE(Solver_t) :: Solver 47 TYPE(Model_t) :: Model 48 REAL(KIND=dp) :: dt 49 LOGICAL :: Transient 50!------------------------------------------------------------------------------ 51 TYPE(ValueList_t), POINTER :: SolverParams 52!------------------------------------------------------------------------------ 53 SolverParams => GetSolverParams() 54 55 IF ( .NOT. ListCheckPresent(SolverParams, 'Bubbles in Global System') ) & 56 CALL ListAddLogical(SolverParams, 'Bubbles in Global System', .FALSE.) 57 58 IF ( .NOT. ListCheckPresent(SolverParams, 'Linear System Solver') ) & 59 CALL ListAddString(SolverParams, 'Linear System Solver', 'Iterative') 60 IF ( .NOT. ListCheckPresent(SolverParams, 'Linear System Iterative Method') ) & 61 CALL ListAddString(SolverParams, 'Linear System Iterative Method', 'GCR') 62 IF ( .NOT. ListCheckPresent(SolverParams, 'Linear System GCR Restart') ) & 63 CALL ListAddInteger(SolverParams, 'Linear System GCR Restart', 50) 64 IF ( .NOT. ListCheckPresent(SolverParams, 'Linear System Max Iterations') ) & 65 CALL ListAddInteger(SolverParams, 'Linear System Max Iterations', 200) 66 IF ( .NOT. ListCheckPresent(SolverParams, 'Linear System Row Equilibration') ) & 67 CALL ListAddLogical(SolverParams, 'Linear System Row Equilibration', .TRUE.) 68 IF ( .NOT. ListCheckPresent(SolverParams, 'Linear System Convergence Tolerance') ) & 69 CALL ListAddConstReal(SolverParams, 'Linear System Convergence Tolerance', 1.0d-6) 70 IF ( .NOT. ListCheckPresent(SolverParams, 'Linear System Base Tolerance') ) & 71 CALL ListAddConstReal(SolverParams, 'Linear System Base Tolerance', 1.0d-3) 72 IF ( .NOT. ListCheckPresent(SolverParams, 'Linear System Relative Tolerance') ) & 73 CALL ListAddConstReal(SolverParams, 'Linear System Relative Tolerance', 1.0d-2) 74 75!------------------------------------------------------------------------------ 76END SUBROUTINE StokesSolver_Init0 77!------------------------------------------------------------------------------ 78 79 80!------------------------------------------------------------------------------ 81SUBROUTINE StokesSolver( Model,Solver,dt,TransientSimulation ) 82!------------------------------------------------------------------------------ 83!****************************************************************************** 84! 85! A parallel solver that uses two-level iterations to solve the discrete 86! Stokes model. Inner iterations can be associated with preconditioning and 87! they provide search directions for the outer iterative method (GCR) applied 88! to the primitive Stokes problem. 89! 90! A key design choice here has been that the inner iterations are performed 91! via calling DefaultSolve routine, so that the full range of standard parallel 92! methods can be applied for this computation and a clean interface is obtained 93! with this piece of code. The outer iterative method is however implemented 94! locally into this solver. 95! 96! ARGUMENTS: 97! 98! TYPE(Model_t) :: Model, 99! INPUT: All model information (mesh, materials, BCs, etc...) 100! 101! TYPE(Solver_t) :: Solver 102! INPUT: Linear & nonlinear equation solver options 103! 104! REAL(KIND=dp) :: dt, 105! INPUT: Timestep size for time dependent simulations 106! 107! LOGICAL :: TransientSimulation 108! INPUT: Steady state or transient simulation 109! 110!****************************************************************************** 111 USE DefUtils 112 USE SolverUtils 113 USE ElementUtils 114 USE MaterialModels 115 USE ElementDescription, ONLY: GetEdgeMap 116 117 IMPLICIT NONE 118!------------------------------------------------------------------------------ 119 TYPE(Solver_t), TARGET :: Solver 120 TYPE(Model_t) :: Model 121 REAL(KIND=dp) :: dt 122 LOGICAL :: TransientSimulation 123!------------------------------------------------------------------------------ 124! Local variables 125!------------------------------------------------------------------------------ 126 LOGICAL :: AllocationsDone = .FALSE., Found, Convect, GotForceBC, NormalTangential, & 127 SkipPowerLaw=.TRUE., Newton = .FALSE. 128 TYPE(Element_t), POINTER :: Element 129 130 REAL(KIND=dp) :: Norm, NonlinearTol, NewtonThreshold, NonLinError, res 131 INTEGER, POINTER :: EdgeMap(:,:) 132 INTEGER :: BrickFaceMap(6,4) 133 INTEGER :: n, m, p, q, nb, nd, t, istat, active, dim, & 134 iter, NonlinearIter, MaxPicardIterations, MidEdgeNodes(12) 135 TYPE(Mesh_t), POINTER :: Mesh 136 TYPE(ValueList_t), POINTER :: BodyForce, Material, BC 137 REAL(KIND=dp), ALLOCATABLE :: STIFF(:,:), LOAD(:,:), FORCE(:), rho(:), mu(:), & 138 Velocity(:,:), LoadVector(:,:), SlipCoeff(:,:), ExtPressure(:) 139 140 TYPE(Nodes_t) :: ElementNodes 141 CHARACTER(LEN=MAX_NAME_LEN) :: NormalTangentialName 142 143 SAVE STIFF, LOAD, FORCE, rho, mu, Velocity, AllocationsDone, ElementNodes, & 144 LoadVector, SlipCoeff, ExtPressure 145 146 147!----------------------------------------------------------------------------- 148! Variables for two-level iterations... 149!----------------------------------------------------------------------------- 150 LOGICAL :: BlockPreconditioning, Parallel, UpdateMatrix, UseTrueResidual, Timing, P2P1 151 LOGICAL :: DoScaling, DoEquilibration, BlockDiagonalA, UseVeloLaplacian, AdaptiveTols 152 CHARACTER(LEN=MAX_NAME_LEN) :: Eq 153 INTEGER :: i, j, k, nlen, Round, MaxIterations, RestartM 154 INTEGER, ALLOCATABLE :: Indexes(:) 155 INTEGER, POINTER :: NodeIndexes(:) 156 REAL(KIND=dp) :: MinTolerance, t0, rt0, st, rst, ct, & 157 BaseTolerance, TargetTol, RelTolerance, PrecondTol 158 TYPE(Solver_t), POINTER :: PressureSolver, VelocitySolver 159 TYPE(Matrix_t), POINTER :: MMatrix, PMatrix, AMatrix 160 REAL(KIND=dp), ALLOCATABLE :: Snew(:), R(:), S(:,:), V(:,:), & 161 ALocal(:,:), PLocal(:,:), Vx(:), Vy(:), Vz(:) 162 REAL(KIND=dp), ALLOCATABLE, TARGET :: Residual(:), RotatedVarValues(:), TempRes(:), TempS(:), & 163 TempMx(:), TempMb(:), TempMr(:), TempRHS(:) 164 REAL(KIND=dp), POINTER :: Mx(:), Mb(:), Mr(:) 165 SAVE Snew, R, S, V, Indexes, ALocal, PLocal, Vx, Vy, Vz 166 REAL(KIND=dp), ALLOCATABLE :: RowScalingA(:), RowScalingP(:) 167 SAVE RowScalingA, RowScalingP 168 SAVE MaxIterations,RestartM 169!------------------------------------------------------------------------------ 170 NonLinError=2.0_dp 171 Solver % Variable % NonLinChange = NonLinError 172 173 Parallel = ParEnv % Pes > 1 174 175 dim = CoordinateSystemDimension() 176 177 P2P1 = ListGetLogical( Solver % Values, 'P2-P1 Approximation', Found) 178 BlockPreconditioning = ListGetLogical( Solver % Values, 'Block Preconditioning', Found) 179 IF ( .NOT. Found) BlockPreconditioning = .FALSE. 180 181 Convect = GetLogical( GetSolverParams(), 'Convective', Found ) 182 IF ( .NOT. Found ) Convect = .FALSE. 183 184 SkipPowerLaw = ListGetLogical( Solver % Values, 'Constant-Viscosity Start', Found) 185 IF ( .NOT. Found) SkipPowerLaw = .TRUE. 186 187 DoScaling=ListGetLogical(Solver % Values, 'Linear System Scaling', Found) 188 DoEquilibration=ListGetLogical(Solver % Values, 'Linear System Row Equilibration', Found) 189 DoScaling = DoScaling .OR. DoEquilibration 190 191 BlockDiagonalA = ListGetLogical( Solver % Values, 'Block Diagonal A', Found) 192 UseVeloLaplacian = ListGetLogical( Solver % Values, 'Use Velocity Laplacian', Found) 193 194 AdaptiveTols = ListGetLogical(Solver % Values, 'Linear System Adaptive Tolerance', Found) 195 IF (AdaptiveTols) THEN 196 TargetTol = ListGetConstReal(Solver % Values, 'Linear System Convergence Tolerance') 197 IF (.NOT. ListCheckPresent(Solver % Values, 'Linear System Relative Tolerance')) THEN 198 RelTolerance = 1.0d-2 199 ELSE 200 RelTolerance = ListGetConstReal(Solver % Values, 'Linear System Relative Tolerance') 201 END IF 202 203 IF (.NOT. ListCheckPresent(Solver % Values, 'Linear System Base Tolerance')) THEN 204 BaseTolerance = 1.0d-4 205 ELSE 206 BaseTolerance = ListGetConstReal(Solver % Values, 'Linear System Base Tolerance') 207 END IF 208 END IF 209 210 !---------------------------------------------------------------------------------- 211 ! Find the preconditioning matrices and solvers generated before simulation 212 !---------------------------------------------------------------------------------- 213 IF (BlockPreconditioning) THEN 214 DO i=1,Model % NumberOfSolvers 215 Eq = ListGetString( Model % Solvers(i) % Values, 'Equation', Found ) 216 IF (Found) THEN 217 nlen = LEN_TRIM(eq) 218 SELECT CASE( Eq(1:nlen) ) 219 CASE ('pressure preconditioning') 220 PressureSolver => Model % Solvers(i) 221 PMatrix => PressureSolver % Matrix 222 CASE ('velocity preconditioning') 223 VelocitySolver => Model % Solvers(i) 224 AMatrix => VelocitySolver % Matrix 225 END SELECT 226 END IF 227 END DO 228 229 IF ( .NOT. ASSOCIATED(PressureSolver) ) & 230 CALL Fatal( 'ParStokes', 'Pressure preconditioning operation missing...' ) 231 232 IF ( .NOT. ASSOCIATED(VelocitySolver) ) & 233 CALL Fatal( 'ParStokes', 'Velocity preconditioning operation missing...' ) 234 235 !--------------------------------------------------------------------------- 236 ! Perform permutation check... 237 !--------------------------------------------------------------------------- 238 n = SIZE(Solver % Variable % Perm) 239 IF ( ANY(Solver % Variable % Perm(1:n) /= PressureSolver % Variable % Perm(1:n)) ) & 240 CALL Fatal( 'ParStokes', & 241 'Nonmatching variable permutations, use Optimize Bandwidth in pressure preconditioning' ) 242 IF ( ANY(Solver % Variable % Perm(1:n) /= VelocitySolver % Variable % Perm(1:n)) ) & 243 CALL Fatal( 'ParStokes',& 244 'Nonmatching variable permutations, use Optimize Bandwidth in velocity preconditioning' ) 245 END IF 246 247 248 !Allocate some permanent storage, this is done first time only: 249 !-------------------------------------------------------------- 250 Mesh => GetMesh() 251 252 IF ( .NOT. AllocationsDone ) THEN 253 p = Mesh % MaxElementDOFs ! just big enough for elemental arrays 254 n = p * (dim+1) 255 m = p * dim 256 ALLOCATE( & 257 FORCE(n), & 258 LOAD(p,dim+1), & 259 STIFF(n,n), & 260 Indexes(p), & 261 rho(p), & 262 mu(p), & 263 Velocity(dim,p), & 264 Vx(p), & 265 Vy(p), & 266 Vz(p), & 267 LoadVector(4,p), & 268 SlipCoeff(3,p), & 269 ExtPressure(p), & 270 ALocal(m,m), & 271 PLocal(p,p), & 272 STAT=istat ) 273 IF ( istat /= 0 ) THEN 274 CALL Fatal( 'ParStokes', 'Memory allocation error.' ) 275 END IF 276 277 IF (BlockPreconditioning) THEN 278 MaxIterations = ListGetInteger(Solver % Values, 'Linear System Max Iterations') 279 ReStartM = ListGetInteger(Solver % Values, 'Linear System GCR Restart', Found) 280 IF ( .NOT. Found) THEN 281 ReStartM = MaxIterations 282 END IF 283 ALLOCATE( & 284 Snew( Solver % Matrix % NumberOfRows ), & ! Search direction variable 285 R( Solver % Matrix % NumberOfRows ), & ! For saving the residual generated by the iterative method 286 S( Solver % Matrix % NumberOfRows, RestartM), & ! For saving search directions 287 V( Solver % Matrix % NumberOfRows, RestartM), & ! For saving the range of coefficient matrix 288 STAT=istat ) 289 IF ( istat /= 0 ) THEN 290 CALL Fatal( 'ParStokes', 'Memory allocation error.' ) 291 END IF 292 END IF 293 294 AllocationsDone = .TRUE. 295 END IF 296 297 !------------------------------------------------------------------- 298 ! System assembly & solve for each system arising from linearization 299 !-------------------------------------------------------------------- 300 NonlinearIter = ListGetInteger( Solver % Values, & 301 'Nonlinear System Max Iterations', minv=0 ) 302 NonlinearTol = ListGetConstReal( Solver % Values, & 303 'Nonlinear System Convergence Tolerance', minv=0.0d0 ) 304 305 IF (BlockPreconditioning) THEN 306 !Solver % Variable % Values = 0.0d0 !Fab: No need to reset to 0 307 Solver % Matrix % RHS = 0.0d0 308 NonlinearIter = NonlinearIter + 1 ! To enable the computation of true nonlinear residual 309 END IF 310 311 Active = GetNOFActive() 312 313 DO iter=1, NonlinearIter 314 !----------------------------------------------------------------------- 315 ! Check whether it is the time to switch to Newton linearization... 316 !----------------------------------------------------------------------- 317 Newton = .FALSE. 318 IF (iter > 2) THEN 319 NewtonThreshold = ListGetConstReal( Solver % Values, & 320 'Nonlinear System Newton After Tolerance', Found ) 321 IF ( Found .AND. (NonLinError < NewtonThreshold) ) THEN 322 Newton = .TRUE. 323 ELSE 324 MaxPicardIterations = ListGetInteger( Solver % Values, & 325 'Nonlinear System Newton After Iterations', Found ) 326 IF ( Found ) THEN 327 IF ( iter > MaxPicardIterations ) THEN 328 Newton = .TRUE. 329 END IF 330 END IF 331 END IF 332 END IF 333 334 IF (Newton) THEN 335 WRITE(Message,'(a)') 'The Newton linearization is used...' 336 CALL Info('ParSolver', Message, Level=4) 337 END IF 338 339 !------------------------------------------------------------------ 340 ! Initialize matrix structures... 341 !------------------------------------------------------------------ 342 CALL DefaultInitialize() 343 344 IF (BlockPreconditioning) THEN 345 !CALL InitializeToZero( AMatrix, AMatrix % RHS ) 346 !CALL InitializeToZero( PMatrix, PMatrix % RHS ) 347 CALL DefaultInitialize(USolver=VelocitySolver) 348 CALL DefaultInitialize(USolver=PressureSolver) 349 END IF 350 351 !------------------------------------------------------------ 352 ! We need to make a copy of the rotated solution vector to 353 ! handle normal-tangential bcs. 354 !------------------------------------------------------------ 355 IF ( BlockPreconditioning .AND. Iter==1 ) THEN 356 n = Solver % Matrix % NumberOfRows 357 ALLOCATE( RotatedVarValues(n) ) 358 RotatedVarValues = 0.0d0 359 END IF 360 361 !--------------------------------------------------------------- 362 ! We always generate the initial guess by adjusting material law 363 ! exponent to obtain linear behaviour 364 !--------------------------------------------------------------- 365 IF ((Iter==1) .AND. (SkipPowerlaw)) THEN !SkipPowerlaw only 1st time 366 SkipPowerlaw = .TRUE. 367 ELSE 368 SkipPowerlaw = .FALSE. 369 END IF 370 371 CALL StartAdvanceOutput( 'StokesSolver', 'Assembly:' ) 372 DO t=1,Active 373 CALL AdvanceOutput(t, Active) 374 375 Element => GetActiveElement(t) 376 n = GetElementNOFNodes() 377 nd = GetElementDOFs( Indexes ) 378 nb = GetElementNOFBDOFs() 379 380 !----------------------------------------------- 381 ! Body forces: 382 !----------------------------------------------- 383 LOAD = 0.0d0 384 BodyForce => GetBodyForce() 385 IF ( ASSOCIATED(BodyForce) ) THEN 386 Load(1:n,1) = GetReal( BodyForce, 'Flow BodyForce 1', Found ) 387 Load(1:n,2) = GetReal( BodyForce, 'Flow BodyForce 2', Found ) 388 IF (dim > 2) & 389 Load(1:n,3) = GetReal( BodyForce, 'Flow BodyForce 3', Found ) 390 END IF 391 392 !----------------------------------------------- 393 ! Material parameters: 394 !----------------------------------------------- 395 Material => GetMaterial() 396 rho(1:n) = GetReal( Material, 'Density' ) 397 mu(1:n) = GetReal( Material, 'Viscosity' ) 398 !------------------------------------------------------- 399 400 !-------------------------------------------------------------------------- 401 ! Get previous elementwise velocity iterate for approximating nonlinearity 402 !-------------------------------------------------------------------------- 403 Vx = 0.0d0 404 Vy = 0.0d0 405 Vz = 0.0d0 406 CALL GetScalarLocalSolution( Vx, ComponentName(Solver % Variable % Name,1) ) 407 CALL GetScalarLocalSolution( Vy, ComponentName(Solver % Variable % Name,2) ) 408 IF( dim > 2 ) THEN 409 CALL GetScalarLocalSolution( Vz, ComponentName(Solver % Variable % Name,3) ) 410 END IF 411 412 !--------------------------------------------------------------------- 413 ! Get element local matrix and rhs vector: 414 !--------------------------------------------------------------------- 415 CALL LocalMatrix( STIFF, ALocal, PLocal, FORCE, LOAD, rho, mu, & 416 Vx, Vy, Vz, Element, n, nd+nb, dim, Convect, P2P1, SkipPowerLaw, Newton, & 417 BlockDiagonalA) 418 419 IF ( nb > 0 ) THEN 420 CALL StaticCondensation( nd, nb, dim, STIFF, FORCE ) 421 END IF 422 423 IF( BlockPreconditioning ) THEN 424 IF (.NOT. BlockDiagonalA) THEN 425 Alocal = 0.0d0 426 DO p=1,nd 427 DO i=1,dim 428 DO q=1,nd 429 DO j=1,dim 430 Alocal( (p-1)*dim+i, (q-1)*dim+j ) = & 431 stiff( (p-1)*(dim+1)+i, (q-1)*(dim+1)+j) 432 END DO 433 END DO 434 END DO 435 END DO 436 ELSE 437 ! The block diagonal approximation of the velocity preconditioning matrix: 438 IF (.NOT. UseVeloLaplacian) THEN 439 Alocal = 0.0d0 440 DO p=1,nd 441 DO i=1,dim 442 DO q=1,nd 443 Alocal( (p-1)*dim+i, (q-1)*dim+i ) = & 444 stiff( (p-1)*(dim+1)+i, (q-1)*(dim+1)+i) 445 END DO 446 END DO 447 END DO 448 END IF 449 END IF 450 END IF 451 452 CALL DefaultUpdateEquations( STIFF, FORCE ) 453 454 IF (BlockPreconditioning) THEN 455 Force = 0.0d0 456 CALL DefaultUpdateEquations( ALocal, FORCE, USolver=VelocitySolver) 457 CALL DefaultUpdateEquations( PLocal, FORCE, USolver=PressureSolver) 458 END IF 459 460 END DO 461 462 CALL DefaultFinishBulkAssembly() 463 464 !------------------------------------------------------------------------------ 465 ! Surface force and slip boundary conditions 466 !------------------------------------------------------------------------------ 467 DO t = 1, GetNOFBoundaryElements() 468 469 Element => GetBoundaryElement(t) 470 IF ( .NOT. ActiveBoundaryElement() ) CYCLE 471 472 n = GetElementNOFNodes() 473 nd = GetElementDOFs( Indexes ) 474 475 IF ( GetElementFamily() == 1 ) CYCLE 476 477 CALL GetElementNodes( ElementNodes ) 478 NodeIndexes => Element % NodeIndexes 479 480 BC => GetBC() 481 IF ( .NOT. ASSOCIATED(BC) ) CYCLE 482 483 GotForceBC = GetLogical( BC, 'Flow Force BC', Found ) 484 IF ( .NOT. Found ) GotForceBC = .TRUE. 485 486 IF ( GotForceBC ) THEN 487 STIFF = 0.0d0 488 FORCE = 0.0d0 489 490 ExtPressure = 0.0d0 491 ExtPressure(1:n) = GetReal( BC, 'Normal Surface Traction', GotForceBC ) 492 IF (.NOT. GotForceBC) ExtPressure(1:n) = GetReal( BC, 'External Pressure', GotForceBC ) 493 494 !------------------------------------------------------------------------------ 495 ! Surface force in given direction BC: \sigma n = F 496 !------------------------------------------------------------------------------ 497 LoadVector = 0.0d0 498 LoadVector(1,1:n) = GetReal( BC, 'Surface Traction 1', Found ) 499 IF (.NOT. Found) LoadVector(1,1:n) = GetReal( BC, 'Pressure 1', Found ) 500 LoadVector(2,1:n) = GetReal( BC, 'Surface Traction 2', Found ) 501 IF (.NOT. Found) LoadVector(2,1:n) = GetReal( BC, 'Pressure 2', Found ) 502 LoadVector(3,1:n) = GetReal( BC, 'Surface Traction 3', Found ) 503 IF (.NOT. Found) LoadVector(3,1:n) = GetReal( BC, 'Pressure 3', Found ) 504 LoadVector(4,1:n) = 0.0d0 505 506 !------------------------------------------------------------------------------ 507 ! Slip boundary condition BC: \sigma n = R_k u_k 508 !------------------------------------------------------------------------------ 509 SlipCoeff = 0.0d0 510 SlipCoeff(1,1:n) = GetReal( BC, 'Slip Coefficient 1', Found ) 511 SlipCoeff(2,1:n) = GetReal( BC, 'Slip Coefficient 2', Found ) 512 SlipCoeff(3,1:n) = GetReal( BC, 'Slip Coefficient 3', Found ) 513 514 515 NormalTangentialName = 'Normal-Tangential' 516 NormalTangentialName = TRIM(NormalTangentialName) // ' ' // & 517 GetVarName(Solver % Variable) 518 519 NormalTangential = GetLogical( BC, & 520 NormalTangentialName, Found ) 521 522 CALL NavierStokesBoundary( STIFF, FORCE, & 523 LoadVector, ExtPressure, SlipCoeff, NormalTangential, & 524 Element, n, nd, ElementNodes ) 525 526 IF (BlockPreconditioning) THEN 527 Alocal = 0.0d0 528 DO p=1,nd 529 DO i=1,dim 530 DO q=1,nd 531 DO j=1,dim 532 Alocal( (p-1)*dim+i, (q-1)*dim+j ) = & 533 stiff( (p-1)*(dim+1)+i, (q-1)*(dim+1)+j) 534 END DO 535 END DO 536 END DO 537 END DO 538 END IF 539 540 CALL DefaultUpdateEquations( STIFF, FORCE ) 541 542 IF (BlockPreconditioning) THEN 543 Force = 0.0d0 544 CALL DefaultUpdateEquations( ALocal, FORCE, USolver=VelocitySolver) 545 END IF 546 547 END IF 548 END DO 549 550 CALL DefaultFinishAssembly() 551 552 CALL DefaultDirichletBCs() 553 554 IF (BlockPreconditioning) THEN 555 ! CALL DefaultFinishAssembly(VelocitySolver) 556 ! CALL DefaultFinishAssembly(PressureSolver) 557 CALL DefaultDirichletBCs(USolver=PressureSolver) 558 CALL DefaultDirichletBCs(USolver=VelocitySolver) 559 END IF 560 561 562 !--------------------------- 563 ! Perform linear solve... 564 !---------------------------- 565 IF (.NOT. BlockPreconditioning) THEN 566 Norm = DefaultSolve() 567 IF ( Solver % Variable % NonlinConverged == 1 ) EXIT 568 ELSE 569 !------------------------------------------------------------------------------- 570 ! This branch performs a locally implemented block preconditioned GCR iteration. 571 ! First, perform linear system scaling by row equilibration and make some 572 ! initializations... 573 !------------------------------------------------------------------------------- 574 IF ( DoScaling ) THEN 575 IF ( DoEquilibration ) THEN 576 CALL RowEquilibration(Solver % Matrix,Solver % Matrix % RHS,Parallel) 577 ELSE 578 CALL ScaleLinearSystem(Solver,Solver % Matrix,Solver % Matrix % RHS) 579 END IF 580 ! store the scaling for the blocks A and Q 581 j = PressureSolver % Matrix % NumberOfRows 582 IF (.NOT. ALLOCATED(RowScalingP)) THEN 583 ALLOCATE(RowScalingP(j)) 584 END IF 585 RowScalingP=1.0_dp 586 IF (.NOT. ALLOCATED(RowScalingA)) THEN 587 ALLOCATE(RowScalingA(AMatrix%NumberOfRows)) 588 END IF 589 RowScalingA=1.0_dp 590 ! for row equilibration we do not scale the A and Q matrices 591 ! so that their symmetry is not destroyed, instead we scale 592 ! the RHS back to the symmetric scaling 593 IF (DoEquilibration) THEN 594 DO i=1,j 595 RowScalingP(i) = 1.0_dp/Solver%Matrix%DiagScaling((dim+1)*i) 596 DO p=1,dim 597 RowScalingA( dim*(i-1)+p ) = & 598 1.0_dp/Solver%Matrix%DiagScaling( (dim+1)*(i-1)+p ) 599 END DO 600 END DO 601 ELSE 602 IF (.NOT. ASSOCIATED(PMatrix%DiagScaling)) THEN 603 ALLOCATE(PMatrix%DiagScaling(j)) 604 END IF 605 PMatrix%DiagScaling(:)=1.0_dp 606 IF (.NOT. ASSOCIATED(AMatrix % DiagScaling)) THEN 607 ALLOCATE(AMatrix%DiagScaling(AMatrix%NumberOfRows)) 608 END IF 609 AMatrix%DiagScaling(:)=1.0_dp 610 DO i=1,j 611 PMatrix%DiagScaling(i) = Solver%Matrix%DiagScaling((dim+1)*i) 612 DO p=1,dim 613 AMatrix%DiagScaling( dim*(i-1)+p ) = & 614 Solver%Matrix%DiagScaling( (dim+1)*(i-1)+p ) 615 END DO 616 END DO 617 CALL ScaleLinearSystem( VelocitySolver, AMatrix, DiagScaling=AMatrix%DiagScaling) 618 CALL ScaleLinearSystem( PressureSolver, PMatrix, DiagScaling=PMatrix%DiagScaling) 619 END IF 620 621 ELSE 622 j = PressureSolver % Matrix % NumberOfRows 623 IF (.NOT. ALLOCATED(RowScalingP)) THEN 624 ALLOCATE(RowScalingP(j)) 625 END IF 626 RowScalingP=1.0_dp 627 IF (.NOT. ALLOCATED(RowScalingA)) THEN 628 ALLOCATE(RowScalingA(AMatrix%NumberOfRows)) 629 END IF 630 RowScalingA=1.0_dp 631 632 IF (.NOT. ASSOCIATED(PMatrix%DiagScaling)) THEN 633 ALLOCATE(PMatrix%DiagScaling(j)) 634 END IF 635 PMatrix%DiagScaling(:)=1.0_dp 636 IF (.NOT. ASSOCIATED(AMatrix % DiagScaling)) THEN 637 ALLOCATE(AMatrix%DiagScaling(AMatrix%NumberOfRows)) 638 END IF 639 AMatrix%DiagScaling(:)=1.0_dp 640 641 IF (.NOT. ASSOCIATED(Solver % Matrix % DiagScaling)) THEN 642 ALLOCATE(Solver % Matrix % DiagScaling(Solver%Matrix%NumberOfRows)) 643 END IF 644 Solver % Matrix % DiagScaling(:)=1.0_dp 645 646 END IF 647 648 !------------------------------------------------------------------------------- 649 ! Switch to using the rotated variables (again) 650 !------------------------------------------------------------------------------- 651 IF (Iter > 1) THEN 652 j = Solver % Matrix % NumberOfRows 653 Solver % Variable % Values(1:j) = RotatedVarValues(1:j) 654 END IF 655 656 IF (Parallel) THEN 657 IF ( .NOT. ASSOCIATED(Solver % Matrix % ParMatrix) ) & 658 CALL ParallelInitMatrix( Solver, Solver % Matrix ) 659 660 n = Solver % Matrix % NumberOfRows 661 IF ( .NOT. ALLOCATED(Residual) ) THEN 662 ALLOCATE( Residual(n) ) 663 Residual = 0.0d0 664 END IF 665 IF ( .NOT. ALLOCATED(TempRes) ) THEN 666 ALLOCATE( TempRes(n) ) 667 TempRes = 0.0d0 668 END IF 669 IF ( .NOT. ALLOCATED(TempS) ) THEN 670 ALLOCATE( TempS(n) ) 671 TempS = 0.0d0 672 END IF 673 IF ( .NOT. ALLOCATED(TempRHS) ) THEN 674 ALLOCATE( TempRHS(n) ) 675 TempRHS = 0.0d0 676 END IF 677 678 679 UpdateMatrix = .TRUE. 680 CALL ParallelInitSolve( Solver % Matrix, Solver % Variable % Values, & 681 Solver % Matrix % RHS, Residual, UpdateMatrix ) 682 683 MMatrix => ParallelMatrix( Solver % Matrix, Mx, Mb, Mr ) 684 n = MMatrix % NumberOfRows 685 686 IF ( .NOT. ALLOCATED(TempMx) ) THEN 687 ALLOCATE( TempMx(n) ) 688 END IF 689 IF ( .NOT. ALLOCATED(TempMb) ) THEN 690 ALLOCATE( TempMb(n) ) 691 END IF 692 IF ( .NOT. ALLOCATED(TempMr) ) THEN 693 ALLOCATE( TempMr(n) ) 694 END IF 695 696 ELSE 697 n = Solver % Matrix % NumberOfRows 698 IF ( .NOT. ALLOCATED(Residual) ) THEN 699 ALLOCATE( Residual(n) ) 700 Residual = 0.0d0 701 END IF 702 IF ( .NOT. ALLOCATED(TempRes) ) THEN 703 ALLOCATE( TempRes(n) ) 704 TempRes = 0.0d0 705 END IF 706 IF ( .NOT. ALLOCATED(TempS) ) THEN 707 ALLOCATE( TempS(n) ) 708 TempS = 0.0d0 709 END IF 710 711 MMatrix => Solver % Matrix 712 Mx => Solver % Variable % Values 713 Mb => Solver % Matrix % RHS 714 Mr => Residual 715 716 END IF 717 718 !------------------------------------------------------------------------------------------ 719 ! Test whether the nonlinear residual is small enough to terminate the nonlinear iteration 720 !------------------------------------------------------------------------------------------ 721 IF (Iter >= 1) THEN 722 !-------------------------------------------------------------------------------------- 723 ! The rotated variables are used to compute the current residual. 724 !-------------------------------------------------------------------------------------- 725 j = Solver % Matrix % NumberOfRows 726 CALL Mymv( Solver % Matrix, Solver % Variable % Values, Residual, .TRUE.) 727 Residual(1:j) = Solver % Matrix % RHS(1:j) - Residual(1:j) 728 729 IF ( Parallel ) THEN 730 DO j=1,SIZE(Mr) 731 Mr(j) = Mb(j) - Mr(j) 732 END DO 733 END IF 734 NonLinError = Mynorm( n, Mr )/Mynorm( n, Mb ) 735 Solver % Variable % NonLinChange = NonLinError 736 737 WRITE(Message,'(a,I4,ES12.3)') 'Residual for nonlinear iterate', & 738 Iter-1, NonLinError 739 CALL Info('StokesSolver', Message, Level=3) 740 741 IF ( NonLinError < NonlinearTol .OR. Iter==NonlinearIter ) THEN 742 DEALLOCATE( Residual ) 743 DEALLOCATE( TempRes ) 744 DEALLOCATE( TempS ) 745 IF (Parallel) THEN 746 DEALLOCATE( TempMx ) 747 DEALLOCATE( TempMr ) 748 DEALLOCATE( TempMb ) 749 DEALLOCATE( TempRHS ) 750 END IF 751 CALL BackRotateNTSystem( Solver % Variable % Values, Solver % Variable % Perm, & 752 Solver % Variable % DOFs ) 753 EXIT 754 END IF 755 !--------------------------------------------------------------------- 756 END IF ! ...Nonlinear convergence check 757 !-------------------------------------------------------------------- 758 759 !------------------------------------------------------------------ 760 ! Initialize the variable containing the new search direction... 761 !------------------------------------------------------------------ 762 Snew = 0.0d0 763 764 IF (AdaptiveTols) THEN 765 !--------------------------------------------------------------------------------------- 766 ! Adapt the linear system convergence tolerance in terms of the current nonlinear error 767 !--------------------------------------------------------------------------------------- 768 IF (Iter >= 1 ) THEN 769 MinTolerance = RelTolerance * NonLinError 770 IF (MinTolerance > BaseTolerance) MinTolerance = BaseTolerance 771 IF (MinTolerance < TargetTol) MinTolerance = TargetTol 772 ELSE 773 MinTolerance = BaseTolerance 774 END IF 775 ELSE 776 MinTolerance = ListGetConstReal( Solver % Values, 'Linear System Convergence Tolerance') 777 END IF 778 779 UseTrueResidual = .TRUE. 780 !---------------------------------------------------------------------- 781 ! Perform solution timing for the actual GCR loop if desired... 782 !--------------------------------------------------------------------- 783 Timing = ListGetLogical(Solver % Values,'Linear System Timing', Found) 784 785 IF( Timing ) THEN 786 t0 = CPUTime(); rt0 = RealTime() 787 END IF 788 789 !--------------------------------------------------------------------------- 790 ! The actual GCR iteration loop... 791 !--------------------------------------------------------------------------- 792 DO Round=1, MaxIterations 793 !----------------------------------------------------------------------- 794 ! Generate a new search direction by applying the preconditioner. That is, 795 ! find an approximation to the error by solving Pe=r, with P the 796 ! preconditioner and r=b-Ax the current residual of the primary system. 797 ! In principle, we can employ the true residual obtained by performing 798 ! the matrix-vector multiplication or utilize the residual vector 799 ! generated by the iterative algorithm. The iterative residual cannot 800 ! be employed currently, since the residual entries corresponding to nodes 801 ! which are not owned by the partition are not received in parallel 802 ! computations... 803 !----------------------------------------------------------------------- 804 IF ( UseTrueResidual ) THEN 805 !------------------------------------------------------------------------- 806 ! Compute the current true residual... 807 !------------------------------------------------------------------------- 808 CALL Mymv( Solver % Matrix, Solver % Variable % Values, Residual, .TRUE. ) 809 j = Solver % Matrix % NumberOfRows 810 Residual(1:j) = Solver % Matrix % RHS(1:j) - Residual(1:j) 811 END IF 812 813 !---------------------------------------------------------------------------- 814 ! Update the preconditioner system rhs by assuming that the preconditioning 815 ! system and the primary system have the same permutation... 816 !-------------------------------------------------------------------------- 817 ! we multiply by the inverse of the row scaling of the global matrix, for 818 ! symmetric scaling we have set the arrays to 1 so nothing happens here. 819 ! For row equilibration this makes sure that the correct system is solved 820 ! without having to scale the A and Q matrices. 821 j = PressureSolver % Matrix % NumberOfRows 822 DO i=1,j 823 PressureSolver % Matrix % RHS(i) = Residual((dim+1)*i) * RowScalingP(i) 824 DO p=1,dim 825 VelocitySolver % Matrix % RHS( dim*(i-1)+p ) = Residual( (dim+1)*(i-1)+p ) & 826 * RowScalingA( dim*(i-1)+p ) 827 END DO 828 END DO 829 830 !-------------------------------------------------------------------------- 831 ! Solve the preconditioning systems sequentially... 832 !-------------------------------------------------------------------------- 833 IF (Round > 1) THEN 834 CALL ListAddLogical(PressureSolver % Values, & 835 'No Precondition Recompute', .TRUE.) 836 CALL ListAddLogical(PressureSolver % Values, & 837 'Linear System Refactorize', .FALSE.) 838 END IF 839 CurrentModel % Solver => PressureSolver 840 CALL Info( 'Pressure Preconditioning', ' ', Level=4 ) 841 CALL Info( 'Pressure Preconditioning', '-------------------------------------', Level=4 ) 842 CALL Info( 'Pressure Preconditioning', 'Performing linear solve', Level=4 ) 843 CALL Info( 'Pressure Preconditioning', '-------------------------------------', Level=4 ) 844 CALL Info( 'Pressure Preconditioning', ' ', Level=4 ) 845 Norm = DefaultSolve(PressureSolver) 846 CurrentModel % Solver => Solver 847 CALL ListAddLogical(PressureSolver % Values, & 848 'No Precondition Recompute', .FALSE.) 849 CALL ListAddLogical(PressureSolver % Values, & 850 'Linear System Refactorize', .TRUE.) 851 852 !----------------------------------------------------------------------------- 853 ! Perform matrix-vector product to produce upper triangular preconditioner... 854 !----------------------------------------------------------------------------- 855 IF (.NOT. Parallel) THEN 856 TempS = 0.0d0 857 j = Solver % Matrix % NumberOfRows/(dim+1) 858 DO i=1,j 859 TempS( i*(dim+1) ) = PressureSolver % Variable % Values(i) 860 END DO 861 CALL Mymv( Solver % Matrix, TempS, TempRes ) 862 DO i=1,j 863 DO p=1,dim 864 VelocitySolver % Matrix % RHS( dim*(i-1)+p ) = & 865 VelocitySolver % Matrix % RHS( dim*(i-1)+p ) - & 866 TempRes( (dim+1)*(i-1)+p ) * RowScalingA( dim*(i-1)+p ) 867 END DO 868 END DO 869 ELSE 870 !-------------------------------------------------------- 871 ! Save the current status of parallel matrix structure 872 !------------------------------------------------------- 873 j = MMatrix % NumberOfRows 874 TempMx(1:j) = Mx(1:j) 875 TempMb(1:j) = Mb(1:j) 876 TempMr(1:j) = Mr(1:j) 877 j=Solver % Matrix % NumberOfRows 878 TempS(1:j) = Solver % Variable % Values(1:j) 879 TempRHS(1:j) = Solver % Matrix % RHS(1:j) 880 881 TempRes = 0.0d0 882 Solver % Matrix % RHS = 0.0d0 883 Solver % Variable % Values = 0.0d0 884 j = Solver % Matrix % NumberOfRows/(dim+1) 885 DO i=1,j 886 Solver % Variable % Values( i*(dim+1) ) = PressureSolver % Variable % Values(i) 887 END DO 888 889 !------------------------------------------------------------------------------- 890 ! Update the vectors of the parallel matrix structure 891 !------------------------------------------------------------------------------- 892 CALL ParallelUpdateSolve( Solver % Matrix, Solver % Variable % Values, TempRes ) 893 !----------------------------------------------------------------------- 894 ! Perform mv-product and insert the result into the argument vectors... 895 !----------------------------------------------------------------------- 896 CALL Mymv( Solver % Matrix, Solver % Variable % Values, TempRes, .TRUE.) 897 898 !------------------------------------------------------------------------- 899 ! Update the velocity solver rhs by adding the result of the gradient mv... 900 !------------------------------------------------------------------------- 901 DO i=1,j 902 DO p=1,dim 903 VelocitySolver % Matrix % RHS( dim*(i-1)+p ) = & 904 VelocitySolver % Matrix % RHS( dim*(i-1)+p ) - & 905 TempRes( (dim+1)*(i-1)+p ) * RowScalingA(dim*(i-1)+p ) 906 907 END DO 908 END DO 909 !----------------------------------------------------------------------------- 910 ! Finally retrieve the original status of the parallel matrix structure.... 911 !---------------------------------------------------------------------------- 912 j=Solver % Matrix % NumberOfRows 913 Solver % Variable % Values(1:j) = TempS(1:j) 914 Solver % Matrix % RHS(1:j) = TempRHS(1:j) 915 916 j = MMatrix % NumberOfRows 917 Mx(1:j) = TempMx(1:j) 918 Mb(1:j) = TempMb(1:j) 919 Mr(1:j) = TempMr(1:j) 920 END IF 921 922 IF (Round > 1) THEN 923 CALL ListAddLogical(VelocitySolver % Values, & 924 'No Precondition Recompute', .TRUE.) 925 CALL ListAddLogical(VelocitySolver % Values, & 926 'Linear System Refactorize', .FALSE.) 927 END IF 928 CurrentModel % Solver => VelocitySolver 929 !---------------------------------------------------------------------------- 930 ! Adapt the preconditioning system convergence tolerance in terms of the 931 ! current convergence tolerance for the primary system 932 !---------------------------------------------------------------------------- 933 !IF (AdaptiveTols) THEN 934 IF (.FALSE.) THEN ! Disabled: Using a fixed tolerance may actually be a better strategy 935 PrecondTol = 1/RelTolerance * MinTolerance 936 IF (PrecondTol > BaseTolerance) PrecondTol = BaseTolerance 937 CALL ListAddConstReal(VelocitySolver % Values, & 938 'Linear System Convergence Tolerance', BaseTolerance) 939 END IF 940 CALL Info( 'Velocity Preconditioning', ' ', Level=4 ) 941 CALL Info( 'Velocity Preconditioning', '-------------------------------------', Level=4 ) 942 CALL Info( 'Velocity Preconditioning', 'Performing linear solve', Level=4 ) 943 CALL Info( 'Velocity Preconditioning', '-------------------------------------', Level=4 ) 944 CALL Info( 'Velocity Preconditioning', ' ', Level=4 ) 945 Norm = DefaultSolve(VelocitySolver) 946 CurrentModel % Solver => Solver 947 CALL ListAddLogical(VelocitySolver % Values, & 948 'No Precondition Recompute', .FALSE.) 949 CALL ListAddLogical(VelocitySolver % Values, & 950 'Linear System Refactorize', .TRUE.) 951 952 !------------------------------------------------------------------------- 953 ! Insert the computed approximation of the error into the search direction 954 ! variable which will be given to the GCR routine... 955 !-------------------------------------------------------------------------- 956 DO t=1, Active 957 Element => GetActiveElement(t) 958 nd = GetElementDOFs( Indexes ) 959 960 DO i=1,nd 961 j = Solver % Variable % Perm( Indexes(i) ) 962 k = PressureSolver % Variable % Perm( Indexes(i) ) 963 Snew( j*(dim+1) ) = PressureSolver % Variable % Values(k) 964 965 k = VelocitySolver % Variable % Perm( Indexes(i) ) 966 DO p=1,dim 967 Snew( (j-1)*(dim+1)+p ) = VelocitySolver % Variable % Values( (j-1)*dim+p ) 968 END DO 969 970 END DO 971 END DO 972 973 !----------------------------------------------------------------------------- 974 ! Perform GCR update for solving the primitive linear system... 975 !----------------------------------------------------------------------------- 976 CALL GCRUpdate(n, Solver % Matrix, MMatrix, Mx, Mb, Mr, Snew, S, V, R, Round, & 977 Norm, RestartM) 978 979 IF (Parallel) & 980 CALL ParallelUpdateResult( Solver % Matrix, Solver % Variable % Values, Residual ) 981 982 IF (Norm < MinTolerance) EXIT 983 984 END DO 985 986 ! unscale solution if column scaling was used 987 IF (DoScaling .AND. (.NOT. DoEquilibration)) THEN 988 CALL BackScaleLinearSystem(Solver,Solver%Matrix,Solver%Matrix%RHS,Solver%Variable%Values) 989 END IF 990 991 IF( Timing ) THEN 992 st = CPUTime() - t0; 993 rst = RealTime() - rt0 994 995 CALL ListAddConstReal(CurrentModel % Simulation,'res: linsys cpu time '& 996 //GetVarName(Solver % Variable),st) 997 CALL ListAddConstReal(CurrentModel % Simulation,'res: linsys real time '& 998 //GetVarName(Solver % Variable),rst) 999 WRITE(Message,'(a,f8.2,f8.2,a)') 'Linear system time (CPU,REAL) for '& 1000 //GetVarName(Solver % Variable)//': ',st,rst,' (s)' 1001 CALL Info('SolveSystem',Message) 1002 1003 IF( ListGetLogical(Solver % Values,'Linear System Timing Cumulative',Found)) THEN 1004 ct = ListGetConstReal(CurrentModel % Simulation,'res: cum linsys cpu time '& 1005 //GetVarName(Solver % Variable),Found) 1006 st = st + ct 1007 ct = ListGetConstReal(CurrentModel % Simulation,'res: cum linsys real time '& 1008 //GetVarName(Solver % Variable),Found) 1009 rst = rst + ct 1010 CALL ListAddConstReal(CurrentModel % Simulation,'res: cum linsys cpu time '& 1011 //GetVarName(Solver % Variable),st) 1012 CALL ListAddConstReal(CurrentModel % Simulation,'res: cum linsys real time '& 1013 //GetVarName(Solver % Variable),rst) 1014 END IF 1015 1016 END IF 1017 1018 DEALLOCATE( Residual ) 1019 DEALLOCATE( TempRes ) 1020 DEALLOCATE( TempS ) 1021 1022 ! Compute Var Loads if needed 1023 IF ( ListGetLogical( Solver % Values,'Calculate Loads', Found ) ) & 1024 CALL ComputeVarLoads(Solver) 1025 1026 1027 n = Solver % Matrix % NumberOfRows 1028 RotatedVarValues(1:n) = Solver % Variable % Values(1:n) 1029 CALL BackRotateNTSystem( Solver % Variable % Values, Solver % Variable % Perm, & 1030 Solver % Variable % DOFs ) 1031 1032 1033 END IF 1034 1035 END DO 1036 1037 IF (BlockPreconditioning) & 1038 DEALLOCATE(RotatedVarValues) 1039 1040 IF (P2P1) THEN 1041 !---------------------------------------------------------------------------------------- 1042 ! Replace the zero pressure solution at the nodes which are not needed in the linear 1043 ! pressure approximation by the interpolated values for right visualization: 1044 !---------------------------------------------------------------------------------------- 1045 DO t=1,Active 1046 ! First the midedge nodes: 1047 Element => GetActiveElement(t) 1048 nd = GetElementDOFs( Indexes ) 1049 k = GetElementFamily() 1050 EdgeMap => GetEdgeMap(k) 1051 SELECT CASE(k) 1052 CASE (3) 1053 MidEdgeNodes(1:3) = (/ 4, 5, 6 /) 1054 CASE (4) 1055 MidEdgeNodes(1:4) = (/ 5, 6, 7, 8 /) 1056 CASE (5) 1057 MidEdgeNodes(1:6) = (/ 5, 6, 7, 8, 9, 10 /) 1058 CASE (6) 1059 MidEdgeNodes(1:8) = (/ 6, 7, 8, 9, 10, 11, 12, 13 /) 1060 CASE (7) 1061 MidEdgeNodes(1:9) = (/ 7, 8, 9, 10, 11, 12, 13, 14, 15 /) 1062 CASE (8) 1063 MidEdgeNodes(1:12) = (/ 9, 10, 11, 12, 17, 18, 19, 20, 13, 14, 15, 16 /) 1064 END SELECT 1065 1066 DO q=1,SIZE(EdgeMap,1) 1067 m = (dim+1) * Solver % Variable % Perm(Indexes(MidEdgeNodes(q))) 1068 i = (dim+1) * Solver % Variable % Perm(Indexes(EdgeMap(q,1))) 1069 j = (dim+1) * Solver % Variable % Perm(Indexes(EdgeMap(q,2))) 1070 Solver % Variable % Values(m) = 0.5d0 * ( Solver % Variable % Values(i) + & 1071 Solver % Variable % Values(j) ) 1072 END DO 1073 1074 ! The pressure at the midface nodes for 409 elements: 1075 IF (k==4 .AND. nd==9) THEN 1076 res = 0.0d0 1077 DO q=1,4 1078 i = (dim+1) * Solver % Variable % Perm(Indexes(q)) 1079 res = res + Solver % Variable % Values(i) 1080 END DO 1081 m = (dim+1) * Solver % Variable % Perm(Indexes(9)) 1082 Solver % Variable % Values(m) = 0.25d0 * res 1083 END IF 1084 1085 ! The pressure at the midpoint and at the midface nodes for 827 elements: 1086 IF (k==8 .AND. nd==27) THEN 1087 BrickFaceMap(1,:) = (/ 1,2,6,5 /) 1088 BrickFaceMap(2,:) = (/ 2,3,7,6 /) 1089 BrickFaceMap(3,:) = (/ 4,3,7,8 /) 1090 BrickFaceMap(4,:) = (/ 1,4,8,5 /) 1091 BrickFaceMap(5,:) = (/ 1,2,3,4 /) 1092 BrickFaceMap(6,:) = (/ 5,6,7,8 /) 1093 DO j=1,6 1094 res = 0.0d0 1095 DO q=1,4 1096 i = (dim+1) * Solver % Variable % Perm(Indexes(BrickFaceMap(j,q))) 1097 res = res + Solver % Variable % Values(i) 1098 END DO 1099 m = (dim+1) * Solver % Variable % Perm(Indexes(20+j)) 1100 Solver % Variable % Values(m) = 0.25d0 * res 1101 END DO 1102 1103 res = 0.0d0 1104 DO q=1,8 1105 i = (dim+1) * Solver % Variable % Perm(Indexes(q)) 1106 res = res + Solver % Variable % Values(i) 1107 END DO 1108 m = (dim+1) * Solver % Variable % Perm(Indexes(27)) 1109 Solver % Variable % Values(m) = 0.125d0 * res 1110 END IF 1111 1112 END DO 1113 END IF 1114 1115CONTAINS 1116 1117 1118!------------------------------------------------------------------------------ 1119 SUBROUTINE GCRUpdate( n, A, M, x, b, r, Snew, S, V, RR, Round, res, MParam) 1120!------------------------------------------------------------------------------ 1121 INTEGER :: n, Round, MParam 1122 TYPE(Matrix_t), POINTER :: A, M 1123 REAL(KIND=dp) :: x(:), b(:), r(:), Snew(:) 1124 REAL(KIND=dp) :: S(:,:), V(:,:), RR(:), res 1125!-------------------------------------------------------------------------------- 1126 REAL(KIND=dp) :: T1(n), T2(n), beta, alpha 1127 INTEGER :: i,j,k 1128!-------------------------------------------------------------------------------- 1129 1130 IF ( Parallel ) CALL ParallelVector(A,Snew) 1131 1132 !---------------------------------------------- 1133 ! Check for restarting 1134 !--------------------------------------------- 1135 IF ( MOD(Round,MParam)==0 ) THEN 1136 j = MParam 1137 ELSE 1138 j = MOD(Round,MParam) 1139 END IF 1140 1141 IF ( j == 1) THEN 1142 CALL Mymv( A, x, r ) 1143 r(1:n) = b(1:n)-r(1:n) 1144 res = MyNorm(n,r)/Mynorm(n,b) 1145 IF (Round==1) THEN 1146 WRITE(Message,'(a,I4,ES12.3)') 'GCR residual for iterate', & 1147 0, res 1148 CALL Info('Outer Iteration',Message,Level=4) 1149 END IF 1150 ELSE 1151 r(1:n) = RR(1:n) 1152 END IF 1153 1154 T1(1:n) = Snew(1:n) 1155 CALL Mymv( A, T1, T2 ) 1156 !-------------------------------------------------------------- 1157 ! Perform the orthogonalisation of the search directions... 1158 !-------------------------------------------------------------- 1159 DO i=1,j-1 1160 beta = Mydot( n, V(1:n,i), T2(1:n) ) 1161 T1(1:n) = T1(1:n) - beta * S(1:n,i) 1162 T2(1:n) = T2(1:n) - beta * V(1:n,i) 1163 END DO 1164 1165 alpha = Mynorm(n,T2) 1166 T1(1:n) = 1.0d0/alpha * T1(1:n) 1167 T2(1:n) = 1.0d0/alpha * T2(1:n) 1168 1169 !------------------------------------------------------------- 1170 ! The update of the solution and save the search data... 1171 !------------------------------------------------------------- 1172 beta = Mydot(n, T2, r) 1173 1174 x(1:n) = x(1:n) + beta * T1(1:n) 1175 r(1:n) = r(1:n) - beta * T2(1:n) 1176 1177 IF ( j /= MParam ) THEN 1178 S(1:n,j) = T1(1:n) 1179 V(1:n,j) = T2(1:n) 1180 END IF 1181 1182 RR(1:n) = r(1:n) 1183 res = MyNorm(n,r)/Mynorm(n,b) 1184 1185 WRITE(Message,'(a,I4,ES12.3)') 'GCR residual for iterate', & 1186 Round, res 1187 CALL Info('Outer Iteration',Message,Level=4) 1188 1189!------------------------------------------- 1190 END SUBROUTINE GCRUpdate 1191!--------------------------------------------- 1192 1193 1194 1195 1196!------------------------------------------------------------------------------ 1197 SUBROUTINE LocalMatrix( STIFF, VeloBlock, Mass, FORCE, LOAD, Nodalrho, & 1198 Nodalmu, Vx, Vy, Vz, Element, n, nd, dim, Convect, P2P1, SkipPowerLaw, & 1199 Newton, DiagonalA) 1200!------------------------------------------------------------------------------ 1201 REAL(KIND=dp), TARGET :: STIFF(:,:), VeloBlock(:,:), Mass(:,:), FORCE(:), LOAD(:,:) 1202 REAL(KIND=dp) :: Nodalmu(:), Nodalrho(:), Vx(:), Vy(:), Vz(:) 1203 INTEGER :: dim, n, nd 1204 TYPE(Element_t), POINTER :: Element 1205 LOGICAL :: Convect, P2P1, SkipPowerLaw, Newton, DiagonalA 1206 !------------------------------------------------------------------------------ 1207 TYPE(Element_t), POINTER :: PressureElement => NULL() 1208 REAL(KIND=dp) :: Basis(nd), dBasisdx(nd,3), & 1209 DetJ, LoadAtIP(dim+1), Velo(dim), Grad(dim,dim), w1, w2, w3, ViscAtIp, RhoAtIP 1210 REAL(KIND=dp) :: LinBasis(nd) 1211 REAL(KIND=dp), POINTER :: A(:,:), F(:), Jac(:,:) 1212 REAL(KIND=dp), TARGET :: JacM(nd*(dim+1),nd*(dim+1)), Sol(nd*(dim+1)) 1213 LOGICAL :: Stat, ViscNewtonLin 1214 INTEGER :: t, i, j, k, l, p, q, nlin 1215 INTEGER :: LinearCode(3:8) = (/ 303,404,504,605,706,808 /) 1216 TYPE(GaussIntegrationPoints_t) :: IP 1217 1218 REAL(KIND=dp) :: mu = 1.0d0, s, & 1219 a1, a2, muder, muder0, Strain(dim,dim) 1220 1221 TYPE(Nodes_t) :: Nodes 1222 SAVE Nodes, PressureElement 1223 !------------------------------------------------------------------------------ 1224 1225 IF (P2P1) THEN 1226 IF ( .NOT. ASSOCIATED(PressureElement) ) PressureElement => AllocateElement() 1227 k = GetElementFamily() 1228 PressureElement % Type => GetElementType( LinearCode(k), .FALSE. ) 1229 nlin = PressureElement % Type % NumberOfNodes 1230 END IF 1231 1232 CALL GetElementNodes( Nodes ) 1233 STIFF = 0.0d0 1234 FORCE = 0.0d0 1235 Mass = 0.0d0 1236 VeloBlock = 0.0d0 1237 IF (Newton) JacM = 0.0d0 1238 ViscNewtonLin = .FALSE. 1239 1240 !---------------------- 1241 ! Numerical integration: 1242 !----------------------- 1243 IP = GaussPoints( Element ) 1244 1245 DO t=1,IP % n 1246 !-------------------------------------------------------------- 1247 ! Basis function values & derivatives at the integration point: 1248 !-------------------------------------------------------------- 1249 stat = ElementInfo( Element, Nodes, IP % U(t), IP % V(t), & 1250 IP % W(t), detJ, Basis, dBasisdx ) 1251 1252 !-------------------------------------------------------------- 1253 ! Linear pressure basis functions for P2-P1 approximation: 1254 !-------------------------------------------------------------- 1255 IF (P2P1) CALL NodalBasisFunctions( nlin, LinBasis, PressureElement, & 1256 IP % U(t), IP % V(t), IP % W(t) ) 1257 1258 s = IP % s(t) * detJ 1259 1260 Grad = 0.0d0 1261 DO i=1,dim 1262 Grad(1,i) = SUM( Vx(1:n) * dBasisdx(1:n,i) ) 1263 Grad(2,i) = SUM( Vy(1:n) * dBasisdx(1:n,i) ) 1264 IF ( DIM > 2 ) Grad(3,i) = SUM( Vz(1:n) * dBasisdx(1:n,i) ) 1265 END DO 1266 1267 !---------------------------------------------- 1268 ! Material parameters at the integration point: 1269 !---------------------------------------------- 1270 ViscAtIP = SUM( Basis(1:n) * Nodalmu(1:n) ) 1271 RhoAtIP = SUM( Basis(1:n) * Nodalrho(1:n) ) 1272 1273 IF ( SkipPowerLaw ) THEN 1274 mu = ViscAtIP 1275 ELSE 1276 IF( ListCheckPresent( Material, 'Viscosity Model' ) ) THEN 1277 mu = EffectiveViscosity( ViscAtIP, RhoAtIp, Vx, Vy, Vz, & 1278 Element, Nodes, n, n, IP % U(t), IP % V(t), & 1279 IP % W(t), muder0, LocalIP=t ) 1280 ViscNewtonLin = Newton .AND. muder0/= 0.0d0 1281 IF ( ViscNewtonLin ) Strain = (Grad+TRANSPOSE(Grad))/2 1282 ELSE 1283 mu = ViscAtIP 1284 END IF 1285 END IF 1286 1287 IF (Convect) THEN 1288 w1 = SUM( Vx(1:n) * Basis(1:n) ) 1289 w2 = SUM( Vy(1:n) * Basis(1:n) ) 1290 IF (dim > 2) w3 = SUM( Vz(1:n) * Basis(1:n) ) 1291 END IF 1292 1293 !-------------------------------------------- 1294 ! The source term at the integration point: 1295 !------------------------------------------ 1296 LoadAtIP(1:dim+1) = MATMUL( Basis(1:n), LOAD(1:n,1:dim+1) ) 1297 1298 !---------------------------------------------------------------------------------- 1299 ! The system matrix with linear pressure approximation 1300 !--------------------------------------------------------------------------------- 1301 DO p=1,nd 1302 DO q=1,nd 1303 i = (dim+1) * (p-1) + 1 1304 j = (dim+1) * (q-1) + 1 1305 A => STIFF(i:i+dim,j:j+dim) 1306 IF ( ViscNewtonLin ) Jac => JacM( i:i+dim,j:j+dim ) 1307 1308 IF ( ViscNewtonLin ) THEN 1309 DO i=1,dim 1310 muder = 4.0d0*muder0*SUM(Strain(i,1:dim)*dBasisdx(q,1:dim)) 1311 DO j=1,dim 1312 a1 = 0.0d0 1313 a2 = 0.0d0 1314 DO k=1,dim 1315 a1 = a1 + Grad(j,k)*dBasisdx(p,k) 1316 a2 = a2 + Grad(k,j)*dBasisdx(p,k) 1317 END DO 1318 Jac(j,i) = Jac(j,i) + s * muder * (a1+a2) 1319 ! Jac(i,i)=Jac(i,i) + s*muder*Grad(i,j)*dBasisdx(p,j) 1320 ! Jac(j,i)=Jac(j,i) + s*muder*Grad(i,j)*dBasisdx(p,i) 1321 END DO 1322 END DO 1323 END IF 1324 1325 DO i=1,dim 1326 DO j = 1,dim 1327 A(i,i) = A(i,i) + s * mu * dBasisdx(q,j) * dBasisdx(p,j) 1328 A(i,j) = A(i,j) + s * mu * dBasisdx(q,i) * dBasisdx(p,j) 1329 END DO 1330 1331 ! Only linear pressure approximation is used... 1332 IF (P2P1) THEN 1333 IF (q <= nlin) & 1334 A(i,dim+1) = A(i,dim+1) - s * LinBasis(q) * dBasisdx(p,i) 1335 IF (p <= nlin) & 1336 A(dim+1,i) = A(dim+1,i) - s * dBasisdx(q,i) * LinBasis(p) 1337 ELSE 1338 IF (q <= n) & 1339 A(i,dim+1) = A(i,dim+1) - s * Basis(q) * dBasisdx(p,i) 1340 1341 IF (p <= n) & 1342 A(dim+1,i) = A(dim+1,i) - s * dBasisdx(q,i) * Basis(p) 1343 END IF 1344 END DO 1345 1346 IF (DiagonalA) THEN 1347 DO i=1,dim 1348 DO j = 1,dim 1349 1350 VeloBlock( dim*(p-1)+i, dim*(q-1)+i ) = VeloBlock( dim*(p-1)+i, dim*(q-1)+i ) + & 1351 s * mu * dBasisdx(q,j) * dBasisdx(p,j) 1352 1353 END DO 1354 END DO 1355 END IF 1356 1357 IF (Convect) THEN 1358 A(1,1) = A(1,1) + s * RhoAtIP * w1 * dBasisdx(q,1) * Basis(p) 1359 A(1,1) = A(1,1) + s * RhoAtIP * w2 * dBasisdx(q,2) * Basis(p) 1360 A(2,2) = A(2,2) + s * RhoAtIP * w1 * dBasisdx(q,1) * Basis(p) 1361 A(2,2) = A(2,2) + s * RhoAtIP * w2 * dBasisdx(q,2) * Basis(p) 1362 1363 IF (dim > 2) THEN 1364 A(1,1) = A(1,1) + s * RhoAtIP * w3 * dBasisdx(q,3) * Basis(p) 1365 A(2,2) = A(2,2) + s * RhoAtIP * w3 * dBasisdx(q,3) * Basis(p) 1366 A(3,3) = A(3,3) + s * RhoAtIP * w1 * dBasisdx(q,1) * Basis(p) 1367 A(3,3) = A(3,3) + s * RhoAtIP * w2 * dBasisdx(q,2) * Basis(p) 1368 A(3,3) = A(3,3) + s * RhoAtIP * w3 * dBasisdx(q,3) * Basis(p) 1369 END IF 1370 1371 END IF 1372 1373 1374 END DO 1375 1376 i = (dim+1) * (p-1) + 1 1377 F => FORCE(i:i+dim) 1378 F = F + s * RhoAtIP * LoadAtIP * Basis(p) 1379 1380 END DO 1381 1382 IF (P2P1) THEN 1383 DO p=1,nlin 1384 DO q=1,nlin 1385 Mass(p,q) = Mass(p,q) - s * 1.0d0/mu * Basis(p) * Basis(q) 1386 END DO 1387 END DO 1388 ELSE 1389 DO p=1,n 1390 DO q=1,n 1391 Mass(p,q) = Mass(p,q) - s * 1.0d0/mu * Basis(p) * Basis(q) 1392 END DO 1393 END DO 1394 END IF 1395 END DO 1396 1397 IF ( ViscNewtonLin ) THEN 1398 SOL=0.0d0 1399 SOL(1:(dim+1)*n:dim+1) = Vx(1:n) 1400 SOL(2:(dim+1)*n:dim+1) = Vy(1:n) 1401 IF ( dim>2 ) SOL(3:(dim+1)*n:dim+1) = Vz(1:n) 1402 p = (dim+1)*n 1403 Stiff(1:p,1:p) = Stiff(1:p,1:p)+JacM(1:p,1:p) 1404 FORCE(1:p)=FORCE(1:p)+MATMUL(JacM(1:p,1:p),SOL(1:p)) 1405 END IF 1406 1407 ! Omit dofs that do not correspond to linear pressure approximation... 1408 IF (P2P1) THEN 1409 DO p = nlin+1,nd 1410 i = (dim+1) * p 1411 FORCE(i) = 0.0d0 1412 STIFF(i,:) = 0.0d0 1413 STIFF(:,i) = 0.0d0 1414 STIFF(i,i) = 1.0d0 1415 Mass(p,p) = 1.0d0 1416 END DO 1417 ELSE 1418 DO p = n+1,nd 1419 i = (dim+1) * p 1420 FORCE(i) = 0.0d0 1421 STIFF(i,:) = 0.0d0 1422 STIFF(:,i) = 0.0d0 1423 STIFF(i,i) = 1.0d0 1424 Mass(p,p) = 1.0d0 1425 END DO 1426 END IF 1427 1428!------------------------------------------------------------------------------ 1429 END SUBROUTINE LocalMatrix 1430!------------------------------------------------------------------------------ 1431 1432 1433!------------------------------------------------------------------------------ 1434 SUBROUTINE StaticCondensation( N, Nb, dim, K, F ) 1435!------------------------------------------------------------------------------ 1436 USE LinearAlgebra 1437 INTEGER :: N, Nb, dim 1438 REAL(KIND=dp) :: K(:,:), F(:), Kbb(nb*dim,nb*dim), & 1439 Kbl(nb*dim,n*(dim+1)),Klb(n*(dim+1),nb*dim),Fb(nb*dim) 1440 1441 INTEGER :: m, i, j, l, p, Cdofs((dim+1)*n), Bdofs(dim*nb) 1442 1443 m = 0 1444 DO p = 1,n 1445 DO i = 1,dim+1 1446 m = m + 1 1447 Cdofs(m) = (dim+1)*(p-1) + i 1448 END DO 1449 END DO 1450 1451 m = 0 1452 DO p = 1,nb 1453 DO i = 1,dim 1454 m = m + 1 1455 Bdofs(m) = (dim+1)*(p-1) + i + n*(dim+1) 1456 END DO 1457 END DO 1458 1459 Kbb = K(Bdofs,Bdofs) 1460 Kbl = K(Bdofs,Cdofs) 1461 Klb = K(Cdofs,Bdofs) 1462 Fb = F(Bdofs) 1463 1464 CALL InvertMatrix( Kbb,nb*dim ) 1465 1466 F(1:(dim+1)*n) = F(1:(dim+1)*n) - MATMUL( Klb, MATMUL( Kbb, Fb ) ) 1467 K(1:(dim+1)*n,1:(dim+1)*n) = & 1468 K(1:(dim+1)*n,1:(dim+1)*n) - MATMUL( Klb, MATMUL( Kbb,Kbl ) ) 1469 1470!------------------------------------------------------------------------------ 1471 END SUBROUTINE StaticCondensation 1472!------------------------------------------------------------------------------ 1473 1474 1475 1476 1477!------------------------------------------------------------------------------ 1478 FUNCTION Mydot( n, x, y ) RESULT(s) 1479!------------------------------------------------------------------------------ 1480 INTEGER :: n 1481 REAL(KIND=dp) :: s,x(:),y(:) 1482!------------------------------------------------------------------------------ 1483 IF ( .NOT. Parallel ) THEN 1484 s = DOT_PRODUCT( x(1:n), y(1:n) ) 1485 ELSE 1486 s = ParallelDot( n, x, y ) 1487 END IF 1488!------------------------------------------------------------------------------ 1489 END FUNCTION Mydot 1490!------------------------------------------------------------------------------ 1491 1492 1493 1494 1495!------------------------------------------------------------------------------ 1496 SUBROUTINE Mymv( A, x, b, Update ) 1497!------------------------------------------------------------------------------ 1498 REAL(KIND=dp) :: x(:), b(:) 1499 LOGICAL, OPTIONAL :: Update 1500 TYPE(Matrix_t), POINTER :: A 1501!------------------------------------------------------------------------------ 1502 IF ( .NOT. Parallel ) THEN 1503 CALL CRS_MatrixVectorMultiply( A, x, b ) 1504 ELSE 1505 IF ( PRESENT( Update ) ) THEN 1506 CALL ParallelMatrixVector( A,x,b,Update ) 1507 ELSE 1508 CALL ParallelMatrixVector( A,x,b ) 1509 END IF 1510 END IF 1511!------------------------------------------------------------------------------ 1512 END SUBROUTINE Mymv 1513!------------------------------------------------------------------------------ 1514 1515 1516!------------------------------------------------------------------------------ 1517 FUNCTION Mynorm( n, x ) RESULT(s) 1518!------------------------------------------------------------------------------ 1519 INTEGER :: n 1520 REAL(KIND=dp) :: s,x(:) 1521!------------------------------------------------------------------------------ 1522 IF ( .NOT. Parallel ) THEN 1523 s = SQRT( DOT_PRODUCT( x(1:n), x(1:n) ) ) 1524 ELSE 1525 s = ParallelNorm( n, x ) 1526 END IF 1527!------------------------------------------------------------------------------ 1528 END FUNCTION Mynorm 1529!------------------------------------------------------------------------------ 1530 1531 1532!------------------------------------------------------------------------------ 1533 SUBROUTINE SetBoundaryConditions( Model, StiffMatrix, & 1534 Name, DOF, NDOFs, Perm, rhs ) 1535!------------------------------------------------------------------------------ 1536! 1537! Set dirichlet boundary condition for given dof 1538! 1539! TYPE(Model_t) :: Model 1540! INPUT: the current model structure 1541! 1542! TYPE(Matrix_t), POINTER :: StiffMatrix 1543! INOUT: The global matrix 1544! 1545! CHARACTER(LEN=*) :: Name 1546! INPUT: name of the dof to be set 1547! 1548! INTEGER :: DOF, NDOFs 1549! INPUT: The order number of the dof and the total number of DOFs for 1550! this equation 1551! 1552! INTEGER :: Perm(:) 1553! INPUT: The node reordering info, this has been generated at the 1554! beginning of the simulation for bandwidth optimization 1555!****************************************************************************** 1556!------------------------------------------------------------------------------ 1557 TYPE(Model_t) :: Model 1558 TYPE(Matrix_t), POINTER :: StiffMatrix 1559 1560 CHARACTER(LEN=*) :: Name 1561 INTEGER :: DOF, NDOFs, Perm(:) 1562 REAL(KIND=dp), OPTIONAL :: rhs(:) 1563!------------------------------------------------------------------------------ 1564 1565 TYPE(Element_t), POINTER :: CurrentElement 1566 INTEGER, POINTER :: NodeIndexes(:) 1567 INTEGER :: i,j,k,n,t,nd 1568 LOGICAL :: GotIt 1569 REAL(KIND=dp) :: Work(Model % MaxElementNodes) 1570 1571!------------------------------------------------------------------------------ 1572 1573 DO t = Model % NumberOfBulkElements + 1, & 1574 Model % NumberOfBulkElements + Model % NumberOfBoundaryElements 1575 1576 CurrentElement => Model % Elements(t) 1577!------------------------------------------------------------------------------ 1578! Set the current element pointer in the model structure to 1579! reflect the element being processed 1580!------------------------------------------------------------------------------ 1581 Model % CurrentElement => Model % Elements(t) 1582!------------------------------------------------------------------------------ 1583 !n = CurrentElement % TYPE % NumberOfNodes 1584 1585 n = GetElementNOFNodes() 1586 nd = GetElementDOFs( Indexes ) 1587 1588 !------------------------------------------------------------------------------ 1589 ! If standard element types are used, there is no need for these modifications 1590 !------------------------------------------------------------------------------ 1591 IF (n==nd) CYCLE 1592 1593 NodeIndexes => CurrentElement % NodeIndexes(1:n) 1594 1595 DO i=1,Model % NumberOfBCs 1596 IF ( CurrentElement % BoundaryInfo % Constraint == & 1597 Model % BCs(i) % Tag ) THEN 1598 1599 Work(1:n) = ListGetReal( Model % BCs(i) % Values, & 1600 Name,n,NodeIndexes, gotIt ) 1601 IF ( gotIt ) THEN 1602 DO j=1,n 1603 k = Perm(NodeIndexes(j)) 1604 IF ( k > 0 ) THEN 1605 k = NDOFs * (k-1) + DOF 1606 CALL ZeroRow( StiffMatrix,k ) 1607 CALL SetMatrixElement( StiffMatrix,k,k, 1.0d0 ) 1608 IF ( PRESENT(rhs) ) rhs(k) = work(j) 1609 END IF 1610 END DO 1611 1612 DO j=n+1,nd 1613 k = Perm(Indexes(j)) 1614 k = NDOFs * (k-1) + DOF 1615 CALL ZeroRow( StiffMatrix,k ) 1616 CALL SetMatrixElement( StiffMatrix,k,k, 1.0d0 ) 1617 IF ( PRESENT(rhs) ) rhs(k) = 0.0d0 1618 END DO 1619 1620 END IF 1621 END IF 1622 END DO 1623 END DO 1624!------------------------------------------------------------------------------ 1625 END SUBROUTINE SetBoundaryConditions 1626!------------------------------------------------------------------------------ 1627 1628 1629 1630 1631 1632 1633 1634!------------------------------------------------------------------------------ 1635 SUBROUTINE NavierStokesBoundary( BoundaryMatrix, BoundaryVector, LoadVector, & 1636 NodalExtPressure, NodalSlipCoeff, NormalTangential, Element, n, nd, Nodes ) 1637 1638!------------------------------------------------------------------------------ 1639!****************************************************************************** 1640! 1641! Return element local matrices and RSH vector for Navier-Stokes-equations 1642! boundary conditions. 1643! 1644! ARGUMENTS: 1645! 1646! REAL(KIND=dp) :: BoundaryMatrix(:,:) 1647! OUTPUT: time derivative coefficient matrix 1648! 1649! REAL(KIND=dp) :: BoundaryVector(:) 1650! OUTPUT: RHS vector 1651! 1652! REAL(KIND=dp) :: LoadVector(:,:) 1653! INPUT: Nodal values force in coordinate directions 1654! 1655! REAL(KIND=dp) :: NodalAlpha(:,:) 1656! INPUT: Nodal values of force in normal direction 1657! 1658! REAL(KIND=dp) :: NodalBeta(:,:) 1659! INPUT: Nodal values of something which will be taken derivative in 1660! tangential direction and added to force... 1661! 1662! 1663! TYPE(Element_t) :: Element 1664! INPUT: Structure describing the element (dimension,nof nodes, 1665! interpolation degree, etc...) 1666! 1667! INTEGER :: n 1668! INPUT: Number of boundary element nodes 1669! 1670! TYPE(Nodes_t) :: Nodes 1671! INPUT: Element node coordinates 1672! 1673!****************************************************************************** 1674!------------------------------------------------------------------------------ 1675 USE ElementUtils 1676 1677 IMPLICIT NONE 1678 1679 REAL(KIND=dp) :: BoundaryMatrix(:,:), BoundaryVector(:), LoadVector(:,:), & 1680 NodalSlipCoeff(:,:), NodalExtPressure(:) 1681 1682 INTEGER :: n, nd, pn 1683 TYPE(Element_t),POINTER :: Element 1684 TYPE(Nodes_t) :: Nodes 1685 LOGICAL :: NormalTangential 1686 1687!------------------------------------------------------------------------------ 1688! Local variables 1689!------------------------------------------------------------------------------ 1690 REAL(KIND=dp) :: Basis(n),dBasisdx(n,3) 1691 REAL(KIND=dp) :: detJ,SlipCoeff 1692 REAL(KIND=dp) :: u,v,w,s 1693 REAL(KIND=dp), POINTER :: U_Integ(:),V_Integ(:),W_Integ(:),S_Integ(:) 1694 REAL(KIND=dp) :: TangentForce(3),Force(3),Normal(3),Tangent(3),Tangent2(3), & 1695 Vect(3), Alpha 1696 1697 REAL(KIND=dp) :: MassFlux 1698 1699 INTEGER :: i,j,k,l,t,q,p,c,dim,N_Integ 1700 1701 LOGICAL :: stat, Found 1702 1703 TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff 1704 1705!------------------------------------------------------------------------------ 1706! NormalTangential = .TRUE. 1707 dim = CoordinateSystemDimension() 1708 c = dim + 1 1709! 1710!------------------------------------------------------------------------------ 1711! Integration stuff 1712!------------------------------------------------------------------------------ 1713 IntegStuff = GaussPoints( element ) 1714 U_Integ => IntegStuff % u 1715 V_Integ => IntegStuff % v 1716 W_Integ => IntegStuff % w 1717 S_Integ => IntegStuff % s 1718 N_Integ = IntegStuff % n 1719 1720!------------------------------------------------------------------------------ 1721! Now we start integrating 1722!------------------------------------------------------------------------------ 1723 DO t=1,N_Integ 1724 1725 u = U_Integ(t) 1726 v = V_Integ(t) 1727 w = W_Integ(t) 1728!------------------------------------------------------------------------------ 1729! Basis function values & derivatives at the integration point 1730!------------------------------------------------------------------------------ 1731 stat = ElementInfo( Element,Nodes,u,v,w,detJ, & 1732 Basis, dBasisdx ) 1733 1734 s = detJ * S_Integ(t) 1735!------------------------------------------------------------------------------ 1736! Add to load: tangential derivative of something 1737!------------------------------------------------------------------------------ 1738 DO i=1,dim 1739 TangentForce(i) = 0.0d0 1740 END DO 1741 1742!------------------------------------------------------------------------------ 1743! Add to load: given force in coordinate directions 1744!------------------------------------------------------------------------------ 1745 Force = 0.0d0 1746 DO i=1,dim 1747 Force(i) = Force(i) + SUM( LoadVector(i,1:n)*Basis ) 1748 END DO 1749 1750!------------------------------------------------------------------------------ 1751! Add to load: given force in normal direction 1752!------------------------------------------------------------------------------ 1753 Normal = NormalVector( Element, Nodes, u,v,.TRUE. ) 1754 1755 Alpha = SUM( NodalExtPressure(1:n) * Basis ) 1756 IF ( NormalTangential ) THEN 1757 Force(1) = Force(1) + Alpha 1758 ELSE 1759 DO i=1,dim 1760 Force(i) = Force(i) + Alpha * Normal(i) 1761 END DO 1762 END IF 1763 1764!------------------------------------------------------------------------------ 1765 1766 Alpha = 0.0d0 1767 MassFlux = 0.0d0 1768 1769!------------------------------------------------------------------------------ 1770 1771 SELECT CASE( Element % TYPE % DIMENSION ) 1772 CASE(1) 1773 Tangent(1) = Normal(2) 1774 Tangent(2) = -Normal(1) 1775 Tangent(3) = 0.0_dp 1776 Tangent2 = 0.0_dp 1777 CASE(2) 1778 CALL TangentDirections( Normal, Tangent, Tangent2 ) 1779 END SELECT 1780 1781 IF ( ANY( NodalSlipCoeff(:,:) /= 0.0d0 ) ) THEN 1782 DO p=1,nd 1783 DO q=1,nd 1784 DO i=1,dim 1785 SlipCoeff = SUM( NodalSlipCoeff(i,1:n) * Basis(1:n) ) 1786 1787 IF ( NormalTangential ) THEN 1788 SELECT CASE(i) 1789 CASE(1) 1790 Vect = Normal 1791 CASE(2) 1792 Vect = Tangent 1793 CASE(3) 1794 Vect = Tangent2 1795 END SELECT 1796 1797 DO j=1,dim 1798 DO k=1,dim 1799 BoundaryMatrix( (p-1)*c+j,(q-1)*c+k ) = & 1800 BoundaryMatrix( (p-1)*c+j,(q-1)*c+k ) + & 1801 s * SlipCoeff * Basis(q) * Basis(p) * Vect(j) * Vect(k) 1802 END DO 1803 END DO 1804 ELSE 1805 BoundaryMatrix( (p-1)*c+i,(q-1)*c+i ) = & 1806 BoundaryMatrix( (p-1)*c+i,(q-1)*c+i ) + & 1807 s * SlipCoeff * Basis(q) * Basis(p) 1808 END IF 1809 END DO 1810 END DO 1811 END DO 1812 END IF 1813 1814 DO q=1,nd 1815 DO i=1,dim 1816 k = (q-1)*c + i 1817 IF ( NormalTangential ) THEN 1818 SELECT CASE(i) 1819 CASE(1) 1820 Vect = Normal 1821 CASE(2) 1822 Vect = Tangent 1823 CASE(3) 1824 Vect = Tangent2 1825 END SELECT 1826 1827 DO j=1,dim 1828 l = (q-1)*c + j 1829 BoundaryVector(l) = BoundaryVector(l) + & 1830 s * Basis(q) * Force(i) * Vect(j) 1831 END DO 1832 ELSE 1833 BoundaryVector(k) = BoundaryVector(k) + s*Basis(q)*Force(i) 1834 END IF 1835 BoundaryVector(k) = BoundaryVector(k) - s * Alpha * dBasisdx(q,i) 1836 BoundaryVector(k) = BoundaryVector(k) + s * TangentForce(i)*Basis(q) 1837 END DO 1838 BoundaryVector(q*c) = BoundaryVector(q*c) + s * MassFlux * Basis(q) 1839 END DO 1840 1841 END DO 1842!------------------------------------------------------------------------------ 1843 END SUBROUTINE NavierStokesBoundary 1844!------------------------------------------------------------------------------ 1845 1846!------------------------------------------------------------------------------ 1847SUBROUTINE ComputeVarLoads(Solver) 1848!------------------------------------------------------------------------------ 1849 1850!------------------------------------------------------------------------------ 1851 TYPE(Solver_t), TARGET :: Solver 1852!------------------------------------------------------------------------------ 1853 TYPE(Matrix_t), POINTER :: Aaid, Projector 1854 TYPE(Variable_t), POINTER :: NodalLoads 1855 1856 REAL(KIND=dp), POINTER CONTIG :: SaveValues(:) 1857 REAL(KIND=dp), ALLOCATABLE :: x(:),TempVector(:), TempRHS(:) 1858 1859 INTEGER :: DOFs 1860 INTEGER :: i,ii,j,l,DOF,This 1861 1862!------------------------------------------------------------------------------ 1863 1864 1865 DOFs=Solver % Variable % DOFs 1866 1867 NodalLoads => VariableGet( Solver % Mesh % Variables, & 1868 GetVarName(Solver % Variable) // ' Loads' ) 1869 1870 Aaid => Solver % Matrix 1871 1872 IF ( ASSOCIATED(NodalLoads) .AND. ASSOCIATED(Aaid % BulkValues) ) THEN 1873 ALLOCATE(x(SIZE(Solver % Variable % Values))) 1874 x(:)=Solver % Variable % Values(:) 1875 ALLOCATE( TempVector(Aaid % NumberOfRows) ) 1876 1877 SaveValues => Aaid % Values 1878 Aaid % Values => Aaid % BulkValues 1879 1880 IF ( ParEnv % PEs > 1 ) THEN 1881 ALLOCATE(TempRHS(SIZE(AAid % BulkRHS))) 1882 TempRHS = Aaid % BulkRHS 1883 CALL ParallelInitSolve( Aaid, x, TempRHS, Tempvector ) 1884 CALL ParallelMatrixVector( Aaid, x, TempVector, .TRUE. ) 1885 ELSE 1886 CALL MatrixVectorMultiply( Aaid, x, TempVector ) 1887 END IF 1888 1889 1890 Aaid % Values => SaveValues 1891 IF ( ParEnv % PEs>1 ) THEN 1892 DO i=1,Aaid % NumberOfRows 1893 IF ( AAid % ParallelInfo % NeighbourList(i) % Neighbours(1) == ParEnv % Mype ) THEN 1894 TempVector(i) = TempVector(i)-TempRHS(i) 1895 ELSE 1896 TempVector(i) = 0 1897 END IF 1898 END DO 1899 DEALLOCATE(TempRHS) 1900 CALL ParallelSumVector( AAid, Tempvector ) 1901 ELSE 1902 TempVector = TempVector - Aaid % BulkRHS 1903 END IF 1904 1905 DO This=1,CurrentModel % NumberOfBCs 1906 Projector=>CurrentModel % BCs(This) % PMatrix 1907 IF (ASSOCIATED(Projector))THEN 1908 DO DOF=1,DOFs 1909 DO i=1,Projector % NumberOfRows 1910 ii = Projector % InvPerm(i) 1911 k = Solver % Variable % Perm(ii) 1912 IF(k<=0) CYCLE 1913 k = DOFs * (k-1) + DOF 1914 TempVector(k)=0 1915 1916 DO l = Projector % Rows(i), Projector % Rows(i+1)-1 1917 IF ( Projector % Cols(l) <= 0 ) CYCLE 1918 m = Solver % Variable % Perm( Projector % Cols(l) ) 1919 IF ( m > 0 ) THEN 1920 m = DOFs * (m-1) + DOF 1921 TempVector(k) = TempVector(k) + Projector % Values(l)*TempVector(m) 1922 END IF 1923 END DO 1924 END DO 1925 END DO 1926 END IF 1927 END DO 1928 1929 DO i=1,SIZE( NodalLoads % Perm ) 1930 IF ( NodalLoads % Perm(i)>0 .AND. Solver % Variable % Perm(i)>0 ) THEN 1931 DO j=1,DOFs 1932 NodalLoads % Values(DOFs*(NodalLoads % Perm(i)-1)+j) = & 1933 TempVector(DOFs*(Solver % Variable % Perm(i)-1)+j) 1934 END DO 1935 END IF 1936 END DO 1937 DEALLOCATE( x,TempVector ) 1938 1939 CALL BackRotateNTSystem(NodalLoads % Values,NodalLoads % Perm,DOFs) 1940 1941 END IF 1942 1943!------------------------------------------------------------------------------ 1944 END SUBROUTINE ComputeVarLoads 1945 1946!------------------------------------------------------------------------------ 1947 1948END SUBROUTINE StokesSolver 1949!------------------------------------------------------------------------------ 1950