RFLU_ModNewtonKrylov.F90

Go to the documentation of this file.

! *********************************************************************
! * Rocstar Simulation Suite                                          *
! * Copyright@2015, Illinois Rocstar LLC. All rights reserved.        *
! *                                                                   *
! * Illinois Rocstar LLC                                              *
! * Champaign, IL                                                     *
! * www.illinoisrocstar.com                                           *
! * sales@illinoisrocstar.com                                         *
! *                                                                   *
! * License: See LICENSE file in top level of distribution package or *
! * http://opensource.org/licenses/NCSA                               *
! *********************************************************************
! *********************************************************************
! * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,   *
! * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   *
! * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          *
! * NONINFRINGEMENT.  IN NO EVENT SHALL THE CONTRIBUTORS OR           *
! * COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER       *
! * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,   *
! * Arising FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE    *
! * USE OR OTHER DEALINGS WITH THE SOFTWARE.                          *
! *********************************************************************
! ******************************************************************************
!
! Purpose: Collection of routines related to Newton-Krylov schemes.
!
! Description: None.
!
! Notes: None. 
!
! ******************************************************************************
!
! $Id: RFLU_ModNewtonKrylov.F90,v 1.16 2008/12/06 08:44:22 mtcampbe Exp $
!
! Copyright: (c) 2005-2006 by the University of Illinois
!
! ******************************************************************************

MODULE rflu_modnewtonkrylov

  USE moddatatypes
  USE modglobal, ONLY: t_global
  USE modparameters
  USE moderror
  USE moddatastruct, ONLY: t_region
  USE modgrid, ONLY: t_grid
  USE modbndpatch, ONLY: t_patch 
  USE modmixture, ONLY: t_mixt,t_mixt_input
  USE modmpi
  
  USE rflu_moddimensions
  USE rflu_modforcesmoments
  USE rflu_modgeometry
  USE rflu_modmpi
  USE rflu_modpatchcoeffs
  USE rflu_modpetscnewtonkrylov  
  USE rflu_modprobes
  USE rflu_modreadwriteflow
  USE rflu_modreadwritegrid
  USE rflu_modreadwritegridspeeds
  USE rflu_modweights  
  
  USE modinterfaces, ONLY: integratesourcetermsmp, &
                           rflu_computegridspeeds, &
                           rflu_computeintegralvalues, & 
                           rflu_checkgridspeeds, &
                           rflu_decideneedbgradface, &
                           rflu_decideprint, & 
                           rflu_decidewrite, & 
                           rflu_explicitmultistage, &
                           rflu_getdeformationwrapper, & 
                           rflu_minimumtimestep, &
                           rflu_movegridwrapper, &
                           rflu_printchangeinfo, &  
                           rflu_printgridinfo, &
                           rflu_printflowinfowrapper, &
                           rflu_printflowinfo, &
                           rflu_printwriteconvergence, &
                           rflu_residualnorm, &
                           rflu_timestepinviscid, &
                           rflu_timestepviscous, &
                           rflu_writerestartinfo, &
                           rflu_writestatsfileoles, &
                           rungekuttamp, &
                           writeprobe, &
                           writetotalmass, &
                           rflu_checkpositivity, &
                           rflu_checkvalidity, &
                           rflu_setdependentvars

  IMPLICIT NONE

  include 'mpif.h'

! DEBUG
!#ifdef GENX
!  INCLUDE 'rocmanf90.h'
!#endif
! END DEBUG

  PRIVATE
  PUBLIC :: rflu_nk_timestepping
  
! ******************************************************************************
! Declarations and definitions
! ******************************************************************************  
      
  CHARACTER(CHRLEN) :: & 
    RCSIdentString = '$RCSfile: RFLU_ModNewtonKrylov.F90,v $ $Revision: 1.16 $' 
              
! ******************************************************************************
! Routines
! ******************************************************************************

  CONTAINS
  
  
  
  


