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), "%D", i);
 82:     PetscSNPrintf(optionname, sizeof(optionname), "-snes_multiblock_%D_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), "%D", 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), "%D", 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) {
223:               ii[nfields*j + k] = rstart + bs*j + fields[k];
224:             }
225:           }
226:           ISCreateGeneral(PetscObjectComm((PetscObject)snes), nslots*nfields, ii, PETSC_OWN_POINTER, &blocks->is);
227:         } else {
228:           ISCreateStride(PetscObjectComm((PetscObject)snes), nslots, rstart+blocks->fields[0], bs, &blocks->is);
229:         }
230:       }
231:       ISSorted(blocks->is, &sorted);
233:       blocks = blocks->next;
234:     }
235:   }

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

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

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

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

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

310:     } else {
311:       KSP  ksp;
312:       char schurprefix[256];

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

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

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

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

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

388: /*
389:   SNESSetFromOptions_Multiblock - Sets various parameters for the SNESMULTIBLOCK method.

391:   Input Parameter:
392: . snes - the SNES context

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

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

420: /*
421:   SNESView_Multiblock - Prints info from the SNESMULTIBLOCK data structure.

423:   Input Parameters:
424: + SNES - the SNES context
425: - viewer - visualization context

427:   Application Interface Routine: SNESView()
428: */
429: static PetscErrorCode SNESView_Multiblock(SNES snes, PetscViewer viewer)
430: {
431:   SNES_Multiblock *mb    = (SNES_Multiblock*) snes->data;
432:   BlockDesc       blocks = mb->blocks;
433:   PetscBool       iascii;

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

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

465: /*
466:   SNESSolve_Multiblock - Solves a nonlinear system with the Multiblock method.

468:   Input Parameters:
469: . snes - the SNES context

471:   Output Parameter:
472: . outits - number of iterations until termination

474:   Application Interface Routine: SNESSolve()
475: */
476: PetscErrorCode SNESSolve_Multiblock(SNES snes)
477: {
478:   SNES_Multiblock *mb = (SNES_Multiblock*) snes->data;
479:   Vec             X, Y, F;
480:   PetscReal       fnorm;
481:   PetscInt        maxits, i;


485:   snes->reason = SNES_CONVERGED_ITERATING;

487:   maxits = snes->max_its;        /* maximum number of iterations */
488:   X      = snes->vec_sol;        /* X^n */
489:   Y      = snes->vec_sol_update; /* \tilde X */
490:   F      = snes->vec_func;       /* residual vector */

492:   VecSetBlockSize(X, mb->bs);
493:   VecSetBlockSize(Y, mb->bs);
494:   VecSetBlockSize(F, mb->bs);
495:   PetscObjectSAWsTakeAccess((PetscObject)snes);
496:   snes->iter = 0;
497:   snes->norm = 0.;
498:   PetscObjectSAWsGrantAccess((PetscObject)snes);

500:   if (!snes->vec_func_init_set) {
501:     SNESComputeFunction(snes, X, F);
502:   } else snes->vec_func_init_set = PETSC_FALSE;

504:   VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
505:   SNESCheckFunctionNorm(snes,fnorm);
506:   PetscObjectSAWsTakeAccess((PetscObject)snes);
507:   snes->norm = fnorm;
508:   PetscObjectSAWsGrantAccess((PetscObject)snes);
509:   SNESLogConvergenceHistory(snes,fnorm,0);
510:   SNESMonitor(snes,0,fnorm);

512:   /* test convergence */
513:   (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
514:   if (snes->reason) return 0;

516:   for (i = 0; i < maxits; i++) {
517:     /* Call general purpose update function */
518:     if (snes->ops->update) {
519:       (*snes->ops->update)(snes, snes->iter);
520:     }
521:     /* Compute X^{new} from subsolves */
522:     if (mb->type == PC_COMPOSITE_ADDITIVE) {
523:       BlockDesc blocks = mb->blocks;

525:       if (mb->defaultblocks) {
526:         /*TODO: Make an array of Vecs for this */
527:         /* VecStrideGatherAll(X, mb->x, INSERT_VALUES); */
528:         while (blocks) {
529:           SNESSolve(blocks->snes, NULL, blocks->x);
530:           blocks = blocks->next;
531:         }
532:         /* VecStrideScatterAll(mb->x, X, INSERT_VALUES); */
533:       } else {
534:         while (blocks) {
535:           VecScatterBegin(blocks->sctx, X, blocks->x, INSERT_VALUES, SCATTER_FORWARD);
536:           VecScatterEnd(blocks->sctx, X, blocks->x, INSERT_VALUES, SCATTER_FORWARD);
537:           SNESSolve(blocks->snes, NULL, blocks->x);
538:           VecScatterBegin(blocks->sctx, blocks->x, X, INSERT_VALUES, SCATTER_REVERSE);
539:           VecScatterEnd(blocks->sctx, blocks->x, X, INSERT_VALUES, SCATTER_REVERSE);
540:           blocks = blocks->next;
541:         }
542:       }
543:     } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "Unsupported or unknown composition", (int) mb->type);
544:     /* Compute F(X^{new}) */
545:     SNESComputeFunction(snes, X, F);
546:     VecNorm(F, NORM_2, &fnorm);
547:     SNESCheckFunctionNorm(snes,fnorm);

549:     if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >=0) {
550:       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
551:       break;
552:     }

554:     /* Monitor convergence */
555:     PetscObjectSAWsTakeAccess((PetscObject)snes);
556:     snes->iter = i+1;
557:     snes->norm = fnorm;
558:     PetscObjectSAWsGrantAccess((PetscObject)snes);
559:     SNESLogConvergenceHistory(snes,snes->norm,0);
560:     SNESMonitor(snes,snes->iter,snes->norm);
561:     /* Test for convergence */
562:     (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
563:     if (snes->reason) break;
564:   }
565:   if (i == maxits) {
566:     PetscInfo(snes, "Maximum number of iterations has been reached: %D\n", maxits);
567:     if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
568:   }
569:   return 0;
570: }

572: PetscErrorCode SNESMultiblockSetFields_Default(SNES snes, const char name[], PetscInt n, const PetscInt fields[])
573: {
574:   SNES_Multiblock *mb = (SNES_Multiblock*) snes->data;
575:   BlockDesc       newblock, next = mb->blocks;
576:   char            prefix[128];
577:   PetscInt        i;

579:   if (mb->defined) {
580:     PetscInfo(snes, "Ignoring new block \"%s\" because the blocks have already been defined\n", name);
581:     return 0;
582:   }
583:   for (i = 0; i < n; ++i) {
586:   }
587:   PetscNew(&newblock);
588:   if (name) {
589:     PetscStrallocpy(name, &newblock->name);
590:   } else {
591:     PetscInt len = floor(log10(mb->numBlocks))+1;

593:     PetscMalloc1(len+1, &newblock->name);
594:     PetscSNPrintf(newblock->name, len, "%s", mb->numBlocks);
595:   }
596:   newblock->nfields = n;

598:   PetscMalloc1(n, &newblock->fields);
599:   PetscArraycpy(newblock->fields, fields, n);

601:   newblock->next = NULL;

603:   SNESCreate(PetscObjectComm((PetscObject)snes), &newblock->snes);
604:   PetscObjectIncrementTabLevel((PetscObject) newblock->snes, (PetscObject) snes, 1);
605:   SNESSetType(newblock->snes, SNESNRICHARDSON);
606:   PetscLogObjectParent((PetscObject) snes, (PetscObject) newblock->snes);
607:   PetscSNPrintf(prefix, sizeof(prefix), "%smultiblock_%s_", ((PetscObject) snes)->prefix ? ((PetscObject) snes)->prefix : "", newblock->name);
608:   SNESSetOptionsPrefix(newblock->snes, prefix);

610:   if (!next) {
611:     mb->blocks         = newblock;
612:     newblock->previous = NULL;
613:   } else {
614:     while (next->next) {
615:       next = next->next;
616:     }
617:     next->next         = newblock;
618:     newblock->previous = next;
619:   }
620:   mb->numBlocks++;
621:   return 0;
622: }

624: PetscErrorCode SNESMultiblockSetIS_Default(SNES snes, const char name[], IS is)
625: {
626:   SNES_Multiblock *mb = (SNES_Multiblock*) snes->data;
627:   BlockDesc       newblock, next = mb->blocks;
628:   char            prefix[128];

630:   if (mb->defined) {
631:     PetscInfo(snes, "Ignoring new block \"%s\" because the blocks have already been defined\n", name);
632:     return 0;
633:   }
634:   PetscNew(&newblock);
635:   if (name) {
636:     PetscStrallocpy(name, &newblock->name);
637:   } else {
638:     PetscInt len = floor(log10(mb->numBlocks))+1;

640:     PetscMalloc1(len+1, &newblock->name);
641:     PetscSNPrintf(newblock->name, len, "%s", mb->numBlocks);
642:   }
643:   newblock->is = is;

645:   PetscObjectReference((PetscObject) is);

647:   newblock->next = NULL;

649:   SNESCreate(PetscObjectComm((PetscObject)snes), &newblock->snes);
650:   PetscObjectIncrementTabLevel((PetscObject) newblock->snes, (PetscObject) snes, 1);
651:   SNESSetType(newblock->snes, SNESNRICHARDSON);
652:   PetscLogObjectParent((PetscObject) snes, (PetscObject) newblock->snes);
653:   PetscSNPrintf(prefix, sizeof(prefix), "%smultiblock_%s_", ((PetscObject) snes)->prefix ? ((PetscObject) snes)->prefix : "", newblock->name);
654:   SNESSetOptionsPrefix(newblock->snes, prefix);

656:   if (!next) {
657:     mb->blocks         = newblock;
658:     newblock->previous = NULL;
659:   } else {
660:     while (next->next) {
661:       next = next->next;
662:     }
663:     next->next         = newblock;
664:     newblock->previous = next;
665:   }
666:   mb->numBlocks++;
667:   return 0;
668: }

670: PetscErrorCode  SNESMultiblockSetBlockSize_Default(SNES snes, PetscInt bs)
671: {
672:   SNES_Multiblock *mb = (SNES_Multiblock*) snes->data;

676:   mb->bs = bs;
677:   return 0;
678: }

680: PetscErrorCode SNESMultiblockGetSubSNES_Default(SNES snes, PetscInt *n, SNES **subsnes)
681: {
682:   SNES_Multiblock *mb    = (SNES_Multiblock*) snes->data;
683:   BlockDesc       blocks = mb->blocks;
684:   PetscInt        cnt    = 0;

686:   PetscMalloc1(mb->numBlocks, subsnes);
687:   while (blocks) {
688:     (*subsnes)[cnt++] = blocks->snes;
689:     blocks            = blocks->next;
690:   }

693:   if (n) *n = mb->numBlocks;
694:   return 0;
695: }

697: PetscErrorCode  SNESMultiblockSetType_Default(SNES snes, PCCompositeType type)
698: {
699:   SNES_Multiblock *mb = (SNES_Multiblock*) snes->data;

701:   mb->type = type;
702:   if (type == PC_COMPOSITE_SCHUR) {
703: #if 1
704:     SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "The Schur composite type is not yet supported");
705: #else
706:     snes->ops->solve = SNESSolve_Multiblock_Schur;
707:     snes->ops->view  = SNESView_Multiblock_Schur;

709:     PetscObjectComposeFunction((PetscObject) snes, "SNESMultiblockGetSubSNES_C", SNESMultiblockGetSubSNES_Schur);
710:     PetscObjectComposeFunction((PetscObject) snes, "SNESMultiblockSchurPrecondition_C", SNESMultiblockSchurPrecondition_Default);
711: #endif
712:   } else {
713:     snes->ops->solve = SNESSolve_Multiblock;
714:     snes->ops->view  = SNESView_Multiblock;

716:     PetscObjectComposeFunction((PetscObject) snes, "SNESMultiblockGetSubSNES_C", SNESMultiblockGetSubSNES_Default);
717:     PetscObjectComposeFunction((PetscObject) snes, "SNESMultiblockSchurPrecondition_C", 0);
718:   }
719:   return 0;
720: }

722: /*@
723:   SNESMultiblockSetFields - Sets the fields for one particular block in the solver

725:   Logically Collective on SNES

727:   Input Parameters:
728: + snes   - the solver
729: . name   - name of this block, if NULL the number of the block is used
730: . n      - the number of fields in this block
731: - fields - the fields in this block

733:   Level: intermediate

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

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

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

747: .seealso: SNESMultiblockGetSubSNES(), SNESMULTIBLOCK, SNESMultiblockSetBlockSize(), SNESMultiblockSetIS()
748: @*/
749: PetscErrorCode SNESMultiblockSetFields(SNES snes, const char name[], PetscInt n, const PetscInt *fields)
750: {
755:   PetscTryMethod(snes, "SNESMultiblockSetFields_C", (SNES, const char[], PetscInt, const PetscInt*), (snes, name, n, fields));
756:   return 0;
757: }

759: /*@
760:   SNESMultiblockSetIS - Sets the global row indices for the block

762:   Logically Collective on SNES

764:   Input Parameters:
765: + snes - the solver context
766: . name - name of this block, if NULL the number of the block is used
767: - is   - the index set that defines the global row indices in this block

769:   Notes:
770:   Use SNESMultiblockSetFields(), for blocks defined by strides.

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

775:   Level: intermediate

777: .seealso: SNESMultiblockGetSubSNES(), SNESMULTIBLOCK, SNESMultiblockSetBlockSize()
778: @*/
779: PetscErrorCode SNESMultiblockSetIS(SNES snes, const char name[], IS is)
780: {
784:   PetscTryMethod(snes, "SNESMultiblockSetIS_C", (SNES, const char[], IS), (snes, name, is));
785:   return 0;
786: }

788: /*@
789:   SNESMultiblockSetType - Sets the type of block combination.

791:   Collective on SNES

793:   Input Parameters:
794: + snes - the solver context
795: - type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE

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

800:   Level: Developer

802: .seealso: PCCompositeSetType()
803: @*/
804: PetscErrorCode SNESMultiblockSetType(SNES snes, PCCompositeType type)
805: {
807:   PetscTryMethod(snes, "SNESMultiblockSetType_C", (SNES, PCCompositeType), (snes, type));
808:   return 0;
809: }

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

814:   Logically Collective on SNES

816:   Input Parameters:
817: + snes - the solver context
818: - bs   - the block size

820:   Level: intermediate

822: .seealso: SNESMultiblockGetSubSNES(), SNESMULTIBLOCK, SNESMultiblockSetFields()
823: @*/
824: PetscErrorCode SNESMultiblockSetBlockSize(SNES snes, PetscInt bs)
825: {
828:   PetscTryMethod(snes, "SNESMultiblockSetBlockSize_C", (SNES, PetscInt), (snes,bs));
829:   return 0;
830: }

832: /*@C
833:   SNESMultiblockGetSubSNES - Gets the SNES contexts for all blocks

835:   Collective on SNES

837:   Input Parameter:
838: . snes - the solver context

840:   Output Parameters:
841: + n       - the number of blocks
842: - subsnes - the array of SNES contexts

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

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

850:   Level: advanced

852: .seealso: SNESMULTIBLOCK
853: @*/
854: PetscErrorCode SNESMultiblockGetSubSNES(SNES snes, PetscInt *n, SNES *subsnes[])
855: {
858:   PetscUseMethod(snes, "SNESMultiblockGetSubSNES_C", (SNES, PetscInt*, SNES **), (snes, n, subsnes));
859:   return 0;
860: }

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

866:   Level: beginner

868: .seealso:  SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESNEWTONTR, SNESNRICHARDSON
869: M*/
870: PETSC_EXTERN PetscErrorCode SNESCreate_Multiblock(SNES snes)
871: {
872:   SNES_Multiblock *mb;

874:   snes->ops->destroy        = SNESDestroy_Multiblock;
875:   snes->ops->setup          = SNESSetUp_Multiblock;
876:   snes->ops->setfromoptions = SNESSetFromOptions_Multiblock;
877:   snes->ops->view           = SNESView_Multiblock;
878:   snes->ops->solve          = SNESSolve_Multiblock;
879:   snes->ops->reset          = SNESReset_Multiblock;

881:   snes->usesksp = PETSC_FALSE;
882:   snes->usesnpc = PETSC_FALSE;

884:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

886:   PetscNewLog(snes,&mb);
887:   snes->data    = (void*) mb;
888:   mb->defined   = PETSC_FALSE;
889:   mb->numBlocks = 0;
890:   mb->bs        = -1;
891:   mb->type      = PC_COMPOSITE_MULTIPLICATIVE;

893:   /* We attach functions so that they can be called on another PC without crashing the program */
894:   PetscObjectComposeFunction((PetscObject) snes, "SNESMultiblockSetFields_C",    SNESMultiblockSetFields_Default);
895:   PetscObjectComposeFunction((PetscObject) snes, "SNESMultiblockSetIS_C",        SNESMultiblockSetIS_Default);
896:   PetscObjectComposeFunction((PetscObject) snes, "SNESMultiblockSetType_C",      SNESMultiblockSetType_Default);
897:   PetscObjectComposeFunction((PetscObject) snes, "SNESMultiblockSetBlockSize_C", SNESMultiblockSetBlockSize_Default);
898:   PetscObjectComposeFunction((PetscObject) snes, "SNESMultiblockGetSubSNES_C",   SNESMultiblockGetSubSNES_Default);
899:   return 0;
900: }