Actual source code: dmlocalsnes.c
1: #include <petsc/private/dmimpl.h>
2: #include <petsc/private/snesimpl.h>
4: typedef struct {
5: PetscErrorCode (*residuallocal)(DM,Vec,Vec,void*);
6: PetscErrorCode (*jacobianlocal)(DM,Vec,Mat,Mat,void*);
7: PetscErrorCode (*boundarylocal)(DM,Vec,void*);
8: void *residuallocalctx;
9: void *jacobianlocalctx;
10: void *boundarylocalctx;
11: } DMSNES_Local;
13: static PetscErrorCode DMSNESDestroy_DMLocal(DMSNES sdm)
14: {
15: PetscFree(sdm->data);
16: sdm->data = NULL;
17: return 0;
18: }
20: static PetscErrorCode DMSNESDuplicate_DMLocal(DMSNES oldsdm,DMSNES sdm)
21: {
22: if (sdm->data != oldsdm->data) {
23: PetscFree(sdm->data);
24: PetscNewLog(sdm,(DMSNES_Local**)&sdm->data);
25: if (oldsdm->data) PetscMemcpy(sdm->data,oldsdm->data,sizeof(DMSNES_Local));
26: }
27: return 0;
28: }
30: static PetscErrorCode DMLocalSNESGetContext(DM dm,DMSNES sdm,DMSNES_Local **dmlocalsnes)
31: {
32: *dmlocalsnes = NULL;
33: if (!sdm->data) {
34: PetscNewLog(dm,(DMSNES_Local**)&sdm->data);
36: sdm->ops->destroy = DMSNESDestroy_DMLocal;
37: sdm->ops->duplicate = DMSNESDuplicate_DMLocal;
38: }
39: *dmlocalsnes = (DMSNES_Local*)sdm->data;
40: return 0;
41: }
43: static PetscErrorCode SNESComputeFunction_DMLocal(SNES snes,Vec X,Vec F,void *ctx)
44: {
45: DMSNES_Local *dmlocalsnes = (DMSNES_Local *) ctx;
46: DM dm;
47: Vec Xloc,Floc;
48: PetscBool transform;
53: SNESGetDM(snes,&dm);
54: DMGetLocalVector(dm,&Xloc);
55: DMGetLocalVector(dm,&Floc);
56: VecZeroEntries(Xloc);
57: VecZeroEntries(Floc);
58: /* Non-conforming routines needs boundary values before G2L */
59: if (dmlocalsnes->boundarylocal) (*dmlocalsnes->boundarylocal)(dm,Xloc,dmlocalsnes->boundarylocalctx);
60: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
61: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
62: /* Need to reset boundary values if we transformed */
63: DMHasBasisTransform(dm, &transform);
64: if (transform && dmlocalsnes->boundarylocal) (*dmlocalsnes->boundarylocal)(dm,Xloc,dmlocalsnes->boundarylocalctx);
65: CHKMEMQ;
66: (*dmlocalsnes->residuallocal)(dm,Xloc,Floc,dmlocalsnes->residuallocalctx);
67: CHKMEMQ;
68: VecZeroEntries(F);
69: DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
70: DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
71: DMRestoreLocalVector(dm,&Floc);
72: DMRestoreLocalVector(dm,&Xloc);
73: {
74: char name[PETSC_MAX_PATH_LEN];
75: char oldname[PETSC_MAX_PATH_LEN];
76: const char *tmp;
77: PetscInt it;
79: SNESGetIterationNumber(snes, &it);
80: PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "Solution, Iterate %d", (int) it);
81: PetscObjectGetName((PetscObject) X, &tmp);
82: PetscStrncpy(oldname, tmp, PETSC_MAX_PATH_LEN-1);
83: PetscObjectSetName((PetscObject) X, name);
84: VecViewFromOptions(X, (PetscObject) snes, "-dmsnes_solution_vec_view");
85: PetscObjectSetName((PetscObject) X, oldname);
86: PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "Residual, Iterate %d", (int) it);
87: PetscObjectSetName((PetscObject) F, name);
88: VecViewFromOptions(F, (PetscObject) snes, "-dmsnes_residual_vec_view");
89: }
90: return 0;
91: }
93: static PetscErrorCode SNESComputeJacobian_DMLocal(SNES snes,Vec X,Mat A,Mat B,void *ctx)
94: {
95: DMSNES_Local *dmlocalsnes = (DMSNES_Local *) ctx;
96: DM dm;
97: Vec Xloc;
98: PetscBool transform;
100: SNESGetDM(snes,&dm);
101: if (dmlocalsnes->jacobianlocal) {
102: DMGetLocalVector(dm,&Xloc);
103: VecZeroEntries(Xloc);
104: /* Non-conforming routines needs boundary values before G2L */
105: if (dmlocalsnes->boundarylocal) (*dmlocalsnes->boundarylocal)(dm,Xloc,dmlocalsnes->boundarylocalctx);
106: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
107: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
108: /* Need to reset boundary values if we transformed */
109: DMHasBasisTransform(dm, &transform);
110: if (transform && dmlocalsnes->boundarylocal) (*dmlocalsnes->boundarylocal)(dm,Xloc,dmlocalsnes->boundarylocalctx);
111: CHKMEMQ;
112: (*dmlocalsnes->jacobianlocal)(dm,Xloc,A,B,dmlocalsnes->jacobianlocalctx);
113: CHKMEMQ;
114: DMRestoreLocalVector(dm,&Xloc);
115: } else {
116: MatFDColoring fdcoloring;
117: PetscObjectQuery((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject*)&fdcoloring);
118: if (!fdcoloring) {
119: ISColoring coloring;
121: DMCreateColoring(dm,dm->coloringtype,&coloring);
122: MatFDColoringCreate(B,coloring,&fdcoloring);
123: ISColoringDestroy(&coloring);
124: switch (dm->coloringtype) {
125: case IS_COLORING_GLOBAL:
126: MatFDColoringSetFunction(fdcoloring,(PetscErrorCode (*)(void))SNESComputeFunction_DMLocal,dmlocalsnes);
127: break;
128: default: SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"No support for coloring type '%s'",ISColoringTypes[dm->coloringtype]);
129: }
130: PetscObjectSetOptionsPrefix((PetscObject)fdcoloring,((PetscObject)dm)->prefix);
131: MatFDColoringSetFromOptions(fdcoloring);
132: MatFDColoringSetUp(B,coloring,fdcoloring);
133: PetscObjectCompose((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject)fdcoloring);
134: PetscObjectDereference((PetscObject)fdcoloring);
136: /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call
137: * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the
138: * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually
139: * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be
140: * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed.
141: */
142: PetscObjectDereference((PetscObject)dm);
143: }
144: MatFDColoringApply(B,fdcoloring,X,snes);
145: }
146: /* This will be redundant if the user called both, but it's too common to forget. */
147: if (A != B) {
148: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
149: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
150: }
151: return 0;
152: }
154: /*@C
155: DMSNESSetFunctionLocal - set a local residual evaluation function. This function is called with local vector
156: containing the local vector information PLUS ghost point information. It should compute a result for all local
157: elements and DMSNES will automatically accumulate the overlapping values.
159: Logically Collective
161: Input Parameters:
162: + dm - DM to associate callback with
163: . func - local residual evaluation
164: - ctx - optional context for local residual evaluation
166: Level: beginner
168: .seealso: DMSNESSetFunction(), DMDASNESSetJacobianLocal(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
169: @*/
170: PetscErrorCode DMSNESSetFunctionLocal(DM dm,PetscErrorCode (*func)(DM,Vec,Vec,void*),void *ctx)
171: {
172: DMSNES sdm;
173: DMSNES_Local *dmlocalsnes;
176: DMGetDMSNESWrite(dm,&sdm);
177: DMLocalSNESGetContext(dm,sdm,&dmlocalsnes);
179: dmlocalsnes->residuallocal = func;
180: dmlocalsnes->residuallocalctx = ctx;
182: DMSNESSetFunction(dm,SNESComputeFunction_DMLocal,dmlocalsnes);
183: if (!sdm->ops->computejacobian) { /* Call us for the Jacobian too, can be overridden by the user. */
184: DMSNESSetJacobian(dm,SNESComputeJacobian_DMLocal,dmlocalsnes);
185: }
186: return 0;
187: }
189: /*@C
190: DMSNESSetBoundaryLocal - set a local boundary value function. This function is called with local vector
191: containing the local vector information PLUS ghost point information. It should insert values into the local
192: vector that do not come from the global vector, such as essential boundary condition data.
194: Logically Collective
196: Input Parameters:
197: + dm - DM to associate callback with
198: . func - local boundary value evaluation
199: - ctx - optional context for local boundary value evaluation
201: Level: intermediate
203: .seealso: DMSNESSetFunctionLocal(), DMDASNESSetJacobianLocal()
204: @*/
205: PetscErrorCode DMSNESSetBoundaryLocal(DM dm,PetscErrorCode (*func)(DM,Vec,void*),void *ctx)
206: {
207: DMSNES sdm;
208: DMSNES_Local *dmlocalsnes;
211: DMGetDMSNESWrite(dm,&sdm);
212: DMLocalSNESGetContext(dm,sdm,&dmlocalsnes);
214: dmlocalsnes->boundarylocal = func;
215: dmlocalsnes->boundarylocalctx = ctx;
217: return 0;
218: }
220: /*@C
221: DMSNESSetJacobianLocal - set a local Jacobian evaluation function
223: Logically Collective
225: Input Parameters:
226: + dm - DM to associate callback with
227: . func - local Jacobian evaluation
228: - ctx - optional context for local Jacobian evaluation
230: Level: beginner
232: .seealso: DMSNESSetJacobian(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
233: @*/
234: PetscErrorCode DMSNESSetJacobianLocal(DM dm,PetscErrorCode (*func)(DM,Vec,Mat,Mat,void*),void *ctx)
235: {
236: DMSNES sdm;
237: DMSNES_Local *dmlocalsnes;
240: DMGetDMSNESWrite(dm,&sdm);
241: DMLocalSNESGetContext(dm,sdm,&dmlocalsnes);
243: dmlocalsnes->jacobianlocal = func;
244: dmlocalsnes->jacobianlocalctx = ctx;
246: DMSNESSetJacobian(dm,SNESComputeJacobian_DMLocal,dmlocalsnes);
247: return 0;
248: }
250: /*@C
251: DMSNESGetFunctionLocal - get the local residual evaluation function information set with DMSNESSetFunctionLocal.
253: Not Collective
255: Input Parameter:
256: . dm - DM with the associated callback
258: Output Parameters:
259: + func - local residual evaluation
260: - ctx - context for local residual evaluation
262: Level: beginner
264: .seealso: DMSNESSetFunction(), DMSNESSetFunctionLocal(), DMDASNESSetJacobianLocal(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
265: @*/
266: PetscErrorCode DMSNESGetFunctionLocal(DM dm,PetscErrorCode (**func)(DM,Vec,Vec,void*),void **ctx)
267: {
268: DMSNES sdm;
269: DMSNES_Local *dmlocalsnes;
272: DMGetDMSNES(dm,&sdm);
273: DMLocalSNESGetContext(dm,sdm,&dmlocalsnes);
274: if (func) *func = dmlocalsnes->residuallocal;
275: if (ctx) *ctx = dmlocalsnes->residuallocalctx;
276: return 0;
277: }
279: /*@C
280: DMSNESGetBoundaryLocal - get the local boundary value function set with DMSNESSetBoundaryLocal.
282: Not Collective
284: Input Parameter:
285: . dm - DM with the associated callback
287: Output Parameters:
288: + func - local boundary value evaluation
289: - ctx - context for local boundary value evaluation
291: Level: intermediate
293: .seealso: DMSNESSetFunctionLocal(), DMSNESSetBoundaryLocal(), DMDASNESSetJacobianLocal()
294: @*/
295: PetscErrorCode DMSNESGetBoundaryLocal(DM dm,PetscErrorCode (**func)(DM,Vec,void*),void **ctx)
296: {
297: DMSNES sdm;
298: DMSNES_Local *dmlocalsnes;
301: DMGetDMSNES(dm,&sdm);
302: DMLocalSNESGetContext(dm,sdm,&dmlocalsnes);
303: if (func) *func = dmlocalsnes->boundarylocal;
304: if (ctx) *ctx = dmlocalsnes->boundarylocalctx;
305: return 0;
306: }
308: /*@C
309: DMSNESGetJacobianLocal - the local Jacobian evaluation function set with DMSNESSetJacobianLocal.
311: Logically Collective
313: Input Parameter:
314: . dm - DM with the associated callback
316: Output Parameters:
317: + func - local Jacobian evaluation
318: - ctx - context for local Jacobian evaluation
320: Level: beginner
322: .seealso: DMSNESSetJacobianLocal(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
323: @*/
324: PetscErrorCode DMSNESGetJacobianLocal(DM dm,PetscErrorCode (**func)(DM,Vec,Mat,Mat,void*),void **ctx)
325: {
326: DMSNES sdm;
327: DMSNES_Local *dmlocalsnes;
330: DMGetDMSNES(dm,&sdm);
331: DMLocalSNESGetContext(dm,sdm,&dmlocalsnes);
332: if (func) *func = dmlocalsnes->jacobianlocal;
333: if (ctx) *ctx = dmlocalsnes->jacobianlocalctx;
334: return 0;
335: }