RFLU_ModPETScNewtonKrylov.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. 
!
! ******************************************************************************

MODULE rflu_modpetscnewtonkrylov

  USE moddatatypes
  USE modglobal, ONLY: t_global
  USE modparameters
  USE moderror
  USE moddatastruct, ONLY: t_region
  USE modgrid, ONLY: t_grid
  USE modmixture, ONLY: t_mixt,t_mixt_input
  USE modmpi

  IMPLICIT NONE

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

  PRIVATE
  PUBLIC :: rflu_petsc_createvectors, & 
            rflu_petsc_createjacobian, &
            rflu_petsc_destroyvectors, & 
            rflu_petsc_destroyjacobian, &
            rflu_petsc_formresidual, &
            rflu_petsc_formresidualfirstorder, &
            rflu_petsc_formjacobian, &
            rflu_petsc_createcoloring, &
            rflu_petsc_createapplicationordering, &
            rflu_petsc_destroyapplicationordering, &
            rflu_petsc_getdtscale

  
! ******************************************************************************
! Declarations and definitions
! ******************************************************************************  
      
  CHARACTER(CHRLEN) :: & 
    RCSIdentString = '$RCSfile: RFLU_ModPETScNewtonKrylov.F90,v $ $Revision: 1.9 $' 
              
! ******************************************************************************
! Routines
! ******************************************************************************

  CONTAINS
  
  
  
  


! ******************************************************************************
!
! Purpose: Create and allocate memory for vectors needed by PETSc
!
! Description: None.
!
! Input: 
!   pRegion             Data for a grid region
!
! Output: 
!   x                   Solution vector (initial solution)
!   r                   Callback vector (zeros)
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_petsc_createvectors(pRegion)

  USE moddatastruct, ONLY: t_region
  
  IMPLICIT NONE

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

#include "finclude/petsc.h"
#include "finclude/petscvec.h"

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

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

  TYPE(t_region), POINTER :: pregion

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

  INTEGER :: ierr, ndof
  petscscalar :: zero = 0.0
  TYPE(t_global),POINTER :: global  

! ******************************************************************************
! Start
! ******************************************************************************

  global => pregion%global
  CALL registerfunction(global,'RFLU_PETSC_CreateVectors',&
  'RFLU_ModPETScNewtonKrylov.F90')

! ******************************************************************************
! Count the number of degrees of freedom
! ******************************************************************************

  ndof = (cv_mixt_ener - cv_mixt_dens + 1) * pregion%grid%nCells

! ******************************************************************************
! Create vectors
! ******************************************************************************

  CALL veccreatempiwitharray(petsc_comm_world,ndof,petsc_decide, &
        pregion%mixt%cv,pregion%x,ierr)
  CALL vecduplicate(pregion%x,pregion%r,ierr)
  CALL vecset(pregion%r,zero,ierr)

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_petsc_createvectors




! ******************************************************************************
!
! Purpose: Create and allocate memory for matrices and solver needed by PETSc
!
! Description: None.
!
! Input: 
!   pRegion        data for grid region
!
! Output: 
!   A              the shell jacobian matrix for matrix-free operations
!   pA             the explicitely formed preconditioning jacobian matrix
!   snes           the PETSc SNES solver context
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_petsc_createjacobian(pRegion)

  USE moddatastruct, ONLY: t_region
  
  IMPLICIT NONE

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

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

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

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

  TYPE(t_region), POINTER :: pregion

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

  INTEGER :: ierr, ndof, gndof, maxitsnes, maxitksp
  ksp :: ksp, subksp
  pc :: pc, subpc
  petscscalar :: snestol, ksptol, snestrtol
  TYPE(t_global), POINTER :: global

  maxitsnes = 1
  maxitksp = 10000
  snestol = 1.0e-32   ! tight to avoid SNES thinking its converged
  ksptol = 1.0e-3     ! loose to facilitate fast solving
  snestrtol = 1.0e-32 ! tight to avoid SNES TR thinking its converged

! ******************************************************************************
! Start
! *****************************************************************************

  global => pregion%global
  CALL registerfunction(global,'RFLU_PETSC_CreateJacobian',&
  'RFLU_ModPETScNewtonKrylov.F90')

! ******************************************************************************
! Create the Application Ordering context
! ******************************************************************************

  IF ( global%nProcAlloc > 1 ) THEN
    CALL rflu_petsc_createapplicationordering(pregion)
  END IF ! global%nProcAlloc

