xref: /petsc/src/snes/tests/ex12f.F90 (revision 77d968b72e8e27b79bcc994c018975de390644ed)

!
!
!  This example demonstrates basic use of the SNES Fortran interface.
!
!
        module ex12fmodule
#include <petsc/finclude/petscsnes.h>
        use petscsnes
        type User
          DM  da
          Vec F
          Vec xl
          MPI_Comm comm
          PetscInt N
        end type User
        save
        type monctx
        PetscInt :: its,lag
        end type monctx
      end module

! ---------------------------------------------------------------------
! ---------------------------------------------------------------------
!  Subroutine FormMonitor
!  This function lets up keep track of the SNES progress at each step
!  In this routine, we determine when the Jacobian is rebuilt with the parameter 'jag'
!
!  Input Parameters:
!    snes    - SNES nonlinear solver context
!    its     - current nonlinear iteration, starting from a call of SNESSolve()
!    norm    - 2-norm of current residual (may be approximate)
!    snesm - monctx designed module (included in Snesmmod)
! ---------------------------------------------------------------------
      subroutine FormMonitor(snes,its,norm,snesm,ierr)
      use ex12fmodule
      implicit none

      SNES ::           snes
      PetscInt ::       its,one,mone
      PetscScalar ::    norm
      type(monctx) ::   snesm
      PetscErrorCode :: ierr

!      write(6,*) ' '
!      write(6,*) '    its ',its,snesm%its,'lag',
!     &            snesm%lag
!      call flush(6)
      if (mod(snesm%its,snesm%lag).eq.0) then
        one = 1
        PetscCall(SNESSetLagJacobian(snes,one,ierr))  ! build jacobian
      else
        mone = -1
        PetscCall(SNESSetLagJacobian(snes,mone,ierr)) ! do NOT build jacobian
      endif
      snesm%its = snesm%its + 1
      end subroutine FormMonitor

!  Note: Any user-defined Fortran routines (such as FormJacobian)
!  MUST be declared as external.
!
!
! Macros to make setting/getting  values into vector clearer.
! The element xx(ib) is the ibth element in the vector indicated by ctx%F
#define xx(ib)  vxx(ixx + (ib))
#define ff(ib)  vff(iff + (ib))
#define F2(ib)  vF2(iF2 + (ib))
      program main
      use ex12fmodule
      implicit none
      type(User) ctx
      PetscMPIInt rank,size
      PetscErrorCode ierr
      PetscInt N,start,end,nn,i
      PetscInt ii,its,i1,i0,i3
      PetscBool  flg
      SNES             snes
      Mat              J
      Vec              x,r,u
      PetscScalar      xp,FF,UU,h
      character*(10)   matrixname
      external         FormJacobian,FormFunction
      external         formmonitor
      type(monctx) :: snesm

      PetscCallA(PetscInitialize(ierr))
      i1 = 1
      i0 = 0
      i3 = 3
      N  = 10
      PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',N,flg,ierr))
      h = 1.0/real(N-1)
      ctx%N = N
      ctx%comm = PETSC_COMM_WORLD

      PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
      PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,size,ierr))

! Set up data structures
      PetscCallA(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,N,i1,i1,PETSC_NULL_INTEGER,ctx%da,ierr))
      PetscCallA(DMSetFromOptions(ctx%da,ierr))
      PetscCallA(DMSetUp(ctx%da,ierr))
      PetscCallA(DMCreateGlobalVector(ctx%da,x,ierr))
      PetscCallA(DMCreateLocalVector(ctx%da,ctx%xl,ierr))

      PetscCallA(PetscObjectSetName(x,'Approximate Solution',ierr))
      PetscCallA(VecDuplicate(x,r,ierr))
      PetscCallA(VecDuplicate(x,ctx%F,ierr))
      PetscCallA(VecDuplicate(x,U,ierr))
      PetscCallA(PetscObjectSetName(U,'Exact Solution',ierr))

      PetscCallA(MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N,N,i3,PETSC_NULL_INTEGER,i0,PETSC_NULL_INTEGER,J,ierr))
      PetscCallA(MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE,ierr))
      PetscCallA(MatGetType(J,matrixname,ierr))

