1 2 3 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/ 4 #include <petsc/private/kspimpl.h> /* This is needed to provide the appropriate PETSC_EXTERN for KSP_Solve_FS ....*/ 5 #include <petscdm.h> 6 7 const char *const PCFieldSplitSchurPreTypes[] = {"SELF","SELFP","A11","USER","FULL","PCFieldSplitSchurPreType","PC_FIELDSPLIT_SCHUR_PRE_",0}; 8 const char *const PCFieldSplitSchurFactTypes[] = {"DIAG","LOWER","UPPER","FULL","PCFieldSplitSchurFactType","PC_FIELDSPLIT_SCHUR_FACT_",0}; 9 10 PetscLogEvent KSP_Solve_FS_0,KSP_Solve_FS_1,KSP_Solve_FS_S,KSP_Solve_FS_U,KSP_Solve_FS_L,KSP_Solve_FS_2,KSP_Solve_FS_3,KSP_Solve_FS_4; 11 12 typedef struct _PC_FieldSplitLink *PC_FieldSplitLink; 13 struct _PC_FieldSplitLink { 14 KSP ksp; 15 Vec x,y,z; 16 char *splitname; 17 PetscInt nfields; 18 PetscInt *fields,*fields_col; 19 VecScatter sctx; 20 IS is,is_col; 21 PC_FieldSplitLink next,previous; 22 PetscLogEvent event; 23 }; 24 25 typedef struct { 26 PCCompositeType type; 27 PetscBool defaultsplit; /* Flag for a system with a set of 'k' scalar fields with the same layout (and bs = k) */ 28 PetscBool splitdefined; /* Flag is set after the splits have been defined, to prevent more splits from being added */ 29 PetscInt bs; /* Block size for IS and Mat structures */ 30 PetscInt nsplits; /* Number of field divisions defined */ 31 Vec *x,*y,w1,w2; 32 Mat *mat; /* The diagonal block for each split */ 33 Mat *pmat; /* The preconditioning diagonal block for each split */ 34 Mat *Afield; /* The rows of the matrix associated with each split */ 35 PetscBool issetup; 36 37 /* Only used when Schur complement preconditioning is used */ 38 Mat B; /* The (0,1) block */ 39 Mat C; /* The (1,0) block */ 40 Mat schur; /* The Schur complement S = A11 - A10 A00^{-1} A01, the KSP here, kspinner, is H_1 in [El08] */ 41 Mat schurp; /* Assembled approximation to S built by MatSchurComplement to be used as a preconditioning matrix when solving with S */ 42 Mat schur_user; /* User-provided preconditioning matrix for the Schur complement */ 43 PCFieldSplitSchurPreType schurpre; /* Determines which preconditioning matrix is used for the Schur complement */ 44 PCFieldSplitSchurFactType schurfactorization; 45 KSP kspschur; /* The solver for S */ 46 KSP kspupper; /* The solver for A in the upper diagonal part of the factorization (H_2 in [El08]) */ 47 PetscScalar schurscale; /* Scaling factor for the Schur complement solution with DIAG factorization */ 48 49 PC_FieldSplitLink head; 50 PetscBool isrestrict; /* indicates PCFieldSplitRestrictIS() has been last called on this object, hack */ 51 PetscBool suboptionsset; /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */ 52 PetscBool dm_splits; /* Whether to use DMCreateFieldDecomposition() whenever possible */ 53 PetscBool diag_use_amat; /* Whether to extract diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */ 54 PetscBool offdiag_use_amat; /* Whether to extract off-diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */ 55 PetscBool detect; /* Whether to form 2-way split by finding zero diagonal entries */ 56 } PC_FieldSplit; 57 58 /* 59 Notes: 60 there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of 61 inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the 62 PC you could change this. 63 */ 64 65 /* This helper is so that setting a user-provided preconditioning matrix is orthogonal to choosing to use it. This way the 66 * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */ 67 static Mat FieldSplitSchurPre(PC_FieldSplit *jac) 68 { 69 switch (jac->schurpre) { 70 case PC_FIELDSPLIT_SCHUR_PRE_SELF: return jac->schur; 71 case PC_FIELDSPLIT_SCHUR_PRE_SELFP: return jac->schurp; 72 case PC_FIELDSPLIT_SCHUR_PRE_A11: return jac->pmat[1]; 73 case PC_FIELDSPLIT_SCHUR_PRE_FULL: /* We calculate this and store it in schur_user */ 74 case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */ 75 default: 76 return jac->schur_user ? jac->schur_user : jac->pmat[1]; 77 } 78 } 79 80 81 #include <petscdraw.h> 82 static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer) 83 { 84 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 85 PetscErrorCode ierr; 86 PetscBool iascii,isdraw; 87 PetscInt i,j; 88 PC_FieldSplitLink ilink = jac->head; 89 90 PetscFunctionBegin; 91 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 92 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 93 if (iascii) { 94 if (jac->bs > 0) { 95 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);CHKERRQ(ierr); 96 } else { 97 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D\n",PCCompositeTypes[jac->type],jac->nsplits);CHKERRQ(ierr); 98 } 99 if (pc->useAmat) { 100 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");CHKERRQ(ierr); 101 } 102 if (jac->diag_use_amat) { 103 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for diagonal blocks\n");CHKERRQ(ierr); 104 } 105 if (jac->offdiag_use_amat) { 106 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for off-diagonal blocks\n");CHKERRQ(ierr); 107 } 108 ierr = PetscViewerASCIIPrintf(viewer," Solver info for each split is in the following KSP objects:\n");CHKERRQ(ierr); 109 for (i=0; i<jac->nsplits; i++) { 110 if (ilink->fields) { 111 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr); 112 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 113 for (j=0; j<ilink->nfields; j++) { 114 if (j > 0) { 115 ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr); 116 } 117 ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr); 118 } 119 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 120 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 121 } else { 122 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);CHKERRQ(ierr); 123 } 124 ierr = KSPView(ilink->ksp,viewer);CHKERRQ(ierr); 125 ilink = ilink->next; 126 } 127 } 128 129 if (isdraw) { 130 PetscDraw draw; 131 PetscReal x,y,w,wd; 132 133 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 134 ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr); 135 w = 2*PetscMin(1.0 - x,x); 136 wd = w/(jac->nsplits + 1); 137 x = x - wd*(jac->nsplits-1)/2.0; 138 for (i=0; i<jac->nsplits; i++) { 139 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr); 140 ierr = KSPView(ilink->ksp,viewer);CHKERRQ(ierr); 141 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 142 x += wd; 143 ilink = ilink->next; 144 } 145 } 146 PetscFunctionReturn(0); 147 } 148 149 static PetscErrorCode PCView_FieldSplit_Schur(PC pc,PetscViewer viewer) 150 { 151 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 152 PetscErrorCode ierr; 153 PetscBool iascii,isdraw; 154 PetscInt i,j; 155 PC_FieldSplitLink ilink = jac->head; 156 MatSchurComplementAinvType atype; 157 158 PetscFunctionBegin; 159 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 160 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 161 if (iascii) { 162 if (jac->bs > 0) { 163 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, blocksize = %D, factorization %s\n",jac->bs,PCFieldSplitSchurFactTypes[jac->schurfactorization]);CHKERRQ(ierr); 164 } else { 165 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, factorization %s\n",PCFieldSplitSchurFactTypes[jac->schurfactorization]);CHKERRQ(ierr); 166 } 167 if (pc->useAmat) { 168 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");CHKERRQ(ierr); 169 } 170 switch (jac->schurpre) { 171 case PC_FIELDSPLIT_SCHUR_PRE_SELF: 172 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from S itself\n");CHKERRQ(ierr); 173 break; 174 case PC_FIELDSPLIT_SCHUR_PRE_SELFP: 175 ierr = MatSchurComplementGetAinvType(jac->schur,&atype);CHKERRQ(ierr); 176 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from Sp, an assembled approximation to S, which uses A00's %sdiagonal's inverse\n",atype == MAT_SCHUR_COMPLEMENT_AINV_DIAG ? "" : (atype == MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG ? "block " : "lumped "));CHKERRQ(ierr);break; 177 case PC_FIELDSPLIT_SCHUR_PRE_A11: 178 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");CHKERRQ(ierr); 179 break; 180 case PC_FIELDSPLIT_SCHUR_PRE_FULL: 181 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from the exact Schur complement\n");CHKERRQ(ierr); 182 break; 183 case PC_FIELDSPLIT_SCHUR_PRE_USER: 184 if (jac->schur_user) { 185 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from user provided matrix\n");CHKERRQ(ierr); 186 } else { 187 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");CHKERRQ(ierr); 188 } 189 break; 190 default: 191 SETERRQ1(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Invalid Schur preconditioning type: %d", jac->schurpre); 192 } 193 ierr = PetscViewerASCIIPrintf(viewer," Split info:\n");CHKERRQ(ierr); 194 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 195 for (i=0; i<jac->nsplits; i++) { 196 if (ilink->fields) { 197 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr); 198 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 199 for (j=0; j<ilink->nfields; j++) { 200 if (j > 0) { 201 ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr); 202 } 203 ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr); 204 } 205 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 206 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 207 } else { 208 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);CHKERRQ(ierr); 209 } 210 ilink = ilink->next; 211 } 212 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for A00 block\n");CHKERRQ(ierr); 213 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 214 if (jac->head) { 215 ierr = KSPView(jac->head->ksp,viewer);CHKERRQ(ierr); 216 } else {ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr);} 217 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 218 if (jac->head && jac->kspupper != jac->head->ksp) { 219 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for upper A00 in upper triangular factor \n");CHKERRQ(ierr); 220 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 221 if (jac->kspupper) {ierr = KSPView(jac->kspupper,viewer);CHKERRQ(ierr);} 222 else {ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr);} 223 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 224 } 225 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for S = A11 - A10 inv(A00) A01 \n");CHKERRQ(ierr); 226 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 227 if (jac->kspschur) { 228 ierr = KSPView(jac->kspschur,viewer);CHKERRQ(ierr); 229 } else { 230 ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr); 231 } 232 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 233 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 234 } else if (isdraw && jac->head) { 235 PetscDraw draw; 236 PetscReal x,y,w,wd,h; 237 PetscInt cnt = 2; 238 char str[32]; 239 240 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 241 ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr); 242 if (jac->kspupper != jac->head->ksp) cnt++; 243 w = 2*PetscMin(1.0 - x,x); 244 wd = w/(cnt + 1); 245 246 ierr = PetscSNPrintf(str,32,"Schur fact. %s",PCFieldSplitSchurFactTypes[jac->schurfactorization]);CHKERRQ(ierr); 247 ierr = PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);CHKERRQ(ierr); 248 y -= h; 249 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_USER && !jac->schur_user) { 250 ierr = PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[PC_FIELDSPLIT_SCHUR_PRE_A11]);CHKERRQ(ierr); 251 } else { 252 ierr = PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[jac->schurpre]);CHKERRQ(ierr); 253 } 254 ierr = PetscDrawStringBoxed(draw,x+wd*(cnt-1)/2.0,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);CHKERRQ(ierr); 255 y -= h; 256 x = x - wd*(cnt-1)/2.0; 257 258 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr); 259 ierr = KSPView(jac->head->ksp,viewer);CHKERRQ(ierr); 260 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 261 if (jac->kspupper != jac->head->ksp) { 262 x += wd; 263 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr); 264 ierr = KSPView(jac->kspupper,viewer);CHKERRQ(ierr); 265 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 266 } 267 x += wd; 268 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr); 269 ierr = KSPView(jac->kspschur,viewer);CHKERRQ(ierr); 270 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 271 } 272 PetscFunctionReturn(0); 273 } 274 275 /* Precondition: jac->bs is set to a meaningful value */ 276 static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc) 277 { 278 PetscErrorCode ierr; 279 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 280 PetscInt i,nfields,*ifields,nfields_col,*ifields_col; 281 PetscBool flg,flg_col; 282 char optionname[128],splitname[8],optionname_col[128]; 283 284 PetscFunctionBegin; 285 ierr = PetscMalloc1(jac->bs,&ifields);CHKERRQ(ierr); 286 ierr = PetscMalloc1(jac->bs,&ifields_col);CHKERRQ(ierr); 287 for (i=0,flg=PETSC_TRUE;; i++) { 288 ierr = PetscSNPrintf(splitname,sizeof(splitname),"%D",i);CHKERRQ(ierr); 289 ierr = PetscSNPrintf(optionname,sizeof(optionname),"-pc_fieldsplit_%D_fields",i);CHKERRQ(ierr); 290 ierr = PetscSNPrintf(optionname_col,sizeof(optionname_col),"-pc_fieldsplit_%D_fields_col",i);CHKERRQ(ierr); 291 nfields = jac->bs; 292 nfields_col = jac->bs; 293 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname,ifields,&nfields,&flg);CHKERRQ(ierr); 294 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname_col,ifields_col,&nfields_col,&flg_col);CHKERRQ(ierr); 295 if (!flg) break; 296 else if (flg && !flg_col) { 297 if (!nfields) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields"); 298 ierr = PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields);CHKERRQ(ierr); 299 } else { 300 if (!nfields || !nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields"); 301 if (nfields != nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Number of row and column fields must match"); 302 ierr = PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields_col);CHKERRQ(ierr); 303 } 304 } 305 if (i > 0) { 306 /* Makes command-line setting of splits take precedence over setting them in code. 307 Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would 308 create new splits, which would probably not be what the user wanted. */ 309 jac->splitdefined = PETSC_TRUE; 310 } 311 ierr = PetscFree(ifields);CHKERRQ(ierr); 312 ierr = PetscFree(ifields_col);CHKERRQ(ierr); 313 PetscFunctionReturn(0); 314 } 315 316 static PetscErrorCode PCFieldSplitSetDefaults(PC pc) 317 { 318 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 319 PetscErrorCode ierr; 320 PC_FieldSplitLink ilink = jac->head; 321 PetscBool fieldsplit_default = PETSC_FALSE,coupling = PETSC_FALSE; 322 PetscInt i; 323 324 PetscFunctionBegin; 325 /* 326 Kinda messy, but at least this now uses DMCreateFieldDecomposition(). 327 Should probably be rewritten. 328 */ 329 if (!ilink) { 330 ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_coupling",&coupling,NULL);CHKERRQ(ierr); 331 if (pc->dm && jac->dm_splits && !jac->detect && !coupling) { 332 PetscInt numFields, f, i, j; 333 char **fieldNames; 334 IS *fields; 335 DM *dms; 336 DM subdm[128]; 337 PetscBool flg; 338 339 ierr = DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms);CHKERRQ(ierr); 340 /* Allow the user to prescribe the splits */ 341 for (i = 0, flg = PETSC_TRUE;; i++) { 342 PetscInt ifields[128]; 343 IS compField; 344 char optionname[128], splitname[8]; 345 PetscInt nfields = numFields; 346 347 ierr = PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%D_fields", i);CHKERRQ(ierr); 348 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg);CHKERRQ(ierr); 349 if (!flg) break; 350 if (numFields > 128) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot currently support %d > 128 fields", numFields); 351 ierr = DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i]);CHKERRQ(ierr); 352 if (nfields == 1) { 353 ierr = PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField);CHKERRQ(ierr); 354 } else { 355 ierr = PetscSNPrintf(splitname, sizeof(splitname), "%D", i);CHKERRQ(ierr); 356 ierr = PCFieldSplitSetIS(pc, splitname, compField);CHKERRQ(ierr); 357 } 358 ierr = ISDestroy(&compField);CHKERRQ(ierr); 359 for (j = 0; j < nfields; ++j) { 360 f = ifields[j]; 361 ierr = PetscFree(fieldNames[f]);CHKERRQ(ierr); 362 ierr = ISDestroy(&fields[f]);CHKERRQ(ierr); 363 } 364 } 365 if (i == 0) { 366 for (f = 0; f < numFields; ++f) { 367 ierr = PCFieldSplitSetIS(pc, fieldNames[f], fields[f]);CHKERRQ(ierr); 368 ierr = PetscFree(fieldNames[f]);CHKERRQ(ierr); 369 ierr = ISDestroy(&fields[f]);CHKERRQ(ierr); 370 } 371 } else { 372 for (j=0; j<numFields; j++) { 373 ierr = DMDestroy(dms+j);CHKERRQ(ierr); 374 } 375 ierr = PetscFree(dms);CHKERRQ(ierr); 376 ierr = PetscMalloc1(i, &dms);CHKERRQ(ierr); 377 for (j = 0; j < i; ++j) dms[j] = subdm[j]; 378 } 379 ierr = PetscFree(fieldNames);CHKERRQ(ierr); 380 ierr = PetscFree(fields);CHKERRQ(ierr); 381 if (dms) { 382 ierr = PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n");CHKERRQ(ierr); 383 for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) { 384 const char *prefix; 385 ierr = PetscObjectGetOptionsPrefix((PetscObject)(ilink->ksp),&prefix);CHKERRQ(ierr); 386 ierr = PetscObjectSetOptionsPrefix((PetscObject)(dms[i]), prefix);CHKERRQ(ierr); 387 ierr = KSPSetDM(ilink->ksp, dms[i]);CHKERRQ(ierr); 388 ierr = KSPSetDMActive(ilink->ksp, PETSC_FALSE);CHKERRQ(ierr); 389 ierr = PetscObjectIncrementTabLevel((PetscObject)dms[i],(PetscObject)ilink->ksp,0);CHKERRQ(ierr); 390 ierr = DMDestroy(&dms[i]);CHKERRQ(ierr); 391 } 392 ierr = PetscFree(dms);CHKERRQ(ierr); 393 } 394 } else { 395 if (jac->bs <= 0) { 396 if (pc->pmat) { 397 ierr = MatGetBlockSize(pc->pmat,&jac->bs);CHKERRQ(ierr); 398 } else jac->bs = 1; 399 } 400 401 if (jac->detect) { 402 IS zerodiags,rest; 403 PetscInt nmin,nmax; 404 405 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 406 ierr = MatFindZeroDiagonals(pc->mat,&zerodiags);CHKERRQ(ierr); 407 ierr = ISComplement(zerodiags,nmin,nmax,&rest);CHKERRQ(ierr); 408 ierr = PCFieldSplitSetIS(pc,"0",rest);CHKERRQ(ierr); 409 ierr = PCFieldSplitSetIS(pc,"1",zerodiags);CHKERRQ(ierr); 410 ierr = ISDestroy(&zerodiags);CHKERRQ(ierr); 411 ierr = ISDestroy(&rest);CHKERRQ(ierr); 412 } else if (coupling) { 413 IS coupling,rest; 414 PetscInt nmin,nmax; 415 416 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 417 ierr = MatFindOffBlockDiagonalEntries(pc->mat,&coupling);CHKERRQ(ierr); 418 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc->mat),nmax-nmin,nmin,1,&rest);CHKERRQ(ierr); 419 ierr = ISSetIdentity(rest);CHKERRQ(ierr); 420 ierr = PCFieldSplitSetIS(pc,"0",rest);CHKERRQ(ierr); 421 ierr = PCFieldSplitSetIS(pc,"1",coupling);CHKERRQ(ierr); 422 ierr = ISDestroy(&coupling);CHKERRQ(ierr); 423 ierr = ISDestroy(&rest);CHKERRQ(ierr); 424 } else { 425 ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&fieldsplit_default,NULL);CHKERRQ(ierr); 426 if (!fieldsplit_default) { 427 /* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit() 428 then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */ 429 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 430 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 431 } 432 if ((fieldsplit_default || !jac->splitdefined) && !jac->isrestrict) { 433 ierr = PetscInfo(pc,"Using default splitting of fields\n");CHKERRQ(ierr); 434 for (i=0; i<jac->bs; i++) { 435 char splitname[8]; 436 ierr = PetscSNPrintf(splitname,sizeof(splitname),"%D",i);CHKERRQ(ierr); 437 ierr = PCFieldSplitSetFields(pc,splitname,1,&i,&i);CHKERRQ(ierr); 438 } 439 jac->defaultsplit = PETSC_TRUE; 440 } 441 } 442 } 443 } else if (jac->nsplits == 1) { 444 if (ilink->is) { 445 IS is2; 446 PetscInt nmin,nmax; 447 448 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 449 ierr = ISComplement(ilink->is,nmin,nmax,&is2);CHKERRQ(ierr); 450 ierr = PCFieldSplitSetIS(pc,"1",is2);CHKERRQ(ierr); 451 ierr = ISDestroy(&is2);CHKERRQ(ierr); 452 } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Must provide at least two sets of fields to PCFieldSplit()"); 453 } 454 455 if (jac->nsplits < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unhandled case, must have at least two fields, not %d", jac->nsplits); 456 PetscFunctionReturn(0); 457 } 458 459 PETSC_EXTERN PetscErrorCode PetscOptionsFindPairPrefix_Private(PetscOptions,const char pre[], const char name[], char *value[], PetscBool *flg); 460 461 static PetscErrorCode PCSetUp_FieldSplit(PC pc) 462 { 463 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 464 PetscErrorCode ierr; 465 PC_FieldSplitLink ilink; 466 PetscInt i,nsplit; 467 PetscBool sorted, sorted_col; 468 469 PetscFunctionBegin; 470 ierr = PCFieldSplitSetDefaults(pc);CHKERRQ(ierr); 471 nsplit = jac->nsplits; 472 ilink = jac->head; 473 474 /* get the matrices for each split */ 475 if (!jac->issetup) { 476 PetscInt rstart,rend,nslots,bs; 477 478 jac->issetup = PETSC_TRUE; 479 480 /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */ 481 if (jac->defaultsplit || !ilink->is) { 482 if (jac->bs <= 0) jac->bs = nsplit; 483 } 484 bs = jac->bs; 485 ierr = MatGetOwnershipRange(pc->pmat,&rstart,&rend);CHKERRQ(ierr); 486 nslots = (rend - rstart)/bs; 487 for (i=0; i<nsplit; i++) { 488 if (jac->defaultsplit) { 489 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+i,nsplit,&ilink->is);CHKERRQ(ierr); 490 ierr = ISDuplicate(ilink->is,&ilink->is_col);CHKERRQ(ierr); 491 } else if (!ilink->is) { 492 if (ilink->nfields > 1) { 493 PetscInt *ii,*jj,j,k,nfields = ilink->nfields,*fields = ilink->fields,*fields_col = ilink->fields_col; 494 ierr = PetscMalloc1(ilink->nfields*nslots,&ii);CHKERRQ(ierr); 495 ierr = PetscMalloc1(ilink->nfields*nslots,&jj);CHKERRQ(ierr); 496 for (j=0; j<nslots; j++) { 497 for (k=0; k<nfields; k++) { 498 ii[nfields*j + k] = rstart + bs*j + fields[k]; 499 jj[nfields*j + k] = rstart + bs*j + fields_col[k]; 500 } 501 } 502 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is);CHKERRQ(ierr); 503 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,jj,PETSC_OWN_POINTER,&ilink->is_col);CHKERRQ(ierr); 504 ierr = ISSetBlockSize(ilink->is, nfields);CHKERRQ(ierr); 505 ierr = ISSetBlockSize(ilink->is_col, nfields);CHKERRQ(ierr); 506 } else { 507 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields[0],bs,&ilink->is);CHKERRQ(ierr); 508 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields_col[0],bs,&ilink->is_col);CHKERRQ(ierr); 509 } 510 } 511 ierr = ISSorted(ilink->is,&sorted);CHKERRQ(ierr); 512 if (ilink->is_col) { ierr = ISSorted(ilink->is_col,&sorted_col);CHKERRQ(ierr); } 513 if (!sorted || !sorted_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Fields must be sorted when creating split"); 514 ilink = ilink->next; 515 } 516 } 517 518 ilink = jac->head; 519 if (!jac->pmat) { 520 Vec xtmp; 521 522 ierr = MatCreateVecs(pc->pmat,&xtmp,NULL);CHKERRQ(ierr); 523 ierr = PetscMalloc1(nsplit,&jac->pmat);CHKERRQ(ierr); 524 ierr = PetscMalloc2(nsplit,&jac->x,nsplit,&jac->y);CHKERRQ(ierr); 525 for (i=0; i<nsplit; i++) { 526 MatNullSpace sp; 527 528 /* Check for preconditioning matrix attached to IS */ 529 ierr = PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &jac->pmat[i]);CHKERRQ(ierr); 530 if (jac->pmat[i]) { 531 ierr = PetscObjectReference((PetscObject) jac->pmat[i]);CHKERRQ(ierr); 532 if (jac->type == PC_COMPOSITE_SCHUR) { 533 jac->schur_user = jac->pmat[i]; 534 535 ierr = PetscObjectReference((PetscObject) jac->schur_user);CHKERRQ(ierr); 536 } 537 } else { 538 const char *prefix; 539 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 540 ierr = KSPGetOptionsPrefix(ilink->ksp,&prefix);CHKERRQ(ierr); 541 ierr = MatSetOptionsPrefix(jac->pmat[i],prefix);CHKERRQ(ierr); 542 ierr = MatViewFromOptions(jac->pmat[i],NULL,"-mat_view");CHKERRQ(ierr); 543 } 544 /* create work vectors for each split */ 545 ierr = MatCreateVecs(jac->pmat[i],&jac->x[i],&jac->y[i]);CHKERRQ(ierr); 546 ilink->x = jac->x[i]; ilink->y = jac->y[i]; ilink->z = NULL; 547 /* compute scatter contexts needed by multiplicative versions and non-default splits */ 548 ierr = VecScatterCreate(xtmp,ilink->is,jac->x[i],NULL,&ilink->sctx);CHKERRQ(ierr); 549 ierr = PetscObjectQuery((PetscObject) ilink->is, "nearnullspace", (PetscObject*) &sp);CHKERRQ(ierr); 550 if (sp) { 551 ierr = MatSetNearNullSpace(jac->pmat[i], sp);CHKERRQ(ierr); 552 } 553 ilink = ilink->next; 554 } 555 ierr = VecDestroy(&xtmp);CHKERRQ(ierr); 556 } else { 557 MatReuse scall; 558 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 559 for (i=0; i<nsplit; i++) { 560 ierr = MatDestroy(&jac->pmat[i]);CHKERRQ(ierr); 561 } 562 scall = MAT_INITIAL_MATRIX; 563 } else scall = MAT_REUSE_MATRIX; 564 565 for (i=0; i<nsplit; i++) { 566 Mat pmat; 567 568 /* Check for preconditioning matrix attached to IS */ 569 ierr = PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &pmat);CHKERRQ(ierr); 570 if (!pmat) { 571 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is_col,scall,&jac->pmat[i]);CHKERRQ(ierr); 572 } 573 ilink = ilink->next; 574 } 575 } 576 if (jac->diag_use_amat) { 577 ilink = jac->head; 578 if (!jac->mat) { 579 ierr = PetscMalloc1(nsplit,&jac->mat);CHKERRQ(ierr); 580 for (i=0; i<nsplit; i++) { 581 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->mat[i]);CHKERRQ(ierr); 582 ilink = ilink->next; 583 } 584 } else { 585 MatReuse scall; 586 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 587 for (i=0; i<nsplit; i++) { 588 ierr = MatDestroy(&jac->mat[i]);CHKERRQ(ierr); 589 } 590 scall = MAT_INITIAL_MATRIX; 591 } else scall = MAT_REUSE_MATRIX; 592 593 for (i=0; i<nsplit; i++) { 594 if (jac->mat[i]) {ierr = MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,scall,&jac->mat[i]);CHKERRQ(ierr);} 595 ilink = ilink->next; 596 } 597 } 598 } else { 599 jac->mat = jac->pmat; 600 } 601 602 /* Check for null space attached to IS */ 603 ilink = jac->head; 604 for (i=0; i<nsplit; i++) { 605 MatNullSpace sp; 606 607 ierr = PetscObjectQuery((PetscObject) ilink->is, "nullspace", (PetscObject*) &sp);CHKERRQ(ierr); 608 if (sp) { 609 ierr = MatSetNullSpace(jac->mat[i], sp);CHKERRQ(ierr); 610 } 611 ilink = ilink->next; 612 } 613 614 if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR) { 615 /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */ 616 /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */ 617 ilink = jac->head; 618 if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) { 619 /* special case need where Afield[0] is not needed and only certain columns of Afield[1] are needed since update is only on those rows of the solution */ 620 if (!jac->Afield) { 621 ierr = PetscCalloc1(nsplit,&jac->Afield);CHKERRQ(ierr); 622 if (jac->offdiag_use_amat) { 623 ierr = MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);CHKERRQ(ierr); 624 } else { 625 ierr = MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);CHKERRQ(ierr); 626 } 627 } else { 628 MatReuse scall; 629 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 630 for (i=0; i<nsplit; i++) { 631 ierr = MatDestroy(&jac->Afield[1]);CHKERRQ(ierr); 632 } 633 scall = MAT_INITIAL_MATRIX; 634 } else scall = MAT_REUSE_MATRIX; 635 636 if (jac->offdiag_use_amat) { 637 ierr = MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,scall,&jac->Afield[1]);CHKERRQ(ierr); 638 } else { 639 ierr = MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,scall,&jac->Afield[1]);CHKERRQ(ierr); 640 } 641 } 642 } else { 643 if (!jac->Afield) { 644 ierr = PetscMalloc1(nsplit,&jac->Afield);CHKERRQ(ierr); 645 for (i=0; i<nsplit; i++) { 646 if (jac->offdiag_use_amat) { 647 ierr = MatCreateSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 648 } else { 649 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 650 } 651 ilink = ilink->next; 652 } 653 } else { 654 MatReuse scall; 655 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 656 for (i=0; i<nsplit; i++) { 657 ierr = MatDestroy(&jac->Afield[i]);CHKERRQ(ierr); 658 } 659 scall = MAT_INITIAL_MATRIX; 660 } else scall = MAT_REUSE_MATRIX; 661 662 for (i=0; i<nsplit; i++) { 663 if (jac->offdiag_use_amat) { 664 ierr = MatCreateSubMatrix(pc->mat,ilink->is,NULL,scall,&jac->Afield[i]);CHKERRQ(ierr); 665 } else { 666 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,NULL,scall,&jac->Afield[i]);CHKERRQ(ierr); 667 } 668 ilink = ilink->next; 669 } 670 } 671 } 672 } 673 674 if (jac->type == PC_COMPOSITE_SCHUR) { 675 IS ccis; 676 PetscBool isspd; 677 PetscInt rstart,rend; 678 char lscname[256]; 679 PetscObject LSC_L; 680 681 if (nsplit != 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields"); 682 683 /* If pc->mat is SPD, don't scale by -1 the Schur complement */ 684 if (jac->schurscale == (PetscScalar)-1.0) { 685 ierr = MatGetOption(pc->pmat,MAT_SPD,&isspd);CHKERRQ(ierr); 686 jac->schurscale = (isspd == PETSC_TRUE) ? 1.0 : -1.0; 687 } 688 689 /* When extracting off-diagonal submatrices, we take complements from this range */ 690 ierr = MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);CHKERRQ(ierr); 691 692 /* need to handle case when one is resetting up the preconditioner */ 693 if (jac->schur) { 694 KSP kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper; 695 696 ierr = MatSchurComplementGetKSP(jac->schur, &kspInner);CHKERRQ(ierr); 697 ilink = jac->head; 698 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 699 if (jac->offdiag_use_amat) { 700 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);CHKERRQ(ierr); 701 } else { 702 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);CHKERRQ(ierr); 703 } 704 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 705 ilink = ilink->next; 706 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 707 if (jac->offdiag_use_amat) { 708 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);CHKERRQ(ierr); 709 } else { 710 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);CHKERRQ(ierr); 711 } 712 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 713 ierr = MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);CHKERRQ(ierr); 714 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) { 715 ierr = MatDestroy(&jac->schurp);CHKERRQ(ierr); 716 ierr = MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);CHKERRQ(ierr); 717 } 718 if (kspA != kspInner) { 719 ierr = KSPSetOperators(kspA,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 720 } 721 if (kspUpper != kspA) { 722 ierr = KSPSetOperators(kspUpper,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 723 } 724 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));CHKERRQ(ierr); 725 } else { 726 const char *Dprefix; 727 char schurprefix[256], schurmatprefix[256]; 728 char schurtestoption[256]; 729 MatNullSpace sp; 730 PetscBool flg; 731 732 /* extract the A01 and A10 matrices */ 733 ilink = jac->head; 734 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 735 if (jac->offdiag_use_amat) { 736 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr); 737 } else { 738 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr); 739 } 740 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 741 ilink = ilink->next; 742 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 743 if (jac->offdiag_use_amat) { 744 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 745 } else { 746 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 747 } 748 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 749 750 /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */ 751 ierr = MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);CHKERRQ(ierr); 752 ierr = MatSetType(jac->schur,MATSCHURCOMPLEMENT);CHKERRQ(ierr); 753 ierr = MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);CHKERRQ(ierr); 754 ierr = PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 755 /* Note that the inner KSP is NOT going to inherit this prefix, and if it did, it would be reset just below. Is that what we want? */ 756 ierr = MatSetOptionsPrefix(jac->schur,schurmatprefix);CHKERRQ(ierr); 757 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr); 758 ierr = MatGetNullSpace(jac->mat[1], &sp);CHKERRQ(ierr); 759 if (sp) { 760 ierr = MatSetNullSpace(jac->schur, sp);CHKERRQ(ierr); 761 } 762 763 ierr = PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);CHKERRQ(ierr); 764 ierr = PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);CHKERRQ(ierr); 765 if (flg) { 766 DM dmInner; 767 KSP kspInner; 768 769 ierr = MatSchurComplementGetKSP(jac->schur, &kspInner);CHKERRQ(ierr); 770 ierr = KSPReset(kspInner);CHKERRQ(ierr); 771 ierr = KSPSetOperators(kspInner,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 772 ierr = PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 773 /* Indent this deeper to emphasize the "inner" nature of this solver. */ 774 ierr = PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject) pc, 2);CHKERRQ(ierr); 775 ierr = PetscObjectIncrementTabLevel((PetscObject)kspInner->pc, (PetscObject) pc, 2);CHKERRQ(ierr); 776 ierr = KSPSetOptionsPrefix(kspInner, schurprefix);CHKERRQ(ierr); 777 778 /* Set DM for new solver */ 779 ierr = KSPGetDM(jac->head->ksp, &dmInner);CHKERRQ(ierr); 780 ierr = KSPSetDM(kspInner, dmInner);CHKERRQ(ierr); 781 ierr = KSPSetDMActive(kspInner, PETSC_FALSE);CHKERRQ(ierr); 782 } else { 783 /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or 784 * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact, 785 * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for 786 * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make 787 * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used 788 * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */ 789 ierr = KSPSetType(jac->head->ksp,KSPGMRES);CHKERRQ(ierr); 790 ierr = MatSchurComplementSetKSP(jac->schur,jac->head->ksp);CHKERRQ(ierr); 791 } 792 ierr = KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 793 ierr = KSPSetFromOptions(jac->head->ksp);CHKERRQ(ierr); 794 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr); 795 796 ierr = PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);CHKERRQ(ierr); 797 ierr = PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);CHKERRQ(ierr); 798 if (flg) { 799 DM dmInner; 800 801 ierr = PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 802 ierr = KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);CHKERRQ(ierr); 803 ierr = KSPSetErrorIfNotConverged(jac->kspupper,pc->erroriffailure);CHKERRQ(ierr); 804 ierr = KSPSetOptionsPrefix(jac->kspupper, schurprefix);CHKERRQ(ierr); 805 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspupper, (PetscObject) pc, 1);CHKERRQ(ierr); 806 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspupper->pc, (PetscObject) pc, 1);CHKERRQ(ierr); 807 ierr = KSPGetDM(jac->head->ksp, &dmInner);CHKERRQ(ierr); 808 ierr = KSPSetDM(jac->kspupper, dmInner);CHKERRQ(ierr); 809 ierr = KSPSetDMActive(jac->kspupper, PETSC_FALSE);CHKERRQ(ierr); 810 ierr = KSPSetFromOptions(jac->kspupper);CHKERRQ(ierr); 811 ierr = KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 812 ierr = VecDuplicate(jac->head->x, &jac->head->z);CHKERRQ(ierr); 813 } else { 814 jac->kspupper = jac->head->ksp; 815 ierr = PetscObjectReference((PetscObject) jac->head->ksp);CHKERRQ(ierr); 816 } 817 818 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) { 819 ierr = MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);CHKERRQ(ierr); 820 } 821 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);CHKERRQ(ierr); 822 ierr = KSPSetErrorIfNotConverged(jac->kspschur,pc->erroriffailure);CHKERRQ(ierr); 823 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);CHKERRQ(ierr); 824 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);CHKERRQ(ierr); 825 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) { 826 PC pcschur; 827 ierr = KSPGetPC(jac->kspschur,&pcschur);CHKERRQ(ierr); 828 ierr = PCSetType(pcschur,PCNONE);CHKERRQ(ierr); 829 /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */ 830 } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) { 831 ierr = MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user);CHKERRQ(ierr); 832 } 833 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));CHKERRQ(ierr); 834 ierr = KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix);CHKERRQ(ierr); 835 ierr = KSPSetOptionsPrefix(jac->kspschur, Dprefix);CHKERRQ(ierr); 836 /* propagate DM */ 837 { 838 DM sdm; 839 ierr = KSPGetDM(jac->head->next->ksp, &sdm);CHKERRQ(ierr); 840 if (sdm) { 841 ierr = KSPSetDM(jac->kspschur, sdm);CHKERRQ(ierr); 842 ierr = KSPSetDMActive(jac->kspschur, PETSC_FALSE);CHKERRQ(ierr); 843 } 844 } 845 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 846 /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */ 847 ierr = KSPSetFromOptions(jac->kspschur);CHKERRQ(ierr); 848 } 849 850 /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */ 851 ierr = PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname);CHKERRQ(ierr); 852 ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 853 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 854 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);CHKERRQ(ierr);} 855 ierr = PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname);CHKERRQ(ierr); 856 ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 857 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 858 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);CHKERRQ(ierr);} 859 } else { 860 /* set up the individual splits' PCs */ 861 i = 0; 862 ilink = jac->head; 863 while (ilink) { 864 ierr = KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i]);CHKERRQ(ierr); 865 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 866 if (!jac->suboptionsset) {ierr = KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr);} 867 i++; 868 ilink = ilink->next; 869 } 870 } 871 872 jac->suboptionsset = PETSC_TRUE; 873 PetscFunctionReturn(0); 874 } 875 876 #define FieldSplitSplitSolveAdd(ilink,xx,yy) \ 877 (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 878 VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 879 PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\ 880 KSPSolve(ilink->ksp,ilink->x,ilink->y) || \ 881 PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\ 882 VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \ 883 VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE)) 884 885 static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y) 886 { 887 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 888 PetscErrorCode ierr; 889 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 890 KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper; 891 892 PetscFunctionBegin; 893 switch (jac->schurfactorization) { 894 case PC_FIELDSPLIT_SCHUR_FACT_DIAG: 895 /* [A00 0; 0 -S], positive definite, suitable for MINRES */ 896 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 897 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 898 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 899 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 900 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 901 ierr = PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 902 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 903 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 904 ierr = PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 905 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 906 ierr = PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 907 ierr = VecScale(ilinkD->y,jac->schurscale);CHKERRQ(ierr); 908 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 909 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 910 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 911 break; 912 case PC_FIELDSPLIT_SCHUR_FACT_LOWER: 913 /* [A00 0; A10 S], suitable for left preconditioning */ 914 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 915 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 916 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 917 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 918 ierr = PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 919 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 920 ierr = VecScale(ilinkD->x,-1.);CHKERRQ(ierr); 921 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 922 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 923 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 924 ierr = PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 925 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 926 ierr = PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 927 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 928 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 929 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 930 break; 931 case PC_FIELDSPLIT_SCHUR_FACT_UPPER: 932 /* [A00 A01; 0 S], suitable for right preconditioning */ 933 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 934 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 935 ierr = PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 936 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 937 ierr = PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 938 ierr = VecScale(ilinkA->x,-1.);CHKERRQ(ierr); 939 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 940 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 941 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 942 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 943 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 944 ierr = PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 945 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 946 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 947 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 948 break; 949 case PC_FIELDSPLIT_SCHUR_FACT_FULL: 950 /* [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 */ 951 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 952 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 953 ierr = PetscLogEventBegin(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 954 ierr = KSPSolve(kspLower,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 955 ierr = PetscLogEventEnd(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 956 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 957 ierr = VecScale(ilinkD->x,-1.0);CHKERRQ(ierr); 958 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 959 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 960 961 ierr = PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 962 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 963 ierr = PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 964 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 965 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 966 967 if (kspUpper == kspA) { 968 ierr = MatMult(jac->B,ilinkD->y,ilinkA->y);CHKERRQ(ierr); 969 ierr = VecAXPY(ilinkA->x,-1.0,ilinkA->y);CHKERRQ(ierr); 970 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 971 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 972 ierr = PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 973 } else { 974 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 975 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 976 ierr = PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 977 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 978 ierr = PetscLogEventBegin(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);CHKERRQ(ierr); 979 ierr = KSPSolve(kspUpper,ilinkA->x,ilinkA->z);CHKERRQ(ierr); 980 ierr = PetscLogEventEnd(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);CHKERRQ(ierr); 981 ierr = VecAXPY(ilinkA->y,-1.0,ilinkA->z);CHKERRQ(ierr); 982 } 983 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 984 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 985 } 986 PetscFunctionReturn(0); 987 } 988 989 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y) 990 { 991 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 992 PetscErrorCode ierr; 993 PC_FieldSplitLink ilink = jac->head; 994 PetscInt cnt,bs; 995 KSPConvergedReason reason; 996 997 PetscFunctionBegin; 998 if (jac->type == PC_COMPOSITE_ADDITIVE) { 999 if (jac->defaultsplit) { 1000 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 1001 if (jac->bs > 0 && bs != jac->bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Blocksize of x vector %D does not match fieldsplit blocksize %D",bs,jac->bs); 1002 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 1003 if (jac->bs > 0 && bs != jac->bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Blocksize of y vector %D does not match fieldsplit blocksize %D",bs,jac->bs); 1004 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 1005 while (ilink) { 1006 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1007 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1008 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1009 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1010 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1011 pc->failedreason = PC_SUBPC_ERROR; 1012 } 1013 ilink = ilink->next; 1014 } 1015 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 1016 } else { 1017 ierr = VecSet(y,0.0);CHKERRQ(ierr); 1018 while (ilink) { 1019 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 1020 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1021 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1022 pc->failedreason = PC_SUBPC_ERROR; 1023 } 1024 ilink = ilink->next; 1025 } 1026 } 1027 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) { 1028 ierr = VecSet(y,0.0);CHKERRQ(ierr); 1029 /* solve on first block for first block variables */ 1030 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1031 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1032 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1033 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1034 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1035 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1036 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1037 pc->failedreason = PC_SUBPC_ERROR; 1038 } 1039 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1040 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1041 1042 /* compute the residual only onto second block variables using first block variables */ 1043 ierr = MatMult(jac->Afield[1],ilink->y,ilink->next->x);CHKERRQ(ierr); 1044 ilink = ilink->next; 1045 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 1046 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1047 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1048 1049 /* solve on second block variables */ 1050 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1051 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1052 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1053 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1054 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1055 pc->failedreason = PC_SUBPC_ERROR; 1056 } 1057 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1058 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1059 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1060 if (!jac->w1) { 1061 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 1062 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 1063 } 1064 ierr = VecSet(y,0.0);CHKERRQ(ierr); 1065 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 1066 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1067 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1068 pc->failedreason = PC_SUBPC_ERROR; 1069 } 1070 cnt = 1; 1071 while (ilink->next) { 1072 ilink = ilink->next; 1073 /* compute the residual only over the part of the vector needed */ 1074 ierr = MatMult(jac->Afield[cnt++],y,ilink->x);CHKERRQ(ierr); 1075 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 1076 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1077 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1078 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1079 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1080 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1081 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1082 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1083 pc->failedreason = PC_SUBPC_ERROR; 1084 } 1085 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1086 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1087 } 1088 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1089 cnt -= 2; 1090 while (ilink->previous) { 1091 ilink = ilink->previous; 1092 /* compute the residual only over the part of the vector needed */ 1093 ierr = MatMult(jac->Afield[cnt--],y,ilink->x);CHKERRQ(ierr); 1094 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 1095 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1096 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1097 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1098 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1099 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1100 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1101 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1102 pc->failedreason = PC_SUBPC_ERROR; 1103 } 1104 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1105 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1106 } 1107 } 1108 } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type); 1109 PetscFunctionReturn(0); 1110 } 1111 1112 #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \ 1113 (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 1114 VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 1115 PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \ 1116 KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \ 1117 PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \ 1118 VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \ 1119 VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE)) 1120 1121 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y) 1122 { 1123 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1124 PetscErrorCode ierr; 1125 PC_FieldSplitLink ilink = jac->head; 1126 PetscInt bs; 1127 KSPConvergedReason reason; 1128 1129 PetscFunctionBegin; 1130 if (jac->type == PC_COMPOSITE_ADDITIVE) { 1131 if (jac->defaultsplit) { 1132 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 1133 if (jac->bs > 0 && bs != jac->bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Blocksize of x vector %D does not match fieldsplit blocksize %D",bs,jac->bs); 1134 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 1135 if (jac->bs > 0 && bs != jac->bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Blocksize of y vector %D does not match fieldsplit blocksize %D",bs,jac->bs); 1136 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 1137 while (ilink) { 1138 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1139 ierr = KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1140 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1141 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1142 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1143 pc->failedreason = PC_SUBPC_ERROR; 1144 } 1145 ilink = ilink->next; 1146 } 1147 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 1148 } else { 1149 ierr = VecSet(y,0.0);CHKERRQ(ierr); 1150 while (ilink) { 1151 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 1152 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1153 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1154 pc->failedreason = PC_SUBPC_ERROR; 1155 } 1156 ilink = ilink->next; 1157 } 1158 } 1159 } else { 1160 if (!jac->w1) { 1161 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 1162 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 1163 } 1164 ierr = VecSet(y,0.0);CHKERRQ(ierr); 1165 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1166 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 1167 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1168 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1169 pc->failedreason = PC_SUBPC_ERROR; 1170 } 1171 while (ilink->next) { 1172 ilink = ilink->next; 1173 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 1174 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 1175 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 1176 } 1177 while (ilink->previous) { 1178 ilink = ilink->previous; 1179 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 1180 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 1181 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 1182 } 1183 } else { 1184 while (ilink->next) { /* get to last entry in linked list */ 1185 ilink = ilink->next; 1186 } 1187 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 1188 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1189 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1190 pc->failedreason = PC_SUBPC_ERROR; 1191 } 1192 while (ilink->previous) { 1193 ilink = ilink->previous; 1194 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 1195 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 1196 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 1197 } 1198 } 1199 } 1200 PetscFunctionReturn(0); 1201 } 1202 1203 static PetscErrorCode PCReset_FieldSplit(PC pc) 1204 { 1205 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1206 PetscErrorCode ierr; 1207 PC_FieldSplitLink ilink = jac->head,next; 1208 1209 PetscFunctionBegin; 1210 while (ilink) { 1211 ierr = KSPDestroy(&ilink->ksp);CHKERRQ(ierr); 1212 ierr = VecDestroy(&ilink->x);CHKERRQ(ierr); 1213 ierr = VecDestroy(&ilink->y);CHKERRQ(ierr); 1214 ierr = VecDestroy(&ilink->z);CHKERRQ(ierr); 1215 ierr = VecScatterDestroy(&ilink->sctx);CHKERRQ(ierr); 1216 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1217 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1218 ierr = PetscFree(ilink->splitname);CHKERRQ(ierr); 1219 ierr = PetscFree(ilink->fields);CHKERRQ(ierr); 1220 ierr = PetscFree(ilink->fields_col);CHKERRQ(ierr); 1221 next = ilink->next; 1222 ierr = PetscFree(ilink);CHKERRQ(ierr); 1223 ilink = next; 1224 } 1225 jac->head = NULL; 1226 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 1227 if (jac->mat && jac->mat != jac->pmat) { 1228 ierr = MatDestroyMatrices(jac->nsplits,&jac->mat);CHKERRQ(ierr); 1229 } else if (jac->mat) { 1230 jac->mat = NULL; 1231 } 1232 if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);} 1233 if (jac->Afield) {ierr = MatDestroyMatrices(jac->nsplits,&jac->Afield);CHKERRQ(ierr);} 1234 jac->nsplits = 0; 1235 ierr = VecDestroy(&jac->w1);CHKERRQ(ierr); 1236 ierr = VecDestroy(&jac->w2);CHKERRQ(ierr); 1237 ierr = MatDestroy(&jac->schur);CHKERRQ(ierr); 1238 ierr = MatDestroy(&jac->schurp);CHKERRQ(ierr); 1239 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1240 ierr = KSPDestroy(&jac->kspschur);CHKERRQ(ierr); 1241 ierr = KSPDestroy(&jac->kspupper);CHKERRQ(ierr); 1242 ierr = MatDestroy(&jac->B);CHKERRQ(ierr); 1243 ierr = MatDestroy(&jac->C);CHKERRQ(ierr); 1244 jac->isrestrict = PETSC_FALSE; 1245 PetscFunctionReturn(0); 1246 } 1247 1248 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 1249 { 1250 PetscErrorCode ierr; 1251 1252 PetscFunctionBegin; 1253 ierr = PCReset_FieldSplit(pc);CHKERRQ(ierr); 1254 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1255 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);CHKERRQ(ierr); 1256 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);CHKERRQ(ierr); 1257 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);CHKERRQ(ierr); 1258 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);CHKERRQ(ierr); 1259 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);CHKERRQ(ierr); 1260 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL);CHKERRQ(ierr); 1261 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL);CHKERRQ(ierr); 1262 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);CHKERRQ(ierr); 1263 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",NULL);CHKERRQ(ierr); 1264 PetscFunctionReturn(0); 1265 } 1266 1267 static PetscErrorCode PCSetFromOptions_FieldSplit(PetscOptionItems *PetscOptionsObject,PC pc) 1268 { 1269 PetscErrorCode ierr; 1270 PetscInt bs; 1271 PetscBool flg; 1272 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1273 PCCompositeType ctype; 1274 1275 PetscFunctionBegin; 1276 ierr = PetscOptionsHead(PetscOptionsObject,"FieldSplit options");CHKERRQ(ierr); 1277 ierr = PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);CHKERRQ(ierr); 1278 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 1279 if (flg) { 1280 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 1281 } 1282 jac->diag_use_amat = pc->useAmat; 1283 ierr = PetscOptionsBool("-pc_fieldsplit_diag_use_amat","Use Amat (not Pmat) to extract diagonal fieldsplit blocks", "PCFieldSplitSetDiagUseAmat",jac->diag_use_amat,&jac->diag_use_amat,NULL);CHKERRQ(ierr); 1284 jac->offdiag_use_amat = pc->useAmat; 1285 ierr = PetscOptionsBool("-pc_fieldsplit_off_diag_use_amat","Use Amat (not Pmat) to extract off-diagonal fieldsplit blocks", "PCFieldSplitSetOffDiagUseAmat",jac->offdiag_use_amat,&jac->offdiag_use_amat,NULL);CHKERRQ(ierr); 1286 ierr = PetscOptionsBool("-pc_fieldsplit_detect_saddle_point","Form 2-way split by detecting zero diagonal entries", "PCFieldSplitSetDetectSaddlePoint",jac->detect,&jac->detect,NULL);CHKERRQ(ierr); 1287 ierr = PCFieldSplitSetDetectSaddlePoint(pc,jac->detect);CHKERRQ(ierr); /* Sets split type and Schur PC type */ 1288 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr); 1289 if (flg) { 1290 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr); 1291 } 1292 /* Only setup fields once */ 1293 if ((jac->bs > 0) && (jac->nsplits == 0)) { 1294 /* only allow user to set fields from command line if bs is already known. 1295 otherwise user can set them in PCFieldSplitSetDefaults() */ 1296 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 1297 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 1298 } 1299 if (jac->type == PC_COMPOSITE_SCHUR) { 1300 ierr = PetscOptionsGetEnum(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);CHKERRQ(ierr); 1301 if (flg) {ierr = PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");CHKERRQ(ierr);} 1302 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,NULL);CHKERRQ(ierr); 1303 ierr = PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);CHKERRQ(ierr); 1304 ierr = PetscOptionsScalar("-pc_fieldsplit_schur_scale","Scale Schur complement","PCFieldSplitSetSchurScale",jac->schurscale,&jac->schurscale,NULL);CHKERRQ(ierr); 1305 } 1306 ierr = PetscOptionsTail();CHKERRQ(ierr); 1307 PetscFunctionReturn(0); 1308 } 1309 1310 /*------------------------------------------------------------------------------------*/ 1311 1312 static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1313 { 1314 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1315 PetscErrorCode ierr; 1316 PC_FieldSplitLink ilink,next = jac->head; 1317 char prefix[128]; 1318 PetscInt i; 1319 1320 PetscFunctionBegin; 1321 if (jac->splitdefined) { 1322 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1323 PetscFunctionReturn(0); 1324 } 1325 for (i=0; i<n; i++) { 1326 if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 1327 if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 1328 } 1329 ierr = PetscNew(&ilink);CHKERRQ(ierr); 1330 if (splitname) { 1331 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1332 } else { 1333 ierr = PetscMalloc1(3,&ilink->splitname);CHKERRQ(ierr); 1334 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr); 1335 } 1336 ilink->event = jac->nsplits < 5 ? KSP_Solve_FS_0 + jac->nsplits : KSP_Solve_FS_0 + 4; /* Any split great than 4 gets logged in the 4th split */ 1337 ierr = PetscMalloc1(n,&ilink->fields);CHKERRQ(ierr); 1338 ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 1339 ierr = PetscMalloc1(n,&ilink->fields_col);CHKERRQ(ierr); 1340 ierr = PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));CHKERRQ(ierr); 1341 1342 ilink->nfields = n; 1343 ilink->next = NULL; 1344 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr); 1345 ierr = KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);CHKERRQ(ierr); 1346 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1347 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1348 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1349 1350 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr); 1351 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1352 1353 if (!next) { 1354 jac->head = ilink; 1355 ilink->previous = NULL; 1356 } else { 1357 while (next->next) { 1358 next = next->next; 1359 } 1360 next->next = ilink; 1361 ilink->previous = next; 1362 } 1363 jac->nsplits++; 1364 PetscFunctionReturn(0); 1365 } 1366 1367 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 1368 { 1369 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1370 PetscErrorCode ierr; 1371 1372 PetscFunctionBegin; 1373 if (!jac->schur) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitGetSubKSP()"); 1374 ierr = PetscMalloc1(jac->nsplits,subksp);CHKERRQ(ierr); 1375 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr); 1376 1377 (*subksp)[1] = jac->kspschur; 1378 if (n) *n = jac->nsplits; 1379 PetscFunctionReturn(0); 1380 } 1381 1382 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1383 { 1384 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1385 PetscErrorCode ierr; 1386 PetscInt cnt = 0; 1387 PC_FieldSplitLink ilink = jac->head; 1388 1389 PetscFunctionBegin; 1390 ierr = PetscMalloc1(jac->nsplits,subksp);CHKERRQ(ierr); 1391 while (ilink) { 1392 (*subksp)[cnt++] = ilink->ksp; 1393 ilink = ilink->next; 1394 } 1395 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); 1396 if (n) *n = jac->nsplits; 1397 PetscFunctionReturn(0); 1398 } 1399 1400 /*@C 1401 PCFieldSplitRestrictIS - Restricts the fieldsplit ISs to be within a given IS. 1402 1403 Input Parameters: 1404 + pc - the preconditioner context 1405 + is - the index set that defines the indices to which the fieldsplit is to be restricted 1406 1407 Level: advanced 1408 1409 @*/ 1410 PetscErrorCode PCFieldSplitRestrictIS(PC pc,IS isy) 1411 { 1412 PetscErrorCode ierr; 1413 1414 PetscFunctionBegin; 1415 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1416 PetscValidHeaderSpecific(isy,IS_CLASSID,2); 1417 ierr = PetscTryMethod(pc,"PCFieldSplitRestrictIS_C",(PC,IS),(pc,isy));CHKERRQ(ierr); 1418 PetscFunctionReturn(0); 1419 } 1420 1421 1422 static PetscErrorCode PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy) 1423 { 1424 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1425 PetscErrorCode ierr; 1426 PC_FieldSplitLink ilink = jac->head, next; 1427 PetscInt localsize,size,sizez,i; 1428 const PetscInt *ind, *indz; 1429 PetscInt *indc, *indcz; 1430 PetscBool flg; 1431 1432 PetscFunctionBegin; 1433 ierr = ISGetLocalSize(isy,&localsize);CHKERRQ(ierr); 1434 ierr = MPI_Scan(&localsize,&size,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)isy));CHKERRQ(ierr); 1435 size -= localsize; 1436 while(ilink) { 1437 IS isrl,isr; 1438 PC subpc; 1439 ierr = ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl);CHKERRQ(ierr); 1440 ierr = ISGetLocalSize(isrl,&localsize);CHKERRQ(ierr); 1441 ierr = PetscMalloc1(localsize,&indc);CHKERRQ(ierr); 1442 ierr = ISGetIndices(isrl,&ind);CHKERRQ(ierr); 1443 ierr = PetscMemcpy(indc,ind,localsize*sizeof(PetscInt));CHKERRQ(ierr); 1444 ierr = ISRestoreIndices(isrl,&ind);CHKERRQ(ierr); 1445 ierr = ISDestroy(&isrl);CHKERRQ(ierr); 1446 for (i=0; i<localsize; i++) *(indc+i) += size; 1447 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)isy),localsize,indc,PETSC_OWN_POINTER,&isr);CHKERRQ(ierr); 1448 ierr = PetscObjectReference((PetscObject)isr);CHKERRQ(ierr); 1449 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1450 ilink->is = isr; 1451 ierr = PetscObjectReference((PetscObject)isr);CHKERRQ(ierr); 1452 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1453 ilink->is_col = isr; 1454 ierr = ISDestroy(&isr);CHKERRQ(ierr); 1455 ierr = KSPGetPC(ilink->ksp, &subpc);CHKERRQ(ierr); 1456 ierr = PetscObjectTypeCompare((PetscObject)subpc,PCFIELDSPLIT,&flg);CHKERRQ(ierr); 1457 if(flg) { 1458 IS iszl,isz; 1459 MPI_Comm comm; 1460 ierr = ISGetLocalSize(ilink->is,&localsize);CHKERRQ(ierr); 1461 comm = PetscObjectComm((PetscObject)ilink->is); 1462 ierr = ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl);CHKERRQ(ierr); 1463 ierr = MPI_Scan(&localsize,&sizez,1,MPIU_INT,MPI_SUM,comm);CHKERRQ(ierr); 1464 sizez -= localsize; 1465 ierr = ISGetLocalSize(iszl,&localsize);CHKERRQ(ierr); 1466 ierr = PetscMalloc1(localsize,&indcz);CHKERRQ(ierr); 1467 ierr = ISGetIndices(iszl,&indz);CHKERRQ(ierr); 1468 ierr = PetscMemcpy(indcz,indz,localsize*sizeof(PetscInt));CHKERRQ(ierr); 1469 ierr = ISRestoreIndices(iszl,&indz);CHKERRQ(ierr); 1470 ierr = ISDestroy(&iszl);CHKERRQ(ierr); 1471 for (i=0; i<localsize; i++) *(indcz+i) += sizez; 1472 ierr = ISCreateGeneral(comm,localsize,indcz,PETSC_OWN_POINTER,&isz);CHKERRQ(ierr); 1473 ierr = PCFieldSplitRestrictIS(subpc,isz);CHKERRQ(ierr); 1474 ierr = ISDestroy(&isz);CHKERRQ(ierr); 1475 } 1476 next = ilink->next; 1477 ilink = next; 1478 } 1479 jac->isrestrict = PETSC_TRUE; 1480 PetscFunctionReturn(0); 1481 } 1482 1483 static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 1484 { 1485 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1486 PetscErrorCode ierr; 1487 PC_FieldSplitLink ilink, next = jac->head; 1488 char prefix[128]; 1489 1490 PetscFunctionBegin; 1491 if (jac->splitdefined) { 1492 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1493 PetscFunctionReturn(0); 1494 } 1495 ierr = PetscNew(&ilink);CHKERRQ(ierr); 1496 if (splitname) { 1497 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1498 } else { 1499 ierr = PetscMalloc1(8,&ilink->splitname);CHKERRQ(ierr); 1500 ierr = PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);CHKERRQ(ierr); 1501 } 1502 ilink->event = jac->nsplits < 5 ? KSP_Solve_FS_0 + jac->nsplits : KSP_Solve_FS_0 + 4; /* Any split great than 4 gets logged in the 4th split */ 1503 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1504 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1505 ilink->is = is; 1506 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1507 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1508 ilink->is_col = is; 1509 ilink->next = NULL; 1510 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr); 1511 ierr = KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);CHKERRQ(ierr); 1512 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1513 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1514 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1515 1516 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr); 1517 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1518 1519 if (!next) { 1520 jac->head = ilink; 1521 ilink->previous = NULL; 1522 } else { 1523 while (next->next) { 1524 next = next->next; 1525 } 1526 next->next = ilink; 1527 ilink->previous = next; 1528 } 1529 jac->nsplits++; 1530 PetscFunctionReturn(0); 1531 } 1532 1533 /*@C 1534 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 1535 1536 Logically Collective on PC 1537 1538 Input Parameters: 1539 + pc - the preconditioner context 1540 . splitname - name of this split, if NULL the number of the split is used 1541 . n - the number of fields in this split 1542 - fields - the fields in this split 1543 1544 Level: intermediate 1545 1546 Notes: 1547 Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 1548 1549 The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block 1550 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 1551 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 1552 where the numbered entries indicate what is in the field. 1553 1554 This function is called once per split (it creates a new split each time). Solve options 1555 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1556 1557 Developer Note: This routine does not actually create the IS representing the split, that is delayed 1558 until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be 1559 available when this routine is called. 1560 1561 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 1562 1563 @*/ 1564 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1565 { 1566 PetscErrorCode ierr; 1567 1568 PetscFunctionBegin; 1569 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1570 PetscValidCharPointer(splitname,2); 1571 if (n < 1) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 1572 PetscValidIntPointer(fields,3); 1573 ierr = PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt*,const PetscInt*),(pc,splitname,n,fields,fields_col));CHKERRQ(ierr); 1574 PetscFunctionReturn(0); 1575 } 1576 1577 /*@ 1578 PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) 1579 1580 Logically Collective on PC 1581 1582 Input Parameters: 1583 + pc - the preconditioner object 1584 - flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from 1585 1586 Options Database: 1587 . -pc_fieldsplit_diag_use_amat 1588 1589 Level: intermediate 1590 1591 .seealso: PCFieldSplitGetDiagUseAmat(), PCFieldSplitSetOffDiagUseAmat(), PCFIELDSPLIT 1592 1593 @*/ 1594 PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg) 1595 { 1596 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1597 PetscBool isfs; 1598 PetscErrorCode ierr; 1599 1600 PetscFunctionBegin; 1601 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1602 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1603 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1604 jac->diag_use_amat = flg; 1605 PetscFunctionReturn(0); 1606 } 1607 1608 /*@ 1609 PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) 1610 1611 Logically Collective on PC 1612 1613 Input Parameters: 1614 . pc - the preconditioner object 1615 1616 Output Parameters: 1617 . flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from 1618 1619 1620 Level: intermediate 1621 1622 .seealso: PCFieldSplitSetDiagUseAmat(), PCFieldSplitGetOffDiagUseAmat(), PCFIELDSPLIT 1623 1624 @*/ 1625 PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg) 1626 { 1627 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1628 PetscBool isfs; 1629 PetscErrorCode ierr; 1630 1631 PetscFunctionBegin; 1632 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1633 PetscValidPointer(flg,2); 1634 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1635 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1636 *flg = jac->diag_use_amat; 1637 PetscFunctionReturn(0); 1638 } 1639 1640 /*@ 1641 PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) 1642 1643 Logically Collective on PC 1644 1645 Input Parameters: 1646 + pc - the preconditioner object 1647 - flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 1648 1649 Options Database: 1650 . -pc_fieldsplit_off_diag_use_amat 1651 1652 Level: intermediate 1653 1654 .seealso: PCFieldSplitGetOffDiagUseAmat(), PCFieldSplitSetDiagUseAmat(), PCFIELDSPLIT 1655 1656 @*/ 1657 PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg) 1658 { 1659 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1660 PetscBool isfs; 1661 PetscErrorCode ierr; 1662 1663 PetscFunctionBegin; 1664 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1665 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1666 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1667 jac->offdiag_use_amat = flg; 1668 PetscFunctionReturn(0); 1669 } 1670 1671 /*@ 1672 PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) 1673 1674 Logically Collective on PC 1675 1676 Input Parameters: 1677 . pc - the preconditioner object 1678 1679 Output Parameters: 1680 . flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 1681 1682 1683 Level: intermediate 1684 1685 .seealso: PCFieldSplitSetOffDiagUseAmat(), PCFieldSplitGetDiagUseAmat(), PCFIELDSPLIT 1686 1687 @*/ 1688 PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg) 1689 { 1690 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1691 PetscBool isfs; 1692 PetscErrorCode ierr; 1693 1694 PetscFunctionBegin; 1695 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1696 PetscValidPointer(flg,2); 1697 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1698 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1699 *flg = jac->offdiag_use_amat; 1700 PetscFunctionReturn(0); 1701 } 1702 1703 1704 1705 /*@C 1706 PCFieldSplitSetIS - Sets the exact elements for field 1707 1708 Logically Collective on PC 1709 1710 Input Parameters: 1711 + pc - the preconditioner context 1712 . splitname - name of this split, if NULL the number of the split is used 1713 - is - the index set that defines the vector elements in this field 1714 1715 1716 Notes: 1717 Use PCFieldSplitSetFields(), for fields defined by strided types. 1718 1719 This function is called once per split (it creates a new split each time). Solve options 1720 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1721 1722 Level: intermediate 1723 1724 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 1725 1726 @*/ 1727 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 1728 { 1729 PetscErrorCode ierr; 1730 1731 PetscFunctionBegin; 1732 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1733 if (splitname) PetscValidCharPointer(splitname,2); 1734 PetscValidHeaderSpecific(is,IS_CLASSID,3); 1735 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr); 1736 PetscFunctionReturn(0); 1737 } 1738 1739 /*@C 1740 PCFieldSplitGetIS - Retrieves the elements for a field as an IS 1741 1742 Logically Collective on PC 1743 1744 Input Parameters: 1745 + pc - the preconditioner context 1746 - splitname - name of this split 1747 1748 Output Parameter: 1749 - is - the index set that defines the vector elements in this field, or NULL if the field is not found 1750 1751 Level: intermediate 1752 1753 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS() 1754 1755 @*/ 1756 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is) 1757 { 1758 PetscErrorCode ierr; 1759 1760 PetscFunctionBegin; 1761 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1762 PetscValidCharPointer(splitname,2); 1763 PetscValidPointer(is,3); 1764 { 1765 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 1766 PC_FieldSplitLink ilink = jac->head; 1767 PetscBool found; 1768 1769 *is = NULL; 1770 while (ilink) { 1771 ierr = PetscStrcmp(ilink->splitname, splitname, &found);CHKERRQ(ierr); 1772 if (found) { 1773 *is = ilink->is; 1774 break; 1775 } 1776 ilink = ilink->next; 1777 } 1778 } 1779 PetscFunctionReturn(0); 1780 } 1781 1782 /*@ 1783 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 1784 fieldsplit preconditioner. If not set the matrix block size is used. 1785 1786 Logically Collective on PC 1787 1788 Input Parameters: 1789 + pc - the preconditioner context 1790 - bs - the block size 1791 1792 Level: intermediate 1793 1794 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 1795 1796 @*/ 1797 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 1798 { 1799 PetscErrorCode ierr; 1800 1801 PetscFunctionBegin; 1802 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1803 PetscValidLogicalCollectiveInt(pc,bs,2); 1804 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr); 1805 PetscFunctionReturn(0); 1806 } 1807 1808 /*@C 1809 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1810 1811 Collective on KSP 1812 1813 Input Parameter: 1814 . pc - the preconditioner context 1815 1816 Output Parameters: 1817 + n - the number of splits 1818 - subksp - the array of KSP contexts 1819 1820 Note: 1821 After PCFieldSplitGetSubKSP() the array of KSPs is to be freed by the user with PetscFree() 1822 (not the KSP just the array that contains them). 1823 1824 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 1825 1826 Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs. 1827 You can call PCFieldSplitGetSubKSP(pc,n,PETSC_NULL_KSP,ierr) to determine how large the 1828 KSP array must be. 1829 1830 1831 Level: advanced 1832 1833 .seealso: PCFIELDSPLIT 1834 @*/ 1835 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1836 { 1837 PetscErrorCode ierr; 1838 1839 PetscFunctionBegin; 1840 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1841 if (n) PetscValidIntPointer(n,2); 1842 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 1843 PetscFunctionReturn(0); 1844 } 1845 1846 /*@ 1847 PCFieldSplitSetSchurPre - Indicates what operator is used to construct the preconditioner for the Schur complement. 1848 A11 matrix. Otherwise no preconditioner is used. 1849 1850 Collective on PC 1851 1852 Input Parameters: 1853 + pc - the preconditioner context 1854 . ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11 (default), PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_SCHUR_PRE_USER 1855 PC_FIELDSPLIT_SCHUR_PRE_SELFP, and PC_FIELDSPLIT_SCHUR_PRE_FULL 1856 - userpre - matrix to use for preconditioning, or NULL 1857 1858 Options Database: 1859 . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11. See notes for meaning of various arguments 1860 1861 Notes: 1862 $ If ptype is 1863 $ a11 then the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner 1864 $ matrix associated with the Schur complement (i.e. A11), not the Schur complement matrix 1865 $ self the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix: 1866 $ The only preconditioner that currently works with this symbolic respresentation matrix object is the PCLSC 1867 $ preconditioner 1868 $ user then the preconditioner for the Schur complement is generated from the user provided matrix (pre argument 1869 $ to this function). 1870 $ selfp then the preconditioning for the Schur complement is generated from an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01 1871 $ This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be 1872 $ lumped before extracting the diagonal using the additional option -fieldsplit_1_mat_schur_complement_ainv_type lump 1873 $ full then the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation computed internally by PFIELDSPLIT (this is expensive) 1874 $ useful mostly as a test that the Schur complement approach can work for your problem 1875 1876 When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense 1877 with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and 1878 -fieldsplit_1_pc_type lsc which uses the least squares commutator to compute a preconditioner for the Schur complement. 1879 1880 Level: intermediate 1881 1882 .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, 1883 MatSchurComplementSetAinvType(), PCLSC 1884 1885 @*/ 1886 PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1887 { 1888 PetscErrorCode ierr; 1889 1890 PetscFunctionBegin; 1891 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1892 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr); 1893 PetscFunctionReturn(0); 1894 } 1895 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */ 1896 1897 /*@ 1898 PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be 1899 preconditioned. See PCFieldSplitSetSchurPre() for details. 1900 1901 Logically Collective on PC 1902 1903 Input Parameters: 1904 . pc - the preconditioner context 1905 1906 Output Parameters: 1907 + ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11, PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_PRE_USER 1908 - userpre - matrix to use for preconditioning (with PC_FIELDSPLIT_PRE_USER), or NULL 1909 1910 Level: intermediate 1911 1912 .seealso: PCFieldSplitSetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC 1913 1914 @*/ 1915 PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre) 1916 { 1917 PetscErrorCode ierr; 1918 1919 PetscFunctionBegin; 1920 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1921 ierr = PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));CHKERRQ(ierr); 1922 PetscFunctionReturn(0); 1923 } 1924 1925 /*@ 1926 PCFieldSplitSchurGetS - extract the MatSchurComplement object used by this PC in case it needs to be configured separately 1927 1928 Not collective 1929 1930 Input Parameter: 1931 . pc - the preconditioner context 1932 1933 Output Parameter: 1934 . S - the Schur complement matrix 1935 1936 Notes: 1937 This matrix should not be destroyed using MatDestroy(); rather, use PCFieldSplitSchurRestoreS(). 1938 1939 Level: advanced 1940 1941 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurRestoreS() 1942 1943 @*/ 1944 PetscErrorCode PCFieldSplitSchurGetS(PC pc,Mat *S) 1945 { 1946 PetscErrorCode ierr; 1947 const char* t; 1948 PetscBool isfs; 1949 PC_FieldSplit *jac; 1950 1951 PetscFunctionBegin; 1952 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1953 ierr = PetscObjectGetType((PetscObject)pc,&t);CHKERRQ(ierr); 1954 ierr = PetscStrcmp(t,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1955 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 1956 jac = (PC_FieldSplit*)pc->data; 1957 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type); 1958 if (S) *S = jac->schur; 1959 PetscFunctionReturn(0); 1960 } 1961 1962 /*@ 1963 PCFieldSplitSchurRestoreS - restores the MatSchurComplement object used by this PC 1964 1965 Not collective 1966 1967 Input Parameters: 1968 + pc - the preconditioner context 1969 . S - the Schur complement matrix 1970 1971 Level: advanced 1972 1973 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurGetS() 1974 1975 @*/ 1976 PetscErrorCode PCFieldSplitSchurRestoreS(PC pc,Mat *S) 1977 { 1978 PetscErrorCode ierr; 1979 const char* t; 1980 PetscBool isfs; 1981 PC_FieldSplit *jac; 1982 1983 PetscFunctionBegin; 1984 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1985 ierr = PetscObjectGetType((PetscObject)pc,&t);CHKERRQ(ierr); 1986 ierr = PetscStrcmp(t,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1987 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 1988 jac = (PC_FieldSplit*)pc->data; 1989 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type); 1990 if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten"); 1991 PetscFunctionReturn(0); 1992 } 1993 1994 1995 static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1996 { 1997 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1998 PetscErrorCode ierr; 1999 2000 PetscFunctionBegin; 2001 jac->schurpre = ptype; 2002 if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) { 2003 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 2004 jac->schur_user = pre; 2005 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 2006 } 2007 PetscFunctionReturn(0); 2008 } 2009 2010 static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre) 2011 { 2012 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2013 2014 PetscFunctionBegin; 2015 *ptype = jac->schurpre; 2016 *pre = jac->schur_user; 2017 PetscFunctionReturn(0); 2018 } 2019 2020 /*@ 2021 PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain in the preconditioner 2022 2023 Collective on PC 2024 2025 Input Parameters: 2026 + pc - the preconditioner context 2027 - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default 2028 2029 Options Database: 2030 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full 2031 2032 2033 Level: intermediate 2034 2035 Notes: 2036 The FULL factorization is 2037 2038 $ (A B) = (1 0) (A 0) (1 Ainv*B) = L D U 2039 $ (C E) (C*Ainv 1) (0 S) (0 1 ) 2040 2041 where S = E - 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, 2042 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). Sign flipping of S can be turned off with PCFieldSplitSetSchurScale(). 2043 2044 $ If A and S are solved exactly 2045 $ *) FULL factorization is a direct solver. 2046 $ *) The preconditioned operator with LOWER or UPPER has all eigenvalues equal to 1 and minimal polynomial of degree 2, so KSPGMRES converges in 2 iterations. 2047 $ *) With DIAG, the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations. 2048 2049 If the iteration count is very low, consider using KSPFGMRES or KSPGCR which can use one less preconditioner 2050 application in this case. Note that the preconditioned operator may be highly non-normal, so such fast convergence may not be observed in practice. 2051 2052 For symmetric problems in which A is positive definite and S is negative definite, DIAG can be used with KSPMINRES. 2053 2054 Note that a flexible method like KSPFGMRES or KSPGCR must be used if the fieldsplit preconditioner is nonlinear (e.g. a few iterations of a Krylov method is used to solve with A or S). 2055 2056 References: 2057 + 1. - Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000). 2058 - 2. - Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001). 2059 2060 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCFieldSplitSetSchurScale() 2061 @*/ 2062 PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype) 2063 { 2064 PetscErrorCode ierr; 2065 2066 PetscFunctionBegin; 2067 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2068 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));CHKERRQ(ierr); 2069 PetscFunctionReturn(0); 2070 } 2071 2072 static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype) 2073 { 2074 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2075 2076 PetscFunctionBegin; 2077 jac->schurfactorization = ftype; 2078 PetscFunctionReturn(0); 2079 } 2080 2081 /*@ 2082 PCFieldSplitSetSchurScale - Controls the sign flip of S for PC_FIELDSPLIT_SCHUR_FACT_DIAG. 2083 2084 Collective on PC 2085 2086 Input Parameters: 2087 + pc - the preconditioner context 2088 - scale - scaling factor for the Schur complement 2089 2090 Options Database: 2091 . -pc_fieldsplit_schur_scale - default is -1.0 2092 2093 Level: intermediate 2094 2095 .seealso: PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurFactType, PCFieldSplitSetSchurScale() 2096 @*/ 2097 PetscErrorCode PCFieldSplitSetSchurScale(PC pc,PetscScalar scale) 2098 { 2099 PetscErrorCode ierr; 2100 2101 PetscFunctionBegin; 2102 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2103 PetscValidLogicalCollectiveScalar(pc,scale,2); 2104 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurScale_C",(PC,PetscScalar),(pc,scale));CHKERRQ(ierr); 2105 PetscFunctionReturn(0); 2106 } 2107 2108 static PetscErrorCode PCFieldSplitSetSchurScale_FieldSplit(PC pc,PetscScalar scale) 2109 { 2110 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2111 2112 PetscFunctionBegin; 2113 jac->schurscale = scale; 2114 PetscFunctionReturn(0); 2115 } 2116 2117 /*@C 2118 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 2119 2120 Collective on KSP 2121 2122 Input Parameter: 2123 . pc - the preconditioner context 2124 2125 Output Parameters: 2126 + A00 - the (0,0) block 2127 . A01 - the (0,1) block 2128 . A10 - the (1,0) block 2129 - A11 - the (1,1) block 2130 2131 Level: advanced 2132 2133 .seealso: PCFIELDSPLIT 2134 @*/ 2135 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 2136 { 2137 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2138 2139 PetscFunctionBegin; 2140 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2141 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 2142 if (A00) *A00 = jac->pmat[0]; 2143 if (A01) *A01 = jac->B; 2144 if (A10) *A10 = jac->C; 2145 if (A11) *A11 = jac->pmat[1]; 2146 PetscFunctionReturn(0); 2147 } 2148 2149 static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 2150 { 2151 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2152 PetscErrorCode ierr; 2153 2154 PetscFunctionBegin; 2155 jac->type = type; 2156 if (type == PC_COMPOSITE_SCHUR) { 2157 pc->ops->apply = PCApply_FieldSplit_Schur; 2158 pc->ops->view = PCView_FieldSplit_Schur; 2159 2160 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 2161 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);CHKERRQ(ierr); 2162 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);CHKERRQ(ierr); 2163 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);CHKERRQ(ierr); 2164 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",PCFieldSplitSetSchurScale_FieldSplit);CHKERRQ(ierr); 2165 2166 } else { 2167 pc->ops->apply = PCApply_FieldSplit; 2168 pc->ops->view = PCView_FieldSplit; 2169 2170 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 2171 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);CHKERRQ(ierr); 2172 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);CHKERRQ(ierr); 2173 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);CHKERRQ(ierr); 2174 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",0);CHKERRQ(ierr); 2175 } 2176 PetscFunctionReturn(0); 2177 } 2178 2179 static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 2180 { 2181 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2182 2183 PetscFunctionBegin; 2184 if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 2185 if (jac->bs > 0 && jac->bs != bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot change fieldsplit blocksize from %D to %D after it has been set",jac->bs,bs); 2186 jac->bs = bs; 2187 PetscFunctionReturn(0); 2188 } 2189 2190 /*@ 2191 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 2192 2193 Collective on PC 2194 2195 Input Parameter: 2196 . pc - the preconditioner context 2197 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 2198 2199 Options Database Key: 2200 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 2201 2202 Level: Intermediate 2203 2204 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 2205 2206 .seealso: PCCompositeSetType() 2207 2208 @*/ 2209 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 2210 { 2211 PetscErrorCode ierr; 2212 2213 PetscFunctionBegin; 2214 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2215 ierr = PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));CHKERRQ(ierr); 2216 PetscFunctionReturn(0); 2217 } 2218 2219 /*@ 2220 PCFieldSplitGetType - Gets the type of fieldsplit preconditioner. 2221 2222 Not collective 2223 2224 Input Parameter: 2225 . pc - the preconditioner context 2226 2227 Output Parameter: 2228 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 2229 2230 Level: Intermediate 2231 2232 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 2233 .seealso: PCCompositeSetType() 2234 @*/ 2235 PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type) 2236 { 2237 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2238 2239 PetscFunctionBegin; 2240 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2241 PetscValidIntPointer(type,2); 2242 *type = jac->type; 2243 PetscFunctionReturn(0); 2244 } 2245 2246 /*@ 2247 PCFieldSplitSetDMSplits - Flags whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 2248 2249 Logically Collective 2250 2251 Input Parameters: 2252 + pc - the preconditioner context 2253 - flg - boolean indicating whether to use field splits defined by the DM 2254 2255 Options Database Key: 2256 . -pc_fieldsplit_dm_splits 2257 2258 Level: Intermediate 2259 2260 .keywords: PC, DM, composite preconditioner, additive, multiplicative 2261 2262 .seealso: PCFieldSplitGetDMSplits() 2263 2264 @*/ 2265 PetscErrorCode PCFieldSplitSetDMSplits(PC pc,PetscBool flg) 2266 { 2267 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2268 PetscBool isfs; 2269 PetscErrorCode ierr; 2270 2271 PetscFunctionBegin; 2272 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2273 PetscValidLogicalCollectiveBool(pc,flg,2); 2274 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 2275 if (isfs) { 2276 jac->dm_splits = flg; 2277 } 2278 PetscFunctionReturn(0); 2279 } 2280 2281 2282 /*@ 2283 PCFieldSplitGetDMSplits - Returns flag indicating whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 2284 2285 Logically Collective 2286 2287 Input Parameter: 2288 . pc - the preconditioner context 2289 2290 Output Parameter: 2291 . flg - boolean indicating whether to use field splits defined by the DM 2292 2293 Level: Intermediate 2294 2295 .keywords: PC, DM, composite preconditioner, additive, multiplicative 2296 2297 .seealso: PCFieldSplitSetDMSplits() 2298 2299 @*/ 2300 PetscErrorCode PCFieldSplitGetDMSplits(PC pc,PetscBool* flg) 2301 { 2302 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2303 PetscBool isfs; 2304 PetscErrorCode ierr; 2305 2306 PetscFunctionBegin; 2307 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2308 PetscValidPointer(flg,2); 2309 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 2310 if (isfs) { 2311 if(flg) *flg = jac->dm_splits; 2312 } 2313 PetscFunctionReturn(0); 2314 } 2315 2316 /*@ 2317 PCFieldSplitGetDetectSaddlePoint - Returns flag indicating whether PCFieldSplit will attempt to automatically determine fields based on zero diagonal entries. 2318 2319 Logically Collective 2320 2321 Input Parameter: 2322 . pc - the preconditioner context 2323 2324 Output Parameter: 2325 . flg - boolean indicating whether to detect fields or not 2326 2327 Level: Intermediate 2328 2329 .seealso: PCFIELDSPLIT, PCFieldSplitSetDetectSaddlePoint() 2330 2331 @*/ 2332 PetscErrorCode PCFieldSplitGetDetectSaddlePoint(PC pc,PetscBool *flg) 2333 { 2334 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2335 2336 PetscFunctionBegin; 2337 *flg = jac->detect; 2338 PetscFunctionReturn(0); 2339 } 2340 2341 /*@ 2342 PCFieldSplitSetDetectSaddlePoint - Sets flag indicating whether PCFieldSplit will attempt to automatically determine fields based on zero diagonal entries. 2343 2344 Logically Collective 2345 2346 Notes: 2347 Also sets the split type to PC_COMPOSITE_SCHUR (see PCFieldSplitSetType()) and the Schur preconditioner type to PC_FIELDSPLIT_SCHUR_PRE_SELF (see PCFieldSplitSetSchurPre()). 2348 2349 Input Parameter: 2350 . pc - the preconditioner context 2351 2352 Output Parameter: 2353 . flg - boolean indicating whether to detect fields or not 2354 2355 Options Database Key: 2356 . -pc_fieldsplit_detect_saddle_point 2357 2358 Level: Intermediate 2359 2360 .seealso: PCFIELDSPLIT, PCFieldSplitSetDetectSaddlePoint(), PCFieldSplitSetType(), PCFieldSplitSetSchurPre() 2361 2362 @*/ 2363 PetscErrorCode PCFieldSplitSetDetectSaddlePoint(PC pc,PetscBool flg) 2364 { 2365 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2366 PetscErrorCode ierr; 2367 2368 PetscFunctionBegin; 2369 jac->detect = flg; 2370 if (jac->detect) { 2371 ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);CHKERRQ(ierr); 2372 ierr = PCFieldSplitSetSchurPre(pc,PC_FIELDSPLIT_SCHUR_PRE_SELF,NULL);CHKERRQ(ierr); 2373 } 2374 PetscFunctionReturn(0); 2375 } 2376 2377 /* -------------------------------------------------------------------------------------*/ 2378 /*MC 2379 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 2380 fields or groups of fields. See the users manual section "Solving Block Matrices" for more details. 2381 2382 To set options on the solvers for each block append -fieldsplit_ to all the PC 2383 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 2384 2385 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 2386 and set the options directly on the resulting KSP object 2387 2388 Level: intermediate 2389 2390 Options Database Keys: 2391 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 2392 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 2393 been supplied explicitly by -pc_fieldsplit_%d_fields 2394 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 2395 . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur> - type of relaxation or factorization splitting 2396 . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11; see PCFieldSplitSetSchurPre() 2397 . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero diagonal and uses Schur complement with no preconditioner as the solver 2398 2399 - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ 2400 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 2401 2402 Notes: 2403 Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 2404 to define a field by an arbitrary collection of entries. 2405 2406 If no fields are set the default is used. The fields are defined by entries strided by bs, 2407 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 2408 if this is not called the block size defaults to the blocksize of the second matrix passed 2409 to KSPSetOperators()/PCSetOperators(). 2410 2411 $ For the Schur complement preconditioner if J = ( A00 A01 ) 2412 $ ( A10 A11 ) 2413 $ the preconditioner using full factorization is 2414 $ ( I -ksp(A00) A01 ) ( inv(A00) 0 ) ( I 0 ) 2415 $ ( 0 I ) ( 0 ksp(S) ) ( -A10 ksp(A00) I ) 2416 where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. S is the Schur complement 2417 $ S = A11 - A10 ksp(A00) A01 2418 which is usually dense and not stored explicitly. The action of ksp(S) is computed using the KSP solver with prefix -fieldsplit_splitname_ (where splitname was given 2419 in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0, 2420 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 2421 A11 is used to construct a preconditioner for S, use PCFieldSplitSetSchurPre() for all the possible ways to construct the preconditioner for S. 2422 2423 The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above, 2424 diag gives 2425 $ ( inv(A00) 0 ) 2426 $ ( 0 -ksp(S) ) 2427 note that slightly counter intuitively there is a negative in front of the ksp(S) so that the preconditioner is positive definite. For SPD matrices J, the sign flip 2428 can be turned off with PCFieldSplitSetSchurScale() or by command line -pc_fieldsplit_schur_scale 1.0. The lower factorization is the inverse of 2429 $ ( A00 0 ) 2430 $ ( A10 S ) 2431 where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of 2432 $ ( A00 A01 ) 2433 $ ( 0 S ) 2434 where again the inverses of A00 and S are applied using KSPs. 2435 2436 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 2437 is used automatically for a second block. 2438 2439 The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1. 2440 Generally it should be used with the AIJ format. 2441 2442 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 2443 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT 2444 inside a smoother resulting in "Distributive Smoothers". 2445 2446 Concepts: physics based preconditioners, block preconditioners 2447 2448 There is a nice discussion of block preconditioners in 2449 2450 [El08] A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations 2451 Howard Elman, V.E. Howle, John Shadid, Robert Shuttleworth, Ray Tuminaro, Journal of Computational Physics 227 (2008) 1790--1808 2452 http://chess.cs.umd.edu/~elman/papers/tax.pdf 2453 2454 The Constrained Pressure Preconditioner (CPR) can be implemented using PCCOMPOSITE with PCGALERKIN. CPR first solves an R A P subsystem, updates the 2455 residual on all variables (PCCompositeSetType(pc,PC_COMPOSITE_MULTIPLICATIVE)), and then applies a simple ILU like preconditioner on all the variables. 2456 2457 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC, 2458 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre(), 2459 MatSchurComplementSetAinvType(), PCFieldSplitSetSchurScale(), 2460 PCFieldSplitSetDetectSaddlePoint() 2461 M*/ 2462 2463 PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc) 2464 { 2465 PetscErrorCode ierr; 2466 PC_FieldSplit *jac; 2467 2468 PetscFunctionBegin; 2469 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr); 2470 2471 jac->bs = -1; 2472 jac->nsplits = 0; 2473 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 2474 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 2475 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL; 2476 jac->schurscale = -1.0; 2477 jac->dm_splits = PETSC_TRUE; 2478 jac->detect = PETSC_FALSE; 2479 2480 pc->data = (void*)jac; 2481 2482 pc->ops->apply = PCApply_FieldSplit; 2483 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 2484 pc->ops->setup = PCSetUp_FieldSplit; 2485 pc->ops->reset = PCReset_FieldSplit; 2486 pc->ops->destroy = PCDestroy_FieldSplit; 2487 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 2488 pc->ops->view = PCView_FieldSplit; 2489 pc->ops->applyrichardson = 0; 2490 2491 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 2492 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 2493 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 2494 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 2495 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 2496 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",PCFieldSplitRestrictIS_FieldSplit);CHKERRQ(ierr); 2497 PetscFunctionReturn(0); 2498 } 2499