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