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: static PetscErrorCode SNESReset_Multiblock(SNES snes)
     27: {
     28:   SNES_Multiblock *mb     = (SNES_Multiblock *)snes->data;
     29:   BlockDesc        blocks = mb->blocks, next;

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

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

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

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

     58:   PetscFunctionBegin;
     59:   PetscCall(SNESReset_Multiblock(snes));
     60:   while (blocks) {
     61:     next = blocks->next;
     62:     PetscCall(SNESDestroy(&blocks->snes));
     63:     PetscCall(PetscFree(blocks->name));
     64:     PetscCall(PetscFree(blocks->fields));
     65:     PetscCall(PetscFree(blocks));
     66:     blocks = next;
     67:   }
     68:   PetscCall(PetscFree(snes->data));
     69:   PetscFunctionReturn(PETSC_SUCCESS);
     70: }

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

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

    102: static PetscErrorCode SNESMultiblockSetDefaults(SNES snes)
    103: {
    104:   SNES_Multiblock *mb     = (SNES_Multiblock *)snes->data;
    105:   BlockDesc        blocks = mb->blocks;
    106:   PetscInt         i;

    108:   PetscFunctionBegin;
    109:   if (!blocks) {
    110:     if (snes->dm) {
    111:       PetscBool dmcomposite;

    113:       PetscCall(PetscObjectTypeCompare((PetscObject)snes->dm, DMCOMPOSITE, &dmcomposite));
    114:       if (dmcomposite) {
    115:         PetscInt nDM;
    116:         IS      *fields;

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

    124:           PetscCall(PetscSNPrintf(name, sizeof(name), "%" PetscInt_FMT, i));
    125:           PetscCall(SNESMultiblockSetIS(snes, name, fields[i]));
    126:           PetscCall(ISDestroy(&fields[i]));
    127:         }
    128:         PetscCall(PetscFree(fields));
    129:       }
    130:     } else {
    131:       PetscBool flg    = PETSC_FALSE;
    132:       PetscBool stokes = PETSC_FALSE;

    134:       if (mb->bs <= 0) {
    135:         if (snes->jacobian_pre) {
    136:           PetscCall(MatGetBlockSize(snes->jacobian_pre, &mb->bs));
    137:         } else mb->bs = 1;
    138:       }

    140:       PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)snes)->prefix, "-snes_multiblock_default", &flg, NULL));
    141:       PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)snes)->prefix, "-snes_multiblock_detect_saddle_point", &stokes, NULL));
    142:       if (stokes) {
    143:         IS       zerodiags, rest;
    144:         PetscInt nmin, nmax;

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

    165:             PetscCall(PetscSNPrintf(name, sizeof(name), "%" PetscInt_FMT, i));
    166:             PetscCall(SNESMultiblockSetFields(snes, name, 1, &i));
    167:           }
    168:           mb->defaultblocks = PETSC_TRUE;
    169:         }
    170:       }
    171:     }
    172:   } else if (mb->numBlocks == 1) {
    173:     if (blocks->is) {
    174:       IS       is2;
    175:       PetscInt nmin, nmax;

    177:       PetscCall(MatGetOwnershipRange(snes->jacobian_pre, &nmin, &nmax));
    178:       PetscCall(ISComplement(blocks->is, nmin, nmax, &is2));
    179:       PetscCall(SNESMultiblockSetIS(snes, "1", is2));
    180:       PetscCall(ISDestroy(&is2));
    181:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least two sets of fields to SNES multiblock");
    182:   }
    183:   PetscCheck(mb->numBlocks >= 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unhandled case, must have at least two blocks");
    184:   PetscFunctionReturn(PETSC_SUCCESS);
    185: }

    187: static PetscErrorCode SNESSetUp_Multiblock(SNES snes)
    188: {
    189:   SNES_Multiblock *mb = (SNES_Multiblock *)snes->data;
    190:   BlockDesc        blocks;
    191:   PetscInt         i, numBlocks;

    193:   PetscFunctionBegin;
    194:   PetscCall(SNESMultiblockSetDefaults(snes));
    195:   numBlocks = mb->numBlocks;
    196:   blocks    = mb->blocks;

    198:   /* Create ISs */
    199:   if (!mb->issetup) {
    200:     PetscInt  ccsize, rstart, rend, nslots, bs;
    201:     PetscBool sorted;

    203:     mb->issetup = PETSC_TRUE;
    204:     bs          = mb->bs;
    205:     PetscCall(MatGetOwnershipRange(snes->jacobian_pre, &rstart, &rend));
    206:     PetscCall(MatGetLocalSize(snes->jacobian_pre, NULL, &ccsize));
    207:     nslots = (rend - rstart) / bs;
    208:     for (i = 0; i < numBlocks; ++i) {
    209:       if (mb->defaultblocks) {
    210:         PetscCall(ISCreateStride(PetscObjectComm((PetscObject)snes), nslots, rstart + i, numBlocks, &blocks->is));
    211:       } else if (!blocks->is) {
    212:         if (blocks->nfields > 1) {
    213:           PetscInt *ii, j, k, nfields = blocks->nfields, *fields = blocks->fields;

    215:           PetscCall(PetscMalloc1(nfields * nslots, &ii));
    216:           for (j = 0; j < nslots; ++j) {
    217:             for (k = 0; k < nfields; ++k) ii[nfields * j + k] = rstart + bs * j + fields[k];
    218:           }
    219:           PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)snes), nslots * nfields, ii, PETSC_OWN_POINTER, &blocks->is));
    220:         } else {
    221:           PetscCall(ISCreateStride(PetscObjectComm((PetscObject)snes), nslots, rstart + blocks->fields[0], bs, &blocks->is));
    222:         }
    223:       }
    224:       PetscCall(ISSorted(blocks->is, &sorted));
    225:       PetscCheck(sorted, PETSC_COMM_SELF, PETSC_ERR_USER, "Fields must be sorted when creating split");
    226:       blocks = blocks->next;
    227:     }
    228:   }

    230: #if 0
    231:   /* Create matrices */
    232:   ilink = jac->head;
    233:   if (!jac->pmat) {
    234:     PetscCall(PetscMalloc1(nsplit,&jac->pmat));
    235:     for (i=0; i<nsplit; i++) {
    236:       PetscCall(MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is,MAT_INITIAL_MATRIX,&jac->pmat[i]));
    237:       ilink = ilink->next;
    238:     }
    239:   } else {
    240:     for (i=0; i<nsplit; i++) {
    241:       PetscCall(MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is,MAT_REUSE_MATRIX,&jac->pmat[i]));
    242:       ilink = ilink->next;
    243:     }
    244:   }
    245:   if (jac->realdiagonal) {
    246:     ilink = jac->head;
    247:     if (!jac->mat) {
    248:       PetscCall(PetscMalloc1(nsplit,&jac->mat));
    249:       for (i=0; i<nsplit; i++) {
    250:         PetscCall(MatCreateSubMatrix(pc->mat,ilink->is,ilink->is,MAT_INITIAL_MATRIX,&jac->mat[i]));
    251:         ilink = ilink->next;
    252:       }
    253:     } else {
    254:       for (i=0; i<nsplit; i++) {
    255:         if (jac->mat[i]) PetscCall(MatCreateSubMatrix(pc->mat,ilink->is,ilink->is,MAT_REUSE_MATRIX,&jac->mat[i]));
    256:         ilink = ilink->next;
    257:       }
    258:     }
    259:   } else jac->mat = jac->pmat;
    260: #endif

    262: #if 0
    263:   if (jac->type != PC_COMPOSITE_ADDITIVE  && jac->type != PC_COMPOSITE_SCHUR) {
    264:     /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */
    265:     ilink = jac->head;
    266:     if (!jac->Afield) {
    267:       PetscCall(PetscMalloc1(nsplit,&jac->Afield));
    268:       for (i=0; i<nsplit; i++) {
    269:         PetscCall(MatCreateSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]));
    270:         ilink = ilink->next;
    271:       }
    272:     } else {
    273:       for (i=0; i<nsplit; i++) {
    274:         PetscCall(MatCreateSubMatrix(pc->mat,ilink->is,NULL,MAT_REUSE_MATRIX,&jac->Afield[i]));
    275:         ilink = ilink->next;
    276:       }
    277:     }
    278:   }
    279: #endif

    281:   if (mb->type == PC_COMPOSITE_SCHUR) {
    282: #if 0
    283:     IS       ccis;
    284:     PetscInt rstart,rend;
    285:     PetscCheck(nsplit == 2,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields");

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

    290:     /* need to handle case when one is resetting up the preconditioner */
    291:     if (jac->schur) {
    292:       ilink = jac->head;
    293:       PetscCall(ISComplement(ilink->is,rstart,rend,&ccis));
    294:       PetscCall(MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B));
    295:       PetscCall(ISDestroy(&ccis));
    296:       ilink = ilink->next;
    297:       PetscCall(ISComplement(ilink->is,rstart,rend,&ccis));
    298:       PetscCall(MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C));
    299:       PetscCall(ISDestroy(&ccis));
    300:       PetscCall(MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->pmat[1]));
    301:       PetscCall(KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),pc->flag));

    303:     } else {
    304:       KSP  ksp;
    305:       char schurprefix[256];

    307:       /* extract the A01 and A10 matrices */
    308:       ilink = jac->head;
    309:       PetscCall(ISComplement(ilink->is,rstart,rend,&ccis));
    310:       PetscCall(MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B));
    311:       PetscCall(ISDestroy(&ccis));
    312:       ilink = ilink->next;
    313:       PetscCall(ISComplement(ilink->is,rstart,rend,&ccis));
    314:       PetscCall(MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C));
    315:       PetscCall(ISDestroy(&ccis));
    316:       /* Use mat[0] (diagonal block of the real matrix) preconditioned by pmat[0] */
    317:       PetscCall(MatCreateSchurComplement(jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1],&jac->schur));
    318:       /* set tabbing and options prefix of KSP inside the MatSchur */
    319:       PetscCall(MatSchurComplementGetKSP(jac->schur,&ksp));
    320:       PetscCall(PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,2));
    321:       PetscCall(PetscSNPrintf(schurprefix,sizeof(schurprefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",jac->head->splitname));
    322:       PetscCall(KSPSetOptionsPrefix(ksp,schurprefix));
    323:       PetscCall(MatSetFromOptions(jac->schur));

    325:       PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur));
    326:       PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1));
    327:       PetscCall(KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac)));
    328:       if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) {
    329:         PC pc;
    330:         PetscCall(KSPGetPC(jac->kspschur,&pc));
    331:         PetscCall(PCSetType(pc,PCNONE));
    332:         /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */
    333:       }
    334:       PetscCall(PetscSNPrintf(schurprefix,sizeof(schurprefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname));
    335:       PetscCall(KSPSetOptionsPrefix(jac->kspschur,schurprefix));
    336:       /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
    337:       PetscCall(KSPSetFromOptions(jac->kspschur));

    339:       PetscCall(PetscMalloc2(2,&jac->x,2,&jac->y));
    340:       PetscCall(MatCreateVecs(jac->pmat[0],&jac->x[0],&jac->y[0]));
    341:       PetscCall(MatCreateVecs(jac->pmat[1],&jac->x[1],&jac->y[1]));
    342:       ilink    = jac->head;
    343:       ilink->x = jac->x[0]; ilink->y = jac->y[0];
    344:       ilink    = ilink->next;
    345:       ilink->x = jac->x[1]; ilink->y = jac->y[1];
    346:     }
    347: #endif
    348:   } else {
    349:     /* Set up the individual SNESs */
    350:     blocks = mb->blocks;
    351:     i      = 0;
    352:     while (blocks) {
    353:       /*TODO: Set these correctly */
    354:       /* PetscCall(SNESSetFunction(blocks->snes, blocks->x, func)); */
    355:       /* PetscCall(SNESSetJacobian(blocks->snes, blocks->x, jac)); */
    356:       PetscCall(VecDuplicate(blocks->snes->vec_sol, &blocks->x));
    357:       /* really want setfromoptions called in SNESSetFromOptions_Multiblock(), but it is not ready yet */
    358:       PetscCall(SNESSetFromOptions(blocks->snes));
    359:       PetscCall(SNESSetUp(blocks->snes));
    360:       blocks = blocks->next;
    361:       i++;
    362:     }
    363:   }

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

    369:     blocks = mb->blocks;
    370:     PetscCall(MatCreateVecs(snes->jacobian_pre, &xtmp, NULL));
    371:     while (blocks) {
    372:       PetscCall(VecScatterCreate(xtmp, blocks->is, blocks->x, NULL, &blocks->sctx));
    373:       blocks = blocks->next;
    374:     }
    375:     PetscCall(VecDestroy(&xtmp));
    376:   }
    377:   PetscFunctionReturn(PETSC_SUCCESS);
    378: }

    380: /*
    381:   SNESSetFromOptions_Multiblock - Sets various parameters for the SNESMULTIBLOCK method.

    383:   Input Parameter:
    384: . snes - the SNES context

    386:   Application Interface Routine: SNESSetFromOptions()
    387: */
    388: static PetscErrorCode SNESSetFromOptions_Multiblock(SNES snes, PetscOptionItems *PetscOptionsObject)
    389: {
    390:   SNES_Multiblock *mb = (SNES_Multiblock *)snes->data;
    391:   PCCompositeType  ctype;
    392:   PetscInt         bs;
    393:   PetscBool        flg;

    395:   PetscFunctionBegin;
    396:   PetscOptionsHeadBegin(PetscOptionsObject, "SNES Multiblock options");
    397:   PetscCall(PetscOptionsInt("-snes_multiblock_block_size", "Blocksize that defines number of fields", "PCFieldSplitSetBlockSize", mb->bs, &bs, &flg));
    398:   if (flg) PetscCall(SNESMultiblockSetBlockSize(snes, bs));
    399:   PetscCall(PetscOptionsEnum("-snes_multiblock_type", "Type of composition", "PCFieldSplitSetType", PCCompositeTypes, (PetscEnum)mb->type, (PetscEnum *)&ctype, &flg));
    400:   if (flg) PetscCall(SNESMultiblockSetType(snes, ctype));
    401:   /* Only setup fields once */
    402:   if ((mb->bs > 0) && (mb->numBlocks == 0)) {
    403:     /* only allow user to set fields from command line if bs is already known, otherwise user can set them in SNESMultiblockSetDefaults() */
    404:     PetscCall(SNESMultiblockSetFieldsRuntime_Private(snes));
    405:     if (mb->defined) PetscCall(PetscInfo(snes, "Blocks defined using the options database\n"));
    406:   }
    407:   PetscOptionsHeadEnd();
    408:   PetscFunctionReturn(PETSC_SUCCESS);
    409: }

    411: static PetscErrorCode SNESView_Multiblock(SNES snes, PetscViewer viewer)
    412: {
    413:   SNES_Multiblock *mb     = (SNES_Multiblock *)snes->data;
    414:   BlockDesc        blocks = mb->blocks;
    415:   PetscBool        iascii;

    417:   PetscFunctionBegin;
    418:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
    419:   if (iascii) {
    420:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Multiblock with %s composition: total blocks = %" PetscInt_FMT ", blocksize = %" PetscInt_FMT "\n", PCCompositeTypes[mb->type], mb->numBlocks, mb->bs));
    421:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Solver info for each split is in the following SNES objects:\n"));
    422:     PetscCall(PetscViewerASCIIPushTab(viewer));
    423:     while (blocks) {
    424:       if (blocks->fields) {
    425:         PetscInt j;

    427:         PetscCall(PetscViewerASCIIPrintf(viewer, "  Block %s Fields ", blocks->name));
    428:         PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
    429:         for (j = 0; j < blocks->nfields; ++j) {
    430:           if (j > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ","));
    431:           PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, blocks->fields[j]));
    432:         }
    433:         PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
    434:         PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
    435:       } else {
    436:         PetscCall(PetscViewerASCIIPrintf(viewer, "  Block %s Defined by IS\n", blocks->name));
    437:       }
    438:       PetscCall(SNESView(blocks->snes, viewer));
    439:       blocks = blocks->next;
    440:     }
    441:     PetscCall(PetscViewerASCIIPopTab(viewer));
    442:   }
    443:   PetscFunctionReturn(PETSC_SUCCESS);
    444: }

    446: static PetscErrorCode SNESSolve_Multiblock(SNES snes)
    447: {
    448:   SNES_Multiblock *mb = (SNES_Multiblock *)snes->data;
    449:   Vec              X, Y, F;
    450:   PetscReal        fnorm;
    451:   PetscInt         maxits, i;

    453:   PetscFunctionBegin;
    454:   PetscCheck(!snes->xl && !snes->xu && !snes->ops->computevariablebounds, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);

    456:   snes->reason = SNES_CONVERGED_ITERATING;

    458:   maxits = snes->max_its;        /* maximum number of iterations */
    459:   X      = snes->vec_sol;        /* X^n */
    460:   Y      = snes->vec_sol_update; /* \tilde X */
    461:   F      = snes->vec_func;       /* residual vector */

    463:   PetscCall(VecSetBlockSize(X, mb->bs));
    464:   PetscCall(VecSetBlockSize(Y, mb->bs));
    465:   PetscCall(VecSetBlockSize(F, mb->bs));
    466:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
    467:   snes->iter = 0;
    468:   snes->norm = 0.;
    469:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));

    471:   if (!snes->vec_func_init_set) {
    472:     PetscCall(SNESComputeFunction(snes, X, F));
    473:   } else snes->vec_func_init_set = PETSC_FALSE;

    475:   PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- ||F||  */
    476:   SNESCheckFunctionNorm(snes, fnorm);
    477:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
    478:   snes->norm = fnorm;
    479:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
    480:   PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0));

    482:   /* test convergence */
    483:   PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm));
    484:   PetscCall(SNESMonitor(snes, 0, fnorm));
    485:   if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);

    487:   for (i = 0; i < maxits; i++) {
    488:     /* Call general purpose update function */
    489:     PetscTryTypeMethod(snes, update, snes->iter);
    490:     /* Compute X^{new} from subsolves */
    491:     if (mb->type == PC_COMPOSITE_ADDITIVE) {
    492:       BlockDesc blocks = mb->blocks;

    494:       if (mb->defaultblocks) {
    495:         /*TODO: Make an array of Vecs for this */
    496:         /* PetscCall(VecStrideGatherAll(X, mb->x, INSERT_VALUES)); */
    497:         while (blocks) {
    498:           PetscCall(SNESSolve(blocks->snes, NULL, blocks->x));
    499:           blocks = blocks->next;
    500:         }
    501:         /* PetscCall(VecStrideScatterAll(mb->x, X, INSERT_VALUES)); */
    502:       } else {
    503:         while (blocks) {
    504:           PetscCall(VecScatterBegin(blocks->sctx, X, blocks->x, INSERT_VALUES, SCATTER_FORWARD));
    505:           PetscCall(VecScatterEnd(blocks->sctx, X, blocks->x, INSERT_VALUES, SCATTER_FORWARD));
    506:           PetscCall(SNESSolve(blocks->snes, NULL, blocks->x));
    507:           PetscCall(VecScatterBegin(blocks->sctx, blocks->x, X, INSERT_VALUES, SCATTER_REVERSE));
    508:           PetscCall(VecScatterEnd(blocks->sctx, blocks->x, X, INSERT_VALUES, SCATTER_REVERSE));
    509:           blocks = blocks->next;
    510:         }
    511:       }
    512:     } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "Unsupported or unknown composition %d", (int)mb->type);
    513:     /* Compute F(X^{new}) */
    514:     PetscCall(SNESComputeFunction(snes, X, F));
    515:     PetscCall(VecNorm(F, NORM_2, &fnorm));
    516:     SNESCheckFunctionNorm(snes, fnorm);

    518:     if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
    519:       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
    520:       break;
    521:     }

    523:     /* Monitor convergence */
    524:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
    525:     snes->iter = i + 1;
    526:     snes->norm = fnorm;
    527:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
    528:     PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
    529:     /* Test for convergence */
    530:     PetscCall(SNESConverged(snes, snes->iter, 0.0, 0.0, fnorm));
    531:     PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
    532:     if (snes->reason) break;
    533:   }
    534:   PetscFunctionReturn(PETSC_SUCCESS);
    535: }

    537: static PetscErrorCode SNESMultiblockSetFields_Default(SNES snes, const char name[], PetscInt n, const PetscInt fields[])
    538: {
    539:   SNES_Multiblock *mb = (SNES_Multiblock *)snes->data;
    540:   BlockDesc        newblock, next = mb->blocks;
    541:   char             prefix[128];
    542:   PetscInt         i;

    544:   PetscFunctionBegin;
    545:   if (mb->defined) {
    546:     PetscCall(PetscInfo(snes, "Ignoring new block \"%s\" because the blocks have already been defined\n", name));
    547:     PetscFunctionReturn(PETSC_SUCCESS);
    548:   }
    549:   for (i = 0; i < n; ++i) {
    550:     PetscCheck(fields[i] < mb->bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field %" PetscInt_FMT " requested but only %" PetscInt_FMT " exist", fields[i], mb->bs);
    551:     PetscCheck(fields[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative field %" PetscInt_FMT " requested", fields[i]);
    552:   }
    553:   PetscCall(PetscNew(&newblock));
    554:   if (name) {
    555:     PetscCall(PetscStrallocpy(name, &newblock->name));
    556:   } else {
    557:     PetscInt len = floor(log10(mb->numBlocks)) + 1;

    559:     PetscCall(PetscMalloc1(len + 1, &newblock->name));
    560:     PetscCall(PetscSNPrintf(newblock->name, len, "%" PetscInt_FMT, mb->numBlocks));
    561:   }
    562:   newblock->nfields = n;

    564:   PetscCall(PetscMalloc1(n, &newblock->fields));
    565:   PetscCall(PetscArraycpy(newblock->fields, fields, n));

    567:   newblock->next = NULL;

    569:   PetscCall(SNESCreate(PetscObjectComm((PetscObject)snes), &newblock->snes));
    570:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)newblock->snes, (PetscObject)snes, 1));
    571:   PetscCall(SNESSetType(newblock->snes, SNESNRICHARDSON));
    572:   PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%smultiblock_%s_", ((PetscObject)snes)->prefix ? ((PetscObject)snes)->prefix : "", newblock->name));
    573:   PetscCall(SNESSetOptionsPrefix(newblock->snes, prefix));

    575:   if (!next) {
    576:     mb->blocks         = newblock;
    577:     newblock->previous = NULL;
    578:   } else {
    579:     while (next->next) next = next->next;
    580:     next->next         = newblock;
    581:     newblock->previous = next;
    582:   }
    583:   mb->numBlocks++;
    584:   PetscFunctionReturn(PETSC_SUCCESS);
    585: }

    587: static PetscErrorCode SNESMultiblockSetIS_Default(SNES snes, const char name[], IS is)
    588: {
    589:   SNES_Multiblock *mb = (SNES_Multiblock *)snes->data;
    590:   BlockDesc        newblock, next = mb->blocks;
    591:   char             prefix[128];

    593:   PetscFunctionBegin;
    594:   if (mb->defined) {
    595:     PetscCall(PetscInfo(snes, "Ignoring new block \"%s\" because the blocks have already been defined\n", name));
    596:     PetscFunctionReturn(PETSC_SUCCESS);
    597:   }
    598:   PetscCall(PetscNew(&newblock));
    599:   if (name) {
    600:     PetscCall(PetscStrallocpy(name, &newblock->name));
    601:   } else {
    602:     PetscInt len = floor(log10(mb->numBlocks)) + 1;

    604:     PetscCall(PetscMalloc1(len + 1, &newblock->name));
    605:     PetscCall(PetscSNPrintf(newblock->name, len, "%" PetscInt_FMT, mb->numBlocks));
    606:   }
    607:   newblock->is = is;

    609:   PetscCall(PetscObjectReference((PetscObject)is));

    611:   newblock->next = NULL;

    613:   PetscCall(SNESCreate(PetscObjectComm((PetscObject)snes), &newblock->snes));
    614:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)newblock->snes, (PetscObject)snes, 1));
    615:   PetscCall(SNESSetType(newblock->snes, SNESNRICHARDSON));
    616:   PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%smultiblock_%s_", ((PetscObject)snes)->prefix ? ((PetscObject)snes)->prefix : "", newblock->name));
    617:   PetscCall(SNESSetOptionsPrefix(newblock->snes, prefix));

    619:   if (!next) {
    620:     mb->blocks         = newblock;
    621:     newblock->previous = NULL;
    622:   } else {
    623:     while (next->next) next = next->next;
    624:     next->next         = newblock;
    625:     newblock->previous = next;
    626:   }
    627:   mb->numBlocks++;
    628:   PetscFunctionReturn(PETSC_SUCCESS);
    629: }

    631: static PetscErrorCode SNESMultiblockSetBlockSize_Default(SNES snes, PetscInt bs)
    632: {
    633:   SNES_Multiblock *mb = (SNES_Multiblock *)snes->data;

    635:   PetscFunctionBegin;
    636:   PetscCheck(bs >= 1, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_OUTOFRANGE, "Blocksize must be positive, you gave %" PetscInt_FMT, bs);
    637:   PetscCheck(mb->bs <= 0 || mb->bs == bs, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Cannot change blocksize from %" PetscInt_FMT " to %" PetscInt_FMT " after it has been set", mb->bs, bs);
    638:   mb->bs = bs;
    639:   PetscFunctionReturn(PETSC_SUCCESS);
    640: }

    642: static PetscErrorCode SNESMultiblockGetSubSNES_Default(SNES snes, PetscInt *n, SNES **subsnes)
    643: {
    644:   SNES_Multiblock *mb     = (SNES_Multiblock *)snes->data;
    645:   BlockDesc        blocks = mb->blocks;
    646:   PetscInt         cnt    = 0;

    648:   PetscFunctionBegin;
    649:   PetscCall(PetscMalloc1(mb->numBlocks, subsnes));
    650:   while (blocks) {
    651:     (*subsnes)[cnt++] = blocks->snes;
    652:     blocks            = blocks->next;
    653:   }
    654:   PetscCheck(cnt == mb->numBlocks, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Corrupt SNESMULTIBLOCK object: number of blocks in linked list %" PetscInt_FMT " does not match number in object %" PetscInt_FMT, cnt, mb->numBlocks);

    656:   if (n) *n = mb->numBlocks;
    657:   PetscFunctionReturn(PETSC_SUCCESS);
    658: }

    660: static PetscErrorCode SNESMultiblockSetType_Default(SNES snes, PCCompositeType type)
    661: {
    662:   SNES_Multiblock *mb = (SNES_Multiblock *)snes->data;

    664:   PetscFunctionBegin;
    665:   mb->type = type;
    666:   if (type == PC_COMPOSITE_SCHUR) {
    667: #if 1
    668:     SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "The Schur composite type is not yet supported");
    669: #else
    670:     snes->ops->solve = SNESSolve_Multiblock_Schur;
    671:     snes->ops->view  = SNESView_Multiblock_Schur;

    673:     PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMultiblockGetSubSNES_C", SNESMultiblockGetSubSNES_Schur));
    674:     PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMultiblockSchurPrecondition_C", SNESMultiblockSchurPrecondition_Default));
    675: #endif
    676:   } else {
    677:     snes->ops->solve = SNESSolve_Multiblock;
    678:     snes->ops->view  = SNESView_Multiblock;

    680:     PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMultiblockGetSubSNES_C", SNESMultiblockGetSubSNES_Default));
    681:     PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMultiblockSchurPrecondition_C", NULL));
    682:   }
    683:   PetscFunctionReturn(PETSC_SUCCESS);
    684: }

    686: /*@
    687:   SNESMultiblockSetFields - Sets the fields for one particular block in a `SNESMULTIBLOCK` solver

    689:   Logically Collective

    691:   Input Parameters:
    692: + snes   - the solver
    693: . name   - name of this block, if `NULL` the number of the block is used
    694: . n      - the number of fields in this block
    695: - fields - the fields in this block

    697:   Level: intermediate

    699:   Notes:
    700:   Use `SNESMultiblockSetIS()` to set a completely general set of row indices as a block.

    702:   The `SNESMultiblockSetFields()` is for defining blocks as a group of strided indices, or fields.
    703:   For example, if the vector block size is three then one can define a block as field 0, or
    704:   1 or 2, or field 0,1 or 0,2 or 1,2 which means
    705:   0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x23x56x8.. x12x45x78x....
    706:   where the numbered entries indicate what is in the block.

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

    711: .seealso: [](ch_snes), `SNES`, `SNESMULTIBLOCK`, `SNESMultiblockGetSubSNES()`, `SNESMultiblockSetBlockSize()`, `SNESMultiblockSetIS()`
    712: @*/
    713: PetscErrorCode SNESMultiblockSetFields(SNES snes, const char name[], PetscInt n, const PetscInt *fields)
    714: {
    715:   PetscFunctionBegin;
    717:   PetscAssertPointer(name, 2);
    718:   PetscCheck(n >= 1, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_OUTOFRANGE, "Provided number of fields %" PetscInt_FMT " in split \"%s\" not positive", n, name);
    719:   PetscAssertPointer(fields, 4);
    720:   PetscTryMethod(snes, "SNESMultiblockSetFields_C", (SNES, const char[], PetscInt, const PetscInt *), (snes, name, n, fields));
    721:   PetscFunctionReturn(PETSC_SUCCESS);
    722: }

    724: /*@
    725:   SNESMultiblockSetIS - Sets the global row indices for one particular block in a `SNESMULTIBLOCK` solver

    727:   Logically Collective

    729:   Input Parameters:
    730: + snes - the solver context
    731: . name - name of this block, if `NULL` the number of the block is used
    732: - is   - the index set that defines the global row indices in this block

    734:   Level: intermediate

    736:   Notes:
    737:   Use `SNESMultiblockSetFields()`, for blocks defined by strides.

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

    742: .seealso: [](ch_snes), `SNES`, `SNESMULTIBLOCK`, `SNESMultiblockGetSubSNES()`, `SNESMultiblockSetBlockSize()`, `SNESMultiblockSetFields()`
    743: @*/
    744: PetscErrorCode SNESMultiblockSetIS(SNES snes, const char name[], IS is)
    745: {
    746:   PetscFunctionBegin;
    748:   PetscAssertPointer(name, 2);
    750:   PetscTryMethod(snes, "SNESMultiblockSetIS_C", (SNES, const char[], IS), (snes, name, is));
    751:   PetscFunctionReturn(PETSC_SUCCESS);
    752: }

    754: /*@
    755:   SNESMultiblockSetType - Sets the type of block combination used for a `SNESMULTIBLOCK` solver

    757:   Logically Collective

    759:   Input Parameters:
    760: + snes - the solver context
    761: - type - `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE` (default), `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`

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

    766:   Level: advanced

    768:   Developer Note:
    769:   This `SNESType` uses `PCCompositeType`, while `SNESCompositeSetType()` uses `SNESCOMPOSITE`, perhaps they should be unified in the future

    771: .seealso: [](ch_snes), `SNES`, `SNESMULTIBLOCK`, `PCCompositeSetType()`, `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`, `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`,
    772:           `PCCompositeType`, `SNESCOMPOSITE`, `SNESCompositeSetType()`
    773: @*/
    774: PetscErrorCode SNESMultiblockSetType(SNES snes, PCCompositeType type)
    775: {
    776:   PetscFunctionBegin;
    778:   PetscTryMethod(snes, "SNESMultiblockSetType_C", (SNES, PCCompositeType), (snes, type));
    779:   PetscFunctionReturn(PETSC_SUCCESS);
    780: }

    782: /*@
    783:   SNESMultiblockSetBlockSize - Sets the block size for structured block division in a `SNESMULTIBLOCK` solver. If not set the matrix block size is used.

    785:   Logically Collective

    787:   Input Parameters:
    788: + snes - the solver context
    789: - bs   - the block size

    791:   Level: intermediate

    793: .seealso: [](ch_snes), `SNES`, `SNESMULTIBLOCK`, `SNESMultiblockGetSubSNES()`, `SNESMultiblockSetFields()`
    794: @*/
    795: PetscErrorCode SNESMultiblockSetBlockSize(SNES snes, PetscInt bs)
    796: {
    797:   PetscFunctionBegin;
    800:   PetscTryMethod(snes, "SNESMultiblockSetBlockSize_C", (SNES, PetscInt), (snes, bs));
    801:   PetscFunctionReturn(PETSC_SUCCESS);
    802: }

    804: /*@C
    805:   SNESMultiblockGetSubSNES - Gets the `SNES` contexts for all blocks in a `SNESMULTIBLOCK` solver.

    807:   Not Collective but each `SNES` obtained is parallel

    809:   Input Parameter:
    810: . snes - the solver context

    812:   Output Parameters:
    813: + n       - the number of blocks
    814: - subsnes - the array of `SNES` contexts

    816:   Level: advanced

    818:   Notes:
    819:   After `SNESMultiblockGetSubSNES()` the array of `SNES`s MUST be freed by the user
    820:   (not each `SNES`, just the array that contains them).

    822:   You must call `SNESSetUp()` before calling `SNESMultiblockGetSubSNES()`.

    824: .seealso: [](ch_snes), `SNES`, `SNESMULTIBLOCK`, `SNESMultiblockSetIS()`, `SNESMultiblockSetFields()`
    825: @*/
    826: PetscErrorCode SNESMultiblockGetSubSNES(SNES snes, PetscInt *n, SNES *subsnes[])
    827: {
    828:   PetscFunctionBegin;
    830:   if (n) PetscAssertPointer(n, 2);
    831:   PetscUseMethod(snes, "SNESMultiblockGetSubSNES_C", (SNES, PetscInt *, SNES **), (snes, n, subsnes));
    832:   PetscFunctionReturn(PETSC_SUCCESS);
    833: }

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

    839:   Level: beginner

    841:   Note:
    842:   This is much like  `PCASM`, `PCBJACOBI`, and `PCFIELDSPLIT` are for linear problems.

    844: .seealso: [](ch_snes), `SNES`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESNEWTONTR`, `SNESNRICHARDSON`, `SNESMultiblockSetType()`,
    845:           `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`, `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `SNESMultiblockSetBlockSize()`,
    846:           `SNESMultiblockGetBlockSize()`, `SNESMultiblockSetFields()`, `SNESMultiblockSetIS()`, `SNESMultiblockGetSubSNES()`, `PCASM`, `PCBJACOBI`
    847: M*/
    848: PETSC_EXTERN PetscErrorCode SNESCreate_Multiblock(SNES snes)
    849: {
    850:   SNES_Multiblock *mb;

    852:   PetscFunctionBegin;
    853:   snes->ops->destroy        = SNESDestroy_Multiblock;
    854:   snes->ops->setup          = SNESSetUp_Multiblock;
    855:   snes->ops->setfromoptions = SNESSetFromOptions_Multiblock;
    856:   snes->ops->view           = SNESView_Multiblock;
    857:   snes->ops->solve          = SNESSolve_Multiblock;
    858:   snes->ops->reset          = SNESReset_Multiblock;

    860:   snes->usesksp = PETSC_FALSE;
    861:   snes->usesnpc = PETSC_FALSE;

    863:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

    865:   PetscCall(SNESParametersInitialize(snes));

    867:   PetscCall(PetscNew(&mb));
    868:   snes->data    = (void *)mb;
    869:   mb->defined   = PETSC_FALSE;
    870:   mb->numBlocks = 0;
    871:   mb->bs        = -1;
    872:   mb->type      = PC_COMPOSITE_MULTIPLICATIVE;

    874:   /* We attach functions so that they can be called on another PC without crashing the program */
    875:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMultiblockSetFields_C", SNESMultiblockSetFields_Default));
    876:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMultiblockSetIS_C", SNESMultiblockSetIS_Default));
    877:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMultiblockSetType_C", SNESMultiblockSetType_Default));
    878:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMultiblockSetBlockSize_C", SNESMultiblockSetBlockSize_Default));
    879:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESMultiblockGetSubSNES_C", SNESMultiblockGetSubSNES_Default));
    880:   PetscFunctionReturn(PETSC_SUCCESS);
    881: }