Actual source code: multiblock.c

      1: #include <petsc/private/snesimpl.h>
      2: #include <petscdmcomposite.h>

      4: typedef struct _BlockDesc *BlockDesc;
      5: struct _BlockDesc {
      6:   char       *name;     /* Block name */
      7:   PetscInt   nfields;   /* If block is defined on a DA, the number of DA fields */
      8:   PetscInt   *fields;   /* If block is defined on a DA, the list of DA fields */
      9:   IS         is;        /* Index sets defining the block */
     10:   VecScatter sctx;      /* Scatter mapping global Vec to blockVec */
     11:   SNES       snes;      /* Solver for this block */
     12:   Vec        x;
     13:   BlockDesc  next, previous;
     14: };

     16: typedef struct {
     17:   PetscBool       issetup;       /* Flag is true after the all ISs and operators have been defined */
     18:   PetscBool       defined;       /* Flag is true after the blocks have been defined, to prevent more blocks from being added */
     19:   PetscBool       defaultblocks; /* Flag is true for a system with a set of 'k' scalar fields with the same layout (and bs = k) */
     20:   PetscInt        numBlocks;     /* Number of blocks (can be fields, domains, etc.) */
     21:   PetscInt        bs;            /* Block size for IS, Vec and Mat structures */
     22:   PCCompositeType type;          /* Solver combination method (additive, multiplicative, etc.) */
     23:   BlockDesc       blocks;        /* Linked list of block descriptors */
     24: } SNES_Multiblock;

     26: PetscErrorCode SNESReset_Multiblock(SNES snes)
     27: {
     28:   SNES_Multiblock *mb    = (SNES_Multiblock*) snes->data;
     29:   BlockDesc       blocks = mb->blocks, next;
     30:   PetscErrorCode  ierr;

     33:   while (blocks) {
     34:     SNESReset(blocks->snes);
     35: #if 0
     36:     VecDestroy(&blocks->x);
     37: #endif
     38:     VecScatterDestroy(&blocks->sctx);
     39:     ISDestroy(&blocks->is);
     40:     next   = blocks->next;
     41:     blocks = next;
     42:   }
     43:   return(0);
     44: }

     46: /*
     47:   SNESDestroy_Multiblock - Destroys the private SNES_Multiblock context that was created with SNESCreate_Multiblock().

     49:   Input Parameter:
     50: . snes - the SNES context

     52:   Application Interface Routine: SNESDestroy()
     53: */
     54: PetscErrorCode SNESDestroy_Multiblock(SNES snes)
     55: {
     56:   SNES_Multiblock *mb    = (SNES_Multiblock*) snes->data;
     57:   BlockDesc       blocks = mb->blocks, next;
     58:   PetscErrorCode  ierr;

     61:   SNESReset_Multiblock(snes);
     62:   while (blocks) {
     63:     next   = blocks->next;
     64:     SNESDestroy(&blocks->snes);
     65:     PetscFree(blocks->name);
     66:     PetscFree(blocks->fields);
     67:     PetscFree(blocks);
     68:     blocks = next;
     69:   }
     70:   PetscFree(snes->data);
     71:   return(0);
     72: }

     74: /* Precondition: blocksize is set to a meaningful value */
     75: static PetscErrorCode SNESMultiblockSetFieldsRuntime_Private(SNES snes)
     76: {
     77:   SNES_Multiblock *mb = (SNES_Multiblock*) snes->data;
     78:   PetscInt        *ifields;
     79:   PetscInt        i, nfields;
     80:   PetscBool       flg = PETSC_TRUE;
     81:   char            optionname[128], name[8];
     82:   PetscErrorCode  ierr;

     85:   PetscMalloc1(mb->bs, &ifields);
     86:   for (i = 0;; ++i) {
     87:     PetscSNPrintf(name, sizeof(name), "%D", i);
     88:     PetscSNPrintf(optionname, sizeof(optionname), "-snes_multiblock_%D_fields", i);
     89:     nfields = mb->bs;
     90:     PetscOptionsGetIntArray(NULL,((PetscObject) snes)->prefix, optionname, ifields, &nfields, &flg);
     91:     if (!flg) break;
     92:     if (!nfields) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot list zero fields");
     93:     SNESMultiblockSetFields(snes, name, nfields, ifields);
     94:   }
     95:   if (i > 0) {
     96:     /* Makes command-line setting of blocks take precedence over setting them in code.
     97:        Otherwise subsequent calls to SNESMultiblockSetIS() or SNESMultiblockSetFields() would
     98:        create new blocks, which would probably not be what the user wanted. */
     99:     mb->defined = PETSC_TRUE;
    100:   }
    101:   PetscFree(ifields);
    102:   return(0);
    103: }

    105: static PetscErrorCode SNESMultiblockSetDefaults(SNES snes)
    106: {
    107:   SNES_Multiblock *mb    = (SNES_Multiblock*) snes->data;
    108:   BlockDesc       blocks = mb->blocks;
    109:   PetscInt        i;
    110:   PetscErrorCode  ierr;

    113:   if (!blocks) {
    114:     if (snes->dm) {
    115:       PetscBool dmcomposite;

    117:       PetscObjectTypeCompare((PetscObject) snes->dm, DMCOMPOSITE, &dmcomposite);
    118:       if (dmcomposite) {
    119:         PetscInt nDM;
    120:         IS       *fields;

    122:         PetscInfo(snes,"Setting up physics based multiblock solver using the embedded DM\n");
    123:         DMCompositeGetNumberDM(snes->dm, &nDM);
    124:         DMCompositeGetGlobalISs(snes->dm, &fields);
    125:         for (i = 0; i < nDM; ++i) {
    126:           char name[8];

    128:           PetscSNPrintf(name, sizeof(name), "%D", i);
    129:           SNESMultiblockSetIS(snes, name, fields[i]);
    130:           ISDestroy(&fields[i]);
    131:         }
    132:         PetscFree(fields);
    133:       }
    134:     } else {
    135:       PetscBool flg    = PETSC_FALSE;
    136:       PetscBool stokes = PETSC_FALSE;

    138:       if (mb->bs <= 0) {
    139:         if (snes->jacobian_pre) {
    140:           MatGetBlockSize(snes->jacobian_pre, &mb->bs);
    141:         } else mb->bs = 1;
    142:       }

    144:       PetscOptionsGetBool(NULL,((PetscObject) snes)->prefix, "-snes_multiblock_default", &flg, NULL);
    145:       PetscOptionsGetBool(NULL,((PetscObject) snes)->prefix, "-snes_multiblock_detect_saddle_point", &stokes, NULL);
    146:       if (stokes) {
    147:         IS       zerodiags, rest;
    148:         PetscInt nmin, nmax;

    150:         MatGetOwnershipRange(snes->jacobian_pre, &nmin, &nmax);
    151:         MatFindZeroDiagonals(snes->jacobian_pre, &zerodiags);
    152:         ISComplement(zerodiags, nmin, nmax, &rest);
    153:         SNESMultiblockSetIS(snes, "0", rest);
    154:         SNESMultiblockSetIS(snes, "1", zerodiags);
    155:         ISDestroy(&zerodiags);
    156:         ISDestroy(&rest);
    157:       } else {
    158:         if (!flg) {
    159:           /* Allow user to set fields from command line, if bs was known at the time of SNESSetFromOptions_Multiblock()
    160:            then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */
    161:           SNESMultiblockSetFieldsRuntime_Private(snes);
    162:           if (mb->defined) {PetscInfo(snes, "Blocks defined using the options database\n");}
    163:         }
    164:         if (flg || !mb->defined) {
    165:           PetscInfo(snes, "Using default splitting of fields\n");
    166:           for (i = 0; i < mb->bs; ++i) {
    167:             char name[8];

    169:             PetscSNPrintf(name, sizeof(name), "%D", i);
    170:             SNESMultiblockSetFields(snes, name, 1, &i);
    171:           }
    172:           mb->defaultblocks = PETSC_TRUE;
    173:         }
    174:       }
    175:     }
    176:   } else if (mb->numBlocks == 1) {
    177:     if (blocks->is) {
    178:       IS       is2;
    179:       PetscInt nmin, nmax;

    181:       MatGetOwnershipRange(snes->jacobian_pre, &nmin, &nmax);
    182:       ISComplement(blocks->is, nmin, nmax, &is2);
    183:       SNESMultiblockSetIS(snes, "1", is2);
    184:       ISDestroy(&is2);
    185:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least two sets of fields to SNES multiblock");
    186:   }
    187:   if (mb->numBlocks < 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unhandled case, must have at least two blocks");
    188:   return(0);
    189: }

    191: /*
    192:    SNESSetUp_Multiblock - Sets up the internal data structures for the later use
    193:    of the SNESMULTIBLOCK nonlinear solver.

    195:    Input Parameters:
    196: +  snes - the SNES context
    197: -  x - the solution vector

    199:    Application Interface Routine: SNESSetUp()
    200: */
    201: PetscErrorCode SNESSetUp_Multiblock(SNES snes)
    202: {
    203:   SNES_Multiblock *mb = (SNES_Multiblock*) snes->data;
    204:   BlockDesc       blocks;
    205:   PetscInt        i, numBlocks;
    206:   PetscErrorCode  ierr;

    209:   SNESMultiblockSetDefaults(snes);
    210:   numBlocks = mb->numBlocks;
    211:   blocks    = mb->blocks;

    213:   /* Create ISs */
    214:   if (!mb->issetup) {
    215:     PetscInt  ccsize, rstart, rend, nslots, bs;
    216:     PetscBool sorted;

    218:     mb->issetup = PETSC_TRUE;
    219:     bs          = mb->bs;
    220:     MatGetOwnershipRange(snes->jacobian_pre, &rstart, &rend);
    221:     MatGetLocalSize(snes->jacobian_pre, NULL, &ccsize);
    222:     nslots      = (rend - rstart)/bs;
    223:     for (i = 0; i < numBlocks; ++i) {
    224:       if (mb->defaultblocks) {
    225:         ISCreateStride(PetscObjectComm((PetscObject)snes), nslots, rstart+i, numBlocks, &blocks->is);
    226:       } else if (!blocks->is) {
    227:         if (blocks->nfields > 1) {
    228:           PetscInt *ii, j, k, nfields = blocks->nfields, *fields = blocks->fields;

    230:           PetscMalloc1(nfields*nslots, &ii);
    231:           for (j = 0; j < nslots; ++j) {
    232:             for (k = 0; k < nfields; ++k) {
    233:               ii[nfields*j + k] = rstart + bs*j + fields[k];
    234:             }
    235:           }
    236:           ISCreateGeneral(PetscObjectComm((PetscObject)snes), nslots*nfields, ii, PETSC_OWN_POINTER, &blocks->is);
    237:         } else {
    238:           ISCreateStride(PetscObjectComm((PetscObject)snes), nslots, rstart+blocks->fields[0], bs, &blocks->is);
    239:         }
    240:       }
    241:       ISSorted(blocks->is, &sorted);
    242:       if (!sorted) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Fields must be sorted when creating split");
    243:       blocks = blocks->next;
    244:     }
    245:   }

    247: #if 0
    248:   /* Create matrices */
    249:   ilink = jac->head;
    250:   if (!jac->pmat) {
    251:     PetscMalloc1(nsplit,&jac->pmat);
    252:     for (i=0; i<nsplit; i++) {
    253:       MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is,MAT_INITIAL_MATRIX,&jac->pmat[i]);
    254:       ilink = ilink->next;
    255:     }
    256:   } else {
    257:     for (i=0; i<nsplit; i++) {
    258:       MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is,MAT_REUSE_MATRIX,&jac->pmat[i]);
    259:       ilink = ilink->next;
    260:     }
    261:   }
    262:   if (jac->realdiagonal) {
    263:     ilink = jac->head;
    264:     if (!jac->mat) {
    265:       PetscMalloc1(nsplit,&jac->mat);
    266:       for (i=0; i<nsplit; i++) {
    267:         MatCreateSubMatrix(pc->mat,ilink->is,ilink->is,MAT_INITIAL_MATRIX,&jac->mat[i]);
    268:         ilink = ilink->next;
    269:       }
    270:     } else {
    271:       for (i=0; i<nsplit; i++) {
    272:         if (jac->mat[i]) {MatCreateSubMatrix(pc->mat,ilink->is,ilink->is,MAT_REUSE_MATRIX,&jac->mat[i]);}
    273:         ilink = ilink->next;
    274:       }
    275:     }
    276:   } else jac->mat = jac->pmat;
    277: #endif

    279: #if 0
    280:   if (jac->type != PC_COMPOSITE_ADDITIVE  && jac->type != PC_COMPOSITE_SCHUR) {
    281:     /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */
    282:     ilink = jac->head;
    283:     if (!jac->Afield) {
    284:       PetscMalloc1(nsplit,&jac->Afield);
    285:       for (i=0; i<nsplit; i++) {
    286:         MatCreateSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);
    287:         ilink = ilink->next;
    288:       }
    289:     } else {
    290:       for (i=0; i<nsplit; i++) {
    291:         MatCreateSubMatrix(pc->mat,ilink->is,NULL,MAT_REUSE_MATRIX,&jac->Afield[i]);
    292:         ilink = ilink->next;
    293:       }
    294:     }
    295:   }
    296: #endif

    298:   if (mb->type == PC_COMPOSITE_SCHUR) {
    299: #if 0
    300:     IS       ccis;
    301:     PetscInt rstart,rend;
    302:     if (nsplit != 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields");

    304:     /* When extracting off-diagonal submatrices, we take complements from this range */
    305:     MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);

    307:     /* need to handle case when one is resetting up the preconditioner */
    308:     if (jac->schur) {
    309:       ilink = jac->head;
    310:       ISComplement(ilink->is,rstart,rend,&ccis);
    311:       MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);
    312:       ISDestroy(&ccis);
    313:       ilink = ilink->next;
    314:       ISComplement(ilink->is,rstart,rend,&ccis);
    315:       MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);
    316:       ISDestroy(&ccis);
    317:       MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->pmat[1]);
    318:       KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),pc->flag);

    320:     } else {
    321:       KSP  ksp;
    322:       char schurprefix[256];

    324:       /* extract the A01 and A10 matrices */
    325:       ilink = jac->head;
    326:       ISComplement(ilink->is,rstart,rend,&ccis);
    327:       MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
    328:       ISDestroy(&ccis);
    329:       ilink = ilink->next;
    330:       ISComplement(ilink->is,rstart,rend,&ccis);
    331:       MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
    332:       ISDestroy(&ccis);
    333:       /* Use mat[0] (diagonal block of the real matrix) preconditioned by pmat[0] */
    334:       MatCreateSchurComplement(jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1],&jac->schur);
    335:       /* set tabbing and options prefix of KSP inside the MatSchur */
    336:       MatSchurComplementGetKSP(jac->schur,&ksp);
    337:       PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,2);
    338:       PetscSNPrintf(schurprefix,sizeof(schurprefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",jac->head->splitname);
    339:       KSPSetOptionsPrefix(ksp,schurprefix);
    340:       MatSetFromOptions(jac->schur);

    342:       KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);
    343:       PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);
    344:       PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);
    345:       KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));
    346:       if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) {
    347:         PC pc;
    348:         KSPGetPC(jac->kspschur,&pc);
    349:         PCSetType(pc,PCNONE);
    350:         /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */
    351:       }
    352:       PetscSNPrintf(schurprefix,sizeof(schurprefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);
    353:       KSPSetOptionsPrefix(jac->kspschur,schurprefix);
    354:       /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
    355:       KSPSetFromOptions(jac->kspschur);

    357:       PetscMalloc2(2,&jac->x,2,&jac->y);
    358:       MatCreateVecs(jac->pmat[0],&jac->x[0],&jac->y[0]);
    359:       MatCreateVecs(jac->pmat[1],&jac->x[1],&jac->y[1]);
    360:       ilink    = jac->head;
    361:       ilink->x = jac->x[0]; ilink->y = jac->y[0];
    362:       ilink    = ilink->next;
    363:       ilink->x = jac->x[1]; ilink->y = jac->y[1];
    364:     }
    365: #endif
    366:   } else {
    367:     /* Set up the individual SNESs */
    368:     blocks = mb->blocks;
    369:     i      = 0;
    370:     while (blocks) {
    371:       /*TODO: Set these correctly */
    372:       /*SNESSetFunction(blocks->snes, blocks->x, func);*/
    373:       /*SNESSetJacobian(blocks->snes, blocks->x, jac);*/
    374:       VecDuplicate(blocks->snes->vec_sol, &blocks->x);
    375:       /* really want setfromoptions called in SNESSetFromOptions_Multiblock(), but it is not ready yet */
    376:       SNESSetFromOptions(blocks->snes);
    377:       SNESSetUp(blocks->snes);
    378:       blocks = blocks->next;
    379:       i++;
    380:     }
    381:   }

    383:   /* Compute scatter contexts needed by multiplicative versions and non-default splits */
    384:   if (!mb->blocks->sctx) {
    385:     Vec xtmp;

    387:     blocks = mb->blocks;
    388:     MatCreateVecs(snes->jacobian_pre, &xtmp, NULL);
    389:     while (blocks) {
    390:       VecScatterCreate(xtmp, blocks->is, blocks->x, NULL, &blocks->sctx);
    391:       blocks = blocks->next;
    392:     }
    393:     VecDestroy(&xtmp);
    394:   }
    395:   return(0);
    396: }

    398: /*
    399:   SNESSetFromOptions_Multiblock - Sets various parameters for the SNESMULTIBLOCK method.

    401:   Input Parameter:
    402: . snes - the SNES context

    404:   Application Interface Routine: SNESSetFromOptions()
    405: */
    406: static PetscErrorCode SNESSetFromOptions_Multiblock(PetscOptionItems *PetscOptionsObject,SNES snes)
    407: {
    408:   SNES_Multiblock *mb = (SNES_Multiblock*) snes->data;
    409:   PCCompositeType ctype;
    410:   PetscInt        bs;
    411:   PetscBool       flg;
    412:   PetscErrorCode  ierr;

    415:   PetscOptionsHead(PetscOptionsObject,"SNES Multiblock options");
    416:   PetscOptionsInt("-snes_multiblock_block_size", "Blocksize that defines number of fields", "PCFieldSplitSetBlockSize", mb->bs, &bs, &flg);
    417:   if (flg) {SNESMultiblockSetBlockSize(snes, bs);}
    418:   PetscOptionsEnum("-snes_multiblock_type", "Type of composition", "PCFieldSplitSetType", PCCompositeTypes, (PetscEnum) mb->type, (PetscEnum*) &ctype, &flg);
    419:   if (flg) {
    420:     SNESMultiblockSetType(snes,ctype);
    421:   }
    422:   /* Only setup fields once */
    423:   if ((mb->bs > 0) && (mb->numBlocks == 0)) {
    424:     /* only allow user to set fields from command line if bs is already known, otherwise user can set them in SNESMultiblockSetDefaults() */
    425:     SNESMultiblockSetFieldsRuntime_Private(snes);
    426:     if (mb->defined) {PetscInfo(snes, "Blocks defined using the options database\n");}
    427:   }
    428:   PetscOptionsTail();
    429:   return(0);
    430: }

    432: /*
    433:   SNESView_Multiblock - Prints info from the SNESMULTIBLOCK data structure.

    435:   Input Parameters:
    436: + SNES - the SNES context
    437: - viewer - visualization context

    439:   Application Interface Routine: SNESView()
    440: */
    441: static PetscErrorCode SNESView_Multiblock(SNES snes, PetscViewer viewer)
    442: {
    443:   SNES_Multiblock *mb    = (SNES_Multiblock*) snes->data;
    444:   BlockDesc       blocks = mb->blocks;
    445:   PetscBool       iascii;
    446:   PetscErrorCode  ierr;

    449:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
    450:   if (iascii) {
    451:     PetscViewerASCIIPrintf(viewer,"  Multiblock with %s composition: total blocks = %D, blocksize = %D\n", PCCompositeTypes[mb->type], mb->numBlocks, mb->bs);
    452:     PetscViewerASCIIPrintf(viewer,"  Solver info for each split is in the following SNES objects:\n");
    453:     PetscViewerASCIIPushTab(viewer);
    454:     while (blocks) {
    455:       if (blocks->fields) {
    456:         PetscInt j;

    458:         PetscViewerASCIIPrintf(viewer, "  Block %s Fields ", blocks->name);
    459:         PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
    460:         for (j = 0; j < blocks->nfields; ++j) {
    461:           if (j > 0) {
    462:             PetscViewerASCIIPrintf(viewer, ",");
    463:           }
    464:           PetscViewerASCIIPrintf(viewer, " %D", blocks->fields[j]);
    465:         }
    466:         PetscViewerASCIIPrintf(viewer, "\n");
    467:         PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
    468:       } else {
    469:         PetscViewerASCIIPrintf(viewer, "  Block %s Defined by IS\n", blocks->name);
    470:       }
    471:       SNESView(blocks->snes, viewer);
    472:       blocks = blocks->next;
    473:     }
    474:     PetscViewerASCIIPopTab(viewer);
    475:   }
    476:   return(0);
    477: }

    479: /*
    480:   SNESSolve_Multiblock - Solves a nonlinear system with the Multiblock method.

    482:   Input Parameters:
    483: . snes - the SNES context

    485:   Output Parameter:
    486: . outits - number of iterations until termination

    488:   Application Interface Routine: SNESSolve()
    489: */
    490: PetscErrorCode SNESSolve_Multiblock(SNES snes)
    491: {
    492:   SNES_Multiblock *mb = (SNES_Multiblock*) snes->data;
    493:   Vec             X, Y, F;
    494:   PetscReal       fnorm;
    495:   PetscInt        maxits, i;
    496:   PetscErrorCode  ierr;

    499:   if (snes->xl || snes->xu || snes->ops->computevariablebounds) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);

    501:   snes->reason = SNES_CONVERGED_ITERATING;

    503:   maxits = snes->max_its;        /* maximum number of iterations */
    504:   X      = snes->vec_sol;        /* X^n */
    505:   Y      = snes->vec_sol_update; /* \tilde X */
    506:   F      = snes->vec_func;       /* residual vector */

    508:   VecSetBlockSize(X, mb->bs);
    509:   VecSetBlockSize(Y, mb->bs);
    510:   VecSetBlockSize(F, mb->bs);
    511:   PetscObjectSAWsTakeAccess((PetscObject)snes);
    512:   snes->iter = 0;
    513:   snes->norm = 0.;
    514:   PetscObjectSAWsGrantAccess((PetscObject)snes);

    516:   if (!snes->vec_func_init_set) {
    517:     SNESComputeFunction(snes, X, F);
    518:   } else snes->vec_func_init_set = PETSC_FALSE;

    520:   VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
    521:   SNESCheckFunctionNorm(snes,fnorm);
    522:   PetscObjectSAWsTakeAccess((PetscObject)snes);
    523:   snes->norm = fnorm;
    524:   PetscObjectSAWsGrantAccess((PetscObject)snes);
    525:   SNESLogConvergenceHistory(snes,fnorm,0);
    526:   SNESMonitor(snes,0,fnorm);

    528:   /* test convergence */
    529:   (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
    530:   if (snes->reason) return(0);

    532:   for (i = 0; i < maxits; i++) {
    533:     /* Call general purpose update function */
    534:     if (snes->ops->update) {
    535:       (*snes->ops->update)(snes, snes->iter);
    536:     }
    537:     /* Compute X^{new} from subsolves */
    538:     if (mb->type == PC_COMPOSITE_ADDITIVE) {
    539:       BlockDesc blocks = mb->blocks;

    541:       if (mb->defaultblocks) {
    542:         /*TODO: Make an array of Vecs for this */
    543:         /*VecStrideGatherAll(X, mb->x, INSERT_VALUES);*/
    544:         while (blocks) {
    545:           SNESSolve(blocks->snes, NULL, blocks->x);
    546:           blocks = blocks->next;
    547:         }
    548:         /*VecStrideScatterAll(mb->x, X, INSERT_VALUES);*/
    549:       } else {
    550:         while (blocks) {
    551:           VecScatterBegin(blocks->sctx, X, blocks->x, INSERT_VALUES, SCATTER_FORWARD);
    552:           VecScatterEnd(blocks->sctx, X, blocks->x, INSERT_VALUES, SCATTER_FORWARD);
    553:           SNESSolve(blocks->snes, NULL, blocks->x);
    554:           VecScatterBegin(blocks->sctx, blocks->x, X, INSERT_VALUES, SCATTER_REVERSE);
    555:           VecScatterEnd(blocks->sctx, blocks->x, X, INSERT_VALUES, SCATTER_REVERSE);
    556:           blocks = blocks->next;
    557:         }
    558:       }
    559:     } else SETERRQ1(PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "Unsupported or unknown composition", (int) mb->type);
    560:     /* Compute F(X^{new}) */
    561:     SNESComputeFunction(snes, X, F);
    562:     VecNorm(F, NORM_2, &fnorm);
    563:     SNESCheckFunctionNorm(snes,fnorm);

    565:     if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >=0) {
    566:       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
    567:       break;
    568:     }

    570:     /* Monitor convergence */
    571:     PetscObjectSAWsTakeAccess((PetscObject)snes);
    572:     snes->iter = i+1;
    573:     snes->norm = fnorm;
    574:     PetscObjectSAWsGrantAccess((PetscObject)snes);
    575:     SNESLogConvergenceHistory(snes,snes->norm,0);
    576:     SNESMonitor(snes,snes->iter,snes->norm);
    577:     /* Test for convergence */
    578:     (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
    579:     if (snes->reason) break;
    580:   }
    581:   if (i == maxits) {
    582:     PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", maxits);
    583:     if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
    584:   }
    585:   return(0);
    586: }

    588: PetscErrorCode SNESMultiblockSetFields_Default(SNES snes, const char name[], PetscInt n, const PetscInt fields[])
    589: {
    590:   SNES_Multiblock *mb = (SNES_Multiblock*) snes->data;
    591:   BlockDesc       newblock, next = mb->blocks;
    592:   char            prefix[128];
    593:   PetscInt        i;
    594:   PetscErrorCode  ierr;

    597:   if (mb->defined) {
    598:     PetscInfo1(snes, "Ignoring new block \"%s\" because the blocks have already been defined\n", name);
    599:     return(0);
    600:   }
    601:   for (i = 0; i < n; ++i) {
    602:     if (fields[i] >= mb->bs) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field %D requested but only %D exist", fields[i], mb->bs);
    603:     if (fields[i] < 0)       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative field %D requested", fields[i]);
    604:   }
    605:   PetscNew(&newblock);
    606:   if (name) {
    607:     PetscStrallocpy(name, &newblock->name);
    608:   } else {
    609:     PetscInt len = floor(log10(mb->numBlocks))+1;

    611:     PetscMalloc1(len+1, &newblock->name);
    612:     PetscSNPrintf(newblock->name, len, "%s", mb->numBlocks);
    613:   }
    614:   newblock->nfields = n;

    616:   PetscMalloc1(n, &newblock->fields);
    617:   PetscArraycpy(newblock->fields, fields, n);

    619:   newblock->next = NULL;

    621:   SNESCreate(PetscObjectComm((PetscObject)snes), &newblock->snes);
    622:   PetscObjectIncrementTabLevel((PetscObject) newblock->snes, (PetscObject) snes, 1);
    623:   SNESSetType(newblock->snes, SNESNRICHARDSON);
    624:   PetscLogObjectParent((PetscObject) snes, (PetscObject) newblock->snes);
    625:   PetscSNPrintf(prefix, sizeof(prefix), "%smultiblock_%s_", ((PetscObject) snes)->prefix ? ((PetscObject) snes)->prefix : "", newblock->name);
    626:   SNESSetOptionsPrefix(newblock->snes, prefix);

    628:   if (!next) {
    629:     mb->blocks         = newblock;
    630:     newblock->previous = NULL;
    631:   } else {
    632:     while (next->next) {
    633:       next = next->next;
    634:     }
    635:     next->next         = newblock;
    636:     newblock->previous = next;
    637:   }
    638:   mb->numBlocks++;
    639:   return(0);
    640: }

    642: PetscErrorCode SNESMultiblockSetIS_Default(SNES snes, const char name[], IS is)
    643: {
    644:   SNES_Multiblock *mb = (SNES_Multiblock*) snes->data;
    645:   BlockDesc       newblock, next = mb->blocks;
    646:   char            prefix[128];
    647:   PetscErrorCode  ierr;

    650:   if (mb->defined) {
    651:     PetscInfo1(snes, "Ignoring new block \"%s\" because the blocks have already been defined\n", name);
    652:     return(0);
    653:   }
    654:   PetscNew(&newblock);
    655:   if (name) {
    656:     PetscStrallocpy(name, &newblock->name);
    657:   } else {
    658:     PetscInt len = floor(log10(mb->numBlocks))+1;

    660:     PetscMalloc1(len+1, &newblock->name);
    661:     PetscSNPrintf(newblock->name, len, "%s", mb->numBlocks);
    662:   }
    663:   newblock->is = is;

    665:   PetscObjectReference((PetscObject) is);

    667:   newblock->next = NULL;

    669:   SNESCreate(PetscObjectComm((PetscObject)snes), &newblock->snes);
    670:   PetscObjectIncrementTabLevel((PetscObject) newblock->snes, (PetscObject) snes, 1);
    671:   SNESSetType(newblock->snes, SNESNRICHARDSON);
    672:   PetscLogObjectParent((PetscObject) snes, (PetscObject) newblock->snes);
    673:   PetscSNPrintf(prefix, sizeof(prefix), "%smultiblock_%s_", ((PetscObject) snes)->prefix ? ((PetscObject) snes)->prefix : "", newblock->name);
    674:   SNESSetOptionsPrefix(newblock->snes, prefix);

    676:   if (!next) {
    677:     mb->blocks         = newblock;
    678:     newblock->previous = NULL;
    679:   } else {
    680:     while (next->next) {
    681:       next = next->next;
    682:     }
    683:     next->next         = newblock;
    684:     newblock->previous = next;
    685:   }
    686:   mb->numBlocks++;
    687:   return(0);
    688: }

    690: PetscErrorCode  SNESMultiblockSetBlockSize_Default(SNES snes, PetscInt bs)
    691: {
    692:   SNES_Multiblock *mb = (SNES_Multiblock*) snes->data;

    695:   if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_OUTOFRANGE, "Blocksize must be positive, you gave %D", bs);
    696:   if (mb->bs > 0 && mb->bs != bs) SETERRQ2(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Cannot change blocksize from %D to %D after it has been set", mb->bs, bs);
    697:   mb->bs = bs;
    698:   return(0);
    699: }

    701: PetscErrorCode SNESMultiblockGetSubSNES_Default(SNES snes, PetscInt *n, SNES **subsnes)
    702: {
    703:   SNES_Multiblock *mb    = (SNES_Multiblock*) snes->data;
    704:   BlockDesc       blocks = mb->blocks;
    705:   PetscInt        cnt    = 0;
    706:   PetscErrorCode  ierr;

    709:   PetscMalloc1(mb->numBlocks, subsnes);
    710:   while (blocks) {
    711:     (*subsnes)[cnt++] = blocks->snes;
    712:     blocks            = blocks->next;
    713:   }
    714:   if (cnt != mb->numBlocks) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Corrupt SNESMULTIBLOCK object: number of blocks in linked list %D does not match number in object %D", cnt, mb->numBlocks);

    716:   if (n) *n = mb->numBlocks;
    717:   return(0);
    718: }

    720: PetscErrorCode  SNESMultiblockSetType_Default(SNES snes, PCCompositeType type)
    721: {
    722:   SNES_Multiblock *mb = (SNES_Multiblock*) snes->data;
    723:   PetscErrorCode  ierr;

    726:   mb->type = type;
    727:   if (type == PC_COMPOSITE_SCHUR) {
    728: #if 1
    729:     SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "The Schur composite type is not yet supported");
    730: #else
    731:     snes->ops->solve = SNESSolve_Multiblock_Schur;
    732:     snes->ops->view  = SNESView_Multiblock_Schur;

    734:     PetscObjectComposeFunction((PetscObject) snes, "SNESMultiblockGetSubSNES_C", SNESMultiblockGetSubSNES_Schur);
    735:     PetscObjectComposeFunction((PetscObject) snes, "SNESMultiblockSchurPrecondition_C", SNESMultiblockSchurPrecondition_Default);
    736: #endif
    737:   } else {
    738:     snes->ops->solve = SNESSolve_Multiblock;
    739:     snes->ops->view  = SNESView_Multiblock;

    741:     PetscObjectComposeFunction((PetscObject) snes, "SNESMultiblockGetSubSNES_C", SNESMultiblockGetSubSNES_Default);
    742:     PetscObjectComposeFunction((PetscObject) snes, "SNESMultiblockSchurPrecondition_C", 0);
    743:   }
    744:   return(0);
    745: }

    747: /*@
    748:   SNESMultiblockSetFields - Sets the fields for one particular block in the solver

    750:   Logically Collective on SNES

    752:   Input Parameters:
    753: + snes   - the solver
    754: . name   - name of this block, if NULL the number of the block is used
    755: . n      - the number of fields in this block
    756: - fields - the fields in this block

    758:   Level: intermediate

    760:   Notes:
    761:     Use SNESMultiblockSetIS() to set a completely general set of row indices as a block.

    763:   The SNESMultiblockSetFields() is for defining blocks as a group of strided indices, or fields.
    764:   For example, if the vector block size is three then one can define a block as field 0, or
    765:   1 or 2, or field 0,1 or 0,2 or 1,2 which means
    766:     0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x....
    767:   where the numbered entries indicate what is in the block.

    769:   This function is called once per block (it creates a new block each time). Solve options
    770:   for this block will be available under the prefix -multiblock_BLOCKNAME_.

    772: .seealso: SNESMultiblockGetSubSNES(), SNESMULTIBLOCK, SNESMultiblockSetBlockSize(), SNESMultiblockSetIS()
    773: @*/
    774: PetscErrorCode SNESMultiblockSetFields(SNES snes, const char name[], PetscInt n, const PetscInt *fields)
    775: {

    781:   if (n < 1) SETERRQ2(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_OUTOFRANGE, "Provided number of fields %D in split \"%s\" not positive", n, name);
    783:   PetscTryMethod(snes, "SNESMultiblockSetFields_C", (SNES, const char[], PetscInt, const PetscInt*), (snes, name, n, fields));
    784:   return(0);
    785: }

    787: /*@
    788:   SNESMultiblockSetIS - Sets the global row indices for the block

    790:   Logically Collective on SNES

    792:   Input Parameters:
    793: + snes - the solver context
    794: . name - name of this block, if NULL the number of the block is used
    795: - is   - the index set that defines the global row indices in this block

    797:   Notes:
    798:   Use SNESMultiblockSetFields(), for blocks defined by strides.

    800:   This function is called once per block (it creates a new block each time). Solve options
    801:   for this block will be available under the prefix -multiblock_BLOCKNAME_.

    803:   Level: intermediate

    805: .seealso: SNESMultiblockGetSubSNES(), SNESMULTIBLOCK, SNESMultiblockSetBlockSize()
    806: @*/
    807: PetscErrorCode SNESMultiblockSetIS(SNES snes, const char name[], IS is)
    808: {

    815:   PetscTryMethod(snes, "SNESMultiblockSetIS_C", (SNES, const char[], IS), (snes, name, is));
    816:   return(0);
    817: }

    819: /*@
    820:   SNESMultiblockSetType - Sets the type of block combination.

    822:   Collective on SNES

    824:   Input Parameters:
    825: + snes - the solver context
    826: - type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE

    828:   Options Database Key:
    829: . -snes_multiblock_type <type: one of multiplicative, additive, symmetric_multiplicative> - Sets block combination type

    831:   Level: Developer

    833: .seealso: PCCompositeSetType()
    834: @*/
    835: PetscErrorCode SNESMultiblockSetType(SNES snes, PCCompositeType type)
    836: {

    841:   PetscTryMethod(snes, "SNESMultiblockSetType_C", (SNES, PCCompositeType), (snes, type));
    842:   return(0);
    843: }

    845: /*@
    846:   SNESMultiblockSetBlockSize - Sets the block size for structured mesh block division. If not set the matrix block size is used.

    848:   Logically Collective on SNES

    850:   Input Parameters:
    851: + snes - the solver context
    852: - bs   - the block size

    854:   Level: intermediate

    856: .seealso: SNESMultiblockGetSubSNES(), SNESMULTIBLOCK, SNESMultiblockSetFields()
    857: @*/
    858: PetscErrorCode SNESMultiblockSetBlockSize(SNES snes, PetscInt bs)
    859: {

    865:   PetscTryMethod(snes, "SNESMultiblockSetBlockSize_C", (SNES, PetscInt), (snes,bs));
    866:   return(0);
    867: }

    869: /*@C
    870:   SNESMultiblockGetSubSNES - Gets the SNES contexts for all blocks

    872:   Collective on SNES

    874:   Input Parameter:
    875: . snes - the solver context

    877:   Output Parameters:
    878: + n       - the number of blocks
    879: - subsnes - the array of SNES contexts

    881:   Note:
    882:   After SNESMultiblockGetSubSNES() the array of SNESs MUST be freed by the user
    883:   (not each SNES, just the array that contains them).

    885:   You must call SNESSetUp() before calling SNESMultiblockGetSubSNES().

    887:   Level: advanced

    889: .seealso: SNESMULTIBLOCK
    890: @*/
    891: PetscErrorCode SNESMultiblockGetSubSNES(SNES snes, PetscInt *n, SNES *subsnes[])
    892: {

    898:   PetscUseMethod(snes, "SNESMultiblockGetSubSNES_C", (SNES, PetscInt*, SNES **), (snes, n, subsnes));
    899:   return(0);
    900: }

    902: /*MC
    903:   SNESMULTIBLOCK - Multiblock nonlinear solver that can use overlapping or nonoverlapping blocks, organized
    904:   additively (Jacobi) or multiplicatively (Gauss-Seidel).

    906:   Level: beginner

    908: .seealso:  SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESNEWTONTR, SNESNRICHARDSON
    909: M*/
    910: PETSC_EXTERN PetscErrorCode SNESCreate_Multiblock(SNES snes)
    911: {
    912:   SNES_Multiblock *mb;
    913:   PetscErrorCode  ierr;

    916:   snes->ops->destroy        = SNESDestroy_Multiblock;
    917:   snes->ops->setup          = SNESSetUp_Multiblock;
    918:   snes->ops->setfromoptions = SNESSetFromOptions_Multiblock;
    919:   snes->ops->view           = SNESView_Multiblock;
    920:   snes->ops->solve          = SNESSolve_Multiblock;
    921:   snes->ops->reset          = SNESReset_Multiblock;

    923:   snes->usesksp = PETSC_FALSE;
    924:   snes->usesnpc = PETSC_FALSE;

    926:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

    928:   PetscNewLog(snes,&mb);
    929:   snes->data    = (void*) mb;
    930:   mb->defined   = PETSC_FALSE;
    931:   mb->numBlocks = 0;
    932:   mb->bs        = -1;
    933:   mb->type      = PC_COMPOSITE_MULTIPLICATIVE;

    935:   /* We attach functions so that they can be called on another PC without crashing the program */
    936:   PetscObjectComposeFunction((PetscObject) snes, "SNESMultiblockSetFields_C",    SNESMultiblockSetFields_Default);
    937:   PetscObjectComposeFunction((PetscObject) snes, "SNESMultiblockSetIS_C",        SNESMultiblockSetIS_Default);
    938:   PetscObjectComposeFunction((PetscObject) snes, "SNESMultiblockSetType_C",      SNESMultiblockSetType_Default);
    939:   PetscObjectComposeFunction((PetscObject) snes, "SNESMultiblockSetBlockSize_C", SNESMultiblockSetBlockSize_Default);
    940:   PetscObjectComposeFunction((PetscObject) snes, "SNESMultiblockGetSubSNES_C",   SNESMultiblockGetSubSNES_Default);
    941:   return(0);
    942: }