Actual source code: ex34.c
slepc-3.8.3 2018-04-03
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: }