Actual source code: dmdasnes.c

  1: #include <petscdmda.h>
  2: #include <petsc/private/dmimpl.h>
  3: #include <petsc/private/snesimpl.h>

  5: /* This structure holds the user-provided DMDA callbacks */
  6: typedef struct {
  7:   PetscErrorCode (*residuallocal)(DMDALocalInfo*,void*,void*,void*);
  8:   PetscErrorCode (*jacobianlocal)(DMDALocalInfo*,void*,Mat,Mat,void*);
  9:   PetscErrorCode (*objectivelocal)(DMDALocalInfo*,void*,PetscReal*,void*);
 10:   void       *residuallocalctx;
 11:   void       *jacobianlocalctx;
 12:   void       *objectivelocalctx;
 13:   InsertMode residuallocalimode;

 15:   /*   For Picard iteration defined locally */
 16:   PetscErrorCode (*rhsplocal)(DMDALocalInfo*,void*,void*,void*);
 17:   PetscErrorCode (*jacobianplocal)(DMDALocalInfo*,void*,Mat,Mat,void*);
 18:   void *picardlocalctx;
 19: } DMSNES_DA;

 21: static PetscErrorCode DMSNESDestroy_DMDA(DMSNES sdm)
 22: {

 26:   PetscFree(sdm->data);
 27:   return(0);
 28: }

 30: static PetscErrorCode DMSNESDuplicate_DMDA(DMSNES oldsdm,DMSNES sdm)
 31: {

 35:   PetscNewLog(sdm,(DMSNES_DA**)&sdm->data);
 36:   if (oldsdm->data) {
 37:     PetscMemcpy(sdm->data,oldsdm->data,sizeof(DMSNES_DA));
 38:   }
 39:   return(0);
 40: }

 42: static PetscErrorCode DMDASNESGetContext(DM dm,DMSNES sdm,DMSNES_DA  **dmdasnes)
 43: {

 47:   *dmdasnes = NULL;
 48:   if (!sdm->data) {
 49:     PetscNewLog(dm,(DMSNES_DA**)&sdm->data);
 50:     sdm->ops->destroy   = DMSNESDestroy_DMDA;
 51:     sdm->ops->duplicate = DMSNESDuplicate_DMDA;
 52:   }
 53:   *dmdasnes = (DMSNES_DA*)sdm->data;
 54:   return(0);
 55: }

 57: static PetscErrorCode SNESComputeFunction_DMDA(SNES snes,Vec X,Vec F,void *ctx)
 58: {
 60:   DM             dm;
 61:   DMSNES_DA      *dmdasnes = (DMSNES_DA*)ctx;
 62:   DMDALocalInfo  info;
 63:   Vec            Xloc;
 64:   void           *x,*f;

 70:   if (!dmdasnes->residuallocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
 71:   SNESGetDM(snes,&dm);
 72:   DMGetLocalVector(dm,&Xloc);
 73:   DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
 74:   DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
 75:   DMDAGetLocalInfo(dm,&info);
 76:   DMDAVecGetArray(dm,Xloc,&x);
 77:   switch (dmdasnes->residuallocalimode) {
 78:   case INSERT_VALUES: {
 79:     DMDAVecGetArray(dm,F,&f);
 80:     PetscLogEventBegin(SNES_FunctionEval,snes,X,F,0);
 81:     CHKMEMQ;
 82:     (*dmdasnes->residuallocal)(&info,x,f,dmdasnes->residuallocalctx);
 83:     CHKMEMQ;
 84:     PetscLogEventEnd(SNES_FunctionEval,snes,X,F,0);
 85:     DMDAVecRestoreArray(dm,F,&f);
 86:   } break;
 87:   case ADD_VALUES: {
 88:     Vec Floc;
 89:     DMGetLocalVector(dm,&Floc);
 90:     VecZeroEntries(Floc);
 91:     DMDAVecGetArray(dm,Floc,&f);
 92:     PetscLogEventBegin(SNES_FunctionEval,snes,X,F,0);
 93:     CHKMEMQ;
 94:     (*dmdasnes->residuallocal)(&info,x,f,dmdasnes->residuallocalctx);
 95:     CHKMEMQ;
 96:     PetscLogEventEnd(SNES_FunctionEval,snes,X,F,0);
 97:     DMDAVecRestoreArray(dm,Floc,&f);
 98:     VecZeroEntries(F);
 99:     DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
100:     DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
101:     DMRestoreLocalVector(dm,&Floc);
102:   } break;
103:   default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdasnes->residuallocalimode);
104:   }
105:   DMDAVecRestoreArray(dm,Xloc,&x);
106:   DMRestoreLocalVector(dm,&Xloc);
107:   if (snes->domainerror) {
108:     VecSetInf(F);
109:   }
110:   return(0);
111: }