! Store right-hand-side of PDE and exact solution
      PetscCallA(VecGetOwnershipRange(x,start,end,ierr))
      xp = h*start
      nn = end - start
      ii = start
      do 10, i=0,nn-1
        FF = 6.0*xp + (xp+1.e-12)**6.e0
        UU = xp*xp*xp
        PetscCallA(VecSetValues(ctx%F,i1,ii,FF,INSERT_VALUES,ierr))
        PetscCallA(VecSetValues(U,i1,ii,UU,INSERT_VALUES,ierr))
        xp = xp + h
        ii = ii + 1
 10   continue
      PetscCallA(VecAssemblyBegin(ctx%F,ierr))
      PetscCallA(VecAssemblyEnd(ctx%F,ierr))
      PetscCallA(VecAssemblyBegin(U,ierr))
      PetscCallA(VecAssemblyEnd(U,ierr))

! Create nonlinear solver
      PetscCallA(SNESCreate(PETSC_COMM_WORLD,snes,ierr))

! Set various routines and options
      PetscCallA(SNESSetFunction(snes,r,FormFunction,ctx,ierr))
      PetscCallA(SNESSetJacobian(snes,J,J,FormJacobian,ctx,ierr))

      snesm%its = 0
      PetscCallA(SNESGetLagJacobian(snes,snesm%lag,ierr))
      PetscCallA(SNESMonitorSet(snes,FormMonitor,snesm,PETSC_NULL_FUNCTION,ierr))
      PetscCallA(SNESSetFromOptions(snes,ierr))

! Solve nonlinear system
      PetscCallA(FormInitialGuess(snes,x,ierr))
      PetscCallA(SNESSolve(snes,PETSC_NULL_VEC,x,ierr))
      PetscCallA(SNESGetIterationNumber(snes,its,ierr))

!  Free work space.  All PETSc objects should be destroyed when they
!  are no longer needed.
      PetscCallA(VecDestroy(x,ierr))
      PetscCallA(VecDestroy(ctx%xl,ierr))
      PetscCallA(VecDestroy(r,ierr))
      PetscCallA(VecDestroy(U,ierr))
      PetscCallA(VecDestroy(ctx%F,ierr))
      PetscCallA(MatDestroy(J,ierr))
      PetscCallA(SNESDestroy(snes,ierr))
      PetscCallA(DMDestroy(ctx%da,ierr))
      PetscCallA(PetscFinalize(ierr))
      end

! --------------------  Evaluate Function F(x) ---------------------

      subroutine FormFunction(snes,x,f,ctx,ierr)
      use ex12fmodule
      implicit none
      SNES             snes
      Vec              x,f
      type(User) ctx
      PetscMPIInt  rank,size,zero
      PetscInt i,s,n
      PetscErrorCode ierr
      PetscOffset      ixx,iff,iF2
      PetscScalar      h,d,vf2(2),vxx(2),vff(2)

      zero = 0
      PetscCallMPI(MPI_Comm_rank(ctx%comm,rank,ierr))
      PetscCallMPI(MPI_Comm_size(ctx%comm,size,ierr))
      h     = 1.0/(real(ctx%N) - 1.0)
      PetscCall(DMGlobalToLocalBegin(ctx%da,x,INSERT_VALUES,ctx%xl,ierr))
      PetscCall(DMGlobalToLocalEnd(ctx%da,x,INSERT_VALUES,ctx%xl,ierr))

      PetscCall(VecGetLocalSize(ctx%xl,n,ierr))
      if (n .gt. 1000) then
        print*, 'Local work array not big enough'
        call MPI_Abort(PETSC_COMM_WORLD,zero,ierr)
      endif

!
! This sets the index ixx so that vxx(ixx+1) is the first local
! element in the vector indicated by ctx%xl.
!
      PetscCall(VecGetArrayRead(ctx%xl,vxx,ixx,ierr))
      PetscCall(VecGetArray(f,vff,iff,ierr))
      PetscCall(VecGetArray(ctx%F,vF2,iF2,ierr))

      d = h*h

