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: PetscNew((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: PetscNew((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:
129: SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "No support for coloring type '%s'", ISColoringTypes[dm->coloringtype]);
130: }
131: PetscObjectSetOptionsPrefix((PetscObject)fdcoloring, ((PetscObject)dm)->prefix);
132: MatFDColoringSetFromOptions(fdcoloring);
133: MatFDColoringSetUp(B, coloring, fdcoloring);
134: PetscObjectCompose((PetscObject)dm, "DMDASNES_FDCOLORING", (PetscObject)fdcoloring);
135: PetscObjectDereference((PetscObject)fdcoloring);
137: /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call
138: * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the
139: * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually
140: * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be
141: * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed.
142: */
143: PetscObjectDereference((PetscObject)dm);
144: }
145: MatFDColoringApply(B, fdcoloring, X, snes);
146: }
147: /* This will be redundant if the user called both, but it's too common to forget. */
148: if (A != B) {
149: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
150: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
151: }
152: return 0;
153: }
155: /*@C
156: DMSNESSetFunctionLocal - set a local residual evaluation function. This function is called with local vector
157: containing the local vector information PLUS ghost point information. It should compute a result for all local
158: elements and `DMSNES` will automatically accumulate the overlapping values.
160: Logically Collective
162: Input Parameters:
163: + dm - `DM` to associate callback with
164: . func - local residual evaluation
165: - ctx - optional context for local residual evaluation
167: Level: advanced
169: .seealso: `DMSNESSetFunction()`, `DMDASNESSetJacobianLocal()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`
170: @*/
171: PetscErrorCode DMSNESSetFunctionLocal(DM dm, PetscErrorCode (*func)(DM, Vec, Vec, void *), void *ctx)
172: {
173: DMSNES sdm;
174: DMSNES_Local *dmlocalsnes;
177: DMGetDMSNESWrite(dm, &sdm);
178: DMLocalSNESGetContext(dm, sdm, &dmlocalsnes);
180: dmlocalsnes->residuallocal = func;
181: dmlocalsnes->residuallocalctx = ctx;
183: DMSNESSetFunction(dm, SNESComputeFunction_DMLocal, dmlocalsnes);
184: if (!sdm->ops->computejacobian) { /* Call us for the Jacobian too, can be overridden by the user. */
185: DMSNESSetJacobian(dm, SNESComputeJacobian_DMLocal, dmlocalsnes);
186: }
187: return 0;
188: }
190: /*@C
191: DMSNESSetBoundaryLocal - set a local boundary value function. This function is called with local vector
192: containing the local vector information PLUS ghost point information. It should insert values into the local
193: vector that do not come from the global vector, such as essential boundary condition data.
195: Logically Collective
197: Input Parameters:
198: + dm - `DM` to associate callback with
199: . func - local boundary value evaluation
200: - ctx - optional context for local boundary value evaluation
202: Level: advanced
204: .seealso: `DMSNESSetFunctionLocal()`, `DMDASNESSetJacobianLocal()`
205: @*/
206: PetscErrorCode DMSNESSetBoundaryLocal(DM dm, PetscErrorCode (*func)(DM, Vec, void *), void *ctx)
207: {
208: DMSNES sdm;
209: DMSNES_Local *dmlocalsnes;
212: DMGetDMSNESWrite(dm, &sdm);
213: DMLocalSNESGetContext(dm, sdm, &dmlocalsnes);
215: dmlocalsnes->boundarylocal = func;
216: dmlocalsnes->boundarylocalctx = ctx;
218: return 0;
219: }
221: /*@C
222: DMSNESSetJacobianLocal - set a local Jacobian evaluation function
224: Logically Collective
226: Input Parameters:
227: + dm - DM to associate callback with
228: . func - local Jacobian evaluation
229: - ctx - optional context for local Jacobian evaluation
231: Level: advanced
233: .seealso: `DMSNESSetJacobian()`, `DMDASNESSetJacobian()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`
234: @*/
235: PetscErrorCode DMSNESSetJacobianLocal(DM dm, PetscErrorCode (*func)(DM, Vec, Mat, Mat, void *), void *ctx)
236: {
237: DMSNES sdm;
238: DMSNES_Local *dmlocalsnes;
241: DMGetDMSNESWrite(dm, &sdm);
242: DMLocalSNESGetContext(dm, sdm, &dmlocalsnes);
244: dmlocalsnes->jacobianlocal = func;
245: dmlocalsnes->jacobianlocalctx = ctx;
247: DMSNESSetJacobian(dm, SNESComputeJacobian_DMLocal, dmlocalsnes);
248: return 0;
249: }
251: /*@C
252: DMSNESGetFunctionLocal - get the local residual evaluation function information set with `DMSNESSetFunctionLocal()`.
254: Not Collective
256: Input Parameter:
257: . dm - `DM` with the associated callback
259: Output Parameters:
260: + func - local residual evaluation
261: - ctx - context for local residual evaluation
263: Level: beginner
265: .seealso: `DMSNESSetFunction()`, `DMSNESSetFunctionLocal()`, `DMDASNESSetJacobianLocal()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`
266: @*/
267: PetscErrorCode DMSNESGetFunctionLocal(DM dm, PetscErrorCode (**func)(DM, Vec, Vec, void *), void **ctx)
268: {
269: DMSNES sdm;
270: DMSNES_Local *dmlocalsnes;
273: DMGetDMSNES(dm, &sdm);
274: DMLocalSNESGetContext(dm, sdm, &dmlocalsnes);
275: if (func) *func = dmlocalsnes->residuallocal;
276: if (ctx) *ctx = dmlocalsnes->residuallocalctx;
277: return 0;
278: }
280: /*@C
281: DMSNESGetBoundaryLocal - get the local boundary value function set with `DMSNESSetBoundaryLocal()`.
283: Not Collective
285: Input Parameter:
286: . dm - `DM` with the associated callback
288: Output Parameters:
289: + func - local boundary value evaluation
290: - ctx - context for local boundary value evaluation
292: Level: intermediate
294: .seealso: `DMSNESSetFunctionLocal()`, `DMSNESSetBoundaryLocal()`, `DMDASNESSetJacobianLocal()`
295: @*/
296: PetscErrorCode DMSNESGetBoundaryLocal(DM dm, PetscErrorCode (**func)(DM, Vec, void *), void **ctx)
297: {
298: DMSNES sdm;
299: DMSNES_Local *dmlocalsnes;
302: DMGetDMSNES(dm, &sdm);
303: DMLocalSNESGetContext(dm, sdm, &dmlocalsnes);
304: if (func) *func = dmlocalsnes->boundarylocal;
305: if (ctx) *ctx = dmlocalsnes->boundarylocalctx;
306: return 0;
307: }
309: /*@C
310: DMSNESGetJacobianLocal - the local Jacobian evaluation function set with `DMSNESSetJacobianLocal()`.
312: Logically Collective
314: Input Parameter:
315: . dm - `DM` with the associated callback
317: Output Parameters:
318: + func - local Jacobian evaluation
319: - ctx - context for local Jacobian evaluation
321: Level: beginner
323: .seealso: `DMSNESSetJacobianLocal()`, `DMDASNESSetJacobian()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`
324: @*/
325: PetscErrorCode DMSNESGetJacobianLocal(DM dm, PetscErrorCode (**func)(DM, Vec, Mat, Mat, void *), void **ctx)
326: {
327: DMSNES sdm;
328: DMSNES_Local *dmlocalsnes;
331: DMGetDMSNES(dm, &sdm);
332: DMLocalSNESGetContext(dm, sdm, &dmlocalsnes);
333: if (func) *func = dmlocalsnes->jacobianlocal;
334: if (ctx) *ctx = dmlocalsnes->jacobianlocalctx;
335: return 0;
336: }