! ******************************************************************************
!
! Purpose: Integrate the governing equations in time.
!
! Description: None.
!
! Input: 
!   dTimeSystem         Total solution time (unsteady flow)
!   dIterSystem         Total number of iterations (steady flow)
!   regions             Data for all grid regions
!
! Output: 
!   regions             Data for all grid regions
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_nk_timestepping(dTimeSystem,dIterSystem,regions)

  USE moddatatypes
  USE modparameters
  USE moderror  
  USE modglobal, ONLY: t_global    
  USE modbndpatch, ONLY: t_patch
  USE moddatastruct, ONLY: t_region
  USE modgrid, ONLY: t_grid  
  USE modmixture, ONLY: t_mixt_input    
  USE modmpi
  
  USE rflu_modpetscnewtonkrylov
  
  IMPLICIT NONE

! ******************************************************************************
! Include PETSc headers
! ******************************************************************************

#include "finclude/petsc.h"
#include "finclude/petscvec.h"
#include "finclude/petscmat.h"
#include "finclude/petscksp.h"
#include "finclude/petscpc.h"
#include "finclude/petscsnes.h"

! ******************************************************************************
! INTERFACE blocks for PETSc SNES callback definition functions
! ******************************************************************************

  INTERFACE
    SUBROUTINE snessetfunction(snes,r,fun,ctx,ierr)
      USE moddatastruct, ONLY: t_region
      snes :: snes
      vec :: r
      EXTERNAL :: fun
      TYPE(t_region),POINTER :: ctx
      INTEGER :: ierr
    END SUBROUTINE
    SUBROUTINE snessetjacobian(snes,A,pA,fun,ctx,ierr)
      USE moddatastruct, ONLY: t_region
      snes :: snes
      mat :: a, pa
      EXTERNAL :: fun
      TYPE(t_region),POINTER :: ctx
      INTEGER :: ierr
    END SUBROUTINE
    SUBROUTINE matfdcoloringsetfunctionsnes(matfd,fun,ctx,ierr)
      USE moddatastruct, ONLY: t_region
      matfdcoloring :: matfd
      EXTERNAL :: fun
      TYPE(t_region),POINTER :: ctx
      INTEGER :: ierr
    END SUBROUTINE
  END INTERFACE

! ******************************************************************************
! Declarations and definitions
! ******************************************************************************  

! ==============================================================================
! Arguments
! ==============================================================================

  LOGICAL :: doprint,doprobe,dowrite
  LOGICAL :: finished = .false., movegrid
  INTEGER :: ditersystem
  REAL(RFREAL) :: dtimesystem
  TYPE(t_region), POINTER :: pregion,regions(:)

! ==============================================================================
! Locals
! ==============================================================================

  INTEGER :: gtemp,ic,ierr,ipatch,ireg,iter,iv,ltemp
  TYPE(t_global), POINTER :: global
  TYPE(t_mixt_input), POINTER :: pmixtinput
  TYPE(t_patch), POINTER :: ppatch
  snesconvergedreason reason
  petsclogdouble time

! DEBUG
!#ifdef GENX
!  DOUBLE PRECISION, DIMENSION(MAN_INTEG_SIZE) :: integ
!#endif
! END DEBUG

! TEMPORARY
  INTEGER :: isub,ivl,ipatch,nsub
  REAL(RFREAL) :: insub, isubterm
  TYPE(t_patch), POINTER :: ppatch
  TYPE(t_grid), POINTER :: pgrid, pgridold, pgridold2
! END TEMPORARY

! ******************************************************************************
! Start, set pointers and variables
! ******************************************************************************

  global  => regions(1)%global
  pregion => regions(1)

  nsub  = 100
  insub = 1.0_rfreal/nsub

  CALL registerfunction(global,'RFLU_NK_TimeStepping',&
  'RFLU_ModNewtonKrylov.F90')
  
! ==============================================================================
! Determine whether have moving grids
! ==============================================================================

  movegrid = .false.

  DO ireg = 1,global%nRegionsLocal
    IF ( regions(ireg)%mixtInput%moveGrid .EQV. .true. ) THEN 
      movegrid = .true.
    END IF ! regions
  END DO ! iReg

