xref: /petsc/src/snes/utils/dmdasnes.c (revision 785e854f82a3c614b452fca2cf5ad4f2afe8bdde)
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   PetscErrorCode (*residuallocal)(DMDALocalInfo*,void*,void*,void*);
8   PetscErrorCode (*jacobianlocal)(DMDALocalInfo*,void*,Mat,Mat,MatStructure*,void*);
9   PetscErrorCode (*objectivelocal)(DMDALocalInfo*,void*,PetscReal*,void*);
10   void       *residuallocalctx;
11   void       *jacobianlocalctx;
12   void       *objectivelocalctx;
13   InsertMode residuallocalimode;
14 
15   /*   For Picard iteration defined locally */
16   PetscErrorCode (*rhsplocal)(DMDALocalInfo*,void*,void*,void*);
17   PetscErrorCode (*jacobianplocal)(DMDALocalInfo*,void*,Mat,Mat,MatStructure*,void*);
18   void *picardlocalctx;
19 } DMSNES_DA;
20 
21 #undef __FUNCT__
22 #define __FUNCT__ "DMSNESDestroy_DMDA"
23 static PetscErrorCode DMSNESDestroy_DMDA(DMSNES sdm)
24 {
25   PetscErrorCode ierr;
26 
27   PetscFunctionBegin;
28   ierr = PetscFree(sdm->data);CHKERRQ(ierr);
29   PetscFunctionReturn(0);
30 }
31 
32 #undef __FUNCT__
33 #define __FUNCT__ "DMSNESDuplicate_DMDA"
34 static PetscErrorCode DMSNESDuplicate_DMDA(DMSNES oldsdm,DMSNES sdm)
35 {
36   PetscErrorCode ierr;
37 
38   PetscFunctionBegin;
39   ierr = PetscNewLog(sdm,DMSNES_DA,&sdm->data);CHKERRQ(ierr);
40   if (oldsdm->data) {
41     ierr = PetscMemcpy(sdm->data,oldsdm->data,sizeof(DMSNES_DA));CHKERRQ(ierr);
42   }
43   PetscFunctionReturn(0);
44 }
45 
46 
47 #undef __FUNCT__
48 #define __FUNCT__ "DMDASNESGetContext"
49 static PetscErrorCode DMDASNESGetContext(DM dm,DMSNES sdm,DMSNES_DA  **dmdasnes)
50 {
51   PetscErrorCode ierr;
52 
53   PetscFunctionBegin;
54   *dmdasnes = NULL;
55   if (!sdm->data) {
56     ierr                = PetscNewLog(dm,DMSNES_DA,&sdm->data);CHKERRQ(ierr);
57     sdm->ops->destroy   = DMSNESDestroy_DMDA;
58     sdm->ops->duplicate = DMSNESDuplicate_DMDA;
59   }
60   *dmdasnes = (DMSNES_DA*)sdm->data;
61   PetscFunctionReturn(0);
62 }
63 
64 #undef __FUNCT__
65 #define __FUNCT__ "SNESComputeFunction_DMDA"
66 /*
67   This function should eventually replace:
68     DMDAComputeFunction() and DMDAComputeFunction1()
69  */
70 static PetscErrorCode SNESComputeFunction_DMDA(SNES snes,Vec X,Vec F,void *ctx)
71 {
72   PetscErrorCode ierr;
73   DM             dm;
74   DMSNES_DA      *dmdasnes = (DMSNES_DA*)ctx;
75   DMDALocalInfo  info;
76   Vec            Xloc;
77   void           *x,*f;
78 
79   PetscFunctionBegin;
80   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
81   PetscValidHeaderSpecific(X,VEC_CLASSID,2);
82   PetscValidHeaderSpecific(F,VEC_CLASSID,3);
83   if (!dmdasnes->residuallocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
84   ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);
85   ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr);
86   ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
87   ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
88   ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);
89   ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr);
90   switch (dmdasnes->residuallocalimode) {
91   case INSERT_VALUES: {
92     ierr = DMDAVecGetArray(dm,F,&f);CHKERRQ(ierr);
93     CHKMEMQ;
94     ierr = (*dmdasnes->residuallocal)(&info,x,f,dmdasnes->residuallocalctx);CHKERRQ(ierr);
95     CHKMEMQ;
96     ierr = DMDAVecRestoreArray(dm,F,&f);CHKERRQ(ierr);
97   } break;
98   case ADD_VALUES: {
99     Vec Floc;
100     ierr = DMGetLocalVector(dm,&Floc);CHKERRQ(ierr);
101     ierr = VecZeroEntries(Floc);CHKERRQ(ierr);
102     ierr = DMDAVecGetArray(dm,Floc,&f);CHKERRQ(ierr);
103     CHKMEMQ;
104     ierr = (*dmdasnes->residuallocal)(&info,x,f,dmdasnes->residuallocalctx);CHKERRQ(ierr);
105     CHKMEMQ;
106     ierr = DMDAVecRestoreArray(dm,Floc,&f);CHKERRQ(ierr);
107     ierr = VecZeroEntries(F);CHKERRQ(ierr);
108     ierr = DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr);
109     ierr = DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr);
110     ierr = DMRestoreLocalVector(dm,&Floc);CHKERRQ(ierr);
111   } break;
112   default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdasnes->residuallocalimode);
113   }
114   ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr);
115   ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr);
116   PetscFunctionReturn(0);
117 }
118 
119 #undef __FUNCT__
120 #define __FUNCT__ "SNESComputeObjective_DMDA"
121 /*
122   This function should eventually replace:
123     DMDAComputeFunction() and DMDAComputeFunction1()
124  */
125 static PetscErrorCode SNESComputeObjective_DMDA(SNES snes,Vec X,PetscReal *ob,void *ctx)
126 {
127   PetscErrorCode ierr;
128   DM             dm;
129   DMSNES_DA      *dmdasnes = (DMSNES_DA*)ctx;
130   DMDALocalInfo  info;
131   Vec            Xloc;
132   void           *x;
133 
134   PetscFunctionBegin;
135   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
136   PetscValidHeaderSpecific(X,VEC_CLASSID,2);
137   PetscValidPointer(ob,3);
138   if (!dmdasnes->objectivelocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
139   ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);
140   ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr);
141   ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
142   ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
143   ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);
144   ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr);
145   CHKMEMQ;
146   ierr = (*dmdasnes->objectivelocal)(&info,x,ob,dmdasnes->objectivelocalctx);CHKERRQ(ierr);
147   CHKMEMQ;
148   ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr);
149   ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr);
150   PetscFunctionReturn(0);
151 }
152 
153 
154 #undef __FUNCT__
155 #define __FUNCT__ "SNESComputeJacobian_DMDA"
156 static PetscErrorCode SNESComputeJacobian_DMDA(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *mstr,void *ctx)
157 {
158   PetscErrorCode ierr;
159   DM             dm;
160   DMSNES_DA      *dmdasnes = (DMSNES_DA*)ctx;
161   DMDALocalInfo  info;
162   Vec            Xloc;
163   void           *x;
164 
165   PetscFunctionBegin;
166   if (!dmdasnes->residuallocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
167   ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);
168 
169   if (dmdasnes->jacobianlocal) {
170     ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr);
171     ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
172     ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
173     ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);
174     ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr);
175     CHKMEMQ;
176     ierr = (*dmdasnes->jacobianlocal)(&info,x,*A,*B,mstr,dmdasnes->jacobianlocalctx);CHKERRQ(ierr);
177     CHKMEMQ;
178     ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr);
179     ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr);
180   } else {
181     MatFDColoring fdcoloring;
182     ierr = PetscObjectQuery((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject*)&fdcoloring);CHKERRQ(ierr);
183     if (!fdcoloring) {
184       ISColoring coloring;
185 
186       ierr = DMCreateColoring(dm,dm->coloringtype,&coloring);CHKERRQ(ierr);
187       ierr = MatFDColoringCreate(*B,coloring,&fdcoloring);CHKERRQ(ierr);
188       switch (dm->coloringtype) {
189       case IS_COLORING_GLOBAL:
190         ierr = MatFDColoringSetFunction(fdcoloring,(PetscErrorCode (*)(void))SNESComputeFunction_DMDA,dmdasnes);CHKERRQ(ierr);
191         break;
192       default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"No support for coloring type '%s'",ISColoringTypes[dm->coloringtype]);
193       }
194       ierr = PetscObjectSetOptionsPrefix((PetscObject)fdcoloring,((PetscObject)dm)->prefix);CHKERRQ(ierr);
195       ierr = MatFDColoringSetFromOptions(fdcoloring);CHKERRQ(ierr);
196       ierr = MatFDColoringSetUp(*B,coloring,fdcoloring);CHKERRQ(ierr);
197       ierr = ISColoringDestroy(&coloring);CHKERRQ(ierr);
198       ierr = PetscObjectCompose((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject)fdcoloring);CHKERRQ(ierr);
199       ierr = PetscObjectDereference((PetscObject)fdcoloring);CHKERRQ(ierr);
200 
201       /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call
202        * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the
203        * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually
204        * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be
205        * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed.
206        */
207       ierr = PetscObjectDereference((PetscObject)dm);CHKERRQ(ierr);
208     }
209     *mstr = SAME_NONZERO_PATTERN;
210     ierr  = MatFDColoringApply(*B,fdcoloring,X,mstr,snes);CHKERRQ(ierr);
211   }
212   /* This will be redundant if the user called both, but it's too common to forget. */
213   if (*A != *B) {
214     ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
215     ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
216   }
217   PetscFunctionReturn(0);
218 }
219 
220 #undef __FUNCT__
221 #define __FUNCT__ "DMDASNESSetFunctionLocal"
222 /*@C
223    DMDASNESSetFunctionLocal - set a local residual evaluation function
224 
225    Logically Collective
226 
227    Input Arguments:
228 +  dm - DM to associate callback with
229 .  imode - INSERT_VALUES if local function computes owned part, ADD_VALUES if it contributes to ghosted part
230 .  func - local residual evaluation
231 -  ctx - optional context for local residual evaluation
232 
233    Calling sequence for func:
234 +  info - DMDALocalInfo defining the subdomain to evaluate the residual on
235 .  x - dimensional pointer to state at which to evaluate residual
236 .  f - dimensional pointer to residual, write the residual here
237 -  ctx - optional context passed above
238 
239    Level: beginner
240 
241 .seealso: DMSNESSetFunction(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
242 @*/
243 PetscErrorCode DMDASNESSetFunctionLocal(DM dm,InsertMode imode,PetscErrorCode (*func)(DMDALocalInfo*,void*,void*,void*),void *ctx)
244 {
245   PetscErrorCode ierr;
246   DMSNES         sdm;
247   DMSNES_DA      *dmdasnes;
248 
249   PetscFunctionBegin;
250   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
251   ierr = DMGetDMSNESWrite(dm,&sdm);CHKERRQ(ierr);
252   ierr = DMDASNESGetContext(dm,sdm,&dmdasnes);CHKERRQ(ierr);
253 
254   dmdasnes->residuallocalimode = imode;
255   dmdasnes->residuallocal      = func;
256   dmdasnes->residuallocalctx   = ctx;
257 
258   ierr = DMSNESSetFunction(dm,SNESComputeFunction_DMDA,dmdasnes);CHKERRQ(ierr);
259   if (!sdm->ops->computejacobian) {  /* Call us for the Jacobian too, can be overridden by the user. */
260     ierr = DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes);CHKERRQ(ierr);
261   }
262   PetscFunctionReturn(0);
263 }
264 
265 #undef __FUNCT__
266 #define __FUNCT__ "DMDASNESSetJacobianLocal"
267 /*@C
268    DMDASNESSetJacobianLocal - set a local Jacobian evaluation function
269 
270    Logically Collective
271 
272    Input Arguments:
273 +  dm - DM to associate callback with
274 .  func - local residual evaluation
275 -  ctx - optional context for local residual evaluation
276 
277    Calling sequence for func:
278 +  info - DMDALocalInfo defining the subdomain to evaluate the residual on
279 .  x - dimensional pointer to state at which to evaluate residual
280 .  f - dimensional pointer to residual, write the residual here
281 -  ctx - optional context passed above
282 
283    Level: beginner
284 
285 .seealso: DMSNESSetJacobian(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
286 @*/
287 PetscErrorCode DMDASNESSetJacobianLocal(DM dm,PetscErrorCode (*func)(DMDALocalInfo*,void*,Mat,Mat,MatStructure*,void*),void *ctx)
288 {
289   PetscErrorCode ierr;
290   DMSNES         sdm;
291   DMSNES_DA      *dmdasnes;
292 
293   PetscFunctionBegin;
294   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
295   ierr = DMGetDMSNESWrite(dm,&sdm);CHKERRQ(ierr);
296   ierr = DMDASNESGetContext(dm,sdm,&dmdasnes);CHKERRQ(ierr);
297 
298   dmdasnes->jacobianlocal    = func;
299   dmdasnes->jacobianlocalctx = ctx;
300 
301   ierr = DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes);CHKERRQ(ierr);
302   PetscFunctionReturn(0);
303 }
304 
305 
306 #undef __FUNCT__
307 #define __FUNCT__ "DMDASNESSetObjectiveLocal"
308 /*@C
309    DMDASNESSetObjectiveLocal - set a local residual evaluation function
310 
311    Logically Collective
312 
313    Input Arguments:
314 +  dm - DM to associate callback with
315 .  func - local objective evaluation
316 -  ctx - optional context for local residual evaluation
317 
318    Calling sequence for func:
319 +  info - DMDALocalInfo defining the subdomain to evaluate the residual on
320 .  x - dimensional pointer to state at which to evaluate residual
321 .  ob - eventual objective value
322 -  ctx - optional context passed above
323 
324    Level: beginner
325 
326 .seealso: DMSNESSetFunction(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
327 @*/
328 PetscErrorCode DMDASNESSetObjectiveLocal(DM dm,DMDASNESObjective func,void *ctx)
329 {
330   PetscErrorCode ierr;
331   DMSNES         sdm;
332   DMSNES_DA      *dmdasnes;
333 
334   PetscFunctionBegin;
335   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
336   ierr = DMGetDMSNESWrite(dm,&sdm);CHKERRQ(ierr);
337   ierr = DMDASNESGetContext(dm,sdm,&dmdasnes);CHKERRQ(ierr);
338 
339   dmdasnes->objectivelocal    = func;
340   dmdasnes->objectivelocalctx = ctx;
341 
342   ierr = DMSNESSetObjective(dm,SNESComputeObjective_DMDA,dmdasnes);CHKERRQ(ierr);
343   PetscFunctionReturn(0);
344 }
345 
346 #undef __FUNCT__
347 #define __FUNCT__ "SNESComputePicard_DMDA"
348 static PetscErrorCode SNESComputePicard_DMDA(SNES snes,Vec X,Vec F,void *ctx)
349 {
350   PetscErrorCode ierr;
351   DM             dm;
352   DMSNES_DA      *dmdasnes = (DMSNES_DA*)ctx;
353   DMDALocalInfo  info;
354   Vec            Xloc;
355   void           *x,*f;
356 
357   PetscFunctionBegin;
358   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
359   PetscValidHeaderSpecific(X,VEC_CLASSID,2);
360   PetscValidHeaderSpecific(F,VEC_CLASSID,3);
361   if (!dmdasnes->rhsplocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
362   ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);
363   ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr);
364   ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
365   ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
366   ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);
367   ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr);
368   switch (dmdasnes->residuallocalimode) {
369   case INSERT_VALUES: {
370     ierr = DMDAVecGetArray(dm,F,&f);CHKERRQ(ierr);
371     CHKMEMQ;
372     ierr = (*dmdasnes->rhsplocal)(&info,x,f,dmdasnes->picardlocalctx);CHKERRQ(ierr);
373     CHKMEMQ;
374     ierr = DMDAVecRestoreArray(dm,F,&f);CHKERRQ(ierr);
375   } break;
376   case ADD_VALUES: {
377     Vec Floc;
378     ierr = DMGetLocalVector(dm,&Floc);CHKERRQ(ierr);
379     ierr = VecZeroEntries(Floc);CHKERRQ(ierr);
380     ierr = DMDAVecGetArray(dm,Floc,&f);CHKERRQ(ierr);
381     CHKMEMQ;
382     ierr = (*dmdasnes->rhsplocal)(&info,x,f,dmdasnes->picardlocalctx);CHKERRQ(ierr);
383     CHKMEMQ;
384     ierr = DMDAVecRestoreArray(dm,Floc,&f);CHKERRQ(ierr);
385     ierr = VecZeroEntries(F);CHKERRQ(ierr);
386     ierr = DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr);
387     ierr = DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr);
388     ierr = DMRestoreLocalVector(dm,&Floc);CHKERRQ(ierr);
389   } break;
390   default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdasnes->residuallocalimode);
391   }
392   ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr);
393   ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr);
394   PetscFunctionReturn(0);
395 }
396 
397 #undef __FUNCT__
398 #define __FUNCT__ "SNESComputePicardJacobian_DMDA"
399 static PetscErrorCode SNESComputePicardJacobian_DMDA(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *mstr,void *ctx)
400 {
401   PetscErrorCode ierr;
402   DM             dm;
403   DMSNES_DA      *dmdasnes = (DMSNES_DA*)ctx;
404   DMDALocalInfo  info;
405   Vec            Xloc;
406   void           *x;
407 
408   PetscFunctionBegin;
409   if (!dmdasnes->jacobianplocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
410   ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);
411 
412   ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr);
413   ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
414   ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
415   ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);
416   ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr);
417   CHKMEMQ;
418   ierr = (*dmdasnes->jacobianplocal)(&info,x,*A,*B,mstr,dmdasnes->picardlocalctx);CHKERRQ(ierr);
419   CHKMEMQ;
420   ierr  = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr);
421   ierr  = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr);
422   *mstr = SAME_NONZERO_PATTERN;
423   if (*A != *B) {
424     ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
425     ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
426   }
427   PetscFunctionReturn(0);
428 }
429 
430 #undef __FUNCT__
431 #define __FUNCT__ "DMDASNESSetPicardLocal"
432 /*@C
433    DMDASNESSetPicardLocal - set a local right hand side and matrix evaluation function for Picard iteration
434 
435    Logically Collective
436 
437    Input Arguments:
438 +  dm - DM to associate callback with
439 .  imode - INSERT_VALUES if local function computes owned part, ADD_VALUES if it contributes to ghosted part
440 .  func - local residual evaluation
441 -  ctx - optional context for local residual evaluation
442 
443    Calling sequence for func:
444 +  info - DMDALocalInfo defining the subdomain to evaluate the residual on
445 .  x - dimensional pointer to state at which to evaluate residual
446 .  f - dimensional pointer to residual, write the residual here
447 -  ctx - optional context passed above
448 
449    Notes:  The user must use
450     extern PetscErrorCode  SNESPicardComputeFunction(SNES,Vec,Vec,void*);
451     extern PetscErrorCode  SNESPicardComputeJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
452     ierr = SNESSetFunction(snes,NULL,SNESPicardComputeFunction,&user);CHKERRQ(ierr);
453     in their code before calling this routine.
454 
455 
456    Level: beginner
457 
458 .seealso: DMSNESSetFunction(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
459 @*/
460 PetscErrorCode DMDASNESSetPicardLocal(DM dm,InsertMode imode,PetscErrorCode (*func)(DMDALocalInfo*,void*,void*,void*),
461                                       PetscErrorCode (*jac)(DMDALocalInfo*,void*,Mat,Mat,MatStructure*,void*),void *ctx)
462 {
463   PetscErrorCode ierr;
464   DMSNES         sdm;
465   DMSNES_DA      *dmdasnes;
466 
467   PetscFunctionBegin;
468   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
469   ierr = DMGetDMSNESWrite(dm,&sdm);CHKERRQ(ierr);
470   ierr = DMDASNESGetContext(dm,sdm,&dmdasnes);CHKERRQ(ierr);
471 
472   dmdasnes->residuallocalimode = imode;
473   dmdasnes->rhsplocal          = func;
474   dmdasnes->jacobianplocal     = jac;
475   dmdasnes->picardlocalctx     = ctx;
476 
477   ierr = DMSNESSetPicard(dm,SNESComputePicard_DMDA,SNESComputePicardJacobian_DMDA,dmdasnes);CHKERRQ(ierr);
478   PetscFunctionReturn(0);
479 }
480