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 : This example implements one of the problems found at
12 : NLEVP: A Collection of Nonlinear Eigenvalue Problems,
13 : The University of Manchester.
14 : The details of the collection can be found at:
15 : [1] T. Betcke et al., "NLEVP: A Collection of Nonlinear Eigenvalue
16 : Problems", ACM Trans. Math. Software 39(2), Article 7, 2013.
17 :
18 : The planar_waveguide problem is a quartic PEP with symmetric matrices,
19 : arising from a finite element solution of the propagation constants in a
20 : six-layer planar waveguide.
21 : */
22 :
23 : static char help[] = "FEM solution of the propagation constants in a six-layer planar waveguide.\n\n"
24 : "The command line options are:\n"
25 : " -n <n>, the dimension of the matrices.\n\n";
26 :
27 : #include <slepcpep.h>
28 :
29 : #define NMAT 5
30 : #define NL 6
31 :
32 2 : int main(int argc,char **argv)
33 : {
34 2 : Mat A[NMAT]; /* problem matrices */
35 2 : PEP pep; /* polynomial eigenproblem solver context */
36 2 : PetscInt n=128,nlocal,k,Istart,Iend,i,j,start_ct,end_ct;
37 2 : PetscReal w=9.92918,a=0.0,b=2.0,h,deltasq;
38 2 : PetscReal nref[NL],K2[NL],q[NL],*md,*supd,*subd;
39 2 : PetscScalar v,alpha;
40 2 : PetscBool terse;
41 :
42 2 : PetscFunctionBeginUser;
43 2 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
44 :
45 2 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
46 2 : n = (n/4)*4;
47 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nPlanar waveguide, n=%" PetscInt_FMT "\n\n",n+1));
48 2 : h = (b-a)/n;
49 2 : nlocal = (n/4)-1;
50 :
51 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
52 : Set waveguide parameters used in construction of matrices
53 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
54 :
55 : /* refractive indices in each layer */
56 2 : nref[0] = 1.5;
57 2 : nref[1] = 1.66;
58 2 : nref[2] = 1.6;
59 2 : nref[3] = 1.53;
60 2 : nref[4] = 1.66;
61 2 : nref[5] = 1.0;
62 :
63 14 : for (i=0;i<NL;i++) K2[i] = w*w*nref[i]*nref[i];
64 2 : deltasq = K2[0] - K2[NL-1];
65 14 : for (i=0;i<NL;i++) q[i] = K2[i] - (K2[0] + K2[NL-1])/2;
66 :
67 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
68 : Compute the polynomial matrices
69 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
70 :
71 : /* initialize matrices */
72 12 : for (i=0;i<NMAT;i++) {
73 10 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A[i]));
74 10 : PetscCall(MatSetSizes(A[i],PETSC_DECIDE,PETSC_DECIDE,n+1,n+1));
75 10 : PetscCall(MatSetFromOptions(A[i]));
76 : }
77 2 : PetscCall(MatGetOwnershipRange(A[0],&Istart,&Iend));
78 :
79 : /* A0 */
80 2 : alpha = (h/6)*(deltasq*deltasq/16);
81 260 : for (i=Istart;i<Iend;i++) {
82 258 : v = 4.0;
83 258 : if (i==0 || i==n) v = 2.0;
84 258 : PetscCall(MatSetValue(A[0],i,i,v*alpha,INSERT_VALUES));
85 258 : if (i>0) PetscCall(MatSetValue(A[0],i,i-1,alpha,INSERT_VALUES));
86 258 : if (i<n) PetscCall(MatSetValue(A[0],i,i+1,alpha,INSERT_VALUES));
87 : }
88 :
89 : /* A1 */
90 2 : if (Istart==0) PetscCall(MatSetValue(A[1],0,0,-deltasq/4,INSERT_VALUES));
91 2 : if (Iend==n+1) PetscCall(MatSetValue(A[1],n,n,deltasq/4,INSERT_VALUES));
92 :
93 : /* A2 */
94 2 : alpha = 1.0/h;
95 260 : for (i=Istart;i<Iend;i++) {
96 258 : v = 2.0;
97 258 : if (i==0 || i==n) v = 1.0;
98 258 : PetscCall(MatSetValue(A[2],i,i,v*alpha,ADD_VALUES));
99 258 : if (i>0) PetscCall(MatSetValue(A[2],i,i-1,-alpha,ADD_VALUES));
100 258 : if (i<n) PetscCall(MatSetValue(A[2],i,i+1,-alpha,ADD_VALUES));
101 : }
102 2 : PetscCall(PetscMalloc3(n+1,&md,n+1,&supd,n+1,&subd));
103 :
104 2 : md[0] = 2.0*q[1];
105 2 : supd[1] = q[1];
106 2 : subd[0] = q[1];
107 :
108 10 : for (k=1;k<=NL-2;k++) {
109 :
110 8 : end_ct = k*(nlocal+1);
111 8 : start_ct = end_ct-nlocal;
112 :
113 256 : for (j=start_ct;j<end_ct;j++) {
114 248 : md[j] = 4*q[k];
115 248 : supd[j+1] = q[k];
116 248 : subd[j] = q[k];
117 : }
118 :
119 8 : if (k < 4) { /* interface points */
120 6 : md[end_ct] = 4*(q[k] + q[k+1])/2.0;
121 6 : supd[end_ct+1] = q[k+1];
122 6 : subd[end_ct] = q[k+1];
123 : }
124 :
125 : }
126 :
127 2 : md[n] = 2*q[NL-2];
128 2 : supd[n] = q[NL-2];
129 2 : subd[n] = q[NL-2];
130 :
131 2 : alpha = -h/6.0;
132 260 : for (i=Istart;i<Iend;i++) {
133 258 : PetscCall(MatSetValue(A[2],i,i,md[i]*alpha,ADD_VALUES));
134 258 : if (i>0) PetscCall(MatSetValue(A[2],i,i-1,subd[i-1]*alpha,ADD_VALUES));
135 258 : if (i<n) PetscCall(MatSetValue(A[2],i,i+1,supd[i+1]*alpha,ADD_VALUES));
136 : }
137 2 : PetscCall(PetscFree3(md,supd,subd));
138 :
139 : /* A3 */
140 2 : if (Istart==0) PetscCall(MatSetValue(A[3],0,0,1.0,INSERT_VALUES));
141 2 : if (Iend==n+1) PetscCall(MatSetValue(A[3],n,n,1.0,INSERT_VALUES));
142 :
143 : /* A4 */
144 2 : alpha = (h/6);
145 260 : for (i=Istart;i<Iend;i++) {
146 258 : v = 4.0;
147 258 : if (i==0 || i==n) v = 2.0;
148 258 : PetscCall(MatSetValue(A[4],i,i,v*alpha,INSERT_VALUES));
149 258 : if (i>0) PetscCall(MatSetValue(A[4],i,i-1,alpha,INSERT_VALUES));
150 258 : if (i<n) PetscCall(MatSetValue(A[4],i,i+1,alpha,INSERT_VALUES));
151 : }
152 :
153 : /* assemble matrices */
154 12 : for (i=0;i<NMAT;i++) PetscCall(MatAssemblyBegin(A[i],MAT_FINAL_ASSEMBLY));
155 12 : for (i=0;i<NMAT;i++) PetscCall(MatAssemblyEnd(A[i],MAT_FINAL_ASSEMBLY));
156 :
157 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158 : Create the eigensolver and solve the problem
159 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
160 :
161 2 : PetscCall(PEPCreate(PETSC_COMM_WORLD,&pep));
162 2 : PetscCall(PEPSetOperators(pep,NMAT,A));
163 2 : PetscCall(PEPSetFromOptions(pep));
164 2 : PetscCall(PEPSolve(pep));
165 :
166 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167 : Display solution and clean up
168 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
169 :
170 : /* show detailed info unless -terse option is given by user */
171 2 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
172 2 : if (terse) PetscCall(PEPErrorView(pep,PEP_ERROR_BACKWARD,NULL));
173 : else {
174 0 : PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
175 0 : PetscCall(PEPConvergedReasonView(pep,PETSC_VIEWER_STDOUT_WORLD));
176 0 : PetscCall(PEPErrorView(pep,PEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD));
177 0 : PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
178 : }
179 2 : PetscCall(PEPDestroy(&pep));
180 12 : for (i=0;i<NMAT;i++) PetscCall(MatDestroy(&A[i]));
181 2 : PetscCall(SlepcFinalize());
182 : return 0;
183 : }
184 :
185 : /*TEST
186 :
187 : test:
188 : suffix: 1
189 : args: -pep_type {{toar linear}} -pep_nev 4 -st_type sinvert -terse
190 : requires: !single
191 :
192 : TEST*/
|
