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