Actual source code: multiblock.c
1: #include <petsc/private/snesimpl.h>
2: #include <petscdmcomposite.h>
4: typedef struct _BlockDesc *BlockDesc;
5: struct _BlockDesc {
6: char *name; /* Block name */
7: PetscInt nfields; /* If block is defined on a DA, the number of DA fields */
8: PetscInt *fields; /* If block is defined on a DA, the list of DA fields */
9: IS is; /* Index sets defining the block */
10: VecScatter sctx; /* Scatter mapping global Vec to blockVec */
11: SNES snes; /* Solver for this block */
12: Vec x;
13: BlockDesc next, previous;
14: };
16: typedef struct {
17: PetscBool issetup; /* Flag is true after the all ISs and operators have been defined */
18: PetscBool defined; /* Flag is true after the blocks have been defined, to prevent more blocks from being added */
19: PetscBool defaultblocks; /* Flag is true for a system with a set of 'k' scalar fields with the same layout (and bs = k) */
20: PetscInt numBlocks; /* Number of blocks (can be fields, domains, etc.) */
21: PetscInt bs; /* Block size for IS, Vec and Mat structures */
22: PCCompositeType type; /* Solver combination method (additive, multiplicative, etc.) */
23: BlockDesc blocks; /* Linked list of block descriptors */
24: } SNES_Multiblock;
26: PetscErrorCode SNESReset_Multiblock(SNES snes)
27: {
28: SNES_Multiblock *mb = (SNES_Multiblock *)snes->data;
29: BlockDesc blocks = mb->blocks, next;
31: while (blocks) {
32: SNESReset(blocks->snes);
33: #if 0
34: VecDestroy(&blocks->x);
35: #endif
36: VecScatterDestroy(&blocks->sctx);
37: ISDestroy(&blocks->is);
38: next = blocks->next;
39: blocks = next;
40: }
41: return 0;
42: }
44: /*
45: SNESDestroy_Multiblock - Destroys the private SNES_Multiblock context that was created with SNESCreate_Multiblock().
47: Input Parameter:
48: . snes - the SNES context
50: Application Interface Routine: SNESDestroy()
51: */
52: PetscErrorCode SNESDestroy_Multiblock(SNES snes)
53: {
54: SNES_Multiblock *mb = (SNES_Multiblock *)snes->data;
55: BlockDesc blocks = mb->blocks, next;
57: SNESReset_Multiblock(snes);
58: while (blocks) {
59: next = blocks->next;
60: SNESDestroy(&blocks->snes);
61: PetscFree(blocks->name);
62: PetscFree(blocks->fields);
63: PetscFree(blocks);
64: blocks = next;
65: }
66: PetscFree(snes->data);
67: return 0;
68: }
70: /* Precondition: blocksize is set to a meaningful value */
71: static PetscErrorCode SNESMultiblockSetFieldsRuntime_Private(SNES snes)
72: {
73: SNES_Multiblock *mb = (SNES_Multiblock *)snes->data;
74: PetscInt *ifields;
75: PetscInt i, nfields;
76: PetscBool flg = PETSC_TRUE;
77: char optionname[128], name[8];
79: PetscMalloc1(mb->bs, &ifields);
80: for (i = 0;; ++i) {
81: PetscSNPrintf(name, sizeof(name), "%" PetscInt_FMT, i);
82: PetscSNPrintf(optionname, sizeof(optionname), "-snes_multiblock_%" PetscInt_FMT "_fields", i);
83: nfields = mb->bs;
84: PetscOptionsGetIntArray(NULL, ((PetscObject)snes)->prefix, optionname, ifields, &nfields, &flg);
85: if (!flg) break;
87: SNESMultiblockSetFields(snes, name, nfields, ifields);
88: }
89: if (i > 0) {
90: /* Makes command-line setting of blocks take precedence over setting them in code.
91: Otherwise subsequent calls to SNESMultiblockSetIS() or SNESMultiblockSetFields() would
92: create new blocks, which would probably not be what the user wanted. */
93: mb->defined = PETSC_TRUE;
94: }
95: PetscFree(ifields);
96: return 0;
97: }
99: static PetscErrorCode SNESMultiblockSetDefaults(SNES snes)
100: {
101: SNES_Multiblock *mb = (SNES_Multiblock *)snes->data;
102: BlockDesc blocks = mb->blocks;
103: PetscInt i;
105: if (!blocks) {
106: if (snes->dm) {
107: PetscBool dmcomposite;
109: PetscObjectTypeCompare((PetscObject)snes->dm, DMCOMPOSITE, &dmcomposite);
110: if (dmcomposite) {
111: PetscInt nDM;
112: IS *fields;
114: PetscInfo(snes, "Setting up physics based multiblock solver using the embedded DM\n");
115: DMCompositeGetNumberDM(snes->dm, &nDM);
116: DMCompositeGetGlobalISs(snes->dm, &fields);
117: for (i = 0; i < nDM; ++i) {
118: char name[8];
120: PetscSNPrintf(name, sizeof(name), "%" PetscInt_FMT, i);
121: SNESMultiblockSetIS(snes, name, fields[i]);
122: ISDestroy(&fields[i]);
123: }
124: PetscFree(fields);
125: }
126: } else {
127: PetscBool flg = PETSC_FALSE;
128: PetscBool stokes = PETSC_FALSE;
130: if (mb->bs <= 0) {
131: if (snes->jacobian_pre) {
132: MatGetBlockSize(snes->jacobian_pre, &mb->bs);
133: } else mb->bs = 1;
134: }
136: PetscOptionsGetBool(NULL, ((PetscObject)snes)->prefix, "-snes_multiblock_default", &flg, NULL);
137: PetscOptionsGetBool(NULL, ((PetscObject)snes)->prefix, "-snes_multiblock_detect_saddle_point", &stokes, NULL);
138: if (stokes) {
139: IS zerodiags, rest;
140: PetscInt nmin, nmax;
142: MatGetOwnershipRange(snes->jacobian_pre, &nmin, &nmax);
143: MatFindZeroDiagonals(snes->jacobian_pre, &zerodiags);
144: ISComplement(zerodiags, nmin, nmax, &rest);
145: SNESMultiblockSetIS(snes, "0", rest);
146: SNESMultiblockSetIS(snes, "1", zerodiags);
147: ISDestroy(&zerodiags);
148: ISDestroy(&rest);
149: } else {
150: if (!flg) {
151: /* Allow user to set fields from command line, if bs was known at the time of SNESSetFromOptions_Multiblock()
152: then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */
153: SNESMultiblockSetFieldsRuntime_Private(snes);
154: if (mb->defined) PetscInfo(snes, "Blocks defined using the options database\n");
155: }
156: if (flg || !mb->defined) {
157: PetscInfo(snes, "Using default splitting of fields\n");
158: for (i = 0; i < mb->bs; ++i) {
159: char name[8];
161: PetscSNPrintf(name, sizeof(name), "%" PetscInt_FMT, i);
162: SNESMultiblockSetFields(snes, name, 1, &i);
163: }
164: mb->defaultblocks = PETSC_TRUE;
165: }
166: }
167: }
168: } else if (mb->numBlocks == 1) {
169: if (blocks->is) {
170: IS is2;
171: PetscInt nmin, nmax;
173: MatGetOwnershipRange(snes->jacobian_pre, &nmin, &nmax);
174: ISComplement(blocks->is, nmin, nmax, &is2);
175: SNESMultiblockSetIS(snes, "1", is2);
176: ISDestroy(&is2);
177: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least two sets of fields to SNES multiblock");
178: }
180: return 0;
181: }
183: /*
184: SNESSetUp_Multiblock - Sets up the internal data structures for the later use
185: of the SNESMULTIBLOCK nonlinear solver.
187: Input Parameters:
188: + snes - the SNES context
189: - x - the solution vector
191: Application Interface Routine: SNESSetUp()
192: */
193: PetscErrorCode SNESSetUp_Multiblock(SNES snes)
194: {
195: SNES_Multiblock *mb = (SNES_Multiblock *)snes->data;
196: BlockDesc blocks;
197: PetscInt i, numBlocks;
199: SNESMultiblockSetDefaults(snes);
200: numBlocks = mb->numBlocks;
201: blocks = mb->blocks;
203: /* Create ISs */
204: if (!mb->issetup) {
205: PetscInt ccsize, rstart, rend, nslots, bs;
206: PetscBool sorted;
208: mb->issetup = PETSC_TRUE;
209: bs = mb->bs;
210: MatGetOwnershipRange(snes->jacobian_pre, &rstart, &rend);
211: MatGetLocalSize(snes->jacobian_pre, NULL, &ccsize);
212: nslots = (rend - rstart) / bs;
213: for (i = 0; i < numBlocks; ++i) {
214: if (mb->defaultblocks) {
215: ISCreateStride(PetscObjectComm((PetscObject)snes), nslots, rstart + i, numBlocks, &blocks->is);
216: } else if (!blocks->is) {
217: if (blocks->nfields > 1) {
218: PetscInt *ii, j, k, nfields = blocks->nfields, *fields = blocks->fields;
220: PetscMalloc1(nfields * nslots, &ii);
221: for (j = 0; j < nslots; ++j) {
222: for (k = 0; k < nfields; ++k) ii[nfields * j + k] = rstart + bs * j + fields[k];
223: }
224: ISCreateGeneral(PetscObjectComm((PetscObject)snes), nslots * nfields, ii, PETSC_OWN_POINTER, &blocks->is);
225: } else {
226: ISCreateStride(PetscObjectComm((PetscObject)snes), nslots, rstart + blocks->fields[0], bs, &blocks->is);
227: }
228: }
229: ISSorted(blocks->is, &sorted);
231: blocks = blocks->next;
232: }
233: }
235: #if 0
236: /* Create matrices */
237: ilink = jac->head;
238: if (!jac->pmat) {
239: PetscMalloc1(nsplit,&jac->pmat);
240: for (i=0; i<nsplit; i++) {
241: MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is,MAT_INITIAL_MATRIX,&jac->pmat[i]);
242: ilink = ilink->next;
243: }
244: } else {
245: for (i=0; i<nsplit; i++) {
246: MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is,MAT_REUSE_MATRIX,&jac->pmat[i]);
247: ilink = ilink->next;
248: }
249: }
250: if (jac->realdiagonal) {
251: ilink = jac->head;
252: if (!jac->mat) {
253: PetscMalloc1(nsplit,&jac->mat);
254: for (i=0; i<nsplit; i++) {
255: MatCreateSubMatrix(pc->mat,ilink->is,ilink->is,MAT_INITIAL_MATRIX,&jac->mat[i]);
256: ilink = ilink->next;
257: }
258: } else {
259: for (i=0; i<nsplit; i++) {
260: if (jac->mat[i]) MatCreateSubMatrix(pc->mat,ilink->is,ilink->is,MAT_REUSE_MATRIX,&jac->mat[i]);
261: ilink = ilink->next;
262: }
263: }
264: } else jac->mat = jac->pmat;
265: #endif
267: #if 0
268: if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR) {
269: /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */
270: ilink = jac->head;
271: if (!jac->Afield) {
272: PetscMalloc1(nsplit,&jac->Afield);
273: for (i=0; i<nsplit; i++) {
274: MatCreateSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);
275: ilink = ilink->next;
276: }
277: } else {
278: for (i=0; i<nsplit; i++) {
279: MatCreateSubMatrix(pc->mat,ilink->is,NULL,MAT_REUSE_MATRIX,&jac->Afield[i]);
280: ilink = ilink->next;
281: }
282: }
283: }
284: #endif
286: if (mb->type == PC_COMPOSITE_SCHUR) {
287: #if 0
288: IS ccis;
289: PetscInt rstart,rend;
292: /* When extracting off-diagonal submatrices, we take complements from this range */
293: MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);
295: /* need to handle case when one is resetting up the preconditioner */
296: if (jac->schur) {
297: ilink = jac->head;
298: ISComplement(ilink->is,rstart,rend,&ccis);
299: MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);
300: ISDestroy(&ccis);
301: ilink = ilink->next;
302: ISComplement(ilink->is,rstart,rend,&ccis);
303: MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);
304: ISDestroy(&ccis);
305: MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->pmat[1]);
306: KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),pc->flag);
308: } else {
309: KSP ksp;
310: char schurprefix[256];
312: /* extract the A01 and A10 matrices */
313: ilink = jac->head;
314: ISComplement(ilink->is,rstart,rend,&ccis);
315: MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);
316: ISDestroy(&ccis);
317: ilink = ilink->next;
318: ISComplement(ilink->is,rstart,rend,&ccis);
319: MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);
320: ISDestroy(&ccis);
321: /* Use mat[0] (diagonal block of the real matrix) preconditioned by pmat[0] */
322: MatCreateSchurComplement(jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1],&jac->schur);
323: /* set tabbing and options prefix of KSP inside the MatSchur */
324: MatSchurComplementGetKSP(jac->schur,&ksp);
325: PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,2);
326: PetscSNPrintf(schurprefix,sizeof(schurprefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",jac->head->splitname);
327: KSPSetOptionsPrefix(ksp,schurprefix);
328: MatSetFromOptions(jac->schur);
330: KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);
331: PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);
332: KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));
333: if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) {
334: PC pc;
335: KSPGetPC(jac->kspschur,&pc);
336: PCSetType(pc,PCNONE);
337: /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */
338: }
339: PetscSNPrintf(schurprefix,sizeof(schurprefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);
340: KSPSetOptionsPrefix(jac->kspschur,schurprefix);
341: /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
342: KSPSetFromOptions(jac->kspschur);
344: PetscMalloc2(2,&jac->x,2,&jac->y);
345: MatCreateVecs(jac->pmat[0],&jac->x[0],&jac->y[0]);
346: MatCreateVecs(jac->pmat[1],&jac->x[1],&jac->y[1]);
347: ilink = jac->head;
348: ilink->x = jac->x[0]; ilink->y = jac->y[0];
349: ilink = ilink->next;
350: ilink->x = jac->x[1]; ilink->y = jac->y[1];
351: }
352: #endif
353: } else {
354: /* Set up the individual SNESs */
355: blocks = mb->blocks;
356: i = 0;
357: while (blocks) {
358: /*TODO: Set these correctly */
359: /* SNESSetFunction(blocks->snes, blocks->x, func); */
360: /* SNESSetJacobian(blocks->snes, blocks->x, jac); */
361: VecDuplicate(blocks->snes->vec_sol, &blocks->x);
362: /* really want setfromoptions called in SNESSetFromOptions_Multiblock(), but it is not ready yet */
363: SNESSetFromOptions(blocks->snes);
364: SNESSetUp(blocks->snes);
365: blocks = blocks->next;
366: i++;
367: }
368: }
370: /* Compute scatter contexts needed by multiplicative versions and non-default splits */
371: if (!mb->blocks->sctx) {
372: Vec xtmp;
374: blocks = mb->blocks;
375: MatCreateVecs(snes->jacobian_pre, &xtmp, NULL);
376: while (blocks) {
377: VecScatterCreate(xtmp, blocks->is, blocks->x, NULL, &blocks->sctx);
378: blocks = blocks->next;
379: }
380: VecDestroy(&xtmp);
381: }
382: return 0;
383: }
385: /*
386: SNESSetFromOptions_Multiblock - Sets various parameters for the SNESMULTIBLOCK method.
388: Input Parameter:
389: . snes - the SNES context
391: Application Interface Routine: SNESSetFromOptions()
392: */
393: static PetscErrorCode SNESSetFromOptions_Multiblock(SNES snes, PetscOptionItems *PetscOptionsObject)
394: {
395: SNES_Multiblock *mb = (SNES_Multiblock *)snes->data;
396: PCCompositeType ctype;
397: PetscInt bs;
398: PetscBool flg;
400: PetscOptionsHeadBegin(PetscOptionsObject, "SNES Multiblock options");
401: PetscOptionsInt("-snes_multiblock_block_size", "Blocksize that defines number of fields", "PCFieldSplitSetBlockSize", mb->bs, &bs, &flg);
402: if (flg) SNESMultiblockSetBlockSize(snes, bs);
403: PetscOptionsEnum("-snes_multiblock_type", "Type of composition", "PCFieldSplitSetType", PCCompositeTypes, (PetscEnum)mb->type, (PetscEnum *)&ctype, &flg);
404: if (flg) SNESMultiblockSetType(snes, ctype);
405: /* Only setup fields once */
406: if ((mb->bs > 0) && (mb->numBlocks == 0)) {
407: /* only allow user to set fields from command line if bs is already known, otherwise user can set them in SNESMultiblockSetDefaults() */
408: SNESMultiblockSetFieldsRuntime_Private(snes);
409: if (mb->defined) PetscInfo(snes, "Blocks defined using the options database\n");
410: }
411: PetscOptionsHeadEnd();
412: return 0;
413: }
415: /*
416: SNESView_Multiblock - Prints info from the SNESMULTIBLOCK data structure.
418: Input Parameters:
419: + SNES - the SNES context
420: - viewer - visualization context
422: Application Interface Routine: SNESView()
423: */
424: static PetscErrorCode SNESView_Multiblock(SNES snes, PetscViewer viewer)
425: {
426: SNES_Multiblock *mb = (SNES_Multiblock *)snes->data;
427: BlockDesc blocks = mb->blocks;
428: PetscBool iascii;
430: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
431: if (iascii) {
432: PetscViewerASCIIPrintf(viewer, " Multiblock with %s composition: total blocks = %" PetscInt_FMT ", blocksize = %" PetscInt_FMT "\n", PCCompositeTypes[mb->type], mb->numBlocks, mb->bs);
433: PetscViewerASCIIPrintf(viewer, " Solver info for each split is in the following SNES objects:\n");
434: PetscViewerASCIIPushTab(viewer);
435: while (blocks) {
436: if (blocks->fields) {
437: PetscInt j;
439: PetscViewerASCIIPrintf(viewer, " Block %s Fields ", blocks->name);
440: PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
441: for (j = 0; j < blocks->nfields; ++j) {
442: if (j > 0) PetscViewerASCIIPrintf(viewer, ",");
443: PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, blocks->fields[j]);
444: }
445: PetscViewerASCIIPrintf(viewer, "\n");
446: PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
447: } else {
448: PetscViewerASCIIPrintf(viewer, " Block %s Defined by IS\n", blocks->name);
449: }
450: SNESView(blocks->snes, viewer);
451: blocks = blocks->next;
452: }
453: PetscViewerASCIIPopTab(viewer);
454: }
455: return 0;
456: }
458: /*
459: SNESSolve_Multiblock - Solves a nonlinear system with the Multiblock method.
461: Input Parameters:
462: . snes - the SNES context
464: Output Parameter:
465: . outits - number of iterations until termination
467: Application Interface Routine: SNESSolve()
468: */
469: PetscErrorCode SNESSolve_Multiblock(SNES snes)
470: {
471: SNES_Multiblock *mb = (SNES_Multiblock *)snes->data;
472: Vec X, Y, F;
473: PetscReal fnorm;
474: PetscInt maxits, i;
478: snes->reason = SNES_CONVERGED_ITERATING;
480: maxits = snes->max_its; /* maximum number of iterations */
481: X = snes->vec_sol; /* X^n */
482: Y = snes->vec_sol_update; /* \tilde X */
483: F = snes->vec_func; /* residual vector */
485: VecSetBlockSize(X, mb->bs);
486: VecSetBlockSize(Y, mb->bs);
487: VecSetBlockSize(F, mb->bs);
488: PetscObjectSAWsTakeAccess((PetscObject)snes);
489: snes->iter = 0;
490: snes->norm = 0.;
491: PetscObjectSAWsGrantAccess((PetscObject)snes);
493: if (!snes->vec_func_init_set) {
494: SNESComputeFunction(snes, X, F);
495: } else snes->vec_func_init_set = PETSC_FALSE;
497: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
498: SNESCheckFunctionNorm(snes, fnorm);
499: PetscObjectSAWsTakeAccess((PetscObject)snes);
500: snes->norm = fnorm;
501: PetscObjectSAWsGrantAccess((PetscObject)snes);
502: SNESLogConvergenceHistory(snes, fnorm, 0);
503: SNESMonitor(snes, 0, fnorm);
505: /* test convergence */
506: PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
507: if (snes->reason) return 0;
509: for (i = 0; i < maxits; i++) {
510: /* Call general purpose update function */
511: PetscTryTypeMethod(snes, update, snes->iter);
512: /* Compute X^{new} from subsolves */
513: if (mb->type == PC_COMPOSITE_ADDITIVE) {
514: BlockDesc blocks = mb->blocks;
516: if (mb->defaultblocks) {
517: /*TODO: Make an array of Vecs for this */
518: /* VecStrideGatherAll(X, mb->x, INSERT_VALUES); */
519: while (blocks) {
520: SNESSolve(blocks->snes, NULL, blocks->x);
521: blocks = blocks->next;
522: }
523: /* VecStrideScatterAll(mb->x, X, INSERT_VALUES); */
524: } else {
525: while (blocks) {
526: VecScatterBegin(blocks->sctx, X, blocks->x, INSERT_VALUES, SCATTER_FORWARD);
527: VecScatterEnd(blocks->sctx, X, blocks->x, INSERT_VALUES, SCATTER_FORWARD);
528: SNESSolve(blocks->snes, NULL, blocks->x);
529: VecScatterBegin(blocks->sctx, blocks->x, X, INSERT_VALUES, SCATTER_REVERSE);
530: VecScatterEnd(blocks->sctx, blocks->x, X, INSERT_VALUES, SCATTER_REVERSE);
531: blocks = blocks->next;
532: }
533: }
534: } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "Unsupported or unknown composition %d", (int)mb->type);
535: /* Compute F(X^{new}) */
536: SNESComputeFunction(snes, X, F);
537: VecNorm(F, NORM_2, &fnorm);
538: SNESCheckFunctionNorm(snes, fnorm);
540: if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
541: snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
542: break;
543: }
545: /* Monitor convergence */
546: PetscObjectSAWsTakeAccess((PetscObject)snes);
547: snes->iter = i + 1;
548: snes->norm = fnorm;
549: PetscObjectSAWsGrantAccess((PetscObject)snes);
550: SNESLogConvergenceHistory(snes, snes->norm, 0);
551: SNESMonitor(snes, snes->iter, snes->norm);
552: /* Test for convergence */
553: PetscUseTypeMethod(snes, converged, snes->iter, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
554: if (snes->reason) break;
555: }
556: if (i == maxits) {
557: PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", maxits);
558: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
559: }
560: return 0;
561: }
563: PetscErrorCode SNESMultiblockSetFields_Default(SNES snes, const char name[], PetscInt n, const PetscInt fields[])
564: {
565: SNES_Multiblock *mb = (SNES_Multiblock *)snes->data;
566: BlockDesc newblock, next = mb->blocks;
567: char prefix[128];
568: PetscInt i;
570: if (mb->defined) {
571: PetscInfo(snes, "Ignoring new block \"%s\" because the blocks have already been defined\n", name);
572: return 0;
573: }
574: for (i = 0; i < n; ++i) {
577: }
578: PetscNew(&newblock);
579: if (name) {
580: PetscStrallocpy(name, &newblock->name);
581: } else {
582: PetscInt len = floor(log10(mb->numBlocks)) + 1;
584: PetscMalloc1(len + 1, &newblock->name);
585: PetscSNPrintf(newblock->name, len, "%" PetscInt_FMT, mb->numBlocks);
586: }
587: newblock->nfields = n;
589: PetscMalloc1(n, &newblock->fields);
590: PetscArraycpy(newblock->fields, fields, n);
592: newblock->next = NULL;
594: SNESCreate(PetscObjectComm((PetscObject)snes), &newblock->snes);
595: PetscObjectIncrementTabLevel((PetscObject)newblock->snes, (PetscObject)snes, 1);
596: SNESSetType(newblock->snes, SNESNRICHARDSON);
597: PetscSNPrintf(prefix, sizeof(prefix), "%smultiblock_%s_", ((PetscObject)snes)->prefix ? ((PetscObject)snes)->prefix : "", newblock->name);
598: SNESSetOptionsPrefix(newblock->snes, prefix);
600: if (!next) {
601: mb->blocks = newblock;
602: newblock->previous = NULL;
603: } else {
604: while (next->next) next = next->next;
605: next->next = newblock;
606: newblock->previous = next;
607: }
608: mb->numBlocks++;
609: return 0;
610: }
612: PetscErrorCode SNESMultiblockSetIS_Default(SNES snes, const char name[], IS is)
613: {
614: SNES_Multiblock *mb = (SNES_Multiblock *)snes->data;
615: BlockDesc newblock, next = mb->blocks;
616: char prefix[128];
618: if (mb->defined) {
619: PetscInfo(snes, "Ignoring new block \"%s\" because the blocks have already been defined\n", name);
620: return 0;
621: }
622: PetscNew(&newblock);
623: if (name) {
624: PetscStrallocpy(name, &newblock->name);
625: } else {
626: PetscInt len = floor(log10(mb->numBlocks)) + 1;
628: PetscMalloc1(len + 1, &newblock->name);
629: PetscSNPrintf(newblock->name, len, "%" PetscInt_FMT, mb->numBlocks);
630: }
631: newblock->is = is;
633: PetscObjectReference((PetscObject)is);
635: newblock->next = NULL;
637: SNESCreate(PetscObjectComm((PetscObject)snes), &newblock->snes);
638: PetscObjectIncrementTabLevel((PetscObject)newblock->snes, (PetscObject)snes, 1);
639: SNESSetType(newblock->snes, SNESNRICHARDSON);
640: PetscSNPrintf(prefix, sizeof(prefix), "%smultiblock_%s_", ((PetscObject)snes)->prefix ? ((PetscObject)snes)->prefix : "", newblock->name);
641: SNESSetOptionsPrefix(newblock->snes, prefix);
643: if (!next) {
644: mb->blocks = newblock;
645: newblock->previous = NULL;
646: } else {
647: while (next->next) next = next->next;
648: next->next = newblock;
649: newblock->previous = next;
650: }
651: mb->numBlocks++;
652: return 0;
653: }
655: PetscErrorCode SNESMultiblockSetBlockSize_Default(SNES snes, PetscInt bs)
656: {
657: SNES_Multiblock *mb = (SNES_Multiblock *)snes->data;
661: mb->bs = bs;
662: return 0;
663: }
665: PetscErrorCode SNESMultiblockGetSubSNES_Default(SNES snes, PetscInt *n, SNES **subsnes)
666: {
667: SNES_Multiblock *mb = (SNES_Multiblock *)snes->data;
668: BlockDesc blocks = mb->blocks;
669: PetscInt cnt = 0;
671: PetscMalloc1(mb->numBlocks, subsnes);
672: while (blocks) {
673: (*subsnes)[cnt++] = blocks->snes;
674: blocks = blocks->next;
675: }
678: if (n) *n = mb->numBlocks;
679: return 0;
680: }
682: PetscErrorCode SNESMultiblockSetType_Default(SNES snes, PCCompositeType type)
683: {
684: SNES_Multiblock *mb = (SNES_Multiblock *)snes->data;
686: mb->type = type;
687: if (type == PC_COMPOSITE_SCHUR) {
688: #if 1
689: SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "The Schur composite type is not yet supported");
690: #else
691: snes->ops->solve = SNESSolve_Multiblock_Schur;
692: snes->ops->view = SNESView_Multiblock_Schur;
694: PetscObjectComposeFunction((PetscObject)snes, "SNESMultiblockGetSubSNES_C", SNESMultiblockGetSubSNES_Schur);
695: PetscObjectComposeFunction((PetscObject)snes, "SNESMultiblockSchurPrecondition_C", SNESMultiblockSchurPrecondition_Default);
696: #endif
697: } else {
698: snes->ops->solve = SNESSolve_Multiblock;
699: snes->ops->view = SNESView_Multiblock;
701: PetscObjectComposeFunction((PetscObject)snes, "SNESMultiblockGetSubSNES_C", SNESMultiblockGetSubSNES_Default);
702: PetscObjectComposeFunction((PetscObject)snes, "SNESMultiblockSchurPrecondition_C", 0);
703: }
704: return 0;
705: }
707: /*@
708: SNESMultiblockSetFields - Sets the fields for one particular block in a `SNESMULTBLOCK` solver
710: Logically Collective
712: Input Parameters:
713: + snes - the solver
714: . name - name of this block, if NULL the number of the block is used
715: . n - the number of fields in this block
716: - fields - the fields in this block
718: Level: intermediate
720: Notes:
721: Use `SNESMultiblockSetIS()` to set a completely general set of row indices as a block.
723: The `SNESMultiblockSetFields()` is for defining blocks as a group of strided indices, or fields.
724: For example, if the vector block size is three then one can define a block as field 0, or
725: 1 or 2, or field 0,1 or 0,2 or 1,2 which means
726: 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x....
727: where the numbered entries indicate what is in the block.
729: This function is called once per block (it creates a new block each time). Solve options
730: for this block will be available under the prefix -multiblock_BLOCKNAME_.
732: .seealso: `SNESMULTBLOCK`, `SNESMultiblockGetSubSNES()`, `SNESMULTIBLOCK`, `SNESMultiblockSetBlockSize()`, `SNESMultiblockSetIS()`
733: @*/
734: PetscErrorCode SNESMultiblockSetFields(SNES snes, const char name[], PetscInt n, const PetscInt *fields)
735: {
740: PetscTryMethod(snes, "SNESMultiblockSetFields_C", (SNES, const char[], PetscInt, const PetscInt *), (snes, name, n, fields));
741: return 0;
742: }
744: /*@
745: SNESMultiblockSetIS - Sets the global row indices for one particular block in a `SNESMULTBLOCK` solver
747: Logically Collective
749: Input Parameters:
750: + snes - the solver context
751: . name - name of this block, if NULL the number of the block is used
752: - is - the index set that defines the global row indices in this block
754: Notes:
755: Use `SNESMultiblockSetFields()`, for blocks defined by strides.
757: This function is called once per block (it creates a new block each time). Solve options
758: for this block will be available under the prefix -multiblock_BLOCKNAME_.
760: Level: intermediate
762: .seealso: `SNESMULTBLOCK`, `SNESMultiblockGetSubSNES()`, `SNESMULTIBLOCK`, `SNESMultiblockSetBlockSize()`
763: @*/
764: PetscErrorCode SNESMultiblockSetIS(SNES snes, const char name[], IS is)
765: {
769: PetscTryMethod(snes, "SNESMultiblockSetIS_C", (SNES, const char[], IS), (snes, name, is));
770: return 0;
771: }
773: /*@
774: SNESMultiblockSetType - Sets the type of block combination used for a `SNESMULTBLOCK` solver
776: Logically Collective
778: Input Parameters:
779: + snes - the solver context
780: - type - `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE` (default), `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`
782: Options Database Key:
783: . -snes_multiblock_type <type: one of multiplicative, additive, symmetric_multiplicative> - Sets block combination type
785: Level: advanced
787: .seealso: `SNESMULTBLOCK`, `PCCompositeSetType()`, `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`, `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`
788: @*/
789: PetscErrorCode SNESMultiblockSetType(SNES snes, PCCompositeType type)
790: {
792: PetscTryMethod(snes, "SNESMultiblockSetType_C", (SNES, PCCompositeType), (snes, type));
793: return 0;
794: }
796: /*@
797: SNESMultiblockSetBlockSize - Sets the block size for structured block division in a `SNESMULTBLOCK` solver. If not set the matrix block size is used.
799: Logically Collective
801: Input Parameters:
802: + snes - the solver context
803: - bs - the block size
805: Level: intermediate
807: .seealso: `SNESMULTBLOCK`, `SNESMultiblockGetSubSNES()`, `SNESMULTIBLOCK`, `SNESMultiblockSetFields()`
808: @*/
809: PetscErrorCode SNESMultiblockSetBlockSize(SNES snes, PetscInt bs)
810: {
813: PetscTryMethod(snes, "SNESMultiblockSetBlockSize_C", (SNES, PetscInt), (snes, bs));
814: return 0;
815: }
817: /*@C
818: SNESMultiblockGetSubSNES - Gets the `SNES` contexts for all blocks in a `SNESMULTBLOCK` solver.
820: Not Collective but each `SNES` obtained is parallel
822: Input Parameter:
823: . snes - the solver context
825: Output Parameters:
826: + n - the number of blocks
827: - subsnes - the array of `SNES` contexts
829: Note:
830: After `SNESMultiblockGetSubSNES()` the array of `SNES`s MUST be freed by the user
831: (not each `SNES`, just the array that contains them).
833: You must call `SNESSetUp()` before calling `SNESMultiblockGetSubSNES()`.
835: Level: advanced
837: .seealso: `SNESMULTBLOCK`, `SNESMultiblockSetIS()`, `SNESMultiblockSetFields()`
838: @*/
839: PetscErrorCode SNESMultiblockGetSubSNES(SNES snes, PetscInt *n, SNES *subsnes[])
840: {
843: PetscUseMethod(snes, "SNESMultiblockGetSubSNES_C", (SNES, PetscInt *, SNES **), (snes, n, subsnes));
844: return 0;
845: }
847: /*MC
848: SNESMULTIBLOCK - Multiblock nonlinear solver that can use overlapping or nonoverlapping blocks, organized
849: additively (Jacobi) or multiplicatively (Gauss-Seidel).
851: Level: beginner
853: .seealso: `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESNEWTONTR`, `SNESNRICHARDSON`, `SNESMultiblockSetType()`,
854: `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`, `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`
855: M*/
856: PETSC_EXTERN PetscErrorCode SNESCreate_Multiblock(SNES snes)
857: {
858: SNES_Multiblock *mb;
860: snes->ops->destroy = SNESDestroy_Multiblock;
861: snes->ops->setup = SNESSetUp_Multiblock;
862: snes->ops->setfromoptions = SNESSetFromOptions_Multiblock;
863: snes->ops->view = SNESView_Multiblock;
864: snes->ops->solve = SNESSolve_Multiblock;
865: snes->ops->reset = SNESReset_Multiblock;
867: snes->usesksp = PETSC_FALSE;
868: snes->usesnpc = PETSC_FALSE;
870: snes->alwayscomputesfinalresidual = PETSC_TRUE;
872: PetscNew(&mb);
873: snes->data = (void *)mb;
874: mb->defined = PETSC_FALSE;
875: mb->numBlocks = 0;
876: mb->bs = -1;
877: mb->type = PC_COMPOSITE_MULTIPLICATIVE;
879: /* We attach functions so that they can be called on another PC without crashing the program */
880: PetscObjectComposeFunction((PetscObject)snes, "SNESMultiblockSetFields_C", SNESMultiblockSetFields_Default);
881: PetscObjectComposeFunction((PetscObject)snes, "SNESMultiblockSetIS_C", SNESMultiblockSetIS_Default);
882: PetscObjectComposeFunction((PetscObject)snes, "SNESMultiblockSetType_C", SNESMultiblockSetType_Default);
883: PetscObjectComposeFunction((PetscObject)snes, "SNESMultiblockSetBlockSize_C", SNESMultiblockSetBlockSize_Default);
884: PetscObjectComposeFunction((PetscObject)snes, "SNESMultiblockGetSubSNES_C", SNESMultiblockGetSubSNES_Default);
885: return 0;
886: }