Actual source code: ex34.c

slepc-3.19.2 2023-09-05
Report Typos and Errors
  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:   PetscFunctionBeginUser;
 63:   PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
 64:   comm = PETSC_COMM_WORLD;
 65:   /* Create a quadrilateral mesh on domain (0,1)x(0,1) */
 66:   PetscCall(CreateSquareMesh(comm,&dm));
 67:   /* Setup basis function */
 68:   PetscCall(SetupDiscretization(dm));
 69:   PetscCall(BoundaryGlobalIndex(dm,"marker",&user.bdis));
 70:   /* Check if we are going to use shell matrices */
 71:   PetscCall(PetscOptionsGetBool(NULL,NULL,"-use_shell_matrix",&use_shell_matrix,NULL));
 72:   if (use_shell_matrix) {
 73:     PetscCall(DMCreateMatrix(dm,&P));
 74:     PetscCall(MatGetLocalSize(P,&m,&n));
 75:     PetscCall(MatGetSize(P,&M,&N));
 76:     PetscCall(MatCreateShell(comm,m,n,M,N,&user,&A));
 77:     PetscCall(MatShellSetOperation(A,MATOP_MULT,(void(*)(void))MatMult_A));
 78:     PetscCall(MatCreateShell(comm,m,n,M,N,&user,&B));
 79:     PetscCall(MatShellSetOperation(B,MATOP_MULT,(void(*)(void))MatMult_B));
 80:   } else {
 81:     PetscCall(DMCreateMatrix(dm,&A));
 82:     PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&B));
 83:   }

 85:   /*
 86:      Compose callback functions and context that will be needed by the solver
 87:   */
 88:   PetscCall(PetscObjectComposeFunction((PetscObject)A,"formFunction",FormFunctionA));
 89:   PetscCall(PetscOptionsGetBool(NULL,NULL,"-form_function_ab",&flg,NULL));
 90:   if (flg) PetscCall(PetscObjectComposeFunction((PetscObject)A,"formFunctionAB",FormFunctionAB));
 91:   PetscCall(PetscObjectComposeFunction((PetscObject)A,"formJacobian",FormJacobianA));
 92:   PetscCall(PetscObjectComposeFunction((PetscObject)B,"formFunction",FormFunctionB));
 93:   PetscCall(PetscContainerCreate(comm,&container));
 94:   PetscCall(PetscContainerSetPointer(container,&user));
 95:   PetscCall(PetscObjectCompose((PetscObject)A,"formFunctionCtx",(PetscObject)container));
 96:   PetscCall(PetscObjectCompose((PetscObject)A,"formJacobianCtx",(PetscObject)container));
 97:   PetscCall(PetscObjectCompose((PetscObject)B,"formFunctionCtx",(PetscObject)container));
 98:   PetscCall(PetscContainerDestroy(&container));

100:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101:                 Create the eigensolver and set various options
102:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

104:   PetscCall(EPSCreate(comm,&eps));
105:   PetscCall(EPSSetOperators(eps,A,B));
106:   PetscCall(EPSSetProblemType(eps,EPS_GNHEP));
107:   /*
108:      Use nonlinear inverse iteration
109:   */
110:   PetscCall(EPSSetType(eps,EPSPOWER));
111:   PetscCall(EPSPowerSetNonlinear(eps,PETSC_TRUE));
112:   /*
113:     Attach DM to SNES
114:   */
115:   PetscCall(EPSPowerGetSNES(eps,&snes));
116:   user.snes = snes;
117:   PetscCall(SNESSetDM(snes,dm));
118:   PetscCall(EPSSetFromOptions(eps));

120:   /* Set a preconditioning matrix to ST */
121:   if (use_shell_matrix) {
122:     PetscCall(EPSGetST(eps,&st));
123:     PetscCall(STSetPreconditionerMat(st,P));
124:   }

126:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127:                       Solve the eigensystem
128:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

130:   PetscCall(EPSSolve(eps));

