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