Line data Source code
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 :
11 : static char help[] = "Illustrates the use of a user-defined stopping test.\n\n"
12 : "This is based on ex22.\n"
13 : "The command line options are:\n"
14 : " -n <n>, where <n> = number of grid subdivisions.\n"
15 : " -tau <tau>, where <tau> is the delay parameter.\n\n";
16 :
17 : /*
18 : Solve parabolic partial differential equation with time delay tau
19 :
20 : u_t = u_xx + a*u(t) + b*u(t-tau)
21 : u(0,t) = u(pi,t) = 0
22 :
23 : with a = 20 and b(x) = -4.1+x*(1-exp(x-pi)).
24 :
25 : Discretization leads to a DDE of dimension n
26 :
27 : -u' = A*u(t) + B*u(t-tau)
28 :
29 : which results in the nonlinear eigenproblem
30 :
31 : (-lambda*I + A + exp(-tau*lambda)*B)*u = 0
32 : */
33 :
34 : #include <slepcnep.h>
35 :
36 : /*
37 : User-defined routines
38 : */
39 : PetscErrorCode MyStoppingTest(NEP,PetscInt,PetscInt,PetscInt,PetscInt,NEPConvergedReason*,void*);
40 :
41 : typedef struct {
42 : PetscInt lastnconv; /* last value of nconv; used in stopping test */
43 : PetscInt nreps; /* number of repetitions of nconv; used in stopping test */
44 : } CTX_DELAY;
45 :
46 1 : int main(int argc,char **argv)
47 : {
48 1 : NEP nep;
49 1 : Mat Id,A,B;
50 1 : FN f1,f2,f3;
51 1 : RG rg;
52 1 : CTX_DELAY *ctx;
53 1 : Mat mats[3];
54 1 : FN funs[3];
55 1 : PetscScalar coeffs[2],b;
56 1 : PetscInt n=128,Istart,Iend,i,mpd;
57 1 : PetscReal tau=0.001,h,a=20,xi;
58 1 : PetscBool terse;
59 1 : PetscViewer viewer;
60 :
61 1 : PetscFunctionBeginUser;
62 1 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
63 1 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
64 1 : PetscCall(PetscOptionsGetReal(NULL,NULL,"-tau",&tau,NULL));
65 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n1-D Delay Eigenproblem, n=%" PetscInt_FMT ", tau=%g\n\n",n,(double)tau));
66 1 : h = PETSC_PI/(PetscReal)(n+1);
67 :
68 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69 : Create nonlinear eigensolver context
70 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
71 :
72 1 : PetscCall(NEPCreate(PETSC_COMM_WORLD,&nep));
73 :
74 : /* Identity matrix */
75 1 : PetscCall(MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,1.0,&Id));
76 1 : PetscCall(MatSetOption(Id,MAT_HERMITIAN,PETSC_TRUE));
77 :
78 : /* A = 1/h^2*tridiag(1,-2,1) + a*I */
79 1 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
80 1 : PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n));
81 1 : PetscCall(MatSetFromOptions(A));
82 1 : PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
83 129 : for (i=Istart;i<Iend;i++) {
84 128 : if (i>0) PetscCall(MatSetValue(A,i,i-1,1.0/(h*h),INSERT_VALUES));
85 128 : if (i<n-1) PetscCall(MatSetValue(A,i,i+1,1.0/(h*h),INSERT_VALUES));
86 128 : PetscCall(MatSetValue(A,i,i,-2.0/(h*h)+a,INSERT_VALUES));
87 : }
88 1 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
89 1 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
90 1 : PetscCall(MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE));
91 :
92 : /* B = diag(b(xi)) */
93 1 : PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
94 1 : PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,n,n));
95 1 : PetscCall(MatSetFromOptions(B));
96 1 : PetscCall(MatGetOwnershipRange(B,&Istart,&Iend));
97 129 : for (i=Istart;i<Iend;i++) {
98 128 : xi = (i+1)*h;
99 128 : b = -4.1+xi*(1.0-PetscExpReal(xi-PETSC_PI));
100 128 : PetscCall(MatSetValues(B,1,&i,1,&i,&b,INSERT_VALUES));
101 : }
102 1 : PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
103 1 : PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
104 1 : PetscCall(MatSetOption(B,MAT_HERMITIAN,PETSC_TRUE));
105 :
106 : /* Functions: f1=-lambda, f2=1.0, f3=exp(-tau*lambda) */
107 1 : PetscCall(FNCreate(PETSC_COMM_WORLD,&f1));
108 1 : PetscCall(FNSetType(f1,FNRATIONAL));
109 1 : coeffs[0] = -1.0; coeffs[1] = 0.0;
110 1 : PetscCall(FNRationalSetNumerator(f1,2,coeffs));
111 :
112 1 : PetscCall(FNCreate(PETSC_COMM_WORLD,&f2));
113 1 : PetscCall(FNSetType(f2,FNRATIONAL));
114 1 : coeffs[0] = 1.0;
115 1 : PetscCall(FNRationalSetNumerator(f2,1,coeffs));
116 :
117 1 : PetscCall(FNCreate(PETSC_COMM_WORLD,&f3));
118 1 : PetscCall(FNSetType(f3,FNEXP));
119 1 : PetscCall(FNSetScale(f3,-tau,1.0));
120 :
121 : /* Set the split operator */
122 1 : mats[0] = A; funs[0] = f2;
123 1 : mats[1] = Id; funs[1] = f1;
124 1 : mats[2] = B; funs[2] = f3;
125 1 : PetscCall(NEPSetSplitOperator(nep,3,mats,funs,SUBSET_NONZERO_PATTERN));
126 :
127 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128 : Customize nonlinear solver; set runtime options
129 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130 :
131 1 : PetscCall(NEPSetType(nep,NEPNLEIGS));
132 1 : PetscCall(NEPGetRG(nep,&rg));
133 1 : PetscCall(RGSetType(rg,RGINTERVAL));
134 : #if defined(PETSC_USE_COMPLEX)
135 : PetscCall(RGIntervalSetEndpoints(rg,5,20,-0.001,0.001));
136 : #else
137 1 : PetscCall(RGIntervalSetEndpoints(rg,5,20,-0.0,0.0));
138 : #endif
139 1 : PetscCall(NEPSetTarget(nep,15.0));
140 1 : PetscCall(NEPSetWhichEigenpairs(nep,NEP_TARGET_MAGNITUDE));
141 :
142 : /*
143 : Set solver options. In particular, we must allocate sufficient
144 : storage for all eigenpairs that may converge (ncv). This is
145 : application-dependent.
146 : */
147 1 : mpd = 40;
148 1 : PetscCall(NEPSetDimensions(nep,2*mpd,3*mpd,mpd));
149 1 : PetscCall(NEPSetTolerances(nep,PETSC_CURRENT,2000));
150 1 : PetscCall(PetscNew(&ctx));
151 1 : ctx->lastnconv = 0;
152 1 : ctx->nreps = 0;
153 1 : PetscCall(NEPSetStoppingTestFunction(nep,MyStoppingTest,(void*)ctx,NULL));
154 :
155 1 : PetscCall(NEPSetFromOptions(nep));
156 :
157 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158 : Solve the eigensystem
159 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
160 :
161 1 : PetscCall(NEPSolve(nep));
162 :
163 : /* show detailed info unless -terse option is given by user */
164 1 : PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer));
165 1 : PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO_DETAIL));
166 1 : PetscCall(NEPConvergedReasonView(nep,viewer));
167 1 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
168 1 : if (!terse) PetscCall(NEPErrorView(nep,NEP_ERROR_BACKWARD,viewer));
169 1 : PetscCall(PetscViewerPopFormat(viewer));
170 :
171 1 : PetscCall(NEPDestroy(&nep));
172 1 : PetscCall(MatDestroy(&Id));
173 1 : PetscCall(MatDestroy(&A));
174 1 : PetscCall(MatDestroy(&B));
175 1 : PetscCall(FNDestroy(&f1));
176 1 : PetscCall(FNDestroy(&f2));
177 1 : PetscCall(FNDestroy(&f3));
178 1 : PetscCall(PetscFree(ctx));
179 1 : PetscCall(SlepcFinalize());
180 : return 0;
181 : }
182 :
183 : /*
184 : Function for user-defined stopping test.
185 :
186 : Ignores the value of nev. It only takes into account the number of
187 : eigenpairs that have converged in recent outer iterations (restarts);
188 : if no new eigenvalues have converged in the last few restarts,
189 : we stop the iteration, assuming that no more eigenvalues are present
190 : inside the region.
191 : */
192 12 : PetscErrorCode MyStoppingTest(NEP nep,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,NEPConvergedReason *reason,void *ptr)
193 : {
194 12 : CTX_DELAY *ctx = (CTX_DELAY*)ptr;
195 :
196 12 : PetscFunctionBeginUser;
197 : /* check usual termination conditions, but ignoring the case nconv>=nev */
198 12 : PetscCall(NEPStoppingBasic(nep,its,max_it,nconv,PETSC_INT_MAX,reason,NULL));
199 12 : if (*reason==NEP_CONVERGED_ITERATING) {
200 : /* check if nconv is the same as before */
201 12 : if (nconv==ctx->lastnconv) ctx->nreps++;
202 : else {
203 1 : ctx->lastnconv = nconv;
204 1 : ctx->nreps = 0;
205 : }
206 : /* check if no eigenvalues converged in last 10 restarts */
207 12 : if (nconv && ctx->nreps>10) *reason = NEP_CONVERGED_USER;
208 : }
209 12 : PetscFunctionReturn(PETSC_SUCCESS);
210 : }
211 :
212 : /*TEST
213 :
214 : test:
215 : suffix: 1
216 : args: -terse
217 :
218 : TEST*/
|
