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: }