Actual source code: ex34.c
slepc-3.10.1 2018-10-23
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2018, 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: PetscDS prob;
197: PetscInt totDim;
198: MPI_Comm comm;
202: /* Create finite element */
203: PetscObjectGetComm((PetscObject)dm,&comm);
204: PetscFECreateDefault(comm,2,1,PETSC_FALSE,NULL,-1,&fe);
205: PetscObjectSetName((PetscObject)fe,"u");
206: DMGetDS(dm,&prob);
207: PetscDSSetDiscretization(prob,0,(PetscObject)fe);
208: DMSetDS(dm,prob);
209: PetscDSGetTotalDimension(prob,&totDim);
210: PetscFEDestroy(&fe);
211: return(0);
212: }
214: PetscErrorCode CreateSquareMesh(MPI_Comm comm,DM *dm)
215: {
216: PetscInt cells[] = {8,8};
217: PetscInt dim = 2;
218: DM pdm;
219: PetscMPIInt size;
223: DMPlexCreateBoxMesh(comm,dim,PETSC_FALSE,cells,NULL,NULL,NULL,PETSC_TRUE,dm);
224: DMSetFromOptions(*dm);
225: DMSetUp(*dm);
226: MPI_Comm_size(comm,&size);
227: if (size > 1) {
228: DMPlexDistribute(*dm,0,NULL,&pdm);
229: DMDestroy(dm);
230: *dm = pdm;
231: }
232: return(0);
233: }
235: PetscErrorCode BoundaryGlobalIndex(DM dm,const char labelname[],IS *bdis)
236: {
237: IS bdpoints;
238: PetscInt nindices,*indices,numDof,offset,npoints,i,j;
239: const PetscInt *bdpoints_indices;
240: DMLabel bdmarker;
241: PetscSection gsection;
245: DMGetDefaultGlobalSection(dm,&gsection);
246: DMGetLabel(dm,labelname,&bdmarker);
247: DMLabelGetStratumIS(bdmarker,1,&bdpoints);
248: ISGetLocalSize(bdpoints,&npoints);
249: ISGetIndices(bdpoints,&bdpoints_indices);
250: nindices = 0;
251: for (i=0;i<npoints;i++) {
252: PetscSectionGetDof(gsection,bdpoints_indices[i],&numDof);
253: if (numDof<=0) continue;
254: nindices += numDof;
255: }
256: PetscCalloc1(nindices,&indices);
257: nindices = 0;
258: for (i=0;i<npoints;i++) {
259: PetscSectionGetDof(gsection,bdpoints_indices[i],&numDof);
260: if (numDof<=0) continue;
261: PetscSectionGetOffset(gsection,bdpoints_indices[i],&offset);
262: for (j=0;j<numDof;j++) indices[nindices++] = offset+j;
263: }
264: ISRestoreIndices(bdpoints,&bdpoints_indices);
265: ISDestroy(&bdpoints);
266: ISCreateGeneral(PetscObjectComm((PetscObject)dm),nindices,indices,PETSC_OWN_POINTER,bdis);
267: return(0);
268: }
270: static PetscErrorCode FormJacobian(SNES snes,Vec X,Mat A,Mat B,void *ctx)
271: {
272: DM dm;
273: Vec Xloc;
277: SNESGetDM(snes,&dm);
278: DMGetLocalVector(dm,&Xloc);
279: VecZeroEntries(Xloc);
280: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
281: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
282: CHKMEMQ;
283: DMPlexSNESComputeJacobianFEM(dm,Xloc,A,B,ctx);
284: CHKMEMQ;
285: DMRestoreLocalVector(dm,&Xloc);
286: if (A != B) {
287: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
288: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
289: }
290: return(0);
291: }
293: PetscErrorCode FormJacobianA(SNES snes,Vec X,Mat A,Mat B,void *ctx)
294: {
296: DM dm;
297: PetscDS prob;
298: AppCtx *userctx = (AppCtx *)ctx;
301: MatSetOption(B,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE);
302: SNESGetDM(snes,&dm);
303: DMGetDS(dm,&prob);
304: PetscDSSetJacobian(prob,0,0,NULL,NULL,NULL,g3_uu);
305: FormJacobian(snes,X,A,B,ctx);
306: MatZeroRowsIS(B,userctx->bdis,1.0,NULL,NULL);
307: return(0);
308: }
310: PetscErrorCode FormJacobianB(SNES snes,Vec X,Mat A,Mat B,void *ctx)
311: {
313: DM dm;
314: PetscDS prob;
315: AppCtx *userctx = (AppCtx *)ctx;
318: MatSetOption(B,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE);
319: SNESGetDM(snes,&dm);
320: DMGetDS(dm,&prob);
321: PetscDSSetJacobian(prob,0,0,g0_uu,NULL,NULL,NULL);
322: FormJacobian(snes,X,A,B,ctx);
323: MatZeroRowsIS(B,userctx->bdis,0.0,NULL,NULL);
324: return(0);
325: }
327: PetscErrorCode FormFunctionAB(SNES snes,Vec x,Vec Ax,Vec Bx,void *ctx)
328: {
332: /*
333: * In real applications, users should have a generic formFunctionAB which
334: * forms Ax and Bx simultaneously for an more efficient calculation.
335: * In this example, we just call FormFunctionA+FormFunctionB to mimic how
336: * to use FormFunctionAB
337: */
338: FormFunctionA(snes,x,Ax,ctx);
339: FormFunctionB(snes,x,Bx,ctx);
340: return(0);
341: }
344: static PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ctx)
345: {
346: DM dm;
347: Vec Xloc,Floc;
351: SNESGetDM(snes,&dm);
352: DMGetLocalVector(dm,&Xloc);
353: DMGetLocalVector(dm,&Floc);
354: VecZeroEntries(Xloc);
355: VecZeroEntries(Floc);
356: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
357: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
358: CHKMEMQ;
359: DMPlexSNESComputeResidualFEM(dm,Xloc,Floc,ctx);
360: CHKMEMQ;
361: VecZeroEntries(F);
362: DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
363: DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
364: DMRestoreLocalVector(dm,&Xloc);
365: DMRestoreLocalVector(dm,&Floc);
366: return(0);
367: }
369: PetscErrorCode FormFunctionA(SNES snes,Vec X,Vec F,void *ctx)
370: {
372: DM dm;
373: PetscDS prob;
374: PetscInt nindices,iStart,iEnd,i;
375: AppCtx *userctx = (AppCtx *)ctx;
376: PetscScalar *array,value;
377: const PetscInt *indices;
378: PetscInt vecstate;
381: SNESGetDM(snes,&dm);
382: DMGetDS(dm,&prob);
383: /* hook functions */
384: PetscDSSetResidual(prob,0,NULL,f1_u);
385: FormFunction(snes,X,F,ctx);
386: /* Boundary condition */
387: VecLockGet(X,&vecstate);
388: if (vecstate>0) {
389: VecLockPop(X);
390: }
391: VecGetOwnershipRange(X,&iStart,&iEnd);
392: VecGetArray(X,&array);
393: ISGetLocalSize(userctx->bdis,&nindices);
394: ISGetIndices(userctx->bdis,&indices);
395: for (i=0;i<nindices;i++) {
396: value = array[indices[i]-iStart] - 0.0;
397: VecSetValue(F,indices[i],value,INSERT_VALUES);
398: }
399: ISRestoreIndices(userctx->bdis,&indices);
400: VecRestoreArray(X,&array);
401: if (vecstate>0) {
402: VecLockPush(X);
403: }
404: VecAssemblyBegin(F);
405: VecAssemblyEnd(F);
406: return(0);
407: }
409: PetscErrorCode FormFunctionB(SNES snes,Vec X,Vec F,void *ctx)
410: {
412: DM dm;
413: PetscDS prob;
414: PetscInt nindices,iStart,iEnd,i;
415: AppCtx *userctx = (AppCtx *)ctx;
416: PetscScalar value;
417: const PetscInt *indices;
420: SNESGetDM(snes,&dm);
421: DMGetDS(dm,&prob);
422: /* hook functions */
423: PetscDSSetResidual(prob,0,f0_u,NULL);
424: FormFunction(snes,X,F,ctx);
425: /* Boundary condition */
426: VecGetOwnershipRange(F,&iStart,&iEnd);
427: ISGetLocalSize(userctx->bdis,&nindices);
428: ISGetIndices(userctx->bdis,&indices);
429: for (i=0;i<nindices;i++) {
430: value = 0.0;
431: VecSetValue(F,indices[i],value,INSERT_VALUES);
432: }
433: ISRestoreIndices(userctx->bdis,&indices);
434: VecAssemblyBegin(F);
435: VecAssemblyEnd(F);
436: return(0);
437: }
439: /*TEST
441: test:
442: suffix: 1
443: args: -petscspace_degree 1 -petscspace_poly_tensor
444: requires: double !complex
446: test:
447: suffix: 2
448: args: -petscspace_degree 1 -petscspace_poly_tensor -eps_power_update
449: requires: double !complex
451: TEST*/