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;

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

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

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

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

 57:   SNESReset_Multiblock(snes);
 58:   while (blocks) {
 59:     next = blocks->next;
 60:     SNESDestroy(&blocks->snes);
 61:     PetscFree(blocks->name);
 62:     PetscFree(blocks->fields);
 63:     PetscFree(blocks);
 64:     blocks = next;
 65:   }
 66:   PetscFree(snes->data);
 67:   return 0;
 68: }

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

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

 99: static PetscErrorCode SNESMultiblockSetDefaults(SNES snes)
100: {
101:   SNES_Multiblock *mb     = (SNES_Multiblock *)snes->data;
102:   BlockDesc        blocks = mb->blocks;
103:   PetscInt         i;

105:   if (!blocks) {
106:     if (snes->dm) {
107:       PetscBool dmcomposite;

109:       PetscObjectTypeCompare((PetscObject)snes->dm, DMCOMPOSITE, &dmcomposite);
110:       if (dmcomposite) {
111:         PetscInt nDM;
112:         IS      *fields;

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

120:           PetscSNPrintf(name, sizeof(name), "%" PetscInt_FMT, i);
121:           SNESMultiblockSetIS(snes, name, fields[i]);
122:           ISDestroy(&fields[i]);
123:         }
124:         PetscFree(fields);
125:       }
126:     } else {
127:       PetscBool flg    = PETSC_FALSE;
128:       PetscBool stokes = PETSC_FALSE;

130:       if (mb->bs <= 0) {
131:         if (snes->jacobian_pre) {
132:           MatGetBlockSize(snes->jacobian_pre, &mb->bs);
133:         } else mb->bs = 1;
134:       }

136:       PetscOptionsGetBool(NULL, ((PetscObject)snes)->prefix, "-snes_multiblock_default", &flg, NULL);
137:       PetscOptionsGetBool(NULL, ((PetscObject)snes)->prefix, "-snes_multiblock_detect_saddle_point", &stokes, NULL);
138:       if (stokes) {
139:         IS       zerodiags, rest;
140:         PetscInt nmin, nmax;

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

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

173:       MatGetOwnershipRange(snes->jacobian_pre, &nmin, &nmax);
174:       ISComplement(blocks->is, nmin, nmax, &is2);
175:       SNESMultiblockSetIS(snes, "1", is2);
176:       ISDestroy(&is2);
177:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least two sets of fields to SNES multiblock");
178:   }
180:   return 0;
181: }

183: /*
184:    SNESSetUp_Multiblock - Sets up the internal data structures for the later use
185:    of the SNESMULTIBLOCK nonlinear solver.

187:    Input Parameters:
188: +  snes - the SNES context
189: -  x - the solution vector

191:    Application Interface Routine: SNESSetUp()
192: */
193: PetscErrorCode SNESSetUp_Multiblock(SNES snes)
194: {
195:   SNES_Multiblock *mb = (SNES_Multiblock *)snes->data;
196:   BlockDesc        blocks;
197:   PetscInt         i, numBlocks;

199:   SNESMultiblockSetDefaults(snes);
200:   numBlocks = mb->numBlocks;
201:   blocks    = mb->blocks;

203:   /* Create ISs */
204:   if (!mb->issetup) {
205:     PetscInt  ccsize, rstart, rend, nslots, bs;
206:     PetscBool sorted;

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

220:           PetscMalloc1(nfields * nslots, &ii);
221:           for (j = 0; j < nslots; ++j) {
222:             for (k = 0; k < nfields; ++k) ii[nfields * j + k] = rstart + bs * j + fields[k];
223:           }
224:           ISCreateGeneral(PetscObjectComm((PetscObject)snes), nslots * nfields, ii, PETSC_OWN_POINTER, &blocks->is);
225:         } else {
226:           ISCreateStride(PetscObjectComm((PetscObject)snes), nslots, rstart + blocks->fields[0], bs, &blocks->is);
227:         }
228:       }
229:       ISSorted(blocks->is, &sorted);
231:       blocks = blocks->next;
232:     }
233:   }

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

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

