1 #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/
2 #include <petscdmcomposite.h>
3
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 };
15
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;
25
SNESReset_Multiblock(SNES snes)26 PetscErrorCode SNESReset_Multiblock(SNES snes)
27 {
28 SNES_Multiblock *mb = (SNES_Multiblock*) snes->data;
29 BlockDesc blocks = mb->blocks, next;
30 PetscErrorCode ierr;
31
32 PetscFunctionBegin;
33 while (blocks) {
34 ierr = SNESReset(blocks->snes);CHKERRQ(ierr);
35 #if 0
36 ierr = VecDestroy(&blocks->x);CHKERRQ(ierr);
37 #endif
38 ierr = VecScatterDestroy(&blocks->sctx);CHKERRQ(ierr);
39 ierr = ISDestroy(&blocks->is);CHKERRQ(ierr);
40 next = blocks->next;
41 blocks = next;
42 }
43 PetscFunctionReturn(0);
44 }
45
46 /*
47 SNESDestroy_Multiblock - Destroys the private SNES_Multiblock context that was created with SNESCreate_Multiblock().
48
49 Input Parameter:
50 . snes - the SNES context
51
52 Application Interface Routine: SNESDestroy()
53 */
SNESDestroy_Multiblock(SNES snes)54 PetscErrorCode SNESDestroy_Multiblock(SNES snes)
55 {
56 SNES_Multiblock *mb = (SNES_Multiblock*) snes->data;
57 BlockDesc blocks = mb->blocks, next;
58 PetscErrorCode ierr;
59
60 PetscFunctionBegin;
61 ierr = SNESReset_Multiblock(snes);CHKERRQ(ierr);
62 while (blocks) {
63 next = blocks->next;
64 ierr = SNESDestroy(&blocks->snes);CHKERRQ(ierr);
65 ierr = PetscFree(blocks->name);CHKERRQ(ierr);
66 ierr = PetscFree(blocks->fields);CHKERRQ(ierr);
67 ierr = PetscFree(blocks);CHKERRQ(ierr);
68 blocks = next;
69 }
70 ierr = PetscFree(snes->data);CHKERRQ(ierr);
71 PetscFunctionReturn(0);
72 }
73
74 /* Precondition: blocksize is set to a meaningful value */
SNESMultiblockSetFieldsRuntime_Private(SNES snes)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;
83
84 PetscFunctionBegin;
85 ierr = PetscMalloc1(mb->bs, &ifields);CHKERRQ(ierr);
86 for (i = 0;; ++i) {
87 ierr = PetscSNPrintf(name, sizeof(name), "%D", i);CHKERRQ(ierr);
88 ierr = PetscSNPrintf(optionname, sizeof(optionname), "-snes_multiblock_%D_fields", i);CHKERRQ(ierr);
89 nfields = mb->bs;
90 ierr = PetscOptionsGetIntArray(NULL,((PetscObject) snes)->prefix, optionname, ifields, &nfields, &flg);CHKERRQ(ierr);
91 if (!flg) break;
92 if (!nfields) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot list zero fields");
93 ierr = SNESMultiblockSetFields(snes, name, nfields, ifields);CHKERRQ(ierr);
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 ierr = PetscFree(ifields);CHKERRQ(ierr);
102 PetscFunctionReturn(0);
103 }
104
SNESMultiblockSetDefaults(SNES snes)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;
111
112 PetscFunctionBegin;
113 if (!blocks) {
114 if (snes->dm) {
115 PetscBool dmcomposite;
116
117 ierr = PetscObjectTypeCompare((PetscObject) snes->dm, DMCOMPOSITE, &dmcomposite);CHKERRQ(ierr);
118 if (dmcomposite) {
119 PetscInt nDM;
120 IS *fields;
121
122 ierr = PetscInfo(snes,"Setting up physics based multiblock solver using the embedded DM\n");CHKERRQ(ierr);
123 ierr = DMCompositeGetNumberDM(snes->dm, &nDM);CHKERRQ(ierr);
124 ierr = DMCompositeGetGlobalISs(snes->dm, &fields);CHKERRQ(ierr);
125 for (i = 0; i < nDM; ++i) {
126 char name[8];
127
128 ierr = PetscSNPrintf(name, sizeof(name), "%D", i);CHKERRQ(ierr);
129 ierr = SNESMultiblockSetIS(snes, name, fields[i]);CHKERRQ(ierr);
130 ierr = ISDestroy(&fields[i]);CHKERRQ(ierr);
131 }
132 ierr = PetscFree(fields);CHKERRQ(ierr);
133 }
134 } else {
135 PetscBool flg = PETSC_FALSE;
136 PetscBool stokes = PETSC_FALSE;
137
138 if (mb->bs <= 0) {
139 if (snes->jacobian_pre) {
140 ierr = MatGetBlockSize(snes->jacobian_pre, &mb->bs);CHKERRQ(ierr);
141 } else mb->bs = 1;
142 }
143
144 ierr = PetscOptionsGetBool(NULL,((PetscObject) snes)->prefix, "-snes_multiblock_default", &flg, NULL);CHKERRQ(ierr);
145 ierr = PetscOptionsGetBool(NULL,((PetscObject) snes)->prefix, "-snes_multiblock_detect_saddle_point", &stokes, NULL);CHKERRQ(ierr);
146 if (stokes) {
147 IS zerodiags, rest;
148 PetscInt nmin, nmax;
149
150 ierr = MatGetOwnershipRange(snes->jacobian_pre, &nmin, &nmax);CHKERRQ(ierr);
151 ierr = MatFindZeroDiagonals(snes->jacobian_pre, &zerodiags);CHKERRQ(ierr);
152 ierr = ISComplement(zerodiags, nmin, nmax, &rest);CHKERRQ(ierr);
153 ierr = SNESMultiblockSetIS(snes, "0", rest);CHKERRQ(ierr);
154 ierr = SNESMultiblockSetIS(snes, "1", zerodiags);CHKERRQ(ierr);
155 ierr = ISDestroy(&zerodiags);CHKERRQ(ierr);
156 ierr = ISDestroy(&rest);CHKERRQ(ierr);
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 ierr = SNESMultiblockSetFieldsRuntime_Private(snes);CHKERRQ(ierr);
162 if (mb->defined) {ierr = PetscInfo(snes, "Blocks defined using the options database\n");CHKERRQ(ierr);}
163 }
164 if (flg || !mb->defined) {
165 ierr = PetscInfo(snes, "Using default splitting of fields\n");CHKERRQ(ierr);
166 for (i = 0; i < mb->bs; ++i) {
167 char name[8];
168
169 ierr = PetscSNPrintf(name, sizeof(name), "%D", i);CHKERRQ(ierr);
170 ierr = SNESMultiblockSetFields(snes, name, 1, &i);CHKERRQ(ierr);
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;
180
181 ierr = MatGetOwnershipRange(snes->jacobian_pre, &nmin, &nmax);CHKERRQ(ierr);
182 ierr = ISComplement(blocks->is, nmin, nmax, &is2);CHKERRQ(ierr);
183 ierr = SNESMultiblockSetIS(snes, "1", is2);CHKERRQ(ierr);
184 ierr = ISDestroy(&is2);CHKERRQ(ierr);
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 PetscFunctionReturn(0);
189 }
190
191 /*
192 SNESSetUp_Multiblock - Sets up the internal data structures for the later use
193 of the SNESMULTIBLOCK nonlinear solver.
194
195 Input Parameters:
196 + snes - the SNES context
197 - x - the solution vector
198
199 Application Interface Routine: SNESSetUp()
200 */
SNESSetUp_Multiblock(SNES snes)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;
207
208 PetscFunctionBegin;
209 ierr = SNESMultiblockSetDefaults(snes);CHKERRQ(ierr);
210 numBlocks = mb->numBlocks;
211 blocks = mb->blocks;
212
213 /* Create ISs */
214 if (!mb->issetup) {
215 PetscInt ccsize, rstart, rend, nslots, bs;
216 PetscBool sorted;
217
218 mb->issetup = PETSC_TRUE;
219 bs = mb->bs;
220 ierr = MatGetOwnershipRange(snes->jacobian_pre, &rstart, &rend);CHKERRQ(ierr);
221 ierr = MatGetLocalSize(snes->jacobian_pre, NULL, &ccsize);CHKERRQ(ierr);
222 nslots = (rend - rstart)/bs;
223 for (i = 0; i < numBlocks; ++i) {
224 if (mb->defaultblocks) {
225 ierr = ISCreateStride(PetscObjectComm((PetscObject)snes), nslots, rstart+i, numBlocks, &blocks->is);CHKERRQ(ierr);
226 } else if (!blocks->is) {
227 if (blocks->nfields > 1) {
228 PetscInt *ii, j, k, nfields = blocks->nfields, *fields = blocks->fields;
229
230 ierr = PetscMalloc1(nfields*nslots, &ii);CHKERRQ(ierr);
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 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)snes), nslots*nfields, ii, PETSC_OWN_POINTER, &blocks->is);CHKERRQ(ierr);
237 } else {
238 ierr = ISCreateStride(PetscObjectComm((PetscObject)snes), nslots, rstart+blocks->fields[0], bs, &blocks->is);CHKERRQ(ierr);
239 }
240 }
241 ierr = ISSorted(blocks->is, &sorted);CHKERRQ(ierr);
242 if (!sorted) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Fields must be sorted when creating split");
243 blocks = blocks->next;
244 }
245 }
246
247 #if 0
248 /* Create matrices */
249 ilink = jac->head;
250 if (!jac->pmat) {
251 ierr = PetscMalloc1(nsplit,&jac->pmat);CHKERRQ(ierr);
252 for (i=0; i<nsplit; i++) {
253 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is,MAT_INITIAL_MATRIX,&jac->pmat[i]);CHKERRQ(ierr);
254 ilink = ilink->next;
255 }
256 } else {
257 for (i=0; i<nsplit; i++) {
258 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is,MAT_REUSE_MATRIX,&jac->pmat[i]);CHKERRQ(ierr);
259 ilink = ilink->next;
260 }
261 }
262 if (jac->realdiagonal) {
263 ilink = jac->head;
264 if (!jac->mat) {
265 ierr = PetscMalloc1(nsplit,&jac->mat);CHKERRQ(ierr);
266 for (i=0; i<nsplit; i++) {
267 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ilink->is,MAT_INITIAL_MATRIX,&jac->mat[i]);CHKERRQ(ierr);
268 ilink = ilink->next;
269 }
270 } else {
271 for (i=0; i<nsplit; i++) {
272 if (jac->mat[i]) {ierr = MatCreateSubMatrix(pc->mat,ilink->is,ilink->is,MAT_REUSE_MATRIX,&jac->mat[i]);CHKERRQ(ierr);}
273 ilink = ilink->next;
274 }
275 }
276 } else jac->mat = jac->pmat;
277 #endif
278
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 ierr = PetscMalloc1(nsplit,&jac->Afield);CHKERRQ(ierr);
285 for (i=0; i<nsplit; i++) {
286 ierr = MatCreateSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);CHKERRQ(ierr);
287 ilink = ilink->next;
288 }
289 } else {
290 for (i=0; i<nsplit; i++) {
291 ierr = MatCreateSubMatrix(pc->mat,ilink->is,NULL,MAT_REUSE_MATRIX,&jac->Afield[i]);CHKERRQ(ierr);
292 ilink = ilink->next;
293 }
294 }
295 }
296 #endif
297
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");
303
304 /* When extracting off-diagonal submatrices, we take complements from this range */
305 ierr = MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);CHKERRQ(ierr);
306
307 /* need to handle case when one is resetting up the preconditioner */
308 if (jac->schur) {
309 ilink = jac->head;
310 ierr = ISComplement(ilink->is,rstart,rend,&ccis);CHKERRQ(ierr);
311 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);CHKERRQ(ierr);
312 ierr = ISDestroy(&ccis);CHKERRQ(ierr);
313 ilink = ilink->next;
314 ierr = ISComplement(ilink->is,rstart,rend,&ccis);CHKERRQ(ierr);
315 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);CHKERRQ(ierr);
316 ierr = ISDestroy(&ccis);CHKERRQ(ierr);
317 ierr = MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->pmat[1]);CHKERRQ(ierr);
318 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),pc->flag);CHKERRQ(ierr);
319
320 } else {
321 KSP ksp;
322 char schurprefix[256];
323
324 /* extract the A01 and A10 matrices */
325 ilink = jac->head;
326 ierr = ISComplement(ilink->is,rstart,rend,&ccis);CHKERRQ(ierr);
327 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr);
328 ierr = ISDestroy(&ccis);CHKERRQ(ierr);
329 ilink = ilink->next;
330 ierr = ISComplement(ilink->is,rstart,rend,&ccis);CHKERRQ(ierr);
331 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr);
332 ierr = ISDestroy(&ccis);CHKERRQ(ierr);
333 /* Use mat[0] (diagonal block of the real matrix) preconditioned by pmat[0] */
334 ierr = MatCreateSchurComplement(jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1],&jac->schur);CHKERRQ(ierr);
335 /* set tabbing and options prefix of KSP inside the MatSchur */
336 ierr = MatSchurComplementGetKSP(jac->schur,&ksp);CHKERRQ(ierr);
337 ierr = PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,2);CHKERRQ(ierr);
338 ierr = PetscSNPrintf(schurprefix,sizeof(schurprefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",jac->head->splitname);CHKERRQ(ierr);
339 ierr = KSPSetOptionsPrefix(ksp,schurprefix);CHKERRQ(ierr);
340 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr);
341
342 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);CHKERRQ(ierr);
343 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);CHKERRQ(ierr);
344 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);CHKERRQ(ierr);
345 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));CHKERRQ(ierr);
346 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) {
347 PC pc;
348 ierr = KSPGetPC(jac->kspschur,&pc);CHKERRQ(ierr);
349 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr);
350 /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */
351 }
352 ierr = PetscSNPrintf(schurprefix,sizeof(schurprefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr);
353 ierr = KSPSetOptionsPrefix(jac->kspschur,schurprefix);CHKERRQ(ierr);
354 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
355 ierr = KSPSetFromOptions(jac->kspschur);CHKERRQ(ierr);
356
357 ierr = PetscMalloc2(2,&jac->x,2,&jac->y);CHKERRQ(ierr);
358 ierr = MatCreateVecs(jac->pmat[0],&jac->x[0],&jac->y[0]);CHKERRQ(ierr);
359 ierr = MatCreateVecs(jac->pmat[1],&jac->x[1],&jac->y[1]);CHKERRQ(ierr);
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 /*ierr = SNESSetFunction(blocks->snes, blocks->x, func);CHKERRQ(ierr);*/
373 /*ierr = SNESSetJacobian(blocks->snes, blocks->x, jac);CHKERRQ(ierr);*/
374 ierr = VecDuplicate(blocks->snes->vec_sol, &blocks->x);CHKERRQ(ierr);
375 /* really want setfromoptions called in SNESSetFromOptions_Multiblock(), but it is not ready yet */
376 ierr = SNESSetFromOptions(blocks->snes);CHKERRQ(ierr);
377 ierr = SNESSetUp(blocks->snes);CHKERRQ(ierr);
378 blocks = blocks->next;
379 i++;
380 }
381 }
382
383 /* Compute scatter contexts needed by multiplicative versions and non-default splits */
384 if (!mb->blocks->sctx) {
385 Vec xtmp;
386
387 blocks = mb->blocks;
388 ierr = MatCreateVecs(snes->jacobian_pre, &xtmp, NULL);CHKERRQ(ierr);
389 while (blocks) {
390 ierr = VecScatterCreate(xtmp, blocks->is, blocks->x, NULL, &blocks->sctx);CHKERRQ(ierr);
391 blocks = blocks->next;
392 }
393 ierr = VecDestroy(&xtmp);CHKERRQ(ierr);
394 }
395 PetscFunctionReturn(0);
396 }
397
398 /*
399 SNESSetFromOptions_Multiblock - Sets various parameters for the SNESMULTIBLOCK method.
400
401 Input Parameter:
402 . snes - the SNES context
403
404 Application Interface Routine: SNESSetFromOptions()
405 */
SNESSetFromOptions_Multiblock(PetscOptionItems * PetscOptionsObject,SNES snes)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;
413
414 PetscFunctionBegin;
415 ierr = PetscOptionsHead(PetscOptionsObject,"SNES Multiblock options");CHKERRQ(ierr);
416 ierr = PetscOptionsInt("-snes_multiblock_block_size", "Blocksize that defines number of fields", "PCFieldSplitSetBlockSize", mb->bs, &bs, &flg);CHKERRQ(ierr);
417 if (flg) {ierr = SNESMultiblockSetBlockSize(snes, bs);CHKERRQ(ierr);}
418 ierr = PetscOptionsEnum("-snes_multiblock_type", "Type of composition", "PCFieldSplitSetType", PCCompositeTypes, (PetscEnum) mb->type, (PetscEnum*) &ctype, &flg);CHKERRQ(ierr);
419 if (flg) {
420 ierr = SNESMultiblockSetType(snes,ctype);CHKERRQ(ierr);
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 ierr = SNESMultiblockSetFieldsRuntime_Private(snes);CHKERRQ(ierr);
426 if (mb->defined) {ierr = PetscInfo(snes, "Blocks defined using the options database\n");CHKERRQ(ierr);}
427 }
428 ierr = PetscOptionsTail();CHKERRQ(ierr);
429 PetscFunctionReturn(0);
430 }
431
432 /*
433 SNESView_Multiblock - Prints info from the SNESMULTIBLOCK data structure.
434
435 Input Parameters:
436 + SNES - the SNES context
437 - viewer - visualization context
438
439 Application Interface Routine: SNESView()
440 */
SNESView_Multiblock(SNES snes,PetscViewer viewer)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;
447
448 PetscFunctionBegin;
449 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
450 if (iascii) {
451 ierr = PetscViewerASCIIPrintf(viewer," Multiblock with %s composition: total blocks = %D, blocksize = %D\n", PCCompositeTypes[mb->type], mb->numBlocks, mb->bs);CHKERRQ(ierr);
452 ierr = PetscViewerASCIIPrintf(viewer," Solver info for each split is in the following SNES objects:\n");CHKERRQ(ierr);
453 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
454 while (blocks) {
455 if (blocks->fields) {
456 PetscInt j;
457
458 ierr = PetscViewerASCIIPrintf(viewer, " Block %s Fields ", blocks->name);CHKERRQ(ierr);
459 ierr = PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);CHKERRQ(ierr);
460 for (j = 0; j < blocks->nfields; ++j) {
461 if (j > 0) {
462 ierr = PetscViewerASCIIPrintf(viewer, ",");CHKERRQ(ierr);
463 }
464 ierr = PetscViewerASCIIPrintf(viewer, " %D", blocks->fields[j]);CHKERRQ(ierr);
465 }
466 ierr = PetscViewerASCIIPrintf(viewer, "\n");CHKERRQ(ierr);
467 ierr = PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);CHKERRQ(ierr);
468 } else {
469 ierr = PetscViewerASCIIPrintf(viewer, " Block %s Defined by IS\n", blocks->name);CHKERRQ(ierr);
470 }
471 ierr = SNESView(blocks->snes, viewer);CHKERRQ(ierr);
472 blocks = blocks->next;
473 }
474 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
475 }
476 PetscFunctionReturn(0);
477 }
478
479 /*
480 SNESSolve_Multiblock - Solves a nonlinear system with the Multiblock method.
481
482 Input Parameters:
483 . snes - the SNES context
484
485 Output Parameter:
486 . outits - number of iterations until termination
487
488 Application Interface Routine: SNESSolve()
489 */
SNESSolve_Multiblock(SNES snes)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;
497
498 PetscFunctionBegin;
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);
500
501 snes->reason = SNES_CONVERGED_ITERATING;
502
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 */
507
508 ierr = VecSetBlockSize(X, mb->bs);CHKERRQ(ierr);
509 ierr = VecSetBlockSize(Y, mb->bs);CHKERRQ(ierr);
510 ierr = VecSetBlockSize(F, mb->bs);CHKERRQ(ierr);
511 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
512 snes->iter = 0;
513 snes->norm = 0.;
514 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
515
516 if (!snes->vec_func_init_set) {
517 ierr = SNESComputeFunction(snes, X, F);CHKERRQ(ierr);
518 } else snes->vec_func_init_set = PETSC_FALSE;
519
520 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */
521 SNESCheckFunctionNorm(snes,fnorm);
522 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
523 snes->norm = fnorm;
524 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
525 ierr = SNESLogConvergenceHistory(snes,fnorm,0);CHKERRQ(ierr);
526 ierr = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr);
527
528 /* test convergence */
529 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
530 if (snes->reason) PetscFunctionReturn(0);
531
532 for (i = 0; i < maxits; i++) {
533 /* Call general purpose update function */
534 if (snes->ops->update) {
535 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr);
536 }
537 /* Compute X^{new} from subsolves */
538 if (mb->type == PC_COMPOSITE_ADDITIVE) {
539 BlockDesc blocks = mb->blocks;
540
541 if (mb->defaultblocks) {
542 /*TODO: Make an array of Vecs for this */
543 /*ierr = VecStrideGatherAll(X, mb->x, INSERT_VALUES);CHKERRQ(ierr);*/
544 while (blocks) {
545 ierr = SNESSolve(blocks->snes, NULL, blocks->x);CHKERRQ(ierr);
546 blocks = blocks->next;
547 }
548 /*ierr = VecStrideScatterAll(mb->x, X, INSERT_VALUES);CHKERRQ(ierr);*/
549 } else {
550 while (blocks) {
551 ierr = VecScatterBegin(blocks->sctx, X, blocks->x, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr);
552 ierr = VecScatterEnd(blocks->sctx, X, blocks->x, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr);
553 ierr = SNESSolve(blocks->snes, NULL, blocks->x);CHKERRQ(ierr);
554 ierr = VecScatterBegin(blocks->sctx, blocks->x, X, INSERT_VALUES, SCATTER_REVERSE);CHKERRQ(ierr);
555 ierr = VecScatterEnd(blocks->sctx, blocks->x, X, INSERT_VALUES, SCATTER_REVERSE);CHKERRQ(ierr);
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 ierr = SNESComputeFunction(snes, X, F);CHKERRQ(ierr);
562 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr);
563 SNESCheckFunctionNorm(snes,fnorm);
564
565 if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >=0) {
566 snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
567 break;
568 }
569
570 /* Monitor convergence */
571 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
572 snes->iter = i+1;
573 snes->norm = fnorm;
574 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
575 ierr = SNESLogConvergenceHistory(snes,snes->norm,0);CHKERRQ(ierr);
576 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr);
577 /* Test for convergence */
578 ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
579 if (snes->reason) break;
580 }
581 if (i == maxits) {
582 ierr = PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", maxits);CHKERRQ(ierr);
583 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
584 }
585 PetscFunctionReturn(0);
586 }
587
SNESMultiblockSetFields_Default(SNES snes,const char name[],PetscInt n,const PetscInt fields[])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;
595
596 PetscFunctionBegin;
597 if (mb->defined) {
598 ierr = PetscInfo1(snes, "Ignoring new block \"%s\" because the blocks have already been defined\n", name);CHKERRQ(ierr);
599 PetscFunctionReturn(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 ierr = PetscNew(&newblock);CHKERRQ(ierr);
606 if (name) {
607 ierr = PetscStrallocpy(name, &newblock->name);CHKERRQ(ierr);
608 } else {
609 PetscInt len = floor(log10(mb->numBlocks))+1;
610
611 ierr = PetscMalloc1(len+1, &newblock->name);CHKERRQ(ierr);
612 ierr = PetscSNPrintf(newblock->name, len, "%s", mb->numBlocks);CHKERRQ(ierr);
613 }
614 newblock->nfields = n;
615
616 ierr = PetscMalloc1(n, &newblock->fields);CHKERRQ(ierr);
617 ierr = PetscArraycpy(newblock->fields, fields, n);CHKERRQ(ierr);
618
619 newblock->next = NULL;
620
621 ierr = SNESCreate(PetscObjectComm((PetscObject)snes), &newblock->snes);CHKERRQ(ierr);
622 ierr = PetscObjectIncrementTabLevel((PetscObject) newblock->snes, (PetscObject) snes, 1);CHKERRQ(ierr);
623 ierr = SNESSetType(newblock->snes, SNESNRICHARDSON);CHKERRQ(ierr);
624 ierr = PetscLogObjectParent((PetscObject) snes, (PetscObject) newblock->snes);CHKERRQ(ierr);
625 ierr = PetscSNPrintf(prefix, sizeof(prefix), "%smultiblock_%s_", ((PetscObject) snes)->prefix ? ((PetscObject) snes)->prefix : "", newblock->name);CHKERRQ(ierr);
626 ierr = SNESSetOptionsPrefix(newblock->snes, prefix);CHKERRQ(ierr);
627
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 PetscFunctionReturn(0);
640 }
641
SNESMultiblockSetIS_Default(SNES snes,const char name[],IS is)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;
648
649 PetscFunctionBegin;
650 if (mb->defined) {
651 ierr = PetscInfo1(snes, "Ignoring new block \"%s\" because the blocks have already been defined\n", name);CHKERRQ(ierr);
652 PetscFunctionReturn(0);
653 }
654 ierr = PetscNew(&newblock);CHKERRQ(ierr);
655 if (name) {
656 ierr = PetscStrallocpy(name, &newblock->name);CHKERRQ(ierr);
657 } else {
658 PetscInt len = floor(log10(mb->numBlocks))+1;
659
660 ierr = PetscMalloc1(len+1, &newblock->name);CHKERRQ(ierr);
661 ierr = PetscSNPrintf(newblock->name, len, "%s", mb->numBlocks);CHKERRQ(ierr);
662 }
663 newblock->is = is;
664
665 ierr = PetscObjectReference((PetscObject) is);CHKERRQ(ierr);
666
667 newblock->next = NULL;
668
669 ierr = SNESCreate(PetscObjectComm((PetscObject)snes), &newblock->snes);CHKERRQ(ierr);
670 ierr = PetscObjectIncrementTabLevel((PetscObject) newblock->snes, (PetscObject) snes, 1);CHKERRQ(ierr);
671 ierr = SNESSetType(newblock->snes, SNESNRICHARDSON);CHKERRQ(ierr);
672 ierr = PetscLogObjectParent((PetscObject) snes, (PetscObject) newblock->snes);CHKERRQ(ierr);
673 ierr = PetscSNPrintf(prefix, sizeof(prefix), "%smultiblock_%s_", ((PetscObject) snes)->prefix ? ((PetscObject) snes)->prefix : "", newblock->name);CHKERRQ(ierr);
674 ierr = SNESSetOptionsPrefix(newblock->snes, prefix);CHKERRQ(ierr);
675
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 PetscFunctionReturn(0);
688 }
689
SNESMultiblockSetBlockSize_Default(SNES snes,PetscInt bs)690 PetscErrorCode SNESMultiblockSetBlockSize_Default(SNES snes, PetscInt bs)
691 {
692 SNES_Multiblock *mb = (SNES_Multiblock*) snes->data;
693
694 PetscFunctionBegin;
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 PetscFunctionReturn(0);
699 }
700
SNESMultiblockGetSubSNES_Default(SNES snes,PetscInt * n,SNES ** subsnes)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;
707
708 PetscFunctionBegin;
709 ierr = PetscMalloc1(mb->numBlocks, subsnes);CHKERRQ(ierr);
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);
715
716 if (n) *n = mb->numBlocks;
717 PetscFunctionReturn(0);
718 }
719
SNESMultiblockSetType_Default(SNES snes,PCCompositeType type)720 PetscErrorCode SNESMultiblockSetType_Default(SNES snes, PCCompositeType type)
721 {
722 SNES_Multiblock *mb = (SNES_Multiblock*) snes->data;
723 PetscErrorCode ierr;
724
725 PetscFunctionBegin;
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;
733
734 ierr = PetscObjectComposeFunction((PetscObject) snes, "SNESMultiblockGetSubSNES_C", SNESMultiblockGetSubSNES_Schur);CHKERRQ(ierr);
735 ierr = PetscObjectComposeFunction((PetscObject) snes, "SNESMultiblockSchurPrecondition_C", SNESMultiblockSchurPrecondition_Default);CHKERRQ(ierr);
736 #endif
737 } else {
738 snes->ops->solve = SNESSolve_Multiblock;
739 snes->ops->view = SNESView_Multiblock;
740
741 ierr = PetscObjectComposeFunction((PetscObject) snes, "SNESMultiblockGetSubSNES_C", SNESMultiblockGetSubSNES_Default);CHKERRQ(ierr);
742 ierr = PetscObjectComposeFunction((PetscObject) snes, "SNESMultiblockSchurPrecondition_C", 0);CHKERRQ(ierr);
743 }
744 PetscFunctionReturn(0);
745 }
746
747 /*@
748 SNESMultiblockSetFields - Sets the fields for one particular block in the solver
749
750 Logically Collective on SNES
751
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
757
758 Level: intermediate
759
760 Notes:
761 Use SNESMultiblockSetIS() to set a completely general set of row indices as a block.
762
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.
768
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_.
771
772 .seealso: SNESMultiblockGetSubSNES(), SNESMULTIBLOCK, SNESMultiblockSetBlockSize(), SNESMultiblockSetIS()
773 @*/
SNESMultiblockSetFields(SNES snes,const char name[],PetscInt n,const PetscInt * fields)774 PetscErrorCode SNESMultiblockSetFields(SNES snes, const char name[], PetscInt n, const PetscInt *fields)
775 {
776 PetscErrorCode ierr;
777
778 PetscFunctionBegin;
779 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
780 PetscValidCharPointer(name, 2);
781 if (n < 1) SETERRQ2(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_OUTOFRANGE, "Provided number of fields %D in split \"%s\" not positive", n, name);
782 PetscValidIntPointer(fields, 3);
783 ierr = PetscTryMethod(snes, "SNESMultiblockSetFields_C", (SNES, const char[], PetscInt, const PetscInt*), (snes, name, n, fields));CHKERRQ(ierr);
784 PetscFunctionReturn(0);
785 }
786
787 /*@
788 SNESMultiblockSetIS - Sets the global row indices for the block
789
790 Logically Collective on SNES
791
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
796
797 Notes:
798 Use SNESMultiblockSetFields(), for blocks defined by strides.
799
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_.
802
803 Level: intermediate
804
805 .seealso: SNESMultiblockGetSubSNES(), SNESMULTIBLOCK, SNESMultiblockSetBlockSize()
806 @*/
SNESMultiblockSetIS(SNES snes,const char name[],IS is)807 PetscErrorCode SNESMultiblockSetIS(SNES snes, const char name[], IS is)
808 {
809 PetscErrorCode ierr;
810
811 PetscFunctionBegin;
812 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
813 PetscValidCharPointer(name, 2);
814 PetscValidHeaderSpecific(is, IS_CLASSID, 3);
815 ierr = PetscTryMethod(snes, "SNESMultiblockSetIS_C", (SNES, const char[], IS), (snes, name, is));CHKERRQ(ierr);
816 PetscFunctionReturn(0);
817 }
818
819 /*@
820 SNESMultiblockSetType - Sets the type of block combination.
821
822 Collective on SNES
823
824 Input Parameters:
825 + snes - the solver context
826 - type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE
827
828 Options Database Key:
829 . -snes_multiblock_type <type: one of multiplicative, additive, symmetric_multiplicative> - Sets block combination type
830
831 Level: Developer
832
833 .seealso: PCCompositeSetType()
834 @*/
SNESMultiblockSetType(SNES snes,PCCompositeType type)835 PetscErrorCode SNESMultiblockSetType(SNES snes, PCCompositeType type)
836 {
837 PetscErrorCode ierr;
838
839 PetscFunctionBegin;
840 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
841 ierr = PetscTryMethod(snes, "SNESMultiblockSetType_C", (SNES, PCCompositeType), (snes, type));CHKERRQ(ierr);
842 PetscFunctionReturn(0);
843 }
844
845 /*@
846 SNESMultiblockSetBlockSize - Sets the block size for structured mesh block division. If not set the matrix block size is used.
847
848 Logically Collective on SNES
849
850 Input Parameters:
851 + snes - the solver context
852 - bs - the block size
853
854 Level: intermediate
855
856 .seealso: SNESMultiblockGetSubSNES(), SNESMULTIBLOCK, SNESMultiblockSetFields()
857 @*/
SNESMultiblockSetBlockSize(SNES snes,PetscInt bs)858 PetscErrorCode SNESMultiblockSetBlockSize(SNES snes, PetscInt bs)
859 {
860 PetscErrorCode ierr;
861
862 PetscFunctionBegin;
863 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
864 PetscValidLogicalCollectiveInt(snes, bs, 2);
865 ierr = PetscTryMethod(snes, "SNESMultiblockSetBlockSize_C", (SNES, PetscInt), (snes,bs));CHKERRQ(ierr);
866 PetscFunctionReturn(0);
867 }
868
869 /*@C
870 SNESMultiblockGetSubSNES - Gets the SNES contexts for all blocks
871
872 Collective on SNES
873
874 Input Parameter:
875 . snes - the solver context
876
877 Output Parameters:
878 + n - the number of blocks
879 - subsnes - the array of SNES contexts
880
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).
884
885 You must call SNESSetUp() before calling SNESMultiblockGetSubSNES().
886
887 Level: advanced
888
889 .seealso: SNESMULTIBLOCK
890 @*/
SNESMultiblockGetSubSNES(SNES snes,PetscInt * n,SNES * subsnes[])891 PetscErrorCode SNESMultiblockGetSubSNES(SNES snes, PetscInt *n, SNES *subsnes[])
892 {
893 PetscErrorCode ierr;
894
895 PetscFunctionBegin;
896 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
897 if (n) PetscValidIntPointer(n, 2);
898 ierr = PetscUseMethod(snes, "SNESMultiblockGetSubSNES_C", (SNES, PetscInt*, SNES **), (snes, n, subsnes));CHKERRQ(ierr);
899 PetscFunctionReturn(0);
900 }
901
902 /*MC
903 SNESMULTIBLOCK - Multiblock nonlinear solver that can use overlapping or nonoverlapping blocks, organized
904 additively (Jacobi) or multiplicatively (Gauss-Seidel).
905
906 Level: beginner
907
908 .seealso: SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESNEWTONTR, SNESNRICHARDSON
909 M*/
SNESCreate_Multiblock(SNES snes)910 PETSC_EXTERN PetscErrorCode SNESCreate_Multiblock(SNES snes)
911 {
912 SNES_Multiblock *mb;
913 PetscErrorCode ierr;
914
915 PetscFunctionBegin;
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;
922
923 snes->usesksp = PETSC_FALSE;
924 snes->usesnpc = PETSC_FALSE;
925
926 snes->alwayscomputesfinalresidual = PETSC_TRUE;
927
928 ierr = PetscNewLog(snes,&mb);CHKERRQ(ierr);
929 snes->data = (void*) mb;
930 mb->defined = PETSC_FALSE;
931 mb->numBlocks = 0;
932 mb->bs = -1;
933 mb->type = PC_COMPOSITE_MULTIPLICATIVE;
934
935 /* We attach functions so that they can be called on another PC without crashing the program */
936 ierr = PetscObjectComposeFunction((PetscObject) snes, "SNESMultiblockSetFields_C", SNESMultiblockSetFields_Default);CHKERRQ(ierr);
937 ierr = PetscObjectComposeFunction((PetscObject) snes, "SNESMultiblockSetIS_C", SNESMultiblockSetIS_Default);CHKERRQ(ierr);
938 ierr = PetscObjectComposeFunction((PetscObject) snes, "SNESMultiblockSetType_C", SNESMultiblockSetType_Default);CHKERRQ(ierr);
939 ierr = PetscObjectComposeFunction((PetscObject) snes, "SNESMultiblockSetBlockSize_C", SNESMultiblockSetBlockSize_Default);CHKERRQ(ierr);
940 ierr = PetscObjectComposeFunction((PetscObject) snes, "SNESMultiblockGetSubSNES_C", SNESMultiblockGetSubSNES_Default);CHKERRQ(ierr);
941 PetscFunctionReturn(0);
942 }
943