Actual source code: ex34.c

slepc-3.8.3 2018-04-03
Report Typos and Errors
  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 is a nonlinear eigenvalue problem. When p=2, it is reduced to a linear Laplace eigenvalue
 12:    problem.

 14:    -\nabla\cdot(|\nabla u|^{p-2} \nabla u) = k |u|^{p-2} u in (0,1)x(0,1),

 16:    u = 0 on the entire boundary.

 18:    The code is implemented based on DMPlex using Q1 FEM on a quadrilateral mesh. In this code, we consider p=3.

 20:    Contributed  by Fande Kong fdkong.jd@gmail.com
 21: */

 23: static char help[] = "Nonlinear inverse iteration for A(x)*x=lambda*B(x)*x.\n\n";


 26: #include <slepceps.h>
 27: #include <petscdmplex.h>
 28: #include <petscds.h>

 30: PetscErrorCode CreateSquareMesh(MPI_Comm,DM*);
 31: PetscErrorCode SetupDiscretization(DM);
 32: PetscErrorCode FormJacobianA(SNES,Vec,Mat,Mat,void*);
 33: PetscErrorCode FormFunctionA(SNES,Vec,Vec,void*);
 34: PetscErrorCode FormJacobianB(SNES,Vec,Mat,Mat,void*);
 35: PetscErrorCode FormFunctionB(SNES,Vec,Vec,void*);
 36: PetscErrorCode BoundaryGlobalIndex(DM,const char*,IS*);

 38: typedef struct {
 39:   IS    bdis; /* global indices for boundary DoFs */
 40: } AppCtx;

 42: int main(int argc,char **argv)
 43: {
 44:   DM             dm;
 45:   MPI_Comm       comm;
 46:   AppCtx         user;
 47:   EPS            eps;  /* eigenproblem solver context */
 48:   EPSType        type;
 49:   Mat            A,B;
 50:   PetscContainer container;
 51:   PetscInt       nev,nconv;
 52:   PetscBool      nonlin;
 53:   SNES           snes;

 56:   SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 57:   comm = PETSC_COMM_WORLD;
 58:   /* Create a quadrilateral mesh on domain (0,1)x(0,1) */
 59:   CreateSquareMesh(comm,&dm);
 60:   /* Setup basis function */
 61:   SetupDiscretization(dm);
 62:   BoundaryGlobalIndex(dm,"marker",&user.bdis);

 64:   DMCreateMatrix(dm,&A);
 65:   MatDuplicate(A,MAT_COPY_VALUES,&B);

 67:   /*
 68:      Compose callback functions and context that will be needed by the solver
 69:   */
 70:   PetscObjectComposeFunction((PetscObject)A,"formFunction",FormFunctionA);
 71:   PetscObjectComposeFunction((PetscObject)A,"formJacobian",FormJacobianA);
 72:   PetscObjectComposeFunction((PetscObject)B,"formFunction",FormFunctionB);
 73:   PetscContainerCreate(comm,&container);
 74:   PetscContainerSetPointer(container,&user);
 75:   PetscObjectCompose((PetscObject)A,"formFunctionCtx",(PetscObject)container);
 76:   PetscObjectCompose((PetscObject)A,"formJacobianCtx",(PetscObject)container);
 77:   PetscObjectCompose((PetscObject)B,"formFunctionCtx",(PetscObject)container);
 78:   PetscContainerDestroy(&container);

 80:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 81:                 Create the eigensolver and set various options
 82:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 84:   EPSCreate(comm,&eps);
 85:   EPSSetOperators(eps,A,B);
 86:   EPSSetProblemType(eps,EPS_GNHEP);
 87:   /*
 88:      Use nonlinear inverse iteration
 89:   */
 90:   EPSSetType(eps,EPSPOWER);
 91:   EPSPowerSetNonlinear(eps,PETSC_TRUE);
 92:   /*
 93:     Attach DM to SNES
 94:   */
 95:   EPSPowerGetSNES(eps,&snes);
 96:   SNESSetDM(snes,dm);
 97:   EPSSetFromOptions(eps);

 99:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100:                       Solve the eigensystem
101:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

103:   EPSSolve(eps);

105:   /*
106:      Optional: Get some information from the solver and display it
107:   */
108:   EPSGetType(eps,&type);
109:   EPSPowerGetNonlinear(eps,&nonlin);
110:   PetscPrintf(comm," Solution method: %s%s\n\n",type,nonlin?" (nonlinear)":"");
111:   EPSGetDimensions(eps,&nev,NULL,NULL);
112:   PetscPrintf(comm," Number of requested eigenvalues: %D\n",nev);

114:   /* print eigenvalue and error */
115:   EPSGetConverged(eps,&nconv);
116:   if (nconv>0) {
117:     PetscScalar   k;
118:     PetscReal     na,nb;
119:     Vec           a,b,eigen;
120:     DMCreateGlobalVector(dm,&a);
121:     VecDuplicate(a,&b);
122:     VecDuplicate(a,&eigen);
123:     EPSGetEigenpair(eps,0,&k,NULL,eigen,NULL);
124:     FormFunctionA(snes,eigen,a,&user);
125:     FormFunctionB(snes,eigen,b,&user);
126:     VecAXPY(a,-k,b);
127:     VecNorm(a,NORM_2,&na);
128:     VecNorm(b,NORM_2,&nb);
129:     PetscPrintf(comm,"k: %g error: %g\n",(double)PetscRealPart(k),(double)na/nb);
130:     VecDestroy(&a);
131:     VecDestroy(&b);
132:     VecDestroy(&eigen);
133:   } else {
134:     PetscPrintf(comm,"Solver did not converge\n");
135:   }

137:   MatDestroy(&A);
138:   MatDestroy(&B);
139:   DMDestroy(&dm);
140:   EPSDestroy(&eps);
141:   ISDestroy(&user.bdis);
142:   SlepcFinalize();
143:   return ierr;
144: }

