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

!
!
!  This example demonstrates basic use of the SNES Fortran interface.
!
!
#include <petsc/finclude/petscsnes.h>
module ex12fmodule
  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) == 0) then
    one = 1
    PetscCall(SNESSetLagJacobian(snes, one, ierr))  ! build jacobian
  else
    mone = -1
    PetscCall(SNESSetLagJacobian(snes, mone, ierr)) ! do NOT build jacobian
  end if
  snesm%its = snesm%its + 1
end subroutine FormMonitor

!  Note: Any user-defined Fortran routines (such as FormJacobian)
!  MUST be declared as external.
!

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_ARRAY, 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_ARRAY, i0, PETSC_NULL_INTEGER_ARRAY, 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 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
  end do
  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
  PetscScalar h, d
  PetscScalar, pointer :: vf2(:), vxx(:), vff(:)

  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 > 1000) then
    print *, 'Local work array not big enough'
    call MPI_Abort(PETSC_COMM_WORLD, zero, ierr)
  end if

  PetscCall(VecGetArrayRead(ctx%xl, vxx, ierr))
  PetscCall(VecGetArray(f, vff, ierr))
  PetscCall(VecGetArray(ctx%F, vF2, 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 == 0) then
    s = 0
    vff(1) = vxx(1)
  else
    s = 1
  end if

  do i = 1, n - 2
    vff(i - s + 1) = d*(vxx(i) - 2.0*vxx(i + 1) + vxx(i + 2)) + vxx(i + 1)*vxx(i + 1) - vF2(i - s + 1)
  end do

  if (rank == size - 1) then
    vff(n - s) = vxx(n) - 1.0
  end if

  PetscCall(VecRestoreArray(f, vff, ierr))
  PetscCall(VecRestoreArrayRead(ctx%xl, vxx, ierr))
  PetscCall(VecRestoreArray(ctx%F, vF2, ierr))
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))
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
  PetscScalar d, A, h
  PetscScalar, pointer :: vxx(:)

  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, ierr))
  PetscCall(VecGetOwnershipRange(x, start, end, ierr))
  n = end - start

  if (rank == 0) then
    A = 1.0
    PetscCall(MatSetValues(jac, i1, [start], i1, [start], [A], INSERT_VALUES, ierr))
    istart = 1
  else
    istart = 0
  end if
  if (rank == 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
  end if
  do 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*vxx(i + 1)
    PetscCall(MatSetValues(jac, i1, [ii], i1, [ii], [A], INSERT_VALUES, ierr))
  end do
  PetscCall(VecRestoreArrayRead(x, vxx, ierr))
  PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr))
  PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr))
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*/