Actual source code: dmlocalsnes.c
petsc-3.11.4 2019-09-28
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: {
55: DM dm;
56: DMSNES_Local *dmlocalsnes = (DMSNES_Local*)ctx;
57: Vec Xloc,Floc;
63: SNESGetDM(snes,&dm);
64: DMGetLocalVector(dm,&Xloc);
65: DMGetLocalVector(dm,&Floc);
66: VecZeroEntries(Xloc);
67: if (dmlocalsnes->boundarylocal) {(*dmlocalsnes->boundarylocal)(dm,Xloc,dmlocalsnes->boundarylocalctx);}
68: VecZeroEntries(Floc);
69: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
70: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
71: CHKMEMQ;
72: (*dmlocalsnes->residuallocal)(dm,Xloc,Floc,dmlocalsnes->residuallocalctx);
73: CHKMEMQ;
74: VecZeroEntries(F);
75: DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
76: DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
77: DMRestoreLocalVector(dm,&Floc);
78: DMRestoreLocalVector(dm,&Xloc);
79: {
80: char name[PETSC_MAX_PATH_LEN];
81: char oldname[PETSC_MAX_PATH_LEN];
82: const char *tmp;
83: PetscInt it;
85: SNESGetIterationNumber(snes, &it);
86: PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "Solution, Iterate %d", (int) it);
87: PetscObjectGetName((PetscObject) X, &tmp);
88: PetscStrncpy(oldname, tmp, PETSC_MAX_PATH_LEN-1);
89: PetscObjectSetName((PetscObject) X, name);
90: VecViewFromOptions(X, (PetscObject) snes, "-dmsnes_solution_vec_view");
91: PetscObjectSetName((PetscObject) X, oldname);
92: PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "Residual, Iterate %d", (int) it);
93: PetscObjectSetName((PetscObject) F, name);
94: VecViewFromOptions(F, (PetscObject) snes, "-dmsnes_residual_vec_view");
95: }
96: return(0);
97: }
99: static PetscErrorCode SNESComputeJacobian_DMLocal(SNES snes,Vec X,Mat A,Mat B,void *ctx)
100: {
102: DM dm;
103: DMSNES_Local *dmlocalsnes = (DMSNES_Local*)ctx;
104: Vec Xloc;
107: SNESGetDM(snes,&dm);
108: if (dmlocalsnes->jacobianlocal) {
109: DMGetLocalVector(dm,&Xloc);
110: VecZeroEntries(Xloc);
111: if (dmlocalsnes->boundarylocal) {(*dmlocalsnes->boundarylocal)(dm,Xloc,dmlocalsnes->boundarylocalctx);}
112: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
113: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
114: CHKMEMQ;
115: (*dmlocalsnes->jacobianlocal)(dm,Xloc,A,B,dmlocalsnes->jacobianlocalctx);
116: CHKMEMQ;
117: DMRestoreLocalVector(dm,&Xloc);
118: } else {
119: MatFDColoring fdcoloring;
120: PetscObjectQuery((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject*)&fdcoloring);
121: if (!fdcoloring) {
122: ISColoring coloring;
124: DMCreateColoring(dm,dm->coloringtype,&coloring);
125: MatFDColoringCreate(B,coloring,&fdcoloring);
126: ISColoringDestroy(&coloring);
127: switch (dm->coloringtype) {
128: case IS_COLORING_GLOBAL:
129: MatFDColoringSetFunction(fdcoloring,(PetscErrorCode (*)(void))SNESComputeFunction_DMLocal,dmlocalsnes);
130: break;
131: default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"No support for coloring type '%s'",ISColoringTypes[dm->coloringtype]);
132: }
133: PetscObjectSetOptionsPrefix((PetscObject)fdcoloring,((PetscObject)dm)->prefix);
134: MatFDColoringSetFromOptions(fdcoloring);
135: MatFDColoringSetUp(B,coloring,fdcoloring);
136: PetscObjectCompose((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject)fdcoloring);
137: PetscObjectDereference((PetscObject)fdcoloring);
139: /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call
140: * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the
141: * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually
142: * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be
143: * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed.
144: */
145: PetscObjectDereference((PetscObject)dm);
146: }
147: MatFDColoringApply(B,fdcoloring,X,snes);
148: }
149: /* This will be redundant if the user called both, but it's too common to forget. */
150: if (A != B) {
151: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
152: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
153: }
154: return(0);
155: }
157: /*@C
158: DMSNESSetFunctionLocal - set a local residual evaluation function. This function is called with local vector
159: containing the local vector information PLUS ghost point information. It should compute a result for all local
160: elements and DMSNES will automatically accumulate the overlapping values.
162: Logically Collective
164: Input Arguments:
165: + dm - DM to associate callback with
166: . func - local residual evaluation
167: - ctx - optional context for local residual evaluation
169: Level: beginner
171: .seealso: DMSNESSetFunction(), DMDASNESSetJacobianLocal(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
172: @*/
173: PetscErrorCode DMSNESSetFunctionLocal(DM dm,PetscErrorCode (*func)(DM,Vec,Vec,void*),void *ctx)
174: {
176: DMSNES sdm;
177: DMSNES_Local *dmlocalsnes;
181: DMGetDMSNESWrite(dm,&sdm);
182: DMLocalSNESGetContext(dm,sdm,&dmlocalsnes);
184: dmlocalsnes->residuallocal = func;
185: dmlocalsnes->residuallocalctx = ctx;
187: DMSNESSetFunction(dm,SNESComputeFunction_DMLocal,dmlocalsnes);
188: if (!sdm->ops->computejacobian) { /* Call us for the Jacobian too, can be overridden by the user. */
189: DMSNESSetJacobian(dm,SNESComputeJacobian_DMLocal,dmlocalsnes);
190: }
191: return(0);
192: }
194: /*@C
195: DMSNESSetBoundaryLocal - set a local boundary value function. This function is called with local vector
196: containing the local vector information PLUS ghost point information. It should insert values into the local
197: vector that do not come from the global vector, such as essential boundary condition data.
199: Logically Collective
201: Input Arguments:
202: + dm - DM to associate callback with
203: . func - local boundary value evaluation
204: - ctx - optional context for local boundary value evaluation
206: Level: intermediate
208: .seealso: DMSNESSetFunctionLocal(), DMDASNESSetJacobianLocal()
209: @*/
210: PetscErrorCode DMSNESSetBoundaryLocal(DM dm,PetscErrorCode (*func)(DM,Vec,void*),void *ctx)
211: {
213: DMSNES sdm;
214: DMSNES_Local *dmlocalsnes;
218: DMGetDMSNESWrite(dm,&sdm);
219: DMLocalSNESGetContext(dm,sdm,&dmlocalsnes);
221: dmlocalsnes->boundarylocal = func;
222: dmlocalsnes->boundarylocalctx = ctx;
224: return(0);
225: }
227: /*@C
228: DMSNESSetJacobianLocal - set a local Jacobian evaluation function
230: Logically Collective
232: Input Arguments:
233: + dm - DM to associate callback with
234: . func - local Jacobian evaluation
235: - ctx - optional context for local Jacobian evaluation
237: Level: beginner
239: .seealso: DMSNESSetJacobian(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
240: @*/
241: PetscErrorCode DMSNESSetJacobianLocal(DM dm,PetscErrorCode (*func)(DM,Vec,Mat,Mat,void*),void *ctx)
242: {
244: DMSNES sdm;
245: DMSNES_Local *dmlocalsnes;
249: DMGetDMSNESWrite(dm,&sdm);
250: DMLocalSNESGetContext(dm,sdm,&dmlocalsnes);
252: dmlocalsnes->jacobianlocal = func;
253: dmlocalsnes->jacobianlocalctx = ctx;
255: DMSNESSetJacobian(dm,SNESComputeJacobian_DMLocal,dmlocalsnes);
256: return(0);
257: }