Actual source code: ex22f90.F90
slepc-3.20.0 2023-09-29
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> ./ex22f90 [-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(MatSetUp(A,ierr))
92: PetscCallA(MatGetOwnershipRange(A,Istart,Iend,ierr))
93: coeffs(1) = 1.0/(h*h)
94: coeffs(2) = -2.0/(h*h)+aa
95: do i=Istart,Iend-1
96: if (i .gt. 0) then
97: PetscCallA(MatSetValue(A,i,i-1,coeffs(1),INSERT_VALUES,ierr))
98: endif
99: if (i .lt. n-1) then
100: PetscCallA(MatSetValue(A,i,i+1,coeffs(1),INSERT_VALUES,ierr))
101: endif
102: PetscCallA(MatSetValue(A,i,i,coeffs(2),INSERT_VALUES,ierr))
103: end do
104: PetscCallA(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr))
105: PetscCallA(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr))
106: PetscCallA(MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE,ierr))
108: ! ** B = diag(bb(xi))
109: PetscCallA(MatCreate(PETSC_COMM_WORLD,B,ierr))
110: PetscCallA(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr))
111: PetscCallA(MatSetFromOptions(B,ierr))
112: PetscCallA(MatSetUp(B,ierr))
113: PetscCallA(MatGetOwnershipRange(B,Istart,Iend,ierr))
114: do i=Istart,Iend-1
115: xi = (i+1)*h
116: bb = -4.1+xi*(1.0-exp(xi-PETSC_PI))
117: PetscCallA(MatSetValue(B,i,i,bb,INSERT_VALUES,ierr))
118: end do
119: PetscCallA(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr))
120: PetscCallA(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr))
121: PetscCallA(MatSetOption(B,MAT_HERMITIAN,PETSC_TRUE,ierr))
123: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124: ! Create problem functions, f1=-lambda, f2=1.0, f3=exp(-tau*lambda)
125: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127: PetscCallA(FNCreate(PETSC_COMM_WORLD,f1,ierr))
128: PetscCallA(FNSetType(f1,FNRATIONAL,ierr))
129: k = 2
130: coeffs(1) = -1.0
131: coeffs(2) = 0.0
132: PetscCallA(FNRationalSetNumerator(f1,k,coeffs,ierr))
134: PetscCallA(FNCreate(PETSC_COMM_WORLD,f2,ierr))
135: PetscCallA(FNSetType(f2,FNRATIONAL,ierr))
136: k = 1
137: coeffs(1) = 1.0
138: PetscCallA(FNRationalSetNumerator(f2,k,coeffs,ierr))
140: PetscCallA(FNCreate(PETSC_COMM_WORLD,f3,ierr))
141: PetscCallA(FNSetType(f3,FNEXP,ierr))
142: scal = -tau
143: PetscCallA(FNSetScale(f3,scal,one,ierr))
145: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146: ! Create the eigensolver and set various options
147: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149: ! ** Create eigensolver context
150: PetscCallA(NEPCreate(PETSC_COMM_WORLD,nep,ierr))
152: ! ** Set the split operator. Note that A is passed first so that
153: ! ** SUBSET_NONZERO_PATTERN can be used
154: k = 3
155: mats(1) = A
156: mats(2) = Id
157: mats(3) = B
158: funs(1) = f2
159: funs(2) = f1
160: funs(3) = f3
161: PetscCallA(NEPSetSplitOperator(nep,k,mats,funs,SUBSET_NONZERO_PATTERN,ierr))
162: PetscCallA(NEPSetProblemType(nep,NEP_GENERAL,ierr))
164: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
165: ! Customize nonlinear solver; set runtime options
166: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
168: tol = 1e-9
169: PetscCallA(NEPSetTolerances(nep,tol,PETSC_DEFAULT_INTEGER,ierr))
170: k = 1
171: PetscCallA(NEPSetDimensions(nep,k,PETSC_DEFAULT_INTEGER,PETSC_DEFAULT_INTEGER,ierr))
172: k = 0
173: PetscCallA(NEPRIISetLagPreconditioner(nep,k,ierr))
175: ! ** Set solver parameters at runtime
176: PetscCallA(NEPSetFromOptions(nep,ierr))
178: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
179: ! Solve the eigensystem
180: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
182: PetscCallA(NEPSolve(nep,ierr))
184: ! ** Optional: Get some information from the solver and display it
185: PetscCallA(NEPGetType(nep,tname,ierr))
186: if (rank .eq. 0) then
187: write(*,120) tname
188: endif
189: 120 format (' Solution method: ',A)
190: PetscCallA(NEPGetDimensions(nep,nev,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
191: if (rank .eq. 0) then
192: write(*,130) nev
193: endif
194: 130 format (' Number of requested eigenvalues:',I4)
196: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
197: ! Display solution and clean up
198: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
200: ! ** show detailed info unless -terse option is given by user
201: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-terse',terse,ierr))
202: if (terse) then
203: PetscCallA(NEPErrorView(nep,PEP_ERROR_RELATIVE,PETSC_NULL_VIEWER,ierr))
204: else
205: PetscCallA(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL,ierr))
206: PetscCallA(NEPConvergedReasonView(nep,PETSC_VIEWER_STDOUT_WORLD,ierr))
207: PetscCallA(NEPErrorView(nep,PEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD,ierr))
208: PetscCallA(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD,ierr))
209: endif
210: PetscCallA(NEPDestroy(nep,ierr))
211: PetscCallA(MatDestroy(Id,ierr))
212: PetscCallA(MatDestroy(A,ierr))
213: PetscCallA(MatDestroy(B,ierr))
214: PetscCallA(FNDestroy(f1,ierr))
215: PetscCallA(FNDestroy(f2,ierr))
216: PetscCallA(FNDestroy(f3,ierr))
217: PetscCallA(SlepcFinalize(ierr))
218: end
220: !/*TEST
221: !
222: ! test:
223: ! suffix: 1
224: ! args: -terse
225: ! requires: !single
226: ! filter: sed -e "s/[+-]0\.0*i//g"
227: !
228: !TEST*/