Actual source code: dmdasnes.c

  1: #include <petscdmda.h>
  2: #include <petsc/private/dmimpl.h>
  3: #include <petsc/private/snesimpl.h>

  5: /* This structure holds the user-provided DMDA callbacks */
  6: typedef struct {
  7:   PetscErrorCode (*residuallocal)(DMDALocalInfo *, void *, void *, void *);
  8:   PetscErrorCode (*jacobianlocal)(DMDALocalInfo *, void *, Mat, Mat, void *);
  9:   PetscErrorCode (*objectivelocal)(DMDALocalInfo *, void *, PetscReal *, void *);

 11:   PetscErrorCode (*residuallocalvec)(DMDALocalInfo *, Vec, Vec, void *);
 12:   PetscErrorCode (*jacobianlocalvec)(DMDALocalInfo *, Vec, Mat, Mat, void *);
 13:   PetscErrorCode (*objectivelocalvec)(DMDALocalInfo *, Vec, PetscReal *, void *);
 14:   void      *residuallocalctx;
 15:   void      *jacobianlocalctx;
 16:   void      *objectivelocalctx;
 17:   InsertMode residuallocalimode;

 19:   /*   For Picard iteration defined locally */
 20:   PetscErrorCode (*rhsplocal)(DMDALocalInfo *, void *, void *, void *);
 21:   PetscErrorCode (*jacobianplocal)(DMDALocalInfo *, void *, Mat, Mat, void *);
 22:   void *picardlocalctx;
 23: } DMSNES_DA;

 25: static PetscErrorCode DMSNESDestroy_DMDA(DMSNES sdm)
 26: {
 27:   PetscFunctionBegin;
 28:   PetscCall(PetscFree(sdm->data));
 29:   PetscFunctionReturn(PETSC_SUCCESS);
 30: }

 32: static PetscErrorCode DMSNESDuplicate_DMDA(DMSNES oldsdm, DMSNES sdm)
 33: {
 34:   PetscFunctionBegin;
 35:   PetscCall(PetscNew((DMSNES_DA **)&sdm->data));
 36:   if (oldsdm->data) PetscCall(PetscMemcpy(sdm->data, oldsdm->data, sizeof(DMSNES_DA)));
 37:   PetscFunctionReturn(PETSC_SUCCESS);
 38: }

 40: static PetscErrorCode DMDASNESGetContext(DM dm, DMSNES sdm, DMSNES_DA **dmdasnes)
 41: {
 42:   PetscFunctionBegin;
 43:   *dmdasnes = NULL;
 44:   if (!sdm->data) {
 45:     PetscCall(PetscNew((DMSNES_DA **)&sdm->data));
 46:     sdm->ops->destroy   = DMSNESDestroy_DMDA;
 47:     sdm->ops->duplicate = DMSNESDuplicate_DMDA;
 48:   }
 49:   *dmdasnes = (DMSNES_DA *)sdm->data;
 50:   PetscFunctionReturn(PETSC_SUCCESS);
 51: }

 53: static PetscErrorCode SNESComputeFunction_DMDA(SNES snes, Vec X, Vec F, void *ctx)
 54: {
 55:   DM            dm;
 56:   DMSNES_DA    *dmdasnes = (DMSNES_DA *)ctx;
 57:   DMDALocalInfo info;
 58:   Vec           Xloc;
 59:   void         *x, *f;

 61:   PetscFunctionBegin;
 65:   PetscCheck(dmdasnes->residuallocal || dmdasnes->residuallocalvec, PetscObjectComm((PetscObject)snes), PETSC_ERR_PLIB, "Corrupt context");
 66:   PetscCall(SNESGetDM(snes, &dm));
 67:   PetscCall(DMGetLocalVector(dm, &Xloc));
 68:   PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, Xloc));
 69:   PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, Xloc));
 70:   PetscCall(DMDAGetLocalInfo(dm, &info));
 71:   switch (dmdasnes->residuallocalimode) {
 72:   case INSERT_VALUES: {
 73:     PetscCall(PetscLogEventBegin(SNES_FunctionEval, snes, X, F, 0));
 74:     if (dmdasnes->residuallocalvec) PetscCallBack("SNES DMDA local callback function", (*dmdasnes->residuallocalvec)(&info, Xloc, F, dmdasnes->residuallocalctx));
 75:     else {
 76:       PetscCall(DMDAVecGetArray(dm, Xloc, &x));
 77:       PetscCall(DMDAVecGetArray(dm, F, &f));
 78:       PetscCallBack("SNES DMDA local callback function", (*dmdasnes->residuallocal)(&info, x, f, dmdasnes->residuallocalctx));
 79:       PetscCall(DMDAVecRestoreArray(dm, Xloc, &x));
 80:       PetscCall(DMDAVecRestoreArray(dm, F, &f));
 81:     }
 82:     PetscCall(PetscLogEventEnd(SNES_FunctionEval, snes, X, F, 0));
 83:   } break;
 84:   case ADD_VALUES: {
 85:     Vec Floc;
 86:     PetscCall(DMGetLocalVector(dm, &Floc));
 87:     PetscCall(VecZeroEntries(Floc));
 88:     PetscCall(PetscLogEventBegin(SNES_FunctionEval, snes, X, F, 0));
 89:     if (dmdasnes->residuallocalvec) PetscCallBack("SNES DMDA local callback function", (*dmdasnes->residuallocalvec)(&info, Xloc, Floc, dmdasnes->residuallocalctx));
 90:     else {
 91:       PetscCall(DMDAVecGetArray(dm, Xloc, &x));
 92:       PetscCall(DMDAVecGetArray(dm, Floc, &f));
 93:       PetscCallBack("SNES DMDA local callback function", (*dmdasnes->residuallocal)(&info, x, f, dmdasnes->residuallocalctx));
 94:       PetscCall(DMDAVecRestoreArray(dm, Xloc, &x));
 95:       PetscCall(DMDAVecRestoreArray(dm, Floc, &f));
 96:     }
 97:     PetscCall(PetscLogEventEnd(SNES_FunctionEval, snes, X, F, 0));
 98:     PetscCall(VecZeroEntries(F));
 99:     PetscCall(DMLocalToGlobalBegin(dm, Floc, ADD_VALUES, F));
100:     PetscCall(DMLocalToGlobalEnd(dm, Floc, ADD_VALUES, F));
101:     PetscCall(DMRestoreLocalVector(dm, &Floc));
102:   } break;
103:   default:
104:     SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_INCOMP, "Cannot use imode=%d", (int)dmdasnes->residuallocalimode);
105:   }
106:   PetscCall(DMRestoreLocalVector(dm, &Xloc));
107:   if (snes->domainerror) PetscCall(VecSetInf(F));
108:   PetscFunctionReturn(PETSC_SUCCESS);
109: }

