Actual source code: multiblock.c

petsc-3.15.0 2021-04-05
Report Typos and Errors
  1: #include <petsc/private/snesimpl.h>
  2: #include <petscdmcomposite.h>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

501:   snes->reason = SNES_CONVERGED_ITERATING;

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

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

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

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

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

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

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

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

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

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

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

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

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

619:   newblock->next = NULL;

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

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

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

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

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

665:   PetscObjectReference((PetscObject) is);

667:   newblock->next = NULL;

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

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

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

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

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

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

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

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

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

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

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

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

750:   Logically Collective on SNES

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

758:   Level: intermediate

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

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

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

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

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

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

790:   Logically Collective on SNES

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

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

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

803:   Level: intermediate

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

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

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

822:   Collective on SNES

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

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

831:   Level: Developer

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

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

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

848:   Logically Collective on SNES

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

854:   Level: intermediate

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

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

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

872:   Collective on SNES

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

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

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

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

887:   Level: advanced

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

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

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

906:   Level: beginner

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

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

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

926:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

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

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