Actual source code: ex34.c

slepc-3.12.2 2020-01-13
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 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 FormFunctionAB(SNES,Vec,Vec,Vec,void*);
 37: PetscErrorCode BoundaryGlobalIndex(DM,const char*,IS*);

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

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

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

 66:   DMCreateMatrix(dm,&A);
 67:   MatDuplicate(A,MAT_COPY_VALUES,&B);

 69:   /*
 70:      Compose callback functions and context that will be needed by the solver
 71:   */
 72:   PetscObjectComposeFunction((PetscObject)A,"formFunction",FormFunctionA);
 73:   PetscOptionsGetBool(NULL,NULL,"-form_function_ab",&flg,NULL);
 74:   if (flg) {
 75:     PetscObjectComposeFunction((PetscObject)A,"formFunctionAB",FormFunctionAB);
 76:   }
 77:   PetscObjectComposeFunction((PetscObject)A,"formJacobian",FormJacobianA);
 78:   PetscObjectComposeFunction((PetscObject)B,"formFunction",FormFunctionB);
 79:   PetscContainerCreate(comm,&container);
 80:   PetscContainerSetPointer(container,&user);
 81:   PetscObjectCompose((PetscObject)A,"formFunctionCtx",(PetscObject)container);
 82:   PetscObjectCompose((PetscObject)A,"formJacobianCtx",(PetscObject)container);
 83:   PetscObjectCompose((PetscObject)B,"formFunctionCtx",(PetscObject)container);
 84:   PetscContainerDestroy(&container);

 86:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 87:                 Create the eigensolver and set various options
 88:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 90:   EPSCreate(comm,&eps);
 91:   EPSSetOperators(eps,A,B);
 92:   EPSSetProblemType(eps,EPS_GNHEP);
 93:   /*
 94:      Use nonlinear inverse iteration
 95:   */
 96:   EPSSetType(eps,EPSPOWER);
 97:   EPSPowerSetNonlinear(eps,PETSC_TRUE);
 98:   /*
 99:     Attach DM to SNES
100:   */
101:   EPSPowerGetSNES(eps,&snes);
102:   SNESSetDM(snes,dm);
103:   EPSSetFromOptions(eps);

105:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106:                       Solve the eigensystem
107:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

109:   EPSSolve(eps);

111:   /*
112:      Optional: Get some information from the solver and display it
113:   */
114:   EPSGetType(eps,&type);
115:   EPSGetTolerances(eps,&tol,NULL);
116:   EPSPowerGetNonlinear(eps,&nonlin);
117:   EPSPowerGetUpdate(eps,&update);
118:   PetscPrintf(comm," Solution method: %s%s\n\n",type,nonlin?(update?" (nonlinear with monolithic update)":" (nonlinear)"):"");
119:   EPSGetDimensions(eps,&nev,NULL,NULL);
120:   PetscPrintf(comm," Number of requested eigenvalues: %D\n",nev);

122:   /* print eigenvalue and error */
123:   EPSGetConverged(eps,&nconv);
124:   if (nconv>0) {
125:     PetscScalar   k;
126:     PetscReal     na,nb;
127:     Vec           a,b,eigen;
128:     DMCreateGlobalVector(dm,&a);
129:     VecDuplicate(a,&b);
130:     VecDuplicate(a,&eigen);
131:     EPSGetEigenpair(eps,0,&k,NULL,eigen,NULL);
132:     FormFunctionA(snes,eigen,a,&user);
133:     FormFunctionB(snes,eigen,b,&user);
134:     VecAXPY(a,-k,b);
135:     VecNorm(a,NORM_2,&na);
136:     VecNorm(b,NORM_2,&nb);
137:     relerr = na/(nb*PetscAbsScalar(k));
138:     if (relerr<10*tol) {
139:       PetscPrintf(comm,"k: %g, relative error below tol\n",(double)PetscRealPart(k));
140:     } else {
141:       PetscPrintf(comm,"k: %g, relative error: %g\n",(double)PetscRealPart(k),(double)relerr);
142:     }
143:     VecDestroy(&a);
144:     VecDestroy(&b);
145:     VecDestroy(&eigen);
146:   } else {
147:     PetscPrintf(comm,"Solver did not converge\n");
148:   }

150:   MatDestroy(&A);
151:   MatDestroy(&B);
152:   DMDestroy(&dm);
153:   EPSDestroy(&eps);
154:   ISDestroy(&user.bdis);
155:   SlepcFinalize();
156:   return ierr;
157: }

159: /* <|u|u, v> */
160: static void f0_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
161:                  const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
162:                  const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
163:                  PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
164: {
165:   PetscScalar cof = PetscAbsScalar(u[0]);

167:   f0[0] = cof*u[0];
168: }

