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