111: static PetscErrorCode SNESComputeObjective_DMDA(SNES snes, Vec X, PetscReal *ob, void *ctx)
112: {
113:   DM            dm;
114:   DMSNES_DA    *dmdasnes = (DMSNES_DA *)ctx;
115:   DMDALocalInfo info;
116:   Vec           Xloc;
117:   void         *x;

119:   PetscFunctionBegin;
122:   PetscAssertPointer(ob, 3);
123:   PetscCheck(dmdasnes->objectivelocal || dmdasnes->objectivelocalvec, PetscObjectComm((PetscObject)snes), PETSC_ERR_PLIB, "Corrupt context");
124:   PetscCall(SNESGetDM(snes, &dm));
125:   PetscCall(DMGetLocalVector(dm, &Xloc));
126:   PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, Xloc));
127:   PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, Xloc));
128:   PetscCall(DMDAGetLocalInfo(dm, &info));
129:   if (dmdasnes->objectivelocalvec) PetscCallBack("SNES DMDA local callback objective", (*dmdasnes->objectivelocalvec)(&info, Xloc, ob, dmdasnes->objectivelocalctx));
130:   else {
131:     PetscCall(DMDAVecGetArray(dm, Xloc, &x));
132:     PetscCallBack("SNES DMDA local callback objective", (*dmdasnes->objectivelocal)(&info, x, ob, dmdasnes->objectivelocalctx));
133:     PetscCall(DMDAVecRestoreArray(dm, Xloc, &x));
134:   }
135:   PetscCall(DMRestoreLocalVector(dm, &Xloc));
136:   PetscFunctionReturn(PETSC_SUCCESS);
137: }

