Actual source code: loaded_string.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, 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 loaded_string problem is a rational eigenvalue problem for the
 19:    finite element model of a loaded vibrating string.
 20:    This example solves the loaded_string problem by first transforming
 21:    it to a quadratic eigenvalue problem.
 22: */

 24: static char help[] = "Finite element model of a loaded vibrating string.\n\n"
 25:   "The command line options are:\n"
 26:   "  -n <n>, dimension of the matrices.\n"
 27:   "  -kappa <kappa>, stiffness of elastic spring.\n"
 28:   "  -mass <m>, mass of the attached load.\n\n";

 30: #include <slepcpep.h>

 32: #define NMAT 3

 34: int main(int argc,char **argv)
 35: {
 36:   Mat            A[3],M;      /* problem matrices */
 37:   PEP            pep;         /* polynomial eigenproblem solver context */
 38:   PetscInt       n=100,Istart,Iend,i;
 39:   PetscBool      terse;
 40:   PetscReal      kappa=1.0,m=1.0;
 41:   PetscScalar    sigma;

 44:   SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;

 46:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 47:   PetscOptionsGetReal(NULL,NULL,"-kappa",&kappa,NULL);
 48:   PetscOptionsGetReal(NULL,NULL,"-mass",&m,NULL);
 49:   sigma = kappa/m;
 50:   PetscPrintf(PETSC_COMM_WORLD,"Loaded vibrating string (QEP), n=%D kappa=%g m=%g\n\n",n,(double)kappa,(double)m);

 52:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 53:      Compute the matrices that define the eigensystem, (k^2*M+k*C+K)x=0
 54:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 55:   /* initialize matrices */
 56:   for (i=0;i<NMAT;i++) {
 57:     MatCreate(PETSC_COMM_WORLD,&A[i]);
 58:     MatSetSizes(A[i],PETSC_DECIDE,PETSC_DECIDE,n,n);
 59:     MatSetFromOptions(A[i]);
 60:     MatSetUp(A[i]);
 61:   }
 62:   MatGetOwnershipRange(A[0],&Istart,&Iend);

 64:   /* A0 */
 65:   for (i=Istart;i<Iend;i++) {
 66:     MatSetValue(A[0],i,i,(i==n-1)?1.0*n:2.0*n,INSERT_VALUES);
 67:     if (i>0) { MatSetValue(A[0],i,i-1,-1.0*n,INSERT_VALUES); }
 68:     if (i<n-1) { MatSetValue(A[0],i,i+1,-1.0*n,INSERT_VALUES); }
 69:   }

 71:   /* A1 */
 72:   for (i=Istart;i<Iend;i++) {
 73:     MatSetValue(A[1],i,i,(i==n-1)?2.0/(6.0*n):4.0/(6.0*n),INSERT_VALUES);
 74:     if (i>0) { MatSetValue(A[1],i,i-1,1.0/(6.0*n),INSERT_VALUES); }
 75:     if (i<n-1) { MatSetValue(A[1],i,i+1,1.0/(6.0*n),INSERT_VALUES); }
 76:   }

 78:   /* A2 */
 79:   if (Istart<=n-1 && n-1<Iend) {
 80:     MatSetValue(A[2],n-1,n-1,kappa,INSERT_VALUES); 
 81:   }

 83:   /* assemble matrices */
 84:   for (i=0;i<NMAT;i++) {
 85:     MatAssemblyBegin(A[i],MAT_FINAL_ASSEMBLY);
 86:   }
 87:   for (i=0;i<NMAT;i++) {
 88:     MatAssemblyEnd(A[i],MAT_FINAL_ASSEMBLY);
 89:   }

 91:   /* build matrices for the QEP */
 92:   MatAXPY(A[2],1.0,A[0],DIFFERENT_NONZERO_PATTERN);
 93:   MatAXPY(A[2],sigma,A[1],SAME_NONZERO_PATTERN);
 94:   MatScale(A[2],-1.0);
 95:   MatScale(A[0],sigma);
 96:   M = A[1];
 97:   A[1] = A[2];
 98:   A[2] = M;

100:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101:                 Create the eigensolver and solve the problem
102:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

104:   PEPCreate(PETSC_COMM_WORLD,&pep);
105:   PEPSetOperators(pep,3,A);
106:   PEPSetFromOptions(pep);
107:   PEPSolve(pep);

109:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110:                     Display solution and clean up
111:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

113:   /* show detailed info unless -terse option is given by user */
114:   PetscOptionsHasName(NULL,NULL,"-terse",&terse);
115:   if (terse) {
116:     PEPErrorView(pep,PEP_ERROR_BACKWARD,NULL);
117:   } else {
118:     PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
119:     PEPReasonView(pep,PETSC_VIEWER_STDOUT_WORLD);
120:     PEPErrorView(pep,PEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD);
121:     PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
122:   }
123:   PEPDestroy(&pep);
124:   for (i=0;i<NMAT;i++) {
125:     MatDestroy(&A[i]);
126:   }
127:   SlepcFinalize();
128:   return ierr;
129: }

131: /*TEST

133:    test:
134:       suffix: 1
135:       args: -pep_hyperbolic -pep_interval 4,900 -pep_type stoar -st_type sinvert -st_pc_type cholesky -terse
136:       requires: !single

138: TEST*/