Actual source code: ex22f.F90
slepc-3.21.1 2024-04-26
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: program main
36: #include <slepc/finclude/slepcnep.h>
37: use slepcnep
38: implicit none
40: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
41: ! Declarations
42: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
43: !
44: ! Variables:
45: ! nep nonlinear eigensolver context
46: ! Id,A,B problem matrices
47: ! f1,f2,f3 functions to define the nonlinear operator
49: Mat Id, A, B, mats(3)
50: FN f1, f2, f3, funs(3)
51: NEP nep
52: NEPType tname
53: PetscScalar one, bb, coeffs(2), scal
54: PetscReal tau, h, aa, xi, tol
55: PetscInt n, i, k, nev, Istart, Iend
56: PetscMPIInt rank
57: PetscErrorCode ierr
58: PetscBool flg, terse
60: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
61: ! Beginning of program
62: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
64: PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER,ierr))
65: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
66: n = 128
67: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr))
68: tau = 0.001
69: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-tau',tau,flg,ierr))
70: if (rank .eq. 0) then
71: write(*,100) n, tau
72: endif
73: 100 format (/'Delay Eigenproblem, n =',I4,', tau =',F6.3)
75: one = 1.0
76: aa = 20.0
77: h = PETSC_PI/real(n+1)
79: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
80: ! Create problem matrices
81: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
83: ! ** Id is the identity matrix
84: PetscCallA(MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,one,Id,ierr))
85: PetscCallA(MatSetOption(Id,MAT_HERMITIAN,PETSC_TRUE,ierr))
87: ! ** A = 1/h^2*tridiag(1,-2,1) + aa*I
88: PetscCallA(MatCreate(PETSC_COMM_WORLD,A,ierr))
89: PetscCallA(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr))
90: PetscCallA(MatSetFromOptions(A,ierr))
91: PetscCallA(MatGetOwnershipRange(A,Istart,Iend,ierr))
92: coeffs(1) = 1.0/(h*h)
93: coeffs(2) = -2.0/(h*h)+aa
94: do i=Istart,Iend-1
95: if (i .gt. 0) then
96: PetscCallA(MatSetValue(A,i,i-1,coeffs(1),INSERT_VALUES,ierr))
97: endif
98: if (i .lt. n-1) then
99: PetscCallA(MatSetValue(A,i,i+1,coeffs(1),INSERT_VALUES,ierr))
100: endif
101: PetscCallA(MatSetValue(A,i,i,coeffs(2),INSERT_VALUES,ierr))
102: end do
103: PetscCallA(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr))
104: PetscCallA(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr))
105: PetscCallA(MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE,ierr))
107: ! ** B = diag(bb(xi))
108: PetscCallA(MatCreate(PETSC_COMM_WORLD,B,ierr))
109: PetscCallA(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr))
110: PetscCallA(MatSetFromOptions(B,ierr))
111: PetscCallA(MatGetOwnershipRange(B,Istart,Iend,ierr))
112: do i=Istart,Iend-1
113: xi = (i+1)*h
114: bb = -4.1+xi*(1.0-exp(xi-PETSC_PI))
115: PetscCallA(MatSetValue(B,i,i,bb,INSERT_VALUES,ierr))
116: end do
117: PetscCallA(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr))
118: PetscCallA(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr))
119: PetscCallA(MatSetOption(B,MAT_HERMITIAN,PETSC_TRUE,ierr))
121: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122: ! Create problem functions, f1=-lambda, f2=1.0, f3=exp(-tau*lambda)
123: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125: PetscCallA(FNCreate(PETSC_COMM_WORLD,f1,ierr))
126: PetscCallA(FNSetType(f1,FNRATIONAL,ierr))
127: k = 2
128: coeffs(1) = -1.0
129: coeffs(2) = 0.0
130: PetscCallA(FNRationalSetNumerator(f1,k,coeffs,ierr))
132: PetscCallA(FNCreate(PETSC_COMM_WORLD,f2,ierr))
133: PetscCallA(FNSetType(f2,FNRATIONAL,ierr))
134: k = 1
135: coeffs(1) = 1.0
136: PetscCallA(FNRationalSetNumerator(f2,k,coeffs,ierr))
138: PetscCallA(FNCreate(PETSC_COMM_WORLD,f3,ierr))
139: PetscCallA(FNSetType(f3,FNEXP,ierr))
140: scal = -tau
141: PetscCallA(FNSetScale(f3,scal,one,ierr))
143: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144: ! Create the eigensolver and set various options
145: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: ! ** Create eigensolver context
148: PetscCallA(NEPCreate(PETSC_COMM_WORLD,nep,ierr))
150: ! ** Set the split operator. Note that A is passed first so that
151: ! ** SUBSET_NONZERO_PATTERN can be used
152: k = 3
153: mats(1) = A
154: mats(2) = Id
155: mats(3) = B
156: funs(1) = f2
157: funs(2) = f1
158: funs(3) = f3
159: PetscCallA(NEPSetSplitOperator(nep,k,mats,funs,SUBSET_NONZERO_PATTERN,ierr))
160: PetscCallA(NEPSetProblemType(nep,NEP_GENERAL,ierr))
162: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
163: ! Customize nonlinear solver; set runtime options
164: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
166: tol = 1e-9
167: PetscCallA(NEPSetTolerances(nep,tol,PETSC_DEFAULT_INTEGER,ierr))
168: k = 1
169: PetscCallA(NEPSetDimensions(nep,k,PETSC_DEFAULT_INTEGER,PETSC_DEFAULT_INTEGER,ierr))
170: k = 0
171: PetscCallA(NEPRIISetLagPreconditioner(nep,k,ierr))
173: ! ** Set solver parameters at runtime
174: PetscCallA(NEPSetFromOptions(nep,ierr))
176: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
177: ! Solve the eigensystem
178: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
180: PetscCallA(NEPSolve(nep,ierr))
182: ! ** Optional: Get some information from the solver and display it
183: PetscCallA(NEPGetType(nep,tname,ierr))
184: if (rank .eq. 0) then
185: write(*,120) tname
186: endif
187: 120 format (' Solution method: ',A)
188: PetscCallA(NEPGetDimensions(nep,nev,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
189: if (rank .eq. 0) then
190: write(*,130) nev
191: endif
192: 130 format (' Number of requested eigenvalues:',I4)
194: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
195: ! Display solution and clean up
196: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
198: ! ** show detailed info unless -terse option is given by user
199: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-terse',terse,ierr))
200: if (terse) then
201: PetscCallA(NEPErrorView(nep,PEP_ERROR_RELATIVE,PETSC_NULL_VIEWER,ierr))
202: else
203: PetscCallA(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL,ierr))
204: PetscCallA(NEPConvergedReasonView(nep,PETSC_VIEWER_STDOUT_WORLD,ierr))
205: PetscCallA(NEPErrorView(nep,PEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD,ierr))
206: PetscCallA(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD,ierr))
207: endif
208: PetscCallA(NEPDestroy(nep,ierr))
209: PetscCallA(MatDestroy(Id,ierr))
210: PetscCallA(MatDestroy(A,ierr))
211: PetscCallA(MatDestroy(B,ierr))
212: PetscCallA(FNDestroy(f1,ierr))
213: PetscCallA(FNDestroy(f2,ierr))
214: PetscCallA(FNDestroy(f3,ierr))
215: PetscCallA(SlepcFinalize(ierr))
216: end
218: !/*TEST
219: !
220: ! test:
221: ! suffix: 1
222: ! args: -terse
223: ! requires: !single
224: ! filter: sed -e "s/[+-]0\.0*i//g"
225: !
226: !TEST*/