Actual source code: dmdasnes.c
petsc-3.5.4 2015-05-23
1: #include <petscdmda.h> /*I "petscdmda.h" I*/
2: #include <petsc-private/dmimpl.h>
3: #include <petsc-private/snesimpl.h> /*I "petscsnes.h" I*/
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,void*);
9: PetscErrorCode (*objectivelocal)(DMDALocalInfo*,void*,PetscReal*,void*);
10: void *residuallocalctx;
11: void *jacobianlocalctx;
12: void *objectivelocalctx;
13: InsertMode residuallocalimode;
15: /* For Picard iteration defined locally */
16: PetscErrorCode (*rhsplocal)(DMDALocalInfo*,void*,void*,void*);
17: PetscErrorCode (*jacobianplocal)(DMDALocalInfo*,void*,Mat,Mat,void*);
18: void *picardlocalctx;
19: } DMSNES_DA;
23: static PetscErrorCode DMSNESDestroy_DMDA(DMSNES sdm)
24: {
28: PetscFree(sdm->data);
29: return(0);
30: }
34: static PetscErrorCode DMSNESDuplicate_DMDA(DMSNES oldsdm,DMSNES sdm)
35: {
39: PetscNewLog(sdm,(DMSNES_DA**)&sdm->data);
40: if (oldsdm->data) {
41: PetscMemcpy(sdm->data,oldsdm->data,sizeof(DMSNES_DA));
42: }
43: return(0);
44: }
49: static PetscErrorCode DMDASNESGetContext(DM dm,DMSNES sdm,DMSNES_DA **dmdasnes)
50: {
54: *dmdasnes = NULL;
55: if (!sdm->data) {
56: PetscNewLog(dm,(DMSNES_DA**)&sdm->data);
57: sdm->ops->destroy = DMSNESDestroy_DMDA;
58: sdm->ops->duplicate = DMSNESDuplicate_DMDA;
59: }
60: *dmdasnes = (DMSNES_DA*)sdm->data;
61: return(0);
62: }
66: static PetscErrorCode SNESComputeFunction_DMDA(SNES snes,Vec X,Vec F,void *ctx)
67: {
69: DM dm;
70: DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx;
71: DMDALocalInfo info;
72: Vec Xloc;
73: void *x,*f;
79: if (!dmdasnes->residuallocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
80: SNESGetDM(snes,&dm);
81: DMGetLocalVector(dm,&Xloc);
82: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
83: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
84: DMDAGetLocalInfo(dm,&info);
85: DMDAVecGetArray(dm,Xloc,&x);
86: switch (dmdasnes->residuallocalimode) {
87: case INSERT_VALUES: {
88: DMDAVecGetArray(dm,F,&f);
89: CHKMEMQ;
90: (*dmdasnes->residuallocal)(&info,x,f,dmdasnes->residuallocalctx);
91: CHKMEMQ;
92: DMDAVecRestoreArray(dm,F,&f);
93: } break;
94: case ADD_VALUES: {
95: Vec Floc;
96: DMGetLocalVector(dm,&Floc);
97: VecZeroEntries(Floc);
98: DMDAVecGetArray(dm,Floc,&f);
99: CHKMEMQ;
100: (*dmdasnes->residuallocal)(&info,x,f,dmdasnes->residuallocalctx);
101: CHKMEMQ;
102: DMDAVecRestoreArray(dm,Floc,&f);
103: VecZeroEntries(F);
104: DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
105: DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
106: DMRestoreLocalVector(dm,&Floc);
107: } break;
108: default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdasnes->residuallocalimode);
109: }
110: DMDAVecRestoreArray(dm,Xloc,&x);
111: DMRestoreLocalVector(dm,&Xloc);
112: return(0);
113: }
117: static PetscErrorCode SNESComputeObjective_DMDA(SNES snes,Vec X,PetscReal *ob,void *ctx)
118: {
120: DM dm;
121: DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx;
122: DMDALocalInfo info;
123: Vec Xloc;
124: void *x;
130: if (!dmdasnes->objectivelocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
131: SNESGetDM(snes,&dm);
132: DMGetLocalVector(dm,&Xloc);
133: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
134: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
135: DMDAGetLocalInfo(dm,&info);
136: DMDAVecGetArray(dm,Xloc,&x);
137: CHKMEMQ;
138: (*dmdasnes->objectivelocal)(&info,x,ob,dmdasnes->objectivelocalctx);
139: CHKMEMQ;
140: DMDAVecRestoreArray(dm,Xloc,&x);
141: DMRestoreLocalVector(dm,&Xloc);
142: return(0);
143: }
148: static PetscErrorCode SNESComputeJacobian_DMDA(SNES snes,Vec X,Mat A,Mat B,void *ctx)
149: {
151: DM dm;
152: DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx;
153: DMDALocalInfo info;
154: Vec Xloc;
155: void *x;
158: if (!dmdasnes->residuallocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
159: SNESGetDM(snes,&dm);
161: if (dmdasnes->jacobianlocal) {
162: DMGetLocalVector(dm,&Xloc);
163: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
164: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
165: DMDAGetLocalInfo(dm,&info);
166: DMDAVecGetArray(dm,Xloc,&x);
167: CHKMEMQ;
168: (*dmdasnes->jacobianlocal)(&info,x,A,B,dmdasnes->jacobianlocalctx);
169: CHKMEMQ;
170: DMDAVecRestoreArray(dm,Xloc,&x);
171: DMRestoreLocalVector(dm,&Xloc);
172: } else {
173: MatFDColoring fdcoloring;
174: PetscObjectQuery((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject*)&fdcoloring);
175: if (!fdcoloring) {
176: ISColoring coloring;
178: DMCreateColoring(dm,dm->coloringtype,&coloring);
179: MatFDColoringCreate(B,coloring,&fdcoloring);
180: switch (dm->coloringtype) {
181: case IS_COLORING_GLOBAL:
182: MatFDColoringSetFunction(fdcoloring,(PetscErrorCode (*)(void))SNESComputeFunction_DMDA,dmdasnes);
183: break;
184: default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"No support for coloring type '%s'",ISColoringTypes[dm->coloringtype]);
185: }
186: PetscObjectSetOptionsPrefix((PetscObject)fdcoloring,((PetscObject)dm)->prefix);
187: MatFDColoringSetFromOptions(fdcoloring);
188: MatFDColoringSetUp(B,coloring,fdcoloring);
189: ISColoringDestroy(&coloring);
190: PetscObjectCompose((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject)fdcoloring);
191: PetscObjectDereference((PetscObject)fdcoloring);
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: PetscObjectDereference((PetscObject)dm);
200: }
201: 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: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
206: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
207: }
208: return(0);
209: }
213: /*@C
214: DMDASNESSetFunctionLocal - set a local residual evaluation function
216: Logically Collective
218: Input Arguments:
219: + dm - DM to associate callback with
220: . imode - INSERT_VALUES if local function computes owned part, ADD_VALUES if it contributes to ghosted part
221: . func - local residual evaluation
222: - ctx - optional context for local residual evaluation
224: Calling sequence for func:
225: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
226: . x - dimensional pointer to state at which to evaluate residual
227: . f - dimensional pointer to residual, write the residual here
228: - ctx - optional context passed above
230: Level: beginner
232: .seealso: DMSNESSetFunction(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
233: @*/
234: PetscErrorCode DMDASNESSetFunctionLocal(DM dm,InsertMode imode,PetscErrorCode (*func)(DMDALocalInfo*,void*,void*,void*),void *ctx)
235: {
237: DMSNES sdm;
238: DMSNES_DA *dmdasnes;
242: DMGetDMSNESWrite(dm,&sdm);
243: DMDASNESGetContext(dm,sdm,&dmdasnes);
245: dmdasnes->residuallocalimode = imode;
246: dmdasnes->residuallocal = func;
247: dmdasnes->residuallocalctx = ctx;
249: DMSNESSetFunction(dm,SNESComputeFunction_DMDA,dmdasnes);
250: if (!sdm->ops->computejacobian) { /* Call us for the Jacobian too, can be overridden by the user. */
251: DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes);
252: }
253: return(0);
254: }
258: /*@C
259: DMDASNESSetJacobianLocal - set a local Jacobian evaluation function
261: Logically Collective
263: Input Arguments:
264: + dm - DM to associate callback with
265: . func - local residual evaluation
266: - ctx - optional context for local residual evaluation
268: Calling sequence for func:
269: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
270: . x - dimensional pointer to state at which to evaluate residual
271: . f - dimensional pointer to residual, write the residual here
272: - ctx - optional context passed above
274: Level: beginner
276: .seealso: DMSNESSetJacobian(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
277: @*/
278: PetscErrorCode DMDASNESSetJacobianLocal(DM dm,PetscErrorCode (*func)(DMDALocalInfo*,void*,Mat,Mat,void*),void *ctx)
279: {
281: DMSNES sdm;
282: DMSNES_DA *dmdasnes;
286: DMGetDMSNESWrite(dm,&sdm);
287: DMDASNESGetContext(dm,sdm,&dmdasnes);
289: dmdasnes->jacobianlocal = func;
290: dmdasnes->jacobianlocalctx = ctx;
292: DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes);
293: return(0);
294: }
299: /*@C
300: DMDASNESSetObjectiveLocal - set a local residual evaluation function
302: Logically Collective
304: Input Arguments:
305: + dm - DM to associate callback with
306: . func - local objective evaluation
307: - ctx - optional context for local residual evaluation
309: Calling sequence for func:
310: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
311: . x - dimensional pointer to state at which to evaluate residual
312: . ob - eventual objective value
313: - ctx - optional context passed above
315: Level: beginner
317: .seealso: DMSNESSetFunction(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
318: @*/
319: PetscErrorCode DMDASNESSetObjectiveLocal(DM dm,DMDASNESObjective func,void *ctx)
320: {
322: DMSNES sdm;
323: DMSNES_DA *dmdasnes;
327: DMGetDMSNESWrite(dm,&sdm);
328: DMDASNESGetContext(dm,sdm,&dmdasnes);
330: dmdasnes->objectivelocal = func;
331: dmdasnes->objectivelocalctx = ctx;
333: DMSNESSetObjective(dm,SNESComputeObjective_DMDA,dmdasnes);
334: return(0);
335: }
339: static PetscErrorCode SNESComputePicard_DMDA(SNES snes,Vec X,Vec F,void *ctx)
340: {
342: DM dm;
343: DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx;
344: DMDALocalInfo info;
345: Vec Xloc;
346: void *x,*f;
352: if (!dmdasnes->rhsplocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
353: SNESGetDM(snes,&dm);
354: DMGetLocalVector(dm,&Xloc);
355: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
356: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
357: DMDAGetLocalInfo(dm,&info);
358: DMDAVecGetArray(dm,Xloc,&x);
359: switch (dmdasnes->residuallocalimode) {
360: case INSERT_VALUES: {
361: DMDAVecGetArray(dm,F,&f);
362: CHKMEMQ;
363: (*dmdasnes->rhsplocal)(&info,x,f,dmdasnes->picardlocalctx);
364: CHKMEMQ;
365: DMDAVecRestoreArray(dm,F,&f);
366: } break;
367: case ADD_VALUES: {
368: Vec Floc;
369: DMGetLocalVector(dm,&Floc);
370: VecZeroEntries(Floc);
371: DMDAVecGetArray(dm,Floc,&f);
372: CHKMEMQ;
373: (*dmdasnes->rhsplocal)(&info,x,f,dmdasnes->picardlocalctx);
374: CHKMEMQ;
375: DMDAVecRestoreArray(dm,Floc,&f);
376: VecZeroEntries(F);
377: DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
378: DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
379: DMRestoreLocalVector(dm,&Floc);
380: } break;
381: default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdasnes->residuallocalimode);
382: }
383: DMDAVecRestoreArray(dm,Xloc,&x);
384: DMRestoreLocalVector(dm,&Xloc);
385: return(0);
386: }
390: static PetscErrorCode SNESComputePicardJacobian_DMDA(SNES snes,Vec X,Mat A,Mat B,void *ctx)
391: {
393: DM dm;
394: DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx;
395: DMDALocalInfo info;
396: Vec Xloc;
397: void *x;
400: if (!dmdasnes->jacobianplocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
401: SNESGetDM(snes,&dm);
403: DMGetLocalVector(dm,&Xloc);
404: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
405: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
406: DMDAGetLocalInfo(dm,&info);
407: DMDAVecGetArray(dm,Xloc,&x);
408: CHKMEMQ;
409: (*dmdasnes->jacobianplocal)(&info,x,A,B,dmdasnes->picardlocalctx);
410: CHKMEMQ;
411: DMDAVecRestoreArray(dm,Xloc,&x);
412: DMRestoreLocalVector(dm,&Xloc);
413: if (A != B) {
414: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
415: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
416: }
417: return(0);
418: }
422: /*@C
423: DMDASNESSetPicardLocal - set a local right hand side and matrix evaluation function for Picard iteration
425: Logically Collective
427: Input Arguments:
428: + dm - DM to associate callback with
429: . imode - INSERT_VALUES if local function computes owned part, ADD_VALUES if it contributes to ghosted part
430: . func - local residual evaluation
431: - ctx - optional context for local residual evaluation
433: Calling sequence for func:
434: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
435: . x - dimensional pointer to state at which to evaluate residual
436: . f - dimensional pointer to residual, write the residual here
437: - ctx - optional context passed above
439: Notes: The user must use
440: extern PetscErrorCode SNESPicardComputeFunction(SNES,Vec,Vec,void*);
441: extern PetscErrorCode SNESPicardComputeJacobian(SNES,Vec,Mat,Mat,MatStructure*,void*);
442: SNESSetFunction(snes,NULL,SNESPicardComputeFunction,&user);
443: in their code before calling this routine.
446: Level: beginner
448: .seealso: DMSNESSetFunction(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
449: @*/
450: PetscErrorCode DMDASNESSetPicardLocal(DM dm,InsertMode imode,PetscErrorCode (*func)(DMDALocalInfo*,void*,void*,void*),
451: PetscErrorCode (*jac)(DMDALocalInfo*,void*,Mat,Mat,void*),void *ctx)
452: {
454: DMSNES sdm;
455: DMSNES_DA *dmdasnes;
459: DMGetDMSNESWrite(dm,&sdm);
460: DMDASNESGetContext(dm,sdm,&dmdasnes);
462: dmdasnes->residuallocalimode = imode;
463: dmdasnes->rhsplocal = func;
464: dmdasnes->jacobianplocal = jac;
465: dmdasnes->picardlocalctx = ctx;
467: DMSNESSetPicard(dm,SNESComputePicard_DMDA,SNESComputePicardJacobian_DMDA,dmdasnes);
468: return(0);
469: }