! ==============================================================================
! Set SNES callback functions
! ==============================================================================

  CALL snessetfunction(pregion%snes,pregion%r, &
       rflu_petsc_formresidual,pregion,ierr)
  CALL snessetjacobian(pregion%snes,pregion%A,pregion%preA, &
       rflu_petsc_formjacobian,pregion,ierr)
  CALL matfdcoloringsetfunctionsnes(pregion%fdcolor, &
       rflu_petsc_formresidualfirstorder,pregion,ierr)  

! ==============================================================================
! Initialize old variables as current variables
! ==============================================================================

  pregion%mixt%cvOld(:,:)   = pregion%mixt%cv(:,:)
  pregion%mixt%cvOld1(:,:)  = pregion%mixt%cv(:,:)
  pregion%mixt%cvOld2(:,:)  = pregion%mixt%cv(:,:)
  pregion%gridOld%vol(:)    = pregion%grid%vol(:)
  pregion%gridOld2%vol(:)   = pregion%grid%vol(:)
  pregion%gridOld%xyz(:,:)  = pregion%grid%xyz(:,:)
  pregion%gridOld2%xyz(:,:) = pregion%grid%xyz(:,:)

! ******************************************************************************
! Loop over iterations/time steps
! ****************************************************************************** 

! DEBUG
  IF ( global%myProcid == 0 ) THEN
    CALL petscgettime(time,ierr)
    print*,'PETSC TIME 1 : ',time
  END IF ! global%myProcid
! END DEBUG

  finished = .false.

  DO

! ==============================================================================
!   Update iteration counter for steady flow.
! ==============================================================================

    IF ( global%flowType == flow_steady ) THEN
      global%currentIter      = global%currentIter      + 1
      global%iterSinceRestart = global%iterSinceRestart + 1
    END IF ! global%flowType

! ==============================================================================
!   Compute time step.
! ==============================================================================

    DO ireg = 1,global%nRegionsLocal  
      pregion => regions(ireg)

      IF ( pregion%mixtInput%flowModel == flow_euler ) THEN  
        CALL rflu_timestepinviscid(pregion)
      ELSE 
        CALL rflu_timestepviscous(pregion)
      END IF ! pRegion
    END DO ! iReg

    global%dtmin = global%dtImposed

! ==============================================================================
!   Move or generate new grid
! ==============================================================================

    IF ( ( global%flowType == flow_unsteady ) .AND.  & 
         ( global%iterSinceRestart > 0 ) ) THEN
      IF ( movegrid .EQV. .true. ) THEN

! ------------------------------------------------------------------------------
!       Get the displacements from GENX (or other source).
! ------------------------------------------------------------------------------

        CALL rflu_getdeformationwrapper(regions)

! ------------------------------------------------------------------------------
!       Store the old grid.
! ------------------------------------------------------------------------------

        DO ireg = 1,global%nRegionsLocal 
          pgrid       => regions(ireg)%grid
          pgridold    => regions(ireg)%gridOld
          pgridold2   => regions(ireg)%gridOld2    
           
          DO ic = 1,pgrid%nCellsTot ! Explicit copy to avoid ASCI White problem
            pgridold2%vol(ic) = pgridold%vol(ic)
          END DO ! ic  
           
          DO iv = 1,pgrid%nVertTot ! Explicit copy to avoid ASCI White problem    
            pgridold2%xyz(xcoord,iv) = pgridold%xyz(xcoord,iv)
            pgridold2%xyz(ycoord,iv) = pgridold%xyz(ycoord,iv)
            pgridold2%xyz(zcoord,iv) = pgridold%xyz(zcoord,iv)            
          END DO ! iv 
        END DO ! iReg
        
        DO ireg = 1,global%nRegionsLocal 
          pgrid    => regions(ireg)%grid
          pgridold => regions(ireg)%gridOld    
          
          DO ic = 1,pgrid%nCellsTot ! Explicit copy to avoid ASCI White problem
            pgridold%vol(ic) = pgrid%vol(ic)
          END DO ! ic  
           
          DO iv = 1,pgrid%nVertTot ! Explicit copy to avoid ASCI White problem    
            pgridold%xyz(xcoord,iv) = pgrid%xyz(xcoord,iv)
            pgridold%xyz(ycoord,iv) = pgrid%xyz(ycoord,iv)
            pgridold%xyz(zcoord,iv) = pgrid%xyz(zcoord,iv)            
          END DO ! iv 
        END DO ! iReg

