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