286:   if (mb->type == PC_COMPOSITE_SCHUR) {
287: #if 0
288:     IS       ccis;
289:     PetscInt rstart,rend;

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

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

308:     } else {
309:       KSP  ksp;
310:       char schurprefix[256];

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

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

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

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

374:     blocks = mb->blocks;
375:     MatCreateVecs(snes->jacobian_pre, &xtmp, NULL);
376:     while (blocks) {
377:       VecScatterCreate(xtmp, blocks->is, blocks->x, NULL, &blocks->sctx);
378:       blocks = blocks->next;
379:     }
380:     VecDestroy(&xtmp);
381:   }
382:   return 0;
383: }

385: /*
386:   SNESSetFromOptions_Multiblock - Sets various parameters for the SNESMULTIBLOCK method.

388:   Input Parameter:
389: . snes - the SNES context

391:   Application Interface Routine: SNESSetFromOptions()
392: */
393: static PetscErrorCode SNESSetFromOptions_Multiblock(SNES snes, PetscOptionItems *PetscOptionsObject)
394: {
395:   SNES_Multiblock *mb = (SNES_Multiblock *)snes->data;
396:   PCCompositeType  ctype;
397:   PetscInt         bs;
398:   PetscBool        flg;

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

415: /*
416:   SNESView_Multiblock - Prints info from the SNESMULTIBLOCK data structure.

418:   Input Parameters:
419: + SNES - the SNES context
420: - viewer - visualization context

422:   Application Interface Routine: SNESView()
423: */
424: static PetscErrorCode SNESView_Multiblock(SNES snes, PetscViewer viewer)
425: {
426:   SNES_Multiblock *mb     = (SNES_Multiblock *)snes->data;
427:   BlockDesc        blocks = mb->blocks;
428:   PetscBool        iascii;

430:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
431:   if (iascii) {
432:     PetscViewerASCIIPrintf(viewer, "  Multiblock with %s composition: total blocks = %" PetscInt_FMT ", blocksize = %" PetscInt_FMT "\n", PCCompositeTypes[mb->type], mb->numBlocks, mb->bs);
433:     PetscViewerASCIIPrintf(viewer, "  Solver info for each split is in the following SNES objects:\n");
434:     PetscViewerASCIIPushTab(viewer);
435:     while (blocks) {
436:       if (blocks->fields) {
437:         PetscInt j;

439:         PetscViewerASCIIPrintf(viewer, "  Block %s Fields ", blocks->name);
440:         PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
441:         for (j = 0; j < blocks->nfields; ++j) {
442:           if (j > 0) PetscViewerASCIIPrintf(viewer, ",");
443:           PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, blocks->fields[j]);
444:         }
445:         PetscViewerASCIIPrintf(viewer, "\n");
446:         PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
447:       } else {
448:         PetscViewerASCIIPrintf(viewer, "  Block %s Defined by IS\n", blocks->name);
449:       }
450:       SNESView(blocks->snes, viewer);
451:       blocks = blocks->next;
452:     }
453:     PetscViewerASCIIPopTab(viewer);
454:   }
455:   return 0;
456: }