! ------------------------------------------------------------------------------
!       We will be substepping Mesquite, so start with a smaller displacement.
! ------------------------------------------------------------------------------

        DO ireg = 1,global%nRegionsLocal 
          pgrid => regions(ireg)%grid    

          DO ipatch = 1,pgrid%nPatches
            ppatch => regions(ireg)%patches(ipatch)
             
            DO ivl = 1,ppatch%nBVert
              ppatch%dXyz(xcoord,ivl) = insub*ppatch%dXyz(xcoord,ivl)
              ppatch%dXyz(ycoord,ivl) = insub*ppatch%dXyz(ycoord,ivl)
              ppatch%dXyz(zcoord,ivl) = insub*ppatch%dXyz(zcoord,ivl)
            END DO ! ivl
          END DO ! iPatch         
        END DO ! iReg

! ------------------------------------------------------------------------------
!       Substep displacements and smoothing to avoid inverting elements.
! ------------------------------------------------------------------------------

        DO isub = 1,nsub
          isubterm = (1.0_rfreal*isub)/(1.0_rfreal*isub-1.0_rfreal)

! ------- Displace the mesh by steadily increasing displacements until the 
!         original displacement is reached.
           
          DO ireg = 1,global%nRegionsLocal 
            pgrid => regions(ireg)%grid    

            DO ipatch = 1,pgrid%nPatches
              ppatch => regions(ireg)%patches(ipatch)
                 
              DO ivl = 1,ppatch%nBVert 
                IF ( isub > 1 ) THEN
                  ppatch%dXyz(xcoord,ivl) = isubterm*ppatch%dXyz(xcoord,ivl)
                  ppatch%dXyz(ycoord,ivl) = isubterm*ppatch%dXyz(ycoord,ivl)
                  ppatch%dXyz(zcoord,ivl) = isubterm*ppatch%dXyz(zcoord,ivl)
                END IF ! iSub
              END DO ! ivl
            END DO ! iPatch         
          END DO ! iReg

! ------- Have Mesquite smooth the grid at each substep. NOTE:  
!         This function will not store the old grid. That must be done before 
!         Mesquite substepping begins.

          CALL rflu_movegridwrapper(regions)
        END DO ! iSub

! ------------------------------------------------------------------------------
!       Finish by generating geometry and grid speeds for the new grid.
! ------------------------------------------------------------------------------

        DO ireg = 1,global%nRegionsLocal
          pregion => regions(ireg)

          CALL rflu_buildgeometry(pregion)
          CALL rflu_computegridspeeds(pregion)

          IF ( global%checkLevel == check_high ) THEN
            CALL rflu_checkgridspeeds(pregion)
          END IF ! global%checkLevel
        END DO ! iReg
      END IF ! moveGrid
    END IF ! global%flowType 

! ==============================================================================
!   Recompute weights
! ==============================================================================

    IF ( global%flowType == flow_unsteady ) THEN
      IF ( movegrid .EQV. .true. ) THEN
        DO ireg=1,global%nRegionsLocal
          pregion => regions(ireg)
          pmixtinput => pregion%mixtInput

          IF ( pmixtinput%spaceOrder > 1 ) THEN
            CALL rflu_computewtsc2cwrapper(pregion,pmixtinput%spaceOrder-1)
          END IF ! pMixtInput%spaceOrder

          IF ( pmixtinput%flowModel == flow_navst ) THEN
            CALL rflu_computewtsf2cwrapper(pregion,pmixtinput%spaceOrder-1)
          END IF ! pMixtInput%flowModel
            
          DO ipatch = 1,pregion%grid%nPatches
            ppatch => pregion%patches(ipatch)

            IF ( rflu_decideneedbgradface(pregion,ppatch) .EQV. .true. ) THEN
              CALL rflu_computewtsbf2cwrapper(pregion,ppatch, &
                                              ppatch%spaceOrder)
            END IF ! RFLU_DecideNeedBGradFace
          END DO ! iPatch  

        END DO ! iReg
      END IF ! moveGrid
    END IF ! global%flowType