146: /* <|u|u, v> */
147: static void f0_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
148:                  const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
149:                  const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
150:                  PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
151: {
152:   PetscScalar cof = PetscAbsScalar(u[0]);

154:   f0[0] = cof*u[0];
155: }

157: /* <|\nabla u| \nabla u, \nabla v> */
158: static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
159:                  const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
160:                  const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
161:                  PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
162: {
163:   PetscInt    d;
164:   PetscScalar cof = 0;
165:   for (d = 0; d < dim; ++d)  cof += u_x[d]*u_x[d];

167:   cof = PetscSqrtScalar(cof);

169:   for (d = 0; d < dim; ++d) f1[d] = u_x[d]*cof;
170: }

172: /* approximate  Jacobian for   <|u|u, v> */
173: static void g0_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
174:                   const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
175:                   const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
176:                   PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
177: {
178:   g0[0] = 1.0*PetscAbsScalar(u[0]);
179: }

181: /* approximate  Jacobian for   <|\nabla u| \nabla u, \nabla v> */
182: static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
183:                   const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
184:                   const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
185:                   PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
186: {
187:   PetscInt d;
188:   for (d = 0; d < dim; ++d) g3[d*dim+d] = 1.0;
189: }

191: PetscErrorCode SetupDiscretization(DM dm)
192: {
193:   PetscFE        fe;
194:   PetscDS        prob;
195:   PetscInt       totDim;

199:   /* Create finite element */
200:   PetscFECreateDefault(dm,2,1,PETSC_FALSE,NULL,-1,&fe);
201:   PetscObjectSetName((PetscObject) fe, "u");
202:   DMGetDS(dm,&prob);
203:   PetscDSSetDiscretization(prob,0,(PetscObject) fe);
204:   DMSetDS(dm,prob);
205:   PetscDSGetTotalDimension(prob,&totDim);
206:   PetscFEDestroy(&fe);
207:   return(0);
208: }