132:   PetscCall(EPSGetConverged(eps,&nconv));
133:   PetscCall(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:     PetscCall(MatCreateVecs(A,&v0,NULL));
140:     PetscCall(EPSGetEigenpair(eps,0,&k,NULL,v0,NULL));
141:     PetscCall(EPSSetInitialSpace(eps,1,&v0));
142:     PetscCall(VecDestroy(&v0));
143:     /* Norm of the previous residual */
144:     PetscCall(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:     PetscCall(SNESSetTolerances(snes,norm0*1.2,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT));
148:     PetscCall(EPSSolve(eps));
149:     /* Number of Newton iterations supposes to be zero */
150:     PetscCall(SNESGetIterationNumber(snes,&nits));
151:     if (nits) PetscCall(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:   PetscCall(EPSGetType(eps,&type));
158:   PetscCall(EPSGetTolerances(eps,&tol,NULL));
159:   PetscCall(EPSPowerGetNonlinear(eps,&nonlin));
160:   PetscCall(EPSPowerGetUpdate(eps,&update));
161:   PetscCall(PetscPrintf(comm," Solution method: %s%s\n\n",type,nonlin?(update?" (nonlinear with monolithic update)":" (nonlinear)"):""));
162:   PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
163:   PetscCall(PetscPrintf(comm," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));

165:   /* print eigenvalue and error */
166:   PetscCall(EPSGetConverged(eps,&nconv));
167:   if (nconv>0) {
168:     PetscScalar   k;
169:     PetscReal     na,nb;
170:     Vec           a,b,eigen;
171:     PetscCall(DMCreateGlobalVector(dm,&a));
172:     PetscCall(VecDuplicate(a,&b));
173:     PetscCall(VecDuplicate(a,&eigen));
174:     PetscCall(EPSGetEigenpair(eps,0,&k,NULL,eigen,NULL));
175:     PetscCall(FormFunctionA(snes,eigen,a,&user));
176:     PetscCall(FormFunctionB(snes,eigen,b,&user));
177:     PetscCall(VecAXPY(a,-k,b));
178:     PetscCall(VecNorm(a,NORM_2,&na));
179:     PetscCall(VecNorm(b,NORM_2,&nb));
180:     relerr = na/(nb*PetscAbsScalar(k));
181:     if (relerr<10*tol) PetscCall(PetscPrintf(comm,"k: %g, relative error below tol\n",(double)PetscRealPart(k)));
182:     else PetscCall(PetscPrintf(comm,"k: %g, relative error: %g\n",(double)PetscRealPart(k),(double)relerr));
183:     PetscCall(VecDestroy(&a));
184:     PetscCall(VecDestroy(&b));
185:     PetscCall(VecDestroy(&eigen));
186:   } else PetscCall(PetscPrintf(comm,"Solver did not converge\n"));

188:   PetscCall(MatDestroy(&A));
189:   PetscCall(MatDestroy(&B));
190:   if (use_shell_matrix) PetscCall(MatDestroy(&P));
191:   PetscCall(DMDestroy(&dm));
192:   PetscCall(EPSDestroy(&eps));
193:   PetscCall(ISDestroy(&user.bdis));
194:   PetscCall(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;

249:   PetscFunctionBeginUser;
250:   /* Create finite element */
251:   PetscCall(PetscObjectGetComm((PetscObject)dm,&comm));
252:   PetscCall(PetscFECreateDefault(comm,2,1,PETSC_FALSE,NULL,-1,&fe));
253:   PetscCall(PetscObjectSetName((PetscObject)fe,"u"));
254:   PetscCall(DMSetField(dm,0,NULL,(PetscObject)fe));
255:   PetscCall(DMCreateDS(dm));
256:   PetscCall(PetscFEDestroy(&fe));
257:   PetscFunctionReturn(PETSC_SUCCESS);
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:   PetscFunctionBegin;
268:   PetscCall(DMPlexCreateBoxMesh(comm,dim,PETSC_FALSE,cells,NULL,NULL,NULL,PETSC_TRUE,dm));
269:   PetscCall(DMSetFromOptions(*dm));
270:   PetscCall(DMSetUp(*dm));
271:   PetscCallMPI(MPI_Comm_size(comm,&size));
272:   if (size > 1) {
273:     PetscCall(DMPlexDistribute(*dm,0,NULL,&pdm));
274:     PetscCall(DMDestroy(dm));
275:     *dm = pdm;
276:   }
277:   PetscFunctionReturn(PETSC_SUCCESS);
278: }

280: PetscErrorCode BoundaryGlobalIndex(DM dm,const char labelname[],IS *bdis)
281: {
282:   IS             bdpoints;
283:   PetscInt       nindices,*indices,numDof,offset,npoints,i,j;
284:   const PetscInt *bdpoints_indices;
285:   DMLabel        bdmarker;
286:   PetscSection   gsection;

288:   PetscFunctionBegin;
289:   PetscCall(DMGetGlobalSection(dm,&gsection));
290:   PetscCall(DMGetLabel(dm,labelname,&bdmarker));
291:   PetscCall(DMLabelGetStratumIS(bdmarker,1,&bdpoints));
292:   PetscCall(ISGetLocalSize(bdpoints,&npoints));
293:   PetscCall(ISGetIndices(bdpoints,&bdpoints_indices));
294:   nindices = 0;
295:   for (i=0;i<npoints;i++) {
296:     PetscCall(PetscSectionGetDof(gsection,bdpoints_indices[i],&numDof));
297:     if (numDof<=0) continue;
298:     nindices += numDof;
299:   }
300:   PetscCall(PetscCalloc1(nindices,&indices));
301:   nindices = 0;
302:   for (i=0;i<npoints;i++) {
303:     PetscCall(PetscSectionGetDof(gsection,bdpoints_indices[i],&numDof));
304:     if (numDof<=0) continue;
305:     PetscCall(PetscSectionGetOffset(gsection,bdpoints_indices[i],&offset));
306:     for (j=0;j<numDof;j++) indices[nindices++] = offset+j;
307:   }
308:   PetscCall(ISRestoreIndices(bdpoints,&bdpoints_indices));
309:   PetscCall(ISDestroy(&bdpoints));
310:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm),nindices,indices,PETSC_OWN_POINTER,bdis));
311:   PetscFunctionReturn(PETSC_SUCCESS);
312: }

314: static PetscErrorCode FormJacobian(SNES snes,Vec X,Mat A,Mat B,void *ctx)
315: {
316:   DM             dm;
317:   Vec            Xloc;

319:   PetscFunctionBegin;
320:   PetscCall(SNESGetDM(snes,&dm));
321:   PetscCall(DMGetLocalVector(dm,&Xloc));
322:   PetscCall(VecZeroEntries(Xloc));
323:   PetscCall(DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc));
324:   PetscCall(DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc));
325:   CHKMEMQ;
326:   PetscCall(DMPlexSNESComputeJacobianFEM(dm,Xloc,A,B,ctx));
327:   if (A!=B) {
328:     PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
329:     PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
330:   }
331:   CHKMEMQ;
332:   PetscCall(DMRestoreLocalVector(dm,&Xloc));
333:   PetscFunctionReturn(PETSC_SUCCESS);
334: }

