xref: /petsc/src/snes/utils/dmdasnes.c (revision 111ade9efc8fd76368da9db8ddd81fbb99846c2d) 
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 = PETSC_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(((PetscObject)snes)->comm,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(((PetscObject)snes)->comm,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(((PetscObject)snes)->comm,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 #undef __FUNCT__
154 #define __FUNCT__ "SNESComputeLocalBlockFunction_DMDA"
155 /*
156   This function should eventually replace:
157     DMDAComputeFunction() and DMDAComputeFunction1()
158  */
159 PetscErrorCode SNESComputeLocalBlockFunction_DMDA(SNES snes,Vec X,Vec F,void *dmctx)
160 {
161   PetscErrorCode ierr;
162   DM             dm = (DM)dmctx; /* the global DMDA context */
163   DMDALocalInfo  info;
164   void           *x,*f;
165   DMSNES_DA      *dmdasnes;
166   DMSNES         sdm;
167 
168   PetscFunctionBegin;
169 
170   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
171   PetscValidHeaderSpecific(X,VEC_CLASSID,2);
172   PetscValidHeaderSpecific(F,VEC_CLASSID,3);
173   PetscValidHeaderSpecific(dm,DM_CLASSID,4);
174 
175   ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr);
176   ierr = DMDASNESGetContext(dm,sdm,&dmdasnes);CHKERRQ(ierr);
177   if (!dmdasnes->residuallocal) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_PLIB,"Corrupt context");
178   ierr = DMDAVecGetArray(dm,X,&x);CHKERRQ(ierr);
179   ierr = DMDAGetLocalBlockInfo(dm,&info);CHKERRQ(ierr);
180   switch (dmdasnes->residuallocalimode) {
181   case INSERT_VALUES: {
182     ierr = DMDAVecGetArray(dm,F,&f);CHKERRQ(ierr);
183     CHKMEMQ;
184     ierr = (*dmdasnes->residuallocal)(&info,x,f,dmdasnes->residuallocalctx);CHKERRQ(ierr);
185     CHKMEMQ;
186     ierr = DMDAVecRestoreArray(dm,F,&f);CHKERRQ(ierr);
187   } break;
188   case ADD_VALUES: {
189     Vec Floc;
190     ierr = DMGetLocalVector(dm,&Floc);CHKERRQ(ierr);
191     ierr = VecZeroEntries(Floc);CHKERRQ(ierr);
192     ierr = DMDAVecGetArray(dm,Floc,&f);CHKERRQ(ierr);
193     CHKMEMQ;
194     ierr = (*dmdasnes->residuallocal)(&info,x,f,dmdasnes->residuallocalctx);CHKERRQ(ierr);
195     CHKMEMQ;
196     ierr = DMDAVecRestoreArray(dm,Floc,&f);CHKERRQ(ierr);
197     ierr = VecZeroEntries(F);CHKERRQ(ierr);
198     ierr = DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr);
199     ierr = DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr);
200     ierr = DMRestoreLocalVector(dm,&Floc);CHKERRQ(ierr);
201   } break;
202   default: SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdasnes->residuallocalimode);
203   }
204   ierr = DMDAVecRestoreArray(dm,X,&x);CHKERRQ(ierr);
205   PetscFunctionReturn(0);
206 }
207 
208 #undef __FUNCT__
209 #define __FUNCT__ "SNESComputeJacobian_DMDA"
210 static PetscErrorCode SNESComputeJacobian_DMDA(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *mstr,void *ctx)
211 {
212   PetscErrorCode ierr;
213   DM             dm;
214   DMSNES_DA      *dmdasnes = (DMSNES_DA *)ctx;
215   DMDALocalInfo  info;
216   Vec            Xloc;
217   void           *x;
218 
219   PetscFunctionBegin;
220   if (!dmdasnes->residuallocal) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_PLIB,"Corrupt context");
221   ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);
222 
223   if (dmdasnes->jacobianlocal) {
224     ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr);
225     ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
226     ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
227     ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);
228     ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr);
229     CHKMEMQ;
230     ierr = (*dmdasnes->jacobianlocal)(&info,x,*A,*B,mstr,dmdasnes->jacobianlocalctx);CHKERRQ(ierr);
231     CHKMEMQ;
232     ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr);
233     ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr);
234   } else {
235     MatFDColoring fdcoloring;
236     ierr = PetscObjectQuery((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject*)&fdcoloring);CHKERRQ(ierr);
237     if (!fdcoloring) {
238       ISColoring     coloring;
239 
240       ierr = DMCreateColoring(dm,dm->coloringtype,dm->mattype,&coloring);CHKERRQ(ierr);
241       ierr = MatFDColoringCreate(*B,coloring,&fdcoloring);CHKERRQ(ierr);
242       ierr = ISColoringDestroy(&coloring);CHKERRQ(ierr);
243       switch (dm->coloringtype) {
244       case IS_COLORING_GLOBAL:
245         ierr = MatFDColoringSetFunction(fdcoloring,(PetscErrorCode(*)(void))SNESComputeFunction_DMDA,dmdasnes);CHKERRQ(ierr);
246         break;
247       default: SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_SUP,"No support for coloring type '%s'",ISColoringTypes[dm->coloringtype]);
248       }
249       ierr = PetscObjectSetOptionsPrefix((PetscObject)fdcoloring,((PetscObject)dm)->prefix);CHKERRQ(ierr);
250       ierr = MatFDColoringSetFromOptions(fdcoloring);CHKERRQ(ierr);
251       ierr = PetscObjectCompose((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject)fdcoloring);CHKERRQ(ierr);
252       ierr = PetscObjectDereference((PetscObject)fdcoloring);CHKERRQ(ierr);
253 
254       /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call
255        * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the
256        * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually
257        * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be
258        * taken when the PetscOList for the Vec inside MatFDColoring is destroyed.
259        */
260       ierr = PetscObjectDereference((PetscObject)dm);CHKERRQ(ierr);
261     }
262     *mstr = SAME_NONZERO_PATTERN;
263     ierr = MatFDColoringApply(*B,fdcoloring,X,mstr,snes);CHKERRQ(ierr);
264   }
265   /* This will be redundant if the user called both, but it's too common to forget. */
266   if (*A != *B) {
267     ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
268     ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
269   }
270   PetscFunctionReturn(0);
271 }
272 
273 #undef __FUNCT__
274 #define __FUNCT__ "SNESComputeLocalBlockJacobian_DMDA"
275 PetscErrorCode SNESComputeLocalBlockJacobian_DMDA(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *mstr,void *dmctx)
276 {
277   PetscErrorCode ierr;
278   DM             dm = (DM)dmctx; /* the global DMDA context */
279   DMSNES_DA      *dmdasnes;
280   DMSNES         sdm;
281   DMDALocalInfo  info;
282   void           *x;
283 
284   PetscFunctionBegin;
285   ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr);
286   ierr = DMDASNESGetContext(dm,sdm,&dmdasnes);CHKERRQ(ierr);
287   if (!dmdasnes->residuallocal) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_PLIB,"Corrupt context");
288   if (dmdasnes->jacobianlocal) {
289     ierr = DMDAVecGetArray(dm,X,&x);CHKERRQ(ierr);
290     CHKMEMQ;
291     ierr = (*dmdasnes->jacobianlocal)(&info,x,*A,*B,mstr,dmdasnes->jacobianlocalctx);CHKERRQ(ierr);
292     CHKMEMQ;
293     ierr = DMDAVecRestoreArray(dm,X,&x);CHKERRQ(ierr);
294   }
295   /* This will be redundant if the user called both, but it's too common to forget. */
296   if (*A != *B) {
297     ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
298     ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
299   }
300   PetscFunctionReturn(0);
301 }
302 
303 #undef __FUNCT__
304 #define __FUNCT__ "DMDASNESSetFunctionLocal"
305 /*@C
306    DMDASNESSetFunctionLocal - set a local residual evaluation function
307 
308    Logically Collective
309 
310    Input Arguments:
311 +  dm - DM to associate callback with
312 .  imode - INSERT_VALUES if local function computes owned part, ADD_VALUES if it contributes to ghosted part
313 .  func - local residual evaluation
314 -  ctx - optional context for local residual evaluation
315 
316    Calling sequence for func:
317 +  info - DMDALocalInfo defining the subdomain to evaluate the residual on
318 .  x - dimensional pointer to state at which to evaluate residual
319 .  f - dimensional pointer to residual, write the residual here
320 -  ctx - optional context passed above
321 
322    Level: beginner
323 
324 .seealso: DMSNESSetFunction(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
325 @*/
326 PetscErrorCode DMDASNESSetFunctionLocal(DM dm,InsertMode imode,PetscErrorCode (*func)(DMDALocalInfo*,void*,void*,void*),void *ctx)
327 {
328   PetscErrorCode ierr;
329   DMSNES         sdm;
330   DMSNES_DA      *dmdasnes;
331 
332   PetscFunctionBegin;
333   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
334   ierr = DMGetDMSNESWrite(dm,&sdm);CHKERRQ(ierr);
335   ierr = DMDASNESGetContext(dm,sdm,&dmdasnes);CHKERRQ(ierr);
336   dmdasnes->residuallocalimode = imode;
337   dmdasnes->residuallocal = func;
338   dmdasnes->residuallocalctx = ctx;
339   ierr = DMSNESSetFunction(dm,SNESComputeFunction_DMDA,dmdasnes);CHKERRQ(ierr);
340   ierr = DMSNESSetBlockFunction(dm,SNESComputeLocalBlockFunction_DMDA,dm);CHKERRQ(ierr);
341   if (!sdm->ops->computejacobian) {  /* Call us for the Jacobian too, can be overridden by the user. */
342     ierr = DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes);CHKERRQ(ierr);
343   }
344   PetscFunctionReturn(0);
345 }
346 
347 #undef __FUNCT__
348 #define __FUNCT__ "DMDASNESSetJacobianLocal"
349 /*@C
350    DMDASNESSetJacobianLocal - set a local residual evaluation function
351 
352    Logically Collective
353 
354    Input Arguments:
355 +  dm - DM to associate callback with
356 .  func - local residual evaluation
357 -  ctx - optional context for local residual evaluation
358 
359    Calling sequence for func:
360 +  info - DMDALocalInfo defining the subdomain to evaluate the residual on
361 .  x - dimensional pointer to state at which to evaluate residual
362 .  f - dimensional pointer to residual, write the residual here
363 -  ctx - optional context passed above
364 
365    Level: beginner
366 
367 .seealso: DMSNESSetJacobian(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
368 @*/
369 PetscErrorCode DMDASNESSetJacobianLocal(DM dm,PetscErrorCode (*func)(DMDALocalInfo*,void*,Mat,Mat,MatStructure*,void*),void *ctx)
370 {
371   PetscErrorCode ierr;
372   DMSNES         sdm;
373   DMSNES_DA      *dmdasnes;
374 
375   PetscFunctionBegin;
376   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
377   ierr = DMGetDMSNESWrite(dm,&sdm);CHKERRQ(ierr);
378   ierr = DMDASNESGetContext(dm,sdm,&dmdasnes);CHKERRQ(ierr);
379   dmdasnes->jacobianlocal = func;
380   dmdasnes->jacobianlocalctx = ctx;
381   ierr = DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes);CHKERRQ(ierr);
382   ierr = DMSNESSetBlockJacobian(dm,SNESComputeLocalBlockJacobian_DMDA,dm);CHKERRQ(ierr);
383   PetscFunctionReturn(0);
384 }
385 
386 
387 #undef __FUNCT__
388 #define __FUNCT__ "DMDASNESSetObjectiveLocal"
389 /*@C
390    DMDASNESSetObjectiveLocal - set a local residual evaluation function
391 
392    Logically Collective
393 
394    Input Arguments:
395 +  dm - DM to associate callback with
396 .  func - local objective evaluation
397 -  ctx - optional context for local residual evaluation
398 
399    Calling sequence for func:
400 +  info - DMDALocalInfo defining the subdomain to evaluate the residual on
401 .  x - dimensional pointer to state at which to evaluate residual
402 .  ob - eventual objective value
403 -  ctx - optional context passed above
404 
405    Level: beginner
406 
407 .seealso: DMSNESSetFunction(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
408 @*/
409 PetscErrorCode DMDASNESSetObjectiveLocal(DM dm,DMDASNESObjective func,void *ctx)
410 {
411   PetscErrorCode ierr;
412   DMSNES         sdm;
413   DMSNES_DA      *dmdasnes;
414 
415   PetscFunctionBegin;
416   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
417   ierr = DMGetDMSNESWrite(dm,&sdm);CHKERRQ(ierr);
418   ierr = DMDASNESGetContext(dm,sdm,&dmdasnes);CHKERRQ(ierr);
419   dmdasnes->objectivelocal = func;
420   dmdasnes->objectivelocalctx = ctx;
421   ierr = DMSNESSetObjective(dm,SNESComputeObjective_DMDA,dmdasnes);CHKERRQ(ierr);
422   PetscFunctionReturn(0);
423 }
424 
425 #undef __FUNCT__
426 #define __FUNCT__ "SNESComputePicard_DMDA"
427 static PetscErrorCode SNESComputePicard_DMDA(SNES snes,Vec X,Vec F,void *ctx)
428 {
429   PetscErrorCode ierr;
430   DM             dm;
431   DMSNES_DA      *dmdasnes = (DMSNES_DA *)ctx;
432   DMDALocalInfo  info;
433   Vec            Xloc;
434   void           *x,*f;
435 
436   PetscFunctionBegin;
437   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
438   PetscValidHeaderSpecific(X,VEC_CLASSID,2);
439   PetscValidHeaderSpecific(F,VEC_CLASSID,3);
440   if (!dmdasnes->rhsplocal) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_PLIB,"Corrupt context");
441   ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);
442   ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr);
443   ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
444   ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
445   ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);
446   ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr);
447   switch (dmdasnes->residuallocalimode) {
448   case INSERT_VALUES: {
449     ierr = DMDAVecGetArray(dm,F,&f);CHKERRQ(ierr);
450     CHKMEMQ;
451     ierr = (*dmdasnes->rhsplocal)(&info,x,f,dmdasnes->picardlocalctx);CHKERRQ(ierr);
452     CHKMEMQ;
453     ierr = DMDAVecRestoreArray(dm,F,&f);CHKERRQ(ierr);
454   } break;
455   case ADD_VALUES: {
456     Vec Floc;
457     ierr = DMGetLocalVector(dm,&Floc);CHKERRQ(ierr);
458     ierr = VecZeroEntries(Floc);CHKERRQ(ierr);
459     ierr = DMDAVecGetArray(dm,Floc,&f);CHKERRQ(ierr);
460     CHKMEMQ;
461     ierr = (*dmdasnes->rhsplocal)(&info,x,f,dmdasnes->picardlocalctx);CHKERRQ(ierr);
462     CHKMEMQ;
463     ierr = DMDAVecRestoreArray(dm,Floc,&f);CHKERRQ(ierr);
464     ierr = VecZeroEntries(F);CHKERRQ(ierr);
465     ierr = DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr);
466     ierr = DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr);
467     ierr = DMRestoreLocalVector(dm,&Floc);CHKERRQ(ierr);
468   } break;
469   default: SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdasnes->residuallocalimode);
470   }
471   ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr);
472   ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr);
473   PetscFunctionReturn(0);
474 }
475 
476 #undef __FUNCT__
477 #define __FUNCT__ "SNESComputePicardJacobian_DMDA"
478 static PetscErrorCode SNESComputePicardJacobian_DMDA(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *mstr,void *ctx)
479 {
480   PetscErrorCode ierr;
481   DM             dm;
482   DMSNES_DA      *dmdasnes = (DMSNES_DA *)ctx;
483   DMDALocalInfo  info;
484   Vec            Xloc;
485   void           *x;
486 
487   PetscFunctionBegin;
488   if (!dmdasnes->jacobianplocal) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_PLIB,"Corrupt context");
489   ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);
490 
491   ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr);
492   ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
493   ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
494   ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);
495   ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr);
496   CHKMEMQ;
497   ierr = (*dmdasnes->jacobianplocal)(&info,x,*A,*B,mstr,dmdasnes->picardlocalctx);CHKERRQ(ierr);
498   CHKMEMQ;
499   ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr);
500   ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr);
501   *mstr = SAME_NONZERO_PATTERN;
502   if (*A != *B) {
503     ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
504     ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
505   }
506   PetscFunctionReturn(0);
507 }
508 
509 #undef __FUNCT__
510 #define __FUNCT__ "DMDASNESSetPicardLocal"
511 /*@C
512    DMDASNESSetPicardLocal - set a local right hand side and matrix evaluation function for Picard iteration
513 
514    Logically Collective
515 
516    Input Arguments:
517 +  dm - DM to associate callback with
518 .  imode - INSERT_VALUES if local function computes owned part, ADD_VALUES if it contributes to ghosted part
519 .  func - local residual evaluation
520 -  ctx - optional context for local residual evaluation
521 
522    Calling sequence for func:
523 +  info - DMDALocalInfo defining the subdomain to evaluate the residual on
524 .  x - dimensional pointer to state at which to evaluate residual
525 .  f - dimensional pointer to residual, write the residual here
526 -  ctx - optional context passed above
527 
528    Notes:  The user must use
529     extern PetscErrorCode  SNESPicardComputeFunction(SNES,Vec,Vec,void*);
530     extern PetscErrorCode  SNESPicardComputeJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
531     ierr = SNESSetFunction(snes,PETSC_NULL,SNESPicardComputeFunction,&user);CHKERRQ(ierr);
532     in their code before calling this routine.
533 
534 
535    Level: beginner
536 
537 .seealso: DMSNESSetFunction(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
538 @*/
539 PetscErrorCode DMDASNESSetPicardLocal(DM dm,InsertMode imode,PetscErrorCode (*func)(DMDALocalInfo*,void*,void*,void*),
540                                         PetscErrorCode (*jac)(DMDALocalInfo*,void*,Mat,Mat,MatStructure*,void*),void *ctx)
541 {
542   PetscErrorCode ierr;
543   DMSNES         sdm;
544   DMSNES_DA      *dmdasnes;
545 
546   PetscFunctionBegin;
547   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
548   ierr = DMGetDMSNESWrite(dm,&sdm);CHKERRQ(ierr);
549   ierr = DMDASNESGetContext(dm,sdm,&dmdasnes);CHKERRQ(ierr);
550   dmdasnes->residuallocalimode = imode;
551   dmdasnes->rhsplocal          = func;
552   dmdasnes->jacobianplocal     = jac;
553   dmdasnes->picardlocalctx     = ctx;
554   ierr = DMSNESSetPicard(dm,SNESComputePicard_DMDA,SNESComputePicardJacobian_DMDA,dmdasnes);CHKERRQ(ierr);
555   PetscFunctionReturn(0);
556 }
557