210: PetscErrorCode CreateSquareMesh(MPI_Comm comm,DM *dm)
211: {
212:   PetscInt       cells[] = {8,8};
213:   PetscInt       dim = 2;
214:   DM             pdm;
215:   PetscMPIInt    size;

219:   DMPlexCreateHexBoxMesh(comm,dim,cells,DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,dm);
220:   DMSetFromOptions(*dm);
221:   DMSetUp(*dm);
222:   MPI_Comm_size(comm,&size);
223:   if (size > 1) {
224:     DMPlexDistribute(*dm,0,NULL,&pdm);
225:     DMDestroy(dm);
226:     *dm = pdm;
227:   }
228:   return(0);
229: }

231: PetscErrorCode BoundaryGlobalIndex(DM dm,const char labelname[],IS *bdis)
232: {
233:   IS             bdpoints;
234:   PetscInt       nindices,*indices,numDof,offset,npoints,i,j;
235:   const PetscInt *bdpoints_indices;
236:   DMLabel        bdmarker;
237:   PetscSection   gsection;

241:   DMGetDefaultGlobalSection(dm,&gsection);
242:   DMGetLabel(dm,labelname,&bdmarker);
243:   DMLabelGetStratumIS(bdmarker,1,&bdpoints);
244:   ISGetLocalSize(bdpoints,&npoints);
245:   ISGetIndices(bdpoints,&bdpoints_indices);
246:   nindices = 0;
247:   for (i=0;i<npoints;i++) {
248:     PetscSectionGetDof(gsection,bdpoints_indices[i],&numDof);
249:     if (numDof<=0) continue;
250:     nindices += numDof;
251:   }
252:   PetscCalloc1(nindices,&indices);
253:   nindices = 0;
254:   for (i=0;i<npoints;i++) {
255:     PetscSectionGetDof(gsection,bdpoints_indices[i],&numDof);
256:     if (numDof<=0) continue;
257:     PetscSectionGetOffset(gsection,bdpoints_indices[i],&offset);
258:     for (j=0;j<numDof;j++) indices[nindices++] = offset+j;
259:   }
260:   ISRestoreIndices(bdpoints,&bdpoints_indices);
261:   ISDestroy(&bdpoints);
262:   ISCreateGeneral(PetscObjectComm((PetscObject)dm),nindices,indices,PETSC_OWN_POINTER,bdis);
263:   return(0);
264: }

266: static PetscErrorCode FormJacobian(SNES snes,Vec X,Mat A,Mat B,void *ctx)
267: {
268:   DM             dm;
269:   Vec            Xloc;

273:   SNESGetDM(snes,&dm);
274:   DMGetLocalVector(dm,&Xloc);
275:   VecZeroEntries(Xloc);
276:   DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
277:   DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
278:   CHKMEMQ;
279:   DMPlexSNESComputeJacobianFEM(dm,Xloc,A,B,ctx);
280:   CHKMEMQ;
281:   DMRestoreLocalVector(dm,&Xloc);
282:   if (A != B) {
283:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
284:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
285:   }
286:   return(0);
287: }

289: PetscErrorCode FormJacobianA(SNES snes,Vec X,Mat A,Mat B,void *ctx)
290: {
292:   DM             dm;
293:   PetscDS        prob;
294:   AppCtx         *userctx = (AppCtx *)ctx;

297:   MatSetOption(B,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE);
298:   SNESGetDM(snes,&dm);
299:   DMGetDS(dm,&prob);
300:   PetscDSSetJacobian(prob,0,0,NULL,NULL,NULL,g3_uu);
301:   FormJacobian(snes,X,A,B,ctx);
302:   MatZeroRowsIS(B,userctx->bdis,1.0,NULL,NULL);
303:   return(0);
304: }

