Actual source code: dmdasnes.c

petsc-3.8.2 2017-11-09
Report Typos and Errors
  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: }


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

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

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

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

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

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


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

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

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

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

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

207: /*@C
208:    DMDASNESSetFunctionLocal - set a local residual evaluation function

210:    Logically Collective

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

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

225:    Level: beginner

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

237:   DMGetDMSNESWrite(dm,&sdm);
238:   DMDASNESGetContext(dm,sdm,&dmdasnes);

240:   dmdasnes->residuallocalimode = imode;
241:   dmdasnes->residuallocal      = func;
242:   dmdasnes->residuallocalctx   = ctx;

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

251: /*@C
252:    DMDASNESSetJacobianLocal - set a local Jacobian evaluation function

254:    Logically Collective

256:    Input Arguments:
257: +  dm - DM to associate callback with
258: .  func - local Jacobian evaluation
259: -  ctx - optional context for local Jacobian evaluation

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

269:    Level: beginner

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

281:   DMGetDMSNESWrite(dm,&sdm);
282:   DMDASNESGetContext(dm,sdm,&dmdasnes);

284:   dmdasnes->jacobianlocal    = func;
285:   dmdasnes->jacobianlocalctx = ctx;

287:   DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes);
288:   return(0);
289: }


292: /*@C
293:    DMDASNESSetObjectiveLocal - set a local residual evaluation function

295:    Logically Collective

297:    Input Arguments:
298: +  dm - DM to associate callback with
299: .  func - local objective evaluation
300: -  ctx - optional context for local residual evaluation

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

308:    Level: beginner

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

320:   DMGetDMSNESWrite(dm,&sdm);
321:   DMDASNESGetContext(dm,sdm,&dmdasnes);

323:   dmdasnes->objectivelocal    = func;
324:   dmdasnes->objectivelocalctx = ctx;

326:   DMSNESSetObjective(dm,SNESComputeObjective_DMDA,dmdasnes);
327:   return(0);
328: }

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

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

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

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

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

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

412:    Logically Collective

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

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

426:    Notes:  The user must use
427:     extern PetscErrorCode  SNESPicardComputeFunction(SNES,Vec,Vec,void*);
428:     extern PetscErrorCode  SNESPicardComputeJacobian(SNES,Vec,Mat,Mat,MatStructure*,void*);
429:     SNESSetFunction(snes,NULL,SNESPicardComputeFunction,&user);
430:     in their code before calling this routine.


433:    Level: beginner

435: .seealso: DMSNESSetFunction(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
436: @*/
437: PetscErrorCode DMDASNESSetPicardLocal(DM dm,InsertMode imode,PetscErrorCode (*func)(DMDALocalInfo*,void*,void*,void*),
438:                                       PetscErrorCode (*jac)(DMDALocalInfo*,void*,Mat,Mat,void*),void *ctx)
439: {
441:   DMSNES         sdm;
442:   DMSNES_DA      *dmdasnes;

446:   DMGetDMSNESWrite(dm,&sdm);
447:   DMDASNESGetContext(dm,sdm,&dmdasnes);

449:   dmdasnes->residuallocalimode = imode;
450:   dmdasnes->rhsplocal          = func;
451:   dmdasnes->jacobianplocal     = jac;
452:   dmdasnes->picardlocalctx     = ctx;

454:   DMSNESSetPicard(dm,SNESComputePicard_DMDA,SNESComputePicardJacobian_DMDA,dmdasnes);
455:   return(0);
456: }