Actual source code: dmdasnes.c
petsc-3.3-p7 2013-05-11
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: void *residuallocalctx;
10: void *jacobianlocalctx;
11: InsertMode residuallocalimode;
12: } DM_DA_SNES;
16: static PetscErrorCode SNESDMDestroy_DMDA(SNESDM sdm)
17: {
21: PetscFree(sdm->data);
22: return(0);
23: }
27: static PetscErrorCode DMDASNESGetContext(DM dm,SNESDM sdm,DM_DA_SNES **dmdasnes)
28: {
32: *dmdasnes = PETSC_NULL;
33: if (!sdm->data) {
34: PetscNewLog(dm,DM_DA_SNES,&sdm->data);
35: sdm->destroy = SNESDMDestroy_DMDA;
36: }
37: *dmdasnes = (DM_DA_SNES*)sdm->data;
38: return(0);
39: }
43: /*
44: This function should eventually replace:
45: DMDAComputeFunction() and DMDAComputeFunction1()
46: */
47: static PetscErrorCode SNESComputeFunction_DMDA(SNES snes,Vec X,Vec F,void *ctx)
48: {
50: DM dm;
51: DM_DA_SNES *dmdasnes = (DM_DA_SNES*)ctx;
52: DMDALocalInfo info;
53: Vec Xloc;
54: void *x,*f;
60: if (!dmdasnes->residuallocal) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_PLIB,"Corrupt context");
61: SNESGetDM(snes,&dm);
62: DMGetLocalVector(dm,&Xloc);
63: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
64: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
65: DMDAGetLocalInfo(dm,&info);
66: DMDAVecGetArray(dm,Xloc,&x);
67: switch (dmdasnes->residuallocalimode) {
68: case INSERT_VALUES: {
69: DMDAVecGetArray(dm,F,&f);
70: CHKMEMQ;
71: (*dmdasnes->residuallocal)(&info,x,f,dmdasnes->residuallocalctx);
72: CHKMEMQ;
73: DMDAVecRestoreArray(dm,F,&f);
74: } break;
75: case ADD_VALUES: {
76: Vec Floc;
77: DMGetLocalVector(dm,&Floc);
78: VecZeroEntries(Floc);
79: DMDAVecGetArray(dm,Floc,&f);
80: CHKMEMQ;
81: (*dmdasnes->residuallocal)(&info,x,f,dmdasnes->residuallocalctx);
82: CHKMEMQ;
83: DMDAVecRestoreArray(dm,Floc,&f);
84: VecZeroEntries(F);
85: DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
86: DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
87: DMRestoreLocalVector(dm,&Floc);
88: } break;
89: default: SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdasnes->residuallocalimode);
90: }
91: DMDAVecRestoreArray(dm,Xloc,&x);
92: DMRestoreLocalVector(dm,&Xloc);
93: return(0);
94: }
98: /*
99: This function should eventually replace:
100: DMComputeJacobian() and DMDAComputeJacobian1()
101: */
102: static PetscErrorCode SNESComputeJacobian_DMDA(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *mstr,void *ctx)
103: {
105: DM dm;
106: DM_DA_SNES *dmdasnes = (DM_DA_SNES*)ctx;
107: DMDALocalInfo info;
108: Vec Xloc;
109: void *x;
112: if (!dmdasnes->residuallocal) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_PLIB,"Corrupt context");
113: SNESGetDM(snes,&dm);
115: if (dmdasnes->jacobianlocal) {
116: DMGetLocalVector(dm,&Xloc);
117: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
118: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
119: DMDAGetLocalInfo(dm,&info);
120: DMDAVecGetArray(dm,Xloc,&x);
121: CHKMEMQ;
122: (*dmdasnes->jacobianlocal)(&info,x,*A,*B,mstr,dmdasnes->jacobianlocalctx);
123: CHKMEMQ;
124: DMDAVecRestoreArray(dm,Xloc,&x);
125: DMRestoreLocalVector(dm,&Xloc);
126: } else {
127: MatFDColoring fdcoloring;
128: PetscObjectQuery((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject*)&fdcoloring);
129: if (!fdcoloring) {
130: ISColoring coloring;
132: DMCreateColoring(dm,dm->coloringtype,MATAIJ,&coloring);
133: MatFDColoringCreate(*B,coloring,&fdcoloring);
134: ISColoringDestroy(&coloring);
135: switch (dm->coloringtype) {
136: case IS_COLORING_GLOBAL:
137: MatFDColoringSetFunction(fdcoloring,(PetscErrorCode(*)(void))SNESComputeFunction_DMDA,dmdasnes);
138: break;
139: default: SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_SUP,"No support for coloring type '%s'",ISColoringTypes[dm->coloringtype]);
140: }
141: PetscObjectSetOptionsPrefix((PetscObject)fdcoloring,((PetscObject)dm)->prefix);
142: MatFDColoringSetFromOptions(fdcoloring);
143: PetscObjectCompose((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject)fdcoloring);
144: PetscObjectDereference((PetscObject)fdcoloring);
146: /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call
147: * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the
148: * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually
149: * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be
150: * taken when the PetscOList for the Vec inside MatFDColoring is destroyed.
151: */
152: PetscObjectDereference((PetscObject)dm);
153: }
154: *mstr = SAME_NONZERO_PATTERN;
155: MatFDColoringApply(*B,fdcoloring,X,mstr,snes);
156: }
157: /* This will be redundant if the user called both, but it's too common to forget. */
158: if (*A != *B) {
159: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
160: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
161: }
162: return(0);
163: }
167: /*@C
168: DMDASNESSetFunctionLocal - set a local residual evaluation function
170: Logically Collective
172: Input Arguments:
173: + dm - DM to associate callback with
174: . func - local residual evaluation
175: - ctx - optional context for local residual evaluation
177: Calling sequence for func:
178: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
179: . imode - INSERT_VALUES if local function computes owned part, ADD_VALUES if it contributes to ghosted part
180: . x - dimensional pointer to state at which to evaluate residual
181: . f - dimensional pointer to residual, write the residual here
182: - ctx - optional context passed above
184: Level: beginner
186: .seealso: DMSNESSetFunction(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
187: @*/
188: PetscErrorCode DMDASNESSetFunctionLocal(DM dm,InsertMode imode,PetscErrorCode (*func)(DMDALocalInfo*,void*,void*,void*),void *ctx)
189: {
191: SNESDM sdm;
192: DM_DA_SNES *dmdasnes;
196: DMSNESGetContext(dm,&sdm);
197: DMDASNESGetContext(dm,sdm,&dmdasnes);
198: dmdasnes->residuallocalimode = imode;
199: dmdasnes->residuallocal = func;
200: dmdasnes->residuallocalctx = ctx;
201: DMSNESSetFunction(dm,SNESComputeFunction_DMDA,dmdasnes);
202: if (!sdm->computejacobian) { /* Call us for the Jacobian too, can be overridden by the user. */
203: DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes);
204: }
205: return(0);
206: }
210: /*@C
211: DMDASNESSetJacobianLocal - set a local residual evaluation function
213: Logically Collective
215: Input Arguments:
216: + dm - DM to associate callback with
217: . func - local residual evaluation
218: - ctx - optional context for local residual evaluation
220: Calling sequence for func:
221: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
222: . x - dimensional pointer to state at which to evaluate residual
223: . f - dimensional pointer to residual, write the residual here
224: - ctx - optional context passed above
226: Level: beginner
228: .seealso: DMSNESSetJacobian(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
229: @*/
230: PetscErrorCode DMDASNESSetJacobianLocal(DM dm,PetscErrorCode (*func)(DMDALocalInfo*,void*,Mat,Mat,MatStructure*,void*),void *ctx)
231: {
233: SNESDM sdm;
234: DM_DA_SNES *dmdasnes;
238: DMSNESGetContext(dm,&sdm);
239: DMDASNESGetContext(dm,sdm,&dmdasnes);
240: dmdasnes->jacobianlocal = func;
241: dmdasnes->jacobianlocalctx = ctx;
242: DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes);
243: return(0);
244: }