module precision_m
use, intrinsic :: iso_fortran_env
  implicit none
  integer,parameter :: DP=REAL64 !kind(1.0D0)
  integer,parameter :: SP=REAL32 !kind(1.0E0)
  integer,parameter :: WP=DP
  integer,parameter :: MSL=100 ! MAX_STR_LENGTH
end module precision_m
module SolverUtils_m
  use precision_m
  implicit none
#include "finclude/petscdef.h"
  type   Solver_t
  ! problem/equation name used as options prefix
    character(len=20) :: solverName
  ! Solver Contexts
    SNES :: snes
    KSP  :: ksp
    PC   :: pc
    MPI_Comm :: comm
    SNESLineSearch :: linesearch
  ! Solver Convergence Information
    SNESConvergedReason :: reasonSNES
  ! Solver Convergence Criteria

  ! Operators
    Vec :: fxn ! F(X), F_i(X), X={x_j} i=1 to numFun, j=1 to numX
    Mat :: Jac ! Jac_ij = dF_i(X)/dx_j
    Vec :: rhs
    Vec :: sol
    Vec :: initGuess ! redundant ? check later

    class(*),pointer :: SolverCtx => null()
   ! Tolorences
    real(WP) :: relTolX
    real(WP) :: absTolX
    !Vec :: absTolVecX
    real(WP) :: absTolF
    real(WP) :: relTolF
    integer :: maxit
    integer :: maxFunEval

    integer :: NumDof !> Size of the problem
    integer :: numFun
    integer :: numX
    logical :: flag_ksp_ksp = .false.
    logical :: flag_pc_pc = .false.
    logical :: flag_vec_fxn = .false.
    logical :: flag_vec_rhs = .false.
    logical :: flag_vec_sol = .false.
    logical :: flag_snes_snes = .false.
    logical :: flag_mat_jac = .false.
    logical :: flag_mat_A = .false.
  ! Flags for Procedure Pointers
    logical :: flag_getFunction = .false.
    logical :: flag_getJacobian = .false.
    logical :: flag_FormPostCheck = .false.
    logical :: flag_FormLineSearch = .false.
    logical :: flag_ComputeInitialGuess = .false.
    logical :: flag_SNESConvergenceTest = .true.
  ! Use this procedures only if the corresponding flags are true
    procedure(ComputeInitialGuess_abs), pointer :: ComputeInitialGuess=> null()
    procedure(getFunction_abs), pointer :: getFunction=> null() ! assign for nonlinear
    procedure(getFunctionWithDM_abs), pointer :: getFunctionWithDM=> null() ! assign for nonlinear
    procedure(getJacobian_abs), pointer :: getJacobian => null()
    procedure(getJacobianWithDM_abs), pointer :: getJacobianWithDM => null()
    procedure(FormPostCheck_abs), pointer :: FormPostCheck => null()
    procedure(FormLineSearch_abs), pointer :: FormLineSearch => null()
    !procedure(SNESConvergenceTest_abs), pointer :: SNESConvergenceTest => null()
  ! These should always be present
  !  procedure(Solve_abs),pointer:: solve => null()
  contains
    procedure :: initwithoutDM => initSolverwithoutDM_Solver
    procedure :: initwithDM => initSolverwithDM_Solver
    generic :: init =>initwithoutDM,initwithDM
    procedure :: solve => solve_solver
    procedure  :: kill => killSolver
    procedure :: setTolerances => setTolerances_Solver
   end type Solver_t

 interface
   subroutine getFunction_abs(this,x,f,ierr)
   import Solver_t
#include "finclude/petscdef.h"
   class(solver_t),intent(inout):: this
   Vec,intent(inout)     :: x
   Vec,intent(inout)  :: f
   PetscErrorCode,intent(inout) :: ierr
   end subroutine getFunction_abs

   subroutine getFunctionWithDM_abs(this,dm,x,f,ierr)
   import Solver_t
#include "finclude/petscdef.h"
   class(solver_t),intent(inout):: this
   DM :: dm
   Vec,intent(inout)     :: x
   Vec,intent(inout)  :: f
   PetscErrorCode,intent(inout) :: ierr
   end subroutine getFunctionWithDM_abs

   subroutine getJacobian_abs(this,x,Jac,pcMat,flag,ierr)
   import Solver_t
