1 2 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/ 3 #include <petscdm.h> 4 5 6 const char *const PCFieldSplitSchurPreTypes[] = {"SELF","SELFP","A11","USER","FULL","PCFieldSplitSchurPreType","PC_FIELDSPLIT_SCHUR_PRE_",0}; 7 const char *const PCFieldSplitSchurFactTypes[] = {"DIAG","LOWER","UPPER","FULL","PCFieldSplitSchurFactType","PC_FIELDSPLIT_SCHUR_FACT_",0}; 8 9 typedef struct _PC_FieldSplitLink *PC_FieldSplitLink; 10 struct _PC_FieldSplitLink { 11 KSP ksp; 12 Vec x,y,z; 13 char *splitname; 14 PetscInt nfields; 15 PetscInt *fields,*fields_col; 16 VecScatter sctx; 17 IS is,is_col; 18 PC_FieldSplitLink next,previous; 19 }; 20 21 typedef struct { 22 PCCompositeType type; 23 PetscBool defaultsplit; /* Flag for a system with a set of 'k' scalar fields with the same layout (and bs = k) */ 24 PetscBool splitdefined; /* Flag is set after the splits have been defined, to prevent more splits from being added */ 25 PetscInt bs; /* Block size for IS and Mat structures */ 26 PetscInt nsplits; /* Number of field divisions defined */ 27 Vec *x,*y,w1,w2; 28 Mat *mat; /* The diagonal block for each split */ 29 Mat *pmat; /* The preconditioning diagonal block for each split */ 30 Mat *Afield; /* The rows of the matrix associated with each split */ 31 PetscBool issetup; 32 33 /* Only used when Schur complement preconditioning is used */ 34 Mat B; /* The (0,1) block */ 35 Mat C; /* The (1,0) block */ 36 Mat schur; /* The Schur complement S = A11 - A10 A00^{-1} A01, the KSP here, kspinner, is H_1 in [El08] */ 37 Mat schurp; /* Assembled approximation to S built by MatSchurComplement to be used as a preconditioning matrix when solving with S */ 38 Mat schur_user; /* User-provided preconditioning matrix for the Schur complement */ 39 PCFieldSplitSchurPreType schurpre; /* Determines which preconditioning matrix is used for the Schur complement */ 40 PCFieldSplitSchurFactType schurfactorization; 41 KSP kspschur; /* The solver for S */ 42 KSP kspupper; /* The solver for A in the upper diagonal part of the factorization (H_2 in [El08]) */ 43 PC_FieldSplitLink head; 44 PetscBool reset; /* indicates PCReset() has been last called on this object, hack */ 45 PetscBool suboptionsset; /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */ 46 PetscBool dm_splits; /* Whether to use DMCreateFieldDecomposition() whenever possible */ 47 PetscBool diag_use_amat; /* Whether to extract diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */ 48 PetscBool offdiag_use_amat; /* Whether to extract off-diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */ 49 } PC_FieldSplit; 50 51 /* 52 Notes: there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of 53 inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the 54 PC you could change this. 55 */ 56 57 /* This helper is so that setting a user-provided preconditioning matrix is orthogonal to choosing to use it. This way the 58 * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */ 59 static Mat FieldSplitSchurPre(PC_FieldSplit *jac) 60 { 61 switch (jac->schurpre) { 62 case PC_FIELDSPLIT_SCHUR_PRE_SELF: return jac->schur; 63 case PC_FIELDSPLIT_SCHUR_PRE_SELFP: return jac->schurp; 64 case PC_FIELDSPLIT_SCHUR_PRE_A11: return jac->pmat[1]; 65 case PC_FIELDSPLIT_SCHUR_PRE_FULL: /* We calculate this and store it in schur_user */ 66 case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */ 67 default: 68 return jac->schur_user ? jac->schur_user : jac->pmat[1]; 69 } 70 } 71 72 73 #include <petscdraw.h> 74 #undef __FUNCT__ 75 #define __FUNCT__ "PCView_FieldSplit" 76 static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer) 77 { 78 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 79 PetscErrorCode ierr; 80 PetscBool iascii,isdraw; 81 PetscInt i,j; 82 PC_FieldSplitLink ilink = jac->head; 83 84 PetscFunctionBegin; 85 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 86 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 87 if (iascii) { 88 if (jac->bs > 0) { 89 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);CHKERRQ(ierr); 90 } else { 91 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D\n",PCCompositeTypes[jac->type],jac->nsplits);CHKERRQ(ierr); 92 } 93 if (pc->useAmat) { 94 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");CHKERRQ(ierr); 95 } 96 if (jac->diag_use_amat) { 97 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for diagonal blocks\n");CHKERRQ(ierr); 98 } 99 if (jac->offdiag_use_amat) { 100 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for off-diagonal blocks\n");CHKERRQ(ierr); 101 } 102 ierr = PetscViewerASCIIPrintf(viewer," Solver info for each split is in the following KSP objects:\n");CHKERRQ(ierr); 103 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 104 for (i=0; i<jac->nsplits; i++) { 105 if (ilink->fields) { 106 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr); 107 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 108 for (j=0; j<ilink->nfields; j++) { 109 if (j > 0) { 110 ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr); 111 } 112 ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr); 113 } 114 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 115 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 116 } else { 117 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);CHKERRQ(ierr); 118 } 119 ierr = KSPView(ilink->ksp,viewer);CHKERRQ(ierr); 120 ilink = ilink->next; 121 } 122 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 123 } 124 125 if (isdraw) { 126 PetscDraw draw; 127 PetscReal x,y,w,wd; 128 129 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 130 ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr); 131 w = 2*PetscMin(1.0 - x,x); 132 wd = w/(jac->nsplits + 1); 133 x = x - wd*(jac->nsplits-1)/2.0; 134 for (i=0; i<jac->nsplits; i++) { 135 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr); 136 ierr = KSPView(ilink->ksp,viewer);CHKERRQ(ierr); 137 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 138 x += wd; 139 ilink = ilink->next; 140 } 141 } 142 PetscFunctionReturn(0); 143 } 144 145 #undef __FUNCT__ 146 #define __FUNCT__ "PCView_FieldSplit_Schur" 147 static PetscErrorCode PCView_FieldSplit_Schur(PC pc,PetscViewer viewer) 148 { 149 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 150 PetscErrorCode ierr; 151 PetscBool iascii,isdraw; 152 PetscInt i,j; 153 PC_FieldSplitLink ilink = jac->head; 154 155 PetscFunctionBegin; 156 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 157 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 158 if (iascii) { 159 if (jac->bs > 0) { 160 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, blocksize = %D, factorization %s\n",jac->bs,PCFieldSplitSchurFactTypes[jac->schurfactorization]);CHKERRQ(ierr); 161 } else { 162 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, factorization %s\n",PCFieldSplitSchurFactTypes[jac->schurfactorization]);CHKERRQ(ierr); 163 } 164 if (pc->useAmat) { 165 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");CHKERRQ(ierr); 166 } 167 switch (jac->schurpre) { 168 case PC_FIELDSPLIT_SCHUR_PRE_SELF: 169 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from S itself\n");CHKERRQ(ierr);break; 170 case PC_FIELDSPLIT_SCHUR_PRE_SELFP: 171 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; 172 case PC_FIELDSPLIT_SCHUR_PRE_A11: 173 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");CHKERRQ(ierr);break; 174 case PC_FIELDSPLIT_SCHUR_PRE_FULL: 175 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from the exact Schur complement\n");CHKERRQ(ierr);break; 176 case PC_FIELDSPLIT_SCHUR_PRE_USER: 177 if (jac->schur_user) { 178 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from user provided matrix\n");CHKERRQ(ierr); 179 } else { 180 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");CHKERRQ(ierr); 181 } 182 break; 183 default: 184 SETERRQ1(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Invalid Schur preconditioning type: %d", jac->schurpre); 185 } 186 ierr = PetscViewerASCIIPrintf(viewer," Split info:\n");CHKERRQ(ierr); 187 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 188 for (i=0; i<jac->nsplits; i++) { 189 if (ilink->fields) { 190 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr); 191 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 192 for (j=0; j<ilink->nfields; j++) { 193 if (j > 0) { 194 ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr); 195 } 196 ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr); 197 } 198 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 199 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 200 } else { 201 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);CHKERRQ(ierr); 202 } 203 ilink = ilink->next; 204 } 205 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for A00 block\n");CHKERRQ(ierr); 206 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 207 if (jac->head) { 208 ierr = KSPView(jac->head->ksp,viewer);CHKERRQ(ierr); 209 } else {ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr);} 210 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 211 if (jac->head && jac->kspupper != jac->head->ksp) { 212 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for upper A00 in upper triangular factor \n");CHKERRQ(ierr); 213 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 214 if (jac->kspupper) {ierr = KSPView(jac->kspupper,viewer);CHKERRQ(ierr);} 215 else {ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr);} 216 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 217 } 218 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for S = A11 - A10 inv(A00) A01 \n");CHKERRQ(ierr); 219 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 220 if (jac->kspschur) { 221 ierr = KSPView(jac->kspschur,viewer);CHKERRQ(ierr); 222 } else { 223 ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr); 224 } 225 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 226 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 227 } else if (isdraw && jac->head) { 228 PetscDraw draw; 229 PetscReal x,y,w,wd,h; 230 PetscInt cnt = 2; 231 char str[32]; 232 233 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 234 ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr); 235 if (jac->kspupper != jac->head->ksp) cnt++; 236 w = 2*PetscMin(1.0 - x,x); 237 wd = w/(cnt + 1); 238 239 ierr = PetscSNPrintf(str,32,"Schur fact. %s",PCFieldSplitSchurFactTypes[jac->schurfactorization]);CHKERRQ(ierr); 240 ierr = PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);CHKERRQ(ierr); 241 y -= h; 242 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_USER && !jac->schur_user) { 243 ierr = PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[PC_FIELDSPLIT_SCHUR_PRE_A11]);CHKERRQ(ierr); 244 } else { 245 ierr = PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[jac->schurpre]);CHKERRQ(ierr); 246 } 247 ierr = PetscDrawStringBoxed(draw,x+wd*(cnt-1)/2.0,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);CHKERRQ(ierr); 248 y -= h; 249 x = x - wd*(cnt-1)/2.0; 250 251 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr); 252 ierr = KSPView(jac->head->ksp,viewer);CHKERRQ(ierr); 253 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 254 if (jac->kspupper != jac->head->ksp) { 255 x += wd; 256 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr); 257 ierr = KSPView(jac->kspupper,viewer);CHKERRQ(ierr); 258 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 259 } 260 x += wd; 261 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr); 262 ierr = KSPView(jac->kspschur,viewer);CHKERRQ(ierr); 263 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 264 } 265 PetscFunctionReturn(0); 266 } 267 268 #undef __FUNCT__ 269 #define __FUNCT__ "PCFieldSplitSetRuntimeSplits_Private" 270 /* Precondition: jac->bs is set to a meaningful value */ 271 static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc) 272 { 273 PetscErrorCode ierr; 274 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 275 PetscInt i,nfields,*ifields,nfields_col,*ifields_col; 276 PetscBool flg,flg_col; 277 char optionname[128],splitname[8],optionname_col[128]; 278 279 PetscFunctionBegin; 280 ierr = PetscMalloc1(jac->bs,&ifields);CHKERRQ(ierr); 281 ierr = PetscMalloc1(jac->bs,&ifields_col);CHKERRQ(ierr); 282 for (i=0,flg=PETSC_TRUE;; i++) { 283 ierr = PetscSNPrintf(splitname,sizeof(splitname),"%D",i);CHKERRQ(ierr); 284 ierr = PetscSNPrintf(optionname,sizeof(optionname),"-pc_fieldsplit_%D_fields",i);CHKERRQ(ierr); 285 ierr = PetscSNPrintf(optionname_col,sizeof(optionname_col),"-pc_fieldsplit_%D_fields_col",i);CHKERRQ(ierr); 286 nfields = jac->bs; 287 nfields_col = jac->bs; 288 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname,ifields,&nfields,&flg);CHKERRQ(ierr); 289 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname_col,ifields_col,&nfields_col,&flg_col);CHKERRQ(ierr); 290 if (!flg) break; 291 else if (flg && !flg_col) { 292 if (!nfields) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields"); 293 ierr = PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields);CHKERRQ(ierr); 294 } else { 295 if (!nfields || !nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields"); 296 if (nfields != nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Number of row and column fields must match"); 297 ierr = PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields_col);CHKERRQ(ierr); 298 } 299 } 300 if (i > 0) { 301 /* Makes command-line setting of splits take precedence over setting them in code. 302 Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would 303 create new splits, which would probably not be what the user wanted. */ 304 jac->splitdefined = PETSC_TRUE; 305 } 306 ierr = PetscFree(ifields);CHKERRQ(ierr); 307 ierr = PetscFree(ifields_col);CHKERRQ(ierr); 308 PetscFunctionReturn(0); 309 } 310 311 #undef __FUNCT__ 312 #define __FUNCT__ "PCFieldSplitSetDefaults" 313 static PetscErrorCode PCFieldSplitSetDefaults(PC pc) 314 { 315 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 316 PetscErrorCode ierr; 317 PC_FieldSplitLink ilink = jac->head; 318 PetscBool fieldsplit_default = PETSC_FALSE,stokes = PETSC_FALSE,coupling = PETSC_FALSE; 319 PetscInt i; 320 321 PetscFunctionBegin; 322 /* 323 Kinda messy, but at least this now uses DMCreateFieldDecomposition() even with jac->reset. 324 Should probably be rewritten. 325 */ 326 if (!ilink || jac->reset) { 327 ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);CHKERRQ(ierr); 328 ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_coupling",&coupling,NULL);CHKERRQ(ierr); 329 if (pc->dm && jac->dm_splits && !stokes && !coupling) { 330 PetscInt numFields, f, i, j; 331 char **fieldNames; 332 IS *fields; 333 DM *dms; 334 DM subdm[128]; 335 PetscBool flg; 336 337 ierr = DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms);CHKERRQ(ierr); 338 /* Allow the user to prescribe the splits */ 339 for (i = 0, flg = PETSC_TRUE;; i++) { 340 PetscInt ifields[128]; 341 IS compField; 342 char optionname[128], splitname[8]; 343 PetscInt nfields = numFields; 344 345 ierr = PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%D_fields", i);CHKERRQ(ierr); 346 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg);CHKERRQ(ierr); 347 if (!flg) break; 348 if (numFields > 128) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot currently support %d > 128 fields", numFields); 349 ierr = DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i]);CHKERRQ(ierr); 350 if (nfields == 1) { 351 ierr = PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField);CHKERRQ(ierr); 352 /* ierr = PetscPrintf(PetscObjectComm((PetscObject)pc), "%s Field Indices:", fieldNames[ifields[0]]);CHKERRQ(ierr); 353 ierr = ISView(compField, NULL);CHKERRQ(ierr); */ 354 } else { 355 ierr = PetscSNPrintf(splitname, sizeof(splitname), "%D", i);CHKERRQ(ierr); 356 ierr = PCFieldSplitSetIS(pc, splitname, compField);CHKERRQ(ierr); 357 /* ierr = PetscPrintf(PetscObjectComm((PetscObject)pc), "%s Field Indices:", splitname);CHKERRQ(ierr); 358 ierr = ISView(compField, NULL);CHKERRQ(ierr); */ 359 } 360 ierr = ISDestroy(&compField);CHKERRQ(ierr); 361 for (j = 0; j < nfields; ++j) { 362 f = ifields[j]; 363 ierr = PetscFree(fieldNames[f]);CHKERRQ(ierr); 364 ierr = ISDestroy(&fields[f]);CHKERRQ(ierr); 365 } 366 } 367 if (i == 0) { 368 for (f = 0; f < numFields; ++f) { 369 ierr = PCFieldSplitSetIS(pc, fieldNames[f], fields[f]);CHKERRQ(ierr); 370 ierr = PetscFree(fieldNames[f]);CHKERRQ(ierr); 371 ierr = ISDestroy(&fields[f]);CHKERRQ(ierr); 372 } 373 } else { 374 for (j=0; j<numFields; j++) { 375 ierr = DMDestroy(dms+j);CHKERRQ(ierr); 376 } 377 ierr = PetscFree(dms);CHKERRQ(ierr); 378 ierr = PetscMalloc1(i, &dms);CHKERRQ(ierr); 379 for (j = 0; j < i; ++j) dms[j] = subdm[j]; 380 } 381 ierr = PetscFree(fieldNames);CHKERRQ(ierr); 382 ierr = PetscFree(fields);CHKERRQ(ierr); 383 if (dms) { 384 ierr = PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n");CHKERRQ(ierr); 385 for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) { 386 const char *prefix; 387 ierr = PetscObjectGetOptionsPrefix((PetscObject)(ilink->ksp),&prefix);CHKERRQ(ierr); 388 ierr = PetscObjectSetOptionsPrefix((PetscObject)(dms[i]), prefix);CHKERRQ(ierr); 389 ierr = KSPSetDM(ilink->ksp, dms[i]);CHKERRQ(ierr); 390 ierr = KSPSetDMActive(ilink->ksp, PETSC_FALSE);CHKERRQ(ierr); 391 ierr = PetscObjectIncrementTabLevel((PetscObject)dms[i],(PetscObject)ilink->ksp,0);CHKERRQ(ierr); 392 ierr = DMDestroy(&dms[i]);CHKERRQ(ierr); 393 } 394 ierr = PetscFree(dms);CHKERRQ(ierr); 395 } 396 } else { 397 if (jac->bs <= 0) { 398 if (pc->pmat) { 399 ierr = MatGetBlockSize(pc->pmat,&jac->bs);CHKERRQ(ierr); 400 } else jac->bs = 1; 401 } 402 403 if (stokes) { 404 IS zerodiags,rest; 405 PetscInt nmin,nmax; 406 407 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 408 ierr = MatFindZeroDiagonals(pc->mat,&zerodiags);CHKERRQ(ierr); 409 ierr = ISComplement(zerodiags,nmin,nmax,&rest);CHKERRQ(ierr); 410 if (jac->reset) { 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 jac->head->is = rest; 429 jac->head->next->is = coupling; 430 } else { 431 ierr = PCFieldSplitSetIS(pc,"0",rest);CHKERRQ(ierr); 432 ierr = PCFieldSplitSetIS(pc,"1",coupling);CHKERRQ(ierr); 433 } 434 ierr = ISDestroy(&coupling);CHKERRQ(ierr); 435 ierr = ISDestroy(&rest);CHKERRQ(ierr); 436 } else { 437 if (jac->reset) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cases not yet handled when PCReset() was used"); 438 ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&fieldsplit_default,NULL);CHKERRQ(ierr); 439 if (!fieldsplit_default) { 440 /* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit() 441 then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */ 442 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 443 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 444 } 445 if (fieldsplit_default || !jac->splitdefined) { 446 ierr = PetscInfo(pc,"Using default splitting of fields\n");CHKERRQ(ierr); 447 for (i=0; i<jac->bs; i++) { 448 char splitname[8]; 449 ierr = PetscSNPrintf(splitname,sizeof(splitname),"%D",i);CHKERRQ(ierr); 450 ierr = PCFieldSplitSetFields(pc,splitname,1,&i,&i);CHKERRQ(ierr); 451 } 452 jac->defaultsplit = PETSC_TRUE; 453 } 454 } 455 } 456 } else if (jac->nsplits == 1) { 457 if (ilink->is) { 458 IS is2; 459 PetscInt nmin,nmax; 460 461 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 462 ierr = ISComplement(ilink->is,nmin,nmax,&is2);CHKERRQ(ierr); 463 ierr = PCFieldSplitSetIS(pc,"1",is2);CHKERRQ(ierr); 464 ierr = ISDestroy(&is2);CHKERRQ(ierr); 465 } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Must provide at least two sets of fields to PCFieldSplit()"); 466 } 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 KSPSolve(ilink->ksp,ilink->x,ilink->y) || \ 836 VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \ 837 VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE)) 838 839 #undef __FUNCT__ 840 #define __FUNCT__ "PCApply_FieldSplit_Schur" 841 static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y) 842 { 843 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 844 PetscErrorCode ierr; 845 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 846 KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper; 847 848 PetscFunctionBegin; 849 switch (jac->schurfactorization) { 850 case PC_FIELDSPLIT_SCHUR_FACT_DIAG: 851 /* [A00 0; 0 -S], positive definite, suitable for MINRES */ 852 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 853 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 854 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 855 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 856 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 857 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 858 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 859 ierr = VecScale(ilinkD->y,-1.);CHKERRQ(ierr); 860 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 861 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 862 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 863 break; 864 case PC_FIELDSPLIT_SCHUR_FACT_LOWER: 865 /* [A00 0; A10 S], suitable for left preconditioning */ 866 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 867 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 868 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 869 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 870 ierr = VecScale(ilinkD->x,-1.);CHKERRQ(ierr); 871 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 872 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 873 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 874 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 875 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 876 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 877 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 878 break; 879 case PC_FIELDSPLIT_SCHUR_FACT_UPPER: 880 /* [A00 A01; 0 S], suitable for right preconditioning */ 881 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 882 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 883 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 884 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 885 ierr = VecScale(ilinkA->x,-1.);CHKERRQ(ierr); 886 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 887 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 888 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 889 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 890 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 891 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 892 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 893 break; 894 case PC_FIELDSPLIT_SCHUR_FACT_FULL: 895 /* [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 */ 896 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 897 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 898 ierr = KSPSolve(kspLower,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 899 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 900 ierr = VecScale(ilinkD->x,-1.0);CHKERRQ(ierr); 901 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 902 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 903 904 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 905 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 906 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 907 908 if (kspUpper == kspA) { 909 ierr = MatMult(jac->B,ilinkD->y,ilinkA->y);CHKERRQ(ierr); 910 ierr = VecAXPY(ilinkA->x,-1.0,ilinkA->y);CHKERRQ(ierr); 911 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 912 } else { 913 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 914 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 915 ierr = KSPSolve(kspUpper,ilinkA->x,ilinkA->z);CHKERRQ(ierr); 916 ierr = VecAXPY(ilinkA->y,-1.0,ilinkA->z);CHKERRQ(ierr); 917 } 918 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 919 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 920 } 921 PetscFunctionReturn(0); 922 } 923 924 #undef __FUNCT__ 925 #define __FUNCT__ "PCApply_FieldSplit" 926 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y) 927 { 928 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 929 PetscErrorCode ierr; 930 PC_FieldSplitLink ilink = jac->head; 931 PetscInt cnt,bs; 932 KSPConvergedReason reason; 933 934 PetscFunctionBegin; 935 if (jac->type == PC_COMPOSITE_ADDITIVE) { 936 if (jac->defaultsplit) { 937 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 938 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); 939 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 940 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); 941 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 942 while (ilink) { 943 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 944 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 945 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 946 pc->failedreason = PC_SUBPC_ERROR; 947 } 948 ilink = ilink->next; 949 } 950 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 951 } else { 952 ierr = VecSet(y,0.0);CHKERRQ(ierr); 953 while (ilink) { 954 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 955 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 956 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 957 pc->failedreason = PC_SUBPC_ERROR; 958 } 959 ilink = ilink->next; 960 } 961 } 962 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) { 963 ierr = VecSet(y,0.0);CHKERRQ(ierr); 964 /* solve on first block for first block variables */ 965 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 966 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 967 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 968 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 969 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 970 pc->failedreason = PC_SUBPC_ERROR; 971 } 972 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 973 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 974 975 /* compute the residual only onto second block variables using first block variables */ 976 ierr = MatMult(jac->Afield[1],ilink->y,ilink->next->x);CHKERRQ(ierr); 977 ilink = ilink->next; 978 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 979 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 980 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 981 982 /* solve on second block variables */ 983 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 984 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 985 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 986 pc->failedreason = PC_SUBPC_ERROR; 987 } 988 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 989 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 990 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 991 if (!jac->w1) { 992 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 993 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 994 } 995 ierr = VecSet(y,0.0);CHKERRQ(ierr); 996 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 997 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 998 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 999 pc->failedreason = PC_SUBPC_ERROR; 1000 } 1001 cnt = 1; 1002 while (ilink->next) { 1003 ilink = ilink->next; 1004 /* compute the residual only over the part of the vector needed */ 1005 ierr = MatMult(jac->Afield[cnt++],y,ilink->x);CHKERRQ(ierr); 1006 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 1007 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1008 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1009 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1010 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1011 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1012 pc->failedreason = PC_SUBPC_ERROR; 1013 } 1014 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1015 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1016 } 1017 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1018 cnt -= 2; 1019 while (ilink->previous) { 1020 ilink = ilink->previous; 1021 /* compute the residual only over the part of the vector needed */ 1022 ierr = MatMult(jac->Afield[cnt--],y,ilink->x);CHKERRQ(ierr); 1023 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 1024 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1025 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1026 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1027 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1028 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1029 pc->failedreason = PC_SUBPC_ERROR; 1030 } 1031 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1032 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1033 } 1034 } 1035 } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type); 1036 PetscFunctionReturn(0); 1037 } 1038 1039 #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \ 1040 (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 1041 VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 1042 KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \ 1043 VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \ 1044 VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE)) 1045 1046 #undef __FUNCT__ 1047 #define __FUNCT__ "PCApplyTranspose_FieldSplit" 1048 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y) 1049 { 1050 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1051 PetscErrorCode ierr; 1052 PC_FieldSplitLink ilink = jac->head; 1053 PetscInt bs; 1054 1055 PetscFunctionBegin; 1056 if (jac->type == PC_COMPOSITE_ADDITIVE) { 1057 if (jac->defaultsplit) { 1058 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 1059 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); 1060 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 1061 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); 1062 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 1063 while (ilink) { 1064 ierr = KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1065 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1066 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1067 pc->failedreason = PC_SUBPC_ERROR; 1068 } 1069 ilink = ilink->next; 1070 } 1071 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 1072 } else { 1073 ierr = VecSet(y,0.0);CHKERRQ(ierr); 1074 while (ilink) { 1075 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 1076 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1077 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1078 pc->failedreason = PC_SUBPC_ERROR; 1079 } 1080 ilink = ilink->next; 1081 } 1082 } 1083 } else { 1084 if (!jac->w1) { 1085 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 1086 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 1087 } 1088 ierr = VecSet(y,0.0);CHKERRQ(ierr); 1089 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1090 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 1091 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1092 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1093 pc->failedreason = PC_SUBPC_ERROR; 1094 } 1095 while (ilink->next) { 1096 ilink = ilink->next; 1097 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 1098 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 1099 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 1100 } 1101 while (ilink->previous) { 1102 ilink = ilink->previous; 1103 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 1104 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 1105 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 1106 } 1107 } else { 1108 while (ilink->next) { /* get to last entry in linked list */ 1109 ilink = ilink->next; 1110 } 1111 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 1112 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1113 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1114 pc->failedreason = PC_SUBPC_ERROR; 1115 } 1116 while (ilink->previous) { 1117 ilink = ilink->previous; 1118 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 1119 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 1120 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 1121 } 1122 } 1123 } 1124 PetscFunctionReturn(0); 1125 } 1126 1127 #undef __FUNCT__ 1128 #define __FUNCT__ "PCReset_FieldSplit" 1129 static PetscErrorCode PCReset_FieldSplit(PC pc) 1130 { 1131 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1132 PetscErrorCode ierr; 1133 PC_FieldSplitLink ilink = jac->head,next; 1134 1135 PetscFunctionBegin; 1136 while (ilink) { 1137 ierr = KSPReset(ilink->ksp);CHKERRQ(ierr); 1138 ierr = VecDestroy(&ilink->x);CHKERRQ(ierr); 1139 ierr = VecDestroy(&ilink->y);CHKERRQ(ierr); 1140 ierr = VecDestroy(&ilink->z);CHKERRQ(ierr); 1141 ierr = VecScatterDestroy(&ilink->sctx);CHKERRQ(ierr); 1142 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1143 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1144 next = ilink->next; 1145 ilink = next; 1146 } 1147 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 1148 if (jac->mat && jac->mat != jac->pmat) { 1149 ierr = MatDestroyMatrices(jac->nsplits,&jac->mat);CHKERRQ(ierr); 1150 } else if (jac->mat) { 1151 jac->mat = NULL; 1152 } 1153 if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);} 1154 if (jac->Afield) {ierr = MatDestroyMatrices(jac->nsplits,&jac->Afield);CHKERRQ(ierr);} 1155 ierr = VecDestroy(&jac->w1);CHKERRQ(ierr); 1156 ierr = VecDestroy(&jac->w2);CHKERRQ(ierr); 1157 ierr = MatDestroy(&jac->schur);CHKERRQ(ierr); 1158 ierr = MatDestroy(&jac->schurp);CHKERRQ(ierr); 1159 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1160 ierr = KSPDestroy(&jac->kspschur);CHKERRQ(ierr); 1161 ierr = KSPDestroy(&jac->kspupper);CHKERRQ(ierr); 1162 ierr = MatDestroy(&jac->B);CHKERRQ(ierr); 1163 ierr = MatDestroy(&jac->C);CHKERRQ(ierr); 1164 jac->reset = PETSC_TRUE; 1165 PetscFunctionReturn(0); 1166 } 1167 1168 #undef __FUNCT__ 1169 #define __FUNCT__ "PCDestroy_FieldSplit" 1170 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 1171 { 1172 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1173 PetscErrorCode ierr; 1174 PC_FieldSplitLink ilink = jac->head,next; 1175 1176 PetscFunctionBegin; 1177 ierr = PCReset_FieldSplit(pc);CHKERRQ(ierr); 1178 while (ilink) { 1179 ierr = KSPDestroy(&ilink->ksp);CHKERRQ(ierr); 1180 next = ilink->next; 1181 ierr = PetscFree(ilink->splitname);CHKERRQ(ierr); 1182 ierr = PetscFree(ilink->fields);CHKERRQ(ierr); 1183 ierr = PetscFree(ilink->fields_col);CHKERRQ(ierr); 1184 ierr = PetscFree(ilink);CHKERRQ(ierr); 1185 ilink = next; 1186 } 1187 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 1188 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1189 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);CHKERRQ(ierr); 1190 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);CHKERRQ(ierr); 1191 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);CHKERRQ(ierr); 1192 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);CHKERRQ(ierr); 1193 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);CHKERRQ(ierr); 1194 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL);CHKERRQ(ierr); 1195 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL);CHKERRQ(ierr); 1196 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);CHKERRQ(ierr); 1197 PetscFunctionReturn(0); 1198 } 1199 1200 #undef __FUNCT__ 1201 #define __FUNCT__ "PCSetFromOptions_FieldSplit" 1202 static PetscErrorCode PCSetFromOptions_FieldSplit(PetscOptionItems *PetscOptionsObject,PC pc) 1203 { 1204 PetscErrorCode ierr; 1205 PetscInt bs; 1206 PetscBool flg,stokes = PETSC_FALSE; 1207 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1208 PCCompositeType ctype; 1209 1210 PetscFunctionBegin; 1211 ierr = PetscOptionsHead(PetscOptionsObject,"FieldSplit options");CHKERRQ(ierr); 1212 ierr = PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);CHKERRQ(ierr); 1213 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 1214 if (flg) { 1215 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 1216 } 1217 jac->diag_use_amat = pc->useAmat; 1218 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); 1219 jac->offdiag_use_amat = pc->useAmat; 1220 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); 1221 /* FIXME: No programmatic equivalent to the following. */ 1222 ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);CHKERRQ(ierr); 1223 if (stokes) { 1224 ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);CHKERRQ(ierr); 1225 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF; 1226 } 1227 1228 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr); 1229 if (flg) { 1230 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr); 1231 } 1232 /* Only setup fields once */ 1233 if ((jac->bs > 0) && (jac->nsplits == 0)) { 1234 /* only allow user to set fields from command line if bs is already known. 1235 otherwise user can set them in PCFieldSplitSetDefaults() */ 1236 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 1237 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 1238 } 1239 if (jac->type == PC_COMPOSITE_SCHUR) { 1240 ierr = PetscOptionsGetEnum(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);CHKERRQ(ierr); 1241 if (flg) {ierr = PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");CHKERRQ(ierr);} 1242 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); 1243 ierr = PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);CHKERRQ(ierr); 1244 } 1245 ierr = PetscOptionsTail();CHKERRQ(ierr); 1246 PetscFunctionReturn(0); 1247 } 1248 1249 /*------------------------------------------------------------------------------------*/ 1250 1251 #undef __FUNCT__ 1252 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 1253 static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1254 { 1255 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1256 PetscErrorCode ierr; 1257 PC_FieldSplitLink ilink,next = jac->head; 1258 char prefix[128]; 1259 PetscInt i; 1260 1261 PetscFunctionBegin; 1262 if (jac->splitdefined) { 1263 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1264 PetscFunctionReturn(0); 1265 } 1266 for (i=0; i<n; i++) { 1267 if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 1268 if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 1269 } 1270 ierr = PetscNew(&ilink);CHKERRQ(ierr); 1271 if (splitname) { 1272 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1273 } else { 1274 ierr = PetscMalloc1(3,&ilink->splitname);CHKERRQ(ierr); 1275 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr); 1276 } 1277 ierr = PetscMalloc1(n,&ilink->fields);CHKERRQ(ierr); 1278 ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 1279 ierr = PetscMalloc1(n,&ilink->fields_col);CHKERRQ(ierr); 1280 ierr = PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));CHKERRQ(ierr); 1281 1282 ilink->nfields = n; 1283 ilink->next = NULL; 1284 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr); 1285 ierr = KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);CHKERRQ(ierr); 1286 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1287 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1288 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1289 1290 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr); 1291 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1292 1293 if (!next) { 1294 jac->head = ilink; 1295 ilink->previous = NULL; 1296 } else { 1297 while (next->next) { 1298 next = next->next; 1299 } 1300 next->next = ilink; 1301 ilink->previous = next; 1302 } 1303 jac->nsplits++; 1304 PetscFunctionReturn(0); 1305 } 1306 1307 #undef __FUNCT__ 1308 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit_Schur" 1309 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 1310 { 1311 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1312 PetscErrorCode ierr; 1313 1314 PetscFunctionBegin; 1315 ierr = PetscMalloc1(jac->nsplits,subksp);CHKERRQ(ierr); 1316 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr); 1317 1318 (*subksp)[1] = jac->kspschur; 1319 if (n) *n = jac->nsplits; 1320 PetscFunctionReturn(0); 1321 } 1322 1323 #undef __FUNCT__ 1324 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 1325 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1326 { 1327 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1328 PetscErrorCode ierr; 1329 PetscInt cnt = 0; 1330 PC_FieldSplitLink ilink = jac->head; 1331 1332 PetscFunctionBegin; 1333 ierr = PetscMalloc1(jac->nsplits,subksp);CHKERRQ(ierr); 1334 while (ilink) { 1335 (*subksp)[cnt++] = ilink->ksp; 1336 ilink = ilink->next; 1337 } 1338 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); 1339 if (n) *n = jac->nsplits; 1340 PetscFunctionReturn(0); 1341 } 1342 1343 #undef __FUNCT__ 1344 #define __FUNCT__ "PCFieldSplitSetIS_FieldSplit" 1345 static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 1346 { 1347 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1348 PetscErrorCode ierr; 1349 PC_FieldSplitLink ilink, next = jac->head; 1350 char prefix[128]; 1351 1352 PetscFunctionBegin; 1353 if (jac->splitdefined) { 1354 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1355 PetscFunctionReturn(0); 1356 } 1357 ierr = PetscNew(&ilink);CHKERRQ(ierr); 1358 if (splitname) { 1359 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1360 } else { 1361 ierr = PetscMalloc1(8,&ilink->splitname);CHKERRQ(ierr); 1362 ierr = PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);CHKERRQ(ierr); 1363 } 1364 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1365 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1366 ilink->is = is; 1367 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1368 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1369 ilink->is_col = is; 1370 ilink->next = NULL; 1371 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr); 1372 ierr = KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);CHKERRQ(ierr); 1373 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1374 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1375 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1376 1377 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr); 1378 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1379 1380 if (!next) { 1381 jac->head = ilink; 1382 ilink->previous = NULL; 1383 } else { 1384 while (next->next) { 1385 next = next->next; 1386 } 1387 next->next = ilink; 1388 ilink->previous = next; 1389 } 1390 jac->nsplits++; 1391 PetscFunctionReturn(0); 1392 } 1393 1394 #undef __FUNCT__ 1395 #define __FUNCT__ "PCFieldSplitSetFields" 1396 /*@ 1397 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 1398 1399 Logically Collective on PC 1400 1401 Input Parameters: 1402 + pc - the preconditioner context 1403 . splitname - name of this split, if NULL the number of the split is used 1404 . n - the number of fields in this split 1405 - fields - the fields in this split 1406 1407 Level: intermediate 1408 1409 Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 1410 1411 The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block 1412 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 1413 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 1414 where the numbered entries indicate what is in the field. 1415 1416 This function is called once per split (it creates a new split each time). Solve options 1417 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1418 1419 Developer Note: This routine does not actually create the IS representing the split, that is delayed 1420 until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be 1421 available when this routine is called. 1422 1423 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 1424 1425 @*/ 1426 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1427 { 1428 PetscErrorCode ierr; 1429 1430 PetscFunctionBegin; 1431 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1432 PetscValidCharPointer(splitname,2); 1433 if (n < 1) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 1434 PetscValidIntPointer(fields,3); 1435 ierr = PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt*,const PetscInt*),(pc,splitname,n,fields,fields_col));CHKERRQ(ierr); 1436 PetscFunctionReturn(0); 1437 } 1438 1439 #undef __FUNCT__ 1440 #define __FUNCT__ "PCFieldSplitSetDiagUseAmat" 1441 /*@ 1442 PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) 1443 1444 Logically Collective on PC 1445 1446 Input Parameters: 1447 + pc - the preconditioner object 1448 - flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from 1449 1450 Options Database: 1451 . -pc_fieldsplit_diag_use_amat 1452 1453 Level: intermediate 1454 1455 .seealso: PCFieldSplitGetDiagUseAmat(), PCFieldSplitSetOffDiagUseAmat(), PCFIELDSPLIT 1456 1457 @*/ 1458 PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg) 1459 { 1460 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1461 PetscBool isfs; 1462 PetscErrorCode ierr; 1463 1464 PetscFunctionBegin; 1465 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1466 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1467 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1468 jac->diag_use_amat = flg; 1469 PetscFunctionReturn(0); 1470 } 1471 1472 #undef __FUNCT__ 1473 #define __FUNCT__ "PCFieldSplitGetDiagUseAmat" 1474 /*@ 1475 PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) 1476 1477 Logically Collective on PC 1478 1479 Input Parameters: 1480 . pc - the preconditioner object 1481 1482 Output Parameters: 1483 . flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from 1484 1485 1486 Level: intermediate 1487 1488 .seealso: PCFieldSplitSetDiagUseAmat(), PCFieldSplitGetOffDiagUseAmat(), PCFIELDSPLIT 1489 1490 @*/ 1491 PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg) 1492 { 1493 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1494 PetscBool isfs; 1495 PetscErrorCode ierr; 1496 1497 PetscFunctionBegin; 1498 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1499 PetscValidPointer(flg,2); 1500 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1501 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1502 *flg = jac->diag_use_amat; 1503 PetscFunctionReturn(0); 1504 } 1505 1506 #undef __FUNCT__ 1507 #define __FUNCT__ "PCFieldSplitSetOffDiagUseAmat" 1508 /*@ 1509 PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) 1510 1511 Logically Collective on PC 1512 1513 Input Parameters: 1514 + pc - the preconditioner object 1515 - flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 1516 1517 Options Database: 1518 . -pc_fieldsplit_off_diag_use_amat 1519 1520 Level: intermediate 1521 1522 .seealso: PCFieldSplitGetOffDiagUseAmat(), PCFieldSplitSetDiagUseAmat(), PCFIELDSPLIT 1523 1524 @*/ 1525 PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg) 1526 { 1527 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1528 PetscBool isfs; 1529 PetscErrorCode ierr; 1530 1531 PetscFunctionBegin; 1532 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1533 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1534 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1535 jac->offdiag_use_amat = flg; 1536 PetscFunctionReturn(0); 1537 } 1538 1539 #undef __FUNCT__ 1540 #define __FUNCT__ "PCFieldSplitGetOffDiagUseAmat" 1541 /*@ 1542 PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) 1543 1544 Logically Collective on PC 1545 1546 Input Parameters: 1547 . pc - the preconditioner object 1548 1549 Output Parameters: 1550 . flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 1551 1552 1553 Level: intermediate 1554 1555 .seealso: PCFieldSplitSetOffDiagUseAmat(), PCFieldSplitGetDiagUseAmat(), PCFIELDSPLIT 1556 1557 @*/ 1558 PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg) 1559 { 1560 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1561 PetscBool isfs; 1562 PetscErrorCode ierr; 1563 1564 PetscFunctionBegin; 1565 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1566 PetscValidPointer(flg,2); 1567 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1568 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1569 *flg = jac->offdiag_use_amat; 1570 PetscFunctionReturn(0); 1571 } 1572 1573 1574 1575 #undef __FUNCT__ 1576 #define __FUNCT__ "PCFieldSplitSetIS" 1577 /*@C 1578 PCFieldSplitSetIS - Sets the exact elements for field 1579 1580 Logically Collective on PC 1581 1582 Input Parameters: 1583 + pc - the preconditioner context 1584 . splitname - name of this split, if NULL the number of the split is used 1585 - is - the index set that defines the vector elements in this field 1586 1587 1588 Notes: 1589 Use PCFieldSplitSetFields(), for fields defined by strided types. 1590 1591 This function is called once per split (it creates a new split each time). Solve options 1592 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1593 1594 Level: intermediate 1595 1596 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 1597 1598 @*/ 1599 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 1600 { 1601 PetscErrorCode ierr; 1602 1603 PetscFunctionBegin; 1604 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1605 if (splitname) PetscValidCharPointer(splitname,2); 1606 PetscValidHeaderSpecific(is,IS_CLASSID,3); 1607 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr); 1608 PetscFunctionReturn(0); 1609 } 1610 1611 #undef __FUNCT__ 1612 #define __FUNCT__ "PCFieldSplitGetIS" 1613 /*@ 1614 PCFieldSplitGetIS - Retrieves the elements for a field as an IS 1615 1616 Logically Collective on PC 1617 1618 Input Parameters: 1619 + pc - the preconditioner context 1620 - splitname - name of this split 1621 1622 Output Parameter: 1623 - is - the index set that defines the vector elements in this field, or NULL if the field is not found 1624 1625 Level: intermediate 1626 1627 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS() 1628 1629 @*/ 1630 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is) 1631 { 1632 PetscErrorCode ierr; 1633 1634 PetscFunctionBegin; 1635 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1636 PetscValidCharPointer(splitname,2); 1637 PetscValidPointer(is,3); 1638 { 1639 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 1640 PC_FieldSplitLink ilink = jac->head; 1641 PetscBool found; 1642 1643 *is = NULL; 1644 while (ilink) { 1645 ierr = PetscStrcmp(ilink->splitname, splitname, &found);CHKERRQ(ierr); 1646 if (found) { 1647 *is = ilink->is; 1648 break; 1649 } 1650 ilink = ilink->next; 1651 } 1652 } 1653 PetscFunctionReturn(0); 1654 } 1655 1656 #undef __FUNCT__ 1657 #define __FUNCT__ "PCFieldSplitSetBlockSize" 1658 /*@ 1659 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 1660 fieldsplit preconditioner. If not set the matrix block size is used. 1661 1662 Logically Collective on PC 1663 1664 Input Parameters: 1665 + pc - the preconditioner context 1666 - bs - the block size 1667 1668 Level: intermediate 1669 1670 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 1671 1672 @*/ 1673 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 1674 { 1675 PetscErrorCode ierr; 1676 1677 PetscFunctionBegin; 1678 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1679 PetscValidLogicalCollectiveInt(pc,bs,2); 1680 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr); 1681 PetscFunctionReturn(0); 1682 } 1683 1684 #undef __FUNCT__ 1685 #define __FUNCT__ "PCFieldSplitGetSubKSP" 1686 /*@C 1687 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1688 1689 Collective on KSP 1690 1691 Input Parameter: 1692 . pc - the preconditioner context 1693 1694 Output Parameters: 1695 + n - the number of splits 1696 - pc - the array of KSP contexts 1697 1698 Note: 1699 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 1700 (not the KSP just the array that contains them). 1701 1702 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 1703 1704 Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs. 1705 You can call PCFieldSplitGetSubKSP(pc,n,NULL_OBJECT,ierr) to determine how large the 1706 KSP array must be. 1707 1708 1709 Level: advanced 1710 1711 .seealso: PCFIELDSPLIT 1712 @*/ 1713 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1714 { 1715 PetscErrorCode ierr; 1716 1717 PetscFunctionBegin; 1718 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1719 if (n) PetscValidIntPointer(n,2); 1720 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 1721 PetscFunctionReturn(0); 1722 } 1723 1724 #undef __FUNCT__ 1725 #define __FUNCT__ "PCFieldSplitSetSchurPre" 1726 /*@ 1727 PCFieldSplitSetSchurPre - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the 1728 A11 matrix. Otherwise no preconditioner is used. 1729 1730 Collective on PC 1731 1732 Input Parameters: 1733 + pc - the preconditioner context 1734 . ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11 (default), PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_PRE_USER 1735 - userpre - matrix to use for preconditioning, or NULL 1736 1737 Options Database: 1738 . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> default is a11 1739 1740 Notes: 1741 $ If ptype is 1742 $ user then the preconditioner for the Schur complement is generated from the user provided matrix (pre argument 1743 $ to this function). 1744 $ a11 then the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner 1745 $ matrix associated with the Schur complement (i.e. A11), not he Schur complement matrix 1746 $ full then the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation computed internally by PFIELDSPLIT (this is expensive) 1747 $ self the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix: 1748 $ The only preconditioner that currently works with this symbolic respresentation matrix object is the PCLSC 1749 $ preconditioner 1750 $ selfp then the preconditioning for the Schur complement is generated from an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01 1751 $ This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be 1752 $ lumped before extracting the diagonal: -fieldsplit_1_mat_schur_complement_ainv_type lump 1753 1754 When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense 1755 with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and 1756 -fieldsplit_1_pc_type lsc which uses the least squares commutator to compute a preconditioner for the Schur complement. 1757 1758 Level: intermediate 1759 1760 .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, 1761 MatSchurComplementSetAinvType(), PCLSC 1762 1763 @*/ 1764 PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1765 { 1766 PetscErrorCode ierr; 1767 1768 PetscFunctionBegin; 1769 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1770 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr); 1771 PetscFunctionReturn(0); 1772 } 1773 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */ 1774 1775 #undef __FUNCT__ 1776 #define __FUNCT__ "PCFieldSplitGetSchurPre" 1777 /*@ 1778 PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be 1779 preconditioned. See PCFieldSplitSetSchurPre() for details. 1780 1781 Logically Collective on PC 1782 1783 Input Parameters: 1784 . pc - the preconditioner context 1785 1786 Output Parameters: 1787 + ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11, PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_PRE_USER 1788 - userpre - matrix to use for preconditioning (with PC_FIELDSPLIT_PRE_USER), or NULL 1789 1790 Level: intermediate 1791 1792 .seealso: PCFieldSplitSetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC 1793 1794 @*/ 1795 PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre) 1796 { 1797 PetscErrorCode ierr; 1798 1799 PetscFunctionBegin; 1800 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1801 ierr = PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));CHKERRQ(ierr); 1802 PetscFunctionReturn(0); 1803 } 1804 1805 #undef __FUNCT__ 1806 #define __FUNCT__ "PCFieldSplitSchurGetS" 1807 /*@ 1808 PCFieldSplitSchurGetS - extract the MatSchurComplement object used by this PC in case it needs to be configured separately 1809 1810 Not collective 1811 1812 Input Parameter: 1813 . pc - the preconditioner context 1814 1815 Output Parameter: 1816 . S - the Schur complement matrix 1817 1818 Notes: 1819 This matrix should not be destroyed using MatDestroy(); rather, use PCFieldSplitSchurRestoreS(). 1820 1821 Level: advanced 1822 1823 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurRestoreS() 1824 1825 @*/ 1826 PetscErrorCode PCFieldSplitSchurGetS(PC pc,Mat *S) 1827 { 1828 PetscErrorCode ierr; 1829 const char* t; 1830 PetscBool isfs; 1831 PC_FieldSplit *jac; 1832 1833 PetscFunctionBegin; 1834 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1835 ierr = PetscObjectGetType((PetscObject)pc,&t);CHKERRQ(ierr); 1836 ierr = PetscStrcmp(t,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1837 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 1838 jac = (PC_FieldSplit*)pc->data; 1839 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type); 1840 if (S) *S = jac->schur; 1841 PetscFunctionReturn(0); 1842 } 1843 1844 #undef __FUNCT__ 1845 #define __FUNCT__ "PCFieldSplitSchurRestoreS" 1846 /*@ 1847 PCFieldSplitSchurRestoreS - restores the MatSchurComplement object used by this PC 1848 1849 Not collective 1850 1851 Input Parameters: 1852 + pc - the preconditioner context 1853 . S - the Schur complement matrix 1854 1855 Level: advanced 1856 1857 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurGetS() 1858 1859 @*/ 1860 PetscErrorCode PCFieldSplitSchurRestoreS(PC pc,Mat *S) 1861 { 1862 PetscErrorCode ierr; 1863 const char* t; 1864 PetscBool isfs; 1865 PC_FieldSplit *jac; 1866 1867 PetscFunctionBegin; 1868 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1869 ierr = PetscObjectGetType((PetscObject)pc,&t);CHKERRQ(ierr); 1870 ierr = PetscStrcmp(t,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1871 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 1872 jac = (PC_FieldSplit*)pc->data; 1873 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type); 1874 if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten"); 1875 PetscFunctionReturn(0); 1876 } 1877 1878 1879 #undef __FUNCT__ 1880 #define __FUNCT__ "PCFieldSplitSetSchurPre_FieldSplit" 1881 static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1882 { 1883 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1884 PetscErrorCode ierr; 1885 1886 PetscFunctionBegin; 1887 jac->schurpre = ptype; 1888 if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) { 1889 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1890 jac->schur_user = pre; 1891 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 1892 } 1893 PetscFunctionReturn(0); 1894 } 1895 1896 #undef __FUNCT__ 1897 #define __FUNCT__ "PCFieldSplitGetSchurPre_FieldSplit" 1898 static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre) 1899 { 1900 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1901 1902 PetscFunctionBegin; 1903 *ptype = jac->schurpre; 1904 *pre = jac->schur_user; 1905 PetscFunctionReturn(0); 1906 } 1907 1908 #undef __FUNCT__ 1909 #define __FUNCT__ "PCFieldSplitSetSchurFactType" 1910 /*@ 1911 PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain 1912 1913 Collective on PC 1914 1915 Input Parameters: 1916 + pc - the preconditioner context 1917 - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default 1918 1919 Options Database: 1920 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full 1921 1922 1923 Level: intermediate 1924 1925 Notes: 1926 The FULL factorization is 1927 1928 $ (A B) = (1 0) (A 0) (1 Ainv*B) 1929 $ (C D) (C*Ainv 1) (0 S) (0 1 ) 1930 1931 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, 1932 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). 1933 1934 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 1935 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 1936 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, 1937 the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations. 1938 1939 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 1940 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). 1941 1942 References: 1943 Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000) pp. 1969-1972. 1944 1945 Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001), pp. 1050-1051. 1946 1947 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType 1948 @*/ 1949 PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype) 1950 { 1951 PetscErrorCode ierr; 1952 1953 PetscFunctionBegin; 1954 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1955 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));CHKERRQ(ierr); 1956 PetscFunctionReturn(0); 1957 } 1958 1959 #undef __FUNCT__ 1960 #define __FUNCT__ "PCFieldSplitSetSchurFactType_FieldSplit" 1961 static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype) 1962 { 1963 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1964 1965 PetscFunctionBegin; 1966 jac->schurfactorization = ftype; 1967 PetscFunctionReturn(0); 1968 } 1969 1970 #undef __FUNCT__ 1971 #define __FUNCT__ "PCFieldSplitGetSchurBlocks" 1972 /*@C 1973 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 1974 1975 Collective on KSP 1976 1977 Input Parameter: 1978 . pc - the preconditioner context 1979 1980 Output Parameters: 1981 + A00 - the (0,0) block 1982 . A01 - the (0,1) block 1983 . A10 - the (1,0) block 1984 - A11 - the (1,1) block 1985 1986 Level: advanced 1987 1988 .seealso: PCFIELDSPLIT 1989 @*/ 1990 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 1991 { 1992 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 1993 1994 PetscFunctionBegin; 1995 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1996 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 1997 if (A00) *A00 = jac->pmat[0]; 1998 if (A01) *A01 = jac->B; 1999 if (A10) *A10 = jac->C; 2000 if (A11) *A11 = jac->pmat[1]; 2001 PetscFunctionReturn(0); 2002 } 2003 2004 #undef __FUNCT__ 2005 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 2006 static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 2007 { 2008 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2009 PetscErrorCode ierr; 2010 2011 PetscFunctionBegin; 2012 jac->type = type; 2013 if (type == PC_COMPOSITE_SCHUR) { 2014 pc->ops->apply = PCApply_FieldSplit_Schur; 2015 pc->ops->view = PCView_FieldSplit_Schur; 2016 2017 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 2018 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);CHKERRQ(ierr); 2019 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);CHKERRQ(ierr); 2020 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);CHKERRQ(ierr); 2021 2022 } else { 2023 pc->ops->apply = PCApply_FieldSplit; 2024 pc->ops->view = PCView_FieldSplit; 2025 2026 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 2027 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);CHKERRQ(ierr); 2028 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);CHKERRQ(ierr); 2029 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);CHKERRQ(ierr); 2030 } 2031 PetscFunctionReturn(0); 2032 } 2033 2034 #undef __FUNCT__ 2035 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 2036 static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 2037 { 2038 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2039 2040 PetscFunctionBegin; 2041 if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 2042 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); 2043 jac->bs = bs; 2044 PetscFunctionReturn(0); 2045 } 2046 2047 #undef __FUNCT__ 2048 #define __FUNCT__ "PCFieldSplitSetType" 2049 /*@ 2050 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 2051 2052 Collective on PC 2053 2054 Input Parameter: 2055 . pc - the preconditioner context 2056 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 2057 2058 Options Database Key: 2059 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 2060 2061 Level: Intermediate 2062 2063 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 2064 2065 .seealso: PCCompositeSetType() 2066 2067 @*/ 2068 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 2069 { 2070 PetscErrorCode ierr; 2071 2072 PetscFunctionBegin; 2073 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2074 ierr = PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));CHKERRQ(ierr); 2075 PetscFunctionReturn(0); 2076 } 2077 2078 #undef __FUNCT__ 2079 #define __FUNCT__ "PCFieldSplitGetType" 2080 /*@ 2081 PCFieldSplitGetType - Gets the type of fieldsplit preconditioner. 2082 2083 Not collective 2084 2085 Input Parameter: 2086 . pc - the preconditioner context 2087 2088 Output Parameter: 2089 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 2090 2091 Level: Intermediate 2092 2093 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 2094 .seealso: PCCompositeSetType() 2095 @*/ 2096 PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type) 2097 { 2098 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2099 2100 PetscFunctionBegin; 2101 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2102 PetscValidIntPointer(type,2); 2103 *type = jac->type; 2104 PetscFunctionReturn(0); 2105 } 2106 2107 #undef __FUNCT__ 2108 #define __FUNCT__ "PCFieldSplitSetDMSplits" 2109 /*@ 2110 PCFieldSplitSetDMSplits - Flags whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 2111 2112 Logically Collective 2113 2114 Input Parameters: 2115 + pc - the preconditioner context 2116 - flg - boolean indicating whether to use field splits defined by the DM 2117 2118 Options Database Key: 2119 . -pc_fieldsplit_dm_splits 2120 2121 Level: Intermediate 2122 2123 .keywords: PC, DM, composite preconditioner, additive, multiplicative 2124 2125 .seealso: PCFieldSplitGetDMSplits() 2126 2127 @*/ 2128 PetscErrorCode PCFieldSplitSetDMSplits(PC pc,PetscBool flg) 2129 { 2130 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2131 PetscBool isfs; 2132 PetscErrorCode ierr; 2133 2134 PetscFunctionBegin; 2135 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2136 PetscValidLogicalCollectiveBool(pc,flg,2); 2137 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 2138 if (isfs) { 2139 jac->dm_splits = flg; 2140 } 2141 PetscFunctionReturn(0); 2142 } 2143 2144 2145 #undef __FUNCT__ 2146 #define __FUNCT__ "PCFieldSplitGetDMSplits" 2147 /*@ 2148 PCFieldSplitGetDMSplits - Returns flag indicating whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 2149 2150 Logically Collective 2151 2152 Input Parameter: 2153 . pc - the preconditioner context 2154 2155 Output Parameter: 2156 . flg - boolean indicating whether to use field splits defined by the DM 2157 2158 Level: Intermediate 2159 2160 .keywords: PC, DM, composite preconditioner, additive, multiplicative 2161 2162 .seealso: PCFieldSplitSetDMSplits() 2163 2164 @*/ 2165 PetscErrorCode PCFieldSplitGetDMSplits(PC pc,PetscBool* flg) 2166 { 2167 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2168 PetscBool isfs; 2169 PetscErrorCode ierr; 2170 2171 PetscFunctionBegin; 2172 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2173 PetscValidPointer(flg,2); 2174 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 2175 if (isfs) { 2176 if(flg) *flg = jac->dm_splits; 2177 } 2178 PetscFunctionReturn(0); 2179 } 2180 2181 /* -------------------------------------------------------------------------------------*/ 2182 /*MC 2183 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 2184 fields or groups of fields. See the users manual section "Solving Block Matrices" for more details. 2185 2186 To set options on the solvers for each block append -fieldsplit_ to all the PC 2187 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 2188 2189 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 2190 and set the options directly on the resulting KSP object 2191 2192 Level: intermediate 2193 2194 Options Database Keys: 2195 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 2196 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 2197 been supplied explicitly by -pc_fieldsplit_%d_fields 2198 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 2199 . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur> - type of relaxation or factorization splitting 2200 . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11; see PCFieldSplitSetSchurPre() 2201 . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver 2202 2203 - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ 2204 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 2205 2206 Notes: 2207 Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 2208 to define a field by an arbitrary collection of entries. 2209 2210 If no fields are set the default is used. The fields are defined by entries strided by bs, 2211 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 2212 if this is not called the block size defaults to the blocksize of the second matrix passed 2213 to KSPSetOperators()/PCSetOperators(). 2214 2215 $ For the Schur complement preconditioner if J = ( A00 A01 ) 2216 $ ( A10 A11 ) 2217 $ the preconditioner using full factorization is 2218 $ ( I -ksp(A00) A01 ) ( inv(A00) 0 ) ( I 0 ) 2219 $ ( 0 I ) ( 0 ksp(S) ) ( -A10 ksp(A00) I ) 2220 where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. S is the Schur complement 2221 $ S = A11 - A10 ksp(A00) A01 2222 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 2223 in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0, 2224 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 2225 A11 is used to construct a preconditioner for S, use PCFieldSplitSetSchurPre() for all the possible ways to construct the preconditioner for S. 2226 2227 The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above, 2228 diag gives 2229 $ ( inv(A00) 0 ) 2230 $ ( 0 -ksp(S) ) 2231 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 2232 $ ( A00 0 ) 2233 $ ( A10 S ) 2234 where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of 2235 $ ( A00 A01 ) 2236 $ ( 0 S ) 2237 where again the inverses of A00 and S are applied using KSPs. 2238 2239 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 2240 is used automatically for a second block. 2241 2242 The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1. 2243 Generally it should be used with the AIJ format. 2244 2245 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 2246 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT 2247 inside a smoother resulting in "Distributive Smoothers". 2248 2249 Concepts: physics based preconditioners, block preconditioners 2250 2251 There is a nice discussion of block preconditioners in 2252 2253 [El08] A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations 2254 Howard Elman, V.E. Howle, John Shadid, Robert Shuttleworth, Ray Tuminaro, Journal of Computational Physics 227 (2008) 1790--1808 2255 http://chess.cs.umd.edu/~elman/papers/tax.pdf 2256 2257 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 2258 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 2259 of backsolving for the saturations in the last step it solves a full coupled (ILU) system for updates to all the variables. 2260 2261 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC, 2262 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre(), 2263 MatSchurComplementSetAinvType() 2264 M*/ 2265 2266 #undef __FUNCT__ 2267 #define __FUNCT__ "PCCreate_FieldSplit" 2268 PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc) 2269 { 2270 PetscErrorCode ierr; 2271 PC_FieldSplit *jac; 2272 2273 PetscFunctionBegin; 2274 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr); 2275 2276 jac->bs = -1; 2277 jac->nsplits = 0; 2278 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 2279 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 2280 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL; 2281 jac->dm_splits = PETSC_TRUE; 2282 2283 pc->data = (void*)jac; 2284 2285 pc->ops->apply = PCApply_FieldSplit; 2286 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 2287 pc->ops->setup = PCSetUp_FieldSplit; 2288 pc->ops->reset = PCReset_FieldSplit; 2289 pc->ops->destroy = PCDestroy_FieldSplit; 2290 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 2291 pc->ops->view = PCView_FieldSplit; 2292 pc->ops->applyrichardson = 0; 2293 2294 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 2295 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 2296 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 2297 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 2298 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 2299 PetscFunctionReturn(0); 2300 } 2301 2302 2303 2304