! ******************************************************************************
! Count the number of degrees of freedom
! ******************************************************************************

  ndof = (cv_mixt_ener - cv_mixt_dens + 1) * pregion%grid%nCells

! ******************************************************************************
! Set PETSc Trust Region Options
! ******************************************************************************

  CALL petscoptionssetvalue('-snes_tr_delta0','1.e64',ierr)

! ******************************************************************************
! Find the global number of dof
! ******************************************************************************

  CALL mpi_allreduce(ndof,gndof,1,mpi_integer,mpi_sum,mpi_comm_world,ierr)

! ******************************************************************************
! Create matrices
! ******************************************************************************

! ==============================================================================
! Create SNES solver context
! ==============================================================================

  CALL snescreate(petsc_comm_world,pregion%snes,ierr)

! ==============================================================================
! Create jacobian matrix and preconditioning matrix
! ==============================================================================

  CALL matcreatesnesmf(pregion%snes,pregion%x,pregion%A,ierr)
  CALL matcreate(petsc_comm_world,pregion%preA,ierr)
  CALL matsetsizes(pregion%preA,ndof,ndof,gndof,gndof,ierr)
  CALL matsetfromoptions(pregion%preA,ierr)

! ==============================================================================
! Set solver parameters
! ==============================================================================

! ------------------------------------------------------------------------------
! Nonlinear solver parameters
! ------------------------------------------------------------------------------

  CALL snessettype(pregion%snes,snestr,ierr)
  CALL snessettolerances(pregion%snes,snestol,petsc_default_double_precision, &
          petsc_default_double_precision,maxitsnes,petsc_default_integer,ierr)
  CALL snessettrustregiontolerance(pregion%snes,snestrtol,ierr)
  CALL snesgetksp(pregion%snes,ksp,ierr)

! ------------------------------------------------------------------------------
! Linear solver & preconditioner parameters
! ------------------------------------------------------------------------------

  CALL kspsettype(ksp,kspgmres,ierr)
  CALL kspgetpc(ksp,pc,ierr)
  CALL pcsettype(pc,pcasm,ierr)
  CALL kspsettolerances(ksp,ksptol,petsc_default_double_precision, &
        petsc_default_double_precision,maxitksp,ierr)
  CALL kspsetinitialguessnonzero(ksp,petsc_true,ierr)
  CALL kspgmressetrestart(ksp,30,ierr)
  CALL kspgmressetpreallocatevectors(ksp,ierr)

! ------------------------------------------------------------------------------
! ASM Sub solver & preconditioner parameters
! ------------------------------------------------------------------------------

  CALL petscoptionssetvalue('-sub_pc_type','ilu',ierr)
  CALL petscoptionssetvalue('-sub_pc_ilu_levels','1',ierr)
  CALL petscoptionssetvalue('-sub_ksp_type','preonly',ierr)

! ------------------------------------------------------------------------------
! Finalize solver parameters
! ------------------------------------------------------------------------------

  CALL snessetfromoptions(pregion%snes,ierr)

! ==============================================================================
! Prefill Jacobian and Create Jacobian Coloring
! ==============================================================================

  IF ( global%myProcid == 0 ) THEN
     WRITE(*,*) 'Please wait while the Jacobian is prefilled and colored.  This may take a while...'
  END IF ! global%myProcid
  CALL rflu_petsc_createcoloring(pregion)

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_petsc_createjacobian

  
! ******************************************************************************
!
! Purpose: Destroy the vectors that were created for PETSc
!
! Description: None.
!
! Input: 
!   pRegion           data for grid region
!
! Output: 
!   none
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_petsc_destroyvectors(pRegion)

  USE moddatastruct, ONLY: t_region  

  IMPLICIT NONE

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

#include "finclude/petsc.h"
#include "finclude/petscvec.h"

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

  INTEGER :: ierr
  TYPE(t_region), POINTER :: pregion
  TYPE(t_global), POINTER :: global

! ******************************************************************************
! Start
! *****************************************************************************

  global => pregion%global
  CALL registerfunction(global,'RFLU_PETSC_DestroyVectors',&
  'RFLU_ModPETScNewtonKrylov.F90')

! ******************************************************************************
! Destroy the vectors
! ******************************************************************************

  CALL vecdestroy(pregion%x,ierr)
  CALL vecdestroy(pregion%r,ierr)

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_petsc_destroyvectors