#include "finclude/petscdef.h"
   class(Solver_t), intent(inout) :: this
   Vec,intent(inout)   :: x
   Mat,intent(inout) :: Jac,pcMat
   MatStructure,intent(inout) :: flag
   PetscErrorCode,intent(inout)  :: ierr
   end subroutine getJacobian_abs
   subroutine getJacobianWithDM_abs(this,dm,x,Jac,pcMat,flag,ierr)
   import Solver_t
#include "finclude/petscdef.h"
   class(Solver_t), intent(inout) :: this
   DM :: dm
   Vec,intent(inout)   :: x
   Mat,intent(inout) :: Jac,pcMat
   MatStructure,intent(inout) :: flag
   PetscErrorCode,intent(inout)  :: ierr
   end subroutine getJacobianWithDM_abs

  subroutine FormPostCheck_abs(this,tX,tY,tW,Changed_Y,Changed_W,ierr)
   import Solver_t
#include "finclude/petscdef.h"
   class(Solver_t), intent(inout) :: this
    Vec  :: tX,tY,tW
    PetscBool :: changed_Y,changed_W
    PetscErrorCode :: ierr
  end subroutine FormPostCheck_abs

  subroutine SNESConvergenceTest_abs(this,it,fnorm,lambda,X,dX,reason,ierr)
    import Solver_t
#include "finclude/petscdef.h"
    class(Solver_t), intent(inout) :: this
    integer, intent(in) :: it
    PetscReal, intent(in) :: fnorm
    PetscReal, intent(in) :: lambda
    Vec :: X
    Vec :: dX
    SNESConvergedReason :: reason
    PetscErrorCode :: ierr
  end subroutine SNESConvergenceTest_abs
  subroutine FormLineSearch_abs(this,tx,tf,tg,ty,tw,tfnorm,tynorm,tgnorm,tflag,ierr)
    import Solver_t
#include "finclude/petscdef.h"
    class(Solver_t) :: this
    Vec :: tx
    Vec :: tf
    Vec :: tg
    Vec :: ty
    Vec :: tw
    PetscReal :: tfnorm
    PetscReal :: tynorm
    PetscReal :: tgnorm
    PetscBool :: tflag
    PetscErrorCode :: ierr
  end subroutine FormLineSearch_abs
  subroutine ComputeInitialGuess_abs(this,tx,ierr)
    import Solver_t
#include "finclude/petsc.h90"
    class(Solver_t) :: this
    Vec :: tx
    PetscErrorCode :: ierr
  end subroutine ComputeInitialGuess_abs

   subroutine initSolver_abs(this,comm,SolverType,numEqation,numUnknown,ierr,EqnCtx)
   import Solver_t
#include "finclude/petscdef.h"
   class(Solver_t), intent(inout)  :: this
   MPI_Comm :: comm
   integer, intent(in)           :: solverType
   integer, intent(in)           :: numEqation
   integer, intent(in)           :: numUnknown
   PetscErrorCode, intent(inout) :: ierr
   class(*),pointer,intent(inout),optional :: EqnCtx
   end subroutine initSolver_abs

   subroutine Solve_abs(this,ierr,opt_solverMode)
   import Solver_t
#include "finclude/petscdef.h"
   class(Solver_t), intent(inout)  :: this
   PetscErrorCode, intent(inout) :: ierr
   integer, optional :: opt_solverMode
   end subroutine Solve_abs
end interface

  contains

  subroutine killSolver(this,ierr)
   implicit none
#include "finclude/petsc.h"

  class(Solver_t) :: this
  PetscErrorCode ierr
   if(this%flag_vec_rhs) call VECDestroy(this%rhs,ierr)
   if(this%flag_vec_fxn) call VECDestroy(this%fxn,ierr)
   if(this%flag_vec_sol) call VECDestroy(this%sol,ierr)

   if(this%flag_mat_jac) call MatDestroy(this%Jac,ierr)

   if(this%flag_pc_pc)   call PCDestroy(this%pc,ierr)
   if(this%flag_ksp_ksp) call KSPDestroy(this%ksp,ierr)
   if(this%flag_snes_snes) call SNESDestroy(this%snes,ierr)

  end subroutine killSolver



