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