170: /* <|\nabla u| \nabla u, \nabla v> */
171: static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
172:                  const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
173:                  const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
174:                  PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
175: {
176:   PetscInt    d;
177:   PetscScalar cof = 0;
178:   for (d = 0; d < dim; ++d)  cof += u_x[d]*u_x[d];

180:   cof = PetscSqrtScalar(cof);

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

185: /* approximate  Jacobian for   <|u|u, v> */
186: static void g0_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
187:                   const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
188:                   const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
189:                   PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
190: {
191:   g0[0] = 1.0*PetscAbsScalar(u[0]);
192: }

194: /* approximate  Jacobian for   <|\nabla u| \nabla u, \nabla v> */
195: static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
196:                   const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
197:                   const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
198:                   PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
199: {
200:   PetscInt d;
201:   for (d = 0; d < dim; ++d) g3[d*dim+d] = 1.0;
202: }

204: PetscErrorCode SetupDiscretization(DM dm)
205: {
206:   PetscFE        fe;
207:   MPI_Comm       comm;

211:   /* Create finite element */
212:   PetscObjectGetComm((PetscObject)dm,&comm);
213:   PetscFECreateDefault(comm,2,1,PETSC_FALSE,NULL,-1,&fe);
214:   PetscObjectSetName((PetscObject)fe,"u");
215:   DMSetField(dm,0,NULL,(PetscObject)fe);
216:   DMCreateDS(dm);
217:   PetscFEDestroy(&fe);
218:   return(0);
219: }

221: PetscErrorCode CreateSquareMesh(MPI_Comm comm,DM *dm)
222: {
223:   PetscInt       cells[] = {8,8};
224:   PetscInt       dim = 2;
225:   DM             pdm;
226:   PetscMPIInt    size;

230:   DMPlexCreateBoxMesh(comm,dim,PETSC_FALSE,cells,NULL,NULL,NULL,PETSC_TRUE,dm);
231:   DMSetFromOptions(*dm);
232:   DMSetUp(*dm);
233:   MPI_Comm_size(comm,&size);
234:   if (size > 1) {
235:     DMPlexDistribute(*dm,0,NULL,&pdm);
236:     DMDestroy(dm);
237:     *dm = pdm;
238:   }
239:   return(0);
240: }

242: PetscErrorCode BoundaryGlobalIndex(DM dm,const char labelname[],IS *bdis)
243: {
244:   IS             bdpoints;
245:   PetscInt       nindices,*indices,numDof,offset,npoints,i,j;
246:   const PetscInt *bdpoints_indices;
247:   DMLabel        bdmarker;
248:   PetscSection   gsection;

252:   DMGetGlobalSection(dm,&gsection);
253:   DMGetLabel(dm,labelname,&bdmarker);
254:   DMLabelGetStratumIS(bdmarker,1,&bdpoints);
255:   ISGetLocalSize(bdpoints,&npoints);
256:   ISGetIndices(bdpoints,&bdpoints_indices);
257:   nindices = 0;
258:   for (i=0;i<npoints;i++) {
259:     PetscSectionGetDof(gsection,bdpoints_indices[i],&numDof);
260:     if (numDof<=0) continue;
261:     nindices += numDof;
262:   }
263:   PetscCalloc1(nindices,&indices);
264:   nindices = 0;
265:   for (i=0;i<npoints;i++) {
266:     PetscSectionGetDof(gsection,bdpoints_indices[i],&numDof);
267:     if (numDof<=0) continue;
268:     PetscSectionGetOffset(gsection,bdpoints_indices[i],&offset);
269:     for (j=0;j<numDof;j++) indices[nindices++] = offset+j;
270:   }
271:   ISRestoreIndices(bdpoints,&bdpoints_indices);
272:   ISDestroy(&bdpoints);
273:   ISCreateGeneral(PetscObjectComm((PetscObject)dm),nindices,indices,PETSC_OWN_POINTER,bdis);
274:   return(0);
275: }

277: static PetscErrorCode FormJacobian(SNES snes,Vec X,Mat A,Mat B,void *ctx)
278: {
279:   DM             dm;
280:   Vec            Xloc;

284:   SNESGetDM(snes,&dm);
285:   DMGetLocalVector(dm,&Xloc);
286:   VecZeroEntries(Xloc);
287:   DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
288:   DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
289:   CHKMEMQ;
290:   DMPlexSNESComputeJacobianFEM(dm,Xloc,A,B,ctx);
291:   CHKMEMQ;
292:   DMRestoreLocalVector(dm,&Xloc);
293:   if (A != B) {
294:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
295:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
296:   }
297:   return(0);
298: }

300: PetscErrorCode FormJacobianA(SNES snes,Vec X,Mat A,Mat B,void *ctx)
301: {
303:   DM             dm;
304:   PetscDS        prob;
305:   AppCtx         *userctx = (AppCtx *)ctx;

308:   MatSetOption(B,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE);
309:   SNESGetDM(snes,&dm);
310:   DMGetDS(dm,&prob);
311:   PetscDSSetJacobian(prob,0,0,NULL,NULL,NULL,g3_uu);
312:   FormJacobian(snes,X,A,B,ctx);
313:   MatZeroRowsIS(B,userctx->bdis,1.0,NULL,NULL);
314:   return(0);
315: }

317: PetscErrorCode FormJacobianB(SNES snes,Vec X,Mat A,Mat B,void *ctx)
318: {
320:   DM             dm;
321:   PetscDS        prob;
322:   AppCtx         *userctx = (AppCtx *)ctx;

325:   MatSetOption(B,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE);
326:   SNESGetDM(snes,&dm);
327:   DMGetDS(dm,&prob);
328:   PetscDSSetJacobian(prob,0,0,g0_uu,NULL,NULL,NULL);
329:   FormJacobian(snes,X,A,B,ctx);
330:   MatZeroRowsIS(B,userctx->bdis,0.0,NULL,NULL);
331:   return(0);
332: }

334: PetscErrorCode FormFunctionAB(SNES snes,Vec x,Vec Ax,Vec Bx,void *ctx)
335: {

339:   /*
340:    * In real applications, users should have a generic formFunctionAB which
341:    * forms Ax and Bx simultaneously for an more efficient calculation.
342:    * In this example, we just call FormFunctionA+FormFunctionB to mimic how
343:    * to use FormFunctionAB
344:    */
345:   FormFunctionA(snes,x,Ax,ctx);
346:   FormFunctionB(snes,x,Bx,ctx);
347:   return(0);
348: }


351: static PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ctx)
352: {
353:   DM             dm;
354:   Vec            Xloc,Floc;

358:   SNESGetDM(snes,&dm);
359:   DMGetLocalVector(dm,&Xloc);
360:   DMGetLocalVector(dm,&Floc);
361:   VecZeroEntries(Xloc);
362:   VecZeroEntries(Floc);
363:   DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
364:   DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
365:   CHKMEMQ;
366:   DMPlexSNESComputeResidualFEM(dm,Xloc,Floc,ctx);
367:   CHKMEMQ;
368:   VecZeroEntries(F);
369:   DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
370:   DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
371:   DMRestoreLocalVector(dm,&Xloc);
372:   DMRestoreLocalVector(dm,&Floc);
373:   return(0);
374: }

