TimeModule.F90

Go to the documentation of this file.

!/*****************************************************************************/
! *
! *  Elmer, A Finite Element Software for Multiphysical Problems
! *
! *  Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland
! * 
! *  This library is free software; you can redistribute it and/or
! *  modify it under the terms of the GNU Lesser General Public
! *  License as published by the Free Software Foundation; either
! *  version 2.1 of the License, or (at your option) any later version.
! *
! *  This library is distributed in the hope that it will be useful,
! *  but WITHOUT ANY WARRANTY; without even the implied warranty of
! *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
! *  Lesser General Public License for more details.
! * 
! *  You should have received a copy of the GNU Lesser General Public
! *  License along with this library (in file ../LGPL-2.1); if not, write 
! *  to the Free Software Foundation, Inc., 51 Franklin Street, 
! *  Fifth Floor, Boston, MA  02110-1301  USA
! *
! *****************************************************************************/
!
!/******************************************************************************
! *
! *  ELMER/FEM Solver main program
! *
! ******************************************************************************
! *
! *  Authors: Juha Ruokolainen
! *  Email:   Juha.Ruokolainen@csc.fi
! *  Web:     http://www.csc.fi/elmer
! *  Address: CSC - IT Center for Science Ltd.
! *           Keilaranta 14
! *           02101 Espoo, Finland 
! *
! *  Original Date: 02 Jun 1997
! *
! *****************************************************************************/






!------------------------------------------------------------------------------
!------------------------------------------------------------------------------
   MODULE timemodule
!------------------------------------------------------------------------------

     USE generalmodule

!------------------------------------------------------------------------------
     IMPLICIT NONE
!------------------------------------------------------------------------------
     REAL(KIND=dp) :: ddt
     INTEGER :: timeleft,cum_timestep
     INTEGER, SAVE ::  stepcount=0, RealTimestep
     LOGICAL :: SteadyStateReached=.FALSE.

     REAL(KIND=dp) :: CumTime, MaxErr, AdaptiveLimit, &
           AdaptiveMinTimestep, AdaptiveMaxTimestep, timePeriod
     INTEGER :: SmallestCount, AdaptiveKeepSmallest, StepControl=-1
     LOGICAL :: AdaptiveTime = .TRUE.

     REAL(KIND=dp) :: newtime, prevtime=0, maxtime, exitcond
     REAL(KIND=dp), ALLOCATABLE :: xx(:,:), xxnrm(:), yynrm(:), PrevXX(:,:,:)

    !Something I (Jess) am adding in order to see if I can access
    !the displacements from this level
    TYPE(valuelist_t), POINTER :: SolverParams
    TYPE(variable_t), POINTER :: StressSol
    CHARACTER(LEN=MAX_NAME_LEN) :: Equation
    REAL(KIND=dp), POINTER :: Displacement(:)
    INTEGER :: BoundaryElementCount, ElementCount, body_id
    INTEGER :: bc_id, nElementNodes, nt, nk
    INTEGER, POINTER :: MyIndex(:), MyNodeIndexes(:), MyPerm(:)
    TYPE(element_t), POINTER :: Element, MyCurrentElement
    TYPE(nodes_t), SAVE :: ElementNodes
    TYPE(mesh_t), POINTER :: MyMesh
    LOGICAL :: IsFSI
    DOUBLE PRECISION :: FinalTime, PreviousTime
    !End of Jess stuff

!------------------------------------------------------------------------------
  CONTAINS
