Actual source code: ex22f.F90

      1: !
      2: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      3: !  SLEPc - Scalable Library for Eigenvalue Problem Computations
      4: !  Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
      5: !
      6: !  This file is part of SLEPc.
      7: !  SLEPc is distributed under a 2-clause BSD license (see LICENSE).
      8: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      9: !
     10: !  Program usage: mpiexec -n <np> ./ex22f [-n <n>] [-tau <tau>] [SLEPc opts]
     11: !
     12: !  Description: Delay differential equation. Fortran90 equivalent of ex22.c
     13: !
     14: !  The command line options are:
     15: !    -n <n>, where <n> = number of grid subdivisions
     16: !    -tau <tau>, where <tau> = delay parameter
     17: !
     18: ! ----------------------------------------------------------------------
     19: !  Solve parabolic partial differential equation with time delay tau
     20: !
     21: !           u_t = u_xx + aa*u(t) + bb*u(t-tau)
     22: !           u(0,t) = u(pi,t) = 0
     23: !
     24: !  with aa = 20 and bb(x) = -4.1+x*(1-exp(x-pi)).
     25: !
     26: !  Discretization leads to a DDE of dimension n
     27: !
     28: !           -u' = A*u(t) + B*u(t-tau)
     29: !
     30: !  which results in the nonlinear eigenproblem
     31: !
     32: !           (-lambda*I + A + exp(-tau*lambda)*B)*u = 0
     33: ! ----------------------------------------------------------------------
     34: !
     35: #include <slepc/finclude/slepcnep.h>
     36: program ex22f
     37:   use slepcnep
     38:   implicit none

     40: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     41: ! Declarations
     42: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     44:   Mat            :: Id, A, B, mats(3)     ! problem matrices
     45:   FN             :: f1, f2, f3, funs(3)   ! functions to define the nonlinear operator
     46:   NEP            :: nep                   ! nonlinear eigensolver context
     47:   NEPType        :: tname
     48:   PetscScalar    :: bb, coeffs(2), scal
     49:   PetscReal      :: tau, h, aa, xi, tol
     50:   PetscInt       :: n, i, nev, Istart, Iend
     51:   PetscMPIInt    :: rank
     52:   PetscErrorCode :: ierr
     53:   PetscBool      :: flg, terse
     54:   PetscScalar, parameter :: one = 1.0

     56: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     57: ! Beginning of program
     58: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     60:   PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER, ierr))
     61:   PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
     62:   n = 128
     63:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))
     64:   tau = 0.001
     65:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-tau', tau, flg, ierr))
     66:   if (rank == 0) then
     67:     write (*, '(/a,i4,a,f6.3)') 'Delay Eigenproblem, n =', n, ', tau =', tau
     68:   end if

     70:   aa = 20.0
     71:   h = PETSC_PI/real(n + 1, PETSC_REAL_KIND)

     73: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     74: ! Create problem matrices
     75: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     77: ! ** Id is the identity matrix
     78:   PetscCallA(MatCreateConstantDiagonal(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n, n, one, Id, ierr))
     79:   PetscCallA(MatSetOption(Id, MAT_HERMITIAN, PETSC_TRUE, ierr))

     81: ! ** A = 1/h^2*tridiag(1,-2,1) + aa*I
     82:   PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
     83:   PetscCallA(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n, ierr))
     84:   PetscCallA(MatSetFromOptions(A, ierr))
     85:   PetscCallA(MatGetOwnershipRange(A, Istart, Iend, ierr))
     86:   coeffs(1) = 1.0/(h*h)
     87:   coeffs(2) = -2.0/(h*h) + aa
     88:   do i = Istart, Iend - 1
     89:     if (i > 0) then
     90:       PetscCallA(MatSetValue(A, i, i - 1, coeffs(1), INSERT_VALUES, ierr))
     91:     end if
     92:     if (i < n - 1) then
     93:       PetscCallA(MatSetValue(A, i, i + 1, coeffs(1), INSERT_VALUES, ierr))
     94:     end if
     95:     PetscCallA(MatSetValue(A, i, i, coeffs(2), INSERT_VALUES, ierr))
     96:   end do
     97:   PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
     98:   PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))
     99:   PetscCallA(MatSetOption(A, MAT_HERMITIAN, PETSC_TRUE, ierr))

    101: ! ** B = diag(bb(xi))
    102:   PetscCallA(MatCreate(PETSC_COMM_WORLD, B, ierr))
    103:   PetscCallA(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, n, n, ierr))
    104:   PetscCallA(MatSetFromOptions(B, ierr))
    105:   PetscCallA(MatGetOwnershipRange(B, Istart, Iend, ierr))
    106:   do i = Istart, Iend - 1
    107:     xi = (i + 1)*h
    108:     bb = -4.1 + xi*(1.0 - exp(xi - PETSC_PI))
    109:     PetscCallA(MatSetValue(B, i, i, bb, INSERT_VALUES, ierr))
    110:   end do
    111:   PetscCallA(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY, ierr))
    112:   PetscCallA(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY, ierr))
    113:   PetscCallA(MatSetOption(B, MAT_HERMITIAN, PETSC_TRUE, ierr))

    115: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    116: ! Create problem functions, f1=-lambda, f2=1.0, f3=exp(-tau*lambda)
    117: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    119:   PetscCallA(FNCreate(PETSC_COMM_WORLD, f1, ierr))
    120:   PetscCallA(FNSetType(f1, FNRATIONAL, ierr))
    121:   coeffs(1) = -1.0
    122:   coeffs(2) = 0.0
    123:   PetscCallA(FNRationalSetNumerator(f1, 2_PETSC_INT_KIND, coeffs, ierr))

    125:   PetscCallA(FNCreate(PETSC_COMM_WORLD, f2, ierr))
    126:   PetscCallA(FNSetType(f2, FNRATIONAL, ierr))
    127:   coeffs(1) = 1.0
    128:   PetscCallA(FNRationalSetNumerator(f2, 1_PETSC_INT_KIND, coeffs, ierr))

    130:   PetscCallA(FNCreate(PETSC_COMM_WORLD, f3, ierr))
    131:   PetscCallA(FNSetType(f3, FNEXP, ierr))
    132:   scal = -tau
    133:   PetscCallA(FNSetScale(f3, scal, one, ierr))

    135: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    136: ! Create the eigensolver and set various options
    137: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    139: ! ** Create eigensolver context
    140:   PetscCallA(NEPCreate(PETSC_COMM_WORLD, nep, ierr))

    142: ! ** Set the split operator. Note that A is passed first so that
    143: ! ** SUBSET_NONZERO_PATTERN can be used
    144:   mats(1) = A
    145:   mats(2) = Id
    146:   mats(3) = B
    147:   funs(1) = f2
    148:   funs(2) = f1
    149:   funs(3) = f3
    150:   PetscCallA(NEPSetSplitOperator(nep, 3_PETSC_INT_KIND, mats, funs, SUBSET_NONZERO_PATTERN, ierr))
    151:   PetscCallA(NEPSetProblemType(nep, NEP_GENERAL, ierr))

    153: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    154: ! Customize nonlinear solver; set runtime options
    155: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    157:   tol = 1e-9
    158:   PetscCallA(NEPSetTolerances(nep, tol, PETSC_CURRENT_INTEGER, ierr))
    159:   PetscCallA(NEPSetDimensions(nep, 1_PETSC_INT_KIND, PETSC_DETERMINE_INTEGER, PETSC_DETERMINE_INTEGER, ierr))
    160:   PetscCallA(NEPRIISetLagPreconditioner(nep, 0_PETSC_INT_KIND, ierr))

    162: ! ** Set solver parameters at runtime
    163:   PetscCallA(NEPSetFromOptions(nep, ierr))

    165: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    166: ! Solve the eigensystem
    167: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    169:   PetscCallA(NEPSolve(nep, ierr))

    171: ! ** Optional: Get some information from the solver and display it
    172:   PetscCallA(NEPGetType(nep, tname, ierr))
    173:   if (rank == 0) then
    174:     write (*, '(a,a)') ' Solution method: ', tname
    175:   end if
    176:   PetscCallA(NEPGetDimensions(nep, nev, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
    177:   if (rank == 0) then
    178:     write (*, '(a,i4)') ' Number of requested eigenvalues:', nev
    179:   end if

    181: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    182: ! Display solution and clean up
    183: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    185: ! ** show detailed info unless -terse option is given by user
    186:   PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-terse', terse, ierr))
    187:   if (terse) then
    188:     PetscCallA(NEPErrorView(nep, NEP_ERROR_RELATIVE, PETSC_NULL_VIEWER, ierr))
    189:   else
    190:     PetscCallA(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL, ierr))
    191:     PetscCallA(NEPConvergedReasonView(nep, PETSC_VIEWER_STDOUT_WORLD, ierr))
    192:     PetscCallA(NEPErrorView(nep, NEP_ERROR_RELATIVE, PETSC_VIEWER_STDOUT_WORLD, ierr))
    193:     PetscCallA(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD, ierr))
    194:   end if
    195:   PetscCallA(NEPDestroy(nep, ierr))
    196:   PetscCallA(MatDestroy(Id, ierr))
    197:   PetscCallA(MatDestroy(A, ierr))
    198:   PetscCallA(MatDestroy(B, ierr))
    199:   PetscCallA(FNDestroy(f1, ierr))
    200:   PetscCallA(FNDestroy(f2, ierr))
    201:   PetscCallA(FNDestroy(f3, ierr))
    202:   PetscCallA(SlepcFinalize(ierr))
    203: end program ex22f

    205: !/*TEST
    206: !
    207: !   test:
    208: !      suffix: 1
    209: !      args: -terse
    210: !      requires: !single
    211: !      filter: sed -e "s/[+-]0\.0*i//g"
    212: !
    213: !TEST*/