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, Mikko Lyly 27! * Email: Juha.Ruokolainen@csc.fi 28! * Web: http://www.csc.fi/elmer 29! * Address: CSC - IT Center for Science Ltd. 30! * Keilaranta 14 31! * 02101 Espoo, Finland 32! * 33! * Original Date: 04 Oct 2000 34! * 35! ****************************************************************************/ 36 37!------------------------------------------------------------------------------ 38!> Solve the complex diffusion-convection-reaction equation. 39!> \ingroup Solvers 40!------------------------------------------------------------------------------ 41SUBROUTINE DCRComplexSolver( Model,Solver,dt,TransientSimulation ) 42!------------------------------------------------------------------------------ 43 USE Types 44 USE Lists 45 USE Adaptive 46 USE Integration 47 USE ElementDescription 48 USE SolverUtils 49 50 IMPLICIT NONE 51!------------------------------------------------------------------------------ 52 TYPE(Solver_t) :: Solver 53 TYPE(Model_t) :: Model 54 55 REAL(KIND=dp) :: dt 56 LOGICAL :: TransientSimulation 57!------------------------------------------------------------------------------ 58! Local variables 59!------------------------------------------------------------------------------ 60 TYPE(Matrix_t),POINTER :: StiffMatrix 61 TYPE(Nodes_t) :: ElementNodes 62 TYPE(Element_t),POINTER :: CurrentElement 63 64 INTEGER, POINTER :: NodeIndexes(:) 65 66 LOGICAL :: AllocationsDone = .FALSE., Bubbles, GotIt, notScalar = .TRUE., stat 67 68 INTEGER, POINTER :: PressurePerm(:) 69 REAL(KIND=dp), POINTER :: Pressure(:), ForceVector(:) 70 71 INTEGER :: iter, i, j, k, n, t, istat, eq, LocalNodes 72 REAL(KIND=dp) :: Norm, PrevNorm, RelativeChange 73 74 TYPE(ValueList_t), POINTER :: Material 75 76 INTEGER :: NonlinearIter 77 REAL(KIND=dp) :: NonlinearTol,s 78 79 REAL(KIND=dp), ALLOCATABLE :: LocalStiffMatrix(:,:), Load(:,:), Work(:), & 80 LocalForce(:), & 81 Amatrix(:,:,:), AvectorReal(:,:), AvectorImag(:,:), AscalarReal(:), & 82 AscalarImag(:), Bvector(:,:), BscalarReal(:), BscalarImag(:) 83 84 CHARACTER(LEN=MAX_NAME_LEN) :: EquationName 85 86 SAVE LocalStiffMatrix, Work, Load, LocalForce, ElementNodes, & 87 AllocationsDone, & 88 Amatrix, AvectorReal, AvectorImag, AscalarReal, AscalarImag, & 89 Bvector, BscalarReal, BscalarImag 90 REAL(KIND=dp) :: at,at0,totat,st,totst,t1 91!------------------------------------------------------------------------------ 92 INTERFACE 93 FUNCTION DCRBoundaryResidual( Model,Edge,Mesh,Quant,Perm,Gnorm ) RESULT(Indicator) 94 USE Types 95 TYPE(Element_t), POINTER :: Edge 96 TYPE(Model_t) :: Model 97 TYPE(Mesh_t), POINTER :: Mesh 98 REAL(KIND=dp) :: Quant(:), Indicator(2), Gnorm 99 INTEGER :: Perm(:) 100 END FUNCTION DCRBoundaryResidual 101 102 FUNCTION DCREdgeResidual( Model,Edge,Mesh,Quant,Perm ) RESULT(Indicator) 103 USE Types 104 TYPE(Element_t), POINTER :: Edge 105 TYPE(Model_t) :: Model 106 TYPE(Mesh_t), POINTER :: Mesh 107 REAL(KIND=dp) :: Quant(:), Indicator(2) 108 INTEGER :: Perm(:) 109 END FUNCTION DCREdgeResidual 110 111 FUNCTION DCRInsideResidual( Model,Element,Mesh,Quant,Perm, Fnorm ) RESULT(Indicator) 112 USE Types 113 TYPE(Element_t), POINTER :: Element 114 TYPE(Model_t) :: Model 115 TYPE(Mesh_t), POINTER :: Mesh 116 REAL(KIND=dp) :: Quant(:), Indicator(2), Fnorm 117 INTEGER :: Perm(:) 118 END FUNCTION DCRInsideResidual 119 END INTERFACE 120 121!------------------------------------------------------------------------------ 122! Get variables needed for solution 123!------------------------------------------------------------------------------ 124 IF ( .NOT.ASSOCIATED( Solver % Matrix ) ) RETURN 125 Solver % Matrix % Complex = .TRUE. 126 127 Pressure => Solver % Variable % Values 128 PressurePerm => Solver % Variable % Perm 129 130 LocalNodes = COUNT( PressurePerm > 0 ) 131 IF ( LocalNodes <= 0 ) RETURN 132 133 StiffMatrix => Solver % Matrix 134 ForceVector => StiffMatrix % RHS 135 136 Norm = Solver % Variable % Norm 137 138!------------------------------------------------------------------------------ 139! Allocate some permanent storage, this is done first time only 140!------------------------------------------------------------------------------ 141 IF ( .NOT. AllocationsDone .OR. Solver % MeshChanged ) THEN 142 N = Solver % Mesh % MaxElementNodes 143 144 IF ( AllocationsDone ) THEN 145 DEALLOCATE( & 146 ElementNodes % x, & 147 ElementNodes % y, & 148 ElementNodes % z, & 149 LocalForce, & 150 Work, & 151 LocalStiffMatrix, & 152 Load, & 153 Amatrix, AvectorReal, AvectorImag, AscalarReal, AscalarImag, & 154 Bvector, BscalarReal, BscalarImag ) 155 END IF 156 157 ALLOCATE( ElementNodes % x( N ), & 158 ElementNodes % y( N ), & 159 ElementNodes % z( N ), & 160 LocalForce( 2*N ), & 161 Work( N ), & 162 LocalStiffMatrix( 2*N,2*N ), & 163 Amatrix(3,3,N), AvectorReal(3,N), AvectorImag(3,N), & 164 AscalarReal(N), AscalarImag(N), Bvector(3,N), & 165 BscalarReal(N), BscalarImag(N), & 166 Load( 2,N ), STAT=istat ) 167 168 IF ( istat /= 0 ) THEN 169 CALL Fatal( 'DCRComplexSolve', 'Memory allocation error.' ) 170 END IF 171 172 AllocationsDone = .TRUE. 173 END IF 174!------------------------------------------------------------------------------ 175 176!------------------------------------------------------------------------------ 177! Do some additional initialization, and go for it 178!------------------------------------------------------------------------------ 179 NonlinearTol = ListGetConstReal( Solver % Values, & 180 'Nonlinear System Convergence Tolerance',GotIt ) 181 182 NonlinearIter = ListGetInteger( Solver % Values, & 183 'Nonlinear System Max Iterations', GotIt ) 184 185 IF ( .NOT.GotIt ) NonlinearIter = 1 186 187 EquationName = ListGetString( Solver % Values, 'Equation' ) 188 189 Bubbles = ListGetLogical( Solver % Values, 'Bubbles', GotIt ) 190 191!------------------------------------------------------------------------------ 192! Iterate over any nonlinearity of material or source 193!------------------------------------------------------------------------------ 194 totat = 0.0d0 195 totst = 0.0d0 196 197 DO iter=1,NonlinearIter 198!------------------------------------------------------------------------------ 199 at = CPUTime() 200 at0 = RealTime() 201 202 CALL Info( 'DCRComplexSolve', ' ', Level=4 ) 203 CALL Info( 'DCRComplexSolve', '-------------------------------------', Level=4 ) 204 WRITE( Message, * ) 'DCRComplex iteration', iter 205 CALL Info( 'DCRComplexSolve', Message, Level=4 ) 206 CALL Info( 'DCRComplexSolve', '-------------------------------------', Level=4 ) 207 CALL Info( 'DCRComplexSolve', ' ', Level=4 ) 208 CALL Info( 'DCRComplexSolve', 'Starting Assmebly', Level=4 ) 209 210 CALL InitializeToZero( StiffMatrix, ForceVector ) 211! 212! Do the bulk assembly: 213! --------------------- 214 215!------------------------------------------------------------------------------ 216 DO t=1,Solver % NumberOfActiveElements 217!------------------------------------------------------------------------------ 218 IF ( RealTime() - at0 > 1.0 ) THEN 219 WRITE(Message,'(a,i3,a)' ) ' Assembly: ', INT(100.0 - 100.0 * & 220 (Solver % NumberOfActiveElements-t) / & 221 (1.0*Solver % NumberOfActiveElements)), ' % done' 222 CALL Info( 'DCRComplexSolve', Message, Level=5 ) 223 224 at0 = RealTime() 225 END IF 226!------------------------------------------------------------------------------ 227! Check if this element belongs to a body where this equation 228! should be computed 229!------------------------------------------------------------------------------ 230 CurrentElement => Solver % Mesh % Elements(Solver % ActiveElements(t)) 231 232! IF ( .NOT. CheckElementEquation( Model, & 233! CurrentElement, EquationName ) ) CYCLE 234!------------------------------------------------------------------------------ 235 Model % CurrentElement => CurrentElement 236 237 n = CurrentElement % Type % NumberOfNodes 238 NodeIndexes => CurrentElement % NodeIndexes 239 240 ElementNodes % x(1:n) = Solver % Mesh % Nodes % x(NodeIndexes) 241 ElementNodes % y(1:n) = Solver % Mesh % Nodes % y(NodeIndexes) 242 ElementNodes % z(1:n) = Solver % Mesh % Nodes % z(NodeIndexes) 243 244!------------------------------------------------------------------------------ 245! Get equation & material parameters 246!------------------------------------------------------------------------------ 247 k = ListGetInteger( Model % Bodies( CurrentElement % & 248 Bodyid ) % Values, 'Material',minv=1,maxv=Model % NumberOfMaterials ) 249 250 Material => Model % Materials(k) % Values 251 252!------------------------------------------------------------------------------ 253! Second and first order time derivative term coefficients on nodes 254!------------------------------------------------------------------------------ 255 CALL InputTensor( Amatrix, notScalar, & 256 'Amatrix', Material, n, NodeIndexes ) 257 258 CALL InputVector( AvectorReal, notScalar, & 259 'Avector 1', Material, n, NodeIndexes ) 260 261 CALL InputVector( AvectorImag, notScalar, & 262 'Avector 2', Material, n, NodeIndexes ) 263 264 AscalarReal(1:n) = ListGetReal( Material, & 265 'Ascalar 1', n, NodeIndexes, GotIt) 266 267 AscalarImag(1:n) = ListGetReal( Material, & 268 'Ascalar 2', n, NodeIndexes, GotIt) 269 270!------------------------------------------------------------------------------ 271! The source term on nodes 272!------------------------------------------------------------------------------ 273 k = ListGetInteger( Model % Bodies( CurrentElement % BodyId ) % & 274 Values, 'Body Force', GotIt, 1, Model % NumberOFBodyForces ) 275 276 Load = 0.0d0 277 IF ( k > 0 ) THEN 278 Load(1,1:n) = ListGetReal( Model % BodyForces(k) % Values, & 279 'Pressure Source 1', n, NodeIndexes, GotIt ) 280 281 Load(2,1:n) = ListGetReal( Model % BodyForces(k) % Values, & 282 'Pressure Source 2', n, NodeIndexes, GotIt ) 283 END IF 284 285!------------------------------------------------------------------------------ 286! Get element local matrix and rhs vector 287!------------------------------------------------------------------------------ 288 CALL LocalMatrix( LocalStiffMatrix, LocalForce, & 289 Load, Bubbles, CurrentElement, n, ElementNodes, & 290 Amatrix, AvectorReal, AvectorImag, AscalarReal, AscalarImag ) 291 292!------------------------------------------------------------------------------ 293! Update global matrix and rhs vector from local matrix & vector 294!------------------------------------------------------------------------------ 295 CALL UpdateGlobalEquations( StiffMatrix, LocalStiffMatrix, & 296 ForceVector, LocalForce, n, Solver % Variable % DOFs, & 297 PressurePerm(NodeIndexes) ) 298!------------------------------------------------------------------------------ 299 END DO 300!------------------------------------------------------------------------------ 301 302! 303! Neumann & Newton BCs: 304! --------------------- 305 306!------------------------------------------------------------------------------ 307 DO t = Solver % Mesh % NumberOfBulkElements + 1, & 308 Solver % Mesh % NumberOfBulkElements + & 309 Solver % Mesh % NumberOfBoundaryElements 310!------------------------------------------------------------------------------ 311 CurrentElement => Solver % Mesh % Elements(t) 312 Model % CurrentElement => CurrentElement 313 314!------------------------------------------------------------------------------ 315! The element type 101 (point element) can only be used 316! to set Dirichlet BCs, so skip em at this stage. 317!------------------------------------------------------------------------------ 318 IF ( CurrentElement % Type % ElementCode == 101 ) CYCLE 319 320!------------------------------------------------------------------------------ 321 DO i=1,Model % NumberOfBCs 322 IF ( CurrentElement % BoundaryInfo % Constraint == & 323 Model % BCs(i) % Tag ) THEN 324!------------------------------------------------------------------------------ 325 n = CurrentElement % Type % NumberOfNodes 326 NodeIndexes => CurrentElement % NodeIndexes 327 328 IF ( ANY( PressurePerm(NodeIndexes) == 0 ) ) CYCLE 329 330 ElementNodes % x(1:n) = Solver % Mesh % Nodes % x(NodeIndexes) 331 ElementNodes % y(1:n) = Solver % Mesh % Nodes % y(NodeIndexes) 332 ElementNodes % z(1:n) = Solver % Mesh % Nodes % z(NodeIndexes) 333 334 Load(1,1:n) = ListGetReal( Model % BCs(i) % Values, & 335 'Wave Flux 1', n, NodeIndexes, GotIt ) 336 337 Load(2,1:n) = ListGetReal( Model % BCs(i) % Values, & 338 'Wave Flux 2', n, NodeIndexes, GotIt ) 339 340 CALL InputVector( Bvector, notScalar, & 341 'Bvector', Model % BCs(i) % Values, n, NodeIndexes ) 342 343 BscalarReal(1:n) = ListGetReal( Model % BCs(i) % Values, & 344 'Bscalar 1', n, NodeIndexes, GotIt) 345 346 BscalarImag(1:n) = ListGetReal( Model % BCs(i) % Values, & 347 'Bscalar 2', n, NodeIndexes, GotIt) 348 349!------------------------------------------------------------------------------ 350! Get element local matrix and rhs vector 351!------------------------------------------------------------------------------ 352 CALL LocalMatrixBoundary( LocalStiffMatrix, LocalForce, & 353 Load, CurrentElement, n, ElementNodes, & 354 Amatrix, AvectorReal, AvectorImag, AscalarReal, AscalarImag, & 355 Bvector, BscalarReal, BscalarImag ) 356 357!------------------------------------------------------------------------------ 358! Update global matrix and rhs vector from local matrix & vector 359!------------------------------------------------------------------------------ 360 CALL UpdateGlobalEquations( StiffMatrix, LocalStiffMatrix, & 361 ForceVector, LocalForce, n, Solver % Variable % DOFs, & 362 PressurePerm(NodeIndexes) ) 363!------------------------------------------------------------------------------ 364 END IF 365 END DO 366!------------------------------------------------------------------------------ 367 END DO 368!------------------------------------------------------------------------------ 369 370 CALL FinishAssembly( Solver, ForceVector ) 371! 372! Dirichlet BCs: 373! -------------- 374 CALL SetDirichletBoundaries( Model, StiffMatrix, ForceVector, & 375 ComponentName(Solver % Variable,1), 1, & 376 Solver % Variable % DOFs, PressurePerm ) 377 378 CALL SetDirichletBoundaries( Model, StiffMatrix, ForceVector, & 379 ComponentName(Solver % Variable,2), 2, & 380 Solver % Variable % DOFs, PressurePerm ) 381 382 CALL Info( 'DCRComplexSolve', 'Assembly done', Level=4 ) 383 384! 385! Solve the system and we are done: 386! --------------------------------- 387 PrevNorm = Norm 388 at = CPUTime() - at 389 st = CPUTime() 390 391 CALL SolveSystem( StiffMatrix, ParMatrix, ForceVector, & 392 Pressure, Norm, Solver % Variable % DOFs, Solver ) 393 394 st = CPUTIme()-st 395 totat = totat + at 396 totst = totst + st 397 WRITE(Message,'(a,i4,a,F8.2,F8.2)') 'iter: ',iter,' Assembly: (s)', at, totat 398 CALL Info( 'DCRComplexSolve', Message, Level=4 ) 399 WRITE(Message,'(a,i4,a,F8.2,F8.2)') 'iter: ',iter,' Solve: (s)', st, totst 400 CALL Info( 'DCRComplexSolve', Message, Level=4 ) 401 402!------------------------------------------------------------------------------ 403 IF ( PrevNorm + Norm /= 0.0d0 ) THEN 404 RelativeChange = 2*ABS(PrevNorm - Norm) / (PrevNorm + Norm) 405 ELSE 406 RelativeChange = 0.0d0 407 END IF 408 409 CALL Info( 'DCRComplexSolve', ' ', Level=4 ) 410 WRITE( Message, * ) 'Result Norm : ',Norm 411 CALL Info( 'DCRComplexSolve', Message, Level=4 ) 412 WRITE( Message, * ) 'Relative Change: ',RelativeChange 413 CALL Info( 'DCRComplexSolve', Message, Level=4 ) 414 415 IF ( RelativeChange < NonlinearTol ) EXIT 416!------------------------------------------------------------------------------ 417 END DO ! of nonlinear iteration 418!------------------------------------------------------------------------------ 419 IF ( ListGetLogical( Solver % Values, 'Adaptive Mesh Refinement', GotIt ) ) & 420 CALL RefineMesh( Model,Solver,Pressure,PressurePerm, & 421 DCRInsideResidual, DCREdgeResidual, DCRBoundaryResidual ) 422 423CONTAINS 424 425!------------------------------------------------------------------------------ 426 SUBROUTINE InputTensor( Tensor, IsScalar, Name, Material, n, NodeIndexes ) 427!------------------------------------------------------------------------------ 428 REAL(KIND=dp) :: Tensor(:,:,:) 429 INTEGER :: n, NodeIndexes(:) 430 LOGICAL :: IsScalar 431 CHARACTER(LEN=*) :: Name 432 TYPE(ValueList_t), POINTER :: Material 433!------------------------------------------------------------------------------ 434 LOGICAL :: FirstTime = .TRUE., stat 435 REAL(KIND=dp), POINTER :: Hwrk(:,:,:) 436 437 SAVE FirstTime, Hwrk 438!------------------------------------------------------------------------------ 439 IF ( FirstTime ) THEN 440 NULLIFY( Hwrk ) 441 FirstTime = .FALSE. 442 END IF 443 444 Tensor = 0.0d0 445 446 CALL ListGetRealArray( Material, Name, Hwrk, n, NodeIndexes, stat ) 447 IsScalar = SIZE(HWrk,1) == 1 .AND. SIZE(HWrk,2) == 1 448 449 IF ( .NOT. stat ) RETURN 450 451 IF ( SIZE(Hwrk,1) == 1 ) THEN 452 453 DO i=1,MIN(3,SIZE(Hwrk,2)) 454 Tensor( i,i,1:n ) = Hwrk( 1,1,1:n ) 455 END DO 456 457 ELSE IF ( SIZE(Hwrk,2) == 1 ) THEN 458 459 DO i=1,MIN(3,SIZE(Hwrk,1)) 460 Tensor(i,i,1:n) = Hwrk(i,1,1:n) 461 END DO 462 463 ELSE 464 465 DO i=1,MIN(3,SIZE(Hwrk,1)) 466 DO j=1,MIN(3,SIZE(Hwrk,2)) 467 Tensor( i,j,1:n ) = Hwrk(i,j,1:n) 468 END DO 469 END DO 470 471 END IF 472!------------------------------------------------------------------------------ 473 END SUBROUTINE InputTensor 474!------------------------------------------------------------------------------ 475 476 477!------------------------------------------------------------------------------ 478 SUBROUTINE InputVector( Tensor, IsScalar, Name, Material, n, NodeIndexes ) 479!------------------------------------------------------------------------------ 480 REAL(KIND=dp) :: Tensor(:,:) 481 INTEGER :: n, NodeIndexes(:) 482 LOGICAL :: IsScalar 483 CHARACTER(LEN=*) :: Name 484 TYPE(ValueList_t), POINTER :: Material 485!------------------------------------------------------------------------------ 486 LOGICAL :: FirstTime = .TRUE., stat 487 REAL(KIND=dp), POINTER :: Hwrk(:,:,:) 488 489 SAVE FirstTime, Hwrk 490!------------------------------------------------------------------------------ 491 IF ( FirstTime ) THEN 492 NULLIFY( Hwrk ) 493 FirstTime = .FALSE. 494 END IF 495 496 Tensor = 0.0d0 497 498 CALL ListGetRealArray( Material, Name, Hwrk, n, NodeIndexes, stat ) 499 IsScalar = SIZE(HWrk,1) == 1 .AND. SIZE(HWrk,2) == 1 500 501 IF ( .NOT. stat ) RETURN 502 503 IF ( SIZE(Hwrk,1) == 1 ) THEN 504 505 DO i=1,MIN(3,SIZE(Hwrk,2)) 506 Tensor( i,1:n ) = Hwrk( 1,1,1:n ) 507 END DO 508 509 ELSE 510 511 DO i=1,MIN(3,SIZE(Hwrk,1)) 512 Tensor( i,1:n ) = Hwrk( i,1,1:n ) 513 END DO 514 515 END IF 516!------------------------------------------------------------------------------ 517 END SUBROUTINE InputVector 518!------------------------------------------------------------------------------ 519 520 521!------------------------------------------------------------------------------ 522 523 SUBROUTINE LocalMatrix( StiffMatrix, Force, Load, Bubbles, Element, n, & 524 Nodes, Amatrix, AvectorReal, AvectorImag, AscalarReal, AscalarImag ) 525!------------------------------------------------------------------------------ 526 REAL(KIND=dp) :: StiffMatrix(:,:), Force(:), Load(:,:), & 527 Amatrix(:,:,:), AvectorReal(:,:), AvectorImag(:,:), & 528 AscalarReal(:), AscalarImag(:) 529 LOGICAL :: Bubbles 530 INTEGER :: n 531 TYPE(Nodes_t) :: Nodes 532 TYPE(Element_t), POINTER :: Element 533!------------------------------------------------------------------------------ 534 REAL(KIND=dp) :: Basis(2*n),dBasisdx(2*n,3) 535 REAL(KIND=dp) :: SqrtElementMetric,U,V,W,S,M,D,L1,L2 536 REAL(KIND=dp) :: DiffCoef(3,3), Velo(3) 537 COMPLEX(KIND=dp) :: LSTIFF(2*n,2*n), LFORCE(2*n), A 538 LOGICAL :: Stat 539 INTEGER :: i,p,q,t,dim, NBasis, CoordSys 540 TYPE(GaussIntegrationPoints_t) :: IntegStuff 541 542 REAL(KIND=dp) :: X,Y,Z,Metric(3,3),SqrtMetric,Symb(3,3,3),dSymb(3,3,3,3) 543 REAL(kind=dp) :: A2(3,3), A1r(3), A1i(3), A0r, A0i 544!------------------------------------------------------------------------------ 545 dim = CoordinateSystemDimension() 546 CoordSys = CurrentCoordinateSystem() 547 548 Metric = 0.0d0 549 Metric(1,1) = 1.0d0 550 Metric(2,2) = 1.0d0 551 Metric(3,3) = 1.0d0 552 553 LSTIFF = 0.0d0 554 LFORCE = 0.0d0 555!------------------------------------------------------------------------------ 556! Numerical integration 557!------------------------------------------------------------------------------ 558 IF ( Bubbles ) THEN 559 IntegStuff = GaussPoints( Element, Element % Type % GaussPoints2 ) 560 NBasis = 2*n 561 ELSE 562 NBasis = n 563 IntegStuff = GaussPoints( Element ) 564 END IF 565!------------------------------------------------------------------------------ 566 DO t=1,IntegStuff % n 567 U = IntegStuff % u(t) 568 V = IntegStuff % v(t) 569 W = IntegStuff % w(t) 570 S = IntegStuff % s(t) 571!------------------------------------------------------------------------------ 572! Basis function values & derivatives at the integration point 573!------------------------------------------------------------------------------ 574 stat = ElementInfo( Element, Nodes, U, V, W, SqrtElementMetric, & 575 Basis, dBasisdx, Bubbles=Bubbles ) 576 577 s = s * SqrtElementMetric 578 IF ( CoordSys /= Cartesian ) THEN 579 X = SUM( Nodes % X(1:n) * Basis(1:n) ) 580 Y = SUM( Nodes % Y(1:n) * Basis(1:n) ) 581 Z = SUM( Nodes % Z(1:n) * Basis(1:n) ) 582 CALL CoordinateSystemInfo( Metric,SqrtMetric,Symb,dSymb,X,Y,Z ) 583 s = s * SqrtMetric 584 END IF 585!------------------------------------------------------------------------------ 586! The source term and the coefficient of the time derivative and 587! diffusion terms at the integration point 588!------------------------------------------------------------------------------ 589! D = WaveNumber * SUM( Damping(1:n) * Basis(1:n) ) 590! M = -WaveNumber**2 591 592 L1 = SUM( Load(1,1:n) * Basis(1:n) ) 593 L2 = SUM( Load(2,1:n) * Basis(1:n) ) 594 595 A2 = 0.0d0 596 A1r = 0.0d0 597 A1i = 0.0d0 598 A0r = 0.0d0 599 A0i = 0.0d0 600 601 A0r = SUM( AscalarReal(1:n) * Basis(1:n) ) 602 A0i = SUM( AscalarImag(1:n) * Basis(1:n) ) 603 do i = 1,dim 604 A1r(i) = A1r(i) + SUM( AvectorReal(i,1:n) * Basis(1:n) ) 605 A1i(i) = A1i(i) + SUM( AvectorImag(i,1:n) * Basis(1:n) ) 606 do j = 1,dim 607 A2(i,j) = A2(i,j) + SUM( Amatrix(i,j,1:n) * Basis(1:n) ) 608 end do 609 end do 610 611 612! Stiffness matrix and load vector 613! -------------------------------- 614 DO p=1,NBasis 615 DO q=1,NBasis 616 A = CMPLX( A0r, A0i,KIND=dp ) * Basis(q) * Basis(p) 617 DO i=1,dim 618 A = A + CMPLX( A1r(i), A1i(i), KIND=dp ) * dBasisdx(q,i) * basis(p) 619 DO j=1,dim 620 DO k = 1,dim 621 A = A + Metric(i,j) * A2(i,k) * dBasisdx(q,k) * dBasisdx(p,j) 622 END DO 623 END DO 624 END DO 625 LSTIFF(p,q) = LSTIFF(p,q) + s*A 626 END DO 627 LFORCE(p) = LFORCE(p) + s * Basis(p) * CMPLX( L1,L2,KIND=dp ) 628 END DO 629 END DO 630!------------------------------------------------------------------------------ 631 632 IF ( Bubbles ) THEN 633 CALL CondensateP( n, n, LSTIFF, LFORCE ) 634 END IF 635 636 DO i=1,n 637 Force( 2*(i-1)+1 ) = REAL( LFORCE(i) ) 638 Force( 2*(i-1)+2 ) = AIMAG( LFORCE(i) ) 639 640 DO j=1,n 641 StiffMatrix( 2*(i-1)+1, 2*(j-1)+1 ) = REAL( LSTIFF(i,j) ) 642 StiffMatrix( 2*(i-1)+1, 2*(j-1)+2 ) = -AIMAG( LSTIFF(i,j) ) 643 StiffMatrix( 2*(i-1)+2, 2*(j-1)+1 ) = AIMAG( LSTIFF(i,j) ) 644 StiffMatrix( 2*(i-1)+2, 2*(j-1)+2 ) = REAL( LSTIFF(i,j) ) 645 END DO 646 END DO 647!------------------------------------------------------------------------------ 648 END SUBROUTINE LocalMatrix 649!------------------------------------------------------------------------------ 650 651 652!------------------------------------------------------------------------------ 653 SUBROUTINE LocalMatrixBoundary( StiffMatrix, Force, & 654 Load, Element, n, Nodes, & 655 Amatrix, AvectorReal, AvectorImag, AscalarReal, AscalarImag, & 656 Bvector, BscalarReal, BscalarImag ) 657!------------------------------------------------------------------------------ 658 REAL(KIND=dp) :: StiffMatrix(:,:),Force(:),Load(:,:) 659! REAL(KIND=dp) :: ConvVelo(:,:) 660 INTEGER :: n 661 REAL(kind=dp) :: Amatrix(:,:,:), AvectorReal(:,:), AvectorImag(:,:), & 662 AscalarReal(:), AscalarImag(:), Bvector(:,:), BscalarReal(:), & 663 BscalarImag(:) 664 TYPE(Nodes_t) :: Nodes 665 TYPE(Element_t), POINTER :: Element 666!------------------------------------------------------------------------------ 667 REAL(KIND=dp) :: SqrtElementMetric,U,V,W,S,L1,L2 668 REAL(KIND=dp) :: Basis(n),dBasisdx(n,3),X,Y,Z 669 REAL(KIND=dp) :: Normal(3), Velo(3), NormVelo, TangVelo(3) 670 REAL(kind=dp) :: A2(3,3), A1r(3), A1i(3), A0r, A0i, B1(3), B0r, B0i, & 671 C1(3), C0 672 COMPLEX(KIND=dp) :: LSTIFF(n,n), LFORCE(n), A 673 LOGICAL :: Stat 674 INTEGER :: i,p,q,t,dim,CoordSys 675 TYPE(GaussIntegrationPoints_t) :: IntegStuff 676!------------------------------------------------------------------------------ 677 dim = CoordinateSystemDimension() 678 CoordSys = CurrentCoordinateSystem() 679 680 LSTIFF = 0.0d0 681 LFORCE = 0.0d0 682!------------------------------------------------------------------------------ 683! Numerical integration 684!------------------------------------------------------------------------------ 685 IntegStuff = GaussPoints( Element ) 686!------------------------------------------------------------------------------ 687 DO t=1,IntegStuff % n 688 U = IntegStuff % u(t) 689 V = IntegStuff % v(t) 690 W = IntegStuff % w(t) 691 S = IntegStuff % s(t) 692!------------------------------------------------------------------------------ 693! Basis function values & derivatives at the integration point 694!------------------------------------------------------------------------------ 695 stat = ElementInfo( Element, Nodes, U, V, W, SqrtElementMetric, & 696 Basis, dBasisdx ) 697 698 s = s * SqrtElementMetric 699 700 Normal = Normalvector(Element, Nodes, U, V, .TRUE.) 701 702 A2 = 0.0d0 703 A1r = 0.0d0 704 A1i = 0.0d0 705 A0r = 0.0d0 706 A0i = 0.0d0 707 708 A0r = SUM( AscalarReal(1:n) * Basis(1:n) ) 709 A0i = SUM( AscalarImag(1:n) * Basis(1:n) ) 710 do i = 1,dim 711 A1r(i) = A1r(i) + SUM( AvectorReal(i,1:n) * Basis(1:n) ) 712 A1i(i) = A1i(i) + SUM( AvectorImag(i,1:n) * Basis(1:n) ) 713 do j = 1,dim 714 A2(i,j) = A2(i,j) + SUM( Amatrix(i,j,1:n) * Basis(1:n) ) 715 end do 716 end do 717 718 B1 = 0.0d0 719 B0r = 0.0d0 720 B0i = 0.0d0 721 722 B0r = SUM( BscalarReal(1:n) * Basis(1:n) ) 723 B0i = SUM( BscalarImag(1:n) * Basis(1:n) ) 724 do i = 1,dim 725 B1(i) = B1(i) + SUM( Bvector(i,1:n) * Basis(1:n) ) 726 end do 727 B1(1:dim) = B1(1:dim) + Normal(1:dim) 728 729 C1 = 0.0d0 730 C0 = 0.0d0 731 do i = 1,dim 732 do j = 1,dim 733 C0 = C0 + Normal(i)*A2(i,j)*Normal(j) 734 end do 735 end do 736 C0 = C0 / SUM( B1(1:dim) * Normal(1:dim) ) 737 C1 = C0*B1 738 do i = 1,dim 739 do j = 1,dim 740 C1(i) = C1(i) - A2(i,j)*Normal(j) 741 end do 742 end do 743 744 745!------------------------------------------------------------------------------ 746 L1 = SUM( Load(1,1:n) * Basis ) 747 L2 = SUM( Load(2,1:n) * Basis ) 748!------------------------------------------------------------------------------ 749 DO p=1,n 750 DO q=1,n 751 A = CMPLX( B0r, B0i, KIND=dp ) * C0 * Basis(p)*Basis(q) 752 A = A + SUM( C1(1:dim) * dBasisdx(q,1:dim) ) * Basis(p) 753 LSTIFF(p,q) = LSTIFF(p,q) + s * A 754 END DO 755 LFORCE(p) = LFORCE(p) + s * Basis(p) * C0 * CMPLX( L1, L2, KIND=dp ) 756 END DO 757!------------------------------------------------------------------------------ 758 END DO 759!------------------------------------------------------------------------------ 760 DO i=1,n 761 Force( 2*(i-1)+1 ) = REAL( LFORCE(i) ) 762 Force( 2*(i-1)+2 ) = AIMAG( LFORCE(i) ) 763 764 DO j=1,n 765 StiffMatrix( 2*(i-1)+1, 2*(j-1)+1 ) = REAL( LSTIFF(i,j) ) 766 StiffMatrix( 2*(i-1)+1, 2*(j-1)+2 ) = -AIMAG( LSTIFF(i,j) ) 767 StiffMatrix( 2*(i-1)+2, 2*(j-1)+1 ) = AIMAG( LSTIFF(i,j) ) 768 StiffMatrix( 2*(i-1)+2, 2*(j-1)+2 ) = REAL( LSTIFF(i,j) ) 769 END DO 770 END DO 771!------------------------------------------------------------------------------ 772 END SUBROUTINE LocalMatrixBoundary 773!------------------------------------------------------------------------------ 774 775!------------------------------------------------------------------------------ 776END SUBROUTINE DCRComplexSolver 777!------------------------------------------------------------------------------ 778 779!------------------------------------------------------------------------------ 780 FUNCTION DCRBoundaryResidual( Model, Edge, Mesh, Quant, Perm,Gnorm ) RESULT( Indicator ) 781!------------------------------------------------------------------------------ 782 USE DefUtils 783 IMPLICIT NONE 784!------------------------------------------------------------------------------ 785 TYPE(Model_t) :: Model 786 INTEGER :: Perm(:) 787 REAL(KIND=dp) :: Quant(:), Indicator(2), Gnorm 788 TYPE( Mesh_t ), POINTER :: Mesh 789 TYPE( Element_t ), POINTER :: Edge 790!------------------------------------------------------------------------------ 791 TYPE(Nodes_t) :: Nodes, EdgeNodes 792 TYPE(Element_t), POINTER :: Element, Bndry 793 794 INTEGER :: i,j,k,n,l,t,DIM,Pn,En 795 LOGICAL :: stat, GotIt, gotWIreal,gotWIimag 796 797 LOGICAL :: notScalar = .TRUE. 798 799 REAL(KIND=dp) :: SqrtMetric, Metric(3,3), Symb(3,3,3), dSymb(3,3,3,3) 800 REAL(KIND=dp) :: Grad(3,3), Normal(3), EdgeLength 801 REAL(KIND=dp) :: Source, ResidualNorm, Area 802 REAL(kind=dp) :: B1(3), B0r, B0i, Greal, Gimag 803 REAL(KIND=dp) :: u, v, w, s, detJ, ResidualReal, ResidualImag, WaveFlux(2) 804 805 REAL(KIND=dp), ALLOCATABLE :: EdgeBasis(:), Basis(:) 806 REAL(KIND=dp), ALLOCATABLE :: dBasisdx(:,:) 807 REAL(KIND=dp), ALLOCATABLE :: Flux(:) 808 REAL(KIND=dp), ALLOCATABLE :: Work(:) 809 REAL(KIND=dp), ALLOCATABLE :: x(:), y(:), z(:) 810 REAL(KIND=dp), ALLOCATABLE :: Pressure(:,:), FluxReal(:), FluxImag(:) 811 REAL(KIND=dp), ALLOCATABLE :: Bvector(:,:), BScalarReal(:), BscalarImag(:) 812 813 TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff 814 815 LOGICAL :: Dirichlet = .FALSE. 816 817! Initialize: 818! ----------- 819 Indicator = 0.0d0 820 Gnorm = 0.0d0 821 822 Metric = 0.0d0 823 DO i=1,3 824 Metric(i,i) = 1.0d0 825 END DO 826 827 SELECT CASE( CurrentCoordinateSystem() ) 828 CASE( AxisSymmetric, CylindricSymmetric ) 829 DIM = 3 830 CASE DEFAULT 831 DIM = CoordinateSystemDimension() 832 END SELECT 833! 834! --------------------------------------------- 835 Element => Edge % BoundaryInfo % Left 836 IF ( .NOT. ASSOCIATED( Element ) ) THEN 837 Element => Edge % BoundaryInfo % Right 838 ELSE IF ( ANY( Perm( Element % NodeIndexes ) <= 0 ) ) THEN 839 Element => Edge % BoundaryInfo % Right 840 END IF 841 842 IF ( .NOT. ASSOCIATED( Element ) ) RETURN 843 IF ( ANY( Perm( Element % NodeIndexes ) <= 0 ) ) RETURN 844 845 En = Edge % TYPE % NumberOfNodes 846 Pn = Element % TYPE % NumberOfNodes 847 848 ALLOCATE( EdgeNodes % x(En), EdgeNodes % y(En), EdgeNodes % z(En) ) 849 850 EdgeNodes % x = Mesh % Nodes % x(Edge % NodeIndexes) 851 EdgeNodes % y = Mesh % Nodes % y(Edge % NodeIndexes) 852 EdgeNodes % z = Mesh % Nodes % z(Edge % NodeIndexes) 853 854 ALLOCATE( Nodes % x(Pn), Nodes % y(Pn), Nodes % z(Pn) ) 855 856 Nodes % x = Mesh % Nodes % x(Element % NodeIndexes) 857 Nodes % y = Mesh % Nodes % y(Element % NodeIndexes) 858 Nodes % z = Mesh % Nodes % z(Element % NodeIndexes) 859 860 ALLOCATE( x(En), y(En), z(En), EdgeBasis(En), Basis(Pn), dBasisdx(Pn,3), & 861 Flux(en), Pressure(2,Pn), FluxReal(En), FluxImag(En), BVector(3,Pn), & 862 BScalarReal(Pn), BScalarImag(Pn) ) 863 864 DO l = 1,En 865 DO k = 1,Pn 866 IF ( Edge % NodeIndexes(l) == Element % NodeIndexes(k) ) THEN 867 x(l) = Element % TYPE % NodeU(k) 868 y(l) = Element % TYPE % NodeV(k) 869 z(l) = Element % TYPE % NodeW(k) 870 EXIT 871 END IF 872 END DO 873 END DO 874 875! 876! Integrate square of residual over boundary element: 877! --------------------------------------------------- 878 Indicator = 0.0d0 879 EdgeLength = 0.0d0 880 ResidualNorm = 0.0d0 881 Gnorm = 0.0d0 882 883 DO j=1,Model % NumberOfBCs 884 IF ( Edge % BoundaryInfo % Constraint /= Model % BCs(j) % Tag ) CYCLE 885! 886! ...given parameters: 887! -------------------- 888 889 CALL InputVector( Bvector, notScalar, 'Bvector', & 890 Model % BCs(j) % Values, Pn, Element % NodeIndexes ) 891 892 BscalarReal(1:Pn) = ListGetReal( Model % BCs(j) % Values, & 893 'Bscalar 1', Pn, Element % NodeIndexes, GotIt) 894 895 BscalarImag(1:Pn) = ListGetReal( Model % BCs(j) % Values, & 896 'Bscalar 2', Pn, Element % NodeIndexes, GotIt) 897 898 FluxReal(1:En) = ListGetReal( Model % BCs(j) % Values, & 899 'Wave Flux 1', En, Edge % NodeIndexes, gotIt ) 900 if( .not.GotIt ) FluxReal(1:En) = 0.0d0 901 902 FluxImag(1:En) = ListGetReal( Model % BCs(j) % Values, & 903 'Wave Flux 2', En, Edge % NodeIndexes, gotIt ) 904 if( .not.GotIt ) FluxImag(1:En) = 0.0d0 905! 906! get material parameters: 907! ------------------------ 908 k = ListGetInteger(Model % Bodies(Element % BodyId) % Values,'Material', & 909 minv=1, maxv=Model % NumberOfMaterials) 910! 911! elementwise nodal solution: 912! --------------------------- 913 Pressure(1,1:Pn) = Quant( 2*Perm(Element % NodeIndexes)-1 ) 914 Pressure(2,1:Pn) = Quant( 2*Perm(Element % NodeIndexes)-0 ) 915 916! do the integration: 917! ------------------- 918 EdgeLength = 0.0d0 919 ResidualNorm = 0.0d0 920 Gnorm = 0.0d0 921 922 IntegStuff = GaussPoints( Edge ) 923 924 DO t=1,IntegStuff % n 925 u = IntegStuff % u(t) 926 v = IntegStuff % v(t) 927 w = IntegStuff % w(t) 928 929 stat = ElementInfo( Edge, EdgeNodes, u, v, w, detJ, & 930 EdgeBasis, dBasisdx ) 931 932 IF ( CurrentCoordinateSystem() == Cartesian ) THEN 933 s = IntegStuff % s(t) * detJ 934 935 ELSE 936 u = SUM( EdgeBasis(1:En) * EdgeNodes % x(1:En) ) 937 v = SUM( EdgeBasis(1:En) * EdgeNodes % y(1:En) ) 938 w = SUM( EdgeBasis(1:En) * EdgeNodes % z(1:En) ) 939 940 CALL CoordinateSystemInfo( Metric, SqrtMetric, & 941 Symb, dSymb, u, v, w ) 942 943 s = IntegStuff % s(t) * detJ * SqrtMetric 944 945 END IF 946 947! 948! Integration point in parent element local 949! coordinates: 950! ----------------------------------------- 951 u = SUM( EdgeBasis(1:En) * x(1:En) ) 952 v = SUM( EdgeBasis(1:En) * y(1:En) ) 953 w = SUM( EdgeBasis(1:En) * z(1:En) ) 954 955 stat = ElementInfo( Element, Nodes, u, v, w, detJ, & 956 Basis, dBasisdx ) 957! 958! Flux and parameters at integration point: 959! ---------------------------------------- 960 961 Normal = NormalVector( Edge, EdgeNodes, u, v, .TRUE. ) 962 963 B1 = 0.0d0 964 B0r = 0.0d0 965 B0i = 0.0d0 966 967 B0r = SUM( BscalarReal(1:Pn) * Basis(1:Pn) ) 968 B0i = SUM( BscalarImag(1:Pn) * Basis(1:Pn) ) 969 do i = 1,dim 970 B1(i) = B1(i) + SUM( Bvector(i,1:Pn) * Basis(1:Pn) ) 971 end do 972 B1(1:dim) = B1(1:dim) + Normal(1:dim) 973 974 WaveFlux = 0.0d0 975 WaveFlux(1) = SUM( FluxReal(1:En) * EdgeBasis(1:En) ) 976 WaveFlux(2) = SUM( FluxImag(1:En) * EdgeBasis(1:En) ) 977 978 ResidualReal = -WaveFlux(1) 979 ResidualImag = -WaveFlux(2) 980 981 Greal = 0.0d0 982 Gimag = 0.0d0 983! 984! flux given by the computed solution, and 985! force norm for scaling the residual: 986! ----------------------------------------- 987 IF ( CurrentCoordinateSystem() == Cartesian ) THEN 988 DO k=1,DIM 989 990 ResidualReal = ResidualReal + & 991 SUM( dBasisdx(1:Pn,k) * Pressure(1,1:Pn) ) * B1(k) 992 993 ResidualImag = ResidualImag + & 994 SUM( dBasisdx(1:Pn,k) * Pressure(2,1:Pn) ) * B1(k) 995 996 GReal = GReal + & 997 SUM( dBasisdx(1:Pn,k) * Pressure(1,1:Pn) ) * B1(k) 998 999 Gimag = Gimag + & 1000 SUM( dBasisdx(1:Pn,k) * Pressure(2,1:Pn) ) * B1(k) 1001 1002 END DO 1003 1004 ResidualReal = ResidualReal + & 1005 SUM( Basis(1:Pn) * Pressure(1,1:Pn) ) * B0r & 1006 -SUM( Basis(1:Pn) * Pressure(2,1:Pn) ) * B0i 1007 1008 ResidualImag = ResidualImag + & 1009 SUM( Basis(1:Pn) * Pressure(1,1:Pn) ) * B0i & 1010 +SUM( Basis(1:Pn) * Pressure(2,1:Pn) ) * B0r 1011 1012 Greal = Greal + & 1013 SUM( Basis(1:Pn) * Pressure(1,1:Pn) ) * B0r & 1014 -SUM( Basis(1:Pn) * Pressure(2,1:Pn) ) * B0i 1015 1016 Gimag = Gimag + & 1017 SUM( Basis(1:Pn) * Pressure(1,1:Pn) ) * B0i & 1018 +SUM( Basis(1:Pn) * Pressure(2,1:Pn) ) * B0r 1019 1020 ELSE 1021! DO k=1,DIM 1022! DO l=1,DIM 1023! Residual = Residual + Metric(k,l) * Conductivity * & 1024! SUM( dBasisdx(1:Pn,k) * Temperature(1:Pn) ) * Normal(l) 1025! 1026! Gnorm = Gnorm + s * (Metric(k,l) * Conductivity * & 1027! SUM(dBasisdx(1:Pn,k) * Temperature(1:Pn) ) * Normal(l))**2 1028! END DO 1029! END DO 1030 END IF 1031 1032 EdgeLength = EdgeLength + s 1033 Gnorm = Gnorm + s * ( Greal **2 + Gimag**2 ) 1034 1035 IF ( .NOT. Dirichlet ) THEN 1036 ResidualNorm = ResidualNorm + s * (ResidualReal**2 + ResidualImag**2) 1037 END IF 1038 1039 END DO 1040 1041 EXIT 1042 1043 END DO 1044 1045 IF ( CoordinateSystemDimension() == 3 ) THEN 1046 EdgeLength = SQRT(EdgeLength) 1047 END IF 1048 1049 Gnorm = EdgeLength * Gnorm 1050 1051 Indicator = EdgeLength * ResidualNorm 1052 1053 DEALLOCATE( Nodes % x, Nodes % y, Nodes % z) 1054 DEALLOCATE( EdgeNodes % x, EdgeNodes % y, EdgeNodes % z) 1055 1056 DEALLOCATE( x, y, z, EdgeBasis, Basis, dBasisdx, & 1057 Flux, Pressure, FluxReal, FluxImag, BVector, & 1058 BScalarReal, BScalarImag ) 1059 1060 1061 1062!------------------------------------------------------------------------------ 1063 1064contains 1065 1066!------------------------------------------------------------------------------ 1067 SUBROUTINE InputVector( Tensor, IsScalar, Name, Material, n, NodeIndexes ) 1068!------------------------------------------------------------------------------ 1069 REAL(KIND=dp) :: Tensor(:,:) 1070 INTEGER :: n, NodeIndexes(:) 1071 LOGICAL :: IsScalar 1072 CHARACTER(LEN=*) :: Name 1073 TYPE(ValueList_t), POINTER :: Material 1074!------------------------------------------------------------------------------ 1075 LOGICAL :: FirstTime = .TRUE., stat 1076 REAL(KIND=dp), POINTER :: Hwrk(:,:,:) 1077 1078 SAVE FirstTime, Hwrk 1079!------------------------------------------------------------------------------ 1080 IF ( FirstTime ) THEN 1081 NULLIFY( Hwrk ) 1082 FirstTime = .FALSE. 1083 END IF 1084 1085 Tensor = 0.0d0 1086 1087 CALL ListGetRealArray( Material, Name, Hwrk, n, NodeIndexes, stat ) 1088 IsScalar = SIZE(HWrk,1) == 1 .AND. SIZE(HWrk,2) == 1 1089 1090 IF ( .NOT. stat ) RETURN 1091 1092 IF ( SIZE(Hwrk,1) == 1 ) THEN 1093 1094 DO i=1,MIN(3,SIZE(Hwrk,2)) 1095 Tensor( i,1:n ) = Hwrk( 1,1,1:n ) 1096 END DO 1097 1098 ELSE 1099 1100 DO i=1,MIN(3,SIZE(Hwrk,1)) 1101 Tensor( i,1:n ) = Hwrk( i,1,1:n ) 1102 END DO 1103 1104 END IF 1105!------------------------------------------------------------------------------ 1106 END SUBROUTINE InputVector 1107!------------------------------------------------------------------------------ 1108 END FUNCTION DCRBoundaryResidual 1109!------------------------------------------------------------------------------ 1110 1111 1112!------------------------------------------------------------------------------ 1113 FUNCTION DCREdgeResidual( Model, Edge, Mesh, Quant, Perm ) RESULT( Indicator ) 1114!------------------------------------------------------------------------------ 1115 USE DefUtils 1116 IMPLICIT NONE 1117 1118 TYPE(Model_t) :: Model 1119 INTEGER :: Perm(:) 1120 REAL(KIND=dp) :: Quant(:), Indicator(2) 1121 TYPE( Mesh_t ), POINTER :: Mesh 1122 TYPE( Element_t ), POINTER :: Edge 1123!------------------------------------------------------------------------------ 1124 1125 TYPE(Nodes_t) :: Nodes, EdgeNodes 1126 TYPE(Element_t), POINTER :: Element, Bndry 1127 1128 INTEGER :: i,j,k,l,n,t,DIM,En,Pn 1129 LOGICAL :: stat, GotIt 1130! REAL(KIND=dp), POINTER :: Hwrk(:,:,:) 1131 1132 REAL(KIND=dp) :: Grad(3,3), Normal(3), EdgeLength, Jump, JumpReal, JumpImag, & 1133 GradReal(3,3),GradImag(3,3) 1134 REAL(KIND=dp) :: u, v, w, s, detJ 1135 REAL(KIND=dp) :: Residual, ResidualNorm, Area 1136 REAL(kind=dp) :: A2(3,3), A1r(3), A1i(3), A0r, A0i 1137 REAL(KIND=dp) :: SqrtMetric, Metric(3,3), Symb(3,3,3), dSymb(3,3,3,3) 1138 1139 REAL(KIND=dp), ALLOCATABLE :: AvectorImag(:,:), AscalarReal(:), AscalarImag(:) 1140 REAL(KIND=dp), ALLOCATABLE :: Flux(:), x(:), y(:), z(:) 1141 REAL(KIND=dp), ALLOCATABLE :: Amatrix(:,:,:) 1142 REAL(KIND=dp), ALLOCATABLE :: EdgeBasis(:) 1143 REAL(KIND=dp), ALLOCATABLE :: AvectorReal(:,:) 1144 REAL(KIND=dp), ALLOCATABLE :: Basis(:), dBasisdx(:,:) 1145 REAL(KIND=dp), ALLOCATABLE :: Temperature(:), Pressure(:,:) 1146 1147 LOGICAL :: notScalar = .TRUE. 1148 1149 TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff 1150 1151! LOGICAL :: First = .TRUE. 1152! SAVE Hwrk, First 1153!------------------------------------------------------------------------------ 1154 1155! Initialize: 1156! ----------- 1157 1158! IF ( First ) THEN 1159! First = .FALSE. 1160! NULLIFY( Hwrk ) 1161! END IF 1162 1163 SELECT CASE( CurrentCoordinateSystem() ) 1164 CASE( AxisSymmetric, CylindricSymmetric ) 1165 DIM = 3 1166 CASE DEFAULT 1167 DIM = CoordinateSystemDimension() 1168 END SELECT 1169 1170 Metric = 0.0d0 1171 DO i = 1,3 1172 Metric(i,i) = 1.0d0 1173 END DO 1174 1175 Grad = 0.0d0 1176 GradReal = 0.0d0 1177 GradImag = 0.0d0 1178! 1179! --------------------------------------------- 1180 1181 Element => Edge % BoundaryInfo % Left 1182 n = Element % TYPE % NumberOfNodes 1183 1184 Element => Edge % BoundaryInfo % Right 1185 n = MAX( n, Element % TYPE % NumberOfNodes ) 1186 1187 ALLOCATE( Nodes % x(n), Nodes % y(n), Nodes % z(n) ) 1188 1189 En = Edge % TYPE % NumberOfNodes 1190 ALLOCATE( EdgeNodes % x(En), EdgeNodes % y(En), EdgeNodes % z(En) ) 1191 1192 EdgeNodes % x = Mesh % Nodes % x(Edge % NodeIndexes) 1193 EdgeNodes % y = Mesh % Nodes % y(Edge % NodeIndexes) 1194 EdgeNodes % z = Mesh % Nodes % z(Edge % NodeIndexes) 1195 1196 ALLOCATE( AvectorImag(3,n), AscalarReal(n), AscalarImag(n), Flux(en), & 1197 x(En), y(En), z(En), AMatrix(3,3,n), EdgeBasis(En), AvectorReal(3,n), & 1198 Basis(n), dBasisdx(n,3), Temperature(n), Pressure(2,n) ) 1199 1200! Integrate square of jump over edge: 1201! ----------------------------------- 1202 ResidualNorm = 0.0d0 1203 EdgeLength = 0.0d0 1204 Indicator = 0.0d0 1205 1206 IntegStuff = GaussPoints( Edge ) 1207 1208 DO t=1,IntegStuff % n 1209 1210 u = IntegStuff % u(t) 1211 v = IntegStuff % v(t) 1212 w = IntegStuff % w(t) 1213 1214 stat = ElementInfo( Edge, EdgeNodes, u, v, w, detJ, & 1215 EdgeBasis, dBasisdx ) 1216 1217 Normal = NormalVector( Edge, EdgeNodes, u, v, .FALSE. ) 1218 1219 IF ( CurrentCoordinateSystem() == Cartesian ) THEN 1220 s = IntegStuff % s(t) * detJ 1221 ELSE 1222 u = SUM( EdgeBasis(1:En) * EdgeNodes % x(1:En) ) 1223 v = SUM( EdgeBasis(1:En) * EdgeNodes % y(1:En) ) 1224 w = SUM( EdgeBasis(1:En) * EdgeNodes % z(1:En) ) 1225 1226 CALL CoordinateSystemInfo( Metric, SqrtMetric, & 1227 Symb, dSymb, u, v, w ) 1228 s = IntegStuff % s(t) * detJ * SqrtMetric 1229 END IF 1230! 1231! Compute flux over the edge as seen by elements 1232! on both sides of the edge: 1233! ---------------------------------------------- 1234 DO i = 1,2 1235 SELECT CASE(i) 1236 CASE(1) 1237 Element => Edge % BoundaryInfo % Left 1238 CASE(2) 1239 Element => Edge % BoundaryInfo % Right 1240 END SELECT 1241! 1242! Can this really happen (maybe it can...) ? 1243! ------------------------------------------- 1244 IF ( ANY( Perm( Element % NodeIndexes ) <= 0 ) ) CYCLE 1245! 1246! Next, get the integration point in parent 1247! local coordinates: 1248! ----------------------------------------- 1249 Pn = Element % TYPE % NumberOfNodes 1250 1251 DO j = 1,En 1252 DO k = 1,Pn 1253 IF ( Edge % NodeIndexes(j) == Element % NodeIndexes(k) ) THEN 1254 x(j) = Element % TYPE % NodeU(k) 1255 y(j) = Element % TYPE % NodeV(k) 1256 z(j) = Element % TYPE % NodeW(k) 1257 EXIT 1258 END IF 1259 END DO 1260 END DO 1261 1262 u = SUM( EdgeBasis(1:En) * x(1:En) ) 1263 v = SUM( EdgeBasis(1:En) * y(1:En) ) 1264 w = SUM( EdgeBasis(1:En) * z(1:En) ) 1265! 1266! Get parent element basis & derivatives at the 1267! integration point: 1268! --------------------------------------------- 1269 Nodes % x(1:Pn) = Mesh % Nodes % x(Element % NodeIndexes) 1270 Nodes % y(1:Pn) = Mesh % Nodes % y(Element % NodeIndexes) 1271 Nodes % z(1:Pn) = Mesh % Nodes % z(Element % NodeIndexes) 1272 1273 stat = ElementInfo( Element, Nodes, u, v, w, detJ, & 1274 Basis, dBasisdx ) 1275! 1276! Material parameters: 1277! -------------------- 1278 k = ListGetInteger( Model % Bodies( & 1279 Element % BodyId) % Values, 'Material', minv=1,maxv=Model % NumberOfMaterials ) 1280 1281 CALL InputTensor( Amatrix, notScalar, & 1282 'Amatrix', Model % Materials(k) % Values, Pn, Element % NodeIndexes ) 1283 1284 CALL InputVector( AvectorReal, notScalar, & 1285 'Avector 1', Model % Materials(k) % Values, Pn, Element % NodeIndexes ) 1286 1287 CALL InputVector( AvectorImag, notScalar, & 1288 'Avector 2', Model % Materials(k) % Values, Pn, Element % NodeIndexes ) 1289 1290 AscalarReal(1:Pn) = ListGetReal( Model % Materials(k) % Values, & 1291 'Ascalar 1', Pn, Element % NodeIndexes, GotIt) 1292 1293 AscalarImag(1:Pn) = ListGetReal( Model % Materials(k) % Values, & 1294 'Ascalar 2', Pn, Element % NodeIndexes, GotIt) 1295 1296 A2 = 0.0d0 1297 A1r = 0.0d0 1298 A1i = 0.0d0 1299 A0r = 0.0d0 1300 A0i = 0.0d0 1301 1302 A0r = SUM( AscalarReal(1:Pn) * Basis(1:Pn) ) 1303 A0i = SUM( AscalarImag(1:Pn) * Basis(1:Pn) ) 1304 do j = 1,dim 1305 A1r(j) = A1r(j) + SUM( AvectorReal(j,1:Pn) * Basis(1:Pn) ) 1306 A1i(j) = A1i(j) + SUM( AvectorImag(j,1:Pn) * Basis(1:Pn) ) 1307 do k = 1,dim 1308 A2(j,k) = A2(j,k) + SUM( Amatrix(j,k,1:Pn) * Basis(1:Pn) ) 1309 end do 1310 end do 1311! 1312! Pressure at element nodal points 1313! --------------------------------- 1314 do k = 1,2 1315 Pressure(k,1:Pn) = Quant( 2*Perm(Element % NodeIndexes)-2+k ) 1316 end do 1317! 1318! Finally, the flux: 1319! ------------------ 1320 DO j=1,DIM 1321 do k = 1,dim 1322 GradReal(j,i) = GradReal(j,i) & 1323 + SUM( dBasisdx(1:Pn,k) * Pressure(1,1:Pn) )*A2(j,k) 1324 GradImag(j,i) = GradImag(j,i) & 1325 + SUM( dBasisdx(1:Pn,k) * Pressure(2,1:Pn) )*A2(j,k) 1326 end do 1327 END DO 1328 1329 END DO 1330 1331! Compute squre of the flux jump: 1332! ------------------------------- 1333 EdgeLength = EdgeLength + s 1334 JumpReal = 0.0d0 1335 JumpImag = 0.0d0 1336 1337 DO k=1,DIM 1338 IF ( CurrentCoordinateSystem() == Cartesian ) THEN 1339 JumpReal = JumpReal + (GradReal(k,1) - GradReal(k,2)) * Normal(k) 1340 JumpImag = JumpImag + (GradImag(k,1) - GradImag(k,2)) * Normal(k) 1341 ELSE 1342! DO l=1,DIM 1343! Jump = Jump + & 1344! Metric(k,l) * (Grad(k,1) - Grad(k,2)) * Normal(l) 1345! END DO 1346 END IF 1347 END DO 1348 ResidualNorm = ResidualNorm + s * ( JumpReal**2 + JumpImag**2 ) 1349 END DO 1350 1351 IF ( CoordinateSystemDimension() == 3 ) THEN 1352 EdgeLength = SQRT(EdgeLength) 1353 END IF 1354 Indicator = EdgeLength * ResidualNorm 1355 1356 DEALLOCATE( Nodes % x, Nodes % y, Nodes % z) 1357 DEALLOCATE( EdgeNodes % x, EdgeNodes % y, EdgeNodes % z) 1358 1359 DEALLOCATE( AvectorImag, AscalarReal, AscalarImag, Flux, & 1360 x, y, z, AMatrix, EdgeBasis, AvectorReal, Basis, dBasisdx,& 1361 Temperature, Pressure ) 1362!------------------------------------------------------------------------------ 1363 1364contains 1365 1366!------------------------------------------------------------------------------ 1367 SUBROUTINE InputTensor( Tensor, IsScalar, Name, Material, n, NodeIndexes ) 1368!------------------------------------------------------------------------------ 1369 REAL(KIND=dp) :: Tensor(:,:,:) 1370 INTEGER :: n, NodeIndexes(:) 1371 LOGICAL :: IsScalar 1372 CHARACTER(LEN=*) :: Name 1373 TYPE(ValueList_t), POINTER :: Material 1374!------------------------------------------------------------------------------ 1375 LOGICAL :: FirstTime = .TRUE., stat 1376 REAL(KIND=dp), POINTER :: Hwrk(:,:,:) 1377 1378 SAVE FirstTime, Hwrk 1379!------------------------------------------------------------------------------ 1380 IF ( FirstTime ) THEN 1381 NULLIFY( Hwrk ) 1382 FirstTime = .FALSE. 1383 END IF 1384 1385 Tensor = 0.0d0 1386 1387 CALL ListGetRealArray( Material, Name, Hwrk, n, NodeIndexes, stat ) 1388 IsScalar = SIZE(HWrk,1) == 1 .AND. SIZE(HWrk,2) == 1 1389 1390 IF ( .NOT. stat ) RETURN 1391 1392 IF ( SIZE(Hwrk,1) == 1 ) THEN 1393 1394 DO i=1,MIN(3,SIZE(Hwrk,2)) 1395 Tensor( i,i,1:n ) = Hwrk( 1,1,1:n ) 1396 END DO 1397 1398 ELSE IF ( SIZE(Hwrk,2) == 1 ) THEN 1399 1400 DO i=1,MIN(3,SIZE(Hwrk,1)) 1401 Tensor(i,i,1:n) = Hwrk(i,1,1:n) 1402 END DO 1403 1404 ELSE 1405 1406 DO i=1,MIN(3,SIZE(Hwrk,1)) 1407 DO j=1,MIN(3,SIZE(Hwrk,2)) 1408 Tensor( i,j,1:n ) = Hwrk(i,j,1:n) 1409 END DO 1410 END DO 1411 1412 END IF 1413!------------------------------------------------------------------------------ 1414 END SUBROUTINE InputTensor 1415!------------------------------------------------------------------------------ 1416 1417 1418!------------------------------------------------------------------------------ 1419 SUBROUTINE InputVector( Tensor, IsScalar, Name, Material, n, NodeIndexes ) 1420!------------------------------------------------------------------------------ 1421 REAL(KIND=dp) :: Tensor(:,:) 1422 INTEGER :: n, NodeIndexes(:) 1423 LOGICAL :: IsScalar 1424 CHARACTER(LEN=*) :: Name 1425 TYPE(ValueList_t), POINTER :: Material 1426!------------------------------------------------------------------------------ 1427 LOGICAL :: FirstTime = .TRUE., stat 1428 REAL(KIND=dp), POINTER :: Hwrk(:,:,:) 1429 1430 SAVE FirstTime, Hwrk 1431!------------------------------------------------------------------------------ 1432 IF ( FirstTime ) THEN 1433 NULLIFY( Hwrk ) 1434 FirstTime = .FALSE. 1435 END IF 1436 1437 Tensor = 0.0d0 1438 1439 CALL ListGetRealArray( Material, Name, Hwrk, n, NodeIndexes, stat ) 1440 IsScalar = SIZE(HWrk,1) == 1 .AND. SIZE(HWrk,2) == 1 1441 1442 IF ( .NOT. stat ) RETURN 1443 1444 IF ( SIZE(Hwrk,1) == 1 ) THEN 1445 1446 DO i=1,MIN(3,SIZE(Hwrk,2)) 1447 Tensor( i,1:n ) = Hwrk( 1,1,1:n ) 1448 END DO 1449 1450 ELSE 1451 1452 DO i=1,MIN(3,SIZE(Hwrk,1)) 1453 Tensor( i,1:n ) = Hwrk( i,1,1:n ) 1454 END DO 1455 1456 END IF 1457!------------------------------------------------------------------------------ 1458 END SUBROUTINE InputVector 1459!------------------------------------------------------------------------------ 1460 1461 1462 END FUNCTION DCREdgeResidual 1463!------------------------------------------------------------------------------ 1464 1465 1466!------------------------------------------------------------------------------ 1467 FUNCTION DCRInsideResidual( Model, Element, Mesh, & 1468 Quant, Perm, Fnorm ) RESULT( Indicator ) 1469!------------------------------------------------------------------------------ 1470 USE DefUtils 1471!------------------------------------------------------------------------------ 1472 IMPLICIT NONE 1473!------------------------------------------------------------------------------ 1474 TYPE(Model_t) :: Model 1475 INTEGER :: Perm(:) 1476 REAL(KIND=dp) :: Quant(:), Indicator(2), Fnorm 1477 TYPE( Mesh_t ), POINTER :: Mesh 1478 TYPE( Element_t ), POINTER :: Element 1479!------------------------------------------------------------------------------ 1480 1481 TYPE(Nodes_t) :: Nodes 1482 1483 INTEGER :: i,j,k,l,n,t,DIM 1484 LOGICAL :: stat, GotIt 1485 TYPE( Variable_t ), POINTER :: Var 1486 1487 REAL(KIND=dp) :: SqrtMetric, Metric(3,3), Symb(3,3,3), dSymb(3,3,3,3) 1488 REAL(KIND=dp) :: Source, Residual, ResidualNorm, Area 1489 REAL(kind=dp) :: A2(3,3), A1r(3), A1i(3), A0r, A0i 1490 REAL(KIND=dp) :: u, v, w, s, detJ, ResidualReal, ResidualImag 1491 1492 REAL(KIND=dp), ALLOCATABLE :: NodalSource(:,:), Basis(:) 1493 REAL(KIND=dp), ALLOCATABLE :: dBasisdx(:,:), ddBasisddx(:,:,:), Pressione(:,:) 1494 REAL(KIND=dp), ALLOCATABLE :: Amatrix(:,:,:), AvectorReal(:,:) 1495 REAL(KIND=dp), ALLOCATABLE :: AvectorImag(:,:), AscalarReal(:), AscalarImag(:) 1496 1497 LOGICAL :: notScalar = .TRUE. 1498 TYPE( ValueList_t ), POINTER :: Material 1499 TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff 1500 1501! LOGICAL :: First = .TRUE. 1502! SAVE Hwrk, First 1503!------------------------------------------------------------------------------ 1504 1505! Initialize: 1506! ----------- 1507 Indicator = 0.0d0 1508 Fnorm = 0.0d0 1509! 1510! Check if this eq. computed in this element: 1511! ------------------------------------------- 1512 IF ( ANY( Perm( Element % NodeIndexes ) <= 0 ) ) RETURN 1513 1514 Metric = 0.0d0 1515 DO i=1,3 1516 Metric(i,i) = 1.0d0 1517 END DO 1518 1519 SELECT CASE( CurrentCoordinateSystem() ) 1520 CASE( AxisSymmetric, CylindricSymmetric ) 1521 DIM = 3 1522 CASE DEFAULT 1523 DIM = CoordinateSystemDimension() 1524 END SELECT 1525! 1526! Element nodal points: 1527! --------------------- 1528 n = Element % TYPE % NumberOfNodes 1529 1530 ALLOCATE( Nodes % x(n), Nodes % y(n), Nodes % z(n) ) 1531 Nodes % x = Mesh % Nodes % x(Element % NodeIndexes) 1532 Nodes % y = Mesh % Nodes % y(Element % NodeIndexes) 1533 Nodes % z = Mesh % Nodes % z(Element % NodeIndexes) 1534 1535 ALLOCATE( NodalSource(2,n), Basis(n), dBasisdx(n,3), ddBasisddx(n,3,3), & 1536 Pressione(2,n), AMatrix(3,3,n), AvectorReal(3,n), AvectorImag(3,n), & 1537 AscalarReal(n), AscalarImag(n) ) 1538! 1539! Elementwise nodal solution: 1540! --------------------------- 1541 Pressione(1,1:n) = Quant( 2*Perm(Element % NodeIndexes)-1 ) 1542 Pressione(2,1:n) = Quant( 2*Perm(Element % NodeIndexes)-0 ) 1543! 1544! Material parameters: sound speed: 1545! --------------------------------- 1546 k = ListGetInteger( Model % Bodies(Element % BodyId) % Values, 'Material', & 1547 minv=1, maxv=Model % NumberOfMaterials ) 1548 1549 Material => Model % Materials(k) % Values 1550 1551 CALL InputTensor( Amatrix, notScalar, & 1552 'Amatrix', Material, n, Element % NodeIndexes ) 1553 1554 CALL InputVector( AvectorReal, notScalar, & 1555 'Avector 1', Material, n, Element % NodeIndexes ) 1556 1557 CALL InputVector( AvectorImag, notScalar, & 1558 'Avector 2', Material, n, Element % NodeIndexes ) 1559 1560 AscalarReal(1:n) = ListGetReal( Material, & 1561 'Ascalar 1', n, Element % NodeIndexes, GotIt) 1562 1563 AscalarImag(1:n) = ListGetReal( Material, & 1564 'Ascalar 2', n, Element % NodeIndexes, GotIt) 1565! 1566! Source term: 1567! ------------ 1568 k = ListGetInteger( & 1569 Model % Bodies(Element % BodyId) % Values,'Body Force',GotIt, & 1570 1, Model % NumberOFBodyForces ) 1571 1572 NodalSource = 0.0d0 1573 IF ( GotIt .AND. k > 0 ) THEN 1574 NodalSource(1,1:n) = ListGetReal( Model % BodyForces(k) % Values, & 1575 'Pressure Source 1', n, Element % NodeIndexes, GotIt ) 1576 1577 NodalSource(2,1:n) = ListGetReal( Model % BodyForces(k) % Values, & 1578 'Pressure Source 2', n, Element % NodeIndexes, GotIt ) 1579 END IF 1580! 1581! Integrate square of residual over element: 1582! ------------------------------------------ 1583 1584 ResidualNorm = 0.0d0 1585 Area = 0.0d0 1586 1587 IntegStuff = GaussPoints( Element ) 1588 1589 DO t=1,IntegStuff % n 1590 u = IntegStuff % u(t) 1591 v = IntegStuff % v(t) 1592 w = IntegStuff % w(t) 1593 1594 stat = ElementInfo( Element, Nodes, u, v, w, detJ, & 1595 Basis, dBasisdx, ddBasisddx, .TRUE. ) 1596 1597 IF ( CurrentCoordinateSystem() == Cartesian ) THEN 1598 s = IntegStuff % s(t) * detJ 1599 1600 ELSE 1601 u = SUM( Basis(1:n) * Nodes % x(1:n) ) 1602 v = SUM( Basis(1:n) * Nodes % y(1:n) ) 1603 w = SUM( Basis(1:n) * Nodes % z(1:n) ) 1604 1605 CALL CoordinateSystemInfo( Metric, SqrtMetric, & 1606 Symb, dSymb, u, v, w ) 1607 s = IntegStuff % s(t) * detJ * SqrtMetric 1608 1609 END IF 1610 1611 ResidualReal = 0.0d0 1612 ResidualImag = 0.0d0 1613 1614 A2 = 0.0d0 1615 A1r = 0.0d0 1616 A1i = 0.0d0 1617 A0r = 0.0d0 1618 A0i = 0.0d0 1619 1620 A0r = SUM( AscalarReal(1:n) * Basis(1:n) ) 1621 A0i = SUM( AscalarImag(1:n) * Basis(1:n) ) 1622 do i = 1,dim 1623 A1r(i) = A1r(i) + SUM( AvectorReal(i,1:n) * Basis(1:n) ) 1624 A1i(i) = A1i(i) + SUM( AvectorImag(i,1:n) * Basis(1:n) ) 1625 do j = 1,dim 1626 A2(i,j) = A2(i,j) + SUM( Amatrix(i,j,1:n) * Basis(1:n) ) 1627 end do 1628 end do 1629 1630 IF ( CurrentCoordinateSystem() == Cartesian ) THEN 1631 1632! diffusion 1633!----------- 1634 do i = 1,dim 1635 do j = 1,dim 1636 ResidualReal = ResidualReal & 1637 - SUM( Pressione(1,1:n) * ddBasisddx(1:n,i,j) ) * A2(i,j) 1638 1639 ResidualImag = ResidualImag & 1640 - SUM( Pressione(2,1:n) * ddBasisddx(1:n,i,j) ) * A2(i,j) 1641 end do 1642 end do 1643 1644! convection 1645!------------ 1646 do i = 1,dim 1647 ResidualReal = ResidualReal + & 1648 SUM( dBasisdx(1:n,i) * Pressione(1,1:n) ) * A1r(i) & 1649 -SUM( dBasisdx(1:n,i) * Pressione(2,1:n) ) * A1i(i) 1650 1651 ResidualImag = ResidualImag + & 1652 SUM( dBasisdx(1:n,i) * Pressione(1,1:n) ) * A1i(i) & 1653 +SUM( dBasisdx(1:n,i) * Pressione(2,1:n) ) * A1r(i) 1654 end do 1655 1656! reaction 1657!---------- 1658 ResidualReal = ResidualReal + & 1659 SUM( Basis(1:n) * Pressione(1,1:n) ) * A0r & 1660 -SUM( Basis(1:n) * Pressione(2,1:n) ) * A0i 1661 1662 ResidualImag = ResidualImag + & 1663 SUM( Basis(1:n) * Pressione(1,1:n) ) * A0i & 1664 +SUM( Basis(1:n) * Pressione(2,1:n) ) * A0r 1665 1666 ELSE 1667 print *,'Only cartesian coordinates implemented at the moment!' 1668 stop 1669 1670! DO j=1,DIM 1671! DO k=1,DIM 1672! 1673! - g^{jk} C_{,k}T_{j}: 1674! --------------------- 1675! 1676! Residual = Residual - Metric(j,k) * & 1677! SUM( Temperature(1:n) * dBasisdx(1:n,j) ) * & 1678! SUM( NodalConductivity(1:n) * dBasisdx(1:n,k) ) 1679 1680! 1681! - g^{jk} C T_{,jk}: 1682! ------------------- 1683! 1684! Residual = Residual - Metric(j,k) * Conductivity * & 1685! SUM( Temperature(1:n) * ddBasisddx(1:n,j,k) ) 1686! 1687! + g^{jk} C {_jk^l} T_{,l}: 1688! --------------------------- 1689! DO l=1,DIM 1690! Residual = Residual + Metric(j,k) * Conductivity * & 1691! Symb(j,k,l) * SUM( Temperature(1:n) * dBasisdx(1:n,l) ) 1692! END DO 1693! END DO 1694! END DO 1695 END IF 1696 1697! 1698! Compute also force norm for scaling the residual: 1699! ------------------------------------------------- 1700 DO i=1,DIM 1701 Fnorm = Fnorm & 1702 + s * ( SUM( NodalSource(1,1:n) * Basis(1:n) ) ) ** 2 & 1703 + s * ( SUM( NodalSource(2,1:n) * Basis(1:n) ) ) ** 2 1704 END DO 1705 1706 Area = Area + s 1707 ResidualNorm = ResidualNorm + s * ( ResidualReal ** 2 + ResidualImag ** 2 ) 1708 END DO 1709 1710 Fnorm = Element % hk**2 * Fnorm 1711 Indicator = Element % hK**2 * ResidualNorm 1712 1713 DEALLOCATE( Nodes % x, Nodes % y, Nodes % z ) 1714 DEALLOCATE( NodalSource, Basis, dBasisdx, ddBasisddx, Pressione, & 1715 AMatrix, AvectorReal, AvectorImag,AscalarReal, AscalarImag ) 1716 1717CONTAINS 1718 1719!------------------------------------------------------------------------------ 1720 SUBROUTINE InputTensor( Tensor, IsScalar, Name, Material, n, NodeIndexes ) 1721!------------------------------------------------------------------------------ 1722 REAL(KIND=dp) :: Tensor(:,:,:) 1723 INTEGER :: n, NodeIndexes(:) 1724 LOGICAL :: IsScalar 1725 CHARACTER(LEN=*) :: Name 1726 TYPE(ValueList_t), POINTER :: Material 1727!------------------------------------------------------------------------------ 1728 LOGICAL :: FirstTime = .TRUE., stat 1729 REAL(KIND=dp), POINTER :: Hwrk(:,:,:) 1730 1731 SAVE FirstTime, Hwrk 1732!------------------------------------------------------------------------------ 1733 IF ( FirstTime ) THEN 1734 NULLIFY( Hwrk ) 1735 FirstTime = .FALSE. 1736 END IF 1737 1738 Tensor = 0.0d0 1739 1740 CALL ListGetRealArray( Material, Name, Hwrk, n, NodeIndexes, stat ) 1741 IsScalar = SIZE(HWrk,1) == 1 .AND. SIZE(HWrk,2) == 1 1742 1743 IF ( .NOT. stat ) RETURN 1744 1745 IF ( SIZE(Hwrk,1) == 1 ) THEN 1746 1747 DO i=1,MIN(3,SIZE(Hwrk,2)) 1748 Tensor( i,i,1:n ) = Hwrk( 1,1,1:n ) 1749 END DO 1750 1751 ELSE IF ( SIZE(Hwrk,2) == 1 ) THEN 1752 1753 DO i=1,MIN(3,SIZE(Hwrk,1)) 1754 Tensor(i,i,1:n) = Hwrk(i,1,1:n) 1755 END DO 1756 1757 ELSE 1758 1759 DO i=1,MIN(3,SIZE(Hwrk,1)) 1760 DO j=1,MIN(3,SIZE(Hwrk,2)) 1761 Tensor( i,j,1:n ) = Hwrk(i,j,1:n) 1762 END DO 1763 END DO 1764 1765 END IF 1766!------------------------------------------------------------------------------ 1767 END SUBROUTINE InputTensor 1768!------------------------------------------------------------------------------ 1769 1770 1771!------------------------------------------------------------------------------ 1772 SUBROUTINE InputVector( Tensor, IsScalar, Name, Material, n, NodeIndexes ) 1773!------------------------------------------------------------------------------ 1774 REAL(KIND=dp) :: Tensor(:,:) 1775 INTEGER :: n, NodeIndexes(:) 1776 LOGICAL :: IsScalar 1777 CHARACTER(LEN=*) :: Name 1778 TYPE(ValueList_t), POINTER :: Material 1779!------------------------------------------------------------------------------ 1780 LOGICAL :: FirstTime = .TRUE., stat 1781 REAL(KIND=dp), POINTER :: Hwrk(:,:,:) 1782 1783 SAVE FirstTime, Hwrk 1784!------------------------------------------------------------------------------ 1785 IF ( FirstTime ) THEN 1786 NULLIFY( Hwrk ) 1787 FirstTime = .FALSE. 1788 END IF 1789 1790 Tensor = 0.0d0 1791 1792 CALL ListGetRealArray( Material, Name, Hwrk, n, NodeIndexes, stat ) 1793 IsScalar = SIZE(HWrk,1) == 1 .AND. SIZE(HWrk,2) == 1 1794 1795 IF ( .NOT. stat ) RETURN 1796 1797 IF ( SIZE(Hwrk,1) == 1 ) THEN 1798 1799 DO i=1,MIN(3,SIZE(Hwrk,2)) 1800 Tensor( i,1:n ) = Hwrk( 1,1,1:n ) 1801 END DO 1802 1803 ELSE 1804 1805 DO i=1,MIN(3,SIZE(Hwrk,1)) 1806 Tensor( i,1:n ) = Hwrk( i,1,1:n ) 1807 END DO 1808 1809 END IF 1810!------------------------------------------------------------------------------ 1811 END SUBROUTINE InputVector 1812!------------------------------------------------------------------------------ 1813 1814!------------------------------------------------------------------------------ 1815 END FUNCTION DCRInsideResidual 1816!------------------------------------------------------------------------------ 1817 1818 1819