458: /*
459:   SNESSolve_Multiblock - Solves a nonlinear system with the Multiblock method.

461:   Input Parameters:
462: . snes - the SNES context

464:   Output Parameter:
465: . outits - number of iterations until termination

467:   Application Interface Routine: SNESSolve()
468: */
469: PetscErrorCode SNESSolve_Multiblock(SNES snes)
470: {
471:   SNES_Multiblock *mb = (SNES_Multiblock *)snes->data;
472:   Vec              X, Y, F;
473:   PetscReal        fnorm;
474:   PetscInt         maxits, i;


478:   snes->reason = SNES_CONVERGED_ITERATING;

480:   maxits = snes->max_its;        /* maximum number of iterations */
481:   X      = snes->vec_sol;        /* X^n */
482:   Y      = snes->vec_sol_update; /* \tilde X */
483:   F      = snes->vec_func;       /* residual vector */

485:   VecSetBlockSize(X, mb->bs);
486:   VecSetBlockSize(Y, mb->bs);
487:   VecSetBlockSize(F, mb->bs);
488:   PetscObjectSAWsTakeAccess((PetscObject)snes);
489:   snes->iter = 0;
490:   snes->norm = 0.;
491:   PetscObjectSAWsGrantAccess((PetscObject)snes);

493:   if (!snes->vec_func_init_set) {
494:     SNESComputeFunction(snes, X, F);
495:   } else snes->vec_func_init_set = PETSC_FALSE;

497:   VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
498:   SNESCheckFunctionNorm(snes, fnorm);
499:   PetscObjectSAWsTakeAccess((PetscObject)snes);
500:   snes->norm = fnorm;
501:   PetscObjectSAWsGrantAccess((PetscObject)snes);
502:   SNESLogConvergenceHistory(snes, fnorm, 0);
503:   SNESMonitor(snes, 0, fnorm);

505:   /* test convergence */
506:   PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
507:   if (snes->reason) return 0;

509:   for (i = 0; i < maxits; i++) {
510:     /* Call general purpose update function */
511:     PetscTryTypeMethod(snes, update, snes->iter);
512:     /* Compute X^{new} from subsolves */
513:     if (mb->type == PC_COMPOSITE_ADDITIVE) {
514:       BlockDesc blocks = mb->blocks;

516:       if (mb->defaultblocks) {
517:         /*TODO: Make an array of Vecs for this */
518:         /* VecStrideGatherAll(X, mb->x, INSERT_VALUES); */
519:         while (blocks) {
520:           SNESSolve(blocks->snes, NULL, blocks->x);
521:           blocks = blocks->next;
522:         }
523:         /* VecStrideScatterAll(mb->x, X, INSERT_VALUES); */
524:       } else {
525:         while (blocks) {
526:           VecScatterBegin(blocks->sctx, X, blocks->x, INSERT_VALUES, SCATTER_FORWARD);
527:           VecScatterEnd(blocks->sctx, X, blocks->x, INSERT_VALUES, SCATTER_FORWARD);
528:           SNESSolve(blocks->snes, NULL, blocks->x);
529:           VecScatterBegin(blocks->sctx, blocks->x, X, INSERT_VALUES, SCATTER_REVERSE);
530:           VecScatterEnd(blocks->sctx, blocks->x, X, INSERT_VALUES, SCATTER_REVERSE);
531:           blocks = blocks->next;
532:         }
533:       }
534:     } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "Unsupported or unknown composition %d", (int)mb->type);
535:     /* Compute F(X^{new}) */
536:     SNESComputeFunction(snes, X, F);
537:     VecNorm(F, NORM_2, &fnorm);
538:     SNESCheckFunctionNorm(snes, fnorm);

540:     if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
541:       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
542:       break;
543:     }

545:     /* Monitor convergence */
546:     PetscObjectSAWsTakeAccess((PetscObject)snes);
547:     snes->iter = i + 1;
548:     snes->norm = fnorm;
549:     PetscObjectSAWsGrantAccess((PetscObject)snes);
550:     SNESLogConvergenceHistory(snes, snes->norm, 0);
551:     SNESMonitor(snes, snes->iter, snes->norm);
552:     /* Test for convergence */
553:     PetscUseTypeMethod(snes, converged, snes->iter, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
554:     if (snes->reason) break;
555:   }
556:   if (i == maxits) {
557:     PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", maxits);
558:     if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
559:   }
560:   return 0;
561: }

