Actual source code: ex69.c
petsc-3.12.4 2020-02-04
Report Typos and Errors

      2: static char help[] = "Tests recovery from domain errors in MatMult() and PCApply()\n\n";

      4: /*
      5:       See src/ksp/ksp/examples/tutorials/ex19.c from which this was copied
      6: */

      8:  #include <petscsnes.h>
      9:  #include <petscdm.h>
     10:  #include <petscdmda.h>

     12: /*
     13:    User-defined routines and data structures
     14: */
     15: typedef struct {
     16:   PetscScalar u,v,omega,temp;
     17: } Field;

     19: PetscErrorCode FormFunctionLocal(DMDALocalInfo*,Field**,Field**,void*);

     21: typedef struct {
     22:   PetscReal   lidvelocity,prandtl,grashof;  /* physical parameters */
     23:   PetscBool   draw_contours;                /* flag - 1 indicates drawing contours */
     24:   PetscBool   errorindomain;
     25:   PetscBool   errorindomainmf;
     26:   SNES        snes;
     27: } AppCtx;

     29: typedef struct {
     30:   Mat Jmf;
     31: } MatShellCtx;

     33: extern PetscErrorCode FormInitialGuess(AppCtx*,DM,Vec);
     34: extern PetscErrorCode MatMult_MyShell(Mat,Vec,Vec);
     35: extern PetscErrorCode MatAssemblyEnd_MyShell(Mat,MatAssemblyType);
     36: extern PetscErrorCode PCApply_MyShell(PC,Vec,Vec);
     37: extern PetscErrorCode SNESComputeJacobian_MyShell(SNES,Vec,Mat,Mat,void*);

     39: int main(int argc,char **argv)
     40: {
     41:   AppCtx         user;                /* user-defined work context */
     42:   PetscInt       mx,my;
     44:   MPI_Comm       comm;
     45:   DM             da;
     46:   Vec            x;
     47:   Mat            J = NULL,Jmf = NULL;
     48:   MatShellCtx    matshellctx;
     49:   PetscInt       mlocal,nlocal;
     50:   PC             pc;
     51:   KSP            ksp;
     52:   PetscBool      errorinmatmult = PETSC_FALSE,errorinpcapply = PETSC_FALSE,errorinpcsetup = PETSC_FALSE;

     54:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
     55:   PetscOptionsGetBool(NULL,NULL,"-error_in_matmult",&errorinmatmult,NULL);
     56:   PetscOptionsGetBool(NULL,NULL,"-error_in_pcapply",&errorinpcapply,NULL);
     57:   PetscOptionsGetBool(NULL,NULL,"-error_in_pcsetup",&errorinpcsetup,NULL);
     58:   user.errorindomain = PETSC_FALSE;
     59:   PetscOptionsGetBool(NULL,NULL,"-error_in_domain",&user.errorindomain,NULL);
     60:   user.errorindomainmf = PETSC_FALSE;
     61:   PetscOptionsGetBool(NULL,NULL,"-error_in_domainmf",&user.errorindomainmf,NULL);

     63:   comm = PETSC_COMM_WORLD;
     64:   SNESCreate(comm,&user.snes);

     66:   /*
     67:       Create distributed array object to manage parallel grid and vectors
     68:       for principal unknowns (x) and governing residuals (f)
     69:   */
     70:   DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,4,4,PETSC_DECIDE,PETSC_DECIDE,4,1,0,0,&da);
     71:   DMSetFromOptions(da);
     72:   DMSetUp(da);
     73:   SNESSetDM(user.snes,da);

     75:   DMDAGetInfo(da,0,&mx,&my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
     76:                      PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
     77:   /*
     78:      Problem parameters (velocity of lid, prandtl, and grashof numbers)
     79:   */
     80:   user.lidvelocity = 1.0/(mx*my);
     81:   user.prandtl     = 1.0;
     82:   user.grashof     = 1.0;

     84:   PetscOptionsGetReal(NULL,NULL,"-lidvelocity",&user.lidvelocity,NULL);
     85:   PetscOptionsGetReal(NULL,NULL,"-prandtl",&user.prandtl,NULL);
     86:   PetscOptionsGetReal(NULL,NULL,"-grashof",&user.grashof,NULL);
     87:   PetscOptionsHasName(NULL,NULL,"-contours",&user.draw_contours);

     89:   DMDASetFieldName(da,0,"x_velocity");
     90:   DMDASetFieldName(da,1,"y_velocity");
     91:   DMDASetFieldName(da,2,"Omega");
     92:   DMDASetFieldName(da,3,"temperature");

     94:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     95:      Create user context, set problem data, create vector data structures.
     96:      Also, compute the initial guess.
     97:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     99:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    100:      Create nonlinear solver context

    102:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    103:   DMSetApplicationContext(da,&user);
    104:   DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,&user);

    106:   if (errorinmatmult) {
    107:     MatCreateSNESMF(user.snes,&Jmf);
    108:     MatSetFromOptions(Jmf);
    109:     MatGetLocalSize(Jmf,&mlocal,&nlocal);
    110:     matshellctx.Jmf = Jmf;
    111:     MatCreateShell(PetscObjectComm((PetscObject)Jmf),mlocal,nlocal,PETSC_DECIDE,PETSC_DECIDE,&matshellctx,&J);
    112:     MatShellSetOperation(J,MATOP_MULT,(void (*)(void))MatMult_MyShell);
    113:     MatShellSetOperation(J,MATOP_ASSEMBLY_END,(void (*)(void))MatAssemblyEnd_MyShell);
    114:     SNESSetJacobian(user.snes,J,J,MatMFFDComputeJacobian,NULL);
    115:   }

    117:   SNESSetFromOptions(user.snes);
    118:   PetscPrintf(comm,"lid velocity = %g, prandtl # = %g, grashof # = %g\n",(double)user.lidvelocity,(double)user.prandtl,(double)user.grashof);

    120:   if (errorinpcapply) {
    121:     SNESGetKSP(user.snes,&ksp);
    122:     KSPGetPC(ksp,&pc);
    123:     PCSetType(pc,PCSHELL);
    124:     PCShellSetApply(pc,PCApply_MyShell);
    125:   }

    127:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    128:      Solve the nonlinear system
    129:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    130:   DMCreateGlobalVector(da,&x);
    131:   FormInitialGuess(&user,da,x);

    133:   if (errorinpcsetup) {
    134:     SNESSetUp(user.snes);
    135:     SNESSetJacobian(user.snes,NULL,NULL,SNESComputeJacobian_MyShell,NULL);
    136:   }
    137:   SNESSolve(user.snes,NULL,x);

    140:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    141:      Free work space.  All PETSc objects should be destroyed when they
    142:      are no longer needed.
    143:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    144:   MatDestroy(&J);
    145:   MatDestroy(&Jmf);
    146:   VecDestroy(&x);
    147:   DMDestroy(&da);
    148:   SNESDestroy(&user.snes);
    149:   PetscFinalize();
    150:   return ierr;
    151: }

    153: /* ------------------------------------------------------------------- */

    156: /*
    157:    FormInitialGuess - Forms initial approximation.

    159:    Input Parameters:
    160:    user - user-defined application context
    161:    X - vector

    163:    Output Parameter:
    164:    X - vector
    165: */
    166: PetscErrorCode FormInitialGuess(AppCtx *user,DM da,Vec X)
    167: {
    168:   PetscInt       i,j,mx,xs,ys,xm,ym;
    170:   PetscReal      grashof,dx;
    171:   Field          **x;

    174:   grashof = user->grashof;

    176:   DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
    177:   dx   = 1.0/(mx-1);

    179:   /*
    180:      Get local grid boundaries (for 2-dimensional DMDA):
    181:        xs, ys   - starting grid indices (no ghost points)
    182:        xm, ym   - widths of local grid (no ghost points)
    183:   */
    184:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

    186:   /*
    187:      Get a pointer to vector data.
    188:        - For default PETSc vectors, VecGetArray() returns a pointer to
    189:          the data array.  Otherwise, the routine is implementation dependent.
    190:        - You MUST call VecRestoreArray() when you no longer need access to
    191:          the array.
    192:   */
    193:   DMDAVecGetArray(da,X,&x);

    195:   /*
    196:      Compute initial guess over the locally owned part of the grid
    197:      Initial condition is motionless fluid and equilibrium temperature
    198:   */
    199:   for (j=ys; j<ys+ym; j++) {
    200:     for (i=xs; i<xs+xm; i++) {
    201:       x[j][i].u     = 0.0;
    202:       x[j][i].v     = 0.0;
    203:       x[j][i].omega = 0.0;
    204:       x[j][i].temp  = (grashof>0)*i*dx;
    205:     }
    206:   }

    208:   /*
    209:      Restore vector
    210:   */
    211:   DMDAVecRestoreArray(da,X,&x);
    212:   return(0);
    213: }

    215: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,Field **x,Field **f,void *ptr)
    216: {
    217:   AppCtx          *user = (AppCtx*)ptr;
    218:   PetscErrorCode  ierr;
    219:   PetscInt        xints,xinte,yints,yinte,i,j;
    220:   PetscReal       hx,hy,dhx,dhy,hxdhy,hydhx;
    221:   PetscReal       grashof,prandtl,lid;
    222:   PetscScalar     u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;
    223:   static PetscInt fail = 0;

    226:   if ((fail++ > 7 && user->errorindomainmf) || (fail++ > 36 && user->errorindomain)){
    227:     PetscMPIInt rank;
    228:     MPI_Comm_rank(PetscObjectComm((PetscObject)user->snes),&rank);
    229:     if (!rank) {
    230:       SNESSetFunctionDomainError(user->snes);
    231:     }
    232:   }
    233:   grashof = user->grashof;
    234:   prandtl = user->prandtl;
    235:   lid     = user->lidvelocity;

    237:   /*
    238:      Define mesh intervals ratios for uniform grid.

    240:      Note: FD formulae below are normalized by multiplying through by
    241:      local volume element (i.e. hx*hy) to obtain coefficients O(1) in two dimensions.

    244:   */
    245:   dhx   = (PetscReal)(info->mx-1);  dhy = (PetscReal)(info->my-1);
    246:   hx    = 1.0/dhx;                   hy = 1.0/dhy;
    247:   hxdhy = hx*dhy;                 hydhx = hy*dhx;

    249:   xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym;

    251:   /* Test whether we are on the bottom edge of the global array */
    252:   if (yints == 0) {
    253:     j     = 0;
    254:     yints = yints + 1;
    255:     /* bottom edge */
    256:     for (i=info->xs; i<info->xs+info->xm; i++) {
    257:       f[j][i].u     = x[j][i].u;
    258:       f[j][i].v     = x[j][i].v;
    259:       f[j][i].omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
    260:       f[j][i].temp  = x[j][i].temp-x[j+1][i].temp;
    261:     }
    262:   }

    264:   /* Test whether we are on the top edge of the global array */
    265:   if (yinte == info->my) {
    266:     j     = info->my - 1;
    267:     yinte = yinte - 1;
    268:     /* top edge */
    269:     for (i=info->xs; i<info->xs+info->xm; i++) {
    270:       f[j][i].u     = x[j][i].u - lid;
    271:       f[j][i].v     = x[j][i].v;
    272:       f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
    273:       f[j][i].temp  = x[j][i].temp-x[j-1][i].temp;
    274:     }
    275:   }

    277:   /* Test whether we are on the left edge of the global array */
    278:   if (xints == 0) {
    279:     i     = 0;
    280:     xints = xints + 1;
    281:     /* left edge */
    282:     for (j=info->ys; j<info->ys+info->ym; j++) {
    283:       f[j][i].u     = x[j][i].u;
    284:       f[j][i].v     = x[j][i].v;
    285:       f[j][i].omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
    286:       f[j][i].temp  = x[j][i].temp;
    287:     }
    288:   }

    290:   /* Test whether we are on the right edge of the global array */
    291:   if (xinte == info->mx) {
    292:     i     = info->mx - 1;
    293:     xinte = xinte - 1;
    294:     /* right edge */
    295:     for (j=info->ys; j<info->ys+info->ym; j++) {
    296:       f[j][i].u     = x[j][i].u;
    297:       f[j][i].v     = x[j][i].v;
    298:       f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
    299:       f[j][i].temp  = x[j][i].temp - (PetscReal)(grashof>0);
    300:     }
    301:   }

    303:   /* Compute over the interior points */
    304:   for (j=yints; j<yinte; j++) {
    305:     for (i=xints; i<xinte; i++) {

    307:       /*
    308:        convective coefficients for upwinding
    309:       */
    310:       vx  = x[j][i].u; avx = PetscAbsScalar(vx);
    311:       vxp = .5*(vx+avx); vxm = .5*(vx-avx);
    312:       vy  = x[j][i].v; avy = PetscAbsScalar(vy);
    313:       vyp = .5*(vy+avy); vym = .5*(vy-avy);

    315:       /* U velocity */
    316:       u         = x[j][i].u;
    317:       uxx       = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
    318:       uyy       = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
    319:       f[j][i].u = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;

    321:       /* V velocity */
    322:       u         = x[j][i].v;
    323:       uxx       = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
    324:       uyy       = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
    325:       f[j][i].v = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;

    327:       /* Omega */
    328:       u             = x[j][i].omega;
    329:       uxx           = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
    330:       uyy           = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
    331:       f[j][i].omega = uxx + uyy + (vxp*(u - x[j][i-1].omega) + vxm*(x[j][i+1].omega - u))*hy +
    332:                       (vyp*(u - x[j-1][i].omega) + vym*(x[j+1][i].omega - u))*hx -
    333:                       .5*grashof*(x[j][i+1].temp - x[j][i-1].temp)*hy;

    335:       /* Temperature */
    336:       u            = x[j][i].temp;
    337:       uxx          = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
    338:       uyy          = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
    339:       f[j][i].temp =  uxx + uyy  + prandtl*((vxp*(u - x[j][i-1].temp) + vxm*(x[j][i+1].temp - u))*hy +
    340:                                             (vyp*(u - x[j-1][i].temp) + vym*(x[j+1][i].temp - u))*hx);
    341:     }
    342:   }

    344:   /*
    345:      Flop count (multiply-adds are counted as 2 operations)
    346:   */
    347:   PetscLogFlops(84.0*info->ym*info->xm);
    348:   return(0);
    349: }

    351: PetscErrorCode MatMult_MyShell(Mat A,Vec x,Vec y)
    352: {
    353:   PetscErrorCode  ierr;
    354:   MatShellCtx     *matshellctx;
    355:   static PetscInt fail = 0;

    358:   MatShellGetContext(A,&matshellctx);
    359:   MatMult(matshellctx->Jmf,x,y);
    360:   if (fail++ > 5) {
    361:     PetscMPIInt rank;
    362:     MPI_Comm_rank(PetscObjectComm((PetscObject)A),&rank);
    363:     if (!rank) {VecSetInf(y);}
    364:   }
    365:   return(0);
    366: }

    368: PetscErrorCode MatAssemblyEnd_MyShell(Mat A,MatAssemblyType tp)
    369: {
    371:   MatShellCtx    *matshellctx;

    374:   MatShellGetContext(A,&matshellctx);
    375:   MatAssemblyEnd(matshellctx->Jmf,tp);
    376:   return(0);
    377: }

    379: PetscErrorCode PCApply_MyShell(PC pc,Vec x,Vec y)
    380: {
    382:   static PetscInt fail = 0;

    385:   VecCopy(x,y);
    386:   if (fail++ > 3) {
    387:     PetscMPIInt rank;
    388:     MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
    389:     if (!rank) {VecSetInf(y);}
    390:   }
    391:   return(0);
    392: }

    394: PETSC_EXTERN PetscErrorCode SNESComputeJacobian_DMDA(SNES,Vec,Mat,Mat,void*);

    396: PetscErrorCode SNESComputeJacobian_MyShell(SNES snes,Vec X,Mat A,Mat B,void *ctx)
    397: {
    398:   static PetscInt fail = 0;
    399:   PetscErrorCode  ierr;

    402:   SNESComputeJacobian_DMDA(snes,X,A,B,ctx);
    403:   if (fail++ > 0) {
    404:     MatZeroEntries(A);
    405:   }
    406:   return(0);
    407: }

    410: /*TEST

    412:    test:
    413:       args: -snes_converged_reason -ksp_converged_reason

    415:    test:
    416:       suffix: 2
    417:       args: -snes_converged_reason -ksp_converged_reason -error_in_matmult

    419:    test:
    420:       suffix: 3
    421:       args: -snes_converged_reason -ksp_converged_reason -error_in_pcapply

    423:    test:
    424:       suffix: 4
    425:       args: -snes_converged_reason -ksp_converged_reason -error_in_pcsetup

    427:    test:
    428:       suffix: 5
    429:       args: -snes_converged_reason -ksp_converged_reason -error_in_pcsetup -pc_type bjacobi

    431:    test:
    432:       suffix: 5_fieldsplit
    433:       args: -snes_converged_reason -ksp_converged_reason -error_in_pcsetup -pc_type fieldsplit
    434:       output_file: output/ex69_5.out

    436:    test:
    437:       suffix: 6
    438:       args: -snes_converged_reason -ksp_converged_reason -error_in_domainmf -snes_mf -pc_type none

    440:    test:
    441:       suffix: 7
    442:       args: -snes_converged_reason -ksp_converged_reason -error_in_domain

    444:    test:
    445:       suffix: 8
    446:       args: -snes_converged_reason -ksp_converged_reason -error_in_domain -snes_mf -pc_type none

    448: TEST*/