113: static PetscErrorCode SNESComputeObjective_DMDA(SNES snes,Vec X,PetscReal *ob,void *ctx)
114: {
116:   DM             dm;
117:   DMSNES_DA      *dmdasnes = (DMSNES_DA*)ctx;
118:   DMDALocalInfo  info;
119:   Vec            Xloc;
120:   void           *x;

126:   if (!dmdasnes->objectivelocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
127:   SNESGetDM(snes,&dm);
128:   DMGetLocalVector(dm,&Xloc);
129:   DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
130:   DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
131:   DMDAGetLocalInfo(dm,&info);
132:   DMDAVecGetArray(dm,Xloc,&x);
133:   CHKMEMQ;
134:   (*dmdasnes->objectivelocal)(&info,x,ob,dmdasnes->objectivelocalctx);
135:   CHKMEMQ;
136:   DMDAVecRestoreArray(dm,Xloc,&x);
137:   DMRestoreLocalVector(dm,&Xloc);
138:   return(0);
139: }

141: /* Routine is called by example, hence must be labeled PETSC_EXTERN */
142: PETSC_EXTERN PetscErrorCode SNESComputeJacobian_DMDA(SNES snes,Vec X,Mat A,Mat B,void *ctx)
143: {
145:   DM             dm;
146:   DMSNES_DA      *dmdasnes = (DMSNES_DA*)ctx;
147:   DMDALocalInfo  info;
148:   Vec            Xloc;
149:   void           *x;

152:   if (!dmdasnes->residuallocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
153:   SNESGetDM(snes,&dm);

155:   if (dmdasnes->jacobianlocal) {
156:     DMGetLocalVector(dm,&Xloc);
157:     DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
158:     DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
159:     DMDAGetLocalInfo(dm,&info);
160:     DMDAVecGetArray(dm,Xloc,&x);
161:     CHKMEMQ;
162:     (*dmdasnes->jacobianlocal)(&info,x,A,B,dmdasnes->jacobianlocalctx);
163:     CHKMEMQ;
164:     DMDAVecRestoreArray(dm,Xloc,&x);
165:     DMRestoreLocalVector(dm,&Xloc);
166:   } else {
167:     MatFDColoring fdcoloring;
168:     PetscObjectQuery((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject*)&fdcoloring);
169:     if (!fdcoloring) {
170:       ISColoring coloring;

172:       DMCreateColoring(dm,dm->coloringtype,&coloring);
173:       MatFDColoringCreate(B,coloring,&fdcoloring);
174:       switch (dm->coloringtype) {
175:       case IS_COLORING_GLOBAL:
176:         MatFDColoringSetFunction(fdcoloring,(PetscErrorCode (*)(void))SNESComputeFunction_DMDA,dmdasnes);
177:         break;
178:       default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"No support for coloring type '%s'",ISColoringTypes[dm->coloringtype]);
179:       }
180:       PetscObjectSetOptionsPrefix((PetscObject)fdcoloring,((PetscObject)dm)->prefix);
181:       MatFDColoringSetFromOptions(fdcoloring);
182:       MatFDColoringSetUp(B,coloring,fdcoloring);
183:       ISColoringDestroy(&coloring);
184:       PetscObjectCompose((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject)fdcoloring);
185:       PetscObjectDereference((PetscObject)fdcoloring);

187:       /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call
188:        * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the
189:        * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually
190:        * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be
191:        * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed.
192:        */
193:       PetscObjectDereference((PetscObject)dm);
194:     }
195:     MatFDColoringApply(B,fdcoloring,X,snes);
196:   }
197:   /* This will be redundant if the user called both, but it's too common to forget. */
198:   if (A != B) {
199:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
200:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
201:   }
202:   return(0);
203: }

205: /*@C
206:    DMDASNESSetFunctionLocal - set a local residual evaluation function

208:    Logically Collective

210:    Input Parameters:
211: +  dm - DM to associate callback with
212: .  imode - INSERT_VALUES if local function computes owned part, ADD_VALUES if it contributes to ghosted part
213: .  func - local residual evaluation
214: -  ctx - optional context for local residual evaluation

216:    Calling sequence:
217:    For PetscErrorCode (*func)(DMDALocalInfo *info,void *x, void *f, void *ctx),
218: +  info - DMDALocalInfo defining the subdomain to evaluate the residual on
219: .  x - dimensional pointer to state at which to evaluate residual (e.g. PetscScalar *x or **x or ***x)
220: .  f - dimensional pointer to residual, write the residual here (e.g. PetscScalar *f or **f or ***f)
221: -  ctx - optional context passed above

223:    Level: beginner

225: .seealso: DMDASNESSetJacobianLocal(), DMSNESSetFunction(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
226: @*/
227: PetscErrorCode DMDASNESSetFunctionLocal(DM dm,InsertMode imode,PetscErrorCode (*func)(DMDALocalInfo*,void*,void*,void*),void *ctx)
228: {
230:   DMSNES         sdm;
231:   DMSNES_DA      *dmdasnes;

235:   DMGetDMSNESWrite(dm,&sdm);
236:   DMDASNESGetContext(dm,sdm,&dmdasnes);

238:   dmdasnes->residuallocalimode = imode;
239:   dmdasnes->residuallocal      = func;
240:   dmdasnes->residuallocalctx   = ctx;

242:   DMSNESSetFunction(dm,SNESComputeFunction_DMDA,dmdasnes);
243:   if (!sdm->ops->computejacobian) {  /* Call us for the Jacobian too, can be overridden by the user. */
244:     DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes);
245:   }
246:   return(0);
247: }

