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