563: PetscErrorCode SNESMultiblockSetFields_Default(SNES snes, const char name[], PetscInt n, const PetscInt fields[])
564: {
565:   SNES_Multiblock *mb = (SNES_Multiblock *)snes->data;
566:   BlockDesc        newblock, next = mb->blocks;
567:   char             prefix[128];
568:   PetscInt         i;

570:   if (mb->defined) {
571:     PetscInfo(snes, "Ignoring new block \"%s\" because the blocks have already been defined\n", name);
572:     return 0;
573:   }
574:   for (i = 0; i < n; ++i) {
577:   }
578:   PetscNew(&newblock);
579:   if (name) {
580:     PetscStrallocpy(name, &newblock->name);
581:   } else {
582:     PetscInt len = floor(log10(mb->numBlocks)) + 1;

584:     PetscMalloc1(len + 1, &newblock->name);
585:     PetscSNPrintf(newblock->name, len, "%" PetscInt_FMT, mb->numBlocks);
586:   }
587:   newblock->nfields = n;

589:   PetscMalloc1(n, &newblock->fields);
590:   PetscArraycpy(newblock->fields, fields, n);

592:   newblock->next = NULL;

594:   SNESCreate(PetscObjectComm((PetscObject)snes), &newblock->snes);
595:   PetscObjectIncrementTabLevel((PetscObject)newblock->snes, (PetscObject)snes, 1);
596:   SNESSetType(newblock->snes, SNESNRICHARDSON);
597:   PetscSNPrintf(prefix, sizeof(prefix), "%smultiblock_%s_", ((PetscObject)snes)->prefix ? ((PetscObject)snes)->prefix : "", newblock->name);
598:   SNESSetOptionsPrefix(newblock->snes, prefix);

600:   if (!next) {
601:     mb->blocks         = newblock;
602:     newblock->previous = NULL;
603:   } else {
604:     while (next->next) next = next->next;
605:     next->next         = newblock;
606:     newblock->previous = next;
607:   }
608:   mb->numBlocks++;
609:   return 0;
610: }

612: PetscErrorCode SNESMultiblockSetIS_Default(SNES snes, const char name[], IS is)
613: {
614:   SNES_Multiblock *mb = (SNES_Multiblock *)snes->data;
615:   BlockDesc        newblock, next = mb->blocks;
616:   char             prefix[128];

618:   if (mb->defined) {
619:     PetscInfo(snes, "Ignoring new block \"%s\" because the blocks have already been defined\n", name);
620:     return 0;
621:   }
622:   PetscNew(&newblock);
623:   if (name) {
624:     PetscStrallocpy(name, &newblock->name);
625:   } else {
626:     PetscInt len = floor(log10(mb->numBlocks)) + 1;

628:     PetscMalloc1(len + 1, &newblock->name);
629:     PetscSNPrintf(newblock->name, len, "%" PetscInt_FMT, mb->numBlocks);
630:   }
631:   newblock->is = is;

633:   PetscObjectReference((PetscObject)is);

635:   newblock->next = NULL;

637:   SNESCreate(PetscObjectComm((PetscObject)snes), &newblock->snes);
638:   PetscObjectIncrementTabLevel((PetscObject)newblock->snes, (PetscObject)snes, 1);
639:   SNESSetType(newblock->snes, SNESNRICHARDSON);
640:   PetscSNPrintf(prefix, sizeof(prefix), "%smultiblock_%s_", ((PetscObject)snes)->prefix ? ((PetscObject)snes)->prefix : "", newblock->name);
641:   SNESSetOptionsPrefix(newblock->snes, prefix);

643:   if (!next) {
644:     mb->blocks         = newblock;
645:     newblock->previous = NULL;
646:   } else {
647:     while (next->next) next = next->next;
648:     next->next         = newblock;
649:     newblock->previous = next;
650:   }
651:   mb->numBlocks++;
652:   return 0;
653: }

655: PetscErrorCode SNESMultiblockSetBlockSize_Default(SNES snes, PetscInt bs)
656: {
657:   SNES_Multiblock *mb = (SNES_Multiblock *)snes->data;

661:   mb->bs = bs;
662:   return 0;
663: }

