Actual source code: ex34.c
slepc-3.11.2 2019-07-30
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;
54: SNES snes;
57: SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
58: comm = PETSC_COMM_WORLD;
59: /* Create a quadrilateral mesh on domain (0,1)x(0,1) */
60: CreateSquareMesh(comm,&dm);
61: /* Setup basis function */
62: SetupDiscretization(dm);
63: BoundaryGlobalIndex(dm,"marker",&user.bdis);
65: DMCreateMatrix(dm,&A);
66: MatDuplicate(A,MAT_COPY_VALUES,&B);
68: /*
69: Compose callback functions and context that will be needed by the solver
70: */
71: PetscObjectComposeFunction((PetscObject)A,"formFunction",FormFunctionA);
72: PetscObjectComposeFunction((PetscObject)A,"formFunctionAB",FormFunctionAB);
73: PetscObjectComposeFunction((PetscObject)A,"formJacobian",FormJacobianA);
74: PetscObjectComposeFunction((PetscObject)B,"formFunction",FormFunctionB);
75: PetscContainerCreate(comm,&container);
76: PetscContainerSetPointer(container,&user);
77: PetscObjectCompose((PetscObject)A,"formFunctionCtx",(PetscObject)container);
78: PetscObjectCompose((PetscObject)A,"formJacobianCtx",(PetscObject)container);
79: PetscObjectCompose((PetscObject)B,"formFunctionCtx",(PetscObject)container);
80: PetscContainerDestroy(&container);
82: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
83: Create the eigensolver and set various options
84: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
86: EPSCreate(comm,&eps);
87: EPSSetOperators(eps,A,B);
88: EPSSetProblemType(eps,EPS_GNHEP);
89: /*
90: Use nonlinear inverse iteration
91: */
92: EPSSetType(eps,EPSPOWER);
93: EPSPowerSetNonlinear(eps,PETSC_TRUE);
94: /*
95: Attach DM to SNES
96: */
97: EPSPowerGetSNES(eps,&snes);
98: SNESSetDM(snes,dm);
99: EPSSetFromOptions(eps);
101: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102: Solve the eigensystem
103: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
105: EPSSolve(eps);
107: /*
108: Optional: Get some information from the solver and display it
109: */
110: EPSGetType(eps,&type);
111: EPSPowerGetNonlinear(eps,&nonlin);
112: PetscPrintf(comm," Solution method: %s%s\n\n",type,nonlin?" (nonlinear)":"");
113: EPSGetDimensions(eps,&nev,NULL,NULL);
114: PetscPrintf(comm," Number of requested eigenvalues: %D\n",nev);
116: /* print eigenvalue and error */
117: EPSGetConverged(eps,&nconv);
118: if (nconv>0) {
119: PetscScalar k;
120: PetscReal na,nb;
121: Vec a,b,eigen;
122: DMCreateGlobalVector(dm,&a);
123: VecDuplicate(a,&b);
124: VecDuplicate(a,&eigen);
125: EPSGetEigenpair(eps,0,&k,NULL,eigen,NULL);
126: FormFunctionA(snes,eigen,a,&user);
127: FormFunctionB(snes,eigen,b,&user);
128: VecAXPY(a,-k,b);
129: VecNorm(a,NORM_2,&na);
130: VecNorm(b,NORM_2,&nb);
131: PetscPrintf(comm,"k: %g error: %g\n",(double)PetscRealPart(k),(double)na/nb);
132: VecDestroy(&a);
133: VecDestroy(&b);
134: VecDestroy(&eigen);
135: } else {
136: PetscPrintf(comm,"Solver did not converge\n");
137: }
139: MatDestroy(&A);
140: MatDestroy(&B);
141: DMDestroy(&dm);
142: EPSDestroy(&eps);
143: ISDestroy(&user.bdis);
144: SlepcFinalize();
145: return ierr;
146: }
148: /* <|u|u, v> */
149: static void f0_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
150: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
151: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
152: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
153: {
154: PetscScalar cof = PetscAbsScalar(u[0]);
156: f0[0] = cof*u[0];
157: }
159: /* <|\nabla u| \nabla u, \nabla v> */
160: static void f1_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 f1[])
164: {
165: PetscInt d;
166: PetscScalar cof = 0;
167: for (d = 0; d < dim; ++d) cof += u_x[d]*u_x[d];
169: cof = PetscSqrtScalar(cof);
171: for (d = 0; d < dim; ++d) f1[d] = u_x[d]*cof;
172: }
174: /* approximate Jacobian for <|u|u, v> */
175: static void g0_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
176: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
177: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
178: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
179: {
180: g0[0] = 1.0*PetscAbsScalar(u[0]);
181: }
183: /* approximate Jacobian for <|\nabla u| \nabla u, \nabla v> */
184: static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
185: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
186: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
187: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
188: {
189: PetscInt d;
190: for (d = 0; d < dim; ++d) g3[d*dim+d] = 1.0;
191: }
193: PetscErrorCode SetupDiscretization(DM dm)
194: {
195: PetscFE fe;
196: MPI_Comm comm;
200: /* Create finite element */
201: PetscObjectGetComm((PetscObject)dm,&comm);
202: PetscFECreateDefault(comm,2,1,PETSC_FALSE,NULL,-1,&fe);
203: PetscObjectSetName((PetscObject)fe,"u");
204: DMSetField(dm,0,NULL,(PetscObject)fe);
205: DMCreateDS(dm);
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: DMPlexCreateBoxMesh(comm,dim,PETSC_FALSE,cells,NULL,NULL,NULL,PETSC_TRUE,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: PetscErrorCode FormFunctionAB(SNES snes,Vec x,Vec Ax,Vec Bx,void *ctx)
324: {
328: /*
329: * In real applications, users should have a generic formFunctionAB which
330: * forms Ax and Bx simultaneously for an more efficient calculation.
331: * In this example, we just call FormFunctionA+FormFunctionB to mimic how
332: * to use FormFunctionAB
333: */
334: FormFunctionA(snes,x,Ax,ctx);
335: FormFunctionB(snes,x,Bx,ctx);
336: return(0);
337: }
340: static PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ctx)
341: {
342: DM dm;
343: Vec Xloc,Floc;
347: SNESGetDM(snes,&dm);
348: DMGetLocalVector(dm,&Xloc);
349: DMGetLocalVector(dm,&Floc);
350: VecZeroEntries(Xloc);
351: VecZeroEntries(Floc);
352: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
353: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
354: CHKMEMQ;
355: DMPlexSNESComputeResidualFEM(dm,Xloc,Floc,ctx);
356: CHKMEMQ;
357: VecZeroEntries(F);
358: DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
359: DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
360: DMRestoreLocalVector(dm,&Xloc);
361: DMRestoreLocalVector(dm,&Floc);
362: return(0);
363: }
365: PetscErrorCode FormFunctionA(SNES snes,Vec X,Vec F,void *ctx)
366: {
368: DM dm;
369: PetscDS prob;
370: PetscInt nindices,iStart,iEnd,i;
371: AppCtx *userctx = (AppCtx *)ctx;
372: PetscScalar *array,value;
373: const PetscInt *indices;
374: PetscInt vecstate;
377: SNESGetDM(snes,&dm);
378: DMGetDS(dm,&prob);
379: /* hook functions */
380: PetscDSSetResidual(prob,0,NULL,f1_u);
381: FormFunction(snes,X,F,ctx);
382: /* Boundary condition */
383: VecLockGet(X,&vecstate);
384: if (vecstate>0) {
385: VecLockReadPop(X);
386: }
387: VecGetOwnershipRange(X,&iStart,&iEnd);
388: VecGetArray(X,&array);
389: ISGetLocalSize(userctx->bdis,&nindices);
390: ISGetIndices(userctx->bdis,&indices);
391: for (i=0;i<nindices;i++) {
392: value = array[indices[i]-iStart] - 0.0;
393: VecSetValue(F,indices[i],value,INSERT_VALUES);
394: }
395: ISRestoreIndices(userctx->bdis,&indices);
396: VecRestoreArray(X,&array);
397: if (vecstate>0) {
398: VecLockReadPush(X);
399: }
400: VecAssemblyBegin(F);
401: VecAssemblyEnd(F);
402: return(0);
403: }
405: PetscErrorCode FormFunctionB(SNES snes,Vec X,Vec F,void *ctx)
406: {
408: DM dm;
409: PetscDS prob;
410: PetscInt nindices,iStart,iEnd,i;
411: AppCtx *userctx = (AppCtx *)ctx;
412: PetscScalar value;
413: const PetscInt *indices;
416: SNESGetDM(snes,&dm);
417: DMGetDS(dm,&prob);
418: /* hook functions */
419: PetscDSSetResidual(prob,0,f0_u,NULL);
420: FormFunction(snes,X,F,ctx);
421: /* Boundary condition */
422: VecGetOwnershipRange(F,&iStart,&iEnd);
423: ISGetLocalSize(userctx->bdis,&nindices);
424: ISGetIndices(userctx->bdis,&indices);
425: for (i=0;i<nindices;i++) {
426: value = 0.0;
427: VecSetValue(F,indices[i],value,INSERT_VALUES);
428: }
429: ISRestoreIndices(userctx->bdis,&indices);
430: VecAssemblyBegin(F);
431: VecAssemblyEnd(F);
432: return(0);
433: }
435: /*TEST
437: test:
438: suffix: 1
439: args: -petscspace_degree 1 -petscspace_poly_tensor
440: requires: double !complex
442: test:
443: suffix: 2
444: args: -petscspace_degree 1 -petscspace_poly_tensor -eps_power_update
445: requires: double !complex
447: TEST*/