Actual source code: ex22.c
slepc-3.6.1 2015-09-03
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2015, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
8: SLEPc is free software: you can redistribute it and/or modify it under the
9: terms of version 3 of the GNU Lesser General Public License as published by
10: the Free Software Foundation.
12: SLEPc is distributed in the hope that it will be useful, but WITHOUT ANY
13: WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14: FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15: more details.
17: You should have received a copy of the GNU Lesser General Public License
18: along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
19: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
20: */
22: static char help[] = "Delay differential equation.\n\n"
23: "The command line options are:\n"
24: " -n <n>, where <n> = number of grid subdivisions.\n"
25: " -tau <tau>, where <tau> is the delay parameter.\n\n";
27: /*
28: Solve parabolic partial differential equation with time delay tau
30: u_t = u_xx + a*u(t) + b*u(t-tau)
31: u(0,t) = u(pi,t) = 0
33: with a = 20 and b(x) = -4.1+x*(1-exp(x-pi)).
35: Discretization leads to a DDE of dimension n
37: -u' = A*u(t) + B*u(t-tau)
39: which results in the nonlinear eigenproblem
41: (-lambda*I + A + exp(-tau*lambda)*B)*u = 0
42: */
44: #include <slepcnep.h>
48: int main(int argc,char **argv)
49: {
50: NEP nep; /* nonlinear eigensolver context */
51: Mat Id,A,B; /* problem matrices */
52: FN f1,f2,f3; /* functions to define the nonlinear operator */
53: Mat mats[3];
54: FN funs[3];
55: NEPType type;
56: PetscScalar coeffs[2],b;
57: PetscInt n=128,nev,Istart,Iend,i,its;
58: PetscReal tau=0.001,h,a=20,xi;
59: PetscBool terse;
62: SlepcInitialize(&argc,&argv,(char*)0,help);
63: PetscOptionsGetInt(NULL,"-n",&n,NULL);
64: PetscOptionsGetReal(NULL,"-tau",&tau,NULL);
65: PetscPrintf(PETSC_COMM_WORLD,"\n1-D Delay Eigenproblem, n=%D, tau=%g\n\n",n,(double)tau);
66: h = PETSC_PI/(PetscReal)(n+1);
68: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69: Create nonlinear eigensolver context
70: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
72: NEPCreate(PETSC_COMM_WORLD,&nep);
74: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
75: Create problem matrices and coefficient functions. Pass them to NEP
76: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
78: /*
79: Identity matrix
80: */
81: MatCreate(PETSC_COMM_WORLD,&Id);
82: MatSetSizes(Id,PETSC_DECIDE,PETSC_DECIDE,n,n);
83: MatSetFromOptions(Id);
84: MatSetUp(Id);
85: MatGetOwnershipRange(Id,&Istart,&Iend);
86: for (i=Istart;i<Iend;i++) {
87: MatSetValue(Id,i,i,1.0,INSERT_VALUES);
88: }
89: MatAssemblyBegin(Id,MAT_FINAL_ASSEMBLY);
90: MatAssemblyEnd(Id,MAT_FINAL_ASSEMBLY);
91: MatSetOption(Id,MAT_HERMITIAN,PETSC_TRUE);
93: /*
94: A = 1/h^2*tridiag(1,-2,1) + a*I
95: */
96: MatCreate(PETSC_COMM_WORLD,&A);
97: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
98: MatSetFromOptions(A);
99: MatSetUp(A);
100: MatGetOwnershipRange(A,&Istart,&Iend);
101: for (i=Istart;i<Iend;i++) {
102: if (i>0) { MatSetValue(A,i,i-1,1.0/(h*h),INSERT_VALUES); }
103: if (i<n-1) { MatSetValue(A,i,i+1,1.0/(h*h),INSERT_VALUES); }
104: MatSetValue(A,i,i,-2.0/(h*h)+a,INSERT_VALUES);
105: }
106: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
107: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
108: MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE);
110: /*
111: B = diag(b(xi))
112: */
113: MatCreate(PETSC_COMM_WORLD,&B);
114: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,n,n);
115: MatSetFromOptions(B);
116: MatSetUp(B);
117: MatGetOwnershipRange(B,&Istart,&Iend);
118: for (i=Istart;i<Iend;i++) {
119: xi = (i+1)*h;
120: b = -4.1+xi*(1.0-PetscExpReal(xi-PETSC_PI));
121: MatSetValues(B,1,&i,1,&i,&b,INSERT_VALUES);
122: }
123: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
124: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
125: MatSetOption(B,MAT_HERMITIAN,PETSC_TRUE);
127: /*
128: Functions: f1=-lambda, f2=1.0, f3=exp(-tau*lambda)
129: */
130: FNCreate(PETSC_COMM_WORLD,&f1);
131: FNSetType(f1,FNRATIONAL);
132: coeffs[0] = -1.0; coeffs[1] = 0.0;
133: FNRationalSetNumerator(f1,2,coeffs);
135: FNCreate(PETSC_COMM_WORLD,&f2);
136: FNSetType(f2,FNRATIONAL);
137: coeffs[0] = 1.0;
138: FNRationalSetNumerator(f2,1,coeffs);
140: FNCreate(PETSC_COMM_WORLD,&f3);
141: FNSetType(f3,FNEXP);
142: FNSetScale(f3,-tau,1.0);
144: /*
145: Set the split operator. Note that A is passed first so that
146: SUBSET_NONZERO_PATTERN can be used
147: */
148: mats[0] = A; funs[0] = f2;
149: mats[1] = Id; funs[1] = f1;
150: mats[2] = B; funs[2] = f3;
151: NEPSetSplitOperator(nep,3,mats,funs,SUBSET_NONZERO_PATTERN);
153: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
154: Customize nonlinear solver; set runtime options
155: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
157: NEPSetTolerances(nep,PETSC_DEFAULT,1e-9,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
158: NEPSetDimensions(nep,1,PETSC_DEFAULT,PETSC_DEFAULT);
159: NEPSetLagPreconditioner(nep,0);
161: /*
162: Set solver parameters at runtime
163: */
164: NEPSetFromOptions(nep);
166: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167: Solve the eigensystem
168: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
170: NEPSolve(nep);
171: NEPGetIterationNumber(nep,&its);
172: PetscPrintf(PETSC_COMM_WORLD," Number of NEP iterations = %D\n\n",its);
174: /*
175: Optional: Get some information from the solver and display it
176: */
177: NEPGetType(nep,&type);
178: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n",type);
179: NEPGetDimensions(nep,&nev,NULL,NULL);
180: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
182: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183: Display solution and clean up
184: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
186: /* show detailed info unless -terse option is given by user */
187: PetscOptionsHasName(NULL,"-terse",&terse);
188: if (terse) {
189: NEPErrorView(nep,NEP_ERROR_RELATIVE,NULL);
190: } else {
191: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
192: NEPReasonView(nep,PETSC_VIEWER_STDOUT_WORLD);
193: NEPErrorView(nep,NEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);
194: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
195: }
196: NEPDestroy(&nep);
197: MatDestroy(&Id);
198: MatDestroy(&A);
199: MatDestroy(&B);
200: FNDestroy(&f1);
201: FNDestroy(&f2);
202: FNDestroy(&f3);
203: SlepcFinalize();
204: return 0;
205: }