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