Actual source code: dmlocalsnes.c

petsc-3.12.0 2019-09-29
Report Typos and Errors
  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: }