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 library is free software; you can redistribute it and/or 8! * modify it under the terms of the GNU Lesser General Public 9! * License as published by the Free Software Foundation; either 10! * version 2.1 of the License, or (at your option) any later version. 11! * 12! * This library 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 GNU 15! * Lesser General Public License for more details. 16! * 17! * You should have received a copy of the GNU Lesser General Public 18! * License along with this library (in file ../LGPL-2.1); if not, write 19! * to the Free Software Foundation, Inc., 51 Franklin Street, 20! * Fifth Floor, Boston, MA 02110-1301 USA 21! * 22! *****************************************************************************/ 23! 24!/****************************************************************************** 25! * 26! * ELMER/FEM Solver main program 27! * 28! ****************************************************************************** 29! * 30! * Authors: Juha Ruokolainen 31! * Email: Juha.Ruokolainen@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: 02 Jun 1997 38! * 39! *****************************************************************************/ 40 41!> \defgroup Solvers Dynamically linked solvers 42 43!> \defgroup UDF Dynamically linked functions 44 45!> \defgroup Programs Utility programs 46 47!> \degroup ElmerLib Elmer library routines 48 49!> \ingroup ElmerLib 50!> \{ 51 52#include "../config.h" 53 54!------------------------------------------------------------------------------ 55!> The main program for Elmer. Solves the equations as defined by the input files. 56!------------------------------------------------------------------------------ 57 SUBROUTINE ElmerSolver(initialize) 58!------------------------------------------------------------------------------ 59 60 USE Lists 61 USE MainUtils 62 63!------------------------------------------------------------------------------ 64 IMPLICIT NONE 65!------------------------------------------------------------------------------ 66 67 INTEGER :: Initialize 68 69!------------------------------------------------------------------------------ 70! Local variables 71!------------------------------------------------------------------------------ 72 73 INTEGER :: i,j,k,n,l,t,k1,k2,iter,Ndeg,istat,nproc,tlen,nthreads 74 CHARACTER(LEN=MAX_STRING_LEN) :: threads, CoordTransform 75 76 REAL(KIND=dp) :: s,dt,dtfunc 77 REAL(KIND=dP), POINTER :: WorkA(:,:,:) => NULL() 78 REAL(KIND=dp), POINTER, SAVE :: sTime(:), sStep(:), sInterval(:), sSize(:), & 79 steadyIt(:),nonlinIt(:),sPrevSizes(:,:),sPeriodic(:),sScan(:),sSweep(:),sPar(:) 80 81 TYPE(Element_t),POINTER :: CurrentElement 82 83 LOGICAL :: GotIt,Transient,Scanning,LastSaved, MeshMode = .FALSE. 84 85 INTEGER :: TimeIntervals,interval,timestep, & 86 TotalTimesteps,SavedSteps,CoupledMaxIter,CoupledMinIter 87 88 INTEGER, POINTER, SAVE :: Timesteps(:),OutputIntervals(:) => NULL(), ActiveSolvers(:) 89 REAL(KIND=dp), POINTER, SAVE :: TimestepSizes(:,:) 90 91 INTEGER(KIND=AddrInt) :: ControlProcedure 92 93 LOGICAL :: InitDirichlet, ExecThis 94 95 TYPE(ElementType_t),POINTER :: elmt 96 97 TYPE(ParEnv_t), POINTER :: ParallelEnv 98 99 CHARACTER(LEN=MAX_NAME_LEN) :: ModelName, eq, ExecCommand, ExtrudedMeshName 100 CHARACTER(LEN=MAX_STRING_LEN) :: OutputFile, PostFile, RestartFile, & 101 OutputName=' ',PostName=' ', When, OptionString 102 103 TYPE(Variable_t), POINTER :: Var 104 TYPE(Mesh_t), POINTER :: Mesh 105 TYPE(Solver_t), POINTER :: Solver 106 107 REAL(KIND=dp) :: CT0,RT0,tt 108 109 LOGICAL :: FirstLoad = .TRUE., FirstTime=.TRUE., Found 110 LOGICAL :: Silent, Version, GotModelName, FinishEarly 111 112 INTEGER :: iargc, NoArgs 113 INTEGER :: iostat, iSweep = 1, OptimIters 114 115 INTEGER :: MeshIndex 116 TYPE(Mesh_t), POINTER :: ExtrudedMesh 117 118 TYPE(Model_t), POINTER, SAVE :: Control 119 CHARACTER(LEN=MAX_NAME_LEN) :: MeshDir, MeshName 120 LOGICAL :: DoControl, GotParams 121 INTEGER :: nr,ni 122 REAL(KIND=dp), ALLOCATABLE :: rpar(:) 123 INTEGER, ALLOCATABLE :: ipar(:) 124 125#ifdef HAVE_TRILINOS 126INTERFACE 127 SUBROUTINE TrilinosCleanup() BIND(C,name='TrilinosCleanup') 128 IMPLICIT NONE 129 END SUBROUTINE TrilinosCleanup 130END INTERFACE 131#endif 132 133 ! Start the watches, store later 134 !-------------------------------- 135 RT0 = RealTime() 136 CT0 = CPUTime() 137 138 ! If parallel execution requested, initialize parallel environment: 139 !------------------------------------------------------------------ 140 IF(FirstTime) ParallelEnv => ParallelInit() 141 142 OutputPE = -1 143 IF( ParEnv % MyPe == 0 ) THEN 144 OutputPE = 0 145 END IF 146 147 IF ( FirstTime ) THEN 148 ! 149 ! Print banner to output: 150 ! ----------------------- 151 NoArgs = COMMAND_ARGUMENT_COUNT() 152 ! Info Level is always true until the model has been read! 153 ! This makes it possible to cast something 154 Silent = .FALSE. 155 Version = .FALSE. 156 157 IF( NoArgs > 0 ) THEN 158 i = 0 159 DO WHILE( i < NoArgs ) 160 i = i + 1 161 CALL GET_COMMAND_ARGUMENT(i, OptionString) 162 IF( OptionString=='-rpar' ) THEN 163 ! Followed by number of paramters + the parameter values 164 i = i + 1 165 CALL GET_COMMAND_ARGUMENT(i, OptionString) 166 READ( OptionString,*) nr 167 ALLOCATE( rpar(nr) ) 168 DO j=1,nr 169 i = i + 1 170 CALL GET_COMMAND_ARGUMENT(i, OptionString) 171 READ( OptionString,*) rpar(j) 172 END DO 173 CALL Info('MAIN','Read '//TRIM(I2S(nr))//' real parameters from command line!') 174 CALL SetRealParametersMATC(nr,rpar) 175 END IF 176 177 IF( OptionString=='-ipar' ) THEN 178 ! Followed by number of paramters + the parameter values 179 i = i + 1 180 CALL GET_COMMAND_ARGUMENT(i, OptionString) 181 READ( OptionString,*) ni 182 ALLOCATE( ipar(nr) ) 183 DO j=1,ni 184 i = i + 1 185 CALL GET_COMMAND_ARGUMENT(i, OptionString) 186 READ( OptionString,*) ipar(j) 187 END DO 188 CALL Info('MAIN','Read '//TRIM(I2S(ni))//' integer parameters from command line!') 189 CALL SetIntegerParametersMATC(ni,ipar) 190 END IF 191 192 Silent = Silent .OR. & 193 ( OptionString=='-s' .OR. OptionString=='--silent' ) 194 Version = Version .OR. & 195 ( OptionString=='-v' .OR. OptionString == '--version' ) 196 END DO 197 END IF 198 199 ! Set number of OpenMP threads 200 nthreads = 1 201 !$ nthreads = omp_get_max_threads() 202 IF (nthreads > 1 ) THEN 203 ! Check if OMP_NUM_THREADS environment variable is set 204 CALL envir( 'OMP_NUM_THREADS', threads, tlen ) 205 206 IF (tlen==0) THEN 207 CALL Info('MAIN','OMP_NUM_THREADS not set. Using only 1 thread per task.',Level=6) 208 nthreads = 1 209 ! Set number of threads to 1 210 !$ CALL omp_set_num_threads(nthreads) 211#ifdef HAVE_MKL 212 CALL mkl_set_num_threads(nthreads) 213#endif 214 END IF 215 END IF 216 217 ParEnv % NumberOfThreads = nthreads 218 219 IF( .NOT. Silent ) THEN 220 CALL Info( 'MAIN', ' ') 221 CALL Info( 'MAIN', '=============================================================') 222 CALL Info( 'MAIN', 'ElmerSolver finite element software, Welcome! ') 223 CALL Info( 'MAIN', 'This program is free software licensed under (L)GPL ') 224 CALL Info( 'MAIN', 'Copyright 1st April 1995 - , CSC - IT Center for Science Ltd.') 225 CALL Info( 'MAIN', 'Webpage http://www.csc.fi/elmer, Email elmeradm@csc.fi ') 226 CALL Info( 'MAIN', 'Version: ' // GetVersion() // ' (Rev: ' // GetRevision() // & 227 ', Compiled: ' // GetCompilationDate() // ')' ) 228 229 IF ( ParEnv % PEs > 1 ) THEN 230 CALL Info( 'MAIN', ' Running in parallel using ' // & 231 TRIM(i2s(ParEnv % PEs)) // ' tasks.') 232 ELSE 233 CALL Info('MAIN', ' Running one task without MPI parallelization.',Level=10) 234 END IF 235 236 ! Print out number of threads in use 237 IF ( nthreads > 1 ) THEN 238 CALL Info('MAIN', ' Running in parallel with ' // & 239 TRIM(i2s(nthreads)) // ' threads per task.') 240 ELSE 241 CALL Info('MAIN', ' Running with just one thread per task.',Level=10) 242 END IF 243 244#ifdef HAVE_FETI4I 245 CALL Info( 'MAIN', ' FETI4I library linked in.') 246#endif 247#ifdef HAVE_HYPRE 248 CALL Info( 'MAIN', ' HYPRE library linked in.') 249#endif 250#ifdef HAVE_TRILINOS 251 CALL Info( 'MAIN', ' Trilinos library linked in.') 252#endif 253#ifdef HAVE_MUMPS 254 CALL Info( 'MAIN', ' MUMPS library linked in.') 255#endif 256#ifdef HAVE_CHOLMOD 257 CALL Info( 'MAIN', ' CHOLMOD library linked in.') 258#endif 259#ifdef HAVE_SUPERLU 260 CALL Info( 'MAIN', ' SUPERLU library linked in.') 261#endif 262#ifdef HAVE_PARDISO 263 CALL Info( 'MAIN', ' PARDISO library linked in.') 264#endif 265#ifdef HAVE_MKL 266 CALL Info( 'MAIN', ' Intel MKL linked in.' ) 267#endif 268#ifdef HAVE_LUA 269 CALL Info( 'MAIN', ' Lua interpreted linked in.' ) 270#endif 271#ifdef HAVE_ZOLTAN 272 CALL Info( 'MAIN', ' Zoltan library linked in.' ) 273#endif 274 CALL Info( 'MAIN', '=============================================================') 275 END IF 276 277 IF( Version ) RETURN 278 279 CALL InitializeElementDescriptions() 280 END IF 281 282 ! Read input file name either as an argument, or from the default file: 283 !---------------------------------------------------------------------- 284 GotModelName = .FALSE. 285 IF ( NoArgs > 0 ) THEN 286 CALL GET_COMMAND_ARGUMENT(1, ModelName) 287 IF( ModelName(1:1) /= '-') THEN 288 GotModelName = .TRUE. 289 IF (NoArgs > 1) CALL GET_COMMAND_ARGUMENT(2, eq) 290 END IF 291 END IF 292 293 IF( .NOT. GotModelName ) THEN 294 OPEN( 1, File='ELMERSOLVER_STARTINFO', STATUS='OLD', IOSTAT=iostat ) 295 IF( iostat /= 0 ) THEN 296 CALL Fatal( 'ElmerSolver', 'Unable to find ELMERSOLVER_STARTINFO, can not execute.' ) 297 END IF 298 READ(1,'(a)') ModelName 299 CLOSE(1) 300 END IF 301 302 ! This sets optionally some internal parameters for doing scanning 303 ! over a parameter space / optimization. 304 !----------------------------------------------------------------- 305 IF( FirstTime ) THEN 306 OPEN( Unit=InFileUnit, Action='Read',File=ModelName,Status='OLD',IOSTAT=iostat) 307 IF( iostat /= 0 ) THEN 308 CALL Fatal( 'ElmerSolver', 'Unable to find input file [' // & 309 TRIM(Modelname) // '], can not execute.' ) 310 END IF 311 ALLOCATE( Control ) 312 ! Read only the "Run Control" section of the sif file. 313 CALL LoadInputFile( Control,InFileUnit,ModelName,MeshDir,MeshName, & 314 .FALSE., .TRUE., ControlOnly = .TRUE.) 315 DoControl = ASSOCIATED( Control % Control ) 316 IF( DoControl ) THEN 317 CALL Info('ElmerSolver','Run Control section active!') 318 OptimIters = ListGetInteger( Control % Control,'Run Control Iterations', Found ) 319 IF(.NOT. Found) OptimIters = 1 320 321 ! If there are no parameters this does nothing 322 CALL ControlParameters(Control % Control,1,GotParams,FinishEarly) 323 ELSE 324 OptimIters = 1 325 END IF 326 END IF 327 328!------------------------------------------------------------------------------ 329! Read element definition file, and initialize element types 330!------------------------------------------------------------------------------ 331 IF ( Initialize==1 ) THEN 332 CALL FreeModel(CurrentModel) 333 FirstLoad=.TRUE. 334 END IF 335 336!------------------------------------------------------------------------------ 337! Read Model and mesh from Elmer mesh data base 338!------------------------------------------------------------------------------ 339 MeshIndex = 0 340 DO WHILE( .TRUE. ) 341 342 IF ( initialize==2 ) GOTO 1 343 344 IF(MeshMode) THEN 345 CALL FreeModel(CurrentModel) 346 MeshIndex = MeshIndex + 1 347 FirstLoad = .TRUE. 348 END IF 349 350 IF ( FirstLoad ) THEN 351 IF( .NOT. Silent ) THEN 352 CALL Info( 'MAIN', ' ') 353 CALL Info( 'MAIN', ' ') 354 CALL Info( 'MAIN', '-------------------------------------') 355 CALL Info( 'MAIN', 'Reading Model: '//TRIM( ModelName) ) 356 END IF 357 358 INQUIRE(Unit=InFileUnit, Opened=GotIt) 359 IF ( gotIt ) CLOSE(inFileUnit) 360 361 ! Here we read the whole model including command file and detault mesh file 362 !--------------------------------------------------------------------------------- 363 OPEN( Unit=InFileUnit, Action='Read',File=ModelName,Status='OLD',IOSTAT=iostat) 364 IF( iostat /= 0 ) THEN 365 CALL Fatal( 'ElmerSolver', 'Unable to find input file [' // & 366 TRIM(Modelname) // '], can not execute.' ) 367 END IF 368 369 370 CurrentModel => LoadModel(ModelName,.FALSE.,ParEnv % PEs,ParEnv % MyPE,MeshIndex) 371 IF(.NOT.ASSOCIATED(CurrentModel)) EXIT 372 373 !---------------------------------------------------------------------------------- 374 ! Set namespace searching mode 375 !---------------------------------------------------------------------------------- 376 CALL SetNamespaceCheck( ListGetLogical( CurrentModel % Simulation, & 377 'Additive namespaces', Found ) ) 378 379 !---------------------------------------------------------------------------------- 380 MeshMode = ListGetLogical( CurrentModel % Simulation, 'Mesh Mode', Found) 381 382 !------------------------------------------------------------------------------ 383 ! Some keywords automatically require other keywords to be set 384 ! We could complain on the missing keywords later on, but sometimes 385 ! it may be just as simple to add them directly. 386 !------------------------------------------------------------------------------ 387 CALL CompleteModelKeywords( ) 388 389 !---------------------------------------------------------------------------------- 390 ! Optionally perform simple extrusion to increase the dimension of the mesh 391 !---------------------------------------------------------------------------------- 392 CALL CreateExtrudedMesh() 393 394 !---------------------------------------------------------------------------------- 395 ! If requested perform coordinate transformation directly after is has been obtained. 396 ! Don't maintain the original mesh. 397 !---------------------------------------------------------------------------------- 398 CoordTransform = ListGetString(CurrentModel % Simulation,'Coordinate Transformation',GotIt) 399 IF( GotIt ) THEN 400 CALL CoordinateTransformation( CurrentModel % Meshes, CoordTransform, & 401 CurrentModel % Simulation, .TRUE. ) 402 END IF 403 404 IF(.NOT. Silent ) THEN 405 CALL Info( 'MAIN', '-------------------------------------') 406 END IF 407 ELSE 408 IF ( Initialize==3 ) THEN 409 INQUIRE(Unit=InFileUnit, Opened=GotIt) 410 IF ( gotIt ) CLOSE(inFileUnit) 411 OPEN( Unit=InFileUnit, Action='Read', & 412 File=ModelName,Status='OLD',IOSTAT=iostat) 413 IF( iostat /= 0 ) THEN 414 CALL Fatal( 'ElmerSolver', 'Unable to find input file [' // & 415 TRIM(Modelname) // '], can not execute.' ) 416 END IF 417 END IF 418 419 IF ( .NOT.ReloadInputFile(CurrentModel) ) EXIT 420 421 Mesh => CurrentModel % Meshes 422 DO WHILE( ASSOCIATED(Mesh) ) 423 Mesh % SavesDone = 0 424 Mesh => Mesh % Next 425 END DO 426 END IF 427 4281 CONTINUE 429 430 CALL ListAddLogical( CurrentModel % Simulation, & 431 'Initialization Phase', .TRUE. ) 432 433 ! Save the start time to the list so that it may be retrieved when necessary 434 ! This could perhaps also be a global variable etc, but this will do for now. 435 !------------------------------------------------------------------------- 436 IF( ListGetLogical( CurrentModel % Simulation,'Simulation Timing',GotIt) ) THEN 437 CALL ListAddConstReal( CurrentModel % Simulation,'cputime0',ct0 ) 438 CALL ListAddConstReal( CurrentModel % Simulation,'realtime0',rt0 ) 439 END IF 440 441!------------------------------------------------------------------------------ 442! Check for transient case 443!------------------------------------------------------------------------------ 444 eq = ListGetString( CurrentModel % Simulation, 'Simulation Type', GotIt ) 445 Scanning = ( eq == 'scanning' ) 446 Transient = ( eq == 'transient' ) 447 448!------------------------------------------------------------------------------ 449! To more conveniently support the use of VTK based visualization there 450! is a hack that recognizes VTU suffix and creates a instance of output solver. 451! Note that this is really quite a dirty hack, and is not a good example. 452!----------------------------------------------------------------------------- 453 CALL AddVtuOutputSolverHack() 454 455! Support easily saving scalars to a file activated by "Scalars File" keyword. 456!----------------------------------------------------------------------------- 457 CALL AddSaveScalarsHack() 458 459!------------------------------------------------------------------------------ 460! Figure out what (flow,heat,stress,...) should be computed, and get 461! memory for the dofs 462!------------------------------------------------------------------------------ 463 CALL AddSolvers() 464 465!------------------------------------------------------------------------------ 466! Time integration and/or steady state steps 467!------------------------------------------------------------------------------ 468 CALL InitializeIntervals() 469 470!------------------------------------------------------------------------------ 471! Add coordinates and simulation time to list of variables so that 472! coordinate dependent parameter computing routines can ask for 473! them... 474!------------------------------------------------------------------------------ 475 IF ( FirstLoad ) CALL AddMeshCoordinatesAndTime() 476 477!------------------------------------------------------------------------------ 478! Get Output File Options 479!------------------------------------------------------------------------------ 480 481 ! Initial Conditions: 482 ! ------------------- 483 IF ( FirstLoad ) THEN 484 CALL SetInitialConditions() 485 486 DO i=1,CurrentModel % NumberOfSolvers 487 Solver => CurrentModel % Solvers(i) 488 IF( ListGetLogical( Solver % Values, 'Initialize Exported Variables', GotIt ) ) THEN 489 CurrentModel % Solver => Solver 490 CALL UpdateExportedVariables( Solver ) 491 END IF 492 END DO 493 END IF 494 495 ! Compute the total number of steps that will be saved to the files 496 ! Particularly look if the last step will be saved, or if it has 497 ! to be saved separately. 498 !------------------------------------------------------------------ 499 CALL CountSavedTimesteps() 500 501 CALL ListAddLogical( CurrentModel % Simulation, & 502 'Initialization Phase', .FALSE. ) 503 504 FirstLoad = .FALSE. 505 IF ( Initialize == 1 ) EXIT 506 507 ! This sets optionally some internal parameters for doing scanning 508 ! over a parameter space / optimization. 509 !----------------------------------------------------------------- 510 DO iSweep = 1, OptimIters 511 sSweep = 1.0_dp * iSweep 512 ! If there are no parameters this does nothing 513 IF( DoControl ) THEN 514 CALL ControlResetMesh(Control % Control, iSweep ) 515 IF( iSweep > 1 ) THEN 516 CALL ControlParameters(Control % Control, & 517 iSweep,GotParams,FinishEarly) 518 IF( FinishEarly ) EXIT 519 Found = ReloadInputFile(CurrentModel,RewindFile=.TRUE.) 520 CALL InitializeIntervals() 521 IF( ListGetLogical( Control % Control,'Reset Initial Conditions',Found ) ) THEN 522 CALL SetInitialConditions() 523 END IF 524 END IF 525 END IF 526 527 !------------------------------------------------------------------------------ 528 ! Here we actually perform the simulation: ExecSimulation does it all .... 529 !------------------------------------------------------------------------------ 530 ExecCommand = ListGetString( CurrentModel % Simulation, & 531 'Control Procedure', GotIt ) 532 IF ( GotIt ) THEN 533 ControlProcedure = GetProcAddr( ExecCommand ) 534 CALL ExecSimulationProc( ControlProcedure, CurrentModel ) 535 ELSE 536 CALL ExecSimulation( TimeIntervals, CoupledMinIter, & 537 CoupledMaxIter, OutputIntervals, Transient, Scanning) 538 END IF 539 540 ! This evaluates the cost function and saves the results of control 541 IF( DoControl ) THEN 542 CALL ControlParameters(Control % Control, & 543 iSweep,GotParams,FinishEarly,.TRUE.) 544 END IF 545 END DO 546 547 ! Comparison to reference is done to enable consistency test underc ctest. 548 !------------------------------------------------------------------------- 549 CALL CompareToReferenceSolution( ) 550 551 IF ( Initialize >= 2 ) EXIT 552 END DO 553 554 555 CALL CompareToReferenceSolution( Finalize = .TRUE. ) 556 557 558#ifdef DEVEL_LISTCOUNTER 559 CALL Info('ElmerSolver','Reporting list counters for code optimization purposes only!') 560 CALL Info('ElmerSolver','If you get these lines with production code undefine > DEVEL_LISTCOUNTER < !') 561 CALL ReportListCounters( CurrentModel ) 562#endif 563 564 565 566!------------------------------------------------------------------------------ 567! THIS IS THE END (...,at last, the end, my friend,...) 568!------------------------------------------------------------------------------ 569 IF ( Initialize /= 1 ) CALL Info( 'ElmerSolver', '*** Elmer Solver: ALL DONE ***',Level=3 ) 570 571 ! This may be used to study problems at the finish 572 IF( ListGetLogical( CurrentModel % Simulation,'Dirty Finish', GotIt ) ) THEN 573 CALL Info('ElmerSolver','Skipping freeing of the Model structure',Level=4) 574 RETURN 575 END IF 576 577 IF ( Initialize <= 0 ) CALL FreeModel(CurrentModel) 578 579#ifdef HAVE_TRILINOS 580 CALL TrilinosCleanup() 581#endif 582 583 IF ( FirstTime ) CALL ParallelFinalize() 584 FirstTime = .FALSE. 585 586 CALL Info('ElmerSolver','The end',Level=3) 587 588 RETURN 589 590 CONTAINS 591 592 593 ! Optionally create extruded mesh on-the-fly. 594 !-------------------------------------------------------------------- 595 SUBROUTINE CreateExtrudedMesh() 596 597 INTEGER :: ExtrudeLayers 598 599 ExtrudeLayers = GetInteger(CurrentModel % Simulation,'Extruded Mesh Levels',Found) - 1 600 IF( .NOT. Found ) THEN 601 ExtrudeLayers = GetInteger(CurrentModel % Simulation,'Extruded Mesh Layers',Found) 602 END IF 603 IF(.NOT. Found ) RETURN 604 605 IF(ExtrudeLayers < 2) THEN 606 CALL Fatal('CreateExtrudedMesh','There must be at least two layers!') 607 END IF 608 609 ExtrudedMeshName = GetString(CurrentModel % Simulation,'Extruded Mesh Name',Found) 610 IF (Found) THEN 611 ExtrudedMesh => MeshExtrude(CurrentModel % Meshes, ExtrudeLayers-1, ExtrudedMeshName) 612 ELSE 613 ExtrudedMesh => MeshExtrude(CurrentModel % Meshes, ExtrudeLayers-1) 614 END IF 615 616 ! Make the solvers point to the extruded mesh, not the original mesh 617 !------------------------------------------------------------------- 618 DO i=1,CurrentModel % NumberOfSolvers 619 IF(ASSOCIATED(CurrentModel % Solvers(i) % Mesh,CurrentModel % Meshes)) & 620 CurrentModel % Solvers(i) % Mesh => ExtrudedMesh 621 END DO 622 ExtrudedMesh % Next => CurrentModel % Meshes % Next 623 CurrentModel % Meshes => ExtrudedMesh 624 625 ! If periodic BC given, compute boundary mesh projector: 626 ! ------------------------------------------------------ 627 DO i = 1,CurrentModel % NumberOfBCs 628 IF(ASSOCIATED(CurrentModel % Bcs(i) % PMatrix)) & 629 CALL FreeMatrix( CurrentModel % BCs(i) % PMatrix ) 630 CurrentModel % BCs(i) % PMatrix => NULL() 631 k = ListGetInteger( CurrentModel % BCs(i) % Values, 'Periodic BC', GotIt ) 632 IF( GotIt ) THEN 633 CurrentModel % BCs(i) % PMatrix => PeriodicProjector( CurrentModel, ExtrudedMesh, i, k ) 634 END IF 635 END DO 636 637 END SUBROUTINE CreateExtrudedMesh 638 639 640 ! Initialize intervals for steady state, transient and scanning types. 641 !---------------------------------------------------------------------- 642 SUBROUTINE InitializeIntervals() 643 644 IF ( Transient .OR. Scanning ) THEN 645 Timesteps => ListGetIntegerArray( CurrentModel % Simulation, & 646 'Timestep Intervals', GotIt ) 647 648 IF ( .NOT.GotIt ) THEN 649 CALL Fatal( ' ', 'Keyword > Timestep Intervals < MUST be ' // & 650 'defined for transient and scanning simulations' ) 651 END IF 652 TimestepSizes => ListGetConstRealArray( CurrentModel % Simulation, & 653 'Timestep Sizes', GotIt ) 654 IF ( .NOT.GotIt ) THEN 655 IF( Scanning .OR. ListCheckPresent( CurrentModel % Simulation,'Timestep Size') ) THEN 656 ALLOCATE(TimestepSizes(SIZE(Timesteps),1)) 657 TimestepSizes = 1.0_dp 658 ELSE 659 CALL Fatal( ' ', 'Keyword [Timestep Sizes] MUST be ' // & 660 'defined for time dependent simulations' ) 661 END IF 662 END IF 663 664 CoupledMaxIter = ListGetInteger( CurrentModel % Simulation, & 665 'Steady State Max Iterations', GotIt, minv=1 ) 666 IF ( .NOT. GotIt ) CoupledMaxIter = 1 667 668 TimeIntervals = SIZE(Timesteps) 669 ELSE 670 ! Steady state 671 !------------------------------------------------------------------------------ 672 IF( .NOT. ASSOCIATED(Timesteps) ) THEN 673 ALLOCATE(Timesteps(1)) 674 END IF 675 Timesteps(1) = ListGetInteger( CurrentModel % Simulation, & 676 'Steady State Max Iterations', GotIt,minv=1 ) 677 IF ( .NOT. GotIt ) Timesteps(1) = 1 678 CoupledMaxIter = 1 679 680 ALLOCATE(TimestepSizes(1,1)) 681 TimestepSizes(1,1) = 1.0_dp 682 TimeIntervals = 1 683 END IF 684 685 686 CoupledMinIter = ListGetInteger( CurrentModel % Simulation, & 687 'Steady State Min Iterations', GotIt ) 688 IF( .NOT. GotIt ) CoupledMinIter = 1 689 690 OutputIntervals => ListGetIntegerArray( CurrentModel % Simulation, & 691 'Output Intervals', GotIt ) 692 IF( GotIt ) THEN 693 IF( SIZE(OutputIntervals) /= SIZE(TimeSteps) ) THEN 694 CALL Fatal('ElmerSolver','> Output Intervals < should have the same size as > Timestep Intervals < !') 695 END IF 696 ELSE 697 IF( .NOT. ASSOCIATED( OutputIntervals ) ) THEN 698 ALLOCATE( OutputIntervals(SIZE(TimeSteps)) ) 699 OutputIntervals = 1 700 END IF 701 END IF 702 703 704 IF ( FirstLoad ) & 705 ALLOCATE( sTime(1), sStep(1), sInterval(1), sSize(1), & 706 steadyIt(1), nonLinit(1), sPrevSizes(1,5), sPeriodic(1), & 707 sPar(1), sScan(1), sSweep(1) ) 708 709 dt = 0._dp 710 711 sTime = 0._dp 712 sStep = 0 713 sPeriodic = 0._dp 714 sScan = 0._dp 715 716 sSize = dt 717 sPrevSizes = 0_dp 718 719 sInterval = 0._dp 720 721 steadyIt = 0 722 nonlinIt = 0 723 sPar = 0 724 725 END SUBROUTINE InitializeIntervals 726 727 728 ! Calculate the number of timesteps for saving. 729 ! This is needed for some legacy formats. 730 !----------------------------------------------------------------- 731 SUBROUTINE CountSavedTimesteps() 732 733 INTEGER :: i 734 735 TotalTimesteps = 0 736 LastSaved = .FALSE. 737 DO i=1,TimeIntervals 738 DO timestep = 1,Timesteps(i) 739 IF( OutputIntervals(i) == 0 ) CYCLE 740 LastSaved = .FALSE. 741 IF ( MOD(Timestep-1, OutputIntervals(i))==0 ) THEN 742 LastSaved = .TRUE. 743 TotalTimesteps = TotalTimesteps + 1 744 END IF 745 END DO 746 END DO 747 748 DO i=1,CurrentModel % NumberOfSolvers 749 Solver => CurrentModel % Solvers(i) 750 IF(.NOT.ASSOCIATED(Solver % Variable)) CYCLE 751 IF(.NOT.ASSOCIATED(Solver % Variable % Values)) CYCLE 752 When = ListGetString( Solver % Values, 'Exec Solver', GotIt ) 753 IF ( GotIt ) THEN 754 IF ( When == 'after simulation' .OR. When == 'after all' ) THEN 755 LastSaved = .FALSE. 756 END IF 757 ELSE 758 IF ( Solver % SolverExecWhen == SOLVER_EXEC_AFTER_ALL ) THEN 759 LastSaved = .FALSE. 760 END IF 761 END IF 762 END DO 763 764 IF ( .NOT.LastSaved ) TotalTimesteps = TotalTimesteps + 1 765 IF( TotalTimesteps == 0 ) TotalTimesteps = 1 766 767 CALL Info('ElmerSolver','Number of timesteps to be saved: '//TRIM(I2S(TotalTimesteps))) 768 769 END SUBROUTINE CountSavedTimesteps 770 771 772 ! The user may request unit tests to be performed. 773 ! This will be done if any reference norm is given. 774 ! The success will be written to file TEST.PASSED as 0/1. 775 !-------------------------------------------------------- 776 SUBROUTINE CompareToReferenceSolution( Finalize ) 777 LOGICAL, OPTIONAL :: Finalize 778 779 INTEGER :: i, j, k, n, solver_id, TestCount=0, PassCount=0, FailCount, Dofs 780 REAL(KIND=dp) :: Norm, RefNorm, Tol, Err, val, refval 781 TYPE(Solver_t), POINTER :: Solver 782 TYPE(Variable_t), POINTER :: Var 783 LOGICAL :: Found, Success = .TRUE., FinalizeOnly, CompareNorm, CompareSolution, AbsoluteErr 784 CHARACTER(LEN=MAX_STRING_LEN) :: PassedMsg 785 786 SAVE TestCount, PassCount 787 788 IF( ParEnv % MyPe > 0 ) RETURN 789 790 ! Write the success to a file for further use e.g. by cmake 791 !---------------------------------------------------------- 792 FinalizeOnly = .FALSE. 793 IF( PRESENT( Finalize ) ) FinalizeOnly = .TRUE. 794 795 IF( FinalizeOnly ) THEN 796 797 ! Nothing tested 798 IF( TestCount == 0 ) RETURN 799 800 Success = ( PassCount == TestCount ) 801 FailCount = TestCount - PassCount 802 803 IF( Success ) THEN 804 CALL Info('CompareToReferenceSolution',& 805 'PASSED all '//TRIM(I2S(TestCount))//' tests!',Level=3) 806 ELSE 807 CALL Warn('CompareToReferenceSolution','FAILED '//TRIM(I2S(FailCount))//& 808 ' tests out of '//TRIM(I2S(TestCount))//'!') 809 END IF 810 811 IF( FinalizeOnly ) THEN 812 IF( ParEnv % MyPe == 0 ) THEN 813 IF( ParEnv % PEs > 1 ) THEN 814 ! Parallel test, add the number of tasks as a suffix 815 WRITE(PassedMsg, '("TEST.PASSED_",I0)') ParEnv % PEs 816 OPEN( 10, FILE = PassedMsg ) 817 ELSE 818 OPEN( 10, FILE = 'TEST.PASSED' ) 819 END IF 820 IF( Success ) THEN 821 WRITE( 10,'(I1)' ) 1 822 ELSE 823 WRITE( 10,'(I1)' ) 0 824 END IF 825 CALL FLUSH( 10 ) 826 CLOSE( 10 ) 827 END IF 828 END IF 829 830 RETURN 831 END IF 832 833 834 DO solver_id=1,CurrentModel % NumberOfSolvers 835 Solver => CurrentModel % Solvers(solver_id) 836 837 RefNorm = ListGetConstReal( Solver % Values,'Reference Norm', CompareNorm ) 838 CompareSolution = ListCheckPrefix( Solver % Values,'Reference Solution') 839 840 IF(.NOT. ( CompareNorm .OR. CompareSolution ) ) CYCLE 841 842 Var => Solver % Variable 843 IF( .NOT. ASSOCIATED( Var ) ) THEN 844 CALL Warn('CompareToReferenceSolution','Variable in Solver '& 845 //TRIM(I2S(i))//' not associated, cannot compare') 846 CYCLE 847 END IF 848 849 TestCount = TestCount + 1 850 Success = .TRUE. 851 852 ! Compare either to existing norm (ensures consistency) 853 ! or to existing solution (may also be used to directly verify) 854 ! Usually only either of these is given but for the sake of completeness 855 ! both may be used at the same time. 856 IF( CompareNorm ) THEN 857 Tol = ListGetConstReal( Solver % Values,'Reference Norm Tolerance', Found ) 858 IF(.NOT. Found ) Tol = 1.0d-5 859 Norm = Var % Norm 860 AbsoluteErr = ListGetLogical( Solver % Values,'Reference Norm Absolute', Found ) 861 Err = ABS( Norm - RefNorm ) 862 863 IF(.NOT. AbsoluteErr ) THEN 864 IF( RefNorm < TINY( RefNorm ) ) THEN 865 CALL Warn('CompareToReferenceSolution','Refenrece norm too small for relative error') 866 AbsoluteErr = .TRUE. 867 ELSE 868 Err = Err / RefNorm 869 END IF 870 END IF 871 872 ! Compare to given reference norm 873 IF( Err > Tol .OR. Err /= Err ) THEN 874 ! Warn only in the main core 875 IF( ParEnv % MyPe == 0 ) THEN 876 WRITE( Message,'(A,I0,A,ES15.8,A,ES15.8)') & 877 'Solver ',solver_id,' FAILED: Norm =',Norm,' RefNorm =',RefNorm 878 CALL Warn('CompareToReferenceSolution',Message) 879 END IF 880 Success = .FALSE. 881 ELSE 882 WRITE( Message,'(A,I0,A,ES15.8,A,ES15.8)') & 883 'Solver ',solver_id,' PASSED: Norm =',Norm,' RefNorm =',RefNorm 884 CALL Info('CompareToReferenceSolution',Message,Level=3) 885 END IF 886 IF( AbsoluteErr ) THEN 887 WRITE( Message,'(A,ES13.6)') 'Absolute Error to reference norm:',Err 888 ELSE 889 WRITE( Message,'(A,ES13.6)') 'Relative Error to reference norm:',Err 890 END IF 891 CALL Info('CompareToReferenceSolution',Message, Level = 3 ) 892 END IF 893 894 IF( CompareSolution ) THEN 895 Tol = ListGetConstReal( Solver % Values,'Reference Solution Tolerance', Found ) 896 IF(.NOT. Found ) Tol = 1.0d-5 897 Dofs = Var % Dofs 898 n = 0 899 RefNorm = 0.0_dp 900 Norm = 0.0_dp 901 Err = 0.0_dp 902 DO i=1,Solver % Mesh % NumberOfNodes 903 j = Var % Perm(i) 904 IF( j == 0 ) CYCLE 905 DO k=1,Dofs 906 IF( Dofs == 1 ) THEN 907 refval = ListGetRealAtNode( Solver % Values,'Reference Solution',i,Found ) 908 ELSE 909 refval = ListGetRealAtNode( Solver % Values,'Reference Solution '//TRIM(I2S(k)),i,Found ) 910 END IF 911 IF( Found ) THEN 912 val = Var % Values( Dofs*(j-1)+k) 913 RefNorm = RefNorm + refval**2 914 Norm = Norm + val**2 915 Err = Err + (refval-val)**2 916 n = n + 1 917 END IF 918 END DO 919 END DO 920 IF( ParEnv % PEs > 1 ) CALL Warn('CompareToReferefenSolution','Not implemented in parallel!') 921 IF( n == 0 ) CALL Fatal('CompareToReferenceSolution','Could not find any reference solution') 922 RefNorm = SQRT( RefNorm / n ) 923 Norm = SQRT( Norm / n ) 924 Err = SQRT( Err / n ) 925 926 IF( Err > Tol ) THEN 927 ! Normally warning is done for every partition but this time it is the same for all 928 IF( ParEnv % MyPe == 0 ) THEN 929 WRITE( Message,'(A,I0,A,ES15.8,A,ES15.8)') & 930 'Solver ',solver_id,' FAILED: Solution = ',Norm,' RefSolution =',RefNorm 931 CALL Warn('CompareToReferenceSolution',Message) 932 WRITE( Message,'(A,ES13.6)') 'Relative Error to reference solution:',Err 933 CALL Info('CompareToReferenceSolution',Message, Level = 3 ) 934 END IF 935 Success = .FALSE. 936 ELSE 937 WRITE( Message,'(A,I0,A,ES15.8,A,ES15.8)') & 938 'Solver ',solver_id,' PASSED: Solution =',Norm,' RefSolution =',RefNorm 939 CALL Info('CompareToReferenceSolution',Message,Level=3) 940 END IF 941 END IF 942 943 IF( Success ) PassCount = PassCount + 1 944 END DO 945 946 947 END SUBROUTINE CompareToReferenceSolution 948 949 950 951 952 ! This is a dirty hack that adds an instance of ResultOutputSolver to the list of Solvers. 953 ! The idea is that it is much easier for the end user to take into use the vtu output this way. 954 ! The solver itself has limited set of parameters needed and is therefore approapriate for this 955 ! kind of hack. It can of course be also added as a regular solver also. 956 !---------------------------------------------------------------------------------------------- 957 SUBROUTINE AddVtuOutputSolverHack() 958 TYPE(Solver_t), POINTER :: ABC(:), PSolver 959 CHARACTER(LEN=MAX_NAME_LEN) :: str 960 INTEGER :: i,j,j2,j3,k,n 961 TYPE(ValueList_t), POINTER :: Params, Simu 962 LOGICAL :: Found, VtuFormat 963 INTEGER :: AllocStat 964 LOGICAL, SAVE :: Visited = .FALSE. 965 966 Simu => CurrentModel % Simulation 967 str = ListGetString( Simu,'Post File',Found) 968 IF(.NOT. Found) RETURN 969 970 k = INDEX( str,'.vtu' ) 971 VtuFormat = ( k /= 0 ) 972 973 IF(.NOT. VtuFormat ) RETURN 974 975 ! No use to create the same solver twice 976 IF( Visited ) RETURN 977 Visited = .TRUE. 978 979 CALL Info('AddVtuOutputSolverHack','Adding ResultOutputSolver to write VTU output in file: '& 980 //TRIM(str(1:k-1))) 981 982 CALL ListRemove( Simu,'Post File') 983 n = CurrentModel % NumberOfSolvers+1 984 ALLOCATE( ABC(n), STAT = AllocStat ) 985 IF( AllocStat /= 0 ) CALL Fatal('AddVtuOutputSolverHack','Allocation error 1') 986 987 CALL Info('AddVtuOutputSolverHack','Increasing number of solver to: '& 988 //TRIM(I2S(n)),Level=8) 989 DO i=1,n-1 990 ! Def_Dofs is the only allocatable structure within Solver_t: 991 IF( ALLOCATED( CurrentModel % Solvers(i) % Def_Dofs ) ) THEN 992 j = SIZE(CurrentModel % Solvers(i) % Def_Dofs,1) 993 j2 = SIZE(CurrentModel % Solvers(i) % Def_Dofs,2) 994 j3 = SIZE(CurrentModel % Solvers(i) % Def_Dofs,3) 995 ALLOCATE( ABC(i) % Def_Dofs(j,j2,j3), STAT = AllocStat ) 996 IF( AllocStat /= 0 ) CALL Fatal('AddVtuOutputSolverHack','Allocation error 2') 997 END IF 998 999 ! Copy the content of the Solver structure 1000 ABC(i) = CurrentModel % Solvers(i) 1001 1002 ! Nullify the old structure since otherwise bad things may happen at deallocation 1003 NULLIFY( CurrentModel % Solvers(i) % ActiveElements ) 1004 NULLIFY( CurrentModel % Solvers(i) % Mesh ) 1005 NULLIFY( CurrentModel % Solvers(i) % BlockMatrix ) 1006 NULLIFY( CurrentModel % Solvers(i) % Matrix ) 1007 NULLIFY( CurrentModel % Solvers(i) % Variable ) 1008 END DO 1009 1010 ! Deallocate the old structure and set the pointer to the new one 1011 DEALLOCATE( CurrentModel % Solvers ) 1012 CurrentModel % Solvers => ABC 1013 CurrentModel % NumberOfSolvers = n 1014 1015 ! Now create the ResultOutputSolver instance on-the-fly 1016 CurrentModel % Solvers(n) % PROCEDURE = 0 1017 NULLIFY( CurrentModel % Solvers(n) % Matrix ) 1018 NULLIFY( CurrentModel % Solvers(n) % BlockMatrix ) 1019 NULLIFY( CurrentModel % Solvers(n) % Variable ) 1020 NULLIFY( CurrentModel % Solvers(n) % ActiveElements ) 1021 CurrentModel % Solvers(n) % NumberOfActiveElements = 0 1022 j = CurrentModel % NumberOfBodies 1023 ALLOCATE( CurrentModel % Solvers(n) % Def_Dofs(10,j,6),STAT=AllocStat) 1024 IF( AllocStat /= 0 ) CALL Fatal('AddVtuOutputSolverHack','Allocation error 3') 1025 CurrentModel % Solvers(n) % Def_Dofs = -1 1026 CurrentModel % Solvers(n) % Def_Dofs(:,:,1) = 1 1027 1028 ! Add some keywords to the list 1029 CurrentModel % Solvers(n) % Values => ListAllocate() 1030 Params => CurrentModel % Solvers(n) % Values 1031 CALL ListAddString( Params,& 1032 'Procedure', 'ResultOutputSolve ResultOutputSolver',.FALSE.) 1033 CALL ListAddString(Params,'Equation','InternalVtuOutputSolver') 1034 CALL ListAddString(Params,'Output Format','vtu') 1035 CALL ListAddString(Params,'Output File Name',str(1:k-1),.FALSE.) 1036 CALL ListAddString(Params,'Exec Solver','after saving') 1037 CALL ListAddLogical(Params,'Save Geometry IDs',.TRUE.) 1038 CALL ListAddLogical(Params,'Check Simulation Keywords',.TRUE.) 1039 1040 ! Add a few often needed keywords also if they are given in simulation section 1041 CALL ListCopyPrefixedKeywords( Simu, Params, 'vtu:' ) 1042 1043 CALL Info('AddVtuOutputSolverHack','Finished appending VTU output solver',Level=12) 1044 1045 END SUBROUTINE AddVtuOutputSolverHack 1046 1047 1048 ! This is a dirty hack that adds an instance of SaveScalars to the list of Solvers. 1049 ! The idea is that it is much easier for the end user to add a basic instance. 1050 !---------------------------------------------------------------------------------------------- 1051 SUBROUTINE AddSaveScalarsHack() 1052 TYPE(Solver_t), POINTER :: ABC(:), PSolver 1053 CHARACTER(LEN=MAX_NAME_LEN) :: str 1054 INTEGER :: i,j,j2,j3,k,n 1055 TYPE(ValueList_t), POINTER :: Params, Simu 1056 LOGICAL :: Found, VtuFormat 1057 INTEGER :: AllocStat 1058 LOGICAL, SAVE :: Visited = .FALSE. 1059 1060 Simu => CurrentModel % Simulation 1061 str = ListGetString( Simu,'Scalars File',Found) 1062 IF(.NOT. Found) RETURN 1063 1064 ! No use to create the same solver twice 1065 IF( Visited ) RETURN 1066 Visited = .TRUE. 1067 1068 CALL Info('AddSaveScalarsHack','Adding SaveScalars solver to write scalars into file: '& 1069 //TRIM(str)) 1070 1071 n = CurrentModel % NumberOfSolvers+1 1072 ALLOCATE( ABC(n), STAT = AllocStat ) 1073 IF( AllocStat /= 0 ) CALL Fatal('AddSaveScalarsHack','Allocation error 1') 1074 1075 CALL Info('AddSaveScalarsHack','Increasing number of solver to: '& 1076 //TRIM(I2S(n)),Level=8) 1077 DO i=1,n-1 1078 ! Def_Dofs is the only allocatable structure within Solver_t: 1079 IF( ALLOCATED( CurrentModel % Solvers(i) % Def_Dofs ) ) THEN 1080 j = SIZE(CurrentModel % Solvers(i) % Def_Dofs,1) 1081 j2 = SIZE(CurrentModel % Solvers(i) % Def_Dofs,2) 1082 j3 = SIZE(CurrentModel % Solvers(i) % Def_Dofs,3) 1083 ALLOCATE( ABC(i) % Def_Dofs(j,j2,j3), STAT = AllocStat ) 1084 IF( AllocStat /= 0 ) CALL Fatal('AddVtuOutputSolverHack','Allocation error 2') 1085 END IF 1086 1087 ! Copy the content of the Solver structure 1088 ABC(i) = CurrentModel % Solvers(i) 1089 1090 ! Nullify the old structure since otherwise bad things may happen at deallocation 1091 NULLIFY( CurrentModel % Solvers(i) % ActiveElements ) 1092 NULLIFY( CurrentModel % Solvers(i) % Mesh ) 1093 NULLIFY( CurrentModel % Solvers(i) % BlockMatrix ) 1094 NULLIFY( CurrentModel % Solvers(i) % Matrix ) 1095 NULLIFY( CurrentModel % Solvers(i) % Variable ) 1096 END DO 1097 1098 ! Deallocate the old structure and set the pointer to the new one 1099 DEALLOCATE( CurrentModel % Solvers ) 1100 CurrentModel % Solvers => ABC 1101 CurrentModel % NumberOfSolvers = n 1102 1103 ! Now create the ResultOutputSolver instance on-the-fly 1104 CurrentModel % Solvers(n) % PROCEDURE = 0 1105 NULLIFY( CurrentModel % Solvers(n) % Matrix ) 1106 NULLIFY( CurrentModel % Solvers(n) % BlockMatrix ) 1107 NULLIFY( CurrentModel % Solvers(n) % Variable ) 1108 NULLIFY( CurrentModel % Solvers(n) % ActiveElements ) 1109 CurrentModel % Solvers(n) % NumberOfActiveElements = 0 1110 j = CurrentModel % NumberOfBodies 1111 ALLOCATE( CurrentModel % Solvers(n) % Def_Dofs(10,j,6),STAT=AllocStat) 1112 IF( AllocStat /= 0 ) CALL Fatal('AddSaveScalarsHack','Allocation error 3') 1113 CurrentModel % Solvers(n) % Def_Dofs = -1 1114 CurrentModel % Solvers(n) % Def_Dofs(:,:,1) = 1 1115 1116 ! Add some keywords to the list 1117 CurrentModel % Solvers(n) % Values => ListAllocate() 1118 Params => CurrentModel % Solvers(n) % Values 1119 CALL ListAddString( Params,'Procedure', 'SaveData SaveScalars',.FALSE.) 1120 CALL ListAddString(Params,'Equation','InternalSaveScalars') 1121 CALL ListAddString(Params,'Filename',TRIM(str),.FALSE.) 1122 CALL ListAddString(Params,'Exec Solver','after saving') 1123 1124 ! Add a few often needed keywords also if they are given in simulation section 1125 CALL ListCopyPrefixedKeywords( Simu, Params, 'scalars:' ) 1126 1127 CALL Info('AddSaveScalarsHack','Finished appending SaveScalars solver',Level=12) 1128 1129 END SUBROUTINE AddSaveScalarsHack 1130 1131 1132 1133!------------------------------------------------------------------------------ 1134!> Adds flags for active solvers. 1135!------------------------------------------------------------------------------ 1136 SUBROUTINE AddSolvers() 1137!------------------------------------------------------------------------------ 1138 INTEGER :: i,j,k,nlen 1139 LOGICAL :: InitSolver, Found 1140!------------------------------------------------------------------------------ 1141 1142 CALL Info('AddSolvers','Setting up '//TRIM(I2S(CurrentModel % NumberOfSolvers))//& 1143 ' solvers',Level=10) 1144 1145 ! This is a hack that sets Equation flags True for the "Active Solvers". 1146 ! The Equation flag is the legacy way of setting a Solver active and is still 1147 ! used internally. 1148 !---------------------------------------------------------------------------- 1149 DO i=1,CurrentModel % NumberOfSolvers 1150 1151 eq = ListGetString( CurrentModel % Solvers(i) % Values,'Equation', Found ) 1152 1153 IF ( Found ) THEN 1154 nlen = LEN_TRIM(eq) 1155 DO j=1,CurrentModel % NumberOFEquations 1156 ActiveSolvers => ListGetIntegerArray( CurrentModel % Equations(j) % Values, & 1157 'Active Solvers', Found ) 1158 IF ( Found ) THEN 1159 DO k=1,SIZE(ActiveSolvers) 1160 IF ( ActiveSolvers(k) == i ) THEN 1161 CALL ListAddLogical( CurrentModel % Equations(j) % Values, eq(1:nlen), .TRUE. ) 1162 EXIT 1163 END IF 1164 END DO 1165 END IF 1166 END DO 1167 END IF 1168 END DO 1169 1170 ! Add the dynamically linked solver to be called later 1171 !--------------------------------------------------------------------- 1172 DO i=1,CurrentModel % NumberOfSolvers 1173 eq = ListGetString( CurrentModel % Solvers(i) % Values,'Equation', Found ) 1174 CALL Info('AddSolvers','Setting up solver '//TRIM(I2S(i))//': '//TRIM(eq),Level=10) 1175 1176 Solver => CurrentModel % Solvers(i) 1177 InitSolver = ListGetLogical( Solver % Values, 'Initialize', Found ) 1178 IF ( Found .AND. InitSolver ) THEN 1179 CALL FreeMatrix( Solver % Matrix ) 1180 CALL ListAddLogical( Solver % Values, 'Initialize', .FALSE. ) 1181 END IF 1182 1183 IF ( Solver % PROCEDURE == 0 .OR. InitSolver ) THEN 1184 IF ( .NOT. ASSOCIATED( Solver % Mesh ) ) THEN 1185 Solver % Mesh => CurrentModel % Meshes 1186 END IF 1187 CurrentModel % Solver => Solver 1188 CALL AddEquationBasics( Solver, eq, Transient ) 1189 CALL AddEquationSolution( Solver, Transient ) 1190 END IF 1191 END DO 1192 1193 CALL Info('AddSolvers','Setting up solvers done',Level=12) 1194 1195!------------------------------------------------------------------------------ 1196 END SUBROUTINE AddSolvers 1197!------------------------------------------------------------------------------ 1198 1199 1200!------------------------------------------------------------------------------ 1201!> Adds coordinate and time variables to the current mesh structure. 1202!------------------------------------------------------------------------------ 1203 SUBROUTINE AddMeshCoordinatesAndTime() 1204!------------------------------------------------------------------------------ 1205 TYPE(Variable_t), POINTER :: DtVar 1206 1207 CALL Info('AddMeshCoordinatesAndTime','Setting mesh coordinates and time',Level=10) 1208 1209 NULLIFY( Solver ) 1210 1211 Mesh => CurrentModel % Meshes 1212 DO WHILE( ASSOCIATED( Mesh ) ) 1213 CALL VariableAdd( Mesh % Variables, Mesh, & 1214 Name='Coordinate 1',DOFs=1,Values=Mesh % Nodes % x ) 1215 1216 CALL VariableAdd(Mesh % Variables,Mesh, & 1217 Name='Coordinate 2',DOFs=1,Values=Mesh % Nodes % y ) 1218 1219 CALL VariableAdd(Mesh % Variables,Mesh, & 1220 Name='Coordinate 3',DOFs=1,Values=Mesh % Nodes % z ) 1221 1222 CALL VariableAdd( Mesh % Variables, Mesh, Name='Time',DOFs=1, Values=sTime ) 1223 CALL VariableAdd( Mesh % Variables, Mesh, Name='Periodic Time',DOFs=1, Values=sPeriodic ) 1224 CALL VariableAdd( Mesh % Variables, Mesh, Name='Timestep', DOFs=1, Values=sStep ) 1225 CALL VariableAdd( Mesh % Variables, Mesh, Name='Timestep size', DOFs=1, Values=sSize ) 1226 CALL VariableAdd( Mesh % Variables, Mesh, Name='Timestep interval', DOFs=1, Values=sInterval ) 1227 1228 ! Save some previous timesteps for variable timestep multistep methods 1229 DtVar => VariableGet( Mesh % Variables, 'Timestep size' ) 1230 DtVar % PrevValues => sPrevSizes 1231 1232 CALL VariableAdd( Mesh % Variables, Mesh, & 1233 Name='nonlin iter', DOFs=1, Values=nonlinIt ) 1234 CALL VariableAdd( Mesh % Variables, Mesh, & 1235 Name='coupled iter', DOFs=1, Values=steadyIt ) 1236 1237 IF( ListCheckPresentAnySolver( CurrentModel,'Scanning Loops') ) THEN 1238 CALL VariableAdd( Mesh % Variables, Mesh, Name='scan', DOFs=1, Values=sScan ) 1239 END IF 1240 1241 IF( DoControl ) THEN 1242 sSweep = 1.0_dp * iSweep 1243 CALL VariableAdd( Mesh % Variables, Mesh, Name='run', DOFs=1, Values=sSweep ) 1244 END IF 1245 1246 sPar(1) = 1.0_dp * ParEnv % MyPe 1247 CALL VariableAdd( Mesh % Variables, Mesh, Name='Partition', DOFs=1, Values=sPar ) 1248 1249 ! Add partition as a elemental field in case we have just one partition 1250 ! and have asked still for partitioning into many. 1251 IF( ParEnv % PEs == 1 .AND. ASSOCIATED( Mesh % Repartition ) ) THEN 1252 BLOCK 1253 REAL(KIND=dp), POINTER :: PartField(:) 1254 INTEGER, POINTER :: PartPerm(:) 1255 INTEGER :: i, n 1256 1257 CALL Info('AddMeshCoordinatesAndTime','Adding partitioning also as a field') 1258 1259 n = Mesh % NumberOfBulkElements 1260 NULLIFY( PartField, PartPerm ) 1261 ALLOCATE( PartField(n), PartPerm(n) ) 1262 DO i=1,n 1263 PartPerm(i) = i 1264 PartField(i) = 1.0_dp * Mesh % RePartition(i) 1265 END DO 1266 1267 CALL VariableAdd( Mesh % Variables, Mesh, Name='PartField',DOFs=1, & 1268 Values=PartField, Perm=PartPerm, TYPE=Variable_on_elements) 1269 END BLOCK 1270 END IF 1271 1272 Mesh => Mesh % Next 1273 1274 END DO 1275!------------------------------------------------------------------------------ 1276 END SUBROUTINE AddMeshCoordinatesAndTime 1277!------------------------------------------------------------------------------ 1278 1279 1280!------------------------------------------------------------------------------ 1281!> Sets initial conditions for the fields. 1282!------------------------------------------------------------------------------ 1283 SUBROUTINE SetInitialConditions() 1284!------------------------------------------------------------------------------ 1285 USE DefUtils 1286 INTEGER :: DOFs 1287 CHARACTER(LEN=MAX_NAME_LEN) :: str 1288 LOGICAL :: Found, NamespaceFound 1289 TYPE(Solver_t), POINTER :: Solver 1290 INTEGER, ALLOCATABLE :: Indexes(:) 1291 REAL(KIND=dp),ALLOCATABLE :: Work(:) 1292 1293 INTEGER :: i,j,k,l,m,vect_dof,real_dof,dim 1294 1295 REAL(KIND=dp) :: nrm(3),t1(3),t2(3),vec(3),tmp(3),udot 1296 TYPE(ValueList_t), POINTER :: BC 1297 TYPE(Nodes_t), SAVE :: Nodes 1298 LOGICAL :: nt_boundary 1299 TYPE(Element_t), POINTER :: Element 1300 TYPE(Variable_t), POINTER :: var, vect_var 1301 LOGICAL :: AnyNameSpace 1302 1303 CALL Info('SetInitialConditions','Setting up initial conditions (if any)',Level=10) 1304 1305 1306 dim = CoordinateSystemDimension() 1307 1308 IF (GetLogical(GetSimulation(),'Restart Before Initial Conditions',Found)) THEN 1309 CALL Restart 1310 CALL InitCond 1311 ELSE 1312 CALL InitCond 1313 CALL Restart 1314 END IF 1315 1316 1317!------------------------------------------------------------------------------ 1318! Make sure that initial values at boundaries are set correctly. 1319! NOTE: This overrides the initial condition setting for field variables!!!! 1320!------------------------------------------------------------------------------- 1321 InitDirichlet = ListGetLogical( CurrentModel % Simulation, & 1322 'Initialize Dirichlet Conditions', GotIt ) 1323 IF ( .NOT. GotIt ) InitDirichlet = .TRUE. 1324 1325 AnyNameSpace = ListCheckPresentAnySolver( CurrentModel,'Namespace') 1326 NamespaceFound = .FALSE. 1327 1328 vect_var => NULL() 1329 IF ( InitDirichlet ) THEN 1330 Mesh => CurrentModel % Meshes 1331 DO WHILE( ASSOCIATED(Mesh) ) 1332 ALLOCATE( Work(Mesh % MaxElementDOFs) ) 1333 CALL SetCurrentMesh( CurrentModel, Mesh ) 1334 1335 DO t = Mesh % NumberOfBulkElements + 1, & 1336 Mesh % NumberOfBulkElements + Mesh % NumberOfBoundaryElements 1337 1338 Element => Mesh % Elements(t) 1339 1340 ! Set also the current element pointer in the model structure to 1341 ! reflect the element being processed: 1342 ! --------------------------------------------------------------- 1343 CurrentModel % CurrentElement => Element 1344 n = Element % TYPE % NumberOfNodes 1345 1346 BC => GetBC() 1347 IF(.NOT.ASSOCIATED(BC)) CYCLE 1348 1349 Var => Mesh % Variables 1350 DO WHILE( ASSOCIATED(Var) ) 1351 Solver => Var % Solver 1352 IF ( .NOT. ASSOCIATED(Solver) ) Solver => CurrentModel % Solver 1353 1354 IF( AnyNameSpace ) THEN 1355 str = ListGetString( Solver % Values, 'Namespace', NamespaceFound ) 1356 IF (NamespaceFound) CALL ListPushNamespace(TRIM(str)) 1357 END IF 1358 1359 IF ( Var % DOFs <= 1 ) THEN 1360 Work(1:n) = GetReal( BC,Var % Name, gotIt ) 1361 IF ( GotIt ) THEN 1362 1363 nt_boundary = .FALSE. 1364 IF ( GetElementFamily() /= 1 ) THEN 1365 k = LEN_TRIM(var % name) 1366 vect_dof = ICHAR(Var % Name(k:k))-ICHAR('0'); 1367 IF ( vect_dof>=1 .AND. vect_dof<= 3 ) THEN 1368 nt_boundary = GetLogical( BC, & 1369 'normal-tangential '//var % name(1:k-2), gotIt) 1370 1371 IF ( nt_boundary ) THEN 1372 nt_boundary = .FALSE. 1373 Vect_var => Mesh % Variables 1374 DO WHILE( ASSOCIATED(Vect_var) ) 1375 IF ( Vect_var % Dofs>=dim ) THEN 1376 DO real_dof=1,Vect_var % Dofs 1377 nt_boundary = ASSOCIATED(Var % Values, & 1378 Vect_var % Values(real_dof::Vect_var % Dofs)) 1379 IF ( nt_boundary ) EXIT 1380 END DO 1381 IF ( nt_boundary ) EXIT 1382 END IF 1383 Vect_var => Vect_var % Next 1384 END DO 1385 END IF 1386 1387 IF ( nt_boundary ) THEN 1388 CALL GetElementNodes(Nodes) 1389 nrm = NormalVector(Element,Nodes,0._dp,0._dp,.TRUE.) 1390 SELECT CASE(Element % TYPE % DIMENSION) 1391 CASE(1) 1392 t1(1) = nrm(2) 1393 t1(2) = -nrm(1) 1394 CASE(2) 1395 CALL TangentDirections(nrm,t1,t2) 1396 END SELECT 1397 1398 SELECT CASE(vect_dof) 1399 CASE(1) 1400 vec = nrm 1401 CASE(2) 1402 vec = t1 1403 CASE(3) 1404 vec = t2 1405 END SELECT 1406 END IF 1407 END IF 1408 END IF 1409 1410 DO j=1,n 1411 IF (Solver % DG) THEN 1412 BLOCK 1413 INTEGER :: i 1414 TYPE(Element_t), POINTER :: P 1415 1416 p => Element % BoundaryInfo % Left 1417 IF(.NOT.ASSOCIATED(p)) p => Element % BoundaryInfo % Right 1418 DO i=1,p % Type % NumberOfNodes 1419 IF(p % NodeIndexes(i) == Element % NodeIndexes(j) ) THEN 1420 k = p % DGIndexes(i); EXIT 1421 END IF 1422 END DO 1423 END BLOCK 1424 ELSE 1425 k = Element % NodeIndexes(j) 1426 END IF 1427 IF ( ASSOCIATED(Var % Perm) ) k = Var % Perm(k) 1428 IF ( k>0 ) THEN 1429 IF ( nt_boundary ) THEN 1430 DO l=1,dim 1431 m = l+real_dof-vect_dof 1432 tmp(l)=Vect_var % Values(Vect_var % Dofs*(k-1)+m) 1433 END DO 1434 udot = SUM(vec(1:dim)*tmp(1:dim)) 1435 tmp(1:dim)=tmp(1:dim)+(work(j)-udot)*vec(1:dim) 1436 DO l=1,dim 1437 m = l+real_dof-vect_dof 1438 Vect_var % Values(Vect_var % Dofs*(k-1)+m)=tmp(l) 1439 END DO 1440 ELSE 1441 Var % Values(k) = Work(j) 1442 END IF 1443 END IF 1444 END DO 1445 END IF 1446 1447 IF ( Transient .AND. Solver % TimeOrder==2 ) THEN 1448 Work(1:n) = GetReal( BC, TRIM(Var % Name) // ' Velocity', GotIt ) 1449 IF ( GotIt ) THEN 1450 DO j=1,n 1451 k = Element % NodeIndexes(j) 1452 IF ( ASSOCIATED(Var % Perm) ) k = Var % Perm(k) 1453 IF ( k>0 ) Var % PrevValues(k,1) = Work(j) 1454 END DO 1455 END IF 1456 Work(1:n) = GetReal( BC, TRIM(Var % Name) // ' Acceleration', GotIt ) 1457 IF ( GotIt ) THEN 1458 DO j=1,n 1459 k = Element % NodeIndexes(j) 1460 IF ( ASSOCIATED(Var % Perm) ) k = Var % Perm(k) 1461 IF ( k>0 ) Var % PrevValues(k,2) = Work(j) 1462 END DO 1463 END IF 1464 END IF 1465 ELSE 1466 CALL ListGetRealArray( BC, & 1467 Var % Name, WorkA, n, Element % NodeIndexes, gotIt ) 1468 IF ( GotIt ) THEN 1469 DO j=1,n 1470 k = Element % NodeIndexes(j) 1471 IF ( ASSOCIATED(Var % Perm) ) k = Var % Perm(k) 1472 IF(k>0) THEN 1473 DO l=1,MIN(SIZE(WorkA,1),Var % DOFs) 1474 Var % Values(Var % DOFs*(k-1)+l) = WorkA(l,1,j) 1475 END DO 1476 END IF 1477 END DO 1478 ELSE 1479 END IF 1480 END IF 1481 1482 IF(NamespaceFound) CALL ListPopNamespace() 1483 Var => Var % Next 1484 END DO 1485 END DO 1486 DEALLOCATE( Work ) 1487 Mesh => Mesh % Next 1488 END DO 1489 END IF 1490!------------------------------------------------------------------------------ 1491 END SUBROUTINE SetInitialConditions 1492!------------------------------------------------------------------------------ 1493 1494 1495!------------------------------------------------------------------------------ 1496 SUBROUTINE InitCond() 1497!------------------------------------------------------------------------------ 1498 USE DefUtils 1499 TYPE(Element_t), POINTER :: Edge 1500 INTEGER :: DOFs,i,j,k,k1,k2,l,n,m,nsize 1501 CHARACTER(LEN=MAX_NAME_LEN) :: str, VarName 1502 LOGICAL :: Found, ThingsToDO, NamespaceFound, AnyNameSpace 1503 TYPE(Solver_t), POINTER :: Solver, CSolver 1504 INTEGER, ALLOCATABLE :: Indexes(:) 1505 REAL(KIND=dp) :: Val 1506 REAL(KIND=dp),ALLOCATABLE :: Work(:) 1507 TYPE(ValueList_t), POINTER :: IC 1508 TYPE(Nodes_t) :: Nodes 1509 TYPE(GaussIntegrationPoints_t) :: IP 1510 REAL(KIND=dp), ALLOCATABLE :: Basis(:) 1511 REAL(KIND=dp) :: DetJ 1512 TYPE(ValueHandle_t) :: LocalSol_h 1513 LOGICAL :: Stat, FoundIC, PrevFoundIC 1514 INTEGER :: VarOrder, PrevBodyId 1515 !------------------------------------------------------------------------------ 1516 1517 AnyNameSpace = ListCheckPresentAnySolver( CurrentModel,'namespace') 1518 NameSpaceFound = .FALSE. 1519 1520 Mesh => CurrentModel % Meshes 1521 DO WHILE( ASSOCIATED( Mesh ) ) 1522 1523 CALL SetCurrentMesh( CurrentModel, Mesh ) 1524 1525 IF( InfoActive( 20 ) ) THEN 1526 PRINT *,'InitCond mesh:',TRIM(Mesh % Name), Mesh % MeshDim, Mesh % NumberOfNodes 1527 Var => Mesh % Variables 1528 DO WHILE( ASSOCIATED(Var) ) 1529 IF( ListCheckPresentAnyIC( CurrentModel, Var % Name ) ) THEN 1530 PRINT *,'InitCond pre range:',TRIM(Var % Name),MINVAL(Var % Values),MAXVAL( Var % Values) 1531 END IF 1532 Var => Var % Next 1533 END DO 1534 END IF 1535 1536 m = Mesh % MaxElementDofs 1537 n = Mesh % MaxElementNodes 1538 ALLOCATE( Indexes(m), Work(m) , Basis(m), Nodes % x(n), Nodes % y(n), Nodes % z(n) ) 1539 1540 ! First set the global variables and check whether there is anything left to do 1541 ThingsToDo = .FALSE. 1542 DO j=1,CurrentModel % NumberOfICs 1543 1544 IC => CurrentModel % ICs(j) % Values 1545 1546 Var => Mesh % Variables 1547 DO WHILE( ASSOCIATED(Var) ) 1548 1549 IF( .NOT. ASSOCIATED( Var % Values ) ) THEN 1550 Var => Var % Next 1551 CYCLE 1552 END IF 1553 1554 IF( SIZE( Var % Values ) == 0 ) THEN 1555 Var => Var % Next 1556 CYCLE 1557 END IF 1558 1559 Solver => Var % Solver 1560 IF ( .NOT. ASSOCIATED(Solver) ) Solver => CurrentModel % Solver 1561 1562 IF( AnyNameSpace ) THEN 1563 str = ListGetString( Solver % Values, 'Namespace', NamespaceFound ) 1564 IF (NamespaceFound) CALL ListPushNamespace(TRIM(str)) 1565 END IF 1566 1567 ! global variable 1568 IF( SIZE( Var % Values ) == Var % DOFs .OR. Var % Type == Variable_global ) THEN 1569 Val = ListGetCReal( IC, Var % Name, GotIt ) 1570 IF( GotIt ) THEN 1571 WRITE( Message,'(A,ES12.3)') 'Initializing global variable: > '& 1572 //TRIM(Var % Name)//' < to :',Val 1573 CALL Info('InitCond', Message,Level=8) 1574 Var % Values = Val 1575 END IF 1576 ELSE 1577 ThingsToDo = ThingsToDo .OR. & 1578 ListCheckPresent( IC, TRIM(Var % Name) ) 1579 ThingsToDo = ThingsToDo .OR. & 1580 ListCheckPresent( IC, TRIM(Var % Name)//' Velocity' ) 1581 ThingsToDo = ThingsToDo .OR. & 1582 ListCheckPresent( IC, TRIM(Var % Name)//' Acceleration' ) 1583 ThingsToDo = ThingsToDo .OR. & 1584 ListCheckPresent( IC, TRIM(Var % Name)//' {e}' ) 1585 END IF 1586 IF (NamespaceFound) CALL ListPopNamespace() 1587 Var => Var % Next 1588 END DO 1589 END DO 1590 1591 1592 ! And now do the ordinary fields 1593 !-------------------------------- 1594 IF( ThingsToDo ) THEN 1595 DO t=1, Mesh % NumberOfBulkElements+Mesh % NumberOfBoundaryElements 1596 1597 CurrentElement => Mesh % Elements(t) 1598 1599 i = CurrentElement % BodyId 1600 IF( i == 0 ) CYCLE 1601 1602 j = ListGetInteger(CurrentModel % Bodies(i) % Values, & 1603 'Initial Condition',GotIt, 1, CurrentModel % NumberOfICs ) 1604 IF ( .NOT. GotIt ) CYCLE 1605 1606 IC => CurrentModel % ICs(j) % Values 1607 CurrentModel % CurrentElement => CurrentElement 1608 n = GetElementNOFNodes() 1609 1610 Var => Mesh % Variables 1611 DO WHILE( ASSOCIATED(Var) ) 1612 1613 IF( .NOT. ASSOCIATED( Var % Values ) ) THEN 1614 Var => Var % Next 1615 CYCLE 1616 END IF 1617 1618 IF( SIZE( Var % Values ) == 0 ) THEN 1619 Var => Var % Next 1620 CYCLE 1621 END IF 1622 1623 IF( t == 1 ) THEN 1624 CALL Info('InitCond','Trying to initialize variable: '//TRIM(Var % Name),Level=20) 1625 END IF 1626 1627 Solver => Var % Solver 1628 IF ( .NOT. ASSOCIATED(Solver) ) THEN 1629 Solver => CurrentModel % Solver 1630 END IF 1631 CSolver => CurrentModel % Solver 1632 CurrentModel % Solver => Solver 1633 1634 IF( AnyNameSpace ) THEN 1635 str = ListGetString( Solver % Values, 'Namespace', NamespaceFound ) 1636 IF (NamespaceFound) CALL ListPushNamespace(TRIM(str)) 1637 END IF 1638 1639 ! global variables were already set 1640 IF( SIZE( Var % Values ) == Var % DOFs .OR. Var % Type == Variable_global ) THEN 1641 CONTINUE 1642 1643 ELSE IF( Var % TYPE == Variable_on_elements ) THEN 1644 IF( Var % DOFs > 1 ) THEN 1645 CALL Fatal('InitCond','Initialization only for scalar elements fields!') 1646 END IF 1647 1648 Work(1:n) = GetReal( IC, Var % Name, GotIt ) 1649 IF ( GotIt ) THEN 1650 k1 = CurrentElement % ElementIndex 1651 IF ( ASSOCIATED(Var % Perm) ) k1 = Var % Perm(k1) 1652 IF ( k1>0 ) Var % Values(k1) = SUM( Work(1:n) ) / n 1653 END IF 1654 1655 ELSE IF( Var % TYPE == Variable_on_gauss_points ) THEN 1656 ! We do this elsewhere in a more efficient manner 1657 CONTINUE 1658 1659 ELSE IF ( Var % DOFs == 1 ) THEN 1660 1661 Work(1:n) = ListGetReal( IC, Var % Name, n, CurrentElement % NodeIndexes, GotIt ) 1662 IF ( GotIt ) THEN 1663 ! Sometimes you may have both DG and bubbles, 1664 ! this way DG always has priority. 1665 IF( Var % TYPE == Variable_on_nodes_on_elements ) THEN 1666 IF(.NOT.ASSOCIATED(CurrentElement % DGIndexes)) GOTO 1 1667 Indexes(1:n) = CurrentElement % DgIndexes(1:n) 1668 ELSE 1669 DOFs = GetElementDOFs( Indexes, USolver=Var % Solver ) 1670 END IF 1671 1672 DO k=1,n 1673 k1 = Indexes(k) 1674 IF ( ASSOCIATED(Var % Perm) ) k1 = Var % Perm(k1) 1675 IF ( k1>0 ) Var % Values(k1) = Work(k) 1676 END DO 1677 16781 CONTINUE 1679 1680 END IF 1681 1682 IF ( Transient .AND. Solver % TimeOrder==2 ) THEN 1683 Work(1:n) = GetReal( IC, TRIM(Var % Name) // ' Velocity', GotIt ) 1684 IF ( GotIt ) THEN 1685 IF( Var % TYPE == Variable_on_nodes_on_elements ) THEN 1686 Indexes(1:n) = CurrentElement % DgIndexes(1:n) 1687 ELSE 1688 DOFs = GetElementDOFs( Indexes, USolver=Var % Solver ) 1689 END IF 1690 1691 DO k=1,n 1692 k1 = Indexes(k) 1693 IF ( ASSOCIATED(Var % Perm) ) k1 = Var % Perm(k1) 1694 IF ( k1>0 ) Var % PrevValues(k1,1) = Work(k) 1695 END DO 1696 END IF 1697 Work(1:n) = GetReal( IC, TRIM(Var % Name) // ' Acceleration', GotIt ) 1698 IF ( GotIt ) THEN 1699 IF( Var % TYPE == Variable_on_nodes_on_elements ) THEN 1700 Indexes(1:n) = CurrentElement % DgIndexes(1:n) 1701 ELSE 1702 DOFs = GetElementDOFs( Indexes, USolver=Var % Solver ) 1703 END IF 1704 1705 DO k=1,n 1706 k1 = Indexes(k) 1707 IF ( ASSOCIATED(Var % Perm) ) k1 = Var % Perm(k1) 1708 IF ( k1>0 ) Var % PrevValues(k1,2) = Work(k) 1709 END DO 1710 END IF 1711 END IF 1712 1713 IF(ASSOCIATED(Mesh % Edges)) THEN 1714 IF ( i<=Mesh % NumberOfBulkElements) THEN 1715 Gotit = ListCheckPresent( IC, TRIM(Var % Name)//' {e}' ) 1716 IF ( Gotit ) THEN 1717 DO k=1,CurrentElement % TYPE % NumberOfedges 1718 Edge => Mesh % Edges(CurrentElement % EdgeIndexes(k)) 1719 l = Var % Perm(CurrentElement % EdgeIndexes(k)+Mesh % NumberOfNodes) 1720 IF ( l>0 ) THEN 1721 CALL VectorElementEdgeDOFs( IC, & 1722 Edge, Edge % TYPE % NumberOfNodes, CurrentElement, n, & 1723 TRIM(Var % Name)//' {e}', Work ) 1724 Var % Values(l) = Work(1) 1725 END IF 1726 END DO 1727 END IF 1728 END IF 1729 END IF 1730 1731 ELSE 1732 CALL ListGetRealArray( IC, & 1733 Var % Name, WorkA, n, CurrentElement % NodeIndexes, gotIt ) 1734 1735 IF ( GotIt ) THEN 1736 DO k=1,n 1737 k1 = Indexes(k) 1738 IF ( ASSOCIATED(Var % Perm) ) k1 = Var % Perm(k1) 1739 IF(k1>0) THEN 1740 DO l=1,MIN(SIZE(WorkA,1),Var % DOFs) 1741 IF ( k1>0 ) Var % Values(Var % DOFs*(k1-1)+l) = WorkA(l,1,k) 1742 END DO 1743 END IF 1744 END DO 1745 END IF 1746 END IF 1747 IF(NamespaceFound) CALL ListPopNamespace() 1748 Var => Var % Next 1749 END DO 1750 CSolver => CurrentModel % Solver 1751 END DO 1752 1753 ! Here we do just the gauss point values for now. 1754 ! It would really make sense to do the ICs in this order since probably 1755 ! there are quite few IC variables to set but quite many elements. 1756 1757 Var => Mesh % Variables 1758 DO WHILE( ASSOCIATED(Var) ) 1759 1760 VarName = Var % Name 1761 Solver => Var % Solver 1762 IF ( .NOT. ASSOCIATED(Solver) ) Solver => CurrentModel % Solver 1763 1764 IF( AnyNameSpace ) THEN 1765 str = ListGetString( Solver % Values, 'Namespace', NamespaceFound ) 1766 IF (NamespaceFound) CALL ListPushNamespace(TRIM(str)) 1767 END IF 1768 1769 VarOrder = -1 1770 DO VarOrder = 0, 2 1771 IF( VarOrder == 0 ) THEN 1772 VarName = Var % Name 1773 ELSE IF( VarOrder == 1 ) THEN 1774 VarName = TRIM( Var % Name )//' Velocity' 1775 ELSE IF( VarOrder == 2 ) THEN 1776 VarName = TRIM( Var % Name )//' Acceleration' 1777 END IF 1778 1779 !CALL ListInitElementKeyword( LocalSol_h,'Initial Condition',VarName, & 1780 ! FoundSomewhere = Found ) 1781 Found = ListCheckPresentAnyIC( CurrentModel, VarName ) 1782 1783 IF( Found .AND. VarOrder > 0 ) THEN 1784 IF ( .NOT. ( Transient .AND. Solver % TimeOrder==2 ) ) THEN 1785 CALL Warn('InitCond','We can only set timederivative for transients') 1786 Found = .FALSE. 1787 END IF 1788 END IF 1789 1790 IF( Found ) THEN 1791 1792 CALL ListInitElementKeyword( LocalSol_h,'Initial Condition',VarName ) 1793 1794 IF( Var % TYPE /= Variable_on_gauss_points ) CYCLE 1795 1796 CALL Info('InitCond','Initializing gauss point field: '//TRIM(VarName),Level=7) 1797 IF( Var % DOFs > 1 ) THEN 1798 CALL Fatal('InitCond','Initialization only for scalar elemental fields!') 1799 END IF 1800 1801100 PrevBodyId = -1 1802 DO t=1, Mesh % NumberOfBulkElements+Mesh % NumberOfBoundaryElements 1803 1804 CurrentElement => Mesh % Elements(t) 1805 1806 i = CurrentElement % BodyId 1807 IF( i == 0 ) CYCLE 1808 1809 IF( i == PrevBodyId ) THEN 1810 FoundIC = PrevFoundIC 1811 ELSE 1812 j = ListGetInteger(CurrentModel % Bodies(i) % Values, & 1813 'Initial Condition',FoundIC, 1, CurrentModel % NumberOfICs ) 1814 IF ( FoundIC ) THEN 1815 IC => CurrentModel % ICs(j) % Values 1816 FoundIC = ListCheckPresent( IC, VarName ) 1817 END IF 1818 PrevFoundIC = FoundIC 1819 END IF 1820 1821 IF( .NOT. FoundIC ) CYCLE 1822 1823 CurrentModel % CurrentElement => CurrentElement 1824 n = GetElementNOFNodes() 1825 1826 k1 = Var % Perm( CurrentElement % ElementIndex ) 1827 k2 = Var % Perm( CurrentElement % ElementIndex + 1 ) 1828 1829 IF( k2- k1 > 0 ) THEN 1830 1831 IP = GaussPointsAdapt( CurrentElement, Solver ) 1832 1833 IF( k2 - k1 /= Ip % n ) THEN 1834 CALL Info('InitCond','Number of Gauss points has changed, redoing permutations!',Level=8) 1835 CALL UpdateIpPerm( Solver, Var % Perm ) 1836 nsize = MAXVAL( Var % Perm ) 1837 1838 CALL Info('InitCond','Total number of new IP dofs: '//TRIM(I2S(nsize))) 1839 1840 IF( SIZE( Var % Values ) /= Var % Dofs * nsize ) THEN 1841 DEALLOCATE( Var % Values ) 1842 ALLOCATE( Var % Values( nsize * Var % Dofs ) ) 1843 END IF 1844 Var % Values = 0.0_dp 1845 GOTO 100 1846 END IF 1847 1848 Nodes % x(1:n) = Mesh % Nodes % x(CurrentElement % NodeIndexes) 1849 Nodes % y(1:n) = Mesh % Nodes % y(CurrentElement % NodeIndexes) 1850 Nodes % z(1:n) = Mesh % Nodes % z(CurrentElement % NodeIndexes) 1851 1852 DO k=1,IP % n 1853 stat = ElementInfo( CurrentElement, Nodes, IP % U(k), IP % V(k), & 1854 IP % W(k), detJ, Basis ) 1855 1856 val = ListGetElementReal( LocalSol_h,Basis,CurrentElement,Found,GaussPoint=k) 1857 1858 IF( VarOrder == 0 ) THEN 1859 Var % Values(k1+k) = val 1860 ELSE 1861 Var % PrevValues(k1+k,VarOrder) = val 1862 END IF 1863 1864 END DO 1865 1866 END IF 1867 END DO 1868 END IF 1869 END DO 1870 IF(NamespaceFound) CALL ListPopNamespace() 1871 Var => Var % Next 1872 END DO 1873 END IF 1874 1875 DEALLOCATE( Indexes, Work, Basis, Nodes % x, Nodes % y, Nodes % z) 1876 1877 IF( InfoActive( 20 ) ) THEN 1878 Var => Mesh % Variables 1879 DO WHILE( ASSOCIATED(Var) ) 1880 IF( ListCheckPresentAnyIC( CurrentModel, Var % Name ) ) THEN 1881 PRINT *,'InitCond post range:',TRIM(Var % Name),& 1882 MINVAL(Var % Values),MAXVAL( Var % Values),SUM(Var % Values)/SIZE(Var % Values) 1883 END IF 1884 Var => Var % Next 1885 END DO 1886 END IF 1887 1888 Mesh => Mesh % Next 1889 END DO 1890 1891 1892 1893!------------------------------------------------------------------------------ 1894 END SUBROUTINE InitCond 1895!------------------------------------------------------------------------------ 1896 1897!------------------------------------------------------------------------------ 1898!> Check if we are restarting are if yes, read in field values. 1899!------------------------------------------------------------------------------ 1900 SUBROUTINE Restart() 1901!------------------------------------------------------------------------------ 1902 USE DefUtils 1903 LOGICAL :: Gotit, DoIt 1904 INTEGER :: i, j, k, l 1905 REAL(KIND=dp) :: StartTime 1906 TYPE(Mesh_t), POINTER :: Mesh, pMesh 1907 TYPE(ValueList_t), POINTER :: RestartList 1908 LOGICAL, ALLOCATABLE :: MeshDone(:) 1909 INTEGER, POINTER :: MeshesToRestart(:) 1910 LOGICAL :: CheckMesh, DoMesh 1911!------------------------------------------------------------------------------ 1912 1913 1914 ! Count the number of meshes first so that we can identify them 1915 j = 0 1916 pMesh => CurrentModel % Meshes 1917 DO WHILE( ASSOCIATED(pMesh) ) 1918 j = j + 1 1919 pMesh => pMesh % Next 1920 END DO 1921 ALLOCATE( MeshDone( j ) ) 1922 MeshDone = .FALSE. 1923 1924 1925 ! Do Solver-mesh specific restart only 1926 !----------------------------------------------------------------- 1927 IF ( ListCheckPresentAnySolver( CurrentModel,'Restart File') ) THEN 1928 DO i=1, CurrentModel % NumberOfSolvers 1929 RestartList => CurrentModel % Solvers(i) % Values 1930 1931 RestartFile = ListGetString( RestartList, 'Restart File', GotIt ) 1932 IF ( GotIt ) THEN 1933 1934 Mesh => CurrentModel % Solvers(i) % Mesh 1935 IF( .NOT. ASSOCIATED(Mesh) ) THEN 1936 CALL Warn('Restart','Solver has no mesh associated!') 1937 CYCLE 1938 END IF 1939 1940 DoIt = .TRUE. 1941 pMesh => CurrentModel % Meshes 1942 j = 0 1943 DO WHILE( ASSOCIATED(pMesh) ) 1944 j = j + 1 1945 pMesh => pMesh % Next 1946 IF( ASSOCIATED( Mesh, pMesh ) ) THEN 1947 IF( MeshDone(j) ) THEN 1948 DoIt = .FALSE. 1949 ELSE 1950 MeshDone(j) = .TRUE. 1951 END IF 1952 END IF 1953 END DO 1954 1955 ! Variables for this mesh already done! 1956 IF(.NOT. DoIt ) CYCLE 1957 1958 CALL Info('Restart','Perfoming solver specific Restart for: '//TRIM(Mesh % Name),Level=6) 1959 IF ( LEN_TRIM(Mesh % Name) > 0 ) THEN 1960 OutputName = TRIM(Mesh % Name) // '/' // TRIM(RestartFile) 1961 ELSE 1962 OutputName = TRIM(RestartFile) 1963 END IF 1964 1965 IF ( ParEnv % PEs > 1 .AND. .NOT. Mesh % SingleMesh ) & 1966 OutputName = TRIM(OutputName) // '.' // TRIM(i2s(ParEnv % MyPe)) 1967 CALL SetCurrentMesh( CurrentModel, Mesh ) 1968 1969 k = ListGetInteger( RestartList,'Restart Position',GotIt, minv=0 ) 1970 CALL LoadRestartFile( OutputName, k, Mesh, SolverId = i ) 1971 1972 StartTime = ListGetConstReal( RestartList ,'Restart Time',GotIt) 1973 IF( GotIt ) THEN 1974 Var => VariableGet( Mesh % Variables, 'Time' ) 1975 IF ( ASSOCIATED( Var ) ) Var % Values(1) = StartTime 1976 END IF 1977 END IF 1978 1979 END DO 1980 END IF 1981 1982 ! Do the standard global restart 1983 !----------------------------------------------------------------- 1984 RestartList => CurrentModel % Simulation 1985 1986 ! We may supress restart from certain meshes. 1987 ! This was initially only related to calving, but no need to limit to that. 1988 l = 0 1989 MeshesToRestart => ListGetIntegerArray(RestartList,& 1990 'Meshes To Restart', CheckMesh ) 1991 1992 RestartFile = ListGetString( RestartList, 'Restart File', GotIt ) 1993 IF ( GotIt ) THEN 1994 k = ListGetInteger( RestartList,'Restart Position',GotIt, minv=0 ) 1995 Mesh => CurrentModel % Meshes 1996 1997 j = 0 1998 DO WHILE( ASSOCIATED(Mesh) ) 1999 j = j + 1 2000 2001 ! Make sure that if a mesh has already been restarted 2002 ! it is not being done again. 2003 IF( MeshDone(j) ) THEN 2004 CALL Info('Restart','Already done mesh: '//TRIM(Mesh % Name)) 2005 Mesh => Mesh % Next 2006 CYCLE 2007 END IF 2008 MeshDone(j) = .TRUE. 2009 2010 CALL Info('Restart','Perfoming global Restart for: '//TRIM(Mesh % Name),Level=6) 2011 2012 IF ( LEN_TRIM(Mesh % Name) > 0 ) THEN 2013 OutputName = TRIM(Mesh % Name) // '/' // TRIM(RestartFile) 2014 ELSE 2015 OutputName = TRIM(RestartFile) 2016 END IF 2017 IF ( ParEnv % PEs > 1 .AND. .NOT. Mesh % SingleMesh ) & 2018 OutputName = TRIM(OutputName) // '.' // TRIM(i2s(ParEnv % MyPe)) 2019 2020 l = l+1 2021 2022 DoMesh = .TRUE. 2023 IF(CheckMesh) THEN 2024 DoMesh = (ANY(MeshesToRestart == l)) 2025 END IF 2026 2027 IF( DoMesh ) THEN 2028 CALL SetCurrentMesh( CurrentModel, Mesh ) 2029 CALL LoadRestartFile( OutputName, k, Mesh ) 2030 2031 StartTime = ListGetConstReal( RestartList ,'Restart Time',GotIt) 2032 IF( GotIt ) THEN 2033 Var => VariableGet( Mesh % Variables, 'Time' ) 2034 IF ( ASSOCIATED( Var ) ) Var % Values(1) = StartTime 2035 END IF 2036 END IF 2037 Mesh => Mesh % Next 2038 END DO 2039 ELSE 2040 StartTime = ListGetConstReal( CurrentModel % Simulation ,'Start Time',GotIt) 2041 Mesh => CurrentModel % Meshes 2042 2043 j = 0 2044 IF( GotIt ) THEN 2045 DO WHILE( ASSOCIATED(Mesh) ) 2046 j = j + 1 2047 2048 ! Make sure that if a mesh has already been restarted 2049 ! it is not being done again. 2050 IF( MeshDone(j) ) THEN 2051 CALL Info('Restart','Already done mesh: '//TRIM(Mesh % Name)) 2052 Mesh => Mesh % Next 2053 CYCLE 2054 END IF 2055 MeshDone(j) = .TRUE. 2056 Var => VariableGet( Mesh % Variables, 'Time' ) 2057 IF ( ASSOCIATED( Var ) ) Var % Values(1) = StartTime 2058 WRITE(Message,*) 'Found "Start Time" (without restart) and setting time-variable to t=',StartTime 2059 CALL INFO("Restart",Message,Level=1) 2060 END DO 2061 END IF 2062 END IF 2063 2064 2065!------------------------------------------------------------------------------ 2066 END SUBROUTINE Restart 2067!------------------------------------------------------------------------------ 2068 2069 2070!------------------------------------------------------------------------------ 2071!> Execute the individual solvers in defined sequence. 2072!------------------------------------------------------------------------------ 2073 SUBROUTINE ExecSimulation(TimeIntervals, CoupledMinIter, & 2074 CoupledMaxIter, OutputIntervals, Transient, Scanning) 2075 IMPLICIT NONE 2076 INTEGER :: TimeIntervals,CoupledMinIter, CoupledMaxIter,OutputIntervals(:) 2077 LOGICAL :: Transient,Scanning 2078!------------------------------------------------------------------------------ 2079 INTEGER :: interval, timestep, i, j, k, n 2080 REAL(KIND=dp) :: dt, ddt, dtfunc, timeleft 2081 INTEGER :: cum_timestep 2082 INTEGER, SAVE :: stepcount=0, RealTimestep 2083 LOGICAL :: ExecThis,SteadyStateReached=.FALSE.,PredCorrControl, & 2084 DivergenceControl, HaveDivergence 2085 REAL(KIND=dp) :: CumTime, MaxErr, AdaptiveLimit, & 2086 AdaptiveMinTimestep, AdaptiveMaxTimestep, timePeriod 2087 INTEGER :: SmallestCount, AdaptiveKeepSmallest, StepControl=-1, nSolvers 2088 LOGICAL :: AdaptiveTime = .TRUE., AdaptiveRough, AdaptiveSmart, Found 2089 INTEGER :: AllocStat 2090 REAL(KIND=dp) :: AdaptiveIncrease, AdaptiveDecrease 2091 TYPE(Solver_t), POINTER :: Solver 2092 TYPE AdaptiveVariables_t 2093 TYPE(Variable_t) :: Var 2094 REAL(KIND=dp) :: Norm 2095 END TYPE AdaptiveVariables_t 2096 TYPE(AdaptiveVariables_t), ALLOCATABLE, SAVE :: AdaptVars(:) 2097 REAL(KIND=dp) :: newtime, prevtime=0, maxtime, exitcond 2098 INTEGER, SAVE :: PrevMeshI = 0 2099 2100 !$OMP PARALLEL 2101 IF(.NOT.GaussPointsInitialized()) CALL GaussPointsInit() 2102 !$OMP END PARALLEL 2103 2104 nSolvers = CurrentModel % NumberOfSolvers 2105 DO i=1,nSolvers 2106 Solver => CurrentModel % Solvers(i) 2107 IF ( Solver % PROCEDURE==0 ) CYCLE 2108 IF ( Solver % SolverExecWhen == SOLVER_EXEC_AHEAD_ALL ) THEN 2109 ! solver to be called prior to time looping can never be transient 2110 dt = 1.0_dp 2111 CALL SolverActivate( CurrentModel,Solver,dt,.FALSE. ) 2112 END IF 2113 END DO 2114 2115 IF( ListGetLogical( CurrentModel % Simulation,'Calculate Mesh Pieces',Found ) ) THEN 2116 CALL CalculateMeshPieces( CurrentModel % Mesh ) 2117 END IF 2118 2119 ! Predictor-Corrector time stepping control 2120 PredCorrControl = ListGetLogical( CurrentModel % Simulation, & 2121 'Predictor-Corrector Control', gotIt) 2122 2123 ! Divergence control 2124 DivergenceControl = ListGetLogical( CurrentModel % Simulation, & 2125 'Convergence Control', gotIt) 2126 2127 AdaptiveTime = ListGetLogical( CurrentModel % Simulation, & 2128 'Adaptive Timestepping', GotIt ) 2129 2130 DO interval = 1, TimeIntervals 2131 stepcount = stepcount + Timesteps(interval) 2132 END DO 2133 2134 dt = 1.0_dp 2135 cum_Timestep = 0 2136 ddt = -1.0_dp 2137 DO interval = 1,TimeIntervals 2138 2139!------------------------------------------------------------------------------ 2140! go through number of timesteps within an interval 2141!------------------------------------------------------------------------------ 2142 timePeriod = ListGetCReal(CurrentModel % Simulation, 'Time Period',gotIt) 2143 IF(.NOT.GotIt) timePeriod = HUGE(timePeriod) 2144 2145 IF(GetNameSpaceCheck()) THEN 2146 IF(Scanning) THEN 2147 CALL ListPushNamespace('scan:') 2148 ELSE IF(Transient) THEN 2149 CALL ListPushNamespace('time:') 2150 ELSE 2151 CALL ListPushNamespace('steady:') 2152 END IF 2153 END IF 2154 2155 RealTimestep = 1 2156 DO timestep = 1,Timesteps(interval) 2157 2158 cum_Timestep = cum_Timestep + 1 2159 sStep(1) = cum_Timestep 2160 2161 IF ( GetNamespaceCheck() ) THEN 2162 IF( Scanning ) THEN 2163 CALL ListPushNamespace('scan '//TRIM(i2s(cum_Timestep))//':') 2164 ELSE IF ( Transient ) THEN 2165 CALL ListPushNamespace('time '//TRIM(i2s(cum_Timestep))//':') 2166 ELSE 2167 CALL ListPushNamespace('steady '//TRIM(i2s(cum_Timestep))//':') 2168 END IF 2169 END IF 2170 2171 IF ( Transient .OR. Scanning ) THEN 2172 dtfunc = ListGetCReal( CurrentModel % Simulation,'Timestep Function',GotIt ) 2173 IF(GotIt) THEN 2174 CALL Warn('ExecSimulation','Obsolete keyword > Timestep Function < , use > Timestep Size < instead') 2175 ELSE 2176 dtfunc = ListGetCReal( CurrentModel % Simulation,'Timestep Size', gotIt) 2177 END IF 2178 IF( GotIt ) THEN 2179 BLOCK 2180 TYPE(Variable_t), POINTER :: tVar 2181 INTEGER :: dtiter 2182 REAL(KIND=dp) :: t0 2183 dtiter = ListGetInteger( CurrentModel % Simulation,'Timestep Size Iterations',Found) 2184 2185 ! If timestep size depends on time i.e. dt=dt(t) we may sometimes want to iterate 2186 ! so that the timestep is consistent with the new time t+dt i.e. dt=dt(t+dt). 2187 IF( dtiter > 0 ) THEN 2188 tVar => VariableGet( CurrentModel % Mesh % Variables, 'time' ) 2189 IF(ASSOCIATED(tVar)) THEN 2190 t0 = tVar % Values(1) 2191 ! Iterate for consistent time + timestep size. 2192 DO i=1,dtiter 2193 tVar % Values(1) = t0 + dtfunc 2194 dtfunc = ListGetCReal( CurrentModel % Simulation,'Timestep Size' ) 2195 END DO 2196 ! Revert back to original time, this will be added later on again. 2197 tVar % Values(1) = t0 2198 END IF 2199 END IF 2200 END BLOCK 2201 END IF 2202 IF(GotIt) THEN 2203 dt = dtfunc 2204 IF(dt < EPSILON(dt) ) THEN 2205 WRITE(Message,'(A,ES12.3)') 'Timestep smaller than epsilon: ',dt 2206 CALL Fatal('ExecSimulation', Message) 2207 END IF 2208 ELSE 2209 dt = TimestepSizes(interval,1) 2210 END IF 2211 END IF 2212 2213!------------------------------------------------------------------------------ 2214 ! Predictor-Corrector time stepping control 2215 IF ( PredCorrControl ) THEN 2216 CALL PredictorCorrectorControl( CurrentModel, dt, timestep ) 2217 END IF 2218 2219!------------------------------------------------------------------------------ 2220 2221 IF( AdaptiveTime ) THEN 2222 AdaptiveLimit = ListGetConstReal( CurrentModel % Simulation, & 2223 'Adaptive Time Error', GotIt ) 2224 IF ( .NOT. GotIt ) THEN 2225 CALL Fatal('ElmerSolver','Adaptive Time Error must be given for ' // & 2226 'adaptive stepping scheme.') 2227 END IF 2228 AdaptiveKeepSmallest = ListGetInteger( CurrentModel % Simulation, & 2229 'Adaptive Keep Smallest', GotIt, minv=0 ) 2230 END IF 2231 2232 IF( AdaptiveTime .OR. DivergenceControl ) THEN 2233 AdaptiveMaxTimestep = ListGetConstReal( CurrentModel % Simulation, & 2234 'Adaptive Max Timestep', GotIt ) 2235 IF ( GotIt ) THEN 2236 AdaptiveMaxTimestep = MIN(AdaptiveMaxTimeStep, dt) 2237 ELSE 2238 AdaptiveMaxTimestep = dt 2239 END IF 2240 2241 AdaptiveMinTimestep = ListGetConstReal( CurrentModel % Simulation, & 2242 'Adaptive Min Timestep', GotIt ) 2243 IF(.NOT. GotIt) AdaptiveMinTimestep = 1.0e-8 * AdaptiveMaxTimestep 2244 2245 AdaptiveIncrease = ListGetConstReal( CurrentModel % Simulation, & 2246 'Adaptive Increase Coefficient', GotIt ) 2247 IF(.NOT. GotIt) AdaptiveIncrease = 2.0_dp 2248 2249 AdaptiveDecrease = ListGetConstReal( CurrentModel % Simulation, & 2250 'Adaptive Decrease Coefficient', GotIt ) 2251 IF(.NOT. GotIt) AdaptiveDecrease = 0.5_dp 2252 2253 AdaptiveRough = ListGetLogical( CurrentModel % Simulation, & 2254 'Adaptive Rough Timestep', GotIt ) 2255 2256 AdaptiveSmart = ListGetLogical( CurrentModel % Simulation, & 2257 'Adaptive Smart Timestep', GotIt ) 2258 END IF 2259 2260 2261!------------------------------------------------------------------------------ 2262 sTime(1) = sTime(1) + dt 2263 sPeriodic(1) = sTime(1) 2264 DO WHILE(sPeriodic(1) > timePeriod) 2265 sPeriodic(1) = sPeriodic(1) - timePeriod 2266 END DO 2267 2268 ! Move the old timesteps one step down the ladder 2269 IF(timestep > 1 .OR. interval > 1) THEN 2270 DO i = SIZE(sPrevSizes,2),2,-1 2271 sPrevSizes(1,i) = sPrevSizes(1,i-1) 2272 END DO 2273 sPrevSizes(1,1) = sSize(1) 2274 END IF 2275 sSize(1) = dt 2276 2277 sInterval(1) = interval 2278 IF (.NOT. Transient ) steadyIt(1) = steadyIt(1) + 1 2279 2280!------------------------------------------------------------------------------ 2281 2282 BLOCK 2283 TYPE(Mesh_t), POINTER :: Mesh 2284 REAL(KIND=dp) :: MeshR 2285 CHARACTER(LEN=MAX_NAME_LEN) :: MeshStr 2286 2287 IF( ListCheckPresent( GetSimulation(), 'Mesh Name Index') ) THEN 2288 IF( Transient ) THEN 2289 CALL Fatal('ExecSimulation','Mesh swapping not supported in transient!') 2290 END IF 2291 2292 ! we cannot have mesh depend on "time" or "timestep" if they are not available as 2293 ! variables. 2294 Mesh => CurrentModel % Meshes 2295 CALL VariableAdd( Mesh % Variables, Mesh, Name='Time',DOFs=1, Values=sTime ) 2296 CALL VariableAdd( Mesh % Variables, Mesh, Name='Timestep', DOFs=1, Values=sStep ) 2297 2298 MeshR = GetCREal( GetSimulation(), 'Mesh Name Index',gotIt ) 2299 i = NINT( MeshR ) 2300 2301 IF( i > 0 .AND. i /= PrevMeshI ) THEN 2302 MeshStr = ListGetString( GetSimulation(),'Mesh Name '//TRIM(I2S(i)),GotIt) 2303 IF( GotIt ) THEN 2304 CALL Info('ExecSimulation','Swapping mesh to: '//TRIM(MeshStr),Level=5) 2305 ELSE 2306 CALL Fatal('ExecSimulation','Could not find >Mesh Name '//TRIM(I2S(i))//'<') 2307 END IF 2308 CALL SwapMesh( CurrentModel, Mesh, MeshStr ) 2309 PrevMeshI = i 2310 END IF 2311 END IF 2312 END BLOCK 2313 2314 2315!------------------------------------------------------------------------------ 2316 IF ( ParEnv % MyPE == 0 ) THEN 2317 CALL Info( 'MAIN', ' ', Level=3 ) 2318 CALL Info( 'MAIN', '-------------------------------------', Level=3 ) 2319 2320 IF ( Transient .OR. Scanning ) THEN 2321 WRITE( Message, * ) 'Time: ',TRIM(i2s(cum_Timestep)),'/', & 2322 TRIM(i2s(stepcount)), sTime(1) 2323 CALL Info( 'MAIN', Message, Level=3 ) 2324 2325 newtime= RealTime() 2326 2327 IF( cum_Timestep > 1 ) THEN 2328 maxtime = ListGetConstReal( CurrentModel % Simulation,'Real Time Max',GotIt) 2329 IF( GotIt ) THEN 2330 WRITE( Message,'(A,F8.3)') 'Fraction of real time left: ',& 2331 1.0_dp-RealTime() / maxtime 2332 END IF 2333 2334 ! Compute estimated time left in seconds 2335 timeleft = (stepcount-(cum_Timestep-1))*(newtime-prevtime) 2336 2337 ! No sense to show too short estimated times 2338 IF( timeleft > 1 ) THEN 2339 IF (timeleft >= 24 * 3600) THEN ! >24 hours 2340 WRITE( Message,'(A)') 'Estimated time left: '//I2S(NINT(timeleft/3600))//' hours' 2341 ELSE IF (timeleft >= 3600) THEN ! 1 to 20 hours 2342 WRITE( Message,'(A,F5.1,A)') 'Estimated time left:',timeleft/3600,' hours' 2343 ELSE IF(timeleft >= 60) THEN ! 1 to 60 minutes 2344 WRITE( Message,'(A,F5.1,A)') 'Estimated time left:',timeleft/60,' minutes' 2345 ELSE ! 1 to 60 seconds 2346 WRITE( Message,'(A,F5.1,A)') 'Estimated time left:',timeleft,' seconds' 2347 END IF 2348 CALL Info( 'MAIN', Message, Level=3 ) 2349 END IF 2350 2351 END IF 2352 prevtime = newtime 2353 ELSE 2354 WRITE( Message, * ) 'Steady state iteration: ',cum_Timestep 2355 CALL Info( 'MAIN', Message, Level=3 ) 2356 END IF 2357 2358 CALL Info( 'MAIN', '-------------------------------------', Level=3 ) 2359 CALL Info( 'MAIN', ' ', Level=3 ) 2360 END IF 2361 2362!------------------------------------------------------------------------------ 2363! Solve any and all governing equations in the system 2364!------------------------------------------------------------------------------ 2365 2366 IF ( Transient .AND. AdaptiveTime ) THEN 2367 2368 IF(.NOT. ALLOCATED( AdaptVars ) ) THEN 2369 ALLOCATE( AdaptVars( nSolvers ), STAT = AllocStat ) 2370 IF( AllocStat /= 0 ) CALL Fatal('ExecSimulation','Allocation error for AdaptVars') 2371 2372 DO i=1,nSolvers 2373 Solver => CurrentModel % Solvers(i) 2374 2375 NULLIFY( AdaptVars(i) % Var % Values ) 2376 NULLIFY( AdaptVars(i) % Var % PrevValues ) 2377 2378 IF( .NOT. ASSOCIATED( Solver % Variable ) ) CYCLE 2379 IF( .NOT. ASSOCIATED( Solver % Variable % Values ) ) CYCLE 2380 CALL Info('ExecSimulation','Allocating adaptive work space for: '//TRIM(I2S(i)),Level=12) 2381 j = SIZE( Solver % Variable % Values ) 2382 ALLOCATE( AdaptVars(i) % Var % Values( j ), STAT=AllocStat ) 2383 IF( AllocStat /= 0 ) CALL Fatal('ExecSimulation','Allocation error AdaptVars Values') 2384 2385 IF( ASSOCIATED( Solver % Variable % PrevValues ) ) THEN 2386 k = SIZE( Solver % Variable % PrevValues, 2 ) 2387 ALLOCATE( AdaptVars(i) % Var % PrevValues( j, k ), STAT=AllocStat) 2388 IF( AllocStat /= 0 ) CALL Fatal('ExecSimulation','Allocation error for AdaptVars PrevValues') 2389 END IF 2390 END DO 2391 END IF 2392 2393 CumTime = 0.0d0 2394 IF ( ddt < 0.0_dp .OR. ddt > AdaptiveMaxTimestep ) ddt = AdaptiveMaxTimestep 2395 2396 s = sTime(1) - dt 2397 SmallestCount = 0 2398 DO WHILE( CumTime < dt-1.0d-12 ) 2399 IF( .NOT. AdaptiveRough ) THEN 2400 ddt = MIN( dt - CumTime, ddt ) 2401 IF( AdaptiveSmart ) THEN 2402 ! If the next timestep will not get us home but the next one would 2403 ! then split the timestep equally into two parts. 2404 IF( dt - CumTime - ddt > 1.0d-12 ) THEN 2405 CALL Info('ExecSimulation','Splitted timestep into two equal parts',Level=12) 2406 ddt = MIN( ddt, ( dt - CumTime ) / 2.0_dp ) 2407 END IF 2408 END IF 2409 2410 END IF 2411 2412 ! Store the initial values before the start of the step 2413 DO i=1,nSolvers 2414 Solver => CurrentModel % Solvers(i) 2415 IF ( .NOT. ASSOCIATED( Solver % Variable ) ) CYCLE 2416 IF ( .NOT. ASSOCIATED( Solver % Variable % Values ) ) CYCLE 2417 AdaptVars(i) % Var % Values = Solver % Variable % Values 2418 AdaptVars(i) % Var % Norm = Solver % Variable % Norm 2419 IF ( ASSOCIATED( Solver % Variable % PrevValues ) ) THEN 2420 AdaptVars(i) % Var % PrevValues = Solver % Variable % PrevValues 2421 END IF 2422 END DO 2423 2424 sTime(1) = s + CumTime + ddt 2425 sSize(1) = ddt 2426 2427 ! Solve with full timestep 2428 CALL SolveEquations( CurrentModel, ddt, Transient, & 2429 CoupledMinIter, CoupledMaxIter, SteadyStateReached, RealTimestep, & 2430 BeforeTime = .TRUE., AtTime = .TRUE., AfterTime = .FALSE.) 2431 2432 ! External adaptive error given 2433 MaxErr = ListGetConstReal( CurrentModel % Simulation, & 2434 'Adaptive Error Measure', GotIt ) 2435 2436 DO i=1,nSolvers 2437 Solver => CurrentModel % Solvers(i) 2438 IF ( .NOT. ASSOCIATED( Solver % Variable ) ) CYCLE 2439 IF ( .NOT. ASSOCIATED( Solver % Variable % Values ) ) CYCLE 2440 Solver % Variable % Values = AdaptVars(i) % Var % Values 2441 AdaptVars(i) % Norm = Solver % Variable % Norm 2442 IF ( ASSOCIATED( Solver % Variable % PrevValues ) ) THEN 2443 Solver % Variable % PrevValues = AdaptVars(i) % Var % PrevValues 2444 END IF 2445 END DO 2446 2447 ! Test the error for half the timestep 2448 sStep(1) = ddt / 2 2449 sTime(1) = s + CumTime + ddt/2 2450 CALL SolveEquations( CurrentModel, ddt/2, Transient, & 2451 CoupledMinIter, CoupledMaxIter, SteadyStateReached, RealTimestep, & 2452 BeforeTime = .TRUE., AtTime = .TRUE., AfterTime = .FALSE.) 2453 2454 sTime(1) = s + CumTime + ddt 2455 CALL SolveEquations( CurrentModel, ddt/2, Transient, & 2456 CoupledMinIter, CoupledMaxIter, SteadyStateReached, RealTimestep, & 2457 BeforeTime = .TRUE., AtTime = .TRUE., AfterTime = .FALSE.) 2458 2459 MaxErr = ABS( MaxErr - ListGetConstReal( CurrentModel % Simulation, & 2460 'Adaptive Error Measure', GotIt ) ) 2461 2462 ! If not measure given then use the maximum change in norm as the measure 2463 IF ( .NOT. GotIt ) THEN 2464 MaxErr = 0.0d0 2465 DO i=1,nSolvers 2466 Solver => CurrentModel % Solvers(i) 2467 IF ( .NOT. ASSOCIATED( Solver % Variable ) ) CYCLE 2468 IF ( .NOT. ASSOCIATED( Solver % Variable % Values ) ) CYCLE 2469 IF ( AdaptVars(i) % norm /= Solver % Variable % Norm ) THEN 2470 Maxerr = MAX(Maxerr,ABS(AdaptVars(i) % norm - Solver % Variable % Norm)/& 2471 AdaptVars(i) % norm ) 2472 END IF 2473 END DO 2474 END IF 2475 2476 ! We have a success no need to redo this step 2477 IF ( MaxErr < AdaptiveLimit .OR. ddt <= AdaptiveMinTimestep ) THEN 2478 CumTime = CumTime + ddt 2479 RealTimestep = RealTimestep+1 2480 IF ( SmallestCount >= AdaptiveKeepSmallest .OR. StepControl > 0 ) THEN 2481 ddt = MIN( AdaptiveIncrease * ddt, AdaptiveMaxTimeStep ) 2482 StepControl = 1 2483 SmallestCount = 0 2484 ELSE 2485 StepControl = 0 2486 SmallestCount = SmallestCount + 1 2487 END IF 2488 2489 ! Finally solve only the postprocessing solver after the timestep has been accepted 2490 CALL SolveEquations( CurrentModel, ddt/2, Transient, & 2491 CoupledMinIter, CoupledMaxIter, SteadyStateReached, RealTimestep, & 2492 BeforeTime = .FALSE., AtTime = .FALSE., AfterTime = .TRUE.) 2493 ELSE 2494 DO i=1,nSolvers 2495 Solver => CurrentModel % Solvers(i) 2496 IF ( .NOT. ASSOCIATED( Solver % Variable ) ) CYCLE 2497 IF ( .NOT. ASSOCIATED( Solver % Variable % Values ) ) CYCLE 2498 Solver % Variable % Norm = AdaptVars(i) % Var % Norm 2499 Solver % Variable % Values = AdaptVars(i) % Var % Values 2500 IF ( ASSOCIATED( Solver % Variable % PrevValues ) ) THEN 2501 Solver % Variable % PrevValues = AdaptVars(i) % Var % PrevValues 2502 END IF 2503 END DO 2504 ddt = AdaptiveDecrease * ddt 2505 StepControl = -1 2506 END IF 2507 2508 WRITE(*,'(a,3e20.12)') 'Adaptive(cum,ddt,err): ', cumtime, ddt, maxerr 2509 END DO 2510 sSize(1) = dt 2511 sTime(1) = s + dt 2512 2513 ELSE IF( DivergenceControl ) THEN 2514 ! This is still tentative 2515 CALL Info('ExecSimulation','Solving equations with divergence control',Level=6) 2516 2517 CumTime = 0.0d0 2518 ddt = AdaptiveIncrease * ddt 2519 IF ( ddt < 0.0_dp .OR. ddt > AdaptiveMaxTimestep ) ddt = AdaptiveMaxTimestep 2520 s = sTime(1) - dt 2521 !StepControl = 0 2522 2523 DO WHILE( CumTime < dt-1.0d-12 ) 2524 2525 IF( .NOT. AdaptiveRough ) THEN 2526 ddt = MIN( dt - CumTime, ddt ) 2527 IF( AdaptiveSmart ) THEN 2528 IF( dt - CumTime - ddt > 1.0d-12 ) THEN 2529 CALL Info('ExecSimulation','Splitted timestep into two equal parts',Level=12) 2530 ddt = MIN( ddt, ( dt - CumTime ) / 2.0_dp ) 2531 END IF 2532 END IF 2533 END IF 2534 2535 sTime(1) = s + CumTime + ddt 2536 sSize(1) = ddt 2537 2538 CALL SolveEquations( CurrentModel, ddt, Transient, & 2539 CoupledMinIter, CoupledMaxIter, SteadyStateReached, RealTimestep, & 2540 BeforeTime = .TRUE., AtTime = .TRUE., AfterTime = .FALSE.) 2541 2542 HaveDivergence = .FALSE. 2543 DO i=1,nSolvers 2544 Solver => CurrentModel % Solvers(i) 2545 IF( ASSOCIATED( Solver % Variable ) ) THEN 2546 IF( Solver % Variable % NonlinConverged > 1 ) THEN 2547 CALL Info('ExecSimulation','Solver '//TRIM(I2S(i))//' has diverged',Level=8) 2548 HaveDivergence = .TRUE. 2549 EXIT 2550 END IF 2551 END IF 2552 END DO 2553 2554 IF( .NOT. HaveDivergence ) THEN 2555 CALL Info('ExecSimulation','No solver has diverged',Level=8) 2556 2557 ! Finally solve only the postprocessing solver after the timestep has been accepted 2558 CALL SolveEquations( CurrentModel, ddt, Transient, & 2559 CoupledMinIter, CoupledMaxIter, SteadyStateReached, RealTimestep, & 2560 BeforeTime = .FALSE., AtTime = .FALSE., AfterTime = .TRUE.) 2561 2562 ! If step control was active for this interval then we can 2563 ! again start to increase timestep, otherwise not 2564 2565 CumTime = CumTime + ddt 2566 RealTimestep = RealTimestep+1 2567 2568 !IF ( StepControl > 0 ) THEN 2569 ddt = AdaptiveIncrease * ddt 2570 IF( ddt > AdaptiveMaxTimestep ) THEN 2571 ddt = AdaptiveMaxTimestep 2572 StepControl = 0 2573 END IF 2574 !END IF 2575 ELSE 2576 IF( ddt < AdaptiveMinTimestep * (1+1.0e-8) ) THEN 2577 CALL Fatal('ExecSimulation','Could not find stable timestep above given minimum') 2578 END IF 2579 2580 CALL Info('ExecSimulation','Reducing timestep due to divergence problems!',Level=6) 2581 2582 ddt = MAX( AdaptiveDecrease * ddt, AdaptiveMinTimestep ) 2583 StepControl = 1 2584 2585 CALL Info('ExecSimulation','Reverting to previous timestep as initial guess',Level=8) 2586 DO i=1,nSolvers 2587 Solver => CurrentModel % Solvers(i) 2588 IF ( ASSOCIATED( Solver % Variable % Values ) ) THEN 2589 IF( ASSOCIATED( Solver % Variable % PrevValues ) ) THEN 2590 Solver % Variable % Values = Solver % Variable % PrevValues(:,1) 2591 Solver % Variable % Norm = Solver % Variable % PrevNorm 2592 END IF 2593 END IF 2594 END DO 2595 END IF 2596 END DO 2597 ELSE 2598 CALL SolveEquations( CurrentModel, dt, Transient, & 2599 CoupledMinIter, CoupledMaxIter, SteadyStateReached, RealTimestep ) 2600 RealTimestep = RealTimestep+1 2601 END IF 2602!------------------------------------------------------------------------------ 2603! Save results to disk, if requested 2604!------------------------------------------------------------------------------ 2605 2606 LastSaved = .FALSE. 2607 2608 IF( OutputIntervals(Interval) /= 0 ) THEN 2609 2610 CALL SaveToPost(0) 2611 k = MOD( Timestep-1, OutputIntervals(Interval) ) 2612 2613 IF ( k == 0 .OR. SteadyStateReached ) THEN 2614 DO i=1,nSolvers 2615 Solver => CurrentModel % Solvers(i) 2616 IF ( Solver % PROCEDURE == 0 ) CYCLE 2617 ExecThis = ( Solver % SolverExecWhen == SOLVER_EXEC_AHEAD_SAVE) 2618 When = ListGetString( Solver % Values, 'Exec Solver', GotIt ) 2619 IF ( GotIt ) ExecThis = ( When == 'before saving') 2620 IF( ExecThis ) CALL SolverActivate( CurrentModel,Solver,dt,Transient ) 2621 END DO 2622 2623 CALL SaveCurrent(Timestep) 2624 LastSaved = .TRUE. 2625 2626 DO i=1,nSolvers 2627 Solver => CurrentModel % Solvers(i) 2628 IF ( Solver % PROCEDURE == 0 ) CYCLE 2629 ExecThis = ( Solver % SolverExecWhen == SOLVER_EXEC_AFTER_SAVE) 2630 When = ListGetString( Solver % Values, 'Exec Solver', GotIt ) 2631 IF ( GotIt ) ExecThis = ( When == 'after saving') 2632 IF( ExecThis ) CALL SolverActivate( CurrentModel,Solver,dt,Transient ) 2633 END DO 2634 END IF 2635 END IF 2636!------------------------------------------------------------------------------ 2637 CALL ListPopNameSpace() 2638!------------------------------------------------------------------------------ 2639 2640 maxtime = ListGetCReal( CurrentModel % Simulation,'Real Time Max',GotIt) 2641 IF( GotIt .AND. RealTime() - RT0 > maxtime ) THEN 2642 CALL Info('ElmerSolver','Reached allowed maximum real time, exiting...') 2643 GOTO 100 2644 END IF 2645 2646 exitcond = ListGetCReal( CurrentModel % Simulation,'Exit Condition',GotIt) 2647 IF( GotIt .AND. exitcond > 0.0_dp ) THEN 2648 CALL Info('ElmerSolver','Found a positive exit condition, exiting...') 2649 GOTO 100 2650 END IF 2651 2652!------------------------------------------------------------------------------ 2653 2654 IF ( SteadyStateReached .AND. .NOT. (Transient .OR. Scanning) ) THEN 2655 IF ( Timestep >= CoupledMinIter ) EXIT 2656 END IF 2657 2658!------------------------------------------------------------------------------ 2659 END DO ! timestep within an iterval 2660!------------------------------------------------------------------------------ 2661 2662!------------------------------------------------------------------------------ 2663 END DO ! timestep intervals, i.e. the simulation 2664!------------------------------------------------------------------------------ 2665 2666100 CONTINUE 2667 2668 CALL ListPopNamespace() 2669 2670 DO i=1,nSolvers 2671 Solver => CurrentModel % Solvers(i) 2672 IF ( Solver % PROCEDURE == 0 ) CYCLE 2673 When = ListGetString( Solver % Values, 'Exec Solver', GotIt ) 2674 IF ( GotIt ) THEN 2675 IF ( When == 'after simulation' .OR. When == 'after all' ) THEN 2676 CALL SolverActivate( CurrentModel,Solver,dt,Transient ) 2677 IF (ASSOCIATED(Solver % Variable % Values) ) LastSaved = .FALSE. 2678 END IF 2679 ELSE 2680 IF ( Solver % SolverExecWhen == SOLVER_EXEC_AFTER_ALL ) THEN 2681 CALL SolverActivate( CurrentModel,Solver,dt,Transient ) 2682 IF (ASSOCIATED(Solver % Variable % Values) ) LastSaved = .FALSE. 2683 END IF 2684 END IF 2685 END DO 2686 2687!------------------------------------------------------------------------------ 2688! Always save the last step to output 2689!----------------------------------------------------------------------------- 2690 IF ( .NOT.LastSaved ) THEN 2691 DO i=1,CurrentModel % NumberOfSolvers 2692 Solver => CurrentModel % Solvers(i) 2693 IF ( Solver % PROCEDURE == 0 ) CYCLE 2694 ExecThis = ( Solver % SolverExecWhen == SOLVER_EXEC_AHEAD_SAVE) 2695 When = ListGetString( Solver % Values, 'Exec Solver', GotIt ) 2696 IF ( GotIt ) ExecThis = ( When == 'before saving') 2697 IF( ExecThis ) CALL SolverActivate( CurrentModel,Solver,dt,Transient ) 2698 END DO 2699 2700 CALL SaveToPost(0) 2701 CALL SaveCurrent(Timestep) 2702 2703 DO i=1,CurrentModel % NumberOfSolvers 2704 Solver => CurrentModel % Solvers(i) 2705 IF ( Solver % PROCEDURE == 0 ) CYCLE 2706 ExecThis = ( Solver % SolverExecWhen == SOLVER_EXEC_AFTER_SAVE) 2707 When = ListGetString( Solver % Values, 'Exec Solver', GotIt ) 2708 IF ( GotIt ) ExecThis = ( When == 'after saving') 2709 IF( ExecThis ) CALL SolverActivate( CurrentModel,Solver,dt,Transient ) 2710 END DO 2711 END IF 2712 2713!------------------------------------------------------------------------------ 2714 END SUBROUTINE ExecSimulation 2715!------------------------------------------------------------------------------ 2716 2717 2718!------------------------------------------------------------------------------ 2719!> Saves current timestep to external files. 2720!------------------------------------------------------------------------------ 2721 SUBROUTINE SaveCurrent( CurrentStep ) 2722!------------------------------------------------------------------------------ 2723 INTEGER :: i, j,k,l,n,q,CurrentStep,nlen 2724 TYPE(Variable_t), POINTER :: Var 2725 LOGICAL :: EigAnal, GotIt 2726 CHARACTER(LEN=MAX_NAME_LEN) :: Simul 2727 LOGICAL :: BinaryOutput, SaveAll 2728 2729 Simul = ListGetString( CurrentModel % Simulation, 'Simulation Type' ) 2730 2731 OutputFile = ListGetString(CurrentModel % Simulation,'Output File',GotIt) 2732 2733 2734 IF ( GotIt ) THEN 2735 i = INDEX( OutputFile,'/') 2736 IF( i > 0 ) THEN 2737 CALL Warn('SaveCurrent','> Output File < for restart should not include directory: '& 2738 //TRIM(OutputFile)) 2739 END IF 2740 2741 IF ( ParEnv % PEs > 1 ) THEN 2742 DO i=1,MAX_NAME_LEN 2743 IF ( OutputFile(i:i) == ' ' ) EXIT 2744 END DO 2745 OutputFile(i:i) = '.' 2746 WRITE( OutputFile(i+1:), '(a)' ) TRIM(i2s(ParEnv % MyPE)) 2747 END IF 2748 2749 BinaryOutput = ListGetLogical( CurrentModel % Simulation,'Binary Output',GotIt ) 2750 IF ( .NOT.GotIt ) BinaryOutput = .FALSE. 2751 2752 SaveAll = .NOT.ListGetLogical( CurrentModel % Simulation,& 2753 'Omit unchanged variables in output',GotIt ) 2754 IF ( .NOT.GotIt ) SaveAll = .TRUE. 2755 2756 Mesh => CurrentModel % Meshes 2757 DO WHILE( ASSOCIATED( Mesh ) ) 2758 IF ( Mesh % OutputActive ) THEN 2759 nlen = LEN_TRIM(Mesh % Name ) 2760 IF ( nlen > 0 ) THEN 2761 OutputName = Mesh % Name(1:nlen) // '/' // TRIM(OutputFile) 2762 ELSE 2763 OutputName = OutputFile 2764 END IF 2765 2766 EigAnal = .FALSE. 2767 DO i=1,CurrentModel % NumberOfSolvers 2768 IF ( ASSOCIATED( CurrentModel % Solvers(i) % Mesh, Mesh ) ) THEN 2769 EigAnal = ListGetLogical( CurrentModel % Solvers(i) % Values, & 2770 'Eigen Analysis', GotIt ) 2771 EigAnal = EigAnal .OR. ListGetLogical( CurrentModel % Solvers(i) % Values, & 2772 'Harmonic Analysis', GotIt ) 2773 2774 IF ( EigAnal ) THEN 2775 Var => CurrentModel % Solvers(i) % Variable 2776 IF ( ASSOCIATED(Var % EigenValues) ) THEN 2777 IF ( TotalTimesteps == 1 ) THEN 2778 DO j=1,CurrentModel % Solvers(i) % NOFEigenValues 2779 IF ( CurrentModel % Solvers(i) % Matrix % COMPLEX ) THEN 2780 2781 n = SIZE(Var % Values)/Var % DOFs 2782 DO k=1,n 2783 DO l=1,Var % DOFs/2 2784 q = Var % DOFs*(k-1) 2785 Var % Values(q+l) = REAL(Var % EigenVectors(j,q/2+l)) 2786 Var % Values(q+l+Var % DOFs/2) = AIMAG(Var % EigenVectors(j,q/2+l)) 2787 END DO 2788 END DO 2789 ELSE 2790 Var % Values = REAL( Var % EigenVectors(j,:) ) 2791 END IF 2792 SavedSteps = SaveResult( OutputName, Mesh, & 2793 j, sTime(1), BinaryOutput, SaveAll ) 2794 END DO 2795 ELSE 2796 j = MIN( CurrentStep, SIZE( Var % EigenVectors,1 ) ) 2797 IF ( CurrentModel % Solvers(i) % Matrix % COMPLEX ) THEN 2798 n = SIZE(Var % Values)/Var % DOFs 2799 DO k=1,n 2800 DO l=1,Var % DOFs/2 2801 q = Var % DOFs*(k-1) 2802 Var % Values(q+l) = REAL(Var % EigenVectors(j,q/2+l)) 2803 Var % Values(q+l+Var % DOFs/2) = AIMAG(Var % EigenVectors(j,q/2+l)) 2804 END DO 2805 END DO 2806 ELSE 2807 Var % Values = REAL(Var % EigenVectors(j,:)) 2808 END IF 2809 SavedSteps = SaveResult( OutputName, Mesh, & 2810 CurrentStep, sTime(1), BinaryOutput, SaveAll ) 2811 END IF 2812 Var % Values = 0.0d0 2813 END IF 2814 END IF 2815 END IF 2816 END DO 2817 2818 IF ( .NOT. EigAnal ) THEN 2819 SavedSteps = SaveResult( OutputName,Mesh, NINT(sStep(1)), & 2820 sTime(1), BinaryOutput, SaveAll ) 2821 END IF 2822 END IF 2823 Mesh => Mesh % Next 2824 END DO 2825 ELSE 2826 Mesh => CurrentModel % Meshes 2827 DO WHILE( ASSOCIATED( Mesh ) ) 2828 IF ( Mesh % OutputActive ) & 2829 Mesh % SavesDone=Mesh % SavesDone+1 2830 Mesh => Mesh % Next 2831 END DO 2832 END IF 2833 CALL SaveToPost(CurrentStep) 2834!------------------------------------------------------------------------------ 2835 END SUBROUTINE SaveCurrent 2836!------------------------------------------------------------------------------ 2837 2838!------------------------------------------------------------------------------ 2839!> Saves results file to post processing file of ElmerPost format, if requested. 2840!------------------------------------------------------------------------------ 2841 SUBROUTINE SaveToPost(CurrentStep) 2842!------------------------------------------------------------------------------ 2843 TYPE(Variable_t), POINTER :: Var 2844 LOGICAL :: EigAnal = .FALSE., Found 2845 INTEGER :: i, j,k,l,n,q,CurrentStep,nlen,timesteps,SavedEigenValues 2846 CHARACTER(LEN=MAX_NAME_LEN) :: Simul, SaveWhich 2847 2848 Simul = ListGetString( CurrentModel % Simulation, 'Simulation Type' ) 2849 2850 OutputFile = ListGetString( CurrentModel % Simulation,'Output File',GotIt ) 2851 IF ( Gotit ) THEN 2852 IF ( ParEnv % PEs > 1 ) THEN 2853 DO i=1,MAX_NAME_LEN 2854 IF ( OutputFile(i:i) == ' ' ) EXIT 2855 END DO 2856 OutputFile(i:i) = '.' 2857 WRITE( OutputFile(i+1:), '(a)' ) TRIM(i2s(ParEnv % MyPE)) 2858 END IF 2859 END IF 2860 2861 PostFile = ListGetString( CurrentModel % Simulation,'Post File',GotIt ) 2862 IF( .NOT. GotIt ) RETURN 2863 2864 IF ( ParEnv % PEs > 1 ) THEN 2865 DO i=1,MAX_NAME_LEN 2866 IF ( PostFile(i:i) == ' ' ) EXIT 2867 END DO 2868 PostFile(i:i) = '.' 2869 WRITE( PostFile(i+1:), '(a)' ) TRIM(i2s(ParEnv % MyPE)) 2870 END IF 2871 2872 ! Loop over all meshes 2873 !--------------------------------------------------- 2874 Mesh => CurrentModel % Meshes 2875 DO WHILE( ASSOCIATED( Mesh ) ) 2876 IF ( CurrentStep == 0 .AND. Mesh % SavesDone > 0) THEN 2877 Mesh => Mesh % Next 2878 CYCLE 2879 END IF 2880 2881 ! Check whether this mesh is active for saving 2882 !-------------------------------------------------- 2883 IF ( Mesh % OutputActive ) THEN 2884 nlen = LEN_TRIM(Mesh % Name) 2885 IF ( nlen==0 .OR. FileNameQualified(OutputFile) ) THEN 2886 OutputName = OutputFile 2887 ELSE 2888 OutputName = Mesh % Name(1:nlen)//'/'//TRIM(OutputFile) 2889 END IF 2890 2891 IF ( nlen==0 .OR. FileNameQualified(PostFile) ) THEN 2892 PostName = PostFile 2893 ELSE 2894 Postname = Mesh % Name(1:nlen)//'/'//TRIM(PostFile) 2895 END IF 2896 2897 IF ( ListGetLogical( CurrentModel % Simulation,'Filename Numbering',GotIt) ) THEN 2898 IF( CurrentStep == 0 ) THEN 2899 PostName = NextFreeFilename(PostName) 2900 ELSE 2901 PostName = NextFreeFilename(PostName,LastExisting = .TRUE. ) 2902 END IF 2903 END IF 2904 2905 ! Set the Mesh pointer in the CurrentModel 2906 CALL SetCurrentMesh( CurrentModel, Mesh ) 2907 2908 ! Use number of timesteps or number of eigenmodes 2909 !------------------------------------------------ 2910 EigAnal = .FALSE. 2911 timesteps = TotalTimeSteps 2912 DO i=1,CurrentModel % NumberOfSolvers 2913 IF (ASSOCIATED(CurrentModel % Solvers(i) % Mesh, Mesh)) THEN 2914 EigAnal = ListGetLogical( CurrentModel % & 2915 Solvers(i) % Values, 'Eigen Analysis', GotIt ) 2916 2917 EigAnal = EigAnal .OR. ListGetLogical( CurrentModel % & 2918 Solvers(i) % Values, 'Harmonic Analysis', GotIt ) 2919 2920 IF ( EigAnal ) timesteps = MAX( timesteps, & 2921 CurrentModel % Solvers(i) % NOFEigenValues ) 2922 END IF 2923 END DO 2924 2925 DO i=1,CurrentModel % NumberOfSolvers 2926 IF (ASSOCIATED(CurrentModel % Solvers(i) % Mesh, Mesh)) THEN 2927 EigAnal = ListGetLogical( CurrentModel % & 2928 Solvers(i) % Values, 'Eigen Analysis', GotIt ) 2929 2930 EigAnal = EigAnal .OR. ListGetLogical( CurrentModel % & 2931 Solvers(i) % Values, 'Harmonic Analysis', GotIt ) 2932 2933 IF ( EigAnal ) THEN 2934 SaveWhich = ListGetString( CurrentModel % Solvers(i) % Values, & 2935 'Eigen and Harmonic Solution Output', Found ) 2936 2937 SavedEigenValues = CurrentModel % Solvers(i) % NOFEigenValues 2938 IF( TotalTimesteps > 1 ) THEN 2939! SavedEiegnValues = MIN( CurrentStep, SIZE( Var % EigenVectors,1 ) ) 2940 END IF 2941 2942 DO j=1, SavedEigenValues 2943 Var => Mesh % Variables 2944 DO WHILE(ASSOCIATED(Var)) 2945 IF ( .NOT. ASSOCIATED(Var % EigenValues) ) THEN 2946 Var => Var % Next 2947 CYCLE 2948 END IF 2949 2950 IF ( CurrentModel % Solvers(i) % Matrix % COMPLEX ) THEN 2951 IF(Var % DOFs==1) THEN 2952 Var => Var % Next 2953 CYCLE 2954 END IF 2955 2956 n = SIZE(Var % Values)/Var % DOFs 2957 DO k=1,n 2958 DO l=1,Var % DOFs/2 2959 q = Var % DOFs*(k-1) 2960 Var % Values(q+l) = REAL(Var % EigenVectors(j,q/2+l)) 2961 Var % Values(q+l+Var % DOFs/2) = AIMAG(Var % EigenVectors(j,q/2+l)) 2962 END DO 2963 END DO 2964 2965 ELSE 2966 SELECT CASE( SaveWhich ) 2967 CASE('real part') 2968 Var % Values = Var % EigenVectors(j,:) 2969 CASE('imag part') 2970 Var % Values = AIMAG(Var % EigenVectors(j,:)) 2971 CASE('abs value') 2972 Var % Values = ABS(Var % EigenVectors(j,:)) 2973 CASE('phase angle') 2974 Var % Values = ATAN2(AIMAG(Var % EigenVectors(j,:)), & 2975 REAL(Var % EigenVectors(j,:))) 2976 CASE DEFAULT 2977 Var % CValues => Var % EigenVectors(j,:) 2978 END SELECT 2979 END IF 2980 Var => Var % Next 2981 END DO 2982 2983 IF ( CurrentStep > 0 ) THEN 2984 IF ( Mesh % SavesDone /= 0 ) THEN 2985 IF( TotalTimeSteps == 1 ) THEN 2986 Mesh % SavesDone = j 2987 ELSE 2988 Mesh % SavesDone = CurrentStep 2989 END IF 2990 END IF 2991 CALL WritePostFile( PostName,OutputName, CurrentModel, & 2992 CurrentModel % Solvers(i) % NOFEigenValues, .TRUE. ) 2993 END IF 2994 END DO 2995 EXIT 2996 END IF 2997 END IF 2998 END DO 2999 3000 ! If this mesh has not been saved, then do so 3001 !---------------------------------------------------------------------------- 3002 IF ( .NOT. EigAnal .OR. CurrentStep == 0 ) THEN 3003 CALL WritePostFile( PostName, OutputName, CurrentModel, timesteps, .TRUE. ) 3004 END IF 3005 3006 Var => Mesh % Variables 3007 DO WHILE(ASSOCIATED(Var)) 3008 IF ( ASSOCIATED(Var % EigenValues) ) THEN 3009 Var % Values = 0._dp 3010 Var % CValues => NULL() 3011 END IF 3012 Var => Var % Next 3013 END DO 3014 END IF 3015 Mesh => Mesh % Next 3016 END DO 3017!------------------------------------------------------------------------------ 3018 END SUBROUTINE SaveToPost 3019!------------------------------------------------------------------------------ 3020 3021!------------------------------------------------------------------------------ 3022 END SUBROUTINE ElmerSolver 3023!------------------------------------------------------------------------------ 3024 3025!> \} ElmerLib 3026