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:   PetscFunctionBegin;
 16:   PetscCall(PetscFree(sdm->data));
 17:   sdm->data = NULL;
 18:   PetscFunctionReturn(PETSC_SUCCESS);
 19: }

 21: static PetscErrorCode DMSNESDuplicate_DMLocal(DMSNES oldsdm, DMSNES sdm)
 22: {
 23:   PetscFunctionBegin;
 24:   if (sdm->data != oldsdm->data) {
 25:     PetscCall(PetscFree(sdm->data));
 26:     PetscCall(PetscNew((DMSNES_Local **)&sdm->data));
 27:     if (oldsdm->data) PetscCall(PetscMemcpy(sdm->data, oldsdm->data, sizeof(DMSNES_Local)));
 28:   }
 29:   PetscFunctionReturn(PETSC_SUCCESS);
 30: }

 32: static PetscErrorCode DMLocalSNESGetContext(DM dm, DMSNES sdm, DMSNES_Local **dmlocalsnes)
 33: {
 34:   PetscFunctionBegin;
 35:   *dmlocalsnes = NULL;
 36:   if (!sdm->data) {
 37:     PetscCall(PetscNew((DMSNES_Local **)&sdm->data));

 39:     sdm->ops->destroy   = DMSNESDestroy_DMLocal;
 40:     sdm->ops->duplicate = DMSNESDuplicate_DMLocal;
 41:   }
 42:   *dmlocalsnes = (DMSNES_Local *)sdm->data;
 43:   PetscFunctionReturn(PETSC_SUCCESS);
 44: }

 46: static PetscErrorCode SNESComputeFunction_DMLocal(SNES snes, Vec X, Vec F, void *ctx)
 47: {
 48:   DMSNES_Local *dmlocalsnes = (DMSNES_Local *)ctx;
 49:   DM            dm;
 50:   Vec           Xloc, Floc;
 51:   PetscBool     transform;

 53:   PetscFunctionBegin;
 57:   PetscCall(SNESGetDM(snes, &dm));
 58:   PetscCall(DMGetLocalVector(dm, &Xloc));
 59:   PetscCall(DMGetLocalVector(dm, &Floc));
 60:   PetscCall(VecZeroEntries(Xloc));
 61:   PetscCall(VecZeroEntries(Floc));
 62:   /* Non-conforming routines needs boundary values before G2L */
 63:   if (dmlocalsnes->boundarylocal) PetscCall((*dmlocalsnes->boundarylocal)(dm, Xloc, dmlocalsnes->boundarylocalctx));
 64:   PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, Xloc));
 65:   PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, Xloc));
 66:   /* Need to reset boundary values if we transformed */
 67:   PetscCall(DMHasBasisTransform(dm, &transform));
 68:   if (transform && dmlocalsnes->boundarylocal) PetscCall((*dmlocalsnes->boundarylocal)(dm, Xloc, dmlocalsnes->boundarylocalctx));
 69:   CHKMEMQ;
 70:   PetscCall((*dmlocalsnes->residuallocal)(dm, Xloc, Floc, dmlocalsnes->residuallocalctx));
 71:   CHKMEMQ;
 72:   PetscCall(VecZeroEntries(F));
 73:   PetscCall(DMLocalToGlobalBegin(dm, Floc, ADD_VALUES, F));
 74:   PetscCall(DMLocalToGlobalEnd(dm, Floc, ADD_VALUES, F));
 75:   PetscCall(DMRestoreLocalVector(dm, &Floc));
 76:   PetscCall(DMRestoreLocalVector(dm, &Xloc));
 77:   {
 78:     char        name[PETSC_MAX_PATH_LEN];
 79:     char        oldname[PETSC_MAX_PATH_LEN];
 80:     const char *tmp;
 81:     PetscInt    it;

 83:     PetscCall(SNESGetIterationNumber(snes, &it));
 84:     PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "Solution, Iterate %d", (int)it));
 85:     PetscCall(PetscObjectGetName((PetscObject)X, &tmp));
 86:     PetscCall(PetscStrncpy(oldname, tmp, PETSC_MAX_PATH_LEN - 1));
 87:     PetscCall(PetscObjectSetName((PetscObject)X, name));
 88:     PetscCall(VecViewFromOptions(X, (PetscObject)snes, "-dmsnes_solution_vec_view"));
 89:     PetscCall(PetscObjectSetName((PetscObject)X, oldname));
 90:     PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "Residual, Iterate %d", (int)it));
 91:     PetscCall(PetscObjectSetName((PetscObject)F, name));
 92:     PetscCall(VecViewFromOptions(F, (PetscObject)snes, "-dmsnes_residual_vec_view"));
 93:   }
 94:   PetscFunctionReturn(PETSC_SUCCESS);
 95: }

 97: static PetscErrorCode SNESComputeJacobian_DMLocal(SNES snes, Vec X, Mat A, Mat B, void *ctx)
 98: {
 99:   DMSNES_Local *dmlocalsnes = (DMSNES_Local *)ctx;
100:   DM            dm;
101:   Vec           Xloc;
102:   PetscBool     transform;

104:   PetscFunctionBegin;
105:   PetscCall(SNESGetDM(snes, &dm));
106:   if (dmlocalsnes->jacobianlocal) {
107:     PetscCall(DMGetLocalVector(dm, &Xloc));
108:     PetscCall(VecZeroEntries(Xloc));
109:     /* Non-conforming routines needs boundary values before G2L */
110:     if (dmlocalsnes->boundarylocal) PetscCall((*dmlocalsnes->boundarylocal)(dm, Xloc, dmlocalsnes->boundarylocalctx));
111:     PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, Xloc));
112:     PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, Xloc));
113:     /* Need to reset boundary values if we transformed */
114:     PetscCall(DMHasBasisTransform(dm, &transform));
115:     if (transform && dmlocalsnes->boundarylocal) PetscCall((*dmlocalsnes->boundarylocal)(dm, Xloc, dmlocalsnes->boundarylocalctx));
116:     CHKMEMQ;
117:     PetscCall((*dmlocalsnes->jacobianlocal)(dm, Xloc, A, B, dmlocalsnes->jacobianlocalctx));
118:     CHKMEMQ;
119:     PetscCall(DMRestoreLocalVector(dm, &Xloc));
120:   } else {
121:     MatFDColoring fdcoloring;
122:     PetscCall(PetscObjectQuery((PetscObject)dm, "DMDASNES_FDCOLORING", (PetscObject *)&fdcoloring));
123:     if (!fdcoloring) {
124:       ISColoring coloring;

126:       PetscCall(DMCreateColoring(dm, dm->coloringtype, &coloring));
127:       PetscCall(MatFDColoringCreate(B, coloring, &fdcoloring));
128:       PetscCall(ISColoringDestroy(&coloring));
129:       switch (dm->coloringtype) {
130:       case IS_COLORING_GLOBAL:
131:         PetscCall(MatFDColoringSetFunction(fdcoloring, (PetscErrorCode(*)(void))SNESComputeFunction_DMLocal, dmlocalsnes));
132:         break;
133:       default:
134:         SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "No support for coloring type '%s'", ISColoringTypes[dm->coloringtype]);
135:       }
136:       PetscCall(PetscObjectSetOptionsPrefix((PetscObject)fdcoloring, ((PetscObject)dm)->prefix));
137:       PetscCall(MatFDColoringSetFromOptions(fdcoloring));
138:       PetscCall(MatFDColoringSetUp(B, coloring, fdcoloring));
139:       PetscCall(PetscObjectCompose((PetscObject)dm, "DMDASNES_FDCOLORING", (PetscObject)fdcoloring));
140:       PetscCall(PetscObjectDereference((PetscObject)fdcoloring));

142:       /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call
143:        * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the
144:        * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually
145:        * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be
146:        * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed.
147:        */
148:       PetscCall(PetscObjectDereference((PetscObject)dm));
149:     }
150:     PetscCall(MatFDColoringApply(B, fdcoloring, X, snes));
151:   }
152:   /* This will be redundant if the user called both, but it's too common to forget. */
153:   if (A != B) {
154:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
155:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
156:   }
157:   PetscFunctionReturn(PETSC_SUCCESS);
158: }