336: PetscErrorCode FormJacobianA(SNES snes,Vec X,Mat A,Mat B,void *ctx)
337: {
338:   DM             dm;
339:   PetscDS        prob;
340:   PetscWeakForm  wf;
341:   AppCtx         *userctx = (AppCtx *)ctx;

343:   PetscFunctionBegin;
344:   PetscCall(MatSetOption(B,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE));
345:   PetscCall(SNESGetDM(snes,&dm));
346:   PetscCall(DMGetDS(dm,&prob));
347:   PetscCall(PetscDSGetWeakForm(prob, &wf));
348:   PetscCall(PetscWeakFormClearIndex(wf, NULL, 0, 0, 0, PETSC_WF_G3, 0));
349:   PetscCall(PetscWeakFormSetIndexJacobian(wf, NULL, 0, 0, 0, 0, 0, NULL, 0, NULL, 0, NULL, 0, g3_uu));
350:   PetscCall(FormJacobian(snes,X,A,B,ctx));
351:   PetscCall(MatZeroRowsIS(B,userctx->bdis,1.0,NULL,NULL));
352:   PetscFunctionReturn(PETSC_SUCCESS);
353: }

355: PetscErrorCode FormJacobianB(SNES snes,Vec X,Mat A,Mat B,void *ctx)
356: {
357:   DM             dm;
358:   PetscDS        prob;
359:   PetscWeakForm  wf;
360:   AppCtx         *userctx = (AppCtx *)ctx;

362:   PetscFunctionBegin;
363:   PetscCall(MatSetOption(B,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE));
364:   PetscCall(SNESGetDM(snes,&dm));
365:   PetscCall(DMGetDS(dm,&prob));
366:   PetscCall(PetscDSGetWeakForm(prob, &wf));
367:   PetscCall(PetscWeakFormClearIndex(wf, NULL, 0, 0, 0, PETSC_WF_G3, 0));
368:   PetscCall(PetscWeakFormSetIndexJacobian(wf, NULL, 0, 0, 0, 0, 0, g0_uu, 0, NULL, 0, NULL, 0, NULL));
369:   PetscCall(FormJacobian(snes,X,A,B,ctx));
370:   PetscCall(MatZeroRowsIS(B,userctx->bdis,0.0,NULL,NULL));
371:   PetscFunctionReturn(PETSC_SUCCESS);
372: }