! ******************************************************************************
!
! Purpose: Destroy the matrices and solver contexts that were created for PETSc
!
! Description: None.
!
! Input:
!   pRegion             data for grid region 
!
! Output: 
!   none
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_petsc_destroyjacobian(pRegion)
  
  USE moddatastruct, ONLY: t_region 

  IMPLICIT NONE

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

#include "finclude/petsc.h"
#include "finclude/petscmat.h"
#include "finclude/petscsnes.h"

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

  INTEGER :: ierr
  TYPE(t_region), POINTER :: pregion
  TYPE(t_global), POINTER :: global

! ******************************************************************************
! Start
! *****************************************************************************

  global => pregion%global
  CALL registerfunction(global,'RFLU_PETSC_DestroyJacobian',&
  'RFLU_ModPETScNewtonKrylov.F90')

! ******************************************************************************
! Destroy the matrices and solver context
! ******************************************************************************
  
  IF ( global%nProcAlloc == 1 ) THEN
    CALL matfdcoloringdestroy(pregion%fdcolor,ierr)
  END IF ! global%nProcAlloc
  CALL matdestroy(pregion%A,ierr)
  CALL matdestroy(pregion%preA,ierr)
  CALL snesdestroy(pregion%snes,ierr)

! ******************************************************************************
! Destroy the Application Ordering context
! ******************************************************************************

  IF ( global%nProcAlloc > 1 ) THEN
    CALL rflu_petsc_destroyapplicationordering(pregion)
  END IF ! global%nProcAlloc

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_petsc_destroyjacobian






! ******************************************************************************
!
! Purpose: Explicitely forms the preconditioning Jacobian matrix for PETSc.
!
! Description: None.
!
! Input: 
!   snes            the PETSc SNES solver context
!   v               the vector at which to compute the jacobian
!   flag            PETSc flag denoting matrix structure
!   dummy           PETSc dummy context
!
! Output:
!   J               matrix-free jacobian matrix (is not modified)
!   pJ              explicitely formed preconditioning jacobian matrix
!   ierr            PETSc error code
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_petsc_formjacobian(snes,v,J,pJ,flag,pRegion,ierr)

  USE moddatastruct, ONLY: t_region
  
  IMPLICIT NONE

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

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

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

  snes :: snes
  mat :: j, pj
  vec :: v
  INTEGER :: ierr, order
  matstructure :: flag
  TYPE(t_region),POINTER :: pregion
  TYPE(t_global), POINTER :: global
  REAL(kind=8) :: val
  INTEGER :: icg, ieq, neq, row
! DEBUG
  INTEGER :: n, i, k, low, high, nz, gnz
  petscoffset xx_i
  petscscalar xx_v(1)
! END DEBUG

! ******************************************************************************
! Start
! ******************************************************************************

  global => pregion%global
  CALL registerfunction(global,'RFLU_PETSC_FormJacobian',&
  'RFLU_ModPETScNewtonKrylov.F90')

! ******************************************************************************
! Compute the preconditioning Jacobian matrix
! ******************************************************************************

  flag = different_nonzero_pattern
  CALL snesdefaultcomputejacobiancolor(snes,v,j,pj,flag, &
      pregion%fdcolor,ierr)

! DEBUG
!
!  IF ( global%currentIter == 2 ) THEN
!   print*,'writing jacobian to file...'
!   if(global%myProcid == 0) then
!     OPEN(unit=123,file='jacobian_computational_0.m')
!   else
!    OPEN(unit=123,file='jacobian_computational_1.m')
!   endif
!   CALL VecGetSize(v,n,ierr)
!   CALL VecGetOwnershipRange(pRegion%x,low,high,ierr)
!   DO i = low, high-1
!     DO k = 1, n
!       CALL MatGetValues(pJ,1,i,1,k-1,val,ierr)
!       if(val /= 0.0) then
!          write(123,*) 'A(',i+1,',',k,') = ',val,';'
!       endif
!     ENDDO
!   ENDDO
!   close(123)
!   print*,'done writing jacobian to file'
!!   CALL VecView(v,PETSC_VIEWER_STDOUT_WORLD,ierr)
!   stop
!  END IF
!
!  IF ( global%myProcid == 0 ) THEN
!    DO i = 1, 5
!      DO k = 1, 5
!        CALL MatGetValues(pJ,1,i-1,1,k-1,val,ierr)
!        print*,'JAC',i,k,val
!      END DO
!    END DO
!  END IF
!  STOP
!
! END DEBUG

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_petsc_formjacobian




