Actual source code: dmdasnes.c
petsc-3.4.5 2014-06-29
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,MatStructure*,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,MatStructure*,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: /*
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: {
73: DM dm;
74: DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx;
75: DMDALocalInfo info;
76: Vec Xloc;
77: void *x,*f;
83: if (!dmdasnes->residuallocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
84: SNESGetDM(snes,&dm);
85: DMGetLocalVector(dm,&Xloc);
86: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
87: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
88: DMDAGetLocalInfo(dm,&info);
89: DMDAVecGetArray(dm,Xloc,&x);
90: switch (dmdasnes->residuallocalimode) {
91: case INSERT_VALUES: {
92: DMDAVecGetArray(dm,F,&f);
93: CHKMEMQ;
94: (*dmdasnes->residuallocal)(&info,x,f,dmdasnes->residuallocalctx);
95: CHKMEMQ;
96: DMDAVecRestoreArray(dm,F,&f);
97: } break;
98: case ADD_VALUES: {
99: Vec Floc;
100: DMGetLocalVector(dm,&Floc);
101: VecZeroEntries(Floc);
102: DMDAVecGetArray(dm,Floc,&f);
103: CHKMEMQ;
104: (*dmdasnes->residuallocal)(&info,x,f,dmdasnes->residuallocalctx);
105: CHKMEMQ;
106: DMDAVecRestoreArray(dm,Floc,&f);
107: VecZeroEntries(F);
108: DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
109: DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
110: DMRestoreLocalVector(dm,&Floc);
111: } break;
112: default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdasnes->residuallocalimode);
113: }
114: DMDAVecRestoreArray(dm,Xloc,&x);
115: DMRestoreLocalVector(dm,&Xloc);
116: return(0);
117: }
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: {
128: DM dm;
129: DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx;
130: DMDALocalInfo info;
131: Vec Xloc;
132: void *x;
138: if (!dmdasnes->objectivelocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
139: SNESGetDM(snes,&dm);
140: DMGetLocalVector(dm,&Xloc);
141: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
142: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
143: DMDAGetLocalInfo(dm,&info);
144: DMDAVecGetArray(dm,Xloc,&x);
145: CHKMEMQ;
146: (*dmdasnes->objectivelocal)(&info,x,ob,dmdasnes->objectivelocalctx);
147: CHKMEMQ;
148: DMDAVecRestoreArray(dm,Xloc,&x);
149: DMRestoreLocalVector(dm,&Xloc);
150: return(0);
151: }
156: static PetscErrorCode SNESComputeJacobian_DMDA(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *mstr,void *ctx)
157: {
159: DM dm;
160: DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx;
161: DMDALocalInfo info;
162: Vec Xloc;
163: void *x;
166: if (!dmdasnes->residuallocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
167: SNESGetDM(snes,&dm);
169: if (dmdasnes->jacobianlocal) {
170: DMGetLocalVector(dm,&Xloc);
171: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
172: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
173: DMDAGetLocalInfo(dm,&info);
174: DMDAVecGetArray(dm,Xloc,&x);
175: CHKMEMQ;
176: (*dmdasnes->jacobianlocal)(&info,x,*A,*B,mstr,dmdasnes->jacobianlocalctx);
177: CHKMEMQ;
178: DMDAVecRestoreArray(dm,Xloc,&x);
179: DMRestoreLocalVector(dm,&Xloc);
180: } else {
181: MatFDColoring fdcoloring;
182: PetscObjectQuery((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject*)&fdcoloring);
183: if (!fdcoloring) {
184: ISColoring coloring;
186: DMCreateColoring(dm,dm->coloringtype,dm->mattype,&coloring);
187: MatFDColoringCreate(*B,coloring,&fdcoloring);
188: ISColoringDestroy(&coloring);
189: switch (dm->coloringtype) {
190: case IS_COLORING_GLOBAL:
191: MatFDColoringSetFunction(fdcoloring,(PetscErrorCode (*)(void))SNESComputeFunction_DMDA,dmdasnes);
192: break;
193: default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"No support for coloring type '%s'",ISColoringTypes[dm->coloringtype]);
194: }
195: PetscObjectSetOptionsPrefix((PetscObject)fdcoloring,((PetscObject)dm)->prefix);
196: MatFDColoringSetFromOptions(fdcoloring);
197: PetscObjectCompose((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject)fdcoloring);
198: PetscObjectDereference((PetscObject)fdcoloring);
200: /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call
201: * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the
202: * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually
203: * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be
204: * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed.
205: */
206: PetscObjectDereference((PetscObject)dm);
207: }
208: *mstr = SAME_NONZERO_PATTERN;
209: MatFDColoringApply(*B,fdcoloring,X,mstr,snes);
210: }
211: /* This will be redundant if the user called both, but it's too common to forget. */
212: if (*A != *B) {
213: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
214: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
215: }
216: return(0);
217: }
221: /*@C
222: DMDASNESSetFunctionLocal - set a local residual evaluation function
224: Logically Collective
226: Input Arguments:
227: + dm - DM to associate callback with
228: . imode - INSERT_VALUES if local function computes owned part, ADD_VALUES if it contributes to ghosted part
229: . func - local residual evaluation
230: - ctx - optional context for local residual evaluation
232: Calling sequence for func:
233: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
234: . x - dimensional pointer to state at which to evaluate residual
235: . f - dimensional pointer to residual, write the residual here
236: - ctx - optional context passed above
238: Level: beginner
240: .seealso: DMSNESSetFunction(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
241: @*/
242: PetscErrorCode DMDASNESSetFunctionLocal(DM dm,InsertMode imode,PetscErrorCode (*func)(DMDALocalInfo*,void*,void*,void*),void *ctx)
243: {
245: DMSNES sdm;
246: DMSNES_DA *dmdasnes;
250: DMGetDMSNESWrite(dm,&sdm);
251: DMDASNESGetContext(dm,sdm,&dmdasnes);
253: dmdasnes->residuallocalimode = imode;
254: dmdasnes->residuallocal = func;
255: dmdasnes->residuallocalctx = ctx;
257: DMSNESSetFunction(dm,SNESComputeFunction_DMDA,dmdasnes);
258: if (!sdm->ops->computejacobian) { /* Call us for the Jacobian too, can be overridden by the user. */
259: DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes);
260: }
261: return(0);
262: }
266: /*@C
267: DMDASNESSetJacobianLocal - set a local Jacobian evaluation function
269: Logically Collective
271: Input Arguments:
272: + dm - DM to associate callback with
273: . func - local residual evaluation
274: - ctx - optional context for local residual evaluation
276: Calling sequence for func:
277: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
278: . x - dimensional pointer to state at which to evaluate residual
279: . f - dimensional pointer to residual, write the residual here
280: - ctx - optional context passed above
282: Level: beginner
284: .seealso: DMSNESSetJacobian(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
285: @*/
286: PetscErrorCode DMDASNESSetJacobianLocal(DM dm,PetscErrorCode (*func)(DMDALocalInfo*,void*,Mat,Mat,MatStructure*,void*),void *ctx)
287: {
289: DMSNES sdm;
290: DMSNES_DA *dmdasnes;
294: DMGetDMSNESWrite(dm,&sdm);
295: DMDASNESGetContext(dm,sdm,&dmdasnes);
297: dmdasnes->jacobianlocal = func;
298: dmdasnes->jacobianlocalctx = ctx;
300: DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes);
301: return(0);
302: }
307: /*@C
308: DMDASNESSetObjectiveLocal - set a local residual evaluation function
310: Logically Collective
312: Input Arguments:
313: + dm - DM to associate callback with
314: . func - local objective evaluation
315: - ctx - optional context for local residual evaluation
317: Calling sequence for func:
318: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
319: . x - dimensional pointer to state at which to evaluate residual
320: . ob - eventual objective value
321: - ctx - optional context passed above
323: Level: beginner
325: .seealso: DMSNESSetFunction(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
326: @*/
327: PetscErrorCode DMDASNESSetObjectiveLocal(DM dm,DMDASNESObjective func,void *ctx)
328: {
330: DMSNES sdm;
331: DMSNES_DA *dmdasnes;
335: DMGetDMSNESWrite(dm,&sdm);
336: DMDASNESGetContext(dm,sdm,&dmdasnes);
338: dmdasnes->objectivelocal = func;
339: dmdasnes->objectivelocalctx = ctx;
341: DMSNESSetObjective(dm,SNESComputeObjective_DMDA,dmdasnes);
342: return(0);
343: }
347: static PetscErrorCode SNESComputePicard_DMDA(SNES snes,Vec X,Vec F,void *ctx)
348: {
350: DM dm;
351: DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx;
352: DMDALocalInfo info;
353: Vec Xloc;
354: void *x,*f;
360: if (!dmdasnes->rhsplocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
361: SNESGetDM(snes,&dm);
362: DMGetLocalVector(dm,&Xloc);
363: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
364: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
365: DMDAGetLocalInfo(dm,&info);
366: DMDAVecGetArray(dm,Xloc,&x);
367: switch (dmdasnes->residuallocalimode) {
368: case INSERT_VALUES: {
369: DMDAVecGetArray(dm,F,&f);
370: CHKMEMQ;
371: (*dmdasnes->rhsplocal)(&info,x,f,dmdasnes->picardlocalctx);
372: CHKMEMQ;
373: DMDAVecRestoreArray(dm,F,&f);
374: } break;
375: case ADD_VALUES: {
376: Vec Floc;
377: DMGetLocalVector(dm,&Floc);
378: VecZeroEntries(Floc);
379: DMDAVecGetArray(dm,Floc,&f);
380: CHKMEMQ;
381: (*dmdasnes->rhsplocal)(&info,x,f,dmdasnes->picardlocalctx);
382: CHKMEMQ;
383: DMDAVecRestoreArray(dm,Floc,&f);
384: VecZeroEntries(F);
385: DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
386: DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
387: DMRestoreLocalVector(dm,&Floc);
388: } break;
389: default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdasnes->residuallocalimode);
390: }
391: DMDAVecRestoreArray(dm,Xloc,&x);
392: DMRestoreLocalVector(dm,&Xloc);
393: return(0);
394: }
398: static PetscErrorCode SNESComputePicardJacobian_DMDA(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *mstr,void *ctx)
399: {
401: DM dm;
402: DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx;
403: DMDALocalInfo info;
404: Vec Xloc;
405: void *x;
408: if (!dmdasnes->jacobianplocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
409: SNESGetDM(snes,&dm);
411: DMGetLocalVector(dm,&Xloc);
412: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
413: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
414: DMDAGetLocalInfo(dm,&info);
415: DMDAVecGetArray(dm,Xloc,&x);
416: CHKMEMQ;
417: (*dmdasnes->jacobianplocal)(&info,x,*A,*B,mstr,dmdasnes->picardlocalctx);
418: CHKMEMQ;
419: DMDAVecRestoreArray(dm,Xloc,&x);
420: DMRestoreLocalVector(dm,&Xloc);
421: *mstr = SAME_NONZERO_PATTERN;
422: if (*A != *B) {
423: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
424: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
425: }
426: return(0);
427: }
431: /*@C
432: DMDASNESSetPicardLocal - set a local right hand side and matrix evaluation function for Picard iteration
434: Logically Collective
436: Input Arguments:
437: + dm - DM to associate callback with
438: . imode - INSERT_VALUES if local function computes owned part, ADD_VALUES if it contributes to ghosted part
439: . func - local residual evaluation
440: - ctx - optional context for local residual evaluation
442: Calling sequence for func:
443: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
444: . x - dimensional pointer to state at which to evaluate residual
445: . f - dimensional pointer to residual, write the residual here
446: - ctx - optional context passed above
448: Notes: The user must use
449: extern PetscErrorCode SNESPicardComputeFunction(SNES,Vec,Vec,void*);
450: extern PetscErrorCode SNESPicardComputeJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
451: SNESSetFunction(snes,NULL,SNESPicardComputeFunction,&user);
452: in their code before calling this routine.
455: Level: beginner
457: .seealso: DMSNESSetFunction(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
458: @*/
459: PetscErrorCode DMDASNESSetPicardLocal(DM dm,InsertMode imode,PetscErrorCode (*func)(DMDALocalInfo*,void*,void*,void*),
460: PetscErrorCode (*jac)(DMDALocalInfo*,void*,Mat,Mat,MatStructure*,void*),void *ctx)
461: {
463: DMSNES sdm;
464: DMSNES_DA *dmdasnes;
468: DMGetDMSNESWrite(dm,&sdm);
469: DMDASNESGetContext(dm,sdm,&dmdasnes);
471: dmdasnes->residuallocalimode = imode;
472: dmdasnes->rhsplocal = func;
473: dmdasnes->jacobianplocal = jac;
474: dmdasnes->picardlocalctx = ctx;
476: DMSNESSetPicard(dm,SNESComputePicard_DMDA,SNESComputePicardJacobian_DMDA,dmdasnes);
477: return(0);
478: }