249: /*@C
250:    DMDASNESSetJacobianLocal - set a local Jacobian evaluation function

252:    Logically Collective

254:    Input Parameters:
255: +  dm - DM to associate callback with
256: .  func - local Jacobian evaluation
257: -  ctx - optional context for local Jacobian evaluation

259:    Calling sequence:
260:    For PetscErrorCode (*func)(DMDALocalInfo *info,void *x,Mat J,Mat M,void *ctx),
261: +  info - DMDALocalInfo defining the subdomain to evaluate the Jacobian at
262: .  x - dimensional pointer to state at which to evaluate Jacobian (e.g. PetscScalar *x or **x or ***x)
263: .  J - Mat object for the Jacobian
264: .  M - Mat object for the Jacobian preconditioner matrix
265: -  ctx - optional context passed above

267:    Level: beginner

269: .seealso: DMDASNESSetFunctionLocal(), DMSNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
270: @*/
271: PetscErrorCode DMDASNESSetJacobianLocal(DM dm,PetscErrorCode (*func)(DMDALocalInfo*,void*,Mat,Mat,void*),void *ctx)
272: {
274:   DMSNES         sdm;
275:   DMSNES_DA      *dmdasnes;

279:   DMGetDMSNESWrite(dm,&sdm);
280:   DMDASNESGetContext(dm,sdm,&dmdasnes);

282:   dmdasnes->jacobianlocal    = func;
283:   dmdasnes->jacobianlocalctx = ctx;

285:   DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes);
286:   return(0);
287: }

289: /*@C
290:    DMDASNESSetObjectiveLocal - set a local residual evaluation function

292:    Logically Collective

294:    Input Parameters:
295: +  dm - DM to associate callback with
296: .  func - local objective evaluation
297: -  ctx - optional context for local residual evaluation

299:    Calling sequence for func:
300: +  info - DMDALocalInfo defining the subdomain to evaluate the residual on
301: .  x - dimensional pointer to state at which to evaluate residual
302: .  ob - eventual objective value
303: -  ctx - optional context passed above

305:    Level: beginner

307: .seealso: DMSNESSetFunction(), DMDASNESSetJacobianLocal(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
308: @*/
309: PetscErrorCode DMDASNESSetObjectiveLocal(DM dm,DMDASNESObjective func,void *ctx)
310: {
312:   DMSNES         sdm;
313:   DMSNES_DA      *dmdasnes;

317:   DMGetDMSNESWrite(dm,&sdm);
318:   DMDASNESGetContext(dm,sdm,&dmdasnes);

320:   dmdasnes->objectivelocal    = func;
321:   dmdasnes->objectivelocalctx = ctx;

323:   DMSNESSetObjective(dm,SNESComputeObjective_DMDA,dmdasnes);
324:   return(0);
325: }