665: PetscErrorCode SNESMultiblockGetSubSNES_Default(SNES snes, PetscInt *n, SNES **subsnes)
666: {
667:   SNES_Multiblock *mb     = (SNES_Multiblock *)snes->data;
668:   BlockDesc        blocks = mb->blocks;
669:   PetscInt         cnt    = 0;

671:   PetscMalloc1(mb->numBlocks, subsnes);
672:   while (blocks) {
673:     (*subsnes)[cnt++] = blocks->snes;
674:     blocks            = blocks->next;
675:   }

678:   if (n) *n = mb->numBlocks;
679:   return 0;
680: }

682: PetscErrorCode SNESMultiblockSetType_Default(SNES snes, PCCompositeType type)
683: {
684:   SNES_Multiblock *mb = (SNES_Multiblock *)snes->data;

686:   mb->type = type;
687:   if (type == PC_COMPOSITE_SCHUR) {
688: #if 1
689:     SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "The Schur composite type is not yet supported");
690: #else
691:     snes->ops->solve = SNESSolve_Multiblock_Schur;
692:     snes->ops->view  = SNESView_Multiblock_Schur;

694:     PetscObjectComposeFunction((PetscObject)snes, "SNESMultiblockGetSubSNES_C", SNESMultiblockGetSubSNES_Schur);
695:     PetscObjectComposeFunction((PetscObject)snes, "SNESMultiblockSchurPrecondition_C", SNESMultiblockSchurPrecondition_Default);
696: #endif
697:   } else {
698:     snes->ops->solve = SNESSolve_Multiblock;
699:     snes->ops->view  = SNESView_Multiblock;

701:     PetscObjectComposeFunction((PetscObject)snes, "SNESMultiblockGetSubSNES_C", SNESMultiblockGetSubSNES_Default);
702:     PetscObjectComposeFunction((PetscObject)snes, "SNESMultiblockSchurPrecondition_C", 0);
703:   }
704:   return 0;
705: }

707: /*@
708:   SNESMultiblockSetFields - Sets the fields for one particular block in a `SNESMULTBLOCK` solver

710:   Logically Collective

712:   Input Parameters:
713: + snes   - the solver
714: . name   - name of this block, if NULL the number of the block is used
715: . n      - the number of fields in this block
716: - fields - the fields in this block

718:   Level: intermediate

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

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

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

732: .seealso: `SNESMULTBLOCK`, `SNESMultiblockGetSubSNES()`, `SNESMULTIBLOCK`, `SNESMultiblockSetBlockSize()`, `SNESMultiblockSetIS()`
733: @*/
734: PetscErrorCode SNESMultiblockSetFields(SNES snes, const char name[], PetscInt n, const PetscInt *fields)
735: {
740:   PetscTryMethod(snes, "SNESMultiblockSetFields_C", (SNES, const char[], PetscInt, const PetscInt *), (snes, name, n, fields));
741:   return 0;
742: }

744: /*@
745:   SNESMultiblockSetIS - Sets the global row indices for one particular block in a `SNESMULTBLOCK` solver

747:   Logically Collective

749:   Input Parameters:
750: + snes - the solver context
751: . name - name of this block, if NULL the number of the block is used
752: - is   - the index set that defines the global row indices in this block

754:   Notes:
755:   Use `SNESMultiblockSetFields()`, for blocks defined by strides.

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

760:   Level: intermediate

762: .seealso: `SNESMULTBLOCK`, `SNESMultiblockGetSubSNES()`, `SNESMULTIBLOCK`, `SNESMultiblockSetBlockSize()`
763: @*/
764: PetscErrorCode SNESMultiblockSetIS(SNES snes, const char name[], IS is)
765: {
769:   PetscTryMethod(snes, "SNESMultiblockSetIS_C", (SNES, const char[], IS), (snes, name, is));
770:   return 0;
771: }

773: /*@
774:   SNESMultiblockSetType - Sets the type of block combination used for a `SNESMULTBLOCK` solver

776:   Logically Collective

778:   Input Parameters:
779: + snes - the solver context
780: - type - `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE` (default), `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`

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

785:   Level: advanced

787: .seealso: `SNESMULTBLOCK`, `PCCompositeSetType()`, `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`, `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`
788: @*/
789: PetscErrorCode SNESMultiblockSetType(SNES snes, PCCompositeType type)
790: {
792:   PetscTryMethod(snes, "SNESMultiblockSetType_C", (SNES, PCCompositeType), (snes, type));
793:   return 0;
794: }

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

