Actual source code: dmdasnes.c
petsc-3.6.4 2016-04-12
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: PetscLogEventBegin(SNES_FunctionEval,snes,X,F,0);
90: CHKMEMQ;
91: (*dmdasnes->residuallocal)(&info,x,f,dmdasnes->residuallocalctx);
92: CHKMEMQ;
93: PetscLogEventEnd(SNES_FunctionEval,snes,X,F,0);
94: DMDAVecRestoreArray(dm,F,&f);
95: } break;
96: case ADD_VALUES: {
97: Vec Floc;
98: DMGetLocalVector(dm,&Floc);
99: VecZeroEntries(Floc);
100: DMDAVecGetArray(dm,Floc,&f);
101: PetscLogEventBegin(SNES_FunctionEval,snes,X,F,0);
102: CHKMEMQ;
103: (*dmdasnes->residuallocal)(&info,x,f,dmdasnes->residuallocalctx);
104: CHKMEMQ;
105: PetscLogEventEnd(SNES_FunctionEval,snes,X,F,0);
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: if (snes->domainerror) {
117: VecSetInf(F);
118: }
119: return(0);
120: }
124: static PetscErrorCode SNESComputeObjective_DMDA(SNES snes,Vec X,PetscReal *ob,void *ctx)
125: {
127: DM dm;
128: DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx;
129: DMDALocalInfo info;
130: Vec Xloc;
131: void *x;
137: if (!dmdasnes->objectivelocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
138: SNESGetDM(snes,&dm);
139: DMGetLocalVector(dm,&Xloc);
140: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
141: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
142: DMDAGetLocalInfo(dm,&info);
143: DMDAVecGetArray(dm,Xloc,&x);
144: CHKMEMQ;
145: (*dmdasnes->objectivelocal)(&info,x,ob,dmdasnes->objectivelocalctx);
146: CHKMEMQ;
147: DMDAVecRestoreArray(dm,Xloc,&x);
148: DMRestoreLocalVector(dm,&Xloc);
149: return(0);
150: }
155: PetscErrorCode SNESComputeJacobian_DMDA(SNES snes,Vec X,Mat A,Mat B,void *ctx)
156: {
158: DM dm;
159: DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx;
160: DMDALocalInfo info;
161: Vec Xloc;
162: void *x;
165: if (!dmdasnes->residuallocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
166: SNESGetDM(snes,&dm);
168: if (dmdasnes->jacobianlocal) {
169: DMGetLocalVector(dm,&Xloc);
170: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
171: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
172: DMDAGetLocalInfo(dm,&info);
173: DMDAVecGetArray(dm,Xloc,&x);
174: CHKMEMQ;
175: (*dmdasnes->jacobianlocal)(&info,x,A,B,dmdasnes->jacobianlocalctx);
176: CHKMEMQ;
177: DMDAVecRestoreArray(dm,Xloc,&x);
178: DMRestoreLocalVector(dm,&Xloc);
179: } else {
180: MatFDColoring fdcoloring;
181: PetscObjectQuery((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject*)&fdcoloring);
182: if (!fdcoloring) {
183: ISColoring coloring;
185: DMCreateColoring(dm,dm->coloringtype,&coloring);
186: MatFDColoringCreate(B,coloring,&fdcoloring);
187: switch (dm->coloringtype) {
188: case IS_COLORING_GLOBAL:
189: MatFDColoringSetFunction(fdcoloring,(PetscErrorCode (*)(void))SNESComputeFunction_DMDA,dmdasnes);
190: break;
191: default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"No support for coloring type '%s'",ISColoringTypes[dm->coloringtype]);
192: }
193: PetscObjectSetOptionsPrefix((PetscObject)fdcoloring,((PetscObject)dm)->prefix);
194: MatFDColoringSetFromOptions(fdcoloring);
195: MatFDColoringSetUp(B,coloring,fdcoloring);
196: ISColoringDestroy(&coloring);
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: MatFDColoringApply(B,fdcoloring,X,snes);
209: }
210: /* This will be redundant if the user called both, but it's too common to forget. */
211: if (A != B) {
212: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
213: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
214: }
215: return(0);
216: }
220: /*@C
221: DMDASNESSetFunctionLocal - set a local residual evaluation function
223: Logically Collective
225: Input Arguments:
226: + dm - DM to associate callback with
227: . imode - INSERT_VALUES if local function computes owned part, ADD_VALUES if it contributes to ghosted part
228: . func - local residual evaluation
229: - ctx - optional context for local residual evaluation
231: Calling sequence:
232: For PetscErrorCode (*func)(DMDALocalInfo *info,void *x, void *f, void *ctx),
233: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
234: . x - dimensional pointer to state at which to evaluate residual (e.g. PetscScalar *x or **x or ***x)
235: . f - dimensional pointer to residual, write the residual here (e.g. PetscScalar *f or **f or ***f)
236: - ctx - optional context passed above
238: Level: beginner
240: .seealso: DMDASNESSetJacobianLocal(), DMSNESSetFunction(), 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 Jacobian evaluation
274: - ctx - optional context for local Jacobian evaluation
276: Calling sequence:
277: For PetscErrorCode (*func)(DMDALocalInfo *info,void *x,Mat J,Mat M,void *ctx),
278: + info - DMDALocalInfo defining the subdomain to evaluate the Jacobian at
279: . x - dimensional pointer to state at which to evaluate Jacobian (e.g. PetscScalar *x or **x or ***x)
280: . J - Mat object for the Jacobian
281: . M - Mat object for the Jacobian preconditioner matrix
282: - ctx - optional context passed above
284: Level: beginner
286: .seealso: DMDASNESSetFunctionLocal(), DMSNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
287: @*/
288: PetscErrorCode DMDASNESSetJacobianLocal(DM dm,PetscErrorCode (*func)(DMDALocalInfo*,void*,Mat,Mat,void*),void *ctx)
289: {
291: DMSNES sdm;
292: DMSNES_DA *dmdasnes;
296: DMGetDMSNESWrite(dm,&sdm);
297: DMDASNESGetContext(dm,sdm,&dmdasnes);
299: dmdasnes->jacobianlocal = func;
300: dmdasnes->jacobianlocalctx = ctx;
302: DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes);
303: return(0);
304: }
309: /*@C
310: DMDASNESSetObjectiveLocal - set a local residual evaluation function
312: Logically Collective
314: Input Arguments:
315: + dm - DM to associate callback with
316: . func - local objective evaluation
317: - ctx - optional context for local residual evaluation
319: Calling sequence for func:
320: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
321: . x - dimensional pointer to state at which to evaluate residual
322: . ob - eventual objective value
323: - ctx - optional context passed above
325: Level: beginner
327: .seealso: DMSNESSetFunction(), DMDASNESSetJacobianLocal(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
328: @*/
329: PetscErrorCode DMDASNESSetObjectiveLocal(DM dm,DMDASNESObjective func,void *ctx)
330: {
332: DMSNES sdm;
333: DMSNES_DA *dmdasnes;
337: DMGetDMSNESWrite(dm,&sdm);
338: DMDASNESGetContext(dm,sdm,&dmdasnes);
340: dmdasnes->objectivelocal = func;
341: dmdasnes->objectivelocalctx = ctx;
343: DMSNESSetObjective(dm,SNESComputeObjective_DMDA,dmdasnes);
344: return(0);
345: }
349: static PetscErrorCode SNESComputePicard_DMDA(SNES snes,Vec X,Vec F,void *ctx)
350: {
352: DM dm;
353: DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx;
354: DMDALocalInfo info;
355: Vec Xloc;
356: void *x,*f;
362: if (!dmdasnes->rhsplocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
363: SNESGetDM(snes,&dm);
364: DMGetLocalVector(dm,&Xloc);
365: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
366: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
367: DMDAGetLocalInfo(dm,&info);
368: DMDAVecGetArray(dm,Xloc,&x);
369: switch (dmdasnes->residuallocalimode) {
370: case INSERT_VALUES: {
371: DMDAVecGetArray(dm,F,&f);
372: CHKMEMQ;
373: (*dmdasnes->rhsplocal)(&info,x,f,dmdasnes->picardlocalctx);
374: CHKMEMQ;
375: DMDAVecRestoreArray(dm,F,&f);
376: } break;
377: case ADD_VALUES: {
378: Vec Floc;
379: DMGetLocalVector(dm,&Floc);
380: VecZeroEntries(Floc);
381: DMDAVecGetArray(dm,Floc,&f);
382: CHKMEMQ;
383: (*dmdasnes->rhsplocal)(&info,x,f,dmdasnes->picardlocalctx);
384: CHKMEMQ;
385: DMDAVecRestoreArray(dm,Floc,&f);
386: VecZeroEntries(F);
387: DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
388: DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
389: DMRestoreLocalVector(dm,&Floc);
390: } break;
391: default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdasnes->residuallocalimode);
392: }
393: DMDAVecRestoreArray(dm,Xloc,&x);
394: DMRestoreLocalVector(dm,&Xloc);
395: return(0);
396: }
400: static PetscErrorCode SNESComputePicardJacobian_DMDA(SNES snes,Vec X,Mat A,Mat B,void *ctx)
401: {
403: DM dm;
404: DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx;
405: DMDALocalInfo info;
406: Vec Xloc;
407: void *x;
410: if (!dmdasnes->jacobianplocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
411: SNESGetDM(snes,&dm);
413: DMGetLocalVector(dm,&Xloc);
414: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
415: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
416: DMDAGetLocalInfo(dm,&info);
417: DMDAVecGetArray(dm,Xloc,&x);
418: CHKMEMQ;
419: (*dmdasnes->jacobianplocal)(&info,x,A,B,dmdasnes->picardlocalctx);
420: CHKMEMQ;
421: DMDAVecRestoreArray(dm,Xloc,&x);
422: DMRestoreLocalVector(dm,&Xloc);
423: if (A != B) {
424: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
425: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
426: }
427: return(0);
428: }
432: /*@C
433: DMDASNESSetPicardLocal - set a local right hand side and matrix evaluation function for Picard iteration
435: Logically Collective
437: Input Arguments:
438: + dm - DM to associate callback with
439: . imode - INSERT_VALUES if local function computes owned part, ADD_VALUES if it contributes to ghosted part
440: . func - local residual evaluation
441: - ctx - optional context for local residual evaluation
443: Calling sequence for func:
444: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
445: . x - dimensional pointer to state at which to evaluate residual
446: . f - dimensional pointer to residual, write the residual here
447: - ctx - optional context passed above
449: Notes: The user must use
450: extern PetscErrorCode SNESPicardComputeFunction(SNES,Vec,Vec,void*);
451: extern PetscErrorCode SNESPicardComputeJacobian(SNES,Vec,Mat,Mat,MatStructure*,void*);
452: SNESSetFunction(snes,NULL,SNESPicardComputeFunction,&user);
453: in their code before calling this routine.
456: Level: beginner
458: .seealso: DMSNESSetFunction(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
459: @*/
460: PetscErrorCode DMDASNESSetPicardLocal(DM dm,InsertMode imode,PetscErrorCode (*func)(DMDALocalInfo*,void*,void*,void*),
461: PetscErrorCode (*jac)(DMDALocalInfo*,void*,Mat,Mat,void*),void *ctx)
462: {
464: DMSNES sdm;
465: DMSNES_DA *dmdasnes;
469: DMGetDMSNESWrite(dm,&sdm);
470: DMDASNESGetContext(dm,sdm,&dmdasnes);
472: dmdasnes->residuallocalimode = imode;
473: dmdasnes->rhsplocal = func;
474: dmdasnes->jacobianplocal = jac;
475: dmdasnes->picardlocalctx = ctx;
477: DMSNESSetPicard(dm,SNESComputePicard_DMDA,SNESComputePicardJacobian_DMDA,dmdasnes);
478: return(0);
479: }