1 /*
2 Defines a SNES that can consist of a collection of SNESes
3 */
4 #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/
5 #include <petscblaslapack.h>
6
7 const char *const SNESCompositeTypes[] = {"ADDITIVE", "MULTIPLICATIVE", "ADDITIVEOPTIMAL", "SNESCompositeType", "SNES_COMPOSITE", NULL};
8
9 typedef struct _SNES_CompositeLink *SNES_CompositeLink;
10 struct _SNES_CompositeLink {
11 SNES snes;
12 PetscReal dmp;
13 Vec X;
14 SNES_CompositeLink next;
15 SNES_CompositeLink previous;
16 };
17
18 typedef struct {
19 SNES_CompositeLink head;
20 PetscInt nsnes;
21 SNESCompositeType type;
22 Vec Xorig;
23 PetscInt innerFailures; /* the number of inner failures we've seen */
24
25 /* context for ADDITIVEOPTIMAL */
26 Vec *Xes, *Fes; /* solution and residual vectors for the subsolvers */
27 PetscReal *fnorms; /* norms of the solutions */
28 PetscScalar *h; /* the matrix formed as q_ij = (rdot_i, rdot_j) */
29 PetscScalar *g; /* the dotproducts of the previous function with the candidate functions */
30 PetscBLASInt n; /* matrix dimension -- nsnes */
31 PetscBLASInt nrhs; /* the number of right-hand sides */
32 PetscBLASInt lda; /* the padded matrix dimension */
33 PetscBLASInt ldb; /* the padded vector dimension */
34 PetscReal *s; /* the singular values */
35 PetscScalar *beta; /* the RHS and combination */
36 PetscReal rcond; /* the exit condition */
37 PetscBLASInt rank; /* the effective rank */
38 PetscScalar *work; /* the work vector */
39 PetscReal *rwork; /* the real work vector used for complex */
40 PetscBLASInt lwork; /* the size of the work vector */
41 PetscBLASInt info; /* the output condition */
42
43 PetscReal rtol; /* restart tolerance for accepting the combination */
44 PetscReal stol; /* restart tolerance for the combination */
45 } SNES_Composite;
46
SNESCompositeApply_Multiplicative(SNES snes,Vec X,Vec B,Vec F,PetscReal * fnorm)47 static PetscErrorCode SNESCompositeApply_Multiplicative(SNES snes, Vec X, Vec B, Vec F, PetscReal *fnorm)
48 {
49 SNES_Composite *jac = (SNES_Composite *)snes->data;
50 SNES_CompositeLink next = jac->head;
51 Vec FSub;
52 SNESConvergedReason reason;
53
54 PetscFunctionBegin;
55 PetscCheck(next, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");
56 if (snes->normschedule == SNES_NORM_ALWAYS) PetscCall(SNESSetInitialFunction(next->snes, F));
57 PetscCall(SNESSolve(next->snes, B, X));
58 PetscCall(SNESGetConvergedReason(next->snes, &reason));
59 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
60 jac->innerFailures++;
61 if (jac->innerFailures >= snes->maxFailures) {
62 snes->reason = SNES_DIVERGED_INNER;
63 PetscFunctionReturn(PETSC_SUCCESS);
64 }
65 }
66
67 while (next->next) {
68 /* only copy the function over in the case where the functions correspond */
69 if (next->snes->npcside == PC_RIGHT && next->snes->normschedule != SNES_NORM_NONE) {
70 PetscCall(SNESGetFunction(next->snes, &FSub, NULL, NULL));
71 next = next->next;
72 PetscCall(SNESSetInitialFunction(next->snes, FSub));
73 } else {
74 next = next->next;
75 }
76 PetscCall(SNESSolve(next->snes, B, X));
77 PetscCall(SNESGetConvergedReason(next->snes, &reason));
78 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
79 jac->innerFailures++;
80 if (jac->innerFailures >= snes->maxFailures) {
81 snes->reason = SNES_DIVERGED_INNER;
82 PetscFunctionReturn(PETSC_SUCCESS);
83 }
84 }
85 }
86 if (next->snes->npcside == PC_RIGHT) {
87 PetscCall(SNESGetFunction(next->snes, &FSub, NULL, NULL));
88 PetscCall(VecCopy(FSub, F));
89 if (fnorm) {
90 if (snes->xl && snes->xu) {
91 PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm));
92 } else {
93 PetscCall(VecNorm(F, NORM_2, fnorm));
94 }
95 SNESCheckFunctionDomainError(snes, *fnorm);
96 }
97 } else if (snes->normschedule == SNES_NORM_ALWAYS) {
98 PetscCall(SNESComputeFunction(snes, X, F));
99 if (fnorm) {
100 if (snes->xl && snes->xu) {
101 PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm));
102 } else {
103 PetscCall(VecNorm(F, NORM_2, fnorm));
104 }
105 SNESCheckFunctionDomainError(snes, *fnorm);
106 }
107 }
108 PetscFunctionReturn(PETSC_SUCCESS);
109 }
110
SNESCompositeApply_Additive(SNES snes,Vec X,Vec B,Vec F,PetscReal * fnorm)111 static PetscErrorCode SNESCompositeApply_Additive(SNES snes, Vec X, Vec B, Vec F, PetscReal *fnorm)
112 {
113 SNES_Composite *jac = (SNES_Composite *)snes->data;
114 SNES_CompositeLink next = jac->head;
115 Vec Y, Xorig;
116 SNESConvergedReason reason;
117
118 PetscFunctionBegin;
119 Y = snes->vec_sol_update;
120 if (!jac->Xorig) PetscCall(VecDuplicate(X, &jac->Xorig));
121 Xorig = jac->Xorig;
122 PetscCall(VecCopy(X, Xorig));
123 PetscCheck(next, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");
124 if (snes->normschedule == SNES_NORM_ALWAYS) {
125 PetscCall(SNESSetInitialFunction(next->snes, F));
126 while (next->next) {
127 next = next->next;
128 PetscCall(SNESSetInitialFunction(next->snes, F));
129 }
130 }
131 next = jac->head;
132 PetscCall(VecCopy(Xorig, Y));
133 PetscCall(SNESSolve(next->snes, B, Y));
134 PetscCall(SNESGetConvergedReason(next->snes, &reason));
135 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
136 jac->innerFailures++;
137 if (jac->innerFailures >= snes->maxFailures) {
138 snes->reason = SNES_DIVERGED_INNER;
139 PetscFunctionReturn(PETSC_SUCCESS);
140 }
141 }
142 PetscCall(VecAXPY(Y, -1.0, Xorig));
143 PetscCall(VecAXPY(X, next->dmp, Y));
144 while (next->next) {
145 next = next->next;
146 PetscCall(VecCopy(Xorig, Y));
147 PetscCall(SNESSolve(next->snes, B, Y));
148 PetscCall(SNESGetConvergedReason(next->snes, &reason));
149 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
150 jac->innerFailures++;
151 if (jac->innerFailures >= snes->maxFailures) {
152 snes->reason = SNES_DIVERGED_INNER;
153 PetscFunctionReturn(PETSC_SUCCESS);
154 }
155 }
156 PetscCall(VecAXPY(Y, -1.0, Xorig));
157 PetscCall(VecAXPY(X, next->dmp, Y));
158 }
159 if (snes->normschedule == SNES_NORM_ALWAYS) {
160 PetscCall(SNESComputeFunction(snes, X, F));
161 if (fnorm) {
162 if (snes->xl && snes->xu) {
163 PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm));
164 } else {
165 PetscCall(VecNorm(F, NORM_2, fnorm));
166 }
167 SNESCheckFunctionDomainError(snes, *fnorm);
168 }
169 }
170 PetscFunctionReturn(PETSC_SUCCESS);
171 }
172
173 /* approximately solve the overdetermined system:
174
175 2*F(x_i)\cdot F(\x_j)\alpha_i = 0
176 \alpha_i = 1
177
178 Which minimizes the L2 norm of the linearization of:
179 ||F(\sum_i \alpha_i*x_i)||^2
180
181 With the constraint that \sum_i\alpha_i = 1
182 Where x_i is the solution from the ith subsolver.
183 */
SNESCompositeApply_AdditiveOptimal(SNES snes,Vec X,Vec B,Vec F,PetscReal * fnorm)184 static PetscErrorCode SNESCompositeApply_AdditiveOptimal(SNES snes, Vec X, Vec B, Vec F, PetscReal *fnorm)
185 {
186 SNES_Composite *jac = (SNES_Composite *)snes->data;
187 SNES_CompositeLink next = jac->head;
188 Vec *Xes = jac->Xes, *Fes = jac->Fes;
189 PetscInt i, j;
190 PetscScalar tot, total, ftf;
191 PetscReal min_fnorm;
192 PetscInt min_i;
193 SNESConvergedReason reason;
194
195 PetscFunctionBegin;
196 PetscCheck(next, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "No composite SNESes supplied via SNESCompositeAddSNES() or -snes_composite_sneses");
197
198 if (snes->normschedule == SNES_NORM_ALWAYS) {
199 next = jac->head;
200 PetscCall(SNESSetInitialFunction(next->snes, F));
201 while (next->next) {
202 next = next->next;
203 PetscCall(SNESSetInitialFunction(next->snes, F));
204 }
205 }
206
207 next = jac->head;
208 i = 0;
209 PetscCall(VecCopy(X, Xes[i]));
210 PetscCall(SNESSolve(next->snes, B, Xes[i]));
211 PetscCall(SNESGetConvergedReason(next->snes, &reason));
212 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
213 jac->innerFailures++;
214 if (jac->innerFailures >= snes->maxFailures) {
215 snes->reason = SNES_DIVERGED_INNER;
216 PetscFunctionReturn(PETSC_SUCCESS);
217 }
218 }
219 while (next->next) {
220 i++;
221 next = next->next;
222 PetscCall(VecCopy(X, Xes[i]));
223 PetscCall(SNESSolve(next->snes, B, Xes[i]));
224 PetscCall(SNESGetConvergedReason(next->snes, &reason));
225 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
226 jac->innerFailures++;
227 if (jac->innerFailures >= snes->maxFailures) {
228 snes->reason = SNES_DIVERGED_INNER;
229 PetscFunctionReturn(PETSC_SUCCESS);
230 }
231 }
232 }
233
234 /* all the solutions are collected; combine optimally */
235 for (i = 0; i < jac->n; i++) {
236 for (j = 0; j < i + 1; j++) PetscCall(VecDotBegin(Fes[i], Fes[j], &jac->h[i + j * jac->n]));
237 PetscCall(VecDotBegin(Fes[i], F, &jac->g[i]));
238 }
239
240 for (i = 0; i < jac->n; i++) {
241 for (j = 0; j < i + 1; j++) {
242 PetscCall(VecDotEnd(Fes[i], Fes[j], &jac->h[i + j * jac->n]));
243 if (i == j) jac->fnorms[i] = PetscSqrtReal(PetscRealPart(jac->h[i + j * jac->n]));
244 }
245 PetscCall(VecDotEnd(Fes[i], F, &jac->g[i]));
246 }
247
248 ftf = (*fnorm) * (*fnorm);
249
250 for (i = 0; i < jac->n; i++) {
251 for (j = i + 1; j < jac->n; j++) jac->h[i + j * jac->n] = jac->h[j + i * jac->n];
252 }
253
254 for (i = 0; i < jac->n; i++) {
255 for (j = 0; j < jac->n; j++) jac->h[i + j * jac->n] = jac->h[i + j * jac->n] - jac->g[j] - jac->g[i] + ftf;
256 jac->beta[i] = ftf - jac->g[i];
257 }
258
259 jac->info = 0;
260 jac->rcond = -1.;
261 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
262 #if defined(PETSC_USE_COMPLEX)
263 PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&jac->n, &jac->n, &jac->nrhs, jac->h, &jac->lda, jac->beta, &jac->lda, jac->s, &jac->rcond, &jac->rank, jac->work, &jac->lwork, jac->rwork, &jac->info));
264 #else
265 PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&jac->n, &jac->n, &jac->nrhs, jac->h, &jac->lda, jac->beta, &jac->lda, jac->s, &jac->rcond, &jac->rank, jac->work, &jac->lwork, &jac->info));
266 #endif
267 PetscCall(PetscFPTrapPop());
268 PetscCheck(jac->info >= 0, PetscObjectComm((PetscObject)snes), PETSC_ERR_LIB, "Bad argument to GELSS");
269 PetscCheck(jac->info <= 0, PetscObjectComm((PetscObject)snes), PETSC_ERR_LIB, "SVD failed to converge");
270 tot = 0.;
271 total = 0.;
272 for (i = 0; i < jac->n; i++) {
273 PetscCheck(!snes->errorifnotconverged || !PetscIsInfOrNanScalar(jac->beta[i]), PetscObjectComm((PetscObject)snes), PETSC_ERR_LIB, "SVD generated inconsistent output");
274 PetscCall(PetscInfo(snes, "%" PetscInt_FMT ": %g\n", i, (double)PetscRealPart(jac->beta[i])));
275 tot += jac->beta[i];
276 total += PetscAbsScalar(jac->beta[i]);
277 }
278 PetscCall(VecScale(X, 1. - tot));
279 PetscCall(VecMAXPY(X, jac->n, jac->beta, Xes));
280 PetscCall(SNESComputeFunction(snes, X, F));
281
282 if (snes->xl && snes->xu) {
283 PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, fnorm));
284 } else {
285 PetscCall(VecNorm(F, NORM_2, fnorm));
286 }
287 SNESCheckFunctionDomainError(snes, *fnorm);
288
289 /* take the minimum-normed candidate if it beats the combination by a factor of rtol or the combination has stagnated */
290 min_fnorm = jac->fnorms[0];
291 min_i = 0;
292 for (i = 0; i < jac->n; i++) {
293 if (jac->fnorms[i] < min_fnorm) {
294 min_fnorm = jac->fnorms[i];
295 min_i = i;
296 }
297 }
298
299 /* stagnation or divergence restart to the solution of the solver that failed the least */
300 if (PetscRealPart(total) < jac->stol || min_fnorm * jac->rtol < *fnorm) {
301 PetscCall(VecCopy(jac->Xes[min_i], X));
302 PetscCall(VecCopy(jac->Fes[min_i], F));
303 *fnorm = min_fnorm;
304 }
305 PetscFunctionReturn(PETSC_SUCCESS);
306 }
307
SNESSetUp_Composite(SNES snes)308 static PetscErrorCode SNESSetUp_Composite(SNES snes)
309 {
310 DM dm;
311 SNES_Composite *jac = (SNES_Composite *)snes->data;
312 SNES_CompositeLink next = jac->head;
313 PetscInt n = 0, i;
314 Vec F;
315
316 PetscFunctionBegin;
317 PetscCall(SNESGetDM(snes, &dm));
318
319 if (snes->ops->computevariablebounds) {
320 /* SNESVI only ever calls computevariablebounds once, so calling it once here is justified */
321 if (!snes->xl) PetscCall(VecDuplicate(snes->vec_sol, &snes->xl));
322 if (!snes->xu) PetscCall(VecDuplicate(snes->vec_sol, &snes->xu));
323 PetscUseTypeMethod(snes, computevariablebounds, snes->xl, snes->xu);
324 }
325
326 while (next) {
327 n++;
328 PetscCall(SNESSetDM(next->snes, dm));
329 PetscCall(SNESSetJacobian(next->snes, snes->jacobian, snes->jacobian_pre, NULL, NULL));
330 PetscCall(SNESSetApplicationContext(next->snes, snes->ctx));
331 if (snes->xl && snes->xu) {
332 if (snes->ops->computevariablebounds) {
333 PetscCall(SNESVISetComputeVariableBounds(next->snes, snes->ops->computevariablebounds));
334 } else {
335 PetscCall(SNESVISetVariableBounds(next->snes, snes->xl, snes->xu));
336 }
337 }
338
339 next = next->next;
340 }
341 jac->nsnes = n;
342 PetscCall(SNESGetFunction(snes, &F, NULL, NULL));
343 if (jac->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
344 PetscCall(VecDuplicateVecs(F, jac->nsnes, &jac->Xes));
345 PetscCall(PetscMalloc1(n, &jac->Fes));
346 PetscCall(PetscMalloc1(n, &jac->fnorms));
347 next = jac->head;
348 i = 0;
349 while (next) {
350 PetscCall(SNESGetFunction(next->snes, &F, NULL, NULL));
351 jac->Fes[i] = F;
352 PetscCall(PetscObjectReference((PetscObject)F));
353 next = next->next;
354 i++;
355 }
356 /* allocate the subspace direct solve area */
357 jac->nrhs = 1;
358 PetscCall(PetscBLASIntCast(jac->nsnes, &jac->lda));
359 PetscCall(PetscBLASIntCast(jac->nsnes, &jac->ldb));
360 PetscCall(PetscBLASIntCast(jac->nsnes, &jac->n));
361
362 PetscCall(PetscMalloc4(jac->n * jac->n, &jac->h, jac->n, &jac->beta, jac->n, &jac->s, jac->n, &jac->g));
363 jac->lwork = 12 * jac->n;
364 #if defined(PETSC_USE_COMPLEX)
365 PetscCall(PetscMalloc1(jac->lwork, &jac->rwork));
366 #endif
367 PetscCall(PetscMalloc1(jac->lwork, &jac->work));
368 }
369 PetscFunctionReturn(PETSC_SUCCESS);
370 }
371
SNESReset_Composite(SNES snes)372 static PetscErrorCode SNESReset_Composite(SNES snes)
373 {
374 SNES_Composite *jac = (SNES_Composite *)snes->data;
375 SNES_CompositeLink next = jac->head;
376
377 PetscFunctionBegin;
378 while (next) {
379 PetscCall(SNESReset(next->snes));
380 next = next->next;
381 }
382 PetscCall(VecDestroy(&jac->Xorig));
383 if (jac->Xes) PetscCall(VecDestroyVecs(jac->nsnes, &jac->Xes));
384 if (jac->Fes) PetscCall(VecDestroyVecs(jac->nsnes, &jac->Fes));
385 PetscCall(PetscFree(jac->fnorms));
386 PetscCall(PetscFree4(jac->h, jac->beta, jac->s, jac->g));
387 PetscCall(PetscFree(jac->work));
388 PetscCall(PetscFree(jac->rwork));
389 PetscFunctionReturn(PETSC_SUCCESS);
390 }
391
SNESDestroy_Composite(SNES snes)392 static PetscErrorCode SNESDestroy_Composite(SNES snes)
393 {
394 SNES_Composite *jac = (SNES_Composite *)snes->data;
395 SNES_CompositeLink next = jac->head, next_tmp;
396
397 PetscFunctionBegin;
398 PetscCall(SNESReset_Composite(snes));
399 while (next) {
400 PetscCall(SNESDestroy(&next->snes));
401 next_tmp = next;
402 next = next->next;
403 PetscCall(PetscFree(next_tmp));
404 }
405 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeSetType_C", NULL));
406 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeAddSNES_C", NULL));
407 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeGetSNES_C", NULL));
408 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeSetDamping_C", NULL));
409 PetscCall(PetscFree(snes->data));
410 PetscFunctionReturn(PETSC_SUCCESS);
411 }
412
SNESSetFromOptions_Composite(SNES snes,PetscOptionItems PetscOptionsObject)413 static PetscErrorCode SNESSetFromOptions_Composite(SNES snes, PetscOptionItems PetscOptionsObject)
414 {
415 SNES_Composite *jac = (SNES_Composite *)snes->data;
416 PetscInt nmax = 8, i;
417 SNES_CompositeLink next;
418 char *sneses[8];
419 PetscReal dmps[8];
420 PetscBool flg;
421
422 PetscFunctionBegin;
423 PetscOptionsHeadBegin(PetscOptionsObject, "Composite preconditioner options");
424 PetscCall(PetscOptionsEnum("-snes_composite_type", "Type of composition", "SNESCompositeSetType", SNESCompositeTypes, (PetscEnum)jac->type, (PetscEnum *)&jac->type, &flg));
425 if (flg) PetscCall(SNESCompositeSetType(snes, jac->type));
426 PetscCall(PetscOptionsStringArray("-snes_composite_sneses", "List of composite solvers", "SNESCompositeAddSNES", sneses, &nmax, &flg));
427 if (flg) {
428 for (i = 0; i < nmax; i++) {
429 PetscCall(SNESCompositeAddSNES(snes, sneses[i]));
430 PetscCall(PetscFree(sneses[i])); /* deallocate string sneses[i], which is allocated in PetscOptionsStringArray() */
431 }
432 }
433 PetscCall(PetscOptionsRealArray("-snes_composite_damping", "Damping of the additive composite solvers", "SNESCompositeSetDamping", dmps, &nmax, &flg));
434 if (flg) {
435 for (i = 0; i < nmax; i++) PetscCall(SNESCompositeSetDamping(snes, i, dmps[i]));
436 }
437 PetscCall(PetscOptionsReal("-snes_composite_stol", "Step tolerance for restart on the additive composite solvers", "", jac->stol, &jac->stol, NULL));
438 PetscCall(PetscOptionsReal("-snes_composite_rtol", "Residual tolerance for the additive composite solvers", "", jac->rtol, &jac->rtol, NULL));
439 PetscOptionsHeadEnd();
440
441 next = jac->head;
442 while (next) {
443 PetscCall(SNESSetFromOptions(next->snes));
444 next = next->next;
445 }
446 PetscFunctionReturn(PETSC_SUCCESS);
447 }
448
SNESView_Composite(SNES snes,PetscViewer viewer)449 static PetscErrorCode SNESView_Composite(SNES snes, PetscViewer viewer)
450 {
451 SNES_Composite *jac = (SNES_Composite *)snes->data;
452 SNES_CompositeLink next = jac->head;
453 PetscBool isascii;
454
455 PetscFunctionBegin;
456 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
457 if (isascii) {
458 PetscCall(PetscViewerASCIIPrintf(viewer, " type - %s\n", SNESCompositeTypes[jac->type]));
459 PetscCall(PetscViewerASCIIPrintf(viewer, " SNESes on composite preconditioner follow\n"));
460 PetscCall(PetscViewerASCIIPrintf(viewer, " ---------------------------------\n"));
461 }
462 if (isascii) PetscCall(PetscViewerASCIIPushTab(viewer));
463 while (next) {
464 PetscCall(SNESView(next->snes, viewer));
465 next = next->next;
466 }
467 if (isascii) {
468 PetscCall(PetscViewerASCIIPopTab(viewer));
469 PetscCall(PetscViewerASCIIPrintf(viewer, " ---------------------------------\n"));
470 }
471 PetscFunctionReturn(PETSC_SUCCESS);
472 }
473
SNESCompositeSetType_Composite(SNES snes,SNESCompositeType type)474 static PetscErrorCode SNESCompositeSetType_Composite(SNES snes, SNESCompositeType type)
475 {
476 SNES_Composite *jac = (SNES_Composite *)snes->data;
477
478 PetscFunctionBegin;
479 jac->type = type;
480 PetscFunctionReturn(PETSC_SUCCESS);
481 }
482
SNESCompositeAddSNES_Composite(SNES snes,SNESType type)483 static PetscErrorCode SNESCompositeAddSNES_Composite(SNES snes, SNESType type)
484 {
485 SNES_Composite *jac;
486 SNES_CompositeLink next, ilink;
487 PetscInt cnt = 0;
488 const char *prefix;
489 char newprefix[20];
490 DM dm;
491
492 PetscFunctionBegin;
493 PetscCall(PetscNew(&ilink));
494 ilink->next = NULL;
495 PetscCall(SNESCreate(PetscObjectComm((PetscObject)snes), &ilink->snes));
496 PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->snes, (PetscObject)snes, 1));
497 PetscCall(SNESGetDM(snes, &dm));
498 PetscCall(SNESSetDM(ilink->snes, dm));
499 PetscCall(SNESSetTolerances(ilink->snes, snes->abstol, snes->rtol, snes->stol, 1, snes->max_funcs));
500 PetscCall(PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)ilink->snes));
501 jac = (SNES_Composite *)snes->data;
502 next = jac->head;
503 if (!next) {
504 jac->head = ilink;
505 ilink->previous = NULL;
506 } else {
507 cnt++;
508 while (next->next) {
509 next = next->next;
510 cnt++;
511 }
512 next->next = ilink;
513 ilink->previous = next;
514 }
515 PetscCall(SNESGetOptionsPrefix(snes, &prefix));
516 PetscCall(SNESSetOptionsPrefix(ilink->snes, prefix));
517 PetscCall(PetscSNPrintf(newprefix, sizeof(newprefix), "sub_%" PetscInt_FMT "_", cnt));
518 PetscCall(SNESAppendOptionsPrefix(ilink->snes, newprefix));
519 PetscCall(SNESSetType(ilink->snes, type));
520 PetscCall(SNESSetNormSchedule(ilink->snes, SNES_NORM_FINAL_ONLY));
521
522 ilink->dmp = 1.0;
523 jac->nsnes++;
524 PetscFunctionReturn(PETSC_SUCCESS);
525 }
526
SNESCompositeGetSNES_Composite(SNES snes,PetscInt n,SNES * subsnes)527 static PetscErrorCode SNESCompositeGetSNES_Composite(SNES snes, PetscInt n, SNES *subsnes)
528 {
529 SNES_Composite *jac;
530 SNES_CompositeLink next;
531 PetscInt i;
532
533 PetscFunctionBegin;
534 jac = (SNES_Composite *)snes->data;
535 next = jac->head;
536 for (i = 0; i < n; i++) {
537 PetscCheck(next->next, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_INCOMP, "Not enough SNESes in composite preconditioner");
538 next = next->next;
539 }
540 *subsnes = next->snes;
541 PetscFunctionReturn(PETSC_SUCCESS);
542 }
543
544 /*@
545 SNESCompositeSetType - Sets the type of composite preconditioner.
546
547 Logically Collective
548
549 Input Parameters:
550 + snes - the preconditioner context
551 - type - `SNES_COMPOSITE_ADDITIVE` (default), `SNES_COMPOSITE_MULTIPLICATIVE`, or `SNES_COMPOSITE_ADDITIVEOPTIMAL`
552
553 Options Database Key:
554 . -snes_composite_type <type: one of multiplicative, additive, additiveoptimal> - Sets composite preconditioner type
555
556 Level: developer
557
558 .seealso: [](ch_snes), `SNES_COMPOSITE_ADDITIVE`, `SNES_COMPOSITE_MULTIPLICATIVE`, `SNESCompositeType`, `SNESCOMPOSITE`, `SNES_COMPOSITE_ADDITIVEOPTIMAL`,
559 `PCCompositeType`
560 @*/
SNESCompositeSetType(SNES snes,SNESCompositeType type)561 PetscErrorCode SNESCompositeSetType(SNES snes, SNESCompositeType type)
562 {
563 PetscFunctionBegin;
564 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
565 PetscValidLogicalCollectiveEnum(snes, type, 2);
566 PetscTryMethod(snes, "SNESCompositeSetType_C", (SNES, SNESCompositeType), (snes, type));
567 PetscFunctionReturn(PETSC_SUCCESS);
568 }
569
570 /*@
571 SNESCompositeAddSNES - Adds another `SNES` to the `SNESCOMPOSITE`
572
573 Collective
574
575 Input Parameters:
576 + snes - the `SNES` context of type `SNESCOMPOSITE`
577 - type - the `SNESType` of the new solver
578
579 Level: developer
580
581 .seealso: [](ch_snes), `SNES`, `SNESCOMPOSITE`, `SNESCompositeGetSNES()`
582 @*/
SNESCompositeAddSNES(SNES snes,SNESType type)583 PetscErrorCode SNESCompositeAddSNES(SNES snes, SNESType type)
584 {
585 PetscFunctionBegin;
586 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
587 PetscTryMethod(snes, "SNESCompositeAddSNES_C", (SNES, SNESType), (snes, type));
588 PetscFunctionReturn(PETSC_SUCCESS);
589 }
590
591 /*@
592 SNESCompositeGetSNES - Gets one of the `SNES` objects in the `SNES` of `SNESType` `SNESCOMPOSITE`
593
594 Not Collective
595
596 Input Parameters:
597 + snes - the `SNES` context
598 - n - the number of the composed `SNES` requested
599
600 Output Parameter:
601 . subsnes - the `SNES` requested
602
603 Level: developer
604
605 .seealso: [](ch_snes), `SNES`, `SNESCOMPOSITE`, `SNESCompositeAddSNES()`, `SNESCompositeGetNumber()`
606 @*/
SNESCompositeGetSNES(SNES snes,PetscInt n,SNES * subsnes)607 PetscErrorCode SNESCompositeGetSNES(SNES snes, PetscInt n, SNES *subsnes)
608 {
609 PetscFunctionBegin;
610 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
611 PetscAssertPointer(subsnes, 3);
612 PetscUseMethod(snes, "SNESCompositeGetSNES_C", (SNES, PetscInt, SNES *), (snes, n, subsnes));
613 PetscFunctionReturn(PETSC_SUCCESS);
614 }
615
616 /*@
617 SNESCompositeGetNumber - Get the number of subsolvers in the `SNESCOMPOSITE`
618
619 Logically Collective
620
621 Input Parameter:
622 . snes - the `SNES` context
623
624 Output Parameter:
625 . n - the number of subsolvers
626
627 Level: developer
628
629 .seealso: [](ch_snes), `SNES`, `SNESCOMPOSITE`, `SNESCompositeAddSNES()`, `SNESCompositeGetSNES()`
630 @*/
SNESCompositeGetNumber(SNES snes,PetscInt * n)631 PetscErrorCode SNESCompositeGetNumber(SNES snes, PetscInt *n)
632 {
633 SNES_Composite *jac;
634 SNES_CompositeLink next;
635
636 PetscFunctionBegin;
637 jac = (SNES_Composite *)snes->data;
638 next = jac->head;
639
640 *n = 0;
641 while (next) {
642 *n = *n + 1;
643 next = next->next;
644 }
645 PetscFunctionReturn(PETSC_SUCCESS);
646 }
647
SNESCompositeSetDamping_Composite(SNES snes,PetscInt n,PetscReal dmp)648 static PetscErrorCode SNESCompositeSetDamping_Composite(SNES snes, PetscInt n, PetscReal dmp)
649 {
650 SNES_Composite *jac;
651 SNES_CompositeLink next;
652 PetscInt i;
653
654 PetscFunctionBegin;
655 jac = (SNES_Composite *)snes->data;
656 next = jac->head;
657 for (i = 0; i < n; i++) {
658 PetscCheck(next->next, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_INCOMP, "Not enough SNESes in composite preconditioner");
659 next = next->next;
660 }
661 next->dmp = dmp;
662 PetscFunctionReturn(PETSC_SUCCESS);
663 }
664
665 /*@
666 SNESCompositeSetDamping - Sets the damping of a subsolver when using `SNES_COMPOSITE_ADDITIVE` with a `SNES` of `SNESType` `SNESCOMPOSITE`
667
668 Not Collective
669
670 Input Parameters:
671 + snes - the `SNES` context
672 . n - the number of the sub-`SNES` object requested
673 - dmp - the damping
674
675 Level: intermediate
676
677 .seealso: [](ch_snes), `SNES`, `SNESCOMPOSITE`, `SNESCompositeAddSNES()`, `SNESCompositeGetSNES()`,
678 `SNES_COMPOSITE_ADDITIVE`, `SNES_COMPOSITE_MULTIPLICATIVE`, `SNESCompositeType`, `SNESCompositeSetType()`
679 @*/
SNESCompositeSetDamping(SNES snes,PetscInt n,PetscReal dmp)680 PetscErrorCode SNESCompositeSetDamping(SNES snes, PetscInt n, PetscReal dmp)
681 {
682 PetscFunctionBegin;
683 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
684 PetscUseMethod(snes, "SNESCompositeSetDamping_C", (SNES, PetscInt, PetscReal), (snes, n, dmp));
685 PetscFunctionReturn(PETSC_SUCCESS);
686 }
687
SNESSolve_Composite(SNES snes)688 static PetscErrorCode SNESSolve_Composite(SNES snes)
689 {
690 Vec F, X, B, Y;
691 PetscInt i;
692 PetscReal fnorm = 0.0, xnorm = 0.0, snorm = 0.0;
693 SNESNormSchedule normtype;
694 SNES_Composite *comp = (SNES_Composite *)snes->data;
695
696 PetscFunctionBegin;
697 X = snes->vec_sol;
698 F = snes->vec_func;
699 B = snes->vec_rhs;
700 Y = snes->vec_sol_update;
701
702 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
703 snes->iter = 0;
704 snes->norm = 0.;
705 comp->innerFailures = 0;
706 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
707 snes->reason = SNES_CONVERGED_ITERATING;
708 PetscCall(SNESGetNormSchedule(snes, &normtype));
709 if (normtype == SNES_NORM_ALWAYS || normtype == SNES_NORM_INITIAL_ONLY || normtype == SNES_NORM_INITIAL_FINAL_ONLY) {
710 if (!snes->vec_func_init_set) PetscCall(SNESComputeFunction(snes, X, F));
711 else snes->vec_func_init_set = PETSC_FALSE;
712
713 if (snes->xl && snes->xu) {
714 PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm));
715 } else {
716 PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- ||F|| */
717 }
718 SNESCheckFunctionDomainError(snes, fnorm);
719 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
720 snes->iter = 0;
721 snes->norm = fnorm;
722 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
723 PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
724
725 /* test convergence */
726 PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm));
727 PetscCall(SNESMonitor(snes, 0, snes->norm));
728 if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
729 } else {
730 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
731 PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
732 PetscCall(SNESMonitor(snes, 0, snes->norm));
733 }
734
735 for (i = 0; i < snes->max_its; i++) {
736 /* Call general purpose update function */
737 PetscTryTypeMethod(snes, update, snes->iter);
738
739 /* Copy the state before modification by application of the composite solver;
740 we will subtract the new state after application */
741 PetscCall(VecCopy(X, Y));
742
743 if (comp->type == SNES_COMPOSITE_ADDITIVE) {
744 PetscCall(SNESCompositeApply_Additive(snes, X, B, F, &fnorm));
745 } else if (comp->type == SNES_COMPOSITE_MULTIPLICATIVE) {
746 PetscCall(SNESCompositeApply_Multiplicative(snes, X, B, F, &fnorm));
747 } else if (comp->type == SNES_COMPOSITE_ADDITIVEOPTIMAL) {
748 PetscCall(SNESCompositeApply_AdditiveOptimal(snes, X, B, F, &fnorm));
749 } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Unsupported SNESComposite type");
750 if (snes->reason < 0) break;
751
752 /* Compute the solution update for convergence testing */
753 PetscCall(VecAYPX(Y, -1.0, X));
754
755 if ((i == snes->max_its - 1) && (normtype == SNES_NORM_INITIAL_FINAL_ONLY || normtype == SNES_NORM_FINAL_ONLY)) {
756 PetscCall(SNESComputeFunction(snes, X, F));
757
758 if (snes->xl && snes->xu) {
759 PetscCall(VecNormBegin(X, NORM_2, &xnorm));
760 PetscCall(VecNormBegin(Y, NORM_2, &snorm));
761 PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm));
762 PetscCall(VecNormEnd(X, NORM_2, &xnorm));
763 PetscCall(VecNormEnd(Y, NORM_2, &snorm));
764 } else {
765 PetscCall(VecNormBegin(F, NORM_2, &fnorm));
766 PetscCall(VecNormBegin(X, NORM_2, &xnorm));
767 PetscCall(VecNormBegin(Y, NORM_2, &snorm));
768
769 PetscCall(VecNormEnd(F, NORM_2, &fnorm));
770 PetscCall(VecNormEnd(X, NORM_2, &xnorm));
771 PetscCall(VecNormEnd(Y, NORM_2, &snorm));
772 }
773 SNESCheckFunctionDomainError(snes, fnorm);
774 } else if (normtype == SNES_NORM_ALWAYS) {
775 PetscCall(VecNormBegin(X, NORM_2, &xnorm));
776 PetscCall(VecNormBegin(Y, NORM_2, &snorm));
777 PetscCall(VecNormEnd(X, NORM_2, &xnorm));
778 PetscCall(VecNormEnd(Y, NORM_2, &snorm));
779 }
780 /* Monitor convergence */
781 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
782 snes->iter = i + 1;
783 snes->norm = fnorm;
784 snes->xnorm = xnorm;
785 snes->ynorm = snorm;
786 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
787 PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
788 /* Test for convergence */
789 PetscCall(SNESConverged(snes, snes->iter, xnorm, snorm, fnorm));
790 PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
791 if (snes->reason) break;
792 }
793 PetscFunctionReturn(PETSC_SUCCESS);
794 }
795
796 /*MC
797 SNESCOMPOSITE - Builds a nonlinear solver/preconditioner by composing together several `SNES` nonlinear solvers {cite}`bruneknepleysmithtu15`
798
799 Options Database Keys:
800 + -snes_composite_type <type: one of multiplicative, additive, additiveoptimal> - Sets composite preconditioner type
801 - -snes_composite_sneses - <snes0,snes1,...> list of `SNES` to compose
802
803 Level: intermediate
804
805 .seealso: [](ch_snes), `SNES`, `SNESCOMPOSITE`, `SNESCompositeAddSNES()`, `SNESCompositeGetSNES()`,
806 `SNES_COMPOSITE_ADDITIVE`, `SNES_COMPOSITE_ADDITIVEOPTIMAL`, `SNES_COMPOSITE_MULTIPLICATIVE`, `SNESCompositeType`,
807 `SNESCompositeSetType()`, `SNESCompositeSetDamping()`, `SNESCompositeGetNumber()`, `PCCompositeType`
808 M*/
809
SNESCreate_Composite(SNES snes)810 PETSC_EXTERN PetscErrorCode SNESCreate_Composite(SNES snes)
811 {
812 SNES_Composite *jac;
813
814 PetscFunctionBegin;
815 PetscCall(PetscNew(&jac));
816
817 snes->ops->solve = SNESSolve_Composite;
818 snes->ops->setup = SNESSetUp_Composite;
819 snes->ops->reset = SNESReset_Composite;
820 snes->ops->destroy = SNESDestroy_Composite;
821 snes->ops->setfromoptions = SNESSetFromOptions_Composite;
822 snes->ops->view = SNESView_Composite;
823
824 snes->usesksp = PETSC_FALSE;
825
826 snes->alwayscomputesfinalresidual = PETSC_FALSE;
827
828 PetscCall(SNESParametersInitialize(snes));
829
830 snes->data = (void *)jac;
831 jac->type = SNES_COMPOSITE_ADDITIVEOPTIMAL;
832 jac->Fes = NULL;
833 jac->Xes = NULL;
834 jac->fnorms = NULL;
835 jac->nsnes = 0;
836 jac->head = NULL;
837 jac->stol = 0.1;
838 jac->rtol = 1.1;
839
840 jac->h = NULL;
841 jac->s = NULL;
842 jac->beta = NULL;
843 jac->work = NULL;
844 jac->rwork = NULL;
845
846 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeSetType_C", SNESCompositeSetType_Composite));
847 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeAddSNES_C", SNESCompositeAddSNES_Composite));
848 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeGetSNES_C", SNESCompositeGetSNES_Composite));
849 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESCompositeSetDamping_C", SNESCompositeSetDamping_Composite));
850 PetscFunctionReturn(PETSC_SUCCESS);
851 }
852