!------------------------------------------------------------------------------

  SUBROUTINE updateloads(global,runs)
     USE testobject
     USE lists
  
     IMPLICIT NONE
  
     include 'comf90.h'
  
     TYPE(t_global), POINTER :: global
     INTEGER :: runs
     REAL(KIND=dp), ALLOCATABLE  :: f(:)
     INTEGER :: myn, currentinterval
     TYPE(valuelist_t), POINTER :: ptr
     
     WRITE(6,*) 'MyTimeModule:UpdateLoads: Updating loads on the structure'
    
     IF (myverbosity > 3) THEN
       WRITE(6,*) 'sTime(1) = ', stime(1)
       WRITE(6,*) 'FinalTime = ', finaltime
       WRITE(6,*) 'PreviousTime = ', previoustime
     END IF
 
     !printing the loads as a check
     IF( myverbosity > 3) THEN
       WRITE(*,*) 'global%NodeLoads:'
       DO t = 1, global%nNodes
           WRITE(*,*) global%NodeLoads(3*(t-1) + 1), global%NodeLoads(3*(t-1) + 2), &
           global%NodeLoads(3*(t-1) + 3)
       END DO
     END IF
     IF( myverbosity > 3) THEN
       WRITE(*,*) 'global%PreviousLoads:'
       DO t = 1, global%nNodes
         DO j =1,3
           WRITE(*,*) global%PreviousLoads(j,t)
         END DO
       END DO
     END IF
  
     !Setting the loads accessible by Elmer to the interpolated value
     !for the current time
     WRITE(*,*) 'MyTimeModule:UpdateLoads: Time Step Summary '
     WRITE(*,*) '    sTime(1)     = ', stime(1)
     WRITE(*,*) '    PreviousTime = ', previoustime
     WRITE(*,*) '    FinalTime    = ', finaltime
     !Original
     DO t = 1, global%nNodes
       j=1
       DO j =1,3
         IF (global%IsLoadFreeNode(t)) THEN
           !Masoud : setting node loads to zero for beam end nodes in the
           !Heron-Turek problem
           currentmodel % NodeLoadsPass(j,t) = 0.d0
           WRITE(*,*) '   Node #',t,' is load free'
           !Masoud End
         ELSE
           currentmodel % NodeLoadsPass(j,t) = global%PreviousLoads(j,t) + &
             (global%NodeLoads(3*(t-1) + j) - global%PreviousLoads(j,t))&
             *(stime(1) - previoustime)/(finaltime - previoustime)
         END IF
       END DO
     END DO
     !Original End

     !Masoud: Sending load update only
     !DO t = 1, global%nNodes
     !  j=1
     !  DO j =1,3
     !    CurrentModel % NodeLoadsPass(j,t) = global%NodeLoads(3*(t-1) + j) - global%PreviousLoads(j,t)
     !  END DO
     !END DO
     !Masoud End
    
 
     !printing the loads as a check
     IF( myverbosity > 3) THEN
       WRITE(*,*) 'CurrentModel % NodeLoadsPass:'
       DO t = 1, global%nNodes
         DO j =1,3
           WRITE(*,*) currentmodel % NodeLoadsPass(j,t)
         END DO
       END DO
     END IF
 
     WRITE(6,*) 'MyTimeModule:UpdateLoads: Finishing'
  
  END SUBROUTINE updateloads
!----------------------------------------------------------------------
SUBROUTINE timestepper(global,runs)

     USE generalmodule
     USE testobject

!------------------------------------------------------------------------------
     IMPLICIT NONE
!------------------------------------------------------------------------------
     include 'comf90.h'

     TYPE(t_global), POINTER :: global
     INTEGER :: runs

     DO interval = 1, timeintervals
        stepcount = stepcount + timesteps(interval)
     END DO

     IF (myverbosity > 3) THEN
       WRITE(*,*) 'In TimeStepper'
       WRITE(*,*) 'TimeIntervals = ', timeintervals 
       WRITE(*,*) 'global%nNodes = ', global%nNodes
     END IF

     cum_timestep = 0
     ddt = 0.0d0
     DO interval = 1,timeintervals

     IF (myverbosity > 3) THEN
       WRITE(*,*) 'interval = ', interval
     END IF 
!------------------------------------------------------------------------------
       IF ( transient .OR. scanning ) THEN
         dt = timestepsizes(interval,1)
       ELSE
         dt = 1
       END IF
!------------------------------------------------------------------------------
!      go trough number of timesteps within an interval
!------------------------------------------------------------------------------
       timeperiod = listgetcreal(currentmodel % Simulation, 'Time Period',gotit)
       IF(.NOT.gotit) timeperiod = huge(timeperiod)


       realtimestep = 1
       WRITE(*,*) 'MyTimeModule:TimeStepper: Timesteps(interval) = ', timesteps(interval)
       DO timestep = 1,timesteps(interval)

         cum_timestep = cum_timestep + 1
         sstep(1) = cum_timestep

         dtfunc = listgetconstreal( currentmodel % Simulation, &
                  'Timestep Function', gotit)
         IF(gotit) THEN
           CALL warn('ExecSimulation','Obsolite keyword > Timestep Function < , use > Timestep Size < instead')
         ELSE   
           dtfunc = listgetcreal( currentmodel % Simulation, &
                  'Timestep Size', gotit)
         END IF
         IF(gotit) dt = dtfunc
       