!
!  Note that the array vxx() was obtained from a ghosted local vector
!  ctx%xl while the array vff() was obtained from the non-ghosted parallel
!  vector F. This is why there is a need for shift variable s. Since vff()
!  does not have locations for the ghost variables we need to index in it
!  slightly different then indexing into vxx(). For example on processor
!  1 (the second processor)
!
!        xx(1)        xx(2)             xx(3)             .....
!      ^^^^^^^        ^^^^^             ^^^^^
!      ghost value   1st local value   2nd local value
!
!                      ff(1)             ff(2)
!                     ^^^^^^^           ^^^^^^^
!                    1st local value   2nd local value
!
       if (rank .eq. 0) then
        s = 0
        ff(1) = xx(1)
      else
        s = 1
      endif

      do 10 i=1,n-2
       ff(i-s+1) = d*(xx(i) - 2.0*xx(i+1) + xx(i+2)) + xx(i+1)*xx(i+1) - F2(i-s+1)
 10   continue

      if (rank .eq. size-1) then
        ff(n-s) = xx(n) - 1.0
      endif

      PetscCall(VecRestoreArray(f,vff,iff,ierr))
      PetscCall(VecRestoreArrayRead(ctx%xl,vxx,ixx,ierr))
      PetscCall(VecRestoreArray(ctx%F,vF2,iF2,ierr))
      return
      end

! --------------------  Form initial approximation -----------------

      subroutine FormInitialGuess(snes,x,ierr)
      use ex12fmodule
      implicit none

      PetscErrorCode   ierr
      Vec              x
      SNES             snes
      PetscScalar      five

      five = .5
      PetscCall(VecSet(x,five,ierr))
      return
      end

! --------------------  Evaluate Jacobian --------------------

      subroutine FormJacobian(snes,x,jac,B,ctx,ierr)
      use ex12fmodule
      implicit none

      SNES             snes
      Vec              x
      Mat              jac,B
      type(User) ctx
      PetscInt  ii,istart,iend
      PetscInt  i,j,n,end,start,i1
      PetscErrorCode ierr
      PetscMPIInt rank,size
      PetscOffset      ixx
      PetscScalar      d,A,h,vxx(2)

      i1 = 1
      h = 1.0/(real(ctx%N) - 1.0)
      d = h*h
      PetscCallMPI(MPI_Comm_rank(ctx%comm,rank,ierr))
      PetscCallMPI(MPI_Comm_size(ctx%comm,size,ierr))

      PetscCall(VecGetArrayRead(x,vxx,ixx,ierr))
      PetscCall(VecGetOwnershipRange(x,start,end,ierr))
      n = end - start

      if (rank .eq. 0) then
        A = 1.0
        PetscCall(MatSetValues(jac,i1,start,i1,start,A,INSERT_VALUES,ierr))
        istart = 1
      else
        istart = 0
      endif
      if (rank .eq. size-1) then
        i = INT(ctx%N-1)
        A = 1.0
        PetscCall(MatSetValues(jac,i1,i,i1,i,A,INSERT_VALUES,ierr))
        iend = n-1
      else
        iend = n
      endif
      do 10 i=istart,iend-1
        ii = i + start
        j = start + i - 1
        PetscCall(MatSetValues(jac,i1,ii,i1,j,d,INSERT_VALUES,ierr))
        j = start + i + 1
        PetscCall(MatSetValues(jac,i1,ii,i1,j,d,INSERT_VALUES,ierr))
        A = -2.0*d + 2.0*xx(i+1)
        PetscCall(MatSetValues(jac,i1,ii,i1,ii,A,INSERT_VALUES,ierr))
 10   continue
      PetscCall(VecRestoreArrayRead(x,vxx,ixx,ierr))
      PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr))
      PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr))
      return
      end

!/*TEST
!
!   test:
!      nsize: 2
!      args: -ksp_gmres_cgs_refinement_type refine_always -n 10 -snes_monitor_short
!      output_file: output/ex12_1.out
!
!TEST*/