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