Actual source code: dmdasnes.c
petsc-3.7.5 2017-01-01
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: /* Routine is called by example, hence must be labeled PETSC_EXTERN */
156: PETSC_EXTERN PetscErrorCode SNESComputeJacobian_DMDA(SNES snes,Vec X,Mat A,Mat B,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,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,&coloring);
187: MatFDColoringCreate(B,coloring,&fdcoloring);
188: switch (dm->coloringtype) {
189: case IS_COLORING_GLOBAL:
190: MatFDColoringSetFunction(fdcoloring,(PetscErrorCode (*)(void))SNESComputeFunction_DMDA,dmdasnes);
191: break;
192: default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"No support for coloring type '%s'",ISColoringTypes[dm->coloringtype]);
193: }
194: PetscObjectSetOptionsPrefix((PetscObject)fdcoloring,((PetscObject)dm)->prefix);
195: MatFDColoringSetFromOptions(fdcoloring);
196: MatFDColoringSetUp(B,coloring,fdcoloring);
197: ISColoringDestroy(&coloring);
198: PetscObjectCompose((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject)fdcoloring);
199: PetscObjectDereference((PetscObject)fdcoloring);
201: /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call
202: * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the
203: * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually
204: * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be
205: * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed.
206: */
207: PetscObjectDereference((PetscObject)dm);
208: }
209: MatFDColoringApply(B,fdcoloring,X,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:
233: For PetscErrorCode (*func)(DMDALocalInfo *info,void *x, void *f, void *ctx),
234: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
235: . x - dimensional pointer to state at which to evaluate residual (e.g. PetscScalar *x or **x or ***x)
236: . f - dimensional pointer to residual, write the residual here (e.g. PetscScalar *f or **f or ***f)
237: - ctx - optional context passed above
239: Level: beginner
241: .seealso: DMDASNESSetJacobianLocal(), DMSNESSetFunction(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
242: @*/
243: PetscErrorCode DMDASNESSetFunctionLocal(DM dm,InsertMode imode,PetscErrorCode (*func)(DMDALocalInfo*,void*,void*,void*),void *ctx)
244: {
246: DMSNES sdm;
247: DMSNES_DA *dmdasnes;
251: DMGetDMSNESWrite(dm,&sdm);
252: DMDASNESGetContext(dm,sdm,&dmdasnes);
254: dmdasnes->residuallocalimode = imode;
255: dmdasnes->residuallocal = func;
256: dmdasnes->residuallocalctx = ctx;
258: DMSNESSetFunction(dm,SNESComputeFunction_DMDA,dmdasnes);
259: if (!sdm->ops->computejacobian) { /* Call us for the Jacobian too, can be overridden by the user. */
260: DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes);
261: }
262: return(0);
263: }
267: /*@C
268: DMDASNESSetJacobianLocal - set a local Jacobian evaluation function
270: Logically Collective
272: Input Arguments:
273: + dm - DM to associate callback with
274: . func - local Jacobian evaluation
275: - ctx - optional context for local Jacobian evaluation
277: Calling sequence:
278: For PetscErrorCode (*func)(DMDALocalInfo *info,void *x,Mat J,Mat M,void *ctx),
279: + info - DMDALocalInfo defining the subdomain to evaluate the Jacobian at
280: . x - dimensional pointer to state at which to evaluate Jacobian (e.g. PetscScalar *x or **x or ***x)
281: . J - Mat object for the Jacobian
282: . M - Mat object for the Jacobian preconditioner matrix
283: - ctx - optional context passed above
285: Level: beginner
287: .seealso: DMDASNESSetFunctionLocal(), DMSNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
288: @*/
289: PetscErrorCode DMDASNESSetJacobianLocal(DM dm,PetscErrorCode (*func)(DMDALocalInfo*,void*,Mat,Mat,void*),void *ctx)
290: {
292: DMSNES sdm;
293: DMSNES_DA *dmdasnes;
297: DMGetDMSNESWrite(dm,&sdm);
298: DMDASNESGetContext(dm,sdm,&dmdasnes);
300: dmdasnes->jacobianlocal = func;
301: dmdasnes->jacobianlocalctx = ctx;
303: DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes);
304: return(0);
305: }
310: /*@C
311: DMDASNESSetObjectiveLocal - set a local residual evaluation function
313: Logically Collective
315: Input Arguments:
316: + dm - DM to associate callback with
317: . func - local objective evaluation
318: - ctx - optional context for local residual evaluation
320: Calling sequence for func:
321: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
322: . x - dimensional pointer to state at which to evaluate residual
323: . ob - eventual objective value
324: - ctx - optional context passed above
326: Level: beginner
328: .seealso: DMSNESSetFunction(), DMDASNESSetJacobianLocal(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
329: @*/
330: PetscErrorCode DMDASNESSetObjectiveLocal(DM dm,DMDASNESObjective func,void *ctx)
331: {
333: DMSNES sdm;
334: DMSNES_DA *dmdasnes;
338: DMGetDMSNESWrite(dm,&sdm);
339: DMDASNESGetContext(dm,sdm,&dmdasnes);
341: dmdasnes->objectivelocal = func;
342: dmdasnes->objectivelocalctx = ctx;
344: DMSNESSetObjective(dm,SNESComputeObjective_DMDA,dmdasnes);
345: return(0);
346: }
350: static PetscErrorCode SNESComputePicard_DMDA(SNES snes,Vec X,Vec F,void *ctx)
351: {
353: DM dm;
354: DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx;
355: DMDALocalInfo info;
356: Vec Xloc;
357: void *x,*f;
363: if (!dmdasnes->rhsplocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
364: SNESGetDM(snes,&dm);
365: DMGetLocalVector(dm,&Xloc);
366: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
367: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
368: DMDAGetLocalInfo(dm,&info);
369: DMDAVecGetArray(dm,Xloc,&x);
370: switch (dmdasnes->residuallocalimode) {
371: case INSERT_VALUES: {
372: DMDAVecGetArray(dm,F,&f);
373: CHKMEMQ;
374: (*dmdasnes->rhsplocal)(&info,x,f,dmdasnes->picardlocalctx);
375: CHKMEMQ;
376: DMDAVecRestoreArray(dm,F,&f);
377: } break;
378: case ADD_VALUES: {
379: Vec Floc;
380: DMGetLocalVector(dm,&Floc);
381: VecZeroEntries(Floc);
382: DMDAVecGetArray(dm,Floc,&f);
383: CHKMEMQ;
384: (*dmdasnes->rhsplocal)(&info,x,f,dmdasnes->picardlocalctx);
385: CHKMEMQ;
386: DMDAVecRestoreArray(dm,Floc,&f);
387: VecZeroEntries(F);
388: DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
389: DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
390: DMRestoreLocalVector(dm,&Floc);
391: } break;
392: default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdasnes->residuallocalimode);
393: }
394: DMDAVecRestoreArray(dm,Xloc,&x);
395: DMRestoreLocalVector(dm,&Xloc);
396: return(0);
397: }
401: static PetscErrorCode SNESComputePicardJacobian_DMDA(SNES snes,Vec X,Mat A,Mat B,void *ctx)
402: {
404: DM dm;
405: DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx;
406: DMDALocalInfo info;
407: Vec Xloc;
408: void *x;
411: if (!dmdasnes->jacobianplocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
412: SNESGetDM(snes,&dm);
414: DMGetLocalVector(dm,&Xloc);
415: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
416: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
417: DMDAGetLocalInfo(dm,&info);
418: DMDAVecGetArray(dm,Xloc,&x);
419: CHKMEMQ;
420: (*dmdasnes->jacobianplocal)(&info,x,A,B,dmdasnes->picardlocalctx);
421: CHKMEMQ;
422: DMDAVecRestoreArray(dm,Xloc,&x);
423: DMRestoreLocalVector(dm,&Xloc);
424: if (A != B) {
425: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
426: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
427: }
428: return(0);
429: }
433: /*@C
434: DMDASNESSetPicardLocal - set a local right hand side and matrix evaluation function for Picard iteration
436: Logically Collective
438: Input Arguments:
439: + dm - DM to associate callback with
440: . imode - INSERT_VALUES if local function computes owned part, ADD_VALUES if it contributes to ghosted part
441: . func - local residual evaluation
442: - ctx - optional context for local residual evaluation
444: Calling sequence for func:
445: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
446: . x - dimensional pointer to state at which to evaluate residual
447: . f - dimensional pointer to residual, write the residual here
448: - ctx - optional context passed above
450: Notes: The user must use
451: extern PetscErrorCode SNESPicardComputeFunction(SNES,Vec,Vec,void*);
452: extern PetscErrorCode SNESPicardComputeJacobian(SNES,Vec,Mat,Mat,MatStructure*,void*);
453: SNESSetFunction(snes,NULL,SNESPicardComputeFunction,&user);
454: in their code before calling this routine.
457: Level: beginner
459: .seealso: DMSNESSetFunction(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
460: @*/
461: PetscErrorCode DMDASNESSetPicardLocal(DM dm,InsertMode imode,PetscErrorCode (*func)(DMDALocalInfo*,void*,void*,void*),
462: PetscErrorCode (*jac)(DMDALocalInfo*,void*,Mat,Mat,void*),void *ctx)
463: {
465: DMSNES sdm;
466: DMSNES_DA *dmdasnes;
470: DMGetDMSNESWrite(dm,&sdm);
471: DMDASNESGetContext(dm,sdm,&dmdasnes);
473: dmdasnes->residuallocalimode = imode;
474: dmdasnes->rhsplocal = func;
475: dmdasnes->jacobianplocal = jac;
476: dmdasnes->picardlocalctx = ctx;
478: DMSNESSetPicard(dm,SNESComputePicard_DMDA,SNESComputePicardJacobian_DMDA,dmdasnes);
479: return(0);
480: }