Actual source code: ex34.c
slepc-3.14.1 2020-12-08
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2020, 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 MatMult_A(Mat A,Vec x,Vec y);
35: PetscErrorCode FormJacobianB(SNES,Vec,Mat,Mat,void*);
36: PetscErrorCode FormFunctionB(SNES,Vec,Vec,void*);
37: PetscErrorCode MatMult_B(Mat A,Vec x,Vec y);
38: PetscErrorCode FormFunctionAB(SNES,Vec,Vec,Vec,void*);
39: PetscErrorCode BoundaryGlobalIndex(DM,const char*,IS*);
41: typedef struct {
42: IS bdis; /* global indices for boundary DoFs */
43: SNES snes;
44: } AppCtx;
46: int main(int argc,char **argv)
47: {
48: DM dm;
49: MPI_Comm comm;
50: AppCtx user;
51: EPS eps; /* eigenproblem solver context */
52: ST st;
53: EPSType type;
54: Mat A,B,P;
55: Vec v0;
56: PetscContainer container;
57: PetscInt nev,nconv,m,n,M,N;
58: PetscBool nonlin,flg=PETSC_FALSE,update;
59: SNES snes;
60: PetscReal tol,relerr;
61: PetscBool use_shell_matrix=PETSC_FALSE,test_init_sol=PETSC_FALSE;
64: SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
65: comm = PETSC_COMM_WORLD;
66: /* Create a quadrilateral mesh on domain (0,1)x(0,1) */
67: CreateSquareMesh(comm,&dm);
68: /* Setup basis function */
69: SetupDiscretization(dm);
70: BoundaryGlobalIndex(dm,"marker",&user.bdis);
71: /* Check if we are going to use shell matrices */
72: PetscOptionsGetBool(NULL,NULL,"-use_shell_matrix",&use_shell_matrix,NULL);
73: if (use_shell_matrix) {
74: DMCreateMatrix(dm,&P);
75: MatGetLocalSize(P,&m,&n);
76: MatGetSize(P,&M,&N);
77: MatCreateShell(comm,m,n,M,N,&user,&A);
78: MatShellSetOperation(A,MATOP_MULT,(void(*)(void))MatMult_A);
79: MatCreateShell(comm,m,n,M,N,&user,&B);
80: MatShellSetOperation(B,MATOP_MULT,(void(*)(void))MatMult_B);
81: } else {
82: DMCreateMatrix(dm,&A);
83: MatDuplicate(A,MAT_COPY_VALUES,&B);
84: }
86: /*
87: Compose callback functions and context that will be needed by the solver
88: */
89: PetscObjectComposeFunction((PetscObject)A,"formFunction",FormFunctionA);
90: PetscOptionsGetBool(NULL,NULL,"-form_function_ab",&flg,NULL);
91: if (flg) {
92: PetscObjectComposeFunction((PetscObject)A,"formFunctionAB",FormFunctionAB);
93: }
94: PetscObjectComposeFunction((PetscObject)A,"formJacobian",FormJacobianA);
95: PetscObjectComposeFunction((PetscObject)B,"formFunction",FormFunctionB);
96: PetscContainerCreate(comm,&container);
97: PetscContainerSetPointer(container,&user);
98: PetscObjectCompose((PetscObject)A,"formFunctionCtx",(PetscObject)container);
99: PetscObjectCompose((PetscObject)A,"formJacobianCtx",(PetscObject)container);
100: PetscObjectCompose((PetscObject)B,"formFunctionCtx",(PetscObject)container);
101: PetscContainerDestroy(&container);
103: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104: Create the eigensolver and set various options
105: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
107: EPSCreate(comm,&eps);
108: EPSSetOperators(eps,A,B);
109: EPSSetProblemType(eps,EPS_GNHEP);
110: /*
111: Use nonlinear inverse iteration
112: */
113: EPSSetType(eps,EPSPOWER);
114: EPSPowerSetNonlinear(eps,PETSC_TRUE);
115: /*
116: Attach DM to SNES
117: */
118: EPSPowerGetSNES(eps,&snes);
119: user.snes = snes;
120: SNESSetDM(snes,dm);
121: EPSSetFromOptions(eps);
123: /* Set a preconditioning matrix to ST */
124: if (use_shell_matrix) {
125: EPSGetST(eps,&st);
126: STPrecondSetMatForPC(st,P);
127: }
129: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130: Solve the eigensystem
131: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
133: EPSSolve(eps);
135: EPSGetConverged(eps,&nconv);
136: PetscOptionsGetBool(NULL,NULL,"-test_init_sol",&test_init_sol,NULL);
137: if (nconv && test_init_sol) {
138: PetscScalar k;
139: PetscReal norm0;
140: PetscInt nits;
142: MatCreateVecs(A,&v0,NULL);
143: EPSGetEigenpair(eps,0,&k,NULL,v0,NULL);
144: EPSSetInitialSpace(eps,1,&v0);
145: VecDestroy(&v0);
146: /* Norm of the previous residual */
147: SNESGetFunctionNorm(snes,&norm0);
148: /* Make the tolerance smaller than the last residual
149: SNES will converge right away if the initial is setup correctly */
150: SNESSetTolerances(snes,norm0*1.2,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
151: EPSSolve(eps);
152: /* Number of Newton iterations supposes to be zero */
153: SNESGetIterationNumber(snes,&nits);
154: if (nits) {
155: PetscPrintf(comm," Number of Newton iterations %D should be zero \n",nits);
156: }
157: }
159: /*
160: Optional: Get some information from the solver and display it
161: */
162: EPSGetType(eps,&type);
163: EPSGetTolerances(eps,&tol,NULL);
164: EPSPowerGetNonlinear(eps,&nonlin);
165: EPSPowerGetUpdate(eps,&update);
166: PetscPrintf(comm," Solution method: %s%s\n\n",type,nonlin?(update?" (nonlinear with monolithic update)":" (nonlinear)"):"");
167: EPSGetDimensions(eps,&nev,NULL,NULL);
168: PetscPrintf(comm," Number of requested eigenvalues: %D\n",nev);
170: /* print eigenvalue and error */
171: EPSGetConverged(eps,&nconv);
172: if (nconv>0) {
173: PetscScalar k;
174: PetscReal na,nb;
175: Vec a,b,eigen;
176: DMCreateGlobalVector(dm,&a);
177: VecDuplicate(a,&b);
178: VecDuplicate(a,&eigen);
179: EPSGetEigenpair(eps,0,&k,NULL,eigen,NULL);
180: FormFunctionA(snes,eigen,a,&user);
181: FormFunctionB(snes,eigen,b,&user);
182: VecAXPY(a,-k,b);
183: VecNorm(a,NORM_2,&na);
184: VecNorm(b,NORM_2,&nb);
185: relerr = na/(nb*PetscAbsScalar(k));
186: if (relerr<10*tol) {
187: PetscPrintf(comm,"k: %g, relative error below tol\n",(double)PetscRealPart(k));
188: } else {
189: PetscPrintf(comm,"k: %g, relative error: %g\n",(double)PetscRealPart(k),(double)relerr);
190: }
191: VecDestroy(&a);
192: VecDestroy(&b);
193: VecDestroy(&eigen);
194: } else {
195: PetscPrintf(comm,"Solver did not converge\n");
196: }
198: MatDestroy(&A);
199: MatDestroy(&B);
200: if (use_shell_matrix) {
201: MatDestroy(&P);
202: }
203: DMDestroy(&dm);
204: EPSDestroy(&eps);
205: ISDestroy(&user.bdis);
206: SlepcFinalize();
207: return ierr;
208: }
210: /* <|u|u, v> */
211: static void f0_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
212: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
213: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
214: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
215: {
216: PetscScalar cof = PetscAbsScalar(u[0]);
218: f0[0] = cof*u[0];
219: }
221: /* <|\nabla u| \nabla u, \nabla v> */
222: static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
223: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
224: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
225: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
226: {
227: PetscInt d;
228: PetscScalar cof = 0;
229: for (d = 0; d < dim; ++d) cof += u_x[d]*u_x[d];
231: cof = PetscSqrtScalar(cof);
233: for (d = 0; d < dim; ++d) f1[d] = u_x[d]*cof;
234: }
236: /* approximate Jacobian for <|u|u, v> */
237: static void g0_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
238: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
239: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
240: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
241: {
242: g0[0] = 1.0*PetscAbsScalar(u[0]);
243: }
245: /* approximate Jacobian for <|\nabla u| \nabla u, \nabla v> */
246: static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
247: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
248: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
249: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
250: {
251: PetscInt d;
253: for (d = 0; d < dim; ++d) g3[d*dim+d] = 1.0;
254: }
256: PetscErrorCode SetupDiscretization(DM dm)
257: {
258: PetscFE fe;
259: MPI_Comm comm;
263: /* Create finite element */
264: PetscObjectGetComm((PetscObject)dm,&comm);
265: PetscFECreateDefault(comm,2,1,PETSC_FALSE,NULL,-1,&fe);
266: PetscObjectSetName((PetscObject)fe,"u");
267: DMSetField(dm,0,NULL,(PetscObject)fe);
268: DMCreateDS(dm);
269: PetscFEDestroy(&fe);
270: return(0);
271: }
273: PetscErrorCode CreateSquareMesh(MPI_Comm comm,DM *dm)
274: {
275: PetscInt cells[] = {8,8};
276: PetscInt dim = 2;
277: DM pdm;
278: PetscMPIInt size;
282: DMPlexCreateBoxMesh(comm,dim,PETSC_FALSE,cells,NULL,NULL,NULL,PETSC_TRUE,dm);
283: DMSetFromOptions(*dm);
284: DMSetUp(*dm);
285: MPI_Comm_size(comm,&size);
286: if (size > 1) {
287: DMPlexDistribute(*dm,0,NULL,&pdm);
288: DMDestroy(dm);
289: *dm = pdm;
290: }
291: return(0);
292: }
294: PetscErrorCode BoundaryGlobalIndex(DM dm,const char labelname[],IS *bdis)
295: {
296: IS bdpoints;
297: PetscInt nindices,*indices,numDof,offset,npoints,i,j;
298: const PetscInt *bdpoints_indices;
299: DMLabel bdmarker;
300: PetscSection gsection;
304: DMGetGlobalSection(dm,&gsection);
305: DMGetLabel(dm,labelname,&bdmarker);
306: DMLabelGetStratumIS(bdmarker,1,&bdpoints);
307: ISGetLocalSize(bdpoints,&npoints);
308: ISGetIndices(bdpoints,&bdpoints_indices);
309: nindices = 0;
310: for (i=0;i<npoints;i++) {
311: PetscSectionGetDof(gsection,bdpoints_indices[i],&numDof);
312: if (numDof<=0) continue;
313: nindices += numDof;
314: }
315: PetscCalloc1(nindices,&indices);
316: nindices = 0;
317: for (i=0;i<npoints;i++) {
318: PetscSectionGetDof(gsection,bdpoints_indices[i],&numDof);
319: if (numDof<=0) continue;
320: PetscSectionGetOffset(gsection,bdpoints_indices[i],&offset);
321: for (j=0;j<numDof;j++) indices[nindices++] = offset+j;
322: }
323: ISRestoreIndices(bdpoints,&bdpoints_indices);
324: ISDestroy(&bdpoints);
325: ISCreateGeneral(PetscObjectComm((PetscObject)dm),nindices,indices,PETSC_OWN_POINTER,bdis);
326: return(0);
327: }
329: static PetscErrorCode FormJacobian(SNES snes,Vec X,Mat A,Mat B,void *ctx)
330: {
331: DM dm;
332: Vec Xloc;
336: SNESGetDM(snes,&dm);
337: DMGetLocalVector(dm,&Xloc);
338: VecZeroEntries(Xloc);
339: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
340: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
341: CHKMEMQ;
342: DMPlexSNESComputeJacobianFEM(dm,Xloc,A,B,ctx);
343: if (A!=B) {
344: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
345: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
346: }
347: CHKMEMQ;
348: DMRestoreLocalVector(dm,&Xloc);
349: return(0);
350: }
352: PetscErrorCode FormJacobianA(SNES snes,Vec X,Mat A,Mat B,void *ctx)
353: {
355: DM dm;
356: PetscDS prob;
357: AppCtx *userctx = (AppCtx *)ctx;
360: MatSetOption(B,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE);
361: SNESGetDM(snes,&dm);
362: DMGetDS(dm,&prob);
363: PetscDSSetJacobian(prob,0,0,NULL,NULL,NULL,g3_uu);
364: FormJacobian(snes,X,A,B,ctx);
365: MatZeroRowsIS(B,userctx->bdis,1.0,NULL,NULL);
366: return(0);
367: }
369: PetscErrorCode FormJacobianB(SNES snes,Vec X,Mat A,Mat B,void *ctx)
370: {
372: DM dm;
373: PetscDS prob;
374: AppCtx *userctx = (AppCtx *)ctx;
377: MatSetOption(B,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE);
378: SNESGetDM(snes,&dm);
379: DMGetDS(dm,&prob);
380: PetscDSSetJacobian(prob,0,0,g0_uu,NULL,NULL,NULL);
381: FormJacobian(snes,X,A,B,ctx);
382: MatZeroRowsIS(B,userctx->bdis,0.0,NULL,NULL);
383: return(0);
384: }
386: PetscErrorCode FormFunctionAB(SNES snes,Vec x,Vec Ax,Vec Bx,void *ctx)
387: {
391: /*
392: * In real applications, users should have a generic formFunctionAB which
393: * forms Ax and Bx simultaneously for an more efficient calculation.
394: * In this example, we just call FormFunctionA+FormFunctionB to mimic how
395: * to use FormFunctionAB
396: */
397: FormFunctionA(snes,x,Ax,ctx);
398: FormFunctionB(snes,x,Bx,ctx);
399: return(0);
400: }
403: static PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ctx)
404: {
405: DM dm;
406: Vec Xloc,Floc;
410: SNESGetDM(snes,&dm);
411: DMGetLocalVector(dm,&Xloc);
412: DMGetLocalVector(dm,&Floc);
413: VecZeroEntries(Xloc);
414: VecZeroEntries(Floc);
415: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
416: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
417: CHKMEMQ;
418: DMPlexSNESComputeResidualFEM(dm,Xloc,Floc,ctx);
419: CHKMEMQ;
420: VecZeroEntries(F);
421: DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
422: DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
423: DMRestoreLocalVector(dm,&Xloc);
424: DMRestoreLocalVector(dm,&Floc);
425: return(0);
426: }
428: PetscErrorCode FormFunctionA(SNES snes,Vec X,Vec F,void *ctx)
429: {
431: DM dm;
432: PetscDS prob;
433: PetscInt nindices,iStart,iEnd,i;
434: AppCtx *userctx = (AppCtx *)ctx;
435: PetscScalar *array,value;
436: const PetscInt *indices;
437: PetscInt vecstate;
440: SNESGetDM(snes,&dm);
441: DMGetDS(dm,&prob);
442: /* hook functions */
443: PetscDSSetResidual(prob,0,NULL,f1_u);
444: FormFunction(snes,X,F,ctx);
445: /* Boundary condition */
446: VecLockGet(X,&vecstate);
447: if (vecstate>0) {
448: VecLockReadPop(X);
449: }
450: VecGetOwnershipRange(X,&iStart,&iEnd);
451: VecGetArray(X,&array);
452: ISGetLocalSize(userctx->bdis,&nindices);
453: ISGetIndices(userctx->bdis,&indices);
454: for (i=0;i<nindices;i++) {
455: value = array[indices[i]-iStart] - 0.0;
456: VecSetValue(F,indices[i],value,INSERT_VALUES);
457: }
458: ISRestoreIndices(userctx->bdis,&indices);
459: VecRestoreArray(X,&array);
460: if (vecstate>0) {
461: VecLockReadPush(X);
462: }
463: VecAssemblyBegin(F);
464: VecAssemblyEnd(F);
465: return(0);
466: }
468: PetscErrorCode MatMult_A(Mat A,Vec x,Vec y)
469: {
471: AppCtx *userctx;
474: MatShellGetContext(A,&userctx);
475: FormFunctionA(userctx->snes,x,y,userctx);
476: return(0);
477: }
479: PetscErrorCode FormFunctionB(SNES snes,Vec X,Vec F,void *ctx)
480: {
482: DM dm;
483: PetscDS prob;
484: PetscInt nindices,iStart,iEnd,i;
485: AppCtx *userctx = (AppCtx *)ctx;
486: PetscScalar value;
487: const PetscInt *indices;
490: SNESGetDM(snes,&dm);
491: DMGetDS(dm,&prob);
492: /* hook functions */
493: PetscDSSetResidual(prob,0,f0_u,NULL);
494: FormFunction(snes,X,F,ctx);
495: /* Boundary condition */
496: VecGetOwnershipRange(F,&iStart,&iEnd);
497: ISGetLocalSize(userctx->bdis,&nindices);
498: ISGetIndices(userctx->bdis,&indices);
499: for (i=0;i<nindices;i++) {
500: value = 0.0;
501: VecSetValue(F,indices[i],value,INSERT_VALUES);
502: }
503: ISRestoreIndices(userctx->bdis,&indices);
504: VecAssemblyBegin(F);
505: VecAssemblyEnd(F);
506: return(0);
507: }
509: PetscErrorCode MatMult_B(Mat B,Vec x,Vec y)
510: {
512: AppCtx *userctx;
515: MatShellGetContext(B,&userctx);
516: FormFunctionB(userctx->snes,x,y,userctx);
517: return(0);
518: }
520: /*TEST
522: testset:
523: requires: double
524: args: -petscspace_degree 1 -petscspace_poly_tensor
525: output_file: output/ex34_1.out
526: test:
527: suffix: 1
528: test:
529: suffix: 2
530: args: -eps_power_update -form_function_ab {{0 1}}
531: filter: sed -e "s/ with monolithic update//"
532: test:
533: suffix: 3
534: args: -use_shell_matrix -eps_power_snes_mf_operator 1
535: test:
536: suffix: 4
537: args: -use_shell_matrix -eps_power_update -init_eps_power_snes_mf_operator 1 -eps_power_snes_mf_operator 1 -form_function_ab {{0 1}}
538: filter: sed -e "s/ with monolithic update//"
539: test:
540: suffix: 5
541: args: -use_shell_matrix -eps_power_update -init_eps_power_snes_mf_operator 1 -eps_power_snes_mf_operator 1 -form_function_ab {{0 1}} -test_init_sol 1
542: filter: sed -e "s/ with monolithic update//"
544: test:
545: suffix: 6
546: args: -use_shell_matrix -eps_power_update -init_eps_power_snes_mf_operator 1 -eps_power_snes_mf_operator 1 -form_function_ab {{0 1}} -eps_monitor_all
547: output_file: output/ex34_6.out
548: filter: sed -e "s/\([+-].*i\)//g" -e "1,3s/[0-9]//g" -e "/[45] EPS/d"
549: TEST*/