799:   Logically Collective

801:   Input Parameters:
802: + snes - the solver context
803: - bs   - the block size

805:   Level: intermediate

807: .seealso: `SNESMULTBLOCK`, `SNESMultiblockGetSubSNES()`, `SNESMULTIBLOCK`, `SNESMultiblockSetFields()`
808: @*/
809: PetscErrorCode SNESMultiblockSetBlockSize(SNES snes, PetscInt bs)
810: {
813:   PetscTryMethod(snes, "SNESMultiblockSetBlockSize_C", (SNES, PetscInt), (snes, bs));
814:   return 0;
815: }

817: /*@C
818:   SNESMultiblockGetSubSNES - Gets the `SNES` contexts for all blocks in a `SNESMULTBLOCK` solver.

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

822:   Input Parameter:
823: . snes - the solver context

825:   Output Parameters:
826: + n       - the number of blocks
827: - subsnes - the array of `SNES` contexts

829:   Note:
830:   After `SNESMultiblockGetSubSNES()` the array of `SNES`s MUST be freed by the user
831:   (not each `SNES`, just the array that contains them).

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

835:   Level: advanced

837: .seealso: `SNESMULTBLOCK`, `SNESMultiblockSetIS()`, `SNESMultiblockSetFields()`
838: @*/
839: PetscErrorCode SNESMultiblockGetSubSNES(SNES snes, PetscInt *n, SNES *subsnes[])
840: {
843:   PetscUseMethod(snes, "SNESMultiblockGetSubSNES_C", (SNES, PetscInt *, SNES **), (snes, n, subsnes));
844:   return 0;
845: }

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

851:   Level: beginner

853: .seealso: `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESNEWTONTR`, `SNESNRICHARDSON`, `SNESMultiblockSetType()`,
854:           `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`, `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`
855: M*/
856: PETSC_EXTERN PetscErrorCode SNESCreate_Multiblock(SNES snes)
857: {
858:   SNES_Multiblock *mb;

860:   snes->ops->destroy        = SNESDestroy_Multiblock;
861:   snes->ops->setup          = SNESSetUp_Multiblock;
862:   snes->ops->setfromoptions = SNESSetFromOptions_Multiblock;
863:   snes->ops->view           = SNESView_Multiblock;
864:   snes->ops->solve          = SNESSolve_Multiblock;
865:   snes->ops->reset          = SNESReset_Multiblock;

867:   snes->usesksp = PETSC_FALSE;
868:   snes->usesnpc = PETSC_FALSE;

870:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

872:   PetscNew(&mb);
873:   snes->data    = (void *)mb;
874:   mb->defined   = PETSC_FALSE;
875:   mb->numBlocks = 0;
876:   mb->bs        = -1;
877:   mb->type      = PC_COMPOSITE_MULTIPLICATIVE;

879:   /* We attach functions so that they can be called on another PC without crashing the program */
880:   PetscObjectComposeFunction((PetscObject)snes, "SNESMultiblockSetFields_C", SNESMultiblockSetFields_Default);
881:   PetscObjectComposeFunction((PetscObject)snes, "SNESMultiblockSetIS_C", SNESMultiblockSetIS_Default);
882:   PetscObjectComposeFunction((PetscObject)snes, "SNESMultiblockSetType_C", SNESMultiblockSetType_Default);
883:   PetscObjectComposeFunction((PetscObject)snes, "SNESMultiblockSetBlockSize_C", SNESMultiblockSetBlockSize_Default);
884:   PetscObjectComposeFunction((PetscObject)snes, "SNESMultiblockGetSubSNES_C", SNESMultiblockGetSubSNES_Default);
885:   return 0;
886: }