1 /* 2 3 */ 4 #include "src/ksp/pc/pcimpl.h" /*I "petscpc.h" I*/ 5 6 typedef struct _PC_FieldSplitLink *PC_FieldSplitLink; 7 struct _PC_FieldSplitLink { 8 KSP ksp; 9 Vec x,y; 10 PetscInt nfields; 11 PetscInt *fields; 12 VecScatter sctx; 13 PC_FieldSplitLink next; 14 }; 15 16 typedef struct { 17 PetscTruth defaultsplit; 18 PetscInt bs; 19 PetscInt nsplits; 20 Vec *x,*y; 21 Mat *pmat; 22 IS *is; 23 PC_FieldSplitLink head; 24 } PC_FieldSplit; 25 26 #undef __FUNCT__ 27 #define __FUNCT__ "PCView_FieldSplit" 28 static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer) 29 { 30 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 31 PetscErrorCode ierr; 32 PetscTruth iascii; 33 PetscInt i,j; 34 PC_FieldSplitLink link = jac->head; 35 36 PetscFunctionBegin; 37 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 38 if (iascii) { 39 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit: total splits = %D",jac->nsplits);CHKERRQ(ierr); 40 ierr = PetscViewerASCIIPrintf(viewer," Solver info for each split is in the following KSP objects:\n");CHKERRQ(ierr); 41 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 42 for (i=0; i<jac->nsplits; i++) { 43 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr); 44 for (j=0; j<link->nfields; j++) { 45 ierr = PetscViewerASCIIPrintf(viewer,"%D \n",link->fields[j]);CHKERRQ(ierr); 46 } 47 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 48 ierr = KSPView(link->ksp,viewer);CHKERRQ(ierr); 49 link = link->next; 50 } 51 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 52 } else { 53 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCFieldSplit",((PetscObject)viewer)->type_name); 54 } 55 PetscFunctionReturn(0); 56 } 57 58 #undef __FUNCT__ 59 #define __FUNCT__ "PCFieldSplitSetDefaults" 60 static PetscErrorCode PCFieldSplitSetDefaults(PC pc) 61 { 62 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 63 PetscErrorCode ierr; 64 PC_FieldSplitLink link = jac->head; 65 PetscInt i; 66 67 PetscFunctionBegin; 68 /* user has not split fields so use default */ 69 if (!link) { 70 ierr = PetscLogInfo(pc,"PCFieldSplitSetDefaults: Using default splitting of fields");CHKERRQ(ierr); 71 if (jac->bs <= 0) { 72 ierr = MatGetBlockSize(pc->pmat,&jac->bs);CHKERRQ(ierr); 73 } 74 for (i=0; i<jac->bs; i++) { 75 ierr = PCFieldSplitSetFields(pc,1,&i);CHKERRQ(ierr); 76 } 77 link = jac->head; 78 jac->defaultsplit = PETSC_TRUE; 79 } 80 PetscFunctionReturn(0); 81 } 82 83 84 #undef __FUNCT__ 85 #define __FUNCT__ "PCSetUp_FieldSplit" 86 static PetscErrorCode PCSetUp_FieldSplit(PC pc) 87 { 88 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 89 PetscErrorCode ierr; 90 PC_FieldSplitLink link; 91 PetscInt i,nsplit; 92 MatStructure flag = pc->flag; 93 94 PetscFunctionBegin; 95 ierr = PCFieldSplitSetDefaults(pc);CHKERRQ(ierr); 96 nsplit = jac->nsplits; 97 link = jac->head; 98 99 /* get the matrices for each split */ 100 if (!jac->is) { 101 PetscInt rstart,rend,nslots,bs; 102 103 ierr = MatGetBlockSize(pc->pmat,&bs);CHKERRQ(ierr); 104 ierr = MatGetOwnershipRange(pc->pmat,&rstart,&rend);CHKERRQ(ierr); 105 nslots = (rend - rstart)/bs; 106 ierr = PetscMalloc(nsplit*sizeof(IS),&jac->is);CHKERRQ(ierr); 107 for (i=0; i<nsplit; i++) { 108 if (jac->defaultsplit) { 109 ierr = ISCreateStride(pc->comm,nslots,rstart+i,nsplit,&jac->is[i]);CHKERRQ(ierr); 110 } else { 111 PetscInt *ii,j,k,nfields = link->nfields,*fields = link->fields; 112 PetscTruth sorted; 113 ierr = PetscMalloc(link->nfields*nslots*sizeof(PetscInt),&ii);CHKERRQ(ierr); 114 for (j=0; j<nslots; j++) { 115 for (k=0; k<nfields; k++) { 116 ii[nfields*j + k] = rstart + bs*j + fields[k]; 117 } 118 } 119 ierr = ISCreateGeneral(pc->comm,nslots*nfields,ii,&jac->is[i]);CHKERRQ(ierr); 120 ierr = ISSorted(jac->is[i],&sorted);CHKERRQ(ierr); 121 if (!sorted) SETERRQ(PETSC_ERR_USER,"Fields must be sorted when creating split"); 122 ierr = PetscFree(ii);CHKERRQ(ierr); 123 link = link->next; 124 } 125 } 126 } 127 128 if (!jac->pmat) { 129 ierr = MatGetSubMatrices(pc->pmat,nsplit,jac->is,jac->is,MAT_INITIAL_MATRIX,&jac->pmat);CHKERRQ(ierr); 130 } else { 131 ierr = MatGetSubMatrices(pc->pmat,nsplit,jac->is,jac->is,MAT_REUSE_MATRIX,&jac->pmat);CHKERRQ(ierr); 132 } 133 134 /* set up the individual PCs */ 135 i = 0; 136 link = jac->head; 137 while (link) { 138 ierr = KSPSetOperators(link->ksp,jac->pmat[i],jac->pmat[i],flag);CHKERRQ(ierr); 139 ierr = KSPSetFromOptions(link->ksp);CHKERRQ(ierr); 140 ierr = KSPSetUp(link->ksp);CHKERRQ(ierr); 141 i++; 142 link = link->next; 143 } 144 145 /* create work vectors for each split */ 146 if (!jac->x) { 147 ierr = PetscMalloc2(nsplit,Vec,&jac->x,nsplit,Vec,&jac->y);CHKERRQ(ierr); 148 link = jac->head; 149 for (i=0; i<nsplit; i++) { 150 Mat A; 151 ierr = KSPGetOperators(link->ksp,PETSC_NULL,&A,PETSC_NULL);CHKERRQ(ierr); 152 ierr = MatGetVecs(A,&link->x,&link->y);CHKERRQ(ierr); 153 jac->x[i] = link->x; 154 jac->y[i] = link->y; 155 link = link->next; 156 } 157 if (!jac->defaultsplit) { 158 Vec xtmp; 159 160 link = jac->head; 161 ierr = MatGetVecs(pc->pmat,&xtmp,PETSC_NULL);CHKERRQ(ierr); 162 for (i=0; i<nsplit; i++) { 163 ierr = VecScatterCreate(xtmp,jac->is[i],jac->x[i],PETSC_NULL,&link->sctx);CHKERRQ(ierr); 164 link = link->next; 165 } 166 ierr = VecDestroy(xtmp);CHKERRQ(ierr); 167 } 168 } 169 PetscFunctionReturn(0); 170 } 171 172 #undef __FUNCT__ 173 #define __FUNCT__ "PCApply_FieldSplit" 174 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y) 175 { 176 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 177 PetscErrorCode ierr; 178 PC_FieldSplitLink link = jac->head; 179 180 PetscFunctionBegin; 181 if (jac->defaultsplit) { 182 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 183 while (link) { 184 ierr = KSPSolve(link->ksp,link->x,link->y);CHKERRQ(ierr); 185 link = link->next; 186 } 187 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 188 } else { 189 PetscScalar zero = 0.0; 190 PetscInt i = 0; 191 192 ierr = VecSet(&zero,y);CHKERRQ(ierr); 193 while (link) { 194 ierr = VecScatterBegin(x,link->x,INSERT_VALUES,SCATTER_FORWARD,link->sctx);CHKERRQ(ierr); 195 ierr = VecScatterEnd(x,link->x,INSERT_VALUES,SCATTER_FORWARD,link->sctx);CHKERRQ(ierr); 196 ierr = KSPSolve(link->ksp,link->x,link->y);CHKERRQ(ierr); 197 ierr = VecScatterBegin(link->y,y,ADD_VALUES,SCATTER_REVERSE,link->sctx);CHKERRQ(ierr); 198 ierr = VecScatterEnd(y,link->y,ADD_VALUES,SCATTER_REVERSE,link->sctx);CHKERRQ(ierr); 199 200 link = link->next; 201 i++; 202 } 203 } 204 PetscFunctionReturn(0); 205 } 206 207 #undef __FUNCT__ 208 #define __FUNCT__ "PCDestroy_FieldSplit" 209 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 210 { 211 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 212 PetscErrorCode ierr; 213 PC_FieldSplitLink link = jac->head,next; 214 215 PetscFunctionBegin; 216 while (link) { 217 ierr = KSPDestroy(link->ksp);CHKERRQ(ierr); 218 if (link->x) {ierr = VecDestroy(link->x);CHKERRQ(ierr);} 219 if (link->y) {ierr = VecDestroy(link->y);CHKERRQ(ierr);} 220 if (link->sctx) {ierr = VecScatterDestroy(link->sctx);CHKERRQ(ierr);} 221 next = link->next; 222 ierr = PetscFree2(link,link->fields);CHKERRQ(ierr); 223 link = next; 224 } 225 if (jac->x) {ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr);} 226 if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);} 227 if (jac->is) { 228 PetscInt i; 229 for (i=0; i<jac->nsplits; i++) {ierr = ISDestroy(jac->is[i]);CHKERRQ(ierr);} 230 ierr = PetscFree(jac->is);CHKERRQ(ierr); 231 } 232 ierr = PetscFree(jac);CHKERRQ(ierr); 233 PetscFunctionReturn(0); 234 } 235 236 #undef __FUNCT__ 237 #define __FUNCT__ "PCSetFromOptions_FieldSplit" 238 static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc) 239 /* This does not call KSPSetFromOptions() on the subksp's, see PCSetFromOptionsBJacobi/ASM() */ 240 { 241 PetscErrorCode ierr; 242 PetscInt i = 0,nfields,fields[12]; 243 PetscTruth flg; 244 char optionname[128]; 245 246 PetscFunctionBegin; 247 ierr = PetscOptionsHead("FieldSplit options");CHKERRQ(ierr); 248 while (PETSC_TRUE) { 249 sprintf(optionname,"-pc_fieldsplit_%d_fields",i++); 250 nfields = 12; 251 ierr = PetscOptionsIntArray(optionname,"Fields in this split","PCFieldSplitSetFields",fields,&nfields,&flg);CHKERRQ(ierr); 252 if (!flg) break; 253 if (!nfields) SETERRQ(PETSC_ERR_USER,"Cannot list zero fields"); 254 ierr = PCFieldSplitSetFields(pc,nfields,fields);CHKERRQ(ierr); 255 } 256 ierr = PetscOptionsTail();CHKERRQ(ierr); 257 PetscFunctionReturn(0); 258 } 259 260 /*------------------------------------------------------------------------------------*/ 261 262 EXTERN_C_BEGIN 263 #undef __FUNCT__ 264 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 265 PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,PetscInt n,PetscInt *fields) 266 { 267 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 268 PetscErrorCode ierr; 269 PC_FieldSplitLink link,next = jac->head; 270 char prefix[128]; 271 272 PetscFunctionBegin; 273 if (n <= 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Negative number of fields requested"); 274 ierr = PetscMalloc2(1,struct _PC_FieldSplitLink,&link,n,PetscInt,&link->fields);CHKERRQ(ierr); 275 ierr = PetscMemcpy(link->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 276 link->nfields = n; 277 link->next = PETSC_NULL; 278 ierr = KSPCreate(pc->comm,&link->ksp);CHKERRQ(ierr); 279 ierr = KSPSetType(link->ksp,KSPPREONLY);CHKERRQ(ierr); 280 281 if (pc->prefix) { 282 sprintf(prefix,"%s_fieldsplit_%d_",pc->prefix,jac->nsplits); 283 } else { 284 sprintf(prefix,"fieldsplit_%d_",jac->nsplits); 285 } 286 ierr = KSPSetOptionsPrefix(link->ksp,prefix);CHKERRQ(ierr); 287 288 if (!next) { 289 jac->head = link; 290 } else { 291 while (next->next) { 292 next = next->next; 293 } 294 next->next = link; 295 } 296 jac->nsplits++; 297 PetscFunctionReturn(0); 298 } 299 EXTERN_C_END 300 301 302 EXTERN_C_BEGIN 303 #undef __FUNCT__ 304 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 305 PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 306 { 307 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 308 PetscErrorCode ierr; 309 PetscInt cnt = 0; 310 PC_FieldSplitLink link; 311 312 PetscFunctionBegin; 313 ierr = PetscMalloc(jac->nsplits*sizeof(KSP*),subksp);CHKERRQ(ierr); 314 while (link) { 315 (*subksp)[cnt++] = link->ksp; 316 link = link->next; 317 } 318 if (cnt != jac->nsplits) SETERRQ2(PETSC_ERR_PLIB,"Corrupt PCFIELDSPLIT object: number splits in linked list %D in object %D",cnt,jac->nsplits); 319 *n = jac->nsplits; 320 PetscFunctionReturn(0); 321 } 322 EXTERN_C_END 323 324 #undef __FUNCT__ 325 #define __FUNCT__ "PCFieldSplitSetFields" 326 /*@ 327 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 328 329 Collective on PC 330 331 Input Parameters: 332 + pc - the preconditioner context 333 . n - the number of fields in this split 334 . fields - the fields in this split 335 336 Level: intermediate 337 338 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT 339 340 @*/ 341 PetscErrorCode PCFieldSplitSetFields(PC pc,PetscInt n, PetscInt *fields) 342 { 343 PetscErrorCode ierr,(*f)(PC,PetscInt,PetscInt *); 344 345 PetscFunctionBegin; 346 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 347 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSetFields_C",(void (**)(void))&f);CHKERRQ(ierr); 348 if (f) { 349 ierr = (*f)(pc,n,fields);CHKERRQ(ierr); 350 } 351 PetscFunctionReturn(0); 352 } 353 354 #undef __FUNCT__ 355 #define __FUNCT__ "PCFieldSplitGetSubKSP" 356 /*@C 357 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 358 359 Collective on KSP 360 361 Input Parameter: 362 . pc - the preconditioner context 363 364 Output Parameters: 365 + n - the number of split 366 - pc - the array of KSP contexts 367 368 Note: 369 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed 370 371 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 372 373 Level: advanced 374 375 .seealso: PCFIELDSPLIT 376 @*/ 377 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 378 { 379 PetscErrorCode ierr,(*f)(PC,PetscInt*,KSP **); 380 381 PetscFunctionBegin; 382 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 383 PetscValidIntPointer(n,2); 384 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",(void (**)(void))&f);CHKERRQ(ierr); 385 if (f) { 386 ierr = (*f)(pc,n,subksp);CHKERRQ(ierr); 387 } else { 388 SETERRQ(PETSC_ERR_ARG_WRONG,"Cannot get subksp for this type of PC"); 389 } 390 PetscFunctionReturn(0); 391 } 392 393 /* -------------------------------------------------------------------------------------*/ 394 /*MC 395 PCFieldSplit - Preconditioner created by combining seperate preconditioners for individual 396 fields or groups of fields 397 398 399 To set options on the solvers for each block append -sub_ to all the PC 400 options database keys. For example, -sub_pc_type ilu -sub_pc_ilu_levels 1 401 402 To set the options on the solvers seperate for each block call PCFieldSplitGetSubKSP() 403 and set the options directly on the resulting KSP object 404 405 Level: intermediate 406 407 Concepts: physics based preconditioners 408 409 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 410 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields() 411 M*/ 412 413 EXTERN_C_BEGIN 414 #undef __FUNCT__ 415 #define __FUNCT__ "PCCreate_FieldSplit" 416 PetscErrorCode PCCreate_FieldSplit(PC pc) 417 { 418 PetscErrorCode ierr; 419 PC_FieldSplit *jac; 420 421 PetscFunctionBegin; 422 ierr = PetscNew(PC_FieldSplit,&jac);CHKERRQ(ierr); 423 PetscLogObjectMemory(pc,sizeof(PC_FieldSplit)); 424 ierr = PetscMemzero(jac,sizeof(PC_FieldSplit));CHKERRQ(ierr); 425 jac->bs = -1; 426 jac->nsplits = 0; 427 pc->data = (void*)jac; 428 429 pc->ops->apply = PCApply_FieldSplit; 430 pc->ops->setup = PCSetUp_FieldSplit; 431 pc->ops->destroy = PCDestroy_FieldSplit; 432 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 433 pc->ops->view = PCView_FieldSplit; 434 pc->ops->applyrichardson = 0; 435 436 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit", 437 PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 438 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","PCFieldSplitSetFields_FieldSplit", 439 PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 440 PetscFunctionReturn(0); 441 } 442 EXTERN_C_END 443 444 445