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