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