Actual source code: multiblock.c
petsc-3.14.6 2021-03-30
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: }