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