139: /* Routine is called by example, hence must be labeled PETSC_EXTERN */
140: PETSC_EXTERN PetscErrorCode SNESComputeJacobian_DMDA(SNES snes, Vec X, Mat A, Mat B, void *ctx)
141: {
142:   DM            dm;
143:   DMSNES_DA    *dmdasnes = (DMSNES_DA *)ctx;
144:   DMDALocalInfo info;
145:   Vec           Xloc;
146:   void         *x;

148:   PetscFunctionBegin;
149:   PetscCheck(dmdasnes->residuallocal || dmdasnes->residuallocalvec, PetscObjectComm((PetscObject)snes), PETSC_ERR_PLIB, "Corrupt context");
150:   PetscCall(SNESGetDM(snes, &dm));

152:   if (dmdasnes->jacobianlocal || dmdasnes->jacobianlocalctx) {
153:     PetscCall(DMGetLocalVector(dm, &Xloc));
154:     PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, Xloc));
155:     PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, Xloc));
156:     PetscCall(DMDAGetLocalInfo(dm, &info));
157:     if (dmdasnes->jacobianlocalvec) PetscCallBack("SNES DMDA local callback Jacobian", (*dmdasnes->jacobianlocalvec)(&info, Xloc, A, B, dmdasnes->jacobianlocalctx));
158:     else {
159:       PetscCall(DMDAVecGetArray(dm, Xloc, &x));
160:       PetscCallBack("SNES DMDA local callback Jacobian", (*dmdasnes->jacobianlocal)(&info, x, A, B, dmdasnes->jacobianlocalctx));
161:       PetscCall(DMDAVecRestoreArray(dm, Xloc, &x));
162:     }
163:     PetscCall(DMRestoreLocalVector(dm, &Xloc));
164:   } else {
165:     MatFDColoring fdcoloring;
166:     PetscCall(PetscObjectQuery((PetscObject)dm, "DMDASNES_FDCOLORING", (PetscObject *)&fdcoloring));
167:     if (!fdcoloring) {
168:       ISColoring coloring;

170:       PetscCall(DMCreateColoring(dm, dm->coloringtype, &coloring));
171:       PetscCall(MatFDColoringCreate(B, coloring, &fdcoloring));
172:       switch (dm->coloringtype) {
173:       case IS_COLORING_GLOBAL:
174:         PetscCall(MatFDColoringSetFunction(fdcoloring, (PetscErrorCode(*)(void))SNESComputeFunction_DMDA, dmdasnes));
175:         break;
176:       default:
177:         SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "No support for coloring type '%s'", ISColoringTypes[dm->coloringtype]);
178:       }
179:       PetscCall(PetscObjectSetOptionsPrefix((PetscObject)fdcoloring, ((PetscObject)dm)->prefix));
180:       PetscCall(MatFDColoringSetFromOptions(fdcoloring));
181:       PetscCall(MatFDColoringSetUp(B, coloring, fdcoloring));
182:       PetscCall(ISColoringDestroy(&coloring));
183:       PetscCall(PetscObjectCompose((PetscObject)dm, "DMDASNES_FDCOLORING", (PetscObject)fdcoloring));
184:       PetscCall(PetscObjectDereference((PetscObject)fdcoloring));

186:       /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call
187:        * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the
188:        * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually
189:        * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be
190:        * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed.
191:        */
192:       PetscCall(PetscObjectDereference((PetscObject)dm));
193:     }
194:     PetscCall(MatFDColoringApply(B, fdcoloring, X, snes));
195:   }
196:   /* This will be redundant if the user called both, but it's too common to forget. */
197:   if (A != B) {
198:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
199:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
200:   }
201:   PetscFunctionReturn(PETSC_SUCCESS);
202: }