!------------------------------------------------------------------------------
         stime(1) = stime(1) + dt
         speriodic(1) = stime(1)
         DO WHILE(speriodic(1) > timeperiod)
           speriodic(1) = speriodic(1) - timeperiod 
         END DO

         ! Move the old timesteps one step down the ladder
         IF(timestep > 1 .OR. interval > 1) THEN
           DO i = SIZE(sprevsizes,2),2,-1
             sprevsizes(1,i) = sprevsizes(1,i-1)
           END DO
           sprevsizes(1,1) = ssize(1)
         END IF 
         ssize(1) = dt

         sinterval(1) = interval
         IF (.NOT. transient ) steadyit(1) = steadyit(1) + 1
!------------------------------------------------------------------------------
         IF ( parenv % MyPE == 0 ) THEN
           CALL info( 'MAIN', ' ', level=3 )
           CALL info( 'MAIN', '-------------------------------------', level=3 )

           IF ( transient .OR. scanning ) THEN
             WRITE( message, * ) 'Time: ',trim(i2s(cum_timestep)),'/', &
                   trim(i2s(stepcount)), stime(1)
             CALL info( 'MAIN', message, level=3 )

             newtime= realtime()

             IF( cum_timestep > 1 ) THEN
               maxtime = listgetconstreal( currentmodel % Simulation,'Real Time Max',gotit)
               IF( gotit ) THEN
                  WRITE( message,'(A,F8.3)') 'Fraction of real time left: ',&
                              1.0_dp-realtime() / maxtime
               ELSE             
                 timeleft = nint((stepcount-(cum_timestep-1))*(newtime-prevtime)/60._dp);
                 IF (timeleft > 120) THEN
                   WRITE( message, *) 'Estimated time left: ', &
                     trim(i2s(timeleft/60)),' hours.'
                 ELSE IF(timeleft > 60) THEN
                   WRITE( message, *) 'Estimated time left: 1 hour ', &
                     trim(i2s(mod(timeleft,60))), ' minutes.'
                 ELSE IF(timeleft >= 1) THEN
                   WRITE( message, *) 'Estimated time left: ', &
                     trim(i2s(timeleft)),' minutes.'
                 ELSE
                   WRITE( message, *) 'Estimated time left: less than a minute.'
                 END IF
               END IF
               CALL info( 'MAIN', message, level=3 )
             END IF
             prevtime = newtime
           ELSE
             WRITE( message, * ) 'Steady state iteration: ',cum_timestep
             CALL info( 'MAIN', message, level=3 )
           END IF

           CALL info( 'MAIN', '-------------------------------------', level=3 )
           CALL info( 'MAIN', ' ', level=3 )
         END IF

!------------------------------------------------------------------------------
!        Call UpdateLoads to the get linearly interpolated load value at 
!        the current time
!------------------------------------------------------------------------------
         IF (myverbosity > 3) THEN
           WRITE(6,*) 'timestepper timestep = ', timestep
           WRITE(6,*) 'sTime(1) = ', stime(1)
           WRITE(6,*) 'dtfunc = ', dtfunc
           WRITE(6,*) 'dt = ', dt
           WRITE(*,*) 'TimeStepper Calling UpdateLoads'
         END IF
         CALL updateloads(global,runs)