374: PetscErrorCode FormFunctionAB(SNES snes,Vec x,Vec Ax,Vec Bx,void *ctx)
375: {
376:   PetscFunctionBegin;
377:   /*
378:    * In real applications, users should have a generic formFunctionAB which
379:    * forms Ax and Bx simultaneously for an more efficient calculation.
380:    * In this example, we just call FormFunctionA+FormFunctionB to mimic how
381:    * to use FormFunctionAB
382:    */
383:   PetscCall(FormFunctionA(snes,x,Ax,ctx));
384:   PetscCall(FormFunctionB(snes,x,Bx,ctx));
385:   PetscFunctionReturn(PETSC_SUCCESS);
386: }

388: static PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ctx)
389: {
390:   DM             dm;
391:   Vec            Xloc,Floc;

393:   PetscFunctionBegin;
394:   PetscCall(SNESGetDM(snes,&dm));
395:   PetscCall(DMGetLocalVector(dm,&Xloc));
396:   PetscCall(DMGetLocalVector(dm,&Floc));
397:   PetscCall(VecZeroEntries(Xloc));
398:   PetscCall(VecZeroEntries(Floc));
399:   PetscCall(DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc));
400:   PetscCall(DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc));
401:   CHKMEMQ;
402:   PetscCall(DMPlexSNESComputeResidualFEM(dm,Xloc,Floc,ctx));
403:   CHKMEMQ;
404:   PetscCall(VecZeroEntries(F));
405:   PetscCall(DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F));
406:   PetscCall(DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F));
407:   PetscCall(DMRestoreLocalVector(dm,&Xloc));
408:   PetscCall(DMRestoreLocalVector(dm,&Floc));
409:   PetscFunctionReturn(PETSC_SUCCESS);
410: }

412: PetscErrorCode FormFunctionA(SNES snes,Vec X,Vec F,void *ctx)
413: {
414:   DM             dm;
415:   PetscDS        prob;
416:   PetscWeakForm  wf;
417:   PetscInt       nindices,iStart,iEnd,i;
418:   AppCtx         *userctx = (AppCtx *)ctx;
419:   PetscScalar    *array,value;
420:   const PetscInt *indices;
421:   PetscInt       vecstate;

423:   PetscFunctionBegin;
424:   PetscCall(SNESGetDM(snes,&dm));
425:   PetscCall(DMGetDS(dm,&prob));
426:   /* hook functions */
427:   PetscCall(PetscDSGetWeakForm(prob, &wf));
428:   PetscCall(PetscWeakFormClearIndex(wf, NULL, 0, 0, 0, PETSC_WF_F0, 0));
429:   PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, 0, 0, 0, NULL, 0, f1_u));
430:   PetscCall(FormFunction(snes,X,F,ctx));
431:   /* Boundary condition */
432:   PetscCall(VecLockGet(X,&vecstate));
433:   if (vecstate>0) PetscCall(VecLockReadPop(X));
434:   PetscCall(VecGetOwnershipRange(X,&iStart,&iEnd));
435:   PetscCall(VecGetArray(X,&array));
436:   PetscCall(ISGetLocalSize(userctx->bdis,&nindices));
437:   PetscCall(ISGetIndices(userctx->bdis,&indices));
438:   for (i=0;i<nindices;i++) {
439:     value = array[indices[i]-iStart] - 0.0;
440:     PetscCall(VecSetValue(F,indices[i],value,INSERT_VALUES));
441:   }
442:   PetscCall(ISRestoreIndices(userctx->bdis,&indices));
443:   PetscCall(VecRestoreArray(X,&array));
444:   if (vecstate>0) PetscCall(VecLockReadPush(X));
445:   PetscCall(VecAssemblyBegin(F));
446:   PetscCall(VecAssemblyEnd(F));
447:   PetscFunctionReturn(PETSC_SUCCESS);
448: }

