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: /* array versions for vector data */
8: DMDASNESFunctionFn *residuallocal;
9: DMDASNESJacobianFn *jacobianlocal;
10: DMDASNESObjectiveFn *objectivelocal;
12: /* Vec version for vector data */
13: DMDASNESFunctionVecFn *residuallocalvec;
14: DMDASNESJacobianVecFn *jacobianlocalvec;
15: DMDASNESObjectiveVecFn *objectivelocalvec;
17: /* user contexts */
18: void *residuallocalctx;
19: void *jacobianlocalctx;
20: void *objectivelocalctx;
21: InsertMode residuallocalimode;
23: /* For Picard iteration defined locally */
24: PetscErrorCode (*rhsplocal)(DMDALocalInfo *, void *, void *, void *);
25: PetscErrorCode (*jacobianplocal)(DMDALocalInfo *, void *, Mat, Mat, void *);
26: void *picardlocalctx;
27: } DMSNES_DA;
29: static PetscErrorCode DMSNESDestroy_DMDA(DMSNES sdm)
30: {
31: PetscFunctionBegin;
32: PetscCall(PetscFree(sdm->data));
33: PetscFunctionReturn(PETSC_SUCCESS);
34: }
36: static PetscErrorCode DMSNESDuplicate_DMDA(DMSNES oldsdm, DMSNES sdm)
37: {
38: PetscFunctionBegin;
39: PetscCall(PetscNew((DMSNES_DA **)&sdm->data));
40: if (oldsdm->data) PetscCall(PetscMemcpy(sdm->data, oldsdm->data, sizeof(DMSNES_DA)));
41: PetscFunctionReturn(PETSC_SUCCESS);
42: }
44: static PetscErrorCode DMDASNESGetContext(DM dm, DMSNES sdm, DMSNES_DA **dmdasnes)
45: {
46: PetscFunctionBegin;
47: *dmdasnes = NULL;
48: if (!sdm->data) {
49: PetscCall(PetscNew((DMSNES_DA **)&sdm->data));
50: sdm->ops->destroy = DMSNESDestroy_DMDA;
51: sdm->ops->duplicate = DMSNESDuplicate_DMDA;
52: }
53: *dmdasnes = (DMSNES_DA *)sdm->data;
54: PetscFunctionReturn(PETSC_SUCCESS);
55: }
57: static PetscErrorCode SNESComputeFunction_DMDA(SNES snes, Vec X, Vec F, void *ctx)
58: {
59: DM dm;
60: DMSNES_DA *dmdasnes = (DMSNES_DA *)ctx;
61: DMDALocalInfo info;
62: Vec Xloc;
63: void *x, *f, *rctx;
65: PetscFunctionBegin;
69: PetscCheck(dmdasnes->residuallocal || dmdasnes->residuallocalvec, PetscObjectComm((PetscObject)snes), PETSC_ERR_PLIB, "Corrupt context");
70: PetscCall(SNESGetDM(snes, &dm));
71: PetscCall(DMGetLocalVector(dm, &Xloc));
72: PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, Xloc));
73: PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, Xloc));
74: PetscCall(DMDAGetLocalInfo(dm, &info));
75: rctx = dmdasnes->residuallocalctx ? dmdasnes->residuallocalctx : snes->ctx;
76: switch (dmdasnes->residuallocalimode) {
77: case INSERT_VALUES: {
78: PetscCall(PetscLogEventBegin(SNES_FunctionEval, snes, X, F, 0));
79: if (dmdasnes->residuallocalvec) PetscCallBack("SNES DMDA local callback function", (*dmdasnes->residuallocalvec)(&info, Xloc, F, rctx));
80: else {
81: PetscCall(DMDAVecGetArray(dm, Xloc, &x));
82: PetscCall(DMDAVecGetArray(dm, F, &f));
83: PetscCallBack("SNES DMDA local callback function", (*dmdasnes->residuallocal)(&info, x, f, rctx));
84: PetscCall(DMDAVecRestoreArray(dm, Xloc, &x));
85: PetscCall(DMDAVecRestoreArray(dm, F, &f));
86: }
87: PetscCall(PetscLogEventEnd(SNES_FunctionEval, snes, X, F, 0));
88: } break;
89: case ADD_VALUES: {
90: Vec Floc;
91: PetscCall(DMGetLocalVector(dm, &Floc));
92: PetscCall(VecZeroEntries(Floc));
93: PetscCall(PetscLogEventBegin(SNES_FunctionEval, snes, X, F, 0));
94: if (dmdasnes->residuallocalvec) PetscCallBack("SNES DMDA local callback function", (*dmdasnes->residuallocalvec)(&info, Xloc, Floc, rctx));
95: else {
96: PetscCall(DMDAVecGetArray(dm, Xloc, &x));
97: PetscCall(DMDAVecGetArray(dm, Floc, &f));
98: PetscCallBack("SNES DMDA local callback function", (*dmdasnes->residuallocal)(&info, x, f, rctx));
99: PetscCall(DMDAVecRestoreArray(dm, Xloc, &x));
100: PetscCall(DMDAVecRestoreArray(dm, Floc, &f));
101: }
102: PetscCall(PetscLogEventEnd(SNES_FunctionEval, snes, X, F, 0));
103: PetscCall(VecZeroEntries(F));
104: PetscCall(DMLocalToGlobalBegin(dm, Floc, ADD_VALUES, F));
105: PetscCall(DMLocalToGlobalEnd(dm, Floc, ADD_VALUES, F));
106: PetscCall(DMRestoreLocalVector(dm, &Floc));
107: } break;
108: default:
109: SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_INCOMP, "Cannot use imode=%d", (int)dmdasnes->residuallocalimode);
110: }
111: PetscCall(DMRestoreLocalVector(dm, &Xloc));
112: PetscCall(VecFlag(F, snes->domainerror));
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: static PetscErrorCode SNESComputeObjective_DMDA(SNES snes, Vec X, PetscReal *ob, void *ctx)
117: {
118: DM dm;
119: DMSNES_DA *dmdasnes = (DMSNES_DA *)ctx;
120: DMDALocalInfo info;
121: Vec Xloc;
122: void *x, *octx;
124: PetscFunctionBegin;
127: PetscAssertPointer(ob, 3);
128: PetscCheck(dmdasnes->objectivelocal || dmdasnes->objectivelocalvec, PetscObjectComm((PetscObject)snes), PETSC_ERR_PLIB, "Corrupt context");
129: PetscCall(SNESGetDM(snes, &dm));
130: PetscCall(DMGetLocalVector(dm, &Xloc));
131: PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, Xloc));
132: PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, Xloc));
133: PetscCall(DMDAGetLocalInfo(dm, &info));
134: octx = dmdasnes->objectivelocalctx ? dmdasnes->objectivelocalctx : snes->ctx;
135: if (dmdasnes->objectivelocalvec) PetscCallBack("SNES DMDA local callback objective", (*dmdasnes->objectivelocalvec)(&info, Xloc, ob, octx));
136: else {
137: PetscCall(DMDAVecGetArray(dm, Xloc, &x));
138: PetscCallBack("SNES DMDA local callback objective", (*dmdasnes->objectivelocal)(&info, x, ob, octx));
139: PetscCall(DMDAVecRestoreArray(dm, Xloc, &x));
140: }
141: PetscCall(DMRestoreLocalVector(dm, &Xloc));
142: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, ob, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)snes)));
143: PetscFunctionReturn(PETSC_SUCCESS);
144: }
146: /* Routine is called by example, hence must be labeled PETSC_EXTERN */
147: PETSC_EXTERN PetscErrorCode SNESComputeJacobian_DMDA(SNES snes, Vec X, Mat A, Mat B, void *ctx)
148: {
149: DM dm;
150: DMSNES_DA *dmdasnes = (DMSNES_DA *)ctx;
151: DMDALocalInfo info;
152: Vec Xloc;
153: void *x, *jctx;
155: PetscFunctionBegin;
156: PetscCheck(dmdasnes->residuallocal || dmdasnes->residuallocalvec, PetscObjectComm((PetscObject)snes), PETSC_ERR_PLIB, "Corrupt context");
157: PetscCall(SNESGetDM(snes, &dm));
158: jctx = dmdasnes->jacobianlocalctx ? dmdasnes->jacobianlocalctx : snes->ctx;
159: if (dmdasnes->jacobianlocal || dmdasnes->jacobianlocalvec) {
160: PetscCall(DMGetLocalVector(dm, &Xloc));
161: PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, Xloc));
162: PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, Xloc));
163: PetscCall(DMDAGetLocalInfo(dm, &info));
164: if (dmdasnes->jacobianlocalvec) PetscCallBack("SNES DMDA local callback Jacobian", (*dmdasnes->jacobianlocalvec)(&info, Xloc, A, B, jctx));
165: else {
166: PetscCall(DMDAVecGetArray(dm, Xloc, &x));
167: PetscCallBack("SNES DMDA local callback Jacobian", (*dmdasnes->jacobianlocal)(&info, x, A, B, jctx));
168: PetscCall(DMDAVecRestoreArray(dm, Xloc, &x));
169: }
170: PetscCall(DMRestoreLocalVector(dm, &Xloc));
171: } else {
172: MatFDColoring fdcoloring;
173: PetscCall(PetscObjectQuery((PetscObject)dm, "DMDASNES_FDCOLORING", (PetscObject *)&fdcoloring));
174: if (!fdcoloring) {
175: ISColoring coloring;
177: PetscCall(DMCreateColoring(dm, dm->coloringtype, &coloring));
178: PetscCall(MatFDColoringCreate(B, coloring, &fdcoloring));
179: switch (dm->coloringtype) {
180: case IS_COLORING_GLOBAL:
181: PetscCall(MatFDColoringSetFunction(fdcoloring, (PetscErrorCode (*)(void))SNESComputeFunction_DMDA, dmdasnes));
182: break;
183: default:
184: SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "No support for coloring type '%s'", ISColoringTypes[dm->coloringtype]);
185: }
186: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)fdcoloring, ((PetscObject)dm)->prefix));
187: PetscCall(MatFDColoringSetFromOptions(fdcoloring));
188: PetscCall(MatFDColoringSetUp(B, coloring, fdcoloring));
189: PetscCall(ISColoringDestroy(&coloring));
190: PetscCall(PetscObjectCompose((PetscObject)dm, "DMDASNES_FDCOLORING", (PetscObject)fdcoloring));
191: PetscCall(PetscObjectDereference((PetscObject)fdcoloring));
193: /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call
194: * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the
195: * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually
196: * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be
197: * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed.
198: */
199: PetscCall(PetscObjectDereference((PetscObject)dm));
200: }
201: PetscCall(MatFDColoringApply(B, fdcoloring, X, snes));
202: }
203: /* This will be redundant if the user called both, but it's too common to forget. */
204: if (A != B) {
205: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
206: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
207: }
208: PetscFunctionReturn(PETSC_SUCCESS);
209: }
211: /*@C
212: DMDASNESSetFunctionLocal - set a local residual evaluation function for use with `DMDA`
214: Logically Collective
216: Input Parameters:
217: + dm - `DM` to associate callback with
218: . imode - `INSERT_VALUES` if local function computes owned part, `ADD_VALUES` if it contributes to ghosted part
219: . func - local residual evaluation
220: - ctx - optional context for local residual evaluation
222: Calling sequence of `func`:
223: + info - `DMDALocalInfo` defining the subdomain to evaluate the residual on
224: . x - dimensional pointer to state at which to evaluate residual (e.g. PetscScalar *x or **x or ***x)
225: . f - dimensional pointer to residual, write the residual here (e.g. PetscScalar *f or **f or ***f)
226: - ctx - optional context passed above
228: Level: beginner
230: .seealso: [](ch_snes), `DMDA`, `DMDASNESSetJacobianLocal()`, `DMSNESSetFunction()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`
231: @*/
232: PetscErrorCode DMDASNESSetFunctionLocal(DM dm, InsertMode imode, PetscErrorCode (*func)(DMDALocalInfo *info, void *x, void *f, void *ctx), void *ctx)
233: {
234: DMSNES sdm;
235: DMSNES_DA *dmdasnes;
237: PetscFunctionBegin;
239: PetscCall(DMGetDMSNESWrite(dm, &sdm));
240: PetscCall(DMDASNESGetContext(dm, sdm, &dmdasnes));
242: dmdasnes->residuallocalimode = imode;
243: dmdasnes->residuallocal = func;
244: dmdasnes->residuallocalctx = ctx;
246: PetscCall(DMSNESSetFunction(dm, SNESComputeFunction_DMDA, dmdasnes));
247: if (!sdm->ops->computejacobian) { /* Call us for the Jacobian too, can be overridden by the user. */
248: PetscCall(DMSNESSetJacobian(dm, SNESComputeJacobian_DMDA, dmdasnes));
249: }
250: PetscFunctionReturn(PETSC_SUCCESS);
251: }
253: /*@C
254: DMDASNESSetFunctionLocalVec - set a local residual evaluation function that operates on a local vector for `DMDA`
256: Logically Collective
258: Input Parameters:
259: + dm - `DM` to associate callback with
260: . imode - `INSERT_VALUES` if local function computes owned part, `ADD_VALUES` if it contributes to ghosted part
261: . func - local residual evaluation
262: - ctx - optional context for local residual evaluation
264: Calling sequence of `func`:
265: + info - `DMDALocalInfo` defining the subdomain to evaluate the residual on
266: . x - state vector at which to evaluate residual
267: . f - residual vector
268: - ctx - optional context passed above
270: Level: beginner
272: .seealso: [](ch_snes), `DMDA`, `DMDASNESSetFunctionLocal()`, `DMDASNESSetJacobianLocalVec()`, `DMSNESSetFunction()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`
273: @*/
274: PetscErrorCode DMDASNESSetFunctionLocalVec(DM dm, InsertMode imode, PetscErrorCode (*func)(DMDALocalInfo *info, Vec x, Vec f, void *ctx), void *ctx)
275: {
276: DMSNES sdm;
277: DMSNES_DA *dmdasnes;
279: PetscFunctionBegin;
281: PetscCall(DMGetDMSNESWrite(dm, &sdm));
282: PetscCall(DMDASNESGetContext(dm, sdm, &dmdasnes));
284: dmdasnes->residuallocalimode = imode;
285: dmdasnes->residuallocalvec = func;
286: dmdasnes->residuallocalctx = ctx;
288: PetscCall(DMSNESSetFunction(dm, SNESComputeFunction_DMDA, dmdasnes));
289: if (!sdm->ops->computejacobian) { /* Call us for the Jacobian too, can be overridden by the user. */
290: PetscCall(DMSNESSetJacobian(dm, SNESComputeJacobian_DMDA, dmdasnes));
291: }
292: PetscFunctionReturn(PETSC_SUCCESS);
293: }
295: /*@C
296: DMDASNESSetJacobianLocal - set a local Jacobian evaluation function for use with `DMDA`
298: Logically Collective
300: Input Parameters:
301: + dm - `DM` to associate callback with
302: . func - local Jacobian evaluation function
303: - ctx - optional context for local Jacobian evaluation
305: Calling sequence of `func`:
306: + info - `DMDALocalInfo` defining the subdomain to evaluate the Jacobian at
307: . x - dimensional pointer to state at which to evaluate Jacobian (e.g. PetscScalar *x or **x or ***x)
308: . J - `Mat` object for the Jacobian
309: . M - `Mat` object for the Jacobian preconditioner matrix, often `J`
310: - ctx - optional context passed above
312: Level: beginner
314: Note:
315: The `J` and `M` matrices are created internally by `DMCreateMatrix()`
317: .seealso: [](ch_snes), `DMDA`, `DMDASNESSetFunctionLocal()`, `DMSNESSetJacobian()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`
318: @*/
319: PetscErrorCode DMDASNESSetJacobianLocal(DM dm, PetscErrorCode (*func)(DMDALocalInfo *info, void *x, Mat J, Mat M, void *ctx), void *ctx)
320: {
321: DMSNES sdm;
322: DMSNES_DA *dmdasnes;
324: PetscFunctionBegin;
326: PetscCall(DMGetDMSNESWrite(dm, &sdm));
327: PetscCall(DMDASNESGetContext(dm, sdm, &dmdasnes));
329: dmdasnes->jacobianlocal = func;
330: dmdasnes->jacobianlocalctx = ctx;
332: PetscCall(DMSNESSetJacobian(dm, SNESComputeJacobian_DMDA, dmdasnes));
333: PetscFunctionReturn(PETSC_SUCCESS);
334: }
336: /*@C
337: DMDASNESSetJacobianLocalVec - set a local Jacobian evaluation function that operates on a local vector with `DMDA`
339: Logically Collective
341: Input Parameters:
342: + dm - `DM` to associate callback with
343: . func - local Jacobian evaluation
344: - ctx - optional context for local Jacobian evaluation
346: Calling sequence of `func`:
347: + info - `DMDALocalInfo` defining the subdomain to evaluate the Jacobian at
348: . x - state vector at which to evaluate Jacobian
349: . J - the Jacobian
350: . M - approximate Jacobian from which the preconditioner will be computed, often `J`
351: - ctx - optional context passed above
353: Level: beginner
355: .seealso: [](ch_snes), `DMDA`, `DMDASNESSetJacobianLocal()`, `DMDASNESSetFunctionLocalVec()`, `DMSNESSetJacobian()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`
356: @*/
357: PetscErrorCode DMDASNESSetJacobianLocalVec(DM dm, PetscErrorCode (*func)(DMDALocalInfo *info, Vec x, Mat J, Mat M, void *), void *ctx)
358: {
359: DMSNES sdm;
360: DMSNES_DA *dmdasnes;
362: PetscFunctionBegin;
364: PetscCall(DMGetDMSNESWrite(dm, &sdm));
365: PetscCall(DMDASNESGetContext(dm, sdm, &dmdasnes));
367: dmdasnes->jacobianlocalvec = func;
368: dmdasnes->jacobianlocalctx = ctx;
370: PetscCall(DMSNESSetJacobian(dm, SNESComputeJacobian_DMDA, dmdasnes));
371: PetscFunctionReturn(PETSC_SUCCESS);
372: }
374: /*@C
375: DMDASNESSetObjectiveLocal - set a local residual evaluation function to used with a `DMDA`
377: Logically Collective
379: Input Parameters:
380: + dm - `DM` to associate callback with
381: . func - local objective evaluation, see `DMDASNESSetObjectiveLocal` for the calling sequence
382: - ctx - optional context for local residual evaluation
384: Calling sequence of `func`:
385: + info - `DMDALocalInfo` defining the subdomain to evaluate the Jacobian at
386: . x - dimensional pointer to state at which to evaluate the objective (e.g. PetscScalar *x or **x or ***x)
387: . obj - returned objective value for the local subdomain
388: - ctx - optional context passed above
390: Level: beginner
392: .seealso: [](ch_snes), `DMDA`, `DMSNESSetFunction()`, `DMDASNESSetJacobianLocal()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`, `DMDASNESObjectiveFn`
393: @*/
394: PetscErrorCode DMDASNESSetObjectiveLocal(DM dm, PetscErrorCode (*func)(DMDALocalInfo *info, void *x, PetscReal *obj, void *), void *ctx)
395: {
396: DMSNES sdm;
397: DMSNES_DA *dmdasnes;
399: PetscFunctionBegin;
401: PetscCall(DMGetDMSNESWrite(dm, &sdm));
402: PetscCall(DMDASNESGetContext(dm, sdm, &dmdasnes));
404: dmdasnes->objectivelocal = func;
405: dmdasnes->objectivelocalctx = ctx;
407: PetscCall(DMSNESSetObjective(dm, SNESComputeObjective_DMDA, dmdasnes));
408: PetscFunctionReturn(PETSC_SUCCESS);
409: }
411: /*@C
412: DMDASNESSetObjectiveLocalVec - set a local residual evaluation function that operates on a local vector with `DMDA`
414: Logically Collective
416: Input Parameters:
417: + dm - `DM` to associate callback with
418: . func - local objective evaluation, see `DMDASNESSetObjectiveLocalVec` for the calling sequence
419: - ctx - optional context for local residual evaluation
421: Calling sequence of `func`:
422: + info - `DMDALocalInfo` defining the subdomain to evaluate the Jacobian at
423: . x - state vector at which to evaluate the objective
424: . obj - returned objective value for the local subdomain
425: - ctx - optional context passed above
427: Level: beginner
429: .seealso: [](ch_snes), `DMDA`, `DMDASNESSetObjectiveLocal()`, `DMSNESSetFunction()`, `DMDASNESSetJacobianLocalVec()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`, `DMDASNESObjectiveVecFn`
430: @*/
431: PetscErrorCode DMDASNESSetObjectiveLocalVec(DM dm, PetscErrorCode (*func)(DMDALocalInfo *info, Vec x, PetscReal *obj, void *), void *ctx)
432: {
433: DMSNES sdm;
434: DMSNES_DA *dmdasnes;
436: PetscFunctionBegin;
438: PetscCall(DMGetDMSNESWrite(dm, &sdm));
439: PetscCall(DMDASNESGetContext(dm, sdm, &dmdasnes));
441: dmdasnes->objectivelocalvec = func;
442: dmdasnes->objectivelocalctx = ctx;
444: PetscCall(DMSNESSetObjective(dm, SNESComputeObjective_DMDA, dmdasnes));
445: PetscFunctionReturn(PETSC_SUCCESS);
446: }
448: static PetscErrorCode SNESComputePicard_DMDA(SNES snes, Vec X, Vec F, void *ctx)
449: {
450: DM dm;
451: DMSNES_DA *dmdasnes = (DMSNES_DA *)ctx;
452: DMDALocalInfo info;
453: Vec Xloc;
454: void *x, *f;
456: PetscFunctionBegin;
460: PetscCheck(dmdasnes->rhsplocal, PetscObjectComm((PetscObject)snes), PETSC_ERR_PLIB, "Corrupt context");
461: PetscCall(SNESGetDM(snes, &dm));
462: PetscCall(DMGetLocalVector(dm, &Xloc));
463: PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, Xloc));
464: PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, Xloc));
465: PetscCall(DMDAGetLocalInfo(dm, &info));
466: PetscCall(DMDAVecGetArray(dm, Xloc, &x));
467: switch (dmdasnes->residuallocalimode) {
468: case INSERT_VALUES: {
469: PetscCall(DMDAVecGetArray(dm, F, &f));
470: PetscCallBack("SNES Picard DMDA local callback function", (*dmdasnes->rhsplocal)(&info, x, f, dmdasnes->picardlocalctx));
471: PetscCall(DMDAVecRestoreArray(dm, F, &f));
472: } break;
473: case ADD_VALUES: {
474: Vec Floc;
475: PetscCall(DMGetLocalVector(dm, &Floc));
476: PetscCall(VecZeroEntries(Floc));
477: PetscCall(DMDAVecGetArray(dm, Floc, &f));
478: PetscCallBack("SNES Picard DMDA local callback function", (*dmdasnes->rhsplocal)(&info, x, f, dmdasnes->picardlocalctx));
479: PetscCall(DMDAVecRestoreArray(dm, Floc, &f));
480: PetscCall(VecZeroEntries(F));
481: PetscCall(DMLocalToGlobalBegin(dm, Floc, ADD_VALUES, F));
482: PetscCall(DMLocalToGlobalEnd(dm, Floc, ADD_VALUES, F));
483: PetscCall(DMRestoreLocalVector(dm, &Floc));
484: } break;
485: default:
486: SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_INCOMP, "Cannot use imode=%d", (int)dmdasnes->residuallocalimode);
487: }
488: PetscCall(DMDAVecRestoreArray(dm, Xloc, &x));
489: PetscCall(DMRestoreLocalVector(dm, &Xloc));
490: PetscFunctionReturn(PETSC_SUCCESS);
491: }
493: static PetscErrorCode SNESComputePicardJacobian_DMDA(SNES snes, Vec X, Mat A, Mat B, void *ctx)
494: {
495: DM dm;
496: DMSNES_DA *dmdasnes = (DMSNES_DA *)ctx;
497: DMDALocalInfo info;
498: Vec Xloc;
499: void *x;
501: PetscFunctionBegin;
502: PetscCheck(dmdasnes->jacobianplocal, PetscObjectComm((PetscObject)snes), PETSC_ERR_PLIB, "Corrupt context");
503: PetscCall(SNESGetDM(snes, &dm));
505: PetscCall(DMGetLocalVector(dm, &Xloc));
506: PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, Xloc));
507: PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, Xloc));
508: PetscCall(DMDAGetLocalInfo(dm, &info));
509: PetscCall(DMDAVecGetArray(dm, Xloc, &x));
510: PetscCallBack("SNES Picard DMDA local callback Jacobian", (*dmdasnes->jacobianplocal)(&info, x, A, B, dmdasnes->picardlocalctx));
511: PetscCall(DMDAVecRestoreArray(dm, Xloc, &x));
512: PetscCall(DMRestoreLocalVector(dm, &Xloc));
513: if (A != B) {
514: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
515: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
516: }
517: PetscFunctionReturn(PETSC_SUCCESS);
518: }
520: /*@C
521: DMDASNESSetPicardLocal - set a local right-hand side and matrix evaluation function for Picard iteration with `DMDA`
523: Logically Collective
525: Input Parameters:
526: + dm - `DM` to associate callback with
527: . imode - `INSERT_VALUES` if local function computes owned part, `ADD_VALUES` if it contributes to ghosted part
528: . func - local residual evaluation
529: . jac - function to compute Jacobian
530: - ctx - optional context for local residual evaluation
532: Calling sequence of `func`:
533: + info - defines the subdomain to evaluate the residual on
534: . x - dimensional pointer to state at which to evaluate residual
535: . f - dimensional pointer to residual, write the residual here
536: - ctx - optional context passed above
538: Calling sequence of `jac`:
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: [](ch_snes), `SNES`, `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: }