! ==============================================================================
!   Relocate probes
! ==============================================================================

    IF ( global%flowType == flow_unsteady ) THEN
      IF ( movegrid .EQV. .true. ) THEN
        IF ( global%nProbes > 0 ) THEN
          DO ireg = 1,global%nRegionsLocal
            pregion => regions(ireg)
            CALL rflu_findprobecells(pregion)
          END DO ! iReg

          CALL rflu_printprobeinfo(global)
        END IF ! global
      END IF ! moveGrid
    END IF ! global%flowType

! ==============================================================================
!   Compute new solution.
! ==============================================================================

    IF ( global%flowType == flow_steady ) THEN

! ------------------------------------------------------------------------------
!   For steady-state problems, simply compute the new solution
! ------------------------------------------------------------------------------

      global%forceX  = 0.0_rfreal
      global%forceY  = 0.0_rfreal
      global%forceZ  = 0.0_rfreal
      global%massIn  = 0.0_rfreal
      global%massOut = 0.0_rfreal
      pregion => regions(1)
      pregion%mixt%cvOld(:,:)  = pregion%mixt%cv(:,:)
      pregion%irkStep = 1
      CALL snessolve(pregion%snes,petsc_null_object,pregion%x,ierr)
      CALL snesgetiterationnumber(pregion%snes,iter,ierr)
!      CALL RFLU_SetDependentVars(pRegion,1,pRegion%grid%nCells)

      CALL rflu_mpi_isendwrapper(pregion)
!        CALL RFLU_MPI_CopyWrapper(regions)
      CALL rflu_setdependentvars(pregion,1,pregion%grid%nCells)
      CALL rflu_mpi_recvwrapper(pregion)
      CALL rflu_mpi_clearrequestwrapper(pregion)
      CALL rflu_setdependentvars(pregion,pregion%grid%nCells+1,pregion%grid%nCellsTot)

    ELSE

! ------------------------------------------------------------------------------
!   For transient problems, first update the iteration variables for the
!   pseudo steady-state problem.
! ------------------------------------------------------------------------------

      pregion => regions(1)

      iter = global%iterSinceRestart
      global%currentIter      = 0
      global%iterSinceRestart = 0

      pregion%mixt%cvOld2(:,:)  = pregion%mixt%cvOld1(:,:)
      pregion%mixt%cvOld1(:,:)  = pregion%mixt%cv(:,:)

      DO

! ------------------------------------------------------------------------------
!   Repeatedly solve the pseudo steady-state problem until the residual drops
!   below the residual tolerance.
! ------------------------------------------------------------------------------

        global%currentIter      = global%currentIter      + 1
        global%iterSinceRestart = global%iterSinceRestart + 1

        DO ireg = 1,global%nRegionsLocal  
          pregion => regions(ireg)
          IF ( pregion%mixtInput%flowModel == flow_euler ) THEN  
            CALL rflu_timestepinviscid(pregion)
          ELSE 
            CALL rflu_timestepviscous(pregion)
          END IF ! pRegion
        END DO ! iReg
        global%dtmin = global%dtImposed

        pregion%mixt%cvOld(:,:)  = pregion%mixt%cv(:,:)

        global%forceX  = 0.0_rfreal
        global%forceY  = 0.0_rfreal
        global%forceZ  = 0.0_rfreal
        global%massIn  = 0.0_rfreal
        global%massOut = 0.0_rfreal

        pregion%irkStep = 1

        CALL snessolve(pregion%snes,petsc_null_object,pregion%x,ierr)
!        CALL RFLU_SetDependentVars(pRegion,1,pRegion%grid%nCells)

        CALL rflu_mpi_isendwrapper(pregion)
!        CALL RFLU_MPI_CopyWrapper(regions)
        CALL rflu_setdependentvars(pregion,1,pregion%grid%nCells)
        CALL rflu_mpi_recvwrapper(pregion)
        CALL rflu_mpi_clearrequestwrapper(pregion)
        CALL rflu_setdependentvars(pregion,pregion%grid%nCells+1,pregion%grid%nCellsTot)
        
        CALL rflu_checkvalidity(pregion)
        CALL rflu_checkpositivity(pregion)

        CALL rflu_residualnorm(regions)    