! ******************************************************************************
!
! Purpose: Returns the residual of the nonlinear equations to PETSc.
!
! Description: None.
!
! Input: 
!   x                  vector at which to compute the residual
!   dummy              PETSc dummy input argument
!
! Output: 
!   f                  residual vector
!   ierr               PETSc error code
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_petsc_formresidual(snes,x,f,pRegion,ierr)

  USE rflu_modresidual
  USE rflu_modmpi

  USE moddatastruct, ONLY: t_region

  USE modinterfaces, ONLY: rflu_setdependentvars, &
                           rflu_checkpositivity

  IMPLICIT NONE

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

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

! ******************************************************************************
! Definitions and declarations
! ******************************************************************************

  snes snes
  vec x, f
  INTEGER ierr, counter, icg, ieq, errorflag
  petscoffset xx_i, ff_i
  petscscalar xx_v(1),ff_v(1),temp
  REAL(KIND=8),ALLOCATABLE,DIMENSION(:,:) :: cv, dv
  TYPE(t_region),POINTER :: pregion
  TYPE(t_global), POINTER :: global
  
! ******************************************************************************
! Start
! ******************************************************************************

  global => pregion%global
  CALL registerfunction(global,'RFLU_PETSC_FormResidual',&
  'RFLU_ModPETScNewtonKrylov.F90')

! ******************************************************************************
! Get the vectors
! ******************************************************************************

  CALL vecgetarray(x,xx_v,xx_i,ierr)
  CALL vecgetarray(f,ff_v,ff_i,ierr)

! ******************************************************************************
! Backup the original values for cv and dv
! ******************************************************************************

  ALLOCATE(cv(cv_mixt_dens:cv_mixt_ener,1:pregion%grid%nCellsTot),stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_allocate,__line__,'cv')
  END IF ! global%errorFlag

  ALLOCATE(dv(dv_mixt_pres:dv_mixt_soun,1:pregion%grid%nCellsTot),stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_allocate,__line__,'dv')
  END IF ! global%errorFlag

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

! ******************************************************************************
! Compute copy the given vector into cv and check for validity
! ******************************************************************************

  counter = 0
  DO icg = 1, pregion%grid%nCells
    DO ieq = cv_mixt_dens, cv_mixt_ener
      counter = counter + 1
      pregion%mixt%cv(ieq,icg) = xx_v(xx_i+counter)
    END DO ! ieq
  END DO ! icg

! ******************************************************************************
! Communicate virtual cell values & compute dependent variables
! ******************************************************************************

  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_checkpositivity(pregion)

! ******************************************************************************
! Compute the residual
! ******************************************************************************

  CALL rflu_res_computeresidual(pregion)

! ******************************************************************************
! Copy the residual array into the PETSc vector
! ******************************************************************************

  counter = 0

  DO icg = 1,pregion%grid%nCells
    DO ieq = cv_mixt_dens,cv_mixt_ener
      counter = counter + 1
      IF ( global%flowType == flow_unsteady ) THEN 
        temp = ( 1.5_rfreal*pregion%mixt%cv(ieq,icg) * pregion%grid%vol(icg) &
               - 2.0_rfreal*pregion%mixt%cvOld1(ieq,icg) * pregion%gridOld%vol(icg) &
               + 0.5_rfreal*pregion%mixt%cvOld2(ieq,icg) * pregion%gridOld2%vol(icg) &
               ) / global%dtImposed
        ff_v(ff_i+counter) = pregion%mixt%rhs(ieq,icg) &
          + pregion%grid%vol(icg)/pregion%dt(icg)*(pregion%mixt%cv(ieq,icg)-cv(ieq,icg)) &
          + temp
      ELSE
        ff_v(ff_i+counter) = pregion%mixt%rhs(ieq,icg) &
          + pregion%grid%vol(icg)/pregion%dt(icg)*(pregion%mixt%cv(ieq,icg)-cv(ieq,icg))
      END IF ! global%flowType
    END DO ! ieq
  END DO ! icg

! ******************************************************************************
! Restore original values for cv and dv
! ******************************************************************************

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

  DEALLOCATE(cv,stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_deallocate,__line__,'cv')
  END IF !global%error

  DEALLOCATE(dv,stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_deallocate,__line__,'dv')
  END IF !global%error

! ******************************************************************************
! Put the vectors back together
! ******************************************************************************

  CALL vecrestorearray(x,xx_v,xx_i,ierr)
  CALL vecrestorearray(f,ff_v,ff_i,ierr)

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_petsc_formresidual








