Actual source code: dmlocalsnes.c
petsc-3.14.6 2021-03-30
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: {
18: PetscFree(sdm->data);
19: sdm->data = NULL;
20: return(0);
21: }
23: static PetscErrorCode DMSNESDuplicate_DMLocal(DMSNES oldsdm,DMSNES sdm)
24: {
28: if (sdm->data != oldsdm->data) {
29: PetscFree(sdm->data);
30: PetscNewLog(sdm,(DMSNES_Local**)&sdm->data);
31: if (oldsdm->data) {PetscMemcpy(sdm->data,oldsdm->data,sizeof(DMSNES_Local));}
32: }
33: return(0);
34: }
36: static PetscErrorCode DMLocalSNESGetContext(DM dm,DMSNES sdm,DMSNES_Local **dmlocalsnes)
37: {
41: *dmlocalsnes = NULL;
42: if (!sdm->data) {
43: PetscNewLog(dm,(DMSNES_Local**)&sdm->data);
45: sdm->ops->destroy = DMSNESDestroy_DMLocal;
46: sdm->ops->duplicate = DMSNESDuplicate_DMLocal;
47: }
48: *dmlocalsnes = (DMSNES_Local*)sdm->data;
49: return(0);
50: }
52: static PetscErrorCode SNESComputeFunction_DMLocal(SNES snes,Vec X,Vec F,void *ctx)
53: {
54: DMSNES_Local *dmlocalsnes = (DMSNES_Local *) ctx;
55: DM dm;
56: Vec Xloc,Floc;
57: PetscBool transform;
64: SNESGetDM(snes,&dm);
65: DMGetLocalVector(dm,&Xloc);
66: DMGetLocalVector(dm,&Floc);
67: VecZeroEntries(Xloc);
68: VecZeroEntries(Floc);
69: /* Non-conforming routines needs boundary values before G2L */
70: if (dmlocalsnes->boundarylocal) {(*dmlocalsnes->boundarylocal)(dm,Xloc,dmlocalsnes->boundarylocalctx);}
71: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
72: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
73: /* Need to reset boundary values if we transformed */
74: DMHasBasisTransform(dm, &transform);
75: if (transform && dmlocalsnes->boundarylocal) {(*dmlocalsnes->boundarylocal)(dm,Xloc,dmlocalsnes->boundarylocalctx);}
76: CHKMEMQ;
77: (*dmlocalsnes->residuallocal)(dm,Xloc,Floc,dmlocalsnes->residuallocalctx);
78: CHKMEMQ;
79: VecZeroEntries(F);
80: DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
81: DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
82: DMRestoreLocalVector(dm,&Floc);
83: DMRestoreLocalVector(dm,&Xloc);
84: {
85: char name[PETSC_MAX_PATH_LEN];
86: char oldname[PETSC_MAX_PATH_LEN];
87: const char *tmp;
88: PetscInt it;
90: SNESGetIterationNumber(snes, &it);
91: PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "Solution, Iterate %d", (int) it);
92: PetscObjectGetName((PetscObject) X, &tmp);
93: PetscStrncpy(oldname, tmp, PETSC_MAX_PATH_LEN-1);
94: PetscObjectSetName((PetscObject) X, name);
95: VecViewFromOptions(X, (PetscObject) snes, "-dmsnes_solution_vec_view");
96: PetscObjectSetName((PetscObject) X, oldname);
97: PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "Residual, Iterate %d", (int) it);
98: PetscObjectSetName((PetscObject) F, name);
99: VecViewFromOptions(F, (PetscObject) snes, "-dmsnes_residual_vec_view");
100: }
101: return(0);
102: }
104: static PetscErrorCode SNESComputeJacobian_DMLocal(SNES snes,Vec X,Mat A,Mat B,void *ctx)
105: {
106: DMSNES_Local *dmlocalsnes = (DMSNES_Local *) ctx;
107: DM dm;
108: Vec Xloc;
109: PetscBool transform;
113: SNESGetDM(snes,&dm);
114: if (dmlocalsnes->jacobianlocal) {
115: DMGetLocalVector(dm,&Xloc);
116: VecZeroEntries(Xloc);
117: /* Non-conforming routines needs boundary values before G2L */
118: if (dmlocalsnes->boundarylocal) {(*dmlocalsnes->boundarylocal)(dm,Xloc,dmlocalsnes->boundarylocalctx);}
119: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
120: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
121: /* Need to reset boundary values if we transformed */
122: DMHasBasisTransform(dm, &transform);
123: if (transform && dmlocalsnes->boundarylocal) {(*dmlocalsnes->boundarylocal)(dm,Xloc,dmlocalsnes->boundarylocalctx);}
124: CHKMEMQ;
125: (*dmlocalsnes->jacobianlocal)(dm,Xloc,A,B,dmlocalsnes->jacobianlocalctx);
126: CHKMEMQ;
127: DMRestoreLocalVector(dm,&Xloc);
128: } else {
129: MatFDColoring fdcoloring;
130: PetscObjectQuery((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject*)&fdcoloring);
131: if (!fdcoloring) {
132: ISColoring coloring;
134: DMCreateColoring(dm,dm->coloringtype,&coloring);
135: MatFDColoringCreate(B,coloring,&fdcoloring);
136: ISColoringDestroy(&coloring);
137: switch (dm->coloringtype) {
138: case IS_COLORING_GLOBAL:
139: MatFDColoringSetFunction(fdcoloring,(PetscErrorCode (*)(void))SNESComputeFunction_DMLocal,dmlocalsnes);
140: break;
141: default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"No support for coloring type '%s'",ISColoringTypes[dm->coloringtype]);
142: }
143: PetscObjectSetOptionsPrefix((PetscObject)fdcoloring,((PetscObject)dm)->prefix);
144: MatFDColoringSetFromOptions(fdcoloring);
145: MatFDColoringSetUp(B,coloring,fdcoloring);
146: PetscObjectCompose((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject)fdcoloring);
147: PetscObjectDereference((PetscObject)fdcoloring);
149: /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call
150: * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the
151: * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually
152: * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be
153: * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed.
154: */
155: PetscObjectDereference((PetscObject)dm);
156: }
157: MatFDColoringApply(B,fdcoloring,X,snes);
158: }
159: /* This will be redundant if the user called both, but it's too common to forget. */
160: if (A != B) {
161: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
162: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
163: }
164: return(0);
165: }
167: /*@C
168: DMSNESSetFunctionLocal - set a local residual evaluation function. This function is called with local vector
169: containing the local vector information PLUS ghost point information. It should compute a result for all local
170: elements and DMSNES will automatically accumulate the overlapping values.
172: Logically Collective
174: Input Arguments:
175: + dm - DM to associate callback with
176: . func - local residual evaluation
177: - ctx - optional context for local residual evaluation
179: Level: beginner
181: .seealso: DMSNESSetFunction(), DMDASNESSetJacobianLocal(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
182: @*/
183: PetscErrorCode DMSNESSetFunctionLocal(DM dm,PetscErrorCode (*func)(DM,Vec,Vec,void*),void *ctx)
184: {
186: DMSNES sdm;
187: DMSNES_Local *dmlocalsnes;
191: DMGetDMSNESWrite(dm,&sdm);
192: DMLocalSNESGetContext(dm,sdm,&dmlocalsnes);
194: dmlocalsnes->residuallocal = func;
195: dmlocalsnes->residuallocalctx = ctx;
197: DMSNESSetFunction(dm,SNESComputeFunction_DMLocal,dmlocalsnes);
198: if (!sdm->ops->computejacobian) { /* Call us for the Jacobian too, can be overridden by the user. */
199: DMSNESSetJacobian(dm,SNESComputeJacobian_DMLocal,dmlocalsnes);
200: }
201: return(0);
202: }
204: /*@C
205: DMSNESSetBoundaryLocal - set a local boundary value function. This function is called with local vector
206: containing the local vector information PLUS ghost point information. It should insert values into the local
207: vector that do not come from the global vector, such as essential boundary condition data.
209: Logically Collective
211: Input Arguments:
212: + dm - DM to associate callback with
213: . func - local boundary value evaluation
214: - ctx - optional context for local boundary value evaluation
216: Level: intermediate
218: .seealso: DMSNESSetFunctionLocal(), DMDASNESSetJacobianLocal()
219: @*/
220: PetscErrorCode DMSNESSetBoundaryLocal(DM dm,PetscErrorCode (*func)(DM,Vec,void*),void *ctx)
221: {
223: DMSNES sdm;
224: DMSNES_Local *dmlocalsnes;
228: DMGetDMSNESWrite(dm,&sdm);
229: DMLocalSNESGetContext(dm,sdm,&dmlocalsnes);
231: dmlocalsnes->boundarylocal = func;
232: dmlocalsnes->boundarylocalctx = ctx;
234: return(0);
235: }
237: /*@C
238: DMSNESSetJacobianLocal - set a local Jacobian evaluation function
240: Logically Collective
242: Input Arguments:
243: + dm - DM to associate callback with
244: . func - local Jacobian evaluation
245: - ctx - optional context for local Jacobian evaluation
247: Level: beginner
249: .seealso: DMSNESSetJacobian(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
250: @*/
251: PetscErrorCode DMSNESSetJacobianLocal(DM dm,PetscErrorCode (*func)(DM,Vec,Mat,Mat,void*),void *ctx)
252: {
254: DMSNES sdm;
255: DMSNES_Local *dmlocalsnes;
259: DMGetDMSNESWrite(dm,&sdm);
260: DMLocalSNESGetContext(dm,sdm,&dmlocalsnes);
262: dmlocalsnes->jacobianlocal = func;
263: dmlocalsnes->jacobianlocalctx = ctx;
265: DMSNESSetJacobian(dm,SNESComputeJacobian_DMLocal,dmlocalsnes);
266: return(0);
267: }
269: /*@C
270: DMSNESGetFunctionLocal - get the local residual evaluation function information set with DMSNESSetFunctionLocal.
272: Not Collective
274: Input Arguments:
275: . dm - DM with the associated callback
277: Output Arguments:
278: + func - local residual evaluation
279: - ctx - context for local residual evaluation
281: Level: beginner
283: .seealso: DMSNESSetFunction(), DMSNESSetFunctionLocal(), DMDASNESSetJacobianLocal(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
284: @*/
285: PetscErrorCode DMSNESGetFunctionLocal(DM dm,PetscErrorCode (**func)(DM,Vec,Vec,void*),void **ctx)
286: {
288: DMSNES sdm;
289: DMSNES_Local *dmlocalsnes;
293: DMGetDMSNES(dm,&sdm);
294: DMLocalSNESGetContext(dm,sdm,&dmlocalsnes);
295: if (func) *func = dmlocalsnes->residuallocal;
296: if (ctx) *ctx = dmlocalsnes->residuallocalctx;
297: return(0);
298: }
300: /*@C
301: DMSNESGetBoundaryLocal - get the local boundary value function set with DMSNESSetBoundaryLocal.
303: Not Collective
305: Input Arguments:
306: . dm - DM with the associated callback
308: Output Arguments:
309: + func - local boundary value evaluation
310: - ctx - context for local boundary value evaluation
312: Level: intermediate
314: .seealso: DMSNESSetFunctionLocal(), DMSNESSetBoundaryLocal(), DMDASNESSetJacobianLocal()
315: @*/
316: PetscErrorCode DMSNESGetBoundaryLocal(DM dm,PetscErrorCode (**func)(DM,Vec,void*),void **ctx)
317: {
319: DMSNES sdm;
320: DMSNES_Local *dmlocalsnes;
324: DMGetDMSNES(dm,&sdm);
325: DMLocalSNESGetContext(dm,sdm,&dmlocalsnes);
326: if (func) *func = dmlocalsnes->boundarylocal;
327: if (ctx) *ctx = dmlocalsnes->boundarylocalctx;
328: return(0);
329: }
331: /*@C
332: DMSNESGetJacobianLocal - the local Jacobian evaluation function set with DMSNESSetJacobianLocal.
334: Logically Collective
336: Input Arguments:
337: . dm - DM with the associated callback
339: Output Arguments:
340: + func - local Jacobian evaluation
341: - ctx - context for local Jacobian evaluation
343: Level: beginner
345: .seealso: DMSNESSetJacobianLocal(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
346: @*/
347: PetscErrorCode DMSNESGetJacobianLocal(DM dm,PetscErrorCode (**func)(DM,Vec,Mat,Mat,void*),void **ctx)
348: {
350: DMSNES sdm;
351: DMSNES_Local *dmlocalsnes;
355: DMGetDMSNES(dm,&sdm);
356: DMLocalSNESGetContext(dm,sdm,&dmlocalsnes);
357: if (func) *func = dmlocalsnes->jacobianlocal;
358: if (ctx) *ctx = dmlocalsnes->jacobianlocalctx;
359: return(0);
360: }