376: PetscErrorCode FormFunctionA(SNES snes,Vec X,Vec F,void *ctx)
377: {
379:   DM             dm;
380:   PetscDS        prob;
381:   PetscInt       nindices,iStart,iEnd,i;
382:   AppCtx         *userctx = (AppCtx *)ctx;
383:   PetscScalar    *array,value;
384:   const PetscInt *indices;
385:   PetscInt       vecstate;

388:   SNESGetDM(snes,&dm);
389:   DMGetDS(dm,&prob);
390:   /* hook functions */
391:   PetscDSSetResidual(prob,0,NULL,f1_u);
392:   FormFunction(snes,X,F,ctx);
393:   /* Boundary condition */
394:   VecLockGet(X,&vecstate);
395:   if (vecstate>0) {
396:     VecLockReadPop(X);
397:   }
398:   VecGetOwnershipRange(X,&iStart,&iEnd);
399:   VecGetArray(X,&array);
400:   ISGetLocalSize(userctx->bdis,&nindices);
401:   ISGetIndices(userctx->bdis,&indices);
402:   for (i=0;i<nindices;i++) {
403:     value = array[indices[i]-iStart] - 0.0;
404:     VecSetValue(F,indices[i],value,INSERT_VALUES);
405:   }
406:   ISRestoreIndices(userctx->bdis,&indices);
407:   VecRestoreArray(X,&array);
408:   if (vecstate>0) {
409:     VecLockReadPush(X);
410:   }
411:   VecAssemblyBegin(F);
412:   VecAssemblyEnd(F);
413:   return(0);
414: }

416: PetscErrorCode FormFunctionB(SNES snes,Vec X,Vec F,void *ctx)
417: {
419:   DM             dm;
420:   PetscDS        prob;
421:   PetscInt       nindices,iStart,iEnd,i;
422:   AppCtx         *userctx = (AppCtx *)ctx;
423:   PetscScalar    value;
424:   const PetscInt *indices;

427:   SNESGetDM(snes,&dm);
428:   DMGetDS(dm,&prob);
429:   /* hook functions */
430:   PetscDSSetResidual(prob,0,f0_u,NULL);
431:   FormFunction(snes,X,F,ctx);
432:   /* Boundary condition */
433:   VecGetOwnershipRange(F,&iStart,&iEnd);
434:   ISGetLocalSize(userctx->bdis,&nindices);
435:   ISGetIndices(userctx->bdis,&indices);
436:   for (i=0;i<nindices;i++) {
437:     value = 0.0;
438:     VecSetValue(F,indices[i],value,INSERT_VALUES);
439:   }
440:   ISRestoreIndices(userctx->bdis,&indices);
441:   VecAssemblyBegin(F);
442:   VecAssemblyEnd(F);
443:   return(0);
444: }

446: /*TEST

448:    testset:
449:       requires: double
450:       args: -petscspace_degree 1 -petscspace_poly_tensor
451:       output_file: output/ex34_1.out
452:       test:
453:          suffix: 1
454:       test:
455:          suffix: 2
456:          args: -eps_power_update -form_function_ab {{0 1}}
457:          filter: sed -e "s/ with monolithic update//"

459: TEST*/