204: /*@C
205:   DMDASNESSetFunctionLocal - set a local residual evaluation function for use with `DMDA`

207:   Logically Collective

209:   Input Parameters:
210: + dm    - `DM` to associate callback with
211: . imode - `INSERT_VALUES` if local function computes owned part, `ADD_VALUES` if it contributes to ghosted part
212: . func  - local residual evaluation
213: - ctx   - optional context for local residual evaluation

215:   Calling sequence of `func`:
216: $   PetscErrorCode func(DMDALocalInfo *info, void *x, void *f, void *ctx)
217: + info - `DMDALocalInfo` defining the subdomain to evaluate the residual on
218: . x    - dimensional pointer to state at which to evaluate residual (e.g. PetscScalar *x or **x or ***x)
219: . f    - dimensional pointer to residual, write the residual here (e.g. PetscScalar *f or **f or ***f)
220: - ctx  - optional context passed above

222:   Level: beginner

224: .seealso: `DMDA`, `DMDASNESSetJacobianLocal()`, `DMSNESSetFunction()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`
225: @*/
226: PetscErrorCode DMDASNESSetFunctionLocal(DM dm, InsertMode imode, PetscErrorCode (*func)(DMDALocalInfo *info, void *x, void *f, void *ctx), void *ctx)
227: {
228:   DMSNES     sdm;
229:   DMSNES_DA *dmdasnes;

231:   PetscFunctionBegin;
233:   PetscCall(DMGetDMSNESWrite(dm, &sdm));
234:   PetscCall(DMDASNESGetContext(dm, sdm, &dmdasnes));

236:   dmdasnes->residuallocalimode = imode;
237:   dmdasnes->residuallocal      = func;
238:   dmdasnes->residuallocalctx   = ctx;

240:   PetscCall(DMSNESSetFunction(dm, SNESComputeFunction_DMDA, dmdasnes));
241:   if (!sdm->ops->computejacobian) { /* Call us for the Jacobian too, can be overridden by the user. */
242:     PetscCall(DMSNESSetJacobian(dm, SNESComputeJacobian_DMDA, dmdasnes));
243:   }
244:   PetscFunctionReturn(PETSC_SUCCESS);
245: }

247: /*@C
248:   DMDASNESSetFunctionLocalVec - set a local residual evaluation function that operates on a local vector for `DMDA`

250:   Logically Collective

252:   Input Parameters:
253: + dm    - `DM` to associate callback with
254: . imode - `INSERT_VALUES` if local function computes owned part, `ADD_VALUES` if it contributes to ghosted part
255: . func  - local residual evaluation
256: - ctx   - optional context for local residual evaluation

258:   Calling sequence of `func`:
259: $   PetscErrorCode func(DMDALocalInfo *info, Vec x, Vec f, void *ctx),
260: + info - `DMDALocalInfo` defining the subdomain to evaluate the residual on
261: . x    - state vector at which to evaluate residual
262: . f    - residual vector
263: - ctx  - optional context passed above

265:   Level: beginner

267: .seealso: `DMDA`, `DMDASNESSetFunctionLocal()`, `DMDASNESSetJacobianLocalVec()`, `DMSNESSetFunction()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`
268: @*/
269: PetscErrorCode DMDASNESSetFunctionLocalVec(DM dm, InsertMode imode, PetscErrorCode (*func)(DMDALocalInfo *info, Vec x, Vec f, void *ctx), void *ctx)
270: {
271:   DMSNES     sdm;
272:   DMSNES_DA *dmdasnes;

274:   PetscFunctionBegin;
276:   PetscCall(DMGetDMSNESWrite(dm, &sdm));
277:   PetscCall(DMDASNESGetContext(dm, sdm, &dmdasnes));

279:   dmdasnes->residuallocalimode = imode;
280:   dmdasnes->residuallocalvec   = func;
281:   dmdasnes->residuallocalctx   = ctx;

283:   PetscCall(DMSNESSetFunction(dm, SNESComputeFunction_DMDA, dmdasnes));
284:   if (!sdm->ops->computejacobian) { /* Call us for the Jacobian too, can be overridden by the user. */
285:     PetscCall(DMSNESSetJacobian(dm, SNESComputeJacobian_DMDA, dmdasnes));
286:   }
287:   PetscFunctionReturn(PETSC_SUCCESS);
288: }