! DEBUG
        IF ( global%myProcid == 0 ) THEN
          print*,'pseudo-steady info:',global%currentIter, &
               global%residual/global%resInit,global%resTol
        END IF ! global%myProcid
! END DEBUG

        IF ( global%residual/global%resInit <= global%resTol ) THEN 
          EXIT
        END IF ! global%residual
      END DO
    END IF ! global%flowType

! ==============================================================================
!   Check for convergence problems
! ==============================================================================

   CALL snesgetconvergedreason(pregion%snes,reason,ierr)
   !print*,'SNES REASON = ',reason
   IF(reason == snes_diverged_function_count) THEN
     global%warnCounter = global%warnCounter + 1
     print*,'*** WARNING ***  SNES FUNCTION COUNT EXCEEDED.  ', &
          'COPYING cvOld INTO cv TO PRESERVE PARTIALLY-CONVERGED SOLUTION.'
     pregion%mixt%cv(:,:) = pregion%mixt%cvOld(:,:)
   ENDIF ! reason
   IF((reason /= snes_diverged_max_it).AND.(reason /= snes_converged_pnorm_relative).AND. &
        (reason /= snes_diverged_function_count)) THEN
     global%warnCounter = global%warnCounter + 1
     print*,'*** WARNING ***  SNES REASON DOES NOT INDICATE CONTINUATION OF SNES STEPPING: ',reason
   ENDIF

! DEBUG
!   CALL RFLU_PrintFlowInfo(pRegion)
! END DEBUG

! ==============================================================================
!   Check validity and positivity
! ==============================================================================

   CALL rflu_checkvalidity(pregion)
   CALL rflu_checkpositivity(pregion)

! ==============================================================================
!   Get statistics
! ==============================================================================

#ifdef STATS
    CALL getstatistics(regions)
#endif

! ==============================================================================
!   Reset global%timeSince* variables. NOTE this must be done right before
!   the update of the various time variables, otherwise calling the routines
!   RFLU_Decide* from within rungeKuttaMP or explicitMultiStage will not work
!   correctly.
! ==============================================================================

    IF ( rflu_decideprint(global) .EQV. .true. ) THEN 
      IF ( global%flowType == flow_unsteady ) THEN    
        IF ( global%iterSinceRestart > 1 ) THEN 
          global%timeSincePrint = 0.0_rfreal
        END IF ! global%iterSinceRestart
      END IF ! global%flowType      
    END IF ! RFLU_DecidePrint

    IF ( rflu_decidewrite(global) .EQV. .true. ) THEN 
      IF ( global%flowType == flow_unsteady ) THEN    
        global%timeSinceWrite = 0.0_rfreal
      END IF ! global%flowType      
    END IF ! RFLU_DecideWrite

    IF ( rflu_decidewriteprobes(global) .EQV. .true. ) THEN 
      IF ( global%flowType == flow_unsteady ) THEN    
        IF ( global%iterSinceRestart > 1 ) THEN 
          global%timeSinceProbe = 0.0_rfreal
        END IF ! global%iterSinceRestart 
      END IF ! global%flowType      
    END IF ! RFLU_DecideWriteProbes  

! ==============================================================================
!   Update times for unsteady flow
! ==============================================================================

    IF ( global%flowType == flow_unsteady ) THEN
      global%currentTime      = global%currentTime      + global%dtImposed    
      global%timeSinceRestart = global%timeSinceRestart + global%dtImposed
      
      global%timeSincePrint = global%timeSincePrint + global%dtImposed
      global%timeSinceWrite = global%timeSinceWrite + global%dtImposed
      global%timeSinceProbe = global%timeSinceProbe + global%dtImposed  
      
      global%iterSinceRestart = iter + 1
    END IF ! global%flowType

! ==============================================================================
!   Decide whether to print/write convergence, data, and probes
! ==============================================================================

    doprint = rflu_decideprint(global)
    dowrite = rflu_decidewrite(global)
    doprobe = rflu_decidewriteprobes(global)

