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! * Subroutine for saving scalar data to files 27! * 28! ****************************************************************************** 29! * 30! * Authors: Peter Råback 31! * Email: Peter.Raback@csc.fi 32! * Web: http://www.csc.fi/elmer 33! * Address: CSC - IT Center for Science Ltd. 34! * Keilaranta 14 35! * 02101 Espoo, Finland 36! * 37! * Original Date: 20 Nov 2001 38! * 39! *****************************************************************************/ 40 41!> \ingroup Solvers 42!> \{ 43 44SUBROUTINE SaveScalars_init( Model,Solver,dt,TransientSimulation ) 45!------------------------------------------------------------------------------ 46 USE DefUtils 47 48 IMPLICIT NONE 49!------------------------------------------------------------------------------ 50 TYPE(Solver_t), TARGET :: Solver 51 TYPE(Model_t) :: Model 52 REAL(KIND=dp) :: dt 53 LOGICAL :: TransientSimulation 54!------------------------------------------------------------------------------ 55! Local variables 56!------------------------------------------------------------------------------ 57 INTEGER :: NormInd, LineInd, MarkerUnit, i 58 LOGICAL :: GotIt, MarkFailed, AvoidFailed 59 CHARACTER(LEN=MAX_NAME_LEN) :: Name 60 61 62 ! If we want to show a pseudonorm add a variable for which the norm 63 ! is associated with. 64 NormInd = ListGetInteger( Solver % Values,'Show Norm Index',GotIt) 65 IF( NormInd > 0 ) THEN 66 Name = ListGetString( Solver % Values, 'Equation',GotIt) 67 IF( .NOT. ListCheckPresent( Solver % Values,'Variable') ) THEN 68 CALL ListAddString( Solver % Values,'Variable',& 69 '-nooutput -global '//TRIM(Name)//'_var') 70 END IF 71 END IF 72 73 IF( ParEnv % MyPe == 0 ) THEN 74 MarkFailed = ListGetLogical( Solver % Values,'Mark Failed Strategy',GotIt) 75 AvoidFailed = ListGetLogical( Solver % Values,'Avoid Failed Strategy',GotIt) 76 IF(.NOT. GotIt) AvoidFailed = MarkFailed 77 78 IF( MarkFailed .OR. AvoidFailed ) THEN 79 LineInd = ListGetInteger( Solver % Values,'Line Marker',GotIt) 80 IF(.NOT. GotIt) THEN 81 CALL Fatal('SaveScalars_init','Failed strategy marked requires > Line Marker <') 82 END IF 83 Name = 'FINISHED_MARKER_'//TRIM(I2S(LineInd)) 84 END IF 85 86 IF( AvoidFailed ) THEN 87 INQUIRE(FILE=TRIM(Name),EXIST=GotIt) 88 IF( GotIt ) THEN 89 OPEN(NEWUNIT=MarkerUnit, FILE=Name) 90 READ(MarkerUnit,*) i 91 IF( i == 0 ) THEN 92 CALL Fatal('SaveScalars_init','Strategy already failed before!') 93 END IF 94 CLOSE(MarkerUnit) 95 END IF 96 END IF 97 98 ! Save a negative status during the execution such that if the 99 ! program terminates the negative status will prevail 100 IF( MarkFailed ) THEN 101 CALL Info('SaveScalars_init','Saving False marker at start') 102 i = 0 103 OPEN(NEWUNIT=MarkerUnit,FILE=Name,STATUS='Unknown') 104 WRITE(MarkerUnit,'(I0)') i 105 CLOSE(MarkerUnit) 106 END IF 107 END IF 108 109 110 111END SUBROUTINE SaveScalars_init 112 113 114 115!------------------------------------------------------------------------------ 116!> This subroutine saves scalar values in ascii format to an external file. 117!> This is a dynamically loaded solver with a standard interface. 118!------------------------------------------------------------------------------ 119SUBROUTINE SaveScalars( Model,Solver,dt,TransientSimulation ) 120!------------------------------------------------------------------------------ 121 USE DefUtils 122 USE Interpolation 123 124 IMPLICIT NONE 125!------------------------------------------------------------------------------ 126 TYPE(Solver_t), TARGET :: Solver 127 TYPE(Model_t) :: Model 128 REAL(KIND=dp) :: dt 129 LOGICAL :: TransientSimulation 130!------------------------------------------------------------------------------ 131! Local variables 132!------------------------------------------------------------------------------ 133 TYPE(Solver_t), POINTER :: ParSolver 134 TYPE(ValueList_t), POINTER :: Params 135 TYPE(ValueListEntry_t), POINTER :: Lst 136 TYPE(Variable_t), POINTER :: Var, OldVar, Var2, Var3 137 TYPE(Mesh_t), POINTER :: Mesh 138 TYPE(Element_t),POINTER :: CurrentElement 139 TYPE(Nodes_t) :: ElementNodes 140 LOGICAL :: MovingMesh, GotCoeff, & 141 GotIt, GotOper, GotParOper, GotVar, GotOldVar, ExactCoordinates, VariablesExist, & 142 ComplexEigenVectors, ComplexEigenValues, IsParallel, ParallelWrite, SaveCVS, & 143 FileAppend, SaveEigenValue, SaveEigenFreq, IsInteger, ParallelReduce, WriteCore, & 144 Hit, SaveToFile, EchoValues, GotAny, BodyOper, BodyForceOper, & 145 MaterialOper, MaskOper, GotMaskName, GotOldOper, ElementalVar, ComponentVar, & 146 Numbering, NodalOper, GotNodalOper, SaveFluxRange, PosOper, NegOper, SaveJson, & 147 StartNewFile 148 LOGICAL, POINTER :: ValuesInteger(:) 149 LOGICAL, ALLOCATABLE :: ActiveBC(:) 150 151 REAL (KIND=DP) :: Minimum, Maximum, AbsMinimum, AbsMaximum, & 152 Mean, Variance, MinDist, x, y, z, Vol, Intmean, intvar, & 153 KineticEnergy, PotentialEnergy, & 154 Coords(3), LocalCoords(3), TempCoordinates(3), Val, Val2, & 155 Change = 0._dp, Norm = 0.0_dp, PrevNorm, ParallelHits, ParallelCands 156 REAL (KIND=DP), ALLOCATABLE :: Values(:), & 157 CoordinateBasis(:), ElementValues(:), BoundaryFluxes(:),BoundaryAreas(:) 158 REAL (KIND=DP), POINTER :: PointCoordinates(:,:), LineCoordinates(:,:), WrkPntr(:) 159 INTEGER, ALLOCATABLE :: BoundaryHits(:) 160 INTEGER, POINTER :: PointIndex(:), NodeIndexes(:), CoordinateIndex(:), CoordinatesElemNo(:) 161 INTEGER, ALLOCATABLE, TARGET :: ClosestIndex(:) 162 CHARACTER(LEN=MAX_NAME_LEN), ALLOCATABLE :: ValueNames(:) 163 CHARACTER(LEN=MAX_NAME_LEN) :: ScalarsFile, ScalarNamesFile, DateStr, & 164 VariableName, OldVariableName, ResultPrefix, Suffix, Oper, Oper0, OldOper0, ParOper, Name, & 165 CoefficientName, ScalarParFile, OutputDirectory, MinOper, MaxOper, & 166 MaskName, OldMaskName, SaveName 167 INTEGER :: i,j,k,l,lpar,q,n,ierr,No,NoPoints,NoCoordinates,NoLines,NumberOfVars,& 168 NoDims, NoDofs, NoOper, NoElements, NoVar, NoValues, PrevNoValues, DIM, & 169 MaxVars, NoEigenValues, Ind, EigenDofs, LineInd, NormInd, CostInd, istat, nlen, & 170 jsonpos 171 INTEGER :: IntVal, FirstInd, LastInd, ScalarsUnit, MarkerUnit, NamesUnit 172 LOGICAL, ALLOCATABLE :: NodeMask(:) 173 REAL (KIND=DP) :: CT, RT 174 175 SAVE :: jsonpos 176 177!------------------------------------------------------------------------------ 178 179 CALL Info('SaveScalars', '-----------------------------------------', Level=4 ) 180 CALL Info('SaveScalars','Saving scalar values of various kinds',Level=4) 181 182 183 Mesh => GetMesh() 184 DIM = CoordinateSystemDimension() 185 Params => GetSolverParams() 186 187 MovingMesh = ListGetLogical(Params,'Moving Mesh',GotIt ) 188 189 FileAppend = ListGetLogical( Params,'File Append',GotIt) 190 191 SaveFluxRange = ListGetLogical( Params,'Save Flux Range',GotIt) 192 IF(.NOT. GotIt) SaveFluxRange = .TRUE. 193 194 ScalarsFile = ListGetString(Params,'Filename',SaveToFile ) 195 IF( SaveToFile ) THEN 196 ! Optionally number files by the number of partitions 197 ! This makes the benchmarking more convenient since each case 198 ! may use the same command file 199 IF(ListGetLogical(Params,'Partition Numbering',GotIt)) THEN 200 i = INDEX( ScalarsFile,'.',.TRUE. ) 201 j = LEN_TRIM(ScalarsFile) 202 IF(i > 0) THEN 203 Suffix = ScalarsFile(i:j) 204 WRITE( ScalarsFile,'(A,I0,A)') & 205 ScalarsFile(1:i-1)//'_',ParEnv % PEs,Suffix(1:j-i+1) 206 ELSE 207 WRITE( ScalarsFile,'(A,I0)') & 208 ScalarsFile(1:j)//'_',ParEnv % PEs 209 END IF 210 END IF 211 212 CALL SolverOutputDirectory( Solver, ScalarsFile, OutputDirectory ) 213 ScalarsFile = TRIM(OutputDirectory)// '/' //TRIM(ScalarsFile) 214 215 Numbering = ListGetLogical(Params,'Filename Numbering',GotIt) 216 217 IF( Numbering ) THEN 218 IF( Solver % TimesVisited > 0 ) THEN 219 ScalarsFile = NextFreeFilename( ScalarsFile, LastExisting = .TRUE. ) 220 ELSE 221 ScalarsFile = NextFreeFilename( ScalarsFile ) 222 END IF 223 END IF 224 225 ScalarNamesFile = TRIM(ScalarsFile) // TRIM(".names") 226 SaveCVS = ListGetLogical(Params,'Live Graph',GotIt) 227 IF(.NOT. GotIt) SaveCVS = ListGetLogical(Params,'CVS Format',GotIt) 228 SaveJson = ListGetLogical(Params,'Json Format',GotIt) 229 END IF 230 231 EchoValues = ListGetLogical( Params,'Echo Values',GotIt) 232 IF(.NOT. GotIt) EchoValues = .NOT. SaveToFile 233 234 ResultPrefix = ListGetString(Params,'Scalars Prefix',GotIt ) 235 IF(.NOT. gotIt) ResultPrefix = 'res:' 236 237 238 IsParallel = .FALSE. 239 WriteCore = .TRUE. 240 ParallelWrite = .FALSE. 241 ParallelReduce = .FALSE. 242 243 IF( ParEnv % PEs > 1 ) THEN 244 IsParallel = .TRUE. 245 ParallelReduce = GetLogical( Params,'Parallel Reduce',GotIt) 246 ParallelWrite = .NOT. ParallelReduce 247 IF( ParEnv % MyPe > 0 .AND. ParallelReduce ) THEN 248 EchoValues = .FALSE. 249 WriteCore = .FALSE. 250 END IF 251 !OutputPE = ParEnv % MYPe 252 END IF 253 254 NoLines = 0 255 LineCoordinates => ListGetConstRealArray(Params,'Polyline Coordinates',gotIt) 256 IF(gotIt) THEN 257 NoLines = SIZE(LineCoordinates,1) / 2 258 NoDims = SIZE(LineCoordinates,2) 259 END IF 260 261 NoPoints = 0 262 PointIndex => ListGetIntegerArray( Params,'Save Points',GotIt) 263 IF ( gotIt ) NoPoints = SIZE(PointIndex) 264 265 NoCoordinates = 0 266 NoElements = 0 267 PointCoordinates => ListGetConstRealArray(Params,'Save Coordinates',gotIt) 268 IF(gotIt) THEN 269 NoDims = SIZE(PointCoordinates,2) 270 ExactCoordinates = ListGetLogical(Params,'Exact Coordinates',GotIt ) 271 272 IF( ParallelReduce .AND. .NOT. ExactCoordinates) THEN 273 CALL Warn('SaveScalars','Only Exact Save Coordinates works in parallel, enforcing...') 274 ExactCoordinates = .TRUE. 275 END IF 276 277 IF(ExactCoordinates) THEN 278 ! Look for the value at the given coordinate point really. 279 NoElements = SIZE(PointCoordinates,1) 280 GotIt = .FALSE. 281 IF( .NOT. MovingMesh ) THEN 282 CoordinatesElemNo => ListGetIntegerArray( Params,'Save Coordinate Elements',GotIt ) 283 END IF 284 IF(.NOT. GotIt ) THEN 285 CALL Info('SaveScalars','Searching for elements containing save coordinates',Level=8) 286 287 ALLOCATE(ClosestIndex(NoElements), STAT=istat) 288 IF( istat /= 0 ) CALL Fatal('SaveScalars','Memory allocation error for CoordinateElemNo') 289 DO j=1,NoElements 290 Coords(1:NoDims) = PointCoordinates(j,1:NoDims) 291 IF(NoDims < 3 ) Coords(NoDims+1:3) = 0.0_dp 292 ClosestIndex(j) = ClosestElementInMesh( Mesh, Coords ) 293 END DO 294 295 CoordinatesElemNo => ClosestIndex 296 IF( .NOT. MovingMesh ) THEN 297 CALL ListAddIntegerArray( Params,'Save Coordinate Elements',& 298 NoElements,ClosestIndex ) 299 END IF 300 END IF 301 ELSE 302 ! Find the indexes of minimum distances 303 NoCoordinates = SIZE(PointCoordinates,1) 304 GotIt = .FALSE. 305 IF( .NOT. MovingMesh ) THEN 306 CoordinateIndex => ListGetIntegerArray( Params,'Save Coordinate Indexes',GotIt ) 307 END IF 308 IF( .NOT. GotIt ) THEN 309 CALL Info('SaveScalars','Searching for closest nodes to coordinates',Level=8) 310 ALLOCATE(ClosestIndex(NoCoordinates), STAT=istat) 311 IF( istat /= 0 ) CALL Fatal('SaveScalars','Memory allocation error for CoordinateIndex') 312 DO j=1,NoCoordinates 313 Coords(1:NoDims) = PointCoordinates(j,1:NoDims) 314 IF(NoDims < 3 ) Coords(NoDims+1:3) = 0.0_dp 315 ClosestIndex(j) = ClosestNodeInMesh( Mesh, Coords ) 316 END DO 317 318 CoordinateIndex => ClosestIndex 319 IF( .NOT. MovingMesh ) THEN 320 CALL ListAddIntegerArray( Params,'Save Coordinate Indexes',& 321 NoCoordinates,ClosestIndex ) 322 END IF 323 END IF 324 END IF 325 END IF 326 327!------------------------------------------------------------------------------ 328 329 n = Mesh % MaxElementNodes 330 ALLOCATE( ElementNodes % x(n), ElementNodes % y(n), ElementNodes % z(n), & 331 ElementValues( n ), CoordinateBasis(n), STAT=istat) 332 IF( istat /= 0 ) CALL Fatal('SaveScalars','Memory allocation error 1') 333 334 n = MAX( Model % NumberOfBodies, MAX(Model % NumberOfBCs, NoLines)) 335 ALLOCATE( BoundaryFluxes(n), BoundaryAreas(n), BoundaryHits(n), STAT=istat ) 336 IF( istat /= 0 ) CALL Fatal('SaveScalars','Memory allocation error 2') 337 338 ALLOCATE( ActiveBC( Model % NumberOfBCs ), STAT=istat ) 339 IF( istat /= 0 ) CALL Fatal('SaveScalars','Memory allocation error 3') 340 341 342 ComplexEigenVectors = ListGetLogical(Params,'Complex Eigen Vectors',GotIt) 343 344 345 346 !------------------------------------------------------------------------------ 347 ! Go through the variables and compute the desired statistical data 348 !------------------------------------------------------------------------------ 349 NoValues = 0 350 GotVar = .TRUE. 351 GotOper = .FALSE. 352 GotOldOper = .FALSE. 353 NULLIFY(OldVar) 354 NoVar = 0 355 MinOper = 'min' 356 MaxOper = 'max' 357 GotNodalOper = .FALSE. 358 OldMaskName = 'save scalars' 359 PosOper = .FALSE. 360 NegOper = .FALSE. 361 362 363 DO WHILE(GotVar .OR. GotOper) 364 365 GotOper = .FALSE. 366 NULLIFY(Var) 367 368 NoVar = NoVar + 1 369 WRITE (Name,'(A,I0)') 'Variable ',NoVar 370 371 VariableName = ListGetString( Params, TRIM(Name), GotVar ) 372 373 374 IF(TRIM(VariableName) == 'cpu time' .OR. TRIM(VariableName) == 'cpu memory') THEN 375 CALL Warn('SaveScalars','This variable should now be invoked as an operator: '//TRIM(VariableName)) 376 CYCLE 377 END IF 378 379 GotOldVar = .FALSE. 380 381 IF(GotVar) THEN 382 Var => VariableGet( Model % Variables, TRIM(VariableName) ) 383 IF ( .NOT. ASSOCIATED( Var ) ) THEN 384 Var => VariableGet( Model % Variables, TRIM(VariableName)//' 1' ) 385 IF( ASSOCIATED( Var ) ) THEN 386 CALL Info('SaveScalars','Treating a component variable: '//TRIM(VariableName),Level=8) 387 ComponentVar = .TRUE. 388 Var2 => VariableGet( Model % Variables, TRIM(VariableName)//' 2' ) 389 Var3 => VariableGet( Model % Variables, TRIM(VariableName)//' 3' ) 390 ELSE 391 CALL Warn('SaveScalars','Requested variable does not exist: '//TRIM(VariableName)) 392 CYCLE 393 END IF 394 ELSE 395 ComponentVar = .FALSE. 396 END IF 397 OldVar => Var 398 OldVariableName = VariableName 399 400 ! A 0D variable cannot really be much operated, hence save it as is 401 !------------------------------------------------------------------- 402 IF(SIZE(Var % Values) == Var % Dofs) THEN 403 IsInteger = .FALSE. 404 IF( VariableName == 'timestep' ) IsInteger = .TRUE. 405 IF( VariableName == 'nonlin iter' ) IsInteger = .TRUE. 406 IF( VariableName == 'coupled iter' ) IsInteger = .TRUE. 407 IF( VariableName == 'run' ) IsInteger = .TRUE. 408 409 IF( Var % Dofs == 1 ) THEN 410 CALL AddToSaveList('value: '//TRIM(VariableName)//' scalar variable', & 411 Var % Values(1), IsInteger ) 412 ELSE 413 DO j=1,Var % DOfs 414 CALL AddToSaveList('value: '//ComponentName(VariableName,j)//' scalar variable', Var % Values(j)) 415 END DO 416 END IF 417 CYCLE 418 END IF 419 420 WRITE (Name,'(A,I0)') 'Nodal Variable ',NoVar 421 NodalOper = ListGetLogical(Params,TRIM(Name),GotNodalOper) 422 ELSE 423 IF(ASSOCIATED(OldVar)) THEN 424 Var => OldVar 425 VariableName = OldVariableName 426 GotOldVar = .TRUE. 427 END IF 428 END IF 429 430 WRITE (Name,'(A,I0)') 'Nodal Variable ',NoVar 431 IF( ListCheckPresent( Params,TRIM(Name) ) ) THEN 432 NodalOper = ListGetLogical(Params,TRIM(Name),GotNodalOper) 433 END IF 434 435 NoOper = NoVar 436 MaskOper = .FALSE. 437 WRITE (Name,'(A,I0)') 'Operator ',NoOper 438 Oper0 = ListGetString(Params,TRIM(Name),GotOper) 439 IF(.NOT. GotOper .AND. GotOldOper ) Oper0 = OldOper0 440 441 442 WRITE (Name,'(A,I0)') 'Mask Name ',NoOper 443 MaskName = ListGetString(Params,TRIM(Name),GotMaskName) 444 IF( GotMaskName ) THEN 445 OldMaskName = MaskName 446 ELSE 447 MaskName = OldMaskName 448 END IF 449 450 451 IF(.NOT. (GotOper .OR. GotVar .OR. GotMaskName ) ) CYCLE 452 453 454 IF( ASSOCIATED( Var ) ) THEN 455 CALL Info('SaveScalars','Treating variable: '//TRIM(VariableName),Level=12) 456 ElementalVar = ( Var % TYPE == Variable_on_nodes_on_elements ) 457 END IF 458 459 IF( GotOper ) THEN 460 CALL Info('SaveScalars','Treating operator: '//TRIM(Oper0),Level=12) 461 OldOper0 = Oper0 462 GotOldOper = .TRUE. 463 ELSE IF( GotOldOper ) THEN 464 Oper0 = OldOper0 465 GotOper = GotOldOper 466 ELSE 467 CALL Info('SaveScalars','No operator given for variable: '//TRIM(VariableName)) 468 CYCLE 469 END IF 470 471 472 BodyOper = .FALSE. 473 BodyForceOper = .FALSE. 474 MaterialOper = .FALSE. 475 nlen = LEN_TRIM(Oper0) 476 IF( Oper0(1:11) == 'body force ') THEN 477 BodyForceOper = .TRUE. 478 Oper = Oper0(12:nlen) 479 ELSE IF( Oper0(1:5) == 'body ') THEN 480 BodyOper = .TRUE. 481 Oper = Oper0(6:nlen) 482 ELSE IF( Oper0(1:9) == 'material ') THEN 483 MaterialOper = .TRUE. 484 Oper = Oper0(10:nlen) 485 ELSE 486 Oper = Oper0 487 END IF 488 MaskOper = ( BodyForceOper .OR. BodyOper .OR. MaterialOper ) 489 IF( MaskOper ) THEN 490 CALL Info('SaveScalars','Operator to be masked: '//TRIM(Oper),Level=12) 491 END IF 492 493 nlen = LEN_TRIM(Oper) 494 PosOper = .FALSE. 495 NegOper = .FALSE. 496 IF( Oper(1:9) == 'positive ') THEN 497 PosOper = .TRUE. 498 Oper = Oper(10:nlen) 499 ELSE IF( Oper(1:9) == 'negative ') THEN 500 NegOper = .TRUE. 501 Oper = Oper(10:nlen) 502 END IF 503 504 505 WRITE (Name,'(A,I0)') 'Coefficient ',NoOper 506 CoefficientName = ListGetString(Params,TRIM(Name),GotCoeff ) 507 508 WRITE (Name,'(A,I0)') 'Parallel Operator ',NoOper 509 ParOper = ListGetString(Params,TRIM(Name),GotParOper) 510 IF(.NOT. GotParOper) ParOper = OperToParOperMap(Oper) 511 512 IF( MaskOper ) THEN 513 GotIt = .FALSE. 514 IF( BodyOper ) THEN 515 GotIt = ListGetLogicalAnyBody( Model, MaskName ) 516 ELSE IF( BodyForceOper ) THEN 517 GotIt = ListGetLogicalAnyBodyForce( Model, MaskName ) 518 ELSE IF( MaterialOper ) THEN 519 GotIt = ListGetLogicalAnyMaterial( Model, MaskName ) 520 END IF 521 IF(.NOT. GotIt ) THEN 522 CALL Warn('SaveScalars','Masked operators require mask: '//TRIM(MaskName)) 523 END IF 524 END IF 525 526 ActiveBC = .FALSE. 527 DO j=1,Model % NumberOfBCs 528 ActiveBC(j) = & 529 ListGetLogical(Model % BCs(j) % Values,'Flux Integrate',gotIt) .OR. & 530 ListGetLogical(Model % BCs(j) % Values, MaskName, gotIt) 531 END DO 532 533 IF ( GotOper ) THEN 534 535 SELECT CASE( Oper ) 536 537 CASE ('partitions') 538 CASE ('partition checksum') 539 CASE ('partition neighbours checksum') 540 CASE ('cpu time') 541 CASE ('wall time') 542 CASE ('cpu memory') 543 CASE ('nodes') 544 CASE ('elements') 545 CASE ('bounding box') 546 CASE ('area') 547 CASE ('volume') 548 CASE ('threads') 549 550 CASE DEFAULT 551 IF( .NOT. (GotVar .OR. GotOldVar ) ) THEN 552 CALL Error('SaveScalars','Operator > '//TRIM(Oper)//' < requires variable') 553 CYCLE 554 END IF 555 END SELECT 556 557 558 ! Set default name for saving 559 IF(GotVar .OR. GotOldVar ) THEN 560 SaveName = TRIM(Oper0)//': '//TRIM(VariableName) 561 IF( GotMaskName ) THEN 562 SaveName = TRIM(SaveName)//' mask '//TRIM(MaskName) 563 END IF 564 IF( GotNodalOper ) THEN 565 IF( NodalOper ) THEN 566 SaveName = TRIM(SaveName)//' nodal' 567 ELSE 568 SaveName = TRIM(SaveName)//' non-nodal' 569 END IF 570 END IF 571 END IF 572 573 574 SELECT CASE(Oper) 575 576 CASE ('partitions') 577 Val = 1.0_dp * ParEnv % PEs 578 SaveName = 'value: number of partitions' 579 CALL AddToSaveList(SaveName,Val,.TRUE.,ParOper) 580 581 CASE ('threads') 582 Val = 1.0_dp * ParEnv % NumberOfThreads 583 SaveName = 'value: number of threads' 584 CALL AddToSaveList(SaveName,Val,.TRUE.,ParOper) 585 586 CASE ('partition checksum') 587 Val = 0.0_dp 588 IF( IsParallel ) THEN 589 Val = 1.0_dp * SUM( 1.0_dp * Mesh % ParallelInfo % GlobalDOFS ) 590 ! Give different partition different weight to create something like a checksum 591 Val = ( ParEnv % MyPe + 1 ) * Val 592 END IF 593 SaveName = 'value: partition checksum' 594 ! Don't use integer as type because it can exceed the bounds, long should be used 595 CALL AddToSaveList(SaveName,Val,.FALSE.,ParOper) 596 597 CASE ('partition neighbours checksum') 598 Val = 0.0_dp 599 IF( IsParallel ) THEN 600 DO j=1,Mesh % NumberOfNodes 601 Val = Val + 1.0_dp * SUM( Mesh % ParallelInfo % NeighbourList(j) % Neighbours ) 602 END DO 603 Val = ( ParEnv % MyPe + 1 ) * Val 604 END IF 605 SaveName = 'value: partition neighbours checksum' 606 ! Don't use integer as type because it can exceed the bounds, long should be used 607 CALL AddToSaveList(SaveName,Val,.FALSE.,ParOper) 608 609 CASE ('cpu time') 610 Val = CPUTime() 611 SaveName = 'value: cpu time (s)' 612 CALL AddToSaveList(SaveName,Val,.FALSE.,ParOper) 613 614 CASE ('wall time') 615 Val = RealTime() 616 SaveName = 'value: real time (s)' 617 CALL AddToSaveList(SaveName,Val,.FALSE.,ParOper) 618 619 CASE ('cpu memory') 620 Val = CPUMemory() 621 SaveName = 'value: maximum memory usage (kb)' 622 CALL AddToSaveList(SaveName,Val,.FALSE.,ParOper) 623 624 CASE ('nodes') 625 Val = 1.0_dp * Solver % Mesh % NumberOfNodes 626 SaveName = TRIM(Oper)//': '//TRIM(VariableName) 627 CALL AddToSaveList(SaveName,Val,.TRUE.,ParOper) 628 629 CASE ('elements') 630 Val = 1.0_dp * Solver % Mesh % NumberOfBulkElements 631 SaveName = TRIM(Oper)//': '//TRIM(VariableName) 632 CALL AddToSaveList(SaveName,Val,.TRUE.,ParOper) 633 634 CASE ('bounding box') 635 Val = MINVAL( Solver % Mesh % Nodes % x ) 636 CALL AddToSaveList(TRIM(Oper)//' min x',Val,ParallelOperator=MinOper) 637 Val = MAXVAL( Solver % Mesh % Nodes % x ) 638 CALL AddToSaveList(TRIM(Oper)//' max x',Val,ParallelOperator=MaxOper) 639 Val = MINVAL( Solver % Mesh % Nodes % y ) 640 CALL AddToSaveList(TRIM(Oper)//' min y',Val,ParallelOperator=MinOper) 641 Val = MAXVAL( Solver % Mesh % Nodes % y ) 642 CALL AddToSaveList(TRIM(Oper)//' max y',Val,ParallelOperator=MaxOper) 643 IF( Solver % Mesh % MeshDim > 2 ) THEN 644 Val = MINVAL( Solver % Mesh % Nodes % z ) 645 CALL AddToSaveList(TRIM(Oper)//' min z',Val,ParallelOperator=MinOper) 646 Val = MAXVAL( Solver % Mesh % Nodes % z ) 647 CALL AddToSaveList(TRIM(Oper)//' max z',Val,ParallelOperator=MaxOper) 648 END IF 649 650 ! Operators that require a variable (except for 'area' and 'volume') 651 !----------------------------------------------------------------- 652 CASE ('norm') 653 Val = Var % Norm 654 CALL AddToSaveList(SaveName,Val,.FALSE.,ParOper) 655 656 CASE ('nonlin change') 657 Val = Var % NonlinChange 658 CALL AddToSaveList(SaveName,Val,.FALSE.,ParOper) 659 660 CASE ('steady state change') 661 Val = Var % SteadyChange 662 CALL AddToSaveList(SaveName,Val,.FALSE.,ParOper) 663 664 CASE ('nonlin iter') 665 Val = Var % NonlinIter 666 CALL AddToSaveList(SaveName,Val,.TRUE.,ParOper) 667 668 CASE ('nonlin converged') 669 Val = 1.0_dp * Var % NonlinConverged 670 CALL AddToSaveList(SaveName,Val,.TRUE.,ParOper) 671 672 CASE ('steady converged') 673 Val = 1.0_dp * Var % SteadyConverged 674 CALL AddToSaveList(SaveName,Val,.TRUE.,ParOper) 675 676 CASE ('dofs') 677 Val = 1.0_dp * SIZE(Var % Values) 678 CALL AddToSaveList(SaveName,Val,.TRUE.,ParOper) 679 680 CASE ('nans') ! number of NaN:s in vector 681 GotIt = .FALSE. 682 j = 0 683 DO i=1,SIZE(Var % Values) 684 IF( Var % Values(i) /= Var % Values(i) ) THEN 685 j = j + 1 686 END IF 687 END DO 688 Val = 1.0_dp * j 689 CALL AddToSaveList(SaveName,Val,.TRUE.,ParOper) 690 691 CASE ('sum','sum abs','mean abs','max','max abs','min','min abs',& 692 'mean','variance','range','sum square','mean square','rms') 693 IF( MaskOper ) CALL CreateNodeMask() 694 IF( GotNodalOper ) THEN 695 Val = VectorStatistics(Var,Oper,NodalOper) 696 ELSE 697 Val = VectorStatistics(Var,Oper) 698 END IF 699 CALL AddToSaveList(SaveName,Val,.FALSE.,ParOper) 700 701 CASE ('deviation') 702 IF( MaskOper ) CALL CreateNodeMask() 703 Val = VectorMeanDeviation(Var,Oper) 704 CALL AddToSaveList(SaveName, Val,.FALSE.,ParOper) 705 706 CASE ('int','int mean','int square','int square mean','int rms','int abs','int abs mean',& 707 'int variance','volume','potential energy','diffusive energy','convective energy') 708 709 IF( MaskOper ) CALL CreateNodeMask() 710 Val = BulkIntegrals(Var, Oper, GotCoeff, CoefficientName) 711 IF(GotCoeff) THEN 712 SaveName = TRIM(SaveName)//' with '//TRIM(CoefficientName) 713 END IF 714 CALL AddToSaveList(SaveName, Val,.FALSE.,ParOper) 715 716 CASE('boundary sum','boundary dofs','boundary max','boundary max abs','boundary min',& 717 'boundary min abs','boundary mean') 718 719 IF( .NOT. ANY( ActiveBC ) ) THEN 720 CALL Error('SaveScalars','No flag > '//TRIM(MaskName)// & 721 ' < active for operator: '// TRIM(Oper)) 722 CYCLE 723 END IF 724 725 BoundaryHits = 0 726 BoundaryFluxes = 0.0_dp 727 CALL BoundaryStatistics(Var, Oper, GotCoeff, & 728 CoefficientName, BoundaryFluxes, BoundaryHits) 729 730 GotAny = .FALSE. 731 DO j=1,Model % NumberOfBCs 732 IF( ActiveBC(j) ) THEN 733 IF( TRIM(Oper) == 'boundary mean' ) THEN 734 IF( IsParallel .AND. ParallelReduce ) THEN 735 CALL Warn('SaveScalars','Operator > boundary mean < not implemented in parallel!') 736 ELSE IF( BoundaryHits(j) > 0 ) THEN 737 BoundaryFluxes(j) = BoundaryFluxes(j) / BoundaryHits(j) 738 END IF 739 END IF 740 WRITE (Name,'(A,A,A,A,I0)') TRIM(Oper),': ',TRIM(VariableName),' over bc ',j 741 CALL AddToSaveList( TRIM(Name), BoundaryFluxes(j),.FALSE.,ParOper) 742 END IF 743 END DO 744 745 CASE ('boundary int','boundary int mean','area','diffusive flux','convective flux') 746 747 IF( .NOT. ANY( ActiveBC ) ) THEN 748 CALL Error('SaveScalars','No flag > '//TRIM(MaskName)// & 749 '< active for operator: '// TRIM(Oper)) 750 ELSE 751 BoundaryHits = 0 752 BoundaryFluxes = 0.0_dp 753 BoundaryAreas = 0.0_dp 754 IF( SaveFluxRange ) THEN 755 Minimum = HUGE(Minimum) 756 Maximum = -HUGE(Maximum) 757 END IF 758 759 CALL BoundaryIntegrals(Var, Oper, GotCoeff, CoefficientName,& 760 BoundaryFluxes,BoundaryAreas,BoundaryHits) 761 762 DO j=1,Model % NumberOfBCs 763 IF( ActiveBC(j) ) THEN 764 IF( TRIM(Oper) == 'boundary int mean' ) THEN 765 IF( IsParallel .AND. ParallelReduce ) THEN 766 CALL Warn('SaveScalars','Operator > boundary int mean < not implemented in parallel!') 767 ELSE IF( BoundaryAreas(j) > 0.0 ) THEN 768 BoundaryFluxes(j) = BoundaryFluxes(j) / BoundaryAreas(j) 769 END IF 770 END IF 771 WRITE (Name,'(A,A,A,A,I0)') TRIM(Oper),': ',TRIM(VariableName),' over bc ',j 772 CALL AddToSaveList( TRIM(Name), BoundaryFluxes(j),.FALSE.,ParOper) 773 END IF 774 END DO 775 776 IF( SaveFluxRange ) THEN 777 IF(TRIM(Oper) == 'diffusive flux' .OR. TRIM(Oper) == 'convective flux') THEN 778 ParOper = 'min' 779 WRITE (Name,'(A,A,A,A)') 'min ',TRIM(Oper),': ',TRIM(VariableName) 780 CALL AddToSaveList( TRIM(Name), Minimum,.FALSE.,ParOper) 781 782 ParOper = 'max' 783 WRITE (Name,'(A,A,A,A)') 'max ',TRIM(Oper),': ',TRIM(VariableName) 784 CALL AddToSaveList( TRIM(Name), Maximum,.FALSE.,ParOper) 785 END IF 786 END IF 787 END IF 788 789 790 IF(NoLines > 0) THEN 791 BoundaryHits = 0 792 BoundaryFluxes = 0.0_dp 793 BoundaryAreas = 0.0_dp 794 795 CALL PolylineIntegrals(Var, Oper, GotCoeff, CoefficientName,& 796 BoundaryFluxes,BoundaryAreas, BoundaryHits) 797 798 DO j=1,NoLines 799 IF( TRIM(Oper) == 'boundary int mean' ) THEN 800 IF( BoundaryHits(j) > 0 ) THEN 801 BoundaryFluxes(j) = BoundaryFluxes(j) / BoundaryAreas(j) 802 END IF 803 END IF 804 WRITE (Name,'(A,A,A,A,I0)') TRIM(Oper),': ',TRIM(VariableName),' over polyline ',j 805 CALL AddToSaveList( TRIM(Name), BoundaryFluxes(j),.FALSE.,ParOper) 806 END DO 807 END IF 808 809 CASE DEFAULT 810 811 WRITE (Message,'(A,A)') 'Unknown operator: ',TRIM(Oper) 812 CALL WARN('SaveScalars',Message) 813 814 END SELECT 815 816 END IF 817 END DO 818 819 IF( NoVar > 0 ) THEN 820 WRITE (Message,'(A,I0,A)') 'Performed ',NoVar-1,' reduction operations' 821 CALL Info('SaveScalars',Message,Level=7) 822 END IF 823 824 825 !------------------------------------------------------------------------------ 826 ! Add eigenvalues on the list if told to 827 !------------------------------------------------------------------------------ 828 829 SaveEigenValue = ListGetLogical( Params, 'Save Eigenvalues', GotIt ) 830 IF(.NOT. GotIt) & 831 SaveEigenValue = ListGetLogical( Params, 'Save Eigen values', GotIt ) 832 833 SaveEigenFreq = ListGetLogical( Params, 'Save Eigenfrequencies', GotIt ) 834 IF(.NOT. GotIt) & 835 SaveEigenFreq = ListGetLogical( Params, 'Save Eigen Frequencies', GotIt ) 836 837 IF ( SaveEigenValue .OR. SaveEigenFreq ) THEN 838 ComplexEigenValues = ListGetLogical(Params,'Complex Eigen Values',GotIt) 839 IF(.NOT. GotIt) & 840 ComplexEigenValues = ListGetLogical(Params,'Complex EigenValues',GotIt) 841 842 l = 0 843 DO i = 1, Model % NumberOfSolvers 844 IF ( Model % Solvers(i) % NOFEigenValues > 0 ) THEN 845 DO k = 1, Model % Solvers(i) % NOFEigenValues 846 847 Val = REAL( Model % Solvers(i) % Variable % EigenValues(k) ) 848 l = l + 1 849 850 IF ( SaveEigenValue ) THEN 851 WRITE( Name, '("eigen: Re value ", I0)' ) k 852 CALL AddToSaveList(TRIM(Name), Val) 853 IF( ComplexEigenValues ) THEN 854 Val2 = AIMAG( Model % Solvers(i) % Variable % EigenValues(k) ) 855 WRITE( Name, '("eigen: Im value ", I0)' ) k 856 CALL AddToSaveList(TRIM(Name), Val2) 857 END IF 858 END IF 859 860 IF ( SaveEigenFreq ) THEN 861 WRITE( Name, '("eigen: frequency ", I0, " [Hz]")' ) k 862 IF ( Val >= 0.0 ) THEN 863 Val = SQRT(Val) / (2*PI) 864 ELSE 865 ! If the eigenvalue is negative take take take a square root of its absolute value and 866 ! return a negative frequency. 867 Val = -SQRT(-Val) / (2*PI) 868 END IF 869 CALL AddToSaveList(TRIM(Name), Val) 870 END IF 871 872 END DO 873 END IF 874 END DO 875 WRITE (Message,'(A,I0,A)') 'Found ',l,' Eigenvalues' 876 CALL Info('SaveScalars',Message) 877 END IF 878 879 !------------------------------------------------------------------------------ 880 ! Get the info at given node points 881 !------------------------------------------------------------------------------ 882 DO k=1,NoPoints+NoCoordinates 883 884 IF(k <= NoPoints) THEN 885 l = PointIndex(k) 886 ELSE 887 l = CoordinateIndex(k-NoPoints) 888 END IF 889 890 Var => Model % Variables 891 892 DO WHILE( ASSOCIATED( Var ) ) 893 894 IF ( .NOT. Var % Output .OR. SIZE(Var % Values) == Var % DOFs ) THEN 895 896 CONTINUE 897 898 ELSE IF (ASSOCIATED (Var % EigenVectors)) THEN 899 900 NoEigenValues = SIZE(Var % EigenValues) 901 EigenDofs = SIZE( Var % EigenVectors(1,:) ) / SIZE( Var % Perm ) 902 903 IF(EigenDofs == Var % DOFs) THEN 904 DO j=1,NoEigenValues 905 DO i=1,Var % DOFs 906 Ind = Var % Perm(l) 907 IF( Ind > 0 ) THEN 908 909 Val = REAL( Var % EigenVectors(j, Var%Dofs*(Ind-1)+i) ) 910 IF(Var % DOFs == 1) THEN 911 WRITE(Name,'("value: Re Eigen ",I0," ",A," at node ",I7)') j,TRIM(Var % Name),l 912 ELSE 913 WRITE(Name,'("value: Re Eigen ",I0," ",A,I2," at node ",I7)') j,TRIM(Var % Name),i,l 914 END IF 915 CALL AddToSaveList( TRIM(Name), Val) 916 917 IF(ComplexEigenVectors) THEN 918 Val2 = AIMAG( Var % EigenVectors(j, Var%Dofs*(Ind-1)+i) ) 919 IF(Var % DOFs == 1) THEN 920 WRITE(Name,'("value: Im Eigen ",I0," ",A," at node ",I7)') j,TRIM(Var % Name),l 921 ELSE 922 WRITE(Name,'("value: Im Eigen ",I0," ",A,I2," at node ",I7)') j,TRIM(Var % Name),i,l 923 END IF 924 CALL AddToSaveList( TRIM(Name), Val2) 925 END IF 926 927 END IF 928 END DO 929 END DO 930 END IF 931 932 ELSE IF( Var % Dofs == 1) THEN 933 ! The variables exist always also as scalars, therefore omit vectors. 934 935 Ind = l 936 IF(ASSOCIATED(Var % Perm)) Ind = Var % Perm(l) 937 938 IF(Ind > 0) THEN 939 WRITE(Name,'("value: ",A," at node ",I7)') TRIM( Var % Name ), l 940 CALL AddToSaveList( TRIM(Name), Var % Values(Ind)) 941 END IF 942 943 END IF 944 945 Var => Var % Next 946 END DO 947 END DO 948 949 IF( NoPoints + NoCoordinates > 0 ) THEN 950 WRITE (Message,'(A,I0,A)') 'Tabulated all field values at ',NoPoints+NoCoordinates,' points' 951 CALL Info('SaveScalars',Message) 952 END IF 953 954 !------------------------------------------------------------------------------ 955 ! Get the info at exact coordinates within elements 956 !------------------------------------------------------------------------------ 957 958 !------------------------------------------------------------------------------ 959 ! For parallel cases set the default to -HUGE and take the max of the values 960 ParOper = 'max' 961 962 DO k=1,NoElements 963 l = CoordinatesElemNo(k) 964 965 lpar = l 966 IF( l > 0 ) THEN 967 CurrentElement => Mesh % Elements(l) 968 IF( IsParallel ) lpar = CurrentElement % GElementIndex 969 END IF 970 971 lpar = NINT( ParallelReduction(1.0_dp * lpar, 2 ) ) 972 IF( lpar == 0 ) CYCLE 973 974 IF( l > 0 ) THEN 975 n = CurrentElement % TYPE % NumberOfNodes 976 977 NodeIndexes => CurrentElement % NodeIndexes 978 979 Coords(1:NoDims) = PointCoordinates(k,1:NoDims) 980 IF(NoDims < 3 ) Coords(NoDims+1:3) = 0.0_dp 981 982 ElementNodes % x(1:n) = Mesh % Nodes % x(NodeIndexes) 983 ElementNodes % y(1:n) = Mesh % Nodes % y(NodeIndexes) 984 ElementNodes % z(1:n) = Mesh % Nodes % z(NodeIndexes) 985 986 Hit = PointInElement( CurrentElement, ElementNodes, & 987 Coords, LocalCoords, GlobalEps = 1.0_dp, LocalEps=0.1_dp ) 988 989 ElementValues(1:n) = 0.0d0 990 CoordinateBasis = 0.0_dp 991 DO q=1,N 992 ElementValues(q) = 1.0d0 993 CoordinateBasis(q) = InterpolateInElement( CurrentElement, ElementValues, & 994 LocalCoords(1), LocalCoords(2), LocalCoords(3) ) 995 ElementValues(q) = 0.0d0 996 END DO 997 END IF 998 999 Var => Model % Variables 1000 DO WHILE( ASSOCIATED( Var ) ) 1001 1002 ElementalVar = ( Var % TYPE == Variable_on_nodes_on_elements ) 1003 1004 IF ( .NOT. Var % Output .OR. SIZE(Var % Values) == Var % DOFs) THEN 1005 CONTINUE 1006 1007 ELSE IF (ASSOCIATED (Var % EigenVectors)) THEN 1008 NoEigenValues = SIZE(Var % EigenValues) 1009 EigenDofs = SIZE( Var % EigenVectors(1,:) ) / SIZE( Var % Perm ) 1010 1011 IF(EigenDofs == Var % DOFs) THEN 1012 DO j=1,NoEigenValues 1013 DO i=1,Var % DOFs 1014 1015 Val = -HUGE(Val) 1016 Val2 = -HUGE(Val) 1017 GotIt = .FALSE. 1018 1019 IF( l > 0 ) THEN 1020 IF( ElementalVar ) THEN 1021 NodeIndexes => CurrentElement % DgIndexes 1022 ELSE 1023 NodeIndexes => CurrentElement % NodeIndexes 1024 END IF 1025 1026 IF( ALL(Var % Perm(NodeIndexes(1:n)) > 0)) THEN 1027 ElementValues(1:n) = REAL( Var % EigenVectors(j,Var%Dofs*(Var % Perm(NodeIndexes(1:n))-1)+i) ) 1028 Val = SUM( CoordinateBasis(1:n) * ElementValues(1:n) ) 1029 1030 IF(ComplexEigenVectors) THEN 1031 ElementValues(1:n) = AIMAG( Var % EigenVectors(j,Var%Dofs*(Var % Perm(NodeIndexes(1:n))-1)+i) ) 1032 Val2 = SUM( CoordinateBasis(1:n) * ElementValues(1:n) ) 1033 END IF 1034 1035 GotIt = .TRUE. 1036 END IF 1037 END IF 1038 1039 IF( GotIt .OR. IsParallel ) THEN 1040 IF(Var % DOFs == 1) THEN 1041 WRITE(Name,'("value: Re Eigen ",I0," ",A," in element ",I0)') j,TRIM(Var % Name),lpar 1042 ELSE 1043 WRITE(Name,'("value: Re Eigen ",I0," ",A," ",I0," in element ",I0)') j,TRIM(Var % Name),i,lpar 1044 END IF 1045 CALL AddToSaveList(TRIM(Name), Val,.FALSE.,ParOper) 1046 1047 IF( ComplexEigenVectors ) THEN 1048 IF(Var % DOFs == 1) THEN 1049 WRITE(Name,'("value: Im Eigen ",I0," ",A," in element ",I0)') j,TRIM(Var % Name),lpar 1050 ELSE 1051 WRITE(Name,'("value: Im Eigen ",I0," ",A," ",I0," in element ",I0)') j,TRIM(Var % Name),i,lpar 1052 END IF 1053 CALL AddToSaveList(TRIM(Name), Val2,.FALSE.,ParOper) 1054 END IF 1055 END IF 1056 END DO 1057 END DO 1058 END IF 1059 1060 ELSE IF(Var % Dofs == 1) THEN 1061 1062 Val = -HUGE( Val ) 1063 GotIt = .FALSE. 1064 1065 IF( l > 0 ) THEN 1066 IF( ElementalVar ) THEN 1067 NodeIndexes => CurrentElement % DgIndexes 1068 ELSE 1069 NodeIndexes => CurrentElement % NodeIndexes 1070 END IF 1071 1072 IF( ASSOCIATED( Var % Perm ) ) THEN 1073 IF( ALL(Var % Perm(NodeIndexes(1:n)) > 0)) THEN 1074 ElementValues(1:n) = Var % Values(Var % Perm(NodeIndexes(1:n))) 1075 Val = SUM( CoordinateBasis(1:n) * ElementValues(1:n) ) 1076 GotIt = .TRUE. 1077 END IF 1078 ELSE 1079 ElementValues(1:n) = Var % Values(NodeIndexes(1:n)) 1080 Val = SUM( CoordinateBasis(1:n) * ElementValues(1:n) ) 1081 GotIt = .TRUE. 1082 END IF 1083 END IF 1084 1085 IF(GotIt .OR. IsParallel) THEN 1086 WRITE(Name,'("value: ",A," in element ",I0)') TRIM(Var % Name),lpar 1087 CALL AddToSaveList(TRIM(Name), Val,.FALSE.,ParOper) 1088 END IF 1089 END IF 1090 1091 Var => Var % Next 1092 1093 END DO 1094 END DO 1095 1096 IF( NoElements > 0 ) THEN 1097 WRITE (Message,'(A,I0,A)') 'Tabulated points within ',NoElements,' elements' 1098 CALL Info('SaveScalars',Message) 1099 END IF 1100 1101 !------------------------------------------------------------------------------ 1102 ! Update time elapsed from start if requested 1103 !------------------------------------------------------------------------------ 1104 IF( ListGetLogical( Model % Simulation,'Simulation Timing',GotIt ) ) THEN 1105 CALL Info('SaveScalars','Adding information on simulation timing',Level=10) 1106 CT = CPUTime() - ListGetConstReal( Model % Simulation,'cputime0',GotIt) 1107 RT = RealTime() - ListGetConstReal( Model % Simulation,'realtime0',GotIt) 1108 CALL AddToSaveList('simulation: cpu time (s)',CT ) 1109 CALL AddToSaveList('simulation: real time (s)',RT ) 1110 END IF 1111 1112 !------------------------------------------------------------------------------ 1113 ! Read the scalars defined in other modules 1114 !------------------------------------------------------------------------------ 1115 Lst => ListHead(Model % Simulation) 1116 l = 0 1117 DO WHILE( ASSOCIATED( Lst ) ) 1118 IF ( Lst % Name(1:4) == TRIM(ResultPrefix) ) THEN 1119 IF ( ASSOCIATED(Lst % Fvalues) ) THEN 1120 CALL AddToSaveList(Lst % Name, Lst % Fvalues(1,1,1)) 1121 l = l + 1 1122 END IF 1123 END IF 1124 Lst => Lst % Next 1125 END DO 1126 IF(l > 0 ) THEN 1127 WRITE (Message,'(A,I0,A)') 'Found ',l,' result scalars in simulation section' 1128 CALL Info('SaveScalars',Message,Level=7) 1129 END IF 1130 1131 !------------------------------------------------------------------------------ 1132 ! Add free form expressions with GetCReal 1133 !------------------------------------------------------------------------------ 1134 NoVar = 0 1135 GotVar = .TRUE. 1136 DO WHILE( GotVar ) 1137 NoVar = NoVar + 1 1138 WRITE (Name,'(A,I0)') 'Expression ',NoVar 1139 Val = ListGetCReal( Params, TRIM(Name), GotVar ) 1140 IF( GotVar ) CALL AddToSaveList(TRIM(Name),Val) 1141 END DO 1142 1143 !------------------------------------------------------------------------------ 1144 ! Add results in Components 1145 !------------------------------------------------------------------------------ 1146 IF( ListGetLogical( Params,'Save Component Results',GotIt ) ) THEN 1147 CALL Info('SaveScalars','Saving results from component',Level=10) 1148 l = 0 1149 DO i = 1, Model % NumberOfComponents 1150 Lst => ListHead( Model % Components(i) % Values ) 1151 DO WHILE( ASSOCIATED( Lst ) ) 1152 IF ( Lst % Name(1:4) == TRIM(ResultPrefix) ) THEN 1153 IF ( ASSOCIATED(Lst % Fvalues) ) THEN 1154 CALL AddToSaveList('component '//TRIM(I2S(i))//': '//TRIM(Lst % Name), Lst % Fvalues(1,1,1)) 1155 l = l + 1 1156 END IF 1157 END IF 1158 Lst => Lst % Next 1159 END DO 1160 END DO 1161 IF(l > 0 ) THEN 1162 WRITE (Message,'(A,I0,A)') 'Found ',l,' result scalars in components section' 1163 CALL Info('SaveScalars',Message,Level=7) 1164 END IF 1165 END IF 1166 1167 1168 !------------------------------------------------------------------------------ 1169 ! If there are no values 1170 !------------------------------------------------------------------------------ 1171 IF( NoValues == 0 ) THEN 1172 CALL Warn('SaveScalars','Found no values to save') 1173 RETURN 1174 ELSE 1175 CALL Info('SaveScalars','Found '//TRIM(I2S(NoValues))//' values to save in total',Level=6) 1176 END IF 1177 1178 1179 1180 !------------------------------------------------------------------------------ 1181 ! Finally save all the scalars into a file 1182 !------------------------------------------------------------------------------ 1183 IF( SaveToFile ) THEN 1184 1185 LineInd = ListGetInteger( Params,'Line Marker',GotIt) 1186 PrevNoValues = ListGetInteger( Params,'Save Scalars Dofs',GotIt) 1187 1188 IF(WriteCore .AND. NoValues /= PrevNoValues) THEN 1189 CALL ListAddInteger( Params,'Save Scalars Dofs',NoValues ) 1190 1191 WRITE( Message, '(A)' ) 'Saving names of values to file: '//TRIM(ScalarNamesFile) 1192 CALL Info( 'SaveScalars', Message, Level=4 ) 1193 1194 IF( Solver % TimesVisited > 0 ) THEN 1195 WRITE ( Message,'(A,I0,A,I0)') 'Number of scalar values differ from previous time:',& 1196 NoValues,' vs. ',PrevNoValues 1197 CALL Warn('SaveScalars',Message) 1198 END IF 1199 1200 IF(ParallelWrite) CALL Info('SaveScalars','Parallel data is written into separate files',Level=6) 1201 IF(ParallelReduce) CALL Info('SaveScalars','Parallel data is reduced into one file',Level=6) 1202 IF(FileAppend) CALL Info('SaveScalars','Data is appended to existing file',Level=6) 1203 1204 OPEN(NEWUNIT=NamesUnit, FILE=ScalarNamesFile) 1205 1206 Message = ListGetString(Model % Simulation,'Comment',GotIt) 1207 IF( GotIt ) THEN 1208 WRITE(NamesUnit,'(A)') TRIM(Message) 1209 END IF 1210 1211 Message = ListGetString(Params,'Comment',GotIt) 1212 IF( GotIt ) THEN 1213 WRITE(NamesUnit,'(A)') TRIM(Message) 1214 END IF 1215 1216 DateStr = GetVersion() 1217 WRITE( NamesUnit,'(A)') 'Elmer version: '//TRIM(DateStr) 1218 1219 DateStr = GetRevision( GotIt ) 1220 IF( GotIt ) THEN 1221 WRITE( NamesUnit,'(A)') 'Elmer revision: '//TRIM(DateStr) 1222 END IF 1223 1224 DateStr = GetCompilationDate( GotIt ) 1225 IF( GotIt ) THEN 1226 WRITE( NamesUnit,'(A)') 'Elmer compilation date: '//TRIM(DateStr) 1227 END IF 1228 1229 DateStr = GetSifName( GotIt ) 1230 IF( GotIt ) THEN 1231 WRITE( NamesUnit,'(A)') 'Solver input file: '//TRIM(DateStr) 1232 END IF 1233 1234 DateStr = FormatDate() 1235 WRITE( NamesUnit,'(A)') 'File started at: '//TRIM(DateStr) 1236 1237 WRITE(NamesUnit,'(A)') ' ' 1238 WRITE(NamesUnit,'(A)') 'Variables in columns of matrix: '//TRIM(ScalarsFile) 1239 IF( LineInd /= 0 ) THEN 1240 i = 1 1241 WRITE(NamesUnit,'(I4,": ",A)') 1,'Line Marker' 1242 ELSE 1243 i = 0 1244 END IF 1245 DO No=1,NoValues 1246 WRITE(NamesUnit,'(I4,": ",A)') No+i,TRIM(ValueNames(No)) 1247 END DO 1248 CLOSE(NamesUnit) 1249 END IF 1250 1251 !------------------------------------------------------------------------------ 1252 ! In parallel case save the data in different files 1253 !------------------------------------------------------------------------------ 1254 1255 WRITE( Message,'(A)' ) 'Saving values to file: '// TRIM(ScalarsFile) 1256 CALL Info( 'SaveScalars', Message, Level=4 ) 1257 1258 IF ( ParallelWrite ) THEN 1259 WRITE( ScalarParFile, '(A,i0)' ) TRIM(ScalarsFile)//'.', ParEnv % MyPE 1260 1261 IF( Solver % TimesVisited > 0 .OR. FileAppend) THEN 1262 OPEN(NEWUNIT=ScalarsUnit, FILE=ScalarParFile,POSITION='APPEND') 1263 ELSE 1264 OPEN(NEWUNIT=ScalarsUnit, FILE=ScalarParFile) 1265 END IF 1266 ELSE IF( WriteCore ) THEN 1267 IF( Solver % TimesVisited > 0 .OR. FileAppend) THEN 1268 OPEN(NEWUNIT=ScalarsUnit, FILE=ScalarsFile,POSITION='APPEND') 1269 ELSE 1270 OPEN(NEWUNIT=ScalarsUnit, FILE=ScalarsFile) 1271 END IF 1272 END IF 1273 1274 1275 IF( WriteCore ) THEN 1276 ! If there are multiple lines it may be a good idea to mark each by an index 1277 IF( LineInd /= 0) THEN 1278 WRITE (ScalarsUnit,'(I6)',advance='no') LineInd 1279 END IF 1280 DO No=1,NoValues-1 1281 IF( ValuesInteger(No) ) THEN 1282 IntVal = NINT( Values(No) ) 1283 WRITE (ScalarsUnit,'(I10)',advance='no') IntVal 1284 ELSE 1285 WRITE (ScalarsUnit,'(ES22.12E3)',advance='no') Values(No) 1286 END IF 1287 END DO 1288 IF( ValuesInteger(NoValues) ) THEN 1289 IntVal = NINT( Values(No) ) 1290 WRITE (ScalarsUnit,'(I10)') IntVal 1291 ELSE 1292 WRITE (ScalarsUnit,'(ES22.12E3)') Values(NoValues) 1293 END IF 1294 CLOSE(ScalarsUnit) 1295 1296 !------------------------------------------------------------------------------ 1297 ! Save comments by line in a metadata file 1298 !------------------------------------------------------------------------------ 1299 IF( Solver % TimesVisited == 0 .AND. FileAppend .AND. LineInd /= 0 ) THEN 1300 Message = ListGetString(Params,'Comment',GotIt) 1301 Name = TRIM(ScalarsFile) // '.' // TRIM("marker") 1302 IF( GotIt ) THEN 1303 OPEN(NEWUNIT=ScalarsUnit, FILE=Name,POSITION='APPEND') 1304 WRITE(ScalarsUnit,'(I6,A,A)') LineInd,': ',TRIM(Message) 1305 CLOSE(ScalarsUnit) 1306 END IF 1307 END IF 1308 1309 !------------------------------------------------------------------------------ 1310 ! Save data in CVS format for, e.g. Livegraph 1311 !------------------------------------------------------------------------------ 1312 1313 IF(SaveCVS .AND. WriteCore ) THEN 1314 ! Save data as comma-separated-values (cvs-file) 1315 IF( Solver % TimesVisited > 0 .OR. FileAppend) THEN 1316 OPEN(NEWUNIT=ScalarsUnit, FILE=TRIM(ScalarsFile)//'.csv',POSITION='APPEND') 1317 ELSE 1318 OPEN(NEWUNIT=ScalarsUnit, FILE=TRIM(ScalarsFile)//'.csv') 1319 DO No=1,NoValues-1 1320 WRITE (ScalarsUnit,'(A)',advance='no') TRIM(ValueNames(No))//"," 1321 END DO 1322 WRITE (ScalarsUnit,'(A)') TRIM(ValueNames(NoValues)) 1323 END IF 1324 1325 DO No=1,NoValues-1 1326 WRITE (ScalarsUnit,'(ES22.12E3,A)',advance='no') Values(No),"," 1327 END DO 1328 WRITE (ScalarsUnit,'(ES22.12E3)') Values(NoValues) 1329 CLOSE(ScalarsUnit) 1330 END IF 1331 1332 IF(SaveJson .AND. WriteCore ) THEN 1333 CALL Info('SaveScalars','Saving data in jason format o file: '& 1334 //TRIM(ScalarsFile)//'.json',Level=6) 1335 1336 IF(ParallelWrite) CALL Warn('SaveScalars','Use "Parallel Reduce=True" with JSON format!') 1337 IF(FileAppend) CALL Info('SaveScalars','Data is appended to existing file',Level=6) 1338 1339 ! Save data in JSON format 1340 StartNewFile = ( Solver % TimesVisited == 0 .AND. .NOT. FileAppend) 1341 1342 IF( StartNewFile ) THEN 1343 OPEN(NEWUNIT=ScalarsUnit, ACCESS='stream',FORM='formatted',& 1344 FILE=TRIM(ScalarsFile)//'.json') 1345 1346 WRITE( ScalarsUnit, '(A)' ) '{' 1347 1348 DateStr = GetVersion() 1349 WRITE( ScalarsUnit,'(A)') ' "elmerver": "'//TRIM(DateStr)//'",' 1350 1351 DateStr = GetRevision( GotIt ) 1352 IF( GotIt ) THEN 1353 WRITE( ScalarsUnit,'(A)') ' "elmerrev": "'//TRIM(DateStr)//'",' 1354 END IF 1355 1356 DateStr = GetCompilationDate( GotIt ) 1357 IF( GotIt ) THEN 1358 WRITE( ScalarsUnit,'(A)') ' "elmercompiled": "'//TRIM(DateStr)//'",' 1359 END IF 1360 1361 DateStr = GetSifName( GotIt ) 1362 IF( GotIt ) THEN 1363 WRITE( ScalarsUnit,'(A)') ' "siffile": "'//TRIM(DateStr)//'",' 1364 END IF 1365 1366 DateStr = FormatDate() 1367 WRITE( ScalarsUnit,'(A)') ' "starttime": "'//TRIM(DateStr)//'",' 1368 1369 WRITE( ScalarsUnit,'(A)') ' "columns": '//TRIM(I2S(NoValues))//',' 1370 1371 WRITE( ScalarsUnit, '(A)',ADVANCE='no' ) ' "names":[' 1372 1373 DO No=1,NoValues 1374 WRITE(ScalarsUnit,'(A)',ADVANCE='no') '"'//TRIM(ValueNames(No))//'"' 1375 IF(No<NoValues) WRITE(ScalarsUnit,'(A)',ADVANCE='no') ',' 1376 END DO 1377 WRITE( ScalarsUnit, '(A)' ) '],' 1378 WRITE( ScalarsUnit, '(A)' ) ' "values":' 1379 WRITE( ScalarsUnit, '(A)' ) ' [' 1380 WRITE( ScalarsUnit, '(A)', ADVANCE='no' ) ' [' 1381 ELSE 1382 OPEN(NEWUNIT=ScalarsUnit, ACCESS='stream',FORM='formatted',& 1383 STATUS='old', POSITION='append',& 1384 FILE=TRIM(ScalarsFile)//'.json') 1385 1386 ! Jump to the start position from previous visit 1387 WRITE( ScalarsUnit,'(A)', POS=jsonpos, ADVANCE='no' ) ' ,[' 1388 END IF 1389 1390 DO No=1,NoValues-1 1391 WRITE (ScalarsUnit,'(ES22.12E3,A)',advance='no') Values(No),"," 1392 END DO 1393 WRITE (ScalarsUnit,'(ES22.12E3,A)') Values(No),"]" 1394 1395 ! Mark the start position for next visit 1396 INQUIRE(ScalarsUnit,POS=jsonpos) 1397 1398 WRITE( ScalarsUnit, '(A)' ) ' ]' 1399 WRITE( ScalarsUnit, '(A)' ) '}' 1400 1401 CLOSE(ScalarsUnit) 1402 END IF 1403 1404 1405 END IF 1406 END IF ! of SaveFile 1407 1408 !------------------------------------------------------------------------------ 1409 ! Echo values if requested, this is the default if no output to file 1410 !------------------------------------------------------------------------------ 1411 IF( EchoValues ) THEN 1412 CALL Info( 'SaveScalars','Showing computed results:') 1413 DO No=1,NoValues 1414 WRITE(Message,'(I4,": ",A,ES22.12E3)') No,TRIM(ValueNames(No)),Values(No) 1415 CALL Info('SaveScalars',Message) 1416 END DO 1417 END IF 1418 1419 !------------------------------------------------------------------------------ 1420 ! Add data to the Simulation block for possible future use 1421 !------------------------------------------------------------------------------ 1422 DO No=1,NoValues 1423 CALL ListAddConstReal(Model % Simulation,TRIM(ValueNames(No)),Values(No)) 1424 END DO 1425 1426 !------------------------------------------------------------------------------ 1427 ! For consistency checks one may print out a value imitating ComputeChange 1428 !------------------------------------------------------------------------------ 1429 NormInd = ListGetInteger( Params,'Show Norm Index',GotIt) 1430 IF( NormInd > 0 .AND. NormInd <= NoValues ) THEN 1431 Norm = Values( NormInd ) 1432 Solver % Variable % Values = Values( NormInd ) 1433 Solver % Variable % Norm = ABS( Values ( NormInd ) ) 1434 1435 Name = ListGetString( Params, 'Equation',GotIt) 1436 IF(.NOT. GotIt) Name = 'SaveScalars' 1437 1438 ! Here the name is ComputeChange in order to get the change also to ElmerGUI 1439 ! albeit in a very dirt style. One could also edit ElmerGUI.... 1440 WRITE( Message, '(a,g15.8,g15.8,a)') & 1441 'SS (ITER=1) (NRM,RELC): (',Norm, Change,& 1442 ' ) :: '//TRIM( Name ) 1443 CALL Info( 'ComputeChange', Message, Level=3 ) 1444 END IF 1445 1446 !------------------------------------------------------------------------------ 1447 ! One may also export desired data as a cost function for FindOptimum solver 1448 !------------------------------------------------------------------------------ 1449 CostInd = ListGetInteger( Params,'Cost Function Index',GotIt) 1450 IF( CostInd > 0 .AND. CostInd <= NoValues ) THEN 1451 CALL ListAddConstReal( Model % Simulation,'Cost Function',Values(CostInd) ) 1452 END IF 1453 1454 1455 IF( ParEnv % MyPe == 0 ) THEN 1456 IF( ListGetLogical( Params,'Mark Failed Strategy',GotIt) ) THEN 1457 LineInd = ListGetInteger( Params,'Line Marker',GotIt) 1458 IF(.NOT. GotIt) THEN 1459 CALL Fatal('SaveScalars','Failed strategy marked requires > Line Marker <') 1460 END IF 1461 1462 CALL Info('SaveScalars','Saving True marker at end') 1463 Name = 'FINISHED_MARKER_'//TRIM(I2S(LineInd)) 1464 i = 1 1465 OPEN(NEWUNIT=MarkerUnit,FILE=Name,STATUS='Unknown') 1466 WRITE(MarkerUnit,'(I0)') i 1467 CLOSE(MarkerUnit) 1468 END IF 1469 END IF 1470 1471 1472 IF( NoElements > 0 ) THEN 1473 DEALLOCATE( ElementNodes % x, ElementNodes % y, ElementNodes % z) 1474 END IF 1475 1476 n = 1 1477 n = NINT( ParallelReduction(1.0_dp*n) ) 1478 1479 CALL Info('SaveScalars', '-----------------------------------------', Level=7 ) 1480 1481 1482!------------------------------------------------------------------------------ 1483 1484CONTAINS 1485 1486 1487 FUNCTION OperToParOperMap( LocalOper ) RESULT ( ParOper ) 1488 1489 CHARACTER(LEN=MAX_NAME_LEN) :: LocalOper, ParOper 1490 1491 1492 ParOper = 'NONE' 1493 IF(.NOT. ParallelReduce ) RETURN 1494 1495 1496 SELECT CASE(LocalOper) 1497 1498 CASE('nodes','elements','dofs','sum','sum square','sum abs','int','int square','int abs','volume',& 1499 'potential energy', 'convective energy','diffusive energy','boundary sum','boundary dofs',& 1500 'boundary int','area','diffusive flux','convective flux','nans','partition checksum',& 1501 'partition neighbours checksum') 1502 ParOper = 'sum' 1503 1504 CASE('max','max abs','boundary max','boundary max abs') 1505 ParOper = 'max' 1506 1507 CASE('min','min abs','boundary min','boundary min abs') 1508 ParOper = 'min' 1509 1510 ! These operators should already be more of less parallel 1511 CASE('partitions','cpu time','wall time','cpu memory','norm','nonlin change','steady state change',& 1512 'nonlin iter','nonlin converged','steady converged','threads') 1513 ParOper = 'none' 1514 1515 CASE DEFAULT 1516 ParOper = ListGetString( Params,'Default Parallel Operator',GotIt) 1517 IF( .NOT. GotIt ) THEN 1518 ParOper = 'none' 1519 CALL Warn('SaveScalars','Reduction not implemented in parallel:'//TRIM(LocalOper)) 1520 END IF 1521 1522 END SELECT 1523 1524 CALL Info('OperToParOperMap',TRIM(LocalOper)//' -> '//TRIM(ParOper),Level=12) 1525 1526 END FUNCTION OperToParOperMap 1527 1528 1529 1530 1531 SUBROUTINE AddToSaveList(Name, Val, ValueIsInteger, ParallelOperator ) 1532 1533 INTEGER :: n 1534 CHARACTER(LEN=*) :: Name 1535 REAL(KIND=dp) :: Val, ParVal 1536 LOGICAL, OPTIONAL :: ValueIsInteger 1537 CHARACTER(LEN=MAX_NAME_LEN), OPTIONAL :: ParallelOperator 1538 !------------------------------------------------------------------------ 1539 CHARACTER(LEN=MAX_NAME_LEN) :: Str, ParOper 1540 REAL(KIND=dp), ALLOCATABLE :: TmpValues(:) 1541 CHARACTER(LEN=MAX_NAME_LEN), ALLOCATABLE :: TmpValueNames(:) 1542 LOGICAL, POINTER :: TmpValuesInteger(:) 1543 TYPE(Variable_t), POINTER :: TargetVar 1544 INTEGER :: MPIOper 1545 INTEGER :: istat 1546 LOGICAL :: GotParOper 1547 1548 SAVE TmpValues, TmpValueNames 1549 1550 ! For the first time allocate some space 1551 IF(.NOT. ALLOCATED(Values)) THEN 1552 n = 20 1553 ALLOCATE( Values(n), ValueNames(n), ValuesInteger(n), STAT=istat ) 1554 IF( istat /= 0 ) CALL Fatal('SaveScalars','Memory allocation error 6') 1555 Values = 0._dp 1556 END IF 1557 1558 n = NoValues 1559 1560 ! If vectors are too small make some more room in a rather dummy way 1561 IF(n >= SIZE(Values) ) THEN 1562 ALLOCATE(TmpValues(n), TmpValueNames(n), TmpValuesInteger(n),STAT=istat) 1563 IF( istat /= 0 ) CALL Fatal('SaveScalars','Memory allocation error 8') 1564 1565 TmpValues(1:n) = Values(1:n) 1566 TmpValuesInteger(1:n) = ValuesInteger(1:n) 1567 TmpValueNames(1:n) = ValueNames(1:n) 1568 DEALLOCATE(Values,ValueNames,ValuesInteger) 1569 1570 ALLOCATE(Values(n+10), ValueNames(n+10), ValuesInteger(n+10), STAT=istat ) 1571 IF( istat /= 0 ) CALL Fatal('SaveScalars','Memory allocation error 9') 1572 Values = 0._dp 1573 1574 Values(1:n) = TmpValues(1:n) 1575 ValuesInteger(1:n) = TmpValuesInteger(1:n) 1576 ValueNames(1:n) = TmpValueNames(1:n) 1577 DEALLOCATE(TmpValues,TmpValueNames,TmpValuesInteger) 1578 END IF 1579 1580 n = n + 1 1581 Values(n) = Val 1582 1583 ! If we have positive or negative operator then revert that back in the save name 1584 IF( PosOper ) THEN 1585 ValueNames(n) = 'positive '//TRIM(Name) 1586 ELSE IF( NegOper ) THEN 1587 ValueNames(n) = 'negative '//TRIM(Name) 1588 ELSE 1589 ValueNames(n) = TRIM(Name) 1590 END IF 1591 1592 IF( PRESENT( ValueIsInteger ) ) THEN 1593 ValuesInteger(n) = ValueIsInteger 1594 ELSE 1595 ValuesInteger(n) = .FALSE. 1596 END IF 1597 NoValues = n 1598 1599 IF( ParallelReduce ) THEN 1600 GotParOper = .FALSE. 1601 IF( PRESENT(ParallelOperator)) THEN 1602 ParOper = ParallelOperator 1603 GotParOper = .TRUE. 1604 ELSE 1605 1606!------------------------------------------------------------------------------------------------ 1607! This is to ensure that parallel operators may be applied to also other scalars than those 1608! computed within this solver with the serial operators. The indexing may not be known in advance 1609! but may be checked by a trial run. Note that the indexing does not often coinsice with the 1610! normal operators but then this loop is never reached. 1611!------------------------------------------------------------------------------------------------ 1612 1613 WRITE (Str,'(A,I0)') 'Parallel Operator ',n 1614 ParOper = ListGetString(Params,TRIM(Str),GotParOper) 1615 END IF 1616 1617 IF( GotParOper ) THEN 1618 SELECT CASE( ParOper ) 1619 1620 CASE('max' ) 1621 MPIOper = MPI_MAX 1622 CASE('min' ) 1623 MPIOper = MPI_MIN 1624 CASE('sum' ) 1625 MPIOper = MPI_SUM 1626 1627 CASE DEFAULT 1628 MPIOper = 0 1629 1630 END SELECT 1631 1632 1633 IF( MPIOper > 0 ) THEN 1634 CALL MPI_ALLREDUCE(Val,ParVal,1,& 1635 MPI_DOUBLE_PRECISION,MPIOper,ELMER_COMM_WORLD,ierr) 1636 Values(n) = ParVal 1637 IF( MPIOper == MPI_MIN ) THEN 1638 WRITE( ValueNames(n),'(A)') TRIM( ValueNames(n) )//' : mpi_min' 1639 ELSE IF( MPIOper == MPI_MAX ) THEN 1640 WRITE( ValueNames(n),'(A)') TRIM( ValueNames(n) )//' : mpi_max' 1641 ELSE IF( MPIOper == MPI_SUM ) THEN 1642 WRITE( ValueNames(n),'(A)') TRIM( ValueNames(n) )//' : mpi_sum' 1643 END IF 1644 END IF 1645 1646 END IF 1647 END IF 1648 1649 !------------------------------------------------------------------------------ 1650 ! If requested, create variable of the result 1651 ! This is performed already here so the variable can be used 1652 ! in subsequent definitions within SaveScalars. 1653 !------------------------------------------------------------------------------ 1654 WRITE (Str,'(A,I0)') 'Target Variable ',n 1655 VariableName = ListGetString( Params, TRIM(Str), GotIt ) 1656 IF( GotIt ) THEN 1657 TargetVar => VariableGet( Model % Variables, TRIM(VariableName) ) 1658 IF(.NOT. ASSOCIATED(TargetVar)) THEN 1659 NULLIFY(WrkPntr) 1660 ALLOCATE(WrkPntr(1),STAT=istat) 1661 IF( istat /= 0 ) CALL Fatal('SaveScalars','Memory allocation error 5') 1662 1663 CALL VariableAdd( Model % Variables, Mesh, Solver, & 1664 TRIM(VariableName), 1, WrkPntr ) 1665 TargetVar => VariableGet( Model % Variables, TRIM(VariableName) ) 1666 END IF 1667 TargetVar % Values(1) = Values(n) 1668 CALL Info('SaveScalars','Defining: '//TRIM(VariableName)//' = '//TRIM(ValueNames(n)),Level=8) 1669 END IF 1670 1671 END SUBROUTINE AddToSaveList 1672 1673 1674 ! Create table for masking nodes in statistical operators. 1675 !------------------------------------------------------------------------- 1676 SUBROUTINE CreateNodeMask( ) 1677 1678 INTEGER :: t 1679 TYPE(Element_t), POINTER :: Element 1680 TYPE(ValueList_t), POINTER :: MaskList 1681 1682 IF(.NOT. MaskOper ) RETURN 1683 1684 CALL Info('SaveScalars','Creating mask for: '//TRIM(MaskName),Level=10) 1685 1686 IF( ALLOCATED( NodeMask ) ) THEN 1687 IF( SIZE( NodeMask ) /= SIZE( Var % Perm ) ) DEALLOCATE( NodeMask ) 1688 END IF 1689 1690 IF(.NOT. ALLOCATED( NodeMask ) ) THEN 1691 ALLOCATE( NodeMask( SIZE( Var % Perm ) ) ) 1692 END IF 1693 NodeMask = .FALSE. 1694 1695 DO t = 1, Mesh % NumberOfBulkElements 1696 Element => Mesh % Elements(t) 1697 n = Element % TYPE % NumberOfNodes 1698 1699 IF( ElementalVar ) THEN 1700 NodeIndexes => Element % DgIndexes 1701 ELSE 1702 NodeIndexes => Element % NodeIndexes 1703 END IF 1704 1705 ! If we are masking operators with correct (body, body force, or material) then do it here 1706 IF( BodyOper ) THEN 1707 MaskList => GetBodyParams( Element ) 1708 IF(.NOT. ASSOCIATED(MaskList)) CYCLE 1709 ELSE IF( BodyForceOper ) THEN 1710 MaskList => GetBodyForce( Element, GotIt ) 1711 IF( .NOT. GotIt ) CYCLE 1712 ELSE IF( MaterialOper ) THEN 1713 MaskList => GetMaterial( Element, GotIt ) 1714 IF(.NOT. GotIt ) CYCLE 1715 ELSE 1716 CALL Fatal('MaskNodes','Unknown mask strategy') 1717 END IF 1718 IF( ListGetLogical( MaskList, MaskName, GotIt ) ) THEN 1719 NodeMask(NodeIndexes(1:n)) = .TRUE. 1720 END IF 1721 END DO 1722 1723 t = COUNT( NodeMask ) 1724 CALL Info('SaveScalars','Created mask of size: '& 1725 //TRIM(I2S(t)),Level=12) 1726 1727 END SUBROUTINE CreateNodeMask 1728 1729 1730 1731 1732 FUNCTION VectorStatistics(Var,OperName,NodalOper) RESULT (operx) 1733 1734 TYPE(Variable_t), POINTER :: Var 1735 CHARACTER(LEN=MAX_NAME_LEN) :: OperName 1736 LOGICAL, OPTIONAL :: NodalOper 1737 REAL(KIND=dp) :: operx 1738 1739 REAL(KIND=dp) :: Minimum, Maximum, AbsMinimum, AbsMaximum, & 1740 Mean, Variance, sumx, sumxx, sumabsx, x, Variance2 1741 INTEGER :: Nonodes, i, j, k, l, NoDofs, sumi 1742 TYPE(NeighbourList_t), POINTER :: nlist(:) 1743 INTEGER, POINTER :: PPerm(:) 1744 1745 CALL Info('SaveScalars','Computing operator: '//TRIM(OperName),Level=12) 1746 1747 sumi = 0 1748 sumx = 0.0 1749 sumxx = 0.0 1750 sumabsx = 0.0 1751 1752 Maximum = -HUGE(x) 1753 Minimum = HUGE(x) 1754 AbsMaximum = -HUGE(x) 1755 AbsMinimum = HUGE(x) 1756 1757 PPerm => Var % Perm 1758 IF( Var % TYPE == Variable_on_gauss_points .OR. & 1759 Var % TYPE == Variable_on_elements ) THEN 1760 NULLIFY( PPerm ) 1761 END IF 1762 1763 1764 NoDofs = Var % Dofs 1765 IF(ASSOCIATED (PPerm)) THEN 1766 Nonodes = SIZE(PPerm) 1767 ELSE 1768 Nonodes = SIZE(Var % Values) / NoDofs 1769 END IF 1770 1771 1772 IF( MaskOper ) THEN 1773 IF( NoNodes > SIZE(NodeMask) ) THEN 1774 CALL Info('SaveScalars','Decreasing operator range to size of mask: '& 1775 //TRIM(I2S(NoNodes))//' vs. '//TRIM(I2S(SIZE(NodeMask))), Level=8) 1776 NoNodes = SIZE(NodeMask) 1777 END IF 1778 END IF 1779 1780 nlist => NULL() 1781 IF( IsParallel ) THEN 1782 IF(ASSOCIATED(Var % Solver)) THEN 1783 IF(ASSOCIATED(Var % Solver % Matrix)) THEN 1784 IF(ASSOCIATED(Var % Solver % Matrix % ParallelInfo)) THEN 1785 nlist => Var % Solver % Matrix % ParallelInfo % NeighbourList 1786 END IF 1787 END IF 1788 END IF 1789 END IF 1790 1791 IF( PRESENT( NodalOper ) ) THEN 1792 IF( NodalOper ) THEN 1793 FirstInd = 1 1794 LastInd = Mesh % NumberOfNodes 1795 ELSE 1796 FirstInd = Mesh % NumberOfNodes + 1 1797 LastInd = SIZE( PPerm ) 1798 END IF 1799 ELSE 1800 FirstInd = 1 1801 LastInd = NoNodes 1802 END IF 1803 1804 1805 DO i=FirstInd,LastInd 1806 IF( MaskOper ) THEN 1807 IF( .NOT. NodeMask(i) ) CYCLE 1808 END IF 1809 1810 j = i 1811 IF(ASSOCIATED(PPerm)) j = Var % Perm(i) 1812 1813 IF(j > 0) THEN 1814 IF( IsParallel .AND. ASSOCIATED(nlist) ) THEN 1815 IF( nlist(j) % Neighbours(1) /= ParEnv % MyPE ) CYCLE 1816 END IF 1817 1818 IF(NoDofs <= 1) THEN 1819 x = Var % Values(j) 1820 ELSE 1821 x = 0.0d0 1822 DO l=1,NoDofs 1823 x = x + Var % Values(NoDofs*(j-1)+l) ** 2 1824 END DO 1825 x = SQRT(x) 1826 END IF 1827 1828 IF( PosOper .AND. x < 0 ) CYCLE 1829 IF( NegOper .AND. x > 0 ) CYCLE 1830 1831 sumi = sumi + 1 1832 sumx = sumx + x 1833 sumxx = sumxx + x*x 1834 sumabsx = sumabsx + ABS( x ) 1835 1836 Maximum = MAX(x,Maximum) 1837 Minimum = MIN(x,Minimum) 1838 IF(ABS(x) > ABS(AbsMaximum) ) AbsMaximum = x 1839 IF(ABS(x) < ABS(AbsMinimum) ) AbsMinimum = x 1840 END IF 1841 END DO 1842 1843 ! If there are no dofs avoid division by zero 1844 IF(sumi == 0) THEN 1845 operx = 0.0d0 1846 RETURN 1847 END IF 1848 1849 Mean = sumx / sumi 1850 1851 Variance2 = sumxx/sumi-Mean*Mean 1852 IF(Variance2 > 0.0d0) THEN 1853 Variance = SQRT(Variance2) 1854 ELSE 1855 Variance = 0.0d0 1856 END IF 1857 1858 SELECT CASE(OperName) 1859 1860 CASE ('sum') 1861 operx = sumx 1862 1863 CASE ('sum square') 1864 operx = sumxx 1865 1866 CASE ('sum abs') 1867 operx = sumabsx 1868 1869 CASE ('max') 1870 operx = Maximum 1871 1872 CASE ('min') 1873 operx = Minimum 1874 1875 CASE ('range') 1876 operx = Maximum - Minimum 1877 1878 CASE ('max abs') 1879 operx = AbsMaximum 1880 1881 CASE ('min abs') 1882 operx = AbsMinimum 1883 1884 CASE ('mean') 1885 operx = Mean 1886 1887 CASE ('mean square') 1888 operx = sumxx / sumi 1889 1890 CASE ('rms') 1891 operx = SQRT( sumxx / sumi ) 1892 1893 CASE ('mean abs') 1894 operx = sumabsx / sumi 1895 1896 CASE ('variance') 1897 operx = Variance 1898 1899 CASE DEFAULT 1900 CALL Warn('SaveScalars','Unknown statistical operator!') 1901 1902 END SELECT 1903 1904 1905 CALL Info('SaveScalars','Finished computing operator',Level=12) 1906 1907 1908 END FUNCTION VectorStatistics 1909 1910 1911 1912 1913!------------------------------------------------------------------------------ 1914 1915 FUNCTION VectorMeanDeviation(Var,OperName) RESULT (Deviation) 1916 1917 TYPE(Variable_t), POINTER :: Var 1918 CHARACTER(LEN=MAX_NAME_LEN) :: OperName 1919 REAL(KIND=dp) :: Mean, Deviation 1920 REAL(KIND=dp) :: sumx, sumdx, x, dx 1921 INTEGER :: Nonodes, i, j, k, NoDofs, sumi 1922 INTEGER, POINTER :: PPerm(:) 1923 1924 NoDofs = Var % Dofs 1925 1926 PPerm => Var % Perm 1927 IF( Var % TYPE == Variable_on_gauss_points .OR. & 1928 Var % TYPE == Variable_on_elements ) THEN 1929 NULLIFY( PPerm ) 1930 END IF 1931 1932 IF(ASSOCIATED (PPerm)) THEN 1933 Nonodes = SIZE(PPerm) 1934 ELSE 1935 Nonodes = SIZE(Var % Values) / NoDofs 1936 END IF 1937 1938 IF( MaskOper ) THEN 1939 IF( NoNodes > SIZE(NodeMask) ) THEN 1940 CALL Info('SaveScalars','Decreasing operator range to size of mask: '& 1941 //TRIM(I2S(NoNodes))//' vs. '//TRIM(I2S(SIZE(NodeMask))) ) 1942 NoNodes = SIZE(NodeMask) 1943 END IF 1944 END IF 1945 1946 sumi = 0 1947 sumx = 0.0 1948 Deviation = 0.0 1949 1950 DO i=1,Nonodes 1951 IF( MaskOper ) THEN 1952 IF( .NOT. NodeMask(i) ) CYCLE 1953 END IF 1954 1955 j = i 1956 1957 IF( IsParallel ) THEN 1958 IF( Mesh % ParallelInfo % NeighbourList(j) % Neighbours(1) /= ParEnv % MyPE ) CYCLE 1959 END IF 1960 1961 IF(ASSOCIATED(PPerm)) j = PPerm(i) 1962 IF(j > 0) THEN 1963 IF(NoDofs <= 1) THEN 1964 x = Var % Values(j) 1965 ELSE 1966 x = 0.0d0 1967 DO k=1,NoDofs 1968 x = x + Var % Values(NoDofs*(j-1)+k) ** 2 1969 END DO 1970 x = SQRT(x) 1971 END IF 1972 1973 IF( PosOper .AND. x < 0 ) CYCLE 1974 IF( NegOper .AND. x > 0 ) CYCLE 1975 1976 sumi = sumi + 1 1977 sumx = sumx + x 1978 END IF 1979 END DO 1980 1981 IF(sumi == 0) RETURN 1982 1983 Mean = sumx / sumi 1984 sumi = 0 1985 sumdx = 0.0 1986 DO i=1,Nonodes 1987 j = i 1988 IF(ASSOCIATED(PPerm)) j = PPerm(i) 1989 IF(j > 0) THEN 1990 IF(NoDofs <= 1) THEN 1991 x = Var % Values(j) 1992 ELSE 1993 x = 0.0d0 1994 DO k=1,NoDofs 1995 x = x + Var % Values(NoDofs*(j-1)+k) ** 2 1996 END DO 1997 x = SQRT(x) 1998 END IF 1999 2000 IF( PosOper .AND. x < 0 ) CYCLE 2001 IF( NegOper .AND. x > 0 ) CYCLE 2002 2003 dx = ABS(x-Mean) 2004 sumi = sumi + 1 2005 sumdx = sumdx + dx 2006 END IF 2007 END DO 2008 Deviation = sumdx / sumi 2009 2010 END FUNCTION VectorMeanDeviation 2011 2012!------------------------------------------------------------------------------ 2013 2014 FUNCTION BulkIntegrals(Var, OperName, GotCoeff, CoeffName) RESULT (operx) 2015 TYPE(Variable_t), POINTER :: Var 2016 CHARACTER(LEN=MAX_NAME_LEN) :: OperName, CoeffName 2017 LOGICAL :: GotCoeff 2018 REAL(KIND=dp) :: operx, vol 2019 2020 INTEGER :: t, hits 2021 TYPE(Element_t), POINTER :: Element 2022 INTEGER, POINTER :: NodeIndexes(:), PermIndexes(:) 2023 REAL(KIND=dp) :: SqrtMetric,Metric(3,3),Symb(3,3,3),dSymb(3,3,3,3) 2024 REAL(KIND=dp) :: Basis(Model % MaxElementNodes), dBasisdx(Model % MaxElementNodes,3) 2025 REAL(KIND=dp) :: EnergyTensor(3,3,Model % MaxElementNodes),& 2026 EnergyCoeff(Model % MaxElementNodes), ElemVals(Model % MaxElementNodes) 2027 REAL(KIND=dp) :: SqrtElementMetric,U,V,W,S,A,L,C(3,3),x,y,z 2028 REAL(KIND=dp) :: func, coeff, integral1, integral2, Grad(3), CoeffGrad(3) 2029 REAL(KIND=DP), POINTER :: Pwrk(:,:,:) 2030 LOGICAL :: Stat 2031 TYPE(ValueList_t), POINTER :: MaskList 2032 2033 INTEGER :: i,j,k,p,q,DIM,NoDofs 2034 2035 TYPE(GaussIntegrationPoints_t) :: IntegStuff 2036 2037 hits = 0 2038 integral1 = 0._dp 2039 integral2 = 0._dp 2040 vol = 0._dp 2041 EnergyCoeff = 1.0d0 2042 2043 NoDofs = Var % Dofs 2044 2045 DIM = CoordinateSystemDimension() 2046 2047 2048 DO t = 1, Mesh % NumberOfBulkElements + Mesh % NumberOfBoundaryElements 2049 2050 IF(t == Mesh % NumberOfBulkElements + 1 .AND. hits > 0) GOTO 10 2051 2052 Element => Mesh % Elements(t) 2053 Model % CurrentElement => Mesh % Elements(t) 2054 n = Element % TYPE % NumberOfNodes 2055 2056 IF ( Element % TYPE % ElementCode == 101 ) CYCLE 2057 2058 IF( ElementalVar ) THEN 2059 PermIndexes => Element % DgIndexes 2060 ELSE 2061 PermIndexes => Element % NodeIndexes 2062 END IF 2063 2064 IF ( ANY(Var % Perm(PermIndexes) == 0 ) ) CYCLE 2065 hits = hits + 1 2066 2067 2068 NodeIndexes => Element % NodeIndexes 2069 ElementNodes % x(1:n) = Mesh % Nodes % x(NodeIndexes(1:n)) 2070 ElementNodes % y(1:n) = Mesh % Nodes % y(NodeIndexes(1:n)) 2071 ElementNodes % z(1:n) = Mesh % Nodes % z(NodeIndexes(1:n)) 2072 2073 2074 ! If we are masking operators with correct (body, body force, or material) then do it here 2075 IF( BodyOper ) THEN 2076 MaskList => GetBodyParams( Element ) 2077 ELSE IF( BodyForceOper ) THEN 2078 MaskList => GetBodyForce( Element, GotIt ) 2079 IF( .NOT. GotIt ) CYCLE 2080 ELSE IF( MaterialOper ) THEN 2081 MaskList => GetMaterial( Element, GotIt ) 2082 IF(.NOT. GotIt ) CYCLE 2083 END IF 2084 IF( MaskOper ) THEN 2085 IF( .NOT. ListGetLogical( MaskList, MaskName, GotIt ) ) CYCLE 2086 END IF 2087 2088 IF( NoDOFs == 1 ) THEN 2089 ElemVals(1:n) = Var % Values(Var % Perm(PermIndexes) ) 2090 ELSE 2091 ElemVals(1:n) = 0.0_dp 2092 DO i=1,NoDOFs 2093 ElemVals(1:n) = ElemVals(1:n) + Var % Values(NoDofs*(Var % Perm(PermIndexes)-1)+i )**2 2094 END DO 2095 ElemVals(1:) = SQRT(ElemVals(1:n)) 2096 END IF 2097 2098 2099 k = ListGetInteger( Model % Bodies( Element % BodyId ) % Values, & 2100 'Material', GotIt, minv=1, maxv=Model % NumberOfMaterials ) 2101 2102 IF( OperName == 'diffusive energy' ) THEN 2103 EnergyTensor = 0.0d0 2104 IF(GotCoeff) THEN 2105 CALL ListGetRealArray( Model % Materials(k) % Values, & 2106 TRIM(CoeffName), Pwrk, n, NodeIndexes ) 2107 2108 IF ( SIZE(Pwrk,1) == 1 ) THEN 2109 DO i=1,3 2110 EnergyTensor( i,i,1:n ) = Pwrk( 1,1,1:n ) 2111 END DO 2112 ELSE IF ( SIZE(Pwrk,2) == 1 ) THEN 2113 DO i=1,MIN(3,SIZE(Pwrk,1)) 2114 EnergyTensor(i,i,1:n) = Pwrk(i,1,1:n) 2115 END DO 2116 ELSE 2117 DO i=1,MIN(3,SIZE(Pwrk,1)) 2118 DO j=1,MIN(3,SIZE(Pwrk,2)) 2119 EnergyTensor( i,j,1:n ) = Pwrk(i,j,1:n) 2120 END DO 2121 END DO 2122 END IF 2123 ELSE 2124 DO i=1,3 2125 EnergyTensor(i,i,1:n) = 1.0d0 2126 END DO 2127 END IF 2128 ELSE 2129 k = ListGetInteger( Model % Bodies( Element % BodyId ) % Values, & 2130 'Material', GotIt, minv=1, maxv=Model % NumberOfMaterials ) 2131 IF(GotCoeff) THEN 2132 EnergyCoeff = ListGetReal( Model % Materials(k) % Values, & 2133 TRIM(CoeffName), n, NodeIndexes(1:n) ) 2134 ELSE 2135 EnergyCoeff(1:n) = 1.0d0 2136 END IF 2137 END IF 2138 2139!------------------------------------------------------------------------------ 2140! Numerical integration 2141!------------------------------------------------------------------------------ 2142 IntegStuff = GaussPoints( Element ) 2143 2144 DO i=1,IntegStuff % n 2145 U = IntegStuff % u(i) 2146 V = IntegStuff % v(i) 2147 W = IntegStuff % w(i) 2148!------------------------------------------------------------------------------ 2149! Basis function values & derivatives at the integration point 2150!------------------------------------------------------------------------------ 2151 stat = ElementInfo( Element,ElementNodes,U,V,W,SqrtElementMetric, & 2152 Basis,dBasisdx ) 2153!------------------------------------------------------------------------------ 2154! Coordinatesystem dependent info 2155!------------------------------------------------------------------------------ 2156 s = 1.0 2157 IF ( CurrentCoordinateSystem() /= Cartesian ) THEN 2158 x = SUM( ElementNodes % x(1:n)*Basis(1:n) ) 2159 y = SUM( ElementNodes % y(1:n)*Basis(1:n) ) 2160 z = SUM( ElementNodes % z(1:n)*Basis(1:n) ) 2161 s = 2*PI 2162 END IF 2163 2164 CALL CoordinateSystemInfo( Metric,SqrtMetric,Symb,dSymb,x,y,z ) 2165 2166 s = s * SqrtMetric * SqrtElementMetric * IntegStuff % s(i) 2167 coeff = SUM( EnergyCoeff(1:n) * Basis(1:n)) 2168 vol = coeff * vol + S 2169 2170 SELECT CASE(OperName) 2171 2172 CASE ('volume') 2173 integral1 = integral1 + coeff * S 2174 2175 CASE ('int','int mean') 2176 func = SUM( ElemVals(1:n) * Basis(1:n) ) 2177 IF( PosOper ) func = MAX( 0.0_dp, func ) 2178 IF( NegOper ) func = MIN( 0.0_dp, func ) 2179 2180 integral1 = integral1 + S * coeff * func 2181 2182 CASE ('int square','int square mean','int rms') 2183 func = SUM( ElemVals(1:n) * Basis(1:n) ) 2184 IF( PosOper ) func = MAX( 0.0_dp, func ) 2185 IF( NegOper ) func = MIN( 0.0_dp, func ) 2186 integral1 = integral1 + S * coeff * func**2 2187 2188 CASE ('int abs','int abs mean') 2189 func = ABS( SUM( ElemVals(1:n) * Basis(1:n) ) ) 2190 IF( PosOper ) func = MAX( 0.0_dp, func ) 2191 IF( NegOper ) func = MIN( 0.0_dp, func ) 2192 integral1 = integral1 + S * coeff * func 2193 2194 CASE ('int variance') 2195 func = SUM( ElemVals(1:n) * Basis(1:n) ) 2196 IF( PosOper ) func = MAX( 0.0_dp, func ) 2197 IF( NegOper ) func = MIN( 0.0_dp, func ) 2198 integral1 = integral1 + S * coeff * func 2199 integral2 = integral2 + S * coeff * func**2 2200 2201 CASE ('diffusive energy') 2202 CoeffGrad = 0.0d0 2203 DO j = 1, DIM 2204 Grad(j) = SUM( dBasisdx(1:n,j) * ElemVals(1:n) ) 2205 DO k = 1, DIM 2206 CoeffGrad(j) = CoeffGrad(j) + SUM( EnergyTensor(j,k,1:n) * Basis(1:n) ) * & 2207 SUM( dBasisdx(1:n,k) * Var % Values(Var % Perm(PermIndexes)) ) 2208 END DO 2209 END DO 2210 2211 integral1 = integral1 + s * SUM( Grad(1:DIM) * CoeffGrad(1:DIM) ) 2212 2213 CASE ('convective energy') 2214 func = SUM( ElemVals(1:n) * Basis(1:n) ) 2215 IF( PosOper ) func = MAX( 0.0_dp, func ) 2216 IF( NegOper ) func = MIN( 0.0_dp, func ) 2217 2218 IF(NoDofs == 1) THEN 2219 func = SUM( ElemVals(1:n) * Basis(1:n) ) 2220 integral1 = integral1 + s * coeff * func**2 2221 ELSE 2222 func = 0.0d0 2223 DO j=1,MIN(DIM,NoDofs) 2224 func = SUM( Var % Values(NoDofs*(Var % Perm(PermIndexes)-1)+j) * Basis(1:n) ) 2225 integral1 = integral1 + s * coeff * func**2 2226 END DO 2227 END IF 2228 2229 CASE ('potential energy') 2230 2231 func = SUM( ElemVals(1:n) * Basis(1:n) ) 2232 IF( PosOper ) func = MAX( 0.0_dp, func ) 2233 IF( NegOper ) func = MIN( 0.0_dp, func ) 2234 2235 integral1 = integral1 + s * coeff * func 2236 2237 CASE DEFAULT 2238 CALL Warn('SaveScalars','Unknown statistical OPERATOR') 2239 2240 END SELECT 2241 2242 END DO 2243 2244 END DO 2245 2246 224710 CONTINUE 2248 2249 operx = 0.0d0 2250 IF(hits == 0) RETURN 2251 2252 SELECT CASE(OperName) 2253 2254 CASE ('volume','int','int abs','int square') 2255 operx = integral1 2256 2257 CASE ('int mean','int square mean','int abs mean') 2258 operx = integral1 / vol 2259 2260 CASE ('int rms') 2261 operx = SQRT( integral1 / vol ) 2262 2263 CASE ('int variance') 2264 operx = SQRT(integral2/vol-(integral1/vol)**2) 2265 2266 CASE ('diffusive energy','convective energy') 2267 operx = 0.5d0 * integral1 2268 2269 CASE ('potential energy') 2270 operx = integral1 2271 2272 END SELECT 2273 2274 2275 END FUNCTION BulkIntegrals 2276!------------------------------------------------------------------------------ 2277 2278 2279 SUBROUTINE BoundaryIntegrals(Var, OperName, GotCoeff, & 2280 CoeffName, fluxes, areas, fluxescomputed) 2281 2282 TYPE(Variable_t), POINTER :: Var 2283 CHARACTER(LEN=MAX_NAME_LEN) :: OperName, CoeffName 2284 LOGICAL :: GotCoeff 2285 REAL(KIND=dp) :: fluxes(:), areas(:) 2286 INTEGER :: fluxescomputed(:) 2287 2288 TYPE(Variable_t), POINTER :: VeloVar 2289 INTEGER :: t, FluxBody, LBody, RBody, NActive 2290 TYPE(Element_t), POINTER :: Element, Parent 2291 TYPE(ValueList_t), POINTER :: Material, BCVal 2292 REAL(KIND=dp) :: SqrtMetric,Metric(3,3),Symb(3,3,3),dSymb(3,3,3,3) 2293 REAL(KIND=dp) :: Basis(Model % MaxElementNodes),dBasisdx(Model % MaxElementNodes,3),& 2294 ParentBasis(Model % MaxElementNodes),& 2295 ParentdBasisdx(Model % MaxElementNodes,3),EnergyTensor(3,3,Model % MaxElementNodes),& 2296 EnergyCoeff(Model % MaxElementNodes) 2297 REAL(KIND=dp) :: SqrtElementMetric,U,V,W,up,vp,wp,S,A,L,C(3,3),x,y,z 2298 REAL(KIND=dp) :: func, coeff, Normal(3), Flow(3), flux 2299 REAL(KIND=DP), POINTER :: Pwrk(:,:,:) => Null() 2300 INTEGER, POINTER :: ParentIndexes(:), PermIndexes(:) 2301 REAL(KIND=dp) :: LocalVectorSolution(3,35), LocalVeloSolution(3,35) 2302 2303 LOGICAL :: Stat, Permutated 2304 INTEGER :: i,j,k,p,q,DIM,bc,NoDofs,pn,hits,istat 2305 TYPE(GaussIntegrationPoints_t) :: IntegStuff 2306 TYPE(Nodes_t) :: ParentNodes 2307 2308 n = Model % MaxElementNodes 2309 2310 DIM = CoordinateSystemDimension() 2311 2312 ALLOCATE(ParentNodes % x(n), ParentNodes % y(n), ParentNodes % z(n), PermIndexes(n), STAT=istat ) 2313 IF( istat /= 0 ) CALL Fatal('BoundaryIntegrals','Memory allocation error') 2314 2315 2316 IF( SaveFluxRange ) THEN 2317 Minimum = HUGE(minimum) 2318 Maximum = -HUGE(maximum) 2319 END IF 2320 2321 IF( OperName == 'area' ) THEN 2322 NoDofs = 1 2323 Permutated = .FALSE. 2324 ELSE 2325 NoDofs = Var % Dofs 2326 Permutated = ASSOCIATED(Var % Perm) 2327 END IF 2328 2329 DO i=1,3 2330 EnergyTensor(i,i,:) = 1.0d0 2331 END DO 2332 2333 SELECT CASE(OperName) 2334 2335 CASE('diffusive flux') 2336 IF(NoDofs /= 1) THEN 2337 CALL Warn('SaveScalars','diffusive flux & NoDofs /= 1?') 2338 RETURN 2339 END IF 2340 2341 CASE ('convective flux') 2342 IF( NoDofs < DIM) THEN 2343 IF( NoDofs == 1 ) THEN 2344 VeloVar => VariableGet( Mesh % Variables,'Flow Solution') 2345 IF( .NOT. ASSOCIATED( VeloVar ) ) THEN 2346 CALL Fatal('SaveScalars','For convective flux of a scalar we need the convective field!') 2347 END IF 2348 ELSE 2349 CALL Warn('SaveScalars','convective flux & NoDofs < DIM?') 2350 RETURN 2351 END IF 2352 END IF 2353 2354 CASE ('area','boundary int','boundary int mean') 2355 2356 CASE DEFAULT 2357 CALL Warn('SaveScalars','Unknown statistical OPERATOR') 2358 2359 END SELECT 2360 2361 fluxes = 0.0_dp 2362 areas = 0.0_dp 2363 coeff = 1.0_dp 2364 2365 2366 NActive = GetNOFBoundaryElements() 2367 2368 DO t = 1, NActive 2369 2370 Element => GetBoundaryElement(t, Var % Solver) 2371 BCVal => GetBC() 2372 IF (.NOT. ASSOCIATED(BCVal) ) CYCLE 2373 2374 IF (NoDOFs > 1) THEN 2375 CALL GetVectorLocalSolution(LocalVectorSolution, & 2376 UElement=Element, USolver=Var % Solver, UVariable=Var) 2377 ELSE 2378 CALL GetScalarLocalSolution(LocalVectorSolution(1,:), & 2379 UElement=Element, USolver=Var % Solver, UVariable=Var) 2380 IF( ComponentVar ) THEN 2381 IF( ASSOCIATED( Var2 ) ) THEN 2382 CALL GetScalarLocalSolution(LocalVectorSolution(2,:), & 2383 UElement=Element, USolver=Var % Solver, UVariable=Var2) 2384 NoDofs = 2 2385 ELSE 2386 LocalVectorSolution(2,:) = 0.0_dp 2387 END IF 2388 IF( ASSOCIATED( Var3 ) ) THEN 2389 CALL GetScalarLocalSolution(LocalVectorSolution(3,:), & 2390 UElement=Element, USolver=Var % Solver, UVariable=Var3) 2391 NoDofs = 3 2392 ELSE 2393 LocalVectorSolution(3,:) = 0.0_dp 2394 END IF 2395 END IF 2396 END IF 2397 2398 Model % CurrentElement => Element 2399 2400 IF ( Element % TYPE % ElementCode == 101 ) CYCLE 2401 2402 n = Element % TYPE % NumberOfNodes 2403 NodeIndexes => Element % NodeIndexes 2404 2405 IF(Permutated) THEN 2406 PermIndexes(1:n) = Var % Perm(NodeIndexes(1:n)) 2407 IF (ANY( PermIndexes(1:n) == 0)) CYCLE 2408 ELSE 2409 PermIndexes(1:n) = NodeIndexes(1:n) 2410 END IF 2411 2412 2413 DO bc=1, Model % NumberOfBCs 2414 2415 IF ( Model % BCs(bc) % Tag /= Element % BoundaryInfo % Constraint ) CYCLE 2416 2417 IF(.NOT. ListGetLogical(Model % BCs(bc) % Values,'Flux Integrate',gotIt ) .AND. & 2418 .NOT. ListGetLogical(Model % BCs(bc) % Values, MaskName, gotIt ) ) CYCLE 2419 2420 ElementNodes % x(1:n) = Mesh % Nodes % x(NodeIndexes(1:n)) 2421 ElementNodes % y(1:n) = Mesh % Nodes % y(NodeIndexes(1:n)) 2422 ElementNodes % z(1:n) = Mesh % Nodes % z(NodeIndexes(1:n)) 2423 2424 2425 2426 SELECT CASE(OperName) 2427 2428 CASE('diffusive flux') 2429 2430 FluxBody = ListGetInteger( Model % BCs(bc) % Values, & 2431 'Flux Integrate Body', gotIt ) 2432 IF ( GotIt ) THEN 2433 IF ( ASSOCIATED( Element % BoundaryInfo % Left ) ) & 2434 Lbody = Element % BoundaryInfo % Left % BodyId 2435 IF ( ASSOCIATED( Element % BoundaryInfo % Right ) ) & 2436 Rbody = Element % BoundaryInfo % Right % BodyId 2437 2438 IF ( Lbody == FluxBody ) THEN 2439 Parent => Element % BoundaryInfo % Left 2440 ELSEIF ( Rbody == FluxBody ) THEN 2441 Parent => Element % BoundaryInfo % Right 2442 ELSE 2443 WRITE( Message, * ) 'No such flux integrate body on bc ', & 2444 Element % BoundaryInfo % Constraint 2445 CALL Fatal( 'SaveScalars', Message ) 2446 END IF 2447 ELSE 2448 Parent => ELement % BoundaryInfo % Left 2449 stat = ASSOCIATED( Parent ) 2450 2451 IF(Permutated) THEN 2452 IF(stat) stat = ALL(Var % Perm(Parent % NodeIndexes) > 0) 2453 2454 IF ( .NOT. stat ) THEN 2455 Parent => ELement % BoundaryInfo % Right 2456 stat = ASSOCIATED( Parent ) 2457 IF(stat) stat = ALL(Var % Perm(Parent % NodeIndexes) > 0) 2458 END IF 2459 END IF 2460 IF ( .NOT. stat ) CALL Fatal( 'SaveScalars',& 2461 'No solution available for specified boundary' ) 2462 END IF 2463 2464 pn = Parent % TYPE % NumberOfNodes 2465 ParentIndexes => Parent % NodeIndexes 2466 i = ListGetInteger( Model % Bodies(Parent % BodyId) % Values, 'Material', & 2467 minv=1, maxv=Model % NumberOFMaterials ) 2468 Material => Model % Materials(i) % Values 2469 2470 ParentNodes % x(1:pn) = Mesh % Nodes % x(ParentIndexes(1:pn)) 2471 ParentNodes % y(1:pn) = Mesh % Nodes % y(ParentIndexes(1:pn)) 2472 ParentNodes % z(1:pn) = Mesh % Nodes % z(ParentIndexes(1:pn)) 2473 2474 EnergyTensor = 0._dp 2475 2476 IF(GotCoeff) THEN 2477 CALL ListGetRealArray( Material, CoeffName, Pwrk, pn, ParentIndexes ) 2478 2479 IF ( SIZE(Pwrk,1) == 1 ) THEN 2480 DO i=1,3 2481 EnergyTensor( i,i,1:pn ) = Pwrk( 1,1,1:pn ) 2482 END DO 2483 ELSE IF ( SIZE(Pwrk,2) == 1 ) THEN 2484 DO i=1,MIN(3,SIZE(Pwrk,1)) 2485 EnergyTensor(i,i,1:pn) = Pwrk(i,1,1:pn) 2486 END DO 2487 ELSE 2488 DO i=1,MIN(3,SIZE(Pwrk,1)) 2489 DO j=1,MIN(3,SIZE(Pwrk,2)) 2490 EnergyTensor( i,j,1:pn ) = Pwrk(i,j,1:pn) 2491 END DO 2492 END DO 2493 END IF 2494 ELSE 2495 DO i=1,3 2496 EnergyTensor(i,i,1:pn) = 1.0d0 2497 END DO 2498 END IF 2499 EnergyCoeff(1:n) = 1.0D00 2500 fluxescomputed(bc) = fluxescomputed(bc) + 1 2501 2502 2503 CASE ('convective flux') 2504 2505 FluxBody = ListGetInteger( Model % BCs(bc) % Values, & 2506 'Flux Integrate Body', gotIt ) 2507 IF ( GotIt ) THEN 2508 IF ( ASSOCIATED( Element % BoundaryInfo % Left ) ) & 2509 Lbody = Element % BoundaryInfo % Left % BodyId 2510 IF ( ASSOCIATED( Element % BoundaryInfo % Right ) ) & 2511 Rbody = Element % BoundaryInfo % Right % BodyId 2512 2513 IF ( Lbody == FluxBody ) THEN 2514 Parent => Element % BoundaryInfo % Left 2515 ELSEIF ( Rbody == FluxBody ) THEN 2516 Parent => Element % BoundaryInfo % Right 2517 ELSE 2518 WRITE( Message, * ) 'No such flux integrate body on bc ', & 2519 Element % BoundaryInfo % Constraint 2520 CALL Fatal( 'SaveScalars', Message ) 2521 END IF 2522 ELSE 2523 Parent => ELement % BoundaryInfo % Left 2524 stat = ASSOCIATED( Parent ) 2525 2526 IF(Permutated) THEN 2527 IF(stat) stat = ALL(Var % Perm(Parent % NodeIndexes) > 0) 2528 2529 IF ( .NOT. stat ) THEN 2530 Parent => ELement % BoundaryInfo % Right 2531 stat = ASSOCIATED( Parent ) 2532 IF(stat) stat = ALL(Var % Perm(Parent % NodeIndexes) > 0) 2533 END IF 2534 END IF 2535 IF ( .NOT. stat ) CALL Fatal( 'SaveScalars',& 2536 'No solution available for specified boundary' ) 2537 END IF 2538 2539 IF( NoDofs == 1 ) THEN 2540 CALL GetVectorLocalSolution(LocalVeloSolution, & 2541 UElement=Element, USolver=VeloVar % Solver, UVariable=VeloVar) 2542 2543 END IF 2544 2545 pn = Parent % TYPE % NumberOfNodes 2546 ParentIndexes => Parent % NodeIndexes 2547 i = ListGetInteger( Model % Bodies(Parent % BodyId) % Values, 'Material', & 2548 minv=1, maxv=Model % NumberOFMaterials ) 2549 Material => Model % Materials(i) % Values 2550 2551 ParentNodes % x(1:pn) = Mesh % Nodes % x(ParentIndexes(1:pn)) 2552 ParentNodes % y(1:pn) = Mesh % Nodes % y(ParentIndexes(1:pn)) 2553 ParentNodes % z(1:pn) = Mesh % Nodes % z(ParentIndexes(1:pn)) 2554 2555 IF(GotCoeff) THEN 2556 EnergyCoeff(1:n) = ListGetReal( Material, CoeffName, n, NodeIndexes ) 2557 ELSE 2558 EnergyCoeff(1:n) = 1.0d0 2559 END IF 2560 fluxescomputed(bc) = fluxescomputed(bc) + 1 2561 2562 CASE ('area','boundary int','boundary int mean') 2563 IF(GotCoeff) THEN 2564 i = ListGetInteger( Model % Bodies(Element % BodyId) % Values, 'Material', & 2565 minv=1, maxv=Model % NumberOFMaterials ) 2566 Material => Model % Materials(i) % Values 2567 EnergyCoeff(1:n) = ListGetReal( Material, CoeffName, n, NodeIndexes ) 2568 ELSE 2569 EnergyCoeff(1:n) = 1.0d0 2570 END IF 2571 fluxescomputed(bc) = fluxescomputed(bc) + 1 2572 2573 2574 END SELECT 2575 2576!------------------------------------------------------------------------------ 2577! Numerical integration 2578!------------------------------------------------------------------------------ 2579 IntegStuff = GaussPoints( Element ) 2580 2581 DO i=1,IntegStuff % n 2582 U = IntegStuff % u(i) 2583 V = IntegStuff % v(i) 2584 W = IntegStuff % w(i) 2585!------------------------------------------------------------------------------ 2586! Basis function values & derivatives at the integration point 2587!------------------------------------------------------------------------------ 2588 stat = ElementInfo( Element,ElementNodes,U,V,W,SqrtElementMetric, & 2589 Basis,dBasisdx ) 2590!------------------------------------------------------------------------------ 2591! Coordinatesystem dependent info 2592!------------------------------------------------------------------------------ 2593 x = SUM( ElementNodes % x(1:n)*Basis(1:n) ) 2594 y = SUM( ElementNodes % y(1:n)*Basis(1:n) ) 2595 z = SUM( ElementNodes % z(1:n)*Basis(1:n) ) 2596 2597 s = 1.0d0 2598 2599 IF(.FALSE.) THEN 2600 IF ( CurrentCoordinateSystem() /= Cartesian ) THEN 2601 s = 2.0d0 * PI 2602 END IF 2603 CALL CoordinateSystemInfo( Metric,SqrtMetric,Symb,dSymb,x,y,z ) 2604 s = s * SqrtMetric * SqrtElementMetric * IntegStuff % s(i) 2605 ELSE 2606 IF ( CurrentCoordinateSystem() /= Cartesian ) THEN 2607 s = 2.0d0 * PI * x 2608 END IF 2609 s = s * SqrtElementMetric * IntegStuff % s(i) 2610 END IF 2611 2612 2613 SELECT CASE(OperName) 2614 2615 CASE ('diffusive flux') 2616 2617 Normal = NormalVector( Element,ElementNodes,U,V,.TRUE. ) 2618 CALL GetParentUVW( Element, n, Parent, pn, Up, Vp, Wp, Basis) 2619 2620 stat = ElementInfo( Parent,ParentNodes,Up,Vp,Wp,SqrtElementMetric, & 2621 ParentBasis,ParentdBasisdx ) 2622 2623 Flow = 0.0d0 2624 DO j = 1, DIM 2625 DO k = 1, DIM 2626 IF(Permutated) THEN 2627 Flow(j) = Flow(j) + SUM( EnergyTensor(j,k,1:pn) * ParentBasis(1:pn) ) * & 2628 SUM( ParentdBasisdx(1:pn,k) * Var % Values(Var % Perm(ParentIndexes(1:pn))) ) 2629 ELSE 2630 Flow(j) = Flow(j) + SUM( EnergyTensor(j,k,1:pn) * ParentBasis(1:pn) ) * & 2631 SUM( ParentdBasisdx(1:pn,k) * Var % Values(ParentIndexes(1:pn)) ) 2632 END IF 2633 END DO 2634 END DO 2635 2636 flux = SUM(Normal(1:DIM) * Flow(1:DIM)) 2637 2638 IF( SaveFluxRange ) THEN 2639 Minimum = MIN(flux,Minimum) 2640 Maximum = MAX(flux,Maximum) 2641 END IF 2642 2643 IF( PosOper ) flux = MAX( flux, 0.0_dp ) 2644 IF( NegOper ) flux = MIN( flux, 0.0_dp ) 2645 2646 fluxes(bc) = fluxes(bc) + s * flux 2647 2648 CASE ('convective flux') 2649 2650 Normal = NormalVector( Element,ElementNodes,u,v,.TRUE. ) 2651 coeff = SUM( EnergyCoeff(1:n) * Basis(1:n)) 2652 2653 IF(NoDofs == 1) THEN 2654 DO j=1,DIM 2655 Flow(j) = coeff * SUM(LocalVectorSolution(1,1:n) * LocalVeloSolution(j,1:n)*Basis(1:n)) 2656 END DO 2657 ELSE 2658 DO j=1,DIM 2659 Flow(j) = coeff * SUM(LocalVectorSolution(j,1:n)*Basis(1:n)) 2660 END DO 2661 END IF 2662 2663 flux = SUM( Normal(1:dim) * Flow(1:dim) ) 2664 2665 IF( SaveFluxRange ) THEN 2666 Minimum = MIN(flux,Minimum) 2667 Maximum = MAX(flux,Maximum) 2668 END IF 2669 2670 ! Enable positive and negative flux integrals 2671 IF( PosOper ) flux = MAX( flux, 0.0_dp ) 2672 IF( NegOper ) flux = MIN( flux, 0.0_dp ) 2673 2674 fluxes(bc) = fluxes(bc) + s * flux 2675 2676 CASE ('boundary int','boundary int mean') 2677 2678 coeff = SUM( EnergyCoeff(1:n) * Basis(1:n)) 2679 2680 IF(NoDofs == 1) THEN 2681 func = SUM( LocalVectorSolution(1,1:n) * Basis(1:n) ) 2682 flux = coeff * func 2683 IF( PosOper ) flux = MAX( flux, 0.0_dp ) 2684 IF( NegOper ) flux = MIN( flux, 0.0_dp ) 2685 fluxes(bc) = fluxes(bc) + s * flux 2686 ELSE 2687 flux = 0.0_dp 2688 DO j=1,NoDofs 2689 flux = flux + SUM( Var % Values(NoDofs*(PermIndexes(1:n)-1)+j) * Basis(1:n) )**2 2690 END DO 2691 flux = coeff * SQRT( flux ) 2692 IF( PosOper ) flux = MAX( flux, 0.0_dp ) 2693 IF( NegOper ) flux = MIN( flux, 0.0_dp ) 2694 fluxes(bc) = fluxes(bc) + s * flux 2695 END IF 2696 2697 CASE ('area') 2698 coeff = SUM( EnergyCoeff(1:n) * Basis(1:n)) 2699 fluxes(bc) = fluxes(bc) + s * coeff 2700 2701 END SELECT 2702 2703 areas(bc) = areas(bc) + s * coeff 2704 2705 END DO 2706 2707 END DO 2708 2709 END DO 2710 2711 DEALLOCATE( ParentNodes % x, ParentNodes % y, ParentNodes % z, PermIndexes ) 2712 2713 END SUBROUTINE BoundaryIntegrals 2714 2715 2716!------------------------------------------------------------------------------ 2717! Take nodal sum over given boundaries 2718!------------------------------------------------------------------------------ 2719 2720 SUBROUTINE BoundaryStatistics(Var, OperName, GotCoeff, & 2721 CoeffName, fluxes, fluxescomputed) 2722 2723 TYPE(Variable_t), POINTER :: Var 2724 CHARACTER(LEN=MAX_NAME_LEN) :: OperName, CoeffName 2725 LOGICAL :: GotCoeff, FindMinMax 2726 REAL(KIND=dp) :: fluxes(:), val 2727 INTEGER :: fluxescomputed(:) 2728 LOGICAL, ALLOCATABLE :: nodescomputed(:) 2729 TYPE(Element_t), POINTER :: Element, Parent 2730 TYPE(ValueList_t), POINTER :: Material 2731 LOGICAL :: Stat, Permutated 2732 INTEGER :: i,j,k,p,q,t,DIM,bc,n,hits,istat 2733 2734 IF( ASSOCIATED( Var % Perm ) ) THEN 2735 n = SIZE( Var % Perm ) 2736 ELSE 2737 n = Mesh % NumberOfNodes 2738 END IF 2739 ALLOCATE(NodesComputed(n),STAT=istat) 2740 IF( istat /= 0 ) CALL Fatal('BoundaryStatistics','Memory allocation error') 2741 2742 NodesComputed = .FALSE. 2743 hits = 0 2744 2745 NoDofs = Var % Dofs 2746 Permutated = ASSOCIATED(Var % Perm) 2747 FindMinMax = .FALSE. 2748 2749 SELECT CASE(OperName) 2750 2751 CASE('boundary sum','boundary dofs','boundary mean') 2752 fluxes = 0.0_dp 2753 2754 CASE('boundary min','boundary min abs') 2755 fluxes = HUGE( val ) 2756 FindMinMax = .TRUE. 2757 2758 CASE('boundary max','boundary max abs') 2759 fluxes = -HUGE( val ) 2760 FindMinMax = .TRUE. 2761 2762 CASE DEFAULT 2763 CALL Warn('SaveScalars','Unknown statistical operator') 2764 2765 END SELECT 2766 2767 2768 DO t = Mesh % NumberOfBulkElements+1, Mesh % NumberOfBulkElements & 2769 + Mesh % NumberOfBoundaryElements 2770 2771 Element => Mesh % Elements(t) 2772 Model % CurrentElement => Mesh % Elements(t) 2773 2774 n = Element % TYPE % NumberOfNodes 2775 NodeIndexes => Element % NodeIndexes 2776 2777 DO bc=1, Model % NumberOfBCs 2778 2779 IF ( Model % BCs(bc) % Tag /= Element % BoundaryInfo % Constraint ) CYCLE 2780 IF(.NOT. ListGetLogical(Model % BCs(bc) % Values, MaskName, gotIt ) ) CYCLE 2781 2782 hits = hits + 1 2783 2784 DO i=1,n 2785 j = NodeIndexes(i) 2786 2787 IF( IsParallel ) THEN 2788 IF( Mesh % ParallelInfo % NeighbourList(j) % Neighbours(1) /= ParEnv % MyPE ) CYCLE 2789 END IF 2790 2791 IF( Permutated ) j = Var % Perm(j) 2792 2793 IF( j == 0) CYCLE 2794 IF( nodescomputed(j) ) CYCLE 2795 2796 IF( NoDofs == 1 ) THEN 2797 val = Var % Values(j) 2798 ELSE 2799 val = 0.0_dp 2800 DO k=1,NoDofs 2801 val = val + Var % Values(NoDofs*(j-1)+k) 2802 END DO 2803 val = SQRT( val ) 2804 END IF 2805 2806 2807 IF(FindMinMax) THEN 2808 SELECT CASE(OperName) 2809 2810 CASE('boundary min') 2811 fluxes(bc) = MIN( val, fluxes(bc) ) 2812 2813 CASE('boundary max') 2814 fluxes(bc) = MAX( val, fluxes(bc) ) 2815 2816 CASE('boundary min abs') 2817 IF(ABS(fluxes(bc)) < ABS( val )) fluxes(bc) = val 2818 2819 CASE('boundary max abs') 2820 IF(ABS(fluxes(bc)) > ABS( val )) fluxes(bc) = val 2821 2822 END SELECT 2823 ELSE 2824 fluxes(bc) = fluxes(bc) + val 2825 END IF 2826 2827 nodescomputed(j) = .TRUE. 2828 fluxescomputed(bc) = fluxescomputed(bc) + 1 2829 END DO 2830 END DO 2831 END DO 2832 2833 2834 SELECT CASE(OperName) 2835 2836 CASE('boundary dofs') 2837 fluxes = 1.0 * fluxescomputed 2838 2839 END SELECT 2840 2841 END SUBROUTINE BoundaryStatistics 2842 2843 2844!------------------------------------------------------------------------------ 2845 2846 2847 SUBROUTINE PolylineIntegrals(Var, OperName, GotCoeff, & 2848 CoeffName, fluxes, areas, fluxescomputed) 2849 2850 TYPE(Variable_t), POINTER :: Var 2851 CHARACTER(LEN=MAX_NAME_LEN) :: OperName, CoeffName 2852 LOGICAL :: GotCoeff 2853 REAL(KIND=dp) :: fluxes(:),areas(:) 2854 INTEGER :: fluxescomputed(:) 2855 2856 INTEGER :: t 2857 TYPE(Element_t), TARGET :: SideElement 2858 TYPE(Element_t), POINTER :: Element, Parent 2859 TYPE(ValueList_t), POINTER :: Material 2860 REAL(KIND=dp) :: SqrtMetric,Metric(3,3),Symb(3,3,3),dSymb(3,3,3,3),LocalCoordinates(3),Point(3) 2861 REAL(KIND=dp) :: Basis(Model % MaxElementNodes),dBasisdx(Model % MaxElementNodes,3),& 2862 ParentBasis(Model % MaxElementNodes),& 2863 ParentdBasisdx(Model % MaxElementNodes,3),EnergyTensor(3,3,Model % MaxElementNodes),& 2864 EnergyCoeff(Model % MaxElementNodes) 2865 REAL(KIND=dp) :: SqrtElementMetric,U,V,W,up,vp,wp,S,A,L,C(3,3),x,y,z,dx,dy,dz,ds,dsmax 2866 REAL(KIND=dp) :: func, coeff, Normal(3), Flow(3), x0, y0, z0, pos(2), flux 2867 REAL(KIND=DP), POINTER :: Pwrk(:,:,:) 2868 INTEGER, POINTER :: ParentIndexes(:), PermIndexes(:), SideIndexes(:), OnLine(:,:) 2869 2870 LOGICAL :: Stat, Permutated, Inside 2871 INTEGER :: i,j,k,p,q,DIM,bc,NoDofs,pn,Line, NoSides, Side, NodeNumber, LineNode(2), istat 2872 TYPE(GaussIntegrationPoints_t) :: IntegStuff 2873 TYPE(Nodes_t) :: ParentNodes, LineNodes, SideNodes 2874 2875 2876 n = Model % MaxElementNodes 2877 DIM = CoordinateSystemDimension() 2878 2879 ALLOCATE(ParentNodes % x(n), ParentNodes % y(n), ParentNodes % z(n), PermIndexes(n), & 2880 SideNodes % x(n), SideNodes % y(n), SideNodes % z(n), SideIndexes(n), & 2881 LineNodes % x(n), LineNodes % y(n), LineNodes % z(n), & 2882 OnLine(Mesh % NumberOfNodes,2), STAT=istat ) 2883 IF( istat /= 0 ) CALL Fatal('PolylineIntegrals','Memory allocation error') 2884 2885 2886 CALL Info('PolylineIntegrals','OPERATOR: '//TRIM(OperName)) 2887 2888 areas = 0.0_dp 2889 coeff = 1.0_dp 2890 2891 NoDofs = Var % Dofs 2892 IF( OperName == 'area' ) THEN 2893 Permutated = .FALSE. 2894 ELSE 2895 Permutated = ASSOCIATED(Var % Perm) 2896 END IF 2897 2898 DO i=1,3 2899 EnergyTensor(i,i,:) = 1.0d0 2900 END DO 2901 2902 SELECT CASE(OperName) 2903 2904 CASE('diffusive flux') 2905 IF(NoDofs /= 1) THEN 2906 CALL Warn('SaveScalars','diffusive flux & NoDofs /= 1?') 2907 RETURN 2908 END IF 2909 2910 CASE ('convective flux') 2911 IF(NoDofs /= 1 .AND. NoDofs < DIM) THEN 2912 CALL Warn('SaveScalars','convective flux & NoDofs < DIM?') 2913 RETURN 2914 END IF 2915 2916 CASE ('area','boundary int','boundary int mean') 2917 2918 CASE DEFAULT 2919 CALL Warn('SaveScalars','Unknown physical OPERATOR') 2920 2921 END SELECT 2922 2923 DIM = CoordinateSystemDimension() 2924 IF(DIM < 3) THEN 2925 LineNodes % z = 0.0d0 2926 SideNodes % z(1:2) = 0.0d0 2927 END IF 2928 2929 fluxes = 0.0d0 2930 SideElement % TYPE => GetElementType( 202, .FALSE.) 2931 SideElement % Bdofs = 0 2932 Element => SideElement 2933 2934! /* Go through the line segments */ 2935 DO Line = 1,NoLines 2936 2937 LineNodes % x(1:2) = LineCoordinates(2*Line-1:2*Line,1) 2938 LineNodes % y(1:2) = LineCoordinates(2*Line-1:2*Line,2) 2939 IF(DIM == 3) LineNodes % z(1:2) = LineCoordinates(2*Line-1:2*Line,3) 2940 OnLine = 0 2941 2942 DO t = 1, Mesh % NumberOfBulkElements 2943 2944 Parent => Mesh % Elements(t) 2945 Model % CurrentElement => Mesh % Elements(t) 2946 2947 NoSides = Parent % TYPE % ElementCode / 100 2948 IF(NoSides < 3 .OR. NoSides > 4) CYCLE 2949 2950 pn = Parent % TYPE % NumberOfNodes 2951 ParentIndexes => Parent % NodeIndexes 2952 2953 ParentNodes % x(1:pn) = Mesh % Nodes % x(ParentIndexes(1:pn)) 2954 ParentNodes % y(1:pn) = Mesh % Nodes % y(ParentIndexes(1:pn)) 2955 ParentNodes % z(1:pn) = Mesh % Nodes % z(ParentIndexes(1:pn)) 2956 2957 IF(Permutated) THEN 2958 PermIndexes(1:pn) = Var % Perm(ParentIndexes(1:pn)) 2959 IF (ANY( PermIndexes(1:pn) == 0)) CYCLE 2960 ELSE 2961 PermIndexes(1:pn) = ParentIndexes(1:pn) 2962 END IF 2963 2964 NodeNumber = 0 2965 2966 DO Side = 1, NoSides 2967 2968 SideIndexes(1) = ParentIndexes(Side) 2969 SideIndexes(2) = ParentIndexes(MOD(Side,NoSides)+1) 2970 2971 SideNodes % x(1:2) = Mesh % Nodes % x(SideIndexes(1:2)) 2972 SideNodes % y(1:2) = Mesh % Nodes % y(SideIndexes(1:2)) 2973 IF(DIM == 3) SideNodes % z(1:2) = Mesh % Nodes % z(SideIndexes(1:2)) 2974 2975 CALL LineIntersectionCoords(SideNodes,LineNodes,Inside,x0,y0,z0,u) 2976 2977 IF(.NOT. Inside) CYCLE 2978 2979 NodeNumber = NodeNumber + 1 2980 ElementNodes % x(NodeNumber) = x0 2981 ElementNodes % y(NodeNumber) = y0 2982 ElementNodes % z(NodeNumber) = z0 2983 pos(NodeNumber) = u 2984 2985 END DO 2986 2987 IF(NodeNumber == 0) CYCLE 2988 2989 IF(NodeNumber > 2) THEN 2990 CALL Warn('PolylineIntergrals','There should not be more than 2 intersections!') 2991 CYCLE 2992 END IF 2993 2994 !--------------------------------------------------------------------------- 2995 ! If there is only one intersection the other end of the node must lie 2996 ! inside the element. Assuming that the line is long compared to the 2997 ! element the correct end of the line segment may be easily deduced. 2998 !--------------------------------------------------------------------------- 2999 IF(NodeNumber == 1) THEN 3000 IF(pos(1) < 0.5d0) THEN 3001 i = 1 3002 pos(2) = 0.0 3003 ELSE 3004 i = 2 3005 pos(2) = 1.0 3006 END IF 3007 x0 = LineNodes % x(i) 3008 y0 = LineNodes % y(i) 3009 z0 = LineNodes % z(i) 3010 3011 ElementNodes % x(2) = LineNodes % x(i) 3012 ElementNodes % y(2) = LineNodes % y(i) 3013 ElementNodes % z(2) = LineNodes % z(i) 3014 END IF 3015 3016 IF(ABS(pos(1)-pos(2)) < 1.0d-8) CYCLE 3017 3018 !----------------------------------------------------------------------------- 3019 ! Change the order of nodes so that the normal always points to the same direction 3020 !----------------------------------------------------------------------------- 3021 IF(pos(1) < pos(2)) THEN 3022 ElementNodes % x(2) = ElementNodes % x(1) 3023 ElementNodes % y(2) = ElementNodes % y(1) 3024 ElementNodes % z(2) = ElementNodes % z(1) 3025 ElementNodes % x(1) = x0 3026 ElementNodes % y(1) = y0 3027 ElementNodes % z(1) = z0 3028 END IF 3029 3030 !-------------------------------------------------------------------------------- 3031 ! The following avoids the cases where the line goes exactly at the element 3032 ! interface and therefore the flux would be computed twice 3033 !-------------------------------------------------------------------------------- 3034 dx = ElementNodes % x(1) - ElementNodes % x(2) 3035 dy = ElementNodes % y(1) - ElementNodes % y(2) 3036 dsmax = SQRT(dx*dx+dy*dy) 3037 LineNode = 0 3038 3039 3040 DO i=1,Parent % TYPE % ElementCode / 100 3041 DO j=1,2 3042 dx = ParentNodes % x(i) - ElementNodes % x(j) 3043 dy = ParentNodes % y(i) - ElementNodes % y(j) 3044 ds = SQRT(dx*dx+dy*dy) 3045 IF(ds < 1.0d-4 * dsmax) LineNode(j) = ParentIndexes(i) 3046 END DO 3047 END DO 3048 3049 IF(ALL(LineNode(1:2) > 0)) THEN 3050 IF(ANY(OnLine(LineNode(1),:) == LineNode(2))) CYCLE 3051 IF(ANY(OnLine(LineNode(2),:) == LineNode(1))) CYCLE 3052 3053 IF(OnLine(LineNode(1),1) == 0) THEN 3054 OnLine(LineNode(1),1) = LineNode(2) 3055 ELSE IF(OnLine(LineNode(1),2) == 0) THEN 3056 OnLine(LineNode(1),2) = LineNode(2) 3057 ELSE 3058 CALL Warn('PolylineIntegrate','This should never happen') 3059 END IF 3060 3061 IF(OnLine(LineNode(2),1) == 0) THEN 3062 OnLine(LineNode(2),1) = LineNode(1) 3063 ELSE IF(OnLine(LineNode(2),2) == 0) THEN 3064 OnLine(LineNode(2),2) = LineNode(1) 3065 ELSE 3066 CALL Warn('PolylineIntegrate','This should never happen') 3067 END IF 3068 END IF 3069 3070 3071 3072 i = ListGetInteger( Model % Bodies(Parent % BodyId) % Values, 'Material', & 3073 minv=1, maxv=Model % NumberOFMaterials ) 3074 Material => Model % Materials(i) % Values 3075 fluxescomputed(Line) = fluxescomputed(Line) + 1 3076 3077 3078 SELECT CASE(OperName) 3079 3080 CASE('diffusive flux') 3081 3082 IF(GotCoeff) THEN 3083 CALL ListGetRealArray( Material, CoeffName, Pwrk, & 3084 pn, ParentIndexes, gotIt ) 3085 3086 EnergyTensor = 0._dp 3087 IF(GotIt) THEN 3088 IF ( SIZE(Pwrk,1) == 1 ) THEN 3089 DO i=1,3 3090 EnergyTensor( i,i,1:pn ) = Pwrk( 1,1,1:pn ) 3091 END DO 3092 ELSE IF ( SIZE(Pwrk,2) == 1 ) THEN 3093 DO i=1,MIN(3,SIZE(Pwrk,1)) 3094 EnergyTensor(i,i,1:pn) = Pwrk(i,1,1:pn) 3095 END DO 3096 ELSE 3097 DO i=1,MIN(3,SIZE(Pwrk,1)) 3098 DO j=1,MIN(3,SIZE(Pwrk,2)) 3099 EnergyTensor( i,j,1:pn ) = Pwrk(i,j,1:pn) 3100 END DO 3101 END DO 3102 END IF 3103 ELSE 3104 DO i=1,3 3105 EnergyTensor(i,i,1:pn) = 1.0d0 3106 END DO 3107 END IF 3108 END IF 3109 3110 CASE ('convective flux','area') 3111 EnergyCoeff(1:n) = ListGetReal( Material, CoeffName, pn, ParentIndexes, gotIt ) 3112 IF(.NOT. GotIt) EnergyCoeff(1:pn) = 1.0d0 3113 3114 END SELECT 3115 3116!------------------------------------------------------------------------------ 3117! Numerical integration 3118!------------------------------------------------------------------------------ 3119 3120 IntegStuff = GaussPoints( Element, 2 ) 3121 3122 DO i=1,IntegStuff % n 3123 U = IntegStuff % u(i) 3124 V = IntegStuff % v(i) 3125 W = IntegStuff % w(i) 3126 3127!------------------------------------------------------------------------------ 3128! Basis function values & derivatives at the integration point 3129!------------------------------------------------------------------------------ 3130 stat = ElementInfo( Element,ElementNodes,U,V,W,SqrtElementMetric, & 3131 Basis,dBasisdx ) 3132 3133 x = SUM( ElementNodes % x(1:n) * Basis(1:n) ) 3134 y = SUM( ElementNodes % y(1:n) * Basis(1:n) ) 3135 z = SUM( ElementNodes % z(1:n) * Basis(1:n) ) 3136 3137!------------------------------------------------------------------------------ 3138! Coordinatesystem dependent info 3139!------------------------------------------------------------------------------ 3140 3141 s = 1.0d0 3142 3143 IF(.FALSE.) THEN 3144 IF(CurrentCoordinateSystem() /= Cartesian ) s = 2.0 * PI 3145 CALL CoordinateSystemInfo( Metric,SqrtMetric,Symb,dSymb,x,y,z ) 3146 s = s * SqrtMetric * SqrtElementMetric * IntegStuff % s(i) 3147 ELSE 3148 IF(CurrentCoordinateSystem() /= Cartesian ) s = 2.0 * PI * x 3149 s = s * SqrtElementMetric * IntegStuff % s(i) 3150 END IF 3151 3152 Normal = NormalVector( Element,ElementNodes,U,V,.FALSE. ) 3153 3154!------------------------------------------------------------------------------ 3155! Because the intersection nodes do not really exist the field variables 3156! must be evaluated using the nodes of the parent element. 3157!------------------------------------------------------------------------------ 3158 3159 Point(1) = x 3160 Point(2) = y 3161 Point(3) = z 3162 3163 IF( .NOT. PointInElement( Parent, ParentNodes, Point, LocalCoordinates ) ) THEN 3164 CALL Warn('PolylineIntegrals','The node should be in the element by construction!') 3165 CYCLE 3166 END IF 3167 3168 Up = LocalCoordinates(1) 3169 Vp = LocalCoordinates(2) 3170 Wp = LocalCoordinates(3) 3171 3172 3173 stat = ElementInfo( Parent,ParentNodes,Up,Vp,Wp,SqrtElementMetric, & 3174 ParentBasis,ParentdBasisdx ) 3175 3176 3177 SELECT CASE(OperName) 3178 3179 CASE ('diffusive flux') 3180 Flow = 0.0d0 3181 DO j = 1, DIM 3182 DO k = 1, DIM 3183 Flow(j) = Flow(j) + SUM( EnergyTensor(j,k,1:pn) * ParentBasis(1:pn) ) * & 3184 SUM( ParentdBasisdx(1:pn,k) * Var % Values(PermIndexes(1:pn)) ) 3185 END DO 3186 END DO 3187 fluxes(Line) = fluxes(Line) + s * SUM(Normal(1:DIM) * Flow(1:DIM)) 3188 3189 3190 CASE ('convective flux') 3191 coeff = SUM( EnergyCoeff(1:pn) * ParentBasis(1:pn)) 3192 IF(NoDofs == 1) THEN 3193 func = SUM( Var % Values(PermIndexes(1:pn)) * ParentBasis(1:pn) ) 3194 fluxes(Line) = fluxes(Line) + s * coeff * func 3195 ELSE 3196 DO j=1,DIM 3197 Flow(j) = coeff * & 3198 SUM( Var % Values(NoDofs*(PermIndexes(1:pn)-1)+j) * ParentBasis(1:pn) ) 3199 END DO 3200 fluxes(Line) = fluxes(Line) + s * coeff * SUM(Normal * Flow) 3201 END IF 3202 3203 CASE ('boundary int','boundary int mean') 3204 coeff = SUM( EnergyCoeff(1:pn) * ParentBasis(1:pn)) 3205 IF(NoDofs == 1) THEN 3206 func = SUM( Var % Values(PermIndexes(1:pn)) * ParentBasis(1:pn) ) 3207 fluxes(Line) = fluxes(Line) + s * coeff * func 3208 ELSE 3209 flux = 0.0_dp 3210 DO j=1,NoDofs 3211 flux = flux + & 3212 SUM( Var % Values(NoDofs*(PermIndexes(1:pn)-1)+j) * ParentBasis(1:pn) )**2 3213 END DO 3214 flux = coeff * SQRT(flux) 3215 fluxes(Line) = fluxes(Line) + s * flux 3216 END IF 3217 3218 CASE ('area') 3219 coeff = SUM( EnergyCoeff(1:pn) * Basis(1:pn)) 3220 fluxes(Line) = fluxes(Line) + s * coeff 3221 3222 END SELECT 3223 3224 areas(Line) = areas(Line) + s * coeff 3225 3226 END DO 3227 3228 END DO 3229 3230 END DO 3231 3232 DEALLOCATE( ParentNodes % x, ParentNodes % y, ParentNodes % z, PermIndexes, & 3233 SideNodes % x, SideNodes % y, SideNodes % z, SideIndexes, & 3234 LineNodes % x, LineNodes % y, LineNodes % z, OnLine) 3235 3236 END SUBROUTINE PolylineIntegrals 3237!------------------------------------------------------------------------------ 3238 3239 3240!------------------------------------------------------------------------------ 3241!> This subroutine tests whether the line segment goes through the current 3242!> face of the element. 3243!------------------------------------------------------------------------------ 3244 SUBROUTINE LineIntersectionCoords(Plane,Line,Inside,x0,y0,z0,frac) 3245 3246 TYPE(Nodes_t) :: Plane, Line 3247 LOGICAL :: Inside 3248 REAL (KIND=dp) :: x0, y0, z0, frac 3249 3250 REAL (KIND=dp) :: A(3,3),B(3),C(3),eps=1.0d-6,detA,absA 3251 3252 Inside = .FALSE. 3253 3254 ! In 2D the intersection is between two lines 3255 A(1,1) = Line % x(2) - Line % x(1) 3256 A(2,1) = Line % y(2) - Line % y(1) 3257 A(1,2) = Plane % x(1) - Plane % x(2) 3258 A(2,2) = Plane % y(1) - Plane % y(2) 3259 3260 detA = A(1,1)*A(2,2)-A(1,2)*A(2,1) 3261 absA = SUM(ABS(A(1,1:2))) * SUM(ABS(A(2,1:2))) 3262 3263 IF(ABS(detA) <= eps * absA + 1.0d-20) RETURN 3264 3265 B(1) = Plane % x(1) - Line % x(1) 3266 B(2) = Plane % y(1) - Line % y(1) 3267 3268 CALL InvertMatrix( A,2 ) 3269 C(1:2) = MATMUL(A(1:2,1:2),B(1:2)) 3270 3271 IF(ANY(C(1:2) < 0.0) .OR. ANY(C(1:2) > 1.0d0)) RETURN 3272 3273 Inside = .TRUE. 3274 frac = C(1) 3275 X0 = Line % x(1) + C(1) * (Line % x(2) - Line % x(1)) 3276 Y0 = Line % y(1) + C(1) * (Line % y(2) - Line % y(1)) 3277 Z0 = Line % z(1) + C(1) * (Line % z(2) - Line % z(1)) 3278 3279 END SUBROUTINE LineIntersectionCoords 3280 3281!------------------------------------------------------------------------------ 3282END SUBROUTINE SaveScalars 3283!------------------------------------------------------------------------------ 3284 3285!> \} 3286