450: PetscErrorCode MatMult_A(Mat A,Vec x,Vec y)
451: {
452:   AppCtx         *userctx;

454:   PetscFunctionBegin;
455:   PetscCall(MatShellGetContext(A,&userctx));
456:   PetscCall(FormFunctionA(userctx->snes,x,y,userctx));
457:   PetscFunctionReturn(PETSC_SUCCESS);
458: }

460: PetscErrorCode FormFunctionB(SNES snes,Vec X,Vec F,void *ctx)
461: {
462:   DM             dm;
463:   PetscDS        prob;
464:   PetscWeakForm  wf;
465:   PetscInt       nindices,iStart,iEnd,i;
466:   AppCtx         *userctx = (AppCtx *)ctx;
467:   PetscScalar    value;
468:   const PetscInt *indices;

470:   PetscFunctionBegin;
471:   PetscCall(SNESGetDM(snes,&dm));
472:   PetscCall(DMGetDS(dm,&prob));
473:   /* hook functions */
474:   PetscCall(PetscDSGetWeakForm(prob, &wf));
475:   PetscCall(PetscWeakFormClearIndex(wf, NULL, 0, 0, 0, PETSC_WF_F1, 0));
476:   PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, 0, 0, 0, f0_u, 0, NULL));
477:   PetscCall(FormFunction(snes,X,F,ctx));
478:   /* Boundary condition */
479:   PetscCall(VecGetOwnershipRange(F,&iStart,&iEnd));
480:   PetscCall(ISGetLocalSize(userctx->bdis,&nindices));
481:   PetscCall(ISGetIndices(userctx->bdis,&indices));
482:   for (i=0;i<nindices;i++) {
483:     value = 0.0;
484:     PetscCall(VecSetValue(F,indices[i],value,INSERT_VALUES));
485:   }
486:   PetscCall(ISRestoreIndices(userctx->bdis,&indices));
487:   PetscCall(VecAssemblyBegin(F));
488:   PetscCall(VecAssemblyEnd(F));
489:   PetscFunctionReturn(PETSC_SUCCESS);
490: }

492: PetscErrorCode MatMult_B(Mat B,Vec x,Vec y)
493: {
494:   AppCtx         *userctx;

496:   PetscFunctionBegin;
497:   PetscCall(MatShellGetContext(B,&userctx));
498:   PetscCall(FormFunctionB(userctx->snes,x,y,userctx));
499:   PetscFunctionReturn(PETSC_SUCCESS);
500: }

502: /*TEST

504:    testset:
505:       requires: double
506:       args: -petscspace_degree 1 -petscspace_poly_tensor -checkfunctionlist 0
507:       output_file: output/ex34_1.out
508:       test:
509:          suffix: 1
510:       test:
511:          suffix: 2
512:          args: -eps_power_update -form_function_ab {{0 1}}
513:          filter: sed -e "s/ with monolithic update//"
514:       test:
515:          suffix: 3
516:          args: -use_shell_matrix -eps_power_snes_mf_operator 1
517:       test:
518:          suffix: 4
519:          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}}
520:          filter: sed -e "s/ with monolithic update//"
521:       test:
522:          suffix: 5
523:          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
524:          filter: sed -e "s/ with monolithic update//"

526:       test:
527:          suffix: 6
528:          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
529:          output_file: output/ex34_6.out
530:          filter: sed -e "s/\([+-].*i\)//g" -e "1,3s/[0-9]//g" -e "/[45] EPS/d"
531: TEST*/