1 #include <petscdmda.h> /*I "petscdmda.h" I*/
2 #include <petsc/private/dmimpl.h>
3 #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/
4
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;
11
12 /* Vec version for vector data */
13 DMDASNESFunctionVecFn *residuallocalvec;
14 DMDASNESJacobianVecFn *jacobianlocalvec;
15 DMDASNESObjectiveVecFn *objectivelocalvec;
16
17 /* user contexts */
18 void *residuallocalctx;
19 void *jacobianlocalctx;
20 void *objectivelocalctx;
21 InsertMode residuallocalimode;
22
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;
28
DMSNESDestroy_DMDA(DMSNES sdm)29 static PetscErrorCode DMSNESDestroy_DMDA(DMSNES sdm)
30 {
31 PetscFunctionBegin;
32 PetscCall(PetscFree(sdm->data));
33 PetscFunctionReturn(PETSC_SUCCESS);
34 }
35
DMSNESDuplicate_DMDA(DMSNES oldsdm,DMSNES sdm)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 }
43
DMDASNESGetContext(DM dm,DMSNES sdm,DMSNES_DA ** dmdasnes)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 }
56
SNESComputeFunction_DMDA(SNES snes,Vec X,Vec F,PetscCtx ctx)57 static PetscErrorCode SNESComputeFunction_DMDA(SNES snes, Vec X, Vec F, PetscCtx ctx)
58 {
59 DM dm;
60 DMSNES_DA *dmdasnes = (DMSNES_DA *)ctx;
61 DMDALocalInfo info;
62 Vec Xloc;
63 void *x, *f, *rctx;
64
65 PetscFunctionBegin;
66 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
67 PetscValidHeaderSpecific(X, VEC_CLASSID, 2);
68 PetscValidHeaderSpecific(F, VEC_CLASSID, 3);
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->functiondomainerror));
113 PetscFunctionReturn(PETSC_SUCCESS);
114 }
115
SNESComputeObjective_DMDA(SNES snes,Vec X,PetscReal * ob,PetscCtx ctx)116 static PetscErrorCode SNESComputeObjective_DMDA(SNES snes, Vec X, PetscReal *ob, PetscCtx ctx)
117 {
118 DM dm;
119 DMSNES_DA *dmdasnes = (DMSNES_DA *)ctx;
120 DMDALocalInfo info;
121 Vec Xloc;
122 void *x, *octx;
123
124 PetscFunctionBegin;
125 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
126 PetscValidHeaderSpecific(X, VEC_CLASSID, 2);
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 }
145
146 /* Routine is called by example, hence must be labeled PETSC_EXTERN */
SNESComputeJacobian_DMDA(SNES snes,Vec X,Mat A,Mat B,PetscCtx ctx)147 PETSC_EXTERN PetscErrorCode SNESComputeJacobian_DMDA(SNES snes, Vec X, Mat A, Mat B, PetscCtx ctx)
148 {
149 DM dm;
150 DMSNES_DA *dmdasnes = (DMSNES_DA *)ctx;
151 DMDALocalInfo info;
152 Vec Xloc;
153 void *x, *jctx;
154
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;
176
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, (MatFDColoringFn *)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));
192
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 }
210
211 /*@C
212 DMDASNESSetFunctionLocal - set a local residual evaluation function for use with `DMDA`
213
214 Logically Collective
215
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
221
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
227
228 Level: beginner
229
230 .seealso: [](ch_snes), `DMDA`, `DMDASNESSetJacobianLocal()`, `DMSNESSetFunction()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`
231 @*/
DMDASNESSetFunctionLocal(DM dm,InsertMode imode,PetscErrorCode (* func)(DMDALocalInfo * info,void * x,void * f,PetscCtx ctx),PetscCtx ctx)232 PetscErrorCode DMDASNESSetFunctionLocal(DM dm, InsertMode imode, PetscErrorCode (*func)(DMDALocalInfo *info, void *x, void *f, PetscCtx ctx), PetscCtx ctx)
233 {
234 DMSNES sdm;
235 DMSNES_DA *dmdasnes;
236
237 PetscFunctionBegin;
238 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
239 PetscCall(DMGetDMSNESWrite(dm, &sdm));
240 PetscCall(DMDASNESGetContext(dm, sdm, &dmdasnes));
241
242 dmdasnes->residuallocalimode = imode;
243 dmdasnes->residuallocal = func;
244 dmdasnes->residuallocalctx = ctx;
245
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 }
252
253 /*@C
254 DMDASNESSetFunctionLocalVec - set a local residual evaluation function that operates on a local vector for `DMDA`
255
256 Logically Collective
257
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
263
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
269
270 Level: beginner
271
272 .seealso: [](ch_snes), `DMDA`, `DMDASNESSetFunctionLocal()`, `DMDASNESSetJacobianLocalVec()`, `DMSNESSetFunction()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`
273 @*/
DMDASNESSetFunctionLocalVec(DM dm,InsertMode imode,PetscErrorCode (* func)(DMDALocalInfo * info,Vec x,Vec f,PetscCtx ctx),PetscCtx ctx)274 PetscErrorCode DMDASNESSetFunctionLocalVec(DM dm, InsertMode imode, PetscErrorCode (*func)(DMDALocalInfo *info, Vec x, Vec f, PetscCtx ctx), PetscCtx ctx)
275 {
276 DMSNES sdm;
277 DMSNES_DA *dmdasnes;
278
279 PetscFunctionBegin;
280 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
281 PetscCall(DMGetDMSNESWrite(dm, &sdm));
282 PetscCall(DMDASNESGetContext(dm, sdm, &dmdasnes));
283
284 dmdasnes->residuallocalimode = imode;
285 dmdasnes->residuallocalvec = func;
286 dmdasnes->residuallocalctx = ctx;
287
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 }
294
295 /*@C
296 DMDASNESSetJacobianLocal - set a local Jacobian evaluation function for use with `DMDA`
297
298 Logically Collective
299
300 Input Parameters:
301 + dm - `DM` to associate callback with
302 . func - local Jacobian evaluation function
303 - ctx - optional context for local Jacobian evaluation
304
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 used to compute the preconditioner often `J`
310 - ctx - optional context passed above
311
312 Level: beginner
313
314 Note:
315 The `J` and `M` matrices are created internally by `DMCreateMatrix()`
316
317 .seealso: [](ch_snes), `DMDA`, `DMDASNESSetFunctionLocal()`, `DMSNESSetJacobian()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`
318 @*/
DMDASNESSetJacobianLocal(DM dm,PetscErrorCode (* func)(DMDALocalInfo * info,void * x,Mat J,Mat M,PetscCtx ctx),PetscCtx ctx)319 PetscErrorCode DMDASNESSetJacobianLocal(DM dm, PetscErrorCode (*func)(DMDALocalInfo *info, void *x, Mat J, Mat M, PetscCtx ctx), PetscCtx ctx)
320 {
321 DMSNES sdm;
322 DMSNES_DA *dmdasnes;
323
324 PetscFunctionBegin;
325 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
326 PetscCall(DMGetDMSNESWrite(dm, &sdm));
327 PetscCall(DMDASNESGetContext(dm, sdm, &dmdasnes));
328
329 dmdasnes->jacobianlocal = func;
330 dmdasnes->jacobianlocalctx = ctx;
331
332 PetscCall(DMSNESSetJacobian(dm, SNESComputeJacobian_DMDA, dmdasnes));
333 PetscFunctionReturn(PETSC_SUCCESS);
334 }
335
336 /*@C
337 DMDASNESSetJacobianLocalVec - set a local Jacobian evaluation function that operates on a local vector with `DMDA`
338
339 Logically Collective
340
341 Input Parameters:
342 + dm - `DM` to associate callback with
343 . func - local Jacobian evaluation
344 - ctx - optional context for local Jacobian evaluation
345
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
352
353 Level: beginner
354
355 .seealso: [](ch_snes), `DMDA`, `DMDASNESSetJacobianLocal()`, `DMDASNESSetFunctionLocalVec()`, `DMSNESSetJacobian()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`
356 @*/
DMDASNESSetJacobianLocalVec(DM dm,PetscErrorCode (* func)(DMDALocalInfo * info,Vec x,Mat J,Mat M,void *),PetscCtx ctx)357 PetscErrorCode DMDASNESSetJacobianLocalVec(DM dm, PetscErrorCode (*func)(DMDALocalInfo *info, Vec x, Mat J, Mat M, void *), PetscCtx ctx)
358 {
359 DMSNES sdm;
360 DMSNES_DA *dmdasnes;
361
362 PetscFunctionBegin;
363 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
364 PetscCall(DMGetDMSNESWrite(dm, &sdm));
365 PetscCall(DMDASNESGetContext(dm, sdm, &dmdasnes));
366
367 dmdasnes->jacobianlocalvec = func;
368 dmdasnes->jacobianlocalctx = ctx;
369
370 PetscCall(DMSNESSetJacobian(dm, SNESComputeJacobian_DMDA, dmdasnes));
371 PetscFunctionReturn(PETSC_SUCCESS);
372 }
373
374 /*@C
375 DMDASNESSetObjectiveLocal - set a local residual evaluation function to used with a `DMDA`
376
377 Logically Collective
378
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
383
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
389
390 Level: beginner
391
392 .seealso: [](ch_snes), `DMDA`, `DMSNESSetFunction()`, `DMDASNESSetJacobianLocal()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`, `DMDASNESObjectiveFn`
393 @*/
DMDASNESSetObjectiveLocal(DM dm,PetscErrorCode (* func)(DMDALocalInfo * info,void * x,PetscReal * obj,void *),PetscCtx ctx)394 PetscErrorCode DMDASNESSetObjectiveLocal(DM dm, PetscErrorCode (*func)(DMDALocalInfo *info, void *x, PetscReal *obj, void *), PetscCtx ctx)
395 {
396 DMSNES sdm;
397 DMSNES_DA *dmdasnes;
398
399 PetscFunctionBegin;
400 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
401 PetscCall(DMGetDMSNESWrite(dm, &sdm));
402 PetscCall(DMDASNESGetContext(dm, sdm, &dmdasnes));
403
404 dmdasnes->objectivelocal = func;
405 dmdasnes->objectivelocalctx = ctx;
406
407 PetscCall(DMSNESSetObjective(dm, SNESComputeObjective_DMDA, dmdasnes));
408 PetscFunctionReturn(PETSC_SUCCESS);
409 }
410
411 /*@C
412 DMDASNESSetObjectiveLocalVec - set a local residual evaluation function that operates on a local vector with `DMDA`
413
414 Logically Collective
415
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
420
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
426
427 Level: beginner
428
429 .seealso: [](ch_snes), `DMDA`, `DMDASNESSetObjectiveLocal()`, `DMSNESSetFunction()`, `DMDASNESSetJacobianLocalVec()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`, `DMDASNESObjectiveVecFn`
430 @*/
DMDASNESSetObjectiveLocalVec(DM dm,PetscErrorCode (* func)(DMDALocalInfo * info,Vec x,PetscReal * obj,void *),PetscCtx ctx)431 PetscErrorCode DMDASNESSetObjectiveLocalVec(DM dm, PetscErrorCode (*func)(DMDALocalInfo *info, Vec x, PetscReal *obj, void *), PetscCtx ctx)
432 {
433 DMSNES sdm;
434 DMSNES_DA *dmdasnes;
435
436 PetscFunctionBegin;
437 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
438 PetscCall(DMGetDMSNESWrite(dm, &sdm));
439 PetscCall(DMDASNESGetContext(dm, sdm, &dmdasnes));
440
441 dmdasnes->objectivelocalvec = func;
442 dmdasnes->objectivelocalctx = ctx;
443
444 PetscCall(DMSNESSetObjective(dm, SNESComputeObjective_DMDA, dmdasnes));
445 PetscFunctionReturn(PETSC_SUCCESS);
446 }
447
SNESComputePicard_DMDA(SNES snes,Vec X,Vec F,PetscCtx ctx)448 static PetscErrorCode SNESComputePicard_DMDA(SNES snes, Vec X, Vec F, PetscCtx ctx)
449 {
450 DM dm;
451 DMSNES_DA *dmdasnes = (DMSNES_DA *)ctx;
452 DMDALocalInfo info;
453 Vec Xloc;
454 void *x, *f;
455
456 PetscFunctionBegin;
457 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
458 PetscValidHeaderSpecific(X, VEC_CLASSID, 2);
459 PetscValidHeaderSpecific(F, VEC_CLASSID, 3);
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 }
492
SNESComputePicardJacobian_DMDA(SNES snes,Vec X,Mat A,Mat B,PetscCtx ctx)493 static PetscErrorCode SNESComputePicardJacobian_DMDA(SNES snes, Vec X, Mat A, Mat B, PetscCtx ctx)
494 {
495 DM dm;
496 DMSNES_DA *dmdasnes = (DMSNES_DA *)ctx;
497 DMDALocalInfo info;
498 Vec Xloc;
499 void *x;
500
501 PetscFunctionBegin;
502 PetscCheck(dmdasnes->jacobianplocal, PetscObjectComm((PetscObject)snes), PETSC_ERR_PLIB, "Corrupt context");
503 PetscCall(SNESGetDM(snes, &dm));
504
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 }
519
520 /*@C
521 DMDASNESSetPicardLocal - set a local right-hand side and matrix evaluation function for Picard iteration with `DMDA`
522
523 Logically Collective
524
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
531
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
537
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
544
545 Level: beginner
546
547 Note:
548 The user must use `SNESSetFunction`(`snes`,`NULL`,`SNESPicardComputeFunction`,&user));
549 in their code before calling this routine.
550
551 .seealso: [](ch_snes), `SNES`, `DMDA`, `DMSNESSetFunction()`, `DMDASNESSetJacobian()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`
552 @*/
DMDASNESSetPicardLocal(DM dm,InsertMode imode,PetscErrorCode (* func)(DMDALocalInfo * info,void * x,void * f,PetscCtx ctx),PetscErrorCode (* jac)(DMDALocalInfo * info,void * x,Mat jac,Mat Jp,PetscCtx ctx),PetscCtx ctx)553 PetscErrorCode DMDASNESSetPicardLocal(DM dm, InsertMode imode, PetscErrorCode (*func)(DMDALocalInfo *info, void *x, void *f, PetscCtx ctx), PetscErrorCode (*jac)(DMDALocalInfo *info, void *x, Mat jac, Mat Jp, PetscCtx ctx), PetscCtx ctx)
554 {
555 DMSNES sdm;
556 DMSNES_DA *dmdasnes;
557
558 PetscFunctionBegin;
559 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
560 PetscCall(DMGetDMSNESWrite(dm, &sdm));
561 PetscCall(DMDASNESGetContext(dm, sdm, &dmdasnes));
562
563 dmdasnes->residuallocalimode = imode;
564 dmdasnes->rhsplocal = func;
565 dmdasnes->jacobianplocal = jac;
566 dmdasnes->picardlocalctx = ctx;
567
568 PetscCall(DMSNESSetPicard(dm, SNESComputePicard_DMDA, SNESComputePicardJacobian_DMDA, dmdasnes));
569 PetscCall(DMSNESSetMFFunction(dm, SNESComputeFunction_DMDA, dmdasnes));
570 PetscFunctionReturn(PETSC_SUCCESS);
571 }
572