Actual source code: test10.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
11: static char help[] = "Tests multiple calls to NEPSolve(). Based on ex22.c.\n\n"
12: "The command line options are:\n"
13: " -n <n>, where <n> = number of grid subdivisions.\n"
14: " -tau <tau>, where <tau> is the delay parameter.\n\n";
16: /*
17: Solve parabolic partial differential equation with time delay tau
19: u_t = u_xx + a*u(t) + b*u(t-tau)
20: u(0,t) = u(pi,t) = 0
22: with a = 20 and b(x) = -4.1+x*(1-exp(x-pi)).
24: Discretization leads to a DDE of dimension n
26: -u' = A*u(t) + B*u(t-tau)
28: which results in the nonlinear eigenproblem
30: (-lambda*I + A + exp(-tau*lambda)*B)*u = 0
31: */
33: #include <slepcnep.h>
35: int main(int argc,char **argv)
36: {
37: NEP nep; /* nonlinear eigensolver context */
38: Mat Id,A,B; /* problem matrices */
39: FN f1,f2,f3; /* functions to define the nonlinear operator */
40: Mat mats[3];
41: FN funs[3];
42: PetscScalar coeffs[2],b;
43: PetscInt n=128,Istart,Iend,i;
44: PetscReal tau=0.001,h,a=20,xi;
47: SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
48: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
49: PetscOptionsGetReal(NULL,NULL,"-tau",&tau,NULL);
50: PetscPrintf(PETSC_COMM_WORLD,"\n1-D Delay Eigenproblem, n=%D, tau=%g\n\n",n,(double)tau);
51: h = PETSC_PI/(PetscReal)(n+1);
53: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
54: Create functions that define the split operator
55: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
57: /* f1=-lambda */
58: FNCreate(PETSC_COMM_WORLD,&f1);
59: FNSetType(f1,FNRATIONAL);
60: coeffs[0] = -1.0; coeffs[1] = 0.0;
61: FNRationalSetNumerator(f1,2,coeffs);
63: /* f2=1.0 */
64: FNCreate(PETSC_COMM_WORLD,&f2);
65: FNSetType(f2,FNRATIONAL);
66: coeffs[0] = 1.0;
67: FNRationalSetNumerator(f2,1,coeffs);
69: /* f3=exp(-tau*lambda) */
70: FNCreate(PETSC_COMM_WORLD,&f3);
71: FNSetType(f3,FNEXP);
72: FNSetScale(f3,-tau,1.0);
74: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
75: Create problem matrices
76: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
78: /* Identity matrix */
79: MatCreate(PETSC_COMM_WORLD,&Id);
80: MatSetSizes(Id,PETSC_DECIDE,PETSC_DECIDE,n,n);
81: MatSetFromOptions(Id);
82: MatSetUp(Id);
83: MatGetOwnershipRange(Id,&Istart,&Iend);
84: for (i=Istart;i<Iend;i++) {
85: MatSetValue(Id,i,i,1.0,INSERT_VALUES);
86: }
87: MatAssemblyBegin(Id,MAT_FINAL_ASSEMBLY);
88: MatAssemblyEnd(Id,MAT_FINAL_ASSEMBLY);
89: MatSetOption(Id,MAT_HERMITIAN,PETSC_TRUE);
91: /* A = 1/h^2*tridiag(1,-2,1) + a*I */
92: MatCreate(PETSC_COMM_WORLD,&A);
93: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
94: MatSetFromOptions(A);
95: MatSetUp(A);
96: MatGetOwnershipRange(A,&Istart,&Iend);
97: for (i=Istart;i<Iend;i++) {
98: if (i>0) { MatSetValue(A,i,i-1,1.0/(h*h),INSERT_VALUES); }
99: if (i<n-1) { MatSetValue(A,i,i+1,1.0/(h*h),INSERT_VALUES); }
100: MatSetValue(A,i,i,-2.0/(h*h)+a,INSERT_VALUES);
101: }
102: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
103: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
104: MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE);
106: /* B = diag(b(xi)) */
107: MatCreate(PETSC_COMM_WORLD,&B);
108: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,n,n);
109: MatSetFromOptions(B);
110: MatSetUp(B);
111: MatGetOwnershipRange(B,&Istart,&Iend);
112: for (i=Istart;i<Iend;i++) {
113: xi = (i+1)*h;
114: b = -4.1+xi*(1.0-PetscExpReal(xi-PETSC_PI));
115: MatSetValue(B,i,i,b,INSERT_VALUES);
116: }
117: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
118: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
119: MatSetOption(B,MAT_HERMITIAN,PETSC_TRUE);
121: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122: Create nonlinear eigensolver and set options
123: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
125: NEPCreate(PETSC_COMM_WORLD,&nep);
126: mats[0] = A; funs[0] = f2;
127: mats[1] = Id; funs[1] = f1;
128: mats[2] = B; funs[2] = f3;
129: NEPSetSplitOperator(nep,3,mats,funs,SUBSET_NONZERO_PATTERN);
130: NEPSetTolerances(nep,1e-9,PETSC_DEFAULT);
131: NEPSetFromOptions(nep);
133: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134: Solve the eigensystem
135: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
137: NEPSolve(nep);
138: NEPErrorView(nep,NEP_ERROR_RELATIVE,NULL);
140: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141: Create problem matrices of size 2*n
142: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
144: MatDestroy(&Id);
145: MatDestroy(&A);
146: MatDestroy(&B);
147: n *= 2;
148: PetscPrintf(PETSC_COMM_WORLD,"\n1-D Delay Eigenproblem, n=%D, tau=%g\n\n",n,(double)tau);
149: h = PETSC_PI/(PetscReal)(n+1);
151: /* Identity matrix */
152: MatCreate(PETSC_COMM_WORLD,&Id);
153: MatSetSizes(Id,PETSC_DECIDE,PETSC_DECIDE,n,n);
154: MatSetFromOptions(Id);
155: MatSetUp(Id);
156: MatGetOwnershipRange(Id,&Istart,&Iend);
157: for (i=Istart;i<Iend;i++) {
158: MatSetValue(Id,i,i,1.0,INSERT_VALUES);
159: }
160: MatAssemblyBegin(Id,MAT_FINAL_ASSEMBLY);
161: MatAssemblyEnd(Id,MAT_FINAL_ASSEMBLY);
162: MatSetOption(Id,MAT_HERMITIAN,PETSC_TRUE);
164: /* A = 1/h^2*tridiag(1,-2,1) + a*I */
165: MatCreate(PETSC_COMM_WORLD,&A);
166: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
167: MatSetFromOptions(A);
168: MatSetUp(A);
169: MatGetOwnershipRange(A,&Istart,&Iend);
170: for (i=Istart;i<Iend;i++) {
171: if (i>0) { MatSetValue(A,i,i-1,1.0/(h*h),INSERT_VALUES); }
172: if (i<n-1) { MatSetValue(A,i,i+1,1.0/(h*h),INSERT_VALUES); }
173: MatSetValue(A,i,i,-2.0/(h*h)+a,INSERT_VALUES);
174: }
175: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
176: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
177: MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE);
179: /* B = diag(b(xi)) */
180: MatCreate(PETSC_COMM_WORLD,&B);
181: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,n,n);
182: MatSetFromOptions(B);
183: MatSetUp(B);
184: MatGetOwnershipRange(B,&Istart,&Iend);
185: for (i=Istart;i<Iend;i++) {
186: xi = (i+1)*h;
187: b = -4.1+xi*(1.0-PetscExpReal(xi-PETSC_PI));
188: MatSetValue(B,i,i,b,INSERT_VALUES);
189: }
190: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
191: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
192: MatSetOption(B,MAT_HERMITIAN,PETSC_TRUE);
194: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
195: Solve again, calling NEPReset() since matrix size has changed
196: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
198: NEPReset(nep); /* if this is omitted, it will be called in NEPSetSplitOperators() */
199: mats[0] = A; funs[0] = f2;
200: mats[1] = Id; funs[1] = f1;
201: mats[2] = B; funs[2] = f3;
202: NEPSetSplitOperator(nep,3,mats,funs,SUBSET_NONZERO_PATTERN);
203: NEPSolve(nep);
204: NEPErrorView(nep,NEP_ERROR_RELATIVE,NULL);
206: NEPDestroy(&nep);
207: MatDestroy(&Id);
208: MatDestroy(&A);
209: MatDestroy(&B);
210: FNDestroy(&f1);
211: FNDestroy(&f2);
212: FNDestroy(&f3);
213: SlepcFinalize();
214: return ierr;
215: }
217: /*TEST
219: testset:
220: nsize: 2
221: requires: !single
222: output_file: output/test10_1.out
223: test:
224: suffix: 1
225: args: -nep_type {{rii narnoldi}} -nep_target 0.55
226: test:
227: suffix: 1_slp
228: args: -nep_type slp -nep_slp_st_pc_type redundant
229: test:
230: suffix: 1_interpol
231: args: -nep_type interpol -rg_type interval -rg_interval_endpoints .5,1,-.1,.1 -nep_target .7 -nep_interpol_st_pc_type redundant
232: test:
233: suffix: 1_narnoldi_sync
234: args: -nep_type narnoldi -ds_parallel synchronized
236: test:
237: suffix: 2
238: args: -nep_nev 2 -nep_type interpol -rg_type interval -rg_interval_endpoints .5,15,-.1,.1 -nep_target .7 -nep_interpol_pep_type jd -nep_interpol_st_pc_type sor
239: requires: !single
241: TEST*/