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[] = "Test the CISS solver with the problem of ex22.\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";
15 :
16 : /*
17 : Solve parabolic partial differential equation with time delay tau
18 :
19 : u_t = u_xx + a*u(t) + b*u(t-tau)
20 : u(0,t) = u(pi,t) = 0
21 :
22 : with a = 20 and b(x) = -4.1+x*(1-exp(x-pi)).
23 :
24 : Discretization leads to a DDE of dimension n
25 :
26 : -u' = A*u(t) + B*u(t-tau)
27 :
28 : which results in the nonlinear eigenproblem
29 :
30 : (-lambda*I + A + exp(-tau*lambda)*B)*u = 0
31 : */
32 :
33 : #include <slepcnep.h>
34 :
35 13 : int main(int argc,char **argv)
36 : {
37 13 : NEP nep;
38 13 : Mat Id,A,B,mats[3];
39 13 : FN f1,f2,f3,funs[3];
40 13 : RG rg;
41 13 : KSP *ksp;
42 13 : PC pc;
43 13 : NEPCISSExtraction ext;
44 13 : PetscScalar coeffs[2],b;
45 13 : PetscInt n=128,Istart,Iend,i,nsolve;
46 13 : PetscReal tau=0.001,h,a=20,xi;
47 13 : PetscBool flg,terse;
48 :
49 13 : PetscFunctionBeginUser;
50 13 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
51 13 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
52 13 : PetscCall(PetscOptionsGetReal(NULL,NULL,"-tau",&tau,NULL));
53 13 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n1-D Delay Eigenproblem, n=%" PetscInt_FMT ", tau=%g\n\n",n,(double)tau));
54 13 : h = PETSC_PI/(PetscReal)(n+1);
55 :
56 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
57 : Create nonlinear eigensolver context
58 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
59 :
60 13 : PetscCall(NEPCreate(PETSC_COMM_WORLD,&nep));
61 :
62 : /* Identity matrix */
63 13 : PetscCall(MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,1.0,&Id));
64 13 : PetscCall(MatSetOption(Id,MAT_HERMITIAN,PETSC_TRUE));
65 :
66 : /* A = 1/h^2*tridiag(1,-2,1) + a*I */
67 13 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
68 13 : PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n));
69 13 : PetscCall(MatSetFromOptions(A));
70 13 : PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
71 1293 : for (i=Istart;i<Iend;i++) {
72 1280 : if (i>0) PetscCall(MatSetValue(A,i,i-1,1.0/(h*h),INSERT_VALUES));
73 1280 : if (i<n-1) PetscCall(MatSetValue(A,i,i+1,1.0/(h*h),INSERT_VALUES));
74 1280 : PetscCall(MatSetValue(A,i,i,-2.0/(h*h)+a,INSERT_VALUES));
75 : }
76 13 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
77 13 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
78 13 : PetscCall(MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE));
79 :
80 : /* B = diag(b(xi)) */
81 13 : PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
82 13 : PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,n,n));
83 13 : PetscCall(MatSetFromOptions(B));
84 13 : PetscCall(MatGetOwnershipRange(B,&Istart,&Iend));
85 1293 : for (i=Istart;i<Iend;i++) {
86 1280 : xi = (i+1)*h;
87 1280 : b = -4.1+xi*(1.0-PetscExpReal(xi-PETSC_PI));
88 1280 : PetscCall(MatSetValues(B,1,&i,1,&i,&b,INSERT_VALUES));
89 : }
90 13 : PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
91 13 : PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
92 13 : PetscCall(MatSetOption(B,MAT_HERMITIAN,PETSC_TRUE));
93 :
94 : /* Functions: f1=-lambda, f2=1.0, f3=exp(-tau*lambda) */
95 13 : PetscCall(FNCreate(PETSC_COMM_WORLD,&f1));
96 13 : PetscCall(FNSetType(f1,FNRATIONAL));
97 13 : coeffs[0] = -1.0; coeffs[1] = 0.0;
98 13 : PetscCall(FNRationalSetNumerator(f1,2,coeffs));
99 :
100 13 : PetscCall(FNCreate(PETSC_COMM_WORLD,&f2));
101 13 : PetscCall(FNSetType(f2,FNRATIONAL));
102 13 : coeffs[0] = 1.0;
103 13 : PetscCall(FNRationalSetNumerator(f2,1,coeffs));
104 :
105 13 : PetscCall(FNCreate(PETSC_COMM_WORLD,&f3));
106 13 : PetscCall(FNSetType(f3,FNEXP));
107 13 : PetscCall(FNSetScale(f3,-tau,1.0));
108 :
109 : /* Set the split operator */
110 13 : mats[0] = A; funs[0] = f2;
111 13 : mats[1] = Id; funs[1] = f1;
112 13 : mats[2] = B; funs[2] = f3;
113 13 : PetscCall(NEPSetSplitOperator(nep,3,mats,funs,SUBSET_NONZERO_PATTERN));
114 :
115 : /* Customize nonlinear solver; set runtime options */
116 13 : PetscCall(NEPSetType(nep,NEPCISS));
117 13 : PetscCall(NEPSetDimensions(nep,1,24,PETSC_DETERMINE));
118 13 : PetscCall(NEPSetTolerances(nep,1e-9,PETSC_CURRENT));
119 13 : PetscCall(NEPGetRG(nep,&rg));
120 13 : PetscCall(RGSetType(rg,RGELLIPSE));
121 13 : PetscCall(RGEllipseSetParameters(rg,10.0,9.5,0.1));
122 13 : PetscCall(NEPCISSSetSizes(nep,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_DETERMINE,1,PETSC_DETERMINE,PETSC_TRUE));
123 13 : PetscCall(NEPCISSGetKSPs(nep,&nsolve,&ksp));
124 221 : for (i=0;i<nsolve;i++) {
125 208 : PetscCall(KSPSetTolerances(ksp[i],1e-12,PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
126 208 : PetscCall(KSPSetType(ksp[i],KSPBCGS));
127 208 : PetscCall(KSPGetPC(ksp[i],&pc));
128 208 : PetscCall(PCSetType(pc,PCSOR));
129 : }
130 13 : PetscCall(NEPSetFromOptions(nep));
131 :
132 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
133 : Solve the eigensystem
134 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
135 :
136 13 : PetscCall(PetscObjectTypeCompare((PetscObject)nep,NEPCISS,&flg));
137 13 : if (flg) {
138 13 : PetscCall(NEPCISSGetExtraction(nep,&ext));
139 13 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Running CISS with %" PetscInt_FMT " KSP solvers (%s extraction)\n",nsolve,NEPCISSExtractions[ext]));
140 : }
141 13 : PetscCall(NEPSolve(nep));
142 :
143 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144 : Display solution and clean up
145 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
146 :
147 : /* show detailed info unless -terse option is given by user */
148 13 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
149 13 : if (terse) PetscCall(NEPErrorView(nep,NEP_ERROR_RELATIVE,NULL));
150 : else {
151 0 : PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
152 0 : PetscCall(NEPConvergedReasonView(nep,PETSC_VIEWER_STDOUT_WORLD));
153 0 : PetscCall(NEPErrorView(nep,NEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
154 0 : PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
155 : }
156 13 : PetscCall(NEPDestroy(&nep));
157 13 : PetscCall(MatDestroy(&Id));
158 13 : PetscCall(MatDestroy(&A));
159 13 : PetscCall(MatDestroy(&B));
160 13 : PetscCall(FNDestroy(&f1));
161 13 : PetscCall(FNDestroy(&f2));
162 13 : PetscCall(FNDestroy(&f3));
163 13 : PetscCall(SlepcFinalize());
164 : return 0;
165 : }
166 :
167 : /*TEST
168 :
169 : build:
170 : requires: complex
171 :
172 : testset:
173 : args: -nep_ciss_extraction {{ritz hankel caa}} -terse
174 : requires: complex !single
175 : output_file: output/test11_1.out
176 : filter: sed -e "s/([A-Z]* extraction)//"
177 : test:
178 : suffix: 1
179 : args: -nep_ciss_delta 1e-10
180 : test:
181 : suffix: 2
182 : nsize: 2
183 : args: -nep_ciss_partitions 2
184 : test:
185 : suffix: 3
186 : args: -nep_ciss_moments 6 -nep_ciss_blocksize 4 -nep_ciss_refine_inner 1 -nep_ciss_refine_blocksize 2
187 :
188 : test:
189 : suffix: 4
190 : args: -terse -nep_view
191 : requires: complex !single
192 : filter: grep -v tolerance | grep -v threshold
193 :
194 : TEST*/
|