!------------------------------------------------------------------------------
!        Solve any and all governing equations in the system
!------------------------------------------------------------------------------
         adaptivetime = listgetlogical( currentmodel % Simulation, &
                  'Adaptive Timestepping', gotit )

         IF ( transient .AND. adaptivetime ) THEN 
            adaptivelimit = listgetconstreal( currentmodel % Simulation, &
                        'Adaptive Time Error', gotit )
 
            IF ( .NOT. gotit ) THEN 
               WRITE( message, * ) 'Adaptive Time Limit must be given for' // &
                        'adaptive stepping scheme.'
               CALL fatal( 'ElmerSolver', message )
            END IF

            adaptivemaxtimestep = listgetconstreal( currentmodel % Simulation, &
                     'Adaptive Max Timestep', gotit )
            IF ( .NOT. gotit ) adaptivemaxtimestep =  dt
            adaptivemaxtimestep =  min(adaptivemaxtimestep, dt)

            adaptivemintimestep = listgetconstreal( currentmodel % Simulation, &
                     'Adaptive Min Timestep', gotit )

            adaptivekeepsmallest = listgetinteger( currentmodel % Simulation, &
                       'Adaptive Keep Smallest', gotit, minv=0  )

            n = currentmodel % NumberOfSolvers
            j = 0
            k = 0
            DO i=1,n
               solver => currentmodel % Solvers(i)
               IF ( ASSOCIATED( solver % Variable  % Values ) ) THEN
                  IF ( ASSOCIATED( solver % Variable % PrevValues ) ) THEN
                     j = max( j, SIZE( solver % Variable % PrevValues,2 ) )
                  END IF
                  k = max( k, SIZE( solver % Variable % Values ) )
               END IF
            END DO
            ALLOCATE( xx(n,k), yynrm(n), xxnrm(n), prevxx( n,k,j ) )

            cumtime = 0.0d0
            IF ( ddt == 0.0d0 .OR. ddt > adaptivemaxtimestep ) ddt = adaptivemaxtimestep

            s = stime(1) - dt
            smallestcount = 0
            DO WHILE( cumtime < dt-1.0d-12 )
               ddt = min( dt - cumtime, ddt )

               DO i=1,currentmodel % NumberOFSolvers
                  solver => currentmodel % Solvers(i)
                  IF ( ASSOCIATED( solver % Variable % Values ) ) THEN
                     n = SIZE( solver % Variable % Values )
                     xx(i,1:n) = solver % Variable % Values
                     xxnrm(i) = solver % Variable % Norm
                     IF ( ASSOCIATED( solver % Variable % PrevValues ) ) THEN
                        DO j=1,SIZE( solver % Variable % PrevValues,2 )
                           prevxx(i,1:n,j) = solver % Variable % PrevValues(:,j)
                        END DO
                     END IF
                  END IF
               END DO

               stime(1) = s + cumtime + ddt
               ssize(1) = ddt
               CALL solveequations( currentmodel, ddt, transient, &
                 coupledminiter, coupledmaxiter, steadystatereached, realtimestep )


               maxerr = listgetconstreal( currentmodel % Simulation, &
                          'Adaptive Error Measure', gotit )

               DO i=1,currentmodel % NumberOFSolvers
                  solver => currentmodel % Solvers(i)
                  IF ( ASSOCIATED( solver % Variable % Values ) ) THEN
                     n = SIZE(solver % Variable % Values)
                     yynrm(i) = solver % Variable % Norm
                     solver % Variable % Values = xx(i,1:n)
                     IF ( ASSOCIATED( solver % Variable % PrevValues ) ) THEN
                        DO j=1,SIZE( solver % Variable % PrevValues,2 )
                           solver % Variable % PrevValues(:,j) = prevxx(i,1:n,j)
                        END DO
                     END IF
                  END IF
               END DO

               sstep(1) = ddt / 2
               stime(1) = s + cumtime + ddt/2
               CALL solveequations( currentmodel, ddt/2, transient, &
                  coupledminiter, coupledmaxiter, steadystatereached, realtimestep )
               stime(1) = s + cumtime + ddt
               CALL solveequations( currentmodel, ddt/2, transient, &
                  coupledminiter, coupledmaxiter, steadystatereached, realtimestep )

               maxerr = abs( maxerr - listgetconstreal( currentmodel % Simulation, &
                           'Adaptive Error Measure', gotit ) )

               IF ( .NOT. gotit ) THEN
                  maxerr = 0.0d0
                  DO i=1,currentmodel % NumberOFSolvers
                     solver => currentmodel % Solvers(i)
                     IF ( ASSOCIATED( solver % Variable % Values ) ) THEN
                        IF ( yynrm(i) /= solver % Variable % Norm ) THEN
                           maxerr = max(maxerr,abs(yynrm(i)-solver % Variable % Norm)/yynrm(i))
                        END IF
                     END IF
                  END DO
               END IF

               IF ( maxerr < adaptivelimit .OR. ddt <= adaptivemintimestep ) THEN
                 cumtime = cumtime + ddt
                 realtimestep = realtimestep+1
                 IF ( smallestcount >= adaptivekeepsmallest .OR. stepcontrol > 0 ) THEN
                    ddt = min( 2*ddt, adaptivemaxtimestep )
                    stepcontrol   = 1
                    smallestcount = 0
                  ELSE
                    stepcontrol   = 0
                    smallestcount = smallestcount + 1
                  END IF
               ELSE
                  DO i=1,currentmodel % NumberOFSolvers
                     solver => currentmodel % Solvers(i)
                     IF ( ASSOCIATED( solver % Variable % Values ) ) THEN
                        n = SIZE(solver % Variable % Values)
                        solver % Variable % Norm = xxnrm(i)
                        solver % Variable % Values = xx(i,1:n)
                        IF ( ASSOCIATED( solver % Variable % PrevValues ) ) THEN
                           DO j=1,SIZE( solver % Variable % PrevValues,2 )
                              solver % Variable % PrevValues(:,j) = prevxx(i,1:n,j)
                           END DO
                        END IF
                     END IF
                  END DO
                  ddt = ddt / 2
                  stepcontrol = -1
               END IF
               WRITE(*,'(a,3e20.12)') 'Adaptive(cum,ddt,err): ', cumtime, ddt, maxerr
            END DO
            ssize(1) = dt
            stime(1) = s + dt
  
            DEALLOCATE( xx, xxnrm, yynrm, prevxx )
         ELSE ! Adaptive timestepping
            !IF (MyVerbosity > 3) THEN
              WRITE(*,*) 'MyTimeModule:TimeStepper: Calling SolveEquations'
              WRITE(*,*) '   Transient = ', transient
            !END IF
            CALL solveequations( currentmodel, dt, transient, &
              coupledminiter, coupledmaxiter, steadystatereached, realtimestep )
            realtimestep = realtimestep+1
         END IF