subroutine FormJacobian_snes(tsnes,tX,tjac,tpreconJac,tflag,this,tierr)
 !use SolverUtils
 implicit none
#include "finclude/petsc.h90"
      SNES     tsnes
      Vec      tx
      Mat          tjac,tpreconJac
      MatStructure tflag
      type(Solver_t), intent(inout) :: this
      PetscErrorCode tierr
    call this%getJacobian(tx,tjac,tpreconJac,tflag,tierr)
end subroutine FormJacobian_snes

  subroutine FunctionJacobianWithDM(dm,X,jac,preconJac,flag,this,ierr)
 !use SolverUtils
    implicit none
#include "finclude/petsc.h90"
      DM     dm
      Vec    X
      Mat    Jac,preconJac
      MatStructure flag
      type(Solver_t), intent(inout) :: this
      PetscErrorCode ierr
     call this%getJacobianWithDM(dm,X,Jac,preconJac,flag,ierr)
  end subroutine FunctionJacobianWithDM

  subroutine FormFunction_snes(tsnes,tx,tf,this,tierr)
  !use SolverUtils
   implicit none
#include "finclude/petsc.h90"
    SNES  :: tsnes
    Vec   :: tx,tf
    type(Solver_t), intent(inout) :: this
    PetscErrorCode tierr
    !  print*,'Calling Get function'
    call this%getFunction(tx,tf,tierr)
    !  print*,'got function'
  end subroutine FormFunction_snes

  subroutine FunctionResidualWithDM(dm,X,F,this,ierr)
   implicit none
#include "finclude/petsc.h90"
   DM :: dm
   Vec :: X
   Vec :: F
   type(Solver_t), intent(inout) :: this
   PetscErrorCode :: ierr
   call this%getFunctionWithDM(dm,X,F,ierr)
  end subroutine FunctionResidualWithDM



  subroutine FormPostCheck_snes(tlinesearch,tX,tY,tW,changed_Y,changed_W,this,ierr)
    implicit none
#include "finclude/petsc.h90"
    SNESLineSearch :: tlinesearch
    Vec  :: tX,tY,tW
    PetscBool :: changed_Y,changed_W
    type(Solver_t) :: this
    PetscErrorCode :: ierr
      call this%FormPostCheck(tX,tY,tW,changed_Y,changed_W,ierr)
  end subroutine FormPostCheck_snes

  subroutine FormLineSearch_snes(tsnes,this,tx,tf,tg,ty,tw,tfnorm,tynorm,tgnorm,tflag,ierr)
    implicit none
#include "finclude/petsc.h90"
    SNES :: tsnes
    type(Solver_t) :: this
    Vec :: tx
    Vec :: tf
    Vec :: tg
    Vec :: ty
    Vec :: tw
    PetscReal :: tfnorm
    PetscReal :: tynorm
    PetscReal :: tgnorm
    PetscBool :: tflag
    PetscErrorCode :: ierr
    call this%FormLineSearch(tx,tf,tg,ty,tw,tfnorm,tynorm,tgnorm,tflag,ierr)
  end subroutine FormLineSearch_snes

subroutine SNESConvergenceTest_snes(snes,it,xnorm,gnorm,fnorm,reason,this, ierr)
  implicit none
#include "finclude/petsc.h90"
  ! IO Variables
  SNES :: snes
  integer :: it
  PetscReal :: xnorm
  PetscReal :: gnorm
  PetscReal :: fnorm
  SNESConvergedReason :: reason
  type(Solver_t) :: this
  PetscErrorCode :: ierr
  ! Local Variables
  Vec :: X, dX
  Vec :: F, W, G
  real(WP) :: lambda
  print*, 'Calling custom test',it,xnorm,fnorm,gnorm
  !call SNESLineSearchGetLambda(this%linesearch, lambda, ierr)
  !call SNESLineSearchGetVecs(this%linesearch, X, F, dX, W, G, ierr)
  !call this%SNESConvergenceTest(it, fnorm, lambda, X, dX, reason, ierr)
  reason = SNES_DIVERGED_FUNCTION_COUNT
