1 2 #include <petsc-private/pcimpl.h> /*I "petscpc.h" I*/ 3 #include <petscdmcomposite.h> /*I "petscdmcomposite.h" I*/ 4 5 const char *const PCFieldSplitSchurPreTypes[] = {"SELF","DIAG","USER","PCFieldSplitSchurPreType","PC_FIELDSPLIT_SCHUR_PRE_",0}; 6 const char *const PCFieldSplitSchurFactTypes[] = {"DIAG","LOWER","UPPER","FULL","PCFieldSplitSchurFactType","PC_FIELDSPLIT_SCHUR_FACT_",0}; 7 8 typedef struct _PC_FieldSplitLink *PC_FieldSplitLink; 9 struct _PC_FieldSplitLink { 10 KSP ksp; 11 Vec x,y; 12 char *splitname; 13 PetscInt nfields; 14 PetscInt *fields,*fields_col; 15 VecScatter sctx; 16 IS is,is_col; 17 PC_FieldSplitLink next,previous; 18 }; 19 20 typedef struct { 21 PCCompositeType type; 22 PetscBool defaultsplit; /* Flag for a system with a set of 'k' scalar fields with the same layout (and bs = k) */ 23 PetscBool splitdefined; /* Flag is set after the splits have been defined, to prevent more splits from being added */ 24 PetscBool realdiagonal; /* Flag to use the diagonal blocks of mat preconditioned by pmat, instead of just pmat */ 25 PetscInt bs; /* Block size for IS and Mat structures */ 26 PetscInt nsplits; /* Number of field divisions defined */ 27 Vec *x,*y,w1,w2; 28 Mat *mat; /* The diagonal block for each split */ 29 Mat *pmat; /* The preconditioning diagonal block for each split */ 30 Mat *Afield; /* The rows of the matrix associated with each split */ 31 PetscBool issetup; 32 /* Only used when Schur complement preconditioning is used */ 33 Mat B; /* The (0,1) block */ 34 Mat C; /* The (1,0) block */ 35 Mat schur; /* The Schur complement S = A11 - A10 A00^{-1} A01 */ 36 Mat schur_user; /* User-provided preconditioning matrix for the Schur complement */ 37 PCFieldSplitSchurPreType schurpre; /* Determines which preconditioning matrix is used for the Schur complement */ 38 PCFieldSplitSchurFactType schurfactorization; 39 KSP kspschur; /* The solver for S */ 40 PC_FieldSplitLink head; 41 PetscBool reset; /* indicates PCReset() has been last called on this object, hack */ 42 PetscBool suboptionsset; /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */ 43 } PC_FieldSplit; 44 45 /* 46 Notes: there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of 47 inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the 48 PC you could change this. 49 */ 50 51 /* This helper is so that setting a user-provided preconditioning matrix is orthogonal to choosing to use it. This way the 52 * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */ 53 static Mat FieldSplitSchurPre(PC_FieldSplit *jac) 54 { 55 switch (jac->schurpre) { 56 case PC_FIELDSPLIT_SCHUR_PRE_SELF: return jac->schur; 57 case PC_FIELDSPLIT_SCHUR_PRE_DIAG: return jac->pmat[1]; 58 case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */ 59 default: 60 return jac->schur_user ? jac->schur_user : jac->pmat[1]; 61 } 62 } 63 64 65 #undef __FUNCT__ 66 #define __FUNCT__ "PCView_FieldSplit" 67 static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer) 68 { 69 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 70 PetscErrorCode ierr; 71 PetscBool iascii; 72 PetscInt i,j; 73 PC_FieldSplitLink ilink = jac->head; 74 75 PetscFunctionBegin; 76 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 77 if (iascii) { 78 if (jac->bs > 0) { 79 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);CHKERRQ(ierr); 80 } else { 81 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D\n",PCCompositeTypes[jac->type],jac->nsplits);CHKERRQ(ierr); 82 } 83 ierr = PetscViewerASCIIPrintf(viewer," Solver info for each split is in the following KSP objects:\n");CHKERRQ(ierr); 84 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 85 for (i=0; i<jac->nsplits; i++) { 86 if (ilink->fields) { 87 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr); 88 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 89 for (j=0; j<ilink->nfields; j++) { 90 if (j > 0) { 91 ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr); 92 } 93 ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr); 94 } 95 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 96 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 97 } else { 98 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);CHKERRQ(ierr); 99 } 100 ierr = KSPView(ilink->ksp,viewer);CHKERRQ(ierr); 101 ilink = ilink->next; 102 } 103 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 104 } else { 105 SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for PCFieldSplit",((PetscObject)viewer)->type_name); 106 } 107 PetscFunctionReturn(0); 108 } 109 110 #undef __FUNCT__ 111 #define __FUNCT__ "PCView_FieldSplit_Schur" 112 static PetscErrorCode PCView_FieldSplit_Schur(PC pc,PetscViewer viewer) 113 { 114 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 115 PetscErrorCode ierr; 116 PetscBool iascii; 117 PetscInt i,j; 118 PC_FieldSplitLink ilink = jac->head; 119 KSP ksp; 120 121 PetscFunctionBegin; 122 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 123 if (iascii) { 124 if (jac->bs > 0) { 125 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, blocksize = %D, factorization %s\n",jac->bs,PCFieldSplitSchurFactTypes[jac->schurfactorization]);CHKERRQ(ierr); 126 } else { 127 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, factorization %s\n",PCFieldSplitSchurFactTypes[jac->schurfactorization]);CHKERRQ(ierr); 128 } 129 ierr = PetscViewerASCIIPrintf(viewer," Split info:\n");CHKERRQ(ierr); 130 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 131 for (i=0; i<jac->nsplits; i++) { 132 if (ilink->fields) { 133 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr); 134 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 135 for (j=0; j<ilink->nfields; j++) { 136 if (j > 0) { 137 ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr); 138 } 139 ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr); 140 } 141 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 142 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 143 } else { 144 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);CHKERRQ(ierr); 145 } 146 ilink = ilink->next; 147 } 148 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for A00 block \n");CHKERRQ(ierr); 149 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 150 if (jac->schur) { 151 ierr = MatSchurComplementGetKSP(jac->schur,&ksp);CHKERRQ(ierr); 152 ierr = KSPView(ksp,viewer);CHKERRQ(ierr); 153 } else { 154 ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr); 155 } 156 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 157 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for S = A11 - A10 inv(A00) A01 \n");CHKERRQ(ierr); 158 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 159 if (jac->kspschur) { 160 ierr = KSPView(jac->kspschur,viewer);CHKERRQ(ierr); 161 } else { 162 ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr); 163 } 164 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 165 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 166 } else { 167 SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for PCFieldSplit",((PetscObject)viewer)->type_name); 168 } 169 PetscFunctionReturn(0); 170 } 171 172 #undef __FUNCT__ 173 #define __FUNCT__ "PCFieldSplitSetRuntimeSplits_Private" 174 /* Precondition: jac->bs is set to a meaningful value */ 175 static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc) 176 { 177 PetscErrorCode ierr; 178 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 179 PetscInt i,nfields,*ifields,nfields_col,*ifields_col; 180 PetscBool flg,flg_col; 181 char optionname[128],splitname[8],optionname_col[128]; 182 183 PetscFunctionBegin; 184 ierr = PetscMalloc(jac->bs*sizeof(PetscInt),&ifields);CHKERRQ(ierr); 185 ierr = PetscMalloc(jac->bs*sizeof(PetscInt),&ifields_col);CHKERRQ(ierr); 186 for (i=0,flg=PETSC_TRUE; ; i++) { 187 ierr = PetscSNPrintf(splitname,sizeof splitname,"%D",i);CHKERRQ(ierr); 188 ierr = PetscSNPrintf(optionname,sizeof optionname,"-pc_fieldsplit_%D_fields",i);CHKERRQ(ierr); 189 ierr = PetscSNPrintf(optionname_col,sizeof optionname_col,"-pc_fieldsplit_%D_fields_col",i);CHKERRQ(ierr); 190 nfields = jac->bs; 191 nfields_col = jac->bs; 192 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->prefix,optionname,ifields,&nfields,&flg);CHKERRQ(ierr); 193 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->prefix,optionname_col,ifields_col,&nfields_col,&flg_col);CHKERRQ(ierr); 194 if (!flg) break; 195 else if (flg && !flg_col) { 196 if (!nfields) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields"); 197 ierr = PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields);CHKERRQ(ierr); 198 } 199 else { 200 if (!nfields || !nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields"); 201 if (nfields != nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Number of row and column fields must match"); 202 ierr = PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields_col);CHKERRQ(ierr); 203 } 204 } 205 if (i > 0) { 206 /* Makes command-line setting of splits take precedence over setting them in code. 207 Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would 208 create new splits, which would probably not be what the user wanted. */ 209 jac->splitdefined = PETSC_TRUE; 210 } 211 ierr = PetscFree(ifields);CHKERRQ(ierr); 212 ierr = PetscFree(ifields_col);CHKERRQ(ierr); 213 PetscFunctionReturn(0); 214 } 215 216 #undef __FUNCT__ 217 #define __FUNCT__ "PCFieldSplitSetDefaults" 218 static PetscErrorCode PCFieldSplitSetDefaults(PC pc) 219 { 220 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 221 PetscErrorCode ierr; 222 PC_FieldSplitLink ilink = jac->head; 223 PetscBool flg = PETSC_FALSE,stokes = PETSC_FALSE; 224 PetscInt i; 225 226 PetscFunctionBegin; 227 if (!ilink) { 228 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,PETSC_NULL);CHKERRQ(ierr); 229 if (pc->dm && !stokes) { 230 PetscInt numFields, f; 231 char **fieldNames; 232 IS *fields; 233 DM *dms; 234 ierr = DMCreateDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms); CHKERRQ(ierr); 235 for(f = 0; f < numFields; ++f) { 236 ierr = PCFieldSplitSetIS(pc, fieldNames[f], fields[f]);CHKERRQ(ierr); 237 ierr = PetscFree(fieldNames[f]);CHKERRQ(ierr); 238 ierr = ISDestroy(&fields[f]);CHKERRQ(ierr); 239 } 240 ierr = PetscFree(fieldNames);CHKERRQ(ierr); 241 ierr = PetscFree(fields);CHKERRQ(ierr); 242 if(dms) { 243 ierr = PetscInfo(pc,"Setting up physics based fieldsplit preconditioner using the embedded DM\n");CHKERRQ(ierr); 244 for (ilink=jac->head,i=0; ilink; ilink=ilink->next,i++) { 245 ierr = KSPSetDM(ilink->ksp, dms[i]);CHKERRQ(ierr); 246 ierr = KSPSetDMActive(ilink->ksp, PETSC_FALSE);CHKERRQ(ierr); 247 ierr = DMDestroy(&dms[i]); CHKERRQ(ierr); 248 } 249 ierr = PetscFree(dms);CHKERRQ(ierr); 250 } 251 } else { 252 if (jac->bs <= 0) { 253 if (pc->pmat) { 254 ierr = MatGetBlockSize(pc->pmat,&jac->bs);CHKERRQ(ierr); 255 } else { 256 jac->bs = 1; 257 } 258 } 259 260 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&flg,PETSC_NULL);CHKERRQ(ierr); 261 if (stokes) { 262 IS zerodiags,rest; 263 PetscInt nmin,nmax; 264 265 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 266 ierr = MatFindZeroDiagonals(pc->mat,&zerodiags);CHKERRQ(ierr); 267 ierr = ISComplement(zerodiags,nmin,nmax,&rest);CHKERRQ(ierr); 268 ierr = PCFieldSplitSetIS(pc,"0",rest);CHKERRQ(ierr); 269 ierr = PCFieldSplitSetIS(pc,"1",zerodiags);CHKERRQ(ierr); 270 ierr = ISDestroy(&zerodiags);CHKERRQ(ierr); 271 ierr = ISDestroy(&rest);CHKERRQ(ierr); 272 } else { 273 if (!flg) { 274 /* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit() 275 then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */ 276 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 277 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 278 } 279 if (flg || !jac->splitdefined) { 280 ierr = PetscInfo(pc,"Using default splitting of fields\n");CHKERRQ(ierr); 281 for (i=0; i<jac->bs; i++) { 282 char splitname[8]; 283 ierr = PetscSNPrintf(splitname,sizeof splitname,"%D",i);CHKERRQ(ierr); 284 ierr = PCFieldSplitSetFields(pc,splitname,1,&i,&i);CHKERRQ(ierr); 285 } 286 jac->defaultsplit = PETSC_TRUE; 287 } 288 } 289 } 290 } else if (jac->nsplits == 1) { 291 if (ilink->is) { 292 IS is2; 293 PetscInt nmin,nmax; 294 295 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 296 ierr = ISComplement(ilink->is,nmin,nmax,&is2);CHKERRQ(ierr); 297 ierr = PCFieldSplitSetIS(pc,"1",is2);CHKERRQ(ierr); 298 ierr = ISDestroy(&is2);CHKERRQ(ierr); 299 } else SETERRQ(((PetscObject) pc)->comm,PETSC_ERR_SUP,"Must provide at least two sets of fields to PCFieldSplit()"); 300 } else if (jac->reset) { 301 /* PCReset() has been called on this PC, ilink exists but all IS and Vec data structures in it must be rebuilt 302 This is basically the !ilink portion of code above copied from above and the allocation of the ilinks removed 303 since they already exist. This should be totally rewritten */ 304 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,PETSC_NULL);CHKERRQ(ierr); 305 if (pc->dm && !stokes) { 306 PetscBool dmcomposite; 307 ierr = PetscObjectTypeCompare((PetscObject)pc->dm,DMCOMPOSITE,&dmcomposite);CHKERRQ(ierr); 308 if (dmcomposite) { 309 PetscInt nDM; 310 IS *fields; 311 DM *dms; 312 ierr = PetscInfo(pc,"Setting up physics based fieldsplit preconditioner using the embedded DM\n");CHKERRQ(ierr); 313 ierr = DMCompositeGetNumberDM(pc->dm,&nDM);CHKERRQ(ierr); 314 ierr = DMCompositeGetGlobalISs(pc->dm,&fields);CHKERRQ(ierr); 315 ierr = PetscMalloc(nDM*sizeof(DM),&dms);CHKERRQ(ierr); 316 ierr = DMCompositeGetEntriesArray(pc->dm,dms);CHKERRQ(ierr); 317 for (i=0; i<nDM; i++) { 318 ierr = KSPSetDM(ilink->ksp,dms[i]);CHKERRQ(ierr); 319 ierr = KSPSetDMActive(ilink->ksp,PETSC_FALSE);CHKERRQ(ierr); 320 ilink->is = fields[i]; 321 ilink = ilink->next; 322 } 323 ierr = PetscFree(fields);CHKERRQ(ierr); 324 ierr = PetscFree(dms);CHKERRQ(ierr); 325 } 326 } else { 327 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&flg,PETSC_NULL);CHKERRQ(ierr); 328 if (stokes) { 329 IS zerodiags,rest; 330 PetscInt nmin,nmax; 331 332 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 333 ierr = MatFindZeroDiagonals(pc->mat,&zerodiags);CHKERRQ(ierr); 334 ierr = ISComplement(zerodiags,nmin,nmax,&rest);CHKERRQ(ierr); 335 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 336 ierr = ISDestroy(&ilink->next->is);CHKERRQ(ierr); 337 ilink->is = rest; 338 ilink->next->is = zerodiags; 339 } else SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Cases not yet handled when PCReset() was used"); 340 } 341 } 342 343 if (jac->nsplits < 2) SETERRQ1(((PetscObject) pc)->comm,PETSC_ERR_PLIB,"Unhandled case, must have at least two fields, not %d", jac->nsplits); 344 PetscFunctionReturn(0); 345 } 346 347 #undef __FUNCT__ 348 #define __FUNCT__ "PCSetUp_FieldSplit" 349 static PetscErrorCode PCSetUp_FieldSplit(PC pc) 350 { 351 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 352 PetscErrorCode ierr; 353 PC_FieldSplitLink ilink; 354 PetscInt i,nsplit; 355 MatStructure flag = pc->flag; 356 PetscBool sorted, sorted_col; 357 358 PetscFunctionBegin; 359 ierr = PCFieldSplitSetDefaults(pc);CHKERRQ(ierr); 360 nsplit = jac->nsplits; 361 ilink = jac->head; 362 363 /* get the matrices for each split */ 364 if (!jac->issetup) { 365 PetscInt rstart,rend,nslots,bs; 366 367 jac->issetup = PETSC_TRUE; 368 369 /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */ 370 if (jac->defaultsplit || !ilink->is) { 371 if (jac->bs <= 0) jac->bs = nsplit; 372 } 373 bs = jac->bs; 374 ierr = MatGetOwnershipRange(pc->pmat,&rstart,&rend);CHKERRQ(ierr); 375 nslots = (rend - rstart)/bs; 376 for (i=0; i<nsplit; i++) { 377 if (jac->defaultsplit) { 378 ierr = ISCreateStride(((PetscObject)pc)->comm,nslots,rstart+i,nsplit,&ilink->is);CHKERRQ(ierr); 379 ierr = ISDuplicate(ilink->is,&ilink->is_col);CHKERRQ(ierr); 380 } else if (!ilink->is) { 381 if (ilink->nfields > 1) { 382 PetscInt *ii,*jj,j,k,nfields = ilink->nfields,*fields = ilink->fields,*fields_col = ilink->fields_col; 383 ierr = PetscMalloc(ilink->nfields*nslots*sizeof(PetscInt),&ii);CHKERRQ(ierr); 384 ierr = PetscMalloc(ilink->nfields*nslots*sizeof(PetscInt),&jj);CHKERRQ(ierr); 385 for (j=0; j<nslots; j++) { 386 for (k=0; k<nfields; k++) { 387 ii[nfields*j + k] = rstart + bs*j + fields[k]; 388 jj[nfields*j + k] = rstart + bs*j + fields_col[k]; 389 } 390 } 391 ierr = ISCreateGeneral(((PetscObject)pc)->comm,nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is);CHKERRQ(ierr); 392 ierr = ISCreateGeneral(((PetscObject)pc)->comm,nslots*nfields,jj,PETSC_OWN_POINTER,&ilink->is_col);CHKERRQ(ierr); 393 } else { 394 ierr = ISCreateStride(((PetscObject)pc)->comm,nslots,rstart+ilink->fields[0],bs,&ilink->is);CHKERRQ(ierr); 395 ierr = ISCreateStride(((PetscObject)pc)->comm,nslots,rstart+ilink->fields_col[0],bs,&ilink->is_col);CHKERRQ(ierr); 396 } 397 } 398 ierr = ISSorted(ilink->is,&sorted);CHKERRQ(ierr); 399 if (ilink->is_col) { ierr = ISSorted(ilink->is_col,&sorted_col);CHKERRQ(ierr); } 400 if (!sorted || !sorted_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Fields must be sorted when creating split"); 401 ilink = ilink->next; 402 } 403 } 404 405 ilink = jac->head; 406 if (!jac->pmat) { 407 Vec xtmp; 408 409 ierr = MatGetVecs(pc->pmat,&xtmp,PETSC_NULL);CHKERRQ(ierr); 410 ierr = PetscMalloc(nsplit*sizeof(Mat),&jac->pmat);CHKERRQ(ierr); 411 ierr = PetscMalloc2(nsplit,Vec,&jac->x,nsplit,Vec,&jac->y);CHKERRQ(ierr); 412 for (i=0; i<nsplit; i++) { 413 MatNullSpace sp; 414 415 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 416 /* create work vectors for each split */ 417 ierr = MatGetVecs(jac->pmat[i],&jac->x[i],&jac->y[i]);CHKERRQ(ierr); 418 ilink->x = jac->x[i]; ilink->y = jac->y[i]; 419 /* compute scatter contexts needed by multiplicative versions and non-default splits */ 420 ierr = VecScatterCreate(xtmp,ilink->is,jac->x[i],PETSC_NULL,&ilink->sctx);CHKERRQ(ierr); 421 /* HACK: Check for the constant null space */ 422 ierr = MatGetNullSpace(pc->pmat, &sp);CHKERRQ(ierr); 423 if (sp) { 424 MatNullSpace subsp; 425 Vec ftmp, gtmp; 426 PetscReal norm; 427 PetscInt N; 428 429 ierr = MatGetVecs(pc->pmat, >mp, PETSC_NULL);CHKERRQ(ierr); 430 ierr = MatGetVecs(jac->pmat[i], &ftmp, PETSC_NULL);CHKERRQ(ierr); 431 ierr = VecGetSize(ftmp, &N);CHKERRQ(ierr); 432 ierr = VecSet(ftmp, 1.0/N);CHKERRQ(ierr); 433 ierr = VecScatterBegin(ilink->sctx, ftmp, gtmp, INSERT_VALUES, SCATTER_REVERSE);CHKERRQ(ierr); 434 ierr = VecScatterEnd(ilink->sctx, ftmp, gtmp, INSERT_VALUES, SCATTER_REVERSE);CHKERRQ(ierr); 435 ierr = MatNullSpaceRemove(sp, gtmp, PETSC_NULL);CHKERRQ(ierr); 436 ierr = VecNorm(gtmp, NORM_2, &norm);CHKERRQ(ierr); 437 if (norm < 1.0e-10) { 438 ierr = MatNullSpaceCreate(((PetscObject)pc)->comm, PETSC_TRUE, 0, PETSC_NULL, &subsp);CHKERRQ(ierr); 439 ierr = MatSetNullSpace(jac->pmat[i], subsp);CHKERRQ(ierr); 440 ierr = MatNullSpaceDestroy(&subsp);CHKERRQ(ierr); 441 } 442 ierr = VecDestroy(&ftmp);CHKERRQ(ierr); 443 ierr = VecDestroy(>mp);CHKERRQ(ierr); 444 } 445 /* Check for null space attached to IS */ 446 ierr = PetscObjectQuery((PetscObject) ilink->is, "nearnullspace", (PetscObject *) &sp);CHKERRQ(ierr); 447 if (sp) { 448 ierr = MatSetNearNullSpace(jac->pmat[i], sp);CHKERRQ(ierr); 449 } 450 ilink = ilink->next; 451 } 452 ierr = VecDestroy(&xtmp);CHKERRQ(ierr); 453 } else { 454 for (i=0; i<nsplit; i++) { 455 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 456 ilink = ilink->next; 457 } 458 } 459 if (jac->realdiagonal) { 460 ilink = jac->head; 461 if (!jac->mat) { 462 ierr = PetscMalloc(nsplit*sizeof(Mat),&jac->mat);CHKERRQ(ierr); 463 for (i=0; i<nsplit; i++) { 464 ierr = MatGetSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->mat[i]);CHKERRQ(ierr); 465 ilink = ilink->next; 466 } 467 } else { 468 for (i=0; i<nsplit; i++) { 469 if (jac->mat[i]) {ierr = MatGetSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->mat[i]);CHKERRQ(ierr);} 470 ilink = ilink->next; 471 } 472 } 473 } else { 474 jac->mat = jac->pmat; 475 } 476 477 if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR) { 478 /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */ 479 ilink = jac->head; 480 if (!jac->Afield) { 481 ierr = PetscMalloc(nsplit*sizeof(Mat),&jac->Afield);CHKERRQ(ierr); 482 for (i=0; i<nsplit; i++) { 483 ierr = MatGetSubMatrix(pc->mat,ilink->is,PETSC_NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 484 ilink = ilink->next; 485 } 486 } else { 487 for (i=0; i<nsplit; i++) { 488 ierr = MatGetSubMatrix(pc->mat,ilink->is,PETSC_NULL,MAT_REUSE_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 489 ilink = ilink->next; 490 } 491 } 492 } 493 494 if (jac->type == PC_COMPOSITE_SCHUR) { 495 IS ccis; 496 PetscInt rstart,rend; 497 char lscname[256]; 498 PetscObject LSC_L; 499 if (nsplit != 2) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields"); 500 501 /* When extracting off-diagonal submatrices, we take complements from this range */ 502 ierr = MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);CHKERRQ(ierr); 503 504 /* need to handle case when one is resetting up the preconditioner */ 505 if (jac->schur) { 506 ilink = jac->head; 507 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 508 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);CHKERRQ(ierr); 509 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 510 ilink = ilink->next; 511 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 512 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);CHKERRQ(ierr); 513 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 514 ierr = MatSchurComplementUpdate(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->pmat[1],pc->flag);CHKERRQ(ierr); 515 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),pc->flag);CHKERRQ(ierr); 516 517 } else { 518 KSP ksp; 519 char schurprefix[256]; 520 MatNullSpace sp; 521 522 /* extract the A01 and A10 matrices */ 523 ilink = jac->head; 524 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 525 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr); 526 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 527 ilink = ilink->next; 528 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 529 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 530 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 531 /* Use mat[0] (diagonal block of the real matrix) preconditioned by pmat[0] */ 532 ierr = MatCreateSchurComplement(jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1],&jac->schur);CHKERRQ(ierr); 533 ierr = MatGetNullSpace(jac->pmat[1], &sp);CHKERRQ(ierr); 534 if (sp) {ierr = MatSetNullSpace(jac->schur, sp);CHKERRQ(ierr);} 535 /* set tabbing and options prefix of KSP inside the MatSchur */ 536 ierr = MatSchurComplementGetKSP(jac->schur,&ksp);CHKERRQ(ierr); 537 ierr = PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,2);CHKERRQ(ierr); 538 ierr = PetscSNPrintf(schurprefix,sizeof schurprefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",jac->head->splitname);CHKERRQ(ierr); 539 ierr = KSPSetOptionsPrefix(ksp,schurprefix);CHKERRQ(ierr); 540 /* Need to call this everytime because new matrix is being created */ 541 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr); 542 ierr = MatSetUp(jac->schur);CHKERRQ(ierr); 543 544 ierr = KSPCreate(((PetscObject)pc)->comm,&jac->kspschur);CHKERRQ(ierr); 545 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);CHKERRQ(ierr); 546 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);CHKERRQ(ierr); 547 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 548 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) { 549 PC pc; 550 ierr = KSPGetPC(jac->kspschur,&pc);CHKERRQ(ierr); 551 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 552 /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */ 553 } 554 ierr = PetscSNPrintf(schurprefix,sizeof schurprefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 555 ierr = KSPSetOptionsPrefix(jac->kspschur,schurprefix);CHKERRQ(ierr); 556 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 557 /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */ 558 ierr = KSPSetFromOptions(jac->kspschur);CHKERRQ(ierr); 559 } 560 561 /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */ 562 ierr = PetscSNPrintf(lscname,sizeof lscname,"%s_LSC_L",ilink->splitname);CHKERRQ(ierr); 563 ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 564 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 565 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);CHKERRQ(ierr);} 566 ierr = PetscSNPrintf(lscname,sizeof lscname,"%s_LSC_Lp",ilink->splitname);CHKERRQ(ierr); 567 ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 568 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 569 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);CHKERRQ(ierr);} 570 } else { 571 /* set up the individual PCs */ 572 i = 0; 573 ilink = jac->head; 574 while (ilink) { 575 ierr = KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i],flag);CHKERRQ(ierr); 576 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 577 if (!jac->suboptionsset) {ierr = KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr);} 578 ierr = KSPSetUp(ilink->ksp);CHKERRQ(ierr); 579 i++; 580 ilink = ilink->next; 581 } 582 } 583 584 jac->suboptionsset = PETSC_TRUE; 585 PetscFunctionReturn(0); 586 } 587 588 #define FieldSplitSplitSolveAdd(ilink,xx,yy) \ 589 (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 590 VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 591 KSPSolve(ilink->ksp,ilink->x,ilink->y) || \ 592 VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \ 593 VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE)) 594 595 #undef __FUNCT__ 596 #define __FUNCT__ "PCApply_FieldSplit_Schur" 597 static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y) 598 { 599 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 600 PetscErrorCode ierr; 601 KSP ksp; 602 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 603 604 PetscFunctionBegin; 605 ierr = MatSchurComplementGetKSP(jac->schur,&ksp);CHKERRQ(ierr); 606 607 switch (jac->schurfactorization) { 608 case PC_FIELDSPLIT_SCHUR_FACT_DIAG: 609 /* [A00 0; 0 -S], positive definite, suitable for MINRES */ 610 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 611 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 612 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 613 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 614 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 615 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 616 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 617 ierr = VecScale(ilinkD->y,-1.);CHKERRQ(ierr); 618 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 619 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 620 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 621 break; 622 case PC_FIELDSPLIT_SCHUR_FACT_LOWER: 623 /* [A00 0; A10 S], suitable for left preconditioning */ 624 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 625 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 626 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 627 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 628 ierr = VecScale(ilinkD->x,-1.);CHKERRQ(ierr); 629 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 630 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 631 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 632 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 633 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 634 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 635 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 636 break; 637 case PC_FIELDSPLIT_SCHUR_FACT_UPPER: 638 /* [A00 A01; 0 S], suitable for right preconditioning */ 639 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 640 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 641 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 642 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 643 ierr = VecScale(ilinkA->x,-1.);CHKERRQ(ierr); 644 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 645 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 646 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 647 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 648 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 649 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 650 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 651 break; 652 case PC_FIELDSPLIT_SCHUR_FACT_FULL: 653 /* [1 0; A10 A00^{-1} 1] [A00 0; 0 S] [1 A00^{-1}A01; 0 1], an exact solve if applied exactly, needs one extra solve with A */ 654 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 655 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 656 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 657 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 658 ierr = VecScale(ilinkD->x,-1.0);CHKERRQ(ierr); 659 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 660 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 661 662 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 663 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 664 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 665 666 ierr = MatMult(jac->B,ilinkD->y,ilinkA->y);CHKERRQ(ierr); 667 ierr = VecAXPY(ilinkA->x,-1.0,ilinkA->y);CHKERRQ(ierr); 668 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 669 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 670 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 671 } 672 PetscFunctionReturn(0); 673 } 674 675 #undef __FUNCT__ 676 #define __FUNCT__ "PCApply_FieldSplit" 677 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y) 678 { 679 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 680 PetscErrorCode ierr; 681 PC_FieldSplitLink ilink = jac->head; 682 PetscInt cnt,bs; 683 684 PetscFunctionBegin; 685 CHKMEMQ; 686 if (jac->type == PC_COMPOSITE_ADDITIVE) { 687 if (jac->defaultsplit) { 688 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 689 if (jac->bs > 0 && bs != jac->bs) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Blocksize of x vector %D does not match fieldsplit blocksize %D",bs,jac->bs); 690 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 691 if (jac->bs > 0 && bs != jac->bs) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Blocksize of y vector %D does not match fieldsplit blocksize %D",bs,jac->bs); 692 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 693 while (ilink) { 694 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 695 ilink = ilink->next; 696 } 697 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 698 } else { 699 ierr = VecSet(y,0.0);CHKERRQ(ierr); 700 while (ilink) { 701 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 702 ilink = ilink->next; 703 } 704 } 705 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 706 if (!jac->w1) { 707 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 708 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 709 } 710 ierr = VecSet(y,0.0);CHKERRQ(ierr); 711 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 712 cnt = 1; 713 while (ilink->next) { 714 ilink = ilink->next; 715 /* compute the residual only over the part of the vector needed */ 716 ierr = MatMult(jac->Afield[cnt++],y,ilink->x);CHKERRQ(ierr); 717 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 718 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 719 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 720 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 721 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 722 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 723 } 724 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 725 cnt -= 2; 726 while (ilink->previous) { 727 ilink = ilink->previous; 728 /* compute the residual only over the part of the vector needed */ 729 ierr = MatMult(jac->Afield[cnt--],y,ilink->x);CHKERRQ(ierr); 730 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 731 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 732 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 733 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 734 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 735 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 736 } 737 } 738 } else SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type); 739 CHKMEMQ; 740 PetscFunctionReturn(0); 741 } 742 743 #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \ 744 (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 745 VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 746 KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \ 747 VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \ 748 VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE)) 749 750 #undef __FUNCT__ 751 #define __FUNCT__ "PCApplyTranspose_FieldSplit" 752 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y) 753 { 754 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 755 PetscErrorCode ierr; 756 PC_FieldSplitLink ilink = jac->head; 757 PetscInt bs; 758 759 PetscFunctionBegin; 760 CHKMEMQ; 761 if (jac->type == PC_COMPOSITE_ADDITIVE) { 762 if (jac->defaultsplit) { 763 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 764 if (jac->bs > 0 && bs != jac->bs) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Blocksize of x vector %D does not match fieldsplit blocksize %D",bs,jac->bs); 765 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 766 if (jac->bs > 0 && bs != jac->bs) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Blocksize of y vector %D does not match fieldsplit blocksize %D",bs,jac->bs); 767 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 768 while (ilink) { 769 ierr = KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 770 ilink = ilink->next; 771 } 772 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 773 } else { 774 ierr = VecSet(y,0.0);CHKERRQ(ierr); 775 while (ilink) { 776 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 777 ilink = ilink->next; 778 } 779 } 780 } else { 781 if (!jac->w1) { 782 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 783 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 784 } 785 ierr = VecSet(y,0.0);CHKERRQ(ierr); 786 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 787 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 788 while (ilink->next) { 789 ilink = ilink->next; 790 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 791 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 792 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 793 } 794 while (ilink->previous) { 795 ilink = ilink->previous; 796 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 797 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 798 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 799 } 800 } else { 801 while (ilink->next) { /* get to last entry in linked list */ 802 ilink = ilink->next; 803 } 804 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 805 while (ilink->previous) { 806 ilink = ilink->previous; 807 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 808 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 809 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 810 } 811 } 812 } 813 CHKMEMQ; 814 PetscFunctionReturn(0); 815 } 816 817 #undef __FUNCT__ 818 #define __FUNCT__ "PCReset_FieldSplit" 819 static PetscErrorCode PCReset_FieldSplit(PC pc) 820 { 821 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 822 PetscErrorCode ierr; 823 PC_FieldSplitLink ilink = jac->head,next; 824 825 PetscFunctionBegin; 826 while (ilink) { 827 ierr = KSPReset(ilink->ksp);CHKERRQ(ierr); 828 ierr = VecDestroy(&ilink->x);CHKERRQ(ierr); 829 ierr = VecDestroy(&ilink->y);CHKERRQ(ierr); 830 ierr = VecScatterDestroy(&ilink->sctx);CHKERRQ(ierr); 831 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 832 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 833 next = ilink->next; 834 ilink = next; 835 } 836 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 837 if (jac->mat && jac->mat != jac->pmat) { 838 ierr = MatDestroyMatrices(jac->nsplits,&jac->mat);CHKERRQ(ierr); 839 } else if (jac->mat) { 840 jac->mat = PETSC_NULL; 841 } 842 if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);} 843 if (jac->Afield) {ierr = MatDestroyMatrices(jac->nsplits,&jac->Afield);CHKERRQ(ierr);} 844 ierr = VecDestroy(&jac->w1);CHKERRQ(ierr); 845 ierr = VecDestroy(&jac->w2);CHKERRQ(ierr); 846 ierr = MatDestroy(&jac->schur);CHKERRQ(ierr); 847 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 848 ierr = KSPDestroy(&jac->kspschur);CHKERRQ(ierr); 849 ierr = MatDestroy(&jac->B);CHKERRQ(ierr); 850 ierr = MatDestroy(&jac->C);CHKERRQ(ierr); 851 jac->reset = PETSC_TRUE; 852 PetscFunctionReturn(0); 853 } 854 855 #undef __FUNCT__ 856 #define __FUNCT__ "PCDestroy_FieldSplit" 857 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 858 { 859 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 860 PetscErrorCode ierr; 861 PC_FieldSplitLink ilink = jac->head,next; 862 863 PetscFunctionBegin; 864 ierr = PCReset_FieldSplit(pc);CHKERRQ(ierr); 865 while (ilink) { 866 ierr = KSPDestroy(&ilink->ksp);CHKERRQ(ierr); 867 next = ilink->next; 868 ierr = PetscFree(ilink->splitname);CHKERRQ(ierr); 869 ierr = PetscFree(ilink->fields);CHKERRQ(ierr); 870 ierr = PetscFree(ilink->fields_col);CHKERRQ(ierr); 871 ierr = PetscFree(ilink);CHKERRQ(ierr); 872 ilink = next; 873 } 874 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 875 ierr = PetscFree(pc->data);CHKERRQ(ierr); 876 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","",PETSC_NULL);CHKERRQ(ierr); 877 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","",PETSC_NULL);CHKERRQ(ierr); 878 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","",PETSC_NULL);CHKERRQ(ierr); 879 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","",PETSC_NULL);CHKERRQ(ierr); 880 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","",PETSC_NULL);CHKERRQ(ierr); 881 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","",PETSC_NULL);CHKERRQ(ierr); 882 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetSchurFactType_C","",PETSC_NULL);CHKERRQ(ierr); 883 PetscFunctionReturn(0); 884 } 885 886 #undef __FUNCT__ 887 #define __FUNCT__ "PCSetFromOptions_FieldSplit" 888 static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc) 889 { 890 PetscErrorCode ierr; 891 PetscInt bs; 892 PetscBool flg,stokes = PETSC_FALSE; 893 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 894 PCCompositeType ctype; 895 DM ddm; 896 char ddm_name[1025]; 897 898 PetscFunctionBegin; 899 ierr = PetscOptionsHead("FieldSplit options");CHKERRQ(ierr); 900 if(pc->dm) { 901 /* Allow the user to request a decomposition DM by name */ 902 ierr = PetscStrncpy(ddm_name, "", 1024); CHKERRQ(ierr); 903 ierr = PetscOptionsString("-pc_fieldsplit_decomposition_dm", "Name of the DM defining the composition", "PCSetDM", ddm_name, ddm_name,1024,&flg); CHKERRQ(ierr); 904 if(flg) { 905 ierr = DMCreateDecompositionDM(pc->dm, ddm_name, &ddm); CHKERRQ(ierr); 906 if(!ddm) { 907 SETERRQ1(((PetscObject)pc)->comm, PETSC_ERR_ARG_WRONGSTATE, "Uknown DM decomposition name %s", ddm_name); 908 } 909 ierr = PetscInfo(pc,"Using decomposition DM defined using options database\n");CHKERRQ(ierr); 910 ierr = PCSetDM(pc,ddm); CHKERRQ(ierr); 911 } 912 } 913 ierr = PetscOptionsBool("-pc_fieldsplit_real_diagonal","Use diagonal blocks of the operator","PCFieldSplitSetRealDiagonal",jac->realdiagonal,&jac->realdiagonal,PETSC_NULL);CHKERRQ(ierr); 914 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 915 if (flg) { 916 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 917 } 918 919 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,PETSC_NULL);CHKERRQ(ierr); 920 if (stokes) { 921 ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);CHKERRQ(ierr); 922 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF; 923 } 924 925 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr); 926 if (flg) { 927 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr); 928 } 929 /* Only setup fields once */ 930 if ((jac->bs > 0) && (jac->nsplits == 0)) { 931 /* only allow user to set fields from command line if bs is already known. 932 otherwise user can set them in PCFieldSplitSetDefaults() */ 933 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 934 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 935 } 936 if (jac->type == PC_COMPOSITE_SCHUR) { 937 ierr = PetscOptionsGetEnum(((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);CHKERRQ(ierr); 938 if (flg) {ierr = PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");CHKERRQ(ierr);} 939 ierr = PetscOptionsEnum("-pc_fieldsplit_schur_fact_type","Which off-diagonal parts of the block factorization to use","PCFieldSplitSetSchurFactType",PCFieldSplitSchurFactTypes,(PetscEnum)jac->schurfactorization,(PetscEnum*)&jac->schurfactorization,PETSC_NULL);CHKERRQ(ierr); 940 ierr = PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSchurPrecondition",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,PETSC_NULL);CHKERRQ(ierr); 941 } 942 ierr = PetscOptionsTail();CHKERRQ(ierr); 943 PetscFunctionReturn(0); 944 } 945 946 /*------------------------------------------------------------------------------------*/ 947 948 EXTERN_C_BEGIN 949 #undef __FUNCT__ 950 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 951 PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 952 { 953 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 954 PetscErrorCode ierr; 955 PC_FieldSplitLink ilink,next = jac->head; 956 char prefix[128]; 957 PetscInt i; 958 959 PetscFunctionBegin; 960 if (jac->splitdefined) { 961 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 962 PetscFunctionReturn(0); 963 } 964 for (i=0; i<n; i++) { 965 if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 966 if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 967 } 968 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 969 if (splitname) { 970 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 971 } else { 972 ierr = PetscMalloc(3*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 973 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr); 974 } 975 ierr = PetscMalloc(n*sizeof(PetscInt),&ilink->fields);CHKERRQ(ierr); 976 ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 977 ierr = PetscMalloc(n*sizeof(PetscInt),&ilink->fields_col);CHKERRQ(ierr); 978 ierr = PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));CHKERRQ(ierr); 979 ilink->nfields = n; 980 ilink->next = PETSC_NULL; 981 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 982 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 983 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 984 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 985 986 ierr = PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 987 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 988 989 if (!next) { 990 jac->head = ilink; 991 ilink->previous = PETSC_NULL; 992 } else { 993 while (next->next) { 994 next = next->next; 995 } 996 next->next = ilink; 997 ilink->previous = next; 998 } 999 jac->nsplits++; 1000 PetscFunctionReturn(0); 1001 } 1002 EXTERN_C_END 1003 1004 EXTERN_C_BEGIN 1005 #undef __FUNCT__ 1006 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit_Schur" 1007 PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 1008 { 1009 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1010 PetscErrorCode ierr; 1011 1012 PetscFunctionBegin; 1013 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 1014 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr); 1015 (*subksp)[1] = jac->kspschur; 1016 if (n) *n = jac->nsplits; 1017 PetscFunctionReturn(0); 1018 } 1019 EXTERN_C_END 1020 1021 EXTERN_C_BEGIN 1022 #undef __FUNCT__ 1023 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 1024 PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1025 { 1026 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1027 PetscErrorCode ierr; 1028 PetscInt cnt = 0; 1029 PC_FieldSplitLink ilink = jac->head; 1030 1031 PetscFunctionBegin; 1032 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 1033 while (ilink) { 1034 (*subksp)[cnt++] = ilink->ksp; 1035 ilink = ilink->next; 1036 } 1037 if (cnt != jac->nsplits) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt PCFIELDSPLIT object: number of splits in linked list %D does not match number in object %D",cnt,jac->nsplits); 1038 if (n) *n = jac->nsplits; 1039 PetscFunctionReturn(0); 1040 } 1041 EXTERN_C_END 1042 1043 EXTERN_C_BEGIN 1044 #undef __FUNCT__ 1045 #define __FUNCT__ "PCFieldSplitSetIS_FieldSplit" 1046 PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 1047 { 1048 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1049 PetscErrorCode ierr; 1050 PC_FieldSplitLink ilink, next = jac->head; 1051 char prefix[128]; 1052 1053 PetscFunctionBegin; 1054 if (jac->splitdefined) { 1055 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1056 PetscFunctionReturn(0); 1057 } 1058 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 1059 if (splitname) { 1060 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1061 } else { 1062 ierr = PetscMalloc(8*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 1063 ierr = PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);CHKERRQ(ierr); 1064 } 1065 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1066 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1067 ilink->is = is; 1068 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1069 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1070 ilink->is_col = is; 1071 ilink->next = PETSC_NULL; 1072 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 1073 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1074 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1075 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1076 1077 ierr = PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 1078 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1079 1080 if (!next) { 1081 jac->head = ilink; 1082 ilink->previous = PETSC_NULL; 1083 } else { 1084 while (next->next) { 1085 next = next->next; 1086 } 1087 next->next = ilink; 1088 ilink->previous = next; 1089 } 1090 jac->nsplits++; 1091 1092 PetscFunctionReturn(0); 1093 } 1094 EXTERN_C_END 1095 1096 #undef __FUNCT__ 1097 #define __FUNCT__ "PCFieldSplitSetFields" 1098 /*@ 1099 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 1100 1101 Logically Collective on PC 1102 1103 Input Parameters: 1104 + pc - the preconditioner context 1105 . splitname - name of this split, if PETSC_NULL the number of the split is used 1106 . n - the number of fields in this split 1107 - fields - the fields in this split 1108 1109 Level: intermediate 1110 1111 Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 1112 1113 The PCFieldSplitSetFields() is for defining fields as a strided blocks. For example, if the block 1114 size is three then one can define a field as 0, or 1 or 2 or 0,1 or 0,2 or 1,2 which mean 1115 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 1116 where the numbered entries indicate what is in the field. 1117 1118 This function is called once per split (it creates a new split each time). Solve options 1119 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1120 1121 Developer Note: This routine does not actually create the IS representing the split, that is delayed 1122 until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be 1123 available when this routine is called. 1124 1125 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 1126 1127 @*/ 1128 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1129 { 1130 PetscErrorCode ierr; 1131 1132 PetscFunctionBegin; 1133 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1134 PetscValidCharPointer(splitname,2); 1135 if (n < 1) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 1136 PetscValidIntPointer(fields,3); 1137 ierr = PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt *,const PetscInt *),(pc,splitname,n,fields,fields_col));CHKERRQ(ierr); 1138 PetscFunctionReturn(0); 1139 } 1140 1141 #undef __FUNCT__ 1142 #define __FUNCT__ "PCFieldSplitSetIS" 1143 /*@ 1144 PCFieldSplitSetIS - Sets the exact elements for field 1145 1146 Logically Collective on PC 1147 1148 Input Parameters: 1149 + pc - the preconditioner context 1150 . splitname - name of this split, if PETSC_NULL the number of the split is used 1151 - is - the index set that defines the vector elements in this field 1152 1153 1154 Notes: 1155 Use PCFieldSplitSetFields(), for fields defined by strided types. 1156 1157 This function is called once per split (it creates a new split each time). Solve options 1158 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1159 1160 Level: intermediate 1161 1162 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 1163 1164 @*/ 1165 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 1166 { 1167 PetscErrorCode ierr; 1168 1169 PetscFunctionBegin; 1170 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1171 if (splitname) PetscValidCharPointer(splitname,2); 1172 PetscValidHeaderSpecific(is,IS_CLASSID,3); 1173 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr); 1174 PetscFunctionReturn(0); 1175 } 1176 1177 #undef __FUNCT__ 1178 #define __FUNCT__ "PCFieldSplitGetIS" 1179 /*@ 1180 PCFieldSplitGetIS - Retrieves the elements for a field as an IS 1181 1182 Logically Collective on PC 1183 1184 Input Parameters: 1185 + pc - the preconditioner context 1186 - splitname - name of this split 1187 1188 Output Parameter: 1189 - is - the index set that defines the vector elements in this field, or PETSC_NULL if the field is not found 1190 1191 Level: intermediate 1192 1193 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS() 1194 1195 @*/ 1196 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is) 1197 { 1198 PetscErrorCode ierr; 1199 1200 PetscFunctionBegin; 1201 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1202 PetscValidCharPointer(splitname,2); 1203 PetscValidPointer(is,3); 1204 { 1205 PC_FieldSplit *jac = (PC_FieldSplit *) pc->data; 1206 PC_FieldSplitLink ilink = jac->head; 1207 PetscBool found; 1208 1209 *is = PETSC_NULL; 1210 while(ilink) { 1211 ierr = PetscStrcmp(ilink->splitname, splitname, &found);CHKERRQ(ierr); 1212 if (found) { 1213 *is = ilink->is; 1214 break; 1215 } 1216 ilink = ilink->next; 1217 } 1218 } 1219 PetscFunctionReturn(0); 1220 } 1221 1222 #undef __FUNCT__ 1223 #define __FUNCT__ "PCFieldSplitSetBlockSize" 1224 /*@ 1225 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 1226 fieldsplit preconditioner. If not set the matrix block size is used. 1227 1228 Logically Collective on PC 1229 1230 Input Parameters: 1231 + pc - the preconditioner context 1232 - bs - the block size 1233 1234 Level: intermediate 1235 1236 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 1237 1238 @*/ 1239 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 1240 { 1241 PetscErrorCode ierr; 1242 1243 PetscFunctionBegin; 1244 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1245 PetscValidLogicalCollectiveInt(pc,bs,2); 1246 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr); 1247 PetscFunctionReturn(0); 1248 } 1249 1250 #undef __FUNCT__ 1251 #define __FUNCT__ "PCFieldSplitGetSubKSP" 1252 /*@C 1253 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1254 1255 Collective on KSP 1256 1257 Input Parameter: 1258 . pc - the preconditioner context 1259 1260 Output Parameters: 1261 + n - the number of splits 1262 - pc - the array of KSP contexts 1263 1264 Note: 1265 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 1266 (not the KSP just the array that contains them). 1267 1268 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 1269 1270 Level: advanced 1271 1272 .seealso: PCFIELDSPLIT 1273 @*/ 1274 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1275 { 1276 PetscErrorCode ierr; 1277 1278 PetscFunctionBegin; 1279 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1280 if (n) PetscValidIntPointer(n,2); 1281 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 1282 PetscFunctionReturn(0); 1283 } 1284 1285 #undef __FUNCT__ 1286 #define __FUNCT__ "PCFieldSplitSchurPrecondition" 1287 /*@ 1288 PCFieldSplitSchurPrecondition - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the 1289 A11 matrix. Otherwise no preconditioner is used. 1290 1291 Collective on PC 1292 1293 Input Parameters: 1294 + pc - the preconditioner context 1295 . ptype - which matrix to use for preconditioning the Schur complement, PC_FIELDSPLIT_SCHUR_PRE_DIAG (diag) is default 1296 - userpre - matrix to use for preconditioning, or PETSC_NULL 1297 1298 Options Database: 1299 . -pc_fieldsplit_schur_precondition <self,user,diag> default is diag 1300 1301 Notes: 1302 $ If ptype is 1303 $ user then the preconditioner for the Schur complement is generated by the provided matrix (pre argument 1304 $ to this function). 1305 $ diag then the preconditioner for the Schur complement is generated by the block diagonal part of the original 1306 $ matrix associated with the Schur complement (i.e. A11) 1307 $ self the preconditioner for the Schur complement is generated from the Schur complement matrix itself: 1308 $ The only preconditioner that currently works directly with the Schur complement matrix object is the PCLSC 1309 $ preconditioner 1310 1311 When solving a saddle point problem, where the A11 block is identically zero, using diag as the ptype only makes sense 1312 with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and 1313 -fieldsplit_1_pc_type lsc which uses the least squares commutator compute a preconditioner for the Schur complement. 1314 1315 Developer Notes: This is a terrible name, gives no good indication of what the function does and should also have Set in 1316 the name since it sets a proceedure to use. 1317 1318 Level: intermediate 1319 1320 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC 1321 1322 @*/ 1323 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1324 { 1325 PetscErrorCode ierr; 1326 1327 PetscFunctionBegin; 1328 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1329 ierr = PetscTryMethod(pc,"PCFieldSplitSchurPrecondition_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr); 1330 PetscFunctionReturn(0); 1331 } 1332 1333 EXTERN_C_BEGIN 1334 #undef __FUNCT__ 1335 #define __FUNCT__ "PCFieldSplitSchurPrecondition_FieldSplit" 1336 PetscErrorCode PCFieldSplitSchurPrecondition_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1337 { 1338 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1339 PetscErrorCode ierr; 1340 1341 PetscFunctionBegin; 1342 jac->schurpre = ptype; 1343 if (pre) { 1344 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1345 jac->schur_user = pre; 1346 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 1347 } 1348 PetscFunctionReturn(0); 1349 } 1350 EXTERN_C_END 1351 1352 #undef __FUNCT__ 1353 #define __FUNCT__ "PCFieldSplitSetSchurFactType" 1354 /*@ 1355 PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain 1356 1357 Collective on PC 1358 1359 Input Parameters: 1360 + pc - the preconditioner context 1361 - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default 1362 1363 Options Database: 1364 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full 1365 1366 1367 Level: intermediate 1368 1369 Notes: 1370 The FULL factorization is 1371 1372 $ (A B) = (1 0) (A 0) (1 Ainv*B) 1373 $ (C D) (C*Ainv 1) (0 S) (0 1 ) 1374 1375 where S = D - C*Ainv*B. In practice, the full factorization is applied via block triangular solves with the grouping L*(D*U). UPPER uses D*U, LOWER uses L*D, 1376 and DIAG is the diagonal part with the sign of S flipped (because this makes the preconditioner positive definite for many formulations, thus allowing the use of KSPMINRES). 1377 1378 If applied exactly, FULL factorization is a direct solver. The preconditioned operator with LOWER or UPPER has all eigenvalues equal to 1 and minimal polynomial 1379 of degree 2, so KSPGMRES converges in 2 iterations. If the iteration count is very low, consider using KSPFGMRES or KSPGCR which can use one less preconditioner 1380 application in this case. Note that the preconditioned operator may be highly non-normal, so such fast convergence may not be observed in practice. With DIAG, 1381 the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations. 1382 1383 For symmetric problems in which A is positive definite and S is negative definite, DIAG can be used with KSPMINRES. Note that a flexible method like KSPFGMRES 1384 or KSPGCR must be used if the fieldsplit preconditioner is nonlinear (e.g. a few iterations of a Krylov method is used inside a split). 1385 1386 References: 1387 Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000) pp. 1969-1972. 1388 1389 Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001), pp. 1050-1051. 1390 1391 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType 1392 @*/ 1393 PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype) 1394 { 1395 PetscErrorCode ierr; 1396 1397 PetscFunctionBegin; 1398 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1399 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));CHKERRQ(ierr); 1400 PetscFunctionReturn(0); 1401 } 1402 1403 #undef __FUNCT__ 1404 #define __FUNCT__ "PCFieldSplitSetSchurFactType_FieldSplit" 1405 PETSC_EXTERN_C PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype) 1406 { 1407 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1408 1409 PetscFunctionBegin; 1410 jac->schurfactorization = ftype; 1411 PetscFunctionReturn(0); 1412 } 1413 1414 #undef __FUNCT__ 1415 #define __FUNCT__ "PCFieldSplitGetSchurBlocks" 1416 /*@C 1417 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 1418 1419 Collective on KSP 1420 1421 Input Parameter: 1422 . pc - the preconditioner context 1423 1424 Output Parameters: 1425 + A00 - the (0,0) block 1426 . A01 - the (0,1) block 1427 . A10 - the (1,0) block 1428 - A11 - the (1,1) block 1429 1430 Level: advanced 1431 1432 .seealso: PCFIELDSPLIT 1433 @*/ 1434 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 1435 { 1436 PC_FieldSplit *jac = (PC_FieldSplit *) pc->data; 1437 1438 PetscFunctionBegin; 1439 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1440 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 1441 if (A00) *A00 = jac->pmat[0]; 1442 if (A01) *A01 = jac->B; 1443 if (A10) *A10 = jac->C; 1444 if (A11) *A11 = jac->pmat[1]; 1445 PetscFunctionReturn(0); 1446 } 1447 1448 EXTERN_C_BEGIN 1449 #undef __FUNCT__ 1450 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 1451 PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 1452 { 1453 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1454 PetscErrorCode ierr; 1455 1456 PetscFunctionBegin; 1457 jac->type = type; 1458 if (type == PC_COMPOSITE_SCHUR) { 1459 pc->ops->apply = PCApply_FieldSplit_Schur; 1460 pc->ops->view = PCView_FieldSplit_Schur; 1461 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit_Schur",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 1462 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","PCFieldSplitSchurPrecondition_FieldSplit",PCFieldSplitSchurPrecondition_FieldSplit);CHKERRQ(ierr); 1463 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetSchurFactType_C","PCFieldSplitSetSchurFactType_FieldSplit",PCFieldSplitSetSchurFactType_FieldSplit);CHKERRQ(ierr); 1464 1465 } else { 1466 pc->ops->apply = PCApply_FieldSplit; 1467 pc->ops->view = PCView_FieldSplit; 1468 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1469 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","",0);CHKERRQ(ierr); 1470 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetSchurFactType_C","",0);CHKERRQ(ierr); 1471 } 1472 PetscFunctionReturn(0); 1473 } 1474 EXTERN_C_END 1475 1476 EXTERN_C_BEGIN 1477 #undef __FUNCT__ 1478 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 1479 PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 1480 { 1481 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1482 1483 PetscFunctionBegin; 1484 if (bs < 1) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 1485 if (jac->bs > 0 && jac->bs != bs) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot change fieldsplit blocksize from %D to %D after it has been set",jac->bs,bs); 1486 jac->bs = bs; 1487 PetscFunctionReturn(0); 1488 } 1489 EXTERN_C_END 1490 1491 #undef __FUNCT__ 1492 #define __FUNCT__ "PCFieldSplitSetType" 1493 /*@ 1494 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 1495 1496 Collective on PC 1497 1498 Input Parameter: 1499 . pc - the preconditioner context 1500 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1501 1502 Options Database Key: 1503 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 1504 1505 Level: Developer 1506 1507 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1508 1509 .seealso: PCCompositeSetType() 1510 1511 @*/ 1512 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 1513 { 1514 PetscErrorCode ierr; 1515 1516 PetscFunctionBegin; 1517 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1518 ierr = PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));CHKERRQ(ierr); 1519 PetscFunctionReturn(0); 1520 } 1521 1522 /* -------------------------------------------------------------------------------------*/ 1523 /*MC 1524 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 1525 fields or groups of fields. See the users manual section "Solving Block Matrices" for more details. 1526 1527 To set options on the solvers for each block append -fieldsplit_ to all the PC 1528 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 1529 1530 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 1531 and set the options directly on the resulting KSP object 1532 1533 Level: intermediate 1534 1535 Options Database Keys: 1536 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 1537 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 1538 been supplied explicitly by -pc_fieldsplit_%d_fields 1539 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 1540 . -pc_fieldsplit_type <additive,multiplicative,schur,symmetric_multiplicative> 1541 . -pc_fieldsplit_schur_precondition <true,false> default is true 1542 . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver 1543 1544 - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ 1545 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 1546 1547 Notes: 1548 Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 1549 to define a field by an arbitrary collection of entries. 1550 1551 If no fields are set the default is used. The fields are defined by entries strided by bs, 1552 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 1553 if this is not called the block size defaults to the blocksize of the second matrix passed 1554 to KSPSetOperators()/PCSetOperators(). 1555 1556 $ For the Schur complement preconditioner if J = ( A00 A01 ) 1557 $ ( A10 A11 ) 1558 $ the preconditioner using full factorization is 1559 $ ( I -A10 ksp(A00) ) ( inv(A00) 0 ) ( I 0 ) 1560 $ ( 0 I ) ( 0 ksp(S) ) ( -A10 ksp(A00) I ) 1561 where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. The action of 1562 ksp(S) is computed using the KSP solver with prefix -fieldsplit_splitname_ (where splitname was given 1563 in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0, 1564 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 1565 A11 is used to construct a preconditioner for S, use PCFieldSplitSchurPrecondition() to turn on or off this 1566 option. You can use the preconditioner PCLSC to precondition the Schur complement with -fieldsplit_1_pc_type lsc. The 1567 factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above, 1568 diag gives 1569 $ ( inv(A00) 0 ) 1570 $ ( 0 -ksp(S) ) 1571 note that slightly counter intuitively there is a negative in front of the ksp(S) so that the preconditioner is positive definite. The lower factorization is the inverse of 1572 $ ( A00 0 ) 1573 $ ( A10 S ) 1574 where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of 1575 $ ( A00 A01 ) 1576 $ ( 0 S ) 1577 where again the inverses of A00 and S are applied using KSPs. 1578 1579 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 1580 is used automatically for a second block. 1581 1582 The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1. 1583 Generally it should be used with the AIJ format. 1584 1585 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 1586 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT 1587 inside a smoother resulting in "Distributive Smoothers". 1588 1589 Concepts: physics based preconditioners, block preconditioners 1590 1591 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC, 1592 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSchurPrecondition() 1593 M*/ 1594 1595 EXTERN_C_BEGIN 1596 #undef __FUNCT__ 1597 #define __FUNCT__ "PCCreate_FieldSplit" 1598 PetscErrorCode PCCreate_FieldSplit(PC pc) 1599 { 1600 PetscErrorCode ierr; 1601 PC_FieldSplit *jac; 1602 1603 PetscFunctionBegin; 1604 ierr = PetscNewLog(pc,PC_FieldSplit,&jac);CHKERRQ(ierr); 1605 jac->bs = -1; 1606 jac->nsplits = 0; 1607 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 1608 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 1609 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL; 1610 1611 pc->data = (void*)jac; 1612 1613 pc->ops->apply = PCApply_FieldSplit; 1614 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 1615 pc->ops->setup = PCSetUp_FieldSplit; 1616 pc->ops->reset = PCReset_FieldSplit; 1617 pc->ops->destroy = PCDestroy_FieldSplit; 1618 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 1619 pc->ops->view = PCView_FieldSplit; 1620 pc->ops->applyrichardson = 0; 1621 1622 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit", 1623 PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1624 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","PCFieldSplitSetFields_FieldSplit", 1625 PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 1626 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","PCFieldSplitSetIS_FieldSplit", 1627 PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 1628 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","PCFieldSplitSetType_FieldSplit", 1629 PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 1630 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","PCFieldSplitSetBlockSize_FieldSplit", 1631 PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 1632 PetscFunctionReturn(0); 1633 } 1634 EXTERN_C_END 1635 1636 1637 1638