! ******************************************************************************
!
! Purpose: Returns the residual of the nonlinear equations to PETSc with first
!          order spacial accuracy.
!
! Description: None.
!
! Input: 
!   x                  vector at which to compute the residual
!   dummy              PETSc dummy input argument
!
! Output: 
!   f                  residual vector
!   ierr               PETSc error code
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_petsc_formresidualfirstorder(snes,x,f,pRegion,ierr)

  USE moddatastruct, ONLY: t_region

  IMPLICIT NONE

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

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

! ******************************************************************************
! Definitions and declarations
! ******************************************************************************

  snes snes
  vec x, f
  TYPE(t_region),POINTER :: pregion
  TYPE(t_global), POINTER :: global
  INTEGER :: ierr, oldspaceorder
  
! ******************************************************************************
! Start
! ******************************************************************************

  global => pregion%global
  CALL registerfunction(global,'RFLU_PETSC_FormResidualFirstOrder',&
  'RFLU_ModPETScNewtonKrylov.F90')

! ==============================================================================
! Set the spacial order of accuracy to first order
! ==============================================================================

  oldspaceorder = pregion%mixtInput%spaceOrder
  pregion%mixtInput%spaceOrder = 1

! ==============================================================================
! Get the residual
! ==============================================================================

  CALL rflu_petsc_formresidual(snes,x,f,pregion,ierr)

! ==============================================================================
! Set the spacial order of accuracy to what it was to start with
! ==============================================================================

  pregion%mixtInput%spaceOrder = oldspaceorder

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_petsc_formresidualfirstorder







! ******************************************************************************
!
! Purpose: Creates the Jacobian matrix coloring contexts needed by PETSc, and
!          preallocates and prefills the Jacobian matrix.
!
! Description: None.
!
! Input:
!   pRegion             data for grid region 
!
! Output: 
!   none
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_petsc_createcoloring(pRegion)

  USE rflu_modcoloring

  USE moddatastruct, ONLY: t_region 
  USE rflu_modresidual, ONLY: rflu_getresidualsupport1, &
                              rflu_getresidualsupport2

  IMPLICIT NONE

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

#include "finclude/petsc.h"
#include "finclude/petscmat.h"
#include "finclude/petscsnes.h"
#include "finclude/petscis.h"
#include "finclude/petscao.h"

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

  TYPE(t_region), POINTER :: pregion
  TYPE(t_global), POINTER :: global
  INTEGER :: i, j, k, ifg, icg, ieq, irs, c1, c2, ierr, gndof, ndof, dof, color, neq
  INTEGER :: low, high, coloroffset, errorflag, rssizemax, rssize
  INTEGER,ALLOCATABLE,DIMENSION(:) :: nnzon, nnzoff, colors, pcolors, rs
  petscscalar :: val
  iscoloring :: iscolor
  INTEGER,ALLOCATABLE,DIMENSION(:) :: rows, cols
  petscscalar,ALLOCATABLE,DIMENSION(:,:) :: vals

! ******************************************************************************
! Start
! ******************************************************************************

  global => pregion%global
  CALL registerfunction(global,'RFLU_PETSC_CreateColoring',&
  'RFLU_ModPETScNewtonKrylov.F90')

! ==============================================================================
! Create the Connectivity Array
! ==============================================================================

  rssizemax = 1000

  ALLOCATE(rs(1:rssizemax),stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_allocate,__line__,'rs')
  END IF ! global%errorFlag

! ==============================================================================
! Preallocate the Jacobian Matrix
! ==============================================================================
 
! ------------------------------------------------------------------------------
! Get misc dimensions
! ------------------------------------------------------------------------------

  neq = cv_mixt_ener - cv_mixt_dens + 1
  CALL vecgetsize(pregion%x,gndof,ierr)
  CALL vecgetlocalsize(pregion%x,ndof,ierr)

! ------------------------------------------------------------------------------
! Allocate preallocation variables
! ------------------------------------------------------------------------------

  ALLOCATE(nnzon(1:ndof),stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_allocate,__line__,'nnzon')
  END IF ! global%errorFlag

  ALLOCATE(nnzoff(1:ndof),stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_allocate,__line__,'nnzoff')
  END IF ! global%errorFlag

! ------------------------------------------------------------------------------
! Calculate number of nonzeros in each row in the Jacobian. NOTE must be careful
! with Fortran and C-style numbering - nnzon uses Fortran numbering...
! ------------------------------------------------------------------------------

  nnzon(:) = neq
  nnzoff(:) = 0

  DO c1 = 1, pregion%grid%nCells