160: /*@C
161:   DMSNESSetFunctionLocal - set a local residual evaluation function. This function is called with local vector
162:   containing the local vector information PLUS ghost point information. It should compute a result for all local
163:   elements and `DMSNES` will automatically accumulate the overlapping values.

165:   Logically Collective

167:   Input Parameters:
168: + dm   - `DM` to associate callback with
169: . func - local residual evaluation
170: - ctx  - optional context for local residual evaluation

172:   Calling sequence of `func`:
173: + dm  - `DM` for the function
174: . x   - vector to state at which to evaluate residual
175: . f   - vector to hold the function evaluation
176: - ctx - optional context passed above

178:   Level: advanced

180: .seealso: [](ch_snes), `DMSNESSetFunction()`, `DMDASNESSetJacobianLocal()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`
181: @*/
182: PetscErrorCode DMSNESSetFunctionLocal(DM dm, PetscErrorCode (*func)(DM dm, Vec x, Vec f, void *ctx), void *ctx)
183: {
184:   DMSNES        sdm;
185:   DMSNES_Local *dmlocalsnes;

187:   PetscFunctionBegin;
189:   PetscCall(DMGetDMSNESWrite(dm, &sdm));
190:   PetscCall(DMLocalSNESGetContext(dm, sdm, &dmlocalsnes));

192:   dmlocalsnes->residuallocal    = func;
193:   dmlocalsnes->residuallocalctx = ctx;

195:   PetscCall(DMSNESSetFunction(dm, SNESComputeFunction_DMLocal, dmlocalsnes));
196:   if (!sdm->ops->computejacobian) { /* Call us for the Jacobian too, can be overridden by the user. */
197:     PetscCall(DMSNESSetJacobian(dm, SNESComputeJacobian_DMLocal, dmlocalsnes));
198:   }
199:   PetscFunctionReturn(PETSC_SUCCESS);
200: }

