Actual source code: gun.c
slepc-3.20.1 2023-11-27
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, 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 https://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: int main(int argc,char **argv)
35: {
36: Mat A[NMAT]; /* problem matrices */
37: FN f[NMAT]; /* functions to define the nonlinear operator */
38: FN ff[2]; /* auxiliary functions to define the nonlinear operator */
39: NEP nep; /* nonlinear eigensolver context */
40: PetscBool terse,flg;
41: const char* string[NMAT]={"-K","-M","-W1","-W2"};
42: char filename[PETSC_MAX_PATH_LEN];
43: PetscScalar numer[2],sigma;
44: PetscInt i;
45: PetscViewer viewer;
47: PetscFunctionBeginUser;
48: PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
50: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"GUN problem\n\n"));
51: #if !defined(PETSC_USE_COMPLEX)
52: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This example requires complex scalars!");
53: #endif
55: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
56: Load the problem matrices
57: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
59: for (i=0;i<NMAT;i++) {
60: PetscCall(PetscOptionsGetString(NULL,NULL,string[i],filename,sizeof(filename),&flg));
61: PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must indicate a filename with the %s option",string[i]);
62: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
63: PetscCall(MatCreate(PETSC_COMM_WORLD,&A[i]));
64: PetscCall(MatSetFromOptions(A[i]));
65: PetscCall(MatLoad(A[i],viewer));
66: PetscCall(PetscViewerDestroy(&viewer));
67: }
69: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
70: Create the problem functions
71: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
73: /* f1=1 */
74: PetscCall(FNCreate(PETSC_COMM_WORLD,&f[0]));
75: PetscCall(FNSetType(f[0],FNRATIONAL));
76: numer[0] = 1.0;
77: PetscCall(FNRationalSetNumerator(f[0],1,numer));
79: /* f2=-lambda */
80: PetscCall(FNCreate(PETSC_COMM_WORLD,&f[1]));
81: PetscCall(FNSetType(f[1],FNRATIONAL));
82: numer[0] = -1.0; numer[1] = 0.0;
83: PetscCall(FNRationalSetNumerator(f[1],2,numer));
85: /* f3=i*sqrt(lambda) */
86: PetscCall(FNCreate(PETSC_COMM_WORLD,&f[2]));
87: PetscCall(FNSetType(f[2],FNSQRT));
88: PetscCall(FNSetScale(f[2],1.0,PETSC_i));
90: /* f4=i*sqrt(lambda-sigma^2) */
91: sigma = SIGMA*SIGMA;
92: PetscCall(FNCreate(PETSC_COMM_WORLD,&ff[0]));
93: PetscCall(FNSetType(ff[0],FNSQRT));
94: PetscCall(FNCreate(PETSC_COMM_WORLD,&ff[1]));
95: PetscCall(FNSetType(ff[1],FNRATIONAL));
96: numer[0] = 1.0; numer[1] = -sigma;
97: PetscCall(FNRationalSetNumerator(ff[1],2,numer));
98: PetscCall(FNCreate(PETSC_COMM_WORLD,&f[3]));
99: PetscCall(FNSetType(f[3],FNCOMBINE));
100: PetscCall(FNCombineSetChildren(f[3],FN_COMBINE_COMPOSE,ff[1],ff[0]));
101: PetscCall(FNSetScale(f[3],1.0,PETSC_i));
103: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104: Create the eigensolver and solve the problem
105: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
107: PetscCall(NEPCreate(PETSC_COMM_WORLD,&nep));
108: PetscCall(NEPSetSplitOperator(nep,4,A,f,UNKNOWN_NONZERO_PATTERN));
109: PetscCall(NEPSetFromOptions(nep));
111: PetscCall(NEPSolve(nep));
113: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: Display solution and clean up
115: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
117: /* show detailed info unless -terse option is given by user */
118: PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
119: if (terse) PetscCall(NEPErrorView(nep,NEP_ERROR_RELATIVE,NULL));
120: else {
121: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
122: PetscCall(NEPConvergedReasonView(nep,PETSC_VIEWER_STDOUT_WORLD));
123: PetscCall(NEPErrorView(nep,NEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
124: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
125: }
126: PetscCall(NEPDestroy(&nep));
127: for (i=0;i<NMAT;i++) {
128: PetscCall(MatDestroy(&A[i]));
129: PetscCall(FNDestroy(&f[i]));
130: }
131: for (i=0;i<2;i++) PetscCall(FNDestroy(&ff[i]));
132: PetscCall(SlepcFinalize());
133: return 0;
134: }
136: /*TEST
138: build:
139: requires: complex
141: test:
142: suffix: 1
143: args: -K ${DATAFILESPATH}/matrices/complex/gun_K.petsc -M ${DATAFILESPATH}/matrices/complex/gun_M.petsc -W1 ${DATAFILESPATH}/matrices/complex/gun_W1.petsc -W2 ${DATAFILESPATH}/matrices/complex/gun_W2.petsc -nep_type nleigs -rg_type polygon -rg_polygon_vertices 12500-1i,120500-1i,120500+30000i,70000+30000i -nep_target 65000 -nep_nev 24 -terse
144: requires: double complex datafilespath !defined(PETSC_USE_64BIT_INDICES) !defined(PETSCTEST_VALGRIND)
145: timeoutfactor: 10
147: TEST*/