! TEMPORARY
!    IF ( pRegion%mixtInput%spaceOrder == DISCR_ORDER_1 ) THEN
      CALL rflu_getresidualsupport1(pregion,c1,rs,rssizemax,rssize)
!    ELSE
!      CALL RFLU_GetResidualSupport2(pRegion,c1,rs,rsSizeMax,rsSize)
!    END IF ! pRegion%mixtInput%spaceOrder
! END TEMPORARY

    DO irs = 1, rssize
      c2 = rs(irs)

      DO ieq = cv_mixt_dens, cv_mixt_ener
        j = neq*(c1-1) + ieq 
        nnzon(j) = nnzon(j) + neq

        IF ( c2 <= pregion%grid%nCells ) THEN
          j = neq*(c2-1) + ieq
          nnzon(j) = nnzon(j) + neq
        ELSE
          j = neq*(c1-1) + ieq 
          nnzoff(j) = nnzoff(j) + neq
        END IF ! c2
      END DO ! ieq
    END DO ! irs
  END DO ! c1

! ------------------------------------------------------------------------------
! Preallocate the matrix
! ------------------------------------------------------------------------------

  IF ( global%nProcAlloc == 1 ) THEN
    CALL matseqaijsetpreallocation(pregion%preA,0,nnzon,ierr)
  ELSE
    CALL matmpiaijsetpreallocation(pregion%preA,0,nnzon,0,nnzoff,ierr)
  END IF ! global%nProcAlloc

! ------------------------------------------------------------------------------
! Deallocate preallocation variables
! ------------------------------------------------------------------------------

  DEALLOCATE(nnzon,stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_deallocate,__line__,'nnzon')
  END IF !global%error

  DEALLOCATE(nnzoff,stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_deallocate,__line__,'nnzoff')
  END IF !global%error

! ==============================================================================
! Pre-Fill the Jacobian Matrix
! ==============================================================================

! ------------------------------------------------------------------------------
! Allocate pre-fill variables
! ------------------------------------------------------------------------------

  ALLOCATE(vals(1:neq,1:neq),stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_allocate,__line__,'vals')
  END IF ! global%errorFlag

  ALLOCATE(rows(1:neq),stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_allocate,__line__,'rows')
  END IF ! global%errorFlag

  ALLOCATE(cols(1:neq),stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_allocate,__line__,'cols')
  END IF ! global%errorFlag

! ------------------------------------------------------------------------------
! Create the array of nonzero values
! ------------------------------------------------------------------------------

  DO i = 1, neq
    DO j = 1, neq
      vals(i,j) = 1.0
    END DO ! j
  END DO ! i

! ------------------------------------------------------------------------------
! Set nonzeros in off-diagonal blocks
! ------------------------------------------------------------------------------

  DO c1 = 1, pregion%grid%nCells
! TEMPORARY
!    IF ( pRegion%mixtInput%spaceOrder == DISCR_ORDER_1 ) THEN
      CALL rflu_getresidualsupport1(pregion,c1,rs,rssizemax,rssize)
!    ELSE
!      CALL RFLU_GetResidualSupport2(pRegion,c1,rs,rsSizeMax,rsSize)
!    END IF ! pRegion%mixtInput%spaceOrder
! END TEMPORARY
    DO irs = 1, rssize
      c2 = rs(irs)
      DO ieq = cv_mixt_dens, cv_mixt_ener
        IF ( global%nProcAlloc == 1 ) THEN
          rows(ieq) = neq*(c1-1) + ieq - 1
          cols(ieq) = neq*(c2-1) + ieq - 1
        ELSE
          rows(ieq) = neq*(pregion%grid%pc2sc(c1)-1) + ieq - 1
          cols(ieq) = neq*(pregion%grid%pc2sc(c2)-1) + ieq - 1
        END IF ! global%nProcAlloc
      END DO ! ieq
      IF ( global%nProcAlloc > 1 ) THEN
        CALL aoapplicationtopetsc(pregion%ao,neq,rows,ierr)
        CALL aoapplicationtopetsc(pregion%ao,neq,cols,ierr)
      END IF ! global%nProcAlloc
      CALL matsetvalues(pregion%preA,neq,rows,neq,cols,vals,insert_values,ierr)
      IF ( c2 <= pregion%grid%nCells ) THEN
        CALL matsetvalues(pregion%preA,neq,cols,neq,rows,vals,insert_values,ierr)
      END IF ! c2
    END DO ! irs
  END DO ! ifg