end subroutine SNESConvergenceTest_snes

! Do not need this, same functionality can be achived by
! SNESKSPSetUseEW(SNES snes,PetscBool flag) flags ==.true.
!subroutine KSPConvergenceTest_snes(ksp, it,rNorm, reason, this, ierr)
!  implicit none
!#include "finclude/petsc.h"
!  ! IO Variables
!  KSP :: ksp
!  integer :: it
!  real(WP) :: rnorm
!  KSPConvergedReason :: reason
!  type(Solver_t) :: this
!  PetscErrorCode :: ierr
!  ! Local Variables
!  real(WP) :: absTol
!  real(WP) :: normF
!  absTol = minval([normF,1E-9*normF**(1.0_WP+0.5_WP)])
!  absTol = maxVal([absTol, this%absTolF])
!  call KSPSetTolerances(ksp,PETSC_DEFAULT_DOUBLE_PRECISION, absTol,  &
!                        PETSC_DEFAULT_DOUBLE_PRECISION,PETSC_DEFAULT_INTEGER,ierr)
!  call KSPDefaultConverged(ksp,it,rnorm,reason,ierr)
!end subroutine KSPConvergenceTest_snes

  subroutine ComputeInitialGuess_snes(tsnes,tx,this,ierr)
    implicit none
#include "finclude/petsc.h"
    SNES :: tsnes
    Vec :: tx
    type(Solver_t) :: this
    PetscErrorCode :: ierr
    call this%ComputeInitialGuess(tx,ierr)
  end subroutine ComputeInitialGuess_snes

  subroutine Solve_solver(this,ierr)
    !use SolverUtils
    implicit none
#include "finclude/petsc.h"
    class(Solver_t), intent(inout)  :: this
    PetscErrorCode, intent(inout) :: ierr
    integer :: solverMode
    PetscInt :: maxPos
    PetscReal :: maxRhs

    !    print*,'entered Nonlinear solve'
        call SNESSolve(this%snes,PETSC_NULL_OBJECT,this%sol,ierr); CHKERRQ(ierr)
        call SNESGetConvergedReason(this%snes,this%reasonSNES,ierr); CHKERRQ(ierr)
        if(this%reasonSNES>0) then
    !      print*,'SNESConvergedReason ', this%reasonSNES
        else
    !      print*,'SNESDivergedReason ', this%reasonSNES
        endif
    !print*,'Exiting solve'
  end subroutine Solve_solver

  subroutine initSolverWithoutDM_Solver(this,comm,solverName,NumEqn,numVar,ierr, &
                               EqnCtx, getFunction, ComputeInitialGuess, &
                               getJacobian, FormLineSearch, &
                               FormPostCheck)
    implicit none