! ==============================================================================
!   Check for end of timestepping
! ==============================================================================
   
   IF ( global%flowType == flow_unsteady ) THEN    
      IF ( global%timeSinceRestart >= dtimesystem ) THEN 
        finished = .true.
      END IF ! global%timeSinceRestart
    ELSE
      IF ( (doprint .EQV. .true.) .OR. (global%iterSinceRestart >= ditersystem) ) THEN 
        CALL rflu_residualnorm(regions)    
        IF ( (global%iterSinceRestart >= ditersystem) .OR. &
             (global%residual/global%resInit <= global%resTol) ) THEN 
          finished = .true.
        END IF ! global%iterSinceRestart      
      END IF ! doPrint
    END IF ! global%flowType

! ==============================================================================
!   Write convergence (file & screen) and total mass (file) history
! ==============================================================================

    IF ( (doprint .EQV. .true.) .OR. (finished .EQV. .true.) ) THEN
      CALL rflu_printwriteconvergence(global)
      
#ifndef GENX      
      IF ( movegrid .EQV. .true. ) THEN
        CALL rflu_computeintegralvalues(regions) 
        CALL writetotalmass(regions)
      END IF ! moveGrid
#endif    
      
      DO ireg = 1,global%nRegionsLocal
        IF ( regions(ireg)%mixtInput%spaceDiscr == discr_opt_les ) THEN 
          CALL rflu_writestatsfileoles(global)
        END IF ! mixtInput
      END DO ! iReg
    END IF ! doPrint

! ==============================================================================
!   Compute forces and mass flow
! ==============================================================================

    IF ( global%forceFlag .EQV. .true. ) THEN 
      IF ( (dowrite .EQV. .true.) .OR. (finished .EQV. .true.) ) THEN 
        DO ireg = 1,global%nRegionsLocal
          pregion => regions(ireg)

          CALL rflu_computelocalforcesmoments(pregion)     
        END DO ! iReg

        CALL rflu_computeglobalforcesmoments(regions)

        pregion => regions(1)

        IF ( pregion%global%myProcid == masterproc ) THEN 
          CALL rflu_printglobalforcesmoments(pregion)
          CALL rflu_writeglobalforcesmoments(pregion)
        END IF ! pRegion
      END IF ! doWrite
    END IF ! global%forceFlag

! ==============================================================================
!   Store probe data
! ==============================================================================

    IF ( global%nProbes > 0 ) THEN
      IF ( doprobe .EQV. .true. ) THEN 
        DO ireg = 1,global%nRegionsLocal
          CALL writeprobe(regions,ireg)
        END DO ! iReg
      END IF ! doProbe
    END IF ! global

! DEBUG
#ifdef GENX
!    CALL RFLU_WriteGridWrapper(pRegion)
!    CALL RFLU_WriteFlowWrapper(pRegion)
!    CALL RFLU_ComputeIntegralValues(regions,integ) 
!    print*,global%myProcid,global%currentTime,'777 777',integ(MAN_INTEG_MASS),integ(MAN_INTEG_VOL)
#endif
! END DEBUG

#ifndef GENX
! ==============================================================================
!   Store flow field (and grid if moving). Write restart info file after 
!   flow (and grid) files so that should those not be written due to reaching
!   the time limit, the restart file will not contain the time level of the 
!   incomplete flow (and grid) files.
! ==============================================================================

    IF ( (dowrite .EQV. .true.) .AND. (finished .EQV. .false.) ) THEN            
      DO ireg=1,global%nRegionsLocal
        pregion => regions(ireg)

        CALL rflu_writedimensionswrapper(pregion,write_dimens_mode_maybe)
                
        IF ( movegrid .EQV. .true. ) THEN        
          CALL rflu_writegridwrapper(pregion)
          CALL rflu_writegridspeedswrapper(pregion)
          
          IF ( global%myProcid == masterproc .AND. &
               global%verbLevel > verbose_none ) THEN           
            CALL rflu_printgridinfo(pregion)
          END IF ! global%myProcid
        END IF ! moveGrid
                       
        CALL rflu_writeflowwrapper(pregion)
        CALL rflu_writepatchcoeffswrapper(pregion)
        
#ifdef PLAG
        CALL plag_writesurfstatswrapper(pregion)