!------------------------------------------------------------------------------
!        Save results to disk, if requested
!------------------------------------------------------------------------------
         lastsaved = .false.
         IF( outputintervals(interval) /= 0 ) THEN

           CALL savetopost(0)
           k = mod( timestep-1, outputintervals(interval) )
           IF ( k == 0 .OR. steadystatereached ) THEN
            
             DO i=1,currentmodel % NumberOfSolvers
               solver => currentmodel % Solvers(i)
               IF ( solver % PROCEDURE == 0 ) cycle
               execthis = ( solver % SolverExecWhen == solver_exec_ahead_save)
               when = listgetstring( solver % Values, 'Exec Solver', gotit )
               IF ( gotit ) execthis = ( when == 'before saving') 
               IF( execthis ) CALL solveractivate( currentmodel,solver,dt,transient )
             END DO 

             CALL savecurrent(timestep)
             lastsaved = .true.

             DO i=1,currentmodel % NumberOfSolvers
               solver => currentmodel % Solvers(i)
               IF ( solver % PROCEDURE == 0 ) cycle
               execthis = ( solver % SolverExecWhen == solver_exec_after_save)
               when = listgetstring( solver % Values, 'Exec Solver', gotit )
               IF ( gotit ) execthis = ( when == 'after saving') 
               IF( execthis ) CALL solveractivate( currentmodel,solver,dt,transient )
             END DO 
           END IF
         END IF
!------------------------------------------------------------------------------

         maxtime = listgetcreal( currentmodel % Simulation,'Real Time Max',gotit)
         IF( gotit .AND. realtime() - rt0 > maxtime ) THEN
            CALL info('ElmerSolver','Reached allowed maximum real time, exiting...')
            goto 100
         END IF

         exitcond = listgetcreal( currentmodel % Simulation,'Exit Condition',gotit)
         IF( gotit .AND. exitcond > 0.0_dp ) THEN
            CALL info('ElmerSolver','Found a positive exit condition, exiting...')
            goto 100
         END IF
         
!------------------------------------------------------------------------------

         IF ( steadystatereached .AND. .NOT. (transient .OR. scanning) ) THEN
            IF ( timestep >= coupledminiter ) THEN
              WRITE(6,*) 'SteadyStateReached, exiting'
              EXIT
            END IF
         END IF

!------------------------------------------------------------------------------
       END DO ! timestep within an iterval
!------------------------------------------------------------------------------

!------------------------------------------------------------------------------
     END DO ! timestep intervals, i.e. the simulation

     IF (myverbosity > 3) THEN
       write(6,*) 'LastSaved = ', lastsaved
     END IF

     runs = 1
      
100  RETURN

END SUBROUTINE timestepper
!------------------------------------------------------------------------------
  END MODULE timemodule
!------------------------------------------------------------------------------