#include "finclude/petsc.h"
    class(Solver_t), intent(inout) :: this
    MPI_Comm            :: comm
    character(len=*), intent(in) :: solverName
    integer, intent(in) :: NumEqn
    integer, intent(in) :: NumVar
    PetscErrorCode :: ierr
    class(*),pointer,intent(inout),optional :: EqnCtx
    procedure(getFunction_abs), optional :: getFunction! assign for nonlinear
    procedure(ComputeInitialGuess_abs), optional :: ComputeInitialGuess
    procedure(getJacobian_abs), optional :: getJacobian
    procedure(FormLineSearch_abs), optional :: FormLineSearch
    procedure(FormPostCheck_abs), optional :: FormPostCheck
    character(len=21) :: optionsPrefix
    PetscReal :: rtol ! relative tolorance
    PetscReal :: abstol ! absolute tolrance
    PetscReal :: dtol
    PetscReal :: solTol
    PetscInt  :: maxits
    PetscInt :: maxfunEval
    ! Local Variables
    PetscInt :: psize ! problem size
    PetscScalar :: pfive !
    PetscInt :: nz
    ! Set solver type

    this%comm = comm

    !if(length(trim(solverName)) > 20 ) then
    this%solverName = trim(solverName) ! this will truncate solverName to 20 characters
    optionsPrefix = trim(this%solverName)//'_'
    ! Set problem size
    this%NumDof = numVar
    this%numFun = numEqn
    this%numX = numVar
    psize = this%numDof
    !

    if(present(EqnCtx)) then
      this%solverCtx => EqnCtx
    end if

    if(present(getFunction)) then
      this%getFunction => getFunction
      this%flag_getFunction = .true.
    else
      print*,'Error :: Must provide get function for a Nonlinear solver'
      stop
    end if

    if(present(ComputeInitialGuess)) then
      this%ComputeInitialGuess => ComputeInitialGuess
      this%flag_ComputeInitialGuess = .true.
    end if
    ! Link the Jacobian Evaluation routine
    ! Genrally optional but here it is must
    if(present(getJacobian)) then
      this%getJacobian => getJacobian
      this%flag_getJacobian = .true.
    end if
    ! Link the FormLineSearch
    ! Optional
    if(present(FormLineSearch)) then
      this%FormLineSearch => FormLineSearch
      this%flag_FormLineSearch = .true.
    end if
    ! Link the FormLineSearch
    ! Optional
    if(present(FormPostCheck)) then
      this%FormPostCheck => FormPostCheck
      this%flag_FormPostCheck = .true.
    end if

        ! Create solution and fxn vectors
    call VecCreateMPI(this%comm,PETSC_DECIDE,this%numx,this%sol,ierr)
    call VecSetOption(this%sol,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE,ierr)
    call VecCreateMPI(this%comm,PETSC_DECIDE,this%numFun,this%fxn,ierr)

    this%flag_vec_sol = .true.
    this%flag_vec_fxn = .true.

    call MatCreate(this%comm,this%Jac,ierr)
    call MatSetSizes(this%Jac,PETSC_DECIDE,PETSC_DECIDE,this%numFun,this%numX,ierr)
    call MatSetType(this%Jac,MATSEQAIJ,ierr)
    !call MatSetType(this%Jac,MATMPIAIJ,ierr)
    call MatSetFromOptions(this%Jac,ierr)
    call MatSetUp(this%Jac,ierr)
    this%flag_mat_Jac = .true.

    !  create snes context
    call SNESCreate(this%comm,this%snes,ierr)
    this%flag_snes_snes = .true.
    !  Set function evaluation routine and vector
    ! use the defualt FormFunction_snes
    if(this%flag_getFunction) then
      call SNESSetFunction(this%snes,this%fxn,FormFunction_snes,this,ierr)
    end if
    !  Set Jacobian matrix data structure and Jacobian evaluation routine
    if(this%flag_getJacobian) then
      call SNESSetJacobian(this%snes,this%Jac,this%Jac,FormJacobian_snes,this,ierr)
    else
    ! use finite difference jacobian
    end if
    !  Set SNES Line Search Method
    if(this%flag_FormLineSearch) then
    !  call SNESLineSearchSet(this%snes,FormLineSearch_snes,this,ierr)
    end if
    ! Set SNES Line Search Post Check
    if(this%flag_FormPostCheck) then
    !  call SNESLineSearchSetPostCheck(this%snes,FormPostCheck_snes,this,ierr)
    end if
    ! Intial Guess
    if(this%flag_ComputeInitialGuess) then
      !call SNESSetComputeInitialGuess(this%snes,ComputeInitialGuess_snes,this,ierr)
      pfive = 0.0_WP
      call VecSet(this%sol,pfive,ierr)
    else
      pfive = 0.0_WP
      call VecSet(this%sol,pfive,ierr)
    end if

   call SNESSetConvergenceTest(this%snes, SNESConvergenceTest_snes, this,PETSC_NULL_FUNCTION, ierr); CHKERRQ(ierr)


    rtol = 1.0E-15!*25.6E-3!PETSC_DEFAULT_DOUBLE_PRECISION
    abstol = 1.0E-6!*25.6E-3!PETSC_DEFAULT_DOUBLE_PRECISION
    soltol = PETSC_DEFAULT_DOUBLE_PRECISION
    maxIts = PETSC_DEFAULT_INTEGER
    maxfunEval = PETSC_DEFAULT_INTEGER
    call SNESSetTolerances(this%snes,abstol,rtol,soltol,maxits,maxfunEval,ierr)
    call SNESSetOptionsPrefix(this%snes,trim(optionsPrefix),ierr)
    call SNESSetFromOptions(this%snes,ierr)
    print*,'Sucessfully intiated Nonlinear Solver for ',trim(solverName)


  end subroutine initSolverWithoutDM_Solver

 ! The DM should already come with field section

  subroutine initSolverWithDM_Solver(this,comm,solverName,meshDM,ierr, &
                               EqnCtx, getFunction, ComputeInitialGuess, &
                               getJacobian, FormLineSearch, &
                               FormPostCheck,SNESConvergenceTest)
    implicit none
