Actual source code: gun.c
slepc-3.8.3 2018-04-03
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2017, 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: */
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.
18: The gun problem arises from model of a radio-frequency gun cavity, with
19: the complex nonlinear function
20: T(lambda) = K-lambda*M+i*lambda^(1/2)*W1+i*(lambda-108.8774^2)^(1/2)*W2
22: Data files can be downloaded from http://slepc.upv.es/datafiles
23: */
25: static char help[] = "Radio-frequency gun cavity.\n\n"
26: "The command line options are:\n"
27: "-K <filename1> -M <filename2> -W1 <filename3> -W2 <filename4>, where filename1,..,filename4 are files containing the matrices in PETSc binary form defining the GUN problem.\n\n";
29: #include <slepcnep.h>
31: #define NMAT 4
32: #define SIGMA 108.8774
34: PetscErrorCode ComputeSingularities(NEP,PetscInt*,PetscScalar*,void*);
36: int main(int argc,char **argv)
37: {
39: Mat A[NMAT]; /* problem matrices */
40: FN f[NMAT]; /* functions to define the nonlinear operator */
41: FN ff[2]; /* auxiliary functions to define the nonlinear operator */
42: NEP nep; /* nonlinear eigensolver context */
43: PetscBool terse,flg;
44: const char* string[NMAT]={"-K","-M","-W1","-W2"};
45: char filename[PETSC_MAX_PATH_LEN];
46: PetscScalar numer[2],sigma;
47: PetscInt i;
48: PetscViewer viewer;
50: SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
52: PetscPrintf(PETSC_COMM_WORLD,"GUN problem\n\n");
53: #if !defined(PETSC_USE_COMPLEX)
54: SETERRQ(PETSC_COMM_WORLD,1,"This example requires complex scalars!");
55: #endif
57: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
58: Load the problem matrices
59: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
61: for (i=0;i<NMAT;i++) {
62: PetscOptionsGetString(NULL,NULL,string[i],filename,PETSC_MAX_PATH_LEN,&flg);
63: if (!flg) SETERRQ1(PETSC_COMM_WORLD,1,"Must indicate a filename with the %s option",string[i]);
64: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
65: MatCreate(PETSC_COMM_WORLD,&A[i]);
66: MatSetFromOptions(A[i]);
67: MatLoad(A[i],viewer);
68: PetscViewerDestroy(&viewer);
69: }
71: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
72: Create the problem functions
73: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
75: /* f1=1 */
76: FNCreate(PETSC_COMM_WORLD,&f[0]);
77: FNSetType(f[0],FNRATIONAL);
78: numer[0] = 1.0;
79: FNRationalSetNumerator(f[0],1,numer);
81: /* f2=-lambda */
82: FNCreate(PETSC_COMM_WORLD,&f[1]);
83: FNSetType(f[1],FNRATIONAL);
84: numer[0] = -1.0; numer[1] = 0.0;
85: FNRationalSetNumerator(f[1],2,numer);
87: /* f3=i*sqrt(lambda) */
88: FNCreate(PETSC_COMM_WORLD,&f[2]);
89: FNSetType(f[2],FNSQRT);
90: FNSetScale(f[2],1.0,PETSC_i);
92: /* f4=i*sqrt(lambda-sigma^2) */
93: sigma = SIGMA*SIGMA;
94: FNCreate(PETSC_COMM_WORLD,&ff[0]);
95: FNSetType(ff[0],FNSQRT);
96: FNCreate(PETSC_COMM_WORLD,&ff[1]);
97: FNSetType(ff[1],FNRATIONAL);
98: numer[0] = 1.0; numer[1] = -sigma;
99: FNRationalSetNumerator(ff[1],2,numer);
100: FNCreate(PETSC_COMM_WORLD,&f[3]);
101: FNSetType(f[3],FNCOMBINE);
102: FNCombineSetChildren(f[3],FN_COMBINE_COMPOSE,ff[1],ff[0]);
103: FNSetScale(f[3],1.0,PETSC_i);
105: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106: Create the eigensolver and solve the problem
107: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109: NEPCreate(PETSC_COMM_WORLD,&nep);
110: NEPSetSplitOperator(nep,4,A,f,DIFFERENT_NONZERO_PATTERN);
111: NEPSetFromOptions(nep);
113: PetscObjectTypeCompare((PetscObject)nep,NEPNLEIGS,&flg);
114: if (flg) {
115: NEPNLEIGSSetSingularitiesFunction(nep,ComputeSingularities,NULL);
116: }
118: NEPSolve(nep);
120: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121: Display solution and clean up
122: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
124: /* show detailed info unless -terse option is given by user */
125: PetscOptionsHasName(NULL,NULL,"-terse",&terse);
126: if (terse) {
127: NEPErrorView(nep,NEP_ERROR_RELATIVE,NULL);
128: } else {
129: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
130: NEPReasonView(nep,PETSC_VIEWER_STDOUT_WORLD);
131: NEPErrorView(nep,NEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);
132: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
133: }
134: NEPDestroy(&nep);
135: for (i=0;i<NMAT;i++) {
136: MatDestroy(&A[i]);
137: FNDestroy(&f[i]);
138: }
139: for (i=0;i<2;i++) {
140: FNDestroy(&ff[i]);
141: }
142: SlepcFinalize();
143: return ierr;
144: }
146: /*
147: ComputeSingularities - Computes maxnp points (at most) in the complex plane where
148: the function T(.) is not analytic.
150: In this case, we discretize the singularity region (-inf,108.8774^2)~(-10e+12,-10e-12+108.8774^2)
151: */
152: PetscErrorCode ComputeSingularities(NEP nep,PetscInt *maxnp,PetscScalar *xi,void *pt)
153: {
154: PetscReal h;
155: PetscInt i;
156: PetscReal sigma,end;
159: sigma = SIGMA*SIGMA;
160: end = PetscLogReal(sigma);
161: h = (12.0+end)/(*maxnp-1);
162: xi[0] = sigma;
163: for (i=1;i<*maxnp;i++) xi[i] = -PetscPowReal(10,h*i)+sigma;
164: return(0);
165: }