290: /*@C
291:   DMDASNESSetJacobianLocal - set a local Jacobian evaluation function for use with `DMDA`

293:   Logically Collective

295:   Input Parameters:
296: + dm   - `DM` to associate callback with
297: . func - local Jacobian evaluation
298: - ctx  - optional context for local Jacobian evaluation

300:   Calling sequence of `func`:
301: $  PetscErrorCode func(DMDALocalInfo *info, void *x, Mat J, Mat M, void *ctx),
302: + info - `DMDALocalInfo` defining the subdomain to evaluate the Jacobian at
303: . x    - dimensional pointer to state at which to evaluate Jacobian (e.g. PetscScalar *x or **x or ***x)
304: . J    - Mat object for the Jacobian
305: . M    - Mat object for the Jacobian preconditioner matrix, often `J`
306: - ctx  - optional context passed above

308:   Level: beginner

310: .seealso: `DMDA`, `DMDASNESSetFunctionLocal()`, `DMSNESSetJacobian()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`
311: @*/
312: PetscErrorCode DMDASNESSetJacobianLocal(DM dm, PetscErrorCode (*func)(DMDALocalInfo *info, void *x, Mat J, Mat M, void *ctx), void *ctx)
313: {
314:   DMSNES     sdm;
315:   DMSNES_DA *dmdasnes;

317:   PetscFunctionBegin;
319:   PetscCall(DMGetDMSNESWrite(dm, &sdm));
320:   PetscCall(DMDASNESGetContext(dm, sdm, &dmdasnes));

322:   dmdasnes->jacobianlocal    = func;
323:   dmdasnes->jacobianlocalctx = ctx;

325:   PetscCall(DMSNESSetJacobian(dm, SNESComputeJacobian_DMDA, dmdasnes));
326:   PetscFunctionReturn(PETSC_SUCCESS);
327: }

329: /*@C
330:   DMDASNESSetJacobianLocalVec - set a local Jacobian evaluation function that operates on a local vector with `DMDA`

332:   Logically Collective

334:   Input Parameters:
335: + dm   - `DM` to associate callback with
336: . func - local Jacobian evaluation
337: - ctx  - optional context for local Jacobian evaluation

339:   Calling sequence of `func`:
340: $  PetscErrorCode (*func)(DMDALocalInfo *info, Vec x, Mat J, Mat M, void *ctx),
341: + info - `DMDALocalInfo` defining the subdomain to evaluate the Jacobian at
342: . x    - state vector at which to evaluate Jacobian
343: . J    - the Jacobian
344: . M    - approximate Jacobian from which the preconditioner will be computed, often `J`
345: - ctx  - optional context passed above

347:   Level: beginner

349: .seealso: `DMDA`, `DMDASNESSetJacobianLocal()`, `DMDASNESSetFunctionLocalVec()`, `DMSNESSetJacobian()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`
350: @*/
351: PetscErrorCode DMDASNESSetJacobianLocalVec(DM dm, PetscErrorCode (*func)(DMDALocalInfo *info, Vec x, Mat J, Mat M, void *), void *ctx)
352: {
353:   DMSNES     sdm;
354:   DMSNES_DA *dmdasnes;

356:   PetscFunctionBegin;
358:   PetscCall(DMGetDMSNESWrite(dm, &sdm));
359:   PetscCall(DMDASNESGetContext(dm, sdm, &dmdasnes));

361:   dmdasnes->jacobianlocalvec = func;
362:   dmdasnes->jacobianlocalctx = ctx;

364:   PetscCall(DMSNESSetJacobian(dm, SNESComputeJacobian_DMDA, dmdasnes));
365:   PetscFunctionReturn(PETSC_SUCCESS);
366: }

368: // PetscClangLinter pragma disable: -fdoc-param-list-func-parameter-documentation
369: /*@C
370:   DMDASNESSetObjectiveLocal - set a local residual evaluation function to used with a `DMDA`

372:   Logically Collective

374:   Input Parameters:
375: + dm   - `DM` to associate callback with
376: . func - local objective evaluation
377: - ctx  - optional context for local residual evaluation

379:   Calling sequence of `func`:
380: $  PetscErrorCode func(DMDALocalInfo *info, void *x, PetscReal obj, void *ctx);
381: +  info - defines the subdomain to evaluate the residual on
382: .  x - dimensional pointer to state at which to evaluate residual
383: .  obj - eventual objective value
384: - ctx - optional context passed above

386:   Level: beginner

388: .seealso: `DMDA`, `DMSNESSetFunction()`, `DMDASNESSetJacobianLocal()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`
389: @*/
390: PetscErrorCode DMDASNESSetObjectiveLocal(DM dm, DMDASNESObjective func, void *ctx)
391: {
392:   DMSNES     sdm;
393:   DMSNES_DA *dmdasnes;

395:   PetscFunctionBegin;
397:   PetscCall(DMGetDMSNESWrite(dm, &sdm));
398:   PetscCall(DMDASNESGetContext(dm, sdm, &dmdasnes));

400:   dmdasnes->objectivelocal    = func;
401:   dmdasnes->objectivelocalctx = ctx;

403:   PetscCall(DMSNESSetObjective(dm, SNESComputeObjective_DMDA, dmdasnes));
404:   PetscFunctionReturn(PETSC_SUCCESS);
405: }