306: PetscErrorCode FormJacobianB(SNES snes,Vec X,Mat A,Mat B,void *ctx)
307: {
309:   DM             dm;
310:   PetscDS        prob;
311:   AppCtx         *userctx = (AppCtx *)ctx;

314:   MatSetOption(B,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE);
315:   SNESGetDM(snes,&dm);
316:   DMGetDS(dm,&prob);
317:   PetscDSSetJacobian(prob,0,0,g0_uu,NULL,NULL,NULL);
318:   FormJacobian(snes,X,A,B,ctx);
319:   MatZeroRowsIS(B,userctx->bdis,0.0,NULL,NULL);
320:   return(0);
321: }

323: static PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ctx)
324: {
325:   DM             dm;
326:   Vec            Xloc,Floc;

330:   SNESGetDM(snes,&dm);
331:   DMGetLocalVector(dm,&Xloc);
332:   DMGetLocalVector(dm,&Floc);
333:   VecZeroEntries(Xloc);
334:   VecZeroEntries(Floc);
335:   DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
336:   DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
337:   CHKMEMQ;
338:   DMPlexSNESComputeResidualFEM(dm,Xloc,Floc,ctx);
339:   CHKMEMQ;
340:   VecZeroEntries(F);
341:   DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
342:   DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
343:   DMRestoreLocalVector(dm,&Xloc);
344:   DMRestoreLocalVector(dm,&Floc);
345:   return(0);
346: }

348: PetscErrorCode FormFunctionA(SNES snes,Vec X,Vec F,void *ctx)
349: {
351:   DM             dm;
352:   PetscDS        prob;
353:   PetscInt       nindices,iStart,iEnd,i;
354:   AppCtx         *userctx = (AppCtx *)ctx;
355:   PetscScalar    *array,value;
356:   const PetscInt *indices;
357:   PetscInt       vecstate;

360:   SNESGetDM(snes,&dm);
361:   DMGetDS(dm,&prob);
362:   /* hook functions */
363:   PetscDSSetResidual(prob,0,NULL,f1_u);
364:   FormFunction(snes,X,F,ctx);
365:   /* Boundary condition */
366:   VecLockGet(X,&vecstate);
367:   if (vecstate>0) {
368:     VecLockPop(X);
369:   }
370:   VecGetOwnershipRange(X,&iStart,&iEnd);
371:   VecGetArray(X,&array);
372:   ISGetLocalSize(userctx->bdis,&nindices);
373:   ISGetIndices(userctx->bdis,&indices);
374:   for (i=0;i<nindices;i++) {
375:     value = array[indices[i]-iStart] - 0.0;
376:     VecSetValue(F,indices[i],value,INSERT_VALUES);
377:   }
378:   ISRestoreIndices(userctx->bdis,&indices);
379:   VecRestoreArray(X,&array);
380:   if (vecstate>0) {
381:     VecLockPush(X);
382:   }
383:   VecAssemblyBegin(F);
384:   VecAssemblyEnd(F);
385:   return(0);
386: }

388: PetscErrorCode FormFunctionB(SNES snes,Vec X,Vec F,void *ctx)
389: {
391:   DM             dm;
392:   PetscDS        prob;
393:   PetscInt       nindices,iStart,iEnd,i;
394:   AppCtx         *userctx = (AppCtx *)ctx;
395:   PetscScalar    value;
396:   const PetscInt *indices;

399:   SNESGetDM(snes,&dm);
400:   DMGetDS(dm,&prob);
401:   /* hook functions */
402:   PetscDSSetResidual(prob,0,f0_u,NULL);
403:   FormFunction(snes,X,F,ctx);
404:   /* Boundary condition */
405:   VecGetOwnershipRange(F,&iStart,&iEnd);
406:   ISGetLocalSize(userctx->bdis,&nindices);
407:   ISGetIndices(userctx->bdis,&indices);
408:   for (i=0;i<nindices;i++) {
409:     value = 0.0;
410:     VecSetValue(F,indices[i],value,INSERT_VALUES);
411:   }
412:   ISRestoreIndices(userctx->bdis,&indices);
413:   VecAssemblyBegin(F);
414:   VecAssemblyEnd(F);
415:   return(0);
416: }