#include "finclude/petsc.h90"
    ! IO Variables
    class(Solver_t), intent(inout) :: this
    MPI_Comm            :: comm
    character(len=*), intent(in) :: solverName
    DM :: meshDM
    PetscErrorCode :: ierr
    class(*),pointer,intent(inout),optional :: EqnCtx
    procedure(getFunctionWithDM_abs), optional :: getFunction! assign for nonlinear
    procedure(ComputeInitialGuess_abs), optional :: ComputeInitialGuess
    procedure(getJacobianWithDM_abs), optional :: getJacobian
    procedure(FormLineSearch_abs), optional :: FormLineSearch
    procedure(FormPostCheck_abs), optional :: FormPostCheck
    procedure(SNESConvergenceTest_abs), optional :: SNESConvergenceTest
    ! Local Variables
    character(len=21) :: optionsPrefix
    KSP :: ksp
    PC :: pc
    PetscReal :: reltolF ! relative tolorance
    PetscReal :: abstolF ! absolute tolrance
    PetscReal :: reltolX ! relative tolorance
    PetscReal :: abstolX ! absolute tolrance
    PetscReal :: dtol
    PetscReal :: solTol
    PetscInt  :: maxit
    PetscInt :: maxfunEval
    PetscScalar :: pfive

    this%comm = comm
    !if(length(trim(solverName)) > 20 ) then
    this%solverName = trim(solverName) ! this will truncate solverName to 20 characters
    optionsPrefix = trim(this%solverName)//'_'

    if(present(EqnCtx)) then
      this%solverCtx => EqnCtx
    end if

    if(present(getFunction)) then
      this%getFunctionWithDM => getFunction
      this%flag_getFunction = .true.
    else
      print*,'Error :: Must provide get function for a Nonlinear solver'
      stop
    end if

    if(present(ComputeInitialGuess)) then
      this%ComputeInitialGuess => ComputeInitialGuess
      this%flag_ComputeInitialGuess = .true.
    end if
    ! Link the Jacobian Evaluation routine
    ! Genrally optional but here it is must
    if(present(getJacobian)) then
      this%getJacobianWithDM => getJacobian
      this%flag_getJacobian = .true.
    end if
    ! Link the FormLineSearch
    ! Optional
    if(present(FormLineSearch)) then
      this%FormLineSearch => FormLineSearch
      this%flag_FormLineSearch = .true.
    end if
    ! Link the FormLineSearch
    ! Optional
    if(present(FormPostCheck)) then
      this%FormPostCheck => FormPostCheck
      this%flag_FormPostCheck = .true.
    end if
    ! Link the SNESConvergenceTest
    if (present(SNESConvergenceTest)) then
      !this%SNESConvergenceTest => SNESConvergenceTest
      this%flag_SNESConvergenceTest = .true.
    end if
    !  create snes context
    call SNESCreate(this%comm,this%snes,ierr)
    this%flag_snes_snes = .true.
    ! Set SNES DM
    call SNESSetDM(this%snes,meshDM,ierr)
    ! Create solution and fxn vectors
    call DMCreateGlobalVector(meshDM,this%sol,ierr)
    call VecSetOption(this%sol,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE,ierr)
    call VecDuplicate(this%sol,this%fxn,ierr)
    this%flag_vec_sol = .true.
    this%flag_vec_fxn = .true.
    ! Create Jacobian Matrix
    call DMSetMatType(meshDM, MATAIJ, ierr)
    call DMCreateMatrix(meshDM, this%Jac,ierr)
    this%flag_mat_Jac = .true.

    !  Set function evaluation routine and vector
    ! use the defualt FormFunction_snes
    if(this%flag_getFunction) then
      call DMSNESSetFunctionLocal(meshDM,FunctionResidualWithDM,this,ierr)
    end if
    !  Set Jacobian matrix data structure and Jacobian evaluation routine
    if(this%flag_getJacobian) then
      call DMSNESSetJacobianLocal(meshDM,FunctionJacobianWithDM,this,ierr)
    else
    ! use finite difference jacobian
    end if
    !  Set SNES Line Search Method
    if(this%flag_FormLineSearch) then
    !  call SNESLineSearchSet(this%snes,FormLineSearch_snes,this,ierr)
    end if
    ! Set SNES Line Search Post Check
    if(this%flag_FormPostCheck) then
      call SNESGetLineSearch(this%snes, this%linesearch, ierr)
      call SNESLineSearchSetPostCheck(this%linesearch,FormPostCheck_snes,this,ierr)
    end if
    ! Set SNES Convergence Test
    !if(this%flag_SNESConvergenceTest) then
      !call SNESGetLineSearch(this%snes, this%linesearch, ierr)
      call SNESSetConvergenceTest(this%snes, SNESConvergenceTest_snes, this,PETSC_NULL_FUNCTION, ierr); CHKERRQ(ierr)
    !end if
    ! Intial Guess
    if(this%flag_ComputeInitialGuess) then
      !call SNESSetComputeInitialGuess(this%snes,ComputeInitialGuess_snes,this,ierr)
      pfive = 0.5
      call VecSet(this%sol,pfive,ierr)
    else
      pfive = 0.5
      call VecSet(this%sol,pfive,ierr)
    end if

    call SNESGetKSP(this%snes, ksp, ierr)
    call KSPSetFromOptions(ksp, ierr)

    call KSPGetPC(ksp, pc, ierr)
    call PCSetFromOptions(pc,ierr)

    reltolF = 1.0E-10!*25.6E-3!PETSC_DEFAULT_DOUBLE_PRECISION
    abstolF = 1.0E-6!*25.6E-3!PETSC_DEFAULT_DOUBLE_PRECISION
    abstolX = 1.0E-3
    relTolX = 1.0E-5
    maxIt = PETSC_DEFAULT_INTEGER
    maxfunEval = PETSC_DEFAULT_INTEGER
    !call SNESSetTolerances(this%snes,abstol,rtol,soltol,maxits,maxfunEval,ierr)
    call this%setTolerances(absTolF,relTolF,absTolX, relTolF, maxIt,maxFunEval)
    !call SNESKSPSetUseEW(this%snes, PETSC_TRUE, ierr)
    call SNESSetOptionsPrefix(this%snes,trim(optionsPrefix),ierr)
    call SNESSetFromOptions(this%snes,ierr)
    print*,'Sucessfully intiated Nonlinear Solver for ',trim(solverName),' model'

  end subroutine initSolverWithDM_Solver

  subroutine setTolerances_Solver(this,absTolF,relTolF,absTolX,relTolX,maxit,maxFunEval)
    implicit none