407: // PetscClangLinter pragma disable: -fdoc-param-list-func-parameter-documentation
408: /*@C
409:   DMDASNESSetObjectiveLocalVec - set a local residual evaluation function that operates on a local vector with `DMDA`

411:   Logically Collective

413:   Input Parameters:
414: + dm   - `DM` to associate callback with
415: . func - local objective evaluation
416: - ctx  - optional context for local residual evaluation

418:   Calling sequence of `func`:
419: $  PetscErrorCode func(DMDALocalInfo *info, Vec x, PetscReal *ob, void *ctx);
420: +  info - defines the subdomain to evaluate the residual on
421: .  x - state vector at which to evaluate residual
422: .  obj - eventual objective value
423: - ctx - optional context passed above

425:   Level: beginner

427: .seealso: `DMDA`, `DMDASNESSetObjectiveLocal()`, `DMSNESSetFunction()`, `DMDASNESSetJacobianLocalVec()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`
428: @*/
429: PetscErrorCode DMDASNESSetObjectiveLocalVec(DM dm, DMDASNESObjectiveVec func, void *ctx)
430: {
431:   DMSNES     sdm;
432:   DMSNES_DA *dmdasnes;

434:   PetscFunctionBegin;
436:   PetscCall(DMGetDMSNESWrite(dm, &sdm));
437:   PetscCall(DMDASNESGetContext(dm, sdm, &dmdasnes));

439:   dmdasnes->objectivelocalvec = func;
440:   dmdasnes->objectivelocalctx = ctx;

442:   PetscCall(DMSNESSetObjective(dm, SNESComputeObjective_DMDA, dmdasnes));
443:   PetscFunctionReturn(PETSC_SUCCESS);
444: }

446: static PetscErrorCode SNESComputePicard_DMDA(SNES snes, Vec X, Vec F, void *ctx)
447: {
448:   DM            dm;
449:   DMSNES_DA    *dmdasnes = (DMSNES_DA *)ctx;
450:   DMDALocalInfo info;
451:   Vec           Xloc;
452:   void         *x, *f;

454:   PetscFunctionBegin;
458:   PetscCheck(dmdasnes->rhsplocal, PetscObjectComm((PetscObject)snes), PETSC_ERR_PLIB, "Corrupt context");
459:   PetscCall(SNESGetDM(snes, &dm));
460:   PetscCall(DMGetLocalVector(dm, &Xloc));
461:   PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, Xloc));
462:   PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, Xloc));
463:   PetscCall(DMDAGetLocalInfo(dm, &info));
464:   PetscCall(DMDAVecGetArray(dm, Xloc, &x));
465:   switch (dmdasnes->residuallocalimode) {
466:   case INSERT_VALUES: {
467:     PetscCall(DMDAVecGetArray(dm, F, &f));
468:     PetscCallBack("SNES Picard DMDA local callback function", (*dmdasnes->rhsplocal)(&info, x, f, dmdasnes->picardlocalctx));
469:     PetscCall(DMDAVecRestoreArray(dm, F, &f));
470:   } break;
471:   case ADD_VALUES: {
472:     Vec Floc;
473:     PetscCall(DMGetLocalVector(dm, &Floc));
474:     PetscCall(VecZeroEntries(Floc));
475:     PetscCall(DMDAVecGetArray(dm, Floc, &f));
476:     PetscCallBack("SNES Picard DMDA local callback function", (*dmdasnes->rhsplocal)(&info, x, f, dmdasnes->picardlocalctx));
477:     PetscCall(DMDAVecRestoreArray(dm, Floc, &f));
478:     PetscCall(VecZeroEntries(F));
479:     PetscCall(DMLocalToGlobalBegin(dm, Floc, ADD_VALUES, F));
480:     PetscCall(DMLocalToGlobalEnd(dm, Floc, ADD_VALUES, F));
481:     PetscCall(DMRestoreLocalVector(dm, &Floc));
482:   } break;
483:   default:
484:     SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_INCOMP, "Cannot use imode=%d", (int)dmdasnes->residuallocalimode);
485:   }
486:   PetscCall(DMDAVecRestoreArray(dm, Xloc, &x));
487:   PetscCall(DMRestoreLocalVector(dm, &Xloc));
488:   PetscFunctionReturn(PETSC_SUCCESS);
489: }

