Actual source code: ex34.c
slepc-3.22.1 2024-10-28
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*);
39: PetscErrorCode FormNorm(SNES,Vec,PetscReal*,void*);
41: typedef struct {
42: IS bdis; /* global indices for boundary DoFs */
43: SNES snes;
44: EPS eps;
45: } AppCtx;
47: int main(int argc,char **argv)
48: {
49: DM dm;
50: MPI_Comm comm;
51: AppCtx user;
52: EPS eps; /* eigenproblem solver context */
53: ST st;
54: EPSType type;
55: Mat A,B,P;
56: Vec v0;
57: PetscContainer container;
58: PetscInt nev,nconv,m,n,M,N;
59: PetscBool nonlin,flg=PETSC_FALSE,update;
60: SNES snes;
61: PetscReal tol,relerr;
62: PetscBool use_shell_matrix=PETSC_FALSE,test_init_sol=PETSC_FALSE,use_custom_norm=PETSC_FALSE,sign_normalization=PETSC_TRUE;
64: PetscFunctionBeginUser;
65: PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
66: comm = PETSC_COMM_WORLD;
67: /* Create a quadrilateral mesh on domain (0,1)x(0,1) */
68: PetscCall(CreateSquareMesh(comm,&dm));
69: /* Setup basis function */
70: PetscCall(SetupDiscretization(dm));
71: PetscCall(BoundaryGlobalIndex(dm,"marker",&user.bdis));
72: /* Check if we are going to use shell matrices */
73: PetscCall(PetscOptionsGetBool(NULL,NULL,"-use_shell_matrix",&use_shell_matrix,NULL));
74: if (use_shell_matrix) {
75: PetscCall(DMCreateMatrix(dm,&P));
76: PetscCall(MatGetLocalSize(P,&m,&n));
77: PetscCall(MatGetSize(P,&M,&N));
78: PetscCall(MatCreateShell(comm,m,n,M,N,&user,&A));
79: PetscCall(MatShellSetOperation(A,MATOP_MULT,(void(*)(void))MatMult_A));
80: PetscCall(MatCreateShell(comm,m,n,M,N,&user,&B));
81: PetscCall(MatShellSetOperation(B,MATOP_MULT,(void(*)(void))MatMult_B));
82: } else {
83: PetscCall(DMCreateMatrix(dm,&A));
84: PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&B));
85: }
86: /* Check whether we should use a custom normalization */
87: PetscCall(PetscOptionsGetBool(NULL,NULL,"-use_custom_norm",&use_custom_norm,NULL));
88: /* Check whether we should normalize Bx by the sign of its first nonzero element */
89: PetscCall(PetscOptionsGetBool(NULL,NULL,"-sign_normalization",&sign_normalization,NULL));
91: /*
92: Compose callback functions and context that will be needed by the solver
93: */
94: PetscCall(PetscObjectComposeFunction((PetscObject)A,"formFunction",FormFunctionA));
95: PetscCall(PetscOptionsGetBool(NULL,NULL,"-form_function_ab",&flg,NULL));
96: if (flg) PetscCall(PetscObjectComposeFunction((PetscObject)A,"formFunctionAB",FormFunctionAB));
97: PetscCall(PetscObjectComposeFunction((PetscObject)A,"formJacobian",FormJacobianA));
98: PetscCall(PetscObjectComposeFunction((PetscObject)B,"formFunction",FormFunctionB));
99: if (use_custom_norm) PetscCall(PetscObjectComposeFunction((PetscObject)B,"formNorm",FormNorm));
100: PetscCall(PetscContainerCreate(comm,&container));
101: PetscCall(PetscContainerSetPointer(container,&user));
102: PetscCall(PetscObjectCompose((PetscObject)A,"formFunctionCtx",(PetscObject)container));
103: PetscCall(PetscObjectCompose((PetscObject)A,"formJacobianCtx",(PetscObject)container));
104: PetscCall(PetscObjectCompose((PetscObject)B,"formFunctionCtx",(PetscObject)container));
105: if (use_custom_norm) PetscCall(PetscObjectCompose((PetscObject)B,"formNormCtx",(PetscObject)container));
106: PetscCall(PetscContainerDestroy(&container));
108: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109: Create the eigensolver and set various options
110: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
112: PetscCall(EPSCreate(comm,&eps));
113: PetscCall(EPSSetOperators(eps,A,B));
114: PetscCall(EPSSetProblemType(eps,EPS_GNHEP));
115: user.eps = eps;
116: /*
117: Use nonlinear inverse iteration
118: */
119: PetscCall(EPSSetType(eps,EPSPOWER));
120: PetscCall(EPSPowerSetNonlinear(eps,PETSC_TRUE));
121: /* Set the Bx sign normalization (or not) */
122: PetscCall(EPSPowerSetSignNormalization(eps,sign_normalization));
123: /*
124: Attach DM to SNES
125: */
126: PetscCall(EPSPowerGetSNES(eps,&snes));
127: user.snes = snes;
128: PetscCall(SNESSetDM(snes,dm));
129: PetscCall(EPSSetFromOptions(eps));
131: /* Set a preconditioning matrix to ST */
132: if (use_shell_matrix) {
133: PetscCall(EPSGetST(eps,&st));
134: PetscCall(STSetPreconditionerMat(st,P));
135: }
137: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138: Solve the eigensystem
139: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
141: PetscCall(EPSSolve(eps));
143: PetscCall(EPSGetConverged(eps,&nconv));
144: PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_init_sol",&test_init_sol,NULL));
145: if (nconv && test_init_sol) {
146: PetscScalar k;
147: PetscReal norm0;
148: PetscInt nits;
150: PetscCall(MatCreateVecs(A,&v0,NULL));
151: PetscCall(EPSGetEigenpair(eps,0,&k,NULL,v0,NULL));
152: PetscCall(EPSSetInitialSpace(eps,1,&v0));
153: PetscCall(VecDestroy(&v0));
154: /* Norm of the previous residual */
155: PetscCall(SNESGetFunctionNorm(snes,&norm0));
156: /* Make the tolerance smaller than the last residual
157: SNES will converge right away if the initial is setup correctly */
158: PetscCall(SNESSetTolerances(snes,norm0*1.2,PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
159: PetscCall(EPSSolve(eps));
160: /* Number of Newton iterations supposes to be zero */
161: PetscCall(SNESGetIterationNumber(snes,&nits));
162: if (nits) PetscCall(PetscPrintf(comm," Number of Newton iterations %" PetscInt_FMT " should be zero \n",nits));
163: }
165: /*
166: Optional: Get some information from the solver and display it
167: */
168: PetscCall(EPSGetType(eps,&type));
169: PetscCall(EPSGetTolerances(eps,&tol,NULL));
170: PetscCall(EPSPowerGetNonlinear(eps,&nonlin));
171: PetscCall(EPSPowerGetUpdate(eps,&update));
172: PetscCall(PetscPrintf(comm," Solution method: %s%s\n\n",type,nonlin?(update?" (nonlinear with monolithic update)":" (nonlinear)"):""));
173: PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
174: PetscCall(PetscPrintf(comm," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
176: /* print eigenvalue and error */
177: PetscCall(EPSGetConverged(eps,&nconv));
178: if (nconv>0) {
179: PetscScalar k;
180: PetscReal na,nb;
181: Vec a,b,eigen;
182: PetscCall(DMCreateGlobalVector(dm,&a));
183: PetscCall(VecDuplicate(a,&b));
184: PetscCall(VecDuplicate(a,&eigen));
185: PetscCall(EPSGetEigenpair(eps,0,&k,NULL,eigen,NULL));
186: PetscCall(FormFunctionA(snes,eigen,a,&user));
187: PetscCall(FormFunctionB(snes,eigen,b,&user));
188: PetscCall(VecAXPY(a,-k,b));
189: PetscCall(VecNorm(a,NORM_2,&na));
190: if (use_custom_norm) PetscCall(FormNorm(snes,b,&nb,&user));
191: else PetscCall(VecNorm(b,NORM_2,&nb));
192: relerr = na/(nb*PetscAbsScalar(k));
193: if (relerr<10*tol) PetscCall(PetscPrintf(comm,"k: %g, relative error below tol\n",(double)PetscRealPart(k)));
194: else PetscCall(PetscPrintf(comm,"k: %g, relative error: %g\n",(double)PetscRealPart(k),(double)relerr));
195: PetscCall(VecDestroy(&a));
196: PetscCall(VecDestroy(&b));
197: PetscCall(VecDestroy(&eigen));
198: } else PetscCall(PetscPrintf(comm,"Solver did not converge\n"));
200: PetscCall(MatDestroy(&A));
201: PetscCall(MatDestroy(&B));
202: if (use_shell_matrix) PetscCall(MatDestroy(&P));
203: PetscCall(DMDestroy(&dm));
204: PetscCall(EPSDestroy(&eps));
205: PetscCall(ISDestroy(&user.bdis));
206: PetscCall(SlepcFinalize());
207: return 0;
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;
261: PetscFunctionBeginUser;
262: /* Create finite element */
263: PetscCall(PetscObjectGetComm((PetscObject)dm,&comm));
264: PetscCall(PetscFECreateDefault(comm,2,1,PETSC_FALSE,NULL,-1,&fe));
265: PetscCall(PetscObjectSetName((PetscObject)fe,"u"));
266: PetscCall(DMSetField(dm,0,NULL,(PetscObject)fe));
267: PetscCall(DMCreateDS(dm));
268: PetscCall(PetscFEDestroy(&fe));
269: PetscFunctionReturn(PETSC_SUCCESS);
270: }
272: PetscErrorCode CreateSquareMesh(MPI_Comm comm,DM *dm)
273: {
274: PetscInt cells[] = {8,8};
275: PetscInt dim = 2;
276: DM pdm;
277: PetscMPIInt size;
279: PetscFunctionBegin;
280: PetscCall(DMPlexCreateBoxMesh(comm,dim,PETSC_FALSE,cells,NULL,NULL,NULL,PETSC_TRUE,0,PETSC_TRUE,dm));
281: PetscCall(DMSetFromOptions(*dm));
282: PetscCall(DMSetUp(*dm));
283: PetscCallMPI(MPI_Comm_size(comm,&size));
284: if (size > 1) {
285: PetscCall(DMPlexDistribute(*dm,0,NULL,&pdm));
286: PetscCall(DMDestroy(dm));
287: *dm = pdm;
288: }
289: PetscFunctionReturn(PETSC_SUCCESS);
290: }
292: PetscErrorCode BoundaryGlobalIndex(DM dm,const char labelname[],IS *bdis)
293: {
294: IS bdpoints;
295: PetscInt nindices,*indices,numDof,offset,npoints,i,j;
296: const PetscInt *bdpoints_indices;
297: DMLabel bdmarker;
298: PetscSection gsection;
300: PetscFunctionBegin;
301: PetscCall(DMGetGlobalSection(dm,&gsection));
302: PetscCall(DMGetLabel(dm,labelname,&bdmarker));
303: PetscCall(DMLabelGetStratumIS(bdmarker,1,&bdpoints));
304: PetscCall(ISGetLocalSize(bdpoints,&npoints));
305: PetscCall(ISGetIndices(bdpoints,&bdpoints_indices));
306: nindices = 0;
307: for (i=0;i<npoints;i++) {
308: PetscCall(PetscSectionGetDof(gsection,bdpoints_indices[i],&numDof));
309: if (numDof<=0) continue;
310: nindices += numDof;
311: }
312: PetscCall(PetscCalloc1(nindices,&indices));
313: nindices = 0;
314: for (i=0;i<npoints;i++) {
315: PetscCall(PetscSectionGetDof(gsection,bdpoints_indices[i],&numDof));
316: if (numDof<=0) continue;
317: PetscCall(PetscSectionGetOffset(gsection,bdpoints_indices[i],&offset));
318: for (j=0;j<numDof;j++) indices[nindices++] = offset+j;
319: }
320: PetscCall(ISRestoreIndices(bdpoints,&bdpoints_indices));
321: PetscCall(ISDestroy(&bdpoints));
322: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm),nindices,indices,PETSC_OWN_POINTER,bdis));
323: PetscFunctionReturn(PETSC_SUCCESS);
324: }
326: static PetscErrorCode FormJacobian(SNES snes,Vec X,Mat A,Mat B,void *ctx)
327: {
328: DM dm;
329: Vec Xloc;
331: PetscFunctionBegin;
332: PetscCall(SNESGetDM(snes,&dm));
333: PetscCall(DMGetLocalVector(dm,&Xloc));
334: PetscCall(VecZeroEntries(Xloc));
335: PetscCall(DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc));
336: PetscCall(DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc));
337: CHKMEMQ;
338: PetscCall(DMPlexSNESComputeJacobianFEM(dm,Xloc,A,B,ctx));
339: if (A!=B) {
340: PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
341: PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
342: }
343: CHKMEMQ;
344: PetscCall(DMRestoreLocalVector(dm,&Xloc));
345: PetscFunctionReturn(PETSC_SUCCESS);
346: }
348: PetscErrorCode FormJacobianA(SNES snes,Vec X,Mat A,Mat B,void *ctx)
349: {
350: DM dm;
351: PetscDS prob;
352: PetscWeakForm wf;
353: AppCtx *userctx = (AppCtx *)ctx;
355: PetscFunctionBegin;
356: PetscCall(MatSetOption(B,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE));
357: PetscCall(SNESGetDM(snes,&dm));
358: PetscCall(DMGetDS(dm,&prob));
359: PetscCall(PetscDSGetWeakForm(prob, &wf));
360: PetscCall(PetscWeakFormClearIndex(wf, NULL, 0, 0, 0, PETSC_WF_G3, 0));
361: PetscCall(PetscWeakFormSetIndexJacobian(wf, NULL, 0, 0, 0, 0, 0, NULL, 0, NULL, 0, NULL, 0, g3_uu));
362: PetscCall(FormJacobian(snes,X,A,B,ctx));
363: PetscCall(MatZeroRowsIS(B,userctx->bdis,1.0,NULL,NULL));
364: PetscFunctionReturn(PETSC_SUCCESS);
365: }
367: PetscErrorCode FormJacobianB(SNES snes,Vec X,Mat A,Mat B,void *ctx)
368: {
369: DM dm;
370: PetscDS prob;
371: PetscWeakForm wf;
372: AppCtx *userctx = (AppCtx *)ctx;
374: PetscFunctionBegin;
375: PetscCall(MatSetOption(B,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE));
376: PetscCall(SNESGetDM(snes,&dm));
377: PetscCall(DMGetDS(dm,&prob));
378: PetscCall(PetscDSGetWeakForm(prob, &wf));
379: PetscCall(PetscWeakFormClearIndex(wf, NULL, 0, 0, 0, PETSC_WF_G3, 0));
380: PetscCall(PetscWeakFormSetIndexJacobian(wf, NULL, 0, 0, 0, 0, 0, g0_uu, 0, NULL, 0, NULL, 0, NULL));
381: PetscCall(FormJacobian(snes,X,A,B,ctx));
382: PetscCall(MatZeroRowsIS(B,userctx->bdis,0.0,NULL,NULL));
383: PetscFunctionReturn(PETSC_SUCCESS);
384: }
386: PetscErrorCode FormFunctionAB(SNES snes,Vec x,Vec Ax,Vec Bx,void *ctx)
387: {
388: PetscFunctionBegin;
389: /*
390: * In real applications, users should have a generic formFunctionAB which
391: * forms Ax and Bx simultaneously for an more efficient calculation.
392: * In this example, we just call FormFunctionA+FormFunctionB to mimic how
393: * to use FormFunctionAB
394: */
395: PetscCall(FormFunctionA(snes,x,Ax,ctx));
396: PetscCall(FormFunctionB(snes,x,Bx,ctx));
397: PetscFunctionReturn(PETSC_SUCCESS);
398: }
400: static PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ctx)
401: {
402: DM dm;
403: Vec Xloc,Floc;
405: PetscFunctionBegin;
406: PetscCall(SNESGetDM(snes,&dm));
407: PetscCall(DMGetLocalVector(dm,&Xloc));
408: PetscCall(DMGetLocalVector(dm,&Floc));
409: PetscCall(VecZeroEntries(Xloc));
410: PetscCall(VecZeroEntries(Floc));
411: PetscCall(DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc));
412: PetscCall(DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc));
413: CHKMEMQ;
414: PetscCall(DMPlexSNESComputeResidualFEM(dm,Xloc,Floc,ctx));
415: CHKMEMQ;
416: PetscCall(VecZeroEntries(F));
417: PetscCall(DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F));
418: PetscCall(DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F));
419: PetscCall(DMRestoreLocalVector(dm,&Xloc));
420: PetscCall(DMRestoreLocalVector(dm,&Floc));
421: PetscFunctionReturn(PETSC_SUCCESS);
422: }
424: PetscErrorCode FormFunctionA(SNES snes,Vec X,Vec F,void *ctx)
425: {
426: DM dm;
427: PetscDS prob;
428: PetscWeakForm wf;
429: PetscInt nindices,iStart,iEnd,i;
430: AppCtx *userctx = (AppCtx *)ctx;
431: PetscScalar *array,value;
432: const PetscInt *indices;
433: PetscInt vecstate;
435: PetscFunctionBegin;
436: PetscCall(SNESGetDM(snes,&dm));
437: PetscCall(DMGetDS(dm,&prob));
438: /* hook functions */
439: PetscCall(PetscDSGetWeakForm(prob, &wf));
440: PetscCall(PetscWeakFormClearIndex(wf, NULL, 0, 0, 0, PETSC_WF_F0, 0));
441: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, 0, 0, 0, NULL, 0, f1_u));
442: PetscCall(FormFunction(snes,X,F,ctx));
443: /* Boundary condition */
444: PetscCall(VecLockGet(X,&vecstate));
445: if (vecstate>0) PetscCall(VecLockReadPop(X));
446: PetscCall(VecGetOwnershipRange(X,&iStart,&iEnd));
447: PetscCall(VecGetArray(X,&array));
448: PetscCall(ISGetLocalSize(userctx->bdis,&nindices));
449: PetscCall(ISGetIndices(userctx->bdis,&indices));
450: for (i=0;i<nindices;i++) {
451: value = array[indices[i]-iStart] - 0.0;
452: PetscCall(VecSetValue(F,indices[i],value,INSERT_VALUES));
453: }
454: PetscCall(ISRestoreIndices(userctx->bdis,&indices));
455: PetscCall(VecRestoreArray(X,&array));
456: if (vecstate>0) PetscCall(VecLockReadPush(X));
457: PetscCall(VecAssemblyBegin(F));
458: PetscCall(VecAssemblyEnd(F));
459: PetscFunctionReturn(PETSC_SUCCESS);
460: }
462: PetscErrorCode FormNorm(SNES snes,Vec Bx,PetscReal *norm,void *ctx)
463: {
464: PetscFunctionBegin;
465: PetscCall(VecNorm(Bx,NORM_2,norm));
466: PetscFunctionReturn(PETSC_SUCCESS);
467: }
469: PetscErrorCode MatMult_A(Mat A,Vec x,Vec y)
470: {
471: AppCtx *userctx;
473: PetscFunctionBegin;
474: PetscCall(MatShellGetContext(A,&userctx));
475: PetscCall(FormFunctionA(userctx->snes,x,y,userctx));
476: PetscFunctionReturn(PETSC_SUCCESS);
477: }
479: PetscErrorCode FormFunctionB(SNES snes,Vec X,Vec F,void *ctx)
480: {
481: DM dm;
482: PetscDS prob;
483: PetscWeakForm wf;
484: PetscInt nindices,iStart,iEnd,i;
485: AppCtx *userctx = (AppCtx *)ctx;
486: PetscScalar value;
487: const PetscInt *indices;
489: PetscFunctionBegin;
490: PetscCall(SNESGetDM(snes,&dm));
491: PetscCall(DMGetDS(dm,&prob));
492: /* hook functions */
493: PetscCall(PetscDSGetWeakForm(prob, &wf));
494: PetscCall(PetscWeakFormClearIndex(wf, NULL, 0, 0, 0, PETSC_WF_F1, 0));
495: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, 0, 0, 0, f0_u, 0, NULL));
496: PetscCall(FormFunction(snes,X,F,ctx));
497: /* Boundary condition */
498: PetscCall(VecGetOwnershipRange(F,&iStart,&iEnd));
499: PetscCall(ISGetLocalSize(userctx->bdis,&nindices));
500: PetscCall(ISGetIndices(userctx->bdis,&indices));
501: for (i=0;i<nindices;i++) {
502: value = 0.0;
503: PetscCall(VecSetValue(F,indices[i],value,INSERT_VALUES));
504: }
505: PetscCall(ISRestoreIndices(userctx->bdis,&indices));
506: PetscCall(VecAssemblyBegin(F));
507: PetscCall(VecAssemblyEnd(F));
508: PetscFunctionReturn(PETSC_SUCCESS);
509: }
511: PetscErrorCode MatMult_B(Mat B,Vec x,Vec y)
512: {
513: AppCtx *userctx;
515: PetscFunctionBegin;
516: PetscCall(MatShellGetContext(B,&userctx));
517: PetscCall(FormFunctionB(userctx->snes,x,y,userctx));
518: PetscFunctionReturn(PETSC_SUCCESS);
519: }
521: /*TEST
523: testset:
524: requires: double
525: args: -petscspace_degree 1 -petscspace_poly_tensor -checkfunctionlist 0
526: output_file: output/ex34_1.out
527: test:
528: suffix: 1
529: test:
530: suffix: 2
531: args: -eps_power_update -form_function_ab {{0 1}}
532: filter: sed -e "s/ with monolithic update//"
533: test:
534: suffix: 3
535: args: -use_shell_matrix -eps_power_snes_mf_operator 1
536: test:
537: suffix: 4
538: 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}}
539: filter: sed -e "s/ with monolithic update//"
540: test:
541: suffix: 5
542: 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
543: filter: sed -e "s/ with monolithic update//"
545: test:
546: suffix: 6
547: 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
548: output_file: output/ex34_6.out
549: filter: sed -e "s/\([+-].*i\)//g" -e "1,3s/[0-9]//g" -e "/[45] EPS/d"
550: test:
551: suffix: 7
552: args: -use_custom_norm -sign_normalization 1 -eps_power_snes_mf_operator 1
553: test:
554: suffix: 8
555: args: -use_custom_norm -sign_normalization 1 -eps_power_update -form_function_ab {{0 1}} -eps_power_snes_mf_operator 1 -init_eps_power_snes_mf_operator 1
556: filter: sed -e "s/ with monolithic update//"
557: test:
558: suffix: 9
559: requires: !complex
560: args: -use_custom_norm {{0 1}} -sign_normalization 0 -eps_power_snes_mf_operator 1
561: test:
562: suffix: 10
563: requires: !complex
564: args: -use_custom_norm {{0 1}} -sign_normalization 0 -eps_power_update -form_function_ab {{0 1}} -eps_power_snes_mf_operator 1 -init_eps_power_snes_mf_operator 1
565: filter: sed -e "s/ with monolithic update//"
566: test:
567: suffix: 11
568: requires: complex
569: args: -use_custom_norm {{0 1}} -sign_normalization 0 -eps_power_snes_type nrichardson -eps_power_snes_atol 1e-12
570: test:
571: suffix: 12
572: requires: complex
573: args: -use_custom_norm {{0 1}} -sign_normalization 0 -eps_power_update -init_eps_power_snes_type nrichardson -init_eps_max_it 2 -eps_power_snes_mf_operator 1
574: filter: sed -e "s/ with monolithic update//"
575: TEST*/