202: /*@C
203:   DMSNESSetBoundaryLocal - set a function to insert, for example, essential boundary conditions into a ghosted solution vector

205:   Logically Collective

207:   Input Parameters:
208: + dm   - `DM` to associate callback with
209: . func - local boundary value evaluation
210: - ctx  - optional context for local boundary value evaluation

212:   Calling sequence of `func`:
213: + dm  - the `DM` context
214: . X   - ghosted solution vector, appropriate locations (such as essential boundary condition nodes) should be filled
215: - ctx - option context passed in `DMSNESSetBoundaryLocal()`

217:   Level: advanced

219: .seealso: [](ch_snes), `DMSNESSetFunctionLocal()`, `DMDASNESSetJacobianLocal()`
220: @*/
221: PetscErrorCode DMSNESSetBoundaryLocal(DM dm, PetscErrorCode (*func)(DM dm, Vec X, void *ctx), void *ctx)
222: {
223:   DMSNES        sdm;
224:   DMSNES_Local *dmlocalsnes;

226:   PetscFunctionBegin;
228:   PetscCall(DMGetDMSNESWrite(dm, &sdm));
229:   PetscCall(DMLocalSNESGetContext(dm, sdm, &dmlocalsnes));

231:   dmlocalsnes->boundarylocal    = func;
232:   dmlocalsnes->boundarylocalctx = ctx;

234:   PetscFunctionReturn(PETSC_SUCCESS);
235: }

237: /*@C
238:   DMSNESSetJacobianLocal - set a local Jacobian evaluation function

240:   Logically Collective

242:   Input Parameters:
243: + dm   - `DM` to associate callback with
244: . func - local Jacobian evaluation
245: - ctx  - optional context for local Jacobian evaluation

247:   Calling sequence of `func`:
248: + dm  - the `DM` context
249: . X   - current solution vector (ghosted or not?)
250: . J   - the Jacobian
251: . Jp  - approximate Jacobian used to compute the preconditioner, often `J`
252: - ctx - a user provided context

254:   Level: advanced

256: .seealso: [](ch_snes), `DMSNESSetJacobian()`, `DMDASNESSetJacobian()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`
257: @*/
258: PetscErrorCode DMSNESSetJacobianLocal(DM dm, PetscErrorCode (*func)(DM dm, Vec X, Mat J, Mat Jp, void *ctx), void *ctx)
259: {
260:   DMSNES        sdm;
261:   DMSNES_Local *dmlocalsnes;

263:   PetscFunctionBegin;
265:   PetscCall(DMGetDMSNESWrite(dm, &sdm));
266:   PetscCall(DMLocalSNESGetContext(dm, sdm, &dmlocalsnes));

268:   dmlocalsnes->jacobianlocal    = func;
269:   dmlocalsnes->jacobianlocalctx = ctx;

271:   PetscCall(DMSNESSetJacobian(dm, SNESComputeJacobian_DMLocal, dmlocalsnes));
272:   PetscFunctionReturn(PETSC_SUCCESS);
273: }