491: static PetscErrorCode SNESComputePicardJacobian_DMDA(SNES snes, Vec X, Mat A, Mat B, void *ctx)
492: {
493:   DM            dm;
494:   DMSNES_DA    *dmdasnes = (DMSNES_DA *)ctx;
495:   DMDALocalInfo info;
496:   Vec           Xloc;
497:   void         *x;

499:   PetscFunctionBegin;
500:   PetscCheck(dmdasnes->jacobianplocal, PetscObjectComm((PetscObject)snes), PETSC_ERR_PLIB, "Corrupt context");
501:   PetscCall(SNESGetDM(snes, &dm));

503:   PetscCall(DMGetLocalVector(dm, &Xloc));
504:   PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, Xloc));
505:   PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, Xloc));
506:   PetscCall(DMDAGetLocalInfo(dm, &info));
507:   PetscCall(DMDAVecGetArray(dm, Xloc, &x));
508:   PetscCallBack("SNES Picard DMDA local callback Jacobian", (*dmdasnes->jacobianplocal)(&info, x, A, B, dmdasnes->picardlocalctx));
509:   PetscCall(DMDAVecRestoreArray(dm, Xloc, &x));
510:   PetscCall(DMRestoreLocalVector(dm, &Xloc));
511:   if (A != B) {
512:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
513:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
514:   }
515:   PetscFunctionReturn(PETSC_SUCCESS);
516: }

518: /*@C
519:   DMDASNESSetPicardLocal - set a local right hand side and matrix evaluation function for Picard iteration with `DMDA`

521:   Logically Collective

523:   Input Parameters:
524: + dm    - `DM` to associate callback with
525: . imode - `INSERT_VALUES` if local function computes owned part, `ADD_VALUES` if it contributes to ghosted part
526: . func  - local residual evaluation
527: . jac   - function to compute Jacobian
528: - ctx   - optional context for local residual evaluation

530:   Calling sequence of `func`:
531: $  PetscErrorCode func(DMDALocalInfo *info, void *x, void *f, void *ctx);
532: + info - defines the subdomain to evaluate the residual on
533: . x    - dimensional pointer to state at which to evaluate residual
534: . f    - dimensional pointer to residual, write the residual here
535: - ctx  - optional context passed above

537:   Calling sequence of `jac`:
538: $  PetscErrorCode jac(DMDALocalInfo *info, void *x, void *f, void *ctx);
539: + info - defines the subdomain to evaluate the residual on
540: . x    - dimensional pointer to state at which to evaluate residual
541: . jac  - the Jacobian
542: . Jp   - approximation to the Jacobian used to compute the preconditioner, often `J`
543: - ctx  - optional context passed above

545:   Level: beginner

547:   Note:
548:   The user must use `SNESSetFunction`(snes,`NULL`,`SNESPicardComputeFunction`,&user));
549:   in their code before calling this routine.

551: .seealso: `DMDA`, `DMSNESSetFunction()`, `DMDASNESSetJacobian()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`
552: @*/
553: PetscErrorCode DMDASNESSetPicardLocal(DM dm, InsertMode imode, PetscErrorCode (*func)(DMDALocalInfo *info, void *x, void *f, void *ctx), PetscErrorCode (*jac)(DMDALocalInfo *info, void *x, Mat jac, Mat Jp, void *ctx), void *ctx)
554: {
555:   DMSNES     sdm;
556:   DMSNES_DA *dmdasnes;

558:   PetscFunctionBegin;
560:   PetscCall(DMGetDMSNESWrite(dm, &sdm));
561:   PetscCall(DMDASNESGetContext(dm, sdm, &dmdasnes));

563:   dmdasnes->residuallocalimode = imode;
564:   dmdasnes->rhsplocal          = func;
565:   dmdasnes->jacobianplocal     = jac;
566:   dmdasnes->picardlocalctx     = ctx;

568:   PetscCall(DMSNESSetPicard(dm, SNESComputePicard_DMDA, SNESComputePicardJacobian_DMDA, dmdasnes));
569:   PetscCall(DMSNESSetMFFunction(dm, SNESComputeFunction_DMDA, dmdasnes));
570:   PetscFunctionReturn(PETSC_SUCCESS);
571: }