#include "finclude/petsc.h"
    ! IO Variables
    class(Solver_t), intent(inout) :: this
    real(WP), intent(in) :: absTolF
    real(WP), intent(in) :: relTolF
    real(WP), intent(in) :: absTolX
    real(WP), intent(in) :: relTolX
    integer, intent(in) :: maxit
    integer, intent(in) :: maxFunEval
    ! Local Variables
    PetscErrorCode :: ierr
    this%absTolF = absTolF
    this%relTolF = relTolF
    this%relTolX = relTolX
    this%absTolX = absTolX
    this%maxit = maxit
    this%maxFunEval = maxFunEVal
    call SNESSetTolerances(this%snes,absTolF,relTolF,relTolX,maxit,maxfunEval,ierr)
  end subroutine setTolerances_Solver
end module SolverUtils_m

#if TestSolver
module SolverTest_m
 use SolverUtils_m
 implicit none


 contains

 subroutine SolverTest()
 use SolverUtils_m
 implicit none
#include "finclude/petsc.h90"

  type(Solver_t),pointer :: test => null()
  class(*),pointer :: solverCtx
  PetscInt :: i2=2
  PetscScalar :: pfive
  PetscErrorCode :: ierr
  PetscScalar, pointer  :: xxv(:)
  allocate(test)
  !test%init => initSolver
  !test%solve => Solve
  test%solverCtx => test
  !test%getFunction =>getFxnTest
  !test%getJacobian => getJacTest
  call test%init(PETSC_COMM_SELF,'test',i2,i2,ierr,solverCtx,&
                 getFunction=getFxnTest,getJacobian=getJacTest)
  call test%solve(ierr)
  call VecGetArrayF90(test%sol,xxv,ierr)
  print*,xxv
  call VecRestoreArrayF90(test%sol,xxv,ierr)
  print*,'Sucessfully Tested solver Module'
 end subroutine SolverTest

  subroutine getFxnTest(this,x,f,ierr)
  use SolverUtils_m
   implicit none
