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