#endif  
        
        IF ( global%myProcid == masterproc .AND. &
             global%verbLevel > verbose_none ) THEN          
! TEMPORARY
!          CALL RFLU_PrintFlowInfoWrapper(pRegion)
! END TEMPORARY
          
          IF ( global%verbLevel > verbose_low ) THEN 
            CALL rflu_printchangeinfo(pregion)
          END IF ! global%verbLevel
        END IF ! global%myProcid        
      END DO ! iReg
      
      CALL rflu_writerestartinfo(global)      
    END IF ! doWrite

#ifdef STATS
! ==============================================================================
!   Output statistics
! ==============================================================================

    IF ( (dowrite .EQV. .true.) .AND. (finished .EQV. .false.) .AND. &
         (global%doStat == active) ) THEN
      IF (global%myProcid==masterproc .AND. &
          global%verbLevel/=verbose_none) THEN
        WRITE(stdout,'(A)') solver_name,'Saving statistics ...'
      ENDIF
      
      DO ireg = 1,global%nRegionsLocal
        pregion => regions(ireg)
        CALL rflu_writestat(pregion)
      END DO ! iReg
    END IF ! iReg
#endif
#endif

! ==============================================================================
!   If run finished, update GENX buffers and exit
! ==============================================================================

    IF ( finished .EQV. .true. ) THEN 
#ifdef GENX
       DO ireg = 1,global%nRegionsLocal
         CALL rflu_putboundaryvalues(regions(ireg))
       END DO ! iReg
       
       global%timeStamp = global%currentTime       
#endif   
! DEBUG
       IF ( global%myProcid == 0 ) THEN
         CALL petscgettime(time,ierr)
         print*,'PETSC TIME 2 : ',time
       END IF ! global%myProcid
! END DEBUG
       EXIT
    END IF ! finished

! ==============================================================================
!   End of timestepping loop
! ==============================================================================

  END DO ! loop over time/iterations

! ******************************************************************************
! End
! ****************************************************************************** 

  CALL deregisterfunction(global)

END SUBROUTINE rflu_nk_timestepping






  
  
  
  
! ******************************************************************************
! End
! ******************************************************************************
  
END MODULE rflu_modnewtonkrylov


! ******************************************************************************
!
! RCS Revision history:
!
! $Log: RFLU_ModNewtonKrylov.F90,v $
! Revision 1.16  2008/12/06 08:44:22  mtcampbe
! Updated license.
!
! Revision 1.15  2008/11/19 22:17:33  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.14  2006/08/19 15:39:10  mparmar
! Added use of RFLU_DecideNeedBGradFace and pPatch%spaceOrder
!
! Revision 1.13  2006/04/07 16:04:02  haselbac
! Adapted to changes in bf2c wts computation
!
! Revision 1.12  2006/04/07 15:19:20  haselbac
! Removed tabs
!
! Revision 1.11  2006/04/07 14:49:41  haselbac
! Adapted to changes in f and bf diff modules and new WENO module
!
! Revision 1.10  2006/03/25 21:54:01  haselbac
! Cosmetics only
!
! Revision 1.9  2006/03/22 14:58:28  hdewey2
! Fixed bug by adding additional virtual cell communication after SNESSolve calls.
!
! Revision 1.8  2006/03/09 14:07:42  haselbac
! Now call wrapper routine for F2C weights
!
! Revision 1.7  2006/02/08 21:23:39  hdewey2
! Added disp ramping and changed grid storage
!
! Revision 1.6  2006/01/06 22:12:44  haselbac
! Added call to cell grad wrapper routine
!
! Revision 1.5  2005/10/25 19:38:47  haselbac
! Reordered modules, added IF on forceFlag
!
! Revision 1.4  2005/09/22 17:15:26  hdewey2
! Added dual timestepping functionality to solve implicit time-dependent problems.
!
! Revision 1.3  2005/08/02 18:19:43  hdewey2
! Added actual procedures
!
! Revision 1.2  2005/05/19 18:19:00  haselbac
! Cosmetics only
!
! Revision 1.1  2005/05/16 21:12:45  haselbac
! Initial revision
!
! ******************************************************************************