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 (the lumped) 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 = PetscDrawBoxedString(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 = PetscDrawBoxedString(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)->prefix,optionname,ifields,&nfields,&flg);CHKERRQ(ierr); 289 ierr = PetscOptionsGetIntArray(((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)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);CHKERRQ(ierr); 328 ierr = PetscOptionsGetBool(((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)->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 if (jac->reset) { 427 jac->head->is = coupling; 428 jac->head->next->is = rest; 429 } else { 430 ierr = PCFieldSplitSetIS(pc,"0",coupling);CHKERRQ(ierr); 431 ierr = PCFieldSplitSetIS(pc,"1",rest);CHKERRQ(ierr); 432 } 433 ierr = ISDestroy(&coupling);CHKERRQ(ierr); 434 ierr = ISDestroy(&rest);CHKERRQ(ierr); 435 } else { 436 if (jac->reset) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cases not yet handled when PCReset() was used"); 437 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&fieldsplit_default,NULL);CHKERRQ(ierr); 438 if (!fieldsplit_default) { 439 /* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit() 440 then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */ 441 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 442 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 443 } 444 if (fieldsplit_default || !jac->splitdefined) { 445 ierr = PetscInfo(pc,"Using default splitting of fields\n");CHKERRQ(ierr); 446 for (i=0; i<jac->bs; i++) { 447 char splitname[8]; 448 ierr = PetscSNPrintf(splitname,sizeof(splitname),"%D",i);CHKERRQ(ierr); 449 ierr = PCFieldSplitSetFields(pc,splitname,1,&i,&i);CHKERRQ(ierr); 450 } 451 jac->defaultsplit = PETSC_TRUE; 452 } 453 } 454 } 455 } else if (jac->nsplits == 1) { 456 if (ilink->is) { 457 IS is2; 458 PetscInt nmin,nmax; 459 460 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 461 ierr = ISComplement(ilink->is,nmin,nmax,&is2);CHKERRQ(ierr); 462 ierr = PCFieldSplitSetIS(pc,"1",is2);CHKERRQ(ierr); 463 ierr = ISDestroy(&is2);CHKERRQ(ierr); 464 } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Must provide at least two sets of fields to PCFieldSplit()"); 465 } 466 467 468 if (jac->nsplits < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unhandled case, must have at least two fields, not %d", jac->nsplits); 469 PetscFunctionReturn(0); 470 } 471 472 PETSC_EXTERN PetscErrorCode PetscOptionsFindPairPrefix_Private(const char pre[], const char name[], char *value[], PetscBool *flg); 473 474 #undef __FUNCT__ 475 #define __FUNCT__ "PCSetUp_FieldSplit" 476 static PetscErrorCode PCSetUp_FieldSplit(PC pc) 477 { 478 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 479 PetscErrorCode ierr; 480 PC_FieldSplitLink ilink; 481 PetscInt i,nsplit; 482 PetscBool sorted, sorted_col; 483 484 PetscFunctionBegin; 485 ierr = PCFieldSplitSetDefaults(pc);CHKERRQ(ierr); 486 nsplit = jac->nsplits; 487 ilink = jac->head; 488 489 /* get the matrices for each split */ 490 if (!jac->issetup) { 491 PetscInt rstart,rend,nslots,bs; 492 493 jac->issetup = PETSC_TRUE; 494 495 /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */ 496 if (jac->defaultsplit || !ilink->is) { 497 if (jac->bs <= 0) jac->bs = nsplit; 498 } 499 bs = jac->bs; 500 ierr = MatGetOwnershipRange(pc->pmat,&rstart,&rend);CHKERRQ(ierr); 501 nslots = (rend - rstart)/bs; 502 for (i=0; i<nsplit; i++) { 503 if (jac->defaultsplit) { 504 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+i,nsplit,&ilink->is);CHKERRQ(ierr); 505 ierr = ISDuplicate(ilink->is,&ilink->is_col);CHKERRQ(ierr); 506 } else if (!ilink->is) { 507 if (ilink->nfields > 1) { 508 PetscInt *ii,*jj,j,k,nfields = ilink->nfields,*fields = ilink->fields,*fields_col = ilink->fields_col; 509 ierr = PetscMalloc1(ilink->nfields*nslots,&ii);CHKERRQ(ierr); 510 ierr = PetscMalloc1(ilink->nfields*nslots,&jj);CHKERRQ(ierr); 511 for (j=0; j<nslots; j++) { 512 for (k=0; k<nfields; k++) { 513 ii[nfields*j + k] = rstart + bs*j + fields[k]; 514 jj[nfields*j + k] = rstart + bs*j + fields_col[k]; 515 } 516 } 517 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is);CHKERRQ(ierr); 518 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,jj,PETSC_OWN_POINTER,&ilink->is_col);CHKERRQ(ierr); 519 } else { 520 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields[0],bs,&ilink->is);CHKERRQ(ierr); 521 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields_col[0],bs,&ilink->is_col);CHKERRQ(ierr); 522 } 523 } 524 ierr = ISSorted(ilink->is,&sorted);CHKERRQ(ierr); 525 if (ilink->is_col) { ierr = ISSorted(ilink->is_col,&sorted_col);CHKERRQ(ierr); } 526 if (!sorted || !sorted_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Fields must be sorted when creating split"); 527 ilink = ilink->next; 528 } 529 } 530 531 ilink = jac->head; 532 if (!jac->pmat) { 533 Vec xtmp; 534 535 ierr = MatCreateVecs(pc->pmat,&xtmp,NULL);CHKERRQ(ierr); 536 ierr = PetscMalloc1(nsplit,&jac->pmat);CHKERRQ(ierr); 537 ierr = PetscMalloc2(nsplit,&jac->x,nsplit,&jac->y);CHKERRQ(ierr); 538 for (i=0; i<nsplit; i++) { 539 MatNullSpace sp; 540 541 /* Check for preconditioning matrix attached to IS */ 542 ierr = PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &jac->pmat[i]);CHKERRQ(ierr); 543 if (jac->pmat[i]) { 544 ierr = PetscObjectReference((PetscObject) jac->pmat[i]);CHKERRQ(ierr); 545 if (jac->type == PC_COMPOSITE_SCHUR) { 546 jac->schur_user = jac->pmat[i]; 547 548 ierr = PetscObjectReference((PetscObject) jac->schur_user);CHKERRQ(ierr); 549 } 550 } else { 551 const char *prefix; 552 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 553 ierr = KSPGetOptionsPrefix(ilink->ksp,&prefix);CHKERRQ(ierr); 554 ierr = MatSetOptionsPrefix(jac->pmat[i],prefix);CHKERRQ(ierr); 555 ierr = MatViewFromOptions(jac->pmat[i],NULL,"-mat_view");CHKERRQ(ierr); 556 } 557 /* create work vectors for each split */ 558 ierr = MatCreateVecs(jac->pmat[i],&jac->x[i],&jac->y[i]);CHKERRQ(ierr); 559 ilink->x = jac->x[i]; ilink->y = jac->y[i]; ilink->z = NULL; 560 /* compute scatter contexts needed by multiplicative versions and non-default splits */ 561 ierr = VecScatterCreate(xtmp,ilink->is,jac->x[i],NULL,&ilink->sctx);CHKERRQ(ierr); 562 /* Check for null space attached to IS */ 563 ierr = PetscObjectQuery((PetscObject) ilink->is, "nullspace", (PetscObject*) &sp);CHKERRQ(ierr); 564 if (sp) { 565 ierr = MatSetNullSpace(jac->pmat[i], sp);CHKERRQ(ierr); 566 } 567 ierr = PetscObjectQuery((PetscObject) ilink->is, "nearnullspace", (PetscObject*) &sp);CHKERRQ(ierr); 568 if (sp) { 569 ierr = MatSetNearNullSpace(jac->pmat[i], sp);CHKERRQ(ierr); 570 } 571 ilink = ilink->next; 572 } 573 ierr = VecDestroy(&xtmp);CHKERRQ(ierr); 574 } else { 575 for (i=0; i<nsplit; i++) { 576 Mat pmat; 577 578 /* Check for preconditioning matrix attached to IS */ 579 ierr = PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &pmat);CHKERRQ(ierr); 580 if (!pmat) { 581 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 582 } 583 ilink = ilink->next; 584 } 585 } 586 if (jac->diag_use_amat) { 587 ilink = jac->head; 588 if (!jac->mat) { 589 ierr = PetscMalloc1(nsplit,&jac->mat);CHKERRQ(ierr); 590 for (i=0; i<nsplit; i++) { 591 ierr = MatGetSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->mat[i]);CHKERRQ(ierr); 592 ilink = ilink->next; 593 } 594 } else { 595 for (i=0; i<nsplit; i++) { 596 if (jac->mat[i]) {ierr = MatGetSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->mat[i]);CHKERRQ(ierr);} 597 ilink = ilink->next; 598 } 599 } 600 } else { 601 jac->mat = jac->pmat; 602 } 603 604 if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR) { 605 /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */ 606 /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */ 607 ilink = jac->head; 608 if (!jac->Afield) { 609 ierr = PetscMalloc1(nsplit,&jac->Afield);CHKERRQ(ierr); 610 for (i=0; i<nsplit; i++) { 611 ierr = MatGetSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 612 ilink = ilink->next; 613 } 614 } else { 615 for (i=0; i<nsplit; i++) { 616 ierr = MatGetSubMatrix(pc->mat,ilink->is,NULL,MAT_REUSE_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 617 ilink = ilink->next; 618 } 619 } 620 } 621 622 if (jac->type == PC_COMPOSITE_SCHUR) { 623 IS ccis; 624 PetscInt rstart,rend; 625 char lscname[256]; 626 PetscObject LSC_L; 627 628 if (nsplit != 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields"); 629 630 /* When extracting off-diagonal submatrices, we take complements from this range */ 631 ierr = MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);CHKERRQ(ierr); 632 633 /* need to handle case when one is resetting up the preconditioner */ 634 if (jac->schur) { 635 KSP kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper; 636 637 ierr = MatSchurComplementGetKSP(jac->schur, &kspInner);CHKERRQ(ierr); 638 ilink = jac->head; 639 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 640 if (jac->offdiag_use_amat) { 641 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);CHKERRQ(ierr); 642 } else { 643 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);CHKERRQ(ierr); 644 } 645 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 646 ilink = ilink->next; 647 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 648 if (jac->offdiag_use_amat) { 649 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);CHKERRQ(ierr); 650 } else { 651 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);CHKERRQ(ierr); 652 } 653 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 654 ierr = MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);CHKERRQ(ierr); 655 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) { 656 ierr = MatDestroy(&jac->schurp);CHKERRQ(ierr); 657 ierr = MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);CHKERRQ(ierr); 658 } 659 if (kspA != kspInner) { 660 ierr = KSPSetOperators(kspA,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 661 } 662 if (kspUpper != kspA) { 663 ierr = KSPSetOperators(kspUpper,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 664 } 665 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));CHKERRQ(ierr); 666 } else { 667 const char *Dprefix; 668 char schurprefix[256], schurmatprefix[256]; 669 char schurtestoption[256]; 670 MatNullSpace sp; 671 PetscBool flg; 672 673 /* extract the A01 and A10 matrices */ 674 ilink = jac->head; 675 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 676 if (jac->offdiag_use_amat) { 677 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr); 678 } else { 679 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr); 680 } 681 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 682 ilink = ilink->next; 683 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 684 if (jac->offdiag_use_amat) { 685 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 686 } else { 687 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 688 } 689 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 690 691 /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */ 692 ierr = MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);CHKERRQ(ierr); 693 ierr = MatSetType(jac->schur,MATSCHURCOMPLEMENT);CHKERRQ(ierr); 694 ierr = MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);CHKERRQ(ierr); 695 ierr = PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 696 /* 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? */ 697 ierr = MatSetOptionsPrefix(jac->schur,schurmatprefix);CHKERRQ(ierr); 698 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr); 699 ierr = MatGetNullSpace(jac->pmat[1], &sp);CHKERRQ(ierr); 700 if (sp) { 701 ierr = MatSetNullSpace(jac->schur, sp);CHKERRQ(ierr); 702 } 703 704 ierr = PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);CHKERRQ(ierr); 705 ierr = PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);CHKERRQ(ierr); 706 if (flg) { 707 DM dmInner; 708 KSP kspInner; 709 710 ierr = MatSchurComplementGetKSP(jac->schur, &kspInner);CHKERRQ(ierr); 711 ierr = PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 712 /* Indent this deeper to emphasize the "inner" nature of this solver. */ 713 ierr = PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject) pc, 2);CHKERRQ(ierr); 714 ierr = KSPSetOptionsPrefix(kspInner, schurprefix);CHKERRQ(ierr); 715 716 /* Set DM for new solver */ 717 ierr = KSPGetDM(jac->head->ksp, &dmInner);CHKERRQ(ierr); 718 ierr = KSPSetDM(kspInner, dmInner);CHKERRQ(ierr); 719 ierr = KSPSetDMActive(kspInner, PETSC_FALSE);CHKERRQ(ierr); 720 } else { 721 /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or 722 * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact, 723 * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for 724 * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make 725 * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used 726 * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */ 727 ierr = KSPSetType(jac->head->ksp,KSPGMRES);CHKERRQ(ierr); 728 ierr = MatSchurComplementSetKSP(jac->schur,jac->head->ksp);CHKERRQ(ierr); 729 } 730 ierr = KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 731 ierr = KSPSetFromOptions(jac->head->ksp);CHKERRQ(ierr); 732 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr); 733 734 ierr = PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);CHKERRQ(ierr); 735 ierr = PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);CHKERRQ(ierr); 736 if (flg) { 737 DM dmInner; 738 739 ierr = PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 740 ierr = KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);CHKERRQ(ierr); 741 ierr = KSPSetOptionsPrefix(jac->kspupper, schurprefix);CHKERRQ(ierr); 742 ierr = KSPGetDM(jac->head->ksp, &dmInner);CHKERRQ(ierr); 743 ierr = KSPSetDM(jac->kspupper, dmInner);CHKERRQ(ierr); 744 ierr = KSPSetDMActive(jac->kspupper, PETSC_FALSE);CHKERRQ(ierr); 745 ierr = KSPSetFromOptions(jac->kspupper);CHKERRQ(ierr); 746 ierr = KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 747 ierr = VecDuplicate(jac->head->x, &jac->head->z);CHKERRQ(ierr); 748 } else { 749 jac->kspupper = jac->head->ksp; 750 ierr = PetscObjectReference((PetscObject) jac->head->ksp);CHKERRQ(ierr); 751 } 752 753 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) { 754 ierr = MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);CHKERRQ(ierr); 755 } 756 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);CHKERRQ(ierr); 757 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);CHKERRQ(ierr); 758 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);CHKERRQ(ierr); 759 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) { 760 PC pcschur; 761 ierr = KSPGetPC(jac->kspschur,&pcschur);CHKERRQ(ierr); 762 ierr = PCSetType(pcschur,PCNONE);CHKERRQ(ierr); 763 /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */ 764 } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) { 765 ierr = MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user);CHKERRQ(ierr); 766 } 767 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));CHKERRQ(ierr); 768 ierr = KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix);CHKERRQ(ierr); 769 ierr = KSPSetOptionsPrefix(jac->kspschur, Dprefix);CHKERRQ(ierr); 770 /* propogate DM */ 771 { 772 DM sdm; 773 ierr = KSPGetDM(jac->head->next->ksp, &sdm);CHKERRQ(ierr); 774 if (sdm) { 775 ierr = KSPSetDM(jac->kspschur, sdm);CHKERRQ(ierr); 776 ierr = KSPSetDMActive(jac->kspschur, PETSC_FALSE);CHKERRQ(ierr); 777 } 778 } 779 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 780 /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */ 781 ierr = KSPSetFromOptions(jac->kspschur);CHKERRQ(ierr); 782 } 783 784 /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */ 785 ierr = PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname);CHKERRQ(ierr); 786 ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 787 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 788 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);CHKERRQ(ierr);} 789 ierr = PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname);CHKERRQ(ierr); 790 ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 791 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 792 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);CHKERRQ(ierr);} 793 } else { 794 /* set up the individual splits' PCs */ 795 i = 0; 796 ilink = jac->head; 797 while (ilink) { 798 ierr = KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i]);CHKERRQ(ierr); 799 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 800 if (!jac->suboptionsset) {ierr = KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr);} 801 i++; 802 ilink = ilink->next; 803 } 804 } 805 806 jac->suboptionsset = PETSC_TRUE; 807 PetscFunctionReturn(0); 808 } 809 810 #define FieldSplitSplitSolveAdd(ilink,xx,yy) \ 811 (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 812 VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 813 KSPSolve(ilink->ksp,ilink->x,ilink->y) || \ 814 VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \ 815 VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE)) 816 817 #undef __FUNCT__ 818 #define __FUNCT__ "PCApply_FieldSplit_Schur" 819 static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y) 820 { 821 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 822 PetscErrorCode ierr; 823 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 824 KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper; 825 826 PetscFunctionBegin; 827 switch (jac->schurfactorization) { 828 case PC_FIELDSPLIT_SCHUR_FACT_DIAG: 829 /* [A00 0; 0 -S], positive definite, suitable for MINRES */ 830 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 831 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 832 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 833 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 834 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 835 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 836 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 837 ierr = VecScale(ilinkD->y,-1.);CHKERRQ(ierr); 838 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 839 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 840 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 841 break; 842 case PC_FIELDSPLIT_SCHUR_FACT_LOWER: 843 /* [A00 0; A10 S], suitable for left preconditioning */ 844 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 845 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 846 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 847 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 848 ierr = VecScale(ilinkD->x,-1.);CHKERRQ(ierr); 849 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 850 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 851 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 852 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 853 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 854 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 855 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 856 break; 857 case PC_FIELDSPLIT_SCHUR_FACT_UPPER: 858 /* [A00 A01; 0 S], suitable for right preconditioning */ 859 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 860 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 861 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 862 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 863 ierr = VecScale(ilinkA->x,-1.);CHKERRQ(ierr); 864 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 865 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 866 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 867 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 868 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 869 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 870 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 871 break; 872 case PC_FIELDSPLIT_SCHUR_FACT_FULL: 873 /* [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 */ 874 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 875 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 876 ierr = KSPSolve(kspLower,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 877 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 878 ierr = VecScale(ilinkD->x,-1.0);CHKERRQ(ierr); 879 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 880 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 881 882 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 883 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 884 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 885 886 if (kspUpper == kspA) { 887 ierr = MatMult(jac->B,ilinkD->y,ilinkA->y);CHKERRQ(ierr); 888 ierr = VecAXPY(ilinkA->x,-1.0,ilinkA->y);CHKERRQ(ierr); 889 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 890 } else { 891 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 892 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 893 ierr = KSPSolve(kspUpper,ilinkA->x,ilinkA->z);CHKERRQ(ierr); 894 ierr = VecAXPY(ilinkA->y,-1.0,ilinkA->z);CHKERRQ(ierr); 895 } 896 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 897 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 898 } 899 PetscFunctionReturn(0); 900 } 901 902 #undef __FUNCT__ 903 #define __FUNCT__ "PCApply_FieldSplit" 904 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y) 905 { 906 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 907 PetscErrorCode ierr; 908 PC_FieldSplitLink ilink = jac->head; 909 PetscInt cnt,bs; 910 911 PetscFunctionBegin; 912 if (jac->type == PC_COMPOSITE_ADDITIVE) { 913 if (jac->defaultsplit) { 914 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 915 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); 916 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 917 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); 918 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 919 while (ilink) { 920 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 921 ilink = ilink->next; 922 } 923 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 924 } else { 925 ierr = VecSet(y,0.0);CHKERRQ(ierr); 926 while (ilink) { 927 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 928 ilink = ilink->next; 929 } 930 } 931 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 932 if (!jac->w1) { 933 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 934 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 935 } 936 ierr = VecSet(y,0.0);CHKERRQ(ierr); 937 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 938 cnt = 1; 939 while (ilink->next) { 940 ilink = ilink->next; 941 /* compute the residual only over the part of the vector needed */ 942 ierr = MatMult(jac->Afield[cnt++],y,ilink->x);CHKERRQ(ierr); 943 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 944 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 945 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 946 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 947 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 948 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 949 } 950 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 951 cnt -= 2; 952 while (ilink->previous) { 953 ilink = ilink->previous; 954 /* compute the residual only over the part of the vector needed */ 955 ierr = MatMult(jac->Afield[cnt--],y,ilink->x);CHKERRQ(ierr); 956 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 957 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 958 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 959 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 960 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 961 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 962 } 963 } 964 } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type); 965 PetscFunctionReturn(0); 966 } 967 968 #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \ 969 (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 970 VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 971 KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \ 972 VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \ 973 VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE)) 974 975 #undef __FUNCT__ 976 #define __FUNCT__ "PCApplyTranspose_FieldSplit" 977 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y) 978 { 979 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 980 PetscErrorCode ierr; 981 PC_FieldSplitLink ilink = jac->head; 982 PetscInt bs; 983 984 PetscFunctionBegin; 985 if (jac->type == PC_COMPOSITE_ADDITIVE) { 986 if (jac->defaultsplit) { 987 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 988 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); 989 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 990 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); 991 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 992 while (ilink) { 993 ierr = KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 994 ilink = ilink->next; 995 } 996 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 997 } else { 998 ierr = VecSet(y,0.0);CHKERRQ(ierr); 999 while (ilink) { 1000 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 1001 ilink = ilink->next; 1002 } 1003 } 1004 } else { 1005 if (!jac->w1) { 1006 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 1007 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 1008 } 1009 ierr = VecSet(y,0.0);CHKERRQ(ierr); 1010 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1011 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 1012 while (ilink->next) { 1013 ilink = ilink->next; 1014 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 1015 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 1016 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 1017 } 1018 while (ilink->previous) { 1019 ilink = ilink->previous; 1020 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 1021 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 1022 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 1023 } 1024 } else { 1025 while (ilink->next) { /* get to last entry in linked list */ 1026 ilink = ilink->next; 1027 } 1028 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 1029 while (ilink->previous) { 1030 ilink = ilink->previous; 1031 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 1032 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 1033 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 1034 } 1035 } 1036 } 1037 PetscFunctionReturn(0); 1038 } 1039 1040 #undef __FUNCT__ 1041 #define __FUNCT__ "PCReset_FieldSplit" 1042 static PetscErrorCode PCReset_FieldSplit(PC pc) 1043 { 1044 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1045 PetscErrorCode ierr; 1046 PC_FieldSplitLink ilink = jac->head,next; 1047 1048 PetscFunctionBegin; 1049 while (ilink) { 1050 ierr = KSPReset(ilink->ksp);CHKERRQ(ierr); 1051 ierr = VecDestroy(&ilink->x);CHKERRQ(ierr); 1052 ierr = VecDestroy(&ilink->y);CHKERRQ(ierr); 1053 ierr = VecDestroy(&ilink->z);CHKERRQ(ierr); 1054 ierr = VecScatterDestroy(&ilink->sctx);CHKERRQ(ierr); 1055 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1056 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1057 next = ilink->next; 1058 ilink = next; 1059 } 1060 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 1061 if (jac->mat && jac->mat != jac->pmat) { 1062 ierr = MatDestroyMatrices(jac->nsplits,&jac->mat);CHKERRQ(ierr); 1063 } else if (jac->mat) { 1064 jac->mat = NULL; 1065 } 1066 if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);} 1067 if (jac->Afield) {ierr = MatDestroyMatrices(jac->nsplits,&jac->Afield);CHKERRQ(ierr);} 1068 ierr = VecDestroy(&jac->w1);CHKERRQ(ierr); 1069 ierr = VecDestroy(&jac->w2);CHKERRQ(ierr); 1070 ierr = MatDestroy(&jac->schur);CHKERRQ(ierr); 1071 ierr = MatDestroy(&jac->schurp);CHKERRQ(ierr); 1072 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1073 ierr = KSPDestroy(&jac->kspschur);CHKERRQ(ierr); 1074 ierr = KSPDestroy(&jac->kspupper);CHKERRQ(ierr); 1075 ierr = MatDestroy(&jac->B);CHKERRQ(ierr); 1076 ierr = MatDestroy(&jac->C);CHKERRQ(ierr); 1077 jac->reset = PETSC_TRUE; 1078 PetscFunctionReturn(0); 1079 } 1080 1081 #undef __FUNCT__ 1082 #define __FUNCT__ "PCDestroy_FieldSplit" 1083 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 1084 { 1085 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1086 PetscErrorCode ierr; 1087 PC_FieldSplitLink ilink = jac->head,next; 1088 1089 PetscFunctionBegin; 1090 ierr = PCReset_FieldSplit(pc);CHKERRQ(ierr); 1091 while (ilink) { 1092 ierr = KSPDestroy(&ilink->ksp);CHKERRQ(ierr); 1093 next = ilink->next; 1094 ierr = PetscFree(ilink->splitname);CHKERRQ(ierr); 1095 ierr = PetscFree(ilink->fields);CHKERRQ(ierr); 1096 ierr = PetscFree(ilink->fields_col);CHKERRQ(ierr); 1097 ierr = PetscFree(ilink);CHKERRQ(ierr); 1098 ilink = next; 1099 } 1100 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 1101 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1102 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);CHKERRQ(ierr); 1103 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);CHKERRQ(ierr); 1104 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);CHKERRQ(ierr); 1105 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);CHKERRQ(ierr); 1106 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);CHKERRQ(ierr); 1107 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL);CHKERRQ(ierr); 1108 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL);CHKERRQ(ierr); 1109 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);CHKERRQ(ierr); 1110 PetscFunctionReturn(0); 1111 } 1112 1113 #undef __FUNCT__ 1114 #define __FUNCT__ "PCSetFromOptions_FieldSplit" 1115 static PetscErrorCode PCSetFromOptions_FieldSplit(PetscOptions *PetscOptionsObject,PC pc) 1116 { 1117 PetscErrorCode ierr; 1118 PetscInt bs; 1119 PetscBool flg,stokes = PETSC_FALSE; 1120 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1121 PCCompositeType ctype; 1122 1123 PetscFunctionBegin; 1124 ierr = PetscOptionsHead(PetscOptionsObject,"FieldSplit options");CHKERRQ(ierr); 1125 ierr = PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);CHKERRQ(ierr); 1126 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 1127 if (flg) { 1128 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 1129 } 1130 jac->diag_use_amat = pc->useAmat; 1131 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); 1132 jac->offdiag_use_amat = pc->useAmat; 1133 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); 1134 /* FIXME: No programmatic equivalent to the following. */ 1135 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);CHKERRQ(ierr); 1136 if (stokes) { 1137 ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);CHKERRQ(ierr); 1138 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF; 1139 } 1140 1141 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr); 1142 if (flg) { 1143 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr); 1144 } 1145 /* Only setup fields once */ 1146 if ((jac->bs > 0) && (jac->nsplits == 0)) { 1147 /* only allow user to set fields from command line if bs is already known. 1148 otherwise user can set them in PCFieldSplitSetDefaults() */ 1149 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 1150 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 1151 } 1152 if (jac->type == PC_COMPOSITE_SCHUR) { 1153 ierr = PetscOptionsGetEnum(((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);CHKERRQ(ierr); 1154 if (flg) {ierr = PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");CHKERRQ(ierr);} 1155 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); 1156 ierr = PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);CHKERRQ(ierr); 1157 } 1158 ierr = PetscOptionsTail();CHKERRQ(ierr); 1159 PetscFunctionReturn(0); 1160 } 1161 1162 /*------------------------------------------------------------------------------------*/ 1163 1164 #undef __FUNCT__ 1165 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 1166 static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1167 { 1168 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1169 PetscErrorCode ierr; 1170 PC_FieldSplitLink ilink,next = jac->head; 1171 char prefix[128]; 1172 PetscInt i; 1173 1174 PetscFunctionBegin; 1175 if (jac->splitdefined) { 1176 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1177 PetscFunctionReturn(0); 1178 } 1179 for (i=0; i<n; i++) { 1180 if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 1181 if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 1182 } 1183 ierr = PetscNew(&ilink);CHKERRQ(ierr); 1184 if (splitname) { 1185 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1186 } else { 1187 ierr = PetscMalloc1(3,&ilink->splitname);CHKERRQ(ierr); 1188 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr); 1189 } 1190 ierr = PetscMalloc1(n,&ilink->fields);CHKERRQ(ierr); 1191 ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 1192 ierr = PetscMalloc1(n,&ilink->fields_col);CHKERRQ(ierr); 1193 ierr = PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));CHKERRQ(ierr); 1194 1195 ilink->nfields = n; 1196 ilink->next = NULL; 1197 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr); 1198 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1199 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1200 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1201 1202 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr); 1203 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1204 1205 if (!next) { 1206 jac->head = ilink; 1207 ilink->previous = NULL; 1208 } else { 1209 while (next->next) { 1210 next = next->next; 1211 } 1212 next->next = ilink; 1213 ilink->previous = next; 1214 } 1215 jac->nsplits++; 1216 PetscFunctionReturn(0); 1217 } 1218 1219 #undef __FUNCT__ 1220 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit_Schur" 1221 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 1222 { 1223 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1224 PetscErrorCode ierr; 1225 1226 PetscFunctionBegin; 1227 ierr = PetscMalloc1(jac->nsplits,subksp);CHKERRQ(ierr); 1228 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr); 1229 1230 (*subksp)[1] = jac->kspschur; 1231 if (n) *n = jac->nsplits; 1232 PetscFunctionReturn(0); 1233 } 1234 1235 #undef __FUNCT__ 1236 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 1237 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1238 { 1239 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1240 PetscErrorCode ierr; 1241 PetscInt cnt = 0; 1242 PC_FieldSplitLink ilink = jac->head; 1243 1244 PetscFunctionBegin; 1245 ierr = PetscMalloc1(jac->nsplits,subksp);CHKERRQ(ierr); 1246 while (ilink) { 1247 (*subksp)[cnt++] = ilink->ksp; 1248 ilink = ilink->next; 1249 } 1250 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); 1251 if (n) *n = jac->nsplits; 1252 PetscFunctionReturn(0); 1253 } 1254 1255 #undef __FUNCT__ 1256 #define __FUNCT__ "PCFieldSplitSetIS_FieldSplit" 1257 static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 1258 { 1259 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1260 PetscErrorCode ierr; 1261 PC_FieldSplitLink ilink, next = jac->head; 1262 char prefix[128]; 1263 1264 PetscFunctionBegin; 1265 if (jac->splitdefined) { 1266 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1267 PetscFunctionReturn(0); 1268 } 1269 ierr = PetscNew(&ilink);CHKERRQ(ierr); 1270 if (splitname) { 1271 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1272 } else { 1273 ierr = PetscMalloc1(8,&ilink->splitname);CHKERRQ(ierr); 1274 ierr = PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);CHKERRQ(ierr); 1275 } 1276 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1277 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1278 ilink->is = is; 1279 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1280 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1281 ilink->is_col = is; 1282 ilink->next = NULL; 1283 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr); 1284 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1285 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1286 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1287 1288 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr); 1289 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1290 1291 if (!next) { 1292 jac->head = ilink; 1293 ilink->previous = NULL; 1294 } else { 1295 while (next->next) { 1296 next = next->next; 1297 } 1298 next->next = ilink; 1299 ilink->previous = next; 1300 } 1301 jac->nsplits++; 1302 PetscFunctionReturn(0); 1303 } 1304 1305 #undef __FUNCT__ 1306 #define __FUNCT__ "PCFieldSplitSetFields" 1307 /*@ 1308 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 1309 1310 Logically Collective on PC 1311 1312 Input Parameters: 1313 + pc - the preconditioner context 1314 . splitname - name of this split, if NULL the number of the split is used 1315 . n - the number of fields in this split 1316 - fields - the fields in this split 1317 1318 Level: intermediate 1319 1320 Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 1321 1322 The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block 1323 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 1324 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 1325 where the numbered entries indicate what is in the field. 1326 1327 This function is called once per split (it creates a new split each time). Solve options 1328 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1329 1330 Developer Note: This routine does not actually create the IS representing the split, that is delayed 1331 until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be 1332 available when this routine is called. 1333 1334 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 1335 1336 @*/ 1337 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1338 { 1339 PetscErrorCode ierr; 1340 1341 PetscFunctionBegin; 1342 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1343 PetscValidCharPointer(splitname,2); 1344 if (n < 1) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 1345 PetscValidIntPointer(fields,3); 1346 ierr = PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt*,const PetscInt*),(pc,splitname,n,fields,fields_col));CHKERRQ(ierr); 1347 PetscFunctionReturn(0); 1348 } 1349 1350 #undef __FUNCT__ 1351 #define __FUNCT__ "PCFieldSplitSetDiagUseAmat" 1352 /*@ 1353 PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) 1354 1355 Logically Collective on PC 1356 1357 Input Parameters: 1358 + pc - the preconditioner object 1359 - flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from 1360 1361 Options Database: 1362 . -pc_fieldsplit_diag_use_amat 1363 1364 Level: intermediate 1365 1366 .seealso: PCFieldSplitGetDiagUseAmat(), PCFieldSplitSetOffDiagUseAmat(), PCFIELDSPLIT 1367 1368 @*/ 1369 PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg) 1370 { 1371 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1372 PetscBool isfs; 1373 PetscErrorCode ierr; 1374 1375 PetscFunctionBegin; 1376 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1377 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1378 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1379 jac->diag_use_amat = flg; 1380 PetscFunctionReturn(0); 1381 } 1382 1383 #undef __FUNCT__ 1384 #define __FUNCT__ "PCFieldSplitGetDiagUseAmat" 1385 /*@ 1386 PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) 1387 1388 Logically Collective on PC 1389 1390 Input Parameters: 1391 . pc - the preconditioner object 1392 1393 Output Parameters: 1394 . flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from 1395 1396 1397 Level: intermediate 1398 1399 .seealso: PCFieldSplitSetDiagUseAmat(), PCFieldSplitGetOffDiagUseAmat(), PCFIELDSPLIT 1400 1401 @*/ 1402 PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg) 1403 { 1404 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1405 PetscBool isfs; 1406 PetscErrorCode ierr; 1407 1408 PetscFunctionBegin; 1409 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1410 PetscValidPointer(flg,2); 1411 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1412 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1413 *flg = jac->diag_use_amat; 1414 PetscFunctionReturn(0); 1415 } 1416 1417 #undef __FUNCT__ 1418 #define __FUNCT__ "PCFieldSplitSetOffDiagUseAmat" 1419 /*@ 1420 PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) 1421 1422 Logically Collective on PC 1423 1424 Input Parameters: 1425 + pc - the preconditioner object 1426 - flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 1427 1428 Options Database: 1429 . -pc_fieldsplit_off_diag_use_amat 1430 1431 Level: intermediate 1432 1433 .seealso: PCFieldSplitGetOffDiagUseAmat(), PCFieldSplitSetDiagUseAmat(), PCFIELDSPLIT 1434 1435 @*/ 1436 PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg) 1437 { 1438 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1439 PetscBool isfs; 1440 PetscErrorCode ierr; 1441 1442 PetscFunctionBegin; 1443 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1444 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1445 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1446 jac->offdiag_use_amat = flg; 1447 PetscFunctionReturn(0); 1448 } 1449 1450 #undef __FUNCT__ 1451 #define __FUNCT__ "PCFieldSplitGetOffDiagUseAmat" 1452 /*@ 1453 PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) 1454 1455 Logically Collective on PC 1456 1457 Input Parameters: 1458 . pc - the preconditioner object 1459 1460 Output Parameters: 1461 . flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 1462 1463 1464 Level: intermediate 1465 1466 .seealso: PCFieldSplitSetOffDiagUseAmat(), PCFieldSplitGetDiagUseAmat(), PCFIELDSPLIT 1467 1468 @*/ 1469 PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg) 1470 { 1471 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1472 PetscBool isfs; 1473 PetscErrorCode ierr; 1474 1475 PetscFunctionBegin; 1476 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1477 PetscValidPointer(flg,2); 1478 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1479 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1480 *flg = jac->offdiag_use_amat; 1481 PetscFunctionReturn(0); 1482 } 1483 1484 1485 1486 #undef __FUNCT__ 1487 #define __FUNCT__ "PCFieldSplitSetIS" 1488 /*@ 1489 PCFieldSplitSetIS - Sets the exact elements for field 1490 1491 Logically Collective on PC 1492 1493 Input Parameters: 1494 + pc - the preconditioner context 1495 . splitname - name of this split, if NULL the number of the split is used 1496 - is - the index set that defines the vector elements in this field 1497 1498 1499 Notes: 1500 Use PCFieldSplitSetFields(), for fields defined by strided types. 1501 1502 This function is called once per split (it creates a new split each time). Solve options 1503 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1504 1505 Level: intermediate 1506 1507 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 1508 1509 @*/ 1510 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 1511 { 1512 PetscErrorCode ierr; 1513 1514 PetscFunctionBegin; 1515 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1516 if (splitname) PetscValidCharPointer(splitname,2); 1517 PetscValidHeaderSpecific(is,IS_CLASSID,3); 1518 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr); 1519 PetscFunctionReturn(0); 1520 } 1521 1522 #undef __FUNCT__ 1523 #define __FUNCT__ "PCFieldSplitGetIS" 1524 /*@ 1525 PCFieldSplitGetIS - Retrieves the elements for a field as an IS 1526 1527 Logically Collective on PC 1528 1529 Input Parameters: 1530 + pc - the preconditioner context 1531 - splitname - name of this split 1532 1533 Output Parameter: 1534 - is - the index set that defines the vector elements in this field, or NULL if the field is not found 1535 1536 Level: intermediate 1537 1538 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS() 1539 1540 @*/ 1541 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is) 1542 { 1543 PetscErrorCode ierr; 1544 1545 PetscFunctionBegin; 1546 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1547 PetscValidCharPointer(splitname,2); 1548 PetscValidPointer(is,3); 1549 { 1550 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 1551 PC_FieldSplitLink ilink = jac->head; 1552 PetscBool found; 1553 1554 *is = NULL; 1555 while (ilink) { 1556 ierr = PetscStrcmp(ilink->splitname, splitname, &found);CHKERRQ(ierr); 1557 if (found) { 1558 *is = ilink->is; 1559 break; 1560 } 1561 ilink = ilink->next; 1562 } 1563 } 1564 PetscFunctionReturn(0); 1565 } 1566 1567 #undef __FUNCT__ 1568 #define __FUNCT__ "PCFieldSplitSetBlockSize" 1569 /*@ 1570 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 1571 fieldsplit preconditioner. If not set the matrix block size is used. 1572 1573 Logically Collective on PC 1574 1575 Input Parameters: 1576 + pc - the preconditioner context 1577 - bs - the block size 1578 1579 Level: intermediate 1580 1581 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 1582 1583 @*/ 1584 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 1585 { 1586 PetscErrorCode ierr; 1587 1588 PetscFunctionBegin; 1589 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1590 PetscValidLogicalCollectiveInt(pc,bs,2); 1591 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr); 1592 PetscFunctionReturn(0); 1593 } 1594 1595 #undef __FUNCT__ 1596 #define __FUNCT__ "PCFieldSplitGetSubKSP" 1597 /*@C 1598 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1599 1600 Collective on KSP 1601 1602 Input Parameter: 1603 . pc - the preconditioner context 1604 1605 Output Parameters: 1606 + n - the number of splits 1607 - pc - the array of KSP contexts 1608 1609 Note: 1610 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 1611 (not the KSP just the array that contains them). 1612 1613 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 1614 1615 Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs. 1616 You can call PCFieldSplitGetSubKSP(pc,n,NULL_OBJECT,ierr) to determine how large the 1617 KSP array must be. 1618 1619 1620 Level: advanced 1621 1622 .seealso: PCFIELDSPLIT 1623 @*/ 1624 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1625 { 1626 PetscErrorCode ierr; 1627 1628 PetscFunctionBegin; 1629 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1630 if (n) PetscValidIntPointer(n,2); 1631 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 1632 PetscFunctionReturn(0); 1633 } 1634 1635 #undef __FUNCT__ 1636 #define __FUNCT__ "PCFieldSplitSetSchurPre" 1637 /*@ 1638 PCFieldSplitSetSchurPre - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the 1639 A11 matrix. Otherwise no preconditioner is used. 1640 1641 Collective on PC 1642 1643 Input Parameters: 1644 + pc - the preconditioner context 1645 . 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 1646 - userpre - matrix to use for preconditioning, or NULL 1647 1648 Options Database: 1649 . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> default is a11 1650 1651 Notes: 1652 $ If ptype is 1653 $ user then the preconditioner for the Schur complement is generated by the provided matrix (pre argument 1654 $ to this function). 1655 $ a11 then the preconditioner for the Schur complement is generated by the block diagonal part of the preconditioner 1656 $ matrix associated with the Schur complement (i.e. A11), not he Schur complement matrix 1657 $ full then the preconditioner uses the exact Schur complement (this is expensive) 1658 $ self the preconditioner for the Schur complement is generated from the Schur complement matrix itself: 1659 $ The only preconditioner that currently works directly with the Schur complement matrix object is the PCLSC 1660 $ preconditioner 1661 $ selfp then the preconditioning matrix is an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01 1662 $ This is only a good preconditioner when diag(A00) is a good preconditioner for A00. 1663 1664 When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense 1665 with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and 1666 -fieldsplit_1_pc_type lsc which uses the least squares commutator to compute a preconditioner for the Schur complement. 1667 1668 Level: intermediate 1669 1670 .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC 1671 1672 @*/ 1673 PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1674 { 1675 PetscErrorCode ierr; 1676 1677 PetscFunctionBegin; 1678 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1679 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr); 1680 PetscFunctionReturn(0); 1681 } 1682 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */ 1683 1684 #undef __FUNCT__ 1685 #define __FUNCT__ "PCFieldSplitGetSchurPre" 1686 /*@ 1687 PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be 1688 preconditioned. See PCFieldSplitSetSchurPre() for details. 1689 1690 Logically Collective on PC 1691 1692 Input Parameters: 1693 . pc - the preconditioner context 1694 1695 Output Parameters: 1696 + ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11, PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_PRE_USER 1697 - userpre - matrix to use for preconditioning (with PC_FIELDSPLIT_PRE_USER), or NULL 1698 1699 Level: intermediate 1700 1701 .seealso: PCFieldSplitSetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC 1702 1703 @*/ 1704 PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre) 1705 { 1706 PetscErrorCode ierr; 1707 1708 PetscFunctionBegin; 1709 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1710 ierr = PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));CHKERRQ(ierr); 1711 PetscFunctionReturn(0); 1712 } 1713 1714 #undef __FUNCT__ 1715 #define __FUNCT__ "PCFieldSplitSchurGetS" 1716 /*@ 1717 PCFieldSplitSchurGetS - extract the MatSchurComplement object used by this PC in case it needs to be configured separately 1718 1719 Not collective 1720 1721 Input Parameter: 1722 . pc - the preconditioner context 1723 1724 Output Parameter: 1725 . S - the Schur complement matrix 1726 1727 Notes: 1728 This matrix should not be destroyed using MatDestroy(); rather, use PCFieldSplitSchurRestoreS(). 1729 1730 Level: advanced 1731 1732 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurRestoreS() 1733 1734 @*/ 1735 PetscErrorCode PCFieldSplitSchurGetS(PC pc,Mat *S) 1736 { 1737 PetscErrorCode ierr; 1738 const char* t; 1739 PetscBool isfs; 1740 PC_FieldSplit *jac; 1741 1742 PetscFunctionBegin; 1743 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1744 ierr = PetscObjectGetType((PetscObject)pc,&t);CHKERRQ(ierr); 1745 ierr = PetscStrcmp(t,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1746 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 1747 jac = (PC_FieldSplit*)pc->data; 1748 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type); 1749 if (S) *S = jac->schur; 1750 PetscFunctionReturn(0); 1751 } 1752 1753 #undef __FUNCT__ 1754 #define __FUNCT__ "PCFieldSplitSchurRestoreS" 1755 /*@ 1756 PCFieldSplitSchurRestoreS - restores the MatSchurComplement object used by this PC 1757 1758 Not collective 1759 1760 Input Parameters: 1761 + pc - the preconditioner context 1762 . S - the Schur complement matrix 1763 1764 Level: advanced 1765 1766 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurGetS() 1767 1768 @*/ 1769 PetscErrorCode PCFieldSplitSchurRestoreS(PC pc,Mat *S) 1770 { 1771 PetscErrorCode ierr; 1772 const char* t; 1773 PetscBool isfs; 1774 PC_FieldSplit *jac; 1775 1776 PetscFunctionBegin; 1777 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1778 ierr = PetscObjectGetType((PetscObject)pc,&t);CHKERRQ(ierr); 1779 ierr = PetscStrcmp(t,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1780 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 1781 jac = (PC_FieldSplit*)pc->data; 1782 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type); 1783 if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten"); 1784 PetscFunctionReturn(0); 1785 } 1786 1787 1788 #undef __FUNCT__ 1789 #define __FUNCT__ "PCFieldSplitSetSchurPre_FieldSplit" 1790 static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1791 { 1792 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1793 PetscErrorCode ierr; 1794 1795 PetscFunctionBegin; 1796 jac->schurpre = ptype; 1797 if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) { 1798 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1799 jac->schur_user = pre; 1800 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 1801 } 1802 PetscFunctionReturn(0); 1803 } 1804 1805 #undef __FUNCT__ 1806 #define __FUNCT__ "PCFieldSplitGetSchurPre_FieldSplit" 1807 static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre) 1808 { 1809 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1810 1811 PetscFunctionBegin; 1812 *ptype = jac->schurpre; 1813 *pre = jac->schur_user; 1814 PetscFunctionReturn(0); 1815 } 1816 1817 #undef __FUNCT__ 1818 #define __FUNCT__ "PCFieldSplitSetSchurFactType" 1819 /*@ 1820 PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain 1821 1822 Collective on PC 1823 1824 Input Parameters: 1825 + pc - the preconditioner context 1826 - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default 1827 1828 Options Database: 1829 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full 1830 1831 1832 Level: intermediate 1833 1834 Notes: 1835 The FULL factorization is 1836 1837 $ (A B) = (1 0) (A 0) (1 Ainv*B) 1838 $ (C D) (C*Ainv 1) (0 S) (0 1 ) 1839 1840 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, 1841 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). 1842 1843 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 1844 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 1845 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, 1846 the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations. 1847 1848 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 1849 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). 1850 1851 References: 1852 Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000) pp. 1969-1972. 1853 1854 Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001), pp. 1050-1051. 1855 1856 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType 1857 @*/ 1858 PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype) 1859 { 1860 PetscErrorCode ierr; 1861 1862 PetscFunctionBegin; 1863 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1864 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));CHKERRQ(ierr); 1865 PetscFunctionReturn(0); 1866 } 1867 1868 #undef __FUNCT__ 1869 #define __FUNCT__ "PCFieldSplitSetSchurFactType_FieldSplit" 1870 static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype) 1871 { 1872 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1873 1874 PetscFunctionBegin; 1875 jac->schurfactorization = ftype; 1876 PetscFunctionReturn(0); 1877 } 1878 1879 #undef __FUNCT__ 1880 #define __FUNCT__ "PCFieldSplitGetSchurBlocks" 1881 /*@C 1882 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 1883 1884 Collective on KSP 1885 1886 Input Parameter: 1887 . pc - the preconditioner context 1888 1889 Output Parameters: 1890 + A00 - the (0,0) block 1891 . A01 - the (0,1) block 1892 . A10 - the (1,0) block 1893 - A11 - the (1,1) block 1894 1895 Level: advanced 1896 1897 .seealso: PCFIELDSPLIT 1898 @*/ 1899 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 1900 { 1901 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 1902 1903 PetscFunctionBegin; 1904 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1905 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 1906 if (A00) *A00 = jac->pmat[0]; 1907 if (A01) *A01 = jac->B; 1908 if (A10) *A10 = jac->C; 1909 if (A11) *A11 = jac->pmat[1]; 1910 PetscFunctionReturn(0); 1911 } 1912 1913 #undef __FUNCT__ 1914 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 1915 static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 1916 { 1917 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1918 PetscErrorCode ierr; 1919 1920 PetscFunctionBegin; 1921 jac->type = type; 1922 if (type == PC_COMPOSITE_SCHUR) { 1923 pc->ops->apply = PCApply_FieldSplit_Schur; 1924 pc->ops->view = PCView_FieldSplit_Schur; 1925 1926 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 1927 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);CHKERRQ(ierr); 1928 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);CHKERRQ(ierr); 1929 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);CHKERRQ(ierr); 1930 1931 } else { 1932 pc->ops->apply = PCApply_FieldSplit; 1933 pc->ops->view = PCView_FieldSplit; 1934 1935 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1936 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);CHKERRQ(ierr); 1937 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);CHKERRQ(ierr); 1938 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);CHKERRQ(ierr); 1939 } 1940 PetscFunctionReturn(0); 1941 } 1942 1943 #undef __FUNCT__ 1944 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 1945 static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 1946 { 1947 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1948 1949 PetscFunctionBegin; 1950 if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 1951 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); 1952 jac->bs = bs; 1953 PetscFunctionReturn(0); 1954 } 1955 1956 #undef __FUNCT__ 1957 #define __FUNCT__ "PCFieldSplitSetType" 1958 /*@ 1959 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 1960 1961 Collective on PC 1962 1963 Input Parameter: 1964 . pc - the preconditioner context 1965 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1966 1967 Options Database Key: 1968 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 1969 1970 Level: Intermediate 1971 1972 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1973 1974 .seealso: PCCompositeSetType() 1975 1976 @*/ 1977 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 1978 { 1979 PetscErrorCode ierr; 1980 1981 PetscFunctionBegin; 1982 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1983 ierr = PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));CHKERRQ(ierr); 1984 PetscFunctionReturn(0); 1985 } 1986 1987 #undef __FUNCT__ 1988 #define __FUNCT__ "PCFieldSplitGetType" 1989 /*@ 1990 PCFieldSplitGetType - Gets the type of fieldsplit preconditioner. 1991 1992 Not collective 1993 1994 Input Parameter: 1995 . pc - the preconditioner context 1996 1997 Output Parameter: 1998 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1999 2000 Level: Intermediate 2001 2002 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 2003 .seealso: PCCompositeSetType() 2004 @*/ 2005 PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type) 2006 { 2007 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2008 2009 PetscFunctionBegin; 2010 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2011 PetscValidIntPointer(type,2); 2012 *type = jac->type; 2013 PetscFunctionReturn(0); 2014 } 2015 2016 #undef __FUNCT__ 2017 #define __FUNCT__ "PCFieldSplitSetDMSplits" 2018 /*@ 2019 PCFieldSplitSetDMSplits - Flags whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 2020 2021 Logically Collective 2022 2023 Input Parameters: 2024 + pc - the preconditioner context 2025 - flg - boolean indicating whether to use field splits defined by the DM 2026 2027 Options Database Key: 2028 . -pc_fieldsplit_dm_splits 2029 2030 Level: Intermediate 2031 2032 .keywords: PC, DM, composite preconditioner, additive, multiplicative 2033 2034 .seealso: PCFieldSplitGetDMSplits() 2035 2036 @*/ 2037 PetscErrorCode PCFieldSplitSetDMSplits(PC pc,PetscBool flg) 2038 { 2039 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2040 PetscBool isfs; 2041 PetscErrorCode ierr; 2042 2043 PetscFunctionBegin; 2044 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2045 PetscValidLogicalCollectiveBool(pc,flg,2); 2046 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 2047 if (isfs) { 2048 jac->dm_splits = flg; 2049 } 2050 PetscFunctionReturn(0); 2051 } 2052 2053 2054 #undef __FUNCT__ 2055 #define __FUNCT__ "PCFieldSplitGetDMSplits" 2056 /*@ 2057 PCFieldSplitGetDMSplits - Returns flag indicating whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 2058 2059 Logically Collective 2060 2061 Input Parameter: 2062 . pc - the preconditioner context 2063 2064 Output Parameter: 2065 . flg - boolean indicating whether to use field splits defined by the DM 2066 2067 Level: Intermediate 2068 2069 .keywords: PC, DM, composite preconditioner, additive, multiplicative 2070 2071 .seealso: PCFieldSplitSetDMSplits() 2072 2073 @*/ 2074 PetscErrorCode PCFieldSplitGetDMSplits(PC pc,PetscBool* flg) 2075 { 2076 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2077 PetscBool isfs; 2078 PetscErrorCode ierr; 2079 2080 PetscFunctionBegin; 2081 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2082 PetscValidPointer(flg,2); 2083 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 2084 if (isfs) { 2085 if(flg) *flg = jac->dm_splits; 2086 } 2087 PetscFunctionReturn(0); 2088 } 2089 2090 /* -------------------------------------------------------------------------------------*/ 2091 /*MC 2092 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 2093 fields or groups of fields. See the users manual section "Solving Block Matrices" for more details. 2094 2095 To set options on the solvers for each block append -fieldsplit_ to all the PC 2096 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 2097 2098 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 2099 and set the options directly on the resulting KSP object 2100 2101 Level: intermediate 2102 2103 Options Database Keys: 2104 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 2105 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 2106 been supplied explicitly by -pc_fieldsplit_%d_fields 2107 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 2108 . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur> - type of relaxation or factorization splitting 2109 . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11 2110 . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver 2111 2112 - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ 2113 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 2114 2115 Notes: 2116 Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 2117 to define a field by an arbitrary collection of entries. 2118 2119 If no fields are set the default is used. The fields are defined by entries strided by bs, 2120 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 2121 if this is not called the block size defaults to the blocksize of the second matrix passed 2122 to KSPSetOperators()/PCSetOperators(). 2123 2124 $ For the Schur complement preconditioner if J = ( A00 A01 ) 2125 $ ( A10 A11 ) 2126 $ the preconditioner using full factorization is 2127 $ ( I -ksp(A00) A01 ) ( inv(A00) 0 ) ( I 0 ) 2128 $ ( 0 I ) ( 0 ksp(S) ) ( -A10 ksp(A00) I ) 2129 where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. S is the Schur complement 2130 $ S = A11 - A10 ksp(A00) A01 2131 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 2132 in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0, 2133 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 2134 A11 is used to construct a preconditioner for S, use PCFieldSplitSetSchurPre() to turn on or off this 2135 option. You can use the preconditioner PCLSC to precondition the Schur complement with -fieldsplit_1_pc_type lsc. 2136 When option -fieldsplit_schur_precondition selfp is given, an approximation to S is assembled -- 2137 Sp = A11 - A10 inv(diag(A00)) A01, which has type AIJ and can be used with a variety of preconditioners 2138 (e.g., -fieldsplit_1_pc_type asm). 2139 The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above, 2140 diag gives 2141 $ ( inv(A00) 0 ) 2142 $ ( 0 -ksp(S) ) 2143 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 2144 $ ( A00 0 ) 2145 $ ( A10 S ) 2146 where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of 2147 $ ( A00 A01 ) 2148 $ ( 0 S ) 2149 where again the inverses of A00 and S are applied using KSPs. 2150 2151 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 2152 is used automatically for a second block. 2153 2154 The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1. 2155 Generally it should be used with the AIJ format. 2156 2157 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 2158 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT 2159 inside a smoother resulting in "Distributive Smoothers". 2160 2161 Concepts: physics based preconditioners, block preconditioners 2162 2163 There is a nice discussion of block preconditioners in 2164 2165 [El08] A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations 2166 Howard Elman, V.E. Howle, John Shadid, Robert Shuttleworth, Ray Tuminaro, Journal of Computational Physics 227 (2008) 1790--1808 2167 http://chess.cs.umd.edu/~elman/papers/tax.pdf 2168 2169 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC, 2170 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre() 2171 M*/ 2172 2173 #undef __FUNCT__ 2174 #define __FUNCT__ "PCCreate_FieldSplit" 2175 PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc) 2176 { 2177 PetscErrorCode ierr; 2178 PC_FieldSplit *jac; 2179 2180 PetscFunctionBegin; 2181 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr); 2182 2183 jac->bs = -1; 2184 jac->nsplits = 0; 2185 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 2186 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 2187 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL; 2188 jac->dm_splits = PETSC_TRUE; 2189 2190 pc->data = (void*)jac; 2191 2192 pc->ops->apply = PCApply_FieldSplit; 2193 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 2194 pc->ops->setup = PCSetUp_FieldSplit; 2195 pc->ops->reset = PCReset_FieldSplit; 2196 pc->ops->destroy = PCDestroy_FieldSplit; 2197 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 2198 pc->ops->view = PCView_FieldSplit; 2199 pc->ops->applyrichardson = 0; 2200 2201 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 2202 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 2203 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 2204 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 2205 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 2206 PetscFunctionReturn(0); 2207 } 2208 2209 2210 2211