#include "finclude/petsc.h90"

   class(Solver_t), intent(inout):: this
   Vec  :: x
   Vec  :: f
   PetscErrorCode  :: ierr
   PetscScalar, pointer  :: xxv(:)
   PetscScalar, pointer :: ffv(:)

   ! print*,'using xxv'
    call VecGetArrayF90(x,xxv,ierr)
    call VecGetArrayF90(f,ffv,ierr)
    !ffv(1) = xxv(1) - 3.0 !xxv(1)*xxv(1) + xxv(1)*xxv(2) - 3.0
    !ffv(2) = xxv(2) - 6.0!xxv(1)*xxv(2)+ xxv(2)*xxv(2) - 6.0
    ffv(1) = xxv(1)*xxv(1) + xxv(1)*xxv(2) - 3.0
    ffv(2) = xxv(1)*xxv(2)+ xxv(2)*xxv(2) - 6.0
  !  print*,'setting ffv',ffv
    call VecRestoreArrayF90(f,ffv,ierr)
    call VecRestoreArrayF90(x,xxv,ierr)

   end subroutine getFxnTest

   subroutine getJacTest(this,x,Jac,pcMat,flag,ierr)
   use SolverUtils_m
   implicit none
#include "finclude/petsc.h90"

   class(Solver_t), intent(inout) :: this
   Vec  :: x
   Mat  :: Jac,pcMAt
   MatStructure  :: flag
   PetscErrorCode  :: ierr
 ! Local Variables
   PetscScalar, pointer   :: xxv(:)
   PetscScalar :: A(4)
   PetscInt :: idx(2),i2
      call VecGetArrayF90(x,xxv,ierr)
      i2=2
      idx(1) = 0
      idx(2) = 1
      !A(1) = 1.0!2.0*xxv(1) + xxv(2)
      !A(2) = 0.0!xxv(1)
      !A(3) = 0.0!xxv(2)
      !A(4) = 1.0!xxv(1) + 2.0*xxv(2)
      A(1) = 2.0*xxv(1) + xxv(2)
      A(2) = xxv(1)
      A(3) = xxv(2)
      A(4) = xxv(1) + 2.0*xxv(2)
      call VecRestoreArrayF90(x,xxv,ierr)
      call MatSetValues(Jac,i2,idx,i2,idx,A,INSERT_VALUES,ierr)
      call MatAssemblyBegin(Jac,MAT_FINAL_ASSEMBLY,ierr)
      call MatAssemblyEnd(Jac,MAT_FINAL_ASSEMBLY,ierr)
      flag = SAME_NONZERO_PATTERN
   end subroutine getJacTest


end module SolverTest_m
program Test
 use SolverTest_m
  implicit none
#include "finclude/petsc.h"

 PetscErrorCode ierr

call PetscInitialize(PETSC_NULL_CHARACTER, ierr)
 print*,'petscIntialiseSucess'
 call SolverTest()

 call PetscFinalize(ierr)
end program Test
#endif