Actual source code: dmdasnes.c
petsc-3.6.1 2015-08-06
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: if (snes->domainerror) {
113: VecSetInf(F);
114: }
115: return(0);
116: }
120: static PetscErrorCode SNESComputeObjective_DMDA(SNES snes,Vec X,PetscReal *ob,void *ctx)
121: {
123: DM dm;
124: DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx;
125: DMDALocalInfo info;
126: Vec Xloc;
127: void *x;
133: if (!dmdasnes->objectivelocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
134: SNESGetDM(snes,&dm);
135: DMGetLocalVector(dm,&Xloc);
136: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
137: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
138: DMDAGetLocalInfo(dm,&info);
139: DMDAVecGetArray(dm,Xloc,&x);
140: CHKMEMQ;
141: (*dmdasnes->objectivelocal)(&info,x,ob,dmdasnes->objectivelocalctx);
142: CHKMEMQ;
143: DMDAVecRestoreArray(dm,Xloc,&x);
144: DMRestoreLocalVector(dm,&Xloc);
145: return(0);
146: }
151: PetscErrorCode SNESComputeJacobian_DMDA(SNES snes,Vec X,Mat A,Mat B,void *ctx)
152: {
154: DM dm;
155: DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx;
156: DMDALocalInfo info;
157: Vec Xloc;
158: void *x;
161: if (!dmdasnes->residuallocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
162: SNESGetDM(snes,&dm);
164: if (dmdasnes->jacobianlocal) {
165: DMGetLocalVector(dm,&Xloc);
166: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
167: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
168: DMDAGetLocalInfo(dm,&info);
169: DMDAVecGetArray(dm,Xloc,&x);
170: CHKMEMQ;
171: (*dmdasnes->jacobianlocal)(&info,x,A,B,dmdasnes->jacobianlocalctx);
172: CHKMEMQ;
173: DMDAVecRestoreArray(dm,Xloc,&x);
174: DMRestoreLocalVector(dm,&Xloc);
175: } else {
176: MatFDColoring fdcoloring;
177: PetscObjectQuery((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject*)&fdcoloring);
178: if (!fdcoloring) {
179: ISColoring coloring;
181: DMCreateColoring(dm,dm->coloringtype,&coloring);
182: MatFDColoringCreate(B,coloring,&fdcoloring);
183: switch (dm->coloringtype) {
184: case IS_COLORING_GLOBAL:
185: MatFDColoringSetFunction(fdcoloring,(PetscErrorCode (*)(void))SNESComputeFunction_DMDA,dmdasnes);
186: break;
187: default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"No support for coloring type '%s'",ISColoringTypes[dm->coloringtype]);
188: }
189: PetscObjectSetOptionsPrefix((PetscObject)fdcoloring,((PetscObject)dm)->prefix);
190: MatFDColoringSetFromOptions(fdcoloring);
191: MatFDColoringSetUp(B,coloring,fdcoloring);
192: ISColoringDestroy(&coloring);
193: PetscObjectCompose((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject)fdcoloring);
194: PetscObjectDereference((PetscObject)fdcoloring);
196: /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call
197: * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the
198: * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually
199: * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be
200: * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed.
201: */
202: PetscObjectDereference((PetscObject)dm);
203: }
204: MatFDColoringApply(B,fdcoloring,X,snes);
205: }
206: /* This will be redundant if the user called both, but it's too common to forget. */
207: if (A != B) {
208: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
209: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
210: }
211: return(0);
212: }
216: /*@C
217: DMDASNESSetFunctionLocal - set a local residual evaluation function
219: Logically Collective
221: Input Arguments:
222: + dm - DM to associate callback with
223: . imode - INSERT_VALUES if local function computes owned part, ADD_VALUES if it contributes to ghosted part
224: . func - local residual evaluation
225: - ctx - optional context for local residual evaluation
227: Calling sequence:
228: For PetscErrorCode (*func)(DMDALocalInfo *info,void *x, void *f, void *ctx),
229: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
230: . x - dimensional pointer to state at which to evaluate residual (e.g. PetscScalar *x or **x or ***x)
231: . f - dimensional pointer to residual, write the residual here (e.g. PetscScalar *f or **f or ***f)
232: - ctx - optional context passed above
234: Level: beginner
236: .seealso: DMDASNESSetJacobianLocal(), DMSNESSetFunction(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
237: @*/
238: PetscErrorCode DMDASNESSetFunctionLocal(DM dm,InsertMode imode,PetscErrorCode (*func)(DMDALocalInfo*,void*,void*,void*),void *ctx)
239: {
241: DMSNES sdm;
242: DMSNES_DA *dmdasnes;
246: DMGetDMSNESWrite(dm,&sdm);
247: DMDASNESGetContext(dm,sdm,&dmdasnes);
249: dmdasnes->residuallocalimode = imode;
250: dmdasnes->residuallocal = func;
251: dmdasnes->residuallocalctx = ctx;
253: DMSNESSetFunction(dm,SNESComputeFunction_DMDA,dmdasnes);
254: if (!sdm->ops->computejacobian) { /* Call us for the Jacobian too, can be overridden by the user. */
255: DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes);
256: }
257: return(0);
258: }
262: /*@C
263: DMDASNESSetJacobianLocal - set a local Jacobian evaluation function
265: Logically Collective
267: Input Arguments:
268: + dm - DM to associate callback with
269: . func - local Jacobian evaluation
270: - ctx - optional context for local Jacobian evaluation
272: Calling sequence:
273: For PetscErrorCode (*func)(DMDALocalInfo *info,void *x,Mat J,Mat M,void *ctx),
274: + info - DMDALocalInfo defining the subdomain to evaluate the Jacobian at
275: . x - dimensional pointer to state at which to evaluate Jacobian (e.g. PetscScalar *x or **x or ***x)
276: . J - Mat object for the Jacobian
277: . M - Mat object for the Jacobian preconditioner matrix
278: - ctx - optional context passed above
280: Level: beginner
282: .seealso: DMDASNESSetFunctionLocal(), DMSNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
283: @*/
284: PetscErrorCode DMDASNESSetJacobianLocal(DM dm,PetscErrorCode (*func)(DMDALocalInfo*,void*,Mat,Mat,void*),void *ctx)
285: {
287: DMSNES sdm;
288: DMSNES_DA *dmdasnes;
292: DMGetDMSNESWrite(dm,&sdm);
293: DMDASNESGetContext(dm,sdm,&dmdasnes);
295: dmdasnes->jacobianlocal = func;
296: dmdasnes->jacobianlocalctx = ctx;
298: DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes);
299: return(0);
300: }
305: /*@C
306: DMDASNESSetObjectiveLocal - set a local residual evaluation function
308: Logically Collective
310: Input Arguments:
311: + dm - DM to associate callback with
312: . func - local objective evaluation
313: - ctx - optional context for local residual evaluation
315: Calling sequence for func:
316: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
317: . x - dimensional pointer to state at which to evaluate residual
318: . ob - eventual objective value
319: - ctx - optional context passed above
321: Level: beginner
323: .seealso: DMSNESSetFunction(), DMDASNESSetJacobianLocal(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
324: @*/
325: PetscErrorCode DMDASNESSetObjectiveLocal(DM dm,DMDASNESObjective func,void *ctx)
326: {
328: DMSNES sdm;
329: DMSNES_DA *dmdasnes;
333: DMGetDMSNESWrite(dm,&sdm);
334: DMDASNESGetContext(dm,sdm,&dmdasnes);
336: dmdasnes->objectivelocal = func;
337: dmdasnes->objectivelocalctx = ctx;
339: DMSNESSetObjective(dm,SNESComputeObjective_DMDA,dmdasnes);
340: return(0);
341: }
345: static PetscErrorCode SNESComputePicard_DMDA(SNES snes,Vec X,Vec F,void *ctx)
346: {
348: DM dm;
349: DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx;
350: DMDALocalInfo info;
351: Vec Xloc;
352: void *x,*f;
358: if (!dmdasnes->rhsplocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
359: SNESGetDM(snes,&dm);
360: DMGetLocalVector(dm,&Xloc);
361: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
362: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
363: DMDAGetLocalInfo(dm,&info);
364: DMDAVecGetArray(dm,Xloc,&x);
365: switch (dmdasnes->residuallocalimode) {
366: case INSERT_VALUES: {
367: DMDAVecGetArray(dm,F,&f);
368: CHKMEMQ;
369: (*dmdasnes->rhsplocal)(&info,x,f,dmdasnes->picardlocalctx);
370: CHKMEMQ;
371: DMDAVecRestoreArray(dm,F,&f);
372: } break;
373: case ADD_VALUES: {
374: Vec Floc;
375: DMGetLocalVector(dm,&Floc);
376: VecZeroEntries(Floc);
377: DMDAVecGetArray(dm,Floc,&f);
378: CHKMEMQ;
379: (*dmdasnes->rhsplocal)(&info,x,f,dmdasnes->picardlocalctx);
380: CHKMEMQ;
381: DMDAVecRestoreArray(dm,Floc,&f);
382: VecZeroEntries(F);
383: DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
384: DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
385: DMRestoreLocalVector(dm,&Floc);
386: } break;
387: default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdasnes->residuallocalimode);
388: }
389: DMDAVecRestoreArray(dm,Xloc,&x);
390: DMRestoreLocalVector(dm,&Xloc);
391: return(0);
392: }
396: static PetscErrorCode SNESComputePicardJacobian_DMDA(SNES snes,Vec X,Mat A,Mat B,void *ctx)
397: {
399: DM dm;
400: DMSNES_DA *dmdasnes = (DMSNES_DA*)ctx;
401: DMDALocalInfo info;
402: Vec Xloc;
403: void *x;
406: if (!dmdasnes->jacobianplocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
407: SNESGetDM(snes,&dm);
409: DMGetLocalVector(dm,&Xloc);
410: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
411: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
412: DMDAGetLocalInfo(dm,&info);
413: DMDAVecGetArray(dm,Xloc,&x);
414: CHKMEMQ;
415: (*dmdasnes->jacobianplocal)(&info,x,A,B,dmdasnes->picardlocalctx);
416: CHKMEMQ;
417: DMDAVecRestoreArray(dm,Xloc,&x);
418: DMRestoreLocalVector(dm,&Xloc);
419: if (A != B) {
420: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
421: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
422: }
423: return(0);
424: }
428: /*@C
429: DMDASNESSetPicardLocal - set a local right hand side and matrix evaluation function for Picard iteration
431: Logically Collective
433: Input Arguments:
434: + dm - DM to associate callback with
435: . imode - INSERT_VALUES if local function computes owned part, ADD_VALUES if it contributes to ghosted part
436: . func - local residual evaluation
437: - ctx - optional context for local residual evaluation
439: Calling sequence for func:
440: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
441: . x - dimensional pointer to state at which to evaluate residual
442: . f - dimensional pointer to residual, write the residual here
443: - ctx - optional context passed above
445: Notes: The user must use
446: extern PetscErrorCode SNESPicardComputeFunction(SNES,Vec,Vec,void*);
447: extern PetscErrorCode SNESPicardComputeJacobian(SNES,Vec,Mat,Mat,MatStructure*,void*);
448: SNESSetFunction(snes,NULL,SNESPicardComputeFunction,&user);
449: in their code before calling this routine.
452: Level: beginner
454: .seealso: DMSNESSetFunction(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
455: @*/
456: PetscErrorCode DMDASNESSetPicardLocal(DM dm,InsertMode imode,PetscErrorCode (*func)(DMDALocalInfo*,void*,void*,void*),
457: PetscErrorCode (*jac)(DMDALocalInfo*,void*,Mat,Mat,void*),void *ctx)
458: {
460: DMSNES sdm;
461: DMSNES_DA *dmdasnes;
465: DMGetDMSNESWrite(dm,&sdm);
466: DMDASNESGetContext(dm,sdm,&dmdasnes);
468: dmdasnes->residuallocalimode = imode;
469: dmdasnes->rhsplocal = func;
470: dmdasnes->jacobianplocal = jac;
471: dmdasnes->picardlocalctx = ctx;
473: DMSNESSetPicard(dm,SNESComputePicard_DMDA,SNESComputePicardJacobian_DMDA,dmdasnes);
474: return(0);
475: }