! ------------------------------------------------------------------------------
! Set nonzeros in on-diagonal blocks
! ------------------------------------------------------------------------------

  DO icg = 1, pregion%grid%nCells
    DO ieq = cv_mixt_dens, cv_mixt_ener
      IF ( global%nProcAlloc == 1 ) THEN
        rows(ieq) = neq*(icg-1) + ieq - 1
      ELSE
        rows(ieq) = neq*(pregion%grid%pc2sc(icg)-1) + ieq - 1
      END IF ! global%nProcAlloc
    END DO ! ieq
    IF ( global%nProcAlloc > 1 ) THEN
      CALL aoapplicationtopetsc(pregion%ao,neq,rows,ierr)
    END IF ! global%nProcAlloc
    CALL matsetvalues(pregion%preA,neq,rows,neq,rows,vals,insert_values,ierr)
  END DO ! icg  

! ------------------------------------------------------------------------------
! Deallocate pre-fill variables
! ------------------------------------------------------------------------------
  
  DEALLOCATE(vals,stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_deallocate,__line__,'vals')
  END IF !global%error

  DEALLOCATE(rows,stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_deallocate,__line__,'rows')
  END IF !global%error

  DEALLOCATE(cols,stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_deallocate,__line__,'cols')
  END IF !global%error

! ------------------------------------------------------------------------------
! Assemble the matrix
! ------------------------------------------------------------------------------

  CALL matassemblybegin(pregion%preA,mat_final_assembly,ierr)
  CALL matassemblyend(pregion%preA,mat_final_assembly,ierr)

! DEBUG
!  print*,'writing jacobian to file...'
!  if(global%myProcid == 0) then
!    OPEN(unit=123,file='jacobian_computational_0.m')
!  else
!    OPEN(unit=123,file='jacobian_computational_1.m')
!  endif
!  CALL VecGetSize(pRegion%x,gndof,ierr)
!  CALL VecGetOwnershipRange(pRegion%x,low,high,ierr)
!  DO i = low, high-1
!    DO k = 1, gndof
!      CALL MatGetValues(pRegion%preA,1,i,1,k-1,val,ierr)
!      if(val /= 0.0) then
!         write(123,*) 'A(',i+1,',',k,') = ',val,';'
!      endif
!    ENDDO
!  ENDDO
!  close(123)
!  print*,'done writing jacobian to file'
!  stop
! END DEBUG

! ==============================================================================
! Create Coloring for Jacobian Evaluation
! ==============================================================================

! ------------------------------------------------------------------------------
! Allocate coloring variables
! ------------------------------------------------------------------------------

  CALL rflu_col_createcoloring(pregion)

  ALLOCATE(colors(1:ndof),stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_allocate,__line__,'colors')
  END IF ! global%errorFlag

! ------------------------------------------------------------------------------
! Read the cell coloring from the .col files
! ------------------------------------------------------------------------------

  CALL rflu_col_readcoloring(pregion)
  
! ------------------------------------------------------------------------------
! Compute the DOF coloring
! ------------------------------------------------------------------------------

  DO icg = 1, pregion%grid%nCells
    DO ieq = cv_mixt_dens, cv_mixt_ener
      colors(neq * (icg-1) + ieq) = neq * (pregion%grid%col(icg) - 1) + ieq - 1
    END DO ! ieq
  END DO ! icg

! ------------------------------------------------------------------------------
! Give PETSc the coloring
! ------------------------------------------------------------------------------

  CALL iscoloringcreate(petsc_comm_world,ndof,colors,iscolor,ierr)
  CALL matfdcoloringcreate(pregion%preA,iscolor,pregion%fdcolor,ierr)
  CALL matfdcoloringsetfromoptions(pregion%fdcolor,ierr)
  CALL iscoloringdestroy(iscolor,ierr)

! ------------------------------------------------------------------------------
! Deallocate the coloring variables
! ------------------------------------------------------------------------------

  CALL rflu_col_destroycoloring(pregion)

  DEALLOCATE(colors,stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_deallocate,__line__,'colors')
  END IF !global%error

! ==============================================================================
! Destroy the Connectivity Array
! ==============================================================================

  DEALLOCATE(rs,stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_deallocate,__line__,'rs')
  END IF !global%error

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_petsc_createcoloring




