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

     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:
    427:     The user must use
    428:     SNESSetFunction(snes,NULL,SNESPicardComputeFunction,&user);
    429:     in their code before calling this routine.

    432:    Level: beginner

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

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

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

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