275: /*@C
276:   DMSNESGetFunctionLocal - get the local residual evaluation function information set with `DMSNESSetFunctionLocal()`.

278:   Not Collective

280:   Input Parameter:
281: . dm - `DM` with the associated callback

283:   Output Parameters:
284: + func - local residual evaluation
285: - ctx  - context for local residual evaluation

287:   Level: beginner

289: .seealso: [](ch_snes), `DMSNESSetFunction()`, `DMSNESSetFunctionLocal()`, `DMDASNESSetJacobianLocal()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`
290: @*/
291: PetscErrorCode DMSNESGetFunctionLocal(DM dm, PetscErrorCode (**func)(DM, Vec, Vec, void *), void **ctx)
292: {
293:   DMSNES        sdm;
294:   DMSNES_Local *dmlocalsnes;

296:   PetscFunctionBegin;
298:   PetscCall(DMGetDMSNES(dm, &sdm));
299:   PetscCall(DMLocalSNESGetContext(dm, sdm, &dmlocalsnes));
300:   if (func) *func = dmlocalsnes->residuallocal;
301:   if (ctx) *ctx = dmlocalsnes->residuallocalctx;
302:   PetscFunctionReturn(PETSC_SUCCESS);
303: }

305: /*@C
306:   DMSNESGetBoundaryLocal - get the local boundary value function set with `DMSNESSetBoundaryLocal()`.

308:   Not Collective

310:   Input Parameter:
311: . dm - `DM` with the associated callback

313:   Output Parameters:
314: + func - local boundary value evaluation
315: - ctx  - context for local boundary value evaluation

317:   Level: intermediate

319: .seealso: [](ch_snes), `DMSNESSetFunctionLocal()`, `DMSNESSetBoundaryLocal()`, `DMDASNESSetJacobianLocal()`
320: @*/
321: PetscErrorCode DMSNESGetBoundaryLocal(DM dm, PetscErrorCode (**func)(DM, Vec, void *), void **ctx)
322: {
323:   DMSNES        sdm;
324:   DMSNES_Local *dmlocalsnes;

326:   PetscFunctionBegin;
328:   PetscCall(DMGetDMSNES(dm, &sdm));
329:   PetscCall(DMLocalSNESGetContext(dm, sdm, &dmlocalsnes));
330:   if (func) *func = dmlocalsnes->boundarylocal;
331:   if (ctx) *ctx = dmlocalsnes->boundarylocalctx;
332:   PetscFunctionReturn(PETSC_SUCCESS);
333: }

335: /*@C
336:   DMSNESGetJacobianLocal - the local Jacobian evaluation function set with `DMSNESSetJacobianLocal()`.

338:   Logically Collective

340:   Input Parameter:
341: . dm - `DM` with the associated callback

343:   Output Parameters:
344: + func - local Jacobian evaluation
345: - ctx  - context for local Jacobian evaluation

347:   Level: beginner

349: .seealso: [](ch_snes), `DMSNESSetJacobianLocal()`, `DMDASNESSetJacobian()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`
350: @*/
351: PetscErrorCode DMSNESGetJacobianLocal(DM dm, PetscErrorCode (**func)(DM, Vec, Mat, Mat, void *), void **ctx)
352: {
353:   DMSNES        sdm;
354:   DMSNES_Local *dmlocalsnes;

356:   PetscFunctionBegin;
358:   PetscCall(DMGetDMSNES(dm, &sdm));
359:   PetscCall(DMLocalSNESGetContext(dm, sdm, &dmlocalsnes));
360:   if (func) *func = dmlocalsnes->jacobianlocal;
361:   if (ctx) *ctx = dmlocalsnes->jacobianlocalctx;
362:   PetscFunctionReturn(PETSC_SUCCESS);
363: }