! ******************************************************************************
!
! Purpose: Creates the PETSc Application Ordering context to map Rocflu global
!          dof indices to PETSc global dof indices.
!
! Description: None.
!
! Input:
!   pRegion             data for grid region 
!
! Output: 
!   none
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_petsc_createapplicationordering(pRegion)
  
  USE moddatastruct, ONLY: t_region 

  IMPLICIT NONE

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

#include "finclude/petsc.h"
#include "finclude/petscvec.h"
#include "finclude/petscao.h"

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

  TYPE(t_region), POINTER :: pregion
  TYPE(t_global), POINTER :: global
  INTEGER,ALLOCATABLE,DIMENSION(:) :: numpetsc, numapp
  INTEGER :: icg, ieq, neq, ipetsc, iapp, low, counter, ierr, errorflag

! ******************************************************************************
! Start
! ******************************************************************************

  global => pregion%global
  CALL registerfunction(global,'RFLU_PETSC_CreateApplicationOrdering',&
  'RFLU_ModPETScNewtonKrylov.F90')
  
! ==============================================================================
! Allocate the numbering arrays
! ============================================================================== 

  neq = cv_mixt_ener - cv_mixt_dens + 1  

  ALLOCATE(numpetsc(1:neq*pregion%grid%nCells),stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_allocate,__line__,'numpetsc')
  END IF ! global%errorFlag

  ALLOCATE(numapp(1:neq*pregion%grid%nCells),stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_allocate,__line__,'numapp')
  END IF ! global%errorFlag

! ==============================================================================
! Fill the application and PETSc numbering arrays
! ==============================================================================

  CALL vecgetownershiprange(pregion%x,low,petsc_null_integer,ierr)
  ipetsc = low - 1
  counter = 0
  DO icg = 1, pregion%grid%nCells
    DO ieq = cv_mixt_dens, cv_mixt_ener
      ipetsc = ipetsc + 1
      counter = counter + 1
      iapp = neq*(pregion%grid%pc2sc(icg) - 1) + ieq - 1
      numapp(counter) = iapp
      numpetsc(counter) = ipetsc
    END DO ! ieq
  END DO ! icg

! ==============================================================================
! Create the Application Ordering context
! ==============================================================================

  CALL aocreatebasic(petsc_comm_world, neq*pregion%grid%nCells, &
       numapp, numpetsc, pregion%ao, ierr)

! ==============================================================================
! Deallocate the numbering arrays
! ============================================================================== 

  DEALLOCATE(numpetsc,stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_deallocate,__line__,'numpetsc')
  END IF !global%error

  DEALLOCATE(numapp,stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_deallocate,__line__,'numapp')
  END IF !global%error

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_petsc_createapplicationordering





! ******************************************************************************
!
! Purpose: Destroys the PETSc Application Ordering context
!
! Description: None.
!
! Input:
!   pRegion             data for grid region 
!
! Output: 
!   none
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_petsc_destroyapplicationordering(pRegion)
  
  USE moddatastruct, ONLY: t_region 

  IMPLICIT NONE

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

#include "finclude/petsc.h"
#include "finclude/petscao.h"

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

  TYPE(t_region), POINTER :: pregion
  TYPE(t_global), POINTER :: global
  INTEGER :: ierr

! ******************************************************************************
! Start
! ******************************************************************************

  global => pregion%global
  CALL registerfunction(global,'RFLU_PETSC_DestroyApplicationOrdering',&
  'RFLU_ModPETScNewtonKrylov.F90')
  
! ==============================================================================
! Destroy the Application Ordering context
! ==============================================================================

  CALL aodestroy(pregion%ao,ierr)

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_petsc_destroyapplicationordering






! ******************************************************************************
!
! Purpose: Provides the dt scaling factor for implicit runs
!
! Description: None.
!
! Input:
!   pRegion             data for grid region 
!
! Output: 
!   none
!
! Notes: None.
!
! ******************************************************************************

SUBROUTINE rflu_petsc_getdtscale(global,dtscale)
  
  USE moddatastruct, ONLY: t_region 

  IMPLICIT NONE

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

  TYPE(t_global), POINTER :: global
  REAL(RFREAL) :: dtscale

! ******************************************************************************
! Start
! ******************************************************************************

  CALL registerfunction(global,'RFLU_PETSC_GetDtScale',&
  'RFLU_ModPETScNewtonKrylov.F90')
  
! ==============================================================================
! Compute the dt scaling
! ==============================================================================

  dtscale = 1.1 ** global%currentIter

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_petsc_getdtscale






! ******************************************************************************
! END
! ******************************************************************************
  
END MODULE rflu_modpetscnewtonkrylov