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 PetscInt bs; 17 PetscInt nsplits; 18 PC_FieldSplitLink head; 19 Vec *x,*y; 20 } PC_FieldSplit; 21 22 #undef __FUNCT__ 23 #define __FUNCT__ "PCView_FieldSplit" 24 static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer) 25 { 26 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 27 PetscErrorCode ierr; 28 PetscTruth iascii; 29 PetscInt i,j; 30 PC_FieldSplitLink link = jac->head; 31 32 PetscFunctionBegin; 33 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 34 if (iascii) { 35 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit: total splits = %D",jac->nsplits);CHKERRQ(ierr); 36 ierr = PetscViewerASCIIPrintf(viewer," Solver info for each split is in the following PC objects:\n");CHKERRQ(ierr); 37 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 38 for (i=0; i<jac->nsplits; i++) { 39 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr); 40 for (j=0; j<link->nfields; j++) { 41 ierr = PetscViewerASCIIPrintf(viewer,"%D \n",link->fields[j]);CHKERRQ(ierr); 42 } 43 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 44 ierr = PCView(link->pc,viewer);CHKERRQ(ierr); 45 link = link->next; 46 } 47 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 48 } else { 49 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCFieldSplit",((PetscObject)viewer)->type_name); 50 } 51 PetscFunctionReturn(0); 52 } 53 54 #undef __FUNCT__ 55 #define __FUNCT__ "PCSetUp_FieldSplit" 56 static PetscErrorCode PCSetUp_FieldSplit(PC pc) 57 { 58 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 59 PetscErrorCode ierr; 60 PC_FieldSplitLink link = jac->head; 61 PetscInt i; 62 63 PetscFunctionBegin; 64 if (jac->bs <= 0) { 65 ierr = MatGetBlockSize(pc->pmat,&jac->bs);CHKERRQ(ierr); 66 } 67 if (!link) { /* user has not split fields so use default */ 68 for (i=0; i<jac->bs; i++) { 69 ierr = PCFieldSplitSetFields(pc,1,&i);CHKERRQ(ierr); 70 } 71 link = jac->head; 72 } 73 while (link) { 74 ierr = PCSetUp(link->pc);CHKERRQ(ierr); 75 link = link->next; 76 } 77 PetscFunctionReturn(0); 78 } 79 80 #undef __FUNCT__ 81 #define __FUNCT__ "PCApply_FieldSplit" 82 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y) 83 { 84 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 85 PetscErrorCode ierr; 86 PC_FieldSplitLink link = jac->head; 87 88 PetscFunctionBegin; 89 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 90 while (link) { 91 ierr = PCApply(link->pc,link->x,link->y);CHKERRQ(ierr); 92 link = link->next; 93 } 94 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 95 PetscFunctionReturn(0); 96 } 97 98 #undef __FUNCT__ 99 #define __FUNCT__ "PCDestroy_FieldSplit" 100 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 101 { 102 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 103 PetscErrorCode ierr; 104 PC_FieldSplitLink link = jac->head; 105 106 PetscFunctionBegin; 107 while (link) { 108 ierr = PCDestroy(link->pc);CHKERRQ(ierr); 109 ierr = PetscFree2(link,link->fields);CHKERRQ(ierr); 110 } 111 if (jac->x) {ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr);} 112 ierr = PetscFree(jac);CHKERRQ(ierr); 113 PetscFunctionReturn(0); 114 } 115 116 #undef __FUNCT__ 117 #define __FUNCT__ "PCSetFromOptions_FieldSplit" 118 static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc) 119 { 120 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 121 PetscErrorCode ierr; 122 PC_FieldSplitLink link = jac->head; 123 124 PetscFunctionBegin; 125 if (!link) { /* user never set fields, so set them ourselves */ 126 link = jac->head; 127 } 128 while (link) { 129 ierr = PCApply(link->pc,link->x,link->y);CHKERRQ(ierr); 130 link = link->next; 131 } 132 PetscFunctionReturn(0); 133 } 134 135 /*------------------------------------------------------------------------------------*/ 136 137 EXTERN_C_BEGIN 138 #undef __FUNCT__ 139 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 140 PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,PetscInt n,PetscInt *fields) 141 { 142 PC_FieldSplit *jac; 143 PetscErrorCode ierr; 144 PetscInt *myfields; 145 PC_FieldSplitLink link,next = jac->head; 146 147 PetscFunctionBegin; 148 if (n <= 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Negative number of fields requested"); 149 ierr = PetscMalloc2(1,struct _PC_FieldSplitLink,&link,n,PetscInt,&myfields);CHKERRQ(ierr); 150 ierr = PetscMemcpy(myfields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 151 link->fields = myfields; 152 link->nfields = n; 153 154 if (!next) { 155 jac->head = link; 156 } else { 157 while (next->next) { 158 next = next->next; 159 } 160 next->next = link; 161 } 162 jac->nsplits++; 163 PetscFunctionReturn(0); 164 } 165 EXTERN_C_END 166 167 168 EXTERN_C_BEGIN 169 #undef __FUNCT__ 170 #define __FUNCT__ "PCFieldSplitGetSubPC_FieldSplit" 171 PetscErrorCode PCFieldSplitGetSubPC_FieldSplit(PC pc,PetscInt *n,PC **subpc) 172 { 173 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 174 PetscErrorCode ierr; 175 PetscInt cnt = 0; 176 PC_FieldSplitLink link; 177 178 PetscFunctionBegin; 179 ierr = PetscMalloc(jac->nsplits*sizeof(PC*),subpc);CHKERRQ(ierr); 180 while (link) { 181 (*subpc)[cnt++] = link->pc; 182 link = link->next; 183 } 184 if (cnt != jac->nsplits) SETERRQ2(PETSC_ERR_PLIB,"Corrupt PCFIELDSPLIT object: number splits in linked list %D in object %D",cnt,jac->nsplits); 185 *n = jac->nsplits; 186 PetscFunctionReturn(0); 187 } 188 EXTERN_C_END 189 190 #undef __FUNCT__ 191 #define __FUNCT__ "PCFieldSplitSetFields" 192 /*@ 193 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 194 195 Collective on PC 196 197 Input Parameters: 198 + pc - the preconditioner context 199 . n - the number of fields in this split 200 . fields - the fields in this split 201 202 Level: intermediate 203 204 .seealso: PCFieldSplitGetSubPC(), PCFIELDSPLIT 205 206 @*/ 207 PetscErrorCode PCFieldSplitSetFields(PC pc,PetscInt n, PetscInt *fields) 208 { 209 PetscErrorCode ierr,(*f)(PC,PetscInt,PetscInt *); 210 211 PetscFunctionBegin; 212 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 213 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSetFields_C",(void (**)(void))&f);CHKERRQ(ierr); 214 if (f) { 215 ierr = (*f)(pc,n,fields);CHKERRQ(ierr); 216 } 217 PetscFunctionReturn(0); 218 } 219 220 #undef __FUNCT__ 221 #define __FUNCT__ "PCFieldSplitGetSubPC" 222 /*@C 223 PCFieldSplitGetSubPC - Gets the PC contexts for all splits 224 225 Collective on PC 226 227 Input Parameter: 228 . pc - the preconditioner context 229 230 Output Parameters: 231 + n - the number of split 232 - pc - the array of PC contexts 233 234 Note: 235 After PCFieldSplitGetSubPC() the array of PCs IS to be freed 236 237 You must call KSPSetUp() before calling PCFieldSplitGetSubPC(). 238 239 Level: advanced 240 241 .keywords: PC, 242 243 .seealso: PCFIELDSPLIT 244 @*/ 245 PetscErrorCode PCFieldSplitGetSubPC(PC pc,PetscInt *n,PC *subpc[]) 246 { 247 PetscErrorCode ierr,(*f)(PC,PetscInt*,PC **); 248 249 PetscFunctionBegin; 250 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 251 PetscValidIntPointer(n,2); 252 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitGetSubPC_C",(void (**)(void))&f);CHKERRQ(ierr); 253 if (f) { 254 ierr = (*f)(pc,n,subpc);CHKERRQ(ierr); 255 } else { 256 SETERRQ(PETSC_ERR_ARG_WRONG,"Cannot get subpc for this type of PC"); 257 } 258 PetscFunctionReturn(0); 259 } 260 261 /* -------------------------------------------------------------------------------------*/ 262 /*MC 263 PCFieldSplit - Preconditioner created by combining seperate preconditioners for individual 264 fields or groups of fields 265 266 267 To set options on the solvers for each block append -sub_ to all the PC 268 options database keys. For example, -sub_pc_type ilu -sub_pc_ilu_levels 1 269 270 To set the options on the solvers seperate for each block call PCFieldSplitGetSubPC() 271 and set the options directly on the resulting PC object 272 273 Level: intermediate 274 275 Concepts: physics based preconditioners 276 277 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 278 PCFieldSplitGetSubPC(), PCFieldSplitSetFields() 279 M*/ 280 281 EXTERN_C_BEGIN 282 #undef __FUNCT__ 283 #define __FUNCT__ "PCCreate_FieldSplit" 284 PetscErrorCode PCCreate_FieldSplit(PC pc) 285 { 286 PetscErrorCode ierr; 287 PC_FieldSplit *jac; 288 289 PetscFunctionBegin; 290 ierr = PetscNew(PC_FieldSplit,&jac);CHKERRQ(ierr); 291 PetscLogObjectMemory(pc,sizeof(PC_FieldSplit)); 292 ierr = PetscMemzero(jac,sizeof(PC_FieldSplit));CHKERRQ(ierr); 293 jac->bs = -1; 294 jac->nsplits = 0; 295 pc->data = (void*)jac; 296 297 pc->ops->apply = PCApply_FieldSplit; 298 pc->ops->setup = PCSetUp_FieldSplit; 299 pc->ops->destroy = PCDestroy_FieldSplit; 300 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 301 pc->ops->view = PCView_FieldSplit; 302 pc->ops->applyrichardson = 0; 303 304 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubPC_C","PCFieldSplitGetSubPC_FieldSplit", 305 PCFieldSplitGetSubPC_FieldSplit);CHKERRQ(ierr); 306 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","PCFieldSplitSetFields_FieldSplit", 307 PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 308 PetscFunctionReturn(0); 309 } 310 EXTERN_C_END 311 312 313