Actual source code: dmlocalsnes.c
1: #include <petsc/private/dmimpl.h>
2: #include <petsc/private/snesimpl.h>
4: typedef struct {
5: PetscErrorCode (*objectivelocal)(DM, Vec, PetscReal *, PetscCtx);
6: PetscErrorCode (*residuallocal)(DM, Vec, Vec, PetscCtx);
7: PetscErrorCode (*jacobianlocal)(DM, Vec, Mat, Mat, PetscCtx);
8: PetscErrorCode (*boundarylocal)(DM, Vec, PetscCtx);
9: PetscCtx objectivelocalctx;
10: PetscCtx residuallocalctx;
11: PetscCtx jacobianlocalctx;
12: PetscCtx boundarylocalctx;
13: } DMSNES_Local;
15: static PetscErrorCode DMSNESDestroy_DMLocal(DMSNES sdm)
16: {
17: PetscFunctionBegin;
18: PetscCall(PetscFree(sdm->data));
19: sdm->data = NULL;
20: PetscFunctionReturn(PETSC_SUCCESS);
21: }
23: static PetscErrorCode DMSNESDuplicate_DMLocal(DMSNES oldsdm, DMSNES sdm)
24: {
25: PetscFunctionBegin;
26: if (sdm->data != oldsdm->data) {
27: PetscCall(PetscFree(sdm->data));
28: PetscCall(PetscNew((DMSNES_Local **)&sdm->data));
29: if (oldsdm->data) PetscCall(PetscMemcpy(sdm->data, oldsdm->data, sizeof(DMSNES_Local)));
30: }
31: PetscFunctionReturn(PETSC_SUCCESS);
32: }
34: static PetscErrorCode DMLocalSNESGetContext(DM dm, DMSNES sdm, DMSNES_Local **dmlocalsnes)
35: {
36: PetscFunctionBegin;
37: *dmlocalsnes = NULL;
38: if (!sdm->data) {
39: PetscCall(PetscNew((DMSNES_Local **)&sdm->data));
41: sdm->ops->destroy = DMSNESDestroy_DMLocal;
42: sdm->ops->duplicate = DMSNESDuplicate_DMLocal;
43: }
44: *dmlocalsnes = (DMSNES_Local *)sdm->data;
45: PetscFunctionReturn(PETSC_SUCCESS);
46: }
48: static PetscErrorCode SNESComputeObjective_DMLocal(SNES snes, Vec X, PetscReal *obj, PetscCtx ctx)
49: {
50: DMSNES_Local *dmlocalsnes = (DMSNES_Local *)ctx;
51: DM dm;
52: Vec Xloc;
53: PetscBool transform;
55: PetscFunctionBegin;
58: PetscCall(SNESGetDM(snes, &dm));
59: PetscCall(DMGetLocalVector(dm, &Xloc));
60: PetscCall(VecZeroEntries(Xloc));
61: /* Non-conforming routines needs boundary values before G2L */
62: if (dmlocalsnes->boundarylocal) PetscCall((*dmlocalsnes->boundarylocal)(dm, Xloc, dmlocalsnes->boundarylocalctx));
63: PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, Xloc));
64: PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, Xloc));
65: /* Need to reset boundary values if we transformed */
66: PetscCall(DMHasBasisTransform(dm, &transform));
67: if (transform && dmlocalsnes->boundarylocal) PetscCall((*dmlocalsnes->boundarylocal)(dm, Xloc, dmlocalsnes->boundarylocalctx));
68: CHKMEMQ;
69: PetscCall((*dmlocalsnes->objectivelocal)(dm, Xloc, obj, dmlocalsnes->objectivelocalctx));
70: CHKMEMQ;
71: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, obj, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)snes)));
72: PetscCall(DMRestoreLocalVector(dm, &Xloc));
73: PetscFunctionReturn(PETSC_SUCCESS);
74: }
76: static PetscErrorCode SNESComputeFunction_DMLocal(SNES snes, Vec X, Vec F, PetscCtx ctx)
77: {
78: DMSNES_Local *dmlocalsnes = (DMSNES_Local *)ctx;
79: DM dm;
80: Vec Xloc, Floc;
81: PetscBool transform;
83: PetscFunctionBegin;
87: PetscCall(SNESGetDM(snes, &dm));
88: PetscCall(DMGetLocalVector(dm, &Xloc));
89: PetscCall(DMGetLocalVector(dm, &Floc));
90: PetscCall(VecZeroEntries(Xloc));
91: PetscCall(VecZeroEntries(Floc));
92: /* Non-conforming routines needs boundary values before G2L */
93: if (dmlocalsnes->boundarylocal) PetscCall((*dmlocalsnes->boundarylocal)(dm, Xloc, dmlocalsnes->boundarylocalctx));
94: PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, Xloc));
95: PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, Xloc));
96: /* Need to reset boundary values if we transformed */
97: PetscCall(DMHasBasisTransform(dm, &transform));
98: if (transform && dmlocalsnes->boundarylocal) PetscCall((*dmlocalsnes->boundarylocal)(dm, Xloc, dmlocalsnes->boundarylocalctx));
99: CHKMEMQ;
100: PetscCall((*dmlocalsnes->residuallocal)(dm, Xloc, Floc, dmlocalsnes->residuallocalctx));
101: CHKMEMQ;
102: PetscCall(VecZeroEntries(F));
103: PetscCall(DMLocalToGlobalBegin(dm, Floc, ADD_VALUES, F));
104: PetscCall(DMLocalToGlobalEnd(dm, Floc, ADD_VALUES, F));
105: PetscCall(DMRestoreLocalVector(dm, &Floc));
106: PetscCall(DMRestoreLocalVector(dm, &Xloc));
107: {
108: char name[PETSC_MAX_PATH_LEN];
109: char oldname[PETSC_MAX_PATH_LEN];
110: const char *tmp;
111: PetscInt it;
113: PetscCall(SNESGetIterationNumber(snes, &it));
114: PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "Solution, Iterate %" PetscInt_FMT, it));
115: PetscCall(PetscObjectGetName((PetscObject)X, &tmp));
116: PetscCall(PetscStrncpy(oldname, tmp, PETSC_MAX_PATH_LEN - 1));
117: PetscCall(PetscObjectSetName((PetscObject)X, name));
118: PetscCall(VecViewFromOptions(X, (PetscObject)snes, "-dmsnes_solution_vec_view"));
119: PetscCall(PetscObjectSetName((PetscObject)X, oldname));
120: PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "Residual, Iterate %" PetscInt_FMT, it));
121: PetscCall(PetscObjectSetName((PetscObject)F, name));
122: PetscCall(VecViewFromOptions(F, (PetscObject)snes, "-dmsnes_residual_vec_view"));
123: }
124: PetscFunctionReturn(PETSC_SUCCESS);
125: }
127: static PetscErrorCode SNESComputeJacobian_DMLocal(SNES snes, Vec X, Mat A, Mat B, PetscCtx ctx)
128: {
129: DMSNES_Local *dmlocalsnes = (DMSNES_Local *)ctx;
130: DM dm;
131: Vec Xloc;
132: PetscBool transform;
134: PetscFunctionBegin;
135: PetscCall(SNESGetDM(snes, &dm));
136: if (dmlocalsnes->jacobianlocal) {
137: PetscCall(DMGetLocalVector(dm, &Xloc));
138: PetscCall(VecZeroEntries(Xloc));
139: /* Non-conforming routines needs boundary values before G2L */
140: if (dmlocalsnes->boundarylocal) PetscCall((*dmlocalsnes->boundarylocal)(dm, Xloc, dmlocalsnes->boundarylocalctx));
141: PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, Xloc));
142: PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, Xloc));
143: /* Need to reset boundary values if we transformed */
144: PetscCall(DMHasBasisTransform(dm, &transform));
145: if (transform && dmlocalsnes->boundarylocal) PetscCall((*dmlocalsnes->boundarylocal)(dm, Xloc, dmlocalsnes->boundarylocalctx));
146: CHKMEMQ;
147: PetscCall((*dmlocalsnes->jacobianlocal)(dm, Xloc, A, B, dmlocalsnes->jacobianlocalctx));
148: CHKMEMQ;
149: PetscCall(DMRestoreLocalVector(dm, &Xloc));
150: } else {
151: MatFDColoring fdcoloring;
152: PetscCall(PetscObjectQuery((PetscObject)dm, "DMDASNES_FDCOLORING", (PetscObject *)&fdcoloring));
153: if (!fdcoloring) {
154: ISColoring coloring;
156: PetscCall(DMCreateColoring(dm, dm->coloringtype, &coloring));
157: PetscCall(MatFDColoringCreate(B, coloring, &fdcoloring));
158: PetscCall(ISColoringDestroy(&coloring));
159: switch (dm->coloringtype) {
160: case IS_COLORING_GLOBAL:
161: PetscCall(MatFDColoringSetFunction(fdcoloring, (MatFDColoringFn *)SNESComputeFunction_DMLocal, dmlocalsnes));
162: break;
163: default:
164: SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "No support for coloring type '%s'", ISColoringTypes[dm->coloringtype]);
165: }
166: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)fdcoloring, ((PetscObject)dm)->prefix));
167: PetscCall(MatFDColoringSetFromOptions(fdcoloring));
168: PetscCall(MatFDColoringSetUp(B, coloring, fdcoloring));
169: PetscCall(PetscObjectCompose((PetscObject)dm, "DMDASNES_FDCOLORING", (PetscObject)fdcoloring));
170: PetscCall(PetscObjectDereference((PetscObject)fdcoloring));
172: /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call
173: * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the
174: * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually
175: * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be
176: * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed.
177: */
178: PetscCall(PetscObjectDereference((PetscObject)dm));
179: }
180: PetscCall(MatFDColoringApply(B, fdcoloring, X, snes));
181: }
182: /* This will be redundant if the user called both, but it's too common to forget. */
183: if (A != B) {
184: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
185: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
186: }
187: PetscFunctionReturn(PETSC_SUCCESS);
188: }
190: /*@C
191: DMSNESSetObjectiveLocal - set a local objective evaluation function. This function is called with local vector
192: containing the local vector information PLUS ghost point information. It should compute a result for all local
193: elements and `DMSNES` will automatically accumulate the overlapping values.
195: Logically Collective
197: Input Parameters:
198: + dm - `DM` to associate callback with
199: . func - local objective evaluation
200: - ctx - optional context for local objective function evaluation
202: Calling sequence of func:
203: + dm - the `DM`
204: . x - the location where the objective is to be evaluated
205: . obj - the value of the objective function
206: - ctx - optional context for the local objective function evaluation
208: Level: advanced
210: .seealso: `DMSNESSetFunctionLocal()`, `DMSNESSetJacobianLocal()`
211: @*/
212: PetscErrorCode DMSNESSetObjectiveLocal(DM dm, PetscErrorCode (*func)(DM dm, Vec x, PetscReal *obj, PetscCtx ctx), PetscCtx ctx)
213: {
214: DMSNES sdm;
215: DMSNES_Local *dmlocalsnes;
217: PetscFunctionBegin;
219: PetscCall(DMGetDMSNESWrite(dm, &sdm));
220: PetscCall(DMLocalSNESGetContext(dm, sdm, &dmlocalsnes));
222: dmlocalsnes->objectivelocal = func;
223: dmlocalsnes->objectivelocalctx = ctx;
225: PetscCall(DMSNESSetObjective(dm, SNESComputeObjective_DMLocal, dmlocalsnes));
226: PetscFunctionReturn(PETSC_SUCCESS);
227: }
229: /*@C
230: DMSNESSetFunctionLocal - set a local residual evaluation function. This function is called with local vector
231: containing the local vector information PLUS ghost point information. It should compute a result for all local
232: elements and `DMSNES` will automatically accumulate the overlapping values.
234: Logically Collective
236: Input Parameters:
237: + dm - `DM` to associate callback with
238: . func - local residual evaluation
239: - ctx - optional context for local residual evaluation
241: Calling sequence of `func`:
242: + dm - `DM` for the function
243: . x - vector to state at which to evaluate residual
244: . f - vector to hold the function evaluation
245: - ctx - optional context passed above
247: Level: advanced
249: .seealso: [](ch_snes), `DMSNESSetFunction()`, `DMSNESSetJacobianLocal()`
250: @*/
251: PetscErrorCode DMSNESSetFunctionLocal(DM dm, PetscErrorCode (*func)(DM dm, Vec x, Vec f, PetscCtx ctx), PetscCtx ctx) PeNSS
252: {
253: DMSNES sdm;
254: DMSNES_Local *dmlocalsnes;
256: PetscFunctionBegin;
258: PetscCall(DMGetDMSNESWrite(dm, &sdm));
259: PetscCall(DMLocalSNESGetContext(dm, sdm, &dmlocalsnes));
261: dmlocalsnes->residuallocal = func;
262: dmlocalsnes->residuallocalctx = ctx;
264: PetscCall(DMSNESSetFunction(dm, SNESComputeFunction_DMLocal, dmlocalsnes));
265: if (!sdm->ops->computejacobian) { /* Call us for the Jacobian too, can be overridden by the user. */
266: PetscCall(DMSNESSetJacobian(dm, SNESComputeJacobian_DMLocal, dmlocalsnes));
267: }
268: PetscFunctionReturn(PETSC_SUCCESS);
269: }
271: /*@C
272: DMSNESSetBoundaryLocal - set a function to insert, for example, essential boundary conditions into a ghosted solution vector
274: Logically Collective
276: Input Parameters:
277: + dm - `DM` to associate callback with
278: . func - local boundary value evaluation
279: - ctx - optional context for local boundary value evaluation
281: Calling sequence of `func`:
282: + dm - the `DM` context
283: . X - ghosted solution vector, appropriate locations (such as essential boundary condition nodes) should be filled
284: - ctx - option context passed in `DMSNESSetBoundaryLocal()`
286: Level: advanced
288: .seealso: [](ch_snes), `DMSNESSetObjectiveLocal()`, `DMSNESSetFunctionLocal()`, `DMSNESSetJacobianLocal()`
289: @*/
290: PetscErrorCode DMSNESSetBoundaryLocal(DM dm, PetscErrorCode (*func)(DM dm, Vec X, PetscCtx ctx), PetscCtx ctx)
291: {
292: DMSNES sdm;
293: DMSNES_Local *dmlocalsnes;
295: PetscFunctionBegin;
297: PetscCall(DMGetDMSNESWrite(dm, &sdm));
298: PetscCall(DMLocalSNESGetContext(dm, sdm, &dmlocalsnes));
300: dmlocalsnes->boundarylocal = func;
301: dmlocalsnes->boundarylocalctx = ctx;
302: PetscFunctionReturn(PETSC_SUCCESS);
303: }
305: /*@C
306: DMSNESSetJacobianLocal - set a local Jacobian evaluation function
308: Logically Collective
310: Input Parameters:
311: + dm - `DM` to associate callback with
312: . func - local Jacobian evaluation
313: - ctx - optional context for local Jacobian evaluation
315: Calling sequence of `func`:
316: + dm - the `DM` context
317: . X - current solution vector (ghosted or not?)
318: . J - the Jacobian
319: . Jp - approximate Jacobian used to compute the preconditioner, often `J`
320: - ctx - a user provided context
322: Level: advanced
324: .seealso: [](ch_snes), `DMSNESSetObjectiveLocal()`, `DMSNESSetFunctionLocal()`, `DMSNESSetBoundaryLocal()`
325: @*/
326: PetscErrorCode DMSNESSetJacobianLocal(DM dm, PetscErrorCode (*func)(DM dm, Vec X, Mat J, Mat Jp, PetscCtx ctx), PetscCtx ctx) PeNSS
327: {
328: DMSNES sdm;
329: DMSNES_Local *dmlocalsnes;
331: PetscFunctionBegin;
333: PetscCall(DMGetDMSNESWrite(dm, &sdm));
334: PetscCall(DMLocalSNESGetContext(dm, sdm, &dmlocalsnes));
336: dmlocalsnes->jacobianlocal = func;
337: dmlocalsnes->jacobianlocalctx = ctx;
339: PetscCall(DMSNESSetJacobian(dm, SNESComputeJacobian_DMLocal, dmlocalsnes));
340: PetscFunctionReturn(PETSC_SUCCESS);
341: }
343: /*@C
344: DMSNESGetObjectiveLocal - get the local objective evaluation function information set with `DMSNESSetObjectiveLocal()`.
346: Not Collective
348: Input Parameter:
349: . dm - `DM` with the associated callback
351: Output Parameters:
352: + func - local objective evaluation
353: - ctx - context for local residual evaluation
355: Calling sequence of func:
356: + dm - the `DM`
357: . x - the location where the objective function is to be evaluated
358: . obj - the value of the objective function
359: - ctx - optional context for the local objective function evaluation
361: Level: beginner
363: .seealso: `DMSNESSetObjective()`, `DMSNESSetObjectiveLocal()`, `DMSNESSetFunctionLocal()`
364: @*/
365: PetscErrorCode DMSNESGetObjectiveLocal(DM dm, PetscErrorCode (**func)(DM dm, Vec x, PetscReal *obj, PetscCtx ctx), PetscCtxRt ctx) PeNSS
366: {
367: DMSNES sdm;
368: DMSNES_Local *dmlocalsnes;
370: PetscFunctionBegin;
372: PetscCall(DMGetDMSNES(dm, &sdm));
373: PetscCall(DMLocalSNESGetContext(dm, sdm, &dmlocalsnes));
374: if (func) *func = dmlocalsnes->objectivelocal;
375: if (ctx) *(void **)ctx = dmlocalsnes->objectivelocalctx;
376: PetscFunctionReturn(PETSC_SUCCESS);
377: }
379: /*@C
380: DMSNESGetFunctionLocal - get the local residual evaluation function information set with `DMSNESSetFunctionLocal()`.
382: Not Collective
384: Input Parameter:
385: . dm - `DM` with the associated callback
387: Output Parameters:
388: + func - local residual evaluation
389: - ctx - context for local residual evaluation
391: Calling sequence of `func`:
392: + dm - `DM` for the function
393: . x - vector to state at which to evaluate residual
394: . f - vector to hold the function evaluation
395: - ctx - optional context passed above
397: Level: beginner
399: .seealso: [](ch_snes), `DMSNESSetFunction()`, `DMSNESSetFunctionLocal()`, `DMSNESSetJacobianLocal()`
400: @*/
401: PetscErrorCode DMSNESGetFunctionLocal(DM dm, PetscErrorCode (**func)(DM dm, Vec x, Vec f, PetscCtx ctx), PetscCtxRt ctx)
402: {
403: DMSNES sdm;
404: DMSNES_Local *dmlocalsnes;
406: PetscFunctionBegin;
408: PetscCall(DMGetDMSNES(dm, &sdm));
409: PetscCall(DMLocalSNESGetContext(dm, sdm, &dmlocalsnes));
410: if (func) *func = dmlocalsnes->residuallocal;
411: if (ctx) *(void **)ctx = dmlocalsnes->residuallocalctx;
412: PetscFunctionReturn(PETSC_SUCCESS);
413: }
415: /*@C
416: DMSNESGetBoundaryLocal - get the local boundary value function set with `DMSNESSetBoundaryLocal()`.
418: Not Collective
420: Input Parameter:
421: . dm - `DM` with the associated callback
423: Output Parameters:
424: + func - local boundary value evaluation
425: - ctx - context for local boundary value evaluation
427: Calling sequence of `func`:
428: + dm - the `DM` context
429: . X - ghosted solution vector, appropriate locations (such as essential boundary condition nodes) should be filled
430: - ctx - option context passed in `DMSNESSetBoundaryLocal()`
432: Level: intermediate
434: .seealso: [](ch_snes), `DMSNESSetFunctionLocal()`, `DMSNESSetBoundaryLocal()`, `DMSNESSetJacobianLocal()`
435: @*/
436: PetscErrorCode DMSNESGetBoundaryLocal(DM dm, PetscErrorCode (**func)(DM dm, Vec X, PetscCtx ctx), PetscCtxRt ctx)
437: {
438: DMSNES sdm;
439: DMSNES_Local *dmlocalsnes;
441: PetscFunctionBegin;
443: PetscCall(DMGetDMSNES(dm, &sdm));
444: PetscCall(DMLocalSNESGetContext(dm, sdm, &dmlocalsnes));
445: if (func) *func = dmlocalsnes->boundarylocal;
446: if (ctx) *(void **)ctx = dmlocalsnes->boundarylocalctx;
447: PetscFunctionReturn(PETSC_SUCCESS);
448: }
450: /*@C
451: DMSNESGetJacobianLocal - the local Jacobian evaluation function set with `DMSNESSetJacobianLocal()`.
453: Logically Collective
455: Input Parameter:
456: . dm - `DM` with the associated callback
458: Output Parameters:
459: + func - local Jacobian evaluation
460: - ctx - context for local Jacobian evaluation
462: Calling sequence of `func`:
463: + dm - the `DM` context
464: . X - current solution vector (ghosted or not?)
465: . J - the Jacobian
466: . Jp - approximate Jacobian used to compute the preconditioner, often `J`
467: - ctx - a user provided context
469: Level: beginner
471: .seealso: [](ch_snes), `DMSNESSetJacobianLocal()`, `DMSNESSetJacobian()`
472: @*/
473: PetscErrorCode DMSNESGetJacobianLocal(DM dm, PetscErrorCode (**func)(DM dm, Vec X, Mat J, Mat Jp, PetscCtx ctx), PetscCtxRt ctx)
474: {
475: DMSNES sdm;
476: DMSNES_Local *dmlocalsnes;
478: PetscFunctionBegin;
480: PetscCall(DMGetDMSNES(dm, &sdm));
481: PetscCall(DMLocalSNESGetContext(dm, sdm, &dmlocalsnes));
482: if (func) *func = dmlocalsnes->jacobianlocal;
483: if (ctx) *(void **)ctx = dmlocalsnes->jacobianlocalctx;
484: PetscFunctionReturn(PETSC_SUCCESS);
485: }