Actual source code: ex34.c

slepc-main 2024-11-09
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*);
 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*/