1 #define PETSCKSP_DLL 2 3 #include "private/pcimpl.h" /*I "petscpc.h" I*/ 4 5 const char *const PCFieldSplitSchurPreTypes[] = {"SELF","DIAG","USER","PCFieldSplitSchurPreType","PC_FIELDSPLIT_SCHUR_PRE_",0}; 6 const char *const PCFieldSplitSchurFactorizationTypes[] = {"DIAG","LOWER","UPPER","FULL","PCFieldSplitSchurFactorizationType","PC_FIELDSPLIT_SCHUR_FACTORIZATION_",0}; 7 8 typedef enum { 9 PC_FIELDSPLIT_SCHUR_FACTORIZATION_DIAG, 10 PC_FIELDSPLIT_SCHUR_FACTORIZATION_LOWER, 11 PC_FIELDSPLIT_SCHUR_FACTORIZATION_UPPER, 12 PC_FIELDSPLIT_SCHUR_FACTORIZATION_FULL 13 } PCFieldSplitSchurFactorizationType; 14 15 typedef struct _PC_FieldSplitLink *PC_FieldSplitLink; 16 struct _PC_FieldSplitLink { 17 KSP ksp; 18 Vec x,y; 19 char *splitname; 20 PetscInt nfields; 21 PetscInt *fields; 22 VecScatter sctx; 23 IS is; 24 PC_FieldSplitLink next,previous; 25 }; 26 27 typedef struct { 28 PCCompositeType type; 29 PetscBool defaultsplit; /* Flag for a system with a set of 'k' scalar fields with the same layout (and bs = k) */ 30 PetscBool splitdefined; /* Flag is set after the splits have been defined, to prevent more splits from being added */ 31 PetscBool realdiagonal; /* Flag to use the diagonal blocks of mat preconditioned by pmat, instead of just pmat */ 32 PetscInt bs; /* Block size for IS and Mat structures */ 33 PetscInt nsplits; /* Number of field divisions defined */ 34 Vec *x,*y,w1,w2; 35 Mat *mat; /* The diagonal block for each split */ 36 Mat *pmat; /* The preconditioning diagonal block for each split */ 37 Mat *Afield; /* The rows of the matrix associated with each split */ 38 PetscBool issetup; 39 /* Only used when Schur complement preconditioning is used */ 40 Mat B; /* The (0,1) block */ 41 Mat C; /* The (1,0) block */ 42 Mat schur; /* The Schur complement S = D - C A^{-1} B */ 43 Mat schur_user; /* User-provided preconditioning matrix for the Schur complement */ 44 PCFieldSplitSchurPreType schurpre; /* Determines which preconditioning matrix is used for the Schur complement */ 45 PCFieldSplitSchurFactorizationType schurfactorization; 46 KSP kspschur; /* The solver for S */ 47 PC_FieldSplitLink head; 48 } PC_FieldSplit; 49 50 /* 51 Notes: there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of 52 inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the 53 PC you could change this. 54 */ 55 56 /* This helper is so that setting a user-provided preconditioning matrix is orthogonal to choosing to use it. This way the 57 * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */ 58 static Mat FieldSplitSchurPre(PC_FieldSplit *jac) 59 { 60 switch (jac->schurpre) { 61 case PC_FIELDSPLIT_SCHUR_PRE_SELF: return jac->schur; 62 case PC_FIELDSPLIT_SCHUR_PRE_DIAG: return jac->pmat[1]; 63 case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */ 64 default: 65 return jac->schur_user ? jac->schur_user : jac->pmat[1]; 66 } 67 } 68 69 70 #undef __FUNCT__ 71 #define __FUNCT__ "PCView_FieldSplit" 72 static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer) 73 { 74 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 75 PetscErrorCode ierr; 76 PetscBool iascii; 77 PetscInt i,j; 78 PC_FieldSplitLink ilink = jac->head; 79 80 PetscFunctionBegin; 81 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 82 if (iascii) { 83 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);CHKERRQ(ierr); 84 ierr = PetscViewerASCIIPrintf(viewer," Solver info for each split is in the following KSP objects:\n");CHKERRQ(ierr); 85 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 86 for (i=0; i<jac->nsplits; i++) { 87 if (ilink->fields) { 88 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr); 89 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 90 for (j=0; j<ilink->nfields; j++) { 91 if (j > 0) { 92 ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr); 93 } 94 ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr); 95 } 96 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 97 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 98 } else { 99 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);CHKERRQ(ierr); 100 } 101 ierr = KSPView(ilink->ksp,viewer);CHKERRQ(ierr); 102 ilink = ilink->next; 103 } 104 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 105 } else { 106 SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for PCFieldSplit",((PetscObject)viewer)->type_name); 107 } 108 PetscFunctionReturn(0); 109 } 110 111 #undef __FUNCT__ 112 #define __FUNCT__ "PCView_FieldSplit_Schur" 113 static PetscErrorCode PCView_FieldSplit_Schur(PC pc,PetscViewer viewer) 114 { 115 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 116 PetscErrorCode ierr; 117 PetscBool iascii; 118 PetscInt i,j; 119 PC_FieldSplitLink ilink = jac->head; 120 KSP ksp; 121 122 PetscFunctionBegin; 123 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 124 if (iascii) { 125 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, blocksize = %D, factorization %s\n",jac->bs,PCFieldSplitSchurFactorizationTypes[jac->schurfactorization]);CHKERRQ(ierr); 126 ierr = PetscViewerASCIIPrintf(viewer," Split info:\n");CHKERRQ(ierr); 127 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 128 for (i=0; i<jac->nsplits; i++) { 129 if (ilink->fields) { 130 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr); 131 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 132 for (j=0; j<ilink->nfields; j++) { 133 if (j > 0) { 134 ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr); 135 } 136 ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr); 137 } 138 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 139 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 140 } else { 141 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);CHKERRQ(ierr); 142 } 143 ilink = ilink->next; 144 } 145 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for A block \n");CHKERRQ(ierr); 146 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 147 if (jac->schur) { 148 ierr = MatSchurComplementGetKSP(jac->schur,&ksp);CHKERRQ(ierr); 149 ierr = KSPView(ksp,viewer);CHKERRQ(ierr); 150 } else { 151 ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr); 152 } 153 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 154 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for S = D - C inv(A) B \n");CHKERRQ(ierr); 155 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 156 if (jac->kspschur) { 157 ierr = KSPView(jac->kspschur,viewer);CHKERRQ(ierr); 158 } else { 159 ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr); 160 } 161 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 162 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 163 } else { 164 SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for PCFieldSplit",((PetscObject)viewer)->type_name); 165 } 166 PetscFunctionReturn(0); 167 } 168 169 #undef __FUNCT__ 170 #define __FUNCT__ "PCFieldSplitSetRuntimeSplits_Private" 171 /* Precondition: jac->bs is set to a meaningful value */ 172 static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc) 173 { 174 PetscErrorCode ierr; 175 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 176 PetscInt i,nfields,*ifields; 177 PetscBool flg; 178 char optionname[128],splitname[8]; 179 180 PetscFunctionBegin; 181 ierr = PetscMalloc(jac->bs*sizeof(PetscInt),&ifields);CHKERRQ(ierr); 182 for (i=0,flg=PETSC_TRUE; ; i++) { 183 ierr = PetscSNPrintf(splitname,sizeof splitname,"%D",i);CHKERRQ(ierr); 184 ierr = PetscSNPrintf(optionname,sizeof optionname,"-pc_fieldsplit_%D_fields",i);CHKERRQ(ierr); 185 nfields = jac->bs; 186 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->prefix,optionname,ifields,&nfields,&flg);CHKERRQ(ierr); 187 if (!flg) break; 188 if (!nfields) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields"); 189 ierr = PCFieldSplitSetFields(pc,splitname,nfields,ifields);CHKERRQ(ierr); 190 } 191 if (i > 0) { 192 /* Makes command-line setting of splits take precedence over setting them in code. 193 Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would 194 create new splits, which would probably not be what the user wanted. */ 195 jac->splitdefined = PETSC_TRUE; 196 } 197 ierr = PetscFree(ifields);CHKERRQ(ierr); 198 PetscFunctionReturn(0); 199 } 200 201 #undef __FUNCT__ 202 #define __FUNCT__ "PCFieldSplitSetDefaults" 203 static PetscErrorCode PCFieldSplitSetDefaults(PC pc) 204 { 205 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 206 PetscErrorCode ierr; 207 PC_FieldSplitLink ilink = jac->head; 208 PetscBool flg = PETSC_FALSE,stokes = PETSC_FALSE; 209 PetscInt i; 210 211 PetscFunctionBegin; 212 if (!ilink) { 213 if (pc->dm) { 214 PetscBool dmcomposite; 215 ierr = PetscTypeCompare((PetscObject)pc->dm,DMCOMPOSITE,&dmcomposite);CHKERRQ(ierr); 216 if (dmcomposite) { 217 PetscInt nDM; 218 IS *fields; 219 ierr = PetscInfo(pc,"Setting up physics based fieldsplit preconditioner using the embedded DM\n");CHKERRQ(ierr); 220 ierr = DMCompositeGetNumberDM(pc->dm,&nDM);CHKERRQ(ierr); 221 ierr = DMCompositeGetGlobalISs(pc->dm,&fields);CHKERRQ(ierr); 222 for (i=0; i<nDM; i++) { 223 char splitname[8]; 224 ierr = PetscSNPrintf(splitname,sizeof splitname,"%D",i);CHKERRQ(ierr); 225 ierr = PCFieldSplitSetIS(pc,splitname,fields[i]);CHKERRQ(ierr); 226 ierr = ISDestroy(fields[i]);CHKERRQ(ierr); 227 } 228 ierr = PetscFree(fields);CHKERRQ(ierr); 229 } 230 } else { 231 if (jac->bs <= 0) { 232 if (pc->pmat) { 233 ierr = MatGetBlockSize(pc->pmat,&jac->bs);CHKERRQ(ierr); 234 } else { 235 jac->bs = 1; 236 } 237 } 238 239 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&flg,PETSC_NULL);CHKERRQ(ierr); 240 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_stokes",&stokes,PETSC_NULL);CHKERRQ(ierr); 241 if (stokes) { 242 IS zerodiags,rest; 243 PetscInt nmin,nmax; 244 245 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 246 ierr = MatFindZeroDiagonals(pc->mat,&zerodiags);CHKERRQ(ierr); 247 ierr = ISComplement(zerodiags,nmin,nmax,&rest);CHKERRQ(ierr); 248 ierr = PCFieldSplitSetIS(pc,"0",rest);CHKERRQ(ierr); 249 ierr = PCFieldSplitSetIS(pc,"1",zerodiags);CHKERRQ(ierr); 250 ierr = ISDestroy(zerodiags);CHKERRQ(ierr); 251 ierr = ISDestroy(rest);CHKERRQ(ierr); 252 } else { 253 if (!flg) { 254 /* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit() 255 then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */ 256 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 257 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 258 } 259 if (flg || !jac->splitdefined) { 260 ierr = PetscInfo(pc,"Using default splitting of fields\n");CHKERRQ(ierr); 261 for (i=0; i<jac->bs; i++) { 262 char splitname[8]; 263 ierr = PetscSNPrintf(splitname,sizeof splitname,"%D",i);CHKERRQ(ierr); 264 ierr = PCFieldSplitSetFields(pc,splitname,1,&i);CHKERRQ(ierr); 265 } 266 jac->defaultsplit = PETSC_TRUE; 267 } 268 } 269 } 270 } else if (jac->nsplits == 1) { 271 if (ilink->is) { 272 IS is2; 273 PetscInt nmin,nmax; 274 275 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 276 ierr = ISComplement(ilink->is,nmin,nmax,&is2);CHKERRQ(ierr); 277 ierr = PCFieldSplitSetIS(pc,"1",is2);CHKERRQ(ierr); 278 ierr = ISDestroy(is2);CHKERRQ(ierr); 279 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Must provide at least two sets of fields to PCFieldSplit()"); 280 } 281 if (jac->nsplits < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unhandled case, must have at least two fields"); 282 PetscFunctionReturn(0); 283 } 284 285 286 #undef __FUNCT__ 287 #define __FUNCT__ "PCSetUp_FieldSplit" 288 static PetscErrorCode PCSetUp_FieldSplit(PC pc) 289 { 290 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 291 PetscErrorCode ierr; 292 PC_FieldSplitLink ilink; 293 PetscInt i,nsplit,ccsize; 294 MatStructure flag = pc->flag; 295 PetscBool sorted; 296 297 PetscFunctionBegin; 298 ierr = PCFieldSplitSetDefaults(pc);CHKERRQ(ierr); 299 nsplit = jac->nsplits; 300 ilink = jac->head; 301 302 /* get the matrices for each split */ 303 if (!jac->issetup) { 304 PetscInt rstart,rend,nslots,bs; 305 306 jac->issetup = PETSC_TRUE; 307 308 /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */ 309 bs = jac->bs; 310 ierr = MatGetOwnershipRange(pc->pmat,&rstart,&rend);CHKERRQ(ierr); 311 ierr = MatGetLocalSize(pc->pmat,PETSC_NULL,&ccsize);CHKERRQ(ierr); 312 nslots = (rend - rstart)/bs; 313 for (i=0; i<nsplit; i++) { 314 if (jac->defaultsplit) { 315 ierr = ISCreateStride(((PetscObject)pc)->comm,nslots,rstart+i,nsplit,&ilink->is);CHKERRQ(ierr); 316 } else if (!ilink->is) { 317 if (ilink->nfields > 1) { 318 PetscInt *ii,j,k,nfields = ilink->nfields,*fields = ilink->fields; 319 ierr = PetscMalloc(ilink->nfields*nslots*sizeof(PetscInt),&ii);CHKERRQ(ierr); 320 for (j=0; j<nslots; j++) { 321 for (k=0; k<nfields; k++) { 322 ii[nfields*j + k] = rstart + bs*j + fields[k]; 323 } 324 } 325 ierr = ISCreateGeneral(((PetscObject)pc)->comm,nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is);CHKERRQ(ierr); 326 } else { 327 ierr = ISCreateStride(((PetscObject)pc)->comm,nslots,rstart+ilink->fields[0],bs,&ilink->is);CHKERRQ(ierr); 328 } 329 } 330 ierr = ISSorted(ilink->is,&sorted);CHKERRQ(ierr); 331 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Fields must be sorted when creating split"); 332 ilink = ilink->next; 333 } 334 } 335 336 ilink = jac->head; 337 if (!jac->pmat) { 338 ierr = PetscMalloc(nsplit*sizeof(Mat),&jac->pmat);CHKERRQ(ierr); 339 for (i=0; i<nsplit; i++) { 340 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ilink->is,MAT_INITIAL_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 341 ilink = ilink->next; 342 } 343 } else { 344 for (i=0; i<nsplit; i++) { 345 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ilink->is,MAT_REUSE_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 346 ilink = ilink->next; 347 } 348 } 349 if (jac->realdiagonal) { 350 ilink = jac->head; 351 if (!jac->mat) { 352 ierr = PetscMalloc(nsplit*sizeof(Mat),&jac->mat);CHKERRQ(ierr); 353 for (i=0; i<nsplit; i++) { 354 ierr = MatGetSubMatrix(pc->mat,ilink->is,ilink->is,MAT_INITIAL_MATRIX,&jac->mat[i]);CHKERRQ(ierr); 355 ilink = ilink->next; 356 } 357 } else { 358 for (i=0; i<nsplit; i++) { 359 ierr = MatGetSubMatrix(pc->mat,ilink->is,ilink->is,MAT_REUSE_MATRIX,&jac->mat[i]);CHKERRQ(ierr); 360 ilink = ilink->next; 361 } 362 } 363 } else { 364 jac->mat = jac->pmat; 365 } 366 367 if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR) { 368 /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */ 369 ilink = jac->head; 370 if (!jac->Afield) { 371 ierr = PetscMalloc(nsplit*sizeof(Mat),&jac->Afield);CHKERRQ(ierr); 372 for (i=0; i<nsplit; i++) { 373 ierr = MatGetSubMatrix(pc->mat,ilink->is,PETSC_NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 374 ilink = ilink->next; 375 } 376 } else { 377 for (i=0; i<nsplit; i++) { 378 ierr = MatGetSubMatrix(pc->mat,ilink->is,PETSC_NULL,MAT_REUSE_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 379 ilink = ilink->next; 380 } 381 } 382 } 383 384 if (jac->type == PC_COMPOSITE_SCHUR) { 385 IS ccis; 386 PetscInt rstart,rend; 387 if (nsplit != 2) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields"); 388 389 /* When extracting off-diagonal submatrices, we take complements from this range */ 390 ierr = MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);CHKERRQ(ierr); 391 392 /* need to handle case when one is resetting up the preconditioner */ 393 if (jac->schur) { 394 ilink = jac->head; 395 ierr = ISComplement(ilink->is,rstart,rend,&ccis);CHKERRQ(ierr); 396 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);CHKERRQ(ierr); 397 ierr = ISDestroy(ccis);CHKERRQ(ierr); 398 ilink = ilink->next; 399 ierr = ISComplement(ilink->is,rstart,rend,&ccis);CHKERRQ(ierr); 400 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);CHKERRQ(ierr); 401 ierr = ISDestroy(ccis);CHKERRQ(ierr); 402 ierr = MatSchurComplementUpdate(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->pmat[1],pc->flag);CHKERRQ(ierr); 403 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),pc->flag);CHKERRQ(ierr); 404 405 } else { 406 KSP ksp; 407 char schurprefix[256]; 408 409 /* extract the B and C matrices */ 410 ilink = jac->head; 411 ierr = ISComplement(ilink->is,rstart,rend,&ccis);CHKERRQ(ierr); 412 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr); 413 ierr = ISDestroy(ccis);CHKERRQ(ierr); 414 ilink = ilink->next; 415 ierr = ISComplement(ilink->is,rstart,rend,&ccis);CHKERRQ(ierr); 416 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 417 ierr = ISDestroy(ccis);CHKERRQ(ierr); 418 /* Better would be to use 'mat[0]' (diagonal block of the real matrix) preconditioned by pmat[0] */ 419 ierr = MatCreateSchurComplement(jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1],&jac->schur);CHKERRQ(ierr); 420 ierr = MatSchurComplementGetKSP(jac->schur,&ksp);CHKERRQ(ierr); 421 ierr = PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,2);CHKERRQ(ierr); 422 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr); 423 424 ierr = KSPCreate(((PetscObject)pc)->comm,&jac->kspschur);CHKERRQ(ierr); 425 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);CHKERRQ(ierr); 426 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);CHKERRQ(ierr); 427 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 428 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) { 429 PC pc; 430 ierr = KSPGetPC(jac->kspschur,&pc);CHKERRQ(ierr); 431 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 432 /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */ 433 } 434 ierr = PetscSNPrintf(schurprefix,sizeof schurprefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 435 ierr = KSPSetOptionsPrefix(jac->kspschur,schurprefix);CHKERRQ(ierr); 436 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 437 ierr = KSPSetFromOptions(jac->kspschur);CHKERRQ(ierr); 438 439 ierr = PetscMalloc2(2,Vec,&jac->x,2,Vec,&jac->y);CHKERRQ(ierr); 440 ierr = MatGetVecs(jac->pmat[0],&jac->x[0],&jac->y[0]);CHKERRQ(ierr); 441 ierr = MatGetVecs(jac->pmat[1],&jac->x[1],&jac->y[1]);CHKERRQ(ierr); 442 ilink = jac->head; 443 ilink->x = jac->x[0]; ilink->y = jac->y[0]; 444 ilink = ilink->next; 445 ilink->x = jac->x[1]; ilink->y = jac->y[1]; 446 } 447 } else { 448 /* set up the individual PCs */ 449 i = 0; 450 ilink = jac->head; 451 while (ilink) { 452 ierr = KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i],flag);CHKERRQ(ierr); 453 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 454 ierr = KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr); 455 ierr = KSPSetUp(ilink->ksp);CHKERRQ(ierr); 456 i++; 457 ilink = ilink->next; 458 } 459 460 /* create work vectors for each split */ 461 if (!jac->x) { 462 ierr = PetscMalloc2(nsplit,Vec,&jac->x,nsplit,Vec,&jac->y);CHKERRQ(ierr); 463 ilink = jac->head; 464 for (i=0; i<nsplit; i++) { 465 Vec *vl,*vr; 466 467 ierr = KSPGetVecs(ilink->ksp,1,&vr,1,&vl);CHKERRQ(ierr); 468 ilink->x = *vr; 469 ilink->y = *vl; 470 ierr = PetscFree(vr);CHKERRQ(ierr); 471 ierr = PetscFree(vl);CHKERRQ(ierr); 472 jac->x[i] = ilink->x; 473 jac->y[i] = ilink->y; 474 ilink = ilink->next; 475 } 476 } 477 } 478 479 480 if (!jac->head->sctx) { 481 Vec xtmp; 482 483 /* compute scatter contexts needed by multiplicative versions and non-default splits */ 484 485 ilink = jac->head; 486 ierr = MatGetVecs(pc->pmat,&xtmp,PETSC_NULL);CHKERRQ(ierr); 487 for (i=0; i<nsplit; i++) { 488 ierr = VecScatterCreate(xtmp,ilink->is,jac->x[i],PETSC_NULL,&ilink->sctx);CHKERRQ(ierr); 489 ilink = ilink->next; 490 } 491 ierr = VecDestroy(xtmp);CHKERRQ(ierr); 492 } 493 PetscFunctionReturn(0); 494 } 495 496 #define FieldSplitSplitSolveAdd(ilink,xx,yy) \ 497 (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 498 VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 499 KSPSolve(ilink->ksp,ilink->x,ilink->y) || \ 500 VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \ 501 VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE)) 502 503 #undef __FUNCT__ 504 #define __FUNCT__ "PCApply_FieldSplit_Schur" 505 static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y) 506 { 507 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 508 PetscErrorCode ierr; 509 KSP ksp; 510 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 511 512 PetscFunctionBegin; 513 ierr = MatSchurComplementGetKSP(jac->schur,&ksp);CHKERRQ(ierr); 514 515 switch (jac->schurfactorization) { 516 case PC_FIELDSPLIT_SCHUR_FACTORIZATION_DIAG: 517 /* [A 0; 0 -S], positive definite, suitable for MINRES */ 518 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 519 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 520 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 521 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 522 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 523 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 524 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 525 ierr = VecScale(ilinkD->y,-1.);CHKERRQ(ierr); 526 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 527 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 528 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 529 break; 530 case PC_FIELDSPLIT_SCHUR_FACTORIZATION_LOWER: 531 /* [A 0; C S], suitable for left preconditioning */ 532 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 533 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 534 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 535 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 536 ierr = VecScale(ilinkD->x,-1.);CHKERRQ(ierr); 537 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 538 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 539 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 540 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 541 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 542 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 543 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 544 break; 545 case PC_FIELDSPLIT_SCHUR_FACTORIZATION_UPPER: 546 /* [A B; 0 S], suitable for right preconditioning */ 547 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 548 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 549 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 550 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 551 ierr = VecScale(ilinkA->x,-1.);CHKERRQ(ierr); 552 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 553 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 554 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 555 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 556 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 557 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 558 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 559 break; 560 case PC_FIELDSPLIT_SCHUR_FACTORIZATION_FULL: 561 /* [1 0; CA^{-1} 1] [A 0; 0 S] [1 A^{-1}B; 0 1], an exact solve if applied exactly, needs one extra solve with A */ 562 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 563 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 564 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 565 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 566 ierr = VecScale(ilinkD->x,-1.0);CHKERRQ(ierr); 567 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 568 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 569 570 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 571 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 572 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 573 574 ierr = MatMult(jac->B,ilinkD->y,ilinkA->y);CHKERRQ(ierr); 575 ierr = VecAXPY(ilinkA->x,-1.0,ilinkA->y);CHKERRQ(ierr); 576 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 577 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 578 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 579 } 580 PetscFunctionReturn(0); 581 } 582 583 #undef __FUNCT__ 584 #define __FUNCT__ "PCApply_FieldSplit" 585 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y) 586 { 587 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 588 PetscErrorCode ierr; 589 PC_FieldSplitLink ilink = jac->head; 590 PetscInt cnt; 591 592 PetscFunctionBegin; 593 CHKMEMQ; 594 ierr = VecSetBlockSize(x,jac->bs);CHKERRQ(ierr); 595 ierr = VecSetBlockSize(y,jac->bs);CHKERRQ(ierr); 596 597 if (jac->type == PC_COMPOSITE_ADDITIVE) { 598 if (jac->defaultsplit) { 599 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 600 while (ilink) { 601 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 602 ilink = ilink->next; 603 } 604 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 605 } else { 606 ierr = VecSet(y,0.0);CHKERRQ(ierr); 607 while (ilink) { 608 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 609 ilink = ilink->next; 610 } 611 } 612 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 613 if (!jac->w1) { 614 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 615 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 616 } 617 ierr = VecSet(y,0.0);CHKERRQ(ierr); 618 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 619 cnt = 1; 620 while (ilink->next) { 621 ilink = ilink->next; 622 /* compute the residual only over the part of the vector needed */ 623 ierr = MatMult(jac->Afield[cnt++],y,ilink->x);CHKERRQ(ierr); 624 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 625 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 626 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 627 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 628 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 629 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 630 } 631 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 632 cnt -= 2; 633 while (ilink->previous) { 634 ilink = ilink->previous; 635 /* compute the residual only over the part of the vector needed */ 636 ierr = MatMult(jac->Afield[cnt--],y,ilink->x);CHKERRQ(ierr); 637 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 638 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 639 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 640 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 641 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 642 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 643 } 644 } 645 } else SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type); 646 CHKMEMQ; 647 PetscFunctionReturn(0); 648 } 649 650 #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \ 651 (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 652 VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 653 KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \ 654 VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \ 655 VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE)) 656 657 #undef __FUNCT__ 658 #define __FUNCT__ "PCApplyTranspose_FieldSplit" 659 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y) 660 { 661 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 662 PetscErrorCode ierr; 663 PC_FieldSplitLink ilink = jac->head; 664 665 PetscFunctionBegin; 666 CHKMEMQ; 667 ierr = VecSetBlockSize(x,jac->bs);CHKERRQ(ierr); 668 ierr = VecSetBlockSize(y,jac->bs);CHKERRQ(ierr); 669 670 if (jac->type == PC_COMPOSITE_ADDITIVE) { 671 if (jac->defaultsplit) { 672 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 673 while (ilink) { 674 ierr = KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 675 ilink = ilink->next; 676 } 677 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 678 } else { 679 ierr = VecSet(y,0.0);CHKERRQ(ierr); 680 while (ilink) { 681 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 682 ilink = ilink->next; 683 } 684 } 685 } else { 686 if (!jac->w1) { 687 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 688 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 689 } 690 ierr = VecSet(y,0.0);CHKERRQ(ierr); 691 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 692 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 693 while (ilink->next) { 694 ilink = ilink->next; 695 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 696 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 697 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 698 } 699 while (ilink->previous) { 700 ilink = ilink->previous; 701 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 702 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 703 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 704 } 705 } else { 706 while (ilink->next) { /* get to last entry in linked list */ 707 ilink = ilink->next; 708 } 709 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 710 while (ilink->previous) { 711 ilink = ilink->previous; 712 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 713 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 714 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 715 } 716 } 717 } 718 CHKMEMQ; 719 PetscFunctionReturn(0); 720 } 721 722 #undef __FUNCT__ 723 #define __FUNCT__ "PCDestroy_FieldSplit" 724 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 725 { 726 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 727 PetscErrorCode ierr; 728 PC_FieldSplitLink ilink = jac->head,next; 729 730 PetscFunctionBegin; 731 while (ilink) { 732 ierr = KSPDestroy(ilink->ksp);CHKERRQ(ierr); 733 if (ilink->x) {ierr = VecDestroy(ilink->x);CHKERRQ(ierr);} 734 if (ilink->y) {ierr = VecDestroy(ilink->y);CHKERRQ(ierr);} 735 if (ilink->sctx) {ierr = VecScatterDestroy(ilink->sctx);CHKERRQ(ierr);} 736 if (ilink->is) {ierr = ISDestroy(ilink->is);CHKERRQ(ierr);} 737 next = ilink->next; 738 ierr = PetscFree(ilink->splitname);CHKERRQ(ierr); 739 ierr = PetscFree(ilink->fields);CHKERRQ(ierr); 740 ierr = PetscFree(ilink);CHKERRQ(ierr); 741 ilink = next; 742 } 743 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 744 if (jac->mat && jac->mat != jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->mat);CHKERRQ(ierr);} 745 if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);} 746 if (jac->Afield) {ierr = MatDestroyMatrices(jac->nsplits,&jac->Afield);CHKERRQ(ierr);} 747 if (jac->w1) {ierr = VecDestroy(jac->w1);CHKERRQ(ierr);} 748 if (jac->w2) {ierr = VecDestroy(jac->w2);CHKERRQ(ierr);} 749 if (jac->schur) {ierr = MatDestroy(jac->schur);CHKERRQ(ierr);} 750 if (jac->schur_user) {ierr = MatDestroy(jac->schur_user);CHKERRQ(ierr);} 751 if (jac->kspschur) {ierr = KSPDestroy(jac->kspschur);CHKERRQ(ierr);} 752 if (jac->B) {ierr = MatDestroy(jac->B);CHKERRQ(ierr);} 753 if (jac->C) {ierr = MatDestroy(jac->C);CHKERRQ(ierr);} 754 ierr = PetscFree(jac);CHKERRQ(ierr); 755 PetscFunctionReturn(0); 756 } 757 758 #undef __FUNCT__ 759 #define __FUNCT__ "PCSetFromOptions_FieldSplit" 760 static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc) 761 { 762 PetscErrorCode ierr; 763 PetscInt bs; 764 PetscBool flg,stokes; 765 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 766 PCCompositeType ctype; 767 768 PetscFunctionBegin; 769 ierr = PetscOptionsHead("FieldSplit options");CHKERRQ(ierr); 770 ierr = PetscOptionsBool("-pc_fieldsplit_real_diagonal","Use diagonal blocks of the operator","PCFieldSplitSetRealDiagonal",jac->realdiagonal,&jac->realdiagonal,PETSC_NULL);CHKERRQ(ierr); 771 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 772 if (flg) { 773 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 774 } 775 776 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_stokes",&stokes,PETSC_NULL);CHKERRQ(ierr); 777 if (stokes) { 778 ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);CHKERRQ(ierr); 779 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF; 780 } 781 782 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr); 783 if (flg) { 784 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr); 785 } 786 787 /* Only setup fields once */ 788 if ((jac->bs > 0) && (jac->nsplits == 0)) { 789 /* only allow user to set fields from command line if bs is already known. 790 otherwise user can set them in PCFieldSplitSetDefaults() */ 791 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 792 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 793 } 794 if (jac->type == PC_COMPOSITE_SCHUR) { 795 ierr = PetscOptionsEnum("-pc_fieldsplit_schur_factorization_type","Factorization to use","None",PCFieldSplitSchurFactorizationTypes,(PetscEnum)jac->schurfactorization,(PetscEnum*)&jac->schurfactorization,PETSC_NULL);CHKERRQ(ierr); 796 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); 797 } 798 ierr = PetscOptionsTail();CHKERRQ(ierr); 799 PetscFunctionReturn(0); 800 } 801 802 /*------------------------------------------------------------------------------------*/ 803 804 EXTERN_C_BEGIN 805 #undef __FUNCT__ 806 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 807 PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields) 808 { 809 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 810 PetscErrorCode ierr; 811 PC_FieldSplitLink ilink,next = jac->head; 812 char prefix[128]; 813 PetscInt i; 814 815 PetscFunctionBegin; 816 if (jac->splitdefined) { 817 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 818 PetscFunctionReturn(0); 819 } 820 for (i=0; i<n; i++) { 821 if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 822 if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 823 } 824 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 825 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 826 ierr = PetscMalloc(n*sizeof(PetscInt),&ilink->fields);CHKERRQ(ierr); 827 ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 828 ilink->nfields = n; 829 ilink->next = PETSC_NULL; 830 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 831 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 832 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 833 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 834 835 ierr = PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",splitname);CHKERRQ(ierr); 836 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 837 838 if (!next) { 839 jac->head = ilink; 840 ilink->previous = PETSC_NULL; 841 } else { 842 while (next->next) { 843 next = next->next; 844 } 845 next->next = ilink; 846 ilink->previous = next; 847 } 848 jac->nsplits++; 849 PetscFunctionReturn(0); 850 } 851 EXTERN_C_END 852 853 EXTERN_C_BEGIN 854 #undef __FUNCT__ 855 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit_Schur" 856 PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 857 { 858 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 859 PetscErrorCode ierr; 860 861 PetscFunctionBegin; 862 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 863 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr); 864 (*subksp)[1] = jac->kspschur; 865 *n = jac->nsplits; 866 PetscFunctionReturn(0); 867 } 868 EXTERN_C_END 869 870 EXTERN_C_BEGIN 871 #undef __FUNCT__ 872 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 873 PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 874 { 875 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 876 PetscErrorCode ierr; 877 PetscInt cnt = 0; 878 PC_FieldSplitLink ilink = jac->head; 879 880 PetscFunctionBegin; 881 ierr = PetscMalloc(jac->nsplits*sizeof(KSP*),subksp);CHKERRQ(ierr); 882 while (ilink) { 883 (*subksp)[cnt++] = ilink->ksp; 884 ilink = ilink->next; 885 } 886 if (cnt != jac->nsplits) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt PCFIELDSPLIT object: number splits in linked list %D in object %D",cnt,jac->nsplits); 887 *n = jac->nsplits; 888 PetscFunctionReturn(0); 889 } 890 EXTERN_C_END 891 892 EXTERN_C_BEGIN 893 #undef __FUNCT__ 894 #define __FUNCT__ "PCFieldSplitSetIS_FieldSplit" 895 PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 896 { 897 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 898 PetscErrorCode ierr; 899 PC_FieldSplitLink ilink, next = jac->head; 900 char prefix[128]; 901 902 PetscFunctionBegin; 903 if (jac->splitdefined) { 904 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 905 PetscFunctionReturn(0); 906 } 907 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 908 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 909 ilink->is = is; 910 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 911 ilink->next = PETSC_NULL; 912 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 913 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 914 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 915 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 916 917 ierr = PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",splitname);CHKERRQ(ierr); 918 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 919 920 if (!next) { 921 jac->head = ilink; 922 ilink->previous = PETSC_NULL; 923 } else { 924 while (next->next) { 925 next = next->next; 926 } 927 next->next = ilink; 928 ilink->previous = next; 929 } 930 jac->nsplits++; 931 932 PetscFunctionReturn(0); 933 } 934 EXTERN_C_END 935 936 #undef __FUNCT__ 937 #define __FUNCT__ "PCFieldSplitSetFields" 938 /*@ 939 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 940 941 Logically Collective on PC 942 943 Input Parameters: 944 + pc - the preconditioner context 945 . splitname - name of this split 946 . n - the number of fields in this split 947 - fields - the fields in this split 948 949 Level: intermediate 950 951 Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 952 953 The PCFieldSplitSetFields() is for defining fields as a strided blocks. For example, if the block 954 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 955 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 956 where the numbered entries indicate what is in the field. 957 958 This function is called once per split (it creates a new split each time). Solve options 959 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 960 961 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 962 963 @*/ 964 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields) 965 { 966 PetscErrorCode ierr; 967 968 PetscFunctionBegin; 969 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 970 PetscValidCharPointer(splitname,2); 971 if (n < 1) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 972 PetscValidIntPointer(fields,3); 973 ierr = PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt *),(pc,splitname,n,fields));CHKERRQ(ierr); 974 PetscFunctionReturn(0); 975 } 976 977 #undef __FUNCT__ 978 #define __FUNCT__ "PCFieldSplitSetIS" 979 /*@ 980 PCFieldSplitSetIS - Sets the exact elements for field 981 982 Logically Collective on PC 983 984 Input Parameters: 985 + pc - the preconditioner context 986 . splitname - name of this split 987 - is - the index set that defines the vector elements in this field 988 989 990 Notes: 991 Use PCFieldSplitSetFields(), for fields defined by strided types. 992 993 This function is called once per split (it creates a new split each time). Solve options 994 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 995 996 Level: intermediate 997 998 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 999 1000 @*/ 1001 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 1002 { 1003 PetscErrorCode ierr; 1004 1005 PetscFunctionBegin; 1006 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1007 PetscValidCharPointer(splitname,2); 1008 PetscValidHeaderSpecific(is,IS_CLASSID,3); 1009 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr); 1010 PetscFunctionReturn(0); 1011 } 1012 1013 #undef __FUNCT__ 1014 #define __FUNCT__ "PCFieldSplitSetBlockSize" 1015 /*@ 1016 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 1017 fieldsplit preconditioner. If not set the matrix block size is used. 1018 1019 Logically Collective on PC 1020 1021 Input Parameters: 1022 + pc - the preconditioner context 1023 - bs - the block size 1024 1025 Level: intermediate 1026 1027 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 1028 1029 @*/ 1030 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 1031 { 1032 PetscErrorCode ierr; 1033 1034 PetscFunctionBegin; 1035 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1036 PetscValidLogicalCollectiveInt(pc,bs,2); 1037 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr); 1038 PetscFunctionReturn(0); 1039 } 1040 1041 #undef __FUNCT__ 1042 #define __FUNCT__ "PCFieldSplitGetSubKSP" 1043 /*@C 1044 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1045 1046 Collective on KSP 1047 1048 Input Parameter: 1049 . pc - the preconditioner context 1050 1051 Output Parameters: 1052 + n - the number of split 1053 - pc - the array of KSP contexts 1054 1055 Note: 1056 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 1057 (not the KSP just the array that contains them). 1058 1059 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 1060 1061 Level: advanced 1062 1063 .seealso: PCFIELDSPLIT 1064 @*/ 1065 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1066 { 1067 PetscErrorCode ierr; 1068 1069 PetscFunctionBegin; 1070 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1071 PetscValidIntPointer(n,2); 1072 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 1073 PetscFunctionReturn(0); 1074 } 1075 1076 #undef __FUNCT__ 1077 #define __FUNCT__ "PCFieldSplitSchurPrecondition" 1078 /*@ 1079 PCFieldSplitSchurPrecondition - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the 1080 D matrix. Otherwise no preconditioner is used. 1081 1082 Collective on PC 1083 1084 Input Parameters: 1085 + pc - the preconditioner context 1086 . ptype - which matrix to use for preconditioning the Schur complement 1087 - userpre - matrix to use for preconditioning, or PETSC_NULL 1088 1089 Notes: 1090 The default is to use the block on the diagonal of the preconditioning matrix. This is D, in the (1,1) position. 1091 There are currently no preconditioners that work directly with the Schur complement so setting 1092 PC_FIELDSPLIT_SCHUR_PRE_SELF is observationally equivalent to -fieldsplit_1_pc_type none. 1093 1094 Options Database: 1095 . -pc_fieldsplit_schur_precondition <self,user,diag> default is diag 1096 1097 Level: intermediate 1098 1099 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType 1100 1101 @*/ 1102 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1103 { 1104 PetscErrorCode ierr; 1105 1106 PetscFunctionBegin; 1107 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1108 ierr = PetscTryMethod(pc,"PCFieldSplitSchurPrecondition_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr); 1109 PetscFunctionReturn(0); 1110 } 1111 1112 EXTERN_C_BEGIN 1113 #undef __FUNCT__ 1114 #define __FUNCT__ "PCFieldSplitSchurPrecondition_FieldSplit" 1115 PetscErrorCode PCFieldSplitSchurPrecondition_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1116 { 1117 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1118 PetscErrorCode ierr; 1119 1120 PetscFunctionBegin; 1121 jac->schurpre = ptype; 1122 if (pre) { 1123 if (jac->schur_user) {ierr = MatDestroy(jac->schur_user);CHKERRQ(ierr);} 1124 jac->schur_user = pre; 1125 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 1126 } 1127 PetscFunctionReturn(0); 1128 } 1129 EXTERN_C_END 1130 1131 #undef __FUNCT__ 1132 #define __FUNCT__ "PCFieldSplitGetSchurBlocks" 1133 /*@C 1134 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 1135 1136 Collective on KSP 1137 1138 Input Parameter: 1139 . pc - the preconditioner context 1140 1141 Output Parameters: 1142 + A - the (0,0) block 1143 . B - the (0,1) block 1144 . C - the (1,0) block 1145 - D - the (1,1) block 1146 1147 Level: advanced 1148 1149 .seealso: PCFIELDSPLIT 1150 @*/ 1151 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A,Mat *B,Mat *C, Mat *D) 1152 { 1153 PC_FieldSplit *jac = (PC_FieldSplit *) pc->data; 1154 1155 PetscFunctionBegin; 1156 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1157 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 1158 if (A) *A = jac->pmat[0]; 1159 if (B) *B = jac->B; 1160 if (C) *C = jac->C; 1161 if (D) *D = jac->pmat[1]; 1162 PetscFunctionReturn(0); 1163 } 1164 1165 EXTERN_C_BEGIN 1166 #undef __FUNCT__ 1167 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 1168 PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 1169 { 1170 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1171 PetscErrorCode ierr; 1172 1173 PetscFunctionBegin; 1174 jac->type = type; 1175 if (type == PC_COMPOSITE_SCHUR) { 1176 pc->ops->apply = PCApply_FieldSplit_Schur; 1177 pc->ops->view = PCView_FieldSplit_Schur; 1178 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit_Schur",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 1179 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","PCFieldSplitSchurPrecondition_FieldSplit",PCFieldSplitSchurPrecondition_FieldSplit);CHKERRQ(ierr); 1180 1181 } else { 1182 pc->ops->apply = PCApply_FieldSplit; 1183 pc->ops->view = PCView_FieldSplit; 1184 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1185 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","",0);CHKERRQ(ierr); 1186 } 1187 PetscFunctionReturn(0); 1188 } 1189 EXTERN_C_END 1190 1191 EXTERN_C_BEGIN 1192 #undef __FUNCT__ 1193 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 1194 PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 1195 { 1196 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1197 1198 PetscFunctionBegin; 1199 if (bs < 1) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 1200 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); 1201 jac->bs = bs; 1202 PetscFunctionReturn(0); 1203 } 1204 EXTERN_C_END 1205 1206 #undef __FUNCT__ 1207 #define __FUNCT__ "PCFieldSplitSetType" 1208 /*@ 1209 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 1210 1211 Collective on PC 1212 1213 Input Parameter: 1214 . pc - the preconditioner context 1215 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1216 1217 Options Database Key: 1218 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 1219 1220 Level: Developer 1221 1222 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1223 1224 .seealso: PCCompositeSetType() 1225 1226 @*/ 1227 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 1228 { 1229 PetscErrorCode ierr; 1230 1231 PetscFunctionBegin; 1232 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1233 ierr = PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));CHKERRQ(ierr); 1234 PetscFunctionReturn(0); 1235 } 1236 1237 /* -------------------------------------------------------------------------------------*/ 1238 /*MC 1239 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 1240 fields or groups of fields 1241 1242 1243 To set options on the solvers for each block append -fieldsplit_ to all the PC 1244 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 1245 1246 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 1247 and set the options directly on the resulting KSP object 1248 1249 Level: intermediate 1250 1251 Options Database Keys: 1252 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 1253 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 1254 been supplied explicitly by -pc_fieldsplit_%d_fields 1255 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 1256 . -pc_fieldsplit_type <additive,multiplicative,schur,symmetric_multiplicative> 1257 . -pc_fieldsplit_schur_precondition <true,false> default is true 1258 . -pc_fieldsplit_stokes - automatically finds rows with zero diagonal and uses Schur complement with no preconditioner as the solver 1259 1260 - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ 1261 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 1262 1263 1264 Notes: use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 1265 to define a field by an arbitrary collection of entries. 1266 1267 If no fields are set the default is used. The fields are defined by entries strided by bs, 1268 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 1269 if this is not called the block size defaults to the blocksize of the second matrix passed 1270 to KSPSetOperators()/PCSetOperators(). 1271 1272 Currently for the multiplicative version, the updated residual needed for the next field 1273 solve is computed via a matrix vector product over the entire array. An optimization would be 1274 to update the residual only for the part of the right hand side associated with the next field 1275 solve. (This would involve more MatGetSubMatrix() calls or some other mechanism to compute the 1276 part of the matrix needed to just update part of the residual). 1277 1278 For the Schur complement preconditioner if J = ( A B ) 1279 ( C D ) 1280 the preconditioner is 1281 (I -B inv(A)) ( inv(A) 0 ) (I 0 ) 1282 (0 I ) ( 0 inv(S) ) (-C inv(A) I ) 1283 where the action of inv(A) is applied using the KSP solver with prefix -fieldsplit_0_. The action of 1284 inv(S) is computed using the KSP solver with prefix -schur_. For PCFieldSplitGetKSP() when field number is 1285 0 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 1286 D is used to construct a preconditioner for S, use PCFieldSplitSchurPrecondition() to turn on or off this 1287 option. 1288 1289 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 1290 is used automatically for a second block. 1291 1292 Concepts: physics based preconditioners, block preconditioners 1293 1294 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners 1295 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSchurPrecondition() 1296 M*/ 1297 1298 EXTERN_C_BEGIN 1299 #undef __FUNCT__ 1300 #define __FUNCT__ "PCCreate_FieldSplit" 1301 PetscErrorCode PCCreate_FieldSplit(PC pc) 1302 { 1303 PetscErrorCode ierr; 1304 PC_FieldSplit *jac; 1305 1306 PetscFunctionBegin; 1307 ierr = PetscNewLog(pc,PC_FieldSplit,&jac);CHKERRQ(ierr); 1308 jac->bs = -1; 1309 jac->nsplits = 0; 1310 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 1311 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 1312 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACTORIZATION_FULL; 1313 1314 pc->data = (void*)jac; 1315 1316 pc->ops->apply = PCApply_FieldSplit; 1317 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 1318 pc->ops->setup = PCSetUp_FieldSplit; 1319 pc->ops->destroy = PCDestroy_FieldSplit; 1320 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 1321 pc->ops->view = PCView_FieldSplit; 1322 pc->ops->applyrichardson = 0; 1323 1324 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit", 1325 PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1326 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","PCFieldSplitSetFields_FieldSplit", 1327 PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 1328 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","PCFieldSplitSetIS_FieldSplit", 1329 PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 1330 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","PCFieldSplitSetType_FieldSplit", 1331 PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 1332 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","PCFieldSplitSetBlockSize_FieldSplit", 1333 PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 1334 PetscFunctionReturn(0); 1335 } 1336 EXTERN_C_END 1337 1338 1339 /*MC 1340 Block_Preconditioners - PETSc provides support for a variety of "block preconditioners", this provides an 1341 overview of these methods. 1342 1343 Consider the solution to ( A B ) (x_1) = (b_1) 1344 ( C D ) (x_2) (b_2) 1345 1346 Important special cases, the Stokes equation: C = B' and D = 0 (A B) (x_1) = (b_1) 1347 B' 0) (x_2) (b_2) 1348 1349 One of the goals of the PCFieldSplit preconditioner in PETSc is to provide a variety of preconditioners 1350 for this block system. 1351 1352 Consider an additional matrix (Ap Bp) 1353 (Cp Dp) where some or all of the entries may be the same as 1354 in the original matrix (for example Ap == A). 1355 1356 In the following, A^ denotes the approximate application of the inverse of A, possibly using Ap in the 1357 approximation. In PETSc this simply means one has called KSPSetOperators(ksp,A,Ap,...) or KSPSetOperators(ksp,Ap,Ap,...) 1358 1359 Block Jacobi: x_1 = A^ b_1 1360 x_2 = D^ b_2 1361 1362 Lower block Gauss-Seidel: x_1 = A^ b_1 1363 x_2 = D^ (b_2 - C x_1) variant x_2 = D^ (b_2 - Cp x_1) 1364 1365 Symmetric Gauss-Seidel: x_1 = x_1 + A^(b_1 - A x_1 - B x_2) variant x_1 = x_1 + A^(b_1 - Ap x_1 - Bp x_2) 1366 Interestingly this form is not actually a symmetric matrix, the symmetric version is 1367 x_1 = A^(b_1 - B x_2) variant x_1 = A^(b_1 - Bp x_2) 1368 1369 Schur complement preconditioner 1370 the preconditioner is 1371 (I -B inv(A)) ( inv(A) 0 ) (I 0 ) 1372 (0 I ) ( 0 inv(S) ) (-C inv(A) I ) 1373 where the action of inv(A) is applied using the KSP solver with prefix -fieldsplit_0_. The action of 1374 inv(S) is computed using the KSP solver with prefix -schur_. For PCFieldSplitGetKSP() when field number is 1375 0 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 1376 D is used to construct a preconditioner for S, use PCFieldSplitSchurPrecondition() to turn on or off this 1377 option. 1378 1379 Level: intermediate 1380 1381 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCFIELDSPLIT 1382 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSchurPrecondition() 1383 M*/ 1384