1!/*****************************************************************************/ 2! * 3! * Elmer, A Finite Element Software for Multiphysical Problems 4! * 5! * Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland 6! * 7! * This program is free software; you can redistribute it and/or 8! * modify it under the terms of the GNU General Public License 9! * as published by the Free Software Foundation; either version 2 10! * of the License, or (at your option) any later version. 11! * 12! * This program is distributed in the hope that it will be useful, 13! * but WITHOUT ANY WARRANTY; without even the implied warranty of 14! * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15! * GNU General Public License for more details. 16! * 17! * You should have received a copy of the GNU General Public License 18! * along with this program (in file fem/GPL-2); if not, write to the 19! * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 20! * Boston, MA 02110-1301, USA. 21! * 22! *****************************************************************************/ 23! 24!/****************************************************************************** 25! * 26! * Authors: Peter Raback & Juha Ruokolainen 27! * Email: Peter.Raback@csc.fi & Juha.Ruokolainen@csc.fi 28! * Web: http://www.csc.fi/elmer 29! * Address: CSC - IT Center for Science Ltd. 30! * Keilaranta 14 31! * 02101 Espoo, Finland 32! * 33! * Original Date: 16.6.2011 34! * 35! *****************************************************************************/ 36 37 38!------------------------------------------------------------------------------- 39!> Subroutine for advecting fields in time using particles to follow them 40!> backwards in time, and taking the field value from the given point. This should overcome 41!> all problems with diffusion. 42!> This is a dynamically loaded solver with a standard interface. 43!> \ingroup Solvers 44!------------------------------------------------------------------------------- 45SUBROUTINE ParticleAdvector( Model,Solver,dt,TransientSimulation ) 46!------------------------------------------------------------------------------ 47 USE DefUtils 48 USE Interpolation 49 USE MeshUtils 50 USE ElementUtils 51 USE ParticleUtils 52 53 IMPLICIT NONE 54!------------------------------------------------------------------------------ 55 TYPE(Solver_t), TARGET :: Solver 56 TYPE(Model_t) :: Model 57 REAL(KIND=dp) :: dt 58 LOGICAL :: TransientSimulation 59!------------------------------------------------------------------------------ 60! Local variables 61!------------------------------------------------------------------------------ 62 63 TYPE(Mesh_t), POINTER :: Mesh 64 TYPE(ValueList_t), POINTER :: Params 65 TYPE(Solver_t), POINTER :: PSolver 66 TYPE(Variable_t), POINTER :: Var, PTimeVar 67 LOGICAL :: GotIt, Debug, Hit, InitLocation, InitTimestep, Found, ParticleInfo, InitAllVelo 68 INTEGER :: i,j,k,n,dim,No,nodims,& 69 ElementIndex, VisitedTimes = 0, nstep, & 70 Status,TimeOrder, PartitionChanges, TimeStepsTaken=0,& 71 ParticleStepsTaken=0, TotParticleStepsTaken, TotNoParticles, & 72 istep,iorder,NoMoving 73 REAL(KIND=dp) :: maxdt, dertime = 0.0, tottime = 0.0 74 CHARACTER(LEN=MAX_NAME_LEN) :: VariableName, IntegMethod 75 TYPE(Particle_t), POINTER :: Particles 76 77 SAVE Nstep, VisitedTimes, TimeOrder, & 78 tottime, TimeStepsTaken, ParticleStepsTaken, ParticleInfo 79 80!------------------------------------------------------------------------------ 81 82 CALL Info('ParticleAdvector','-----------------------------------------', Level=4 ) 83 CALL Info('ParticleAdvector','Advecting fields using particle tracking',Level=4) 84 85 Particles => GlobalParticles 86 VisitedTimes = VisitedTimes + 1 87 88 Params => GetSolverParams() 89 Mesh => Solver % Mesh 90 PSolver => Solver 91 DIM = CoordinateSystemDimension() 92 93 maxdt = 0.0_dp 94 istep = 1 95 iorder = 1 96 97 InitAllVelo = .TRUE. 98 99 ! Do some initalialization: allocate space, check fields 100 !------------------------------------------------------------------------ 101 IF( VisitedTimes == 1 ) THEN 102 TimeOrder = GetInteger( Params,'Time Order',GotIt) 103 CALL SetParticlePreliminaries( Particles, dim, TimeOrder ) 104 Nstep = GetInteger( Params,'Max Timestep Intervals',Found) 105 IF(.NOT. Found) Nstep = 1 106 ParticleInfo = GetLogical( Params,'Particle Info',Found) 107 END IF 108 109 ! Initialize particles always since we just advance -dt each time 110 !------------------------------------------------------------------------- 111 IF( VisitedTimes == 1 .OR. GetLogical( Params,'Reinitialize Particles',GotIt ) ) THEN 112 CALL InitializeParticles( Particles, SaveOrigin = .TRUE. ) 113 CALL ReleaseWaitingParticles(Particles) 114 Particles % Status = PARTICLE_LOCATED 115 ELSE 116 ! in case the velocity field is changed update also the particle velocities 117 CALL SetParticleVelocities(InitAllVelo) 118 InitAllVelo = .FALSE. 119 END IF 120 121 IF( VisitedTimes == 1 ) THEN 122 IF( GetLogical( Params,'Particle Time',Found) ) THEN 123 CALL ParticleVariableCreate( Particles,'particle time') 124 END IF 125 CALL ParticleVariableCreate( Particles,'particle distance') 126 ELSE 127 PtimeVar => ParticleVariableGet( Particles, 'particle time' ) 128 IF( ASSOCIATED( PTimeVar ) ) THEN 129 PTimeVar % Values = 0.0_dp 130 END IF 131 END IF 132 133 ! Freeze particles that are known not to move (e.g. no-slip wall) 134 !---------------------------------------------------------------- 135 CALL SetFixedParticles( ) 136 137 ! We are advecting backwards in time! 138 !---------------------------------------------------- 139 InitTimestep = .TRUE. 140 IF( Particles % RK2 ) THEN 141 iorder = 2 142 ELSE 143 iorder = 1 144 END IF 145 146! CALL ParticleStatusCount( Particles ) 147 148 149 DO i=1,nstep 150 151 ! Get the timestep size, initialize at 1st round 152 !-------------------------------------------------------------- 153 maxdt = GetParticleTimeStep( Particles, InitTimestep ) 154 155 ! If size of timestep goes to zero then no more steps are needed 156 !--------------------------------------------------------------- 157 IF( ABS( maxdt ) < TINY( maxdt ) ) THEN 158 WRITE (Message,'(A,I0)') 'Number of steps used: ',i-1 159 CALL Info('ParticleAdvector',Message,Level=6) 160 EXIT 161 END IF 162 163 dertime = dertime + maxdt 164 tottime = tottime + maxdt 165 166 TimeStepsTaken = TimeStepsTaken + 1 167 ParticleStepsTaken = ParticleStepsTaken + Particles % NumberOfParticles 168 169 ! If there are periodic BCs apply them just before locating the particles 170 !------------------------------------------------------------------------ 171 ! CALL ParticleBoxPeriodic( Particles ) 172 173 174 ! Loop over possible runge-kutta steps 175 !------------------------------------------------------------------------ 176 177 DO istep = 1, iorder 178 179 ! Set velocity field at points. This is not done the first time. 180 ! This has to be here since the velocity field could have changed 181 ! between calls. 182 !------------------------------------------------------------------ 183 !IF( .NOT. InitTimestep ) CALL SetParticleVelocities() 184 185 ! Advance the coordinate, r = r0 + v * dt 186 ! either with 2nd order Runge-Kutta, or with 1st order explicit scheme. 187 !------------------------------------------------------------------ 188 CALL ParticleAdvanceTimestep( Particles, istep ) 189 190 IF( InfoActive( 20 ) ) THEN 191 CALL ParticleStatusCount( Particles ) 192 END IF 193 194 ! Find the elements (and only the elements) in which the particles are in. 195 !------------------------------------------------------------------------ 196 CALL LocateParticles( Particles ) 197 198 CALL SetParticleVelocities(InitAllVelo) 199 InitAllVelo = .FALSE. 200 201 ! Integrate over the particle path (\int f(r) ds or \int f(r) dt ) 202 !------------------------------------------------------------------ 203 CALL ParticlePathIntegral( Particles, istep ) 204 205 InitTimestep = .FALSE. 206 END DO 207 208 NoMoving = Particles % NumberOfMovingParticles 209 NoMoving = NINT( ParallelReduction( 1.0_dp * NoMoving ) ) 210 WRITE (Message,'(A,I0,A,I0,A)') 'Timestep ',i,' with ',NoMoving,' moving particles' 211 CALL Info('ParticleAdvector',Message,Level=6) 212 213 ! Freeze particles that are too old 214 !---------------------------------------------------------------- 215 CALL SetRetiredParticles( ) 216 217 IF( InfoActive( 15 ) ) THEN 218 CALL ParticleInformation(Particles, ParticleStepsTaken, & 219 TimeStepsTaken, tottime ) 220 END IF 221 222 END DO 223 224 ! Set the advected field giving the final locations of the particles backward in time 225 !------------------------------------------------------------------------------------ 226 CALL SetAdvectedField() 227 228 ! In the end show some statistical info 229 !--------------------------------------------------------------- 230 IF( ParticleInfo ) THEN 231 CALL ParticleInformation(Particles, ParticleStepsTaken, & 232 TimeStepsTaken, tottime ) 233 END IF 234 235 CALL Info('ParticleAdvector','All done',Level=4) 236 CALL Info('ParticleAdvector', '-----------------------------------------', Level=4 ) 237 238 239CONTAINS 240 241 !> Eliminate particles that sit an fixed boundaries. 242 !------------------------------------------------------------------------- 243 SUBROUTINE SetFixedParticles( ) 244 245 INTEGER :: i,j,k 246 TYPE( ValueList_t), POINTER :: BC, BodyForce 247 TYPE( Element_t ), POINTER :: Element 248 REAL(KIND=dp), ALLOCATABLE :: PTime(:), PCond(:) 249 REAL(KIND=dp) :: PTimeConst 250 LOGICAL :: Found, GotTime, SomeBodyForce, SomeBC 251 INTEGER :: FixedCount 252 TYPE(ValueList_t), POINTER :: Params 253 254 255 Params => GetSolverParams() 256 257 SomeBC = ListCheckPresentAnyBC( Model,'Particle Fixed Condition') 258 SomeBodyForce = ListCheckPresentAnyBodyForce( Model,'Particle Fixed Condition') 259 260 IF( .NOT. (SomeBC .OR. SomeBodyForce ) ) RETURN 261 262 i = Model % MaxElementNodes 263 ALLOCATE( PTime(i), PCond(i) ) 264 265 PtimeVar => ParticleVariableGet( Particles, 'particle time' ) 266 GotTime = ASSOCIATED( PTimeVar ) 267 268 PTimeConst = ListGetCReal( Params,'Fixed Particle Time',Found) 269 FixedCount = 0 270 271 IF( SomeBC ) THEN 272 DO j=Mesh % NumberOfBulkElements+1,& 273 Mesh % NumberOfBulkElements + Mesh % NumberOfBoundaryElements 274 Element => Mesh % Elements(j) 275 Model % CurrentElement => Element 276 BC => GetBC( Element ) 277 n = GetElementNOFNodes() 278 279 PCond(1:n) = GetReal( BC,'Particle Fixed Condition',Found) 280 IF(.NOT. Found ) CYCLE 281 282 PTime(1:n) = GetReal( BC,'Particle Time',Found) 283 284 DO i=1,n 285 IF( PCond(i) < 0.0_dp ) CYCLE 286 287 k = Element % NodeIndexes(i) 288 IF( Particles % Status( k ) == PARTICLE_FIXEDCOORD ) CYCLE 289 290 Particles % Status( k ) = PARTICLE_FIXEDCOORD 291 FixedCount = FixedCount + 1 292 IF(.NOT. GotTime ) CYCLE 293 294 IF( Found ) THEN 295 PTimeVar % Values( k ) = PTime(i) 296 ELSE 297 PTimeVar % Values( k ) = PTimeConst 298 END IF 299 END DO 300 END DO 301 END IF 302 303 304 IF( SomeBodyForce ) THEN 305 DO j=1, Mesh % NumberOfBulkElements 306 Element => Mesh % Elements(j) 307 Model % CurrentElement => Element 308 BodyForce => GetBodyForce( Element ) 309 n = GetElementNOFNodes() 310 311 PCond(1:n) = GetReal( BodyForce,'Particle Fixed Condition',Found) 312 IF(.NOT. Found ) CYCLE 313 314 PTime(1:n) = GetReal( BodyForce,'Particle Time',Found) 315 316 DO i=1,n 317 IF( PCond(i) < 0.0_dp ) CYCLE 318 319 k = Element % NodeIndexes(i) 320 IF( Particles % Status( k ) == PARTICLE_FIXEDCOORD ) CYCLE 321 322 Particles % Status( k ) = PARTICLE_FIXEDCOORD 323 FixedCount = FixedCount + 1 324 IF(.NOT. GotTime ) CYCLE 325 326 IF( Found ) THEN 327 PTimeVar % Values( k ) = PTime(i) 328 ELSE 329 PTimeVar % Values( k ) = PTimeConst 330 END IF 331 END DO 332 END DO 333 END IF 334 DEALLOCATE( PTime, PCond ) 335 336 CALL Info('SetFixedParticles','Number of new fixed particles: '& 337 //TRIM(I2S(FixedCount)),Level=7) 338 339 END SUBROUTINE SetFixedParticles 340 341 342 !> Eliminate particles that are too old sit an fixed boundaries. 343 !------------------------------------------------------------------------- 344 SUBROUTINE SetRetiredParticles( ) 345 346 INTEGER :: i,j,k 347 REAL(KIND=dp) :: MaxIntegTime 348 LOGICAL :: Found, GotMaxIntegTime 349 INTEGER :: FixedCount 350 TYPE(ValueList_t), POINTER :: Params 351 352 Params => GetSolverParams() 353 354 MaxIntegTime = ListGetCReal( Params,'Max Integration Time',GotMaxIntegTime ) 355 356 IF( .NOT. GotMaxIntegTime ) RETURN 357 358 PtimeVar => ParticleVariableGet( Particles, 'particle time' ) 359 IF(.NOT. ASSOCIATED( PTimeVar ) ) THEN 360 CALL Fatal('SetRetiredParticles','Cannot retire particles without time!') 361 END IF 362 363 FixedCount = 0 364 DO k=1,SIZE( PTimeVar % Values ) 365 IF( Particles % Status( k ) == PARTICLE_FIXEDCOORD ) CYCLE 366 IF( PTimeVar % Values( k ) >= MaxIntegTime ) THEN 367 FixedCount = FixedCount + 1 368 Particles % Status( k ) = PARTICLE_FIXEDCOORD 369 END IF 370 END DO 371 372 IF( FixedCount > 0 ) THEN 373 CALL Info('SetRetiredParticles','Number of new retired particles: '& 374 //TRIM(I2S(FixedCount)),Level=7) 375 END IF 376 377 END SUBROUTINE SetRetiredParticles 378 379 380 !------------------------------------------------------------------------ 381 !> Compute field values at the given points in the FE mesh. 382 !------------------------------------------------------------------------- 383 SUBROUTINE SetParticleVelocities( FirstStep ) 384 LOGICAL :: FirstStep 385 386 TYPE(Element_t), POINTER :: BulkElement 387 INTEGER :: No, Status 388 REAL(KIND=dp) :: Coord(3),Velo(3),GradVelo(3,3) 389 TYPE(Element_t), POINTER :: BulkElement2 390 TYPE(Mesh_t), POINTER :: Mesh 391 TYPE(Valuelist_t), POINTER :: Params 392 REAL(KIND=dp) :: VeloAtPoint(3), GradVeloAtPoint(3,3),dtime 393 LOGICAL :: Stat, UseGradVelo, Visited = .FALSE. 394 INTEGER :: i,j,k,l,n,dim,TimeOrder, NewLost(3), OldLost, FixedLost 395 INTEGER, POINTER :: NodeIndexes(:), FieldPerm(:),FieldIndexes(:) 396 REAL(KIND=dp) :: SqrtElementMetric, Weight, Speed, SpeedMin 397 REAL(KIND=dp), POINTER :: Basis(:), dBasisdx(:,:), Coordinate(:,:), Velocity(:,:) 398 LOGICAL :: GotIt, SkipZeroTime 399 CHARACTER(LEN=MAX_NAME_LEN) :: VariableName 400 TYPE(Variable_t), POINTER :: VeloVar 401 TYPE(Variable_t), POINTER :: DtVar 402 403 404 SAVE :: Visited, Mesh, dim, Basis, dBasisdx, Params, VeloVar, UseGradvelo, DtVar, & 405 SpeedMin, NewLost 406 407 IF( .NOT. Visited ) THEN 408 Mesh => GetMesh() 409 dim = Mesh % MeshDim 410 n = Mesh % MaxElementNodes 411 ALLOCATE( Basis(n), dBasisdx(n, 3) ) 412 413 Params => GetSolverParams() 414 415 VariableName = ListGetString(Params,'Velocity Variable Name',Found) 416 IF(.NOT. Found) VariableName = 'flow solution' 417 VeloVar => VariableGet( Mesh % Variables, TRIM(VariableName) ) 418 IF(.NOT. ASSOCIATED( VeloVar ) ) THEN 419 CALL Fatal('ParticleFieldInteraction','Velocity field variable does not exist: '//TRIM(VariableName)) 420 END IF 421 UseGradVelo = GetLogical( Params,'Velocity Gradient Correction',Found) 422 423 IF(UseGradVelo .AND. Particles % RK2 ) THEN 424 CALL Warn('ParticlePathIntegral','Quadratic source correction incompatibe with Runge-Kutta') 425 UseGradVelo = .FALSE. 426 END IF 427 428 IF( .NOT. Particles % DtConstant ) THEN 429 DtVar => ParticleVariableGet( Particles,'particle dt') 430 IF( .NOT. ASSOCIATED( DtVar ) ) THEN 431 CALL Fatal('SetParticleVelocities','Required field > particle dt < not present!') 432 END IF 433 END IF 434 435 SpeedMin = ListGetConstReal( Params,'Particle Min Speed',Found) 436 IF(.NOT. Found) SpeedMin = EPSILON( SpeedMin ) 437 438 NewLost = 0 439 440 Visited = .TRUE. 441 END IF 442 443 Coordinate => Particles % Coordinate 444 Velocity => Particles % Velocity 445 Coord = 0.0_dp 446 Velo = 0.0_dp 447 448 OldLost = 0 449 FixedLost = 0 450 451 IF( Particles % DtConstant ) THEN 452 dtime = Particles % DtSign * Particles % dtime 453 END IF 454 455 SkipZeroTime = .NOT. ( Particles % DtConstant .OR. FirstStep ) 456 457 DO No = 1, Particles % NumberOfParticles 458 Status = GetParticleStatus( Particles, No ) 459 IF( Status >= PARTICLE_LOST .OR. & 460 Status <= PARTICLE_INITIATED .OR. & 461 Status == PARTICLE_FIXEDCOORD .OR. & 462 Status == PARTICLE_WALLBOUNDARY ) THEN 463 OldLost = OldLost + 1 464 CYCLE 465 END IF 466 467 ElementIndex = GetParticleElement( Particles, No ) 468 IF( ElementIndex == 0 ) THEN 469 Particles % Status(No) = PARTICLE_LOST 470 NewLost(1) = NewLost(1) + 1 471 CYCLE 472 END IF 473 474 ! If the particle has not moved then it cannot have 475 ! any change in the velocity. 476 IF( SkipZeroTime ) THEN 477 IF( ABS( DtVar % Values(No) ) < TINY( dtime ) ) CYCLE 478 END IF 479 480 481 BulkElement => Mesh % Elements( ElementIndex ) 482 483 Coord(1:dim) = Coordinate( No, 1:dim ) 484 485 !------------------------------------------------------------------------- 486 ! Get velocity from mesh 487 !------------------------------------------------------------------------- 488 IF( UseGradVelo ) THEN 489 stat = ParticleElementInfo( BulkElement, Coord, & 490 SqrtElementMetric, Basis, dBasisdx ) 491 ELSE 492 stat = ParticleElementInfo( BulkElement, Coord, & 493 SqrtElementMetric, Basis ) 494 END IF 495 496 IF(.NOT. Stat ) THEN 497 Particles % Status(No) = PARTICLE_LOST 498 NewLost(2) = NewLost(2) + 1 499 CYCLE 500 END IF 501 502 IF( UseGradVelo ) THEN 503 CALL GetVectorFieldInMesh(VeloVar,BulkElement, Basis, VeloAtPoint, & 504 dBasisdx, GradVeloAtPoint ) 505 IF( .NOT. Particles % DtConstant ) THEN 506 dtime = Particles % DtSign * DtVar % Values(No) 507 END IF 508 DO i=1,dim 509 Velo(i) = VeloAtPoint(i) + & 510 0.5_dp * SUM( GradVeloAtPoint(i,1:dim) * VeloAtPoint(1:dim) ) * dtime 511 END DO 512 ELSE 513 CALL GetVectorFieldInMesh(VeloVar, BulkElement, Basis, VeloAtPoint ) 514 Velo(1:dim) = VeloAtPoint(1:dim) 515 END IF 516 517 Speed = SQRT( SUM( Velo(1:dim) ** 2 ) ) 518 IF( Speed < SpeedMin ) THEN 519 Particles % Status(No) = PARTICLE_FIXEDCOORD 520 Velocity( No, 1:dim ) = 0.0_dp 521 FixedLost = FixedLost + 1 522 ELSE 523 Velocity( No, 1:dim ) = Velo(1:dim) 524 END IF 525 END DO 526 527 IF( .FALSE. ) THEN 528 IF( NewLost(1) > 0 ) CALL Warn('SetParticleVelocities','Some particles could not be located') 529 PRINT *,'Total number of particles:',Particles % NumberOfParticles 530 PRINT *,'Passive particles:',OldLost 531 PRINT *,'New lost particles:',NewLost 532 PRINT *,'New fixed velo particles:',FixedLost 533 END IF 534 535 536 END SUBROUTINE SetParticleVelocities 537 538 539 540 541 !------------------------------------------------------------------------ 542 !> Compute field values at the given points in the FE mesh. 543 !------------------------------------------------------------------------- 544 SUBROUTINE SetAdvectedField() 545 546 TYPE(Element_t), POINTER :: BulkElement 547 INTEGER :: No, Status, NoParticles 548 REAL(KIND=dp) :: dtime, Coord(3),Velo(3),val,vals(10) 549 TYPE(Mesh_t), POINTER :: Mesh 550 TYPE(Valuelist_t), POINTER :: Params 551 LOGICAL :: Stat, Visited = .FALSE. 552 INTEGER :: i,j,k,l,n,nsize,dim,wallcount,NoVar,NoNorm,dofs,maxdim,VarType 553 INTEGER, POINTER :: NodeIndexes(:), PPerm(:) 554 REAL(KIND=dp) :: SqrtElementMetric, Norm, PrevNorm = 0.0_dp, Change 555 REAL(KIND=dp), POINTER :: Basis(:) 556 LOGICAL :: GotIt, Difference,Cumulative,Derivative,GotVar,GotRes,GotOper,Debug,& 557 UsePerm,InternalVariable,Initiated, Parallel 558 CHARACTER(LEN=MAX_NAME_LEN) :: VariableName, ResultName, OperName, Name 559 TYPE(Variable_t), POINTER :: TargetVar, ResultVar, DataVar, Var 560 TYPE(Variable_t), POINTER :: ParticleVar 561 REAL(KIND=dp), POINTER :: TmpValues(:), NodeValues(:), NewValues(:) 562 INTEGER, POINTER :: TmpPerm(:), UnitPerm(:) 563 REAL(KIND=dp) :: x,y,z 564 565 566 SAVE :: Visited, PrevNorm, UnitPerm 567 568 CALL Info('ParticleAdvector','Setting the advected fields',Level=10) 569 570 571 Mesh => GetMesh() 572 dim = Mesh % MeshDim 573 n = Mesh % MaxElementNodes 574 ALLOCATE( Basis(n) ) 575 Coord = 0.0_dp 576 Velo = 0.0_dp 577 578 Params => GetSolverParams() 579 NoNorm = GetInteger( Params,'Norm Variable Index',GotIt) 580 NoParticles = Particles % NumberOfParticles 581 maxdim = 0 582 583 584 DataVar => VariableGet( Mesh % Variables,'AdvectorData') 585 IF( ASSOCIATED( DataVar ) ) THEN 586 nsize = SIZE( DataVar % Values ) 587 DataVar % Output = .FALSE. 588 ELSE 589 nsize = Mesh % NumberOfNodes 590 END IF 591 592 Parallel = ( ParEnv % PEs > 1 ) 593 594 Initiated = .FALSE. 595100 NoVar = 0 596 DO WHILE(.TRUE.) 597 NoVar = NoVar + 1 598 599 WRITE (Name,'(A,I0)') 'Variable ',NoVar 600 VariableName = GetString( Params,Name,GotVar) 601 IF(.NOT. GotVar ) EXIT 602 603 CALL Info('ParticleAdvector','Setting field for variable: '//TRIM(VariableName),Level=15) 604 605 ! Get the target variables 606 ! Variables starting with 'particle' as associated with particles 607 !---------------------------------------------------------------- 608 TargetVar => NULL() 609 610 IF( VariableName == 'particle coordinate' .OR. & 611 VariableName == 'particle velocity' .OR. & 612 VariableName == 'particle force') THEN 613 614 dofs = dim 615 InternalVariable = .TRUE. 616 maxdim = MAX( dim, maxdim ) 617 ELSE IF( SEQL(VariableName, 'particle') ) THEN 618 dofs = 1 619 InternalVariable = .TRUE. 620 maxdim = MAX( 1, maxdim ) 621 ELSE 622 TargetVar => VariableGet( Mesh % Variables, TRIM(VariableName) ) 623 IF( .NOT. ASSOCIATED(TargetVar)) THEN 624 CALL Fatal('ParticleAdvector','Could not obtain target variable:'& 625 //TRIM(VariableName)) 626 CYCLE 627 END IF 628 dofs = TargetVar % dofs 629 IF( dofs > 1 ) THEN 630 CALL Fatal('ParticleAdvector','Advection implemented so far only for scalars') 631 END IF 632 InternalVariable = .FALSE. 633 maxdim = MAX( dofs, maxdim ) 634 END IF 635 636 IF(.NOT. Initiated ) CYCLE 637 638 639 ! For internal variables the target name is the field name 640 !--------------------------------------------------------- 641 Difference = .FALSE. 642 Derivative = .FALSE. 643 Cumulative = .FALSE. 644 645 WRITE (Name,'(A,I0)') 'Operator ',NoVar 646 OperName = GetString( Params,Name,GotOper) 647 648 IF( GotOper ) THEN 649 IF( OperName == 'difference' ) THEN 650 Difference = .TRUE. 651 ELSE IF( OperName == 'derivative') THEN 652 Derivative = .TRUE. 653 ELSE IF( OperName == 'cumulative' .OR. OperName == 'age' ) THEN 654 Cumulative = .TRUE. 655 ELSE 656 CALL Fatal('SetAdvectedField','Unknown operator: '//TRIM(OperName)) 657 END IF 658 ELSE 659 OperName = 'adv' 660 END IF 661 662 CALL Info('ParticleAdvector','Using operator for variable: '//TRIM(OperName),Level=15) 663 664 WRITE (Name,'(A,I0)') 'Result Variable ',NoVar 665 ResultName = GetString( Params,Name,GotRes) 666 IF( .NOT. GotRes ) THEN 667 ResultName = TRIM(OperName)//' '//TRIM(VariableName) 668 END IF 669 670 671 ! Create variables if they do not exist 672 !--------------------------------------------------------- 673 674 ResultVar => VariableGet( Mesh % Variables, TRIM(ResultName) ) 675 IF( ASSOCIATED(ResultVar) ) THEN 676 IF( ASSOCIATED( DataVar) ) THEN 677 IF( DataVar % TYPE /= ResultVar % TYPE ) THEN 678 CALL Fatal('ParticleAdvector','ResultVar is of wrong type, use new name for result variable!') 679 END IF 680 IF( SIZE( DataVar % Values ) /= SIZE( ResultVar % Values) ) THEN 681 CALL Fatal('ParticleAdvector','ResultVar is of wrong size, use new name for result variable!') 682 END IF 683 END IF 684 CALL Info('ParticleAdvector','Found a pre-existing result variable: '//TRIM(ResultName),Level=20) 685 ELSE 686 GotIt = .FALSE. 687 PPerm => NULL() 688 IF( ASSOCIATED( DataVar ) ) THEN 689 CALL Info('ParticleAdvector','Using non-nodal given permutation for data',Level=15) 690 PPerm => DataVar % Perm 691 VarType = DataVar % TYPE 692 ELSE IF( ASSOCIATED( TargetVar ) ) THEN 693 CALL Info('ParticleAdvector','Using inherited permutation for data',Level=15) 694 PPerm => TargetVar % Perm 695 VarType = TargetVar % TYPE 696 END IF 697 698 IF( .NOT. ASSOCIATED( PPerm ) ) THEN 699 IF(.NOT. ASSOCIATED(UnitPerm)) THEN 700 CALL Info('ParticleAdvector','Creating unity permutation for data',Level=15) 701 ALLOCATE( UnitPerm(nsize) ) 702 DO i=1,nsize 703 UnitPerm(i) = i 704 END DO 705 END IF 706 PPerm => UnitPerm 707 VarType = 0 708 END IF 709 710 CALL VariableAddVector( Mesh % Variables,Mesh,PSolver,ResultName,dofs,& 711 Perm = PPerm, VarType = VarType ) 712 713 IF( dofs == 1 ) THEN 714 CALL Info('ParticleAdvector','Created a scalar variable: '//TRIM(ResultName) ) 715 ELSE 716 CALL Info('ParticleAdvector','Created a vector variable: '//TRIM(ResultName) ) 717 END IF 718 ResultVar => VariableGet( Mesh % Variables, TRIM(ResultName)) 719 IF(.NOT. ASSOCIATED(ResultVar)) CALL Fatal('ParticleAdvector','Problems in VariableAdd') 720 END IF 721 722 723 ! Finally, set the values 724 !--------------------------------------------------------- 725 IF( InternalVariable ) THEN 726 CALL Info('ParticleAdvector','Setting particle variable to fields',Level=15) 727 728 IF( VariableName == 'particle coordinate') THEN 729 IF( ResultVar % Dofs /= dim ) THEN 730 CALL Fatal('ParticleAdvector','Variable should have dim dofs: '//TRIM(VariableName)) 731 END IF 732 DO i=1,NoParticles 733 DO j=1,dim 734 NewValues(dim*(i-1)+j) = Particles % Coordinate(i,j) 735 END DO 736 END DO 737 ELSE IF( VariableName == 'particle coordinate_abs') THEN 738 DO i=1,NoParticles 739 val = 0.0_dp 740 DO j=1,dim 741 val = val + Particles % Coordinate(i,j)**2 742 END DO 743 NewValues(i) = SQRT( val ) 744 END DO 745 ELSE IF( VariableName == 'particle velocity') THEN 746 IF( ResultVar % Dofs /= dim ) THEN 747 CALL Fatal('ParticleAdvector','Variable should have dim dofs: '//TRIM(VariableName)) 748 END IF 749 DO i=1,NoParticles 750 DO j=1,dim 751 NewValues(dim*(i-1)+j) = Particles % Velocity(i,j) 752 END DO 753 END DO 754 755 ELSE IF( VariableName == 'particle velocity_abs') THEN 756 DO i=1,NoParticles 757 val = 0.0_dp 758 DO j=1,dim 759 val = val + Particles % Velocity(i,j)**2 760 END DO 761 NewValues(i) = SQRT( val ) 762 END DO 763 764 ELSE IF( VariableName == 'particle force') THEN 765 IF( ResultVar % Dofs /= dim ) THEN 766 CALL Fatal('ParticleAdvector','Variable should have dim dofs: '//TRIM(VariableName)) 767 END IF 768 DO i=1,NoParticles 769 DO j=1,dim 770 NewValues(dim*(i-1)+j) = Particles % Force(i,j) 771 END DO 772 END DO 773 774 ELSE IF( VariableName == 'particle status') THEN 775 DO i=1,NoParticles 776 NewValues(i) = 1.0_dp * Particles % Status(i) 777 END DO 778 779 ELSE IF( VariableName == 'particle number') THEN 780 DO i=1,NoParticles 781 NewValues(i) = 1.0_dp * i 782 END DO 783 784 ELSE IF( VariableName == 'particle index') THEN 785 DO i=1,NoParticles 786 NewValues(i) = 1.0_dp * Particles % NodeIndex(i) 787 END DO 788 789 ELSE IF( SEQL(VariableName, 'particle') ) THEN 790 ParticleVar => ParticleVariableGet( Particles, VariableName ) 791 IF( ASSOCIATED( ParticleVar ) ) THEN 792 NewValues = ParticleVar % Values(1:SIZE(NewValues)) 793 ELSE 794 CALL Warn('ParticleAdvector','Field does not exist: '//TRIM(VariableName)) 795 END IF 796 END IF 797 798 ELSE 799 CALL Info('ParticleAdvector','Setting field variable to advected fields',Level=15) 800 801 DO i = 1, NoParticles 802 Status = GetParticleStatus( Particles, i ) 803 804 IF( Status >= PARTICLE_LOST ) CYCLE 805 IF( Status <= PARTICLE_INITIATED ) CYCLE 806 807 ElementIndex = GetParticleElement( Particles, i ) 808 809 IF( ElementIndex == 0 ) CYCLE 810 811 BulkElement => Mesh % Elements( ElementIndex ) 812 813 Coord(1:dim) = Particles % Coordinate(i, 1:dim) 814 Velo(1:dim) = Particles % Velocity(i, 1:dim) 815 816 stat = ParticleElementInfo( BulkElement, Coord, SqrtElementMetric, Basis ) 817 IF(.NOT. stat) CYCLE 818 819 IF( dofs == 1 ) THEN 820 CALL GetScalarFieldInMesh(TargetVar, BulkElement, Basis, val ) 821 NewValues( i ) = val 822 ELSE 823 CALL GetVectorFieldInMesh(TargetVar, BulkElement, Basis, vals ) 824 DO j=1,dofs 825 NewValues( dofs*(i-1)+j ) = vals(j) 826 END DO 827 END IF 828 END DO 829 END IF 830 831 ! In a serial case the nodes and particles are directly associated. 832 ! In a parallel case we need to transfer the values from particles in 833 ! different partitions to nodes. 834 !--------------------------------------------------------------------- 835 IF( Parallel ) THEN 836 NodeValues = 0._dp 837 CALL ParticleAdvectParallel( Particles, NewValues, NodeValues, dofs ) 838 END IF 839 840 ! Finally move the nodal values to the target variable 841 !--------------------------------------------------------------------- 842 IF( ASSOCIATED( DataVar ) ) THEN 843 IF( Difference .OR. Derivative ) THEN 844 ResultVar % Values = NodeValues - TargetVar % Values 845 ELSE IF( Cumulative ) THEN 846 ResultVar % Values = NodeValues + ResultVar % Values 847 ELSE 848 ResultVar % Values = NodeValues 849 END IF 850 ELSE 851 DO j=1,nsize 852 k = j 853 IF( ASSOCIATED( ResultVar % Perm ) ) k = ResultVar % Perm( k ) 854 IF( k == 0 ) CYCLE 855 IF( Difference .OR. Derivative ) THEN 856 ResultVar % Values( k ) = NodeValues( j ) - TargetVar % Values( k ) 857 ELSE IF( Cumulative ) THEN 858 ResultVar % Values( k ) = NodeValues( j ) + ResultVar % Values( k ) 859 ELSE 860 ResultVar % Values( k ) = NodeValues( j ) 861 END IF 862 END DO 863 END IF 864 865 IF( Derivative ) ResultVar % Values = ResultVar % Values / dertime 866 867 BLOCK 868 INTEGER :: t, LocalPerm(10) 869 REAL(KIND=DP) :: cval 870 TYPE(Element_t), POINTER :: Element 871 REAL(KIND=dp) :: DgScale 872 LOGICAL :: GotScale 873 874 IF( ResultVar % TYPE == variable_on_nodes_on_elements ) THEN 875 876 DGScale = ListGetCReal( Solver % Values,'DG Nodes Scale',GotScale ) 877 IF(.NOT. GotScale ) DgScale = 1.0 / SQRT( 3.0_dp ) 878 GotScale = ( ABS( DGScale - 1.0_dp ) > TINY( DgScale ) ) 879 880 IF( GotScale ) THEN 881 CALL Info('ParticleAdvector','Expanding shrinked DG field',Level=12) 882 DO t=1, Mesh % NumberOfBulkElements 883 Element => Mesh % Elements(t) 884 n = Element % TYPE % NumberOfNodes 885 LocalPerm(1:n) = ResultVar % Perm( Element % DGIndexes ) 886 IF( ANY( LocalPerm(1:n) == 0) ) CYCLE 887 cval = SUM( ResultVar % Values( LocalPerm(1:n) ) ) / n 888 DO i=1,n 889 j = LocalPerm(i) 890 ResultVar % Values(j) = cval + ( ResultVar % Values(j) - cval ) * ( 1.0_dp / DgScale ) 891 END DO 892 END DO 893 END IF 894 END IF 895 END BLOCK 896 897 ! To allow computation of change in the standard manner the Variable 898 ! is set to point to the one of interest. This is mainly used in the 899 ! tests, or could be used in for convergence monitoring also. 900 !--------------------------------------------------------------- 901 IF( NoVar == NoNorm ) THEN 902 n = SIZE( ResultVar % Values ) 903 Norm = SQRT( SUM( ResultVar % Values ** 2) / n ) 904 Change = 2.0 * ABS( Norm-PrevNorm ) / ( Norm + PrevNorm ) 905 PrevNorm = Norm 906 907 Solver % Variable % Norm = Norm 908 Solver % Variable % NonlinChange = Change 909 Solver % Variable % Values = Norm 910 911 ! Here the name is ComputeChange in order to get the change also to ElmerGUI 912 ! albeit in a very dirt style. One could also edit ElmerGUI.... 913 WRITE( Message, '(a,g15.8,g15.8,a)') & 914 'SS (ITER=1) (NRM,RELC): (',Norm, Change,' ) :: '//TRIM( ResultName ) 915 CALL Info( 'ComputeChange', Message, Level=3 ) 916 END IF 917 END DO 918 919 ! Allocate the local new temporal values 920 IF(.NOT. Initiated ) THEN 921 CALL Info('ParticleAdvector','Allocating for temporal value vectors',Level=15) 922 NoVar = NoVar - 1 923 IF( NoVar < 1 ) THEN 924 CALL Fatal('ParticleAdvector','No target and result variables exist!') 925 END IF 926 ALLOCATE( NewValues( maxdim * Particles % NumberOfParticles ) ) 927 NewValues = 0.0_dp 928 IF( Parallel ) THEN 929 ALLOCATE( NodeValues( maxdim * nsize ) ) 930 NodeValues = 0.0_dp 931 ELSE 932 NodeValues => NewValues 933 END IF 934 Initiated = .TRUE. 935 GOTO 100 936 END IF 937 938 DEALLOCATE( NewValues ) 939 IF( Parallel ) DEALLOCATE( NodeValues ) 940 Visited = .TRUE. 941 942 END SUBROUTINE SetAdvectedField 943 944!------------------------------------------------------------------------------ 945END SUBROUTINE ParticleAdvector 946!------------------------------------------------------------------------------ 947 948 949!------------------------------------------------------------------------------ 950!> Initialization for the primary solver: ParticleAdvector 951!> \ingroup Solvers 952!------------------------------------------------------------------------------ 953SUBROUTINE ParticleAdvector_Init( Model,Solver,dt,TransientSimulation ) 954 955 USE DefUtils 956 USE Interpolation 957 USE MeshUtils 958 USE ElementUtils 959 USE ParticleUtils 960 961 IMPLICIT NONE 962!------------------------------------------------------------------------------ 963 TYPE(Solver_t), TARGET :: Solver 964 TYPE(Model_t) :: Model 965 REAL(KIND=dp) :: dt 966 LOGICAL :: TransientSimulation 967!------------------------------------------------------------------------------ 968! Local variables 969!------------------------ 970 971 TYPE(ValueList_t), POINTER :: Params 972 LOGICAL :: Found, AdvectElemental, AdvectDG, AdvectIp 973 INTEGER :: NormInd 974 975 Params => GetSolverParams() 976 977 ! These are default setting that make the operation of the advection solver 978 ! possible. There should always be one passive particle for each active node. 979 !--------------------------------------------------------------------------- 980 AdvectElemental = ListGetLogical( Params,'Advect Elemental',Found) 981 AdvectDG = ListGetLogical( Params,'Advect DG',Found) 982 AdvectIp = ListGetLogical( Params,'Advect Ip',Found) 983 984 IF( AdvectElemental .OR. AdvectDg .OR. AdvectIp ) THEN 985 IF( AdvectElemental ) THEN 986 CALL ListAddString( Params,'Exported Variable 1','-elem AdvectorData') 987 ELSE IF( AdvectDg ) THEN 988 CALL ListAddString( Params,'Exported Variable 1','-dg AdvectorData') 989 ELSE 990 CALL ListAddString( Params,'Exported Variable 1','-ip AdvectorData') 991 END IF 992 CALL ListAddString( Params,'Coordinate Initialization Method','advector') 993 CALL ListAddString( Params,'Velocity Initialization Method','advector') 994 ELSE 995 CALL ListAddString( Params,'Coordinate Initialization Method','nodal ordered') 996 CALL ListAddString( Params,'Velocity Initialization Method','nodal velocity') 997 CALL ListAddConstReal( Params,'Particle Node Fraction',1.0_dp) 998 END IF 999 1000 CALL ListAddInteger( Params,'Time Order',0 ) 1001 CALL ListAddNewLogical( Params,'Particle Accurate At Face',.FALSE.) 1002 CALL ListAddLogical( Params,'Particle Dt Negative',.TRUE.) 1003 CALL ListAddLogical( Params,'Particle Fix Frozen',.TRUE.) 1004 1005 ! If we want to show a pseudonorm add a variable for which the norm 1006 ! is associated with. 1007 NormInd = ListGetInteger( Params,'Norm Variable Index',Found) 1008 IF( NormInd > 0 ) THEN 1009 IF( .NOT. ListCheckPresent( Params,'Variable') ) THEN 1010 CALL ListAddString( Solver % Values,'Variable','-nooutput -global particleadvector_var') 1011 END IF 1012 END IF 1013 1014 CALL ListAddNewLogical( Params,'No Matrix',.TRUE.) 1015 1016END SUBROUTINE ParticleAdvector_Init 1017