327: static PetscErrorCode SNESComputePicard_DMDA(SNES snes,Vec X,Vec F,void *ctx)
328: {
330:   DM             dm;
331:   DMSNES_DA      *dmdasnes = (DMSNES_DA*)ctx;
332:   DMDALocalInfo  info;
333:   Vec            Xloc;
334:   void           *x,*f;

340:   if (!dmdasnes->rhsplocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
341:   SNESGetDM(snes,&dm);
342:   DMGetLocalVector(dm,&Xloc);
343:   DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
344:   DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
345:   DMDAGetLocalInfo(dm,&info);
346:   DMDAVecGetArray(dm,Xloc,&x);
347:   switch (dmdasnes->residuallocalimode) {
348:   case INSERT_VALUES: {
349:     DMDAVecGetArray(dm,F,&f);
350:     CHKMEMQ;
351:     (*dmdasnes->rhsplocal)(&info,x,f,dmdasnes->picardlocalctx);
352:     CHKMEMQ;
353:     DMDAVecRestoreArray(dm,F,&f);
354:   } break;
355:   case ADD_VALUES: {
356:     Vec Floc;
357:     DMGetLocalVector(dm,&Floc);
358:     VecZeroEntries(Floc);
359:     DMDAVecGetArray(dm,Floc,&f);
360:     CHKMEMQ;
361:     (*dmdasnes->rhsplocal)(&info,x,f,dmdasnes->picardlocalctx);
362:     CHKMEMQ;
363:     DMDAVecRestoreArray(dm,Floc,&f);
364:     VecZeroEntries(F);
365:     DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
366:     DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
367:     DMRestoreLocalVector(dm,&Floc);
368:   } break;
369:   default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdasnes->residuallocalimode);
370:   }
371:   DMDAVecRestoreArray(dm,Xloc,&x);
372:   DMRestoreLocalVector(dm,&Xloc);
373:   return(0);
374: }

376: static PetscErrorCode SNESComputePicardJacobian_DMDA(SNES snes,Vec X,Mat A,Mat B,void *ctx)
377: {
379:   DM             dm;
380:   DMSNES_DA      *dmdasnes = (DMSNES_DA*)ctx;
381:   DMDALocalInfo  info;
382:   Vec            Xloc;
383:   void           *x;

386:   if (!dmdasnes->jacobianplocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
387:   SNESGetDM(snes,&dm);

389:   DMGetLocalVector(dm,&Xloc);
390:   DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
391:   DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
392:   DMDAGetLocalInfo(dm,&info);
393:   DMDAVecGetArray(dm,Xloc,&x);
394:   CHKMEMQ;
395:   (*dmdasnes->jacobianplocal)(&info,x,A,B,dmdasnes->picardlocalctx);
396:   CHKMEMQ;
397:   DMDAVecRestoreArray(dm,Xloc,&x);
398:   DMRestoreLocalVector(dm,&Xloc);
399:   if (A != B) {
400:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
401:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
402:   }
403:   return(0);
404: }

406: /*@C
407:    DMDASNESSetPicardLocal - set a local right hand side and matrix evaluation function for Picard iteration

409:    Logically Collective

411:    Input Parameters:
412: +  dm - DM to associate callback with
413: .  imode - INSERT_VALUES if local function computes owned part, ADD_VALUES if it contributes to ghosted part
414: .  func - local residual evaluation
415: -  ctx - optional context for local residual evaluation

417:    Calling sequence for func:
418: +  info - DMDALocalInfo defining the subdomain to evaluate the residual on
419: .  x - dimensional pointer to state at which to evaluate residual
420: .  f - dimensional pointer to residual, write the residual here
421: -  ctx - optional context passed above

423:    Notes:
424:     The user must use
425:     SNESSetFunction(snes,NULL,SNESPicardComputeFunction,&user);
426:     in their code before calling this routine.

428:    Level: beginner

430: .seealso: DMSNESSetFunction(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
431: @*/
432: PetscErrorCode DMDASNESSetPicardLocal(DM dm,InsertMode imode,PetscErrorCode (*func)(DMDALocalInfo*,void*,void*,void*),
433:                                       PetscErrorCode (*jac)(DMDALocalInfo*,void*,Mat,Mat,void*),void *ctx)
434: {
436:   DMSNES         sdm;
437:   DMSNES_DA      *dmdasnes;

441:   DMGetDMSNESWrite(dm,&sdm);
442:   DMDASNESGetContext(dm,sdm,&dmdasnes);

444:   dmdasnes->residuallocalimode = imode;
445:   dmdasnes->rhsplocal          = func;
446:   dmdasnes->jacobianplocal     = jac;
447:   dmdasnes->picardlocalctx     = ctx;

449:   DMSNESSetPicard(dm,SNESComputePicard_DMDA,SNESComputePicardJacobian_DMDA,dmdasnes);
450:   DMSNESSetMFFunction(dm,SNESComputeFunction_DMDA,dmdasnes);
451:   return(0);
452: }