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